source: project/release/3/nemo/trunk/nemo-matlab.scm @ 12554

Last change on this file since 12554 was 12554, checked in by Ivan Raikov, 11 years ago

Many fixes to the matlab code generator.

File size: 27.6 KB
Line 
1;;       
2;;
3;; An extension for translating NEMO models to Matlab code.
4;;
5;; Copyright 2008 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(require-extension syntax-case)
22(require-extension matchable)
23(require-extension strictly-pretty)
24(require-extension environments)
25(require-extension nemo-core)
26(require-extension srfi-1)
27(require-extension srfi-4)
28(require-extension srfi-13)
29(require-extension srfi-14)
30(require-extension utils)
31(require-extension lolevel)
32(require-extension varsubst)
33(require-extension digraph)
34(require-extension datatype)
35
36(define-extension nemo-matlab)
37
38(declare
39 (lambda-lift)
40 (export  nemo:matlab-translator))
41
42(register-feature! 'nemo-matlab)
43
44(define (lookup-def k lst . rest)
45  (let-optionals rest ((default #f))
46      (let ((kv (assoc k lst)))
47        (if (not kv) default
48            (match kv ((k v) v) (else (cdr kv)))))))
49
50(define (matlab-name s)
51  (let ((cs (string->list (->string s))))
52    (let loop ((lst (list)) (cs cs))
53      (if (null? cs) (string->symbol (list->string (reverse lst)))
54          (let* ((c (car cs))
55                 (c1 (cond ((or (char-alphabetic? c) (char-numeric? c) (char=? c #\_)) c)
56                           (else #\_))))
57            (loop (cons c1 lst) (cdr cs)))))))
58                           
59                 
60
61(define (matlab-state-name n s)
62  (matlab-name (s+ n s)))
63
64(define (enum-bnds expr ax)
65  (match expr 
66         (('if . es)        (fold enum-bnds ax es))
67         (('let bnds body)  (enum-bnds body (append (map car bnds) (fold enum-bnds ax (map cadr bnds)))))
68         ((s . es)          (if (symbol? s)  (fold enum-bnds ax es) ax))
69         (else ax)))
70
71
72(define (enum-freevars expr bnds ax)
73  (match expr 
74         (('if . es) 
75          (fold (lambda (x ax) (enum-freevars x bnds ax)) ax es))
76         (('let bnds body) 
77          (let ((bnds1 (append (map first bnds) bnds)))
78            (enum-freevars body bnds1 (fold (lambda (x ax) (enum-freevars x bnds ax)) ax (map second bnds)))))
79         ((s . es)    (if (symbol? s)  (fold (lambda (x ax) (enum-freevars x bnds ax)) ax es) ax))
80         (id          (if (and (symbol? id) (not (member id bnds)))  (cons id ax) ax))))
81
82(define (rhsvars rhs)
83  (enum-freevars rhs (list) (list)))
84
85(define (rhsexpr expr)
86  (match expr 
87         (('if . es)  `(if . ,(map (lambda (x) (rhsexpr x)) es)))
88         (('pow x y)  (if (and (integer? y)  (positive? y))
89                          (if (> y 1)  (let ((tmp (gensym "x")))
90                                         `(let ((,tmp ,x)) (* . ,(list-tabulate (inexact->exact y) (lambda (i) tmp)))))
91                              x)
92                            expr))
93         ((s . es)    (if (symbol? s)  (cons s (map (lambda (x) (rhsexpr x)) es)) expr))
94         (id          (if (symbol? id) (matlab-name id) id))))
95
96
97;;; Procedures for string concatenation and pretty-printing
98
99(define (s+ . lst)    (string-concatenate (map ->string lst)))
100(define (sw+ lst)     (string-intersperse (filter-map (lambda (x) (and x (->string x))) lst) " "))
101(define (sl\ p lst)   (string-intersperse (map ->string lst) p))
102(define nl "\n")
103
104(define (spaces n)  (list->string (list-tabulate n (lambda (x) #\space))))
105
106(define (ppf indent . lst)
107  (let ((sp (spaces indent)))
108    (for-each (lambda (x)
109                (and x (match x 
110                              ((i . x1) (if (and (number? i) (positive? i))
111                                            (for-each (lambda (x) (ppf (+ indent i) x)) x1)
112                                            (print sp (sw+ x))))
113                              (else   (print sp (if (list? x) (sw+ x) x))))))
114              lst)))
115
116
117(define-syntax pp
118  (syntax-rules ()
119    ((pp indent val ...) (ppf indent (quasiquote val) ...))))
120
121(define (ifthen/MATLAB c e1 e2)
122  (doc:nest 
123   2 (doc:connect
124      (doc:connect (doc:group (doc:connect (doc:text "if") c))
125                   (doc:connect (doc:nest 2 e1)
126                                (doc:nest 2 (doc:connect (doc:text "else") e2))))
127      (doc:text "endif"))))
128
129
130(define (letblk/MATLAB e1 e2)
131  (cond ((equal? e1 (doc:empty)) (doc:group (doc:nest 2 e2)))
132        ((equal? e2 (doc:empty)) (doc:group (doc:nest 2 e1)))
133        (else (doc:connect (doc:group (doc:concat (list (doc:nest 2 e1) (doc:text ";"))))
134                           (doc:group (doc:nest 2 e2))))))
135       
136
137(define group/MATLAB   (doc:block 2 (doc:text "(") (doc:text ")")))
138(define block/MATLAB   (doc:block 2 (doc:empty) (doc:empty)))
139
140(define (format-op/MATLAB indent op args)
141  (let ((op1 (doc:text (->string op))))
142    (if (null? args) op1
143        (match args
144               ((x)      (doc:concat (list op1 x)))
145               ((x y)    (doc:concat (intersperse (list x op1 y) (doc:space))))
146               ((x y z)  (doc:concat (intersperse (list x op1 y op1 z) (doc:space))))
147               (lst      (let* ((n   (length lst))
148                                (n/2 (inexact->exact (round (/ n 2)))))
149                           (doc:concat 
150                            (intersperse 
151                             (list (format-op/MATLAB indent op (take lst n/2 )) op1 
152                                   (format-op/MATLAB indent op (drop lst n/2 )))
153                             (doc:space)))))))))
154
155(define (format-fncall/MATLAB indent op args)
156  (let ((op1 (doc:text (->string op))))
157    (doc:cons op1 (group/MATLAB ((doc:list indent identity (lambda () (doc:text ", "))) args)))))
158
159(define matlab-builtin-consts
160  `())
161
162(define matlab-ops
163  `(+ - * / > < <= >= = ^))
164
165(define builtin-fns
166  `(+ - * / pow neg abs atan asin acos sin cos exp ln
167      sqrt tan cosh sinh tanh hypot gamma lgamma log10 log2 log1p ldexp cube
168      > < <= >= = and or round ceiling floor max min
169      fpvector-ref))
170
171(define (sum lst)
172  (if (null? lst) lst
173      (match lst
174             ((x)   x)
175             ((x y) `(+ ,x ,y))
176             ((x y . rest) `(+ (+ ,x ,y) ,(sum rest)))
177             ((x . rest) `(+ ,x ,(sum rest))))))
178
179
180(define (subst-term t subst k)
181  (match t
182         (('if c t e)
183          `(if ,(k c subst) ,(k t subst) ,(k e subst)))
184         (('let bs e)
185          (k `(let ,(map (lambda (b) `(,(car b) ,(k (cadr b) subst))) bs) ,(k e subst)) subst))
186         ((f . es)
187          (cons (k f subst) (map (lambda (e) (k e subst)) es)))
188         ((? symbol? )  (lookup-def t subst t))
189         ((? atom? ) t)))
190
191(define (binding? t) (and (list? t) (eq? 'let (car t)) (cdr t)))
192
193(define (bind ks vs e) `(let ,(zip ks vs) ,e))
194
195(define (name-normalize expr)
196  (match expr 
197         (('if c t e)  `(if ,(name-normalize c) ,(name-normalize t) ,(name-normalize e)))
198         (('let bs e)
199          `(let ,(map (lambda (b) `(,(car b) ,(name-normalize (cadr b)))) bs) ,(name-normalize e)))
200         ((f . es) 
201          (cons f (map name-normalize es)))
202         ((? symbol? ) (matlab-name expr))
203         ((? atom? ) expr)))
204
205(define (if-convert expr)
206  (match expr 
207         (('if c t e) 
208          (let ((r (gensym "if")))
209            `(let ((,r (if ,(if-convert c) ,(if-convert t) ,(if-convert e)))) 
210               ,r)))
211         (('let bs e)
212          `(let ,(map (lambda (b) `(,(car b) ,(if-convert (cadr b)))) bs) ,(if-convert e)))
213         ((f . es)
214          (cons f (map if-convert es)))
215         ((? atom? ) expr)))
216
217         
218(define (let-enum expr ax)
219  (match expr
220         (('let ((x ('if c t e))) y)
221          (let ((ax (fold let-enum ax (list c ))))
222            (if (eq? x y)  (append ax (list (list x `(if ,c ,t ,e)))) ax)))
223
224         (('let bnds body)  (let-enum body (append ax bnds)))
225
226         (('if c t e)  (let-enum ax c))
227
228         ((f . es)  (fold let-enum ax es))
229
230         (else ax)))
231
232
233(define (let-elim expr)
234  (match expr
235         (('let ((x ('if c t e))) y)
236          (if (eq? x y)  y expr))
237
238         (('let bnds body) (let-elim body))
239
240         (('if c t e)  `(if ,(let-elim c) ,(let-lift t) ,(let-lift e)))
241
242         ((f . es)  `(,f . ,(map let-elim es)))
243
244         (else expr)))
245 
246
247(define (let-lift expr)
248  (let ((bnds (let-enum expr (list))))
249    (if (null? bnds) expr
250        `(let ,(map (lambda (b) (list (car b) (let-elim (cadr b)))) bnds) ,(let-elim expr)))))
251
252(define (canonicalize-expr/MATLAB expr)
253  (let ((subst-convert  (subst-driver (lambda (x) (and (symbol? x) x)) binding? identity bind subst-term)))
254    (let* ((expr1 (if-convert expr))
255           (expr2 (subst-convert expr1 subst-empty))
256           (expr3 (let-lift expr2))
257           (expr4 (name-normalize expr3)))
258      expr4)))
259
260
261(define (format-expr/MATLAB indent expr . rest) 
262  (let-optionals rest ((rv #f))
263   (let ((indent+ (+ 2 indent)))
264    (match expr
265       (('let bindings body)
266        (letblk/MATLAB
267         (fold-right 
268           (lambda (x ax)
269             (letblk/MATLAB
270              (match (second x)
271                     (('if c t e)
272                      (ifthen/MATLAB
273                       (group/MATLAB (format-expr/MATLAB indent c))
274                       (block/MATLAB (format-expr/MATLAB indent t (first x)))
275                       (block/MATLAB (format-expr/MATLAB indent e (first x)))))
276                     (else
277                      (format-op/MATLAB indent+ " = "
278                                       (list (format-expr/MATLAB indent (first x) )
279                                             (format-expr/MATLAB indent (second x))))))
280              ax))
281           (doc:empty) bindings)
282         (let ((body1 (doc:nest indent (format-expr/MATLAB indent body))))
283           (if rv  (format-op/MATLAB indent " = " (list (format-expr/MATLAB indent+ rv ) body1))
284               body1))))
285       
286       (('if . rest) (error 'format-expr/MATLAB "invalid if statement " expr))
287
288       ((op . rest) 
289       (let ((op (case op ((pow)  '^)  (else op))))
290         (let ((fe
291                (if (member op matlab-ops)
292                    (let ((mdiv?  (any (lambda (x) (match x (('* . _) #t) (('/ . _) #t) (else #f))) rest))
293                          (mul?   (any (lambda (x) (match x (('* . _) #t) (else #f))) rest))
294                          (plmin? (any (lambda (x) (match x (('+ . _) #t) (('- . _) #t) (else #f))) rest)))
295                      (case op
296                        ((/) 
297                         (format-op/MATLAB indent op 
298                                          (map (lambda (x) 
299                                                 (let ((fx (format-expr/MATLAB indent+ x)))
300                                                   (if (or (symbol? x) (number? x)) fx
301                                                       (if (or mul? plmin?) (group/MATLAB fx) fx)))) rest)))
302                        ((*) 
303                         (format-op/MATLAB indent op 
304                                          (map (lambda (x) 
305                                                 (let ((fx (format-expr/MATLAB indent+ x)))
306                                                   (if (or (symbol? x) (number? x)) fx
307                                                       (if plmin? (group/MATLAB fx) fx)))) rest)))
308                       
309                        ((^) 
310                         (format-op/MATLAB indent op 
311                                          (map (lambda (x) 
312                                                 (let ((fx (format-expr/MATLAB indent+ x)))
313                                                   (if (or (symbol? x)  (number? x)) fx
314                                                       (if (or mdiv? plmin?) (group/MATLAB fx) fx)))) rest)))
315                       
316                        (else
317                         (format-op/MATLAB indent op 
318                                          (map (lambda (x) 
319                                                 (let ((fx (format-expr/MATLAB indent+ x))) fx)) rest)))))
320                   
321                    (let ((op (case op ((neg) '-) (else op))))
322                      (format-fncall/MATLAB indent op (map (lambda (x) (format-expr/MATLAB indent+ x)) rest))))))
323           (if rv 
324               (format-op/MATLAB indent " = " (list (format-expr/MATLAB indent+ rv ) fe))
325               fe))))
326     
327      (else  (let ((fe (doc:text (->string expr))))
328               (if rv 
329                   (format-op/MATLAB indent " = " (list (format-expr/MATLAB indent+ rv ) fe))
330                   fe)))))))
331               
332
333         
334(define (expr->string/MATLAB x . rest)
335  (let-optionals rest ((rv #f) (width 72))
336    (sdoc->string (doc:format width (format-expr/MATLAB 2 x rv)))))
337 
338(define (make-define-fn table? min-v max-v with)
339  (lambda (indent n proc)
340    (let ((lst (procedure-data proc))
341          (indent+ (+ 2 indent)))
342      (let ((retval   (matlab-name (gensym "retval")))
343            (rt       (lookup-def 'rt lst))
344            (formals  (lookup-def 'formals lst))
345            (vars     (lookup-def 'vars lst))
346            (body     (lookup-def 'body lst)))
347        (pp indent ,nl (function ,retval = ,n (,(sl\ ", " vars)) ))
348        (let* ((body1 (canonicalize-expr/MATLAB (rhsexpr body)))
349               (lbs   (enum-bnds body1 (list))))
350          (pp indent+ ,(expr->string/MATLAB body1 retval))
351          (pp indent endfunction))
352          ))))
353
354(define (define-state indent n)
355  (pp indent (,n)))
356
357
358(define (state-eqs n initial open transitions power)
359  (let* ((subst-convert  (subst-driver (lambda (x) (and (symbol? x) x)) binding? identity bind subst-term))
360         (g          (make-digraph n (string-append (->string n) " transitions graph")))
361         (add-node!  (g 'add-node!))
362         (add-edge!  (g 'add-edge!))
363         (out-edges  (g 'out-edges))
364         (in-edges   (g 'in-edges))
365         (node-info  (g 'node-info))
366         (node-list  (let loop ((lst (list)) (tlst transitions))
367                       (if (null? tlst)  (delete-duplicates lst eq?)
368                           (match (car tlst) 
369                                  (('-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr)
370                                   (loop (cons* s0 s1 lst) (cdr tlst)))
371                                  (((and (? symbol?) s0) '-> (and (? symbol? s1)) rate-expr)
372                                   (loop (cons* s0 s1 lst) (cdr tlst)))
373                                  (('<-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr1 rate-expr2)
374                                   (loop (cons* s0 s1 lst) (cdr tlst)))
375                                  (((and (? symbol?) s0) 'M-> (and (? symbol? s1)) rate-expr1 rate-expr2)
376                                   (loop (cons* s0 s1 lst) (cdr tlst)))
377                                  (else
378                                   (nemo:error 'nemo:matlab-state-eqs ": invalid transition equation " 
379                                                  (car tlst) " in state complex " n))
380                                  (else (loop lst (cdr tlst)))))))
381         (node-ids      (list-tabulate (length node-list) identity))
382         (name->id-map  (zip node-list node-ids))
383         (node-subs     (fold (lambda (s ax) (subst-extend s (matlab-state-name n s) ax)) subst-empty node-list)))
384    ;; insert state nodes in the dependency graph
385    (for-each (lambda (i n) (add-node! i n)) node-ids node-list)
386    (let* ((nodes  ((g 'nodes)))
387           (snode   (find (lambda (s) (not (eq? (second s) open))) nodes))
388           (snex   `(- 1 ,(sum (filter-map (lambda (s) (and (not (= (first s) (first snode))) (second s))) nodes))))
389           (add-tredge (lambda (s0 s1 rexpr1 rexpr2)
390                         (let ((i   (car (alist-ref s0 name->id-map)))
391                               (j   (car (alist-ref s1 name->id-map)))
392                               (x0  (if (eq? s0 (second snode)) snex s0))
393                               (x1  (if (eq? s1 (second snode)) snex s1)))
394                           (add-edge! (list i j `(* ,(subst-convert x0 node-subs) 
395                                                    ,(subst-convert rexpr1 node-subs))))
396                           (if rexpr2 (add-edge! (list j i `(* ,(subst-convert x1 node-subs) 
397                                                               ,(subst-convert rexpr2 node-subs)))))))))
398      ;; create rate edges in the graph
399      (for-each (lambda (e) 
400                  (match e
401                         (('-> s0 s1 rexpr)  (add-tredge s0 s1 rexpr #f))
402                         ((s0 '-> s1 rexpr)  (add-tredge s0 s1 rexpr #f))
403                         (('<-> s0 s1 rexpr1 rexpr2)  (add-tredge s0 s1 rexpr1 rexpr2))
404                         ((s0 '<-> s1 rexpr1 rexpr2)  (add-tredge s0 s1 rexpr1 rexpr2))
405                         ))
406                transitions)
407
408      ;; generate differential equations for each state in the transitions system
409      (let ((eqs    (fold (lambda (s ax) 
410                            (if (= (first snode) (first s) ) ax
411                                (let* ((out   (out-edges (first s)))
412                                       (in    (in-edges (first s)))
413                                       (open? (eq? (second s) open))
414                                       (name  (matlab-name (lookup-def (second s) node-subs))))
415                                  (let* ((rhs1  (cond ((and (not (null? out)) (not (null? in)))
416                                                       `(+ (neg ,(sum (map third out)))
417                                                           ,(sum (map third in))))
418                                                      ((and (not (null? out)) (null? in))
419                                                       `(neg ,(sum (map third out))))
420                                                      ((and (null? out) (not (null? in)))
421                                                       (sum (map third in)))))
422                                         (fbody  (rhsexpr rhs1))
423                                         (fbody1 (canonicalize-expr/MATLAB fbody)))
424                                    (cons (list name  fbody1) ax)))))
425                          (list) nodes)))
426        eqs))))
427           
428       
429
430
431(define (state-init n init)
432  (let* ((init  (rhsexpr init))
433         (init1 (canonicalize-expr/MATLAB init)))
434    (list (matlab-name n) init1)))
435
436
437
438(define (asgn-eq n rhs)
439  (let* ((fbody   (rhsexpr rhs))
440         (fbody1  (canonicalize-expr/MATLAB fbody)))
441    (list (matlab-name n) fbody1)))
442
443
444(define (stcomp-eq n open transitions)
445  (list (matlab-name n) (matlab-name (matlab-state-name n open))))
446
447
448(define (poset->asgn-eq-defs poset sys)
449  (fold-right
450   (lambda (lst ax)
451     (fold  (lambda (x ax) 
452              (match-let (((i . n)  x))
453                         (let ((en (environment-ref sys n)))
454                           (if (nemo:quantity? en)
455                               (cases nemo:quantity en
456                                      (ASGN  (name value rhs) (cons (asgn-eq name rhs) ax))
457                                      (else  ax))
458                               ax))))
459            ax lst))
460   (list) poset))
461
462
463(define (poset->state-eq-defs poset sys)
464  (fold-right
465   (lambda (lst ax)
466     (fold  (lambda (x ax) 
467              (match-let (((i . n)  x))
468                         (let ((en (environment-ref sys n)))
469                           (if (nemo:quantity? en)
470                               (cases nemo:quantity en
471                                      (TSCOMP  (name initial open transitions conserve power) 
472                                               (append (state-eqs name initial open transitions power) ax))
473                                      (else  ax))
474                               ax))))
475            ax lst))
476   (list) poset))
477
478
479
480
481(define (poset->stcomp-eq-defs poset sys)
482  (fold-right
483   (lambda (lst ax)
484     (fold  (lambda (x ax) 
485              (match-let (((i . n)  x))
486                         (let ((en (environment-ref sys n)))
487                           (if (nemo:quantity? en)
488                               (cases nemo:quantity en
489                                      (TSCOMP  (name initial open transitions conserve power) 
490                                               (cons (stcomp-eq name open transitions) ax))
491                                      (else  ax))
492                               ax))))
493            ax lst))
494   (list) poset))
495
496
497(define (poset->state-init-defs poset sys)
498  (fold-right
499   (lambda (lst ax)
500     (fold  (lambda (x ax) 
501              (match-let (((i . n)  x))
502                         (let ((en (environment-ref sys n)))
503                           (if (nemo:quantity? en)
504                               (cases nemo:quantity en
505                                      (TSCOMP  (name initial open transitions conserve power)
506                                               (if (nemo:rhs? initial)
507                                                   (cons* (state-init name initial) 
508                                                          (state-init (matlab-state-name name open) name) ax) 
509                                                   ax))
510                                      (else  ax))
511                               ax))))
512            ax lst))
513   (list) poset))
514
515
516(define (find-locals defs)
517  (concatenate (map (lambda (def) (match def (('let bnds _) (map first bnds)) (else (list)))) defs)))
518
519
520(define (state-power sys n)
521  (let ((en (environment-ref sys n)))
522    (if (nemo:quantity? en)
523        (cases nemo:quantity en
524               (TSCOMP  (name initial open transitions conserve power)  power)
525               (else  #f))  #f)))
526
527
528(define (bucket-partition p lst)
529  (let loop ((lst lst) (ax (list)))
530    (if (null? lst) ax
531        (let ((x (car lst)))
532          (let bkt-loop ((old-bkts ax) (new-bkts (list)))
533            (if (null? old-bkts) (loop (cdr lst) (cons (list x) new-bkts))
534                (if (p x (caar old-bkts ))
535                    (loop (cdr lst) (append (cdr old-bkts) (cons (cons x (car old-bkts)) new-bkts)))
536                    (bkt-loop (cdr old-bkts) (cons (car old-bkts) new-bkts)))))))))
537
538
539(define (collect-epools sys)
540   (match-let ((($ nemo:quantity 'DISPATCH  dis) (environment-ref sys (nemo-intern 'dispatch))))
541      (let recur ((comp-name (nemo-intern 'toplevel)) (ax (list)))
542        (let* ((comp-symbols   ((dis 'component-symbols)   sys comp-name))
543               (subcomps       ((dis 'component-subcomps)  sys comp-name)))
544          (fold recur 
545                (fold (lambda (sym ax)
546                        (let ((en (environment-ref sys sym)))
547                          (match en
548                                 ((or (('decaying 'pool)  ('name (? symbol? ion)) . alst)
549                                      (('decaying-pool)   ('name (? symbol? ion)) . alst))
550                                  (cons (list ion alst) ax))
551                                 (else ax)))) ax comp-symbols)
552                (map third subcomps))))))
553
554
555(define (nemo:matlab-translator sys . rest)
556  (define (cid x)  (second x))
557  (define (cn x)   (first x))
558  (let-optionals rest ((method 'cnexp) (table? #f) (min-v -100) (max-v 100) (step 0.5) )
559  (match-let ((($ nemo:quantity 'DISPATCH  dis) (environment-ref sys (nemo-intern 'dispatch))))
560    (let ((imports  ((dis 'imports)  sys))
561          (exports  ((dis 'exports)  sys)))
562      (let* ((indent      0)
563             (indent+     (+ 2 indent ))
564             (eval-const  (dis 'eval-const))
565             (sysname     (matlab-name ((dis 'sysname) sys)))
566             (deps*       ((dis 'depgraph*) sys))
567             (consts      ((dis 'consts)  sys))
568             (asgns       ((dis 'asgns)   sys))
569             (states      ((dis 'states)  sys))
570             (stcomps     ((dis 'stcomps) sys))
571             (defuns      ((dis 'defuns)  sys))
572             (components  ((dis 'components) sys))
573             (ionchs      (filter-map (match-lambda ((name 'ion-channel id) (list name id)) (else #f)) components))
574             (epools      (collect-epools sys)))
575
576
577        (match-let (((state-list asgn-list g) deps*))
578         (let* ((poset            (vector->list ((dis 'depgraph->bfs-dist-poset) g)))
579                (asgn-eq-defs     (poset->asgn-eq-defs poset sys))
580                (perm-ions (fold (lambda (ionch ax) 
581                                    (let* ((subcomps ((dis 'component-subcomps) sys (cid ionch)))
582                                           (perm      (lookup-def 'permeating-substance subcomps)))
583                                      (if perm 
584                                          (case (cn perm)
585                                            ((non-specific)   
586                                             (let* ((erev (car ((dis 'component-exports) sys (cid perm))))
587                                                    (i    (matlab-name 'i))
588                                                    (e    (matlab-name 'e)))
589                                               (cons `(,(cn perm) ,i ,e ,erev) ax)))
590                                            (else (let* ((erev (car ((dis 'component-exports) sys (cid perm))))
591                                                         (i    (matlab-name (s+ 'i (cn perm))))
592                                                         (e    (matlab-name (s+ 'e (cn perm)))))
593                                                    (cons `(,(cn perm) ,i ,e ,erev) ax))))
594                                          ax)))
595                                  (list) ionchs))
596               (acc-ions (fold (lambda (ionch ax) 
597                                  (let* ((subcomps ((dis 'component-subcomps) sys (cid ionch)))
598                                         (acc   (lookup-def 'accumulating-substance subcomps))
599                                         (i     (and acc (matlab-name (s+ 'i (cn acc)))))
600                                         (in    (and acc (matlab-name (s+ (cn acc) 'i))))
601                                         (out   (and acc (matlab-name (s+ (cn acc) 'o)))))
602                                    (if acc  (cons `(,(cn acc) ,i ,in ,out) ax) ax)))
603                                (list) ionchs))
604               (pool-ions (map (lambda (ep)
605                                  (let ((ion (car ep)))
606                                    `(,(matlab-name ion) ,(matlab-name (s+ 'i ion)) ,(matlab-name (s+ ion 'i)))))
607                               epools)))
608
609           (for-each
610            (lambda (a)
611              (let ((acc-ion   (car a)))
612                (if (assoc acc-ion perm-ions)
613                    (nemo:error 'nemo:matlab-translator 
614                                ": ion species " acc-ion " cannot be declared as both accumulating and permeating"))))
615            acc-ions)
616
617           (for-each
618            (lambda (p)
619              (let ((pool-ion  (car p)))
620                (if (assoc pool-ion perm-ions)
621                    (nemo:error 'nemo:matlab-translator 
622                                ": ion species " pool-ion " cannot be declared as both pool and permeating"))))
623            pool-ions)
624
625           (let* ((with      (inexact->exact (round (/ (abs (- max-v min-v)) step))))
626                  (define-fn (make-define-fn table? min-v max-v with)))
627             (for-each (lambda (fndef) 
628                         (if (not (member (car fndef) builtin-fns))
629                             (apply define-fn (cons indent fndef)))) 
630                       defuns))
631
632           (pp indent ,nl (function dy = ,sysname (,(sl\ ", " `(dy t)))))
633           
634           (if (not (null? exports)) (pp indent+ (global ,(sl\ ", " exports))))
635
636           (let* ((const-defs (filter-map
637                               (lambda (nv)
638                                 (and (not (member (first nv) matlab-builtin-consts))
639                                      (let ((v1 (canonicalize-expr/MATLAB (second nv))))
640                                        (list (first nv) v1))))
641                                   consts))
642                  (locals  (find-locals (map second const-defs))))
643             (for-each (lambda (def)
644                         (let ((n (first def)) (b (second def)))
645                           (pp indent+ ,(expr->string/MATLAB b n)))) const-defs)
646             (for-each (lambda (ep)
647                         (let* ((ep-name     (first ep))
648                                (ep-props    (second ep))
649                                (init-expr   (lookup-def 'initial ep-props))
650                                (temp-expr   (lookup-def 'temp-adj ep-props))
651                                (beta-expr   (lookup-def 'beta ep-props))
652                                (depth-expr  (lookup-def 'depth ep-props))
653                                (init-name   (matlab-name (s+ ep-name '-init)))
654                                (temp-name   (matlab-name (s+ ep-name '-temp-adj)))
655                                (beta-name   (matlab-name (s+ ep-name '-beta)))
656                                (depth-name  (matlab-name (s+ ep-name '-depth))))
657                           (if (or (not beta-expr) (not depth-expr) (not init-expr))
658                               (nemo:error 'nemo:matlab-translator 
659                                           ": ion pool " ep-name " requires initial value, depth and beta parameters"))
660                           (let ((temp-val  (and temp-expr (eval-const sys temp-expr)))
661                                 (init-val  (eval-const sys init-expr))
662                                 (beta-val  (eval-const sys beta-expr))
663                                 (depth-val (eval-const sys depth-expr)))
664                             (pp indent+ ,(expr->string/MATLAB init-val init-name))
665                             (pp indent+ ,(expr->string/MATLAB temp-val temp-name))
666                             (pp indent+ ,(expr->string/MATLAB beta-val beta-name))
667                             (pp indent+ ,(expr->string/MATLAB depth-val depth-name)))))
668                       epools))
669
670
671
672           (if (not (null? asgns))
673               (for-each (lambda (def)
674                           (let ((n (matlab-name (first def)) )
675                                 (b (second def)))
676                             (pp indent+ ,(expr->string/MATLAB b n)))) asgn-eq-defs))
677
678           (let* ((i-eqs (filter-map
679                          (lambda (ionch) 
680
681                            (let* ((label     (first ionch))
682                                   (n         (second ionch))
683                                   (subcomps  ((dis 'component-subcomps) sys n))
684                                   (acc       (lookup-def 'accumulating-substance subcomps))
685                                   (perm      (lookup-def 'permeating-substance subcomps))
686                                   (permqs    (and perm ((dis 'component-exports) sys (cid perm))))
687                                   (pore      (lookup-def 'pore subcomps))
688                                   (gate      (lookup-def 'gate subcomps))
689                                   (sts       (and gate ((dis 'component-exports) sys (cid gate)))))
690
691                              (cond ((and perm pore gate)
692                                     (case (cn perm)
693                                       ((non-specific)
694                                        (let* ((i     (matlab-name 'i))
695                                               (e     (car permqs))
696                                               (gmax  (car ((dis 'component-exports) sys (cid pore))))
697                                               (pwrs  (map (lambda (n) (state-power sys n)) sts))
698                                               (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
699                                               (gion  `(* ,gmax ,@sptms)))
700                                          (list i e gion)))
701                                       (else
702                                        (let* ((i     (matlab-name (s+ 'i (cn perm))))
703                                               (e     (matlab-name (s+ 'e (cn perm))))
704                                               (gmax  (car ((dis 'component-exports) sys (cid pore))))
705                                               (pwrs  (map (lambda (n) (state-power sys n)) sts))
706                                               (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
707                                               (gion  `(* ,gmax ,@sptms)))
708                                          (list i e gion)))))
709
710                                    ((and perm pore)
711                                     (case (cn perm)
712                                       ((non-specific)
713                                        (let* ((i     (matlab-name 'i))
714                                               (e     (car permqs))
715                                               (gmax  (car ((dis 'component-exports) sys (cid pore)))))
716                                          (list i e gmax)))
717                                       (else
718                                        (nemo:error 'nemo:matlab-translator ": invalid ion channel definition " label))))
719
720                                      ((and acc pore gate)
721                                       (let* ((i     (matlab-name (s+ 'i (cn acc))))
722                                              (gmax  (car ((dis 'component-exports) sys (cid pore))))
723                                              (pwrs  (map (lambda (n) (state-power sys n)) sts))
724                                              (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
725                                              (gion  `(* ,gmax ,@sptms)))
726                                         (list i #f gion)))
727                                      (else (nemo:error 'nemo:matlab-translator ": invalid ion channel definition " label))
728                                      )))
729                          ionchs))
730
731                  (i-bkts (bucket-partition (lambda (x y) (eq? (car x) (car y))) i-eqs))
732
733                  (i-eqs  (fold (lambda (b ax) 
734                                  (match b 
735                                         ((and ps ((i e gion) . rst)) 
736                                          (let* ((sum   (if e (sum (map (lambda (b) `(* ,(third b) (- v ,(second b)))) 
737                                                                        ps))
738                                                            (sum (map third ps))))
739                                                 (sum0  (rhsexpr sum))
740                                                 (sum1  (canonicalize-expr/MATLAB sum0)))
741                                            (cons (list i sum1) ax)))
742                                           
743                                         ((i e gion)
744                                          (let* ((expr0  (rhsexpr (if e `(* ,gion (- v ,e)) gion)))
745                                                 (expr1  (canonicalize-expr/MATLAB expr0)))
746                                              (cons (list i expr1) ax)))
747                                         
748                                         (else ax)))
749                                (list) i-bkts))
750
751                  (state-eq-defs    (reverse (poset->state-eq-defs poset sys)))
752
753                  (stcomp-eq-defs   (poset->stcomp-eq-defs poset sys))
754
755                  (pool-eq-defs
756                   (map (lambda (ep)
757                          (let* ((ep-name     (first ep))
758                                 (pool-ion    (assoc ep-name pool-ions))
759                                 (i-name      (second pool-ion))
760                                 (init-name   (matlab-name (s+ ep-name '-init)))
761                                 (temp-name   (matlab-name (s+ ep-name '-temp-adj)))
762                                 (beta-name   (matlab-name (s+ ep-name '-beta)))
763                                 (depth-name  (matlab-name (s+ ep-name '-depth)))
764                                 (rhs         `(let ((F 96485.0))
765                                                 (- (/ (neg ,i-name) (* 2 F ,init-name ,depth-name)) 
766                                                    (* ,beta-name ,ep-name . 
767                                                       ,(if temp-name (list temp-name) (list)))))))
768                            `(,(s+ ep-name)  ,rhs)))
769                        epools))
770
771                  (rate-eq-defs (append pool-eq-defs state-eq-defs))
772                 
773                  (state-index-map  (let ((acc (fold (lambda (def ax)
774                                                       (let ((st-name (first def)))
775                                                         (list (+ 1 (first ax)) 
776                                                               (cons `(,st-name ,(first ax)) (second ax)))))
777                                                     (list 0 (list)) rate-eq-defs)))
778                                    (second acc)))
779                  )
780
781
782             (pp indent+ (dy = zeros(length(y),1)))
783
784             (for-each (lambda (def)
785                         (let ((n (first def)) )
786                           (pp indent+ (,n = ,(s+ "y(" (lookup-def n state-index-map) ")"))))) rate-eq-defs)
787
788             (for-each (lambda (def)
789                         (let ((n (first def)) (b (second def)))
790                           (pp indent+ ,(expr->string/MATLAB b n)))) stcomp-eq-defs)
791
792             (for-each (lambda (def)
793                         (let ((n (s+ "dy(" (lookup-def (first def) state-index-map) ")")) (b (second def)))
794                           (pp indent+ ,(expr->string/MATLAB b n)))) rate-eq-defs)
795
796             (for-each (lambda (pool-ion)
797                         (pp indent+ (,(third pool-ion) = ,(first pool-ion)))) pool-ions)
798
799             (for-each (lambda (def) 
800                         (pp indent+ ,(expr->string/MATLAB (second def) (first def)))) i-eqs)
801             
802             ))))))))
Note: See TracBrowser for help on using the repository browser.