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 

4 

5from pyrocko import gf, util, moment_tensor as mtm 

6from pyrocko.guts import String, Float, Dict, StringChoice 

7 

8from grond.meta import Forbidden, expand_template, Parameter, \ 

9 has_get_plot_classes 

10 

11from ..base import Problem, ProblemConfig 

12 

13guts_prefix = 'grond' 

14logger = logging.getLogger('grond.problems.cmt.problem') 

15km = 1e3 

16as_km = dict(scale_factor=km, scale_unit='km') 

17 

18 

19def as_arr(mat_or_arr): 

20 try: 

21 return mat_or_arr.A 

22 except AttributeError: 

23 return mat_or_arr 

24 

25 

26class MTType(StringChoice): 

27 choices = ['full', 'deviatoric', 'dc'] 

28 

29 

30class STFType(StringChoice): 

31 choices = ['HalfSinusoidSTF', 'ResonatorSTF'] 

32 

33 cls = { 

34 'HalfSinusoidSTF': gf.HalfSinusoidSTF, 

35 'ResonatorSTF': gf.ResonatorSTF} 

36 

37 @classmethod 

38 def base_stf(cls, name): 

39 return cls.cls[name]() 

40 

41 

42class CMTProblemConfig(ProblemConfig): 

43 

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 

49 def get_problem(self, event, target_groups, targets): 

50 self.check_deprecations() 

51 

52 if event.depth is None: 

53 event.depth = 0. 

54 

55 base_source = gf.MTSource.from_pyrocko_event(event) 

56 

57 stf = STFType.base_stf(self.stf_type) 

58 stf.duration = event.duration or 0.0 

59 

60 base_source.stf = stf 

61 

62 subs = dict( 

63 event_name=event.name, 

64 event_time=util.time_to_str(event.time)) 

65 

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) 

76 

77 return problem 

78 

79 

80@has_get_plot_classes 

81class CMTProblem(Problem): 

82 

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$')] 

95 

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')]} 

102 

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}$')] 

112 

113 distance_min = Float.T(default=0.0) 

114 mt_type = MTType.T(default='full') 

115 stf_type = STFType.T(default='HalfSinusoidSTF') 

116 

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) 

123 

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]) 

128 

129 return self._base_stf.clone(**d_stf) 

130 

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) 

135 

136 m0 = mtm.magnitude_to_moment(d.magnitude) 

137 m6 = rm6 * m0 

138 

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])) 

144 

145 source = self.base_source.clone(m6=m6, stf=self.get_stf(d), **p) 

146 return source 

147 

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] 

152 

153 if pname not in self.dependant_names: 

154 raise KeyError(pname) 

155 

156 mt = self.base_source.pyrocko_moment_tensor() 

157 

158 sdrs_ref = mt.both_strike_dip_rake() 

159 

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) 

170 

171 cache[k] = mt, res, sdrs 

172 

173 mt, res, sdrs = cache[k] 

174 

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] 

186 

187 return y 

188 

189 def pack_stf(self, stf): 

190 return [ 

191 stf[p.name] for p in self.problem_parameters_stf[self.stf_type]] 

192 

193 def pack(self, source): 

194 m6 = source.m6 

195 mt = source.pyrocko_moment_tensor() 

196 rm6 = m6 / mt.scalar_moment() 

197 

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) 

205 

206 return x 

207 

208 def event_to_x(self, event): 

209 bs = self.base_source 

210 source = bs.__class__.from_pyrocko_event(event) 

211 

212 if source.stf is None: 

213 rstatem = self.get_rstate_manager() 

214 rstate = rstatem.get_rstate(name='event2x') 

215 

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) 

221 

222 return self.source_to_x(source) 

223 

224 def random_uniform(self, xbounds, rstate, fixed_magnitude=None): 

225 

226 x = num.zeros(self.nparameters) 

227 for i in range(self.nparameters): 

228 x[i] = rstate.uniform(xbounds[i, 0], xbounds[i, 1]) 

229 

230 if fixed_magnitude is not None: 

231 x[4] = fixed_magnitude 

232 

233 x[5:11] = mtm.random_m6(x=rstate.random_sample(6)) 

234 

235 return x.tolist() 

236 

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) 

241 

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 

247 

248 elif self.mt_type == 'dc': 

249 mt = mtm.MomentTensor(m=m9) 

250 m9 = mt.standard_decomposition()[1][2] 

251 

252 m0_unscaled = math.sqrt(num.sum(as_arr(m9)**2)) / math.sqrt(2.) 

253 

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) 

258 

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() 

263 

264 return x 

265 

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.)] 

276 

277 return out 

278 

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 

286 

287 

288__all__ = ''' 

289 CMTProblem 

290 CMTProblemConfig 

291'''.split()