2015年8月1日土曜日

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

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

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

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

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

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

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

この2種類のファイルから,Input_sizeInput_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

approx

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

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

test_sizeを小さくするとこんな感じに.(学習回数1,000,000回,test_size=10

approx2

また,学習回数が少ない場合など.(学習回数10,000回,test_size=100

approx3

学習回数が少ないと,与えた点との近似率が落ちています.(点を通っていない)

学習回数が更に少ない場合.(学習回数1,000回,test_size=100

approx3

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

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

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

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

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);
}

10 件のコメント:

Unknown さんのコメント...

こんばんは、最近ニューラルネットワークに興味を持ちここに行き着いたものです。
上記のソースを拝借し、gcc neural.c -lm でコンパイルにかけたのですが、fopen_s,fscan_sについて.text0x**定義されていない参照です、とエラーが出て先に進めません。原因は何でしょうか?

hamachan さんのコメント...

こんばんは.
こんなヘンピなところに興味を持ってくださりありがとうございます(笑)

こちらの記載不足ですが,fopen_sとfscanf_sはVisual studio用コンパイラの記載となっています.(うろ覚えですが)

gccでコンパイルされる際には,
if(fopen_s(&fp, "FILENAME", "w") != 0)
となっている部分を,
if((fp = fopen("FILENAME", "W")) == NULL)
のように書き換え,さらに
fscanf_sをfscanfに書き換えていただければコンパイルが通るのではないかと思います.

注意すべきは,アドレスとして渡していたfpが関数の外に出ていることと,
fp = fopen(~)の部分全体がカッコで囲まれていること,
比較が != 0ではなく == NULLになっているところあたりかと思います.

1年以上前の汚らしいコードですが,こちらの環境では上記で動作確認できましたので,何かあれば気軽にご連絡くださいませ

Hamachi さんのコメント...

上記でfopenの中が大文字のWになっていますが,正しくは小文字のwです.

Naoki Odajima さんのコメント...

先ほど質問させていただいたものです。
早いご返答ありがとうございます!おかげで解決できました!

Hamachi さんのコメント...

わかりづらいところばかりですが何かのご一助になれば幸いです(^-^)

NN さんのコメント...

いつも参考にさせて頂いています。
ありがとうございます。
void Test()以降の処理時間はHamachiさんのPCでどれくらいですか?
コードそのまま実行すると1時間以上結果が出ずどこかマズいのか困ってます。
1e6を小さくしても変わりません。
アドバイスお願いします!

Hamachi さんのコメント...

NNさん
こんばんは.いつもありがとうございます.

void TEST()以降の計算ですが,この関数が呼ばれるのは学習が終わった後なので,
この先の処理には殆ど時間はかからないように思います.

試しに自分の環境で実行してみたところ,学習も含めた全体の実行時間が3.218秒でした.

実行環境は以下のとおりです.(システムプロパティより)

OS: Windows7 Home Premium (64 bit)
CPU: Intel(R) Core(TM)i7-4710MQ CPU @ 2.50GHz 2.50Ghz
メモリ: 8.00 GB

1時間も貴重な時間を棒に振ってしまい申し訳ない限りです……
もっとわかりやすいコードを書かねば……
実行環境等,記載していただけるとヒント程度なら得られるかもしれません.

参考までに,コンソールで出力される実行結果は次のような感じです.
グラフについては,実行後に作成される"test_out.dat"をExcelに貼り付けて出しています.

----
a=
0.000000
1.000000
2.000000
3.000000

rho=
0.510000
0.530000
0.520000
0.540000

v=
0.014035 -0.079638
-0.062505 0.093683
0.054643 0.026487
0.016984 0.016752
-0.064397 0.028739
0.087671 -0.066326
-0.086114 -0.041295
-0.052513 0.072643
-0.061016 -0.034367
0.055437 0.036454

k=1000000
0.000000
0.000000
0.000000
0.000000

O_out=
0.000000 1.000000
1.000000 3.000000
2.000000 2.000000
3.000000 4.000000
max=3.000000,min=0.000000

Process returned 0 (0x0) execution time : 3.218 s
Press any key to continue.

Hamachi さんのコメント...

蛇足ですが,今回の実行環境はCode Blocks 16.01です(今インストールしたので特別な設定はしていないです).
コード作成時はVisual Studioだったように思います.

匿名 さんのコメント...

こんなにすぐにコメントありがとうございます!
確認したら、a[i][j]が入ってませんでした。自分のdatファイルの作り方がまずかったようです。本当に申し訳ありません、修正後に実行したら瞬足で終了しました!
よくよく確認してから書き込みます。ありがとうございました!

Hamachi さんのコメント...

動きましたか!良かったです.
いま見返してみるとdatファイルについては少し記述が不親切だったかもしれませんね.
今後も(もし機会があれば)お気軽にどうぞ!

コメントを投稿