Changeset 30942 in project


Ignore:
Timestamp:
05/30/14 07:35:18 (6 years ago)
Author:
Ivan Raikov
Message:

signal-diagram/flsim: synchronizing changes related to adaptive time step solvers

Location:
release/4
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • release/4/flsim/trunk/flsim.scm

    r30939 r30942  
    317317    result))
    318318 
    319 (define (prelude/scheme #!key (solver 'rk4b) (integral-index 0))
     319(define (prelude/scheme #!key (solver 'rk4b) (random #f) (integral-index 0))
    320320
    321321`(
     
    469469         (V:Vec     (lst)
    470470                    (let ((n (length lst)))
    471                       (list "([" (intersperse (map (lambda (v) (value->ML v)) lst) ", ") "])")))
     471                      (list "(Vector.fromList [" (intersperse (map (lambda (v) (value->ML v)) lst) ", ") "])")))
    472472         (V:Sub     (index v)
    473                     (list "List.nth (" (value->ML v) ", " index ")"))
     473                    (list "Unsafe.Vector.sub (" (value->ML v) ", " index ")"))
    474474         (V:Fn      (args body)
    475475                    (list "(fn (" (intersperse (map name/ML args) ",") ") => "
     
    480480                                     (else #f)))
    481481                           (op (if infix? (list "(op " name ")") name)))
    482                       (let ((name (name/ML name)))
     482                      (let ((name1 (name/ML name)))
    483483                        (cond ((null? args)
    484484                               (case name
    485                                  ((NONE)  (list name))
    486                                  (else    (list name "()"))))
     485                                 ((NONE)  (list name1))
     486                                 (else    (list name1 "()"))))
    487487                              ((null? (cdr args))
    488488                               (list "(" op " " (value->ML (car args)) ")"))
     
    500500
    501501
    502 (define (prelude/ML  #!key (solver 'rk4b) )
     502(define (prelude/ML  #!key (solver 'rk4b) (random #f))
    503503`(
    504504 #<<EOF
     
    509509open Math
    510510open RungeKutta
    511 open RandomMTZig
    512511
    513512datatype ('b,'c) trs = TRSA of 'b | TRSB of 'c
     
    533532end
    534533
     534fun 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
    535542exception EmptySignal
     543
    536544val neg = (op ~)
    537545val swap = fn (x,v) => (case v of NONE => x | SOME v => v)
     
    539547val signalOf = fn (v) => (case v of NONE => raise EmptySignal | SOME v => v)
    540548val heaviside = fn (v) => (if Real.< (v, 0.0) then 0.0 else 1.0)
     549EOF
     550
     551,(if random
     552#<<EOF
    541553
    542554fun RandomInit () = RandomMTZig.fromEntropy()
     
    566578
    567579EOF
     580"")
    568581
    569582,(if (not solver)
    570583     `(("(* dummy solver; returns only the computed derivatives *)")
    571        ("fun integrate (f,x: real,y: real list,h,i) = (f (x,y))" ,nl)
     584       ("fun integrate (f,x: real,y: real vector,h,i) = (f (x,y))" ,nl)
    572585       )
    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)
     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)
    577588       . ,(case solver 
    578589            ;; adaptive solvers
    579590            ((rkhe rkbs rkf45 rkck rkdp rkf78 rkv65)
    580591             `(
     592               ("val " ,solver ": (real vector) stepper2 = make_" ,solver "()" ,nl)
     593               ("fun make_stepper (deriv) = " ,solver " (scaler,summer,deriv)" ,nl)
    581594#<<EOF
     595
     596
    582597
    583598val tol = Real.Math.pow (10.0, ~6.0)
    584599val lb = 0.5 * tol
    585600val ub = 1.0 * tol
    586 
     601val maxit = 100
    587602
    588603datatype ('a, 'b) either = Left of 'a | Right of 'b
     
    591606fun predictor tol (h,ys) =
    592607  let open Real
    593       val e = foldl (fn (y,ax) => Real.+ ((abs y),ax)) 0.0 ys
     608      val e = Vector.foldl (fn (y,ax) => Real.+ ((abs y),ax)) 0.0 ys
    594609  in
    595610      if e < lb
     
    601616
    602617
     618exception ConvergenceError
     619
    603620fun integrate (f,x,ys,h,i) =
    604621  let open Real
     622      val xn = x+h
    605623      val stepper = make_stepper f
    606       fun recur (hs) =
     624      fun recur (hs,x,ys,it) =
    607625          let
    608               val (stn, etn) = (stepper hs) (x,ys)
     626              val hf = xn-x
    609627          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)
     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 (h,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)
    616638          end
    617639  in
    618       recur h
     640      recur (h,x,ys,0)
    619641  end
     642
    620643
    621644EOF
     
    624647            (else
    625648             `(
    626                ("fun integrate (f,x: real,y: real list,h,i) = ((make_stepper f) h) (x,y)" ,nl)
     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)
    627652               ))
    628653            ))
  • release/4/signal-diagram/trunk/examples/Izhikevich03.scm

    r28566 r30942  
    1010(define (Izhikevich03:construct t V U k1 k2 k3 theta a b c d spike tspike Isyn )
    1111
    12   `((f . ,(let* ((subthreshold-eq         (ODE h t
     12  `((f . ,(let* ((subthreshold-eq         (ODE (variable h) t
    1313                                               (V  (+ (* k1 V V)  (* k2 V) k3 (neg U) Isyn))
    1414                                               (U  (* a (- (* b V) U)))))
     
    6060
    6161
    62 (define (codegen name model #!key (language 'scheme) (pre #t) (post #t) (solver 'rk3))
     62(define (codegen name model #!key (language 'scheme) (random #f) (pre #t) (post #t) (solver 'rkdp))
    6363  (let* ((f         (alist-ref 'f model))
    6464         (initial   (alist-ref 'initial model)))
     
    6666      ((octave Octave) (codegen/Octave name (construct f) initial: initial pre: pre solver: solver))
    6767      ((scheme Scheme) (codegen/scheme name (construct f) initial: initial pre: pre solver: solver))
    68       ((ML)     (codegen/ML name (construct f) initial: initial pre: pre post: post solver: solver)))
     68      ((ML)     (codegen/ML name (construct f) initial: initial random: random pre: pre post: post solver: solver)))
    6969    ))
    7070
     
    161161
    162162
    163 (define models `((RS  . ,RS) (IB . ,IB) (CH . ,CH)
    164                  (FS . ,FS)  (RZ . ,RZ) (LTS . ,LTS) ))
     163(define models `((RS  . ,RS) (IB . ,IB) (CH . ,CH) (FS . ,FS)  (RZ . ,RZ) (LTS . ,LTS) ))
    165164
    166165 
     166(with-output-to-file "Izhikevich03_solver.sml"
     167  (lambda ()
     168    (let recur ((models models) (pre #t))
     169      (if (pair? models)
     170          (let* ((name.model (car models))
     171                 (name  (car name.model))
     172                 (model (cdr name.model)))
     173            (codegen name model language: 'ML solver: 'rkdp random: #f pre: pre post: (null? (cdr models)))
     174            (recur (cdr models) #f)
     175            )))
     176    ))
     177
    167178(with-output-to-file "Izhikevich03_solver.m"
    168179  (lambda ()
     
    190201
    191202
    192  
    193 (with-output-to-file "Izhikevich03_solver.sml"
    194   (lambda ()
    195     (let recur ((models models) (pre #t))
    196       (if (pair? models)
    197           (let* ((name.model (car models))
    198                  (name  (car name.model))
    199                  (model (cdr name.model)))
    200             (codegen name model language: 'ML pre: pre post: (null? (cdr models)))
    201             (recur (cdr models) #f)
    202             )))
    203     ))
    204 
    205 
     203
  • release/4/signal-diagram/trunk/examples/Izhikevich03_run.mlb

    r23938 r30942  
    1 (* Compile with mlton -mlb-path-var 'RK_LIB $(CHICKEN_HOME)/share/chicken/flsim/sml-lib' *)
     1(* Compile with mlton -mlb-path-var 'RK_LIB $(CHICKEN_HOME)/share/chicken/flsim/sml-lib/rk' *)
    22
    33$(SML_LIB)/basis/basis.mlb
    4 $(RK_LIB)/rk/rk.mlb
     4$(SML_LIB)/basis/unsafe.mlb
     5$(RK_LIB)/rk.mlb
     6
    57local
    68        Izhikevich03_solver.sml
  • release/4/signal-diagram/trunk/examples/Izhikevich03_run.sml

    r23917 r30942  
    4848
    4949val _ = (start (Model.RS_initial,240.0,Model.RS,make_pulse (40.0,200.0,0.0,10.0)))
     50(*
    5051val _ = (start (Model.IB_initial,240.0,Model.IB,make_pulse (40.0,200.0,0.0,10.0)))
    5152val _ = (start (Model.CH_initial,240.0,Model.CH,make_pulse (40.0,200.0,0.0,10.0)))
    5253val _ = (start (Model.FS_initial,240.0,Model.FS,make_pulse (40.0,200.0,0.0,10.0)))
    5354val _ = (start (Model.LTS_initial,240.0,Model.LTS,make_pulse (40.0,200.0,0.0,10.0)))
     55*)
  • release/4/signal-diagram/trunk/examples/Morris-Lecar81.scm

    r30940 r30942  
    1515 
    1616
    17   `((f . ,(DAE h t
     17  `((f . ,(DAE (variable h) t
    1818               (minf (v) (* 0.5 (+ 1 (tanh (/ (- v v1) v2)))))
    1919               (winf (v) (* 0.5 (+ 1 (tanh (/ (- v v3) v4)))))
     
    6464                     ((scheme Scheme) codegen/scheme)
    6565                     ((ML) codegen/ML))))
    66       (codegen name (construct f adaptive-integral-method:
    67                                (case solver
    68                                      ((rkhe rkbs rkf45 rkck rkdp rkf78 rkv65) #t)
    69                                      (else #f)))
     66      (codegen name (construct f)
    7067               initial: initial solver: solver ))
    7168    ))
  • release/4/signal-diagram/trunk/examples/Morris-Lecar81_run.mlb

    r23938 r30942  
    1 (* Compile with mlton -mlb-path-var 'RK_LIB $(CHICKEN_HOME)/share/chicken/flsim/sml-lib' *)
     1(* Compile with mlton -mlb-path-var 'RK_LIB $(CHICKEN_HOME)/share/chicken/flsim/sml-lib/rk' *)
    22
    33$(SML_LIB)/basis/basis.mlb
    4 $(RK_LIB)/rk/rk.mlb
     4$(SML_LIB)/basis/unsafe.mlb
     5$(RK_LIB)/rk.mlb
     6
    57local
    68        Morris-Lecar81_solver.sml
  • release/4/signal-diagram/trunk/signal-diagram-dynamics.scm

    r30940 r30942  
    8080          (else  expr)))
    8181
    82   (let ((dqs (filter-map (lambda (x) (cond ((= 2 (length x)) (car x))
     82  (let ((varh (case (car h)
     83                ((variable) #t)
     84                (else #f)))
     85        (hname (cadr h))
     86
     87        (dqs (filter-map (lambda (x) (cond ((= 2 (length x)) (car x))
    8388                                            (else #f)))
    8489                          eqs))
     
    125130                 eqs)))
    126131
    127       (let ((du (ACTUATE dqs (INTEGRALH indep dqs h dfs))))
     132      (let ((du (ACTUATE (if varh (cons hname dqs) dqs)
     133                         (INTEGRAL indep dqs h dfs))))
    128134
    129135        (make-relation (zip aqs ads afs) ;;(SEQUENCE du (make-assign-system `((,indep (+ ,indep ,h))))))
    130                        (SEQUENCE du (ACTUATE (list indep h) (PURE (make-function (list indep h) `(+ ,indep ,h)))))
    131                        )
     136                       (UNION du (ACTUATE (list indep)
     137                                          (PURE (make-function
     138                                                 (list indep hname)
     139                                                 `(+ ,indep ,hname))))))
    132140      ))
    133141    ))
     
    142150 
    143151(define (make-ode-system h indep eqs)
    144   (let ((deps (map car eqs))
     152  (let ((varh (case (car h)
     153                ((variable) #t)
     154                (else #f)))
     155        (hname (cadr h))
     156        (deps (map car eqs))
    145157        (rhss (map cadr eqs)))
     158
    146159    (let ((fs (map
    147160               (lambda (rhs)                   
     
    150163               rhss)))
    151164
    152       (let ((u (cond ((or (null? fs) (null? deps)) (error 'make-ode-system "empty list of equations"))
    153                      ((null? (cdr fs))
    154                       (let ((d (car deps)))
    155                         (ACTUATE (list d) (INTEGRALH indep (list d) h (list (car fs))))
    156                         ))
    157                      (else
    158                       (ACTUATE deps (INTEGRALH indep deps h fs)))
    159                      )))
    160 
    161         (SEQUENCE u (ACTUATE (list indep h) (PURE (make-function (list indep h) `(+ ,indep ,h)))))
     165      (let ((u (cond
     166                ((or (null? fs) (null? deps))
     167                 (error 'make-ode-system "empty list of equations"))
     168
     169                ((null? (cdr fs))
     170                 (let ((d (car deps)))
     171                   (ACTUATE (if varh (list hname d) (list d))
     172                            (INTEGRAL indep (list d) h (list (car fs))))
     173                   ))
     174
     175                (else
     176                 (ACTUATE (if varh (cons hname deps) deps)
     177                          (INTEGRAL indep deps h fs)))
     178                )))
     179
     180        (UNION u (ACTUATE (list indep)
     181                          (PURE (make-function
     182                                 (list indep hname)
     183                                 `(+ ,indep ,hname)))))
    162184
    163185        ))
  • release/4/signal-diagram/trunk/signal-diagram.scm

    r30940 r30942  
    181181       (function? (caddr r))))
    182182
     183(define (hspec? h)
     184  (and (list? h)
     185       (case (car h)
     186         ((variable fixed var fix) #t)
     187         (else #f))
     188       (symbol? (cadr h))))
    183189
    184190(define-datatype diagram diagram?
     
    200206  (TRANSIENT    (f diagram?) (g diagram?) (e symbol?) )
    201207  (ON           (f diagram?) (e symbol?) )
    202   (INTEGRAL     (i symbol?) (d symbol-list?) (h (lambda (x) (or (symbol? x) (number? x))))
    203                 (f function-list?))
     208  (INTEGRAL     (i symbol?)
     209                (d symbol-list?) (h hspec?)
     210                (f function-list?) )
    204211  )
    205212
     
    255262       (lambda (s) (if (prim? fd) outputs (list name)))
    256263       ;; in
    257        (lambda (s) s)
     264       (lambda (s) (if (function? fd)
     265                       (lset-intersection eq? (function-formals fd) s)
     266                       s))
    258267       ;; out
    259268       (lambda (s) (if (prim? fd) outputs (list name))))
     
    598607               (rv  (gensym 'sequence))
    599608               )
     609
     610
    600611
    601612          (make-codegen
     
    12461257
    12471258
    1248 (define (sf-integral0 x ys h fs method)
    1249 
    1250   (let* ((xn     (gensym (string->symbol (s+ x "+h"))))
    1251          (yis    (if (and (adaptive-integral) (symbol? h))
    1252                      (list-tabulate (length ys) (lambda (i) (+ 1 i)))
    1253                      (list-tabulate (length ys) (lambda (i) i))))
     1259(define (sf-integral x ys h fs)
     1260
     1261  (let* ((varh   (case (car h)
     1262                   ((variable) #t)
     1263                   (else #f)))
     1264         (hname  (cadr h))
     1265         (xn     (gensym (string->symbol (s+ x "+h"))))
     1266         (yis    (list-tabulate (length ys) (lambda (i) i)))
    12541267         (yns    (map (lambda (y) (gensym (string->symbol (s+ y "(" xn ")")))) ys))
    12551268         (ynvs   (map (lambda (yn) (gensym (string->symbol (s+ yn "v")))) yns))
     
    12681281
    12691282        ;; gen
    1270         (lambda (s) (if (and (adaptive-integral) (symbol? h))
    1271                         (cons h yns) yns))
     1283        (lambda (s) (if varh (cons hname yns) yns))
    12721284       
    12731285        ;; kill
     
    12801292                               (lset-union eq?
    12811293                                           (concatenate fs-formals)
    1282                                            (append (if (symbol? h) (list h) '())
     1294                                           (append (list hname)
    12831295                                                   (cons x ys))))))
    12841296            x))
    12851297       
    12861298        ;; out
    1287         (lambda (s) (if (and (adaptive-integral) (symbol? h))
    1288                         (cons h yns) yns))
     1299        (lambda (s) (if varh (cons hname yns) yns))
    12891300        )
    12901301       
     
    13041315
    13051316                  (tstep     (if (symbol? h)
    1306                                  (select-signal 'integral1 h env)
    1307                                  (V:C h)))
     1317                                 (select-signal 'integral1 hname env)
     1318                                 (V:C hname)))
    13081319                 
    13091320                  (fenv      (list->cgenenv 'integral2 (concatenate fs-formals)
     
    13131324                  )
    13141325
     1326
    13151327             (make-codegen
    13161328
    13171329              rv2
    13181330
    1319               ((lambda (env) (if (and (adaptive-integral) (symbol? h))
    1320                                  (cons (cons h rv2) env) env))
     1331              ((lambda (env) (if varh (cons (cons hname rv2) env) env))
    13211332               (map (lambda (s) (cons s rv2)) yns))
    13221333
     
    13411352                                   idxv
    13421353                                   )))
    1343 
    1344 
    1345                 (B:Val rv2
    1346                        (V:Rec ((lambda (flds)
    1347                                  (if (and (adaptive-integral) (symbol? h))
    1348                                      (cons `(,h ,(V:Sub 0 (V:Var rv1))) flds)
    1349                                      flds))
    1350                                (map (lambda (yn yi) `(,yn ,(V:Sub yi (V:Var rv1)) )) yns yis)) ))
    13511354                )
     1355
     1356                (if varh
     1357                    (let ((stn1 (gensym 'stn)))
     1358                      (list
     1359                       (B:Val stn1 (V:Sel 'stn (V:Var rv1)))
     1360                       (B:Val rv2
     1361                              (V:Rec
     1362                               (cons `(,hname ,(V:Sel 'tstep (V:Var rv1)))
     1363                                     (map (lambda (yn yi) `(,yn ,(V:Sub yi (V:Var stn1)))) yns yis))))))
     1364                     (list (B:Val rv2
     1365                                  (V:Rec (map (lambda (yn yi) `(,yn ,(V:Sub yi (V:Var rv1)) )) yns yis)) )))
    13521366               ))
    13531367             ))
     
    13681382
    13691383
    1370 (define (sf-integralh x y h f)
    1371   (sf-integral0 x y h f))
    1372 
    1373 
    1374 (define (construct d #!key (method 'rk3))
     1384
     1385(define (construct d)
    13751386  (integral-index 0)
    13761387  (dynvector-clear! integral-events 0)
    1377   (construct1 d method))
    1378 
    1379 
    1380 (define (construct1 d method)
     1388  (construct1 d))
     1389
     1390
     1391(define (construct1 d)
    13811392  (let recur ((d d))
    13821393    (cases diagram d
     
    13941405           (TRANSIENT (f g e)          (sf-transient (recur f) (recur g) e))
    13951406           (ON (f e)                   (sf-on (recur f) e))
    1396            (INTEGRAL  (x ys h fs)      (sf-integralh method x ys h fs))
     1407           (INTEGRAL (x ys h fs)       (sf-integral x ys h fs))
    13971408           )))
    13981409
     
    15491560      )))
    15501561
    1551 (define (codegen/ML name f #!key (initial #f) (pre #t) (post #t) (solver 'rk4b))
     1562(define (codegen/ML name f #!key (initial #f) (random #f) (pre #t) (post #t) (solver 'rk4b))
    15521563
    15531564  (if (and solver (not (member solver '(rkfe rk3 rk4a rk4b rkhe rkbs rkf45 rkck rkdp rkf78 rkv65))))
    15541565      (error 'codegen/ML "unknown solver" solver))
     1566
    15551567
    15561568  (let ((dfe (sfarrow-dfe f)))
     
    15651577           )
    15661578
    1567       (if pre (print-fragments (prelude/ML solver: solver)))
     1579      (if pre (print-fragments (prelude/ML solver: solver random: random)))
    15681580
    15691581      (print-fragments (list (map (lambda (x)
Note: See TracChangeset for help on using the changeset viewer.