C言語でニューラルネット（関数近似）

C言語で関数近似を行うためのニューラルネットワーク（誤差逆伝播法，B ack P ropagation N eural N etwork）のプログラムです．

プログラム中で，Input_sizeは学習用データの入力の長さ，Input_dimは学習用データの次元，test_sizeは検証用データの入力の長さです．

たとえば，以下のような点の集合を関数で近似したい場合を考えます．

このとき，学習用データは次のように与えられます．

data.dat	rho.dat
0	1
1	3
2	2
3	4

というわけで，上のように実行ファイルの場所に，data.datとrho.datを作成します．
見て分かる通り，data.datがx軸，rho.datがy軸を表しています（一次元の場合）．

この2種類のファイルから，Input_sizeとInput_dimを設定します．
この場合，Input_sizeは，データが横軸に4つなので，Input_size=4とします．
また，入力は1次元なのでInput_dim=1です．

また，test_sizeは，学習が終わったあとのデータ点の数で，x軸の0~3を何分割するかを決める滑らかさの値です．
理想的な場合学習が起こった場合では，test_size=4とすると，上の図の点と同じ位置に点が打たれるはずです．
ここでは，どのように近似されたかを詳細に見るために，test_size=100とします．

これらのデータから実行した結果が以下の図．（学習回数1,000,000回，test_size=100）

このように，学習データとして与えた点を全て通るような曲線が描かれています．
これが，ニューラルネットワークによる関数近似です．

先ほど与えた点の値と，現在のニューラルネットワークの出力の誤差がどんどん小さくなっていくことで，与えた点を通る曲線が描かれています．
この曲線を拡大してみると，100個の点からなる直線の集合です（=test_sizeの値）．

test_sizeを小さくするとこんな感じに．（学習回数1,000,000回，test_size=10）

また，学習回数が少ない場合など．（学習回数10,000回，test_size=100）

学習回数が少ないと，与えた点との近似率が落ちています．（点を通っていない）

学習回数が更に少ない場合．（学習回数1,000回，test_size=100）

ここまでくると，もはや点を通ってすらいないです．

このように，ニューラルネットワークでの関数近似は，与えられた点と現在のニューラルネットワークによる出力の誤差が小さくなるように繰り返し学習することで，与えられた点を通る出力を得られるというものです．

というわけで以下リスト．

プログラム内の変数の詳細についてまとめておくと以下のようになっています．

label	details
Input_size	入力の長さ（大きさ）
Input_dim	入力の次元
Hidden_dim	中間層の数
alpha	シグモイド関数の傾き（後述）
eta	学習係数（他サイト参照）
gain	ゲイン（後述）
test_size	学習後の検証用データの量（大きさ）
th	未使用

ここで，ニューラルネットワークに使用しているシグモイド関数は，0~1の範囲の実数しか取ることができないので，このままでは0~1の出力の範囲の点の近似しかできません．

このため，シグモイド関数にgainを掛けることで，擬似的にシグモイド関数の出力範囲を広げています．
たとえば，gain=100の場合では，シグモイド関数の出力範囲は-50~50の実数となります．
一応，負の値も取れるように，負方向にgainの半分だけずらしています．

これらの処理により，ニューラルネットワークで-50~50までの範囲の点を近似することができるようになっています．

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

#define Input_size 4
#define Input_dim 1
#define Hidden_dim 10
#define alpha 10
#define eta 0.8
#define gain 100
#define test_size 100
#define th 0.5

#define INPUT_NAME "data.dat"
#define RHO_NAME "rho.dat"

double a[Input_size][Input_dim] = { 0 }, rho[Input_size] = { 0 };
double v[Hidden_dim][Input_dim + 1] = { 0 }, w[Hidden_dim + 1] = { 0 };
double I[Input_dim + 1] = { 0 }, H[Hidden_dim + 1] = { 0 };
double E[Input_size] = { 0 };
double O = 0;
double O_out[Input_size][Input_dim + 1] = { 0 }, test_out[test_size][Input_dim + 1];
double dk, dj[Hidden_dim] = { 0 };

void Loading();
void Init();
void Display();
void Test();
void Save();

int main(void)
{
    int i, j, k, t;
    double temp, E_max = 0;

    Loading();
    Init();

    Display();

    //main loop
    k = 0;
    while (1){
        k++;

        for (t = 0; t < Input_size; t++){
            for (i = 0; i < Input_dim; i++){
                I[i + 1] = a[t][i];
            }

            for (i = 0; i < Hidden_dim; i++){
                temp = 0;
                for (j = 0; j < Input_dim + 1; j++){
                    temp += I[j] * v[i][j];
                }
                H[i + 1] = 1 / (1 + exp(-alpha*temp));
            }

            temp = 0;
            for (i = 0; i < Hidden_dim + 1; i++){
                temp += H[i] * w[i];
            }

            O = 1 / (1 + exp(-alpha*temp));

            E[t] = 1 / 2.0*pow(gain*(rho[t] - O), 2);

            dk = -(rho[t] - O)*O*(1 - O);

            for (i = 0; i < Hidden_dim; i++){
                dj[i] = dk*w[i + 1] * H[i + 1] * (1 - H[i + 1]);
            }

            for (i = 0; i < Hidden_dim + 1; i++){
                w[i] -= eta*dk*H[i];
            }

            for (i = 0; i < Hidden_dim; i++){
                for (j = 0; j < Input_dim + 1; j++){
                    v[i][j] = v[i][j] - eta*dj[i] * I[j];
                }
            }

            for (i = 0; i < Input_dim + 1; i++){
                if (i == Input_dim){
                    O_out[t][i] = O*gain - gain / 2.0;
                }
                else{
                    O_out[t][i] = a[t][i];
                }
            }
        }
        /*
        E_max=0;
        for (i=0;i<Input_size;i++){
        if(fabs(E[i])>E_max){
        E_max=fabs(E[i]);
        }
        }

        if(E_max<1e-5){
        printf("\nk=%d\nE_max=%f\n",k,E_max);
        for (i=0;i<Input_size;i++){
        printf("%f\n",E[i]);
        }
        printf("\nO_out=\n");
        for (i=0;i<Input_size;i++){
        for (j=0;j<Input_dim+1;j++){
        printf("%f\t",O_out[i][j]);
        }
        printf("\n");
        }
        break;
        }
        */
        if (k >= 1e6){
            printf("\nk=%d\n", k);
            for (i = 0; i < Input_size; i++){
                printf("%f\n", E[i]);
            }
            printf("\nO_out=\n");
            for (i = 0; i < Input_size; i++){
                for (j = 0; j < Input_dim + 1; j++){
                    printf("%f\t", O_out[i][j]);
                }
                printf("\n");
            }
            break;
        }
    }

    Test();
    Save();
    return 0;
}

