Simple artifical neural network (2d multilayered)

/*
  Переписка genann на структурный лад.

  Полезная информация:
    Сигмоидная функция активации,
    Алгоритм обратного распространения ошибки (используется в network_train):
        http://www.aiportal.ru/articles/neural-networks/back-propagation.html
        https://ru.coursera.org/learn/vvedenie-mashinnoe-obuchenie/lecture/SjSlD/nieironnyie-sieti-vviedieniie
        https://ru.coursera.org/learn/vvedenie-mashinnoe-obuchenie/lecture/Lmbqi/nieironnyie-sieti-mietod-obratnogho-rasprostranieniia-oshibki
*/

/*
  TODO
  
  1) рефакторинг кода: уж слишком длинные выражения в циклах обучения сети
  2) сохранение в файл, извлечение состояния сети
*/

#include <stdio.h>  // printf
#include <stdlib.h> // rand, srand, malloc, calloc
#include <time.h>   // time
#include <math.h>   // exp


typedef struct network network;
typedef struct layer layer;
typedef struct neuron neuron;


struct neuron {
  double* weight;
  double output;
  double delta;
};


struct layer {
  neuron* neuron;
};


struct network {
  int error;
  int inputs, hidden_layers, neurons, outputs;
  layer* hidden;
  layer output;
};


/* Создает и возвращает нейрон с заданным кол-вом весов */
neuron neuron_create(int weights) {
  neuron neuron;

  /* calloc, в отличии от malloc, заполняет выделенную область нулями */
  /* Т.е. веса изначально будут равны нулю */
  neuron.weight = calloc(weights, sizeof(double));

  /* output и delta вычисляются позже, им не нужно сейчас присваивать значение */

  return neuron;
}


/* Создает и возвращает слой нейронов */
layer layer_create(int neurons, int weights) {
  layer layer;

  /* Выделяем место под массив из neurons нейронов */
  layer.neuron = malloc(neurons * sizeof(neuron));

  /* Заполняем массив нейронов новыми нейронами */
  /* ВНИМАНИЕ! У каждого нейрона есть один дополнительный вес с нулевым индексом,
     который всегда умножается на -1.0 */
  int n;
  for (n = 0; n < neurons; n++) {
    layer.neuron[n] = neuron_create(weights + 1);
  }

  return layer;
}


/* Возвращает случайное число от 0 до 1 */
double random_double() {
  /* Статическая переменная сохранит свое значение даже после выхода из функции */
  static int initalized = 0;

  /* Используем это для инициализации рандома (seed random) только при первом запуске функции */
  if (!initalized) {
    srand(time(NULL));
    initalized = 1;
  }

  return (double) rand() / RAND_MAX;
}


/* Рандомизирует веса нейронов сети. Вызывается в конце network_init */
void network_randomize(network* net) {
  int h, n, w;

  /* Скрытые слои (если есть) */
  for (h = 0; h < net->hidden_layers; h++) {
    for (n = 0; n < net->neurons; n++) {
      for (w = 0; w < 1 + (h == 0 ? net->inputs : net-> neurons); w++) {
        /* Задаем весам случайные значения от -0.5 до 0.5 */
        net->hidden[h].neuron[n].weight[w] = random_double() - 0.5;
      }
    }
  }

  /* Выходной слой */
  for (n = 0; n < net->outputs; n++) {
    for (w = 0; w < 1 + (net->hidden_layers ? net->neurons : net->inputs); w++) {
      net->output.neuron[n].weight[w] = random_double() - 0.5;
    }
  }
}


/* Создает и возвращает сеть с заданными параметрами */
network network_init(int inputs, int hidden_layers, int neurons, int outputs) {
  network net;

  /* Проверка входных параметров */
  net.error = (
    inputs < 1
    || hidden_layers < 0
    || (hidden_layers > 0 && neurons < 1)
    || outputs < 1
  );

  if (!net.error) {
    net.inputs = inputs;
    net.hidden_layers = hidden_layers;
    net.neurons = neurons;
    net.outputs = outputs;

    /* Выделяем память под массив скрытых слоёв (если задано 0, то и выделится 0) */
    net.hidden = malloc(hidden_layers * sizeof(layer));

    /* Заполняем массив слоёв */
    int h;
    for (h = 0; h < hidden_layers; h++) {
      net.hidden[h] = layer_create(neurons, h == 0 ? inputs : neurons);
    }

    /* Выходной слой */
    net.output = layer_create(outputs, neurons);

    network_randomize(&net);
  }

  return net;
}


/* Сигмоидная функция активации */
double act_sigmoid(double a) {
  return 1.0 / (1 + exp(-a));
}


/* Фидфорвард сети */
void network_run(network* net, double* input) {
  int h, n, w;

  /* Скрытые слои (если есть) */
  for (h = 0; h < net->hidden_layers; h++) {
    int first_hidden = h == 0;

    for (n = 0; n < net->neurons; n++) {
      double sum = 0;

      for (w = 0; w < 1 + (first_hidden ? net->inputs : net-> neurons); w++) {
        if (w == 0) {
          sum += net->hidden[h].neuron[n].weight[w] * -1.0;
        } else if (first_hidden) {
          sum += net->hidden[h].neuron[n].weight[w] * input[w];
        } else {
          sum += net->hidden[h].neuron[n].weight[w] * net->hidden[h-1].neuron[w].output;
        }
      }
      net->hidden[h].neuron[n].output = act_sigmoid(sum);
    }
  }

  /* Выходной слой */
  for (n = 0; n < net->outputs; n++) {
    double sum = 0;

    for (w = 0; w < 1 + (net->hidden_layers ? net->neurons : net->inputs); w++) {
      if (w == 0) {
        sum += net->output.neuron[n].weight[w] * -1.0;
      } else if (net->hidden_layers) {
        sum += net->output.neuron[n].weight[w] * net->hidden[net->hidden_layers-1].neuron[w].output;
      } else {
        sum += net->output.neuron[n].weight[w] * input[w];
      }
    }

    net->output.neuron[n].output = act_sigmoid(sum);
  }
}


