Coverage for /usr/local/lib/python3.11/dist-packages/grond/problems/plot.py: 14%

1import math 

2import logging 

3import random 

4 

5import numpy as num 

6from matplotlib import cm, patches, colors as mcolors 

7 

8from pyrocko.guts import Tuple, Float, Int, String, List, Bool, StringChoice 

9 

10from pyrocko.plot import mpl_margins, mpl_color, mpl_init, mpl_graph_color 

11from pyrocko.plot import beachball, hudson 

12 

13from grond.plot.config import PlotConfig 

14from grond.plot.section import SectionPlotConfig, SectionPlot 

15from grond.plot.collection import PlotItem 

16from grond import meta, core, stats 

17from matplotlib import pyplot as plt 

18 

19logger = logging.getLogger('grond.problem.plot') 

20 

21guts_prefix = 'grond' 

22 

23 

24def cluster_label(icluster, perc): 

25 if icluster == -1: 

26 return 'Unclust. (%.0f%%)' % perc 

27 else: 

28 return 'Cluster %i (%.0f%%)' % (icluster, perc) 

29 

30 

31def cluster_color(icluster): 

32 if icluster == -1: 

33 return mpl_color('aluminium3') 

34 else: 

35 return mpl_graph_color(icluster) 

36 

37 

38def fixlim(lo, hi): 

39 if lo == hi: 

40 return lo - 1.0, hi + 1.0 

41 else: 

42 return lo, hi 

43 

44 

45def eigh_sorted(mat): 

46 evals, evecs = num.linalg.eigh(mat) 

47 iorder = num.argsort(evals) 

48 return evals[iorder], evecs[:, iorder] 

49 

50 

51class JointparPlot(PlotConfig): 

52 ''' 

53 Source problem parameter's tradeoff plots. 

54 ''' 

55 

56 name = 'jointpar' 

57 size_cm = Tuple.T(2, Float.T(), default=(20., 20.)) 

58 misfit_cutoff = Float.T(optional=True) 

59 ibootstrap = Int.T(optional=True) 

60 color_parameter = String.T(default='misfit') 

61 exclude = List.T(String.T()) 

62 include = List.T(String.T()) 

63 show_ellipses = Bool.T(default=False) 

64 nsubplots = Int.T(default=6) 

65 show_ticks = Bool.T(default=False) 

66 show_reference = Bool.T(default=True) 

67 

68 def make(self, environ): 

69 cm = environ.get_plot_collection_manager() 

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

71 optimiser = environ.get_optimiser() 

72 sref = 'Dark gray boxes mark reference solution.' \ 

73 if self.show_reference else '' 

74 

75 mpl_init(fontsize=self.font_size) 

76 cm.create_group_mpl( 

77 self, 

78 self.draw_figures(history, optimiser), 

79 title=u'Jointpar Plot', 

80 section='solution', 

81 feather_icon='crosshair', 

82 description=u''' 

83Source problem parameter's scatter plots, to evaluate the resolution of source 

84parameters and possible trade-offs between pairs of model parameters. 

85 

86A subset of model solutions (from harvest) is shown in two dimensions for all 

87possible parameter pairs as points. The point color indicates the misfit for 

88the model solution with cold colors (blue) for high misfit models and warm 

89colors (red) for low misfit models. The plot ranges are defined by the given 

90parameter bounds and shows the model space of the optimsation. %s''' % sref) 

91 

92 def draw_figures(self, history, optimiser): 

93 

94 color_parameter = self.color_parameter 

95 exclude = self.exclude 

96 include = self.include 

97 nsubplots = self.nsubplots 

98 figsize = self.size_inch 

99 ibootstrap = 0 if self.ibootstrap is None else self.ibootstrap 

100 misfit_cutoff = self.misfit_cutoff 

101 show_ellipses = self.show_ellipses 

102 msize = 1.5 

103 cmap = 'coolwarm' 

104 

105 problem = history.problem 

106 if not problem: 

107 return [] 

108 

109 models = history.models 

110 

