# include	<iostream.h>
# include	<stdlib.h>
# include	<stdio.h>
# include	<math.h>
# include	<ctype.h>
# include	<time.h>
# define N 10
# define NN N*N
# define G(x) ((1.0+ tanh(x/u0))/2.0)  /* threshold funtion */
void scities();               /* select city position */
void sinit();                  /* select initial neural states */
void cstates() ;               /* calculate neural states */
void dstates();               /* display neural states */
int recheck();           
int m=15;
int tm, aa;
clock_t start0,finish0;//計時
double v[NN],      /* neuron output*/
       v1[14000],
       u[NN],      /* neurons input */
       dd[NN],     /* traverlling distance */
       t[NN],
       xx[N], yy[N], /* coordinate of cities */
       e,           /* energy */
       f,
       sub=0.00001;
double a=0.5,
       b=0.500,
       c=0.200,
       d=0.5,
       u0=0.03,
       h=0.05,
       l[NN],
       pi=3.1415926;
int right=0;
int wrong=0;
double bestd;
FILE *fp,*fopen();
void main()
{
      int i;
      double f1;
      fp=fopen("result.txt","w");
      i=0;
	  bestd=1000;
	  //預置一批S型激勵函數值,以留作后面使用
      f1=0-0.07;
      do
	  {   
		  i++;
          f1+=sub;
          v1[i]=G(f1);
	  } while((v1[i]<=0.999)&&(i<=13999));

	  srand(time(NULL));
	 //產生十個城市的位置分布
	 scities();

	 //產生50組神經網絡
     start0=clock();
 	 for(i=1;i<=100;i++)
	 {   tm=0;
         aa=i*10;
//         cout<<"組號"<<i;         
	
		 sinit();	  //產生神經元的初始狀態
 		 f=0;
         do
		 { //進行神經元狀態的迭代運算
			 cstates();
             //判斷前后兩次的能量函數值是否很接近,若很接近,則結束運算
			 if(fabs(e-f)<1e-20) 
				 break;
 			 //檢查旅行路線的合法性
			 if( recheck())
				 break;
             f=e;
		 }while(tm<1000);  //迭代超過1000次未達到穩定狀態,則求解失敗
 		 //計算在最終穩定狀態下合法旅行路徑的長度,并顯示所有神經元的狀態
		 dstates();
	 }
	 finish0=clock();
	 cout<<(double)(finish0-start0)/CLK_TCK<<"秒"<<endl;
//	 fclose(fp);
	 cout<<"right:"<<right<<",wrong:"<<wrong<<endl;
	 cout<<"bestd:"<<bestd<<endl;
}
void scities()
{   int i=0;
    int j;
    double h[N],o,w,oo;
/*	for(i=0;i<N;i++)
	{
		xx[i]=rand()/(float)32767;
		yy[i]=rand()/(float)32767;
	}*/
	//給出每個城市的坐標
    xx[0]=0.4;   yy[0]=0.4493;
	xx[1]=0.2493;yy[1]=0.1463;
	xx[2]=0.1707;yy[2]=0.2293;
	xx[3]=0.2293;yy[3]=0.7610;
	xx[4]=0.5171;yy[4]=0.9414;
	xx[5]=0.8732;yy[5]=0.6536;
	xx[6]=0.6878;yy[6]=0.5219;
	xx[7]=0.8488;yy[7]=0.3609;
	xx[8]=0.6683;yy[8]=0.2536;
	xx[9]=0.6195;yy[9]=0.2643;

	//計算城市之間的距離
    for(i=0;i<N;i++)
		for(j=0;j<N;j++)
		{
			if(i==j) continue;
            dd[i*N+j]=hypot(xx[i]-xx[j],yy[i]-yy[j]);
        }

     //根據坐標先初始化下神經元的閾值
     for(i=0;i<N;i++)
	 {	 o=(yy[i]-0.5)/(xx[i]-0.5);
         h[i]=atan(o);
         oo=hypot(xx[i]-0.5,yy[i]-0.5);
         for(j=0;j<N;j++)
		 {
			 w=h[i]+(j-1)*2*pi/(float)N;
             l[i*N+j]=cos(w)*oo;
         }
       }
	 /*
	for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			l[i*N+j]=rand()/(float)32767;
*/
	
}
void sinit()//初始化神經元
{
     int i,i1;
     double u00=0-u0*log(N-1)/2.0;
	 //初始化神經元參數
     for(i=0;i<aa;i++) 
	 {
		 t[0]=rand()/(float)32767;
	 }
     for(i=aa;i<aa+NN;i++) 
	 {
		 t[i-aa]=rand()/(float)32767;
	 }
	 //求神經元的初始狀態
     for(i=0;i<NN;i++)
     {
         u[i]=u00+0.001*(t[i]*2-1)+0.002*l[i];
         i1=(int)(u[i]*100000.0+0.5)+7000;
         if(i1 > 13907) v[i]=v1[13907];
         if(i1<=1) v[i]=v1[1];
         if(i1>1 && i1<=13907)
			 v[i]=v1[i1];
      }
}
void cstates()//能量計算
{
	int i1,i,j,x,y,x0,y0;
	double z,k,z1;
    e=0.0;k=0;
    for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			k+=v[i*N+j];
	//求能量函數值
	e=0.0;
	for(x=0;x<N;x++)
	{
		x0=x*N;
		for(i=0;i<N;i++)
		{
			if(i==j) continue;
			e+=v[x0+i]*v[x0+j];
		}
	}
	for(i=0;i<N;i++)
		for(x=0;x<N;x++)
		{
			x0=x* N;
			for(y=0;y<N;y++)
			{
				if (x==y) continue;
				e+=v[x0+i]*v[y*N+i];
			}
		}
	for(x=0;x<N;x++)
	{
		x0=x*N;
        for(y=0;y<N;y++)
		{
			if (y==x) continue;
            y0=y*N;
            for(i=0;i<N;i++)
			{
				if (i==0)
                   e+=v[x0]*dd[x0+y]*(v[y0+1]+v[y0+N-1]);
                else 
					if(i==N-1)
						e+=v[x0+i]*dd[x0+y]*(v[y0+N-2]+v[y0]);
					else
						e+=v[x0+i] * dd[x0+y]*(v[y0+i-1]+v[y0+i+1]);
              }
         }
    }
    e+=(e*a+c*(k-N)*(k-N))/2.0;
    //計算 duxi/dt 
    for(x=0;x<N;x++)
	{
		x0=x*N;
        for(i=0;i<N;i++)
		{
			z=0-c*(k-m);
			for(j=0;j<N;j++)
			{
				if(i==j) continue;
                z-=v[x0+j];
            }
            for(y=0;y<N;y++)
			{
				if(x==y) continue;
				z-=v[y*N+i];
			}
            u[x0+i]+=h*z;
            //再求神經元的狀態
			z1=u[x0+i]* 100000.0+0.5;
            i1=(int) z1 +7000;
            if(i1>13907) v[x0+i]=v1[13907];
            if(i1<=1) v[x0+i]=v1[1];
            if(i1>1 && i1<=13907) v[x0+i]=v1[i1];
		}
	}
	tm+=1;
}
void dstates()
{
     int i, j, x0;
     double dis;
//     fprintf(fp,"iterations=%d   e=%f",tm, e);
//     cout<<"  iterations="<<tm<<"  e="<<e<<"    ";
     if(recheck())
	 {
//		 cout<<" right path \n";
 //        fprintf(fp, "  right path\n");
		 right++;
		//求最終的路徑長度
		dis=0;
		for (i=0;i<N;i++)
		  for (j=0;j<N;j++)
		  {
			  if (v[i*N+j]>0.99)
				  dis=dis+dd[i*N+j];
		  }
 //       fprintf(fp,"distance = %f \n",dis);
//		cout<<"最短距離是"<<dis<<endl;
		if(dis<bestd)
			bestd=dis;
        /* ovput the result of neuron satrix */
        for(i=0;i<N;i++)
		{
			x0=i*N;
//			for(j=0;j<N;j++)
//				fprintf(fp,"%3.0f%s",v[x0+j],",   ");
//			fprintf(fp,"\n");
		}
			fprintf(fp,"\n\n");
	 }
     else
	 {
//		 fprintf(fp, " wrong path \n\n");
//         cout<<"wrong path \n";
		 wrong++;
	 }
}
int recheck()
{
	int i,j,x0;
    double k;
	/*neuron's State must access 0 or 1 */
    for(i=0;i<NN;i++)
		if((v[i]>0.01) && (v[i]<0.99))
			return 0;
	/*every row have and only have one 1 */
	for(i=0;i<N;i++)
	{
		k=0.0;
        x0 =i*N;
        for(j=0;j<N;j++)
			k+=v[x0+j];
        if((k-1.0)>0.1)
			return 0;
	}
    /* every column have and only have one 1*/
    for(i=0;i<N;i++)
	{
		k=0.0;
        for(j=0;j<N;j++)
			k+=v[j*N+i];
        if((k -1.0)>0.1)
            return 0;
   }
   return 1;
}
	/*****end******/