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

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

9ML-toolkit: bringing octave and scheme backends up to date; changed Izhikevich FS example to use heaviside function

File size: 11.7 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(define (chicken-value v)
37  (cond
38   ((pair? v)
39    (case (car v)
40      ((realsig)   (chicken-value (caddr v)))
41      ((realconst) (chicken-value (cadr v)))
42      ((generator) (sprintf "~A ()" (cadr v)))
43      ((random)    (sprintf "~A ()" (cadr v)))
44      ((neg)       (sprintf "(- (~A))" (chicken-value (cadr v))))
45      ((+ - * / >= <= > <) 
46       (sprintf "(~A ~A ~A)"
47                (car v) 
48                (chicken-value (cadr v)) 
49                (chicken-value (caddr v))))
50      ((log ln sin cos cosh tanh exp)
51       (sprintf "(~A ~A)"
52                (car v) (chicken-value (cadr v)) ))
53      (else (error 'chicken-value "invalid value" v))))
54   ((boolean? v)  (if v "#t" "#f"))
55   (else (sprintf "~A" v))))
56
57(define chicken-run
58#<<EOF
59
60(define-syntax run
61  (syntax-rules ()
62    ((_ f indep dep events end input parameters)
63     (let ((nstate input))
64       (printf "# ~A " indep)
65       (for-each (lambda (x) (printf "~A " x)) dep)
66       (printf "~%")
67       (let recur ((nstate nstate))
68         (let ((ival (alist-ref indep nstate)))
69           (printf "~A " ival)
70           (for-each (lambda (x)
71                       (let ((v (alist-ref x nstate)))
72                         (printf "~A " (if (boolean? v) (or (and v 1) 0) v))))
73                     dep)
74           (printf "~%")
75           (if (> ival end)
76               (print "# All done!")
77               (recur (append (f nstate) parameters))))))
78     )))
79EOF
80)
81
82
83(define chicken-run/cvode
84#<<EOF
85
86
87(define-syntax run
88  (syntax-rules ()
89    ((_ f h indep dep events end input parameters)
90     (let* ((nstate (make-parameter input))
91
92            (numstates (filter (lambda (y) (number? (alist-ref y input))) dep))
93           
94            (update-nstate (lambda (old new)
95                             (old (fold (lambda (kv ax) (alist-update (car kv) (cdr kv) ax)) (old) new))))
96
97            (vector->nstate (lambda (yvec) 
98                              (fold (lambda (y yv nstate) (cons (cons y yv) nstate))
99                                    '()
100                                    numstates (f64vector->list yvec))))
101            (nstate->vector (lambda (nstate) 
102                              (list->f64vector (map (lambda (y) (alist-ref y nstate)) numstates))))
103
104            (fwrap (lambda (t yvec nstate)
105                     (let ((nstate1 (cons (cons indep t) (vector->nstate yvec))))
106                       (nstate nstate1)
107                       (let ((nstate2 (f (append nstate1 parameters))))
108                         (nstate->vector nstate2)
109                         ))
110                     ))
111
112            (fcall (lambda (t yvec nstate)
113                     (let ((nstate1 (cons (cons indep t) (vector->nstate yvec))))
114                       (nstate nstate1)
115                       (let ((nstate2 (f (append nstate1 parameters))))
116                         nstate2
117                         ))
118                     ))
119
120            (event-detect (lambda (t yy nstate)
121                            (let ((Vpeak (alist-ref 'Vpeak parameters))
122                                  (V (f64vector-ref yy 2)))
123                              (s32vector (floor (- V Vpeak))))))
124
125
126            (solver (cvode-create-solver
127                     (alist-ref indep input) 
128                     (nstate->vector input)
129                     fwrap
130                     abstol: 1e-4
131                     reltol: 1e-4
132                     events: (make-s32vector (length events) 0)
133                     event-fn: event-detect
134                     user-data: nstate
135                     ))
136
137            (solvewrap (lambda (solver tnext)
138                         (let ((status (cvode-solve solver tnext)))
139                           (if (negative? status) (error "CVODE solver error" status))
140                           (if (zero? status)
141                               (cons (cons indep tnext) (vector->nstate (cvode-yy solver)))
142                               (let ((yy (cvode-yy solver)) (tt (cvode-t solver)))
143                                 (nstate (f (append (fcall tt yy nstate) parameters)))
144                                 (cvode-reinit-solver solver tnext (nstate->vector (nstate)))
145                                 (cons (cons indep tt) (vector->nstate (cvode-yy solver)))
146                                 )))
147                           ))
148            )
149       (printf "# ~A " indep)
150       (for-each (lambda (x) (printf "~A " x)) dep)
151       (printf "~%")
152       (let recur ((hv (alist-ref h parameters)))
153         (let ((ival (alist-ref indep (nstate))))
154           (printf "~A " ival)
155           (for-each (lambda (x)
156                       (let ((v (alist-ref x (nstate))))
157                         (printf "~A " (if (boolean? v) (or (and v 1) 0) v))))
158                     dep)
159           (printf "~%")
160           (if (> ival end)
161               (print "# All done!")
162               (let ((nstate1 (solvewrap solver (+ hv ival))))
163                 (recur hv)))
164           ))
165       ))
166    ))
167
168EOF
169)
170
171(define (ivp-chicken prefix ivp-id ivar dvars pvars events start end ic sd
172                     #!key (solver 'rk3))
173
174  (let* ((dir (or (pathname-directory prefix) "."))
175         (solver-path (make-pathname dir (conc ivp-id "_solver.scm")))
176         (run-path    (make-pathname dir (sprintf "~A_run.scm" ivp-id)))
177         (exec-path   (make-pathname dir (sprintf "~A_run" ivp-id)))
178         (log-path    (make-pathname dir (sprintf "~A_~A.log" (pathname-file prefix) ivp-id)))
179         (csc-path    (make-pathname (program-path) "csc")))
180   
181    (make 
182        (
183         (solver-path (prefix)
184                      (with-output-to-file solver-path
185                        (lambda () (codegen/scheme ivp-id sd solver: solver))))
186         
187         (run-path (prefix)
188                   (with-output-to-file run-path
189                     (lambda () 
190                       (print-fragments
191                        (list
192                         (sprintf "(include \"~A_solver.scm\")~%~%" ivp-id)
193                         chicken-run nl
194                         (sprintf "(define initial (quote ~A))~%~%" 
195                                  (cons (cons ivar start) 
196                                        (map (lambda (x)
197                                               (let ((n (car x))
198                                                     (v (chicken-value (cdr x))))
199                                                 (cons n v)))
200                                             ic)))
201                         (sprintf "(define parameters (quote ~A))~%~%" 
202                                  (map (lambda (x) (assoc x ic)) pvars))
203                         (sprintf "(run ~A (quote ~A) (quote ~A) (quote ~A) ~A initial parameters)~%~%" 
204                                  ivp-id ivar dvars events end)
205                         ))
206                       ))
207                   )
208         
209         (exec-path (run-path solver-path)
210                    (run (,csc-path -w -I ,dir -b -S -d0 -O3 -disable-interrupts ,run-path)))
211         
212         (log-path (exec-path) (run (,exec-path > ,log-path)))
213         )
214     
215      (list log-path) )
216    ))
217
218
219(define (ivp-chicken/cvode prefix ivp-id hvar ivar dvars pvars events start end ic sd)
220
221  (let* ((dir (or (pathname-directory prefix) "."))
222         (solver-path (make-pathname dir (conc ivp-id "_solver.scm")))
223         (run-path    (make-pathname dir (sprintf "~A_run.scm" ivp-id)))
224         (exec-path   (make-pathname dir (sprintf "~A_run" ivp-id)))
225         (log-path    (make-pathname dir (sprintf "~A_~A.log" (pathname-file prefix) ivp-id)))
226         (csc-path    (make-pathname (program-path) "csc")))
227   
228    (make 
229        (
230         (solver-path (prefix)
231                      (with-output-to-file solver-path
232                        (lambda () (codegen/scheme ivp-id sd solver: 'cvode))))
233         
234         (run-path (prefix)
235                   (with-output-to-file run-path
236                     (lambda () 
237                       (print-fragments
238                        (list
239                         (sprintf "(include \"~A_solver.scm\")~%~%" ivp-id)
240                         chicken-run/cvode nl
241                         (sprintf "(define initial (quote ~A))~%~%" 
242                                  (cons (cons ivar start) 
243                                        (map (lambda (x)
244                                               (let ((n (car x))
245                                                     (v (chicken-value (cdr x))))
246                                                 (cons n v)))
247                                             ic)))
248                         (sprintf "(define parameters (quote ~A))~%~%" 
249                                  (map (lambda (x) (assoc x ic)) pvars))
250                         (sprintf "(run ~A (quote ~A) (quote ~A) (quote ~A) (quote ~A) ~A initial parameters)~%~%" 
251                                  ivp-id hvar ivar dvars events end)
252                         )))))
253         
254         (exec-path (run-path solver-path)
255                    (run (,csc-path -w -I ,dir -b -S -d0 -O3 -disable-interrupts ,run-path)))
256         
257         (log-path (exec-path) (run (,exec-path > ,log-path)))
258         )
259     
260      (list log-path) )
261    ))
262
263
264(define (ivp-chicken-codegen prefix ivp-id ivar dvars pvars events ic sd
265                     #!key (solver 'rk3))
266
267  (let* ((dir (or (pathname-directory prefix) "."))
268         (solver-path (make-pathname dir (conc ivp-id "_solver.scm")))
269         (run-path    (make-pathname dir (sprintf "~A_run.scm" ivp-id)))
270         )
271   
272    (make 
273        (
274         (solver-path (prefix)
275                      (with-output-to-file solver-path
276                        (lambda () (codegen/scheme ivp-id sd solver: solver))))
277         
278         (run-path (prefix)
279                   (with-output-to-file run-path
280                     (lambda () 
281                       (print-fragments
282                        (list
283                         (sprintf "(include \"~A_solver.scm\")~%~%" ivp-id)
284                         chicken-run nl
285                         (sprintf "(define initial (quote ~A))~%~%" 
286                                  (cons (cons ivar 0.0) 
287                                        (map (lambda (x)
288                                               (let ((n (car x))
289                                                     (v (chicken-value (cdr x))))
290                                                 (cons n v)))
291                                             ic)))
292                         (sprintf "(define parameters (quote ~A))~%~%" 
293                                  (map (lambda (x) (assoc x ic)) pvars))
294                         ))
295                       ))
296                   )
297         
298         )
299     
300      (list solver-path run-path) )
301    ))
302
303
304(define (ivp-chicken-codegen/cvode prefix ivp-id hvar ivar dvars pvars events ic sd)
305
306  (let* ((dir (or (pathname-directory prefix) "."))
307         (solver-path (make-pathname dir (conc ivp-id "_solver.scm")))
308         (run-path    (make-pathname dir (sprintf "~A_run.scm" ivp-id)))
309         )
310
311    (make 
312        (
313         (solver-path (prefix)
314                      (with-output-to-file solver-path
315                        (lambda () (codegen/scheme ivp-id sd solver: 'cvode))))
316         
317         (run-path (prefix)
318                   (with-output-to-file run-path
319                     (lambda () 
320                       (print-fragments
321                        (list
322                         (sprintf "(include \"~A_solver.scm\")~%~%" ivp-id)
323                         chicken-run/cvode nl
324                         (sprintf "(define initial (quote ~A))~%~%" 
325                                  (cons (cons ivar 0.0) 
326                                        (map (lambda (x)
327                                               (let ((n (car x))
328                                                     (v (chicken-value (cdr x))))
329                                                 (cons n v)))
330                                             ic)))
331                         (sprintf "(define parameters (quote ~A))~%~%" 
332                                  (map (lambda (x) (assoc x ic)) pvars))
333                         )))))
334         )
335     
336      (list solver-path run-path) )
337    ))
338
339
340)
Note: See TracBrowser for help on using the repository browser.