// $Id: arm.hoc,v 1.252 2012/07/21 21:55:29 samn Exp $ 

DOPE_INTF6 = 1
EDOPE_INTF6 = 1
IDOPE_INTF6 = 0
ESTDP_INTF6 = ISTDP_INTF6 = 0
declare("verbosearm",0)

{CTYP.o(DP).s="P" CTYP.o(IML).s="ILM" CTYP.o(ISL).s="ILS"} // labels for cells

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

declare("aid","o[1]") // finitializehandler for arm
declare("nqaupd","o[1]") // NQS for arm updates (not all involving position changes)

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("armSeg",2) // number of segments
declare("armLen","d[2]") // length of each arm segment -- fixed throughout sim
declare("MLen","d[4]") // length of each muscle -- varies throughout sim
declare("MFctr",1)
declare("minMLen","d[4]","maxMLen","d[4]") // min/max muscle lengths
declare("rotfctr","d[2]") // multiplies diff in muscle group activation into rotation angle
declare("rotNorm",0) // normalize rotation by distance from origin?
declare("minang","d[2]","maxang","d[2]") // min/max angles for each joint
declare("aDT",5) // how often to update the whole arm apparatus
declare("amoveDT",50) // 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",50) // how often to update the muscle spindles (DP cells)
declare("DPlag",25) // 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)?
maxplastt_INTF6=rlDT + 50
declare("minRLerrchangeLTP",0.01) // minimum error change visible to RL algorithm for LTP (units in Cartesian coordinates)
declare("minRLerrchangeLTD",0.01) // minimum error change visible to RL algorithm for LTD (units in Cartesian coordinates)
declare("lrec","o[2]") // recording from ES, EM
declare("vEM",new Vector())
declare("armAng","d[2]","armPos",new p2d(),"ePos",new p2d(),"tPos",new p2d())
declare("tAng","d[2]") // target angle
declare("sAng","d[2]") // starting angle
declare("APHASE",0,"AMOVE",0,"AHOLD",1,"ARESET",2,"HoldDur",400,"MoveDur",300,"ResetDur",100)
declare("DoLearn",0,"DoReset",2,"DoAnim",1,"DoRDM",2)
declare("RLMode",3)  // reinforcement learning mode (0=none,1=reward,2=punishment,3=reward+punishment)
declare("ATYP",new List())
declare("SUBPHASE",0)
declare("ardm",new Random(),"aseed",913015) // for random set of arm pos during training
declare("pnq","o[2]")
declare("DiffRDEM",0) // whether to take difference from last activity when reading out EM activity
declare("nqsy","o[2]") // for recording synaptic weights over time --  only saves L5 cell synapses
declare("syDT",0) // dt for recording synaptic weights into nqsy -- only used when syDT>0
declare("lastsysave",0) // when was the nqsy synaptic weight recording (in ms)?
declare("shlock",0) // whether to lock shoulder
declare("ellock",0) // whether to lock elbow
declare("targid",11)  // target ID for target to set up
sprint(tstr,"o[%d][%d]",CTYPi+1,CTYPi+1)
declare("lssyavg",tstr) // list of average synaptic weights vs time for a particular population
declare("XYERR",0) // whether to use diff btwn targ and arm angle for error -- defaults to 0==cartesian error
declare("ANGERR",1) // whether to use diff btwn targ and arm angle for error -- defaults to 0==cartesian error
declare("errTY",XYERR) // 
declare("COMBERR",0) // whether to check angular error in combination with cartesian to ensure good path to targ
declare("HoldStill",0) // whether to hold the arm still - useful for debugging EM output at a given position
// declare("AdaptLearn",0) // whether to modulate learning level by distance from target
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",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",1) // reset EMNoiseRate to sgrhzEE @ start of run ?
declare("EMNoiseRateMax",3e3) // rate of noise inputs to EM cells
declare("MLenIndepArmLen",1) // whether MLen values are independent of arm segment length
declare("maxang0",135,"maxang1",135) // maxangs

