Motor cortex microcircuit simulation based on brain activity mapping (Chadderdon et al. 2014)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:146949
"... We developed a computational model based primarily on a unified set of brain activity mapping studies of mouse M1. The simulation consisted of 775 spiking neurons of 10 cell types with detailed population-to-population connectivity. Static analysis of connectivity with graph-theoretic tools revealed that the corticostriatal population showed strong centrality, suggesting that would provide a network hub. ... By demonstrating the effectiveness of combined static and dynamic analysis, our results show how static brain maps can be related to the results of brain activity mapping."
Reference:
1 . Chadderdon GL, Mohan A, Suter BA, Neymotin SA, Kerr CC, Francis JT, Shepherd GM, Lytton WW (2014) Motor cortex microcircuit simulation based on brain activity mapping. Neural Comput 26:1239-62 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Neocortex;
Cell Type(s): Neocortex V1 pyramidal corticothalamic L6 cell; Neocortex M1 pyramidal intratelencephalic L2-5 cell; Neocortex fast spiking (FS) interneuron; Neocortex spiking regular (RS) neuron; Neocortex spiking low threshold (LTS) neuron;
Channel(s):
Gap Junctions:
Receptor(s): GabaA; AMPA; NMDA;
Gene(s):
Transmitter(s): Gaba; Glutamate;
Simulation Environment: NEURON;
Model Concept(s): Oscillations; Laminar Connectivity;
Implementer(s): Lytton, William [billl at neurosim.downstate.edu]; Neymotin, Sam [samn at neurosim.downstate.edu]; Shepherd, Gordon MG [g-shepherd at northwestern.edu]; Chadderdon, George [gchadder3 at gmail.com]; Kerr, Cliff [cliffk at neurosim.downstate.edu];
Search NeuronDB for information about:  Neocortex V1 pyramidal corticothalamic L6 cell; Neocortex M1 pyramidal intratelencephalic L2-5 cell; GabaA; AMPA; NMDA; Gaba; Glutamate;
/
src
README
infot.mod *
intf6.mod *
intfsw.mod *
matrix.mod
misc.mod *
nstim.mod *
staley.mod *
stats.mod *
vecst.mod *
boxes.hoc *
col.hoc
declist.hoc *
decmat.hoc *
decnqs.hoc *
decvec.hoc *
default.hoc *
drline.hoc *
filtutils.hoc *
gcelldata.hoc
gmgs102.nqs
grvec.hoc *
infot.hoc *
init.hoc
intfsw.hoc *
labels.hoc *
load.py
local.hoc *
main.hoc
misc.h *
miscfuncs.py
network.hoc
neuroplot.py *
nload.hoc
nqs.hoc *
nqsnet.hoc
nrnoc.hoc *
params.hoc
run.hoc
samutils.hoc *
saveoutput.hoc
saveweights.hoc
setup.hoc *
simctrl.hoc *
spkts.hoc *
staley.hoc *
stats.hoc *
stdgui.hoc *
syncode.hoc *
updown.hoc *
wdmaps2.nqs
xgetargs.hoc *
                            
/* $Id: samutils.hoc,v 1.59 2009/12/11 03:56:04 samn Exp $ */

use_samn_utils=0
if(name_declared("INSTALLED_sn")){
  use_samn_utils=1
  install_sn()
} else {
  printf("Warning: couldn't install samnutils.mod\n")
}

//whether to print debug info
mydbg=0


//////////////////////////////////////////////////////////////////////////////////////
//color consts
black=1
red=2
blue=3
green=4
//////////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////////
//time functions



//////////////////////////////////////////////////////////////////////////////////////


//////////////////////////////////////////////////////////////////////////////////////
//drawing functions
//plot and mark $o1 in color $2 using optional xinc = $3
proc plotit () { local xi localobj vv
  vv=$o1 clr=$2
  if(numarg()>2)xi=$3 else xi=1
  vv.plot(g,xi,clr,4)
  vv.mark(g,xi,"O",12,clr,1)
}

//plot $o1 with error bars in $o2, using xvec $o3, and color $4
proc plotite () { local clr localobj vv,ve,vx
  vv=$o1 ve=$o2 vx=$o3 clr=$4
  vv.ploterr(g,vx,ve,15,clr,4)
  vv.plot(g,vx,clr,4)
  vv.mark(g,vx,"O",15,clr,1)
}
//////////////////////////////////////////////////////////////////////////////////////


///////////////////////////////////////////////////////////
//string functions

func strlen(){ localobj strobj
  strobj = new StringFunctions()
  if(argtype(1)==1){
    return strobj.len($o1.s)
  } else if(argtype(1)==2){
    return strobj.len($s1)
  }
  printf("strlen ERRA: invalid argtype\n")
  return 0
}

func strcpy(){
  if(argtype(1)==1 && argtype(2)==1){
    sprint($o1.s,$o2.s)
    return 1
  } else if(argtype(1)==1 && argtype(2)==2){
    sprint($o1.s,$s2)
    return 1
  }
  return 0
}
///////////////////////////////////////////////////////////

//utility functions

func MIN(){
  if($1 <= $2) return $1
  return $2
}

func MAX(){
  if($1 >= $2) return $1
  return $2
}

//return 1 iff int($1) is a power of 2
//$o1 can be list or vector, returns
//true iff their size,count is powof2
func powof2(){ local val

  if(argtype(1)==0){
    val = int($1)
  } else if(isobj($o1,"Vector")){
    val = $o1.size
  } else if(isobj($01,"List")){
    val = $o1.count
  } else {
    printf("powof2 ERRA: invalid arg type\n")
    return -1
  }

  if(val==1) return 1

  if(val%2==1 || val<=0) return 0

  while(val>=2){
    if(val%2==1) return 0
    val /= 2
  }

  return 1
}

//return closest powof2 >= $1
//$o1 can be list or vector, returns
//closest powof2 >= their size,count
func ceilpowof2(){ local ceilval,val

  if(argtype(1)==0){
    val = $1
  } else if(isobj($o1,"Vector")){
    val = $o1.size
  } else if(isobj($o1,"List")){
    val = $o1.count
  } else {
    printf("ceilpowof2 ERRA: invalid arg type\n")
    return -1
  }

  ceilval = 1

  while(ceilval < val) ceilval *= 2

  return ceilval
}

//return closest powof2 <= $1
//$o1 can be list or vector, returns
//closest powof2 <= their size,count
func floorpowof2(){ local floorval,val

  if(argtype(1)==0){
    val = $1
  } else if(isobj($o1,"Vector")){
    val = $o1.size
  } else if(isobj($o1,"List")){
    val = $o1.count
  } else {
    printf("floorpowof2 ERRA: invalid arg type\n")
    return -1
  }

  floorval = val

  while(!powof2(floorval) && floorval >= 2) floorval -= 1

  while(floorval > val) floorval /= 2

  return floorval
}

