source: project/release/4/nemo/trunk/extensions/nemo-vclamp.scm @ 27021

Last change on this file since 27021 was 27021, checked in by Ivan Raikov, 9 years ago

nemo: bug fixes in hh and vclamp extension modules

File size: 14.1 KB
Line 
1;;       
2;;
3;; An extension for generating NEURON HOC voltage clamp scripts from NEMO models.
4;;
5;; Copyright 2008-2011 Ivan Raikov and the Okinawa Institute of Science and Technology
6;;
7;; This program is free software: you can redistribute it and/or
8;; modify it under the terms of the GNU General Public License as
9;; published by the Free Software Foundation, either version 3 of the
10;; License, or (at your option) any later version.
11;;
12;; This program is distributed in the hope that it will be useful, but
13;; WITHOUT ANY WARRANTY; without even the implied warranty of
14;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15;; General Public License for more details.
16;;
17;; A full copy of the GPL license can be found at
18;; <http://www.gnu.org/licenses/>.
19;;
20
21(module nemo-vclamp
22
23        (nemo:vclamp-translator)
24
25        (import scheme chicken utils data-structures lolevel srfi-1 srfi-13 )
26
27        (require-extension lolevel datatype matchable strictly-pretty 
28                           environments varsubst datatype 
29                           nemo-core nemo-utils nemo-gate-complex)
30
31(define-syntax pp
32  (syntax-rules ()
33    ((pp indent val ...) (ppf indent (quasiquote val) ...))))
34
35(define (hoc-cname sys v)
36  (s+ (hoc-name v) "_" (hoc-name sys)))
37
38(define (hoc-name s)
39  (let ((cs (string->list (->string s))))
40    (let loop ((lst (list)) (cs cs))
41      (if (null? cs) (string->symbol (list->string (reverse lst)))
42          (let* ((c (car cs))
43                 (c1 (cond ((or (char-alphabetic? c) (char-numeric? c) (char=? c #\_)) c)
44                           (else #\_))))
45            (loop (cons c1 lst) (cdr cs)))))))
46                           
47(define (matlab-name s)
48  (let ((cs (string->list (->string s))))
49    (let loop ((lst (list)) (cs cs))
50      (if (null? cs) (string->symbol (list->string (reverse lst)))
51          (let* ((c (car cs))
52                 (c1 (cond ((or (char-alphabetic? c) (char-numeric? c) (char=? c #\_)) c)
53                           (else #\_))))
54            (loop (cons c1 lst) (cdr cs)))))))
55           
56(define (const-val v)
57  (and (nemo:quantity? v)
58       (cases nemo:quantity v
59              (CONST (name value) value)
60              (else (error 'const-val "invalid constant" v)))))
61
62(define (print-hoc-prelude indent)
63
64(pp indent (#<<EOF
65
66vgraphbox.intercept(0)
67
68vec_sizes = tstop/dt + 1        // recorded traces are all this size
69
70soma vce=new VClamp(0)
71
72dt = 0.001
73
74proc vcrun() {
75       
76        vce.dur[0]=vchdur
77        vce.dur[1]=vcbdur
78        vce.dur[2]=vchdur
79
80        tstop=vce.dur[0]+vce.dur[1]+vce.dur[2]
81       
82        vce.amp[0]=vchold
83        vce.amp[2]=vchold
84       
85        for j=0, vcsteps-1 {
86                x=vcbase+j*vcincrement
87                vce.amp[1]=x
88
89                init()
90
91                run()
92                $o2.line(g[$1], dt)
93
94                for i=0,$o2.size()-1 {
95                    $o4.printf("%g %g\n", $o3.x[i], $o2.x[i])
96                    }
97
98                if (stoppedrun()) {
99                        break
100                }
101        }
102}
103EOF
104))
105)
106
107(define (nemo:vclamp-translator sys . rest)
108
109  (define (cid x)  (second x))
110  (define (cn x)   (first x))
111 
112  (let-optionals rest ((target 'hoc) (filename #f))
113
114    (cases nemo:quantity (environment-ref sys (nemo-intern 'dispatch)) 
115
116         (DISPATCH  (dis) 
117
118          (let ((imports  ((dis 'imports)  sys))
119                (exports  ((dis 'exports)  sys)))
120           
121            (let* ((indent        0)
122                   (indent+       (+ 2 indent ))
123                   (eval-const    (dis 'eval-const))
124                   
125                   (sysname       (hoc-name ((dis 'sysname) sys)))
126                   (filename      (or filename (s+ sysname "_vclamp" (case target ((matlab octave) ".m") (else ".ses")))))
127
128                   (consts        ((dis 'consts)  sys))
129                   (asgns         ((dis 'asgns)   sys))
130                   (states        ((dis 'states)  sys))
131                   (reactions     ((dis 'reactions) sys))
132                   (rates         ((dis 'rates) sys))
133                   (defuns        ((dis 'defuns)  sys))
134                   (components    ((dis 'components) sys))
135                   
136                   (vc-components (filter-map (lambda (x) (and (eq? 'voltage-clamp (second x)) x)) components))
137                   (g             (match-let (((state-list asgn-list g) ((dis 'depgraph*) sys))) g))
138                   (poset         (vector->list ((dis 'depgraph->bfs-dist-poset) g)))
139                   
140                   (gate-complex-info    (nemo:gate-complex-query sys))
141                   (gate-complexes       (lookup-def 'gate-complexes gate-complex-info))
142                   
143                   (i-gates       (lookup-def 'i-gates gate-complex-info))
144                   
145                   )
146             
147              (if (not (null? vc-components))
148
149                 
150                  (case target 
151                   
152                    ((hoc)
153
154                     (with-output-to-file filename
155                       (lambda ()
156                         (let* ((vcn   (length vc-components))
157                                (hboxn (inexact->exact (ceiling (/ vcn 3)))))
158                       
159                           (pp indent
160                               ("load_file(\"nrngui.hoc\")")
161                               ("objref vce     // voltage clamp)")
162                               ("objref g[" ,vcn "] // graph objects"))
163
164                           (pp indent
165                               ("objref vgraphbox, hgraphbox[" ,hboxn "]")
166                               ("vgraphbox=new VBox()")
167                               ("vgraphbox.intercept(1)"))
168                           
169                           (let recur ((a 0) (b (min (- vcn 1) 2)) (hbi 0))
170                             (pp indent
171                                 ("hgraphbox[" ,hbi "]=new HBox()")
172                                 ("hgraphbox[" ,hbi "].intercept(1)"))
173                             
174                             (if (< a b)
175                                 (pp indent ("for i=" ,a "," ,b " {"))
176                                 (pp indent ("i = " ,a)))
177                             (pp (if (< a b) indent+ indent)
178                                 ("g[i]=new Graph()")
179                                 ("g[i].exec_menu(\"Keep Lines\")"))
180                             (if (< a b) (pp indent ("}")))
181                             (pp indent
182                                 ("hgraphbox[" ,hbi "].intercept(0)")
183                                 ("hgraphbox[" ,hbi "].map()"))
184                             (if (> hboxn (+ 1 hbi))
185                                 (recur (+ b 1) (min (- vcn 1) (+ b 3)) (+ 1 hbi)))
186                             )
187
188                           (print-hoc-prelude indent)
189
190                           (for-each
191                            (lambda (comp i) 
192                             
193                              (let ((name (first comp))
194                                    (sym  (third comp)))
195                               
196                                (match-let (((vchold vcbase vcinc vcsteps vchdur vcbdur) ((dis 'component-exports) sys sym)))
197
198                                    (pp indent ("print \"generating " ,name "\""))
199
200                                    (pp indent
201                                        ("vchold"      = ,(s+ "soma." (hoc-cname sysname vchold)))
202                                        ("vcbase"      = ,(s+ "soma." (hoc-cname sysname vcbase)))
203                                        ("vcincrement" = ,(s+ "soma." (hoc-cname sysname vcinc)))
204                                        ("vcsteps"     = ,(s+ "soma." (hoc-cname sysname vcsteps)))
205                                        ("vchdur"      = ,(s+ "soma." (hoc-cname sysname vchdur)))
206                                        ("vcbdur"      = ,(s+ "soma." (hoc-cname sysname vcbdur)))
207                                        )
208                                   
209                                    (pp indent
210                                        ("objref " ,name)
211                                        (,name = "new Vector(vec_sizes)")
212                                        (,(s+ name ".record") ,(s+ "(&soma.i_" name "_" sysname "( 0.5 ))"))
213                                        ("objref tlog")
214                                        ("tlog = new Vector(vec_sizes,0)")
215                                        ("tlog.record (&t)"))
216                                   
217                                    (pp indent ("objref logfile")
218                                        ("logfile = new File()")
219                                        ("logfile.wopen (" ,(s+ "\"" name ".dat\"") ")"))
220                                   
221                                    (pp indent ("vcrun(" ,i ", " ,name  ", tlog, logfile)")
222                                        ("g[" ,i "].label(.5,.85,\"" ,name "\")")
223                                        "logfile.close()")
224                                    )))
225                            vc-components 
226                            (list-tabulate (length vc-components) (lambda (i) i)))
227                           
228                           (if (> vcn 1) (pp indent ("for i=0, " ,(- vcn 1) " {")) (pp indent "i=0"))
229                           (pp (if (> vcn 1) indent+ indent) ("g[i].exec_menu(\"View = plot\")"))
230                           (if (> vcn 1) (pp indent ("}")))
231                           (pp indent ("vgraphbox.map()"))
232                           
233                           ))
234                       ))
235
236          ((matlab) 
237           (with-output-to-file filename
238            (lambda ()
239              (let* ((str          (lambda (x) (s+ "\"" x "\"")))
240                     (defs         (s+ sysname "_defs"))
241                     (init         (s+ sysname "_init"))
242                     (current-vars (map (lambda (comp) (s+ "i_" (first comp))) vc-components))
243                     (vcn          (length vc-components))
244                     (hboxn        (inexact->exact (ceiling (/ vcn 3))))
245                     )
246
247           (pp indent (,defs = ,(str (s+ sysname ".m")) #\;)
248               (autoload (,(str sysname ) #\, ,defs ) #\;)
249               (autoload (,(str init) #\, ,defs ) #\;))
250
251           (pp indent (,init #\;))
252           (pp indent (y0 = ,init (-65) #\;))
253
254
255           (pp indent (global N #\;)
256                      (N = length(y0) #\;))
257
258           (pp indent (function is = ,(s+ sysname "_currents_odepkg (y)")))
259           (pp indent+ (global N #\;)
260                       (global ,@(intersperse current-vars #\space) #\;))
261
262           (pp indent+ (t = y(1) #\;)
263                       (s = y(2:N+1) #\;)
264                       (y1 = ,sysname (t #\, s) #\;))
265
266           (pp indent+ (is = #\[ ,@(cons 't current-vars) #\] #\;))
267
268           (pp indent (endfunction))
269
270           (pp indent (function is = ,(s+ sysname "_currents_lsode (y)")))
271           (pp indent+ (global N #\;)
272                       (global ,@(intersperse current-vars #\space) #\;))
273
274           (pp indent+ (t = y(1) #\;)
275                       (s = y(2:N+1) #\;)
276                       (y1 = ,sysname (s #\, t) #\;))
277
278           (pp indent+ (is = #\[ ,@(cons 't current-vars) #\] #\;))
279
280           (pp indent (endfunction))
281
282           (pp indent (function ys = ,(s+ sysname "_vclamp_run1_odepkg") (N #\, dt #\, vchold #\, vchdur #\, vcbase #\, vcdur)))
283
284           (pp indent+ (global reltol abstol))
285
286           (pp indent+ (P = odeset (,(str "RelTol") #\, reltol #\, 
287                                    ,(str "AbsTol") #\, abstol #\, 
288                                    ,(str "MaxStep") #\, 1 #\, 
289                                    ,(str "InitialStep") #\, dt) #\;))
290
291           (pp indent+ (t0 = 0.0 #\;)
292                       (t1 = t0+vchdur #\;)
293                       (t2 = t1 #\;)
294                       (t3 = t2+vcdur #\;)
295                       (t4 = t3 #\;)
296                       (t5 = t4+vchdur #\;))
297
298           (pp indent+ (vinit = vchold #\;)
299                       (y0 = ,(s+ sysname "_init") (vinit))
300                       (sol = ode2r (,(s+ "@" sysname) #\, #\[ t0 t1 #\] #\, y0 #\, P) #\;)
301                       (ys = #\[ sol.x sol.y #\] #\;))
302
303           (pp indent+ (vinit = vcbase #\;)
304               (y0 = sol.y(size(sol.y)(1) #\, :)#\' #\; y0(1) = vinit)
305               (sol = ode2r(,(s+ "@" sysname) #\, #\[ t2 t3 #\] #\, y0 #\, P) #\;)
306               (ys = vertcat (ys #\, #\[ sol.x sol.y #\]) #\;))
307
308           (pp indent+ (vinit = vchold #\;)
309                       (y0 = sol.y(size(sol.y)(1) #\, :) #\' #\; y0(1) = vinit)
310                       (sol = ode2r(,(s+ "@" sysname) #\, #\[ t4 t5 #\] #\, y0 #\, P) #\;)
311                       (ys = vertcat (ys #\, #\[ sol.x sol.y #\]) #\;))
312
313           (pp indent (endfunction #\;))
314
315           (pp indent (function ys = ,(s+ sysname "_vclamp_run1_lsode") (N #\, dt #\, vchold #\, vchdur #\, vcbase #\, vcdur)))
316
317           (pp indent+ (global reltol abstol))
318
319           (pp indent+ (lsode_options (,(str "absolute tolerance") #\, abstol) #\;)
320                       (lsode_options (,(str "relative tolerance") #\, reltol) #\;)
321                       (lsode_options (,(str "integration method") #\, "bdf") #\;)
322                       (lsode_options (,(str "initial step size") #\, dt) #\;))
323
324           (pp indent+ (t0 = 0.0 #\;)
325                       (t1 = t0+vchdur #\;)
326                       (thold0 = linspace(t0 #\, t1 #\, (t1-t0)/dt) #\;)
327                       (t2 = t1 #\;)
328                       (t3 = t2+vcdur #\;)
329                       (tclamp = linspace(t2 #\, t3 #\, (t3-t2)/dt) #\;)
330                       (t4 = t3 #\;)
331                       (t5 = t4+vchdur #\;)
332                       (thold1 = linspace(t4 #\, t5 #\, (t5-t4)/dt) #\;)
333                       (vinit = vchold #\;))
334
335            (pp indent+ (y0 = ,(s+ sysname "_init") (vinit))
336                        (y = lsode(,(s+ "@" sysname) #\, y0 #\, thold0) #\;)
337                        (ys = #\[ "thold0'" y #\] #\;))
338
339            (pp indent+ (vinit = vcbase #\;)
340                        (y0 = y(size(y)(1) #\, :) #\' #\; y0(1) = vinit)
341                        (y = lsode(,(s+ "@" sysname) #\, y0 #\, tclamp) #\;)
342                        (ys = vertcat (ys #\, #\[ "tclamp'" y #\]) #\;))
343
344            (pp indent+ (vinit = vchold #\;)
345                        (y0 = y(size(y)(1) #\, :) #\' #\; y0(1) = vinit)
346                        (y = lsode(,(s+ "@" sysname) #\, y0 #\, thold1) #\;))
347
348            (pp indent+ (ys = vertcat (ys #\, #\[ "thold1'" y #\]) #\;))
349
350            (pp indent (endfunction))
351
352            (pp indent (function ilog = Khaliq03_vclamp (N #\, dt #\, vchold #\, vchdur #\, vcbase #\, vcdur #\, vcinc #\, vcsteps)))
353
354            (pp indent+ (vc = vcbase #\;)
355                        (ys = cell(1 #\, vcsteps) #\;))
356
357            (pp indent+ (for i = 1:vcsteps))
358            (pp (+ 2 indent+) 
359                (y = ,(s+ sysname "_vclamp_run1_odepkg") (N #\, dt #\, vchold #\, vchdur #\, vc #\, vcdur) #\;)
360                (ys(i) = y #\;)
361                (vc = vc + vcinc #\;))
362            (pp indent+ (endfor #\;))
363 
364            (pp indent+ (ilog = cell(vcsteps #\, 1) #\;))
365            (pp indent+ (for i = 1:vcsteps))
366            (pp (+ 2 indent+)
367                (ilogv = #\[#\] #\;)
368                (n = size(ys "{i}" )(1) #\;))
369            (pp (+ 2 indent+) (for j = 1:n))
370            (pp (+ 4 indent+) 
371                (next = Khaliq03_currents_odepkg (ys "{i}"(j #\, :)) #\;)
372                (ilogv = vertcat (ilogv #\, next) #\;))
373            (pp (+ 2 indent+) (endfor #\; ))
374            (pp (+ 2 indent+) (ilog "{i}" = ilogv #\;))
375            (pp indent+ (endfor #\;))
376     
377            (pp indent (endfunction))
378
379            (pp indent (global reltol abstol))
380
381            (pp indent (reltol = 1e-6 #\;))
382            (pp indent (abstol = 1e-6 #\;))
383            (pp indent (dt = 0.001 #\;))
384
385           (for-each
386            (lambda (comp i) 
387              (let* ((name (first comp))
388                     (sym  (third comp))
389                     (env  ((dis 'component-env) sys sym)))
390
391                (match-let (((vchold vcbase vcinc vcsteps vchdur vcbdur) ((dis 'component-exports) sys sym)))
392                           
393                   (pp indent ("##" ,(s+ "i_" name) plot))
394                   (pp indent (,(s+ "i_" name "_index") = ,(+ 1 i) #\;))
395
396                   (pp indent (global
397                               ;; vchold,vchdur,vcbase,vcdur,vcinc,vcsteps
398                               ,(matlab-name vchold) 
399                               ,(matlab-name vchdur) ,(matlab-name vcbase) 
400                               ,(matlab-name vcbdur) ,(matlab-name vcinc) 
401                               ,(matlab-name vcsteps)) )
402
403                   (pp indent (,(s+ sysname "_ilog") = ,(s+ sysname "_vclamp") 
404                               ;; vchold,vchdur,vcbase,vcdur,vcinc,vcsteps
405                               (N #\, dt #\, ,(matlab-name vchold) #\, 
406                                  ,(matlab-name vchdur) #\, ,(matlab-name vcbase) #\,
407                                  ,(matlab-name vcbdur) #\, ,(matlab-name vcinc) #\, 
408                                  ,(matlab-name vcsteps)) #\;))
409
410                   (pp indent (subplot(,hboxn #\, 3 #\, ,i) #\;))
411                   (pp indent (plot( ,@(intersperse
412                                        (concatenate
413                                         (list-tabulate (inexact->exact (const-val (environment-ref env vcsteps)))
414                                          (lambda (i) 
415                                           `(,(s+ sysname "_ilog{" (+ 1 i) "}(:,1)") 
416                                             ,(s+ (s+ sysname "_ilog{" (+ 1 i) "}")
417                                                  (s+ "(:,i_" name "_index)"))))))
418                                        #\, )
419                                    )
420                                   #\;))
421
422
423                   (pp indent (,(s+ name "_log") = vertcat ,(intersperse 
424                                       (list-tabulate (inexact->exact (const-val (environment-ref env vcsteps)))
425                                         (lambda (i) 
426                                           (s+ #\[ 
427                                               (s+ sysname "_ilog{" (+ 1 i) "}(:,1)") #\,
428                                               (s+ sysname "_ilog{" (+ 1 i) "}")
429                                               (s+ "(:,i_" name "_index)")
430                                               #\])))
431                                       #\, )
432                                     #\;))
433
434
435                   (pp indent (save -ascii ,(str (s+ name ".dat")) ,(s+ name "_log") #\;))
436
437                   )))
438                vc-components (list-tabulate vcn (lambda (i) (+ 1 i))))
439          ))))
440  )))))
441)))
442)
Note: See TracBrowser for help on using the repository browser.