path: root/src/lib/src/neuralnet.c

#include <neuralnet/neuralnet.h>

#include "activation.h"
#include "neuralnet_impl.h"
#include <neuralnet/matrix.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;
}