Some more commits

Author: BenLauwens

src/continuous.jl

@@ -14,42 +14,41 @@ function advance_time(var::Variable, t::Float64)
   var.x.coeffs[1]
 end
 
-function evaluate(var::Variable, t::Float64)
+function evaluate(var::Variable, t::Float64=now(environment(var)))
   evaluate(var.x, t - var.t)
 end
 
 struct ZeroCrossing <: AbstractEvent
 
 end
 
-macro zerocrossing(expr::Expr)
+macro trigger(expr::Expr)
   expr.head != :call && error("Expression is not a function call!")
   nothing
 end
 
 struct Continuous <: AbstractProcess
   bev :: BaseEvent
-  integrator :: Integrator
   vars :: Vector{Variable}
-  function Continuous(env::Environment, integrator::Integrator)
-    cont = new(BaseEvent(env), integrator, Vector{Variable}())
-    t = now(env)
-    x = initial_values(integrator, t)
-    for (i, x₀) in enumerate(x)
-      push!(cont.vars, Variable(env, i, x₀, t))
-    end
-    for var in cont.vars
-      @callback step(var, cont, integrator)
-      schedule(var)
-    end
-    cont
+  function Continuous(env::Environment)
+    new(BaseEvent(env), Vector{Variable}())
   end
 end
 
 function Continuous{I<:Integrator}(model::Model, env::Environment, ::Type{I},
     x₀::Vector{Float64}, p::Vector{Float64}=Float64[]; args...)
+  cont = Continuous(env)
+  t = now(env)
   integrator = I(model, now(env), x₀, p; args...)
-  Continuous(env, integrator)
+  x = initial_values(integrator, t)
+  for (i, x₀) in enumerate(x)
+    push!(cont.vars, Variable(env, i, x₀, t))
+  end
+  for var in cont.vars
+    @callback step(var, cont, integrator)
+    schedule(var)
+  end
+  cont
 end
 
 macro continuous(expr::Expr)

src/odes/QSS.jl

@@ -11,10 +11,10 @@ struct QSS{T} <: Integrator
   function QSS{T}(model::Model, t::Float64, x₀::Vector{Float64}, p::Vector{Float64};
                   order::Number=4, Δrel::Float64=1e-6, Δabs::Float64=1e-6) where T
     qss = new(UInt8(order), model, p, Vector{Taylor1}(), Vector{Float64}(), Δrel, Δabs)
-    t₀ = t + Taylor1(Float64, qss.order+1)
+    t₀ = t + Taylor1(Float64, order)
     for q₀ in x₀
       push!(qss.t, t)
-      push!(qss.q, q₀ + Taylor1(zeros(Float64, qss.order+1)))
+      push!(qss.q, q₀ + Taylor1(zeros(Float64, order+1)))
     end
     for i in 1:qss.order-1
       for (j, q₀) in enumerate(x₀)
@@ -32,7 +32,7 @@ function Continuous(model::Model, env::Environment, x₀::Vector{Float64}, p::Ve
 end
 
 function initial_values(qss::QSS, t::Float64)
-  t₀ = t + Taylor1(Float64, qss.order+1)
+  t₀ = t + Taylor1(Float64, qss.order+0)
   x₀ = Vector{Taylor1}()
   for (i, f) in enumerate(qss.model.f)
     push!(x₀, integrate(f(t₀, qss.q, qss.p), qss.q[i].coeffs[1]))
@@ -44,7 +44,7 @@ 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+1)
+  t₀ = t + Taylor1(Float64, qss.order+0)
   x₀ = advance_time(var, t)
   update_quantized_state(qss, cont.vars, i, t)
   Δt = compute_next_time(var.x, max(qss.Δrel*x₀, qss.Δabs))
@@ -63,7 +63,7 @@ function step(var::Variable, cont::Continuous, qss::QSS)
 end
 
 function advance_time(qss::QSS, i::Int, t::Float64)
-  qss.q[i] = evaluate(qss.q[i], t - qss.t[i] + Taylor1(qss.q[i].order))
+  qss.q[i] = evaluate(qss.q[i], t - qss.t[i] + Taylor1(qss.order+0))
   qss.t[i] = t
   qss.q[i].coeffs[1]
 end
