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

14 件のコメント:

Unknown さんのコメント...

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

Hamachi さんのコメント...

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

こちらの記載不足ですが,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です.

Unknown さんのコメント...

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

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ファイルについては少し記述が不親切だったかもしれませんね.
今後も(もし機会があれば)お気軽にどうぞ!

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

質問失礼します。
コンパイルは通りましたが、ファイルの読み込みがうまくいきません。
他の方のコメントから、自分の作成したdatファイルに問題があると思うのですが
datファイルの具体的な書き方について教えていただけるとうれしいです。

Hamachi さんのコメント...

匿名さん
はじめまして.
datファイルの作り方なのですが,まずテキストファイルが"data.dat"と"rho.dat"の2つ必要です.
"data.dat"はx軸に関するファイルで,今回の例ですと
0
1
2
3
のような改行で区切った数値です.もう片方の"rho.dat"はy軸に関するファイルで,先程と同様に
1
3
2
4
のような内容となります.
お使いになっている環境がわかりませんが,環境によってはこれら2つのファイルを置く場所が特殊な場合がありますので(Debugフォルダの中など),
1. 心当たりのある場所に複数コピーして試していただくか,
2. 適当なファイルを書き込みタイプでfopenして場所を確認して頂く
ことが良いかもしれません.

CodeBlocksで実行した環境を以下のリンクにアップロードしましたので,一部でも参考になれば幸いです.
記事中のコードはVisualStudioで実行したときのもので,以下のリンクとはfopen_sなどのファイル処理部分が少し変わっています.

1週間で消えますので,消えてしまった場合はご一報ください.
http://fast-uploader.com/file/7096377823868/

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

ありがとうございます。
おかげさまで正常に動作しました。
どうやらファイルの置き場所が悪かったようです。
今後もわからないことがあれば質問させていただくつもりです。
お手数をおかけしますがまたよろしくお願いいたします。

Hamachi さんのコメント...

匿名さん
解決されましたか!よかったです.
ファイルの置き場所は自分で作っていても間違うことが多いので難しいですね.
あまりお役に立てることは多くありませんが,お気軽にご質問ください.

ブログは移転しまして今後は(https://hamachannel.hatenablog.com/)が更新されます.なお,こちらを削除することはありませんのでご安心ください(?)

コメントを投稿