// CA1 heteroassociative memory network: Storage and recall
// CA1 PCs, BCs, AACs, BSCs and OLMs (using moderate cell models)
// EC, CA3 (excitatory) and Septal (inhibitory) inputs
// Cycle is: Recall-Storage-Recall etc
// Parallel code adapted from Hines' ran3par.hoc
// VCU & BPG 8-1-09


{load_file("nrngui.hoc")}  // load the GUI and standard run libraries

{load_file("pyramidal_cell4.hoc")}
{load_file("basket_cell17S.hoc")}
{load_file("axoaxonic_cell17S.hoc")}
{load_file("bistratified_cell13S.hoc")}
{load_file("olm_cell2.hoc")}
{load_file("stim_cell.hoc")}
{load_file("burst_cell.hoc")}

{load_file("ranstream.hoc")}  // to give each cell its own sequence generator

objref pc
pc = new ParallelContext()

STARTDEL = 100	// msecs
THETA = 250	// msecs (4 Hz)
GAMMA =25   // (50 Hz)  25	// msecs (40 Hz)
ECCA3DEL = 9	// msecs

if (name_declared("SPATT") == 0) {execute("SPATT= 20")}  // number of active cells per pattern


npcell = 100	
nbcell = 2	
naacell = 1		
nbscell = 1	
nolm = 	1
nCA3 = 100
nEC = SPATT
nSEP = 10

ncell = npcell+nbcell+naacell+nbscell+nolm	// total number of cells
nstim = nCA3+nEC+nSEP		// total number of inputs
ntot = ncell+nstim

// gid ordering: 
//	PCs:0..npcell-1
//	BCs:npcell..npcell+nbcell-1
//	etc
// indices of first cell of each type in list "cells"
iPC = 0
iBC = npcell
iAAC = npcell+nbcell
iBSC = npcell+nbcell+naacell
iOLM = npcell+nbcell+naacell+nbscell
iCA3 = npcell+nbcell+naacell+nbscell+nolm
iEC = npcell+nbcell+naacell+nbscell+nolm+nCA3
iSEP = npcell+nbcell+naacell+nbscell+nolm+nCA3+nEC


//////////////////////////////////////////////////////////////
// Steps 2 and 3 are to create the cells and connect the cells
//////////////////////////////////////////////////////////////

C_P = 1  // probability of excitatory connections received by each CA1 PC
         // from CA3 inputs (1 gives full connectivity)
         
NPATT = 0	// number of stored patterns
CFRAC = 1	// fraction of active cells in cue
iNPPC=1		// index of a non-pattern PC (1st patt in 5 patterns)
CPATT=1
ECPATT=1

strdef FCONN, FSTORE	// file name of connection weights and patterns

//  variables set in the batch file


if (name_declared("MOLT_THETA") == 0) {execute("MOLT_THETA = 100")}	// multiplicator of THETA parameter
SIMDUR = STARTDEL + (THETA*MOLT_THETA)	// simulation duration (msecs)
STORE_NEW_WEIGHTS = 1 //  enables the weight matrix storage, the matrix will be the input of the Recall phase 

if (name_declared("NSTORE") == 0) {execute("NSTORE = 1")}	// indicates the number of patterns to store; 
if (name_declared("SETPATT") == 0) {execute("SETPATT = 1")} //indicates the set of patterns to be used; default is 1 (setpatt1)
if (name_declared("ORT")==0){execute("ORT=0")}  // indicates if orthogonal patterns are presented
if (name_declared("NORT")==0){execute("NORT=1")} // /indicates the set of orthogonal patterns to be used;
if (name_declared("ECWGT") == 0) {execute("ECWGT = 0.001")} //EC versuc PC;
if (name_declared("CLWGT") == 0) {execute("CLWGT = 0.00045")} //cue weight; default is 0.
if (name_declared("CHWGT") == 0) {execute("CHWGT = 0.000738")} //cue weight  control case
//if (name_declared("CHWGT") == 0) {execute("CHWGT = 0.001197")} //cue weight  creb case
if (name_declared("STDPPFAC") == 0) {execute("STDPPFAC = 0.6")} //potentiation factor
if (name_declared("STDPDFAC") == 0) {execute("STDPDFAC = 0.8")} //depression factor
if (name_declared("inputrand") == 0) {execute("inputrand = 6")}
if (name_declared("CREB") == 0) {execute("CREB=0")} // indicates if CREB case runs
if (name_declared("N_PRINT_INFO_CELLS") == 0) {execute("N_PRINT_INFO_CELLS = 2")} 



sprint(FCONN,"Weights/wgtsN%dS%dP0.dat",npcell,SPATT)

if (ORT==0) {
    
    sprint(FSTORE,"Weights/setpatt%d/pattsN%dS%dP%d.dat",SETPATT,npcell,SPATT,NSTORE)

} else {
  sprint(FSTORE,"Weights/setpatt%d/pattsN%dS%dP%do%d.dat",SETPATT,npcell,SPATT,NSTORE,NORT)
}


NDIM=NSTORE


// varibales for STDP window
STDPDMed=-22
STDPDVar=8
STDPPTAU=10

// variables for CA1 template
vinit=-65
somakap =0.007
somakad =0.007
somacaL=0.001
ghsoma = 0.00004 
somacaT=0.0005


// output directory
strdef running_proc,position,position2,s_sys, path, pathNORT, sfSCOPE
if (ORT==0) {
	sprint(path, "MT_%d_NS_%d_SPATT_%d_STDPP_%.2f_STDPD_%.2f_SET_%d_INPR_%d", MOLT_THETA, NSTORE, SPATT, STDPPFAC, STDPDFAC, SETPATT, inputrand)
} else {
	sprint(path, "MT_%d_NS_%d_SPATT_%d_STDPP_%.2f_STDPD_%.2f_SET_%d_NORT_%d_INPR_%d", MOLT_THETA, NSTORE, SPATT, STDPPFAC, STDPDFAC, SETPATT, NORT, inputrand)
}

sprint(running_proc,"bpattrun2")

	if (CREB==0) {
		sprint(position,"Results/%s/%s", running_proc, path)
		sprint(position2,"Results/bpattrun/%s_RECALL",path)
	} else {
		
		sprint(position,"Results/%s/%s_CREB", running_proc, path)
		sprint(position2,"Results/bpattrun/%s_RECALL_CREB", path)
	}


