#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "grid.h"

#define Steps 400

int main(){
  int i, j;
  bool control_mean;
  FILE *W_score, *W_average, *W_aux;

  // Random Sequence Initialization
  srand(1); 

  // Making dirs
  make_dirs();
  sprintf(filename,"%s%cEvolGridScores.dat",foldername, 47);
  W_score = fopen(filename,"w");
  sprintf(filename,"%s%cEvolMeanGridScores.dat",foldername, 47);
  W_average = fopen(filename,"w");
  
  // Setting place cells in a lattice
  set_input();
  sprintf(filename,"%s%cGrid_PlaceUnits.dat",foldername, 47);
  W_aux=fopen(filename,"w");
  for(i=0;i<numhippunit;i++){    
    fprintf(W_aux,"%d %f %f\n", i+1, x_place_cell[i], y_place_cell[i]);
    } // end for i
  fclose(W_aux);

  // Initializating dynamics (input, adaptation variables, and output + weights from Hipp to EC + initial location)
  initialize();
  sprintf(filename,"%s%cProps.dat",foldername, 47);
  W_aux=fopen(filename,"w");
  for(i=0;i<numecunit;i++){
    fprintf(W_aux,"%d %f %f %f %f %f %f\n",i+1,(180./M_PI)*ConjCell_pref_HD[i],ConjCell_X[i],ConjCell_Y[i],b1[i],nu[i],rho[i]);
	}
  fclose(W_aux);
  sprintf(filename,"%s%cLat_weights.dat",foldername, 47);
  W_aux=fopen(filename,"w+");
  for(i=0;i<numecunit;i++){
    for(j=0;j<numecunit;j++) fprintf(W_aux,"%.3f ",weight_lateral[i][j]);
    fprintf(W_aux,"\n");
    }
  fclose(W_aux);

  // Evolution
  for(idstep=0;idstep<plotstep*Steps;idstep++){
    if(idstep%10000==0) printf("Iteration %d\n",idstep);

    // Update Rat Position
    update_position();

    // Updating neurons according to current position
    control_mean = (idstep > 15);
    update_neurons(control_mean);

    // Storing history for some variables
    idhis = idstep%lenhis;
    speed_list[idhis] = speed; // nevertheless, speed is constant in this simulation
    x_list[idhis] = x;
    y_list[idhis] = y;
    theta_list[idhis] = theta; // movement direction
    mean_list[idhis] = mean_network_act;
    sparsity_list[idhis]= sparsity;
    
    // Once output of the system is defined, learning in connections is applied
    update_weights();
    update_mean();

    // Update positional maps of activity of individual neurons
    // (according to a smoothing procedure, which takes into account ~100 historical values)
    update_maps();
      
    // Creates auto-correlograms, separated by plotstep steps 
    if(idstep%plotstep==0){	  
	  printf("Step %d out of %d: \n" , idstep, plotstep*Steps);
      compute_autocorrelation(auto_corre, map_pos);
      compute_gridness(auto_corre,grid_score,0);

      fprintf(W_score,"%d ",idstep);
      for(i=0;i<numecunit;i++)
	    fprintf(W_score,"%.3f ",grid_score[i]);
      fprintf(W_score,"\n");
      fflush(W_score);
	  
      compute_statistics(grid_score);
	  fprintf(W_average,"%d %.3f %.3f\n",idstep,mean,sqrt(var));      
      fflush(W_average);

      save_all();

	  } // end if idstep     

    } // end for idstep

  fclose(W_score);
  fclose(W_average);
  
  return(1);
}

// makes directories, basically one for each neuron in mec
void make_dirs(){
  int i;
  char command[500];

  sprintf(foldername,"Simulation-%d",1);
  sprintf(command,"mkdir %s", foldername);
  system(command);

  for(i=0;i<numecunit;i++){
    sprintf(command,"mkdir %s%ccell_%d", foldername, 47, i+1);
	system(command);
    }

  return;
}

//generates center coordinates of input maps that are like hippocampal place cells arranged evenly across the environment
void set_input(){
  int i, j, total, n;
  float r, spacing, radius, ang, coef_a0, coef_a1, coef_b0, coef_b1, coef_c0, coef_c1;

  total = 0;
 
  // Primer elemento a radio nulo
  i = 0; total++;
  radius = 0.;
  ang = 0.;
  x = radius*cos(ang);
  y = radius*sin(ang);
  r = map_perturb*(-.5+ uniform());
  x_place_cell[i] = x + environmentsize/2. + r;
  r = map_perturb*(-.5+ uniform());
  y_place_cell[i] = y + environmentsize/2. + r;
  
  spacing = 2.*M_PI*(environmentsize/2.) / (float)ext_sidehippgrid;
  
  for(i=1;i<=Nrings;i++){
    radius = (i/(float)Nrings)*(environmentsize/2.);
    n = (int) floor(2.*M_PI*radius/spacing + 0.5);
	ang = (2.*M_PI/(float)n) * uniform();
	for(j=0;j<n;j++){
	  ang += (2.*M_PI/(float)n);
      x = radius*cos(ang);
      y = radius*sin(ang);
      r = map_perturb*(-.5+ uniform());
      x_place_cell[total] = x + environmentsize/2. + r;
      r = map_perturb*(-.5+ uniform());
      y_place_cell[total] = y + environmentsize/2. + r;
	  total++;
	  }
	}

  // Una capa mas por afuera del circulo
  radius = ((Nrings+1)/(float)Nrings)*(environmentsize/2.);
  n = (int) floor(2.*M_PI*radius/spacing + 0.5);
  ang = (2.*M_PI/(float)n) * uniform();
  for(j=0;j<n;j++){
    ang += (2.*M_PI/(float)n);
    x = radius*cos(ang);
    y = radius*sin(ang);
    r = map_perturb*(-.5+ uniform());
    x_place_cell[total] = x + environmentsize/2. + r;
    r = map_perturb*(-.5+ uniform());
    y_place_cell[total] = y + environmentsize/2. + r;
    total++;
    }

  numhippunit = total;
  nconnhipp2ec = total;

  printf("Hippocampal units: %d\n",numhippunit);

  // Setting time constants in the network
  coef_a0 = (exp(b1_min/B1) - (numecunit/2)*exp(b1_max/B1)) / (exp(b1_min/B1) - exp(b1_max/B1));
  coef_a1 = (float)((numecunit/2)-1) / (exp(-b1_max/B1) - exp(-b1_min/B1));
  coef_b0 = 1.0;
  coef_b1 = -0.8*exp(b1_min/B1);
  coef_c0 = 0.2;
  coef_c1 = 0.2*exp(b1_min/B1);
  for(i=0;i<numecunit/2;i++){
    b1[i] = B1*log(coef_a1/((float)(i+1)-coef_a0));
    b1[numecunit-1-i] = b1[i];

    b2[i] = b1[i]/3.;
    b2[numecunit-1-i] = b2[i];

    nu[i] = coef_b0 + coef_b1*exp(-b1[i]/B1);
    nu[numecunit-1-i] = nu[i];

    rho[i] = coef_c0 + coef_c1*exp(-b1[i]/B1);
    rho[numecunit-1-i] = rho[i];
	}

  return;
}

