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

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

Added some more detailed error messages; restructured hh example.

File size: 25.2 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* ((poset            (vector->list ((dis 'depgraph->bfs-dist-poset) g)))
461                (asgn-eq-defs     (poset->asgn-eq-defs poset sys))
462                (perm-ions (fold (lambda (ionch ax) 
463                                    (let* ((subcomps ((dis 'component-subcomps) sys (cid ionch)))
464                                           (perm      (lookup-def 'permeating-substance subcomps)))
465                                      (if perm 
466                                          (case (cn perm)
467                                            ((non-specific)   
468                                             (let* ((erev (car ((dis 'component-exports) sys (cid perm))))
469                                                    (i    (matlab-name 'i))
470                                                    (e    (matlab-name 'e)))
471                                               (cons `(,(cn perm) ,i ,e ,erev) ax)))
472                                            (else (let* ((erev (car ((dis 'component-exports) sys (cid perm))))
473                                                         (i    (matlab-name (s+ 'i (cn perm))))
474                                                         (e    (matlab-name (s+ 'e (cn perm)))))
475                                                    (cons `(,(cn perm) ,i ,e ,erev) ax))))
476                                          ax)))
477                                  (list) ionchs))
478               (acc-ions (fold (lambda (ionch ax) 
479                                  (let* ((subcomps ((dis 'component-subcomps) sys (cid ionch)))
480                                         (acc   (lookup-def 'accumulating-substance subcomps))
481                                         (i     (and acc (matlab-name (s+ 'i (cn acc)))))
482                                         (in    (and acc (matlab-name (s+ (cn acc) 'i))))
483                                         (out   (and acc (matlab-name (s+ (cn acc) 'o)))))
484                                    (if acc  (cons `(,(cn acc) ,i ,in ,out) ax) ax)))
485                                (list) ionchs))
486               (pool-ions (map (lambda (ep)
487                                  (let ((ion (car ep)))
488                                    `(,(matlab-name ion) ,(matlab-name (s+ 'i ion)) ,(matlab-name (s+ ion 'i)))))
489                               epools)))
490
491           (for-each
492            (lambda (a)
493              (let ((acc-ion   (car a)))
494                (if (assoc acc-ion perm-ions)
495                    (nemo:error 'nemo:matlab-translator 
496                                ": ion species " acc-ion " cannot be declared as both accumulating and permeating"))))
497            acc-ions)
498
499           (for-each
500            (lambda (p)
501              (let ((pool-ion  (car p)))
502                (if (assoc pool-ion perm-ions)
503                    (nemo:error 'nemo:matlab-translator 
504                                ": ion species " pool-ion " cannot be declared as both pool and permeating"))))
505            pool-ions)
506
507           (let* ((with      (inexact->exact (round (/ (abs (- max-v min-v)) step))))
508                  (define-fn (make-define-fn table? min-v max-v with)))
509             (for-each (lambda (fndef) 
510                         (if (not (member (car fndef) builtin-fns))
511                             (apply define-fn (cons indent fndef)))) 
512                       defuns))
513
514           (pp indent ,nl)
515
516           (let* ((const-defs (filter-map
517                               (lambda (nv)
518                                 (and (not (member (first nv) matlab-builtin-consts))
519                                      (let ((v1 (canonicalize-expr/MATLAB (second nv))))
520                                        (list (first nv) v1))))
521                                   consts)))
522
523             (for-each (lambda (def)
524                         (let ((n (first def)) (b (second def)))
525                           (pp indent ,(expr->string/MATLAB b n)))) const-defs)
526
527             (for-each (lambda (ep)
528                         (let* ((ep-name     (first ep))
529                                (ep-props    (second ep))
530                                (init-expr   (lookup-def 'initial ep-props))
531                                (temp-expr   (lookup-def 'temp-adj ep-props))
532                                (beta-expr   (lookup-def 'beta ep-props))
533                                (depth-expr  (lookup-def 'depth ep-props))
534                                (init-name   (matlab-name (s+ ep-name '-init)))
535                                (temp-name   (matlab-name (s+ ep-name '-temp-adj)))
536                                (beta-name   (matlab-name (s+ ep-name '-beta)))
537                                (depth-name  (matlab-name (s+ ep-name '-depth))))
538                           (if (or (not beta-expr) (not depth-expr) (not init-expr))
539                               (nemo:error 'nemo:matlab-translator 
540                                           ": ion pool " ep-name " requires initial value, depth and beta parameters"))
541                           (let ((temp-val  (and temp-expr (eval-const sys temp-expr)))
542                                 (init-val  (eval-const sys init-expr))
543                                 (beta-val  (eval-const sys beta-expr))
544                                 (depth-val (eval-const sys depth-expr)))
545                             (pp indent 
546                                 ,(expr->string/MATLAB init-val init-name)
547                                 ,(expr->string/MATLAB temp-val temp-name)
548                                 ,(expr->string/MATLAB beta-val beta-name)
549                                 ,(expr->string/MATLAB depth-val depth-name)))))
550                       epools))
551
552             (pp indent ,nl (function dy = ,sysname (,(sl\ ", " `(dy t)))))
553
554             (if (not (null? exports)) (pp indent+ (global ,(sl\ ", " exports))))
555
556             (if (not (null? asgns))
557                 (for-each (lambda (def)
558                             (let ((n (matlab-name (first def)) )
559                                   (b (second def)))
560                               (pp indent+ ,(expr->string/MATLAB b n)))) asgn-eq-defs))
561
562           (let* ((i-eqs (filter-map
563                          (lambda (ionch) 
564
565                            (let* ((label     (first ionch))
566                                   (n         (second ionch))
567                                   (subcomps  ((dis 'component-subcomps) sys n))
568                                   (acc       (lookup-def 'accumulating-substance subcomps))
569                                   (perm      (lookup-def 'permeating-substance subcomps))
570                                   (permqs    (and perm ((dis 'component-exports) sys (cid perm))))
571                                   (pore      (lookup-def 'pore subcomps))
572                                   (gate      (lookup-def 'gate subcomps))
573                                   (sts       (and gate ((dis 'component-exports) sys (cid gate)))))
574
575                              (cond ((and perm pore gate)
576                                     (case (cn perm)
577                                       ((non-specific)
578                                        (let* ((i     (matlab-name 'i))
579                                               (e     (car permqs))
580                                               (gmax  (car ((dis 'component-exports) sys (cid pore))))
581                                               (pwrs  (map (lambda (n) (state-power sys n)) sts))
582                                               (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
583                                               (gion  `(* ,gmax ,@sptms)))
584                                          (list i e gion)))
585                                       (else
586                                        (let* ((i     (matlab-name (s+ 'i (cn perm))))
587                                               (e     (matlab-name (s+ 'e (cn perm))))
588                                               (gmax  (car ((dis 'component-exports) sys (cid pore))))
589                                               (pwrs  (map (lambda (n) (state-power sys n)) sts))
590                                               (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
591                                               (gion  `(* ,gmax ,@sptms)))
592                                          (list i e gion)))))
593
594                                    ((and perm pore)
595                                     (case (cn perm)
596                                       ((non-specific)
597                                        (let* ((i     (matlab-name 'i))
598                                               (e     (car permqs))
599                                               (gmax  (car ((dis 'component-exports) sys (cid pore)))))
600                                          (list i e gmax)))
601                                       (else
602                                        (nemo:error 'nemo:matlab-translator ": invalid ion channel definition " label))))
603
604                                      ((and acc pore gate)
605                                       (let* ((i     (matlab-name (s+ 'i (cn acc))))
606                                              (gmax  (car ((dis 'component-exports) sys (cid pore))))
607                                              (pwrs  (map (lambda (n) (state-power sys n)) sts))
608                                              (sptms (map (lambda (st pwr) `(pow ,st ,pwr)) sts pwrs))
609                                              (gion  `(* ,gmax ,@sptms)))
610                                         (list i #f gion)))
611                                      (else (nemo:error 'nemo:matlab-translator ": invalid ion channel definition " label))
612                                      )))
613                          ionchs))
614
615                  (i-bkts (bucket-partition (lambda (x y) (eq? (car x) (car y))) i-eqs))
616
617                  (i-eqs  (fold (lambda (b ax) 
618                                  (match b 
619                                         ((and ps ((i e gion) . rst)) 
620                                          (let* ((sum   (if e (sum (map (lambda (b) `(* ,(third b) (- v ,(second b)))) 
621                                                                        ps))
622                                                            (sum (map third ps))))
623                                                 (sum0  (rhsexpr/MATLAB sum))
624                                                 (sum1  (canonicalize-expr/MATLAB sum0)))
625                                            (cons (list i sum1) ax)))
626                                           
627                                         ((i e gion)
628                                          (let* ((expr0  (rhsexpr/MATLAB (if e `(* ,gion (- v ,e)) gion)))
629                                                 (expr1  (canonicalize-expr/MATLAB expr0)))
630                                              (cons (list i expr1) ax)))
631                                         
632                                         (else ax)))
633                                (list) i-bkts))
634
635                  (state-eq-defs    (reverse (poset->state-eq-defs poset sys)))
636
637                  (stcomp-eq-defs   (poset->stcomp-eq-defs poset sys))
638
639                  (pool-eq-defs
640                   (map (lambda (ep)
641                          (let* ((ep-name     (first ep))
642                                 (pool-ion    (assoc ep-name pool-ions))
643                                 (i-name      (second pool-ion))
644                                 (init-name   (matlab-name (s+ ep-name '-init)))
645                                 (temp-name   (matlab-name (s+ ep-name '-temp-adj)))
646                                 (beta-name   (matlab-name (s+ ep-name '-beta)))
647                                 (depth-name  (matlab-name (s+ ep-name '-depth)))
648                                 (rhs         `(let ((F 96485.0))
649                                                 (- (/ (neg ,i-name) (* 2 F ,init-name ,depth-name)) 
650                                                    (* ,beta-name ,ep-name . 
651                                                       ,(if temp-name (list temp-name) (list)))))))
652                            `(,(s+ ep-name)  ,rhs)))
653                        epools))
654
655                  (rate-eq-defs (append pool-eq-defs state-eq-defs))
656                 
657                  (state-index-map  (let ((acc (fold (lambda (def ax)
658                                                       (let ((st-name (first def)))
659                                                         (list (+ 1 (first ax)) 
660                                                               (cons `(,st-name ,(first ax)) (second ax)))))
661                                                     (list 0 (list)) rate-eq-defs)))
662                                    (second acc)))
663                  )
664
665
666
667             (for-each (lambda (def)
668                         (let ((n (first def)) )
669                           (pp indent+ (,n = ,(s+ "y(" (lookup-def n state-index-map) ")"))))) rate-eq-defs)
670
671             (for-each (lambda (def)
672                         (let ((n (first def)) (b (second def)))
673                           (pp indent+ ,(expr->string/MATLAB b n)))) stcomp-eq-defs)
674
675             (pp indent+ (dy = zeros(length(y) #\, 1)))
676
677             (for-each (lambda (def)
678                         (let ((n (s+ "dy(" (lookup-def (first def) state-index-map) ")")) (b (second def)))
679                           (pp indent+ ,(expr->string/MATLAB b n)))) rate-eq-defs)
680
681             (for-each (lambda (pool-ion)
682                         (pp indent+ (,(third pool-ion) = ,(first pool-ion)))) pool-ions)
683
684             (for-each (lambda (def) 
685                         (pp indent+ ,(expr->string/MATLAB (second def) (first def)))) i-eqs)
686             
687             (pp indent ,nl (endfunction))
688
689             (pp indent ,nl (function y0 = ,(s+ sysname "_init") (v)))
690
691             (if (not (null? exports)) (pp indent+ (global ,(sl\ ", " exports))))
692
693             (let* ((init-defs         (poset->state-init-defs poset sys)))
694;;                  (init-eq-defs      (poset->state-init-eq-defs poset sys)))
695
696               (for-each (lambda (def)
697                           (let ((n (matlab-name (first def)) )
698                                 (b (second def)))
699                             (pp indent+ ,(expr->string/MATLAB b n)))) asgn-eq-defs)
700
701               (for-each (lambda (def)
702                           (let* ((n  (first def)) 
703                                  (ni (lookup-def n state-index-map))
704                                  (b (second def)))
705                             (pp indent+ ,(expr->string/MATLAB b (if ni (s+ "y0(" ni ")") n))))) 
706                         init-defs)
707
708#|
709               (for-each
710                (lambda (x)
711                  (let ((lineqs  (second x)))
712                    (for-each (lambda (eq)
713                                (let ((val  (first eq))
714                                      (expr (third eq)))
715                                  (pp indent+ ,(lineq->string/NMODL expr val))))
716                              lineqs)))
717                init-eq-defs)
718|#
719
720               )
721             (pp indent ,nl (endfunction))
722             
723             
724             ))))))))
Note: See TracBrowser for help on using the repository browser.