if(pc.id==0) {
	sprint(s_sys,"mkdir -p %s", position)
	system(s_sys)  // to run under Windows write 'WinExec(s_sys)' instead 'system(sys)'

	sprint(s_sys,"mkdir -p %s", position2)
	system(s_sys)  // to run under Windows write 'WinExec(s_sys)' instead 'system(sys)'

	
}

pc.barrier()

if (CREB==0) {

     redsAHP=1
     redmAHP=1
     tauca=3.6
     thna3=-28
     thna3dend=-28
     shift=0    

} else {

     redsAHP=0.36
     redmAHP=0.48
     tauca=3.6
     thna3=-28
     thna3dend=-28
     shift=0
}




// Simple connectivity
CA3_PC = nCA3  // # of connections received by each PC from CA3 cells (excit)
CA3_BC = nCA3  // # of connections received by each BC from CA3 cells (excit)
CA3_AAC = nCA3  // # of connections received by each BC from CA3 cells (excit)
CA3_BSC = nCA3  // # of connections received by each BC from CA3 cells (excit)
EC_PC = nEC  // # of connections received by each PC from EC cells (excit)
EC_BC = nEC  // # of connections received by each BC from EC cells (excit)
EC_AAC = nEC  // # of connections received by each AAC from EC cells (excit)

SEP_BC = nSEP  // # of connections received by each basket cell from septum (inhib)
SEP_AAC = nSEP  // # of connections received by each AAC cell from septum (inhib)
SEP_BSC = nSEP  // # of connections received by each BSC cell from septum (inhib)
SEP_OLM = nSEP  // # of connections received by each OLM cell from septum (inhib)

PC_PC = 1  // # of connections received by each PC from other PCs (excit)
PC_BC = npcell  // # of connections received by each basket cell from PCs (excit)
PC_BSC = npcell  // # of connections received by each bistratified cell from PCs (excit)
PC_AAC = npcell  // # of connections received by each bistratified cell from PCs (excit)
PC_OLM = npcell  // # of connections received by each OLM cell from PCs (excit)

BC_PC = 2  // # of connections received by each PC from basket cells (inhib)
BC_BSC = 2  // # of connections received by each BSC from basket cells (inhib)
BC_OLM = 2  // # of connections received by each OLM from basket cells (inhib)
BC_BC = 1  // # of connections received by each BC from other BCs (inhib)
AAC_PC = 1  // # of connections received by each PC from axoaxonic cells (inhib)
BSC_PC = 1  // # of connections received by each PC from bistratified cells (inhib)
BSC_BC = 1  // # of connections received by each BC from bistratified cells (inhib)
OLM_PC = 1  // # of connections received by each PC from OLM cells (inhib)
OLM_BC = 1  // # of connections received by each basket cell from OLM cells (inhib)

Pcell2Pcell_weight = 0.001
Pcell2Pcell_delay = 1

Bcell2Pcell_weight =0.03 //0.02
Bcell2Pcell_delay = 1
Pcell2Bcell_weight =0.0005 
Pcell2Bcell_delay = 1
Bcell2Bcell_weight = 0.001
Bcell2Bcell_delay = 1
Bcell2BScell_weight =0.02
Bcell2BScell_delay = 1
Bcell2OLMcell_weight = 0.0
Bcell2OLMcell_delay = 1

AAcell2Pcell_weight = 0.05 //0.04
AAcell2Pcell_delay = 1
Pcell2AAcell_weight =0.0005
Pcell2AAcell_delay = 1

BScell2Pcell_weight =0.003 //0.002	
BScell2Pcell_delay = 1
BScell2Pcell_GABAB_weight=0.0004 //0.0004
BScell2Pcell_delay = 1
Pcell2BScell_weight =0.0005
Pcell2BScell_delay = 1
BScell2Bcell_weight = 0.01
BScell2Bcell_delay = 1

OLMcell2Pcell_weight =0.05//0.04
OLMcell2Pcell_GABAB_weight =0.0004
OLMcell2Pcell_delay = 1
OLMcell2Bcell_weight = 0.01
OLMcell2Bcell_delay = 1
Pcell2OLMcell_weight =0.00005
Pcell2OLMcell_delay = 1
OLMcell2Bcell_weight = 0.0
OLMcell2Bcell_delay = 1

// Synapse indices
// onto CA1 PCs
E_EC = 0	// EC AMPA excit to medium SLM (2 of)
E_CA3 = 2	// CA3 AMPA excit to medium SR
EN_CA3 = 3	// CA3 NMDA excit to medium SR
EM_CA3 = 23	// CA3 modifiable (STDP) AMPA excit to medium SR
E_PC = 4	// CA1 recurrent AMPA excit to proximal SR
I_BC = 5	// ff&fb inhib via BCs to soma
I_AAC = 6	// ff&fb inhib via AACs to axon initial segment
I_BSC = 11	// ff&fb inhib via BSCs to SR med (12 of: 6 GABAA, 6 GABAB)
I_OLM = 7	// fb inhib via OLMs to SLM (4 of: 2 GABAA, 2 GABAB)

// onto INs (BC, AAC, BSC)
EI_EC = 0	// EC AMPA excit (2 of; not onto BSC)
EI_CA3 = 2	// CA3 AMPA excit (4 of)
EI_PC = 6	// CA1 PC AMPA excit (2 of)
II_SAME = 8	// inhib from neighbouring INs (BC->BC; BSC->BSC)
II_OPP = 9	// inhib from other INs (BSC->BC; BC->BSC)
II_SEP = 10	// inhib from septum (4 of: 2 GABAA, 2 GABAB)

// onto INs (OLM)
EO_PC = 0	// CA1 PC AMPA excit (2 of)
IO_IN = 2	// inhib from INs and septum (2 of: 1 GABAA, 1 GABAB)

// Septal inhibition
SEPNUM = 100000	// number of SEP spikes
SEPSTART = STARTDEL+(THETA/12)	// time of first SEP spike
SEPINT = 20	// SEP spike ISI (during burst)
SEPNOISE = 0.4	// SEP ISI noise
SEPBINT = 2*THETA/3	// SEP interburst interval
SEPBLEN = THETA/3	// SEP burst length
SEPWGT = 0.02	// SEP weight to BCs and AACs
SEPWGTL = 0.0002 //0.0002	// SEP weight to BSCs and OLMs
SEPWGTL2 = 0.0002 //0.0002	// SEP weight to BSCs and OLMs
SEPDEL = 1	// SEP delay

// Background excitation (not used)
ENUM = 0	// number of spikes
ESTART = 0	// time of first spike
EINT = 200	// spike ISI
ENOISE = 1	// ISI noise
EWGT = 0.001	// excitatory weights (AMPA)
ENWGT = 0.002	// excitatory weights (NMDA)
EDEL = 1	// delay (msecs)

