From fb9979cf0f776496b101d5ac8108199413a1bcbf Mon Sep 17 00:00:00 2001 From: Thibaut Date: Wed, 27 Jul 2022 17:38:30 -0700 Subject: [PATCH] example cleanup, fixing access to properties --- SimPEG/regularization/pgi.py | 30 +++++++++---------- examples/10-pgi/plot_inv_0_PGI_Linear_1D.py | 6 ++-- ...1_PGI_Linear_1D_joint_WithRelationships.py | 11 ------- ...t_inv_1_joint_pf_pgi_full_info_tutorial.py | 9 ++++-- ...lot_inv_2_joint_pf_pgi_no_info_tutorial.py | 5 +++- 5 files changed, 27 insertions(+), 34 deletions(-) diff --git a/SimPEG/regularization/pgi.py b/SimPEG/regularization/pgi.py index 8feaf29756..79590835b6 100644 --- a/SimPEG/regularization/pgi.py +++ b/SimPEG/regularization/pgi.py @@ -71,7 +71,7 @@ def __init__( self.approx_eval = approx_eval self.approx_hessian = approx_hessian self.non_linear_relationships = non_linear_relationships - self.gmm = gmm + self._gmm = copy.deepcopy(gmm) self.wiresmap = wiresmap self.maplist = maplist @@ -127,11 +127,11 @@ def shape(self): def membership(self, m): modellist = self.wiresmap * m model = np.c_[[a * b for a, b in zip(self.maplist, modellist)]].T - return self.gmm.predict(model) # mkvc(m, numDims=2)) + return self.gmm.predict(model) def compute_quasi_geology_model(self): # used once mref is built - mreflist = self.wiresmap * self.mref + mreflist = self.wiresmap * self.reference_model mrefarray = np.c_[[a * b for a, b in zip(self.maplist, mreflist)]].T return np.c_[ [((mrefarray - mean) ** 2).sum(axis=1) for mean in self.gmm.means_] @@ -203,12 +203,12 @@ def __call__(self, m, external_weights=True): W = Identity() if getattr(self, "mref", None) is None: - self.mref = mkvc(self.gmm.means_[self.membership(m)]) + self.reference_model = mkvc(self.gmm.means_[self.membership(m)]) if self.approx_eval: membership = self.compute_quasi_geology_model() dm = self.wiresmap * (m) - dmref = self.wiresmap * (self.mref) + dmref = self.wiresmap * (self.reference_model) dmm = np.c_[[a * b for a, b in zip(self.maplist, dm)]].T if self.non_linear_relationships: dmm = np.r_[ @@ -276,12 +276,12 @@ def __call__(self, m, external_weights=True): def deriv(self, m): if getattr(self, "mref", None) is None: - self.mref = mkvc(self.gmm.means_[self.membership(m)]) + self.reference_model = mkvc(self.gmm.means_[self.membership(m)]) membership = self.compute_quasi_geology_model() modellist = self.wiresmap * m dmmodel = np.c_[[a * b for a, b in zip(self.maplist, modellist)]].T - mreflist = self.wiresmap * self.mref + mreflist = self.wiresmap * self.reference_model mD = [a.deriv(b) for a, b in zip(self.maplist, modellist)] mD = sp.block_diag(mD) @@ -467,7 +467,7 @@ def deriv(self, m): def deriv2(self, m, v=None): if getattr(self, "mref", None) is None: - self.mref = mkvc(self.gmm.means_[self.membership(m)]) + self.reference_model = mkvc(self.gmm.means_[self.membership(m)]) if self.approx_hessian: # we approximate it with the covariance of the cluster @@ -712,6 +712,8 @@ def __init__( self._wiresmap = wiresmap self._maplist = maplist self.regularization_mesh = mesh + self.gmmref = copy.deepcopy(gmmref) + self.gmmref.order_clusters_GM_weight() objfcts = [ PGIsmallness( @@ -771,13 +773,11 @@ def alpha_pgi(self, value): @property def gmm(self): - return self._gmm + return self.objfcts[0].gmm @gmm.setter def gmm(self, gm): - if gm is not None: - self._gmm = copy.deepcopy(gm) - self.objfcts[0].gmm = self.gmm + self.objfcts[0].gmm = copy.deepcopy(gm) def membership(self, m): return self.objfcts[0].membership(m) @@ -833,7 +833,7 @@ def reference_model_in_smooth(self, value: bool): @property def reference_model(self) -> np.ndarray: """Reference physical property model""" - return self._reference_model + return self.objfcts[0].reference_model @reference_model.setter def reference_model(self, values: np.ndarray | float): @@ -844,11 +844,9 @@ def reference_model(self, values: np.ndarray | float): for fct in self.objfcts: fct.reference_model = values - self._reference_model = values - mref = deprecate_property( reference_model, "mref", - "0.x.0", + "reference_model", error=False, ) diff --git a/examples/10-pgi/plot_inv_0_PGI_Linear_1D.py b/examples/10-pgi/plot_inv_0_PGI_Linear_1D.py index 07f9996f97..b733d613fa 100644 --- a/examples/10-pgi/plot_inv_0_PGI_Linear_1D.py +++ b/examples/10-pgi/plot_inv_0_PGI_Linear_1D.py @@ -114,7 +114,7 @@ def g(k): ) # Optimization -opt = optimization.ProjectedGNCG(maxIter=10, maxIterCG=50, tolCG=1e-4) +opt = optimization.ProjectedGNCG(maxIter=20, maxIterCG=50, tolCG=1e-4) opt.remember("xc") # Setup new inverse problem @@ -122,7 +122,7 @@ def g(k): # directives Alphas = directives.AlphasSmoothEstimate_ByEig(alpha0_ratio=10.0, verbose=True) -beta = directives.BetaEstimate_ByEig(beta0_ratio=1e-6) +beta = directives.BetaEstimate_ByEig(beta0_ratio=1e-8) betaIt = directives.PGI_BetaAlphaSchedule( verbose=True, coolingFactor=2.0, @@ -159,7 +159,7 @@ def g(k): axes[2].plot(mesh.cell_centers_x, mtrue, color="black", linewidth=3) axes[2].plot(mesh.cell_centers_x, mnormal, color="blue") axes[2].plot(mesh.cell_centers_x, mcluster, "r-") -axes[2].plot(mesh.cell_centers_x, invProb.reg.objfcts[0].mref, "r--") +axes[2].plot(mesh.cell_centers_x, invProb.reg.objfcts[0].reference_model, "r--") axes[2].legend(("True Model", "L2 Model", "PGI Model", "Learned Mref")) axes[2].set_ylim([-2, 2]) diff --git a/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py b/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py index 1be4efc858..b5b2ab392d 100644 --- a/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py +++ b/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py @@ -333,13 +333,6 @@ def g(k): axes[3].set_xlabel("Property 1") axes[3].set_ylabel("Property 2") -# fig.suptitle( -# 'Doodling with Mapping: one mapping per identified rock unit\n' + -# 'Joint inversion of 1D Linear Problems ' + -# 'with nonlinear petrophysical relationships', -# fontsize=24 -# ) - axes[4].set_axis_off() axes[4].text( 0.5 * (left + right), @@ -353,8 +346,6 @@ def g(k): ) axes[5].plot(mesh.cell_centers_x, wires.m1 * mcluster_no_map, "b.-", ms=5, marker="v") -# axes[5].plot(mesh.cell_centers_x, wires.m1 * reg_simple_no_map.objfcts[0].mref, 'g--') - axes[5].plot(mesh.cell_centers_x, wires.m1 * m, "k--") axes[5].set_title("Problem 1") axes[5].legend(["Recovered Model", "True Model"], loc=1) @@ -362,8 +353,6 @@ def g(k): axes[5].set_ylabel("Property 1") axes[6].plot(mesh.cell_centers_x, wires.m2 * mcluster_no_map, "r.-", ms=5, marker="v") -# axes[6].plot(mesh.cell_centers_x, wires.m2 * reg_simple_no_map.objfcts[0].mref, 'g--') - axes[6].plot(mesh.cell_centers_x, wires.m2 * m, "k--") axes[6].set_title("Problem 2") axes[6].legend(["Recovered Model", "True Model"], loc=1) diff --git a/tutorials/13-joint_inversion/plot_inv_1_joint_pf_pgi_full_info_tutorial.py b/tutorials/13-joint_inversion/plot_inv_1_joint_pf_pgi_full_info_tutorial.py index b7d8afa234..4f3b4769ce 100644 --- a/tutorials/13-joint_inversion/plot_inv_1_joint_pf_pgi_full_info_tutorial.py +++ b/tutorials/13-joint_inversion/plot_inv_1_joint_pf_pgi_full_info_tutorial.py @@ -334,7 +334,10 @@ alpha_xx=0.0, alpha_yy=0.0, alpha_zz=0.0, - weights_list=[wr_grav, wr_mag], # weights each phys. prop. by correct sensW + # use the classification of the initial model (here, all background unit) + # as initial reference model + reference_model=utils.mkvc(gmmref.means_[gmmref.predict(m0.reshape(actvMap.nP,-1))]), + weights_list=[wr_grav, wr_mag], # weights each phys. prop. by correct sensW ) ######################################################################### @@ -347,8 +350,8 @@ # ratio to use for each phys prop. smoothness in each direction; # roughly the ratio of the order of magnitude of each phys. prop. alpha0_ratio = np.r_[ - 1e-2 * np.ones(len(reg.objfcts[1].objfcts[1:])), - 1e-2 * 100.0 * np.ones(len(reg.objfcts[2].objfcts[1:])), + 1e-4 * np.ones(len(reg.objfcts[1].objfcts[1:])), + 1e-4 * 100.0 * np.ones(len(reg.objfcts[2].objfcts[1:])), ] Alphas = directives.AlphasSmoothEstimate_ByEig(alpha0_ratio=alpha0_ratio, verbose=True) # initialize beta and beta/alpha_s schedule diff --git a/tutorials/13-joint_inversion/plot_inv_2_joint_pf_pgi_no_info_tutorial.py b/tutorials/13-joint_inversion/plot_inv_2_joint_pf_pgi_no_info_tutorial.py index 760659b2f5..2024f3ef23 100644 --- a/tutorials/13-joint_inversion/plot_inv_2_joint_pf_pgi_no_info_tutorial.py +++ b/tutorials/13-joint_inversion/plot_inv_2_joint_pf_pgi_no_info_tutorial.py @@ -354,7 +354,10 @@ alpha_xx=0.0, alpha_yy=0.0, alpha_zz=0.0, - weights_list=[wr_grav, wr_mag], # weights each phys. prop. by each sensW + # use the classification of the initial model (here, all background unit) + # as initial reference model + reference_model=utils.mkvc(gmmref.means_[gmmref.predict(m0.reshape(actvMap.nP,-1))]), + weights_list=[wr_grav, wr_mag], # weights each phys. prop. by correct sensW ) # Directives