Motor system model with reinforcement learning drives virtual arm (Dura-Bernal et al 2017)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:194897
"We implemented a model of the motor system with the following components: dorsal premotor cortex (PMd), primary motor cortex (M1), spinal cord and musculoskeletal arm (Figure 1). PMd modulated M1 to select the target to reach, M1 excited the descending spinal cord neurons that drove the arm muscles, and received arm proprioceptive feedback (information about the arm position) via the ascending spinal cord neurons. The large-scale model of M1 consisted of 6,208 spiking Izhikevich model neurons [37] of four types: regular-firing and bursting pyramidal neurons, and fast-spiking and low-threshold-spiking interneurons. These were distributed across cortical layers 2/3, 5A, 5B and 6, with cell properties, proportions, locations, connectivity, weights and delays drawn primarily from mammalian experimental data [38], [39], and described in detail in previous work [29]. The network included 486,491 connections, with synapses modeling properties of four different receptors ..."
Reference:
1 . Dura-Bernal S, Neymotin SA, Kerr CC, Sivagnanam S, Majumdar A, Francis JT, Lytton WW (2017) Evolutionary algorithm optimization of biological learning parameters in a biomimetic neuroprosthesis. IBM Journal of Research and Development (Computational Neuroscience special issue) 61(2/3):6:1-6:14
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism:
Cell Type(s): Abstract Izhikevich neuron;
Channel(s):
Gap Junctions:
Receptor(s): GabaA; GabaB; NMDA; AMPA;
Gene(s):
Transmitter(s): Glutamate; Gaba;
Simulation Environment: NEURON; Python;
Model Concept(s): Learning; Reinforcement Learning; Reward-modulated STDP; STDP; Motor control; Sensory processing;
Implementer(s): Dura-Bernal, Salvador [salvadordura at gmail.com]; Kerr, Cliff [cliffk at neurosim.downstate.edu];
Search NeuronDB for information about:  GabaA; GabaB; AMPA; NMDA; Gaba; Glutamate;
// $Id: msarm.hoc,v 1.252 2012/07/21 21:55:29 samn Exp $ 
// Last revision: 5/25/13 (sald)
// Modification of the code for simple virtual 2D arm to interface with musculoskeletal 2D arm

// Set up Python object to allow interface with arminterface.py module.
objref pyobj
pyobj = new PythonObject()

objref gsyn // for graph with synaptic changes
objref ls // for power spectra

declare("verbosearm",0)

//* templates - need to write
begintemplate p2d
double x[1],y[1]
public x,y
proc init () {
  if(numarg()==2) {x=$1 y=$2}
}
endtemplate p2d

//* variables/declares

// arm variables
declare("aid","o[1]") // finitializehandler for arm
declare("nqaupd","o[1]") // NQS for arm updates (not all involving position changes)
declare("nqa","o[1]") // NQS for recording arm joint angles/positions in time
declare("nqMCMD","o[1]") // NQS for recording motor commands during sim (to draw figure)
declare("armAng","d[2]","ePos",new p2d(), "armPos",new p2d(),"tPos",new p2d()) // arm angles; elbow position; end-effector position; target position
declare("tAng","d[2]") // target angle
declare("sAng","d[2]") // starting angle
declare("sPos","d[2]") // starting angle
declare("minang","d[2]","maxang","d[2]") // min/max angles for each joint
declare("armLen","d[2]") // length of each arm segment -- fixed throughout sim

// muscle variables
declare("nqE","o[1]") // has E assignment to muscle groups
declare("nqDP","o[1]") // has DP assignment to muscle groups
declare("MTYP","o[1]") // has muscle names
declare("MLen","d[4]") // length of each muscle -- varies throughout sim
declare("minMLen","d[4]","maxMLen","d[4]") // min/max muscle lengths
declare("splitEM", 1) // whether to read only from half of the EM population
declare("DiffRDEM",0) // whether to take difference from last activity when reading out EM activity
declare("musExc","d[4]") // output muscle excitation
declare("vEMmax",30) // max num of spikes of output EM population (depends on mcmdspkwd) - used to calculate normalized muscle excitation (musExc) 
declare("minDPrate",0.00001,"maxDPrate",100)  // min and max firing rate of DP NSLOCs (tuned to different muscle lengths)
declare("DPnoise",0.0) // noise of NSLOC input to DP cells
declare("DPoverlap", 1) // whether to have overlap in the encoding of muscle lengths (~population coding)

// time variables
declare("aDT",10) // how often to update the whole arm apparatus
declare("amoveDT",10) // how often to update the arm's position
declare("mcmdspkwd",50) // motor command spike window; how wide to make the counting window (in ms)
declare("EMlag",50) // lag between EM commands (at neuromuscular junction) and arm update
declare("lastmovetime",0) // when was the last arm move (in ms)?
declare("spdDT",10) // how often to update the muscle spindles (DP cells)
declare("DPlag",50) // lag between arm movement and DP update
declare("lastspdupdate",0) // when was the last spindle update (in ms)?
declare("rlDT",50) // how often to check RL status
declare("lastrlupdate",0) // when was the last RL update (in ms)?
declare("muscleNoiseChangeDT",500) // how often to alternate noise to muscle groups
declare("randomMuscleDTmax",1000) // max random delay to alternate muscle groups (changes every time) (-1 = off = use fixed muscleNoiseChangeDT)
declare("lastMuscleNoiseChange",0) // when was the last RL update (in ms)?
declare("syDT",0) // dt for recording synaptic weights into nqsy -- only used when syDT>0
declare("lastSynScaling", 0) // last time syns scaling occurred
declare("synScalingDT", 1000) // how often to apply synaptic scaling

// learning variables
//maxplastt_INTF6=rlDT + 50 // spike time difference over which to turn on eligibility
declare("minRLerrchangeLTP",0.001) // minimum error change visible to RL algorithm for LTP (units in Cartesian coordinates)
declare("minRLerrchangeLTD",0.001) // minimum error change visible to RL algorithm for LTD (units in Cartesian coordinates)
declare("DoLearn",4,"DoReset",2,"DoRDM",2, "DoAnim",0) // learn; reset; use random targets for training; animate arm
 declare("invertAxis",0) // invert axis when plotting arm trajectory and target
declare("RLMode",3)  // reinforcement learning mode (0=none,1=reward,2=punishment,3=reward+punishment)
declare("targid",2)  // target ID for target to set up
declare("XYERR",0) // whether to use cartesian error
declare("ANGERR",1) // whether to use diff btwn targ and arm angle for error 
declare("TRAJERR",2) // whether to use diff btwn targ and arm angle for error 
declare("COMBERR",0) // whether to use diff btwn targ and arm angle for error 
declare("errTY",XYERR) // type of error - either XYERR or ANGERR
declare("centerOutTask", 1) // use targets from center-out reaching task (1->L=2*monkey; 2->L=1.5*monkey)

// declare("AdaptLearn",0) // whether to modulate learning level by distance from target

// recording/visualization variables
declare("lrec","o[2]") // recording from ES, EM
declare("vEM",new Vector())
declare("pnq","o[2]") // NQS for motor command map
declare("nqsy","o[2]") // for recording synaptic weights over time --  only saves L5 cell synapses
declare("lastsysave",0) // when was the nqsy synaptic weight recording (in ms)?
sprint(tstr,"o[%d][%d]",CTYPi+1,CTYPi+1)
declare("lssyavg",tstr) // list of average synaptic weights vs time for a particular population
declare("syCTYP",1) // 1=only ES->EM ; 2 = all connections

// babble noise variables
declare("lem","o[1]") // List of EM AM2 'noise' NetStims to allow modulation
declare("AdaptNoise",0) // whether to adapt noise
declare("LTDCount",0) // number of recent LTD periods - used for noise adaptation
declare("StuckCount",2) // number of periods where arm doesn't improve and should adapt noise
declare("EMNoiseRate",250)//sgrhzEE) // rate of noise inputs to EM cells
declare("EMNoiseRateInc",50) // rate by which to increase noise rate when arm gets stuck
declare("EMNoiseRateDec",25) // rate by which to decrease noise rate when arm gets stuck
declare("ResetEMNoiseRate",0) // reset EMNoiseRate to sgrhzEE @ start of run ?
declare("EMNoiseRateMax",3e3) // rate of noise inputs to EM cells
declare("EMNoiseRateMin",50) // rate of noise inputs to EM cells
declare("EMNoiseMuscleGroup", 0) // alternate muscle group (-1=no alternation; initial muscle: 0=shext; 1=shflex ; 2=elext; 3=elflex)  or pattern in expSeq
declare("exploreTot", 32) // max number of exploratory sequences to use from expSeq; if =4, activate only 1 muscle at a time
declare("expSeq","d[32]") // sequence of muscles activated during exploratory movements (in pairs)

/* --------------------------------------------------*/
/* DP/NSLOC FUNCTIONS */
/* ---------------------------------------------------*/

//* assignDP -  set DP cell muscle length range responsiveness. store info in nqDP 
proc assignDP () { local ii,jj,shmlenrangewid,elmlenrangewid,shminmlen,shmaxmlen,elminmlen,elmaxmlen localobj xo
  {nqsdel(nqDP) nqDP = new NQS("col","id","ty","mid","mids","MLmin","MLmax") nqDP.strdec("mids")}

  // Calculate the width of each subrange.
  shmlenrangewid = (maxMLen[0]-minMLen[0]) / (cedp.count() / NMUSCLES) // muscle length range / (total num of dp units / number of muscles) 
  elmlenrangewid = (maxMLen[2]-minMLen[2]) / (cedp.count() / NMUSCLES)

  // Initialize the subrange bounds to the first subrange.
  shminmlen = minMLen[0]
  shmaxmlen = minMLen[0]+shmlenrangewid
  elminmlen = minMLen[2]
  elmaxmlen = minMLen[2]+elmlenrangewid


  // Loop over the different subranges...
  for ii=0,cedp.count() / NMUSCLES - 1 {
    
    // For the shoulder cells...
    for jj=0,1 {
    xo=cedp.o(ii*NMUSCLES+jj)
          
    xo.start    = 1  // if set to -1 they never spike; set to >0 to start spiking   
      xo.number = maxDPrate*tstop // set to max spikes (eg. 1000 Hz * 30 sec)
      xo.noise = DPnoise // noise level - if 0 spike times are deterministic
      xo.interval = 1/minDPrate // set base firing rate
      
      // Set the muscle length range in the particular cell.
      if (DPoverlap) {
    xo.mlenmin = shminmlen-shmlenrangewid
    xo.mlenmax = shmaxmlen+shmlenrangewid

    // Add the muscle length bounds to the DP NQS table.
    nqDP.append(1,xo.id,DP,xo.zloc,MTYP.o(xo.zloc).s,shminmlen-shmlenrangewid,shmaxmlen+shmlenrangewid)
    }  else {
    xo.mlenmin = shminmlen
    xo.mlenmax = shmaxmlen

    // Add the muscle length bounds to the DP NQS table.
    nqDP.append(1,xo.id,DP,xo.zloc,MTYP.o(xo.zloc).s,shminmlen,shmaxmlen)
    }
    }
    // For the elbow cells...
    for jj=2,3 {
    xo=cedp.o(ii*NMUSCLES+jj)
        
      xo.start    = 1  // if set to -1 they never spike; set to >0 to start spiking   
      xo.number = maxDPrate*tstop // set to max spikes (eg. 1000 Hz * 30 sec)
      xo.noise = DPnoise // noise level - if 0 spike times are deterministic
      xo.interval = 1/minDPrate // set base firing rate
      
      if (DPoverlap) {
    xo.mlenmin = elminmlen-elmlenrangewid
    xo.mlenmax = elmaxmlen+elmlenrangewid

    // Add the muscle length bounds to the DP NQS table.
    nqDP.append(1,xo.id,DP,xo.zloc,MTYP.o(xo.zloc).s,elminmlen-elmlenrangewid,elmaxmlen+elmlenrangewid)
    }  else {
        // Set the muscle length range in the particular cell.
        xo.mlenmin = elminmlen
        xo.mlenmax = elmaxmlen

        // Add the muscle length bounds to the DP NQS table.
        nqDP.append(1,xo.id,DP,xo.zloc,MTYP.o(xo.zloc).s,elminmlen,elmaxmlen)
    }
    }    

    // Move to the next subrange.
    shminmlen += shmlenrangewid
    shmaxmlen += shmlenrangewid
    elminmlen += elmlenrangewid
    elmaxmlen += elmlenrangewid
  }
}

//* updateDP - update the DP cell drive 
proc updateDP () { local mlen, mlentime,nqarow,nqcnum localobj xo
  // The first argument is the time at which to read muscle lengths.  If none 
  // is provided, the default is the current time.
  if (numarg() > 0) {
    mlentime = $1 
    nqa.verbose = 0
    nqarow = nqa.select("t","<=",mlentime) - 1
    nqa.tog()
    nqa.verbose = 1
  } else mlentime = t

  nqcnum = nqaupd.fi("spupd") //Remember that we update the DP cells in nqaupd.
  nqaupd.v[nqcnum].x(nqaupd.v.size-1) = 1
  
  // Read arm angles from Python interface to MSM
  for i = 0, 1 {
  armAng[i] = pyobj.axf.getAngleReceived(i)
  }
  
  // Read muscle lengths from Python interface to MSM
  for i = 0, 3 {
  MLen[i] = pyobj.axf.getLengthReceived(i)
  }

  // Update the interval (firing rate) of the DP/NSLOC units according to muscle length tuning curve
  nqcnum = nqa.fi("ML0")
  for ltr(xo,cedp) {
    // select mlen value as a function of t and zloc
    if (mlentime == t) {  // use current muscle length 
    mlen = MLen[xo.zloc]
  } else {  // use muscle length stored in nqa (previous t)
    mlen = nqa.v[nqcnum+xo.zloc].x(nqarow)
  } 
    
  // update value of DP NSLOCs based on whether mlen falls in its range [mlenmin, mlenmax}  
  if ((mlen >= xo.mlenmin) && (mlen <= xo.mlenmax)) {
    xo.interval = 1000.0/maxDPrate  // interval = inverse of max firing rate (in ms)
  } else {
    xo.interval = 1000.0/minDPrate  // interval = inverse of min firing rate (in ms)
    }
  
  /*if (t <600) {
    xo.interval = 1000.0/maxDPrate
  } else if (t >600 && t<1200) {
    xo.interval = 1000.0/minDPrate
  } else if (t >1200) {
    xo.interval = 1000.0/maxDPrate
  }*/
  
  }
}



/* --------------------------------------------------*/
/* EM / ES FUNCTIONS */
/* ---------------------------------------------------*/

//* GetEMType(cell id) - returns type of motor unit - only used for EM cells
// lookup in MTYP List to see string representation
func GetEMType () { 
    if($1 < col.ix[EM] || $1 > col.ixe[EM]) return -1 // not an EM cell? return -1
  return $1 % 4  // otherwise, return type
}

//* assignEM() - create nq with list of EM cells and associated muscle
proc assignEM () { local i,ct,mid,a localobj vty
  {nqsdel(nqE) nqE = new NQS("col","id","ty","mid","mids") nqE.strdec("mids")}
  a=allocvecs(vty)
  vty.append(EM)

  if (splitEM) {
    nqE.clear(col[0].numc[EM]/2)
    for vtr(&ct,vty) for i=col[0].ix[ct]+col[0].numc[EM]/2,col[0].ixe[ct] {
      mid = i%4
      nqE.append(1,i,ct,mid,MTYP.o(mid).s)
  } 
  } else {
  nqE.clear(col[0].numc[EM])
    for vtr(&ct,vty) for i=col[0].ix[ct],col[0].ixe[ct] {
      mid = i%4
      nqE.append(1,i,ct,mid,MTYP.o(mid).s)
    }
  }
  dealloc(a)
}

