1 files changed, 228 insertions, 0 deletions

diff --git a/src/lib/src/neuralnet.c b/src/lib/src/neuralnet.c
new file mode 100644
index 0000000..cac611a
--- /dev/null
+++ b/src/lib/src/neuralnet.c
@@ -0,0 +1,228 @@
+#include <neuralnet/neuralnet.h>
+
+#include <neuralnet/matrix.h>
+#include "activation.h"
+#include "neuralnet_impl.h"
+
+#include <assert.h>
+#include <stdlib.h>
+
+nnNeuralNetwork* nnMakeNet(int num_layers, const int* layer_sizes, const nnActivation* activations) {
+  assert(num_layers > 0);
+  assert(layer_sizes);
+  assert(activations);
+
+  nnNeuralNetwork* net = calloc(1, sizeof(nnNeuralNetwork));
+  if (net == 0) {
+    return 0;
+  }
+
+  net->num_layers = num_layers;
+
+  net->weights = calloc(num_layers, sizeof(nnMatrix));
+  net->biases  = calloc(num_layers, sizeof(nnMatrix));
+  net->activations = calloc(num_layers, sizeof(nnActivation));
+  if ( (net->weights == 0) || (net->biases == 0) || (net->activations == 0) ) {
+    nnDeleteNet(&net);
+    return 0;
+  }
+
+  for (int l = 0; l < num_layers; ++l) {
+    // layer_sizes = { input layer size, first hidden layer size, ...}
+    const int layer_input_size  = layer_sizes[l];
+    const int layer_output_size = layer_sizes[l+1];
+
+    // We store the transpose of the weight matrix as written in textbooks.
+    // Our vectors are row vectors and the matrices row-major.
+    const int rows = layer_input_size;
+    const int cols = layer_output_size;
+
+    net->weights[l] = nnMatrixMake(rows, cols);
+    net->biases[l]  = nnMatrixMake(1, cols);
+    net->activations[l] = activations[l];
+  }
+
+  return net;
+}
+
+void nnDeleteNet(nnNeuralNetwork** net) {
+  if ( (!net) || (!(*net)) ) {
+    return;
+  }
+  if ((*net)->weights != 0) {
+    for (int l = 0; l < (*net)->num_layers; ++l) {
+      nnMatrixDel(&(*net)->weights[l]);
+    }
+    free((*net)->weights);
+    (*net)->weights = 0;
+  }
+  if ((*net)->biases != 0) {
+    for (int l = 0; l < (*net)->num_layers; ++l) {
+      nnMatrixDel(&(*net)->biases[l]);
+    }
+    free((*net)->biases);
+    (*net)->biases = 0;
+  }
+  if ((*net)->activations) {
+    free((*net)->activations);
+    (*net)->activations = 0;
+  }
+  free(*net);
+  *net = 0;
+}
+
+void nnSetWeights(nnNeuralNetwork* net, const R* weights) {
+  assert(net);
+  assert(weights);
+
+  for (int l = 0; l < net->num_layers; ++l) {
+    nnMatrix* layer_weights = &net->weights[l];
+    R* layer_values = layer_weights->values;
+
+    for (int j = 0; j < layer_weights->rows * layer_weights->cols; ++j) {
+      *layer_values++ = *weights++;
+    }
+  }
+}
+
+void nnSetBiases(nnNeuralNetwork* net, const R* biases) {
+  assert(net);
+  assert(biases);
+
+  for (int l = 0; l < net->num_layers; ++l) {
+    nnMatrix* layer_biases = &net->biases[l];
+    R* layer_values = layer_biases->values;
+
+    for (int j = 0; j < layer_biases->rows * layer_biases->cols; ++j) {
+      *layer_values++ = *biases++;
+    }
+  }
+}
+
+void nnQuery(const nnNeuralNetwork* net, nnQueryObject* query, const nnMatrix* input) {
+  assert(net);
+  assert(query);
+  assert(input);
+  assert(net->num_layers == query->num_layers);
+  assert(input->rows <= query->network_outputs->rows);
+  assert(input->cols == nnNetInputSize(net));
+
+  for (int i = 0; i < input->rows; ++i) {
+    // Not mutating the input, but we need the cast to borrow.
+    nnMatrix input_vector = nnMatrixBorrowRows((nnMatrix*)input, i, 1);
+
+    for (int l = 0; l < net->num_layers; ++l) {
+      const nnMatrix* layer_weights = &net->weights[l];
+      const nnMatrix* layer_biases  = &net->biases[l];
+      // Y^T = (W*X)^T = X^T*W^T
+      //
+      // TODO: If we had a row-row matrix multiplication, we could compute:
+      //   Y^T = W ** X^T
+      // The row-row multiplication could be more cache-friendly. We just need
+      // to store W as is, without transposing.
+      // We could also rewrite the original Mul function to go row x row,
+      // decomposing the multiplication. Preserving the original meaning of Mul
+      // makes everything clearer.
+      nnMatrix output_vector = nnMatrixBorrowRows(&query->layer_outputs[l], i, 1);
+      nnMatrixMul(&input_vector, layer_weights, &output_vector);
+      nnMatrixAddRow(&output_vector, layer_biases, &output_vector);
+
+      switch (net->activations[l]) {
+        case nnIdentity:
+          break;  // Nothing to do for the identity function.
+        case nnSigmoid:
+          sigmoid_array(output_vector.values, output_vector.values, output_vector.cols);
+          break;
+        case nnRelu:
+          relu_array(output_vector.values, output_vector.values, output_vector.cols);
+          break;
+        default:
+          assert(0);
+      }
+
+      input_vector = output_vector;  // Borrow.
+    }
+  }
+}
+
+void nnQueryArray(const nnNeuralNetwork* net, nnQueryObject* query, const R* input, R* output) {
+  assert(net);
+  assert(query);
+  assert(input);
+  assert(output);
+  assert(net->num_layers > 0);
+
+  nnMatrix input_vector = nnMatrixMake(net->weights[0].cols, 1);
+  nnMatrixInit(&input_vector, input);
+  nnQuery(net, query, &input_vector);
+  nnMatrixRowToArray(query->network_outputs, 0, output);
+}
+
+nnQueryObject* nnMakeQueryObject(const nnNeuralNetwork* net, int num_inputs) {
+  assert(net);
+  assert(num_inputs > 0);
+  assert(net->num_layers > 0);
+
+  nnQueryObject* query = calloc(1, sizeof(nnQueryObject));
+  if (!query) {
+    return 0;
+  }
+
+  query->num_layers = net->num_layers;
+
+  // Allocate the intermediate layer output matrices.
+  query->layer_outputs = calloc(net->num_layers, sizeof(nnMatrix));
+  if (!query->layer_outputs) {
+    free(query);
+    return 0;
+  }
+  for (int l = 0; l < net->num_layers; ++l) {
+      const nnMatrix* layer_weights = &net->weights[l];
+      const int layer_output_size = nnLayerOutputSize(layer_weights);
+      query->layer_outputs[l] = nnMatrixMake(num_inputs, layer_output_size);
+  }
+  query->network_outputs = &query->layer_outputs[net->num_layers - 1];
+
+  return query;
+}
+
+void nnDeleteQueryObject(nnQueryObject** query) {
+  if ( (!query) || (!(*query)) ) {
+    return;
+  }
+  if ((*query)->layer_outputs != 0) {
+    for (int l = 0; l < (*query)->num_layers; ++l) {
+      nnMatrixDel(&(*query)->layer_outputs[l]);
+    }
+  }
+  free((*query)->layer_outputs);
+  free(*query);
+  *query = 0;
+}
+
+const nnMatrix* nnNetOutputs(const nnQueryObject* query) {
+  assert(query);
+  return query->network_outputs;
+}
+
+int nnNetInputSize(const nnNeuralNetwork* net) {
+  assert(net);
+  assert(net->num_layers > 0);
+  return net->weights[0].rows;
+}
+
+int nnNetOutputSize(const nnNeuralNetwork* net) {
+  assert(net);
+  assert(net->num_layers > 0);
+  return net->weights[net->num_layers - 1].cols;
+}
+
+int nnLayerInputSize(const nnMatrix* weights) {
+  assert(weights);
+  return weights->rows;
+}
+
+int nnLayerOutputSize(const nnMatrix* weights) {
+  assert(weights);
+  return weights->cols;
+}