/* Simple spiking neuron according to Izhikevich 
   for the dentate granule cell

Nearest-neighbour additive STDP
BCM sliding threshold
max and min weight limit

Copyright by Lubica Benuskova

*/

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<sys/time.h>

typedef struct _neuron *neuron;

neuron make_neuron_struct();
void init_neuron_memory(neuron n);
neuron make_neuron();
void spontaneous_input(int time, neuron single_cell);
void test_input(int time, neuron single_cell);
void Med50(int time, neuron single_cell);
void Lat50(int time, neuron single_cell);
void cycle_cell(int time, neuron single_cell);
void stdp_train_cell(neuron single_cell);
void bcm_threshold(int time, neuron single_cell);
void test_sample(int time, neuron single_cell);
void run_network(int time);
void getTime(long *z);

#define SIMULATION_PERIOD 25200000  // total # of time steps in the simulation
#define SPONT_PERIOD       3600000  // spontaneous period before MPP HFS
#define MED_ONSET          5400000  // onset of med HFS
#define LAT_ONSET         19800000  // onset of lat HFS
#define MIN30              1800000  // 30 min in ms
#define MIN10               600000    // 10 min in ms
#define MIN1                 60000   // 1 min in ms

#define TIME_STEP                1    // time step in ms
#define SAMPLE_PERIOD        60000    // time step between samples
#define minutes              60000.0  // to get time in minutes
#define TEST_PERIOD          10000    // interval between test pulses in ms

#define TEST_INTENSITY         150   // number of input fibers
#define TRAIN_INTENSITY        250   // number of input fibers
#define AP                      55.0   // amplitude of action potential
#define P_HFS                    0.400   // probability of mpp firing during HFS
#define P_theta                  0.008 // probability of pp synchronous spont. firing
#define P_pp                     0.0001 // probability of pp asynchronous spont. firing
#define MAX                  10000 // max number of presynaptic spikes between 2 postsynaptic
#define CW                       5  // correlation window for inputs in ms

#define TAU                  60000    // Integration period for theta_M
#define BOLTZMAN                 0  	// should a Boltzman average be used
#define theta_M0             2000.0  // thetaM is scaled by this parameter 
#define MAX_LOOP               1    // number of measurments

//FILE *f_gc_voltage,*f_lpp_input,*f_mpp_input;
FILE *f_mpp_weight,*f_lpp_weight;
//FILE *f_pre_spikes;
FILE *f_theta_M,*f_mpp_ltp,*f_mpp_ltd;
/**************************************************************************
		Variables
***************************************************************************/

/* weight initialization. (Initial bias) */
float MPP_WEIGHT_INIT = 0.05;
float LPP_WEIGHT_INIT = 0.05;
float firing_theta = 24.0;  // excitation threshold in mV

/* STDP rule parameters */
float tau_ltp = 20;  // LTP tau in ms 
float tau_ltd = 100; // LTD tau in ms 
float A0_ltp = 0.02;    // LTP relative change x current weight 
float A0_ltd = 0.01;    // LTD relative change x current weight 

/* Izhikevich's neuron parameters values for RS cell*/
float coeff1 = 0.04;
float coeff2 = 5;
float coeff3 = 140;
float coeff4 = 1;

float a = 0.02;
float b = 0.2;
float c = -69.0;
float d = 2;

/* general variables */
int mpp_input, lpp_input, t1_gc, t2_gc, index_lpp, index_mpp;
int t_lpp[MAX],t_mpp[MAX];
int last_mpp_test, last_lpp_test, last_sample;
int N_MPP,N_LPP; // number of input fibers
int loop, k, burst, train, impulse;
float input, rnd, gc_av_frequency, mpp_thetaM[MAX],lpp_thetaM[MAX];
float av_mpp_weight[422],av_lpp_weight[422],av_thetaM[422];
/**************************************************************************
		End of Variables and Parameters
***************************************************************************/

/***********************************************************************
		Network structures
***********************************************************************/

struct _neuron
{
/* variables describing the neuron  */
        float v;      // total membrane potential in mV 
        float u;      // membrane recovery variable
        int output;   // binary output
/* weights of neuron */
        float mpp_weight; // mpp input weight 
        float lpp_weight; // lpp input weight 
/* BCM averaging stuff */
        float theta_M;	// dynamic modification LTD-LTP threshold
        float mem;      // average of activity
        float pf; 	// used to calculate mem
        float boltzfactor;	// used to calculate mem
}; /*end struct _neuron */

