1import numpy as num
2import logging
4from pyrocko import trace
5from pyrocko.guts import (Object, Dict, String, Float, Bool, Int)
7logger = logging.getLogger('grond.synthetic_tests')
9guts_prefix = 'grond'
12class SyntheticWaveformNotAvailable(Exception):
13 pass
16class SyntheticTest(Object):
17 inject_solution = Bool.T(default=False)
18 respect_data_availability = Bool.T(default=False)
19 real_noise_scale = Float.T(default=0.0)
20 white_noise_scale = Float.T(default=0.0)
21 relative_white_noise_scale = Float.T(default=0.0)
22 random_response_scale = Float.T(default=0.0)
23 real_noise_toffset = Float.T(default=-3600.)
24 random_seed = Int.T(optional=True)
25 x = Dict.T(String.T(), Float.T())
27 def __init__(self, **kwargs):
28 Object.__init__(self, **kwargs)
29 self._problem = None
30 self._synthetics = None
32 def set_problem(self, problem):
33 self._problem = problem
34 self._synthetics = None
36 def get_problem(self):
37 if self._problem is None:
38 raise SyntheticWaveformNotAvailable(
39 'SyntheticTest.set_problem() has not been called yet')
41 return self._problem
43 def get_x(self):
44 problem = self.get_problem()
45 if self.x:
46 x = problem.preconstrain(
47 problem.get_parameter_array(self.x))
49 else:
50 x = problem.preconstrain(
51 problem.pack(
52 problem.base_source))
54 return x
56 def get_synthetics(self):
57 problem = self.get_problem()
58 if self._synthetics is None:
59 x = self.get_x()
60 results = problem.forward(x)
61 synthetics = {}
62 for iresult, result in enumerate(results):
63 tr = result.trace.pyrocko_trace()
64 tfade = tr.tmax - tr.tmin
65 tr_orig = tr.copy()
66 tr.extend(tr.tmin - tfade, tr.tmax + tfade)
67 rstate = num.random.RandomState(
68 (self.random_seed or 0) + iresult)
70 if self.random_response_scale != 0:
71 tf = RandomResponse(scale=self.random_response_scale)
72 tf.set_random_state(rstate)
73 tr = tr.transfer(
74 tfade=tfade,
75 transfer_function=tf)
77 if self.white_noise_scale != 0.0:
78 u = rstate.normal(
79 scale=self.white_noise_scale,
80 size=tr.data_len())
82 tr.ydata += u
84 if self.relative_white_noise_scale != 0.0:
85 u = rstate.normal(
86 scale=self.relative_white_noise_scale * num.std(
87 tr_orig.ydata),
88 size=tr.data_len())
90 tr.ydata += u
92 synthetics[result.trace.codes] = tr
94 self._synthetics = synthetics
96 return self._synthetics
98 def get_waveform(self, nslc, tmin, tmax, tfade=0., freqlimits=None):
99 synthetics = self.get_synthetics()
100 if nslc not in synthetics:
101 s = 'no synthetics for %s available' % '.'.join(nslc)
102 logger.warning(s)
103 from grond.dataset import NotFound
104 raise NotFound(s)
106 tr = synthetics[nslc]
107 tr.extend(tmin - tfade * 2.0, tmax + tfade * 2.0)
109 tr = tr.transfer(
110 tfade=tfade,
111 freqlimits=freqlimits)
113 tr.chop(tmin, tmax)
114 return tr
117class RandomResponse(trace.FrequencyResponse):
119 scale = Float.T(default=0.0)
121 def set_random_state(self, rstate):
122 self._rstate = rstate
124 def evaluate(self, freqs):
125 n = freqs.size
126 return 1.0 + (
127 self._rstate.normal(scale=self.scale, size=n) +
128 0.0J * self._rstate.normal(scale=self.scale, size=n))
131__all__ = '''
132SyntheticTest
133SyntheticWaveformNotAvailable
134'''.split()