Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Multipliable dimarrays #58

Closed
wants to merge 60 commits into from
Closed
Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
60 commits
Select commit Hold shift + click to select a range
ec5e0ee
separate out DimEstimate code
ggebbie May 28, 2024
d62e679
start making tests w or w/o units
ggebbie May 28, 2024
24610e5
split tests into `Estimate` component
ggebbie May 29, 2024
0d63fe6
rm old test file
ggebbie May 29, 2024
8c58bf9
uncomment
ggebbie May 29, 2024
592c4a3
extend `(\)`
ggebbie May 29, 2024
5ed6164
specialize sigma to work with `AbstractDimArray`
ggebbie May 29, 2024
bfbf383
wip: deprecate Unitful Vectors
ggebbie May 30, 2024
54346dc
compat for DimensionalData 0.27
ggebbie May 30, 2024
9f0dc3a
poly test with + w/o units, no DimArray
ggebbie May 30, 2024
e55b09a
`test_estimate` passes with/w/o units
ggebbie May 30, 2024
e038ffb
first wm exmple: Estimate w/ DimArray
ggebbie May 30, 2024
b169f67
first wm example with units
ggebbie May 30, 2024
81d1b0b
update `Missing` struct field type
ggebbie May 31, 2024
7cd087a
minor change for ULA convenience
ggebbie Jun 13, 2024
944916a
new source file for Dimensional stuff
ggebbie Jun 17, 2024
1693ee6
add new blockdim files
ggebbie Jun 17, 2024
d2c2036
simple `convolve` returns `DimArray`
ggebbie Jun 17, 2024
47d6637
start dispatching `standard_error`
ggebbie Jun 17, 2024
867ff0d
split unitful material
ggebbie Jun 18, 2024
7e6c656
split src to under/over-determined files
ggebbie Jun 18, 2024
ef96c1d
`NamedTuple` conveniences in sep. file
ggebbie Jun 18, 2024
0f7250e
ignore auto saved files in git
ggebbie Jun 18, 2024
8f9a2da
simple `combine` `Estimate`s
ggebbie Jun 18, 2024
580a738
update overdetermined/mixed signals test
ggebbie Jun 18, 2024
e459bf4
wip: mixed signals fails w/units
ggebbie Jun 18, 2024
521d056
mixed signals test passes w w/o units
ggebbie Jun 19, 2024
44590df
use diag-regression branch: UnitfulLinearAlgebra
ggebbie Jun 20, 2024
5f48639
remove `update`
ggebbie Jun 20, 2024
e1bde2b
standard error for `DimArray`s
ggebbie Jun 20, 2024
6e4ae0b
rename blockdim -> dimensional_data
ggebbie Jun 20, 2024
8c059e9
update `show` for `DimArray`s
ggebbie Jun 20, 2024
d45eda4
bugfix: `observematrix`
ggebbie Jun 20, 2024
3234a56
wip: 1st water-mass test
ggebbie Jun 21, 2024
74b16ef
first water-mass test returns obs
ggebbie Jun 21, 2024
1c5d65f
`algebraic_object` placeholder name
ggebbie Jun 22, 2024
ac9439e
wip: uncertainty calc. algebraic DimArray
ggebbie Jun 22, 2024
2ebc2e3
shorten func names that don't pirate
ggebbie Jun 22, 2024
80c2abf
`observematrix` -> `convolve` via dispatch
ggebbie Jun 22, 2024
9e8eb54
passes: case false,false,false, no units
ggebbie Jun 22, 2024
8fa4bbd
`standard_error` with units
ggebbie Jun 22, 2024
8ae8b28
bugfix: overdetermined problems
ggebbie Jun 22, 2024
9ce60fc
passes: f,f,f with units, first dimdata test
ggebbie Jun 22, 2024
8d45a7f
`show` DimArray more general for state vars
ggebbie Jun 22, 2024
23f80ad
`DimArray` for `coeffs`
ggebbie Jun 24, 2024
4c9fb1b
bugfix: years |> yrs
ggebbie Jun 24, 2024
f254f6d
statevars tests pass
ggebbie Jun 24, 2024
d5ae083
t,t,t w/o units passes
ggebbie Jun 24, 2024
6ddeae4
wip: units with DimArray + statevecs
ggebbie Jun 24, 2024
cdffc3c
wip: wrap DimArray with UnitfulMatrix
ggebbie Jun 26, 2024
4ee32f2
wip: unitfulmatrix with dimarray in ULA
ggebbie Jun 26, 2024
8c7eb7b
move code to `MultipliableDimArrays.jl`
ggebbie Jun 27, 2024
c3b3a66
add `MultipliableDimArrays`
ggebbie Jun 27, 2024
2aa30cc
update first water-mass test for MDA
ggebbie Jun 27, 2024
1bdc8c9
subset of tests w/MultipiableDimArrays
ggebbie Jun 27, 2024
d8c99dc
remove code that is duplicated in MDA.jl
ggebbie Jun 27, 2024
5a63cb9
bump project
ggebbie Jun 27, 2024
2e666cd
add Manifest for unregistered MDA
ggebbie Jun 27, 2024
64b7d2f
next step: marked
ggebbie Jun 28, 2024
e99360d
use GitHub available packages
ggebbie Jun 29, 2024
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
Prev Previous commit
Next Next commit
update Missing struct field type
  • Loading branch information
ggebbie committed May 31, 2024
commit 81d1b0bbaff9a8e0d173b2e1c924d46ce9780da0
4 changes: 2 additions & 2 deletions src/BLUEs.jl
Original file line number Diff line number Diff line change
Expand Up @@ -71,8 +71,8 @@ struct UnderdeterminedProblem
y :: AbstractVector
E :: AbstractMatrix
Cnn :: AbstractMatrix
Cxx :: Union{AbstractMatrix}
x₀ :: Union{AbstractVector,AbstractDimArray,Missing}
Cxx :: Union{AbstractMatrix, Missing}
x₀ :: Union{AbstractVector, AbstractDimArray, Missing}
end

