/***********************************************************/
/*         Computer simulation of on-line learninig        */
/*                   (perceptron traininig)                */
/*                                    J. Inoue             */
/***********************************************************/
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#define size 5000           /* 入力の次元 */
#define a 6.0              /* 教師機械の内部パラメータ */ 
#define SEEDINITIAL 341
#define SEEDINITIAL_S 953
float J0[size];            /* 教師機械の結合 */
float J[size];             /* 生徒機械の結合*/
float x[size];             /* 例題 */
int OutTeach;              /* 教師機械の出力 */
int OutStudent;            /* 生徒機械の出力 */
int Delta;             /* 生徒がアップデートするか否かを決めるファクター */
float hT;              /* 教師機械の内部ポテンシャル */
float hS;              /* 生徒機械の内部ポテンシャル */ 
float Norm;            /* 生徒機械の結合の長さの規格化因子 */
float L;               /* 秩序変数（結合の長さ) */
float Q;               /* 秩序変数 (重なり) */  
float ll; 
/**********************************************************************/
/*    Generating randomnumbers by Numerical Recepies in C (ran3)      */
/**********************************************************************/
#define MBIG 1000000000
#define MSEED 161803398
#define MZ 0
#define FAC (1.0/MBIG)
float ran3(long *idum)
{
	static int inext,inextp;
	static long ma[56];
	static int iff=0;
	long mj,mk;
	int i,ii,k;

	if (*idum < 0 || iff == 0) {
		iff=1;
		mj=MSEED-(*idum < 0 ? -*idum : *idum);
		mj %= MBIG;
		ma[55]=mj;
		mk=1;
		for (i=1;i<=54;i++) {
			ii=(21*i) % 55;
			ma[ii]=mk;
			mk=mj-mk;
			if (mk < MZ) mk += MBIG;
			mj=ma[ii];
		}
		for (k=1;k<=4;k++)
			for (i=1;i<=55;i++) {
				ma[i] -= ma[1+(i+30) % 55];
				if (ma[i] < MZ) ma[i] += MBIG;
			}
		inext=0;
		inextp=31;
		*idum=1;
	}
	if (++inext == 56) inext=1;
	if (++inextp == 56) inextp=1;
	mj=ma[inext]-ma[inextp];
	if (mj < MZ) mj += MBIG;
	ma[inext]=mj;
	return mj*FAC;
}
#undef MBIG
#undef MSEED
#undef MZ
#undef FAC
/* (C) Copr. 1986-92 Numerical Recipes Software 1+5-5i. */
/************************************************************/
/*              Generating random examples                  */
/************************************************************/
void GenerateRandomPattern(int k)
{
  long *idum,m;
  int i;
  idum = &m;
  *idum = k;
    for(i=0; i < size; i++){
      x[i] =  2.0*ran3(&m)-1.0;
    }
    for(i=0, Norm=0; i < size; i++){
      Norm = Norm + x[i]*x[i];
    }    
    for(i=0; i < size; i++){
      x[i] = x[i]/sqrt(Norm);
    }
}
/*********************************************************/
/*       Set initial Teacher's weight vector             */
/*********************************************************/
void SetInitial()
{
  long *idum,m; 
  int i,j,k;
  float J0norm;
  idum = &m;
  *idum = -SEEDINITIAL;
  for(i=0; i < size; i++){
    J0[i] = 2.0*ran3(&m) -1.0;
  }
  for(j=0, J0norm=0; j < size; j++){
    J0norm = J0norm + J0[j]*J0[j];
  }
  for(k=0; k < size; k++){
    J0[k] = J0[k]*sqrt(size)/sqrt(J0norm);
  }
}
/******************************************************/
/*               Set initial Student                  */
/******************************************************/
void SetInitialStudent(float d, float r)
{
  int i,j,k;
  float Jnorm;
  long *idum,*idum,m,m2;
  idum = &m;
  idum = &m2;
  *idum = -SEEDINITIAL_S;
  *idum2 = -SEEDINITIAL;
  for(i=0; i < (int)(r*size); i++){
    J[i] = 2.0*ran3(&m2) -1.0;
  }
  
  for(i=(int)(r*size); i < size; i++){
     J[i] = 2.0*ran3(&m) -1.0;
     }
  
  for(j=0, Jnorm=0; j < size; j++){
    Jnorm = Jnorm + J[j]*J[j];
  }
  for(k=0; k < size; k++){
    J[k] = J[k]*(d)*sqrt(size)/sqrt(Jnorm);
  }
}
/******************************************************/
/*           Calculating local fields                */
/*****************************************************/
/*------ 教師の内部ポテンシャルを計算　 -------------*/
void LocalFieldT()
{
  int i;
    for(i=0,hT=0; i < size; i++){
      hT = hT+ J0[i]*x[i];
    }
}
/*----- 生徒の内部ポテンシャルを計算 ----------*/
void LocalFieldS()
{
  int i;
    for(i=0,hS=0; i < size; i++){
      hS = hS + J[i]*x[i]*sqrt(ll)/sqrt(size);
    }
}
/*********************************************************/
/*         Calculating Teacher's OutPut                  */
/*********************************************************/
void TeacherOut()
{
       if((hT > 0 && hT  < a) || (hT < -a)){
      OutTeach = 1;
    }else{
      OutTeach= -1;
    }
}
/*******************************************************/
/*      Calculating Student's OutPut                   */
/*******************************************************/
void StudentOut()
{   
    if(hS > 0){
      OutStudent = 1;
    }else{
      OutStudent = -1;
    }
}
/*************************************************************/
/*    Condition for updating the weight of the Student       */
/*************************************************************/
void CalDelta()
{
    if(OutStudent*OutTeach < 0){
      Delta = 1;
    }else{
      Delta = 0;
    }
}
/*****************************************************/
/*     Update rule of the Student's weight           */
/*****************************************************/
void Update()
{
  int i;
    for(i=0; i < size; i++){
      J[i] = J[i] - Delta*OutStudent*x[i];
    }
}
/**************************************************************/
/*    Calculating the length of the Student's weight vector   */
/**************************************************************/
void CalLength()
{
  int i;  
  for(i=0,ll=0; i < size; i++){
    ll = ll + J[i]*J[i];
  }
  L = sqrt(ll)/sqrt(size);
}
/******************************************************************/
/*  Calculating the overlap bewteen the Teacher and the Student  */
/******************************************************************/
void CalOverlap()
{
  int i;
  for(i=0,Q=0; i < size; i++){
    Q = Q + (J0[i]*J[i]/sqrt(ll))/sqrt(size);
  }
}
/********************************************************/
/*                   Main program                       */
/*******************************************************/
main()
{
  FILE *pf;
  int i,imax,c;
  float d,r;
  imax = 200000;
  d = 2.00;  /* 生徒機械の長さの初期値設定 */
  r = 0.90; /* 生徒機械と教師機械の結合の重なりの初期値設定 */
  
  SetInitial();  
  SetInitialStudent(d,r); 
  CalLength();
  CalOverlap();
  if((pf = fopen("on-b.dat","wt")) != NULL){
    fprintf(pf,"\n %f %f", L,Q);
  }
  fclose(pf);
  if((pf = fopen("on-b.dat","wt")) != NULL){
    for(i=0; i < imax; i++){
      GenerateRandomPattern(-i);
      LocalFieldT();
      LocalFieldS();
      TeacherOut();
      StudentOut();
      CalDelta();
      Update();
      CalLength();
      CalOverlap();
      /* 出力ファイルが大きくなりすぎるので, 5000ステップ毎にモニタリング */
      /* 全ての点をプロットしたければ, fmod(i,5000)をfmod(i,1)に修正する */
      if(fmod(i,5000)==0){fprintf(pf,"\n %f %f",L,Q);}
    }
  }
  fclose(pf);
}