Coverage for /usr/local/lib/python3.11/dist-packages/grond/core.py: 51%

546 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-10-26 16:25 +0000

1from __future__ import print_function 

2 

3import sys 

4import logging 

5import time 

6import copy 

7import shutil 

8import glob 

9import math 

10import os 

11import numpy as num 

12from contextlib import contextmanager 

13 

14from pyrocko.guts import Object, String, Float, List 

15from pyrocko import gf, trace, guts, util, weeding 

16from pyrocko import parimap, model, marker as pmarker 

17 

18from .dataset import NotFound, InvalidObject 

19from .problems.base import Problem, load_problem_info_and_data, \ 

20 load_problem_data, ProblemDataNotAvailable 

21 

22from .optimisers.base import BadProblem 

23from .targets.waveform.target import WaveformMisfitResult 

24from .targets.base import dump_misfit_result_collection, \ 

25 MisfitResultCollection, MisfitResult, MisfitResultError 

26from .meta import expand_template, GrondError, selected 

27from .environment import Environment 

28from .monitor import GrondMonitor 

29from .config import get_global_config 

30 

31logger = logging.getLogger('grond.core') 

32guts_prefix = 'grond' 

33op = os.path 

34 

35 

36class RingBuffer(num.ndarray): 

37 def __new__(cls, *args, **kwargs): 

38 cls = num.ndarray.__new__(cls, *args, **kwargs) 

39 cls.fill(0.) 

40 return cls 

41 

42 def __init__(self, *args, **kwargs): 

43 self.pos = 0 

44 

45 def put(self, value): 

46 self[self.pos] = value 

47 self.pos += 1 

48 self.pos %= self.size 

49 

50 

51def mahalanobis_distance(xs, mx, cov): 

52 imask = num.diag(cov) != 0. 

53 icov = num.linalg.inv(cov[imask, :][:, imask]) 

54 temp = xs[:, imask] - mx[imask] 

55 return num.sqrt(num.sum(temp * num.dot(icov, temp.T).T, axis=1)) 

56 

57 

58@contextmanager 

59def lock_rundir(rundir): 

60 statefn = op.join(rundir, '.running') 

61 if op.exists(statefn): 

62 raise EnvironmentError('file %s already exists!' % statefn) 

63 try: 

64 with open(statefn, 'w') as f: 

65 f.write('') 

66 yield True 

67 finally: 

68 os.remove(statefn) 

69 

70 

71class DirectoryAlreadyExists(GrondError): 

72 pass 

73 

74 

75def weed(origin, targets, limit, neighborhood=3): 

76 

77 azimuths = num.zeros(len(targets)) 

78 dists = num.zeros(len(targets)) 

79 for i, target in enumerate(targets): 

80 _, azimuths[i] = target.azibazi_to(origin) 

81 dists[i] = target.distance_to(origin) 

82 

83 badnesses = num.ones(len(targets), dtype=float) 

84 deleted, meandists_kept = weeding.weed( 

85 azimuths, dists, badnesses, 

86 nwanted=limit, 

87 neighborhood=neighborhood) 

88 

89 targets_weeded = [ 

90 target for (delete, target) in zip(deleted, targets) if not delete] 

91 

92 return targets_weeded, meandists_kept, deleted 

93 

94 

95def sarr(a): 

96 return ' '.join('%15g' % x for x in a) 

97 

98 

99def forward(env, show='filtered'): 

100 payload = [] 

101 if env.have_rundir(): 

102 env.setup_modelling() 

103 history = env.get_history(subset='harvest') 

104 xbest = history.get_best_model() 

105 problem = env.get_problem() 

106 ds = env.get_dataset() 

107 payload.append((ds, problem, xbest)) 

108 

109 else: 

110 for event_name in env.get_selected_event_names(): 

111 env.set_current_event_name(event_name) 

112 env.setup_modelling() 

113 problem = env.get_problem() 

114 ds = env.get_dataset() 

115 xref = problem.preconstrain(problem.get_reference_model()) 

116 payload.append((ds, problem, xref)) 

117 

118 all_trs = [] 

119 events = [] 

120 stations = {} 

121 for (ds, problem, x) in payload: 

122 results = problem.evaluate(x) 

123 

124 event = problem.get_source(x).pyrocko_event() 

125 events.append(event) 

126 

127 for result in results: 

