A Simple NN

// A simple neural network framework

// TODO: Add back propagation

// Basic assumptions:
//  1. Fully connected NN only.
//  2. Each neuron in the output layer has only one output, so output layer has as many neurons as the outputs.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>



float ti [] = {
0, 1, 2, 3, 4, 5
};

float to [] = {
0, 2, 4, 6, 8, 10
};



typedef struct {
    int numRows;
    int numCols;
    float *elements;
} Matrix;

#define MAT_AT(m, i, j) ((m)->elements[(i)*(m)->numCols+(j)])



typedef struct {
    int numInputs; // Number of inputs to the NN
    int numOutputs; // Number of output neurons in output layer
    int numHiddenLayers; // Number of hidden layers. Can be 0
    int *pnWidth; // Width of each hidden layer saved as a pointer to a list of width. numHiddenLayers is the number of layers.

    Matrix matInputs; // Single value for each input i.e. x1, x2.. etc. Can be represented by a single row [x1, x2...]
    Matrix *matLayerWeights; // Weights of layers.
    Matrix *matLayerBiases; // Biases of layers
    Matrix *matLayerOutputs; // Intermediate values from Layers after passing in values from previous layers
} NN;


void matrix_print(Matrix *m);

void matrix_fill_zero(Matrix *m)
{
    for(int i=0; i < m->numRows; i++) {
        for(int j=0; j < m->numCols; ++j){
            int k = i*m->numCols + j;
                m->elements[k] = 0;
        }
    }
}


Matrix *matrix_create(int r, int c)
{
    Matrix *m = NULL;

    m = malloc(sizeof(m));
    assert(m != NULL);

    m->numRows = r;
    m->numCols = c;
    m->elements = NULL;
    m->elements = malloc(sizeof(m->elements)*r*c);
    assert(m->elements != NULL);

    matrix_fill_zero(m);

    return m;
}



void matrix_dot(Matrix *dst, Matrix *a, Matrix *b)
{
    float n, sum;

    assert(a->numCols == b->numRows);
    assert(dst->numRows == a->numRows);
    assert(dst->numCols == b->numCols);

    n = a->numCols;
    for(int i=0; i<a->numRows; ++i) {
        for(int j=0; j<b->numCols; ++j) {
            sum = 0;
            for(int k=0; k<n; k++) {
                sum += MAT_AT(a, i, k) * MAT_AT(b, k, j);
            }
            MAT_AT(dst, i, j) = sum;
        }
    }
}


void matrix_add(Matrix *dst, Matrix *a, Matrix *b)
{
    float sum;

    assert((a->numCols == b->numCols) && (b->numCols == dst->numCols));
    assert((a->numRows == b->numRows) && (b->numRows == dst->numRows));

    for(int i=0; i<a->numRows; ++i) {
        for(int j=0; j<a->numCols; ++j) {
            sum = MAT_AT(a, i, j) + MAT_AT(b, i, j);
            MAT_AT(dst, i, j) = sum;
        }
    }
}



void matrix_print(Matrix *m)
{
    printf("matrix_print: Rows: %d Cols: %d\n", m->numRows, m->numCols);
    for(int i=0; i<m->numRows; i++) {
        for(int j=0; j<m->numCols; ++j){
            printf("%f ", MAT_AT(m, i, j));
        }
        printf("\n");
    }
}



void matrix_randomize(Matrix *m, float low, float high)
{
    for(int i=0; i<m->numRows; i++) {
        for(int j=0; j<m->numCols; j++) {
            MAT_AT(m, i, j) = rand()/(float)RAND_MAX;
        }
    }
}