//* updateMLen - update muscle lengths based on current angles
proc updateMLen () {
  if(MLenIndepArmLen) { // is muscle length independent of arm segment length?
   if(NMUSCLES==2) {
      // Shoulder extensor muscle length is upper-arm length under maximum flex, zero under minimum flex.
      MLen[0] = (armAng[0] - minang[0]) / (maxang[0] - minang[0])
      // Elbow extensor muscle length is forearm length under maximum flex, zero under minimum flex.
      MLen[1] = (armAng[1] - minang[1]) / (maxang[1] - minang[1])
      MLen[2] = MLen[3] = 0 // not using these
    } else if(NMUSCLES==4) { // redundant muscles when using 4 muscles for 2 joints      
      // Shoulder extensor muscle length is upper-arm length under maximum flex, zero under minimum flex.
      MLen[0] = (armAng[0] - minang[0]) / (maxang[0] - minang[0])      
      // Shoulder flexor muscle length is upper-arm length under minimum flex, zero under maximum flex.
      MLen[1] = 1 - MLen[0]      
      // Elbow extensor muscle length is forearm length under maximum flex, zero under minimum flex.
      MLen[2] = (armAng[1] - minang[1]) / (maxang[1] - minang[1])      
      // Elbow flexor muscle length is forearm length under minimum flex, zero under maximum flex.
      MLen[3] = 1 - MLen[2]
    }
  } else { // muscle lengths (MLen) not independent of arm segment length (armLen)
    if(NMUSCLES==2) {
      // Shoulder extensor muscle length is upper-arm length under maximum flex, zero under minimum flex.
      MLen[0] = armLen[0] * (armAng[0] - minang[0]) / (maxang[0] - minang[0])
      // Elbow extensor muscle length is forearm length under maximum flex, zero under minimum flex.
      MLen[1] = armLen[1] * (armAng[1] - minang[1]) / (maxang[1] - minang[1])
      MLen[2] = MLen[3] = 0 // not using these
    } else if(NMUSCLES==4) { // redundant muscles when using 4 muscles for 2 joints
      
      // Shoulder extensor muscle length is upper-arm length under maximum flex, zero under minimum flex.
      MLen[0] = armLen[0] * (armAng[0] - minang[0]) / (maxang[0] - minang[0])
      
      // Shoulder flexor muscle length is upper-arm length under minimum flex, zero under maximum flex.
      MLen[1] = armLen[0] - MLen[0]
      
      // Elbow extensor muscle length is forearm length under maximum flex, zero under minimum flex.
      MLen[2] = armLen[1] * (armAng[1] - minang[1]) / (maxang[1] - minang[1])
      
      // Elbow flexor muscle length is forearm length under minimum flex, zero under maximum flex.
      MLen[3] = armLen[1] - MLen[2]
    }
  }
}

//* resetarm - put arm in starting position
proc resetarm () { local dr0,dr1
  armAng[0] = sAng[0]
  armAng[1] = sAng[1]

  dr0 = deg2rad(armAng[0]) // convert to radians
  dr1 = deg2rad(armAng[1])

  ePos.x = armLen[0] * cos(dr0) // end of elbow
  ePos.y = armLen[0] * sin(dr0)

  armPos.x = ePos.x + armLen[1] * cos(dr0 + dr1) // wrist=arm position
  armPos.y = ePos.y + armLen[1] * sin(dr0 + dr1)

  vEM.fill(0)
  updateMLen()
}

//* initarmang - init arm angles - checks lock variables
proc initarmang () {
  if(shlock) {  // shoulder-lock configuration
    minang[0] = -15.1
    maxang[0] = -14.9
    minang[1] = 0
    maxang[1] = 135
    sAng[0] = -15
    sAng[1] = 90 
  } else if(ellock) { // elbow-lock straight-arm configuration
    minang[0] = -45
    maxang[0] = 135
    minang[1] = -0.1
    maxang[1] = 0.1
    sAng[0] = 45
    sAng[1] = 0 
  } else {
    minang[0] = -45
    maxang[0] = maxang0
    minang[1] = 0
    maxang[1] = maxang1
    sAng[0] = minang[0]
    sAng[1] = minang[1]
  }
}

