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, Int
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')
48 nthreads = Int.T(default=1)
50 def get_problem(self, event, target_groups, targets):
51 if event.depth is None:
52 event.depth = 0.
54 base_source = gf.MTSource.from_pyrocko_event(event)
56 stf = STFType.base_stf(self.stf_type)
57 stf.duration = event.duration or 0.0
59 base_source.stf = stf
61 subs = dict(
62 event_name=event.name,
63 event_time=util.time_to_str(event.time))
65 problem = CMTProblem(
66 name=expand_template(self.name_template, subs),
67 base_source=base_source,
68 target_groups=target_groups,
69 targets=targets,
70 ranges=self.ranges,
71 distance_min=self.distance_min,
72 mt_type=self.mt_type,
73 stf_type=self.stf_type,
74 norm_exponent=self.norm_exponent,
75 nthreads=self.nthreads)
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 random_uniform(self, xbounds, rstate, fixed_magnitude=None):
210 x = num.zeros(self.nparameters)
211 for i in range(self.nparameters):
212 x[i] = rstate.uniform(xbounds[i, 0], xbounds[i, 1])
214 if fixed_magnitude is not None:
215 x[4] = fixed_magnitude
217 x[5:11] = mtm.random_m6(x=rstate.random_sample(6))
219 return x.tolist()
221 def preconstrain(self, x):
222 d = self.get_parameter_dict(x)
223 m6 = num.array([d.rmnn, d.rmee, d.rmdd, d.rmne, d.rmnd, d.rmed],
224 dtype=float)
226 m9 = mtm.symmat6(*m6)
227 if self.mt_type == 'deviatoric':
228 trace_m = num.trace(m9)
229 m_iso = num.diag([trace_m / 3., trace_m / 3., trace_m / 3.])
230 m9 -= m_iso
232 elif self.mt_type == 'dc':
233 mt = mtm.MomentTensor(m=m9)
234 m9 = mt.standard_decomposition()[1][2]
236 m0_unscaled = math.sqrt(num.sum(as_arr(m9)**2)) / math.sqrt(2.)
238 m9 /= m0_unscaled
239 m6 = mtm.to6(m9)
240 d.rmnn, d.rmee, d.rmdd, d.rmne, d.rmnd, d.rmed = m6
241 x = self.get_parameter_array(d)
243 source = self.get_source(x)
244 for t in self.waveform_targets:
245 if (self.distance_min > num.asarray(t.distance_to(source))).any():
246 raise Forbidden()
248 return x
250 def get_dependant_bounds(self):
251 out = [
252 (0., 360.),
253 (0., 90.),
254 (-180., 180.),
255 (0., 360.),
256 (0., 90.),
257 (-180., 180.),
258 (-1., 1.),
259 (-1., 1.)]
261 return out
263 @classmethod
264 def get_plot_classes(cls):
265 from .. import plot
266 plots = super(CMTProblem, cls).get_plot_classes()
267 plots.extend([plot.HudsonPlot, plot.MTDecompositionPlot,
268 plot.MTLocationPlot, plot.MTFuzzyPlot])
269 return plots
272__all__ = '''
273 CMTProblem
274 CMTProblemConfig
275'''.split()