jacobian

a basic keras-like neural network library for c++/python
Log | Files | Refs | README

bpnn.cpp (11814B)

      1 #include "bpnn.hpp"
      2 #include "utils.hpp"
      3 #include <random>
      4 
      5 namespace Jacobian {
      6 Layer::Layer(int batch_sz, int nodes)
      7 {
      8     contents = Eigen::MatrixXf(batch_sz, nodes);
      9     dZ = Eigen::MatrixXf(batch_sz, nodes);
     10     int datalen = batch_sz*nodes;
     11     for (int i = 0; i < datalen; i++) {
     12         contents(static_cast<int>(i / nodes),i%nodes) = 0;
     13         dZ(static_cast<int>(i / nodes),i%nodes) = 0;
     14     }
     15     bias = Eigen::MatrixXf (batch_sz, nodes);
     16     for (int i = 0; i < nodes; i++) {
     17         for (int j = 0; j < batch_sz; j++) bias(j, i) = 0;
     18     }
     19 }
     20 
     21 
     22 void Layer::init_weights(Layer next)
     23 {
     24     v = Eigen::MatrixXf (contents.cols(), next.contents.cols());
     25     m =  Eigen::MatrixXf (contents.cols(), next.contents.cols());
     26     weights = Eigen::MatrixXf (contents.cols(), next.contents.cols());
     27     int nodes = weights.cols();
     28     int n = contents.cols() + next.contents.cols();
     29     std::normal_distribution<float> d(0,sqrt(1.0/n));
     30     for (int i = 0; i < (weights.rows()*weights.cols()); i++) {
     31         std::random_device rd;
     32         std::mt19937 gen(rd());
     33         weights(static_cast<int>(i / nodes), i%nodes) = d(gen);
     34         v(static_cast<int>(i / nodes), i%nodes) = 0;
     35         m(static_cast<int>(i / nodes), i%nodes) = 0;
     36     }
     37 }
     38 
     39 Network::Network(const char* path, int batch_sz, float learn_rate, float bias_rate,
     40                  Regularization regularization, float l, float ratio, bool early_exit, float cutoff)
     41     :batch_size(batch_sz), learning_rate(learn_rate), bias_lr(bias_rate), reg_type(regularization),
     42      lambda(l), early_stop(early_exit), threshold(cutoff)
     43 {
     44     Expects(batch_size > 0 && learning_rate > 0 &&
     45             bias_rate > 0 && l >= 0 && ratio >= 0 && ratio <= 1);
     46     int total_instances = prep_file(path, SHUFFLED_PATH);
     47     val_instances = split_file(SHUFFLED_PATH, total_instances, ratio);
     48     prep(TRAIN_PATH, TRAIN_BIN_PATH);
     49     prep(VAL_PATH, VAL_BIN_PATH);
     50     data = open(TRAIN_BIN_PATH, O_RDONLY | O_NONBLOCK);
     51     val_data = open(VAL_BIN_PATH, O_RDONLY | O_NONBLOCK);
     52     instances = total_instances - val_instances;
     53     decay = [](float& learning_rate) -> void {};
     54     update = [](Layer& layer, const Eigen::MatrixXf delta, const float learning_rate) {
     55         layer.weights = (learning_rate * delta);
     56     };
     57     // File descriptors are nonnegative integers and open() returns -1 on failure.
     58     Ensures(batch_size < instances && data > 0 && val_data > 0);
     59 }
     60 
     61 Network::~Network()
     62 {
     63     close(data);
     64     close(val_data);
     65 }
     66 
     67 void Network::add_layer(int nodes, std::function<float(float)> activation,
     68                         std::function<float(float)> activation_deriv)
     69 {
     70     Expects(nodes > 0);
     71     length++;
     72     layers.emplace_back(batch_size, nodes);
     73     layers[length-1].activation = activation;
     74     layers[length-1].activation_deriv = activation_deriv;
     75 }
     76 
     77 void Network::initialize()
     78 {
     79     Expects(length > 1);
     80     labels = new Eigen::MatrixXf (batch_size,layers[length-1].contents.cols());
     81     for (int i = 0; i < length-1; i++) layers[i].init_weights(layers[i+1]);
     82 }
     83 
     84 void Network::set_activation(int index, std::function<float(float)> custom,
     85                              std::function<float(float)> custom_deriv)
     86 {
     87     Expects(index >= 0 && index < length);
     88     layers[index].activation = custom;
     89     layers[index].activation_deriv = custom_deriv;
     90 }
     91 
     92 void Network::softmax()
     93 {
     94     for (int i = 0; i < layers[length-1].contents.rows(); i++) {
     95         Eigen::MatrixXf m = layers[length-1].contents.block(i,0,1,layers[length-1].contents.cols());
     96         Eigen::MatrixXf::Index maxRow, maxCol;
     97         float max = m.maxCoeff(&maxRow, &maxCol);
     98         m = (m.array() - max).matrix();
     99         float sum = 0;
    100         for (int j = 0; j < layers[length-1].contents.cols(); j++) {
    101             sum += exp(m(0,j));
    102         }
    103         for (int j = 0; j < layers[length-1].contents.cols(); j++) {
    104             m(0,j) = exp(m(0,j))/sum;
    105         }
    106         layers[length-1].contents.block(i,0,1,layers[length-1].contents.cols()) = m;
    107     }
    108 }
    109 
    110 void Network::feedforward()
    111 {
    112     for (int i = 0; i < length-1; i++) {
    113         for (int j = 0; j < layers[i].contents.rows(); j++) {
    114             for (int k = 0; k < layers[i].contents.cols(); k++) {
    115                 layers[i].dZ(j,k) = layers[i].activation_deriv(layers[i].contents(j,k));
    116                 layers[i].contents(j,k) = layers[i].activation(layers[i].contents(j,k));
    117             }
    118         }
    119         layers[i+1].contents = layers[i].contents * layers[i].weights;
    120         layers[i+1].contents += layers[i+1].bias;
    121     }
    122     for (int j = 0; j < layers[length-1].contents.rows(); j++) {
    123         for (int k = 0; k < layers[length-1].contents.cols(); k++) {
    124             layers[length-1].dZ(j,k) = layers[length-1].activation_deriv(layers[length-1].contents(j,k));
    125             layers[length-1].contents(j,k) = layers[length-1].activation(layers[length-1].contents(j,k));
    126         }
    127     }
    128     softmax();
    129 }
    130 
    131 std::function<void(float&)> decays::step(float a_0, float k)
    132 {
    133     return [a_0, k](float& learning_rate) -> void {
    134         learning_rate = a_0 * learning_rate/k;
    135     };
    136 }
    137 
    138 std::function<void(float&)> decays::exponential(float a_0, float k)
    139 {
    140     int epochs = 0;
    141     return [a_0, k, epochs](float& learning_rate) mutable -> void {
    142         learning_rate = a_0 * exp(-k * epochs);
    143         epochs++;
    144     };
    145 }
    146 
    147 std::function<void(float&)> decays::fractional(float a_0, float k)
    148 {
    149     int epochs = 0;
    150     return [a_0, k, epochs](float& learning_rate) mutable -> void {
    151         learning_rate = a_0 / (1+(k * epochs));
    152         epochs++;
    153     };
    154 }
    155 
    156 std::function<void(float&)> decays::linear(int max_ep)
    157 {
    158     int epochs = 0;
    159     return [max_ep, epochs](float& learning_rate) mutable -> void {
    160         learning_rate = 1 - epochs/max_ep;
    161         epochs++;
    162     };
    163 }
    164 
    165 
    166 void Network::init_decay(std::function<void(float&)> f)
    167 {
    168     decay = f;
    169 }
    170 
    171 void Network::list_net()
    172 {
    173     Expects(length > 1);
    174     std::cout << "-----------------------\nINPUT LAYER (LAYER 0)\n"
    175               <<  "\n\n\u001b[31mACTIVATIONS:\x1B[0;37m\n" << layers[0].contents
    176               << "\n\n\u001b[31mWEIGHTS:\x1B[0;37m\n" << layers[0].weights
    177               << "\n\n\u001b[31mBIASES:\x1B[0;37m\n" << layers[0].bias << "\n\n\n";
    178     for (int i = 1; i < length-1; i++) {
    179         std::cout << "-----------------------\nLAYER " << i
    180                   << "\n\n\u001b[31mACTIVATIONS:\x1B[0;37m\n" << layers[i].contents
    181                   << "\n\n\u001b[31mBIASES:\x1B[0;37m\n" << layers[i].bias
    182                   << "\n\n\u001b[31mWEIGHTS:\x1B[0;37m\n" << layers[i].weights << "\n\n\n";
    183     }
    184     std::cout << "-----------------------\nOUTPUT LAYER (LAYER " << length-1
    185               <<"\n\n\u001b[31mACTIVATIONS:\x1B[0;37m\n" << layers[length-1].contents
    186               << "\n\n\u001b[31BIASES:\x1B[0;37m\n" << layers[length-1].bias <<  "\n\n\n";
    187 }
    188 
    189 float Network::cost()
    190 {
    191     float sum = 0;
    192     float reg = 0; // Regularization term
    193     for (int i = 0; i < layers[length-1].contents.rows(); i++) {
    194         float tempsum = 0;
    195         for (int j = 0; j < layers[length-1].contents.cols(); j++) {
    196             float truth;
    197             if (j==(*labels)(i,0)) truth = 1;
    198             else truth = 0;
    199             if (layers[length-1].contents(i,j) == 0) layers[length-1].contents(i,j) += 0.00001;
    200             tempsum += truth * log(layers[length-1].contents(i,j));
    201         }
    202         sum-=tempsum;
    203     }
    204     for (unsigned long i = 0; i < layers.size()-1; i++) {
    205         if (reg_type == Regularization::L2) {
    206             reg += layers[i].weights.cwiseProduct(layers[i].weights).sum();
    207         } else if (reg_type == Regularization::L1) {
    208             reg += (layers[i].weights.array().abs().matrix()).sum();
    209         }
    210     }
    211     return ((1.0/batch_size) * sum) + (1/2*lambda*reg);
    212 }
    213 
    214 float Network::accuracy()
    215 {
    216     float correct = 0;
    217     for (int i = 0; i < layers[length-1].contents.rows(); i++) {
    218         float ans = -INFINITY;
    219         float index = -1;
    220         for (int j = 0; j < layers[length-1].contents.cols(); j++) {
    221             if (layers[length-1].contents(i, j) > ans) {
    222                 ans = layers[length-1].contents(i, j);
    223                 index = j;
    224             }
    225         }
    226         if ((*labels)(i, 0) == index) correct += 1;
    227     }
    228     return (1.0/batch_size) * correct;
    229 }
    230 
    231 Eigen::MatrixXf l1_deriv(Eigen::MatrixXf m)
    232 {
    233     Eigen::MatrixXf r(m.rows(), m.cols());
    234     for (int i = 0; i < m.rows(); i++) {
    235         for (int j = 0; j < m.cols(); j++) {
    236             if (m(i,j) == 0) r(i,j) = 0;
    237             else r(i,j) = 1;
    238         }
    239     }
    240     return r;
    241 }
    242 
    243 Eigen::MatrixXf Network::backpropagate()
    244 {
    245     std::vector<Eigen::MatrixXf> gradients;
    246     std::vector<Eigen::MatrixXf> deltas;
    247     Eigen::MatrixXf error (layers[length-1].contents.rows(), layers[length-1].contents.cols());
    248     for (int i = 0; i < error.rows(); i++) {
    249         for (int j = 0; j < error.cols(); j++) {
    250             float truth;
    251             if (j==(*labels)(i,0)) truth = 1;
    252             else truth = 0;
    253             error(i,j) = layers[length-1].contents(i,j) - truth;
    254         }
    255     }
    256     gradients.push_back(error);
    257     deltas.push_back(layers[length-2].contents.transpose() * gradients[0]);
    258     int counter = 1;
    259     for (int i = length-2; i >= 1; i--) {
    260         // TODO: Add nesterov momentum
    261         // (*layers[i].weights-((learning_rate * *layers[i].weights) + (0.9 * *layers[i].v))).transpose()
    262         //grad_calc(gradients, counter, i)
    263         gradients.push_back((gradients[counter-1] * layers[i].weights.transpose()).cwiseProduct(layers[i].dZ));
    264         deltas.push_back(layers[i-1].contents.transpose() * gradients[counter]);
    265         counter++;
    266     }
    267     for (int i = 0; i < length-1; i++) {
    268         update(layers[length-2-i], deltas[i], learning_rate);
    269         if (reg_type == Regularization::L2) {
    270           layers[length-2-i].weights -= ((lambda/batch_size) * (layers[length-2-i].weights));
    271         } else if (reg_type == Regularization::L1) {
    272           layers[length-2-i].weights -= ((lambda/(2*batch_size)) * l1_deriv(layers[length-2-i].weights));
    273         }
    274         layers[length-1-i].bias -= bias_lr * gradients[i];
    275     }
    276     return gradients.back();
    277 }
    278 
    279 #include "data.cpp"
    280 
    281 void Network::validate(const char* path)
    282 {
    283     if (val_instances == 0) return;
    284     float costsum = 0;
    285     float accsum = 0;
    286     for (int i = 0; i <= val_instances-batch_size; i+=batch_size) {
    287         next_batch(val_data);
    288         feedforward();
    289         costsum += cost();
    290         accsum += accuracy();
    291     }
    292     val_acc = 1.0/(static_cast<float>(val_instances/batch_size)) * accsum;
    293     val_cost = 1.0/(static_cast<float>(val_instances/batch_size)) * costsum;
    294     val_data = open(VAL_BIN_PATH, O_RDONLY | O_NONBLOCK);
    295     Ensures(lseek(val_data, 0, SEEK_CUR) == 0);
    296 }
    297 
    298 void Network::interactive_next_batch()
    299 {
    300     if (batches < instances/batch_size-batch_size) next_batch(data);
    301     else {
    302         batches = 0;
    303         data = open(TRAIN_BIN_PATH, O_RDONLY | O_NONBLOCK);
    304         decay(learning_rate);
    305     }
    306     batches++;
    307 }
    308 
    309 void Network::train()
    310 {
    311     float cost_sum = 0;
    312     float acc_sum = 0;
    313     for (int i = 0; i <= instances - batch_size; i += batch_size) {
    314         if (i != instances - batch_size)
    315             next_batch(data);
    316         feedforward();
    317         backpropagate();
    318         cost_sum += cost();
    319         acc_sum += accuracy();
    320         batches++;
    321     }
    322     epoch_acc =
    323         1.0 / (static_cast<float>(instances / batch_size)) * acc_sum;
    324     epoch_cost =
    325         1.0 / (static_cast<float>(instances / batch_size)) * cost_sum;
    326     validate(VAL_PATH);
    327     if (silenced == false)
    328         printf("Epoch %i complete - cost %f - acc %f - val_cost %f - val_acc %f\n",
    329                epochs, epoch_cost, epoch_acc, val_cost, val_acc);
    330     batches = 1;
    331     data = open(TRAIN_BIN_PATH, O_RDONLY | O_NONBLOCK);
    332     decay(learning_rate);
    333     epochs++;
    334     Ensures(lseek(data, 0, SEEK_CUR) == 0);
    335 }
    336 }