Context Navigation

source: project/release/4/nemo/trunk/nemo-matlab.scm @ 25870

Last change on this file since 25870 was 25870, checked in by Ivan Raikov, 9 years ago
nemo: updated copyright year
File size: 34.8 KB

Line
1	;;
2	;;
3	;; An extension for generating Matlab/Octave code from NEMO models.
4	;;
5	;; Copyright 2008-2012 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	(module nemo-matlab
22
23	(nemo:matlab-translator
24	nemo:octave-translator)
25
26	(import scheme chicken utils data-structures lolevel ports srfi-1 srfi-13)
27
28	(require-extension lolevel matchable strictly-pretty environments
29	varsubst datatype nemo-core nemo-utils
30	nemo-gate-complex)
31
32	(define matlab-builtin-consts
33	`())
34
35	(define matlab-ops
36	`(+ - * / > < <= >= = ^))
37
38	(define builtin-fns
39	`(+ - * / pow neg abs atan asin acos sin cos exp ln
40	sqrt tan cosh sinh tanh hypot gamma lgamma log10 log2 log1p ldexp cube
41	> < <= >= = and or round ceiling floor max min
42	fpvector-ref))
43
44	(define (matlab-name s)
45	(let ((cs (string->list (->string s))))
46	(let loop ((lst (list)) (cs cs))
47	(if (null? cs) (string->symbol (list->string (reverse lst)))
48	(let* ((c (car cs))
49	(c1 (cond ((or (char-alphabetic? c) (char-numeric? c) (char=? c #\_)) c)
50	(else #\_))))
51	(loop (cons c1 lst) (cdr cs)))))))
52
53	(define (rhsexpr/MATLAB expr)
54	(match expr
55	(('if . es) `(if . ,(map (lambda (x) (rhsexpr/MATLAB x)) es)))
56	(('pow x y) (if (and (integer? y) (positive? y))
57	(if (> y 1) (let ((tmp (gensym "x")))
58	`(let ((,tmp ,x)) (* . ,(list-tabulate (inexact->exact y) (lambda (i) tmp)))))
59	x)
60	expr))
61	((s . es) (if (symbol? s) (cons (if (member s builtin-fns) s (matlab-name s)) (map (lambda (x) (rhsexpr/MATLAB x)) es)) expr))
62	(id (if (symbol? id) (matlab-name id) id))))
63
64	(define (matlab-state-name n s)
65	(matlab-name (s+ n s)))
66
67	(define-syntax pp
68	(syntax-rules ()
69	((pp indent val ...) (ppf indent (quasiquote val) ...))))
70
71
72	(define group/MATLAB (doc:block 2 (doc:text "(") (doc:text ")")))
73	(define block/MATLAB (doc:block 2 (doc:empty) (doc:empty)))
74	(define (stmt/MATLAB x)
75	(match x
76	(($ doc 'DocCons _ ($ doc 'DocText ";")) x)
77	(else (doc:cons x (doc:text ";")))))
78
79
80	(define (ifthen/MATLAB c e1 e2)
81	(doc:nest
82	2 (doc:connect
83	(doc:connect (doc:group (doc:connect (doc:text "if") c))
84	(doc:connect (doc:nest 2 e1)
85	(doc:nest 2 (doc:connect (doc:text "else") e2))))
86	(doc:text "end"))))
87
88	(define (letblk/MATLAB e1 e2)
89	(cond ((equal? e1 (doc:empty)) (doc:group (doc:nest 2 e2)))
90	((equal? e2 (doc:empty)) (doc:group (doc:nest 2 e1)))
91	(else (doc:connect (doc:group (doc:nest 2 (stmt/MATLAB e1)))
92	(doc:group (doc:nest 2 e2))))))
93
94
95	(define (format-op/MATLAB indent op args)
96	(let ((op1 (doc:text (->string op))))
97	(if (null? args) op1
98	(match args
99	((x) (doc:concat (list op1 x)))
100	((x y) (doc:concat (intersperse (list x op1 y) (doc:space))))
101	((x y z) (doc:concat (intersperse (list x op1 y op1 z) (doc:space))))
102	(lst (let* ((n (length lst))
103	(n/2 (inexact->exact (round (/ n 2)))))
104	(doc:concat
105	(intersperse
106	(list (format-op/MATLAB indent op (take lst n/2 )) op1
107	(format-op/MATLAB indent op (drop lst n/2 )))
108	(doc:space)))))))))
109
110	(define (format-fncall/MATLAB indent op args)
111	(let ((op1 (doc:text (->string op))))
112	(doc:cons op1 (group/MATLAB ((doc:list indent identity (lambda () (doc:text ", "))) args)))))
113
114	(define (name-normalize expr)
115	(match expr
116	(('if c t e) `(if ,(name-normalize c) ,(name-normalize t) ,(name-normalize e)))
117	(('let bs e)
118	`(let ,(map (lambda (b) `(,(car b) ,(name-normalize (cadr b)))) bs) ,(name-normalize e)))
119	((f . es)
120	(cons f (map name-normalize es)))
121	((? symbol? ) (matlab-name expr))
122	((? atom? ) expr)))
123
124
125	(define (canonicalize-expr/MATLAB expr)
126	(let ((subst-convert (subst-driver (lambda (x) (and (symbol? x) x)) nemo:binding? identity nemo:bind nemo:subst-term)))
127	(let* ((expr1 (if-convert expr))
128	(expr2 (subst-convert expr1 subst-empty))
129	(expr3 (let-lift expr2))
130	(expr4 (name-normalize expr3)))
131	expr4)))
132
133
134	(define (format-expr/MATLAB indent expr . rest)
135	(let-optionals rest ((rv #f))
136	(let ((indent+ (+ 2 indent)))
137	(match expr
138	(('let bindings body)
139	(letblk/MATLAB
140	(fold-right
141	(lambda (x ax)
142	(letblk/MATLAB
143	(match (second x)
144	(('if c t e)
145	(ifthen/MATLAB
146	(group/MATLAB (format-expr/MATLAB indent c))
147	(block/MATLAB (format-expr/MATLAB indent t (first x)))
148	(block/MATLAB (format-expr/MATLAB indent e (first x)))))
149	(else
150	(stmt/MATLAB
151	(format-op/MATLAB indent+ " = "
152	(list (format-expr/MATLAB indent (first x) )
153	(format-expr/MATLAB indent (second x)))))))
154	ax))
155	(doc:empty) bindings)
156	(match body
157	(('let _ _) (format-expr/MATLAB indent body rv))
158	(else
159	(let ((body1 (doc:nest indent (format-expr/MATLAB indent body))))
160	(if rv (stmt/MATLAB (format-op/MATLAB indent " = " (list (format-expr/MATLAB indent+ rv ) body1)))
161	body1))))))
162
163	(('if . rest) (error 'format-expr/MATLAB "invalid if statement " expr))
164
165	((op . rest)
166	(let ((op (case op ((pow) '^) ((ln) 'log) (else op))))
167	(let ((fe
168	(if (member op matlab-ops)
169	(let ((mdiv? (any (lambda (x) (match x (('* . _) #t) (('/ . _) #t) (else #f))) rest))
170	(mul? (any (lambda (x) (match x (('* . _) #t) (else #f))) rest))
171	(plmin? (any (lambda (x) (match x (('+ . _) #t) (('- . _) #t) (else #f))) rest)))
172	(case op
173	((/)
174	(format-op/MATLAB indent op
175	(map (lambda (x)
176	(let ((fx (format-expr/MATLAB indent+ x)))
177	(if (or (symbol? x) (number? x)) fx
178	(if (or mul? plmin?) (group/MATLAB fx) fx)))) rest)))
179	((*)
180	(format-op/MATLAB indent op
181	(map (lambda (x)
182	(let ((fx (format-expr/MATLAB indent+ x)))
183	(if (or (symbol? x) (number? x)) fx
184	(if plmin? (group/MATLAB fx) fx)))) rest)))
185
186	((^)
187	(format-op/MATLAB indent op
188	(map (lambda (x)
189	(let ((fx (format-expr/MATLAB indent+ x)))
190	(if (or (symbol? x) (number? x)) fx
191	(if (or mdiv? plmin?) (group/MATLAB fx) fx)))) rest)))
192
193	(else
194	(format-op/MATLAB indent op
195	(map (lambda (x)
196	(let ((fx (format-expr/MATLAB indent+ x))) fx)) rest)))))
197
198	(let ((op (case op ((neg) '-) (else op))))
199	(format-fncall/MATLAB indent op (map (lambda (x) (format-expr/MATLAB indent+ x)) rest))))))
200	(if rv
201	(stmt/MATLAB (format-op/MATLAB indent " = " (list (format-expr/MATLAB indent+ rv ) fe)))
202	fe))))
203
204	(else (let ((fe (doc:text (->string expr))))
205	(if rv
206	(stmt/MATLAB (format-op/MATLAB indent " = " (list (format-expr/MATLAB indent+ rv ) fe)))
207	fe)))))))
208
209
210	(define (expr->string/MATLAB x . rest)
211	(let-optionals rest ((rv #f) (width 72))
212	(sdoc->string (doc:format width (format-expr/MATLAB 2 x rv)))))
213
214
215	(define (expeuler dt name rhs)
216	(define (isname? x) (equal? x name))
217	(let ((res
218	(match rhs
219	((or ('- A ('* B (and x (? isname?))))
220	('+ ('neg ('* B (and x (? isname?)))) A))
221	(let ((xexp (string->symbol (s+ x 'exp))))
222	`(let ((,xexp (exp (* (neg ,B) ,dt))))
223	(+ (* ,x ,xexp) (* (- 1 ,xexp) (/ ,A ,B))))))
224
225	((or ('- A ('* (and x (? isname?)) . B))
226	('+ ('neg ('* (and x (? isname?)) . B)) A))
227	(let ((xexp (string->symbol (s+ x 'exp)))
228	(B1 (if (null? (cdr B)) (car B) `(* ,@B))))
229	`(let ((,xexp (exp (* (neg ,B1) ,dt))))
230	(+ (* ,x ,xexp) (* (- 1 ,xexp) (/ ,A ,B1))))))
231
232	(('+ ('neg ('* (and x1 (? isname?)) Alpha))
233	('* ('- 1 (and x2 (? isname?))) Beta))
234	(let ((A Alpha)
235	(B `(+ ,Alpha ,Beta)))
236	(let ((xexp (string->symbol (s+ x1 'exp))))
237	`(let ((,xexp (exp (* (neg ,B) ,dt))))
238	(+ (* ,x1 ,xexp) (* (- 1 ,xexp) (/ ,A ,B)))))))
239
240	(('let bnds body)
241	`(let ,bnds ,(expeuler dt name body)))
242
243	(else (nemo:error 'nemo:expeuler ": unable to rewrite equation " rhs
244	"in exponential Euler form")))))
245
246	res))
247
248
249
250	(define (state-init n init)
251	(let* ((init (rhsexpr/MATLAB init))
252	(init1 (canonicalize-expr/MATLAB init)))
253	(list (matlab-name n) init1)))
254
255
256	(define (asgn-eq n rhs)
257	(let* ((fbody (rhsexpr/MATLAB rhs))
258	(fbody1 (canonicalize-expr/MATLAB fbody)))
259	(list (matlab-name n) fbody1)))
260
261
262	(define (reaction-eq n open transitions conserve)
263	(match-let (((g cnode node-subs) (transitions-graph n open transitions conserve matlab-state-name)))
264	(let ((nodes ((g 'nodes))))
265	; (if (< 2 (length nodes))
266	; (list (matlab-name n) `(abs (/ ,(matlab-state-name n open)
267	; ,(sum (filter-map (lambda (x) (and (not (= (first x) (first cnode)))
268	; (matlab-state-name n (second x)))) nodes)))))
269	(list (matlab-name n) (matlab-state-name n open)))))
270
271
272	(define (reaction-transition-eqs n initial open transitions conserve power method)
273	(match-let (((g cnode node-subs) (transitions-graph n open transitions conserve matlab-state-name)))
274	(let* ((out-edges (g 'out-edges))
275	(in-edges (g 'in-edges))
276	(nodes ((g 'nodes))))
277	;; generate differential equations for each state in the transitions system
278	(let ((eqs (fold (lambda (s ax)
279	(if (and cnode (= (first cnode) (first s) )) ax
280	(let* ((out (out-edges (first s)))
281	(in (in-edges (first s)))
282	(open? (eq? (second s) open))
283	(name (matlab-name (lookup-def (second s) node-subs))))
284	(let* ((rhs1 (cond ((and (not (null? out)) (not (null? in)))
285	`(+ (neg ,(sum (map third out)))
286	,(sum (map third in))))
287	((and (not (null? out)) (null? in))
288	`(neg ,(sum (map third out))))
289	((and (null? out) (not (null? in)))
290	(sum (map third in)))))
291	(fbody0 (rhsexpr/MATLAB rhs1))
292	(fbody1 (case method
293	((expeuler) (canonicalize-expr/MATLAB (expeuler 'dt name fbody0)))
294	(else (canonicalize-expr/MATLAB fbody0)))))
295	(cons (list name fbody1) ax))
296	)))
297	(list) nodes)))
298	eqs))))
299
300
301
302
303	(define (poset->asgn-eq-defs poset sys)
304	(fold-right
305	(lambda (lst ax)
306	(fold (lambda (x ax)
307	(match-let (((i . n) x))
308	(let ((en (environment-ref sys n)))
309	(if (nemo:quantity? en)
310	(cases nemo:quantity en
311	(ASGN (name value rhs) (cons (asgn-eq name rhs) ax))
312	(else ax))
313	ax))))
314	ax lst))
315	(list) poset))
316
317	(define (poset->rate-eq-defs poset sys method)
318	(fold-right
319	(lambda (lst ax)
320	(fold (lambda (x ax)
321	(match-let (((i . n) x))
322	(let ((en (environment-ref sys n)))
323	(if (nemo:quantity? en)
324	(cases nemo:quantity en
325
326	(REACTION (name initial open transitions conserve power)
327	(append (reaction-transition-eqs name initial open transitions
328	conserve power method) ax))
329
330	(RATE (name initial rhs power)
331	(let ((fbody0 (rhsexpr/MATLAB rhs))
332	(dy (matlab-name name) ))
333	(case method
334	((expeuler)
335	(cons (list dy (canonicalize-expr/MATLAB (expeuler 'dt name fbody0)))
336	ax))
337	(else
338	(cons (list dy (canonicalize-expr/MATLAB fbody0)) ax)))))
339
340	(else ax))
341	ax))))
342	ax lst))
343	(list) poset))
344
345
346	(define (poset->reaction-eq-defs poset sys)
347	(fold-right
348	(lambda (lst ax)
349	(fold (lambda (x ax)
350	(match-let (((i . n) x))
351	(let ((en (environment-ref sys n)))
352	(if (nemo:quantity? en)
353	(cases nemo:quantity en
354	(REACTION (name initial open transitions conserve power)
355	(cons (reaction-eq name open transitions conserve) ax))
356	(else ax))
357	ax))))
358	ax lst))
359	(list) poset))
360
361
362	(define (poset->init-defs poset sys)
363	(fold-right
364	(lambda (lst ax)
365	(fold-right
366	(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	(REACTION (name initial open transitions conserve power)
372	(if (nemo:rhs? initial)
373	(cons* (state-init name initial)
374	(state-init (matlab-state-name name open) name) ax)
375	ax))
376
377	(RATE (name initial rhs power)
378	(if (and initial (nemo:rhs? initial))
379	(cons (state-init name initial) ax)
380	ax))
381
382	(else ax))
383	ax))))
384	ax lst))
385	(list) poset))
386
387
388	(define (poset->state-conserve-eq-defs poset sys)
389	(fold-right
390	(lambda (lst ax)
391	(fold (lambda (x ax)
392	(match-let (((i . n) x))
393	(let ((en (environment-ref sys n)))
394	(if (nemo:quantity? en)
395	(cases nemo:quantity en
396	(REACTION (name initial open transitions conserve power)
397	(if (and (list? conserve) (every nemo:conseq? conserve))
398	(cons (state-conseqs (matlab-name name) transitions conserve
399	matlab-state-name) ax)
400	ax))
401	(else ax))
402	ax))))
403	ax lst))
404	(list) poset))
405
406
407	(define (rate/reaction-power sys n)
408	(let ((en (environment-ref sys n)))
409	(if (nemo:quantity? en)
410	(cases nemo:quantity en
411	(REACTION (name initial open transitions conserve power) power)
412	(RATE (name initial rhs power) power)
413	(else #f)) #f)))
414
415
416	(define (bucket-partition p lst)
417	(let loop ((lst lst) (ax (list)))
418	(if (null? lst) ax
419	(let ((x (car lst)))
420	(let bkt-loop ((old-bkts ax) (new-bkts (list)))
421	(if (null? old-bkts) (loop (cdr lst) (cons (list x) new-bkts))
422	(if (p x (caar old-bkts ))
423	(loop (cdr lst) (append (cdr old-bkts) (cons (cons x (car old-bkts)) new-bkts)))
424	(bkt-loop (cdr old-bkts) (cons (car old-bkts) new-bkts)))))))))
425
426
427	(define (make-define-fn)
428	(lambda (indent n proc)
429	(let ((lst (procedure-data proc))
430	(indent+ (+ 2 indent)))
431	(let ((retval (matlab-name (gensym "retval")))
432	(rt (lookup-def 'rt lst))
433	(formals (lookup-def 'formals lst))
434	(vars (lookup-def 'vars lst))
435	(body (lookup-def 'body lst)))
436	(pp indent ,nl (function ,retval = ,(matlab-name n) (,(slp ", " vars)) ))
437	(let* ((body0 (rhsexpr/MATLAB body))
438	(body1 (canonicalize-expr/MATLAB body0))
439	(lbs (enum-bnds body1 (list))))
440	(pp indent+ ,(expr->string/MATLAB body1 retval))
441	(pp indent end))
442	))))
443
444
445	(define (output-dy sysname method globals state-index-map
446	rate-eq-defs reaction-eq-defs asgn-eq-defs
447	pool-ions mcap i-eqs v-eq indent indent+)
448
449	(pp indent ,nl (function dy = ,sysname (,(slp ", " (case method
450	((lsode) `(y t))
451	((odepkg) `(t y))
452	(else `(y t)))))))
453
454	(if (not (null? globals)) (pp indent+ (global ,(slp " " globals))))
455
456	(if (member 'v globals)
457	(let ((vi (lookup-def 'v state-index-map)))
458	(if vi (pp indent+ ,(expr->string/MATLAB (s+ "y(" vi ")") 'v)))))
459
460	(for-each (lambda (def)
461	(let ((n (first def)) )
462	(pp indent+
463	,(s+ "% rate eq " n)
464	,(expr->string/MATLAB
465	(s+ "y(" (lookup-def n state-index-map) ")") n))))
466	rate-eq-defs)
467
468	(let* ((eqs
469	(append
470
471	asgn-eq-defs
472	reaction-eq-defs
473
474	(map (lambda (pool-ion)
475	(let ((n (third pool-ion))
476	(b (first pool-ion)))
477	(list n b)))
478	pool-ions)))
479
480	(eq-dag
481	(map (lambda (def)
482	(cons (first def) (enum-freevars (second def) '() '())))
483	eqs))
484	(eq-order
485	(reverse
486	(topological-sort eq-dag
487	(lambda (x y) (string=? (->string x) (->string y)))))))
488	(for-each (lambda (n)
489	(let ((b (lookup-def n eqs)))
490
491	(if b (pp indent+
492	,(s+ "% equation for " n)
493	,(expr->string/MATLAB b (matlab-name n))))))
494	eq-order))
495
496
497
498	(pp indent+ ,(expr->string/MATLAB `(zeros (length y) 1) 'dy))
499
500	(for-each (lambda (def)
501	(let ((n (s+ "dy(" (lookup-def (first def) state-index-map) ")"))
502	(b (second def)))
503	(pp indent+ ,(s+ "% state " (first def))
504	,(expr->string/MATLAB b n))))
505	rate-eq-defs)
506
507	(for-each (lambda (def)
508	(pp indent+ ,(expr->string/MATLAB (second def) (first def))))
509	i-eqs)
510
511	(if v-eq
512	(let ((vi (lookup-def 'v state-index-map)))
513	(if vi
514	(pp indent+ ,(expr->string/MATLAB (second v-eq) (s+ "dy(" vi ")"))))))
515
516	(pp indent ,nl (end))
517	)
518
519	(define (output-steadystate sysname method globals steady-state-index-map pool-ions
520	const-defs init-eq-defs asgn-eq-defs rate-eq-defs reaction-eq-defs
521	indent indent+)
522
523	(if (not (null? steady-state-index-map))
524	(begin
525	(pp indent ,nl (function y = ,(s+ sysname "_steadystate") (x)))
526
527	(if (not (null? globals)) (pp indent+ (global ,(slp " " globals))))
528
529	(for-each
530	(lambda (def)
531	(let* ((n (first def))
532	(ni (lookup-def n steady-state-index-map)))
533	(if ni (pp indent+ ,(expr->string/MATLAB (s+ "x(" ni ")") n)))
534	))
535	rate-eq-defs)
536
537	(pp indent+ ,(expr->string/MATLAB `(zeros ,(length steady-state-index-map) 1) 'y))
538
539
540	(let* ((init-eqs
541	(append
542
543	const-defs
544	asgn-eq-defs
545	init-eq-defs
546
547	(map (lambda (pool-ion)
548	(let ((n (third pool-ion))
549	(b (first pool-ion)))
550	(list n b)))
551	pool-ions)))
552
553	(init-dag
554	(map (lambda (def)
555	(cons (first def) (enum-freevars (second def) '() '())))
556	init-eqs))
557	(init-order
558	(reverse
559	(topological-sort init-dag
560	(lambda (x y) (string=? (->string x) (->string y)))))))
561	(for-each (lambda (n)
562	(let ((b (lookup-def n init-eqs)))
563	(if b (pp indent+ ,(expr->string/MATLAB b (matlab-name n))))))
564	init-order))
565
566	(for-each
567	(lambda (def)
568	(let* ((n (first def))
569	(ni (lookup-def n steady-state-index-map))
570	(b (second def)))
571	(if ni (pp indent+ ,(s+ "% state " n)
572	,(expr->string/MATLAB b (s+ "y(" ni ")"))))
573	))
574	rate-eq-defs)
575
576	(pp indent ,nl (end))))
577	)
578
579
580	(define (output-print-state sysname state-index-map indent indent+)
581	(pp indent ,nl (function ,(s+ sysname "_print_state") (y)))
582
583	(let ((lst (sort (map (lambda (x) (cons (->string (car x)) (cdr x))) state-index-map)
584	(lambda (x y) (string<? (car x) (car y))))))
585	(for-each (lambda (p)
586	(let ((n (first p)) (i (second p)))
587	(pp indent+ (,n " = " "y(" ,i ")"))))
588	lst))
589
590	(pp indent ,nl (end))
591	)
592
593
594	(define (output-init sysname globals state-index-map steady-state-index-map
595	const-defs asgn-eq-defs init-eq-defs rate-eq-defs
596	reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
597
598	(pp indent ,nl (function y0 = ,(s+ sysname "_init") (Vinit)))
599
600	(if (not (null? globals)) (pp indent+ (global ,(slp " " globals))))
601
602	(pp indent+ ,(expr->string/MATLAB `(zeros ,(length state-index-map) 1) 'y0))
603
604	(if (member 'v globals)
605	(let ((vi (lookup-def 'v state-index-map)))
606	(pp indent+ ,(expr->string/MATLAB `Vinit 'v))
607	(if vi (pp indent+ ,(expr->string/MATLAB 'v (s+ "y0(" vi ")") )))))
608
609
610	(let* ((init-eqs
611	(append
612
613	const-defs
614	asgn-eq-defs
615	init-eq-defs
616
617	(map (lambda (pool-ion)
618	(let ((n (third pool-ion))
619	(b (first pool-ion)))
620	(list n b)))
621	pool-ions)))
622
623	(init-dag
624	(map (lambda (def)
625	(cons (first def) (enum-freevars (second def) '() '())))
626	init-eqs))
627	(init-order
628	(reverse
629	(topological-sort init-dag
630	(lambda (x y) (string=? (->string x) (->string y)))))))
631
632	(for-each (lambda (n)
633	(let ((b (lookup-def n init-eqs)))
634	(if b (pp indent+ ,(expr->string/MATLAB b (matlab-name n))))))
635	init-order))
636
637	(for-each (lambda (def)
638	(let* ((n (first def))
639	(ni (lookup-def n state-index-map)))
640	(if ni (pp indent+ ,(expr->string/MATLAB n (s+ "y0(" ni ")"))))))
641	init-eq-defs)
642
643	(if (not (null? steady-state-index-map))
644	(begin
645	(pp indent+ ,(expr->string/MATLAB `(zeros ,(length steady-state-index-map) 1) 'xs)
646	,(expr->string/MATLAB
647	`(fsolve ,(s+ #\' sysname "_steadystate" #\' ) xs)
648	"[ys,fval,info]"))
649
650
651	(for-each
652	(lambda (def)
653	(let* ((n (first def))
654	(si (lookup-def n steady-state-index-map))
655	(ni (lookup-def n state-index-map)))
656	(if si (begin
657	(pp indent+ ,(expr->string/MATLAB (s+ "ys(" si ")") n))
658	(pp indent+ ,(expr->string/MATLAB (s+ "ys(" si ")") (s+ "y0(" ni ")")))))))
659	rate-eq-defs)
660
661
662	))
663
664	(for-each (lambda (def)
665	(let ((n (first def)) (b (second def)))
666	(if (not (lookup-def n init-eq-defs))
667	(pp indent+ ,(expr->string/MATLAB b n)))))
668	reaction-eq-defs)
669
670	(for-each (lambda (def)
671	(pp indent+ ,(expr->string/MATLAB (second def) (first def))))
672	i-eqs)
673
674	(pp indent ,nl (end))
675
676	)
677
678	(define (matlab-translator1 sys . rest)
679	(define (cid x) (second x))
680	(define (cn x) (first x))
681	(let-optionals rest ((mode 'multiple) (method 'lsode) (filename #f))
682	(match-let ((($ nemo:quantity 'DISPATCH dis) (environment-ref sys (nemo-intern 'dispatch))))
683	(let ((imports ((dis 'imports) sys))
684	(exports ((dis 'exports) sys)))
685	(let* ((indent 0)
686	(indent+ (+ 2 indent ))
687	(eval-const (dis 'eval-const))
688	(sysname (matlab-name ((dis 'sysname) sys)))
689	(prefix sysname)
690	(filename (or filename (s+ sysname ".m")))
691	(deps* ((dis 'depgraph*) sys))
692	(consts ((dis 'consts) sys))
693	(asgns ((dis 'asgns) sys))
694	(states ((dis 'states) sys))
695	(reactions ((dis 'reactions) sys))
696	(defuns ((dis 'defuns) sys))
697	(components ((dis 'components) sys))
698
699	(g (match-let (((state-list asgn-list g) ((dis 'depgraph*) sys))) g))
700	(poset (vector->list ((dis 'depgraph->bfs-dist-poset) g)))
701
702	(const-defs (filter-map
703	(lambda (nv)
704	(and (not (member (first nv) matlab-builtin-consts))
705	(let ((v1 (canonicalize-expr/MATLAB (second nv))))
706	(list (matlab-name (first nv)) v1))))
707	consts))
708
709	(gate-complex-info (nemo:gate-complex-query sys))
710
711	(gate-complexes (lookup-def 'gate-complexes gate-complex-info))
712	(perm-ions (map (match-lambda ((comp i e erev) `(,comp ,(matlab-name i) ,(matlab-name e) ,erev)))
713	(lookup-def 'perm-ions gate-complex-info)))
714	(acc-ions (map (match-lambda ((comp i in out) `(,comp ,@(map matlab-name (list i in out)))))
715	(lookup-def 'acc-ions gate-complex-info)))
716	(epools (lookup-def 'pool-ions gate-complex-info))
717	(pool-ions (map (lambda (lst) (map matlab-name lst)) epools))
718
719	(i-gates (lookup-def 'i-gates gate-complex-info))
720
721	(capcomp (any (match-lambda ((name 'membrane-capacitance id) (list name id)) (else #f)) components))
722	(mcap (and capcomp (car ((dis 'component-exports) sys (cid capcomp)))))
723
724
725	(i-eqs (filter-map
726	(lambda (gate-complex)
727
728	(let* ((label (first gate-complex))
729	(n (second gate-complex))
730	(subcomps ((dis 'component-subcomps) sys n))
731	(acc (lookup-def 'accumulating-substance subcomps))
732	(perm (lookup-def 'permeating-ion subcomps))
733	(permqs (and perm ((dis 'component-exports) sys (cid perm))))
734	(pore (lookup-def 'pore subcomps))
735	(permeability (lookup-def 'permeability subcomps))
736	(gate (lookup-def 'gate subcomps))
737	(sts (and gate ((dis 'component-exports) sys (cid gate)))))
738
739	(if (not (or pore permeability))
740	(nemo:error 'nemo:matlab-translator ": ion channel definition " label
741	"lacks any pore or permeability components"))
742
743
744	(cond ((and perm permeability gate)
745	(let* ((i (matlab-name (s+ 'i (cn perm))))
746	(pmax (car ((dis 'component-exports) sys (cid permeability))))
747	(pwrs (map (lambda (n) (rate/reaction-power sys n)) sts))
748	(sptms (map (lambda (st pwr) (if pwr `(pow ,st ,pwr) st)) sts pwrs))
749	(gion `(* ,pmax ,@sptms)))
750	(list i #f gion (matlab-name (s+ 'i_ label) ))))
751
752	((and perm pore gate)
753	(case (cn perm)
754	((non-specific)
755	(let* ((i (matlab-name 'i))
756	(e (car permqs))
757	(gmax (car ((dis 'component-exports) sys (cid pore))))
758	(pwrs (map (lambda (n) (rate/reaction-power sys n)) sts))
759	(sptms (map (lambda (st pwr) (if pwr `(pow ,st ,pwr) st)) sts pwrs))
760	(gion `(* ,gmax ,@sptms)))
761	(list i e gion (matlab-name (s+ 'i_ label) ))))
762
763	(else
764	(let* ((i (matlab-name (s+ 'i (cn perm))))
765	(e (car permqs))
766	(gmax (car ((dis 'component-exports) sys (cid pore))))
767	(pwrs (map (lambda (n) (rate/reaction-power sys n)) sts))
768	(sptms (map (lambda (st pwr) (if pwr `(pow ,st ,pwr) st)) sts pwrs))
769	(gion `(* ,gmax ,@sptms)))
770	(list i e gion (matlab-name (s+ 'i_ label) ))))))
771
772	((and perm pore)
773	(case (cn perm)
774	((non-specific)
775	(let* ((i (matlab-name 'i))
776	(e (car permqs))
777	(gmax (car ((dis 'component-exports) sys (cid pore)))))
778	(list i e gmax (matlab-name (s+ 'i_ label) ))))
779	(else
780	(nemo:error 'nemo:matlab-translator ": invalid ion channel definition " label))))
781
782	((and acc pore gate)
783	(let* ((i (matlab-name (s+ 'i (cn acc))))
784	(gmax (car ((dis 'component-exports) sys (cid pore))))
785	(pwrs (map (lambda (n) (rate/reaction-power sys n)) sts))
786	(sptms (map (lambda (st pwr) (if pwr `(pow ,st ,pwr) st)) sts pwrs))
787	(gion `(* ,gmax ,@sptms)))
788	(list i #f gion (matlab-name (s+ 'i_ label) ))))
789	(else (nemo:error 'nemo:matlab-translator ": invalid ion channel definition " label))
790	)))
791	gate-complexes))
792
793	(i-names (delete-duplicates (map first i-eqs)))
794
795	(i-eqs (fold (lambda (i-gate ax)
796	(let ((i-gate-var (first i-gate)))
797	(cons (list (matlab-name 'i) #f i-gate-var (s+ 'i_ (second i-gate))) ax)))
798	i-eqs i-gates))
799
800	(i-bkts (bucket-partition (lambda (x y) (eq? (car x) (car y))) i-eqs))
801
802	(i-eqs (fold (lambda (b ax)
803	(match b
804	((and ps ((i e gion ii) . rst))
805	(let loop ((ps ps) (summands (list)) (eqs (list)))
806	(if (null? ps)
807
808	(let* ((sum0 (sum summands))
809	(sum1 (rhsexpr/MATLAB sum0))
810	(sum2 (canonicalize-expr/MATLAB sum1)))
811	(append eqs (list (list i sum2)) ax))
812
813	(match-let (((i e gion ii) (car ps)))
814	(loop (cdr ps)
815	(cons ii summands)
816	(let* ((expr0 (rhsexpr/MATLAB (if e `(* ,gion (- v ,e)) gion)))
817	(expr1 (canonicalize-expr/MATLAB expr0)))
818	(cons (list ii expr1) eqs)))))))
819
820	((i e gion ii)
821	(let* ((expr0 (rhsexpr/MATLAB (if e `(* ,gion (- v ,e)) gion)))
822	(expr1 (canonicalize-expr/MATLAB expr0)))
823	(cons (list i expr1) ax)))
824
825	(else ax)))
826	(list) i-bkts))
827
828	(asgn-eq-defs (poset->asgn-eq-defs poset sys))
829
830	(rate-eq-defs (reverse (poset->rate-eq-defs poset sys method)))
831
832	(reaction-eq-defs (poset->reaction-eq-defs poset sys))
833
834	(init-eq-defs (poset->init-defs poset sys))
835
836	(conserve-eq-defs (map (lambda (eq) (list 0 `(- ,(second eq) ,(first eq))))
837	(poset->state-conserve-eq-defs poset sys)))
838
839	(globals (map matlab-name
840	(delete-duplicates (append
841	exports
842	(map second perm-ions)
843	(map third perm-ions)
844	(map second acc-ions)
845	(map third acc-ions)
846	(map fourth acc-ions)
847	(map second pool-ions)
848	(map third pool-ions)
849	(map first imports)
850	(map first const-defs)
851	(map first asgn-eq-defs)
852	(map (lambda (gate-complex) (matlab-name (s+ 'i_ (first gate-complex)))) gate-complexes )
853	(map (lambda (i-gate) (matlab-name (s+ 'i_ (second i-gate)))) i-gates )
854
855	))))
856	(v-eq (if (and mcap (member 'v globals))
857	(list 'v (rhsexpr/MATLAB `(/ (neg ,(sum i-names)) ,mcap)))
858	(list 'v 0.0)))
859
860	(state-index-map (let ((acc (fold (lambda (def ax)
861	(let ((st-name (first def)))
862	(list (+ 1 (first ax))
863	(cons `(,st-name ,(first ax)) (second ax)))))
864	(list 1 (list))
865	(if (member 'v globals)
866	(cons (list 'v) rate-eq-defs)
867	rate-eq-defs)
868	)))
869	(second acc)))
870
871	(steady-state-index-map (let ((acc (fold (lambda (def ax)
872	(let ((st-name (first def)))
873	(if (not (alist-ref st-name init-eq-defs))
874	(list (+ 1 (first ax))
875	(cons `(,st-name ,(first ax)) (second ax)))
876	ax)))
877	(list 1 (list))
878	rate-eq-defs)))
879	(second acc)))
880
881	(dfenv
882	(map (lambda (x) (let ((n (first x)))
883	(list n (matlab-name (s+ "d_" n )))))
884	defuns))
885
886	(output-globals
887	(lambda () (if (not (null? globals)) (pp indent (global ,(slp " " globals))))))
888	)
889
890
891
892	(for-each
893	(lambda (a)
894	(let ((acc-ion (car a)))
895	(if (assoc acc-ion perm-ions)
896	(nemo:error 'nemo:matlab-translator
897	": ion species " acc-ion " cannot be declared as both accumulating and permeating"))))
898	acc-ions)
899
900	(for-each
901	(lambda (p)
902	(let ((pool-ion (car p)))
903	(if (assoc pool-ion perm-ions)
904	(nemo:error 'nemo:matlab-translator
905	": ion species " pool-ion " cannot be declared as both pool and permeating"))))
906	pool-ions)
907
908	(let ((output (case mode
909	((single) (open-output-file filename))
910	(else #f))))
911
912	(if output (with-output-to-port output (lambda () (output-globals))))
913
914	;; derivative function
915	(let ((output1 (or output (open-output-file (s+ prefix "_dy.m")))))
916	(with-output-to-port output1
917	(lambda ()
918	(output-dy sysname method globals state-index-map
919	rate-eq-defs reaction-eq-defs asgn-eq-defs
920	pool-ions mcap i-eqs v-eq
921	indent indent+)))
922	(if (not output) (close-output-port output1)))
923
924
925	;; steady state solver
926	(let ((output1 (or output (open-output-file (s+ prefix "_steadystate.m")))))
927	(with-output-to-port output1
928	(lambda ()
929	(output-steadystate sysname method globals steady-state-index-map pool-ions
930	const-defs init-eq-defs asgn-eq-defs rate-eq-defs reaction-eq-defs
931	indent indent+)))
932	(if (not output) (close-output-port output1)))
933
934	;; initial values function
935	(let ((output1 (or output (open-output-file (s+ prefix "_init.m")))))
936	(with-output-to-port output1
937	(lambda ()
938	(output-init sysname globals state-index-map steady-state-index-map
939	const-defs asgn-eq-defs init-eq-defs rate-eq-defs
940	reaction-eq-defs i-eqs pool-ions perm-ions indent indent+)
941	(pp indent ,nl)))
942	(if (not output) (close-output-port output1)))
943
944	;; user-defined functions
945	(let* (;;(with (inexact->exact (round (/ (abs (- max-v min-v)) step))))
946	(define-fn (make-define-fn)))
947	(for-each (lambda (fndef)
948	(if (not (member (car fndef) builtin-fns))
949	(let ((output1 (or output (open-output-file (s+ (matlab-name (car fndef)) ".m")))))
950	(with-output-to-port output1
951	(lambda ()
952	(apply define-fn (cons indent fndef))
953	(pp indent ,nl)))
954	(if (not output) (close-output-port output1)))))
955	defuns))
956
957	(if output (close-output-port output)))
958
959	)))))
960
961
962	(define (nemo:matlab-translator sys . rest)
963	(apply matlab-translator1 (cons sys (cons 'multiple (cons 'odepkg rest)))))
964
965	(define (nemo:octave-translator sys . rest)
966	(apply matlab-translator1 (cons sys (cons 'single rest))))
967
968	#\|
969
970	(define (make-define-dfn fenv)
971
972	(define (D vars body)
973	(let loop ((vars vars) (exprs (list)))
974	(if (null? vars) (if (pair? (cdr exprs))
975	`(+ . ,(reverse exprs)) (car exprs))
976	(let ((de (let ((de (differentiate fenv (car vars) body)))
977	(let loop ((e (simplify de)) (de de))
978	(if (equal? de e) de
979	(loop (simplify e) e))))))
980	(loop (cdr vars) (cons de exprs))))))
981
982	(lambda (indent n proc)
983	(let ((lst (procedure-data proc))
984	(indent+ (+ 2 indent)))
985	(let ((retval (matlab-name (gensym "retval")))
986	(rt (lookup-def 'rt lst))
987	(formals (lookup-def 'formals lst))
988	(vars (lookup-def 'vars lst))
989	(body (lookup-def 'body lst)))
990	(pp indent ,nl (function ,retval = ,(s+ "d_" (matlab-name n) ) (,(slp ", " vars)) ))
991	(let* ((dbody (D vars body))
992	(dbody (canonicalize-expr/MATLAB (rhsexpr/MATLAB dbody))))
993	(pp indent+ ,(expr->string/MATLAB dbody retval))
994	(pp indent end))
995	)))
996	)
997
998	(define (output-jac sysname method globals fenv state-index-map
999	rate-eq-defs reaction-eq-defs asgn-eq-defs
1000	mcap i-eqs v-eq
1001	defuns indent indent+)
1002
1003
1004	(define (D var expr)
1005	(let* ((dexpr (distribute expr))
1006	(de (differentiate fenv var dexpr))
1007	(se (let loop ((e (simplify de)) (de de))
1008	(if (equal? de e) de
1009	(loop (simplify e) e)))))
1010	se))
1011
1012	(pp indent ,nl (function jac = ,(s+ sysname "_jac")
1013	(,(slp ", " `(t y)))))
1014
1015	(if (not (null? globals)) (pp indent+ (global ,(slp " " globals))))
1016
1017	(if (member 'v globals)
1018	(let ((vi (lookup-def 'v state-index-map)))
1019	(if vi (pp indent+ ,(expr->string/MATLAB (s+ "y(" vi ")") 'v)))))
1020
1021	(for-each (lambda (def)
1022	(let ((name (first def)))
1023	(pp indent+ ,(expr->string/MATLAB
1024	(s+ "y(" (lookup-def name state-index-map) ")") name))))
1025	rate-eq-defs)
1026
1027	(for-each (lambda (def)
1028	(let ((n (first def)) (b (second def)))
1029	(pp indent+
1030	,(expr->string/MATLAB b n))))
1031	reaction-eq-defs)
1032
1033	(pp indent+ ,(expr->string/MATLAB `(zeros (length y) (length y)) 'jac))
1034
1035
1036	(for-each
1037	(lambda (def)
1038
1039	(let* ((name (first def))
1040	(rhs (second def))
1041	(fv (enum-freevars rhs '() '() ))
1042	(drvs (map (lambda (st)
1043	(let ((var (first st))
1044	(expr `(let ,(or (and v-eq (list v-eq)) '())
1045	(let ,rate-eq-defs
1046	(let ,reaction-eq-defs
1047	(let ,(filter-map (lambda (x)
1048	(or (assoc x i-eqs)
1049	(assoc x asgn-eq-defs)))
1050	fv))
1051	,rhs)))))
1052	(D var expr)))
1053	state-index-map))
1054	(vs (map (lambda (st)
1055	(let ((ni (lookup-def name state-index-map))
1056	(i (second st)))
1057	(s+ "jac(" (slp "," (list ni i)) ")")))
1058	state-index-map))
1059	(ndrvs (map (lambda (x) (canonicalize-expr/MATLAB (rhsexpr/MATLAB x))) drvs)))
1060
1061	(pp indent+ ,(s+ "% state " name))
1062	(for-each (lambda (ndrv v)
1063	(pp indent+ ,(expr->string/MATLAB ndrv v)))
1064	ndrvs vs)
1065	))
1066
1067	(if (and mcap (member 'v globals))
1068	(let ((vi (lookup-def 'v state-index-map)))
1069	(if vi (cons
1070	(list 'v `(/ (neg ,(sum (map first i-eqs))) ,mcap))
1071	rate-eq-defs)
1072	rate-eq-defs))
1073	rate-eq-defs))
1074
1075	(pp indent ,nl (end))
1076
1077	)
1078
1079	\|#
1080
1081	)

Note: See TracBrowser for help on using the repository browser.

Download in other formats: