From cfd7369f3acddd650515f923501a1695c1f2fb71 Mon Sep 17 00:00:00 2001 From: Thomas Wutzler Date: Mon, 28 Apr 2025 10:31:50 +0200 Subject: [PATCH 1/7] test matrix version and neglogden on GPU --- src/ComponentArrayInterpreter.jl | 47 ++++++++++++-- src/DoubleMM/f_doubleMM.jl | 53 ++++++++++++---- src/HybridVariationalInference.jl | 2 +- src/logden_normal.jl | 3 +- test/test_ComponentArrayInterpreter.jl | 6 ++ test/test_doubleMM.jl | 85 ++++++++++++++++++++++++-- 6 files changed, 171 insertions(+), 25 deletions(-) diff --git a/src/ComponentArrayInterpreter.jl b/src/ComponentArrayInterpreter.jl index 46cd5f2..6f89b35 100644 --- a/src/ComponentArrayInterpreter.jl +++ b/src/ComponentArrayInterpreter.jl @@ -5,9 +5,12 @@ Interface for Type that implements - `as_ca(::AbstractArray, interpreter) -> ComponentArray` +- `CompoentArrays.getaxes(interpeter)` - `Base.length(interpreter) -> Int` When called on a vector, forwards to `as_ca`. + +There is a default implementation for Base.length based on CompoentArrays.getaxes. """ abstract type AbstractComponentArrayInterpreter end @@ -18,6 +21,11 @@ Returns a ComponentArray with underlying data `v`. """ function as_ca end +function Base.length(cai::AbstractComponentArrayInterpreter) + prod(_axis_length.(CA.getaxes(cai))) +end + + (interpreter::AbstractComponentArrayInterpreter)(v::AbstractArray) = as_ca(v, interpreter) """ @@ -36,9 +44,13 @@ function as_ca(v::AbstractArray, ::StaticComponentArrayInterpreter{AX}) where {A CA.ComponentArray(vr, AX) end -function Base.length(::StaticComponentArrayInterpreter{AX}) where {AX} - #sum(length, typeof(AX).parameters[1]) - prod(_axis_length.(AX)) +# function Base.length(::StaticComponentArrayInterpreter{AX}) where {AX} +# #sum(length, typeof(AX).parameters[1]) +# prod(_axis_length.(AX)) +# end + +function CA.getaxes(int::StaticComponentArrayInterpreter{AX}) where {AX} + AX end get_concrete(cai::StaticComponentArrayInterpreter) = cai @@ -63,10 +75,11 @@ function as_ca(v::AbstractArray, cai::ComponentArrayInterpreter) CA.ComponentArray(vr, cai.axes) end -function Base.length(cai::ComponentArrayInterpreter) - prod(_axis_length.(cai.axes)) +function CA.getaxes(cai::ComponentArrayInterpreter) + cai.axes end + get_concrete(cai::ComponentArrayInterpreter) = StaticComponentArrayInterpreter{cai.axes}() @@ -120,6 +133,10 @@ function ComponentArrayInterpreter( ca::CA.AbstractComponentArray, n_dims::NTuple{N,<:Integer}) where N ComponentArrayInterpreter(CA.getaxes(ca), n_dims) end +function ComponentArrayInterpreter( + cai::AbstractComponentArrayInterpreter, n_dims::NTuple{N,<:Integer}) where N + ComponentArrayInterpreter(CA.getaxes(cai), n_dims) +end function ComponentArrayInterpreter( axes::NTuple{M, <:CA.AbstractAxis}, n_dims::NTuple{N,<:Integer}) where {M,N} axes_ext = (axes..., map(n_dim -> CA.Axis(i=1:n_dim), n_dims)...) @@ -131,12 +148,17 @@ function ComponentArrayInterpreter( n_dims::NTuple{N,<:Integer}, ca::CA.AbstractComponentArray) where N ComponentArrayInterpreter(n_dims, CA.getaxes(ca)) end +function ComponentArrayInterpreter( + n_dims::NTuple{N,<:Integer}, cai::AbstractComponentArrayInterpreter) where N + ComponentArrayInterpreter(n_dims, CA.getaxes(cai)) +end function ComponentArrayInterpreter( n_dims::NTuple{N,<:Integer}, axes::NTuple{M, <:CA.AbstractAxis}) where {N,M} axes_ext = (map(n_dim -> CA.Axis(i=1:n_dim), n_dims)..., axes...) ComponentArrayInterpreter(axes_ext) end + # ambuiguity with two empty Tuples (edge prob that does not make sense) # Empty ComponentVector with no other array dimensions -> empty componentVector function ComponentArrayInterpreter(n_dims1::Tuple{}, n_dims2::Tuple{}) @@ -156,6 +178,8 @@ _axis_length(::CA.FlatAxis) = 0 _axis_length(::CA.UnitRange) = 0 """ + flatten1(cv::CA.ComponentVector) + Removes the highest level of keys. Keeps the reference to the underlying data, but changes the axis. If first-level vector has no sub-names, an error (Aguement Error tuple must be non-empty) @@ -174,3 +198,16 @@ function flatten1(cv::CA.ComponentVector) CA.ComponentVector(cv, first(CA.getaxes(cv_new))) end end + + +""" + get_positions(cai::AbstractComponentArrayInterpreter) + +Create a NamedTuple of integer indices for each component. +Assumes that interpreter results in a one-dimensional array, i.e. in a ComponentVector. +""" +function get_positions(cai::AbstractComponentArrayInterpreter) + @assert length(CA.getaxes(cai)) == 1 + cv = cai(1:length(cai)) + (; (k => cv[k] for k in keys(cv))... ) +end diff --git a/src/DoubleMM/f_doubleMM.jl b/src/DoubleMM/f_doubleMM.jl index 18d1f4e..b68f1fe 100644 --- a/src/DoubleMM/f_doubleMM.jl +++ b/src/DoubleMM/f_doubleMM.jl @@ -19,34 +19,61 @@ function f_doubleMM(θ::AbstractVector, x, intθ) θc = intθ(θ) #using ComponentArrays: ComponentArrays as CA #r0, r1, K1, K2 = θc[(:r0, :r1, :K1, :K2)] # does not work on Zygote+GPU - r0 = θc[:r0] - r1 = θc[:r1] - K1 = θc[:K1] - K2 = θc[:K2] + (r0, r1, K1, K2) = map((:r0, :r1, :K1, :K2)) do par + # vector will be repeated when broadcasted by a matrix + CA.getdata(θc[par]) + end + # r0 = θc[:r0] + # r1 = θc[:r1] + # K1 = θc[:K1] + # K2 = θc[:K2] y = r0 .+ r1 .* x.S1 ./ (K1 .+ x.S1) .* x.S2 ./ (K2 .+ x.S2) end return (y) end -function f_doubleMM(θ::AbstractMatrix, x::NTuple{N, AbstractMatrix}, intθ) where N +function f_doubleMM(θ::AbstractMatrix, x::NamedTuple, intθ::HVI.AbstractComponentArrayInterpreter) # provide θ for n_row sites # provide x.S1 as Matrix n_site x n_obs # extract parameters not depending on order, i.e whether they are in θP or θM θc = intθ(θ) - #using ComponentArrays: ComponentArrays as CA - #r0, r1, K1, K2 = θc[(:r0, :r1, :K1, :K2)] # does not work on Zygote+GPU @assert size(x.S1,1) == size(θ,1) # same number of sites @assert size(x.S1) == size(x.S2) # same number of observations #@assert length(x.s2 == n_obs) - r0 = θc[:,:r0] # vector will be repeated when broadcasted by a matrix - r1 = θc[:,:r1] - K1 = θc[:,:K1] - K2 = θc[:,:K2] + # problems on AD on GPU with indexing CA may be related to printing result, use ";" + (r0, r1, K1, K2) = map((:r0, :r1, :K1, :K2)) do par + # vector will be repeated when broadcasted by a matrix + CA.getdata(θc[:,par]) + end + # r0 = CA.getdata(θc[:,:r0]) # vector will be repeated when broadcasted by a matrix + # r1 = CA.getdata(θc[:,:r1]) + # K1 = CA.getdata(θc[:,:K1]) + # K2 = CA.getdata(θc[:,:K2]) y = r0 .+ r1 .* x.S1 ./ (K1 .+ x.S1) .* x.S2 ./ (K2 .+ x.S2) - return (y) + return (y) end - +# function f_doubleMM(θ::AbstractMatrix, x::NamedTuple, θpos::NamedTuple) +# # provide θ for n_row sites +# # provide x.S1 as Matrix n_site x n_obs +# # extract parameters not depending on order, i.e whether they are in θP or θM +# @assert size(x.S1,1) == size(θ,1) # same number of sites +# @assert size(x.S1) == size(x.S2) # same number of observations +# (r0, r1, K1, K2) = map((:r0, :r1, :K1, :K2)) do par +# # vector will be repeated when broadcasted by a matrix +# CA.getdata(θ[:,θpos[par]]) +# end +# # r0 = CA.getdata(θ[:,θpos.r0]) # vector will be repeated when broadcasted by a matrix +# # r1 = CA.getdata(θ[:,θpos.r1]) +# # K1 = CA.getdata(θ[:,θpos.K1]) +# # K2 = CA.getdata(θ[:,θpos.K2]) +# #y = r0 .+ r1 +# #y = x.S1 + x.S2 +# #y = (K1 .+ x.S1) +# #y = r1 .* x.S1 ./ (K1 .+ x.S1) +# y = r0 .+ r1 .* x.S1 ./ (K1 .+ x.S1) .* x.S2 ./ (K2 .+ x.S2) +# return (y) +# end function HVI.get_hybridproblem_par_templates(::DoubleMMCase; scenario::NTuple = ()) if (:omit_r0 ∈ scenario) diff --git a/src/HybridVariationalInference.jl b/src/HybridVariationalInference.jl index efb2f3f..e8ad83a 100644 --- a/src/HybridVariationalInference.jl +++ b/src/HybridVariationalInference.jl @@ -20,7 +20,7 @@ using Distributions, DistributionFits using StaticArrays: StaticArrays as SA using Functors -export ComponentArrayInterpreter, flatten1, get_concrete +export ComponentArrayInterpreter, flatten1, get_concrete, get_positions include("ComponentArrayInterpreter.jl") export AbstractModelApplicator, construct_ChainsApplicator diff --git a/src/logden_normal.jl b/src/logden_normal.jl index 50a8509..e25fafb 100644 --- a/src/logden_normal.jl +++ b/src/logden_normal.jl @@ -13,7 +13,8 @@ a low uncertainty estimate and means closer to the observations to help an initial fit. The obtained parameters then can be used as starting values for a the proper fit with `σfac=1.0`. """ -function neg_logden_indep_normal(obs::AbstractArray, μ::AbstractArray, logσ2::AbstractArray; σfac=1.0) +function neg_logden_indep_normal(obs::AbstractArray, μ::AbstractArray, logσ2::AbstractArray{ET}; + σfac=one(ET)) where ET # log of independent Normal distributions # estimate independent uncertainty of each θM, rather than full covariance #nlogL = sum(σfac .* log.(σs) .+ 1 / 2 .* abs2.((obs .- μ) ./ σs)) diff --git a/test/test_ComponentArrayInterpreter.jl b/test/test_ComponentArrayInterpreter.jl index f392dae..f3b8fda 100644 --- a/test/test_ComponentArrayInterpreter.jl +++ b/test/test_ComponentArrayInterpreter.jl @@ -57,6 +57,8 @@ end; testm(mm) mmc = get_concrete(mm) testm(mmc) + mmi = ComponentArrayInterpreter(mv, (n_col,)) # construct on interpreter itself + testm(mmi) # n_z = 3 mm = ComponentArrayInterpreter(cv, (n_col, n_z)) @@ -67,6 +69,8 @@ end; end testm(mm) testm(get_concrete(mm)) + mmi = ComponentArrayInterpreter(mv, (n_col, n_z)) # construct on interpreter itself + testm(mmi) # n_row = 3 mm = ComponentArrayInterpreter((n_row,), cv) @@ -77,6 +81,8 @@ end; end testm(mm) testm(get_concrete(mm)) + mm = ComponentArrayInterpreter((n_row,), mv) # construct on interpreter itself + testm(mmi) end; @testset "empty ComponentVector" begin diff --git a/test/test_doubleMM.jl b/test/test_doubleMM.jl index 6b42803..2f59843 100644 --- a/test/test_doubleMM.jl +++ b/test/test_doubleMM.jl @@ -14,6 +14,14 @@ import Zygote using OptimizationOptimisers using MLDataDevices +using CUDA: CUDA +using Flux +using GPUArraysCore + +gdev = gpu_device() +cdev = cpu_device() + + const prob = DoubleMM.DoubleMMCase() scenario = (:default,) #using Flux @@ -32,6 +40,8 @@ rng = StableRNG(111) (; xM, n_site, θP_true, θMs_true, xP, y_global_true, y_true, y_global_o, y_o, y_unc ) = gen_hybridproblem_synthetic(rng, prob; scenario); i_sites = 1:n_site +fneglogden = get_hybridproblem_neg_logden_obs(prob; scenario) + @testset "gen_hybridproblem_synthetic" begin @test isapprox( @@ -45,6 +55,76 @@ i_sites = 1:n_site @test gen2.y_o == y_o end +@testset "f_doubleMM_Matrix" begin + is = repeat(axes(θP_true,1)', n_site) + θvec = CA.ComponentVector(P = θP_true, Ms = θMs_true) + xPM = map(xP1s -> repeat(xP1s', n_site), xP[1]) + #θ = hcat(θP_true[is], θMs_true') + intθ1 = get_concrete(ComponentArrayInterpreter(vcat(θP_true, θMs_true[:,1]))) + #θpos = get_positions(intθ1) + intθ = get_concrete(ComponentArrayInterpreter((n_site,), intθ1)) + fcost = (θvec, xPM, y_o, y_unc) -> begin + θ = hcat(CA.getdata(θvec.P[is]), CA.getdata(θvec.Ms')) + y = HVI.DoubleMM.f_doubleMM(θ, xPM, intθ) + #y = HVI.DoubleMM.f_doubleMM(θ, xPM, θpos) + fneglogden(y_o, y', y_unc) + end + cost = fcost(θvec, xPM, y_o, y_unc) + ygrad = Zygote.gradient(θv -> fcost(θv, xPM), θvec)[1]; + if gdev isa MLDataDevices.AbstractGPUDevice + # θg = gdev(θ) + # xPMg = gdev(xPM) + # yg = HVI.DoubleMM.f_doubleMM(θg, xPMg, intθ); + θvecg = gdev(θvec); + xPMg = gdev(xPM) + y_og = gdev(y_o) + y_uncg = gdev(y_unc) + costg = fcost(θvecg, xPMg, y_og, y_uncg); + ygradg = Zygote.gradient(θv -> fcost(θv, xPMg, y_og, y_uncg), θvecg)[1]; # erros without ";" + @test ygradg isa CA.ComponentArray + @test CA.getdata(ygradg) isa GPUArraysCore.AbstractGPUArray + ygradgc = HVI.apply_preserve_axes(cdev, ygradg) # can print the cpu version + # ygradgc.P .- ygrad.P + # ygradgc.Ms + end +end + +@testset "neg_logden_obs Matrix" begin + is = repeat(axes(θP_true,1)', n_site) + θvec = CA.ComponentVector(P = θP_true, Ms = θMs_true) + xPM = map(xP1s -> repeat(xP1s', n_site), xP[1]) + #θ = hcat(θP_true[is], θMs_true') + intθ1 = get_concrete(ComponentArrayInterpreter(vcat(θP_true, θMs_true[:,1]))) + #θpos = get_positions(intθ1) + intθ = get_concrete(ComponentArrayInterpreter((n_site,), intθ1)) + fl = (θvec, xPM) -> begin + θ = hcat(CA.getdata(θvec.P[is]), CA.getdata(θvec.Ms')) + y = HVI.DoubleMM.f_doubleMM(θ, xPM, intθ) + #y = HVI.DoubleMM.f_doubleMM(θ, xPM, θpos) + + end + y = fy(θvec, xPM) + y_exp = applyf(HVI.DoubleMM.f_doubleMM, θMs_true, θP_true, Vector{eltype(θP_true)}(undef, 0), xP, intθ1) + @test y == y_exp' + ygrad = Zygote.gradient(θv -> sum(fy(θv, xPM)), θvec)[1]; + if gdev isa MLDataDevices.AbstractGPUDevice + # θg = gdev(θ) + # xPMg = gdev(xPM) + # yg = HVI.DoubleMM.f_doubleMM(θg, xPMg, intθ); + θvecg = gdev(θvec); + xPMg = gdev(xPM) + yg = fy(θvecg, xPMg); + @test cdev(yg) == y_exp' + ygradg = Zygote.gradient(θv -> sum(fy(θv, xPMg)), θvecg)[1]; # erros without ";" + @test ygradg isa CA.ComponentArray + @test CA.getdata(ygradg) isa GPUArraysCore.AbstractGPUArray + ygradgc = HVI.apply_preserve_axes(cdev, ygradg) # can print the cpu version + # ygradgc.P .- ygrad.P + # ygradgc.Ms + end +end + + @testset "loss_g" begin g, ϕg0 = get_hybridproblem_MLapplicator(rng, prob; scenario) (; transP, transM) = get_hybridproblem_transforms(prob; scenario) @@ -141,12 +221,7 @@ end end end -using CUDA: CUDA -using Flux -using GPUArraysCore -gdev = gpu_device() -cdev = cpu_device() if gdev isa MLDataDevices.AbstractGPUDevice @testset "transfer NormalScalingModelApplicator to gpu" begin scenario = (:use_Flux,) From 05dfc1f781c79ff684009221617fefec6bbe6651 Mon Sep 17 00:00:00 2001 From: Thomas Wutzler Date: Tue, 29 Apr 2025 11:31:33 +0200 Subject: [PATCH 2/7] Implement Bijector Exp that also works on AD on GPU --- src/bijectors_exp.jl | 17 +++++++++ test/test_bijectors_utils.jl | 68 ++++++++++++++++++++++++++++++++++++ 2 files changed, 85 insertions(+) create mode 100644 src/bijectors_exp.jl create mode 100644 test/test_bijectors_utils.jl diff --git a/src/bijectors_exp.jl b/src/bijectors_exp.jl new file mode 100644 index 0000000..6213440 --- /dev/null +++ b/src/bijectors_exp.jl @@ -0,0 +1,17 @@ +struct Exp <: Bijector +end + +#Functors.@functor Exp + +transform(b::Exp, x) = exp.(x) # note the broadcast +transform(ib::Inverse{<:Exp}, y) = log.(y) + +# `logabsdetjac` +logabsdetjac(b::Exp, x) = sum(x) + +`with_logabsdet_jacobian` +function Bijectors.with_logabsdet_jacobian(b::Exp, x) + return transform(b, x), logabsdetjac(b, x) +end + +is_monotonically_increasing(::Exp) = true diff --git a/test/test_bijectors_utils.jl b/test/test_bijectors_utils.jl new file mode 100644 index 0000000..c666cac --- /dev/null +++ b/test/test_bijectors_utils.jl @@ -0,0 +1,68 @@ +using Test +using HybridVariationalInference +using HybridVariationalInference: HybridVariationalInference as CP + +using Bijectors + +using MLDataDevices +import CUDA, cuDNN +using Zygote + + +x = [0.1, 0.2, 0.3, 0.4] +gdev = gpu_device() +cdev = cpu_device() + +function trans(x, b) + y, logjac = Bijectors.with_logabsdet_jacobian(b, x) + sum(y .+ logjac) +end + +b2 = elementwise(exp) +b2s = Stacked((b2,b2),(1:3,4:4)) +b3 = HybridVariationalInference.Exp() +b3s = Stacked((b3,b3), (1:3,4:4)) + + +y = trans(x, b2) +dy = Zygote.gradient(x -> trans(x,b2), x) + + +@testset "elementwise exp" begin + ys = trans(x,b2s) + @test ys == y + Zygote.gradient(x -> trans(x,b2s), x) +end; + +@testset "Exp" begin + ye = trans(x, b3) + dye = Zygote.gradient(x -> trans(x,b3), x) + @test ye == y + @test dye == dy + ys = trans(x,b3s) + dys = Zygote.gradient(x -> trans(x,b2s), x) + @test dys == dy +end; + + +if gdev isa MLDataDevices.AbstractGPUDevice + xd = gdev(x) + @testset "elementwise exp gpu" begin + ys = trans(xd,b2) + @test ys ≈ y + @test_broken Zygote.gradient(x -> trans(x,b2), xd) + @test_broken Zygote.gradient(x -> trans(x,b2s), xd) + end; + + @testset "Exp" begin + ye = trans(xd, b3) + dye = Zygote.gradient(x -> trans(x,b3), xd) + @test ye ≈ y + @test all(cdev(dye) .≈ dy) + ys = trans(xd,b3s) + dys = Zygote.gradient(x -> trans(x,b3s), xd) + @test ys ≈ y + @test all(cdev(dys) .≈ dy) + end; +end + From 56adfbf7ad4541279553a4643e3c63a54c6050fd Mon Sep 17 00:00:00 2001 From: Thomas Wutzler Date: Wed, 30 Apr 2025 11:04:43 +0200 Subject: [PATCH 3/7] implement broadcasted doubleMM f_PBM move n_batch from solver to problem remove problem responsibilty of putting dataloader to gpu but do it in fitting and predicting functions fix error in solver of properly updating problem.phi_unc --- dev/doubleMM.jl | 33 +++++----- dev/negLogDensity.pdf | Bin 10348 -> 0 bytes dev/r1_density.pdf | Bin 10198 -> 0 bytes dev/ys_density.pdf | Bin 11863 -> 0 bytes src/AbstractHybridProblem.jl | 47 ++++++++++---- src/DoubleMM/f_doubleMM.jl | 60 ++++++++++++++---- src/HybridProblem.jl | 71 +++++++++++---------- src/HybridSolver.jl | 79 +++++++++++++----------- src/HybridVariationalInference.jl | 6 +- src/elbo.jl | 20 +++--- src/gf.jl | 15 ++++- test/runtests.jl | 2 + test/test_HybridProblem.jl | 93 ++++++++++++++++------------ test/test_doubleMM.jl | 99 ++++++++++++++++-------------- test/test_elbo.jl | 18 +++--- test/test_sample_zeta.jl | 2 +- 16 files changed, 331 insertions(+), 214 deletions(-) delete mode 100644 dev/negLogDensity.pdf delete mode 100644 dev/r1_density.pdf delete mode 100644 dev/ys_density.pdf diff --git a/dev/doubleMM.jl b/dev/doubleMM.jl index 8691365..a534183 100644 --- a/dev/doubleMM.jl +++ b/dev/doubleMM.jl @@ -28,35 +28,34 @@ gdev = :use_gpu ∈ scenario ? gpu_device() : identity cdev = gdev isa MLDataDevices.AbstractGPUDevice ? cpu_device() : identity #------ setup synthetic data and training data loader +prob0_ = HybridProblem(DoubleMM.DoubleMMCase(); scenario); (; xM, n_site, θP_true, θMs_true, xP, y_global_true, y_true, y_global_o, y_o, y_unc -) = gen_hybridproblem_synthetic(rng, DoubleMM.DoubleMMCase(); scenario); -#n_site = get_hybridproblem_n_site(DoubleMM.DoubleMMCase(); scenario) +) = gen_hybridproblem_synthetic(rng, prob0_; scenario); +n_site, n_batch = get_hybridproblem_n_site_and_batch(prob0_; scenario) ζP_true, ζMs_true = log.(θP_true), log.(θMs_true) i_sites = 1:n_site -xM_cpu = xM; -xM = xM_cpu |> gdev; -get_train_loader = (; n_batch, kwargs...) -> MLUtils.DataLoader( +n_site, n_batch = get_hybridproblem_n_site_and_batch(prob0_; scenario) +train_dataloader = MLUtils.DataLoader( (xM, xP, y_o, y_unc, 1:n_site); batchsize = n_batch, partial = false) σ_o = exp.(y_unc[:, 1] / 2) - # assign the train_loader, otherwise it eatch time creates another version of synthetic data -prob0 = HVI.update(HybridProblem(DoubleMM.DoubleMMCase(); scenario); get_train_loader) +prob0 = HVI.update(prob0_; train_dataloader); #tmp = HVI.get_hybridproblem_ϕunc(prob0; scenario) #------- pointwise hybrid model fit -solver_point = HybridPointSolver(; alg = OptimizationOptimisers.Adam(0.01), n_batch = 30) +solver_point = HybridPointSolver(; alg = OptimizationOptimisers.Adam(0.01)) #solver_point = HybridPointSolver(; alg = Adam(0.01), n_batch = 30) #solver_point = HybridPointSolver(; alg = Adam(0.01), n_batch = 10) #solver_point = HybridPointSolver(; alg = Adam(), n_batch = 200) -n_batches_in_epoch = n_site ÷ solver_point.n_batch +n_batches_in_epoch = n_site ÷ n_batch n_epoch = 80 (; ϕ, resopt, probo) = solve(prob0, solver_point; scenario, rng, callback = callback_loss(n_batches_in_epoch * 10), maxiters = n_batches_in_epoch * n_epoch); # update the problem with optimized parameters prob0o = probo; -y_pred_global, y_pred, θMs = gf(prob0o, xM, xP; scenario); +y_pred_global, y_pred, θMs = gf(prob0o, scenario); plt = scatterplot(θMs_true[1, :], θMs[1, :]); lineplot!(plt, 0, 1) scatterplot(θMs_true[2, :], θMs[2, :]) @@ -149,10 +148,10 @@ probh = prob0o # start from point optimized to infer uncertainty #probh = prob1o # start from point optimized to infer uncertainty #probh = prob0 # start from no information solver_post = HybridPosteriorSolver(; - alg = OptimizationOptimisers.Adam(0.01), n_batch = min(50, n_site), n_MC = 3) + alg = OptimizationOptimisers.Adam(0.01), n_MC = 3) #solver_point = HybridPointSolver(; alg = Adam(), n_batch = 200) -n_batches_in_epoch = n_site ÷ solver_post.n_batch -n_epoch = 80 +n_batches_in_epoch = n_site ÷ n_batch +n_epoch = 40 (; ϕ, θP, resopt, interpreters, probo) = solve(probh, solver_post; scenario, rng, callback = callback_loss(n_batches_in_epoch * 5), maxiters = n_batches_in_epoch * n_epoch, @@ -213,6 +212,7 @@ end n_sample_pred = 400 (; θ, y, entropy_ζ) = predict_gf(rng, prob2o_indep, xM, xP; scenario, n_sample_pred); (θ2_indep, y2_indep) = (θ, y) + #(θ2_indep, y2_indep) = (θ2, y2) # workaround to use covarK2 when loading failed end () -> begin # otpimize using LUX @@ -246,7 +246,7 @@ exp.(ϕunc_VI.coef_logσ2_logMs[1, :]) # test predicting correct obs-uncertainty of predictive posterior n_sample_pred = 400 -(; θ, y, entropy_ζ) = predict_gf(rng, prob2o, xM, xP; scenario, n_sample_pred); +(; θ, y, entropy_ζ) = predict_gf(rng, prob2o; scenario, n_sample_pred); (θ2, y2) = (θ, y) size(y) # n_obs x n_site, n_sample_pred size(θ) # n_θP + n_site * n_θM x n_sample @@ -506,12 +506,13 @@ chain = sample(model, NUTS(), MCMCThreads(), ceil(Integer,n_sample_NUTS/n_thread using JLD2 fname = "intermediate/doubleMM_chain_zeta_$(last(scenario)).jld2" jldsave(fname, false, IOStream; chain) - chain = load(fname, "chain"; iotype = IOStream) + chain = load(fname, "chain"; iotype = IOStream); end #ζi = first(eachrow(Array(chain))) +f_allsites = get_hybridproblem_PBmodel(prob0; scenario, use_all_sites = true) ζs = mapreduce(ζi -> transposeMs(ζi, intm_PMs_gen, true), hcat, eachrow(Array(chain))); -(; θ, y) = HVI.predict_ζf(ζs, f, xP, trans_PMs_gen, intm_PMs_gen); +(; θ, y) = HVI.predict_ζf(ζs, f_allsites, xP, trans_PMs_gen, intm_PMs_gen); (ζs_hmc, θ_hmc, y_hmc) = (ζs, θ, y); diff --git a/dev/negLogDensity.pdf b/dev/negLogDensity.pdf deleted file mode 100644 index dec12262e1427cc44abb960069189e6610c2ee83..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 10348 zcmdsdbzD^4+O~8H3?0Id0unRBkRsjP(mm205(<(cAt@+Q64D?@N=XVxC?Oz?;D9tp z3Bo%AVw~rB&N;vLd;j_vX7;RmueJ7C_geeDuIt{zqAo4V1?CpOWvN;#J;H?oK|p5< zJ6sVFAdfQ4$=bsP$cHd#;sSv{9$8xl511R`=wR*vlZIJ3TfuNeMRDCd++gO8xZbI2 zDJn6{H3VTB-^e^^>1D|aoL_IYxDWxwFF%BnFT5B_(=U6^OS2)2y}hv$A77|5LeJaJ zC2qMk-Nd<((dQgSez>#vX?yK@&7EUAt9D+|e*XE9XU+3V zTl*;STVIkzA<}uD7H0VOUeKCGYpAhmV~10iM)cn7ax-cYnJ39J`Mx%=d2bl)>8EYy z@YA5v0KS0JH}l1YqLdY+QeKtEo`;SMI(7k;S4e|Y5-CbeKHoJ7yF^a-cHphrXX#L- zD(>9u9wH-M96I|v1sHZ|zP546O9sBsil+KPiEy<7sqE(Z=MsTjb`R39(cQjSgn<_HC&-p~BTEMg z327yd1XuYsE@z0H90)>F-B+!v)EzQ=8CFM3UiB2jfX%WQ^GF_*!H*-bBO_^T#J&w} zP^2lMEc>d$#jH7FWB3#DQph(7bH7OS+4tsmtX?0>+<4Z2lanC=Oo8K%KcB}rEKl3% zH^#?JWcjvst$LsGVG1A)wqog$n%461QH{69D&?fBSbk9*prs#YIcZdzB z67$$&ttj5Z89B-Pn4nlQr4#Rt8f-z;?pn;DHUycvrKR!ivS#ZZ=D3OC-A7)!Pnd#V z%QKp@`#r*wHDFXpVl&7hvafY2SXTL5?R)$!1wgd21 z&vv{lLFs^&a%J|5Lt)1u#}xoqY_7m|9o^%KN7D;dGrYCwxT%C@s;Km+v*{Y>m&w6E zr6ZKDXrkhd(rVi#X8QQz_U`Xdx`pomL{8rCkCQ{mD0}o@$zbx}@vwc4wL*n@(n~N0 zP0BF=Vubfo0;jJMqZ<pc$Ac{C3)g1Ybh-b$Sk8 z!dEn(jX54o1yIZ2j$I|OKS(#JLa&7CGgXk4hWWgn4Au2z_V5SEv$5{g3%CkPQ%`0k zA{35PJ5?W`XB>K=;}TKIzm~X~#W`MPSy7U^3%Eg^Nnm*&O&G(L6WF52X!!sqaKxuu zpjl96qCvv9KFapR18GO)3Kk23v8L+lA-ivVg)Hk#KB(9f19aNP3VIY|fLs=@se&6- zL-=*FodK+mGOs-li(Q!PTzW=$v|N@cnQ4r#INWN4ZJeu^H99a(e zEodhujb6T^yaej)Jaa965aW^i&MOvhY6~4{8XI&+oZHp?qukUGrDL42?9xkzdS&$b zCpyeN9_szULxtqf1dZ*ak?i7(aCCi>qdoaJ3q~OgtFmWKo=WVtO$rmY?CkQdI$c%R zrOI8i6El$ArLxt8JXx+}=S!BHl+PlVgg>C>^_H!Q-r~BoS^V7Y+cSyp?1L%^zHFDs zaSW%ViW!#`4F%pkQrwPYR`|R$SYV?{Hp`4o_9#pB6FN6L`%c$Ytd5{EE;8Pt`}tm+ z%@|U(Kt^#Gr}Ev>E)sM*T)B=t-SiK1>etw+bj&)(_e1J^0k?#3X0ZL3R-<1Cv`%9m zH4|6wxL?_Ni=&aTsgf_;tJ-QR%RRg`5qm=Jy8JGZ_FCplgw=DD`wy++7e>FVDzj!5 z@zzLD**nyYKI04|SjL8pN6OQ%;i1K-CzXZCC|ZN4bT)c=t%lR=InugaSxSMFCZ+du zP-&)yCxKIxQ9K0s<2Wtb6Y7tNqj>yNy%|Ah!qT^7Ux-P7jIIf(V#}5@MQekwS1Z^1 zLS{3o_!tC};zq6orM0@iK+KbF1Ik2~v)3_0T5LHz!$8a|hd8aN@t7v8)I_S$@+2W~ z#G7_#=BsxulZ4clU12npXccje;={g8P8$&#$%rl=8mr9H+O|th9sfs&~i3^q( z(Gn#k<&^k=#==va$wYw4mg!`zC~=WA*kxSq79hZT2Vhj1yFZ$F4yq#q{z^IVZBVVe z9FztdFWd`%z3v4>(0=U-E4rm9>AfeE_!>noTg=|O*@nHad5>0$p@C5vL;p zxq|)huLniX^S0jzF~m$0*XvPJC43pe_o%o0s=^@wD4{7PI%{6lHomx<8<1DWK*jDi zZkT8lw2+gS(pg9Ps4}>O&Mwg)C3l{Xxnw@q2T!$%#&3p|BgNTHWx)6&bqn3DycCZb z5wmySawecmeQjm%fue)ll%*EHO=Dc%tY&aHD>yY+>xvc)nU8vspLC&O7G3MGoRghG z496=XJAk%=og*zepO)6}z2fHDjAfszGr~?>>mTNSeVjBKOrAu}E^Y+c$M^~R@klh| z?5&_!(+7In~==4oDx+`6+YT#*|K|w#n>HXH#p*OHMeunebo0Bm$3QI<4X#6ef zS7x_qvW?pwr$^T0#fo+x77C{d$KJpo5RNjWtP~wTIcyC+0EZp}dx`sCyd6_VL-)6g zr7kCqbZ@yAZ5?c{zs(zfmvm21o2l3UAOOCsU96c5C1<)0u?#(vsznOzRJ$JCqxG#< z8bY40aEUgCryGqbnQzj%sNyX~Wenw_ot#7$^~IN1-#jWy-WgTWh2y{}l3lB3dR_NV zc3^NSIonEjW6#^snPwqB1q(LCwh}Otx?u-TvNFci?-B1sZ&PoG8ds5~(&17Wd%)*q zl$iBzx`L)SMEnIs`DVU`+c(fdMNAesqQg{NcHKS@kNNIIe4V;>u+ZIEMRN-_MA+lC zm(d%5qKCPD2^Z#M^@B4Zz9HEbl42o`2$BWih4LXQAd9>|$=DhFlJWMC)AB%2FQAAB z(t!BmF$W}2R zojiWK;w>%^7YHbDK?;$BK$^fn!9Pu4L@l8U6QV1`|Bq>bK|kjOK(6|Ty*~2DTpI_#0#^D@~Bd}z*gyoa{XbSk^&ZuDg=11gs~a#~II<#aUtOU8TC zQOCxv+T85OQ`4!b?uMd@jlwInc7&y2kz#{Srm0pw@$p>a6}j6oPm#Rv?daf;J2h?Y z3Xz}~7aU_yxLl9_L9{iU1MSl;)|WQ?p}8;K=&OBrx@)n;86*R-z6DdO5G3!6*ArB- zx_kYdtA!rFU*uE;XBE_wa>m?GF6N>8^=d(Bp<%(j+jpk6RxMj%*47us*9SnXM*!5`sG*?#`Qmg77jdp*G2nPwp7Yi{h6g6aY zPu=dZ-wNTO$iQfhBzXfg(kAkg7Z$RU)lKa5J;ek^h}^Hl#@d1&2;Hg6WYnFb_UxsI zem9Zz05bwisbJ{ zmm|MZ)?+>v<1hRWO|NVz-&VMoRnTRsZ;fUapaZp95zl>>CZ9c;A#lpkh@mUrb1%{Z zY~_gVecj819C)EJdd{Q{9oaf$r&^V)Yb1*1zoD>;&b z8>tG?_e9f+_)k+DpKn%LY=#RSl6Fkc$GsU1Uq`PkFyY22>~g%+e6MX$u(v!wLH`c> za?lJ*+G;t|+NDyr$wqCaR5qC@4H{FyaQlE+44o&7pIVr}_H@tPtrFqGrNF7hWQE&e zM?2%#AGlU_m*%2s`;r{B93HGuNa91K_xb9aUC=ydHb`mA`o6&{f?%`nyCwrIKk5=j z;|OF{rKAeH6wFU!vb8RwamZn9Q(26OxJ6nvjd#M&PnJ}yHHgjRH(p(02f8Wh!6>1- z^tCSGwfl=Z-?s(G3IjtnirKlm@$T^SC@Ss8xs2BAjX!+Rm`9(o<}aO5ZN7&^`1JcT z3|X{dk^wEH{Vm+sVX;B97WgT{aW9||yfpijenejPiWKjg+=U*7FN@2vUv~?ArZ%@! zp>EI|UCiJk3thXwdi^|43meM>O8NYzS^(y9JF)0(*S$i+^iPw zea9;B5w{roG5=>O%{F&j3;<{?H>$zAiv2oGIo24(qDY4msc1$T$$$}4%UX8(=<4?a zqCz3+tLfNTx8veddFoTFeT_p*iYVwvb?GO!?yo*KiZInKI(%B9YgHCh)-swevU_c! z!YA(qQNZA2O1CUQb-0D1L~KU=YEg>t;`b+SJl$8F1iVY$@;O4s;Ni4_cP&PJD%OHl zoMPuhGQ+or+~+&2TfP(1qc1BTx{3hF7l#m$kbbtN7LU<5gqS zA8?Hn@;YTS4DTJa_sdp3%@q!|*Ufphy%RYLW79ucr)Hs*YfZ*8l@j$9PFl2h2zWTi zGwE&BlL(hSKD?s^>5;}ZpX>|7e?MHD^Wx=pe#=a`nno~}gjgy#2C6iK6_Y5dT9?CZ zs$nS*%e6+xI6>QBQ63wy4zbGTc!l3?^Pb9rIG3)kRp#9wIezq~9%g|yef{ez;+zL& z%$JR?z1jROS0q}ZPx8J zQyqsAEhb+Njs}mWZ!`*HI#$Myn+wf^4tBaY2aZAp_xU>nhGlb3OrL|A=cPKkN#MHYk18ie|gCG&9Ym2SXC~DB4xhvcY3(1OXS}$hKtk~ z0`b8@0s=od!v%T!EoBDg|F@h0L4|(i%#7m)oPt64of|$8wDnRLA^h)|L!h1%myfc` zZ(LGIy=)u8_72Ti?0TuNE!HMa|GvIz>iqDH+vY~Ft@K+HZ5-bzje=lJta#h_$>4BB z1%Z#d)ja{%2{cG;)@xv)1$h4CDol}tlSUsZyQdPK4r~H*5<3jtVmn>p>dUkS3Qnj$ z5h@eu>PJ*Rx-639M>W3!3goQtnGu?#A9akW7VChMOB6Pg*0Nl`t~YPSKrn}z#r{^8 zf4T-W=%^+SYS=&RPpAOxc~;2@4ws$irHsAt@+yBX=Albh6v~cVKmCJ1s#ujYDzwk^ zLZ|XRd(_+1{WLwJ#bpYn$-OM)7P{?Z?lLAHCuLp58mZQ?(2Q5tFqb2ixhAcx%h>0}wu8HOvJ|3G)}1-}Li{!B9b<@mT`}4|4}wO9>}y2N)1^ zW|45Wgdw=O5LDo-h#WMs%mo1po*AUfUF2c5);1nz-JF@UJYbHxK*96&$b$Su0Xc0% zf)|FkCXxkf{zlmU-bd&C{Xu=tOp@mAurtl@-#H6O>ZNSV-GIE1v*ya?KgwXp58D67 zdH}&cr3a6U6Cwp|>tqe&xoPVp;pA?6VL0pPrmdBSjXT1nNV)O%FJ4|=AVdHuT=@9- z|MEV!{(VcNRR99{hkpOL&BYny^$`7oeq9+kj&m#IN59CUAYz2Sj|~Jwjv3kRSAArk zNXywaLF{5@9Vc64ZXO6ePu>5fypU79*zFfW>!(Qj+x-3?`bJ7Y8o9E_`vml-%yHBJ zD+NhVHtlHD&9*7Kgp2M<+1clF85S;Uy^vu>z1qgata1Qp$<@EZke6Aa6Q21zDW3)Y z8TbuHWhAW;cPdR$scmXUbYHY#pv^B4o&8xNH~C@H=$=o~QFOx8w#FlhgE`xhRT66E z5Oq)f#dwd`sy*m&-S?FB2MJK8PrOz0$%Jc$a9CC*-r{xzH)hTGo6HR0F)P~&0*t5j zyf~uK(u+1f-KrW_&3jZ(WNd7?LT!56?%32h=>fVnqump|c1H4CpHt&>V}cdeOO&H1 zl=mjlf)qCVhDeO9@n%=ZaLu8nD1#UvtSFgjxBR9T%|NFJTXxMPwR?@M>x%ms`}jHO@rhc zQM=%87^^hj7(bGL$q$b(4L1^(U}ByADA|k!S$_jU$i?~@#dt5G7zDWm{$*VuzXbGu zdTwJJJ7CuiB-_dU3i!C{B;aXY%uQ>z`enByM8rXTM-W=Iwv-*heo8#3q`k|M&ty@z^DWs}GSU5qN6#FrJ5Oj$ z4Q17JYUG%o(sGaMF_T;+HEzp&wqs;q=5hM%OQx)ISy*8`tVXkBP={T&@6&@BE4-kl%A9znt_ZD*PTW{skQV8>b2PoCxpZa%X7CNCmrYYRc6Zk$Pu0pxJq`eW- zylm}OfJs>er{S?^>PC(2(`%W_0Li{5NsY^ge#Idnu~p0|?FnoatTf%^5Sgor5iE_0 z_uTZUMn=BKevSL)R#o{i-SQGB1x;OCC2+5G{Y8MS*M!aE_OUcRdK$_%QNj+T{m4gf<^kVtBtpFr~<7&wLNuMXR)4jWz zUsnHt0RtMf(AeJ}u8=A`Yja2W)j&#$`zr}U*HcAU0=l1q^InjQ!la4~PfI+?3;Ab_ zb1NV`(vK}m@cvoX#G^f`p6OhHj$o(G^(dwym{*ciL;)A*c-inktb1Rl{|d^(mOT%L zc$1TUjV^S>M)=sfkSpADhVblH0+JudrogjWRh8viUUpD|fCt=-zFFkEP%b`pHF-{% z1xtLbJI;O}O%_dI4)g_2*bZ1ExArL_w*{rf=z!fSSZCPbsi3-RY|*EkYQkHqN@ADn zsW*^3$r<~cmWI$8Ytg#!k(QY6>Ov`mc-%``59dCvQmjPWxcV_sB~~?-H&0<|(X0GO zf};UiQDao; z=BB78us%jsCsXFEH>1oXeGU5?z_}2Mi2F+r^ixAZkYw-zaDMPI$Sb*_f+jU7|&r6wW&&zZ?{J}+2N?-Hs$ zq?Dj3LM?T|1DLmt>kU?BMn&B$nsg&hf5%L>&o>XzTf>Hni)p$OQ!e|d=MUZHX>E?c2qfQfRcCja_gAiHjvAO^s_byGQyOG9 zYK)}}NsDvo^qn=>M=$wC;a41?@Yb0adt+UT)pD8PSo=`BwrCR8+4H!rV3OJ21~2dG z;lM~hfAO`^Z=;HJ2@Mq!;cj>gDmV8GI(yKlB5Owc;4cePO|{AbHvY!V7m5YpW(eOe zH~)Ru{ztJVt*omlDgReomvnZp`UxQa4@wX55Q2*rB!E74z`gfWZk*AJNrb|JeQc~ZTJz3j(j6ayn;f16>Vp+;P*8_Ji7Z89L~vx+Fu@*@&01_ zsuxLxv&Ur+-t(RMAAY)}6Q^V!0wN9FKZ&Sk*TGG>??ukU2B)DO>n_u-PtGXElxX(Rol=v&s0`e36fkAiaJ_}&3MmInU2 zTn$0B)4Z*?t;=DeFRxgug@^*;IaveBXHiPp{F#bzg(NF*8VYnHvkzNp(bY`Mt~D7; z?b9t@-qs1IxJ8@E*+%qv+%gOA6|gGnw$x^P4My|5nrwsTxH?teK4_+?y}1N-pS51> zwaImReRBYhTTx8z_5t4M7oFBfo}sCUIl{@e&z2~z89R;tHLjd_`9d!U2qG&WOT2$w zAw5mltux18sO%qi<==(!Gp=Rp3qwfzvnB}g2?j!dkn>=m>Fn$Q%y!>|zz#yZzWGBFCngsalF&5)JaxjN6 z2C*}+-wd6&H|onPaj+)^Cg!E3T0q8M9p&y}?&jg`hVTIx1mc6>varag$>RPGL6xa< diff --git a/dev/r1_density.pdf b/dev/r1_density.pdf deleted file mode 100644 index a05680ac22b85258e9e4971d14548cba2179c1f4..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 10198 zcmdsdbzD^48ZM=DD=9E^3Cs+QDBaQxA`(OAFi1*whk$elD1wx9BaKK2A{`O}(jo{5 zcLu~b-*>)qe)rtJF2l_1wbr}V-s^qWexK)Adl*$Dr8z*H{8)^Y8^y<1JOCiT-q;FD zSQx-{7h-1)vjFhItJJUn005V?r40-Mg@4-^!61?l6MIt#)}1?8PB19M$QH{rZA0%( z)V*5V>DiK+H(ztMWJmmo2uwku1OQQzex8_)7VNfYx0a1i52C)Ul^bdYs!MCTuH2%2 zyH5){h#2qhZpXY(OG7=KkvR$B=(DGOrmUX!?O-+dTk4q6>HJtIH}%E_(7SZ(36__> zq4y1LhGU#2==)%|M+W8%nNXeH8&{x}hQm1!FZW;6%>V*qkRVsbUtGo zenrmk3o>E6S$i<^X{_lnlO9`|O44)Q03i`!rP+pN{DY~@`R}lUYlgxnHb=W0uYjy6F?&dO&~32QD!qXKuM#1exBD<1M4AJ-o)G%advEp`;5=RkPN%*U3R4zC~1 zXXel$E8_z~VMRmRsl9av%{2#4UX)=__u~rOypA`am-rZ5#l2CTRN4HRB8&AP%)4sl zcs{hOWfO0RQ7c=7a5;}OcGzV$kbHg%nB9kO|E%6sj44T*y6ur$QHHgh!!2#)k)t=L zgyp8Ww0yWBQrHFDE>C@rs^vJ%%nkh*{JmO8WCxmfHaKP5ABqh>v#D7FsS_j(8>p~% z`57FWqJ-Yfj7z-M&48O>d(BiQstn51b*0`oqJh(7YNUYqS*Tx_;tE-rRQ2`NyFfNG zl$4J_2_Yh;52`S!^gC-!f+gwM)oMp9?XC?5P$VHyYLJdDe$0=JP3XH`c$2L+DQ4q) z0M83EUDiD`vChQUO$TBtU9`b&H`ZllCjD5!$~=;CB~6=*cC}(P@-Z$@K9%@bR-du_|{*NfbJ)k{?S9RkAPlGm|N{Qa|QW zTwX4UZnII{tiysLEJ)Z#Zsyl|mf%E!lgdiAfzr8)tk4 zN$3|l6eK2N_iJ<;s2}h9C`IYp2WX(01`zqh;{{s|=U~1V?2H=pv%`gC*Z35w4~B-~ z^<6_ArZ={2pY=@a$Wdf)kytE9nuv z-B_r|8r%VWq)-C17ChanI?!!rpJzIFYjJG0k;sRkZy5|;$oT-wt`Dn?tc2V~XB?lB z#q9P5t7>YaMnU=HUg>blS9%nLn&6eTyNbqj4JY_!vXG->$)!vCt@F?`#aGaNdfT(x=lp$9Sy!2TKAcG}gjrLR{d>8J|K zIeHlVs)Iy75KKKkkzV3wJ2T9ON~}lrp3N4Vh>ma`?l^6Bj#KtWHcG2C=Aw{B#@D2u zX4ozw-!z?8I`yF=6^s$-pJkrHcSwLScnjCQ>yl7m*VcAd-CxYac zcBfQC^m>D^z%goS2EPAyWcSEti^@$CJv6n;s$`3^JzBR8=+ol~C+EBf{9bHM;djcS z4OT!!WEq8AZXNgs)-8RU^K#I_N8dy*tgIyQIbOl6kdd@q#MjF8q(3tg{m?^WveG4X z`at*1LF9>C+r-0NbSZw1cCk}3YxpI0PPa71H49Ywi&uJmgsS-|OBi+M$6;YgiMi_6$fjVQ_&{x%Z$}n5u zB=Y3_Hzx;g{wiS=5!Xs`6!blRsYXNcc6$a2{04li{Zmu=Q?VzvcKmNN9u>Fa+x2nK zgk+99S@8YV6H!T%oH_C#x(}ZTXTJM{Kttkv?@5-l=}7E134W|QK<^KsFrWyOGa+Fm z@*C}tp7Il#2)0$(69S6_`W=7c_(b|0`uaAVu8>SWGv9lmG$(p*L6kdP##yt;G+)*l zhkH*{Wvs!reaaa_#LLm_0WQl2uGGmigq57e6Uid($92Tb5uFd1y>{-ei=JRq-x_!G zrpWZHqW|{P&O^N-R)pP6L1s6&eX=h7gmzc@tvl+D`xct4f|XM!nUtfx8ZzUI$t07F zGYbWYqr2sMX0&Kluk}igmANL;DwlDF^{phD4YtIUk*yB5>7DFpO>5sClJT0daZ}GirE3TUOfTI)kG=pYNs1}2S&hC>m&8nRl)lj|frXNtclNeB5c^>?=~sE<3|kZ zfl&9&S=<1z6v8Ww8Qm`< zS7!35sc}2aM6Q}Q3`NV-c4x(c$;h)Se%w)lt`XRvsx4>!v5{&n@MJ^!xK*Vk-g%>; zP+3la`epSYo3N&P(TERgRtu2Ab}No$_34n-`WLj^=k#7hYL9tt-4oQp@PGL#Pjbv& zk8}TJrnju4wYt6S+Dpa8%7SXX{rCAZ6C?Vm?$+un`l&}Ul?8_Tv7aXjnlkVAXEIe< z%YwJ6b58}mg=gZ+mCeW4NhhTW1jz$FuM~b$(9R$9Ag3>?PJX{LVf-#ef%J7>z{;kD zLBpvS9#h4w625U_$M}1iebn#1 z7n9cBdqewSY@`a?P+~pYrRjOa^bv);VZ)(R-p1OPics@+=cPhJ_P}ntfEKIUT=us7 zhn&-xLT#%Mo1*ZzAWKsA4zX>ID35gtVXjg2wx-SXNat+5jgxXpZZvb3#eS}?Z?;uw zlv}F0s9AB7H6*0f){e7{ZZ0*xqerTSF4ea%M>gk!B z(|K*=w3TuCLa8U%%ry2})AXuNg)FBB!f$z1o=sDGc)rKGCa9^RwV*4wvY|9v?=Kx8 zgPT}57&sf3Y0=PRT1(OCll!ighQIA8-pbyU&JDhG8U}W11WR z+cpY5axYlII(vxhdO8$Udw9{LGSaOYYMTgE#;~otZ>b+Il}l*zre$AG$53ue+PD=t zmERes#hdQ7yg9Zmz1xQLGA!e}#>TAjKppCeikNoV0IuivcU>;_2i8H%3CmWoI<+nu zV5906sdGouQ!h(D0HZ#>(Tqy^NIfv~eg93<`YxYP@hS9(D~k>X1q)(l`h&2-j}TNB z!EX`YaE1%!2lF8e5Jm2vO!%DpO1Z*h)M0Sm3lJ7YRKWkZjDYY0fJ;IQ00MBC7{Nc` zZ1}PbblC>v1cA?+&kugH0RLveduhS@V~h(sE@^u^*l%Yv;{b910sNPY9B~n593Xg= zz@Jq>4iG@_vI;&F{Qt+aK)|2#0&r=#IY0nh%EngeFx&H4{-FFI&}CmPB>?EWMcNYT z1OxE=IR0lHz^x9P3tJ9bLP5GJfNbzrkDm|!!{wl#4hQ}6_)qNtLg_Ntn-@S5P_8E9yLa*ClYSqhYF&kdmlts5pr^j?s}a zbdCE%AhWE@=mKtT=Q5D9t(p5f zl8G_zKvOTkW;=zp*!hI;ui_y7+YF)=Jb3Xc_)**Q+72^GK_`FjZkxsV`+@X33~-tDdZcLy)fGDNiFn zQ+u3p(fEf7Rbi0$u|>@C^Zl!1!gE|VXx}@8jJgEGe^Pk=UFE5usQy4>c_e`<=05ct zwd{5PU(~yj=UXP*hN_43?UGC2$;USYx&VbbAWv8O=EM51Y)o%@p|b zT{G%f-ACiHPEmI9Wbo6aKr*AB*iTKvYd*kU@{TE-XrR5<3d?w`#Vq{Q*Q`W*lSuzL zmbbn3OKw+;g|A9=xuUmiu2{7}*;UYq1-&+UTlE)vb|%m-v*dm4(P5qq9S-w)%QTVhcXr1SOp zP&m?D4HH$0%={p~Ga)&>k8&W{PGEs|ZD4*a#(g`Or6A{4N}^}1i9C>P@Uf+d-+s}P z@dxP~>bdt;sO6yb+c4fkyv3r~S}iP@7zfK`p|7T9bj`I^|dQx`2S=)o!QIGU}OHXDdE zg&j`}zQ^%)3UH(`8nT%?N`4Sju;72^i2(G@nbBsf&X*ND+yvZYA}fW09Jygi?kJgV zJ{g%ggVE(jv^x)GUzO>epC>Y{0;jyj}=*pg zQAtg4%DJ6wC#ie$h!xf^yYX#Or0=N*)_2{+$3gKbPQYlt!Hk@=4#aO4V8Cp!0|Rk z6PBwvRcOcXqpV{0XxXWd8aIg+54s{U=ic;JqfwS@f!lrs{4MJ$_o75Jjjw3QN7g4C zJLi9+8>x|hMRVhMWnQ_#{64xQQl*{F?aKbZIjhk6dHo7Jk{jBZcZJAr9z6W4Vo@?Q zBV>S=)7TiO9$Y=DCoQ(zsJl3bcg4MSNZ4V6<4`ylgRV;~Is-c{b8ZEma0D-UKL$%oLV_tQYvuBs*dO8nwK)3&e2krSY3+@Qk%$q$rDAYl#=nkr_o4Hvy_ z&G{Sbx!mD`V+V+rTkt0?UebWycD_J-|A;^LfSi{a$(&?ds*#&OcjI7WxelTjbS2xp-=u%W-W z!2Gz7>V~NPb^H1g?RTG~V*P@hLk0Hx`*Liv5AOym;IP|d?0ddOM@H)K7$hsY)wX;r zP|9Nbz5Z#o6EnNJb@Ki*C(1EftYat7m()f{2iKEw%_75Tt_KaJ zt5!h9#jNi#53YKfj}ix`Mbt%zRX~%9Aju>{?=_bz7&wy^jK7X7#uh0Ir*`#O%1r_J z7fZ$s`nS#I3BUEjF2*hlLsjrUo){b`P+#Mv`LwCNcnYBL9iV}H6kz<(44ql)jfCDVR0|0UJasxmh z9sz*Cc?UU|k&UH^n4P%|1OPm*5pyzuzzLin5C3@)anXn}2N)!9ULj%RAPccHw}71w zb6%wmgV<^T1TOj`3bL04#I)gCun_n;5mZU-H&*xeKDrq14}NrBC2r&dIVYC?owE?k zPr||o3g8BxcfV`&qYMIFOzY2h06;(aDVLNTe3RJH&K$s{ZD}WF=VWuFvj73SK=}K|Il=r0ZVd#0!Qda)D8SAAv-QHB8&L!A zgRn)^oG%CXYPHw2vqbEY!!r%M@9zNvAqKr%)t9XMCu9HHo%TP-u7sQ!2>koLYiR_(>k|B zuWs;3jL8x0z=?!!-tj?@;U>5UQq(9cUoHxX6lN0*Gs2H%3z_K0r;Ek_mx93MJD8sY zYHsB(GUvJH=PoChCTg&;B6DaWjlV#u(n7W=K^ZZ&fdctKK;L`qI1fo;-nMstdi7Ce z`XsY9XL@?t+o`@@u&&-az}bIwGn3$y2Izfes1TEqS@fW6Yb(p@RQ#~H)w79)-}~lvM>fjao%-Nvn8_l?Nx0{b)-UoA^PN`vEJifSu+eK|9PRn@9xSZiFluofs>=@*K zu0OF?%2JEkET7$_YKDk*mRunpgWy-}99R*&I;$QGn?*6K7i@AaIpiEXGYIj=Cu|bV zYB5F%T!7Ic*;^$5{w6gcB6FFVz&w|!35-}5{}M#-FX{B39+w!zv|6_ z9ES}aV)zajKtBj2ZXn;#SHcaqxwaBhow*Q%$5bkIbT)93*aRz~3dA|`f4Twk$6q8{ zsVS8rbHK@66xxp3gLg!3FYREohjMj+?dh{#3CP*VsLk=t41OyQtr z33L0|r0(qCU<0v5K;1cx!DoT!@h>=W9`*}yuMN-#!f}8{5T17hKmY@S->~iS-un~( z|D5-E9{j(U@Fg4pxGv|wA$vZ|1yY+D!4VgpW)T>5?tirlp#PQF&r`tfi5+~&FMs2a z;D3IE`bT0%%;aBoF&?Q|%RUynmTo`RC24{G;eLuucuHVFwy&0+o+bV}m}Vb9ZIWd~ zOOa6_H(1KKJcr05MlUGkN?^FOs+CT1ph(=X*R^Y64<55T*|Uuw?yIIC8;VYzpYhlf zJOE)is)gjRyo$S70P6UbxS!~`-t=}=*W;$@`_Ji;j&uPm3uuMZlKUeCd;2evh+S6( z$dp$;_U}n8cJjCxrf}L%;kFIVW@$Ginv#~Y8~G4P(Wtt6#-6oGB>pxz zv2pdAM^QjPOeI4~dpwIV6J-wxSc+LbjImMP2dYaxI=UwPCH4qfSuv4ef&@%KRS{M6 z{rY;d)Z6mmq{Wl=59z!#lw@5ILN>+Yn_oriDPK{ve6P9l*dBSwaqbJfPE6mDgOQbAwdhc&9$tLyr`l zuwzioo><|jZK}^q9BB4aAxE`-h1ta2KkZlTMw4$``7jtj#YwHd^0ba${QdQ5 z(42Z@McK|nD;`|$P|im89Flz=4qjGe*_%>}COGPz_8verMl~Tev_)sg9!NN^Ce?$3 zoQL8{2gIyklXj1*ocxi2@mr0w``k53#8yepIs-{lH)HZ|P~w}TFPXncP#1CkyjTn- z9DgXOgBgsa5F-{lu5v<9fnJ5~%Jm|xuEE)lB9U=cg$q(I>*IlE`v5n~`S|zXbga*jy?) z`28jDOHY6=WybF}=r6|lpFz9Ut{p%SEH7m4gDI_uQ{ z$N6gfE-U&xq$WYH<*cN6w^AmQNBBxnWMY(sD8+WzL`JX2b%ralA|kX4r=WxxgACLM zybEBR4U8305jBTMqu{o9DDey-QkK1#5%6J9+Zd#%*m3U^owS0GY}G^M`8yY`*Ue#9 zd`Wg3Raos7yzZ7aM+{8US9I7}DGamTYm6ZaNRM^sbf42ZKzn&a>QNLX_r{(OV{22L z$z+xGMB|%AP2m)zv-e4D!4!j+1@<%MZv&&={YC6!M`QA}@eSpZAy8~uMeVP8oxP~! z;nkxaE3aOp8LF3fZ~cv%FSQ5U&0yZ2vJLz@!TkkS|5mn1-qlhQm;EcIi`(0n{sfQz zL9@a4oH^jC5~0}y;bIG(lK)?5HZbo64*eIJ?Hm|>M~c69(8YNFq}l$nr~W~({n)^V zW7@As0M27|0mQ!m@qZF*|CJQbHVR1fj_DF#po|qt@2-!! z!m5dta`?=ebi9r2X5O{-Y5i5&_BpH$St5mlN&cZxH5AmnkBT#tL4K|~*kFAz?lnfV z9?tzRl%XX)EJ@z=fUiiB_Kym8e96#o@0e@!S^~3X=7Gv1@>Rs7Sa-}Dh>C#y>;Ao@A&<0tLpa-X8^#HZye2izedDq`|3+kWtJ z+VkyW2f=&gnT&&!#qvRh5D5oK`N?>WGefB#P|hq`f=fbhcT&)byeDj)%SaFj^m~u-SY4~ zhH8Dmrk%9w|7%`3_wuDy;1@s`AWGbSjgXES#O&PhmpU4*(9Y@R-^AgZYgxKO0NfX` zf)h^=01N04*gF!!_!x+mnLW+)vad_W))Qg{}j-(UUY Y1T%udT%m9m2nqsuuoxMol%=u$4?QrS%>V!Z diff --git a/dev/ys_density.pdf b/dev/ys_density.pdf deleted file mode 100644 index 4ca94db5bdb275c65ce7e817f29bf447c203e72e..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 11863 zcmdsdbzD?U|2Cb{NDE8pg48aJbR#0&9ZSO!D@ZreEhY`por07gU4oQ=AR#HOG`zbg z>iyjJ^ZcIA`+NU-IqdG4IWu!+zBA{0UDtOGvx>AVJBWiHhq?M|$q5b*fC~UMv&9h> z25{bkI9kE10le@mH5>o{z$t5E4}-YCf9*|S5NU`x)B=JdDvIL@bAgyT;5>Z!)u6Y4 zx$|+#vi4#bPaf`9CmPAsc*c`t#Quh zcP=2CZOM3;m}N9}bFlALg3}Q zeBK`w-%<`4{;nL$<6WX(=wscdOu#!oX40WjM!@%Jevve#*{kwloku&!54Q|hPX>-s zgzWJiq>l4%J(+#u_0`BVhkq+s`{HuLrNtOr?qD?6^cJ z>YlMw{l|l5M|S_y-JDQ0PPvPtO?=zQm%X2lP5jP#4YTxSa8Lk7UysLYYm3$#3UeNL zLqE_MTr4M&I)ROTI_J>Z?IVYcmZ#inS1e?osAAE?xFUzPWElu%}7sB8=y)KAL$`9l! z$d4Cm>KK|Gl0;d3l$fGtu5 zr@xeG zIw7KGgE_p$LnD$l0NcbX86d=mr{DsOlr}A#V31MG4ZzouOm5Nqirdff$a*4kO4vba zCXOjXm4lF3?rdn?Qferzj2l~}nn~fsfdO`ElCqSUrb3(vr%+@_U7273vj=eyxYFMi z>FsO9!1>JzrUfi)yc@T07zN1q;uiy@qp0!U;>T|!E#`&CnWT10g^x(*>6P#+B-E+^ zST|%CtT3as-V2RCN~<(0nCK4XBPa)dFDLuR9-Ye+7sF22`o!o?p~4QCzghiAA5s|X zh6D_N*#grwnBnzhfD)S^RR8%H498tk-=+rdHKjF$?3&4MU<%v6u5($hfyi|)$XnCL$c^cOFMkQ&^S)(f0qdxz$Kzsx3O?~O0HGA~ek@s;pZLSsG zP#do_qwf=lf5adE#t8cp50C{W6wN0w>7;j?t?#4fzT2D=^&V}Yzc1xU1q1Ef!2}%^ zlo22+HSpDhu`W?G8Oz4*^=H?##v}Uhik><5ceQ1=F}&{Y(#sCK0c>DiLmkf5qReEC z778lH>_<+HF);SPVD$oNk&TBawZsa>HS|q@qm}stJbeH#4f@#UgyFi#I;Nx}rx2nC zz!r)L*1Jzb0rp6|c67|FW?cB>19YjaWQ@CGUZq*D(0WQZkUnS`vAl)at71tOA`w}2 z&aVnKyG1g6mDi8W&^qxqR3Vy^8=U??5qKxNnA68XTzYoUJJ0`gmF#{_MHb_tie87X zA$LVJ(#6`wdw0UU2ty7D1&tJZtkg9g)qd8)qNz_W&?31UKA&rKNGqL%$_UEp2D2@3 z<0j9U_)1#`3_KEeaV>x{O*QzbrVm9PkX*rI{R0!w#*?-6-hG*pl^8bcN!p@eS+w%y zJ~U8)kN5HCAraIV2-T1IYu6UdWO*M&>!I z=je?Oc_ay>HbDua?4>}-B_#=^kdJWyWlDo2=?C60*6ZH9v@d{6?{+fPmK;>#*MZ2; zxL|26-9&X2=4=}E&;-i{xdV$qLtw>mVqB?DRQXI>^B8O@EI&B7I(4c;;TDrQwh=KC ztvOD$UGJFLi(@M@%Kfo7@w5poxJ^2Dv{f;;hkZZn#&6-u7zW(daq^CXg$(cOuxytP zlf;VbO#JQZ%y_;$(IT*ru9Rq^0sVR#@&gYPzIN-Cz3B2!Xri94WXYS2O_7MI^jf>js> ziLl#a9(-gY>S+m+8bW8o?@Cx--l@^_f$m_PCg4p~I|fx8M=^CXX;v$LxoskhB8DGj z9ba1}9*Rs{Id|P4<33}!Xx99-VC(LApZI*L4B}0q@9dKb*EN}A%f@@rmNDnIvAn*2 z3lhZg>Ii32vZG&9ry(I~Xj<{nh-$li*V;|rPQb|jMM8RKAy{=(X)8mrDx4yR)3c#8 zwX#howvO&RW)#gZM1~k$FoJ|n18)^0M_@_}bp@@_Sx?gHZ4)2&d}^zbLb*r13K_b1 z=ruF~0QMt!3+oP4`HBYNgSY>v96ON6I+8#o|dAlb>_y2tsI#~jSnN#L`GT< zJzAo34J|sG;wDx>T7d%@R?t>HL&1%EO>)=~u`~i$Y0_^~gEsResa>gHoKYFK>T`RZ zhupq3o0vy^%RcMi9+3Yil}W-JV|K)B<|b8eY-aj3>fk}^IEkIs1E6V$*p3K@@#R?s zNfY)<kWGdZy%T4pwBto)R13B}nX0iq@bpwFUL*L_CPJD(M^T-a`PU~;43 z-;WxMEq$~RJQ%&(#4U=};(yk0*yZ8oj?~E&zhQ^59)MZ{9<)*u$9PfjA=@97)deDS24pwJ51{)?m z#X&AsB5Eey1Rn-EQ}enMm;#{pOlTG?n0||cI3kTEE`N(n00Yxv=Ci~NAOOGXju#)d z0N}t;*wRJZzPEm3069Et92Hp@gvIC|sh9b9B?Y-&h&K8jM%Np{!L-v@P_9AherW_3 zb}PX-W~Y!k=vdCKVXonQ^J97kB+oLjiM8tjL1&D60&jmO(@_?pP!q~W z88I#7aAZAiVv_c3J-3b%C)=|w#3SyC@Vbwp33M<{M4c!7wrwF8FOKobxmE$Ck`C{| ztpE>#if&}ZvNwkJnacV!Qkk4c15vSyAk=W1w|u!N$A@FB9D|VpH(=}&i@8PruG#J_ z49rET0^dddNZN2nmUj@^?&Nck^ruX5{o)y5K|mGW^ac%~T&y{<_A{n9hYm_r8m72n zdFMOw>D3+fQG;7P6mq5EZp0HGuE(&I^(x;Y1CC`KyD|i%9K8y6Bb&&Txr4jmHk^od zS7lHzo%do%8YLXMp*x%HJ+I!u8NGu|n!~s-wKHC!F8Z0T(n>|N%`Jgqq0(x=>`BlL zHhqKvN8&D}5Wk*v(w6k|nA3r|$Ehu|>8|`~En<^XiH#!v;p>9LtHIaU{iuF(9xoMC z=iaNm&iQ7hZZ{;o)y(kXwDeos4oaTUa(`*>7Qgs-0PWEh|9b9bpzkq0v%~8L)H3kl|pas8?E+iCTm?^d_n z&{jX$=he>tN^baQ?498_rm?CT0k(~WwM&`O&`3b|jvjRfV&#rKJFCOGJ4JuJ+$H01?dw_pI<+IOfvtC4H3@?S z7t(~E!}w=f%R`q>Su`C@U$hww<-U3-8H5Qsq0DlYp6`4v(e(^nnv+(|6M(+viR)GxKU%tM%*It;GM0m8zF4EP_XDHl8e;FOX8fB>B4rtojL zTEDsmy1K^20RmrMzC8J(1o)2dH$c{6x`}ua}i%6mQau|1dt8?Zt=_Jf0`V`_p4{WE&eMu0z$+~{(W!+F@E(2 zX0Xd)AO-`H`b)qhMAJhfY5clh8!;URs%=M2giv}QuOx+K)N^c@G$C`aVbVKR9FB0D z7=1u=elST8Ig@R+nLLx?ZOLh7>^veF;ot&N#usL6@;GBvpwaI>`>i%xC&^<9BDsFv zvpxq+C!4!o2cL2oIfODtR(+G>$_0d{QbA&Esjqi(s|$0fLacSR3@>Y_Zonh+$O8)Hfs%fY zu0YM&Xw9pCIGp+=CC+(Lv2ufQ)lSPEd<@~3tz&MPHN)pAep$}paAql)5FXv0@DLmu z9FIkOm&`2BlUV9Y0$N2_6w_wAcexn5j0m-P5&9aKx=4jD|iefaW zn+xOy7nds(j7iiCB281R^miiTJz9)?VnjG~8dk)os0fyO<-mG^l8|U@SnvivNhz3v zb1vh9nl%z1xc*Ex-7>5i#Xa|Mbg&V3UwLQ<-=1Q$l4Fd49qc6K>v?TEqasV0hDHH` zA0vI&O6b0}5t80DFSq1H+|RDIDI4pH}s5<POQZW%34o-&r zgM*vtwSu-LU*&x%_;2I8qpaFH^((&l-#=fZ`(iV z3GKBAt?@7(dVEqthk8uA!3q4toNXFmrz+E_o4S%_pm|}JgU;ht;%V-s0~)L8Uko}> zJ{oL2MbZJ@r2Dw%pAjsAI@8obuA+!+)b;d`+E@IH9-o9RBs{i>bQrUWd`;i~bzyO! zQg4vb^kZctDwOc0+f>1KUTlnh5g<5Ew0gbrxYzi6(+(zpT<<-EQvD_be+nrK}AvbiWklaI7&FnPYZVypA)zA#636$sa z)Sw*%infSsItzSdntruC(=YozjN-cpx@KAvBEcfGhVxqJV;`9cim8@==T zysj{DW1lj#7qq-HHu%=Lew&35b~etY(oj5;Xs5OrXUA~N3 ztk|Cqr#OF9LR_hgjdxZb*N5EeRdJA#WU=!QC4+L*(32-@6XF+j=hno1ai|lwm#AAR z9u8=LJiN*0S$gNSwZ3#jN45~4761{6S-FDF_QZG%(qiR z7PEy#@8QnfuxuE4`^p5#JY(>hcb!vr6Zv?C>Fw*RLr*n~klKY6325_7gcjGUDW7KF zk27N#1;}NKr>-611MQ!-0d2Hi(@N+Sia{~;aNgY!o1hQgQ#rJ4or=TcL9MZDv!#E@ zPPgJ~-vuP-q<*M5;_V~I&?6U6?I^N0TvE`dmlzYe)wbk$81Oc01*51NO((HZm{K!U z``+kXsifr(-N-t=vZmZek@CuPqSY5KTab=Z(H*mo&y&vsJ0A&AlkTsrC8}mpFJx`g z6Ek~FA!Ey|eq1Br#Lfx6b^B~(vRymE_r1>9lx4b-N_9CyYHFi_<`|pj>M6r^;tp+T zpsKEQSI<;(8Y)b%h;>jz<_VLW#r(tVK=cRORyfV2c{Yc^TAO~&VR$0qrJL+`kmN!F z2to>7Z64W+2La4|K1F^c*~O=R@0=Xn(T+Eigm&&yXN>g~Z|kgp=3qz0r&dLtQO5YF z%;~+(jk0HbCIe4H5){e}VrM3;SbZ1jP3r!5_|P ze(j7Hswq0Q0SISYBXnOycJVRkJZ4kWtl!C`8~ut*L!qX~D?6ck_#N^T__)4o8F>r! z+Pi$O&C@!`9a^iA1sunHzIm{5mW03pWIr>(M6U9a>=E=8eHcUdiG!~5+M5?LBjP*p zWKVb<%e*S~DCy9UNTKh|XwJ~xjB`_SkL|F;1Tu^>NF8v=3@f}S->!4flj^?Nve%;< zcSy7FiyCAHbP1O{RJ5$A3vM-`IYpOIPiE&Cs77f#tG#KL{4(W;Hq&~oQ<}Y3L)sNa znVaxMF_jZru>)-C(`NEcv3DgRx9ed%`lWc(wQGX(@MFB@bk@ zPdS9+$n6DeW}jozhi9BZ9d~b`yD_UIns#gJc3*(T(FytrUa{vt%phV`uJ|$!|E~!| zkWxf8Vgi5gN(6V6hPax$*f_zUEVduK@6xiS4mS2403f^|08tWfSrl+tmgPsA_NG>@SIiXNzN94dAwVC%&dVhT zU(K0dQUBNVu9q;5=IpJ{$z4M06UFW(R`= zE;FP|o#Y`lR@Sh~W-hbTVGsu`fWVLXh=lxA0x@j(P96k)PXwb?`-6P{vyFbV_ZJ#_ znI&oJ3c2L}|DCZA1YOG7)CIr|zN~)F^k*6b{z>8gG9LiYFB!op;|Sktv~jcoaBAB) zN;taOTxDD~q-|pXvv!4>6wHs18UK86b94W^&dbXO_@~c5*Zbe(BJRTV`(9kX&#!t4 z27`b0R{-9>KYPUmK=hawUI$V3@(TPMaUk+87YckGLp2?35ZkZ-(2rf(KXeqr6j$s0 zN;dtHT>rS0`xmt$r67%%Tf~#X^;gm2pbAn9l%QzZQ-3|zcF!qH^qZ8eT^_rB!K!+h z3=_(&Hby3;LvU-3t~Y&dda-6$`m5wT=H)MdV{E0-lt!GHl)H*;Gkc;3qHhMY*V~JLd6uU@_=tC0oHQA2Q9wa|DNtJ9?GZlhxhZY z3)PNbGp|qe;&cTyW`6cHnjOStx@RLmY&i4XjWrrIwQzgEvU>7**NZmc64W-{u5jTnHdIxLu!ISun_#XMVv_*K z=KPJ%Rnj=iNl6m8>NjApFlv$NhMd!Fy(c13-$2Lc8`Q@PA4x!DM<*EiTS+Sru?}D4 zEQb8dfAC|7$@=BTxUc*e7_kKYWnRI*u)az!55(ok6 zp>s^QsDG=e?3M%z+bg!djuj0rK|Sj3&l$qVwR$4!O^5`T?WPYAXk9%4(gpb@aGX+bw1gt+BXAlHwl;7M1&3D&c+3lOiKhrJrR-t>aY5h#P{9HYgf=jb zUzgM$I62uv91!q$iGA=j9T;_eb1D8=%hx#||DrxMzifF|Pm!U}*RUCjRM1 z|J~Rxz2yJJps!#IzAIUtZ$RzfP(+ho$9ug_5YOC`sL@Z(0my}e(BnkL*&mnQRua=79eeAQ(pS(5& z4?sB1YGIkc$^^E2(7U7LZ^_>4P5rC7UTmsk8!6Jx48g#8v_e|xZy)maz7;$ndAK}C zp}f2>uqU(F!Q)~4oC7*_&4a{QFsF~SFX9t7u-$@qT2{gF!$dSyqw4O(?etY*$^K`_ zjjQnYfP#bLs+pd*Cj!lGQg@SqWp3S#U~atY@1jdNI=UviAAjspUG*{59Es~Ws*1Q$ z!1uSCWqvkpQ`S%0CsKIns3|(5gzQTuH@}NFP*+m5oYjdYL6Mi7=k^(O;`(+*UEVgh zT?_D)=yZ6K{;;Xg(zF`uSR<0+z)s#C(>257RV#r#9xxJv@pcPpcBpEv-wTsH%xLkC zj!dfAJY>vXq7*dVsXdQ8?eG%BSiB?jU@X}euirQtbKm3 z2S-#6Jeb&r>E*IcED;o6OXf@ZD!H9@=vH27!x%j}Pt;=Lz(AP7OQAVy?|b!w&!4;2 zOXxda+=V2fc_~1@2g)c+D_L{4CLou|XEc6Z=fWlZ*t)VjFz1|f@}08hQx5;TAjgi) zD8@pFTe4I{K0DXxs{Wx^_krfX4dj^C?=brWqw@jPF0{Li%M(MvH#lhZm$T~mCC9GM zfacVztIBuWYRx!U zQ}R%uzk^ujZ_@8^R#5sFn)PeEG~%wiPh$JzfzIHQX|}jmG}HuEm`hd#kJZIIHx^64 zM3ZjPI@phJ6yqe~Csj@ft1xRYA95DFEUeFq#5)t^m&GA3iD%R%@d03&1U4;7)RyQ%@`(s39V%hUe-aE%6uOgyY^-{?jap!#7 z9DyD{w&SeA;yCYnucA3>aE7t!oujSdFpEKB97S+SyiFfMyCSPmZ+*?P3qkaQLx5tmi@4ilaQ!y3hf=jQY{avr4 z2bD6icGPRRvf!n$da2*mKbZMSv%t*^=KXEve|FpdR_;mP(^8X^|GV6igxXvD0+9cO z)`M>;uyb?qBeb3%T%^HW^8X912j;!(+kd0=E`i}sq`0D&h;{g%aP5aS`>pl49-yw)`pTpl#3@sHeLHU1 z6y<#M<9nJ+>sKjc3)N3T9m{cf8VEVErG)i)o->^1v9; zl{bqH?>J4q1+Qx)6ujLsGIJFblz(;uI+UlN$P;IuO0JjW7Vfx0em;6T f_doubleMM(θ) # omit x_site drivers +# par_templates = get_hybridproblem_par_templates(prob; scenario) +# intθ, θFix = setup_PBMpar_interpreter(par_templates.θP, par_templates.θM, θall) +# let θFix = gdev(θFix), intθ = get_concrete(intθ) +# function f_doubleMM_with_global(θP::AbstractVector, θMs::AbstractMatrix, xP) +# pred_sites = applyf(f_doubleMM, θMs, θP, θFix, xP, intθ) +# pred_global = eltype(pred_sites)[] +# return pred_global, pred_sites +# end +# end +# end + function HVI.get_hybridproblem_PBmodel(prob::DoubleMMCase; scenario::NTuple = (), + use_all_sites = false, gdev = :f_on_gpu ∈ scenario ? gpu_device() : identity, ) + n_site, n_batch = get_hybridproblem_n_site_and_batch(prob; scenario) + n_site_batch = use_all_sites ? n_site : n_batch #fsite = (θ, x_site) -> f_doubleMM(θ) # omit x_site drivers par_templates = get_hybridproblem_par_templates(prob; scenario) - intθ, θFix = setup_PBMpar_interpreter(par_templates.θP, par_templates.θM, θall) - let θFix = gdev(θFix), intθ = get_concrete(intθ) - function f_doubleMM_with_global(θP::AbstractVector, θMs::AbstractMatrix, x) - pred_sites = applyf(f_doubleMM, θMs, θP, θFix, x, intθ) + intθ1, θFix1 = setup_PBMpar_interpreter(par_templates.θP, par_templates.θM, θall) + θFix = repeat(θFix1', n_site_batch) + intθ = get_concrete(ComponentArrayInterpreter((n_site_batch,), intθ1)) + isP = repeat(axes(par_templates.θP,1)', n_site_batch) + let θFix = θFix, θFix_dev = gdev(θFix), intθ = get_concrete(intθ), isP=isP, n_site_batch=n_site_batch + function f_doubleMM_with_global(θP::AbstractVector, θMs::AbstractMatrix, xP) + @assert length(xP) == n_site_batch + @assert size(θMs,2) == n_site_batch + # convert vector of tuples to tuple of matricesByRows + # need to supply xP as vectorOfTuples to work with DataLoader + # k = first(keys(xP[1])) + xPM = (; zip(keys(xP[1]), map(keys(xP[1])) do k + #stack(map(r -> r[k], xP))' + stack(map(r -> r[k], xP); dims = 1) + end)...) + #xPM = map(transpose, xPM1) + # make sure the same order of columns as in intθ + θFixd = (θP isa GPUArraysCore.AbstractGPUVector) ? θFix_dev : θFix + θ = hcat(CA.getdata(θP[isP]), CA.getdata(θMs)', θFixd) + pred_sites = f_doubleMM(θ, xPM, intθ)' pred_global = eltype(pred_sites)[] return pred_global, pred_sites end @@ -157,6 +191,7 @@ function HVI.get_hybridproblem_PBmodel(prob::DoubleMMCase; scenario::NTuple = () end end + function HVI.get_hybridproblem_neg_logden_obs(::DoubleMMCase; scenario::NTuple = ()) neg_logden_indep_normal end @@ -173,25 +208,28 @@ const xP_S2 = Float32[1.0, 3.0, 4.0, 5.0, 5.0, 5.0, 5.0, 5.0] # const xP_S2 = Float32[1.0, 2.0, 3.0, 4.0, 5.0, 5.0, 5.0, 5.0] HVI.get_hybridproblem_n_covar(prob::DoubleMMCase; scenario) = 5 -function HVI.get_hybridproblem_n_site(prob::DoubleMMCase; scenario) +function HVI.get_hybridproblem_n_site_and_batch(prob::DoubleMMCase; scenario) + n_batch = 20 + n_site = 800 if (:few_sites ∈ scenario) - return(100) + n_site = 100 elseif (:sites20 ∈ scenario) - return(20) + n_site = 20 end - 800 + (n_site, n_batch) end function HVI.get_hybridproblem_train_dataloader(prob::DoubleMMCase; scenario = (), - n_batch, rng::AbstractRNG = StableRNG(111), kwargs... + rng::AbstractRNG = StableRNG(111), kwargs... ) + n_site, n_batch = get_hybridproblem_n_site_and_batch(prob; scenario) construct_dataloader_from_synthetic(rng, prob; scenario, n_batch, kwargs...) end function HVI.gen_hybridproblem_synthetic(rng::AbstractRNG, prob::DoubleMMCase; scenario = ()) n_covar_pc = 2 - n_site = get_hybridproblem_n_site(prob; scenario) + n_site, n_batch = get_hybridproblem_n_site_and_batch(prob; scenario) n_covar = get_hybridproblem_n_covar(prob; scenario) n_θM = length(θM) FloatType = get_hybridproblem_float_type(prob; scenario) @@ -201,7 +239,7 @@ function HVI.gen_hybridproblem_synthetic(rng::AbstractRNG, prob::DoubleMMCase; int_θMs_sites = ComponentArrayInterpreter(θM, (n_site,)) # normalize to be distributed around the prescribed true values θMs_true = int_θMs_sites(scale_centered_at(θMs_true0, θM, FloatType(0.1))) - f = get_hybridproblem_PBmodel(prob; scenario, gdev=identity) + f = get_hybridproblem_PBmodel(prob; scenario, gdev=identity, use_all_sites = true) xP = fill((; S1 = xP_S1, S2 = xP_S2), n_site) θP = par_templates.θP y_global_true, y_true = f(θP, θMs_true, xP) diff --git a/src/HybridProblem.jl b/src/HybridProblem.jl index 625c3e3..5be6a91 100644 --- a/src/HybridProblem.jl +++ b/src/HybridProblem.jl @@ -1,70 +1,71 @@ struct HybridProblem <: AbstractHybridProblem θP::Any θM::Any - f::Any + f_batch::Any + f_allsites::Any g::Any ϕg::Any ϕunc::Any priors::Any py::Any - transP::Any transM::Any + transP::Any cor_ends::Any # = (P=(1,),M=(1,)) - get_train_loader::Any + train_dataloader::Any n_covar::Int n_site::Int + n_batch::Int pbm_covars::NTuple - # inner constructor to constrain the types + #inner constructor to constrain the types function HybridProblem( θP::CA.ComponentVector, θM::CA.ComponentVector, g::AbstractModelApplicator, ϕg::AbstractVector, ϕunc::CA.ComponentVector, - f::Function, + f_batch::Function, + f_allsites::Function, priors::AbstractDict, py::Function, transM::Union{Function, Bijectors.Transform}, transP::Union{Function, Bijectors.Transform}, # return a function that constructs the trainloader based on n_batch - get_train_loader::Function, + train_dataloader::MLUtils.DataLoader, n_covar::Int, n_site::Int, + n_batch::Int, cor_ends::NamedTuple = (P = [length(θP)], M = [length(θM)]), - pbm_covars::NTuple{N,Symbol} = (), + pbm_covars::NTuple{N,Symbol} = () ) where N new( - θP, θM, f, g, ϕg, ϕunc, priors, py, transM, transP, cor_ends, get_train_loader, - n_covar, n_site, pbm_covars) + θP, θM, f_batch, f_allsites, g, ϕg, ϕunc, priors, py, transM, transP, cor_ends, + train_dataloader, n_covar, n_site, n_batch, pbm_covars) end end function HybridProblem(θP::CA.ComponentVector, θM::CA.ComponentVector, # note no ϕg argument and g_chain unconstrained - g_chain, f::Function, + g_chain, f_batch::Function, args...; rng = Random.default_rng(), kwargs...) # dispatches on type of g_chain g, ϕg = construct_ChainsApplicator(rng, g_chain, eltype(θM)) - HybridProblem(θP, θM, g, ϕg, f, args...; kwargs...) + HybridProblem(θP, θM, g, ϕg, f_batch, args...; kwargs...) end function HybridProblem(prob::AbstractHybridProblem; scenario = ()) (; θP, θM) = get_hybridproblem_par_templates(prob; scenario) g, ϕg = get_hybridproblem_MLapplicator(prob; scenario) ϕunc = get_hybridproblem_ϕunc(prob; scenario) - f = get_hybridproblem_PBmodel(prob; scenario) + f_batch = get_hybridproblem_PBmodel(prob; scenario, use_all_sites = false) + f_allsites = get_hybridproblem_PBmodel(prob; scenario, use_all_sites = true) py = get_hybridproblem_neg_logden_obs(prob; scenario) (; transP, transM) = get_hybridproblem_transforms(prob; scenario) - get_train_loader = let prob = prob, scenario = scenario - function inner_get_train_loader(;kwargs...) - get_hybridproblem_train_dataloader(prob; scenario, kwargs...) - end - end + train_dataloader = get_hybridproblem_train_dataloader(prob; scenario) cor_ends = get_hybridproblem_cor_ends(prob; scenario) pbm_covars = get_hybridproblem_pbmpar_covars(prob; scenario) priors = get_hybridproblem_priors(prob; scenario) n_covar = get_hybridproblem_n_covar(prob; scenario) - n_site = get_hybridproblem_n_site(prob; scenario) - HybridProblem(θP, θM, g, ϕg, ϕunc, f, priors, py, transP, transM, get_train_loader, - n_covar, n_site, cor_ends, pbm_covars) + n_site, n_batch = get_hybridproblem_n_site_and_batch(prob; scenario) + HybridProblem(θP, θM, g, ϕg, ϕunc, f_batch, f_allsites, priors, py, transM, transP, train_dataloader, + n_covar, n_site, n_batch, cor_ends, pbm_covars) end function update(prob::HybridProblem; @@ -73,19 +74,23 @@ function update(prob::HybridProblem; g::AbstractModelApplicator = prob.g, ϕg::AbstractVector = prob.ϕg, ϕunc::CA.ComponentVector = prob.ϕunc, - f::Function = prob.f, + f_batch::Function = prob.f_batch, + f_allsites::Function = prob.f_allsites, priors::AbstractDict = prob.priors, py::Function = prob.py, - transM::Union{Function, Bijectors.Transform} = prob.transM, - transP::Union{Function, Bijectors.Transform} = prob.transP, + # transM::Union{Function, Bijectors.Transform} = prob.transM, + # transP::Union{Function, Bijectors.Transform} = prob.transP, + transM = prob.transM, + transP = prob.transP, cor_ends::NamedTuple = prob.cor_ends, pbm_covars::NTuple{N,Symbol} = prob.pbm_covars, - get_train_loader::Function = prob.get_train_loader, + train_dataloader::MLUtils.DataLoader = prob.train_dataloader, n_covar::Integer = prob.n_covar, - n_site::Integer = prob.n_site + n_site::Integer = prob.n_site, + n_batch::Integer = prob.n_batch, ) where N - HybridProblem(θP, θM, g, ϕg, ϕunc, f, priors, py, transP, transM, get_train_loader, - n_covar, n_site, cor_ends, pbm_covars) + HybridProblem(θP, θM, g, ϕg, ϕunc, f_batch, f_allsites, priors, py, transM, transP, + train_dataloader, n_covar, n_site, n_batch, cor_ends, pbm_covars) end function get_hybridproblem_par_templates(prob::HybridProblem; scenario::NTuple = ()) @@ -110,16 +115,16 @@ end # (; n_covar=prob.n_covar, n_batch=prob.n_batch, n_θM, n_θP) # end -function get_hybridproblem_PBmodel(prob::HybridProblem; scenario::NTuple = ()) - prob.f +function get_hybridproblem_PBmodel(prob::HybridProblem; scenario::NTuple = (), use_all_sites=false) + use_all_sites ? prob.f_allsites : prob.f_batch end function get_hybridproblem_MLapplicator(prob::HybridProblem; scenario::NTuple = ()) prob.g, prob.ϕg end -function get_hybridproblem_train_dataloader(prob::HybridProblem; kwargs...) - return prob.get_train_loader(;kwargs...) +function get_hybridproblem_train_dataloader(prob::HybridProblem; scenario = ()) + prob.train_dataloader end function get_hybridproblem_cor_ends(prob::HybridProblem; scenario = ()) @@ -131,8 +136,8 @@ end function get_hybridproblem_n_covar(prob::HybridProblem; scenario = ()) prob.n_covar end -function get_hybridproblem_n_site(prob::HybridProblem; scenario = ()) - prob.n_site +function get_hybridproblem_n_site_and_batch(prob::HybridProblem; scenario = ()) + prob.n_site, prob.n_batch end function get_hybridproblem_priors(prob::HybridProblem; scenario = ()) diff --git a/src/HybridSolver.jl b/src/HybridSolver.jl index 6e00147..60eb1dc 100644 --- a/src/HybridSolver.jl +++ b/src/HybridSolver.jl @@ -2,11 +2,9 @@ abstract type AbstractHybridSolver end struct HybridPointSolver{A} <: AbstractHybridSolver alg::A - n_batch::Int end -HybridPointSolver(; alg, n_batch = 10) = HybridPointSolver(alg, n_batch) -#HybridPointSolver(; alg = Adam(0.02), n_batch = 10) = HybridPointSolver(alg,n_batch) +HybridPointSolver(; alg) = HybridPointSolver(alg) function CommonSolve.solve(prob::AbstractHybridProblem, solver::HybridPointSolver; scenario, rng = Random.default_rng(), @@ -21,32 +19,32 @@ function CommonSolve.solve(prob::AbstractHybridProblem, solver::HybridPointSolve ϕg = 1:length(ϕg0), ϕP = par_templates.θP)) #ϕ0_cpu = vcat(ϕg0, par_templates.θP .* FT(0.9)) # slightly disturb θP_true ϕ0_cpu = vcat(ϕg0, apply_preserve_axes(inverse(transP),par_templates.θP)) + train_loader = get_hybridproblem_train_dataloader(prob; scenario) if gdev isa MLDataDevices.AbstractGPUDevice ϕ0_dev = gdev(ϕ0_cpu) g_dev = gdev(g) + train_loader_dev = gdev_hybridproblem_dataloader(train_loader; scenario, gdev) else ϕ0_dev = ϕ0_cpu g_dev = g + train_loader_dev = train_loader end - train_loader = get_hybridproblem_train_dataloader( - prob; scenario, n_batch = solver.n_batch) - f = get_hybridproblem_PBmodel(prob; scenario) + f = get_hybridproblem_PBmodel(prob; scenario, use_all_sites = false) y_global_o = FT[] # TODO pbm_covars = get_hybridproblem_pbmpar_covars(prob; scenario) #intP = ComponentArrayInterpreter(par_templates.θP) loss_gf = get_loss_gf(g_dev, transM, transP, f, y_global_o, intϕ; cdev, pbm_covars) # call loss function once - l1 = loss_gf(ϕ0_dev, first(train_loader)...)[1] + l1 = loss_gf(ϕ0_dev, first(train_loader_dev)...)[1] # and gradient - # xMg, xP, y_o, y_unc = first(train_loader) + # xMg, xP, y_o, y_unc = first(train_loader_dev) # gr1 = Zygote.gradient( # p -> loss_gf(p, xMg, xP, y_o, y_unc)[1], # ϕ0_dev) - # data1 = first(train_loader) - # Zygote.gradient(ϕ0_dev -> loss_gf(ϕ0_dev, data1...)[1], ϕ0_dev) +165 # Zygote.gradient(ϕ0_dev -> loss_gf(ϕ0_dev, data1...)[1], ϕ0_dev) optf = Optimization.OptimizationFunction((ϕ, data) -> loss_gf(ϕ, data...)[1], Optimization.AutoZygote()) - optprob = OptimizationProblem(optf, CA.getdata(ϕ0_dev), train_loader) + optprob = OptimizationProblem(optf, CA.getdata(ϕ0_dev), train_loader_dev) res = Optimization.solve(optprob, solver.alg; kwargs...) ϕ = intϕ(res.u) θP = cpu_ca(apply_preserve_axes(transP, cpu_ca(ϕ).ϕP)) @@ -56,19 +54,17 @@ end struct HybridPosteriorSolver{A} <: AbstractHybridSolver alg::A - n_batch::Int n_MC::Int n_MC_cap::Int end -function HybridPosteriorSolver(; alg, n_batch = 10, n_MC = 12, n_MC_cap = n_MC) - HybridPosteriorSolver(alg, n_batch, n_MC, n_MC_cap) +function HybridPosteriorSolver(; alg, n_MC = 12, n_MC_cap = n_MC) + HybridPosteriorSolver(alg, n_MC, n_MC_cap) end function update(solver::HybridPosteriorSolver; alg = solver.alg, - n_batch = solver.n_batch, n_MC = solver.n_MC, n_MC_cap = n_MC) - HybridPosteriorSolver(alg, n_batch, n_MC, n_MC_cap) + HybridPosteriorSolver(alg, n_MC, n_MC_cap) end function CommonSolve.solve(prob::AbstractHybridProblem, solver::HybridPosteriorSolver; @@ -84,34 +80,37 @@ function CommonSolve.solve(prob::AbstractHybridProblem, solver::HybridPosteriorS ϕunc0 = get_hybridproblem_ϕunc(prob; scenario) (; transP, transM) = get_hybridproblem_transforms(prob; scenario) pbm_covars = get_hybridproblem_pbmpar_covars(prob; scenario) + n_site, n_batch = get_hybridproblem_n_site_and_batch(prob; scenario) (; ϕ, transPMs_batch, interpreters, get_transPMs, get_ca_int_PMs) = init_hybrid_params( - θP, θM, cor_ends, ϕg0, solver.n_batch; transP, transM, ϕunc0) + θP, θM, cor_ends, ϕg0, n_batch; transP, transM, ϕunc0) + train_loader = get_hybridproblem_train_dataloader(prob; scenario) if gdev isa MLDataDevices.AbstractGPUDevice ϕ0_dev = gdev(ϕ) g_dev = gdev(g) # zygote fails if gdev is a CPUDevice, although should be non-op + train_loader_dev = gdev_hybridproblem_dataloader(train_loader; scenario, gdev) else ϕ0_dev = ϕ g_dev = g + train_loader_dev = train_loader end - train_loader = get_hybridproblem_train_dataloader(prob; scenario, solver.n_batch) - f = get_hybridproblem_PBmodel(prob; scenario) + f = get_hybridproblem_PBmodel(prob; scenario, use_all_sites = false) py = get_hybridproblem_neg_logden_obs(prob; scenario) priors_θ_mean = construct_priors_θ_mean( prob, ϕ0_dev.ϕg, keys(θM), θP, θmean_quant, g_dev, transM, transP; - scenario, get_ca_int_PMs, cdev, pbm_covars) + scenario, get_ca_int_PMs, gdev, cdev, pbm_covars) y_global_o = Float32[] # TODO loss_elbo = get_loss_elbo( g_dev, transPMs_batch, f, py, y_global_o, interpreters; solver.n_MC, solver.n_MC_cap, cor_ends, priors_θ_mean, cdev, pbm_covars, θP) # test loss function once - l0 = loss_elbo(ϕ0_dev, rng, first(train_loader)...) + l0 = loss_elbo(ϕ0_dev, rng, first(train_loader_dev)...) optf = Optimization.OptimizationFunction((ϕ, data) -> loss_elbo(ϕ, rng, data...)[1], Optimization.AutoZygote()) - optprob = OptimizationProblem(optf, CA.getdata(ϕ0_dev), train_loader) + optprob = OptimizationProblem(optf, CA.getdata(ϕ0_dev), train_loader_dev) res = Optimization.solve(optprob, solver.alg; kwargs...) ϕc = interpreters.μP_ϕg_unc(res.u) θP = cpu_ca(apply_preserve_axes(transP, ϕc.μP)) - probo = update(prob; ϕg = cpu_ca(ϕ).ϕg, θP = θP, ϕunc = cpu_ca(ϕ).unc); + probo = update(prob; ϕg = cpu_ca(ϕc).ϕg, θP = θP, ϕunc = cpu_ca(ϕc).unc); (; ϕ = ϕc, θP, resopt = res, interpreters, probo) end @@ -157,14 +156,15 @@ end """ Compute the components of the elbo for given initial conditions of the problems -for the first batch of the trainloader, whose `n_batch` defaults to all sites. +for the first batch of the trainloader. """ function compute_elbo_components( prob::AbstractHybridProblem, solver::HybridPosteriorSolver; scenario, rng = Random.default_rng(), gdev = gpu_device(), θmean_quant = 0.0, - n_batch = get_hybridproblem_n_site(prob; scenario), + use_all_sites = false, kwargs...) + n_site, n_batch = get_hybridproblem_n_site_and_batch(prob; scenario) par_templates = get_hybridproblem_par_templates(prob; scenario) (; θP, θM) = par_templates cor_ends = get_hybridproblem_cor_ends(prob; scenario) @@ -172,21 +172,24 @@ function compute_elbo_components( ϕunc0 = get_hybridproblem_ϕunc(prob; scenario) (; transP, transM) = get_hybridproblem_transforms(prob; scenario) (; ϕ, transPMs_batch, interpreters, get_transPMs, get_ca_int_PMs) = init_hybrid_params( - θP, θM, cor_ends, ϕg0, n_batch; transP, transM, ϕunc0) + θP, θM, cor_ends, ϕg0, n_batch; transP, transM, ϕunc0) + train_loader = get_hybridproblem_train_dataloader(prob; scenario) if gdev isa MLDataDevices.AbstractGPUDevice ϕ0_dev = gdev(ϕ) g_dev = gdev(g) # zygote fails if gdev is a CPUDevice, although should be non-op + train_loader_dev = gdev_hybridproblem_dataloader(train_loader; scenario, gdev) else ϕ0_dev = ϕ g_dev = g + train_loader_dev = train_loader end - train_loader = get_hybridproblem_train_dataloader(prob; scenario, n_batch) - f = get_hybridproblem_PBmodel(prob; scenario) + f = get_hybridproblem_PBmodel(prob; scenario, use_all_sites) py = get_hybridproblem_neg_logden_obs(prob; scenario) priors_θ_mean = construct_priors_θ_mean( prob, ϕ0_dev.ϕg, keys(θM), θP, θmean_quant, g_dev, transM; - scenario, get_ca_int_PMs) - xM, xP, y_o, y_unc, i_sites = first(train_loader) + scenario, get_ca_int_PMs, gdev, cdev, pbm_covars) + # TODO replace train_loader.data by proper function that pulls all the data + xM, xP, y_o, y_unc, i_sites = use_all_sites ? train_loader_dev.data : first(train_loader_dev) neg_elbo_gtf_components( rng, ϕ0_dev, g_dev, transPMs_batch, f, py, xM, xP, y_o, y_unc, i_sites, interpreters; solver.n_MC, solver.n_MC_cap, cor_ends, priors_θ_mean) @@ -197,17 +200,21 @@ In order to let mean of θ stay close to initial point parameter estimates construct a prior on mean θ to a Normal around initial prediction. """ function construct_priors_θ_mean(prob, ϕg, keysθM, θP, θmean_quant, g_dev, transM, transP; - scenario, get_ca_int_PMs, cdev, pbm_covars) + scenario, get_ca_int_PMs, gdev, cdev, pbm_covars) iszero(θmean_quant) ? [] : begin - n_site = get_hybridproblem_n_site(prob; scenario) - all_loader = get_hybridproblem_train_dataloader(prob; scenario, n_batch = n_site) - xM_all = first(all_loader)[1] - #Main.@infiltrate_main + n_site, n_batch = get_hybridproblem_n_site_and_batch(prob; scenario) + # all_loader = MLUtils.DataLoader( + # get_hybridproblem_train_dataloader(prob; scenario).data, batchsize = n_site) + # xM_all = first(all_loader)[1] + is_gpu = :use_gpu ∈ scenario + xM_all_cpu = get_hybridproblem_train_dataloader(prob; scenario).data[1] + xM_all = is_gpu ? gdev(xM_all_cpu) : xM_all_cpu ζP = apply_preserve_axes(inverse(transP), θP) pbm_covar_indices = get_pbm_covar_indices(θP, pbm_covars) xMP_all = _append_each_covars(xM_all, CA.getdata(ζP), pbm_covar_indices) - θMs = gtrans(g_dev, transM, xMP_all, CA.getdata(ϕg); cdev) + #Main.@infiltrate_main + θMs = gtrans(g_dev, transM, xMP_all, CA.getdata(ϕg); cdev = cpu_device()) priors_dict = get_hybridproblem_priors(prob; scenario) priorsP = [priors_dict[k] for k in keys(θP)] priors_θP_mean = map(priorsP, θP) do priorsP, θPi diff --git a/src/HybridVariationalInference.jl b/src/HybridVariationalInference.jl index e8ad83a..a7fa4e9 100644 --- a/src/HybridVariationalInference.jl +++ b/src/HybridVariationalInference.jl @@ -20,6 +20,9 @@ using Distributions, DistributionFits using StaticArrays: StaticArrays as SA using Functors +#export Exp +include("bijectors_exp.jl") + export ComponentArrayInterpreter, flatten1, get_concrete, get_positions include("ComponentArrayInterpreter.jl") @@ -38,13 +41,14 @@ export AbstractHybridProblem, get_hybridproblem_MLapplicator, get_hybridproblem_ get_hybridproblem_train_dataloader, get_hybridproblem_neg_logden_obs, get_hybridproblem_n_covar, - get_hybridproblem_n_site, + get_hybridproblem_n_site_and_batch, get_hybridproblem_cor_ends, get_hybridproblem_priors, get_hybridproblem_pbmpar_covars, #update, gen_cov_pred, construct_dataloader_from_synthetic, + gdev_hybridproblem_dataloader, setup_PBMpar_interpreter include("AbstractHybridProblem.jl") diff --git a/src/elbo.jl b/src/elbo.jl index bc9b3ad..dfa69ea 100644 --- a/src/elbo.jl +++ b/src/elbo.jl @@ -63,6 +63,8 @@ function neg_elbo_gtf_components(rng, ϕ::AbstractVector, g, transPMs, f, py, #nLmean_θP = map((d, θi) -> -logpdf(d, θi), CA.getdata(priors_θ_mean.P), mean_θP) #workaround for Zygote failing on `priors_θ_mean.P` iθ = CA.ComponentArray(1:length(priors_θ_mean), CA.getaxes(priors_θ_mean)) + # need to apply different dist to each entry in θP and mean_θMs -> @allowscalar + # but does not work with Zygote nLmean_θP = map((d, θi) -> -logpdf(d, θi), priors_θ_mean[CA.getdata(iθ.P)], mean_θP) θMss = map(θ -> interpreters.PMs(θ).Ms, θs) |> stack mean_θMs = mean(θMss; dims = (3))[:, :, 1] @@ -131,9 +133,10 @@ Prediction function for hybrid model. Returns an NamedTuple with entries - `y`: Array `(n_obs, n_site, n_sample_pred)` of model predictions. """ function predict_gf(rng, prob::AbstractHybridProblem; scenario, kwargs...) - n_batch = get_hybridproblem_n_site(prob; scenario) - data = first(get_hybridproblem_train_dataloader(prob; scenario, n_batch)) - predict_gf(rng, prob, data[1], data[2]; scenario, kwargs...) + dl = get_hybridproblem_train_dataloader(prob; scenario) + dl_dev = gdev_hybridproblem_dataloader(dl; scenario) + xM, xP = dl_dev.data[1:2] + predict_gf(rng, prob, xM, xP; scenario, kwargs...) end function predict_gf(rng, prob::AbstractHybridProblem, xM::AbstractMatrix, xP; scenario, @@ -141,19 +144,22 @@ function predict_gf(rng, prob::AbstractHybridProblem, xM::AbstractMatrix, xP; gdev = :use_gpu ∈ scenario ? gpu_device() : identity, cdev = gdev isa MLDataDevices.AbstractGPUDevice ? cpu_device() : identity ) - n_site = length(xP) - @assert size(xM, 2) == n_site + n_site, n_batch = get_hybridproblem_n_site_and_batch(prob; scenario) + is_predict_batch = (n_batch == length(xP)) + n_site_pred = is_predict_batch ? n_batch : n_site + @assert length(xP) == n_site_pred + @assert size(xM, 2) == n_site_pred + f = get_hybridproblem_PBmodel(prob; scenario, use_all_sites = !is_predict_batch) par_templates = get_hybridproblem_par_templates(prob; scenario) (; θP, θM) = par_templates cor_ends = get_hybridproblem_cor_ends(prob; scenario) g, ϕg0 = get_hybridproblem_MLapplicator(prob; scenario) ϕunc0 = get_hybridproblem_ϕunc(prob; scenario) (; transP, transM) = get_hybridproblem_transforms(prob; scenario) - f = get_hybridproblem_PBmodel(prob; scenario) pbm_covars = get_hybridproblem_pbmpar_covars(prob; scenario) pbm_covar_indices = get_pbm_covar_indices(θP, pbm_covars) (; ϕ, transPMs_batch, interpreters, get_transPMs, get_ca_int_PMs) = init_hybrid_params( - θP, θM, cor_ends, ϕg0, n_site; transP, transM, ϕunc0) + θP, θM, cor_ends, ϕg0, n_site_pred; transP, transM, ϕunc0) g_dev, ϕ_dev = gdev(g), gdev(ϕ) predict_gf(rng, g_dev, f, ϕ_dev, xM, xP, interpreters; get_transPMs, get_ca_int_PMs, n_sample_pred, cdev, cor_ends, pbm_covar_indices) diff --git a/src/gf.jl b/src/gf.jl index 81e129a..58868b5 100644 --- a/src/gf.jl +++ b/src/gf.jl @@ -19,13 +19,24 @@ end """ composition f ∘ transM ∘ g: mechanistic model after machine learning parameter prediction """ -function gf(prob::AbstractHybridProblem, xM, xP, args...; +function gf(prob::AbstractHybridProblem, args...; scenario = (), kwargs...) + train_loader = get_hybridproblem_train_dataloader(prob; scenario) + train_loader_dev = gdev_hybridproblem_dataloader(train_loader; scenario) + xM, xP = train_loader_dev.data[1:2] + gf(prob, xM, xP, args...; kwargs...) +end +function gf(prob::AbstractHybridProblem, xM::AbstractMatrix, xP::AbstractVector, args...; scenario = (), gdev = :use_gpu ∈ scenario ? gpu_device() : identity, cdev = gdev isa MLDataDevices.AbstractGPUDevice ? cpu_device() : identity, kwargs...) g, ϕg = get_hybridproblem_MLapplicator(prob; scenario) - f = get_hybridproblem_PBmodel(prob; scenario) + n_site, n_batch = get_hybridproblem_n_site_and_batch(prob; scenario) + is_predict_batch = (n_batch == length(xP)) + n_site_pred = is_predict_batch ? n_batch : n_site + @assert length(xP) == n_site_pred + @assert size(xM, 2) == n_site_pred + f = get_hybridproblem_PBmodel(prob; scenario, use_all_sites = !is_predict_batch) (; θP, θM) = get_hybridproblem_par_templates(prob; scenario) (; transP, transM) = get_hybridproblem_transforms(prob; scenario) intP = ComponentArrayInterpreter(θP) diff --git a/test/runtests.jl b/test/runtests.jl index 45f7450..2c87aa2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,8 @@ const GROUP = get(ENV, "GROUP", "All") # defined in in CI.yml @time begin if GROUP == "All" || GROUP == "Basic" + #@safetestset "test" include("test/test_bijectors_utils.jl") + @time @safetestset "test_bijectors_utils" include("test_bijectors_utils.jl") #@safetestset "test" include("test/test_ComponentArrayInterpreter.jl") @time @safetestset "test_ComponentArrayInterpreter" include("test_ComponentArrayInterpreter.jl") #@safetestset "test" include("test/test_ModelApplicator.jl") diff --git a/test/test_HybridProblem.jl b/test/test_HybridProblem.jl index 36939ea..1fc3916 100644 --- a/test/test_HybridProblem.jl +++ b/test/test_HybridProblem.jl @@ -1,6 +1,6 @@ using Test using HybridVariationalInference -using HybridVariationalInference: HybridVariationalInference as HVI +using HybridVariationalInference: HybridVariationalInference as CP using StableRNGs using Random using Statistics @@ -46,6 +46,8 @@ construct_problem = (;scenario=(:default,)) -> begin end n_out = length(θM) rng = StableRNG(111) + # n_batch = 10 + n_site, n_batch = get_hybridproblem_n_site_and_batch(CP.DoubleMM.DoubleMMCase(); scenario) # dependency on DeoubleMMCase -> take care of changes in covariates (; xM, n_site, θP_true, θMs_true, xP, y_global_true, y_true, y_global_o, y_o, y_unc ) = gen_hybridproblem_synthetic(rng, DoubleMM.DoubleMMCase()) @@ -62,13 +64,13 @@ construct_problem = (;scenario=(:default,)) -> begin # g, ϕg = construct_SimpleChainsApplicator(g_chain) # py = neg_logden_indep_normal - n_batch = 10 i_sites = 1:n_site - get_train_loader = let xM = xM, xP = xP, y_o = y_o, y_unc = y_unc, i_sites = i_sites - function inner_get_train_loader(; n_batch, kwargs...) - MLUtils.DataLoader((xM, xP, y_o, y_unc, i_sites), batchsize=n_batch, partial=false) - end - end + # get_train_loader = let xM = xM, xP = xP, y_o = y_o, y_unc = y_unc, i_sites = i_sites + # function inner_get_train_loader(; n_batch, kwargs...) + # MLUtils.DataLoader((xM, xP, y_o, y_unc, i_sites), batchsize=n_batch, partial=false) + # end + # end + train_dataloader = MLUtils.DataLoader((xM, xP, y_o, y_unc, i_sites), batchsize=n_batch, partial=false) θall = vcat(θP, θM) priors_dict = Dict{Symbol, Distribution}(keys(θall) .=> fit.(LogNormal, θall, QuantilePoint.(θall .* 3, 0.95))) priors_dict[:r1] = fit(Normal, θall.r1, qp_uu(3 * θall.r1)) # not transformed to log-scale @@ -79,8 +81,10 @@ construct_problem = (;scenario=(:default,)) -> begin #g_chain_scaled = app ϕunc0 = init_hybrid_ϕunc(cor_ends, zero(FT)) pbm_covars = (:covarK2 ∈ scenario) ? (:K2,) : () - HybridProblem(θP, θM, g_chain_scaled, ϕg0, ϕunc0, f_doubleMM_with_global, priors_dict, py, - transM, transP, get_train_loader, n_covar, n_site, cor_ends, pbm_covars) + HybridProblem(θP, θM, g_chain_scaled, ϕg0, ϕunc0, + f_doubleMM_with_global, f_doubleMM_with_global, priors_dict, py, + transM, transP, train_dataloader, n_covar, n_site, n_batch, + cor_ends, pbm_covars) end test_without_flux = (scenario) -> begin @@ -104,9 +108,10 @@ test_without_flux = (scenario) -> begin #----------- fit g and θP to y_o rng = StableRNG(111) g, ϕg0 = get_hybridproblem_MLapplicator(prob; scenario) - train_loader = get_hybridproblem_train_dataloader(prob; n_batch=10, scenario) + n_site, n_batch = get_hybridproblem_n_site_and_batch(prob; scenario) + train_loader = get_hybridproblem_train_dataloader(prob; scenario) (xM, xP, y_o, y_unc, i_sites) = first(train_loader) - f = get_hybridproblem_PBmodel(prob; scenario) + f = get_hybridproblem_PBmodel(prob; scenario, use_all_sites = false) par_templates = get_hybridproblem_par_templates(prob; scenario) #f(par_templates.θP, hcat(par_templates.θM, par_templates.θM), xP[1:2]) (; transM, transP) = get_hybridproblem_transforms(prob; scenario) @@ -120,7 +125,7 @@ test_without_flux = (scenario) -> begin y_global_o = Float64[] loss_gf = get_loss_gf(g, transM, transP, f, y_global_o, intϕ; pbm_covars) l1 = loss_gf(p0, first(train_loader)...) - tld = train_loader.data + tld = first(train_loader) gr = Zygote.gradient(p -> loss_gf(p, tld...)[1], CA.getdata(p0)) @test gr[1] isa Vector @@ -147,15 +152,15 @@ using GPUArraysCore import Flux gdev = gpu_device() -#methods(HVI.vec2uutri) +#methods(CP.vec2uutri) test_with_flux = (scenario) -> begin prob = probc = construct_problem(;scenario); @testset "HybridPointSolver" begin rng = StableRNG(111) - solver = HybridPointSolver(; alg=Adam(0.02), n_batch=11) - (; ϕ, resopt) = solve(prob, solver; scenario, rng, + solver = HybridPointSolver(; alg=Adam(0.02)) + (; ϕ, resopt, probo) = solve(prob, solver; scenario, rng, #callback = callback_loss(100), maxiters = 1200 #maxiters = 1200 #maxiters = 20 @@ -164,13 +169,17 @@ test_with_flux = (scenario) -> begin #gpu_handler = NullGPUDataHandler ) (; θP) = get_hybridproblem_par_templates(prob; scenario) - @test ϕ.ϕP.r0 < 1.5 * θP.r0 + θPo = (() -> begin + (; θP) = get_hybridproblem_par_templates(probo; scenario); + θP + end)() + @test θPo.r0 < 1.5 * θP.r0 @test ϕ.ϕP.K2 < 1.5 * log(θP.K2) end; @testset "HybridPosteriorSolver" begin rng = StableRNG(111) - solver = HybridPosteriorSolver(; alg=Adam(0.02), n_batch=11, n_MC=3) + solver = HybridPosteriorSolver(; alg=Adam(0.02), n_MC=3) (; ϕ, θP, resopt) = solve(prob, solver; scenario, rng, #callback = callback_loss(100), maxiters = 1200, #maxiters = 20 # too small so that it yields error @@ -189,9 +198,11 @@ test_with_flux = (scenario) -> begin @testset "HybridPosteriorSolver gpu" begin scenf = (scenario..., :use_Flux, :use_gpu, :omit_r0) rng = StableRNG(111) - prob = probg = HybridProblem(DoubleMM.DoubleMMCase(); scenario = scenf) - solver = HybridPosteriorSolver(; alg=Adam(0.02), n_batch=11, n_MC=3) - n_batches_in_epoch = get_hybridproblem_n_site(prob; scenario) ÷ solver.n_batch + # here using DoubleMMCase() directly rather than construct_problem + prob = probg = HybridProblem(DoubleMM.DoubleMMCase(); scenario = scenf); + solver = HybridPosteriorSolver(; alg=Adam(0.02), n_MC=3) + n_site, n_batch = get_hybridproblem_n_site_and_batch(prob; scenario = scenf) + n_batches_in_epoch = n_site ÷ n_batch (; ϕ, θP, resopt) = solve(prob, solver; scenario = scenf, rng, maxiters = 37, # smallest value by trial and error #maxiters = 20 # too small so that it yields error @@ -200,7 +211,11 @@ test_with_flux = (scenario) -> begin @test CA.getdata(ϕ) isa GPUArraysCore.AbstractGPUVector #@test cdev(ϕ.unc.ρsM)[1] > 0 # too few iterations in test -> may fail # - solver = HybridPosteriorSolver(; alg=Adam(0.02), n_batch=11, n_MC=3) + solver = HybridPosteriorSolver(; alg=Adam(0.02), n_MC=3) + (; ϕ, θP, resopt, probo) = solve(prob, solver; scenario = scenf, + maxiters = 37, + ); + @test cdev(ϕ.unc.ρsM)[1] > 0 test_correlation = () -> begin n_epoch = 100 # requires (; ϕ, θP, resopt, probo) = solve(prob, solver; scenario = scenf, @@ -221,24 +236,24 @@ test_with_flux = (scenario) -> begin end end; - # does not work with general Bijector: - # @testset "HybridPosteriorSolver also f on gpu" begin - # scenario = (:use_Flux, :use_gpu, :omit_r0, :f_on_gpu) - # rng = StableRNG(111) - # probg = HybridProblem(DoubleMM.DoubleMMCase(); scenario) - # prob = HVI.update(probg) - # #prob = HVI.update(probg, transM = identity, transP = identity) - # solver = HybridPosteriorSolver(; alg=Adam(0.02), n_batch=11, n_MC=3) - # n_batches_in_epoch = get_hybridproblem_n_site(prob; scenario) ÷ solver.n_batch - # (; ϕ, θP, resopt) = solve(prob, solver; scenario, rng, - # maxiters = 37, # smallest value by trial and error - # #maxiters = 20 # too small so that it yields error - # #θmean_quant = 0.01, # TODO make possible on gpu - # cdev = identity # do not move ζ to cpu # TODO infer in solve from scenario - # ) - # @test CA.getdata(ϕ) isa GPUArraysCore.AbstractGPUVector - # @test cdev(ϕ.unc.ρsM)[1] > 0 - # end; + @testset "HybridPosteriorSolver also f on gpu" begin + scenf = (scenario..., :use_Flux, :use_gpu, :omit_r0, :f_on_gpu) + rng = StableRNG(111) + probg = HybridProblem(DoubleMM.DoubleMMCase(); scenario = scenf); + # TODO replace exp by Exp Transformer + prob = CP.update(probg, transM = identity, transP = identity); + solver = HybridPosteriorSolver(; alg=Adam(0.02), n_MC=3) + n_site, n_batch = get_hybridproblem_n_site_and_batch(prob; scenario = scenf) + n_batches_in_epoch = n_site ÷ n_batch + (; ϕ, θP, resopt, probo) = solve(prob, solver; scenario = scenf, rng, + maxiters = 37, # smallest value by trial and error + #maxiters = 20 # too small so that it yields error + #θmean_quant = 0.01, # TODO make possible on gpu + cdev = identity # do not move ζ to cpu # TODO infer in solve from scenario + ); + @test CA.getdata(ϕ) isa GPUArraysCore.AbstractGPUVector + # @test cdev(ϕ.unc.ρsM)[1] > 0 # too few iterations + end; end # if gdev isa MLDataDevices.AbstractGPUDevice end # test_with flux diff --git a/test/test_doubleMM.jl b/test/test_doubleMM.jl index 2f59843..6a1b97b 100644 --- a/test/test_doubleMM.jl +++ b/test/test_doubleMM.jl @@ -21,8 +21,7 @@ using GPUArraysCore gdev = gpu_device() cdev = cpu_device() - -const prob = DoubleMM.DoubleMMCase() +prob = DoubleMM.DoubleMMCase() scenario = (:default,) #using Flux #scenario = (:use_Flux,) @@ -36,18 +35,19 @@ par_templates = get_hybridproblem_par_templates(prob; scenario) @test quantile(priors[:K2], 0.95) ≈ θall.K2 * 3 # fitted in f_doubleMM end -rng = StableRNG(111) +rng = StableRNG(111) # make sure to be the same as when constructing train_dataloader (; xM, n_site, θP_true, θMs_true, xP, y_global_true, y_true, y_global_o, y_o, y_unc ) = gen_hybridproblem_synthetic(rng, prob; scenario); i_sites = 1:n_site fneglogden = get_hybridproblem_neg_logden_obs(prob; scenario) - @testset "gen_hybridproblem_synthetic" begin @test isapprox( vec(mean(CA.getdata(θMs_true); dims = 2)), CA.getdata(par_templates.θM), rtol = 0.02) @test isapprox(vec(std(CA.getdata(θMs_true); dims = 2)), CA.getdata(par_templates.θM) .* 0.1, rtol = 0.02) + @test size(xP) == (n_site,) + @test size(y_o) == (8, n_site) # test same results for same rng rng2 = StableRNG(111) @@ -56,31 +56,32 @@ fneglogden = get_hybridproblem_neg_logden_obs(prob; scenario) end @testset "f_doubleMM_Matrix" begin - is = repeat(axes(θP_true,1)', n_site) + is = repeat(axes(θP_true, 1)', n_site) θvec = CA.ComponentVector(P = θP_true, Ms = θMs_true) - xPM = map(xP1s -> repeat(xP1s', n_site), xP[1]) + xPM = map(xP1s -> repeat(xP1s', n_site), xP[1]) #θ = hcat(θP_true[is], θMs_true') - intθ1 = get_concrete(ComponentArrayInterpreter(vcat(θP_true, θMs_true[:,1]))) + intθ1 = get_concrete(ComponentArrayInterpreter(vcat(θP_true, θMs_true[:, 1]))) #θpos = get_positions(intθ1) intθ = get_concrete(ComponentArrayInterpreter((n_site,), intθ1)) - fcost = (θvec, xPM, y_o, y_unc) -> begin + fy = (θvec, xPM) -> begin θ = hcat(CA.getdata(θvec.P[is]), CA.getdata(θvec.Ms')) y = HVI.DoubleMM.f_doubleMM(θ, xPM, intθ) #y = HVI.DoubleMM.f_doubleMM(θ, xPM, θpos) - fneglogden(y_o, y', y_unc) end - cost = fcost(θvec, xPM, y_o, y_unc) - ygrad = Zygote.gradient(θv -> fcost(θv, xPM), θvec)[1]; + y = fy(θvec, xPM) + y_exp = applyf(HVI.DoubleMM.f_doubleMM, θMs_true, θP_true, + Vector{eltype(θP_true)}(undef, 0), xP, intθ1) + @test y == y_exp' + ygrad = Zygote.gradient(θv -> sum(fy(θv, xPM)), θvec)[1] if gdev isa MLDataDevices.AbstractGPUDevice # θg = gdev(θ) # xPMg = gdev(xPM) # yg = HVI.DoubleMM.f_doubleMM(θg, xPMg, intθ); - θvecg = gdev(θvec); + θvecg = gdev(θvec) xPMg = gdev(xPM) - y_og = gdev(y_o) - y_uncg = gdev(y_unc) - costg = fcost(θvecg, xPMg, y_og, y_uncg); - ygradg = Zygote.gradient(θv -> fcost(θv, xPMg, y_og, y_uncg), θvecg)[1]; # erros without ";" + yg = fy(θvecg, xPMg) + @test cdev(yg) == y_exp' + ygradg = Zygote.gradient(θv -> sum(fy(θv, xPMg)), θvecg)[1] # erros without ";" @test ygradg isa CA.ComponentArray @test CA.getdata(ygradg) isa GPUArraysCore.AbstractGPUArray ygradgc = HVI.apply_preserve_axes(cdev, ygradg) # can print the cpu version @@ -90,32 +91,32 @@ end end @testset "neg_logden_obs Matrix" begin - is = repeat(axes(θP_true,1)', n_site) + is = repeat(axes(θP_true, 1)', n_site) θvec = CA.ComponentVector(P = θP_true, Ms = θMs_true) - xPM = map(xP1s -> repeat(xP1s', n_site), xP[1]) + xPM = map(xP1s -> repeat(xP1s', n_site), xP[1]) #θ = hcat(θP_true[is], θMs_true') - intθ1 = get_concrete(ComponentArrayInterpreter(vcat(θP_true, θMs_true[:,1]))) + intθ1 = get_concrete(ComponentArrayInterpreter(vcat(θP_true, θMs_true[:, 1]))) #θpos = get_positions(intθ1) intθ = get_concrete(ComponentArrayInterpreter((n_site,), intθ1)) - fl = (θvec, xPM) -> begin + fcost = (θvec, xPM, y_o, y_unc) -> begin θ = hcat(CA.getdata(θvec.P[is]), CA.getdata(θvec.Ms')) y = HVI.DoubleMM.f_doubleMM(θ, xPM, intθ) #y = HVI.DoubleMM.f_doubleMM(θ, xPM, θpos) - + fneglogden(y_o, y', y_unc) end - y = fy(θvec, xPM) - y_exp = applyf(HVI.DoubleMM.f_doubleMM, θMs_true, θP_true, Vector{eltype(θP_true)}(undef, 0), xP, intθ1) - @test y == y_exp' - ygrad = Zygote.gradient(θv -> sum(fy(θv, xPM)), θvec)[1]; + cost = fcost(θvec, xPM, y_o, y_unc) + ygrad = Zygote.gradient(θv -> fcost(θv, xPM, y_o, y_unc), θvec)[1] if gdev isa MLDataDevices.AbstractGPUDevice # θg = gdev(θ) # xPMg = gdev(xPM) # yg = HVI.DoubleMM.f_doubleMM(θg, xPMg, intθ); - θvecg = gdev(θvec); + θvecg = gdev(θvec) xPMg = gdev(xPM) - yg = fy(θvecg, xPMg); - @test cdev(yg) == y_exp' - ygradg = Zygote.gradient(θv -> sum(fy(θv, xPMg)), θvecg)[1]; # erros without ";" + y_og = gdev(y_o) + y_uncg = gdev(y_unc) + costg = fcost(θvecg, xPMg, y_og, y_uncg) + @test costg ≈ cost + ygradg = Zygote.gradient(θv -> fcost(θv, xPMg, y_og, y_uncg), θvecg)[1] # erros without ";" @test ygradg isa CA.ComponentArray @test CA.getdata(ygradg) isa GPUArraysCore.AbstractGPUArray ygradgc = HVI.apply_preserve_axes(cdev, ygradg) # can print the cpu version @@ -124,7 +125,6 @@ end end end - @testset "loss_g" begin g, ϕg0 = get_hybridproblem_MLapplicator(rng, prob; scenario) (; transP, transM) = get_hybridproblem_transforms(prob; scenario) @@ -169,23 +169,28 @@ end #----------- fit g and θP to y_o (without uncertainty, without transforming θP) g, ϕg0 = get_hybridproblem_MLapplicator(prob; scenario) (; transP, transM) = get_hybridproblem_transforms(prob; scenario) - f = get_hybridproblem_PBmodel(prob; scenario) + n_site, n_site_batch = get_hybridproblem_n_site_and_batch(prob; scenario) + f = get_hybridproblem_PBmodel(prob; scenario, use_all_sites = false) + f2 = get_hybridproblem_PBmodel(prob; scenario, use_all_sites = true) intϕ = ComponentArrayInterpreter(CA.ComponentVector( ϕg = 1:length(ϕg0), ϕP = par_templates.θP)) - p = p0 = vcat(ϕg0, HVI.apply_preserve_axes(inverse(transP), par_templates.θP) .- + p = p0 = vcat(ϕg0, + HVI.apply_preserve_axes(inverse(transP), par_templates.θP) .- convert(eltype(ϕg0), 0.1)) # slightly disturb θP_true #p = p0 = vcat(ϕg_opt1, par_templates.θP); # almost true # Pass the site-data for the batches as separate vectors wrapped in a tuple - n_batch = 10 - train_loader = MLUtils.DataLoader((xM, xP, y_o, y_unc, i_sites), batchsize = n_batch) - # get_hybridproblem_train_dataloader recreates synthetic data different θ_true - train_loader2 = get_hybridproblem_train_dataloader(prob; scenario, n_batch = n_site) - pbm_covars = get_hybridproblem_pbmpar_covars(prob; scenario) + # train_loader = MLUtils.DataLoader( + # (xM, xP, y_o, y_unc, i_sites), batchsize = n_site_batch) + train_loader = get_hybridproblem_train_dataloader(prob; scenario) + @assert train_loader.data == (xM, xP, y_o, y_unc, i_sites) + train_lodaer2 = MLUtils.DataLoader(train_loader.data, batchsize = n_site) + pbm_covars = get_hybridproblem_pbmpar_covars(prob; scenario) #loss_gf = get_loss_gf(g, transM, f, y_global_o, intϕ; gdev = identity) loss_gf = get_loss_gf(g, transM, transP, f, y_global_o, intϕ; pbm_covars) + loss_gf2 = get_loss_gf(g, transM, transP, f2, y_global_o, intϕ; pbm_covars) l1 = loss_gf(p0, first(train_loader)...)[1] (xM_batch, xP_batch, y_o_batch, y_unc_batch, i_sites_batch) = first(train_loader) Zygote.gradient( @@ -200,7 +205,7 @@ end #optprob, Adam(0.02), callback = callback_loss(100), maxiters = 5000); optprob, Adam(0.02), maxiters = 1000) - l1, y_pred_global, y_pred, θMs_pred, θP_pred = loss_gf(res.u, train_loader.data...) + l1, y_pred_global, y_pred, θMs_pred, θP_pred = loss_gf2(res.u, train_loader.data...) #l1, y_pred_global, y_pred, θMs_pred = loss_gf(p0, xM, xP, y_o, y_unc); θMs_pred = CA.ComponentArray(θMs_pred, CA.getaxes(θMs_true)) #TODO @test isapprox(par_templates.θP, intϕ(res.u).ϕP, rtol = 0.15) @@ -221,17 +226,17 @@ end end end - if gdev isa MLDataDevices.AbstractGPUDevice + scenario = (:use_Flux,) + g, ϕg0 = get_hybridproblem_MLapplicator(rng, prob; scenario) + ϕg0_gpu = gdev(ϕg0) + xM_gpu = gdev(xM) + g_gpu = gdev(g) + @testset "transfer NormalScalingModelApplicator to gpu" begin - scenario = (:use_Flux,) - g, ϕg0 = get_hybridproblem_MLapplicator(rng, prob; scenario) - ϕg = gdev(ϕg0) - xM_gpu = gdev(xM) - g_gpu = gdev(g) @test g_gpu.μ isa GPUArraysCore.AbstractGPUArray - y_gpu = g_gpu(xM_gpu, ϕg) + y_gpu = g_gpu(xM_gpu, ϕg0_gpu) y = g(xM, ϕg0) @test cdev(y_gpu) ≈ y - end; -end + end +end \ No newline at end of file diff --git a/test/test_elbo.jl b/test/test_elbo.jl index af85a73..0ee9aad 100644 --- a/test/test_elbo.jl +++ b/test/test_elbo.jl @@ -36,15 +36,17 @@ test_scenario = (scenario) -> begin #θsite_true = get_hybridproblem_par_templates(prob; scenario) + n_covar = 5 + n_site, n_batch = get_hybridproblem_n_site_and_batch(prob; scenario) + (; xM, n_site, θP_true, θMs_true, xP, y_global_true, y_true, y_global_o, y_o, y_unc + ) = gen_hybridproblem_synthetic(rng, prob; scenario); + g, ϕg0 = get_hybridproblem_MLapplicator(prob; scenario); - f = get_hybridproblem_PBmodel(prob; scenario) + f = get_hybridproblem_PBmodel(prob; scenario, use_all_sites = false) + f_pred = get_hybridproblem_PBmodel(prob; scenario, use_all_sites = true) - n_covar = 5 - n_batch = 10 n_θM, n_θP = values(map(length, get_hybridproblem_par_templates(prob; scenario))) - (; xM, n_site, θP_true, θMs_true, xP, y_global_true, y_true, y_global_o, y_o, y_unc - ) = gen_hybridproblem_synthetic(rng, prob; scenario); py = neg_logden_indep_normal @@ -149,7 +151,7 @@ test_scenario = (scenario) -> begin trans_PMs_gen = get_transPMs(n_site) @test length(intm_PMs_gen) == 402 @test trans_PMs_gen.length_in == 402 - (; θ, y) = predict_gf(rng, g, f, ϕ_ini, xM, xP, map(get_concrete, interpreters); + (; θ, y) = predict_gf(rng, g, f_pred, ϕ_ini, xM, xP, map(get_concrete, interpreters); get_transPMs, get_ca_int_PMs, n_sample_pred, cor_ends, pbm_covar_indices) @test θ isa CA.ComponentMatrix @test θ[:, 1].P.r0 > 0 @@ -162,7 +164,7 @@ test_scenario = (scenario) -> begin n_sample_pred = 200 ϕ = ggdev(CA.getdata(ϕ_ini)) xMg = ggdev(xM) - (; θ, y) = predict_gf(rng, g_gpu, f, ϕ, xMg, xP, map(get_concrete, interpreters); + (; θ, y) = predict_gf(rng, g_gpu, f_pred, ϕ, xMg, xP, map(get_concrete, interpreters); get_transPMs, get_ca_int_PMs, n_sample_pred, cor_ends, pbm_covar_indices) @test θ isa CA.ComponentMatrix # only ML parameters are on gpu @test θ[:, 1].P.r0 > 0 @@ -182,7 +184,7 @@ test_scenario = (scenario) -> begin # n_sample_pred = 200 # ϕ = ggdev(CA.getdata(ϕ_ini)) # xMg = ggdev(xM) - # (; θ, y) = predict_gf(rng, g_gpu, f, ϕ, xMg, ggdev(xP), map(get_concrete, interpreters); + # (; θ, y) = predict_gf(rng, g_gpu, f_pred, ϕ, xMg, ggdev(xP), map(get_concrete, interpreters); # get_ca_int_PMs, n_sample_pred, cor_ends, pbm_covar_indices, # get_transPMs = get_transPMs_ident, # cdev = identity); # keep on gpu diff --git a/test/test_sample_zeta.jl b/test/test_sample_zeta.jl index 16719d2..16413c9 100644 --- a/test/test_sample_zeta.jl +++ b/test/test_sample_zeta.jl @@ -65,7 +65,7 @@ end if ggdev isa MLDataDevices.AbstractGPUDevice @testset "sample_ζ_norm0 gpu" begin # sample only n_batch of 50 - n_batch = 40 + n_site, n_batch = get_hybridproblem_n_site_and_batch(prob; scenario) ϕb = CA.ComponentVector(P = ϕ_cpu.P, Ms = ϕ_cpu.Ms[:,1:n_batch], unc = ϕ_cpu.unc) intb = ComponentArrayInterpreter(ϕb) ϕ = ggdev(CA.getdata(ϕb)) From 7cd7c9ee3311e7a1fef12d6a668e97f6009daba6 Mon Sep 17 00:00:00 2001 From: Thomas Wutzler Date: Wed, 30 Apr 2025 12:07:28 +0200 Subject: [PATCH 4/7] replace transformations by Exp because exlementwise(exp) failed on AD on GPU test on GPU --- src/DoubleMM/f_doubleMM.jl | 14 ++++++++------ src/HybridVariationalInference.jl | 2 +- src/bijectors_exp.jl | 17 ----------------- src/bijectors_utils.jl | 21 +++++++++++++++++++++ test/test_HybridProblem.jl | 8 +++++--- test/test_bijectors_utils.jl | 5 +++++ 6 files changed, 40 insertions(+), 27 deletions(-) delete mode 100644 src/bijectors_exp.jl create mode 100644 src/bijectors_utils.jl diff --git a/src/DoubleMM/f_doubleMM.jl b/src/DoubleMM/f_doubleMM.jl index 564a56a..1436495 100644 --- a/src/DoubleMM/f_doubleMM.jl +++ b/src/DoubleMM/f_doubleMM.jl @@ -6,10 +6,10 @@ const θall = vcat(θP, θM) const θP_nor0 = θP[(:K2,)] -const transP = elementwise(exp) -const transM = elementwise(exp) +# const transP = elementwise(exp) +# const transM = elementwise(exp) -const transMS = Stacked(elementwise(identity), elementwise(exp)) +# const transMS = Stacked(elementwise(identity), elementwise(exp)) const int_θdoubleMM = ComponentArrayInterpreter(flatten1(CA.ComponentVector(; θP, θM))) @@ -114,15 +114,17 @@ end function HVI.get_hybridproblem_transforms(prob::DoubleMMCase; scenario::NTuple = ()) + _θP, _θM = get_hybridproblem_par_templates(prob; scenario) if (:stackedMS ∈ scenario) - return ((; transP, transM = transMS)) + return (; transP = Stacked((HVI.Exp(),),(1:length(_θP),)), + transM = Stacked((identity,HVI.Exp(),),(1:1, 2:length(_θM),))) elseif (:transIdent ∈ scenario) # identity transformations, should AD on GPU - _θP, _θM = get_hybridproblem_par_templates(prob; scenario) return (; transP = Stacked((identity,),(1:length(_θP),)), transM = Stacked((identity,),(1:length(_θM),))) end - (; transP, transM) + (; transP = Stacked((HVI.Exp(),),(1:length(_θP),)), + transM = Stacked((HVI.Exp(),),(1:length(_θM),))) end # function HVI.get_hybridproblem_sizes(::DoubleMMCase; scenario = ()) diff --git a/src/HybridVariationalInference.jl b/src/HybridVariationalInference.jl index a7fa4e9..095e19d 100644 --- a/src/HybridVariationalInference.jl +++ b/src/HybridVariationalInference.jl @@ -21,7 +21,7 @@ using StaticArrays: StaticArrays as SA using Functors #export Exp -include("bijectors_exp.jl") +include("bijectors_utils.jl") export ComponentArrayInterpreter, flatten1, get_concrete, get_positions include("ComponentArrayInterpreter.jl") diff --git a/src/bijectors_exp.jl b/src/bijectors_exp.jl deleted file mode 100644 index 6213440..0000000 --- a/src/bijectors_exp.jl +++ /dev/null @@ -1,17 +0,0 @@ -struct Exp <: Bijector -end - -#Functors.@functor Exp - -transform(b::Exp, x) = exp.(x) # note the broadcast -transform(ib::Inverse{<:Exp}, y) = log.(y) - -# `logabsdetjac` -logabsdetjac(b::Exp, x) = sum(x) - -`with_logabsdet_jacobian` -function Bijectors.with_logabsdet_jacobian(b::Exp, x) - return transform(b, x), logabsdetjac(b, x) -end - -is_monotonically_increasing(::Exp) = true diff --git a/src/bijectors_utils.jl b/src/bijectors_utils.jl new file mode 100644 index 0000000..3cc70ee --- /dev/null +++ b/src/bijectors_utils.jl @@ -0,0 +1,21 @@ +struct Exp <: Bijector +end + +#Functors.@functor Exp +Bijectors.transform(b::Exp, x) = exp.(x) # note the broadcast +Bijectors.transform(ib::Inverse{<:Exp}, y) = log.(y) + +# `logabsdetjac` +Bijectors.logabsdetjac(b::Exp, x) = sum(x) + +`with_logabsdet_jacobian` +function Bijectors.with_logabsdet_jacobian(b::Exp, x) + return exp.(x), sum(x) +end +# function Bijectors.with_logabsdet_jacobian(ib::Inverse{<:Exp}, y) +# x = transform(ib, y) +# return x, -logabsdetjac(inverse(ib), x) +# end + + +Bijectors.is_monotonically_increasing(::Exp) = true diff --git a/test/test_HybridProblem.jl b/test/test_HybridProblem.jl index 1fc3916..4f004bd 100644 --- a/test/test_HybridProblem.jl +++ b/test/test_HybridProblem.jl @@ -199,6 +199,7 @@ test_with_flux = (scenario) -> begin scenf = (scenario..., :use_Flux, :use_gpu, :omit_r0) rng = StableRNG(111) # here using DoubleMMCase() directly rather than construct_problem + #(;transP, transM) = get_hybridproblem_transforms(DoubleMM.DoubleMMCase(); scenario = scenf) prob = probg = HybridProblem(DoubleMM.DoubleMMCase(); scenario = scenf); solver = HybridPosteriorSolver(; alg=Adam(0.02), n_MC=3) n_site, n_batch = get_hybridproblem_n_site_and_batch(prob; scenario = scenf) @@ -216,13 +217,15 @@ test_with_flux = (scenario) -> begin maxiters = 37, ); @test cdev(ϕ.unc.ρsM)[1] > 0 + @test probo.ϕunc == cdev(ϕ.unc) test_correlation = () -> begin - n_epoch = 100 # requires + n_epoch = 20 # requires (; ϕ, θP, resopt, probo) = solve(prob, solver; scenario = scenf, maxiters = n_batches_in_epoch * n_epoch, callback = callback_loss(n_batches_in_epoch*5) ); @test cdev(ϕ.unc.ρsM)[1] > 0 + @test probo.ϕunc == cdev(ϕ.unc) # predict using problem and its associated dataloader (; θ, y, entropy_ζ) = predict_gf(rng, probo; scenario = scenf, n_sample_pred = 200); mean_θ = CA.ComponentVector(mean(CA.getdata(θ); dims = 2)[:, 1], CA.getaxes(θ[:, 1])[1]) @@ -240,8 +243,7 @@ test_with_flux = (scenario) -> begin scenf = (scenario..., :use_Flux, :use_gpu, :omit_r0, :f_on_gpu) rng = StableRNG(111) probg = HybridProblem(DoubleMM.DoubleMMCase(); scenario = scenf); - # TODO replace exp by Exp Transformer - prob = CP.update(probg, transM = identity, transP = identity); + #prob = CP.update(probg, transM = identity, transP = identity); solver = HybridPosteriorSolver(; alg=Adam(0.02), n_MC=3) n_site, n_batch = get_hybridproblem_n_site_and_batch(prob; scenario = scenf) n_batches_in_epoch = n_site ÷ n_batch diff --git a/test/test_bijectors_utils.jl b/test/test_bijectors_utils.jl index c666cac..351ee86 100644 --- a/test/test_bijectors_utils.jl +++ b/test/test_bijectors_utils.jl @@ -22,6 +22,7 @@ b2 = elementwise(exp) b2s = Stacked((b2,b2),(1:3,4:4)) b3 = HybridVariationalInference.Exp() b3s = Stacked((b3,b3), (1:3,4:4)) +#b3s = Stacked((b3,),(1:4,)) y = trans(x, b2) @@ -35,6 +36,10 @@ dy = Zygote.gradient(x -> trans(x,b2), x) end; @testset "Exp" begin + y1 = b3(x) + y2 = b3s(x) + @test all(inverse(b3)(y2) .≈ x) + @test all(inverse(b3s)(y2) .≈ x) ye = trans(x, b3) dye = Zygote.gradient(x -> trans(x,b3), x) @test ye == y From d75482184ec6701e4f94e51dc1da5a3fb6b59c04 Mon Sep 17 00:00:00 2001 From: Thomas Wutzler Date: Wed, 30 Apr 2025 12:23:53 +0200 Subject: [PATCH 5/7] remove n_site from gen_hybridproblem_synthetic --- dev/doubleMM.jl | 4 ++-- src/DoubleMM/f_doubleMM.jl | 1 - test/test_HybridProblem.jl | 2 +- test/test_doubleMM.jl | 3 ++- test/test_elbo.jl | 2 +- test/test_sample_zeta.jl | 2 +- 6 files changed, 7 insertions(+), 7 deletions(-) diff --git a/dev/doubleMM.jl b/dev/doubleMM.jl index a534183..d443ee8 100644 --- a/dev/doubleMM.jl +++ b/dev/doubleMM.jl @@ -29,8 +29,8 @@ cdev = gdev isa MLDataDevices.AbstractGPUDevice ? cpu_device() : identity #------ setup synthetic data and training data loader prob0_ = HybridProblem(DoubleMM.DoubleMMCase(); scenario); -(; xM, n_site, θP_true, θMs_true, xP, y_global_true, y_true, y_global_o, y_o, y_unc -) = gen_hybridproblem_synthetic(rng, prob0_; scenario); +(; xM, θP_true, θMs_true, xP, y_global_true, y_true, y_global_o, y_o, y_unc +) = gen_hybridproblem_synthetic(rng, DoubleMM.DoubleMMCase(); scenario); n_site, n_batch = get_hybridproblem_n_site_and_batch(prob0_; scenario) ζP_true, ζMs_true = log.(θP_true), log.(θMs_true) i_sites = 1:n_site diff --git a/src/DoubleMM/f_doubleMM.jl b/src/DoubleMM/f_doubleMM.jl index 1436495..a5ef56a 100644 --- a/src/DoubleMM/f_doubleMM.jl +++ b/src/DoubleMM/f_doubleMM.jl @@ -253,7 +253,6 @@ function HVI.gen_hybridproblem_synthetic(rng::AbstractRNG, prob::DoubleMMCase; y_o = y_true .+ randn(rng, FloatType, size(y_true)) .* σ_o (; xM, - n_site, θP_true = θP, θMs_true, xP, diff --git a/test/test_HybridProblem.jl b/test/test_HybridProblem.jl index 4f004bd..e9d717f 100644 --- a/test/test_HybridProblem.jl +++ b/test/test_HybridProblem.jl @@ -49,7 +49,7 @@ construct_problem = (;scenario=(:default,)) -> begin # n_batch = 10 n_site, n_batch = get_hybridproblem_n_site_and_batch(CP.DoubleMM.DoubleMMCase(); scenario) # dependency on DeoubleMMCase -> take care of changes in covariates - (; xM, n_site, θP_true, θMs_true, xP, y_global_true, y_true, y_global_o, y_o, y_unc + (; xM, θP_true, θMs_true, xP, y_global_true, y_true, y_global_o, y_o, y_unc ) = gen_hybridproblem_synthetic(rng, DoubleMM.DoubleMMCase()) n_covar = size(xM,1) n_input = (:covarK2 ∈ scenario) ? n_covar +1 : n_covar diff --git a/test/test_doubleMM.jl b/test/test_doubleMM.jl index 6a1b97b..0fec9b0 100644 --- a/test/test_doubleMM.jl +++ b/test/test_doubleMM.jl @@ -36,8 +36,9 @@ par_templates = get_hybridproblem_par_templates(prob; scenario) end rng = StableRNG(111) # make sure to be the same as when constructing train_dataloader -(; xM, n_site, θP_true, θMs_true, xP, y_global_true, y_true, y_global_o, y_o, y_unc +(; xM, θP_true, θMs_true, xP, y_global_true, y_true, y_global_o, y_o, y_unc ) = gen_hybridproblem_synthetic(rng, prob; scenario); +n_site, n_batch = get_hybridproblem_n_site_and_batch(prob; scenario) i_sites = 1:n_site fneglogden = get_hybridproblem_neg_logden_obs(prob; scenario) diff --git a/test/test_elbo.jl b/test/test_elbo.jl index 0ee9aad..3c1cf12 100644 --- a/test/test_elbo.jl +++ b/test/test_elbo.jl @@ -38,7 +38,7 @@ test_scenario = (scenario) -> begin #θsite_true = get_hybridproblem_par_templates(prob; scenario) n_covar = 5 n_site, n_batch = get_hybridproblem_n_site_and_batch(prob; scenario) - (; xM, n_site, θP_true, θMs_true, xP, y_global_true, y_true, y_global_o, y_o, y_unc + (; xM, θP_true, θMs_true, xP, y_global_true, y_true, y_global_o, y_o, y_unc ) = gen_hybridproblem_synthetic(rng, prob; scenario); g, ϕg0 = get_hybridproblem_MLapplicator(prob; scenario); diff --git a/test/test_sample_zeta.jl b/test/test_sample_zeta.jl index 16413c9..eb0db4a 100644 --- a/test/test_sample_zeta.jl +++ b/test/test_sample_zeta.jl @@ -23,7 +23,7 @@ scenario = (:default,) n_θM, n_θP = length.(values(get_hybridproblem_par_templates(prob; scenario))) -(; xM, n_site, θP_true, θMs_true, xP, y_global_true, y_true, y_global_o, y_o +(; xM, θP_true, θMs_true, xP, y_global_true, y_true, y_global_o, y_o ) = gen_hybridproblem_synthetic(rng, prob; scenario) FT = get_hybridproblem_float_type(prob; scenario) From 705fb806dec267ff85890cefde03d55bb63ca733 Mon Sep 17 00:00:00 2001 From: Thomas Wutzler Date: Wed, 30 Apr 2025 13:17:01 +0200 Subject: [PATCH 6/7] typos --- src/ComponentArrayInterpreter.jl | 4 ++-- test/test_doubleMM.jl | 5 ++--- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/src/ComponentArrayInterpreter.jl b/src/ComponentArrayInterpreter.jl index 6f89b35..fe5bb2c 100644 --- a/src/ComponentArrayInterpreter.jl +++ b/src/ComponentArrayInterpreter.jl @@ -5,12 +5,12 @@ Interface for Type that implements - `as_ca(::AbstractArray, interpreter) -> ComponentArray` -- `CompoentArrays.getaxes(interpeter)` +- `ComponentArrays.getaxes(interpreter)` - `Base.length(interpreter) -> Int` When called on a vector, forwards to `as_ca`. -There is a default implementation for Base.length based on CompoentArrays.getaxes. +There is a default implementation for Base.length based on ComponentArrays.getaxes. """ abstract type AbstractComponentArrayInterpreter end diff --git a/test/test_doubleMM.jl b/test/test_doubleMM.jl index 0fec9b0..e3e84ad 100644 --- a/test/test_doubleMM.jl +++ b/test/test_doubleMM.jl @@ -82,7 +82,7 @@ end xPMg = gdev(xPM) yg = fy(θvecg, xPMg) @test cdev(yg) == y_exp' - ygradg = Zygote.gradient(θv -> sum(fy(θv, xPMg)), θvecg)[1] # erros without ";" + ygradg = Zygote.gradient(θv -> sum(fy(θv, xPMg)), θvecg)[1] # errors without ";" @test ygradg isa CA.ComponentArray @test CA.getdata(ygradg) isa GPUArraysCore.AbstractGPUArray ygradgc = HVI.apply_preserve_axes(cdev, ygradg) # can print the cpu version @@ -117,7 +117,7 @@ end y_uncg = gdev(y_unc) costg = fcost(θvecg, xPMg, y_og, y_uncg) @test costg ≈ cost - ygradg = Zygote.gradient(θv -> fcost(θv, xPMg, y_og, y_uncg), θvecg)[1] # erros without ";" + ygradg = Zygote.gradient(θv -> fcost(θv, xPMg, y_og, y_uncg), θvecg)[1] # errors without ";" @test ygradg isa CA.ComponentArray @test CA.getdata(ygradg) isa GPUArraysCore.AbstractGPUArray ygradgc = HVI.apply_preserve_axes(cdev, ygradg) # can print the cpu version @@ -186,7 +186,6 @@ end # (xM, xP, y_o, y_unc, i_sites), batchsize = n_site_batch) train_loader = get_hybridproblem_train_dataloader(prob; scenario) @assert train_loader.data == (xM, xP, y_o, y_unc, i_sites) - train_lodaer2 = MLUtils.DataLoader(train_loader.data, batchsize = n_site) pbm_covars = get_hybridproblem_pbmpar_covars(prob; scenario) #loss_gf = get_loss_gf(g, transM, f, y_global_o, intϕ; gdev = identity) From d4594e7880bdd04b04012d24ef07e70c8d0cee9a Mon Sep 17 00:00:00 2001 From: Thomas Wutzler Date: Wed, 30 Apr 2025 14:39:00 +0200 Subject: [PATCH 7/7] fix test where n_site is not provided any more --- test/test_sample_zeta.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/test_sample_zeta.jl b/test/test_sample_zeta.jl index eb0db4a..4749c11 100644 --- a/test/test_sample_zeta.jl +++ b/test/test_sample_zeta.jl @@ -25,6 +25,8 @@ n_θM, n_θP = length.(values(get_hybridproblem_par_templates(prob; scenario))) (; xM, θP_true, θMs_true, xP, y_global_true, y_true, y_global_o, y_o ) = gen_hybridproblem_synthetic(rng, prob; scenario) +n_site, n_batch = get_hybridproblem_n_site_and_batch(prob; scenario) + FT = get_hybridproblem_float_type(prob; scenario)