Changeset 30972 in project


Ignore:
Timestamp:
06/06/14 06:13:11 (6 years ago)
Author:
Ivan Raikov
Message:

signal-diagram: bug fixes in integral variable propagation

Location:
release/4/signal-diagram/trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • release/4/signal-diagram/trunk/examples/Izhikevich03.scm

    r30964 r30972  
    1616                 (threshold-detect        (ASSIGN (spike (- V theta))))
    1717
    18                  (refractory-eq           (ASSIGN (t (+ t h) )
     18                 (refractory-eq           (ASSIGN (t t)
    1919                                                  (h h)
    2020                                                  (V c)
  • release/4/signal-diagram/trunk/signal-diagram-dynamics.scm

    r30970 r30972  
    130130                 eqs)))
    131131
    132       (let ((du (ACTUATE (if varh (append dqs (list hname)) dqs)
     132      (let ((du (ACTUATE (if varh (append (cons indep dqs) (list hname)) dqs)
    133133                         (INTEGRAL indep dqs h dfs))))
    134134
    135135        (make-relation (zip aqs ads afs) ;;(SEQUENCE du (make-assign-system `((,indep (+ ,indep ,h))))))
    136                        (UNION du (ACTUATE (list indep)
    137                                           (PURE (make-function
    138                                                  (list indep hname)
    139                                                  `(+ ,indep ,hname))))))
     136                       du )
    140137      ))
    141138    ))
     
    169166                ((null? (cdr fs))
    170167                 (let ((d (car deps)))
    171                    (ACTUATE (if varh (list d hname) (list d))
     168                   (ACTUATE (if varh (list indep d hname) (list d))
    172169                            (INTEGRAL indep (list d) h (list (car fs))))
    173170                   ))
    174171
    175172                (else
    176                  (ACTUATE (if varh (append deps (list hname)) deps)
     173                 (ACTUATE (if varh
     174                              (append (cons indep deps) (list hname))
     175                              (cons indep deps))
    177176                          (INTEGRAL indep deps h fs)))
    178177                )))
    179178
    180         (UNION u (ACTUATE (list indep)
    181                           (PURE (make-function
    182                                  (list indep hname)
    183                                  `(+ ,indep ,hname)))))
    184 
     179        u
    185180        ))
    186181    ))
  • release/4/signal-diagram/trunk/signal-diagram.scm

    r30970 r30972  
    13021302        (lambda (s) (lset-union eq?
    13031303                               (or (and ee (ee-gen s)) '())
    1304                                (if varh (cons hname yns) yns)))
     1304                               (if varh (cons hname (cons x yns))
     1305                                   (cons x yns))))
    13051306       
    13061307        ;; kill
     
    13201321       
    13211322        ;; out
    1322         (lambda (s) (append yns
     1323        (lambda (s) (append (cons x yns)
    13231324                            (or (and varh (list hname)) '())
    13241325                            (or (and ee (ee-out s)) '())
     
    13561357
    13571358              ((lambda (env) (if varh (cons (cons hname rv2) env) env))
    1358                (map (lambda (s) (cons s rv2)) yns))
     1359               (cons (cons x rv2) (map (lambda (s) (cons s rv2)) yns)))
    13591360
    13601361              (append
     
    14221423
    14231424                (if varh
     1425
    14241426                    (let* ((ysn (gensym 'ysn))
     1427                           (xn  (gensym 'xn))
    14251428                           (retflds (cons `(,hname ,(V:Sel 'h (V:Var rv1)))
    1426                                           (map (lambda (yn yi) `(,yn ,(V:Sub yi (V:Var ysn)))) yns yis))))
     1429                                          (cons `(,x ,(V:Var xn))
     1430                                                (map (lambda (yn yi) `(,yn ,(V:Sub yi (V:Var ysn))))
     1431                                                     yns yis)))))
    14271432                      (list
    14281433                       (B:Val ysn (V:Sel 'ysn (V:Var rv1)))
     1434                       (B:Val xn  (V:Sel 'xn (V:Var rv1)))
    14291435                       (B:Val rv2 (V:Rec retflds))))
     1436
    14301437                    (let* ((ysn (gensym 'ysn))
    1431                            (retflds (map (lambda (yn yi) `(,yn ,(V:Sub yi (V:Var rv1)) )) yns yis)))
    1432                       (list (B:Val rv2 (V:Rec retflds)))))
     1438                           (xn  (gensym 'xn))
     1439                           (retflds (cons `(,x ,(V:Var xn))
     1440                                          (map (lambda (yn yi) `(,yn ,(V:Sub yi (V:Var rv1)) )) yns yis))))
     1441
     1442                      (list
     1443                       (B:Val ysn (V:Sel 'ysn (V:Var rv1)))
     1444                       (B:Val xn  (V:Sel 'xn  (V:Var rv1)))
     1445                       (B:Val rv2 (V:Rec retflds)))))
    14331446                ))
    14341447             ))
     
    15241537    ))
    15251538
     1539 
     1540(define (prelude/scheme #!key (solver 'rk4b) (random #f) (integral-index 0))
     1541
     1542`(
     1543  ,(case solver
     1544     ((cvode) `("(use sundials random-mtzig mathh)" ,nl))
     1545     (else    `("(use runge-kutta random-mtzig mathh)" ,nl)))
     1546 #<<EOF
     1547
     1548;; Variant types
     1549
     1550(define-syntax define-datatype
     1551  (syntax-rules ()
     1552    [(_ type (name field ...) ...)
     1553     (begin
     1554       (define-constructors type ((name field ...) ...)))]))
     1555
     1556
     1557(define-syntax define-constructors
     1558  (syntax-rules ()
     1559    [(define-constructors type ((name field ...) ...))
     1560     (define-constructors type ((name field ...) ...) (name ...))]
     1561    [(define-constructors type ((name field ...) ...) names)
     1562     (begin
     1563       (define-constructor type (name field ...) names)
     1564       ...)]))
     1565
     1566
     1567(define-syntax define-constructor
     1568  (syntax-rules ()
     1569    [(_ type (name field ...) names)
     1570     (define (name field ...)
     1571       (cons 'type
     1572             (lambda names
     1573               (name field ...))))]))
     1574
     1575
     1576(define-syntax cases
     1577  (syntax-rules ()
     1578    [(_ type x [(name field ...) exp]
     1579          ...)
     1580     ((cdr x) (lambda (field ...) exp)
     1581              ...)]))
     1582
     1583(define-datatype trs (TRSA a) (TRSB b))
     1584(define-datatype trc (TRC f fk e))
     1585
     1586(define (tsCase fa fb x) (cases trs x ((TRSA a) (fa a)) ((TRSB b) (fb b))))
     1587(define (trfOf x)  (cases trc x ((TRC f fk e) f)))
     1588(define (trfkOf x) (cases trc x ((TRC f fk e) fk)))
     1589(define (treOf x)  (cases trc x ((TRC f fk e) e)))
     1590
     1591(define-datatype option (NONE) (SOME a))
     1592(define (swap x v) (cases option v ((NONE) x) ((SOME v) v)))
     1593
     1594(define false #f)
     1595(define true  #t)
     1596
     1597(define equal equal?)
     1598
     1599(define (signalOf v) (if (not v) (error 'signalOf "empty signal" v) v))
     1600
     1601(define (heaviside x) (if (negative? x) 0 1))
     1602
     1603EOF
     1604
     1605,(if (positive? integral-index)
     1606     (case solver
     1607       ((cvode)
     1608        (list "(define integral-solvers (make-vector " integral-index " #f))" nl))
     1609       (else '())
     1610       ) '())
     1611
     1612,(if (not solver)
     1613     `((";; dummy solver; returns only the computed derivatives")
     1614       ("(define (integral f x y h i) (f x y))" ,nl)
     1615       )
     1616     (case solver
     1617       ((cvode)
     1618        `(
     1619          ("(define (integral f x y h i) (f x y))" ,nl)
     1620          ))
     1621
     1622       (else
     1623        `(("(define (scaler x a) (map (lambda (k) (fp* x k)) a))" ,nl)
     1624          ("(define (summer a b) (map fp+ a b))" ,nl)
     1625          ("(define " ,solver " (make-" ,solver "))" ,nl)
     1626          ("(define (make_stepper deriv) (" ,solver " scaler summer deriv))" ,nl)
     1627          ("(define (integral f x y h i) (((make_stepper f) h) x y))" ,nl)
     1628          ))
     1629       ))
     1630))
     1631
     1632
     1633
     1634(define (prelude/ML  #!key (solver 'rk4b) (random #f))
     1635`(
     1636 #<<EOF
     1637structure Model =
     1638struct
     1639
     1640open Real
     1641open Math
     1642open RungeKutta
     1643
     1644datatype ('b,'c) trs = TRSA of 'b | TRSB of 'c
     1645datatype ('a,'b,'c) trc = TRC of ((('a -> (('b,'c) trs))) *
     1646                                  (('a -> (('b,'c) trs))) *
     1647                                  ((('b,'c) trs) -> bool))
     1648         
     1649fun tsCase (fa,fb,x) = case x of TRSA a => (fa a) | TRSB b => (fb b)
     1650fun trfOf x = case x of TRC (f,fk,e) => f
     1651fun trfkOf x = case x of TRC (f,fk,e) => fk
     1652fun treOf x = case x of TRC (f,fk,e) => e
     1653
     1654fun putStrLn str =
     1655  (TextIO.output (TextIO.stdOut, str);
     1656   TextIO.output (TextIO.stdOut, "\n"))
     1657
     1658fun putStr str = (TextIO.output (TextIO.stdOut, str))
     1659
     1660fun showReal n =
     1661let open StringCvt
     1662in
     1663(if n < 0.0 then "-" else "") ^ (fmt (FIX (SOME 12)) (abs n))
     1664end
     1665
     1666fun vmap2 f (v1,v2) =
     1667    let
     1668        val n = Vector.length v1
     1669    in
     1670        Vector.tabulate (n, fn (i) => f (Unsafe.Vector.sub (v1,i),
     1671                                        Unsafe.Vector.sub (v2,i)))
     1672    end
     1673
     1674exception EmptySignal
     1675
     1676val neg = (op ~)
     1677val swap = fn (x,v) => (case v of NONE => x | SOME v => v)
     1678val equal = fn (x,y) => (x = y)
     1679val signalOf = fn (v) => (case v of NONE => raise EmptySignal | SOME v => v)
     1680val heaviside = fn (v) => (if Real.< (v, 0.0) then 0.0 else 1.0)
     1681EOF
     1682
     1683,(if random
     1684#<<EOF
     1685
     1686fun RandomInit () = RandomMTZig.fromEntropy()
     1687
     1688val RandomState = RandomInit ()
     1689
     1690fun random_uniform () = RandomMTZig.randUniform RandomState
     1691
     1692
     1693fun PoissonInit () =
     1694 let
     1695     val zt = RandomMTZig.ztnew()
     1696     val st = RandomMTZig.fromEntropy()
     1697 in
     1698     {st=st,zt=zt}
     1699 end
     1700
     1701fun PoissonStep (rate,t,h,st,zt) =
     1702 let
     1703    val rv     = RandomMTZig.randPoisson (rate*0.001*h,st,zt)
     1704    val spike' = Real.> (rv,0.0)
     1705    val spikeCount' = if spike' then rv else 0.0
     1706  in
     1707     {t=t+h,spike=spike',spikeCount=spikeCount',st=st,zt=zt}
     1708  end
     1709
     1710
     1711EOF
     1712"")
     1713
     1714,(if (not solver)
     1715     `(("(* dummy solver; returns only the computed derivatives *)")
     1716       ("fun integral (f,x: real,y: real vector,h,i) = (f (x,y))" ,nl)
     1717       )
     1718     `(("val summer = fn (a,b) => (vmap2 (fn (x,y) => x+y) (a,b))" ,nl)
     1719       ("val scaler = fn(a,lst) => (Vector.map (fn (x) => a*x) lst)" ,nl)
     1720       . ,(case solver 
     1721            ;; adaptive solvers
     1722            ((rkhe rkbs rkf45 rkck rkdp rkf78 rkv65)
     1723             `(
     1724               ("val " ,solver ": (real vector) stepper2 = make_" ,solver "()" ,nl)
     1725               ("fun make_stepper (deriv) = " ,solver " (scaler,summer,deriv)" ,nl)
     1726               ("val cerkdp: (real vector) stepper3 = make_cerkdp()" ,nl)
     1727               ("fun make_estepper (deriv) = cerkdp (scaler,summer,deriv)" ,nl)
     1728#<<EOF
     1729
     1730val tol = Real.Math.pow (10.0, ~7.0)
     1731val lb = 0.5 * tol
     1732val ub = 0.9 * tol
     1733
     1734datatype ('a, 'b) either = Left of 'a | Right of 'b
     1735
     1736
     1737fun predictor tol (h,ys) =
     1738  let open Real
     1739      val e = Vector.foldl (fn (y,ax) => Real.+ ((abs y),ax)) 0.0 ys
     1740  in
     1741      if e < lb
     1742      then Right (1.414*h)      (* step too small, accept but grow *)
     1743      else (if e < ub
     1744            then Right h        (* step just right *)
     1745            else Left (0.5*h))  (* step too large, reject and shrink *)
     1746  end
     1747
     1748
     1749exception ConvergenceError
     1750
     1751
     1752fun secant tol f fg0 guess1 guess0 =
     1753    let open Real
     1754        val fg1 = f guess1
     1755        val newGuess = guess1 - fg1 * (guess1 - guess0) / (fg1 - fg0)
     1756        val err =  abs (newGuess - guess1)
     1757    in
     1758        if (err < tol)
     1759        then newGuess
     1760        else secant tol f fg1 newGuess guess1
     1761    end
     1762
     1763
     1764datatype 'a result = Next of 'a | Root of 'a
     1765
     1766
     1767fun esolver (stepper,evtest) (x,ys,h) =
     1768    let open Real
     1769        val (ys',e,finterp) = stepper h (x,ys)
     1770    in
     1771        case predictor tol (h,e) of
     1772            Right h' =>
     1773            if (evtest (ys') >= 0.0)
     1774            then (let
     1775                     val theta   = secant tol (evtest o finterp) (evtest ys) 1.0 0.0
     1776                     val ys''    = finterp (theta+tol)
     1777                 in
     1778                     Root (x+(theta+tol)*h,ys'',h')
     1779                 end)
     1780            else Next (x+h,ys',h')
     1781          | Left h'  =>
     1782            esolver (stepper,evtest) (x,ys,h')
     1783    end
     1784
     1785
     1786fun eintegral (f,x,ys,evtest,h,i) =
     1787    case esolver (make_estepper f,evtest) (x,ys,h) of
     1788        Next (xn,ysn,h') =>
     1789        ({xn=xn,h=h',ysn=ysn})
     1790      | Root (xn,ysn,h') =>
     1791        ({xn=xn,ysn=ysn,h=h'})
     1792
     1793fun solver stepper (x,ys,h) =
     1794    let open Real
     1795        val (ys',e) = stepper h (x,ys)
     1796    in
     1797        case predictor tol (h,e) of
     1798            Right h' =>
     1799            (x+h,ys',h')
     1800          | Left h'  =>
     1801            solver (stepper) (x,ys,h')
     1802    end
     1803
     1804fun integral (f,x,ys,evtest,h,i) =
     1805    let
     1806        val (xn,ysn,h') = solver (make_stepper f) (x,ys,h)
     1807    in
     1808        {xn=xn,ysn=ysn,h=h'}
     1809    end
     1810
     1811
     1812
     1813EOF
     1814
     1815               ))
     1816            (else
     1817             `(
     1818               ("val " ,solver ": (real vector) stepper1 = make_" ,solver "()" ,nl)
     1819               ("fun make_stepper (deriv) = " ,solver " (scaler,summer,deriv)" ,nl)
     1820               ("fun integral (f,x: real,y: real vector,h,i) = ((make_stepper f) h) (x,y)" ,nl)
     1821               ))
     1822            ))
     1823       )
     1824))
     1825
     1826
     1827 
     1828(define (prelude/Octave #!key (solver 'lsode))
     1829`(
     1830#<<EOF
     1831
     1832function res = TRSA(v)
     1833  res = struct ("TRSA",v); 
     1834end
     1835
     1836function res = TRSB(v)
     1837  res = struct ("TRSB",v); 
     1838end
     1839
     1840function res = tsCase(fa,fb,v)
     1841   if (isfield (v, "TRSA"))
     1842     res = fa(getfield(v,"TRSA"));
     1843   else
     1844     res = fb(getfield(v,"TRSB"));
     1845   endif
     1846end
     1847
     1848function res = TRSC(f,fk,e)
     1849  res = struct ("TRSC",[f,fk,e]); 
     1850end
     1851
     1852function res = trfOf(x)
     1853   res = getfield(x,"TRC")(1);
     1854end
     1855
     1856function res = trfkOf(x)
     1857   res = getfield(x,"TRC")(2);
     1858end
     1859
     1860function res = treOf(x)
     1861   res = getfield(x,"TRC")(3);
     1862end
     1863
     1864function res = NONE()
     1865  res = struct (); 
     1866end
     1867
     1868function res = SOME(v)
     1869  res = struct ("SOME",v); 
     1870end
     1871
     1872function res = swap (x,v)
     1873   if (isfield (v, "SOME"))
     1874     res = getfield(v,"SOME");
     1875   else
     1876     res = x;
     1877   endif
     1878end
     1879
     1880function res = equal (x, y)
     1881   res = (x==y);
     1882end
     1883
     1884function res = signalOf (v)
     1885  if (not(v))
     1886   error ("empty signal")
     1887 else
     1888   res = v;
     1889 endif
     1890end
     1891
     1892function res = neg (x)
     1893  res = -x;
     1894end
     1895
     1896function res = ifv (b,x,y)
     1897  if (b)
     1898    res = x;
     1899  else
     1900    res = y;
     1901  endif
     1902end
     1903
     1904EOF
     1905
     1906))
    15261907
    15271908(define (codegen/Octave name f #!key (initial #f) (pre #t) (solver #f))
  • release/4/signal-diagram/trunk/signal-diagram.setup

    r30957 r30972  
    8282    ))
    8383
     84
     85(make (
     86       ((dynld-name "runge-kutta") ("runge-kutta.scm" )
     87        (compile -O2 -S -s runge-kutta.scm -j runge-kutta))
     88
     89       ((dynld-name "runge-kutta.import") ("runge-kutta.import.scm")
     90        (compile -O2 -S -s runge-kutta.import.scm))
     91       )
     92
     93  (list (dynld-name "runge-kutta")
     94        (dynld-name "runge-kutta.import"))
     95  )
     96
     97(install-extension
     98
     99  ; Name of your extension:
     100  'runge-kutta
     101
     102  ; Files to install for your extension:
     103  `(,(dynld-name "runge-kutta") ,(dynld-name "runge-kutta.import") )
     104
     105  ; Assoc list with properties for your extension:
     106  `((version ,version)
     107    ))
     108
     109
     110;; From setup-header.scm by Kon Lovett
     111
     112(define (installation-chicken-home)
     113  (if (not (installation-prefix)) (chicken-home)
     114    (make-pathname `(,(installation-prefix) "share") "chicken") ) )
     115
     116;;; Constants & Procedures
     117
     118(define SHARED-DIR (installation-chicken-home))
     119(define SIGNAL-DIAGRAM-DIR (make-pathname SHARED-DIR "signal-diagram"))
     120
     121;; File Copy Operations
     122
     123(define (copy-file-to-signal-diagram-dir fn)
     124  (let ([fn (->string fn)])
     125    (copy-file fn (make-pathname SIGNAL-DIAGRAM-DIR fn)) ) )
     126
     127(copy-file-to-signal-diagram-dir "sml-lib/rk/rk.sml")
     128(copy-file-to-signal-diagram-dir "sml-lib/rk/rk.mlb")
     129
     130(copy-file-to-signal-diagram-dir "sml-lib/randmtzig/randmtziglib.c")
     131(copy-file-to-signal-diagram-dir "sml-lib/randmtzig/randmtzig.sml")
     132(copy-file-to-signal-diagram-dir "sml-lib/randmtzig/randmtzig.mlb")
     133
     134(copy-file-to-signal-diagram-dir "sml-lib/tensor/DynArray.sml")
     135(copy-file-to-signal-diagram-dir "sml-lib/tensor/tensor.sml")
     136(copy-file-to-signal-diagram-dir "sml-lib/tensor/sparse.sml")
     137(copy-file-to-signal-diagram-dir "sml-lib/tensor/tensor.mlb")
     138(copy-file-to-signal-diagram-dir "sml-lib/tensor/sparse.mlb")
     139
     140
     141(copy-file-to-signal-diagram-dir "octave-lib/ode15s.m")
Note: See TracChangeset for help on using the changeset viewer.