Skip to content

Commit

Permalink
Merge pull request #115 from geoscixyz/fix_linear_inversion_app
Browse files Browse the repository at this point in the history 
fix linear inversion app
  • Loading branch information
@sgkang
sgkang authored Jul 26, 2018
2 parents 0b50c0b + c36deaf commit e905862
Show file tree
Hide file tree
Showing 4 changed files with 45 additions and 44 deletions.
2 changes: 1 addition & 1 deletion .bumpversion.cfg
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
[bumpversion]
current_version = 0.0.29
current_version = 0.0.30
files = em_examples/__init__.py setup.py

83 changes: 42 additions & 41 deletions em_examples/LinearInversion.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,8 +65,8 @@ def set_G(
M=100,
p=-0.25,
q=0.25,
k1=1,
kn=60,
j1=1,
jn=60,
):
"""
Parameters
Expand All @@ -79,7 +79,7 @@ def set_G(
self.N=N
self.M=M
self._mesh=Mesh.TensorMesh([M])
jk=np.linspace(k1, kn, N)
jk=np.linspace(j1, jn, N)
self._G=np.zeros((N, self.mesh.nC), dtype=float, order='C')

def g(k):
Expand All @@ -98,19 +98,20 @@ def plot_G(
M=100,
p=-0.25,
q=0.25,
k1=1,
kn=60,
j1=1,
jn=60,
scale='log',
fixed=False,
ymin=-0.001,
ymax=0.011,
scale='log'
ymax=0.011
):
self.set_G(
N=N,
M=M,
p=p,
q=q,
k1=k1,
kn=kn,
j1=j1,
jn=jn,

)

Expand All @@ -125,7 +126,8 @@ def plot_G(
ax2 = plt.subplot(gs1[0, 3:])

ax1.plot(self.mesh.vectorCCx, self.G.T)
ax1.set_ylim(ymin, ymax)
if fixed:
ax1.set_ylim(ymin, ymax)
ax1.set_xlabel("x")
ax1.set_ylabel("g(x)")

Expand All @@ -139,10 +141,6 @@ def plot_G(
ax2.xaxis.set_major_formatter(plt.NullFormatter())
ax2.xaxis.set_minor_formatter(plt.NullFormatter())

# ax2.yaxis.set_major_locator(plt.NullLocator())
# ax2.yaxis.set_minor_locator(plt.NullLocator())
# ax2.yaxis.set_major_formatter(plt.NullFormatter())
# ax2.yaxis.set_minor_formatter(plt.NullFormatter())
plt.tight_layout()
plt.show()

Expand All @@ -168,13 +166,13 @@ def plot_model(
self,
m_background=0.,
m1=1.,
m2=-1.,
m1_center=0.2,
dm1=0.2,
m2=-1.,
m2_center=0.5,
sigma_2=1.,
option="model",
add_noise=False,
add_noise=True,
percentage =10,
floor=1e-1,
):
Expand Down Expand Up @@ -240,8 +238,9 @@ def plot_model(
yerr=self.uncertainty*visualization_factor,
color='k'
)
ax.plot(self.jk, self.data, 'ko')
else:
ax.plot(self.jk, self.data, 'k')
ax.plot(self.jk, self.data, 'ko-')

ax.set_title('Data')
ax.set_xlabel("$k_j$")
Expand Down Expand Up @@ -383,7 +382,7 @@ def plot_inversion(
if run:
axes[0].plot(self.mesh.vectorCCx, self.model[i_plot])
axes[0].set_ylim([-2.5, 2.5])
axes[1].plot(self.jk, self.data, 'k')
axes[1].plot(self.jk, self.data, 'ko')
if run:
axes[1].plot(self.jk, self.pred[i_plot], 'bx')
axes[1].legend(("Observed", "Predicted"))
Expand Down Expand Up @@ -413,13 +412,13 @@ def plot_inversion(

if i_iteration == 0:
i_iteration = 1
axes[2].plot(np.arange(len(self.save.phi_d))[i_iteration-1]+1, self.save.phi_d[i_iteration-1], 'go', ms=10)
axes[2].plot(np.arange(len(self.save.phi_d))[i_iteration-1]+1, self.save.phi_d[i_iteration-1]*2, 'go', ms=10)

ax_1 = axes[2].twinx()
axes[2].semilogy(np.arange(len(self.save.phi_d))+1, self.save.phi_d, 'k-', lw=2)
axes[2].semilogy(np.arange(len(self.save.phi_d))+1, self.save.phi_d*2, 'k-', lw=2)
if self.save.i_target is not None:
axes[2].plot(np.arange(len(self.save.phi_d))[self.save.i_target]+1, self.save.phi_d[self.save.i_target], 'k*', ms=10)
axes[2].plot(np.r_[axes[2].get_xlim()[0], axes[2].get_xlim()[1]], np.ones(2)*self.save.target_misfit, 'k:')
axes[2].plot(np.arange(len(self.save.phi_d))[self.save.i_target]+1, self.save.phi_d[self.save.i_target]*2, 'k*', ms=10)
axes[2].plot(np.r_[axes[2].get_xlim()[0], axes[2].get_xlim()[1]], np.ones(2)*self.save.target_misfit*2, 'k:')

ax_1.semilogy(np.arange(len(self.save.phi_d))+1, self.save.phi_m, 'r', lw=2)
axes[2].set_xlabel("Iteration")
Expand All @@ -440,14 +439,14 @@ def plot_inversion(

if i_iteration == 0:
i_iteration = 1
axes[2].plot(self.save.phi_m[i_iteration-1], self.save.phi_d[i_iteration-1], 'go', ms=10)
axes[2].plot(self.save.phi_m[i_iteration-1], self.save.phi_d[i_iteration-1]*2, 'go', ms=10)

axes[2].plot(self.save.phi_m, self.save.phi_d, 'k-', lw=2)
axes[2].plot(self.save.phi_m, self.save.phi_d*2, 'k-', lw=2)
axes[2].set_xlim(np.hstack(self.save.phi_m).min(), np.hstack(self.save.phi_m).max())
axes[2].set_xlabel("$\phi_m$", fontsize=14)
axes[2].set_ylabel("$\phi_d$", fontsize=14)
if self.save.i_target is not None:
axes[2].plot(self.save.phi_m[self.save.i_target], self.save.phi_d[self.save.i_target], 'k*', ms=10)
axes[2].plot(self.save.phi_m[self.save.i_target], self.save.phi_d[self.save.i_target]*2., 'k*', ms=10)
axes[2].set_title('Tikhonov curve')
plt.tight_layout()

Expand All @@ -456,53 +455,55 @@ def interact_plot_G(self):
self.plot_G,
N=IntSlider(min=1, max=100, step=1, value=20, continuous_update=False),
M=IntSlider(min=1, max=100, step=1, value=100, continuous_update=False),
p =FloatSlider(min=-1, max=0, step=0.05, value=-0.25, continuous_update=False),
p =FloatSlider(min=-1, max=0, step=0.05, value=-0.15, continuous_update=False),
q=FloatSlider(min=0, max=1, step=0.05, value=0.25, continuous_update=False),
k1 =FloatText(value=1.),
kn=FloatText(value=19.),
ymin=FloatText(value=-0.005),
ymax=FloatText(value=0.011),
j1 =FloatText(value=1.),
jn=FloatText(value=19.),
scale=ToggleButtons(
options=["linear", "log"], value="log"
),
fixed=False,
ymin=FloatText(value=-0.005),
ymax=FloatText(value=0.011),

)
return Q

def interact_plot_model(self):
Q=interact(
self.plot_model,
m_background=FloatSlider(
min=-2, max=2, step=0.05, value=0., continuous_update=False
min=-2, max=2, step=0.05, value=0., continuous_update=False, description="m$_{background}$",
),
m1=FloatSlider(
min=-2, max=2, step=0.05, value=1., continuous_update=False
min=-2, max=2, step=0.05, value=1., continuous_update=False, description="m1",
),
m2=FloatSlider(
min=-2, max=2, step=0.05, value=2., continuous_update=False
min=-2, max=2, step=0.05, value=2., continuous_update=False, description="m2",
),
m1_center=FloatSlider(
min=-2, max=2, step=0.05, value=0.2, continuous_update=False
min=-2, max=2, step=0.05, value=0.2, continuous_update=False, description="m1$_{center}$",
),
dm1 =FloatSlider(
min=0, max=0.5, step=0.05, value=0.2, continuous_update=False
min=0, max=0.5, step=0.05, value=0.2, continuous_update=False, description="m1_{width}",
),
m2_center=FloatSlider(
min=-2, max=2, step=0.05, value=0.75, continuous_update=False
min=-2, max=2, step=0.05, value=0.75, continuous_update=False, description="m2$_{center}$",
),
sigma_2=FloatSlider(
min=0.01, max=0.1, step=0.01, value=0.07, continuous_update=False
min=0.01, max=0.1, step=0.01, value=0.07, continuous_update=False, description="m2$_{sigma}$",
),
option=SelectMultiple(
options=["kernel", "model", "data"],
value=["model"],
description='option'
),
percentage=FloatText(value=20),
floor=FloatText(value=0.001),
percentage=FloatText(value=5),
floor=FloatText(value=0.02),
)
return Q

def interact_plot_inversion(self, maxIter=20):
def interact_plot_inversion(self, maxIter=30):
Q = interact(
self.plot_inversion,
maxIter=IntText(value=maxIter),
Expand All @@ -514,7 +515,7 @@ def interact_plot_inversion(self, maxIter=20):
beta0_ratio=FloatText(value=100),
coolingFactor=FloatSlider(min=0.1, max=10, step=1, value=2, continuous_update=False),
coolingRate=IntSlider(min=1, max=10, step=1, value=1, continuous_update=False),
alpha_s=FloatText(value=1),
alpha_s=FloatText(value=1e-10),
alpha_x=FloatText(value=0),
run = True,
target = False,
Expand Down
2 changes: 1 addition & 1 deletion em_examples/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@
from matplotlib import rcParams
rcParams['font.size'] = 16

__version__ = '0.0.29'
__version__ = '0.0.30'
__author__ = 'GeoScixyz developers'
__license__ = 'MIT'
__copyright__ = 'Copyright 2017 GeoScixyz developers'
2 changes: 1 addition & 1 deletion setup.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@

setup(
name = 'em_examples',
version = '0.0.29',
version = '0.0.30',
packages = find_packages(),
install_requires = [
'future',
Expand Down

0 comments on commit e905862

Please sign in to comment.