Commit ccbb734c authored by Sebastian Heimann's avatar Sebastian Heimann

more flexibility in modelling

- MisfitTarget and gf.Target are now conceptually independent
- MisfitTarget can produce more than one misfit (impl. incomplete)
parent 975d5b99
......@@ -6,7 +6,7 @@ license-file=LICENSE
universal=1
[build_ext]
inplace=1
inplace=0
[nosetests]
verbosity=2
......
......@@ -19,6 +19,7 @@ from .problems.base import ProblemConfig, Problem, \
from .optimizers.base import OptimizerConfig, BadProblem
from .targets.base import TargetGroup
from .targets.waveform.target import WaveformMisfitResult
from .analysers.base import AnalyserConfig
from .meta import Path, HasPaths, expand_template, GrondError, Forbidden
......@@ -130,7 +131,8 @@ class Config(HasPaths):
def get_problem(self, event):
targets = self.get_targets(event)
problem = self.problem_config.get_problem(event, targets)
problem = self.problem_config.get_problem(
event, self.target_groups, targets)
self.setup_modelling_environment(problem)
return problem
......@@ -254,7 +256,7 @@ def forward(rundir_or_config_path, event_names):
events.append(event)
for result in results:
if not isinstance(result, gf.SeismosizerError):
if isinstance(result, WaveformMisfitResult):
result.filtered_obs.set_codes(location='ob')
result.filtered_syn.set_codes(location='sy')
all_trs.append(result.filtered_obs)
......@@ -452,6 +454,7 @@ def check(
backazimuth=target.
get_backazimuth_for_waveform(),
debug=True)
except NotFound as e:
logger.warn(str(e))
continue
......@@ -498,7 +501,7 @@ def check(
yabsmaxs = []
for results in results_list:
result = results[itarget]
if not isinstance(result, gf.SeismosizerError):
if isinstance(result, WaveformMisfitResult):
yabsmaxs.append(
num.max(num.abs(
result.filtered_obs.get_ydata())))
......@@ -512,7 +515,7 @@ def check(
ii = 0
for results in results_list:
result = results[itarget]
if not isinstance(result, gf.SeismosizerError):
if isinstance(result, WaveformMisfitResult):
if fig is None:
fig = plt.figure()
axes = fig.add_subplot(1, 1, 1)
......@@ -521,7 +524,8 @@ def check(
xdata = result.filtered_obs.get_xdata()
ydata = result.filtered_obs.get_ydata() / yabsmax
axes.plot(xdata, ydata*0.5 + 3.5, color='black')
axes.plot(
xdata, ydata*0.5 + 3.5, color='black')
color = plot.mpl_graph_color(ii)
......
......@@ -442,8 +442,8 @@ class HighScoreOptimizer(Optimizer):
isbad_mask != isbad_mask_new):
errmess = [
'problem %s: inconsistency in data availability' %
problem.name]
'problem %s: inconsistency in data availability at iteration %i' %
(problem.name, iiter)]
for target, isbad_new, isbad in zip(
problem.targets, isbad_mask_new, isbad_mask):
......@@ -452,7 +452,7 @@ class HighScoreOptimizer(Optimizer):
errmess.append(' %s, %s -> %s' % (
target.string_id(), isbad, isbad_new))
raise BadProblem(errmess)
raise BadProblem('\n'.join(errmess))
isbad_mask = isbad_mask_new
......
......@@ -2086,6 +2086,10 @@ def plot_result(dirname, plotnames_want,
fns.extend(
save_figs(figs, plot_dirname, plotname, formats, dpi))
for fig in figs:
plt.close(fig)
if 7 != len({
'fits',
'fits_statics',
......@@ -2114,6 +2118,9 @@ def plot_result(dirname, plotnames_want,
fns.extend(
save_figs(figs, plot_dirname, plotname, formats, dpi))
for fig in figs:
plt.close(fig)
for plotname in ['jointpar', 'hudson', 'solution', 'location']:
if plotname in plotnames_want:
figs = plot_dispatch[plotname](history, optimizer, plt)
......@@ -2121,6 +2128,9 @@ def plot_result(dirname, plotnames_want,
fns.extend(
save_figs(figs, plot_dirname, plotname, formats, dpi))
for fig in figs:
plt.close(fig)
if not save:
plt.show()
......
......@@ -11,7 +11,8 @@ from pyrocko import gf, util, guts
from pyrocko.guts import Object, String, Bool, List, Dict, Int
from ..meta import ADict, Parameter, GrondError, xjoin
from ..targets import WaveformMisfitTarget, SatelliteMisfitTarget
from ..targets import MisfitTarget, TargetGroup, WaveformMisfitTarget, \
SatelliteMisfitTarget
guts_prefix = 'grond'
......@@ -37,8 +38,8 @@ class Problem(Object):
apply_balancing_weights = Bool.T(default=True)
norm_exponent = Int.T(default=2)
base_source = gf.Source.T(optional=True)
targets = List.T()
targets = List.T(MisfitTarget.T())
target_groups = List.T(TargetGroup.T())
def __init__(self, **kwargs):
Object.__init__(self, **kwargs)
......@@ -212,10 +213,10 @@ class Problem(Object):
def get_target_weights(self):
if self._target_weights is None:
self._target_weights = num.array(
self._target_weights = num.concatenate(
[target.get_combined_weight(
apply_balancing_weights=self.apply_balancing_weights)
for target in self.targets], dtype=num.float)
for target in self.targets])
return self._target_weights
......@@ -352,26 +353,47 @@ class Problem(Object):
return self._family_mask
def evaluate(self, x, dump_data=False):
raise NotImplementedError()
def forward(self, x):
def evaluate(self, x, mask=None, result_mode='sparse'):
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)
for target in self.targets:
target.set_result_mode(result_mode)
modelling_targets = []
t2m_map = {}
for itarget, target in enumerate(self.targets):
t2m_map[target] = target.prepare_modelling()
if mask is None or mask[itarget]:
modelling_targets.extend(t2m_map[target])
resp = engine.process(source, modelling_targets)
modelling_results = list(resp.results_list[0])
imt = 0
imisfit = 0
misfits = num.zeros((self.nmisfits, 2))
misfits.fill(None)
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),)))
for itarget, target in enumerate(self.targets):
nmt_this = len(t2m_map[target])
if mask is None or mask[itarget]:
misfits[imisfit:imisfit+target.nmisfits, :], result = \
target.finalize_modelling(
modelling_results[imt:imt+nmt_this])
results.append(None)
imt += nmt_this
else:
results.append(result)
result = gf.SeismosizerError(
'target was excluded from modelling')
return results
results.append(result)
imisfit += target.nmisfits
if result_mode == 'full':
return misfits, results
else:
return misfits
class InvalidRundir(Exception):
......
......@@ -21,7 +21,7 @@ class CMTProblemConfig(ProblemConfig):
distance_min = Float.T(default=0.0)
mt_type = StringChoice.T(choices=['full', 'deviatoric'])
def get_problem(self, event, targets):
def get_problem(self, event, target_groups, targets):
if event.depth is None:
event.depth = 0.
......@@ -36,6 +36,7 @@ class CMTProblemConfig(ProblemConfig):
name=expand_template(self.name_template, subs),
apply_balancing_weights=self.apply_balancing_weights,
base_source=base_source,
target_groups=target_groups,
targets=targets,
ranges=self.ranges,
distance_min=self.distance_min,
......@@ -206,70 +207,6 @@ class CMTProblem(Problem):
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))
misfits = num.array(data, dtype=num.float)
if result_mode == 'full':
return misfits, results
else:
return misfits
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
__all__ = '''
CMTProblem
......
......@@ -126,52 +126,6 @@ class DoubleDCProblem(Problem):
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()
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
__all__ = '''
DoubleDCProblem
......
......@@ -99,70 +99,6 @@ class RectangularProblem(Problem):
# 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
__all__ = '''
RectangularProblem
......
import copy
import numpy as num
from pyrocko import gf
from pyrocko.guts import Object, Float
......@@ -7,17 +9,11 @@ from pyrocko.guts import Object, Float
guts_prefix = 'grond'
class MisfitConfig(Object):
pass
class TargetGroup(Object):
normalisation_family = gf.StringID.T(optional=True)
path = gf.StringID.T(optional=True)
weight = Float.T(default=1.0)
misfit_config = MisfitConfig.T(optional=True)
interpolation = gf.InterpolationMethod.T()
store_id = gf.StringID.T(optional=True)
......@@ -32,9 +28,8 @@ class TargetAnalysisResult(Object):
balancing_weight = Float.T()
class MisfitResult(gf.Result):
misfit_value = Float.T()
misfit_norm = Float.T()
class MisfitResult(Object):
pass
class MisfitTarget(Object):
......@@ -44,19 +39,15 @@ class MisfitTarget(Object):
help='Relative weight of this target')
analysis_result = TargetAnalysisResult.T(
optional=True)
misfit_config = MisfitConfig.T(
optional=True,
help='Configuration object how the misfit is calculated.')
normalisation_family = gf.StringID.T(
optional=True,
help='Normalisation family of this misfitTarget')
help='Normalisation family of this misfit target')
path = gf.StringID.T(
help='A path identifier used for plotting')
def __init__(self, **kwargs):
Object.__init__(self, **kwargs)
self.parameters = []
self.nmisfits = 0
self._ds = None
self._result_mode = 'sparse'
......@@ -70,6 +61,10 @@ class MisfitTarget(Object):
def get_dataset(self):
return self._ds
@property
def nmisfits(self):
return 1
@property
def nparameters(self):
return len(self._target_parameters)
......@@ -84,12 +79,7 @@ class MisfitTarget(Object):
@property
def target_ranges(self):
if self._target_ranges is None:
self._target_ranges = self.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
return {}
def set_parameter_values(self, model):
for i, p in enumerate(self.parameters):
......@@ -105,12 +95,17 @@ class MisfitTarget(Object):
raise NotImplementedError()
def get_combined_weight(self, apply_balancing_weights=False):
return 1.0
return num.ones(1, dtype=num.float)
def init_modelling(self):
return []
def finalize_modelling(self, results):
raise NotImplemented('must be overloaded in subclass')
__all__ = '''
TargetGroup
MisfitConfig
MisfitTarget
MisfitResult
TargetAnalysisResult
......
import logging
from pyrocko import gf
from pyrocko.guts import Object
from ..base import MisfitTarget, MisfitConfig, MisfitResult, TargetGroup
from ..base import MisfitTarget, MisfitResult, TargetGroup
guts_prefix = 'grond'
logger = logging.getLogger('grond.targets.gnss.target')
......@@ -16,7 +17,7 @@ class GNSSMisfitResult(MisfitResult):
pass
class GNSSMisfitConfig(MisfitConfig):
class GNSSMisfitConfig(Object):
pass
......
......@@ -2,18 +2,28 @@ import logging
import numpy as num
from pyrocko import gf
from pyrocko.guts import String, Bool, Dict, List
from pyrocko.guts import String, Bool, Dict, List, Object, Float
from grond.meta import Parameter
from ..base import MisfitTarget, MisfitConfig, MisfitResult, TargetGroup
from ..base import MisfitTarget, MisfitResult, TargetGroup
guts_prefix = 'grond'
logger = logging.getLogger('grond.targets.satellite.target')
class SatelliteMisfitConfig(Object):
use_weight_focal = Bool.T(default=False)
optimize_orbital_ramp = Bool.T(default=True)
ranges = Dict.T(String.T(), gf.Range.T(),
default={'offset': '-0.5 .. 0.5',
'ramp_north': '-1e-4 .. 1e-4',
'ramp_east': '-1e-4 .. 1e-4'})
class SatelliteTargetGroup(TargetGroup):
kite_scenes = List.T(optional=True)
misfit_config = SatelliteMisfitConfig.T()
def get_targets(self, ds, event, default_path):
logger.debug('Selecting satellite targets...')
......@@ -24,11 +34,6 @@ class SatelliteTargetGroup(TargetGroup):
'*all' not in self.kite_scenes:
continue
if not isinstance(self.misfit_config,
SatelliteMisfitConfig):
raise AttributeError('misfit_config must be of type'
' SatelliteMisfitConfig')
qt = scene.quadtree
lats = num.empty(qt.nleaves)
......@@ -61,27 +66,21 @@ class SatelliteTargetGroup(TargetGroup):
return targets
class SatelliteMisfitResult(MisfitResult):
class SatelliteMisfitResult(gf.Result, MisfitResult):
misfit_value = Float.T()
misfit_norm = Float.T()
statics_syn = Dict.T(optional=True)
statics_obs = Dict.T(optional=True)
class SatelliteMisfitConfig(MisfitConfig):
use_weight_focal = Bool.T(default=False)
optimize_orbital_ramp = Bool.T(default=True)
ranges = Dict.T(String.T(), gf.Range.T(),
default={'offset': '-0.5 .. 0.5',
'ramp_north': '-1e-4 .. 1e-4',
'ramp_east': '-1e-4 .. 1e-4'})
class SatelliteMisfitTarget(MisfitTarget, gf.SatelliteTarget):
class SatelliteMisfitTarget(gf.SatelliteTarget, MisfitTarget):
scene_id = String.T()
available_parameters = [
Parameter('offset', 'm'),
Parameter('ramp_north', 'm/m'),
Parameter('ramp_east', 'm/m'),
]
misfit_config = SatelliteMisfitConfig.T()
def __init__(self, *args, **kwargs):
gf.SatelliteTarget.__init__(self, *args, **kwargs)
......@@ -93,6 +92,15 @@ class SatelliteMisfitTarget(MisfitTarget, gf.SatelliteTarget):
self.parameter_values = {}
@property
def target_ranges(self):
if self._target_ranges is None:
self._target_ranges = self.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