//Author: Aaron Miller
//Last Modified: 3/02/2010

#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string>
#include <iostream>
#include <fstream>
#include <sstream>
using namespace std;
#include "synfireGrowth.h"
#include "ran1.h"

//random number seed
long seed = -100000;

double DT=.1; //timestep (ms)
double INV_DT=10; //(1/DT)

int runid =1;  //identifies the particular run
bool LOAD = false; //1 means load data
bool ILOAD = false;
int trials = 200000, trial_t = 2000, *train_lab; //default # of trials, length of trial (ms)

//Save strings
string loadpath, iloadpath, synapse_path="syn/syn", roster_path="roster/roster", volt_e_path="volt/volt_ensemble"; 
string volt_t_path="volt/volt_track_",conn_path="connect/connect", r_conn_path = "read/net_struct";
string dat_suff = ".dat", txt_suff = ".txt", sim_base=".", sim_lab;
int synapse_save = 100, roster_save = 1000, volt_save = 1000, conn_save = 1000; //save data intervals

//Tracking variables during trial
ofstream track_volt; //stream that tracks volt and conds during trial
int volt_post, volt_pre=0; //contains the label of the neuron whose voltage is being tracked, it postsynaptic to volt_pre

double t = 0.0; //current time in ms
int *whospike, whocount = 0;//tracks labels of neurons that spike during current step, length of whospike[]
int group_rank[SIZE]; //contains group ranking of each neuron after chain network forms

//Spontaneous activity defaults
double exfreq=40, infreq=200, examp=1.3, inamp=.1, global_i = .3, inh_d=0, leak = -85.0;

//Synapses defaults
int NSS = SIZE, tempNSS=10; //max # of supersynapses allowed
double act = .2, sup =.4, cap = .6, frac = .1, isynmax = .3, eq_syn = .3;
double syndec = .99999;
int conn_type=1;
bool plasticity=true;
double window=200; //history time window size (ms)

//Training defaults
int ntrain = 10, ntrg = 1; //number of training neurons
bool man_tt=false; //manual training time bool (false if training occurs at t=0)
double *train_times;
double training_f = 1.5; //training spike frequency (ms)^(-1)
double training_amp = .7; //training spike strength
double training_t = 8.0; //training duration in ms

//Stats
bool stats_on=true;
int sc=0; //total spike counter
double av = 0;
time_t rock, roll;

//-----------------------------------------------------------------------------------------------------

//Neuron object
class neuron {

 private:
  double nvolt, gexh, ginh, *spkhist;//membrane potential; excitatory, inhibitory conductances; spike times
  int latcount, refcount, spkcount, label;//latent, refractory period counters; spike counter; neuron label
  double spexfreq, spinfreq, LEAKREV;//spontaneous excitatory, inhibitory frequencies
  double spexamp, spinamp;//amplitudes of excitatory, inhibitory excitations 
  double global_inh; //amplitude of global inhibition
  
  //Integrate & Fire model parameters
  const static double INDECAY = 3.0, EXDECAY = 5.0, MEMDECAY=20.0;//inhibitory, excitatory, membrane time constants in ms 
  const static double SPKTHRES = -50.0, RESET = -80.0;//spike threshold, membrane reset potentials in mV
  const static double INHIBREV = -75.0;//inhibitory reversal potentials in mV
  const static double REFTIME = 25.0, LATTIME = 2.0;//refractory, latency time intervals

 public:
  neuron (int lab, double e, double i, double ea, double ia, double g);
  void set_sp_freq(double e, double i) {spinfreq = i*.001; spexfreq = e*.001;}; //convert from Hz to (ms)^(-1)
  void resetval(); 
  void resetspike() {if(spkcount!=0) {spkcount = 0; delete[] spkhist;}};
  int get_spkhist(double* & spikes) {if (spkcount!=0) spikes=spkhist; else spikes=NULL; return spkcount;};
  int get_label() {return label;}; //return label
  double get_pvt_var(char code); //return private variables
  bool updatevolt (); //4th order RK algorithm updates membrane potentials and conductances
  void neur_dyn(bool no_volt);
  void recspike(double t); //record spike times in spkhist arrays
  void excite_inhibit(double amp, char p);//increment conductances
};

//constuctor
neuron::neuron (int lab, double e, double i, double ea, double ia, double g){

  set_sp_freq(e, i);
  spexamp=ea;
  spinamp=ia;
  global_inh = g;
  label = lab;
  latcount = 0;
  refcount = 0;
  spkcount = 0;
  resetval();
  LEAKREV = leak;

}

void neuron::neur_dyn(bool no_volt){

  //update membrane potential,conductances with 4th order RK
  double c1[3], c2[3], c3[3], c4[3];
  c1[0]=-DT*gexh/EXDECAY;
  c1[1]=-DT*ginh/INDECAY;
  if (no_volt==false) c1[2] = (DT/MEMDECAY)*((LEAKREV-nvolt)-gexh*nvolt+ginh*(INHIBREV-nvolt));
  else c1[2]=0;
  
  c2[0] = -DT*(gexh+c1[0]/2.0)/EXDECAY;
  c2[1] = -DT*(ginh+c1[1]/2.0)/INDECAY;
  if (no_volt==false) c2[2] = (DT/MEMDECAY)*(LEAKREV-(nvolt+(c1[2]/2.0))-(gexh+c1[0]/2.0)*(nvolt+(c1[2]/2.0))+(ginh+c1[1]/2.0)*(INHIBREV-(nvolt+(c1[2]/2.0))));
  else c2[2]=0;
  
  c3[0] = -DT*(gexh+c2[0]/2.0)/EXDECAY;
  c3[1] = -DT*(ginh+c2[1]/2.0)/INDECAY;
  if (no_volt==false) c3[2] = (DT/MEMDECAY)*(LEAKREV-(nvolt+(c2[2]/2.0))-(gexh+c2[0]/2.0)*(nvolt+(c2[2]/2.0))+(ginh+c2[1]/2.0)*(INHIBREV-(nvolt+(c2[2]/2.0))));
  else c3[2]=0;
  
  c4[0] = -DT*(gexh+c3[0])/EXDECAY;
  c4[1] = -DT*(ginh+c3[1])/INDECAY;
  if (no_volt==false) c4[2] = (DT/MEMDECAY)*(LEAKREV-(nvolt+c3[2])-(gexh+c3[0])*(nvolt+(c3[2]))+(ginh+c3[1])*(INHIBREV-(nvolt+c3[2])));
  else c4[2]=0;
  
  nvolt += (c1[2]+2*c2[2]+2*c3[2]+c4[2])/6.0;
  gexh += (c1[0]+2*c2[0]+2*c3[0]+c4[0])/6.0;
  ginh += (c1[1]+2*c2[1]+2*c3[1]+c4[1])/6.0;
}

