Simplest ANN implementation in pure C without libraries

lib_math_vectors.c

#include "lib_math_vectors.h"

// Simplifies math power function
double pow2(double x) { return pow(x, 2); }

// Sigmoidal function
double sigmoid(double x) { return 1.0 / (1.0 + exp(-x)); }

// ReLU function
double relu(double x) { return fmax(0.0, x); }

// Signum function
double sign(double x) {
    if (x > 0.0) { return 1.0; }
    if (x < 0.0) { return -1.0; }
    return 0.0;
}

// Allocates zero matrix of specified num of rows
matrix_t* matrix_init(uint16_t num)
{
    matrix_t* mat = (matrix_t*) malloc(sizeof(matrix_t));
    mat->num_rows = num;
    mat->row_vectors = (vector_t**) malloc(num * sizeof(vector_t*));
    return mat;
}

// Initializes matrix and fills it with specified vectors
matrix_t* matrix_create(uint16_t num_rows, ...)
{
    matrix_t* mat = matrix_init(num_rows);
    va_list valist;
    va_start(valist, num_rows);
    for(uint16_t i = 0; i < num_rows; i++) {
        mat->row_vectors[i] = va_arg(valist, vector_t*);
        assert((i == 0 || mat->row_vectors[0]->len == mat->row_vectors[i]->len) && "Vector appended to matrix has incorrect size!");
    }
    va_end(valist);
    return mat;
}

// Initializes matrix of specified number of rows and cols, can be filled by const or randomly
matrix_t* matrix_create_sized(uint16_t num_rows, uint16_t num_cols, bool random, double c)
{
    matrix_t* mat = matrix_init(num_rows);
    for(uint16_t i = 0; i < num_rows; i++) {
        mat->row_vectors[i] = vector_init(num_cols, 0);
        if (c != 0 || random) {
            for (uint16_t j = 0; j < num_cols; j++) {
                if (random) {
                    mat->row_vectors[i]->values[j] = ((double) rand()) / RAND_MAX;
                } else {
                    mat->row_vectors[i]->values[j] = c;
                }
            }
        }
    }
    return mat;
}

// Frees space allocated by matrix
void matrix_destroy(matrix_t* m)
{
    if (m == NULL) {
        return;
    }
    for(int i = 0; i < m->num_rows; i++) {
        vector_destroy(m->row_vectors[i]);
    }
    free(m->row_vectors);
    free(m);
}

// Copies matrix
matrix_t* matrix_copy(matrix_t* m)
{
    matrix_t* res = matrix_init(m->num_rows);
    for (uint16_t i = 0; i < res->num_rows; i++) {
        res->row_vectors[i] = vector_copy(m->row_vectors[i]);
    }
    return res;
}

// Adds vector to each row of matrix
matrix_t* matrix_add_vector_rowwise(matrix_t* m, bool inplace, bool substract, vector_t* v)
{
    matrix_t* res;
    if (inplace) {
        res = m;
    } else {
        res = matrix_copy(m);
    }
    for (int i = 0; i < res->num_rows; i++) {
        vectors_add(true, substract, 2, res->row_vectors[i], v);
    }
    return res;
}

// Adds matrices together
matrix_t* matrices_add(bool to_first_inplace, bool substract, uint16_t num, ...)
{
    matrix_t* res;
    va_list valist;
    va_start(valist, num);
    matrix_t* mat_0 = va_arg(valist, matrix_t*);
    if (to_first_inplace) {
        res = mat_0;
    } else {
        res = matrix_copy(mat_0);
    }
    matrix_t* mat_temp;
    for(uint16_t i = 1; i < num; i++) {
        mat_temp = va_arg(valist, matrix_t*);
        assert(mat_0->num_rows == mat_temp->num_rows && mat_0->row_vectors[0]->len == mat_temp->row_vectors[0]->len);
        for (uint16_t i = 0; i < res->num_rows; i++) {
            matrix_add_vector_rowwise(res, true, substract, mat_temp->row_vectors[i]);
        }
    }
    va_end(valist);
    return res;
}

