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

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

Place the main function at the beginning of the file.

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