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

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

nemo: more uniform code generaton for octave and nmodl

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