Commit b55c490a authored by Marius Isken's avatar Marius Isken
Browse files

Cleaning up // rebase

parent e0dc12e1
from core import * # noqa
from cmt import * # noqa
from double_dc import * # noqa
from dataset import * # noqa
from .problems import * # noqa
......@@ -617,6 +617,43 @@ tautoshift**2 / tautoshift_max**2``
return result
class MisfitSatelliteTarget(gf.SatelliteTarget):
scene_id = guts.String.T()
super_group = gf.StringID.T()
group = gf.StringID.T()
def __init__(self, *args, **kwargs):
gf.SatelliteTarget.__init__(self, *args, **kwargs)
self._ds = None
def set_dataset(self, ds):
self._ds = ds
def get_dataset(self):
return self._ds
def post_process(self, engine, source, statics):
scene = self._ds.get_kite_scene(self.scene_id)
stat_obs = scene.quadtree.leaf_medians
stat_syn = statics['displacement.los']
mf = misfit_static(
stat_obs,
stat_syn,
exponent=2)
wf = mf * scene.covariance.weight_matrix_focal
return wf
def misfit_static(stat_obs, stat_syn, exponent):
mv = num.abs(stat_obs - stat_syn)**exponent
result = MisfitResult(
misfit_value=mv,
misfit_norm=exponent)
return result
def _process(tr, tmin, tmax, taper, domain):
tr_proc = _extend_extract(tr, tmin, tmax)
tr_proc.taper(taper)
......@@ -1004,6 +1041,33 @@ class TargetConfig(Object):
origin = event
targets = []
for scene in ds.get_kite_scenes():
qt = scene.quadtree
lats = num.empty(qt.nleafs)
lons = num.empty(qt.nleafs)
lats.fill(qt.frame.llLat)
lons.fill(qt.frame.llLon)
east_shifts = qt.leaf_focal_points()[:, 0]
north_shifts = qt.leaf_focal_points()[:, 1]
tsnapshot = scene.meta.time_slave - scene.meta.time_master
target = MisfitSatelliteTarget(
quantity='displacement',
scene_id=scene.meta.scene_id,
lats=lats,
lons=lons,
east_shifts=east_shifts,
north_shifts=north_shifts,
tsnapshot=tsnapshot,
interpolation=self.interpolation,
store_id=self.store_id,
super_group=self.super_group,
group=self.group or default_group)
for st in ds.get_stations():
for cha in self.channels:
if ds.is_blacklisted((st.nsl() + (cha,))):
......@@ -1066,6 +1130,16 @@ class TargetConfig(Object):
return targets
class SatelliteTargetConfig(Object):
super_group = gf.StringID.T(default='', optional=True)
group = gf.StringID.T(optional=True)
limit = Int.T(optional=True)
inner_misfit_config = InnerMisfitConfig.T()
interpolation = gf.InterpolationMethod.T()
store_id = gf.StringID.T()
weight = Float.T(default=1.0)
class AnalyserConfig(Object):
niter = Int.T(default=1000)
......
......@@ -76,7 +76,7 @@ class Dataset(object):
self.apply_correction_delays = True
self.apply_correction_factors = True
self.clip_handling = 'by_nsl'
self.static_scenes = []
self.kite_scenes = []
self.synthetic_test = None
self._picks = None
self._cache = {}
......@@ -226,7 +226,7 @@ class Dataset(object):
raise ImportError('Module kite could not be imported, '
'please install from http://pyrocko.org')
for file in filenames:
self.static_scenes.append(read(file))
self.kite_scenes.append(read(file))
def is_blacklisted(self, obj):
try:
......@@ -319,6 +319,17 @@ class Dataset(object):
if not self.is_blacklisted(self.stations[k])
and self.is_whitelisted(self.stations[k])]
def get_kite_displacement(self, scene_id=None):
if scene_id is None:
if len(self.kite_scene) == 0:
raise AttributeError('No kite displacements defined')
return self.kite_scenes[0]
else:
for scene in self.kite_scenes:
if scene.meta.scene_id is scene_id:
return scene
raise AttributeError('No kite scene with id %s defined' % scene_id)
def get_response(self, obj):
if (self.responses is None or len(self.responses) == 0) \
and (self.responses_stationxml is None
......@@ -649,9 +660,6 @@ class Dataset(object):
cache[nslc, tmin, tmax] = e
raise
def get_displacement(self):
return
def get_events(self, magmin=None, event_names=None):
evs = []
for ev in self.events:
......
from cmt import * # noqa
from double_dc import * # noqa
from rectangular import * # noqa
import math
import logging
import numpy as num
from pyrocko import gf, moment_tensor as mtm
from pyrocko.guts import Float, String, Dict, List, Int
from grond import core
guts_prefix = 'grond'
logger = logging.getLogger('grond.cmt')
km = 1000.
as_km = dict(scale_factor=km, scale_unit='km')
class RectangularProblem(core.Problem):
parameters = [
core.Parameter('north_shift', 'm', label='Northing', **as_km),
core.Parameter('east_shift', 'm', label='Easting', **as_km),
core.Parameter('depth', 'm', label='Depth', **as_km),
core.Parameter('length', 'm', label='Length', **as_km),
core.Parameter('width', 'm', label='Width', **as_km),
core.Parameter('dip', 'deg', label='Dip'),
core.Parameter('strike', 'deg', label='Strike'),
core.Parameter('rake', 'deg', label='Rake'),
core.Parameter('slip', 'm', label='Slip'),
]
dependants = []
targets = List.T(gf.Target.T())
ranges = Dict.T(String.T(), gf.Range.T())
distance_min = Float.T(default=0.0)
nbootstrap = Int.T(default=10)
def pack(self, source):
return self.parameter_array(source)
def unpack(self, x):
source = self.base_source.clone(**self.parameter_dict(x))
return source
def extract(self, xs, i):
if xs.ndim == 1:
return self.extract(xs[num.newaxis, :], i)[0]
return xs[:, i]
def random_uniform(self, xbounds):
x = num.zeros(self.nparameters)
for i in [0, 1, 2, 3, 4, 11]:
x[i] = num.random.uniform(xbounds[i, 0], xbounds[i, 1])
x[5:11] = mtm.random_m6()
return x.tolist()
def preconstrain(self, x):
source = self.unpack(x)
if any(self.distance_min > source.distance_to(t)
for t in self.targets):
raise core.Forbidden()
return x
def bounds(self):
out = []
for p in self.parameters:
r = self.ranges[p.name]
out.append((r.start, r.stop))
return out
def dependant_bounds(self):
print('dependent_bounds')
out = [
(-1., 1.),
(-1., 1.)]
return out
def evaluate(self, x, result_mode='sparse', mask=None):
source = self.unpack(x)
engine = self.get_engine()
for target in self.targets:
target.set_result_mode(result_mode)
if mask is not None:
assert len(mask) == len(self.targets)
targets_ok = [
target for (target, ok) in zip(self.targets, mask) if ok]
else:
targets_ok = self.targets
resp = engine.process(source, targets_ok)
if mask is not None:
ires_ok = 0
results = []
for target, ok in zip(self.targets, mask):
if ok:
results.append(resp.results_list[0][ires_ok])
ires_ok += 1
else:
results.append(
gf.SeismosizerError(
'skipped because of previous failure'))
else:
results = list(resp.results_list[0])
data = []
for target, result in zip(self.targets, results):
if isinstance(result, gf.SeismosizerError):
logger.debug(
'%s.%s.%s.%s: %s' % (target.codes + (str(result),)))
data.append((None, None))
else:
data.append((result.misfit_value, result.misfit_norm))
ms, ns = num.array(data, dtype=num.float).T
if result_mode == 'full':
return ms, ns, results
else:
return ms, ns
def forward(self, x):
source = self.unpack(x)
engine = self.get_engine()
plain_targets = [target.get_plain_target() for target in self.targets]
resp = engine.process(source, plain_targets)
results = []
for target, result in zip(self.targets, resp.results_list[0]):
if isinstance(result, gf.SeismosizerError):
logger.debug(
'%s.%s.%s.%s: %s' % (target.codes + (str(result),)))
results.append(None)
else:
results.append(result)
return results
def get_target_weights(self):
if self._target_weights is None:
self._target_weights = num.array(
[target.get_combined_weight(
apply_balancing_weights=self.apply_balancing_weights)
for target in self.targets], dtype=num.float)
return self._target_weights
def inter_group_weights(self, ns):
group, ngroups = self.get_group_mask()
ws = num.zeros(self.ntargets)
for igroup in xrange(ngroups):
mask = group == igroup
ws[mask] = 1. / num.sqrt(num.nansum(ns[mask]**2))
return ws
def inter_group_weights2(self, ns):
group, ngroups = self.get_group_mask()
ws = num.zeros(ns.shape)
for igroup in xrange(ngroups):
mask = group == igroup
ws[:, mask] = (1. / num.sqrt(
num.nansum(ns[:, mask]**2, axis=1)))[:, num.newaxis]
return ws
def bootstrap_misfit(self, ms, ns, ibootstrap=None):
w = self.get_bootstrap_weights(ibootstrap) * \
self.get_target_weights() * self.inter_group_weights(ns)
if ibootstrap is None:
return num.sqrt(
num.nansum((w*ms[num.newaxis, :])**2, axis=1) /
num.nansum((w*ns[num.newaxis, :])**2, axis=1))
else:
return num.sqrt(num.nansum((w*ms)**2) / num.nansum((w*ns)**2))
def bootstrap_misfits(self, misfits, ibootstrap):
w = self.get_bootstrap_weights(ibootstrap)[num.newaxis, :] * \
self.get_target_weights()[num.newaxis, :] * \
self.inter_group_weights2(misfits[:, :, 1])
bms = num.sqrt(num.nansum((w*misfits[:, :, 0])**2, axis=1) /
num.nansum((w*misfits[:, :, 1])**2, axis=1))
return bms
def global_misfit(self, ms, ns):
ws = self.get_target_weights() * self.inter_group_weights(ns)
m = num.sqrt(num.nansum((ws*ms)**2) / num.nansum((ws*ns)**2))
return m
def global_misfits(self, misfits):
ws = self.get_target_weights()[num.newaxis, :] * \
self.inter_group_weights2(misfits[:, :, 1])
gms = num.sqrt(num.nansum((ws*misfits[:, :, 0])**2, axis=1) /
num.nansum((ws*misfits[:, :, 1])**2, axis=1))
return gms
def global_contributions(self, misfits):
ws = self.get_target_weights()[num.newaxis, :] * \
self.inter_group_weights2(misfits[:, :, 1])
gcms = (ws*misfits[:, :, 0])**2 / \
num.nansum((ws*misfits[:, :, 1])**2, axis=1)[:, num.newaxis]
return gcms
class RectangularProblemConfig(core.ProblemConfig):
ranges = Dict.T(String.T(), gf.Range.T())
name = String.T()
distance_min = Float.T(default=0.0)
nbootstrap = Int.T(default=0)
apply_balancing_weights = False
def get_problem(self, rect_source, targets):
base_source = gf.RecangularSource.load(rect_source)
problem = RectangularProblem(
name=core.expand_template(self.name_template, self.name),
apply_balancing_weights=self.apply_balancing_weights,
base_source=base_source,
targets=targets,
ranges=self.ranges,
distance_min=self.distance_min,
nbootstrap=self.nbootstrap,
mt_type=self.mt_type)
return problem
__all__ = [
'RectangularProblem',
'RectangularProblemConfig',
]
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment