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

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

nemo: updated copyright year

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