Skip to content

GitLab

Commit 38bb0322 authored by Marius Isken's avatar Marius Isken
Browse files

refactoring problems

parent c54c953e
Hide whitespace changes
Inline Side-by-side
setup.py
View file @ 38bb0322
......@@ -15,7 +15,7 @@ setup(
version='0.1',
author='Sebastian Heimann',
author_email='sebastian.heimann@gfz-potsdam.de',
packages=['grond', 'grond.baraddur'],
packages=['grond', 'grond.baraddur', 'grond.problems'],
scripts=['apps/grond'],
package_dir={'grond': 'src'},
package_data={'grond': ['baraddur/templates/*.html',
......
src/core.py
View file @ 38bb0322
......@@ -14,7 +14,7 @@ from pyrocko import orthodrome as od, gf, trace, guts, util, weeding
from pyrocko import parimap, model, marker as pmarker
from .dataset import DatasetConfig, NotFound
from .problems import ProblemConfig, Problem
from .problems.base import ProblemConfig, Problem
from .targets import TargetAnalysisResult, TargetConfig
from .meta import (Path, HasPaths, expand_template, xjoin, GrondError,
Forbidden)
......
src/problems/__init__.py 0 → 100644
View file @ 38bb0322
from .cmt import CMTProblem, CMTProblemConfig # noqa
from .rectangular import RectangularProblem, RectangularProblemConfig # noqa
from .double_dc import DoubleDCProblem, DoubleDCProblemConfig # noqa
src/problems.py → src/problems/base.py
View file @ 38bb0322
......@@ -2,14 +2,12 @@ import numpy as num
import copy
import logging
import os.path as op
import math
from pyrocko import gf, util, guts, moment_tensor as mtm
from pyrocko.guts import (Object, String, Bool, List, Float, Dict, Int,
StringChoice)
from pyrocko import gf, util, guts
from pyrocko.guts import Object, String, Bool, List, Dict, Int
from .meta import Forbidden, expand_template, ADict, Parameter, GrondError
from .targets import MisfitTarget, MisfitSatelliteTarget
from ..meta import ADict, Parameter, GrondError
from ..targets import MisfitTarget, MisfitSatelliteTarget
guts_prefix = 'grond'
......@@ -360,366 +358,6 @@ class Problem(Object):
def evaluate(self, x, dump_data=False):
raise NotImplementedError()
class CMTProblem(Problem):
problem_parameters = [
Parameter('time', 's', label='Time'),
Parameter('north_shift', 'm', label='Northing', **as_km),
Parameter('east_shift', 'm', label='Easting', **as_km),
Parameter('depth', 'm', label='Depth', **as_km),
Parameter('magnitude', label='Magnitude'),
Parameter('rmnn', label='$m_{nn} / M_0$'),
Parameter('rmee', label='$m_{ee} / M_0$'),
Parameter('rmdd', label='$m_{dd} / M_0$'),
Parameter('rmne', label='$m_{ne} / M_0$'),
Parameter('rmnd', label='$m_{nd} / M_0$'),
Parameter('rmed', label='$m_{ed} / M_0$'),
Parameter('duration', 's', label='Duration')]
dependants = [
Parameter('rel_moment_iso', label='$M_{0}^{ISO}/M_{0}$'),
Parameter('rel_moment_clvd', label='$M_{0}^{CLVD}/M_{0}$')]
distance_min = Float.T(default=0.0)
nbootstrap = Int.T(default=10)
mt_type = StringChoice.T(
default='full', choices=['full', 'deviatoric', 'dc'])
def get_source(self, x):
d = self.get_parameter_dict(x)
rm6 = num.array([d.rmnn, d.rmee, d.rmdd, d.rmne, d.rmnd, d.rmed],
dtype=num.float)
m0 = mtm.magnitude_to_moment(d.magnitude)
m6 = rm6 * m0
p = {}
for k in self.base_source.keys():
if k in d:
p[k] = float(
self.ranges[k].make_relative(self.base_source[k], d[k]))
stf = gf.HalfSinusoidSTF(duration=float(d.duration))
source = self.base_source.clone(m6=m6, stf=stf, **p)
return source
def make_dependant(self, xs, pname):
if xs.ndim == 1:
return self.make_dependant(xs[num.newaxis, :], pname)[0]
if pname in ('rel_moment_iso', 'rel_moment_clvd'):
y = num.zeros(xs.shape[0])
for i, x in enumerate(xs):
source = self.get_source(x)
mt = source.pyrocko_moment_tensor()
res = mt.standard_decomposition()
if pname == 'rel_moment_iso':
ratio_iso, m_iso = res[0][1:3]
y[i] = ratio_iso * num.sign(m_iso[0, 0])
else:
ratio_clvd, m_clvd = res[2][1:3]
evals, evecs = mtm.eigh_check(m_clvd)
ii = num.argmax(num.abs(evals))
y[i] = ratio_clvd * num.sign(evals[ii])
return y
else:
raise KeyError(pname)
def extract(self, xs, i):
if xs.ndim == 1:
return self.extract(xs[num.newaxis, :], i)[0]
if i < self.nparameters:
return xs[:, i]
else:
return self.make_dependant(
xs, self.dependants[i-self.nparameters].name)
def pack(self, source):
m6 = source.m6
mt = source.pyrocko_moment_tensor()
rm6 = m6 / mt.scalar_moment()
x = num.array([
source.time - self.base_source.time,
source.north_shift,
source.east_shift,
source.depth,
mt.moment_magnitude(),
] + rm6.tolist() + [source.stf.duration], dtype=num.float)
return x
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):
d = self.get_parameter_dict(x)
m6 = num.array([d.rmnn, d.rmee, d.rmdd, d.rmne, d.rmnd, d.rmed],
dtype=num.float)
m9 = mtm.symmat6(*m6)
if self.mt_type == 'deviatoric':
trace_m = num.trace(m9)
m_iso = num.diag([trace_m / 3., trace_m / 3., trace_m / 3.])
m9 -= m_iso
elif self.mt_type == 'dc':
mt = mtm.MomentTensor(m=m9)
m9 = mt.standard_decomposition()[1][2]
m0_unscaled = math.sqrt(num.sum(m9.A**2)) / math.sqrt(2.)
m9 /= m0_unscaled
m6 = mtm.to6(m9)
d.rmnn, d.rmee, d.rmdd, d.rmne, d.rmnd, d.rmed = m6
x = self.get_parameter_array(d)
source = self.get_source(x)
if any(self.distance_min > source.distance_to(t)
for t in self.targets):
raise Forbidden()
return x
def get_dependant_bounds(self):
out = [
(-1., 1.),
(-1., 1.)]
return out
def evaluate(self, x, result_mode='sparse', mask=None):
source = self.get_source(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.get_source(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
class CMTProblemConfig(ProblemConfig):
ranges = Dict.T(String.T(), gf.Range.T())
distance_min = Float.T(default=0.0)
nbootstrap = Int.T(default=10)
mt_type = StringChoice.T(choices=['full', 'deviatoric'])
def get_problem(self, event, targets):
if event.depth is None:
event.depth = 0.
base_source = gf.MTSource.from_pyrocko_event(event)
base_source.stf = gf.HalfSinusoidSTF(duration=event.duration or 0.0)
subs = dict(
event_name=event.name,
event_time=util.time_to_str(event.time))
problem = CMTProblem(
name=expand_template(self.name_template, subs),
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,
norm_exponent=self.norm_exponent)
return problem
class DoubleDCProblem(Problem):
problem_parameters = [
Parameter('time', 's', label='Time'),
Parameter('north_shift', 'm', label='Northing', **as_km),
Parameter('east_shift', 'm', label='Easting', **as_km),
Parameter('depth', 'm', label='Depth', **as_km),
Parameter('magnitude', label='Magnitude'),
Parameter('strike1', 'deg', label='Strike 1'),
Parameter('dip1', 'deg', label='Dip 1'),
Parameter('rake1', 'deg', label='Rake 1'),
Parameter('strike2', 'deg', label='Strike 2'),
Parameter('dip2', 'deg', label='Dip 2'),
Parameter('rake2', 'deg', label='Rake 2'),
Parameter('delta_time', 's', label='$\Delta$ Time'),
Parameter('delta_depth', 'm', label='$\Delta$ Depth'),
Parameter('azimuth', 'deg', label='Azimuth'),
Parameter('distance', 'm', label='Distance'),
Parameter('mix', label='Mix'),
Parameter('duration1', 's', label='Duration 1'),
Parameter('duration2', 's', label='Duration 2')]
dependants = []
distance_min = Float.T(default=0.0)
nbootstrap = Int.T(default=100)
def get_source(self, x):
d = self.get_parameter_dict(x)
p = {}
for k in self.base_source.keys():
if k in d:
p[k] = float(
self.ranges[k].make_relative(self.base_source[k], d[k]))
stf1 = gf.HalfSinusoidSTF(duration=float(d.duration1))
stf2 = gf.HalfSinusoidSTF(duration=float(d.duration2))
source = self.base_source.clone(stf1=stf1, stf2=stf2, **p)
return source
def make_dependant(self, xs, pname):
if xs.ndim == 1:
return self.make_dependant(xs[num.newaxis, :], pname)[0]
raise KeyError(pname)
def extract(self, xs, i):
if xs.ndim == 1:
return self.extract(xs[num.newaxis, :], i)[0]
if i < self.nparameters:
return xs[:, i]
else:
return self.make_dependant(
xs, self.dependants[i-self.nparameters].name)
def pack(self, source):
x = num.array([
source.time - self.base_source.time,
source.north_shift,
source.east_shift,
source.depth,
source.magnitude,
source.strike1,
source.dip1,
source.rake1,
source.strike2,
source.dip2,
source.rake2,
source.delta_time,
source.delta_depth,
source.azimuth,
source.distance,
source.mix,
source.stf1.duration,
source.stf2.duration], dtype=num.float)
return x
def random_uniform(self, xbounds):
x = num.zeros(self.nparameters)
for i in xrange(self.nparameters):
x[i] = num.random.uniform(xbounds[i, 0], xbounds[i, 1])
return x.tolist()
def preconstrain(self, x):
source = self.get_source(x)
if any(self.distance_min > source.distance_to(t)
for t in self.targets):
raise Forbidden()
return num.array(x, dtype=num.float)
def evaluate(self, x, result_mode='sparse'):
source = self.get_source(x)
engine = self.get_engine()
for target in self.targets:
target.set_result_mode(result_mode)
resp = engine.process(source, self.targets)
data = []
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),)))
data.append((None, None))
results.append(result)
else:
data.append((result.misfit_value, result.misfit_norm))
results.append(result)
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.get_source(x)
engine = self.get_engine()
......@@ -737,202 +375,3 @@ class DoubleDCProblem(Problem):
results.append(result)
return results
class DoubleDCProblemConfig(ProblemConfig):
ranges = Dict.T(String.T(), gf.Range.T())
distance_min = Float.T(default=0.0)
nbootstrap = Int.T(default=100)
def get_problem(self, event, targets):
if event.depth is None:
event.depth = 0.
base_source = gf.DoubleDCSource.from_pyrocko_event(event)
base_source.stf = gf.HalfSinusoidSTF(duration=event.duration or 0.0)
subs = dict(
event_name=event.name,
event_time=util.time_to_str(event.time))
problem = DoubleDCProblem(
name=expand_template(self.name_template, subs),
apply_balancing_weights=self.apply_balancing_weights,
base_source=base_source,
targets=targets,
ranges=self.ranges,
distance_min=self.distance_min,
nbootstrap=self.nbootstrap,
norm_exponent=self.norm_exponent)
return problem
class RectangularProblem(Problem):
# nucleation_x
# nucleation_y
# time
# stf
problem_parameters = [
Parameter('north_shift', 'm', label='Northing', **as_km),
Parameter('east_shift', 'm', label='Easting', **as_km),
Parameter('depth', 'm', label='Depth', **as_km),
Parameter('length', 'm', label='Length', **as_km),
Parameter('width', 'm', label='Width', **as_km),
Parameter('dip', 'deg', label='Dip'),
Parameter('strike', 'deg', label='Strike'),
Parameter('rake', 'deg', label='Rake'),
Parameter('slip', 'm', label='Slip'),
]
problem_waveform_parameters = [
Parameter('nucleation_x', 'offset', label='Nucleation X'),
Parameter('nucleation_y', 'offset', label='Nucleation Y'),
Parameter('time', 's', label='Time'),
]
dependants = []
distance_min = Float.T(default=0.0)
nbootstrap = Int.T(default=10)
def pack(self, source):
arr = self.get_parameter_array(source)
for ip, p in enumerate(self.parameters):
if p.name == 'time':
arr[ip] -= self.base_source.time
return arr
def get_source(self, x):
d = self.get_parameter_dict(x)
if 'time' in d.keys():
d['time'] += self.base_source['time']
source = self.base_source.clone(**d)
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 xrange(self.nparameters):
x[i] = num.random.uniform(xbounds[i, 0], xbounds[i, 1])
return x.tolist()
def preconstrain(self, x):
# source = self.get_source(x)
# if any(self.distance_min > source.distance_to(t)
# for t in self.targets):
# raise Forbidden()
return x
def evaluate(self, x, result_mode='sparse', mask=None, nprocs=0):
source = self.get_source(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
self.set_target_parameter_values(x)
resp = engine.process(source, targets_ok, nprocs=nprocs)
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, 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, nprocs=nprocs)
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
class RectangularProblemConfig(ProblemConfig):
ranges = Dict.T(String.T(), gf.Range.T())
apply_balancing_weights = Bool.T(default=False)
decimation_factor = Int.T(default=1)
nbootstrap = Int.T(default=10)
distance_min = Float.T(default=0.)
def get_problem(self, event, targets):
base_source = gf.RectangularSource(
lat=event.lat,
lon=event.lon,
time=event.time,
depth=event.depth,
anchor='top',
decimation_factor=self.decimation_factor,
)
problem = RectangularProblem(
name=expand_template(self.name_template, event.name),
apply_balancing_weights=self.apply_balancing_weights,
base_source=base_source,
nbootstrap=self.nbootstrap,
distance_min=self.distance_min,
targets=targets,
ranges=self.ranges,
norm_exponent=self.norm_exponent)