// EC excitation
ECNUM = 100000	// number of EC spikes
ECSTART = STARTDEL	// time of first EC spike
ECINT = GAMMA	// EC spike ISI
ECNOISE = 0.2	// EC ISI noise
//ECWGT = 0.0005	// 0.001// EC weight to PCs
ECWGT0=0.0
ECDEL = 1	// EC delay
ECIWGT =0.00015	// excitatory weights to INs
CA3IWGT =0.00015     //0.00015	// excitatory weights to INs
CA3IWGT2 =0.00015     //0.00015	// excitatory weights to INs

EIDEL = 1	// delay (msecs)

// Cue (CA3) excitation
CNUM = 100000	// number of cue spikes
CSTART = STARTDEL+ECCA3DEL	// time of first cue spike
CINT = GAMMA	// cue spike ISI
CNOISE = 0.2	// cue ISI noise
//CHWGT = 0.0015	// cue weight
//CLWGT = 0  //0.0005	// unlearnt weight (usually 0)
CNWGT = CLWGT //0.0005  // 0.0005	// excitatory weights (NMDA)
CDEL = 1	// cue delay

// STDP configuration


AMPASUPP = 0.4	// fraction of AMPA during storage
STDPTHRESH = -55	// voltage threshold for STDP
STDPSTART = STARTDEL+(THETA/2)	// STDP starts at same time as EC input
STDPINT = THETA/2	// STDP interburst (recall) interval
STDPLEN = THETA/2	// STDP burst (storage) length


connect_random_low_start_ = 1  // low seed for mcell_ran4_init()

objref cells, nclist, ncslist, ncelist, wslist, wtlist, wsEClist, wtEClist  // cells will be a List that holds 
  // all instances of network cells that exist on this host
  // nclist will hold all NetCon instances that exist on this host
  // and connect network spike sources to targets on this host (nclist)
  // ncslist holds NetConns from input spike sources (NetStims)
objref ranlist  // for RandomStreams on this host

objref stims, stimlist, cuelist, EClist	// phasic and tonic cell stimulation

objref gidvec  // to associate gid and position in cells List
  // useful for setting up connections and reporting connectivity


// Make the network
proc mknet() {

   mcell_ran4_init(connect_random_low_start_)
  nclist = new List()
    
  print "Make cells..."
  mkcells()  // create the cells
 
  print "Make inputs..."
  mkinputs()  // create the CA3, EC and septal inputs
    
  print "Connect cells..."
  

  //print "   EC..."
  // EC to BC
  connectcells(nbcell, iBC, nEC, iEC, EC_BC, EI_EC, EI_EC+1, EIDEL, ECIWGT)
  // EC to AAC
  connectcells(naacell, iAAC, nEC, iEC, EC_AAC, EI_EC, EI_EC+1, EIDEL, ECIWGT)

  //print "   CA3..."
  // CA3 to BC
  connectcells(nbcell, iBC, nCA3, iCA3, CA3_BC, EI_CA3, EI_CA3+3, EIDEL, CA3IWGT)
  // CA3 to AAC
  connectcells(naacell, iAAC, nCA3, iCA3, CA3_AAC, EI_CA3, EI_CA3+3, EIDEL, CA3IWGT)
  // CA3 to BSC
  connectcells(nbscell, iBSC, nCA3, iCA3, CA3_BSC, EI_CA3, EI_CA3+3, EIDEL, CA3IWGT2)

  //print "   septum..."
  // SEP to BC
  connectcells(nbcell, iBC, nSEP, iSEP, SEP_BC, II_SEP, II_SEP+1, SEPDEL, SEPWGT)
  // SEP to AAC
  connectcells(naacell, iAAC, nSEP, iSEP, SEP_AAC, II_SEP, II_SEP+1, SEPDEL, SEPWGT)
  // SEP to BSC
  connectcells(nbscell, iBSC, nSEP, iSEP, SEP_BSC, II_SEP, II_SEP+1, SEPDEL, SEPWGTL)
  // SEP to OLM
  connectcells(nolm, iOLM, nSEP, iSEP, SEP_OLM, IO_IN, IO_IN, SEPDEL, SEPWGTL2)

  //print "   PC..."
  // PC to PC
  connectcells(npcell, iPC, npcell, iPC, PC_PC, E_PC, E_PC, Pcell2Pcell_delay, Pcell2Pcell_weight)
  // PC to BC
  connectcells(nbcell, iBC, npcell, iPC, PC_BC, EI_PC, EI_PC+1, Pcell2Bcell_delay, Pcell2Bcell_weight)
  // PC to AAC
  connectcells(naacell, iAAC, npcell, iPC, PC_AAC, EI_PC, EI_PC+1, Pcell2AAcell_delay, Pcell2AAcell_weight)
  // PC to BSC
  connectcells(nbscell, iBSC, npcell, iPC, PC_BSC, EI_PC, EI_PC+1, Pcell2BScell_delay, Pcell2BScell_weight)
  // PC to OLM
  connectcells(nolm, iOLM, npcell, iPC, PC_OLM, EO_PC, EO_PC+1, Pcell2OLMcell_delay, Pcell2OLMcell_weight)

  //print "   INs..."
  // BC to PC
  connectcells(npcell, iPC, nbcell, iBC, BC_PC, I_BC, I_BC, Bcell2Pcell_delay, Bcell2Pcell_weight)
  // BC to BC
  connectcells(nbcell, iBC, nbcell, iBC, BC_BC, II_SAME, II_SAME, Bcell2Bcell_delay, Bcell2Bcell_weight)
  // BC to BSC
  connectcells(nbscell, iBSC, nbcell, iBC, BC_BSC, II_OPP, II_OPP, Bcell2BScell_delay, Bcell2BScell_weight)
  // BC to OLM
  //connectcells(nolm, iOLM, nbcell, iBC, BC_OLM, IO_IN, IO_IN, Bcell2OLMcell_delay, Bcell2OLMcell_weight)
  // AAC to PC
  connectcells(npcell, iPC, naacell, iAAC, AAC_PC, I_AAC, I_AAC, AAcell2Pcell_delay, AAcell2Pcell_weight)
  // BSC to PC
  connectcells(npcell, iPC, nbscell, iBSC, BSC_PC, I_BSC, I_BSC+5, BScell2Pcell_delay, BScell2Pcell_weight)
  connectcells(npcell, iPC, nbscell, iBSC, BSC_PC, I_BSC+6, I_BSC+11, BScell2Pcell_delay, BScell2Pcell_GABAB_weight)
  // BSC to BSC
  //connectcells(nbscell, iBSC, nbscell, iBSC, BSC_BSC, II_SAME, II_SAME, BScell2BScell_delay, BScell2BScell_weight)
  // BSC to BC
  connectcells(nbcell, iBC, nbscell, iBSC, BSC_PC, II_OPP, II_OPP, BScell2Bcell_delay, BScell2Bcell_weight)
  // OLM to PC
  connectcells(npcell, iPC, nolm, iOLM, OLM_PC, I_OLM, I_OLM+1, OLMcell2Pcell_delay, OLMcell2Pcell_weight)
  connectcells(npcell, iPC, nolm, iOLM, OLM_PC, I_OLM+2, I_OLM+3, OLMcell2Pcell_delay, OLMcell2Pcell_GABAB_weight)
  // OLM to BC
  //connectcells(nbcell, iBC, nolm, iOLM, OLM_BC, II_OPP, II_OPP, OLMcell2Bcell_delay, OLMcell2Bcell_weight)

  

  //print "Connect inputs..."
  // EC input to PCs
     connectEC(E_EC,2)
    connectEC_2(FSTORE, $3, NDIM)	// store new pattern
  // CA3 input to PCs
     connectCA3(FCONN, $2, EM_CA3, EN_CA3)	// with modifiable synapses
  
}




