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

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

nemo: bug fixes and error handling

File size: 35.4 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 srfi-69)
27       
28        (require-extension lolevel matchable strictly-pretty 
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      ))
43
44(define (matlab-name s)
45  (let ((cs (string->list (->string s))))
46    (let loop ((lst (list)) (cs cs))
47      (if (null? cs) (string->symbol (list->string (reverse lst)))
48          (let* ((c (car cs))
49                 (c1 (cond ((or (char-alphabetic? c) (char-numeric? c) (char=? c #\_)) c)
50                           (else #\_))))
51            (loop (cons c1 lst) (cdr cs)))))))
52           
53(define (rhsexpr/MATLAB expr)
54  (match expr 
55         (('if . es)  `(if . ,(map (lambda (x) (rhsexpr/MATLAB x)) es)))
56         (('pow x y)  (if (and (integer? y)  (positive? y))
57                          (if (> y 1)  (let ((tmp (gensym "x")))
58                                         `(let ((,tmp ,x)) (* . ,(list-tabulate (inexact->exact y) (lambda (i) tmp)))))
59                              x)
60                            (if (and (number? y) (zero? y)) 1.0 expr)))
61         ((s . es)    (if (symbol? s)  (cons (if (member s builtin-fns) s (matlab-name s)) (map (lambda (x) (rhsexpr/MATLAB x)) es)) expr))
62         (id          (if (symbol? id) (matlab-name id) id))))
63
64(define (matlab-state-name n s)
65  (matlab-name (s+ n s)))
66
67(define-syntax pp
68  (syntax-rules ()
69    ((pp indent val ...) (ppf indent (quasiquote val) ...))))
70
71
72(define group/MATLAB   (doc:block 2 (doc:text "(") (doc:text ")")))
73(define block/MATLAB   (doc:block 2 (doc:empty) (doc:empty)))
74(define (stmt/MATLAB x) 
75  (match x
76         (($ doc 'DocCons _ ($ doc 'DocText ";")) x)
77         (else  (doc:cons x (doc:text ";")))))
78
79
80(define (ifthen/MATLAB c e1 e2)
81  (doc:nest 
82   2 (doc:connect
83      (doc:connect (doc:group (doc:connect (doc:text "if") c))
84                   (doc:connect (doc:nest 2 e1)
85                                (doc:nest 2 (doc:connect (doc:text "else") e2))))
86      (doc:text "end"))))
87
88(define (letblk/MATLAB e1 e2)
89  (cond ((equal? e1 (doc:empty)) (doc:group (doc:nest 2 e2)))
90        ((equal? e2 (doc:empty)) (doc:group (doc:nest 2 e1)))
91        (else (doc:connect (doc:group (doc:nest 2 (stmt/MATLAB e1)))
92                           (doc:group (doc:nest 2 e2))))))
93       
94
95(define (format-op/MATLAB indent op args)
96  (let ((op1 (doc:text (->string op))))
97    (if (null? args) op1
98        (match args
99               ((x)      (doc:concat (list op1 x)))
100               ((x y)    (doc:concat (intersperse (list x op1 y) (doc:space))))
101               ((x y z)  (doc:concat (intersperse (list x op1 y op1 z) (doc:space))))
102               (lst      (let* ((n   (length lst))
103                                (n/2 (inexact->exact (round (/ n 2)))))
104                           (doc:concat 
105                            (intersperse 
106                             (list (format-op/MATLAB indent op (take lst n/2 )) op1 
107                                   (format-op/MATLAB indent op (drop lst n/2 )))
108                             (doc:space)))))))))
109
110(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 (hash-table-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 (hash-table-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 (hash-table-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 (hash-table-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 (hash-table-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 (hash-table-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 #t) (dirname "."))
682  (match-let ((($ nemo:quantity 'DISPATCH  dis) (hash-table-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      (->string sysname))
690             (deps*       ((dis 'depgraph*) sys))
691             (consts      ((dis 'consts)  sys))
692             (asgns       ((dis 'asgns)   sys))
693             (states      ((dis 'states)  sys))
694             (reactions   ((dis 'reactions) sys))
695             (defuns      ((dis 'defuns)  sys))
696             (components  ((dis 'components) sys))
697             
698             (g             (match-let (((state-list asgn-list g) ((dis 'depgraph*) sys))) g))
699             (poset         (vector->list ((dis 'depgraph->bfs-dist-poset) g)))
700
701             (const-defs       (filter-map
702                                (lambda (nv)
703                                  (and (not (member (first nv) matlab-builtin-consts))
704                                       (let ((v1 (canonicalize-expr/MATLAB (second nv))))
705                                         (list (matlab-name (first nv)) v1))))
706                                consts))
707             
708             (gate-complex-info    (nemo:gate-complex-query sys))
709
710             (gate-complexes       (lookup-def 'gate-complexes gate-complex-info))
711             (perm-ions     (map (match-lambda ((comp i e erev) `(,comp ,(matlab-name i) ,(matlab-name e) ,erev)))
712                                 (lookup-def 'perm-ions gate-complex-info)))
713             (acc-ions      (map (match-lambda ((comp i in out) `(,comp ,@(map matlab-name (list i in out)))))
714                                 (lookup-def 'acc-ions gate-complex-info)))
715             (epools        (lookup-def 'pool-ions gate-complex-info))
716             (pool-ions     (map (lambda (lst) (map matlab-name lst)) epools))
717
718             (i-gates       (lookup-def 'i-gates gate-complex-info))
719
720             (capcomp       (any (match-lambda ((name 'membrane-capacitance id) (list name id)) (else #f)) components))
721             (mcap          (and capcomp (car ((dis 'component-exports) sys (cid capcomp)))))
722               
723
724             (i-eqs (filter-map
725                     (lambda (gate-complex) 
726                       
727                       (let* ((label             (first gate-complex))
728                              (n                 (second gate-complex))
729                              (subcomps          ((dis 'component-subcomps) sys n))
730                              (acc               (lookup-def 'accumulating-substance subcomps))
731                              (perm              (lookup-def 'permeating-ion subcomps))
732                              (permqs            (and perm ((dis 'component-exports) sys (cid perm))))
733                              (pore              (lookup-def 'pore subcomps))
734                              (permeability      (lookup-def 'permeability subcomps))
735                              (gate              (lookup-def 'gate subcomps))
736                              (sts               (and gate ((dis 'component-exports) sys (cid gate)))))
737
738                         (if (null? permqs)
739                             (nemo:error 'nemo:matlab-translator ": ion channel definition " label
740                                         "permeating-ion component lacks exported quantities"))
741                         (if (null? sts)
742                             (nemo:error 'nemo:matlab-translator ": ion channel definition " label
743                                         "gate component lacks exported quantities"))
744                         
745                         (if (not (or pore permeability))
746                             (nemo:error 'nemo:matlab-translator ": ion channel definition " label
747                                         "lacks any pore or permeability components"))
748
749
750                         (cond ((and perm permeability gate)
751                                     (let* ((i     (matlab-name (s+ 'i (cn perm))))
752                                            (pmax  (car ((dis 'component-exports) sys (cid permeability))))
753                                            (pwrs  (map (lambda (n) (rate/reaction-power sys n)) sts))
754                                            (sptms (map (lambda (st pwr) (if pwr `(pow ,st ,pwr) st)) sts pwrs))
755                                            (gion  `(* ,pmax ,@sptms)))
756                                         (list i #f gion (matlab-name (s+ 'i_ label) ))))
757
758                               ((and perm pore gate)
759                                (case (cn perm)
760                                  ((non-specific)
761                                   (let* ((i     (matlab-name 'i))
762                                          (e     (car permqs))
763                                          (gmax  (car ((dis 'component-exports) sys (cid pore))))
764                                          (pwrs  (map (lambda (n) (rate/reaction-power sys n)) sts))
765                                          (sptms (map (lambda (st pwr) (if pwr `(pow ,st ,pwr) st)) sts pwrs))
766                                          (gion  `(* ,gmax ,@sptms)))
767                                     (list i e gion  (matlab-name (s+ 'i_ label) ))))
768
769                                  (else
770                                   (let* ((i     (matlab-name (s+ 'i (cn perm))))
771                                          (e     (car permqs))
772                                          (gmax  (car ((dis 'component-exports) sys (cid pore))))
773                                          (pwrs  (map (lambda (n) (rate/reaction-power sys n)) sts))
774                                          (sptms (map (lambda (st pwr) (if pwr `(pow ,st ,pwr) st)) sts pwrs))
775                                          (gion  `(* ,gmax ,@sptms)))
776                                     (list i e gion (matlab-name (s+ 'i_ label) ))))))
777                               
778                               ((and perm pore)
779                                (case (cn perm)
780                                  ((non-specific)
781                                   (let* ((i     (matlab-name 'i))
782                                          (e     (car permqs))
783                                          (gmax  (car ((dis 'component-exports) sys (cid pore)))))
784                                     (list i e gmax (matlab-name (s+ 'i_ label) ))))
785                                  (else
786                                   (nemo:error 'nemo:matlab-translator ": invalid ion channel definition " label))))
787                               
788                               ((and acc pore gate)
789                                (let* ((i     (matlab-name (s+ 'i (cn acc))))
790                                       (gmax  (car ((dis 'component-exports) sys (cid pore))))
791                                       (pwrs  (map (lambda (n) (rate/reaction-power sys n)) sts))
792                                       (sptms (map (lambda (st pwr) (if pwr `(pow ,st ,pwr) st)) sts pwrs))
793                                       (gion  `(* ,gmax ,@sptms)))
794                                  (list i #f gion  (matlab-name (s+ 'i_ label) ))))
795                               (else (nemo:error 'nemo:matlab-translator ": invalid ion channel definition " label))
796                               )))
797                     gate-complexes))
798
799             (i-names (delete-duplicates (map first i-eqs)))
800               
801             (i-eqs  (fold  (lambda (i-gate ax) 
802                              (let ((i-gate-var (first i-gate)))
803                                (cons (list (matlab-name 'i) #f i-gate-var (s+ 'i_ (second i-gate))) ax)))
804                            i-eqs i-gates))
805
806             (i-bkts (bucket-partition (lambda (x y) (eq? (car x) (car y))) i-eqs))
807             
808             (i-eqs  (fold (lambda (b ax) 
809                             (match b 
810                                    ((and ps ((i e gion ii) . rst)) 
811                                     (let loop ((ps ps) (summands (list)) (eqs (list)))
812                                       (if (null? ps)
813                                           
814                                           (let* ((sum0  (sum summands))
815                                                  (sum1  (rhsexpr/MATLAB sum0))
816                                                  (sum2  (canonicalize-expr/MATLAB sum1)))
817                                             (append eqs (list (list i sum2)) ax))
818                                           
819                                           (match-let (((i e gion ii) (car ps)))
820                                                      (loop (cdr ps) 
821                                                            (cons ii summands) 
822                                                            (let* ((expr0 (rhsexpr/MATLAB (if e `(* ,gion (- v ,e)) gion)))
823                                                                   (expr1 (canonicalize-expr/MATLAB expr0)))
824                                                              (cons (list ii expr1) eqs)))))))
825                                   
826                                    ((i e gion ii)
827                                     (let* ((expr0  (rhsexpr/MATLAB (if e `(* ,gion (- v ,e)) gion)))
828                                            (expr1  (canonicalize-expr/MATLAB expr0)))
829                                       (cons (list i expr1) ax)))
830                                   
831                                    (else ax)))
832                           (list) i-bkts))
833
834             (asgn-eq-defs     (poset->asgn-eq-defs poset sys))
835             
836             (rate-eq-defs       (reverse (poset->rate-eq-defs poset sys method)))
837             
838             (reaction-eq-defs   (poset->reaction-eq-defs poset sys))
839             
840             (init-eq-defs       (poset->init-defs poset sys))
841             
842             (conserve-eq-defs   (map (lambda (eq) (list 0 `(- ,(second eq) ,(first eq)))) 
843                                      (poset->state-conserve-eq-defs poset sys)))
844             
845             (globals          (map matlab-name 
846                                    (delete-duplicates (append
847                                                        exports
848                                                        (map second  perm-ions)
849                                                        (map third   perm-ions)
850                                                        (map second  acc-ions)
851                                                        (map third   acc-ions)
852                                                        (map fourth  acc-ions)
853                                                        (map second  pool-ions)
854                                                        (map third   pool-ions)
855                                                        (map first   imports) 
856                                                        (map first   const-defs)
857                                                        (map first   asgn-eq-defs)
858                                                        (map (lambda (gate-complex) (matlab-name (s+ 'i_ (first gate-complex)))) gate-complexes )
859                                                        (map (lambda (i-gate) (matlab-name (s+ 'i_ (second i-gate)))) i-gates )
860
861                                                        ))))
862             (v-eq    (if (and mcap (member 'v globals))
863                          (list 'v (rhsexpr/MATLAB `(/ (neg ,(sum i-names)) ,mcap)))
864                          (list 'v 0.0)))
865             
866             (state-index-map  (let ((acc (fold (lambda (def ax)
867                                                  (let ((st-name (first def)))
868                                                    (list (+ 1 (first ax)) 
869                                                          (cons `(,st-name ,(first ax)) (second ax)))))
870                                                (list 1 (list)) 
871                                                (if (member 'v globals)
872                                                    (cons (list 'v) rate-eq-defs)
873                                                    rate-eq-defs)
874                                                )))
875                                 (second acc)))
876             
877             (steady-state-index-map  (let ((acc (fold (lambda (def ax)
878                                                         (let ((st-name (first def)))
879                                                           (if (not (alist-ref st-name init-eq-defs))
880                                                               (list (+ 1 (first ax)) 
881                                                                     (cons `(,st-name ,(first ax)) (second ax)))
882                                                               ax)))
883                                                       (list 1 (list)) 
884                                                       rate-eq-defs)))
885                                        (second acc)))
886             
887             (dfenv 
888              (map (lambda (x) (let ((n (first x)))
889                                 (list n (matlab-name (s+ "d_" n )))))
890                   defuns))
891
892             (output-globals 
893              (lambda () (if (not (null? globals)) (pp indent (global ,(slp " " globals))))))
894             )
895       
896       
897       
898        (for-each
899         (lambda (a)
900           (let ((acc-ion   (car a)))
901             (if (assoc acc-ion perm-ions)
902                 (nemo:error 'nemo:matlab-translator 
903                             ": ion species " acc-ion " cannot be declared as both accumulating and permeating"))))
904         acc-ions)
905
906           (for-each
907            (lambda (p)
908              (let ((pool-ion  (car p)))
909                (if (assoc pool-ion perm-ions)
910                    (nemo:error 'nemo:matlab-translator 
911                                ": ion species " pool-ion " cannot be declared as both pool and permeating"))))
912            pool-ions)
913
914           (let ((output  (case mode
915                            ((single)  (open-output-file (make-output-fname dirname prefix ".m" filename)))
916                            (else #f))))
917
918             (if output (with-output-to-port output (lambda () (output-globals))))
919             
920             ;; derivative function
921             (let ((output1 (or output (open-output-file (make-output-fname dirname prefix "_dy.m")))))
922               (with-output-to-port output1
923                 (lambda ()
924                   (output-dy sysname method globals state-index-map 
925                              rate-eq-defs reaction-eq-defs asgn-eq-defs
926                              pool-ions mcap i-eqs v-eq
927                              indent indent+)))
928               (if (not output) (close-output-port output1)))
929                 
930
931             ;; steady state solver
932             (let ((output1 (or output (open-output-file (make-output-fname dirname prefix "_steadystate.m")))))
933               (with-output-to-port output1
934                 (lambda ()
935                   (output-steadystate sysname method globals steady-state-index-map pool-ions
936                                       const-defs init-eq-defs asgn-eq-defs rate-eq-defs reaction-eq-defs 
937                                       indent indent+)))
938               (if (not output) (close-output-port output1)))
939
940             ;; initial values function
941             (let ((output1 (or output (open-output-file (make-output-fname dirname prefix "_init.m")))))
942               (with-output-to-port output1
943                 (lambda ()
944                   (output-init sysname globals state-index-map steady-state-index-map 
945                                const-defs asgn-eq-defs init-eq-defs rate-eq-defs 
946                                reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
947                   (pp indent ,nl)))
948               (if (not output) (close-output-port output1)))
949
950           ;; user-defined functions
951           (let* (;;(with       (inexact->exact (round (/ (abs (- max-v min-v)) step))))
952                  (define-fn  (make-define-fn)))
953             (for-each (lambda (fndef) 
954                         (if (not (member (car fndef) builtin-fns))
955                             (let ((output1 (or output (open-output-file (make-output-fname dirname (matlab-name (car fndef)) ".m")))))
956                               (with-output-to-port output1
957                                 (lambda ()
958                                   (apply define-fn (cons indent fndef))
959                                   (pp indent ,nl)))
960                               (if (not output) (close-output-port output1)))))
961                       defuns))
962         
963           (if output (close-output-port output)))
964
965           )))))
966
967
968(define (nemo:matlab-translator sys . rest)
969  (apply matlab-translator1 (cons sys (cons 'multiple (cons 'odepkg rest)))))
970
971(define (nemo:octave-translator sys . rest)
972  (apply matlab-translator1 (cons sys (cons 'single rest))))
973
974#|
975 
976(define (make-define-dfn fenv)
977
978  (define (D vars body)
979    (let loop ((vars vars) (exprs (list)))
980      (if (null? vars) (if (pair? (cdr exprs))
981                           `(+ . ,(reverse exprs))  (car exprs))
982          (let ((de (let ((de (differentiate fenv (car vars) body)))
983                       (let loop ((e (simplify de)) (de de))
984                         (if (equal? de e) de
985                             (loop (simplify e) e))))))
986            (loop (cdr vars) (cons de exprs))))))
987
988  (lambda (indent n proc)
989    (let ((lst (procedure-data proc))
990          (indent+ (+ 2 indent)))
991      (let ((retval   (matlab-name (gensym "retval")))
992            (rt       (lookup-def 'rt lst))
993            (formals  (lookup-def 'formals lst))
994            (vars     (lookup-def 'vars lst))
995            (body     (lookup-def 'body lst)))
996        (pp indent ,nl (function ,retval = ,(s+ "d_" (matlab-name n) ) (,(slp ", " vars)) ))
997        (let* ((dbody (D vars body))
998               (dbody (canonicalize-expr/MATLAB (rhsexpr/MATLAB dbody))))
999          (pp indent+ ,(expr->string/MATLAB dbody retval))
1000          (pp indent end))
1001          )))
1002)
1003
1004(define (output-jac sysname method globals fenv state-index-map
1005                    rate-eq-defs reaction-eq-defs asgn-eq-defs
1006                    mcap i-eqs v-eq
1007                    defuns indent indent+)
1008
1009
1010  (define (D var expr)
1011    (let* ((dexpr  (distribute expr))
1012           (de (differentiate fenv var dexpr))
1013           (se (let loop ((e (simplify de)) (de de))
1014                 (if (equal? de e) de
1015                     (loop (simplify e) e)))))
1016      se))
1017   
1018  (pp indent ,nl (function jac = ,(s+ sysname "_jac")
1019                           (,(slp ", " `(t y)))))
1020
1021  (if (not (null? globals)) (pp indent+ (global ,(slp " " globals))))
1022 
1023  (if (member 'v globals)
1024      (let ((vi (lookup-def 'v state-index-map)))
1025        (if vi (pp indent+ ,(expr->string/MATLAB (s+ "y(" vi ")") 'v)))))
1026 
1027  (for-each (lambda (def)
1028              (let ((name (first def)))
1029                (pp indent+ ,(expr->string/MATLAB
1030                              (s+ "y(" (lookup-def name state-index-map) ")") name))))
1031            rate-eq-defs)
1032
1033  (for-each (lambda (def)
1034              (let ((n (first def)) (b (second def)))
1035                (pp indent+
1036                    ,(expr->string/MATLAB b n))))
1037            reaction-eq-defs)
1038 
1039  (pp indent+ ,(expr->string/MATLAB `(zeros (length y) (length y)) 'jac))
1040
1041
1042  (for-each
1043   (lambda (def)
1044     
1045     (let* ((name  (first def))
1046            (rhs   (second def))
1047            (fv    (enum-freevars rhs '() '() ))
1048            (drvs  (map (lambda (st)
1049                          (let ((var (first st))
1050                                (expr `(let ,(or (and v-eq (list v-eq)) '())
1051                                         (let ,rate-eq-defs
1052                                           (let ,reaction-eq-defs
1053                                             (let ,(filter-map (lambda (x)
1054                                                                 (or (assoc x i-eqs)
1055                                                                       (assoc x asgn-eq-defs)))
1056                                                               fv))
1057                                             ,rhs)))))
1058                            (D var expr)))
1059                        state-index-map))
1060            (vs (map (lambda (st)
1061                       (let ((ni (lookup-def name state-index-map))
1062                             (i (second st)))
1063                         (s+ "jac(" (slp "," (list ni i)) ")")))
1064                     state-index-map))
1065            (ndrvs  (map (lambda (x) (canonicalize-expr/MATLAB (rhsexpr/MATLAB x))) drvs)))
1066       
1067       (pp indent+ ,(s+ "% state " name))
1068       (for-each (lambda (ndrv v)
1069                   (pp indent+ ,(expr->string/MATLAB ndrv v)))
1070                 ndrvs vs)
1071       ))
1072           
1073     (if (and mcap (member 'v globals))
1074         (let ((vi (lookup-def 'v state-index-map)))
1075           (if vi (cons
1076                   (list 'v `(/ (neg ,(sum (map first i-eqs))) ,mcap))
1077                   rate-eq-defs)
1078               rate-eq-defs))
1079         rate-eq-defs))
1080 
1081  (pp indent ,nl (end))
1082 
1083)
1084
1085|#
1086
1087)
Note: See TracBrowser for help on using the repository browser.