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

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

Factoring out the code for creating graph representation of the kinetic equations.

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