From 006cd7f75354cdb1fbe6764c75ebf4b7b349b1e1 Mon Sep 17 00:00:00 2001 From: b-r-hamilton Date: Mon, 27 Mar 2023 20:22:06 -0400 Subject: [PATCH 1/8] new testset started, methods for x with two state var --- test/runtests.jl | 20 ++++++++++++++++++++ test/test_functions.jl | 11 +++++++++++ 2 files changed, 31 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index 7e14e38..ac3383c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -293,6 +293,26 @@ include("test_functions.jl") @test cost(x̃,problem) < 5e-2 # no noise in ob end + @testset "source water inversion: obs at one time, many surface regions, with circulation lag, TWO STATE VARIABLES" begin + + @dim YearCE "years Common Era" + @dim SurfaceRegion "surface location" + @dim InteriorLocation "interior location" + @dim StateVariable "state variable" + + #yr = u"yr" + nτ = 5 # how much of a lag is possible? + lags = (0:(nτ-1))yr + surfaceregions = [:NATL,:ANT,:SUBANT] + years = (1990:2000)yr + n = length(surfaceregions) + + M = source_water_matrix_with_lag(surfaceregions,lags,years) + + statevariables = [:θ, :d18O] + x = source_water_solution(surfaceregions,years, statevariables) + end + @testset "source water inversion: obs TIMESERIES, many surface regions, with circulation lag" begin @dim YearCE "years Common Era" diff --git a/test/test_functions.jl b/test/test_functions.jl index 5c018b7..1653212 100644 --- a/test/test_functions.jl +++ b/test/test_functions.jl @@ -106,6 +106,17 @@ function source_water_solution(surfaceregions,years) return x end +function source_water_solution(surfaceregions, years, statevar) + yr = u"yr" + K = u"K" + permil = u"permille" + m = length(years) + n = length(surfaceregions) + mat = cat(randn(m, n, 1)permil, randn(m, n, 1)K; dims = 3) + x = DimArray(mat, (Ti(years), SurfaceRegion(surfaceregions), StateVariable(statevar))) + return x +end + """ function source_water_matrix_vector_pair_with_lag(M) From c2600fc864af18fd92aa918a7fae6cd87e917744 Mon Sep 17 00:00:00 2001 From: b-r-hamilton Date: Tue, 28 Mar 2023 23:50:44 -0400 Subject: [PATCH 2/8] two state variable single obs test seems to be working --- src/BLUEs.jl | 15 ++++++++++++++- test/runtests.jl | 42 +++++++++++++++++++++++++++++++++++++++--- test/test_functions.jl | 2 +- 3 files changed, 54 insertions(+), 5 deletions(-) diff --git a/src/BLUEs.jl b/src/BLUEs.jl index 31b89b4..a25f608 100644 --- a/src/BLUEs.jl +++ b/src/BLUEs.jl @@ -2,6 +2,7 @@ module BLUEs using LinearAlgebra, Statistics, Unitful, UnitfulLinearAlgebra, Measurements using DimensionalData +using Revise export Estimate, DimEstimate, OverdeterminedProblem, UnderdeterminedProblem export solve, show, cost, datacost, controlcost @@ -559,6 +560,19 @@ function convolve(x::AbstractDimArray,E::AbstractDimArray) lags = first(dims(E)) return sum([E[ii,:] ⋅ x[Near(tnow-ll),:] for (ii,ll) in enumerate(lags)]) end +""" +function convolve(x::AbstractDimArray, E::AbstractDimArray, statevars::Vector{Symbol} + + This is a silly method - the info in statevars is a dimension in x + I can't figure out how to have one convolve method for a 2D AbstractDimArray + and another convolve method for a 3D AbstractDimArray + AbstractDimArray{Any, 3} does not work, even though I really want it to +""" +function convolve(x::AbstractDimArray,E::AbstractDimArray, statevars::Vector{Symbol}) + mat = [convolve(x[:,:,At(s)], E) for s in statevars] + #return DimArray(mat, statevars) + return mat +end function convolve(x::AbstractDimArray,E::AbstractDimArray,t::Number) lags = first(dims(E)) return sum([E[ii,:] ⋅ x[Near(t-ll),:] for (ii,ll) in enumerate(lags)]) @@ -626,5 +640,4 @@ function addcontrol!(x::AbstractDimArray,u) end return x end - end # module diff --git a/test/runtests.jl b/test/runtests.jl index ac3383c..13e33b6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -260,14 +260,15 @@ include("test_functions.jl") y = predictobs(convolve,x,M) ## invert for y for x̃ - # Given, M and y. Make first guess for x. + # Given, M and y. Make first guess for x.x # add adjustment # DimArray is good enough. This is an array, not necessarily a matrix. x₀ = DimArray(zeros(size(x))K,(Ti(years),last(dims(M)))) # probe to get E matrix. Use function convolve E = impulseresponse(convolve,x₀,M) - + + @test (E*UnitfulMatrix(vec(x₀)))[1] .== ustrip(convolve(x₀, M)) # Does E matrix work properly? ỹ = E*UnitfulMatrix(vec(x)) @test isapprox(y,getindexqty(ỹ,1)) @@ -280,7 +281,8 @@ include("test_functions.jl") # unitful scalars in UnitfulLinearAlgebra σₙ = 0.01 - σₓ = 100.0 + σₓ_d18O = 100.0 + σₓ_θ = 100.0 Cnn = UnitfulMatrix(Diagonal(fill(σₙ,length(y))),fill(unit.(y).^1,length(y)),fill(unit.(y).^-1,length(y)),exact=false) @@ -311,6 +313,40 @@ include("test_functions.jl") statevariables = [:θ, :d18O] x = source_water_solution(surfaceregions,years, statevariables) + + #convolve(x, M, s) gives us x variables propagated to location + #here is a function to linearly combine those two propagated value into a single value + coeffs = UnitfulMatrix(reshape(rand(2), (2,1)), [unit(1.0), K], [K*permil^-1]) + combine(y, coeffs) = UnitfulMatrix(transpose(y)) * coeffs + y = combine(convolve(x, M, statevariables), coeffs) + #let's make one function to give to impulseresponse + prop_and_combine(x, M, statevariables, coeffs) = getindexqty(combine(convolve(x, M, statevariables), coeffs), 1, 1) + @test getindexqty(y,1,1) .== prop_and_combine(x, M, statevariables, coeffs) + + x₀ = copy(x).*0 + y₀ = convolve(x₀, M, statevariables) #seems to work + E = impulseresponse(prop_and_combine,x₀,M, statevariables, coeffs) + #does it work? + @test getindexqty(E*UnitfulMatrix(vec(x₀)),1) == prop_and_combine(x₀, M, statevariables, coeffs) + + #E matrix tests from prior example + ỹ = E*UnitfulMatrix(vec(x)) + @test isapprox(getindexqty(y,1,1),getindexqty(ỹ,1)) + + #had to get a little crafty to make this work... + x̂ = E\UnitfulMatrix(parent(y), [permil], [unit(1.0)]) + @test isapprox(getindexqty(y,1,1),getindexqty(E*x̂,1,1)) + + σₙ = 0.01 + σₓ = 100.0 + + Cnn = UnitfulMatrix(Diagonal(fill(σₙ,length(y))),fill(unit(getindexqty(y, 1,1)).^1,length(y)),fill(unit(getindexqty(y,1,1)).^-1,length(y)),exact=false) + + Cxx = UnitfulMatrix(Diagonal(fill(σₓ,length(x₀))),unit.(x₀)[:],unit.(x₀)[:].^-1,exact=false) + problem = UnderdeterminedProblem(UnitfulMatrix([getindexqty(y,1,1)]), E, Cnn, Cxx, x₀) + x̃ = solve(problem) + @test within(getindexqty(y, 1,1), getindexqty(E*x̃.v, 1), 3σₙ) + @test cost(x̃,problem) < 5e-2 end @testset "source water inversion: obs TIMESERIES, many surface regions, with circulation lag" begin diff --git a/test/test_functions.jl b/test/test_functions.jl index 1653212..a1a66fb 100644 --- a/test/test_functions.jl +++ b/test/test_functions.jl @@ -112,7 +112,7 @@ function source_water_solution(surfaceregions, years, statevar) permil = u"permille" m = length(years) n = length(surfaceregions) - mat = cat(randn(m, n, 1)permil, randn(m, n, 1)K; dims = 3) + mat = cat(randn(m, n, 1)K, randn(m, n, 1)permil; dims = 3) x = DimArray(mat, (Ti(years), SurfaceRegion(surfaceregions), StateVariable(statevar))) return x end From ad1cae0b9ada08949f2c69a2b8be98950aa9a1bf Mon Sep 17 00:00:00 2001 From: b-r-hamilton Date: Wed, 29 Mar 2023 21:01:02 -0400 Subject: [PATCH 3/8] timeseries, two state variables working? --- src/BLUEs.jl | 33 +++++++++++++++++++----- test/runtests.jl | 65 +++++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 91 insertions(+), 7 deletions(-) diff --git a/src/BLUEs.jl b/src/BLUEs.jl index a25f608..3223fda 100644 --- a/src/BLUEs.jl +++ b/src/BLUEs.jl @@ -474,12 +474,12 @@ function impulseresponse(funk::Function,x,args...) y₀ = funk(x,args...) #Eunits = expectedunits(y₀,x) #Eu = Quantity.(zeros(length(y₀),length(x)),Eunits) - - Ep = zeros(length(y₀),length(x)) + Ep = zeros(length(y₀), length(x)) + for rr in eachindex(x) #u = zeros(length(x)).*unit.(x)[:] if length(x) > 1 - u = Quantity.(zeros(length(x)),unit.(vec(x))) + u = Quantity.(zeros(length(x)), unit.(vec(x))) else u = Quantity(0.0,unit(x)) end @@ -490,7 +490,9 @@ function impulseresponse(funk::Function,x,args...) y = funk(x₁,args...) #println(y) if y isa AbstractDimArray - Ep[:,rr] .= vec(parent((y - y₀)/Δu)) + #Ep[:,rr] .= vec(parent((y - y₀)/Δu)) + Ep[:, rr] .= ustrip.(vec(parent((y - y₀)/Δu))) + else tmp = ustrip.((y - y₀)/Δu) # getting ugly around here @@ -501,7 +503,8 @@ function impulseresponse(funk::Function,x,args...) end end end - + + # This function could use vcat to be cleaner (but maybe slower) # Note: Need to add ability to return sparse matrix #return E = UnitfulMatrix(Ep,unit.(vec(y₀)),unit.(vec(x))) @@ -562,7 +565,6 @@ function convolve(x::AbstractDimArray,E::AbstractDimArray) end """ function convolve(x::AbstractDimArray, E::AbstractDimArray, statevars::Vector{Symbol} - This is a silly method - the info in statevars is a dimension in x I can't figure out how to have one convolve method for a 2D AbstractDimArray and another convolve method for a 3D AbstractDimArray @@ -577,6 +579,25 @@ function convolve(x::AbstractDimArray,E::AbstractDimArray,t::Number) lags = first(dims(E)) return sum([E[ii,:] ⋅ x[Near(t-ll),:] for (ii,ll) in enumerate(lags)]) end +# two more silly convolve function +function convolve(x::AbstractDimArray, E::AbstractDimArray, t::Number, statevars::Vector{Symbol}) + mat = [convolve(x[:,:,At(s)], E, t) for s in statevars] + return mat +end +#don't handle the ndims(M) == 3 case here but I'll get back to it +function convolve(x::AbstractDimArray, M::AbstractDimArray, Tx::Union{Ti, Vector}, statevars::Vector{Symbol}) + if ndims(M) == 2 + return DimArray([convolve(x,M,Tx[tt], statevars) for (tt, yy) in enumerate(Tx)], Tx) + elseif ndims(M) == 3 + error("some code should go here") + else + error("M has wrong number of dims") + end + + +end + + function convolve(x::AbstractDimArray,M::AbstractDimArray,Tx::Union{Ti,Vector}) if ndims(M) == 2 return DimArray([convolve(x,M,Tx[tt]) for (tt,yy) in enumerate(Tx)],Tx) diff --git a/test/runtests.jl b/test/runtests.jl index 13e33b6..1043014 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -346,9 +346,72 @@ include("test_functions.jl") problem = UnderdeterminedProblem(UnitfulMatrix([getindexqty(y,1,1)]), E, Cnn, Cxx, x₀) x̃ = solve(problem) @test within(getindexqty(y, 1,1), getindexqty(E*x̃.v, 1), 3σₙ) - @test cost(x̃,problem) < 5e-2 + @test cost(x̃,problem) < 5e-2 end + @testset "source water inversion: obs TIMESERIES, many surface regions, with circulation lag, TWO STATE VARIABLES" begin + + @dim YearCE "years Common Era" + @dim SurfaceRegion "surface location" + @dim InteriorLocation "interior location" + @dim StateVariable "state variable" + + yr = u"yr" + nτ = 5 # how much of a lag is possible? + lags = (0:(nτ-1))yr + surfaceregions = [:NATL,:ANT,:SUBANT] + years = (1990:2000)yr + n = length(surfaceregions) + + M = source_water_matrix_with_lag(surfaceregions,lags,years) + + statevariables = [:θ, :d18O] + x = source_water_solution(surfaceregions,years, statevariables) + + Tx = first(dims(x)) + x₀ = copy(x).*0 + + y_ = convolve(x, M, Tx, statevariables) + coeffs = UnitfulMatrix(reshape(rand(2), (2,1)), [unit(1.0), K], [K*permil^-1]) + #function to linearly combine state variable to one variable + combine(y_, coeffs) = [getindexqty(UnitfulMatrix(transpose(y__))*coeffs,1,1) for y__ in y_] + #test! + combine(y_, coeffs) + + #one function to hand to impulseresponse + prop_and_combine(x, M, Tx, statevariables, coeffs) = combine(convolve(x, M, Tx, statevariables), coeffs) + y = prop_and_combine(x, M, Tx, statevariables, coeffs) + y₀ = prop_and_combine(x₀, M, Tx, statevariables, coeffs) + + E = impulseresponse(prop_and_combine, x₀, M, Tx, statevariables, coeffs) + + # Does E matrix work properly? + ỹ = E*UnitfulMatrix(vec(x)) + for jj in eachindex(vec(y)) + @test isapprox.(vec(y)[jj],getindexqty(ỹ,jj)) + end + + x̂ = E\UnitfulMatrix(vec(y)) + for jj in eachindex(vec(y)) + @test isapprox.(vec(y)[jj],getindexqty(E*x̂,jj)) + end + + σₙ = 0.01 + σₓ = 100.0 + Cnn = UnitfulMatrix(Diagonal(fill(σₙ,length(y))),vec(unit.(y)).^1,vec(unit.(y)).^-1,exact=true) + + Cxx = UnitfulMatrix(Diagonal(fill(σₓ,length(x₀))),vec(unit.(x₀)),vec(unit.(x₀)).^-1,exact=true) + + problem = UnderdeterminedProblem(UnitfulMatrix(vec(y)),E,Cnn,Cxx,x₀) + x̃ = solve(problem) + for jj in eachindex(vec(y)) + @test within(y[jj],getindexqty(E*x̃.v,jj),3σₙ) # within 3-sigma + end + + @test cost(x̃,problem) < 0.5 # ad-hoc choice + + end + @testset "source water inversion: obs TIMESERIES, many surface regions, with circulation lag" begin @dim YearCE "years Common Era" From 25c1934ae1cef76894d19b5b351d0deb42f5785d Mon Sep 17 00:00:00 2001 From: b-r-hamilton Date: Wed, 29 Mar 2023 21:04:51 -0400 Subject: [PATCH 4/8] slightly cleaned up second runtest --- test/runtests.jl | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 1043014..2d04a43 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -371,18 +371,16 @@ include("test_functions.jl") Tx = first(dims(x)) x₀ = copy(x).*0 - y_ = convolve(x, M, Tx, statevariables) + #choose random linear coefficients coeffs = UnitfulMatrix(reshape(rand(2), (2,1)), [unit(1.0), K], [K*permil^-1]) #function to linearly combine state variable to one variable combine(y_, coeffs) = [getindexqty(UnitfulMatrix(transpose(y__))*coeffs,1,1) for y__ in y_] - #test! - combine(y_, coeffs) #one function to hand to impulseresponse prop_and_combine(x, M, Tx, statevariables, coeffs) = combine(convolve(x, M, Tx, statevariables), coeffs) - y = prop_and_combine(x, M, Tx, statevariables, coeffs) - y₀ = prop_and_combine(x₀, M, Tx, statevariables, coeffs) + #synthetic timeseries + y = prop_and_combine(x, M, Tx, statevariables, coeffs) E = impulseresponse(prop_and_combine, x₀, M, Tx, statevariables, coeffs) # Does E matrix work properly? @@ -390,18 +388,18 @@ include("test_functions.jl") for jj in eachindex(vec(y)) @test isapprox.(vec(y)[jj],getindexqty(ỹ,jj)) end - x̂ = E\UnitfulMatrix(vec(y)) for jj in eachindex(vec(y)) @test isapprox.(vec(y)[jj],getindexqty(E*x̂,jj)) end + #generate Cnn and Cxx matrices σₙ = 0.01 σₓ = 100.0 Cnn = UnitfulMatrix(Diagonal(fill(σₙ,length(y))),vec(unit.(y)).^1,vec(unit.(y)).^-1,exact=true) - Cxx = UnitfulMatrix(Diagonal(fill(σₓ,length(x₀))),vec(unit.(x₀)),vec(unit.(x₀)).^-1,exact=true) + #create problem and solve problem = UnderdeterminedProblem(UnitfulMatrix(vec(y)),E,Cnn,Cxx,x₀) x̃ = solve(problem) for jj in eachindex(vec(y)) From 3bf58a48973064088f9e429a686752dc1f86d547 Mon Sep 17 00:00:00 2001 From: b-r-hamilton Date: Thu, 30 Mar 2023 21:49:10 -0400 Subject: [PATCH 5/8] exchanged prop + combine method for adapted `convolve` --- src/BLUEs.jl | 18 +++++++++--------- test/runtests.jl | 35 +++++++++++++---------------------- 2 files changed, 22 insertions(+), 31 deletions(-) diff --git a/src/BLUEs.jl b/src/BLUEs.jl index 3223fda..cd39a3a 100644 --- a/src/BLUEs.jl +++ b/src/BLUEs.jl @@ -564,17 +564,17 @@ function convolve(x::AbstractDimArray,E::AbstractDimArray) return sum([E[ii,:] ⋅ x[Near(tnow-ll),:] for (ii,ll) in enumerate(lags)]) end """ -function convolve(x::AbstractDimArray, E::AbstractDimArray, statevars::Vector{Symbol} - This is a silly method - the info in statevars is a dimension in x - I can't figure out how to have one convolve method for a 2D AbstractDimArray - and another convolve method for a 3D AbstractDimArray - AbstractDimArray{Any, 3} does not work, even though I really want it to +function convolve(x::AbstractDimArray, E::AbstractDimArray, coeffs::UnitfulMatrix} + the `coeffs` argument signifies that x is a 3D array (i.e. >1 state variables) + + this function both convolves, and linearly combines the propagated state variables """ -function convolve(x::AbstractDimArray,E::AbstractDimArray, statevars::Vector{Symbol}) - mat = [convolve(x[:,:,At(s)], E) for s in statevars] - #return DimArray(mat, statevars) - return mat +function convolve(x::AbstractDimArray,E::AbstractDimArray, coeffs::UnitfulMatrix) + statevars = x.dims[3] + mat = UnitfulMatrix(transpose([convolve(x[:,:,At(s)], E) for s in statevars])) * coeffs + return getindexqty(mat, 1,1) end + function convolve(x::AbstractDimArray,E::AbstractDimArray,t::Number) lags = first(dims(E)) return sum([E[ii,:] ⋅ x[Near(t-ll),:] for (ii,ll) in enumerate(lags)]) diff --git a/test/runtests.jl b/test/runtests.jl index 2d04a43..b7d1b58 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -302,7 +302,7 @@ include("test_functions.jl") @dim InteriorLocation "interior location" @dim StateVariable "state variable" - #yr = u"yr" + yr = u"yr" nτ = 5 # how much of a lag is possible? lags = (0:(nτ-1))yr surfaceregions = [:NATL,:ANT,:SUBANT] @@ -314,38 +314,29 @@ include("test_functions.jl") statevariables = [:θ, :d18O] x = source_water_solution(surfaceregions,years, statevariables) - #convolve(x, M, s) gives us x variables propagated to location - #here is a function to linearly combine those two propagated value into a single value coeffs = UnitfulMatrix(reshape(rand(2), (2,1)), [unit(1.0), K], [K*permil^-1]) - combine(y, coeffs) = UnitfulMatrix(transpose(y)) * coeffs - y = combine(convolve(x, M, statevariables), coeffs) - #let's make one function to give to impulseresponse - prop_and_combine(x, M, statevariables, coeffs) = getindexqty(combine(convolve(x, M, statevariables), coeffs), 1, 1) - @test getindexqty(y,1,1) .== prop_and_combine(x, M, statevariables, coeffs) - + y = convolve(x, M, coeffs) + x₀ = copy(x).*0 - y₀ = convolve(x₀, M, statevariables) #seems to work - E = impulseresponse(prop_and_combine,x₀,M, statevariables, coeffs) - #does it work? - @test getindexqty(E*UnitfulMatrix(vec(x₀)),1) == prop_and_combine(x₀, M, statevariables, coeffs) + E = impulseresponse(convolve,x₀,M,coeffs) #E matrix tests from prior example ỹ = E*UnitfulMatrix(vec(x)) - @test isapprox(getindexqty(y,1,1),getindexqty(ỹ,1)) + @test isapprox(y,getindexqty(ỹ,1)) #had to get a little crafty to make this work... - x̂ = E\UnitfulMatrix(parent(y), [permil], [unit(1.0)]) - @test isapprox(getindexqty(y,1,1),getindexqty(E*x̂,1,1)) + x̂ = E\UnitfulMatrix(reshape([ustrip.(y)], (1,1)), [permil], [unit(1.0)]) + @test isapprox(y,getindexqty(E*x̂,1,1)) σₙ = 0.01 σₓ = 100.0 - Cnn = UnitfulMatrix(Diagonal(fill(σₙ,length(y))),fill(unit(getindexqty(y, 1,1)).^1,length(y)),fill(unit(getindexqty(y,1,1)).^-1,length(y)),exact=false) + Cnn = UnitfulMatrix(Diagonal(fill(σₙ.^2,length(y))),fill(unit(y).^1,length(y)),fill(unit(y).^-1,length(y)),exact=false) - Cxx = UnitfulMatrix(Diagonal(fill(σₓ,length(x₀))),unit.(x₀)[:],unit.(x₀)[:].^-1,exact=false) - problem = UnderdeterminedProblem(UnitfulMatrix([getindexqty(y,1,1)]), E, Cnn, Cxx, x₀) + Cxx = UnitfulMatrix(Diagonal(fill(σₓ.^2,length(x₀))),unit.(x₀)[:],unit.(x₀)[:].^-1,exact=false) + problem = UnderdeterminedProblem(UnitfulMatrix([y]), E, Cnn, Cxx, x₀) x̃ = solve(problem) - @test within(getindexqty(y, 1,1), getindexqty(E*x̃.v, 1), 3σₙ) + @test within(y, getindexqty(E*x̃.v, 1), 3σₙ) @test cost(x̃,problem) < 5e-2 end @@ -396,8 +387,8 @@ include("test_functions.jl") #generate Cnn and Cxx matrices σₙ = 0.01 σₓ = 100.0 - Cnn = UnitfulMatrix(Diagonal(fill(σₙ,length(y))),vec(unit.(y)).^1,vec(unit.(y)).^-1,exact=true) - Cxx = UnitfulMatrix(Diagonal(fill(σₓ,length(x₀))),vec(unit.(x₀)),vec(unit.(x₀)).^-1,exact=true) + Cnn = UnitfulMatrix(Diagonal(fill(σₙ.^2,length(y))),vec(unit.(y)).^1,vec(unit.(y)).^-1,exact=true) + Cxx = UnitfulMatrix(Diagonal(fill(σₓ.^2,length(x₀))),vec(unit.(x₀)),vec(unit.(x₀)).^-1,exact=true) #create problem and solve problem = UnderdeterminedProblem(UnitfulMatrix(vec(y)),E,Cnn,Cxx,x₀) From 217d0fc42f3c9985a5b400e37383a8dbc0e6fea1 Mon Sep 17 00:00:00 2001 From: b-r-hamilton Date: Thu, 30 Mar 2023 22:01:53 -0400 Subject: [PATCH 6/8] `convolve` updated for timeseries case --- src/BLUEs.jl | 15 +++++++++------ test/runtests.jl | 18 ++++++------------ 2 files changed, 15 insertions(+), 18 deletions(-) diff --git a/src/BLUEs.jl b/src/BLUEs.jl index cd39a3a..3fc7d20 100644 --- a/src/BLUEs.jl +++ b/src/BLUEs.jl @@ -579,15 +579,18 @@ function convolve(x::AbstractDimArray,E::AbstractDimArray,t::Number) lags = first(dims(E)) return sum([E[ii,:] ⋅ x[Near(t-ll),:] for (ii,ll) in enumerate(lags)]) end -# two more silly convolve function -function convolve(x::AbstractDimArray, E::AbstractDimArray, t::Number, statevars::Vector{Symbol}) - mat = [convolve(x[:,:,At(s)], E, t) for s in statevars] - return mat + +#coeffs signifies that x is 3D +function convolve(x::AbstractDimArray, E::AbstractDimArray, t::Number, coeffs::UnitfulMatrix) + statevars = x.dims[3] + mat = UnitfulMatrix(transpose([convolve(x[:,:,At(s)], E, t) for s in statevars]))*coeffs + return getindexqty(mat, 1,1) end + #don't handle the ndims(M) == 3 case here but I'll get back to it -function convolve(x::AbstractDimArray, M::AbstractDimArray, Tx::Union{Ti, Vector}, statevars::Vector{Symbol}) +function convolve(x::AbstractDimArray, M::AbstractDimArray, Tx::Union{Ti, Vector}, coeffs::UnitfulMatrix) if ndims(M) == 2 - return DimArray([convolve(x,M,Tx[tt], statevars) for (tt, yy) in enumerate(Tx)], Tx) + return DimArray([convolve(x,M,Tx[tt], coeffs) for (tt, yy) in enumerate(Tx)], Tx) elseif ndims(M) == 3 error("some code should go here") else diff --git a/test/runtests.jl b/test/runtests.jl index b7d1b58..fbe4695 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -281,12 +281,11 @@ include("test_functions.jl") # unitful scalars in UnitfulLinearAlgebra σₙ = 0.01 - σₓ_d18O = 100.0 - σₓ_θ = 100.0 + σₓ = 100.0 - Cnn = UnitfulMatrix(Diagonal(fill(σₙ,length(y))),fill(unit.(y).^1,length(y)),fill(unit.(y).^-1,length(y)),exact=false) + Cnn = UnitfulMatrix(Diagonal(fill(σₙ^2,length(y))),fill(unit.(y).^1,length(y)),fill(unit.(y).^-1,length(y)),exact=false) - Cxx = UnitfulMatrix(Diagonal(fill(σₓ,length(x₀))),unit.(x₀)[:],unit.(x₀)[:].^-1,exact=false) + Cxx = UnitfulMatrix(Diagonal(fill(σₓ^2,length(x₀))),unit.(x₀)[:],unit.(x₀)[:].^-1,exact=false) problem = UnderdeterminedProblem(UnitfulMatrix([y]),E,Cnn,Cxx,x₀) x̃ = solve(problem) @@ -364,15 +363,10 @@ include("test_functions.jl") #choose random linear coefficients coeffs = UnitfulMatrix(reshape(rand(2), (2,1)), [unit(1.0), K], [K*permil^-1]) - #function to linearly combine state variable to one variable - combine(y_, coeffs) = [getindexqty(UnitfulMatrix(transpose(y__))*coeffs,1,1) for y__ in y_] - - #one function to hand to impulseresponse - prop_and_combine(x, M, Tx, statevariables, coeffs) = combine(convolve(x, M, Tx, statevariables), coeffs) - + #synthetic timeseries - y = prop_and_combine(x, M, Tx, statevariables, coeffs) - E = impulseresponse(prop_and_combine, x₀, M, Tx, statevariables, coeffs) + y = convolve(x, M, Tx, coeffs) + E = impulseresponse(convolve, x₀, M, Tx,coeffs) # Does E matrix work properly? ỹ = E*UnitfulMatrix(vec(x)) From d55ff8dcc1e169b258735bb2bbb5c4afe7177f2a Mon Sep 17 00:00:00 2001 From: b-r-hamilton Date: Thu, 30 Mar 2023 22:08:13 -0400 Subject: [PATCH 7/8] `show` method for 3D x (needs some help) --- src/BLUEs.jl | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/src/BLUEs.jl b/src/BLUEs.jl index 3fc7d20..2555995 100644 --- a/src/BLUEs.jl +++ b/src/BLUEs.jl @@ -102,6 +102,19 @@ function show(io::IO, mime::MIME{Symbol("text/plain")}, x::Union{DimEstimate,Est # show(io, mime, F.V) end +function show(io::IO, mime::MIME{Symbol("text/plain")}, x::DimArray{Quantity{Float64}, 3}) + summary(io, x); println(io) + statevars = x.dims[3] + for (i, s) in enumerate(statevars) + if i != 1 + println() + end + + println(io, "State Variable " * string(i) * ": " * string(s)) + show(io, mime, x[:,:,At(s)]) + end +end + """ function getproperty(x::Estimate, d::Symbol) From 1ceb3b01df70c3fd5a36a380cc5899f61a4c2478 Mon Sep 17 00:00:00 2001 From: b-r-hamilton Date: Thu, 13 Apr 2023 16:10:08 -0400 Subject: [PATCH 8/8] datacost + controlcost checks --- test/runtests.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index fbe4695..c0f056a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -336,7 +336,8 @@ include("test_functions.jl") problem = UnderdeterminedProblem(UnitfulMatrix([y]), E, Cnn, Cxx, x₀) x̃ = solve(problem) @test within(y, getindexqty(E*x̃.v, 1), 3σₙ) - @test cost(x̃,problem) < 5e-2 + @test cost(x̃,problem) < 5e-2 + @test cost(x̃, problem) == datacost(x̃, problem) + controlcost(x̃, problem) end @testset "source water inversion: obs TIMESERIES, many surface regions, with circulation lag, TWO STATE VARIABLES" begin @@ -448,6 +449,9 @@ include("test_functions.jl") # no noise in obs but some control penalty @test cost(x̃,problem) < 0.5 # ad-hoc choice + @test cost(x̃, problem) == datacost(x̃, problem) + controlcost(x̃, problem) + + end @testset "source water inversion: MANY OBS TIMESERIES, many surface regions, with circulation lag" begin