1 | |
---|
2 | ;; |
---|
3 | ;; Definition and code generators for a simple applicative language |
---|
4 | ;; for numerical simulation. |
---|
5 | ;; |
---|
6 | ;; Copyright 2010-2013 Ivan Raikov and the Okinawa Institute of |
---|
7 | ;; Science and Technology. |
---|
8 | ;; |
---|
9 | ;; This program is free software: you can redistribute it and/or |
---|
10 | ;; modify it under the terms of the GNU General Public License as |
---|
11 | ;; published by the Free Software Foundation, either version 3 of the |
---|
12 | ;; License, or (at your option) any later version. |
---|
13 | ;; |
---|
14 | ;; This program is distributed in the hope that it will be useful, but |
---|
15 | ;; WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
16 | ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
---|
17 | ;; General Public License for more details. |
---|
18 | ;; |
---|
19 | ;; A full copy of the GPL license can be found at |
---|
20 | ;; <http://www.gnu.org/licenses/>. |
---|
21 | ;; |
---|
22 | |
---|
23 | (module flsim |
---|
24 | |
---|
25 | ( |
---|
26 | value? |
---|
27 | V:C V:Var V:Rec V:Sel V:Vec V:Sub V:Fn V:Op V:Ifv |
---|
28 | |
---|
29 | expr? |
---|
30 | E:Ife E:Let E:Ret E:Noop |
---|
31 | |
---|
32 | binding? B:Val |
---|
33 | |
---|
34 | typenv |
---|
35 | |
---|
36 | name/scheme prelude/scheme expr->scheme value->scheme binding->scheme |
---|
37 | name/ML prelude/ML expr->ML value->ML binding->ML |
---|
38 | name/Octave prelude/Octave expr->Octave value->Octave binding->Octave |
---|
39 | ) |
---|
40 | |
---|
41 | (import scheme chicken) |
---|
42 | |
---|
43 | (require-extension extras data-structures srfi-1 datatype) |
---|
44 | |
---|
45 | (define nl "\n") |
---|
46 | |
---|
47 | |
---|
48 | (define-datatype value value? |
---|
49 | (V:C (v (lambda (x) (or (symbol? x) (number? x) )))) |
---|
50 | (V:Var (name symbol?)) |
---|
51 | (V:Rec (flds (lambda (xs) (and (pair? xs) (every (lambda (x) (and (symbol? (car x)) (value? (cadr x)))) xs))))) |
---|
52 | (V:Sel (field (lambda (x) (or (symbol? x) (integer? x)))) (v value?)) |
---|
53 | (V:Vec (vals (lambda (x) (every value? x)))) |
---|
54 | (V:Sub (index (lambda (x) (and (integer? x) (or (zero? x) (positive? x))))) (v value?)) |
---|
55 | (V:Fn (args list?) (body expr?)) |
---|
56 | (V:Op (name symbol?) (args (lambda (x) (every value? x)))) |
---|
57 | (V:Ifv (test value?) (ift value?) (iff value?)) |
---|
58 | ) |
---|
59 | |
---|
60 | |
---|
61 | (define-record-printer (value x out) |
---|
62 | (fprintf out "~A" |
---|
63 | (cases value x |
---|
64 | (V:C (v) (sprintf "V:C ~A" v)) |
---|
65 | (V:Var (n) (sprintf "V:Var ~A " n)) |
---|
66 | (V:Rec (lst) (sprintf "V:Rec ~A" lst)) |
---|
67 | (V:Sel (f v) (sprintf "V:Sel ~A ~A" f v)) |
---|
68 | (V:Vec (lst) (sprintf "V:Vec ~A" lst)) |
---|
69 | (V:Sub (i v) (sprintf "V:Sub ~A ~A" i v)) |
---|
70 | (V:Ifv (test tv fv) "V:Ifv ~A ~A ~A" test tv fv) |
---|
71 | (V:Fn (args body) (sprintf "V:Fn ~A = ~A" args body)) |
---|
72 | (V:Op (name args) (sprintf "(V:Op ~A ~A)" name args))))) |
---|
73 | |
---|
74 | |
---|
75 | (define-datatype binding binding? |
---|
76 | (B:Val (name symbol?) (v value?)) |
---|
77 | ) |
---|
78 | |
---|
79 | |
---|
80 | (define-record-printer (binding x out) |
---|
81 | (fprintf out "~A" |
---|
82 | (cases binding x |
---|
83 | (B:Val (name v) (sprintf "B:Val ~A = ~A" name v)) |
---|
84 | ))) |
---|
85 | |
---|
86 | |
---|
87 | (define-datatype expr expr? |
---|
88 | (E:Ife (test value?) (ift expr?) (iff expr?)) |
---|
89 | (E:Let (bnds (lambda (x) (every binding? x))) (body expr?)) |
---|
90 | (E:Ret (v value?)) |
---|
91 | (E:Noop) |
---|
92 | ) |
---|
93 | |
---|
94 | |
---|
95 | (define-record-printer (expr x out) |
---|
96 | (fprintf out "~A" |
---|
97 | (cases expr x |
---|
98 | (E:Ife (test ift iff) (sprintf "E:Ife ~A ~A ~A" test ift iff)) |
---|
99 | (E:Let (bnds body) (sprintf "E:Let ( ~A ) ~A" bnds body)) |
---|
100 | (E:Ret (v) (sprintf "E:Ret ~A" v)) |
---|
101 | (E:Noop () (sprintf "E:Noop")) |
---|
102 | ))) |
---|
103 | |
---|
104 | |
---|
105 | ;; Builds a type environment for an expression list |
---|
106 | (define (typenv bindings) |
---|
107 | |
---|
108 | (define (numeric-type? x) |
---|
109 | (case x ((numeric) #t) (else #f))) |
---|
110 | (define (boolean-type? x) |
---|
111 | (case x ((bool) #t) (else #f))) |
---|
112 | (define (type-equal? ty1 ty2) |
---|
113 | (case (car ty1) |
---|
114 | ((numeric boolean unit) |
---|
115 | (eq? (car ty1) (car ty2))) |
---|
116 | ((record) |
---|
117 | (case (car ty2) |
---|
118 | ((record) |
---|
119 | (lset= (lambda (x y) (and (eq? (car x) (car y)) |
---|
120 | (type-equal? (cdr x) (cdr y)))) |
---|
121 | (cdr ty1) (cdr ty2))) |
---|
122 | (else #f))) |
---|
123 | ((vector) |
---|
124 | (case (car ty2) |
---|
125 | ((vector) |
---|
126 | (and (= (length (cdr ty1)) (length (cdr ty2))) |
---|
127 | (every type-equal? (cdr ty1) (cdr ty2)))) |
---|
128 | (else #f))) |
---|
129 | ((function) |
---|
130 | (case (car ty2) |
---|
131 | ((function) |
---|
132 | (and (equal? (cadr ty1) (cadr ty2)) |
---|
133 | (equal? (caddr ty1) (caddr ty2)))))) |
---|
134 | )) |
---|
135 | |
---|
136 | (define (value-type v env) |
---|
137 | (cases value v |
---|
138 | (V:C (v) (if (number? v) 'number |
---|
139 | (case v |
---|
140 | ((false true) 'bool) |
---|
141 | (else (error 'value-type "unknown constant" v))))) |
---|
142 | (V:Var (s) (alist-ref s env)) |
---|
143 | (V:Rec (flds) |
---|
144 | `(record . ,(map (lambda (f) (cons (car f) (value-type (cadr f) env))) flds))) |
---|
145 | (V:Sel (field u) |
---|
146 | (let ((ut (value-type u env))) |
---|
147 | (cond ((and (eq? (car ut) 'record) (symbol? field)) |
---|
148 | (alist-ref field (cdr ut))) |
---|
149 | (else (error 'value-type "invalid select value" field u)) |
---|
150 | ))) |
---|
151 | (V:Vec (vs) `(vector . ,(map (lambda (x) (value-type x env)) vs))) |
---|
152 | (V:Sub (index u) |
---|
153 | (let ((ut (value-type u env))) |
---|
154 | (cond |
---|
155 | ((and (eq? (car ut) 'vector) (integer? index)) |
---|
156 | (let ((limit (length (cdr ut)))) |
---|
157 | (if (not (or (zero? index) (positive? index))) |
---|
158 | (error 'value-type "invalid vector index" index u)) |
---|
159 | (if (< index limit) |
---|
160 | (error 'value-type "vector index out of range" index u)) |
---|
161 | (list-ref (cdr ut) index) |
---|
162 | )) |
---|
163 | (else (error 'value-type "invalid sub value" index u)) |
---|
164 | ))) |
---|
165 | (V:Fn (args body) `(function ,args ,body)) |
---|
166 | (V:Op (name args) |
---|
167 | (let ((f (alist-ref name env))) |
---|
168 | (if f |
---|
169 | (apply-type f args env) |
---|
170 | (case name |
---|
171 | ((+ - * /) |
---|
172 | (if (every numeric-type? (map (lambda (x) (value-type x env)) args)) |
---|
173 | 'number |
---|
174 | (error 'value-type "non-numeric arguments to arithmetic operator" |
---|
175 | name args))) |
---|
176 | ((>= > < <=) |
---|
177 | (if (every numeric-type? (map (lambda (x) (value-type x env)) args)) |
---|
178 | 'bool |
---|
179 | (error 'value-type "non-numeric arguments to comparison operator" |
---|
180 | name args))) |
---|
181 | |
---|
182 | ((TRSA TRSB) |
---|
183 | (let ((ty (value-type (car args) env))) |
---|
184 | `(record (,name . ,ty)))) |
---|
185 | |
---|
186 | (else (error 'value-type "unknown operator" name)) |
---|
187 | )) |
---|
188 | )) |
---|
189 | (V:Ifv (test ift iff) |
---|
190 | (if (boolean-type? (value-type test env)) |
---|
191 | (let ((ty1 (value-type ift env)) |
---|
192 | (ty2 (value-type iff env))) |
---|
193 | (if (type-equal? ty1 ty2) ty1 |
---|
194 | (error 'value-type "if value branches have different types" |
---|
195 | test ift ty1 iff ty2))) |
---|
196 | (error 'value-type "if value has non-boolean test value"))) |
---|
197 | )) |
---|
198 | |
---|
199 | (define (binding-env b env) |
---|
200 | (cases binding b |
---|
201 | (B:Val (name v) |
---|
202 | (let ((vt (value-type v env))) |
---|
203 | (cons `(name . vt) env))) |
---|
204 | )) |
---|
205 | |
---|
206 | (define (expr-type e env) |
---|
207 | (cases expr e |
---|
208 | (E:Ife (test ift iff) |
---|
209 | (if (boolean-type? (value-type test env)) |
---|
210 | (let ((ty1 (expr-type ift env)) |
---|
211 | (ty2 (expr-type iff env))) |
---|
212 | (if (type-equal? ty1 ty2) ty1 |
---|
213 | (error 'expr-type "if expression branches have different types" |
---|
214 | test ift ty1 iff ty2))) |
---|
215 | (error 'expr-type "if expression has non-boolean test value"))) |
---|
216 | (E:Let (bnds body) |
---|
217 | (let ((env1 (fold binding-env env bnds))) |
---|
218 | (expr-type e env1))) |
---|
219 | (E:Ret (v) (value-type v env)) |
---|
220 | (E:Noop () 'unit))) |
---|
221 | |
---|
222 | (define (apply-type f args env) |
---|
223 | (let ((formals (cadr f)) |
---|
224 | (body (caddr f))) |
---|
225 | (let ((env (map (lambda (formal arg) `(,formal . ,(value-type arg env))) formals args))) |
---|
226 | (expr-type body env)))) |
---|
227 | |
---|
228 | (fold binding-env '() bindings) |
---|
229 | |
---|
230 | ) |
---|
231 | |
---|
232 | |
---|
233 | |
---|
234 | (define (name/scheme s) |
---|
235 | (let ((cs (string->list (->string s)))) |
---|
236 | (let loop ((lst (list)) (cs cs)) |
---|
237 | (if (null? cs) (string->symbol (list->string (reverse lst))) |
---|
238 | (let* ((c (car cs)) |
---|
239 | (c1 (cond ((or (char-alphabetic? c) (char-numeric? c) |
---|
240 | (char=? c #\_) (char=? c #\-)) c) |
---|
241 | (else #\-)))) |
---|
242 | (loop (cons c1 lst) (cdr cs))))))) |
---|
243 | |
---|
244 | (define (binding->scheme x) |
---|
245 | (cases binding x |
---|
246 | (B:Val (name v) |
---|
247 | (list "(" (name/scheme name) " " (value->scheme v) ")" nl)))) |
---|
248 | |
---|
249 | |
---|
250 | (define (expr->scheme x . rest) |
---|
251 | (let-optionals rest ((bnd? #t)) |
---|
252 | (cases expr x |
---|
253 | |
---|
254 | (E:Ife (test ift iff) |
---|
255 | (list "(cond (" (value->scheme test) " " nl |
---|
256 | " " (expr->scheme ift ) ")" nl |
---|
257 | "(else " (expr->scheme iff) "))" nl)) |
---|
258 | |
---|
259 | (E:Let (bnds body) |
---|
260 | (list "(let* (" nl |
---|
261 | (map (lambda (x) (binding->scheme x)) bnds) nl |
---|
262 | ") " nl |
---|
263 | (expr->scheme body #f) nl |
---|
264 | ")" nl)) |
---|
265 | |
---|
266 | (E:Ret (v) (value->scheme v)) |
---|
267 | |
---|
268 | (E:Noop () (list "(void)")) |
---|
269 | ))) |
---|
270 | |
---|
271 | |
---|
272 | (define (value->scheme v) |
---|
273 | (let ((result |
---|
274 | (cases value v |
---|
275 | (V:C (v) v) |
---|
276 | (V:Var (name) (name/scheme name)) |
---|
277 | (V:Rec (lst) |
---|
278 | (list "`(" (intersperse (map (lambda (nv) (list "(" (name/scheme (car nv)) " . ," |
---|
279 | (value->scheme (cadr nv)) ")")) lst) " ") ")")) |
---|
280 | (V:Sel (field v) |
---|
281 | (if (number? field) |
---|
282 | (list "(list-ref " (value->scheme v) " " (- field 1) ")") |
---|
283 | (list "(alist-ref '" (name/scheme field) " " (value->scheme v) ")"))) |
---|
284 | (V:Vec (lst) |
---|
285 | (list "(list " (intersperse (map value->scheme lst) " ") ")")) |
---|
286 | (V:Sub (index v) |
---|
287 | (list "(list-ref " (value->scheme v) " " index ")")) |
---|
288 | (V:Fn (args body) |
---|
289 | (list "(lambda (" (intersperse (map name/scheme args) " ") ") " |
---|
290 | (expr->scheme body #f) ")")) |
---|
291 | (V:Op (name args) |
---|
292 | (let* ((fp? (case name |
---|
293 | ((+ - * / >= > < <= neg) #t) |
---|
294 | (else #f))) |
---|
295 | (op (if fp? (conc "fp" name) name))) |
---|
296 | (cond ((null? args) |
---|
297 | (case name |
---|
298 | ((NONE) (list "#f")) |
---|
299 | (else (list "(" name ")")))) |
---|
300 | |
---|
301 | (fp? |
---|
302 | (if (pair? (cdr args)) |
---|
303 | (fold-right (lambda (x ax) (list "(" op " " (value->scheme x) " " ax ")")) |
---|
304 | (list "(" op " " (value->scheme (car args)) " " (value->scheme (cadr args)) ")") |
---|
305 | (cddr args)) |
---|
306 | (list "(" op " " (value->scheme (car args)) ")") |
---|
307 | )) |
---|
308 | |
---|
309 | (else |
---|
310 | (list "(" op " " (intersperse (map value->scheme args) " ") ")"))))) |
---|
311 | (V:Ifv (test ift iff) |
---|
312 | (list "(if " (value->scheme test) " " |
---|
313 | (value->scheme ift) " " |
---|
314 | (value->scheme iff) ")")) |
---|
315 | |
---|
316 | ))) |
---|
317 | result)) |
---|
318 | |
---|
319 | (define (prelude/scheme #!key (solver 'rk4b) (random #f) (integral-index 0)) |
---|
320 | |
---|
321 | `( |
---|
322 | ,(case solver |
---|
323 | ((cvode) `("(use sundials random-mtzig mathh)" ,nl)) |
---|
324 | (else `("(use runge-kutta random-mtzig mathh)" ,nl))) |
---|
325 | #<<EOF |
---|
326 | |
---|
327 | ;; Variant types |
---|
328 | |
---|
329 | (define-syntax define-datatype |
---|
330 | (syntax-rules () |
---|
331 | [(_ type (name field ...) ...) |
---|
332 | (begin |
---|
333 | (define-constructors type ((name field ...) ...)))])) |
---|
334 | |
---|
335 | |
---|
336 | (define-syntax define-constructors |
---|
337 | (syntax-rules () |
---|
338 | [(define-constructors type ((name field ...) ...)) |
---|
339 | (define-constructors type ((name field ...) ...) (name ...))] |
---|
340 | [(define-constructors type ((name field ...) ...) names) |
---|
341 | (begin |
---|
342 | (define-constructor type (name field ...) names) |
---|
343 | ...)])) |
---|
344 | |
---|
345 | |
---|
346 | (define-syntax define-constructor |
---|
347 | (syntax-rules () |
---|
348 | [(_ type (name field ...) names) |
---|
349 | (define (name field ...) |
---|
350 | (cons 'type |
---|
351 | (lambda names |
---|
352 | (name field ...))))])) |
---|
353 | |
---|
354 | |
---|
355 | (define-syntax cases |
---|
356 | (syntax-rules () |
---|
357 | [(_ type x [(name field ...) exp] |
---|
358 | ...) |
---|
359 | ((cdr x) (lambda (field ...) exp) |
---|
360 | ...)])) |
---|
361 | |
---|
362 | (define-datatype trs (TRSA a) (TRSB b)) |
---|
363 | (define-datatype trc (TRC f fk e)) |
---|
364 | |
---|
365 | (define (tsCase fa fb x) (cases trs x ((TRSA a) (fa a)) ((TRSB b) (fb b)))) |
---|
366 | (define (trfOf x) (cases trc x ((TRC f fk e) f))) |
---|
367 | (define (trfkOf x) (cases trc x ((TRC f fk e) fk))) |
---|
368 | (define (treOf x) (cases trc x ((TRC f fk e) e))) |
---|
369 | |
---|
370 | (define-datatype option (NONE) (SOME a)) |
---|
371 | (define (swap x v) (cases option v ((NONE) x) ((SOME v) v))) |
---|
372 | |
---|
373 | (define false #f) |
---|
374 | (define true #t) |
---|
375 | |
---|
376 | (define equal equal?) |
---|
377 | |
---|
378 | (define (signalOf v) (if (not v) (error 'signalOf "empty signal" v) v)) |
---|
379 | |
---|
380 | (define (heaviside x) (if (negative? x) 0 1)) |
---|
381 | |
---|
382 | EOF |
---|
383 | |
---|
384 | ,(if (positive? integral-index) |
---|
385 | (case solver |
---|
386 | ((cvode) |
---|
387 | (list "(define integral-solvers (make-vector " integral-index " #f))" nl)) |
---|
388 | (else '()) |
---|
389 | ) '()) |
---|
390 | |
---|
391 | ,(if (not solver) |
---|
392 | `((";; dummy solver; returns only the computed derivatives") |
---|
393 | ("(define (integrate f x y h i) (f x y))" ,nl) |
---|
394 | ) |
---|
395 | (case solver |
---|
396 | ((cvode) |
---|
397 | `( |
---|
398 | ("(define (integrate f x y h i) (f x y))" ,nl) |
---|
399 | )) |
---|
400 | |
---|
401 | (else |
---|
402 | `(("(define (scaler x a) (map (lambda (k) (fp* x k)) a))" ,nl) |
---|
403 | ("(define (summer a b) (map fp+ a b))" ,nl) |
---|
404 | ("(define " ,solver " (make-" ,solver "))" ,nl) |
---|
405 | ("(define (make_stepper deriv) (" ,solver " scaler summer deriv))" ,nl) |
---|
406 | ("(define (integrate f x y h i) (((make_stepper f) h) x y))" ,nl) |
---|
407 | )) |
---|
408 | )) |
---|
409 | )) |
---|
410 | |
---|
411 | (define (name/ML s) |
---|
412 | (let ((cs (string->list (->string s)))) |
---|
413 | (let loop ((lst (list)) (cs cs)) |
---|
414 | (cond ((null? cs) (string->symbol (list->string (reverse lst)))) |
---|
415 | ((null? (cdr cs)) |
---|
416 | (let ((c (car cs))) |
---|
417 | (if (or (char-alphabetic? c) (char-numeric? c)) |
---|
418 | (loop (cons c lst) (cdr cs)) |
---|
419 | (loop (append (reverse (string->list (->string (gensym 't)))) lst) (cdr cs)) |
---|
420 | ))) |
---|
421 | (else |
---|
422 | (let* ((c (car cs)) |
---|
423 | (c1 (cond ((or (char-alphabetic? c) (char-numeric? c) |
---|
424 | (char=? c #\_) (char=? c #\#)) c) |
---|
425 | (else #\_)))) |
---|
426 | (loop (cons c1 lst) (cdr cs)))))))) |
---|
427 | |
---|
428 | |
---|
429 | (define (binding->ML x . rest) |
---|
430 | (let-optionals rest ((bnd? #t)) |
---|
431 | (cases binding x |
---|
432 | (B:Val (name v) |
---|
433 | (list "val " (name/ML name) " = " (value->ML v) nl)) |
---|
434 | ))) |
---|
435 | |
---|
436 | (define (expr->ML x . rest) |
---|
437 | (let-optionals rest ((bnd? #t)) |
---|
438 | (cases expr x |
---|
439 | |
---|
440 | (E:Ife (test ift iff) |
---|
441 | (list "if (" (value->ML test) ") " nl |
---|
442 | "then " (expr->ML ift ) nl |
---|
443 | "else " (expr->ML iff) nl)) |
---|
444 | |
---|
445 | (E:Let (bnds body) |
---|
446 | (list "let " nl |
---|
447 | (map (lambda (x) (binding->ML x #t)) bnds) nl |
---|
448 | "in " nl |
---|
449 | (expr->ML body #f) nl |
---|
450 | "end" nl)) |
---|
451 | |
---|
452 | (E:Ret (v) (value->ML v)) |
---|
453 | |
---|
454 | (E:Noop () (list "()")) |
---|
455 | ))) |
---|
456 | |
---|
457 | |
---|
458 | (define (value->ML v) |
---|
459 | (cases value v |
---|
460 | (V:C (v) (cond ((and (number? v) (negative? v)) |
---|
461 | (list "~" (abs v))) |
---|
462 | (else v))) |
---|
463 | (V:Var (name) (name/ML name)) |
---|
464 | (V:Rec (lst) |
---|
465 | (list "{" (intersperse (map (lambda (nv) (list (name/ML (first nv)) " = " |
---|
466 | (value->ML (cadr nv)))) lst) ", ") "}")) |
---|
467 | (V:Sel (field v) |
---|
468 | (list "(#" (name/ML field) "(" (value->ML v) "))")) |
---|
469 | (V:Vec (lst) |
---|
470 | (let ((n (length lst))) |
---|
471 | (list "(Vector.fromList [" (intersperse (map (lambda (v) (value->ML v)) lst) ", ") "])"))) |
---|
472 | (V:Sub (index v) |
---|
473 | (list "Unsafe.Vector.sub (" (value->ML v) ", " index ")")) |
---|
474 | (V:Fn (args body) |
---|
475 | (list "(fn (" (intersperse (map name/ML args) ",") ") => " |
---|
476 | (expr->ML body #f) ")")) |
---|
477 | (V:Op (name args) |
---|
478 | (let* ((infix? (case name |
---|
479 | ((+ - * / >= > < <=) #t) |
---|
480 | (else #f))) |
---|
481 | (op (if infix? (list "(op " name ")") name))) |
---|
482 | (let ((name1 (name/ML name))) |
---|
483 | (cond ((null? args) |
---|
484 | (case name |
---|
485 | ((NONE) (list name1)) |
---|
486 | (else (list name1 "()")))) |
---|
487 | ((null? (cdr args)) |
---|
488 | (list "(" op " " (value->ML (car args)) ")")) |
---|
489 | ((and infix? (pair? (cddr args))) |
---|
490 | (list "(foldr " op "(" (value->ML (V:Op name (list (car args) (cadr args)))) ")" |
---|
491 | "[" (intersperse (map value->ML (cddr args)) ",") "])")) |
---|
492 | (else |
---|
493 | (list "(" op "(" (intersperse (map value->ML args) ",") "))")))))) |
---|
494 | (V:Ifv (test ift iff) |
---|
495 | (list "(if (" (value->ML test) ") " |
---|
496 | "then " (value->ML ift ) " " |
---|
497 | "else " (value->ML iff) ")")) |
---|
498 | |
---|
499 | )) |
---|
500 | |
---|
501 | |
---|
502 | (define (prelude/ML #!key (solver 'rk4b) (random #f)) |
---|
503 | `( |
---|
504 | #<<EOF |
---|
505 | structure Model = |
---|
506 | struct |
---|
507 | |
---|
508 | open Real |
---|
509 | open Math |
---|
510 | open RungeKutta |
---|
511 | |
---|
512 | datatype ('b,'c) trs = TRSA of 'b | TRSB of 'c |
---|
513 | datatype ('a,'b,'c) trc = TRC of ((('a -> (('b,'c) trs))) * |
---|
514 | (('a -> (('b,'c) trs))) * |
---|
515 | ((('b,'c) trs) -> bool)) |
---|
516 | |
---|
517 | fun tsCase (fa,fb,x) = case x of TRSA a => (fa a) | TRSB b => (fb b) |
---|
518 | fun trfOf x = case x of TRC (f,fk,e) => f |
---|
519 | fun trfkOf x = case x of TRC (f,fk,e) => fk |
---|
520 | fun treOf x = case x of TRC (f,fk,e) => e |
---|
521 | |
---|
522 | fun putStrLn str = |
---|
523 | (TextIO.output (TextIO.stdOut, str); |
---|
524 | TextIO.output (TextIO.stdOut, "\n")) |
---|
525 | |
---|
526 | fun putStr str = (TextIO.output (TextIO.stdOut, str)) |
---|
527 | |
---|
528 | fun showReal n = |
---|
529 | let open StringCvt |
---|
530 | in |
---|
531 | (if n < 0.0 then "-" else "") ^ (fmt (FIX (SOME 12)) (abs n)) |
---|
532 | end |
---|
533 | |
---|
534 | fun vmap2 f (v1,v2) = |
---|
535 | let |
---|
536 | val n = Vector.length v1 |
---|
537 | in |
---|
538 | Vector.tabulate (n, fn (i) => f (Unsafe.Vector.sub (v1,i), |
---|
539 | Unsafe.Vector.sub (v2,i))) |
---|
540 | end |
---|
541 | |
---|
542 | exception EmptySignal |
---|
543 | |
---|
544 | val neg = (op ~) |
---|
545 | val swap = fn (x,v) => (case v of NONE => x | SOME v => v) |
---|
546 | val equal = fn (x,y) => (x = y) |
---|
547 | val signalOf = fn (v) => (case v of NONE => raise EmptySignal | SOME v => v) |
---|
548 | val heaviside = fn (v) => (if Real.< (v, 0.0) then 0.0 else 1.0) |
---|
549 | EOF |
---|
550 | |
---|
551 | ,(if random |
---|
552 | #<<EOF |
---|
553 | |
---|
554 | fun RandomInit () = RandomMTZig.fromEntropy() |
---|
555 | |
---|
556 | val RandomState = RandomInit () |
---|
557 | |
---|
558 | fun random_uniform () = RandomMTZig.randUniform RandomState |
---|
559 | |
---|
560 | |
---|
561 | fun PoissonInit () = |
---|
562 | let |
---|
563 | val zt = RandomMTZig.ztnew() |
---|
564 | val st = RandomMTZig.fromEntropy() |
---|
565 | in |
---|
566 | {st=st,zt=zt} |
---|
567 | end |
---|
568 | |
---|
569 | fun PoissonStep (rate,t,h,st,zt) = |
---|
570 | let |
---|
571 | val rv = RandomMTZig.randPoisson (rate*0.001*h,st,zt) |
---|
572 | val spike' = Real.> (rv,0.0) |
---|
573 | val spikeCount' = if spike' then rv else 0.0 |
---|
574 | in |
---|
575 | {t=t+h,spike=spike',spikeCount=spikeCount',st=st,zt=zt} |
---|
576 | end |
---|
577 | |
---|
578 | |
---|
579 | EOF |
---|
580 | "") |
---|
581 | |
---|
582 | ,(if (not solver) |
---|
583 | `(("(* dummy solver; returns only the computed derivatives *)") |
---|
584 | ("fun integrate (f,x: real,y: real vector,h,i) = (f (x,y))" ,nl) |
---|
585 | ) |
---|
586 | `(("val summer = fn (a,b) => (vmap2 (fn (x,y) => x+y) (a,b))" ,nl) |
---|
587 | ("val scaler = fn(a,lst) => (Vector.map (fn (x) => a*x) lst)" ,nl) |
---|
588 | . ,(case solver |
---|
589 | ;; adaptive solvers |
---|
590 | ((rkhe rkbs rkf45 rkck rkdp rkf78 rkv65) |
---|
591 | `( |
---|
592 | ("val " ,solver ": (real vector) stepper2 = make_" ,solver "()" ,nl) |
---|
593 | ("fun make_stepper (deriv) = " ,solver " (scaler,summer,deriv)" ,nl) |
---|
594 | #<<EOF |
---|
595 | |
---|
596 | |
---|
597 | |
---|
598 | val tol = Real.Math.pow (10.0, ~6.0) |
---|
599 | val lb = 0.5 * tol |
---|
600 | val ub = 1.0 * tol |
---|
601 | val maxit = 100 |
---|
602 | |
---|
603 | datatype ('a, 'b) either = Left of 'a | Right of 'b |
---|
604 | |
---|
605 | |
---|
606 | fun predictor tol (h,ys) = |
---|
607 | let open Real |
---|
608 | val e = Vector.foldl (fn (y,ax) => Real.+ ((abs y),ax)) 0.0 ys |
---|
609 | in |
---|
610 | if e < lb |
---|
611 | then Right (1.414*h) (* step too small, accept but grow *) |
---|
612 | else (if e < ub |
---|
613 | then Right h (* step just right *) |
---|
614 | else Left (0.5*h)) (* step too large, reject and shrink *) |
---|
615 | end |
---|
616 | |
---|
617 | |
---|
618 | exception ConvergenceError |
---|
619 | |
---|
620 | fun integrate (f,x,ys,h,i) = |
---|
621 | let open Real |
---|
622 | val xn = x+h |
---|
623 | val stepper = make_stepper f |
---|
624 | fun recur (hs,x,ys,it) = |
---|
625 | let |
---|
626 | val hf = xn-x |
---|
627 | in |
---|
628 | if Real.>= (hs, hf) |
---|
629 | then (let val (stn, etn) = (stepper hf) (x,ys) |
---|
630 | in {tstep=hs, stn=stn} end) |
---|
631 | else (let val (stn, etn) = (stepper hs) (x,ys) |
---|
632 | in case predictor tol (hs,etn) of |
---|
633 | Left bad => (if Int.<(it,maxit) |
---|
634 | then recur (bad,x,ys,Int.+(it,1)) |
---|
635 | else raise ConvergenceError) |
---|
636 | | Right good => recur (good,x+hs,stn,0) |
---|
637 | end) |
---|
638 | end |
---|
639 | in |
---|
640 | recur (h,x,ys,0) |
---|
641 | end |
---|
642 | |
---|
643 | |
---|
644 | EOF |
---|
645 | |
---|
646 | )) |
---|
647 | (else |
---|
648 | `( |
---|
649 | ("val " ,solver ": (real vector) stepper1 = make_" ,solver "()" ,nl) |
---|
650 | ("fun make_stepper (deriv) = " ,solver " (scaler,summer,deriv)" ,nl) |
---|
651 | ("fun integrate (f,x: real,y: real vector,h,i) = ((make_stepper f) h) (x,y)" ,nl) |
---|
652 | )) |
---|
653 | )) |
---|
654 | ) |
---|
655 | )) |
---|
656 | |
---|
657 | |
---|
658 | (define (name/Octave s) |
---|
659 | (let ((cs (string->list (->string s)))) |
---|
660 | (let loop ((lst (list)) (cs cs)) |
---|
661 | (cond ((null? cs) (string->symbol (list->string (reverse lst)))) |
---|
662 | ((null? (cdr cs)) |
---|
663 | (let ((c (car cs))) |
---|
664 | (if (or (char-alphabetic? c) (char-numeric? c)) |
---|
665 | (loop (cons c lst) (cdr cs)) |
---|
666 | (loop (append (reverse (string->list (->string (gensym 't)))) lst) (cdr cs)) |
---|
667 | ))) |
---|
668 | (else |
---|
669 | (let* ((c (car cs)) |
---|
670 | (c1 (cond ((or (char-alphabetic? c) (char-numeric? c) |
---|
671 | (char=? c #\_)) c) |
---|
672 | (else #\_)))) |
---|
673 | (loop (cons c1 lst) (cdr cs)))))))) |
---|
674 | |
---|
675 | |
---|
676 | (define (binding->Octave x . rest) |
---|
677 | (let-optionals rest ((inner? #f)) |
---|
678 | (cases binding x |
---|
679 | (B:Val (name v) |
---|
680 | (if inner? |
---|
681 | (list (name/Octave name) " = " (value->Octave v) ) |
---|
682 | (list (name/Octave name) " = " (value->Octave v) ";" nl)))))) |
---|
683 | |
---|
684 | (define (expr->Octave x . rest) |
---|
685 | (let-optionals rest ((inner? #f)) |
---|
686 | (cases expr x |
---|
687 | |
---|
688 | (E:Ife (test ift iff) |
---|
689 | (list "if (" (value->Octave test) ") " nl |
---|
690 | (expr->Octave ift ) nl |
---|
691 | " else " (expr->Octave iff) " endif" nl)) |
---|
692 | |
---|
693 | (E:Let (bnds body) |
---|
694 | (if inner? |
---|
695 | (list #\{ |
---|
696 | (intersperse (append (map (lambda (x) (list "(" (binding->Octave x #t) ")")) bnds) |
---|
697 | (list "(" (expr->Octave body #t) ")")) |
---|
698 | ", ") |
---|
699 | #\} |
---|
700 | "{" (+ 1 (length bnds)) "}" nl) |
---|
701 | (list (map (lambda (x) (list (binding->Octave x) nl)) bnds) nl |
---|
702 | (expr->Octave body) nl))) |
---|
703 | |
---|
704 | (E:Ret (v) (value->Octave v)) |
---|
705 | |
---|
706 | (E:Noop () (list)) |
---|
707 | ))) |
---|
708 | |
---|
709 | |
---|
710 | (define (value->Octave v) |
---|
711 | (cases value v |
---|
712 | (V:C (v) v) |
---|
713 | (V:Var (name) (name/Octave name)) |
---|
714 | (V:Rec (lst) |
---|
715 | (list "struct (" (intersperse (map (lambda (nv) (list #\" (name/Octave (first nv)) #\" ", " |
---|
716 | (value->Octave (cadr nv)))) lst) ", ") ")")) |
---|
717 | (V:Sel (field v) |
---|
718 | (list (value->Octave v) "." (name/Octave field))) |
---|
719 | (V:Vec (lst) |
---|
720 | (let ((n (length lst))) |
---|
721 | (list "([" (intersperse (map (lambda (v) (value->Octave v)) lst) ", ") "])"))) |
---|
722 | (V:Sub (index v) |
---|
723 | (list (value->Octave v) "((" index ")+1" ")")) |
---|
724 | (V:Fn (args body) |
---|
725 | (list "(@ (" (intersperse (map name/Octave args) ",") ") " |
---|
726 | (expr->Octave body #t) ")")) |
---|
727 | (V:Op (name args) |
---|
728 | (let* ((infix? (case name |
---|
729 | ((+ - * / >= > < <=) #t) |
---|
730 | (else #f))) |
---|
731 | (op name)) |
---|
732 | (cond ((null? args) |
---|
733 | (case name |
---|
734 | ((NONE) (list name)) |
---|
735 | (else (list name "()")))) |
---|
736 | ((null? (cdr args)) |
---|
737 | (list op "(" (value->Octave (car args)) ")")) |
---|
738 | ((and infix? (null? (cddr args))) |
---|
739 | (list "(" (value->Octave (car args)) ")" op "(" (value->Octave (cadr args)) ")")) |
---|
740 | (infix? |
---|
741 | (let ((op (case op ((+) 'sum) ((*) 'prod) (else op)))) |
---|
742 | (list op "([" (intersperse (map value->Octave args) ",") "])"))) |
---|
743 | (else |
---|
744 | (list op "(" (intersperse (map value->Octave args) ",") ")"))))) |
---|
745 | (V:Ifv (test ift iff) |
---|
746 | (list "(ifv (" (value->Octave test) ", " |
---|
747 | (value->Octave ift ) ", " |
---|
748 | (value->Octave iff) "))")) |
---|
749 | |
---|
750 | )) |
---|
751 | |
---|
752 | |
---|
753 | (define (prelude/Octave #!key (solver 'lsode)) |
---|
754 | `( |
---|
755 | #<<EOF |
---|
756 | |
---|
757 | function res = TRSA(v) |
---|
758 | res = struct ("TRSA",v); |
---|
759 | end |
---|
760 | |
---|
761 | function res = TRSB(v) |
---|
762 | res = struct ("TRSB",v); |
---|
763 | end |
---|
764 | |
---|
765 | function res = tsCase(fa,fb,v) |
---|
766 | if (isfield (v, "TRSA")) |
---|
767 | res = fa(getfield(v,"TRSA")); |
---|
768 | else |
---|
769 | res = fb(getfield(v,"TRSB")); |
---|
770 | endif |
---|
771 | end |
---|
772 | |
---|
773 | function res = TRSC(f,fk,e) |
---|
774 | res = struct ("TRSC",[f,fk,e]); |
---|
775 | end |
---|
776 | |
---|
777 | function res = trfOf(x) |
---|
778 | res = getfield(x,"TRC")(1); |
---|
779 | end |
---|
780 | |
---|
781 | function res = trfkOf(x) |
---|
782 | res = getfield(x,"TRC")(2); |
---|
783 | end |
---|
784 | |
---|
785 | function res = treOf(x) |
---|
786 | res = getfield(x,"TRC")(3); |
---|
787 | end |
---|
788 | |
---|
789 | function res = NONE() |
---|
790 | res = struct (); |
---|
791 | end |
---|
792 | |
---|
793 | function res = SOME(v) |
---|
794 | res = struct ("SOME",v); |
---|
795 | end |
---|
796 | |
---|
797 | function res = swap (x,v) |
---|
798 | if (isfield (v, "SOME")) |
---|
799 | res = getfield(v,"SOME"); |
---|
800 | else |
---|
801 | res = x; |
---|
802 | endif |
---|
803 | end |
---|
804 | |
---|
805 | function res = equal (x, y) |
---|
806 | res = (x==y); |
---|
807 | end |
---|
808 | |
---|
809 | function res = signalOf (v) |
---|
810 | if (not(v)) |
---|
811 | error ("empty signal") |
---|
812 | else |
---|
813 | res = v; |
---|
814 | endif |
---|
815 | end |
---|
816 | |
---|
817 | function res = neg (x) |
---|
818 | res = -x; |
---|
819 | end |
---|
820 | |
---|
821 | function res = ifv (b,x,y) |
---|
822 | if (b) |
---|
823 | res = x; |
---|
824 | else |
---|
825 | res = y; |
---|
826 | endif |
---|
827 | end |
---|
828 | |
---|
829 | EOF |
---|
830 | ,(if (not solver) |
---|
831 | `("function res = integrate (f,x,y,h,i)" ,nl |
---|
832 | " res = f(x,y); " ,nl end ,nl |
---|
833 | ,nl) |
---|
834 | `(( |
---|
835 | #<<EOF |
---|
836 | |
---|
837 | EOF |
---|
838 | ) |
---|
839 | ("global " ,solver #\; ,nl) |
---|
840 | (,solver = #\@ ,(case solver |
---|
841 | ((rk3) 'ode23) ((rk4b rk4a) 'ode45) ((rkf45) 'ode54) |
---|
842 | (else 'lsode)) #\; ,nl) |
---|
843 | |
---|
844 | ,(cond |
---|
845 | |
---|
846 | ((member solver '(rk3 rk4b rk4a rkf45)) |
---|
847 | |
---|
848 | `( |
---|
849 | #<<EOF |
---|
850 | global reltol abstol dt; |
---|
851 | reltol = 0.1; |
---|
852 | abstol = 0.01; |
---|
853 | dt = 0.001; |
---|
854 | |
---|
855 | global P; |
---|
856 | P = odeset ("RelTol" , reltol , "AbsTol" , abstol , "MaxStep" , 1 , "InitialStep" , dt) ; |
---|
857 | |
---|
858 | EOF |
---|
859 | "function res = integrate (f,x,y,h,i)" ,nl |
---|
860 | " global P " ,solver ";" ,nl |
---|
861 | " [t,y] = " ,solver "(f,[x,x+h],y,P); " ,nl |
---|
862 | " res = y(2);" ,nl end ,nl |
---|
863 | )) |
---|
864 | (else |
---|
865 | |
---|
866 | `( |
---|
867 | #<<EOF |
---|
868 | lsode_options("relative tolerance",1e-4); |
---|
869 | lsode_options("absolute tolerance",1e-4); |
---|
870 | |
---|
871 | |
---|
872 | EOF |
---|
873 | "function res = integrate (f,x,y0,h,i)" ,nl |
---|
874 | " global " ,solver ";" ,nl |
---|
875 | " y = " ,solver "(@(yvec,t) f(t,yvec),y0,[x,x+h]); " ,nl |
---|
876 | " res = y(2);" ,nl end ,nl |
---|
877 | ) |
---|
878 | )) |
---|
879 | )) |
---|
880 | )) |
---|
881 | |
---|
882 | |
---|
883 | |
---|
884 | ) |
---|