source: project/ode/trunk/extensions/ode-stalin.scm @ 8365

Last change on this file since 8365 was 8365, checked in by Ivan Raikov, 13 years ago

Added stalin-related extension.

File size: 12.4 KB
Line 
1;;
2;;
3;; An extension for translating ODE systems to Stalin Scheme programs.
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 ode)
22(require-extension srfi-1)
23(require-extension srfi-4)
24(require-extension srfi-13)
25(require-extension srfi-14)
26(require-extension environments)
27(require-extension runcmd)
28(require-extension utils)
29(require-extension lolevel)
30
31
32(define (lookup-def k lst . rest)
33  (let-optionals rest ((default #f))
34      (let ((kv (assoc k lst)))
35        (if (not kv) default
36            (match kv ((k v) v) (else (cdr kv)))))))
37
38(define nl "\n")
39
40(define (s+ x . lst)
41  (apply string-append (map ->string (cons x lst))))
42
43(define ($ . rest)
44  (string->symbol (apply s+ rest)))
45
46(define (fnoutf explist)
47  (for-each (lambda (str) (printf "~a\n" str)) explist))
48
49(define-macro (outf . explist)
50  `(fnoutf (list ,@(map (lambda (x) (list 'quasiquote x)) explist))))
51
52
53(define (enumprocs expr ax)
54  (match expr 
55         (('if . es)  (fold enumprocs ax es))
56         ((s . es)    (if (symbol? s)  (cons s (fold enumprocs ax es)) ax))
57         (else ax)))
58
59(define (enumvars expr ax)
60  (match expr 
61         (('if . es)  (fold enumvars ax es))
62         ((s . es)    (if (symbol? s)  (fold enumvars ax es) ax))
63         (id          (if (symbol? id) (cons id ax) ax))))
64
65(define (rhsvars rhs)
66  (if (rhs-procedure? rhs)
67      (rhs-procedure-vars rhs)
68      (enumvars rhs (list))))
69
70(define (rhsexpr dyenv expr)
71  (match expr 
72         (('if . es)  `(if . ,(map (lambda (x) (rhsexpr dyenv x)) es)))
73         (('pow x y)  (if (integer? y) `(let ((x ,x)) (* . ,(list-tabulate y (lambda (i) 'x))))
74                          `(expt ,x ,y)))
75         ((s . es)    (if (symbol? s)  (cons s (map (lambda (x) (rhsexpr dyenv x)) es))))
76         (id          (if (symbol? id) 
77                          (let ((p (lookup-def id dyenv)))
78                            (if p `(e^ %e ,(car p)) id))
79                          id))))
80
81(define (rhsbody dyenv rhs)
82  (if (rhs-procedure? rhs) (rhs-procedure-proc rhs)
83      (rhsexpr dyenv rhs)))
84 
85
86(define-record kinds const fn asgn vasgn state)
87
88(define (define-const k v)
89  (outf (define ,k ,v)))
90
91(define (define-fn dyenv n proc)
92  (let ((lst (procedure-data proc))) 
93    (let ((rt       (lookup-def 'rt lst))
94          (formals  (lookup-def 'formals lst))
95          (vars     (lookup-def 'vars lst))
96          (body     (lookup-def 'body lst))
97          (consts   (lookup-def 'consts lst)))
98      (if (and vars body) (outf (define (,n . ,vars) ,body))))))
99
100(define (define-state dyenv n init d rhs)
101  (match-let (((i t) (lookup-def n dyenv)))
102             (outf (vector-set! solve-env ,i (vector ,init 0.0 0.0 0.0))
103                   (vector-set! eval-env  ,i ,init))))
104
105(define (define-asgn dyenv n v rhs)
106  (match-let (((i t) (lookup-def n dyenv)))
107             (outf  (vector-set! solve-env ,i (vector ,v 0.0 0.0))
108                    (vector-set! eval-env  ,i ,v))))
109
110(define (define-vasgn dyenv n v rhs)
111  (match-let (((i t) (lookup-def n dyenv)))
112             (outf  (vector-set! solve-env ,i (list ,v))
113                    (vector-set! eval-env  ,i ,v))))
114
115
116(define (define-state-rhs dyenv n i d rhs)
117  (let ((fname  ($ n "-rhs"))
118        (vars   (rhsvars rhs))
119        (fbody  (rhsbody dyenv rhs)))
120    (if (not (procedure? fbody))
121        (outf (define (,fname %e) ,fbody)))))
122
123(define (define-asgn-rhs dyenv n v rhs)
124  (let ((fname  ($ n "-rhs"))
125        (vars   (rhsvars rhs))
126        (fbody  (rhsbody dyenv rhs)))
127    (if (not (procedure? fbody))
128        (outf (define (,fname %e) ,fbody)))))
129 
130(define (define-vasgn-rhs dyenv n v rhs)
131  (let ((fname  ($ n "-rhs"))
132        (vars   (rhsvars rhs))
133        (fbody  (rhsbody dyenv rhs)))
134    (if (not (procedure? fbody))
135        (outf (define (,fname %e) ,fbody)))))
136
137
138(define (ode:stalin-translator sys solver start stop initial-h . rest)
139 (let-optionals rest ((fptype "double"))
140  (match-let
141    ((($ ode:quantity 'SYSNAME  sysname)   
142      (environment-ref sys (ode-intern 'name)))
143     (($ ode:quantity 'ODECORE  ode)       
144      (environment-ref sys (ode-intern 'odecore)))
145     (($ ode:quantity 'TIMESTEP dt _)       
146      (environment-ref sys (ode-intern 'timestep)))
147     (($ ode:quantity 'INDEP indep _)       
148      (environment-ref sys (ode-intern 'indep)))
149     (($ ode:quantity 'PRTLIST prtitems prtevery prtfrom) 
150      (environment-ref sys (ode-intern 'prtlist)))
151     (($ ode:quantity 'GUARDS guards) 
152      (environment-ref sys (ode-intern 'guards))))
153     
154                         
155    (let* ((sfname  (string-append (symbol->string sysname) "_stalin.scm"))
156           (deps*   ((ode 'depgraph*) sys))
157           (sks     (make-kinds (list) (list) (list) (list) (list)))
158           (dyenv   (cadr
159                     (fold
160                      (lambda (x ax)
161                        (match (cons (environment-ref sys x) ax)
162                               (((? procedure? proc) i dyenv)
163                                (kinds-fn-set! sks (cons (list x proc) (kinds-fn sks)))
164                                ax)
165                               ((($ ode:quantity 'CONST  n v) i dyenv)
166                                (kinds-const-set! sks (cons (list n v) (kinds-const sks)))
167                                ax)
168                               ((($ ode:quantity 'ASGN   n v rhs _) i dyenv)
169                                (kinds-asgn-set! sks (cons (list n v rhs) (kinds-asgn sks)))
170                                (list (+ 1 i) (alist-update! n (list i 'a) dyenv)))
171                               ((($ ode:quantity 'VASGN  n v rhs _) i dyenv)
172                                (kinds-vasgn-set! sks (cons (list n v rhs) (kinds-vasgn sks)))
173                                (list (+ 1 i) (alist-update! n (list i 'va) dyenv)))
174                               ((($ ode:quantity 'STATE  n init v d rhs _) i dyenv)
175                                (kinds-state-set! sks (cons (list n init d rhs) (kinds-state sks)))
176                                (list (+ 1 i) (alist-update! n (list i d) dyenv)))
177                               (else ax)))
178                      (list 2 (list `(,indep 0 'q) `(,dt 1 'a))) 
179                      (environment-symbols sys)))))
180      (print "dyenv = " dyenv)
181      (match-let (((_ _ g) deps*))
182         (let ((poset     (vector->list ((ode 'depgraph->bfs-dist-poset) g)))
183               (N_dyenv   (length dyenv)))
184           (print "poset = " poset)
185           (with-output-to-file sfname
186             (lambda ()
187               (outf 
188                ,nl
189                ;; TODO: definitions of multi-argument print and error
190                (define error panic) 
191                (define pow expt)
192                (define (neg x) (- x))
193                ,nl
194                (include "\"stalin-box.scm\"")
195                (include "\"abm4.scm\"")
196                (include "\"rkf45.scm\"")
197                (include "\"euler.scm\"")
198                (include "\"stalin-solvers.scm\"")
199                ,nl)
200               (outf
201                ";; Definitions of basic datastructures and global variables"
202                ;; (define-record qsym s i t)
203                (define qsym    (primitive-procedure make-structure qsym 3))
204                (define qsym-s  (primitive-procedure structure-ref qsym 0))
205                (define qsym-i  (primitive-procedure structure-ref qsym 1))
206                (define qsym-t  (primitive-procedure structure-ref qsym 2))
207                (define e^ vector-ref)
208                (define N ,N_dyenv)
209                (define solve-env  (make-vector N))
210                (define eval-env   (make-vector N))
211                (define indep      (qsym (quote ,indep) 0 #t))
212                (define timestep   (qsym (quote ,dt) 1 #f))
213                )
214
215               (outf ";;" 
216                     ";; System constants" 
217                     ";;")
218               (for-each (lambda (c) (apply define-const c)) (kinds-const sks))
219               (outf ";;" 
220                     ";; User-defined functions" 
221                     ";;")
222               (for-each (lambda (st) (apply define-fn (cons dyenv st))) (kinds-fn sks))
223               (outf ";;" 
224                     ";; System states" 
225                     ";;")
226               (outf (vector-set! solve-env ,(car (lookup-def indep dyenv)) (vector ,start 0.0 0.0 0.0)))
227               (for-each (lambda (st) (apply define-state (cons dyenv st))) (kinds-state sks))
228               (outf ";;" 
229                     ";; Assigned quantities" 
230                     ";;")
231               (for-each (lambda (st) (apply define-asgn (cons dyenv st))) (kinds-asgn sks))
232               (for-each (lambda (st) (apply define-vasgn (cons dyenv st))) (kinds-vasgn sks))
233               (outf ";;" 
234                     ";; Equation definitions" 
235                     ";;")
236               (for-each (lambda (st) (apply define-state-rhs (cons dyenv st))) (kinds-state sks))
237               (for-each (lambda (st) (apply define-asgn-rhs  (cons dyenv st))) (kinds-asgn sks))
238               (for-each (lambda (st) (apply define-vasgn-rhs (cons dyenv st))) (kinds-vasgn sks))
239               (outf ";;"
240                     ";; Sets that contain the names of states and assigned quantities"
241                     ";;"
242                     (define states  (list . ,(map (lambda (x) 
243                                                     (let ((sym (car x)))
244                                                       (match-let (((i t) (lookup-def sym dyenv)))
245                                                         `(qsym (quote ,sym) ,i (quote ,t)))))
246                                                   (kinds-state sks))))
247                     (define asgns   (list . ,(map (lambda (x) 
248                                                     (let ((sym (car x)))
249                                                       (match-let (((i t) (lookup-def sym dyenv)))
250                                                          `(qsym (quote ,sym) ,i (quote ,t)))))
251                                                   (append (kinds-asgn sks) (kinds-vasgn sks))))))
252               (outf  ";;"
253                     ";; Evaluation order"
254                     ";;"
255                     (define eval-poset 
256                       (list
257                        . ,(let loop ((poset poset) (new-poset (list)) (i 0))
258                             (if (null? poset) (reverse new-poset)
259                                 (let ((lst1
260                                        (let inner-loop ((lst (car poset)) (new-lst (list)))
261                                          (if (null? lst) new-lst
262                                              (let* ((x (car lst)) (q (cdr x)) 
263                                                     (qrhs ($ q "-rhs")))
264                                                (match-let (((i t) (lookup-def q dyenv)))
265                                                  (let ((qq `(list (quote ,t) (qsym (quote ,q) ,i (quote ,t)) ,qrhs)))
266                                                    (inner-loop (cdr lst) (cons qq new-lst)))))))))
267                                   (loop (cdr poset) (cons `(list ,i . ,lst1) new-poset) (+ 1 i))))))))
268               (outf ";;"
269                     ";; Guards"
270                     ";;"
271                     (define guards (list . ,(map (lambda (expr) `(lambda (eval-env) ,(rhsexpr dyenv expr))) guards)))
272                    )
273               (outf  (define order                     ,(length (kinds-state sks)))
274                      (define svec-value-idx            0)
275                      (define svec-abserr-idx           1)
276                      (define svec-relerr-idx           2)
277                      (define svec-deriv-idx            3)
278                      (define (is-state? q)             (or (eq? (qsym-t q) 'q) (eq? (qsym-t q) 'd)))
279                      (define (is-asgn? q)              (not (or (eq? (qsym-t q) 'q) (eq? (qsym-t q) 'd))))
280                      (define (eval-poset-foreach proc) 
281                        (for-each (lambda (x) 
282                                    (let ((order (car x))
283                                          (lst   (cdr x)))
284                                      (proc order lst)))
285                                  eval-poset))
286                      (define (eval-poset-fold proc initial)
287                        (fold (lambda (x ax)
288                                (let ((order (car x))
289                                      (lst   (cdr x)))
290                                  (proc order ax lst))) 
291                              initial eval-poset))
292                      (define (eval-env-ref q)          (vector-ref eval-env (qsym-i q)))
293                      (define (eval-env-set! q v)       (vector-set! eval-env (qsym-i q) v))
294                      (define (solve-env-ref q)         (vector-ref solve-env (qsym-i q)))
295                      (define (solve-env-set! q i v)    (vector-set! (vector-ref solve-env (qsym-i q)) i v))
296                      (define (solve-env-state-foreach proc)
297                        (for-each (lambda (x) (apply proc x)) 
298                                  (map (lambda (q) (list q (vector-ref solve-env (qsym-i q)))) states)))
299                      (define (solve-env-state-fold proc init)
300                        (fold (lambda (x ax) (apply proc (append x (list ax))))
301                              init  (map (lambda (q) (list q (vector-ref solve-env (qsym-i q)))) states)))
302                      (define (solve-env-asgn-fold proc init)
303                        (fold (lambda (x ax) (apply proc (append x (list ax))))
304                              init  (map (lambda (q) (list q (vector-ref solve-env (qsym-i q)))) asgns)))
305                      (define (eval-rhs q qrhs . rest)  (qrhs eval-env))
306                      (define (eval-guard gfun)         (gfun eval-env))
307                      (define hmax-factor               2)
308                      (define hmin-factor               0.5)
309                      (define hscale-factor             0.9)
310                      (define  hmin                     1e-16)
311                      (define  hmax                     1.0)
312                      (define relmax                    ,(ode 'relmax))
313                      (define relmin                    ,(ode 'relmin))
314                      (define absmax                    ,(ode 'absmax))
315                      (define absmin                    ,(ode 'absmin))
316                      (define prti        0)
317                      (define prtevery    ,prtevery)
318                      (define prtfrom     ,prtfrom)
319                      (define prtitems    (list
320                                           . ,(map (lambda (x)
321                                                     (match-let ((( tag sym ) x))
322                                                        (match-let (((i t) (lookup-def sym dyenv)))
323                                                          `(list (quote ,tag )
324                                                                 (qsym (quote ,sym) ,i (quote ,t))))))
325                                                   prtitems)))
326
327                      (include "\"stalin-numerr.scm\"")
328                      (include "\"stalin-qprint.scm\"")
329                      (include "\"stalin-dispatch.scm\"") 
330                      )
331               (outf ";;"
332                     ";; Main simulation loop"
333                     ";;"
334                     ((,solver error make-vector vector-ref vector-set! vector?)
335                      ode-runtime indep ,start ,stop timestep ,initial-h guards))
336               ))))))))
337
338
339;; stalin -copt "-O3 -finline-functions-called-once"
340;; -no-clone-size-limit -d -k -On -I ../extensions -I ..
341;; file.scm
Note: See TracBrowser for help on using the repository browser.