forked from simpeg/simpeg
-
Notifications
You must be signed in to change notification settings - Fork 1
/
plot_inversion_linear_irls.py
151 lines (115 loc) · 3.92 KB
/
plot_inversion_linear_irls.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
"""
Inversion: Linear: IRLS
=======================
Here we go over the basics of creating a linear problem and inversion.
"""
from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt
from SimPEG import Mesh
from SimPEG import Problem
from SimPEG import Survey
from SimPEG import DataMisfit
from SimPEG import Directives
from SimPEG import Optimization
from SimPEG import Regularization
from SimPEG import InvProblem
from SimPEG import Inversion
from SimPEG import Maps
def run(N=100, plotIt=True):
np.random.seed(1)
std_noise = 1e-2
mesh = Mesh.TensorMesh([N])
m0 = np.ones(mesh.nC) * 1e-4
mref = np.zeros(mesh.nC)
nk = 20
jk = np.linspace(1., 60., nk)
p = -0.25
q = 0.25
def g(k):
return (
np.exp(p*jk[k]*mesh.vectorCCx) *
np.cos(np.pi*q*jk[k]*mesh.vectorCCx)
)
G = np.empty((nk, mesh.nC))
for i in range(nk):
G[i, :] = g(i)
mtrue = np.zeros(mesh.nC)
mtrue[mesh.vectorCCx > 0.3] = 1.
mtrue[mesh.vectorCCx > 0.45] = -0.5
mtrue[mesh.vectorCCx > 0.6] = 0
prob = Problem.LinearProblem(mesh, G=G)
survey = Survey.LinearSurvey()
survey.pair(prob)
survey.dobs = prob.fields(mtrue) + std_noise * np.random.randn(nk)
wd = np.ones(nk) * std_noise
# Distance weighting
wr = np.sum(prob.getJ(m0)**2., axis=0)**0.5
wr = wr/np.max(wr)
dmis = DataMisfit.l2_DataMisfit(survey)
dmis.W = 1./wd
betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e0)
# Creat reduced identity map
idenMap = Maps.IdentityMap(nP=mesh.nC)
reg = Regularization.Sparse(mesh, mapping=idenMap)
reg.mref = mref
reg.cell_weights = wr
reg.norms = np.c_[0., 0., 2., 2.]
reg.mref = np.zeros(mesh.nC)
opt = Optimization.ProjectedGNCG(
maxIter=100, lower=-2., upper=2.,
maxIterLS=20, maxIterCG=10, tolCG=1e-3
)
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
update_Jacobi = Directives.UpdatePreconditioner()
# Set the IRLS directive, penalize the lowest 25 percentile of model values
# Start with an l2-l2, then switch to lp-norms
IRLS = Directives.Update_IRLS(
maxIRLSiter=40, minGNiter=1, f_min_change=1e-4)
saveDict = Directives.SaveOutputEveryIteration(save_txt=False)
inv = Inversion.BaseInversion(
invProb,
directiveList=[IRLS, betaest, update_Jacobi, saveDict]
)
# Run inversion
mrec = inv.run(m0)
print("Final misfit:" + str(invProb.dmisfit(mrec)))
if plotIt:
fig, axes = plt.subplots(2, 2, figsize=(12*1.2, 8*1.2))
for i in range(prob.G.shape[0]):
axes[0, 0].plot(prob.G[i, :])
axes[0, 0].set_title('Columns of matrix G')
axes[0, 1].plot(mesh.vectorCCx, mtrue, 'b-')
axes[0, 1].plot(mesh.vectorCCx, invProb.l2model, 'r-')
# axes[0, 1].legend(('True Model', 'Recovered Model'))
axes[0, 1].set_ylim(-1.0, 1.25)
axes[0, 1].plot(mesh.vectorCCx, mrec, 'k-', lw=2)
axes[0, 1].legend(
(
'True Model',
'Smooth l2-l2',
'Sparse norms: {0}'.format(*reg.norms)
),
fontsize=12
)
axes[1, 1].plot(saveDict.phi_d, 'k', lw=2)
twin = axes[1, 1].twinx()
twin.plot(saveDict.phi_m, 'k--', lw=2)
axes[1, 1].plot(
np.r_[IRLS.iterStart, IRLS.iterStart],
np.r_[0, np.max(saveDict.phi_d)], 'k:'
)
axes[1, 1].text(
IRLS.iterStart, 0.,
'IRLS Start', va='bottom', ha='center',
rotation='vertical', size=12,
bbox={'facecolor': 'white'}
)
axes[1, 1].set_ylabel('$\phi_d$', size=16, rotation=0)
axes[1, 1].set_xlabel('Iterations', size=14)
axes[1, 0].axis('off')
twin.set_ylabel('$\phi_m$', size=16, rotation=0)
return prob, survey, mesh, mrec
if __name__ == '__main__':
run()
plt.show()