// Multiplies vector by matrix from left, returns new vector
vector_t* matrix_dot_vector(matrix_t* m, vector_t* v)
{
    assert(m->row_vectors[0]->len == v->len);
    vector_t* temp = vector_copy(v);
    vector_t* res = vector_init(m->num_rows, 0);
    for (uint16_t i = 0; i < m->num_rows; i++) {
        res->values[i] = vector_mult_vector(m->row_vectors[i], temp);
    }
    vector_destroy(temp);
    return res;
}

// Multiplies matrix by matrix, returns new matrix
matrix_t* matrix_dot_matrix(matrix_t* m1, matrix_t* m2)
{
    assert(m1->row_vectors[0]->len == m2->num_rows);
    matrix_t* m2T = matrix_transpose(m2);
    matrix_t* resT = matrix_init(m2T->num_rows);
    for (int i = 0; i < m2T->num_rows; i++) {
        resT->row_vectors[i] = matrix_dot_vector(m1, m2T->row_vectors[i]);
    }
    matrix_t* res = matrix_transpose(resT);
    matrix_destroy(m2T);
    matrix_destroy(resT);
    return res;
}

// Creates transponed matrix
matrix_t* matrix_transpose(matrix_t* m)
{
    uint16_t width = m->row_vectors[0]->len;
    matrix_t* res = matrix_init(width);
    for (uint16_t i = 0; i < res->num_rows; i++) {
        res->row_vectors[i] = vector_init(m->num_rows, 0);
    }
    for (uint16_t i = 0; i < m->num_rows; i++) {
        for (uint16_t j = 0; j < width; j++) {
            res->row_vectors[j]->values[i] = m->row_vectors[i]->values[j];
        }
    }
    return res;
}

// Applies function to each element of provided matrix
matrix_t* matrix_elwise_apply(matrix_t* m, bool inplace, double (*f)(double))
{
    matrix_t* res;
    if (inplace) {
        res = m;
    } else {
        res = matrix_copy(m);
    }
    for (uint16_t i = 0; i < m->num_rows; i++) {
        vector_elwise_apply(res->row_vectors[i], true, f);
    }
    return res;
}

// Mutliplies matrix by matrix element-wise
matrix_t* matrix_mult_matrix(matrix_t* m, bool inplace, matrix_t* m2)
{
    assert(m->num_rows == m2->num_rows && m->row_vectors[0]->len == m2->row_vectors[0]->len);
    matrix_t* res;
    if (inplace) {
        res = m;
    } else {
        res = matrix_copy(m);
    }
    for (uint16_t i = 0; i < m->num_rows; i++) {
        for (uint16_t j = 0; j < m->row_vectors[0]->len; j++) {
            res->row_vectors[i]->values[j] *= m2->row_vectors[i]->values[j];
        }
    }
    return res;
}

// Performs softmax function on each row vector of provided matrix
matrix_t* matrix_row_softmax(matrix_t* m, bool inplace) 
{
    matrix_t* res;
    if (inplace) {
        res = m;
    } else {
        res = matrix_copy(m);
    }
    for (uint16_t i = 0; i < m->num_rows; i++) {
        vector_softmax(res->row_vectors[i], true);
    }
    return res;
}

// Sums matrix over specified axis: 0 = rows (top to bottom), 1 = columns (left to right)
vector_t* matrix_sum(matrix_t* m, uint8_t axis) 
{
    vector_t* res = NULL;
    if (axis == 0) {
        matrix_t* mT = matrix_transpose(m);
        vector_t* ones = vector_init(m->num_rows, 1.0);
        res = matrix_dot_vector(mT, ones);
        vector_destroy(ones);
        matrix_destroy(mT);
    } else if (axis == 1) {
        res = vector_init(m->num_rows, 0);
        for (int i = 0; i < m->num_rows; i++) {
            res->values[i] = vector_sum(m->row_vectors[i]);
        }
    }
    return res;
}

// Multiplies matrix values by specified constant number
matrix_t* matrix_mult_const(matrix_t* m, bool inplace, double c)
{
    matrix_t* res;
    if (inplace) {
        res = m;
    } else {
        res = matrix_copy(m);
    }
    for (uint16_t i = 0; i < m->num_rows; i++) {
        vector_mult_const(m->row_vectors[i], true, c);
    }
    return res;
}

