source: project/release/3/nemo/trunk/nemo-nmodl.scm @ 13050

Last change on this file since 13050 was 13050, checked in by Ivan Raikov, 11 years ago

Fixes to the matlab backend; better error messages in nmodl backend.

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