/* Обучает сеть 1 раз */
void network_train(network* net, double* input, double* desired_output, double learning_rate) {
  /* Фидфорвардим сеть. Вычисляются выходные значения. */
  network_run(net, input);

  int h, n, w;

  /* Вычисляем разность между выходом сети и ожидаемым значением (сигнал ошибки) */
  /* Первый этап случая 1 */
  for (n = 0; n < net->outputs; n++) {
    double desired = desired_output[n];
    double output = net->output.neuron[n].output;
    net->output.neuron[n].delta = (desired - output) * output * (1.0 - output);
  }

  /* Распространяем сигнал ошибки к началу сети */
  /* Кстати, этот цикл будет пропущен если hidden_layers == 0 */
  /* http://www.aiportal.ru/images/articles/NN_back_propagation_algorithm_02.png
   * Первый этап случая 2.
   * forward_weight это wq-k на картинке,
   * forward_delta это дельта k-того нейрона на след. слое
   * delta это сумма дельт нейронов след. слоя */
  for (h = net->hidden_layers - 1; h >= 0; h--) {
    int last_hidden = h == net->hidden_layers - 1;

    for (n = 0; n < net->neurons; n++) {
      double delta_sum = 0;

      for (w = 0; w < (last_hidden ? net->outputs : net->neurons); w++) {
        double forward_delta = last_hidden ? net->output.neuron[w].delta : net->hidden[h+1].neuron[w].delta;
        double forward_weight = last_hidden ? net->output.neuron[w].weight[n] : net->hidden[h+1].neuron[w].weight[n];
        delta_sum += forward_delta * forward_weight;
      }

      double output = net->hidden[h].neuron[n].output;
      net->hidden[h].neuron[n].delta = output * (1.0-output) * delta_sum;
    }
  }

  /* Тренируем выходной слой */
  /* Второй этап случая 1 */
  for (n = 0; n < net->outputs; n++) {
    for (w = 0; w < 1 + (net->hidden_layers ? net->neurons : net->inputs); w++) {
      if (w == 0) {
        net->output.neuron[n].weight[w] += learning_rate * net->output.neuron[n].delta * -1.0;
      } else if (net->hidden_layers) {
        net->output.neuron[n].weight[w] += learning_rate * net->output.neuron[n].delta * net->hidden[net->hidden_layers-1].neuron[w].output;
      } else {
        net->output.neuron[n].weight[w] += learning_rate * net->output.neuron[n].delta * input[w];
      }
    }
  }

  /* Тренируем скрытые слои */
  /* Второй этап случая 2 */
  for (h = 0; h < net->hidden_layers; h++) {
    for (n = 0; n < net->neurons; n++) {
      for (w = 0; w < 1 + (h == 0 ? net->inputs : net->neurons); w++) {
        if (w == 0) {
          net->hidden[h].neuron[n].weight[w] += net->hidden[h].neuron[n].delta * -1.0;
        } else if (h == 0) {
          net->hidden[h].neuron[n].weight[w] += net->hidden[h].neuron[n].delta * input[w];
        } else {
          net->hidden[h].neuron[n].weight[w] += net->hidden[h].neuron[n].delta * net->hidden[h-1].neuron[w].output;
        }
      }
    }
  }
}


/* Отдает входные данные на обработку сетью, выводит результат */
void test(network* net, double* input, double* desired_output) {
  network_run(net, input);
  printf("Output for {%1.f, %1.f} is %1.f (should be %1.f)\n", input[0], input[1], net->output.neuron[0].output, desired_output[0]);
}


int main() {
  /* Будем учить нашу сеть выполнять функцию XOR */

  /* Создаем сеть с 2 входами
     1 скрытым слоем из 2 нейронов
     и 1 выходом */
  network net = network_init(2, 1, 2, 1);
  if (net.error) return 1;

  /* Для XOR в принципе достаточно иметь лишь выходной слой нейронов */
  /* Но мы покажем, что сеть может работеть и с несколькими слоями с различным кол-вом нейронов внутри */
  /* Кстати, при увеличении кол-ва скрытых слоев придется подкрутить и настройки обучения */

  double input[4][2] = {{0, 0}, {0, 1}, {1, 0}, {1, 1}};
  double desired_output[4][1] = {{0}, {1}, {1}, {0}};
  double learning_rate = 3;

  /* Тренируем сеть 300 раз */
  int i;
  for (i = 0; i < 300; i++) {
    network_train(&net, input[0], desired_output[0], learning_rate);
    network_train(&net, input[1], desired_output[1], learning_rate);
    network_train(&net, input[2], desired_output[2], learning_rate);
    network_train(&net, input[3], desired_output[3], learning_rate);
  }

  /* Проверяем результат обучения */
  for (i = 0; i < 4; i++) {
    test(&net, input[i], desired_output[i]);
  }

  return 0;
}