//* assignDP - set DP cell muscle length range responsiveness. store info in nqDP
proc assignDP () { local ii,jj,shmlenrangewid,elmlenrangewid,shminmlen,shmaxmlen,elminmlen,elmaxmlen,shinc,elinc localobj xo
  {nqsdel(nqDP) nqDP = new NQS("col","id","ty","mid","mids","MLmin","MLmax") nqDP.strdec("mids")}
  if(MLenIndepArmLen) {
    shinc = 1.0 / (cedp.count() / NMUSCLES) //increments for centers of muscle-length tuning
    elinc = 1.0 / (cedp.count() / NMUSCLES)

    shmlenrangewid = 1.0 / (cedp.count() / NMUSCLES) //range of muscle lengths over which DP cells responsive
    elmlenrangewid = 1.0 / (cedp.count() / NMUSCLES)
  } else {
    shinc = armLen[0] / (cedp.count() / NMUSCLES) //increments for centers of muscle-length tuning
    elinc = armLen[1] / (cedp.count() / NMUSCLES)

    shmlenrangewid = armLen[0] / (cedp.count() / NMUSCLES) //range of muscle lengths over which DP cells responsive
    elmlenrangewid = armLen[1] / (cedp.count() / NMUSCLES)
  }
  shminmlen = 0
  shmaxmlen = shmlenrangewid
  elminmlen = 0
  elmaxmlen = elmlenrangewid
  for ii=0,cedp.count() / NMUSCLES - 1 {
    if(NMUSCLES==2) {
      {xo=cedp.o(ii*NMUSCLES+0) xo.setmlenrange(shminmlen,shmaxmlen)}
      nqDP.append(1,xo.id,DP,xo.zloc,MTYP.o(xo.zloc).s,shminmlen,shmaxmlen)

      {xo=cedp.o(ii*NMUSCLES+1) xo.setmlenrange(elminmlen,elmaxmlen)}
      nqDP.append(1,xo.id,DP,xo.zloc,MTYP.o(xo.zloc).s,elminmlen,elmaxmlen)
    } else if(NMUSCLES==4) {
      for jj=0,1 {
        {xo=cedp.o(ii*NMUSCLES+jj) xo.setmlenrange(shminmlen,shmaxmlen)}
        nqDP.append(1,xo.id,DP,xo.zloc,MTYP.o(xo.zloc).s,shminmlen,shmaxmlen)
      }
      for jj=2,3 {
        {xo=cedp.o(ii*NMUSCLES+jj) xo.setmlenrange(elminmlen,elmaxmlen)}
        nqDP.append(1,xo.id,DP,xo.zloc,MTYP.o(xo.zloc).s,elminmlen,elmaxmlen)
      }    
    } else print "assignDP ERRA: invalid NMUSCLES = ", NMUSCLES
    shminmlen += shinc
    shmaxmlen += shinc
    elminmlen += elinc
    elmaxmlen += elinc
  }
}

//* initarm - init arm location/params
proc initarm () { local ii
  armLen[0] = 1 // shoulder to forearm
  armLen[1] = 2 // forearm to wrist
  for ii=0,1 rotfctr[ii] = 1 // max of 1 degrees
  initarmang() // initialize arm angles, checks lock variables
  assignDP() // assign DP responsiveness, store info in nqDP
  resetarm()
}

//* mkmyTYP - set up names of muscle groups and phases
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"))
  ATYP=new List()
  ATYP.append(new String("Move"))
  ATYP.append(new String("Hold"))
  ATYP.append(new String("Reset"))
}

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

//* 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()
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)
  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) vEM.fill(0) ldx=1
  for i=0,nqE.v.size-1 {    
    idx=nqE.v[1].x(i)-col[0].ix[EM]
    vEM.x(nqE.v[3].x(i)) += lrec[ldx].o(1).o(idx).size
    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
    }
  }
}

//* getarmX(shoulder angle in radians, elbow angle in radians) -
// gets the end effector x coordinate
func getarmX () { local ex
  ex = armLen[0] * cos($1) // end of elbow
  return ex + armLen[1] * cos($1+$2) // wrist=arm position
}