void Loading(){
    FILE *fp;
    int i, j;

    //Input a
    if (fopen_s(&fp, INPUT_NAME, "r") != 0){
        printf("INPUT FILE open error!\n");
        exit(-1);
    }

    for (i = 0; i < Input_size; i++){
        for (j = 0; j < Input_dim; j++){
            fscanf_s(fp, "%lf ", &a[i][j]);
        }
    }
    fclose(fp);

    //Input rho 
    if (fopen_s(&fp, RHO_NAME, "r") != 0){
        printf("RHO FILE open error!\n");
        exit(-1);
    }

    for (i = 0; i < Input_size; i++){
        fscanf_s(fp, "%lf ", &rho[i]);
    }
    fclose(fp);
}

void Init(){
    int i, j;
    for (i = 0; i < Input_size; i++){
        rho[i] = (gain / 2 + rho[i]) / gain;
    }

    I[0] = 1;
    H[0] = 1;

    srand((unsigned int)time(NULL));

    for (i = 0; i < Hidden_dim; i++){
        for (j = 0; j < Input_dim + 1; j++){
            v[i][j] = 0.2*rand() / (double)RAND_MAX - 0.1;
        }
    }

    for (i = 0; i < Hidden_dim + 1; i++){
        w[i] = 0.2*rand() / (double)RAND_MAX - 0.1;
    }
}

void Display(){
    int i, j;

    printf("a=\n");
    for (i = 0; i < Input_size; i++){
        for (j = 0; j < Input_dim; j++){
            printf("%f", a[i][j]);
        }
        printf("\n");
    }

    printf("\nrho=\n");
    for (i = 0; i < Input_size; i++){
        printf("%f\n", rho[i]);
    }

    printf("\nv=\n");
    for (i = 0; i < Hidden_dim; i++){
        for (j = 0; j < Input_dim + 1; j++){
            printf("%f\t", v[i][j]);
        }
        printf("\n");
    }
}

void Test(){
    int i, j, t;
    double max[Input_dim] = { 0 }, min[Input_dim] = { 0 };
    double p;
    double test[test_size + 1][Input_dim];
    double temp;

    for (i = 0; i < Input_dim; i++){
        max[i] = a[0][i];
        min[i] = a[0][i];
    }

    for (i = 0; i < Input_size; i++){
        for (j = 0; j < Input_dim; j++){
            if (max[j] < a[i][j]){
                max[j] = a[i][j];
            }
            if (min[j] > a[i][j]){
                min[j] = a[i][j];
            }
        }
    }

    printf("max=%f,min=%f\n", max[0], min[0]);

    for (i = 0; i < Input_dim; i++){
        j = 0;
        for (p = min[0]; p <= max[0]; p += (max[0] - min[0]) / test_size){
            test[j][i] = p;
            j++;
        }
    }

    for (t = 0; t < test_size + 1; t++){
        for (i = 0; i < Input_dim; i++){
            I[i + 1] = test[t][i];
        }

        for (i = 0; i < Hidden_dim; i++){
            temp = 0;
            for (j = 0; j < Input_dim + 1; j++){
                temp += I[j] * v[i][j];
            }
            H[i + 1] = 1 / (1 + exp(-alpha*temp));
        }

        temp = 0;
        for (i = 0; i < Hidden_dim + 1; i++){
            temp += H[i] * w[i];
        }

        O = 1 / (1 + exp(-alpha*temp));

        for (i = 0; i < Input_dim + 1; i++){
            if (i == Input_dim){
                test_out[t][i] = O*gain - gain / 2.0;
            }
            else{
                test_out[t][i] = test[t][i];
            }
        }
    }
}

void Save(){
    int i, j;
    FILE *fp;

    if (fopen_s(&fp, "./O_out.dat", "w") != 0){
        printf("O_out.dat open Error!\n");
        exit(-1);
    }

    for (i = 0; i < Input_size; i++){
        for (j = 0; j < Input_dim; j++){
            fprintf(fp, "%f\t%f\n", a[i][j], rho[i] * gain - gain / 2.0);
        }
    }

    fclose(fp);

    if (fopen_s(&fp, "./test_out.dat", "w") != 0){
        exit(-1);
    }

    for (i = 0; i < test_size + 1; i++){
        for (j = 0; j < Input_dim + 1; j++){
            fprintf(fp, "%f\t", test_out[i][j]);
        }
        fprintf(fp, "\n");
    }
    fclose(fp);
}