/************************************************************************
			Neuron Stuff
*************************************************************************/

/********* This function creates a new neuron initializes it to 0 ******/
                   /* It is added to the neuron list.*/
neuron make_neuron_struct()
{
    neuron n;
 
    n = (neuron) malloc(sizeof(struct _neuron));

    n->v = c;
    n->u = b * c; 
    n->output = 0;	

    n->mpp_weight = MPP_WEIGHT_INIT;
    n->lpp_weight = LPP_WEIGHT_INIT; 
   
    n->theta_M = 1.0; 
    n->mem = 0.0;
    n->pf = 1.0;
    n->boltzfactor = 0.0; 

    return n;
} /*end neuron make_neuron_struct()*/

void init_neuron_memory(neuron n)
{
    n->boltzfactor = exp( -((double) TIME_STEP / (double) TAU));
    if(BOLTZMAN == 0)
	n->pf = 1.0 / (1 - n->boltzfactor);
    else
	n->pf = 1.0;
} /* end void init_neuron(n)*/

/******** This function creates and initializes a new neuron **********/
neuron make_neuron()
{
    neuron n;

    n = make_neuron_struct();
    init_neuron_memory(n);

    return n;
} /* end neuron make_neuron()*/

/************************************************************************
			End Neuron Stuff
*************************************************************************/

void main()
{ 
  long z;
  getTime(&z);
  srandom(z);

   for(loop = 0; loop <= 420; ++loop) 
    {av_mpp_weight[loop] = 0.0; av_lpp_weight[loop] = 0.0; av_thetaM[loop] = 0.0;}  

   for(loop = 0; loop < MAX_LOOP; ++loop) run_network(SIMULATION_PERIOD);

  /* f_mpp_input = fopen("mpp_pulse_times.dat","w");
  f_lpp_input = fopen("lpp_pulse_times.dat","w");
  f_gc_voltage = fopen("gc_voltage.dat","w"); 
  f_pre_spikes = fopen("pre_spikes.dat","w");
  f_mpp_ltp = fopen("mpp_ltp.dat","w");
  f_mpp_ltd = fopen("mpp_ltd.dat","w");*/

  f_mpp_weight = fopen("mpp_weight.dat","w");
  f_lpp_weight = fopen("lpp_weight.dat","w");
  f_theta_M = fopen("theta_M.dat","w");
 
  for(loop = 0; loop <= 420; ++loop) {
    fprintf(f_mpp_weight,"%d\t%6.3f\n",loop-60,av_mpp_weight[loop] / MAX_LOOP);
    fprintf(f_lpp_weight,"%d\t%6.3f\n",loop-60,av_lpp_weight[loop] / MAX_LOOP);
    fprintf(f_theta_M,"%d\t%f\n",loop-60,av_thetaM[loop] / MAX_LOOP);
  }   

  /*fclose(f_mpp_input);
  fclose(f_lpp_input);
  fclose(f_gc_voltage);
  fclose(f_pre_spikes);
  fclose(f_mpp_ltp);
  fclose(f_mpp_ltd);*/

  fclose(f_mpp_weight);
  fclose(f_lpp_weight);
  fclose(f_theta_M);

      
} /* end void main()*/

/************************************************************************
			Network Stuff
*************************************************************************/