bool neuron::updatevolt(){

  bool isspike = false;

  //spontaneous excitation and inhibition
  if (ran1(&seed)<DT*spexfreq) excite_inhibit(spexamp*ran1(&seed), 'e');
  if (ran1(&seed)<DT*spinfreq) excite_inhibit(spinamp*ran1(&seed), 'i');
  
  if (latcount<1 && refcount<1) {//if neuron isn't in latent period before spike

    neur_dyn(false);
    
    //neuron goes into latency before spike if neuron potential is greater than theshold & is not in refractory state
    if (nvolt>=SPKTHRES) {
      nvolt = 0; //spike
      latcount=(int)(LATTIME*INV_DT); //sets latency counter 
    }
    
  }
  else{

    if(refcount>=1){
      refcount -= 1; //update refractory period counter
    } 

    if(latcount>=1){//if neuron is in latency period before spike
      latcount -= 1; //update latency counter
      if(nvolt==0){
	nvolt = RESET;
      }
      if (latcount<1) {//if latent time period ends on this timestep, neuron spikes
	latcount=0; //reset counter (should already be zero, but just in case)
	refcount=(int)(REFTIME*INV_DT); //set refractory time counter
	isspike = true;
      }
    }

    neur_dyn(true);
   
  }

  return isspike;
}

void neuron::excite_inhibit(double amp, char p){
  //p=='e' means excitatory, p=='i' means inhibitory
  switch (p) {
  case 'e': gexh+=amp;
    break;
  case 'i': ginh+=amp;
    break;
  default: cout<<"neuron.excite called with invalid char arg"<<endl; exit(1);
    break;
  }
}

double neuron::get_pvt_var(char code){
  //code = 'v' means membrane volt, code = 'e' means excitatory cond, code = 'i' means inhibitory cond.
  switch (code){
  case 'e': return gexh;
    break;
  case 'i': return ginh;
    break;
  case 'v': return nvolt;
    break;
  default: cout<<"neuron.get_pvt_var called with invalid char arg"<<endl; exit(1);
    break;
  }
}

void neuron::recspike(double t){

  if (spkcount!=0){
    double *temp = new double[spkcount]; //allocate temp array
    if (temp == NULL){
      cerr << "Memory allocation error in neuron.recspike" <<endl;
      exit(1);
    }
    //copy previous times to temp
    for (int i=0; i<spkcount; i++){
      temp[i]=spkhist[i];
    }

    delete[] spkhist; //deallocate old array
    spkhist = new double[spkcount + 1]; //allocate new (larger) array
    if (spkhist == NULL){
      cerr << "Memory allocation error in neuron.recspike" <<endl;
      exit(1);
    }
    //copy from temp to new array
    for (int i=0; i<spkcount; i++){
      spkhist[i] = temp[i];
    }
    delete [] temp; //deallocate temp array
  }
  else {
    spkhist = new double[1];
  }
  
  spkhist[spkcount] = t; //store new time
  spkcount++; //increment counter
}

void neuron::resetval(){ //randomly selects initial conditions before the beginning of each trial (see README)

  static double gexh_avg = .5*EXDECAY*spexamp*spexfreq;
  static double ginh_avg = .5*INDECAY*spinamp*spinfreq;

  nvolt = ran1(&seed)*(-55+80)-80;
  //select a conductance
  gexh =2*gexh_avg*ran1(&seed); 
  ginh =2*ginh_avg*ran1(&seed);
}

//----------------------------------------------------------------------------------------------------------------

//Synapses Object
class synapses {

 private:
  
  double G[SIZE][SIZE]; //synaptic strength matrix
  int actcount[SIZE], supcount[SIZE], NSS[SIZE]; //arrays tracking numbers of active and super cns of each neuron in the network
  int *actsyn[SIZE], *supsyn[SIZE]; //arrays containing the postsynaptic cns of each neuron in the network
  double actthres, supthres, synmax; //active, super thresholds, synapse cap
  double syndec, GLTP; //synaptic decay
  //Synaptic plasticity parameters
  const static double ALTP = .01, ALTD = .0105; //ltp parameter, ltp amplitude, ltd parameter
  const static double INV_DECTLTP = .05, INV_DECTLTD = .05; //inverses of ltp, ltd decay times in inv(ms)
  const static double POTFUNT = 5.0, DEPFUNT = 5.25;

 public:
  synapses (double fract_act, double glob, double a, double s, double m, double d, int form_opt); //initialize new network;
  synapses (ifstream& synfile, double a, double s, double m, double d); //build synapses object from file
  void activate (char p, int pre, int post); //activate synapse G[pre][post]
  void deactivate (char p, int pre, int post); //deactivate synapse G[pre][post]
  double getsynstr(int pre, int post) {return G[pre][post];}; //return synaptic strength  
  int getpost(char p, int pre, int* & post); //get postsynaptic neuron labels
  void synaptic_plasticity(int spiker, double t, neuron **net); //update synaptic strengths after network spikes
  void checkthres(double new_syn,int pre,int post,char p,char q); //check if synapse crosses threshold
  double potfun(double time, int spsc, double * hist, char p); //potentiation/depression function 
  int count_syns(char p); //count the active or super synapses
  void synaptic_decay(); //scale down synapses by syndec
  void ordersyn(char p, int k); //order synapse list sequentially by label (for analysis purposes)
  void write_syns(ofstream& outfile); //write synaptic matrix to outfile
  void write_groups(ofstream& readfile, ofstream& datfile, int *group_rank, char p); //write connectivity to outfile
  void rank_groups(int *group_one, int size_group_one, int *group_rank, char p); //determine group structure of network
  void afferent(int * a, char p); //create afferent connection matrix
  int write_aff(ofstream& outfile, int * a, int post); //print afferent neuron labels of 'post' to outfile, given afferent connection matrix a
  void shortestpath(int * group, int group_size, int group_rank, int * res, char p); //compute shortest path to training neurons
  int findmode(int * vec, int ma); //determine group rank
  double getNSS(int lab) {return NSS[lab];};
};

