Commit 88fa2aaa authored by Marius Isken's avatar Marius Isken
Browse files

TagetParameters // Leveling plane params for SatelliteMisfitTarget

parent c3cf0951
......@@ -13,7 +13,7 @@ from pyrocko.guts import (load, Object, String, Float, Int, Bool, List,
from pyrocko import orthodrome as od, gf, trace, guts, util
from pyrocko import parimap, model, marker as pmarker
from .dataset import DatasetConfig
from .dataset import DatasetConfig, NotFound
from .problems import ProblemConfig, Problem
from .targets import TargetAnalysisResult, TargetConfig
from .meta import (Path, HasPaths, expand_template, xjoin, GrondError,
......@@ -249,7 +249,10 @@ def analyse(problem, niter=1000, show_progress=False):
return
wtargets = []
for target in problem.targets:
if not problem.has_waveforms:
return
for target in problem.waveform_targets:
wtarget = copy.copy(target)
wtarget.flip_norm = True
wtarget.weight = 1.0
......@@ -317,14 +320,14 @@ def excentricity_compensated_choice(xs, sbx, factor):
inonflat = num.where(sbx != 0.0)[0]
scale = num.zeros_like(sbx)
scale[inonflat] = 1.0 / (sbx[inonflat] * (factor if factor != 0. else 1.0))
#distances_all = math.sqrt(num.sum(
# ((xs[num.newaxis, :, :] - xs[:, num.newaxis, :]) *
# scale[num.newaxis, num.newaxis, :])**2, axis=2))
# distances_all = math.sqrt(num.sum(
# ((xs[num.newaxis, :, :] - xs[:, num.newaxis, :]) *
# scale[num.newaxis, num.newaxis, :])**2, axis=2))
distances_sqr_all = num.sum(
((xs[num.newaxis, :, :] - xs[:, num.newaxis, :]) *
scale[num.newaxis, num.newaxis, :])**2, axis=2)
probabilities = 1.0 / num.sum(distances_sqr_all < 1.0, axis=1)
print num.sort(num.sum(distances_sqr_all < 1.0, axis=1))
# print num.sort(num.sum(distances_sqr_all < 1.0, axis=1))
probabilities /= num.sum(probabilities)
r = num.random.random()
ichoice = num.searchsorted(num.cumsum(probabilities), r)
......@@ -339,7 +342,7 @@ def select_most_excentric(xcandidates, xs, sbx, factor):
distances_sqr_all = num.sum(
((xcandidates[num.newaxis, :, :] - xs[:, num.newaxis, :]) *
scale[num.newaxis, num.newaxis, :])**2, axis=2)
#print num.sort(num.sum(distances_sqr_all < 1.0, axis=0))
# print num.sort(num.sum(distances_sqr_all < 1.0, axis=0))
ichoice = num.argmin(num.sum(distances_sqr_all < 1.0, axis=0))
return xcandidates[ichoice]
......@@ -355,11 +358,11 @@ def solve(problem,
chain_length_factor=8.0,
xs_inject=None,
sampler_distribution='multivariate_normal',
compensate_excentricity=True,
compensate_excentricity=False,
status=()):
xbounds = num.array(problem.get_parameter_bounds(), dtype=num.float)
npar = xbounds.shape[0]
npar = problem.nparameters
nlinks_cap = int(round(chain_length_factor * npar + 1))
chains_m = num.zeros((1 + problem.nbootstrap, nlinks_cap), num.float)
......@@ -383,7 +386,7 @@ def solve(problem,
isbad_mask = None
accept_sum = num.zeros(1 + problem.nbootstrap, dtype=num.int)
accept_hist = num.zeros(niter, dtype=num.int)
pnames = [p.name for p in problem.parameters]
pnames = problem.parameter_names
while iiter < niter:
......@@ -832,7 +835,8 @@ def check(
check_problem(problem)
xbounds = num.array(problem.get_parameter_bounds(), dtype=num.float)
xbounds = num.array(problem.get_parameter_bounds(),
dtype=num.float)
results_list = []
......@@ -907,9 +911,10 @@ def check(
tfade=tfade,
freqlimits=freqlimits,
deltat=deltat,
backazimuth=target.get_backazimuth_for_waveform(),
backazimuth=target.
get_backazimuth_for_waveform(),
debug=True)
except dataset.NotFound, e:
except NotFound, e:
logger.warn(str(e))
continue
......
......@@ -51,7 +51,7 @@ class ADict(dict):
class Parameter(Object):
name = String.T()
name__ = String.T()
unit = String.T(optional=True)
scale_factor = Float.T(default=1., optional=True)
scale_unit = String.T(optional=True)
......@@ -63,7 +63,8 @@ class Parameter(Object):
if len(args) >= 2:
kwargs['unit'] = args[1]
self.target = kwargs.pop('target', None)
self.groups = [None]
self._name = None
Object.__init__(self, **kwargs)
......@@ -76,6 +77,25 @@ class Parameter(Object):
return ' '.join(l)
def set_groups(self, groups):
if not isinstance(groups, list):
raise AttributeError('Groups must be a list of strings.')
self.groups = groups
def _get_name(self):
if None not in self.groups:
return '%s:%s' % (':'.join(self.groups), self._name)
return self._name
def _set_name(self, value):
self._name = value
name = property(_get_name, _set_name)
@property
def name_nogroups(self):
return self._name
def get_value_label(self, value, format='%(value)g%(unit)s'):
value = self.scaled(value)
unit = self.get_unit_suffix()
......
......@@ -33,9 +33,7 @@ class Problem(Object):
targets = List.T()
def __init__(self, **kwargs):
Object.__init__(self, **kwargs)
self.init_satellite_target_leveling()
self._bootstrap_weights = None
self._target_weights = None
......@@ -52,14 +50,6 @@ class Problem(Object):
o._target_weights = None
return o
@property
def parameters(self):
ps = self.problem_parameters
for target in self.targets:
ps.extend(target.target_parameters)
return ps
def set_target_parameter_values(self, x):
i = len(self.problem_parameters)
for target in self.targets:
......@@ -67,16 +57,16 @@ class Problem(Object):
target.set_parameter_values(x[i:i+n])
i += n
def get_parameter_dict(self, x):
return ADict((p.name, v) for p, v in zip(self.parameters, x))
def get_parameter_dict(self, model, group=None):
params = []
for ip, p in enumerate(self.parameters):
if group in p.groups:
params.append((p.name, model[ip]))
return ADict(params)
def get_parameter_array(self, d):
return num.array([d[p.name] for p in self.parameters], dtype=num.float)
@property
def parameter_names(self):
return [p.name for p in self.combined]
def dump_problem_info(self, dirname):
fn = op.join(dirname, 'problem.yaml')
util.ensuredirs(fn)
......@@ -98,6 +88,17 @@ class Problem(Object):
pnames = [p.name for p in self.combined]
return pnames.index(name)
@property
def parameters(self):
target_parameters = []
for target in self.targets:
target_parameters.extend(target.target_parameters)
return self.problem_parameters + target_parameters
@property
def parameter_names(self):
return [p.name for p in self.combined]
@property
def nparameters(self):
return len(self.parameters)
......@@ -143,6 +144,9 @@ class Problem(Object):
def set_engine(self, engine):
self._engine = engine
def random_uniform(self, xbounds):
raise NotImplementedError()
def make_bootstrap_weights(self, nbootstrap):
ntargets = len(self.targets)
ws = num.zeros((nbootstrap, ntargets))
......@@ -204,9 +208,8 @@ class Problem(Object):
for target in self.targets:
for p in target.target_parameters:
r = target.ranges[p.name]
r = target.target_ranges[p.name]
out.append((r.start, r.stop))
return out
def get_dependant_bounds(self):
......@@ -678,6 +681,11 @@ class DoubleDCProblemConfig(ProblemConfig):
class RectangularProblem(Problem):
# nucleation_x
# nucleation_y
# rupture_velocity
# t0
# stf
problem_parameters = [
Parameter('north_shift', 'm', label='Northing', **as_km),
......@@ -699,7 +707,8 @@ class RectangularProblem(Problem):
return self.get_parameter_array(source)
def get_source(self, x):
source = self.base_source.clone(**self.get_parameter_dict(x))
source = self.base_source.clone(
**self.get_parameter_dict(x, group=None))
return source
def extract(self, xs, i):
......@@ -729,8 +738,6 @@ class RectangularProblem(Problem):
for target in self.targets:
target.set_result_mode(result_mode)
self.set_satellite_scene_levels(x)
if mask is not None:
assert len(mask) == len(self.targets)
targets_ok = [
......@@ -738,6 +745,7 @@ class RectangularProblem(Problem):
else:
targets_ok = self.targets
self.set_target_parameter_values(x)
resp = engine.process(source, targets_ok, nprocs=nprocs)
if mask is not None:
......@@ -770,12 +778,12 @@ class RectangularProblem(Problem):
else:
return ms, ns
def forward(self, x):
def forward(self, x, nprocs=0):
source = self.get_source(x)
engine = self.get_engine()
plain_targets = [target.get_plain_target() for target in self.targets]
resp = engine.process(source, plain_targets)
resp = engine.process(source, plain_targets, nprocs=nprocs)
results = []
for target, result in zip(self.targets, resp.results_list[0]):
if isinstance(result, gf.SeismosizerError):
......
import logging
import math
import numpy as num
import copy
from pyrocko import gf, trace, weeding
from pyrocko.guts import (Object, String, Float, Bool, Int, List, Dict,
......@@ -8,6 +9,7 @@ from pyrocko.guts import (Object, String, Float, Bool, Int, List, Dict,
from pyrocko.guts_array import Array
from .dataset import NotFound
from .meta import Parameter
guts_prefix = 'grond'
logger = logging.getLogger('grond.target')
......@@ -57,6 +59,50 @@ class TraceSpectrum(Object):
return self.fmin + num.arange(self.ydata.size) * self.deltaf
class GrondTarget(object):
def set_dataset(self, ds):
self._ds = ds
def get_dataset(self):
return self._ds
@property
def id(self):
return self.scene_id
@property
def nparameters(self):
return len(self._target_parameters)
@property
def target_parameters(self):
if self._target_parameters is None:
self._target_parameters = copy.copy(self.parameters)
for p in self._target_parameters:
p.set_groups([self.id])
return self._target_parameters
@property
def target_ranges(self):
if self._target_ranges is None:
self._target_ranges = self.inner_misfit_config.ranges.copy()
for k in self._target_ranges.keys():
self._target_ranges['%s:%s' % (self.id, k)] =\
self._target_ranges.pop(k)
return self._target_ranges
def set_parameter_values(self, model):
for i, p in enumerate(self.parameters):
self.parameter_values[p.name_nogroups] = model[i]
def set_result_mode(self, result_mode):
pass
def string_id(self):
return '.'.join([self.super_group, self.group, self.id])
class InnerMisfitConfig(Object):
fmin = Float.T()
fmax = Float.T()
......@@ -95,6 +141,8 @@ class InnerMisfitConfig(Object):
'``autoshift_penalty_max * normalization_factor * tautoshift**2 '
'/ tautoshift_max**2``')
ranges = {}
def get_full_frequency_range(self):
return self.fmin / self.ffactor, self.fmax * self.ffactor
......@@ -115,46 +163,37 @@ class MisfitResult(gf.Result):
cc = Trace.T(optional=True)
class MisfitTarget(gf.Target):
class MisfitTarget(gf.Target, GrondTarget):
misfit_config = InnerMisfitConfig.T()
flip_norm = Bool.T(default=False)
manual_weight = Float.T(default=1.0)
analysis_result = TargetAnalysisResult.T(optional=True)
super_group = gf.StringID.T()
group = gf.StringID.T()
analysis_result = TargetAnalysisResult.T(optional=True)
parameters = []
def __init__(self, **kwargs):
gf.Target.__init__(self, **kwargs)
self._ds = None
self._result_mode = 'sparse'
self._target_parameters = None
self._target_ranges = None
def string_id(self):
return '.'.join(x for x in (
self.super_group, self.group) + self.codes if x)
@property
def id(self):
return self.codes
def get_plain_target(self):
d = dict(
(k, getattr(self, k)) for k in gf.Target.T.propnames)
return gf.Target(**d)
def get_dataset(self):
return self._ds
def set_dataset(self, ds):
self._ds = ds
def set_result_mode(self, result_mode):
self._result_mode = result_mode
def get_combined_weight(self, apply_balancing_weights):
w = self.manual_weight
if apply_balancing_weights:
w *= self.get_balancing_weight()
return w
def get_balancing_weight(self):
......@@ -468,53 +507,36 @@ def _process(tr, tmin, tmax, taper, domain):
class InnerSatelliteMisfitConfig(Object):
use_weight_focal = Bool.T(default=False)
optimize_leveling_plane = Bool.T(default=True)
ranges = Dict.T(String.T(), gf.Range.T(),
default={'waterlevel': '-0.5 .. 0.5',
'ramp_north': '-1e-4 .. 1e-4',
'ramp_east': '-1e-4 .. 1e-4'})
default={'waterlevel': '-0.5 .. 0.5',
'ramp_north': '-1e-4 .. 1e-4',
'ramp_east': '-1e-4 .. 1e-4'})
class MisfitSatelliteTarget(gf.SatelliteTarget):
class MisfitSatelliteTarget(gf.SatelliteTarget, GrondTarget):
scene_id = String.T()
super_group = gf.StringID.T()
inner_misfit_config = InnerSatelliteMisfitConfig.T()
manual_weight = Float.T(default=1.0)
group = gf.StringID.T()
target_parameters = [
...]
ranges = Dict.T(String.T(), gf.Range.T())
parameters = [
Parameter('waterlevel', 'm'),
Parameter('ramp_north', 'm/m'),
Parameter('ramp_east', 'm/m'),
]
def __init__(self, *args, **kwargs):
gf.SatelliteTarget.__init__(self, *args, **kwargs)
self.waterlevel = 0.
self.ramp_north = 0.
self.ramp_east = 0.
self._ds = None
self._distance_cache = None
def set_dataset(self, ds):
self._ds = ds
@property
def id(self):
return self.scene_id
def get_dataset(self):
return self._ds
if not self.inner_misfit_config.optimize_leveling_plane:
self.parameters = []
def set_result_mode(self, result_mode):
pass
def string_id(self):
return '.'.join([self.super_group, self.group, self.scene_id])
self._ds = None
def set_scene_levels(self, waterlevel, ramp_north, ramp_east):
self.waterlevel = waterlevel
self.ramp_north = ramp_north
self.ramp_east = ramp_east
self.parameter_values = {}
self._target_parameters = None
self._target_ranges = None
def post_process(self, engine, source, statics):
scene = self._ds.get_kite_scene(self.scene_id)
......@@ -523,12 +545,16 @@ class MisfitSatelliteTarget(gf.SatelliteTarget):
stat_obs = scene.quadtree.leaf_medians
stat_syn = statics['displacement.los']
stat_level = num.zeros_like(stat_obs)
stat_level.fill(self.waterlevel)
stat_level += (quadtree.leaf_center_distance[:, 0] * self.ramp_east)
stat_level += (quadtree.leaf_center_distance[:, 1] * self.ramp_north)
if self.inner_misfit_config.optimize_leveling_plane:
stat_level = num.zeros_like(stat_obs)
stat_level.fill(self.parameter_values['waterlevel'])
stat_level += (quadtree.leaf_center_distance[:, 0]
* self.parameter_values['ramp_east'])
stat_level += (quadtree.leaf_center_distance[:, 1]
* self.parameter_values['ramp_north'])
stat_syn += stat_level
res = num.abs(stat_obs - (stat_syn + stat_level))
res = num.abs(stat_obs - stat_syn)
misfit_value = num.sqrt(
num.sum((res * scene.covariance.weight_vector)**2))
......@@ -542,11 +568,6 @@ class MisfitSatelliteTarget(gf.SatelliteTarget):
return 1.
class GrondTarget(MisfitSatelliteTarget, MisfitTarget):
def __init__(self, args):
super(MisfitSatelliteTarget, self).__init__()
class TargetConfig(Object):
super_group = gf.StringID.T(default='', optional=True)
......
Supports Markdown
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