SDESystem

Construction of SDESystems

A SDESystem is represented by the state function

\[ dx = f(x, u, t) dt + h(x, u, t)dW \]

where $t$ is the time, $x \in R^n$ is the value of state, $u \in R^p$ is the value of the input. $W$ is the Wiener process of the system. The output function is defined by

\[ y = g(x, u, t)\]

where $y$ is the value of output at time $t$.

As an example consider a system with the following stochastic differential equation

\[ \begin{array}{l} dx = -x dt - x dW \end{array}\]

and the following output equation

\[y = x\]

The state function statefunc and the output function outputfunc is defined as follows.

julia> using Causal # hide

julia> f(dx, x, u, t) = (dx[1] = -x[1])
f (generic function with 1 method)

julia> h(dx, x, u, t) = (dx[1] = -x[1])
h (generic function with 1 method)

The state function statefunc is the tuple of drift and diffusion functions

julia> statefunc = (f, h)
(Main.ex-sde_system_ex.f, Main.ex-sde_system_ex.h)

The output function outputfunc is defined as,

julia> g(x, u, t) = x
g (generic function with 1 method)

Note that the in drift function f and diffusion function g, the vector dx is mutated while in the output function g no mutation is done, but the output value is generated instead.

From the definition of drift function f and the diffusion function g, it is seen that the system does not have any input, that is, the input of the system is nothing. Since all the state variables are taken as outputs, the system needs an output bus of length 1. Thus,

julia> input = nothing

julia> output = Outport(1)
1-element Outport{Outpin{Float64}}:
 Outpin(eltype:Float64, isbound:false)

At this point, we are ready to construct the system ds.

julia> ds = SDESystem(righthandside=statefunc, readout=g, state=[1.], input=input, output=output)
ERROR: UndefKeywordError: keyword argument drift not assigned

Basic Operation of SDESystems

The basic operation of a SDESystem is the same as those of other dynamical systems. When triggered from its trigger link, a SDESystem reads its time t from its trigger link, reads its input value from its input, solves its state equation, which is a stochastic differential equation, computes its output and writes its computed output to its output bus.

In this section, we continue with the system ds constructed in the previous section. To make ds drivable, we need to launch it.

julia> iport, trg, hnd = Inport(1), Outpin(), Inpin{Bool}()
(Inport(numpins:1, eltype:Inpin{Float64}), Outpin(eltype:Float64, isbound:false), Inpin(eltype:Bool, isbound:false))

julia> connect!(ds.output, iport)
ERROR: UndefVarError: ds not defined

julia> connect!(trg, ds.trigger)
ERROR: UndefVarError: ds not defined

julia> connect!(ds.handshake, hnd)
ERROR: UndefVarError: ds not defined

julia> task = launch(ds)
ERROR: UndefVarError: ds not defined

julia> task2 = @async while true
           all(take!(iport) .=== NaN) && break
           end
Task (failed) @0x00007f1fd69d6710
MethodError: no method matching take!(::Missing)
Closest candidates are:
  take!(!Matched::Distributed.WorkerPool) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/Distributed/src/workerpool.jl:136
  take!(!Matched::IOStream) at iostream.jl:419
  take!(!Matched::Base.GenericIOBuffer{Array{UInt8,1}}) at iobuffer.jl:385
  ...
take!(::Inpin{Float64}) at /home/runner/work/Causal.jl/Causal.jl/src/connections/pin.jl:111
_broadcast_getindex_evalf at ./broadcast.jl:648 [inlined]
_broadcast_getindex at ./broadcast.jl:621 [inlined]
getindex at ./broadcast.jl:575 [inlined]
macro expansion at ./broadcast.jl:932 [inlined]
macro expansion at ./simdloop.jl:77 [inlined]
copyto! at ./broadcast.jl:931 [inlined]
copyto! at ./broadcast.jl:886 [inlined]
copy at ./broadcast.jl:862 [inlined]
materialize at ./broadcast.jl:837 [inlined]
take!(::Inport{Inpin{Float64}}) at /home/runner/work/Causal.jl/Causal.jl/src/connections/port.jl:186
macro expansion at ./none:2 [inlined]
(::Main.ex-sde_system_ex.var"#1#2")() at ./task.jl:356

