From 81c44a150303570bbd1baf4c698a715043cffecd Mon Sep 17 00:00:00 2001 From: slwu89 <10673535+slwu89@users.noreply.github.com> Date: Wed, 3 Aug 2022 15:50:54 -0700 Subject: [PATCH 1/6] (not working) attempt to add vector wires for instantaneous continuous machines --- src/dwd_dynam.jl | 60 ++++++++++++++++++++++++++++++++++++++++------- test/dwd_dynam.jl | 16 +++++++++++++ 2 files changed, 67 insertions(+), 9 deletions(-) diff --git a/src/dwd_dynam.jl b/src/dwd_dynam.jl index aeef108..562baa0 100644 --- a/src/dwd_dynam.jl +++ b/src/dwd_dynam.jl @@ -70,16 +70,49 @@ dependency_pairs(interface::InstantaneousDirectedInterface) = map(apex(dependenc legs(dependency(interface))[2](i) => legs(dependency(interface))[1](i) end |> sort +# InstantaneousDirectedVectorInterface +struct InstantaneousDirectedVectorInterface{T,N} <: AbstractDirectedInterface{T} + input_ports::Vector + output_ports::Vector + dependency::Span # P_in <- R -> P_out +end + +InstantaneousDirectedVectorInterface{T,N}(input_ports::AbstractVector, output_ports::AbstractVector, dependency_pairs::AbstractVector) where {T,N} = +InstantaneousDirectedVectorInterface{T,N}(input_ports, output_ports, + Span( + FinFunction(Array{Int}(last.(dependency_pairs)), length(dependency_pairs), length(input_ports)), + FinFunction(Array{Int}(first.(dependency_pairs)), length(dependency_pairs), length(output_ports)) + )) + +InstantaneousDirectedVectorInterface{T,N}(ninputs::Int, noutputs::Int, dependency) where {T,N} = + InstantaneousDirectedVectorInterface{T,N}(1:ninputs, 1:noutputs, dependency) + +InstantaneousDirectedVectorInterface{T,N}(input_ports::AbstractVector, output_ports::AbstractVector, ::Nothing) where {T,N} = + InstantaneousDirectedVectorInterface{T,N}(input_ports, output_ports, vcat( + map(1:length(input_ports)) do i + map(1:length(output_ports)) do j + j => i + end + end...) +) + +dependency(interface::InstantaneousDirectedVectorInterface) = interface.dependency +dependency_pairs(interface::InstantaneousDirectedVectorInterface) = map(apex(dependency(interface))) do i + legs(dependency(interface))[2](i) => legs(dependency(interface))[1](i) +end |> sort + +# methods for directed interfaces input_ports(interface::AbstractDirectedInterface) = interface.input_ports output_ports(interface::AbstractDirectedInterface) = interface.output_ports ninputs(interface::AbstractDirectedInterface) = length(input_ports(interface)) noutputs(interface::AbstractDirectedInterface) = length(output_ports(interface)) ndims(::DirectedVectorInterface{T, N}) where {T,N} = N - +ndims(::InstantaneousDirectedVectorInterface{T, N}) where {T,N} = N zero(::Type{I}) where {T, I<:AbstractDirectedInterface{T}} = zero(T) zero(::Type{DirectedVectorInterface{T,N}}) where {T,N} = zeros(T,N) +zero(::Type{InstantaneousDirectedVectorInterface{T,N}}) where {T,N} = zeros(T,N) ### Dynamics abstract type AbstractDirectedSystem{T} end @@ -144,7 +177,7 @@ The readout function `r` may depend on the state, parameters and time, so it mus If it is left out, then ``r=u``. """ const ContinuousMachine{T,I} = Machine{T, I, ContinuousDirectedSystem{T}} -const InstantaneousContinuousMachine{T} = Machine{T, InstantaneousDirectedInterface{T}, ContinuousDirectedSystem{T}} +const InstantaneousContinuousMachine{T,I} = Machine{T, I, ContinuousDirectedSystem{T}} ContinuousMachine{T}(interface::I, system::ContinuousDirectedSystem{T}) where {T, I <: AbstractDirectedInterface} = ContinuousMachine{T, I}(interface, system) @@ -164,13 +197,16 @@ ContinuousMachine{T,I}(ninputs::Int, nstates::Int, dynamics) where {T,I} = ContinuousMachine{T,I}(ninputs, nstates, noutputs, dynamics, readout, dependency) where {T, I<:InstantaneousDirectedInterface{T}} = ContinuousMachine{T, I}(I(ninputs, noutputs, dependency), ContinuousDirectedSystem{T}(nstates, dynamics, readout)) -InstantaneousContinuousMachine{T}(ninputs, nstates, noutputs, dynamics, readout, dependency) where T = - ContinuousMachine{T}(InstantaneousDirectedInterface{T}(ninputs, noutputs, dependency), ContinuousDirectedSystem{T}(nstates, dynamics, readout)) +InstantaneousContinuousMachine{T}(ninputs, nstates, noutputs, dynamics, readout, dependency) where {T, I <: AbstractDirectedInterface} = + ContinuousMachine{T, I}(InstantaneousDirectedInterface{T}(ninputs, noutputs, dependency), ContinuousDirectedSystem{T}(nstates, dynamics, readout)) + +InstantaneousContinuousMachine{T,N}(ninputs, nstates, noutputs, dynamics, readout, dependency) where {T,N} = + ContinuousMachine{T, InstantaneousDirectedVectorInterface{T,N}}(InstantaneousDirectedVectorInterface{T,N}(ninputs, noutputs, dependency), ContinuousDirectedSystem{T}(nstates, dynamics, readout)) -InstantaneousContinuousMachine{T}(f::Function, ninputs::Int, noutputs::Int, dependency = nothing) where T = - InstantaneousContinuousMachine{T}(ninputs, 0, noutputs, (u,x,p,t)->T[], (u,x,p,t)->f(x), dependency) +InstantaneousContinuousMachine{T, I}(f::Function, ninputs::Int, noutputs::Int, dependency = nothing) where {T,I} = + InstantaneousContinuousMachine{T, I}(ninputs, 0, noutputs, (u,x,p,t)->T[], (u,x,p,t)->f(x), dependency) -InstantaneousContinuousMachine{T}(m::ContinuousMachine{T, I}) where {T, I<:DirectedInterface{T}} = +InstantaneousContinuousMachine{T, I}(m::ContinuousMachine{T, I}) where {T, I<:DirectedInterface{T}} = ContinuousMachine{T}(InstantaneousDirectedInterface{T}(input_ports(m), output_ports(m), []), ContinuousDirectedSystem{T}(nstates(m), dynamics(m), (u,x,p,t) -> readout(m, u, p, t)) ) @@ -186,7 +222,7 @@ The readout function `r` may depend on the state, history, parameters and time, If it is left out, then ``r=u``. """ const DelayMachine{T,I} = Machine{T, I, DelayDirectedSystem{T}} -const InstantaneousDelayMachine{T} = Machine{T, InstantaneousDirectedInterface{T}, DelayDirectedSystem{T}} +const InstantaneousDelayMachine{T, I} = Machine{T, I, DelayDirectedSystem{T}} DelayMachine{T}(interface::I, system::DelayDirectedSystem{T}) where {T, I<:AbstractDirectedInterface} = DelayMachine{T, I}(interface, system) @@ -209,6 +245,9 @@ DelayMachine{T,I}(ninputs, nstates, noutputs, dynamics, readout, dependency) whe InstantaneousDelayMachine{T}(ninputs, nstates, noutputs, dynamics, readout, dependency) where T = DelayMachine{T}(InstantaneousDirectedInterface{T}(ninputs, noutputs, dependency), DelayDirectedSystem{T}(nstates, dynamics, readout)) +InstantaneousDelayMachine{T,N}(ninputs, nstates, noutputs, dynamics, readout, dependency) where {T,N} = + DelayMachine{T}(InstantaneousDirectedVectorInterface{T,N}(ninputs, noutputs, dependency), DelayDirectedSystem{T}(nstates, dynamics, readout)) + InstantaneousDelayMachine{T}(f::Function, ninputs::Int, noutputs::Int, dependency = nothing) where T = InstantaneousDelayMachine{T}(ninputs, 0, noutputs, (u,x,p,t)->T[], (u,x,p,t)->f(x), dependency) @@ -273,11 +312,14 @@ readout(f::AbstractMachine, u::FinDomFunction, args...) = readout(f, collect(u), readout(m::Machine{T,I,S}, u::AbstractVector, x::AbstractVector, p=nothing, t=0) where {T, I<:InstantaneousDirectedInterface{T}, S} = readout(m)(u,x,p,t) +readout(m::Machine{T,I,S}, u::AbstractVector, x::AbstractVector, p=nothing, t=0) where {T, I<:InstantaneousDirectedVectorInterface, S} = + readout(m)(u,x,p,t) readout(f::DelayMachine{T,I}, u::AbstractVector, h=nothing, p = nothing, t = 0) where {T, I<:Union{DirectedInterface, DirectedVectorInterface}}= readout(f)(u, h, p, t) readout(m::DelayMachine{T,I}, u::AbstractVector, x::AbstractVector, h=nothing, p=nothing, t=0) where {T, I<:InstantaneousDirectedInterface} = readout(m)(u,x,h,p,t) - +readout(m::DelayMachine{T,I}, u::AbstractVector, x::AbstractVector, h=nothing, p=nothing, t=0) where {T, I<:InstantaneousDirectedVectorInterface} = + readout(m)(u,x,h,p,t) """ eval_dynamics(m::AbstractMachine, u::AbstractVector, x:AbstractVector{F}, p, t) where {F<:Function} eval_dynamics(m::AbstractMachine{T}, u::AbstractVector, x:AbstractVector{T}, p, t) where T diff --git a/test/dwd_dynam.jl b/test/dwd_dynam.jl index 8729fe3..3d7a268 100644 --- a/test/dwd_dynam.jl +++ b/test/dwd_dynam.jl @@ -310,6 +310,22 @@ end @test readout(m, vcat(u1, u2, u3), nothing, 0) == [u1 + u2, u2, u3] @test eval_dynamics(m, vcat(u1, u2, u3), [x1, x2], nothing, 0) == vcat(x1, u2 + x1 + 2*u2 + u3, x2) end + + @testset "VectorInterface for InstantaneousContinuousMachine" begin + m1 = InstantaneousContinuousMachine{Float64,3}( + 1, 3, 1, + (u, x, p, t) -> x*1.5, + (u, x, p, t) -> u+x, + [1 => 1] + ) + + x1 = [1,2,3] + u1 = [4,5,6] + + @test readout(m1, u1, x1, nothing, 0) == u1+x1 + @test eval_dynamics(m1, u1, [x1], nothing, 0) == x1*1.5 + + end end From 9ebea218f1f92aa8c4930dff4d15daea0205c876 Mon Sep 17 00:00:00 2001 From: slwu89 <10673535+slwu89@users.noreply.github.com> Date: Wed, 3 Aug 2022 17:11:36 -0700 Subject: [PATCH 2/6] fix errors and check InstantaneousContinuousMachine with vector readout works with single machine --- test/dwd_dynam.jl | 30 ++++++++++++++++-------------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/test/dwd_dynam.jl b/test/dwd_dynam.jl index 3d7a268..1edfb91 100644 --- a/test/dwd_dynam.jl +++ b/test/dwd_dynam.jl @@ -312,21 +312,23 @@ end end @testset "VectorInterface for InstantaneousContinuousMachine" begin - m1 = InstantaneousContinuousMachine{Float64,3}( + m = InstantaneousContinuousMachine{Float64,3}( 1, 3, 1, - (u, x, p, t) -> x*1.5, - (u, x, p, t) -> u+x, + (u, x, p, t) -> u*1.5, # dynamics + (u, x, p, t) -> u+x, # readout [1 => 1] ) - - x1 = [1,2,3] - u1 = [4,5,6] - - @test readout(m1, u1, x1, nothing, 0) == u1+x1 - @test eval_dynamics(m1, u1, [x1], nothing, 0) == x1*1.5 - + + x1 = Float64.([1,2,3]) + u1 = Float64.([4,5,6]) + + @test readout(m, u1, x1, nothing, 0) == u1+x1 + @test eval_dynamics(m, u1, [x1], nothing, 0) == u1*1.5 + + # compare to analytic soln + prob = ODEProblem(m, u1, [x1], (0.0, 0.05)) + sol = solve(prob, Tsit5()) + + @test sol.u[end] ≈ [x*exp(1.5*0.05) for x in u1] end -end - - - +end \ No newline at end of file From 13693d3fcce78f60b2f1d7c12aa38c2dc8202077 Mon Sep 17 00:00:00 2001 From: slwu89 <10673535+slwu89@users.noreply.github.com> Date: Wed, 3 Aug 2022 17:50:26 -0700 Subject: [PATCH 3/6] test composition works --- examples/Ross-Macdonald.jl | 5 ++++- src/dwd_dynam.jl | 18 +++++++++--------- test/dwd_dynam.jl | 24 +++++++++++++++++++++++- 3 files changed, 36 insertions(+), 11 deletions(-) diff --git a/examples/Ross-Macdonald.jl b/examples/Ross-Macdonald.jl index 8d6552d..f738c20 100644 --- a/examples/Ross-Macdonald.jl +++ b/examples/Ross-Macdonald.jl @@ -166,7 +166,7 @@ add_wires!(rmb, Pair[ # These two terms are sent from the bloodmeal machine to the mosquito and human machines # via their input ports. Then the dynamical system filling the mosquito machine is -# $\dot{Z} = a*\kappa*(e^{-gn} - Z) - gZ$ and $\dot{X} = bEIR(1-X) - rX$ is the dynamical system +# $\dot{Z} = a\kappa (e^{-gn} - Z) - gZ$ and $\dot{X} = bEIR(1-X) - rX$ is the dynamical system # filling the human machine. to_graphviz(rmb) @@ -270,6 +270,9 @@ plot(sol, label=["non-infectious mosquito population" "infectious mosquito popul xlabel = "time", ylabel = "proportion infectious", color = ["magenta" "red" "blue"] ) +N = length(sol) +plot!(sol.t, fill(X̄, N), label = "human equilibrium", ls = :dash, lw = 2, color = "blue") +plot!(sol.t, fill(Z̄, N), label = "infectious mosquito equilibrium", ls = :dash, lw = 2, color = "red") # While the equilibrium points of the two models are identical, they exhibit different dynamical behavior # before settling down to equilibrium. Because models are often used to examine how the system may respond diff --git a/src/dwd_dynam.jl b/src/dwd_dynam.jl index 562baa0..76be5c2 100644 --- a/src/dwd_dynam.jl +++ b/src/dwd_dynam.jl @@ -197,8 +197,8 @@ ContinuousMachine{T,I}(ninputs::Int, nstates::Int, dynamics) where {T,I} = ContinuousMachine{T,I}(ninputs, nstates, noutputs, dynamics, readout, dependency) where {T, I<:InstantaneousDirectedInterface{T}} = ContinuousMachine{T, I}(I(ninputs, noutputs, dependency), ContinuousDirectedSystem{T}(nstates, dynamics, readout)) -InstantaneousContinuousMachine{T}(ninputs, nstates, noutputs, dynamics, readout, dependency) where {T, I <: AbstractDirectedInterface} = - ContinuousMachine{T, I}(InstantaneousDirectedInterface{T}(ninputs, noutputs, dependency), ContinuousDirectedSystem{T}(nstates, dynamics, readout)) +InstantaneousContinuousMachine{T}(ninputs, nstates, noutputs, dynamics, readout, dependency) where {T} = + ContinuousMachine{T, InstantaneousDirectedInterface{T}}(InstantaneousDirectedInterface{T}(ninputs, noutputs, dependency), ContinuousDirectedSystem{T}(nstates, dynamics, readout)) InstantaneousContinuousMachine{T,N}(ninputs, nstates, noutputs, dynamics, readout, dependency) where {T,N} = ContinuousMachine{T, InstantaneousDirectedVectorInterface{T,N}}(InstantaneousDirectedVectorInterface{T,N}(ninputs, noutputs, dependency), ContinuousDirectedSystem{T}(nstates, dynamics, readout)) @@ -206,7 +206,7 @@ InstantaneousContinuousMachine{T,N}(ninputs, nstates, noutputs, dynamics, readou InstantaneousContinuousMachine{T, I}(f::Function, ninputs::Int, noutputs::Int, dependency = nothing) where {T,I} = InstantaneousContinuousMachine{T, I}(ninputs, 0, noutputs, (u,x,p,t)->T[], (u,x,p,t)->f(x), dependency) -InstantaneousContinuousMachine{T, I}(m::ContinuousMachine{T, I}) where {T, I<:DirectedInterface{T}} = +InstantaneousContinuousMachine{T}(m::ContinuousMachine{T, I}) where {T, I<:DirectedInterface{T}} = ContinuousMachine{T}(InstantaneousDirectedInterface{T}(input_ports(m), output_ports(m), []), ContinuousDirectedSystem{T}(nstates(m), dynamics(m), (u,x,p,t) -> readout(m, u, p, t)) ) @@ -242,16 +242,16 @@ DelayMachine{T,I}(ninputs::Int, nstates::Int, dynamics) where {T,I} = DelayMachine{T,I}(ninputs, nstates, noutputs, dynamics, readout, dependency) where {T, I<:InstantaneousDirectedInterface{T}} = DelayMachine{T, I}(I(ninputs, noutputs, dependency), DelayDirectedSystem{T}(nstates, dynamics, readout)) -InstantaneousDelayMachine{T}(ninputs, nstates, noutputs, dynamics, readout, dependency) where T = - DelayMachine{T}(InstantaneousDirectedInterface{T}(ninputs, noutputs, dependency), DelayDirectedSystem{T}(nstates, dynamics, readout)) +InstantaneousDelayMachine{T}(ninputs, nstates, noutputs, dynamics, readout, dependency) where {T} = + DelayMachine{T, InstantaneousDirectedInterface{T}}(InstantaneousDirectedInterface{T}(ninputs, noutputs, dependency), DelayDirectedSystem{T}(nstates, dynamics, readout)) InstantaneousDelayMachine{T,N}(ninputs, nstates, noutputs, dynamics, readout, dependency) where {T,N} = - DelayMachine{T}(InstantaneousDirectedVectorInterface{T,N}(ninputs, noutputs, dependency), DelayDirectedSystem{T}(nstates, dynamics, readout)) + DelayMachine{T,InstantaneousDirectedVectorInterface{T,N}}(InstantaneousDirectedVectorInterface{T,N}(ninputs, noutputs, dependency), DelayDirectedSystem{T}(nstates, dynamics, readout)) -InstantaneousDelayMachine{T}(f::Function, ninputs::Int, noutputs::Int, dependency = nothing) where T = - InstantaneousDelayMachine{T}(ninputs, 0, noutputs, (u,x,p,t)->T[], (u,x,p,t)->f(x), dependency) +InstantaneousDelayMachine{T, I}(f::Function, ninputs::Int, noutputs::Int, dependency = nothing) where {T,I} = + InstantaneousDelayMachine{T, I}(ninputs, 0, noutputs, (u,x,p,t)->T[], (u,x,p,t)->f(x), dependency) -InstantaneousDelayMachine(m::DelayMachine{T, I}) where {T, I<:DirectedInterface{T}} = +InstantaneousDelayMachine{T}(m::DelayMachine{T, I}) where {T, I<:DirectedInterface{T}} = DelayMachine{T}(InstantaneousDirectedInterface{T}(input_ports(m), output_ports(m), []), DelayDirectedSystem{T}(nstates(m), dynamics(m), (u,x,p,t) -> readout(m, u, p, t)) ) diff --git a/test/dwd_dynam.jl b/test/dwd_dynam.jl index 1edfb91..09bad7e 100644 --- a/test/dwd_dynam.jl +++ b/test/dwd_dynam.jl @@ -330,5 +330,27 @@ end sol = solve(prob, Tsit5()) @test sol.u[end] ≈ [x*exp(1.5*0.05) for x in u1] + + # check it composes + m1 = InstantaneousContinuousMachine{Float64,3}(2, 3, 1, + (u, x, p, t) -> x[1] + x[2], + (u,p,t) -> [u]) + m2 = InstantaneousContinuousMachine{Float64, 3}(1, 3, 2, + (u, x, p, t) -> x[1] + u, + (u,p,t) -> [u, 2*u]) + m3 = InstantaneousContinuousMachine{Float64, 3}(1, 3, 1, + (u, x, p, t) -> x[1], + (u,p,t) -> [u]) + m = oapply(d_big, [m1, m2, m3]) + + u1 = 1:3; u2 = 4:6; u3 = 7:9; + x1 = ones(3); x2 = fill(0.5, 3) + + @test readout(m, vcat(u1, u2, u3), nothing, 0) == [u1 + u2, u2, u3] + @test eval_dynamics(m, vcat(u1, u2, u3), [x1, x2], nothing, 0) == vcat(x1, u2 + x1 + 2*u2 + u3, x2) + end -end \ No newline at end of file +end + + + From 52009ef10882954fe22ba9dba21609651e405ab2 Mon Sep 17 00:00:00 2001 From: slwu89 <10673535+slwu89@users.noreply.github.com> Date: Thu, 4 Aug 2022 15:16:17 -0700 Subject: [PATCH 4/6] add test for InstantaneousDelayMachine with vectors on wires --- src/dwd_dynam.jl | 17 ++++++++--------- test/dwd_dynam.jl | 35 ++++++++++++++++++++++++++++++++--- 2 files changed, 40 insertions(+), 12 deletions(-) diff --git a/src/dwd_dynam.jl b/src/dwd_dynam.jl index 76be5c2..3665504 100644 --- a/src/dwd_dynam.jl +++ b/src/dwd_dynam.jl @@ -39,8 +39,11 @@ end DirectedVectorInterface{T,N}(ninputs::Int, noutputs::Int) where {T,N} = DirectedVectorInterface{T,N}(1:ninputs, 1:noutputs) +# Instantaneous interfaces +abstract type AbstractInstantaneousDirectedInterface{T} <: AbstractDirectedInterface{T} end + # InstantaneousDirectedInterface -struct InstantaneousDirectedInterface{T} <: AbstractDirectedInterface{T} +struct InstantaneousDirectedInterface{T} <: AbstractInstantaneousDirectedInterface{T} input_ports::Vector output_ports::Vector dependency::Span # P_in <- R -> P_out @@ -65,13 +68,8 @@ InstantaneousDirectedInterface{T}(input_ports::AbstractVector, output_ports::Abs end...) ) -dependency(interface::InstantaneousDirectedInterface) = interface.dependency -dependency_pairs(interface::InstantaneousDirectedInterface) = map(apex(dependency(interface))) do i - legs(dependency(interface))[2](i) => legs(dependency(interface))[1](i) -end |> sort - # InstantaneousDirectedVectorInterface -struct InstantaneousDirectedVectorInterface{T,N} <: AbstractDirectedInterface{T} +struct InstantaneousDirectedVectorInterface{T,N} <: AbstractInstantaneousDirectedInterface{T} input_ports::Vector output_ports::Vector dependency::Span # P_in <- R -> P_out @@ -96,8 +94,9 @@ InstantaneousDirectedVectorInterface{T,N}(input_ports::AbstractVector, output_po end...) ) -dependency(interface::InstantaneousDirectedVectorInterface) = interface.dependency -dependency_pairs(interface::InstantaneousDirectedVectorInterface) = map(apex(dependency(interface))) do i +# get dependency +dependency(interface::AbstractInstantaneousDirectedInterface) = interface.dependency +dependency_pairs(interface::AbstractInstantaneousDirectedInterface) = map(apex(dependency(interface))) do i legs(dependency(interface))[2](i) => legs(dependency(interface))[1](i) end |> sort diff --git a/test/dwd_dynam.jl b/test/dwd_dynam.jl index 09bad7e..f584b34 100644 --- a/test/dwd_dynam.jl +++ b/test/dwd_dynam.jl @@ -350,7 +350,36 @@ end @test eval_dynamics(m, vcat(u1, u2, u3), [x1, x2], nothing, 0) == vcat(x1, u2 + x1 + 2*u2 + u3, x2) end -end - - + @testset "VectorInterface for InstantaneousDelayMachine" begin + m = InstantaneousDelayMachine{Float64,3}( + 1, 3, 1, # ninputs, nstates, noutputs + (u, x, h, p, t) -> -h(p, t-1), # dynamics + (u, x, h, p, t) -> u+x, # readout + [1 => 1] # output -> input dependency + ) + + x1 = Float64.([1,2,3]) + u1 = Float64.([4,5,6]) + hist(p,t) = u1 + + @test readout(m, u1, x1, hist, nothing, 0) == u1+x1 + @test eval_dynamics(m, u1, [x1], hist, nothing, 0) == -u1 + + prob = DDEProblem(m, u1, [x1], hist, (0.0, 3.0), nothing) + sol = solve(prob,MethodOfSteps(Tsit5()),abstol=1e-12,reltol=1e-12) + + # DiffEq.jl solution + prob1 = DDEProblem((du,u,h,p,t) -> begin + x_lag = -h(p, t-1) + du[1] = x_lag[1] + du[2] = x_lag[2] + du[3] = x_lag[3] + end, u1, hist, (0.0, 3.0), nothing) + + sol1 = solve(prob1,MethodOfSteps(Tsit5()),abstol=1e-12,reltol=1e-12) + + @test sol.u[end] ≈ sol1.u[end] + + end +end \ No newline at end of file From 710b78cd573af57de973f54959f10b037e6298b4 Mon Sep 17 00:00:00 2001 From: slwu89 <10673535+slwu89@users.noreply.github.com> Date: Tue, 16 Aug 2022 21:12:08 -0700 Subject: [PATCH 5/6] ensure single-typed ports are default --- src/dwd_dynam.jl | 48 ++++++++++++++++++++++++++++++++++------------- test/dwd_dynam.jl | 8 ++++---- 2 files changed, 39 insertions(+), 17 deletions(-) diff --git a/src/dwd_dynam.jl b/src/dwd_dynam.jl index 3665504..d64bddd 100644 --- a/src/dwd_dynam.jl +++ b/src/dwd_dynam.jl @@ -199,11 +199,20 @@ ContinuousMachine{T,I}(ninputs, nstates, noutputs, dynamics, readout, dependency InstantaneousContinuousMachine{T}(ninputs, nstates, noutputs, dynamics, readout, dependency) where {T} = ContinuousMachine{T, InstantaneousDirectedInterface{T}}(InstantaneousDirectedInterface{T}(ninputs, noutputs, dependency), ContinuousDirectedSystem{T}(nstates, dynamics, readout)) +InstantaneousContinuousMachine{T,I}(ninputs, nstates, noutputs, dynamics, readout, dependency) where {T,I<:InstantaneousDirectedInterface{T}} = + ContinuousMachine{T,I}(I(ninputs, noutputs, dependency), ContinuousDirectedSystem{T}(nstates, dynamics, readout)) + InstantaneousContinuousMachine{T,N}(ninputs, nstates, noutputs, dynamics, readout, dependency) where {T,N} = ContinuousMachine{T, InstantaneousDirectedVectorInterface{T,N}}(InstantaneousDirectedVectorInterface{T,N}(ninputs, noutputs, dependency), ContinuousDirectedSystem{T}(nstates, dynamics, readout)) -InstantaneousContinuousMachine{T, I}(f::Function, ninputs::Int, noutputs::Int, dependency = nothing) where {T,I} = - InstantaneousContinuousMachine{T, I}(ninputs, 0, noutputs, (u,x,p,t)->T[], (u,x,p,t)->f(x), dependency) +InstantaneousContinuousMachine{T}(f::Function, ninputs::Int, noutputs::Int, dependency = nothing) where {T} = + ContinuousMachine{T,InstantaneousDirectedInterface{T}}(ninputs, 0, noutputs, (u,x,p,t)->f(x), (u,x,p,t)->T[], dependency) + +InstantaneousContinuousMachine{T,I}(f::Function, ninputs::Int, noutputs::Int, dependency = nothing) where {T,I<:AbstractInstantaneousDirectedInterface{T}} = + ContinuousMachine{T,I}(ninputs, 0, noutputs, (u,x,p,t)->f(x), (u,x,p,t)->T[], dependency) + +InstantaneousContinuousMachine{T,N}(f::Function, ninputs::Int, noutputs::Int, dependency = nothing) where {T,N} = + InstantaneousContinuousMachine{T,N}(ninputs, 0, noutputs, (u,x,p,t)->f(x), (u,x,p,t)->T[], dependency) InstantaneousContinuousMachine{T}(m::ContinuousMachine{T, I}) where {T, I<:DirectedInterface{T}} = ContinuousMachine{T}(InstantaneousDirectedInterface{T}(input_ports(m), output_ports(m), []), @@ -221,7 +230,7 @@ The readout function `r` may depend on the state, history, parameters and time, If it is left out, then ``r=u``. """ const DelayMachine{T,I} = Machine{T, I, DelayDirectedSystem{T}} -const InstantaneousDelayMachine{T, I} = Machine{T, I, DelayDirectedSystem{T}} +const InstantaneousDelayMachine{T,I} = Machine{T, I, DelayDirectedSystem{T}} DelayMachine{T}(interface::I, system::DelayDirectedSystem{T}) where {T, I<:AbstractDirectedInterface} = DelayMachine{T, I}(interface, system) @@ -244,15 +253,24 @@ DelayMachine{T,I}(ninputs, nstates, noutputs, dynamics, readout, dependency) whe InstantaneousDelayMachine{T}(ninputs, nstates, noutputs, dynamics, readout, dependency) where {T} = DelayMachine{T, InstantaneousDirectedInterface{T}}(InstantaneousDirectedInterface{T}(ninputs, noutputs, dependency), DelayDirectedSystem{T}(nstates, dynamics, readout)) +InstantaneousDelayMachine{T,I}(ninputs, nstates, noutputs, dynamics, readout, dependency) where {T,I<:InstantaneousDirectedInterface{T}} = + DelayMachine{T,I}(I(ninputs, noutputs, dependency), DelayDirectedSystem{T}(nstates, dynamics, readout)) + InstantaneousDelayMachine{T,N}(ninputs, nstates, noutputs, dynamics, readout, dependency) where {T,N} = - DelayMachine{T,InstantaneousDirectedVectorInterface{T,N}}(InstantaneousDirectedVectorInterface{T,N}(ninputs, noutputs, dependency), DelayDirectedSystem{T}(nstates, dynamics, readout)) + DelayMachine{T, InstantaneousDirectedVectorInterface{T,N}}(InstantaneousDirectedVectorInterface{T,N}(ninputs, noutputs, dependency), DelayDirectedSystem{T}(nstates, dynamics, readout)) -InstantaneousDelayMachine{T, I}(f::Function, ninputs::Int, noutputs::Int, dependency = nothing) where {T,I} = - InstantaneousDelayMachine{T, I}(ninputs, 0, noutputs, (u,x,p,t)->T[], (u,x,p,t)->f(x), dependency) +InstantaneousDelayMachine{T}(f::Function, ninputs::Int, noutputs::Int, dependency = nothing) where {T} = + DelayMachine{T,InstantaneousDirectedInterface{T}}(ninputs, 0, noutputs, (u,x,h,p,t)->f(x), (u,x,h,p,t)->T[], dependency) + +InstantaneousDelayMachine{T,I}(f::Function, ninputs::Int, noutputs::Int, dependency = nothing) where {T,I<:AbstractInstantaneousDirectedInterface{T}} = + DelayMachine{T,I}(ninputs, 0, noutputs, (u,x,h,p,t)->f(x), (u,x,h,p,t)->T[], dependency) + +InstantaneousDelayMachine{T,N}(f::Function, ninputs::Int, noutputs::Int, dependency = nothing) where {T,N} = + InstantaneousDelayMachine{T,N}(ninputs, 0, noutputs, (u,x,h,p,t)->f(x), (u,x,h,p,t)->T[], dependency) InstantaneousDelayMachine{T}(m::DelayMachine{T, I}) where {T, I<:DirectedInterface{T}} = DelayMachine{T}(InstantaneousDirectedInterface{T}(input_ports(m), output_ports(m), []), - DelayDirectedSystem{T}(nstates(m), dynamics(m), (u,x,p,t) -> readout(m, u, p, t)) + DelayDirectedSystem{T}(nstates(m), dynamics(m), (u,x,h,p,t) -> readout(m, u, p, t)) ) @@ -309,16 +327,20 @@ eltype(::AbstractMachine{T}) where T = T readout(f::AbstractMachine, u::AbstractVector, p = nothing, t = 0) = readout(f)(u, p, t) readout(f::AbstractMachine, u::FinDomFunction, args...) = readout(f, collect(u), args...) -readout(m::Machine{T,I,S}, u::AbstractVector, x::AbstractVector, p=nothing, t=0) where {T, I<:InstantaneousDirectedInterface{T}, S} = - readout(m)(u,x,p,t) -readout(m::Machine{T,I,S}, u::AbstractVector, x::AbstractVector, p=nothing, t=0) where {T, I<:InstantaneousDirectedVectorInterface, S} = +readout(m::Machine{T,I,S}, u::AbstractVector, x::AbstractVector, p=nothing, t=0) where {T, I<:AbstractInstantaneousDirectedInterface{T}, S} = readout(m)(u,x,p,t) +# readout(m::Machine{T,I,S}, u::AbstractVector, x::AbstractVector, p=nothing, t=0) where {T, I<:InstantaneousDirectedInterface{T}, S} = +# readout(m)(u,x,p,t) +# readout(m::Machine{T,I,S}, u::AbstractVector, x::AbstractVector, p=nothing, t=0) where {T, I<:InstantaneousDirectedVectorInterface, S} = +# readout(m)(u,x,p,t) readout(f::DelayMachine{T,I}, u::AbstractVector, h=nothing, p = nothing, t = 0) where {T, I<:Union{DirectedInterface, DirectedVectorInterface}}= readout(f)(u, h, p, t) -readout(m::DelayMachine{T,I}, u::AbstractVector, x::AbstractVector, h=nothing, p=nothing, t=0) where {T, I<:InstantaneousDirectedInterface} = - readout(m)(u,x,h,p,t) -readout(m::DelayMachine{T,I}, u::AbstractVector, x::AbstractVector, h=nothing, p=nothing, t=0) where {T, I<:InstantaneousDirectedVectorInterface} = +readout(m::DelayMachine{T,I}, u::AbstractVector, x::AbstractVector, h=nothing, p=nothing, t=0) where {T, I<:AbstractInstantaneousDirectedInterface} = readout(m)(u,x,h,p,t) +# readout(m::DelayMachine{T,I}, u::AbstractVector, x::AbstractVector, h=nothing, p=nothing, t=0) where {T, I<:InstantaneousDirectedInterface} = +# readout(m)(u,x,h,p,t) +# readout(m::DelayMachine{T,I}, u::AbstractVector, x::AbstractVector, h=nothing, p=nothing, t=0) where {T, I<:InstantaneousDirectedVectorInterface} = +# readout(m)(u,x,h,p,t) """ eval_dynamics(m::AbstractMachine, u::AbstractVector, x:AbstractVector{F}, p, t) where {F<:Function} eval_dynamics(m::AbstractMachine{T}, u::AbstractVector, x:AbstractVector{T}, p, t) where T diff --git a/test/dwd_dynam.jl b/test/dwd_dynam.jl index f584b34..2ebd895 100644 --- a/test/dwd_dynam.jl +++ b/test/dwd_dynam.jl @@ -315,14 +315,14 @@ end m = InstantaneousContinuousMachine{Float64,3}( 1, 3, 1, (u, x, p, t) -> u*1.5, # dynamics - (u, x, p, t) -> u+x, # readout + (u, x, p, t) -> u+x[1], # readout [1 => 1] ) x1 = Float64.([1,2,3]) u1 = Float64.([4,5,6]) - @test readout(m, u1, x1, nothing, 0) == u1+x1 + @test readout(m, u1, [x1], nothing, 0) == u1+x1 @test eval_dynamics(m, u1, [x1], nothing, 0) == u1*1.5 # compare to analytic soln @@ -355,7 +355,7 @@ end m = InstantaneousDelayMachine{Float64,3}( 1, 3, 1, # ninputs, nstates, noutputs (u, x, h, p, t) -> -h(p, t-1), # dynamics - (u, x, h, p, t) -> u+x, # readout + (u, x, h, p, t) -> u+x[1], # readout [1 => 1] # output -> input dependency ) @@ -363,7 +363,7 @@ end u1 = Float64.([4,5,6]) hist(p,t) = u1 - @test readout(m, u1, x1, hist, nothing, 0) == u1+x1 + @test readout(m, u1, [x1], hist, nothing, 0) == u1+x1 @test eval_dynamics(m, u1, [x1], hist, nothing, 0) == -u1 prob = DDEProblem(m, u1, [x1], hist, (0.0, 3.0), nothing) From d09c7316f30aaeef5e5d1d6cc8ecbf8c404447ab Mon Sep 17 00:00:00 2001 From: slwu89 <10673535+slwu89@users.noreply.github.com> Date: Tue, 16 Aug 2022 21:15:17 -0700 Subject: [PATCH 6/6] remove commented out code --- src/dwd_dynam.jl | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/src/dwd_dynam.jl b/src/dwd_dynam.jl index d64bddd..79d9b3f 100644 --- a/src/dwd_dynam.jl +++ b/src/dwd_dynam.jl @@ -329,18 +329,11 @@ readout(f::AbstractMachine, u::FinDomFunction, args...) = readout(f, collect(u), readout(m::Machine{T,I,S}, u::AbstractVector, x::AbstractVector, p=nothing, t=0) where {T, I<:AbstractInstantaneousDirectedInterface{T}, S} = readout(m)(u,x,p,t) -# readout(m::Machine{T,I,S}, u::AbstractVector, x::AbstractVector, p=nothing, t=0) where {T, I<:InstantaneousDirectedInterface{T}, S} = -# readout(m)(u,x,p,t) -# readout(m::Machine{T,I,S}, u::AbstractVector, x::AbstractVector, p=nothing, t=0) where {T, I<:InstantaneousDirectedVectorInterface, S} = -# readout(m)(u,x,p,t) readout(f::DelayMachine{T,I}, u::AbstractVector, h=nothing, p = nothing, t = 0) where {T, I<:Union{DirectedInterface, DirectedVectorInterface}}= readout(f)(u, h, p, t) readout(m::DelayMachine{T,I}, u::AbstractVector, x::AbstractVector, h=nothing, p=nothing, t=0) where {T, I<:AbstractInstantaneousDirectedInterface} = readout(m)(u,x,h,p,t) -# readout(m::DelayMachine{T,I}, u::AbstractVector, x::AbstractVector, h=nothing, p=nothing, t=0) where {T, I<:InstantaneousDirectedInterface} = -# readout(m)(u,x,h,p,t) -# readout(m::DelayMachine{T,I}, u::AbstractVector, x::AbstractVector, h=nothing, p=nothing, t=0) where {T, I<:InstantaneousDirectedVectorInterface} = -# readout(m)(u,x,h,p,t) + """ eval_dynamics(m::AbstractMachine, u::AbstractVector, x:AbstractVector{F}, p, t) where {F<:Function} eval_dynamics(m::AbstractMachine{T}, u::AbstractVector, x:AbstractVector{T}, p, t) where T