//initialize variables
void initialize(){
  int i,j,k,l;
  float sum_weight, angle, dist;

  speed = constspeed;
  threshold = 0.;

  for(i=0;i<numecunit;i++){
    act_ec[i] = 0.;       // entorhinal neurons state
    mean_act_ec[i] = 0.;
    active_ec[i] = 0.;    // adaptation variable
    inactive_ec[i] = 0.;  // adaptation variable
    input2ec[i] = 0.;     // activation of previous variables 
    for(k=0;k<delay;k++){
      delayed_act_ec[i][k] = 0.;
	  } // end for k
    for(k=0;k<numbinmap;k++){
      map_ang[i][k] = 0.;
      for(l=0;l<numbinmap;l++){
        map_pos[i][k][l] = 0.;
		} // end for l
	  } // end for k
	} // end for i

  for(j=0;j<numhippunit;j++){
    act_hipp[j] = 0.;     //hippocampal neurons state
    mean_act_hipp[j] = 0.;
    } // end for j
  
  for(i=0;i<numecunit;i++){
    sum_weight = 0.;
    for(k=0;k<nconnhipp2ec;k++){
	  weight_hipp2ec[i][k] = 0.9 + 0.1*uniform(); //random weights (from hc to ec) initialization
	  sum_weight +=  weight_hipp2ec[i][k]*weight_hipp2ec[i][k];
      } // end for k
      
    for(k=0;k<nconnhipp2ec;k++)//normalization
	  weight_hipp2ec[i][k] *= sqrt(1./sum_weight);
    
    } // end for i

  // Lateral connectivity - HD cell + Virtual associated cell
  for(i=0;i<numecunit;i++){
    ConjCell_pref_HD[i] = 2.*M_PI*uniform();
    do{
      ConjCell_X[i] = environmentsize * uniform();
      ConjCell_Y[i] = environmentsize * uniform();
      }
	while(isin(ConjCell_X[i],ConjCell_Y[i])==0);
	}

  // Set weights
  for(i=0;i<numecunit;i++){
    sum_weight = 0.;
    for(k=0;k<numecunit;k++){
      if(k!=i && (fabs(k-i)<=numside || fabs(k+numecunit-i)<=numside || fabs(k-numecunit-i)<=numside)){

	    angle = atan2(ConjCell_Y[i]-ConjCell_Y[k], ConjCell_X[i]-ConjCell_X[k]);
		angle = modulus(angle,2.*M_PI);

        dist = sqrt( pow(ConjCell_X[i]-ConjCell_X[k]-expect_length*cos(angle),2) + pow(ConjCell_Y[i]-ConjCell_Y[k]-expect_length*sin(angle),2) );
  	    weight_lateral[i][k] = fdir(angle,ConjCell_pref_HD[k],k)*fdir(angle,ConjCell_pref_HD[i],i)*exp(-dist*dist/(2.*spreadlateral*spreadlateral)) - kappa;
        if (weight_lateral[i][k] < 0.) weight_lateral[i][k] = 0.;
		sum_weight +=  weight_lateral[i][k]*weight_lateral[i][k];
  	    }
  	  else
  	    weight_lateral[i][k] = 0.;
	  } // end for k

    if(sum_weight>0.){
	  for(k=0;k<numecunit;k++)
  	    weight_lateral[i][k] *= sqrt(1./sum_weight);
      }

	} // end for i
  
  // Initial position
  x = -1.;
  y = -1.;
  while(isin(x,y)==0){ // While the position is outside, keep doing. Until position is in the box
    x = environmentsize * uniform();
    y = environmentsize * uniform();
    }

  return;
}

