Skip to content
New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Jump to bottom

Laplace approximation to posterior #20

Merged
merged 29 commits into from
Apr 24, 2024
Merged

Laplace approximation to posterior #20

Show file tree
Hide file tree
Changes from 20 commits
Commits
Show all changes
29 commits
Select commit Hold shift + click to select a range
b5e53af
laplace approximation
Apr 10, 2024
e705e03
update __del__ method
Apr 10, 2024
0ac0944
minor changes
Apr 11, 2024
68a5c03
minor comment fixes
Apr 11, 2024
c65c361
changes to del and solve method
Apr 11, 2024
2f14024
formatting fix
Apr 11, 2024
db41420
prior and posterior samples
Apr 12, 2024
d06de44
change posterior mean, formatting
Apr 12, 2024
f23f48d
test for LowRankHessain added
Apr 12, 2024
e4d6e46
test for lowRankHessian
Apr 12, 2024
6a30e9d
posterior samples, fix low rank Hessian test
Apr 12, 2024
07b5348
changes to writing out prior, posterior samples, llaplaceApproximatio…
Apr 15, 2024
0a4ec14
removing commented code
Apr 15, 2024
847b2f7
removing code for Variational Regularization prior
Apr 15, 2024
b096580
linter fix
Apr 15, 2024
bcfe06f
fixing solve method for cgsolver
Apr 16, 2024
fb0cf34
adding tests for sampling and preconditioner
Apr 16, 2024
ee5edbb
misfit argument fix
Apr 16, 2024
085d5bf
return type changed
Apr 16, 2024
a8ecee4
sampler check not working
Apr 16, 2024
b44b0ae
checking sample test
Apr 17, 2024
ecf7c85
rename sample to mult
Apr 17, 2024
81d7a9a
time series function plot not working
Apr 17, 2024
3ff241a
changing noise level, sampler
Apr 19, 2024
c31e622
writing samples in xdmf format working
Apr 19, 2024
0cb7b58
Modified script so that it saves as time series.
uvilla Apr 19, 2024
6f96898
Rename unused variables
uvilla Apr 19, 2024
1026f11
prior and posterior sample time series plots
Apr 19, 2024
4ef37ef
spelling fix
Apr 24, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 15 additions & 0 deletions .github/workflows/CI_testing.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,21 @@ jobs:
cd ./hippylibX/test &&
mpirun -n 2 python3 test_eigendecomposition.py

- name: low rank Hessian testing
run: |
cd ./hippylibX/test &&
mpirun -n 2 python3 test_lowRankHessian.py

- name: low rank Hessian precondtioner testing
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Spelling: preconditioner

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Spelling changed. Thanks for pointing it out.

run: |
cd ./hippylibX/test &&
mpirun -n 2 python3 test_lowRankHessian_preconditioner.py

- name: sample testing
run: |
cd ./hippylibX/test &&
mpirun -n 2 python3 test_sampling.py

- name: run serial check
run: |
cd ./hippylibX/test &&
Expand Down
60 changes: 47 additions & 13 deletions example/poisson_dirichlet_example.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,6 @@
from matplotlib import pyplot as plt
from typing import Sequence, Dict


sys.path.append(os.environ.get("HIPPYLIBX_BASE_DIR", "../"))
import hippylibX as hpx

Expand Down Expand Up @@ -184,7 +183,7 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]:

with dlx.io.VTXWriter(
msh.comm,
"poisson_Dirichlet_BiLaplacian_prior_np{0:d}_Prior.bp".format(nproc),
"poisson_Dirichlet_BiLaplacian_prior_np{0:d}.bp".format(nproc),
[m_fun, m_true_fun, u_map_fun, u_true_fun, d_fun],
) as vtx:
vtx.write(0.0)
Expand Down Expand Up @@ -215,17 +214,51 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]:

hpx.parRandom.normal(1.0, Omega)

d, U = hpx.doublePassG(
Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1, check=False
)
d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1)

eigen_decomposition_results = {
"A": Hmisfit.mat,
"B": prior.R,
"k": k,
"d": d,
"U": U,
}
# generating prior and posterior samples
lap_aprx = hpx.LaplaceApproximator(prior, d, U)
lap_aprx.mean = prior.generate_parameter(0)
lap_aprx.mean.array[:] = x[hpx.PARAMETER].array[:]

m_prior = prior.generate_parameter(0)
m_post = prior.generate_parameter(0)

noise = prior.generate_parameter("noise")

num_samples_generate = 5

prior_samples = []
posterior_samples = []
for i in range(num_samples_generate):
hpx.parRandom.normal(1.0, noise)
lap_aprx.sample(noise, m_prior, m_post)
prior_sample = hpx.vector2Function(
m_prior, Vh[hpx.PARAMETER], name=f"prior_sample_{i}"
)
posterior_sample = hpx.vector2Function(
m_post, Vh[hpx.PARAMETER], name=f"posterior_sample_{i}"
)
prior_samples.append(prior_sample)
posterior_samples.append(posterior_sample)