void run_network(int t)
{
    int time, i; 
    neuron single_cell;
    single_cell = make_neuron();
            
    for(i = 0; i < MAX; ++i) 
    {t_mpp[i] = 0; t_lpp[i] = 0; mpp_thetaM[i] = 0.0; lpp_thetaM[i] = 0.0;}
    gc_av_frequency = 0;
    last_mpp_test =0;
    last_lpp_test =0;
    t1_gc = 0;
    t2_gc = 0;
    index_lpp = 0;
    index_mpp = 0;
    last_sample = 0;
    burst = 1;
    train = 1;
    impulse = 1;
    mpp_input =0;
    lpp_input = 0;
    k = 0;

    N_MPP = TRAIN_INTENSITY; N_LPP = TRAIN_INTENSITY;

    time = 0;
    while(time <= t)
    {
     if(
     (time < MED_ONSET) || 
     ((time > (MED_ONSET+MIN10)) && (time < LAT_ONSET)) || 
     (time > (LAT_ONSET+MIN10))
     ) 
        spontaneous_input(time, single_cell);

     if(
     ((time > SPONT_PERIOD) && (time < MED_ONSET)) || 
     ((time > (MED_ONSET+MIN10)) && (time < LAT_ONSET)) ||
     (time > (LAT_ONSET+MIN10))
     ) 
        test_input(time, single_cell);

     if((time >= MED_ONSET) && (time <= (MED_ONSET+MIN10))) Med50(time, single_cell);
       
     if((time >= LAT_ONSET) && (time <= (LAT_ONSET+MIN10))) Lat50(time, single_cell);

     cycle_cell(time, single_cell);
     
     if((time == t2_gc) && (t1_gc > 0)) stdp_train_cell(single_cell); 

     bcm_threshold(time, single_cell);
     test_sample(time, single_cell); 
     time = time + TIME_STEP;      
    } /* while */ 

    gc_av_frequency = gc_av_frequency / ((float)SIMULATION_PERIOD / 1000.0);
    printf("Average output frequency\t");
    printf("%f\n",gc_av_frequency);
        
}  /* end  void run_network()*/ 

/************************************************************************
This function generates a spontaneous input
************************************************************************/
void spontaneous_input(int time, neuron n)
{
    float rnd;
    int ro;
    
    mpp_input = 0; lpp_input = 0;

// correlated inputs
    rnd = (float)(random())/(float)RAND_MAX;   
    if(P_theta > rnd) 
       {
        mpp_input = 1;  N_MPP = TEST_INTENSITY;
        lpp_input = 1;  N_LPP = TEST_INTENSITY;
        //fprintf(f_mpp_input,"%f\t%d\n",(float)time/minutes,mpp_input); 
        //fprintf(f_lpp_input,"%f\t%d\n",(float)time/minutes,lpp_input); 

        if(t1_gc > 0)
        { if((index_mpp < MAX) && (index_lpp < MAX) && (time > t2_gc))
	  { t_mpp[index_mpp] = time;
	    t_lpp[index_lpp] = time; 
            mpp_thetaM[index_mpp] = n->theta_M; 
            lpp_thetaM[index_lpp] = n->theta_M;            
            index_mpp = index_mpp + 1;
            index_lpp = index_lpp + 1;
          }
	}
       }

// asynchronous spontaneous LPP input 
   rnd = (float)(random())/(float)RAND_MAX;   
   if((lpp_input == 0) && (P_pp > rnd) )
      {
       lpp_input = 1;  N_LPP = TEST_INTENSITY;
       //fprintf(f_lpp_input,"%f\t%d\n",(float)time/minutes,lpp_input);       

       if(t1_gc > 0)
        { if((index_lpp < MAX) && (time > t2_gc))
	  {t_lpp[index_lpp] = time; 
           lpp_thetaM[index_lpp] = n->theta_M;
           index_lpp = index_lpp + 1;
          }
	}
       }  
 // asynchronous spontaneous MPP input
   rnd = (float)(random())/(float)RAND_MAX;   
   if((mpp_input == 0) && (P_pp > rnd)) 
      {
       mpp_input = 1;  N_MPP = TEST_INTENSITY;
       //fprintf(f_mpp_input,"%f\t%d\n",(float)time/minutes,mpp_input);       

       if(t1_gc > 0)
        { if((index_mpp < MAX) && (time > t2_gc))
	  {t_mpp[index_mpp] = time; 
           mpp_thetaM[index_mpp] = n->theta_M;
           index_mpp = index_mpp + 1;
          }
	}
       }                

}  /* end void spontaneous_input(time) */


/************************************************************************
This function generates a test input
************************************************************************/
void test_input(int time, neuron n)
{
 
  if((time >= (SPONT_PERIOD+10000)) && (time - last_mpp_test >= 2*TEST_PERIOD))
   {
    last_mpp_test = time;
    if(mpp_input == 0)
    {
     mpp_input = 1; N_MPP = TEST_INTENSITY;
     //fprintf(f_mpp_input,"%f\t%d\n",(float)time/minutes,mpp_input); 
    
        if(t1_gc > 0)
        { if((index_mpp < MAX) && (time > t2_gc))
	  { t_mpp[index_mpp] = time;
            mpp_thetaM[index_mpp] = n->theta_M;           
            index_mpp = index_mpp + 1;
          }
	}
     }
   }

  if((time >= (SPONT_PERIOD+20000)) && (time - last_lpp_test >= 2*TEST_PERIOD))
  {
   last_lpp_test = time;
   if(lpp_input == 0)
    {
     lpp_input = 1; N_LPP = TEST_INTENSITY;
     //fprintf(f_lpp_input,"%f\t%d\n",(float)time/minutes,lpp_input); 
    
        if(t1_gc > 0)
        { if((index_lpp < MAX) && (time > t2_gc))
	  { t_lpp[index_lpp] = time;
            lpp_thetaM[index_lpp] = n->theta_M;           
            index_lpp = index_lpp + 1;
          }
	}
    }
  }

} /* end void test_input(time, n) */

