// $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)
//print MLen[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(140)
// 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 arm length
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 arm length
} 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 arm length
} 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 arm length
} 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 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) // commented out because got error when running main.hoc by itself (filestem not found) -- fix!!
//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