111 exclude = list(exclude) 

112 bounds = problem.get_combined_bounds() 

113 for ipar in range(problem.ncombined): 

114 par = problem.combined[ipar] 

115 lo, hi = bounds[ipar] 

116 if lo == hi: 

117 exclude.append(par.name) 

118 

119 xref = problem.get_reference_model() 

120 

121 isort = history.get_sorted_misfits_idx(chain=ibootstrap)[::-1] 

122 models = history.get_sorted_models(chain=ibootstrap)[::-1] 

123 nmodels = history.nmodels 

124 

125 gms = history.get_sorted_misfits(chain=ibootstrap)[::-1] 

126 if misfit_cutoff is not None: 

127 ibest = gms < misfit_cutoff 

128 gms = gms[ibest] 

129 models = models[ibest] 

130 

131 kwargs = {} 

132 

133 if color_parameter == 'dist': 

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

135 cov = num.cov(models.T) 

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

137 icolor = meta.ordersort(mdists) 

138 

139 elif color_parameter == 'misfit': 

140 iorder = num.arange(nmodels) 

141 icolor = iorder 

142 

143 elif color_parameter in problem.parameter_names: 

144 ind = problem.name_to_index(color_parameter) 

145 icolor = problem.extract(models, ind) 

146 

147 elif color_parameter in history.attribute_names: 

148 icolor = history.get_attribute(color_parameter)[isort] 

149 icolor_need = num.unique(icolor) 

150 

151 colors = [] 

152 for i in range(icolor_need[-1]+1): 

153 colors.append(mpl_graph_color(i)) 

154 

155 cmap = mcolors.ListedColormap(colors) 

156 cmap.set_under(mpl_color('aluminium3')) 

157 kwargs.update(dict(vmin=0, vmax=icolor_need[-1])) 

158 else: 

159 raise meta.GrondError( 

160 'Invalid color_parameter: %s' % color_parameter) 

161 

162 smap = {} 

163 iselected = 0 

164 for ipar in range(problem.ncombined): 

165 par = problem.combined[ipar] 

166 if exclude and par.name in exclude or \ 

167 include and par.name not in include: 

168 continue 

169 

170 smap[iselected] = ipar 

171 iselected += 1 

172 

173 nselected = iselected 

174 

175 if nselected < 2: 

176 logger.warning( 

177 'Cannot draw joinpar figures with less than two parameters ' 

178 'selected.') 

179 return [] 

180 

181 nfig = (nselected - 2) // nsubplots + 1 

182 

183 figs = [] 

184 for ifig in range(nfig): 

185 figs_row = [] 

186 for jfig in range(nfig): 

187 if ifig >= jfig: 

188 item = PlotItem(name='fig_%i_%i' % (ifig, jfig)) 

189 item.attributes['parameters'] = [] 

190 figs_row.append((item, plt.figure(figsize=figsize))) 

191 else: 

192 figs_row.append(None) 

193 

194 figs.append(figs_row) 

195 

196 for iselected in range(nselected): 

197 ipar = smap[iselected] 

198 ypar = problem.combined[ipar] 

199 for jselected in range(iselected): 

200 jpar = smap[jselected] 

201 xpar = problem.combined[jpar] 

202 

203 ixg = (iselected - 1) 

204 iyg = jselected 

205 

206 ix = ixg % nsubplots 

207 iy = iyg % nsubplots 

208 

209 ifig = ixg // nsubplots 

210 jfig = iyg // nsubplots 

211 

212 aind = (nsubplots, nsubplots, (ix * nsubplots) + iy + 1) 

213 

214 item, fig = figs[ifig][jfig] 

215 

216 tlist = item.attributes['parameters'] 

217 if xpar.name not in tlist: 

218 tlist.append(xpar.name) 

219 if ypar.name not in tlist: 

220 tlist.append(ypar.name) 

221 

222 axes = fig.add_subplot(*aind) 

223 

224 axes.axvline(0., color=mpl_color('aluminium3'), lw=0.5) 

