source: project/release/4/nemo/tags/5.1/nemo-nest.scm @ 25871

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

nemo release 5.1

File size: 26.7 KB
Line 
1;;       
2;;
3;; An extension for translating NEMO models to NEST code.
4;;
5;; Copyright 2011-2012 Ivan Raikov and the Okinawa Institute of Science and Technology
6;;
7;; This program is free software: you can redistribute it and/or
8;; modify it under the terms of the GNU General Public License as
9;; published by the Free Software Foundation, either version 3 of the
10;; License, or (at your option) any later version.
11;;
12;; This program is distributed in the hope that it will be useful, but
13;; WITHOUT ANY WARRANTY; without even the implied warranty of
14;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15;; General Public License for more details.
16;;
17;; A full copy of the GPL license can be found at
18;; <http://www.gnu.org/licenses/>.
19;;
20
21(module nemo-nest
22
23        (nemo:nest-translator)
24
25        (import scheme chicken utils data-structures lolevel ports srfi-1 srfi-13)
26       
27        (require-extension lolevel matchable strictly-pretty environments
28                           varsubst datatype nemo-core nemo-utils
29                           nemo-gate-complex)
30
31(define nest-builtin-consts
32  `())
33
34(define C++-ops
35  `(+ - * / > < <= >= = ^))
36
37(define builtin-fns
38  `(+ - * / pow neg abs atan asin acos sin cos exp ln
39      sqrt tan cosh sinh tanh hypot gamma lgamma log10 log2 log1p ldexp cube
40      > < <= >= = and or round ceiling floor max min
41      fpvector-ref))
42
43(define (nest-name s)
44  (let ((cs (string->list (->string s))))
45    (let loop ((lst (list)) (cs cs))
46      (if (null? cs) (string->symbol (list->string (reverse lst)))
47          (let* ((c (car cs))
48                 (c1 (cond ((or (char-alphabetic? c) (char-numeric? c) (char=? c #\_)) c)
49                           (else #\_))))
50            (loop (cons c1 lst) (cdr cs)))))))
51           
52(define (rhsexpr/C++ expr)
53  (match expr 
54         (('if . es)  `(if . ,(map (lambda (x) (rhsexpr/C++ x)) es)))
55         (('pow x y)  (if (and (integer? y)  (positive? y))
56                          (if (> y 1)  (let ((tmp (gensym "x")))
57                                         `(let ((,tmp ,x)) (* . ,(list-tabulate (inexact->exact y) (lambda (i) tmp)))))
58                              x)
59                            expr))
60         ((s . es)    (if (symbol? s)  (cons (if (member s builtin-fns) s (nest-name s)) (map (lambda (x) (rhsexpr/C++ x)) es)) expr))
61         (id          (if (symbol? id) (nest-name id) id))))
62
63(define (nest-state-name n s)
64  (nest-name (s+ n s)))
65
66(define-syntax pp
67  (syntax-rules ()
68    ((pp indent val ...) (ppf indent (quasiquote val) ...))))
69
70
71(define group/C++   (doc:block 2 (doc:text "(") (doc:text ")")))
72(define block/C++   (doc:block 2 (doc:empty) (doc:empty)))
73(define (stmt/C++ x) 
74  (match x
75         (($ doc 'DocCons _ ($ doc 'DocText ";")) x)
76         (else  (doc:cons x (doc:text ";")))))
77
78
79(define (ifthen/C++ c e1 e2)
80  (doc:nest 
81   2 (doc:connect
82      (doc:connect (doc:group (doc:connect (doc:text "if") c))
83                   (doc:connect (doc:nest 2 e1)
84                                (doc:nest 2 (doc:connect (doc:text "else") e2))))
85      (doc:text "end"))))
86
87(define (letblk/C++ e1 e2)
88  (cond ((equal? e1 (doc:empty)) (doc:group (doc:nest 2 e2)))
89        ((equal? e2 (doc:empty)) (doc:group (doc:nest 2 e1)))
90        (else (doc:connect (doc:group (doc:nest 2 (stmt/C++ e1)))
91                           (doc:group (doc:nest 2 e2))))))
92       
93
94(define (format-op/C++ indent op args)
95  (let ((op1 (doc:text (->string op))))
96    (if (null? args) op1
97        (match args
98               ((x)      (doc:concat (list op1 x)))
99               ((x y)    (doc:concat (intersperse (list x op1 y) (doc:space))))
100               ((x y z)  (doc:concat (intersperse (list x op1 y op1 z) (doc:space))))
101               (lst      (let* ((n   (length lst))
102                                (n/2 (inexact->exact (round (/ n 2)))))
103                           (doc:concat 
104                            (intersperse 
105                             (list (format-op/C++ indent op (take lst n/2 )) op1 
106                                   (format-op/C++ indent op (drop lst n/2 )))
107                             (doc:space)))))))))
108
109(define (format-fncall/C++ indent op args)
110  (let ((op1 (doc:text (->string op))))
111    (doc:cons op1 (group/C++ ((doc:list indent identity (lambda () (doc:text ", "))) args)))))
112
113(define (name-normalize expr)
114  (match expr 
115         (('if c t e)  `(if ,(name-normalize c) ,(name-normalize t) ,(name-normalize e)))
116         (('let bs e)
117          `(let ,(map (lambda (b) `(,(car b) ,(name-normalize (cadr b)))) bs) ,(name-normalize e)))
118         ((f . es) 
119          (cons f (map name-normalize es)))
120         ((? symbol? ) (nest-name expr))
121         ((? atom? ) expr)))
122
123
124(define (canonicalize-expr/C++ expr)
125  (let ((subst-convert  (subst-driver (lambda (x) (and (symbol? x) x)) nemo:binding? identity nemo:bind nemo:subst-term)))
126    (let* ((expr1 (if-convert expr))
127           (expr2 (subst-convert expr1 subst-empty))
128           (expr3 (let-lift expr2))
129           (expr4 (name-normalize expr3)))
130      expr4)))
131
132
133(define (format-expr/C++ indent expr . rest) 
134  (let-optionals rest ((rv #f))
135   (let ((indent+ (+ 2 indent)))
136    (match expr
137       (('let bindings body)
138        (letblk/C++
139         (fold-right 
140           (lambda (x ax)
141             (letblk/C++
142              (match (second x)
143                     (('if c t e)
144                      (ifthen/C++
145                       (group/C++ (format-expr/C++ indent c))
146                       (block/C++ (format-expr/C++ indent t (first x)))
147                       (block/C++ (format-expr/C++ indent e (first x)))))
148                     (else
149                      (stmt/C++
150                       (format-op/C++ indent+ " = "
151                                         (list (format-expr/C++ indent (first x) )
152                                               (format-expr/C++ indent (second x)))))))
153              ax))
154           (doc:empty) bindings)
155         (match body
156                (('let _ _) (format-expr/C++ indent body rv))
157                (else
158                 (let ((body1 (doc:nest indent (format-expr/C++ indent body))))
159                   (if rv (stmt/C++ (format-op/C++ indent " = " (list (format-expr/C++ indent+ rv ) body1)))
160                       body1))))))
161       
162       (('if . rest) (error 'format-expr/C++ "invalid if statement " expr))
163
164       ((op . rest) 
165       (let ((op (case op ((pow)  '^)  ((ln) 'log) (else op))))
166         (let ((fe
167                (if (member op C++-ops)
168                    (let ((mdiv?  (any (lambda (x) (match x (('* . _) #t) (('/ . _) #t) (else #f))) rest))
169                          (mul?   (any (lambda (x) (match x (('* . _) #t) (else #f))) rest))
170                          (plmin? (any (lambda (x) (match x (('+ . _) #t) (('- . _) #t) (else #f))) rest)))
171                      (case op
172                        ((/) 
173                         (format-op/C++ indent op 
174                                          (map (lambda (x) 
175                                                 (let ((fx (format-expr/C++ indent+ x)))
176                                                   (if (or (symbol? x) (number? x)) fx
177                                                       (if (or mul? plmin?) (group/C++ fx) fx)))) rest)))
178                        ((*) 
179                         (format-op/C++ indent op 
180                                          (map (lambda (x) 
181                                                 (let ((fx (format-expr/C++ indent+ x)))
182                                                   (if (or (symbol? x) (number? x)) fx
183                                                       (if plmin? (group/C++ fx) fx)))) rest)))
184                       
185                        ((^) 
186                         (format-op/C++ indent op 
187                                          (map (lambda (x) 
188                                                 (let ((fx (format-expr/C++ indent+ x)))
189                                                   (if (or (symbol? x)  (number? x)) fx
190                                                       (if (or mdiv? plmin?) (group/C++ fx) fx)))) rest)))
191                       
192                        (else
193                         (format-op/C++ indent op 
194                                          (map (lambda (x) 
195                                                 (let ((fx (format-expr/C++ indent+ x))) fx)) rest)))))
196                   
197                    (let ((op (case op ((neg) '-) (else op))))
198                      (format-fncall/C++ indent op (map (lambda (x) (format-expr/C++ indent+ x)) rest))))))
199           (if rv 
200               (stmt/C++ (format-op/C++ indent " = " (list (format-expr/C++ indent+ rv ) fe)))
201               fe))))
202     
203      (else  (let ((fe (doc:text (->string expr))))
204               (if rv 
205                   (stmt/C++ (format-op/C++ indent " = " (list (format-expr/C++ indent+ rv ) fe)))
206                   fe)))))))
207               
208         
209(define (expr->string/C++ x . rest)
210  (let-optionals rest ((rv #f) (width 72))
211    (sdoc->string (doc:format width (format-expr/C++ 2 x rv)))))
212
213
214
215(define (state-init n init)
216  (let* ((init  (rhsexpr/C++ init))
217         (init1 (canonicalize-expr/C++ init)))
218    (list (nest-name n) init1)))
219
220
221(define (asgn-eq n rhs)
222  (let* ((fbody   (rhsexpr/C++ rhs))
223         (fbody1  (canonicalize-expr/C++ fbody)))
224    (list (nest-name n) fbody1)))
225
226
227(define (reaction-eq n open transitions conserve)
228  (match-let (((g cnode node-subs)  (transitions-graph n open transitions conserve nest-state-name)))
229    (let ((nodes ((g 'nodes))))
230      (list (nest-name n) (nest-state-name n open)))))
231
232
233(define (reaction-transition-eqs n initial open transitions conserve power)
234  (match-let (((g cnode node-subs)  (transitions-graph n open transitions conserve nest-state-name)))
235     (let* ((out-edges  (g 'out-edges))
236            (in-edges   (g 'in-edges))
237            (nodes      ((g 'nodes))))
238       ;; generate differential equations for each state in the transitions system
239      (let ((eqs    (fold (lambda (s ax) 
240                            (if (= (first cnode) (first s) ) ax
241                                (let* ((out   (out-edges (first s)))
242                                       (in    (in-edges (first s)))
243                                       (open? (eq? (second s) open))
244                                       (name  (nest-name (lookup-def (second s) node-subs))))
245                                  (let* ((rhs1  (cond ((and (not (null? out)) (not (null? in)))
246                                                       `(+ (neg ,(sum (map third out)))
247                                                           ,(sum (map third in))))
248                                                      ((and (not (null? out)) (null? in))
249                                                       `(neg ,(sum (map third out))))
250                                                      ((and (null? out) (not (null? in)))
251                                                       (sum (map third in)))))
252                                         (fbody0 (rhsexpr/C++ rhs1)))
253                                   
254                                    (cons (list name (canonicalize-expr/C++ fbody0)) ax)
255                                    ))))
256                          (list) nodes)))
257        eqs))))
258           
259       
260
261(define (poset->asgn-eq-defs poset sys)
262  (fold-right
263   (lambda (lst ax)
264     (fold  (lambda (x ax) 
265              (match-let (((i . n)  x))
266                         (let ((en (environment-ref sys n)))
267                           (if (nemo:quantity? en)
268                               (cases nemo:quantity en
269                                      (ASGN  (name value rhs) (cons (asgn-eq name rhs) ax))
270                                      (else  ax))
271                               ax))))
272            ax lst))
273   (list) poset))
274
275(define (poset->rate-eq-defs poset sys)
276  (fold-right
277   (lambda (lst ax)
278     (fold  (lambda (x ax) 
279              (match-let (((i . n)  x))
280                         (let ((en (environment-ref sys n)))
281                           (if (nemo:quantity? en)
282                               (cases nemo:quantity en
283
284                                      (REACTION  (name initial open transitions conserve power) 
285                                                 (append (reaction-transition-eqs name initial open transitions 
286                                                                                  conserve power) ax))
287                                     
288                                      (RATE (name initial rhs power)
289                                            (let ((fbody0  (rhsexpr/C++ rhs))
290                                                  (dy      (nest-name name) ))
291                                              (cons (list dy (canonicalize-expr/C++ fbody0)) ax)))
292
293                                      (else  ax))
294                               ax))))
295            ax lst))
296   (list) poset))
297
298
299(define (poset->reaction-eq-defs poset sys)
300  (fold-right
301   (lambda (lst ax)
302     (fold  (lambda (x ax) 
303              (match-let (((i . n)  x))
304                         (let ((en (environment-ref sys n)))
305                           (if (nemo:quantity? en)
306                               (cases nemo:quantity en
307                                      (REACTION  (name initial open transitions conserve power) 
308                                                 (cons (reaction-eq name open transitions conserve) ax))
309                                      (else  ax))
310                               ax))))
311            ax lst))
312   (list) poset))
313
314
315(define (poset->init-defs poset sys)
316  (fold-right
317   (lambda (lst ax)
318     (fold-right
319      (lambda (x ax) 
320              (match-let (((i . n)  x))
321                         (let ((en (environment-ref sys n)))
322                           (if (nemo:quantity? en)
323                               (cases nemo:quantity en
324                                      (REACTION  (name initial open transitions conserve power)
325                                                 (if (nemo:rhs? initial)
326                                                     (cons* (state-init name initial) 
327                                                            (state-init (nest-state-name name open) name) ax) 
328                                                     ax))
329
330                                      (RATE  (name initial rhs power)
331                                             (if (and initial (nemo:rhs? initial))
332                                                 (cons (state-init name initial) ax) 
333                                                 ax))
334
335                                      (else  ax))
336                               ax))))
337            ax lst))
338   (list) poset))
339
340
341(define (poset->state-conserve-eq-defs poset sys)
342  (fold-right
343   (lambda (lst ax)
344     (fold  (lambda (x ax) 
345              (match-let (((i . n)  x))
346                         (let ((en (environment-ref sys n)))
347                           (if (nemo:quantity? en)
348                               (cases nemo:quantity en
349                                      (REACTION (name initial open transitions conserve power)
350                                                (if (and (list? conserve) (every nemo:conseq? conserve))
351                                                    (cons (state-conseqs (nest-name name) transitions conserve
352                                                                        nest-state-name) ax) 
353                                                    ax))
354                                      (else  ax))
355                               ax))))
356            ax lst))
357   (list) poset))
358
359
360(define (reaction-power sys n)
361  (let ((en (environment-ref sys n)))
362    (if (nemo:quantity? en)
363        (cases nemo:quantity en
364               (REACTION  (name initial open transitions conserve power)  power)
365               (else  #f))  #f)))
366
367
368(define (bucket-partition p lst)
369  (let loop ((lst lst) (ax (list)))
370    (if (null? lst) ax
371        (let ((x (car lst)))
372          (let bkt-loop ((old-bkts ax) (new-bkts (list)))
373            (if (null? old-bkts) (loop (cdr lst) (cons (list x) new-bkts))
374                (if (p x (caar old-bkts ))
375                    (loop (cdr lst) (append (cdr old-bkts) (cons (cons x (car old-bkts)) new-bkts)))
376                    (bkt-loop (cdr old-bkts) (cons (car old-bkts) new-bkts)))))))))
377
378
379(define (make-define-fn )
380  (lambda (indent n proc)
381    (let ((lst (procedure-data proc))
382          (indent+ (+ 2 indent)))
383      (let ((rt       (lookup-def 'rt lst))
384            (formals  (lookup-def 'formals lst))
385            (vars     (lookup-def 'vars lst))
386            (body     (lookup-def 'body lst))
387            (rv       (gensym 'rv)))
388        (pp indent ,nl (,rt ,(nest-name n) (,(slp ", " vars)) "{" ))
389        (let* ((body0 (rhsexpr/C++ body))
390               (body1 (canonicalize-expr/C++ body0))
391               (lbs   (enum-bnds body1 (list))))
392          (if (not (null? lbs)) (pp indent+ (double ,(slp ", " lbs))))
393          (pp indent+ ,(expr->string/C++ body1 (nest-name rv)))
394          (pp indent+ ,(s+ "return " rv ";"))
395          )
396        (pp indent "}"))) 
397    ))
398
399
400(define (output-dy sysname parameters state-index-map 
401                   rate-eq-defs reaction-eq-defs asgn-eq-defs
402                   pool-ions mcap i-eqs v-eq indent indent+)
403
404  (pp indent ,nl (int ,(s+ sysname "_dynamics")
405                      (,(slp ", " `("double t"
406                                    "const double y[]"
407                                    "double f[]"
408                                    "void* pnode"
409                                    )))
410                      #\{
411                      ))
412
413  (pp indent+ ,nl
414      ("// S is shorthand for the type that describes the model state ")
415      (typedef ,(s+ sysname "::state_t") "S;")
416
417      ,nl
418
419      ("// cast the node ptr to an object of the proper type")
420      ("assert(pnode);")
421      ("const " ,sysname "& node =  *(reinterpret_cast<" ,sysname "*>(pnode));")
422
423      ,nl
424     
425      ("// y[] must be the state vector supplied by the integrator, ")
426      ("// not the state vector in the node, node.S_.y[]. ")
427      ,nl
428      )
429
430  ;; TODO parameters
431           
432  (for-each (lambda (def)
433              (let ((n (first def)) )
434                (pp indent+
435                    ,(s+ "// rate equation " n)
436                    ,(expr->string/C++ 
437                      (s+ "y[" (lookup-def n state-index-map) "]") n))))
438            rate-eq-defs)
439
440  (let* ((eqs 
441          (append
442           
443           asgn-eq-defs
444           reaction-eq-defs
445                 
446           (map (lambda (pool-ion)
447                  (let ((n (third pool-ion))
448                        (b (first pool-ion)))
449                    (list n b)))
450                pool-ions)))
451         
452         (eq-dag 
453          (map (lambda (def)
454                 (cons (first def) (enum-freevars (second def) '() '())))
455               eqs))
456         (eq-order
457          (reverse
458           (topological-sort eq-dag 
459                             (lambda (x y) (string=? (->string x) (->string y)))))))
460          (for-each (lambda (n)
461                      (let ((b  (lookup-def n eqs)))
462                        (if b (pp indent+ 
463                                  ,(s+ "// equation for " n)
464                                  ,(expr->string/C++ b (nest-name n))))))
465                    eq-order))
466 
467  (for-each (lambda (def)
468              (let ((n (s+ "f[" (lookup-def (first def) state-index-map) "]")) 
469                    (b (second def)))
470                (pp indent+ ,(s+ "// state " (first def))
471                    ,(expr->string/C++ b n))))
472            rate-eq-defs)
473 
474  (for-each (lambda (def) 
475              (pp indent+ ,(expr->string/C++ (second def) (first def)))) 
476            i-eqs)
477
478  (if v-eq
479      (let ((vi (lookup-def 'v state-index-map)))
480        (if vi
481            (pp indent+ ,(expr->string/C++ (second v-eq) (s+ "f[" vi "]"))))))
482 
483  (pp indent+ ,nl ("return GSL_SUCCESS;"))
484  (pp indent  #\})
485)
486
487#|
488(define (output-init sysname globals state-index-map steady-state-index-map
489                     const-defs asgn-eq-defs init-eq-defs rate-eq-defs
490                     reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
491
492    (pp indent ,nl (function y0 = ,(s+ sysname "_init") (Vinit)))
493   
494    (if (not (null? globals)) (pp indent+ (global ,(slp " " globals))))
495   
496    (pp indent+ ,(expr->string/C++ `(zeros ,(length state-index-map) 1) 'y0))
497   
498    (if (member 'v globals)
499        (let ((vi (lookup-def 'v state-index-map)))
500          (pp indent+ ,(expr->string/C++ `Vinit 'v))
501          (if vi (pp indent+ ,(expr->string/C++ 'v (s+ "y0(" vi ")") )))))
502
503
504    (let* ((init-eqs
505            (append
506             
507             const-defs
508             asgn-eq-defs
509             init-eq-defs
510             
511             (map (lambda (pool-ion)
512                    (let ((n (third pool-ion))
513                          (b (first pool-ion)))
514                      (list n b)))
515                  pool-ions)))
516
517           (init-dag
518            (map (lambda (def)
519                   (cons (first def) (enum-freevars (second def) '() '())))
520                 init-eqs))
521           (init-order
522            (reverse
523             (topological-sort init-dag
524               (lambda (x y) (string=? (->string x) (->string y)))))))
525
526      (for-each (lambda (n)
527                  (let ((b  (lookup-def n init-eqs)))
528                    (if b (pp indent+ ,(expr->string/C++ b (nest-name n))))))
529                init-order))
530
531    (for-each (lambda (def)
532                (let* ((n  (first def))
533                       (ni (lookup-def n state-index-map)))
534                  (if ni (pp indent+ ,(expr->string/C++ n  (s+ "y0(" ni ")"))))))
535              init-eq-defs)
536
537    (if (not (null? steady-state-index-map))
538        (begin
539          (pp indent+ ,(expr->string/C++ `(zeros ,(length steady-state-index-map) 1) 'xs)
540              ,(expr->string/C++
541                `(fsolve ,(s+ #\' sysname "_steadystate" #\' ) xs)
542                "[ys,fval,info]"))
543         
544         
545          (for-each
546           (lambda (def)
547             (let* ((n (first def))
548                    (si (lookup-def n steady-state-index-map))
549                    (ni (lookup-def n state-index-map)))
550               (if si (begin
551                        (pp indent+ ,(expr->string/C++ (s+ "ys(" si ")") n))
552                        (pp indent+ ,(expr->string/C++ (s+ "ys(" si ")") (s+ "y0(" ni ")")))))))
553           rate-eq-defs)
554
555
556          ))
557
558    (for-each (lambda (def)
559                (let ((n (first def)) (b (second def)))
560                  (if (not (lookup-def n init-eq-defs))
561                      (pp indent+ ,(expr->string/C++ b n)))))
562            reaction-eq-defs)
563
564    (for-each (lambda (def)
565                (pp indent+ ,(expr->string/C++ (second def) (first def))))
566              i-eqs)
567   
568    (pp indent ,nl (end))
569    (pp indent ,nl (}))
570
571    )
572|#
573
574(define (nemo:nest-translator sys . rest)
575
576  (define (cid x)  (second x))
577  (define (cn x)   (first x))
578
579  (let-optionals rest ((filename #f))
580  (match-let ((($ nemo:quantity 'DISPATCH  dis) (environment-ref sys (nemo-intern 'dispatch))))
581    (let ((imports  ((dis 'imports)  sys))
582          (exports  ((dis 'exports)  sys)))
583      (let* ((indent      0)
584             (indent+     (+ 2 indent ))
585             (eval-const  (dis 'eval-const))
586             (sysname     (nest-name ((dis 'sysname) sys)))
587             (prefix      sysname)
588             (filename    (or filename (s+ sysname ".cpp")))
589             (deps*       ((dis 'depgraph*) sys))
590             (consts      ((dis 'consts)  sys))
591             (asgns       ((dis 'asgns)   sys))
592             (states      ((dis 'states)  sys))
593             (reactions   ((dis 'reactions) sys))
594             (defuns      ((dis 'defuns)  sys))
595             (components  ((dis 'components) sys))
596             
597             (g             (match-let (((state-list asgn-list g) ((dis 'depgraph*) sys))) g))
598             (poset         (vector->list ((dis 'depgraph->bfs-dist-poset) g)))
599
600             (const-defs       (filter-map
601                                (lambda (nv)
602                                  (and (not (member (first nv) nest-builtin-consts))
603                                       (let ((v1 (canonicalize-expr/C++ (second nv))))
604                                         (list (nest-name (first nv)) v1))))
605                                consts))
606             
607             (gate-complex-info    (nemo:gate-complex-query sys))
608
609             (gate-complexes       (lookup-def 'gate-complexes gate-complex-info))
610             (perm-ions     (map (match-lambda ((comp i e erev) `(,comp ,(nest-name i) ,(nest-name e) ,erev)))
611                                 (lookup-def 'perm-ions gate-complex-info)))
612             (acc-ions      (map (match-lambda ((comp i in out) `(,comp ,@(map nest-name (list i in out)))))
613                                 (lookup-def 'acc-ions gate-complex-info)))
614             (epools        (lookup-def 'pool-ions gate-complex-info))
615             (pool-ions     (map (lambda (lst) (map nest-name lst)) epools))
616
617             (i-gates       (lookup-def 'i-gates gate-complex-info))
618
619             (capcomp       (any (match-lambda ((name 'membrane-capacitance id) (list name id)) (else #f)) components))
620             (mcap          (and capcomp (car ((dis 'component-exports) sys (cid capcomp)))))
621               
622
623             (i-eqs (filter-map
624                     (lambda (gate-complex) 
625                       
626                       (let* ((label             (first gate-complex))
627                              (n                 (second gate-complex))
628                              (subcomps          ((dis 'component-subcomps) sys n))
629                              (acc               (lookup-def 'accumulating-substance subcomps))
630                              (perm              (lookup-def 'permeating-ion subcomps))
631                              (permqs            (and perm ((dis 'component-exports) sys (cid perm))))
632                              (pore              (lookup-def 'pore subcomps))
633                              (permeability      (lookup-def 'permeability subcomps))
634                              (gate              (lookup-def 'gate subcomps))
635                              (sts               (and gate ((dis 'component-exports) sys (cid gate)))))
636                         
637                         (if (not (or pore permeability))
638                             (nemo:error 'nemo:nest-translator ": ion channel definition " label
639                                         "lacks any pore or permeability components"))
640
641
642                         (cond ((and perm permeability gate)
643                                     (let* ((i     (nest-name (s+ 'i (cn perm))))
644                                            (pmax  (car ((dis 'component-exports) sys (cid permeability))))
645                                            (pwrs  (map (lambda (n) (reaction-power sys n)) sts))
646                                            (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
647                                            (gion  `(* ,pmax ,@sptms)))
648                                         (list i #f gion (nest-name (s+ 'i_ label) ))))
649
650                               ((and perm pore gate)
651                                (case (cn perm)
652                                  ((non-specific)
653                                   (let* ((i     (nest-name 'i))
654                                          (e     (car permqs))
655                                          (gmax  (car ((dis 'component-exports) sys (cid pore))))
656                                          (pwrs  (map (lambda (n) (reaction-power sys n)) sts))
657                                          (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
658                                          (gion  `(* ,gmax ,@sptms)))
659                                     (list i e gion  (nest-name (s+ 'i_ label) ))))
660
661                                  (else
662                                   (let* ((i     (nest-name (s+ 'i (cn perm))))
663                                          (e     (car permqs))
664                                          (gmax  (car ((dis 'component-exports) sys (cid pore))))
665                                          (pwrs  (map (lambda (n) (reaction-power sys n)) sts))
666                                          (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
667                                          (gion  `(* ,gmax ,@sptms)))
668                                     (list i e gion (nest-name (s+ 'i_ label) ))))))
669                               
670                               ((and perm pore)
671                                (case (cn perm)
672                                  ((non-specific)
673                                   (let* ((i     (nest-name 'i))
674                                          (e     (car permqs))
675                                          (gmax  (car ((dis 'component-exports) sys (cid pore)))))
676                                     (list i e gmax (nest-name (s+ 'i_ label) ))))
677                                  (else
678                                   (nemo:error 'nemo:nest-translator ": invalid ion channel definition " label))))
679                               
680                               ((and acc pore gate)
681                                (let* ((i     (nest-name (s+ 'i (cn acc))))
682                                       (gmax  (car ((dis 'component-exports) sys (cid pore))))
683                                       (pwrs  (map (lambda (n) (reaction-power sys n)) sts))
684                                       (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
685                                       (gion  `(* ,gmax ,@sptms)))
686                                  (list i #f gion  (nest-name (s+ 'i_ label) ))))
687                               (else (nemo:error 'nemo:nest-translator ": invalid ion channel definition " label))
688                               )))
689                     gate-complexes))
690
691             (i-names (delete-duplicates (map first i-eqs)))
692               
693             (i-eqs  (fold  (lambda (i-gate ax) 
694                              (let ((i-gate-var (first i-gate)))
695                                (cons (list (nest-name 'i) #f i-gate-var (s+ 'i_ (second i-gate))) ax)))
696                            i-eqs i-gates))
697
698             (i-bkts (bucket-partition (lambda (x y) (eq? (car x) (car y))) i-eqs))
699             
700             (i-eqs  (fold (lambda (b ax) 
701                             (match b 
702                                    ((and ps ((i e gion ii) . rst)) 
703                                     (let loop ((ps ps) (summands (list)) (eqs (list)))
704                                       (if (null? ps)
705                                           
706                                           (let* ((sum0  (sum summands))
707                                                  (sum1  (rhsexpr/C++ sum0))
708                                                  (sum2  (canonicalize-expr/C++ sum1)))
709                                             (append eqs (list (list i sum2)) ax))
710                                           
711                                           (match-let (((i e gion ii) (car ps)))
712                                                      (loop (cdr ps) 
713                                                            (cons ii summands) 
714                                                            (let* ((expr0 (rhsexpr/C++ (if e `(* ,gion (- v ,e)) gion)))
715                                                                   (expr1 (canonicalize-expr/C++ expr0)))
716                                                              (cons (list ii expr1) eqs)))))))
717                                   
718                                    ((i e gion ii)
719                                     (let* ((expr0  (rhsexpr/C++ (if e `(* ,gion (- v ,e)) gion)))
720                                            (expr1  (canonicalize-expr/C++ expr0)))
721                                       (cons (list i expr1) ax)))
722                                   
723                                    (else ax)))
724                           (list) i-bkts))
725
726             (asgn-eq-defs     (poset->asgn-eq-defs poset sys))
727             
728             (rate-eq-defs       (reverse (poset->rate-eq-defs poset sys)))
729             
730             (reaction-eq-defs   (poset->reaction-eq-defs poset sys))
731             
732             (init-eq-defs       (poset->init-defs poset sys))
733             
734             (conserve-eq-defs   (map (lambda (eq) (list 0 `(- ,(second eq) ,(first eq)))) 
735                                      (poset->state-conserve-eq-defs poset sys)))
736             
737             (globals          (map nest-name 
738                                    (delete-duplicates (append
739                                                        exports
740                                                        (map second  perm-ions)
741                                                        (map third   perm-ions)
742                                                        (map second  acc-ions)
743                                                        (map third   acc-ions)
744                                                        (map fourth  acc-ions)
745                                                        (map second  pool-ions)
746                                                        (map third   pool-ions)
747                                                        (map first   imports) 
748                                                        (map first   const-defs)
749                                                        (map first   asgn-eq-defs)
750                                                        (map (lambda (gate-complex) (nest-name (s+ 'i_ (first gate-complex)))) gate-complexes )
751                                                        (map (lambda (i-gate) (nest-name (s+ 'i_ (second i-gate)))) i-gates )
752                                                        ))))
753
754             (parameters (map nest-name (map first const-defs)))
755
756             (v-eq    (if (and mcap (member 'v globals))
757                          (list 'v (rhsexpr/C++ `(/ (neg ,(sum i-names)) ,mcap)))
758                          (list 'v 0.0)))
759             
760             (state-index-map  (let ((acc (fold (lambda (def ax)
761                                                  (let ((st-name (first def)))
762                                                    (list (+ 1 (first ax)) 
763                                                          (cons `(,st-name ,(first ax)) (second ax)))))
764                                                (list 1 (list)) 
765                                                (if (member 'v globals)
766                                                    (cons (list 'v) rate-eq-defs)
767                                                    rate-eq-defs)
768                                                )))
769                                 (second acc)))
770             
771             (steady-state-index-map  (let ((acc (fold (lambda (def ax)
772                                                         (let ((st-name (first def)))
773                                                           (if (not (alist-ref st-name init-eq-defs))
774                                                               (list (+ 1 (first ax)) 
775                                                                     (cons `(,st-name ,(first ax)) (second ax)))
776                                                               ax)))
777                                                       (list 1 (list)) 
778                                                       rate-eq-defs)))
779                                        (second acc)))
780             
781             (dfenv 
782              (map (lambda (x) (let ((n (first x)))
783                                 (list n (nest-name (s+ "d_" n )))))
784                   defuns))
785
786             )
787       
788       
789       
790        (for-each
791         (lambda (a)
792           (let ((acc-ion   (car a)))
793             (if (assoc acc-ion perm-ions)
794                 (nemo:error 'nemo:nest-translator 
795                             ": ion species " acc-ion " cannot be declared as both accumulating and permeating"))))
796         acc-ions)
797
798           (for-each
799            (lambda (p)
800              (let ((pool-ion  (car p)))
801                (if (assoc pool-ion perm-ions)
802                    (nemo:error 'nemo:nest-translator 
803                                ": ion species " pool-ion " cannot be declared as both pool and permeating"))))
804            pool-ions)
805
806           (let ((output  (open-output-file filename)))
807             
808             ;; derivative function
809             (with-output-to-port output
810               (lambda ()
811                 (output-dy sysname globals state-index-map 
812                            rate-eq-defs reaction-eq-defs asgn-eq-defs
813                            pool-ions mcap i-eqs v-eq
814                            indent indent+)
815                 ))
816
817
818#|
819
820             ;; initial values function
821             (with-output-to-port output
822               (lambda ()
823                 (output-init sysname globals state-index-map steady-state-index-map
824                              const-defs asgn-eq-defs init-eq-defs rate-eq-defs
825                              reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
826                 (pp indent ,nl)
827                 ))
828               
829|#
830
831           ;; user-defined functions
832           (let ((define-fn  (make-define-fn)))
833             (for-each (lambda (fndef) 
834                         (if (not (member (car fndef) builtin-fns))
835                             (with-output-to-port output
836                               (lambda ()
837                                 (apply define-fn (cons indent fndef))
838                                 (pp indent ,nl)))
839                             ))
840                       defuns))
841         
842           (close-output-port output))
843
844           )))
845  ))
846
847
848)
Note: See TracBrowser for help on using the repository browser.