jacobian

Unnamed repository; edit this file 'description' to name the repository.
Log | Files | Refs | README

commit 17ace7c4da8fd6e11435e2ebac925bd6ccf44f4e
parent 9fc32c95da5344e8b735b0fecae7e632d47c7aa6
Author: David Freifeld <freifeld.david@gmail.com>
Date:   Fri,  9 Apr 2021 20:41:42 -0700

Move optimizers to util.hpp, consolidate namespaces, fix python demo

Diffstat:
Mexample.cpp | 44++++++++++++++++++++++----------------------
Mexample.py | 28++--------------------------
Msrc/bpnn.cpp | 29+++++++++++++++++------------
Msrc/bpnn.hpp | 40++++++++++++++++++----------------------
Dsrc/optimizers.cpp | 50--------------------------------------------------
Msrc/pybind.cpp | 45+++++++++++++++++++++++----------------------
Msrc/utils.cpp | 89+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------
Msrc/utils.hpp | 24+++++++++++++++++++++---
8 files changed, 172 insertions(+), 177 deletions(-)

diff --git a/example.cpp b/example.cpp @@ -10,36 +10,36 @@ // #include <indicators/block_progress_bar.hpp> // using namespace indicators; -#include "./src/bpnn.hpp" -#include "./src/utils.hpp" +#include "src/bpnn.hpp" +#include "src/utils.hpp" #include "unistd.h" #include <ctime> #include <chrono> double bench(int batch_sz, int epochs) { - auto start = std::chrono::high_resolution_clock::now(); - Network net ("./data_banknote_authentication.txt", batch_sz, 0.0155, 0.03, L2, 0, 0.9); - net.add_layer(4, linear, linear_deriv); - net.add_layer(5, lecun_tanh, lecun_tanh_deriv); - net.add_layer(2, linear, linear_deriv); - // net.init_optimizer(optimizers::demon(0.1, 50)); - net.initialize(); - for (int i = 0; i < epochs; i++) { - net.train(); - } - auto end = std::chrono::high_resolution_clock::now(); - return std::chrono::duration_cast<std::chrono::nanoseconds>(end - start).count() / pow(10,9); + auto start = std::chrono::high_resolution_clock::now(); + Jacobian::Network net ("./data_banknote_authentication.txt", batch_sz, 0.0155, 0.03, Jacobian::Regularization::L2, 0, 0.9); + net.add_layer(4, Jacobian::activations::linear, Jacobian::activations::linear_deriv); + net.add_layer(5, Jacobian::activations::lecun_tanh, Jacobian::activations::lecun_tanh_deriv); + net.add_layer(2, Jacobian::activations::linear, Jacobian::activations::linear_deriv); + net.init_optimizer(Jacobian::optimizers::momentum(0.1)); + net.initialize(); + for (int i = 0; i < epochs; i++) { + net.train(); + } + auto end = std::chrono::high_resolution_clock::now(); + return std::chrono::duration_cast<std::chrono::nanoseconds>(end - start).count() / pow(10,9); } int main(int argc, char** argv) { - if (argc < 2) { - std::cout << "Invalid command! Either pass a special option or pass two integers - batch_size and epochs (in that order)." << "\n"; - exit(1); - } - else { - sleep(strtol(argv[3], NULL, 10)); - std::cout << bench(strtol(argv[1], NULL, 10), strtol(argv[2], NULL, 10)) << "\n"; - } + if (argc < 2) { + std::cout << "Invalid command! Either pass a special option or pass two integers - batch_size and epochs (in that order)." << "\n"; + exit(1); + } + else { + sleep(strtol(argv[3], NULL, 10)); + std::cout << bench(strtol(argv[1], NULL, 10), strtol(argv[2], NULL, 10)) << "\n"; + } } diff --git a/example.py b/example.py @@ -10,29 +10,5 @@ net.add_layer(2, jcb.activations.linear, jcb.activations.linear_deriv) net.init_optimizer(jcb.optimizers.momentum(0.1)) net.init_decay(jcb.decays.exponential(1, 0.5)) net.initialize() -fig, axs = plt.subplots() -accuracies = [] -costs = [] -def update(frame): - net.feedforward() - net.backpropagate() - for i in range(3): - plt.subplot(3, 4, 4*i + 1) - plt.imshow(net.layers[i].get_contents()) - plt.colorbar() - plt.subplot(3, 4, 4*i + 2) - plt.imshow(net.layers[i].get_weights()) - plt.colorbar() - plt.subplot(3, 4, 4*i + 3) - plt.imshow(net.layers[i].get_bias()) - plt.colorbar() - accuracies.append(net.accuracy()) - costs.append(net.cost()) - plt.subplot(3,4,4) - plt.plot(costs[1:], color="blue") - plt.subplot(3,4,8) - plt.plot(accuracies[1:], color="orange") - time.sleep(0.01) - net.next_batch() -ani = FuncAnimation(fig, update, interval=1) -plt.show() +for i in range(50): + net.train() diff --git a/src/bpnn.cpp b/src/bpnn.cpp @@ -9,6 +9,7 @@ #include "utils.hpp" #include <random> +namespace Jacobian { Layer::Layer(int batch_sz, int nodes) { contents = Eigen::MatrixXf(batch_sz, nodes); @@ -205,8 +206,8 @@ float Network::cost() sum-=tempsum; } for (unsigned long i = 0; i < layers.size()-1; i++) { - if (reg_type == L2) reg += layers[i].weights.cwiseProduct(layers[i].weights).sum(); - else if (reg_type == L1) reg += (layers[i].weights.array().abs().matrix()).sum(); + if (reg_type == Regularization::L2) reg += layers[i].weights.cwiseProduct(layers[i].weights).sum(); + else if (reg_type == Regularization::L1) reg += (layers[i].weights.array().abs().matrix()).sum(); } return ((1.0/batch_size) * sum) + (1/2*lambda*reg); } @@ -266,8 +267,8 @@ Eigen::MatrixXf Network::backpropagate() } for (int i = 0; i < length-1; i++) { update(layers[length-2-i], deltas[i], learning_rate); - if (reg_type == L2) layers[length-2-i].weights -= ((lambda/batch_size) * (layers[length-2-i].weights)); - else if (reg_type == L1) layers[length-2-i].weights -= ((lambda/(2*batch_size)) * l1_deriv(layers[length-2-i].weights)); + if (reg_type == Regularization::L2) layers[length-2-i].weights -= ((lambda/batch_size) * (layers[length-2-i].weights)); + else if (reg_type == Regularization::L1) layers[length-2-i].weights -= ((lambda/(2*batch_size)) * l1_deriv(layers[length-2-i].weights)); layers[length-1-i].bias -= bias_lr * gradients[i]; } return gradients.back(); @@ -292,8 +293,6 @@ void Network::validate(const char* path) Ensures(lseek(val_data, 0, SEEK_CUR) == 0); } -#include "optimizers.cpp" - void Network::interactive_next_batch() { if (batches < instances/batch_size-batch_size) next_batch(data); @@ -309,21 +308,27 @@ void Network::train() { float cost_sum = 0; float acc_sum = 0; - for (int i = 0; i <= instances-batch_size; i+=batch_size) { - if (i != instances-batch_size) next_batch(data); + for (int i = 0; i <= instances - batch_size; i += batch_size) { + if (i != instances - batch_size) + next_batch(data); feedforward(); backpropagate(); cost_sum += cost(); acc_sum += accuracy(); batches++; } - epoch_acc = 1.0/(static_cast<float>(instances/batch_size)) * acc_sum; - epoch_cost = 1.0/(static_cast<float>(instances/batch_size)) * cost_sum; + epoch_acc = + 1.0 / (static_cast<float>(instances / batch_size)) * acc_sum; + epoch_cost = + 1.0 / (static_cast<float>(instances / batch_size)) * cost_sum; validate(VAL_PATH); - if (silenced == false) printf("Epoch %i complete - cost %f - acc %f - val_cost %f - val_acc %f\n", epochs, epoch_cost, epoch_acc, val_cost, val_acc); - batches=1; + if (silenced == false) + printf("Epoch %i complete - cost %f - acc %f - val_cost %f - val_acc %f\n", + epochs, epoch_cost, epoch_acc, val_cost, val_acc); + batches = 1; data = open(TRAIN_BIN_PATH, O_RDONLY | O_NONBLOCK); decay(learning_rate); epochs++; Ensures(lseek(data, 0, SEEK_CUR) == 0); } +} diff --git a/src/bpnn.hpp b/src/bpnn.hpp @@ -13,9 +13,10 @@ #include <fcntl.h> #include <unistd.h> +namespace Jacobian { #define BUFFER_SIZE 600*1024 #define LARGE_BUF 600*1024*15 -enum Regularization {L1, L2}; +enum class Regularization {L1, L2}; class Layer { public: @@ -73,7 +74,10 @@ public: ~Network(); void add_layer(int nodes, std::function<float(float)> activation, std::function<float(float)> activation_deriv); void initialize(); - void init_optimizer(std::function<void(Layer&, Eigen::MatrixXf, float)> f); + void init_optimizer(std::function<void(Layer &, Eigen::MatrixXf, float)> f) + { + update = f; + }; void init_decay(std::function<void(float&)> f); void set_activation(int index, std::function<float(float)> custom, std::function<float(float)> custom_deriv); void feedforward(); @@ -88,34 +92,25 @@ public: float get_acc() {return epoch_acc;} float get_val_acc() {return val_acc;} float get_cost() {return epoch_cost;} - float get_val_cost() {return val_cost;} + float get_val_cost() + { + return val_cost; + } }; -int prep_file(const char* path, const char* out_path); -int split_file(const char* path, int lines, float ratio); +int prep_file(const char *path, const char *out_path); +int split_file(const char *path, int lines, float ratio); -void prep(const char* rname, const char* wname); -void compress(const char* rname, const char* wname); +void prep(const char *rname, const char *wname); +void compress(const char *rname, const char *wname); Eigen::MatrixXf l1_deriv(Eigen::MatrixXf m); -namespace optimizers { -std::function<void(Layer&, Eigen::MatrixXf, float)> momentum(float beta); -std::function<void(Layer&, Eigen::MatrixXf, float)> demon(float beta_init, int max_ep); -std::function<void(Layer&, Eigen::MatrixXf, float)> adam(float beta1, float beta2, float epsilon); -std::function<void(Layer&, Eigen::MatrixXf, float)> adamax(float beta1, float beta2, float epsilon); -} - -namespace decays { -std::function<void(float&)> step(float a_0, float k); -std::function<void(float&)> exponential(float a_0, float k); -std::function<void(float&)> fractional(float a_0, float k); -std::function<void(float&)> linear(int max_ep); -} - #define MAXLINE 1024 #if (!RECKLESS) -#define checknan(x, loc) if(x==INFINITY || x==NAN || x == -INFINITY) throw ValueError("Detected NaN in operation", loc) +#define checknan(x, loc) \ + if (x == INFINITY || x == NAN || x == -INFINITY) \ + throw ValueError("Detected NaN in operation", loc) #define Expects(cond) assert(cond); #define Ensures(cond) assert(cond); #else @@ -130,4 +125,5 @@ std::function<void(float&)> linear(int max_ep); #define VAL_BIN_PATH "./test.bin" #define TRAIN_BIN_PATH "./train.bin" +} #endif /* MODULE_H */ diff --git a/src/optimizers.cpp b/src/optimizers.cpp @@ -1,50 +0,0 @@ -// -// optimizers.cpp -// Jacobian -// -// Created by David Freifeld -// - - -std::function<void(Layer&, Eigen::MatrixXf, float)> optimizers::momentum(float beta) { - return [beta](Layer& layer, const Eigen::MatrixXf delta, const float learning_rate) { - layer.weights -= (beta * layer.m) + (learning_rate * delta); - layer.m = (learning_rate * delta); - }; -} - -std::function<void(Layer&, Eigen::MatrixXf, float)> optimizers::demon(float beta, int max_ep) { - float beta_init = beta; - float prev_epoch = -1; - float epochs = 0; - return [max_ep, epochs, beta_init, beta](Layer& layer, const Eigen::MatrixXf delta, const float learning_rate) mutable { - beta = beta_init * (1-(epochs/max_ep)) / ((beta_init * (1-(epochs/max_ep))) + (1-beta_init)); - layer.weights -= (beta * layer.m) + (learning_rate * delta); - layer.m = (learning_rate * delta); - epochs++; - }; -} - -std::function<void(Layer&, Eigen::MatrixXf, float)> optimizers::adam(float beta1, float beta2, float epsilon) { - return [beta1, beta2, epsilon](Layer& layer, const Eigen::MatrixXf delta, const float learning_rate) { - layer.m = (beta1 * layer.m) + ((1-beta1)*delta); - layer.v = (beta2 * layer.v) + (1-beta2)*(delta.cwiseProduct(delta)); - layer.weights -= learning_rate * - ((layer.v.cwiseSqrt()).array()+epsilon).pow(-1).cwiseProduct(layer.m.array()).matrix(); - }; -} - -std::function<void(Layer&, Eigen::MatrixXf, float)> optimizers::adamax(float beta1, float beta2, float epsilon) { - return [beta1, beta2, epsilon](Layer& layer, const Eigen::MatrixXf delta, const float learning_rate) { - layer.m = (beta1 * layer.m) + ((1-beta1)*delta); - if ((beta2 * layer.v).sum() > delta.array().abs().sum()) layer.v = (beta2 * layer.v); - else layer.v = delta.array().abs().matrix(); - layer.weights -= learning_rate * - (layer.v.array().pow(-1).cwiseProduct(layer.m.array())).matrix(); - }; -} - -void Network::init_optimizer(std::function<void(Layer&, Eigen::MatrixXf, float)> f) -{ - update = f; -} diff --git a/src/pybind.cpp b/src/pybind.cpp @@ -12,6 +12,7 @@ #include "bpnn.hpp" #include "utils.hpp" +using namespace Jacobian; namespace py = pybind11; PYBIND11_MODULE(_jacobian, m) @@ -59,28 +60,28 @@ PYBIND11_MODULE(_jacobian, m) .def("get_val_acc", &Network::get_val_acc) .def_readonly("layers", &Network::layers); auto a = m.def_submodule("activations", "Submodule supplying built-in activation functions."); - a.def("linear", &linear, py::arg("x")); - a.def("linear_deriv", &linear_deriv, py::arg("x")); - a.def("sigmoid", &sigmoid, py::arg("x")); - a.def("sigmoid_deriv", &sigmoid_deriv, py::arg("x")); - a.def("lecun_tanh", &lecun_tanh, py::arg("x")); - a.def("lecun_tanh_deriv", &lecun_tanh_deriv, py::arg("x")); - a.def("softplus", &softplus, py::arg("x")); - a.def("softplus_deriv", &softplus_deriv, py::arg("x")); - a.def("inverse_logit", &inverse_logit, py::arg("x")); - a.def("inverse_logit_deriv", &inverse_logit_deriv, py::arg("x")); - a.def("cloglog", &cloglog, py::arg("x")); - a.def("cloglog_deriv", &cloglog_deriv, py::arg("x")); - a.def("bipolar", &bipolar, py::arg("x")); - a.def("bipolar_deriv", &bipolar_deriv, py::arg("x")); - a.def("step", &step, py::arg("x")); - a.def("step_deriv", &step_deriv, py::arg("x")); - a.def("hard_tanh", &hard_tanh, py::arg("x")); - a.def("hard_tanh_deriv", &hard_tanh_deriv, py::arg("x")); - a.def("leaky_relu", &leaky_relu, py::arg("x")); - a.def("leaky_relu_deriv", &leaky_relu_deriv, py::arg("x")); - a.def("relu", (rectifier(linear)), py::arg("x")); - a.def("relu_deriv", (rectifier(linear_deriv)), py::arg("x")); + a.def("linear", &activations::linear, py::arg("x")); + a.def("linear_deriv", &activations::linear_deriv, py::arg("x")); + a.def("sigmoid", &activations::sigmoid, py::arg("x")); + a.def("sigmoid_deriv", &activations::sigmoid_deriv, py::arg("x")); + a.def("lecun_tanh", &activations::lecun_tanh, py::arg("x")); + a.def("lecun_tanh_deriv", &activations::lecun_tanh_deriv, py::arg("x")); + a.def("softplus", &activations::softplus, py::arg("x")); + a.def("softplus_deriv", &activations::softplus_deriv, py::arg("x")); + a.def("inverse_logit", &activations::inverse_logit, py::arg("x")); + a.def("inverse_logit_deriv", &activations::inverse_logit_deriv, py::arg("x")); + a.def("cloglog", &activations::cloglog, py::arg("x")); + a.def("cloglog_deriv", &activations::cloglog_deriv, py::arg("x")); + a.def("bipolar", &activations::bipolar, py::arg("x")); + a.def("bipolar_deriv", &activations::bipolar_deriv, py::arg("x")); + a.def("step", &activations::step, py::arg("x")); + a.def("step_deriv", &activations::step_deriv, py::arg("x")); + a.def("hard_tanh", &activations::hard_tanh, py::arg("x")); + a.def("hard_tanh_deriv", &activations::hard_tanh_deriv, py::arg("x")); + a.def("leaky_relu", &activations::leaky_relu, py::arg("x")); + a.def("leaky_relu_deriv", &activations::leaky_relu_deriv, py::arg("x")); + a.def("relu", (activations::rectifier(activations::linear)), py::arg("x")); + a.def("relu_deriv", (activations::rectifier(activations::linear_deriv)), py::arg("x")); auto o = m.def_submodule("optimizers", "Submodule supplying built-in gradient descent optimizers."); o.def("momentum", &optimizers::momentum, py::arg("beta")); o.def("demon", &optimizers::demon, py::arg("beta"), py::arg("max_ep")); diff --git a/src/utils.cpp b/src/utils.cpp @@ -18,10 +18,10 @@ #include <sys/stat.h> #include <Eigen/Dense> -// A bunch of hardcoded activation functions. Avoids much of the slowness of custom functions. -// Although the std::function makes it not the fastest way, the functionality is worth it. -// Yes, these functions may be a frustrating to read but they're just equations and I want to conserve space. +#include "utils.hpp" +namespace Jacobian { +namespace activations { float sigmoid(float x) {return 1.0/(1+exp(-x));} float sigmoid_deriv(float x) {return 1.0/(1+exp(-x)) * (1 - 1.0/(1+exp(-x)));} @@ -42,16 +42,16 @@ float cloglog_deriv(float x) {return exp(x-exp(x));} float step(float x) { - if (x > 0) return 1; - else return 0; + if (x > 0) return 1; + else return 0; } float step_deriv(float x) {return 0;} float bipolar(float x) { - if (x > 0) return 1; - else if (x == 0) return 0; - else return -1; + if (x > 0) return 1; + else if (x == 0) return 0; + else return -1; } float bipolar_deriv(float x) {return 0;} @@ -61,28 +61,77 @@ float bipolar_sigmoid_deriv(float x) {return (2*exp(x))/(pow(exp(x)+1,2));} float hard_tanh(float x) {return fmax(-1, fmin(1,x));} float hard_tanh_deriv(float x) { - if (-1 < x && x < 1) return 1; - else return 0; + if (-1 < x && x < 1) return 1; + else return 0; } float leaky_relu(float x) { - if (x > 0) return x; - else return 0.01 * x; + if (x > 0) return x; + else return 0.01 * x; } float leaky_relu_deriv(float x) { - if (x > 0) return 1; - else return 0.01; + if (x > 0) return 1; + else return 0.01; } std::function<float(float)> rectifier(float (*activation)(float)) { - auto rectified = [activation](float x) -> float - { - if (x > 0) return (*activation)(x); - else return 0; - }; - return rectified; + auto rectified = [activation](float x) -> float { + if (x > 0) + return (*activation)(x); + else + return 0; + }; + return rectified; +} +} // namespace activations + +namespace optimizers { +std::function<void(Layer&, Eigen::MatrixXf, float)> momentum(float beta) { + return [beta](Layer& layer, const Eigen::MatrixXf delta, const float learning_rate) { + layer.weights -= (beta * layer.m) + (learning_rate * delta); + layer.m = (learning_rate * delta); + }; +} + +std::function<void(Layer&, Eigen::MatrixXf, float)> demon(float beta, int max_ep) { + float beta_init = beta; + float prev_epoch = -1; + float epochs = 0; + return [max_ep, epochs, beta_init, beta](Layer& layer, const Eigen::MatrixXf delta, const float learning_rate) mutable { + beta = beta_init * (1-(epochs/max_ep)) / ((beta_init * (1-(epochs/max_ep))) + (1-beta_init)); + layer.weights -= (beta * layer.m) + (learning_rate * delta); + layer.m = (learning_rate * delta); + epochs++; + }; +} + +std::function<void(Layer&, Eigen::MatrixXf, float)> adam(float beta1, float beta2, float epsilon) { + return [beta1, beta2, epsilon](Layer& layer, const Eigen::MatrixXf delta, const float learning_rate) { + layer.m = (beta1 * layer.m) + ((1-beta1)*delta); + layer.v = (beta2 * layer.v) + (1-beta2)*(delta.cwiseProduct(delta)); + layer.weights -= learning_rate * + ((layer.v.cwiseSqrt()).array()+epsilon).pow(-1).cwiseProduct(layer.m.array()).matrix(); + }; +} + +std::function<void(Layer&, Eigen::MatrixXf, float)> adamax(float beta1, float beta2, float epsilon) { + return [beta1, beta2, epsilon](Layer &layer, + const Eigen::MatrixXf delta, + const float learning_rate) { + layer.m = (beta1 * layer.m) + ((1 - beta1) * delta); + if ((beta2 * layer.v).sum() > delta.array().abs().sum()) + layer.v = (beta2 * layer.v); + else + layer.v = delta.array().abs().matrix(); + layer.weights -= + learning_rate * + (layer.v.array().pow(-1).cwiseProduct(layer.m.array())) + .matrix(); + }; +} +} } diff --git a/src/utils.hpp b/src/utils.hpp @@ -9,9 +9,11 @@ #define UTILS_H #include <functional> -#include <immintrin.h> +#include <Eigen/Dense> +#include "bpnn.hpp" -// A zoo of activation functions. +namespace Jacobian { +namespace activations { float sigmoid(float x); float sigmoid_deriv(float x); float linear(float x); @@ -37,8 +39,24 @@ float bipolar_sigmoid_deriv(float x); float leaky_relu(float x); float leaky_relu_deriv(float x); std::function<float(float)> rectifier(float (*activation)(float)); +} + +namespace optimizers { +std::function<void(Layer&, Eigen::MatrixXf, float)> momentum(float beta); +std::function<void(Layer&, Eigen::MatrixXf, float)> demon(float beta_init, int max_ep); +std::function<void(Layer&, Eigen::MatrixXf, float)> adam(float beta1, float beta2, float epsilon); +std::function<void(Layer&, Eigen::MatrixXf, float)> adamax(float beta1, float beta2, float epsilon); +} + +namespace decays { +std::function<void(float &)> step(float a_0, float k); +std::function<void(float &)> exponential(float a_0, float k); +std::function<void(float &)> fractional(float a_0, float k); +std::function<void(float &)> linear(int max_ep); +} Eigen::MatrixXf strassen_mul(Eigen::MatrixXf a, Eigen::MatrixXf b); -#endif /* MODULE_H */ +} +#endif /* MODULE_H */