Skip to content

example cleanup, fixing access to properties #23

New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Merged
domfournier merged 1 commit into from
Jul 28, 2022
Merged

example cleanup, fixing access to properties #23

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
30 changes: 14 additions & 16 deletions SimPEG/regularization/pgi.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Copy link
Copy Markdown

@domfournier domfournier Jul 28, 2022

Choose a reason for hiding this comment

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

Should we do some checks on what the input is? Hence using the setter that does the copy as well?

Copy link
Copy Markdown
Author

@thibaut-kobold thibaut-kobold Jul 28, 2022 •
edited
Loading

Choose a reason for hiding this comment

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

We could, what do you have in mind?

Some thoughts / insights:

Admittedly, this option is not often use. Its main goal is to be able to restart a PGI where the GMM is being learnt (thus the difference between the reference GMM and the GMM model). gmmref is really the most important here.

The reason I do the deepcopy is that (remembering things from several years ago) I was having issues when pointing to the original object, so it was best to have only internal to the regularization (plus, it is potentially getting modified every iteration, something that the GMM class was not originally designed for, thus the creation of a new object that overwrite the previous).

Copy link
Copy Markdown

@domfournier domfournier Jul 28, 2022 •
edited
Loading

Choose a reason for hiding this comment

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

Yep for sure. I just think we should check for at least the type and shape of the input to make sure the code can fail fast with a nice message if the input doesn't make sense. We should probably do some typing of the return value as well. I don't even know myself what the input/output looks like.

Copy link
Copy Markdown
Author

@thibaut-kobold thibaut-kobold Jul 28, 2022

Choose a reason for hiding this comment

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

gmm is an object of class "WeightedGaussianMixture" (or one of its children class). We would need to have a validate method somewhere then.
We can think of something. This is a very quick first pass. I just checked that tests and examples were running. Up to you if you want me to have all in a single PR or if you want to merge this one and I create others as I go.

self.wiresmap = wiresmap
self.maplist = maplist

Expand Down Expand Up @@ -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_]
Expand Down Expand Up @@ -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_[
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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(
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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):
Expand All @@ -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,
)
6 changes: 3 additions & 3 deletions examples/10-pgi/plot_inv_0_PGI_Linear_1D.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -114,15 +114,15 @@ 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
invProb = inverse_problem.BaseInvProblem(dmis, reg, opt)

# 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,
Expand Down Expand Up @@ -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])
Expand Down
11 changes: 0 additions & 11 deletions examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Expand All @@ -353,17 +346,13 @@ 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)
axes[5].set_xlabel("X")
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)
Expand Down
9 changes: 6 additions & 3 deletions tutorials/13-joint_inversion/plot_inv_1_joint_pf_pgi_full_info_tutorial.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
)

#########################################################################
Expand All @@ -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
Expand Down
5 changes: 4 additions & 1 deletion tutorials/13-joint_inversion/plot_inv_2_joint_pf_pgi_no_info_tutorial.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down