Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

pass vectors along wires for instantaneous machines #112

Merged
merged 6 commits into from
Aug 17, 2022
Merged

pass vectors along wires for instantaneous machines #112

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
81c44a1
(not working) attempt to add vector wires for instantaneous continuou…
slwu89 Aug 3, 2022
9ebea21
fix errors and check InstantaneousContinuousMachine with vector reado…
slwu89 Aug 4, 2022
13693d3
test composition works
slwu89 Aug 4, 2022
52009ef
add test for InstantaneousDelayMachine with vectors on wires
slwu89 Aug 4, 2022
710b78c
ensure single-typed ports are default
slwu89 Aug 17, 2022
d09c731
remove commented out code
slwu89 Aug 17, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 4 additions & 1 deletion examples/Ross-Macdonald.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
73 changes: 57 additions & 16 deletions src/dwd_dynam.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -65,21 +68,50 @@ InstantaneousDirectedInterface{T}(input_ports::AbstractVector, output_ports::Abs
end...)
)

dependency(interface::InstantaneousDirectedInterface) = interface.dependency
dependency_pairs(interface::InstantaneousDirectedInterface) = map(apex(dependency(interface))) do i
# InstantaneousDirectedVectorInterface
struct InstantaneousDirectedVectorInterface{T,N} <: AbstractInstantaneousDirectedInterface{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...)
)

# 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

# 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
Expand Down Expand Up @@ -144,7 +176,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}}
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Check here that the given interface I is a subtype of AbstractInstantaneousDirectedInterface{T}

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@slibkind I tried to implement something like const InstantaneousContinuousMachine{T,I<: AbstractInstantaneousDirectedInterface{T}} = Machine{T, I, ContinuousDirectedSystem{T}} but it ran into problems with the constructors as implemented. That is, the constructor typed as {T,N} complained when N was given as an Int64 for the dimensions of the vectors on the ports, because it expected an interface. So I did the checking in the constructors.

I think this is a general problem we have to live with because when I tried to do some type checking for regular continuous machines as const ContinuousMachine{T,I <: AbstractDirectedInterface{T}} = Machine{T, I, ContinuousDirectedSystem{T}} the constructor typed as {T,N} also failed with the type error ERROR: TypeError: in Type, in I, expected I<:AlgebraicDynamics.DWDDynam.AbstractDirectedInterface{Float64}, got a value of type Int64.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I agree that having the vector length as part of the type parameter has led to some clunky issues. I made note of this in Issue #113.


ContinuousMachine{T}(interface::I, system::ContinuousDirectedSystem{T}) where {T, I <: AbstractDirectedInterface} =
ContinuousMachine{T, I}(interface, system)
Expand All @@ -164,11 +196,14 @@ 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} =
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))

InstantaneousContinuousMachine{T}(f::Function, ninputs::Int, noutputs::Int, dependency = nothing) where T =
slibkind marked this conversation as resolved.
Show resolved Hide resolved
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}} =
ContinuousMachine{T}(InstantaneousDirectedInterface{T}(input_ports(m), output_ports(m), []),
Expand All @@ -186,7 +221,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}}
slibkind marked this conversation as resolved.
Show resolved Hide resolved
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)
Expand All @@ -206,13 +241,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}}(InstantaneousDirectedVectorInterface{T,N}(ninputs, noutputs, dependency), DelayDirectedSystem{T}(nstates, dynamics, readout))

InstantaneousDelayMachine{T}(f::Function, ninputs::Int, noutputs::Int, dependency = nothing) where T =
slibkind marked this conversation as resolved.
Show resolved Hide resolved
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))
)
Expand Down Expand Up @@ -273,11 +311,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} =
slibkind marked this conversation as resolved.
Show resolved Hide resolved
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)
slibkind marked this conversation as resolved.
Show resolved Hide resolved

""" 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
Expand Down
71 changes: 70 additions & 1 deletion test/dwd_dynam.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -310,7 +310,76 @@ 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
end

@testset "VectorInterface for InstantaneousContinuousMachine" begin
m = InstantaneousContinuousMachine{Float64,3}(
1, 3, 1,
(u, x, p, t) -> u*1.5, # dynamics
(u, x, p, t) -> u+x, # readout
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm surprised this works, since u is a vector of length 3 and x is a vector (of length 1) of vectors of length 3.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fixed.

[1 => 1]
)

x1 = Float64.([1,2,3])
u1 = Float64.([4,5,6])

@test readout(m, u1, x1, nothing, 0) == u1+x1
slibkind marked this conversation as resolved.
Show resolved Hide resolved
@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]

# 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

@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