// Prints out matrix
void matrix_print(matrix_t* m)
{
    uint16_t width = m->row_vectors[0]->len;
    printf("Matrix of %d x %d: [ \n", m->num_rows, width);
    for(int i = 0; i < m->num_rows; i++) {
        printf("\t[ ");
        for(int j = 0; j < width; j++) {
            printf("%f ", m->row_vectors[i]->values[j]);
        }
        printf("]\n");
    }
    printf("]\n");
}

// Allocates zero vector of specified length
vector_t* vector_init(uint16_t len, double init_value)
{
    vector_t* vec = (vector_t*) malloc(sizeof(vector_t));
    vec->len = len;
    vec->values = (double*) calloc(len, sizeof(double));
    if (init_value != 0) {
        for (int i = 0; i < len; i++) {
            vec->values[i] = init_value;
        }
    }
    return vec;
}

// Initializes vector and fills it with specified values
vector_t* vector_create(uint16_t len, ...)
{
    vector_t* vec = vector_init(len, 0);
    va_list valist;
    va_start(valist, len);
    for(uint16_t i = 0; i < len; i++) {
        vec->values[i] = va_arg(valist, double);
    }
    va_end(valist);
    return vec;
}

// Frees memory of vector
void vector_destroy(vector_t* vec)
{
    if (vec == NULL) {
        return;
    }
    free(vec->values);
    free(vec);
}

// Applies provided function on provided vector element-wise
vector_t* vector_elwise_apply(vector_t* v, bool inplace, double (*f)(double))
{
    vector_t* res;
    if (inplace) {
        res = v;
    } else {
        res = vector_init(v->len, 0);
    }
    for (uint16_t i = 0; i < res->len; i++) {
        res->values[i] = f(v->values[i]);
    }
    return res;
}

// Sums all values of provided vector
double vector_sum(vector_t* v)
{
    double sum = 0;
    for (uint16_t i = 0; i < v->len; i++) {
        sum += v->values[i];
    }
    return sum;
}

// Copies provided vector
vector_t* vector_copy(vector_t* v)
{
    vector_t* res = vector_init(v->len, 0);
    for (uint16_t i = 0; i < res->len; i++) {
        res->values[i] = v->values[i];
    }
    return res;
}

// Adds vectors together, optionally substracts from first vector
vector_t* vectors_add(bool to_first_inplace, bool substract, uint16_t num, ...)
{
    vector_t* res;
    va_list valist;
    va_start(valist, num);
    vector_t* vec_0 = va_arg(valist, vector_t*);
    if (to_first_inplace) {
        res = vec_0;
    } else {
        res = vector_copy(vec_0);
    }
    vector_t* vec_temp;
    for(uint16_t i = 1; i < num; i++) {
        vec_temp = va_arg(valist, vector_t*);
        assert(vec_temp->len == vec_0->len);
        for (uint16_t i = 0; i < res->len; i++) {
            if (substract) {
                res->values[i] -= vec_temp->values[i];
            } else {
                res->values[i] += vec_temp->values[i];
            }
        }
    }
    va_end(valist);
    return res;
}

// Multiplies vector values by specified constant number
vector_t* vector_mult_const(vector_t* v, bool inplace, double c)
{
    vector_t* res;
    if (inplace) {
        res = v;
    } else {
        res = vector_copy(v);
    }
    for (uint16_t i = 0; i < res->len; i++) {
        res->values[i] *= c;
    }
    return res;
}

// Multiplies vector by another vector, standard scalar multiplication
double vector_mult_vector(vector_t* v1, vector_t* v2) 
{
    assert(v1->len == v2->len);
    double res = 0;
    for (uint16_t i = 0; i < v1->len; i++) {
        res += v1->values[i] * v2->values[i];
    }
    assert("Result of vector multiplication is NaN!" && (res == res));
    return res;
}

// Calculates Eucleidian vector distance
double vector_length(vector_t* v)
{
    return sqrt(vector_mult_vector(v, v));
}

// Performs softmax function on provided vector
vector_t* vector_softmax(vector_t* v, bool inplace) 
{
    vector_t* res = vector_elwise_apply(v, inplace, &exp);
    vector_mult_const(res, true, 1./vector_sum(res));
    return res;
}