When launched, ds can be driven. For this, either of the syntax put!(ds.trigger, t) or drive(ds, t) can be used.

julia> put!(trg, 1.)
ERROR: MethodError: no method matching iterate(::Missing)
Closest candidates are:
  iterate(!Matched::Combinatorics.Combinations) at /home/runner/.julia/packages/Combinatorics/Udg6X/src/combinations.jl:13
  iterate(!Matched::Combinatorics.Combinations, !Matched::Any) at /home/runner/.julia/packages/Combinatorics/Udg6X/src/combinations.jl:13
  iterate(!Matched::SimpleTraits.GenerateTypeVars{:upcase}) at /home/runner/.julia/packages/SimpleTraits/1wYi7/src/SimpleTraits.jl:445
  ...

After this command, ds reads its time t from its trigger link, solves its state function and computes its output. The calculated output value is written to the buffer of output. To signal that, the step is takes with success, ds writes true to its handshake link. To further drive ds, this handshake link must be read. For this either of the syntax, take!(ds.handshake) or approve!(ds) can be used

julia> hnd.link
missing

julia> take!(hnd)
ERROR: MethodError: no method matching take!(::Missing)
Closest candidates are:
  take!(!Matched::Distributed.WorkerPool) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/Distributed/src/workerpool.jl:136
  take!(!Matched::IOStream) at iostream.jl:419
  take!(!Matched::Base.GenericIOBuffer{Array{UInt8,1}}) at iobuffer.jl:385
  ...

At this point, we can further drive ds.

julia> for t in 2. : 10.
           put!(trg, t)
           take!(hnd)
       end
ERROR: MethodError: no method matching iterate(::Missing)
Closest candidates are:
  iterate(!Matched::Combinatorics.Combinations) at /home/runner/.julia/packages/Combinatorics/Udg6X/src/combinations.jl:13
  iterate(!Matched::Combinatorics.Combinations, !Matched::Any) at /home/runner/.julia/packages/Combinatorics/Udg6X/src/combinations.jl:13
  iterate(!Matched::SimpleTraits.GenerateTypeVars{:upcase}) at /home/runner/.julia/packages/SimpleTraits/1wYi7/src/SimpleTraits.jl:445
  ...

Note that during the evolution, the output of ds is written into the buffers of output bus.

julia> iport[1].link.buffer
ERROR: type Missing has no field buffer
Warning

The values of the output is written into buffers if the output of the systems is not nothing.

When we launched ds, we constructed a task whose state is running which implies that the ds can be drivable. As long as this task is running, ds can be drivable.

Warning

The state of the task is different from running in case an exception is thrown.

To terminate the task securely, we need to terminate ds securely. To do that, can use terminate!(ds).

julia> put!(trg, NaN)
ERROR: MethodError: no method matching iterate(::Missing)
Closest candidates are:
  iterate(!Matched::Combinatorics.Combinations) at /home/runner/.julia/packages/Combinatorics/Udg6X/src/combinations.jl:13
  iterate(!Matched::Combinatorics.Combinations, !Matched::Any) at /home/runner/.julia/packages/Combinatorics/Udg6X/src/combinations.jl:13
  iterate(!Matched::SimpleTraits.GenerateTypeVars{:upcase}) at /home/runner/.julia/packages/SimpleTraits/1wYi7/src/SimpleTraits.jl:445
  ...

julia> put!(ds.output, [NaN])
ERROR: UndefVarError: ds not defined

Note that the task is terminated without a hassle.

julia> task
ERROR: UndefVarError: task not defined