@@ -86,7 +86,7 @@ function compute_next_time(x::Taylor1, Δq::Float64)
 end
 
 function recompute_next_time(::QSS{non_stiff}, x::Taylor1, q::Taylor1, Δq::Float64)
-  p = x.coeffs - q.coeffs
+  p = (x-q).coeffs
   p[1] -= Δq
   neg = roots(p)
   p[1] += 2Δq

src/odes/macros.jl

@@ -1,6 +1,6 @@
 function check_dependencies(ex, deps::Array{Bool}, i::Int, x::Symbol)
   if ex isa Expr
-    if ex.head == :ref && x == ex.args[1]
+    if ex.head == :ref && ex.args[1] == x
       deps[i, ex.args[2]] = true
     else
       for arg in ex.args
@@ -15,13 +15,14 @@ macro model(expr::Expr)
   args = getArguments(expr)
   func_name = shift!(args)
   f_vec = Vector{Expr}()
+  dx = expr.args[1].args[5]
   for ex in expr.args[2].args
-    ex isa Expr && ex.head == Symbol("=") && push!(f_vec, ex)
+    ex isa Expr && ex.head == Symbol("=") && ex.args[1] isa Expr && ex.args[1].head == :ref && ex.args[1].args[1] == dx && push!(f_vec, ex)
   end
   n = length(f_vec)
   deps = zeros(Bool, n, n)
   for ex in expr.args[2].args
-    ex isa Expr && ex.head == Symbol("=") && check_dependencies(ex.args[2], deps, ex.args[1].args[2], expr.args[1].args[3])
+    ex isa Expr && ex.head == Symbol("=") && ex.args[1] isa Expr && ex.args[1].head == :ref && ex.args[1].args[1] == dx && check_dependencies(ex.args[2], deps, ex.args[1].args[2], expr.args[1].args[3])
   end
   esc(:(function $func_name()
       f = Array{Function}($n)
@@ -30,7 +31,7 @@ macro model(expr::Expr)
     end))
 end
 
-macro trigger(expr::Expr)
+macro zerocrossing(expr::Expr)
   expr.head != :function && error("Expression is not a function definition!")
   nothing
 end

test/continuous.jl

@@ -1,27 +1,46 @@
 using SimJulia
 
-@model function diffeq(t, x, p, dq)
+@model function diffeq(t, x, p, dx)
   dx[1] = 0.01*x[2]
   dx[2] = p[1]-100.0*x[1]-100.0*x[2]
 end
 
-@trigger function less_prey(t, x, p, res)
-  res = x[1] < x[2]
-  if res
-    p[1] = 0.0
-  end
-end
-
 function report(sim::Simulation, cont::Continuous)
   while true
-    t = now(sim)
-    println(t, " ", evaluate(cont.vars[1], t), " ", evaluate(cont.vars[2], t))
+    println(now(sim), " ", evaluate(cont.vars[1]), " ", evaluate(cont.vars[2]))
     yield(Timeout(sim, 1.0))
   end
 end
 
 sim = Simulation()
-cont = @continuous diffeq(sim, [0.0, 20.0], [2020.0]; stiff=false, order=5)
+cont = @continuous diffeq(sim, [0.0, 20.0], [2020.0]; stiff=false, order=4)
 @process report(sim, cont)
-zc = @zerocrossing less_prey(cont)
 run(sim, 71)
+
+@model function bouncing_ball(t, x, p, dx)
+  dx[1] = x[3]
+  dx[3] = zero(t)
+  dx[2] = x[4]
+  dx[4] = -9.81 + p[1]*(2.0e4/0.1*(x[2]-(3.0-ceil(evaluate(x[1])/0.2)*0.2)))
+end
+
+@zerocrossing function contact(t, x, p, res)
+  yf = 3.0 - ceil(x[1]/0.2)*0.2
+  res = x[4] < 0.0 && x[2] - yf < 0.05
+  if res
+    p[1] = 1.0
+  end
+end
+
+@zerocrossing function up(t, x, p, res)
+  yf = 3.0 - ceil(x[1]/0.2)*0.2
+  res = x[4] > 0.0 && x[2] - yf > 0.05
+  if res
+    p[1] = 0.0
+  end
+end
+
+sim = Simulation()
+cont = @continuous bouncing_ball(sim, [0.0, 3.05, 1.0, 0.0], [0.0])
+@process report(sim, cont)
+run(sim, 5)