//update the position with a step v. If the new position is out of the box, choose another one
void update_position(){
  float theta0, x0, y0, spread_turn0;
  float par, coef_a, coef_b, coef_c, xt, yt, alfa, theta0r;
  int i, count = 0;
	
  spread_turn0 = spreadturn;

  do{
    theta0 = theta + spread_turn0 * normal();
    x0 = x + speed * cos(theta0);
    y0 = y + speed * sin(theta0);
    
    if(isin(x0,y0)==0){
      coef_a = 1.;
	  coef_b = 2.0*cos(theta0)*(x-environmentsize/2.0) + 2.0*sin(theta0)*(y-environmentsize/2.0);
      coef_c = (x-environmentsize/2.0)*(x-environmentsize/2.0) + (y-environmentsize/2.0)*(y-environmentsize/2.0) - environmentsize*environmentsize/4.0;
	  par = (-coef_b + sqrt(coef_b*coef_b-4.0*coef_a*coef_c))/(2.0*coef_a);
	  	  
      xt = x + par * cos(theta0);
      yt = y + par * sin(theta0);

      alfa = atan2(-(xt-environmentsize/2.0),(yt-environmentsize/2.0));
      
	  // reflejamos
	  theta0r = 2.0*alfa - theta0;
	  theta0 = theta0r;
      x0 = xt + (speed-par) * cos(theta0);
      y0 = yt + (speed-par) * sin(theta0);
	  }

    //spread_turn0 *=1.1;
    //if(spread_turn0> 2.*M_PI) spread_turn0 = 2.*M_PI;
    
    //count++; 
    //if(count>50){
    //  printf("Stuck near a wall !!!!!!\n");
    //  count = 0;
    //  }

    } //end do
  while(isin(x0,y0)==0); // up to position is inside (at least, one updating)

  // Updating global variables
  theta = modulus(theta0,2.*M_PI);
  x = x0;
  y = y0;
  
  x_map = (int) floor(x*numbinmap/environmentsize);	// Pixeling the position (corresponding bin)
  y_map = (int) floor(y*numbinmap/environmentsize);
  theta_map = (int) floor(theta*numbinmap/(2.*M_PI));

  return;
}

//main part of the program, updates all mEC neurons by integrating feedforward activity
void update_neurons(bool control_mean){
  int i,j,k;
  float input, factor;

  for(j=0;j<numhippunit;j++){
    act_hipp[j]= fhc(x - x_place_cell[j],y - y_place_cell[j]);     //activation in hippocampus (time t)
    }

  for(i=0;i<numecunit;i++){
    // Updating input to EC (from Hippocampus - time t)
    input = 0.;
    for(k=0;k<nconnhipp2ec;k++)
      input += weight_hipp2ec[i][k]*act_hipp[k];

    input2ec[i] = input;

    // Updating input to EC (from lateral connectivity - time t)
	input = 0.;
	for(k=0;k<numecunit;k++)
	  input += weight_lateral[i][k]*delayed_act_ec[k][0];

    //factor = (float)idstep/(float)(plotstep*Steps/4.);
    //if (factor > 1.) factor = 1.0;
    factor = 1.0;
    input2ec[i] += rho[i]*factor*input;
    input2ec[i] *= fdir(theta, ConjCell_pref_HD[i],i);

    // Updating adaptation variables, including selection of most/less active cells (time t+1)
    active_ec[i] += b1[i]*(input2ec[i]-inactive_ec[i]-active_ec[i]);  //adaptation variable (activation)
    inactive_ec[i] += b2[i]*(input2ec[i]-inactive_ec[i]);             //adaptation variable (inactivation)
    
    if(i==0){
      min_active_ec = active_ec[0];  // EC neuron less active (initial arbitrary selection)
      max_active_ec = active_ec[0];  // EC neuron more active (initial arbitrary selection)
       }
    else{
      if(active_ec[i] < min_active_ec) min_active_ec = active_ec[i];
      if(active_ec[i] > max_active_ec) max_active_ec = active_ec[i];
      }

	} // end for i

  networkoutput(control_mean);

  return;
}

void networkoutput(bool control_mean){
  int i, j;
  
  if(!control_mean){
    ec_output();
    return; 
    } // get out without any control
  
  n_iter = 0;

  bool fixed = false;
  bool go_on = true;
  
  if(threshold > max_active_ec - 0.01) //choose a good initial threshold
    threshold = max_active_ec - sparsitylevel * (max_active_ec - min_active_ec);
  
  do{
    fix_sparsity();                   //fix sparsity (but maybe not, depending on "nummaxiter").
    fixed = is_fixed();

    if(!fixed &&  n_iter <= nummaxiter){
      comput_gain_gradient();        //fix mean by gradient descent
      adapt_gain();
      }

    go_on = ((!fixed) && (n_iter <= nummaxiter)); // to close the loop, but maybe not fixed!
    } // end do
  while(go_on);

  ec_output(); // includes the last update

  // Updates delayed activity
  for(i=0;i<numecunit;i++){
     for(j=0;j<delay-1;j++){
       delayed_act_ec[i][j] = delayed_act_ec[i][j+1];
	   }
     delayed_act_ec[i][delay-1] = act_ec[i];
	 }

  return;
}

void fix_sparsity(){
  float cur_spar_error;

  bool fixed = find_threshold_range(max_active_ec);  // sums one iteration
  bool go_on = (!fixed);

  while(go_on){
    ec_output(); // With the new threshold selected in find_threshold_range: (lower+upper)/2
    n_iter++;

    cur_spar_error = fabs(sparsity - sparsitylevel);
    go_on =  ((cur_spar_error > tolerancemean * sparsitylevel)  && n_iter <= nummaxiter);

    if(go_on){  //new threshold
	  
      if(sparsity < sparsitylevel){  //threshold is too high, move to a lower value
        upper_threshold = threshold;
        threshold = (lower_threshold + threshold) / 2.0;
        }
      else{
        lower_threshold = threshold;
        threshold = (upper_threshold + threshold) / 2.0;
        } // end else

      } // end if go_on

    // A go condition can be reached by iterations, with a poor performance in cur_spar_error.
    } // end while 

  inc_threshold = upper_threshold - lower_threshold;
  inc_threshold = Max(inc_threshold, 0.005);
}

