付 録B Ctype気温予測学習プログラムリスト

/*一日の気温だけを入力、出力１日分とゴミ六日分 c_type*/
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <local/bgi.h>

#define Iunits    14              /*    各ユニットの個数    */
#define Hunits    200
#define Ounits    14
#define Tunits    2

#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 hkeep[Hunits];             /*    houtをとっておく  */
 
double drand48();
double Eta,Alpha;
double data[5000];               /*気温データ配列*/

FILE *fpkion; /*気温、湿度ファイルポインタ*/

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


main(argc,argv)
     int   argc;
     char *argv[];    /*メインへのいん数引渡し 実際にはkaisuumaxを決める*/

{
  FILE  *fp,*fp2;
  int gd,gm,cc,ave,j;
  long   kaisuu;
  long   kaisuumax;
  unsigned int   i,dn,k,kk;
  char  name[64],kion[64],a[64],graph[64];
  double averror,dummy,error,serror,max_syokichi,mim_syokichi,x;
  Eta=atof(argv[1]);
  Alpha=atof(argv[2]);

   if (argc != 4){
    printf(" Usage : kion \n" );
    exit(-1);
  }

  kaisuumax=3;      /*学習年数*/
  kaisuumax =kaisuumax*365+1;

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

/*気温データファイルの指定*/
  strcpy(kion,argv[3]);

  error = serror = 0;

/*気温データファイルのオープン*/
  fpkion   = fopen(kion,"r");
  if (fpkion==NULL){
    printf("can't open the file\n");
    exit(1);
  }

/*気温データを配列に代入*/
  for (i=0 ; i<kaisuumax*2 ; i++){
    fscanf(fpkion,"%lf\n",&x);
    data[i]=x;
  }
  fclose(fpkion);

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

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

for (k=1; k<=kk ; k++){
/*初期値入力*/
  for (i=0 ; i<Tunits ; i++){
    iout[i] = data[i];
    }

/*一定の学習回数で定数を変化*/
if(k==3000)Eta=0.00004;
if(k==8000)Eta=0.00002;
if(k==16200)Eta=0.00001;

/*六日分空回りさせる*/
for ( kaisuu = 1 ; kaisuu < 7 ; kaisuu++ ){
  tsignal[0]=data[kaisuu*2];
  tsignal[1]=data[kaisuu*2+1];
  foward_propagation();
  for (i =Tunits ; i < Iunits; i++){
    iout[i] = oout[i];
  }
  iout[0] = tsignal[0];
  iout[1] = tsignal[1]; 
}

for ( kaisuu = 7 ; kaisuu < kaisuumax ; kaisuu++ ){

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

  foward_propagation();
  back_propagation();

  for(i=0;i<Hunits;i++){
  hkeep[i]=hout[i];
  }
  
  for (i =Tunits ; i < Iunits; i++){
  iout[i] = oout[i];
  }
  iout[0] = tsignal[0];
  iout[1] = tsignal[1];

}

/*500回ごとに各ニューロンの重みの記録をする*/
if (k%500 == 0){
  dn=k/100;
  sprintf(name,"c%d_E%0.5f_A%0.1f.data",dn,Eta,Alpha);
  printf("500data OK\n");

  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 ; i<Tunits ; i++){
    iout[i] = data[i];
    }
serror=0;
/*六日分空回りさせる*/
for ( kaisuu = 1 ; kaisuu < 7 ; kaisuu++ ){
  tsignal[0]=data[kaisuu*2];
  tsignal[1]=data[kaisuu*2+1];

  foward_propagation();

  for (i =Tunits ; i < Iunits; i++){
    iout[i] = oout[i];
  }
  iout[0] = tsignal[0];
  iout[1] = tsignal[1]; 
}
for ( kaisuu = 7 ; kaisuu < kaisuumax ; kaisuu++ ){
  tsignal[0]=data[kaisuu*2];
  tsignal[1]=data[kaisuu*2+1];

  foward_propagation();  

  error=0;
  for (i = 0; i < Tunits; i++){
    dummy = tsignal[i] - oout[i];
    error += (dummy * dummy);
  }
  error  =sqrt(error);
  serror += error;  
  for (i =Tunits ; i < Iunits; i++){
    iout[i] = oout[i];
  }
  iout[0] = tsignal[0];
  iout[1] = tsignal[1];
}
averror=serror/kaisuumax;
printf( "%d\tavERROR = %5.10f\n",k,averror);
fprintf(fp,"%d\t%f\n",k,averror);
}

}
  fclose(fp);
  printf("\nkaisuu = %d\n",k);
  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];
    }
    /*出力は0~1ではないので伝達関数にはかけない*/
    oout[i] = sum + shikiio[i];
  }
}

void back_propagation()
{
  int    i,j;
  double dih[Hunits],dho[Ounits];           /* 学習信号 */
  double sum;

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

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