/************************************************************************
This function generates a MPP CS consisting of 50 HFS trains
************************************************************************/
void Med50(int time, neuron n)
{   
  float rnd;

  mpp_input = 0; lpp_input = 0;

// 10 HFS bursts of 5 trains to MPP
if( 
(time >= (MED_ONSET + (burst-1)*MIN1 + (train-1)*1000 + (train-1)*25)) && 
(time <= (MED_ONSET + (burst-1)*MIN1 + (train-1)*1000 + (train-1)*25 + 25)) 
)
{
     rnd = (float)(random())/(float)RAND_MAX;   
     if(P_HFS > rnd)
     { 
      mpp_input = 1; N_MPP = TRAIN_INTENSITY;
      //fprintf(f_mpp_input,"%f\t%d\n",(float)time/minutes,mpp_input); 
    
        if(t1_gc > 0)
        { if((index_mpp < MAX) && (time > t2_gc))
	  { t_mpp[index_mpp] = time;
            mpp_thetaM[index_mpp] = n->theta_M;           
            index_mpp = index_mpp + 1;
          }
	}
      }

       if(impulse == 26) {train = train+1; impulse = 1;}
       if(train == 6) {burst = burst+1; train = 1;} 
       impulse = impulse + 1;

 // spontaneous LPP input during the bursts
   rnd = (float)(random())/(float)RAND_MAX;   
    if((P_theta+P_pp) > rnd) 
      {
       lpp_input = 1; N_LPP = TEST_INTENSITY;
       //fprintf(f_lpp_input,"%f\t%d\n",(float)time/minutes,lpp_input);       

       if(t1_gc > 0)
        { if((index_lpp < MAX) && (time > t2_gc))
	  {t_lpp[index_lpp] = time; 
           lpp_thetaM[index_lpp] = n->theta_M;
           index_lpp = index_lpp + 1;
          }
	}
       }              

}
else
{
// spontaneous LPP input in between the bursts
   rnd = (float)(random())/(float)RAND_MAX;   
    if((P_theta+P_pp) > rnd) 
      {
       lpp_input = 1; N_LPP = TEST_INTENSITY;
       //fprintf(f_lpp_input,"%f\t%d\n",(float)time/minutes,lpp_input);       

       if(t1_gc > 0)
        { if((index_lpp < MAX) && (time > t2_gc))
	  {t_lpp[index_lpp] = time; 
           lpp_thetaM[index_lpp] = n->theta_M;
           index_lpp = index_lpp + 1;
          }
	}
       }  
 // spontaneous MPP input in between the bursts
   rnd = (float)(random())/(float)RAND_MAX;   
    if((P_theta+P_pp) > rnd) 
      {
       mpp_input = 1; N_MPP = TEST_INTENSITY;
       //fprintf(f_mpp_input,"%f\t%d\n",(float)time/minutes,mpp_input);       

       if(t1_gc > 0)
        { if((index_mpp < MAX) && (time > t2_gc))
	  {t_mpp[index_mpp] = time; 
           mpp_thetaM[index_mpp] = n->theta_M;
           index_mpp = index_mpp + 1;
          }
	}
       }                
     
    test_input(time, n);        
}

      if(time == (MED_ONSET+MIN10)) 
      {
      last_mpp_test = MED_ONSET+MIN10-10000;
      last_lpp_test = MED_ONSET+MIN10;
      burst = 1; train = 1; impulse = 1;
      }

} /* end Med50(time, n) */


