// $Id: nload.hoc,v 1.189 2011/02/16 15:46:37 samn Exp $ 

strdef strvecs,strnclnq //vecs file, vq nqs file, nqin nqs file, nqfld nqs file
strdef strv
strv="09apr17.18"
sprint(tstr,"o[%d]",numcols)
declare("nqtt",tstr,"nqte",tstr,"nqtea",tstr,"nqktau",tstr,"nclnq",tstr,"vLFP",tstr,"glfp","o[1]")
declare("vspudx",new Vector(),"vspudy",new Vector())
declare("nqCO",tstr,"nqisi",tstr,"nqCTY",tstr)
if(name_declared("snq")) {
  nqsdel(snq)
  objref snq[numcols]
} else declare("snq",tstr)
if(name_declared("fnq")) {
  nqsdel(fnq)
  objref fnq[numcols]
} else declare("fnq",tstr)
if(name_declared("vq")) {
  nqsdel(vq)
  objref vq[numcols]
} else declare("vq",tstr)
declare("CDX",0)
sprint(tstr,"o[%d][10]",numcols)
declare("nqIN",tstr)
declare("skipsnq",1) // iff==1 dont get snq and isinq
declare("binsz",100) //100ms
{pflg=1 ers=0}
declare("htmin",0)//sgrdel //min time in histogram
declare("htmax",mytstop)  //max time in histogram
for i=0,numcols-1 vq[i]=lcol.o(i).cstim.vq

//* VQIsiNQS(vq,sy) - get an NQS with ISI from external inputs stored in vq
// only uses synapses type specified by sy
obfunc VQIsiNQS() { local idx,tt,jdx,sy,a localobj nqisi,vid,vd,vt,vq
  vq=$o1 sy=$2 a=allocvecs(vid,vd)
  {nqisi=new NQS("ind","visi") nqisi.odec("visi")}
  {vq.verbose=nqisi.verbose=0 vq.tog("DB")}
  vid.resize(vq.v.size) vq.getcol("ind").uniq(vid)
  for vtr(&idx,vid) if(vq.select("ind",idx,"sy",sy)>2) {
    vt=vq.getcol("time")
    vd.resize(0) vd.deriv(vt,1,1)
    nqisi.append(idx,vd)
  }
  {vq.tog("DB") vq.verbose=nqisi.verbose=1 dealloc(a)}
  return nqisi
}

//* getLV(visi) -- get LV (local variation)
func getLV () { local i,j,s,n,d localobj vec
  vec=$o1 s=0
  if(vec.size < 2) {printf("getLV ERRA: input vec too small!\n") return 0}
  for i=0,vec.size-2 {
    n = (vec.x(i)-vec.x(i+1))^2
    d = (vec.x(i)+vec.x(i+1))^2
    s += 3*n/d
  }
  return s / (vec.size-1)
}

//* addCVcol(nqisi) - add column to input NQS with CV -- input NQS has to have a visi objref col
// containing a Vector of interspike intervals
proc addCVcol () {  local i,c,m,a localobj nqi,vec
  a=allocvecs(vec) nqi=$o1 nqi.tog("DB")
  if((c=nqi.fi("cv")) == -1){ nqi.resize("cv") nqi.pad() c=nqi.m-1}
  for i=0,nqi.v.size-1 {
    vec.resize(0) vec.copy(nqi.get("visi",i).o)
    m=vec.mean()
    if(m) nqi.v[nqi.m-1].x(i) = vec.stdev / m else nqi.v[nqi.m-1].x(i)=0
  }
  dealloc(a)
}

//* addLVcol(nqisi) - add column to input NQS with LV (local variability of interspike interval) --
// input NQS has to have a visi objref col containing a Vector of interspike intervals
// LV is 1 for poisson ISI, 0 for regular
proc addLVcol () {  local i,c,m,a localobj nqi,vec
  a=allocvecs(vec) nqi=$o1 nqi.tog("DB")
  if((c=nqi.fi("lv")) == -1){ nqi.resize("lv") nqi.pad() c=nqi.m-1}
  for i=0,nqi.v.size-1 {
    {vec.resize(0) vec.copy(nqi.get("visi",i).o)}
    nqi.v[nqi.m-1].x(i) = vec.lvar()
  }
  dealloc(a)
}

