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

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

Added skeleton for matlab code generator.

File size: 35.2 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 (s\ p . lst)  (string-intersperse (map ->string lst) p))
102(define (sl\ p lst)   (string-intersperse (map ->string lst) p))
103(define nl "\n")
104
105(define (spaces n)  (list->string (list-tabulate n (lambda (x) #\space))))
106
107(define (ppf indent . lst)
108  (let ((sp (spaces indent)))
109    (for-each (lambda (x)
110                (and x (match x 
111                              ((i . x1) (if (and (number? i) (positive? i))
112                                            (for-each (lambda (x) (ppf (+ indent i) x)) x1)
113                                            (print sp (sw+ x))))
114                              (else   (print sp (if (list? x) (sw+ x) x))))))
115              lst)))
116
117
118(define-syntax pp
119  (syntax-rules ()
120    ((pp indent val ...) (ppf indent (quasiquote val) ...))))
121
122(define (ifthen/MATLAB c e1 e2)
123  (doc:nest 
124   2 (doc:connect
125      (doc:connect (doc:group (doc:connect (doc:text "if") c))
126                   (doc:connect (doc:nest 2 e1)
127                                (doc:nest 2 (doc:connect (doc:text "else") e2))))
128      (doc:text "endif"))))
129
130
131(define (letblk/MATLAB e1 e2)
132  (cond ((equal? e1 (doc:empty)) (doc:group (doc:nest 2 e2)))
133        ((equal? e2 (doc:empty)) (doc:group (doc:nest 2 e1)))
134        (else (doc:connect (doc:group (doc:concat (list (doc:nest 2 e1) (doc:text ";"))))
135                           (doc:group (doc:nest 2 e2))))))
136       
137
138(define group/MATLAB   (doc:block 2 (doc:text "(") (doc:text ")")))
139(define block/MATLAB   (doc:block 2 (doc:empty) (doc:empty)))
140
141(define (format-op/MATLAB indent op args)
142  (let ((op1 (doc:text (->string op))))
143    (if (null? args) op1
144        (match args
145               ((x)      (doc:concat (list op1 x)))
146               ((x y)    (doc:concat (intersperse (list x op1 y) (doc:space))))
147               ((x y z)  (doc:concat (intersperse (list x op1 y op1 z) (doc:space))))
148               (lst      (let* ((n   (length lst))
149                                (n/2 (inexact->exact (round (/ n 2)))))
150                           (doc:concat 
151                            (intersperse 
152                             (list (format-op/MATLAB indent op (take lst n/2 )) op1 
153                                   (format-op/MATLAB indent op (drop lst n/2 )))
154                             (doc:space)))))))))
155
156(define (format-fncall/MATLAB indent op args)
157  (let ((op1 (doc:text (->string op))))
158    (doc:cons op1 (group/MATLAB ((doc:list indent identity (lambda () (doc:text ", "))) args)))))
159
160(define matlab-builtin-consts
161  `())
162
163(define matlab-ops
164  `(+ - * / > < <= >= = ^))
165
166(define builtin-fns
167  `(+ - * / pow neg abs atan asin acos sin cos exp ln
168      sqrt tan cosh sinh tanh hypot gamma lgamma log10 log2 log1p ldexp cube
169      > < <= >= = and or round ceiling floor max min
170      fpvector-ref))
171
172(define (sum lst)
173  (if (null? lst) lst
174      (match lst
175             ((x)   x)
176             ((x y) `(+ ,x ,y))
177             ((x y . rest) `(+ (+ ,x ,y) ,(sum rest)))
178             ((x . rest) `(+ ,x ,(sum rest))))))
179
180
181(define (subst-term t subst k)
182  (match t
183         (('if c t e)
184          `(if ,(k c subst) ,(k t subst) ,(k e subst)))
185         (('let bs e)
186          (k `(let ,(map (lambda (b) `(,(car b) ,(k (cadr b) subst))) bs) ,(k e subst)) subst))
187         ((f . es)
188          (cons (k f subst) (map (lambda (e) (k e subst)) es)))
189         ((? symbol? )  (lookup-def t subst t))
190         ((? atom? ) t)))
191
192(define (binding? t) (and (list? t) (eq? 'let (car t)) (cdr t)))
193
194(define (bind ks vs e) `(let ,(zip ks vs) ,e))
195
196(define (name-normalize expr)
197  (match expr 
198         (('if c t e)  `(if ,(name-normalize c) ,(name-normalize t) ,(name-normalize e)))
199         (('let bs e)
200          `(let ,(map (lambda (b) `(,(car b) ,(name-normalize (cadr b)))) bs) ,(name-normalize e)))
201         ((f . es) 
202          (cons f (map name-normalize es)))
203         ((? symbol? ) (matlab-name expr))
204         ((? atom? ) expr)))
205
206(define (if-convert expr)
207  (match expr 
208         (('if c t e) 
209          (let ((r (gensym "if")))
210            `(let ((,r (if ,(if-convert c) ,(if-convert t) ,(if-convert e)))) 
211               ,r)))
212         (('let bs e)
213          `(let ,(map (lambda (b) `(,(car b) ,(if-convert (cadr b)))) bs) ,(if-convert e)))
214         ((f . es)
215          (cons f (map if-convert es)))
216         ((? atom? ) expr)))
217
218         
219(define (let-enum expr ax)
220  (match expr
221         (('let ((x ('if c t e))) y)
222          (let ((ax (fold let-enum ax (list c ))))
223            (if (eq? x y)  (append ax (list (list x `(if ,c ,t ,e)))) ax)))
224
225         (('let bnds body)  (let-enum body (append ax bnds)))
226
227         (('if c t e)  (let-enum ax c))
228
229         ((f . es)  (fold let-enum ax es))
230
231         (else ax)))
232
233
234(define (let-elim expr)
235  (match expr
236         (('let ((x ('if c t e))) y)
237          (if (eq? x y)  y expr))
238
239         (('let bnds body) (let-elim body))
240
241         (('if c t e)  `(if ,(let-elim c) ,(let-lift t) ,(let-lift e)))
242
243         ((f . es)  `(,f . ,(map let-elim es)))
244
245         (else expr)))
246 
247
248(define (let-lift expr)
249  (let ((bnds (let-enum expr (list))))
250    (if (null? bnds) expr
251        `(let ,(map (lambda (b) (list (car b) (let-elim (cadr b)))) bnds) ,(let-elim expr)))))
252
253(define (canonicalize-expr/MATLAB expr)
254  (let ((subst-convert  (subst-driver (lambda (x) (and (symbol? x) x)) binding? identity bind subst-term)))
255    (let* ((expr1 (if-convert expr))
256           (expr2 (subst-convert expr1 subst-empty))
257           (expr3 (let-lift expr2))
258           (expr4 (name-normalize expr3)))
259      expr4)))
260
261
262(define (format-expr/MATLAB indent expr . rest) 
263  (let-optionals rest ((rv #f))
264   (let ((indent+ (+ 2 indent)))
265    (match expr
266       (('let bindings body)
267        (letblk/MATLAB
268         (fold-right 
269           (lambda (x ax)
270             (letblk/MATLAB
271              (match (second x)
272                     (('if c t e)
273                      (ifthen/MATLAB
274                       (group/MATLAB (format-expr/MATLAB indent c))
275                       (block/MATLAB (format-expr/MATLAB indent t (first x)))
276                       (block/MATLAB (format-expr/MATLAB indent e (first x)))))
277                     (else
278                      (format-op/MATLAB indent+ " = "
279                                       (list (format-expr/MATLAB indent (first x) )
280                                             (format-expr/MATLAB indent (second x))))))
281              ax))
282           (doc:empty) bindings)
283         (let ((body1 (doc:nest indent (format-expr/MATLAB indent body))))
284           (if rv  (format-op/MATLAB indent " = " (list (format-expr/MATLAB indent+ rv ) body1))
285               body1))))
286       
287       (('if . rest) (error 'format-expr/MATLAB "invalid if statement " expr))
288
289       ((op . rest) 
290       (let ((op (case op ((pow)  '^)  (else op))))
291         (let ((fe
292                (if (member op matlab-ops)
293                    (let ((mdiv?  (any (lambda (x) (match x (('* . _) #t) (('/ . _) #t) (else #f))) rest))
294                          (mul?   (any (lambda (x) (match x (('* . _) #t) (else #f))) rest))
295                          (plmin? (any (lambda (x) (match x (('+ . _) #t) (('- . _) #t) (else #f))) rest)))
296                      (case op
297                        ((/) 
298                         (format-op/MATLAB indent op 
299                                          (map (lambda (x) 
300                                                 (let ((fx (format-expr/MATLAB indent+ x)))
301                                                   (if (or (symbol? x) (number? x)) fx
302                                                       (if (or mul? plmin?) (group/MATLAB fx) fx)))) rest)))
303                        ((*) 
304                         (format-op/MATLAB indent op 
305                                          (map (lambda (x) 
306                                                 (let ((fx (format-expr/MATLAB indent+ x)))
307                                                   (if (or (symbol? x) (number? x)) fx
308                                                       (if plmin? (group/MATLAB fx) fx)))) rest)))
309                       
310                        ((^) 
311                         (format-op/MATLAB indent op 
312                                          (map (lambda (x) 
313                                                 (let ((fx (format-expr/MATLAB indent+ x)))
314                                                   (if (or (symbol? x)  (number? x)) fx
315                                                       (if (or mdiv? plmin?) (group/MATLAB fx) fx)))) rest)))
316                       
317                        (else
318                         (format-op/MATLAB indent op 
319                                          (map (lambda (x) 
320                                                 (let ((fx (format-expr/MATLAB indent+ x))) fx)) rest)))))
321                   
322                    (let ((op (case op ((neg) '-) (else op))))
323                      (format-fncall/MATLAB indent op (map (lambda (x) (format-expr/MATLAB indent+ x)) rest))))))
324           (if rv 
325               (format-op/MATLAB indent " = " (list (format-expr/MATLAB indent+ rv ) fe))
326               fe))))
327     
328      (else  (let ((fe (doc:text (->string expr))))
329               (if rv 
330                   (format-op/MATLAB indent " = " (list (format-expr/MATLAB indent+ rv ) fe))
331                   fe)))))))
332               
333
334         
335(define (expr->string/MATLAB x . rest)
336  (let-optionals rest ((rv #f) (width 72))
337    (sdoc->string (doc:format width (format-expr/MATLAB 2 x rv)))))
338 
339
340(define (format-lineq/MATLAB indent expr . rest) 
341  (let-optionals rest ((rv #f))
342   (let ((indent+ (+ 2 indent)))
343    (match expr
344       (('let bindings body)
345        (letblk/MATLAB
346         (fold-right 
347           (lambda (x ax)
348             (letblk/MATLAB
349              (match (second x)
350                     (('if c t e)
351                      (ifthen/MATLAB
352                       (group/MATLAB (format-lineq/MATLAB indent c))
353                       (block/MATLAB (format-lineq/MATLAB indent t (first x)))
354                       (block/MATLAB (format-lineq/MATLAB indent e (first x)))))
355                     (else
356                      (format-op/MATLAB indent+ " = "
357                                       (list (format-lineq/MATLAB indent (first x) )
358                                             (format-lineq/MATLAB indent (second x))))))
359              ax))
360           (doc:empty) bindings)
361         (let ((body1 (doc:nest indent (format-lineq/MATLAB indent body))))
362           (if rv  (format-op/MATLAB indent " = " (list (format-lineq/MATLAB indent+ rv ) body1))
363               body1))))
364       
365       (('if . rest) (error 'format-lineq/MATLAB "invalid if statement " expr))
366
367       ((op . rest) 
368        (let ((op (case op ((pow)  '^) ((abs) 'fabs) (else op))))
369          (let ((fe
370                (if (member op matlab-ops)
371                    (let ((mdiv?  (any (lambda (x) (match x (('* . _) #t) (('/ . _) #t) (else #f))) rest))
372                          (mul?   (any (lambda (x) (match x (('* . _) #t) (else #f))) rest))
373                          (plmin? (any (lambda (x) (match x (('+ . _) #t) (('- . _) #t) (else #f))) rest)))
374                      (case op
375                        ((/) 
376                         (format-op/MATLAB indent op 
377                                          (map (lambda (x) 
378                                                 (let ((fx (format-lineq/MATLAB indent+ x)))
379                                                   (if (or (symbol? x) (number? x)) fx
380                                                       (if (or mul? plmin?) (group/MATLAB fx) fx)))) rest)))
381                        ((*) 
382                         (format-op/MATLAB indent op 
383                                          (map (lambda (x) 
384                                                 (let ((fx (format-lineq/MATLAB indent+ x)))
385                                                   (if (or (symbol? x) (number? x)) fx
386                                                       (if plmin? (group/MATLAB fx) fx)))) rest)))
387                       
388                        ((^) 
389                         (format-op/MATLAB indent op 
390                                          (map (lambda (x) 
391                                                 (let ((fx (format-lineq/MATLAB indent+ x)))
392                                                   (if (or (symbol? x)  (number? x)) fx
393                                                       (if (or mdiv? plmin?) (group/MATLAB fx) fx)))) rest)))
394                       
395                        (else
396                         (format-op/MATLAB indent op 
397                                          (map (lambda (x) 
398                                                 (let ((fx (format-lineq/MATLAB indent+ x))) fx)) rest)))))
399                   
400                    (case op
401                      ((neg) (format-op/MATLAB indent '* (map (lambda (x) (format-lineq/MATLAB indent+ x)) 
402                                                             (cons "(-1)" rest))))
403                      (else  (format-fncall/MATLAB indent op (map (lambda (x) (format-lineq/MATLAB indent+ x)) 
404                                                                 rest)))))))
405
406           (if rv (format-op/MATLAB indent " = " (list (format-lineq/MATLAB indent+ rv ) fe)) fe))))
407     
408      (else  (let ((fe (doc:text (->string expr))))
409               (if rv 
410                   (format-op/MATLAB indent " = " (list (format-lineq/MATLAB indent+ rv ) fe))
411                   fe)))))))
412               
413
414         
415(define (lineq->string/MATLAB x val . rest)
416  (let-optionals rest ((width 72))
417    (s+ "~ " (sdoc->string (doc:format width (format-lineq/MATLAB 2 x #f))) 
418        " = " (number->string val))))
419 
420
421(define (make-define-fn table? min-v max-v with depend)
422  (lambda (indent n proc)
423    (let ((lst (procedure-data proc))
424          (indent+ (+ 2 indent)))
425      (let ((retval   (matlab-name (gensym "retval")))
426            (rt       (lookup-def 'rt lst))
427            (formals  (lookup-def 'formals lst))
428            (vars     (lookup-def 'vars lst))
429            (body     (lookup-def 'body lst)))
430        (pp indent ,nl (function ,retval = ,n (,(sl\ ", " vars)) ))
431        (let* ((body1 (canonicalize-expr/MATLAB (rhsexpr body)))
432               (lbs   (enum-bnds body1 (list))))
433;         (if (and table? min-v max-v with)
434;             (match vars
435;                    (('v)  (pp indent+ (TABLE ,@(if depend `(DEPEND ,depend) `(""))
436;                                              FROM ,min-v TO ,max-v WITH ,with)))
437;                    (else  (void))))
438          (pp indent+ ,(expr->string/MATLAB body1 retval))
439          (pp indent endfunction))
440          ))))
441
442(define (define-state indent n)
443  (pp indent (,n)))
444
445
446(define (state-eqs n initial open transitions power)
447  (let* ((subst-convert  (subst-driver (lambda (x) (and (symbol? x) x)) binding? identity bind subst-term))
448         (g          (make-digraph n (string-append (->string n) " transitions graph")))
449         (add-node!  (g 'add-node!))
450         (add-edge!  (g 'add-edge!))
451         (out-edges  (g 'out-edges))
452         (in-edges   (g 'in-edges))
453         (node-info  (g 'node-info))
454         (node-list  (let loop ((lst (list)) (tlst transitions))
455                       (if (null? tlst)  (delete-duplicates lst eq?)
456                           (match (car tlst) 
457                                  (('-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr)
458                                   (loop (cons* s0 s1 lst) (cdr tlst)))
459                                  (((and (? symbol?) s0) '-> (and (? symbol? s1)) rate-expr)
460                                   (loop (cons* s0 s1 lst) (cdr tlst)))
461                                  (('<-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr1 rate-expr2)
462                                   (loop (cons* s0 s1 lst) (cdr tlst)))
463                                  (((and (? symbol?) s0) 'M-> (and (? symbol? s1)) rate-expr1 rate-expr2)
464                                   (loop (cons* s0 s1 lst) (cdr tlst)))
465                                  (else
466                                   (nemo:error 'nemo:matlab-state-eqs ": invalid transition equation " 
467                                                  (car tlst) " in state complex " n))
468                                  (else (loop lst (cdr tlst)))))))
469         (node-ids      (list-tabulate (length node-list) identity))
470         (name->id-map  (zip node-list node-ids))
471         (node-subs     (fold (lambda (s ax) (subst-extend s (matlab-state-name n s) ax)) subst-empty node-list)))
472    ;; insert state nodes in the dependency graph
473    (for-each (lambda (i n) (add-node! i n)) node-ids node-list)
474    (let* ((nodes  ((g 'nodes)))
475           (snode   (find (lambda (s) (not (eq? (second s) open))) nodes))
476           (snex   `(- 1 ,(sum (filter-map (lambda (s) (and (not (= (first s) (first snode))) (second s))) nodes))))
477           (add-tredge (lambda (s0 s1 rexpr1 rexpr2)
478                         (let ((i   (car (alist-ref s0 name->id-map)))
479                               (j   (car (alist-ref s1 name->id-map)))
480                               (x0  (if (eq? s0 (second snode)) snex s0))
481                               (x1  (if (eq? s1 (second snode)) snex s1)))
482                           (add-edge! (list i j `(* ,(subst-convert x0 node-subs) 
483                                                    ,(subst-convert rexpr1 node-subs))))
484                           (if rexpr2 (add-edge! (list j i `(* ,(subst-convert x1 node-subs) 
485                                                               ,(subst-convert rexpr2 node-subs)))))))))
486      ;; create rate edges in the graph
487      (for-each (lambda (e) 
488                  (match e
489                         (('-> s0 s1 rexpr)  (add-tredge s0 s1 rexpr #f))
490                         ((s0 '-> s1 rexpr)  (add-tredge s0 s1 rexpr #f))
491                         (('<-> s0 s1 rexpr1 rexpr2)  (add-tredge s0 s1 rexpr1 rexpr2))
492                         ((s0 '<-> s1 rexpr1 rexpr2)  (add-tredge s0 s1 rexpr1 rexpr2))
493                         ))
494                transitions)
495
496      ;; generate differential equations for each state in the transitions system
497      (let ((eqs    (fold (lambda (s ax) 
498                            (if (= (first snode) (first s) ) ax
499                                (let* ((out   (out-edges (first s)))
500                                       (in    (in-edges (first s)))
501                                       (open? (eq? (second s) open))
502                                       (name  (matlab-name (lookup-def (second s) node-subs))))
503                                  (let* ((rhs1  (cond ((and (not (null? out)) (not (null? in)))
504                                                       `(+ (neg ,(sum (map third out)))
505                                                           ,(sum (map third in))))
506                                                      ((and (not (null? out)) (null? in))
507                                                       `(neg ,(sum (map third out))))
508                                                      ((and (null? out) (not (null? in)))
509                                                       (sum (map third in)))))
510                                         (fbody  (rhsexpr rhs1))
511                                         (fbody1 (canonicalize-expr/MATLAB fbody)))
512                                    (cons (list (s+ name "'") fbody1) ax)))))
513                          (list) nodes)))
514        eqs))))
515           
516       
517
518(define (kstate-eqs n initial open transitions power)
519  (let* ((subst-convert  (subst-driver (lambda (x) (and (symbol? x) x)) binding? identity bind subst-term))
520         (state-list     (let loop ((lst (list)) (tlst transitions))
521                           (if (null? tlst)  (delete-duplicates lst eq?)
522                               (match (car tlst) 
523                                      (('-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr)
524                                       (loop (cons* s0 s1 lst) (cdr tlst)))
525                                      (((and (? symbol?) s0) '-> (and (? symbol? s1)) rate-expr)
526                                       (loop (cons* s0 s1 lst) (cdr tlst)))
527                                      (('<-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr1 rate-expr2)
528                                       (loop (cons* s0 s1 lst) (cdr tlst)))
529                                      (((and (? symbol?) s0) 'M-> (and (? symbol? s1)) rate-expr1 rate-expr2)
530                                       (loop (cons* s0 s1 lst) (cdr tlst)))
531                                      (else
532                                       (nemo:error 'nemo:matlab-kstate-eqs ": invalid transition equation " 
533                                                   (car tlst) " in state complex " n))
534                                      (else (loop lst (cdr tlst)))))))
535         (state-subs     (fold (lambda (s ax) (subst-extend s (matlab-state-name n s) ax)) subst-empty state-list)))
536    ;; generate kinetic equations for each edge in the transitions system
537    (map
538     (lambda (e) 
539       (match e
540              (('-> s0 s1 rexpr)
541               (let ((i  (lookup-def s0 state-subs))
542                     (j  (lookup-def s1 state-subs)))
543                 `(-> ,i ,j ,(subst-convert rexpr state-subs))))
544             
545              ((s0 '-> s1 rexpr)
546               (let ((i  (lookup-def s0 state-subs))
547                     (j  (lookup-def s1 state-subs)))
548                 `(-> ,i ,j ,(subst-convert rexpr state-subs))))
549             
550              (('<-> s0 s1 rexpr1 rexpr2)
551               (let ((i  (lookup-def s0 state-subs))
552                     (j  (lookup-def s1 state-subs)))
553                 `(<-> ,i ,j ,(subst-convert rexpr1 state-subs) ,(subst-convert rexpr2 state-subs))))
554             
555              ((s0 '<-> s1 rexpr1 rexpr2)
556               (let ((i  (lookup-def s0 state-subs))
557                     (j  (lookup-def s1 state-subs)))
558                 `(<-> ,i ,j ,(subst-convert rexpr1 state-subs) ,(subst-convert rexpr2 state-subs))))
559             
560                 
561              (else (nemo:error 'nemo:matlab-kstate-eqs ": invalid transition equation " 
562                                e " in state complex " n))))
563     transitions)))
564       
565
566
567(define (state-init n init)
568  (let* ((init  (rhsexpr init))
569         (init1 (canonicalize-expr/MATLAB init)))
570    (list (matlab-name n) init1)))
571
572
573(define (state-init-eq n transitions init)
574  (let* ((subst-convert  (subst-driver (lambda (x) (and (symbol? x) x)) binding? identity bind subst-term))
575         (state-list     (let loop ((lst (list)) (tlst transitions))
576                           (if (null? tlst)  (delete-duplicates lst eq?)
577                               (match (car tlst) 
578                                      (('-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr)
579                                       (loop (cons* s0 s1 lst) (cdr tlst)))
580                                      (((and (? symbol?) s0) '-> (and (? symbol? s1)) rate-expr)
581                                       (loop (cons* s0 s1 lst) (cdr tlst)))
582                                      (('<-> (and (? symbol?) s0) (and (? symbol?) s1) rate-expr1 rate-expr2)
583                                       (loop (cons* s0 s1 lst) (cdr tlst)))
584                                      (((and (? symbol?) s0) 'M-> (and (? symbol? s1)) rate-expr1 rate-expr2)
585                                       (loop (cons* s0 s1 lst) (cdr tlst)))
586                                      (else
587                                       (nemo:error 'nemo:state-init-eq ": invalid transition equation " 
588                                                   (car tlst) " in state complex " n))
589                                      (else (loop lst (cdr tlst)))))))
590         (state-subs     (fold (lambda (s ax) (subst-extend s (matlab-state-name n s) ax)) subst-empty state-list))
591         (init1          (map (lambda (lineq) (match lineq ((i '= . expr) `(,i = . ,(subst-convert expr state-subs)))))
592                              init)))
593    (list (matlab-name n) init1)))
594
595
596(define (asgn-eq n rhs)
597  (let* ((fbody   (rhsexpr rhs))
598         (fbody1  (canonicalize-expr/MATLAB fbody)))
599    (list (matlab-name n) fbody1)))
600
601
602(define (stcomp-eq n open transitions)
603  (list (matlab-name n) (matlab-name (matlab-state-name n open))))
604
605
606(define (poset->asgn-eq-defs poset sys)
607  (fold-right
608   (lambda (lst ax)
609     (fold  (lambda (x ax) 
610              (match-let (((i . n)  x))
611                         (let ((en (environment-ref sys n)))
612                           (if (nemo:quantity? en)
613                               (cases nemo:quantity en
614                                      (ASGN  (name value rhs) (cons (asgn-eq name rhs) ax))
615                                      (else  ax))
616                               ax))))
617            ax lst))
618   (list) poset))
619
620
621(define (poset->state-eq-defs poset sys kinetic)
622  (fold-right
623   (lambda (lst ax)
624     (fold  (lambda (x ax) 
625              (match-let (((i . n)  x))
626                         (let ((en (environment-ref sys n)))
627                           (if (and (not (member n kinetic)) (nemo:quantity? en))
628                               (cases nemo:quantity en
629                                      (TSCOMP  (name initial open transitions power) 
630                                               (append (state-eqs name initial open transitions power) ax))
631                                      (else  ax))
632                               ax))))
633            ax lst))
634   (list) poset))
635
636
637(define (poset->kstate-eq-defs poset sys kinetic)
638  (fold-right
639   (lambda (lst ax)
640     (fold  (lambda (x ax) 
641              (match-let (((i . n)  x))
642                         (let ((en (environment-ref sys n)))
643                           (if (and (member n kinetic) (nemo:quantity? en))
644                               (cases nemo:quantity en
645                                      (TSCOMP  (name initial open transitions power) 
646                                               (append (kstate-eqs name initial open transitions power) ax))
647                                      (else  ax))
648                               ax))))
649            ax lst))
650   (list) poset))
651
652
653(define (poset->stcomp-eq-defs poset sys)
654  (fold-right
655   (lambda (lst ax)
656     (fold  (lambda (x ax) 
657              (match-let (((i . n)  x))
658                         (let ((en (environment-ref sys n)))
659                           (if (nemo:quantity? en)
660                               (cases nemo:quantity en
661                                      (TSCOMP  (name initial open transitions power) 
662                                               (cons (stcomp-eq name open transitions) ax))
663                                      (else  ax))
664                               ax))))
665            ax lst))
666   (list) poset))
667
668
669(define (poset->state-init-defs poset sys)
670  (fold-right
671   (lambda (lst ax)
672     (fold  (lambda (x ax) 
673              (match-let (((i . n)  x))
674                         (let ((en (environment-ref sys n)))
675                           (if (nemo:quantity? en)
676                               (cases nemo:quantity en
677                                      (TSCOMP  (name initial open transitions power)
678                                               (if (nemo:rhs? initial)
679                                                   (cons* (state-init name initial) 
680                                                          (state-init (matlab-state-name name open) name) ax) 
681                                                   ax))
682                                      (else  ax))
683                               ax))))
684            ax lst))
685   (list) poset))
686
687
688(define (poset->state-init-eq-defs poset sys)
689  (fold-right
690   (lambda (lst ax)
691     (fold  (lambda (x ax) 
692              (match-let (((i . n)  x))
693                         (let ((en (environment-ref sys n)))
694                           (if (nemo:quantity? en)
695                               (cases nemo:quantity en
696                                      (TSCOMP (name initial open transitions power)
697                                              (if (and (list? initial) (every nemo:lineq? initial))
698                                                  (cons (state-init-eq name transitions initial) ax) 
699                                                  ax))
700                                      (else  ax))
701                               ax))))
702            ax lst))
703   (list) poset))
704
705
706(define (find-locals defs)
707  (concatenate (map (lambda (def) (match def (('let bnds _) (map first bnds)) (else (list)))) defs)))
708
709
710(define (state-power sys n)
711  (let ((en (environment-ref sys n)))
712    (if (nemo:quantity? en)
713        (cases nemo:quantity en
714               (TSCOMP  (name initial open transitions power)  power)
715               (else  #f))  #f)))
716
717
718(define (bucket-partition p lst)
719  (let loop ((lst lst) (ax (list)))
720    (if (null? lst) ax
721        (let ((x (car lst)))
722          (let bkt-loop ((old-bkts ax) (new-bkts (list)))
723            (if (null? old-bkts) (loop (cdr lst) (cons (list x) new-bkts))
724                (if (p x (caar old-bkts ))
725                    (loop (cdr lst) (append (cdr old-bkts) (cons (cons x (car old-bkts)) new-bkts)))
726                    (bkt-loop (cdr old-bkts) (cons (car old-bkts) new-bkts)))))))))
727
728
729(define (collect-epools sys)
730   (match-let ((($ nemo:quantity 'DISPATCH  dis) (environment-ref sys (nemo-intern 'dispatch))))
731      (let recur ((comp-name (nemo-intern 'toplevel)) (ax (list)))
732        (let* ((comp-symbols   ((dis 'component-symbols)   sys comp-name))
733               (subcomps       ((dis 'component-subcomps)  sys comp-name)))
734          (fold recur 
735                (fold (lambda (sym ax)
736                        (let ((en (environment-ref sys sym)))
737                          (match en
738                                 ((or (('decaying 'pool)  ('name (? symbol? ion)) . alst)
739                                      (('decaying-pool)   ('name (? symbol? ion)) . alst))
740                                  (cons (list ion alst) ax))
741                                 (else ax)))) ax comp-symbols)
742                (map third subcomps))))))
743
744
745(define (nemo:matlab-translator sys . rest)
746  (define (cid x)  (second x))
747  (define (cn x)   (first x))
748  (let-optionals rest ((method 'cnexp) (table? #f) (min-v -100) (max-v 100) (step 0.5) 
749                       (depend #f)  (kinetic (list)) )
750  (match-let ((($ nemo:quantity 'DISPATCH  dis) (environment-ref sys (nemo-intern 'dispatch))))
751    (let ((imports  ((dis 'imports)  sys))
752          (exports  ((dis 'exports)  sys)))
753      (let* ((indent      0)
754             (indent+     (+ 2 indent ))
755             (eval-const  (dis 'eval-const))
756             (sysname     (matlab-name ((dis 'sysname) sys)))
757             (deps*       ((dis 'depgraph*) sys))
758             (consts      ((dis 'consts)  sys))
759             (asgns       ((dis 'asgns)   sys))
760             (states      ((dis 'states)  sys))
761             (stcomps     ((dis 'stcomps) sys))
762             (defuns      ((dis 'defuns)  sys))
763             (components  ((dis 'components) sys))
764             (ionchs      (filter-map (match-lambda ((name 'ion-channel id) (list name id)) (else #f)) components))
765             (epools      (collect-epools sys)))
766
767
768        (match-let (((state-list asgn-list g) deps*))
769         (let* ((poset          (vector->list ((dis 'depgraph->bfs-dist-poset) g)))
770                (asgn-eq-defs   (poset->asgn-eq-defs poset sys))
771                (perm-ions (fold (lambda (ionch ax) 
772                                    (let* ((subcomps ((dis 'component-subcomps) sys (cid ionch)))
773                                           (perm      (lookup-def 'permeating-substance subcomps)))
774                                      (if perm 
775                                          (case (cn perm)
776                                            ((non-specific)   
777                                             (let* ((erev (car ((dis 'component-exports) sys (cid perm))))
778                                                    (i    (matlab-name 'i))
779                                                    (e    (matlab-name 'e)))
780                                               (cons `(,(cn perm) ,i ,e ,erev) ax)))
781                                            (else (let* ((erev (car ((dis 'component-exports) sys (cid perm))))
782                                                         (i    (matlab-name (s+ 'i (cn perm))))
783                                                         (e    (matlab-name (s+ 'e (cn perm)))))
784                                                    (cons `(,(cn perm) ,i ,e ,erev) ax))))
785                                          ax)))
786                                  (list) ionchs))
787               (acc-ions (fold (lambda (ionch ax) 
788                                  (let* ((subcomps ((dis 'component-subcomps) sys (cid ionch)))
789                                         (acc   (lookup-def 'accumulating-substance subcomps))
790                                         (i     (and acc (matlab-name (s+ 'i (cn acc)))))
791                                         (in    (and acc (matlab-name (s+ (cn acc) 'i))))
792                                         (out   (and acc (matlab-name (s+ (cn acc) 'o)))))
793                                    (if acc  (cons `(,(cn acc) ,i ,in ,out) ax) ax)))
794                                (list) ionchs))
795               (pool-ions (map (lambda (ep)
796                                  (let ((ion (car ep)))
797                                    `(,(matlab-name ion) ,(matlab-name (s+ 'i ion)) ,(matlab-name (s+ ion 'i)))))
798                               epools))
799               (has-kinetic? (or (not (null? (filter (lambda (x) (member (car x) kinetic)) states)))))
800               (has-ode?     (or (not (null? (filter (lambda (x) (not (member (car x) kinetic))) states)))
801                                 (not (null? pool-ions)))))
802
803
804           (for-each
805            (lambda (a)
806              (let ((acc-ion   (car a)))
807                (if (assoc acc-ion perm-ions)
808                    (nemo:error 'nemo:matlab-translator 
809                                ": ion species " acc-ion " cannot be declared as both accumulating and permeating"))))
810            acc-ions)
811
812           (for-each
813            (lambda (p)
814              (let ((pool-ion  (car p)))
815                (if (assoc pool-ion perm-ions)
816                    (nemo:error 'nemo:matlab-translator 
817                                ": ion species " pool-ion " cannot be declared as both pool and permeating"))))
818            pool-ions)
819
820           (let* ((with      (inexact->exact (round (/ (abs (- max-v min-v)) step))))
821                  (define-fn (make-define-fn table? min-v max-v with depend)))
822             (for-each (lambda (fndef) 
823                         (if (not (member (car fndef) builtin-fns))
824                             (apply define-fn (cons indent fndef)))) 
825                       defuns))
826
827           (pp indent ,nl (dy = ,sysname (,(sl\ ", " `(dy t)))))
828           
829           (if (not (null? exports)) (pp indent+ (global ,(sl\ ", " exports))))
830
831           (let* ((const-defs (filter-map
832                               (lambda (nv)
833                                 (and (not (member (first nv) matlab-builtin-consts))
834                                      (let ((v1 (canonicalize-expr/MATLAB (second nv))))
835                                        (list (first nv) v1))))
836                                   consts))
837                  (locals  (find-locals (map second const-defs))))
838             (for-each (lambda (def)
839                         (let ((n (first def)) (b (second def)))
840                           (pp indent+ ,(expr->string/MATLAB b n)))) const-defs)
841             (for-each (lambda (ep)
842                         (let* ((ep-name     (first ep))
843                                (ep-props    (second ep))
844                                (init-expr   (lookup-def 'initial ep-props))
845                                (temp-expr   (lookup-def 'temp-adj ep-props))
846                                (beta-expr   (lookup-def 'beta ep-props))
847                                (depth-expr  (lookup-def 'depth ep-props))
848                                (init-name   (matlab-name (s+ ep-name '-init)))
849                                (temp-name   (matlab-name (s+ ep-name '-temp-adj)))
850                                (beta-name   (matlab-name (s+ ep-name '-beta)))
851                                (depth-name  (matlab-name (s+ ep-name '-depth))))
852                           (if (or (not beta-expr) (not depth-expr) (not init-expr))
853                               (nemo:error 'nemo:matlab-translator 
854                                           ": ion pool " ep-name " requires initial value, depth and beta parameters"))
855                           (let ((temp-val  (and temp-expr (eval-const sys temp-expr)))
856                                 (init-val  (eval-const sys init-expr))
857                                 (beta-val  (eval-const sys beta-expr))
858                                 (depth-val (eval-const sys depth-expr)))
859                             (pp indent+ ,(expr->string/MATLAB init-val init-name))
860                             (pp indent+ ,(expr->string/MATLAB temp-val temp-name))
861                             (pp indent+ ,(expr->string/MATLAB beta-val beta-name))
862                             (pp indent+ ,(expr->string/MATLAB depth-val depth-name)))))
863                       epools))
864
865
866
867           (if (not (null? asgns))
868               (for-each (lambda (def)
869                           (let ((n (matlab-name (first def)) )
870                                 (b (second def)))
871                             (pp indent+ ,(expr->string/MATLAB b n)))) asgn-eq-defs))
872
873#|
874           (let* ((i-eqs (filter-map
875                          (lambda (ionch)
876
877                            (let* ((label     (first ionch))
878                                   (n         (second ionch))
879                                   (subcomps  ((dis 'component-subcomps) sys n))
880                                   (acc       (lookup-def 'accumulating-substance subcomps))
881                                   (perm      (lookup-def 'permeating-substance subcomps))
882                                   (permqs    (and perm ((dis 'component-exports) sys (cid perm))))
883                                   (pore      (lookup-def 'pore subcomps))
884                                   (gate      (lookup-def 'gate subcomps))
885                                   (sts       (and gate ((dis 'component-exports) sys (cid gate)))))
886
887                              (cond ((and perm pore gate)
888                                     (case (cn perm)
889                                       ((non-specific)
890                                        (let* ((i     (matlab-name 'i))
891                                               (e     (car permqs))
892                                               (gmax  (car ((dis 'component-exports) sys (cid pore))))
893                                               (pwrs  (map (lambda (n) (state-power sys n)) sts))
894                                               (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
895                                               (gion  `(* ,gmax ,@sptms)))
896                                          (list i e gion)))
897                                       (else
898                                        (let* ((i     (matlab-name (s+ 'i (cn perm))))
899                                               (e     (matlab-name (s+ 'e (cn perm))))
900                                               (gmax  (car ((dis 'component-exports) sys (cid pore))))
901                                               (pwrs  (map (lambda (n) (state-power sys n)) sts))
902                                               (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
903                                               (gion  `(* ,gmax ,@sptms)))
904                                          (list i e gion)))))
905
906                                    ((and perm pore)
907                                     (case (cn perm)
908                                       ((non-specific)
909                                        (let* ((i     (matlab-name 'i))
910                                               (e     (car permqs))
911                                               (gmax  (car ((dis 'component-exports) sys (cid pore)))))
912                                          (list i e gmax)))
913                                       (else
914                                        (nemo:error 'nemo:matlab-translator ": invalid ion channel definition " label))))
915
916                                      ((and acc pore gate)
917                                       (let* ((i     (matlab-name (s+ 'i (cn acc))))
918                                              (gmax  (car ((dis 'component-exports) sys (cid pore))))
919                                              (pwrs  (map (lambda (n) (state-power sys n)) sts))
920                                              (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
921                                              (gion  `(* ,gmax ,@sptms)))
922                                         (list i #f gion)))
923                                      (else (nemo:error 'nemo:matlab-translator ": invalid ion channel definition " label))
924                                      )))
925                          ionchs))
926                  (i-bkts (bucket-partition (lambda (x y) (eq? (car x) (car y))) i-eqs))
927                  (i-eqs  (fold (lambda (b ax)
928                                  (match b
929                                         ((and ps ((i e gion) . rst)) 
930                                          (let* ((sum   (if e (sum (map (lambda (b) `(* ,(third b) (- v ,(second b))))
931                                                                        ps))
932                                                            (sum (map third ps))))
933                                                 (sum0  (rhsexpr sum))
934                                                 (sum1  (canonicalize-expr/MATLAB sum0)))
935                                            (cons (list i sum1) ax)))
936                                           
937                                         ((i e gion)
938                                          (let* ((expr0  (rhsexpr (if e `(* ,gion (- v ,e)) gion)))
939                                                 (expr1  (canonicalize-expr/MATLAB expr0)))
940                                              (cons (list i expr1) ax)))
941                                         
942                                         (else ax)))
943                                (list) i-bkts))
944                  (locals (find-locals (map second i-eqs))))
945             (if (not (null? locals))    (pp indent+ (LOCAL ,(sl\ ", " locals))))
946             (if (not (null? asgns))     (pp indent+ (rates ())))
947             (if has-ode?
948                 (if (not method) (pp indent+ (SOLVE states))
949                     (pp indent+ (SOLVE states METHOD ,method))))
950             (if has-kinetic?   (pp indent+  (SOLVE kstates METHOD sparse)))
951             (if (not (null? stcomps))   (pp indent+ (stcomps ())))
952             (if (not (null? pool-ions)) (pp indent+ (pools ())))
953             (for-each (lambda (p) (pp indent+ ,(expr->string/MATLAB (second p) (first p)))) i-eqs)
954             (pp indent "}"))
955           
956           (if has-ode?
957               (begin
958                 (pp indent ,nl (DERIVATIVE states "{"))
959                 (let* ((state-eq-defs  (reverse (poset->state-eq-defs poset sys kinetic)))
960                        (pool-eq-defs
961                         (map (lambda (ep)
962                                (let* ((ep-name     (first ep))
963                                       (pool-ion    (assoc ep-name pool-ions))
964                                       (i-name      (second pool-ion))
965                                       (init-name   (matlab-name (s+ ep-name '-init)))
966                                       (temp-name   (matlab-name (s+ ep-name '-temp-adj)))
967                                       (beta-name   (matlab-name (s+ ep-name '-beta)))
968                                       (depth-name  (matlab-name (s+ ep-name '-depth)))
969                                       (rhs         `(let ((F 96485.0))
970                                                       (- (/ (neg ,i-name) (* 2 F ,init-name ,depth-name))
971                                                          (* ,beta-name ,ep-name .
972                                                             ,(if temp-name (list temp-name) (list)))))))
973                                  `(,(s+ ep-name "'") ,rhs)))
974                              epools))
975                        (eq-defs (append pool-eq-defs state-eq-defs))
976                        (locals  (find-locals (map second eq-defs))))
977                   (if (not (null? locals)) (pp indent+ (LOCAL ,(sl\ ", " locals))))
978                   (for-each (lambda (def)
979                               (let ((n (first def)) (b (second def)))
980                                 (pp indent+ ,(expr->string/MATLAB b n)))) eq-defs))
981                   
982                 (pp indent "}")))
983           
984|#
985
986)))))))
Note: See TracBrowser for help on using the repository browser.