/**********************  INIT  *************************************/
void init() {
 INT i,j,k;

  dt    = 0.1;						// Integration time step [ms].
  sdt   = sqrt(0.1);				// Square root of the dt [ms^0.5].

  eta   = 1.;						// Learning rate

  // (1) Model neuron parameters (i.e. exponential I&F)
  th    = 20.;						// [mV] - peak value for a spike
  C     = 281.;     					// Membrane capacitance [pF]
  g_L   = 30.;						// Membrane conductance [nS] (i.e. taum = 16.6667 ms) 
  E_L   = -70.6;						// Resting membrane potential [mV]
  V_T   = -50.4;                    // Excitability threshold [mV].
  H     = -70.6;                    // Reset potential [mV].
  Delta_T = 2.;						// Excitability slope [mV]
  Tarp  = 2.;						// Absolute refractory period [ms].

 // (2) Model neuron parameters (i.e. spike-freq. adaptation of the exp I&F)
  tau_w   = 144.;					// [ms] - decay time constant for adaptation variable
  a       = 4.;						// [nS] - voltage-dependence of adaptation variable
  b       = 0.0805;					// [nA] - spike-dependence of adaptation variable

 // (3) Synaptic coupling (and plasticity) parameters
  maxweight	= 5.;					// was 1.2
  tau_Isyn= 5.;						// Decay time constant of individual EPSC [ms]
  //GA    = 200.;						// Maximal synaptic weight (for each synapse, same value) - was 200 for facilitating and depressing synapses
  fU    = 0.1;						// Probability of release (for each synapse, same value) 
  fD    = 100.;	// was 5					// [ms] recovery from short-term depression (for each synapse, same value)
  fF    = 900.; // was 5000					// [ms] recovery from short-term facilitation (for each synapse, same value)
 
  dU    = 0.8;						// Probability of release (for each synapse, same value) 
  dD    = 900.;	// was 500					// [ms] recovery from short-term depression (for each synapse, same value)
  dF    = 100.;	// was 5					// [ms] recovery from short-term facilitation (for each synapse, same value)

  rate  = 0.;  
  Trate = 5000.;                     // [ms] temporal window to compute mean firing frequencies over
    
 // -- VISUAL CORTEX, MINIMAL MODEL (Pfister & Gerstner, 2006) ------------------------------
  taur2    = 101.;           // Decay time constant of presynaptic indicators (tau_x)
  tauo2    = 114.;           // Decay time constant of postsynaptic indicators (tau_y)
  taur1    = 16.8;           // Decay time constant of presynaptic indicators (tau^+)
  tauo1    = 33.7;           // Decay time constant of postsynaptic indicators (tau^-)
	
  A2p      = 0;				// Amplitude of weight change, in case of a pre-post
  A3p      = 6.5e-3;		// Amplitude of weight change (Triplet-term), in case of a pre-post
  A2m      = 7.1e-3;		// Amplitude of weight change, in case of a post-pre
  A3m      = 0;				// Amplitude of weight change (Triplet-term), in case of a post-pre
 //-------------------------------------------------------------------------------------------	
  tmp_Ge		=	exp(-dt/tau_Isyn);

  tmp_r1		=	exp(-dt/taur1);
  tmp_r2		=	exp(-dt/taur2);
  tmp_o1		=	exp(-dt/tauo1);
  tmp_o2		=	exp(-dt/tauo2);
  
  tmp_u1     =   (dt / C) * g_L * E_L;
  tmp_u2     =   (dt / C) * g_L;
  tmp_u3     =   (dt / C) * g_L * Delta_T;
  tmp_u4     =   (dt / C);
  
//  mu    = 0.;
//  sigma = 0.;
  taui  = 5.;

  tmpnoise1 = (1. - dt/taui);
  //tmpnoise2 = dt*mu/taui;
  tmpnoise3 = sigma * sqrt(2 * dt/taui);


for(i=0;i<N;i++) {
	to[i]    = -9999.;                  // Initialization of the last firing times (to -infty)
    to_old[i]= -99999.;                 // Initialization of the last firing times (to -infty)
    u[i]     = E_L +(V_T - E_L) * 0*drand49(); // Random initialization of the membrane voltages.
    w[i]     = 0.;						// Spike-frequency adaptation variable [pA]
    Ge[i]    = 0.;
    mu[i]    = Io * (1 + 1 * gauss());
    //mu[i]    = Io;
    Iext[i]  = mu[i];
	r1[i]    = 0.;						// Presynaptic indicator, one for each neuron
	r2[i]    = 0.;						// Presynaptic indicator, one for each neuron
	o1[i]    = 0.;						// Postsynaptic indicator, one for each neuron
	o2[i]    = 0.;						// Postsynaptic indicator, one for each neuron
    
    rates[i] = 0.;
    }
    
    for(i=0; i<N; i++)
     for(j=0; j<N; j++) {
     
         if ((i<(N/2)) & (j<(N/2)))
             W[i][j] = drand49() * maxweight;
         else if ((i>(N/2)) & (j>(N/2)))
             W[i][j] = drand49() * maxweight;
         else
             W[i][j] = 0.1* drand49() * maxweight;
       
     //    W[i][j] = drand49() * maxweight;
      
        if (CC[j][i] != 0) // The current neuron may receive an EPSC..
            if (CC[j][i] == 1)   // Facilitatory excitatory synapse
                    {
                     Gr[j][i] = 1.;
                     Gu[j][i] = fU;
                    }
            if (CC[j][i] == -1)   // Depressing excitatory synapse
                    {
                     Gr[j][i] = 1.;
                     Gu[j][i] = dU;
                    }
     }

 Nspiked    = 0;			// ???		
 return;
}// end init()