//* readoutEM - read activity from EM - store results in vEM 
proc readoutEM () { local i,idx,ldx,sz
  vEM.resize(MTYP.count)  // resize to the number of muscle groups
  vEM.fill(0) 
  ldx=1 // column?
  for i=0,nqE.v.size-1 {  // for all EM units    
    idx=nqE.v[1].x(i)-col[0].ix[EM]  // obtain the relative idx value
    vEM.x(nqE.v[3].x(i)) += lrec[ldx].o(1).o(idx).size  // add all recorded spikes from that unit to corresponding count
    lrec[ldx].o(1).o(idx).resize(0) // reset to 0
  }
  if(DiffRDEM) {
    sz=nqa.v.size
    if(sz) {
      vEM.x(0) -= nqa.v[5].x(sz-1)
      vEM.x(1) -= nqa.v[6].x(sz-1)
      vEM.x(2) -= nqa.v[7].x(sz-1)
      vEM.x(3) -= nqa.v[8].x(sz-1)
      for i=0,3 if(vEM.x(i) < 0) vEM.x(i) = 0
    }
  }
}

//* recE - set up recording from EM cells
proc recE () {
  //lrec[0] = mkrecl(col[0],ES)
  lrec[1] = mkrecl(col[0],EM)
}

//* svsywts - save synaptic weights from ES cells
proc svsywts () { local i,a localobj xo,vwg,vtau,vinc,vmaxw,vidx,vt,vpre
  // print "svsywts: t " , t
  a=allocvecs(vidx,vwg,vtau,vinc,vmaxw,vt,vpre)
  if (syCTYP == 1) { //record only ES synapses
    for ltr(xo,col[0].ce) {
      if(xo.type==ES) {
        vrsz(xo.getdvi(vidx),vwg,vtau,vinc,vmaxw,vt,vpre)
        xo.getplast(vwg,vtau,vinc,vmaxw) // get current weight gains    
        vt.fill(t) //time
        vpre.fill(xo.id) //presynaptic id
        nqsy.v[0].append(vpre)
        nqsy.v[1].append(vidx)
        nqsy.v[2].append(vwg)
        nqsy.v[3].append(vt)
      }
    }
  } else (syCTYP == 2) { //record all synapses
    for ltr(xo,col[0].ce) {
        vrsz(xo.getdvi(vidx),vwg,vtau,vinc,vmaxw,vt,vpre)
        xo.getplast(vwg,vtau,vinc,vmaxw) // get current weight gains    
        vt.fill(t) //time
        vpre.fill(xo.id) //presynaptic id
        nqsy.v[0].append(vpre)
        nqsy.v[1].append(vidx)
        nqsy.v[2].append(vwg)
        nqsy.v[3].append(vt)
    } 
  }
  
  
  dealloc(a)
}

proc synScaling () {local i,a,wtSumOriginal,wtSum,scaleFactor localobj xo,vidx,delay_vec,prob_vec,wt1vec,wt1vecAbs,wt2vec,wt2vecAbs,distalsyns,wgain,vwgain,vplasttau,vplastinc,vplastmaxw
  // allocate vectors
  a=allocvecs(vidx,delay_vec,prob_vec,wt1vec,wt1vecAbs,wt2vec,wt2vecAbs,distalsyns,wgain,vwgain,vplasttau,vplastinc,vplastmaxw)
  
  wtSumOriginal=0
  wtSum=0
  scaleFactor=0

  // loop over all cells
  for ltr(xo,col[0].ce) {
    // set vector length based on number of postsynaptic connections
    vrsz(xo.getdvi(vidx),delay_vec,prob_vec,wt1vec,wt1vecAbs,wt2vec,wt2vecAbs,distalsyns,wgain)

    // read all weight gains and vwt1- getdvi
    xo.getdvi(vidx,delay_vec,prob_vec,wt1vec,wt2vec,distalsyns,wgain)

    // get absolute weight of wt1
    wt1vecAbs.abs(wt1vec)

    // get absolute weight of wt1
    wt2vecAbs.abs(wt2vec)

    if (wgain.sum() > 0) {
      // sum original weights
      wtSumOriginal=wtSumOriginal+wt1vecAbs.sum()+wt2vecAbs.sum()

      // sum current weights
      wtSum=wtSum+wt1vecAbs.mul(wgain).sum()+wt2vecAbs.mul(wgain).sum()
    }
  }
  
  // calculate normalizing constant (scaling factor)
  scaleFactor = 1.0 * wtSumOriginal / wtSum
  //print "Scaling Factor: ", scaleFactor
  
  if (scaleFactor < 1) {
    // set all weight gains via setplast
    for ltr(xo,col[0].ce) {
      //print "pre:", xo.gid
      // set vector length based on number of postsynaptic connections
      vrsz(xo.getdvi(vidx),vwgain,vplasttau,vplastinc,vplastmaxw)

      // read wgains
      xo.getplast(vwgain,vplasttau,vplastinc,vplastmaxw)

      // scale synaptic weight gain
      vwgain.mul(scaleFactor)

      // update scaled wgains
      xo.setplast(vwgain,vplasttau,vplastinc,vplastmaxw)
    }
  }
  dealloc(a)
}

  

/* --------------------------------------------------*/
/*  ARM/MUSCLE FUNCTIONS */
/* ---------------------------------------------------*/

//* initarm - init arm length, excitation, DP responsiveness, nrnpython axf - combine with initarmCB??
proc initarm () { 
 
  // arm dimensions - multiplied by 0.5 to match Alex's model dimensions 
  armLen[0] = 0.4634 - 0.173 // elbow - shoulder from MSM; previously, 0.5*(0.55) // shoulder to elbow (check if ms = wam length !)
  
  if (centerOutTask >= 3) {
    armLen[1] = 0.7169 - 0.4634 + 0.18 // radioulnar - elbow from MSM; previously, 0.5*(0.3+0.06+0.095) // elbow to center of hand = elbow to wrist + wrist length + hand length/2 (check if ms = wam legth)
  } else {
    armLen[1] = 0.7169 - 0.4634 // radioulnar - elbow from MSM; previously, 0.5*(0.3+0.06+0.095) // elbow to center of hand = elbow to wrist + wrist length + hand length/2 (check if ms = wam legth)
  }


  // min and max joint angles (in degrees)
  minang[0] = deg2rad(-10) //shoulder
  maxang[0] = deg2rad(110)
  minang[1] = deg2rad(0) //elbow
  maxang[1] = deg2rad(135)
  
  // if using normalizedFiberLength
  //minMLen[0] = 0.4  maxMLen[0] = 1.6  // define min and max muscle length values
  //minMLen[1] = 0.4  maxMLen[1] = 1.6
  //minMLen[2] = 0.4  maxMLen[2] = 1.6
  //minMLen[3] = 0.4  maxMLen[3] = 1.6
  
  // if using fiberLength
  minMLen[0] = 0.03  maxMLen[0] = 0.23 // define min and max muscle length values
  minMLen[1] = 0.03  maxMLen[1] = 0.23
  minMLen[2] = 0.03  maxMLen[2] = 0.23
  minMLen[3] = 0.03  maxMLen[3] = 0.23
  
  musExc[0] = 0 // shoulder extensor (deltoid)
  musExc[1] = 0 // shoulder flexor (pectoralis)
  musExc[2] = 0 // elbow extensor (triceps)
  musExc[3] = 0 // elbow flexor (biceps+brachialis)
  
  // create exploratory pattern of movements as pairs of muscle groups that are activated together (-1=no muscle)
  expSeq[0] = 0 expSeq[1] = -1
  expSeq[2] = 1 expSeq[3] = -1
  expSeq[4] = 2 expSeq[5] = -1
  expSeq[6] = 3 expSeq[7] = -1
  expSeq[8] = 0 expSeq[9] = 2
  expSeq[10] = 1 expSeq[11] = 3
  expSeq[12] = 0 expSeq[13] = 3
  expSeq[14] = 1 expSeq[15] = 2
  expSeq[16] = 3 expSeq[17] = -1
  expSeq[18] = 1 expSeq[19] = -1
  expSeq[20] = 2 expSeq[21] = -1
  expSeq[22] = 0 expSeq[23] = -1
  expSeq[24] = 1 expSeq[25] = 2
  expSeq[26] = 1 expSeq[27] = 3
  expSeq[28] = 0 expSeq[29] = 2
  expSeq[30] = 0 expSeq[31] = 3
  
  calculateCenterPos(sAng0,sAng1) // calculate center position (used for TRAJERR)

  assignDP() // assign DP responsiveness, store info in nqDP
  
  setTargByID(targid) // set target

  // Initialize the socket interface to the virtual or real (robotic) arm.
  nrnpython("import arminterface_pipe as axf")
  
  // run setup() in python- pyobj.axf.setup()
  sprint(tstr,"axf.setup(%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f)",tstop/1e3, aDT,sAng[0], sAng[1], tPos.x, tPos.y, damping, shExtGain, shFlexGain, elExtGain, elFlexGain)
  nrnpython(tstr)
} 

//* initArmCB - initialize nqaupd, nqa and arm callbacks 
proc initArmCB () { local i
  // Set up the muscle commands NQS table.
  nqsdel(nqaupd)
  nqaupd = new NQS("t","shext","shflex","elext","elflex","reinf","mvmt","spupd")

  // Set up the arm NQS table.
  {nqsdel(nqa) nqa=new NQS("t","ang0","ang1","x","y","shext","shflex","elext","elflex","ex","ey","phase")}
  nqa.resize("ML0","ML1","ML2","ML3","subphase","errxy","errang")
  
  // Clear out the current muscle command array vEM.
  if(vEM==nil) vEM = new Vector(MTYP.count) else {vEM.resize(MTYP.count) vEM.fill(0)}

  // initialize arm callback
  cvode.event(aDT,"updateArm()")
  
  // Do RL with dopamine
  if(DoLearn==4) {
    LearnON()
  } else LearnOFF()
  
  // create nqsy
  if(syDT>0) {nqsdel(nqsy) nqsy=new NQS("id1","id2","wg","t") cvode.event(syDT,"svsywts()")}

  // initialize last variables to 0
  LTDCount = lastmovetime = lastspdupdate = lastrlupdate = lastsysave = lastmuscle = lastMuscleNoiseChange = lastSynScaling = 0
  
  // reset EM noise rate 
  if(ResetEMNoiseRate && EMNoiseRate!=sgrhzEE) SetEMNoiseRate(EMNoiseRate=sgrhzEE)
}

//* mkmyTYP - set up names of muscle groups
proc mkmyTYP () {
  MTYP=new List()
  MTYP.append(new String("shext"))
  MTYP.append(new String("shflex"))
  MTYP.append(new String("elext"))
  MTYP.append(new String("elflex"))
}

//* updateArm - update all arm apparatus 
// This includes muscle commands, joints, muscle spindles, and reinforcement learning signals.
proc updateArm () { local errxy,errang,nqcnum,ii,tt localobj randomMuscleDT
  readoutEM()//Read M1 muscle group command spike counts since updateArm() call. These go into vEM.

  if (t == aDT) { // If we are at the beginning of the simulation...
    nqaupd.append(0,0,0,0,0,0,0,0) // Set up an arm update entry for t=0.
    // Set up an arm position entry for t=0.
    errxy = getArmErr(XYERR)   // gets error in position
    if (errTY==ANGERR) {
      errang = getArmErr(ANGERR) // gets error in angle
    } else if (errTY==TRAJERR) {
      errang = getArmErr(TRAJERR) // gets error in angle
    }

    nqa.append(0,armAng[0],armAng[1],armPos.x,armPos.y,vEM.x(0),vEM.x(1),vEM.x(2),vEM.x(3),ePos.x,ePos.y,0,MLen[0],MLen[1],MLen[2],MLen[3],0,errxy,errang)
    updateDP() // Set up the initial DP cell activity.
  }
  // Add vEM to the muscle commands part of the arm update table, nqaupd.
  nqaupd.append(t,vEM.x(0),vEM.x(1),vEM.x(2),vEM.x(3),0,-1,0)

  // If it's time for another arm motion...
  if ((t - lastmovetime >= amoveDT)) {
    if (t >= mcmdspkwd + EMlag) {
      // Get all of the motor command entries in a window mcmdspkwd ms wide, back
      // EMlag ms from the current time t.
      {nqaupd.verbose=0 nqaupd.select("t","[)",t-mcmdspkwd-EMlag,t-EMlag)}
      
      // Set vEM to the sum of each of the relevant motor command entries.
      vEM.x(0) = nqaupd.getcol("shext").sum()
      vEM.x(1) = nqaupd.getcol("shflex").sum()
      vEM.x(2) = nqaupd.getcol("elext").sum()
      vEM.x(3) = nqaupd.getcol("elflex").sum()
      
      musExc[0] = vEM.x(0) / vEMmax 
      musExc[1] = vEM.x(1) / vEMmax 
      musExc[2] = vEM.x(2) / vEMmax 
      musExc[3] = vEM.x(3) / vEMmax
   // if not enough time to collect spikes from sim (mcmdspikwindow) then use 0 excitation
   } else {
      musExc[0] = 0
      musExc[1] = 0
      musExc[2] = 0
      musExc[3] = 0
    }
    // Send a packet out to the virtual or real arm with neural excitations
    sprint(tstr,"axf.sendAndReceiveDataPackets(%f,%f,%f,%f,%f)",t, musExc[0], musExc[1], musExc[2], musExc[3])
    nrnpython(tstr)
    
    // Update the arm position
    ePos.x = armLen[0] * cos(armAng[0]) // end of elbow
    ePos.y = armLen[0] * sin(armAng[0])
    armPos.x = ePos.x + armLen[1] * cos(armAng[0] + armAng[1]) // wrist=arm position
    armPos.y = ePos.y + armLen[1] * sin(armAng[0] + armAng[1])
  
    // Get errors
    errxy = getArmErr(XYERR)   // gets error in position
    
    if (errTY==ANGERR) {
      errang = getArmErr(ANGERR) // gets error in angle
    } else if (errTY==TRAJERR) {
      errang = getArmErr(TRAJERR) // gets error in angle
    }
      
    // Set the database back.
    {nqaupd.tog() nqaupd.verbose=1}

    if (DoAnim) drarm() // Draw the arm configuration.
    lastmovetime = t    // Remember the time of the event.

    // Add the nqa row number to the nqaupd table.
    nqcnum = nqaupd.fi("mvmt")
    nqaupd.v[nqcnum].x(nqaupd.v.size - 1) = nqa.v.size - 1   
  }
  
  //If it's time for a muscle spindle cell update...
  if ((t - lastspdupdate >= spdDT) && (t >= DPlag)){
    updateDP(t-DPlag) // Update DP cell activity.    
    lastspdupdate = t // Remember the time of the event.
  }
  
  // If it's time for an RL update...
  if ((DoLearn==4) && (t - lastrlupdate > rlDT)) {
  
    // update EM muscle group that receives noise as input
    if (EMNoiseMuscleGroup > -1 && (t-lastMuscleNoiseChange >= muscleNoiseChangeDT)) {
      //print "muscle group=", EMNoiseMuscleGroup
      lastMuscleNoiseChange = t
      EMNoiseMuscleGroup += 2  // increase exploratory sequence index (seee expSeq)
      if (EMNoiseMuscleGroup >= exploreTot) { // wrap around
        EMNoiseMuscleGroup = 0
      }
      // Use random delays to alternate the muscle that receives noise during exploratory movements
      if (randomMuscleDTmax>0) {
        randomMuscleDT = new Random()
        randomMuscleDT.ACG(inputseed+iepoch+t) //initialize random number generator
        randomMuscleDT.uniform(0,1) 
        muscleNoiseChangeDT= randomMuscleDT.uniform(100,randomMuscleDTmax)
        EMNoiseMuscleGroup = int(randomMuscleDT.uniform(0,exploreTot/2))*2
        print "muscleNoiseChangeDT= ",muscleNoiseChangeDT
        print "muscle group=", EMNoiseMuscleGroup
      }
    }     
    if (verbosearm) print "Babble noise EM Muscle Group: ", EMNoiseMuscleGroup
      
    RLUpdate()       // Do an RL update.    
    lastrlupdate = t // Remember the time of the event.
  }
  
  // If it's time for synaptic scaling
  if ((DoLearn==4) && (t - lastSynScaling >= synScalingDT)) {
    synScaling()
    lastSynScaling = t
  }

  // Add a position entry to nqa.
  //print "shAng: ", armAng[0]
  //print "elAng: ", armAng[1]
  nqa.append(t,armAng[0],armAng[1],armPos.x,armPos.y,vEM.x(0),vEM.x(1),vEM.x(2),vEM.x(3),ePos.x,ePos.y,0,MLen[0],MLen[1],MLen[2],MLen[3],0,errxy,errang)
  
  if (syDT > 0 && t - lastsysave >= syDT) { // time for nqsy synaptic save
    svsywts()      // Do a synaptic save.
    lastsysave = t // Remember the time of the event.
  }  
  
  cvode.event(t+aDT,"updateArm()") // Post the next updateArm() event.

  // IF end, call python function to close sockets, save data and plot arm graphs (pass sim length as argument)
  if (t >= tstop-aDT/2) { 
   //sprint(tstr,"axf.closeSavePlot(%f, %f)",t/1000, aDT)
   //nrnpython(tstr)
   print "Finished sim..."
   if (syDT > 0 && lastsysave < t) { // save synaptic weights in last step
      svsywts()      // Do a synaptic save.
    }
  }
} 

