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

Last change on this file since 12685 was 12685, checked in by Ivan Raikov, 13 years ago

Bug fixes.

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