225 axes.axhline(0., color=mpl_color('aluminium3'), lw=0.5) 

226 for spine in axes.spines.values(): 

227 spine.set_edgecolor(mpl_color('aluminium5')) 

228 spine.set_linewidth(0.5) 

229 

230 xmin, xmax = fixlim(*xpar.scaled(bounds[jpar])) 

231 ymin, ymax = fixlim(*ypar.scaled(bounds[ipar])) 

232 

233 if ix == 0 or jselected + 1 == iselected: 

234 for (xpos, xoff, x) in [ 

235 (0.0, 10., xmin), 

236 (1.0, -10., xmax)]: 

237 

238 axes.annotate( 

239 '%.3g%s' % (x, xpar.get_unit_suffix()), 

240 xy=(xpos, 1.05), 

241 xycoords='axes fraction', 

242 xytext=(xoff, 5.), 

243 textcoords='offset points', 

244 verticalalignment='bottom', 

245 horizontalalignment='left', 

246 rotation=45.) 

247 

248 if iy == nsubplots - 1 or jselected + 1 == iselected: 

249 for (ypos, yoff, y) in [ 

250 (0., 10., ymin), 

251 (1.0, -10., ymax)]: 

252 

253 axes.annotate( 

254 '%.3g%s' % (y, ypar.get_unit_suffix()), 

255 xy=(1.0, ypos), 

256 xycoords='axes fraction', 

257 xytext=(5., yoff), 

258 textcoords='offset points', 

259 verticalalignment='bottom', 

260 horizontalalignment='left', 

261 rotation=45.) 

262 

263 axes.set_xlim(xmin, xmax) 

264 axes.set_ylim(ymin, ymax) 

265 

266 if not self.show_ticks: 

267 axes.get_xaxis().set_ticks([]) 

268 axes.get_yaxis().set_ticks([]) 

269 else: 

270 axes.tick_params(length=4, which='both') 

271 axes.get_yaxis().set_ticklabels([]) 

272 axes.get_xaxis().set_ticklabels([]) 

273 

274 if iselected == nselected - 1 or ix == nsubplots - 1: 

275 axes.annotate( 

276 xpar.get_label(with_unit=False), 

277 xy=(0.5, -0.05), 

278 xycoords='axes fraction', 

279 verticalalignment='top', 

280 horizontalalignment='right', 

281 rotation=45.) 

282 

283 if iy == 0: 

284 axes.annotate( 

285 ypar.get_label(with_unit=False), 

286 xy=(-0.05, 0.5), 

287 xycoords='axes fraction', 

288 verticalalignment='top', 

289 horizontalalignment='right', 

290 rotation=45.) 

291 

292 fx = problem.extract(models, jpar) 

293 fy = problem.extract(models, ipar) 

294 

295 axes.scatter( 

296 xpar.scaled(fx), 

297 ypar.scaled(fy), 

298 c=icolor, 

299 s=msize, alpha=0.5, cmap=cmap, edgecolors='none', **kwargs) 

300 

301 if show_ellipses: 

302 cov = num.cov((xpar.scaled(fx), ypar.scaled(fy))) 

303 evals, evecs = eigh_sorted(cov) 

304 evals = num.sqrt(evals) 

305 ell = patches.Ellipse( 

306 xy=( 

307 num.mean(xpar.scaled(fx)), 

308 num.mean(ypar.scaled(fy))), 

309 width=evals[0] * 2, 

310 height=evals[1] * 2, 

311 angle=num.rad2deg( 

312 num.arctan2(evecs[1][0], evecs[0][0]))) 

313 

314 ell.set_facecolor('none') 

315 axes.add_artist(ell) 

316 

317 if self.show_reference: 

318 fx = problem.extract(xref, jpar) 

319 fy = problem.extract(xref, ipar) 

320 

321 ref_color = mpl_color('aluminium6') 

322 ref_color_light = 'none' 

323 axes.plot( 

324 xpar.scaled(fx), ypar.scaled(fy), 's', 

325 mew=1.5, ms=5, mfc=ref_color_light, mec=ref_color) 

