-
Notifications
You must be signed in to change notification settings - Fork 0
/
eTrace-p.hoc
4386 lines (3885 loc) · 149 KB
/
eTrace-p.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
/* ================================================================================
Operations on Traces
================================================================================ */
eTrace_loaded_ = 1 //- USE load_file( "eTrace-p.hoc", "eTrace_loaded" ) to avoid multiple loadings
//================================================================================
// Initialize to steady state by iterating in negative time (see NEURON Book).
begintemplate initss
public cvode, v_init, idebug, t0, dt0, t1, dt1, t2, dt2, T, V, G, dV, sd
objref T, V, G, cvode
proc init(){
cvode = $o1
v_init = -65
if( numarg() > 1 ) v_init = $2
idebug = 0
t2 = 1e9 // start at time -t2-t1-t0. finish at t=0
dt2 = 1e8 // iterate each with dt2, dt1 and dt0 respectively
t1 = 1e7 //100
dt1 = 1e3//0.1
t0 = 15//30
dt0 = 0.01//0.1
// dV = maximal variation in t0 window
// sd = SD in t0 window
}
public init_steady_state
proc init_steady_state(){ local dtsav, was_active
if( idebug > 1 ) {
sprint( tstr, "%s", G )
if( 0 == strcmp( tstr, "NULLobject" )) G = new Graph()
printf( "initss t %g dt %g tstop %g cvode %g\n", t, dt, tstop, cvode.active() )
}
it0 = -1 * t0
it1 = -1*( t0 + t1 )
it2 = -1*( t0 + t1 + t2 )
finitialize( v_init )
dtsav = dt // record dt and cvode state
was_active = cvode.active()
t = it2
dt = dt2
cvode.active( 0 )
while( t < it1 ) {
if( idebug > 1 ) {
printf( "s1 %g %g\n", t, v )
G.mark( t, v )
}
fadvance()
}
/*
if( idebug > 1 ) printf( "initss Done s1 t %g v %g\n", t, v )
dt = dt1
cvode.active( 1 )
cvode.atol(1.e-3)
while( t < it0 ) {
if( idebug > 1 ) {
printf( "s2 %g %g\n", t, v )
G.mark( t, v )
}
fadvance()
}
if( idebug > 1 ) printf( "initss Done s2 t %g v %g\n", t, v )
*/
T = new Vector()
V = new Vector()
cvode.active( 0 )
dt = dt0
t = it0
while( t < 0 ){
if( idebug > 1 ) printf( "s3 %g %g\n", t, v )
T.append( t )
V.append( v )
fadvance()
}
if( idebug > 1 ) V.line( G, T )
dV = V.max() - V.min()
sd = V.stdev()/sqrt( dt0 )
if( idebug > 0 ) printf( "initss Done s3 t %g v %g dV %g SD %g\n", t, v, dV, sd )
cvode.active( was_active ) // restore cvode and dt state
dt = dtsav
t = 0
if( cvode.active() ) {
cvode.re_init()
} else {
fcurrent()
}
frecord_init()
}
endtemplate initss
//================================================================================
begintemplate Stat
public id, n, x, avg, sd, cv, min, max
strdef id
proc init(){
if( numarg()>0 ) id = $s1
n = avg = S = cv = 0
min = 1e60
max = -1e60
}
public push
proc push(){ local avgm1
x = $1
n += 1
if( x > max ) max = x
if( x < min ) min = x
if( n==1 ){
avg = x
S = 0
} else {
avgm1 = avg
avg = avgm1 + (x - avgm1)/n
S = S + (x - avgm1) * (x - avg)
}
if( n>1 ) {
sd = sqrt( S/(n-1) )
if( avg != 0.0 ) cv = sd / avg
}
}
endtemplate Stat
//================================================================================
begintemplate str_obj
public s1, s2, s3, s4, s5, sres, aN, aS, idebug
strdef s1, s2, s3, s4, s5, sres, stmp
objref this, aN, aS, SF
proc init(){
if( numarg() > 0 ) s1 = $s1
if( numarg() > 1 ) s2 = $s2
if( numarg() > 2 ) s3 = $s3
if( numarg() > 3 ) s4 = $s4
if( numarg() > 4 ) s5 = $s5
idebug = 0
SF = new StringFunctions()
}
//----------------------------------------------------------------------
public copy
obfunc copy(){ local i localobj sN
sN = new str_obj()
sN.s1 = s1
sN.s2 = s2
sN.s3 = s3
sN.s4 = s4
sN.s5 = s5
return sN
}
// === split( s1 [, seps] ) returns a list of str_objs using seps to split s1
// default seps=[ \t\n\r] is whitespace
public split
obfunc split(){ local i localobj so, sf, sL
so = new str_obj( $s1 )
so.s2 = "[ \t\n\r]" // default seps
if( numarg()>1 ) so.s2 = $s2
sf = new StringFunctions()
sL = new List()
while( sf.head( so.s1, so.s2, so.s3 ) >= 0 ){
if(sf.len(so.s3)>0 ) sL.append( new str_obj( so.s3 ))
sf.tail( so.s1, so.s2, so.s1 )
}
if( sf.len(so.s1) > 0 ) sL.append( new str_obj( so.s1 ))
return sL
}
public nsystem // remove nl from system result
func nsystem(){ local sret localobj sf_
sf_ = new StringFunctions()
sret = system( $s1, $s2 )
if( sf_.substr( $s2, "\n" ) > 0 ) sf_.head( $s2, "\n", $s2 )
return sret
}
public substr // ( "substr" [,i] ) return value of SF.substr( si, "substr" ). i=1 Default.
func substr(){
if( numarg()==1 ) return SF.substr( s1, $s1 )
sprint( stmp, "stmp = s%d", $2 )
execute( stmp, this )
return SF.substr( stmp, $s1 )
}
public system // ( format, sprint-parameters ... )
proc system(){ local i
stmp = ""
aN = new Vector( numarg() ) // for numeric args
aS = new List() // for string args
for i=2, numarg(){
aS.append( new str_obj(""))
if( argtype(i) == 0 ) {
aN.x[i-2] = $i sprint( stmp, "%s, aN.x[%d]", stmp, i-2 )
}
if( argtype(i) == 2 ) {
aS.o(i-2).s1 = $si sprint( stmp, "%s, aS.o(%d).s1", stmp, i-2 )
}
if( idebug ) printf( "system stmp >>%s<<\n", stmp )
}
sprint( stmp, "sprint( stmp, \"%s\" %s )", $s1, stmp )
if( idebug ) printf( "stmp >>%s<<\n", stmp )
execute( stmp, this )
if( idebug ) printf( "stmp >>%s<<\n", stmp )
sprint( stmp, "system( \"%s\", %s.sres )", stmp, this )
if( idebug ) printf( "stmp >>%s<<\n", stmp )
execute( stmp )
// system( stmp )
}
public sPrint // ( [i,] format, sprint-parameters ... )
// Analogous to sprint. An sprint function; places result in $si (Default=1)
obfunc sPrint(){ local i, nda localobj str
str = new str_obj()
if( argtype(1) == 0 ){ // numeric first arg. Use to set si
sprint( str.s1, "s%d", $1 )
str.s2 = $s2 // format
nda = 3
}
if( argtype(1) == 2 ){ // string first arg. Use as format and set to s1
str.s1 = "s1"
str.s2 = $s1 // format
nda = 2
}
aN = new Vector( numarg() ) // for numeric args
aS = new List() // for string args
sprint( str.s3, "sprint( %s, \"%s\"", str.s1, str.s2 )
for i=nda, numarg(){
aS.append( new str_obj( "" ))
if( argtype(i) == 0 ) {
aN.x[i-nda] = $i
sprint( str.s3, "%s, aN.x[%d]", str.s3, i-nda )
}
if( argtype(i) == 2 ) {
aS.o(i-nda).s1 = $si
sprint( str.s3, "%s, aS.o(%d).s1", str.s3, i-nda )
}
}
sprint( str.s3, "%s )", str.s3 )
if( idebug ) printf( "sprint command >>%s<<\n", str.s3 )
execute( str.s3, this )
return this
}
endtemplate str_obj
//================================================================================
begintemplate jQueue
objref q
public add, sum, sum2, pos_count
//--------------------------------------------------------------------------------
proc init(){ //- ( size )
size = $1
q = new Vector( size, 0 )
i = -1
sum = 0
sum2 = 0
pos_count = 0
}
//--------------------------------------------------------------------------------
proc add( ){ //- ( number )
a = $1
i = (i + 1) % size
sum = sum + a - q.x[i]
sum2 += a^2 - q.x[i]^2
if( a > 0 ) pos_count += 1
if( q.x[i] > 0 ) pos_count -= 1
q.x[i] = a
}
endtemplate jQueue
//================================================================================
begintemplate eTrace
objref this
strdef tstr1, tstr2
public col_t, col_v, col_sd, col_n, row_n, name, short_name
public vec_t, vec_v, vec_i, vec_sd, vec_dv, vec_ddv, vec_tmp, vec_ddv_sp
public plot, addplot, plotG, plot_color, plot_brush, plot_recenter, load_file, smooth_exp, insert_ends
public set_dt_sample, dt_sample, load_Vectors, save_file
public stim_del, stim_dur, stim_amp, has_sd, plot_sd, plot_i, has_i
public zc_cache, idebug
strdef header, name, short_name, buf, tok
objref vec_t, vec_v, vec_i, vec_sd, vec_dv, vec_ddv, vec_tmp, vec_ddv_sp, edf, plotG, g2, g3
strdef zc_cache, seg_cache, VdV_cache
public seg_nspikes, eAvg_n, tmp, sf_cutoff
//--------------------------------------------------------------------------------
proc init(){ //- ( )
col_t = 1 //- Column at which to find time in file
col_v = 2 //- Column at which to find voltage in file
col_sd = 3 //- Column at which to find voltage SD in file
col_n = 2 //- Number of columns in file
has_sd = 0 //- files have no SDs by default
plot_color = 1
plot_brush = 1
plot_recenter = 1
plot_sd = 0 //- 0-no plot; 1 plot on zero; 2- plot +- trace
plot_i = 0 //- 0-no plot; 1 plot
vec_t = new Vector() //- Time vector
vec_v = new Vector() //- Voltage vector
vec_i = new Vector() //- current vector (atf files)
has_i = 0 //- 1 if it has current
vec_sd = new Vector() //- SD of Voltage vector
vec_dv = new Vector() //- First derivative
vec_ddv = new Vector() //- Second derivative
vec_tmp = new Vector() //- For various uses
// Defaults. stim_xxx are read from TraceSet files
stim_del = 20 //- stimulation delay
stim_dur = 500 //- stimulation duration
stim_amp = 0.8 //- stimulation amplitude
zc_dt = 0.4 //- zc window size (ms); all points on left >=0 and on righ <= 0
zc_rms_th = 0 //- default rms threshold (factor) to detect zero crossings (zc)
zc_cache = "" //- zc cache flag
sf_cutoff = 200 // cutoff for variable_sf
seg_cache = "" //- spk cache flag
seg_nspikes = -1 //- Num of spikes found by seg_look
VdV_cache = "" //- VdV cache flag
idebug = 0
htf_ntno = 0 // number or traces in htf loaded file
eAvg_n = 0 // to count number of traces for eAvg
tmp = 0 // temp var
}
//--------------------------------------------------------------------------------
public copy
obfunc copy(){ localobj eTc
eTc = new eTrace()
eTc.name = name
eTc.dt_sample = dt_sample
eTc.col_t=col_t
eTc.col_v=col_v
eTc.col_sd=col_sd
eTc.col_n=col_n
eTc.plot_color=plot_color
eTc.plot_brush=plot_brush
eTc.plot_recenter=plot_recenter
eTc.plot_sd=plot_sd
eTc.vec_v=vec_v.c
eTc.vec_t=vec_t.c
eTc.vec_i=vec_i.c eTc.has_i=has_i
eTc.vec_v=vec_v.c
eTc.vec_sd=vec_sd.c
eTc.vec_dv=vec_dv.c
eTc.vec_ddv=vec_ddv.c
eTc.stim_del=stim_del
eTc.stim_dur=stim_dur
eTc.stim_amp=stim_amp
eTc.zc_dt=zc_dt
eTc.zc_rms_th=zc_rms_th
eTc.zc_cache=zc_cache
eTc.seg_cache=seg_cache
eTc.idebug=idebug
return eTc
}
//--------------------------------------------------------------------------------
public systemLF
func systemLF(){ local sret localobj sf_
sf_ = new StringFunctions()
sret = system( $s1, $s2 )
if( sf_.substr( $s2, "\n" ) > 0 ) sf_.head( $s2, "\n", $s2 )
return sret
}
//--------------------------------------------------------------------------------
proc load_file(){ localobj sf //- ( name )
edf = new File()
name = $s1
if( numarg()>= 2 ) short_name = $s2
sf = new StringFunctions()
edf.ropen( name )
comm_n = 0 //- Count comment lines; lines with '#' as first char
go_on = 1
if( idebug ) printf( "load_file: comm_n %d col_n %d col_t %d col_v %d col_sd %d has_sd %d\n", \
comm_n, col_n, col_t, col_v, col_sd, has_sd )
while( go_on ){
edf.gets( header )
if( sf.head( header, "#", tok ) != 0 ) {
go_on = 0
} else {
if( sf.tail( header, "#nrn ", tok ) > -1 ){
sf.left( tok, sf.len(tok)-1 )
sprint( tok, "%s", tok )
if( idebug ) printf( "Executing %s\n", tok )
execute( tok, this )
}
comm_n += 1
}
}
if( idebug ) printf( "load_file: comm_n %d col_n %d col_t %d col_v %d col_sd %d has_sd %d\n", \
comm_n, col_n, col_t, col_v, col_sd, has_sd )
edf.ropen( name ) //- reopen file
for i=1, comm_n edf.gets( header ) //- read comment lines (ignored)
vec_t.scanf( edf, col_t, col_n )
// vec_t.mul( 1000 ) //- Convert to milliseconds
edf.ropen( name ) //- reopen file
for i=1, comm_n edf.gets( header ) //- read comment lines (ignored)
vec_v.scanf( edf, col_v, col_n )
if( has_sd ){
edf.ropen( name ) //- reopen file
for i=1, comm_n edf.gets( header ) //- read comment lines (ignored)
vec_sd.scanf( edf, col_sd, col_n )
}
set_dt_sample()
clear_cache()
edf.close()
}
//--------------------------------------------------------------------------------
public load_atf //( filename, V_colN [, I_colN ] )
obfunc load_atf(){ local colV, colI, i, ivS, bS localobj edf, str, tB
has_i = 0
name = $s1
colV = $2 // tno for voltage readings
colI = -1
if( numarg() > 2 ) colI = $3
str = new str_obj()
edf = new File()
edf.ropen( name )
if( ! edf.isopen()){
printf( "\n\n \t\tload_atf ERROR: >>%s<< not found\n\n", name )
return
}
edf.scanstr( str.s1 )
if( strcmp( "ATF", str.s1 ) != 0 ) {
printf( "\n\n \t\tload_atf ERROR: Not ATF format.\n\n" )
return
}
edf.scanstr( str.s2 )
ncoms = edf.scanvar()
ncols = edf.scanvar()
if ( idebug > 1 ) printf( "ncoms %g ncols %g\n", ncoms, ncols )
if( colV > ncols || colI > ncols ){
printf( "\n\n \t\tload_atf ERROR: trial colV=%g or colI=%g not found in >>%s<<\n\n", colV, colI, name )
return
}
vec_t = new Vector()
vec_v = new Vector()
if( colI > 0 ) vec_i = new Vector()
bS = 0
tB = new Vector( ncols )
while( ! edf.eof ){
if( bS==0 || ( bS <= vec_t.size ) ){
if( bS==0 ) bS = 2^14
if( idebug > 1 ) printf( "bS %g vec_t.size %g\n", bS, vec_t.size )
bS = 2 * bS
vec_t.buffer_size( bS )
vec_v.buffer_size( bS )
if( colI > 0 ) vec_i.buffer_size( bS )
}
tB.scanf( edf, ncols )
vec_t.append( tB.x[0]*1000 )
vec_v.append( tB.x[colV-1] )
if( colI > 0 ) vec_i.append( tB.x[colI-1] )
}
set_dt_sample()
clear_cache()
edf.close()
return this
}
//--------------------------------------------------------------------------------
// Use HOME environment variable to resolve file name if not begin with "/"
public load_atf_home
obfunc load_atf_home(){ localobj st, sf
sf = new StringFunctions()
st = new str_obj( $s1 )
if( sf.substr( $s1, "/" ) != 0 ){
systemLF( "echo $HOME", st.s1 )
sprint( st.s1, "%s/%s", st.s1, $s1 )
}
if( numarg() >= 3 ) load_atf( st.s1, $2, $3 )
if( numarg() < 3 ) load_atf( st.s1, $2 )
return this
}
//--------------------------------------------------------------------------------
// Use HOME environment variable to resolve file name if not begin with "/"
// take colNs as text file:c1:c2
public load_atf_home_S
proc load_atf_home_S(){ local c1, c2 localobj str, sf
sf = new StringFunctions()
str = new str_obj( $s1 )
sf.tail( str.s1, ":", str.s2 )
sf.head( str.s1, ":", str.s1 )
c2 = -1
if( sscanf( str.s2, "%d:%d", &c1, &c2 ) < 1 ){
printf( "\n\n\tload_atf_home_S: ERROR wrong parameters. expecting fname:x:y got >>%s<<\n", $s1 )
return
}
load_atf_home( str.s1, c1, c2 )
}
//--------------------------------------------------------------------------------
public save_atf
proc save_atf(){ local i, ym, yM, im, iM localobj str, eF, hiS // ( eT, fname [,comm] )
str = new str_obj( $s1 )
if( numarg() > 1 ) str.s2=$s2
ym = im = 1e60
yM = iM = -1e60
for i=0, vec_v.size-1 {
if( vec_v.x[i] > yM ) yM = vec_v.x[i]
if( vec_v.x[i] < ym ) ym = vec_v.x[i]
if( has_i ){
if( vec_i.x[i] > iM ) iM = vec_i.x[i]
if( vec_i.x[i] < im ) im = vec_i.x[i]
}
}
hiS = new str_obj()
if( has_i ){
sprint( hiS.s1, ",%g", iM )
sprint( hiS.s2, ",%g", im )
sprint( hiS.s3, "\t\"Im\"" )
sprint( hiS.s4, "\t\"Trace #1 (nA)\"" )
sprint( hiS.s5, ",Im" )
}
eF = new File()
eF.wopen( str.s1 )
eF.printf( "ATF\t1.0\n7\t%d\n", 2+has_i )
eF.printf( "\"AcquisitionMode=Episodic Stimulation\"\n" )
eF.printf( "\"Comment=%s Created by save_atf jose@kiubo.net c2006\"\n", str.s2 )
eF.printf( "\"YTop=%g%s\"\n", yM, hiS.s1 )
eF.printf( "\"YBottom=%g%s\"\n", ym, hiS.s2 )
eF.printf( "\"SweepStartTimesMS=0.000\"\n" )
eF.printf( "\"SignalsExported=Vm%s\"\n", hiS.s5 )
eF.printf( "\"Signals=\"\t\"Vm\"%s\n", hiS.s3 )
eF.printf( "\"Time (ms)\"\t\"Trace #1 (mV)\"%s\n", hiS.s4 )
for i=0, vec_v.size-1 {
eF.printf( "%g\t%g", vec_t.x[i], vec_v.x[i] )
if( has_i ) eF.printf( "\t%g", vec_i.x[i] )
eF.printf( "\n" )
}
eF.close()
}
//--------------------------------------------------------------------------------
public load_atf_all
obfunc load_atf_all(){ local i, bS localobj edf, str, vL, tB
str = new str_obj( $s1 )
vL = new List()
edf = new File()
edf.ropen( str.s1 )
if( ! edf.isopen()){
printf( "\n\n \t\tload_atf ERROR: >>%s<< not found\n\n", str.s1 )
return
}
edf.scanstr( str.s2 )
if( strcmp( "ATF", str.s2 ) != 0 ) {
printf( "\n\n \t\tload_atf ERROR: Not ATF format.\n\n" )
return
}
edf.scanstr( str.s2 )
ncoms = edf.scanvar()
ncols = edf.scanvar()
for i=0, ncols-1 vL.append( new Vector() )
bS = 0
tB = new Vector( ncols )
while( ! edf.eof ){
if( bS==0 || ( bS <= vec_t.size ) ){
if( bS==0 ) bS = 2^14
if( idebug > 1 ) printf( "bS %g vec_t.size %g\n", bS, vec_t.size )
bS = 2 * bS
for i=0, ncols-1 vL.o(i).buffer_size( bS )
}
tB.scanf( edf, ncols )
for i=0, ncols-1 vL.o(i).append( tB.x[i] )
}
edf.close()
return vL
}
//--------------------------------------------------------------------------------
public save_atf_all_htf
proc save_atf_all_htf(){ local i localobj vL, str, eF, aV
vL = $o1
str = new str_obj( $s2 )
eF = new File()
eF.wopen( str.s1 )
if( ! eF.isopen()){
printf( "\n\n \t\tsave_atf_all_htf ERROR cant write to >>%s<<\n\n", str.s1 )
return
}
aV = new Vector()
aV.buffer_size( vL.count * vL.o(0).size + 100 )
aV.append( vL.o(0).size )
aV.append( vL.count )
for i=0, vL.count-1 aV.append( vL.o(i) )
aV.vwrite( eF, 4 )
eF.close()
}
//--------------------------------------------------------------------------------
public save_htf_1p0
proc save_htf_1p0(){ local i, nc localobj eL, str, sf, eF, aV
eL = $o1 // list of eTraces (all same time index)
str = new str_obj( $s2 ) // filename in str.s1
if( numarg()>2 ) str.s2 = $s3 // optional comment
for i=0, eL.count-1 {
if( eL.o(0).vec_v.size != eL.o(i).vec_v.size ) \
printf( "\n\n\tsave_htf_1p0: WARNING diff size between 0 and %d\n\n", i )
}
if( vec_v.size != vec_t.size ){
printf( "\n\n \t\tsave_htf_1p0 ERROR: vec_t.size %d != vec_v.size %d\n", vec_t.size, vec_v.size )
return
}
sf = new StringFunctions()
nc = 0
str.s3 = str.s2
while( sf.tail( str.s3, "\n", str.s3 )!= -1 ) nc += 1 // count comments
if( idebug ) printf( "save_htf_1p0: nc %d\n", nc )
eF = new File()
eF.wopen( str.s1 )
if( ! eF.isopen()){
printf( "\n\n \t\tsave_htf_1p0 ERROR cant write to >>%s<<\n\n", str.s1 )
return
}
eF.printf( "iHTF1.0 %d\n%s\n", nc+1, str.s2 )
aV = new Vector()
aV.buffer_size( (eL.count + 1) * eL.o(0).vec_t.size + 100 ) // Add one for time index vector
aV.append( eL.o(0).vec_t.size ) // x[0] size
aV.append( eL.count ) // x[1] number of trials
aV.append( eL.o(0).vec_t ) // time index
for i=0, eL.count-1 aV.append( eL.o(i).vec_v ) // the vec_v parts
aV.vwrite( eF, 4 )
eF.close()
}
//--------------------------------------------------------------------------------
public load_htf_1p0, htf_ntno
//( filename, tno [, Root-env-var] )
obfunc load_htf_1p0(){ local tno, ncom, i, nsize localobj eF, str, aV
has_i = 0
str = new str_obj( $s1 )
tno = $2 // trial number to load [1..htf_ntno] (time vector in zero)
if( numarg()>2 ){
sprint( str.s4, "echo $%s", $s3 ) // prepend EnvVar to filename
systemLF( str.s4, str.s5 )
//sprint( str.s1, "%s/%s", str.s5, str.s1 )
if( idebug ) printf( "load_htf_1p0: opening >>%s<<\n", str.s1 )
}
sprint( name, "%s:%d", str.s1, tno )
eF = new File()
aV = new Vector()
eF.ropen( str.s1 )
if( ! eF.isopen()){
printf( "\n\n \t\tload_htf_1p0 ERROR: >>%s<< not found\n\n", str.s1 )
return
}
eF.gets( str.s4 )
if( 1 != sscanf( str.s4, "iHTF1.0 %d", &ncom )){
printf( "\n\n \t\tload_htf_1p0 ERROR: expecting iHTF1.0 got >>%s<< \n\n", str.s4 )
return
}
for i=1, ncom eF.gets( str.s4 ) // skip comments
aV.vread( eF )
nsize = aV.x[0]
htf_ntno = aV.x[1]
if( tno > htf_ntno ){
printf( "\n\n \t\tload_atf ERROR: trial tno=%g not found in >>%s<<\n\n", tno, str.s1 )
return
}
if( idebug ) printf( "load_htf_1p0: aV.size %d nsize %d tno %d\n", aV.size, nsize, tno )
vec_t = aV.c( 2, 2+nsize-1 )
vec_v = aV.c( 2 + (tno)*nsize, 2 + (tno)*nsize + nsize -1 )
set_dt_sample()
clear_cache()
eF.close()
return this
}
public load_htf_1p0c // ( tno:filename [, Env-var] )
obfunc load_htf_1p0c(){ local tno localobj str
str = new str_obj( $s1 )
sscanf( str.s1, "%d:%s", &tno, str.s2 )
if( numarg()>1 ) { return load_htf_1p0( str.s2, tno, $s2 )
}else{ return load_htf_1p0( str.s2, tno ) }
}
//--------------------------------------------------------------------------------
public atf_2_htf
proc atf_2_htf(){ local i localobj vL
vL = load_atf_all( $s1 )
save_atf_all_htf( vL, $s2 )
}
//--------------------------------------------------------------------------------
public load_htf //( filename, V_colN [, I_colN ] )
obfunc load_htf(){ local colV, colI, i, nsize, ncols localobj edf, str, aV
has_i = 0
name = $s1
colV = $2 // tno for voltage readings
colI = -1
if( numarg() > 2 ) colI = $3
str = new str_obj()
edf = new File()
aV = new Vector()
edf.ropen( name )
if( ! edf.isopen()){
printf( "\n\n \t\tload_atf ERROR: >>%s<< not found\n\n", name )
return
}
aV.vread( edf )
nsize = aV.x[0]
ncols = aV.x[1]
if( colV > ncols || colI > ncols ){
printf( "\n\n \t\tload_atf ERROR: trial colV=%g or colI=%g not found in >>%s<<\n\n", colV, colI, name )
return
}
vec_t = aV.c( 2, 2+nsize-1 ).mul( 1000 )
vec_v = aV.c( 2 + (colV-1)*nsize, 2 + (colV-1)*nsize + nsize -1 )
if( colI > 0 ) vec_i = aV.c( 2 + (colI-1)*nsize, 2 + (colI-1)*nsize + nsize -1 )
set_dt_sample()
clear_cache()
edf.close()
return this
}
//--------------------------------------------------------------------------------
// Use HOME environment variable to resolve file name if not begin with "/"
public load_htf_home
obfunc load_htf_home(){ localobj st, sf
sf = new StringFunctions()
st = new str_obj( $s1 )
/*
if( sf.substr( $s1, "/" ) != 0 ){
systemLF( "echo $HOME", st.s1 )
sprint( st.s1, "%s/%s", st.s1, $s1 )
}
if( numarg() >= 3 ) load_htf( st.s1, $2, $3 )
if( numarg() < 3 ) load_htf( st.s1, $2 )
*/
load_htf($s1)
return this
}
//--------------------------------------------------------------------------------
// Use HOME environment variable to resolve file name if not begin with "/"
// take colNs as text file:c1:c2
public load_htf_home_S
proc load_htf_home_S(){ local c1, c2 localobj str, sf
sf = new StringFunctions()
str = new str_obj( $s1 )
sf.tail( str.s1, ":", str.s2 )
sf.head( str.s1, ":", str.s1 )
c2 = -1
if( sscanf( str.s2, "%d:%d", &c1, &c2 ) < 1 ){
printf( "\n\n\tload_htf_home_S: ERROR wrong parameters. expecting fname:x:y got >>%s<<\n", $s1 )
return
}
load_htf_home( str.s1, c1, c2 )
}
//--------------------------------------------------------------------------------
proc load_Vectors(){ //- load_Vectors( vec_t, vec_v )
vec_v = $o1.c
vec_t = $o2.c
set_dt_sample()
clear_cache()
}
//--------------------------------------------------------------------------------
// Use cNeuro environment variable to make code more portable
public load_cNeuro
proc load_cNeuro(){
systemLF( "echo $cNeuro", tstr1 )
sprint( tstr1, "%s/%s", tstr1, $s1 )
if( numarg() >= 2 ) load_file( tstr1, $s2 )
if( numarg() < 2 ) load_file( tstr1 )
}
//--------------------------------------------------------------------------------
// Use HOME environment variable to resolve file name if not begin with "/"
public load_file_home
proc load_file_home(){ localobj st, sf
sf = new StringFunctions()
st = new str_obj( $s1 )
if( sf.substr( $s1, "/" ) != 0 ){
systemLF( "echo $HOME", st.s1 )
sprint( st.s1, "%s/%s", st.s1, $s1 )
}
if( numarg() >= 2 ) load_file( st.s1, $s2 )
if( numarg() < 2 ) load_file( st.s1 )
}
//--------------------------------------------------------------------------------
public load_wcp
obfunc load_wcp(){ local tno, chn, ljp localobj str
str = new str_obj( $s1 )
tno = $2
chn = $3
ljp = $4
sprint( str.s2, "Emes5 %s %d %d -xWcp %d -dcOff -show_Trace > tmp_wcp.txt", str.s1, tno, chn, ljp )
if( idebug ) printf( "load_wcp: executing >>%s<<\n", str.s2 )
system( str.s2 )
load_file( "tmp_wcp.txt" )
sprint( name, "%s-t=%d-ch=%d-ljp=%d", str.s1, tno, chn, ljp )
return this
}
//--------------------------------------------------------------------------------
public load_wcp_all
obfunc load_wcp_all(){ local i, nrec, nchan, bS localobj edf, str, vL, tB
str = new str_obj( $s1 )
vL = new List()
edf = new File()
edf.ropen( str.s1 )
if( ! edf.isopen()){
printf( "\n\n \t\tload_wcp_all ERROR: >>%s<< not found\n\n", str.s1 )
return vL
}
sprint( str.s2, "egrep -a 'NR=|NC=' %s | sort | tr '\n\r' ' '", str.s1 )
system( str.s2, str.s3 )
sscanf( str.s3, "NC=%d NR=%d", &nchan, &nrec )
if( idebug ) printf( "load_wcp_all: NC %g NR %g\n", nchan, nrec )
edf.close()
// for ic=0, nchan vL.append( new List() )
for ic=0, nchan-1 {
for ir=1, nrec {
load_wcp( str.s1, ir, ic, 0 )
if( ir==1 && ic==0 ) vL.append( vec_t.c ) // place time as first member
vL.append( vec_v.c )
}
}
return vL
}
//--------------------------------------------------------------------------------
proc save_file(){ localobj edf //- ( name )
edf = new File()
name = $s1
edf.wopen( name )
edf.printf( "#nrn col_n=2 col_t=1 col_v=2 stim_del=%g stim_dur=%g stim_amp=%g\n", \
stim_del, stim_dur, stim_amp )
if( numarg() > 1 ) edf.printf( $s2 )
for i=0, vec_t.size()-1 edf.printf( "%.10g %.10g\n", vec_t.x[i], vec_v.x[i] )
edf.close()
}
//--------------------------------------------------------------------------------
public parse_name, cell_id, ms_id, pA_id, cond_id, num_id
strdef cell_id, ms_id, pA_id, cond_id, num_id
proc parse_name(){ localobj sf
sf = new StringFunctions()
sf.tail( name, "blk/|ACSF/", ms_id )
sf.head( ms_id, "/", cell_id )
sf.tail( ms_id, "/", ms_id )
sf.head( ms_id, ".txt", ms_id )
sf.tail( ms_id, "_", cond_id )
sf.head( ms_id, "_", ms_id )
sf.tail( ms_id, "ms", pA_id )
sf.head( ms_id, "\\+|-", ms_id )
sf.tail( cond_id, "_", num_id )
sf.head( cond_id, "_", cond_id )
}
//--------------------------------------------------------------------------------
public set_new_dt
obfunc set_new_dt(){ local newdt localobj taux
newdt = $1
taux = new Vector()
taux.indgen( vec_t.x[0], vec_t.x[vec_t.size-1], newdt )
vec_v.interpolate( taux, vec_t )
vec_t = taux.c
dt_sample = newdt //taux.x[1] - taux.x[0]
return this
}
//--------------------------------------------------------------------------------
// dt_sample = -1 if sampling freq is not uniform; otherwise dt_sample = sample interval
//
proc set_dt_sample(){ local i
dt_sample = vec_t.x[1]-vec_t.x[0]
for i=1, vec_t.size()-1 {
if( vec_t.x[i]-vec_t.x[i-1] == dt_sample ) continue
dt_sample = -1
break
}
}
//--------------------------------------------------------------------------------
obfunc plot() {
// if( object_id( plotG ) == 0 ) {
plotG = new Graph()
plotG.family( 1 ) //- To keep lines
if( numarg() >= 1 ) plot_color = $1
if( numarg() >= 2 ) plot_brush = $2
plotG.color( plot_color )
plotG.brush( plot_brush )
vec_v.line( plotG, vec_t )
plotG.label( name )
if( has_sd && plot_sd>0 ) {
plot_color += 1
plotG.color( plot_color )
if( plot_sd == 1 ) vec_sd.line( plotG, vec_t )
if( plot_sd==2 ) vec_sd.c.add(vec_v).line( plotG, vec_t )
}
if( has_i && plot_i ) vec_i.line( plotG, vec_t )
recenter()
return this
}
//--------------------------------------------------------------------------------
proc addplot() {
if( numarg() > 1 ) {
plot_color += 1
if( plot_color%10 == 0 || plot_color%10 == 8 ) plot_color += 1
}
if( object_id( plotG )==0 ) { plotG = new Graph() }
plotG.color( plot_color )
plotG.brush( plot_brush )
$o1.vec_v.line( plotG, $o1.vec_t )
plotG.label( $o1.name )
if( $o1.has_i && plot_i ) $o1.vec_i.line( plotG, $o1.vec_t )
recenter()
}
//--------------------------------------------------------------------------------
public addplot_c
proc addplot_c() { local col localobj eT
eT = $o1
if( numarg() > 1 ) {
plot_color = $2
}
if( numarg() > 2 ) {
plot_brush = $3
}
if( object_id( plotG )==0 ) { plotG = new Graph() }
plotG.color( plot_color )
plotG.brush( plot_brush )
eT.vec_v.line( plotG, eT.vec_t )
plotG.label( eT.name )
if( eT.has_i && plot_i ) eT.vec_i.line( plotG, eT.vec_t )
recenter()
}
//--------------------------------------------------------------------------------
public recenter
proc recenter() {
if( plot_recenter ) plotG.exec_menu( "View = plot" )
}
//--------------------------------------------------------------------------------
// plot v,dv,ddv according to option XXX (X=1 plot, X=0 no plot)
public plot_ddv
obfunc plot_ddv(){ local i, opt, col, vm, dvm, ddvm localobj pG
pG = $o1
opt = $2
{ col=1 if( numarg()>2 ) col = $3 }
{ vm=0.1 dvm=0.01 ddvm=0.001 if( numarg()>3 ){ vm=$4 dvm= $5 ddvm=$6 }}
if( object_id( pG )==0 ) { pG = new Graph() }
pG.color( col )
if( int(opt/100)%10 ) vec_v.c.mul(vm).line( pG, vec_t )
if( int(opt/10)%10 ) vec_dv.c.mul(dvm).line( pG, vec_t )
if( int(opt/1)%10 ) vec_ddv.c.mul(ddvm).line( pG, vec_t )
pG.label( name )
pG.exec_menu( "View = plot" )
return pG
}