//find a threshold range to bound the sparsity
//return true if sparsity criteria is satisfied
bool find_threshold_range(float max_act){
  int direction;
  float cur_spar_error;
  
  ec_output();
  n_iter++;

  cur_spar_error = fabs(sparsity - sparsitylevel);
  if(cur_spar_error < tolerancemean * sparsitylevel)
    return true;

  if(sparsity < sparsitylevel){
    direction = -1;
    search_threshold(direction, upper_threshold, lower_threshold); //search for a lower threshold
    }
  else{
    direction = 1;
    search_threshold(direction, lower_threshold, upper_threshold); //search for a higher threshold
    }

  upper_threshold = Min(upper_threshold, max_act);

  threshold = (lower_threshold + upper_threshold) / 2.0;//save threshold for bisect search.
     
  return false; //need to do bisection further.
}

void ec_output(){
  int i;
  float act;

  mean_network_act = 0;
  mean_square_network_act = 0;

  for(i=0;i<numecunit;i++){
    act = active_ec[i] - threshold;
    
    if(act < 0.)
      act_ec[i]=0.;
    
    else
      act_ec[i]= firingconst * atan(gain * act);//state of ec unit. maximal rate is 1.

    mean_network_act += act_ec[i];
    mean_square_network_act += act_ec[i]*act_ec[i];
    } // end for

  mean_network_act = mean_network_act / (float)numecunit;  //average activity in EC: <r>
  if(mean_square_network_act > 1e-6)
    mean_square_network_act = mean_square_network_act / (float)numecunit;  //average dispersion in EC <r^2>
  else 
    mean_square_network_act = 1e-6;

  sparsity = mean_network_act * mean_network_act / mean_square_network_act;

  return;
}

//search for threshold range, return to root.
//When direction = 1, ascending search, and find [back_bound, ahead_bound].
//When direction = -1, descending search, and find [ahead_bound, back_bound]
void search_threshold(int direction, float &back_bound, float &ahead_bound){
  bool go_on = true;
  bool search_low = (direction == -1);

  inc_threshold = Max(inc_threshold, 0.005);  // inc_threshold is posteriorly updated according to upper/lower values; actually, this line doesn't care (it is already done before)
 
  back_bound = threshold; // current threshold is good as a back bound.
  
  while(go_on){
    threshold += direction * inc_threshold;

    ec_output();          // re-calculate network global state with this new threshold
    n_iter++;

    if((sparsity > sparsitylevel && search_low) || (sparsity < sparsitylevel && !search_low)){  //ahead bound found.
      ahead_bound = threshold; 
      go_on = false;
      }
    else{ //update back bound.
      back_bound = threshold;  
      inc_threshold *= 4;
      } 
    } // end while

  return;
}

bool is_fixed(){
  bool fixed = false;

  float cur_spar_error = fabs(sparsity - sparsitylevel);
  float cur_mean_error = fabs(mean_network_act - meanactivitylevel);

  fixed  = ((cur_spar_error < tolerancemean * sparsitylevel) && (cur_mean_error < tolerancemean * meanactivitylevel));
 
  return fixed;
}

void comput_gain_gradient(){
  int i;
  float dif, d_gain, denorm;
  
  d_mean_gain = 0.;

  for(i=0;i<numecunit;i++){
    dif = active_ec[i] - threshold;

    d_gain = 0.;
    if(dif > 0.){
      denorm = 1. / (1. + gain * gain *dif * dif);
	  d_gain = firingconst * dif * denorm;
	  }
      
    d_mean_gain += d_gain;
    } // end for
  
  d_mean_gain/= (float)numecunit;

  return;
}

void adapt_gain(){
  float f_mean, inc_gain, len, new_gain;

  f_mean = mean_network_act - meanactivitylevel;
  inc_gain = 0;
  
  if(fabs(d_mean_gain) > 1e-7){      // if mean gradient obtained before is significant
    inc_gain = f_mean / d_mean_gain; // increment in g
    len = fabs(inc_gain);
    if(len > 1000)
      inc_gain *= 1000.0 / len;       // limited to 1000 (+ or -)
    }
  else{  // mean gradient too small; not useful to update gain
    if(fabs(f_mean) > tolerancemean * meanactivitylevel){
      if(f_mean > 0) //mean too high, reduce gain
	    inc_gain = gain/(1.0+n_iter);
      else
	    inc_gain = -gain/(1.0+n_iter);
      }
    }

  new_gain = gain - inc_gain;

  if(new_gain < mingain || new_gain > maxgain){ //keep in safe range.
    new_gain = Max(new_gain, mingain);
    new_gain = Min(new_gain, maxgain);
    }

  gain = new_gain;

  return;
}

//learning weights
void update_weights(){
  int i,k;
  float sum_weight, lrweight;

  lrweight = lrweight1 + (lrweight2-lrweight1)*idstep/((float)plotstep*Steps/4.);
  if (lrweight < lrweight2) lrweight = lrweight2;

  for(i=0;i<numecunit;i++){
    sum_weight = 0.;
	
    for(k=0;k<nconnhipp2ec;k++){
      weight_hipp2ec[i][k] += lrweight * (act_ec[i]*act_hipp[k] - mean_act_ec[i] * mean_act_hipp[k]);
      
      if(weight_hipp2ec[i][k]<0.)
	    weight_hipp2ec[i][k]=0.;
	  else
	    sum_weight += weight_hipp2ec[i][k] * weight_hipp2ec[i][k]; 
	  } // end for k
      
    if(sum_weight > 1e-20){
	  for(k=0;k<nconnhipp2ec;k++)//normalization
	    weight_hipp2ec[i][k] *= sqrt(1./sum_weight);
      }
    else{
	  printf(" Weight diverge %e!\n", sum_weight);
	  for(k=0;k<nconnhipp2ec;k++)
	    weight_hipp2ec[i][k] = 1.0/sqrt((float)nconnhipp2ec);	  
      } // end else
    
    } // end for i

  return;
}

