basic distribution and getters based on GRF.jl

Project.toml

format

add plot to tests

fix global_rng consistency

improved interface and seperated coeff vs function sampling

constrained sampling from distribution and testing

sample from prior test

changed transforms to be more consistent and compatible with func dists

removed dimensional ambiguities in ndims for functions and tested c->u and u->c

get_logpdf -> logpdf

2D examples, unit tests done

docstrings

codecov

docstrings

format

docstring?

docstrings

Darcy example

removed additional sampling in KL

adjusted tests for user dist, and new func construction

removed unnecessary functions

plot bugfix

fixed bug runtest

removed rng from build without coeffs, more robust tests

rm rng

codecov

Project.toml

@@ -7,6 +7,7 @@ version = "0.14.2"
 Convex = "f65535da-76fb-5f13-bab9-19810c17039a"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+GaussianRandomFields = "e4b2fa32-6e09-5554-b718-106ed5adafe9"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
@@ -33,6 +34,7 @@ julia = "1.5"
 [extras]
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 
 [targets]
-test = ["StableRNGs", "Test"]
+test = ["StableRNGs", "Test", "Plots"]

docs/src/API/ParameterDistributions.md

@@ -4,33 +4,60 @@
 CurrentModule = EnsembleKalmanProcesses.ParameterDistributions
 ```
 
+## ParameterDistributionTypes
 ```@docs
-ParameterDistribution
-constrained_gaussian
 Parameterized
 Samples
+VectorOfParameterized
+```
+
+## Constraints
+```@docs
 Constraint
 no_constraint
 bounded_below
 bounded_above
 bounded
 length(c::CType) where {CType <: ConstraintType}
 size(c::CType) where {CType <: ConstraintType}
+```
+
+## ParameterDistributions
+
+```@docs
+ParameterDistribution
+constrained_gaussian
 n_samples
 get_name
 get_dimensions
 get_n_samples
-get_all_constraints
-get_bounds
+get_all_constraints(::ParameterDistribution)
 get_constraint_type
+get_bounds
 batch
 get_distribution
 sample
-get_logpdf
+logpdf
 mean
 var
 cov
