Browse Source

add waveform piggyback stuff

whitelist
Sebastian Heimann 2 years ago
parent
commit
b99af2f1e1
10 changed files with 200 additions and 18 deletions
  1. +1
    -0
      setup.py
  2. +21
    -4
      src/problems/base.py
  3. +1
    -0
      src/targets/__init__.py
  4. +1
    -4
      src/targets/base.py
  5. +3
    -1
      src/targets/gnss_campaign/target.py
  6. +3
    -1
      src/targets/satellite/target.py
  7. +26
    -6
      src/targets/waveform/target.py
  8. +4
    -2
      src/targets/waveform_phase_ratio/target.py
  9. +1
    -0
      src/targets/waveform_piggyback/__init__.py
  10. +139
    -0
      src/targets/waveform_piggyback/target.py

+ 1
- 0
setup.py View File

@ -15,6 +15,7 @@ setup(
'grond.targets',
'grond.targets.waveform',
'grond.targets.waveform_phase_ratio',
'grond.targets.waveform_piggyback',
'grond.targets.satellite',
'grond.targets.gnss_campaign',
'grond.problems',


+ 21
- 4
src/problems/base.py View File

@ -14,7 +14,7 @@ import os
import time
from pyrocko import gf, util, guts
from pyrocko.guts import Object, String, Bool, List, Dict, Int
from pyrocko.guts import Object, String, List, Dict, Int
from ..meta import ADict, Parameter, GrondError, xjoin, Forbidden
from ..targets import MisfitResult, MisfitTarget, TargetGroup, \
@ -403,12 +403,29 @@ class Problem(Object):
modelling_targets = []
t2m_map = {}
for itarget, target in enumerate(targets):
t2m_map[target] = target.prepare_modelling(engine, source)
t2m_map[target] = target.prepare_modelling(engine, source, targets)
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])
u2m_map = {}
for imtarget, mtarget in enumerate(modelling_targets):
if mtarget not in u2m_map:
u2m_map[mtarget] = []
u2m_map[mtarget].append(imtarget)
modelling_targets_unique = list(u2m_map.keys())
resp = engine.process(source, modelling_targets_unique)
modelling_results_unique = list(resp.results_list[0])
modelling_results = [None] * len(modelling_targets)
for mtarget, mresult in zip(
modelling_targets_unique, modelling_results_unique):
for itarget in u2m_map[mtarget]:
modelling_results[itarget] = mresult
imt = 0
results = []


+ 1
- 0
src/targets/__init__.py View File

@ -1,5 +1,6 @@
from .base import * # noqa
from .waveform import * # noqa
from .waveform_phase_ratio import * # noqa
from .waveform_piggyback import * # noqa
from .satellite import * # noqa
from .gnss_campaign import * # noqa

+ 1
- 4
src/targets/base.py View File