// mean activity in each neuron (averaging ~20 time steps)
void update_mean(){
  int i;

  for(i=0;i<numhippunit;i++)
    mean_act_hipp[i] += lrmean * (act_hipp[i] - mean_act_hipp[i]);
	  
  for(i=0;i<numecunit;i++)
    mean_act_ec[i] += lrmean * (act_ec[i] - mean_act_ec[i]);

  return;
}

//updates the map of each mec cell by averaging the last few times the rat passed through this place
void update_maps(){
  int i, j, pos_x, pos_y;
  
  for(i=0;i<numecunit;i++){
    map_pos[i][x_map][y_map] += lrpositionmap*(act_ec[i]-map_pos[i][x_map][y_map]);
    map_ang[i][theta_map] += lranglemap*(act_ec[i]-map_ang[i][theta_map]);
    }

  return;
}

//obtain auto-correlograms 
void compute_autocorrelation(float ac[numecunit][2*numbinmap-1][2*numbinmap-1], float map[numecunit][numbinmap][numbinmap]){
  int sideY = numbinmap, sideX = numbinmap;
  int i, xOverlap, yOverlap, nOverlap, iXMin, iXMax, iX, iY;
  float aux, aux1, aux2, aux3, aux4, aux5, act1, act2;

  for(i=0;i<numecunit;i++)
    for(yOverlap = 1; yOverlap< sideY+1; yOverlap++)
      for(xOverlap =1; xOverlap < 2 * sideX; xOverlap++){
		nOverlap = 0;

		aux=0.;
		if(xOverlap > sideX){
          nOverlap = yOverlap *(2* sideX- xOverlap);
          iXMin = xOverlap - sideX;  // when xOverlap = sideX+1 :: all pixels except one are available, corresponds to tauX = (+/-) 1
          iXMax = sideX;             // when xOverlap = 2*sideX-1 :: just one pixel is available, corresponds to tauX = (+/-)(numbinmap-1) 
		  }
		else{
          nOverlap = xOverlap * yOverlap;
          iXMin = 0;         
          iXMax = xOverlap;  // when xOverlap = sideX :: all pixels are available, corresponds to tauX = 0
          }                  // when xOverlap = 1 :: just one pixel is available, corresponds to tauX = (+/-)(numbinmap-1)

		if(nOverlap > 6){
          aux1=0.; aux2=0.; aux3=0.; aux4=0.; aux5=0.;
          for(iY =0; iY < yOverlap; iY++)        // which pixels to span
            for(iX = iXMin; iX < iXMax; iX++){   // which pixels to span
			  act1 = map[i][iY][iX];
			  act2 = map[i][sideY-yOverlap+iY][sideX-xOverlap+iX];

			  aux1 += act1 * act2;
			  aux2 += act1;
			  aux3 += act2;
			  aux4 += act1 * act1;
			  aux5 += act2 * act2;
              } // end for iX

          // all sums have been calculated, for a given (tauX,tauY), whenever statistics is adequate (n>6)
          aux = (nOverlap*aux4 - aux2*aux2)*(nOverlap*aux5 - aux3*aux3);
          if(aux<1e-20){ // if variance is negligible across any two locations separated by (tauX,tauY) 
            aux = 0;
            }
          else
            aux = (nOverlap*aux1 - aux2*aux3)/sqrt(aux);

          } // end if(nOverlap > 6)

		ac[i][yOverlap-1][xOverlap-1] = aux;
		ac[i][2*sideY-yOverlap-1][2*sideX-xOverlap-1] = aux;
		} // end for xOverlap

  return;
}

void compute_gridness(float ac[numecunit][2*numbinmap-1][2*numbinmap-1], float score[numecunit], int flag_print){
  int i;
  float UpperLimit, rmax_aux, aux;
  FILE *W1;

  assess_rmin(ac);
  //for(i=0;i<numecunit;i++) rmin[i] = 10.;

//  printf("  rmin\n");

  UpperLimit = ((numbinmap)*environmentsize/(float)numbinmap);

  for(i=0;i<numecunit;i++){
//    printf("    %d: %f -",i+1,rmin[i]);

    if(flag_print){
      sprintf(filename,"%s%ccell_%d%cScoreVsRadius-%d.dat",foldername, 47, i+1, 47, idstep);
      W1=fopen(filename,"w+");
      }

    rmax_aux = rmin[i] + environmentsize/10.; // Starts with an offset
    if(rmax_aux<Rmax_min) rmax_aux = Rmax_min;
	score[i] = -10.;
    rmax[i] = rmax_aux;
    while(rmax_aux<UpperLimit-environmentsize/10.){
      aux = assess_corr(i,ac,rmin[i],rmax_aux,0);
	  if(flag_print) fprintf(W1,"%f %f\n",rmax_aux,aux);
	  if(aux>score[i]){
	    score[i] = aux;
	    rmax[i] = rmax_aux;
	    }
      rmax_aux += environmentsize/(float)numbinmap;
  	  } // end while

//    printf(" %f\n",rmax[i]);

    if(flag_print) fclose(W1);

    // Once rmax[i] is determined, we print correspondind maps (if asked)
    aux = assess_corr(i,ac,rmin[i],rmax[i],flag_print);
    
    } // end for i

  return;
}

