source: project/release/4/nemo/trunk/nemo-matlab.scm @ 28033

Last change on this file since 28033 was 28033, checked in by Ivan Raikov, 8 years ago

nemo: added several variants of membrane potential equation

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