1import numpy as num
2import logging
4from pyrocko import gf, util
5from pyrocko.guts import String, Float, Dict, StringChoice, Int
7from grond.meta import Forbidden, expand_template, Parameter, \
8 has_get_plot_classes
10from ..base import Problem, ProblemConfig
12guts_prefix = 'grond'
13logger = logging.getLogger('grond.problems.double_sf.problem')
14km = 1e3
15as_km = dict(scale_factor=km, scale_unit='km')
18class DoubleSFProblemConfig(ProblemConfig):
20 ranges = Dict.T(String.T(), gf.Range.T())
21 distance_min = Float.T(default=0.0)
22 nthreads = Int.T(default=1)
23 force_directions = StringChoice.T(
24 choices=('off', 'unidirectional', 'counterdirectional'),
25 default='off')
27 def get_problem(self, event, target_groups, targets):
28 if event.depth is None:
29 event.depth = 0.
31 base_source = gf.DoubleSFSource.from_pyrocko_event(event)
33 base_source.stf1 = gf.HalfSinusoidSTF(duration=event.duration or 0.0)
34 base_source.stf2 = gf.HalfSinusoidSTF(duration=event.duration or 0.0)
36 subs = dict(
37 event_name=event.name,
38 event_time=util.time_to_str(event.time))
40 problem = DoubleSFProblem(
41 name=expand_template(self.name_template, subs),
42 base_source=base_source,
43 target_groups=target_groups,
44 targets=targets,
45 ranges=self.ranges,
46 distance_min=self.distance_min,
47 norm_exponent=self.norm_exponent,
48 nthreads=self.nthreads,
49 force_directions=self.force_directions)
51 return problem
54@has_get_plot_classes
55class DoubleSFProblem(Problem):
57 problem_parameters = [
58 Parameter('time', 's', label='Time'),
59 Parameter('north_shift', 'm', label='Northing', **as_km),
60 Parameter('east_shift', 'm', label='Easting', **as_km),
61 Parameter('depth', 'm', label='Depth', **as_km),
62 Parameter('force', 'N', label='$||F||$'),
63 Parameter('rfn1', '', label='$rF_{n1}$'),
64 Parameter('rfe1', '', label='$rF_{e1}$'),
65 Parameter('rfd1', '', label='$rF_{d1}$'),
66 Parameter('rfn2', '', label='$rF_{n2}$'),
67 Parameter('rfe2', '', label='$rF_{e2}$'),
68 Parameter('rfd2', '', label='$rF_{d2}$'),
69 Parameter('delta_time', 's', label='$\\Delta$ Time'),
70 Parameter('delta_depth', 'm', label='$\\Delta$ Depth'),
71 Parameter('azimuth', 'deg', label='Azimuth'),
72 Parameter('distance', 'm', label='Distance'),
73 Parameter('mix', label='Mix'),
74 Parameter('duration1', 's', label='Duration 1'),
75 Parameter('duration2', 's', label='Duration 2')]
77 dependants = [
78 Parameter('fn1', 'N', label='$F_{n1}$'),
79 Parameter('fe1', 'N', label='$F_{e1}$'),
80 Parameter('fd1', 'N', label='$F_{d1}$'),
81 Parameter('fn2', 'N', label='$F_{n2}$'),
82 Parameter('fe2', 'N', label='$F_{e2}$'),
83 Parameter('fd2', 'N', label='$F_{d2}$')]
85 distance_min = Float.T(default=0.0)
86 force_directions = String.T()
88 def __init__(self, **kwargs):
89 Problem.__init__(self, **kwargs)
90 self.deps_cache = {}
92 def get_source(self, x):
93 d = self.get_parameter_dict(x)
95 p = {}
96 for k in self.base_source.keys():
97 if k in d:
98 p[k] = float(
99 self.ranges[k].make_relative(self.base_source[k], d[k]))
101 stf1 = gf.HalfSinusoidSTF(duration=float(d.duration1))
102 stf2 = gf.HalfSinusoidSTF(duration=float(d.duration2))
104 source = self.base_source.clone(stf1=stf1, stf2=stf2, **p)
105 return source
107 def make_dependant(self, xs, pname):
108 cache = self.deps_cache
109 if xs.ndim == 1:
110 return self.make_dependant(xs[num.newaxis, :], pname)[0]
112 if pname not in self.dependant_names:
113 raise KeyError(pname)
115 y = num.zeros(xs.shape[0])
116 for i, x in enumerate(xs):
117 k = tuple(x.tolist())
118 if k not in cache:
119 source = self.get_source(x)
120 cache[k] = source
122 source = cache[k]
124 y[i] = getattr(source, pname)
126 return y
128 def pack(self, source):
129 arr = self.get_parameter_array(source)
130 for ip, p in enumerate(self.parameters):
131 if p.name == 'time':
132 arr[ip] -= self.base_source.time
133 if p.name == 'duration1':
134 arr[ip] = source.stf1.duration if source.stf1 else 0.0
135 if p.name == 'duration2':
136 arr[ip] = source.stf2.duration if source.stf2 else 0.0
137 return arr
139 def random_uniform(self, xbounds, rstate, fixed_magnitude=None):
141 x = num.zeros(self.nparameters)
142 for i in range(self.nparameters):
143 x[i] = rstate.uniform(xbounds[i, 0], xbounds[i, 1])
145 return x.tolist()
147 def preconstrain(self, x):
148 source = self.get_source(x)
149 if any(self.distance_min > source.distance_to(t)
150 for t in self.targets):
151 raise Forbidden()
153 if self.force_directions == 'unidirectional':
154 idx_rfn2 = self.get_parameter_index('rfn2')
155 idx_rfe2 = self.get_parameter_index('rfe2')
156 idx_rfd2 = self.get_parameter_index('rfd2')
158 x[idx_rfn2] = source.rfn1
159 x[idx_rfe2] = source.rfe1
160 x[idx_rfd2] = source.rfd1
162 elif self.force_directions == 'counterdirectional':
163 idx_rfn2 = self.get_parameter_index('rfn2')
164 idx_rfe2 = self.get_parameter_index('rfe2')
165 idx_rfd2 = self.get_parameter_index('rfd2')
167 x[idx_rfn2] = -source.rfn1
168 x[idx_rfe2] = -source.rfe1
169 x[idx_rfd2] = -source.rfd1
171 return num.array(x, dtype=num.float)
173 def get_dependant_bounds(self):
174 range_force = self.ranges['force']
175 out = [
176 (-range_force.stop, range_force.stop),
177 (-range_force.stop, range_force.stop),
178 (-range_force.stop, range_force.stop),
179 (-range_force.stop, range_force.stop),
180 (-range_force.stop, range_force.stop),
181 (-range_force.stop, range_force.stop)]
183 return out
185 @classmethod
186 def get_plot_classes(cls):
187 from . import plot
188 from ..singleforce import plot as sfplot
189 plots = super(DoubleSFProblem, cls).get_plot_classes()
190 plots.extend([
191 sfplot.SFLocationPlot,
192 plot.SFForcePlot,
193 plot.DoubleSFDecompositionPlot])
194 return plots
197class IndependentDoubleSFProblemConfig(ProblemConfig):
199 ranges = Dict.T(String.T(), gf.Range.T())
200 distance_min = Float.T(default=0.0)
201 nthreads = Int.T(default=1)
202 force_directions = StringChoice.T(
203 choices=('off', 'unidirectional', 'counterdirectional'),
204 default='off')
206 def get_problem(self, event, target_groups, targets):
207 if event.depth is None:
208 event.depth = 0.
210 source = gf.SFSource.from_pyrocko_event(event)
211 source.stf = gf.HalfSinusoidSTF(duration=event.duration or 0.0)
213 base_source = gf.CombiSFSource(
214 name=source.name,
215 subsources=[
216 source.clone(name=''),
217 source.clone(name='')])
219 source.lat, source.lon = event.effective_latlon
221 subs = dict(
222 event_name=event.name,
223 event_time=util.time_to_str(event.time))
225 problem = IndependentDoubleSFProblem(
226 name=expand_template(self.name_template, subs),
227 base_source=base_source,
228 target_groups=target_groups,
229 targets=targets,
230 ranges=self.ranges,
231 distance_min=self.distance_min,
232 norm_exponent=self.norm_exponent,
233 nthreads=self.nthreads,
234 force_directions=self.force_directions)
236 return problem
239@has_get_plot_classes
240class IndependentDoubleSFProblem(Problem):
242 problem_parameters = [
243 Parameter('time1', 's', label='Time'),
244 Parameter('north_shift1', 'm', label='Northing', **as_km),
245 Parameter('east_shift1', 'm', label='Easting', **as_km),
246 Parameter('depth1', 'm', label='Depth', **as_km),
247 Parameter('time2', 's', label='Time'),
248 Parameter('north_shift2', 'm', label='Northing', **as_km),
249 Parameter('east_shift2', 'm', label='Easting', **as_km),
250 Parameter('depth2', 'm', label='Depth', **as_km),
251 Parameter('fn1', 'N', label='$F_{n1}$'),
252 Parameter('fe1', 'N', label='$F_{e1}$'),
253 Parameter('fd1', 'N', label='$F_{d1}$'),
254 Parameter('fn2', 'N', label='$F_{n2}$'),
255 Parameter('fe2', 'N', label='$F_{e2}$'),
256 Parameter('fd2', 'N', label='$F_{d2}$'),
257 Parameter('duration1', 's', label='Duration 1'),
258 Parameter('duration2', 's', label='Duration 2')]
260 dependants = [
261 Parameter('force1', 'N', label='$||F_1||$'),
262 Parameter('force2', 'N', label='$||F_2||$')]
264 distance_min = Float.T(default=0.0)
265 force_directions = String.T()
267 def __init__(self, **kwargs):
268 Problem.__init__(self, **kwargs)
269 self.deps_cache = {}
271 def get_source(self, x):
272 d = self.get_parameter_dict(x)
274 sources = []
276 for i, subsource in enumerate(self.base_source.subsources):
277 p = {}
279 for k in subsource.keys():
280 key = '%s%d' % (k, i+1)
281 if key in d:
282 p[k] = float(
283 self.ranges[key].make_relative(subsource[k], d[key]))
285 sources.append(self.base_source.subsources[i].clone(**p))
287 sources[0].stf = gf.HalfSinusoidSTF(duration=float(d.duration1))
288 sources[1].stf = gf.HalfSinusoidSTF(duration=float(d.duration2))
290 return self.base_source.clone(subsources=sources)
292 def make_dependant(self, xs, pname):
293 cache = self.deps_cache
294 if xs.ndim == 1:
295 return self.make_dependant(xs[num.newaxis, :], pname)[0]
297 if pname not in self.dependant_names:
298 raise KeyError(pname)
300 y = num.zeros(xs.shape[0])
301 for i, x in enumerate(xs):
302 k = tuple(x.tolist())
303 if k not in cache:
304 source = self.get_source(x)
305 cache[k] = source
307 source = cache[k]
309 y[i] = getattr(source.subsources[int(pname[-1]) - 1], pname[:-1])
311 return y
313 def pack(self, source):
314 arr = self.get_parameter_array(source)
315 subsrcs = source.subsources
316 for ip, p in enumerate(self.parameters):
317 # if p.name == 'time':
318 # arr[ip] -= self.base_source.time
319 if p.name == 'duration1':
320 arr[ip] = subsrcs[0].stf.duration if subsrcs[0].stf else 0.0
321 if p.name == 'duration2':
322 arr[ip] = subsrcs[1].stf.duration if subsrcs[1].stf else 0.0
323 return arr
325 def random_uniform(self, xbounds, rstate, fixed_magnitude=None):
327 x = num.zeros(self.nparameters)
328 for i in range(self.nparameters):
329 x[i] = rstate.uniform(xbounds[i, 0], xbounds[i, 1])
331 return x.tolist()
333 def preconstrain(self, x):
334 source = self.get_source(x)
335 if any(self.distance_min > source.distance_to(t)
336 for t in self.targets):
337 raise Forbidden()
339 if self.force_directions == 'unidirectional':
340 force1 = source.subsources[0].force
341 force2 = source.subsources[1].force
343 ratio = force2 / force1
345 idx_fn2 = self.get_parameter_index('fn2')
346 idx_fe2 = self.get_parameter_index('fe2')
347 idx_fd2 = self.get_parameter_index('fd2')
349 x[idx_fn2] = source.subsources[0].fn * ratio
350 x[idx_fe2] = source.subsources[0].fe * ratio
351 x[idx_fd2] = source.subsources[0].fd * ratio
353 elif self.force_directions == 'counterdirectional':
354 force1 = source.subsources[0].force
355 force2 = source.subsources[1].force
357 ratio = force2 / force1
359 idx_fn2 = self.get_parameter_index('fn2')
360 idx_fe2 = self.get_parameter_index('fe2')
361 idx_fd2 = self.get_parameter_index('fd2')
363 x[idx_fn2] = -source.subsources[0].fn * ratio
364 x[idx_fe2] = -source.subsources[0].fe * ratio
365 x[idx_fd2] = -source.subsources[0].fd * ratio
367 return num.array(x, dtype=num.float)
369 def get_dependant_bounds(self):
370 range_start = num.min([
371 (self.ranges['f{}1'.format(f)].start,
372 self.ranges['f{}2'.format(f)].start) for f in 'n e d'.split()],
373 axis=0)
375 range_stop = num.max([
376 (self.ranges['f{}1'.format(f)].stop,
377 self.ranges['f{}2'.format(f)].stop) for f in 'n e d'.split()],
378 axis=0)
380 force_range = (
381 -num.linalg.norm(range_start),
382 num.linalg.norm(range_stop))
384 out = [force_range, force_range]
386 return out
388 @classmethod
389 def get_plot_classes(cls):
390 from . import plot
391 plots = super(IndependentDoubleSFProblem, cls).get_plot_classes()
392 plots.extend([
393 plot.SFForcePlot,
394 plot.DoubleSFDecompositionPlot])
395 return plots
398__all__ = '''
399 DoubleSFProblem
400 DoubleSFProblemConfig
401 IndependentDoubleSFProblem
402 IndependentDoubleSFProblemConfig
403'''.split()