source: project/release/4/flsim/trunk/flsim.scm @ 30939

Last change on this file since 30939 was 30939, checked in by Ivan Raikov, 6 years ago

flsim: support for variable timestep integration

File size: 24.8 KB
Line 
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) (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
382EOF
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 "([" (intersperse (map (lambda (v) (value->ML v)) lst) ", ") "])")))
472         (V:Sub     (index v) 
473                    (list "List.nth (" (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 ((name (name/ML name)))
483                        (cond ((null? args) 
484                               (case name 
485                                 ((NONE)  (list name))
486                                 (else    (list name "()"))))
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) )
503`(
504 #<<EOF
505structure Model = 
506struct
507
508open Real
509open Math
510open RungeKutta
511open RandomMTZig
512
513datatype ('b,'c) trs = TRSA of 'b | TRSB of 'c
514datatype ('a,'b,'c) trc = TRC of ((('a -> (('b,'c) trs))) * 
515                                  (('a -> (('b,'c) trs))) * 
516                                  ((('b,'c) trs) -> bool))
517         
518fun tsCase (fa,fb,x) = case x of TRSA a => (fa a) | TRSB b => (fb b)
519fun trfOf x = case x of TRC (f,fk,e) => f
520fun trfkOf x = case x of TRC (f,fk,e) => fk
521fun treOf x = case x of TRC (f,fk,e) => e
522
523fun putStrLn str = 
524  (TextIO.output (TextIO.stdOut, str);
525   TextIO.output (TextIO.stdOut, "\n"))
526
527fun putStr str = (TextIO.output (TextIO.stdOut, str))
528
529fun showReal n = 
530let open StringCvt
531in
532(if n < 0.0 then "-" else "") ^ (fmt (FIX (SOME 12)) (abs n))
533end
534
535exception EmptySignal
536val neg = (op ~)
537val swap = fn (x,v) => (case v of NONE => x | SOME v => v) 
538val equal = fn (x,y) => (x = y) 
539val signalOf = fn (v) => (case v of NONE => raise EmptySignal | SOME v => v) 
540val heaviside = fn (v) => (if Real.< (v, 0.0) then 0.0 else 1.0)
541
542fun RandomInit () = RandomMTZig.fromEntropy()
543
544val RandomState = RandomInit ()
545
546fun random_uniform () = RandomMTZig.randUniform RandomState
547
548
549fun PoissonInit () =
550 let
551     val zt = RandomMTZig.ztnew()
552     val st = RandomMTZig.fromEntropy()
553 in
554     {st=st,zt=zt}
555 end
556
557fun PoissonStep (rate,t,h,st,zt) =
558 let
559    val rv     = RandomMTZig.randPoisson (rate*0.001*h,st,zt) 
560    val spike' = Real.> (rv,0.0)
561    val spikeCount' = if spike' then rv else 0.0
562  in
563     {t=t+h,spike=spike',spikeCount=spikeCount',st=st,zt=zt}
564  end
565
566
567EOF
568
569,(if (not solver)
570     `(("(* dummy solver; returns only the computed derivatives *)")
571       ("fun integrate (f,x: real,y: real list,h,i) = (f (x,y))" ,nl)
572       )
573     `(("val summer = fn (a,b) => (ListPair.map (fn (x,y) => x+y) (a,b))" ,nl)
574       ("val scaler = fn(a,lst) => (map (fn (x) => a*x) lst)" ,nl)
575       ("val " ,solver ": (real list) stepper1 = make_" ,solver "()" ,nl)
576       ("fun make_stepper (deriv) = " ,solver " (scaler,summer,deriv)" ,nl)
577       . ,(case solver 
578            ;; adaptive solvers
579            ((rkhe rkbs rkf45 rkck rkdp rkf78 rkv65)
580             `(
581#<<EOF
582
583val tol = Real.Math.pow (10.0, ~6.0)
584val lb = 0.5 * tol
585val ub = 1.0 * tol
586
587
588datatype ('a, 'b) either = Left of 'a | Right of 'b
589
590
591fun predictor tol (h,ys) =
592  let open Real
593      val e = foldl (fn (y,ax) => Real.+ ((abs y),ax)) 0.0 ys
594  in 
595      if e < lb 
596      then Right (1.414*h)              (* step too small, accept but grow *)
597      else (if e < ub 
598            then Right h                (* step just right *)
599            else Left (0.5*h))          (* step too large, reject and shrink *)
600  end
601
602
603fun integrate (f,x,ys,h,i) =
604  let open Real
605      val stepper = make_stepper f
606      fun recur (hs) =
607          let
608              val (stn, etn) = (stepper hs) (x,ys)
609          in
610              if hs < 1.0 andalso Real.== (hs, hf)
611              then (hs, xn, stn)
612              else
613                  (case predictor tol (h,etn) of
614                       Left bad => recur bad
615                     | Right good => good::stn)
616          end
617  in
618      recur h
619  end
620
621EOF
622
623               ))
624            (else
625             `(
626               ("fun integrate (f,x: real,y: real list,h,i) = ((make_stepper f) h) (x,y)" ,nl)
627               ))
628            ))
629       )
630))
631
632
633(define (name/Octave s)
634  (let ((cs (string->list (->string s))))
635    (let loop ((lst (list)) (cs cs))
636      (cond ((null? cs) (string->symbol (list->string (reverse lst))))
637            ((null? (cdr cs))
638             (let ((c (car cs)))
639               (if (or (char-alphabetic? c) (char-numeric? c))
640                   (loop (cons c lst) (cdr cs))
641                   (loop (append (reverse (string->list (->string (gensym 't)))) lst) (cdr cs))
642                   )))
643            (else
644             (let* ((c (car cs))
645                    (c1 (cond ((or (char-alphabetic? c) (char-numeric? c) 
646                                   (char=? c #\_)) c)
647                              (else #\_))))
648               (loop (cons c1 lst) (cdr cs))))))))
649
650
651(define (binding->Octave x . rest)
652  (let-optionals rest ((inner? #f))
653    (cases binding x
654         (B:Val     (name v)
655                    (if inner?
656                        (list (name/Octave name) " = " (value->Octave v) )
657                        (list (name/Octave name) " = " (value->Octave v) ";" nl))))))
658
659(define (expr->Octave x . rest)
660  (let-optionals rest ((inner? #f))
661    (cases expr x
662
663         (E:Ife     (test ift iff)
664                    (list "if (" (value->Octave test) ") " nl
665                          (expr->Octave ift ) nl
666                          " else " (expr->Octave iff) " endif" nl))
667
668         (E:Let     (bnds body)
669                    (if inner?
670                        (list #\{ 
671                              (intersperse (append (map (lambda (x) (list "(" (binding->Octave x #t) ")")) bnds)
672                                                   (list "(" (expr->Octave body #t) ")"))
673                                           ", ")
674                              #\}
675                              "{" (+ 1 (length bnds)) "}" nl)
676                        (list (map (lambda (x) (list (binding->Octave x) nl)) bnds) nl
677                              (expr->Octave body) nl)))
678                         
679         (E:Ret     (v)  (value->Octave v))
680                   
681         (E:Noop    () (list))
682         )))
683
684
685(define (value->Octave v)
686  (cases value v
687         (V:C       (v) v)
688         (V:Var     (name) (name/Octave name))
689         (V:Rec     (lst) 
690                    (list "struct (" (intersperse (map (lambda (nv) (list #\" (name/Octave (first nv)) #\" ", " 
691                                                                          (value->Octave (cadr nv)))) lst) ", ") ")"))
692         (V:Sel     (field v) 
693                    (list (value->Octave v) "." (name/Octave field)))
694         (V:Vec     (lst) 
695                    (let ((n (length lst)))
696                      (list "([" (intersperse (map (lambda (v) (value->Octave v)) lst) ", ") "])")))
697         (V:Sub     (index v) 
698                    (list (value->Octave v) "((" index ")+1" ")"))
699         (V:Fn      (args body) 
700                    (list "(@ (" (intersperse (map name/Octave args) ",") ") "
701                          (expr->Octave body #t) ")"))
702         (V:Op    (name args)
703                    (let* ((infix? (case name
704                                     ((+ - * / >= > < <=)  #t)
705                                     (else #f)))
706                           (op name))
707                      (cond ((null? args) 
708                             (case name 
709                               ((NONE)  (list name))
710                               (else    (list name "()"))))
711                            ((null? (cdr args))
712                             (list op "(" (value->Octave (car args)) ")"))
713                            ((and infix? (null? (cddr args)))
714                             (list "(" (value->Octave (car args)) ")" op "(" (value->Octave (cadr args)) ")"))
715                            (infix?
716                             (let ((op (case op ((+) 'sum) ((*) 'prod) (else op))))
717                               (list  op "([" (intersperse (map value->Octave args) ",") "])")))
718                            (else
719                             (list  op "(" (intersperse (map value->Octave args) ",") ")")))))
720         (V:Ifv     (test ift iff)
721                    (list "(ifv (" (value->Octave test) ", " 
722                          (value->Octave ift ) ", " 
723                          (value->Octave iff) "))"))
724
725         ))
726
727 
728(define (prelude/Octave #!key (solver 'lsode))
729`(
730#<<EOF
731
732function res = TRSA(v)
733  res = struct ("TRSA",v)
734end
735
736function res = TRSB(v)
737  res = struct ("TRSB",v)
738end
739
740function res = tsCase(fa,fb,v)
741   if (isfield (v, "TRSA"))
742     res = fa(getfield(v,"TRSA"));
743   else
744     res = fb(getfield(v,"TRSB"));
745   endif
746end
747
748function res = TRSC(f,fk,e)
749  res = struct ("TRSC",[f,fk,e])
750end
751
752function res = trfOf(x)
753   res = getfield(x,"TRC")(1);
754end
755
756function res = trfkOf(x)
757   res = getfield(x,"TRC")(2);
758end
759
760function res = treOf(x)
761   res = getfield(x,"TRC")(3);
762end
763
764function res = NONE()
765  res = struct ()
766end
767
768function res = SOME(v)
769  res = struct ("SOME",v)
770end
771
772function res = swap (x,v)
773   if (isfield (v, "SOME"))
774     res = getfield(v,"SOME");
775   else
776     res = x;
777   endif
778end
779
780function res = equal (x, y)
781   res = (x==y);
782end
783
784function res = signalOf (v)
785  if (not(v))
786   error ("empty signal")
787 else
788   res = v;
789 endif
790end
791
792function res = neg (x)
793  res = -x;
794end
795
796function res = ifv (b,x,y)
797  if (b) 
798    res = x;
799  else
800    res = y;
801  endif
802end
803
804EOF
805,(if (not solver)
806     `("function res = integrate (f,x,y,h,i)" ,nl 
807       "  res = f(x,y); " ,nl end ,nl
808       ,nl)
809     `((
810#<<EOF
811
812EOF
813)
814       ("global " ,solver #\; ,nl)
815       (,solver = #\@ ,(case solver
816                         ((rk3) 'ode23) ((rk4b rk4a) 'ode45) ((rkf45) 'ode54) 
817                         (else 'lsode)) #\; ,nl)
818
819       ,(cond
820
821         ((member solver '(rk3 rk4b rk4a rkf45))
822
823               `(
824#<<EOF
825global reltol abstol dt;
826reltol = 0.1;
827abstol = 0.01;
828dt = 0.001;
829
830global P;
831P = odeset ("RelTol" , reltol , "AbsTol" , abstol , "MaxStep" , 1 , "InitialStep" , dt) ;
832
833EOF
834       "function res = integrate (f,x,y,h,i)" ,nl 
835       "  global P " ,solver ";" ,nl
836       "  [t,y] = " ,solver "(f,[x,x+h],y,P); " ,nl 
837       "  res = y(2);" ,nl end ,nl
838       ))
839              (else
840
841               `(
842#<<EOF
843lsode_options("relative tolerance",1e-4);
844lsode_options("absolute tolerance",1e-4);
845
846
847EOF
848       "function res = integrate (f,x,y0,h,i)" ,nl 
849       "  global " ,solver ";" ,nl
850       "  y = " ,solver "(@(yvec,t) f(t,yvec),y0,[x,x+h]); " ,nl 
851       "  res = y(2);" ,nl end ,nl
852       )
853   ))
854 ))
855))
856
857         
858
859)
Note: See TracBrowser for help on using the repository browser.