void assess_rmin(float ac[numecunit][2*numbinmap-1][2*numbinmap-1]){
  int i, j, k, n;
  float Dx, Dy;
  float UpperLimit, BinWidth;
  int Position, flag1, flag2, flag3;
  
  UpperLimit = sqrt(2)*((numbinmap)*environmentsize/(float)numbinmap);
  BinWidth = UpperLimit/(float)Nbins;
  
  for(i=0;i<numecunit;i++){ // loop over all EC cells

    for(j=0;j<(2*numbinmap-1);j++){   // loop over X in autocorrelogram
      for(k=0;k<(2*numbinmap-1);k++){ // loop over Y in autocorrelogram
	    Dx = ((j+1) - numbinmap)*environmentsize/(float)numbinmap;
	    Dy = ((k+1) - numbinmap)*environmentsize/(float)numbinmap;
	    Dist[j][k] = sqrt(pow(Dx,2)+pow(Dy,2));
        Corr[j][k] = ac[i][j][k];
		} // end for k 
      } // end for j

    for(j=0;j<Nbins;j++){
      BinData[i][j] = 0.;
      NData[j] = 0;
      }

    // Ordering data
    for(j=0;j<(2*numbinmap-1);j++){
      for(k=0;k<(2*numbinmap-1);k++){
        Position = (int) floor(Dist[j][k]/BinWidth);
        if(Position<Nbins)
		  BinData[i][Position] += Corr[j][k];
        else printf(" Error in determining Rmin\n");
		NData[Position]++;
		} // end for k 
      } // end for j

    rmin_aux[i] = 0;
    corr_aux[i] = 1.1;
    flag1 = 1;
    flag2 = 1;

    for(j=0;j<Nbins;j++){
      if(NData[j]!=0){
        BinData[i][j] /= (float)(NData[j]);
        if(flag1 && BinData[i][j]<corr_aux[i]){ // always decreasing up to minimum
          rmin_aux[i] = (j+0.5)*BinWidth;
          corr_aux[i] = BinData[i][j];
          
          rmin_extreme[i] = rmin_aux[i];        // stores the minimum value
		  if(flag2) rmin[i] = rmin_aux[i];      // by default; it is not, if 0.2 is reached before
		  if(flag2 && corr_aux[i]<0.2){
		    // interpolate
		    flag3 = 1;
		    n = 1;
			while(flag3 && (j-n) >= 0){
			  if(NData[j-n]==0){
			    n++;
				}
			  else flag3 = 0;
			  }
			// Interpolation
			rmin[i] = rmin_aux[i] + ((0.2-corr_aux[i])/(BinData[i][j-n]-corr_aux[i]))*((j-n+0.5)*BinWidth-rmin_aux[i]);
			flag2 = 0;
			} // end if (flag2)
		  } // end if (flag1)
		else
		  flag1 = 0; // leaves the previous loop (the first minimum was already caught) 
		}
  	  }

    } // end i
  
  return;
}

