Coverage for /usr/local/lib/python3.11/dist-packages/grond/problems/cmt/problem.py: 72%
156 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-10-26 16:25 +0000
1import numpy as num
2import math
3import logging
5from pyrocko import gf, util, moment_tensor as mtm
6from pyrocko.guts import String, Float, Dict, StringChoice
8from grond.meta import Forbidden, expand_template, Parameter, \
9 has_get_plot_classes
11from ..base import Problem, ProblemConfig
13guts_prefix = 'grond'
14logger = logging.getLogger('grond.problems.cmt.problem')
15km = 1e3
16as_km = dict(scale_factor=km, scale_unit='km')
19def as_arr(mat_or_arr):
20 try:
21 return mat_or_arr.A
22 except AttributeError:
23 return mat_or_arr
26class MTType(StringChoice):
27 choices = ['full', 'deviatoric', 'dc']
30class STFType(StringChoice):
31 choices = ['HalfSinusoidSTF', 'ResonatorSTF']
33 cls = {
34 'HalfSinusoidSTF': gf.HalfSinusoidSTF,
35 'ResonatorSTF': gf.ResonatorSTF}
37 @classmethod
38 def base_stf(cls, name):
39 return cls.cls[name]()
42class CMTProblemConfig(ProblemConfig):
44 ranges = Dict.T(String.T(), gf.Range.T())
45 distance_min = Float.T(default=0.0)
46 mt_type = MTType.T(default='full')
47 stf_type = STFType.T(default='HalfSinusoidSTF')
49 def get_problem(self, event, target_groups, targets):
50 self.check_deprecations()
52 if event.depth is None:
53 event.depth = 0.
55 base_source = gf.MTSource.from_pyrocko_event(event)
57 stf = STFType.base_stf(self.stf_type)
58 stf.duration = event.duration or 0.0
60 base_source.stf = stf
62 subs = dict(
63 event_name=event.name,
64 event_time=util.time_to_str(event.time))
66 problem = CMTProblem(
67 name=expand_template(self.name_template, subs),
68 base_source=base_source,
69 target_groups=target_groups,
70 targets=targets,
71 ranges=self.ranges,
72 distance_min=self.distance_min,
73 mt_type=self.mt_type,
74 stf_type=self.stf_type,
75 norm_exponent=self.norm_exponent)
77 return problem
80@has_get_plot_classes
81class CMTProblem(Problem):
83 problem_parameters = [
84 Parameter('time', 's', label='Time'),
85 Parameter('north_shift', 'm', label='Northing', **as_km),
86 Parameter('east_shift', 'm', label='Easting', **as_km),
87 Parameter('depth', 'm', label='Depth', **as_km),
88 Parameter('magnitude', label='Magnitude'),
89 Parameter('rmnn', label='$m_{nn} / M_0$'),
90 Parameter('rmee', label='$m_{ee} / M_0$'),
91 Parameter('rmdd', label='$m_{dd} / M_0$'),
92 Parameter('rmne', label='$m_{ne} / M_0$'),
93 Parameter('rmnd', label='$m_{nd} / M_0$'),
94 Parameter('rmed', label='$m_{ed} / M_0$')]
96 problem_parameters_stf = {
97 'HalfSinusoidSTF': [
98 Parameter('duration', 's', label='Duration')],
99 'ResonatorSTF': [
100 Parameter('duration', 's', label='Duration'),
101 Parameter('frequency', 'Hz', label='Frequency')]}
103 dependants = [
104 Parameter('strike1', u'\u00b0', label='Strike 1'),
105 Parameter('dip1', u'\u00b0', label='Dip 1'),
106 Parameter('rake1', u'\u00b0', label='Rake 1'),
107 Parameter('strike2', u'\u00b0', label='Strike 2'),
108 Parameter('dip2', u'\u00b0', label='Dip 2'),
109 Parameter('rake2', u'\u00b0', label='Rake 2'),
110 Parameter('rel_moment_iso', label='$M_{0}^{ISO}/M_{0}$'),
111 Parameter('rel_moment_clvd', label='$M_{0}^{CLVD}/M_{0}$')]
113 distance_min = Float.T(default=0.0)
114 mt_type = MTType.T(default='full')
115 stf_type = STFType.T(default='HalfSinusoidSTF')
117 def __init__(self, **kwargs):
118 Problem.__init__(self, **kwargs)
119 self.deps_cache = {}
120 self.problem_parameters = self.problem_parameters \
121 + self.problem_parameters_stf[self.stf_type]
122 self._base_stf = STFType.base_stf(self.stf_type)
124 def get_stf(self, d):
125 d_stf = {}
126 for p in self.problem_parameters_stf[self.stf_type]:
127 d_stf[p.name] = float(d[p.name])
129 return self._base_stf.clone(**d_stf)
131 def get_source(self, x):
132 d = self.get_parameter_dict(x)
133 rm6 = num.array([d.rmnn, d.rmee, d.rmdd, d.rmne, d.rmnd, d.rmed],
134 dtype=float)
136 m0 = mtm.magnitude_to_moment(d.magnitude)
137 m6 = rm6 * m0
139 p = {}
140 for k in self.base_source.keys():
141 if k in d:
142 p[k] = float(
143 self.ranges[k].make_relative(self.base_source[k], d[k]))
145 source = self.base_source.clone(m6=m6, stf=self.get_stf(d), **p)
146 return source
148 def make_dependant(self, xs, pname):
149 cache = self.deps_cache
150 if xs.ndim == 1:
151 return self.make_dependant(xs[num.newaxis, :], pname)[0]
153 if pname not in self.dependant_names:
154 raise KeyError(pname)
156 mt = self.base_source.pyrocko_moment_tensor()
158 sdrs_ref = mt.both_strike_dip_rake()
160 y = num.zeros(xs.shape[0])
161 for i, x in enumerate(xs):
162 k = tuple(x.tolist())
163 if k not in cache:
164 source = self.get_source(x)
165 mt = source.pyrocko_moment_tensor()
166 res = mt.standard_decomposition()
167 sdrs = mt.both_strike_dip_rake()
168 if sdrs_ref:
169 sdrs = mtm.order_like(sdrs, sdrs_ref)
171 cache[k] = mt, res, sdrs
173 mt, res, sdrs = cache[k]
175 if pname == 'rel_moment_iso':
176 ratio_iso, m_iso = res[0][1:3]
177 y[i] = ratio_iso * num.sign(m_iso[0, 0])
178 elif pname == 'rel_moment_clvd':
179 ratio_clvd, m_clvd = res[2][1:3]
180 evals, evecs = mtm.eigh_check(m_clvd)
181 ii = num.argmax(num.abs(evals))
182 y[i] = ratio_clvd * num.sign(evals[ii])
183 else:
184 isdr = {'strike': 0, 'dip': 1, 'rake': 2}[pname[:-1]]
185 y[i] = sdrs[int(pname[-1])-1][isdr]
187 return y
189 def pack_stf(self, stf):
190 return [
191 stf[p.name] for p in self.problem_parameters_stf[self.stf_type]]
193 def pack(self, source):
194 m6 = source.m6
195 mt = source.pyrocko_moment_tensor()
196 rm6 = m6 / mt.scalar_moment()
198 x = num.array([
199 source.time - self.base_source.time,
200 source.north_shift,
201 source.east_shift,
202 source.depth,
203 mt.moment_magnitude(),
204 ] + rm6.tolist() + self.pack_stf(source.stf), dtype=float)
206 return x
208 def event_to_x(self, event):
209 bs = self.base_source
210 source = bs.__class__.from_pyrocko_event(event)
212 if source.stf is None:
213 rstatem = self.get_rstate_manager()
214 rstate = rstatem.get_rstate(name='event2x')
216 source.stf = bs.stf
217 source.stf.duration = rstate.uniform(
218 low=self.ranges['time'].start,
219 high=self.ranges['time'].stop,
220 size=1)
222 return self.source_to_x(source)
224 def random_uniform(self, xbounds, rstate, fixed_magnitude=None):
226 x = num.zeros(self.nparameters)
227 for i in range(self.nparameters):
228 x[i] = rstate.uniform(xbounds[i, 0], xbounds[i, 1])
230 if fixed_magnitude is not None:
231 x[4] = fixed_magnitude
233 x[5:11] = mtm.random_m6(x=rstate.random_sample(6))
235 return x.tolist()
237 def preconstrain(self, x, optimizer=False):
238 d = self.get_parameter_dict(x)
239 m6 = num.array([d.rmnn, d.rmee, d.rmdd, d.rmne, d.rmnd, d.rmed],
240 dtype=float)
242 m9 = mtm.symmat6(*m6)
243 if self.mt_type == 'deviatoric':
244 trace_m = num.trace(m9)
245 m_iso = num.diag([trace_m / 3., trace_m / 3., trace_m / 3.])
246 m9 -= m_iso
248 elif self.mt_type == 'dc':
249 mt = mtm.MomentTensor(m=m9)
250 m9 = mt.standard_decomposition()[1][2]
252 m0_unscaled = math.sqrt(num.sum(as_arr(m9)**2)) / math.sqrt(2.)
254 m9 /= m0_unscaled
255 m6 = mtm.to6(m9)
256 d.rmnn, d.rmee, d.rmdd, d.rmne, d.rmnd, d.rmed = m6
257 x = self.get_parameter_array(d)
259 source = self.get_source(x)
260 for t in self.waveform_targets:
261 if (self.distance_min > num.asarray(t.distance_to(source))).any():
262 raise Forbidden()
264 return x
266 def get_dependant_bounds(self):
267 out = [
268 (0., 360.),
269 (0., 90.),
270 (-180., 180.),
271 (0., 360.),
272 (0., 90.),
273 (-180., 180.),
274 (-1., 1.),
275 (-1., 1.)]
277 return out
279 @classmethod
280 def get_plot_classes(cls):
281 from .. import plot
282 plots = super(CMTProblem, cls).get_plot_classes()
283 plots.extend([plot.HudsonPlot, plot.MTDecompositionPlot,
284 plot.MTLocationPlot, plot.MTFuzzyPlot])
285 return plots
288__all__ = '''
289 CMTProblem
290 CMTProblemConfig
291'''.split()