//* mkaid - setup handlers for periodic arm update
proc mkaid () {
  aid = new FInitializeHandler(1,"initArmCB()") 
}


/* --------------------------------------------------*/
/* LEARNING/NOISE FUNCTIONS */
/* ---------------------------------------------------*/

//* getArmErr([errtype]) - gets current error in position
// 1st arg is optional - if left out, uses global errTY
func getArmErr () { local dsq,ety, xtmp, ytmp, rotAng, xrot,yrot localobj tAngs
  if(numarg()>0) ety=$1 else ety=errTY
  if(ety == XYERR) {
    // error as Cartesian distance between hand position and target position
    return sqrt((armPos.x - tPos.x)^2 + (armPos.y - tPos.y)^2)
  } else if(ety == ANGERR) {
    // error as squared difference between joint and target angles
    return sqrt((armAng[0] - tAng[0])^2 + (armAng[1] - tAng[1])^2)
  } else if(ety == TRAJERR) {
    // error as squared difference between current position and ideal (straight line) trajectory line
    xtmp = armPos.x-sPos[0]
    ytmp = armPos.y-sPos[1]

    //print "armPos.x:",armPos.x,"ArmPos.y",armPos.y
    //print "xtemp:",xtmp,"ytemp:",ytmp
    //set angle of rotation as a function of target
    tAngs = new Vector()
    tAngs.append(0, -PI, -PI/2, PI/2, -PI/4, -3*PI/4, PI/4, 3*PI/4)
    rotAng=tAngs.x[targid]

    // rotate trajectory to position over positive x-axis (=target 0, center-right)
    //xrot = xtmp*cos(rotAng) - ytmp*sin(rotAng)
    yrot = xtmp*sin(rotAng) + ytmp*cos(rotAng)
    //print "yrot(error):", yrot
    // retrun deviation from ideal trajectory
    return abs(yrot)
  } 
}
  
//* LearnON - turn on learning
proc LearnON () {
  plaststartT_INTF6 = 0
  plastendT_INTF6 = tstop * 2
}

//* LearnOFF - turn off learning
proc LearnOFF () {
  plaststartT_INTF6 = tstop * 2
  plastendT_INTF6 = tstop * 3
}

//* RLUpdate - reinforcement learning updateupdate - MODIFY
proc RLUpdate () { local err1xy,err2xy,err1ang,err2ang,err1,err2,nqcnum,combRL
  // If we have at least 2 nqa entries...
  if (nqa.v.size > 1) {
    // Find the distance error of the 2nd-to-last hand position (in 
    // nqa) from the target hand position.
    err1xy = nqa.getcol("errxy").x(nqa.v.size-2)

    // Find the distance error of the last hand position (in 
    // nqa) from the target hand position.
    err2xy = nqa.getcol("errxy").x(nqa.v.size-1)

    // Find the distance error of the 2nd-to-last hand position (in 
    // nqa) from the target hand position.
    err1ang = nqa.getcol("errang").x(nqa.v.size-2)

    // Find the distance error of the last hand position (in 
    // nqa) from the target hand position.
    err2ang = nqa.getcol("errang").x(nqa.v.size-1)

    if(XYERR == errTY) {
      err1 = err1xy
      err2 = err2xy
    } else if (ANGERR == errTY) {
      err1 = err1ang
      err2 = err2ang
    } else if (TRAJERR == errTY) {
      err1 = err1ang // store traj error in same variable as ang
      err2 = err2ang
    }
    
    if(verbosearm) print "err1 : ", err1, " err2 " , err2
    if ((DoLearn == 4) && (RLMode > 0)) { // RL with dopamine

      combRL=0
    if(COMBERR==4) {
      if( err1xy - err2xy >= minRLerrchangeLTP && err1ang >= err2ang ) { //dev
        combRL = 1 // do LTP
      } else if( err2xy - err1xy >= minRLerrchangeLTD && err1ang < err2ang ) { //deviation change
         combRL = 2 // do LTD
      }
    } else if (COMBERR==2) {
      if( err1xy - err2xy >= minRLerrchangeLTP && abs(err1ang) < 0.03) { // midev 
        combRL = 1 // do LTP
      } else if( err2xy - err1xy >= minRLerrchangeLTD && abs(err1ang) > 0.03) { // min deviation
        combRL = 2 // do LTD
      }
    } else if (COMBERR==3) {
      if( err1xy - err2xy >= minRLerrchangeLTP && abs(err1ang) < 0.02) { // midev2
        combRL = 1 // do LTP
      } else if( err1ang<err2ang && abs(err1ang) > 0.02) { // min deviation2
        combRL = 2 // do LTD
        }
    } else if (COMBERR==1) { // alternate XYERR and TRAJERR 
      if (errTY == XYERR) {
        errTY == TRAJERR
      } else if (errTY = TRAJERR) {
        errTY == XYERR
      }
    }
    
      nqcnum = nqaupd.fi("reinf")
      if (verbosearm) {print "err1: ", err1, " err2: " , err2, " err1-err2:  ", err1-err2, " minRLerrchangeLTP: ", minRLerrchangeLTP,  " minRLerrchangeLTD: ", minRLerrchangeLTD, " LTDCount: ", LTDCount}
      // Reward
      if (((err1 - err2 >= minRLerrchangeLTP) && ((RLMode == 1) || (RLMode == 3)) && COMBERR==0) || combRL==1) {
        // if(AdaptLearn) EPOTW_INTF6 = (err1 - err2) / err1 // else EPOTW_INTF6=1
        col[0].ce.o(0).dopelearn(1)                // LTP
        if (verbosearm) print "LTP!"
        nqaupd.v[nqcnum].x(nqaupd.v.size - 1) = 1  // put in nqaupd under "reinf"
        LTDCount = 0
        // print "LTDCount is now 0 , A"
        if(AdaptNoise) { // drop noise towards baseline
          EMNoiseRate -= EMNoiseRateDec
          if(EMNoiseRate < EMNoiseRateMin) EMNoiseRate = EMNoiseRateMin
          if(verbosearm) print "setting EMNoiseRate to ", EMNoiseRate
          //SetEMNoiseRate(EMNoiseRate)
        }
      // Punishment
      } else if (combRL==2 || ((err2 - err1 >= minRLerrchangeLTD) && ((RLMode == 2) || (RLMode == 3)) && COMBERR==0)) {
        // if(AdaptLearn) EDEPW_INTF6 = (err2 - err1) / err1 // else EDEPW_INTF6=1
        col[0].ce.o(0).dopelearn(-1)               // LTD
        if (verbosearm) print "LTD!"
        nqaupd.v[nqcnum].x(nqaupd.v.size - 1) = -1 // put in nqaupd
        LTDCount += 1
      } else if(err2>0 && err2==err1) LTDCount += 1
      
      // Adapt noise
      if(AdaptNoise && LTDCount >= StuckCount) { // if arm gets stuck, increase noise
        EMNoiseRate += EMNoiseRateInc
        if(EMNoiseRate > EMNoiseRateMax) EMNoiseRate = EMNoiseRateMax
        if(verbosearm) print "setting EMNoiseRate to ", EMNoiseRate
        //SetEMNoiseRate(EMNoiseRate)
        LTDCount = 0 // reset counter to 0
        // print "LTDCount is now 0 , B"
      }
      SetEMNoiseRate(EMNoiseRate)
      if (verbosearm) print "setting EMNoiseRate to ", EMNoiseRate
    }
  }
}

//* get a List of Lists with the noise NetStims for EM cells
obfunc LEMNoise () { local i,sy localobj nc,ncl,nsl,ls
  ncl=col.cstim.ncl // list of netconnections = lcstim.o(0).ncl
  nsl=col.cstim.nsl // list of netstims = lcstim.o(0).nsl
  ls=new List()  // list of lists with noise netstims
  for case(&sy,AM2,NM2,GA,GA2) ls.append(new List())  // for each synapse type add new list
  for ltr(nc,ncl,&i) {  // iterate through netconnections
    if(nc.syn.type == EM) {  // if EM
      if(nc.weight[AM2]>0) {  
        ls.o(0).append(nsl.o(i))  // add Netstim to list
      } else if(nc.weight[NM2]>0) {
        ls.o(1).append(nsl.o(i))
      } else if(nc.weight[GA]>0) {
        ls.o(2).append(nsl.o(i))
      } else if(nc.weight[GA2]>0) {
        ls.o(3).append(nsl.o(i))
      }
    }
  }
  return ls
}

//* SetEMNoiseRate(rate,index) - set rate of EM cell noise
// index: 0==AM2, 1==NM2, 2==GA, 3==GA2
proc SetEMNoiseRate () { local i,rate,idx localobj ns
  //splitEM=0
  rate=$1 if(numarg()>1) idx=$2 else idx=0  
  //ncl=col.cstim.ncl // list of netconnections = lcstim.o(0).ncl
  if (verbosearm) print "rate is ", rate, "idx is ", idx
  if(lem==nil) lem=LEMNoise() // get list of lists with NetStims
  for ltr(ns,lem.o(idx), &i) {  // iterate through netstims
  if ((splitEM && i >= col[0].numc[EM]/2) || splitEM==0) { // if splitEM only apply noise to half of EM 
    isplit = i
    if (splitEM) {isplit = i-col[0].numc[EM]/2}
    if (EMNoiseMuscleGroup == -1 || (nqE.v(3).x(isplit) == expSeq[EMNoiseMuscleGroup]) || (nqE.v(3).x(isplit) == expSeq[EMNoiseMuscleGroup+1])) { // compare muscle id of cell with current muscle group
      if(rate <= 0) {
        ns.number = 0
      } else {
        ns.interval = 1e3 / rate
        ns.number = tstop * rate / 1e3
        //print "on:", i , " interval:", ns.interval
      } 
    } else {  // set to 0 the rest of muscle groups
      ns.interval = 1e3 / EMNoiseRateMin
      //print "off:", i, " interval:", ns.interval
    }
  }
  }
  //splitEM=1
} 


/* --------------------------------------------------*/
/* DRAW ARM FUNCTIONS */
/* ---------------------------------------------------*/

//* drtarg - draw target location
proc drtarg () { local xsz,clr
  xsz = 0.05
  if(numarg()>0) clr=$1 else clr=1
  if (invertAxis) {
  drline(-tPos.x-xsz,-tPos.y-xsz,-tPos.x+xsz,-tPos.y+xsz,g,clr,4) // draw an x for the target
  drline(-tPos.x-xsz,-tPos.y+xsz,-tPos.x+xsz,-tPos.y-xsz,g,clr,4)
  } else {
    drline(tPos.x-xsz, tPos.y-xsz, tPos.x+xsz, tPos.y+xsz,g,clr,4) // draw an x for the target
    drline(tPos.x-xsz, tPos.y+xsz, tPos.x+xsz, tPos.y-xsz,g,clr,4)
  }
}

//* drarm([nqa,row from nqa,erase]) - draw arm location 
proc drarm () { local idx,shx, shy,ex,ey,wx,wy,ang0,ang1,ers,ln,xsz,she,shf,ele,elf localobj nqa,s
  ers=1 xsz=0.15
  shx = 0 //0.173 ; use shoulder joint as reference
  shy = 0  
  ln=armLen[0]+armLen[1]
  g.size(-ln,ln,-ln,ln)
  s=new String()
  // update using data from nqa
  if(numarg()>=2) {
    {nqa=$o1 idx=$2}
    if(idx < 0 || idx >= nqa.v.size) {printf("drarm ERRA: invalid index: %d,%d\n",idx,nqa.v.size) return}
    ex = nqa.v[9].x(idx)
    ey = nqa.v[10].x(idx)
    wx = nqa.v[3].x(idx)
    wy = nqa.v[4].x(idx)
    ang0 = nqa.v[1].x(idx)
    ang1 = nqa.v[2].x(idx)
    if(numarg()>=3) ers=$3 // if erase argument present
    if(DoLearn) {
      sprint(s.s,"t=%g: %g %g %s",nqa.v[0].x(idx),ang0,ang1)
    } else {
      sprint(s.s,"t=%g: %g %g",nqa.v[0].x(idx),ang0,ang1)
    }
    // if updated from current ePos and armPos
  } else {
    if(numarg()>=1) ers=$1 // if erase argument present
    if (invertAxis) {
    ex = -(shx+ePos.x)
    ey = -(shy+ePos.y)
    wx = -(shx+armPos.x)
    wy = -(shy+armPos.y)
    } else {
    ex = shx+ePos.x
    ey = shy+ePos.y
    wx = shx+armPos.x
    wy = shy+armPos.y 
  }
    ang0 = armAng[0]
    ang1 = armAng[1]
    if(0 && DoLearn) {
      sprint(s.s,"t=%g: %g %g %s",t,ang0,ang1) //,ATYP.o(APHASE).s)
    } else {
      sprint(s.s,"t=%g: %g %g",t,ang0,ang1)
    }
  }
  // if(verbose) print "ex=",ex,"ey=",ey,"ang0=",ang0,"wx=",wx,"wy=",wy,"ang1=",ang1
  if(ers) g.erase_all()
  drtarg() // draw an x for the target
  drline(shx,shy,ex,ey,g,2,4) // draw arm
  drline(ex,ey,wx,wy,g,3,4)

  if(DoAnim > 1) {
    if(NMUSCLES==2) { // draw the muscle lengths
      drline(3.1,0,3.1,MLen[0],g,2,3)
      drline(3.2,0,3.2,MLen[1],g,3,3)
    } else {
      drline(3.1,0,3.1,MLen[0],g,2,3)
      drline(3.2,0,3.2,MLen[2],g,3,3)
    }
  }
  // draw motor commands
  if(DoAnim > 2 && nqaupd!=nil) {
    {nqaupd.verbose=0 nqaupd.select("t","[)",t-mcmdspkwd-EMlag,t-EMlag)}
    she = nqaupd.getcol("shext").sum() / ( col.numc[EM] / 4 )
    shf = nqaupd.getcol("shflex").sum()  / ( col.numc[EM] / 4 )
    ele = nqaupd.getcol("elext").sum()   / ( col.numc[EM] / 4 )
    elf = nqaupd.getcol("elflex").sum()  / ( col.numc[EM] / 4 )
    drline(0,-3.1,she,-3.1,g,2,3)
    drline(0,-3.1,-shf,-3.1,g,2,3)
    drline(0,-3.2,ele,-3.2,g,3,3)
    drline(0,-3.2,-elf,-3.2,g,3,3)
    {nqaupd.tog("DB") nqaupd.verbose=1}
  }

  g.label(0.55,0.95,s.s)
  g.flush()
  doNotify()
}