128 if isinstance(result, WaveformMisfitResult): 

129 if show == 'filtered': 

130 result.filtered_obs.set_codes(location='ob') 

131 result.filtered_syn.set_codes(location='sy') 

132 all_trs.append(result.filtered_obs) 

133 all_trs.append(result.filtered_syn) 

134 elif show == 'processed': 

135 result.processed_obs.set_codes(location='ob') 

136 result.processed_syn.set_codes(location='sy') 

137 all_trs.append(result.processed_obs) 

138 all_trs.append(result.processed_syn) 

139 else: 

140 raise ValueError('Invalid argument for show: %s' % show) 

141 

142 for station in ds.get_stations(): 

143 stations[station.nsl()] = station 

144 

145 markers = [] 

146 for ev in events: 

147 markers.append(pmarker.EventMarker(ev)) 

148 

149 trace.snuffle(all_trs, markers=markers, stations=list(stations.values())) 

150 

151 

152def harvest( 

153 rundir, problem=None, nbest=10, force=False, weed=0, 

154 export_fits=[]): 

155 

156 env = Environment([rundir]) 

157 optimiser = env.get_optimiser() 

158 nchains = env.get_optimiser().nchains 

159 

160 if problem is None: 

161 problem, xs, misfits, bootstrap_misfits, _ = \ 

162 load_problem_info_and_data(rundir, nchains=nchains) 

163 else: 

164 xs, misfits, bootstrap_misfits, _ = \ 

165 load_problem_data(rundir, problem, nchains=nchains) 

166 

167 logger.info('Harvesting problem "%s"...' % problem.name) 

168 

169 dumpdir = op.join(rundir, 'harvest') 

170 if op.exists(dumpdir): 

171 if force: 

172 shutil.rmtree(dumpdir) 

173 else: 

174 raise DirectoryAlreadyExists( 

175 'Harvest directory already exists: %s' % dumpdir) 

176 

177 util.ensuredir(dumpdir) 

178 

179 ibests_list = [] 

180 ibests = [] 

181 gms = bootstrap_misfits[:, 0] 

182 isort = num.argsort(gms) 

183 

184 ibests_list.append(isort[:nbest]) 

185 

186 if weed != 3: 

187 for ibootstrap in range(optimiser.nbootstrap): 

188 bms = bootstrap_misfits[:, ibootstrap] 

189 isort = num.argsort(bms) 

190 ibests_list.append(isort[:nbest]) 

191 ibests.append(isort[0]) 

192 

193 if weed: 

194 mean_gm_best = num.median(gms[ibests]) 

195 std_gm_best = num.std(gms[ibests]) 

196 ibad = set() 

197 

198 for ibootstrap, ibest in enumerate(ibests): 

199 if gms[ibest] > mean_gm_best + std_gm_best: 

200 ibad.add(ibootstrap) 

201 

202 ibests_list = [ 

203 ibests_ for (ibootstrap, ibests_) in enumerate(ibests_list) 

204 if ibootstrap not in ibad] 

205 

206 ibests = num.concatenate(ibests_list) 

207 

208 if weed == 2: 

209 ibests = ibests[gms[ibests] < mean_gm_best] 

210 

211 for i in ibests: 

212 problem.dump_problem_data(dumpdir, xs[i], misfits[i, :, :]) 

213 

214 if export_fits: 

215 env.setup_modelling() 

216 problem = env.get_problem() 

217 history = env.get_history(subset='harvest') 

218 

219 for what in export_fits: 

220 if what == 'best': 

221 models = [history.get_best_model()] 

222 elif what == 'mean': 

223 models = [history.get_mean_model()] 

224 elif what == 'ensemble': 

225 models = history.models 

226 else: 

227 raise GrondError( 

228 'Invalid option for harvest\'s export_fits argument: %s' 

229 % what) 

230 

231 results = [] 

232 for x in models: 

233 results.append([ 

234 (result if isinstance(result, MisfitResult) 

235 else MisfitResultError(message=str(result))) for 

236 result in problem.evaluate(x)]) 

237 

238 emr = MisfitResultCollection(results=results) 

239 

240 dump_misfit_result_collection( 

241 emr, 

242 op.join(dumpdir, 'fits-%s.yaml' % what)) 

243 

244 logger.info('Done harvesting problem "%s".' % problem.name) 

245 

246 

