// $Id: network.hoc,v 1.415 2012/07/19 00:25:46 samn Exp $
//* Numbers and connectivity params
DOPE_INTF6=1
EDOPE_INTF6=1
ESTDP_INTF6 = ISTDP_INTF6 = 0
declare("scale",4)
declare("colW",1,"colH",1)
declare("numcols",colW*colH)
declare("dbgcols",0) // whether to debug columns by making them have the same wiring and inputs
declare("colr",2) // maximal trans-column projection distance; 0 within col; 1 next col etc
declare("colnq","o[5]","lcol",new List())
declare("autotune",0) //
declare("useSTDP",1) //
{sprint(tstr,"o[%d]",numcols) declare("col",tstr)}
{sprint(tstr,"o[%d][%d]",colH,colW) declare("gcol",tstr)} // 2D column grid
double div[numcols][numcols][CTYPi][CTYPi]//div[i][j]==# of outputs from type i->j
double wmat[numcols][numcols][CTYPi][CTYPi][STYPi] // wmat[i][j][k]==weight from type i->j for synapse k
double delm[numcols][numcols][CTYPi][CTYPi]//avg. delay from type i->j
double deld[numcols][numcols][CTYPi][CTYPi]//delay variance from type i->j
double conv[numcols][numcols][CTYPi][CTYPi]
dosetpmat=name_declared("pmat")==0
{sprint(tstr,"d[%d][%d][%d][%d]",numcols,numcols,CTYPi,CTYPi) declare("pmat",tstr)}
double synloc[CTYPi][CTYPi]//location of synapses
declare("cedp",new List()) // list for DP cells
declare("ncdp",new List()) // list of NetCons from DP -> X cells
declare("useEV",0) // whether to use EV cells for target select
//* wiring & related column variables
declare("wirety",4) // type of wiring to use, 0=fixed divergence, 1=swire, 2=swirecut, 3=swirecutfl, 4=fixed convergence (convwire)
// NB: wirty==4 (convwire) still has spatial wiring from DP -> ES
declare("colside",30) // column diameter in micrometers
declare("checkers",0) // whether to arrange cells in checkerboard pattern
declare("gridpos",0) // whether to arrange INTF6 cells on a 2D grid
declare("DPgridpos",0) // whether to arrange DP cells on a 2D grid
declare("DPpad",0) // whether to pad DP cells when they're separated
declare("shareI",0) // whether to share I cells between S,M cells
declare("cxinc",3,"cyinc",3,"crad",1)
declare("EEsq",(colside/2)^2,"EIsq",(colside/2)^2,"IEsq",(colside/1)^2,"IIsq",(colside/2)^2) // params for swirecut
double dels[CTYPi][CTYPi]//[STYPi] //stdev of delays
double delv[CTYPi][CTYPi]//variance of delays
double vcond[CTYPi][CTYPi]//[STYPi] //conduction velocities
declare("layerzvar",25) // range of z location of cells in a layer (in micrometers)
declare("colxydist",50) // distance between adjacent columns (x,y)
declare("vlayerz",new Vector(CTYPi))
declare("pseed",4321)
declare("dvcond","d[2]") //intra layer conduction delays
declare("maxsfall",0.001) // max fall-off in prob of connection @ opposite side of column (used by swire)
declare("slambda",colside/3) // spatial length constant for probability of connections, used in swirecut
dvcond[0] = 1 //if presynaptic cell is excit, conduction delay == 1.0 m/s
dvcond[1] = 0.4 //if presynaptic cell is inhib, conduction delay == 0.4 m/s
declare("hval",0)
//* other variables
declare("mknetnqss",1) // whether COLUMN should make connsnq,cellsnq
declare("fc",0.55)
declare("EEGain",2,"EIGain",8.5,"IEGain",3,"IIGain",3,"ELTSGain",0.5)
declare("NMAMR",0.1,"EENMGain",1)
declare("EIGainInterC0",0.125,"EIGainInterC1",1,"EEGainInterC0",1,"EEGainInterC1",0.25)
double EIGainInterC[2],EEGainInterC[2]
EIGainInterC[0] = EIGainInterC0 // 0->1
EIGainInterC[1] = EIGainInterC1
EEGainInterC[0] = EEGainInterC0
EEGainInterC[1] = EEGainInterC1
declare("TCNM",0) // whether TC->NM2 on
declare("NMUSCLES",2) // number of muscle groups in model
declare("sepDPpos",1) // whether to position DP cells in each group separately
declare("DPESGain",15) // 11.89875, weight gain from DP -> ES cells
if(autotune) {
seadsetting_INTF6 = 2
wsetting_INTF6 = 1 // use sywv during sim
} else if(useSTDP) {
seadsetting_INTF6 = 3 // use plasticity
} else {
seadsetting_INTF6 = 2 // fixed weights
}
declare("wnqstr","") // intracol wiring NQS (single column)
//* prdiv() - print div
proc prdiv () { local ii,jj,kk,ll
for ii=0,numcols-1 for jj=0,numcols-1 for kk=0,CTYPi-1 for ll=0,CTYPi-1 if(div[ii][jj][kk][ll]) {
printf("div[%d][%d][%s][%s]=%g\n",ii,jj,CTYP.o(kk).s,CTYP.o(ll).s,div[ii][jj][kk][ll])
}
}
// %con (con/pre) = %div (div/post)
DEAD_DIV_INTF6=0
declare("jcn",1)
declare("disinhib",0) //iff==1 , turn off inhibition, by setting wmat[I%d][...]==0 in inhiboff()
declare("pmatscale",0.75*6.0/scale) // scale for pmat - allows keeping it fixed while changing # of cells in network
declare("wmatscale",1) // scale for wmat - called after setwmat
batch_flag=declare("dstr",datestr,"setdviPT",NORM)
declare("params","not batch","ofile",output_file)
declare("dvseed",534023) // seed for wiring
dosetcpercol=name_declared("cpercol")==0 // whether to set values in cpercol or use user-supplied values
{sprint(tstr,"d[%d][%d]",numcols,CTYPi) declare("cpercol",tstr)} // cells of a specific type per column
{sprint(tstr,"o[%d]",numcols) declare("vcpercol",tstr)}
for i=0,numcols-1 vcpercol[i]=new Vector(CTYPi)
declare("DPESFeedForward",0.1) // level of feedforward connectivity from DP -> ES
declare("EMRecur",0.05) // level of recurrence between EM cells
declare("ESRecur",0.05)// level of recurrence between ES cells
declare("EMESFeedBack",0.017)// level of feedback from EM -> ES
declare("ESEMFeedForward",0.08)// level of feedforward connectivity from ES -> EM
declare("ESIMFeedForward",0.0)
declare("EMISFeedBack",0.0)
declare("ESEMGain",1.0) // weight gain from ES -> EM cells
sprint(tstr,"d[%d][%d]",CTYPi,CTYPi)
if(!(i=name_declared("delmscale"))) {
declare("delmscale",tstr) // scale delm values by this #
}
if(!i) {
for i=0,CTYPi-1 for j=0,CTYPi-1 delmscale[i][j]=1
}
//** mklayerz -- get average z location of a layer, based on frana and tononi
// returns value in micrometers. 0 is at top (layer 1), max val
// is at layer 6 ... actual #s don't matter, only distances between
// the layers matters...
proc mkvlayerz () { local i,vl
vlayerz.indgen(0,CTYPi-1,1)
vlayerz.apply("layer")
for vtr(&vl,vlayerz,&i) {
if(vl == 2 || vl == 3){
vlayerz.x(i)=1740 // average of 1540+1940 from frana
} else if(vl == 4) {
vlayerz.x(i) = (1740+1130) / 2
} else if(vl == 5 ){
vlayerz.x(i) = 1130
} else if(vl == 6 ){
vlayerz.x(i) = 488
} else vlayerz.x(i) = 200
}
}
//* setcpercol - set # of cells per column
proc setcpercol () { local i,j // (/u/samn/vcsim/notebook.dol_1:24562)(notebook.dol_1:24492)
if(dosetcpercol) { // if user didn't supply values (default), set # of cells of a type per column
cpercol[0][ES] = int(48 * scale) // sensory layer
cpercol[0][ISL] = 5 * scale
cpercol[0][IS] = 11 * scale
cpercol[0][EM] = int(48 * scale) // motor layer
cpercol[0][IML] = 5 * scale
cpercol[0][IM] = 11 * scale
if(shareI) for i=0,CTYPi-1 if(ice(i)) cpercol[0][i] *= shareI // if shareI!=0, it scales # of I cells
if(useEV) cpercol[0][EV] = int(48*scale)
cpercol[0][DP] = cpercol[0][ES] // arm-cell/spindle layer
for i=0,CTYPi-1 if(cpercol[0][i]) cpercol[0][i]=int(cpercol[0][i]+0.5) // round
//store values in vcpercol Vectors
for i=0,numcols-1 {
{vcpercol[i].resize(CTYPi) vcpercol[i].fill(0)} // store the values in a vector
for j=0,CTYPi-1 vcpercol[i].x(j)=cpercol[i][j]
}
}
}
//* setpmat()
proc setpmat () { local pre,po,ii,jj,kk,ll
if(!dosetpmat) return // if pmat setup by user (in notebook), then don't reset its values
for ii=0,numcols-1 for jj=0,numcols-1 for kk=0,CTYPi-1 for ll=0,CTYPi-1 pmat[ii][jj][kk][ll]=0
// intracol connections
ii=0
pmat[0][0][DP][ES]=DPESFeedForward // 0.1, DP = drawspike cells (muscle/arm layer)
// sensory layer
pmat[ii][ii][ES][EM]=ESEMFeedForward
pmat[ii][ii][ES][IM]=pmat[ii][ii][ES][IML]=ESIMFeedForward
pmat[ii][ii][ES][ES]=ESRecur
pmat[ii][ii][ES][ISL]=0.51
pmat[ii][ii][ES][IS]=0.43
pmat[ii][ii][ISL][ES]=0.35
pmat[ii][ii][ISL][ISL]=0.09
pmat[ii][ii][ISL][IS]=0.53
pmat[ii][ii][IS][ES]=0.44
pmat[ii][ii][IS][ISL]=0.34
pmat[ii][ii][IS][IS]=0.62
// motor layer
pmat[ii][ii][EM][EM]=EMRecur
pmat[ii][ii][EM][ES]=EMESFeedBack
pmat[ii][ii][EM][IS]=pmat[ii][ii][EM][ISL]=EMISFeedBack
pmat[ii][ii][EM][IML]=0.51
pmat[ii][ii][EM][IM]=0.43
pmat[ii][ii][IML][EM]=0.35
pmat[ii][ii][IML][IML]=0.09
pmat[ii][ii][IML][IM]=0.53
pmat[ii][ii][IM][EM]=0.44
pmat[ii][ii][IM][IML]=0.34
pmat[ii][ii][IM][IM]=0.62
if(shareI) { // share interneurons between S,M cells?
pmat[ii][ii][EM][ISL]=0.51
pmat[ii][ii][EM][IS]=0.43
pmat[ii][ii][ISL][EM]=0.35
pmat[ii][ii][IS][EM]=0.44
pmat[ii][ii][ES][IML]=0.51
pmat[ii][ii][ES][IM]=0.43
pmat[ii][ii][IML][ES]=0.35
pmat[ii][ii][IM][ES]=0.44
pmat[ii][ii][ISL][IML]=0.09
pmat[ii][ii][ISL][IM]=0.53
pmat[ii][ii][IS][IML]=0.34
pmat[ii][ii][IS][IM]=0.62
pmat[ii][ii][IML][ISL]=0.09
pmat[ii][ii][IML][IS]=0.53
pmat[ii][ii][IM][ISL]=0.34
pmat[ii][ii][IM][IS]=0.62
}
if(useEV) pmat[ii][ii][EV][EM]=0.25
}
//* scalepmat(fctr) - multiply values in pmat by fctr
proc scalepmat () { local fctr,from,to,c1,c2
fctr=$1
for c1=0,numcols-1 for c2=0,numcols-1 {
for from=0,CTYPi-1 for to=0,CTYPi-1 pmat[c1][c2][from][to] *= fctr
}
}
//* pmat2nq - return an NQS with info in pmat
obfunc pmat2nq () { local i,j,k,l localobj nqpmat
nqpmat=new NQS("col1","col2","froms","tos","from","to","pij")
{nqpmat.strdec("froms") nqpmat.strdec("tos")}
for i=0,numcols-1 for j=0,numcols-1 for k=0,CTYPi-1 for l=0,CTYPi-1 if(pmat[i][j][k][l]) {
nqpmat.append(i,j,CTYP.o(k).s,CTYP.o(l).s,k,l,pmat[i][j][k][l])
}
return nqpmat
}
//* nq2pmat - load NQS ($o1) into pmat
proc nq2pmat () { local i,j,k,l localobj nq,vf,vto,vc1,vc2,vpij
{nq=$o1 nq.tog("DB") vc1=nq.getcol("col1") vc2=nq.getcol("col2")}
{vf=nq.getcol("from") vto=nq.getcol("to") vpij=nq.getcol("pij")}
for i=0,numcols-1 for j=0,numcols-1 for k=0,CTYPi-1 for l=0,CTYPi-1 pmat[i][j][k][l]=0
for i=0,vf.size-1 pmat[vc1.x(i)][vc2.x(i)][vf.x(i)][vto.x(i)]=vpij.x(i)
print "loaded " , nq , " into pmat"
}
//* synapse locations DEND SOMA AXON
proc setsynloc () { local from,to
for from=0,CTYPi-1 for to=0,CTYPi-1 {
if(ice(from)) {
if(IsLTS(from)) {
synloc[from][to]=DEND // distal [GA2] - from LTS
} else {
synloc[from][to]=SOMA // proximal [GA] - from FS
}
} else {
synloc[from][to]=DEND // E always distal. use AM2,NM2
}
}
}
//* setdelmats -- setup delm,deld
proc setdelmats () { local from,to,ii,jj,fcm,fcd
for ii=0,numcols-1 for jj=0,numcols-1 {
if(ii==jj) {
fcm=fcd=1
} else {fcm=3 fcd=6}
for from=0,CTYPi-1 for to=0,CTYPi-1 {
if(synloc[from][to]==DEND) {
delm[ii][jj][from][to]=4 * delmscale[from][to] * fcm
deld[ii][jj][from][to]=1 * fcd
} else {
delm[ii][jj][from][to]=2.0 * delmscale[from][to] * fcm
deld[ii][jj][from][to]=0.2 * fcd
}
}
}
// delm[0][0][DP][ES] = 8
// deld[0][0][DP][ES] = 4
}
//* weight params
//** delay all 2+/-0.02 within column for now
proc setwmat () { local from,to,sy,gn,c,ii,jj
for ii=0,numcols-1 for jj=0,numcols-1 for from=0,CTYPi-1 for to=0,CTYPi-1 for sy=0,STYPi-1 wmat[ii][jj][from][to][sy]=0
ii=0
// wire DP directly to ES
wmat[0][0][DP][ES][AM2]= DPESGain / EEGain
// sensory layer
wmat[ii][ii][ES][ES][AM2]=0.66
wmat[ii][ii][ES][ISL][AM2]=0.23 // * 1.375 // * 1.375 * 1.25
wmat[ii][ii][ES][IS][AM2]=0.23 // * 1.25 * 1.25
wmat[ii][ii][ES][EM][AM2]=0.88 * ESEMGain
wmat[ii][ii][IS][ES][GA]=1.5
wmat[ii][ii][IS][ISL][GA]=1.5
wmat[ii][ii][IS][IS][GA]=1.5
wmat[ii][ii][ISL][ES][GA2]=0.83
wmat[ii][ii][ISL][ISL][GA2]=1.5
wmat[ii][ii][ISL][IS][GA2]=1.5
// motor layer
wmat[ii][ii][EM][EM][AM2]=0.66 * 0.9
wmat[ii][ii][EM][IML][AM2]=0.23 //* 0.8 // * 1.8 * 1.25
wmat[ii][ii][EM][IM][AM2]=0.23 //* 0.7 // * 1.5 * 1.25
wmat[ii][ii][EM][ES][AM2]=0.24
wmat[ii][ii][IM][EM][GA]=1.5 * 2
wmat[ii][ii][IM][IM][GA]=1.5
wmat[ii][ii][IM][IML][GA]=1.5
wmat[ii][ii][IML][EM][GA2]=0.83 * 2
wmat[ii][ii][IML][IM][GA2]=1.5
wmat[ii][ii][IML][IML][GA2]=1.5
if(shareI) { // share interneurons between S,M cells?
wmat[ii][ii][ES][IML][AM2]=0.23
wmat[ii][ii][ES][IM][AM2]=0.23
wmat[ii][ii][IM][ES][GA]=1.5
wmat[ii][ii][IML][ES][GA2]=0.83
wmat[ii][ii][EM][ISL][AM2]=0.23
wmat[ii][ii][EM][IS][AM2]=0.23
wmat[ii][ii][IS][EM][GA]=1.5
wmat[ii][ii][ISL][EM][GA2]=0.83
wmat[ii][ii][IS][IML][GA]=1.5
wmat[ii][ii][IS][IM][GA]=1.5
wmat[ii][ii][ISL][IML][GA2]=1.5
wmat[ii][ii][ISL][IM][GA2]=1.5
wmat[ii][ii][IM][ISL][GA]=1.5
wmat[ii][ii][IM][IS][GA]=1.5
wmat[ii][ii][IML][ISL][GA2]=1.5
wmat[ii][ii][IML][IS][GA2]=1.5
}
if(ESIMFeedForward) {
wmat[ii][ii][ES][IML][AM2]=0.23
wmat[ii][ii][ES][IM][AM2]=0.23
}
if(EMISFeedBack) {
wmat[ii][ii][EM][ISL][AM2]=0.23
wmat[ii][ii][EM][IS][AM2]=0.23
}
if(useEV) wmat[ii][ii][EV][EM][AM2]=0.5
//set NMDA weights
for ii=0,numcols-1 for jj=0,numcols-1 {
for from=0,CTYPi-1 for to=0,CTYPi-1 wmat[ii][jj][from][to][NM2]=NMAMR*wmat[ii][jj][from][to][AM2]
}
//gain control
for ii=0,numcols-1 for jj=0,numcols-1 for from=0,CTYPi-1 for to=0,CTYPi-1 for sy=AM,GA2 if(wmat[ii][jj][from][to][sy]>0) {
if(ice(from)) {
if(ice(to)) {
gn = IIGain
} else {
gn = IEGain
}
if(IsLTS(from) && !IsLTS(to)) gn *= 0.5
} else {
if(ice(to)) {
gn = EIGain
if(IsLTS(to)) gn *= ELTSGain
} else {
gn = EEGain
if(sy==NM || sy==NM2) gn *= EENMGain // E->E NMDA gain
}
}
wmat[ii][jj][from][to][sy] *= gn
}
}
//* scalewmat(fctr) - multiply values in wmat by fctr
proc scalewmat () { local fctr,from,to,c1,c2,sy
fctr=$1
for c1=0,numcols-1 for c2=0,numcols-1 {
for from=0,CTYPi-1 {
if(from==DP || IsTHAL(from)) continue
for to=0,CTYPi-1 {
if(IsTHAL(to)) continue
for sy=0,STYPi-1 wmat[c1][c2][from][to][sy] *= fctr
}
}
}
}
//* %con (con/pre) = %div (div/post)
{sp = new NQS() cp = new NQS()}
//** gethublims(col,hubtype,hubfactor,numhubs,mode)
// get a matrix of size CTYPi X CTYPi, specifying div with mat.x(hubtype,othertype)
// and conv with mat.x(othertype,hubtype)
// hubtype = type of hub. hubfactor = desired ratio of hub div/conv vs non-hub div/conv
// numhubs = # of hubs. col = COLUMN for which to set hubs.
// mode == 0 <-- hub div(conv) is set to hubfactor * original div(conv)
// mode == 1 <-- hub div(conv) is set so that final hub div = hubfactor * final non_hub div (same for conv)
// formula is based on: m / ((N-H*m) / (C-H)) = F , and then solving for m
// m = div for the hubs, F = desired ratio of final hub div to final non-hub div
// N = # of synapses (links), C = total # of postsynaptic cells (including hubs) , H = # of hubs
// similarly done for conv , but replace N with appropriate values
// (/u/samn/intfcol/notebook.dol_1:21933)
obfunc gethublims () { local ct,mode,from,to,lim,nc,nhubs,fctr localobj col,mat
{col=$o1 ct=$2 fctr=$3 nhubs=$4 mode=$5 mat=new Matrix(CTYPi,CTYPi)}
for to=0,CTYPi-1 if(col.numc[to] && col.div[ct][to]) {
{nc=col.numc[to] if(ct==to)nc-=1} // deduct for self-link
if(mode==0) {
lim = int( 0.5 + col.div[ct][to]*fctr )
} else {
lim = int( 0.5 + col.div[ct][to]*col.numc[ct]*fctr/(col.numc[ct]-nhubs+fctr*nhubs) )
}
mat.x(ct,to) = MINxy(lim, nc) // at most div to all postsynaptic cells
}
for from=0,CTYPi-1 if(col.numc[from] && col.div[from][ct]) {
{nc=col.numc[from] if(ct==from)nc-=1} // deduct for self-link
if(mode==0) {
lim = int( 0.5 + col.conv[from][ct]*fctr )
} else {
lim = int( 0.5 + col.div[from][ct]*col.numc[from]*fctr/(col.numc[ct]-nhubs+fctr*nhubs) )
}
mat.x(from,ct) = MAXxy(MINxy(lim, nc),1) // at most conv from all presynaptic cells, but at least 1
}
return mat
}
//** addhubs(column,cell-type,numhubs,scaling factor,skipI[,seed,allowz,hubmode,verbose])
// add hubs to the network by stealing wires from other neurons
// $o1 == column
// $2 == cell type of hub
// $3 == number of hubs to add
// $4 == scaling factor (should be > 1.0) for conv,div of hub
// $5 == skip div/conv of I cells
// $6 == seed - optional
// $7 == allowz - whether to allow pulling all links from/to another cell
// $8 == hubmode - which mode to use for gethublims (see above)
// $9 == verbose - optional
// function returns a Vector containing the ids of the cells selected as hubs (within column ids)
obfunc addhubs () { local a,ct,fctr,nhubs,idx,jdx,lseed,hubid,szorig,cursz,preid,poid,lim,skipI,to,from,vrb,changed,allowz,hmode\
localobj col,ce,vin,vout,nq,vd,vc,vdd,vdt,vddt,vpicked,vhubid,vw1,vw2,vsyn,vprob,vsynt,vtmp,vdsz,vcsz,mhlim
col=$o1 ct=$2 nhubs=$3 fctr=$4 skipI=$5
if(numarg()>5) lseed=$6 else lseed=1234
if(numarg()>6) allowz=$7 else allowz=1
if(numarg()>7) hmode=$8 else hmode=0
if(numarg()>8) vrb=$9 else vrb=0
{ce=col.ce hashseed_stats(lseed,lseed,lseed)}
a=allocvecs(vin,vout,vd,vc,vdd,vdt,vddt,vpicked,vw1,vw2,vsyn,vprob,vsynt,vtmp,vdsz,vcsz)
vrsz(col.allcells,vin,vout,vd,vc,vdd,vdt,vddt,vpicked,vw1,vw2,vsyn,vprob,vsynt,vdsz,vcsz,vtmp)
mhlim=gethublims(col,ct,fctr,nhubs,hmode)
//vin,vout = input/output markers. vd,vc = div/conv.
//vdd div/conv delays, vdt div/conv temp. vddt=div/conv delay temp
//vpicked=which cells already picked as hubs
vhubid=new Vector()
{vhubid.indgen(col.ix[ct],col.ixe[ct],1) vhubid.shuffle() vhubid.resize(nhubs)}
if(vrb) vlk(vhubid)
for idx=0,vhubid.size-1 vpicked.x(vhubid.x(idx))=1
for idx=0,vhubid.size-1 { hubid=vhubid.x(idx)
if(vrb) printf("hub%d id = %d\n",idx+1,hubid)
{ce.o(hubid).getdvi(vd,vdd,vw1,vw2,vprob,vsyn) ce.o(hubid).getconv(vc)}//IDs of post/presynaptic cells
{ce.o(hubid).getconv(1.2,vcsz) vdsz.resize(CTYPi) vdsz.fill(0)}//counts of post/pre types
for jdx=0,vd.size-1 vdsz.x(ce.o(vd.x(jdx)).type)+=1
{vout.fill(0) vin.fill(0)} //init as 0
for jdx=0,vd.size-1 vout.x(vd.x(jdx))=1 //mark current postsynaptic cells
for jdx=0,vc.size-1 vin.x(vc.x(jdx))=1 //mark current presynaptic cells
for to=0,CTYPi-1 if(col.numc[to] && col.div[ct][to] && (!skipI || !ice(to))) {
cursz=szorig=vdsz.x(to) // update divergence
if(vrb) print "\torig div -> " , CTYP.o(to).s, " = " , szorig
{lim=mhlim.x(ct,to) changed=1}
while(cursz<lim && changed==1) { changed=0
for(preid=col.ix[ct];preid<=col.ixe[ct] && cursz<lim;preid+=1) {// pick same presynaptic type
if(vpicked.x(preid)) continue //dont take from other hubs
ce.o(preid).getdvi(vdt,vddt,vw1,vw2,vprob,vsynt)
vtmp.fill(0)
for jdx=0,vdt.size-1 vtmp.x(ce.o(vdt.x(jdx)).type)+=1
if(!allowz && vtmp.x(to)<=1)continue//dont want to turn div of another cell to 0
for jdx=0,vdt.size-1 { poid=vdt.x(jdx) // go thru postsynaptic cells looking for target type
if(ce.o(poid).type==to && poid!=hubid && vout.x(poid)==0) { cursz+=1
{vd.append(poid) vdd.append(vddt.x(jdx)) vsyn.append(vsynt.x(jdx))}
{vdt.remove(jdx) vddt.remove(jdx) vsynt.remove(jdx)}
ce.o(preid).setdvi(vdt,vddt,vsynt,1) // update presynaptic cell
vout.x(poid)=changed=1 // this cell synapses on poid
break
}
}
}
}
if(vrb) print "\tnew div -> " , CTYP.o(to).s, " = " , cursz
}
ce.o(hubid).setdvi(vd,vdd,vsyn,1) // update hub dvi
for from=0,CTYPi-1 if(col.numc[from] && col.div[from][ct] && (!skipI || !ice(from))) {
cursz=szorig=vcsz.x(from) // update convergence
{lim=mhlim.x(from,ct) changed=1}
if(vrb) print "\torig conv <- ", CTYP.o(from).s, " = " , szorig
while(cursz<lim && changed==1) { changed=0
for(preid=col.ix[from];preid<=col.ixe[from]&&cursz<lim;preid+=1) {
if(preid==hubid || vin.x(preid)) continue // don't make self or double-connects
ce.o(preid).getdvi(vdt,vddt,vw1,vw2,vprob,vsynt)
for jdx=0,vdt.size-1{
poid = vdt.x(jdx)
if(vpicked.x(poid)) continue // don't take wires from other hubs
if(ce.o(poid).type==ct){ ce.o(poid).getconv(1.2,vtmp)
if(allowz || vtmp.x(from)>1) { // make sure not to remove all inputs of a type to a cell
vdt.x( jdx ) = hubid // reassign input to hub
ce.o(preid).setdvi(vdt,vddt,vsynt,1) // reset presynaptic cell's div
vin.x( preid ) = changed = 1 // mark input
cursz += 1
break
}
}
}
}
}
if(vrb) print "\tnew conv <- " , CTYP.o(from).s, " = " , cursz
}
}
ce.o(0).finishdvir()
{dealloc(a) return vhubid}
}
//* mkcolnqs - make an nqs with current pmat,wmat,delm,deld info for use by both COLUMNs(AREAS) for wiring
// "dist" represents distance between columns: dist==0 for intra-COLUMN setup, dist>0 for INTER-COLUMN setup
proc mkcolnqs () { local from,to,sy,i,j,d localobj froms,tos,sys
for i=0,numcols-1 {
{nqsdel(colnq[i]) colnq[i]=new NQS("froms","tos","sys","from","to","sy","w","pij","delm","deld","loc","dist")}
colnq[i].strdec("froms","tos","sys")
d=0 //intracol
for from=0,CTYPi-1 { froms=CTYP.o(from)
for to=0,CTYPi-1 { tos=CTYP.o(to)
if(pmat[i][i][from][to]>0) for sy=0,STYPi-1 if(wmat[i][i][from][to][sy]>0) { sys=STYP.o(sy)
colnq[i].append(froms.s,tos.s,sys.s,from,to,sy,wmat[i][i][from][to][sy],pmat[i][i][from][to],delm[i][i][from][to],deld[i][i][from][to],synloc[from][to],d)
}
}
}
for j=0,numcols-1 if(i!=j) { d=1 //intercol
for from=0,CTYPi-1 { froms=CTYP.o(from)
for to=0,CTYPi-1 { tos=CTYP.o(to)
if(pmat[i][j][from][to]>0) for sy=0,STYPi-1 if(wmat[i][j][from][to][sy]>0) { sys=STYP.o(sy)
colnq[i].append(froms.s,tos.s,sys.s,from,to,sy,wmat[i][j][from][to][sy],pmat[i][j][from][to],delm[i][j][from][to],deld[i][j][from][to],synloc[from][to],d)
}
}
}
}
}
}
//** swirecut (col,EEsq,EIsq,IEsq,IIsq[,seed,slambda]) - spatial wiring: wires the column using pmat and cell positions
// (wiring probability effected by distance btwn cells)
// seed is random # seed
// lambda specifies length-constant for spatially-dependent fall-off in wiring probability
// at one lambda away, probability of connect is value in pmat
proc swirecut () { local x,y,z,ii,jj,a,del,prid,poid,prty,poty,dv,lseed,h,prob,slambda,dsq,dist,slambdasq,\
EEsq,EIsq,IEsq,IIsq,ic1,ic2,pdx\
localobj col,ce,vidx,vdel,vdist,vwt1,vwt2,vtmp,opr,opo,st,rdm,vprob
a=allocvecs(vidx,vdel,vtmp,vdist,vwt1,vwt2,vprob) z=0
col=$o1 ce=col.ce EEsq=$2 EIsq=$3 IEsq=$4 IIsq=$5
if (argtype(6)==0) lseed=$6 else lseed=1234
if(numarg()>6) slambda=$7 else slambda=10
if(slambda<=0){
printf("swirecut WARN: invalid slambda=%g, setting slambda to %g\n",colside/3)
slambda=colside/3
}
slambdasq = slambda^2 // using squared distance
vrsz(1e4,vidx,vdel,vdist,vtmp)
rdm=new Random() rdm.ACG(lseed) // initdivrnd(lseed)//initialize random # generator
rdm.uniform(0,1) vprob.resize(ce.count^2) vprob.setrand(rdm) pdx=0
// Create the connectivity NQS table.
if(col.mknetnqss) {nqsdel(col.connsnq) col.connsnq=new NQS("id1","id2","del","dist","wt1","wt2")}
for y=0,ce.count-1 { opr=ce.o(y)
vrsz(0,vidx,vdel,vdist,vwt1,vwt2)
prid=opr.id prty=opr.type ic1=ice(opr.type)
for poty=0,CTYPi-1 if (col.numc[poty]!=0 && (h=col.pmat[prty][poty])>0) {
for poid=col.ix[poty],col.ixe[poty] if(prid!=poid) { // go thru postsynaptic cells
opo = ce.o(poid) ic2=ice(ce.o(poid).type)
dx = opr.xloc - opo.xloc
dy = opr.yloc - opo.yloc
dsq = dx^2 + dy^2
if(ic1 && ic2 && dsq > IIsq) { // check max dists
continue
} else if(!ic1 && !ic2 && dsq > EEsq) {
continue
} else if(!ic1 && ic2 && dsq > EIsq) {
continue
} else if(ic1 && !ic2 && dsq > IEsq) continue
if(dsq <= slambdasq) {
prob=h
} else {
prob = h * exp(1 - sqrt(dsq)/slambda) // probability of connect
// print slambda,dx,dy,dsq,h,prob
}
if( prob >= vprob.x(pdx) ) { // rdm.uniform(0,1)
del = rdm.uniform(col.delm[prty][poty]-col.deld[prty][poty],col.delm[prty][poty]+col.deld[prty][poty])
{vidx.append(poid) vdel.append(del)}
if(synloc[prty][poty]==DEND) vdist.append(1) else vdist.append(0)
if(col.mknetnqss || wsetting_INTF6==1) {
if(col.syty1[prty][poty]>=0) vwt1.append(col.wmat[prty][poty][col.syty1[prty][poty]]) else vwt1.append(0)
if(col.syty2[prty][poty]>=0) vwt2.append(col.wmat[prty][poty][col.syty2[prty][poty]]) else vwt2.append(0)
}
}
pdx += 1
}
}
if(vidx.size>0) {
if(wsetting_INTF6==1) opr.setdvi(vidx,vdel,vdist,1,vwt1,vwt2) else opr.setdvi(vidx,vdel,vdist,1)
if(col.mknetnqss) {
for ii=0,vidx.size-1 col.connsnq.append(prid,vidx.x(ii),vdel.x(ii),vdist.x(ii),vwt1.x(ii),vwt2.x(ii))
}
}
}
dealloc(a)
if(verbose) printf("\n")
}
//* conprob - calculate connection probability using global slambda and hval
func conprob () {
if ($1<=slambda) return hval else return hval*exp(1-$1/slambda) // probability of connect
}
//** swirecutfl (col,EEsq,EIsq,IEsq,IIsq[,seed]) - spatial wiring: wires the column using pmat and cell positions
// (wiring probability effected by distance btwn cells)
// seed is random # seed
// the global slambda specifies length-constant for spatially-dependent fall-off in wiring probability
// at one slambda away, probability of connect is value in pmat - NB: slambda is a global!
proc swirecutfl () { local x,y,z,ii,jj,a,del,prid,poid,prty,poty,dv,lseed,h,\
EEsq,EIsq,IEsq,IIsq,ic1,ic2,pdx,n,nn,maxd,rsz\
localobj col,ce,vidx,vdel,vdist,vwt1,vwt2,vtmp,vpre,opr,opo,st,rdm,vprob,oq,ol,mmaxd
oq=new NQS("id1","id2","del","dist","wt1","wt2","dist0","prob") oq.verbose=0 ol=new List() mmaxd=new Matrix(2,2)
a=allocvecs(vidx,vdel,vtmp,vdist,vwt1,vwt2,vprob,vpre) z=0
{ol.append(vpre) ol.append(vidx) ol.append(vdel) ol.append(vdist) ol.append(vwt1) ol.append(vwt2)}
col=$o1 ce=col.ce EEsq=$2 EIsq=$3 IEsq=$4 IIsq=$5
if (argtype(6)==0) lseed=$6 else lseed=1234
vrsz(1e4,vidx,vdel,vdist,vtmp,vpre)
rdm=new Random() rdm.ACG(lseed) // initdivrnd(lseed)//initialize random # generator
rdm.uniform(0,1) rsz=ce.count^2 vprob.resize(rsz) vprob.setrand(rdm) pdx=0
{mmaxd.x(0,0)=sqrt(EEsq) mmaxd.x(0,1)=sqrt(EIsq) mmaxd.x(1,0)=sqrt(IEsq) mmaxd.x(1,1)=sqrt(IIsq)}
// Create the connectivity NQS table.
if(col.mknetnqss) {nqsdel(col.connsnq) col.connsnq=new NQS("id1","id2","del","dist","wt1","wt2")}
for y=0,ce.count-1 { opr=ce.o(y)
vrsz(0,vidx,vdel,vdist,vwt1,vwt2)
prid=opr.id prty=opr.type ic1=ice(opr.type)
for poty=0,CTYPi-1 if (col.numc[poty]!=0 && (h=hval=col.pmat[prty][poty])>0) {
ic2 = ice(poty)
oq.clear()
if ((n=opr.floc(opr.xloc,opr.yloc,1e9,oq.getcol("id2"),oq.getcol("dist"),mmaxd.x(ic1,ic2),poty))==0) {
if(verbose) printf("swirecutfl: No possibilities for connections found from %d\n",opr.id)
} else {
oq.getcol("dist0").copy(oq.getcol("dist"))
oq.pad(n)
oq.getcol("dist0").apply("conprob") // turn distances into probabilities
if (pdx+n>=rsz) { rdm.uniform(0,1) vprob.setrand(rdm) pdx=0 }
oq.getcol("prob").copy(vprob,pdx,pdx+n-1) // pick out some random probabilities to compare to
pdx+=n
oq.calc("<dist0>.sub(<prob>)")
if ((nn=oq.select("dist0",">",0,"id2","!=",prid))==0) { // looks for dist_prob > prob
if(verbose) printf("swirecutfl: No connections found from %d out of %d assayed\n",opr.id,n)
} else {
oq.cpout()
//oq.fill("id1",prid)
rdm.uniform(col.delm[prty][poty]-col.deld[prty][poty],col.delm[prty][poty]+col.deld[prty][poty])
oq.getcol("del").setrand(rdm)
if (col.syty1[prty][poty]>=0) oq.getcol("wt1").fill(col.wmat[prty][poty][col.syty1[prty][poty]]) else oq.getcol("wt1").fill(0)
if (col.syty2[prty][poty]>=0) oq.getcol("wt2").fill(col.wmat[prty][poty][col.syty2[prty][poty]]) else oq.getcol("wt2").fill(0)
{vwt1.append(oq.getcol("wt1")) vwt2.append(oq.getcol("wt2"))}
{vidx.append(oq.getcol("id2")) vdel.append(oq.getcol("del"))}
vdist.resize(vdist.size+nn)
if(synloc[prty][poty]==DEND) vdist.fill(1,vdist.size-nn,vdist.size-1) else vdist.fill(0,vdist.size-nn,vdist.size-1)
}
}
}
if(vidx.size>0) {
if(wsetting_INTF6==1) opr.setdvi(vidx,vdel,vdist,1,vwt1,vwt2) else opr.setdvi(vidx,vdel,vdist,1)
if(col.mknetnqss){vpre.resize(vidx.size) vpre.fill(prid) col.connsnq.append(ol)}
// for ii=0,vidx.size-1 col.connsnq.append(prid,vidx.x(ii),vdel.x(ii),vdist.x(ii),vwt1.x(ii),vwt2.x(ii))
} else print "empty div", y
}
col.ce.o(0).finishdvir()
dealloc(a)
nqsdel(oq)
if(verbose) printf("\n")
}
//* swireDPs(seed) - wire the DP (muscle proprioceptive) cells - they're only in COLUMN[0] (sensory)
// NB: this function relies on the global slambda, via the conprob function call
proc swireDPs () { local x,y,z,ii,jj,a,n,del,prid,poid,prty,poty,dv,lseed,prob,ic1,ic2,pdx,rsz,cdx\
localobj cl,lce,vidx,vdel,vdend,vwt1,vwt2,vprob,vtmp,opr,opo,st,nc,oq,ol,rdm,vpoty
//prid,poid,vdel.x[ii],vdist.x[ii],vwt1.x[ii],vwt2.x[ii]
oq=new NQS("id1","id2","del","dist","wt1","wt2","dist0","prob") oq.verbose=0
a=allocvecs(vprob,vpoty) ol=new List()
for ii=0,5 ol.append(oq.v[ii]) // will use ol for the appends
vidx=oq.getcol("id1") // vdel,vdend,vwt1,vwt2,vtmp,vprob
cdx=0
cl=col[cdx] lce=cl.ce
if (argtype(1)==0) lseed=$1 else lseed=1234
rdm=new Random() rdm.ACG(lseed) // initdivrnd(lseed)//initialize random # generator
rdm.uniform(0,1) vprob.resize(rsz=lce.count^2) vprob.setrand(rdm) pdx=0
vpoty.append(ES)
//** wire in the DPs
prty=DP
for vtr(&poty,vpoty) {
hval=pmat[cdx][cdx][prty][poty]
for ltr(opr,cedp,&ii) { prid=opr.id
oq.clear()
if ((n=col[cdx].ce.o(0).floc(opr.xloc,opr.yloc,1e9,oq.getcol("id2"),oq.getcol("dist"),4*slambda,poty))==0) {
printf("swireDPs: No possibilities for connections found from %d\n",opr.id)
} else {
oq.getcol("dist0").copy(oq.getcol("dist"))
oq.pad(n)
oq.getcol("dist0").apply("conprob") // turn distances into probabilities
if (pdx+n>=rsz) { vprob.setrand(rdm) pdx=0 }
oq.getcol("prob").copy(vprob,pdx,pdx+n-1) // pick out some random probabilities to compare to
pdx+=n
oq.calc("<dist0>.sub(<prob>)")
if ((nn=oq.select("dist0",">",0))==0) { // looks for dist_prob > prob
printf("swireDPs: No connections found from %d out of %d assayed\n",opr.id,n)
} else {
oq.cpout()
oq.fill("id1",prid)
rdm.uniform(cl.delm[prty][poty]-cl.deld[prty][poty],cl.delm[prty][poty]+cl.deld[prty][poty])
oq.getcol("del").setrand(rdm)
if (cl.syty1[prty][poty]>=0) oq.getcol("wt1").fill(cl.wmat[prty][poty][cl.syty1[prty][poty]]) else oq.getcol("wt1").fill(0)
if (cl.syty2[prty][poty]>=0) oq.getcol("wt2").fill(cl.wmat[prty][poty][cl.syty2[prty][poty]]) else oq.getcol("wt2").fill(0)
if (cl.mknetnqss) cl.connsnq.append(ol)
for vtr(&poid,oq.out.getcol("id2"),&jj) {
ncl.append(nc=new NetCon(opr.drspk,lce.o(poid)))
nc.delay=oq.getcol("del").x[jj]
nc.weight[cl.syty1[prty][poty]]=oq.getcol("wt1").x[jj]
nc.weight[cl.syty2[prty][poty]]=oq.getcol("wt2").x[jj]
ncdp.append(nc) // save to ncdp too
}
}
}
}
}
dealloc(a)
nqsdel(oq)
}
//* wireDPs(seed) - wire the DP (muscle proprioceptive) cells - they're only in COLUMN[0] (sensory)
// this uses fixed convergence
proc wireDPs () { local x,y,z,ii,jj,a,n,del,prid,poid,prty,poty,dv,lseed,cdx,cv,cnt\
localobj cl,lce,vidx,vdist,vdel,vwt1,vwt2,opr,opo,st,nc,oq,ol,rdm,vpoty,v1,v2,v3,v4,v5,v6,v7,l
a=allocvecs(vidx,vdel,vdist,vpoty,vwt1,vwt2,v1,v2,v3,v4,v5,v6,v7) ol=new List() l=new List() l.append(v1) l.append(v4)
cdx=0
cl=col[cdx] lce=cl.ce
if (argtype(1)==0) lseed=$1 else lseed=1234
rdm=new Random() rdm.ACG(lseed) // initdivrnd(lseed)//initialize random # generator
rdm.uniform(0,1)
vpoty.append(ES)
//** wire in the DPs
prty=DP
for vtr(&poty,vpoty) {
if((cv = int(.5+cl.pmat[prty][poty]*cl.numc[prty])) < 1) continue // convergence
for poid = cl.ix[poty],cl.ixe[poty] {
vrsz(0,vidx,vdel,vdist,vwt1,vwt2)
vrsz(MAXxy(2*cv,4*cl.numc[prty]),v1,v2,v3)
{rdm.discunif(0,cl.numc[prty]-1) v3.setrand(rdm)}
if(prty==poty) {v1.cull(v3,prid) v3.copy(v1) v2.resize(v3.size)} // get rid of self connect
if( (cnt=v3.uniq(l,1)) > cv){ cnt = cv } // puts unsorted uniq vals into v4
vrsz(cnt,v2,v4,v5,v6,v7) // v4 has prids
rdm.uniform(cl.delm[prty][poty]-cl.deld[prty][poty],cl.delm[prty][poty]+cl.deld[prty][poty])
v2.setrand(rdm)
v2.abs() // no negative delays
if(cl.synloc[prty][poty]==DEND) v5.fill(1) else v5.fill(0)
vidx.append(v4) vdel.append(v2) vdist.append(v5) // vidx.append(v1) will give sorted bug
if(cl.syty1[prty][poty]>=0) v6.fill(cl.wmat[prty][poty][cl.syty1[prty][poty]]) else v6.fill(0)
if(cl.syty2[prty][poty]>=0) v7.fill(cl.wmat[prty][poty][cl.syty2[prty][poty]]) else v7.fill(0)
vwt1.append(v6) vwt2.append(v7)
if (cl.mknetnqss) {
for ii=0,vidx.size-1 cl.connsnq.append(vidx.x(ii)+cl.ix[prty],poid,vdel.x(ii),vdist.x(ii),vwt1.x(ii),vwt2.x(ii))
}
for vtr(&prid,vidx,&jj) { opr=cedp.o(prid)
ncl.append(nc=new NetCon(opr.drspk,lce.o(poid)))
nc.delay=vdel.x(jj)
nc.weight[cl.syty1[prty][poty]]=vwt1.x(jj)
nc.weight[cl.syty2[prty][poty]]=vwt2.x(jj)
ncdp.append(nc) // save to ncdp too
}
}
}
dealloc(a)
}
//* checkerpos(col,rad,xinc,yinc,seed) - arrange cells in checkerboard pattern
proc checkerpos () { local i,j,x,y,csz,xinc,yinc,seed,rad,a localobj col,ce,rdm,xo,nqc
col=$o1 csz=col.cside rad=$2 xinc=$3 yinc=$4 rdm=new Random()
if(numarg()>4) seed=$5 else seed=1234 rdm.ACG(seed) ce=col.ce nqc=new NQS("x","y")
for(y=0;y<csz;y+=yinc) {
if(y%2==0) {
for(x=0;x<csz;x+=xinc) nqc.append(x,y)
} else {
for(x=xinc/2;x<csz;x+=xinc) nqc.append(x,y)
}
}
if(verbose) {gvmarkflag=1 g.erase nqc.gr("y","x")}
for ltr(xo,ce,&i) {
r = rdm.discunif(0,nqc.v.size-1)
x = nqc.v[0].x(r) + rdm.uniform(-rad,rad)
y = nqc.v[1].x(r) + rdm.uniform(-rad,rad)
while(x<0 || y<0 || x>=csz || y>=csz) {
r = rdm.uniform(0,nqc.v.size-1)
x = nqc.v[0].x(r) + rdm.uniform(-rad,rad)
y = nqc.v[1].x(r) + rdm.uniform(-rad,rad)
}
{xo.xloc=x xo.yloc=y}
if(col.cellsnq!=nil) {col.cellsnq.v[5].x(i)=xo.xloc col.cellsnq.v[6].x(i)=xo.yloc}
}
nqsdel(nqc)
}
//* closestsqrt(number) - return the closest integer square root
func closestsqrt () { local n,s,si,sa,sb
n = $1
s = sqrt(n)
si = int(s)
sa = si*si
sb = si*(si+2)+1
if (sb-n < n-sa) return sqrt(sb) else return si
}
//* setgridpos(col,cside) - arrange cells on a 2D grid
proc setgridpos () { local done,odd,i,j,ct,x,y,csz,xinc,yinc,a localobj col,ce,xo,nqc
col=$o1 csz=$2
ce=col.ce nqc=new NQS("id","x","y")
for ct=0,CTYPi-1 if(ct!=DP && col.numc[ct]) {
xinc = yinc = csz / (closestsqrt(col.numc[ct]) - 1)
x = y = odd = done = 0
for i=col.ix[ct],col.ixe[ct] {
xo = ce.o(i)
xo.xloc = x
xo.yloc = y
if(verbose) nqc.append(xo.id,x,y)
if(col.cellsnq!=nil) {col.cellsnq.v[5].x(i)=xo.xloc col.cellsnq.v[6].x(i)=xo.yloc}
if(!done) x += xinc
if(x > csz ) {
odd = !odd
if(odd) x = xinc / 2 else x = 0
if(!done) y += yinc
if(y > csz) {
y = csz / 2
x = csz / 2
done = 1
}
}
}
}
if(verbose) {gvmarkflag=1 g.erase nqc.gr("y","x")}
nqsdel(nqc)
}
//* convwire(col,[sead]) - wire the column with fixed convergence
proc convwire () { local x,y,z,ii,jj,a,del,prid,prty,poid,poty,cv,cvt,dv,dvt,lseed,dvrand,szo,mkt\
localobj v1,v2,v3,v4,v5,v6,v7,vidx,vdel,vdist,vwt1,vwt2,vtmp,opr,opo,st,l,rdm,col,ce
a=allocvecs(v1,v2,v3,v4,v5,v6,v7,vidx,vdel,vtmp,vdist,vwt1,vwt2) z=0 l=new List() l.append(v1) l.append(v4)
col=$o1 ce=col.ce
if (numarg()>1) lseed=$2 else lseed=1234
vrsz(1e4,vidx,vdel,vdist,vtmp)
{rdm=new Random() rdm.MCellRan4(lseed,lseed)}
// Create the connectivity NQS table.
nqsdel(col.connsnq)
col.connsnq=new NQS("id1","id2","del","dist","wt1","wt2")
for y=0,ce.count-1 { opo=ce.o(y) // go thru the postsynaptic cells
vrsz(0,vidx,vdel,vdist,vwt1,vwt2)
poid=opo.id
poty=opo.type
for prty=0,CTYPi-1 if (col.numc[prty]!=0 && (cv=int(col.conv[prty][poty]))>0){//go thru presynaptic types
vrsz(MAXxy(2*cv,4*col.numc[prty]),v1,v2,v3)
{rdm.discunif(col.ix[prty],col.ixe[prty]) v3.setrand(rdm)}
// 2023-02-28: the use of an unitialised value of prid was causing inconsistent results between different versions of NEURON
prid = 0
if(prty==poty) {v1.cull(v3,prid) v3.copy(v1) v2.resize(v3.size)} // get rid of self connect -- NB: prid wasn't set!! bug--fixit (poid==prid anyway)
if( (cnt=v3.uniq(l,1)) > cv){ cnt = cv } // puts unsorted uniq vals into v4
vrsz(cnt,v2,v4,v5,v6,v7) // v4 has prids
rdm.uniform(col.delm[prty][poty]-col.deld[prty][poty],col.delm[prty][poty]+col.deld[prty][poty])
v2.setrand(rdm)
v2.abs() // no negative delays
if(col.synloc[prty][poty]==DEND) v5.fill(1) else v5.fill(0)
vidx.append(v4) vdel.append(v2) vdist.append(v5) // vidx.append(v1) will give sorted bug
if(col.syty1[prty][poty]>=0) v6.fill(col.wmat[prty][poty][col.syty1[prty][poty]]) else v6.fill(0)
if(col.syty2[prty][poty]>=0) v7.fill(col.wmat[prty][poty][col.syty2[prty][poty]]) else v7.fill(0)
vwt1.append(v6) vwt2.append(v7)
}//end prty
if(vidx.size>0) {
for ii=0,vidx.size-1 col.connsnq.append(vidx.x(ii),poid,vdel.x(ii),vdist.x(ii),vwt1.x(ii),vwt2.x(ii))
}
}//end cell loop
mkt=col.mknetnqss
col.mknetnqss=0
col.wirenq(col.connsnq) // wire from the NQS just created
col.mknetnqss=mkt
dealloc(a)
if(verbose) printf("\n")
}
//* mkcols - make the COLUMNs
proc mkcols () { local id,x,y,seed,svnum localobj nq
id=nextGID_INTF6=0
for y=0,colH-1 for x=0,colW-1 {
if(dbgcols)seed=dvseed else seed=(id+1)*dvseed
svnum=vcpercol[id].x[DP] vcpercol[id].x[DP]=0 // ==== REMOVE DPs so can create separately ====
lcol.append(gcol[y][x]=new COLUMN(id,vcpercol[id],colnq[id],seed,x,y,setdviPT,mknetnqss,1))
print "made COLUMN, wiring"
// vcpercol[id].x[DP]=svnum // replace
col[id]=gcol[y][x]
col[id].verbose=verbose_INTF6
nextGID_INTF6 += cpercol[0][DP] // inc by # of DPs
if(strlen(wnqstr)>0 && id==1) { // wire from NQS
nq=new NQS(wnqstr)
col[id].wirenq(nq)
nqsdel(nq)
} else if(wirety>0) {
col[id].setcellpos(vlayerz,layerzvar,pseed*(id+1),colside) // cell positions used for spatial-wiring
if(checkers) checkerpos(col[id],crad,cxinc,cyinc,pseed*(id+1)) // arrange cells on checker-like grid
if(gridpos) setgridpos(col[id],colside) // arrange cells on 2D grid
if(wirety==1) {
col[id].swire(col[id].wseed,maxsfall)
} else if(wirety==2) {
swirecut(col[id],EEsq,EIsq,IEsq,IIsq,col[id].wseed,slambda)
} else if(wirety==3) {
swirecutfl(col[id],EEsq,EIsq,IEsq,IIsq,col[id].wseed) // uses intf6::floc for speedup
} else if(wirety==4) {
convwire(col[id],col[id].wseed) // random convergence-based wiring (no spatial dependence)
}
} else col[id].wire(col[id].wseed) // random divergence-based wiring (no spatial dependence)
id+=1
}
}
//* hexgridnq(width,height,number of elements) - makes an nqs with hexagonal grid arrangement
// tries to place 'cells' w. equidistant separation. if fills up, moves to center.
obfunc hexgridnq () { local done,odd,i,w,h,sz,n,xinc,yinc,x,y localobj nq
w=$1 h=$2 n=$3
nq=new NQS("i","x","y")
xinc = w / (closestsqrt(n) - 1)
yinc = h / (closestsqrt(n) - 1)
x = y = odd = done = 0
for i=0,n-1 {
nq.append(i,x,y)
if(!done) x += xinc
if(x > w ) {
odd = !odd
if(odd) x = xinc / 2 else x = 0
if(!done) y += yinc
if(y > h) {
y = h / 2
x = w / 2
done = 1
}
}
}
return nq
}
//** mkDPs() create and locate the Stimulators -- these are DRSPK cells rather than INTF6
proc mkDPs () { local ii,x,y,c,ndp,suty,rw,cl,rs,cs,pad,cdx,xinc,yinc\
localobj cellsnq,xo,rdm,ce,nqh
cdx=0
ndp=col[cdx].numc[DP]=cpercol[cdx][DP] // will be outside of col now
col[cdx].ix[DP]=col[cdx].allcells col[cdx].ixe[DP]=col[cdx].allcells+ndp-1
ce=col[cdx].ce gid=0
if(NMUSCLES==4) rs=cs=2 else {rs=1 cs=2} // # rows and columns for DP spatial layout
for ii=0,col[cdx].numc[DP]-1 cedp.append(new DPC()) // create the DPs
rdm=new Random() rdm.ACG(pseed+DP)
if(DPpad) pad=colside/cs/8 else pad=0
if (col[cdx].cellsnq==nil) c=-1 else {
cellsnq=col[cdx].cellsnq cellsnq.tog("DB") c=cellsnq.fi("xloc")
gid=cellsnq.applf(".max","gid")+1
if(gid!=col[cdx].ix[DP]) {printf("mkDPS ERRA %d\n",gid) return }
}
if (gid==0) gid=col[cdx].ix[DP]
//** assign DP identity
if(DPgridpos) {
if(sepDPpos) {
nqh=hexgridnq(colside/2,colside/2,cedp.count()/NMUSCLES)
} else nqh=hexgridnq(colside,colside,cedp.count())
for ltr(xo,cedp,&ii) {
xo.id=gid xo.type=DP suty=xo.subtype=int(xo.id%NMUSCLES) // which muscle it belongs to
pad = 0
rw=int(suty/cs) cl=suty%cs // xloc,yloc
// print suty,rw,cl,cl*colside/cs,(cl+1)*colside/cs,rw*colside/cs,(rw+1)*colside/rs
if(sepDPpos) { // each set of DP cells is spatially separated
xo.xloc=x= cl*colside/2 + nqh.getcol("x").x((xo.id-col.ix[DP])/NMUSCLES)
xo.yloc=y= rw*colside/2 + nqh.getcol("y").x((xo.id-col.ix[DP])/NMUSCLES)
} else { // dont segregate DP cells
xo.xloc=x=nqh.getcol("x").x(ii)
xo.yloc=y=nqh.getcol("y").x(ii)
}
xo.zloc=z=suty
ce.o(col[cdx].ix[ES]+ii).zloc=ce.o(col[cdx].ix[EM]+ii).zloc=z//same group
// add DPs to NQS; gid(0) id(1) col(2) ty(3) ice(4) xloc(5) yloc(6) zloc(7)
// gid id col ty ice x,y,z
if(c!=-1) col[cdx].cellsnq.append(gid,xo.id,0,xo.type,0,x,y,suty) // use zloc to store muscle #
gid+=1
}
nqsdel(nqh)
} else for ltr(xo,cedp,&ii) {
xo.id=gid xo.type=DP suty=xo.subtype=int(xo.id%NMUSCLES) // which muscle it belongs to
rw=int(suty/cs) cl=suty%cs // xloc,yloc
// print suty,rw,cl,cl*colside/cs,(cl+1)*colside/cs,rw*colside/cs,(rw+1)*colside/rs
if(sepDPpos) { // each set of DP cells is spatially separated
xo.xloc=x=rdm.uniform(cl*colside/cs+pad,(cl+1)*colside/cs-pad)
xo.yloc=y=rdm.uniform(rw*colside/rs+pad,(rw+1)*colside/rs-pad)
} else { // dont segregate DP cells
xo.xloc=x=rdm.uniform(0,colside)
xo.yloc=y=rdm.uniform(0,colside)
}
xo.zloc=z=suty
ce.o(col[cdx].ix[ES]+ii).zloc=ce.o(col[cdx].ix[EM]+ii).zloc=z//same group
// add DPs to NQS; gid(0) id(1) col(2) ty(3) ice(4) xloc(5) yloc(6) zloc(7)
// gid id col ty ice x,y,z
if(c!=-1) col[cdx].cellsnq.append(gid,xo.id,0,xo.type,0,x,y,suty) // use zloc to store muscle #
gid+=1
}
}
//* wirecols - setup inter-COLUMN connectivity with NetCon
proc wirecols () { local i,j
if(numcols==1) return
for case(&i,0,1) if(EEGainInterC[i] || EIGainInterC[i]) {
j = 1 - i
gcol[0][i].wire2col(gcol[0][j],colnq[i],1,ncl) // unidirectional wiring
}
}
//* intercoloff(c0,c1) - turn off all weights between COLUMNs
proc intercoloff () { local i,c0,c1 localobj xo
c0=$1 c1=$2
for ltr(xo,ncl) if(isojt(xo.pre,col.ce.o(0)) && isojt(xo.syn,col.ce.o(0))) {
if(xo.pre.col!=c0 || xo.syn.col!=c1) continue
for i=0,6 xo.weight(i)=0
}
}
//* intercolmul(c0,c1,from,to,sy,w)
proc intercolsyw () { local c0,c1,from,to,sy,w localobj xo
c0=$1 c1=$2 from=$3 to=$4 sy=$5 w=$6
for ltr(xo,ncl) if(isojt(xo.pre,col.ce.o(0)) && isojt(xo.syn,col.ce.o(0))) {
if(xo.pre.col!=c0 || xo.syn.col!=c1) continue
if(xo.pre.type==from && xo.syn.type==to) xo.weight(sy)=w
}
}
//** drdp() DP locations
proc drdp () { local cdx localobj nq
nq=col[0].cellsnq gvmarkflag=1 ge(0)
nq.marksym="o"
for ii=0,5 if (nq.select("ty",DP,"zloc",ii)) nq.gr("yloc","xloc",0,ii+1,8)
}
//* drcelllocs - draw location of cells
proc drcelllocs () { local i,gvt,ty
gvt=gvmarkflag
{gvmarkflag=1 col[CDX].cellsnq.marksym="O"}
if(numarg()>0) {
if(col[CDX].cellsnq.select("ty",$1)) col[CDX].cellsnq.gr("yloc","xloc",0,1,3)
} else for i=0,1 if(col[CDX].cellsnq.select("ice",i)) col[CDX].cellsnq.gr("yloc","xloc",0,i+2,3)
gvmarkflag=gvt
}
//* drcelldiv(cellid[,skipI,skipE]) - draw outputs from a cell
proc drcelldiv () { local id1,id2,x1,y1,x2,y2,ic1,ic2,clr1,clr2,gvt,skipI,skipE localobj cnq,nq
gvt=gvmarkflag gvmarkflag=1
cnq=col[CDX].connsnq
nq=col[CDX].cellsnq
id1=$1 nq.verbose=0
if(!nq.select(-1,"id",id1)) return
if(numarg()>1) skipI=$2 else skipI=0
if(numarg()>2) skipE=$3 else skipE=0
x1=nq.v[5].x(nq.ind.x(0))
y1=nq.v[6].x(nq.ind.x(0))
ic1=nq.v[4].x(nq.ind.x(0))
if(ic1) clr1=3 else clr1=2
nq.select("id",id1)
nq.gr("yloc","xloc",0,clr1,15)
cnq.select(-1,"id1",id1)
cnq.verbose=0
for vtr(&i,cnq.ind) {
id2 = cnq.v[1].x(i)
if(nq.select("id",id2)) {
x2=nq.fetch("xloc")
y2=nq.fetch("yloc")
ic2=nq.fetch("ice")
if(ic2) {
clr2=3
if(skipI) continue
} else {
clr2=2
if(skipE) continue
}
nq.gr("yloc","xloc",0,clr2,8)
drline(x1,y1,x2,y2,g,clr1,7)
}
}
cnq.verbose=nq.verbose=1
gvmarkflag=gvt
}
//* drcellconv(cellid) - draw inputs to a cell
proc drcellconv () { local id1,id2,x1,y1,x2,y2,ic1,ic2,clr1,clr2,gvt localobj cnq,nq
gvt=gvmarkflag gvmarkflag=1
cnq=col[CDX].connsnq
nq=col[CDX].cellsnq
id2=$1 nq.verbose=0
nq.select(-1,"id",id2)
x2=nq.v[5].x(nq.ind.x(0))
y2=nq.v[6].x(nq.ind.x(0))
ic2=nq.v[4].x(nq.ind.x(0))
if(ic2) clr2=3 else clr2=2
nq.select("id",id2)
nq.gr("yloc","xloc",0,clr2,15)
cnq.select(-1,"id2",id2)
cnq.verbose=0
for vtr(&i,cnq.ind) {
id1 = cnq.v[0].x(i)
if(nq.select("id",id1)) {
x1=nq.fetch("xloc")
y1=nq.fetch("yloc")
ic1=nq.fetch("ice")
if(ic1) clr1=3 else clr1=2
nq.gr("yloc","xloc",0,clr1,8)
drline(x1,y1,x2,y2,g,clr1,7)
}
}
cnq.verbose=nq.verbose=1
gvmarkflag=gvt
}
//* mkwgnq() - return nqs with intracolumnar weight info
obfunc mkwgnq () { local a,i,j,id1,id2,wg,ic,ty1,ty2,wt1,wt2,del,dist,c\
localobj mycel,vwt1,vwt2,myv,myv1,myv2,myv3,myv4,myv5,myv6,wgnq,col
a=allocvecs(vwt1,vwt2,myv,myv1,myv2,myv3,myv4,myv5,myv6)
{nqsdel(wgnq) wgnq=new NQS("col","id1","id2","wg","ice","ty1","ty2","wt1","wt2","del","dist")}
for ltr(col,lcol,&c) for i=0,col.allcells-1 {
mycel=col.ce.o(i)
mycel.getdvi(myv,myv1,myv2,myv3,myv4,myv5,myv6)
id1=i
ic=ice(mycel.type)
ty1=mycel.type
vrsz(mycel.getdvi,vwt1,vwt2)
if(wsetting_INTF6==1) {
mycel.getsywv(vwt1,vwt2)
} else {
vwt1.copy(myv3)
vwt2.copy(myv4)
}
for j=0,myv6.size-1 {
id2=myv.x(j)
ty2=col.ce.o(id2).type
wg = myv6.x(j)
wt1 = wg * vwt1.x(j)
wt2 = vwt2.x(j)
del = myv1.x(j)
dist = myv5.x(j)
wgnq.append(c,id1,id2,wg,ic,ty1,ty2,wt1,wt2,del,dist)
}
}
dealloc(a)
return wgnq
}
//* getavgwgnq(wgnq) - get avg/stderr weight gains btwn populations -
// wgnq is nqs from mkwgnq function
obfunc getavgwgnq () { local i,ty1,ty2,a localobj vpre,vpo,nqw,nqa
nqw=$o1
a=allocvecs(vpre,vpo)
vpre.append(ES,ES, ES,ES,EM,EM,EM, EM)
vpo.append(ES,IS,ISL,EM,ES,EM,IM,IML)
nqa=new NQS("ty1","ty2","wgavg","wgstderr")
nqw.verbose=0
for i=0,vpre.size-1 {
ty1=vpre.x(i)
ty2=vpo.x(i)
if(nqw.select("ty1",ty1,"ty2",ty2)) {
nqa.append(ty1,ty2,nqw.getcol("wg").mean(),nqw.getcol("wg").stderr())
}
}
nqw.verbose=1
dealloc(a)
return nqa
}
//* setNCWeight(NetCon List, synapse, weight) - sets the weights of the NetCons to specified value
proc setNCWeight () { local i,sy,w localobj nc,ncl
ncl=$o1 sy=$2 w=$3
for ltr(nc,ncl,&i) nc.weight[sy]=w
}
//* setDPESWeightAM2(weight) - set weight of DP -> ES AM2 connections
proc setDPESWeightAM2 () {
setNCWeight(ncdp,AM2,$1)
}
//* setDPESWeightNM2(weight) - set weight of DP -> ES NM2 connections
proc setDPESWeightNM2 () {
setNCWeight(ncdp,NM2,$1)
}
//* function calls to setup network
//** z location of each layer in microns
mkvlayerz()
//** # of cells per column
setcpercol()
//** setup pmat
if(name_declared("nqpmat")==2) { // read pmat from NQS if available, else set to default
if(nqpmat!=nil) nq2pmat(nqpmat) else setpmat()
} else setpmat()
if(pmatscale!=1) scalepmat(pmatscale)
print "pmat[ES][IML]=",pmat[0][0][ES][IML],"ESIMFeedForward=",ESIMFeedForward
//** setup synapse locations,delays,wmat
setsynloc()
setdelmats()
setwmat() // new KMJ version
if(wmatscale!=1) scalewmat(wmatscale)
scrsz=50*1e3
double scr[scrsz]
//** make cells, columns, wire columns
{mkcolnqs() printf("finished mkcolnqs\n")}
{mkcols() printf("finished mkcols\n")}
wirecols(1) // only 1 col at present
{mkDPs() printf("finished mkDPs\n")}
if(wirety==4) {wireDPs(col[0].wseed) printf("finished wireDPs\n")} else{swireDPs(col[0].wseed) printf("finished swireDPs\n")}