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

Last change on this file since 12685 was 12685, checked in by Ivan Raikov, 13 years ago

Bug fixes.

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