//* IsiNQS() - makes an NQS with ISI,CV,LV information for each cell, uses globals: CDX,col,snq!!!
obfunc IsiNQS () { local idx,tt,jdx,a localobj nqisi,vid,vd,vt,vty,vic
  if(snq[CDX]==nil)snq[CDX]=SpikeNQS(vit[CDX])
  {nqisi=new NQS("gid","id","type","ice","visi") nqisi.odec("visi")}
  snq[CDX].verbose=nqisi.verbose=0
  a=allocvecs(vid,vd)
  for idx=0,col[CDX].allcells-1 {
    if(snq[CDX].select("id",idx)>2) {
      vt=snq[CDX].getcol("t")
      vty=snq[CDX].getcol("type")
      vic=snq[CDX].getcol("ice")
      vd.resize(0) vd.deriv(vt,1,1) // ISI
      nqisi.append(col[CDX].ce.o(idx).gid,idx,vty.x(0),vic.x(0),vd)
    }
  }
  addCVcol(nqisi) // coefficient of variation
  addLVcol(nqisi) // local variation
  snq[CDX].tog("DB")
  dealloc(a)
  snq[CDX].verbose=nqisi.verbose=1
  return nqisi
}

//* SpikeNQS (time vec, id vec, cell list) returns NQS containing gid,id,type,ice,t<-spike times
obfunc SpikeNQS (){ local idx,gcol,tycol,icecol,idcol,tcol,colcol localobj vec,tvec,nq,ce,c
  {tvec=$o1 vec=$o2 ce=$o3.ce}
  nq=new NQS("col","gid","id","type","ice","t") 
  {gcol=nq.fi("gid") tycol=nq.fi("type") icecol=nq.fi("ice") idcol=nq.fi("id") tcol=nq.fi("t") colcol=nq.fi("col")}
  {nq.v[idcol].copy(vec) nq.v[tcol].copy(tvec) nq.pad()}
  for idx=0,vec.size-1 { c=ce.o(vec.x(idx))
    nq.v[colcol].x(idx)=c.col
    nq.v[gcol].x(idx)=c.gid()
    nq.v[tycol].x(idx)=c.type
  }
  {nq.v[icecol].copy(nq.v[tycol]) nq.v[icecol].apply("ice")}
  return nq
}

//* FreqNQS(snq,ce,[,duration])
//get an NQS with rate of firing of each cell
//$o1 = spikenqs, $o2 = cell list, $3 = duration (starts from 0)
obfunc FreqNQS () { local i,a,ns,tt localobj snq,fnq,vi,ce
  a=allocvecs(vi) fnq=new NQS("col","gid","id","type","ice","freq") snq=$o1 ce=$o2.ce
  if(numarg()>2) tt=$3 else tt=tstop
  snq.tog("DB") vi.indgen(0,$o2.allcells-1,1)
  for vtr(&i,vi) fnq.append(ce.o(i).col,ce.o(i).gid,i,ce.o(i).type,ice(ce.o(i).type),1e3*snq.select(-1,"id",i)/tt)
  dealloc(a) return fnq
}

//* clrMyNQs - delete nqCO[CDX], nqIN[CDX]
proc clrMyNQs () { local i
  nqsdel(nqCO[CDX])
  for i=0,9 nqsdel(nqIN[CDX][i])
  nqsdel(snq[CDX])
  nqsdel(nqisi[CDX])
  nqsdel(nqCTY[CDX])
}

//* rrounddI
func rrounddI () {
  return int($1+0.5)
}

//** getnqlfpcor(nqs with columns as LFP, window size in ms)
obfunc getnqlfpcor () { local winsz,tt,et,maxt,maxi,i,j,a localobj v0,v1,nqlfp,nqo
  nqlfp=$o1 a=allocvecs(v0,v1) winsz=$2
  {maxt=nqlfp.v.size*vdt_INTF6-1e3 maxi=nqlfp.v.size-1 nqo=new NQS("c0","c1","tt","r")}
  {nqo.v.resize((maxt/winsz)*nqlfp.m*(nqlfp.m-1)/2) nqo.pad() nqo.clear()}
  for(tt=1e3;tt+winsz<=maxt;tt+=winsz) {
    {et=MINxy((tt+winsz)/vdt_INTF6,maxi) vrsz(0,v0,v1)}
    for i=0,nqlfp.m-1 { v0.copy(nqlfp.v[i],tt/vdt_INTF6,et)
      for j=i+1,nqlfp.m-1 { v1.copy(nqlfp.v[j],tt/vdt_INTF6,et)
        nqo.append(i,j,tt,v0.pcorrel(v1))
      }
    }
  }
  {dealloc(a) return nqo}
}

//** getlfppco(nqs from getnqlfpcor)
obfunc getlfppco () { local t1,t2,csz,rc,i,j,a localobj mc,v0,v1,nqc
  {a=allocvecs(v0,v1) nqc=$o1 nqc.tog("DB") nqc.verbose=0}
  {v0.resize(nqc.v.size) nqc.getcol("c0").uniq(v0)}
  {csz=(1+v0.size)*(v0.size)/2 rc=nqc.fi("r") vrsz(0,v0,v1)}
  {mc=new Matrix(nqc.v.size/csz,nqc.v.size/csz) i=j=0}
  for(t1=0;t1<nqc.v.size;t1+=csz) { v0.copy(nqc.v[rc],t1,t1+csz-1)
    {mc.x(i,i)=1 j=i+1}
    for(t2=t1+csz;t2<nqc.v.size;t2+=csz) { v1.copy(nqc.v[rc],t2,t2+csz-1)
      mc.x(i,j)=mc.x(j,i)=v0.pcorrel(v1)
      j+=1
    }
    i+=1
  }
  {nqc.verbose=1 dealloc(a) return mc}
}