247def cluster(rundir, clustering, metric): 

248 env = Environment([rundir]) 

249 history = env.get_history(subset='harvest') 

250 problem = history.problem 

251 models = history.models 

252 

253 events = [problem.get_source(model).pyrocko_event() for model in models] 

254 

255 from grond.clustering import metrics 

256 

257 if metric not in metrics.metrics: 

258 raise GrondError('Unknown metric: %s' % metric) 

259 

260 mat = metrics.compute_similarity_matrix(events, metric) 

261 

262 clusters = clustering.perform(mat) 

263 

264 labels = num.sort(num.unique(clusters)) 

265 bins = num.concatenate((labels, [labels[-1]+1])) 

266 ns = num.histogram(clusters, bins)[0] 

267 

268 history.set_attribute('cluster', clusters) 

269 

270 for i in range(labels.size): 

271 if labels[i] == -1: 

272 logging.info( 

273 'Number of unclustered events: %5i' % ns[i]) 

274 else: 

275 logging.info( 

276 'Number of events in cluster %i: %5i' % (labels[i], ns[i])) 

277 

278 

279def get_event_names(config): 

280 return config.get_event_names() 

281 

282 

283def check_problem(problem): 

284 if len(problem.targets) == 0: 

285 raise GrondError('No targets available') 

286 

287 

288def check( 

289 config, 

290 event_names=None, 

291 target_string_ids=None, 

292 show_waveforms=False, 

293 n_random_synthetics=10, 

294 stations_used_path=None): 

295 

296 markers = [] 

297 stations_used = {} 

298 erroneous = [] 

299 for ievent, event_name in enumerate(event_names): 

300 ds = config.get_dataset(event_name) 

301 event = ds.get_event() 

302 trs_all = [] 

303 try: 

304 problem = config.get_problem(event) 

305 

306 _, nfamilies = problem.get_family_mask() 

307 logger.info('Problem: %s' % problem.name) 

308 logger.info('Number of target families: %i' % nfamilies) 

309 logger.info('Number of targets (total): %i' % len(problem.targets)) 

310 

311 if target_string_ids: 

312 problem.targets = [ 

313 target for target in problem.targets 

314 if util.match_nslc(target_string_ids, target.string_id())] 

315 

316 logger.info( 

317 'Number of targets (selected): %i' % len(problem.targets)) 

318 

319 check_problem(problem) 

320 

321 results_list = [] 

322 sources = [] 

323 if n_random_synthetics == 0: 

324 x = problem.preconstrain(problem.get_reference_model()) 

325 sources.append(problem.base_source) 

326 results = problem.evaluate(x) 

327 results_list.append(results) 

328 

329 else: 

330 for i in range(n_random_synthetics): 

331 x = problem.get_random_model() 

332 sources.append(problem.get_source(x)) 

333 results = problem.evaluate(x) 

334 results_list.append(results) 

335 

336 if show_waveforms: 

337 engine = config.engine_config.get_engine() 

338 times = [] 

339 tdata = [] 

340 for target in problem.targets: 

341 tobs_shift_group = [] 

342 tcuts = [] 

343 for source in sources: 

344 tmin_fit, tmax_fit, tfade, tfade_taper = \ 

345 target.get_taper_params(engine, source) 

346 

347 times.extend((tmin_fit-tfade*2., tmax_fit+tfade*2.)) 

348 

349 tobs, tsyn = target.get_pick_shift(engine, source) 

350 if None not in (tobs, tsyn): 

351 tobs_shift = tobs - tsyn 

352 else: 

353 tobs_shift = 0.0 

354 

355 tcuts.append(target.get_cutout_timespan( 

356 tmin_fit+tobs_shift, tmax_fit+tobs_shift, tfade)) 

357 

358 tobs_shift_group.append(tobs_shift) 

359 

360 tcuts = num.array(tcuts, dtype=float) 

361 

362 tdata.append(( 

363 tfade, 

364 num.mean(tobs_shift_group), 

365 (num.min(tcuts[:, 0]), num.max(tcuts[:, 1])))) 

366 

367 tmin = min(times) 

368 tmax = max(times) 

369 

370 tmax += (tmax-tmin)*2 

371 

372 for (tfade, tobs_shift, tcut), target in zip( 

373 tdata, problem.targets): 

374 

375 store = engine.get_store(target.store_id) 

376 