326 

327 figs_flat = [] 

328 for figs_row in figs: 

329 figs_flat.extend( 

330 item_fig for item_fig in figs_row if item_fig is not None) 

331 

332 return figs_flat 

333 

334 

335class HistogramPlot(PlotConfig): 

336 ''' 

337 Histograms or Gaussian kernel densities (default) of all parameters 

338 (marginal distributions of model parameters). 

339 

340 The histograms (by default shown as Gaussian kernel densities) show (red 

341 curved solid line) the distributions of the parameters (marginals) along 

342 with some characteristics: The red solid vertical line gives the median of 

343 the distribution and the dashed red vertical line the mean value. Dark gray 

344 vertical lines show reference values (given in the event.txt file). The 

345 overlapping red-shaded areas show the 68% confidence intervals (innermost 

346 area), the 90% confidence intervals (middle area) and the minimum and 

347 maximum values (widest area). The plot ranges are defined by the given 

348 parameter bounds and show the model space. Well resolved model parameters 

349 show peaked distributions. 

350 ''' 

351 

352 name = 'histogram' 

353 size_cm = Tuple.T(2, Float.T(), default=(12.5, 7.5)) 

354 exclude = List.T(String.T()) 

355 include = List.T(String.T()) 

356 method = StringChoice.T( 

357 choices=['gaussian_kde', 'histogram'], 

358 default='gaussian_kde') 

359 show_reference = Bool.T(default=True) 

360 

361 def make(self, environ): 

362 cm = environ.get_plot_collection_manager() 

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

364 

365 mpl_init(fontsize=self.font_size) 

366 cm.create_group_mpl( 

367 self, 

368 self.draw_figures(history), 

369 title=u'Histogram', 

370 section='solution', 

371 feather_icon='bar-chart-2', 

372 description=u''' 

373Distribution of the problem's parameters. 

374 

375The histograms are shown either as Gaussian kernel densities (red curved solid 

376line) or as bar plots the distributions of the parameters (marginals) along 

377with some characteristics: 

378 

379The red solid vertical line gives the median of the distribution and the dashed 

380red vertical line the mean value. Dark gray vertical lines show reference 

381parameter values if given in the event.txt file. The overlapping red-shaded 

382areas show the 68% confidence intervals (innermost area), the 90% confidence 

383intervals (middle area) and the minimum and maximum values (widest area). The 

384plot ranges are defined by the given parameter bounds and show the model 

385space. 

386''') 

387 

388 def draw_figures(self, history): 

389 

390 import scipy.stats 

391 from grond.core import make_stats 

392 

393 exclude = self.exclude 

394 include = self.include 

395 figsize = self.size_inch 

396 fontsize = self.font_size 

397 method = self.method 

398 ref_color = mpl_color('aluminium6') 

399 stats_color = mpl_color('scarletred2') 

400 bar_color = mpl_color('scarletred1') 

401 stats_color3 = mpl_color('scarletred3') 

402 

403 problem = history.problem 

404 

405 models = history.models 

406 

407 bounds = problem.get_combined_bounds() 

408 exclude = list(exclude) 

409 for ipar in range(problem.ncombined): 

410 par = problem.combined[ipar] 

411 vmin, vmax = bounds[ipar] 

412 if vmin == vmax: 

413 exclude.append(par.name) 

414 

415 xref = problem.get_reference_model() 

416 

417 smap = {} 

418 iselected = 0 

419 for ipar in range(problem.ncombined): 

420 par = problem.combined[ipar] 

421 if exclude and par.name in exclude or \ 

422 include and par.name not in include: 

423 continue 

424 

425 smap[iselected] = ipar 

426 iselected += 1 

427 

428 nselected = iselected 

429 del iselected 

430 

431 pnames = [ 

432 problem.combined[smap[iselected]].name 

433 for iselected in range(nselected)] 

434 

435 rstats = make_stats(problem, models, 

436 history.get_primary_chain_misfits(), 

437 pnames=pnames) 

438 