// Prints out vector (length and content)
void vector_print(vector_t* vec)
{
    printf("Vector of %d values: [ ", vec->len);
    for(int i = 0; i < vec->len; i++) {
        printf("%f ", vec->values[i]);
    }
    printf("]\n");
}

lib_math_vectors.h

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include <stdarg.h>
#include <stdbool.h>
#include <assert.h>

double pow2(double x);
double sigmoid(double x);
double relu(double x);
double sign(double x);

typedef struct {
    double* values;
    uint16_t len;
} vector_t;

vector_t* vector_init(uint16_t len, double init_value);
vector_t* vector_create(uint16_t len, ...);
void vector_destroy(vector_t* vec);
vector_t* vector_elwise_apply(vector_t* v, bool inplace, double (*f)(double));
double vector_sum(vector_t* v);
vector_t* vector_copy(vector_t* v);
vector_t* vectors_add(bool to_first_inplace, bool substract, uint16_t num, ...);
vector_t* vector_mult_const(vector_t* v, bool inplace, double c);
double vector_mult_vector(vector_t* v1, vector_t* v2);
vector_t* vector_softmax(vector_t* v, bool inplace);
double vector_length(vector_t* v);
void vector_print(vector_t* vec);

typedef struct {
    vector_t** row_vectors;
    uint16_t num_rows;
} matrix_t;

matrix_t* matrix_init(uint16_t num);
matrix_t* matrix_create(uint16_t num_rows, ...);
matrix_t* matrix_create_sized(uint16_t num_rows, uint16_t num_cols, bool random, double c);
void matrix_destroy(matrix_t* m);
matrix_t* matrix_copy(matrix_t* m);
vector_t* matrix_dot_vector(matrix_t* m, vector_t* v);
void matrix_print(matrix_t* m);
matrix_t* matrix_transpose(matrix_t* m);
matrix_t* matrix_dot_matrix(matrix_t* m1, matrix_t* m2);
matrix_t* matrix_add_vector_rowwise(matrix_t* m, bool inplace, bool substract, vector_t* v);
matrix_t* matrix_elwise_apply(matrix_t* m, bool inplace, double (*f)(double));
matrix_t* matrix_row_softmax(matrix_t* m, bool inplace);
matrix_t* matrices_add(bool to_first_inplace, bool substract, uint16_t num, ...);
vector_t* matrix_sum(matrix_t* m, uint8_t axis);
matrix_t* matrix_mult_const(matrix_t* m, bool inplace, double c);
matrix_t* matrix_mult_matrix(matrix_t* m, bool inplace, matrix_t* m2);

lib_nn_layers.c

#include "lib_nn_layers.h"

/* DEFINITIONS OF ACTIVATION FUNCTIONS */

// Derivative of tanh function for matrices: 1-f^2(x)
matrix_t* tanh_deriv(matrix_t* prev_prop)
{
    matrix_t* temp = matrix_elwise_apply(prev_prop, false, *pow2);
    matrix_t* temp2 = matrix_create_sized(prev_prop->num_rows, prev_prop->row_vectors[0]->len, false, 1.0);
    matrices_add(true, true, 2, temp2, temp);
    matrix_destroy(temp);
    return temp2;
}
// tanh activation function with derivative 1-tanh^2
const activation_function_t fa_tanh = { .calculate = *tanh, .calculate_derivative = *tanh_deriv };

// Derivative of sigmoidal function for matrices: f(x)*(1-f(x))
matrix_t* sigmoid_deriv(matrix_t* prev_prop)
{
    matrix_t* temp = matrix_create_sized(prev_prop->num_rows, prev_prop->row_vectors[0]->len, false, 1.0);
    matrices_add(true, true, 2, temp, prev_prop);
    matrix_mult_matrix(temp, true, prev_prop);
    return temp;
}
// sigmoid activation function with derivative sigmoid*(1-sigmoid)
const activation_function_t fa_sigmoid = { .calculate = *sigmoid, .calculate_derivative = *sigmoid_deriv };

// Derivative of ReLU function for matrices: sign(f(x))
matrix_t* relu_deriv(matrix_t* prev_prop)
{
    return matrix_elwise_apply(prev_prop, false, *sign);;
}
// ReLU activation function with derivative
const activation_function_t fa_relu = { .calculate = *relu, .calculate_derivative = *relu_deriv };