synapses::synapses (double fract_act, double glob, double a, double s, double m, double d, int form_opt){
  
  actthres = a;
  supthres = s;
  synmax = m;
  syndec = d;
  GLTP = eq_syn;
  for(int i=0; i<SIZE; i++){
    NSS[i]=tempNSS;
  }

  switch(form_opt){

  case 1://randomly activate fract_act of all synapses

    for (int i=0; i<SIZE; i++){
      //cout<<i<<endl;
      actcount[i]=0;
      supcount[i]=0;
      for (int j=0; j<SIZE; j++){
	//cout<<"we are"<<endl;
        if (i!=j){
	  if (ran1(&seed)<=fract_act){
	    G[i][j]=actthres+(isynmax-actthres)*ran1(&seed);
	    activate('a',i,j);
	  }
	
	  else{
	    G[i][j]=actthres*ran1(&seed);
	  }
	  if(G[i][j]>synmax) G[i][j]=synmax;
	}
	else G[i][j]=0; //self synapses not allowed
      }
    
    }
    break;

  case 2:

    for (int i=0; i<SIZE; i++){
      actcount[i]=0;
      supcount[i]=0;
      for (int j=0; j<SIZE; j++){
	G[i][j]=glob;
      }
    }
    break;

  default:
    cout<<"Invalid form opt"<<endl;
    cout<<"add code to set an inital connectivity"<<endl;
    exit(0);
    break;
  }
  
  //int dum=0;
  //for (int i=0; i<SIZE; i++){
  //  dum += actcount[i];
  //}
  //cout<<"init act count = "<<dum<<endl;
  
}

synapses::synapses(ifstream& synfile, double a, double s, double m, double d) {
  
  actthres = a;
  supthres = s;
  synmax = m;
  syndec = d;
  for (int i=0; i<SIZE; i++){
    actcount[i]=0;
    supcount[i]=0;
    for (int j=0; j<SIZE; j++){
      synfile >> G[i][j];
      if (!synfile) {
	cerr << "load failed " <<i<<" "<<j<< endl;
	exit(1);
      }
      if (i==j) G[i][j]=0; //do not allow self connections
      if (G[i][j]>=actthres){
	activate('a',i,j);
	if (G[i][j]>=supthres){
	  activate('s',i,j);
	  if (G[i][j]>=synmax){
	    cerr<<"Saved matrix contains elements larger than synmax."<<endl;
	    exit(1);
	  }
	}
      }
    }
    if(supcount[i]>tempNSS){
      NSS[i]=supcount[i];
    }
    else{NSS[i]=tempNSS;}
  }
}

void synapses::activate (char p, int pre, int post) {
  //p == 'a' for active, p == 's' for super
  int *temp, *l, *m;
  switch (p){
  case 'a': 
    l = &actcount[pre]; //l points to # of act synapses of pre
    if ((*l)!=0){//if old actsyn exists
      temp = new int[(*l)]; //make temp array
      if (temp == NULL){
	cerr << "Memory allocation error in neuron.activate" <<endl;
	exit(1);
      }
      for (int i=0;i<(*l);i++){
	temp[i]=actsyn[pre][i]; //store current values of actsyn[pre] in temp
      }
      delete[] actsyn[pre]; //deallocate actsyn[pre]
    }
    actsyn[pre] = new int[(*l)+1]; //allocate larger actsyn[pre]
    if (actsyn[pre] == NULL){
      cerr << "Memory allocation error in neuron.activate" <<endl;
      exit(1);
    }
    m = actsyn[pre]; //m points to new array
    break;
    
  case 's':
    l = &supcount[pre]; //l points to # of super synapses of pre
    if ((*l)!=0){//if old actsyn exists
      temp = new int[(*l)]; //make temp array
      if (temp == NULL){
	cerr << "Memory allocation error in neuron.activate" <<endl;
	exit(1);
      }
      for (int i=0;i<(*l);i++){
	temp[i]=supsyn[pre][i]; //store current values of supsyn[pre] in temp
      }
      delete[] supsyn[pre]; //deallocate supsyn
    }
    supsyn[pre] = new int[(*l)+1]; //allocate larger supsyn[pre]
    if (supsyn[pre] == NULL){
      cerr << "Memory allocation error in neuron.activate" <<endl;
      exit(1);
    }
    m = supsyn[pre]; //m points to new array
    break;

  default: cout<<"synapses.activate called with invalid char arg"<<endl;
    break;
  }

  if((*l)!=0){//if temp exists
    //copy old vals from temp into new array
    for (int i=0; i<(*l); i++){
      m[i]=temp[i];
    }
    delete [] temp;//deallocate temp
  }
  
  m[(*l)]=post;//place new synapse in new array
  (*l)++;//increment counter
}

void synapses::deactivate(char p, int pre, int post) { 
  //p == 'a' active thres, p == 's'  for sup thres
  int *temp, *l, *m;
  switch (p){
  case 'a': 
    l = &actcount[pre]; //l points to # of act synapses of pre
    if ((*l)!=1){//if the last post is not being deactivated
      temp = new int[(*l)]; //make temp array
      if (temp == NULL){
	cerr << "Memory allocation error in neuron.deactivate" <<endl;
	exit(1);
      }
      for (int i=0;i<(*l);i++){
	temp[i]=actsyn[pre][i]; //store current values of actsyn[pre] in temp
      }
    }
    delete[] actsyn[pre]; //deallocate actsyn[pre]
    if ((*l)!=1){//if the last post is not being deactivated
      actsyn[pre] = new int[(*l)-1]; //allocate smaller actsyn[pre]
      if (actsyn[pre] == NULL){
	cerr << "Memory allocation error in neuron.deactivate" <<endl;
	exit(1);
      }
      m = actsyn[pre]; //m points to new array
    }
    break;
    
  case 's':
    l = &supcount[pre]; //l points to # of super synapses of pre
    if ((*l)!=1){//if the last post is not being deactivated
      temp = new int[(*l)]; //make temp array
      if (temp == NULL){
	cerr << "Memory allocation error in neuron.deactivate" <<endl;
	exit(1);
      }
      for (int i=0;i<(*l);i++){
	temp[i]=supsyn[pre][i]; //store current values of supsyn[pre] in temp
      }
    }
    delete[] supsyn[pre]; //deallocate supsyn
    if ((*l)!=1){//if the last post is not being deactivated
      supsyn[pre] = new int[(*l)-1]; //allocate smaller supsyn[pre]
      if (supsyn[pre] == NULL){
	cerr << "Memory allocation error in neuron.deactivate" <<endl;
	exit(1);
      }
      m = supsyn[pre]; //m points to new array
    }
    break;

  default: cout<<"synapses.deactivate called with invalid char arg"<<endl;
    break;
  }
  
  if ((*l)!=1){//if the last post is not being deactivated
    //don't include post when placing values in temp back into supsyn[pre], so move the last value into its place
    for(int i=0; i<(*l); i++){
      if(temp[i]==post){
	temp[i]=temp[(*l)-1];
	break;
      }
    }
    (*l)--; //decrease counter
    
    //copy old vals from temp into new array
    for (int i=0; i<(*l); i++){
      m[i]=temp[i];
    }
    
    delete [] temp;//deallocate temp
  }
  else{//if the last post is being deactivated
    (*l)--;
  }
}