int read_connectivity() {		// UPDATED JULY 2011 !!!!!!!!!!!!!!!!!!!
	// ..reads 'connectivity.dat' (binary) 
	// short **CC; <-- global variable
	// INT N;      <-- global variable
	INT i,j,k,out,MM;					// Useful temporary variables (e.g. indexes, return values, etc.)
	FILE *fopen(), *fp;					// Needed for file I/O operations...
	char fname[20];						// This will contain the file name for connectivity...
	double Ncells, tmp;					// This will contain the number of cells defined in the file...
	
	sprintf(fname, "connectivity.dat");	// This is the file generated by the MATLAB GUI [FIXED NAME]
    //NOTE:
    // This file has been generated in MATLAB by the following statements:
    //MM    = length(find(C(:)~=0));  % Number of non-zero elements of C(i,j)
    //Ctemp = C';
    //Ctemp = Ctemp(:);   % A single vector, encoding 'C' row-wise
    //tmp   = find(Ctemp ~= 0); // Only the non-zero elements of C are considered to be written on disk
    //fname = 'connectivity.dat';
    //fp    = fopen(fname, 'wb');
    //fwrite(fp, Ncell_exc, 'double');  
    //fwrite(fp, MM, 'double');  
    //for kk=1:length(tmp)
    //    fwrite(fp, tmp(kk) * sign(Ctemp(tmp(kk))), 'double');
    //end
    //fclose(fp);
	//------------------------------------------------------------------------------------------- 
	printf("Opening connectivity file <%s> (generated by MATLAB-GUI)...\n", fname);
	fp = fopen(fname, "r");				// Try to open file 'fname' for reading...
	if (fp==NULL) {  printf("Error: unable to open the file..\n");  return -1; }
	
	out = fread(&Ncells, sizeof(double), 1, fp); 	// I expect the first element to be 'Ncells'..
	if (out!=1) {  printf("Error: corrupted or wrong format for the connectivity file - first element is NOT Ncells!\n");  return -1; }
	
	N  = (INT) Ncells;								// First element stored as the global variable 'N'.
	printf("Found connectivity information for %d cells!\n", (int) N);
	
	out = fread(&tmp, sizeof(double), 1, fp); 	// I expect the second element to be 'MM'..
	if (out!=1) {  printf("Error: corrupted or wrong format for the connectivity file! - second element is NOT MM\n");  return -1; }
	MM = (INT) tmp;
	
	//------------------------------------------------------------------------------------------- 
	printf("Allocating memory [CC and W] for this simulation..\n");
	CC  = calloc(N, sizeof(short*));			// Reserving memory for the N x N int connectivity matrix
	if (CC==NULL) {  printf("Error: unable to allocate memory [CC]..\n"); fflush(NULL); return -1; }

	W  = calloc(N, sizeof(double*));				// Reserving memory for the N x N double synaptic matrix
	if (W==NULL) {  printf("Error: unable to allocate memory {W]..\n"); fflush(NULL); return -1; }
	
	for (i=0; i<N; i++) {
		CC[i] = calloc(N, sizeof(short)); // Reserving memory for the N x N int connectivity matrix
		if (CC[i]==NULL) {  printf("Error: unable to allocate memory [CCi]..\n"); fflush(NULL); return -1; }  

		W[i] = calloc(N, sizeof(double)); // Reserving memory for the N x N double synaptic matrix
		if (W[i]==NULL) {  printf("Error: unable to allocate memory [Wi]..\n"); fflush(NULL); return -1; }  
		
		//out = fread(CC[i], sizeof(short), N, fp);
		//if (out!=N) {  printf("Error: corrupted or wrong format for the connectivity file!\n", out);  return -1; }  
	}

    for (k=0; k<MM; k++) {
     out = fread(&tmp, sizeof(double), 1, fp); 	// I expect the second element to be 'MM'..
	 if (out!=1) {  printf("Error: corrupted or wrong format for the connectivity file!\n"); fflush(NULL); return -1; }
      
        
      i = (INT) ((fabs(tmp)-1)/N);  
      j = (INT) fmod(fabs(tmp)-1,N); 
      CC[i][j] = (short) ((tmp>0) ? 1 : -1);
    }
        printf("\n");
	// PRINT THE CONNECTIVITY MATRIX FOR DEBUGGING PURPOUSES
	/*
	for (i=0;i<N;i++) {
	 for (j=0; j<N; j++)
	 	printf("%d ", CC[i][j]);
	 printf("\n");
	}
	*/
	printf("Connectivity file <%s> acquired successfully!\n\n\n", fname);
	fclose(fp);									// File is closed here.
	//------------------------------------------------------------------------------------------- 
	
	return 0;
} // end read_connectivity()