/************************************************************************
This function generates 50 HFS LPP trains 
************************************************************************/
void Lat50(int time, neuron n)
{   
  float rnd;

  mpp_input = 0; lpp_input = 0;

// 10 HFS bursts of 5 trains to LPP
if( 
(time >= (LAT_ONSET + (burst-1)*MIN1 + (train-1)*1000 + (train-1)*25)) && 
(time <= (LAT_ONSET + (burst-1)*MIN1 + (train-1)*1000 + (train-1)*25 + 25)) 
)
{
     rnd = (float)(random())/(float)RAND_MAX;   
     if(P_HFS > rnd)
     { 
      lpp_input = 1; N_LPP = TRAIN_INTENSITY;
      //fprintf(f_lpp_input,"%f\t%d\n",(float)time/minutes,lpp_input); 
    
        if(t1_gc > 0)
        { if((index_lpp < MAX) && (time > t2_gc))
	  { t_lpp[index_lpp] = time;
            lpp_thetaM[index_lpp] = n->theta_M;           
            index_lpp = index_lpp + 1;
          }
	}
      }

       if(impulse == 26) {train = train+1; impulse = 1;}
       if(train == 6) {burst = burst+1; train = 1;} 
       impulse = impulse + 1;

 // spontaneous MPP input during the bursts
   rnd = (float)(random())/(float)RAND_MAX;   
    if((P_theta+P_pp) > rnd) 
      {
       mpp_input = 1;  N_MPP = TEST_INTENSITY;
       //fprintf(f_mpp_input,"%f\t%d\n",(float)time/minutes,mpp_input);       

       if(t1_gc > 0)
        { if((index_mpp < MAX) && (time > t2_gc))
	  {t_mpp[index_mpp] = time; 
           mpp_thetaM[index_mpp] = n->theta_M;
           index_mpp = index_mpp + 1;
          }
	}
       }              

}
else
{
// spontaneous LPP input in between the bursts
   rnd = (float)(random())/(float)RAND_MAX;   
    if((P_theta+P_pp) > rnd) 
      {
       lpp_input = 1;  N_LPP = TEST_INTENSITY;
       //fprintf(f_lpp_input,"%f\t%d\n",(float)time/minutes,lpp_input);       

       if(t1_gc > 0)
        { if((index_lpp < MAX) && (time > t2_gc))
	  {t_lpp[index_lpp] = time; 
           lpp_thetaM[index_lpp] = n->theta_M;
           index_lpp = index_lpp + 1;
          }
	}
       }  
 // spontaneous MPP input in between the bursts
   rnd = (float)(random())/(float)RAND_MAX;   
    if((P_theta+P_pp) > rnd) 
      {
       mpp_input = 1;  N_MPP = TEST_INTENSITY;
       //fprintf(f_mpp_input,"%f\t%d\n",(float)time/minutes,mpp_input);       

       if(t1_gc > 0)
        { if((index_mpp < MAX) && (time > t2_gc))
	  {t_mpp[index_mpp] = time; 
           mpp_thetaM[index_mpp] = n->theta_M;
           index_mpp = index_mpp + 1;
          }
	}
       }                
     
    test_input(time, n);        
}

      if(time == (LAT_ONSET+MIN10)) 
      {
      last_mpp_test = LAT_ONSET+MIN10-10000;
      last_lpp_test = LAT_ONSET+MIN10;
      }


} /* end Lat50(time, n) */


/**********************************************************************
This function updates the neuron output according to its input
**********************************************************************/
void cycle_cell(int time, neuron n)
{   
  n->output = 0;

   if(n->v >= firing_theta) 
     { 
       n->v = c; n->u = n->u + d;    
     }
   else 
     {input =  mpp_input * n->mpp_weight * N_MPP + lpp_input *  n->lpp_weight * N_LPP;
      n->v = n->v + coeff1*(n->v * n->v) + coeff2*n->v + coeff3 - coeff4 * n->u + input;

      if(n->v >= firing_theta) {
          n->v = AP; 
          n->output = 1;
          if(t1_gc == 0) {t1_gc = time; t2_gc = time;}
          else {t1_gc = t2_gc; t2_gc = time;}  
          gc_av_frequency = gc_av_frequency + 1;             
          }
      n->u = n->u + a*(b * n->v - n->u); 
     }
         
   //fprintf(f_pre_spikes,"%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\n",t1_gc,t2_gc,mpp_input,index_mpp,t_mpp[index_mpp-1],lpp_input,index_lpp,t_lpp[index_lpp-1]); 

                                                                       
}  /* end void cycle_cell(time, n) */

