// $Id: parnqsnet.hoc,v 1.2 2006/02/11 13:34:17 hines Exp $
// load_file("nqsnet.hoc")
objref sp[5],nq[5]
datestr="06feb04"
scale=1 // redefined in network.hoc
double pmat[CTYPi][CTYPi][STYPi],numc[CTYPi],ix[CTYPi],wmat[CTYPi][CTYPi][STYPi]
obfunc cobj(){} // return cell given gid (defined in network.hoc)
//* routines
//** styp() sets synapse type based on presynaptic cell
func styp () { local pr,po
pr=$1 po=$2
if (pr==IN && po==IN) { return GA
} else if (pr==IN) { return IH
} else if (pr==SU) { return EX
} else if (pr==SM) { return AM
} else printf("styp ERR %s->%s not classified",CTYP.object(pr).s,CTYP.object(po).s)
}
for ii=0,2 {sp[ii]=new NQS("PR","PO","SID","WT","DEL") sp[ii].clear(1e5*scale)}
batch_flag=1
for ii=0,2 sp[ii].mo(1)
batch_flag=0
// mkspmat() uses sp and sp[1], returns sp
proc mkspmat () { local a,sn,ii,jj,n,x,y,z,z1 localobj vr,vo
x=$1 y=$2 z=$3
rdm.ACG(903342)
a=allocvecs(vr,vo)
sn=pmat[x][y][z]*numc[x]*numc[y] // number of synapses
for ii=0,1 sp[ii].clear
sp.mo(2)
sp.pad(1.5*sn) // leave room
rdm.discunif(ix[x],ix[x]+numc[x]-1) // randomly chosen presynaptic cells
PRv.setrand(rdm)
rdm.discunif(ix[y],ix[y]+numc[y]-1) // randomly chosen postsynaptic cells
POv.setrand(rdm)
if (z==EX) z1=AM else z1=z
sytyp(c[y],z1,vr) // will want to replace this since got rid of c[]
SIDv.vfill(vr)
rdm.normal(wmat[x][y][z1],(0.5*wmat[x][y][z1])^2) // calc random weights
WTv.setrand(rdm) // set the weight column
rdm.normal(2.2,0.1) // random delays
DELv.setrand(rdm)
sp.elimrepeats("PR","PO","SID") // so don't have 2 connections between same cells
sp.shuffle() // randomize again
sp.pad(sn) // reduce the length down to desired
if (z==EX) { // now have to repeat everything, also will need for IX
z1=NM
sp[1].mo(2) // reset the PRv,POv,... pointers
sp[1].cp(sp) // this will be the same as prior
rdm.normal(wmat[x][y][z1],(0.5*wmat[x][y][z1])^2)
WTv.setrand(rdm)
sytyp(c[y],z1,vo) // vo and vr represent list indices for cell synlists
vo.sub(vr) // difference between the NMDA synnums and AMPA synnums
if (! vo.ismono(0)){ printf("mkspmat ERRA\n") return }
SIDv.add(vo.x[0])
sp.grow(sp[1])
}
sp.mo(2) // reset the PRv,POv,... pointers
WTv.negwrap() // take any neg weights and make them positive
dealloc(a)
}
// mkpomat() only make the array for a single postsynaptic cell
proc mkpomat () { local a,sn,ii,jj,n,x,y,yid,z,z1 localobj vr,vo
x=$1 y=$2 z=$3 yid=$4
// rdm.ACG(903342)
// We want to have a postsynaptic gid specific stream but streams from
// different gids should be statistically independent.
// A stream is random but we want to be able to reproduce it at will
// Assume this is the only place where this is required
mcell_ran4_init(yid) // for independent sims add offset > i*allcells
rdm.MCellRan4(1) //avoid 0 want to always start at known seed
a=allocvecs(vr,vo)
sn=pmat[x][y][z]*numc[x]*1 // number of synapses
if (sn<1) { // treat sn<1 as a prob of connecting rather than a density
if (rdm.uniform(0,1)<sn) sn=1 else sn=0 // flip coin
} else if (sn<5) sn=round(sn) // possibly round it up
for ii=0,1 sp[ii].clear
if (sn==0) { dealloc(a) return }
sp.mo(2)
sp.pad(1.5*sn) // leave room
rdm.discunif(ix[x],ix[x]+numc[x]-1) // randomly chosen presynaptic cells
PRv.setrand(rdm)
POv.fill(yid)
if (z==EX) z1=AM else z1=z
sytyp(cobj(yid),z1,vr) // brief syncode.hoc routine translates synaptic types
SIDv.vfill(vr)
rdm.normal(wmat[x][y][z1],(0.5*wmat[x][y][z1])^2) // calc random weights
WTv.setrand(rdm) // set the weight column
rdm.normal(2.2,0.1) // random delays
DELv.setrand(rdm)
sp.elimrepeats("PR","PO") // so don't have 2 connections between same cells
sp.shuffle() // randomize again
sp.pad(sn) // reduce the length down to desired
if (z==EX) { // now have to repeat everything, also will need for IX
z1=NM
sp[1].mo(2) // reset the PRv,POv,... pointers
sp[1].cp(sp) // this will be the same as prior
rdm.normal(wmat[x][y][z1],(0.5*wmat[x][y][z1])^2)
WTv.setrand(rdm)
sytyp(cobj(yid),z1,vo) // vo and vr represent list indices for cell synlists
vo.sub(vr) // difference between the NMDA synnums and AMPA synnums
if (! vo.ismono(0)){ printf("mkspmat ERRA\n") return }
SIDv.add(vo.x[0])
sp.grow(sp[1])
}
sp.mo(2) // reset the PRv,POv,... pointers
WTv.negwrap() // take any neg weights and make them positive
dealloc(a)
}
proc mkall () { local x,y,z localobj fname
fname=new String()
for case(&y,FP,TP,B5) { // FP,TP,B5
sp[2].clear
for case(&x,SM,FP,TP,B5) for case(&z,AM,GA,GB,EX,E2) {
if (pmat[x][y][z]!=0) {
mkspmat(x,y,z)
sp[2].grow(sp)
}
}
sprint(fname.s,"data/%s%sS%02d.nqs",datestr,CTYP.o(y).s,scale)
sp[2].sv(fname.s)
}
}
// nqs2txt()
proc nqs2txt () { local ii,last,lnum,n,fl1 localobj f,fn
nq=new NQS("GID","LINE","NUML","TELL")
fn=new String()
fl1=lnum=n=0
f=new File()
sprint(tstr,"data/%sS%02d.syns",datestr,scale)
f.wopen(tstr)
for case(&y,FP) { // FP,TP,B5 need proper order
sprint(fn.s,"data/%s%sS%02d.nqs",datestr,CTYP.o(y).s,scale)
sp.rd(fn.s)
sp.sort("SID") sp.sort("PR") sp.sort("PO")
if (fl1==0) {
last=sp.v[1].x[0]
nq.append(last,lnum+1,-1,f.tell)
fl1=1
}
for ii=0,sp.size(1)-1 {
PR=sp.v[0].x[ii] PO=sp.v[1].x[ii] SID=sp.v[2].x[ii] WT=sp.v[3].x[ii] DEL=sp.v[4].x[ii]
if (PO==last) {
lnum+=1 n+=1
f.printf("%d %d %d %g %g\n",PO,PR,SID,WT,DEL)
} else {
nq.set("NUML",-1,n)
n=0
if (PO==last+1) {
nq.append(PO,lnum+1,-1,f.tell+1)
} else {
printf("%d is not a target\n")
nq.append(PO,lnum+1,0,f.tell+1)
}
last=PO
ii-=1 // do this one again
}
}
nq.set("NUML",-1,n)
}
sprint(tstr,"data/%sS%02d.index",datestr,scale)
f.wopen(tstr)
for ii=0,nq.size(1)-1 {
f.printf("%d %d %d %d\n",nq.v[0].x[ii],nq.v[1].x[ii],nq.v[2].x[ii],nq.v[3].x[ii])
}
f.close
}
// nqssplit() // divide into multiple smaller files
proc nqssplit () { local ii,last,lnum,n,fl1 localobj f,fn
fn=new String()
fl1=lnum=n=0
f=new File()
for case(&y,FP,TP,B5) {
sprint(fn.s,"data/%s%sS%02d.nqs",datestr,CTYP.o(y).s,scale)
sp.rd(fn.s)
sp.sort("PO")
for (ii=ix[y];ii<ix[y]+numc[y];ii+=200) {
sp.select("PO","[]",ii,ii+200-1)
sprint(tstr,"data/S%dzip/%s_%d-%d_S%02d.nqs",scale,datestr,ii,ii+200-1,scale)
sp.out.sv(tstr)
}
}
}
proc ncl2mat () { local ii,wt,del,pc,pr,po,prid,poid,ps,sy,syn localobj Xo
if (!isassigned(sp)) {
sp=new NQS("PRID","POID","PR","PO","SYNID","DEL","WT0","THRESH") }
for ii=0,ncl.count-1 {
Xo=ncl.o(ii)
pc=id(Xo.pre) prid=id(Xo.precell) sy=id(Xo.syn) poid=id(Xo.postcell)
pr=ojtnum(Xo.precell) po=objnum(Xo.postcell)
if (!isojt(Xo.postcell,nil)) syn=Xo.postcell.synlist.index(Xo.syn)
if (poid==-1 && sy==-1) { // nil targ
} else if (poid==-1) { // ACELL targ -- not implemented
} else if (prid==-1) { // ACELL src
sp.append(pc,poid,ojtnum(Xo.pre),po,syn,Xo.delay,Xo.weight,Xo.threshold)
} else {
sp.append(prid,poid,pr, po,syn,Xo.delay,Xo.weight,Xo.threshold)
}
}
}
//* fill everything in and save it
// mksfc() // synfire chain
// ncl2mat()
// mkall()
// go to global id's -- can leave fp alone since starts at 0
proc rel2glbl () {
sp.select("PRID",TP) sp.spr("<PR>.add(idtp)") sp.delect()
sp.select("POID",TP) sp.spr("<PO>.add(idtp)") sp.delect()
sp.select("PRID",B5) sp.spr("<PR>.add(idb5)") sp.delect()
sp.select("POID",B5) sp.spr("<PO>.add(idb5)") sp.delect()
sp.select("PRID",SM) sp.spr("<PR>.add(idstim)") sp.delect()
}