void print_output() {
    if (fmod(t, Trate)<dt) {
        rate = 1000. * rate / (N * Trate);
        printf("Mean firing rate: %.0f Hz\r", rate);
        frate = fopen("rates.x", "w");
        for (i=0;i<N;i++) {
            fprintf(frate, "%f\n", 1000. * rates[i] / Trate);
            rates[i] = 0.;
        }
        fclose(frate);
    }
    
    if (fmod(t, 10*Trate)<dt) {
        gfile  = fopen("G.x","w");
        for (i=0;i<N;i++) {
            for (j=0;j<N;j++)
                fprintf(gfile, "%f ", W[i][j]);
            fprintf(gfile, "\n");
        }
        fclose(gfile);
    }

/*    
    //  if(fmod(t,20*dt)<=dt)   // Every 20 simulation steps, the membrane voltage is wrote on disk.
    fprintf(sample,"%f %f %f %f %f %f %f\n",t, u[0], Ge[0], u[1], Ge[1], u[2], Ge[2]);
    //if(fmod(t,100*dt)<dt) fprintf(gfile,"%f %f %f %f %f %f %f\n",t, W[0][1], W[1][0], W[3][5], W[5][3], W[7][4], W[4][7]);
 */
    return;
}


void log_weights() {
	INT i,j,k,out,MM;					// Useful temporary variables (e.g. indexes, return values, etc.)
	FILE *fopen();
	char fname[20];						// This will contain the file name for connectivity...
	
	sprintf(fname, "W_%.0f.x", t);	

	weightmatrix = fopen(fname, "w");				 // Open output file for writing...
    for (i=0; i<N; i++) {
     for (j=0; j<N; j++) 
	fprintf(weightmatrix, "%f ", W[i][j]);
	fprintf(weightmatrix, "\n");
	}
return;
} // end log_weights()




