付録A 学習プログラム（ニューラルネットワーク）

/* 学習プログラム（ニューラルネットワーク） : main.c */

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

#define Iunits  20            /* 各ユニットの個数 */
#define Hunits  50
#define Ounits  2
#define SMAX    256

#define Rlow         -0.30
#define Rhigh         0.30
#define A             1.4
#define B             0.3

#define fout(x)   ( 1 / ( 1 + exp( -(x) ) ))                /* シグモイド関数 */
#define urand()   ( drand48() * (Rhigh - Rlow) + Rlow )

double wih[Hunits][Iunits];   /* 結合荷重(重み) */
double who[Ounits][Hunits];

double dwih[Hunits][Iunits];  /*  重みの変化量  */
double dwho[Ounits][Hunits];

double shikiih[Hunits];       /*    しきい値    */
double shikiio[Ounits];

double dshikiih[Hunits];      /* しきい値の変化量 */
double dshikiio[Ounits];

double iout[Iunits];          /*    各層の出力  */ 
double hout[Hunits];
double oout[Ounits];

double tsignal[Ounits];       /*     教師信号    */

double drand48();
double Eta,Alpha;
double ts[SMAX];              /* musicデータを配列 */
int    tsmax;

FILE *fpmusic;                /* musicファイルのポインタ */

  void foward_propagation();
  void back_propagation();
  void initialize();

main( int argc , char *argv[] )
{
  FILE  *fp,*fp2;
  int  ct;
  long kaisuu;
  char music[SMAX],name[SMAX];
  unsigned int dn,i,j,k,kk;
  double averror,dummy,error,serror;

  if( argc !=4 ){
    printf(" Usage : Eta Alpha filename \n");
    exit(1);
  }

  Eta = atof(argv[1]);
  Alpha = atof(argv[2]);

  kk = 150000;              /* 学習回数 */

  /* musicファイルの指定 */
  strcpy(music,argv[3]);

  /* ニューロンの初期化 */
  initialize();

  /* musicファイルのオープン */
  fpmusic = fopen(music,"r");
  if( fpmusic == NULL ){
    printf("Can't open the file \n");
    exit(1);
  }

  /* データの読み込みとデータ数のカウント */
  ct = 0;
  while( fscanf(fpmusic,"%lf %lf",&ts[ct],&ts[ct+1]) == 2 )
    ct += 2;

  /* 教師信号数を記憶 */
  tsmax = ct/2;

  /* musicファイルのクローズ */
  fclose(fpmusic);

  error = serror = 0;

  /*学習中における誤差変化ファイルのオープン*/
  sprintf(name,"Berror_E%0.4f_A%0.2f.data",Eta,Alpha);
  fp = fopen(name,"w");
  if (fp == NULL){
    fprintf(stderr,"error\n");
    exit(1);
  }

  for( k = 1 ; k <= kk ; k++ ){
    /* 初期値入力 */
    for( i = 0, j = tsmax * 2 - 20 ; i < 20 ; i++, j++ ){
      iout[i] = ts[j];
    }
    for( kaisuu = 0 ; kaisuu < tsmax ; kaisuu++ ){

      tsignal[0] = ts[kaisuu*2];
      tsignal[1] = ts[kaisuu*2+1];

      foward_propagation();

      back_propagation();

      for( i = 0 ; i < 18 ; i++ ){
        iout[i] = iout[i+2];
      }
      iout[18] = tsignal[0];
      iout[19] = tsignal[1];

    }

    /* 最期に各ニューロンの重みの記録をする */
    if( k == kk ){
      dn=k/100;
      sprintf(name,"b%d_E%0.4f_A%0.2f.data",dn,Eta,Alpha);
      fp2 = fopen(name,"w");
      if (fp2 == NULL){
        fprintf(stderr,"error\n");
        exit(1);
      }
      /*最終的な各ユニット、その他のデータをファイルに書き込む*/
      fwrite(wih,sizeof(double),Iunits * Hunits,fp2);
      fwrite(who,sizeof(double),Hunits * Ounits,fp2);
      fwrite(shikiih,sizeof(double),Hunits,fp2);
      fwrite(shikiio,sizeof(double),Ounits,fp2);
      fclose(fp2);
    }

    /*50回ごとにネットワークに予測をさせ誤差を記録*/
    if( k%50  == 0 ){
      for( i = 0, j = tsmax * 2 - 20 ; i < 20 ; i++, j++ ){
        iout[i] = ts[j];
      }
      serror = 0;
      for( kaisuu = 0 ; kaisuu < tsmax ; kaisuu++ ){
        tsignal[0] = ts[kaisuu*2];
        tsignal[1] = ts[kaisuu*2+1];
        foward_propagation();
        error = 0;
        for ( i = 0 ; i < Ounits ; i++ ){
          dummy = tsignal[i] - oout[i];
          error += (dummy * dummy);
        }
        error = sqrt(error);
        serror += error;
        for( i = 0 ; i < 18 ; i++ ){
          iout[i] = iout[i+2];
        }
        iout[18] = tsignal[0];
        iout[19] = tsignal[1];
      }
      averror = serror/tsmax;
      fprintf(fp,"%d\t%f\n",k,averror);
      printf("%d\t%f\n",k,averror);
    }
  }
  fclose(fp);
  putchar('G'-'@');
}