julia> task2
Task (failed) @0x00007f1fd69d6710
MethodError: no method matching take!(::Missing)
Closest candidates are:
  take!(!Matched::Distributed.WorkerPool) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/Distributed/src/workerpool.jl:136
  take!(!Matched::IOStream) at iostream.jl:419
  take!(!Matched::Base.GenericIOBuffer{Array{UInt8,1}}) at iobuffer.jl:385
  ...
take!(::Inpin{Float64}) at /home/runner/work/Causal.jl/Causal.jl/src/connections/pin.jl:111
_broadcast_getindex_evalf at ./broadcast.jl:648 [inlined]
_broadcast_getindex at ./broadcast.jl:621 [inlined]
getindex at ./broadcast.jl:575 [inlined]
macro expansion at ./broadcast.jl:932 [inlined]
macro expansion at ./simdloop.jl:77 [inlined]
copyto! at ./broadcast.jl:931 [inlined]
copyto! at ./broadcast.jl:886 [inlined]
copy at ./broadcast.jl:862 [inlined]
materialize at ./broadcast.jl:837 [inlined]
take!(::Inport{Inpin{Float64}}) at /home/runner/work/Causal.jl/Causal.jl/src/connections/port.jl:186
macro expansion at ./none:2 [inlined]
(::Main.ex-sde_system_ex.var"#1#2")() at ./task.jl:356

Full API

Causal.@def_sde_systemMacro
@def_sde_system ex

where ex is the expression to define to define a new AbstractSDESystem component type. The usage is as follows:

@def_sde_system mutable struct MySDESystem{T1,T2,T3,...,TN,OP,RH,RO,ST,IP,OP} <: AbstractSDESystem
    param1::T1 = param1_default                 # optional field 
    param2::T2 = param2_default                 # optional field 
    param3::T3 = param3_default                 # optional field
        ⋮
    paramN::TN = paramN_default                 # optional field 
    drift::DR = drift_function                  # mandatory field
    diffusion::DF = diffusion_function          # mandatory field
    readout::RO = readout_functtion             # mandatory field
    state::ST = state_default                   # mandatory field
    input::IP = input_default                   # mandatory field
    output::OP = output_default                 # mandatory field
end

Here, MySDESystem has N parameters. MySDESystem is represented by the drift, diffusion and readout function. state, input and output is the initial state, input port and output port of MySDESystem.

Warning

drift must have the signature

