多層パーセプトロンのC言語による実装

多層パーセプトロンのC言語による実装を示します。

この実装は、下記のWebページにあった実装をもとにしています。
http://www.ccodechamp.com/c-program-of-multilayer-perceptron-net-using-backpropagation/

主な変更点は以下のとおりです。
1. ReLUを追加。
2. 最適化アルゴリズムをAdaGradに変更（[An overview of gradient descent optimization algorithms](http://sebastianruder.com/optimizing-gradient-descent/index.html)）。
3. 乱数をSFMTに変更（[SIMD-oriented Fast Mersenne Twister (SFMT)](http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/)）。
4. いくつかの固定長の配列を動的配列に変更。
5. loadNet関数を削除。

```c:

#include <stdio.h>
#include <math.h>
#include "SFMT-src-1.4.1/SFMT.h"

#define DTYPE float
#define M_PI 3.14159265358979323846

//================================================================
// http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/
sfmt_t *sfmt;
unsigned int initRand(unsigned int nSeed)
{
  sfmt = (sfmt_t *) malloc(sizeof(sfmt_t));
  if (sfmt == NULL) return 1;
  sfmt_init_gen_rand(sfmt, nSeed);
  return 0;
}
DTYPE genrand_real1() // [0,1]
{ return (DTYPE) sfmt_genrand_real1(sfmt); }
DTYPE genrand_real2() // [0,1)
{ return (DTYPE) sfmt_genrand_real2(sfmt); }
DTYPE genrand_normal() // Box-Muller
{ return (DTYPE) sqrt(- 2.0 * log(1.0 - sfmt_genrand_real2(sfmt)))
    * cos(2.0 * M_PI * (1.0 - sfmt_genrand_real2(sfmt))); }
//================================================================

struct neuron{
  DTYPE *output;
  char axonFamily;
  DTYPE *weights;
  DTYPE *weightGrad;
  DTYPE *weightHistoricalGrad;
  DTYPE *bias;
  DTYPE *biasGrad;
  DTYPE *biasHistoricalGrad;
  int noOfInputs;
  DTYPE *inputs;
  DTYPE *error;
};

struct layer {
  int noOfNeurons;
  struct neuron *neurons;
};

struct neuralNet {
  int noOfInputs;
  DTYPE *inputs;
  DTYPE *outputs;
  int noOfLayers;
  struct layer *layers;
  int noOfBatchChanges;
  int noOfIterations;
} theNet;

DTYPE weightInit() { return 0.3 * genrand_normal(0); }
DTYPE biasInit() { return 0.1; }
DTYPE getRand() { return 2.0 * M_PI * (2.0 * genrand_real1(0) - 1.0); }

int totalNoOfNeurons;
int totalNoOfWeights;

struct neuron *netNeurons;
DTYPE *netInputs;
DTYPE *netNeuronOutputs;
DTYPE *netErrors;
struct layer *netLayers;
DTYPE *netWeights;
DTYPE *netWeightGrad;
DTYPE *netWeightHistoricalGrad;
DTYPE *netBiases;
DTYPE *netBiasGrad;
DTYPE *netBiasHistoricalGrad;

void createNet(int noOfLayers, int *noOfNeurons, int *noOfInputs, char *axonFamilies, int initWeights) {

  int i, j, counter, counter2, counter3, counter4;

  theNet.noOfLayers = noOfLayers;
  netLayers = (struct layer *) malloc(sizeof(struct layer) * theNet.noOfLayers);
  theNet.layers = netLayers;

  theNet.noOfInputs = noOfInputs[0];
  netInputs = (DTYPE *) malloc(sizeof(DTYPE) * theNet.noOfInputs);
  theNet.inputs = netInputs;

  totalNoOfNeurons = 0;
  for(i = 0; i < theNet.noOfLayers; i++) {
    totalNoOfNeurons += noOfNeurons[i];
  }
  netNeurons = (struct neuron *) malloc(sizeof(struct neuron) * totalNoOfNeurons);
  netNeuronOutputs = (DTYPE *) malloc(sizeof(DTYPE) * totalNoOfNeurons);
  netErrors = (DTYPE *) malloc(sizeof(DTYPE) * totalNoOfNeurons);
  netBiases = (DTYPE *) malloc(sizeof(DTYPE) * totalNoOfNeurons);
  netBiasGrad = (DTYPE *) malloc(sizeof(DTYPE) * totalNoOfNeurons);
  netBiasHistoricalGrad = (DTYPE *) malloc(sizeof(DTYPE) * totalNoOfNeurons);
  for(i = 0; i < totalNoOfNeurons; i++) { netNeuronOutputs[i] = 0.0; }

  totalNoOfWeights = 0;
  for(i = 0; i < theNet.noOfLayers; i++) {
    totalNoOfWeights += noOfInputs[i] * noOfNeurons[i];
  }
  netWeights = (DTYPE *) malloc(sizeof(DTYPE) * totalNoOfWeights);
  netWeightGrad = (DTYPE *) malloc(sizeof(DTYPE) * totalNoOfWeights);
  netWeightHistoricalGrad = (DTYPE *) malloc(sizeof(DTYPE) * totalNoOfWeights);

  counter = 0;
  counter2 = 0;
  counter3 = 0;
  counter4 = 0;
  for(i = 0; i < theNet.noOfLayers; i++) {
    for(j = 0; j < noOfNeurons[i]; j++) {
      if(i == theNet.noOfLayers-1 && j == 0) { // beginning of the output layer
        theNet.outputs = &netNeuronOutputs[counter];
      }
      netNeurons[counter].output = &netNeuronOutputs[counter];
      netNeurons[counter].axonFamily = axonFamilies[i];
      netNeurons[counter].weights = &netWeights[counter2];
      netNeurons[counter].weightGrad = &netWeightGrad[counter2];
      netNeurons[counter].weightHistoricalGrad = &netWeightHistoricalGrad[counter2];
      netNeurons[counter].bias = &netBiases[counter];
      netNeurons[counter].biasGrad = &netBiasGrad[counter];
      netNeurons[counter].biasHistoricalGrad = &netBiasHistoricalGrad[counter];
      netNeurons[counter].noOfInputs = noOfInputs[i];
      if (i == 0) {
        netNeurons[counter].inputs = netInputs;
      }
      else {
        netNeurons[counter].inputs = &netNeuronOutputs[counter3];
      }
      netNeurons[counter].error = &netErrors[counter];
      counter2 += noOfInputs[i];
      counter++;
    }
    netLayers[i].noOfNeurons = noOfNeurons[i];
    netLayers[i].neurons = &netNeurons[counter4];
    if (i > 0) {
      counter3 += noOfNeurons[i-1];
    }
    counter4 += noOfNeurons[i];
  }

  if (initWeights == 1) {
    for(i = 0; i < totalNoOfNeurons; i++) { netBiases[i] = biasInit(); }
    for(i = 0; i < totalNoOfNeurons; i++) { netBiasHistoricalGrad[i] = 0.0; }
    for(i = 0; i < totalNoOfWeights; i++) { netWeights[i] = weightInit(); }
    for(i = 0; i < totalNoOfWeights; i++) { netWeightHistoricalGrad[i] = 0.0; }
  }

  for(i = 0; i < totalNoOfNeurons; i++) { netBiasGrad[i] = 0.0; }
  for(i = 0; i < totalNoOfWeights; i++) { netWeightGrad[i] = 0.0; }
  theNet.noOfBatchChanges = 0;
  theNet.noOfIterations = 0;

}

void feedNetInputs(DTYPE *inputs) {

  int i;
  for ( i = 0; i < theNet.noOfInputs; i++ ) {
    netInputs[i] = inputs[i];
  }

}

void updateNeuronOutput(struct neuron *myNeuron) {

  DTYPE activation = 0.0;
  int i;

  for (i = 0; i < myNeuron->noOfInputs; i++) {
    activation += myNeuron->inputs[i] * myNeuron->weights[i];
  }
  activation += *(myNeuron->bias);
  DTYPE temp;
  switch (myNeuron->axonFamily) {
  case 'g': // sigmoid
    temp = - activation;
    if (temp > 45.0) { *(myNeuron->output) = 0.0; }
    else if (temp < -45.0) { *(myNeuron->output) = 1.0; }
    else { *(myNeuron->output) = 1.0 / (1.0 + exp(temp)); }
    break;
  case 't': // tanh
    temp = - activation;
    if (temp > 45.0) { *(myNeuron->output) = -1.0; }
    else if (temp < -45.0) { *(myNeuron->output) = 1.0; }
    else { *(myNeuron->output) = (2.0 / (1.0 + exp(temp))) - 1.0; }
    break;
  case 'r': // ReLU
    if (activation > 0.0) { *(myNeuron->output) = activation; }
    else { *(myNeuron->output) = 0.0; }
    break;
  default:
    break;
  }

}

void updateNetOutput() {

  int i, j;
  for(i = 0; i < theNet.noOfLayers; i++) {
    for(j = 0; j < theNet.layers[i].noOfNeurons; j++) {
      updateNeuronOutput(&(theNet.layers[i].neurons[j]));
    }
  }

}

DTYPE derivative(struct neuron *myNeuron) {

  DTYPE temp;
  switch (myNeuron->axonFamily) {
  case 'g': // sigmoid
    temp = *(myNeuron->output) * (1.0 - *(myNeuron->output)); break;
  case 't': // tanh
    temp = 1.0 - 0.5 * *(myNeuron->output) * *(myNeuron->output); break;
  case 'r': // ReLU
    if (*(myNeuron->output) > 0.0) { temp = 1.0; }
    else { temp = 0.0; }
    break;
  default:
    temp = 0; break;
  }

  return temp;

}

void trainNet(DTYPE *outputTargets) {

  int i,j,k;
  DTYPE temp;
  struct layer *currLayer, *nextLayer;

  // calculate errors
  for(i = theNet.noOfLayers - 1; i >= 0; i--) {
    currLayer = &theNet.layers[i];
    if (i == theNet.noOfLayers - 1 ) { // output layer
      for (j = 0; j < currLayer->noOfNeurons; j++ ) {
        *(currLayer->neurons[j].error)
          = derivative(&currLayer->neurons[j]) * (outputTargets[j] - *(currLayer->neurons[j].output));
      }
    }
    else { // other layers
      nextLayer = &theNet.layers[i+1];
      for (j = 0; j < currLayer->noOfNeurons; j++ ) {
        temp = 0.0;
        for (k = 0; k < nextLayer->noOfNeurons; k++ ) {
          temp += *(nextLayer->neurons[k].error) * nextLayer->neurons[k].weights[j];
        }
        *(currLayer->neurons[j].error) = derivative(&currLayer->neurons[j]) * temp;
      }
    }
  }

  // update weights n biases
  for(i = theNet.noOfLayers - 1; i >= 0; i--) {
    currLayer = &theNet.layers[i];
    for (j = 0; j < currLayer->noOfNeurons; j++) {
      *(currLayer->neurons[j].biasGrad) += *(currLayer->neurons[j].error);
      // weights
      for(k = 0; k < currLayer->neurons[j].noOfInputs; k++) {
        currLayer->neurons[j].weightGrad[k]
          += *(currLayer->neurons[j].error) * currLayer->neurons[j].inputs[k];
      }
    }
  }

  theNet.noOfBatchChanges++;

}

void setGradient(DTYPE alpha) {

  int i,j,k;
  struct layer *currLayer;

  for(i = theNet.noOfLayers - 1; i >= 0; i--) {
    currLayer = &theNet.layers[i];
    for (j = 0; j < currLayer->noOfNeurons; j++) {
      // biases
      *(currLayer->neurons[j].biasGrad) /= theNet.noOfBatchChanges;
      if (theNet.noOfIterations == 0) {
        *(currLayer->neurons[j].biasHistoricalGrad)
          = *(currLayer->neurons[j].biasGrad) * *(currLayer->neurons[j].biasGrad);
      }
      else {
        *(currLayer->neurons[j].biasHistoricalGrad)
          = alpha * *(currLayer->neurons[j].biasHistoricalGrad)
          + (1.0 - alpha) * *(currLayer->neurons[j].biasGrad) * *(currLayer->neurons[j].biasGrad);
      }
      // weights
      for (k = 0; k < currLayer->neurons[j].noOfInputs; k++) {
        currLayer->neurons[j].weightGrad[k] /= theNet.noOfBatchChanges;
        if (theNet.noOfIterations == 0) {
          currLayer->neurons[j].weightHistoricalGrad[k]
            = currLayer->neurons[j].weightGrad[k] * currLayer->neurons[j].weightGrad[k];
        }
        else {
          currLayer->neurons[j].weightHistoricalGrad[k]
            = alpha * currLayer->neurons[j].weightHistoricalGrad[k]
            + (1.0 - alpha) * currLayer->neurons[j].weightGrad[k] * currLayer->neurons[j].weightGrad[k];
        }
      }
    }
  }

}

void update(DTYPE stepsize) {

  int i,j,k;
  struct layer *currLayer;
  DTYPE fudgeFactor = 1e-8;

  for(i = theNet.noOfLayers - 1; i >= 0; i--) {
    currLayer = &theNet.layers[i];
    for (j = 0; j < currLayer->noOfNeurons; j++) {
      // biases
      *(currLayer->neurons[j].bias)
        += stepsize * (*(currLayer->neurons[j].biasGrad)
                       / sqrt(fudgeFactor + *(currLayer->neurons[j].biasHistoricalGrad)));
      *(currLayer->neurons[j].biasGrad) = 0.0;
      // weights
      for (k = 0; k < currLayer->neurons[j].noOfInputs; k++) {
        currLayer->neurons[j].weights[k]
          += stepsize * (currLayer->neurons[j].weightGrad[k]
                         / (sqrt(fudgeFactor + currLayer->neurons[j].weightHistoricalGrad[k])));
        currLayer->neurons[j].weightGrad[k] = 0.0;
      }
    }
  }

  theNet.noOfBatchChanges = 0;
  theNet.noOfIterations++;

}

DTYPE *getOutputs() {

  return theNet.outputs;

}

int main() {

  if (initRand(0)) return 1;

  /* determine layer paramaters */
  int noOfLayers = 4; // input layer excluded
  int noOfNeurons[] = {20,20,20,2};
  int noOfInputs[] = {2,20,20,20};
  char axonFamilies[] = {'r','r','r','t'};

  DTYPE *inputs;
  DTYPE *outputTargets;

  inputs = (DTYPE *) malloc(sizeof(DTYPE) * noOfInputs[0]);
  outputTargets = (DTYPE *) malloc(sizeof(DTYPE) * noOfNeurons[noOfLayers - 1]);

  createNet(noOfLayers, noOfNeurons, noOfInputs, axonFamilies, 1);

  /* train it using batch method */
  int i;
  DTYPE tempTotal;
  int counter = 0;
  for(i = 0; i < 10000000; i++) {
    inputs[0] = getRand();
    inputs[1] = getRand();
    tempTotal = inputs[0] + inputs[1];
    feedNetInputs(inputs);
    updateNetOutput();
    outputTargets[0] = (DTYPE) sin(tempTotal);
    outputTargets[1] = (DTYPE) cos(tempTotal);
    trainNet(outputTargets);
    counter++;
    if (counter == 50) { setGradient(0.9); update(0.001); counter = 0; }
  }

  /* test it */
  DTYPE *outputs;
  printf("Sin Target \t Output \t Cos Target \t Output\n");
  printf("---------- \t -------- \t ---------- \t --------\n");
  for(i = 0; i < 50; i++) {
    inputs[0] = getRand();
    inputs[1] = getRand();
    tempTotal = inputs[0] + inputs[1];
    feedNetInputs(inputs);
    updateNetOutput();
    outputs = getOutputs();
    printf( "%f \t %f \t %f \t %f \n", sin(tempTotal), outputs[0], cos(tempTotal), outputs[1]);
  }

  return 0;

}
```

実行結果は以下のような感じになると思います。

```shell-session
$ ./mlp
Sin Target       Output          Cos Target      Output
----------       --------        ----------      --------
0.640082         0.746855        -0.768307       -0.812494
0.069101         0.079599        -0.997610       -0.988040
-0.985831        -0.984480       -0.167741       -0.157381
0.521136         0.532716        0.853474        0.864914
-0.097435        -0.109023       0.995242        0.994873
-0.114213        -0.112507       0.993456        0.992926
-0.917489        -0.925177       -0.397762       -0.380658
-0.046450        -0.051414       0.998921        0.996574
-0.692018        -0.679392       -0.721880       -0.727636
0.409927         0.424983        -0.912118       -0.924427
-0.916550        -0.913544       0.399921        0.423543
-0.824584        -0.721065       0.565740        0.638374
0.023286         0.022497        0.999729        0.967530
-0.650853        -0.650613       0.759204        0.782016
0.972549         0.942901        -0.232696       -0.201588
0.014787         -0.000262       -0.999891       -0.962320
-0.969422        -0.968556       0.245399        0.293858
0.338011         0.355989        -0.941142       -0.948886
-0.981046        -0.981085       -0.193777       -0.182955
0.520488         0.611798        0.853869        0.867387
0.958533         0.969545        -0.284982       -0.269451
-0.978963        -0.977859       -0.204036       -0.193433
-0.905485        -0.877378       0.424378        0.436468
-0.438091        -0.431923       -0.898930       -0.885023
-0.985005        -0.986799       0.172525        0.187228
-0.400110        -0.335487       -0.916467       -0.900038
-0.248896        -0.262504       0.968530        0.979315
0.443793         0.505794        -0.896129       -0.895520
-0.182260        -0.181806       -0.983250       -0.978943
0.322289         0.332658        0.946641        0.947475
-0.291678        -0.259068       0.956516        0.947655
0.102971         0.131531        -0.994684       -0.988836
-0.779955        -0.766913       0.625836        0.646452
0.906765         0.914458        0.421636        0.311500
0.877212         0.894052        0.480104        0.521253
-0.653435        -0.650393       0.756982        0.782935
-0.727792        -0.727932       0.685798        0.689874
-0.085042        -0.070495       -0.996377       -0.988692
0.124316         0.127393        0.992243        0.992283
-0.996733        -0.990036       -0.080773       -0.043595
-0.472917        -0.456472       -0.881107       -0.884031
-0.594023        -0.583654       0.804448        0.834040
0.980313         0.967607        0.197448        0.239313
-0.590370        -0.562181       0.807133        0.805516
-0.289410        -0.270091       -0.957205       -0.961424
-0.555648        -0.559167       0.831418        0.860375
0.642193         0.643480        0.766543        0.785267
0.486129         0.552793        -0.873887       -0.881352
-0.491824        -0.466695       0.870694        0.875552
-0.974766        -0.975019       -0.223228       -0.213617
```