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

Merge branch 'master' of gitext.gfz-potsdam.de:heimann/grond

parents 0d8e2b83 46276a26
......@@ -298,7 +298,17 @@ def command_harvest(args):
help='take NEACH best samples from each chain (default: 10)')
parser.add_option(
'--weed', dest='weed', type=int, default=0,
help='weed out bootstrap samples with bad global performance')
help='weed out bootstrap samples with bad global performance. '
'0: no weeding (default), '
'1: only bootstrap chains where all NEACH best samples '
'global misfit is less than the global average misfit of all '
'NEACH best in all chains plus one standard deviation are '
'included in the harvest ensemble, '
'2: same as 1 but additionally individual samples are '
'removed if their global misfit is greater than the global '
'average misfit of all NEACH best in all chains, '
'3: harvesting is done on the global chain only, bootstrap '
'chains are excluded')
parser, options, args = cl_parse('harvest', args, setup)
if len(args) != 1:
......
......@@ -265,12 +265,37 @@ class InnerMisfitConfig(Object):
fmin = Float.T()
fmax = Float.T()
ffactor = Float.T(default=1.5)
tmin = gf.Timing.T()
tmax = gf.Timing.T()
tfade = Float.T(optional=True)
pick_synthetic_traveltime = gf.Timing.T(optional=True)
pick_phasename = String.T(optional=True)
domain = DomainChoice.T(default='time_domain')
tmin = gf.Timing.T(
help='Start of main time window used for waveform fitting.')
tmax = gf.Timing.T(
help='End of main time window used for waveform fitting.')
tfade = Float.T(
optional=True,
help='Decay time of taper prepended and appended to main time window '
'used for waveform fitting [s].')
pick_synthetic_traveltime = gf.Timing.T(
optional=True,
help='Synthetic phase arrival definition for alignment of observed '
'and synthetic traces.')
pick_phasename = String.T(
optional=True,
help='Name of picked phase for alignment of observed and synthetic '
'traces.')
domain = DomainChoice.T(
default='time_domain',
help='Type of data characteristic to be fitted.\n\nAvailable choices '
'are: %s' % ', '.join("``'%s'``" % s
for s in DomainChoice.choices))
tautoshift_max = Float.T(
default=0.0,
help='If non-zero, allow synthetic and observed traces to be shifted '
'against each other by up to +/- the given value [s].')
autoshift_penalty_max = Float.T(
default=0.0,
help='If non-zero, a penalty misfit is added for non-zero shift '
'values.\n\nThe penalty value is computed as '
'``autoshift_penalty_max * normalization_factor * tautoshift**2 '
'/ tautoshift_max**2``')
def get_full_frequency_range(self):
return self.fmin / self.ffactor, self.fmax * self.ffactor
......@@ -466,7 +491,9 @@ class MisfitTarget(gf.Target):
domain=config.domain,
exponent=2,
flip=self.flip_norm,
result_mode=self._result_mode)
result_mode=self._result_mode,
tautoshift_max=config.tautoshift_max,
autoshift_penalty_max=config.autoshift_penalty_max)
mr.tobs_shift = float(tobs_shift)
mr.tsyn_pick = float_or_none(tsyn)
......@@ -479,7 +506,8 @@ class MisfitTarget(gf.Target):
def misfit(
tr_obs, tr_syn, taper, domain, exponent, flip, result_mode='sparse'):
tr_obs, tr_syn, taper, domain, exponent, tautoshift_max,
autoshift_penalty_max, flip, result_mode='sparse'):
'''
Calculate misfit between observed and synthetic trace.
......@@ -490,6 +518,12 @@ def misfit(
:py:class:`pyrocko.trace.Taper`
:param domain: how to calculate difference, see :py:class:`DomainChoice`
:param exponent: exponent of Lx type norms
:param tautoshift_max: if non-zero, return lowest misfit when traces are
allowed to shift against each other by up to +/- ``tautoshift_max``
:param autoshift_penalty_max: if non-zero, a penalty misfit is added for
for non-zero shift values. The penalty value is
``autoshift_penalty_max * normalization_factor * \
tautoshift**2 / tautoshift_max**2``
:param flip: ``bool``, if set to ``True``, normalization factor is
computed against *tr_syn* rather than *tr_obs*
:param result_mode: ``'full'``, include traces and spectra or ``'sparse'``,
......@@ -506,12 +540,39 @@ def misfit(
cc_shift = None
ctr = None
deltat = tr_proc_obs.deltat
if domain in ('time_domain', 'envelope', 'absolute'):
a, b = tr_proc_syn.ydata, tr_proc_obs.ydata
if flip:
b, a = a, b
m, n = trace.Lx_norm(a, b, norm=exponent)
nshift_max = max(0, min(a.size-1,
int(math.floor(tautoshift_max / deltat))))
if nshift_max == 0:
m, n = trace.Lx_norm(a, b, norm=exponent)
else:
mns = []
for ishift in xrange(-nshift_max, nshift_max+1):
if ishift < 0:
a_cut = a[-ishift:]
b_cut = b[:ishift]
elif ishift == 0:
a_cut = a
b_cut = b
elif ishift > 0:
a_cut = a[:-ishift]
b_cut = b[ishift:]
mns.append(trace.Lx_norm(a_cut, b_cut, norm=exponent))
ms, ns = num.array(mns).T
iarg = num.argmin(ms)
tshift = (iarg-nshift_max)*deltat
m, n = ms[iarg], ns[iarg]
m += autoshift_penalty_max * n * tshift**2 / tautoshift_max**2
elif domain == 'cc_max_norm':
......
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