Update to cope with changes from TalorSeries v0.5

Author: BenLauwens

src/continuous.jl

@@ -3,15 +3,15 @@ abstract type ContinuousProcess <: AbstractProcess end
 mutable struct Variable <: AbstractEvent
   bev :: BaseEvent
   id :: UInt
-  x :: Taylor1
+  x :: Taylor1{Float64}
   t :: Float64
-  function Variable(env::Environment, id::Int, x::Taylor1, t::Float64)
+  function Variable(env::Environment, id::Int, x::Taylor1{Float64}, t::Float64)
     new(BaseEvent(env), id, x, t)
   end
 end
 
 function advance_time(var::Variable, t::Float64=now(environment(var)))
-  var.x = evaluate(var.x, t - var.t + Taylor1(var.x.order))
+  var.x = evaluate(var.x, t - var.t + Taylor1(Float64, var.x.order))
   var.t = t
   var.x[1]
 end
@@ -67,9 +67,3 @@ macro continuous(expr::Expr)
   end
   esc(:(Continuous($func(), $(args...); $(params...))))
 end
-
-function schedule(var::Variable, cont::Continuous, integrator::Integrator, Δt::Float64=0.0)
-  var.bev = BaseEvent(environment(var))
-  @callback step(var, cont, integrator)
-  schedule(var, Δt)
-end

src/odes/QSS.jl

@@ -17,51 +17,52 @@ struct QSS{T} <: Integrator
       push!(qss.q, q₀ + Taylor1(zeros(Float64, order+1)))
     end
     for i in 1:order-1
+      q = deepcopy(qss.q)
       for (j, q₀) in enumerate(x₀)
-        qss.q[j] = integrate(model.f[j](t₀, qss.q, p), q₀)
+        qss.q[j] = integrate(model.f[j](t₀, q, p), q₀)
       end
     end
-    println(qss.q)
     qss
   end
 end
 
 function Continuous(model::Model, env::Environment, x₀::Vector{Float64}, p::Vector{Float64}=Float64[];
-                    stiff::Bool=false, order::Number=4, Δrel::Float64=1e-6, Δabs::Float64=1e-6)
-  stiff ? Continuous(model, env, QSS{SimJulia.stiff}, x₀, p; order=order, Δrel=Δrel, Δabs=Δabs) :
-          Continuous(model, env, QSS{SimJulia.non_stiff}, x₀, p; order=order, Δrel=Δrel, Δabs=Δabs)
+                    stiff::Bool=false, order::Number=4, Δrel::Number=1e-6, Δabs::Number=1e-6)
+  stiff ? Continuous(model, env, QSS{SimJulia.stiff}, x₀, p; order=order, Δrel=float(Δrel), Δabs=float(Δabs)) :
+          Continuous(model, env, QSS{SimJulia.non_stiff}, x₀, p; order=order, Δrel=float(Δrel), Δabs=float(Δabs))
 end
 
 function initial_values(qss::QSS, t::Float64)
-  t₀ = t + Taylor1(Float64, qss.order+0)
+  t₀ = t + Taylor1(Float64, qss.order+1)
   x₀ = Vector{Taylor1{Float64}}()
   for (i, f) in enumerate(qss.model.f)
     push!(x₀, integrate(f(t₀, qss.q, qss.p), qss.q[i][1]))
   end
-  println(x₀)
   x₀
 end
 
 function step(var::Variable, cont::Continuous, qss::QSS)
   t = now(environment(var))
   n = length(qss.model.f)
   i = var.id
-  t₀ = t + Taylor1(Float64, qss.order+0)
+  t₀ = t + Taylor1(Float64, qss.order+1)
   x₀ = advance_time(var, t)
-  update_quantized_state(qss, cont.vars, i, t)
-  println(x₀)
+  update_quantized_state(qss, var, t)
   Δt = compute_next_time(var.x, max(qss.Δrel*x₀, qss.Δabs))
-  schedule(var, cont, qss, Δt)
+  reset(var)
+  schedule(var, Δt)
   for j in filter(j->qss.model.deps[j,i], 1:n)
     dep = cont.vars[j]
     x₀ = evaluate(dep.x, t - dep.t)
     dep.t = t
