Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
3 changes: 0 additions & 3 deletions .JuliaFormatter.toml

This file was deleted.

19 changes: 19 additions & 0 deletions .github/workflows/FormatCheck.yml
Original file line number Diff line number Diff line change
@@ -0,0 +1,19 @@
name: format-check

on:
push:
branches:
- 'master'
- 'main'
- 'release-'
tags: '*'
pull_request:

jobs:
runic:
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v4
- uses: fredrikekre/runic-action@v1
with:
version: '1'
4 changes: 2 additions & 2 deletions docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -17,10 +17,10 @@ makedocs(;
"Examples" => [
"examples/DiffEqFlux.md",
"examples/adaptive_control.md",
"examples/ODE_jac.md",
"examples/coulomb_control.md",
"examples/ODE_jac.md"
],
"API" => "api.md"
"API" => "api.md",
],
repo = GitHub("SciML/ComponentArrays.jl"),
sitename = "ComponentArrays.jl",
Expand Down
26 changes: 16 additions & 10 deletions examples/DiffEqFlux_example.jl
Original file line number Diff line number Diff line change
Expand Up @@ -23,12 +23,12 @@ tspan2 = (0.0f0, 25.0f0)
# Make truth data
function trueODEfunc(du, u, p, t)
true_A = [-0.1 2.0; -2.0 -0.1]
du .= ((u .^ 3)'true_A)'
return du .= ((u .^ 3)'true_A)'
end

t = range(tspan[1], tspan[2], length = datasize)
# t = Float32.(vcat(range(0.0, 0.9, length=10), 10 .^ range(log10(tspan[1]+1), log10(tspan[2]), length=datasize-10)))
t = Float32.([0; 10 .^ range(log10(tspan[1] + 0.01), log10(tspan[2]), length = datasize-1)])
t = Float32.([0; 10 .^ range(log10(tspan[1] + 0.01), log10(tspan[2]), length = datasize - 1)])
prob = ODEProblem(trueODEfunc, u0, tspan2)
ode_sol = solve(prob, Tsit5())
ode_data = Array(ode_sol(t)) + MeasurementNoise(0.1)
Expand All @@ -39,7 +39,7 @@ neural_layer(in, out) = ComponentArray{Float32}(W = glorot_uniform(out, in), b =
# Dense neural layer function
dense(layer, activation = identity) = u -> activation.(layer.W * u + layer.b)

# Neural ODE function
# Neural ODE function
dudt(u, p, t) = u .^ 3 |> dense(p.L1, σ) |> dense(p.L2)

prob = ODEProblem(dudt, u0, tspan2)
Expand All @@ -56,7 +56,7 @@ full_sol(θ) = solve(prob, Tsit5(), u0 = θ.u, p = θ.p)
function loss_n_ode(θ)
pred = predict_n_ode(θ)

loss = sum(abs2, ode_data .- pred)/datasize + 0.1*(sum(abs, θ.p)/length(θ.p))
loss = sum(abs2, ode_data .- pred) / datasize + 0.1 * (sum(abs, θ.p) / length(θ.p))
return loss, pred
end
loss_n_ode(θ)
Expand All @@ -83,12 +83,18 @@ cb = function (θ, loss, pred; doplot = false)
plot!(pl_3, ode_sol, vars = (1, 2), label = "truth")
scatter!(pl_3, pred[1, :], pred[2, :], label = "predicted data")
scatter!(pl_3, ode_data[1, :], ode_data[2, :], label = "measured data")
plot!(pl_3, hcat(pred[1, :], ode_data[1, :])', hcat(pred[2, :], ode_data[2, :])',
label = false, color = :lightgray, legend = :bottomright)

display(plot(
plot(pl_1, pl_2, layout = (2, 1), size = (400, 500)), pl_3, layout = (1, 2), size = (
950, 500)))
plot!(
pl_3, hcat(pred[1, :], ode_data[1, :])', hcat(pred[2, :], ode_data[2, :])',
label = false, color = :lightgray, legend = :bottomright
)

display(
plot(
plot(pl_1, pl_2, layout = (2, 1), size = (400, 500)), pl_3, layout = (1, 2), size = (
950, 500,
)
)
)
# frame(anim)
return false
end
Expand Down
18 changes: 9 additions & 9 deletions examples/ODE_example.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,13 +9,13 @@ function lorenz!(D, u, p, t; f = 0.0)
@unpack σ, ρ, β = p
@unpack x, y, z = u

D.x = σ*(y - x)
D.y = x*(ρ - z) - y - f
D.z = x*y - β*z
D.x = σ * (y - x)
D.y = x * (ρ - z) - y - f
D.z = x * y - β * z
return nothing
end

lorenz_p = (σ = 10.0, ρ = 28.0, β = 8/3)
lorenz_p = (σ = 10.0, ρ = 28.0, β = 8 / 3)
lorenz_ic = ComponentArray(x = 0.0, y = 0.0, z = 0.0)
lorenz_prob = ODEProblem(lorenz!, lorenz_ic, tspan, lorenz_p)

Expand All @@ -24,12 +24,12 @@ function lotka!(D, u, p, t; f = 0.0)
@unpack α, β, γ, δ = p
@unpack x, y = u

D.x = α*x - β*x*y + f
D.y = -γ*y + δ*x*y
D.x = α * x - β * x * y + f
D.y = -γ * y + δ * x * y
return nothing
end

lotka_p = (α = 2/3, β = 4/3, γ = 1.0, δ = 1.0)
lotka_p = (α = 2 / 3, β = 4 / 3, γ = 1.0, δ = 1.0)
lotka_ic = ComponentArray(x = 1.0, y = 1.0)
lotka_prob = ODEProblem(lotka!, lotka_ic, tspan, lotka_p)

Expand All @@ -38,8 +38,8 @@ function composed!(D, u, p, t)
c = p.c #coupling parameter
@unpack lorenz, lotka = u

lorenz!(D.lorenz, lorenz, p.lorenz, t, f = c*lotka.x)
lotka!(D.lotka, lotka, p.lotka, t, f = c*lorenz.x)
lorenz!(D.lorenz, lorenz, p.lorenz, t, f = c * lotka.x)
lotka!(D.lotka, lotka, p.lotka, t, f = c * lorenz.x)
return nothing
end

Expand Down
26 changes: 13 additions & 13 deletions examples/ODE_jac_example.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,9 +9,9 @@ function lorenz!(D, u, p, t; f = 0.0)
@unpack σ, ρ, β = p
@unpack x, y, z = u

D.x = σ*(y - x)
D.y = x*(ρ - z) - y - f
D.z = x*y - β*z
D.x = σ * (y - x)
D.y = x * (ρ - z) - y - f
D.z = x * y - β * z
return nothing
end
function lorenz_jac!(D, u, p, t)
Expand All @@ -31,7 +31,7 @@ function lorenz_jac!(D, u, p, t)
return nothing
end

lorenz_p = (σ = 10.0, ρ = 28.0, β = 8/3)
lorenz_p = (σ = 10.0, ρ = 28.0, β = 8 / 3)
lorenz_ic = ComponentArray(x = 0.0, y = 0.0, z = 0.0)
lorenz_fun = ODEFunction(lorenz!, jac = lorenz_jac!)
lorenz_prob = ODEProblem(lorenz_fun, lorenz_ic, tspan, lorenz_p)
Expand All @@ -41,23 +41,23 @@ function lotka!(D, u, p, t; f = 0.0)
@unpack α, β, γ, δ = p
@unpack x, y = u

D.x = α*x - β*x*y + f
D.y = -γ*y + δ*x*y
D.x = α * x - β * x * y + f
D.y = -γ * y + δ * x * y
return nothing
end
function lotka_jac!(D, u, p, t)
@unpack α, β, γ, δ = p
@unpack x, y = u

D[:x, :x] = α - β*y
D[:x, :y] = -β*x
D[:x, :x] = α - β * y
D[:x, :y] = -β * x

D[:y, :x] = δ*y
D[:y, :y] = -γ + δ*x
D[:y, :x] = δ * y
D[:y, :y] = -γ + δ * x
return nothing
end

lotka_p = (α = 2/3, β = 4/3, γ = 1.0, δ = 1.0)
lotka_p = (α = 2 / 3, β = 4 / 3, γ = 1.0, δ = 1.0)
lotka_ic = ComponentArray(x = 1.0, y = 1.0)
lotka_fun = ODEFunction(lotka!, jac = lotka_jac!)
lotka_prob = ODEProblem(lotka_fun, lotka_ic, tspan, lotka_p)
Expand All @@ -67,8 +67,8 @@ function composed!(D, u, p, t)
c = p.c #coupling parameter
@unpack lorenz, lotka = u

lorenz!(D.lorenz, lorenz, p.lorenz, t, f = c*lotka.x)
lotka!(D.lotka, lotka, p.lotka, t, f = c*lorenz.x)
lorenz!(D.lorenz, lorenz, p.lorenz, t, f = c * lotka.x)
lotka!(D.lotka, lotka, p.lotka, t, f = c * lorenz.x)
return nothing
end
function composed_jac!(D, u, p, t)
Expand Down
50 changes: 30 additions & 20 deletions examples/adaptive_control_example.jl
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,9 @@ maybe_apply(f, x, p, t) = f

function apply_inputs(func; kwargs...)
simfun(dx, x, p, t) = func(
dx, x, p, t; map(f->maybe_apply(f, x, p, t), (; kwargs...))...)
simfun(x, p, t) = func(x, p, t; map(f->maybe_apply(f, x, p, t), (; kwargs...))...)
dx, x, p, t; map(f -> maybe_apply(f, x, p, t), (; kwargs...))...
)
simfun(x, p, t) = func(x, p, t; map(f -> maybe_apply(f, x, p, t), (; kwargs...))...)
return simfun
end

Expand All @@ -27,7 +28,7 @@ SISO_simulator(P::TransferFunction) = SISO_simulator(ss(P))
function SISO_simulator(P::AbstractStateSpace)
@unpack A, B, C, D = P

if size(D)!=(1, 1)
if size(D) != (1, 1)
error("This is not a SISO system")
end

Expand All @@ -37,8 +38,8 @@ function SISO_simulator(P::AbstractStateSpace)
DD = D[1, 1]

return function sim!(dx, x, p, t; u = 0.0)
dx .= A*x + BB*u
return CC*x + DD*u
dx .= A * x + BB * u
return CC * x + DD * u
end
end

Expand All @@ -60,13 +61,13 @@ nominal_sim! = SISO_simulator(nominal_plant)

# To test robustness to uncertainty, we'll also include unmodeled dynamics with an entirely
# different structure than our nominal plant model.
unmodeled_dynamics = 229/(s^2 + 30s + 229)
unmodeled_dynamics = 229 / (s^2 + 30s + 229)
truth_plant = nominal_plant * unmodeled_dynamics
truth_sim! = SISO_simulator(truth_plant)

# We'll make a first-order sensor as well so we can add noise to our measurement
τ = 0.005
sensor_plant = 1 / (τ*s + 1)
sensor_plant = 1 / (τ * s + 1)
sensor_sim! = SISO_simulator(sensor_plant)

## Derivative functions
Expand All @@ -76,7 +77,7 @@ control(θ, w) = θ'w

# We'll use a simple gradient descent adaptation law
function adapt!(Dθ, θ, γ, t; e, w)
Dθ .= -γ*e*w
Dθ .= -γ * e * w
return nothing
end

Expand All @@ -98,7 +99,7 @@ function feedback_sys!(D, vars, p, t; ym, r, n)
return yp
end
# Now the full system takes in an input signal `r`, feeds it through the reference model,
# and feeds the output of the reference model `ym` and the input signal to `feedback_sys`.
# and feeds the output of the reference model `ym` and the input signal to `feedback_sys`.
function system!(D, vars, p, t; r = 0.0, n = 0.0)
@unpack reference_model, feedback_loop = vars

Expand All @@ -122,31 +123,33 @@ sensor_ic = zeros(1)
θ_est_ic = ComponentArray(θr = 0.0, θy = 0.0)

## Set up and run Simulation
function simulate(plant_fun, plant_ic;
function simulate(
plant_fun, plant_ic;
tspan = tspan,
input_signal = input_signal,
adapt_gain = 1.5,
noise_param = nothing,
deterministic_noise = 0.0)
deterministic_noise = 0.0
)
noise(D, vars, p, t) = (D.feedback_loop.sensor[1] = noise_param)

# Truth control parameters
θ_truth = (r = bm/bp, y = (ap-am)/bp)
θ_truth = (r = bm / bp, y = (ap - am) / bp)

# Initial conditions
ic = ComponentArray(
reference_model = ref_ic,
feedback_loop = (
parameter_estimates = θ_est_ic,
sensor = sensor_ic,
plant_model = plant_ic
plant_model = plant_ic,
)
)

# Model parameters
p = (
gamma = adapt_gain,
plant_fun = plant_fun
plant_fun = plant_fun,
)

sim_fun = apply_inputs(system!; r = input_signal, n = deterministic_noise)
Expand All @@ -172,9 +175,14 @@ function simulate(plant_fun, plant_ic;
)

# Parameter estimate tracking
bottom = plot(sol,
vars = Symbol.([
"feedback_loop.parameter_estimates.θr", "feedback_loop.parameter_estimates.θy"]))
bottom = plot(
sol,
vars = Symbol.(
[
"feedback_loop.parameter_estimates.θr", "feedback_loop.parameter_estimates.θy",
]
)
)
plot!(
bottom,
[tspan...], [θ_truth.r θ_truth.y; θ_truth.r θ_truth.y],
Expand All @@ -184,8 +192,10 @@ function simulate(plant_fun, plant_ic;
)

# Combine both plots
plot(top, bottom, layout = (2, 1), size = (800, 800))
return plot(top, bottom, layout = (2, 1), size = (800, 800))
end

simulate(truth_sim!, truth_ic; input_signal = 2.0,
deterministic_noise = (x, p, t)->0.5sin(16.1t), noise_param = nothing)
simulate(
truth_sim!, truth_ic; input_signal = 2.0,
deterministic_noise = (x, p, t) -> 0.5sin(16.1t), noise_param = nothing
)
Loading
Loading