439 for iselected in range(nselected): 

440 ipar = smap[iselected] 

441 par = problem.combined[ipar] 

442 vs = problem.extract(models, ipar) 

443 vmin, vmax = bounds[ipar] 

444 

445 fig = plt.figure(figsize=figsize) 

446 labelpos = mpl_margins( 

447 fig, nw=1, nh=1, w=7., bottom=5., top=1, units=fontsize) 

448 

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

450 labelpos(axes, 2.5, 2.0) 

451 axes.set_xlabel(par.get_label()) 

452 axes.set_ylabel('PDF') 

453 axes.set_xlim(*fixlim(*par.scaled((vmin, vmax)))) 

454 

455 if method == 'gaussian_kde': 

456 try: 

457 kde = scipy.stats.gaussian_kde(vs) 

458 except Exception: 

459 logger.warning( 

460 'Cannot create plot histogram with gaussian_kde: ' 

461 'possibly all samples have the same value.') 

462 continue 

463 

464 vps = num.linspace(vmin, vmax, 600) 

465 pps = kde(vps) 

466 

467 axes.plot( 

468 par.scaled(vps), par.inv_scaled(pps), color=stats_color) 

469 

470 elif method == 'histogram': 

471 pps, edges = num.histogram( 

472 vs, 

473 bins=num.linspace(vmin, vmax, num=40), 

474 density=True) 

475 vps = 0.5 * (edges[:-1] + edges[1:]) 

476 

477 axes.bar(par.scaled(vps), par.inv_scaled(pps), 

478 par.scaled(2.*(vps - edges[:-1])), 

479 color=bar_color) 

480 

481 pstats = rstats.parameter_stats_list[iselected] 

482 

483 axes.axvspan( 

484 par.scaled(pstats.minimum), 

485 par.scaled(pstats.maximum), 

486 color=stats_color, alpha=0.1) 

487 axes.axvspan( 

488 par.scaled(pstats.percentile16), 

489 par.scaled(pstats.percentile84), 

490 color=stats_color, alpha=0.1) 

491 axes.axvspan( 

492 par.scaled(pstats.percentile5), 

493 par.scaled(pstats.percentile95), 

494 color=stats_color, alpha=0.1) 

495 

496 axes.axvline( 

497 par.scaled(pstats.median), 

498 color=stats_color3, alpha=0.5) 

499 axes.axvline( 

500 par.scaled(pstats.mean), 

501 color=stats_color3, ls=':', alpha=0.5) 

502 

503 if self.show_reference: 

504 axes.axvline( 

505 par.scaled(problem.extract(xref, ipar)), 

506 color=ref_color) 

507 

508 item = PlotItem(name=par.name) 

509 item.attributes['parameters'] = [par.name] 

510 yield item, fig 

511 

512 

513class MTDecompositionPlot(PlotConfig): 

514 ''' 

515 Moment tensor decomposition plot. 

516 ''' 

517 

518 name = 'mt_decomposition' 

519 size_cm = Tuple.T(2, Float.T(), default=(15., 5.)) 

520 cluster_attribute = meta.StringID.T( 

521 optional=True, 

522 help='name of attribute to use as cluster IDs') 

523 show_reference = Bool.T(default=True) 

524 

525 def make(self, environ): 

526 cm = environ.get_plot_collection_manager() 

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

528 mpl_init(fontsize=self.font_size) 

529 cm.create_group_mpl( 

530 self, 

531 self.draw_figures(history), 

532 title=u'MT Decomposition', 

533 section='solution', 

534 feather_icon='sun', 

535 description=u''' 

536Moment tensor decomposition of the best-fitting solution into isotropic, 

537deviatoric and best double couple components. 

538 

539Shown are the ensemble best, the ensemble mean%s and, if available, a reference 

540mechanism. The symbol size indicates the relative strength of the components. 

541The inversion result is consistent and stable if ensemble mean and ensemble 

542best have similar symbol size and patterns. 

543''' % (', cluster results' if self.cluster_attribute else '')) 