//** gotolfp(nqlfp,nqspud,tmin,tmax,offy)
func gotolfp () { local a,b,c,i,offy,tmin,tmax localobj nqlfp,nqspud,vtmpx,vtmpy,vsub
  if(glfp==nil)glfp=new Graph()
  a=allocvecs(vtmpx,vtmpy,vsub)
  {nqlfp=$o1 nqspud=$o2 tmin=$3 tmax=$4 offy=$5}
  if(vLFP==nil) {
    for i=0,nqlfp.m-1 {
      vLFP[i]=new Vector()
      vLFP[i].plot(glfp,vdt_INTF6)
    }
  }
  glfp.erase
  vsub.resize(nqlfp.m)
  for i=0,nqlfp.m-1 { // draw lfps
    {vLFP[i].resize(0) vLFP[i].copy(nqlfp.v[i],tmin/vdt_INTF6,tmax/vdt_INTF6)}
    {vsub.x(i)=vLFP[i].min vLFP[i].add(offy*i-vsub.x(i)) vLFP[i].line(glfp,vdt_INTF6)}
  }
  if(nqspud!=nil) { // draw peak of bumps
    {nqspud.verbose=0 vrsz(0,vspudx,vspudy)}
    for i=0,nqlfp.m-1 if(nqspud.select("LOC","[]",tmin,tmax,"COL",i)) {
      vrsz(0,vtmpx,vtmpy)
      {vtmpy.copy(nqspud.getcol("PEAK")) vtmpy.add(i*offy-nqlfp.v[i].min)}
      {vtmpx.copy(nqspud.getcol("LOC")) vtmpx.sub(tmin)}
      {vspudx.append(vtmpx) vspudy.append(vtmpy)}
    }
    vspudy.mark(glfp,vspudx,"O",4,2,1)
    {nqspud.verbose=1 dealloc(a)}
  }
  {dealloc(a) glfp.exec_menu("View = plot")}
  return 1
}

//** make a single NQS out of nqLFP, using last full LFP in each
obfunc getnqlfp () { local i localobj nqlfp 
  nqlfp=new NQS(numcols)
  for i=0,numcols-1 nqlfp.v[i].copy(nqLFP[i].v[nqLFP[i].m-1])
  return nqlfp
}

//** getnqspud(nqs with columns as LFP)
obfunc getnqspud () { local i,off,x,a localobj vec,nqu,nqin,nqtmp
  {a=allocvecs(vec) test_do_graphs=0 nqin=$o1 printStep=vdt_INTF6}
  if(numarg()>1) off=$2 else off=1e3
  for i=0,nqin.m-1 {
    vec.resize(0) vec.copy(nqin.v[i],off/vdt_INTF6,(tstop-off)/vdt_INTF6)
    nqtmp=testupdown(vec,10,0)
    {nqtmp.resize("COL") nqtmp.pad() nqtmp.v[nqtmp.m-1].fill(i)}
    if(nqu==nil) {
      {nqu=new NQS() nqu.cp(nqtmp)}
    } else nqu.append(nqtmp)
    nqsdel(nqtmp)
  }  
  for case(&x,0,5) nqu.v[x].add(off) // add offset so have sim times
  {dealloc(a) return nqu}
}

//* initnqCTY(col id) - makes an NQS with summed spike-train activity by type, requires nqCO to be setup already
// nqCO[cid].v[celltype] is summed activity of the given celltype
proc initnqCTY () { local i,cid,ty
  {cid=$1 nqsdel(nqCTY[cid]) nqCTY[cid]=new NQS(CTYPi)}
  for i=0,nqCO[cid].m-1 {
    ty=col[cid].ce.o(i).type
    if(nqCTY[cid].v[ty].size==0) nqCTY[cid].v[ty].resize(nqCO[cid].v[i].size)
    nqCTY[cid].v[ty].add(nqCO[cid].v[i])
  }
}

