#define environmentsize 200. //size of square room
#define constspeed 0.4 //speed of the rat
#define spreadturn 0.2 //angular head turn dispersion
#define numbinmap 50 //pixels per dimension in a map
#define map_perturb 0. //randomness in place fields.
#define lrpositionmap 0.01 //speed of map update
#define lranglemap 0.1*lrpositionmap
#define plotstep 250000
#define lenhis 250000 //remember past steps.
#define numecunit 600 //number of ec units
#define numside 50
#define Nrings 30
#define ext_sidehippgrid 190
#define numhippunit_aux Nrings*ext_sidehippgrid //number of place units
#define nconnhipp2ec_aux numhippunit_aux
#define spreadplacefield 5.0 //place field standard deviation
//#define b1 0.100 //adaptation parameters
//#define b2 b1/3.
#define b1_min 0.015
#define b1_max 0.300
#define B1 0.06293
#define b3 0.01
#define b4 0.1
#define lrmean 0.05 //averaging factor for firing rate (takes ~1000 steps).
#define lrweight1 0.005 //speed of weight learning
#define lrweight2 0.001 //speed of weight learning
#define sparsitylevel 0.3 //desired sparsity
#define meanactivitylevel 0.1 //desired mean
#define mingain 0.01
#define maxgain 50000.0
#define nummaxiter 1000 //max step in sparsity and mean control
#define tolerancemean 0.1 //relative error in spr mean control
#define firingconst 2.0/M_PI
#define c0 0.2
//define nu 0.8
#define spreadlateral 15. //lateral connectivity standard deviation
#define kappa 0.05
#define expect_length 10.
#define delay 25
//#define rho 0.2
#define Nbins 50
#define Nbins_histo 80
#define CorrSup 2.0
#define CorrInf -2.0
#define r_avoid 20.0 // to prevent close peaks in autocorrelogram
#define Rmax_min 45.0
#define Max(a, b) ((a) > (b) ? (a) : (b))
#define Min(a, b) ((a) < (b) ? (a) : (b))
// Functions
void initialize();
float fhc(float, float); // place field
float fdir(float angle, float pref_angle, int n); // head direction selectivity
void update_position();
void update_neurons(bool control_mean);
void networkoutput(bool control_mean);
void ec_output();
void fix_sparsity();
bool find_threshold_range(float max_act);
void search_threshold(int direction, float &back_bound, float &ahead_bound);
bool is_fixed();
void comput_gain_gradient();
void adapt_gain();
void update_maps();
void update_weights();
void update_mean();
void set_input();
void compute_autocorrelation(float ac[numecunit][2*numbinmap-1][2*numbinmap-1], float map[numecunit][numbinmap][numbinmap]);
void compute_gridness(float ac[numecunit][2*numbinmap-1][2*numbinmap-1], float score[numecunit], int flag_print);
void assess_rmin(float ac[numecunit][2*numbinmap-1][2*numbinmap-1]);
float assess_corr(int i, float ac[numecunit][2*numbinmap-1][2*numbinmap-1], float r1, float r2, int print);
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);
void compute_statistics(float score[numecunit]);
int isin(float, float); //is this position inside the box?
void make_dirs();
void save_all();
float normal(); //random normal generator
float uniform();
float modulus(float, float);
// Global variables
float x_place_cell[numhippunit_aux],y_place_cell[numhippunit_aux]; //centers of input distributions for HC (space)
float weight_hipp2ec[numecunit][nconnhipp2ec_aux]; //weights from hippocampus to entorhinal cortex
float act_ec[numecunit], act_hipp[numhippunit_aux]; //states of entorhinal and hippocampal neurons
float active_ec[numecunit],inactive_ec[numecunit],input2ec[numecunit]; //adaptation variables for MEC
float mean_act_ec[numecunit],mean_act_hipp[numhippunit_aux]; //temporal mean of EC and HC
float mean_network_act, mean_square_network_act, sparsity; //mean, square mean, sparsity EC activity
float threshold = 0, gain = 2; //threshold, gain of the trnasfer function.
float max_active_ec, min_active_ec; //range of active_ec.
int n_iter;
float inc_threshold = 0.005;
float upper_threshold = 1, lower_threshold = -1; //threshold bound.
float d_mean_gain; //mean gradient of activity w.r.t. gain.
float b1[numecunit], b2[numecunit], nu[numecunit], rho[numecunit];
float map_pos[numecunit][numbinmap][numbinmap]; //spatial map of EC cells
float map_ang[numecunit][numbinmap]; //angular map of EC cells
float auto_corre[numecunit][2*numbinmap-1][2*numbinmap-1]; //autocorrelatoin
float ac_rot60[2*numbinmap-1][2*numbinmap-1];
float ac_rot120[2*numbinmap-1][2*numbinmap-1];
float ac_rot30[2*numbinmap-1][2*numbinmap-1];
float ac_rot90[2*numbinmap-1][2*numbinmap-1];
float ac_rot150[2*numbinmap-1][2*numbinmap-1];
float theta, x, y; //head direction and position
float speed;
int x_map, y_map, theta_map; //variables to construct maps
float mean_list[lenhis]; //average activity
float sparsity_list[lenhis]; //sparsity
float x_list[lenhis], y_list[lenhis], theta_list[lenhis], speed_list[lenhis];
char foldername[500], filename[500];
int idstep;
int idhis;
float ConjCell_pref_HD[numecunit], ConjCell_X[numecunit], ConjCell_Y[numecunit];
float weight_lateral[numecunit][numecunit];
float delayed_act_ec[numecunit][delay];
float rmin[numecunit], rmin_extreme[numecunit], rmax[numecunit], grid_score[numecunit];
float Dist[2*numbinmap-1][2*numbinmap-1], Corr[2*numbinmap-1][2*numbinmap-1], rmin_aux[numecunit], corr_aux[numecunit];
float BinData[numecunit][Nbins];
int NData[Nbins];
int numhippunit, nconnhipp2ec;
float mean, var, Histo[Nbins_histo];