//* drxytraj - draws the x,y position from nqa
proc drxytraj () { local gvt
  g=new Graph()
  g.size(0,tstop,0,2)
  {gvt=gvmarkflag gvmarkflag=0 g.erase if(nqa==nil) return}
  drarm()
  nqa.gr("y","x",0,2,1)
  g.exec_menu("View = plot")
  gvmarkflag=gvt
}

//* nqa2gif(nqa[,inc]) - saves arm locations in nqa as gif
// inc specifies how many frames to skip
proc nqa2gif () { local i,j,inc localobj nqa,s
  if(!FileExists("gif/wg")) system("mkdir gif/wg")
  nqa=$o1 s=new String() j=0
  if(numarg()>1) inc=$2 else inc=1
  for(i=0;i<nqa.v.size;i+=inc){
    drarm(nqa,i)
    sprint(s.s,"xcalc2gif gif/wg/%010d_nqa.gif",j)
    system(s.s)
    j+=1
  }  
}

//* animnqa(nqa[,startidx,endidx,delays]) - animates contents of nqa (arm position info)
proc animnqa () { local i,is,ie,del localobj nqa
  nqa=$o1
  if(numarg()>1) is=$2 else is=0
  if(numarg()>2) ie=$3 else ie=nqa.v.size-1
  if(numarg()>3) del=$4 else del=0.25e9
  for i=is,ie {
    drarm(nqa,i)
    sleepfor(0,del)
  }
}

//* whirlgif(outputfile,framesglob[,delframes]) - makes a moving gif using whirlgif
// delframes specifies whether to delete the gif frames after
func whirlgif () { local del localobj s
  if(!FileExists("/usr/site/pkgs/download/whirlgif304/whirlgif")) {
    print "whirlgif ERR0: can't find whirlgif binary!"
    return 0
  }
  if(numarg()>2) del=$3 else del=1
  s=new String2()
  sprint(s.s,"/usr/site/pkgs/download/whirlgif304/whirlgif -globmap -o %s %s",$s1,$s2)
  system(s.s)
  if(del) {system("rm gif/wg/*.gif") printf("whirlgif WARN: deleted %s !\n",$s2)}
  return 1
}

//* whirlnqa(nqa,file[,inc,delframes]) - nqa has arm trajectories
// delframes specifies whether to delete the gif frames after
func whirlnqa () { local i,del,inc localobj nqa
  nqa=$o1
  if(numarg()>3)inc=$3 else inc=1
  if(numarg()>2)del=$4 else del=1
  nqa2gif(nqa,inc)
  return whirlgif($s2, "gif/wg/*.gif",del)
}


/* --------------------------------------------------*/
/* TRAINING/TESTING FUNCTIONS */
/* ---------------------------------------------------*/

//* calculateCenterPos(angle0, angle1)
proc calculateCenterPos () {local dr0,dr1,ex,ey
  
  dr0 = $1 
  dr1 = $2
  //print "dr0", dr0, "dr1", dr1

  ex = armLen[0] * cos(dr0) // end of elbow
  ey = armLen[0] * sin(dr0)
  //print "ex", ex, "ey", ey

  sPos[0] = ex + armLen[1] * cos(dr0 + dr1) // target position for end of arm in x,y
  sPos[1] = ey + armLen[1] * sin(dr0 + dr1)
  //print "sPos0", sPos[0], "sPos1", sPos[1]
 }

//* setTarg(angle0,angle1) - set target angle
proc setTarg () { local dr0,dr1,ex,ey
  tAng[0] = deg2rad($1)
  tAng[1] = deg2rad($2)

  dr0 = tAng[0] // convert to radians
  dr1 = tAng[1]

  ex = armLen[0] * cos(dr0) // end of elbow
  ey = armLen[0] * sin(dr0)

  tPos.x = ex + armLen[1] * cos(dr0 + dr1) // target position for end of arm in x,y
  tPos.y = ey + armLen[1] * sin(dr0 + dr1)
}

//* setTargByID(targid) - sets target by code $1 - MODIFY ?
proc setTargByID () { local tid
  /* Center out targets in joint angles
   * [[ 0.50551768,  1.38993546],
       [ 0.78660115,  1.5914252 ],
       [ 0.86059412,  1.0805086 ],
       [ 0.41088479,  1.88770502],
       [ 0.71390705,  1.13083334],
       [ 0.87861959,  1.28225856],
       [ 0.36748288,  1.69689411],
       [ 0.6047993 ,  1.8411986 ]])
       * 
       * ([[ 0.45519364,  1.16293505],
       [ 0.9911765 ,  1.57515715],
       [ 1.29295714,  0.20699215],
       [ 0.20234833,  2.19736336],
       [ 0.95196038,  0.47728848],
       [ 1.16275143,  0.92036049],
       [ 0.12116826,  1.78758638],
       [ 0.67440973,  2.09186399]])

*/
  tid=$1
  // targets at 4cm * 2 * human(wrist)/monkey hand size
  if (centerOutTask==1) { 
    if (tid == 0) {
      setTarg(rad2deg(0.45519364),  rad2deg(1.16293505)) // right
    } else if (tid == 1) {
      setTarg(rad2deg(0.9911765),  rad2deg(1.57515715)) // left
    } else if (tid == 2) {
      setTarg(rad2deg(1.29295714),  rad2deg(0.20699215)) // top
    } else if (tid == 3) {
      setTarg( rad2deg(0.20234833),  rad2deg(2.19736336)) // bottom
    } else if (tid == 4) {
      setTarg(rad2deg(0.95196038),  rad2deg(0.47728848)) // right-top
    } else if (tid == 5) {
      setTarg(rad2deg(1.16275143),  rad2deg(0.92036049)) // left-top
    } else if (tid == 6) {
      setTarg(rad2deg(0.12116826),  rad2deg(1.78758638)) // left-bottom
    } else if (tid == 7) {
      setTarg(rad2deg(0.67440973) ,  rad2deg(2.09186399)) // right-bottom
    } 
  // targets at 4cm * 1.5 * human(wrist)/monkey hand size
  } else if (centerOutTask==2) {
    if (tid == 0) {
      setTarg(rad2deg(0.47152101),  rad2deg(1.2885104)) // right
    } else if (tid == 1) {
      setTarg(rad2deg(0.88510267),  rad2deg(1.59300214)) // left
    } else if (tid == 2) {
      setTarg(rad2deg(1.01646978),  rad2deg(0.77409105)) // top
    } else if (tid == 3) {
      setTarg( rad2deg(0.30815398),  rad2deg(2.047338)) // bottom
    } else if (tid == 4) {
      setTarg(rad2deg(0.80110947),  rad2deg(0.86659837)) // right-top
    } else if (tid == 5) {
      setTarg(rad2deg(1.01492911),  rad2deg(1.12032497)) // left-top
    } else if (tid == 6) {
      setTarg(rad2deg(0.24298669),  rad2deg(1.75186615)) // left-bottom
    } else if (tid == 7) {
      setTarg(rad2deg(0.62773373) ,  rad2deg(1.97372577)) // right-bottom
    } 
  // targets at 4cm * 1 * human(finger tip)/monkey hand size
  } else if (centerOutTask==3) {
    if (tid == 0) {
      setTarg(rad2deg(0.45971385),  rad2deg(1.49043484)) // right
    } else if (tid == 1) {
      setTarg(rad2deg(0.83499696),  rad2deg(1.48894409)) // left
    } else if (tid == 2) {
      setTarg(rad2deg(0.93937918),  rad2deg(1.04029084)) // top
    } else if (tid == 3) {
      setTarg( rad2deg(0.33628155),  rad2deg(1.92189911)) // bottom
    } else if (tid == 4) {
      setTarg(rad2deg(0.73274572),  rad2deg(1.18116696)) // right-top
    } else if (tid == 5) {
      setTarg(rad2deg(0.96793499),  rad2deg(1.18003431)) // left-top
    } else if (tid == 6) {
      setTarg(rad2deg(0.28174105),  rad2deg(1.79193016)) // left-bottom
    } else if (tid == 7) {
      setTarg(rad2deg(0.5882491) ,  rad2deg(1.79085662)) // right-bottom
    }  
  // targets at 4cm * 1.5 * human(finger tip)/monkey hand size
  } else if (centerOutTask==4) {
    if (tid == 0) {
      setTarg(rad2deg(0.40384714),  rad2deg(1.44023095)) // right
    } else if (tid == 1) {
      setTarg(rad2deg(0.95883405),  rad2deg(1.43798272)) // left
    } else if (tid == 2) {
      setTarg(rad2deg(1.15383275),  rad2deg(0.69393338)) // top
    } else if (tid == 3) {
      setTarg( rad2deg(0.19010291),  rad2deg(2.10147419)) // bottom
    } 
  } else {
    if (tid == 0) {
      setTarg(-15,135) // extreme flexion target
    } else if (tid == 1) {
      setTarg(-15,105)
    } else if (tid == 2) {
      setTarg(-15,75)  // normal target
    } else if (tid == 3) {
      setTarg(-15,35)
    } else if (tid == 4) {
      setTarg(-15,0)   // extreme extension target
    } else if (tid == 5) {
      setTarg(90,90)
    } else if (tid == 6) {
      setTarg(0,120)
    } else if (tid == 7) {
      setTarg(120,0)
    } else if (tid == 8) {
      setTarg(120,90)
    } else if (tid == 9) {
      setTarg(120,120)
    } else if (tid == 10) {
      setTarg(rad2deg(minang[0]),rad2deg(minang[1]))
    } else if (tid == 11) {
      setTarg(rad2deg(maxang[0]),rad2deg(maxang[1]))
    } else if (tid == 12) {
      setTarg(90,80)
    }  else if (tid == 13) {
      setTarg(90,70)
    } else if (tid == 14) {
      setTarg(90,100)
    } else if (tid == 15) {
      setTarg(100,80)
    } else if (tid == 16) {
      setTarg(100,70)   
    } else if (tid == 17) {
      setTarg(95,95)   
    } else if (tid == 18) {
      setTarg(95,0)   
      } else if (tid == 19) {
      setTarg(-10,95)   
    } else print "setTargByID ERRA: unknown tid == ", tid , " ! "
  }
}

//* NRTrain(noisemin,noisemax,noisedec,dur) - runs training by gradually decreasing noise to EM cells
proc NRTrain () { local i,j,dur,noisemax,noisemin,noisedec
  resetplast_INTF6=0
  noisemin=$1 noisemax=$2 noisedec=$3 dur=$4
  tstop = dur
  for(EMNoiseRate=noisemax;EMNoiseRate>=noisemin;EMNoiseRate-=noisedec) {
    SetEMNoiseRate(sgrhzEE=EMNoiseRate,0)
    for i=0,1 {
      if(i==0) {
        for j=0,1 sAng[j]=minang[j]
        print "EMNoiseRate = ", EMNoiseRate, ", start @ minang"
      } else {
        for j=0,1 sAng[j]=maxang[j]
        print "EMNoiseRate = ", EMNoiseRate, " start @ maxang"
      }
      run()
      pravgrates()
    }
  }
}

//* RandTrain(num locations, num iterations @ each location, duration for each iter of each location [,seed]) 
// performs training with random starting positions
proc RandTrain () { local nlocs,i,j,dur,niter,se localobj rdm
  nlocs=$1 niter=$2 dur=$3 if(numarg()>2)se=$3 else se=213951 resetplast_INTF6=0
  rdm=new Random() rdm.ACG(se)
  tstop = dur
  for i=0,nlocs-1 {
    print "location number " , i
    for j=0,1 sAng[j] = rdm.discunif(minang[j],maxang[j])
    for j=0,niter-1 run()
  }
}

//* IterTrain(number of iterations, number of increments, duration for each starting position[,savew,incvrse])
proc IterTrain () { local nangs,i,j,k,dur,niter,savew,itmp,a localobj vinc,str,nqw,vrse
  a=allocvecs(vinc,vrse) vinc.resize(2)
  niter=$1 nangs=$2 dur=$3 resetplast_INTF6=0
  if(numarg()>3) savew=$4 else savew=0
  if(numarg()>4) incvrse=$5 else incvrse=1
  if(incvrse){vrse.copy(col.cstim.vrse) itmp=initrands initrands=1}
  str=new String2()
  for i=0,1 vinc.x(i) = (maxang[i]-minang[i])/nangs
  for i=0,niter-1 {
    print "iteration number " , i
    for j=0,1 sAng[j] = minang[j]
    if(incvrse) col.cstim.vrse.copy(vrse) // reset rand seeds to 1st
    for j=0,nangs {
      tstop = dur
      run()
      for k=0,1 { // increment starting position
        sAng[k] += vinc.x(k)
        if(sAng[k] > maxang[k]) sAng[k]=maxang[k]
      }
      if(incvrse) col.cstim.vrse.add(1) // increment random seeds for new stream of external inputs
    }
    if(savew) if(i%savew==0) { // save weights
      {sprint(str.s,"data/%s_IterTrain_plastnq_iter_%d_.nqs",strv,i) nqw=getplastnq(col[0]) nqw.sv(str.s) nqsdel(nqw)}
    }    
  }
  if(incvrse) {col.cstim.vrse.copy(vrse) initrands=itmp} // restore rand seeds
  dealloc(a)
}

//* IterTrain2D(number of iterations, number of increments, duration for each starting position)
proc IterTrain2D () { local nangs,i,j,k,l,dur,niter,savew,a localobj vinc,str,nqw
  a=allocvecs(vinc) vinc.resize(2)
  niter=$1 nangs=$2 dur=$3 resetplast_INTF6=0
  if(numarg()>3) savew=$4 else savew=0
  str=new String2()
  for i=0,1 vinc.x(i) = (maxang[i]-minang[i])/nangs
  for i=0,niter-1 { sAng[0]=minang[0]
    l=0 // subiteration counter
    for j=0,nangs { sAng[1]=minang[1]
      for k=0,nangs {
        print "train: iter ", i, ", subiter ", l
        tstop = dur
        initarm()
        run()
        sAng[1] += vinc.x(1) // inc elbow angle
        if(sAng[1] > maxang[1]) sAng[1]=maxang[1]
        l += 1 // subiteration
      }
      sAng[0] += vinc.x(0) // inc shoulder angle
      if(sAng[0] > maxang[0]) sAng[0]=maxang[0]
    }
    if(savew) if(i%savew==0) { // save weights
      {sprint(str.s,"data/%s_IterTrain2D_plastnq_iter_%d_.nqs",strv,i) nqw=getplastnq(col[0]) nqw.sv(str.s) nqsdel(nqw)}
    }    
  }
  dealloc(a)
}

