Reinforcement learning of targeted movement (Chadderdon et al. 2012)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:144538
"Sensorimotor control has traditionally been considered from a control theory perspective, without relation to neurobiology. In contrast, here we utilized a spiking-neuron model of motor cortex and trained it to perform a simple movement task, which consisted of rotating a single-joint “forearm” to a target. Learning was based on a reinforcement mechanism analogous to that of the dopamine system. This provided a global reward or punishment signal in response to decreasing or increasing distance from hand to target, respectively. Output was partially driven by Poisson motor babbling, creating stochastic movements that could then be shaped by learning. The virtual forearm consisted of a single segment rotated around an elbow joint, controlled by flexor and extensor muscles. ..."
Reference:
1 . Chadderdon GL, Neymotin SA, Kerr CC, Lytton WW (2012) Reinforcement learning of targeted movement in a spiking neuronal model of motor cortex. PLoS One 7:e47251 [PubMed]
Citations  Citation Browser
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 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): Dopamine; Gaba; Glutamate;
Simulation Environment: NEURON;
Model Concept(s): Simplified Models; Synaptic Plasticity; Long-term Synaptic Plasticity; Reinforcement Learning; Reward-modulated STDP;
Implementer(s): Neymotin, Sam [Samuel.Neymotin at nki.rfmh.org]; Chadderdon, George [gchadder3 at gmail.com];
Search NeuronDB for information about:  GabaA; AMPA; NMDA; Dopamine; Gaba; Glutamate;
/
arm1d
README
drspk.mod *
infot.mod *
intf6_.mod *
intfsw.mod *
misc.mod *
nstim.mod *
stats.mod *
updown.mod *
vecst.mod *
arm.hoc
basestdp.hoc
col.hoc *
colors.hoc *
declist.hoc *
decmat.hoc *
decnqs.hoc *
decvec.hoc *
default.hoc *
drline.hoc *
filtutils.hoc *
geom.hoc
grvec.hoc *
hinton.hoc *
infot.hoc *
init.hoc
intfsw.hoc *
labels.hoc *
local.hoc *
misc.h *
mosinit.hoc
network.hoc
nload.hoc
nqs.hoc *
nqsnet.hoc *
nrnoc.hoc *
params.hoc
run.hoc
samutils.hoc *
sense.hoc *
setup.hoc *
sim.hoc
simctrl.hoc *
stats.hoc *
stim.hoc
syncode.hoc *
units.hoc *
xgetargs.hoc *
                            
// $Id: network.hoc,v 1.352 2012/03/08 17:56:22 samn Exp $

//* Numbers and connectivity params

DOPE_INTF6=1
EDOPE_INTF6=1
ESTDP_INTF6 = ISTDP_INTF6 = 0

declare("scale",2)  // Sam had this on 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

//* swire & related column variables
declare("colside",30) // column diameter in micrometers
declare("swire",0) // whether to use 'spatial' wiring, 2 means use swirecut, 3 means use swirecutfl  (Sam had set to 3)
declare("checkers",0) // whether to arrange cells in checkerboard pattern
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("lambda",30,"hval",0)

//* other variables
declare("mknetnqss",1) // whether COLUMN should make connsnq,cellsnq
declare("fc",0.55)
declare("EEGain",fc*4*15/11,"EIGain",fc*15,"IEGain",fc*4*15/11,"IIGain",fc*4*15/11)
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",4) // number of muscle groups in model
declare("sepDPpos",1) // whether to position DP cells in each group separately

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",1) // 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.15) // level of feedforward connectivity from DP -> ES
// declare("DPESFeedForward",0.1) // level of feedforward connectivity from DP -> ES
declare("EMRecur",0.0) //declare("EMRecur",0.07) // level of recurrence between EM cells
declare("ESRecur",0.0)//declare("ESRecur",0.191)// level of recurrence between ES cells
declare("EMESFeedBack",0.0)//declare("EMESFeedBack",0.017)// level of feedback from EM -> ES
declare("ESEMFeedForward",0.08)// level of feedforward connectivity from ES -> EM

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 

    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][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][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
}

//* 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
      }
    }
  }
}

//* 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]= 1.75 * 5.01 / EEGain // 1.67 

  // sensory layer
  wmat[ii][ii][ES][ES][AM2]=0.66
  wmat[ii][ii][ES][ISL][AM2]=0.23
  wmat[ii][ii][ES][IS][AM2]=0.23
  wmat[ii][ii][ES][EM][AM2]=0.88 * 2

  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.71
  wmat[ii][ii][EM][IML][AM2]=0.23
  wmat[ii][ii][EM][IM][AM2]=0.23
  wmat[ii][ii][EM][ES][AM2]=0.24

  wmat[ii][ii][IM][EM][GA]=1.5
  wmat[ii][ii][IM][IM][GA]=1.5
  wmat[ii][ii][IM][IML][GA]=1.5

  wmat[ii][ii][IML][EM][GA2]=0.83
  wmat[ii][ii][IML][IM][GA2]=1.5
  wmat[ii][ii][IML][IML][GA2]=1.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 *= 0.5
      } 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
func conprob () { // calculate connection probability using global slambda and hval
  if ($1<=slambda) return hval else return hval*exp(1-$1/slambda) // probability of connect
}

//** swirecutfl (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
// slambda specifies length-constant for spatially-dependent fall-off in wiring probability
// at one slambda away, probability of connect is value in pmat
proc swirecutfl () { 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,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  
  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,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))
    }
  }
  col.ce.o(0).finishdvir()
  dealloc(a)
  nqsdel(oq)
  if(verbose) printf("\n")
}

//* swireDPs(seed,slambda,) - wire the DP (muscle proprioceptive) cells - they're only in COLUMN[0] (sensory)
proc swireDPs () { local x,y,z,ii,jj,a,n,del,prid,poid,prty,poty,dv,lseed,prob,dsq,dist,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  
  if(argtype(2)==0) slambda=$2 else slambda=10
  if(slambda<=0){slambda=colside/3 printf("swireDPs WARN: invalid lambda=%g, setting lambda to %g\n",slambda)}
  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]
          }
        }
      }
    }
  }
  dealloc(a)
  nqsdel(oq)
}

//* 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)
}

//* 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(swire>0) { // spatial-wiring
      col[id].setcellpos(vlayerz,layerzvar,pseed*(id+1),colside)
      if(checkers) checkerpos(col[id],crad,cxinc,cyinc,pseed*(id+1))
      if(swire==1) {
        col[id].swire(col[id].wseed,maxsfall)
      } else if(swire==2) {
        swirecut(col[id],EEsq,EIsq,IEsq,IIsq,col[id].wseed,slambda)
      } else if(swire==3) {
        swirecutfl(col[id],EEsq,EIsq,IEsq,IIsq,col[id].wseed,slambda) // uses intf6::floc for speedup
      }
    } else col[id].wire(col[id].wseed) // random wiring (no spatial dependence)
    id+=1
  }
}

//** 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 localobj cellsnq,xo,rdm,ce
  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 rs=2 cs=2 // 2 rows and 2 columns
  for ii=0,col[cdx].numc[DP]-1 cedp.append(new DPC()) // create the DPs
  rdm=new Random() rdm.ACG(pseed+DP)
  pad=colside/cs/8
  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
  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)
      if(ic) wt2 = vwt2.x(j) else wt2 = wg * 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
}

//* 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)

//** 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")}
{swireDPs(col[0].wseed,lambda) printf("finished swireDPs\n")}