float assess_corr(int i, float ac[numecunit][2*numbinmap-1][2*numbinmap-1], float r1, float r2, int print){
  int j, k, n;
  bool flag1;
  float Dx, Dy, Radius;
  float m0, m30, m60, m90, m120, m150;
  float v0, v30, v60, v90, v120, v150;
  float cross30, cross60, cross90, cross120, cross150;
  float corr30, corr60, corr90, corr120, corr150, corr_even, corr_odd;
  float aux1, aux2, corr;
  FILE *W0,*W30,*W60,*W90,*W120,*W150;
  
  if(print){
    sprintf(filename,"%s%ccell_%d%cOriginalMap-%d.dat",foldername, 47, i+1, 47, idstep);
    W0=fopen(filename,"w+");
    sprintf(filename,"%s%ccell_%d%cRotatedMap_30degrees-%d.dat",foldername, 47, i+1, 47, idstep);
    W30=fopen(filename,"w+");
    sprintf(filename,"%s%ccell_%d%cRotatedMap_60degrees-%d.dat",foldername, 47, i+1, 47, idstep);
    W60=fopen(filename,"w+");
    sprintf(filename,"%s%ccell_%d%cRotatedMap_90degrees-%d.dat",foldername, 47, i+1, 47, idstep);
    W90=fopen(filename,"w+");
    sprintf(filename,"%s%ccell_%d%cRotatedMap_120degrees-%d.dat",foldername, 47, i+1, 47, idstep);
    W120=fopen(filename,"w+");
    sprintf(filename,"%s%ccell_%d%cRotatedMap_150degrees-%d.dat",foldername, 47, i+1, 47, idstep);
    W150=fopen(filename,"w+");
	}


  n = 0;
  m0 = 0.; v0 =0.;
  m30 = 0.; v30 = 0.; cross30 = 0.; corr30 = 0.;
  m60 = 0.; v60 = 0.; cross60 = 0.; corr60 = 0.;
  m90 = 0.; v90 = 0.; cross90 = 0.; corr90 = 0.;
  m120 = 0.; v120 = 0.; cross120 = 0.; corr120 = 0.;
  m150 = 0.; v150 = 0.; cross150 = 0.; corr150 = 0.;
  
  rotate(i,ac,ac_rot30,30.*M_PI/180.0,r1,r2);
  rotate(i,ac,ac_rot60,60.*M_PI/180.0,r1,r2);
  rotate(i,ac,ac_rot90,90.*M_PI/180.0,r1,r2);
  rotate(i,ac,ac_rot120,120.*M_PI/180.0,r1,r2);
  rotate(i,ac,ac_rot150,150.*M_PI/180.0,r1,r2);

  for(j=0;j<2*numbinmap-1;j++){
	for(k=0;k<2*numbinmap-1;k++){
      // Truncated map within the anulus
      Dx = ((j+1) - numbinmap)*environmentsize/(float)numbinmap;
      Dy = ((k+1) - numbinmap)*environmentsize/(float)numbinmap;
      Radius = sqrt(pow(Dx,2)+pow(Dy,2));
      // Anulus condition
	  flag1 = (Radius>=r1 && Radius<r2);
	  
      if(flag1){
		m0 += ac[i][j][k];
		m30 += ac_rot30[j][k];
		m60 += ac_rot60[j][k];
		m90 += ac_rot90[j][k];
		m120 += ac_rot120[j][k];
		m150 += ac_rot150[j][k];
		n++;

        if(print){
          fprintf(W0,"%.3f ", ac[i][j][k]);
          fprintf(W30,"%.3f ", ac_rot30[j][k]);
          fprintf(W60,"%.3f ", ac_rot60[j][k]);
          fprintf(W90,"%.3f ", ac_rot90[j][k]);
          fprintf(W120,"%.3f ", ac_rot120[j][k]);
          fprintf(W150,"%.3f ", ac_rot150[j][k]);          
		  } // end if print
		} // end if
	  
	  else{
        if(print){
  	      fprintf(W0,"%.3f ", -10.);
	      fprintf(W30,"%.3f ", -10.);
	      fprintf(W60,"%.3f ", -10.);
	      fprintf(W90,"%.3f ", -10.);
	      fprintf(W120,"%.3f ", -10.);
	      fprintf(W150,"%.3f ", -10.);
		  }	// end if print
		} // end else
	  }
    
	if(print){
      fprintf(W0,"\n ");
      fprintf(W30,"\n ");
      fprintf(W60,"\n ");
      fprintf(W90,"\n ");
      fprintf(W120,"\n ");
      fprintf(W150,"\n ");
	  }
	
	}

  m0 /= (float)n; 
  m30 /= (float)n; 
  m60 /= (float)n; 
  m90 /= (float)n; 
  m120 /= (float)n; 
  m150 /= (float)n; 

  for(j=0;j<2*numbinmap-1;j++){
	for(k=0;k<2*numbinmap-1;k++){
      // Truncated map within the anulus
      Dx = ((j+1) - numbinmap)*environmentsize/(float)numbinmap;
      Dy = ((k+1) - numbinmap)*environmentsize/(float)numbinmap;
      Radius = sqrt(pow(Dx,2)+pow(Dy,2));
      // Anulus condition
	  flag1 = (Radius>=r1 && Radius<r2);
	  
      if(flag1){
		v0 += pow(ac[i][j][k]-m0,2);
		v30 += pow(ac_rot30[j][k]-m30,2);
		v60 += pow(ac_rot60[j][k]-m60,2);
		v90 += pow(ac_rot90[j][k]-m90,2);
		v120 += pow(ac_rot120[j][k]-m120,2);
		v150 += pow(ac_rot150[j][k]-m150,2);

        cross30 += (ac[i][j][k]-m0)*(ac_rot30[j][k]-m30);
        cross60 += (ac[i][j][k]-m0)*(ac_rot60[j][k]-m60);
        cross90 += (ac[i][j][k]-m0)*(ac_rot90[j][k]-m90);
        cross120 += (ac[i][j][k]-m0)*(ac_rot120[j][k]-m120);
        cross150 += (ac[i][j][k]-m0)*(ac_rot150[j][k]-m150);
		} // end if
	  }
    }

  corr30 = cross30/sqrt(v0*v30);
  corr60 = cross60/sqrt(v0*v60);
  corr90 = cross90/sqrt(v0*v90);
  corr120 = cross120/sqrt(v0*v120);
  corr150 = cross150/sqrt(v0*v150);

  aux1 = 0.; aux2 = 0.;
  if(corr60<=corr120) aux1 = corr60;
  else aux1 = corr120;
  if(corr30>=corr90 && corr30>=corr150) aux2 = corr30;
  if(corr90>=corr30 && corr90>=corr150) aux2 = corr90;
  if(corr150>=corr30 && corr150>=corr90) aux2 = corr150;
//  corr = aux1 - aux2;
  corr_even = 0.5*(corr60+corr120);
  corr_odd = (1./3.)*(corr30+corr90+corr150);
  corr = corr_even - corr_odd;

  if(print){
    fclose(W0);
    fclose(W30);
    fclose(W60);
    fclose(W90);
    fclose(W120);
    fclose(W150);
    }
  
  return(corr);
}

void rotate(int i, float ac[numecunit][2*numbinmap-1][2*numbinmap-1], float ac_rot[2*numbinmap-1][2*numbinmap-1], float alpha, float r1, float r2){
  bool flag1;
  int j, k, j1, j2, k1, k2;
  float Dx, Dy, Radius, Dx_, Dy_, Dx_rot, Dy_rot, Q11, Q12, Q21, Q22;

  for(j=0;j<2*numbinmap-1;j++){
	for(k=0;k<2*numbinmap-1;k++){
      // Truncated map within the anulus
      Dx = ((j+1) - numbinmap)*environmentsize/(float)numbinmap;
      Dy = ((k+1) - numbinmap)*environmentsize/(float)numbinmap;
      Radius = sqrt(pow(Dx,2)+pow(Dy,2));

      // Anulus condition
	  flag1 = (Radius>=r1 && Radius<r2);

      if(flag1){
  	    // Rotated coordinates (with _ )
	    Dx_ = ((j+1) - numbinmap)*cos(alpha) + ((k+1) - numbinmap)*sin(alpha);
	    Dy_ = -((j+1) - numbinmap)*sin(alpha) + ((k+1) - numbinmap)*cos(alpha);
	    Dx_rot = floor(Dx_);
        Dy_rot = floor(Dy_);

        // Bilinear interpolation
        j1 = (int)Dx_rot + (numbinmap-1);
	    j2 = (int)Dx_rot + (numbinmap-1) + 1;
        k1 = (int)Dy_rot + (numbinmap-1);
	    k2 = (int)Dy_rot + (numbinmap-1) + 1;
      
        Q11 = ac[i][k1][j1];
        Q12 = ac[i][k2][j1];
        Q21 = ac[i][k1][j2];
        Q22 = ac[i][k2][j2];
	  
  	    ac_rot[k][j] = (1./(1.)*(1.)) * ( Q11*(Dx_rot+1-Dx_)*(Dy_rot+1-Dy_) + Q21*(Dx_-Dx_rot)*(Dy_rot+1-Dy_) + Q12*(Dx_rot+1-Dx_)*(Dy_-Dy_rot) + Q22*(Dx_-Dx_rot)*(Dy_-Dy_rot) );
        } // end if
	  }
    }

  return;
}

