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

Last change on this file since 13040 was 13040, checked in by Ivan Raikov, 12 years ago

Save.

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