//* getarmY(shoulder angle in radians, elbow angle in radians) -
// gets the end effector y coordinate
func getarmY () { local ey
  ey = armLen[0] * sin($1)
  return ey + armLen[1] * sin($1+$2)
}

//* clampval(pointer to value, min value, max value) - clamps value to min,max
proc clampval () { local val
  $&1 = MAXxy($&1,$2)
  $&1 = MINxy($&1,$3)
}

//* rotArm(angle0,angle1) - rotate arm by angle0,angle1
proc rotArm () { local dsh,delb,dr0,dr1,i
  dsh=$1 delb=$2
  armAng[0] += dsh // inc angles
  armAng[1] += delb
  for i=0,1 clampval(&armAng[i],minang[i],maxang[i])

  dr0 = deg2rad(armAng[0]) // convert to radians
  dr1 = deg2rad(armAng[1])

  ePos.x = armLen[0] * cos(dr0) // end of elbow
  ePos.y = armLen[0] * sin(dr0)

  armPos.x = ePos.x + armLen[1] * cos(dr0 + dr1) // wrist=arm position
  armPos.y = ePos.y + armLen[1] * sin(dr0 + dr1)

  updateMLen() // update the muscle lengths

//  if(verbosearm<0) print "x,y=",armPos.x,armPos.y
}

//* rotArmTo(angle0,angle1) - rotate the arm to angle0,angle1
proc rotArmTo () { local dsh,delb
  dsh = $1 - armAng[0]
  delb = $2 - armAng[1]
  rotArm(dsh,delb)
}

//* getArmErr([errtype]) - gets current error in position
// 1st arg is optional - if left out, uses global errTY
func getArmErr () { local dsq,ety
  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)
  } 
}

//* holdArmPos - hold arm position at target location
proc holdArmPos () { local dsh,delb,ph,errxy,errang
  if(numarg()>=2) {
    dsh = $1 - armAng[0]
    delb = $2 - armAng[1]
    if(numarg()>2) ph=$3 else ph=AMOVE
  } else {
    ph = AHOLD
    dsh = tAng[0] - armAng[0]
    delb = tAng[1] - armAng[1]
    if(APHASE == AMOVE) return // last event on queue for hold gets skipped
  }

  rotArm(dsh,delb)

  errxy = getArmErr(XYERR)   // gets error in position
  errang = getArmErr(ANGERR) // gets error in angle

  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,ph,MLen[0],MLen[1],MLen[2],MLen[3],SUBPHASE,errxy,errang)
  if(DoAnim) drarm()
}

//* setArmPos - set/save arm position using new joint angles, based on M1 activity in vEM
proc setArmPos () { local dsh,delb,ph,errxy,errang,fctr
  if (numarg()>0) ph=$1 else ph=AMOVE

  if(rotNorm) { // normalize rotation dx,dy by distance from origin
    if(armPos.x || armPos.y) fctr = 1 / sqrt(armPos.x^2 + armPos.y^2) else fctr = 1
    dsh = rotfctr[0] * fctr * (vEM.x(1) - vEM.x(0)) // should use a sigmoid (?)
    delb = rotfctr[1] * fctr * (vEM.x(3) - vEM.x(2))
  } else {
    dsh = rotfctr[0] * (vEM.x(1) - vEM.x(0)) // should use a sigmoid (?)
    delb = rotfctr[1] * (vEM.x(3) - vEM.x(2))
  }

  if(!HoldStill) rotArm(dsh,delb)

  errxy = getArmErr(XYERR)   // gets error in position
  errang = getArmErr(ANGERR) // gets error in angle

  // Add a position entry to nqa.
  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,APHASE,MLen[0],MLen[1],MLen[2],MLen[3],SUBPHASE,errxy,errang)
}

//* updateDP - update the DP cell drive
proc updateDP () { local 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

  // Loop over all DP cells and update them.
  nqcnum = nqa.fi("ML0")
  for ltr(xo,cedp) {
    if (mlentime == t) {
      xo.updatedrive()//Call cell update, having it use the current MLen[zloc].
    } else {
      xo.updatedrive(nqa.v[nqcnum+xo.zloc].x(nqarow))//Call cell update passing in appropriate muscle length.
    }
  }
}