//discrete wavelet transform using db4
//$o1 = input vector
//$2 = floor ? (if true use floor of pow of 2 of size, else ceil)
//$3 = inverse transform
//returns vector of 2^n dwt coefficients from d4t
obfunc dwt(){ local sz,inv,flr localobj d4t,vt

  d4t = new Vector()

  if(!INSTALLED_d4t_wavelet){
    printf("dwt ERRA: d4t not installed!\n")
    return d4t
  }

  if(!isobj($o1,"Vector")){
    printf("dwt ERRB: arg 1 must be vector!\n")
    return d4t
  }

  if(numarg() > 1) flr = $2 else flr = 1
  if(numarg() > 2) inv = $3 else inv = 1

  if(powof2($o1)){
    vt = $o1
    sz = $o1.size
  } else {
    if(flr) sz = floorpowof2($o1) else sz = ceilpowof2($o1)
    vt = new Vector(sz)
    vt.samp($o1,sz)
  }
 
  d4t = new Vector(sz)

  d4t.d4t(vt,inv,0)

  return d4t
}

//zeroes a level of dwt coefficients (for filtering)
//levels are 0 - (n-1) where 2^n is size of
//dwt coeffs vec
//$o1 = vector of dwt coefficients
//$2 = level to zero out or vec of levels to zero out
func zerodwtlevel(){ local idx,jdx,maxdx,offset,level localobj vlev

  if(!isobj($o1,"Vector")){
    printf("zerodwtlevel ERRA: arg 1 must be vector!\n")
    return 0
  }

  if(argtype(2)==0){
    level = $2
  
    offset = 2^level

    maxdx = 2*offset

    for(idx=offset;idx<maxdx;idx+=1) $o1.x(idx) = 0

  } else {
    vlev = $o2    
    for(idx=0;idx<vlev.size;idx+=1){
      zerodwtlevel($o1,vlev.x(idx))
    }
  }


  return 1
}

//upsample input using steps
//$o1 = input vec
//$2 = output size
obfunc stepupsamp(){ local idx,jdx,kdx,fctr localobj vtmp
  vtmp=new Vector($2)

  if(!isobj($o1,"Vector")){
    printf("stepupsamp ERRA: arg 1 must be Vector!\n")
    return vtmp
  }

  if($o1.size == 0){
    printf("stepupsamp ERRB: arg 1 Vector size is 0!\n")
    return vtmp
  }

  fctr = $2 / $o1.size

  if(fctr < 1){
    printf("stepupsamp ERRC: downsampling not allowed!\n")
    return vtmp
  }
  
  kdx = 0
  for(idx=0;idx<$o1.size;idx+=1){
    for(jdx=0;jdx<fctr;jdx+=1){
      vtmp.x(kdx) = $o1.x(idx)
      kdx += 1
    }
  }

  return vtmp
}

//returns dwt coefficient levels in list of vectors
//$o1 = vector of dwt coefficients
//$2 = resample to size [optional]
//$3 = use step resample
obfunc dwtlevels(){ local nlevels,rsz,idx,stp localobj lv,vtmp,vtmp2

  lv = new List()

  if(!INSTALLED_d4t_wavelet){
    printf("dwtlevels ERRA: d4t not installed!\n")
    return lv
  }

  if(!isobj($o1,"Vector")){
    printf("dwtlevels ERRB: arg 1 must be vector!\n")
    return lv
  }

  if(numarg() > 1) rsz = $2 else rsz = -1
  if(numarg() > 2) stp = $3 else stp = 0

  nlevels = log($o1.size) / log(2)

  for(idx=0;idx<nlevels;idx+=1){
    vtmp = new Vector(2^idx)
    vtmp.d4tlncoeffs($o1,idx)
    if(rsz > 0){
      if(stp){
        vtmp2 = stepupsamp(vtmp,rsz)
      } else {
        vtmp2 = new Vector(rsz)
        vtmp2.samp(vtmp,rsz)
      }
      lv.append(vtmp2)
    } else {
      lv.append(vtmp)
    }
  }

  return lv
}

//takes list of vectors , returns vector with avgs
obfunc GetListAvg(){ local idx,jdx,mx localobj vavg
  mx=0
  for idx=0,$o1.count-1 if($o1.o(idx).size>mx) mx = $o1.o(idx).size
  if(mx==0) return nil else vavg = new Vector(mx)
  for idx=0,$o1.count-1 {
    for jdx=0,$o1.o(idx).size-1{
      vavg.x(jdx) = vavg.x(jdx) + $o1.o(idx).x(jdx)
    }
  }
  vavg.div($o1.count) //
  return vavg
}

//intervals(TRAIN,OUTPUT)
//$o1 spikes/voltages
//$o2 interval vector
func intervals () {
  if ($o1.size<=1) { printf("%s size <2 in intervals()\n",$o1) return 0}
  $o2.deriv($o1,1,1)
  return $o2.size
}

//return derivative of input vector
//output vector has same size as input
//$o1 = input vector
//vp = output vector
obfunc Deriv(){ localobj vp
  vp=new Vector($o1.size)
  if($o1.size < 2){
    printf("Deriv ERRAA: input vec size < 2!\n")
    return vp
  }
  vp.deriv($o1,1,2)
  return vp
}
 
//return list of input vector ($o1),
//and 1st & 2nd derivs 
obfunc Traj(){ localobj lv,vp,vpp
  lv = new List()
  vp = Deriv($o1)
  vpp = Deriv(vp)
  lv.append($o1)
  lv.append(vp)
  lv.append(vpp)
  return lv
}
 
//plot list.o(1) vs list.o(0)
//and list.o(2) vs list.o(1)
//return list of graphs
obfunc PlotTraj(){ localobj lv,vp,vpp,lg
  lg = new List()
  lv = $o1
  vp = lv.o(1)
  vpp = lv.o(2)
  vp.label("v'")
  vpp.label("v''")
  lg.append(new Graph())
  vp.plot(lg.o(lg.count-1),lv.o(0))
  lg.append(new Graph())
  vpp.plot(lg.o(lg.count-1),vp)
  return lg
}

//sets up recording of spike times
//$o1 = cell
//$o2 = printlist obj
//$s3 = string
//$4 = spike thresh
proc RecordCell(){
  if(numarg()==3){
    $o2 = new_printlist_nc($o1,0,$s3)
  } else if(numarg()==4){
    $o2 = new_printlist_nc($o1,0,$s3,$4)
  }
}

//$o1 = vec
//$2 = size
//allocates memory for vector but resizes it to 0
proc VecMalloc(){
  $o1 = new Vector()
  $o1.resize($2)
  $o1.resize(0)
}

proc PrintBurstInfo(){
  if($o1.size < 4){
    printf("invalid size for burst info vector\n")
    return
  }
  printf("num bursts = %d\n",$o1.x(0))
  printf("spikes per burst = %d\n",$o1.x(1))
  printf("inter-spike time = %g\n",$o1.x(2))
  printf("inter-burst time = %g\n",$o1.x(3))
}