544 

545 def draw_figures(self, history): 

546 

547 fontsize = self.font_size 

548 

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

550 axes = fig.add_subplot(1, 1, 1, aspect=1.0) 

551 fig.subplots_adjust(left=0., right=1., bottom=0., top=1.) 

552 

553 problem = history.problem 

554 models = history.models 

555 

556 if models.size == 0: 

557 logger.warning('Empty models vector.') 

558 return [] 

559 

560 # gms = problem.combine_misfits(history.misfits) 

561 # isort = num.argsort(gms) 

562 # iorder = num.empty_like(isort) 

563 # iorder[isort] = num.arange(iorder.size)[::-1] 

564 

565 ref_source = problem.base_source 

566 

567 mean_source = stats.get_mean_source( 

568 problem, history.models) 

569 

570 best_source = history.get_best_source() 

571 

572 nlines_max = int(round(self.size_cm[1] / 5. * 4. - 1.0)) 

573 

574 if self.cluster_attribute: 

575 cluster_sources = history.mean_sources_by_cluster( 

576 self.cluster_attribute) 

577 else: 

578 cluster_sources = [] 

579 

580 def get_deco(source): 

581 mt = source.pyrocko_moment_tensor() 

582 return mt.standard_decomposition() 

583 

584 lines = [] 

585 lines.append( 

586 ('Ensemble best', get_deco(best_source), mpl_color('aluminium5'))) 

587 

588 lines.append( 

589 ('Ensemble mean', get_deco(mean_source), mpl_color('aluminium5'))) 

590 

591 for (icluster, perc, source) in cluster_sources: 

592 if len(lines) < nlines_max - int(self.show_reference): 

593 lines.append( 

594 (cluster_label(icluster, perc), 

595 get_deco(source), 

596 cluster_color(icluster))) 

597 else: 

598 logger.warning( 

599 'Skipping display of cluster %i because figure height is ' 

600 'too small. Figure height should be at least %g cm.' % ( 

601 icluster, (3 + len(cluster_sources) 

602 + int(self.show_reference)) * 5/4.)) 

603 

604 if self.show_reference: 

605 lines.append( 

606 ('Reference', get_deco(ref_source), mpl_color('aluminium3'))) 

607 

608 moment_full_max = max(deco[-1][0] for (_, deco, _) in lines) 

609 

610 for xpos, label in [ 

611 (0., 'Full'), 

612 (2., 'Isotropic'), 

613 (4., 'Deviatoric'), 

614 (6., 'CLVD'), 

615 (8., 'DC')]: 

616 

617 axes.annotate( 

618 label, 

619 xy=(1 + xpos, nlines_max), 

620 xycoords='data', 

621 xytext=(0., 0.), 

622 textcoords='offset points', 

623 ha='center', 

624 va='center', 

625 color='black', 

626 fontsize=fontsize) 

627 

628 for i, (label, deco, color_t) in enumerate(lines): 

629 ypos = nlines_max - i - 1.0 

630 

631 [(moment_iso, ratio_iso, m_iso), 

632 (moment_dc, ratio_dc, m_dc), 

633 (moment_clvd, ratio_clvd, m_clvd), 

634 (moment_devi, ratio_devi, m_devi), 

635 (moment_full, ratio_full, m_full)] = deco 

636 

637 size0 = moment_full / moment_full_max 

638 

639 axes.annotate( 

640 label, 

641 xy=(-2., ypos), 

642 xycoords='data', 

643 xytext=(0., 0.), 

644 textcoords='offset points', 

645 ha='left', 

646 va='center', 

647 color='black', 

648 fontsize=fontsize) 

649 

650 for xpos, mt_part, ratio, ops in [ 

651 (0., m_full, ratio_full, '-'), 

652 (2., m_iso, ratio_iso, '='), 

653 (4., m_devi, ratio_devi, '='), 

654 (6., m_clvd, ratio_clvd, '+'), 

655 (8., m_dc, ratio_dc, None)]: 

656 

657 if ratio > 1e-4: 