// creates the cells and appends them to a List called cells
// argument is the number of cells to be created
proc mkcells() {local i,j  localobj cell, nc, nil
  cells = new List()
  ranlist = new List()
  gidvec = new Vector()
  // each host gets every nhost'th cell,
  // starting from the id of the host
  // and continuing until no more cells are left
  for (i=pc.id; i < ntot; i += pc.nhost) {
    if (i < iBC) {
      cell = new PyramidalCell(redsAHP,redmAHP,v_init,tauca,thna3,thna3dend,somakap,somakad,ghsoma,somacaT,somacaL,shift)
    } else if (i < iAAC) {
      cell = new BasketCell()	// BC
    } else if (i < iBSC) {
      cell = new AACell()	// AAC
    } else if (i < iOLM) {
      cell = new BistratifiedCell()	// BSC
    } else if (i < iCA3) {
      cell = new OLMCell()	// OLM
    } else if (i < iEC) {
      cell = new StimCell()	// CA3 input
    } else if (i < iSEP) {
      cell = new StimCell()	// EC input
    } else {
      cell = new BurstCell()	// Septal input
    }
    cells.append(cell)
    pc.set_gid2node(i, pc.id)  // associate gid i with this host
    nc = cell.connect2target(nil)  // attach spike detector to cell
    pc.cell(i, nc)  // associate gid i with spike detector
    ranlist.append(new RandomStream(inputrand*iSEP+i))  // ranlist.o(i) corresponds to
    gidvec.append(i)
  }
}

// Procedure for CA1 templeate

proc caT_insert() {
    
     for (x) if (x>0 && x<1)  {  
        xdist = distance(x)
       insert cat  
        if (xdist < 300) {
           gcatbar_cat(x) = $1*(1+xdist/60)

        } else {
           gcatbar_cat(x) = $1*6     
    }
}
}

proc caL_insert() {

     cal_distance=200
     for (x) if (x>0 && x<1) {  
         xdist = distance(x)
         insert calH
         mytau_calH=$2
         if (xdist < cal_distance) {
           gcalbar_calH(x) = $1 *(1-xdist/cal_distance)

        } else {
           gcalbar_calH(x) = 0       
    }

                    
     }
}

     

proc A_h_insert(){

     ghend=$1*7
     dhalf=280
     steep=50
     KMA=3
   
	insert h
        ghdbar_h=0
 	insert kap
   	  gkabar_kap = 0
	insert kad
   	 gkabar_kad = 0


	for (x) if (x>0 && x<1) {  
       	xdist=distance(x)

          	ghdbar_h(x)= $1 + (ghend - $1)/(1.0 + exp((dhalf-xdist)/steep))
             
     		 if (xdist < 100){
	         gkabar_kap(x) = $2*(1+KMA*xdist/100)
               vhalfl_h=-73

              	   }else{

                        vhalfl_h=-81
                       gkabar_kad(x) = $3*(1+KMA*xdist/100)
                      // print secname(), " ",xdist, gkabar_kad(x)
                 		         }
       

}
}


proc acc_dist(){local i

 for i=0, cells.count-1 {
    gid = gidvec.x[i]  // id of cell
    if (gid >= 0 && gid < iBC) {    
      
       access cells.object(i).soma
       distance(0,x)

       access cells.object(i).radTprox
        A_h_insert(ghsoma,somakap,somakad)
        caT_insert(somacaT)
        caL_insert(somacaL,tauca)

     access cells.object(i).radTmed
        A_h_insert(ghsoma,somakap,somakad)
        caT_insert(somacaT)
         caL_insert(somacaL,tauca)

      access cells.object(i).radTdist
        A_h_insert(ghsoma,somakap,somakad)
        caT_insert(somacaT)
         caL_insert(somacaL,tauca)

 access cells.object(i).lm_thick1
        A_h_insert(ghsoma,somakap,somakad)
  access cells.object(i).lm_medium1
        A_h_insert(ghsoma,somakap,somakad)
 access cells.object(i).lm_thin1
        A_h_insert(ghsoma,somakap,somakad)

access cells.object(i).lm_thick2
        A_h_insert(ghsoma,somakap,somakad)
 access cells.object(i).lm_medium2
        A_h_insert(ghsoma,somakap,somakad)
 access cells.object(i).lm_thin2
        A_h_insert(ghsoma,somakap,somakad)
  


  }

}

}


// reports distribution of cells across hosts
proc report_gidvecs() { local i, rank
  pc.barrier()  // wait for all hosts to get to this point
  if (pc.id==0) printf("\ngidvecs on the various hosts\n")
  for rank=0, pc.nhost-1 {  // host 0 first, then 1, 2, etc.
    if (rank==pc.id) {
      print "host ", pc.id
      gidvec.printf()
    }
    pc.barrier()  // wait for all hosts to get to this point
  }
}


