source: project/release/4/9ML-toolkit/trunk/ivp-chicken.scm @ 29732

Last change on this file since 29732 was 29732, checked in by Ivan Raikov, 8 years ago

9ML-toolkit: bug fixes in propagation of initial states of primitives

File size: 12.2 KB
Line 
1;;
2;;  NineML IVP code generator for Chicken Scheme.
3;;
4;;
5;; Copyright 2010-2013 Ivan Raikov and the Okinawa Institute of
6;; Science and Technology.
7;;
8;; This program is free software: you can redistribute it and/or
9;; modify it under the terms of the GNU General Public License as
10;; published by the Free Software Foundation, either version 3 of the
11;; License, or (at your option) any later version.
12;;
13;; This program is distributed in the hope that it will be useful, but
14;; WITHOUT ANY WARRANTY; without even the implied warranty of
15;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
16;; General Public License for more details.
17;;
18;; A full copy of the GPL license can be found at
19;; <http://www.gnu.org/licenses/>.
20;;
21
22
23(module 9ML-ivp-chicken
24
25        (ivp-chicken ivp-chicken/cvode ivp-chicken-codegen ivp-chicken-codegen/cvode)
26
27        (import scheme chicken )
28
29(import (only files make-pathname pathname-directory pathname-file)
30        (only data-structures conc))
31(require-extension make datatype signal-diagram 9ML-eval setup-api)
32
33
34(define nl "\n")
35       
36
37(define chicken-run
38#<<EOF
39
40(define-syntax run
41  (syntax-rules ()
42    ((_ f indep dep events end input parameters)
43     (let ((nstate input))
44       (printf "# ~A " indep)
45       (for-each (lambda (x) (printf "~A " x)) dep)
46       (printf "~%")
47       (let recur ((nstate nstate))
48         (let ((ival (alist-ref indep nstate)))
49           (printf "~A " ival)
50           (for-each (lambda (x)
51                       (let ((v (alist-ref x nstate)))
52                         (printf "~A " (if (boolean? v) (or (and v 1) 0) v))))
53                     dep)
54           (printf "~%")
55           (if (> ival end)
56               (print "# All done!")
57               (recur (append (f nstate) parameters))))))
58     )))
59EOF
60)
61
62
63(define chicken-run/cvode
64#<<EOF
65
66
67(define-syntax run
68  (syntax-rules ()
69    ((_ f h indep dep events end input parameters)
70     (let* ((nstate (make-parameter input))
71
72            (numstates (filter (lambda (y) (number? (alist-ref y input))) dep))
73           
74            (update-nstate (lambda (old new)
75                             (old (fold (lambda (kv ax) (alist-update (car kv) (cdr kv) ax)) (old) new))))
76
77            (vector->nstate (lambda (yvec) 
78                              (fold (lambda (y yv nstate) (cons (cons y yv) nstate))
79                                    '()
80                                    numstates (f64vector->list yvec))))
81            (nstate->vector (lambda (nstate) 
82                              (list->f64vector (map (lambda (y) (alist-ref y nstate)) numstates))))
83
84            (fwrap (lambda (t yvec nstate)
85                     (let ((nstate1 (cons (cons indep t) (vector->nstate yvec))))
86                       (nstate nstate1)
87                       (let ((nstate2 (f (append nstate1 parameters))))
88                         (nstate->vector nstate2)
89                         ))
90                     ))
91
92            (fcall (lambda (t yvec nstate)
93                     (let ((nstate1 (cons (cons indep t) (vector->nstate yvec))))
94                       (nstate nstate1)
95                       (let ((nstate2 (f (append nstate1 parameters))))
96                         nstate2
97                         ))
98                     ))
99
100            (event-detect (lambda (t yy nstate)
101                            (let ((Vpeak (alist-ref 'Vpeak parameters))
102                                  (V (f64vector-ref yy 2)))
103                              (s32vector (floor (- V Vpeak))))))
104
105
106            (solver (cvode-create-solver
107                     (alist-ref indep input) 
108                     (nstate->vector input)
109                     fwrap
110                     abstol: 1e-4
111                     reltol: 1e-4
112                     events: (make-s32vector (length events) 0)
113                     event-fn: event-detect
114                     user-data: nstate
115                     ))
116
117            (solvewrap (lambda (solver tnext)
118                         (let ((status (cvode-solve solver tnext)))
119                           (if (negative? status) (error "CVODE solver error" status))
120                           (if (zero? status)
121                               (cons (cons indep tnext) (vector->nstate (cvode-yy solver)))
122                               (let ((yy (cvode-yy solver)) (tt (cvode-t solver)))
123                                 (nstate (f (append (fcall tt yy nstate) parameters)))
124                                 (cvode-reinit-solver solver tnext (nstate->vector (nstate)))
125                                 (cons (cons indep tt) (vector->nstate (cvode-yy solver)))
126                                 )))
127                           ))
128            )
129       (printf "# ~A " indep)
130       (for-each (lambda (x) (printf "~A " x)) dep)
131       (printf "~%")
132       (let recur ((hv (alist-ref h parameters)))
133         (let ((ival (alist-ref indep (nstate))))
134           (printf "~A " ival)
135           (for-each (lambda (x)
136                       (let ((v (alist-ref x (nstate))))
137                         (printf "~A " (if (boolean? v) (or (and v 1) 0) v))))
138                     dep)
139           (printf "~%")
140           (if (> ival end)
141               (print "# All done!")
142               (let ((nstate1 (solvewrap solver (+ hv ival))))
143                 (recur hv)))
144           ))
145       ))
146    ))
147
148EOF
149)
150
151(define (ivp-chicken prefix ivp-id ivar dvars pvars events start end ic sd
152                     #!key (solver 'rk3))
153
154  (let* ((dir (or (pathname-directory prefix) "."))
155         (solver-path (make-pathname dir (conc ivp-id "_solver.scm")))
156         (run-path    (make-pathname dir (sprintf "~A_run.scm" ivp-id)))
157         (exec-path   (make-pathname dir (sprintf "~A_run" ivp-id)))
158         (log-path    (make-pathname dir (sprintf "~A_~A.log" (pathname-file prefix) ivp-id)))
159         (csc-path    (make-pathname (program-path) "csc")))
160   
161    (make 
162        (
163         (solver-path (prefix)
164                      (with-output-to-file solver-path
165                        (lambda () (codegen/scheme ivp-id sd solver: solver))))
166         
167         (run-path (prefix)
168                   (with-output-to-file run-path
169                     (lambda () 
170                       (print-fragments
171                        (list
172                         (sprintf "(include \"~A_solver.scm\")~%~%" ivp-id)
173                         chicken-run nl
174                         (sprintf "(define initial (quote ~A))~%~%" 
175                                  (cons (cons ivar start) 
176                                        (map (lambda (x)
177                                               (let ((n (car x))
178                                                     (v (cdr x)))
179                                                 (if (pair? v)
180                                                     (case (car v)
181                                                       ((generator) (sprintf "(~A)" (cadr v)))
182                                                       (else (error 'ivp-chicken-codegen "invalid initial value" v)))
183                                                     x)))
184                                             ic)))
185                         (sprintf "(define parameters (quote ~A))~%~%" 
186                                  (map (lambda (x) (assoc x ic)) pvars))
187                         (sprintf "(run ~A (quote ~A) (quote ~A) (quote ~A) ~A initial parameters)~%~%" 
188                                  ivp-id ivar dvars events end)
189                         ))
190                       ))
191                   )
192         
193         (exec-path (run-path solver-path)
194                    (run (,csc-path -w -I ,dir -b -S -d0 -O3 -disable-interrupts ,run-path)))
195         
196         (log-path (exec-path) (run (,exec-path > ,log-path)))
197         )
198     
199      (list log-path) )
200    ))
201
202
203(define (ivp-chicken/cvode prefix ivp-id hvar ivar dvars pvars events start end ic sd)
204
205  (let* ((dir (or (pathname-directory prefix) "."))
206         (solver-path (make-pathname dir (conc ivp-id "_solver.scm")))
207         (run-path    (make-pathname dir (sprintf "~A_run.scm" ivp-id)))
208         (exec-path   (make-pathname dir (sprintf "~A_run" ivp-id)))
209         (log-path    (make-pathname dir (sprintf "~A_~A.log" (pathname-file prefix) ivp-id)))
210         (csc-path    (make-pathname (program-path) "csc")))
211   
212    (make 
213        (
214         (solver-path (prefix)
215                      (with-output-to-file solver-path
216                        (lambda () (codegen/scheme ivp-id sd solver: 'cvode))))
217         
218         (run-path (prefix)
219                   (with-output-to-file run-path
220                     (lambda () 
221                       (print-fragments
222                        (list
223                         (sprintf "(include \"~A_solver.scm\")~%~%" ivp-id)
224                         chicken-run/cvode nl
225                         (sprintf "(define initial (quote ~A))~%~%" 
226                                  (cons (cons ivar start) 
227                                        (map (lambda (x)
228                                               (let ((n (car x))
229                                                     (v (cdr x)))
230                                                 (if (pair? v)
231                                                     (case (car v)
232                                                       ((generator) (sprintf "(~A)" (cadr v)))
233                                                       (else (error 'ivp-chicken-codegen "invalid initial value" v)))
234                                                     x)))
235                                             ic)))
236                         (sprintf "(define parameters (quote ~A))~%~%" 
237                                  (map (lambda (x) (assoc x ic)) pvars))
238                         (sprintf "(run ~A (quote ~A) (quote ~A) (quote ~A) (quote ~A) ~A initial parameters)~%~%" 
239                                  ivp-id hvar ivar dvars events end)
240                         )))))
241         
242         (exec-path (run-path solver-path)
243                    (run (,csc-path -w -I ,dir -b -S -d0 -O3 -disable-interrupts ,run-path)))
244         
245         (log-path (exec-path) (run (,exec-path > ,log-path)))
246         )
247     
248      (list log-path) )
249    ))
250
251
252(define (ivp-chicken-codegen prefix ivp-id ivar dvars pvars events ic sd
253                     #!key (solver 'rk3))
254
255  (let* ((dir (or (pathname-directory prefix) "."))
256         (solver-path (make-pathname dir (conc ivp-id "_solver.scm")))
257         (run-path    (make-pathname dir (sprintf "~A_run.scm" ivp-id)))
258         )
259   
260    (make 
261        (
262         (solver-path (prefix)
263                      (with-output-to-file solver-path
264                        (lambda () (codegen/scheme ivp-id sd solver: solver))))
265         
266         (run-path (prefix)
267                   (with-output-to-file run-path
268                     (lambda () 
269                       (print-fragments
270                        (list
271                         (sprintf "(include \"~A_solver.scm\")~%~%" ivp-id)
272                         chicken-run nl
273                         (sprintf "(define initial (quote ~A))~%~%" 
274                                  (cons (cons ivar 0.0) 
275                                        (map (lambda (x)
276                                               (let ((n (car x))
277                                                     (v (cdr x)))
278                                                 (if (pair? v)
279                                                     (case (car v)
280                                                       ((generator) (sprintf "(~A)" (cadr v)))
281                                                       (else (error 'ivp-chicken-codegen "invalid initial value" v)))
282                                                     x)))
283                                             ic)))
284                         (sprintf "(define parameters (quote ~A))~%~%" 
285                                  (map (lambda (x) (assoc x ic)) pvars))
286                         ))
287                       ))
288                   )
289         
290         )
291     
292      (list solver-path run-path) )
293    ))
294
295
296(define (ivp-chicken-codegen/cvode prefix ivp-id hvar ivar dvars pvars events ic sd)
297
298  (let* ((dir (or (pathname-directory prefix) "."))
299         (solver-path (make-pathname dir (conc ivp-id "_solver.scm")))
300         (run-path    (make-pathname dir (sprintf "~A_run.scm" ivp-id)))
301         )
302
303    (make 
304        (
305         (solver-path (prefix)
306                      (with-output-to-file solver-path
307                        (lambda () (codegen/scheme ivp-id sd solver: 'cvode))))
308         
309         (run-path (prefix)
310                   (with-output-to-file run-path
311                     (lambda () 
312                       (print-fragments
313                        (list
314                         (sprintf "(include \"~A_solver.scm\")~%~%" ivp-id)
315                         chicken-run/cvode nl
316                         (sprintf "(define initial (quote ~A))~%~%" 
317                                  (cons (cons ivar 0.0) 
318                                        (map (lambda (x)
319                                               (let ((n (car x))
320                                                     (v (cdr x)))
321                                                 (if (pair? v)
322                                                     (case (car v)
323                                                       ((generator) (sprintf "(~A)" (cadr v)))
324                                                       (else (error 'ivp-chicken-codegen "invalid initial value" v)))
325                                                     x)))
326                                             ic)))
327                         (sprintf "(define parameters (quote ~A))~%~%" 
328                                  (map (lambda (x) (assoc x ic)) pvars))
329                         )))))
330         )
331     
332      (list solver-path run-path) )
333    ))
334
335
336)
Note: See TracBrowser for help on using the repository browser.