//$o1 = spike times
//$o2 = output vector[num_bursts,spikes_per_burst,inter_spike_time,inter_burst_time]
func GetBurstInfo(){ local i,lg_thresh,sm_thresh,spikes_per_burst localobj v_inter_spike,v_inter_burst,v_intervals

  lg_thresh = 20   //min time btwn bursts
  sm_thresh = 3.55  //max time btwn spikes within a burst

  if($o1.size < 2){
    $o2 = new Vector()
    $o2.resize(4)
    $o2.x(0)=$o2.x(1)=$o2.x(2)=$o2.x(3)=0
    printf("need at least two spikes to get burst info\n")
    return 0
  }

  v_intervals = new Vector() intervals($o1,v_intervals)
  //printf("here are the intervals:\n")
  //v_intervals.printf

  VecMalloc(v_inter_spike,v_intervals.size)
  VecMalloc(v_inter_burst,v_intervals.size)

  for(i=0;i<v_intervals.size;i+=1){
    if(v_intervals.x(i) >= lg_thresh){
      //printf("\tinter-burst time = %g\n",v_intervals.x(i))
      v_inter_burst.append(v_intervals.x(i))
    } else if(v_intervals.x(i) <= sm_thresh) {
      //printf("\tinter-spike time = %g\n",v_intervals.x(i))
      v_inter_spike.append(v_intervals.x(i))
    }
  }

  $o2 = new Vector() $o2.resize(4)

  //if no valid inter-spike / inter-burst found, return 0
  if(v_inter_spike.size == 0 && v_inter_burst.size == 0) {
    printf("\tno valid inter-spike / inter-burst found\n")
    return 0
  }

  //if( ((v_inter_burst.max-v_inter_burst.min) / v_inter_burst.max) >= 0.3){
    //printf("\tno consistent inter-burst time\n")
  //  return 0
  //}

  $o2.x(0) = v_inter_burst.size() + 1 //# of bursts
  $o2.x(1) = $o1.size / $o2.x(0)  //# of spikes per burst, assumes they are equally distributed

  if(v_inter_spike.size > 0) { $o2.x(2) = v_inter_spike.median() }  //inter-spike time
  if(v_inter_burst.size > 0) { $o2.x(3) = v_inter_burst.median() }  //inter-burst time

  //if no valid inter-spike / inter-burst found, return 0
  //if($o2.x(2) == 0 && $o2.x(3) == 0) {
  //  printf("\tno valid inter-spike / inter-burst found\n")
  //  return 0
  //}

  return 1
}

//$o1 = stim
//$2 = min amp
//$3 = max amp
//$4 = amp inc
//$o5 = vector of # of bursts
//$o6 = vector of # of spikes per burst
//$o7 = inter-spike times
//$o8 = inter-burst times
//$o9 = cell
func GetBurstInfoVecs(){ local curamp,minamp,maxamp,ampinc localobj cell_rec,v_burst
   RecordCell($o9,cell_rec,"cell.soma.v")
   minamp = $2 maxamp = $3 ampinc = $4
   $o5 = new Vector()
   $o6 = new Vector()
   $o7 = new Vector()
   $o8 = new Vector()
   v_burst = new Vector()
   $o1.del = 40
   $o1.dur = 400
   tstop = 500
   for(curamp=minamp;curamp<=maxamp;curamp+=ampinc){
     $o1.amp = curamp
     printf("getting burst info for amp = %g\n",curamp)
     run()
     GetBurstInfo(cell_rec.tvec,v_burst)
     $o5.append(v_burst.x(0))
     $o6.append(v_burst.x(1))
     $o7.append(v_burst.x(2))
     $o8.append(v_burst.x(3))
   }
   return 1
}

//$o1 = graph
//$o2 = vector of # of bursts
//$o3 = vector of # of spikes per burst
//$o4 = inter-spike times
//$o5 = inter-burst times
//$o6 = amp
//$7 , whether to label
//$8 , whether to make new graph
proc PlotBurstInfoVecs(){
  if($8){
    $o1 = new Graph()
  }

  if($7){
    $o2.label("num bursts")
    $o3.label("spikes per burst")
    $o4.label("inter-spike time")
    $o5.label("inter-burst time")
  }

  $o3.plot($o1,$o6,1,1) //color brush
  $o5.plot($o1,$o6,2,1)
  $o2.plot($o1,$o6,3,1)
  $o4.plot($o1,$o6,4,1)

//  $o3.mark($o1,$o6,"+",6,1,1) //style size color brush
//  $o5.mark($o1,$o6,"t",6,2,1)
//  $o2.mark($o1,$o6,"s",6,3,1)
//  $o4.mark($o1,$o6,"o",6,4,1)

  $o1.flush()
}

//$o1 = voltage vector (x-axis=time, y-axis=voltage)
//$o2 = output spike times
//$3 = threshold for spike
//$4 = minimum interspike time
proc GetSpikeTimes() { local ii,v,dipthresh
  $o2 = new Vector($o1.size)
  $o2.resize(0)
  dipthresh = 1.5 * $3
  for(ii=0;ii<$o1.size;ii+=1) {
    v = $o1.x(ii)
    if(v >= $3) { //is v > threshold?
      if($o2.size>0) { //make sure at least $4 time has passed
        if(dt*ii-$o2.x($o2.size-1) < $4) {
          continue
        }
      }
      while( ii+1<$o1.size) { //get peak of spike
        if( $o1.x(ii) > $o1.x(ii+1) ) {
          break
        }
        ii += 1
      }
      $o2.append(dt*ii) //store spike location

      while(ii<$o1.size) { //must dip down sufficiently
        if($o1.x(ii) <= dipthresh) {
          break
        }
        ii += 1
      }
    }
  }
}

//just moves 1.25 ms to left & right of spikes
//returns vector with 1 if its in spike time window
//and 0 otherwise
//$o1 = voltage vector
//$o2 = spike times
//$3 = 1/2 time window around spikes [optional, default=1.25ms]
obfunc GetSpikeWindow(){ local leftv,rightv,meanv,dthresh,left,right,w,idx,jdx,mint,st,twin2\
                    localobj vspikes,vvolt,vns

  vvolt = $o1
  vspikes = $o2

  if(numarg() > 2) twin2 = $3 else twin2 = 1.25

  vns = new Vector(vvolt.size)
//  vns.copy(vvolt)

  for(idx = 0; idx < $o2.size; idx+=1){

    st = vspikes.x(idx)

    left = st/dt - twin2/dt
    if(left < 0) left = 0
    right = st/dt + twin2/dt
    if(right >= vvolt.size) right = vvolt.size-1
 
    w = right - left
    if(w == 0) w = 0

    leftv = vvolt.x(left)
    rightv = vvolt.x(right)

    for(jdx=left;jdx<=right;jdx+=1){
      vns.x(jdx) = 1      
//      vns.x(jdx) = ((w-(jdx-left))*leftv)/w + ((w-(right-jdx))*rightv)/w
    }   
  }

  return vns
}



//$o1 = voltage vector
//$o2 = spike times
obfunc CutOutSpikes(){ local leftv,rightv,meanv,dthresh,left,right,w,idx,jdx,mint,st\
                    localobj vspikes,vvolt,vns,vderiv

  vvolt = $o1
  vspikes = $o2

  vderiv = Deriv(vvolt)

  meanv = vderiv.mean
  dthresh = meanv+vderiv.stdev/4

  printf("dthresh = %f\n",dthresh)

  vns = new Vector(vvolt.size)
  vns.copy(vvolt)

  for(idx = 0; idx < $o2.size; idx+=1){

    st = vspikes.x(idx)

    left = st/dt - 5
    if(left < 0) left = 0
    right = st/dt + 5
    if(right >= vderiv.size) right = vderiv.size-1
 
    while(left > 0){
      if(vderiv.x(left) < dthresh){
        break
      }
      left-=1
    }

    while(right < vvolt.size-1){

      if(abs(vderiv.x(right)) < dthresh){
        break
      }
      right+=1
    }

    w = right - left
    if(w == 0) w = 0

    leftv = vvolt.x(left)
    rightv = vvolt.x(right)

    for(jdx=left;jdx<=right;jdx+=1){
//      vns.x(jdx) = 60      
      vns.x(jdx) = ((w-(jdx-left))*leftv)/w + ((w-(right-jdx))*rightv)/w
    }   
  }

  return vns
}

//$o1 = voltage vector (x-axis=time, y-axis=voltage)
//$o2 = spike time vector
//$o3 = graph obj
proc PlotSpikeTimes(){ local idx
  $o1.plot($o3,dt)
  for(idx=0;idx<$o2.size;idx+=1){
    $o3.mark($o2.x(idx),$o1.x($o2.x(idx)/dt),1,4,red,1)
  }
}

//get interspike interval times
//$o1 = spike train
//$o2 = output ISI start times
//$o3 = output ISI end times
func GetISITimes(){ local idx,mindist,tmp
  mindist = 1.25 //min dist after & before spike to go
  $o2 = new Vector()
  $o3 = new Vector()
  if($o1.size < 2){
    printf("need at least two spikes to get ISI\n")
    return 0
  }
  for(idx=0;idx<$o1.size-1;idx+=1){
    $o2.append($o1.x(idx)+mindist)
    $o3.append($o1.x(idx+1)-mindist)
    //swap if incorrect order due to mindist requirement
    if($o3.x($o3.size-1) <= $o2.x($o2.size-1)){
      tmp = $o3.x($o3.size-1)
      $o3.x($o3.size-1) = $o2.x($o2.size-1) + dt
      $o2.x($o2.size-1) = tmp
    }
  }
  return 1
}

//get trajectory v,v',v''
//$o1 = voltage vector
//$2 = start time
//$3 = end time
//$o4 = output copy of v
//$o5 = 1st deriv (v')
//$o6 = 2nd deriv (v'')
func GetTrajectory(){ local startt,endt
  $o4 = new Vector()  $o5 = new Vector() $o6 = new Vector()

  startt = $2
  endt = $3

  if(startt >= endt){
    printf("GetTrajectory ERRA: invalid time range (%f,%f)!\n",$2,$3)
    return 0
  }

  endt = endt / dt
  startt = startt / dt

  $o4.copy($o1,startt,endt) //copy v
  if($o4.size < 2){
    $o4.resize(0)
    printf("GetTrajectory ERRB: time range %f ms too small!\n",(endt-startt))
    return 0
  }

  $o5 = Deriv($o4)
  $o6 = Deriv($o5)
 
  return 1
}

//get all trajectories from start/end time vecs
//$o1 = voltage vector
//$o2 = start times
//$o3 = end times
//$o4 = list of v
//$o5 = list of v'
//$o6 = list of v''
func GetTrajectories(){ local idx,allbad localobj v1,v2,v3
  allbad = 1
  for(idx=0;idx<$o2.size;idx+=1){
    v1 = new Vector() v2 = new Vector() v3 = new Vector()
    if(GetTrajectory($o1,$o2.x(idx),$o3.x(idx),v1,v2,v3)){
      $o4.append(v1) $o5.append(v2) $o6.append(v3)
      allbad = 0
    } 
  }
  return !allbad
}

//$o1 = graph1
//$o2 = graph2
//$o3 = v
//$o4 = v'
//$o5 = v''
func PlotTrajectory(){
  if(!isassigned($o1)) $o1 = new Graph()
  $o4.label("v'")
  $o4.plot($o1,$o3)

  if(!isassigned($o2)) $o2 = new Graph()
  $o5.label("v''")
  $o5.plot($o2,$o4)

  return 1
}

//$o1 = list of v
//$o2 = list of v'
//$o3 = list of v''
//$o4 = list of graphs
func PlotTrajectories(){ local idx localobj mygv,mygvv
  for(idx=0;idx<$o1.count;idx+=1){
    mygv=new Graph()
    mygvv=new Graph()
    PlotTrajectory(mygv,mygvv,$o1.o(idx),$o2.o(idx),$o3.o(idx))
    $o4.append(mygv)
    $o4.append(mygvv)
  }
  return 1
}

//euclid dist btwn 2 points in 3D space
func EuclidDist(){ local x1,y1,z1,x2,y2,z2
  x1=$1 y1=$2 z1=$3 x2=$4 y2=$5 z2=$6
  return sqrt( (x1-x2)^2 + (y1-y2)^2 + (z1-z2)^2 )
}

//$o1 = vec of x1
//$o2 = vec of y1
//$o3 = vec of z1
//$o4 = vec of x2
//$o5 = vec of y2
//$o6 = vec of z2
func VEuclidDist(){ localobj l1,l2
  if(!use_samn_utils){
    printf("VEuclidDist ERRA: samnutils.mod not available\n")
    return -1
  }
  l1=new List() l2=new List()
  l1.append($o1) l1.append($o2) l1.append($o3)
  l2.append($o4) l2.append($o5) l2.append($o6)
  return LDist_sn(l1,l2)//in /usr/site/nrniv/local/mod/samnutils.mod
}

//$o1 = vec of x1
//$o2 = vec of y1
//$o3 = vec of z1
//$o4 = vec of x2
//$o5 = vec of y2
//$o6 = vec of z2
func VSQDist(){ localobj l1,l2
  if(!use_samn_utils){
    printf("VSQDist ERRA: samnutils.mod not available\n")
    return -1
  }
  l1=new List() l2=new List()
  l1.append($o1) l1.append($o2) l1.append($o3)
  l2.append($o4) l2.append($o5) l2.append($o6)
  return LDist_sn(l1,l2,SQDIFF_sn)//in /usr/site/nrniv/local/mod/samnutils.mod
}

//resample a vec to new size using linear interpolation
//$o1 = vec
//$2 = new size
func Resample(){ local newsz,idxdest,idxsrc,val,fctr,frac,last,lastset localobj vtmp
  vtmp = new Vector($2)
  fctr = $o1.size / $2
  vtmp.x(0) = $o1.x(0)
  idxsrc = fctr
  for(idxdest=1;idxdest<$2-1;idxdest+=1){
    idxsrc = idxdest * fctr
    frac = idxsrc - int(idxsrc)
    idxsrc = int(idxsrc)
    if(idxsrc+1>=$o1.size){
      vtmp.x(idxdest) = $o1.x(idxsrc)
      continue
    }
    val = (1-frac) * $o1.x(idxsrc) + frac * $o1.x(idxsrc+1)
    vtmp.x(idxdest) = val
  }
  vtmp.x($2-1) = $o1.x($o1.size-1)
  $o1.resize($2)
  $o1.copy(vtmp)
  return 1
}

