/* Program calculates average attraction radius (H0) as a function of
   M(number of patterns in a learning set) and D (Cii=D*Cii) using 
 flood-fill neuroproccessing algorithm */

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

#define MAXLONG -((1<<31)+1)  /* For random number generator */
#define MaxN 10000

#define N 200      /* The size of a net            */

#define ALLOWABLE_H  N*0.001 /* (i.e. only 1% of patterns may be retrieved
			       incompletely (i.e. with H=1) */
#define Nimpl 10  /* No of Noise Implementaions  for each amount H0*/
#define Nsets 5   /* No of sets of prototypes used */
#define SETSTART 0/* The learning set to start     */
int dH=N/100;


		       
int  M, m, i, j, it;


int V[N], Y[N], 
    ebuf[N], eN, eNlast, e, Pcyc;
float S[N],S0[N], C[N][N], E;
int H0, H0min, H0max;
float Have, cyc, it_ave;

/* for findC and testC */
  char  pattern[N+1], filename[20], s[4];
  FILE *f;
/* for iterate() */
  int  i1, i2, lasti1, lasti2; 
/* for testC() */
  int initH0, H;
  
/* for main */
  int init;
  float D;
  float HH0,HHave,Ccyc,Iit_all;
  FILE *fout;
  char fname[30];


void findC()
{
  float X,Y,Z;

  f=fopen(filename, "r");
  fgets(pattern, MaxN, f);  /* first comment line */

  for(i=0;i<N;i++)
    for(j=0;j<N;j++)
      C[i][j]=0;

  for(m=0;m<M;m++)
    {
      fgets(pattern, MaxN, f); 
      for(i=0;i<N;i++)
	V[i]='*' - 1 - pattern[i];
      
      for(E=N,i=0;i<N;i++)
	{
	  for(S[i]=0,j=0;j<N;j++)
	    S[i]+=C[j][i]*V[j];
	  E-=S[i]*V[i];
	}
      for(i=0;i<N;i++)
	for(j=0;j<N;j++)
	  C[i][j]+=(V[i]-S[i])*(V[j]-S[j])/E;
    }

  for(X=0,Y=0,Z=0, i=0;i<N;i++)
    {
      for(j=0;j<N;j++)
	{
	  Y+=C[i][j]*C[i][j];  Z+=fabs(C[i][j]);
	}
      X+=C[i][i];
    }
  
#if 0  /* to output average values of weights Cii and Cij */
  X=X/N; Y=Y/N/N; Z=(Z/N/N)*(Z/N/N); 
  printf("\nC:  <Cii> = %.5f,\t   <Cij^> = %.5f  vs  <|Cij|>^=%.5f\n", X,Y,Z);
  X=(float)M/N; Y=(float)M*(N-M)/N/N/N;  
  printf("vs   M/N  = %.5f,\tM(N-M)/N3 = %.5f \n", X,Y);
#endif

  fclose(f);
}     

void invert(int Noise)
{
  int kk;

  for(i=0;i<N;i++)  
    Y[i]=-1;

  for(i=0;i<Noise; )  
    {  
      kk=random() % N;
      if (Y[kk] != 0)
	{
	  Y[kk]=0;
	  i++;
	}
    }
  for(i=0;i<N;i++)  
    {
    if (Y[i] == 0)
      Y[i]=-V[i];
    else
      Y[i]=V[i];
    }
}

int iterate()
{
  for(it=0;  ;it++)
    {
      for(e=0;e<eN;e++)
	if (ebuf[e]<0)
	  for(j=0;j<N;j++)
	    S[j]-=2*C[-ebuf[e]-1][j];
	else
	  for(j=0;j<N;j++)
	      S[j]+=2*C[ebuf[e]-1][j];
      
      eNlast=eN;
      eN=0;Pcyc=0;
      for(i=0;i<N;i++)
	{
	  if (S[i]<0)
	    if (Y[i]>0)
	    { 
	      Y[i]=-1; 
	      if (ebuf[eN]==i+1)  /* for check of a cycle */
		Pcyc++;
	      ebuf[eN]=-i-1;eN++;
	    }
	    else ;
	  else 
	    if (Y[i]<0)
	    {
	      Y[i]=+1; 
	      if (ebuf[eN]==-i-1)  /* for check of a cycle */
		Pcyc++;
	      ebuf[eN]=+i+1;eN++;
	    }
	}
      /* printf("-%i", eN); */
      if (eN==0)
	return (0);
      
      if ((eNlast=eN)&&(Pcyc==eN)) /* check  if it is a cycle */ 
	{
	  /*  printf("%i ", eN); */
	  return (-eN);
	}	
    } /* end of iterations */
}

