#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#define N 10
#define NN N*N
#define G(x) ((1.0+tanh(x/u0))/2.0)  //threshold function
void scities();  //select city position
void sinit();  //select initial neural states
void cstates();  //caculate neural states
void dstates();  //display neural states
int recheck();
int m=15;
int tm,aa;
float v[100],v1[14000],u[100],dd[100],t[100],xx[10],yy[10],e,f,sub=0.00001;
double a=0.5,b=0.5,c=0.2,d=0.5,u0=0.02,h=0.01,l[100],pi=3.1415926;
FILE *fp,*fopen();

main()
{
  int i,j,in,peng;
  float f1;
  fp=fopen("result.dat","w");
  i=0;
  f1=0-0.07;
  do{
    i++;
    f1+=sub;
    v1[i]=G(f1);
  }while((v1[i]<=0.999) && (i<=13999));
  scities();
  for(i=1;i<=50;i++)
  {
    tm=0;
    aa=i*10;
    printf("%d",i);
    sinit();
    f=0;
    do
    {
      cstates();
      if(fabs(e-f)<1e-20)
        break;
      f=e;
    }while(tm<1000);
    dstates();
  }
}

void scities()
{
  int i,j;
  double h[N],o,w,oo;
  //get the coordinate of cities using random data
  //cites coordinates given by Hopfield-Tank
  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]);
    }
  }
  //caculate initial bias
  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;
    }
  }
}

void sinit()
{
  int i,j,i1;
  float u00=0-u0*log(N-1)/2.0;
  //get initial neuron's state
  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>13908) v[i]=v1[13908];
    if(i1<=1) v[i]=v1[1];
    if(i1>1 && i1<=13908)
      v[i]=v1[i1];
  }
}

void cstates()
{
  int i1,i,j,q,x,r,y,x0,y0,z0;
  float z,k,e1,z1;
  e=0;
  k=0;
  for(i=0;i<N;i++)
    for(j=0;j<N;j++)
      k+=v[i*N+j];
    //caculate energy function
    e=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;
      //caculate 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+j];
          }
          u[x0+i]+=h*z;
          z1=u[x0+i]*100000.0+0.5;
          i1=(int)z1+7000;
          if(i1>13908) v[x0+i]=v1[13908];
          if(i1<=1) v[x0+i]=v1[1];
          if(i1>1 && i1<=13908) v[x0+i]=v1[i1];
        }
      }
      tm+=1;
}

void dstates()
{
  int i,j,x0;
  float dis;
  fprintf(fp,"iterations=%d e=%f",tm,e);
  printf("iterations=%d e=%f",tm,e);
  if(recheck())
  {
    printf("right path\n");
    fprintf(fp,"right path\n");
    for(j=0;j<N;j++)
    {
      if(v[j*N]>=0.99)
      {
        x0=j;
        break;
      }
    }
    dis=0;
    for(i=0;i<N;i++)
    {
      //?????????????/
      if(i==0)
      {
        dis+=dd[x0*N];
        break;
      }
      for(j=0;j<=N;j++)
      {
        if(v[j*N+i]>=0.99)
        {
          dis+=dd[x0*N+j];
          x0=j;
        }
        break;
      }
    }
    fprintf(fp,"distance = &f\n",dis);
    //output the result of neuron matrix
    for(i=0;i<N;i++)
      x0=i*N;
    for(j=0;j<N;j++)
      fprintf(fp,"%3.1f",v[x0+j]);
    fprintf(fp,"\n");
  }
  else
  {
    fprintf(fp,"wrong path\n");
    printf("wrong path\n");
  }
}

int recheck()
{
  int i,j,x0;
  float k;
  for(i=0;i<NN;i++)
    if((v[i]>0.01)&&(v[i]<0.99))
      return 0;
  //neuron's state must access 0 or 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;
  }
  //ervey row 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;
  }
  //ervey column have and only have one 1
  return 1;
}