int synapses::getpost(char p, int pre, int* & post){
  
  int res = -1;
  switch(p){
  case 'a':
    post = actsyn[pre];
    res=actcount[pre];
    break;
  case 's':
    post = supsyn[pre];
    res=supcount[pre];
    break;
  default: cout<<"synapses.getpost called with invalid char arg"<<endl;
    exit(1);
  }
  return res;
}

void synapses::synaptic_plasticity(int spiker, double t, neuron **net){

  double * spk_times;
  int spk_count;
  double tempPot, tempDep, GPot, GDep;
  for (int k = 0; k<SIZE; k++){
    if(k!=spiker){
      
      GPot = G[k][spiker];
      GDep = G[spiker][k];
      spk_count = net[k]->get_spkhist(spk_times);
      if(spk_count!=0 && spk_times[spk_count-1]+window>=t){
	tempPot = GPot+ALTP*GLTP*potfun(t, spk_count, spk_times, 'p'); //potentiation
	tempDep = GDep*(1-ALTD*potfun(t, spk_count, spk_times, 'd')); //depression
	
	if(tempPot>synmax) tempPot=synmax;
	if(tempDep<0) tempDep=0;
	
	//Potentiate G[k][spiker]
	if(supcount[k]<NSS[k] || (supcount[k]==NSS[k] && GPot>=supthres)){
	  checkthres(tempPot,k,spiker,'a','p');
	  checkthres(tempPot,k,spiker,'s','p');
	  G[k][spiker] = tempPot;
	}
	//Depress G[spiker][k]
	if(supcount[spiker]<NSS[spiker] || (supcount[spiker]==NSS[spiker] && GDep>=supthres)){
	  checkthres(tempDep,spiker,k,'a','d');
	  checkthres(tempDep,spiker,k,'s','d');
	  G[spiker][k] = tempDep;
	}
      }
    }
  }
}

double synapses::potfun(double time, int spsc, double * hist, char p) {
  //p=='p' for potentiation, p=='d' for depression
   double res=0.0;
   double delt,pwt,inv_dect;
   double a=0.0;
   
   switch(p){
   case 'p': 
     pwt=POTFUNT;
     inv_dect=INV_DECTLTP;
     break;
   case 'd':
     pwt = DEPFUNT;
     inv_dect=INV_DECTLTD;
     break;
   default: cout<<"Invalid input to synapses.potfun"<<endl;
     break;
   }
   
   for(int i=spsc-1; i>=0; i--){
     delt = time - hist[i];
     //if(p=='d'){
     //cout<<delt<<endl;
     //}
     if(delt<=window){
       if (delt <= pwt){
	 a = delt/pwt;
       }
       else {
	 a = exp(-(delt - pwt)*inv_dect);
       }
       res += a;
     }
     else break;
   }
  return res;
}  

void synapses::checkthres(double new_syn, int pre, int post, char p, char q){
  //p=='a'/p=='s', check if G crossed active/super threshold 
  //q == 'p' for potentiation, q == 'd' for depression
  double thres;
  switch(p){
  case 'a': thres = actthres;
    break;
  case 's': thres = supthres;
    break;
  default: cout<<"Invalid input for p in synapses.checkthres"<<endl;
    break;
  }

  switch(q){
  case 'p':
    if(G[pre][post]<thres && new_syn>=thres) activate(p, pre, post);
    break;
  case 'd':
    if(G[pre][post]>=thres && new_syn<thres) deactivate(p, pre, post);
    break;
  default: cout<<"Invalid input for q in synapses.checkthres"<<endl;
    break;
  }
}

void synapses::synaptic_decay(){

  for (int i=0; i<SIZE; i++){
    for (int j=0; j<SIZE; j++){
      checkthres(G[i][j]*syndec,i,j,'a','d');
      checkthres(G[i][j]*syndec,i,j,'s','d');
      G[i][j]*=syndec;
    }
  }
}

int synapses::count_syns(char p){
  
  int sum=0;
  switch(p){
  case 'a':
    for (int i=0; i<SIZE; i++){
      sum += actcount[i];
    }
    break;
  case 's':
    for (int i=0; i<SIZE; i++){
      sum += supcount[i];
    }
    break;
  default:cout<<"Invalid input to synapses.count_syns"<<endl;
    break;
  }

  return sum;
}

void synapses::write_syns(ofstream& outfile){

  for (int i=0; i<SIZE; i++){
    for (int j=0; j<SIZE; j++){
      outfile<<G[i][j]<<" ";
    }
    outfile<<endl;
  }
}

