Batch solving

src/AlgebraicDynamics.jl

@@ -8,6 +8,5 @@ using Catlab.Programs
 include("uwd_dynam.jl")
 include("dwd_dynam.jl")
 include("cpg_dynam.jl")
-include("trajectories.jl")
 
 end # module

src/dwd_dynam.jl

@@ -43,8 +43,8 @@ noutputs(interface::AbstractDirectedInterface) = length(output_ports(interface))
 
 ndims(::DirectedVectorInterface{T, N}) where {T,N} = N
 
-unit(::Type{I}, dims) where {T, I<:AbstractDirectedInterface{T}} = zeros(T, dims)
-unit(::Type{DirectedVectorInterface{T,N}}, dims) where {T, N} = fill(zeros(T, N), dims)
+unit(::Type{I}, shape) where {T, I<:AbstractDirectedInterface{T}} = fill(zero(T), shape)
+unit(::Type{DirectedVectorInterface{T,N}}, shape) where {T, N} = fill(zeros(T, N), shape)
 
 
 ### Dynamics
@@ -188,30 +188,30 @@ show(io::IO, vf::DiscreteMachine) = print("DiscreteMachine(ℝ^$(nstates(vf)) ×
 eltype(::AbstractMachine{T}) where T = T
 
 
-readout(f::DelayMachine, u::AbstractVector, h = nothing, p = nothing, t = 0) = readout(f)(u, h, p, 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(f::DelayMachine, u, h = nothing, p = nothing, t = 0) = readout(f)(collect(u), h, p, t)
+readout(f::AbstractMachine, u, p = nothing, t = 0) = readout(f)(collect(u), p, t)
 
 """    eval_dynamics(m::AbstractMachine, u::AbstractVector, xs:AbstractVector, p, t)
 
 Evaluates the dynamics of the machine `m` at state `u`, parameters `p`, and time `t`. The exogenous variables are set by `xs` which may either be a collection of functions ``x(t)`` or a collection of constant values. 
 
 The length of `xs` must equal the number of inputs to `m`.
 """
+
 eval_dynamics(f::DelayMachine, u, xs, h, p=nothing, t=0) = begin
-    ninputs(f) == length(xs) || error("$xs must have length $(ninputs(f)) to set the exogenous variables.")
+    #ninputs(f) == length(xs) || error("$xs must have length $(ninputs(f)) to set the exogenous variables.")
     dynamics(f)(collect(u), collect(xs), h, p, t)
 end
 
 eval_dynamics(f::AbstractMachine, u, xs, p=nothing, t=0) = begin
-    ninputs(f) == length(xs) || error("$xs must have length $(ninputs(f)) to set the exogenous variables.")
+    #ninputs(f) == length(xs) || error("$xs must have length $(ninputs(f)) to set the exogenous variables.")
     dynamics(f)(collect(u), collect(xs), p, t)
 end
 
 # eval_dynamics(f::AbstractMachine, u::S, xs::T, args...) where {S,T <: Union{FinDomFunction, AbstractVector}} =
 #     eval_dynamics(f, collect(u), collect(xs), args...)
 
-eval_dynamics(f::AbstractMachine, u::AbstractVector, xs::AbstractVector{T}, p=nothing, t=0) where T <: Function =
+eval_dynamics(f::AbstractMachine, u, xs::AbstractVector{T}, p=nothing, t=0) where T <: Function =
     eval_dynamics(f, u, [x(t) for x in xs], p, t)
 
 """    euler_approx(m::ContinuousMachine, h)
@@ -246,18 +246,18 @@ euler_approx(fs::AbstractDict{S, M}, args...) where {S, M<:ContinuousMachine} =
 
 Constructs an ODEProblem from the vector field defined by `(u,p,t) -> m.dynamics(u, x, p, t)`. The exogenous variables are determined by `xs`.
 """
-ODEProblem(m::ContinuousMachine{T}, u0::AbstractVector, xs::AbstractVector, tspan, p=nothing; kwargs...)  where T= 
+ODEProblem(m::ContinuousMachine{T}, u0, xs, tspan, p=nothing; kwargs...)  where T= 
     ODEProblem((u,p,t) -> eval_dynamics(m, u, xs, p, t), u0, tspan, p; kwargs...)
 
-ODEProblem(m::ContinuousMachine{T}, u0::AbstractVector, x::Union{T, Function}, tspan, p=nothing; kwargs...) where T= 
+ODEProblem(m::ContinuousMachine{T}, u0, x::Union{T, Function}, tspan, p=nothing; kwargs...) where T= 
     ODEProblem(m, u0, collect(repeated(x, ninputs(m))), tspan, p; kwargs...)
 
-ODEProblem(m::ContinuousMachine{T}, u0::AbstractVector, tspan, p=nothing; kwargs...) where T = 
+ODEProblem(m::ContinuousMachine{T}, u0, tspan::Tuple, p=nothing; kwargs...) where T = 
     ODEProblem(m, u0, T[], tspan, p; kwargs...)
 
 """    DDEProblem(m::DelayMachine, u0::Vector, xs::Vector, h::Function, tspan, p = nothing; kwargs...)
 """
-DDEProblem(m::DelayMachine, u0::AbstractVector, xs::AbstractVector, hist, tspan, params=nothing; kwargs...) = 
+DDEProblem(m::DelayMachine, u0, xs::AbstractVector, hist, tspan, params=nothing; kwargs...) = 
     DDEProblem((u,h,p,t) -> eval_dynamics(m, u, xs, h, p, t), u0, hist, tspan, params; kwargs...)
 
 
@@ -266,10 +266,10 @@ DDEProblem(m::DelayMachine, u0::AbstractVector, xs::AbstractVector, hist, tspan,
 Constructs an DiscreteDynamicalSystem from the equation of motion defined by 
 `(u,p,t) -> m.dynamics(u, x, p, t)`. The exogenous variables are determined by `xs`. Pass `nothing` in place of `p` if your system does not have parameters.
 """
-DiscreteProblem(m::DiscreteMachine, u0::AbstractVector, xs::AbstractVector, tspan, p; kwargs...) = 
+DiscreteProblem(m::DiscreteMachine, u0, xs::AbstractVector, tspan, p; kwargs...) = 
     DiscreteProblem((u,p,t) -> eval_dynamics(m, u, xs, p, t), u0, tspan, p; kwargs...)
 
-DiscreteProblem(m::DiscreteMachine, u0::AbstractVector, x, tspan, p; kwargs...) = 
+DiscreteProblem(m::DiscreteMachine, u0, x, tspan, p; kwargs...) = 
   DiscreteProblem(m, u0, collect(repeated(x, ninputs(m))), tspan, p; kwargs...)
 
 DiscreteProblem(m::DiscreteMachine{T}, u0, tspan, p; kwargs...) where T = 
@@ -347,26 +347,27 @@ end
 
 function induced_dynamics(d::WiringDiagram, ms::Vector{M}, S, Inputs) where {T,I, M<:AbstractMachine{T,I}}
 
-    function v(u::AbstractVector, xs::AbstractVector, p, t::Real)  
+    function v(u, xs, p, t)  
         states = destruct(S, u) # a list of the states by box
         readouts = get_readouts(ms, states, p, t)
-        readins = unit(I, length(apex(Inputs)))
+        readin_shape = u isa AbstractVector ? length(apex(Inputs)) : (length(apex(Inputs)), size(u,2))
+        readins = unit(I, readin_shape)
         fill_readins!(readins, d, Inputs, readouts, xs)
 
         reduce(vcat, map(enumerate(destruct(Inputs, readins))) do (i,x)
             eval_dynamics(ms[i], states[i], x, p, t)
         end)
     end 
-
 end
 
 function induced_dynamics(d::WiringDiagram, ms::Vector{M}, S, Inputs) where {T,I, M<:DelayMachine{T,I}}
 
-    function v(u::AbstractVector, xs::AbstractVector, h, p, t::Real) 
+    function v(u, xs, h, p, t) 
         states = destruct(S, u) # a list of the states by box
         hists = destruct(S, h)
         readouts = get_readouts(ms, states, hists, p, t)
-        readins = unit(I, length(apex(Inputs)))
+        readin_shape = u isa AbstractVector ? length(apex(Inputs)) : (length(apex(Inputs)), size(u,2))
+        readins = unit(I, readin_shape)
         fill_readins!(readins, d, Inputs, readouts, xs)
 
         reduce(vcat, map(enumerate(destruct(Inputs, readins))) do (i,x)
@@ -376,20 +377,22 @@ function induced_dynamics(d::WiringDiagram, ms::Vector{M}, S, Inputs) where {T,I
 end 
 
 function induced_readout(d::WiringDiagram, ms::Vector{M}, S) where {T, I, M<:AbstractMachine{T,I}}
-    function r(u::AbstractVector, p, t)
+    function r(u, p, t)
         states = destruct(S, u)
         readouts = get_readouts(ms, states, p, t)
-        outputs = unit(I, length(output_ports(d)))
+        readout_shape = u isa AbstractVector ? length(output_ports(d)) : (length(output_ports(d)), size(u,2))
+        outputs = unit(I, readout_shape)
         fill_outputs!(outputs, d, readouts)
     end
 end
 
 function induced_readout(d::WiringDiagram, ms::Vector{M}, S) where {T, I, M<:DelayMachine{T,I}}
-    function r(u::AbstractVector, h, p, t)
+    function r(u, h, p, t)
         states = destruct(S, u)
         hists = destruct(S, h)
         readouts = get_readouts(ms, states, hists, p, t)
-        outputs = unit(I, length(output_ports(d)))
+        readout_shape = u isa AbstractVector ? length(output_ports(d)) : (length(output_ports(d)), size(u,2))
+        outputs = unit(I, readout_shape)
         fill_outputs!(outputs, d, readouts)
     end
 end
@@ -406,10 +409,13 @@ function fills(m::AbstractMachine, d::WiringDiagram, b::Int)
 end
 
 
-destruct(C::Colimit, xs::FinDomFunction) = map(1:length(C)) do i 
-    collect(compose(legs(C)[i], xs))
+destruct(C::Colimit, xs::AbstractVector) = map(1:length(C)) do i 
+  xs[legs(C)[i].func]
+end
+
+destruct(C::Colimit, xs::AbstractMatrix) = map(1:length(C)) do i 
+  xs[legs(C)[i].func, :]
 end
-destruct(C::Colimit, xs::AbstractVector) = destruct(C, FinDomFunction(xs))
 
 destruct(C::Colimit, h) = map(1:length(C)) do i 
     (p,t) -> destruct(C, h(p,t))[i]
@@ -425,12 +431,11 @@ end
 
 
 function fill_readins!(readins, d::WiringDiagram, Inputs::Colimit, readouts, xs) 
-
     for w in wires(d, :Wire)
-        readins[legs(Inputs)[w.target.box](w.target.port)] += readouts[w.source.box][w.source.port]
+        readins[legs(Inputs)[w.target.box](w.target.port), :] += readouts[w.source.box][w.source.port, :]
     end
     for w in wires(d, :InWire)
-        readins[legs(Inputs)[w.target.box](w.target.port)] += xs[w.source.port]
+        readins[legs(Inputs)[w.target.box](w.target.port), :] += xs[w.source.port, :]
     end
 
     return readins
@@ -440,7 +445,7 @@ end
 function fill_outputs!(outs, d::WiringDiagram, readouts) 
 
     for w in wires(d, :OutWire)
-        outs[w.target.port] += readouts[w.source.box][w.source.port]
+        outs[w.target.port, :] += readouts[w.source.box][w.source.port, :]
     end
 
     return outs

src/uwd_dynam.jl

@@ -133,7 +133,7 @@ DelayResourceSharer{T}(nports, nstates, dynamics, portmap) where {T}=
   DelayResourceSharer{T}(UndirectedInterface{T}(nports), DelayUndirectedSystem{T}(nstates, dynamics, portmap))
 
 DelayResourceSharer{T}(nstates::Int, dynamics::Function) where T = 
-    DelayResourceSharer{T}(nstates,nstates, dynamics, 1:nstates)
+    DelayResourceSharer{T}(nstates, nstates, dynamics, 1:nstates)
 
 """ 
 
@@ -153,24 +153,24 @@ DiscreteResourceSharer{T}(nports, nstates, dynamics, portmap) where {T}=
   DiscreteResourceSharer{T}(UndirectedInterface{T}(nports), DiscreteUndirectedSystem{T}(nstates, dynamics, portmap))
 
 DiscreteResourceSharer{T}(nstates::Int, dynamics::Function) where T = 
-    DiscreteResourceSharer{T}(nstates,nstates, dynamics, 1:nstates)
+    DiscreteResourceSharer{T}(nstates, nstates, dynamics, 1:nstates)
 
 """    eval_dynamics(r::AbstractResourceSharer, u::AbstractVector, p, t)
 
 Evaluates the dynamics of the resource sharer `r` at state `u`, parameters `p`, and time `t`.
 
 Omitting `t` and `p` is allowed if the dynamics of `r` does not depend on them.
 """
-eval_dynamics(r::DelayResourceSharer, u::AbstractVector, h, p, t::Real) = dynamics(r)(u, h, p, t)
-eval_dynamics!(du, r::DelayResourceSharer, u::AbstractVector, h, p, t::Real) = begin
+eval_dynamics(r::DelayResourceSharer, u, h, p, t) = dynamics(r)(u, h, p, t)
+eval_dynamics!(du, r::DelayResourceSharer, u, h, p, t) = begin
     du .= eval_dynamics(r, u, h, p, t)
 end
-eval_dynamics(r::AbstractResourceSharer, u::AbstractVector, p, t::Real) = dynamics(r)(u, p, t)
-eval_dynamics!(du, r::AbstractResourceSharer, u::AbstractVector, p, t::Real) = begin
+eval_dynamics(r::AbstractResourceSharer, u, p, t) = dynamics(r)(u, p, t)
+eval_dynamics!(du, r::AbstractResourceSharer, u, p, t) = begin
     du .= eval_dynamics(r, u, p, t)
 end
-eval_dynamics(r::AbstractResourceSharer, u::AbstractVector) = eval_dynamics(r, u, [], 0)
-eval_dynamics(r::AbstractResourceSharer, u::AbstractVector, p) = eval_dynamics(r, u, p, 0)
+eval_dynamics(r::AbstractResourceSharer, u) = eval_dynamics(r, u, [], 0)
+eval_dynamics(r::AbstractResourceSharer, u, p) = eval_dynamics(r, u, p, 0)
 
 show(io::IO, vf::ContinuousResourceSharer) = print("ContinuousResourceSharer(ℝ^$(nstates(vf)) → ℝ^$(nstates(vf))) with $(nports(vf)) exposed ports")
 show(io::IO, vf::DelayResourceSharer) = print("DelayResourceSharer(ℝ^$(nstates(vf)) → ℝ^$(nstates(vf))) with $(nports(vf)) exposed ports")
@@ -181,7 +181,7 @@ eltype(r::AbstractResourceSharer{T}) where T = T
 
 Transforms a continuous resource sharer `r` into a discrete resource sharer via Euler's method with step size `h`. If the dynamics of `r` is given by ``\\dot{u}(t) = f(u(t),p,t)`` the the dynamics of the new discrete system is given by the update rule ``u_{n+1} = u_n + h f(u_n, p, t)``.
 """
-euler_approx(f::ContinuousResourceSharer{T}, h::Float64) where T = DiscreteResourceSharer{T}(
+euler_approx(f::ContinuousResourceSharer{T}, h::Real) where T = DiscreteResourceSharer{T}(
         nports(f), nstates(f), 
         (u, p, t) -> u + h*eval_dynamics(f, u, p, t),
         portmap(f)
@@ -207,28 +207,28 @@ euler_approx(fs::AbstractDict{S, R}, args...) where {S, T, R<:ContinuousResource
 
 Constructs an ODEProblem from the vector field defined by `r.dynamics(u,p,t)`.
 """
-ODEProblem(r::ContinuousResourceSharer, u0::AbstractVector, tspan, p=nothing; kwargs...) = 
+ODEProblem(r::ContinuousResourceSharer, u0, tspan, p=nothing; kwargs...) = 
     ODEProblem(dynamics(r), u0, tspan, p; kwargs...)
 
 """    DDEProblem(r::DelayResourceSharer, u0::Vector, h, tspan)
 
 Constructs an DDEProblem from the vector field defined by `r.dynamics(u,h,p,t)`.
 """
-DDEProblem(r::DelayResourceSharer, u0::AbstractVector, h, tspan, p=nothing; kwargs...) = 
+DDEProblem(r::DelayResourceSharer, u0, h, tspan, p=nothing; kwargs...) = 
     DDEProblem(dynamics(r), u0, h, tspan, p; kwargs...)
 
 """    DiscreteProblem(r::DiscreteResourceSharer, u0::Vector, p)
 
 Constructs a DiscreteDynamicalSystem from the equation of motion `r.dynamics(u,p,t)`.  Pass `nothing` in place of `p` if your system does not have parameters.
 """
-DiscreteProblem(r::DiscreteResourceSharer, u0::AbstractVector, tspan, p=nothing; kwargs...) =
+DiscreteProblem(r::DiscreteResourceSharer, u0, tspan, p=nothing; kwargs...) =
     DiscreteProblem(dynamics(r), u0, tspan, p; kwargs...)
 
 """    trajectory(r::DiscreteResourceSharer, u0::AbstractVector, p, nsteps::Int; dt::Int = 1)
 
 Evolves the resouce sharer `r` for `nsteps` times with step size `dt`, initial condition `u0`, and parameters `p`.
 """
-function trajectory(r::DiscreteResourceSharer, u0::AbstractVector, p, T::Int; dt::Int= 1)
+function trajectory(r::DiscreteResourceSharer, u0, p, T::Int; dt::Int= 1)
   prob = DiscreteProblem(r, u0, (0, T), p)
   sol = solve(problem, FunctionMap(); dt = dt)
   return sol.xs

test/dwd_dynam.jl

@@ -42,7 +42,7 @@ add_wires!(d_big, Pair[
 
 @testset "ODE Problems" begin
   # Identity 
-  uf(u, x, p, t) = [x[1] - u[1]]
+  uf(u, x, p, t) = x - u
   rf(u, args...) = u
   mf = ContinuousMachine{Float64}(1,1,1, uf, rf)
 
@@ -53,7 +53,7 @@ add_wires!(d_big, Pair[
   @test eval_dynamics(m_id, [x0], [p0]) == [p0 - x0]
   @test readout(m_id, [x0]) == [x0]
 
-  # unfed parameter 
+  # compose 
   m12 = oapply(d12, [mf, mf])
 
   x0 = -1
@@ -62,6 +62,11 @@ add_wires!(d_big, Pair[
   @test eval_dynamics(m12, [x0, y0], [p0]) == [p0 - x0, x0 - y0]
   @test readout(m12, [x0,y0]) == [y0]
 
+  # test batch 
+  batch_size = 10
+  us = reshape(1:(2*batch_size), 2, batch_size)
+  xs = reshape(1:batch_size, 1, batch_size)
+  @test eval_dynamics(m12, us, xs) == (hcat(collect(1:batch_size) .- collect(us[1,:]), collect(us[1,:]) - collect(us[2,:])) |> transpose)
 
   # break and back together
   m = oapply(d_copymerge, Dict(:f => mf))

test/trajectories.jl

@@ -17,7 +17,7 @@ approx_equal(u, v) = abs(maximum(u - v)) < .1
 dx(x) = [1 - x[1]^2, 2*x[1]-x[2]]
 dy(y) = [1 - y[1]^2]
 
-r = ContinuousResourceSharer{Real}(2, (u,p,t) -> dx(u))
+r = ContinuousResourceSharer{Float64}(2, (u,p,t) -> dx(u))
 
 u0 = [-1.0, -2.0]
 tspan = (0.0, 100.0)
@@ -43,7 +43,7 @@ t = solve(dds, FunctionMap())
 
 
 
-dr = DiscreteResourceSharer{Real}(1, (u,p,t) -> dy(u))
+dr = DiscreteResourceSharer{Float64}(1, (u,p,t) -> dy(u))
 u0 = [1.0]
 dds = DiscreteProblem(dr, u0, tspan, nothing)
 t = solve(dds, FunctionMap())
@@ -109,9 +109,9 @@ dotr(u,p,t) = p[1]*u
 dotrf(u,p,t) = [-p[2]*u[1]*u[2], p[3]*u[1]*u[2]]
 dotf(u,p,t) = -p[4]*u
 
-r = ContinuousResourceSharer{Real}(1, dotr)
-rf_pred = ContinuousResourceSharer{Real}(2, dotrf)
-f = ContinuousResourceSharer{Real}(1, dotf)
+r = ContinuousResourceSharer{Float64}(1, dotr)
+rf_pred = ContinuousResourceSharer{Float64}(2, dotrf)
+f = ContinuousResourceSharer{Float64}(1, dotf)
 
 
 rf_pattern = UWD(0)
@@ -142,8 +142,8 @@ end)
 dotr(u, x, p, t) = [p[1]*u[1] - p[2]*u[1]*x[1]]
 dotf(u, x, p, t) = [p[3]*u[1]*x[1] - p[4]*u[1]]
 
-rmachine = ContinuousMachine{Real}(1,1,1, dotr, (r,p,t) -> r)
-fmachine = ContinuousMachine{Real}(1,1,1, dotf, (f,p,t) -> f)
+rmachine = ContinuousMachine{Float64}(1,1,1, dotr, (r,p,t) -> r)
+fmachine = ContinuousMachine{Float64}(1,1,1, dotf, (f,p,t) -> f)
 
 rf_pattern = WiringDiagram([],[])
 boxr = add_box!(rf_pattern, Box(nothing, [nothing], [nothing]))
@@ -180,9 +180,9 @@ dotfish(f, x, p, t) = [p[1]*f[1] - p[2]*x[1]*f[1]]
 dotFISH(F, x, p, t) = [p[3]*x[1]*F[1] - p[4]*F[1] - p[5]*x[2]*F[1]]
 dotsharks(s, x, p, t) = [-p[7]*s[1] + p[6]*s[1]*x[1]]
 
-fish   = ContinuousMachine{Real}(1,1,1, dotfish,   (f,p,t) ->f)
-FISH   = ContinuousMachine{Real}(2,1,2, dotFISH,   (F,p,t)->[F[1], F[1]])
-sharks = ContinuousMachine{Real}(1,1,1, dotsharks, (s,p,t)->s)
+fish   = ContinuousMachine{Float64}(1,1,1, dotfish,   (f,p,t) ->f)
+FISH   = ContinuousMachine{Float64}(2,1,2, dotFISH,   (F,p,t)->[F[1], F[1]])
+sharks = ContinuousMachine{Float64}(1,1,1, dotsharks, (s,p,t)->s)
 
 ocean_pat = WiringDiagram([], [])
 boxf = add_box!(ocean_pat, Box(nothing, [nothing], [nothing]))