load_problem_data, ProblemDataNotAvailable
MisfitResultCollection, MisfitResult, MisfitResultError
cls = num.ndarray.__new__(cls, *args, **kwargs) cls.fill(0.) return cls
self.pos = 0
self[self.pos] = value self.pos += 1 self.pos %= self.size
imask = num.diag(cov) != 0. icov = num.linalg.inv(cov[imask, :][:, imask]) temp = xs[:, imask] - mx[imask] return num.sqrt(num.sum(temp * num.dot(icov, temp.T).T, axis=1))
def lock_rundir(rundir): raise EnvironmentError('file %s already exists!' % statefn) finally:
azimuths = num.zeros(len(targets)) dists = num.zeros(len(targets)) for i, target in enumerate(targets): _, azimuths[i] = target.azibazi_to(origin) dists[i] = target.distance_to(origin)
badnesses = num.ones(len(targets), dtype=float) deleted, meandists_kept = weeding.weed( azimuths, dists, badnesses, nwanted=limit, neighborhood=neighborhood)
targets_weeded = [ target for (delete, target) in zip(deleted, targets) if not delete]
return targets_weeded, meandists_kept, deleted
return ' '.join('%15g' % x for x in a)
payload = [] if env.have_rundir(): env.setup_modelling() history = env.get_history(subset='harvest') xbest = history.get_best_model() problem = env.get_problem() ds = env.get_dataset() payload.append((ds, problem, xbest))
else: for event_name in env.get_selected_event_names(): env.set_current_event_name(event_name) env.setup_modelling() problem = env.get_problem() ds = env.get_dataset() xref = problem.preconstrain(problem.get_reference_model()) payload.append((ds, problem, xref))
all_trs = [] events = [] stations = {} for (ds, problem, x) in payload: results = problem.evaluate(x)
event = problem.get_source(x).pyrocko_event() events.append(event)
for result in results: if isinstance(result, WaveformMisfitResult): if show == 'filtered': result.filtered_obs.set_codes(location='ob') result.filtered_syn.set_codes(location='sy') all_trs.append(result.filtered_obs) all_trs.append(result.filtered_syn) elif show == 'processed': result.processed_obs.set_codes(location='ob') result.processed_syn.set_codes(location='sy') all_trs.append(result.processed_obs) all_trs.append(result.processed_syn) else: raise ValueError('Invalid argument for show: %s' % show)
for station in ds.get_stations(): stations[station.nsl()] = station
markers = [] for ev in events: markers.append(pmarker.EventMarker(ev))
trace.snuffle(all_trs, markers=markers, stations=list(stations.values()))
rundir, problem=None, nbest=10, force=False, weed=0, export_fits=[]):
load_problem_info_and_data(rundir, nchains=nchains) else: load_problem_data(rundir, problem, nchains=nchains)
else: raise DirectoryAlreadyExists( 'Harvest directory already exists: %s' % dumpdir)
mean_gm_best = num.median(gms[ibests]) std_gm_best = num.std(gms[ibests]) ibad = set()
for ibootstrap, ibest in enumerate(ibests): if gms[ibest] > mean_gm_best + std_gm_best: ibad.add(ibootstrap)
ibests_list = [ ibests_ for (ibootstrap, ibests_) in enumerate(ibests_list) if ibootstrap not in ibad]
ibests = ibests[gms[ibests] < mean_gm_best]
elif what == 'ensemble': models = history.models else: raise GrondError( 'Invalid option for harvest\'s export_fits argument: %s' % what)
(result if isinstance(result, MisfitResult) else MisfitResultError(message=str(result))) for result in problem.evaluate(x)])
emr, op.join(dumpdir, 'fits-%s.yaml' % what))
raise GrondError('Unknown metric: %s' % metric)
'Number of unclustered events: %5i' % ns[i]) else: 'Number of events in cluster %i: %5i' % (labels[i], ns[i]))
raise GrondError('No targets available')
config, event_names=None, target_string_ids=None, show_waveforms=False, n_random_synthetics=10, stations_used_path=None):
problem.targets = [ target for target in problem.targets if util.match_nslc(target_string_ids, target.string_id())]
'Number of targets (selected): %i' % len(problem.targets))
x = problem.preconstrain(problem.get_reference_model()) sources.append(problem.base_source) results = problem.evaluate(x) results_list.append(results)
else:
engine = config.engine_config.get_engine() times = [] tdata = [] for target in problem.targets: tobs_shift_group = [] tcuts = [] for source in sources: tmin_fit, tmax_fit, tfade, tfade_taper = \ target.get_taper_params(engine, source)
times.extend((tmin_fit-tfade*2., tmax_fit+tfade*2.))
tobs, tsyn = target.get_pick_shift(engine, source) if None not in (tobs, tsyn): tobs_shift = tobs - tsyn else: tobs_shift = 0.0
tcuts.append(target.get_cutout_timespan( tmin_fit+tobs_shift, tmax_fit+tobs_shift, tfade))
tobs_shift_group.append(tobs_shift)
tcuts = num.array(tcuts, dtype=num.float)
tdata.append(( tfade, num.mean(tobs_shift_group), (num.min(tcuts[:, 0]), num.max(tcuts[:, 1]))))
tmin = min(times) tmax = max(times)
tmax += (tmax-tmin)*2
for (tfade, tobs_shift, tcut), target in zip( tdata, problem.targets):
store = engine.get_store(target.store_id)
deltat = store.config.deltat
freqlimits = list(target.get_freqlimits()) freqlimits[2] = 0.45/deltat freqlimits[3] = 0.5/deltat freqlimits = tuple(freqlimits)
try: trs_projected, trs_restituted, trs_raw, _ = \ ds.get_waveform( target.codes, tmin=tmin+tobs_shift, tmax=tmax+tobs_shift, tfade=tfade, freqlimits=freqlimits, deltat=deltat, backazimuth=target. get_backazimuth_for_waveform(), debug=True)
except NotFound as e: logger.warn(str(e)) continue
trs_projected = copy.deepcopy(trs_projected) trs_restituted = copy.deepcopy(trs_restituted) trs_raw = copy.deepcopy(trs_raw)
for trx in trs_projected + trs_restituted + trs_raw: trx.shift(-tobs_shift) trx.set_codes( network='', station=target.string_id(), location='')
for trx in trs_projected: trx.set_codes(location=trx.location + '2_proj')
for trx in trs_restituted: trx.set_codes(location=trx.location + '1_rest')
for trx in trs_raw: trx.set_codes(location=trx.location + '0_raw')
trs_all.extend(trs_projected) trs_all.extend(trs_restituted) trs_all.extend(trs_raw)
for source in sources: tmin_fit, tmax_fit, tfade, tfade_taper = \ target.get_taper_params(engine, source)
markers.append(pmarker.Marker( nslc_ids=[('', target.string_id(), '*_proj', '*')], tmin=tmin_fit, tmax=tmax_fit))
markers.append(pmarker.Marker( nslc_ids=[('', target.string_id(), '*_raw', '*')], tmin=tcut[0]-tobs_shift, tmax=tcut[1]-tobs_shift, kind=1))
else:
else: sok = 'not used (%i/%i ok)' % (nok, len(results_list))
(target.string_id() + ':', sok)))
except GrondError as e: logger.error('Event %i, "%s": %s' % ( ievent, event.name or util.time_to_str(event.time), str(e)))
erroneous.append(event)
trace.snuffle(trs_all, stations=ds.get_stations(), markers=markers)
raise GrondError( 'Check failed for events: %s' % ', '.join(ev.name for ev in erroneous))
force=False, preserve=False, nparallel=1, status='state', nthreads=0):
status, nparallel, nthreads)
process_event, range(environment.nevents_selected), [id(g_data)] * nevents, nprocs=nparallel):
g_state[g_data_id]
problem.base_source = problem.get_source(synt.get_x())
config.rundir_template, dict(problem_name=problem.name))
if preserve: nold_rundirs = len(glob.glob(rundir + '*')) shutil.move(rundir, rundir+'-old-%d' % (nold_rundirs)) elif force: shutil.rmtree(rundir) else: logger.warn('Skipping problem "%s": rundir already exists: %s' % (problem.name, rundir)) return
'Starting event %i / %i' % (ievent+1, nevents))
xs_inject = synt.get_x()[num.newaxis, :]
from .optimisers import highscore if not isinstance(optimiser, highscore.HighScoreOptimiser): raise GrondError( 'Optimiser does not support injections.')
optimiser.sampler_phases[0:0] = [ highscore.InjectionSamplerPhase(xs_inject=xs_inject)]
problem, rundir=rundir)
except BadProblem as e: logger.error(str(e))
except GrondError as e: logger.error(str(e))
finally:
'Stop %i / %i (%g min)' % (ievent+1, nevents, (tstop - tstart)/60.))
'Done with problem "%s", rundir is "%s".' % (problem.name, rundir))
return dict( (self.name+'.' + k, getattr(self, k)) for k in self.T.propnames if k != 'name')
d = {} for ps in self.parameter_stats_list: d.update(ps.get_values_dict()) return d
vs, [0., 5., 16., 50., 84., 95., 100.]))
pname, mean, std, best, mi, p5, p16, median, p84, p95, ma)
vs = [data.get(k, None) for k in ( 'north_shift.std', 'east_shift.std', 'depth.std')]
if None not in vs: data['location_uncertainty'] = math.sqrt(sum(v**2 for v in vs)) types['location_uncertainty'] = float
pname_to_pindex = dict( (p.name, i) for (i, p) in enumerate(rs.parameter_stats_list))
values = [] headers = [] for x in fmt: pname, qname = x.split('.') pindex = pname_to_pindex[pname] values.append(getattr(rs.parameter_stats_list[pindex], qname)) headers.append(x)
return ' '.join('%16.7g' % v for v in values)
what, rundirs, type=None, pnames=None, filename=None, selection=None):
pnames_clean = [pname.split('.')[0] for pname in pnames] shortform = all(len(pname.split('.')) == 2 for pname in pnames) else:
raise GrondError('Invalid argument combination: what=%s, type=%s' % ( repr(what), repr(type)))
raise GrondError('Invalid argument combination: what=%s, pnames=%s' % ( repr(what), repr(pnames)))
raise GrondError( 'Invalid argument combination: what=%s, type=%s, pnames=%s' % ( repr(what), repr(type), repr(pnames)))
else:
print('#', ' '.join(['%16s' % x for x in pnames]), file=out)
print(' ', ' '.join( '%16.7g' % problem.extract(x, i) for i in indices), '%16.7g' % gm, file=out)
source = problem.get_source(x) guts.dump(source, stream=out)
else: raise GrondError('Invalid argument: type=%s' % repr(type))
except ProblemDataNotAvailable as e: logger.error( 'Harvest not available (Did the run succeed?): %s' % str(e)) continue
rs = make_stats( problem, models, history.get_primary_chain_misfits())
data = dict(tags=info.tags) types = dict(tags=(list, str))
for k, v in rs.get_values_dict().items(): data[k] = v types[k] = float
try_add_location_uncertainty(data, types)
if not selected(selection, data=data, types=types): continue
else:
pnames_take = pnames_clean or \ problem.parameter_names[:problem.nparameters]
indices = num.array( [problem.name_to_index(pname) for pname in pnames_take])
if type == 'vector' and what in ('best', 'mean', 'ensemble'): extra = ['global_misfit'] else: extra = []
new_header = '# ' + ' '.join( '%16s' % x for x in pnames_take + extra)
if type == 'vector' and header != new_header: print(new_header, file=out)
header = new_header else:
history.get_primary_chain_misfits(), pnames_clean)
print(' ', format_stats(rs, pnames), file=out) else:
else: raise GrondError('Invalid argument: what=%s' % repr(what))
DirectoryAlreadyExists forward harvest cluster go get_event_names check export '''.split() |