//$o1 = stim
//$o2 = voltage vector (x-axis = time , y-axis = volt)
//$3 = min amp
//$4 = max amp
//$5 = amp inc
//$6 = threshold for spike
//$o7 = amp vec
//$o8 = f vec
//$o9 = cell
proc GetIFVecs() { local minamp,maxamp,ampinc,totalt,curamp,thresh
  minamp = $3 maxamp = $4 ampinc = $5 thresh = $6
  $o7 = new Vector()
  $o8 = new Vector()
  totalt = $o1.dur / 1000
  for(curamp = minamp; curamp <= maxamp; curamp += ampinc) {
    $o1.amp = curamp
    RecordCell($o9,$o2,"soma.v",$6)
    run()
    $o7.append(curamp)
    $o8.append($o2.tvec.size * (1/totalt))
  }
}

proc PlotIFVec(){ localobj g
  g = new Graph()
  $o1.label("Hz")
  $o2.label("nA")
  $o1.plot(g,$o2)
}

//$s1 = output file name
//$o2 = amp vec
//$o3 = freq vec
func SaveIFCurve(){ localobj f
  f = new File()
  f.wopen($s1)
  if(!f.isopen()) {
    printf("couldn't open %s for writing\n",$s1)
    return 0
  }
  if(!$o2.vwrite(f)){
    printf("couldn't write amp vec to %s\n",$s1)
    return 0
  }
  if(!$o3.vwrite(f)){
    printf("couldn't write freq vec to %s\n",$s1)
    return 0
  }
  f.close()
  return 1
}

//$s1 = input file name
//$o2 = amp vec
//$o3 = freq vec
func ReadIFCurve(){ localobj f
  f = new File()
  f.ropen($s1)
  if(!f.isopen()) {
    printf("couldn't open %s for reading\n",$s1)
    return 0
  }
  $o2 = new Vector()
  if(!$o2.vread(f)){
    printf("couldn't read amp vec from %s\n",$s1)
    return 0
  }
  $o3 = new Vector()
  if(!$o3.vread(f)){
    printf("couldn't read freq vec from %s\n",$s1)
    return 0
  }
  f.close()
  return 1
}

//writes vecs to file
//$s1 = path to file
//$o2 ... $on = vectors to be written to file
func WriteVecs(){ local i localobj f
  f = new File()
  f.wopen($s1)
  if(!f.isopen()) {
    printf("couldn't open %s for writing\n",$s1)
    return 0
  }
  for(i=2;i<=numarg();i+=1){
    printf("writing vector %d\n",i-1)
    if(!$oi.vwrite(f)){
      printf("couldn't write vector %d to %s\n",i-1,$s1)
      return 0
    }
  }
  f.close()
  return 1
}

//reads vecs from file
//$s1 = path to file
//$o2 ... $on = vectors to be read from file
func ReadVecs(){ local i localobj f
  f = new File()
  f.ropen($s1)
  if(!f.isopen()) {
    printf("couldn't open %s for reading\n",$s1)
    return 0
  }
  for(i=2;i<=numarg();i+=1){
    printf("reading vector %d\n",i-1)
    if(!$oi.vread(f)){
      printf("couldn't read vector %d from %s\n",i-1,$s1)
      return 0
    }
  }
  f.close()
  return 1
}

//$s1 = path to file
//$2 = index of vector to read
//$o3 = vector to read to
func ReadVecX(){ local idx localobj f
  f = new File()
  f.ropen($s1)
  if(!f.isopen()) {
    printf("couldn't open %s for reading\n",$s1)
    return 0
  }
  if(!isassigned($o3)) $o3 = new Vector()
  for(idx=0;idx<=$2;idx+=1){
    if(f.eof()){
      f.close()
      return 0
    }
    $o3.vread(f)
  }
  f.close()
  return 1
}

//tests whether file exists by opening in 'r' mode
func FileExists(){ localobj f
  f = new File()
  f.ropen($s1)
  if(f.isopen()){
    f.close()
    return 1
  }
  return 0
}

//reads vecs from file
//$s1 = path to file
//$o2 ... $on = vectors to be read from file
/*func Read2DVec(){ local i localobj f,cv
  f = new File()
  f.ropen($s1)
  if(!f.isopen()) {
    printf("couldn't open %s for reading\n",$s1)
    return 0
  }

  //first vector just has # of vectors in file
  cv = new Vector()
  if(!cv.vread(f)){
    printf("couldn't get # vectors in file\n")
    return 0
  }
  printf("there are %d vecs in file\n",cv.size)
  
  for(i=0;i<cv.size;i+=1){
    printf("reading vector %d\n",i+1)
    $3[i] = new Vector()
    if(!$3[i].vread(f)){
      printf("couldn't read vector %d from %s\n",i+1,$s1)
      return 0
    }
  }
  f.close()
  return 1
}*/

//$s1 = file path
//$o2 = vector of # of bursts
//$o3 = vector of # of spikes per burst
//$o4 = inter-spike times
//$o5 = inter-burst times
//$o6 = amp
func WriteBurstInfoVecs(){
  return WriteVecs($s1,$o2,$o3,$o4,$o5,$o6)
}

//$s1 = file path
//$o2 = vector of # of bursts
//$o3 = vector of # of spikes per burst
//$o4 = inter-spike times
//$o5 = inter-burst times
//$o6 = amp
func ReadBurstInfoVecs(){
  return ReadVecs($s1,$o2,$o3,$o4,$o5,$o6)
}

///////////////////////
//random related funcs

my_ran_ind = 0

objref r_uniform

proc InitRand(){
  r_uniform = new Random(1234567891011121314)
  my_ran_ind = r_uniform.MCellRan4()
  r_uniform.uniform(0,1)
}

func uniform(){
  return r_uniform.repick()
}

//$1 - mean
//$2 - stddev
white_noise_N = 20
sqrt_val = sqrt(12/white_noise_N)
func GetWhiteNoise(){ local X,N,U,i,mu,std
  mu=$1
  std=$2

  N=white_noise_N
  X=0

  for(i=0;i<N;i+=1){
     U = uniform()
     X +=  U
  }
	
  // for uniform randoms in [0,1], mu = 0.5 and var = 1/12 
  // adjust X so mu = 0 and var = 1 
	
  X  -=  N/2                // set mean to 0 
  X *= sqrt_val           // adjust variance to 1 X=X*sqrt(12/N)

  return mu + std*X
}

InitRand()