// sets the CA3, EC and Septal background inputs
proc mkinputs() {local i localobj stim, rs
  for i=0, cells.count-1 {
    gid = gidvec.x[i]	// id of cell
    if (gid >= iCA3 && gid < ntot-nSEP-nEC) {	// appropriate target cell
    // set background activity for excitatory inputs
    stim = cells.object(i).stim
    stim.number = ENUM
    stim.start = ESTART
    stim.interval = EINT
    stim.noise = ENOISE
    }
    if (gid >= iEC && gid < ntot-nSEP) {	// appropriate target cell
    // set background activity for excitatory inputs
    stim = cells.object(i).stim
    stim.number = ENUM
    stim.start = ESTART
    stim.interval = EINT
    stim.noise = ENOISE
    }
    if (gid >= iSEP && gid < ntot) {	// appropriate target cell
    // set background activity for septum
    stim = cells.object(i).stim
    rs = ranlist.object(i)
    stim.number = SEPNUM
    stim.start = SEPSTART
    stim.interval = SEPINT
    stim.noise = SEPNOISE
    stim.burstint = SEPBINT
    stim.burstlen = SEPBLEN
    // Use the gid-specific random generator so random streams are
    // independent of where and how many stims there are.
    stim.noiseFromRandom(rs.r)
    rs.r.negexp(1)
    rs.start()
    }
  }
}


// Target cells receive "convergence" number of inputs from
// the pool of source cells (only one input per source cell at most)
// ("convergence" not reached if no. of sources < convergence)
// connectcells(number of targets, first target cell, 
//		number of source cells, first source cell, 
//		convergence, first synapse,
//		last synapse, connection delay, weight)
// appends the NetCons to a List called nclist

proc connectcells() {local i, j, gid, nsyn, r  localobj syn, nc, rs, u
  // initialize the pseudorandom number generator
  u = new Vector($3)  // for sampling without replacement
 // print "NUMERO ELEMENTI: ", nclist.count, " per il processore ", pc.id
  for i=0, cells.count-1 {	// loop over possible target cells
    gid = gidvec.x[i]	// id of cell
    if (gid >= $2 && gid < $1+$2) {	// appropriate target cell
      rs = ranlist.object(i)  // RandomStream for cells.object(i)
      rs.start()
      rs.r.discunif($4, $4+$3-1)  // return source cell index
      u.fill(0)  // u.x[i]==1 means spike source i has already been chosen
      nsyn = 0
      while (nsyn < $5 && nsyn < $3) {
        r = rs.repick()
        //print "connect rs", r
        // no self-connection and only one connection from any source
        if (r != gidvec.x[i]) if (u.x[r-$4] == 0) {
          // target synapses
          for j = $6, $7 {
            // set up connection from source to target
            syn = cells.object(i).pre_list.object(j)
            nc = pc.gid_connect(r, syn)
            nclist.append(nc)
            nc.delay = $8
            nc.weight = $9
          }
          u.x[r-$4] = 1
          nsyn += 1
        }
      }
    }
  }
}

// connects the EC input layer to PC cells
// read active PCs from pattern file
// all-to-all connectivity between EC and PC pattern cells
// appends the PC NetCons to a List called ncslist
proc connectEC() {local i, j, k, gid localobj syn, src, nc, fp, target
  ncelist = new List()
  wsEClist = new Vector()
  wtEClist = new Vector()

 for i=0, cells.count-1 {	// loop over possible target cells
    gid = gidvec.x[i]	// id of cell
    if (gid >= iPC && gid < npcell+iPC) {	// appropriate target cell
    	if (!pc.gid_exists(gid)) { continue }
		  print "cell ", gid
               target=cells.object(i)
      // target all synapses

      for k = $1, $1+$2-1 {
        syn = target.pre_list.object(k)	// excitatory synapse
        // create pattern stimulus
        for j=0, nEC-1 {
          // set up connection from source to target
          nc = pc.gid_connect(j+iEC, syn)
          nc.delay = ECDEL
          nc.weight = ECWGT0
          ncelist.append(nc)
          wsEClist.append(j+iEC)
	     		wtEClist.append(gid)

          }
      }
    }
  }
  
}



proc connectEC_2() {local i, j,gid, ncue, index  localobj cue, cstim, syn, src, fp, target, file,mm2
  
         
                 fp = new File($s1)
                 fp.ropen()
                 cue = new Vector(npcell)
                 cue.scanf(fp, $2, $3)	// read pattern
                 fp.close()
                 ncue = 0
             
                // find active cells in pattern
                for i=0, cue.size()-1 {
                  if (!pc.gid_exists(i+iPC)) { continue }
                     index=int(i/pc.nhost)
                     //    print "INDEX ", index, " PC.ID ", pc.id, " i ", i

                   for j=0, 2*nEC-1 {
                   			//print "WEIGHT BEFORE :", ncelist.object(index*2*SPATT+j).weight, " i: ", i, " j: ", j 
                       ncelist.object(index*2*SPATT+j).weight = ECWGT0
                   			//print "WEIGHT AFTER :", ncelist.object(index*2*SPATT+j).weight, " i: ", i, " j: ", j
                       if (cue.x[i] == 1) {
                                   // print "Pattern cell ", i
                       			// set up connection from source to target
          			        ncelist.object(index*2*SPATT+j).weight = ECWGT
					      //print "weightEC",  ncelist.object(index*2*SPATT+j).weight[0], "target", i  	

                                  }
                      
                   }
               }
  
}




// connects the CA3 input layer to output cells (PCs and INs)
// read PC connections from a file, with connections to
// a target being a column with index i for target cell i
// appends the PC NetCons to a List called ncslist