-transform_constrained_to_unconstrained
-transform_unconstrained_to_constrained
+transform_constrained_to_unconstrained(::ParameterDistribution, ::AbstractVector)
+transform_constrained_to_unconstrained(::ParameterDistribution, ::AbstractMatrix)
+transform_constrained_to_unconstrained(::ParameterDistribution, ::Dict)
+transform_unconstrained_to_constrained(::ParameterDistribution, ::AbstractVector)
+transform_unconstrained_to_constrained(::ParameterDistribution, ::AbstractMatrix)
+transform_unconstrained_to_constrained(::ParameterDistribution, ::Dict)
 ```
 
+## FunctionParameterDistributions
+
+```@docs
+GaussianRandomFieldsPackage
+GaussianRandomFieldInterface
+ndims(grfi::GaussianRandomFieldInterface)
+get_all_constraints(grfi::GaussianRandomFieldInterface)
+transform_constrained_to_unconstrained(::GaussianRandomFieldInterface, ::AbstractVector, ::AbstractVector{FT}) where {FT <: Real}
+transform_constrained_to_unconstrained(::GaussianRandomFieldInterface, ::AbstractVector, ::AbstractMatrix{FT}) where {FT <: Real}
+transform_unconstrained_to_constrained(::GaussianRandomFieldInterface, ::AbstractVector, ::AbstractVector{FT}) where {FT <: Real}
+transform_unconstrained_to_constrained(::GaussianRandomFieldInterface, ::AbstractVector, ::AbstractMatrix{FT}) where {FT <: Real}
+```

examples/Darcy/GModel.jl

@@ -0,0 +1,219 @@
+##################
+# Copied on 3/16/23 and modified from 
+# https://github.com/Zhengyu-Huang/InverseProblems.jl/blob/master/Fluid/Darcy-2D.jl
+##################
+
+using JLD2
+using Statistics
+using LinearAlgebra
+using Distributions
+using Random
+using SparseArrays
+
+
+
+
+mutable struct Setup_Param{FT <: AbstractFloat, IT <: Int}
+    # physics
+    N::IT         # number of grid points for both x and y directions (including both ends)
+    Δx::FT
+    xx::AbstractVector{FT}  # uniform grid [a, a+Δx, a+2Δx ... b] (in each dimension)
+
+    #for source term
+    f_2d::AbstractMatrix{FT}
+
+    κ::AbstractMatrix{FT}
+
+    # observation locations is tensor product x_locs × y_locs
+    x_locs::AbstractVector{IT}
+    y_locs::AbstractVector{IT}
+
+    N_y::IT
+end
+
+
+function Setup_Param(
+    xx::AbstractVector{FT},
+    obs_ΔN::IT,
+    κ::AbstractMatrix;
+    seed::IT = 123,
+) where {FT <: AbstractFloat, IT <: Int}
+
+    N = length(xx)
+    Δx = xx[2] - xx[1]
+
+    #  logκ_2d, φ, λ, θ_ref = generate_θ_KL(xx, N_KL, d, τ, seed=seed)
+    f_2d = compute_f_2d(xx)
+
+    x_locs = Array(obs_ΔN:obs_ΔN:(N - obs_ΔN))
+    y_locs = Array(obs_ΔN:obs_ΔN:(N - obs_ΔN))
+    N_y = length(x_locs) * length(y_locs)
+
+    Setup_Param(N, Δx, xx, f_2d, κ, x_locs, y_locs, N_y)
+end
+
+
+
+#=
+A hardcoding source function, 
+which assumes the computational domain is
+[0 1]×[0 1]
+f(x,y) = f(y),
+which dependes only on y
+=#
+function compute_f_2d(yy::AbstractVector{FT}) where {FT <: AbstractFloat}
+    N = length(yy)
+    f_2d = zeros(FT, N, N)
+    for i in 1:N
+        if (yy[i] <= 4 / 6)
+            f_2d[:, i] .= 1000.0
+        elseif (yy[i] >= 4 / 6 && yy[i] <= 5 / 6)
+            f_2d[:, i] .= 2000.0
+        elseif (yy[i] >= 5 / 6)
+            f_2d[:, i] .= 3000.0
+        end
+    end
+    return f_2d
+end
+
+"""
+    run_G_ensemble(darcy,κs::AbstractMatrix)
+
+Computes the forward map `G` (`solve_Darcy_2D` followed by `compute_obs`) over an ensemble of `κ`'s, stored flat as columns of `κs`
+"""
+function run_G_ensemble(darcy, κs::AbstractMatrix)
+    N_ens = size(κs, 2) # ens size
+    nd = darcy.N_y #num obs
+    g_ens = zeros(nd, N_ens)
+    for i in 1:N_ens
+        # run the model with the current parameters, i.e., map θ to G(θ)
+        κ_i = reshape(κs[:, i], darcy.N, darcy.N) # unflatten
+        h_i = solve_Darcy_2D(darcy, κ_i) # run model
+        g_ens[:, i] = compute_obs(darcy, h_i) # observe solution
+    end
+    return g_ens
+end
+
+
+
+#=
+    return the unknow index for the grid point
+    Since zero-Dirichlet boundary conditions are imposed on  
+    all four edges, the freedoms are only on interior points
+=#
+function ind(darcy::Setup_Param{FT, IT}, ix::IT, iy::IT) where {FT <: AbstractFloat, IT <: Int}
+    return (ix - 1) + (iy - 2) * (darcy.N - 2)
+end
+
+function ind_all(darcy::Setup_Param{FT, IT}, ix::IT, iy::IT) where {FT <: AbstractFloat, IT <: Int}
+    return ix + (iy - 1) * darcy.N
+end
+
+#=
+    solve Darcy equation with finite difference method:
+    -∇(κ∇h) = f
+    with Dirichlet boundary condition, h=0 on ∂Ω
+=#
+function solve_Darcy_2D(darcy::Setup_Param{FT, IT}, κ_2d::AbstractMatrix{FT}) where {FT <: AbstractFloat, IT <: Int}
+    Δx, N = darcy.Δx, darcy.N
+
+    indx = IT[]
+    indy = IT[]
+    vals = FT[]
+
+    f_2d = darcy.f_2d
+
+    𝓒 = Δx^2
+    for iy in 2:(N - 1)
+        for ix in 2:(N - 1)
+
+            ixy = ind(darcy, ix, iy)
+
+            #top
+            if iy == N - 1
+                #ft = -(κ_2d[ix, iy] + κ_2d[ix, iy+1])/2.0 * (0 - h_2d[ix,iy])
+                push!(indx, ixy)
+                push!(indy, ixy)
+                push!(vals, (κ_2d[ix, iy] + κ_2d[ix, iy + 1]) / 2.0 / 𝓒)
+            else
+                #ft = -(κ_2d[ix, iy] + κ_2d[ix, iy+1])/2.0 * (h_2d[ix,iy+1] - h_2d[ix,iy])
+                append!(indx, [ixy, ixy])
+                append!(indy, [ixy, ind(darcy, ix, iy + 1)])
+                append!(
+                    vals,
+                    [(κ_2d[ix, iy] + κ_2d[ix, iy + 1]) / 2.0 / 𝓒, -(κ_2d[ix, iy] + κ_2d[ix, iy + 1]) / 2.0 / 𝓒],
+                )
+            end
+
+            #bottom
+            if iy == 2
+                #fb = (κ_2d[ix, iy] + κ_2d[ix, iy-1])/2.0 * (h_2d[ix,iy] - 0)
+                push!(indx, ixy)
+                push!(indy, ixy)
+                push!(vals, (κ_2d[ix, iy] + κ_2d[ix, iy - 1]) / 2.0 / 𝓒)
+            else
+                #fb = (κ_2d[ix, iy] + κ_2d[ix, iy-1])/2.0 * (h_2d[ix,iy] - h_2d[ix,iy-1])
+                append!(indx, [ixy, ixy])
+                append!(indy, [ixy, ind(darcy, ix, iy - 1)])
+                append!(
+                    vals,
+                    [(κ_2d[ix, iy] + κ_2d[ix, iy - 1]) / 2.0 / 𝓒, -(κ_2d[ix, iy] + κ_2d[ix, iy - 1]) / 2.0 / 𝓒],
+                )
+            end
+
+            #right
+            if ix == N - 1
+                #fr = -(κ_2d[ix, iy] + κ_2d[ix+1, iy])/2.0 * (0 - h_2d[ix,iy])
+                push!(indx, ixy)
+                push!(indy, ixy)
+                push!(vals, (κ_2d[ix, iy] + κ_2d[ix + 1, iy]) / 2.0 / 𝓒)
+            else
+                #fr = -(κ_2d[ix, iy] + κ_2d[ix+1, iy])/2.0 * (h_2d[ix+1,iy] - h_2d[ix,iy])
+                append!(indx, [ixy, ixy])
+                append!(indy, [ixy, ind(darcy, ix + 1, iy)])
+                append!(
+                    vals,
+                    [(κ_2d[ix, iy] + κ_2d[ix + 1, iy]) / 2.0 / 𝓒, -(κ_2d[ix, iy] + κ_2d[ix + 1, iy]) / 2.0 / 𝓒],
+                )
+            end
+
+            #left
+            if ix == 2
+                #fl = (κ_2d[ix, iy] + κ_2d[ix-1, iy])/2.0 * (h_2d[ix,iy] - 0)
+                push!(indx, ixy)
+                push!(indy, ixy)
+                push!(vals, (κ_2d[ix, iy] + κ_2d[ix - 1, iy]) / 2.0 / 𝓒)
+            else
+                #fl = (κ_2d[ix, iy] + κ_2d[ix-1, iy])/2.0 * (h_2d[ix,iy] - h_2d[ix-1,iy])
+                append!(indx, [ixy, ixy])
+                append!(indy, [ixy, ind(darcy, ix - 1, iy)])
+                append!(
+                    vals,
+                    [(κ_2d[ix, iy] + κ_2d[ix - 1, iy]) / 2.0 / 𝓒, -(κ_2d[ix, iy] + κ_2d[ix - 1, iy]) / 2.0 / 𝓒],
+                )
+            end
+
+        end
+    end
+
+
+
+    df = sparse(indx, indy, vals, (N - 2)^2, (N - 2)^2)
+    # Multithread does not support sparse matrix solver
+    h = df \ (f_2d[2:(N - 1), 2:(N - 1)])[:]
+
+    h_2d = zeros(FT, N, N)
+    h_2d[2:(N - 1), 2:(N - 1)] .= reshape(h, N - 2, N - 2)
+
+    return h_2d
+end
+
+#=
+Compute observation values
+=#
+function compute_obs(darcy::Setup_Param{FT, IT}, h_2d::AbstractMatrix{FT}) where {FT <: AbstractFloat, IT <: Int}
+    # X---X(1)---X(2) ... X(obs_N)---X
+    obs_2d = h_2d[darcy.x_locs, darcy.y_locs]
+
+    return obs_2d[:]
+end

examples/Darcy/Project.toml

@@ -0,0 +1,10 @@
+[deps]
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+EnsembleKalmanProcesses = "aa8a2aa5-91d8-4396-bcef-d4f2ec43552d"
+GaussianRandomFields = "e4b2fa32-6e09-5554-b718-106ed5adafe9"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"