-
Notifications
You must be signed in to change notification settings - Fork 2
/
msarm.hoc
2527 lines (2256 loc) · 94.1 KB
/
msarm.hoc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
// $Id: msarm.hoc,v 1.252 2012/07/21 21:55:29 samn Exp $
// Last revision: 5/25/13 (sald)
// Modification of the code for simple virtual 2D arm to interface with musculoskeletal 2D arm
// Set up Python object to allow interface with arminterface.py module.
objref pyobj
pyobj = new PythonObject()
objref gsyn // for graph with synaptic changes
objref ls // for power spectra
declare("verbosearm",0)
//* templates - need to write
begintemplate p2d
double x[1],y[1]
public x,y
proc init () {
if(numarg()==2) {x=$1 y=$2}
}
endtemplate p2d
//* variables/declares
// arm variables
declare("aid","o[1]") // finitializehandler for arm
declare("nqaupd","o[1]") // NQS for arm updates (not all involving position changes)
declare("nqa","o[1]") // NQS for recording arm joint angles/positions in time
declare("nqMCMD","o[1]") // NQS for recording motor commands during sim (to draw figure)
declare("armAng","d[2]","ePos",new p2d(), "armPos",new p2d(),"tPos",new p2d()) // arm angles; elbow position; end-effector position; target position
declare("tAng","d[2]") // target angle
declare("sAng","d[2]") // starting angle
declare("sPos","d[2]") // starting angle
declare("minang","d[2]","maxang","d[2]") // min/max angles for each joint
declare("armLen","d[2]") // length of each arm segment -- fixed throughout sim
// muscle variables
declare("nqE","o[1]") // has E assignment to muscle groups
declare("nqDP","o[1]") // has DP assignment to muscle groups
declare("MTYP","o[1]") // has muscle names
declare("MLen","d[4]") // length of each muscle -- varies throughout sim
declare("minMLen","d[4]","maxMLen","d[4]") // min/max muscle lengths
declare("splitEM", 1) // whether to read only from half of the EM population
declare("DiffRDEM",0) // whether to take difference from last activity when reading out EM activity
declare("musExc","d[4]") // output muscle excitation
declare("vEMmax",30) // max num of spikes of output EM population (depends on mcmdspkwd) - used to calculate normalized muscle excitation (musExc)
declare("minDPrate",0.00001,"maxDPrate",100) // min and max firing rate of DP NSLOCs (tuned to different muscle lengths)
declare("DPnoise",0.0) // noise of NSLOC input to DP cells
declare("DPoverlap", 1) // whether to have overlap in the encoding of muscle lengths (~population coding)
// time variables
declare("aDT",10) // how often to update the whole arm apparatus
declare("amoveDT",10) // how often to update the arm's position
declare("mcmdspkwd",50) // motor command spike window; how wide to make the counting window (in ms)
declare("EMlag",50) // lag between EM commands (at neuromuscular junction) and arm update
declare("lastmovetime",0) // when was the last arm move (in ms)?
declare("spdDT",10) // how often to update the muscle spindles (DP cells)
declare("DPlag",50) // lag between arm movement and DP update
declare("lastspdupdate",0) // when was the last spindle update (in ms)?
declare("rlDT",50) // how often to check RL status
declare("lastrlupdate",0) // when was the last RL update (in ms)?
declare("muscleNoiseChangeDT",500) // how often to alternate noise to muscle groups
declare("randomMuscleDTmax",1000) // max random delay to alternate muscle groups (changes every time) (-1 = off = use fixed muscleNoiseChangeDT)
declare("lastMuscleNoiseChange",0) // when was the last RL update (in ms)?
declare("syDT",0) // dt for recording synaptic weights into nqsy -- only used when syDT>0
declare("lastSynScaling", 0) // last time syns scaling occurred
declare("synScalingDT", 1000) // how often to apply synaptic scaling
// learning variables
//maxplastt_INTF6=rlDT + 50 // spike time difference over which to turn on eligibility
declare("minRLerrchangeLTP",0.001) // minimum error change visible to RL algorithm for LTP (units in Cartesian coordinates)
declare("minRLerrchangeLTD",0.001) // minimum error change visible to RL algorithm for LTD (units in Cartesian coordinates)
declare("DoLearn",4,"DoReset",2,"DoRDM",2, "DoAnim",0) // learn; reset; use random targets for training; animate arm
declare("invertAxis",0) // invert axis when plotting arm trajectory and target
declare("RLMode",3) // reinforcement learning mode (0=none,1=reward,2=punishment,3=reward+punishment)
declare("targid",2) // target ID for target to set up
declare("XYERR",0) // whether to use cartesian error
declare("ANGERR",1) // whether to use diff btwn targ and arm angle for error
declare("TRAJERR",2) // whether to use diff btwn targ and arm angle for error
declare("COMBERR",0) // whether to use diff btwn targ and arm angle for error
declare("errTY",XYERR) // type of error - either XYERR or ANGERR
declare("centerOutTask", 1) // use targets from center-out reaching task (1->L=2*monkey; 2->L=1.5*monkey)
// declare("AdaptLearn",0) // whether to modulate learning level by distance from target
// recording/visualization variables
declare("lrec","o[2]") // recording from ES, EM
declare("vEM",new Vector())
declare("pnq","o[2]") // NQS for motor command map
declare("nqsy","o[2]") // for recording synaptic weights over time -- only saves L5 cell synapses
declare("lastsysave",0) // when was the nqsy synaptic weight recording (in ms)?
sprint(tstr,"o[%d][%d]",CTYPi+1,CTYPi+1)
declare("lssyavg",tstr) // list of average synaptic weights vs time for a particular population
declare("syCTYP",1) // 1=only ES->EM ; 2 = all connections
// babble noise variables
declare("lem","o[1]") // List of EM AM2 'noise' NetStims to allow modulation
declare("AdaptNoise",0) // whether to adapt noise
declare("LTDCount",0) // number of recent LTD periods - used for noise adaptation
declare("StuckCount",2) // number of periods where arm doesn't improve and should adapt noise
declare("EMNoiseRate",250)//sgrhzEE) // rate of noise inputs to EM cells
declare("EMNoiseRateInc",50) // rate by which to increase noise rate when arm gets stuck
declare("EMNoiseRateDec",25) // rate by which to decrease noise rate when arm gets stuck
declare("ResetEMNoiseRate",0) // reset EMNoiseRate to sgrhzEE @ start of run ?
declare("EMNoiseRateMax",3e3) // rate of noise inputs to EM cells
declare("EMNoiseRateMin",50) // rate of noise inputs to EM cells
declare("EMNoiseMuscleGroup", 0) // alternate muscle group (-1=no alternation; initial muscle: 0=shext; 1=shflex ; 2=elext; 3=elflex) or pattern in expSeq
declare("exploreTot", 32) // max number of exploratory sequences to use from expSeq; if =4, activate only 1 muscle at a time
declare("expSeq","d[32]") // sequence of muscles activated during exploratory movements (in pairs)
/* --------------------------------------------------*/
/* DP/NSLOC FUNCTIONS */
/* ---------------------------------------------------*/
//* assignDP - set DP cell muscle length range responsiveness. store info in nqDP
proc assignDP () { local ii,jj,shmlenrangewid,elmlenrangewid,shminmlen,shmaxmlen,elminmlen,elmaxmlen localobj xo
{nqsdel(nqDP) nqDP = new NQS("col","id","ty","mid","mids","MLmin","MLmax") nqDP.strdec("mids")}
// Calculate the width of each subrange.
shmlenrangewid = (maxMLen[0]-minMLen[0]) / (cedp.count() / NMUSCLES) // muscle length range / (total num of dp units / number of muscles)
elmlenrangewid = (maxMLen[2]-minMLen[2]) / (cedp.count() / NMUSCLES)
// Initialize the subrange bounds to the first subrange.
shminmlen = minMLen[0]
shmaxmlen = minMLen[0]+shmlenrangewid
elminmlen = minMLen[2]
elmaxmlen = minMLen[2]+elmlenrangewid
// Loop over the different subranges...
for ii=0,cedp.count() / NMUSCLES - 1 {
// For the shoulder cells...
for jj=0,1 {
xo=cedp.o(ii*NMUSCLES+jj)
xo.start = 1 // if set to -1 they never spike; set to >0 to start spiking
xo.number = maxDPrate*tstop // set to max spikes (eg. 1000 Hz * 30 sec)
xo.noise = DPnoise // noise level - if 0 spike times are deterministic
xo.interval = 1/minDPrate // set base firing rate
// Set the muscle length range in the particular cell.
if (DPoverlap) {
xo.mlenmin = shminmlen-shmlenrangewid
xo.mlenmax = shmaxmlen+shmlenrangewid
// Add the muscle length bounds to the DP NQS table.
nqDP.append(1,xo.id,DP,xo.zloc,MTYP.o(xo.zloc).s,shminmlen-shmlenrangewid,shmaxmlen+shmlenrangewid)
} else {
xo.mlenmin = shminmlen
xo.mlenmax = shmaxmlen
// Add the muscle length bounds to the DP NQS table.
nqDP.append(1,xo.id,DP,xo.zloc,MTYP.o(xo.zloc).s,shminmlen,shmaxmlen)
}
}
// For the elbow cells...
for jj=2,3 {
xo=cedp.o(ii*NMUSCLES+jj)
xo.start = 1 // if set to -1 they never spike; set to >0 to start spiking
xo.number = maxDPrate*tstop // set to max spikes (eg. 1000 Hz * 30 sec)
xo.noise = DPnoise // noise level - if 0 spike times are deterministic
xo.interval = 1/minDPrate // set base firing rate
if (DPoverlap) {
xo.mlenmin = elminmlen-elmlenrangewid
xo.mlenmax = elmaxmlen+elmlenrangewid
// Add the muscle length bounds to the DP NQS table.
nqDP.append(1,xo.id,DP,xo.zloc,MTYP.o(xo.zloc).s,elminmlen-elmlenrangewid,elmaxmlen+elmlenrangewid)
} else {
// Set the muscle length range in the particular cell.
xo.mlenmin = elminmlen
xo.mlenmax = elmaxmlen
// Add the muscle length bounds to the DP NQS table.
nqDP.append(1,xo.id,DP,xo.zloc,MTYP.o(xo.zloc).s,elminmlen,elmaxmlen)
}
}
// Move to the next subrange.
shminmlen += shmlenrangewid
shmaxmlen += shmlenrangewid
elminmlen += elmlenrangewid
elmaxmlen += elmlenrangewid
}
}
//* updateDP - update the DP cell drive
proc updateDP () { local mlen, mlentime,nqarow,nqcnum localobj xo
// The first argument is the time at which to read muscle lengths. If none
// is provided, the default is the current time.
if (numarg() > 0) {
mlentime = $1
nqa.verbose = 0
nqarow = nqa.select("t","<=",mlentime) - 1
nqa.tog()
nqa.verbose = 1
} else mlentime = t
nqcnum = nqaupd.fi("spupd") //Remember that we update the DP cells in nqaupd.
nqaupd.v[nqcnum].x(nqaupd.v.size-1) = 1
// Read arm angles from Python interface to MSM
for i = 0, 1 {
armAng[i] = pyobj.axf.getAngleReceived(i)
}
// Read muscle lengths from Python interface to MSM
for i = 0, 3 {
MLen[i] = pyobj.axf.getLengthReceived(i)
}
// Update the interval (firing rate) of the DP/NSLOC units according to muscle length tuning curve
nqcnum = nqa.fi("ML0")
for ltr(xo,cedp) {
// select mlen value as a function of t and zloc
if (mlentime == t) { // use current muscle length
mlen = MLen[xo.zloc]
} else { // use muscle length stored in nqa (previous t)
mlen = nqa.v[nqcnum+xo.zloc].x(nqarow)
}
// update value of DP NSLOCs based on whether mlen falls in its range [mlenmin, mlenmax}
if ((mlen >= xo.mlenmin) && (mlen <= xo.mlenmax)) {
xo.interval = 1000.0/maxDPrate // interval = inverse of max firing rate (in ms)
} else {
xo.interval = 1000.0/minDPrate // interval = inverse of min firing rate (in ms)
}
/*if (t <600) {
xo.interval = 1000.0/maxDPrate
} else if (t >600 && t<1200) {
xo.interval = 1000.0/minDPrate
} else if (t >1200) {
xo.interval = 1000.0/maxDPrate
}*/
}
}
/* --------------------------------------------------*/
/* EM / ES FUNCTIONS */
/* ---------------------------------------------------*/
//* GetEMType(cell id) - returns type of motor unit - only used for EM cells
// lookup in MTYP List to see string representation
func GetEMType () {
if($1 < col.ix[EM] || $1 > col.ixe[EM]) return -1 // not an EM cell? return -1
return $1 % 4 // otherwise, return type
}
//* assignEM() - create nq with list of EM cells and associated muscle
proc assignEM () { local i,ct,mid,a localobj vty
{nqsdel(nqE) nqE = new NQS("col","id","ty","mid","mids") nqE.strdec("mids")}
a=allocvecs(vty)
vty.append(EM)
if (splitEM) {
nqE.clear(col[0].numc[EM]/2)
for vtr(&ct,vty) for i=col[0].ix[ct]+col[0].numc[EM]/2,col[0].ixe[ct] {
mid = i%4
nqE.append(1,i,ct,mid,MTYP.o(mid).s)
}
} else {
nqE.clear(col[0].numc[EM])
for vtr(&ct,vty) for i=col[0].ix[ct],col[0].ixe[ct] {
mid = i%4
nqE.append(1,i,ct,mid,MTYP.o(mid).s)
}
}
dealloc(a)
}
//* readoutEM - read activity from EM - store results in vEM
proc readoutEM () { local i,idx,ldx,sz
vEM.resize(MTYP.count) // resize to the number of muscle groups
vEM.fill(0)
ldx=1 // column?
for i=0,nqE.v.size-1 { // for all EM units
idx=nqE.v[1].x(i)-col[0].ix[EM] // obtain the relative idx value
vEM.x(nqE.v[3].x(i)) += lrec[ldx].o(1).o(idx).size // add all recorded spikes from that unit to corresponding count
lrec[ldx].o(1).o(idx).resize(0) // reset to 0
}
if(DiffRDEM) {
sz=nqa.v.size
if(sz) {
vEM.x(0) -= nqa.v[5].x(sz-1)
vEM.x(1) -= nqa.v[6].x(sz-1)
vEM.x(2) -= nqa.v[7].x(sz-1)
vEM.x(3) -= nqa.v[8].x(sz-1)
for i=0,3 if(vEM.x(i) < 0) vEM.x(i) = 0
}
}
}
//* recE - set up recording from EM cells
proc recE () {
//lrec[0] = mkrecl(col[0],ES)
lrec[1] = mkrecl(col[0],EM)
}
//* svsywts - save synaptic weights from ES cells
proc svsywts () { local i,a localobj xo,vwg,vtau,vinc,vmaxw,vidx,vt,vpre
// print "svsywts: t " , t
a=allocvecs(vidx,vwg,vtau,vinc,vmaxw,vt,vpre)
if (syCTYP == 1) { //record only ES synapses
for ltr(xo,col[0].ce) {
if(xo.type==ES) {
vrsz(xo.getdvi(vidx),vwg,vtau,vinc,vmaxw,vt,vpre)
xo.getplast(vwg,vtau,vinc,vmaxw) // get current weight gains
vt.fill(t) //time
vpre.fill(xo.id) //presynaptic id
nqsy.v[0].append(vpre)
nqsy.v[1].append(vidx)
nqsy.v[2].append(vwg)
nqsy.v[3].append(vt)
}
}
} else (syCTYP == 2) { //record all synapses
for ltr(xo,col[0].ce) {
vrsz(xo.getdvi(vidx),vwg,vtau,vinc,vmaxw,vt,vpre)
xo.getplast(vwg,vtau,vinc,vmaxw) // get current weight gains
vt.fill(t) //time
vpre.fill(xo.id) //presynaptic id
nqsy.v[0].append(vpre)
nqsy.v[1].append(vidx)
nqsy.v[2].append(vwg)
nqsy.v[3].append(vt)
}
}
dealloc(a)
}
proc synScaling () {local i,a,wtSumOriginal,wtSum,scaleFactor localobj xo,vidx,delay_vec,prob_vec,wt1vec,wt1vecAbs,wt2vec,wt2vecAbs,distalsyns,wgain,vwgain,vplasttau,vplastinc,vplastmaxw
// allocate vectors
a=allocvecs(vidx,delay_vec,prob_vec,wt1vec,wt1vecAbs,wt2vec,wt2vecAbs,distalsyns,wgain,vwgain,vplasttau,vplastinc,vplastmaxw)
wtSumOriginal=0
wtSum=0
scaleFactor=0
// loop over all cells
for ltr(xo,col[0].ce) {
// set vector length based on number of postsynaptic connections
vrsz(xo.getdvi(vidx),delay_vec,prob_vec,wt1vec,wt1vecAbs,wt2vec,wt2vecAbs,distalsyns,wgain)
// read all weight gains and vwt1- getdvi
xo.getdvi(vidx,delay_vec,prob_vec,wt1vec,wt2vec,distalsyns,wgain)
// get absolute weight of wt1
wt1vecAbs.abs(wt1vec)
// get absolute weight of wt1
wt2vecAbs.abs(wt2vec)
if (wgain.sum() > 0) {
// sum original weights
wtSumOriginal=wtSumOriginal+wt1vecAbs.sum()+wt2vecAbs.sum()
// sum current weights
wtSum=wtSum+wt1vecAbs.mul(wgain).sum()+wt2vecAbs.mul(wgain).sum()
}
}
// calculate normalizing constant (scaling factor)
scaleFactor = 1.0 * wtSumOriginal / wtSum
//print "Scaling Factor: ", scaleFactor
if (scaleFactor < 1) {
// set all weight gains via setplast
for ltr(xo,col[0].ce) {
//print "pre:", xo.gid
// set vector length based on number of postsynaptic connections
vrsz(xo.getdvi(vidx),vwgain,vplasttau,vplastinc,vplastmaxw)
// read wgains
xo.getplast(vwgain,vplasttau,vplastinc,vplastmaxw)
// scale synaptic weight gain
vwgain.mul(scaleFactor)
// update scaled wgains
xo.setplast(vwgain,vplasttau,vplastinc,vplastmaxw)
}
}
dealloc(a)
}
/* --------------------------------------------------*/
/* ARM/MUSCLE FUNCTIONS */
/* ---------------------------------------------------*/
//* initarm - init arm length, excitation, DP responsiveness, nrnpython axf - combine with initarmCB??
proc initarm () {
// arm dimensions - multiplied by 0.5 to match Alex's model dimensions
armLen[0] = 0.4634 - 0.173 // elbow - shoulder from MSM; previously, 0.5*(0.55) // shoulder to elbow (check if ms = wam length !)
if (centerOutTask >= 3) {
armLen[1] = 0.7169 - 0.4634 + 0.18 // radioulnar - elbow from MSM; previously, 0.5*(0.3+0.06+0.095) // elbow to center of hand = elbow to wrist + wrist length + hand length/2 (check if ms = wam legth)
} else {
armLen[1] = 0.7169 - 0.4634 // radioulnar - elbow from MSM; previously, 0.5*(0.3+0.06+0.095) // elbow to center of hand = elbow to wrist + wrist length + hand length/2 (check if ms = wam legth)
}
// min and max joint angles (in degrees)
minang[0] = deg2rad(-10) //shoulder
maxang[0] = deg2rad(110)
minang[1] = deg2rad(0) //elbow
maxang[1] = deg2rad(135)
// if using normalizedFiberLength
//minMLen[0] = 0.4 maxMLen[0] = 1.6 // define min and max muscle length values
//minMLen[1] = 0.4 maxMLen[1] = 1.6
//minMLen[2] = 0.4 maxMLen[2] = 1.6
//minMLen[3] = 0.4 maxMLen[3] = 1.6
// if using fiberLength
minMLen[0] = 0.03 maxMLen[0] = 0.23 // define min and max muscle length values
minMLen[1] = 0.03 maxMLen[1] = 0.23
minMLen[2] = 0.03 maxMLen[2] = 0.23
minMLen[3] = 0.03 maxMLen[3] = 0.23
musExc[0] = 0 // shoulder extensor (deltoid)
musExc[1] = 0 // shoulder flexor (pectoralis)
musExc[2] = 0 // elbow extensor (triceps)
musExc[3] = 0 // elbow flexor (biceps+brachialis)
// create exploratory pattern of movements as pairs of muscle groups that are activated together (-1=no muscle)
expSeq[0] = 0 expSeq[1] = -1
expSeq[2] = 1 expSeq[3] = -1
expSeq[4] = 2 expSeq[5] = -1
expSeq[6] = 3 expSeq[7] = -1
expSeq[8] = 0 expSeq[9] = 2
expSeq[10] = 1 expSeq[11] = 3
expSeq[12] = 0 expSeq[13] = 3
expSeq[14] = 1 expSeq[15] = 2
expSeq[16] = 3 expSeq[17] = -1
expSeq[18] = 1 expSeq[19] = -1
expSeq[20] = 2 expSeq[21] = -1
expSeq[22] = 0 expSeq[23] = -1
expSeq[24] = 1 expSeq[25] = 2
expSeq[26] = 1 expSeq[27] = 3
expSeq[28] = 0 expSeq[29] = 2
expSeq[30] = 0 expSeq[31] = 3
calculateCenterPos(sAng0,sAng1) // calculate center position (used for TRAJERR)
assignDP() // assign DP responsiveness, store info in nqDP
setTargByID(targid) // set target
// Initialize the socket interface to the virtual or real (robotic) arm.
nrnpython("import arminterface_pipe as axf")
// run setup() in python- pyobj.axf.setup()
sprint(tstr,"axf.setup(%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f)",tstop/1e3, aDT,sAng[0], sAng[1], tPos.x, tPos.y, damping, shExtGain, shFlexGain, elExtGain, elFlexGain)
nrnpython(tstr)
}
//* initArmCB - initialize nqaupd, nqa and arm callbacks
proc initArmCB () { local i
// Set up the muscle commands NQS table.
nqsdel(nqaupd)
nqaupd = new NQS("t","shext","shflex","elext","elflex","reinf","mvmt","spupd")
// Set up the arm NQS table.
{nqsdel(nqa) nqa=new NQS("t","ang0","ang1","x","y","shext","shflex","elext","elflex","ex","ey","phase")}
nqa.resize("ML0","ML1","ML2","ML3","subphase","errxy","errang")
// Clear out the current muscle command array vEM.
if(vEM==nil) vEM = new Vector(MTYP.count) else {vEM.resize(MTYP.count) vEM.fill(0)}
// initialize arm callback
cvode.event(aDT,"updateArm()")
// Do RL with dopamine
if(DoLearn==4) {
LearnON()
} else LearnOFF()
// create nqsy
if(syDT>0) {nqsdel(nqsy) nqsy=new NQS("id1","id2","wg","t") cvode.event(syDT,"svsywts()")}
// initialize last variables to 0
LTDCount = lastmovetime = lastspdupdate = lastrlupdate = lastsysave = lastmuscle = lastMuscleNoiseChange = lastSynScaling = 0
// reset EM noise rate
if(ResetEMNoiseRate && EMNoiseRate!=sgrhzEE) SetEMNoiseRate(EMNoiseRate=sgrhzEE)
}
//* mkmyTYP - set up names of muscle groups
proc mkmyTYP () {
MTYP=new List()
MTYP.append(new String("shext"))
MTYP.append(new String("shflex"))
MTYP.append(new String("elext"))
MTYP.append(new String("elflex"))
}
//* updateArm - update all arm apparatus
// This includes muscle commands, joints, muscle spindles, and reinforcement learning signals.
proc updateArm () { local errxy,errang,nqcnum,ii,tt localobj randomMuscleDT
readoutEM()//Read M1 muscle group command spike counts since updateArm() call. These go into vEM.
if (t == aDT) { // If we are at the beginning of the simulation...
nqaupd.append(0,0,0,0,0,0,0,0) // Set up an arm update entry for t=0.
// Set up an arm position entry for t=0.
errxy = getArmErr(XYERR) // gets error in position
if (errTY==ANGERR) {
errang = getArmErr(ANGERR) // gets error in angle
} else if (errTY==TRAJERR) {
errang = getArmErr(TRAJERR) // gets error in angle
}
nqa.append(0,armAng[0],armAng[1],armPos.x,armPos.y,vEM.x(0),vEM.x(1),vEM.x(2),vEM.x(3),ePos.x,ePos.y,0,MLen[0],MLen[1],MLen[2],MLen[3],0,errxy,errang)
updateDP() // Set up the initial DP cell activity.
}
// Add vEM to the muscle commands part of the arm update table, nqaupd.
nqaupd.append(t,vEM.x(0),vEM.x(1),vEM.x(2),vEM.x(3),0,-1,0)
// If it's time for another arm motion...
if ((t - lastmovetime >= amoveDT)) {
if (t >= mcmdspkwd + EMlag) {
// Get all of the motor command entries in a window mcmdspkwd ms wide, back
// EMlag ms from the current time t.
{nqaupd.verbose=0 nqaupd.select("t","[)",t-mcmdspkwd-EMlag,t-EMlag)}
// Set vEM to the sum of each of the relevant motor command entries.
vEM.x(0) = nqaupd.getcol("shext").sum()
vEM.x(1) = nqaupd.getcol("shflex").sum()
vEM.x(2) = nqaupd.getcol("elext").sum()
vEM.x(3) = nqaupd.getcol("elflex").sum()
musExc[0] = vEM.x(0) / vEMmax
musExc[1] = vEM.x(1) / vEMmax
musExc[2] = vEM.x(2) / vEMmax
musExc[3] = vEM.x(3) / vEMmax
// if not enough time to collect spikes from sim (mcmdspikwindow) then use 0 excitation
} else {
musExc[0] = 0
musExc[1] = 0
musExc[2] = 0
musExc[3] = 0
}
// Send a packet out to the virtual or real arm with neural excitations
sprint(tstr,"axf.sendAndReceiveDataPackets(%f,%f,%f,%f,%f)",t, musExc[0], musExc[1], musExc[2], musExc[3])
nrnpython(tstr)
// Update the arm position
ePos.x = armLen[0] * cos(armAng[0]) // end of elbow
ePos.y = armLen[0] * sin(armAng[0])
armPos.x = ePos.x + armLen[1] * cos(armAng[0] + armAng[1]) // wrist=arm position
armPos.y = ePos.y + armLen[1] * sin(armAng[0] + armAng[1])
// Get errors
errxy = getArmErr(XYERR) // gets error in position
if (errTY==ANGERR) {
errang = getArmErr(ANGERR) // gets error in angle
} else if (errTY==TRAJERR) {
errang = getArmErr(TRAJERR) // gets error in angle
}
// Set the database back.
{nqaupd.tog() nqaupd.verbose=1}
if (DoAnim) drarm() // Draw the arm configuration.
lastmovetime = t // Remember the time of the event.
// Add the nqa row number to the nqaupd table.
nqcnum = nqaupd.fi("mvmt")
nqaupd.v[nqcnum].x(nqaupd.v.size - 1) = nqa.v.size - 1
}
//If it's time for a muscle spindle cell update...
if ((t - lastspdupdate >= spdDT) && (t >= DPlag)){
updateDP(t-DPlag) // Update DP cell activity.
lastspdupdate = t // Remember the time of the event.
}
// If it's time for an RL update...
if ((DoLearn==4) && (t - lastrlupdate > rlDT)) {
// update EM muscle group that receives noise as input
if (EMNoiseMuscleGroup > -1 && (t-lastMuscleNoiseChange >= muscleNoiseChangeDT)) {
//print "muscle group=", EMNoiseMuscleGroup
lastMuscleNoiseChange = t
EMNoiseMuscleGroup += 2 // increase exploratory sequence index (seee expSeq)
if (EMNoiseMuscleGroup >= exploreTot) { // wrap around
EMNoiseMuscleGroup = 0
}
// Use random delays to alternate the muscle that receives noise during exploratory movements
if (randomMuscleDTmax>0) {
randomMuscleDT = new Random()
randomMuscleDT.ACG(inputseed+iepoch+t) //initialize random number generator
randomMuscleDT.uniform(0,1)
muscleNoiseChangeDT= randomMuscleDT.uniform(100,randomMuscleDTmax)
EMNoiseMuscleGroup = int(randomMuscleDT.uniform(0,exploreTot/2))*2
print "muscleNoiseChangeDT= ",muscleNoiseChangeDT
print "muscle group=", EMNoiseMuscleGroup
}
}
if (verbosearm) print "Babble noise EM Muscle Group: ", EMNoiseMuscleGroup
RLUpdate() // Do an RL update.
lastrlupdate = t // Remember the time of the event.
}
// If it's time for synaptic scaling
if ((DoLearn==4) && (t - lastSynScaling >= synScalingDT)) {
synScaling()
lastSynScaling = t
}
// Add a position entry to nqa.
//print "shAng: ", armAng[0]
//print "elAng: ", armAng[1]
nqa.append(t,armAng[0],armAng[1],armPos.x,armPos.y,vEM.x(0),vEM.x(1),vEM.x(2),vEM.x(3),ePos.x,ePos.y,0,MLen[0],MLen[1],MLen[2],MLen[3],0,errxy,errang)
if (syDT > 0 && t - lastsysave >= syDT) { // time for nqsy synaptic save
svsywts() // Do a synaptic save.
lastsysave = t // Remember the time of the event.
}
cvode.event(t+aDT,"updateArm()") // Post the next updateArm() event.
// IF end, call python function to close sockets, save data and plot arm graphs (pass sim length as argument)
if (t >= tstop-aDT/2) {
//sprint(tstr,"axf.closeSavePlot(%f, %f)",t/1000, aDT)
//nrnpython(tstr)
print "Finished sim..."
if (syDT > 0 && lastsysave < t) { // save synaptic weights in last step
svsywts() // Do a synaptic save.
}
}
}
//* mkaid - setup handlers for periodic arm update
proc mkaid () {
aid = new FInitializeHandler(1,"initArmCB()")
}
/* --------------------------------------------------*/
/* LEARNING/NOISE FUNCTIONS */
/* ---------------------------------------------------*/
//* getArmErr([errtype]) - gets current error in position
// 1st arg is optional - if left out, uses global errTY
func getArmErr () { local dsq,ety, xtmp, ytmp, rotAng, xrot,yrot localobj tAngs
if(numarg()>0) ety=$1 else ety=errTY
if(ety == XYERR) {
// error as Cartesian distance between hand position and target position
return sqrt((armPos.x - tPos.x)^2 + (armPos.y - tPos.y)^2)
} else if(ety == ANGERR) {
// error as squared difference between joint and target angles
return sqrt((armAng[0] - tAng[0])^2 + (armAng[1] - tAng[1])^2)
} else if(ety == TRAJERR) {
// error as squared difference between current position and ideal (straight line) trajectory line
xtmp = armPos.x-sPos[0]
ytmp = armPos.y-sPos[1]
//print "armPos.x:",armPos.x,"ArmPos.y",armPos.y
//print "xtemp:",xtmp,"ytemp:",ytmp
//set angle of rotation as a function of target
tAngs = new Vector()
tAngs.append(0, -PI, -PI/2, PI/2, -PI/4, -3*PI/4, PI/4, 3*PI/4)
rotAng=tAngs.x[targid]
// rotate trajectory to position over positive x-axis (=target 0, center-right)
//xrot = xtmp*cos(rotAng) - ytmp*sin(rotAng)
yrot = xtmp*sin(rotAng) + ytmp*cos(rotAng)
//print "yrot(error):", yrot
// retrun deviation from ideal trajectory
return abs(yrot)
}
}
//* LearnON - turn on learning
proc LearnON () {
plaststartT_INTF6 = 0
plastendT_INTF6 = tstop * 2
}
//* LearnOFF - turn off learning
proc LearnOFF () {
plaststartT_INTF6 = tstop * 2
plastendT_INTF6 = tstop * 3
}
//* RLUpdate - reinforcement learning updateupdate - MODIFY
proc RLUpdate () { local err1xy,err2xy,err1ang,err2ang,err1,err2,nqcnum,combRL
// If we have at least 2 nqa entries...
if (nqa.v.size > 1) {
// Find the distance error of the 2nd-to-last hand position (in
// nqa) from the target hand position.
err1xy = nqa.getcol("errxy").x(nqa.v.size-2)
// Find the distance error of the last hand position (in
// nqa) from the target hand position.
err2xy = nqa.getcol("errxy").x(nqa.v.size-1)
// Find the distance error of the 2nd-to-last hand position (in
// nqa) from the target hand position.
err1ang = nqa.getcol("errang").x(nqa.v.size-2)
// Find the distance error of the last hand position (in
// nqa) from the target hand position.
err2ang = nqa.getcol("errang").x(nqa.v.size-1)
if(XYERR == errTY) {
err1 = err1xy
err2 = err2xy
} else if (ANGERR == errTY) {
err1 = err1ang
err2 = err2ang
} else if (TRAJERR == errTY) {
err1 = err1ang // store traj error in same variable as ang
err2 = err2ang
}
if(verbosearm) print "err1 : ", err1, " err2 " , err2
if ((DoLearn == 4) && (RLMode > 0)) { // RL with dopamine
combRL=0
if(COMBERR==4) {
if( err1xy - err2xy >= minRLerrchangeLTP && err1ang >= err2ang ) { //dev
combRL = 1 // do LTP
} else if( err2xy - err1xy >= minRLerrchangeLTD && err1ang < err2ang ) { //deviation change
combRL = 2 // do LTD
}
} else if (COMBERR==2) {
if( err1xy - err2xy >= minRLerrchangeLTP && abs(err1ang) < 0.03) { // midev
combRL = 1 // do LTP
} else if( err2xy - err1xy >= minRLerrchangeLTD && abs(err1ang) > 0.03) { // min deviation
combRL = 2 // do LTD
}
} else if (COMBERR==3) {
if( err1xy - err2xy >= minRLerrchangeLTP && abs(err1ang) < 0.02) { // midev2
combRL = 1 // do LTP
} else if( err1ang<err2ang && abs(err1ang) > 0.02) { // min deviation2
combRL = 2 // do LTD
}
} else if (COMBERR==1) { // alternate XYERR and TRAJERR
if (errTY == XYERR) {
errTY == TRAJERR
} else if (errTY = TRAJERR) {
errTY == XYERR
}
}
nqcnum = nqaupd.fi("reinf")
if (verbosearm) {print "err1: ", err1, " err2: " , err2, " err1-err2: ", err1-err2, " minRLerrchangeLTP: ", minRLerrchangeLTP, " minRLerrchangeLTD: ", minRLerrchangeLTD, " LTDCount: ", LTDCount}
// Reward
if (((err1 - err2 >= minRLerrchangeLTP) && ((RLMode == 1) || (RLMode == 3)) && COMBERR==0) || combRL==1) {
// if(AdaptLearn) EPOTW_INTF6 = (err1 - err2) / err1 // else EPOTW_INTF6=1
col[0].ce.o(0).dopelearn(1) // LTP
if (verbosearm) print "LTP!"
nqaupd.v[nqcnum].x(nqaupd.v.size - 1) = 1 // put in nqaupd under "reinf"
LTDCount = 0
// print "LTDCount is now 0 , A"
if(AdaptNoise) { // drop noise towards baseline
EMNoiseRate -= EMNoiseRateDec
if(EMNoiseRate < EMNoiseRateMin) EMNoiseRate = EMNoiseRateMin
if(verbosearm) print "setting EMNoiseRate to ", EMNoiseRate
//SetEMNoiseRate(EMNoiseRate)
}
// Punishment
} else if (combRL==2 || ((err2 - err1 >= minRLerrchangeLTD) && ((RLMode == 2) || (RLMode == 3)) && COMBERR==0)) {
// if(AdaptLearn) EDEPW_INTF6 = (err2 - err1) / err1 // else EDEPW_INTF6=1
col[0].ce.o(0).dopelearn(-1) // LTD
if (verbosearm) print "LTD!"
nqaupd.v[nqcnum].x(nqaupd.v.size - 1) = -1 // put in nqaupd
LTDCount += 1
} else if(err2>0 && err2==err1) LTDCount += 1
// Adapt noise
if(AdaptNoise && LTDCount >= StuckCount) { // if arm gets stuck, increase noise
EMNoiseRate += EMNoiseRateInc
if(EMNoiseRate > EMNoiseRateMax) EMNoiseRate = EMNoiseRateMax
if(verbosearm) print "setting EMNoiseRate to ", EMNoiseRate
//SetEMNoiseRate(EMNoiseRate)
LTDCount = 0 // reset counter to 0
// print "LTDCount is now 0 , B"
}
SetEMNoiseRate(EMNoiseRate)
if (verbosearm) print "setting EMNoiseRate to ", EMNoiseRate
}
}
}
//* get a List of Lists with the noise NetStims for EM cells
obfunc LEMNoise () { local i,sy localobj nc,ncl,nsl,ls
ncl=col.cstim.ncl // list of netconnections = lcstim.o(0).ncl
nsl=col.cstim.nsl // list of netstims = lcstim.o(0).nsl
ls=new List() // list of lists with noise netstims
for case(&sy,AM2,NM2,GA,GA2) ls.append(new List()) // for each synapse type add new list
for ltr(nc,ncl,&i) { // iterate through netconnections
if(nc.syn.type == EM) { // if EM
if(nc.weight[AM2]>0) {
ls.o(0).append(nsl.o(i)) // add Netstim to list
} else if(nc.weight[NM2]>0) {
ls.o(1).append(nsl.o(i))
} else if(nc.weight[GA]>0) {
ls.o(2).append(nsl.o(i))
} else if(nc.weight[GA2]>0) {
ls.o(3).append(nsl.o(i))
}
}
}
return ls
}
//* SetEMNoiseRate(rate,index) - set rate of EM cell noise
// index: 0==AM2, 1==NM2, 2==GA, 3==GA2
proc SetEMNoiseRate () { local i,rate,idx localobj ns
//splitEM=0
rate=$1 if(numarg()>1) idx=$2 else idx=0
//ncl=col.cstim.ncl // list of netconnections = lcstim.o(0).ncl
if (verbosearm) print "rate is ", rate, "idx is ", idx
if(lem==nil) lem=LEMNoise() // get list of lists with NetStims
for ltr(ns,lem.o(idx), &i) { // iterate through netstims
if ((splitEM && i >= col[0].numc[EM]/2) || splitEM==0) { // if splitEM only apply noise to half of EM
isplit = i
if (splitEM) {isplit = i-col[0].numc[EM]/2}
if (EMNoiseMuscleGroup == -1 || (nqE.v(3).x(isplit) == expSeq[EMNoiseMuscleGroup]) || (nqE.v(3).x(isplit) == expSeq[EMNoiseMuscleGroup+1])) { // compare muscle id of cell with current muscle group
if(rate <= 0) {
ns.number = 0
} else {
ns.interval = 1e3 / rate
ns.number = tstop * rate / 1e3
//print "on:", i , " interval:", ns.interval
}
} else { // set to 0 the rest of muscle groups
ns.interval = 1e3 / EMNoiseRateMin
//print "off:", i, " interval:", ns.interval
}
}
}
//splitEM=1
}
/* --------------------------------------------------*/
/* DRAW ARM FUNCTIONS */
/* ---------------------------------------------------*/
//* drtarg - draw target location
proc drtarg () { local xsz,clr
xsz = 0.05
if(numarg()>0) clr=$1 else clr=1
if (invertAxis) {
drline(-tPos.x-xsz,-tPos.y-xsz,-tPos.x+xsz,-tPos.y+xsz,g,clr,4) // draw an x for the target
drline(-tPos.x-xsz,-tPos.y+xsz,-tPos.x+xsz,-tPos.y-xsz,g,clr,4)
} else {
drline(tPos.x-xsz, tPos.y-xsz, tPos.x+xsz, tPos.y+xsz,g,clr,4) // draw an x for the target
drline(tPos.x-xsz, tPos.y+xsz, tPos.x+xsz, tPos.y-xsz,g,clr,4)
}
}
//* drarm([nqa,row from nqa,erase]) - draw arm location
proc drarm () { local idx,shx, shy,ex,ey,wx,wy,ang0,ang1,ers,ln,xsz,she,shf,ele,elf localobj nqa,s
ers=1 xsz=0.15
shx = 0 //0.173 ; use shoulder joint as reference
shy = 0
ln=armLen[0]+armLen[1]
g.size(-ln,ln,-ln,ln)
s=new String()
// update using data from nqa
if(numarg()>=2) {
{nqa=$o1 idx=$2}
if(idx < 0 || idx >= nqa.v.size) {printf("drarm ERRA: invalid index: %d,%d\n",idx,nqa.v.size) return}
ex = nqa.v[9].x(idx)
ey = nqa.v[10].x(idx)
wx = nqa.v[3].x(idx)
wy = nqa.v[4].x(idx)
ang0 = nqa.v[1].x(idx)
ang1 = nqa.v[2].x(idx)
if(numarg()>=3) ers=$3 // if erase argument present
if(DoLearn) {
sprint(s.s,"t=%g: %g %g %s",nqa.v[0].x(idx),ang0,ang1)
} else {
sprint(s.s,"t=%g: %g %g",nqa.v[0].x(idx),ang0,ang1)
}
// if updated from current ePos and armPos
} else {
if(numarg()>=1) ers=$1 // if erase argument present
if (invertAxis) {
ex = -(shx+ePos.x)
ey = -(shy+ePos.y)
wx = -(shx+armPos.x)
wy = -(shy+armPos.y)
} else {
ex = shx+ePos.x
ey = shy+ePos.y
wx = shx+armPos.x
wy = shy+armPos.y
}
ang0 = armAng[0]
ang1 = armAng[1]
if(0 && DoLearn) {
sprint(s.s,"t=%g: %g %g %s",t,ang0,ang1) //,ATYP.o(APHASE).s)
} else {
sprint(s.s,"t=%g: %g %g",t,ang0,ang1)
}
}
// if(verbose) print "ex=",ex,"ey=",ey,"ang0=",ang0,"wx=",wx,"wy=",wy,"ang1=",ang1
if(ers) g.erase_all()
drtarg() // draw an x for the target
drline(shx,shy,ex,ey,g,2,4) // draw arm
drline(ex,ey,wx,wy,g,3,4)
if(DoAnim > 1) {
if(NMUSCLES==2) { // draw the muscle lengths
drline(3.1,0,3.1,MLen[0],g,2,3)
drline(3.2,0,3.2,MLen[1],g,3,3)
} else {
drline(3.1,0,3.1,MLen[0],g,2,3)
drline(3.2,0,3.2,MLen[2],g,3,3)
}
}
// draw motor commands
if(DoAnim > 2 && nqaupd!=nil) {
{nqaupd.verbose=0 nqaupd.select("t","[)",t-mcmdspkwd-EMlag,t-EMlag)}
she = nqaupd.getcol("shext").sum() / ( col.numc[EM] / 4 )
shf = nqaupd.getcol("shflex").sum() / ( col.numc[EM] / 4 )
ele = nqaupd.getcol("elext").sum() / ( col.numc[EM] / 4 )
elf = nqaupd.getcol("elflex").sum() / ( col.numc[EM] / 4 )
drline(0,-3.1,she,-3.1,g,2,3)
drline(0,-3.1,-shf,-3.1,g,2,3)
drline(0,-3.2,ele,-3.2,g,3,3)
drline(0,-3.2,-elf,-3.2,g,3,3)
{nqaupd.tog("DB") nqaupd.verbose=1}
}
g.label(0.55,0.95,s.s)
g.flush()
doNotify()
}
//* drxytraj - draws the x,y position from nqa
proc drxytraj () { local gvt
g=new Graph()
g.size(0,tstop,0,2)
{gvt=gvmarkflag gvmarkflag=0 g.erase if(nqa==nil) return}
drarm()
nqa.gr("y","x",0,2,1)
g.exec_menu("View = plot")
gvmarkflag=gvt
}
//* nqa2gif(nqa[,inc]) - saves arm locations in nqa as gif
// inc specifies how many frames to skip
proc nqa2gif () { local i,j,inc localobj nqa,s
if(!FileExists("gif/wg")) system("mkdir gif/wg")
nqa=$o1 s=new String() j=0
if(numarg()>1) inc=$2 else inc=1
for(i=0;i<nqa.v.size;i+=inc){
drarm(nqa,i)
sprint(s.s,"xcalc2gif gif/wg/%010d_nqa.gif",j)
system(s.s)
j+=1
}
}
//* animnqa(nqa[,startidx,endidx,delays]) - animates contents of nqa (arm position info)
proc animnqa () { local i,is,ie,del localobj nqa
nqa=$o1
if(numarg()>1) is=$2 else is=0
if(numarg()>2) ie=$3 else ie=nqa.v.size-1
if(numarg()>3) del=$4 else del=0.25e9
for i=is,ie {
drarm(nqa,i)
sleepfor(0,del)
}
}
//* whirlgif(outputfile,framesglob[,delframes]) - makes a moving gif using whirlgif
// delframes specifies whether to delete the gif frames after
func whirlgif () { local del localobj s
if(!FileExists("/usr/site/pkgs/download/whirlgif304/whirlgif")) {
print "whirlgif ERR0: can't find whirlgif binary!"
return 0
}
if(numarg()>2) del=$3 else del=1
s=new String2()
sprint(s.s,"/usr/site/pkgs/download/whirlgif304/whirlgif -globmap -o %s %s",$s1,$s2)
system(s.s)
if(del) {system("rm gif/wg/*.gif") printf("whirlgif WARN: deleted %s !\n",$s2)}
return 1
}