//* IterTest(number of iterations, number of increments, duration for each starting position[,control,incvrse,svspks])
// NB: WHEN control==1, the weights will be reset to baseline, SO USE WITH CARE!!!!!!!!!!!!!
// incvrse allows better control of sequence of random inputs - increments them by 1 for each iteration and then
// restores them to initial value at the end. this way, control and trained networks can get same random streams
// but also have variability in the streams (if number of iterations > 1).
obfunc IterTest () { local ii,nangs,i,j,k,dur,niter,dltmp,ctl,incvrse,itmp,svspks,a localobj vinc,nqo,vrse,str
  a=allocvecs(vinc,vrse) vinc.resize(2)
  niter=$1 nangs=$2 dur=$3
  if(numarg()>3) ctl=$4 else ctl=0
  if(numarg()>4) incvrse=$5 else incvrse=1
  if(numarg()>5) svspks=$6 else svspks=0
  if(incvrse){vrse.copy(col.cstim.vrse) itmp=initrands initrands=1}
  str=new String()
  if(ctl) {
    print "IterTest WARNING: reset weights to baseline!"
    resetplast_INTF6=1
  } else resetplast_INTF6=0
  dltmp=DoLearn DoLearn=0 // turn off learning
  for i=0,1 vinc.x(i) = (maxang[i]-minang[i])/nangs
  for i=0,niter-1 {
    for j=0,1 sAng[j] = minang[j]
    if(incvrse) col.cstim.vrse.copy(vrse) 
    for j=0,nangs {
      print "test: iter ", i, ", subiter ", j
      tstop = dur
      run()
      {nqa.resize("iter") nqa.pad() nqa.v[nqa.m-1].fill(i)}
      {nqa.resize("subiter") nqa.pad() nqa.v[nqa.m-1].fill(j)}
      {nqa.resize("sAng0") nqa.pad() nqa.v[nqa.m-1].fill(sAng[0])}
      {nqa.resize("sAng1") nqa.pad() nqa.v[nqa.m-1].fill(sAng[1])}
      if(nqo==nil) {nqo=new NQS() nqo.cp(nqa)} else {nqo.append(nqa)}

      if(svspks) { // save spikes?
        {CDX=0 mksnq() CDX=0} // make the NQS with spikes (snq)
        if(ctl) {
          sprint(str.s,"data/%s_iter_%d_subiter_%d_snq_control_A5.nqs",strv,i,j)
        } else {
          sprint(str.s,"data/%s_iter_%d_subiter_%d_snq_A5.nqs",strv,i,j)
        }
        snq.sv(str.s)
      }

      for k=0,1 { // increment starting position
        sAng[k] += vinc.x(k)
        if(sAng[k] > maxang[k]) sAng[k]=maxang[k]
      }
      if(incvrse) col.cstim.vrse.add(1) // increment random seeds for new stream of external inputs
    }
  }
  if(incvrse) {col.cstim.vrse.copy(vrse) initrands=itmp} // restore rand seeds
  {DoLearn=dltmp dealloc(a) return nqo}
}

//* AddCountCol(nqa) - add a column of spike counts to nqa. assumes arm in same position in the nqa.
proc AddCountCol () { local i,j,tcind,mint,maxt,nc,a localobj nqa,vcnt,vrt
  {a=allocvecs(vcnt,vrt) nqa=$o1}
  nc = col[0].cellsnq.size() // get the number of cells from cellsnq
  mksnq() // make the NQS with spike info
  {nqa.tog("DB") nqa.resize("vcnt","vrt") nqa.odec("vcnt") nqa.odec("vrt") nqa.pad()}
  snq.verbose=0
  vcnt.resize(nc)
  vrt.resize(nc)
  // Remember the column index for "t" in nqa.
  tcind = nqa.fi("t")  
  // Loop over all the nqa entries...
  for i = 0, nqa.v.size-1 {
    // If we are in the first entry (t=0)...
    if (i == 0) {
      // Set spike counts and rates to all zeros.
      vcnt.resize(0)
      vcnt.resize(nc)
      vrt.resize(0)
      vrt.resize(nc)
    } else {
      // The minimum time is the previous nqa entry time, and the max. is the current entry.
      mint = nqa.v[tcind].x(i-1)
      maxt = nqa.v[tcind].x(i)     
      // Loop over all cell indices...
      for j = 0, nc - 1 {
        // Count only those spikes from the right time interval and cell id.
        vcnt.x(j) = snq.select("id",j,"t","(]",mint,maxt)
        // Calculate the rate in spikes / s.
        vrt.x(j) = vcnt.x(j) * 1e3 / (maxt-mint)
      }
    }
    // Set the vcnt and vrt entries for this nqa entry.
    nqa.set("vcnt",i,vcnt)
    nqa.set("vrt",i,vrt)
  }
  snq.verbose=1
  dealloc(a)
}

//* IterTest2D(number of iterations, number of increments, duration for each starting position[,control,savecounts,savespikes])
// NB: WHEN control==1, the weights will be reset to baseline, SO USE WITH CARE!!!!!!!!!!!!!
// savecounts flag specifies whether to add a column to the output NQS that has # of spikes of each
// cell at each position (a Vector column). when savespikes is true, a list is returned with two objects. the first
// is the nqa and the second is a NQS with spikes for each subiteration.
obfunc IterTest2D () { local ii,nangs,i,j,k,l,dur,niter,dltmp,ctl,savec,savespks,a localobj vinc,nqo
  a=allocvecs(vinc) vinc.resize(2)
  niter=$1 nangs=$2 dur=$3
  if(numarg()>3) ctl=$4 else ctl=0
  if(numarg()>4) savec=$5 else savec=0
  if(numarg()>5) savespks=$6 else savespks=0
  if(ctl) {
    print "IterTest2D WARNING: reset weights to baseline!"
    resetplast_INTF6=1
  } else resetplast_INTF6=0
  dltmp=DoLearn DoLearn=0 // turn off learning
  for i=0,1 vinc.x(i) = (maxang[i]-minang[i])/nangs
  for i=0,niter-1 { sAng[0]=minang[0]
    l=0 // subiteration counter
    for j=0,nangs { sAng[1]=minang[1]
      for k=0,nangs {
        print "test: iter ", i, ", subiter ", l
        tstop = dur
        initarm()
        run()
        {nqa.resize("iter") nqa.pad() nqa.v[nqa.m-1].fill(i)}
        {nqa.resize("subiter") nqa.pad() nqa.v[nqa.m-1].fill(l)}
        {nqa.resize("sAng0") nqa.pad() nqa.v[nqa.m-1].fill(sAng[0])}
        {nqa.resize("sAng1") nqa.pad() nqa.v[nqa.m-1].fill(sAng[1])}
        if(savec) AddCountCol(nqa)
        if(nqo==nil) {nqo=new NQS() nqo.cp(nqa)} else {nqo.append(nqa)}

        sAng[1] += vinc.x(1) // inc elbow angle
        if(sAng[1] > maxang[1]) sAng[1]=maxang[1]

        l += 1 // subiteration
      }
      sAng[0] += vinc.x(0) // inc shoulder angle
      if(sAng[0] > maxang[0]) sAng[0]=maxang[0]
    }
  }
  {DoLearn=dltmp dealloc(a) return nqo}
}

//* Eval(nq[,noiserate,getall]) - evaluate performance at 256 different positions - meant for short intervals to see if generates
// correct motor commands - stores results in $o1 (nq) and returns the same nqs. first time Eval is called can
// pass in a nil object. default noise rate is 0. if getall==1, it will store output at each location, regardless
// of whether error was increasing or decreasing for that interval
obfunc Eval () { local i,j,k,err0,err1,nr,getall,x0,y0,x1,y1,a localobj vinc,nq
  a=allocvecs(vinc)
  nq=$o1
  if(nq==nil) nq=new NQS("sAng0","sAng1","err0","err1","x0","y0","x1","y1","sid","derr") else nq.clear()
  if(numarg()>1) nr=$2 else nr=0
  if(numarg()>2) getall=$3 else getall=0
  EMNoiseRate = sgrhzEE = nr
  SetEMNoiseRate(EMNoiseRate,0) // set noise to EM cells
  DoLearn=0 // turn off learning
  vinc.resize(2)
  for i=0,1 vinc.x(i) = (maxang[i]-minang[i])/15
  sAng[0]=minang[0]
  k=0
  for i=0,15 { sAng[1]=minang[1]
    for j=0,15 { rotArmTo(sAng[0],sAng[1])
      {err0=getArmErr() x0=armPos.x y0=armPos.y}
      run()
      {err1=getArmErr() x1=armPos.x y1=armPos.y}
      if(verbosearm) print "Eval:",sAng[0]," , ",sAng[1]," err0: ",err0,", err1: ",err1
      if(getall || (err1 >= err0 && err0 > 0)) nq.append(sAng[0],sAng[1],err0,err1,x0,y0,x1,y1,k,err1-err0)
      sAng[1] += vinc.x(1) // inc starting elbow angle
      if(sAng[1] > maxang[1]) sAng[1]=maxang[1]
      k += 1
    }
    sAng[0] += vinc.x(0) // inc starting shoulder angle
    if(sAng[0] > maxang[0]) sAng[0]=maxang[0]
  }
  dealloc(a)
  return nq
}

//* TrainAndEval([noiseratelearn,noiserateeval,durtrain,dureval,maxi,savew,evall,dureval2]) - does training from 256  positions
// and saving weights and evaluations after each iteration.
// noiseratelearn is noise rate for learning mode and noiserateeval is noise rate for 
// evaluation mode. noiseratelearn default is 200, noiserateeval default is 0. this procedure is meant for
// training using short durations (durlearn) for each iteration (default is 400 ms). dureval is the duration
// for eavaluation starting from each location. maxi is max # of iterations
// of train+eval. default is maxi of 0, which means keeps going until perfect learning attained (non-advisable,
// as it may never terminate). when savew>0 the intermediate weights will be saved every savew sessions. the global strv
// variable will then be used to determine the filenames to save to. trall determines whether to train on all starting positions
// or just those that have incorrect evaluations on last iteration.
obfunc TrainAndEval () { local i,j,durlearn,dureval,noiseratelearn,noiserateeval,maxi,savew,trall,dureval2\
                        localobj nqe,nqo,nqtmp,str,nqw
  if(numarg()>0) noiseratelearn=$1 else noiseratelearn=200
  if(numarg()>1) noiserateeval=$2 else noiserateeval=0
  if(numarg()>2) durlearn=$3 else durlearn=500
  if(numarg()>3) dureval=$4 else dureval=100
  if(numarg()>4) maxi=$5 else maxi=200
  if(numarg()>5) savew=$6 else savew=25
  if(numarg()>6) trall=$7 else trall=1
  if(numarg()>7) dureval2=$8 else dureval2=30e3
  {str=new String2() i=0}
  tstop=1
  nqtmp=Eval(nqtmp,noiserateeval,1)//learning is turned off in Eval and noiserate is set to noiserateeval
  {nqe=new NQS() nqe.cp(nqtmp) nqe.resize("session") nqe.pad() nqe.v[nqe.m-1].fill(-1)}
  i+=1 // 1 is first training
  print "did pre-eval"  
  while((i<=maxi || maxi<=0) && nqe.v.size>0) {
    print "session ", i , " nqe.v.size = " , nqe.v.size
    DoLearn = 4
    tstop=durlearn
    EMNoiseRate = sgrhzEE = noiseratelearn
    SetEMNoiseRate(EMNoiseRate,0)
    for j=0,nqtmp.v.size-1 if(trall || (nqtmp.v[3].x(j)>=nqtmp.v[2].x(j) && nqtmp.v[2].x(j)>0)) {
      sAng[0] = nqtmp.v[0].x(j)
      sAng[1] = nqtmp.v[1].x(j)
      run()
    }
    if(savew) if(i%savew==0 || i==1) {// save weights & do evaluations
      print "Evaluating progress..."
      tstop=dureval
      {nqsdel(nqtmp) nqtmp=Eval(nqtmp,noiserateeval,1)}
      {nqtmp.resize("session") nqtmp.pad() nqtmp.v[nqtmp.m-1].fill(i) nqe.append(nqtmp)}

      if(dureval2>0) { // do long-term evaluation and save output (trajectory and spikes)
        tstop=dureval2

        print "minang long-term eval..."
        {sAng[0]=minang[0] sAng[1]=minang[1] run()}
        {sprint(str.s,"data/%s_minang_nqa_session_%d_.nqs",strv,i) nqa.sv(str.s)}
        {CDX=0 mksnq() CDX=0 sprint(str.s,"data/%s_minang_snq_session_%d_.nqs",strv,i) snq.sv(str.s)}      
        
        print "maxang long-term eval..."
        {sAng[0]=maxang[0] sAng[1]=maxang[1] run()}
        {sprint(str.s,"data/%s_maxang_nqa_session_%d_.nqs",strv,i) nqa.sv(str.s)}
        {CDX=0 mksnq() CDX=0 sprint(str.s,"data/%s_maxang_snq_session_%d_.nqs",strv,i) snq.sv(str.s)}      
      }
      {sprint(str.s,"data/%s_plastnq_session_%d_.nqs",strv,i) nqw=getplastnq(col[0]) nqw.sv(str.s) nqsdel(nqw)}//save weights
    }
    i += 1
  }
  return nqe
}

//* MultTargTrain(nqc,iters,dur) - performs multiple target training. starting position is middle angle. 
// nqc must have 2 columns: targid and vrse. targid is target id. vrse is corresponding vector of random
// seeds that are copied over to col.cstim.vrse before the run. initrands will be set to 1 so that the
// Random objects are initialized properly.
proc MultTargTrain () { local i,j,dur,iters localobj nqc,str
  {initrands=1 nqc=$o1 iters=$2 dur=tstop=$3 resetplast_INTF6=0}
  for i=0,1 sAng[i] = minang[i] + (maxang[i]-minang[i]) / 2 // start in the middle
  for i=0,iters-1 {
    for j=0,nqc.v.size-1 { 
      setTargByID(targid=nqc.v.x(j))
      col.cstim.vrse.copy(nqc.get("vrse",j).o)
      print "train: targid is " , targid, " iter " , i
      run()
    }
  }
}

//* multtargtrainsv(strv,iters,dur)
proc multtargtrainsv () { local dur,iters,a localobj nqc,str,vorig,vtmp
  a=allocvecs(vorig,vtmp)
  vorig.copy(col.cstim.vrse) // copy the current random seeds
  str=new String2()
  sprint(str.s,"data/%s_",$s1) // base path for output files
  {nqc=new NQS("targid","vrse") nqc.odec("vrse")}
  nqc.append(10,col.cstim.vrse) // targ at minima in both angles
  vtmp.copy(col.cstim.vrse) // this assumes col.cstim.vrse wasn't messed with
  {vtmp.add(10) nqc.append(11,vtmp)} // targ at maxima in both angles
  iters=$2 dur=tstop=$3
  MultTargTrain(nqc,iters,dur) // does the training
  sprint(str.t,"%s_multtarg_plastnq_B1.nqs",str.s)
  {nq=getplastnq(col[0]) nq.sv(str.t) nqsdel(nq)}
  col.cstim.vrse.copy(vorig) // restore the origina random seeds
  {sprint(str.t,"%s_multtarg_nqs_B2.nqs",str.s) nqc.sv(str.t)}
  {dealloc(a) nqsdel(nqc)}
}