#include("dim_estimate.jl")
Expand Down
53 changes: 35 additions & 18 deletions test/test_dimestimate.jl
Original file line number Diff line number Diff line change
Expand Up @@ -7,10 +7,10 @@
@dim StateVariable "state variable"
yr = u"yr"
surfaceregions = [:NATL,:ANT,:SUBANT]
n = length(surfaceregions)
N = length(surfaceregions)
years = (1990:2000)yr
statevariables = [:θ, :δ¹⁸O]
m = 5 # Interior Locations with obs
M = 5 # Interior Locations with obs

@testset "E matrix = UnitfulDimMatrix (many surface regions, one state variable, multiple obs locations but no obs timeseries, no circulation lag)" begin

Expand All @@ -19,15 +19,15 @@
# 2) Obs at one end time
# 3) No circulation lag
σₙ = 1.0
E,x = random_source_water_matrix_vector_pair(m)
E,x = random_source_water_matrix_vector_pair(M)
if use_units
Cnndims = (first(dims(E)),first(dims(E)))
Cnn⁻¹ = Diagonal(fill(σₙ^-1,m),unitrange(E).^-1,unitrange(E),Cnndims,exact=true)
Cnn⁻¹ = Diagonal(fill(σₙ^-1,M),unitrange(E).^-1,unitrange(E),Cnndims,exact=true)
else
E = DimArray(parent(E),dims(E))
x = DimArray(parent(x),dims(x))
Cnndims = (first(dims(E)),first(dims(E)))
Cnn⁻¹ = DimArray(Diagonal(fill(σₙ^-1,m)),Cnndims)
Cnn⁻¹ = DimArray(Diagonal(fill(σₙ^-1,M)),Cnndims)
end

# Run model to predict interior location temperature
Expand Down Expand Up @@ -63,29 +63,40 @@

#(statevars,timeseries,lag) = cases[1] # for interactive use
for (statevars,timeseries,lag) in cases
println("statevars,timeseries,lag = ",statevars,timeseries,lag)
println("statevars,timeseries,lag = ",statevars, " ", timeseries, " ", lag)

#define constants
lag ? nτ = 5 : nτ = 1
lags = (0:(nτ-1))yr
σₙ = 0.01
σₓ = 100.0

# take all years or just the end year
timeseries ? yrs = years : yrs = [years[end]]

M = source_water_matrix_with_lag(surfaceregions,lags)

# Step 1: get synthetic solution
if !statevars

x = source_water_solution(surfaceregions,yrs)

# DimArray is good enough. This is an array, not necessarily a matrix.
x₀ = DimArray(zeros(size(x))K,(Ti(yrs),last(dims(M))))

if use_units
# DimArray is good enough. This is an array, not necessarily a matrix.
x₀ = DimArray(zeros(size(x))K,(Ti(yrs),last(dims(M))))
else
x = DimArray(ustrip.(parent(x)), dims(x))
x₀ = DimArray(zeros(size(x)),(Ti(yrs),last(dims(M))))
end

elseif statevars

x = source_water_solution(surfaceregions,years, statevariables)
coeffs = UnitfulMatrix(rand(2,1), [NoUnits, K], [K/permil])
x₀ = copy(x).*0
if use_units
else
end

else
error("no case")
Expand All @@ -105,7 +116,7 @@
global predict(x) = convolve(x,M)
end

# probe to get E matrix. Use function convolve
# probe to get E matrix. Use a convolution.
E = impulseresponse(predict,x₀)

# Given, M and x. Make synthetic data for observations.
Expand All @@ -115,12 +126,14 @@

# test compatibility
@test first(E*UnitfulMatrix(vec(x₀))) .== first(predict(x₀))
@test first(E*vec(x₀)) .== first(predict(x₀)) # not necessary to transfer x₀ to UnitfulMatrix

# Julia doesn't have method to make scalar into vector
!(y isa AbstractVector) && (y = [y])

# Does E matrix work properly?
ỹ = E*UnitfulMatrix(vec(x))
ỹ = E*UnitfulMatrix(vec(x)) # not necessary
ỹ = E*vec(x)
for jj in eachindex(y)
@test isapprox(vec(y)[jj],vec(ỹ)[jj])
end
Expand All @@ -130,12 +143,16 @@
end

# now in a position to use BLUEs to solve
σₙ = 0.01
σₓ = 100.0

Cnn = Diagonal(fill(σₙ^2,length(y)),unit.(y),unit.(y).^-1)
Cxx = Diagonal(fill(σₓ^2,length(x₀)),vec(unit.(x₀)),vec(unit.(x₀)).^-1)
problem = UnderdeterminedProblem(UnitfulMatrix(y),E,Cnn,Cxx,x₀)
if use_units
Cnn = Diagonal(fill(σₙ^2,length(y)),unit.(y),unit.(y).^-1)
Cxx = Diagonal(fill(σₓ^2,length(x₀)),vec(unit.(x₀)),vec(unit.(x₀)).^-1)
else
Cnn = Diagonal(fill(σₙ^2,length(y)))
Cxx = Diagonal(fill(σₓ^2,length(x₀)))
end

#problem = UnderdeterminedProblem(UnitfulMatrix(y),E,Cnn,Cxx,x₀)
problem = UnderdeterminedProblem(y,E,Cnn,Cxx,x₀)

# when x₀ is a DimArray, then x̃ is a DimEstimate
x̃ = solve(problem) # ::DimEstimate
Expand Down