/* NN functions definitions */
// TODO put NN functions here

lib_nn_layers.h

#include "lib_math_vectors.h"

// Lambda functions implementation for C language (works only with GCC)
#define __LAMBDA__(return_type, body_and_args) ({ return_type __fn__ body_and_args __fn__; })

typedef struct {
    double (*calculate)(double);
    matrix_t* (*calculate_derivative)(matrix_t*);
} activation_function_t;

typedef struct {
    matrix_t* weights_matrix;
    vector_t* b_vector;
    activation_function_t act_fun;
} nn_layer_t;

typedef struct {
    nn_layer_t** layers;
    uint16_t num_layers;
    double epsilon;
    double reg_lambda;
    activation_function_t act_fun;
} neural_network_t;

extern const activation_function_t fa_tanh;
extern const activation_function_t fa_sigmoid;
extern const activation_function_t fa_relu;

nn_multilayer.c

// Compile: gcc lib_math_vectors.c lib_nn_layers.c nn_multilayer.c -g -lm -o nn_multilayer.exe

#include "lib_nn_layers.h"


// Creates NN layer with specified input dimension, number of neurons and activation function
nn_layer_t* nn_layer_create(int in_dim, int num_neurons, activation_function_t act_fun)
{
    nn_layer_t* l = (nn_layer_t*) malloc(sizeof(nn_layer_t));
    l->weights_matrix = matrix_create_sized(num_neurons, in_dim, true, 0);
    l->b_vector = vector_init(num_neurons, 0);
    l->act_fun = act_fun;
    return l;
}

// Frees space allocated by NN layer
void nn_layer_destroy(nn_layer_t* l)
{
    matrix_destroy(l->weights_matrix);
    vector_destroy(l->b_vector);
    free(l);
}

// Performs forward propagation with NN layer on input vector
vector_t* nn_layer_forward_prop(nn_layer_t* l, vector_t* v)
{
    vector_t* res = matrix_dot_vector(l->weights_matrix, v);
    vectors_add(true, false, 2, res, l->b_vector);
    vector_elwise_apply(res, true, l->act_fun.calculate);
    return res;
}

// Performs backward propagation with NN layer, returns new delta matrix for previous layer
matrix_t* nn_layer_backward_prop(nn_layer_t* l, matrix_t* delta, double reg_lambda, double epsilon, matrix_t* prev_prop)
{
    matrix_t* temp;
    
    // Calculate new weight matrix from delta
    temp = matrix_transpose(delta);
    matrix_t* dW = matrix_dot_matrix(temp, prev_prop);
    matrix_destroy(temp);
    matrix_mult_const(dW, true, epsilon);
    // + Apply regularization to weight matrix
    temp = matrix_mult_const(l->weights_matrix, false, reg_lambda);
    matrices_add(true, false, 2, dW, temp);
    matrix_destroy(temp);
    matrices_add(true, true, 2, l->weights_matrix, dW);
    matrix_destroy(dW);
    
    // Calculate new b vector from delta
    vector_t* db = matrix_sum(delta, 0);
    vectors_add(true, true, 2, l->b_vector, db);
    vector_destroy(db);

    // Calculate new delta matrix
    matrix_t* delta_new = matrix_dot_matrix(delta, l->weights_matrix);
    temp = l->act_fun.calculate_derivative(prev_prop);
    matrix_mult_matrix(delta_new, true, temp);
    matrix_destroy(temp);

    return delta_new;
}

// Creates new neural network with specified parameters: input_dim, epsilon, 
// regularization lambda, activation function, number of layers and their sizes
// - Last layer is the output layer
neural_network_t* neural_network_create(uint16_t input_dim, double epsilon, double reg_lambda, activation_function_t act_fun, uint16_t num_layers, ...)
{
    neural_network_t* nn = (neural_network_t*) malloc(sizeof(neural_network_t));
    nn->epsilon = epsilon;
    nn->reg_lambda = reg_lambda;
    nn->num_layers = num_layers;
    nn->layers = (nn_layer_t**) malloc(sizeof(nn_layer_t*) * num_layers);
    va_list valist;
    va_start(valist, num_layers);
    uint16_t num_neurons;
    uint16_t last_dim = input_dim;
    for(uint16_t i = 0; i < num_layers; i++) {
        num_neurons = va_arg(valist, int);
        if (i == num_layers-1) {
            act_fun = (activation_function_t) { .calculate = __LAMBDA__(double, (double x) { return x; }), .calculate_derivative = __LAMBDA__(matrix_t*, (matrix_t* m) { return matrix_create_sized(m->num_rows, m->row_vectors[0]->len, false, 1.0); }) };
        }
        nn->layers[i] = nn_layer_create(last_dim, num_neurons, act_fun);
        last_dim = num_neurons;
    }
    va_end(valist);
    return nn;
}

