// FLEXINPUT
// This is a collection of scripts for adding flexible input to
// intfcol. It uses the lfpstim function that Bill
// wrote, plus an adaptation of "sgrcells" from col.hoc.
// Version: 2011may19 by cliffk


// POISADD
// This function adds an arbitrary Poisson input to a particular
// population or populations of cells. It calls poistim followed
// by stimadd.
// Usage:
//   poisadd(signal,timei,timef,freq,cellpop,cellprct,cellwt,whichsy)
// where
//   signal describes the probability of a spike at a given time (e.g. a 10K-element sine wave)
//   timei is the start time of the stimulus (in ms, e.g. 2e3)
//   timef is the end time of the stimulus (in ms, e.g. 5e3)
//   freq is the number of spikes (in Hz, e.g. 10) (note: signal.size() must be greater than (timef-timei)*freq!)
//   pop is a vector of cell populations (e.g. [E2,E4,E5])
//   cellwt is the weight given to each spike  (e.g. 1e9)
//   whichsy is the synapse used (e.g. AM2)
// Version: 2011may20

proc poisadd () { local timei,timef,freq,cellprct,cellwt,whichsy,npops,cellstart,cellfinish,pickthiscell localobj pickcell,signal,cellpop,spkoutput
  pickcell=new Random() //
  signal=$o1 // Signal to base the Poisson spike train on
  timei=$2 // Start time of signal
  timef=$3 // End time of signal
  freq=$4 // Frequency/rate of the signal
  cellpop=$o5 // Cell populations to add signal to
  cellprct=$6 // Percent of cells to stimulate in each population
  cellwt=$7 // Weight of each synapse
  whichsy=$8 // Type of each synapse
  npops=cellpop.size() // Number of cell populations
  count=0
  for h=0,numcols-1 { // Loop over columns
    for i=0,npops-1 { // Loop over each cell population
	   	cellstart=col[h].ix[cellpop.x[i]] // Starting cell index
	  	cellfinish=col[h].ixe[cellpop.x[i]] // Finishing cell index
	  	for cellid=cellstart,cellfinish { // Loop over each cell in the population
	  		pickthiscell=100*pickcell.repick() // Whether or not to pick this cell
	  		if(cellprct>pickthiscell) { // Pick out cellprct percent of cells
	  			thisseed=7829*cellid+24091*i+251 // Create a pseudorandom seed
	  			spkoutput=poistim(signal,timei,timef,freq,thisseed) // Calculate Poisson train
	  			stimadd(spkoutput,cellid,cellwt,whichsy)
	  		}
	  	}
  	}
/*  col[h].cstim.pushspks() // Test -- stim wasn't having any effect before*/
  }
}


// POISTIM -- arbitrary Poisson generator
//** spktimevec = poistim(signal,timei,timef,freq)
// signal is vector giving the input signal - eg LFP
// timei gives the initial time time of the signal
// timef gives the final time of the signal, thus timespan is timef-timei
// freq gives the target freq for the spike train -- this is approximate
// Example: 
// objref signal, spktimevec
// signal=new Vector()
// signal.indgen(0.1,0.9,0.001)
// spktimevec=poistim(signal,10,5)
// spktimevec.size() = 50
// Note: the number of points in "signal" must be equal to or greater than the number of spikes!
// Version: 2011may20
obfunc poistim () { local a,timei,timef,thisseed localobj signal,v1,v2,vt
  signal=$o1 timei=$2 timef=$3 freq=$4 thisseed=$5 // Handle input arguments: signal
  a=allocvecs(v1,v2) // Allocate vectors
  vt=new Vector(signal.size) // but ((timef-timei)*freq) is number of spikes desired in period
  vt.setrnd(4,thisseed) // seed for 0-1
  v1.copy(signal) v1.inv()
  vt.mul(v1) // scale the intervals by the signal
  vt.mul((timef-timei)/vt.sum)
  vt.integral() // turn intervals into times
  v1.resize((timef-timei)*freq/1e3)
  v1.setrnd(6,0,vt.size-1,thisseed) // rand unique indices; to cull to get only (maxt*freq/1e3)
  v2.index(vt,v1) // pick the times randomly
  vt.copy(v2)
  vt.add(timei) // Add start time
  dealloc(a)
  return vt
}



// STIMADD -- add stimulus to the input list for a single cell
// This function, based on sgrcells, adds an arbitrary
// stimulus to the rest of the input NQS table vq.
// Usage:
//	 stimadd(times,cellid,cellwt,whichsy)
// where
//	 times is a length-N vector of spike times (e.g. 0, 1.34, 2.53, 7.34, 7.45)
//	 cellid is the cell ID (e.g. 142)
//   cellwt is the synaptic weight (e.g. 1e9)
//   whichsy is the synapse type (e.g. AMPA)
// Version: 2011may20
proc stimadd () { local cellid,cellwt,whichsy,npts,ii,foo localobj times,vqtmp
   if (vq==nil) vq=new NQS("ind","time","cellwt","whichsy") // Initialize NQS to store spikes
   vqtmp=new NQS("ind","time","cellwt","whichsy")
   times=$o1 // Incoming spike times (e.g. 0, 1.34, 2.53, 7.34, 7.45)
   cellid=$2 // Cell ID (e.g. 142)
   cellwt=$3 // Synaptic weights (e.g. 1e9)
   whichsy=$4 // Synapse type (e.g. AMPA)
   npts=times.size() // Find the number of points
   for ii=0,3 vqtmp.v[ii]=new Vector(npts) // Initialize vectors
   vqtmp.v[0].fill(cellid) // Assign the cell ID
   vqtmp.v[1]=times // Assign the times to the second column
   vqtmp.v[2].fill(cellwt) // Assign weights
   vqtmp.v[3].fill(whichsy) // Assign synapse type
   vqtmp.pad() // Shouldn't be necessary, but it is -- make sure all columns are the same size
   vq.append(vqtmp) // Append to original array -- won't take effect until pushspks() call, however
   nqsdel(vqtmp) // Garbage collection
 }