//* getnqCTYnt - returns an nqs of nTE between summed counts per type (nqCTY[...]), across/within columns
obfunc getnqCTYnt () { local i,c1,c2,nt,nshuf,from,to,a localobj nq,vtmp
  a=allocvecs(vtmp) vtmp.resize(35) nshuf=30
  {nq=new NQS("col1","col2","fs","ts","from","to","ic1","ic2","nte","samec") nq.strdec("fs") nq.strdec("ts")}
  for c1=0,numcols-1 for c2=0,numcols-1 {
    for from=0,CTYPi-1 for to=0,CTYPi-1 {
      if(from==to && c1==c2) continue
      if(nqCTY[c1].v[from].size<1 || nqCTY[c2].v[to].size<1) continue
      nt = nqCTY[c1].v[from].tentropspks(nqCTY[c2].v[to],vtmp,nshuf)
      nq.append(c1,c2,CTYP.o(from).s,CTYP.o(to).s,from,to,ice(from),ice(to),nt,c1==c2)
    }
  }
  {dealloc(a) return nq}
}

//* getnTEmap(nqs from getnqCTYnt, col1, col2, type1, type2, map of types to matrix coordinates)
// returns a matrix where element x,y is nTE from row y -> column x
obfunc getnTEmap () { local ty1,ty2,c1,c2,cid,from,to localobj mc,vty1,vty2,vm,nq
  nq=$o1 c1=$2 c2=$3 vty1=$o4 vty2=$o5 vm=$o6 nq.verbose=0
  mc=new Matrix(vty1.size,vty2.size)
  for vtr(&from,vty1) for vtr(&to,vty2) if(nq.select("col1",c1,"col2",c2,"from",from,"to",to)) {    
    mc.x(vm.x(from),vm.x(to)) = MAXxy(nq.fetch("nte"),0)
  } else mc.x(vm.x(from),vm.x(to))=0
  nq.verbose=1
  return mc
}

//* colnTEmap(nq from getnqCTYnt, list for output matrices,inter flag)
// gets COLUMN nTE map, uses all cell types and columns
// returns average of list for output matrices
// iff inter == 0, then gets INTRA-column map, iff inter == 1, gets INTER-column map
obfunc colnTEmap () { local i,j,cnt,inter,a localobj mc,vty1,vty2,vm,lc,nq
  {a=allocvecs(vty1,vty2,vm) vrsz(0,vty1,vty2)}
  nq=$o1 lc=$o2 if(lc==nil){lc=new List()} inter=$3
  {vty1.append(E2,I2,I2L,E4,I4,I4L,E5R,E5B,I5,I5L,E6,I6,I6L) vty2.copy(vty1)}
  {vm.resize(CTYPi) vm.fill(0) cnt=0}
  for vtr(&i,vty1) {vm.x(i)=cnt cnt+=1} // map cell types to index in Matrix
  if(inter) {
    for i=0,numcols-1 for j=0,numcols-1 if(i!=j) lc.append(getnTEmap(nq,i,j,vty1,vty2,vm))
  } else for i=0,numcols-1 lc.append(getnTEmap(nq,i,i,vty1,vty2,vm))
  mc=avgmat(lc) // average Matrix
  {dealloc(a) return mc}
}

//* initMyNQs - init binned value NQSs using raster info from vit[CDX] and input info from vq
// nqCO[CDX] = nqs with counts/time for each cell
// nqIN[CDX][sy] = nqs with counts/time for external inputs to synapse type sy @ all the cells
proc initMyNQs () { local ct,idx,sy,a localobj vidx,vt
  a=allocvecs(vidx,vt) vq[CDX].verbose=0
  if(nqCO[CDX]==nil) nqCO[CDX]=new NQS(col[CDX].allcells) else nqCO[CDX].clear() //NQS for the binned spike-trains
  vit[CDX].vec.ihist(vit[CDX].tvec,nqCO[CDX].vl,htmin,htmax,binsz)
  for case(&sy,AM,NM,GA,AM2,NM2,GA2) if(vq[CDX].select("sy",sy)) {
    flag_stats=1 //allow vq[CDX].v[1] to be unsorted  
    if(nqIN[CDX][sy]==nil) nqIN[CDX][sy]=new NQS(vq[CDX].v[0].max+1) else nqIN[CDX][sy].clear() //NQS for the binned input events
    {vidx.copy(vq[CDX].getcol("ind")) vt.copy(vq[CDX].getcol("time"))}
    vidx.ihist(vt,nqIN[CDX][sy].vl,htmin,htmax,binsz)  
  }
  if(!skipsnq) {
    {nqsdel(snq[CDX]) snq[CDX]=SpikeNQS(vit[CDX].tvec,vit[CDX].vec,lcol.o(CDX))}
    {nqsdel(fnq[CDX]) fnq[CDX]=FreqNQS(snq[CDX],lcol.o(CDX))}
  }
  initnqCTY(CDX)
  {vq[CDX].verbose=1 dealloc(a)}
}

//* initAllMyNQs() - initMyNQs for all numcols
proc initAllMyNQs () { local cdx
  cdx=CDX
  for CDX=0,numcols-1 initMyNQs()
  CDX=cdx
}