//init rand vector of size $1
//with mean = $2 , stddev = $3
//$o4 is vector
proc WhiteNoiseVec(){ local mu,std,inc,idx
  mu=$2 std=$3
  $o4 = new Vector($1)
  for(idx=0;idx<$o4.size;idx+=1){
    $o4.x(idx) = GetWhiteNoise(mu,std)
  }
}

proc SaveListOfVecs(){ local idx
  for(idx=0;idx<$o1.count;idx+=1){
    savevec($o1.o(idx))
  }
}

proc Smooth() { local idx,w,val,jdx,cnt localobj myv,myvtmp
  myv=$o1
  if(numarg()>1) w=$2 else w=5
  myvtmp=new Vector(myv.size)
  myvtmp.copy(myv)
  for(idx=w;idx<myv.size-w;idx+=1){
    cnt=0
    val=0
    for(jdx=idx-w;jdx<idx+w;jdx+=1){
      val += myv.x(jdx)
      cnt += 1
    }
    if(cnt>0) myvtmp.x(idx) = val / cnt
  }
  myv.copy(myvtmp)
}

proc Smooth2(){ local idx,w,val,jdx,cnt,wght,taus localobj myv,myvtmp
  myv=$o1
  w = taus/2
  wght=0
  myvtmp=new Vector(myv.size)
  if(numarg()>1) taus=3/$2 else taus=3/dt
//  myvtmp.copy(myv)
  for(idx=0;idx<2*w;idx+=1) wght += exp(-abs(idx-w)/taus)
  for(idx=w;idx<myv.size-w;idx+=1){
    cnt=0
    val=0
    for(jdx=idx-w;jdx<idx+w;jdx+=1){
      val += myv.x(jdx) * exp(-abs(jdx-idx)/taus)
      cnt += 1
    }
//    myvtmp.x(idx)=val
    if(cnt>0) myvtmp.x(idx) = val / wght//cnt
//    if(cnt>0) myvtmp.x(idx) = val / exp(1)
  }
  myv.copy(myvtmp)
}

//does smoothing of 1d vector using exp kernel
//$o1 = input vector , (also output)
//$2 = 1/2 width of kernel
//$3 = tstep
proc ExpSmooth(){ local idx,w,val,jdx,wght,myt localobj myv,myvtmp
  myv = $o1
  if(numarg()>1) w = $2 else w = 4
  if(numarg()>2) tstep = $3 else tstep=1
  wght = 0
  myt=0
  myvtmp=new Vector(myv.size)
  myt = tstep
  for(idx=0;idx<w;idx+=1){  
    wght += exp(-abs(tstep))
    myt += tstep
  }
  wght = 2*wght + 1
  if(mydbg) printf("weight=%f\n",mydbg)
  for(idx=w;idx<myv.size-w;idx+=1){
    val=0
    myt=w*tstep
    for(jdx=idx-w;jdx<idx;jdx+=1){
      val += myv.x(jdx) * exp(-myt)
      myt-=tstep
    }
    val += myv.x(jdx)
    myt=tstep
    for(;jdx<idx+w;jdx+=1){
      val += myv.x(jdx) * exp(-myt)
      myt+=tstep
    }
    myvtmp.x(idx)=val/wght
  }
  myv.copy(myvtmp)
}

//fix points from EEG trace where voltage diff
//between sample and median is > thresh since it's probably
//a problem with recording equipment ... assumes only 1 point
//in a row will be an outlier (since it replaces point with
//average of surrounding two neighbors)
//$o1 = vector
//$2 = thresh
//$3 = whether to only fix endpoints
//returns -1 on error, 1 if fixed any, 0 otherwise
func FixEEGOutliers(){ local idx,thresh,med,ept,fixedany localobj vtmp

  if(numarg() < 1){
    printf("FixEEGOutliers ERRA: incorrect # of args\n")
    return -1
  }

  if(!isobj($o1,"Vector")){
    printf("FixEEGOutliers ERRB: require vector arg\n")
    return -1
  }

  vtmp = $o1

  med = $o1.median

  thresh = 20000

  fixedany = 0

  if(numarg() > 1) if(argtype(2)==0) thresh = $2

  ept = 0//default don't ONLY fix endpoints
  if(numarg() > 2) if(argtype(3)==0) ept = $3

  if(mydbg) printf("FixEEGOutliers: thresh = %d\n",thresh)

  for(idx=0;idx<vtmp.size;idx+=1){
    if(abs(vtmp.x(idx)-med) >=thresh){

      fixedany = 1

      printf("FixEEGOutliers: fixing val=%f at idx=%d, ",vtmp.x(idx),idx)

      //take average of surrounding two points
      if(idx + 1 < vtmp.size && idx > 0){
        vtmp.x(idx) = ( vtmp.x(idx-1) + vtmp.x(idx+1) ) / 2
      } else if(idx > 0){ //take value from left
        vtmp.x(idx) = vtmp.x(idx - 1)
      } else if(idx + 1 < vtmp.size){ //take value from right
        vtmp.x(idx) = vtmp.x(idx+1)
      }
      printf("new val=%f, med=%f\n",vtmp.x(idx),med)
    }
    if(ept && idx==0) idx = vtmp.size - 2
  }

  return fixedany
}

///////////////////////////////////////////////////////////////////////////
//////////////////////// filter functions /////////////////////////////////
//$o1 = filter vector
//$2 = size
//$3 = tau
//$4 = inc
proc AlphaFilter(){ local idx,inc,tau,sz,myt
  sz = $2
  tau = $3
  myt = inc = $4
  $o1 = new Vector(sz)
  printf("tau=%f inc=%f\n",tau,inc)
  for(idx=0;idx<sz;idx+=1){
    $o1.x(idx) = myt * exp(-myt/tau)
    myt += inc
  }
  $o1.div($o1.sum)
}

//$o1 = filter vector
//$2 = size of ~ 1/2 filter
//$3 = tau
//$4 = inc
proc SymmetricAlphaFilter(){ local idx,inc,tau,sz localobj vtmp,vtmpr
  sz=$2
  tau=$3
  inc=$4

  vtmp=new Vector()
  AlphaFilter(vtmp,sz,tau,inc)

  vtmpr=new Vector()
  vtmpr.copy(vtmp,vtmp.max_ind,vtmp.size-1)

  $o1 = new Vector()
  vtmpr.reverse()
  $o1.copy(vtmpr)
  vtmpr.reverse()
  $o1.append(vtmpr)

  $o1.div($o1.sum)
}

//compute sample correlation
//$o1 = vector # 1
//$o2 = vector # 2
//must have same size
func cor(){ local idx,m1,m2,sum
  if($o1.size!=$o2.size){
    printf("cross_cor err: input vecs not same size: %d & %d\n",$o1.size,$o2.size)
    return 0
  } else if($o1.size<2){
    printf("cross_cor err: input vecs must have size > 1\n")
    return 0
  }
  m1 = $o1.mean m2 = $o2.mean
  for(idx=0;idx<$o1.size;idx+=1){
    sum += ($o1.x(idx)-m1)*($o2.x(idx)-m2)
  }
  sum = sum / (($o1.size-1)*$o1.stdev*$o2.stdev)
  return sum
}