/********************************************************************** 
                   STDP training
**********************************************************************/
void stdp_train_cell(neuron n)
{
  float ltp, ltd, A_ltp, A_ltd;
  int k;

  index_lpp = 0; index_mpp = 0; 

  while(index_mpp < MAX)
    {if(t_mpp[index_mpp] == 0) index_mpp = MAX;  
     else 
      {
       {
        if(mpp_thetaM[index_mpp]!=0) A_ltp = A0_ltp * (1.0/mpp_thetaM[index_mpp]);
	A_ltd = A0_ltd * (mpp_thetaM[index_mpp]);
	 
        if(t2_gc > t_mpp[index_mpp])
        ltp = A_ltp * exp(-(t2_gc - t_mpp[index_mpp]) / tau_ltp);
        if(t1_gc < t_mpp[index_mpp])
        ltd = A_ltd * exp((t1_gc - t_mpp[index_mpp]) / tau_ltd);

        //if(ltp>0) fprintf(f_mpp_ltp,"%d\t%6.3f\n",(t2_gc-t_mpp[index_mpp]),ltp);
        //if(ltd>0) fprintf(f_mpp_ltd,"%d\t%6.3f\n",(t1_gc-t_mpp[index_mpp]),-ltd);

        n->mpp_weight = n->mpp_weight * (1 + ltp-ltd);   
        if(n->mpp_weight < 0.03) n->mpp_weight = 0.03;
        if(n->mpp_weight > 0.08) n->mpp_weight = 0.08;   
       }
        index_mpp = index_mpp + 1;
      }
    }

  while(index_lpp < MAX)
    {if(t_lpp[index_lpp] == 0) index_lpp = MAX;
     else
      {     
       {
	if(lpp_thetaM[index_lpp]!=0) A_ltp = A0_ltp * (1.0/lpp_thetaM[index_lpp]);
	A_ltd = A0_ltd * (lpp_thetaM[index_lpp]);
     
        if(t2_gc > t_lpp[index_lpp])
        ltp = A_ltp * exp(-(t2_gc - t_lpp[index_lpp]) / tau_ltp);
        if(t1_gc < t_lpp[index_lpp])
        ltd =  A_ltd * exp((t1_gc - t_lpp[index_lpp]) / tau_ltd); 

        n->lpp_weight = n->lpp_weight * (1 + ltp-ltd);
        if(n->lpp_weight < 0.03) n->lpp_weight = 0.03;
        if(n->lpp_weight > 0.08) n->lpp_weight = 0.08;
        }
        index_lpp = index_lpp + 1;
       } 
      }

  for(k = 0; k < MAX; ++k) 
  {t_mpp[k] = 0; t_lpp[k] = 0; mpp_thetaM[k] = 0.0; lpp_thetaM[k] = 0.0;}
  index_mpp = index_lpp = 0; 
  t1_gc = t2_gc;

} /* end void stdp_train_cell(n) */


/********************************************************************** 
                   BCM training
**********************************************************************/
void bcm_threshold(int time, neuron n)
{

  n->mem = (n->mem * (n->pf - 1.0)  + (n->output * n->output)) / n->pf;
  n->pf = n->pf * n->boltzfactor + 1.0;
  n->theta_M = theta_M0 * n->mem;

} /* end void bcm_threshold(time, n) */

/************************************************************************
This function will test the cell every so often
************************************************************************/
void test_sample(int time, neuron n)
{  float percent_mpp_w, percent_lpp_w;
   
  //if(n->output > 0) fprintf(f_gc_voltage,"%f\t%1d\n",(float)time/minutes,1);

   if((time == 0) || (time - last_sample >= SAMPLE_PERIOD))
   {
    last_sample = time;  
    percent_mpp_w = 100*(n->mpp_weight-0.05)/0.05;
    percent_lpp_w = 100*(n->lpp_weight-0.05)/0.05;
    av_mpp_weight[k] =  av_mpp_weight[k] + percent_mpp_w;
    av_lpp_weight[k] =  av_lpp_weight[k] + percent_lpp_w;
    av_thetaM[k] = av_thetaM[k] + n->theta_M;
    k = k + 1;

   }

} /* end void test_sample(time, n) */

/***********************************************************************
This function returns the current time in miniseconds
***********************************************************************/
void getTime(long *z)
{
 struct timeval tp;
 struct timezone tzp;
  
 gettimeofday(&tp,&tzp);
 *z = tp.tv_usec;
}