with dlx.io.VTXWriter(
msh.comm,
"poisson_Dirichlet_prior_Bilaplacian_samples_prior_np{0:d}.bp".format(nproc),
prior_samples,
) as vtx:
vtx.write(0.0)

with dlx.io.VTXWriter(
msh.comm,
"poisson_Dirichlet_prior_Bilaplacian_samples_posterior_np{0:d}.bp".format(
nproc
),
posterior_samples,
) as vtx:
vtx.write(0.0)

eigen_decomposition_results = {"A": Hmisfit, "B": prior, "k": k, "d": d, "U": U}

final_results = {
"data_misfit_True": data_misfit_True,
Expand All @@ -235,15 +268,16 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]:
}

return final_results
#######################################


if __name__ == "__main__":
nx = 64
ny = 64
noise_variance = 1e-4
prior_param = {"gamma": 0.03, "delta": 0.3}

final_results = run_inversion(nx, ny, noise_variance, prior_param)

k, d = (
final_results["eigen_decomposition_results"]["k"],
final_results["eigen_decomposition_results"]["d"],
Expand Down
12 changes: 2 additions & 10 deletions example/poisson_dirichlet_example_reg.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -225,17 +225,9 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]:

hpx.parRandom.normal(1.0, Omega)

d, U = hpx.doublePassG(
Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1, check=False
)
d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1)

eigen_decomposition_results = {
"A": Hmisfit.mat,
"B": prior.R,
"k": k,
"d": d,
"U": U,
}
eigen_decomposition_results = {"A": Hmisfit, "B": prior, "k": k, "d": d, "U": U}