/*各ユニットに乱数で重みを振り当てる*/
void initialize()
{
  int i,j;

  for ( i = 0 ; i < Hunits ; i++ ){
    for ( j = 0 ; j < Iunits ; j++ ){
      wih[i][j] = urand();
    }
    shikiih[i] = urand();
  }

  for ( i = 0 ; i < Ounits ; i++ ){
    for ( j = 0; j < Hunits; j++ ){
      who[i][j] = urand();
    }
    shikiio[i] = urand();
  }
}

/* フォワードプロパゲーション */
void foward_propagation()
{
  int i,j;
  double sum;

  for( i = 0 ; i < Hunits ; i++ ){
    sum = 0;
    for( j = 0 ; j < Iunits ; j++ ){
      sum += wih[i][j] * iout[j];
    }
    hout[i] = fout(sum + shikiih[i] );
  }
  for ( i = 0; i < Ounits ; i++ ){
    sum = 0;
    for ( j = 0; j < Hunits ; j++ ){
      sum += who[i][j] * hout[j];
    }
    oout[i] = fout(sum + shikiio[i]);
  }
}

/* バックプロパゲーション */
void back_propagation()
{
  int    i,j;
  double dih[Hunits],dho[Ounits];           /* 学習信号 */
  double sum;

  for ( i = 0; i < Ounits; i++ ){
    /*出力の伝達関数はy=xという形になっているので伝達関数の微分は1ということになる*/
    /*出力層の学習信号の計算*/
    dho[i] = oout[i] * ( 1 - oout[i] ) * ( tsignal[i] - oout[i] );

  }
  /*出力の重み補正*/
  for ( i = 0; i < Hunits ; i++ ){
    for ( sum = 0, j = 0 ; j < Ounits ; j++ ){
      dwho[j][i] = ( Eta * dho[j] * hout[i] ) + ( Alpha * dwho[j][i] );
      sum += ( dho[j] * who[j][i] );
      who[j][i] += dwho[j][i];
      
    }
    /*中間層の教師信号の計算*/
    dih[i] = hout[i] * ( 1 - hout[i] ) * sum;
  }
  /*出力のしきい値補正*/
  for ( i = 0 ; i < Ounits ; i++ ){
    dshikiio[i] = ( Eta * dho[i] ) + ( Alpha * dshikiio[i] );
    shikiio[i] += dshikiio[i];
  }
  /*中間層の重み補正*/
  for ( i = 0 ; i < Iunits ; i++ ){
    for ( j = 0 ; j < Hunits ; j++ ){
      dwih[j][i] = ( Eta * dih[j] * iout[i] ) + ( Alpha * dwih[j][i] );
      wih[j][i] += dwih[j][i];
    }
  }
  /*中間層のしきい値補正*/
  for ( i = 0 ; i < Hunits ; i++ ){
    dshikiih[i] = ( Eta * dih[i] ) + ( Alpha * dshikiih[i] );
    shikiih[i] += dshikiih[i];
  }
}