init function for QSS

Author: BenLauwens

src/SimJulia.jl

@@ -11,6 +11,7 @@ module SimJulia
 
   import Base.run, Base.now, Base.isless, Base.show, Base.interrupt, Base.yield
   import Base.(&), Base.(|)
+  import TaylorSeries.integrate
 
   export AbstractEvent
   export value, state, environment
@@ -44,6 +45,9 @@ module SimJulia
   include("resources/base.jl")
   include("resources/containers.jl")
   include("resources/stores.jl")
+  include("odes/base.jl")
   include("odes/macros.jl")
+  include("odes/QSS.jl")
   include("continuous.jl")
+  include("odes/commons.jl")
 end

src/continuous.jl

@@ -1,31 +1,21 @@
+mutable struct Variable <: AbstractEvent
+  bev :: BaseEvent
+  id :: UInt
+  f :: Function
+  x :: Taylor1
+  t :: Float64
+  function Variable(env::Environment, id::Int, f::Vector{Function}, x₀::Taylor1)
+    new(BaseEvent(env), UInt(id), f[id], x₀, now(env))
+  end
+end
+
 struct Continuous <: AbstractProcess
   bev :: BaseEvent
-  f :: Vector{Function}
-  order :: UInt8
-  stiff :: Bool
-  t :: Vector{Float64}
-  q :: Vector{Taylor1}
-  x :: Vector{Taylor1}
+  vars :: Vector{Variable}
   p :: Vector{Float64}
-  function Continuous(func::Function, env::Environment, q::Vector{Float64}, p::Vector{Float64}=Float64[]; order::Number=4, stiff::Bool=false)
-    cont = new(BaseEvent(env), func(), UInt8(order), stiff, Vector{Float64}(), Vector{Taylor1}(), Vector{Taylor1}(), p)
-    zero_taylor = 0*Taylor1(Float64, order+1)
-    for q₀ in q
-      push!(cont.t, now(env))
-      push!(cont.q, q₀ + zero_taylor)
-    end
-    t₀ = now(env) + Taylor1(Float64, order+1)
-    for (i, q₀) in enumerate(q)
-      push!(cont.x, integrate(cont.f[i](t₀, cont.q, cont.p), q₀))
-    end
-    for i in 1:order-1
-      for (i, xᵢ) in enumerate(cont.x)
-        cont.q[i] = copy(cont.x[i])
-      end
-      for (i, qᵢ) in enumerate(cont.q)
-        cont.x[i] = integrate(cont.f[i](t₀, cont.q, cont.p), q[i])
-      end
-    end
+  function Continuous(model::Function, env::Environment, x₀::Vector{Float64}, p::Vector{Float64}=Float64[]; integrator::Integrator=QSS())
+    cont = new(BaseEvent(env), Vector{Variable}(), p)
+    init(env, cont, integrator, model, x₀)
     cont
   end
 end

src/odes/QSS.jl

@@ -0,0 +1,16 @@
+abstract type Quantizer end
+
+struct DirectQuantizer <: Quantizer end
+
+struct ImplicitQuantizer <: Quantizer end
+
+struct QSS <: Integrator
+  order :: UInt8
+  quantizer :: Quantizer
+  q :: Vector{Taylor1}
+  t :: Vector{Float64}
+  function QSS(;order::Number=4, stiff::Bool=false)
+    quantizer = stiff ? ImplicitQuantizer() : DirectQuantizer()
+    new(UInt8(order), quantizer, Vector{Taylor1}(), Vector{Float64}())
+  end
+end

src/odes/base.jl

@@ -0,0 +1 @@
+abstract type Integrator end

src/odes/commons.jl

@@ -0,0 +1,36 @@
+function init(env::Environment, cont::Continuous, integrator::QSS, model::Function, x₀::Vector{Float64})
+  n = length(x₀)
+  t = now(env)
+  f, deps = model()
+  zero_taylor = 0*Taylor1(Float64, integrator.order+1)
+  for (i, q₀) in enumerate(x₀)
+    push!(integrator.t, t)
+    push!(integrator.q, q₀ + zero_taylor)
+  end
+  t₀ = t + Taylor1(Float64, integrator.order+1)
+  x_vec = Vector{Taylor1}()
+  for (i, q₀) in enumerate(x₀)
+    push!(x_vec, integrate(f[i](t₀, integrator.q, cont.p), q₀))
+  end
+  for i in 1:integrator.order-1
+    for (j, x) in enumerate(x_vec)
+      integrator.q[j] = copy(x)
+    end
+    for (j, q₀) in enumerate(x₀)
+      x_vec[j] = integrate(f[j](t₀, integrator.q, cont.p), q₀)
+    end
+  end
+  for (i, x) in enumerate(x_vec)
+    var = Variable(env, i, f, x)
+    push!(cont.vars, var)
+    @callback step(var, cont, integrator)
+    schedule(var)
+  end
+end
+
+function step(var::Variable, cont::Continuous, integrator::QSS)
+  i = var.id
+  env = environment(var)
+  t = now(env)
+  println(integrator.q)
+end

src/odes/macros.jl

@@ -1,3 +1,15 @@
+function check_dependencies(ex, deps::Array{Bool}, i::Int)
+  if ex isa Expr
+    if ex.head == :ref
+      deps[i, ex.args[2]] = true
+    else
+      for arg in ex.args
+        check_dependencies(arg, deps, i)
+      end
+    end
+  end
+end
+
 macro model(expr::Expr)
   expr.head != :function && error("Expression is not a function definition!")
   args = getArguments(expr)
@@ -6,9 +18,16 @@ macro model(expr::Expr)
   for ex in expr.args[2].args
     ex isa Expr && ex.head == Symbol("=") && push!(f_vec, ex)
   end
+  n = length(f_vec)
+  deps = zeros(Bool, n, n)
+  for ex in expr.args[2].args
+    if ex isa Expr && ex.head == Symbol("=")
+      check_dependencies(ex.args[2], deps, ex.args[1].args[2])
+    end
+  end
   esc(:(function $func_name()
       f = Array{Function}(length($f_vec))
       $((:(f[$(f_vec[i].args[1].args[2])] = (t::TaylorSeries.Taylor1, q::Vector{TaylorSeries.Taylor1}, p::Vector{Float64})->$(f_vec[i].args[2])) for i in 1:length(:($f_vec)))...)
-      f
+      f, $deps
     end))
 end

test/continuous.jl

@@ -6,7 +6,8 @@ using SimJulia
 end
 
 sim = Simulation()
-cont = @continuous stiffeq(sim, [0.0, 20.0], [2020.0]; order=6, stiff=true)
+cont = @continuous stiffeq(sim, [0.0, 20.0], [2020.0])
 run(sim, 500.0)
-println(cont.x)
-println(cont.q)
+for var in cont.vars
+  println(var.x)
+end