//* clrAllMyNQs() - clrMyNQs for all numcols
proc clrAllMyNQs () { local cdx
  cdx=CDX
  for CDX=0,numcols-1 clrMyNQs()
  CDX=cdx
}
 
//* getTEKTaunq() gets transfer entropy and ktau from inputs->cells
obfunc getTEKTaunq () { local sy,loc,idx,vbt localobj nqt,vtmp
  {nqt=new NQS("id","sy","loc","te","nte","ktau") initMyNQs()}
  {vbt=verbose_infot verbose_infot=0}
  for case(&sy,AM,NM,GA,AM2,NM2,GA2) if(nqIN[CDX][sy]!=nil) {
    if(sy<AM2)loc=SOMA else loc=DEND
    for idx=0,col[CDX].allcells-1 { vtmp=normte(nqIN[CDX][sy].v[idx],nqCO[CDX].v[idx],30)
      if(idx%100==0) printf("%s:idx=%d\n",STYP.o(sy).s,idx)
      nqt.append(idx,sy,loc,vtmp.x(0),vtmp.x(2),nqIN[CDX][sy].v[idx].kcorrel(nqCO[CDX].v[idx],1))
    }
  }
  for idx=nqt.m-2,nqt.m-1 nqt.v[idx].unnan() //get rid of nans
  {verbose_infot=vbt return nqt}
}
 
// //* getTEnqAllsub -- 
// // gets nTE from all inputs of a cell type -> all outputs of a specific cell type
// obfunc getTEnqAllsub () { local a,row,ct,sy,id,poid,pls,prid,cts,tyt,\
//                        idcol,sycol,INTFidcol,pulsecol,loccol\
//                        localobj nqt,vtmp
// cts=$1
// nqin.tog("DB")  nqt=new NQS("id","type","prid","poid","pulse","sy","same","loc","te","nte")
// initMyNQs()
// {idcol=nqin.fi("id")  sycol=nqin.fi("sy")  INTFidcol=nqin.fi("INTFid")}
// {pulsecol=nqin.fi("pulse")  loccol=nqin.fi("loc")}
// {nqin.select(-1,"type",cts)  nqin.verbose=0 pls=0}
// printf("\nid:")
// for id=ix[cts],ixe[cts] {
//   ct = ce.o(id).type
//   if(id%20==0) printf("\n")
//   printf("%d ",id)
//   for vtr(&row,nqin.ind) {
//     {prid=nqin.v[INTFidcol].x(row)      poid=nqin.v[idcol].x(row)}
//     {sy=nqin.v[sycol].x(row)      loc=nqin.v[loccol].x(row)}
//     vtmp = normte(nqIN[CDX].v[prid],nqCO[CDX].v[id],30)
//     nqt.append(id,ct,prid,poid,pls,sy,poid==id,loc,vtmp.x(0),vtmp.x(2))
//   }
// }
// printf("\n")
// for id=0,nqt.m-1 nqt.v[id].unnan()
// nqin.verbose=1
// return nqt
// }
// 
// 
// //* ccntevsbinsz(prid,poid,vbins,nq)
// obfunc ccntevsbinsz () { local prid,poid localobj vbins,vtmp,nq
// prid=$1 poid=$2 vbins=$o3
// nq=$o4
// if(nq==nil) nq=new NQS("prid","poid","binsz","te","nte","prty","poty") else nq.tog("DB")
// for vtr(&binsz,vbins) {
//   print "binsz" , binsz
//   clrMyNQs()
//   initMyNQs()
//   vtmp = normte(nqCO[CDX].v[prid],nqCO[CDX].v[poid],30)
//   nq.append(prid,poid,binsz,vtmp.x(0),vtmp.x(2),ce.o(prid).type,ce.o(poid).type)
// }
// return nq
// }
// //* tyntevsbinsz(Vec of pre-type, Vec of po-type, Vec of bin sizes, nq output)
// obfunc tyntevsbinsz () { local prid,poid,prty,poty,verb localobj vbins,vtmp,nq,vprty,vpoty
// vprty=$o1 vpoty=$o2 vbins=$o3 nq=$o4
// if(numarg()>4)verb=$5 else verb=0
// if(nq==nil) nq=new NQS("prid","poid","binsz","te","nte","prty","poty") else nq.tog("DB")
// nclnq[CDX].verbose=0
// for vtr(&binsz,vbins) {
//   print "binsz " , binsz
//   clrMyNQs()
//   initMyNQs()
//   for vtr(&prty,vprty) for vtr(&poty,vpoty) {
//     if(prty==poty || div[prty][poty][0]<1) continue
//     print "prty : " , CTYP.o(prty).s, " , poty : ", CTYP.o(poty).s
//     for prid=ix[prty],ixe[prty] {
//       if(nclnq[CDX].select("prid",prid,"poty",poty)) for nclnq[CDX].qt(&poid,"poid") {
//         if(verb) print "prid " , prid , " poid " , poid
//         vtmp=normte(nqCO[CDX].v[prid],nqCO[CDX].v[poid],30)
//         nq.append(prid,poid,binsz,vtmp.x(0),vtmp.x(2),prty,poty)
//       }
//     }
//   }
// }
// nclnq[CDX].verbose=1
// return nq
// }