658 try: 

659 beachball.plot_beachball_mpl( 

660 mt_part, axes, 

661 beachball_type='full', 

662 position=(1. + xpos, ypos), 

663 size=0.9 * size0 * math.sqrt(ratio), 

664 size_units='data', 

665 color_t=color_t, 

666 linewidth=1.0) 

667 

668 except beachball.BeachballError as e: 

669 logger.warning(str(e)) 

670 

671 axes.annotate( 

672 'ERROR', 

673 xy=(1. + xpos, ypos), 

674 ha='center', 

675 va='center', 

676 color='red', 

677 fontsize=fontsize) 

678 

679 else: 

680 axes.annotate( 

681 'N/A', 

682 xy=(1. + xpos, ypos), 

683 ha='center', 

684 va='center', 

685 color='black', 

686 fontsize=fontsize) 

687 

688 if ops is not None: 

689 axes.annotate( 

690 ops, 

691 xy=(2. + xpos, ypos), 

692 ha='center', 

693 va='center', 

694 color='black', 

695 fontsize=fontsize) 

696 

697 axes.axison = False 

698 axes.set_xlim(-2.25, 9.75) 

699 axes.set_ylim(-0.5, nlines_max+0.5) 

700 

701 item = PlotItem(name='main') 

702 return [[item, fig]] 

703 

704 

705class MTLocationPlot(SectionPlotConfig): 

706 ''' MT location plot of the best solutions in three cross-sections. ''' 

707 name = 'location_mt' 

708 beachball_type = StringChoice.T( 

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

710 default='dc') 

711 normalisation_gamma = Float.T( 

712 default=3., 

713 help='Normalisation of colors and alpha as :math:`x^\\gamma`.' 

714 'A linear colormap/alpha with :math:`\\gamma=1`.') 

715 

716 def make(self, environ): 

717 environ.setup_modelling() 

718 cm = environ.get_plot_collection_manager() 

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

720 mpl_init(fontsize=self.font_size) 

721 self._to_be_closed = [] 

722 cm.create_group_mpl( 

723 self, 

724 self.draw_figures(history), 

725 title=u'MT Location', 

726 section='solution', 

727 feather_icon='target', 

728 description=u''' 

729Location plot of the ensemble of best solutions in three cross-sections. 

730 

731The coordinate range is defined by the search space given in the config file. 

732Symbols show best double-couple mechanisms, and colors indicate low (red) and 

733high (blue) misfit. 

734''') 

735 for obj in self._to_be_closed: 

736 obj.close() 

737 

738 def draw_figures(self, history, color_p_axis=False): 

739 from matplotlib import colors 

740 

741 color = 'black' 

742 fontsize = self.font_size 

743 markersize = fontsize * 1.5 

744 beachballsize_small = markersize * 0.5 

745 beachball_type = self.beachball_type 

746 

747 problem = history.problem 

748 sp = SectionPlot(config=self) 

749 self._to_be_closed.append(sp) 

750 

751 fig = sp.fig 

752 axes_en = sp.axes_xy 

753 axes_dn = sp.axes_zy 

754 axes_ed = sp.axes_xz 

755 

756 bounds = problem.get_combined_bounds() 

757 

758 models = history.get_sorted_primary_models()[::-1] 

759 

760 iorder = num.arange(history.nmodels) 

761 

762 for parname, set_label, set_lim in [ 

763 ['east_shift', sp.set_xlabel, sp.set_xlim], 

764 ['north_shift', sp.set_ylabel, sp.set_ylim], 

765 ['depth', sp.set_zlabel, sp.set_zlim]]: 

766 

767 ipar = problem.name_to_index(parname) 

768 par = problem.combined[ipar] 

769 set_label(par.get_label()) 

770 xmin, xmax = fixlim(*par.scaled(bounds[ipar])) 

771 set_lim(xmin, xmax) 

772 

773 if 'volume_change' in problem.parameter_names: 

774 volumes = models[:, problem.name_to_index('volume_change')] 

775 volume_max = volumes.max()