377 deltat = store.config.deltat 

378 

379 freqlimits = list(target.get_freqlimits()) 

380 freqlimits[2] = 0.45/deltat 

381 freqlimits[3] = 0.5/deltat 

382 freqlimits = tuple(freqlimits) 

383 

384 try: 

385 trs_projected, trs_restituted, trs_raw, _ = \ 

386 ds.get_waveform( 

387 target.codes, 

388 tmin=tmin+tobs_shift, 

389 tmax=tmax+tobs_shift, 

390 tfade=tfade, 

391 freqlimits=freqlimits, 

392 deltat=deltat, 

393 backazimuth=target. 

394 get_backazimuth_for_waveform(), 

395 debug=True) 

396 

397 except NotFound as e: 

398 logger.warning(str(e)) 

399 continue 

400 

401 trs_projected = copy.deepcopy(trs_projected) 

402 trs_restituted = copy.deepcopy(trs_restituted) 

403 trs_raw = copy.deepcopy(trs_raw) 

404 

405 for trx in trs_projected + trs_restituted + trs_raw: 

406 trx.shift(-tobs_shift) 

407 trx.set_codes( 

408 network='', 

409 station=target.string_id(), 

410 location='') 

411 

412 for trx in trs_projected: 

413 trx.set_codes(location=trx.location + '2_proj') 

414 

415 for trx in trs_restituted: 

416 trx.set_codes(location=trx.location + '1_rest') 

417 

418 for trx in trs_raw: 

419 trx.set_codes(location=trx.location + '0_raw') 

420 

421 trs_all.extend(trs_projected) 

422 trs_all.extend(trs_restituted) 

423 trs_all.extend(trs_raw) 

424 

425 for source in sources: 

426 tmin_fit, tmax_fit, tfade, tfade_taper = \ 

427 target.get_taper_params(engine, source) 

428 

429 markers.append(pmarker.Marker( 

430 nslc_ids=[('', target.string_id(), '*_proj', '*')], 

431 tmin=tmin_fit, tmax=tmax_fit)) 

432 

433 markers.append(pmarker.Marker( 

434 nslc_ids=[('', target.string_id(), '*_raw', '*')], 

435 tmin=tcut[0]-tobs_shift, tmax=tcut[1]-tobs_shift, 

436 kind=1)) 

437 

438 else: 

439 for itarget, target in enumerate(problem.targets): 

440 

441 nok = 0 

442 for results in results_list: 

443 result = results[itarget] 

444 if not isinstance(result, gf.SeismosizerError): 

445 nok += 1 

446 

447 if nok == 0: 

448 sok = 'not used' 

449 elif nok == len(results_list): 

450 sok = 'ok' 

451 try: 

452 s = ds.get_station(target) 

453 stations_used[s.nsl()] = s 

454 except (NotFound, InvalidObject): 

455 pass 

456 else: 

457 sok = 'not used (%i/%i ok)' % (nok, len(results_list)) 

458 

459 logger.info('%-40s %s' % ( 

460 (target.string_id() + ':', sok))) 

461 

462 except GrondError as e: 

463 logger.error('Event %i, "%s": %s' % ( 

464 ievent, 

465 event.name or util.time_to_str(event.time), 

466 str(e))) 

467 

468 erroneous.append(event) 

469 

470 if show_waveforms: 

471 trace.snuffle(trs_all, stations=ds.get_stations(), markers=markers) 

472 

473 if stations_used_path: 

474 stations = list(stations_used.values()) 

475 stations.sort(key=lambda s: s.nsl()) 

476 model.dump_stations(stations, stations_used_path) 

477 

478 if erroneous: 

479 raise GrondError( 

480 'Check failed for events: %s' 

481 % ', '.join(ev.name for ev in erroneous)) 

482 

483 

484g_state = {} 

485 

486 

487def go(environment, force=False, preserve=False): 

488 

489 global_config = get_global_config() 

490 

491 g_data = (environment, force, preserve, global_config, False) 

492 g_state[id(g_data)] = g_data 

493 

494 nevents = environment.nevents_selected 

495 for x in parimap.parimap( 

496 process_event, 

497 range(environment.nevents_selected), 

498 [id(g_data)] * nevents, 

499 nprocs=global_config.nparallel): 

500 

501 pass 

502 

503 

504def continue_run(environment, force=False, preserve=False): 

505 