proc connectCA3() {local i, j, cp, gid  localobj src, syn, synN, nc, fc, rs, conns, rc
  cp = $2	// connection probability
  mcell_ran4_init(connect_random_low_start_)
  conns = new Vector(nCA3)  // connection weights
  rc = new Vector(nCA3)  // random physical connectivity
  ncslist = new List()
  wslist = new Vector()
  wtlist = new Vector()
  //ncilist = new List()

  // inputs to PCs determined by weight matrix
  for i=0, cells.count-1 {	// loop over possible target cells
    gid = gidvec.x[i]	// id of cell
    if (gid >= iPC && gid < npcell+iPC) {	// appropriate target cell
		  print "cell ", gid
		  syn = cells.object(i).pre_list.object($3)	// AMPA synapse with STDP
		  syn.wmax = CHWGT
		  syn.wmin = CLWGT
		  syn.d = STDPDFAC	// depression
		  syn.p = STDPPFAC	// potentiation
		  syn.ptau = STDPPTAU // potentiation time
              syn.dM = STDPDMed 
              syn.dV = STDPDVar 
		  syn.gscale = AMPASUPP	// fraction of AMPA during storage
		  syn.thresh = STDPTHRESH	// threshold for postsynaptic voltage detection
		  syn.gbdel = STDPSTART
		  syn.gbint = STDPINT
		  syn.gblen = STDPLEN

              synN = cells.object(i).pre_list.object($4)	// NMDA synapse
		  rs = ranlist.object(i)  // the corresponding RandomStream
		  rs.start()
		  rs.r.uniform(0, 1)  // return integer in range 0..1
		  rc.setrand(rs.r)	// generate random connectivity

		  // open connections file
		  fc = new File($s1)
		  fc.ropen()
		  conns.scanf(fc, gid+1, nCA3)	// read incoming weights for cell gid
		  fc.close()
		  for j=0, nCA3-1 {
		    // only connection if physical connection exists
		    if (rc.x[j] <= cp) {
		      //print "   src ", j
		      // set up connection from source to target NMDA synapse
		      nc = pc.gid_connect(j+iCA3, synN)
		      ncslist.append(nc)
			wslist.append(j+iCA3)
		      wtlist.append(gid)
		      nc.delay = CDEL
		      nc.weight = CNWGT	// NMDA weight same for all connections

		     // set up connection from source to AMPA target
			      nc = pc.gid_connect(j+iCA3, syn)
			      ncslist.append(nc)
				wslist.append(j+iCA3)
				wtlist.append(gid)
			      nc.delay = CDEL
			      nc.weight = CLWGT	// unlearned weight
			      				
						
		       }
		  }
    }
  }
}




proc connectCA3_2() {local i, j, k, gid, index  localobj src, syn, nc, fc, conns
	conns = new Vector(nCA3)  // connection weights 
	 // inputs to PCs determined by weight matrix
	for i=0, cells.count-1 {	// loop over possible target cells
		gid = gidvec.x[i]	// id of cell
		if (gid >= iPC && gid < npcell+iPC) {	// appropriate target cell
			print "cell ", gid
		  // open connections file
		  fc = new File($s1)
		  fc.ropen()
		  conns.scanf(fc, gid+1, nCA3)	// read incoming weights for cell gid
		  fc.close()
			k=0
			index=int(gid/pc.nhost)
			//print "ncslist.count ", ncslist.count, " nCA3 ", nCA3
			for(j=1; j<2*nCA3; j=j+2){
			// set up connection from source to target
				ncslist.object(index*2*npcell+j).weight[0]=conns.x[k]
				//print "weight",  ncslist.object(index*200+j).weight[0], " s ", wslist.x[index*2*npcell+j], " t ", wtlist.x[index*2*npcell+j], " gid ", gid, " j ", j, "INDICE ", index*2*npcell+j, " INDEX ", index, " i/pc.nhost ", i/pc.nhost, " i ", i
				k=k+1  	
			}	// for j	            
		}  // if gid
	}  // for i
}




proc cancstim() {local i, j, old_cstim, new_cstim localobj cstim
	
	
	for i=0, cells.count-1 {
		gid = gidvec.x[i]	// id of cell
		if (gid >= iCA3 && gid < nCA3+iCA3) {	// appropriate target cell
			if(!pc.gid_exists(gid)) { continue }   
			cstim = cells.object(i).stim
			old_cstim=cstim.number
			cstim.number = 0
			new_cstim=cells.object(i).stim.number
			
		}
		if (gid >= iEC && gid < nEC+iEC) {	// appropriate target cell
			if(!pc.gid_exists(gid)) { continue }   
			cstim = cells.object(i).stim
			old_cstim=cstim.number
			cstim.number = 0
			new_cstim=cells.object(i).stim.number
			
		}
	}
	
}

strdef s, s_system
objref m, v1, v2


proc aweights() {local i, j, simul, count localobj matrix, f, f3, f4, f_T
	
	simul=$1
	
	f4 = new File()
	sprint(s,"%s/EC_WEIGHTf%d.dat",$s2,pc.id)
	f4.wopen(s)

	
	//PRINT EC synpases
	for(j=1; j<ncelist.count; j=j+1){
		f4.printf("%d\t%d\t%d\t%g\n",SIMDUR,wsEClist.x[j],wtEClist.x[j],ncelist.o(j).weight[0])
	}
	f4.close()

      f3 = new File()
	sprint(s,"%s/AMPA_WEIGHTf%d.dat",$s2,pc.id)
	f3.wopen(s)
	
	//PRINT only AMPA synapses
	for(j=1; j<ncslist.count; j=j+2){
		f3.printf("%d\t%d\t%d\t%g\t%g\t%g\n",SIMDUR,wslist.x[j],wtlist.x[j],(ncslist.o(j).weight[0]+ncslist.o(j).weight[1]),ncslist.o(j).weight[1],ncslist.o(j).weight[0])  
	}
	f3.close()
	

	pc.barrier()
	if (pc.id==0){
		//AMPA, AMPA_WEIGHT, EC_WEIGHT
		//system("cat AMPAf*.dat > AMPAtemp_f.dat")
		sprint(s_system,"cat %s/AMPA_WEIGHTf*.dat > %s/AMPA_WEIGHTtemp_f.dat",$s2,$s2)
		system(s_system)
		sprint(s_system,"cat %s/EC_WEIGHTf*.dat > %s/EC_WEIGHTtemp_f.dat",$s2,$s2)
		system(s_system)
		/*sprint(s_system,"sort -n -k 2 -k 3 AMPAtemp_f.dat > AMPAweights_%d_simul_%d.dat",MOLT_THETA,simul)
		system(s_system)*/
		sprint(s_system,"sort -n -k 2 -k 3 %s/AMPA_WEIGHTtemp_f.dat > %s/AMPA_WEIGHTweights_%d_simul_%d.dat",$s2,$s2,MOLT_THETA,simul)
		system(s_system)
		sprint(s_system,"sort -n -k 2 -k 3 %s/EC_WEIGHTtemp_f.dat > %s/EC_WEIGHTweights_%d_simul_%d.dat",$s2,$s2,MOLT_THETA,simul)
		system(s_system)
		//system("rm AMPAf*.dat")
		sprint(s_system,"rm %s/AMPA_WEIGHTf*.dat",$s2)
		system(s_system)
		sprint(s_system,"rm %s/EC_WEIGHTf*.dat",$s2)
		system(s_system)
		//system("rm AMPAtemp_f.dat")
		sprint(s_system,"rm %s/AMPA_WEIGHTtemp_f.dat",$s2)
		system(s_system)
		sprint(s_system,"rm %s/EC_WEIGHTtemp_f.dat",$s2)
		system(s_system)

		/*sprint(s_system,"mv AMPAweights*.dat %s",$s2 )
		system(s_system)*/
		/*sprint(s_system,"mv AMPA_WEIGHTweights*.dat %s",$s2 )
		system(s_system)*/
		//sprint(s_system,"mv EC_WEIGHTweights*.dat %s",$s2 )
		//system(s_system)
		/*sprint(s_system,"cat %s/AMPAweights_%d_simul_%d.dat > %s/AMPAweights_simul_%d.dat",$s2,MOLT_THETA,simul,$s2,simul)
		system(s_system)*/
		/*sprint(s_system,"cat %s/AMPA_WEIGHTweights_%d_simul_%d.dat > %s/AMPA_WEIGHTweights_simul_%d.dat",$s2,MOLT_THETA,simul,$s2,simul)
		system(s_system)*/
		/*sprint(s_system,"cat %s/EC_WEIGHTweights_%d_simul_%d.dat > %s/EC_WEIGHTweights_simul_%d.dat",$s2,MOLT_THETA,simul,$s2,simul)
		system(s_system)*/

				
		if (STORE_NEW_WEIGHTS == 1) {
			f = new File()
			sprint(s,"%s/AMPA_WEIGHTweights_%d_simul_%d.dat",$s2,MOLT_THETA,simul)
			f.ropen(s)
			m = new Matrix()
			m.scanf(f, npcell*nCA3, 6)
			f.close()
			v2 = m.getcol(3)
						
                  sprint(s,"%s/wgtsN%dS%dP%db_bpattrun2.dat",$s2,npcell,SPATT,simul )
                  
			print "Making file"
			f.wopen(s)
			for i=0, nCA3-1 {  // CA3
				for j=0, npcell-1 {    //CA1
					f.printf("%g ",v2.x[i*nCA3+j])
				}
				f.printf("\n")
			}
			f.close()
			print "End making file"
     
    sprint(s_system,"cp %s/wgtsN%dS%dP%db_bpattrun2.dat %s/wgtsN%dS%dP%db_bpattrun2.dat",$s2,npcell,SPATT,simul,$s3,npcell,SPATT,simul)
              system(s_system)


			
		}
	}
}