/*
int read_pars() {		// UPDATED JULY 2011 !!!!!!!!!!!!!!!!!!!
	//
	// ..reads 'pars.dat' (binary)
	//
	//[Ncell_exc, Prob_ee, Prob_loop, Prob_ff, Prob_dd, Prob_fd, Prob_f, Prob_d, fixseed, seed];
	//
	// INT N;      <-- global
	
	FILE *fopen(), *fp;					// Needed for file I/O operations...
	char fname[20];						// This will contain the file name for parameters...
	INT i, out;							// Useful temporary variables (e.g. indexes, return values, etc.) 
	double par;
	
	sprintf(fname, "pars.dat"); // This is the file generated by the MATLAB GUI [FIXED NAME]
	// NOTE: 
	//This file has been created by the following MATLAB-statements
    //
	//fname = 'pars.dat';
	//fp    = fopen(fname, 'wb');
	//all_pars = [Ncell_exc, Prob_ee, Prob_loop, Prob_ff, Prob_dd, Prob_fd, Prob_f, Prob_d, fixseed, seed];
	//fwrite(fp, pars, 'double');  
	//fclose(fp);
	//
	//We aim at extracting only two elements: the 9th "fixseed" and the 10th "seed".

	
	//------------------------------------------------------------------------------------------- 
	printf("Opening parameter file <%s>..\n", fname);
	fp = fopen(fname, "r");				// Try to open file 'fname' for reading...
	if (fp==NULL) {  printf("Error: unable to open the file..\n");  return -1; }
	
	// The first 8 elementE_Ls are inessential at run-time for the simulation...
	for (i=0; i<8; i++) {							// Let's go through each of them...
		out = fread(&par, sizeof(double), 1, fp);	// and simply ignore it, while the file-access pointer moves forward...
		if (out!=1) {  printf("Error: unable to read connectivity file headers..\n");  return -1; }
	}
    
	out = fread(&par, sizeof(double), 1, fp);  fixseed = par;
	if (out!=1) {  printf("Error: unable to read connectivity file headers..\n");  return -1; }
	out = fread(&par, sizeof(double), 1, fp);  seed = par;
	if (out!=1) {  printf("Error: unable to read connectivity file headers..\n");  return -1; }
	
	printf("Parameter file <%s> acquired successfully!\n\n\n", fname);
	fclose(fp);									// File is closed here.
	//------------------------------------------------------------------------------------------- 
	
	return 0;
} // end read_pars()
*/

/*

void manage_synapses(INT i) {	// UPDATED JULY 2011 !!!!!!!!!!!!!!!!!!!
	// Neuron i-th is examined, ready to receive inputs from its presynaptic peers (j=1..N), if any...
	double delta;
	for(j=0;j<N;j++)     // Now, all the connected presynaptic neurons are checked for (recent) spiking activity.
	  if (CC[i][j] != 0) // A connection exists from neuron j-th (pre) to neuron i-th (post)!
         if ((t-(to[j]+latency))<dt && (t-(to[j]+latency))>=0) // Indeed, presynaptic neuron j-th fired - after latency of EPSPs...
            { 
            delta   = (to[j]-to_old[j]);
            if (CC[i][j] == 1)   // Facilitatory exciting synapse
                    {
                     Gr[i][j] = (Gr[i][j] * (1 - Gu[i][j]) - 1) * exp(-delta/fD) + 1;
                     Gu[i][j] = Gu[i][j] * (1 - fU) * exp(-delta/fF) + fU;
                    }
            else // (CC[i][j] == -1)   // Depressing exciting synapse
                    {
					Gr[i][j] = (Gr[i][j] * (1 - Gu[i][j]) - 1) * exp(-delta/dD) + 1;
					Gu[i][j] = Gu[i][j] * (1 - dU) * exp(-delta/dF) + dU;
                    }
			Ge[i]  += GA *Gr[j][i] * Gu[j][i];
            }
return;
} // end manage_synapses()




void manage_stdp(INT j) {		// UPDATED JULY 2011 !!!!!!!!!!!!!!!!!!!
   // Neuron j-th just fired...
   INT i;
   for(i=0;i<N;i++)    // Now, all the connected neurons, postsynaptic to the j-th must be updated..
        {
			W[i][j] -= eta * o1[i] * abs(CC[i][j]) * (A2m + A3m * r2[j]); // the j-th neuron is presynaptic to all (connected) generic neurons i-th
			W[j][i] += eta * r1[i] * abs(CC[j][i]) * (A2p + A3p * o2[j]); // the j-th neuron is postsynaptic to all (connected) generic neurons i-th
		} 
	
   r1[j]++;			// update the pre/post-synaptic detectors of the neuron j-th, which just fired
   r2[j]++;			// update the pre/post-synaptic detectors of the neuron j-th, which just fired 
   o1[j]++;			// update the pre/post-synaptic detectors of the neuron j-th, which just fired
   o2[j]++;			// update the pre/post-synaptic detectors of the neuron j-th, which just fired
	return;
} // end manage_stdp()
*/