//plot clusters in 2 dimensions
//$o1 = list of vecs
//$o2 = classes
//$o3 = graph
//$4 = dim1
//$5 = dim2
func plot_kmeans(){ local idx
  if($4 >= $o1.count || $5 >= $o1.count){
     printf("idx out of bounds %d %d %d",$4,$5,$o1.count)
     return 0
  }
  if(mydbg){
    printf("$o1.o(%d).size=%d\n",$4,$o1.o($4).size)
    printf("$o1.o(%d).size=%d\n",$5,$o1.o($5).size)
  }
  if(mydbg){
    for(idx=0;idx<$o2.size;idx+=1){
      if(idx >= $o2.size){
        printf("idx out of $o2 bounds %d %d\n",idx,$o2.size)
        continue
      } else if(idx >= $o1.o($4).size || idx >= $o1.o($5).size){
        printf("idx out of $o1.o($4) or $o1.o($5) %d %d %d\n",idx,$o1.o($4).size,$o1.o($5).size)
        continue
      }
      $o3.mark($o1.o($4).x(idx),$o1.o($5).x(idx),"o",4,$o2.x(idx)+2,1)
    }
  } else {
    for(idx=0;idx<$o2.size;idx+=1){
      if(idx >= $o2.size){
        printf("idx out of $o2 bounds %d %d\n",idx,$o2.size)
        continue
      } else if(idx >= $o1.o($4).size || idx >= $o1.o($5).size){
        printf("idx out of $o1.o($4) or $o1.o($5) %d %d %d\n",idx,$o1.o($4).size,$o1.o($5).size)
        continue
      }
      $o3.mark($o1.o($4).x(idx),$o1.o($5).x(idx),"o",4,$o2.x(idx)+2,1)
    }
  }
  return 1
}


///////////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////
//The return value is
//0 for numbers, 
//1 for objref, 
//2 for strdef, and
//3 for pointers to numbers.
//const for return codes from argtype()
NUMBER_T  = 0
OBJREF_T  = 1
STRING_T  = 2
POINTER_T = 3
/////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////
/////////////////// fitness functions for ga /////////////////

//$o1 = spike train 1
//$o2 = spike train 2
func SpikeTrainFitness(){ local thresh,fitness,tmp,idx1,idx2,spiket1,spiket2,foundm,matchspikes
  ///////////////////////////////////////////////////////////////////
  //calculate fitness
  //
  //this will determine the "fitness", a measure of the similarity
  //between two spike trains

  thresh = 1.5         //max time between "matching" spikes
  fitness = 0
  
  tmp = abs($o1.size-$o2.size)
  fitness -= (tmp*tmp*tmp) //penalty for diff # of spikes

  //find # of matching spikes from set 1 to set 2
  matchspikes = 0  
  for(idx1 = 0; idx1 < $o1.size ; idx1 += 1){
    spiket1 = $o1.x(idx1)
    foundm = 0
    for(idx2 = 0; idx2 < $o2.size && !foundm; idx2 += 1){
      spiket2 = $o2.x(idx2)
      if(abs(spiket1-spiket2)<=thresh){
        matchspikes += 1
        foundm = 1
        break
      }
    }
  }

  //matching spikes add to fitness
  fitness += (matchspikes*matchspikes)
  tmp = $o1.size - matchspikes  
  //unmatching spikes sub from fitness
  fitness -= (tmp*tmp)
  
  //find # of matching spikes from set 2 to set 1
  matchspikes = 0  
  for(idx2 = 0; idx2 < $o2.size ; idx2 += 1){
    spiket2 = $o2.x(idx2)
    foundm = 0
    for(idx1 = 0; idx1 < $o1.size && !foundm; idx1 += 1){
      spiket1 = $o1.x(idx1)
      if(abs(spiket1-spiket2)<=thresh){
        matchspikes += 1
        foundm = 1
        break
      }
    }
  }  
  
  fitness += (matchspikes*matchspikes)
  tmp = $o2.size - matchspikes
  fitness -= (tmp*tmp)

  return fitness
  ///////////////////////////////////////////////////////////////////
}


//$o1 = spike train 1
//$o2 = spike train 2
//$3 = cost per unit time to move spike
maxspikes=512
double dtab[maxspikes][maxspikes]
func SpikeTrainEditDist(){ local cost,nspi,nspj,row,col
  nspi = $o1.size
  nspj = $o2.size
  cost = $3
  
  if(cost==0){
    return abs(nspi-nspj)
  } else if(cost >= 9e24) {
    return nspi+nspj
  }

  if(nspi==0){
    return nspj
  } else if(nspj==0){
    return nspi
  }
  
  for(row=0;row<nspi+1;row+=1){
    for(col=0;col<nspj+1;col+=1){
      dtab[row][col]=0
    }
  }

  for(row=0;row<nspi+1;row+=1){
    dtab[row][0] = row
  }

  for(col=0;col<nspj+1;col+=1){
    dtab[0][col] = col
  }

  for(row=1;row<nspi+1;row+=1){
    for(col=1;col<nspj+1;col+=1){
      dtab[row][col]=MIN(MIN(dtab[row-1][col]+1,dtab[row][col-1]+1),dtab[row-1][col-1]+cost*abs($o1.x(row-1)-$o2.x(col-1)))
    }
  }
  return dtab[nspi][nspj]   
}

//distance between params of two mainen cells
//$o1 = cell 1
//$o2 = cell 2
func MainenDist(){ local dist
  dist = abs($o1.rho-$o2.rho)+\
         abs($o1.kappa-$o2.kappa)+\
         abs($o1.dend.gmax_naz-$o2.dend.gmax_naz)+\
         abs($o1.dend.gmax_Nca-$o2.dend.gmax_Nca)+\
         abs($o1.dend.gmax_km-$o2.dend.gmax_km)+\
         abs($o1.dend.gmax_kca-$o2.dend.gmax_kca)
  return dist
}

//$o1 = spike train 1 - "guess" spike train - this tries to reach target
//$o2 = spike train 2 - target spike train - this is the end result we want
func SpikeTrainFitness2(){ local thresh,fitness,tmp,idx1,idx2,spiket1,spiket2,matchspikes

  //this function determines the "fitness" of how close
  //spike train $o1 is to spike train $o2

  thresh = 0.5         //max time between "matching" spikes

  if($o1.size==0){
    return -$o2.size
  }
  
  //find # of matching spikes from set 1 to set 2
  matchspikes = 0  
  for(idx1 = 0; idx1 < $o1.size ; idx1 += 1){
    spiket1 = $o1.x(idx1)
    for(idx2 = 0; idx2 < $o2.size ; idx2 += 1){
      spiket2 = $o2.x(idx2)
      if(abs(spiket1-spiket2)<=thresh){
        matchspikes += 1
        break
      }
    }
  }

  return (matchspikes/$o1.size) - abs($o1.size-$o2.size) / $o1.size
}

