source: project/release/4/nemo/trunk/nemo-nmodl.scm @ 27113

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

nemo: bug fixes and error handling

File size: 37.2 KB
Line 
1;;       
2;;
3;; An extension for translating NEMO models to NMODL descriptions.
4;;
5;; Copyright 2008-2012 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-nmodl
22
23        (nemo:nmodl-translator)
24
25        (import scheme chicken utils data-structures lolevel srfi-1 srfi-13 srfi-69)
26
27        (require-extension lolevel datatype matchable strictly-pretty 
28                           varsubst datatype 
29                           nemo-core nemo-utils nemo-gate-complex)
30
31(define nmodl-builtin-consts
32  `(celsius diam))
33
34(define nmodl-ops
35  `(+ - * / > < <= >= = ^))
36
37(define builtin-fns
38  `(+ - * / pow neg abs atan asin acos sin cos exp ln
39      sqrt tan cosh sinh tanh hypot gamma lgamma log10 log2 log1p ldexp cube
40      > < <= >= = and or round ceiling floor max min
41      ))
42
43
44(define (nmodl-name s)
45  (let ((cs (string->list (->string s))))
46    (let loop ((lst (list)) (cs cs))
47      (if (null? cs) (string->symbol (list->string (reverse lst)))
48          (let* ((c (car cs))
49                 (c1 (cond ((or (char-alphabetic? c) (char-numeric? c) (char=? c #\_)) c)
50                           (else #\_))))
51            (loop (cons c1 lst) (cdr cs)))))))
52                           
53
54(define (nmodl-state-name n s)
55  (nmodl-name (if n (s+ n s) s)))
56
57(define (rhsvars rhs)
58  (enum-freevars rhs (list) (list)))
59
60(define (rhsexpr/NMODL expr)
61  (match expr 
62         (('if . es)  `(if . ,(map (lambda (x) (rhsexpr/NMODL x)) es)))
63         (('let bnds body) `(let ,(map (lambda (x) (list (car x) (rhsexpr/NMODL (cadr x)))) bnds) ,(rhsexpr/NMODL body)))
64         (('pow x y)  (if (and (integer? y)  (positive? y))
65                          (if (> y 1)  (let ((tmp (gensym "x")))
66                                         `(let ((,tmp ,x)) (* . ,(list-tabulate (inexact->exact y) (lambda (i) tmp)))))
67                              x)
68                          (if (and (number? y) (zero? y)) 1.0 expr)))
69         ((s . es)    (if (symbol? s)   (cons (if (member s builtin-fns) s (nmodl-name s)) (map (lambda (x) (rhsexpr/NMODL x)) es)) expr))
70         (id          (if (symbol? id) (nmodl-name id) id))))
71
72
73(define-syntax pp
74  (syntax-rules ()
75    ((pp indent val ...) (ppf indent (quasiquote val) ...))))
76
77
78(define (letblk/NMODL e1 e2)
79  (cond ((equal? e1 (doc:empty)) (doc:group (doc:nest 2 e2)))
80        ((equal? e2 (doc:empty)) (doc:group (doc:nest 2 e1)))
81        (else (doc:connect (doc:group (doc:nest 2 e1))
82                           (doc:group (doc:nest 2 e2))))))
83       
84(define ifthen/NMODL  (doc:ifthen 0 (doc:text "if") (doc:text "") (doc:text "else")))
85(define group/NMODL   (doc:block 2 (doc:text "(") (doc:text ")")))
86(define block/NMODL   (doc:block 2 (doc:text "{") (doc:text "}")))
87(define binop/NMODL   (doc:binop 2))
88
89(define (format-op/NMODL indent op args)
90  (let ((op1 (doc:text (->string op))))
91    (let ((res
92           (if (null? args) op1
93               (match args
94                      ((x)           (doc:connect op1 x))
95                      ((x y)         (binop/NMODL x op1 y))
96                      ((x y z)       (binop/NMODL x op1 (binop/NMODL y op1 z)))
97                      (lst           (let* ((n   (length lst))
98                                            (n/2 (inexact->exact (round (/ n 2)))))
99                                       (binop/NMODL (format-op/NMODL indent op (take lst n/2 )) op1 
100                                                    (format-op/NMODL indent op (drop lst n/2 )))))))))
101      res)))
102
103
104(define (format-conseq-op/NMODL indent op args)
105  (let ((op1 (doc:text (->string op))))
106    (if (null? args) op1
107        (match args
108               ((x)      (doc:concat (list op1 x)))
109               ((x y)    (doc:concat (intersperse (list x op1 y) (doc:space))))
110               ((x y z)  (doc:concat (intersperse (list x op1 y op1 z) (doc:space))))
111               (lst      (let* ((n   (length lst))
112                                (n/2 (inexact->exact (round (/ n 2)))))
113                           (doc:concat 
114                            (intersperse 
115                             (list (format-conseq-op/NMODL indent op (take lst n/2 )) op1 
116                                   (format-conseq-op/NMODL indent op (drop lst n/2 )))
117                             (doc:space)))))))))
118
119(define (format-fncall/NMODL indent op args)
120  (let ((op1 (doc:text (->string op))))
121    (doc:cons op1 (group/NMODL ((doc:list indent identity (lambda () (doc:text ", "))) args)))))
122
123(define (name-normalize expr)
124  (match expr 
125         (('if c t e)  `(if ,(name-normalize c) ,(name-normalize t) ,(name-normalize e)))
126         (('let bs e)
127          `(let ,(map (lambda (b) `(,(car b) ,(name-normalize (cadr b)))) bs) ,(name-normalize e)))
128         ((f . es) 
129          (cons f (map name-normalize es)))
130         ((? symbol? ) (nmodl-name expr))
131         ((? atom? ) expr)))
132
133
134(define (canonicalize-expr/NMODL expr)
135  (let ((subst-convert  (subst-driver (lambda (x) (and (symbol? x) x)) 
136                                      nemo:binding? identity nemo:bind 
137                                      nemo:subst-term)))
138    (let* ((expr1 (if-convert expr))
139           (expr2 (subst-convert expr1 subst-empty))
140           (expr3 (let-lift expr2))
141           (expr4 (name-normalize expr3)))
142      expr4)))
143
144
145(define (format-expr/NMODL indent expr . rest) 
146  (let-optionals rest ((rv #f))
147   (let ((indent+ (+ 2 indent)))
148    (match expr
149       (('let bindings body)
150        (letblk/NMODL
151         (fold-right 
152           (lambda (x ax)
153             (let ((res
154                    (letblk/NMODL
155                     (match (second x)
156                            (('if c t e)
157                             (ifthen/NMODL
158                              (group/NMODL (format-expr/NMODL indent c))
159                              (block/NMODL (format-expr/NMODL indent t (first x)))
160                              (block/NMODL (format-expr/NMODL indent e (first x)))))
161                            (else
162                             (format-op/NMODL indent+ " = "
163                                              (list (format-expr/NMODL indent (first x) )
164                                                    (format-expr/NMODL indent (second x))))))
165                     ax)))
166               res
167               ))
168           (doc:empty) bindings)
169         (match body
170                (('let _ _) (format-expr/NMODL indent body rv))
171                (else
172                 (let ((body1 (doc:nest indent (format-expr/NMODL indent body))))
173                   (if rv  (format-op/NMODL indent " = " (list (format-expr/NMODL indent+ rv ) body1))
174                       body1))))))
175       
176       (('if . rest) (error 'format-expr/NMODL "invalid if statement " expr))
177
178       ((op . rest) 
179       (let ((op (case op ((pow)  '^) ((abs) 'fabs) ((ln) 'log) (else op))))
180         (let ((fe
181                (if (member op nmodl-ops)
182                    (let ((mdiv?  (any (lambda (x) (match x (('* . _) #t) (('/ . _) #t) (else #f))) rest))
183                          (mul?   (any (lambda (x) (match x (('* . _) #t) (else #f))) rest))
184                          (plmin? (any (lambda (x) (match x (('+ . _) #t) (('- . _) #t) (else #f))) rest)))
185                      (case op
186                        ((/) 
187                         (format-op/NMODL indent op 
188                                          (map (lambda (x) 
189                                                 (let ((fx (format-expr/NMODL indent+ x)))
190                                                   (if (or (symbol? x) (number? x)) fx
191                                                       (if (or mul? plmin?) (group/NMODL fx) fx)))) rest)))
192                        ((*) 
193                         (format-op/NMODL indent op 
194                                          (map (lambda (x) 
195                                                 (let ((fx (format-expr/NMODL indent+ x)))
196                                                   (if (or (symbol? x) (number? x)) fx
197                                                       (if plmin? (group/NMODL fx) fx)))) rest)))
198                       
199                        ((^) 
200                         (format-op/NMODL indent op 
201                                          (map (lambda (x) 
202                                                 (let ((fx (format-expr/NMODL indent+ x)))
203                                                   (if (or (symbol? x)  (number? x)) fx
204                                                       (if (or mdiv? plmin?) (group/NMODL fx) fx)))) rest)))
205                       
206                        (else
207                         (format-op/NMODL indent op 
208                                          (map (lambda (x) 
209                                                 (let ((fx (format-expr/NMODL indent+ x))) fx)) rest)))))
210                   
211                    (let ((op (case op ((neg) '-) (else op))))
212                      (format-fncall/NMODL indent op (map (lambda (x) (format-expr/NMODL indent+ x)) rest))))))
213           (if rv 
214               (format-op/NMODL indent " = " (list (format-expr/NMODL indent+ rv ) fe))
215               fe))))
216     
217      (else  (let ((fe (doc:text (->string expr))))
218               (if rv 
219                   (format-op/NMODL indent " = " (list (format-expr/NMODL indent+ rv ) fe))
220                   fe)))))))
221               
222
223         
224(define (expr->string/NMODL x . rest)
225  (let-optionals rest ((rv #f) (width 72))
226    (sdoc->string (doc:format width (format-expr/NMODL 2 x rv)))))
227 
228
229(define (format-conseq/NMODL indent expr . rest) 
230  (let-optionals rest ((rv #f))
231   (let ((indent+ (+ 2 indent)))
232    (match expr
233       (('let bindings body)
234        (letblk/NMODL
235         (fold-right 
236           (lambda (x ax)
237             (letblk/NMODL
238              (match (second x)
239                     (('if c t e)
240                      (ifthen/NMODL
241                       (group/NMODL (format-conseq/NMODL indent c))
242                       (block/NMODL (format-conseq/NMODL indent t (first x)))
243                       (block/NMODL (format-conseq/NMODL indent e (first x)))))
244                     (else
245                      (format-conseq-op/NMODL indent+ " = "
246                                       (list (format-conseq/NMODL indent (first x) )
247                                             (format-conseq/NMODL indent (second x))))))
248              ax))
249           (doc:empty) bindings)
250         (let ((body1 (doc:nest indent (format-conseq/NMODL indent body))))
251           (if rv  (format-conseq-op/NMODL indent " = " (list (format-conseq/NMODL indent+ rv ) body1))
252               body1))))
253       
254       (('if . rest) (error 'format-conseq/NMODL "invalid if statement " expr))
255
256       ((op . rest) 
257        (let ((op (case op ((pow)  '^) ((abs) 'fabs) (else op))))
258          (let ((fe
259                (if (member op nmodl-ops)
260                    (let ((mdiv?  (any (lambda (x) (match x (('* . _) #t) (('/ . _) #t) (else #f))) rest))
261                          (mul?   (any (lambda (x) (match x (('* . _) #t) (else #f))) rest))
262                          (plmin? (any (lambda (x) (match x (('+ . _) #t) (('- . _) #t) (else #f))) rest)))
263                      (case op
264                        ((/) 
265                         (format-conseq-op/NMODL indent op 
266                                          (map (lambda (x) 
267                                                 (let ((fx (format-conseq/NMODL indent+ x)))
268                                                   (if (or (symbol? x) (number? x)) fx
269                                                       (if (or mul? plmin?) (group/NMODL fx) fx)))) rest)))
270                        ((*) 
271                         (format-conseq-op/NMODL indent op 
272                                          (map (lambda (x) 
273                                                 (let ((fx (format-conseq/NMODL indent+ x)))
274                                                   (if (or (symbol? x) (number? x)) fx
275                                                       (if plmin? (group/NMODL fx) fx)))) rest)))
276                       
277                        ((^) 
278                         (format-conseq-op/NMODL indent op 
279                                          (map (lambda (x) 
280                                                 (let ((fx (format-conseq/NMODL indent+ x)))
281                                                   (if (or (symbol? x)  (number? x)) fx
282                                                       (if (or mdiv? plmin?) (group/NMODL fx) fx)))) rest)))
283                       
284                        (else
285                         (format-conseq-op/NMODL indent op 
286                                          (map (lambda (x) 
287                                                 (let ((fx (format-conseq/NMODL indent+ x))) fx)) rest)))))
288                   
289                    (case op
290                      ((neg) (format-conseq-op/NMODL indent '* (map (lambda (x) (format-conseq/NMODL indent+ x)) 
291                                                             (cons "(-1)" rest))))
292                      (else  (format-fncall/NMODL indent op (map (lambda (x) (format-conseq/NMODL indent+ x)) 
293                                                                 rest)))))))
294
295           (if rv (format-conseq-op/NMODL indent " = " (list (format-conseq/NMODL indent+ rv ) fe)) fe))))
296     
297      (else  (let ((fe (doc:text (->string expr))))
298               (if rv 
299                   (format-conseq-op/NMODL indent " = " (list (format-conseq/NMODL indent+ rv ) fe))
300                   fe)))))))
301               
302
303(define (conserve-conseq->string/NMODL x val . rest)
304  (let-optionals rest ((width 72))
305    (s+ "CONSERVE " (sdoc->string (doc:format width (format-conseq/NMODL 2 x #f))) 
306        " = " (number->string val))))
307 
308
309(define (make-define-fn table? min-v max-v with depend)
310  (lambda (indent n proc)
311    (let ((lst (procedure-data proc))
312          (indent+ (+ 2 indent)))
313      (let ((rt       (lookup-def 'rt lst))
314            (formals  (lookup-def 'formals lst))
315            (vars     (lookup-def 'vars lst))
316            (body     (lookup-def 'body lst)))
317        (pp indent ,nl (FUNCTION ,(nmodl-name n) (,(slp ", " vars)) "{" ))
318        (let* ((body0 (rhsexpr/NMODL body))
319               (body1 (canonicalize-expr/NMODL body0))
320               (lbs   (enum-bnds body1 (list))))
321          (if (not (null? lbs)) (pp indent+ (LOCAL ,(slp ", " lbs))))
322          (if (and table? min-v max-v with)
323              (match vars
324                     (('v)  (pp indent+ (TABLE ,@(if depend `(DEPEND ,depend) `(""))
325                                               FROM ,min-v TO ,max-v WITH ,with)))
326                     (else  (void))))
327          (pp indent+ ,(expr->string/NMODL body1 (nmodl-name n))))
328        (pp indent "}"))) 
329    ))
330
331(define (expeuler dt name rhs)
332  (define (isname? x) (equal? x name))
333  (let ((res 
334         (match rhs
335                ((or ('- A ('* B (and x (? isname?))))
336                     ('+ ('neg ('* B (and x (? isname?)))) A))
337                 (let ((xexp (string->symbol (s+ x 'exp))))
338                   `(let ((,xexp (exp (* (neg ,B) ,dt))))
339                      (+ (* ,x ,xexp)  (* (- 1 ,xexp) (/ ,A ,B))))))
340
341                ((or ('- A ('* (and x (? isname?)) . B))
342                     ('+ ('neg ('* (and x (? isname?)) . B)) A))
343                 (let ((xexp (string->symbol (s+ x 'exp)))
344                       (B1   (if (null? (cdr B)) (car B) `(* ,@B))))
345                   `(let ((,xexp (exp (* (neg ,B1) ,dt))))
346                      (+ (* ,x ,xexp)  (* (- 1 ,xexp) (/ ,A ,B1))))))
347               
348                (('+ ('neg ('* (and x1 (? isname?)) Alpha))
349                     ('* ('- 1 (and x2 (? isname?))) Beta))
350                 (let ((A  Alpha)
351                       (B  `(+ ,Alpha ,Beta)))
352                   (let ((xexp (string->symbol (s+ x1 'exp))))
353                     `(let ((,xexp (exp (* (neg ,B) ,dt))))
354                        (+ (* ,x1 ,xexp) (* (- 1 ,xexp) (/ ,A ,B)))))))
355               
356                (('let bnds body)
357                 `(let ,bnds ,(expeuler dt name body)))
358               
359                (else (nemo:error 'nemo:expeuler ": unable to rewrite equation " rhs 
360                                  "in exponential Euler form")))))
361
362    res))
363 
364
365(define (reaction-transition-eqs n initial open transitions conserve power method)
366  (match-let (((g cnode node-subs)  (transitions-graph n open transitions conserve nmodl-state-name)))
367     (let* ((out-edges  (g 'out-edges))
368            (in-edges   (g 'in-edges))
369            (nodes      ((g 'nodes))))
370       ;; generate differential equations for each state in the transitions system
371      (let ((eqs    (fold (lambda (s ax) 
372                            (if (and cnode (= (first cnode) (first s) )) ax
373                                (let* ((out   (out-edges (first s)))
374                                       (in    (in-edges (first s)))
375                                       (open? (eq? (second s) open))
376                                       (name  (nmodl-name (lookup-def (second s) node-subs))))
377                                  (let* ((rhs1  (cond ((and (not (null? out)) (not (null? in)))
378                                                       `(+ (neg ,(sum (map third out)))
379                                                           ,(sum (map third in))))
380                                                      ((and (not (null? out)) (null? in))
381                                                       `(neg ,(sum (map third out))))
382                                                      ((and (null? out) (not (null? in)))
383                                                       (sum (map third in)))))
384                                         (fbody0 (rhsexpr/NMODL rhs1))
385                                         (fbody1 (case method
386                                                   ((expeuler) (canonicalize-expr/NMODL (expeuler 'dt name fbody0)))
387                                                   (else       (canonicalize-expr/NMODL fbody0)))))
388                                    (cons (list name fbody1) ax))
389                                  )))
390                          (list) nodes)))
391        eqs))))
392
393
394           
395       
396
397(define (reaction-keqs n initial open transitions power)
398  (let* ((subst-convert  (subst-driver (lambda (x) (and (symbol? x) x)) 
399                                       nemo:binding? identity nemo:bind nemo:subst-term))
400         (state-list     (let loop ((lst (list)) (tlst transitions))
401                           (if (null? tlst)  (delete-duplicates lst eq?)
402                               (match (car tlst) 
403                                      (('-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr)
404                                       (loop (cons* s0 s1 lst) (cdr tlst)))
405                                      (((and (? symbol?) s0) '-> (and (? symbol? s1)) rate-expr)
406                                       (loop (cons* s0 s1 lst) (cdr tlst)))
407                                      (('<-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr1 rate-expr2)
408                                       (loop (cons* s0 s1 lst) (cdr tlst)))
409                                      (((and (? symbol?) s0) 'M-> (and (? symbol? s1)) rate-expr1 rate-expr2)
410                                       (loop (cons* s0 s1 lst) (cdr tlst)))
411                                      (else
412                                       (nemo:error 'nemo:nmodl-reaction-keqs ": invalid transition equation " 
413                                                   (car tlst) " in state complex " n))
414                                      (else (loop lst (cdr tlst)))))))
415         (state-subs     (fold (lambda (s ax) (subst-extend s (nmodl-state-name n s) ax)) subst-empty state-list)))
416    ;; generate kinetic equations for each edge in the transitions system
417    (list n 
418          (map
419           (lambda (e) 
420             (match e
421                    (('-> s0 s1 rexpr)
422                     (let ((i  (lookup-def s0 state-subs))
423                           (j  (lookup-def s1 state-subs)))
424                       `(-> ,i ,j ,(canonicalize-expr/NMODL 
425                                    (subst-convert rexpr state-subs)))))
426                   
427                    ((s0 '-> s1 rexpr)
428                     (let ((i  (lookup-def s0 state-subs))
429                           (j  (lookup-def s1 state-subs)))
430                       `(-> ,i ,j ,(canonicalize-expr/NMODL 
431                                    (subst-convert rexpr state-subs)))))
432                   
433                    (('<-> s0 s1 rexpr1 rexpr2)
434                     (let ((i  (lookup-def s0 state-subs))
435                           (j  (lookup-def s1 state-subs)))
436                       `(<-> ,i ,j 
437                             ,(canonicalize-expr/NMODL (subst-convert rexpr1 state-subs))
438                             ,(canonicalize-expr/NMODL (subst-convert rexpr2 state-subs)))))
439                   
440                    ((s0 '<-> s1 rexpr1 rexpr2)
441                     (let ((i  (lookup-def s0 state-subs))
442                           (j  (lookup-def s1 state-subs)))
443                       `(<-> ,i ,j 
444                             ,(canonicalize-expr/NMODL (subst-convert rexpr1 state-subs))
445                             ,(canonicalize-expr/NMODL (subst-convert rexpr2 state-subs)))))
446                   
447                   
448                    (else (nemo:error 'nemo:nmodl-reaction-keqs ": invalid transition equation " 
449                                      e " in state complex " n))))
450           transitions))))
451       
452
453
454(define (state-init n init)
455  (let* ((init  (rhsexpr/NMODL init))
456         (init1 (canonicalize-expr/NMODL init)))
457    (list (nmodl-name n) init1)))
458
459
460(define (asgn-eq n rhs)
461  (let* ((fbody   (rhsexpr/NMODL rhs))
462         (fbody1  (canonicalize-expr/NMODL fbody)))
463    (list (nmodl-name n) fbody1)))
464
465
466(define (reaction-eq n open transitions conserve)
467  (list (nmodl-name n) (nmodl-state-name n open)))
468
469(define (poset->reaction-eq-defs poset sys kinetic)
470  (fold-right
471   (lambda (lst ax)
472     (fold  (lambda (x ax) 
473              (match-let (((i . n)  x))
474                         (let ((en (hash-table-ref sys n)))
475                           (if (nemo:quantity? en)
476                               (cases nemo:quantity en
477                                      (REACTION  (name initial open transitions conserve power) 
478                                                 (cons (reaction-eq name open transitions conserve) ax))
479
480                                      (else  ax))
481                               ax))))
482            ax lst))
483   (list) poset))
484
485
486(define (poset->asgn-eq-defs poset sys)
487  (fold-right
488   (lambda (lst ax)
489     (fold  (lambda (x ax) 
490              (match-let (((i . n)  x))
491                         (let ((en (hash-table-ref sys n)))
492                           (if (nemo:quantity? en)
493                               (cases nemo:quantity en
494                                      (ASGN  (name value rhs) (cons (asgn-eq name rhs) ax))
495                                      (else  ax))
496                               ax))))
497            ax lst))
498   (list) poset))
499
500
501(define (poset->rate-eq-defs poset sys kinetic method)
502  (fold-right
503   (lambda (lst ax)
504     (fold  (lambda (x ax) 
505              (match-let (((i . n)  x))
506                         (let ((en (hash-table-ref sys n)))
507                           (if (and (not (member n kinetic)) (nemo:quantity? en))
508                               (cases nemo:quantity en
509                                      (REACTION  (name initial open transitions conserve power) 
510                                                 (append (reaction-transition-eqs name initial open transitions 
511                                                                                  conserve power method) ax))
512                                     
513                                      (RATE (name initial rhs power)
514                                            (let ((fbody0  (rhsexpr/NMODL rhs))
515                                                  (dy      (nmodl-name name )))
516                                              (case method
517                                                ((expeuler) 
518                                                 (cons (list dy (canonicalize-expr/NMODL (expeuler 'dt name fbody0))) 
519                                                       ax))
520                                                (else
521                                                 (cons (list dy (canonicalize-expr/NMODL fbody0)) ax)))))
522
523                                      (else  ax))
524                               ax))))
525            ax lst))
526   (list) poset))
527
528
529(define (poset->kinetic-eq-defs poset sys kinetic)
530  (fold-right
531   (lambda (lst ax)
532     (fold  (lambda (x ax) 
533              (match-let (((i . n)  x))
534                         (let ((en (hash-table-ref sys n)))
535                           (if (and (member n kinetic) (nemo:quantity? en))
536                               (cases nemo:quantity en
537                                      (REACTION  (name initial open transitions conserve power) 
538                                                 (cons (reaction-keqs name initial open transitions power) ax))
539                                      (else  ax))
540                               ax))))
541            ax lst))
542   (list) poset))
543
544
545(define (poset->state-init-defs poset sys)
546  (fold-right
547   (lambda (lst ax)
548     (fold  (lambda (x ax) 
549              (match-let (((i . n)  x))
550                         (let ((en (hash-table-ref sys n)))
551                           (if (nemo:quantity? en)
552                               (cases nemo:quantity en
553                                      (REACTION  (name initial open transitions conserve power)
554                                                 (if (nemo:rhs? initial)
555                                                     (cons* (state-init name initial) 
556                                                            (state-init (nmodl-state-name name open) name) ax) 
557                                                     ax))
558                                      (RATE  (name initial rhs power)
559                                             (if (nemo:rhs? initial)
560                                                 (cons (state-init name initial) ax) 
561                                                 ax))
562                                      (else  ax))
563                               ax))))
564            ax lst))
565   (list) poset))
566
567(define (poset->state-conserve-eq-defs poset sys)
568  (fold-right
569   (lambda (lst ax)
570     (fold  (lambda (x ax) 
571              (match-let (((i . n)  x))
572                         (let ((en (hash-table-ref sys n)))
573                           (if (nemo:quantity? en)
574                               (cases nemo:quantity en
575                                      (REACTION (name initial open transitions conserve power)
576                                                (if (and (list? conserve) (every nemo:conseq? conserve))
577                                                    (cons (state-conseqs (nmodl-name name) transitions conserve
578                                                                        nmodl-state-name) ax) 
579                                                    ax))
580                                      (else  ax))
581                               ax))))
582            ax lst))
583   (list) poset))
584
585
586(define (find-locals defs)
587  (concatenate (map (lambda (def) 
588                      (match def 
589                             (('let bnds body) 
590                              (let ((bexprs (map second bnds)))
591                                (concatenate (list (map first bnds) 
592                                                   (find-locals bexprs )
593                                                   (find-locals (list body))))))
594                             (('if c t e)      (append (find-locals (list t)) (find-locals (list e))))
595                             ((s . rest)       (find-locals rest))
596                             (else (list)))) 
597                    defs)))
598
599
600(define (rate/reaction-power sys n)
601  (let ((en (hash-table-ref sys n)))
602    (if (nemo:quantity? en)
603        (cases nemo:quantity en
604               (REACTION  (name initial open transitions conserve power) 
605                          power)
606               (RATE    (name initial rhs power)
607                        power)
608               (else  #f)) 
609        #f)))
610
611
612(define (bucket-partition p lst)
613  (let loop ((lst lst) (ax (list)))
614    (if (null? lst) ax
615        (let ((x (car lst)))
616          (let bkt-loop ((old-bkts ax) (new-bkts (list)))
617            (if (null? old-bkts) (loop (cdr lst) (cons (list x) new-bkts))
618                (if (p x (caar old-bkts ))
619                    (loop (cdr lst) (append (cdr old-bkts) (cons (cons x (car old-bkts)) new-bkts)))
620                    (bkt-loop (cdr old-bkts) (cons (car old-bkts) new-bkts)))))))))
621
622
623
624(define (nemo:nmodl-translator sys . rest)
625  (define (cid x)  (second x))
626  (define (cn x)   (first x))
627  (let-optionals rest ((method 'cnexp) (table? #f) (min-v -100) (max-v 100) (step 0.5) 
628                       (depend #f)  (kinetic (list)) (linear? #f))
629  (match-let ((($ nemo:quantity 'DISPATCH  dis) (hash-table-ref sys (nemo-intern 'dispatch))))
630    (let ((imports  ((dis 'imports)  sys))
631          (exports  ((dis 'exports)  sys)))
632      (let* ((indent        0)
633             (indent+       (+ 2 indent ))
634             (table-with    (and table? (inexact->exact (round (/ (abs (- max-v min-v)) step)))))
635             (sysname       (nmodl-name ((dis 'sysname) sys)))
636             (consts        ((dis 'consts)  sys))
637             (asgns         ((dis 'asgns)   sys))
638             (states        ((dis 'states)  sys))
639             (kinetic       (or kinetic '()))
640             (kinetic       (delete-duplicates
641                             (cond ((eq? kinetic 'all) (filter-map first states))
642                                   ((symbol? kinetic) 
643                                    (let ((sk (->string kinetic)))
644                                      (filter-map (lambda (s) (and s (and (string-suffix? sk (->string s)) s)) )
645                                                  (map first states))))
646                                   (else
647                                    (let ((kinetic (map ->string kinetic))
648                                          (ss      (map (compose ->string first) states)))
649                                      (concatenate
650                                       (map (lambda (sk)
651                                              (filter-map (lambda (s) (and (string-suffix? sk s) s))
652                                                          ss))
653                                            kinetic)))))))
654             (reactions     ((dis 'reactions) sys))
655             (rates         ((dis 'rates) sys))
656             (defuns        ((dis 'defuns)  sys))
657             (components    ((dis 'components) sys))
658             (g             (match-let (((state-list asgn-list g) ((dis 'depgraph*) sys))) g))
659             (poset         (vector->list ((dis 'depgraph->bfs-dist-poset) g)))
660             
661             (gate-complex-info    (nemo:gate-complex-query sys))
662             (gate-complexes       (lookup-def 'gate-complexes gate-complex-info))
663             (perm-ions     (map (match-lambda ((comp i e erev) `(,comp ,(nmodl-name i) ,(nmodl-name e) ,erev)))
664                                 (lookup-def 'perm-ions gate-complex-info)))
665             (acc-ions      (map (match-lambda ((comp i in out) `(,comp ,@(map nmodl-name (list i in out)))))
666                                 (lookup-def 'acc-ions gate-complex-info)))
667             (epools        (lookup-def 'pool-ions gate-complex-info))
668             (pool-ions     (map (lambda (lst) (map nmodl-name lst)) epools))
669
670             (i-gates       (lookup-def 'i-gates gate-complex-info))
671
672             (has-kinetic?  (or (not (null? (filter (lambda (x) (member (car x) kinetic)) states)))))
673             (has-ode?      (or (not (null? (filter (lambda (x) (not (member (car x) kinetic))) states)))
674                                (not (null? pool-ions))))
675
676             (asgn-eq-defs       (poset->asgn-eq-defs poset sys))
677             (reaction-eq-defs   (poset->reaction-eq-defs poset sys kinetic)) 
678             (rate-eq-defs       (reverse (poset->rate-eq-defs poset sys kinetic method)))
679             (kstate-eq-defs     (poset->kinetic-eq-defs poset sys kinetic))
680             (conserve-eq-defs   (poset->state-conserve-eq-defs poset sys))
681             (state-init-defs    (poset->state-init-defs poset sys))
682
683             )
684
685           (pp indent ,nl (TITLE ,sysname))
686           
687           (pp indent ,nl (NEURON "{"))
688           (let recur ((exports exports))
689             (if (not (null? exports)) 
690                 (begin
691                   (pp indent+ (RANGE ,(slp ", " (map nmodl-name (take exports (min 10 (length exports)))))))
692                   (recur (drop exports (min 10 (length exports)))))))
693
694           (let ((currents (append (map (lambda (gate-complex) (nmodl-name (s+ 'i_ (first gate-complex)))) gate-complexes )
695                                   (map (lambda (i-gate) (nmodl-name (s+ 'i_ (second i-gate)))) i-gates ))))
696             (if (not (null? currents)) (pp indent+ (RANGE ,(slp ", " currents)))))
697
698
699           (for-each (lambda (x)
700                       (case (first x)
701                         ((non-specific) 
702                          (pp indent+ (RANGE ,(third x))
703                              (NONSPECIFIC_CURRENT ,(second x))))
704                         (else
705                          (cond ((fourth x)
706                                 (pp indent+ (RANGE ,(second x))
707                                     (USEION ,(first x) READ ,(third x) WRITE ,(second x))))
708                                (else (pp indent+ (RANGE ,(second x))))))))
709                     (delete-duplicates perm-ions (lambda (x y) (eq? (car x) (car y)))))
710
711           
712           (if (null? acc-ions)
713               (for-each (lambda (pool-ion)
714                           (pp indent+ (RANGE ,(slp ", " (list (second pool-ion) (third pool-ion))))
715                               (USEION ,(fourth pool-ion) 
716                                       READ  ,(slp ", " (list (second pool-ion)))
717                                       WRITE ,(slp ", " (list (third pool-ion ))))))
718                         pool-ions)
719               (for-each (lambda (acc-ion)
720                           (let ((pool-ion (assoc (first acc-ion) pool-ions)))
721                             (if pool-ion
722                                 (pp indent+ (RANGE ,(second acc-ion))
723                                     (USEION ,(first acc-ion) 
724                                             READ  ,(slp ", " (list (third acc-ion) (fourth acc-ion) (second pool-ion)))
725                                             WRITE ,(slp ", " (list (second acc-ion) (third pool-ion )))))
726                                 (pp indent+ (RANGE ,(second acc-ion))
727                                     (USEION ,(first acc-ion) 
728                                             READ ,(slp ", "  (list (third acc-ion) (fourth acc-ion) ))
729                                             WRITE ,(second acc-ion))))))
730                         (delete-duplicates acc-ions (lambda (x y) (eq? (car x) (car y))))))
731
732           (let* ((const-names   (map first consts))
733                  (is-const?     (lambda (x) (member x const-names)))
734                  (range-consts  (delete-duplicates 
735                                  (fold (lambda (def ax) 
736                                          (let* ((rhs   (second def))
737                                                 (vars  (rhsvars rhs)))
738                                            (append (filter is-const? vars) ax)))
739                                        (list) asgn-eq-defs ))))
740             (if (not (null? range-consts)) (pp indent+ (RANGE ,(slp ", " range-consts)))))
741           
742           
743           (pp indent "}")
744           
745           (let* ((define-fn (make-define-fn table? min-v max-v table-with depend)))
746             (for-each (lambda (fndef) 
747                         (if (not (member (car fndef) builtin-fns))
748                             (apply define-fn (cons indent fndef)))) 
749                       defuns))
750
751           (let* ((parameter-defs 
752                   (filter-map
753                    (lambda (nv)
754                      (and (not (member (first nv) nmodl-builtin-consts))
755                           (let ((v1 (canonicalize-expr/NMODL (second nv))))
756                             (list (first nv) v1))))
757                    consts))
758                  (parameter-locals  (find-locals (map second parameter-defs)))
759                  (state-defs 
760                   (append
761                    (map (lambda (st)
762                           (if (pair? st) (nmodl-state-name (first st) (second st)) 
763                               (nmodl-name st)))
764                         states)
765                    (map nmodl-name reactions)))
766                  (assigned-defs
767                   (filter-map
768                    (lambda (x) 
769                      (let ((x1 (nmodl-name x)))
770                        (and (not (or (member x1 state-defs) (assoc x1 parameter-defs)))
771                             x1)))
772                    (delete-duplicates
773                     (append asgns
774                             (map first imports) 
775                             (map second perm-ions) (map third perm-ions)
776                             (map second acc-ions)  (map fourth acc-ions)
777                             (map second pool-ions) (map third pool-ions)
778                             (map (lambda (gate-complex) (nmodl-name (s+ 'i_ (first gate-complex)))) gate-complexes )
779                             (map (lambda (i-gate) (nmodl-name (s+ 'i_ (second i-gate)))) i-gates )
780                             )))))
781             
782             (pp indent ,nl (PARAMETER "{"))
783             (if (not (null? parameter-locals)) (pp indent+ (LOCAL ,(slp ", " parameter-locals))))
784             (for-each (lambda (def)
785                         (let ((n (nmodl-name (first def))) (b (second def)))
786                           (pp indent+ ,(expr->string/NMODL b n)))) parameter-defs)
787             (case method  ((expeuler)  (pp indent+ dt)))
788             (pp indent "}")
789             (pp indent ,nl (STATE "{"))
790             (for-each (lambda (x) (pp indent+ ,x)) state-defs)
791             (pp indent "}")
792             (pp indent ,nl (ASSIGNED "{"))
793             (for-each (lambda (x) (pp indent+ ,x)) assigned-defs)
794             (pp indent "}"))
795
796           (if (not (null? asgns))
797               (begin
798                 (pp indent ,nl (PROCEDURE asgns () "{"))
799                 (let ((locals    (find-locals (map second asgn-eq-defs))) )
800                   (if (not (null? locals)) (pp indent+ (LOCAL ,(slp ", " locals)))))
801#|
802                 This seems to cause a segmentation fault in nrnoc:
803
804                 (if (and table? min-v max-v table-with)
805                     (pp indent+ (TABLE ,(slp ", " (map first asgn-eq-defs))
806                                  ,@(if depend `(DEPEND ,depend) `(""))
807                                  FROM ,min-v TO ,max-v WITH ,table-with)))
808|#
809
810                 (for-each (lambda (def)
811                             (let ((n (nmodl-name (first def)) )
812                                   (b (second def)))
813                               (pp indent+ ,(expr->string/NMODL b n)))) asgn-eq-defs)
814                 (pp indent "}")))
815             
816           (if (not (null? reactions))
817               (begin
818                 (pp indent ,nl (PROCEDURE reactions () "{"))
819                 (let ((locals    (find-locals (map second reaction-eq-defs))) )
820                   (if (not (null? locals)) (pp indent+ (LOCAL ,(slp ", " locals))))
821                   (for-each (lambda (def)
822                               (let ((n (nmodl-name (first def))) (b (second def)))
823                                 (pp indent+ ,(expr->string/NMODL b n)))) 
824                             reaction-eq-defs))
825                 (pp indent "}")))
826
827           (if (not (null? pool-ions))
828               (begin
829                 (pp indent ,nl (PROCEDURE pools () "{"))
830                 (for-each (lambda (pool-ion)
831                             (pp indent+ (,(third pool-ion) = ,(first pool-ion))))
832                           pool-ions)
833                 (pp indent "}")))
834
835           (pp indent ,nl (BREAKPOINT "{"))
836           (let* ((i-eqs (filter-map
837                          (lambda (gate-complex) 
838
839                            (let* ((label             (first gate-complex))
840                                   (n                 (second gate-complex))
841                                   (subcomps          ((dis 'component-subcomps) sys n))
842                                   (acc               (lookup-def 'accumulating-substance subcomps))
843                                   (perm              (lookup-def 'permeating-ion subcomps))
844                                   (permqs            (and perm ((dis 'component-exports) sys (cid perm))))
845                                   (pore              (lookup-def 'pore subcomps))
846                                   (permeability      (lookup-def 'permeability subcomps))
847                                   (gate              (lookup-def 'gate subcomps))
848                                   (sts               (and gate ((dis 'component-exports) sys (cid gate)))))
849
850                              (if (null? permqs)
851                                  (nemo:error 'nemo:nmodl-translator ": ion channel definition " label
852                                              "permeating-ion component lacks exported quantities"))
853                              (if (null? sts)
854                                  (nemo:error 'nemo:nmodl-translator ": ion channel definition " label
855                                              "gate component lacks exported quantities"))
856
857                              (if (not (or pore permeability))
858                                  (nemo:error 'nemo:nmodl-translator ": ion channel definition " label
859                                              "lacks any pore or permeability components"))
860 
861                              (cond ((and perm permeability gate)
862                                     (let* ((i     (nmodl-name (s+ 'i (cn perm))))
863                                            (pmax  (car ((dis 'component-exports) sys (cid permeability))))
864                                            (pwrs  (map (lambda (n) (rate/reaction-power sys n)) sts))
865                                            (sptms (map (lambda (st pwr) (if pwr `(pow ,st ,pwr) st)) sts pwrs))
866                                            (gion  `(* ,pmax ,@sptms)))
867                                         (list i #f gion (nmodl-name (s+ 'i_ label) ))))
868
869                                    ((and perm pore gate)
870
871                                     (case (cn perm)
872                                       ((non-specific)
873                                        (let* ((i     (nmodl-name 'i))
874                                               (e     (car permqs))
875                                               (gmax  (car ((dis 'component-exports) sys (cid pore))))
876                                               (pwrs  (map (lambda (n) (rate/reaction-power sys n)) sts))
877                                               (sptms (map (lambda (st pwr) (if pwr `(pow ,st ,pwr) st)) sts pwrs))
878                                               (gion  `(* ,gmax ,@sptms)))
879                                          (list i e gion (nmodl-name (s+ 'i_ label) ))))
880
881                                       (else
882                                        (let* ((i     (nmodl-name (s+ 'i (cn perm))))
883                                               (e     (nmodl-name (car permqs)))
884                                               (gmax  (car ((dis 'component-exports) sys (cid pore))))
885                                               (pwrs  (map (lambda (n) (rate/reaction-power sys n)) sts))
886                                               (sptms (map (lambda (st pwr) (if pwr `(pow ,st ,pwr) st)) sts pwrs))
887                                               (gion  `(* ,gmax ,@sptms)))
888
889                                          (list i e gion (nmodl-name (s+ 'i_ label)))))))
890
891                                    ((and perm pore)
892                                     (case (cn perm)
893                                       ((non-specific)
894                                        (let* ((i     (nmodl-name 'i))
895                                               (e     (car permqs))
896                                               (gmax  (car ((dis 'component-exports) sys (cid pore)))))
897                                          (list i e gmax (nmodl-name (s+ 'i_ label)))))
898                                       (else
899                                        (nemo:error 'nemo:nmodl-translator ": ion channel definition " label
900                                                    (s+ "(" n ")")
901                                                    "lacks gate component"))))
902                                   
903                                    ((and acc pore gate)
904                                     (let* ((i     (nmodl-name (s+ 'i (cn acc))))
905                                            (gmax  (car ((dis 'component-exports) sys (cid pore))))
906                                            (pwrs  (map (lambda (n) (rate/reaction-power sys n)) sts))
907                                            (sptms (map (lambda (st pwr) (if pwr `(pow ,st ,pwr) st)) sts pwrs))
908                                            (gion  `(* ,gmax ,@sptms)))
909                                       (list i #f gion (nmodl-name (s+ 'i_ label) ))))
910                                   
911                                    (else (nemo:error 'nemo:nmodl-translator ": invalid ion channel definition " 
912                                                      label))
913                                    )))
914                          gate-complexes))
915
916                  (i-eqs  (fold  (lambda (i-gate ax) 
917                                   (let ((i-gate-var (first i-gate)))
918                                     (cons (list (nmodl-name 'i) #f i-gate-var (s+ 'i_ (second i-gate)) ) ax)))
919                                 i-eqs i-gates))
920
921                  (i-bkts (bucket-partition (lambda (x y) (eq? (car x) (car y))) i-eqs))
922
923                  (i-eqs  (fold (lambda (b ax) 
924                                  (match b 
925                                         ((and ps ((i e gion ii) . rst))
926                                          (let loop ((ps ps) (summands (list)) (eqs (list)))
927                                            (if (null? ps)
928
929                                                (let* ((sum0  (sum summands))
930                                                       (sum1  (rhsexpr/NMODL sum0))
931                                                       (sum2  (canonicalize-expr/NMODL sum1)))
932                                                  (append eqs (list (list i sum2)) ax))
933
934                                                (match-let (((i e gion ii) (car ps)))
935
936                                                  (loop (cdr ps) 
937                                                        (cons ii summands) 
938                                                        (let* ((expr0 (rhsexpr/NMODL (if e `(* ,gion (- v ,e)) gion)))
939                                                               (expr1 (canonicalize-expr/NMODL expr0)))
940                                                          (cons (list ii expr1) eqs)))))))
941                                           
942                                         ((i e gion ii)
943                                          (let* ((expr0  (rhsexpr/NMODL (if e `(* ,gion (- v ,e)) gion)))
944                                                 (expr1  (canonicalize-expr/NMODL expr0)))
945                                            (cons (list i expr1) ax)))
946                                         
947                                         (else ax)))
948                                (list) i-bkts))
949                 
950                  (locals (find-locals (map second i-eqs))))
951
952             (if (not (null? locals))    (pp indent+ (LOCAL ,(slp ", " locals))))
953             (if has-ode?
954                 (case method
955                   ((#f expeuler)  (pp indent+ (SOLVE states)))
956                   (else           (pp indent+ (SOLVE states METHOD ,method)))))
957             (if has-kinetic?   (pp indent+  (SOLVE kstates METHOD sparse)))
958             (if (not (null? reactions))   (pp indent+ (reactions ())))
959             (if (not (null? pool-ions)) (pp indent+ (pools ())))
960             (for-each (lambda (p) (pp indent+ ,(expr->string/NMODL (second p) (first p)))) i-eqs)
961             (pp indent "}"))
962           
963           (if has-ode?
964               (let ((locals   (find-locals (map second rate-eq-defs))))
965                 (case method
966                   ((expeuler) (pp indent ,nl (PROCEDURE states () "{")))
967                   (else       (pp indent ,nl (DERIVATIVE states "{"))))
968                 (if (not (null? locals)) (pp indent+ (LOCAL ,(slp ", " locals))))
969                 (if (not (null? asgns))     (pp indent+ (asgns ())))
970                 (let ((prime (case method 
971                                ((expeuler) identity)
972                                (else  (lambda (x) (s+ x "'"))))))
973                   (for-each (lambda (def)
974                               (let ((n (prime (first def)))
975                                     (b (second def)))
976                                 (pp indent+ ,(expr->string/NMODL b n))))
977                             rate-eq-defs))
978                 (pp indent "}")))
979
980           (if has-kinetic?
981               (begin
982                 (pp indent ,nl (KINETIC kstates "{"))
983                 (let* ((exprs             (map second kstate-eq-defs))
984                        (locals            (concatenate 
985                                            (find-locals
986                                             (append (map (lambda (x) (map fourth x)) exprs)
987                                                     (map (lambda (x) (map fifth x)) exprs))))))
988                   (if (not (null? locals)) (pp indent+ (LOCAL ,(slp ", " locals))))
989                   (if (not (null? asgns))     (pp indent+ (asgns ())))
990                   (for-each
991                    (lambda (def)
992                      (let* ((n             (first def))
993                             (eqs           (second def))
994                             (conserve-eqs  (lookup-def (nmodl-name n) conserve-eq-defs)))
995                        (for-each
996                         (lambda (eq)
997                           (match eq
998                                  (('-> s0 s1 rexpr) 
999                                   (pp indent+ (~ ,s0 -> ,s1 (,(expr->string/NMODL rexpr)))))
1000                                  (('<-> s0 s1 rexpr1 rexpr2) 
1001                                   (pp indent+ (~ ,s0 <-> ,s1 (,(expr->string/NMODL rexpr1) #\,
1002                                                               ,(expr->string/NMODL rexpr2) 
1003                                                               ))))
1004                                  ))
1005                           eqs)
1006                        (if conserve-eqs
1007                            (for-each (lambda (eq) 
1008                                        (let ((val  (first eq))
1009                                              (expr (third eq)))
1010                                          (pp indent+ ,(conserve-conseq->string/NMODL expr val))))
1011                                    conserve-eqs))
1012                        ))
1013                    kstate-eq-defs))
1014                 (pp indent "}")))
1015           
1016           
1017           (let ((locals (concatenate (find-locals (map second state-init-defs)))) )
1018               (pp indent ,nl (INITIAL "{"))
1019               (if (not (null? locals)) (pp indent+ (LOCAL ,(slp ", " locals))))
1020               (if (not (null? asgns))  (pp indent+ (asgns ())))
1021               (for-each (lambda (def)
1022                           (let ((n (first def)) (b (second def)))
1023                             (pp indent+ ,(expr->string/NMODL b n)))) state-init-defs)
1024
1025               (if has-kinetic?
1026                   (pp indent+ (SOLVE kstates STEADYSTATE sparse)))
1027               
1028               (pp indent "}")
1029
1030               (pp indent ,nl (PROCEDURE print_state () "{"))
1031
1032               (let ((lst (sort (map (compose ->string first) rate-eq-defs) string<?)))
1033                 (for-each (lambda (x)
1034                             (pp indent+ (printf (,(s+ #\" x " = %g\\n"  #\") ", " ,x ))))
1035                           lst))
1036
1037               (pp indent "}")
1038
1039               ))))))
1040)
Note: See TracBrowser for help on using the repository browser.