506 global_config = get_global_config() 

507 

508 g_data = (environment, force, preserve, global_config, True) 

509 g_state[id(g_data)] = g_data 

510 

511 nevents = environment.nevents_selected 

512 for x in parimap.parimap( 

513 process_event, 

514 range(environment.nevents_selected), 

515 [id(g_data)] * nevents, 

516 nprocs=global_config.nparallel): 

517 

518 pass 

519 

520 

521def process_event(ievent, g_data_id): 

522 

523 env, force, preserve, global_config, continue_run = g_state[g_data_id] 

524 

525 config = env.get_config() 

526 event_name = env.get_selected_event_names()[ievent] 

527 nevents = env.nevents_selected 

528 tstart = time.time() 

529 

530 ds = config.get_dataset(event_name) 

531 event = ds.get_event() 

532 problem = config.get_problem(event) 

533 

534 synt = ds.synthetic_test 

535 if synt: 

536 problem.base_source = problem.get_source(synt.get_x()) 

537 

538 check_problem(problem) 

539 

540 rundir = expand_template( 

541 config.rundir_template, 

542 dict(problem_name=problem.name)) 

543 env.set_rundir_path(rundir) 

544 nold_rundirs = len(glob.glob(rundir + '*')) 

545 

546 if op.exists(rundir) and not continue_run: 

547 if preserve: 

548 shutil.move(rundir, rundir+'-old-%d' % nold_rundirs) 

549 elif force: 

550 shutil.rmtree(rundir) 

551 else: 

552 logger.warning('Skipping problem "%s": rundir already exists: %s', 

553 problem.name, rundir) 

554 return 

555 

556 if op.exists(rundir) and continue_run: 

557 logger.info( 

558 'Continuing event %i / %i', ievent + 1, nevents) 

559 

560 env_old = Environment(rundir) 

561 

562 history = env_old.get_history() 

563 targets = env_old.get_problem().targets 

564 for target in targets: 

565 target.set_dataset(ds) 

566 

567 problem.targets = targets 

568 

569 if preserve: 

570 shutil.copytree(rundir, rundir+'-old-%d' % nold_rundirs) 

571 

572 elif not op.exists(rundir) and continue_run: 

573 logger.warning('Cannot find rundir %s to continue...', rundir) 

574 return 

575 

576 else: 

577 logger.info( 

578 'Starting event %i / %i', ievent+1, nevents) 

579 history = None 

580 

581 util.ensuredir(rundir) 

582 logger.info('Rundir: %s', rundir) 

583 

584 optimiser = config.optimiser_config.get_optimiser() 

585 

586 if not continue_run: 

587 logger.info('Analysing problem "%s".', problem.name) 

588 for analyser_conf in config.analyser_configs: 

589 analyser = analyser_conf.get_analyser() 

590 analyser.analyse(problem, ds) 

591 

592 basepath = config.get_basepath() 

593 config.change_basepath(rundir) 

594 guts.dump(config, filename=op.join(rundir, 'config.yaml')) 

595 config.change_basepath(basepath) 

596 optimiser.init_bootstraps(problem) 

597 

598 problem.dump_problem_info(rundir) 

599 

600 monitor = None 

601 

602 xs_inject = None 

603 synt = ds.synthetic_test 

604 if synt and synt.inject_solution: 

605 xs_inject = synt.get_x()[num.newaxis, :] 

606 

607 try: 

608 if xs_inject is not None: 

609 from .optimisers import highscore 

610 if not isinstance(optimiser, highscore.HighScoreOptimiser): 

611 raise GrondError( 

612 'Optimiser does not support injections.') 

613 

614 optimiser.sampler_phases[0:0] = [ 

615 highscore.InjectionSamplerPhase(xs_inject=xs_inject)] 

616 

617 gconf = get_global_config() 

618 with lock_rundir(rundir): 

619 if gconf.status == 'state': 

620 monitor = GrondMonitor.watch(rundir) 

621 optimiser.optimise( 

622 problem, 

623 rundir=rundir, 

624 history=history) 

625 

626 harvest(rundir, problem, force=True) 

627 

628 except BadProblem as e: 

629 logger.error(str(e)) 

630 

631 except GrondError as e: 

632 logger.error(str(e)) 

633 

634 finally: 

635 if monitor: 

636 monitor.terminate() 

637 