//* binarize(vec) - set all values > 0 to 1, everything else to 0
proc binarize () { local i,sz localobj vin
  vin=$o1 sz=vin.size()-1
  for i=0,sz if(vin.x(i)>0) vin.x(i)=1 else vin.x(i)=0
}

//* mknclnq - make an NQS with connectivity info
proc mknclnq () { local idx,jdx,ty1,ty2,sy,loc,ic1,ic2,lt1,gid1,gid2,c1,c2,c localobj adj,ce,nc
  if(nclnq!=nil) return
  nclnq=new NQS("prid","poid","prty","poty","sy","loc","col1","col2","gid1","gid2","w")
  for c=0,numcols-1 { // intra-column wiring for all the columns
    {ce=col[c].ce adj=AdjList(0,ce.count-1,0,ce)}
    for idx=0,adj.count-1 { ty1=ce.o(idx).type
      gid1=ce.o(idx).gid
      if(!ice(ty1)) {
        for vtr(&jdx,adj.o(idx)) { ty2=ce.o(jdx).type
          nclnq.append(idx,jdx,ty1,ty2,AM2,synloc[ty1][ty2],c,c,gid1,ce.o(jdx).gid,wmat[ty1][ty2][AM2][0])
          nclnq.append(idx,jdx,ty1,ty2,NM2,synloc[ty1][ty2],c,c,gid1,ce.o(jdx).gid,wmat[ty1][ty2][NM2][0])
        }
      } else {
        if(IsLTS(ty1)) sy=GA2 else sy=GA
        for vtr(&jdx,adj.o(idx)) { ty2=ce.o(jdx).type
          nclnq.append(idx,jdx,ty1,ty2,sy,synloc[ty1][ty2],c,c,gid1,ce.o(jdx).gid,wmat[ty1][ty2][sy][0])
        }
      }
    }
  }
  for ltr(nc,ncl) { // inter-column wiring for all the columns
    if(!(isojt(nc.pre,INTF6[0]) && isojt(nc.syn,INTF6[0]))) continue // only looking btwn INTF6 cells
    {idx=nc.pre.id ty1=nc.pre.type gid1=nc.pre.gid c1=nc.pre.col}
    {jdx=nc.syn.id ty2=nc.syn.type gid2=nc.syn.gid c2=nc.syn.col}
    if(!ice(ty1)) {
      nclnq.append(idx,jdx,ty1,ty2,AM2,synloc[ty1][ty2],c1,c2,gid1,gid2,nc.weight(AM2))
      nclnq.append(idx,jdx,ty1,ty2,NM2,synloc[ty1][ty2],c1,c2,gid1,gid2,nc.weight(NM2))
    } else {
      if(IsLTS(ty1)) sy=GA2 else sy=GA
      nclnq.append(idx,jdx,ty1,ty2,sy,synloc[ty1][ty2],c1,c2,gid1,gid2,nc.weight(sy))
    }
  }
}

//* FullAdjList(makenclnq,skipI) - make adjacency list that has intercol/intracol info
// uses global identifier (gid) of each cell (not intracol id)
obfunc FullAdjList () { local i,j,sz,skipI,mknq,a localobj adj,vt,vg1,vg2,vex
  if(nclnq==nil)mknq=1 else {
    if(numarg()>0)mknq=$1 else mknq=1
  }
  if(mknq) {nqsdel(nclnq) mknclnq()} // has all the connectivity info
  nclnq.tog("DB")
  if(numarg()>1)skipI=$2 else skipI=1 //skip inhib cells by default
  a=allocvecs(vex,vt)
  {adj=new List() vg1=nclnq.getcol("gid1") vg2=nclnq.getcol("gid2")}
  sz=MAXxy(vg1.max,vg2.max)
  for i=0,sz adj.append(new Vector())
  if(skipI) { vex.resize(sz+1) // skip inhib cells
    for i=0,sz vex.x(i)=(INTF6[i].flag("inh")==0)
    for i=0,vg1.size-1 if(vex.x(vg1.x(i)) && vex.x(vg2.x(i))) {
      adj.o(vg1.x(i)).append(vg2.x(i))
    }
  } else for i=0,vg1.size-1 adj.o(vg1.x(i)).append(vg2.x(i))
  for i=0,sz if(adj.o(i).size) { // make sure all are unique, since nclnq has connection info for each synapse type
    {vt.resize(adj.o(i).size) adj.o(i).uniq(vt)}
    {adj.o(i).resize(0) vt.sort() adj.o(i).copy(vt)}
  }
  dealloc(a)
  return adj
}