//returns min ISI distance btwn 2 sets of ISI trajectories
//$o1 = list of v1
//$o2 = list of v1'
//$o3 = list of v1''
//$o4 = list of v2
//$o5 = list of v2'
//$o6 = list of v2''
func ISIDistance(){ local idx1,idx2,dist,mindist,newsz localobj tmp1,tmp2,tmp3,tmp4,tmp5,tmp6
  dist = mindist = 9e308
  for(idx1=0;idx1<$o1.count;idx1+=1){
    for(idx2=0;idx2<$o4.count;idx2+=1){
      if($o1.o(idx1).size > $o4.o(idx2).size){ //resample $o4,5,6.o(idx) to size of $o1.o(idx)
        newsz = $o1.o(idx1).size
        tmp4 = new Vector(newsz) tmp4.samp($o4.o(idx2),newsz) 
        tmp5 = new Vector(newsz) tmp5.samp($o5.o(idx2),newsz)
        tmp6 = new Vector(newsz) tmp6.samp($o6.o(idx2),newsz)
        dist = VSQDist($o1.o(idx1),$o2.o(idx1),$o3.o(idx1),tmp4,tmp5,tmp6)
      } else if($o1.o(idx1).size < $o4.o(idx2).size){
        newsz = $o4.o(idx2).size
        tmp1 = new Vector(newsz) tmp1.samp($o1.o(idx1),newsz)
        tmp2 = new Vector(newsz) tmp2.samp($o2.o(idx1),newsz)
        tmp3 = new Vector(newsz) tmp3.samp($o3.o(idx1),newsz)
        dist = VSQDist(tmp1,tmp2,tmp3,$o4.o(idx2),$o5.o(idx2),$o6.o(idx2))
      } else{
        dist = VSQDist($o1.o(idx1),$o2.o(idx1),$o3.o(idx1),$o4.o(idx2),$o5.o(idx2),$o6.o(idx2))
      }
      if(dist != -1 && dist < mindist){
        mindist = dist
        if(mydbg){
          printf("idx1=%d idx2=%d mindist_so_far=%f\n",idx1,idx2,mindist)
        }
      }
    }
  }
  return mindist
}

//returns sum of min distance between each ISI trajectory in first set to 2nd set
//$o1 = list of v1
//$o2 = list of v1'
//$o3 = list of v1''
//$o4 = list of v2
//$o5 = list of v2'
//$o6 = list of v2''
//$&7 = best idx1 (optional)
//$&8 = best idx2 (optional)
func ISIDistance2(){ local idx1,idx2,dist,mindist,newsz,totaldist,savebestidx localobj tmp1,tmp2,tmp3,tmp4,tmp5,tmp6
  if(numarg()==8) savebestidx = 1 else savebestidx = 0 //
  totaldist = 0
  for(idx1=0;idx1<$o1.count;idx1+=1){
    dist = -1.0
    mindist = 9e25
    for(idx2=0;idx2<$o4.count;idx2+=1){
      if($o1.o(idx1).size > $o4.o(idx2).size){ //resample $o4,5,6.o(idx) to size of $o1.o(idx)
        newsz = $o1.o(idx1).size
        tmp4 = new Vector(newsz) tmp4.samp($o4.o(idx2),newsz) 
        tmp5 = new Vector(newsz) tmp5.samp($o5.o(idx2),newsz)
        tmp6 = new Vector(newsz) tmp6.samp($o6.o(idx2),newsz)
        dist = VSQDist($o1.o(idx1),$o2.o(idx1),$o3.o(idx1),tmp4,tmp5,tmp6)
      } else if($o1.o(idx1).size < $o4.o(idx2).size){
        newsz = $o4.o(idx2).size
        tmp1 = new Vector(newsz) tmp1.samp($o1.o(idx1),newsz)
        tmp2 = new Vector(newsz) tmp2.samp($o2.o(idx1),newsz)
        tmp3 = new Vector(newsz) tmp3.samp($o3.o(idx1),newsz)
        dist = VSQDist(tmp1,tmp2,tmp3,$o4.o(idx2),$o5.o(idx2),$o6.o(idx2))
      } else{
        dist = VSQDist($o1.o(idx1),$o2.o(idx1),$o3.o(idx1),$o4.o(idx2),$o5.o(idx2),$o6.o(idx2))
      }
      if(dist != -1.0 && dist < mindist){
        mindist = dist
        if(mydbg){
          printf("idx1=%d idx2=%d mindist_so_far=%f\n",idx1,idx2,mindist)
        }
        if(savebestidx){
          $&7 = idx1 $&8 = idx2
        }
      }
    }
    if(mindist < 9e25) totaldist += mindist
  }
  return totaldist
}

//////////////////////////////////////////////////////////////

///some templates...

//////////////////////////////////////////////////////////////////////////////////////
//template for storing inter-spike interval information
begintemplate isi_info

external GetISITimes,GetTrajectories,GetTrajectory
public start_times,end_times,lv1,lv2,lv3
objref start_times,end_times,lv1,lv2,lv3 //start,end isi times , v,v',v'' vectors
public setup

//$o1 = spike train
//$o2 = voltage vector
func setup(){ localobj spikes,voltvec
  if(numarg()!=2){
    return 0
  }
  spikes = $o1 voltvec = $o2
  if(!GetISITimes(spikes,start_times,end_times)) return 0
  lv1=new List() lv2=new List() lv3=new List()
  if(!GetTrajectories(voltvec,start_times,end_times,lv1,lv2,lv3)) return 0
  return 1
}

//$o1 = spike train
//$o2 = voltage vector
proc init(){

  start_times=new Vector()
  end_times=new Vector()
  lv1=new List()
  lv2=new List()
  lv3=new List()

  if(numarg()==0){
    return
  } else if(numarg()==2){
    setup($o1,$o2)
  }
}

endtemplate isi_info
//////////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////////
//template for storing spike-related info , including inter-spike interval info
begintemplate spike_info

external GetSpikeTimes,GetISITimes,GetTrajectories,GetTrajectory
public isi,spikes,setup
objref isi,spikes

//$o1 = voltage vector
//$2 = spike threshold
//$3 = min interspike time
//$4 = whether to initialize isi info
func setup(){ local sz
  if(numarg()!=4){
    printf("invalid args to spike_info.setup\n")
    return 0
  }
  spikes = new Vector($o1.size)
  sz = spikes.GetSpikeTimes($o1,$2,$3)
  spikes.resize(sz)
  if($4) return isi.setup(spikes,$o1)
  return 1
}

//either give 0 or 4 arguments
//$o1 = voltage vector
//$2 = spike threshold
//$3 = min interspike time
//$4 = whether to initialize isi info
proc init(){ local thresh,mintimediff,getisi localobj voltvec
  thresh = -20
  mintimediff = 2
  isi=new isi_info()
  spikes=new Vector()
  getisi = 1
  if(numarg()==0){
    return
  } else if(numarg()==4){
    voltvec = $o1
    thresh = $2
    mintimediff = $3
    getisi = $4
  }
  setup(voltvec,thresh,mintimediff,getisi)
}

endtemplate spike_info

///
//////////////////////////////////////////////////////////////////////////////////////

Loading data, please wait...