/
rodesystems.jl
142 lines (124 loc) · 4.76 KB
/
rodesystems.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
# This file includes RODESystems
import DifferentialEquations: RandomEM, RODEProblem
import UUIDs: uuid4
"""
@def_rode_system ex
where `ex` is the expression to define to define a new AbstractRODESystem component type. The usage is as follows:
```julia
@def_rode_system mutable struct MyRODESystem{T1,T2,T3,...,TN,OP,RH,RO,ST,IP,OP} <: AbstractRODESystem
param1::T1 = param1_default # optional field
param2::T2 = param2_default # optional field
param3::T3 = param3_default # optional field
⋮
paramN::TN = paramN_default # optional field
righthandside::RH = righthandside_function # mandatory field
readout::RO = readout_function # mandatory field
state::ST = state_default # mandatory field
input::IP = input_default # mandatory field
output::OP = output_default # mandatory field
end
```
Here, `MyRODESystem` has `N` parameters. `MyRODESystem` is represented by the `righthandside` and `readout` function. `state`, `input` and `output` is the initial state, input port and output port of `MyRODESystem`.
!!! warning
`righthandside` must have the signature
```julia
function righthandside((dx, x, u, t, W, args...; kwargs...)
dx .= .... # update dx
end
```
and `readout` must have the signature
```julia
function readout(x, u, t)
y = ...
return y
end
```
!!! warning
New RODE system must be a subtype of `AbstractRODESystem` to function properly.
# Example
```julia
julia> @def_rode_system mutable struct MySystem{RH, RO, IP, OP} <: AbstractRODESystem
A::Matrix{Float64} = [2. 0.; 0 -2]
righthandside::RH = (dx, x, u, t, W) -> (dx .= A * x * W)
readout::RO = (x, u, t) -> x
state::Vector{Float64} = rand(2)
input::IP = nothing
output::OP = Outport(2)
end
julia> ds = MySystem();
```
"""
macro def_rode_system(ex)
checksyntax(ex, :AbstractRODESystem)
appendcommonex!(ex)
foreach(nex -> appendex!(ex, nex), [
:( alg::$ALG_TYPE_SYMBOL = Causal.RandomEM() ),
:( integrator::$INTEGRATOR_TYPE_SYMBOL = Causal.construct_integrator(
Causal.RODEProblem, input, righthandside, state, t, modelargs, solverargs; alg=alg, modelkwargs=modelkwargs,
solverkwargs = :dt in keys(solverkwargs) ? solverkwargs :
(; zip((keys(solverkwargs)..., :dt), (values(solverkwargs)..., 0.01))...),
numtaps=3) )
])
quote
Base.@kwdef $ex
end |> esc
end
##### Define RODE sytem library
"""
$TYPEDEF
Constructs a generic RODE system
# Fields
$TYPEDFIELDS
"""
@def_rode_system mutable struct RODESystem{RH,
RO,
ST <: AbstractVector{<:Real},
IP <: Union{<:Inport, <:Nothing},
OP <: Union{<:Outport,<:Nothing}} <: AbstractRODESystem
"Right-hand-side function"
righthandside::RH
"Readout function"
readout::RO
"State"
state::ST
"Input. Expected to be an `Inport` or `Nothing`"
input::IP
"Output port"
output::OP
end
"""
$TYPEDEF
Constructs a `MultiplicativeNoiseLinearSystem` with the dynamics
```math
\\begin{array}{l}
\\dot{x} = A x W
\\end{array}
where `W` is the noise process.
```
# Fields
$TYPEDFIELDS
"""
@def_rode_system mutable struct MultiplicativeNoiseLinearSystem{T1 <: AbstractMatrix{<:Real},
RH,
RO,
ST <: AbstractVector{<:Real},
IP <: Union{<:Inport, <:Nothing},
OP <: Union{<:Outport,<:Nothing}} <: AbstractRODESystem
"A"
A::T1 = [2. 0.; 0 -2]
"Right-hand-side function"
righthandside::RH = (dx, x, u, t, W) -> (dx .= A * x * W)
"Readout function"
readout::RO = (x, u, t) -> x
"State"
state::ST = rand(2)
"Input. Expected to be an `Inport` or `Nothing`"
input::IP = nothing
"Output port"
output::OP = Outport(2)
end
##### Pretty printing
show(io::IO, ds::RODESystem) = print(io, "RODESystem(righthandside:$(ds.righthandside), ",
"readout:$(ds.readout), state:$(ds.state), t:$(ds.t), input:$(ds.input), output:$(ds.output))")
show(io::IO, ds::MultiplicativeNoiseLinearSystem) = print(io,
"MultiplicativeNoiseLinearSystem(A:$(ds.A), state:$(ds.state), t:$(ds.t), input:$(ds.input), output:$(ds.output))")