//* getsy(from,to) - get synapse type , either returns AM2 or GA,GA2 , never returns NMDA
func getsy () { local from,to
  from=$1 to=$2
  if(ice(from)) {
    if(IsLTS(from)) return GA2 else return GA
  } else return AM2
}
 
//* getvpow(nqf,minf,maxf)
//get the power in a range of frequencies from nqf, the matspecgram spectrogram nqs from nqfld
obfunc getvpow () { local a,minf,maxf,i localobj vtmp,vout,nqf,vf
  nqf=$o1 minf=$2 maxf=$3
  a=allocvecs(vtmp,vf) vout=new Vector()
  {vtmp=nqf.get("pow",0).o print vtmp vout.resize(vtmp.size) vout.fill(0)}
  {vf.resize(nqf.v.size)  nqf.getcol("f").uniq(vf)}
  for i=0,vf.size-1 {
    if(vf.x(i)>=minf && vf.x(i)<=maxf) {
      vtmp=nqf.get("pow",i).o
      vout.add(vtmp)
    }
    if(vf.x(i)>maxf) break
  }
  dealloc(a)
  return vout
}

//* getpfcormat(nq1[,nq2]) -- nq1,nq2 are from matspecgram, getpfcormat gets power fluctuation correlation matrix
// can get nq1,nq2 from nqfld+matspecgram like this: nq1=matspecgram(nqLFP.v,sampr,100,0,1)
// nq2 is optional
obfunc getpfcormat () { local r1,r2,sz,a localobj mc,v1,v2,nq1,nq2
  a=allocvecs(v1,v2) 
  if(numarg()==1){nq1=nq2=$o1} else {nq1=$o1 nq2=$o2}
  {sz=MAXxy(nq1.v.size,nq2.v.size) mc=new Matrix(sz,sz)}
  if(nq1==nq2) {
    for r1=0,nq1.v.size-1 { v1.copy(nq1.get("pow",r1).o)
      mc.x(r1,r1)=1
      for r2=r1+1,nq2.v.size-1 { v2.copy(nq2.get("pow",r2).o)
        mc.x(r1,r2) = mc.x(r2,r1) = v1.pcorrel(v2)
      }
    }
  } else {
    for r1=0,nq1.v.size-1 { v1.copy(nq1.get("pow",r1).o)
      for r2=0,nq2.v.size-1 { v2.copy(nq2.get("pow",r2).o)
        mc.x(r1,r2) = v1.pcorrel(v2)
      }
    }
  }
  dealloc(a)
  return mc
}

//* getsyncarea(nqf,mc,minf,maxf,[th]) -- nqf is from matspecgram, mc is from getpfcormat
// gets the 'synchronization area' -- the ratio of pixels in matrix with r > th / that area of the matrix
// frequencies/area used are in range [minf,maxf]
func getsyncarea () { local r1,c1,sa,cnt,minf,maxf,th localobj nqf,mc
  nqf=$o1 mc=$o2 minf=$3 maxf=$4 sa=cnt=0
  if(numarg()>4)th=$5 else th=0.2
  for r1=0,mc.nrow-1 if(nqf.v[0].x(r1)>=minf && nqf.v[0].x(r1)<=maxf) {
    for c1=0,r1-1 if(nqf.v[0].x(c1)>=minf && nqf.v[0].x(c1)<=maxf) {
      if(mc.x(r1,c1)>th) sa+=1
      cnt+=1
    }
  }
  if(cnt>0) return sa/cnt else return 0
}
 
//* FileExists(path) - 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
}

