/* program for measure 1 
 mpicc hopfield_measure1.c -lm -o hopfield_measure1c
 mpirun -np 2 hopfield_measure1c N M filename */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include <malloc.h>
#include <time.h>

#define rnd() ( (double)rand() / 0x7fff )

int N,M;
char filename[30];

double **distance;
double **unit,Acoef,Bcoef,Dcoef;
double **city_pos;

int myrank,size;

void initialize(){
  char buf[256],*buf2;
  FILE *fp;
  int i,j;
  double dtotal;

  i=0;
  fp = fopen(filename,"r");
  if ( fp == NULL ){
    exit(1);
  }
  for(i=0;i<N;i++){
    if ( fgets(buf,sizeof(buf),fp) == 0 ){
      break; 
    }
    // sscanf is not work
    buf2 = strtok(buf,",");
    city_pos[i][0] = atof(buf2);
    buf2 = strtok(NULL,",");
    city_pos[i][1] = atof(buf2);
  }

  /* printf("city_pos[1][0]=%lf, [1][1]=%lf\n", city_pos[1][0], city_pos[1][1]); */
  fclose(fp);

  Acoef = 1.0;
  Bcoef = 1.0;
  Dcoef = 2.0;

  for(i=0;i<N;i++){
    for(j=0;j<N;j++){
      distance[i][j] = 
	sqrt ( (city_pos[i][0] - city_pos[j][0]) * (city_pos[i][0] - city_pos[j][0]) +
	       (city_pos[i][1] - city_pos[j][1]) * (city_pos[i][1] - city_pos[j][1]) );
      dtotal += distance[i][j];
    }
  }

  for(i=0;i<N;i++){
    for(j=0;j<N;j++){
      distance[i][j] = 10.0 * distance[i][j] / dtotal;
      unit[i][j] = rnd();
    }
  }

  
}
  
void update_state(){
  int i,j,n,m,jm,jp;
  double un,unitin;
  double aterm,bterm,dterm;

  for (i=0;i<N;i++){
    for(j=0;j<N;j++){
      aterm = bterm = dterm = 0.0;
      for (n=0;n<N;n++) aterm += unit[i][n];
      aterm = -Acoef*(aterm - unit[i][j]);

      for (m=0;m<N;m++) bterm += unit[m][j];
      bterm = -Bcoef*(bterm-unit[i][j]);

      if(j-1 == -1) jm = N-1; else jm=j-1;
      if(j+1 == N) jp=0; else jp = j+1;
      for(m=0;m<N;m++){
	dterm += distance[i][m] * (unit[m][jp] + unit[m][jm]);
      }
      dterm = -Dcoef*dterm;

      unitin = aterm + bterm + dterm + Acoef + Bcoef;
      unit[i][j] = 0.5 * (1.0 + tanh(unitin/0.5) );
    }
  }
}

double energy(){
  int i,j,m,jp,jm;
  double term1,term2,term3;

  term1 = term2 = term3 = 0.0;

  for(i=0;i<N;i++){
    for(j=0;j<N;j++){
      for(m=0;m<N;m++){
	if(m!=j) term1 += unit[i][j]*unit[i][m];
      }
    }
  }
  term1 = 0.5 * Acoef * term1;

  for (j=0;j<N;j++){
    for(i=0;i<N;i++){
      for(m=0;m<N;m++){
	if(m!=i) term2 += unit[i][j]*unit[m][j];
      }
    }
  }
  term2 = 0.5 * Bcoef * term2;

  for(i=0;i<N;i++){
    for(j=0;j<N;j++){
      if(j-1==-1) jm = N-1; else jm=j-1;
      if(j+1==N) jp=0; else jp=j+1;
      for(m=0;m<N;m++) term3 += distance[i][m]*unit[i][j]*(unit[m][jp] + unit[m][jm]);
    }
  }
  term3 = 0.5 * Dcoef * term3;

  return (term1+term2+term3);
}

    

main(int argc,char *argv[]){
  int i;
  double elapsed_time1,elapsed_time2;
  MPI_Status status;
  MPI_Init(&argc,&argv);
  MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
  N = atoi(argv[1]);
  M = atoi(argv[2]);
  strcpy(filename,argv[3]);
  
  distance = (double **)malloc(N * sizeof(double *));
  unit = (double **)malloc(N * sizeof(double *));
  city_pos = (double **)malloc(N * sizeof(double *));
  for(i=0;i<N;i++){
    distance[i] = (double *)malloc(N * sizeof(double));
    unit[i] = (double *)malloc(N * sizeof(double));
    city_pos[i] = (double *)malloc(2 * sizeof(double));
  }

  initialize();

  elapsed_time1 = (double)clock();
  for(i=0;i<M;i++){
    update_state();
    energy();
  }
  elapsed_time2 = (double)clock();
  if (myrank==0) printf("elapsed time:%lf second\n",(elapsed_time2-elapsed_time1)/CLOCKS_PER_SEC);

  for(i=0;i<N;i++){
    free(distance[i]);
    free(unit[i]);
    free(city_pos[i]);
  }
  free(distance);
  free(unit);
  free(city_pos);
  

  MPI_Finalize();
}