638 tstop = time.time() 

639 logger.info( 

640 'Stop %i / %i (%g min)', ievent+1, nevents, (tstop - tstart)/60.) 

641 

642 logger.info( 

643 'Done with problem "%s", rundir is "%s".', problem.name, rundir) 

644 

645 

646class ParameterStats(Object): 

647 name = String.T() 

648 mean = Float.T() 

649 std = Float.T() 

650 best = Float.T() 

651 minimum = Float.T() 

652 percentile5 = Float.T() 

653 percentile16 = Float.T() 

654 median = Float.T() 

655 percentile84 = Float.T() 

656 percentile95 = Float.T() 

657 maximum = Float.T() 

658 

659 def __init__(self, *args, **kwargs): 

660 kwargs.update(zip(self.T.propnames, args)) 

661 Object.__init__(self, **kwargs) 

662 

663 def get_values_dict(self): 

664 return dict( 

665 (self.name+'.' + k, getattr(self, k)) 

666 for k in self.T.propnames 

667 if k != 'name') 

668 

669 

670class ResultStats(Object): 

671 problem = Problem.T() 

672 parameter_stats_list = List.T(ParameterStats.T()) 

673 

674 def get_values_dict(self): 

675 d = {} 

676 for ps in self.parameter_stats_list: 

677 d.update(ps.get_values_dict()) 

678 return d 

679 

680 

681def make_stats(problem, models, gms, pnames=None): 

682 ibest = num.argmin(gms) 

683 rs = ResultStats(problem=problem) 

684 if pnames is None: 

685 pnames = problem.parameter_names 

686 

687 for pname in pnames: 

688 iparam = problem.name_to_index(pname) 

689 vs = problem.extract(models, iparam) 

690 mi, p5, p16, median, p84, p95, ma = map(float, num.percentile( 

691 vs, [0., 5., 16., 50., 84., 95., 100.])) 

692 

693 mean = float(num.mean(vs)) 

694 std = float(num.std(vs)) 

695 best = float(vs[ibest]) 

696 s = ParameterStats( 

697 pname, mean, std, best, mi, p5, p16, median, p84, p95, ma) 

698 

699 rs.parameter_stats_list.append(s) 

700 

701 return rs 

702 

703 

704def try_add_location_uncertainty(data, types): 

705 vs = [data.get(k, None) for k in ( 

706 'north_shift.std', 'east_shift.std', 'depth.std')] 

707 

708 if None not in vs: 

709 data['location_uncertainty'] = math.sqrt(sum(v**2 for v in vs)) 

710 types['location_uncertainty'] = float 

711 

712 

713def format_stats(rs, fmt): 

714 pname_to_pindex = dict( 

715 (p.name, i) for (i, p) in enumerate(rs.parameter_stats_list)) 

716 

717 values = [] 

718 headers = [] 

719 for x in fmt: 

720 if x == 'problem.name': 

721 headers.append(x) 

722 values.append('%-16s' % rs.problem.name) 

723 else: 

724 pname, qname = x.split('.') 

725 pindex = pname_to_pindex[pname] 

726 values.append( 

727 '%16.7g' % getattr(rs.parameter_stats_list[pindex], qname)) 

728 headers.append(x) 

729 

730 return ' '.join(values) 

731 

732 

733def export( 

734 what, rundirs, type=None, pnames=None, filename=None, selection=None, 

735 effective_lat_lon=False): 

736 

737 if pnames is not None: 

738 pnames_clean = [ 

739 pname.split('.')[0] for pname in pnames 

740 if not pname.startswith('problem.')] 

741 shortform = all(len(pname.split('.')) == 2 for pname in pnames) 

742 else: 

743 pnames_clean = None 

744 shortform = False 

745 

746 if what == 'stats' and type is not None: 

747 raise GrondError('Invalid argument combination: what=%s, type=%s' % ( 

748 repr(what), repr(type))) 

749 

750 if what != 'stats' and shortform: 

751 raise GrondError('Invalid argument combination: what=%s, pnames=%s' % ( 

752 repr(what), repr(pnames))) 

753 

754 if what != 'stats' and type != 'vector' and pnames is not None: 

755 raise GrondError( 

756 'Invalid argument combination: what=%s, type=%s, pnames=%s' % ( 

757 repr(what), repr(type), repr(pnames))) 

758 

759 if filename is None: 