NN NN_create(int inputs, int outputs, int layers, int *width)
{
    NN nn, *tnn;
    Matrix  *mw, *mb, *mo;

    nn.numInputs = inputs;
    nn.numOutputs = outputs;
    nn.numHiddenLayers = layers;
    nn.pnWidth = width;

    if (layers < 0) { // Can't have negative number of layers
        assert(layers >= 0);
        return nn;

    } else if (layers == 0) { // No hidden layer. Only input and output layer. So basically only as many neurons as the output layer.
        nn.matInputs = *matrix_create(1, inputs);

        // As there are no hidden layers, there is only going to be one layer i.e. the output layer
        mw = matrix_create(inputs, outputs);
        mb = matrix_create(1, outputs);
        mo = matrix_create(1, outputs);

        nn.matLayerWeights = mw;
        nn.matLayerBiases = mb;
        nn.matLayerOutputs = mo;
        nn.pnWidth = NULL; //Width should be NULL for neural networks with no hidden layers

    } else { // Hidden layers are 1 or more than 1.
        Matrix *m;
        m = matrix_create(1, inputs);
        nn.matInputs = *m;

        if (layers == 1) {
            // Layer 1 is hidden and the next is the output layer. So 2 layers.
            nn.matLayerWeights = malloc(sizeof(nn.matLayerWeights)*(layers + outputs));
            nn.matLayerBiases = malloc(sizeof(nn.matLayerBiases)*(layers + outputs));
            nn.matLayerOutputs = malloc(sizeof(nn.matLayerOutputs)*(layers + outputs));

            mw = matrix_create(inputs, *width);
            mb = matrix_create(1, *width);
            mo = matrix_create(1, *width);
            nn.matLayerWeights[0] = *mw;
            nn.matLayerBiases[0] = *mb;
            nn.matLayerOutputs[0] =  *mo;

            mw = matrix_create(*width, outputs);
            mb = matrix_create(1, outputs);
            mo = matrix_create(1, outputs);
            nn.matLayerWeights[1] = *mw;
            nn.matLayerBiases[1] = *mb;
            nn.matLayerOutputs[1] =  *mo;
        } else {

            nn.matLayerWeights = malloc(sizeof(nn.matLayerWeights)*(layers + 10));
            nn.matLayerBiases = malloc(sizeof(nn.matLayerBiases)*(layers + 10));
            nn.matLayerOutputs = malloc(sizeof(nn.matLayerOutputs)*(layers + 10));


                // First hidden layer's number of rows is equal to the number of columns of input layer
                mw = matrix_create(inputs, *width);
                mb = matrix_create(1, *width);
                mo = matrix_create(1, *width);

                nn.matLayerWeights[0] = *mw;
                nn.matLayerBiases[0] = *mb;
                nn.matLayerOutputs[0] =  *mo;

            for (int i=1; i<layers; i++) {
                mw = matrix_create(width[i-1], width[i]);
                mb = matrix_create(1, width[i]);
                mo = matrix_create(1, width[i]);

                nn.matLayerWeights[i] = *mw;
                nn.matLayerBiases[i] = *mb;
                nn.matLayerOutputs[i] =  *mo;

            }

            // Now the last layer i.e. output layer
            mw = matrix_create(width[layers-1], outputs);
            mb = matrix_create(1, outputs);
            mo = matrix_create(1, outputs);

            nn.matLayerWeights[layers] = *mw;
            nn.matLayerBiases[layers] = *mb;
            nn.matLayerOutputs[layers] = *mo;
        }
    }
    return nn;
}



void NN_print(NN nn)
{
    printf("NN_Print: Printing Neural network\nInputs: %d Outputs: %d Layers: %d Width: ", nn.numInputs, nn.numOutputs, nn.numHiddenLayers);

    if (nn.numHiddenLayers == 0)
        printf("0\n");
    else
        printf("%d \n", nn.numHiddenLayers);

    printf("NN_Print: Input layer: ****\n");
    matrix_print(&nn.matInputs);

    for(int i=0; i<nn.numHiddenLayers+1; ++i) {
        printf("NN_Print: Printing layer weights: %d***\n", i);
        matrix_print(&nn.matLayerWeights[i]);

        printf("NN_Print: Printing biases: %d***\n", i);
        matrix_print(&nn.matLayerBiases[i]);

        printf("NN_Print: Printing intermediate: %d***\n", i);
        matrix_print(&nn.matLayerOutputs[i]);

        printf("****************************\n");
    }
}



void NN_randomize(NN nn)
{
    // We randomize weights and biases of the hidden and output layers
    for(int i=0; i<nn.numHiddenLayers+1; i++) {
        matrix_randomize(&nn.matLayerWeights[i], 0, 1);
        matrix_randomize(&nn.matLayerBiases[i], 0, 1);
    }
}


// Function to feed values to NN. As many inputs as numInputs at a time.
void NN_feed(NN nn, float *inputs, float inputLen)
{
    assert(inputLen == nn.numInputs);

    Matrix *m = &nn.matInputs;

    for(int i=0; i<inputLen; i++) {
        MAT_AT(m, 0, i) = inputs[i];
    }

    if (nn.numHiddenLayers == 0) {

        // No hidden layers. Output can be directly computed.
        matrix_dot(nn.matLayerOutputs, &nn.matInputs, nn.matLayerWeights);
        matrix_add(nn.matLayerOutputs, nn.matLayerOutputs, nn.matLayerBiases);

    } else if (nn.numHiddenLayers == 1) { // 1 Hidden Layer

        matrix_dot(&nn.matLayerOutputs[0], &nn.matInputs, &nn.matLayerWeights[0]);
        matrix_add(&nn.matLayerOutputs[0], &nn.matLayerOutputs[0], &nn.matLayerBiases[0]);

        matrix_dot(&nn.matLayerOutputs[1], &nn.matLayerOutputs[0], &nn.matLayerWeights[1]);
        matrix_add(&nn.matLayerOutputs[1], &nn.matLayerOutputs[1], &nn.matLayerBiases[1]);

    } else { // More than 1 Hidden Layers

        matrix_dot(&nn.matLayerOutputs[0], &nn.matInputs, &nn.matLayerWeights[0]);
        matrix_add(&nn.matLayerOutputs[0], &nn.matLayerOutputs[0], &nn.matLayerBiases[0]);

        for(int i=1; i<nn.numHiddenLayers + 1; i++) {
            matrix_dot(&nn.matLayerOutputs[i], &nn.matLayerOutputs[i-1], &nn.matLayerWeights[i]);
            matrix_add(&nn.matLayerOutputs[i], &nn.matLayerOutputs[i], &nn.matLayerBiases[i]);
        }
    }
}



