From b5e53af2ed0c061192d5ea7a60c22ead2a58dff3 Mon Sep 17 00:00:00 2001 From: V-Rang Date: Wed, 10 Apr 2024 18:10:42 -0500 Subject: [PATCH 01/29] laplace approximation --- hippylibX/algorithms/__init__.py | 1 + hippylibX/algorithms/lowRankOperator.py | 120 ++++++++++ hippylibX/algorithms/multivector.py | 12 +- hippylibX/modeling/__init__.py | 5 + hippylibX/modeling/laplaceApproximation.py | 264 +++++++++++++++++++++ 5 files changed, 395 insertions(+), 7 deletions(-) create mode 100644 hippylibX/algorithms/lowRankOperator.py create mode 100644 hippylibX/modeling/laplaceApproximation.py diff --git a/hippylibX/algorithms/__init__.py b/hippylibX/algorithms/__init__.py index 7649993a..f198b04c 100644 --- a/hippylibX/algorithms/__init__.py +++ b/hippylibX/algorithms/__init__.py @@ -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 diff --git a/hippylibX/algorithms/lowRankOperator.py b/hippylibX/algorithms/lowRankOperator.py new file mode 100644 index 00000000..7debc550 --- /dev/null +++ b/hippylibX/algorithms/lowRankOperator.py @@ -0,0 +1,120 @@ +# Copyright (c) 2016-2018, The University of Texas at Austin +# & University of California--Merced. +# Copyright (c) 2019-2020, The University of Texas at Austin +# University of California--Merced, Washington University in St. Louis. +# +# All Rights reserved. +# See file COPYRIGHT for details. +# +# This file is part of the hIPPYlib library. For more information and source code +# availability see https://hippylib.github.io. +# +# hIPPYlib is free software; you can redistribute it and/or modify it under the +# terms of the GNU General Public License (as published by the Free +# Software Foundation) version 2.0 dated June 1991. + +import petsc4py.PETSc +from .multivector import MultiVector, MatMvMult +import numpy as np +import petsc4py + + +class LowRankOperator: + """ + This class model the action of a low rank operator :math:`A = U D U^T`. + Here :math:`D` is a diagonal matrix, and the columns of are orthonormal + in some weighted inner-product. + """ + + def __init__( + self, d: np.array, U: MultiVector, createVecLeft=None, createVecRight=None + ): + """ + Construct the low rank operator given :code:`d` and :code:`U`. + """ + self.d = d + self.U = U + self.createVecLeft = createVecLeft + self.createVecRight = createVecRight + + def __del__(self) -> None: + for i in range(self.U.nvec): + self.U[i].destroy() + + def mult(self, x: petsc4py.PETSc.Vec, y: petsc4py.PETSc.Vec) -> None: + """ + Compute :math:`y = Ax = U D U^T x` + """ + Utx = self.U.dot(x) + dUtx = self.d * Utx # elementwise mult + self.U.reduce(y, dUtx) + + def solve(self, sol: petsc4py.PETSc.Vec, rhs: petsc4py.PETSc.Vec) -> None: + """ + Compute :math:`\mbox{sol} = U D^-1 U^T x` + """ + Utr = self.U.dot(rhs) + dinvUtr = Utr / self.d + self.U.reduce(sol, dinvUtr) + + def get_diagonal(self, diag: petsc4py.PETSc.Vec) -> None: + """ + Compute the diagonal of :code:`A`. + """ + diag.scale(0.0) + tmp = self.U[0].duplicate() + for i in range(self.U.nvec): + tmp.scale(0.0) + tmp.axpy(1.0, self.U[i]) + tmp.pointwiseMult(tmp, self.U[i]) + diag.axpy(self.d[i], tmp) + + def trace(self, W=None) -> float: + """ + Compute the trace of :code:`A`. + If the weight :code:`W` is given, compute the trace of :math:`W^{1/2} A W^{1/2}`. + This is equivalent to :math:`\mbox{tr}_W(A) = \sum_i \lambda_i`, + where :math:`\lambda_i` are the generalized eigenvalues of + :math:`A x = \lambda W^{-1} x`. + + .. note:: If :math:`U` is a :math:`W`-orthogonal matrix then :math:`\mbox{tr}_W(A) = \sum_i D(i,i)`. + """ + if W is None: + tmp = self.U[0].duplicate() + tmp.scale(0.0) + self.U.reduce(tmp, np.sqrt(self.d)) + tr = tmp.dot(tmp) + else: + WU = MultiVector.createFromVec(self.U[0], self.U.nvec) + MatMvMult(W, self.U, WU) + diagWUtU = np.zeros_like(self.d) + for i in range(self.d.shape[0]): + diagWUtU[i] = WU[i].dot(self.U[i]) + tr = np.sum(self.d * diagWUtU) + + return tr + + def trace2(self, W=None) -> float: + """ + Compute the trace of :math:`A A` (Note this is the square of Frobenius norm, since :math:`A` is symmetic). + If the weight :code:`W` is provided, it will compute the trace of :math:`(AW)^2`. + + This is equivalent to :math:`\mbox{tr}_W(A) = \sum_i \lambda_i^2`, + where :math:`\lambda_i` are the generalized eigenvalues of + :math:`A x = \lambda W^{-1} x`. + + .. note:: If :math:`U` is a :math:`W`-orthogonal matrix then :math:`\mbox{tr}_W(A) = \sum_i D(i,i)^2`. + """ + if W is None: + UtU = self.U.dot(self.U) + dUtU = self.d[:, None] * UtU # diag(d)*UtU. + tr2 = np.sum(dUtU * dUtU) + else: + WU = MultiVector.createFromVec(self.U[0], self.U.nvec) + MatMvMult(W, self.U, WU) + WU = np.zeros(self.U.shape, dtype=self.U.dtype) + UtWU = self.U.dot(WU) + dUtWU = self.d[:, None] * UtWU # diag(d)*UtU. + tr2 = np.power(np.linalg.norm(dUtWU), 2) + + return tr2 diff --git a/hippylibX/algorithms/multivector.py b/hippylibX/algorithms/multivector.py index 67828d1a..9889c5e3 100644 --- a/hippylibX/algorithms/multivector.py +++ b/hippylibX/algorithms/multivector.py @@ -56,18 +56,15 @@ def dot(self, v: Union[petsc4py.PETSc.Vec, Type["MultiVector"]]) -> np.array: return return_values - def reduce(self, alpha: np.array) -> petsc4py.PETSc.Vec: - return_vec = self[0].duplicate() - return_vec.scale(0.0) + def reduce(self, y: petsc4py.PETSc.Vec, alpha: np.array) -> None: + y.scale(0.0) for i in range(self.nvec): - return_vec.axpy(alpha[i], self[i]) - return return_vec + y.axpy(alpha[i], self[i]) def axpy(self, alpha: Union[float, np.array], Y: Type["MultiVector"]) -> None: if isinstance(alpha, float): for i in range(self.nvec): self[i].axpy(alpha, Y[i]) - else: for i in range(self.nvec): self[i].axpy(alpha[i], Y[i]) @@ -154,5 +151,6 @@ def MvDSmatMult(X: MultiVector, A: np.array, Y: MultiVector) -> None: ), "Y Number of vecs incompatible with number of cols in A" for j in range(Y.nvec): Y[j].scale(0.0) - reduced_vec = X.reduce(A[:, j].flatten()) + reduced_vec = X[0].duplicate() + X.reduce(reduced_vec, A[:, j].flatten()) Y.data[j] = X.data[0].duplicate(reduced_vec.getArray()) diff --git a/hippylibX/modeling/__init__.py b/hippylibX/modeling/__init__.py index 4493c511..57a5493e 100644 --- a/hippylibX/modeling/__init__.py +++ b/hippylibX/modeling/__init__.py @@ -6,3 +6,8 @@ from .Regularization import H1TikhonvFunctional, VariationalRegularization # noqa from .variables import STATE, PARAMETER, ADJOINT, NVAR # noqa from .reducedHessian import ReducedHessian # noqa +from .laplaceApproximation import ( + LowRankHessian, # noqa + LowRankPosteriorSampler, # noqa + LaplaceApproximator, # noqa +) diff --git a/hippylibX/modeling/laplaceApproximation.py b/hippylibX/modeling/laplaceApproximation.py new file mode 100644 index 00000000..ae0b66f6 --- /dev/null +++ b/hippylibX/modeling/laplaceApproximation.py @@ -0,0 +1,264 @@ +# Copyright (c) 2016-2018, The University of Texas at Austin +# & University of California--Merced. +# Copyright (c) 2019-2020, The University of Texas at Austin +# University of California--Merced, Washington University in St. Louis. +# +# All Rights reserved. +# See file COPYRIGHT for details. +# +# This file is part of the hIPPYlib library. For more information and source code +# availability see https://hippylib.github.io. +# +# hIPPYlib is free software; you can redistribute it and/or modify it under the +# terms of the GNU General Public License (as published by the Free +# Software Foundation) version 2.0 dated June 1991. + +from ..algorithms.lowRankOperator import LowRankOperator +from ..algorithms.multivector import MultiVector +import numpy as np +import petsc4py +import dolfinx as dlx + + +def not_implemented(func): + def wrapper(*args, **kwargs): + raise NotImplementedError(f"{func.__name__} Function has not been implemented.") + + return wrapper + + +class LowRankHessian: + """ + Operator that represents the action of the low rank approximation + of the Hessian and of its inverse. + """ + + def __init__( + self, + prior, + d: np.array, + U: MultiVector, + createVecLeft=None, + createVecRight=None, + ): + self.prior = prior + self.LowRankH = LowRankOperator(d, U) + dsolve = d / (np.ones(d.shape, dtype=d.dtype) + d) + self.LowRankHinv = LowRankOperator(dsolve, U) + + self.createVecLeft = createVecLeft + self.createVecRight = createVecRight + + self.help = self.prior.R.createVecRight() + self.help1 = self.prior.R.createVecLeft() + + def __del__(self) -> None: + for i in range(self.U.nvec): + self.U[i].destroy() + + def mult(self, x: petsc4py.PETSc.Vec, y: petsc4py.PETSc.Vec) -> None: + self.prior.R.mult(x, y) + self.LowRankH.mult(y, self.help) + self.prior.R.mult(self.help, self.help1) + y.axpy(1, self.help1) + + def solve(self, sol: petsc4py.PETSc.Vec, rhs: petsc4py.PETSc.Vec) -> None: + # self.prior.Rsolver.solve(sol, rhs) + self.prior.Rsolver.solve(rhs, sol) + self.LowRankHinv.mult(rhs, self.help) + sol.axpy(-1, self.help) + + +class LowRankPosteriorSampler: + """ + Object to sample from the low rank approximation + of the posterior. + + .. math:: y = ( I - U S U^{T}) x, + + where + + :math:`S = I - (I + D)^{-1/2}, x \\sim \\mathcal{N}(0, R^{-1}).` + """ + + def __init__( + self, + prior, + d: np.array, + U: MultiVector, + createVecLeft=None, + createVecRight=None, + ): + self.prior = prior + self.createVecLeft = createVecLeft + self.createVecRight = createVecRight + + ones = np.ones(d.shape, dtype=d.dtype) + self.d = ones - np.power(ones + d, -0.5) + self.lrsqrt = LowRankOperator(self.d, U) + self.help = self.prior.R.createVecLeft() + + def __del__(self) -> None: + for i in range(self.U.nvec): + self.U[i].destroy() + + def sample(self, noise: dlx.la.Vector, s: dlx.la.Vector): + temp_petsc_vec_noise = dlx.la.create_petsc_vector_wrap(noise) + temp_petsc_vec_s = dlx.la.create_petsc_vector_wrap(s) + + self.prior.R.mult(temp_petsc_vec_noise, self.help) + self.lrsqrt.mult(self.help, temp_petsc_vec_s) + temp_petsc_vec_s.axpy(-1.0, temp_petsc_vec_noise) + temp_petsc_vec_s.scale(-1.0) + temp_petsc_vec_noise.destroy() + temp_petsc_vec_s.destroy() + + +class LaplaceApproximator: + """ + Class for the low rank Gaussian Approximation of the Posterior. + This class provides functionality for approximate Hessian + apply, solve, and Gaussian sampling based on the low rank + factorization of the Hessian. + + In particular if :math:`d` and :math:`U` are the dominant eigenpairs of + :math:`H_{\\mbox{misfit}} U[:,i] = d[i] R U[:,i]` + then we have: + + - low rank Hessian apply: :math:`y = ( R + RU D U^{T}) x.` + - low rank Hessian solve: :math:`y = (R^-1 - U (I + D^{-1})^{-1} U^T) x.` + - low rank Hessian Gaussian sampling: :math:`y = ( I - U S U^{T}) x`, where :math:`S = I - (I + D)^{-1/2}` and :math:`x \\sim \\mathcal{N}(0, R^{-1}).` + """ + + def __init__(self, prior, d: np.array, U: MultiVector, mean=None): + """ + Construct the Gaussian approximation of the posterior. + Input: + - :code:`prior`: the prior mode. + - :code:`d`: the dominant generalized eigenvalues of the Hessian misfit. + - :code:`U`: the dominant generalized eigenvector of the Hessian misfit :math:`U^T R U = I.` + - :code:`mean`: the MAP point. + """ + self.prior = prior + self.d = d + self.U = U + self.Hlr = LowRankHessian(prior, d, U) + self.sampler = LowRankPosteriorSampler(self.prior, self.d, self.U) + self.mean = None + + def __del__(self) -> None: + for i in range(self.U.nvec): + self.U[i].destroy() + + def cost(self, m: dlx.la.Vector) -> float: + if self.mean is None: + dm = m + else: + dm = m - self.mean + temp_petsc_vec_dm = dlx.la.create_petsc_vector_wrap(dm) + self.Hlr.mult(temp_petsc_vec_dm, self.Hlr.help1) + return_value = 0.5 * self.Hlr.dot(self.Hlr.help1) + temp_petsc_vec_dm.destroy() + return return_value + + def sample(self, *args, **kwargs): + """ + possible calls: + + 1) :code:`sample(s_prior, s_post, add_mean=True)` + + Given a prior sample :code:`s_prior` compute a sample :code:`s_post` from the posterior. + + - :code:`s_prior` is a sample from the prior centered at 0 (input). + - :code:`s_post` is a sample from the posterior (output). + - if :code:`add_mean=True` (default) then the samples will be centered at the map point. + + 2) :code:`sample(noise, s_prior, s_post, add_mean=True)` + + Given :code:`noise` :math:`\\sim \\mathcal{N}(0, I)` compute a sample :code:`s_prior` from the prior and + :code:`s_post` from the posterior. + + - :code:`noise` is a realization of white noise (input). + - :code:`s_prior` is a sample from the prior (output). + - :code:`s_post` is a sample from the posterior. + - if :code:`add_mean=True` (default) then the prior and posterior samples will be centered at the respective means. + + """ + add_mean = True + for name, value in kwargs.items(): + if name == "add_mean": + add_mean = value + else: + raise NameError(name) + + if len(args) == 2: + self._sample_given_prior(args[0], args[1]) + if add_mean: + args[1].axpy(1.0, self.mean) + elif len(args) == 3: + self._sample_given_white_noise(args[0], args[1], args[2]) + if add_mean: + args[1].axpy(1.0, self.prior.mean) + args[2].axpy(1.0, self.mean) + else: + raise NameError("Invalid number of parameters in Posterior::sample") + + def _sample_given_white_noise( + self, noise: dlx.la.Vector, s_prior: dlx.la.Vector, s_post: dlx.la.Vector + ): + self.prior.sample(noise, s_prior, add_mean=False) + self.sampler.sample(s_prior, s_post) + + def _sample_given_prior(self, s_prior: dlx.la.Vector, s_post: dlx.la.Vector): + self.sampler.sample(s_prior, s_post) + + @not_implemented + def trace(self, **kwargs): + """ + Compute/estimate the trace of the posterior, prior distribution + and the trace of the data informed correction. + + See :code:`_Prior.trace` for more details. + """ + pr_trace = self.prior.trace(**kwargs) + corr_trace = self.trace_update() + post_trace = pr_trace - corr_trace + return post_trace, pr_trace, corr_trace + + @not_implemented + def trace_update(self): + return self.Hlr.LowRankHinv.trace(self.prior.M) + + @not_implemented + def pointwise_variance(self, **kwargs): + """ + Compute/estimate the pointwise variance of the posterior, prior distribution + and the pointwise variance reduction informed by the data. + + See :code:`_Prior.pointwise_variance` for more details. + """ + pr_pointwise_variance = self.prior.pointwise_variance(**kwargs) + # correction_pointwise_variance = Vector(self.prior.R.mpi_comm()) + correction_pointwise_variance = None + self.init_vector(correction_pointwise_variance, 0) + self.Hlr.LowRankHinv.get_diagonal(correction_pointwise_variance) + post_pointwise_variance = pr_pointwise_variance - correction_pointwise_variance + return ( + post_pointwise_variance, + pr_pointwise_variance, + correction_pointwise_variance, + ) + + def klDistanceFromPrior(self, sub_comp=False): + dplus1 = self.d + np.ones_like(self.d) + + c_logdet = 0.5 * np.sum(np.log(dplus1)) + c_trace = -0.5 * np.sum(self.d / dplus1) + c_shift = self.prior.cost(self.mean) + + kld = c_logdet + c_trace + c_shift + + if sub_comp: + return kld, c_logdet, c_trace, c_shift + else: + return kld From e705e03dc69b0b14a05035537872b5be3a835665 Mon Sep 17 00:00:00 2001 From: V-Rang Date: Wed, 10 Apr 2024 18:16:07 -0500 Subject: [PATCH 02/29] update __del__ method --- hippylibX/modeling/laplaceApproximation.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/hippylibX/modeling/laplaceApproximation.py b/hippylibX/modeling/laplaceApproximation.py index ae0b66f6..89683104 100644 --- a/hippylibX/modeling/laplaceApproximation.py +++ b/hippylibX/modeling/laplaceApproximation.py @@ -55,6 +55,9 @@ def __init__( def __del__(self) -> None: for i in range(self.U.nvec): self.U[i].destroy() + + self.help.destroy() + self.help1.destroy() def mult(self, x: petsc4py.PETSc.Vec, y: petsc4py.PETSc.Vec) -> None: self.prior.R.mult(x, y) @@ -63,7 +66,6 @@ def mult(self, x: petsc4py.PETSc.Vec, y: petsc4py.PETSc.Vec) -> None: y.axpy(1, self.help1) def solve(self, sol: petsc4py.PETSc.Vec, rhs: petsc4py.PETSc.Vec) -> None: - # self.prior.Rsolver.solve(sol, rhs) self.prior.Rsolver.solve(rhs, sol) self.LowRankHinv.mult(rhs, self.help) sol.axpy(-1, self.help) @@ -101,7 +103,9 @@ def __init__( def __del__(self) -> None: for i in range(self.U.nvec): self.U[i].destroy() - + + self.help.destroy() + def sample(self, noise: dlx.la.Vector, s: dlx.la.Vector): temp_petsc_vec_noise = dlx.la.create_petsc_vector_wrap(noise) temp_petsc_vec_s = dlx.la.create_petsc_vector_wrap(s) From 0ac0944e2cca6e502e904c07e03e74d866c25d26 Mon Sep 17 00:00:00 2001 From: V-Rang Date: Wed, 10 Apr 2024 21:01:30 -0500 Subject: [PATCH 03/29] minor changes --- hippylibX/algorithms/lowRankOperator.py | 25 ---------------------- hippylibX/modeling/laplaceApproximation.py | 16 ++++++++------ 2 files changed, 9 insertions(+), 32 deletions(-) diff --git a/hippylibX/algorithms/lowRankOperator.py b/hippylibX/algorithms/lowRankOperator.py index 7debc550..9e18a699 100644 --- a/hippylibX/algorithms/lowRankOperator.py +++ b/hippylibX/algorithms/lowRankOperator.py @@ -93,28 +93,3 @@ def trace(self, W=None) -> float: tr = np.sum(self.d * diagWUtU) return tr - - def trace2(self, W=None) -> float: - """ - Compute the trace of :math:`A A` (Note this is the square of Frobenius norm, since :math:`A` is symmetic). - If the weight :code:`W` is provided, it will compute the trace of :math:`(AW)^2`. - - This is equivalent to :math:`\mbox{tr}_W(A) = \sum_i \lambda_i^2`, - where :math:`\lambda_i` are the generalized eigenvalues of - :math:`A x = \lambda W^{-1} x`. - - .. note:: If :math:`U` is a :math:`W`-orthogonal matrix then :math:`\mbox{tr}_W(A) = \sum_i D(i,i)^2`. - """ - if W is None: - UtU = self.U.dot(self.U) - dUtU = self.d[:, None] * UtU # diag(d)*UtU. - tr2 = np.sum(dUtU * dUtU) - else: - WU = MultiVector.createFromVec(self.U[0], self.U.nvec) - MatMvMult(W, self.U, WU) - WU = np.zeros(self.U.shape, dtype=self.U.dtype) - UtWU = self.U.dot(WU) - dUtWU = self.d[:, None] * UtWU # diag(d)*UtU. - tr2 = np.power(np.linalg.norm(dUtWU), 2) - - return tr2 diff --git a/hippylibX/modeling/laplaceApproximation.py b/hippylibX/modeling/laplaceApproximation.py index 89683104..74c18abe 100644 --- a/hippylibX/modeling/laplaceApproximation.py +++ b/hippylibX/modeling/laplaceApproximation.py @@ -41,10 +41,11 @@ def __init__( createVecLeft=None, createVecRight=None, ): + self.U = U self.prior = prior - self.LowRankH = LowRankOperator(d, U) + self.LowRankH = LowRankOperator(d, self.U) dsolve = d / (np.ones(d.shape, dtype=d.dtype) + d) - self.LowRankHinv = LowRankOperator(dsolve, U) + self.LowRankHinv = LowRankOperator(dsolve, self.U) self.createVecLeft = createVecLeft self.createVecRight = createVecRight @@ -55,7 +56,7 @@ def __init__( def __del__(self) -> None: for i in range(self.U.nvec): self.U[i].destroy() - + self.help.destroy() self.help1.destroy() @@ -91,21 +92,22 @@ def __init__( createVecLeft=None, createVecRight=None, ): + self.U = U self.prior = prior self.createVecLeft = createVecLeft self.createVecRight = createVecRight ones = np.ones(d.shape, dtype=d.dtype) self.d = ones - np.power(ones + d, -0.5) - self.lrsqrt = LowRankOperator(self.d, U) + self.lrsqrt = LowRankOperator(self.d, self.U) self.help = self.prior.R.createVecLeft() def __del__(self) -> None: for i in range(self.U.nvec): self.U[i].destroy() - + self.help.destroy() - + def sample(self, noise: dlx.la.Vector, s: dlx.la.Vector): temp_petsc_vec_noise = dlx.la.create_petsc_vector_wrap(noise) temp_petsc_vec_s = dlx.la.create_petsc_vector_wrap(s) @@ -161,7 +163,7 @@ def cost(self, m: dlx.la.Vector) -> float: dm = m - self.mean temp_petsc_vec_dm = dlx.la.create_petsc_vector_wrap(dm) self.Hlr.mult(temp_petsc_vec_dm, self.Hlr.help1) - return_value = 0.5 * self.Hlr.dot(self.Hlr.help1) + return_value = 0.5 * self.Hlr.help1.dot(temp_petsc_vec_dm) temp_petsc_vec_dm.destroy() return return_value From 68a5c03fc15908c26d52570b84f7844b39d2b07b Mon Sep 17 00:00:00 2001 From: V-Rang Date: Thu, 11 Apr 2024 08:28:23 -0500 Subject: [PATCH 04/29] minor comment fixes --- example/poisson_dirichlet_example.py | 4 +--- example/poisson_dirichlet_example_reg.py | 4 +--- example/poisson_example.py | 4 +--- example/poisson_example_reg.py | 4 +--- example/sfsi_toy_gaussian.py | 4 +--- example/sfsi_toy_gaussian_reg.py | 4 +--- hippylibX/algorithms/randomizedEigensolver.py | 1 - 7 files changed, 6 insertions(+), 19 deletions(-) diff --git a/example/poisson_dirichlet_example.py b/example/poisson_dirichlet_example.py index af19e6d2..211c4783 100644 --- a/example/poisson_dirichlet_example.py +++ b/example/poisson_dirichlet_example.py @@ -215,9 +215,7 @@ 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, diff --git a/example/poisson_dirichlet_example_reg.py b/example/poisson_dirichlet_example_reg.py index b11d7439..a5a641d2 100644 --- a/example/poisson_dirichlet_example_reg.py +++ b/example/poisson_dirichlet_example_reg.py @@ -225,9 +225,7 @@ 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, diff --git a/example/poisson_example.py b/example/poisson_example.py index 162d575d..18364e43 100644 --- a/example/poisson_example.py +++ b/example/poisson_example.py @@ -195,9 +195,7 @@ 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, diff --git a/example/poisson_example_reg.py b/example/poisson_example_reg.py index 0da1e975..e17c7a0e 100644 --- a/example/poisson_example_reg.py +++ b/example/poisson_example_reg.py @@ -198,9 +198,7 @@ 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, diff --git a/example/sfsi_toy_gaussian.py b/example/sfsi_toy_gaussian.py index 742dbbbd..1522482c 100644 --- a/example/sfsi_toy_gaussian.py +++ b/example/sfsi_toy_gaussian.py @@ -228,9 +228,7 @@ 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, diff --git a/example/sfsi_toy_gaussian_reg.py b/example/sfsi_toy_gaussian_reg.py index 05e829d3..2f56affa 100644 --- a/example/sfsi_toy_gaussian_reg.py +++ b/example/sfsi_toy_gaussian_reg.py @@ -229,9 +229,7 @@ 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, diff --git a/hippylibX/algorithms/randomizedEigensolver.py b/hippylibX/algorithms/randomizedEigensolver.py index 5f4e8d24..4af6c539 100644 --- a/hippylibX/algorithms/randomizedEigensolver.py +++ b/hippylibX/algorithms/randomizedEigensolver.py @@ -13,7 +13,6 @@ def doublePassG( Omega: MultiVector, k: int, s=1, - check=False, ) -> tuple[np.array, MultiVector]: nvec = Omega.nvec assert nvec >= k From c65c361ae158b44778f2e89f0fa7c0bf6485c6e2 Mon Sep 17 00:00:00 2001 From: V-Rang Date: Thu, 11 Apr 2024 11:12:41 -0500 Subject: [PATCH 05/29] changes to del and solve method --- hippylibX/algorithms/lowRankOperator.py | 12 +++--------- hippylibX/algorithms/multivector.py | 9 +-------- hippylibX/modeling/laplaceApproximation.py | 12 +----------- 3 files changed, 5 insertions(+), 28 deletions(-) diff --git a/hippylibX/algorithms/lowRankOperator.py b/hippylibX/algorithms/lowRankOperator.py index 9e18a699..d005b331 100644 --- a/hippylibX/algorithms/lowRankOperator.py +++ b/hippylibX/algorithms/lowRankOperator.py @@ -37,10 +37,6 @@ def __init__( self.createVecLeft = createVecLeft self.createVecRight = createVecRight - def __del__(self) -> None: - for i in range(self.U.nvec): - self.U[i].destroy() - def mult(self, x: petsc4py.PETSc.Vec, y: petsc4py.PETSc.Vec) -> None: """ Compute :math:`y = Ax = U D U^T x` @@ -49,7 +45,7 @@ def mult(self, x: petsc4py.PETSc.Vec, y: petsc4py.PETSc.Vec) -> None: dUtx = self.d * Utx # elementwise mult self.U.reduce(y, dUtx) - def solve(self, sol: petsc4py.PETSc.Vec, rhs: petsc4py.PETSc.Vec) -> None: + def solve(self, rhs: petsc4py.PETSc.Vec, sol: petsc4py.PETSc.Vec) -> None: """ Compute :math:`\mbox{sol} = U D^-1 U^T x` """ @@ -63,10 +59,8 @@ def get_diagonal(self, diag: petsc4py.PETSc.Vec) -> None: """ diag.scale(0.0) tmp = self.U[0].duplicate() - for i in range(self.U.nvec): - tmp.scale(0.0) - tmp.axpy(1.0, self.U[i]) - tmp.pointwiseMult(tmp, self.U[i]) + for i in range(self.U.nvec): + tmp.pointwiseMult(self.U[i], self.U[i]) diag.axpy(self.d[i], tmp) def trace(self, W=None) -> float: diff --git a/hippylibX/algorithms/multivector.py b/hippylibX/algorithms/multivector.py index 9889c5e3..75e37e20 100644 --- a/hippylibX/algorithms/multivector.py +++ b/hippylibX/algorithms/multivector.py @@ -27,10 +27,6 @@ def createFromMultiVec(cls, mv: Type["MultiVector"]) -> Type["MultiVector"]: return mv_copy - def __del__(self) -> None: - for d in self.data: - d.destroy() - def __getitem__(self, k: int) -> petsc4py.PETSc.Vec: return self.data[k] @@ -57,7 +53,6 @@ def dot(self, v: Union[petsc4py.PETSc.Vec, Type["MultiVector"]]) -> np.array: return return_values def reduce(self, y: petsc4py.PETSc.Vec, alpha: np.array) -> None: - y.scale(0.0) for i in range(self.nvec): y.axpy(alpha[i], self[i]) @@ -151,6 +146,4 @@ def MvDSmatMult(X: MultiVector, A: np.array, Y: MultiVector) -> None: ), "Y Number of vecs incompatible with number of cols in A" for j in range(Y.nvec): Y[j].scale(0.0) - reduced_vec = X[0].duplicate() - X.reduce(reduced_vec, A[:, j].flatten()) - Y.data[j] = X.data[0].duplicate(reduced_vec.getArray()) + X.reduce(Y[j],A[:,j].flatten()) \ No newline at end of file diff --git a/hippylibX/modeling/laplaceApproximation.py b/hippylibX/modeling/laplaceApproximation.py index 74c18abe..95c48123 100644 --- a/hippylibX/modeling/laplaceApproximation.py +++ b/hippylibX/modeling/laplaceApproximation.py @@ -54,9 +54,6 @@ def __init__( self.help1 = self.prior.R.createVecLeft() def __del__(self) -> None: - for i in range(self.U.nvec): - self.U[i].destroy() - self.help.destroy() self.help1.destroy() @@ -66,7 +63,7 @@ def mult(self, x: petsc4py.PETSc.Vec, y: petsc4py.PETSc.Vec) -> None: self.prior.R.mult(self.help, self.help1) y.axpy(1, self.help1) - def solve(self, sol: petsc4py.PETSc.Vec, rhs: petsc4py.PETSc.Vec) -> None: + def solve(self, rhs: petsc4py.PETSc.Vec, sol: petsc4py.PETSc.Vec) -> None: self.prior.Rsolver.solve(rhs, sol) self.LowRankHinv.mult(rhs, self.help) sol.axpy(-1, self.help) @@ -103,9 +100,6 @@ def __init__( self.help = self.prior.R.createVecLeft() def __del__(self) -> None: - for i in range(self.U.nvec): - self.U[i].destroy() - self.help.destroy() def sample(self, noise: dlx.la.Vector, s: dlx.la.Vector): @@ -152,10 +146,6 @@ def __init__(self, prior, d: np.array, U: MultiVector, mean=None): self.sampler = LowRankPosteriorSampler(self.prior, self.d, self.U) self.mean = None - def __del__(self) -> None: - for i in range(self.U.nvec): - self.U[i].destroy() - def cost(self, m: dlx.la.Vector) -> float: if self.mean is None: dm = m From 2f14024ad509ecccf1a03958a1119c6a870804c5 Mon Sep 17 00:00:00 2001 From: V-Rang Date: Thu, 11 Apr 2024 11:14:53 -0500 Subject: [PATCH 06/29] formatting fix --- hippylibX/algorithms/lowRankOperator.py | 2 +- hippylibX/algorithms/multivector.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/hippylibX/algorithms/lowRankOperator.py b/hippylibX/algorithms/lowRankOperator.py index d005b331..c93f4965 100644 --- a/hippylibX/algorithms/lowRankOperator.py +++ b/hippylibX/algorithms/lowRankOperator.py @@ -59,7 +59,7 @@ def get_diagonal(self, diag: petsc4py.PETSc.Vec) -> None: """ diag.scale(0.0) tmp = self.U[0].duplicate() - for i in range(self.U.nvec): + for i in range(self.U.nvec): tmp.pointwiseMult(self.U[i], self.U[i]) diag.axpy(self.d[i], tmp) diff --git a/hippylibX/algorithms/multivector.py b/hippylibX/algorithms/multivector.py index 75e37e20..c01ab402 100644 --- a/hippylibX/algorithms/multivector.py +++ b/hippylibX/algorithms/multivector.py @@ -146,4 +146,4 @@ def MvDSmatMult(X: MultiVector, A: np.array, Y: MultiVector) -> None: ), "Y Number of vecs incompatible with number of cols in A" for j in range(Y.nvec): Y[j].scale(0.0) - X.reduce(Y[j],A[:,j].flatten()) \ No newline at end of file + X.reduce(Y[j], A[:, j].flatten()) From db414205b63afc504a08beae97f0f07afd655384 Mon Sep 17 00:00:00 2001 From: V-Rang Date: Thu, 11 Apr 2024 21:40:50 -0500 Subject: [PATCH 07/29] prior and posterior samples --- example/poisson_dirichlet_example.py | 33 ++++++++++++++++ example/sfsi_toy_gaussian.py | 44 +++++++++++++++++++--- hippylibX/algorithms/lowRankOperator.py | 2 + hippylibX/modeling/Regularization.py | 3 ++ hippylibX/modeling/laplaceApproximation.py | 28 ++++++++++++-- 5 files changed, 100 insertions(+), 10 deletions(-) diff --git a/example/poisson_dirichlet_example.py b/example/poisson_dirichlet_example.py index 211c4783..80401aba 100644 --- a/example/poisson_dirichlet_example.py +++ b/example/poisson_dirichlet_example.py @@ -217,6 +217,39 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1) + lap_aprx = hpx.LaplaceApproximator(prior, d, U) + lap_aprx.mean = prior.generate_parameter(0) + lap_aprx.mean.array[:] = prior_mean.array[:] + + noise = prior.generate_parameter("noise") + hpx.parRandom.normal(1.0, noise) + + m0 = prior.generate_parameter(0) + m_post = prior.generate_parameter(0) + + lap_aprx.sample(noise, m0, m_post) + + true_param = hpx.vector2Function(m_true, Vh[hpx.PARAMETER]) + with dlx.io.XDMFFile( + msh.comm, "dirichlet_Poisson_true_para_np{0:d}_X.xdmf".format(nproc), "w" + ) as file: + file.write_mesh(msh) + file.write_function(true_param) + + prior_sample = hpx.vector2Function(m0, Vh[hpx.PARAMETER]) + with dlx.io.XDMFFile( + msh.comm, "dirichlet_Poisson_prior_sample_np{0:d}_X.xdmf".format(nproc), "w" + ) as file: + file.write_mesh(msh) + file.write_function(prior_sample) + + posterior_sample = hpx.vector2Function(m_post, Vh[hpx.PARAMETER]) + with dlx.io.XDMFFile( + msh.comm, "dirichlet_Poisson_posterior_sample_np{0:d}_X.xdmf".format(nproc), "w" + ) as file: + file.write_mesh(msh) + file.write_function(posterior_sample) + eigen_decomposition_results = { "A": Hmisfit.mat, "B": prior.R, diff --git a/example/sfsi_toy_gaussian.py b/example/sfsi_toy_gaussian.py index 1522482c..3d7a594d 100644 --- a/example/sfsi_toy_gaussian.py +++ b/example/sfsi_toy_gaussian.py @@ -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) @@ -230,6 +225,43 @@ def run_inversion( d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1) + lap_aprx = hpx.LaplaceApproximator( + prior, + d, + U, + ) + lap_aprx.mean = prior.generate_parameter(0) + lap_aprx.mean.array[:] = prior_mean.array[:] + + noise = prior.generate_parameter("noise") + hpx.parRandom.normal(1.0, noise) + + m0 = prior.generate_parameter(0) + m_post = prior.generate_parameter(0) + + lap_aprx.sample(noise, m0, m_post) + + true_param = hpx.vector2Function(m_true, Vh[hpx.PARAMETER]) + with dlx.io.XDMFFile( + msh.comm, "qpact_true_para_np{0:d}_X.xdmf".format(nproc), "w" + ) as file: + file.write_mesh(msh) + file.write_function(true_param) + + prior_sample = hpx.vector2Function(m0, Vh[hpx.PARAMETER]) + with dlx.io.XDMFFile( + msh.comm, "qpact_prior_sample_np{0:d}_X.xdmf".format(nproc), "w" + ) as file: + file.write_mesh(msh) + file.write_function(prior_sample) + + posterior_sample = hpx.vector2Function(m_post, Vh[hpx.PARAMETER]) + with dlx.io.XDMFFile( + msh.comm, "qpact_posterior_sample_np{0:d}_X.xdmf".format(nproc), "w" + ) as file: + file.write_mesh(msh) + file.write_function(posterior_sample) + eigen_decomposition_results = { "A": Hmisfit.mat, "B": prior.R, @@ -246,7 +278,7 @@ def run_inversion( } return final_results - ####################################### + ###################################### if __name__ == "__main__": diff --git a/hippylibX/algorithms/lowRankOperator.py b/hippylibX/algorithms/lowRankOperator.py index c93f4965..803e7503 100644 --- a/hippylibX/algorithms/lowRankOperator.py +++ b/hippylibX/algorithms/lowRankOperator.py @@ -43,6 +43,7 @@ def mult(self, x: petsc4py.PETSc.Vec, y: petsc4py.PETSc.Vec) -> None: """ Utx = self.U.dot(x) dUtx = self.d * Utx # elementwise mult + y.scale(0.0) self.U.reduce(y, dUtx) def solve(self, rhs: petsc4py.PETSc.Vec, sol: petsc4py.PETSc.Vec) -> None: @@ -51,6 +52,7 @@ def solve(self, rhs: petsc4py.PETSc.Vec, sol: petsc4py.PETSc.Vec) -> None: """ Utr = self.U.dot(rhs) dinvUtr = Utr / self.d + sol.scale(0.0) self.U.reduce(sol, dinvUtr) def get_diagonal(self, diag: petsc4py.PETSc.Vec) -> None: diff --git a/hippylibX/modeling/Regularization.py b/hippylibX/modeling/Regularization.py index ffea59a2..37408d61 100644 --- a/hippylibX/modeling/Regularization.py +++ b/hippylibX/modeling/Regularization.py @@ -77,6 +77,9 @@ def __del__(self): self.Msolver.destroy() self.M.destroy() + def generate_parameter(self, dim: int) -> dlx.la.Vector: + return dlx.la.Vector(self.Vh.dofmap.index_map) + def _createsolver(self, petsc_options: dict) -> petsc4py.PETSc.KSP: ksp = petsc4py.PETSc.KSP().create(self.Vh.mesh.comm) problem_prefix = f"dolfinx_solve_{id(self)}" diff --git a/hippylibX/modeling/laplaceApproximation.py b/hippylibX/modeling/laplaceApproximation.py index 95c48123..2e2504f4 100644 --- a/hippylibX/modeling/laplaceApproximation.py +++ b/hippylibX/modeling/laplaceApproximation.py @@ -144,7 +144,10 @@ def __init__(self, prior, d: np.array, U: MultiVector, mean=None): self.U = U self.Hlr = LowRankHessian(prior, d, U) self.sampler = LowRankPosteriorSampler(self.prior, self.d, self.U) - self.mean = None + self.mean = mean + + if self.mean is None: + self.mean = self.prior.generate_parameter(0) def cost(self, m: dlx.la.Vector) -> float: if self.mean is None: @@ -190,12 +193,29 @@ def sample(self, *args, **kwargs): if len(args) == 2: self._sample_given_prior(args[0], args[1]) if add_mean: - args[1].axpy(1.0, self.mean) + temp_petsc_vec_arg1 = dlx.la.create_petsc_vector_wrap(args[1]) + temp_petsc_vec_mean = dlx.la.create_petsc_vector_wrap(self.mean) + temp_petsc_vec_arg1.axpy(1.0, temp_petsc_vec_mean) + temp_petsc_vec_arg1.destroy() + temp_petsc_vec_mean.destroy() + elif len(args) == 3: self._sample_given_white_noise(args[0], args[1], args[2]) if add_mean: - args[1].axpy(1.0, self.prior.mean) - args[2].axpy(1.0, self.mean) + temp_petsc_vec_arg1 = dlx.la.create_petsc_vector_wrap(args[1]) + temp_petsc_vec_prior_mean = dlx.la.create_petsc_vector_wrap( + self.prior.mean + ) + temp_petsc_vec_arg2 = dlx.la.create_petsc_vector_wrap(args[2]) + temp_petsc_vec_mean = dlx.la.create_petsc_vector_wrap(self.mean) + + temp_petsc_vec_arg1.axpy(1.0, temp_petsc_vec_prior_mean) + temp_petsc_vec_arg2.axpy(1.0, temp_petsc_vec_mean) + + temp_petsc_vec_arg1.destroy() + temp_petsc_vec_prior_mean.destroy() + temp_petsc_vec_arg2.destroy() + temp_petsc_vec_mean.destroy() else: raise NameError("Invalid number of parameters in Posterior::sample") From d06de44cda1dfef6c240cd91ff7c0cef90f874bd Mon Sep 17 00:00:00 2001 From: V-Rang Date: Fri, 12 Apr 2024 06:53:41 -0500 Subject: [PATCH 08/29] change posterior mean, formatting --- example/poisson_dirichlet_example.py | 8 ++++---- example/sfsi_toy_gaussian.py | 14 +++++--------- 2 files changed, 9 insertions(+), 13 deletions(-) diff --git a/example/poisson_dirichlet_example.py b/example/poisson_dirichlet_example.py index 80401aba..fb8c35e9 100644 --- a/example/poisson_dirichlet_example.py +++ b/example/poisson_dirichlet_example.py @@ -219,15 +219,15 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: lap_aprx = hpx.LaplaceApproximator(prior, d, U) lap_aprx.mean = prior.generate_parameter(0) - lap_aprx.mean.array[:] = prior_mean.array[:] + lap_aprx.mean.array[:] = x[hpx.PARAMETER].array[:] noise = prior.generate_parameter("noise") hpx.parRandom.normal(1.0, noise) - m0 = prior.generate_parameter(0) + m_prior = prior.generate_parameter(0) m_post = prior.generate_parameter(0) - lap_aprx.sample(noise, m0, m_post) + lap_aprx.sample(noise, m_prior, m_post) true_param = hpx.vector2Function(m_true, Vh[hpx.PARAMETER]) with dlx.io.XDMFFile( @@ -236,7 +236,7 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: file.write_mesh(msh) file.write_function(true_param) - prior_sample = hpx.vector2Function(m0, Vh[hpx.PARAMETER]) + prior_sample = hpx.vector2Function(m_prior, Vh[hpx.PARAMETER]) with dlx.io.XDMFFile( msh.comm, "dirichlet_Poisson_prior_sample_np{0:d}_X.xdmf".format(nproc), "w" ) as file: diff --git a/example/sfsi_toy_gaussian.py b/example/sfsi_toy_gaussian.py index 3d7a594d..70be97da 100644 --- a/example/sfsi_toy_gaussian.py +++ b/example/sfsi_toy_gaussian.py @@ -225,21 +225,17 @@ def run_inversion( d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1) - lap_aprx = hpx.LaplaceApproximator( - prior, - d, - U, - ) + lap_aprx = hpx.LaplaceApproximator(prior, d, U) lap_aprx.mean = prior.generate_parameter(0) - lap_aprx.mean.array[:] = prior_mean.array[:] + lap_aprx.mean.array[:] = x[hpx.PARAMETER].array[:] noise = prior.generate_parameter("noise") hpx.parRandom.normal(1.0, noise) - m0 = prior.generate_parameter(0) + m_prior = prior.generate_parameter(0) m_post = prior.generate_parameter(0) - lap_aprx.sample(noise, m0, m_post) + lap_aprx.sample(noise, m_prior, m_post) true_param = hpx.vector2Function(m_true, Vh[hpx.PARAMETER]) with dlx.io.XDMFFile( @@ -248,7 +244,7 @@ def run_inversion( file.write_mesh(msh) file.write_function(true_param) - prior_sample = hpx.vector2Function(m0, Vh[hpx.PARAMETER]) + prior_sample = hpx.vector2Function(m_prior, Vh[hpx.PARAMETER]) with dlx.io.XDMFFile( msh.comm, "qpact_prior_sample_np{0:d}_X.xdmf".format(nproc), "w" ) as file: From f23f48dd81851dee21025dc27603220ab2d5c6da Mon Sep 17 00:00:00 2001 From: V-Rang Date: Fri, 12 Apr 2024 12:02:19 -0500 Subject: [PATCH 09/29] test for LowRankHessain added --- example/poisson_dirichlet_example.py | 9 +-- example/poisson_dirichlet_example_reg.py | 2 +- example/poisson_example.py | 2 +- example/poisson_example_reg.py | 2 +- example/sfsi_toy_gaussian.py | 2 +- example/sfsi_toy_gaussian_reg.py | 2 +- hippylibX/test/test_eigendecomposition.py | 2 +- hippylibX/test/test_lowRankHessian.py | 92 +++++++++++++++++++++++ 8 files changed, 99 insertions(+), 14 deletions(-) create mode 100644 hippylibX/test/test_lowRankHessian.py diff --git a/example/poisson_dirichlet_example.py b/example/poisson_dirichlet_example.py index fb8c35e9..3bda7c58 100644 --- a/example/poisson_dirichlet_example.py +++ b/example/poisson_dirichlet_example.py @@ -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 @@ -250,13 +249,7 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: file.write_mesh(msh) file.write_function(posterior_sample) - eigen_decomposition_results = { - "A": Hmisfit.mat, - "B": prior.R, - "k": k, - "d": d, - "U": U, - } + eigen_decomposition_results = {"A": Hmisfit.mat, "B": prior, "k": k, "d": d, "U": U} final_results = { "data_misfit_True": data_misfit_True, diff --git a/example/poisson_dirichlet_example_reg.py b/example/poisson_dirichlet_example_reg.py index a5a641d2..239c9406 100644 --- a/example/poisson_dirichlet_example_reg.py +++ b/example/poisson_dirichlet_example_reg.py @@ -229,7 +229,7 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: eigen_decomposition_results = { "A": Hmisfit.mat, - "B": prior.R, + "B": prior, "k": k, "d": d, "U": U, diff --git a/example/poisson_example.py b/example/poisson_example.py index 18364e43..1545611e 100644 --- a/example/poisson_example.py +++ b/example/poisson_example.py @@ -199,7 +199,7 @@ def run_inversion( eigen_decomposition_results = { "A": Hmisfit.mat, - "B": prior.R, + "B": prior, "k": k, "d": d, "U": U, diff --git a/example/poisson_example_reg.py b/example/poisson_example_reg.py index e17c7a0e..9443ec61 100644 --- a/example/poisson_example_reg.py +++ b/example/poisson_example_reg.py @@ -202,7 +202,7 @@ def run_inversion( eigen_decomposition_results = { "A": Hmisfit.mat, - "B": prior.R, + "B": prior, "k": k, "d": d, "U": U, diff --git a/example/sfsi_toy_gaussian.py b/example/sfsi_toy_gaussian.py index 70be97da..9d9aa137 100644 --- a/example/sfsi_toy_gaussian.py +++ b/example/sfsi_toy_gaussian.py @@ -260,7 +260,7 @@ def run_inversion( eigen_decomposition_results = { "A": Hmisfit.mat, - "B": prior.R, + "B": prior, "k": k, "d": d, "U": U, diff --git a/example/sfsi_toy_gaussian_reg.py b/example/sfsi_toy_gaussian_reg.py index 2f56affa..24b2316f 100644 --- a/example/sfsi_toy_gaussian_reg.py +++ b/example/sfsi_toy_gaussian_reg.py @@ -233,7 +233,7 @@ def run_inversion( eigen_decomposition_results = { "A": Hmisfit.mat, - "B": prior.R, + "B": prior, "k": k, "d": d, "U": U, diff --git a/hippylibX/test/test_eigendecomposition.py b/hippylibX/test/test_eigendecomposition.py index c05dbd4a..f7a1ed5d 100644 --- a/hippylibX/test/test_eigendecomposition.py +++ b/hippylibX/test/test_eigendecomposition.py @@ -83,7 +83,7 @@ def test_eigendecomposition(self): out["eigen_decomposition_results"]["d"], out["eigen_decomposition_results"]["U"], ) - result = check_g(A, B, U, d) + result = check_g(A, B.R, U, d) check_output(self, result) diff --git a/hippylibX/test/test_lowRankHessian.py b/hippylibX/test/test_lowRankHessian.py new file mode 100644 index 00000000..adfe80f5 --- /dev/null +++ b/hippylibX/test/test_lowRankHessian.py @@ -0,0 +1,92 @@ +# Test to check if Hlr^(-1).Hlr == Hlr.Hlr(-1) == Identity +# Hlr -> hpx.LowRankHessian + +# x ,y z as petsc4py Vectors: +# y = Hlr(x) +# z = Hlr^(-1)y => z == x + +# y = Hlr(-1)x +# z = Hlr(y) => z == x + +import unittest +import sys +import os +import numpy as np +import dolfinx as dlx +import dolfinx.fem.petsc +import petsc4py +from typing import Any + +sys.path.append(os.path.abspath("../..")) + +import hippylibX as hpx + +sys.path.append(os.path.abspath("../../example")) +from example import poisson_dirichlet_example + + +def low_Rank_Hessian_mult_solve( + prior: Any, d: np.array, U: hpx.MultiVector +) -> tuple[float, float]: + Hlr = hpx.LowRankHessian(prior, d, U) + vec1 = dlx.la.vector(prior.Vh.dofmap.index_map) + vec2 = dlx.la.vector(prior.Vh.dofmap.index_map) + vec3 = dlx.la.vector(prior.Vh.dofmap.index_map) + + hpx.parRandom.normal(1.0, vec1) + temp_petsc_vec1 = dlx.la.create_petsc_vector_wrap(vec1) + temp_petsc_vec2 = dlx.la.create_petsc_vector_wrap(vec2) + temp_petsc_vec3 = dlx.la.create_petsc_vector_wrap(vec3) + + Hlr.mult(temp_petsc_vec1, temp_petsc_vec2) + Hlr.solve(temp_petsc_vec2, temp_petsc_vec3) + temp_petsc_vec3.axpy(-1.0, temp_petsc_vec1) + + value1 = temp_petsc_vec3.norm(petsc4py.PETSc.NormType.N2) + + hpx.parRandom.normal(1.0, vec1) + + Hlr.solve(temp_petsc_vec1, temp_petsc_vec2) + Hlr.mult(temp_petsc_vec2, temp_petsc_vec3) + temp_petsc_vec3.axpy(-1.0, temp_petsc_vec1) + + value2 = temp_petsc_vec3.norm(petsc4py.PETSc.NormType.N2) + + return value1, value2 + + +def check_output(self, result: tuple[float, float]): + self.assertLessEqual( + np.abs(result[0]), + 1e-6, + "Hlr(-1).Hlr not equal to identity", + ) + + self.assertLessEqual( + np.abs(result[1]), + 1e-6, + "Hlr.Hlr(-1) not equal to identity", + ) + + +class Testing_Execution(unittest.TestCase): + def test_lowRankHessian(self): + hpx.parRandom.replay() + nx = 64 + ny = 64 + noise_variance = 1e-4 + prior_param = {"gamma": 0.1, "delta": 1.0} + out = poisson_dirichlet_example.run_inversion( + nx, ny, noise_variance, prior_param + ) + prior, d, U = ( + out["eigen_decomposition_results"]["B"], + out["eigen_decomposition_results"]["d"], + out["eigen_decomposition_results"]["U"], + ) + result = low_Rank_Hessian_mult_solve(prior, d, U) + check_output(self, result) + + +if __name__ == "__main__": + unittest.main() From e4d6e4600bdbdaa8a39e949777c7559c7e7925b0 Mon Sep 17 00:00:00 2001 From: V-Rang Date: Fri, 12 Apr 2024 12:05:13 -0500 Subject: [PATCH 10/29] test for lowRankHessian --- .github/workflows/CI_testing.yml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.github/workflows/CI_testing.yml b/.github/workflows/CI_testing.yml index 70f69762..b5233cbf 100644 --- a/.github/workflows/CI_testing.yml +++ b/.github/workflows/CI_testing.yml @@ -35,6 +35,11 @@ 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: run serial check run: | cd ./hippylibX/test && From 6a30e9df90b628fbbade3eb557ac00a5af4b1cf6 Mon Sep 17 00:00:00 2001 From: V-Rang Date: Fri, 12 Apr 2024 14:42:28 -0500 Subject: [PATCH 11/29] posterior samples, fix low rank Hessian test --- example/poisson_dirichlet_example.py | 58 +++++++++++++++++++++---- example/sfsi_toy_gaussian.py | 62 ++++++++++++++++++++------- hippylibX/test/test_lowRankHessian.py | 4 ++ 3 files changed, 99 insertions(+), 25 deletions(-) diff --git a/example/poisson_dirichlet_example.py b/example/poisson_dirichlet_example.py index 3bda7c58..a1813831 100644 --- a/example/poisson_dirichlet_example.py +++ b/example/poisson_dirichlet_example.py @@ -220,31 +220,67 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: 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") + + ###################################################### hpx.parRandom.normal(1.0, noise) + lap_aprx.sample(noise, m_prior, m_post) - m_prior = prior.generate_parameter(0) - m_post = prior.generate_parameter(0) + prior_sample = hpx.vector2Function(m_prior, Vh[hpx.PARAMETER]) + with dlx.io.XDMFFile( + msh.comm, "dirichlet_Poisson_prior_sample_1_np{0:d}_X.xdmf".format(nproc), "w" + ) as file: + file.write_mesh(msh) + file.write_function(prior_sample) + posterior_sample = hpx.vector2Function(m_post, Vh[hpx.PARAMETER]) + with dlx.io.XDMFFile( + msh.comm, + "dirichlet_Poisson_posterior_sample_1_np{0:d}_X.xdmf".format(nproc), + "w", + ) as file: + file.write_mesh(msh) + file.write_function(posterior_sample) + + ###################################################### + hpx.parRandom.normal(1.0, noise) lap_aprx.sample(noise, m_prior, m_post) - true_param = hpx.vector2Function(m_true, Vh[hpx.PARAMETER]) + prior_sample = hpx.vector2Function(m_prior, Vh[hpx.PARAMETER]) with dlx.io.XDMFFile( - msh.comm, "dirichlet_Poisson_true_para_np{0:d}_X.xdmf".format(nproc), "w" + msh.comm, "dirichlet_Poisson_prior_sample_2_np{0:d}_X.xdmf".format(nproc), "w" ) as file: file.write_mesh(msh) - file.write_function(true_param) + file.write_function(prior_sample) + + posterior_sample = hpx.vector2Function(m_post, Vh[hpx.PARAMETER]) + with dlx.io.XDMFFile( + msh.comm, + "dirichlet_Poisson_posterior_sample_2_np{0:d}_X.xdmf".format(nproc), + "w", + ) as file: + file.write_mesh(msh) + file.write_function(posterior_sample) + + ###################################################### + hpx.parRandom.normal(1.0, noise) + lap_aprx.sample(noise, m_prior, m_post) prior_sample = hpx.vector2Function(m_prior, Vh[hpx.PARAMETER]) with dlx.io.XDMFFile( - msh.comm, "dirichlet_Poisson_prior_sample_np{0:d}_X.xdmf".format(nproc), "w" + msh.comm, "dirichlet_Poisson_prior_sample_3_np{0:d}_X.xdmf".format(nproc), "w" ) as file: file.write_mesh(msh) file.write_function(prior_sample) posterior_sample = hpx.vector2Function(m_post, Vh[hpx.PARAMETER]) with dlx.io.XDMFFile( - msh.comm, "dirichlet_Poisson_posterior_sample_np{0:d}_X.xdmf".format(nproc), "w" + msh.comm, + "dirichlet_Poisson_posterior_sample_3_np{0:d}_X.xdmf".format(nproc), + "w", ) as file: file.write_mesh(msh) file.write_function(posterior_sample) @@ -265,8 +301,12 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: if __name__ == "__main__": nx = 64 ny = 64 - noise_variance = 1e-4 - prior_param = {"gamma": 0.03, "delta": 0.3} + # noise_variance = 1e-4 + # prior_param = {"gamma": 0.03, "delta": 0.3} + + noise_variance = (1 / 100) * 1e-4 + prior_param = {"gamma": 10 * 0.03, "delta": 10 * 0.3} + final_results = run_inversion(nx, ny, noise_variance, prior_param) k, d = ( final_results["eigen_decomposition_results"]["k"], diff --git a/example/sfsi_toy_gaussian.py b/example/sfsi_toy_gaussian.py index 9d9aa137..122f44c8 100644 --- a/example/sfsi_toy_gaussian.py +++ b/example/sfsi_toy_gaussian.py @@ -229,42 +229,66 @@ def run_inversion( 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") + + ###################################################### hpx.parRandom.normal(1.0, noise) + lap_aprx.sample(noise, m_prior, m_post) - m_prior = prior.generate_parameter(0) - m_post = prior.generate_parameter(0) + prior_sample = hpx.vector2Function(m_prior, Vh[hpx.PARAMETER]) + with dlx.io.XDMFFile( + msh.comm, "pact_prior_sample_1_np{0:d}_X.xdmf".format(nproc), "w" + ) as file: + file.write_mesh(msh) + file.write_function(prior_sample) + + posterior_sample = hpx.vector2Function(m_post, Vh[hpx.PARAMETER]) + with dlx.io.XDMFFile( + msh.comm, "pact_posterior_sample_1_np{0:d}_X.xdmf".format(nproc), "w" + ) as file: + file.write_mesh(msh) + file.write_function(posterior_sample) + ###################################################### + hpx.parRandom.normal(1.0, noise) lap_aprx.sample(noise, m_prior, m_post) - true_param = hpx.vector2Function(m_true, Vh[hpx.PARAMETER]) + prior_sample = hpx.vector2Function(m_prior, Vh[hpx.PARAMETER]) with dlx.io.XDMFFile( - msh.comm, "qpact_true_para_np{0:d}_X.xdmf".format(nproc), "w" + msh.comm, "pact_prior_sample_2_np{0:d}_X.xdmf".format(nproc), "w" ) as file: file.write_mesh(msh) - file.write_function(true_param) + file.write_function(prior_sample) + + posterior_sample = hpx.vector2Function(m_post, Vh[hpx.PARAMETER]) + with dlx.io.XDMFFile( + msh.comm, "pact_posterior_sample_2_np{0:d}_X.xdmf".format(nproc), "w" + ) as file: + file.write_mesh(msh) + file.write_function(posterior_sample) + + ###################################################### + hpx.parRandom.normal(1.0, noise) + lap_aprx.sample(noise, m_prior, m_post) prior_sample = hpx.vector2Function(m_prior, Vh[hpx.PARAMETER]) with dlx.io.XDMFFile( - msh.comm, "qpact_prior_sample_np{0:d}_X.xdmf".format(nproc), "w" + msh.comm, "pact_prior_sample_3_np{0:d}_X.xdmf".format(nproc), "w" ) as file: file.write_mesh(msh) file.write_function(prior_sample) posterior_sample = hpx.vector2Function(m_post, Vh[hpx.PARAMETER]) with dlx.io.XDMFFile( - msh.comm, "qpact_posterior_sample_np{0:d}_X.xdmf".format(nproc), "w" + msh.comm, "pact_posterior_sample_3_np{0:d}_X.xdmf".format(nproc), "w" ) as file: file.write_mesh(msh) file.write_function(posterior_sample) - eigen_decomposition_results = { - "A": Hmisfit.mat, - "B": prior, - "k": k, - "d": d, - "U": U, - } + eigen_decomposition_results = {"A": Hmisfit.mat, "B": prior, "k": k, "d": d, "U": U} final_results = { "data_misfit_True": data_misfit_True, @@ -273,6 +297,8 @@ def run_inversion( "eigen_decomposition_results": eigen_decomposition_results, } + print(comm.rank, ":", misfit.cost(x)) + return final_results ###################################### @@ -280,8 +306,12 @@ def run_inversion( if __name__ == "__main__": nx = 64 ny = 64 - noise_variance = 1e-6 - prior_param = {"gamma": 0.040, "delta": 0.8} + # noise_variance = 1e-6 + # prior_param = {"gamma": 0.040, "delta": 0.8} + + noise_variance = (1 / 100) * 1e-6 + prior_param = {"gamma": 10 * 0.040, "delta": 10 * 0.8} + mesh_filename = "./meshes/circle.xdmf" final_results = run_inversion(mesh_filename, nx, ny, noise_variance, prior_param) k, d = ( diff --git a/hippylibX/test/test_lowRankHessian.py b/hippylibX/test/test_lowRankHessian.py index adfe80f5..47b7c91c 100644 --- a/hippylibX/test/test_lowRankHessian.py +++ b/hippylibX/test/test_lowRankHessian.py @@ -52,6 +52,10 @@ def low_Rank_Hessian_mult_solve( value2 = temp_petsc_vec3.norm(petsc4py.PETSc.NormType.N2) + temp_petsc_vec1.destroy() + temp_petsc_vec2.destroy() + temp_petsc_vec3.destroy() + return value1, value2 From 07b5348810ebdaa167cd907f43d675f6a2fae5fa Mon Sep 17 00:00:00 2001 From: V-Rang Date: Mon, 15 Apr 2024 14:25:39 -0500 Subject: [PATCH 12/29] changes to writing out prior, posterior samples, llaplaceApproximation modification --- example/poisson_dirichlet_example.py | 91 +++++++-------------- example/poisson_dirichlet_example_reg.py | 52 ++++++++++-- example/poisson_example.py | 48 +++++++++-- example/poisson_example_reg.py | 44 ++++++++++ example/sfsi_toy_gaussian.py | 93 ++++++++-------------- example/sfsi_toy_gaussian_reg.py | 50 ++++++++++-- hippylibX/algorithms/multivector.py | 4 + hippylibX/modeling/laplaceApproximation.py | 56 +++++-------- 8 files changed, 259 insertions(+), 179 deletions(-) diff --git a/example/poisson_dirichlet_example.py b/example/poisson_dirichlet_example.py index a1813831..6433a925 100644 --- a/example/poisson_dirichlet_example.py +++ b/example/poisson_dirichlet_example.py @@ -183,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) @@ -216,6 +216,7 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1) + # 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[:] @@ -225,65 +226,37 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: noise = prior.generate_parameter("noise") - ###################################################### - hpx.parRandom.normal(1.0, noise) - lap_aprx.sample(noise, m_prior, m_post) - - prior_sample = hpx.vector2Function(m_prior, Vh[hpx.PARAMETER]) - with dlx.io.XDMFFile( - msh.comm, "dirichlet_Poisson_prior_sample_1_np{0:d}_X.xdmf".format(nproc), "w" - ) as file: - file.write_mesh(msh) - file.write_function(prior_sample) + num_samples_generate = 5 - posterior_sample = hpx.vector2Function(m_post, Vh[hpx.PARAMETER]) - with dlx.io.XDMFFile( - msh.comm, - "dirichlet_Poisson_posterior_sample_1_np{0:d}_X.xdmf".format(nproc), - "w", - ) as file: - file.write_mesh(msh) - file.write_function(posterior_sample) - - ###################################################### - hpx.parRandom.normal(1.0, noise) - lap_aprx.sample(noise, m_prior, m_post) - - prior_sample = hpx.vector2Function(m_prior, Vh[hpx.PARAMETER]) - with dlx.io.XDMFFile( - msh.comm, "dirichlet_Poisson_prior_sample_2_np{0:d}_X.xdmf".format(nproc), "w" - ) as file: - file.write_mesh(msh) - file.write_function(prior_sample) + 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) - posterior_sample = hpx.vector2Function(m_post, Vh[hpx.PARAMETER]) - with dlx.io.XDMFFile( + with dlx.io.VTXWriter( msh.comm, - "dirichlet_Poisson_posterior_sample_2_np{0:d}_X.xdmf".format(nproc), - "w", - ) as file: - file.write_mesh(msh) - file.write_function(posterior_sample) - - ###################################################### - hpx.parRandom.normal(1.0, noise) - lap_aprx.sample(noise, m_prior, m_post) - - prior_sample = hpx.vector2Function(m_prior, Vh[hpx.PARAMETER]) - with dlx.io.XDMFFile( - msh.comm, "dirichlet_Poisson_prior_sample_3_np{0:d}_X.xdmf".format(nproc), "w" - ) as file: - file.write_mesh(msh) - file.write_function(prior_sample) + "poisson_Dirichlet_prior_Bilaplacian_samples_prior_np{0:d}.bp".format(nproc), + prior_samples, + ) as vtx: + vtx.write(0.0) - posterior_sample = hpx.vector2Function(m_post, Vh[hpx.PARAMETER]) - with dlx.io.XDMFFile( + with dlx.io.VTXWriter( msh.comm, - "dirichlet_Poisson_posterior_sample_3_np{0:d}_X.xdmf".format(nproc), - "w", - ) as file: - file.write_mesh(msh) - file.write_function(posterior_sample) + "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.mat, "B": prior, "k": k, "d": d, "U": U} @@ -295,17 +268,13 @@ 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} - - noise_variance = (1 / 100) * 1e-4 - prior_param = {"gamma": 10 * 0.03, "delta": 10 * 0.3} + noise_variance = 1e-4 + prior_param = {"gamma": 0.03, "delta": 0.3} final_results = run_inversion(nx, ny, noise_variance, prior_param) k, d = ( diff --git a/example/poisson_dirichlet_example_reg.py b/example/poisson_dirichlet_example_reg.py index 239c9406..c80bb3d5 100644 --- a/example/poisson_dirichlet_example_reg.py +++ b/example/poisson_dirichlet_example_reg.py @@ -227,13 +227,51 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1) - eigen_decomposition_results = { - "A": Hmisfit.mat, - "B": prior, - "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_Variational_Regularization_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_Variational_Regularization_samples_posterior_np{0:d}.bp".format( + nproc + ), + posterior_samples, + ) as vtx: + vtx.write(0.0) + + eigen_decomposition_results = {"A": Hmisfit.mat, "B": prior, "k": k, "d": d, "U": U} final_results = { "data_misfit_True": data_misfit_True, diff --git a/example/poisson_example.py b/example/poisson_example.py index 1545611e..c3148978 100644 --- a/example/poisson_example.py +++ b/example/poisson_example.py @@ -197,13 +197,47 @@ def run_inversion( d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1) - eigen_decomposition_results = { - "A": Hmisfit.mat, - "B": prior, - "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.mat, "B": prior, "k": k, "d": d, "U": U} final_results = { "data_misfit_True": data_misfit_True, diff --git a/example/poisson_example_reg.py b/example/poisson_example_reg.py index 9443ec61..8599f213 100644 --- a/example/poisson_example_reg.py +++ b/example/poisson_example_reg.py @@ -200,6 +200,50 @@ def run_inversion( d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1) + # 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_Variational_Regularization_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_Variational_Regularization_samples_posterior_np{0:d}.bp".format( + nproc + ), + posterior_samples, + ) as vtx: + vtx.write(0.0) + eigen_decomposition_results = { "A": Hmisfit.mat, "B": prior, diff --git a/example/sfsi_toy_gaussian.py b/example/sfsi_toy_gaussian.py index 122f44c8..1b568b4d 100644 --- a/example/sfsi_toy_gaussian.py +++ b/example/sfsi_toy_gaussian.py @@ -225,6 +225,7 @@ def run_inversion( d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1) + # 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[:] @@ -234,59 +235,35 @@ def run_inversion( noise = prior.generate_parameter("noise") - ###################################################### - hpx.parRandom.normal(1.0, noise) - lap_aprx.sample(noise, m_prior, m_post) - - prior_sample = hpx.vector2Function(m_prior, Vh[hpx.PARAMETER]) - with dlx.io.XDMFFile( - msh.comm, "pact_prior_sample_1_np{0:d}_X.xdmf".format(nproc), "w" - ) as file: - file.write_mesh(msh) - file.write_function(prior_sample) - - posterior_sample = hpx.vector2Function(m_post, Vh[hpx.PARAMETER]) - with dlx.io.XDMFFile( - msh.comm, "pact_posterior_sample_1_np{0:d}_X.xdmf".format(nproc), "w" - ) as file: - file.write_mesh(msh) - file.write_function(posterior_sample) - - ###################################################### - hpx.parRandom.normal(1.0, noise) - lap_aprx.sample(noise, m_prior, m_post) - - prior_sample = hpx.vector2Function(m_prior, Vh[hpx.PARAMETER]) - with dlx.io.XDMFFile( - msh.comm, "pact_prior_sample_2_np{0:d}_X.xdmf".format(nproc), "w" - ) as file: - file.write_mesh(msh) - file.write_function(prior_sample) - - posterior_sample = hpx.vector2Function(m_post, Vh[hpx.PARAMETER]) - with dlx.io.XDMFFile( - msh.comm, "pact_posterior_sample_2_np{0:d}_X.xdmf".format(nproc), "w" - ) as file: - file.write_mesh(msh) - file.write_function(posterior_sample) - - ###################################################### - hpx.parRandom.normal(1.0, noise) - lap_aprx.sample(noise, m_prior, m_post) - - prior_sample = hpx.vector2Function(m_prior, Vh[hpx.PARAMETER]) - with dlx.io.XDMFFile( - msh.comm, "pact_prior_sample_3_np{0:d}_X.xdmf".format(nproc), "w" - ) as file: - file.write_mesh(msh) - file.write_function(prior_sample) - - posterior_sample = hpx.vector2Function(m_post, Vh[hpx.PARAMETER]) - with dlx.io.XDMFFile( - msh.comm, "pact_posterior_sample_3_np{0:d}_X.xdmf".format(nproc), "w" - ) as file: - file.write_mesh(msh) - file.write_function(posterior_sample) + 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.mat, "B": prior, "k": k, "d": d, "U": U} @@ -297,8 +274,6 @@ def run_inversion( "eigen_decomposition_results": eigen_decomposition_results, } - print(comm.rank, ":", misfit.cost(x)) - return final_results ###################################### @@ -306,12 +281,8 @@ def run_inversion( if __name__ == "__main__": nx = 64 ny = 64 - # noise_variance = 1e-6 - # prior_param = {"gamma": 0.040, "delta": 0.8} - - noise_variance = (1 / 100) * 1e-6 - prior_param = {"gamma": 10 * 0.040, "delta": 10 * 0.8} - + noise_variance = 1e-6 + prior_param = {"gamma": 0.040, "delta": 0.8} mesh_filename = "./meshes/circle.xdmf" final_results = run_inversion(mesh_filename, nx, ny, noise_variance, prior_param) k, d = ( diff --git a/example/sfsi_toy_gaussian_reg.py b/example/sfsi_toy_gaussian_reg.py index 24b2316f..35310ebd 100644 --- a/example/sfsi_toy_gaussian_reg.py +++ b/example/sfsi_toy_gaussian_reg.py @@ -231,13 +231,49 @@ def run_inversion( d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1) - eigen_decomposition_results = { - "A": Hmisfit.mat, - "B": prior, - "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_Variational_Regularization_samples_prior_np{0:d}.bp".format(nproc), + prior_samples, + ) as vtx: + vtx.write(0.0) + + with dlx.io.VTXWriter( + msh.comm, + "qpact_prior_Variational_Regularization_samples_posterior_np{0:d}.bp".format( + nproc + ), + posterior_samples, + ) as vtx: + vtx.write(0.0) + + eigen_decomposition_results = {"A": Hmisfit.mat, "B": prior, "k": k, "d": d, "U": U} final_results = { "data_misfit_True": data_misfit_True, diff --git a/hippylibX/algorithms/multivector.py b/hippylibX/algorithms/multivector.py index c01ab402..b732b150 100644 --- a/hippylibX/algorithms/multivector.py +++ b/hippylibX/algorithms/multivector.py @@ -27,6 +27,10 @@ def createFromMultiVec(cls, mv: Type["MultiVector"]) -> Type["MultiVector"]: return mv_copy + def __del__(self) -> None: + for d in self.data: + d.destroy() + def __getitem__(self, k: int) -> petsc4py.PETSc.Vec: return self.data[k] diff --git a/hippylibX/modeling/laplaceApproximation.py b/hippylibX/modeling/laplaceApproximation.py index 2e2504f4..a7f08518 100644 --- a/hippylibX/modeling/laplaceApproximation.py +++ b/hippylibX/modeling/laplaceApproximation.py @@ -13,6 +13,7 @@ # terms of the GNU General Public License (as published by the Free # Software Foundation) version 2.0 dated June 1991. +import petsc4py.PETSc from ..algorithms.lowRankOperator import LowRankOperator from ..algorithms.multivector import MultiVector import numpy as np @@ -33,23 +34,13 @@ class LowRankHessian: of the Hessian and of its inverse. """ - def __init__( - self, - prior, - d: np.array, - U: MultiVector, - createVecLeft=None, - createVecRight=None, - ): + def __init__(self, prior, d: np.array, U: MultiVector): self.U = U self.prior = prior self.LowRankH = LowRankOperator(d, self.U) dsolve = d / (np.ones(d.shape, dtype=d.dtype) + d) self.LowRankHinv = LowRankOperator(dsolve, self.U) - self.createVecLeft = createVecLeft - self.createVecRight = createVecRight - self.help = self.prior.R.createVecRight() self.help1 = self.prior.R.createVecLeft() @@ -57,6 +48,12 @@ def __del__(self) -> None: self.help.destroy() self.help1.destroy() + def createVecRight(self) -> petsc4py.PETSc.Vec: + return self.prior.R.createVecRight() + + def createVecLeft(self) -> petsc4py.PETSc.Vec: + return self.prior.R.createVecLeft() + def mult(self, x: petsc4py.PETSc.Vec, y: petsc4py.PETSc.Vec) -> None: self.prior.R.mult(x, y) self.LowRankH.mult(y, self.help) @@ -86,13 +83,9 @@ def __init__( prior, d: np.array, U: MultiVector, - createVecLeft=None, - createVecRight=None, ): self.U = U self.prior = prior - self.createVecLeft = createVecLeft - self.createVecRight = createVecRight ones = np.ones(d.shape, dtype=d.dtype) self.d = ones - np.power(ones + d, -0.5) @@ -102,14 +95,21 @@ def __init__( def __del__(self) -> None: self.help.destroy() + def createVecRight(self) -> petsc4py.PETSc.Vec: + return self.prior.R.createVecRight() + + def createVecLeft(self) -> petsc4py.PETSc.Vec: + return self.prior.R.createVecLeft() + def sample(self, noise: dlx.la.Vector, s: dlx.la.Vector): temp_petsc_vec_noise = dlx.la.create_petsc_vector_wrap(noise) temp_petsc_vec_s = dlx.la.create_petsc_vector_wrap(s) self.prior.R.mult(temp_petsc_vec_noise, self.help) self.lrsqrt.mult(self.help, temp_petsc_vec_s) - temp_petsc_vec_s.axpy(-1.0, temp_petsc_vec_noise) - temp_petsc_vec_s.scale(-1.0) + # temp_petsc_vec_s.axpy(-1.0, temp_petsc_vec_noise) + # temp_petsc_vec_s.scale(-1.0) + temp_petsc_vec_s.axpby(-1.0, 1.0, temp_petsc_vec_noise) temp_petsc_vec_noise.destroy() temp_petsc_vec_s.destroy() @@ -193,29 +193,13 @@ def sample(self, *args, **kwargs): if len(args) == 2: self._sample_given_prior(args[0], args[1]) if add_mean: - temp_petsc_vec_arg1 = dlx.la.create_petsc_vector_wrap(args[1]) - temp_petsc_vec_mean = dlx.la.create_petsc_vector_wrap(self.mean) - temp_petsc_vec_arg1.axpy(1.0, temp_petsc_vec_mean) - temp_petsc_vec_arg1.destroy() - temp_petsc_vec_mean.destroy() + args[1].array[:] += self.mean.array elif len(args) == 3: self._sample_given_white_noise(args[0], args[1], args[2]) if add_mean: - temp_petsc_vec_arg1 = dlx.la.create_petsc_vector_wrap(args[1]) - temp_petsc_vec_prior_mean = dlx.la.create_petsc_vector_wrap( - self.prior.mean - ) - temp_petsc_vec_arg2 = dlx.la.create_petsc_vector_wrap(args[2]) - temp_petsc_vec_mean = dlx.la.create_petsc_vector_wrap(self.mean) - - temp_petsc_vec_arg1.axpy(1.0, temp_petsc_vec_prior_mean) - temp_petsc_vec_arg2.axpy(1.0, temp_petsc_vec_mean) - - temp_petsc_vec_arg1.destroy() - temp_petsc_vec_prior_mean.destroy() - temp_petsc_vec_arg2.destroy() - temp_petsc_vec_mean.destroy() + args[1].array[:] += self.prior.mean.array + args[2].array[:] += self.mean.array else: raise NameError("Invalid number of parameters in Posterior::sample") From 0a4ec14239030cf99c4b6471ad5620f3d8bf9dab Mon Sep 17 00:00:00 2001 From: V-Rang Date: Mon, 15 Apr 2024 14:30:02 -0500 Subject: [PATCH 13/29] removing commented code --- hippylibX/modeling/laplaceApproximation.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/hippylibX/modeling/laplaceApproximation.py b/hippylibX/modeling/laplaceApproximation.py index a7f08518..69fde63e 100644 --- a/hippylibX/modeling/laplaceApproximation.py +++ b/hippylibX/modeling/laplaceApproximation.py @@ -107,8 +107,6 @@ def sample(self, noise: dlx.la.Vector, s: dlx.la.Vector): self.prior.R.mult(temp_petsc_vec_noise, self.help) self.lrsqrt.mult(self.help, temp_petsc_vec_s) - # temp_petsc_vec_s.axpy(-1.0, temp_petsc_vec_noise) - # temp_petsc_vec_s.scale(-1.0) temp_petsc_vec_s.axpby(-1.0, 1.0, temp_petsc_vec_noise) temp_petsc_vec_noise.destroy() temp_petsc_vec_s.destroy() From 847b2f7858aff50c1cc9d57482698dec8a08dff3 Mon Sep 17 00:00:00 2001 From: V-Rang Date: Mon, 15 Apr 2024 14:39:24 -0500 Subject: [PATCH 14/29] removing code for Variational Regularization prior --- example/poisson_dirichlet_example_reg.py | 45 +----------------------- example/poisson_example_reg.py | 44 ----------------------- example/sfsi_toy_gaussian_reg.py | 42 ---------------------- 3 files changed, 1 insertion(+), 130 deletions(-) diff --git a/example/poisson_dirichlet_example_reg.py b/example/poisson_dirichlet_example_reg.py index c80bb3d5..df17ba21 100644 --- a/example/poisson_dirichlet_example_reg.py +++ b/example/poisson_dirichlet_example_reg.py @@ -226,50 +226,7 @@ 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) - - # 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_Variational_Regularization_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_Variational_Regularization_samples_posterior_np{0:d}.bp".format( - nproc - ), - posterior_samples, - ) as vtx: - vtx.write(0.0) + eigen_decomposition_results = {"A": Hmisfit.mat, "B": prior, "k": k, "d": d, "U": U} diff --git a/example/poisson_example_reg.py b/example/poisson_example_reg.py index 8599f213..9443ec61 100644 --- a/example/poisson_example_reg.py +++ b/example/poisson_example_reg.py @@ -200,50 +200,6 @@ def run_inversion( d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1) - # 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_Variational_Regularization_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_Variational_Regularization_samples_posterior_np{0:d}.bp".format( - nproc - ), - posterior_samples, - ) as vtx: - vtx.write(0.0) - eigen_decomposition_results = { "A": Hmisfit.mat, "B": prior, diff --git a/example/sfsi_toy_gaussian_reg.py b/example/sfsi_toy_gaussian_reg.py index 35310ebd..8f7c8dce 100644 --- a/example/sfsi_toy_gaussian_reg.py +++ b/example/sfsi_toy_gaussian_reg.py @@ -231,48 +231,6 @@ def run_inversion( d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1) - # 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_Variational_Regularization_samples_prior_np{0:d}.bp".format(nproc), - prior_samples, - ) as vtx: - vtx.write(0.0) - - with dlx.io.VTXWriter( - msh.comm, - "qpact_prior_Variational_Regularization_samples_posterior_np{0:d}.bp".format( - nproc - ), - posterior_samples, - ) as vtx: - vtx.write(0.0) - eigen_decomposition_results = {"A": Hmisfit.mat, "B": prior, "k": k, "d": d, "U": U} final_results = { From b096580cea5348b966a8ddc9e027c1ee6d584426 Mon Sep 17 00:00:00 2001 From: V-Rang Date: Mon, 15 Apr 2024 14:42:34 -0500 Subject: [PATCH 15/29] linter fix --- example/poisson_dirichlet_example_reg.py | 1 - 1 file changed, 1 deletion(-) diff --git a/example/poisson_dirichlet_example_reg.py b/example/poisson_dirichlet_example_reg.py index df17ba21..aa6b7201 100644 --- a/example/poisson_dirichlet_example_reg.py +++ b/example/poisson_dirichlet_example_reg.py @@ -226,7 +226,6 @@ 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) - eigen_decomposition_results = {"A": Hmisfit.mat, "B": prior, "k": k, "d": d, "U": U} From bcfe06fb5f4cedefa67edeab4f2fae7952e66e4e Mon Sep 17 00:00:00 2001 From: V-Rang Date: Tue, 16 Apr 2024 13:58:51 -0500 Subject: [PATCH 16/29] fixing solve method for cgsolver --- example/poisson_dirichlet_example.py | 35 ++++++++++++++++++++++++ hippylibX/algorithms/NewtonCG.py | 5 ++-- hippylibX/algorithms/cgsolverSteihaug.py | 6 ++-- 3 files changed, 41 insertions(+), 5 deletions(-) diff --git a/example/poisson_dirichlet_example.py b/example/poisson_dirichlet_example.py index 6433a925..91798f0b 100644 --- a/example/poisson_dirichlet_example.py +++ b/example/poisson_dirichlet_example.py @@ -267,6 +267,40 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: "eigen_decomposition_results": eigen_decomposition_results, } + # Hlr = hpx.LowRankHessian(prior, d, U) + # matr = final_results['eigen_decomposition_results']['A'] + # vec1 = dlx.la.vector(prior.Vh.dofmap.index_map) + # hpx.parRandom.normal(1.0, vec1) + + # vec1 = dlx.fem.Function(Vh_m) + # vec1.interpolate( + # lambda x: np.log(2 + 7 * (((x[0] - 0.5) ** 2 + (x[1] - 0.5) ** 2) ** 0.5 > 0.2)) + # ) + # vec1.x.scatter_forward() + # vec1 = vec1.x + + # model.evalGradientParameter(x,vec1) + # vec1.array[:] *= -1. + + # print(vec1.array.min(),":",vec1.array.max()) + + # need to manufacture a vec1 that will not blow up the solver + # vec1 = dlx.la.vector(prior.Vh.dofmap.index_map) + # hpx.parRandom.normal(1.0, vec1) + + # vec2 = dlx.la.vector(prior.Vh.dofmap.index_map) + + # parameters = hpx.algorithms.cgsolverSteihaug.CGSolverSteihaug_ParameterList() + # solver = hpx.algorithms.cgsolverSteihaug.CGSolverSteihaug(parameters, comm) + # solver.set_operator(matr) + + # # solver.set_preconditioner(prior.Rsolver) + # solver.set_preconditioner(Hlr) + + # solver.solve(vec1, vec2) + # print(vec2.array.min(),":",vec2.array.max()) + # print(solver.iter) + return final_results @@ -277,6 +311,7 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: 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"], diff --git a/hippylibX/algorithms/NewtonCG.py b/hippylibX/algorithms/NewtonCG.py index be345a55..1406bde8 100644 --- a/hippylibX/algorithms/NewtonCG.py +++ b/hippylibX/algorithms/NewtonCG.py @@ -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 @@ -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) diff --git a/hippylibX/algorithms/cgsolverSteihaug.py b/hippylibX/algorithms/cgsolverSteihaug.py index 10720413..d378a50d 100644 --- a/hippylibX/algorithms/cgsolverSteihaug.py +++ b/hippylibX/algorithms/cgsolverSteihaug.py @@ -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` """ From fb0cf34784046f96cf15d56532a87f979b0920e4 Mon Sep 17 00:00:00 2001 From: V-Rang Date: Tue, 16 Apr 2024 16:50:51 -0500 Subject: [PATCH 17/29] adding tests for sampling and preconditioner --- .github/workflows/CI_testing.yml | 10 +++ example/poisson_dirichlet_example.py | 36 +------- example/poisson_dirichlet_example_reg.py | 2 +- example/poisson_example.py | 2 +- example/poisson_example_reg.py | 2 +- example/sfsi_toy_gaussian.py | 2 +- example/sfsi_toy_gaussian_reg.py | 2 +- hippylibX/test/test_eigendecomposition.py | 2 +- .../test_lowRankHessian_preconditioner.py | 64 ++++++++++++++ hippylibX/test/test_sampling.py | 83 +++++++++++++++++++ 10 files changed, 164 insertions(+), 41 deletions(-) create mode 100644 hippylibX/test/test_lowRankHessian_preconditioner.py create mode 100644 hippylibX/test/test_sampling.py diff --git a/.github/workflows/CI_testing.yml b/.github/workflows/CI_testing.yml index b5233cbf..718bb6b7 100644 --- a/.github/workflows/CI_testing.yml +++ b/.github/workflows/CI_testing.yml @@ -40,6 +40,16 @@ jobs: cd ./hippylibX/test && mpirun -n 2 python3 test_lowRankHessian.py + - name: low rank Hessian precondtioner testing + 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 && diff --git a/example/poisson_dirichlet_example.py b/example/poisson_dirichlet_example.py index 91798f0b..92e64618 100644 --- a/example/poisson_dirichlet_example.py +++ b/example/poisson_dirichlet_example.py @@ -258,7 +258,7 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: ) as vtx: vtx.write(0.0) - eigen_decomposition_results = {"A": Hmisfit.mat, "B": prior, "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, @@ -267,40 +267,6 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: "eigen_decomposition_results": eigen_decomposition_results, } - # Hlr = hpx.LowRankHessian(prior, d, U) - # matr = final_results['eigen_decomposition_results']['A'] - # vec1 = dlx.la.vector(prior.Vh.dofmap.index_map) - # hpx.parRandom.normal(1.0, vec1) - - # vec1 = dlx.fem.Function(Vh_m) - # vec1.interpolate( - # lambda x: np.log(2 + 7 * (((x[0] - 0.5) ** 2 + (x[1] - 0.5) ** 2) ** 0.5 > 0.2)) - # ) - # vec1.x.scatter_forward() - # vec1 = vec1.x - - # model.evalGradientParameter(x,vec1) - # vec1.array[:] *= -1. - - # print(vec1.array.min(),":",vec1.array.max()) - - # need to manufacture a vec1 that will not blow up the solver - # vec1 = dlx.la.vector(prior.Vh.dofmap.index_map) - # hpx.parRandom.normal(1.0, vec1) - - # vec2 = dlx.la.vector(prior.Vh.dofmap.index_map) - - # parameters = hpx.algorithms.cgsolverSteihaug.CGSolverSteihaug_ParameterList() - # solver = hpx.algorithms.cgsolverSteihaug.CGSolverSteihaug(parameters, comm) - # solver.set_operator(matr) - - # # solver.set_preconditioner(prior.Rsolver) - # solver.set_preconditioner(Hlr) - - # solver.solve(vec1, vec2) - # print(vec2.array.min(),":",vec2.array.max()) - # print(solver.iter) - return final_results diff --git a/example/poisson_dirichlet_example_reg.py b/example/poisson_dirichlet_example_reg.py index aa6b7201..2b9e927d 100644 --- a/example/poisson_dirichlet_example_reg.py +++ b/example/poisson_dirichlet_example_reg.py @@ -227,7 +227,7 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1) - eigen_decomposition_results = {"A": Hmisfit.mat, "B": prior, "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, diff --git a/example/poisson_example.py b/example/poisson_example.py index c3148978..0ef0e7b2 100644 --- a/example/poisson_example.py +++ b/example/poisson_example.py @@ -237,7 +237,7 @@ def run_inversion( ) as vtx: vtx.write(0.0) - eigen_decomposition_results = {"A": Hmisfit.mat, "B": prior, "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, diff --git a/example/poisson_example_reg.py b/example/poisson_example_reg.py index 9443ec61..52b3d004 100644 --- a/example/poisson_example_reg.py +++ b/example/poisson_example_reg.py @@ -198,7 +198,7 @@ def run_inversion( hpx.parRandom.normal(1.0, Omega) - d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1) + d, U = hpx.doublePassG(Hmisfit, prior.R, prior.Rsolver, Omega, k, s=1) eigen_decomposition_results = { "A": Hmisfit.mat, diff --git a/example/sfsi_toy_gaussian.py b/example/sfsi_toy_gaussian.py index 1b568b4d..7d52dbb2 100644 --- a/example/sfsi_toy_gaussian.py +++ b/example/sfsi_toy_gaussian.py @@ -223,7 +223,7 @@ def run_inversion( hpx.parRandom.normal(1.0, Omega) - d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1) + d, U = hpx.doublePassG(Hmisfit, prior.R, prior.Rsolver, Omega, k, s=1) # generating prior and posterior samples lap_aprx = hpx.LaplaceApproximator(prior, d, U) diff --git a/example/sfsi_toy_gaussian_reg.py b/example/sfsi_toy_gaussian_reg.py index 8f7c8dce..9a74503a 100644 --- a/example/sfsi_toy_gaussian_reg.py +++ b/example/sfsi_toy_gaussian_reg.py @@ -231,7 +231,7 @@ def run_inversion( d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1) - eigen_decomposition_results = {"A": Hmisfit.mat, "B": prior, "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, diff --git a/hippylibX/test/test_eigendecomposition.py b/hippylibX/test/test_eigendecomposition.py index f7a1ed5d..be7553c5 100644 --- a/hippylibX/test/test_eigendecomposition.py +++ b/hippylibX/test/test_eigendecomposition.py @@ -83,7 +83,7 @@ def test_eigendecomposition(self): out["eigen_decomposition_results"]["d"], out["eigen_decomposition_results"]["U"], ) - result = check_g(A, B.R, U, d) + result = check_g(A.mat, B.R, U, d) check_output(self, result) diff --git a/hippylibX/test/test_lowRankHessian_preconditioner.py b/hippylibX/test/test_lowRankHessian_preconditioner.py new file mode 100644 index 00000000..2974ca1c --- /dev/null +++ b/hippylibX/test/test_lowRankHessian_preconditioner.py @@ -0,0 +1,64 @@ +# Low rank Hessian should be a very good preconditioner for the data +# misfit Hessian. + +import unittest +import sys +import os +import dolfinx as dlx + +sys.path.append(os.path.abspath("../..")) + +import hippylibX as hpx + +sys.path.append(os.path.abspath("../../example")) +from example import poisson_dirichlet_example + + +class Testing_Execution(unittest.TestCase): + def test_lowRankHessian_precondtioner(self): + hpx.parRandom.replay() + nx = 64 + ny = 64 + noise_variance = 1e-4 + prior_param = {"gamma": 0.1, "delta": 1.0} + out = poisson_dirichlet_example.run_inversion( + nx, ny, noise_variance, prior_param + ) + Hmisfit, prior, d, U = ( + out["eigen_decomposition_results"]["A"], + out["eigen_decomposition_results"]["B"], + out["eigen_decomposition_results"]["d"], + out["eigen_decomposition_results"]["U"], + ) + lowRank_Hessian = hpx.LowRankHessian(prior, d, U) + + Hmisfit.misfit_only = False + matr = Hmisfit.mat + + vec1 = dlx.la.vector(prior.Vh.dofmap.index_map) + hpx.parRandom.normal(1.0, vec1) + + vec2 = dlx.la.vector(prior.Vh.dofmap.index_map) + parameters = hpx.algorithms.cgsolverSteihaug.CGSolverSteihaug_ParameterList() + solver = hpx.algorithms.cgsolverSteihaug.CGSolverSteihaug( + parameters, prior.Vh.mesh.comm + ) + solver.set_operator(matr) + + solver.set_preconditioner(prior.Rsolver) + solver.solve(vec1, vec2) + Rsolver_iter_count = solver.iter + + solver.set_preconditioner(lowRank_Hessian) + solver.solve(vec1, vec2) + Hlr_iter_count = solver.iter + + self.assertLess( + Hlr_iter_count, + Rsolver_iter_count, + "lowRank Hessian not better pc than prior.Rsolver", + ) + + +if __name__ == "__main__": + unittest.main() diff --git a/hippylibX/test/test_sampling.py b/hippylibX/test/test_sampling.py new file mode 100644 index 00000000..3a5da8d8 --- /dev/null +++ b/hippylibX/test/test_sampling.py @@ -0,0 +1,83 @@ +# testing the lowRankHessian.mult method from modeling/laplaceApproximation.py +# x, y petsc vectors +# h1 = Hlr(x) +# h2 = Hlr(y) +# result1 = h2.h1 (inner product) + +# h1 = Hlr(x) +# h2 = Hlr(h1) +# result2 = y.dot(h2) (inner product) + +# result1 == result2 + +import unittest +import sys +import os +import numpy as np +import dolfinx as dlx + +sys.path.append(os.path.abspath("../..")) + +import hippylibX as hpx + +sys.path.append(os.path.abspath("../../example")) +from example import poisson_dirichlet_example + + +class Testing_Execution(unittest.TestCase): + def test_sampling(self): + hpx.parRandom.replay() + + nx = 64 + ny = 64 + noise_variance = 1e-4 + prior_param = {"gamma": 0.1, "delta": 1.0} + out = poisson_dirichlet_example.run_inversion( + nx, ny, noise_variance, prior_param + ) + prior, d, U = ( + out["eigen_decomposition_results"]["B"], + out["eigen_decomposition_results"]["d"], + out["eigen_decomposition_results"]["U"], + ) + + Hlr = hpx.LowRankHessian(prior, d, U) + vec1 = dlx.la.vector(prior.Vh.dofmap.index_map) + hpx.parRandom.normal(1.0, vec1) + + vec2 = dlx.la.vector(prior.Vh.dofmap.index_map) + hpx.parRandom.normal(1.0, vec2) + + temp_petsc_vec1 = dlx.la.create_petsc_vector_wrap(vec1) + temp_petsc_vec2 = dlx.la.create_petsc_vector_wrap(vec2) + + help1 = dlx.la.create_petsc_vector( + prior.Vh.dofmap.index_map, prior.Vh.dofmap.bs + ) + help2 = dlx.la.create_petsc_vector( + prior.Vh.dofmap.index_map, prior.Vh.dofmap.bs + ) + + Hlr.mult(temp_petsc_vec1, help1) + Hlr.mult(temp_petsc_vec2, help2) + + result_1 = help2.dot(help1) + + Hlr.mult(temp_petsc_vec1, help1) + Hlr.mult(help1, help2) + result_2 = temp_petsc_vec2.dot(help2) + + temp_petsc_vec1.destroy() + temp_petsc_vec2.destroy() + help1.destroy() + help2.destroy() + + self.assertLessEqual( + np.abs(result_1 - result_2), + 1e-3, + "lowRankHessian.mult failed", + ) + + +if __name__ == "__main__": + unittest.main() From ee5edbb51351a5d67cb89db389dde343ef0dc55f Mon Sep 17 00:00:00 2001 From: V-Rang Date: Tue, 16 Apr 2024 17:09:36 -0500 Subject: [PATCH 18/29] misfit argument fix --- example/poisson_example_reg.py | 2 +- example/sfsi_toy_gaussian.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/example/poisson_example_reg.py b/example/poisson_example_reg.py index 52b3d004..9443ec61 100644 --- a/example/poisson_example_reg.py +++ b/example/poisson_example_reg.py @@ -198,7 +198,7 @@ def run_inversion( hpx.parRandom.normal(1.0, Omega) - d, U = hpx.doublePassG(Hmisfit, prior.R, prior.Rsolver, Omega, k, s=1) + d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1) eigen_decomposition_results = { "A": Hmisfit.mat, diff --git a/example/sfsi_toy_gaussian.py b/example/sfsi_toy_gaussian.py index 7d52dbb2..1b568b4d 100644 --- a/example/sfsi_toy_gaussian.py +++ b/example/sfsi_toy_gaussian.py @@ -223,7 +223,7 @@ def run_inversion( hpx.parRandom.normal(1.0, Omega) - d, U = hpx.doublePassG(Hmisfit, prior.R, prior.Rsolver, Omega, k, s=1) + d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1) # generating prior and posterior samples lap_aprx = hpx.LaplaceApproximator(prior, d, U) From 085d5bfe1ec861db7c2af060a83eb0bcdd52fbab Mon Sep 17 00:00:00 2001 From: V-Rang Date: Tue, 16 Apr 2024 17:15:26 -0500 Subject: [PATCH 19/29] return type changed --- example/poisson_example_reg.py | 8 +------- example/sfsi_toy_gaussian.py | 2 +- 2 files changed, 2 insertions(+), 8 deletions(-) diff --git a/example/poisson_example_reg.py b/example/poisson_example_reg.py index 9443ec61..6b8034aa 100644 --- a/example/poisson_example_reg.py +++ b/example/poisson_example_reg.py @@ -200,13 +200,7 @@ def run_inversion( d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1) - eigen_decomposition_results = { - "A": Hmisfit.mat, - "B": prior, - "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, diff --git a/example/sfsi_toy_gaussian.py b/example/sfsi_toy_gaussian.py index 1b568b4d..5217507e 100644 --- a/example/sfsi_toy_gaussian.py +++ b/example/sfsi_toy_gaussian.py @@ -265,7 +265,7 @@ def run_inversion( ) as vtx: vtx.write(0.0) - eigen_decomposition_results = {"A": Hmisfit.mat, "B": prior, "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, From a8ecee490ea5270570b0aba43a22af152fa52600 Mon Sep 17 00:00:00 2001 From: V-Rang Date: Tue, 16 Apr 2024 18:42:25 -0500 Subject: [PATCH 20/29] sampler check not working --- hippylibX/test/test_sampling.py | 58 ++++++++++++++++----------------- 1 file changed, 29 insertions(+), 29 deletions(-) diff --git a/hippylibX/test/test_sampling.py b/hippylibX/test/test_sampling.py index 3a5da8d8..d4456cce 100644 --- a/hippylibX/test/test_sampling.py +++ b/hippylibX/test/test_sampling.py @@ -13,7 +13,6 @@ import unittest import sys import os -import numpy as np import dolfinx as dlx sys.path.append(os.path.abspath("../..")) @@ -41,42 +40,43 @@ def test_sampling(self): out["eigen_decomposition_results"]["U"], ) - Hlr = hpx.LowRankHessian(prior, d, U) - vec1 = dlx.la.vector(prior.Vh.dofmap.index_map) - hpx.parRandom.normal(1.0, vec1) - - vec2 = dlx.la.vector(prior.Vh.dofmap.index_map) - hpx.parRandom.normal(1.0, vec2) + x = dlx.la.vector(prior.Vh.dofmap.index_map) + hpx.parRandom.normal(1.0, x) - temp_petsc_vec1 = dlx.la.create_petsc_vector_wrap(vec1) - temp_petsc_vec2 = dlx.la.create_petsc_vector_wrap(vec2) + y = dlx.la.vector(prior.Vh.dofmap.index_map) + hpx.parRandom.normal(1.0, y) - help1 = dlx.la.create_petsc_vector( - prior.Vh.dofmap.index_map, prior.Vh.dofmap.bs - ) - help2 = dlx.la.create_petsc_vector( - prior.Vh.dofmap.index_map, prior.Vh.dofmap.bs + sx, sy = ( + dlx.la.vector(prior.Vh.dofmap.index_map), + dlx.la.vector(prior.Vh.dofmap.index_map), ) - Hlr.mult(temp_petsc_vec1, help1) - Hlr.mult(temp_petsc_vec2, help2) + PS_lr = hpx.LowRankPosteriorSampler(prior, d, U) + PS_lr.sample(x, sx) + PS_lr.sample(y, sy) + + result_1 = dlx.cpp.la.inner_product(sy._cpp_object, sx._cpp_object) - result_1 = help2.dot(help1) + Hlr = hpx.LowRankHessian(prior, d, U) + Hx = dlx.la.create_petsc_vector(prior.Vh.dofmap.index_map, prior.Vh.dofmap.bs) - Hlr.mult(temp_petsc_vec1, help1) - Hlr.mult(help1, help2) - result_2 = temp_petsc_vec2.dot(help2) + temp_petsc_vec_x = dlx.la.create_petsc_vector_wrap(x) + temp_petsc_vec_y = dlx.la.create_petsc_vector_wrap(y) - temp_petsc_vec1.destroy() - temp_petsc_vec2.destroy() - help1.destroy() - help2.destroy() + Hlr.mult(temp_petsc_vec_x, Hx) + result_2 = temp_petsc_vec_y.dot(Hx) - self.assertLessEqual( - np.abs(result_1 - result_2), - 1e-3, - "lowRankHessian.mult failed", - ) + temp_petsc_vec_x.destroy() + temp_petsc_vec_y.destroy() + Hx.destroy() + + print(result_1, result_2) + + # self.assertLessEqual( + # np.abs(result_1 - result_2), + # 1e-3, + # "lowRankHessian.mult failed", + # ) if __name__ == "__main__": From b44b0ae499fdc01ac0190742fecccc52f979a42f Mon Sep 17 00:00:00 2001 From: V-Rang Date: Wed, 17 Apr 2024 11:23:09 -0500 Subject: [PATCH 21/29] checking sample test --- hippylibX/test/test_sampling.py | 67 ++++++++++++++++----------------- 1 file changed, 33 insertions(+), 34 deletions(-) diff --git a/hippylibX/test/test_sampling.py b/hippylibX/test/test_sampling.py index d4456cce..273f88ae 100644 --- a/hippylibX/test/test_sampling.py +++ b/hippylibX/test/test_sampling.py @@ -1,19 +1,18 @@ -# testing the lowRankHessian.mult method from modeling/laplaceApproximation.py -# x, y petsc vectors -# h1 = Hlr(x) -# h2 = Hlr(y) -# result1 = h2.h1 (inner product) - -# h1 = Hlr(x) -# h2 = Hlr(h1) -# result2 = y.dot(h2) (inner product) - -# result1 == result2 +# testing the lowRankPosteriorSampler.sample method from modeling/laplaceApproximation.py +# x : random vector +# prior.Rsolver.solve(x, x1) +# Hlrsampler.sample(x1, x2) +# Hlrsampler.sample(x2, x3) +# Hlr.solve(x, y) +# assert x3 == y import unittest import sys import os import dolfinx as dlx +import petsc4py +import dolfinx.fem.petsc +import numpy as np sys.path.append(os.path.abspath("../..")) @@ -41,42 +40,42 @@ def test_sampling(self): ) x = dlx.la.vector(prior.Vh.dofmap.index_map) - hpx.parRandom.normal(1.0, x) - + x1 = dlx.la.vector(prior.Vh.dofmap.index_map) + x2 = dlx.la.vector(prior.Vh.dofmap.index_map) + x3 = dlx.la.vector(prior.Vh.dofmap.index_map) y = dlx.la.vector(prior.Vh.dofmap.index_map) - hpx.parRandom.normal(1.0, y) - sx, sy = ( - dlx.la.vector(prior.Vh.dofmap.index_map), - dlx.la.vector(prior.Vh.dofmap.index_map), - ) + hpx.parRandom.normal(1.0, x) - PS_lr = hpx.LowRankPosteriorSampler(prior, d, U) - PS_lr.sample(x, sx) - PS_lr.sample(y, sy) + temp_petsc_vec_x = dlx.la.create_petsc_vector_wrap(x) + temp_petsc_vec_x1 = dlx.la.create_petsc_vector_wrap(x1) + prior.Rsolver.solve(temp_petsc_vec_x, temp_petsc_vec_x1) - result_1 = dlx.cpp.la.inner_product(sy._cpp_object, sx._cpp_object) + PS_lr = hpx.LowRankPosteriorSampler(prior, d, U) + PS_lr.sample(x1, x2) + PS_lr.sample(x2, x3) Hlr = hpx.LowRankHessian(prior, d, U) - Hx = dlx.la.create_petsc_vector(prior.Vh.dofmap.index_map, prior.Vh.dofmap.bs) - temp_petsc_vec_x = dlx.la.create_petsc_vector_wrap(x) temp_petsc_vec_y = dlx.la.create_petsc_vector_wrap(y) + Hlr.solve(temp_petsc_vec_x, temp_petsc_vec_y) - Hlr.mult(temp_petsc_vec_x, Hx) - result_2 = temp_petsc_vec_y.dot(Hx) + # assert x3 == y + temp_petsc_vec_x3 = dlx.la.create_petsc_vector_wrap(x3) + temp_petsc_vec_y.axpy(-1.0, temp_petsc_vec_x3) + + value = temp_petsc_vec_y.norm(petsc4py.PETSc.NormType.N2) temp_petsc_vec_x.destroy() + temp_petsc_vec_x1.destroy() + temp_petsc_vec_x3.destroy() temp_petsc_vec_y.destroy() - Hx.destroy() - - print(result_1, result_2) - # self.assertLessEqual( - # np.abs(result_1 - result_2), - # 1e-3, - # "lowRankHessian.mult failed", - # ) + self.assertLessEqual( + np.abs(value), + 1e-6, + "lowRankPosteriorSampler.sample failed", + ) if __name__ == "__main__": From ecf7c855eeeaf3b724c90663e6f45af3486ee9a8 Mon Sep 17 00:00:00 2001 From: V-Rang Date: Wed, 17 Apr 2024 11:33:58 -0500 Subject: [PATCH 22/29] rename sample to mult --- hippylibX/modeling/laplaceApproximation.py | 6 +++--- hippylibX/test/test_sampling.py | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/hippylibX/modeling/laplaceApproximation.py b/hippylibX/modeling/laplaceApproximation.py index 69fde63e..5e999c43 100644 --- a/hippylibX/modeling/laplaceApproximation.py +++ b/hippylibX/modeling/laplaceApproximation.py @@ -101,7 +101,7 @@ def createVecRight(self) -> petsc4py.PETSc.Vec: def createVecLeft(self) -> petsc4py.PETSc.Vec: return self.prior.R.createVecLeft() - def sample(self, noise: dlx.la.Vector, s: dlx.la.Vector): + def mult(self, noise: dlx.la.Vector, s: dlx.la.Vector): temp_petsc_vec_noise = dlx.la.create_petsc_vector_wrap(noise) temp_petsc_vec_s = dlx.la.create_petsc_vector_wrap(s) @@ -205,10 +205,10 @@ def _sample_given_white_noise( self, noise: dlx.la.Vector, s_prior: dlx.la.Vector, s_post: dlx.la.Vector ): self.prior.sample(noise, s_prior, add_mean=False) - self.sampler.sample(s_prior, s_post) + self.sampler.mult(s_prior, s_post) def _sample_given_prior(self, s_prior: dlx.la.Vector, s_post: dlx.la.Vector): - self.sampler.sample(s_prior, s_post) + self.sampler.mult(s_prior, s_post) @not_implemented def trace(self, **kwargs): diff --git a/hippylibX/test/test_sampling.py b/hippylibX/test/test_sampling.py index 273f88ae..fdf427f7 100644 --- a/hippylibX/test/test_sampling.py +++ b/hippylibX/test/test_sampling.py @@ -52,8 +52,8 @@ def test_sampling(self): prior.Rsolver.solve(temp_petsc_vec_x, temp_petsc_vec_x1) PS_lr = hpx.LowRankPosteriorSampler(prior, d, U) - PS_lr.sample(x1, x2) - PS_lr.sample(x2, x3) + PS_lr.mult(x1, x2) + PS_lr.mult(x2, x3) Hlr = hpx.LowRankHessian(prior, d, U) From 81d7a9af3c5891916bf56c5d8eda498bf2af81e3 Mon Sep 17 00:00:00 2001 From: V-Rang Date: Wed, 17 Apr 2024 13:57:05 -0500 Subject: [PATCH 23/29] time series function plot not working --- example/poisson_dirichlet_example.py | 44 +++++++++++++++++++++------- 1 file changed, 33 insertions(+), 11 deletions(-) diff --git a/example/poisson_dirichlet_example.py b/example/poisson_dirichlet_example.py index 92e64618..d55ddea7 100644 --- a/example/poisson_dirichlet_example.py +++ b/example/poisson_dirichlet_example.py @@ -242,21 +242,43 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: prior_samples.append(prior_sample) posterior_samples.append(posterior_sample) - with dlx.io.VTXWriter( + with dlx.io.XDMFFile( 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( + "poisson_Dirichlet_prior_Bilaplacian_samples_prior_np{0:d}.xdmf".format(nproc), + "w", + ) as file: + file.write_mesh(msh) + for i, sample in enumerate(prior_samples): + time_value = i / (num_samples_generate - 1) + file.write_function(sample, time_value) + + with dlx.io.XDMFFile( msh.comm, - "poisson_Dirichlet_prior_Bilaplacian_samples_posterior_np{0:d}.bp".format( + "poisson_Dirichlet_prior_Bilaplacian_samples_posterior_np{0:d}.xdmf".format( nproc ), - posterior_samples, - ) as vtx: - vtx.write(0.0) + "w", + ) as file: + file.write_mesh(msh) + for i, sample in enumerate(posterior_samples): + time_value = i / (num_samples_generate - 1) + file.write_function(sample, time_value) + + # 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} From 3ff241ae40cf4e0459885c9685216f5052bbf8a8 Mon Sep 17 00:00:00 2001 From: V-Rang Date: Thu, 18 Apr 2024 21:01:19 -0500 Subject: [PATCH 24/29] changing noise level, sampler --- example/poisson_dirichlet_example.py | 46 +++++++----------------- example/poisson_dirichlet_example_reg.py | 2 +- example/poisson_example.py | 5 +-- example/poisson_example_reg.py | 4 +-- 4 files changed, 18 insertions(+), 39 deletions(-) diff --git a/example/poisson_dirichlet_example.py b/example/poisson_dirichlet_example.py index d55ddea7..8e70bae7 100644 --- a/example/poisson_dirichlet_example.py +++ b/example/poisson_dirichlet_example.py @@ -105,7 +105,7 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: misfit_form = PoissonMisfitForm(d, noise_variance) misfit = hpx.NonGaussianContinuousMisfit(Vh, misfit_form, [bc0]) prior_mean = dlx.fem.Function(Vh_m) - prior_mean.x.array[:] = 0.01 + prior_mean.x.array[:] = np.log(2) prior_mean = prior_mean.x prior = hpx.BiLaplacianPrior( @@ -242,43 +242,21 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: prior_samples.append(prior_sample) posterior_samples.append(posterior_sample) - with dlx.io.XDMFFile( + with dlx.io.VTXWriter( msh.comm, - "poisson_Dirichlet_prior_Bilaplacian_samples_prior_np{0:d}.xdmf".format(nproc), - "w", - ) as file: - file.write_mesh(msh) - for i, sample in enumerate(prior_samples): - time_value = i / (num_samples_generate - 1) - file.write_function(sample, time_value) - - with dlx.io.XDMFFile( + "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}.xdmf".format( + "poisson_Dirichlet_prior_Bilaplacian_samples_posterior_np{0:d}.bp".format( nproc ), - "w", - ) as file: - file.write_mesh(msh) - for i, sample in enumerate(posterior_samples): - time_value = i / (num_samples_generate - 1) - file.write_function(sample, time_value) - - # 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) + posterior_samples, + ) as vtx: + vtx.write(0.0) eigen_decomposition_results = {"A": Hmisfit, "B": prior, "k": k, "d": d, "U": U} diff --git a/example/poisson_dirichlet_example_reg.py b/example/poisson_dirichlet_example_reg.py index 2b9e927d..50a8649b 100644 --- a/example/poisson_dirichlet_example_reg.py +++ b/example/poisson_dirichlet_example_reg.py @@ -109,7 +109,7 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: misfit = hpx.NonGaussianContinuousMisfit(Vh, misfit_form, [bc0]) prior_mean = dlx.fem.Function(Vh_m) - prior_mean.x.array[:] = 0.01 + prior_mean.x.array[:] = np.log(2) prior_gamma = prior_param["gamma"] prior_delta = prior_param["delta"] diff --git a/example/poisson_example.py b/example/poisson_example.py index 0ef0e7b2..905887eb 100644 --- a/example/poisson_example.py +++ b/example/poisson_example.py @@ -87,7 +87,7 @@ def run_inversion( misfit_form = PoissonMisfitForm(d, noise_variance) misfit = hpx.NonGaussianContinuousMisfit(Vh, misfit_form) prior_mean = dlx.fem.Function(Vh_m) - prior_mean.x.array[:] = 0.01 + prior_mean.x.array[:] = np.log(2) prior_mean = prior_mean.x prior = hpx.BiLaplacianPrior( @@ -254,9 +254,10 @@ def run_inversion( if __name__ == "__main__": nx = 64 ny = 64 - noise_variance = 1e-4 + noise_variance = 1e-6 prior_param = {"gamma": 0.07, "delta": 0.7} final_results = run_inversion(nx, ny, noise_variance, prior_param) + k, d = ( final_results["eigen_decomposition_results"]["k"], final_results["eigen_decomposition_results"]["d"], diff --git a/example/poisson_example_reg.py b/example/poisson_example_reg.py index 6b8034aa..a658b948 100644 --- a/example/poisson_example_reg.py +++ b/example/poisson_example_reg.py @@ -85,7 +85,7 @@ def run_inversion( misfit = hpx.NonGaussianContinuousMisfit(Vh, misfit_form) prior_mean = dlx.fem.Function(Vh_m) - prior_mean.x.array[:] = 0.01 + prior_mean.x.array[:] = np.log(2) prior_gamma = prior_param["gamma"] prior_delta = prior_param["delta"] @@ -216,7 +216,7 @@ def run_inversion( if __name__ == "__main__": nx = 64 ny = 64 - noise_variance = 1e-4 + noise_variance = 1e-6 prior_param = {"gamma": 0.02, "delta": 0.2} final_results = run_inversion(nx, ny, noise_variance, prior_param) k, d = ( From c31e622ca1c3c8476c8125bc0633412a4da5f42f Mon Sep 17 00:00:00 2001 From: V-Rang Date: Fri, 19 Apr 2024 08:35:22 -0500 Subject: [PATCH 25/29] writing samples in xdmf format working --- example/poisson_dirichlet_example.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/example/poisson_dirichlet_example.py b/example/poisson_dirichlet_example.py index 8e70bae7..19e212ee 100644 --- a/example/poisson_dirichlet_example.py +++ b/example/poisson_dirichlet_example.py @@ -242,6 +242,29 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: prior_samples.append(prior_sample) posterior_samples.append(posterior_sample) + ############################################ + # time series in XDMFFile format: + with dlx.io.XDMFFile( + msh.comm, + "poisson_Dirichlet_prior_Bilaplacian_samples_prior_np{0:d}.xdmf".format(nproc), + "w", + ) as file: + file.write_mesh(msh) + for sample in prior_samples: + file.write_function(sample) + + with dlx.io.XDMFFile( + msh.comm, + "poisson_Dirichlet_prior_Bilaplacian_samples_posterior_np{0:d}.xdmf".format( + nproc + ), + "w", + ) as file: + file.write_mesh(msh) + for sample in posterior_samples: + file.write_function(sample) + ############################################ + with dlx.io.VTXWriter( msh.comm, "poisson_Dirichlet_prior_Bilaplacian_samples_prior_np{0:d}.bp".format(nproc), From 0cb7b58ee2e701b4ff45196076b3e927fe730396 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Fri, 19 Apr 2024 10:02:05 -0500 Subject: [PATCH 26/29] Modified script so that it saves as time series. --- example/poisson_dirichlet_example.py | 80 ++++++++++------------------ 1 file changed, 29 insertions(+), 51 deletions(-) diff --git a/example/poisson_dirichlet_example.py b/example/poisson_dirichlet_example.py index 19e212ee..b12d5124 100644 --- a/example/poisson_dirichlet_example.py +++ b/example/poisson_dirichlet_example.py @@ -227,59 +227,37 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: 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) - + use_vtx = False + prior_sample = dlx.fem.Function(Vh[hpx.PARAMETER], name="prior_sample") + posterior_sample = dlx.fem.Function(Vh[hpx.PARAMETER], name="posterior_sample") + if use_vtx: + with dlx.io.VTXWriter( + msh.comm, + "poisson_Dirichlet_prior_Bilaplacian_samples_np{0:d}.bp".format(nproc), + [prior_sample, posterior_sample] + ) as vtx: + for i in range(num_samples_generate): + hpx.parRandom.normal(1.0, noise) + lap_aprx.sample(noise, prior_sample.x, posterior_sample.x) + prior_sample.x.scatter_forward() + posterior_sample.x.scatter_forward() + vtx.write(float(i)) + else: ############################################ - # time series in XDMFFile format: - with dlx.io.XDMFFile( - msh.comm, - "poisson_Dirichlet_prior_Bilaplacian_samples_prior_np{0:d}.xdmf".format(nproc), - "w", - ) as file: - file.write_mesh(msh) - for sample in prior_samples: - file.write_function(sample) - - with dlx.io.XDMFFile( - msh.comm, - "poisson_Dirichlet_prior_Bilaplacian_samples_posterior_np{0:d}.xdmf".format( - nproc - ), - "w", - ) as file: - file.write_mesh(msh) - for sample in posterior_samples: - file.write_function(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.XDMFFile( + msh.comm, + "poisson_Dirichlet_prior_Bilaplacian_samples_np{0:d}.xdmf".format(nproc), + "w", + ) as file: + file.write_mesh(msh) + for i in range(num_samples_generate): + hpx.parRandom.normal(1.0, noise) + lap_aprx.sample(noise, prior_sample.x, posterior_sample.x) + prior_sample.x.scatter_forward() + posterior_sample.x.scatter_forward() + file.write_function(prior_sample, float(i)) + file.write_function(posterior_sample, float(i)) - 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} From 6f9689874166bd8403ffd730150a25304981c701 Mon Sep 17 00:00:00 2001 From: Umberto Villa Date: Fri, 19 Apr 2024 10:33:17 -0500 Subject: [PATCH 27/29] Rename unused variables --- example/poisson_dirichlet_example.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/example/poisson_dirichlet_example.py b/example/poisson_dirichlet_example.py index b12d5124..26dc415e 100644 --- a/example/poisson_dirichlet_example.py +++ b/example/poisson_dirichlet_example.py @@ -221,8 +221,6 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: 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") From 1026f1130d6bbacfb6da78de6b11a5626f789395 Mon Sep 17 00:00:00 2001 From: V-Rang Date: Fri, 19 Apr 2024 12:01:44 -0500 Subject: [PATCH 28/29] prior and posterior sample time series plots --- example/poisson_dirichlet_example.py | 12 +++--- example/poisson_example.py | 61 ++++++++++++++-------------- example/sfsi_toy_gaussian.py | 61 ++++++++++++++-------------- 3 files changed, 65 insertions(+), 69 deletions(-) diff --git a/example/poisson_dirichlet_example.py b/example/poisson_dirichlet_example.py index 26dc415e..b62d3d69 100644 --- a/example/poisson_dirichlet_example.py +++ b/example/poisson_dirichlet_example.py @@ -221,18 +221,17 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: lap_aprx.mean = prior.generate_parameter(0) lap_aprx.mean.array[:] = x[hpx.PARAMETER].array[:] - noise = prior.generate_parameter("noise") num_samples_generate = 5 - use_vtx = False + use_vtx = False prior_sample = dlx.fem.Function(Vh[hpx.PARAMETER], name="prior_sample") - posterior_sample = dlx.fem.Function(Vh[hpx.PARAMETER], name="posterior_sample") + posterior_sample = dlx.fem.Function(Vh[hpx.PARAMETER], name="posterior_sample") if use_vtx: with dlx.io.VTXWriter( msh.comm, "poisson_Dirichlet_prior_Bilaplacian_samples_np{0:d}.bp".format(nproc), - [prior_sample, posterior_sample] + [prior_sample, posterior_sample], ) as vtx: for i in range(num_samples_generate): hpx.parRandom.normal(1.0, noise) @@ -240,8 +239,8 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: prior_sample.x.scatter_forward() posterior_sample.x.scatter_forward() vtx.write(float(i)) - else: - ############################################ + else: + ############################################ with dlx.io.XDMFFile( msh.comm, "poisson_Dirichlet_prior_Bilaplacian_samples_np{0:d}.xdmf".format(nproc), @@ -256,7 +255,6 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]: file.write_function(prior_sample, float(i)) file.write_function(posterior_sample, float(i)) - eigen_decomposition_results = {"A": Hmisfit, "B": prior, "k": k, "d": d, "U": U} final_results = { diff --git a/example/poisson_example.py b/example/poisson_example.py index 905887eb..01177d37 100644 --- a/example/poisson_example.py +++ b/example/poisson_example.py @@ -202,40 +202,39 @@ def run_inversion( 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) + use_vtx = False + prior_sample = dlx.fem.Function(Vh[hpx.PARAMETER], name="prior_sample") + posterior_sample = dlx.fem.Function(Vh[hpx.PARAMETER], name="posterior_sample") + if use_vtx: + with dlx.io.VTXWriter( + msh.comm, + "poisson_Robin_prior_Bilaplacian_samples_np{0:d}.bp".format(nproc), + [prior_sample, posterior_sample], + ) as vtx: + for i in range(num_samples_generate): + hpx.parRandom.normal(1.0, noise) + lap_aprx.sample(noise, prior_sample.x, posterior_sample.x) + prior_sample.x.scatter_forward() + posterior_sample.x.scatter_forward() + vtx.write(float(i)) + else: + ############################################ + with dlx.io.XDMFFile( + msh.comm, + "poisson_Robin_prior_Bilaplacian_samples_np{0:d}.xdmf".format(nproc), + "w", + ) as file: + file.write_mesh(msh) + for i in range(num_samples_generate): + hpx.parRandom.normal(1.0, noise) + lap_aprx.sample(noise, prior_sample.x, posterior_sample.x) + prior_sample.x.scatter_forward() + posterior_sample.x.scatter_forward() + file.write_function(prior_sample, float(i)) + file.write_function(posterior_sample, float(i)) eigen_decomposition_results = {"A": Hmisfit, "B": prior, "k": k, "d": d, "U": U} diff --git a/example/sfsi_toy_gaussian.py b/example/sfsi_toy_gaussian.py index 5217507e..e63bbb20 100644 --- a/example/sfsi_toy_gaussian.py +++ b/example/sfsi_toy_gaussian.py @@ -230,40 +230,39 @@ def run_inversion( 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) + use_vtx = False + prior_sample = dlx.fem.Function(Vh[hpx.PARAMETER], name="prior_sample") + posterior_sample = dlx.fem.Function(Vh[hpx.PARAMETER], name="posterior_sample") + if use_vtx: + with dlx.io.VTXWriter( + msh.comm, + "pact_prior_Bilaplacian_samples_np{0:d}.bp".format(nproc), + [prior_sample, posterior_sample], + ) as vtx: + for i in range(num_samples_generate): + hpx.parRandom.normal(1.0, noise) + lap_aprx.sample(noise, prior_sample.x, posterior_sample.x) + prior_sample.x.scatter_forward() + posterior_sample.x.scatter_forward() + vtx.write(float(i)) + else: + ############################################ + with dlx.io.XDMFFile( + msh.comm, + "pact_prior_Bilaplacian_samples_np{0:d}.xdmf".format(nproc), + "w", + ) as file: + file.write_mesh(msh) + for i in range(num_samples_generate): + hpx.parRandom.normal(1.0, noise) + lap_aprx.sample(noise, prior_sample.x, posterior_sample.x) + prior_sample.x.scatter_forward() + posterior_sample.x.scatter_forward() + file.write_function(prior_sample, float(i)) + file.write_function(posterior_sample, float(i)) eigen_decomposition_results = {"A": Hmisfit, "B": prior, "k": k, "d": d, "U": U} From 4ef37ef647f6c0cd5c2924508c6c08a3b01ab908 Mon Sep 17 00:00:00 2001 From: V-Rang Date: Wed, 24 Apr 2024 08:30:28 -0500 Subject: [PATCH 29/29] spelling fix --- .github/workflows/CI_testing.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI_testing.yml b/.github/workflows/CI_testing.yml index 718bb6b7..3fc778f3 100644 --- a/.github/workflows/CI_testing.yml +++ b/.github/workflows/CI_testing.yml @@ -40,7 +40,7 @@ jobs: cd ./hippylibX/test && mpirun -n 2 python3 test_lowRankHessian.py - - name: low rank Hessian precondtioner testing + - name: low rank Hessian preconditioner testing run: | cd ./hippylibX/test && mpirun -n 2 python3 test_lowRankHessian_preconditioner.py