//* resetArmPos - put arm back in starting position, updating nqa and vEM settings
proc resetArmPos () {
  if(DoReset == 1) {
    holdArmPos(sAng[0],sAng[1],ARESET)
  } else if(DoReset==2) {
    holdArmPos(ardm.discunif(minang[0],maxang[0]),ardm.discunif(minang[1],maxang[1]),ARESET)
  }
}

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

//* svsywts - save synaptic weights from L5 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)
  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)
    }
  }
  dealloc(a)
}

//* RLUpdate - reinforcement learning updateupdate
proc RLUpdate () { local err1xy,err2xy,err1ang,err2ang,err1,err2,nqcnum,cflip
  // 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 {
      err1 = err1ang
      err2 = err2ang
    }


    if(verbosearm) print "err1 : ", err1, " err2 " , err2
    if ((DoLearn == 4) && (RLMode > 0)) { // RL with dopamine

      cflip=0
      if(COMBERR) {
        if( err1xy - err2xy >= minRLerrchangeLTP && err1ang < err2ang ) {
          cflip = 1 // do LTD
        } //else if( err2xy - err1xy >= minRLerrchangeLTD && err1ang > err2ang ) {
          //cflip = 1 // do LTP ??
          // }
      }

      nqcnum = nqaupd.fi("reinf")
      // print "err1 ", err1, " err2 " , err2, "LTDCount", LTDCount
      if ((err1 - err2 >= minRLerrchangeLTP) && ((RLMode == 1) || (RLMode == 3)) && !cflip) {
        // if(AdaptLearn) EPOTW_INTF6 = (err1 - err2) / err1 // else EPOTW_INTF6=1
        col[0].ce.o(0).dopelearn(1)                // LTP
        nqaupd.v[nqcnum].x(nqaupd.v.size - 1) = 1  // put in nqaupd
        LTDCount = 0
        // print "LTDCount is now 0 , A"
        if(AdaptNoise) { // drop noise towards baseline
          EMNoiseRate -= EMNoiseRateDec
          if(EMNoiseRate < 0) EMNoiseRate = 0
          if(verbosearm) print "setting EMNoiseRate to ", EMNoiseRate
          SetEMNoiseRate(EMNoiseRate)
        }
      } else if (cflip || ((err2 - err1 >= minRLerrchangeLTD) && ((RLMode == 2) || (RLMode == 3)))) {
        // if(AdaptLearn) EDEPW_INTF6 = (err2 - err1) / err1 // else EDEPW_INTF6=1
        col[0].ce.o(0).dopelearn(-1)               // LTD
        nqaupd.v[nqcnum].x(nqaupd.v.size - 1) = -1 // put in nqaupd
        LTDCount += 1
      } else if(err2>0 && err2==err1) LTDCount += 1
      if(AdaptNoise && LTDCount >= StuckCount) { // if arm gets stuck, increase noise
        EMNoiseRate += EMNoiseRateInc
        print "now EMNoiseRate is ", EMNoiseRate
        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"
      }
    }
  }
}

//* resetPhase - check which movement phase (move or hold) and set plasticity params
proc resetPhase () { local md
  SUBPHASE = 0
  if(APHASE==AMOVE) { // phase was move, now set to hold/learn
    APHASE = AHOLD
    cvode.event(t+HoldDur,"resetPhase()")
    LearnON()
  } else if(APHASE==AHOLD) { // phase was hold/learn
    LearnOFF() // turn off learning
    if(DoReset) {
      APHASE = ARESET // now set to reset
      cvode.event(t+ResetDur,"resetPhase()")
    } else {
      APHASE = AMOVE // now set to move
      if(DoRDM==2) {
        md = MAXxy(aDT,aDT*int(ardm.negexp(MoveDur/aDT)))
        cvode.event(t+md,"resetPhase()")
      } else {
        cvode.event(t+MoveDur,"resetPhase()")
      }
    }
  } else if(APHASE==ARESET) { // phase waws reset, now set to move
    APHASE = AMOVE
    if(DoRDM==2) {
      md = MAXxy(aDT,aDT*int(ardm.negexp(MoveDur/aDT)))
      cvode.event(t+md,"resetPhase()")
    } else {
      cvode.event(t+MoveDur,"resetPhase()")
    }
  }
  print "t : ", t , " set phase to ", ATYP.o(APHASE).s
}

