Coverage for /usr/local/lib/python3.11/dist-packages/grond/problems/cmt/problem.py: 76%
152 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-11-01 12:39 +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 MTRandomType(StringChoice):
31 choices = ['uniform', 'use_bounds']
34class STFType(StringChoice):
35 choices = ['HalfSinusoidSTF', 'ResonatorSTF']
37 cls = {
38 'HalfSinusoidSTF': gf.HalfSinusoidSTF,
39 'ResonatorSTF': gf.ResonatorSTF}
41 @classmethod
42 def base_stf(cls, name):
43 return cls.cls[name]()
46class CMTProblemConfig(ProblemConfig):
48 ranges = Dict.T(String.T(), gf.Range.T())
49 distance_min = Float.T(default=0.0)
50 mt_type = MTType.T(default='full')
51 mt_random = MTRandomType.T(default='uniform')
52 stf_type = STFType.T(default='HalfSinusoidSTF')
54 def get_problem(self, event, target_groups, targets):
55 self.check_deprecations()
57 if event.depth is None:
58 event.depth = 0.
60 base_source = gf.MTSource.from_pyrocko_event(event)
62 stf = STFType.base_stf(self.stf_type)
63 stf.duration = event.duration or 0.0
65 base_source.stf = stf
67 subs = dict(
68 event_name=event.name,
69 event_time=util.time_to_str(event.time))
71 problem = CMTProblem(
72 name=expand_template(self.name_template, subs),
73 base_source=base_source,
74 target_groups=target_groups,
75 targets=targets,
76 ranges=self.ranges,
77 distance_min=self.distance_min,
78 mt_type=self.mt_type,
79 mt_random=self.mt_random,
80 stf_type=self.stf_type,
81 norm_exponent=self.norm_exponent)
83 return problem
86@has_get_plot_classes
87class CMTProblem(Problem):
89 problem_parameters = [
90 Parameter('time', 's', label='Time'),
91 Parameter('north_shift', 'm', label='Northing', **as_km),
92 Parameter('east_shift', 'm', label='Easting', **as_km),
93 Parameter('depth', 'm', label='Depth', **as_km),
94 Parameter('magnitude', label='Magnitude'),
95 Parameter('rmnn', label='$m_{nn} / M_0$'),
96 Parameter('rmee', label='$m_{ee} / M_0$'),
97 Parameter('rmdd', label='$m_{dd} / M_0$'),
98 Parameter('rmne', label='$m_{ne} / M_0$'),
99 Parameter('rmnd', label='$m_{nd} / M_0$'),
100 Parameter('rmed', label='$m_{ed} / M_0$')]
102 problem_parameters_stf = {
103 'HalfSinusoidSTF': [
104 Parameter('duration', 's', label='Duration')],
105 'ResonatorSTF': [
106 Parameter('duration', 's', label='Duration'),
107 Parameter('frequency', 'Hz', label='Frequency')]}
109 dependants = [
110 Parameter('strike1', u'\u00b0', label='Strike 1'),
111 Parameter('dip1', u'\u00b0', label='Dip 1'),
112 Parameter('rake1', u'\u00b0', label='Rake 1'),
113 Parameter('strike2', u'\u00b0', label='Strike 2'),
114 Parameter('dip2', u'\u00b0', label='Dip 2'),
115 Parameter('rake2', u'\u00b0', label='Rake 2'),
116 Parameter('rel_moment_iso', label='$M_{0}^{ISO}/M_{0}$'),
117 Parameter('rel_moment_clvd', label='$M_{0}^{CLVD}/M_{0}$')]
119 distance_min = Float.T(default=0.0)
120 mt_type = MTType.T(default='full')
121 mt_random = MTRandomType.T(default='uniform')
122 stf_type = STFType.T(default='HalfSinusoidSTF')
124 def __init__(self, **kwargs):
125 Problem.__init__(self, **kwargs)
126 self.deps_cache = {}
127 self.problem_parameters = self.problem_parameters \
128 + self.problem_parameters_stf[self.stf_type]
129 self._base_stf = STFType.base_stf(self.stf_type)
131 def get_stf(self, d):
132 d_stf = {}
133 for p in self.problem_parameters_stf[self.stf_type]:
134 d_stf[p.name] = float(d[p.name])
136 return self._base_stf.clone(**d_stf)
138 def get_source(self, x):
139 d = self.get_parameter_dict(x)
140 rm6 = num.array([d.rmnn, d.rmee, d.rmdd, d.rmne, d.rmnd, d.rmed],
141 dtype=float)
143 m0 = mtm.magnitude_to_moment(d.magnitude)
144 m6 = rm6 * m0
146 p = {}
147 for k in self.base_source.keys():
148 if k in d:
149 p[k] = float(
150 self.ranges[k].make_relative(self.base_source[k], d[k]))
152 source = self.base_source.clone(m6=m6, stf=self.get_stf(d), **p)
153 return source
155 def make_dependant(self, xs, pname):
156 cache = self.deps_cache
157 if xs.ndim == 1:
158 return self.make_dependant(xs[num.newaxis, :], pname)[0]
160 if pname not in self.dependant_names:
161 raise KeyError(pname)
163 mt = self.base_source.pyrocko_moment_tensor()
165 sdrs_ref = mt.both_strike_dip_rake()
167 y = num.zeros(xs.shape[0])
168 for i, x in enumerate(xs):
169 k = tuple(x.tolist())
170 if k not in cache:
171 source = self.get_source(x)
172 mt = source.pyrocko_moment_tensor()
173 res = mt.standard_decomposition()
174 sdrs = mt.both_strike_dip_rake()
175 if sdrs_ref:
176 sdrs = mtm.order_like(sdrs, sdrs_ref)
178 cache[k] = mt, res, sdrs
180 mt, res, sdrs = cache[k]
182 if pname == 'rel_moment_iso':
183 ratio_iso, m_iso = res[0][1:3]
184 y[i] = ratio_iso * num.sign(m_iso[0, 0])
185 elif pname == 'rel_moment_clvd':
186 ratio_clvd, m_clvd = res[2][1:3]
187 evals, evecs = mtm.eigh_check(m_clvd)
188 ii = num.argmax(num.abs(evals))
189 y[i] = ratio_clvd * num.sign(evals[ii])
190 else:
191 isdr = {'strike': 0, 'dip': 1, 'rake': 2}[pname[:-1]]
192 y[i] = sdrs[int(pname[-1])-1][isdr]
194 return y
196 def pack_stf(self, stf):
197 return [
198 stf[p.name] for p in self.problem_parameters_stf[self.stf_type]]
200 def pack(self, source):
201 m6 = source.m6
202 mt = source.pyrocko_moment_tensor()
203 rm6 = m6 / mt.scalar_moment()
205 x = num.array([
206 source.time - self.base_source.time,
207 source.north_shift,
208 source.east_shift,
209 source.depth,
210 mt.moment_magnitude(),
211 ] + rm6.tolist() + self.pack_stf(source.stf), dtype=float)
213 return x
215 def random_uniform(self, xbounds, rstate, fixed_magnitude=None):
217 x = num.zeros(self.nparameters)
218 for i in range(self.nparameters):
219 x[i] = rstate.uniform(xbounds[i, 0], xbounds[i, 1])
221 if fixed_magnitude is not None:
222 x[4] = fixed_magnitude
224 if self.mt_random == 'uniform':
225 x[5:11] = mtm.random_m6(x=rstate.random_sample(6))
227 return x.tolist()
229 def preconstrain(self, x):
230 d = self.get_parameter_dict(x)
231 m6 = num.array([d.rmnn, d.rmee, d.rmdd, d.rmne, d.rmnd, d.rmed],
232 dtype=float)
234 m9 = mtm.symmat6(*m6)
235 if self.mt_type == 'deviatoric':
236 trace_m = num.trace(m9)
237 m_iso = num.diag([trace_m / 3., trace_m / 3., trace_m / 3.])
238 m9 -= m_iso
240 elif self.mt_type == 'dc':
241 mt = mtm.MomentTensor(m=m9)
242 m9 = mt.standard_decomposition()[1][2]
244 m0_unscaled = math.sqrt(num.sum(as_arr(m9)**2)) / math.sqrt(2.)
246 m9 /= m0_unscaled
247 m6 = mtm.to6(m9)
248 d.rmnn, d.rmee, d.rmdd, d.rmne, d.rmnd, d.rmed = m6
249 x = self.get_parameter_array(d)
251 source = self.get_source(x)
252 for t in self.waveform_targets:
253 if (self.distance_min > num.asarray(t.distance_to(source))).any():
254 raise Forbidden()
256 return x
258 def get_dependant_bounds(self):
259 out = [
260 (0., 360.),
261 (0., 90.),
262 (-180., 180.),
263 (0., 360.),
264 (0., 90.),
265 (-180., 180.),
266 (-1., 1.),
267 (-1., 1.)]
269 return out
271 @classmethod
272 def get_plot_classes(cls):
273 from .. import plot
274 plots = super(CMTProblem, cls).get_plot_classes()
275 plots.extend([plot.HudsonPlot, plot.MTDecompositionPlot,
276 plot.MTLocationPlot, plot.MTFuzzyPlot])
277 return plots
280__all__ = '''
281 CMTProblem
282 CMTProblemConfig
283'''.split()