function drift((dx, x, u, t, args...; kwargs...)
    dx .= .... # update dx
end

and diffusion must have the signature

function diffusion((dx, x, u, t, args...; kwargs...)
    dx .= .... # update dx
end

and readout must have the signature

function readout(x, u, t)
    y = ...
    return y
end
Warning

New SDE system must be a subtype of AbstractSDESystem to function properly.

Example

julia> @def_sde_system mutable struct MySDESystem{DR, DF, RO, IP, OP} <: AbstractSDESystem
           η::Float64 = 1.
           drift::DR = (dx, x, u, t) -> (dx .= x)
           diffusion::DF = (dx, x, u, t, η=η) -> (dx .= η)
           readout::RO = (x, u, t) -> x 
           state::Vector{Float64} = rand(2) 
           input::IP = nothing 
           output::OP = Outport(2)
       end

julia> ds = MySDESystem();
source
Causal.SDESystemType
mutable struct SDESystem{DR, DF, RO, ST<:(AbstractArray{var"#s301",1} where var"#s301"<:Real), IP<:(Union{var"#s300", var"#s299"} where var"#s299"<:Nothing where var"#s300"<:Inport), OP<:(Union{var"#s298", var"#s297"} where var"#s297"<:Nothing where var"#s298"<:Outport), var"253", var"254", var"255", Symbol, var"256", Float64, var"257", var"258", var"259", var"260", var"261", var"262"} <: AbstractSDESystem

Constructs a SDE system.

Fields

  • drift::Any

    Drift function

  • diffusion::Any

    diffusion function

  • readout::Any

    Readout function

  • state::AbstractArray{var"#s301",1} where var"#s301"<:Real

    State

  • input::Union{var"#s300", var"#s299"} where var"#s299"<:Nothing where var"#s300"<:Inport

    Input. Expected to be an Inport or Nothing

  • output::Union{var"#s298", var"#s297"} where var"#s297"<:Nothing where var"#s298"<:Outport

    Output port

  • trigger::Any

  • handshake::Any

  • callbacks::Any

  • name::Any

  • id::Any

  • t::Any

  • modelargs::Any

  • modelkwargs::Any

  • solverargs::Any

  • solverkwargs::Any

  • alg::Any

  • integrator::Any

source
Causal.NoisyLorenzSystemType
mutable struct NoisyLorenzSystem{T1<:Real, T2<:Real, T3<:Real, T4<:Real, T5<:Real, DR, DF, RO, ST<:(AbstractArray{var"#s301",1} where var"#s301"<:Real), IP<:(Union{var"#s300", var"#s299"} where var"#s299"<:Nothing where var"#s300"<:Inport), OP<:(Union{var"#s298", var"#s297"} where var"#s297"<:Nothing where var"#s298"<:Outport), var"253", var"254", var"255", Symbol, var"256", Float64, var"257", var"258", var"259", var"260", var"261", var"262"} <: AbstractSDESystem

Constructs a noisy Lorenz system

Fields

  • σ::Real

    σ

  • β::Real

    β

  • ρ::Real

    ρ

  • η::Real

    η

  • γ::Real

    γ

  • drift::Any

    Drift function

  • diffusion::Any

    Diffusion function

  • readout::Any

    Readout function

  • state::AbstractArray{var"#s301",1} where var"#s301"<:Real

    State

  • input::Union{var"#s300", var"#s299"} where var"#s299"<:Nothing where var"#s300"<:Inport

    Input. Expected to be an Inport or Nothing

  • output::Union{var"#s298", var"#s297"} where var"#s297"<:Nothing where var"#s298"<:Outport

    Output port

  • trigger::Any

  • handshake::Any

  • callbacks::Any

  • name::Any

  • id::Any

  • t::Any

  • modelargs::Any

  • modelkwargs::Any

  • solverargs::Any

  • solverkwargs::Any

  • alg::Any

  • integrator::Any

source
Causal.ForcedNoisyLorenzSystemType
mutable struct ForcedNoisyLorenzSystem{T1<:Real, T2<:Real, T3<:Real, T4<:Real, CM<:(AbstractArray{var"#s301",2} where var"#s301"<:Real), T5<:Real, DR, DF, RO, ST<:(AbstractArray{var"#s300",1} where var"#s300"<:Real), IP<:(Union{var"#s299", var"#s298"} where var"#s298"<:Nothing where var"#s299"<:Inport), OP<:(Union{var"#s297", var"#s296"} where var"#s296"<:Nothing where var"#s297"<:Outport), var"253", var"254", var"255", Symbol, var"256", Float64, var"257", var"258", var"259", var"260", var"261", var"262"} <: AbstractSDESystem

Constructs a noisy Lorenz system

Fields

  • σ::Real

    σ

  • β::Real

    β

  • ρ::Real

    ρ

  • η::Real

    η

  • cplmat::AbstractArray{var"#s301",2} where var"#s301"<:Real

    Input coupling matrix. Expected to be a diagonal matrix.

  • γ::Real

    γ

  • drift::Any

    Drift function

  • diffusion::Any

    Diffusion function

  • readout::Any

    Readout function

  • state::AbstractArray{var"#s300",1} where var"#s300"<:Real

    State

  • input::Union{var"#s299", var"#s298"} where var"#s298"<:Nothing where var"#s299"<:Inport

    Input. Expected to be an Inport or Nothing

  • output::Union{var"#s297", var"#s296"} where var"#s296"<:Nothing where var"#s297"<:Outport

    Output port

  • trigger::Any

  • handshake::Any

  • callbacks::Any

  • name::Any

  • id::Any

  • t::Any

  • modelargs::Any

  • modelkwargs::Any

  • solverargs::Any

  • solverkwargs::Any

  • alg::Any

  • integrator::Any

source