Coverage for /usr/local/lib/python3.11/dist-packages/grond/targets/waveform_phase_ratio/plot.py: 21%

1import numpy as num 

2from matplotlib import pyplot as plt 

3from matplotlib import cm 

4 

5from pyrocko.guts import Tuple, Float, String, Int 

6from pyrocko.plot import mpl_init, mpl_margins 

7from pyrocko import gf 

8 

9from grond import core, meta 

10from grond.plot.config import PlotConfig 

11from grond.plot.collection import PlotItem 

12 

13from .target import PhaseRatioTarget 

14 

15guts_prefix = 'grond' 

16 

17 

18def S(text): 

19 return ' '.join(text.split()) 

20 

21 

22class FitsPhaseRatioPlot(PlotConfig): 

23 '''Plot showing the phase ratio fits for the best model.''' 

24 

25 name = 'fits_phase_ratio' 

26 

27 size_cm = Tuple.T( 

28 2, Float.T(), 

29 default=(15., 7.5), 

30 help='Width and length of the figure in [cm]') 

31 

32 misfit_cutoff = Float.T( 

33 optional=True, 

34 help='Plot fits for models up to this misfit value.') 

35 

36 color_parameter = String.T( 

37 default='misfit', 

38 help='Choice of value to color, options: "misfit" (default), ' 

39 '"dist" or any source parameter name.') 

40 

41 istride_ensemble = Int.T( 

42 default=10, 

43 help='Stride value N to choose every Nth model from the solution ' 

44 'ensemble.') 

45 

46 font_size_title = Float.T( 

47 default=10, 

48 help='Font size of title [pt]') 

49 

50 def make(self, environ): 

51 cm = environ.get_plot_collection_manager() 

52 mpl_init(fontsize=self.font_size) 

53 

54 environ.setup_modelling() 

55 ds = environ.get_dataset() 

56 history = environ.get_history(subset='harvest') 

57 

58 scolor = { 

59 'misfit': S(''' 

60 The synthetic markers are colored according to their 

61 respective global (non-bootstrapped) misfit value. Red 

62 indicates better fit, blue worse.'''), 

63 

64 'dist': S(''' 

65 The synthetic markers are colored according to their 

66 Mahalanobis distance from the mean solution.''') 

67 

68 }.get(self.color_parameter, S(''' 

69 The synthetic markers are colored according to source 

70 parameter "%s".''' % self.color_parameter)) 

71 

72 cm.create_group_mpl( 

73 self, self.draw_figures(ds, history), 

74 title=u'Fits Phase Ratios', 

75 section='fits', 

76 feather_icon='activity', 

77 description=u''' 

78Observed (black markers) and synthetic waveform amplitude phase ratio estimates 

79(spectral average ratio; colored markers) at different stations for every Nth 

80model in the bootstrap solution ensemble (N=%i). 

81 

82%s 

83 

84The frequency range used to estimate spectral averages is not shown (see config 

85file). Optimal solutions show good agreement between black and colored markers. 

86''' % (self.istride_ensemble, scolor)) 

87 

88 def draw_figures(self, ds, history): 

89 problem = history.problem 

90 

91 for target in problem.targets: 

92 target.set_dataset(ds) 

93 

94 targets = [ 

95 t for t in problem.targets if isinstance(t, PhaseRatioTarget)] 

96 

97 tpaths = sorted(set(t.path for t in targets)) 

98 

99 for tpath in tpaths: 

100 for (item, fig) in self.draw_figure(ds, history, tpath): 

101 yield item, fig 

102 

103 def draw_figure(self, ds, history, tpath): 

104 problem = history.problem 

105 color_parameter = self.color_parameter 

106 misfit_cutoff = self.misfit_cutoff 

107 fontsize = self.font_size 

108 targets = [ 

109 t for t in problem.targets 

110 if isinstance(t, PhaseRatioTarget) and t.path == tpath] 

111 

112 gms = history.get_sorted_misfits(chain=0)[::-1] 

113 models = history.get_sorted_models(chain=0)[::-1] 

114 

115 if misfit_cutoff is not None: 

116 ibest = gms < misfit_cutoff 

117 gms = gms[ibest] 