760 out = sys.stdout 

761 else: 

762 out = open(filename, 'w') 

763 

764 if type is None: 

765 type = 'event' 

766 

767 if shortform: 

768 print('#', ' '.join(['%16s' % x for x in pnames]), file=out) 

769 

770 def dump(x, gm, indices): 

771 if type == 'vector': 

772 print(' ', ' '.join( 

773 '%16.7g' % problem.extract(x, i) for i in indices), 

774 '%16.7g' % gm, file=out) 

775 

776 elif type == 'source': 

777 source = problem.get_source(x) 

778 if effective_lat_lon: 

779 source.set_origin(*source.effective_latlon) 

780 guts.dump(source, stream=out) 

781 

782 elif type == 'event': 

783 ev = problem.get_source(x).pyrocko_event() 

784 if effective_lat_lon: 

785 ev.set_origin(*ev.effective_latlon) 

786 

787 model.dump_events([ev], stream=out) 

788 

789 elif type == 'event-yaml': 

790 ev = problem.get_source(x).pyrocko_event() 

791 if effective_lat_lon: 

792 ev.set_origin(*ev.effective_latlon) 

793 guts.dump_all([ev], stream=out) 

794 

795 else: 

796 raise GrondError('Invalid argument: type=%s' % repr(type)) 

797 

798 header = None 

799 for rundir in rundirs: 

800 env = Environment(rundir) 

801 info = env.get_run_info() 

802 

803 try: 

804 history = env.get_history(subset='harvest') 

805 except ProblemDataNotAvailable as e: 

806 logger.error( 

807 'Harvest not available (Did the run succeed?): %s' % str(e)) 

808 continue 

809 

810 problem = history.problem 

811 models = history.models 

812 misfits = history.get_primary_chain_misfits() 

813 

814 config = env.get_config() 

815 engine = config.engine_config.get_engine() 

816 problem.set_engine(engine) 

817 

818 if selection: 

819 rs = make_stats( 

820 problem, models, 

821 history.get_primary_chain_misfits()) 

822 

823 data = dict(tags=info.tags) 

824 types = dict(tags=(list, str)) 

825 

826 for k, v in rs.get_values_dict().items(): 

827 data[k] = v 

828 types[k] = float 

829 

830 try_add_location_uncertainty(data, types) 

831 

832 if not selected(selection, data=data, types=types): 

833 continue 

834 

835 else: 

836 rs = None 

837 

838 if type == 'vector': 

839 pnames_take = pnames_clean or \ 

840 problem.parameter_names[:problem.nparameters] 

841 

842 indices = num.array( 

843 [problem.name_to_index(pname) for pname in pnames_take]) 

844 

845 if type == 'vector' and what in ('best', 'mean', 'ensemble'): 

846 extra = ['global_misfit'] 

847 else: 

848 extra = [] 

849 

850 new_header = '# ' + ' '.join( 

851 '%16s' % x for x in pnames_take + extra) 

852 

853 if type == 'vector' and header != new_header: 

854 print(new_header, file=out) 

855 

856 header = new_header 

857 else: 

858 indices = None 

859 

860 if what == 'best': 

861 x_best = history.get_best_model() 

862 gm_best = history.get_best_misfit() 

863 dump(x_best, gm_best, indices) 

864 

865 elif what == 'mean': 

866 x_mean = history.get_mean_model() 

867 gm_mean = history.get_mean_misfit(chain=0) 

868 dump(x_mean, gm_mean, indices) 

869 

870 elif what == 'ensemble': 

871 isort = num.argsort(misfits) 

872 for i in isort: 

873 dump(models[i], misfits[i], indices) 

874 

875 elif what == 'stats': 

876 if not rs: 

877 rs = make_stats(problem, models, 

878 history.get_primary_chain_misfits(), 

879 pnames_clean) 

880 

881 if shortform: 

882 print(' ', format_stats(rs, pnames), file=out) 

883 else: 

884 print(rs, file=out) 

885 

886 else: 

887 raise GrondError('Invalid argument: what=%s' % repr(what)) 

888 

889 if out is not sys.stdout: 

890 out.close() 

891 

892 

893__all__ = ''' 

894 DirectoryAlreadyExists 

895 forward 

896 harvest 

897 cluster 

898 go 

899 continue_run 

900 get_event_names 

901 check 

902 export 

903'''.split()