mknet(FCONN, C_P,ECPATT)  // go ahead and create the net!
acc_dist()


//////////////////////////////////////////////////
// Instrumentation, i.e. stimulation and recording
//////////////////////////////////////////////////

// setup activity in EC stims
proc mkEC() {local i, necs localobj cstim, rs
  EClist = new Vector()
  necs = 0
  for i=0, cells.count-1 {
    gid = gidvec.x[i]	// id of cell
    if (gid >= iEC && gid < iEC+nEC) {	// appropriate target cell
        // create cue stimulus
        cstim = cells.object(i).stim
    	rs = ranlist.object(i)
        cstim.number = ECNUM
        cstim.start = ECSTART
        cstim.interval = ECINT
        cstim.noise = ECNOISE
        // Use the gid-specific random generator so random streams are
        // independent of where and how many stims there are.
        cstim.noiseFromRandom(rs.r)
        rs.r.normal(0, 1)
        rs.start()
        EClist.append(i)
        necs += 1
    }
  }
}


objref cue, fp

// setup activity pattern in input cue stims
proc mkcue() {local i, j, ncue localobj cstim, target, rs
  cuelist = new Vector()
  // open patterns file
  fp = new File($s1)
  fp.ropen()
  cue = new Vector(nCA3)
  cue.scanf(fp, $2, $4)	// read pattern
//  cue.printf()
  fp.close()
  ncue = 0
  // find active cells in pattern
  for i=0, cue.size()-1 {
    if (!pc.gid_exists(i+iCA3)) { continue }
    if (ncue <= SPATT*$3) { 	// fraction of active cells in cue
      if (cue.x[i] == 1) {
        print "Cue cell ", i
        cstim = pc.gid2cell(i+iCA3)
        for j=0, cells.count-1 {
          if (gidvec.x[j] == i+iCA3) {break}	// find cell index
        }
    	rs = ranlist.object(j)
        // create cue stimulus
        //cstim = target.stim
        cstim.number = CNUM
        cstim.start = CSTART
        cstim.interval = CINT
        cstim.noise = CNOISE
        // Use the gid-specific random generator so random streams are
        // independent of where and how many stims there are.
        cstim.noiseFromRandom(rs.r)
        rs.r.normal(0, 1)
        rs.start()
        cuelist.append(i)
        ncue += 1
      }
    }
  }
//  print "cue size ", ncue
}

// remove activity pattern in input cue stims
proc erasecue() {local i, j localobj cstim
  for i=0, cue.size()-1 {
    if (!pc.gid_exists(i+iCA3)) { continue }
    cstim = pc.gid2cell(i+iCA3)
    cstim.number = 0
  }
}




//mkcue(FSTORE, CPATT, CFRAC, NSTORE)	// cue from new pattern
//mkEC()



// Spike recording
objref tvec, idvec  // will be Vectors that record all spike times (tvec)
        // and the corresponding id numbers of the cells that spiked (idvec)
proc spikerecord() {local i  localobj nc, nil
  tvec = new Vector()
  idvec = new Vector()
  for i=0, cells.count-1 {
    nc = cells.object(i).connect2target(nil)
    nc.record(tvec, idvec, nc.srcgid)
    // the Vector will continue to record spike times
    // even after the NetCon has been destroyed
  }
}

spikerecord()


// Record cell voltage traces
objref pvsoma, pvsr, pvslm  // Vectors that record voltages from pattern PC
objref npvsoma, npvsr, npvslm  // Vectors that record voltages from non-pattern PC
objref vBC, vAAC, vBSC, vOLM  // Vectors that record voltages from INs

objref list_pvsoma, list_pvsr, list_pvslm