//* MultTargTest(nqc,dur[,control,svspks]) - does testing of mult target sim
// nqc has the random seeds used as a cue. returns nqs with trajectories
// for the different cues. this function assumes arm starts in midpoint of
// both of its angles
obfunc MultTargTest () { local i,dur,ctl,itmp,dltmp,svspks,a localobj nq,nqc,vec,vrse,nqo,str
  a=allocvecs(vrse,vec)
  nqc=$o1 dur=tstop=$2
  if(numarg()>2) ctl=$3 else ctl=0
  if(numarg()>3) svspks=$4 else svspks=0
  str=new String2()
  if(ctl) {
    print "MultTargTest WARNING: reset weights to baseline!"
    resetplast_INTF6=1
  } else resetplast_INTF6=0
  vrse.copy(col.cstim.vrse) // backup the seeds
  dltmp=DoLearn DoLearn=0 // turn off learning
  itmp=initrands initrands=1
  for i=0,1 sAng[i] = minang[i] + (maxang[i]-minang[i]) / 2 // start in the middle
  for i=0,nqc.v.size-1 {
    setTargByID(targid=nqc.v.x(i))
    print "targ " , targid 
    col.cstim.vrse.copy(nqc.get("vrse",i).o)
    run()
    {nqa.resize("targid") nqa.pad() nqa.v[nqa.m-1].fill(targid)}
    {nqa.resize("tAng0") nqa.pad() nqa.v[nqa.m-1].fill(tAng[0])}
    {nqa.resize("tAng1") nqa.pad() nqa.v[nqa.m-1].fill(tAng[1])}
    {nqa.resize("subiter") nqa.pad() nqa.v[nqa.m-1].fill(i)}
    if(nqo==nil) {nqo=new NQS() nqo.cp(nqa)} else {nqo.append(nqa)}

    if(svspks) { // save spikes?
      {CDX=0 mksnq() CDX=0} // make the NQS with spikes (snq)
      if(ctl) {
        sprint(str.s,"data/%s_multtarg_subiter_%d_snq_control_B5.nqs",strv,i)
      } else {
        sprint(str.s,"data/%s_multtarg_subiter_%d_snq_B5.nqs",strv,i)
      }
      snq.sv(str.s)
    }
  }
  {DoLearn=dltmp initrands=itmp col.cstim.vrse.copy(vrse) dealloc(a) return nqo}
}

//* multtargtestsv(simstr[,dur,twod,skipc,svspks]) - run mult targ testing and save output
// iff skipc==1 , skip control. 
// see MultTargTest for more details. svspks - save the snq spike NQS for each subiteration of MultTargTest? uses
// strv and iter,subiter,etc. information for the output filename.
proc multtargtestsv () { local iters,nl,dur,twod,c,skipc,svspks localobj str,nq,nqc
  DoLearn = syDT = 0
  str=new String2()
  sprint(str.s,"data/%s",$s1)
  sprint(str.t,"%s__multtarg_plastnq_B1.nqs",str.s)
  {nq=new NQS(str.t) setplastnq(nq,col[0]) nqsdel(nq)}// this loads the learned weights 
  print "loaded weights from ", str.t
  {sprint(str.t,"%s__multtarg_nqs_B2.nqs",str.s) nqc=new NQS(str.t)}
  print "loaded vrse from " , str.t
  if(numarg()>1) dur=$2 else dur=30e3
  if(numarg()>2) twod=$3 else twod=0
  if(numarg()>3) skipc=$4 else skipc=0
  if(numarg()>4) svspks=$5 else svspks=0
  for c=0,1 {
    if(skipc && c==1) continue
    if(twod) {
      //nq=IterTest2D(iters,nl,dur,c)
      //if(HoldStill) addhyperrcols(nq) // add the hypothetical moves if running with HoldStill==1
      //if(c==0) sprint(str.t,"%s_itertest2D_A3.nqs",str.s) else sprint(str.t,"%s_itertest2D_control_A4.nqs",str.s)
    } else {
      nq=MultTargTest(nqc,dur,c,svspks)
      if(c==0) sprint(str.t,"%s_MultTargTest1D_B3.nqs",str.s) else sprint(str.t,"%s_MultTargTest1D_control_B4.nqs",str.s)
    }
    nq.sv(str.t)
    print "saved ", str.t
    nqsdel(nq)
  }
  nqsdel(nqc)
}

//* addhyperrcols(nqs from IterTest or IterTest2D) - adds hypothetical moves/errors from 
// a IterTest with HoldStill==1, to see what the output moves would be
// with fixed positions and then to calculate the errors
proc addhyperrcols () { local i,s0,s1,errsh,errel,errp,sx,sy,nx,ny,err0,err1,dx,dy localobj nqo
  nqo=$o1 nqo.tog("DB")
  if(nqo.fi("dy")==-1) {
    nqo.resize("dirsh","direl","errsh","errel","errp","errxy","dx","dy")
    nqo.pad()
  }
  for i=0,nqo.v.size-1 {
    dirsh = nqo.getcol("shflex").x(i) - nqo.getcol("shext").x(i) // shoulder rot command
    direl = nqo.getcol("elflex").x(i) - nqo.getcol("elext").x(i) // elbow rot command
    s0 = nqo.getcol("sAng0").x(i) // starting angle
    s1 = nqo.getcol("sAng1").x(i)
    if(abs(s0 + dirsh - tAng[0] ) < abs(s0 - tAng[0])) {
      errsh = -abs(dirsh) // angular error was reduced
    } else errsh = abs(dirsh) // angular error was increased
    if(abs(s1 + direl - tAng[1] ) < abs(s1 - tAng[1])) {
      errel = -abs(direl) // angular error was reduced
    } else errel = abs(direl) // angular error was increased
    errp = errsh + errel
    
    sx = armLen[0] * cos(s0) + armLen[1] * cos(s0+s1) // get starting x
    sy = armLen[0] * sin(s0) + armLen[1] * sin(s0+s1) // get starting y
    err0 = sqrt( (sx-tPos.x)^2 + (sy-tPos.y)^2) // starting error
    
    nx = armLen[0] * cos((s0+dirsh)) + armLen[1] * cos((s0+dirsh) + (s1+direl)) // get next x
    ny = armLen[0] * sin((s0+dirsh)) + armLen[1] * sin((s0+dirsh)+(s1+direl)) // get next y
    err1 = sqrt( (nx-tPos.x)^2 + (ny-tPos.y)^2) // error after hypothetical move
    
    if(err1<err0) errxy=-abs(err1-err0) else errxy=abs(err1+err0)

    dx = nx - sx // delta-x 
    dy = ny - sy // delta-y
    
    nqo.getcol("dirsh").x(i)=dirsh // hypothetical rotation command for shoulder
    nqo.getcol("direl").x(i)=direl // hypothetical rotation command for elbow
    nqo.getcol("errsh").x(i)=errsh // hypothetical error for shoulder
    nqo.getcol("errel").x(i)=errel // hypothetical error for elbow
    nqo.getcol("errp").x(i)=errp // total hypothetical error - useless?
    nqo.getcol("errxy").x(i)=errxy // hypothetical cartesian error
    nqo.getcol("dx").x(i)=dx // hypothetical move in x 
    nqo.getcol("dy").x(i)=dy // hypothetical move in y
  }
}


/* --------------------------------------------------*/
/* ANALYSIS FUNCTIONS */
/* ---------------------------------------------------*/

//* getnqrf([celltype,colid]) - make an NQS with approx celltype RFs -- default is for DP cells - INTERESTING!!
obfunc getnqrf () { local cdx,i,tt,tdx,id,ang0,ang1,mid,ml,ct localobj nqrf
  nqrf=new NQS("i","id","ty","ang0","ang1","ML0","ML1","ML2","ML3","t")
  if(numarg()>0)ct=$1 else ct=DP
  if(numarg()>1)cdx=$2 else cdx=0
  if(ct==DP) nqrf.resize("mid","ml")
  if(snq[cdx]==nil) {i=CDX CDX=cdx mksnq() CDX=i}
  nqrf.clear(snq[cdx].select("type",DP))
  snq[cdx].tog("DB")
  nqa.tog("DB")
  nqa.verbose=snq[cdx].verbose=0
  tdx=0     // starting time
  tt=snq[cdx].getcol("t").x(0) // get first spike time
  for i=1,nqa.v.size-1 {    // loop through all arm updates
    if(i%10==0)printf(".")    
    while(tt < nqa.getcol("t").x(i) && tdx < snq[cdx].v.size) { // while spike time tt is within time of current arm update
      if(snq[cdx].getcol("type").x(tdx)==ct) {    // if the spike was from a cell type ct (eg. DP)
        id=snq[cdx].getcol("id").x(tdx)
        ang0=nqa.getcol("ang0").x(i-1)
        ang1=nqa.getcol("ang1").x(i-1)
        //      cedp.o(id).subtype
        if(ct==DP) {
          mid=cedp.o(id-col[cdx].ix[ct]).zloc
          ml=nqa.v[mid+12].x(i-1)
          nqrf.append(i-1,id,ct,ang0,ang1,nqa.v[12].x(i-1),nqa.v[13].x(i-1),nqa.v[14].x(i-1),nqa.v[15].x(i-1),tt,mid,ml)
        } else {
          nqrf.append(i-1,id,ct,ang0,ang1,nqa.v[12].x(i-1),nqa.v[13].x(i-1),nqa.v[14].x(i-1),nqa.v[15].x(i-1),tt)
        }
      }      
      tdx += 1
      if(tdx<snq[cdx].v.size) tt=snq[cdx].getcol("t").x(tdx) else break
    }
  }
  nqa.verbose=snq[cdx].verbose=1
  return nqrf
}

//* getcellrf(nqa from IterTest2D with cell counts @ each position, cell id[,normalize]) -
// returns an nqs with spike counts per position for a given cell
// normalize divides by max # of spikes at a location so they're all between 0 and 1
obfunc getcellrf () { local id,x,y,r,cnt,nrm,a localobj nq,vx,vy,vcnt,nqrf,ls,mc
  nq=$o1 id=$2 if(numarg()>2)nrm=$3 else nrm=1
  a=allocvecs(vx,vy,vcnt)
  vrsz(nq.getcol("x").uniq,vx,vy)
  nq.getcol("x").uniq(vx)
  nq.getcol("y").uniq(vy)
  nqrf=new NQS("x","y","cnt")
  for i=0,nq.v.size-1 {
    x = nq.getcol("x").x(i)
    y = nq.getcol("y").x(i)
    cnt=nq.get("vcnt",i).o.x(id)
    if(nqrf.select(-1,"x",x,"y",y)) {
      nqrf.v[2].x(nqrf.ind.x(0)) += cnt
    } else nqrf.append(x,y,cnt)
  }
  dealloc(a)
  if(nrm) nqrf.getcol("cnt").div(nqrf.getcol("cnt").max())
  return nqrf
}

//* drcellrf(nqs from getcellrf[,sz]) - draws the count per position using data from the nqs
// green dots are starting x,y . blue is ending. red line height is # spikes @ x,y.
// sz is the height,width of markers for the x,y positions
// the nqs ($o1) is obtained from getcellrf
proc drcellrf () { local i,x,y,n,sz localobj nqr
  nqr=$o1
  if(numarg()>1)sz=$2 else sz=0.01
  for i=0,nqr.v.size-1 {
    x = nqr.v[0].x(i)
    y = nqr.v[1].x(i)
    n = nqr.v[2].x(i)
    drline(x,y,x,y+n,g,2,3)
    drline(x-sz,y-sz,x+sz,y-sz,g,4,3)
    drline(x-sz,y+n+sz,x+sz,y+n+sz,g,3,3)
  }
}

//* drerrRL - draw error and reinforcement signal
proc drerrRL () { local gvt
  g=new Graph()
  g.size(0,tstop,-1,3.5)
  gvt=gvmarkflag
  gvmarkflag=1
  nqaupd.gr("reinf","t",0,2,12) // reinforcement signal in red
  nqa.gr("errxy","t",0,3,12) // cartesian error in blue
  nqa.gr("errang","t",0,4,12) // angular error in green
  gvmarkflag=0
  nqaupd.gr("reinf","t",0,2,4) // reinforcement signal in red
  nqa.gr("errxy","t",0,3,4) // cartesian error in blue
  nqa.gr("errang","t",0,4,12) // angular error in green
  gvmarkflag=gvt
}

//* drelbowtrajectory([linestyle,erase]) -- show the elbows trajectory vs. the target angle
proc drelbowtrajectory () { local ln,ers  
  if(numarg()>0) ln=$1 else ln=1
  if(numarg()>1) ers=$2 else ers=0
  if(ers) g.erase_all()
  drline(0,minang[1],t,minang[1],g,1,5)
  drline(0,maxang[1],t,maxang[1],g,1,5)
  nqa.gr("ang1","t",0,3,ln)
  drline(0,tAng[1],tstop,tAng[1],g,3,5)
}

//* drshouldertrajectory([linestyle,erase]) -- show the shoulder's trajectory vs. the target angle
proc drshouldertrajectory () { local ln,ers
  if(numarg()>0) ln=$1 else ln=1
  if(numarg()>1) ers=$2 else ers=0
  if(ers) g.erase_all()
  drline(0,minang[0],t,minang[0],g,9,5)
  drline(0,maxang[0],t,maxang[0],g,9,5)
  nqa.gr("ang0","t",0,2,ln)
  drline(0,tAng[0],tstop,tAng[0],g,2,5)
}

//* addnqacol2nqaupd
proc addnqacol2nqaupd () { local nqcnum,nqcnum2,nqcnum3,nqaind,ii
  nqaupd.resize($s1)                    // add column to look at
  nqcnum = nqaupd.fi($s1)
  nqcnum2 = nqaupd.fi("mvmt")
  nqcnum3 = nqa.fi($s1)
  nqaupd.v[nqcnum].resize(nqaupd.v.size)  // pad the column with zeros
  nqaind = -1
  for ii=0,nqaupd.v.size-1 {
    if ((nqaupd.v[nqcnum2].x(ii) != -1) && (nqaupd.v[nqcnum2].x(ii) != nqaind)) {
      nqaind = nqaupd.v[nqcnum2].x(ii)
    }
    nqaupd.v[nqcnum].x(ii) = nqa.v[nqcnum3].x(nqaind)
  }
}

//* getavgsyvst(nqsy,ty1,ty2[,CDX]) - get average synaptic weight vs time
obfunc getavgsyvst () { local i,ty1,ty2,cdx localobj ls,vt,vwg,nqsy
  nqsy=$o1 ty1=$2 ty2=$3 if(numarg()>3) cdx=$4 else cdx=0
  ls=new List()
  ls.append(vt=new Vector(nqsy.v.size))
  ls.append(vwg=new Vector(nqsy.v.size))
  vwg.resize(0)
  nqsy.tog("DB")
  nqsy.getcol("t").uniq(vt)
  nqsy.verbose=0
  for vtr(&i,vt) if(nqsy.select("t",i,"id1","[]",col[cdx].ix[ty1],col[cdx].ixe[ty1],"id2","[]",col[cdx].ix[ty2],col[cdx].ixe[ty2])) {
    vwg.append(nqsy.getcol("wg").mean)
  }
  nqsy.verbose=1
  return ls // return list with vector of times and average weights
}

//* mkavgsyvst(nqsy)
proc mkavgsyvst () { local i,j localobj str //,nqsy
  str=new String() 
  //nqsy=$o1
  
  if (syCTYP == 1) {
    for case(&i,ES) for case(&j,EM) if(pmat[0][0][i][j]) {
      print CTYP.o(i).s,"->",CTYP.o(j).s
      lssyavg[i][j] = getavgsyvst(nqsy,i,j)
      sprint(str.s,"%s->%s",CTYP.o(i).s,CTYP.o(j).s)
      lssyavg[i][j].o(1).label(str.s)
    } 
  } else if (syCTYP == 2) {
    for case(&i,ES,EM) for case(&j,ES,IS,ISL,EM,IM,IML) if(pmat[0][0][i][j]) {
      print CTYP.o(i).s,"->",CTYP.o(j).s
      lssyavg[i][j] = getavgsyvst(nqsy,i,j)
      sprint(str.s,"%s->%s",CTYP.o(i).s,CTYP.o(j).s)
      lssyavg[i][j].o(1).label(str.s)
    } 
  }
}