@ -26,9 +26,6 @@ class TargetGroup(Object):
default=1.0,
help='Additional manual weight of the target group')
interpolation = gf.InterpolationMethod.T()
store_id = gf.StringID.T(optional=True)
def get_targets(self, ds, event, default_path):
raise NotImplementedError()
@ -116,7 +113,7 @@ class MisfitTarget(Object):
def get_combined_weight(self):
return num.ones(1, dtype=num.float)
def prepare_modelling(self, engine, source):
def prepare_modelling(self, engine, source, targets):
return []
def finalize_modelling(


+ 3
- 1
src/targets/gnss_campaign/target.py View File

@ -33,6 +33,8 @@ class GNSSCampaignTargetGroup(TargetGroup):
optional=True,
help='List of individual campaign names (`name` in `gnss.yaml` files).')
misfit_config = GNSSCampaignMisfitConfig.T()
interpolation = gf.InterpolationMethod.T()
store_id = gf.StringID.T(optional=True)
def get_targets(self, ds, event, default_path):
logger.debug('Selecting GNSS targets...')
@ -165,7 +167,7 @@ class GNSSCampaignMisfitTarget(gf.GNSSCampaignTarget, MisfitTarget):
def get_combined_weight(self):
return num.array([self.manual_weight])
def prepare_modelling(self, engine, source):
def prepare_modelling(self, engine, source, targets):
return [self]
def finalize_modelling(


+ 3
- 1
src/targets/satellite/target.py View File

@ -48,6 +48,8 @@ class SatelliteTargetGroup(TargetGroup):
misfit_config = SatelliteMisfitConfig.T(
help='Carries the settings of the objective function for these targets'
)
interpolation = gf.InterpolationMethod.T()
store_id = gf.StringID.T(optional=True)
def get_targets(self, ds, event, default_path):
logger.debug('Selecting satellite targets...')
@ -190,7 +192,7 @@ class SatelliteMisfitTarget(gf.SatelliteTarget, MisfitTarget):
def get_combined_weight(self):
return num.array([self.manual_weight], dtype=num.float)
def prepare_modelling(self, engine, source):
def prepare_modelling(self, engine, source, targets):
return [self]
def finalize_modelling(


+ 26
- 6
src/targets/waveform/target.py View File

@ -11,6 +11,7 @@ from grond.dataset import NotFound
from grond.analysers.base import AnalyserResult
from ..base import (MisfitConfig, MisfitTarget, MisfitResult, TargetGroup)
from ..waveform_piggyback.target import WaveformPiggybackSubresult
guts_prefix = 'grond'
logger = logging.getLogger('grond.targets.waveform.target')
@ -31,8 +32,8 @@ class Trace(Object):
class WaveformMisfitConfig(MisfitConfig):
fmin = Float.T(help='minimum frequency of bandpass filter' )
fmax = Float.T(help='maximum frequency of bandpass filter' )
fmin = Float.T(help='minimum frequency of bandpass filter')
fmax = Float.T(help='maximum frequency of bandpass filter')
ffactor = Float.T(default=1.5)
tmin = gf.Timing.T(
optional=True,
@ -109,6 +110,8 @@ class WaveformTargetGroup(TargetGroup):
optional=True,
help='set channels to include, e.g. \[\'Z\',\'T\'\]')
misfit_config = WaveformMisfitConfig.T()
interpolation = gf.InterpolationMethod.T()
store_id = gf.StringID.T(optional=True)
def get_targets(self, ds, event, default_path):
logger.debug('Selecting waveform targets...')
@ -221,6 +224,8 @@ class WaveformMisfitResult(gf.Result, MisfitResult):
tshift = Float.T(optional=True)
cc = Trace.T(optional=True)
piggyback_subresults = List.T(WaveformPiggybackSubresult.T())
class WaveformMisfitTarget(gf.Target, MisfitTarget):
flip_norm = Bool.T(default=False)
@ -229,9 +234,10 @@ class WaveformMisfitTarget(gf.Target, MisfitTarget):
def __init__(self, **kwargs):
gf.Target.__init__(self, **kwargs)
MisfitTarget.__init__(self, **kwargs)
self._piggyback_subtargets = []
def string_id(self):
return '.'.join(x for x in (self.path,) + self.codes if x)
return '.'.join(x for x in (self.path,) + self.codes)
@classmethod
def get_plot_classes(cls):
@ -379,7 +385,10 @@ class WaveformMisfitTarget(gf.Target, MisfitTarget):
flip=self.flip_norm,
result_mode=self._result_mode,
tautoshift_max=config.tautoshift_max,
autoshift_penalty_max=config.autoshift_penalty_max)
autoshift_penalty_max=config.autoshift_penalty_max,
subtargets=self._piggyback_subtargets)
self._piggyback_subtargets = []
mr.tobs_shift = float(tobs_shift)
mr.tsyn_pick = float_or_none(tsyn)
@ -390,7 +399,7 @@ class WaveformMisfitTarget(gf.Target, MisfitTarget):
logger.debug(str(e))
raise gf.SeismosizerError('no waveform data, %s' % str(e))
def prepare_modelling(self, engine, source):
def prepare_modelling(self, engine, source, targets):
return [self]
def finalize_modelling(
@ -403,10 +412,13 @@ class WaveformMisfitTarget(gf.Target, MisfitTarget):
(k, getattr(self, k)) for k in gf.Target.T.propnames)
return [gf.Target(**d)]
def add_piggyback_subtarget(self, subtarget):
self._piggyback_subtargets.append(subtarget)
def misfit(
tr_obs, tr_syn, taper, domain, exponent, tautoshift_max,
autoshift_penalty_max, flip, result_mode='sparse'):
autoshift_penalty_max, flip, result_mode='sparse', subtargets=[]):
'''
Calculate misfit between observed and synthetic trace.
@ -438,6 +450,12 @@ tautoshift**2 / tautoshift_max**2``
tr_proc_obs, trspec_proc_obs = _process(tr_obs, tmin, tmax, taper, domain)
tr_proc_syn, trspec_proc_syn = _process(tr_syn, tmin, tmax, taper, domain)
piggyback_results = []
for subtarget in subtargets:
piggyback_results.append(
subtarget.evaluate(
tr_proc_obs, trspec_proc_obs, tr_proc_syn, trspec_proc_syn))
tshift = None
ctr = None
deltat = tr_proc_obs.deltat
@ -529,6 +547,8 @@ tautoshift**2 / tautoshift_max**2``
else:
assert False
result.piggyback_subresults = piggyback_results
return result


+ 4
- 2
src/targets/waveform_phase_ratio/target.py View File

@ -28,6 +28,8 @@ class PhaseRatioTargetGroup(TargetGroup):
depth_max = Float.T(optional=True)
measure_a = fm.FeatureMeasure.T()
measure_b = fm.FeatureMeasure.T()
interpolation = gf.InterpolationMethod.T()
store_id = gf.StringID.T(optional=True)
def get_targets(self, ds, event, default_path):
logger.debug('Selecting phase ratio targets...')
@ -125,9 +127,9 @@ class PhaseRatioTarget(gf.Location, MisfitTarget):
return '.'.join(x for x in (self.path,) + self.codes if x)
def get_plain_targets(self, engine, source):
return self.prepare_modelling(engine, source)
return self.prepare_modelling(engine, source, None)
def prepare_modelling(self, engine, source):
def prepare_modelling(self, engine, source, targets):
modelling_targets = []
for measure in [self.measure_a, self.measure_b]:
modelling_targets.extend(measure.get_modelling_targets(


+ 1
- 0
src/targets/waveform_piggyback/__init__.py View File

@ -0,0 +1 @@
from .target import * # noqa

+ 139
- 0
src/targets/waveform_piggyback/target.py View File

@ -0,0 +1,139 @@
import logging
import numpy as num
from pyrocko.guts import Float, Object, Int
from pyrocko import gf
from ..base import (
MisfitTarget, TargetGroup, MisfitResult)
guts_prefix = 'grond'
logger = logging.getLogger('grond.targets.waveform_piggyback.target')
def log_exclude(target, reason):
logger.debug('excluding potential target %s: %s' % (
target.string_id(), reason))
class WaveformPiggybackTargetGroup(TargetGroup):
associated_path = gf.StringID.T()
def get_targets(self, ds, event, default_path):
logger.debug('Selecting waveform piggyback targets...')
target = WaveformPiggybackTarget(
path=self.path,
associated_path=self.associated_path)
return [target]
class WaveformPiggybackSubresult(Object):
piggy_id = Int.T()
amplitude_obs = Float.T()
amplitude_syn = Float.T()
class WaveformPiggybackSubtarget(Object):
piggy_id = Int.T()
def evaluate(
self, tr_proc_obs, trspec_proc_obs, tr_proc_syn, trspec_proc_syn):
a_obs = num.sqrt(num.mean(num.abs(tr_proc_obs.ydata**2)))
a_syn = num.sqrt(num.mean(num.abs(tr_proc_syn.ydata**2)))
res = WaveformPiggybackSubresult(
piggy_id=self.piggy_id,
amplitude_obs=float(a_obs),
amplitude_syn=float(a_syn))
return res
class WaveformPiggybackMisfitResult(MisfitResult):
pass
class WaveformPiggybackTarget(MisfitTarget):
associated_path = gf.StringID.T()
_next_piggy_id = 0
def __init__(self, **kwargs):
MisfitTarget.__init__(self, **kwargs)
self.piggy_ids = set()
@classmethod
def new_piggy_id(cls):
piggy_id = cls._next_piggy_id
cls._next_piggy_id += 1
return piggy_id
def string_id(self):
return self.path
def get_plain_targets(self, engine, source):
return []
def distance_to(self, other):
return num.array([], dtype=float)
def prepare_modelling(self, engine, source, targets):
from ..waveform.target import WaveformMisfitTarget
relevant = []
for target in targets:
if isinstance(target, WaveformMisfitTarget) \
and target.path == self.associated_path:
piggy_id = self.new_piggy_id()
self.piggy_ids.add(piggy_id)
target.add_piggyback_subtarget(
WaveformPiggybackSubtarget(
piggy_id=piggy_id))
relevant.append(target)
return relevant
def finalize_modelling(
self, engine, source, modelling_targets, modelling_results):
from ..waveform.target import WaveformMisfitResult
amps_obs = []
amps_syn = []
for mtarget, mresult in zip(modelling_targets, modelling_results):
if isinstance(mresult, WaveformMisfitResult):
for sr in mresult.piggyback_subresults:
try:
self.piggy_ids.remove(sr.piggy_id)
amps_obs.append(sr.amplitude_obs)
amps_syn.append(sr.amplitude_syn)
except KeyError:
pass
amps_obs = num.array(amps_obs)
amps_syn = num.array(amps_syn)
amp_obs = num.median(amps_obs[num.isfinite(amps_obs)])
amp_syn = num.median(amps_syn[num.isfinite(amps_syn)])
m = num.log(amp_obs / amp_syn)
result = WaveformPiggybackMisfitResult(
misfits=num.array([[m, 1.]], dtype=num.float))
return result
__all__ = '''
WaveformPiggybackTargetGroup
WaveformPiggybackTarget
WaveformPiggybackMisfitResult
WaveformPiggybackSubtarget
WaveformPiggybackSubresult
'''.split()

Loading…
Cancel
Save