-
Notifications
You must be signed in to change notification settings - Fork 33
Expand file tree
/
Copy paththeory.py
More file actions
232 lines (201 loc) · 8.23 KB
/
theory.py
File metadata and controls
232 lines (201 loc) · 8.23 KB
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
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
import warnings
import numpy as np
from srfpython.Herrmann.Herrmann import HerrmannCallerBasis
from srfpython.HerrMet.parameterizers import Parameterizer
from srfpython.HerrMet.datacoders import Datacoder
"""
Parameterizer : object to convert an array of parameters into smth that can be understood by Herrmann.dispersion
Theory : object to convert a model array into a data array
Datacoder : object to convert output from Herrmann.dispersion to an array of data
a model array (m)
| ^ |
| | |
| mod96file -----> parameterizer --------> m_apr, CM
| (apriori) | |
| | v
| depth model
| (ztop, vp, vs, rh)
| |
theory Herrmann.HerrmannCaller.disperse
(forward problem) |
(Frechet derivatives) v
| dispersion data
| (waves, types, modes, freqs, values, (dvalues))
| ^ |
| | |
| surf96file -----> datacoder --------> d_obs, CD => logRHOD
| (target) | |
v | v
a data array (d)
"""
class Theory(object):
def __init__(self, parameterizer, datacoder, h=0.005, ddc=0.005):
if not isinstance(parameterizer, Parameterizer):
raise TypeError(type(parameterizer))
if not isinstance(datacoder, Datacoder):
raise TypeError(type(datacoder))
self.parameterizer, self.datacoder = parameterizer, datacoder
self.herrmanncaller = HerrmannCallerBasis(
waves=datacoder.waves, types=datacoder.types,
modes=datacoder.modes, freqs=datacoder.freqs,
h=h, ddc=ddc)
# the output of self.herrmanncaller.disperse must match excatcly the datacoder order
assert (self.herrmanncaller.waves == datacoder.waves).all()
assert (self.herrmanncaller.types == datacoder.types).all()
assert (self.herrmanncaller.modes == datacoder.modes).all()
assert (self.herrmanncaller.freqs == datacoder.freqs).all()
def __call__(self, m):
"""solves the forward problem"""
# recover model from parameterized array (m)
ztop, vp, vs, rh = self.parameterizer.inv(m)
# call Herrmann's codes
values = self.herrmanncaller.disperse(ztop, vp, vs, rh)
# convert and return dispersion data to coded array (d)
return self.datacoder(values)
def frechet_derivatives(self, m, gm=None):
"""
:param m: the model near which to compute frechet derivatives
:param gm: if the dispersion data is known already (otherwise, I compute it using self.__call__(m)
:return fd: the frechet derivatives
"""
# compute the data for the current model if not provided
if gm is None:
gm = self.__call__(m)
if np.isnan(gm).any():
# probably because some freq points are above the cut off period
warnings.warn('g(m) contains nans')
# raise ValueError('g(m) contains nans')
if np.isinf(gm).any():
raise ValueError('g(m) contains infs')
# shift the current model parameters by a delta
deltam = self.parameterizer.frechet_deltas()
if not np.all(deltam > 0.):
raise ValueError('deltam must be positive')
fd = np.zeros((len(gm), len(m)), float)
for j in range(len(m)):
mj = m.copy()
mj[j] += deltam[j]
gmj = self.__call__(mj)
if np.isnan(gmj).any() or np.isinf(gmj).any():
# probably because some freq points are above the cut off period
warnings.warn('g(m+dm) contains nans of infs')
# raise ValueError('g(m+dm) contains nans of infs')
fd[:, j] = (gmj - gm) / deltam[j]
fd[np.isnan(fd)] = 0.
return fd
if __name__ == '__main__':
from srfpython.standalone.asciifile import AsciiFile_fromstring
from srfpython.HerrMet.parameterizers import Parameterizer_mZVSVPvsRHvp
from srfpython.HerrMet.datacoders import Datacoder_log, makedatacoder
if False:
parameter_file_string = """
#met NLAYER = 9
#met TYPE = 'mZVSPRRH'
#met PRIORTYPE = "DVPDVSDRHDPR"
#met DVSMIN = -0.5
#met DVSMAX = +1.5
#met DVPMIN = -0.5
#met DVPMAX = +1.5
#met DPRMIN = -1.0
#met DPRMAX = +0.0
#met DRHMIN = +0.0
#met DRHMAX = +1.0
#fld KEY VINF VSUP
#unt [] [] []
#fmt %s %f %f
-Z1 -0.200000 -0.05
-Z2 -0.400000 -0.2
-Z3 -0.800000 -0.400000
-Z4 -1.000000 -0.800000
-Z5 -1.500000 -1.000000
-Z6 -2.000000 -1.500000
-Z7 -2.500000 -2.000000
-Z8 -3.000000 -2.500000
VS0 0.2 1.500000
VS1 0.2 1.700000
VS2 0.5 2.000000
VS3 0.5 2.500000
VS4 1.0 3.000000
VS5 1.0 3.250000
VS6 1.3 3.500000
VS7 1.3 3.550000
VS8 3.2 3.600000
PR0 1.800000 3.6
PR1 1.800000 2.3
PR2 1.800000 2.3
PR3 1.800000 2.3
PR4 1.700000 2.3
PR5 1.700000 2.3
PR6 1.600000 2.0
PR7 1.600000 2.0
PR8 1.600000 1.9
RH0 1.80 2.3
RH1 1.80 2.3
RH2 1.90 2.4
RH3 2.00 2.4
RH4 2.00 2.5
RH5 2.10 2.5
RH6 2.30 2.7
RH7 2.40 2.7
RH8 2.50 2.7
"""
else:
parameter_file_string = """
# BROCHER2005
#met RHvp = 'def RHvp(VP):return 1.6612 * VP - 0.4721 * VP ** 2.0 + 0.0671 * VP ** 3.0 - 0.0043 * VP ** 4.0 + 0.000106 * VP ** 5.0'
#met NLAYER = 10
#met TYPE = 'mZVSVPvsRHvp'
#met VPvs = 'def VPvs(VS):return 0.9409 + 2.0947 * VS - 0.8206 * VS ** 2.0 + 0.2683 * VS ** 3.0 - 0.0251 * VS ** 4.0'
#fld KEY VINF VSUP
#unt - - -
#fmt %5s %16f %16f
-Z1 -0.333333 -0.001000
-Z2 -0.666667 -0.333333
-Z3 -1.000000 -0.666667
-Z4 -1.333333 -1.000000
-Z5 -1.666667 -1.333333
-Z6 -2.000000 -1.666667
-Z7 -2.333333 -2.000000
-Z8 -2.666667 -2.333333
-Z9 -3.000000 -2.666667
VS0 1.000000 1.100000
VS1 1.100000 1.200000
VS2 1.200000 1.300000
VS3 1.300000 1.400000
VS4 1.400000 1.500000
VS5 1.500000 1.600000
VS6 1.600000 2.000000
VS7 2.000000 2.500000
VS8 2.500000 3.000000
VS9 3.000000 3.1
"""
dispersion_file_string = """SURF96 R U T 0 5.000000 2.368112 0.1
SURF96 R U T 0 4.429334 2.289762 0.1
SURF96 R U T 0 3.923800 2.184182 0.1
SURF96 R U T 0 3.475964 2.028462 0.1
SURF96 R U T 0 3.079241 1.787782 0.1
SURF96 R U T 0 2.727797 1.457564 0.1
SURF96 R U T 0 2.416465 1.165715 0.1
SURF96 R U T 0 2.140666 1.051432 0.1
SURF96 R U T 0 1.896345 1.051072 0.1
SURF96 R U T 0 1.679909 1.085069 0.1
SURF96 R U T 0 1.488176 1.126465 0.1
SURF96 R U T 0 1.318325 1.168252 0.1
SURF96 R U T 0 1.167861 1.206051 0.1
SURF96 R U T 0 1.034569 1.234272 0.1
SURF96 R U T 0 0.916490 1.246419 0.1
SURF96 R U T 0 0.811888 1.235000 0.1
SURF96 R U T 0 0.719225 1.192478 0.1
SURF96 R U T 0 0.637137 1.115804 0.1
SURF96 R U T 0 0.564419 1.017582 0.1
SURF96 R U T 0 0.500000 0.929962 0.1
"""
A = AsciiFile_fromstring(parameter_file_string)
# p = Parameterizer_mZVSPRRH(A)
p = Parameterizer_mZVSVPvsRHvp(A)
d = makedatacoder(surf96filename=dispersion_file_string, which=Datacoder_log)
t = Theory(parameterizer=p, datacoder=d)
fd = t.frechet_derivatives(m=p.MMEAN, gm=None)
import matplotlib.pyplot as plt
plt.colorbar(plt.imshow(fd))
plt.show()