Coverage for /usr/local/lib/python3.11/dist-packages/grond/analysers/target_balancing/analyser.py: 90%

1import copy 

2import time 

3import logging 

4import numpy as num 

5from pyrocko.guts import Int, Float, Bool 

6from pyrocko import gf 

7 

8from grond.meta import Forbidden, GrondError 

9 

10from ..base import Analyser, AnalyserConfig, AnalyserResult 

11 

12logger = logging.getLogger('grond.analysers.target_balancer') 

13 

14 

15guts_prefix = 'grond' 

16 

17 

18class TargetBalancingAnalyser(Analyser): 

19 """ Estimating target weights that balance the signal amplitudes. 

20 

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. 

26 

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 """ 

32 

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 

38 

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: 

44 

45 logger.info( 

46 'Target balancing for "%s" at %i/%i.' % ( 

47 problem.name, 

48 iiter, niter)) 

49 

50 self._tlog_last = t 

51 

52 def analyse(self, problem, ds): 

53 if self.niter == 0: 

54 return 

55 

56 wtargets = [] 

57 if not problem.has_waveforms: 

58 return 

59 

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) 

65 

66 wproblem = problem.copy() 

67 if hasattr(problem, 'cmt_problem'): 

68 if problem.cmt_problem is not None: 

69 wproblem = problem.cmt_problem.copy() 

70 

71 wproblem.targets = wtargets 

72 

73 xbounds = wproblem.get_parameter_bounds() 

74 

75 misfits = num.zeros((self.niter, wproblem.ntargets, 2)) 

76 rstate = self.get_rstate(problem) 

77 

78 isbad_mask = None 

79 

80 self._tlog_last = 0 

81 for iiter in range(self.niter): 

82 self.log_progress(problem, iiter, self.niter) 

83 while True: 

84 if self.use_reference_magnitude: 

85 try: 

86 fixed_magnitude = wproblem.base_source.get_magnitude() 

87 except gf.DerivedMagnitudeError: 

88 raise GrondError( 

89 'Cannot use use_reference_magnitude for this type ' 

90 'of source model.') 

91 else: 

92 fixed_magnitude = None 

93 

94 x = wproblem.random_uniform( 

95 xbounds, rstate, fixed_magnitude=fixed_magnitude) 

96 

97 try: 

98 x = wproblem.preconstrain(x) 

99 break 

100 

101 except Forbidden: 

102 pass 

103 

104 if isbad_mask is not None and num.any(isbad_mask): 

105 isok_mask = num.logical_not(isbad_mask) 

106 else: 

107 isok_mask = None 

108 misfits[iiter, :, :] = wproblem.misfits(x, mask=isok_mask) 

109 

110 isbad_mask = num.isnan(misfits[iiter, :, 1]) 

111 

112 mean_ms = num.mean(misfits[:, :, 0], axis=0) 

113 mean_ps = num.mean(misfits[:, :, 1], axis=0) 

114 

115 weights = 1. / mean_ps 

116 families, nfamilies = wproblem.get_family_mask() 

117 

118 for ifamily in range(nfamilies): 

119 weights[families == ifamily] /= ( 

120 num.nansum(weights[families == ifamily]) / 

121 num.nansum(num.isfinite(weights[families == ifamily]))) 

122 

123 if self.cutoff is not None: 

124 weights[mean_ms / mean_ps > self.cutoff] = 0.0 

125 

126 for weight, target in zip(weights, problem.waveform_targets): 

127 target.analyser_results['target_balancing'] = \ 

128 TargetBalancingAnalyserResult(weight=float(weight)) 

129 

130 for itarget, target in enumerate(problem.waveform_targets): 

131 logger.info(( 

132 'Balancing analysis for target "%s":\n' 

133 ' m/p: %g\n' 

134 ' weight: %g\n' 

135 ) % ( 

136 target.string_id(), 

137 mean_ms[itarget] / mean_ps[itarget], 

138 weights[itarget])) 

139 

140 

141class TargetBalancingAnalyserResult(AnalyserResult): 

142 weight = Float.T() 

143 

144 

145class TargetBalancingAnalyserConfig(AnalyserConfig): 

146 """Configuration parameters of the target balancing.""" 

147 niterations = Int.T( 

148 default=1000, 

149 help='Number of random forward models for mean phase amplitude ' 

150 'estimation') 

151 

152 use_reference_magnitude = Bool.T( 

153 default=False, 

154 help='Fix magnitude of random sources to the magnitude of the ' 

155 'reference event.') 

156 

157 cutoff = Float.T( 

158 optional=True, 

159 help='Remove targets where ratio m/p > cutoff, where m is the misfit ' 

160 'between synthetics and observations and p is the misfit between ' 

161 'synthetics and zero-traces. Magnitude should be fixed to use ' 

162 'this.') 

163 

164 def get_analyser(self): 

165 return TargetBalancingAnalyser( 

166 niter=self.niterations, 

167 use_reference_magnitude=self.use_reference_magnitude, 

168 cutoff=self.cutoff) 

169 

170 

171__all__ = ''' 

172 TargetBalancingAnalyser 

173 TargetBalancingAnalyserConfig 

174'''.split()