void compute_statistics(float score[numecunit]){
  int i, Position;
  float BinWidth;
 
  mean = 0.;
  for(i=0;i<numecunit;i++)
    mean += score[i];
  mean /= (float)numecunit;
      
  var = 0.;
  for(i=0;i<numecunit;i++)
    var += pow(score[i]-mean,2);
  var /= (float)(numecunit-1);

  for(i=0;i<Nbins_histo;i++){
    Histo[i] = 0.;
    }
    
  BinWidth = (CorrSup-CorrInf)/(float)Nbins_histo;
  for(i=0;i<numecunit;i++){
    Position = (int) floor((score[i]-CorrInf)/BinWidth);
    Histo[Position]++;
	}
  for(i=0;i<Nbins_histo;i++){
    Histo[i] /= (numecunit*BinWidth);
    }

  return;
}

//saves maps, autocorrelations and more
void save_all(){
  int i,j,k,l;
  float BinWidth;
  FILE *disp;
  
  for(i=0;i<numecunit;i++){
    sprintf(filename,"%s%ccell_%d%cweight-%d.dat",foldername, 47, i+1, 47, idstep);
    disp=fopen(filename,"w+");
    for(k=0;k<numhippunit;k++){
      fprintf(disp,"%d %.3f\n",k+1,weight_hipp2ec[i][k]);
      }
    fclose(disp);
		
    sprintf(filename,"%s%ccell_%d%cmap-%d.dat",foldername, 47, i+1, 47, idstep);
    disp=fopen(filename,"w+");
    for(l=0;l<numbinmap;l++){
      for(k=0;k<numbinmap;k++) fprintf(disp,"%.3f ",map_pos[i][k][l]); 
      fprintf(disp,"\n ");
      }
    fclose(disp);

    sprintf(filename,"%s%ccell_%d%cautoc-%d.dat",foldername, 47, i+1, 47, idstep);
    disp=fopen(filename,"w+");
    for(l=0;l<2*numbinmap-1;l++){
      for(k=0;k<2*numbinmap-1;k++) fprintf(disp,"%.3f ",auto_corre[i][l][k]);
      fprintf(disp,"\n ");
      }
    fclose(disp);

    sprintf(filename,"%s%ccell_%d%cmapang-%d.dat",foldername, 47, i+1, 47, idstep);
    disp=fopen(filename,"w+");
    for(k=0;k<numbinmap;k++) fprintf(disp,"%.3f ",map_ang[i][k]);
    fprintf(disp,"\n ");
    fclose(disp);
    } // end for

  BinWidth = (CorrSup-CorrInf)/(float)Nbins_histo;

  sprintf(filename,"%s%cGridScore_Histogram-%d.dat",foldername, 47,idstep);
  disp=fopen(filename,"w+");
  for(l=0;l<Nbins_histo;l++) fprintf(disp,"%f %.3f\n",CorrInf+(l+0.5)*BinWidth,Histo[l]);
  fclose(disp);

  sprintf(filename,"%s%cpath-x.dat",foldername, 47);
  disp=fopen(filename,"w+");
  for(l=0;l<lenhis;l++) fprintf(disp,"%d %.3f\n",l,x_list[l]);
  fclose(disp);

  sprintf(filename,"%s%cpath-y.dat",foldername, 47);
  disp=fopen(filename,"w+");
  for(l=0;l<lenhis;l++) fprintf(disp,"%d %.3f\n",l,y_list[l]);
  fclose(disp); 

  sprintf(filename,"%s%cpath-theta.dat",foldername, 47);
  disp=fopen(filename,"w+");
  for(l=0;l<lenhis;l++) fprintf(disp,"%d %.3f\n",l,theta_list[l]);
  fclose(disp);

  return;
}

//activation of a simple place cell like input neuron
float fhc(float x0, float y0){
  return exp(-(x0*x0+y0*y0)/(2.*spreadplacefield*spreadplacefield));
}

//activation of a simple place cell like input neuron
float fdir(float angle, float pref_angle, int n){
  return ( c0 + (1.-c0) * exp( nu[n]*( cos(pref_angle-angle) - 1. ) ) );
}

//Is this position inside or outside the box?
int isin(float x0, float y0){
  float radius;

  // circle of diameter environmentsize
  radius = sqrt(pow(x0-environmentsize/2.,2) + pow(y0-environmentsize/2.,2));  
  
  if(radius<=environmentsize/2.0) return 1;
  else return 0;
}

// the remainder (modulus b), even for negative values.
float modulus(float a, float b){
  while(a <0.) a += b;
  return fmod(a,b);
}
// normal distribution
float normal(){
  return sqrt(-2.*log((rand()+1.)/((float)RAND_MAX+1.)))*sin(2.*M_PI*rand()/((float)RAND_MAX+1.));
}

//uniform distribution in [0 1).
float uniform(){
  return ((float)rand())/((float)RAND_MAX + 1.0);
}