float NN_cost(NN nn)
{
    // Calculate the cost of NN
    // Input all the rows of the input data and compare with output

    int n = sizeof(ti)/sizeof(ti[0]);
    float cost = 0.0f, y = 0.0f;


    for (int i=0; i<n; i++) {
        float d = 0.0f;
        NN_feed(nn, &ti[i], 1);
        for(int j=0; j<nn.numOutputs; j++) {
            y = MAT_AT(&nn.matLayerOutputs[nn.numHiddenLayers], 0, j);
            d = to[i*nn.numOutputs+j] - y;
            cost += d*d;
        }
    }

    return cost/n;
}



void NN_finite_difference(NN nn, NN g, float eps)
{
    // Here we calculate the difference in costs of original NN and changing one weight by epsilon.
    // nn is the original neural network, g is the neural network that contains the finite differences
    // Later we'll use NN_learn_with_rate to make changes to the original neural network

    // We'll cycle through the layers to the output layer and change one weight at a time and calculate the cost
    // We have to do both weights and biases

    float saved = 0;
    float cost = NN_cost(nn); // Save this cost as this will be required later
    for(int i=0; i<nn.numHiddenLayers+1; i++) {

        // First for weights
        for(int j=0; j<nn.matLayerWeights[i].numRows; j++) {
            for(int k=0; k<nn.matLayerWeights[i].numRows; k++) {
                saved = MAT_AT(&nn.matLayerWeights[i], j, k);
                MAT_AT(&nn.matLayerWeights[i], j, k) += eps;
                MAT_AT(&g.matLayerWeights[i], j, k) = (NN_cost(nn) - cost)/eps;
                MAT_AT(&nn.matLayerWeights[i], j, k) = saved;
            }
        }

        // Now for biases
        for(int j=0; j<nn.matLayerBiases[i].numRows; j++) {
            for(int k=0; k<nn.matLayerBiases[i].numRows; k++) {
                saved = MAT_AT(&nn.matLayerBiases[i], j, k);
                MAT_AT(&nn.matLayerBiases[i], j, k) += eps;
                MAT_AT(&g.matLayerBiases[i], j, k) = (NN_cost(nn) - cost)/eps;
                MAT_AT(&nn.matLayerBiases[i], j, k) = saved;
            }
        }
    }
}



void NN_learn_with_rate(NN nn, NN g, float rate)
{
    for(int i=0; i<nn.numHiddenLayers+1; i++) {

        // First for weights
        for(int j=0; j<nn.matLayerWeights[i].numRows; j++) {
            for(int k=0; k<nn.matLayerWeights[i].numRows; k++) {
                MAT_AT(&nn.matLayerWeights[i], j, k) -= rate * MAT_AT(&g.matLayerWeights[i], j, k);
            }
        }

        // Now for biases
        for(int j=0; j<nn.matLayerBiases[i].numRows; j++) {
            for(int k=0; k<nn.matLayerBiases[i].numRows; k++) {
                MAT_AT(&nn.matLayerBiases[i], j, k) -= rate * MAT_AT(&g.matLayerBiases[i], j, k);
            }
        }
    }
}



int main() {

    int input = 1;
    int numberOfHiddenLayers = 1;
    int widthOfHiddenLayers[] = {2, 3, 4};
    int output = 1;
    float eps = 1e-3, rate = 1e-3;

    int n = sizeof(ti)/sizeof(ti[0]);

    srand(69);

    NN nn = NN_create(input, output, numberOfHiddenLayers, widthOfHiddenLayers);

    NN g = NN_create(input, output, numberOfHiddenLayers, widthOfHiddenLayers);

    NN_randomize(nn);

printf("Initial Cost :%f\n", NN_cost(nn));

    for(int i=0; i<10000; ++i) {
//printf("Cost %d :%f\n", i, NN_cost(nn));
NN_finite_difference(nn, g, eps);
NN_learn_with_rate(nn, g, rate);
    }

printf("Final Cost :%f\n", NN_cost(nn));

for(int i=0; i<n; i++) {
  NN_feed(nn, &ti[i], 1);
  printf("ti: %f to: %f NN Output: %f\n", ti[i], to[i], MAT_AT(&nn.matLayerOutputs[nn.numHiddenLayers], 0, 0));
}

NN_print(nn);
    // TODO: Free the allocated memory!
    return 0;
}