-    for k in filter(k->qss.model.deps[j,k], 1:n)
+    advance_time(qss, j, t)
+    for k in filter(k->qss.model.deps[j,k] && k!=j, 1:n)
       advance_time(qss, k, t)
     end
     dep.x = integrate(qss.model.f[j](t₀, qss.q, qss.p), x₀)
     Δt = recompute_next_time(qss, dep.x, qss.q[j], max(qss.Δrel*x₀, qss.Δabs))
-    schedule(dep, cont, qss, Δt)
+    reset(dep)
+    schedule(dep, Δt)
   end
 end
 
@@ -71,8 +72,9 @@ function advance_time(qss::QSS, i::Int, t::Float64)
   qss.q[i][1]
 end
 
-function update_quantized_state(qss::QSS{non_stiff}, vars::Vector{Variable}, i::UInt, t::Float64)
-  qss.q[i] = copy(vars[i].x)
+function update_quantized_state(qss::QSS{non_stiff}, var::Variable, t::Float64)
+  i = var.id
+  qss.q[i] = deepcopy(var.x)
   qss.q[i][end] = 0.0
   qss.t[i] = t
 end
@@ -88,7 +90,7 @@ function compute_next_time(x::Taylor1, Δq::Float64)
   (abs(Δq/x[end]))^(1.0/x.order)
 end
 
-function recompute_next_time(::QSS{non_stiff}, x::Taylor1, q::Taylor1, Δq::Float64)
+function recompute_next_time(::QSS{non_stiff}, x::Taylor1{Float64}, q::Taylor1{Float64}, Δq::Float64)
   p = (x-q).coeffs
   p[1] -= Δq
   neg = roots(p)

test/continuous.jl

@@ -1,21 +1,33 @@
 using SimJulia
 
+@model function simple(t, x, p, dx)
+  dx[1] = x[1]
+end
+
+#sim = Simulation()
+#cont = @continuous simple(sim, [1.0]; stiff=false, order=4, Δrel=1e-16, Δabs=1e-6)
+#run(sim, 10)
+
 @model function diffeq(t, x, p, dx)
   dx[1] = p[2]+0.01*x[2]
   dx[2] = p[1]-100.0*x[1]-100.0*x[2]
 end
 
 function report(sim::Simulation, cont::Continuous)
+  λ, P = eig([0.0 0.01; -100.0 -100.0])
   while true
-    println(now(sim), " ", evaluate(cont.vars[1]), " ", evaluate(cont.vars[2]))
+    t = now(sim)
+    xs = [evaluate(cont.vars[1]), evaluate(cont.vars[2])]
+    xe = P*diagm(exp.(λ*t))*inv(P)*[0.0; 20.0]+P*diagm((exp.(λ*t)-1)./λ)*inv(P)*[0.0;2020.0]
+    println(t, " ", xs[1], " ", abs(xe[1]-xs[1]), " ", xs[2], " ", abs(xe[2]-xs[2]))
     yield(Timeout(sim, 1.0))
   end
 end
 
 sim = Simulation()
-cont = @continuous diffeq(sim, [0.0, 20.0], [2020.0, 0.0]; stiff=false, order=4)
+cont = @continuous diffeq(sim, [0.0, 20.0], [2020.0, 0.0]; stiff=false, order=4, Δrel=1e-16, Δabs=1e-6)
 @process report(sim, cont)
-run(sim, 71)
+@time run(sim, 500)
 
 @model function bouncing_ball(t, x, p, dx)
   g = 9.81
@@ -32,7 +44,7 @@ run(sim, 71)
   dx[2] = -g - p[1]*f/m
 end
 
-sim = Simulation()
-cont = @continuous bouncing_ball(sim, [2.0, 0.0], [0.0])
-@process report(sim, cont)
-run(sim, 5)
+#sim = Simulation()
+#cont = @continuous bouncing_ball(sim, [2.0, 0.0], [0.0])
+#@process report(sim, cont)
+#run(sim, 5)