//* updateArm - update all arm apparatus
// This includes muscle commands, joints, muscle spindles, and reinforcement learning signals.
proc updateArm () { local errxy,errang,nqcnum,ii,tt
  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
    errang = getArmErr(ANGERR) // 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,APHASE,MLen[0],MLen[1],MLen[2],MLen[3],SUBPHASE,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) && (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()
    
    {nqaupd.tog() nqaupd.verbose=1}// Set the database back.

    setArmPos()         // Set the arm position.
    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 ((t - lastspdupdate >= spdDT) && (t >= DPlag)){//If it's time for a muscle spindle cell update...
    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)) {    
    RLUpdate()       // Do an RL update.    
    lastrlupdate = t // Remember the time of the event.
  }
  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.
}

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

  resetarm() // reset
  cvode.event(aDT,"updateArm()")
  if(DoLearn==4) {
    LearnON()
  } else LearnOFF()
  APHASE=AMOVE
  if(DoRDM) ardm.ACG(aseed)
  if(DoLearn==1) cvode.event(MoveDur,"resetPhase()")    
  SUBPHASE=0
  if(syDT>0) {nqsdel(nqsy) nqsy=new NQS("id1","id2","wg","t") cvode.event(syDT,"svsywts()")}

  LTDCount = lastmovetime = lastspdupdate = lastrlupdate = lastsysave = 0
  if(ResetEMNoiseRate && EMNoiseRate!=sgrhzEE) SetEMNoiseRate(EMNoiseRate=sgrhzEE)
}

//* mkaid - setup handlers for EM readout -> arm motion
proc mkaid () {
  aid = new FInitializeHandler(1,"initArmCB()") 
}

//* drtarg - draw target location
proc drtarg () { local xsz,clr
  xsz = 0.15
  if(numarg()>0) clr=$1 else clr=1
  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,ex,ey,wx,wy,ang0,ang1,ers,ln,xsz,she,shf,ele,elf localobj nqa,s
  ers=1 xsz=0.15
  ln=armLen[0]+armLen[1]
  g.size(-ln,ln,-ln,ln)
  s=new String()
  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(DoLearn) {
      sprint(s.s,"t=%g: %g %g %s",nqa.v[0].x(idx),ang0,ang1,ATYP.o(nqa.v[11].x(idx)).s)
    } else {
      sprint(s.s,"t=%g: %g %g",nqa.v[0].x(idx),ang0,ang1)
    }
  } else {
    if(numarg()>=1)ers=$1
    ex = ePos.x
    ey = ePos.y
    wx = armPos.x
    wy = 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(0,0,ex,ey,g,1,4) // draw arm
  drline(ex,ey,wx,wy,g,9,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
  {gvt=gvmarkflag gvmarkflag=0 g.erase if(nqa==nil) return}
  {rotArmTo(tAng[0],tAng[1]) 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)
}

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

  dr0 = deg2rad(tAng[0]) // convert to radians
  dr1 = deg2rad(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
proc setTargByID () { local tid
  tid=$1
  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(130,90)
  } else if (tid == 10) {
    setTarg(minang[0],minang[1])
  } else if (tid == 11) {
    setTarg(maxang[0],maxang[1])
  } 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
        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
        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=getarmX(deg2rad(s0),deg2rad(s1)) // starting x
    sy=getarmY(deg2rad(s0),deg2rad(s1)) // starting y
    err0 = sqrt( (sx-tPos.x)^2 + (sy-tPos.y)^2) // starting error
    
    nx=getarmX(deg2rad(s0+dirsh),deg2rad(s1+direl)) // new x
    ny=getarmY(deg2rad(s0+dirsh),deg2rad(s1+direl)) // new 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
  }
}