final_results = {
"data_misfit_True": data_misfit_True,
Expand Down
52 changes: 42 additions & 10 deletions example/poisson_example.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -195,17 +195,49 @@ def run_inversion(

hpx.parRandom.normal(1.0, Omega)

d, U = hpx.doublePassG(
Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1, check=False
)
d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1)

eigen_decomposition_results = {
"A": Hmisfit.mat,
"B": prior.R,
"k": k,
"d": d,
"U": U,
}
# generating prior and posterior samples
lap_aprx = hpx.LaplaceApproximator(prior, d, U)
lap_aprx.mean = prior.generate_parameter(0)
lap_aprx.mean.array[:] = x[hpx.PARAMETER].array[:]

m_prior = prior.generate_parameter(0)
m_post = prior.generate_parameter(0)

noise = prior.generate_parameter("noise")

num_samples_generate = 5

prior_samples = []
posterior_samples = []
for i in range(num_samples_generate):
hpx.parRandom.normal(1.0, noise)
lap_aprx.sample(noise, m_prior, m_post)
prior_sample = hpx.vector2Function(
m_prior, Vh[hpx.PARAMETER], name=f"prior_sample_{i}"
)
posterior_sample = hpx.vector2Function(
m_post, Vh[hpx.PARAMETER], name=f"posterior_sample_{i}"
)
prior_samples.append(prior_sample)
posterior_samples.append(posterior_sample)

with dlx.io.VTXWriter(
msh.comm,
"poisson_Robin_prior_Bilaplacian_samples_prior_np{0:d}.bp".format(nproc),
prior_samples,
) as vtx:
vtx.write(0.0)

with dlx.io.VTXWriter(
msh.comm,
"poisson_Robin_prior_Bilaplacian_samples_posterior_np{0:d}.bp".format(nproc),
posterior_samples,
) as vtx:
vtx.write(0.0)

eigen_decomposition_results = {"A": Hmisfit, "B": prior, "k": k, "d": d, "U": U}

final_results = {
"data_misfit_True": data_misfit_True,
Expand Down
12 changes: 2 additions & 10 deletions example/poisson_example_reg.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -198,17 +198,9 @@ def run_inversion(

hpx.parRandom.normal(1.0, Omega)

d, U = hpx.doublePassG(
Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1, check=False
)
d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1)

eigen_decomposition_results = {
"A": Hmisfit.mat,
"B": prior.R,
"k": k,
"d": d,
"U": U,
}
eigen_decomposition_results = {"A": Hmisfit, "B": prior, "k": k, "d": d, "U": U}

final_results = {
"data_misfit_True": data_misfit_True,
Expand Down
59 changes: 43 additions & 16 deletions example/sfsi_toy_gaussian.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -205,11 +205,6 @@ def run_inversion(
optimizer_results["optimizer"] = True
else:
optimizer_results["optimizer"] = False
final_results = {
"data_misfit_True": data_misfit_True,
"data_misfit_False": data_misfit_False,
"optimizer_results": optimizer_results,
}

Hmisfit = hpx.ReducedHessian(model, misfit_only=True)

Expand All @@ -228,17 +223,49 @@ def run_inversion(

hpx.parRandom.normal(1.0, Omega)

d, U = hpx.doublePassG(
Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1, check=False
)
d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1)

eigen_decomposition_results = {
"A": Hmisfit.mat,
"B": prior.R,
"k": k,
"d": d,
"U": U,
}
# generating prior and posterior samples
lap_aprx = hpx.LaplaceApproximator(prior, d, U)
lap_aprx.mean = prior.generate_parameter(0)
lap_aprx.mean.array[:] = x[hpx.PARAMETER].array[:]

m_prior = prior.generate_parameter(0)
m_post = prior.generate_parameter(0)

noise = prior.generate_parameter("noise")

num_samples_generate = 5

prior_samples = []
posterior_samples = []
for i in range(num_samples_generate):
hpx.parRandom.normal(1.0, noise)
lap_aprx.sample(noise, m_prior, m_post)
prior_sample = hpx.vector2Function(
m_prior, Vh[hpx.PARAMETER], name=f"prior_sample_{i}"
)
posterior_sample = hpx.vector2Function(
m_post, Vh[hpx.PARAMETER], name=f"posterior_sample_{i}"
)
prior_samples.append(prior_sample)
posterior_samples.append(posterior_sample)

with dlx.io.VTXWriter(
msh.comm,
"qpact_prior_Bilaplacian_samples_prior_np{0:d}.bp".format(nproc),
prior_samples,
) as vtx:
vtx.write(0.0)

with dlx.io.VTXWriter(
msh.comm,
"qpact_prior_Bilaplacian_samples_posterior_np{0:d}.bp".format(nproc),
posterior_samples,
) as vtx:
vtx.write(0.0)

eigen_decomposition_results = {"A": Hmisfit, "B": prior, "k": k, "d": d, "U": U}

final_results = {
"data_misfit_True": data_misfit_True,
Expand All @@ -248,7 +275,7 @@ def run_inversion(
}

return final_results
#######################################
######################################


if __name__ == "__main__":
Expand Down
12 changes: 2 additions & 10 deletions example/sfsi_toy_gaussian_reg.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -229,17 +229,9 @@ def run_inversion(

hpx.parRandom.normal(1.0, Omega)

d, U = hpx.doublePassG(
Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1, check=False
)
d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1)

eigen_decomposition_results = {
"A": Hmisfit.mat,
"B": prior.R,
"k": k,
"d": d,
"U": U,
}
eigen_decomposition_results = {"A": Hmisfit, "B": prior, "k": k, "d": d, "U": U}

final_results = {
"data_misfit_True": data_misfit_True,
Expand Down
5 changes: 3 additions & 2 deletions hippylibX/algorithms/NewtonCG.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -240,7 +240,7 @@ def _solve_ls(self, x: list) -> list:
solver.parameters["print_level"] = print_level - 1
mg_neg = self.model.generate_vector(PARAMETER)
mg_neg.array[:] = -1 * mg.array[:]
solver.solve(mhat, mg_neg)
solver.solve(mg_neg, mhat)
self.total_cg_iter += HessApply.ncalls
alpha = 1.0
descent = 0
Expand Down Expand Up @@ -374,7 +374,8 @@ def _solve_tr(self, x):
solver.parameters["zero_initial_guess"] = True
solver.parameters["print_level"] = print_level - 1

solver.solve(mhat, -mg)
# solver.solve(mhat, -mg)
solver.solve(-mg, mhat)
self.total_cg_iter += HessApply.ncalls
if self.it == 1:
self.model.prior.R.mult(mhat, R_mhat)
Expand Down
1 change: 1 addition & 0 deletions hippylibX/algorithms/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,3 +2,4 @@
from .NewtonCG import ReducedSpaceNewtonCG, ReducedSpaceNewtonCG_ParameterList # noqa
from .multivector import MultiVector, MatMvMult, MatMvTranspmult, MvDSmatMult # noqa
from .randomizedEigensolver import doublePassG # noqa
from .lowRankOperator import LowRankOperator # noqa
6 changes: 3 additions & 3 deletions hippylibX/algorithms/cgsolverSteihaug.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -160,14 +160,14 @@ def update_x_with_TR(self, x, alpha, d):
x.axpy(tau * alpha, d)
return True

def solve(self, x: dlx.la.Vector, b: dlx.la.Vector) -> None:
def solve(self, b: dlx.la.Vector, x: dlx.la.Vector) -> None:
temp_petsc_vec_x = dlx.la.create_petsc_vector_wrap(x)
temp_petsc_vec_b = dlx.la.create_petsc_vector_wrap(b)
self.solve_petsc(temp_petsc_vec_x, temp_petsc_vec_b)
self.solve_petsc(temp_petsc_vec_b, temp_petsc_vec_x)
temp_petsc_vec_x.destroy()
temp_petsc_vec_b.destroy()

def solve_petsc(self, x: petsc4py.PETSc.Vec, b: petsc4py.PETSc.Vec) -> None:
def solve_petsc(self, b: petsc4py.PETSc.Vec, x: petsc4py.PETSc.Vec) -> None:
"""
Solve the linear system :math:`Ax = b`
"""
Expand Down
Loading
Loading