1import copy
2import time
3import logging
4import numpy as num
5from pyrocko.guts import Int, Float, Bool
6from pyrocko import gf
8from grond.meta import Forbidden, GrondError
10from ..base import Analyser, AnalyserConfig, AnalyserResult
12logger = logging.getLogger('grond.analysers.target_balancer')
15guts_prefix = 'grond'
18class TargetBalancingAnalyser(Analyser):
19 """ Estimating target weights that balance the signal amplitudes.
21 Signal amplitudes depend on the source-receiver distance, on the
22 phase type and the taper used. Large signals have in general
23 a higher contribution to the misfit than smaller signals,
24 without carrying more information. With this function, weights are
25 estimated that shall balance the phase contributions.
27 The weight estimation is based on synthetic waveforms that stem from
28 a given number of random forward models. The inverse of the mean
29 synthetic signal amplitudes gives the balancing weight. This is
30 described as adaptive station weighting in Heimann (2011).
31 """
33 def __init__(self, niter, use_reference_magnitude, cutoff):
34 Analyser.__init__(self)
35 self.niter = niter
36 self.use_reference_magnitude = use_reference_magnitude
37 self.cutoff = cutoff
39 def log_progress(self, problem, iiter, niter):
40 t = time.time()
41 if self._tlog_last < t - 10. \
42 or iiter == 0 \
43 or iiter == niter - 1:
45 logger.info(
46 'Target balancing for "%s" at %i/%i.' % (
47 problem.name,
48 iiter, niter))
50 self._tlog_last = t
52 def analyse(self, problem, ds):
53 if self.niter == 0:
54 return
56 wtargets = []
57 if not problem.has_waveforms:
58 return
60 for target in problem.waveform_targets:
61 wtarget = copy.copy(target)
62 wtarget.flip_norm = True
63 wtarget.weight = 1.0
64 wtargets.append(wtarget)
66 wproblem = problem.copy()
67 wproblem.targets = wtargets
69 xbounds = wproblem.get_parameter_bounds()
71 misfits = num.zeros((self.niter, wproblem.ntargets, 2))
72 rstate = num.random.RandomState(123)
74 isbad_mask = None
76 self._tlog_last = 0
77 for iiter in range(self.niter):
78 self.log_progress(problem, iiter, self.niter)
79 while True:
80 if self.use_reference_magnitude:
81 try:
82 fixed_magnitude = wproblem.base_source.get_magnitude()
83 except gf.DerivedMagnitudeError:
84 raise GrondError(
85 'Cannot use use_reference_magnitude for this type '
86 'of source model.')
87 else:
88 fixed_magnitude = None
90 x = wproblem.random_uniform(
91 xbounds, rstate, fixed_magnitude=fixed_magnitude)
93 try:
94 x = wproblem.preconstrain(x)
95 break
97 except Forbidden:
98 pass
100 if isbad_mask is not None and num.any(isbad_mask):
101 isok_mask = num.logical_not(isbad_mask)
102 else:
103 isok_mask = None
104 misfits[iiter, :, :] = wproblem.misfits(x, mask=isok_mask)
106 isbad_mask = num.isnan(misfits[iiter, :, 1])
108 mean_ms = num.mean(misfits[:, :, 0], axis=0)
110 mean_ps = num.mean(misfits[:, :, 1], axis=0)
112 weights = 1. / mean_ps
113 families, nfamilies = wproblem.get_family_mask()
115 for ifamily in range(nfamilies):
116 weights[families == ifamily] /= (
117 num.nansum(weights[families == ifamily]) /
118 num.nansum(num.isfinite(weights[families == ifamily])))
120 if self.cutoff is not None:
121 weights[mean_ms / mean_ps > self.cutoff] = 0.0
123 for weight, target in zip(weights, problem.waveform_targets):
124 target.analyser_results['target_balancing'] = \
125 TargetBalancingAnalyserResult(weight=float(weight))
127 for itarget, target in enumerate(problem.waveform_targets):
128 logger.info((
129 'Balancing analysis for target "%s":\n'
130 ' m/p: %g\n'
131 ' weight: %g\n'
132 ) % (
133 target.string_id(),
134 mean_ms[itarget] / mean_ps[itarget],
135 weights[itarget]))
138class TargetBalancingAnalyserResult(AnalyserResult):
139 weight = Float.T()
142class TargetBalancingAnalyserConfig(AnalyserConfig):
143 """Configuration parameters of the target balancing."""
144 niterations = Int.T(
145 default=1000,
146 help='Number of random forward models for mean phase amplitude '
147 'estimation')
149 use_reference_magnitude = Bool.T(
150 default=False,
151 help='Fix magnitude of random sources to the magnitude of the '
152 'reference event.')
154 cutoff = Float.T(
155 optional=True,
156 help='Remove targets where ratio m/p > cutoff, where m is the misfit '
157 'between synthetics and observations and p is the misfit between '
158 'synthetics and zero-traces. Magnitude should be fixed to use '
159 'this.')
161 def get_analyser(self):
162 return TargetBalancingAnalyser(
163 niter=self.niterations,
164 use_reference_magnitude=self.use_reference_magnitude,
165 cutoff=self.cutoff)
168__all__ = '''
169 TargetBalancingAnalyser
170 TargetBalancingAnalyserConfig
171'''.split()