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

Last change on this file since 23596 was 23596, checked in by Ivan Raikov, 10 years ago

nemo: removing lambda-lift declarations

File size: 14.0 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 (nemo:vclamp-translator sys . rest)
63
64  (define (cid x)  (second x))
65  (define (cn x)   (first x))
66 
67  (let-optionals rest ((target 'hoc) (filename #f))
68
69    (cases nemo:quantity (environment-ref sys (nemo-intern 'dispatch)) 
70
71         (DISPATCH  (dis) 
72
73          (let ((imports  ((dis 'imports)  sys))
74                (exports  ((dis 'exports)  sys)))
75           
76            (let* ((indent        0)
77                   (indent+       (+ 2 indent ))
78                   (eval-const    (dis 'eval-const))
79                   
80                   (sysname       (hoc-name ((dis 'sysname) sys)))
81                   (filename      (or filename (s+ sysname "_vclamp" (case target ((matlab octave) ".m") (else ".ses")))))
82
83                   (consts        ((dis 'consts)  sys))
84                   (asgns         ((dis 'asgns)   sys))
85                   (states        ((dis 'states)  sys))
86                   (reactions     ((dis 'reactions) sys))
87                   (rates         ((dis 'rates) sys))
88                   (defuns        ((dis 'defuns)  sys))
89                   (components    ((dis 'components) sys))
90                   
91                   (vc-components (filter-map (lambda (x) (and (eq? 'voltage-clamp (second x)) x)) components))
92                   (g             (match-let (((state-list asgn-list g) ((dis 'depgraph*) sys))) g))
93                   (poset         (vector->list ((dis 'depgraph->bfs-dist-poset) g)))
94                   
95                   (gate-complex-info    (nemo:gate-complex-query sys))
96                   (gate-complexes       (lookup-def 'gate-complexes gate-complex-info))
97                   
98                   (i-gates       (lookup-def 'i-gates gate-complex-info))
99                   
100                   )
101             
102              (if (not (null? vc-components))
103                 
104                  (case target 
105                   
106                    ((hoc)
107
108                     (with-output-to-file filename
109                       (lambda ()
110                         (let* ((vcn   (length vc-components))
111                                (hboxn (inexact->exact (ceiling (/ vcn 3)))))
112                       
113                           (pp indent
114                               ("load_file(\"nrngui.hoc\")")
115                               ("objref vce     // voltage clamp)")
116                               ("objref g[" ,vcn "] // graph objects"))
117
118                           (pp indent
119                               ("objref vgraphbox, hgraphbox[" ,hboxn "]")
120                               ("vgraphbox=new VBox()")
121                               ("vgraphbox.intercept(1)"))
122                           
123                           (let recur ((a 0) (b (min (- vcn 1) 2)) (hbi 0))
124                             (pp indent
125                                 ("hgraphbox[" ,hbi "]=new HBox()")
126                                 ("hgraphbox[" ,hbi "].intercept(1)"))
127                             
128                             (if (< a b)
129                                 (pp indent ("for i=" ,a "," ,b " {"))
130                                 (pp indent ("i = " ,a)))
131                             (pp (if (< a b) indent+ indent)
132                                 ("g[i]=new Graph()")
133                                 ("g[i].exec_menu(\"Keep Lines\")"))
134                             (if (< a b) (pp indent ("}")))
135                             (pp indent
136                                 ("hgraphbox[" ,hbi "].intercept(0)")
137                                 ("hgraphbox[" ,hbi "].map()"))
138                             (if (> hboxn (+ 1 hbi))
139                                 (recur (+ b 1) (min (- vcn 1) (+ b 3)) (+ 1 hbi)))
140                             )
141
142(pp indent (#<<EOF
143
144vgraphbox.intercept(0)
145
146vec_sizes = tstop/dt + 1        // recorded traces are all this size
147
148soma vce=new VClamp(0)
149
150dt = 0.001
151
152proc vcrun() {
153       
154        vce.dur[0]=vchdur
155        vce.dur[1]=vcbdur
156        vce.dur[2]=vchdur
157
158        tstop=vce.dur[0]+vce.dur[1]+vce.dur[2]
159       
160        vce.amp[0]=vchold
161        vce.amp[2]=vchold
162       
163        for j=0, vcsteps-1 {
164                x=vcbase+j*vcincrement
165                vce.amp[1]=x
166
167                init()
168
169                run()
170                $o2.line(g[$1], dt)
171
172                for i=0,$o2.size()-1 {
173                    $o4.printf("%g %g\n", $o3.x[i], $o2.x[i])
174                    }
175
176                if (stoppedrun()) {
177                        break
178                }
179        }
180}
181
182
183EOF
184))
185                    (for-each (lambda (comp i) 
186
187                       (let ((name (first comp))
188                             (sym  (third comp)))
189
190                         (match-let (((vchold vcbase vcinc vcsteps vchdur vcbdur) ((dis 'component-exports) sys sym)))
191
192                                    (pp indent ("print \"generating " ,name "\""))
193
194                                    (pp indent
195                                        ("vchold"      = ,(s+ "soma." (hoc-cname sysname vchold)))
196                                        ("vcbase"      = ,(s+ "soma." (hoc-cname sysname vcbase)))
197                                        ("vcincrement" = ,(s+ "soma." (hoc-cname sysname vcinc)))
198                                        ("vcsteps"     = ,(s+ "soma." (hoc-cname sysname vcsteps)))
199                                        ("vchdur"      = ,(s+ "soma." (hoc-cname sysname vchdur)))
200                                        ("vcbdur"      = ,(s+ "soma." (hoc-cname sysname vcbdur)))
201                                        )
202                                   
203                                    (pp indent
204                                        ("objref " ,name)
205                                        (,name = "new Vector(vec_sizes)")
206                                        (,(s+ name ".record") ,(s+ "(&soma.i_" name "_" sysname "( 0.5 ))"))
207                                        ("objref tlog")
208                                        ("tlog = new Vector(vec_sizes,0)")
209                                        ("tlog.record (&t)"))
210                                   
211                                    (pp indent ("objref logfile")
212                                        ("logfile = new File()")
213                                        ("logfile.wopen (" ,(s+ "\"" name ".dat\"") ")"))
214                                   
215                                    (pp indent ("vcrun(" ,i ", " ,name  ", tlog, logfile)")
216                                        ("g[" ,i "].label(.5,.85,\"" ,name "\")")
217                                        "logfile.close()")
218                                    )))
219                     vc-components (list-tabulate (length vc-components) (lambda (i) i)))
220           
221           (if (> vcn 1) (pp indent ("for i=0, " ,(- vcn 1) " {")) (pp indent "i=0"))
222           (pp (if (> vcn 1) indent+ indent) ("g[i].exec_menu(\"View = plot\")"))
223           (if (> vcn 1) (pp indent ("}")))
224           (pp indent ("vgraphbox.map()"))
225
226           ))
227          ))
228
229          ((matlab) 
230           (with-output-to-file filename
231            (lambda ()
232              (let* ((str          (lambda (x) (s+ "\"" x "\"")))
233                     (defs         (s+ sysname "_defs"))
234                     (init         (s+ sysname "_init"))
235                     (current-vars (map (lambda (comp) (s+ "i_" (first comp))) vc-components))
236                     (vcn          (length vc-components))
237                     (hboxn        (inexact->exact (ceiling (/ vcn 3))))
238                     )
239
240           (pp indent (,defs = ,(str (s+ sysname ".m")) #\;)
241               (autoload (,(str sysname ) #\, ,defs ) #\;)
242               (autoload (,(str init) #\, ,defs ) #\;))
243
244           (pp indent (,init #\;))
245           (pp indent (y0 = ,init (-65) #\;))
246
247
248           (pp indent (global N #\;)
249                      (N = length(y0) #\;))
250
251           (pp indent (function is = ,(s+ sysname "_currents_odepkg (y)")))
252           (pp indent+ (global N #\;)
253                       (global ,@(intersperse current-vars #\space) #\;))
254
255           (pp indent+ (t = y(1) #\;)
256                       (s = y(2:N+1) #\;)
257                       (y1 = ,sysname (t #\, s) #\;))
258
259           (pp indent+ (is = #\[ ,@(cons 't current-vars) #\] #\;))
260
261           (pp indent (endfunction))
262
263           (pp indent (function is = ,(s+ sysname "_currents_lsode (y)")))
264           (pp indent+ (global N #\;)
265                       (global ,@(intersperse current-vars #\space) #\;))
266
267           (pp indent+ (t = y(1) #\;)
268                       (s = y(2:N+1) #\;)
269                       (y1 = ,sysname (s #\, t) #\;))
270
271           (pp indent+ (is = #\[ ,@(cons 't current-vars) #\] #\;))
272
273           (pp indent (endfunction))
274
275           (pp indent (function ys = ,(s+ sysname "_vclamp_run1_odepkg") (N #\, dt #\, vchold #\, vchdur #\, vcbase #\, vcdur)))
276
277           (pp indent+ (global reltol abstol))
278
279           (pp indent+ (P = odeset (,(str "RelTol") #\, reltol #\, 
280                                    ,(str "AbsTol") #\, abstol #\, 
281                                    ,(str "MaxStep") #\, 1 #\, 
282                                    ,(str "InitialStep") #\, dt) #\;))
283
284           (pp indent+ (t0 = 0.0 #\;)
285                       (t1 = t0+vchdur #\;)
286                       (t2 = t1 #\;)
287                       (t3 = t2+vcdur #\;)
288                       (t4 = t3 #\;)
289                       (t5 = t4+vchdur #\;))
290
291           (pp indent+ (vinit = vchold #\;)
292                       (y0 = ,(s+ sysname "_init") (vinit))
293                       (sol = ode2r (,(s+ "@" sysname) #\, #\[ t0 t1 #\] #\, y0 #\, P) #\;)
294                       (ys = #\[ sol.x sol.y #\] #\;))
295
296           (pp indent+ (vinit = vcbase #\;)
297               (y0 = sol.y(size(sol.y)(1) #\, :)#\' #\; y0(1) = vinit)
298               (sol = ode2r(,(s+ "@" sysname) #\, #\[ t2 t3 #\] #\, y0 #\, P) #\;)
299               (ys = vertcat (ys #\, #\[ sol.x sol.y #\]) #\;))
300
301           (pp indent+ (vinit = vchold #\;)
302                       (y0 = sol.y(size(sol.y)(1) #\, :) #\' #\; y0(1) = vinit)
303                       (sol = ode2r(,(s+ "@" sysname) #\, #\[ t4 t5 #\] #\, y0 #\, P) #\;)
304                       (ys = vertcat (ys #\, #\[ sol.x sol.y #\]) #\;))
305
306           (pp indent (endfunction #\;))
307
308           (pp indent (function ys = ,(s+ sysname "_vclamp_run1_lsode") (N #\, dt #\, vchold #\, vchdur #\, vcbase #\, vcdur)))
309
310           (pp indent+ (global reltol abstol))
311
312           (pp indent+ (lsode_options (,(str "absolute tolerance") #\, abstol) #\;)
313                       (lsode_options (,(str "relative tolerance") #\, reltol) #\;)
314                       (lsode_options (,(str "integration method") #\, "bdf") #\;)
315                       (lsode_options (,(str "initial step size") #\, dt) #\;))
316
317           (pp indent+ (t0 = 0.0 #\;)
318                       (t1 = t0+vchdur #\;)
319                       (thold0 = linspace(t0 #\, t1 #\, (t1-t0)/dt) #\;)
320                       (t2 = t1 #\;)
321                       (t3 = t2+vcdur #\;)
322                       (tclamp = linspace(t2 #\, t3 #\, (t3-t2)/dt) #\;)
323                       (t4 = t3 #\;)
324                       (t5 = t4+vchdur #\;)
325                       (thold1 = linspace(t4 #\, t5 #\, (t5-t4)/dt) #\;)
326                       (vinit = vchold #\;))
327
328            (pp indent+ (y0 = ,(s+ sysname "_init") (vinit))
329                        (y = lsode(,(s+ "@" sysname) #\, y0 #\, thold0) #\;)
330                        (ys = #\[ "thold0'" y #\] #\;))
331
332            (pp indent+ (vinit = vcbase #\;)
333                        (y0 = y(size(y)(1) #\, :) #\' #\; y0(1) = vinit)
334                        (y = lsode(,(s+ "@" sysname) #\, y0 #\, tclamp) #\;)
335                        (ys = vertcat (ys #\, #\[ "tclamp'" y #\]) #\;))
336
337            (pp indent+ (vinit = vchold #\;)
338                        (y0 = y(size(y)(1) #\, :) #\' #\; y0(1) = vinit)
339                        (y = lsode(,(s+ "@" sysname) #\, y0 #\, thold1) #\;))
340
341            (pp indent+ (ys = vertcat (ys #\, #\[ "thold1'" y #\]) #\;))
342
343            (pp indent (endfunction))
344
345            (pp indent (function ilog = Khaliq03_vclamp (N #\, dt #\, vchold #\, vchdur #\, vcbase #\, vcdur #\, vcinc #\, vcsteps)))
346
347            (pp indent+ (vc = vcbase #\;)
348                        (ys = cell(1 #\, vcsteps) #\;))
349
350            (pp indent+ (for i = 1:vcsteps))
351            (pp (+ 2 indent+) 
352                (y = ,(s+ sysname "_vclamp_run1_odepkg") (N #\, dt #\, vchold #\, vchdur #\, vc #\, vcdur) #\;)
353                (ys(i) = y #\;)
354                (vc = vc + vcinc #\;))
355            (pp indent+ (endfor #\;))
356 
357            (pp indent+ (ilog = cell(vcsteps #\, 1) #\;))
358            (pp indent+ (for i = 1:vcsteps))
359            (pp (+ 2 indent+)
360                (ilogv = #\[#\] #\;)
361                (n = size(ys "{i}" )(1) #\;))
362            (pp (+ 2 indent+) (for j = 1:n))
363            (pp (+ 4 indent+) 
364                (next = Khaliq03_currents_odepkg (ys "{i}"(j #\, :)) #\;)
365                (ilogv = vertcat (ilogv #\, next) #\;))
366            (pp (+ 2 indent+) (endfor #\; ))
367            (pp (+ 2 indent+) (ilog "{i}" = ilogv #\;))
368            (pp indent+ (endfor #\;))
369     
370            (pp indent (endfunction))
371
372            (pp indent (global reltol abstol))
373
374            (pp indent (reltol = 1e-6 #\;))
375            (pp indent (abstol = 1e-6 #\;))
376            (pp indent (dt = 0.001 #\;))
377
378           (for-each
379            (lambda (comp i) 
380              (let* ((name (first comp))
381                     (sym  (third comp))
382                     (env  ((dis 'component-env) sys sym)))
383
384                (match-let (((vchold vcbase vcinc vcsteps vchdur vcbdur) ((dis 'component-exports) sys sym)))
385                           
386                   (pp indent ("##" ,(s+ "i_" name) plot))
387                   (pp indent (,(s+ "i_" name "_index") = ,(+ 1 i) #\;))
388
389                   (pp indent (global
390                               ;; vchold,vchdur,vcbase,vcdur,vcinc,vcsteps
391                               ,(matlab-name vchold) 
392                               ,(matlab-name vchdur) ,(matlab-name vcbase) 
393                               ,(matlab-name vcbdur) ,(matlab-name vcinc) 
394                               ,(matlab-name vcsteps)) )
395
396                   (pp indent (,(s+ sysname "_ilog") = ,(s+ sysname "_vclamp") 
397                               ;; vchold,vchdur,vcbase,vcdur,vcinc,vcsteps
398                               (N #\, dt #\, ,(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 (subplot(,hboxn #\, 3 #\, ,i) #\;))
404                   (pp indent (plot( ,@(intersperse
405                                        (concatenate
406                                         (list-tabulate (inexact->exact (const-val (environment-ref env vcsteps)))
407                                          (lambda (i) 
408                                           `(,(s+ sysname "_ilog{" (+ 1 i) "}(:,1)") 
409                                             ,(s+ (s+ sysname "_ilog{" (+ 1 i) "}")
410                                                  (s+ "(:,i_" name "_index)"))))))
411                                        #\, )
412                                    )
413                                   #\;))
414
415
416                   (pp indent (,(s+ name "_log") = vertcat ,(intersperse 
417                                       (list-tabulate (inexact->exact (const-val (environment-ref env vcsteps)))
418                                         (lambda (i) 
419                                           (s+ #\[ 
420                                               (s+ sysname "_ilog{" (+ 1 i) "}(:,1)") #\,
421                                               (s+ sysname "_ilog{" (+ 1 i) "}")
422                                               (s+ "(:,i_" name "_index)")
423                                               #\])))
424                                       #\, )
425                                     #\;))
426
427
428                   (pp indent (save -ascii ,(str (s+ name ".dat")) ,(s+ name "_log") #\;))
429
430                   )))
431                vc-components (list-tabulate vcn (lambda (i) (+ 1 i))))
432          ))))
433  )))))
434)))
435)
Note: See TracBrowser for help on using the repository browser.