From 5d7ee8c35a905e37dd11847a78dac39674737558 Mon Sep 17 00:00:00 2001 From: thorek1 Date: Tue, 3 Oct 2023 16:43:18 +0200 Subject: [PATCH] first steps with filtering package --- Project.toml | 4 ++ test/figure_out_filtering.jl | 103 +++++++++++++++++++++++++++++++++++ 2 files changed, 107 insertions(+) create mode 100644 test/figure_out_filtering.jl diff --git a/Project.toml b/Project.toml index 3c8a71a0..3785515b 100644 --- a/Project.toml +++ b/Project.toml @@ -6,9 +6,11 @@ version = "0.1.29" [deps] AbstractDifferentiation = "c29ec348-61ec-40c8-8164-b8c60e9d9f3d" AxisKeys = "94b1ba4f-4ee9-5380-92f1-94cde586c3c5" +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" BlockTriangularForm = "adeb47b7-70bf-415a-bb24-c358563e873a" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" DynarePreprocessor_jll = "23afba7c-24e5-5ee2-bc2c-b42e07f0492a" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" @@ -18,6 +20,7 @@ Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearOperators = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125" +LowLevelParticleFilters = "d9d29d28-c116-5dba-9239-57a5fe23875b" MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" @@ -29,6 +32,7 @@ RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" SpeedMapping = "f1835b91-879b-4a3f-a438-e4baacf14412" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Subscripts = "2b7f82d5-8785-4f63-971e-f18ddbeb808e" SymPyPythonCall = "bc8888f7-b21e-4b7c-a06a-5d9c9496438c" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" diff --git a/test/figure_out_filtering.jl b/test/figure_out_filtering.jl new file mode 100644 index 00000000..1ef03cea --- /dev/null +++ b/test/figure_out_filtering.jl @@ -0,0 +1,103 @@ +using MacroModelling, LowLevelParticleFilters, Distributions +import LinearAlgebra as ℒ + +include("models/FS2000.jl") + +simulations = get_simulation(m) + +# (x,u,p,t) = (state, input, parameters, time) +# state_update₁ = function(state::Vector{Float64}, shock::Vector{Float64}) sol_mat * [state[𝓂.timings.past_not_future_and_mixed_idx]; shock] end + +# dynamics(state, input, parameters, time) = m.solution.perturbation.first_order.solution_matrix * [state[m.timings.past_not_future_and_mixed_idx]; input] +dynamics(state, input, parameters, time) = m.solution.perturbation.first_order.state_update(state,input) +measurement(state, input, parameters, time) = m.solution.perturbation.first_order.state_update(state,input)#ℒ.I * state + +covariance = get_covariance(m) |> collect +# covariance = sparse(covariance) +# droptol!(covariance,eps()) +# covariance = sparse(Symmetric(covariance)) +# isposdef(covariance) +ukf = UnscentedKalmanFilter(dynamics, + measurement, + eye(m.timings.nVars), + eye(m.timings.nVars), + MvNormal(diag(covariance)), + nu = m.timings.nExo, + ny = m.timings.nVars) + + +x,u,y = LowLevelParticleFilters.simulate(ukf,10,MvNormal(ones(m.timings.nExo))) # Simuate trajectory using the model in the filter + + +using LinearAlgebra +covariance +isposdef((covariance+covariance')/2) +factorize(Symmetric(covariance)) + + + + +using Statistics, LinearAlgebra, StaticArrays +m = randn(3) +S = randn(3,3) +S = S'S +xs = LowLevelParticleFilters.sigmapoints(m, S) +X = reduce(hcat, xs) +@test vec(mean(X, dims=2)) ≈ m +@test Statistics.cov(X, dims=2) ≈ S + +m = [1,2] +S = [3. 1; 1 4] +xs = LowLevelParticleFilters.sigmapoints(m, S) +X = reduce(hcat, xs) +# @test vec(mean(X, dims=2)) ≈ m +# @test cov(X, dims=2) ≈ S + + + + +eye(n) = Matrix{Float64}(I,n,n) +mvnormal(d::Int, σ::Real) = MvNormal(LinearAlgebra.Diagonal(fill(float(σ) ^ 2, d))) +mvnormal(μ::AbstractVector{<:Real}, σ::Real) = MvNormal(μ, float(σ) ^ 2 * I) + +nx = 2 # Dinemsion of state +nu = 2 # Dinemsion of input +ny = 2 # Dinemsion of measurements + +d0 = mvnormal(randn(nx),2.0) # Initial state Distribution +du = mvnormal(2,1) # Control input distribution + +# Define random linenar state-space system +Tr = randn(nx,nx) +A = SA[0.99 0.1; 0 0.2] +B = @SMatrix randn(nx,nu) +C = SMatrix{ny,ny}(eye(ny)) +# C = SMatrix{p,n}([1 1]) + +dynamics(x,u,p,t) = A*x .+ B*u +measurement(x,u,p,t) = C*x + +# sss = SS(m,derivatives=false) |> collect +# I(m.timings.nVars)[m.timings.past_not_future_and_mixed_idx,:] * [sss...,0,0] + +T = 200 # Number of time steps +kf = KalmanFilter(A, B, C, 0, eye(nx), eye(ny), d0) +ukf = UnscentedKalmanFilter(dynamics, measurement, eye(nx), eye(ny), d0; ny, nu) +x,u,y = LowLevelParticleFilters.simulate(kf,T,du) # Simuate trajectory using the model in the filter + +using SparseArrays, LinearOperators, BenchmarkTools +n = 20 +p = 10 +x = sprand(p,n^3,.1) +kk = sprand(n,n,.2) +kkLinOp = kk' |> LinearOperator + +@benchmark kron(kron(kron(kkLinOp, kkLinOp), kkLinOp), I(p)) * vec(x) +@benchmark x * kron(kron(kk,kk),kk) + +sparse(reshape(kron(kron(kron(kk', kk'), kk'), I(p)) * vec(x),p,n^3)) + +findnz(sparse(reshape(kron(I(p),kron(kron(kk,kk),kk)) * vec(x),p,n^3))) + +x * kron(kron(kk,kk),kk) +findnz(x * kron(kron(kk,kk),kk))