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

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

nemo: more uniform code generaton for octave and nmodl

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