//* loadmydat() - load in data from sim - this func depends on strv being set
// $s1 == vecs file, [$s2 == nclnq file]
func loadmydat () { local cdx,idx,a localobj co,vi,str,vtmp,svq
  {a=allocvecs(vtmp) strvecs=$s1 str=new String2() svq=new String()}
  clrAllMyNQs() // free histogram NQSs
  printf("reading in %s vecs file\n",strvecs)
  gvnew(strvecs) // load the Vectors
  for ltr(co,lcol,&cdx) {
    {vit[cdx].vec.resize(0) vit[cdx].tvec.resize(0)}
    {sprint(str.s,"%s_SPKS",co.name) sprint(str.t,"%s_LFP",co.name)}
    for ltr(vi,panobj.printlist,&idx) {
      if(!strcmp(vi.name,str.s)) {
        print "reading in spikes for ",co.name, " from vitem ", vi.name
        panobj.rv_readvec(idx,vit[cdx].tvec,vit[cdx].vec) //read in raster
      }
    }
    sprint(svq.s,"data/%s_%s_vq.nqs",strv,co.name)
    if(FileExists(svq.s)) {
      print "reading in vq for ", co.name, " from ", svq.s
      vq[cdx].rd(svq.s)
    } else print "loadmydat WARNING: file ", svq.s , " doesn't exist."
    sprint(str.s,"data/%s_%s_LFP.nqs",strv,col[cdx].name)
    if(FileExists(str.s)) {
      {nqsdel(nqLFP[cdx]) nqLFP[cdx]=new NQS(str.s)}
    } else print "loadmydat WARNING: file ", str.s , " doesn't exist."
  }
  if(numarg()>1) {
    {strnclnq=$s2 printf("reading in %s nclnq nqs file\n",strnclnq)}
    if(nclnq==nil) nclnq=new NQS(strnclnq) else nclnq.rd(strnclnq)
  }
  {dealloc(a) return 1}
}
 
//* pravgrates([duration,allcolumns,strdef output]) - print average rates of cell types after run
//print average firing rates from vit[...], assumes sgrdur,numc same as when ran sim
//when allcolumns==0 (default) uses vit[CDX] (CDX is global index for COLUMN)
proc pravgrates () { local idx,ct,dur,cdx,sidx,eidx,all,a localobj vspk,vnc,str
  if(numarg()>0)dur=$1 else dur=tstop
  if(numarg()>1)all=$2 else all=0   
  {a=allocvecs(vspk,vnc) vrsz(CTYPi+1,vspk,vnc) idx=ct=0 sidx=eidx=CDX str=new String2()}
  if(all){sidx=0 eidx=numcols-1}
  for cdx=sidx,eidx {
    for vtr(&idx,vit[cdx].vec) vspk.x(col[cdx].ce.o(idx).type)+=1
    for ct=0,CTYPi-1 vnc.x(ct)+=col[cdx].numc[ct]
  }
  for ct=0,CTYPi-1 if(vnc.x(ct)) {
    sprint(str.t,"avg rate for %s: %gHz\n",CTYP.o(ct).s,(1e3*vspk.x(ct)/dur)/vnc.x(ct))
    strcat(str.s,str.t)
  }
  print str.s
  if(numarg()>2) $s3=str.s
}
//* prctcellsfired - get % of cells that fired after a run
func prctcellsfired () { local a,i,pr,ct localobj vf
  a=allocvecs(vf) 
  if(numarg()>0) {
    ct=$1
    if(ct<0 || ct>=CTYPi) {
      print "invalid ct = ", ct
      return 0
    }
    vrsz(numc[ct],vf) vf.fill(0)
    for vtr(&i,vit[CDX].vec) if(col[CDX].ce.o(i).type==ct) vf.x(i-ix[ct])=1
    pr=vf.count(1)/numc[ct]
  } else {
    vrsz(col[CDX].allcells,vf) vf.fill(0)
    for vtr(&i,vit[CDX].vec) vf.x(i)=1  
    pr=vf.count(1)/col[CDX].allcells
  }
  dealloc(a)
  return pr
}
 
//* minrunsv() minimal run/save of data , just save vecs, vq, nclnq
// $1 - optional , specifices whether to save nclnq
proc minrunsv () { local svc,cdx,svq localobj st,co
  {st=new String() time("run()")}
  printf("\nran sim for %g\n",tstop)
  sprint(st.s,"data/%s_.vecs",strv)
  {panobj.pvplist(st.s,strv) printf("saved vecs\n")}
  for cdx=0,numcols-1 { // save NQS with LFPs
    sprint(tstr,"data/%s_%s_LFP.nqs",strv,col[cdx].name)
    {nqLFP[cdx].tog("DB") nqLFP[cdx].sv(tstr)}
  }
  if(numarg()>0) svc=$1 else svc=1
  if(numarg()>1) svq=$2 else svq=1
  if(svc) {
    sprint(st.s,"data/%s_nclnq.nqs",strv)
    {mknclnq() nclnq.tog("DB") nclnq.sv(st.s)}
    printf("saved connectivity\n")
  }
  if(svq) for ltr(co,lcol,&cdx) {
    sprint(st.s,"data/%s_%s_vq.nqs",strv,co.name)
    {vq[cdx].tog("DB") vq[cdx].sv(st.s) printf("saved inputs for %s\n",co.name)}
  }
}

//$s1=version string, loads the data created from minrunsv - should set strv first
proc loadminrundat () { localobj st
  {st=new String2() sprint(st.s,"data/%s_.vecs",$s1) sprint(st.t,"data/%s_nclnq.nqs",$s1)}
  if(FileExists(st.t)) loadmydat(st.s,st.t) else loadmydat(st.s)
}