118 models = models[ibest] 

119 

120 gms = gms[::self.istride_ensemble] 

121 models = models[::self.istride_ensemble] 

122 

123 nmodels = models.shape[0] 

124 if color_parameter == 'dist': 

125 mx = num.mean(models, axis=0) 

126 cov = num.cov(models.T) 

127 mdists = core.mahalanobis_distance(models, mx, cov) 

128 icolor = meta.ordersort(mdists) 

129 

130 elif color_parameter == 'misfit': 

131 iorder = num.arange(nmodels) 

132 icolor = iorder 

133 

134 elif color_parameter in problem.parameter_names: 

135 ind = problem.name_to_index(color_parameter) 

136 icolor = problem.extract(models, ind) 

137 

138 from matplotlib import colors 

139 cmap = cm.ScalarMappable( 

140 norm=colors.Normalize(vmin=num.min(icolor), vmax=num.max(icolor)), 

141 cmap=plt.get_cmap('coolwarm')) 

142 

143 imodel_to_color = [] 

144 for imodel in range(nmodels): 

145 imodel_to_color.append(cmap.to_rgba(icolor[imodel])) 

146 

147 data = [] 

148 for imodel in range(nmodels): 

149 model = models[imodel, :] 

150 

151 # source = problem.get_source(model) 

152 results = problem.evaluate(model, targets=targets) 

153 

154 for target, result in zip(targets, results): 

155 if isinstance(result, gf.SeismosizerError): 

156 continue 

157 

158 if not isinstance(target, PhaseRatioTarget): 

159 continue 

160 

161 a_obs = result.a_obs 

162 b_obs = result.b_obs 

163 a_syn = result.a_syn 

164 b_syn = result.b_syn 

165 

166 r_obs = a_obs / (a_obs + b_obs) 

167 r_syn = a_syn / (a_syn + b_syn) 

168 

169 data.append(('.'.join(target.codes), imodel, r_obs, r_syn)) 

170 

171 fontsize = self.font_size 

172 

173 item = PlotItem( 

174 name='fig_%s' % tpath) 

175 

176 item.attributes['targets'] = [ 

177 t.string_id() for t in targets] 

178 

179 fig = plt.figure(figsize=self.size_inch) 

180 

181 labelpos = mpl_margins( 

182 fig, nw=1, nh=1, left=7., right=1., bottom=10., top=3, 

183 units=fontsize) 

184 

185 axes = fig.add_subplot(1, 1, 1) 

186 labelpos(axes, 2.5, 2.0) 

187 

188 labels = sorted(set(x[0] for x in data)) 

189 

190 ntargets = len(labels) 

191 string_id_to_itarget = dict((x, i) for (i, x) in enumerate(labels)) 

192 

193 itargets = num.array([string_id_to_itarget[x[0]] for x in data]) 

194 

195 imodels = num.array([x[1] for x in data], dtype=int).T 

196 r_obs, r_syn = num.array([x[2:] for x in data]).T 

197 

198 r_obs_median = num.zeros(ntargets) 

199 for itarget in range(ntargets): 

200 r_obs_median[itarget] = num.median(r_obs[itargets == itarget]) 

201 

202 iorder = meta.ordersort(r_obs_median) 

203 

204 for imodel in range(nmodels): 

205 mask = imodels == imodel 

206 axes.plot( 

207 iorder[itargets[mask]], r_obs[mask], '_', 

208 ms=20., 

209 zorder=-10, 

210 alpha=0.5, 

211 color='black') 

212 axes.plot( 

213 iorder[itargets[mask]], r_syn[mask], '_', 

214 ms=10., 

215 alpha=0.5, 

216 color=imodel_to_color[imodel]) 

217 

218 axes.set_yscale('log') 

219 axes.set_ylabel('Ratio') 

220 

221 axes.set_xticks( 

222 iorder[num.arange(ntargets)]) 

223 axes.set_xticklabels(labels, rotation='vertical') 

224 

225 fig.suptitle(tpath, fontsize=self.font_size_title) 

226 

227 yield item, fig 

228 

229 

230def get_plot_classes(): 

231 return [ 

232 FitsPhaseRatioPlot, 

233 ]