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

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

Save.

File size: 26.1 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)  '^)  ((ln) 'log) (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->init-defs poset sys)
383  (fold-right
384   (lambda (lst ax)
385     (fold-right
386      (lambda (x ax) 
387              (match-let (((i . n)  x))
388                         (let ((en (environment-ref sys n)))
389                           (if (nemo:quantity? en)
390                               (cases nemo:quantity en
391                                      (REACTION  (name initial open transitions conserve power)
392                                                 (if (nemo:rhs? initial)
393                                                     (cons* (state-init name initial) 
394                                                            (state-init (matlab-state-name name open) name) ax) 
395                                                     ax))
396
397                                      (RATE  (name initial rhs)
398                                             (if (and initial (nemo:rhs? initial))
399                                                 (cons (state-init name initial) ax) 
400                                                 ax))
401
402                                      (else  ax))
403                               ax))))
404            ax lst))
405   (list) poset))
406
407
408(define (poset->state-conserve-eq-defs poset sys)
409  (fold-right
410   (lambda (lst ax)
411     (fold  (lambda (x ax) 
412              (match-let (((i . n)  x))
413                         (let ((en (environment-ref sys n)))
414                           (if (nemo:quantity? en)
415                               (cases nemo:quantity en
416                                      (REACTION (name initial open transitions conserve power)
417                                                (if (and (list? conserve) (every nemo:lineq? conserve))
418                                                    (cons (state-lineqs (matlab-name name) transitions conserve
419                                                                        matlab-state-name) ax) 
420                                                    ax))
421                                      (else  ax))
422                               ax))))
423            ax lst))
424   (list) poset))
425
426
427(define (reaction-power sys n)
428  (let ((en (environment-ref sys n)))
429    (if (nemo:quantity? en)
430        (cases nemo:quantity en
431               (REACTION  (name initial open transitions conserve power)  power)
432               (else  #f))  #f)))
433
434
435(define (bucket-partition p lst)
436  (let loop ((lst lst) (ax (list)))
437    (if (null? lst) ax
438        (let ((x (car lst)))
439          (let bkt-loop ((old-bkts ax) (new-bkts (list)))
440            (if (null? old-bkts) (loop (cdr lst) (cons (list x) new-bkts))
441                (if (p x (caar old-bkts ))
442                    (loop (cdr lst) (append (cdr old-bkts) (cons (cons x (car old-bkts)) new-bkts)))
443                    (bkt-loop (cdr old-bkts) (cons (car old-bkts) new-bkts)))))))))
444
445
446(define (nemo:matlab-translator sys . rest)
447  (define (cid x)  (second x))
448  (define (cn x)   (first x))
449  (let-optionals rest ((method 'lsode) (table? #f) (min-v -100) (max-v 100) (step 0.5) )
450  (match-let ((($ nemo:quantity 'DISPATCH  dis) (environment-ref sys (nemo-intern 'dispatch))))
451    (let ((imports  ((dis 'imports)  sys))
452          (exports  ((dis 'exports)  sys)))
453      (let* ((indent      0)
454             (indent+     (+ 2 indent ))
455             (eval-const  (dis 'eval-const))
456             (sysname     (matlab-name ((dis 'sysname) sys)))
457             (deps*       ((dis 'depgraph*) sys))
458             (consts      ((dis 'consts)  sys))
459             (asgns       ((dis 'asgns)   sys))
460             (states      ((dis 'states)  sys))
461             (reactions   ((dis 'reactions) sys))
462             (defuns      ((dis 'defuns)  sys))
463             (components  ((dis 'components) sys))
464             
465             (g             (match-let (((state-list asgn-list g) ((dis 'depgraph*) sys))) g))
466             (poset         (vector->list ((dis 'depgraph->bfs-dist-poset) g)))
467
468             (const-defs       (filter-map
469                                (lambda (nv)
470                                  (and (not (member (first nv) matlab-builtin-consts))
471                                       (let ((v1 (canonicalize-expr/MATLAB (second nv))))
472                                         (list (matlab-name (first nv)) v1))))
473                                consts))
474             
475             (ionch-info    (nemo:ionch-query sys))
476             (ionchs        (lookup-def 'ion-channels ionch-info))
477             (perm-ions     (map (match-lambda ((comp i e erev) `(,comp ,(matlab-name i) ,(matlab-name e) ,erev)))
478                                 (lookup-def 'perm-ions ionch-info)))
479             (acc-ions      (map (match-lambda ((comp i in out) `(,comp ,@(map matlab-name (list i in out)))))
480                                 (lookup-def 'acc-ions ionch-info)))
481             (epools        (lookup-def 'pool-ions ionch-info))
482             (pool-ions     (map (lambda (lst) (map matlab-name lst)) epools))
483
484             (i-gates       (lookup-def 'i-gates ionch-info))
485
486             (capcomp       (any (match-lambda ((name 'membrane-capacitance id) (list name id)) (else #f)) components))
487             (mcap          (and capcomp (car ((dis 'component-exports) sys (cid capcomp)))))
488               
489
490             (i-eqs0 (filter-map
491                     (lambda (ionch) 
492                       
493                       (let* ((label     (first ionch))
494                              (n         (second ionch))
495                              (subcomps  ((dis 'component-subcomps) sys n))
496                              (acc       (lookup-def 'accumulating-substance subcomps))
497                              (perm      (lookup-def 'permeating-substance subcomps))
498                              (permqs    (and perm ((dis 'component-exports) sys (cid perm))))
499                              (pore      (lookup-def 'pore subcomps))
500                              (gate      (lookup-def 'gate subcomps))
501                              (sts       (and gate ((dis 'component-exports) sys (cid gate)))))
502                         
503                         (cond ((and perm pore gate)
504                                (case (cn perm)
505                                  ((non-specific)
506                                   (let* ((i     (matlab-name 'i))
507                                          (e     (car permqs))
508                                          (gmax  (car ((dis 'component-exports) sys (cid pore))))
509                                          (pwrs  (map (lambda (n) (reaction-power sys n)) sts))
510                                          (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
511                                          (gion  `(* ,gmax ,@sptms)))
512                                     (list i e gion)))
513                                  (else
514                                   (let* ((i     (matlab-name (s+ 'i (cn perm))))
515                                          (e     (car permqs))
516                                          (gmax  (car ((dis 'component-exports) sys (cid pore))))
517                                          (pwrs  (map (lambda (n) (reaction-power sys n)) sts))
518                                          (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
519                                          (gion  `(* ,gmax ,@sptms)))
520                                     (list i e gion)))))
521                               
522                               ((and perm pore)
523                                (case (cn perm)
524                                  ((non-specific)
525                                   (let* ((i     (matlab-name 'i))
526                                          (e     (car permqs))
527                                          (gmax  (car ((dis 'component-exports) sys (cid pore)))))
528                                     (list i e gmax)))
529                                  (else
530                                   (nemo:error 'nemo:matlab-translator ": invalid ion channel definition " label))))
531                               
532                               ((and acc pore gate)
533                                (let* ((i     (matlab-name (s+ 'i (cn acc))))
534                                       (gmax  (car ((dis 'component-exports) sys (cid pore))))
535                                       (pwrs  (map (lambda (n) (reaction-power sys n)) sts))
536                                       (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
537                                       (gion  `(* ,gmax ,@sptms)))
538                                  (list i #f gion)))
539                               (else (nemo:error 'nemo:matlab-translator ": invalid ion channel definition " label))
540                               )))
541                     ionchs))
542               
543             (i-bkts (bucket-partition (lambda (x y) (eq? (car x) (car y))) i-eqs0))
544             
545             (i-eqs  (fold (lambda (b ax) 
546                             (match b 
547                                    ((and ps ((i e gion) . rst)) 
548                                     (let* ((sum   (if e (sum (map (lambda (b) `(* ,(third b) (- v ,(second b)))) 
549                                                                   ps))
550                                                       (sum (map third ps))))
551                                            (sum0  (rhsexpr/MATLAB sum))
552                                            (sum1  (canonicalize-expr/MATLAB sum0)))
553                                       (cons (list i sum1) ax)))
554                                   
555                                    ((i e gion)
556                                     (let* ((expr0  (rhsexpr/MATLAB (if e `(* ,gion (- v ,e)) gion)))
557                                            (expr1  (canonicalize-expr/MATLAB expr0)))
558                                       (cons (list i expr1) ax)))
559                                   
560                                    (else ax)))
561                           (list) i-bkts))
562
563             (asgn-eq-defs     (poset->asgn-eq-defs poset sys))
564             
565             
566             (rate-eq-defs       (reverse (poset->rate-eq-defs poset sys method)))
567             
568             (reaction-eq-defs   (poset->reaction-eq-defs poset sys))
569             
570             (init-eq-defs       (poset->init-defs poset sys))
571             
572             (conserve-eq-defs   (map (lambda (eq) (list 0 `(- ,(second eq) ,(first eq)))) 
573                                      (poset->state-conserve-eq-defs poset sys)))
574             
575             (pool-ion-i          (map second pool-ions))
576             
577             
578             (state-index-map  (let ((acc (fold (lambda (def ax)
579                                                  (let ((st-name (first def)))
580                                                    (list (+ 1 (first ax)) 
581                                                          (cons `(,st-name ,(first ax)) (second ax)))))
582                                                (list 1 (list)) 
583                                                (cons (list 'v) rate-eq-defs))))
584                                 
585                                 (second acc)))
586             
587             (steady-state-index-map  (let ((acc (fold (lambda (def ax)
588                                                         (let ((st-name (first def)))
589                                                           (if (not (alist-ref st-name init-eq-defs))
590                                                               (list (+ 1 (first ax)) 
591                                                                     (cons `(,st-name ,(first ax)) (second ax)))
592                                                               ax)))
593                                                       (list 1 (list)) 
594                                                       rate-eq-defs)))
595                                        (second acc)))
596             
597             (globals          (map matlab-name 
598                                    (delete-duplicates (append
599                                                        exports
600                                                        (map second  perm-ions)
601                                                        (map third   perm-ions)
602                                                        (map second  acc-ions)
603                                                        (map third   acc-ions)
604                                                        (map fourth  acc-ions)
605                                                        (map second  pool-ions)
606                                                        (map third   pool-ions)
607                                                        (map first   imports) 
608                                                        (map first   const-defs)))))
609
610               )
611
612           (for-each
613            (lambda (a)
614              (let ((acc-ion   (car a)))
615                (if (assoc acc-ion perm-ions)
616                    (nemo:error 'nemo:matlab-translator 
617                                ": ion species " acc-ion " cannot be declared as both accumulating and permeating"))))
618            acc-ions)
619
620           (for-each
621            (lambda (p)
622              (let ((pool-ion  (car p)))
623                (if (assoc pool-ion perm-ions)
624                    (nemo:error 'nemo:matlab-translator 
625                                ": ion species " pool-ion " cannot be declared as both pool and permeating"))))
626            pool-ions)
627
628           (if (not (null? globals)) (pp indent (global ,(sl\ " " globals))))
629
630           (pp indent ,nl (function dy = ,sysname (,(sl\ ", " (case method
631                                                                ((lsode)  `(y t))
632                                                                ((odepkg) `(t y))
633                                                                (else     `(y t)))))))
634
635           (if (not (null? globals)) (pp indent+ (global ,(sl\ " " globals))))
636           
637           (if (member 'v globals)
638               (let ((vi (lookup-def 'v state-index-map)))
639                 (if vi (pp indent+ ,(expr->string/MATLAB (s+ "y(" vi ")") 'v)))))
640
641           
642           (for-each (lambda (def)
643                       (let ((n (first def)) )
644                         (pp indent+ ,(expr->string/MATLAB 
645                                       (s+ "y(" (lookup-def n state-index-map) ")") n))))
646                     rate-eq-defs)
647
648           (for-each (lambda (def)
649                       (let ((n (matlab-name (first def)) )
650                             (b (second def)))
651                         (pp indent+ ,(expr->string/MATLAB b n)))) 
652                     asgn-eq-defs)
653           
654           (for-each (lambda (def)
655                       (let ((n (first def)) (b (second def)))
656                         (pp indent+ ,(expr->string/MATLAB b n)))) 
657                     reaction-eq-defs)
658
659           (for-each (lambda (pool-ion)
660                       (pp indent+ ,(expr->string/MATLAB (first pool-ion) (third pool-ion) )))
661                     pool-ions)
662
663         (pp indent+ ,(expr->string/MATLAB `(zeros (length y) 1) 'dy))
664         
665
666         (for-each (lambda (def)
667                     (let ((n (s+ "dy(" (lookup-def (first def) state-index-map) ")")) 
668                           (b (second def)))
669                       (pp indent+ ,(s+ "# state " (first def))
670                           ,(expr->string/MATLAB b n))))
671                   rate-eq-defs)
672         
673         (for-each (lambda (def) 
674                     (pp indent+ ,(expr->string/MATLAB (second def) (first def)))) 
675                   i-eqs)
676         
677         (if (and mcap (member 'v globals))
678             (let ((vi (lookup-def 'v state-index-map)))
679               (if vi (pp indent+  ,(expr->string/MATLAB `(/ (neg ,(sum (map first i-eqs))) ,mcap) 
680                                                         (s+ "dy(" vi ")"))))))
681         
682         (pp indent ,nl (endfunction))
683
684         (if (not (null? steady-state-index-map))
685             (begin
686               (pp indent ,nl (function y = ,(s+ sysname "_steadystate") (x)))
687               
688               (if (not (null? globals)) (pp indent+ (global ,(sl\ " " globals))))
689               
690               (for-each
691                (lambda (def)
692                  (let* ((n   (first def)) 
693                         (ni  (lookup-def n steady-state-index-map)))
694                    (if ni (pp indent+ ,(expr->string/MATLAB (s+ "x(" ni ")") n)))
695                    ))
696                rate-eq-defs)
697               
698               (pp indent+ ,(expr->string/MATLAB `(zeros ,(length steady-state-index-map) 1) 'y))
699               
700               (for-each (lambda (def)
701                           (let ((n (matlab-name (first def)) )
702                                 (b (second def)))
703                             (pp indent+ ,(expr->string/MATLAB b n)))) asgn-eq-defs)
704               
705               (for-each
706                (lambda (def)
707                  (let* ((n   (first def)) 
708                         (ni  (lookup-def n steady-state-index-map))
709                         (b   (second def)))
710                    (if ni (pp indent+ ,(s+ "# state " n)
711                               ,(expr->string/MATLAB b (s+ "y(" ni ")"))))
712                    ))
713                rate-eq-defs)
714               
715               (pp indent ,nl (endfunction))))
716             
717
718         (pp indent ,nl (function  ,(s+ sysname "_print_state") (y)))
719
720         (let ((lst (sort (map (lambda (x) (cons (->string (car x)) (cdr x))) state-index-map)
721                          (lambda (x y) (string<? (car x) (car y))))))
722           (for-each (lambda (p)
723                       (let ((n (first p)) (i (second p)))
724                         (pp indent+ (,n " = " "y(" ,i ")"))))
725                     lst))
726
727         (pp indent ,nl (endfunction))
728
729         
730         (pp indent ,nl (function y0 = ,(s+ sysname "_init") (Vinit)))
731         
732         (if (not (null? globals)) (pp indent+ (global ,(sl\ " " globals))))
733         
734         (pp indent+ ,(expr->string/MATLAB `(zeros ,(length state-index-map) 1) 'y0))
735         
736         (if (member 'v globals)
737             (let ((vi (lookup-def 'v state-index-map)))
738               (pp indent+ ,(expr->string/MATLAB `Vinit 'v))
739               (if vi (pp indent+ ,(expr->string/MATLAB 'v (s+ "y0(" vi ")") )))))
740         
741         (for-each (lambda (def)
742                     (let ((n (first def)) (b (second def)))
743                       (pp indent+ ,(expr->string/MATLAB b n)))) 
744                   const-defs)
745         
746         (for-each (lambda (def)
747                     (let* ((n  (first def)) 
748                            (ni (lookup-def n state-index-map))
749                            (b (second def)))
750                       (pp indent+ ,(expr->string/MATLAB b n))
751                       (if ni (pp indent+ ,(expr->string/MATLAB n  (s+ "y0(" ni ")"))))))
752                   init-eq-defs)
753         
754         (for-each (lambda (pool-ion)
755                     (pp indent+ ,(expr->string/MATLAB (first pool-ion) (third pool-ion) )))
756                   pool-ions)
757
758         (if (not (null? steady-state-index-map))
759             (begin
760               (for-each (lambda (def) 
761                                 (if (member (first def) pool-ion-i)
762                                     (pp indent+ ,(expr->string/MATLAB (second def) (first def)))))
763                         i-eqs)
764
765               (pp indent+ ,(expr->string/MATLAB `(zeros ,(length steady-state-index-map) 1) 'xs)
766                   ,(expr->string/MATLAB 
767                     `(fsolve ,(s+ #\" sysname "_steadystate" #\" ) xs) 
768                     "[ys,fval,info]"))
769               
770
771               (for-each
772                (lambda (def)
773                  (let* ((n (first def)) 
774                         (si (lookup-def n steady-state-index-map))
775                         (ni (lookup-def n state-index-map)))
776                    (if si (pp indent+ ,(expr->string/MATLAB (s+ "ys(" si ")") (s+ "y0(" ni ")"))))))
777                rate-eq-defs)
778               ))
779
780
781         (for-each (lambda (def) 
782                     (if (not (member (first def) pool-ion-i))
783                         (pp indent+ ,(expr->string/MATLAB (second def) (first def)))))
784                   i-eqs)
785
786         (pp indent ,nl (endfunction))
787         
788         (pp indent ,nl)
789         
790         (let* ((with       (inexact->exact (round (/ (abs (- max-v min-v)) step))))
791                (define-fn  (make-define-fn table? min-v max-v with)))
792           (for-each (lambda (fndef) 
793                       (if (not (member (car fndef) builtin-fns))
794                           (apply define-fn (cons indent fndef)))) 
795                     defuns))
796         )))))
Note: See TracBrowser for help on using the repository browser.