-
Notifications
You must be signed in to change notification settings - Fork 0
/
updown.hoc
167 lines (146 loc) · 5.73 KB
/
updown.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
// $Id: updown.hoc,v 1.121 2013/02/15 22:42:34 samn Exp $
print "Loading updown.hoc..."
declare("printStep",0.1)
CREEP_UPDOWN=1
NOV_UPDOWN=1
declare("test_do_graphs", 1)
declare("drawlr", 1, "drawth", 1,"drawscale",0,"shapesz",15,"flipit",0)
declare("minthreshlx", 0, "maxthreshlx",9e3*printStep)
declare("black",1,"red",2,"blue",3,"green",4)
declare("vthg",new Vector()) // global, stores threshold lines
// EXAMPLE USAGE
// {gvnew(-2) panobj.read_vfile("/u/billl/nrniv/sync/data/05nov02.501/v05nov02.1000")}
// panobj.rv(7)
// vec.copy(panobj.vrtmp)
// vec.mul(-1)
// gg(vec,printStep)
// objref output
// output=updown(vec)
// output.pr(10)
// updown(vec[,nq,#CUTS,LOGCUT,MIN]) -- use updown() to find spikes
// LOC(0) PEAK(1) WIDTH(2) BASE(3) HEIGHT(4) START(5) SLICES(6) SHARP(7) INDEX(8) FILE(9)
// other options
declare("pos_updown",1) // set to 1 to move whole curve up above 0
declare("maxp_updown",0.95) // draw top sample at 95% of max
declare("minp_updown",0.05) // draw bottom sample at 5% of max
declare("over_updown",1) // turn over and try again if nothing found
declare("allover_updown",0) // turn over and add these locs (not debugged)
declare("verbose_updown",0) // give messages, can also turn on DEBUG_VECST for messages from updn()
declare("dynamic_th",0)
declare("padsz_updown",5000) // default padsize for updown function
declare("bumpcap_updown",2e4) // default capacity for number of bumps
//** updown - calls Vector updown function and returns NQS with spikes/bumps
obfunc updown () { local a,ii,npts,logflag,min,x,sz,midp localobj bq,cq,v1,v2,v3,bb,tl,vtmp
if (verbose_updown) printf("MAXTIME appears to be %g (printStep=%g)\n",$o1.size*printStep,printStep)
npts=10 // n sample locations
tl=new List()
bq=new NQS("LOC","PEAK","WIDTH","BASE","HEIGHT","START","SLICES","SHARP","INDEX","FILE","NESTED")
if (numarg()>1) cq=$o2 else cq=new NQS()
if (cq.m!=10) {
cq.verbose = 0
cq.resize(0)
cq.resize("LOC","PEAK","WIDTH","BASE","HEIGHT","START","SLICES","SHARP","INDEX","FILE","NESTED") }
cq.verbose = 1
if (numarg()>2) npts=$3
if (numarg()>3) logflag=$4 else logflag=0
if (numarg()>4) {
if (npts!=$o5.size) printf("Correcting npts from %d to %d\n",npts,npts=$o5.size)
if ($o5.ismono(1)) $o5.reverse
if (! $o5.ismono(-1)) {printf("updown: Arg 5 (%s) must be monotonic\n",$o5) return}
}
bq.listvecs(bb)
bq.pad(padsz_updown)
a=allocvecs(npts+3,bumpcap_updown)
v1=mso[a+0] v2=mso[a+1] v3=mso[a+2]
for ii=0,npts-1 tl.append(mso[ii+a+3])
v1.copy($o1)
if (pos_updown) {
min=v1.min
v1.sub(min) // make it uniformly positive
} else min=0
if (numarg()>4) v2.copy($o5) else {
if(logflag==2){ //"double-log" spacing
midp = v1.max * 0.5 //50% of max is center of log axis in vertical direction
// 1/2 of lines btwn max and midpoint
vtmp=new Vector()
vtmp.indgen(2,2+npts/2-1,1)
vtmp.log()
vtmp.scale(maxp_updown*v1.max,midp)
v2.append(vtmp)
// 1/2 of lines btwn min and midpoint
vtmp.indgen(2,2+npts/2-1,1)
vtmp.log()
vtmp.scale(minp_updown*v1.max,midp)
v2.append(vtmp)
//make sure they're monotonically increasing
v2.sort()
v2.reverse()
} else {
v2.indgen(2,2+npts-1,1) // sampling at npts points, start at 2 to avoid log(1)=0
if (logflag) v2.log() // log sampling
v2.scale(-maxp_updown*v1.max,-minp_updown*v1.max) v2.mul(-1)
}
}
if(dynamic_th) { //allocate memory so updown can use th
v2.resize(v1.max()+1)
v2.resize(0)
printf("dynamic_th\n")
}
vthg.copy(v2) // <--------- whats this for? copies threshold lines to global vthg
//v2 = thresh
v1.updown(v2,tl,bb)
if (pos_updown) { bq.v[1].add(min) bq.v[3].add(min) }
cq.append(bq)
sz=bq.size(1)
if (allover_updown) { // do it upside down as well
v1.mul(-1) // v2 will be upside-down
if (pos_updown) {min=v1.min v1.sub(min)}
v1.updown(v2,tl,bb)
bq.v[8].add(sz) bq.v[4].mul(-1) // turn HEIGHT upside down
cq.append(bq)
} else if (over_updown && sz==0) { // turn it over an try again
print "updown() checking upside-down"
v1.mul(-1) // v2 will be upside-down
v1.updown(v2,tl,bb)
}
for case(&x,0,2,5) cq.v[x].mul(printStep)
nqsdel(bq)
dealloc(a)
return cq
}
//** drupdown(voltage,vector with updown output,vth - vector of threshold slices) - draw output from updown
proc drupdown () { local idx,x,left,right localobj vec,output,vth, vslice
vec=$o1 output=$o2 vth=$o3
vslice = new Vector()
if(g==nil) gg()
gg(vec,printStep,black,3) //5 is for thick lines
// output.v[1].mark(g,output.v[0],"o",shapesz,red,1) //draw circles at top of peaks // Sam's line
output.v[1].add(vec.min()).mark(g,output.v[0],"o",shapesz,red,1) //draw circles at top of peaks
if(drawlr) for idx=0,output.v.size-1 {
left = output.v[5].x(idx)
right = left+output.v[2].x(idx)
g.mark(left, vec.x(left/printStep), "t",shapesz,blue ,1)
g.mark(right,vec.x(right/printStep),"s",shapesz,green,1)
}
// draw horizontal lines for slices
// if(drawth) for vtr(&x,vth) drline(minthreshlx,x,maxthreshlx,x,red,1) // Sam's line
vslice = vth.add(vec.min())
if(drawth) for vtr(&x,vth) drline(minthreshlx,x,maxthreshlx,x,red,1)
}
//** testupdown($o1=input vector,$2=num slices,$3=logarithmic-spaced slices,$o4=thresholds slices locations)
// runs updown on input vector and returns an NQS with detected spikes/bumps
obfunc testupdown (){ local idx,left,right,x,a localobj vtmp,output,vec
{a=allocvecs(vec) output = new NQS() vec.copy($o1)}
if(flipit) vec.mul(-1)
vec.sub(vec.min)
if(numarg()>3){
updown(vec,output,$2,$3,$o4)
} else if(numarg()>2){
updown(vec,output,$2,$3)
} else{
updown(vec,output)
}
if(test_do_graphs) drupdown(vec,output,vthg)
dealloc(a)
return output
}