void synapses::write_groups(ofstream& readfile, ofstream& datfile, int* group_rank, char p){

  int c, *post, *throwaway, max_rank=0;
  for (int k=0; k<SIZE; k++){
    if (group_rank[k]>max_rank){
      max_rank = group_rank[k];//compute max_rank
    }
  }

  //Readable file
  switch (p){
  case 'a':
    readfile<<"##Active synapse network##"<<endl;
    break;
  case 's':
    readfile<<"##Super synapse network##"<<endl;
    break;
  default:cout<<"synapses.write_groups called with invalid char"<<endl;
    break;
  }
  for (int i=1; i<=max_rank; i++){
    readfile<<"GROUP "<<i<<":"<<endl;
    for (int j=0; j<SIZE; j++){
      if (group_rank[j]==i){
	readfile<<j<<": ";
	ordersyn('s',j);
	c = getpost(p,j,post);
	for (int k=0; k<c; k++){
	  readfile<<post[k]<<"(G"<<group_rank[post[k]]<<") ";
	}
	readfile<<endl;
      }
    }
    readfile<<endl;
  }
  
  //Data file
  for (int i=0; i<SIZE; i++){
    if (group_rank[i]!=0){
      c = getpost(p,i,post);
      for (int j=0; j<c; j++){
	datfile<<i<<" "<<group_rank[i]<<" "<<post[j]<<" "<<group_rank[post[j]]<<" "<<(getpost(p,post[j],throwaway)==tempNSS)<<" "<<G[i][post[j]]<<endl;
      }
    }
  }
}

  void synapses::rank_groups(int *group_one, int size_group_one, int *group_rank, char p){

  int aff[SIZE][SIZE], short_rank[SIZE];
  //Initialize these arrays
  for (int i=0; i<SIZE; i++){
    short_rank[i]=-1;
    group_rank[i]=-1;
    for (int j=0; j<SIZE; j++){
      aff[i][j]=0;
    }
  }
  
  afferent(&aff[0][0], p); //aff[c][d] contains 1 if G[d][c] is a supersynapse
  //short_rank[i] contains the length of the shortest path via synapses of type 'p' from neuron i to a neuron in group_one
  shortestpath(group_one,size_group_one,0,short_rank,p); 
  
  for(int i=0; i<SIZE; i++){
    for (int j=0; j<SIZE; j++){
      aff[i][j]*=(short_rank[j]+1); //rank the afferent connections
    }
    //compute group_rank by finding mode of the short_rank of all afferent neurons
    group_rank[i] = findmode(aff[i],SIZE/size_group_one)+1; 
  }

  for (int i=0; i<size_group_one; i++){
    group_rank[group_one[i]]=1; //set group_one's member's group_rank to one
  }
}
	  
void synapses::afferent(int * a, char p){
  //p=='a' for active connections, p=='s' for super connections

  //aff[c][d] contains 1 if G[d][c] is a supersynapse
  switch (p){
  case 'a':
    for (int i=0; i<SIZE; i++){
      for (int j=0; j<actcount[i]; j++){
	*(a+SIZE*actsyn[i][j]+i)=1;
      }
    }
    break;
  case 's':
    for (int i=0; i<SIZE; i++){
      for (int j=0; j<supcount[i]; j++){
	*(a+SIZE*supsyn[i][j]+i)=1;
      }
    }
    break;
  default: cout<<"Invalid input to synapses.afferent."<<endl;
    break;
  }
}

int synapses::write_aff(ofstream& outfile, int * a, int post){
  int y=0; //tracks total # of afferent neurons
  for (int i=0; i<SIZE; i++){
    if (*(a+SIZE*post+i)!=0){
      y++;
      outfile<<i<<" ";
    }
  }
  outfile<<endl;
  return y;
}

void synapses::ordersyn(char p, int k){
  //p is 'a' for active, 's' for sup
  int count,temp;
  int * me;

  temp=0;  
  
  switch(p){
  case 'a':
    count = actcount[k];
    me = actsyn[k];
    break;
  case 's':
    count = supcount[k];
    me = supsyn[k];
    break;
  default:cout<<"Invalid input to synapses.ordersyn"<<endl;
    break;
  }
  
  for (int i=0; i<(count-1); i++){
    for (int j=0; j<(count-1); j++){
      if(me[j]>me[j+1]){
	temp = me[j];
	me[j]=me[j+1];
	me[j+1]=temp;
      }
    }
  } 
}

void synapses::shortestpath(int * group, int group_size, int group_rank, int * res, char p){
  //group is a set of <group_size> neurons which are a shortest distance of <group_rank> from group_one
  //p=='a' means follow all active connections, p=='s' means follow only superconnections
  for (int i=0; i<group_size; i++){
    if ((res[group[i]]>group_rank || res[group[i]]==-1)){//shortest_path passed through "res" must be initialized to -1
      res[group[i]]=group_rank;//change current value of shortest_path if current path is shorter
      
      switch(p){
      case 's':
	if (supcount[group[i]]!=0){
	  shortestpath(supsyn[group[i]],supcount[group[i]],group_rank+1,res,p);
	}
	break;
	
      case 'a':
	if (actcount[group[i]]!=0){
	  shortestpath(actsyn[group[i]],actcount[group[i]],group_rank+1,res,p);
	}
	break;

      default: cout<<"Invalid char input to synapses.shortestpath."<<endl;
	break;
      }
    }
  }
}

int synapses::findmode(int * vec, int ma){
//determines group rank of a neuron whose afferent connections are contains in vec (length SIZE)
  int res[ma];
  for (int i=0; i<ma; i++){
    res[i]=0;
  }
  int mode = -1;
  int count = -1;
  for (int i=0; i<SIZE; i++){
    if(vec[i]!=0){
      res[vec[i]-1]++;
      if (res[vec[i]-1]>count){
	count = res[vec[i]-1];
	mode = vec[i];
      }
    }    
  }

  //check for tie, goes to lowest group number
  int check=0;
  for (int i=0; i<ma; i++){
    check += (int)(res[i]==count);
  }

  if (check>=2){
    for (int i=0; i<ma; i++){
      if (res[i]==count){
	mode = i+1;
	break;
      }
    }
  }

  return mode;
}

//-------------------------------------------------------------------------------------------------