//* getnqrf([celltype,colid]) - make an NQS with approx celltype RFs -- default is for DP cells
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 
  tt=snq[cdx].getcol("t").x(0)
  for i=1,nqa.v.size-1 {
  if(i%10==0)printf(".")
    while(tt < nqa.getcol("t").x(i) && tdx < snq[cdx].v.size) {
      if(snq[cdx].getcol("type").x(tdx)==ct) {
        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)
  }
}

//
// Analysis functions (adding these here to avoid putting them always in notebook)
//

//* drerrRL - draw error and reinforcement signal
proc drerrRL () { local gvt
  gvt=gvmarkflag
  gvmarkflag=1
  nqaupd.gr("reinf","t",0,2,12) // reinforcement signal in red
  nqa.gr("err","t",0,3,12) // cartesian error in blue
  gvmarkflag=0
  nqaupd.gr("reinf","t",0,2,4) // reinforcement signal in red
  nqa.gr("err","t",0,3,4) // cartesian error in blue
  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],t,tAng[1],g,3,2)
}

//* 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],t,tAng[0],g,2,2)
}

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

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

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

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

//* DPatang(shoulder angle,elbow angle) - get Vector of DP cell IDs activated at the particular angles
obfunc DPatang () { local i,a0,a1 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
  rotArmTo(a0,a1) // rotate arm to target angle so can check muscle length activations
  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 MCMDatang ()  { local a,id1,id2,i,a0,a1,fctr localobj vidm,vcmd,vem
  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
  if(rotNorm) { // normalize rotation dx,dy by distance from origin
    if(armPos.x || armPos.y) fctr = 1 / sqrt(armPos.x^2 + armPos.y^2) else fctr = 1
    vcmd.x(0) = rotfctr[0] * fctr * (vcmd.x(1) - vcmd.x(0)) // should use a sigmoid (?)
    vcmd.x(1) = rotfctr[1] * fctr * (vcmd.x(3) - vcmd.x(2))
  } else {
    vcmd.x(0) = rotfctr[0] * (vcmd.x(1) - vcmd.x(0)) // should use a sigmoid (?)
    vcmd.x(1) = rotfctr[1] * (vcmd.x(3) - vcmd.x(2))
  }
  vcmd.resize(2)
  return vcmd
}

//* MCMDmapnq([angle increment]) - get an NQS with motor command for each angle pair
// just uses static connectivity from DP -> ES -> EM
obfunc MCMDmapnq () { local a0,a1,inc,x0,y0,x1,y1 localobj vcmd,nq
  if(numarg()>0) inc=$1 else inc=5
  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) {
      vcmd=MCMDatang(a0,a1)
      x0=armPos.x
      y0=armPos.y
      rotArm(vcmd.x(0),vcmd.x(1))
      x1=armPos.x
      y1=armPos.y
      nq.append(a0,a1,x0,y0,x1,y1,vcmd.x(0),vcmd.x(1))
    }
  }
  return nq
}

//* drMCMDmapnq(nqs from MCMDmapnq) - draws the NQS
proc drMCMDmapnq () { local x0,y0,x1,y1,i localobj nq
  nq=$o1
  nq.verbose=0
  nq.tog("DB")
  gvmarkflag=1
  nq.gr("y0","x0",0,4,8)
  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,8)
  nq.verbose=1
}

//* 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 {
    print "CDX:",CDX
    mksnq()
    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)
  }
}

//* 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 nsl=col.cstim.nsl
  ls=new List()
  for case(&sy,AM2,NM2,GA,GA2) ls.append(new List())
  for ltr(nc,ncl,&i) {
    if(nc.syn.type == EM) {
      if(nc.weight[AM2]>0) {
        ls.o(0).append(nsl.o(i))
      } 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
  rate=$1 if(numarg()>1) idx=$2 else idx=0  
  // print "rate is ", rate, "idx is ", idx
  if(lem==nil) lem=LEMNoise()
  for ltr(ns,lem.o(idx)) {
    if(rate <= 0) {
      ns.number = 0
    } else {
      ns.interval = 1e3 / rate
      ns.number = tstop * rate / 1e3
    }
  }
}

//* func calls

mkmyTYP()
initarm()
assignEM()
recE()
mkaid()
setTargByID(targid)
// if(DoAnim) gg()
LearnOFF() // start with no learning