//* plotsyvst - plot average synaptic weight change vs time
proc plotavgsyvst () { local i,j,k localobj str
  g=new Graph()
  g.size(0,tstop,0,4)
  k=0
  if (syCTYP == 1) {
    for case(&i,ES) for case(&j,EM) if(pmat[0][0][i][j]) {
      sprint(tstr,"%s->%s",CTYP.o(i).s,CTYP.o(j).s)
      lssyavg[i][j].o(1).label(tstr)
      lssyavg[i][j].o(1).plot(g,lssyavg[i][j].o(0),k+1,1)
      k+=1
    } 
  } else if (syCTYP == 2) {
    for case(&i,ES,EM) for case(&j,ES,IS,ISL,EM,IM,IML,DP) if(pmat[0][0][i][j]) {
      sprint(tstr,"%s->%s",CTYP.o(i).s,CTYP.o(j).s)
      lssyavg[i][j].o(1).label(tstr)
      lssyavg[i][j].o(1).plot(g,lssyavg[i][j].o(0),k+1,1)
      k+=1
      // if run out of colors, change to dotted lines
      if (k == 9) {
        k =0
        gvmarkflag=1
      }
    }
  } 
}


//* drnqEval(nqs from Eval function[, scale]) - draws vectors from Eval, red for lines that decrease error
// gray for lines that increased it
proc drnqEval () { local eidx,i,dx,dy,x0,y0,x1,y1,xsz,scale,e0,e1,a localobj vx1,vy1,vx2,vy2,nq
  nq=$o1 if(numarg()>1)scale=$2 else scale=1
  {a=allocvecs(vx1,vx2,vy1,vy2) vrsz(0,vx1,vx2,vy1,vy2)}
  {rotArmTo(tAng[0],tAng[1]) drarm()}
  {nq.verbose=0 nq.tog("DB") eidx=nq.v.size-1}
  for i=0,eidx  {
    x0 = nq.getcol("x0").x(i)
    y0 = nq.getcol("y0").x(i)
    x1 = nq.getcol("x1").x(i) 
    y1 = nq.getcol("y1").x(i) 
    dx = (x1 - x0) * scale
    dy = (y1 - y0) * scale
    e0 = nq.getcol("err0").x(i)
    e1 = nq.getcol("err1").x(i)
    if(e0 > e1) {
      drline(x0,y0,x0+dx,y0+dy,g,2,3)
    } else drline(x0,y0,x0+dx,y0+dy,g,9,3)
    {vx1.append(x0) vy1.append(y0) vx2.append(x1) vy2.append(y1)}
  }
  nq.verbose=1
  vy1.mark(g,vx1,"O",6,4,1) // start  (green)
  vy2.mark(g,vx2,"O",6,3,1) // end    (blue)
  xsz=0.1
  drline(tPos.x-xsz,tPos.y-xsz,tPos.x+xsz,tPos.y+xsz,g,1,4) // draw an x for the target
  drline(tPos.x-xsz,tPos.y+xsz,tPos.x+xsz,tPos.y-xsz,g,1,4)
  dealloc(a)
}


/* --------------------------------------------------------------------------*/
/* ANALYSIS FUNCTIONS -- NOT ADAPTED TO MS ARM*/
/* --------------------------------------------------------------------------*/

//* drxyvecfield(nqo[,scale]) - draw xy trajectory from IterTest or IterTest2D NQS
// that had addhyperrefcol called on it - this NQS should be generated
// post learning with HoldStill set to 1. scale sets scaling for vectors (default==1)
proc drxyvecfield () { local eidx,i,x0,y0,x1,y1,xsz,scale,err0,err1,clr,a localobj vx1,vy1,vx2,vy2,nqo
  g=new Graph()
  ln=armLen[0]+armLen[1]
  g.size(-ln,ln,-ln,ln)
  nqo=$o1 if(numarg()>1)scale=$2 else scale=1
  {a=allocvecs(vx1,vx2,vy1,vy2) vrsz(0,vx1,vx2,vy1,vy2)}
  //{rotArmTo(tAng[0],tAng[1]) drarm()}
  {nqo.verbose=0 nqo.tog("DB") eidx=nqo.getcol("subiter").max()}
  for i=0,eidx if(nqo.select("subiter",i)) {
    x0 = nqo.getcol("x").mean()
    y0 = nqo.getcol("y").mean()
    err0 = sqrt((x0-tPos.x)^2+(y0-tPos.y)^2)
    x1 = x0 + nqo.getcol("dx").mean() // without scale for err calc
    y1 = y0 + nqo.getcol("dy").mean()
    err1 = sqrt((x1-tPos.x)^2+(y1-tPos.y)^2)
    if(err1<=err0) clr=2 else clr=9
    x1 = x0 + nqo.getcol("dx").mean() * scale // with scale for drawing
    y1 = y0 + nqo.getcol("dy").mean() * scale
    drline(x0,y0,x1,y1,g,clr,3)
    {vx1.append(x0) vy1.append(y0) vx2.append(x1) vy2.append(y1)}
  }
  nqo.verbose=1
  vy1.mark(g,vx1,"O",6,4,1) // start  (green)
  vy2.mark(g,vx2,"O",6,3,1) // end    (blue)
  xsz=0.01
  drline(tPos.x-xsz,tPos.y-xsz,tPos.x+xsz,tPos.y+xsz,g,1,4) // draw an x for the target
  drline(tPos.x-xsz,tPos.y+xsz,tPos.x+xsz,tPos.y-xsz,g,1,4)
  dealloc(a)
}

//* drrotvecfield(nqo[,xsz,scale]) - draw angular trajectory from IterTest or IterTest2D NQS
// that had addhyperrefcol called on it - this NQS should be generated
// post learning with HoldStill set to 1. xsz is size of target
proc drrotvecfield () { local eidx,i,s0,s1,n0,n1,xsz,scale,a localobj vx1,vy1,vx2,vy2,nqo
  {nqo=$o1 a=allocvecs(vx1,vx2,vy1,vy2) vrsz(0,vx1,vx2,vy1,vy2)}
  {nqo.verbose=0 nqo.tog("DB") eidx=nqo.getcol("subiter").max()}
  if(numarg()>1)xsz=$2 else xsz=5
  if(numarg()>2)scale=$3 else scale=1
  for i=0,eidx if(nqo.select("subiter",i)) {
    s0 = nqo.getcol("sAng0").mean()
    s1 = nqo.getcol("sAng1").mean()
    n0 = s0 + nqo.getcol("dirsh").mean() * scale
    n1 = s1 + nqo.getcol("direl").mean() * scale
    drline(s0,s1,n0,n1,g,2,3)
    {vx1.append(s0) vy1.append(s1) vx2.append(n0) vy2.append(n1)}
  }
  nqo.verbose=1
  vy1.mark(g,vx1,"O",6,4,1) // start  (green)
  vy2.mark(g,vx2,"O",6,3,1) // end    (blue)
  if(numarg()>1) xsz=$2 else xsz=5
  drline(tAng[0]-xsz,tAng[1]-xsz,tAng[0]+xsz,tAng[1]+xsz,g,1,4) // draw an x for the target
  drline(tAng[0]-xsz,tAng[1]+xsz,tAng[0]+xsz,tAng[1]-xsz,g,1,4)
  dealloc(a)
}

//* DPatang(shoulder angle,elbow angle) - get Vector of DP cell IDs activated at the particular angles -
obfunc DPatang () { local i,a0,a1,rel_shang, rel_elang localobj vid,c
  vid=new Vector()
  a0=$1 a1=$2 // a0 is shoulder, a1 is elbow angle
  if(verbosearm) print "WARN: rotated to ", a0, a1
  
  //cannot rotate arm so need way of mapping joint angles to muscle lengths
  //rotArmTo(a0,a1) // rotate arm to target angle so can check muscle length activations
  rel_shang= (a0-minang[0])/(maxang[0]-minang[0]) // relative shoulder angle (0-1)
  rel_elang= (a1-minang[1])/(maxang[1]-minang[1]) // relative shoulder angle (0-1)
 
  MLen[0] = minMLen[0] + (maxMLen[0]-minMLen[0])*(1-rel_shang)  // sh ext
  MLen[1] = minMLen[1] + (maxMLen[1]-minMLen[1])*rel_shang  // sh flex
  MLen[2] = minMLen[2] + (maxMLen[2]-minMLen[2])*(1-rel_elang)  // el ext
  MLen[3] = minMLen[3] + (maxMLen[3]-minMLen[3])*rel_elang  // el flex
  
  for ltr(c,cedp) { // go thru each DP cell
    if(MLen[c.zloc] >= c.mlenmin && MLen[c.zloc] <= c.mlenmax) { // check each DP cell's muscle activation
      vid.append(c.id) // save the cell's id
    }
  }
  return vid // return the Vector of IDs
}

//* ESatang(shoulder angle,elbow angle) - get Vector of ES cell IDs activated at the particular angles
obfunc ESatang () { local a,id1,id2,i,a0,a1 localobj vid,vidu,vDP,c
  a=allocvecs(vid) vidu=new Vector() 
  a0=$1 a1=$2
  vDP = DPatang(a0,a1) // first get DP cells activated at the angle, then see which ES cells they project to
  col.connsnq.verbose=0
  for vtr(&id1,vDP) if(col.connsnq.select("id1",id1,"id2","[]",col.ix[ES],col.ixe[ES])) {
    vid.append(col.connsnq.getcol("id2"))
  }
  col.connsnq.verbose=1
  vidu.resize(vid.size)
  vid.uniq(vidu) // get rid of duplicate ES cell IDs
  dealloc(a)
  return vidu // return the IDs
}

//* EMatang(shoulder angle,elbow angle) - get Vector of EM cell IDs activated at the particular angles
obfunc EMatang () { local a,id1,id2,i,a0,a1 localobj vid,vidu,vES,c
  a=allocvecs(vid) vidu=new Vector() 
  a0=$1 a1=$2
  vES = ESatang(a0,a1) // first get the ES cells activated at the angle, then see which EM cells they project to
  col.connsnq.verbose=0
  for vtr(&id1,vES) { // go thru presynaptic ES cells & check which EM cells they project to
    if(col.connsnq.select("id1",id1,"id2","[]",col.ix[EM],col.ixe[EM])){//use connectivity NQS for projections
      vid.append(col.connsnq.getcol("id2")) // save the IDs
    }
  }
  col.connsnq.verbose=1
  vidu.resize(vid.size)
  vid.uniq(vidu) // get rid of duplicate EM cell IDs
  dealloc(a)
  return vidu // return the IDs
}

//* MCMDatang(shoulder angle,elbow angle) - motor command at angle - gets motor command based on which EM cells activated
obfunc MSMCMDatang ()  { local a,id1,id2,i,a0,a1,fctr localobj vidm,vcmd,vem
  scaling = 0.05 // excitation to angle change factor
  vcmd=new Vector(4)
  a0=$1 a1=$2
  vidm=EMatang(a0,a1) // get which EM cells activated at the angle
  for vtr(&id1,vidm) vcmd.x(nqE.v[3].x(id1-col.ix[EM]))+=1
 
  vcmd.x(0) = scaling*(vcmd.x(1)/vEMmax - vcmd.x(0)/vEMmax)  
  vcmd.x(1) = scaling*(vcmd.x(3)/vEMmax - vcmd.x(2)/vEMmax)
  
  vcmd.resize(2)
  return vcmd
}

//* MSMCMDmapnq2([angle increment]) - get an NQS with motor command for each angle pair
// just uses static connectivity from DP -> ES -> EM
// adapted to muscskel arm: x1 and y1 calculated approximately
obfunc MSMCMDmapnq2 () { local a0,a1,inc,x0,y0,x1,y1 localobj vcmd,nq
  if(numarg()>0) inc=$1 else inc=0.2
  nq=new NQS("a0","a1","x0","y0","x1","y1","dsh","del")
  for(a0=0;a0<=maxang[0];a0+=inc) {
    for(a1=0;a1<=maxang[1];a1+=inc) {
    // calcualate initial position at ang a0,a1 
    ePos.x = armLen[0] * cos(a0) // end of elbow
    ePos.y = armLen[0] * sin(a0)
    armPos.x = ePos.x + armLen[1] * cos(a0 + a1) // wrist=arm position
    armPos.y = ePos.y + armLen[1] * sin(a0 + a1)
      x0 = armPos.x  
      y0 = armPos.y
      
      // obtain EM command at ang a0,a1
      // command already approximated as joint angle change in radians
      vcmd=MSMCMDatang(a0,a1)  
      
      // calculate end position after EM command 
      // only approx: doesn't take into account complex internal dynamis of muscskel arm
      a0cmd = a0 + vcmd.x(0)
      a1cmd =a1 + vcmd.x(1)
      ePos.x = armLen[0] * cos(a0cmd) // end of elbow
    ePos.y = armLen[0] * sin(a0cmd)
    armPos.x = ePos.x + armLen[1] * cos(a0cmd + a1cmd) // wrist=arm position
    armPos.y = ePos.y + armLen[1] * sin(a0cmd + a1cmd)
      x1 = armPos.x  
      y1 = armPos.y
      
      nq.append(a0,a1,x0,y0,x1,y1,vcmd.x(0),vcmd.x(1))
    }
  }
  return nq
}

//* MSMCMDmapnq([angle increment]) - get an NQS with motor command for each angle pair
// just uses static connectivity from DP -> ES -> EM
// this function has been adapted from the musculoskeletal arm -- it uses the nqa table to find the appropriate values
obfunc MSMCMDmapnq () { local a0,a1,inc,x0,y0,x1,y1,tmax,interval localobj vcmd,nq
  if(numarg()>0) inc=$1 else inc=0.1
  angInterval = inc/10
  tInterval = 10
  nq=new NQS("a0","a1","x0","y0","x1","y1")
  nqa.verbose = 0
  for(a0=0;a0<=maxang[0];a0+=inc) {
    for(a1=0;a1<=maxang[1];a1+=inc) {
    if (nqa.select("ang0","[]",a0-angInterval, a0+angInterval,"ang1","[]",a1-angInterval,a1+angInterval)) { // find rows in nqa that match desired angles
      if (nqa.size() > 1) {
      tmax=nqa.stat("t","max") // find row with latest time
      nqa.select("t","==",tmax)
      } else tmax = nqa.fetch("t")  
      x0=nqa.fetch("x")
      y0=nqa.fetch("y")
      nqa.select("t","==",tmax+tInterval)  
      x1=nqa.fetch("x")
      y1=nqa.fetch("y")
      nq.append(a0,a1,x0,y0,x1,y1)
    }
    }
  }
  return nq
}

//* drMCMDmapnq(nqs from MCMDmapnq) - draws the NQS
proc drMCMDmapnq () { local x0,y0,x1,y1,i localobj nq
  g=new Graph()
  ln=armLen[0]+armLen[1]
  g.size(-ln,ln,-ln,ln)
  nq=$o1
  nq.verbose=0
  nq.tog("DB")
  gvmarkflag=1
  nq.gr("y0","x0",0,4,4)
  for i=0,nq.v.size-1 {
    x0 = nq.getcol("x0").x(i)
    y0 = nq.getcol("y0").x(i)
    x1 = nq.getcol("x1").x(i)
    y1 = nq.getcol("y1").x(i)
    drline(x0,y0,x1,y1,g,2,3)
  }
  nq.gr("y1","x1",0,3,4)
  nq.verbose=1
}


/* --------------------------------------------------*/
/* RUN AND SAVE FUNCTIONS */
/* ---------------------------------------------------*/

//* run and init nqs objects
proc myrun () { local i localobj xo
  run()    // Do the normal run.
  // For all of the columns, make spike NQS tables, snq[].
  for CDX=0,numcols-1 {
    {nqsdel(snq[CDX]) snq[CDX]=SpikeNQS(vit[CDX].tvec,vit[CDX].vec,col[CDX])}
    snq[CDX].marksym="O"
    print "CDX:",CDX
    pravgrates()  // Show average firing rates for each column.
  }
  CDX=0 // Load the analysis tables for column 1.
  if(numarg()>0) initMyNQs()
  // make all of the usual analysis tables.
  // initAllMyNQs()
  // make all of the desired synaptic weights.
  if(syDT) mkavgsyvst(nqsy)
}