void testC()
{
  for(Have=0,cyc=0, H0=dH; Have<=ALLOWABLE_H; H0+=dH)  
    {
      f=fopen(filename, "r");
      fgets(pattern, MaxN, f);  /* skip first line, as  it's comments*/

      Have=0; cyc=0; it_ave=0; 
      for(m=0;m<M;m++)
	{
	  fgets(pattern, MaxN, f); 
	  for(i=0;i<N;i++)
	    V[i]='*' - 1 - pattern[i];
      
	  for (initH0=0; initH0<Nimpl; initH0++)
	    {
	      srandom( (long)initH0 );   /* e.g. 23th one for all
					    vectors */
	      invert(H0); 
	      
	      for(eN=0,i=0;i<N;i++)
		if (Y[i]<0) 
		  {
		    ebuf[eN]=-i-1; eN++; 
		  }
	      /* printf(" m=%i:eN=%i", m, eN); */

	      for(i=0;i<N;i++)
		S[i]=S0[i];

	      if ( iterate() < 0 ) /* i.e it was a cycle */
		cyc++;

	      for(H=0,i=0;i<N;i++) 
		if (Y[i]!=V[i])
		  H++;

	      Have += H;  it_ave += it;
	    }	 
	}
      Have = Have/M/Nimpl; cyc=cyc/Nimpl; it_ave=it_ave/M/Nimpl;
      fclose(f);
    }
  
}
        
main()
{
  printf("\nN=%i (Allowable <Hfinal>=%.2f ).   Sets: %i+%i, #noises %i", N, ALLOWABLE_H, SETSTART, Nsets, Nimpl ); 

  for (M=(int)(0.6*N);M>=(int)(0.4*N);M-=0.2*N) /*(int)(0.2*N)*/
    {
      printf("\nM=%i", M);
      printf("\n D \t H0max \t it(H0max) \tH0min \tH0max  H0max+1 gives( H ; cyc )   \n");
      HH0=22;
      for (D=0.00;D<1.01;D+=0.05)
	{
	  if ( (HH0<=1+dH)&&(D>0.1) ) break; /* to finish earlier */
	  HH0=0.0; Iit_all=0.0; HHave=0.0; Ccyc=0.0; H0max=0;H0min=N;
	  for (init=0;init<Nsets;init++)
	    {     
	      strcpy(filename, "set100.");   
	      sprintf(s, "%i", init+SETSTART); strcat(filename, s); 
	      
	      findC();

	      for(i=0;i<N;i++)
		C[i][i]=D*C[i][i];

	      for(i=0;i<N;i++)
		for(S0[i]=0,j=0;j<N;j++)
		  S0[i]+=C[j][i];
    
	      testC(); H0=H0-1-dH; /* 1) as "for" adds +1, 2) as last
				   iteration is when H0=Hattr+dH */
	      HH0+=H0; Iit_all+=it_ave; HHave+=Have; Ccyc+=cyc; 
	      if (H0>H0max)
		H0max=H0;
	      else if (H0<H0min)
		H0min=H0;
	    }

	  HH0=HH0/Nsets; Iit_all=Iit_all/Nsets; 
	  HHave=HHave/Nsets; Ccyc=Ccyc/Nsets;	      
     
	  printf("\n %.2f \t %.3f \t  %.3f  \t%i \t%i \t( %.2f; %.2f)",
		 D, HH0, Iit_all, H0min, H0max, HHave, Ccyc );     
	} /* D */
    } /* M */
}