proc vrecord2() {local i, gid 

list_pvsoma = new List()
list_pvsr = new List()
list_pvslm = new List()

  for i=0, cells.count-1 {	// loop over possible target cells
    gid = gidvec.x[i]	// id of cell
    if (gid< N_PRINT_INFO_CELLS) {
      pvsoma = new Vector()
      pvsr = new Vector()
      pvslm = new Vector()
      pvsoma.record(&cells.object(i).soma.v(0.5))
      pvsr.record(&cells.object(i).radTmed.v(0.5))
      pvslm.record(&cells.object(i).lm_thick1.v(0.5))
      list_pvsoma.append(pvsoma)
	list_pvsr.append(pvsr)
	list_pvslm.append(pvslm)
    }
    if (gid==iNPPC) {
      npvsoma = new Vector()
      npvsr = new Vector()
      npvslm = new Vector()
      npvsoma.record(&cells.object(i).soma.v(0.5))
      npvsr.record(&cells.object(i).radTmed.v(0.5))
      npvslm.record(&cells.object(i).lm_thick1.v(0.5))
    }
    if (gid==iBC) {
      vBC = new Vector()
      vBC.record(&cells.object(i).soma.v(0.5))
    }
    if (gid==iAAC) {
      vAAC = new Vector()
      vAAC.record(&cells.object(i).soma.v(0.5))
    }
    if (gid==iBSC) {
      vBSC = new Vector()
      vBSC.record(&cells.object(i).soma.v(0.5))
    }
    if (gid==iOLM) {
      vOLM = new Vector()
      vOLM.record(&cells.object(i).soma.v(0.5))
    }
  }
}

vrecord2()


////////////////////////////
// Simulation control
////////////////////////////

strdef fstem//, s_sys,position,position2
tstop = SIMDUR
celsius = 34


proc bpattrun2() { local i
 	
for i=1, NSTORE {
           
		if (pc.id==0) {print "Cue pattern ", i} // print header once
		        
		sprint(fstem, "%s/HAM_P%dR%d",position,NPATT,i)
		if(i!=1) {
			  sprint(FCONN,"%s/wgtsN%dS%dP%db_bpattrun2.dat" ,position,npcell,SPATT, i-1)
			   // EC input to PCs
			  connectEC_2(FSTORE, i, NDIM)    // store new pattern
			  // CA3 input to PCs
                    connectCA3_2(FCONN)    // with modifiable synapses
             
            
		}
		mkEC()
            {mkcue(FSTORE, i, CFRAC, NDIM)}    // cue from stored pattern
		pc.barrier()  // wait for all hosts to get to this point
		//{pc.set_maxstep(10)}
		{pc.set_maxstep(1)}
		stdinit()
		{pc.psolve(tstop)}
		aweights(i,position,position2)
		pc.barrier()
            spikeout()
            vout2()
           //erasecue()
          cancstim()
   
 }
}


////////////////////////////
// Report simulation results
////////////////////////////

objref fo
strdef fno

proc spikeout() { local i  
  pc.barrier()  // wait for all hosts to get to this point
  sprint(fno,"%s_spt_proc%d.dat", fstem, pc.id)
  fo = new File(fno)
  fo.wopen()
  for i=0, tvec.size-1 fo.printf("%g\t %d\n", tvec.x[i], idvec.x[i])
  fo.close()
  pc.barrier()
  if (pc.id==0){
  	sprint(s,"cat %s_spt_proc* > %s_spt_temp.dat", fstem, fstem)
  	system(s)
		sprint(s,"sort -n -k 1 -k 2 %s_spt_temp.dat > %s_spt.dat", fstem, fstem)
		system(s)
		sprint(s,"rm %s*_spt_proc*", fstem)
		system(s)
		sprint(s,"rm %s*_spt_temp*", fstem)
		system(s)
	}
  if (pc.id==0) printf("\ntime\t cell\n")  // print header once
  for rank=0, pc.nhost-1 {  // host 0 first, then 1, 2, etc.
    if (rank==pc.id) {
      for i=0, tvec.size-1 {
        printf("%g\t %d\n", tvec.x[i], idvec.x[i])
      }
      fo.close()
    }
    pc.barrier()  // wait for all hosts to get to this point
  }
}



proc vout2() { local i, j, gid
  
  for j=0, cells.count-1 {	// loop over possible target cells
    gid = gidvec.x[j]	// id of cell
    if (gid< N_PRINT_INFO_CELLS) {
      sprint(fno,"%s_pvsoma_%d.dat", fstem,gid)
      fo = new File(fno)
      fo.wopen()
      for i=0, pvsoma.size-1 {
        fo.printf("%g\n", list_pvsoma.o(j).x[i])
      }
      fo.close()
      sprint(fno,"%s_pvsr_%d.dat", fstem,gid)
      fo = new File(fno)
      fo.wopen()
      for i=0, pvsr.size-1 {
      	fo.printf("%g\n", list_pvsr.o(j).x[i])
      }
      fo.close()
      sprint(fno,"%s_pvslm_%d.dat", fstem,gid)
      fo = new File(fno)
      fo.wopen()
      for i=0, pvslm.size-1 {
      	fo.printf("%g\n", list_pvslm.o(j).x[i])
      }
      fo.close()
    }
    if (gid==iNPPC) {
      sprint(fno,"%s_npvsoma.dat", fstem)
      fo = new File(fno)
      fo.wopen()
      for i=0, npvsoma.size-1 {
        fo.printf("%g\n", npvsoma.x[i])
      }
      fo.close()
      sprint(fno,"%s_npvsr.dat", fstem)
      fo = new File(fno)
      fo.wopen()
      for i=0, npvsr.size-1 {
        fo.printf("%g\n", npvsr.x[i])
      }
      fo.close()
      sprint(fno,"%s_npvslm.dat", fstem)
      fo = new File(fno)
      fo.wopen()
      for i=0, npvslm.size-1 {
        fo.printf("%g\n", npvslm.x[i])
      }
      fo.close()
    }
    if (gid==iBC) {
      sprint(fno,"%s_BC.dat", fstem)
      fo = new File(fno)
      fo.wopen()
      for i=0, vBC.size-1 {
        fo.printf("%g\n", vBC.x[i])
      }
      fo.close()
    }
    if (gid==iAAC) {
      sprint(fno,"%s_AAC.dat", fstem)
      fo = new File(fno)
      fo.wopen()
      for i=0, vAAC.size-1 {
        fo.printf("%g\n", vAAC.x[i])
      }
      fo.close()
    }
    if (gid==iBSC) {
      sprint(fno,"%s_BSC.dat", fstem)
      fo = new File(fno)
      fo.wopen()
      for i=0, vBSC.size-1 {
        fo.printf("%g\n", vBSC.x[i])
      }
      fo.close()
    }
    if (gid==iOLM) {
      sprint(fno,"%s_OLM.dat", fstem)
      fo = new File(fno)
      fo.wopen()
      for i=0, vOLM.size-1 {
        fo.printf("%g\n", vOLM.x[i])
      }
      fo.close()
    }
  }
}


bpattrun2()
{pc.runworker()}
{pc.done()}

quit()