//* mysv - save output after myrun
proc mysv () { localobj s,nq
  s = new String()
  CDX=0
  // Save the snq for column 1.
  sprint(s.s,"data/%s_%s_snq.nqs",$s1,col[CDX].name)
  {snq[CDX].tog("DB") snq[CDX].sv(s.s)}
  // Save the global LFP for column 1.
  sprint(s.s,"data/%s_%s_LFP.nqs",$s1,col[CDX].name)
  {nqLFP[CDX].tog("DB") nqLFP[CDX].sv(s.s)}
  // Save the nqa table.
  {sprint(s.s,"data/%s_nqa.nqs",$s1) nqa.tog("DB") nqa.sv(s.s)}
  // Save an nq for the average synaptic weight changes.
  if(syDT) {
    nq = new NQS("t","EStoEMavgwtgn")
    nq.setcol("t",lssyavg[ES][EM].o(0))
    nq.setcol("EStoEMavgwtgn",lssyavg[ES][EM].o(1))
    sprint(s.s,"data/%s_EStoEMavgwtgn.nqs",$s1)
    {nq.tog("DB") nq.sv(s.s)}
    nqsdel(nq)
  }
/*
  {sprint(s.s,"/u/samn/intfstdp/data/%s_snq.nqs",$s1) snq.tog("DB") snq.sv(s.s)}
  if(nqrat!=nil) {
    {sprint(s.s,"/u/samn/intfstdp/data/%s_nqrat.nqs",$s1) nqrat.tog("DB") nqrat.sv(s.s)}
  }
  if(mynqp.size>0) {
    {sprint(s.s,"/u/samn/intfstdp/data/%s_mynqp.nqs",$s1) mynqp.tog("DB") mynqp.sv(s.s)}
  }
  {nqsdel(wgnq) wgnq=mkwgnq(col) sprint(s.s,"/u/samn/intfstdp/data/%s_wgnq.nqs",strv) wgnq.sv(s.s)} */
}

//* myrunsv(simstr) - run & save output
proc myrunsv () { 
  myrun()
  mysv($s1)
}

//* mytrainrunsv(simstr,type[,iters,numlevels/numlocations,dur,seed,savew]) - run with training and save weights
// type specifies training method: -2 == NRTrain,  -1 == RandTrain, 0==IterTrain, 1==IterTrain2D
// when savew is specified, save weights every savew iteration in IterTrain
// seed specifies random seed for use in RandTrain
// when using type==-2:NRTrain $s1==strv,$3==noisemin,$4==noisemax,$5==noisedec,$6==dur
proc mytrainrunsv () { local iters,nl,dur,ty,se,noisemin,noisemax,noisedec,savew localobj str
  str=new String2()
  ty=$2 // type of training
  if(ty!=-2) {
    if(numarg()>2) iters=$3 else iters=100 // # of iterations at each location
    if(numarg()>3) nl=$4 else nl=15  // # of levels for iterative training or # of locations for random training
    if(numarg()>4) dur=$5 else dur=10e3 // duration for each subiteration
    if(numarg()>5) se=$6 else se=213951 // random seed for RandTrain - determines starting locations
    if(numarg()>6) savew=$7 else savew=0
  } else {
    noisemin=$3
    noisemax=$4
    noisedec=$5
    dur=$6
  }
  sprint(str.s,"data/%s_",$s1) // base path for output files
  if(ty == -2) {
    NRTrain(noisemin,noisemax,noisedec,dur)
  } else if(ty == -1) {
    RandTrain(nl,iters,dur,se) // training with randomized start positions
  } else if(ty==1) {
    IterTrain2D(iters,nl,dur) // 2D iterative train
  } else IterTrain(iters,nl,dur,savew) // 1D iterative train
  {nqsdel(nq[0]) nq[0]=mkwgnq(col[0]) sprint(str.t,"%s_wgnq_A1.nqs",str.s) nq[0].sv(str.t)} // save output
  print "saved ", str.t
  {nqsdel(nq[1]) nq[1]=getplastnq(col[0]) sprint(str.t,"%s_plastnq_A2.nqs",str.s)  nq[1].sv(str.t)}
  print "saved ", str.t
}

//* mytestrunsv(simstr[,iters,numlevels,dur,twod,skipc,incvrse,svspks]) - run with training and save weights
// iff skipc==1 , skip control. incvrse controls the random seeds used in IterTest (not used in IterTest2D).
// see IterTest for more details. svspks - save the snq spike NQS for each subiteration of IterTest? uses
// strv and iter,subiter,etc. information for the output filename.
proc mytestrunsv () { local iters,nl,dur,twod,c,skipc,incvrse,svspks localobj str,nq
  DoLearn = syDT = 0
  str=new String2()
  sprint(str.s,"data/%s",$s1)
  sprint(str.t,"%s__plastnq_A2.nqs",str.s)
  nq=new NQS(str.t)
  setplastnq(nq,col[0]) // this loads the learned weights 
  nqsdel(nq)
  print "loaded weights from ", str.t
  if(numarg()>1) iters=$2 else iters=1
  if(numarg()>2) nl=$3 else nl=15
  if(numarg()>3) dur=$4 else dur=10e3
  if(numarg()>4) twod=$5 else twod=0
  if(numarg()>5) skipc=$6 else skipc=0
  if(numarg()>6) incvrse=$7 else incvrse=0
  if(numarg()>7) svspks=$8 else svspks=0
  for c=0,1 {
    if(skipc && c==1) continue
    if(twod) {
      nq=IterTest2D(iters,nl,dur,c)
      if(HoldStill) addhyperrcols(nq) // add the hypothetical moves if running with HoldStill==1
      if(c==0) sprint(str.t,"%s_itertest2D_A3.nqs",str.s) else sprint(str.t,"%s_itertest2D_control_A4.nqs",str.s)
    } else {
      nq=IterTest(iters,nl,dur,c,incvrse,svspks)
      if(c==0) sprint(str.t,"%s_itertest1D_A3.nqs",str.s) else sprint(str.t,"%s_itertest1D_control_A4.nqs",str.s)
    }
    nq.sv(str.t)
    print "saved ", str.t
    nqsdel(nq)
  }
}

// save plasticity
proc saveplast () {localobj nq
   strdef str 
   nq=getplastnq(col[0]) 
   sprint(str,"%s_plastnq.nqs",filestem)  
   nq.sv(str)
 }

proc saveweights () {localobj nq
    strdef outcon,outsyn 
    sprint(outcon,"%s_train-con.nqs",filestem) // Store the connectivity information (NQS version)
    sprint(outsyn,"%s_train-syn.nqs",filestem) // Store the synaptic changes (NQS version)
    print "Saving nqs with cell connectivity..."
    col[0].connsnq.sv(outcon)  // assume 1 column for the time being (0)
    print "Saving nqs with synaptic changes ..."
    nqsy.sv(outsyn)
  }

//* loadplastwam(simstr) - load weights to run sim and send trajectory in real.time to WAM
// eg.  call using: loadplastwam("12jun14.06_inputseed_4321_dvseed_534023_targid_11_")
proc loadplast () {localobj str,nq
  load_file("load.hoc")
  DoLearn = syDT = 0
  str=new String2()
  sprint(str.s,"data/%s",$s1)
  sprint(str.t,"%s_plastnq.nqs",str.s)
  nq=new NQS(str.t)
  setplastnq(nq,col[0]) // this loads the learned weights 
  nqsdel(nq)
  print "loaded weights from ", str.t
}

// loadnqa(batchsim, param1, param2, iseed, wseed)
// eg. loadnqa(13sep12_sim1, 150, 2, 1235, 120456)

obfunc loadnqa () { localobj str,nqa
  str = new String2()
  sprint(str.s,"data/%s",$s1)
  sprint(str.t,"%s/p1-%s_p2-%s_i-%s_w-%s_test-nqa.nqs",str.s, $s2, $s3, $s4,$s5)
  nqa = new NQS(str.t)
  return nqa
}

proc loadWeightsPrincipe() {localobj nqtemp, filename
  // load weights of trained network
  nqtemp = new NQS($s1)
  setplastnq(nqtemp,col[0])
  nqsdel(nqtemp)
  }
  

// plotTraj(tstop, targid) - plots angular and xy trajectory -- need to load nqa first
proc plotTraj () {
  // set arm lengths -- needed for plotting target
  armLen[0] = 0.4634 - 0.173 // elbow - shoulder from MSM; previously, 0.5*(0.55) // shoulder to elbow (check if ms = wam length !)
  armLen[1] = 0.7169 - 0.4634 // radioulnar - elbow from MSM; previously, 0.5*(0.3+0.06+0.095) // elbow to center of hand = elbow to wrist + wrist length + hand length/2 (check if ms = wam legth)

  // set input variables
  tstop = $1
  setTargByID($2)
  
  // draw elbow and shoulder traj vs time
  g=new Graph()
  g.size(0,tstop,-1,3.3)
  drshouldertrajectory()
  drelbowtrajectory()
  
  // draw x-y traj
  g=new Graph()
  ln=armLen[0]+armLen[1]
  g.size(-ln,ln,-ln,ln)
  drxytraj()
}

/* --------------------------------------------------*/
/* MS ARM FUNCTIONS (SALVA) */
/* ---------------------------------------------------*/

//* initBeforeRun ()  - required initialization prior to running msarm sim
proc initBeforeRun () {
  // initialization
  sAng[0] = sAng0
  sAng[1] = sAng1
  
  initarm()
// ?  mkaid()
}

func afterRun () {
  // close sockets
  sprint(tstr,"axf.closeSavePlot(%f, %f, '%s')",t/1000, aDT, filestem)
  //sprint(tstr,"axf.closeSavePlot(%f, %f)",t/1000, aDT)
  nrnpython(tstr)
  print "Closing sockets..."

  return pyobj.axf.getPacketLoss() // return 1 if there has been packet loss -- so can discard results
}


//* msrun () - musculoskeletal arm model single run
proc msrun () { local i localobj xo   
  initBeforeRun()
  run()    // Do the normal run.
  
  CDX=0 // select column1
  pravgrates()  // Show average firing rates for each column.
  initMyNQs() // Load the analysis tables 
  
  // plot raster
  gg() drit(0, tstop) grlines()
  
  // plot error and RL signal
  drerrRL()
  
  // if syDT>0, record and plot all of the desired synaptic weights.
  if(syDT) {
    mkavgsyvst(nqsy)
    plotavgsyvst()
  }
  
  // plot motor command map (not obtained properly -- includes noise)
  //drMCMDmapnq(MSMCMDmapnq())
  
}

//* singleTrain () 
proc singleTrain () { local i localobj xo 
  //****************
  // Train stage
  //****************
  
  // sim params
  tstop=mytstop=htmax= 10 *1e3 // sim time
  syDT = 20 // dt for recording synaptic weights into nqsy -- only used when syDT>0
  syCTYP = 2 // 1=only ES->EM ; 2 = all connections
   
   // set joint angles starting post
  sAng0 = 0.62 // starting shoulder and elbow angles
  sAng1 = 1.53 //1.57
  
  // set target
  targid = 9
  
  // set training params
  // learning params
  minRLerrchangeLTP = 0.001// minimum error change visible to RL algorithm for LTP (units in Cartesian coordinates)
  minRLerrchangeLTD = 0.001 // minimum error change visible to RL algorithm for LTD (units in Cartesian coordinates)
  DoLearn= 4 // learning mode
  errTY= ANGERR//XYERR // type of error - either XYERR or ANGERR
  
  // noise params
  EMNoiseRate=300 //sgrhzEE) // rate of noise inputs to EM cells
  EMNoiseMuscleGroup= 0 // alternate muscle group (-1=no alternation; initial muscle: 0=shext; 1=shflex ; 2=elext; 3=elflex) 
  muscleNoiseChangeDT = 500 // how often to alternate stim to different muscles 
  AdaptNoise=0 // whether to adapt noise
  StuckCount=2 // number of periods where arm doesn't improve and should adapt noise
  EMNoiseRateInc=50 // rate by which to increase noise rate when arm gets stuck
  EMNoiseRateDec=25 // rate by which to decrease noise rate when arm gets stuck
  EMNoiseRateMax=3e3 // rate of noise inputs to EM cells
  EMNoiseRateMin=50 // rate of noise inputs to EM cells
  
  // run sim
  initBeforeRun()
  run()    // Do the normal run.
  //iters=1 nl=1
  //IterTrain2D(iters,nl,tstop)
  
  // show results
  CDX=0 // select column1
  pravgrates()  // Show average firing rates for each column.
  initMyNQs() // Load the analysis tables 
   
  // plot raster
  gg() drit(0, tstop) grlines()
  
  // plot error and RL signal
  drerrRL()
  
  // if syDT>0, record and plot all of the desired synaptic weights.
  if(syDT) {
    mkavgsyvst(nqsy)
    plotavgsyvst()
  }
  
  type = 1

}

 //* singleTest ()  
proc singleTest () {local i localobj xo
  //****************
  // Test stage
  //****************
  
  // sim params
  tstop=mytstop=htmax= 10 *1e3 // sim time
  syDT = 0 // dt for recording synaptic weights into nqsy -- only used when syDT>0
  syCTYP = 1 // 1=only ES->EM ; 2 = all connections
   
   // set joint angles starting post
  sAng0 = 0.62 // starting shoulder and elbow angles
  sAng1 = 1.53 //1.57
  
  // set training params 
  // learning params (learning turned off)
  DoLearn= 0 // learning mode
  //errTY= XYERR // type of error - either XYERR or ANGERR
  
  // noise params (reduced noise after training)
  EMNoiseRate=0 //sgrhzEE) // rate of noise inputs to EM cells
  EMNoiseMuscleGroup= -1 // alternate muscle group (-1=no alternation; initial muscle: 0=shext; 1=shflex ; 2=elext; 3=elflex) 
  muscleNoiseChangeDT = 0 // how often to alternate stim to different muscles 
  AdaptNoise=0 // whether to adapt noise
  StuckCount=2 // number of periods where arm doesn't improve and should adapt noise
  EMNoiseRateInc=50 // rate by which to increase noise rate when arm gets stuck
  EMNoiseRateDec=25 // rate by which to decrease noise rate when arm gets stuck
  EMNoiseRateMax=3e3 // rate of noise inputs to EM cells
  EMNoiseRateMin=50 // rate of noise inputs to EM cells
  
  // run sim
  initBeforeRun()
  run()    // Do the normal run.
  // {iters=1 nl=8  nqiter2d=IterTest2D(iters,nl,tstop)} 
  //addhyperrcols(nqiter2d) -- not working with msarm
  
  // show results
  CDX=0 // select column1
  pravgrates()  // Show average firing rates for each column.
  initMyNQs() // Load the analysis tables 
  
  // plot raster
  gg() drit(0, tstop/1e3) grlines()
  
  // plot error and RL signal
  drerrRL()
  
  // draw elbow and shoulder traj vs time
  g=new Graph()
  g.size(0,tstop,-1,4)
  drshouldertrajectory()
  drelbowtrajectory()
  
  // draw x-y traj
  g=new Graph()
  ln=armLen[0]+armLen[1]
  g.size(-ln,ln,-ln,ln)
  drxytraj()
  
  // draw vector field - not working with msarm
  // drxyvecfield(nqiter2d)
  
  // draw motor command map - not working yet -- need to mapping between joint angles and muscle lengths
  // drMCMDmapnq(MSMCMDmapnq2())
  
  // if syDT>0, record and plot all of the desired synaptic weights.
  if(syDT) {
    mkavgsyvst(nqsy)
    plotavgsyvst()
  }
}


/* --------------------------------------------------*/
/* FUNCTION CALLS (ONLY ONCE)*/
/* ---------------------------------------------------*/
mkmyTYP()
assignEM()
recE()
mkaid()
if(DoAnim) gg()
LearnOFF() // start with no learning


Loading data, please wait...