int main(int argc, char *argv[]) {

  //set arguments passed in command line to override defaults
  for (int i=1; i<argc; i++){
    
    string commnd = argv[i];
    if (commnd[0] != '-') continue;
    else{
      switch (commnd[1]){
	
	//"DEFAULT" SIMUATIONS/USER-DEFINED FLAGS (in addition to the absolute default:1000 neuron network) 
      case '2': //'-2' sets the default parameters for a typical 200 neuron network
	exfreq = 30, examp=1.88, syndec = .99992;
	break;
      case '4':
	exfreq = 33, examp=1.57, syndec = .99994;
	break;
	//OUTPUT & ANALYSIS FLAGS
      case 'd': // '-d' outputs various statistics after each trial for diagnostic purposes (i.e. spk count, avg volt, etc)
	stats_on = 1;
	break;
      case 't': // '-t <int>' saves synapse data every <int> trials (default is 1000)
	synapse_save = atoi(argv[i+1]);
	roster_save = synapse_save;
	break;
      case 'v':// '-v <int1> <int2>' tracks membrane voltage of neuron <int2> and one of its postsynaptic neurons 
	       // every <int1> trials (default is int1=1000, int2=0 (training neuron))
	volt_save = atoi(argv[i+1]);
	if (argv[i+2][0]!='-') volt_pre = atoi(argv[i+2]);
	break;
      case 'w':// '-w <int>' saves connectivity info in readable format and also in a format suitable for analysis (default is 1000)
	conn_save = atoi(argv[i+1]);
	break;
	//PARAMETER FLAGS
      case '#': // '-# <int>' sets the ID number for the run (default is 1)
	runid = atoi(argv[i+1]);
	break;
      case 'A':// '-A <int>' sets the numbers of trials before termination of the run (default is 200000)
	trials = atoi(argv[i+1]);
	break;
      case 'B':// '-B <int>' sets the number of training neurons (default is 10)
	ntrain = atoi(argv[i+1]);
	break;
      case 'C':// '-C <int>' sets the number of ms per trial (default is 2000);
	trial_t = atoi(argv[i+1]);
	break;
      case 'D':// '-D <int>' sets the max number of supersynapses a neuron may have
	tempNSS = atoi(argv[i+1]);
	break;
      case 'F':
	ntrg = atoi(argv[i+1]);
	train_times = new double[ntrg+1];
	man_tt=true;
	for(int j=0; j<ntrg; j++){
	  train_times[j]=atof(argv[i+j+2]);
	}
	train_times[ntrg]=trial_t+1000;
	break;
      case 'J'://load inhibitory weights from file
	ILOAD=true;
	iloadpath = argv[i+1];
	break;
      case 'L':// synapses are static
	plasticity=false;
	break;
      case 'R'://set leak reversal (default is -85)
	leak = atof(argv[i+1]);
	leak = leak*(-1);
	break;
      case 'X': //base directory
	sim_base = argv[i+1];
	break;
      case 'Y': //numeric directory label
	sim_lab = argv[i+1];
	break;
      case 'a':// sets connectivity type (see README)
	conn_type = atoi(argv[i+1]);
	break;
      case 'c':// '-c <double>' sets the decay rate of the synapses (default is .999996)
	syndec = atof(argv[i+1]);
	break;
      case 'f':// '-f <double>' sets the fraction of synapses initially active (default is .1)
	frac = atof(argv[i+1]);
	if (frac>=1 || frac<0){
	  cerr << "Command line input -f is invalid, must be <1 and >=0" << endl;
	  exit(1);
	}
	break;
      case 'g': //sets GLTP
	eq_syn = atof(argv[i+1]);
	break;
      case 'i':// '-i <double>' sets the amplitude of the global inhibition (default is .3)
	global_i = atof(argv[i+1]);
	break;
      case 'j':// sets the inhibition delay
	inh_d = atof(argv[i+1]);
	break;
      case 'l': // '-l <string>' loads synapse data from path <string>
	LOAD = 1;
	loadpath = argv[i+1];
	break;
      case 'm': // '-m <double>' sets the excitatory spontaneous frequency (default is 40Hz)
	exfreq = atof(argv[i+1]);
	break;
      case 'n':// '-n <double>' sets the inhibitory spontaneous frequency (default is 200Hz)
	infreq = atof(argv[i+1]);
	break;
      case 'o':// '-o <double>' sets the amplitude of the spontaneous excitation (default is 1.3)
	examp = atof(argv[i+1]);
	break;
      case 'p':// '-p <double>' sets the amplitude of the spontaneous inhibition (default is .1)
	inamp = atof(argv[i+1]);
	break;
      case 'q':// '-q <double>' sets the active syn threshold (default is .2)
	act = atof(argv[i+1]);
	break;
      case 'r':// '-r <double>' sets the super syn threshold (default is .4)
	sup = atof(argv[i+1]);
	break;
      case 's':// '-s <double>' sets the synapse maximum (default is .6)
	cap = atof(argv[i+1]);
	break;
      case 'u'://set maximum inhibitory synaptic strength
	isynmax = atof(argv[i+1]);
	break;
      case 'x': // '-x <double>' sets spike rate of training excitation (default is 1500 Hz)
	training_f = atof(argv[i+1])*.001;//convert to (ms)^(-1)
	break;
      case 'y': // '-y <double>' sets amplitude of training excitation (default is .7)
	training_amp = atof(argv[i+1]);
	break;
      case 'z': // '-z <double>' sets the training excitation duration (default is 8 ms)
	training_t = atof(argv[i+1]);
	break;
      default:
	cout<<"Warning: command line flag "<< commnd <<" is not recognized."<<endl;
	break;
		      
      }
    }
  }
  
  //Create neurons and synapses
  synapses *connectivity;
  synapses *inhibition_strength;
  neuron *network[SIZE];
  int group_s = (int) SIZE/ntrg;
  
  //Seed random number generator & check RAND_MAX
  seed = time(NULL)*(-1);
  //ran1(&seed);

  //initialize synapse object
  if (LOAD == false){
    //cout<<"loading"<<endl;
    connectivity = new synapses(frac,0,act,sup,cap,syndec,conn_type);
    if (connectivity == NULL){
      cerr << "Memory allocation error while creating synapses object" <<endl;
      exit(1);
    }
  }
  
  if (LOAD == true){
    ifstream infile;
    infile.open(loadpath.c_str());
    if (!infile){
      cerr<<"Failed to open "<<loadpath<<" containing synapse matrix"<<endl;
      exit(1);
    }
    connectivity = new synapses(infile,act,sup,cap,syndec);
    if (connectivity == NULL){
      cerr << "Memory allocation error while creating synapses object" <<endl;
      exit(1);
    }
  }

  //inhibition
  if(ILOAD==false){
    inhibition_strength = new synapses(1,global_i,1,1,1,1,2);
    if (inhibition_strength == NULL){
      cerr << "Memory allocation error while creating inhibition object" <<endl;
      exit(1);
    }
  }

  if(ILOAD==true){
    ifstream infile;
    infile.open(iloadpath.c_str());
    if (!infile){
      cerr<<"Failed to open "<<iloadpath<<" containing synapse matrix"<<endl;
      exit(1);
    }
    inhibition_strength = new synapses(infile,1000,1000,1000,1);
    if (inhibition_strength == NULL){
      cerr << "Memory allocation error while creating synapses object" <<endl;
      exit(1);
    }
  }
  
  //initialize neuron objects
  for (int i=0; i<SIZE; i++){
    network[i]= new neuron(i,exfreq,infreq,examp,inamp,global_i);
    if (network[i] == NULL){
      cerr << "Memory allocation error while creating neuron object" <<endl;
      exit(1);
    }
  }
  
  if(man_tt==false){
    train_times = new double[2];
    train_times[0]=0.0;
    train_times[1]=trial_t+1000;
  }

  //initialize array of training neuron labels
  train_lab = new int[ntrain*ntrg];
  if (train_lab == NULL){
    cerr<<"Could not allocate memory for train_lab in main."<<endl;
    exit(1);
  }
  int tc=0;
  for(int j=0; j<ntrg; j++){
    for (int i=0; i<ntrain; i++){
      train_lab[tc]=j*group_s+i;
      tc++;
    }
  }
  
  //set up inhibition delay
  int dsteps = 1 + (int) INV_DT*inh_d;
  int inh[dsteps];
  for(int i=0; i<dsteps; i++){
    inh[i] = 0;
  }

 //Write simulation parameters in logfile before simulation starts
  stringstream id_num;
  id_num << runid;
  string fname = sim_base+sim_lab+"/log"+id_num.str()+".txt";
  
  ofstream writelog;
  writelog.open(fname.c_str());

  writelog<<"Run started with command line:"<<endl;
  for (int i=0; i<argc; i++){
    writelog<<argv[i]<<" ";
  }
  writelog<<endl;
  writelog<<"SIZE = "<<SIZE<<" :: Max # of supersynapses = "<<tempNSS<<endl;
  writelog<<"Trials = "<<trials<<" :: Trial length  = "<<trial_t<<"ms"<<endl;
  writelog<<"Timestep = "<<DT<<"ms"<<endl;
  writelog<<endl<<"Neuron class:"<<endl;
  writelog<<"Sp. excitatory freq. = "<<exfreq<<"Hz :: Sp. inhibitory freq = "<<infreq<<"Hz"<<endl;
  writelog<<"Sp. excitatory amp. = "<<examp<<" :: Sp. inhibitory amp = "<<inamp<<endl;
  writelog<<"Ext current = "<<(leak+85)<<endl;
  if(ILOAD=true){
    writelog<<"Inh loaded from "<<iloadpath;
  }
  else{
    writelog<<"Global inh. amp. = "<<global_i;
  }
  writelog<<":: Inh delay = "<<inh_d<<endl;

  writelog<<endl<<"Synapses class:"<<endl;
  writelog<<"Act. thres. = "<<act<<" :: Sup. thres. = "<<sup<<" :: Max. syn. str. = "<<cap<<endl;
  writelog<<"Synapse decay factor = "<<syndec<<" :: Frac. initially active = "<<frac<<endl;
  writelog<<"GLTP = "<<eq_syn<<endl;
  writelog<<endl<<"Training:"<<endl;
  writelog<<"Training neurons = ";
  for(int i=0; i<ntrain; i++){
    writelog<<train_lab[i]<<" ";
  }
  writelog<<endl<<"Training time =  ";
  writelog<<train_times<<" ";
  writelog<<endl<<"Training duration = "<<training_t<<"ms"<<" :: Training freq = "<<training_f*1000<<"Hz :: Training amp = "<<training_amp<<endl;
  writelog<<endl<<"Data output and analysis:"<<endl;
  if (LOAD ==1){
    writelog<<"Synapse data loaded from "<<loadpath<<endl;
  }
  else writelog<<"Synapses initialized with option "<<conn_type<<endl;

  writelog<<"Synapse data saved every "<<synapse_save<<" trials"<<endl;
  writelog<<"Spike data saved every "<<roster_save<<" trials"<<endl;
  writelog<<"Neuron microvariables tracked and saved every "<<volt_save<<" trials"<<endl;
  writelog<<"Connection data saved every "<<conn_save<<" trials"<<endl;

  writelog.close();

  int trial_steps = (int) trial_t*INV_DT;
  
  for (int a=0; a<=trials; a++){//****Trial Loop****//
    rock = time(NULL);//start timer
    int train_time_counter=0;
    int train_group_lab = 0;
    
    //Save membrane potential data for a pre TN and a post 
    //(FORMAT: <pre.nvolt> <pre.gexh> <pre.ginh> <post.nvolt> <post.gexh> <post.ginh> <t>)
    if ((a%volt_save)==0){
      stringstream id_num, trial_num, id_pre, id_post;
      id_num << runid;
      trial_num << a;

      //select neuron to track (which isn't receiving training spikes)
      int *post;
      int c = connectivity->getpost('a',volt_pre,post);
      for (int i=0; i<c; i++){
	if (post[i]>=ntrain){
	  volt_post = post[i];
	  break;
	}
      }
      id_pre << volt_pre;
      id_post << volt_post;
      
      string fname = sim_base+sim_lab+"/"+volt_t_path + id_pre.str() + "-" + id_post.str() + "_" + trial_num.str() + "r" + id_num.str() + dat_suff;
      track_volt.open(fname.c_str()); 
    }
    
    //At t=0, on with the trial!
    for (int i=0; i<trial_steps; i++) {//****Timestep Loop****//
      //Excite training neurons
	
      if (t>=train_times[train_time_counter]){//****Training Loop****//
	for(int j=train_group_lab*ntrain; j<ntrain*(train_group_lab+1); j++){
	  if(ran1(&seed)<training_f*DT){
	    network[train_lab[j]]->excite_inhibit(training_amp, 'e');
	    //cout<<"excited "<<train_lab[j]<<" at "<<t<<endl;
	  }
	}
   
	if(t>=(train_times[train_time_counter]+training_t)){
	  //cout<<train_group_lab<<" "<<t<<endl;
	  train_time_counter++;
	  train_group_lab = 1;//rand()%ntrg;
	}
      }
      
      //Track voltages and conductances
      if ((a%volt_save)==0){
	track_volt<<network[volt_pre]->get_pvt_var('v')<<" "<<network[volt_pre]->get_pvt_var('e')<<" ";
	track_volt<<network[volt_pre]->get_pvt_var('i')<<" "<<network[volt_post]->get_pvt_var('v')<<" ";
	track_volt<<network[volt_post]->get_pvt_var('e')<<" "<<network[volt_post]->get_pvt_var('i')<<" "<<t<<endl;
      }
      
      //Update membrane potentials first, keep track of who spikes
      for (int j=0; j<SIZE; j++) {//****Membrane Potential Loop****
	if (network[j]->updatevolt()){//if neuron j spikes this timestep, increase the length of whospike, store label 
	  if (whocount!=0){
	    int *temp = new int[whocount];  //temp array while whospike is destroyed
	    if (temp==NULL){
	      cerr<<"Too many spikes this timestep, could not allocate memory for whospike"<<endl;
	      exit(1);
	    }
	    
	    for (int k=0; k<whocount; k++){
	      temp[k]=whospike[k]; //put spikes in temp
	    }
	    delete[] whospike;
	    whocount++;
	    whospike = new int[whocount]; //create longer whospike
	    if (whospike==NULL){
	      cerr<<"Too many spikes this timestep, could not allocate memory for whospike"<<endl;
	      exit(1);
	    }
	    
	    for (int k=0; k<(whocount-1); k++){
	      whospike[k]=temp[k]; //put spikes back into whospike
	    }
	    delete[] temp; //destroy temp
	  }
	  else{ //if this is the first spike this timestep
	    whocount++;
	    whospike = new int[whocount];
	    if (whospike==NULL){
	      cerr<<"Too many spikes this timestep, could not allocate memory for whospike"<<endl;
	      exit(1);
	    }
	  }
	  whospike[whocount-1]=j; //store spiking neuron in longer array
	}
      }

      for (int j=0; j<whocount; j++){//****Spike Loop****
	int *send_to;//pointer to array containing post neurons
	int send_count;//number of post neurons receiving spike

	//cout<<whospike[j]<<" spikes!!! at "<<t<<endl;
	network[whospike[j]]->recspike(t); //record time of spike
	sc++; //keep track of total number of spikes this trial

	//Send spikes
	if (connectivity->getpost('s',whospike[j],send_to)==connectivity->getNSS(whospike[j])){ //check to see if spiking neuron is saturated
	  for(int k = 0; k<connectivity->getNSS(whospike[j]); k++){
	    network[send_to[k]]->excite_inhibit(connectivity->getsynstr(whospike[j],send_to[k]),'e'); //send spikes along super synapses
	  }
	}
	else {//spiking neuron isn't saturated, send spikes along active connections
	  send_count = connectivity->getpost('a',whospike[j],send_to);
	  for(int k = 0; k<send_count; k++){
	    network[send_to[k]]->excite_inhibit(connectivity->getsynstr(whospike[j],send_to[k]),'e');
	  }
	}

	//Plasticity and remodeling
	if(plasticity==true){
	  connectivity->synaptic_plasticity(whospike[j],t,network);
	}
	
      }//****Spike Loop****
      
	//inhibition
    inh[dsteps-1]=whocount;

    for (int z=0; z<whocount; z++){

      for (int j=0; j<SIZE; j++){
	network[j]->excite_inhibit(inhibition_strength->getsynstr(whospike[z],j),'i');
      }
      //cout<<t<<" "<<inh[0]*global_i<<endl;
      
    }
    
    for(int i=0; i<dsteps-1; i++){
      inh[i]=inh[i+1];
    }
    
    if (whocount!=0) {
      delete [] whospike; //reset whocount and destroy whospike at end of timestep
      whocount = 0;
    }
    t += DT; //increment t
    
    }//****Timestep Loop****
    
    t=0.0; //reset timer
    seed = time(NULL)*(-1); //reseed
    ran1(&seed);
    
    //stop tracking volt when trial is done
    if ((a%volt_save)==0){
      track_volt.close();
    }
    
    if (plasticity==true){
      connectivity->synaptic_decay(); //synapses decay after each trial
    }

    if (stats_on==1){//(FORMAT: <trial> <spike total> <av. volt> <runtime> <# of active connections>)
      for (int i=0; i<SIZE; i++){
	av += network[i]->get_pvt_var('v');
      }
      roll = time(NULL);
      cout << a <<" "<< sc <<" "<< (av/SIZE) <<" "<<(roll-rock)<<" "<<connectivity->count_syns('a')<<endl;

      sc = 0;
      av = 0;
    }
    
    //SAVE DATA FILES AND WEIGHTS
    //Connectivity  files (FORMAT: <pre> <pre_group#> <post> <post_group#> <sat(1)/unsat(0)> <G[pre][post]>)
    if((a%conn_save)==0){
      /*
      stringstream id_num,trial_num;
      id_num << runid;
      trial_num << a;
      string fname = r_conn_path + trial_num.str() + "r" + id_num.str() + txt_suff;
      string fname2 = conn_path + trial_num.str() + "r" + id_num.str() + dat_suff;
      ofstream write_conn_read;
      ofstream write_conn_dat;
      write_conn_read.open(fname.c_str());
      write_conn_dat.open(fname2.c_str());
      connectivity->rank_groups(train_lab,ntrain,group_rank,'s');//group structure
      connectivity->write_groups(write_conn_read,write_conn_dat,group_rank,'s');
      write_conn_read.close();
      write_conn_dat.close();
      */
    }
    

    //micro-variable (end of trial) data file (FORMAT: <n.label> <membrane volt> <exc.cond.> <inh.cond.>)
    if((a%volt_save)==0){
      stringstream id_num, trial_num;
      id_num << runid;
      trial_num << a;
      string fname = sim_base+sim_lab+"/"+volt_e_path + trial_num.str() + "r" + id_num.str() + dat_suff;

      ofstream write_volt;
      write_volt.open(fname.c_str()); 
      for (int i=0; i<SIZE; i++){
	write_volt<<i<<" "<<network[i]->get_pvt_var('v')<<" "<<network[i]->get_pvt_var('e')<<" ";
	write_volt<<network[i]->get_pvt_var('i')<<endl;
      }
      write_volt.close();
    }

    //Save syn weights
    if((a%synapse_save)==0){
      stringstream id_num, trial_num;
      id_num << runid;
      trial_num << a;
      string fname = sim_base+sim_lab+"/"+synapse_path + trial_num.str() + "r" + id_num.str() + dat_suff;

      ofstream write_syn;
      write_syn.open(fname.c_str()); 
      connectivity->write_syns(write_syn);
      write_syn.close();
    }

    //Save spike times for roster plot (FORMAT: <n.label> <group rank (0 if not part of chain)> <spike t>)
    if((a%roster_save)==0){
      stringstream id_num, trial_num;
      id_num << runid;
      trial_num << a;
      string fname = sim_base+sim_lab+"/"+roster_path + trial_num.str() + "r" + id_num.str() + dat_suff;

      //if((a%synapse_save)!=0){//group_rank is required to make a meaningful roster plot
      //	connectivity->rank_groups(train_lab,ntrain,group_rank,'s');
      //}
      ofstream write_roster;
      double *times;
      write_roster.open(fname.c_str()); 
      for (int i=0; i<SIZE; i++){
	for (int j=0; j<network[i]->get_spkhist(times); j++){
	  write_roster<<i<<" "<<times[j]<<endl;
	}
      }
      write_roster.close();
    }
    
    //Reset neuron values for next trial
    for (int i=0; i<SIZE; i++){
      network[i]->resetspike();
      network[i]->resetval();
    }
  
  }//end of simulation
 
  return 0;
}//end of main