// Destroys neural network and frees allocated space
void neural_network_destroy(neural_network_t* nn)
{
    for(uint16_t i = 0; i < nn->num_layers; i++) {
        nn_layer_destroy(nn->layers[i]);
    }
    free(nn->layers);
    free(nn);
}

// Trains model of neural network, returns final loss
double neural_network_train_model(neural_network_t* nn, matrix_t* input_set, matrix_t* target_set, uint16_t num_passes)
{
    double model_loss = 0; // TODO calc Loss function

    matrix_t* temp;
    matrix_t** fwd_prop_outputs = (matrix_t**) malloc(sizeof(matrix_t*) * (nn->num_layers + 1));
    fwd_prop_outputs[0] = input_set;
    
    for (uint16_t i = 0; i < num_passes; i++) {
        // Forward propagation
        for (uint16_t n = 0; n < nn->num_layers; n++) {
            fwd_prop_outputs[n+1] = matrix_init(input_set->num_rows);
            for (uint16_t j = 0; j < input_set->num_rows; j++) {
                fwd_prop_outputs[n+1]->row_vectors[j] = nn_layer_forward_prop(nn->layers[n], fwd_prop_outputs[n]->row_vectors[j]);
            }
        }
        matrix_row_softmax(fwd_prop_outputs[nn->num_layers], true);

        // Backward propagation
        matrices_add(true, true, 2, fwd_prop_outputs[nn->num_layers], target_set);
        matrix_t* delta = fwd_prop_outputs[nn->num_layers];
        matrix_t* delta_new;
        for (uint16_t n = nn->num_layers-1; n >= 0; n--) {
            delta_new = nn_layer_backward_prop(nn->layers[n], delta, nn->reg_lambda, nn->epsilon, fwd_prop_outputs[n]);
            matrix_destroy(delta);
            if (n == 0) {
				matrix_destroy(delta_new);
				break;
			}
			delta = delta_new;
        }

        for (uint16_t n = 1; n < nn->num_layers; n++) {
            matrix_destroy(fwd_prop_outputs[n]);
        }

        printf("Training passed iteration #%d\n", i);
    }

    free(fwd_prop_outputs);

    return model_loss;
}

int main(int argc, char *argv[])
{
    matrix_t* input_set = matrix_create(10,
    vector_create(2, 1., 1.),
    vector_create(2, 3., 5.),
    vector_create(2, 2., 0.),
    vector_create(2, 1., 5.),
    vector_create(2, 1., 3.),
    vector_create(2, 10., 1.),
    vector_create(2, 30., 5.),
    vector_create(2, 20., 0.),
    vector_create(2, 10., 5.),
    vector_create(2, 10., 3.)
    );

    matrix_t* labels = matrix_create(10,
    vector_create(2, 0., 1.),
    vector_create(2, 0., 1.),
    vector_create(2, 1., 0.),
    vector_create(2, 1., 0.),
    vector_create(2, 1., 0.),
    vector_create(2, 1., 0.),
    vector_create(2, 1., 0.),
    vector_create(2, 1., 0.),
    vector_create(2, 1., 0.),
    vector_create(2, 0., 1.)
    );

    printf("Creating new NN.\n");
    neural_network_t* nn_test = neural_network_create(2, 0.01, 0.01, fa_tanh, 2, 10, 2);
    printf("Training NN.\n");
    double final_loss = neural_network_train_model(nn_test, input_set, labels, 70);
    printf("NN trained with loss = %f\n", final_loss);
    neural_network_destroy(nn_test);

    matrix_destroy(input_set);
    matrix_destroy(labels);
}