base.py 11.4 KB
Newer Older
Marius Isken's avatar
Marius Isken committed
1
2
3
4
5
6
7
8
9
import numpy as num
import copy
import logging
import os.path as op

from pyrocko import gf, util, guts
from pyrocko.guts import Object, String, Bool, List, Dict, Int

from ..meta import ADict, Parameter, GrondError
Marius Isken's avatar
Marius Isken committed
10
from ..targets import WaveformMisfitTarget, SatelliteMisfitTarget
Marius Isken's avatar
Marius Isken committed
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139


guts_prefix = 'grond'
logger = logging.getLogger('grond')
km = 1e3
as_km = dict(scale_factor=km, scale_unit='km')


class ProblemConfig(Object):
    name_template = String.T()
    apply_balancing_weights = Bool.T(default=True)
    norm_exponent = Int.T(default=2)


class Problem(Object):
    name = String.T()
    ranges = Dict.T(String.T(), gf.Range.T())
    dependants = List.T(Parameter.T())
    apply_balancing_weights = Bool.T(default=True)
    norm_exponent = Int.T(default=2)
    base_source = gf.Source.T(optional=True)

    targets = List.T()

    def __init__(self, **kwargs):
        Object.__init__(self, **kwargs)

        self._bootstrap_weights = None
        self._target_weights = None
        self._engine = None
        self._group_mask = None

        if hasattr(self, 'problem_waveform_parameters') and self.has_waveforms:
            self.problem_parameters =\
                self.problem_parameters + self.problem_waveform_parameters

        logger.name = self.__class__.__name__

    def get_engine(self):
        return self._engine

    def copy(self):
        o = copy.copy(self)
        o._bootstrap_weights = None
        o._target_weights = None
        return o

    def set_target_parameter_values(self, x):
        i = len(self.problem_parameters)
        for target in self.targets:
            n = len(target.target_parameters)
            target.set_parameter_values(x[i:i+n])
            i += n

    def get_parameter_dict(self, model, group=None):
        params = []
        for ip, p in enumerate(self.parameters):
            if group in p.groups:
                params.append((p.name, model[ip]))
        return ADict(params)

    def get_parameter_array(self, d):
        arr = num.zeros(self.nparameters, dtype=num.float)
        for ip, p in enumerate(self.parameters):
            if p.name in d.keys():
                arr[ip] = d[p.name]
        return arr

    def dump_problem_info(self, dirname):
        fn = op.join(dirname, 'problem.yaml')
        util.ensuredirs(fn)
        guts.dump(self, filename=fn)

    def dump_problem_data(
            self, dirname, x, ms, ns,
            accept=None, ibootstrap_choice=None, ibase=None):

        fn = op.join(dirname, 'models')
        if not isinstance(x, num.ndarray):
            x = num.array(x)
        with open(fn, 'ab') as f:
            x.astype('<f8').tofile(f)

        fn = op.join(dirname, 'misfits')
        with open(fn, 'ab') as f:
            ms.astype('<f8').tofile(f)
            ns.astype('<f8').tofile(f)

        if None not in (ibootstrap_choice, ibase):
            fn = op.join(dirname, 'choices')
            with open(fn, 'ab') as f:
                num.array((ibootstrap_choice, ibase), dtype='<i8').tofile(f)

        if accept is not None:
            fn = op.join(dirname, 'accepted')
            with open(fn, 'ab') as f:
                accept.astype('<i1').tofile(f)

    def name_to_index(self, name):
        pnames = [p.name for p in self.combined]
        return pnames.index(name)

    @property
    def parameters(self):
        target_parameters = []
        for target in self.targets:
            target_parameters.extend(target.target_parameters)
        return self.problem_parameters + target_parameters

    @property
    def parameter_names(self):
        return [p.name for p in self.combined]

    @property
    def nparameters(self):
        return len(self.parameters)

    @property
    def ntargets(self):
        return len(self.targets)

    @property
    def ntargets_waveform(self):
        return len(self.waveform_targets)

    @property
    def ntargets_static(self):
        return len(self.satellite_targets)

Marius Isken's avatar
Marius Isken committed
140
141
142
143
144
145
146
    @property
    def nmisfits(self):
        nmisfits = 0
        for target in self.targets:
            nmisfits += target.nmisfits
        return nmisfits

Marius Isken's avatar
Marius Isken committed
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
    @property
    def ndependants(self):
        return len(self.dependants)

    @property
    def ncombined(self):
        return len(self.parameters) + len(self.dependants)

    @property
    def combined(self):
        return self.parameters + self.dependants

    @property
    def satellite_targets(self):
        return [t for t in self.targets
Marius Isken's avatar
Marius Isken committed
162
                if isinstance(t, SatelliteMisfitTarget)]
Marius Isken's avatar
Marius Isken committed
163
164
165
166

    @property
    def waveform_targets(self):
        return [t for t in self.targets
Marius Isken's avatar
Marius Isken committed
167
                if isinstance(t, WaveformMisfitTarget)]
Marius Isken's avatar
Marius Isken committed
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207

    @property
    def has_statics(self):
        if self.satellite_targets:
            return True
        return False

    @property
    def has_waveforms(self):
        if self.waveform_targets:
            return True
        return False

    def set_engine(self, engine):
        self._engine = engine

    def random_uniform(self, xbounds):
        raise NotImplementedError()

    def make_bootstrap_weights(self, nbootstrap, type='classic'):
        ntargets = self.ntargets
        ws = num.zeros((nbootstrap, ntargets))
        rstate = num.random.RandomState(23)
        for ibootstrap in xrange(nbootstrap):
            if type == 'classic':
                ii = rstate.randint(0, ntargets, size=self.ntargets)
                ws[ibootstrap, :] = num.histogram(
                    ii, ntargets, (-0.5, ntargets - 0.5))[0]
            elif type == 'bayesian':
                f = rstate.uniform(0., 1., size=self.ntargets+1)
                f[0] = 0.
                f[-1] = 1.
                f = num.sort(f)
                g = f[1:] - f[:-1]
                ws[ibootstrap, :] = g * ntargets
            else:
                assert False

        return ws

Sebastian Heimann's avatar
wip    
Sebastian Heimann committed
208
    def get_bootstrap_weights(self, nbootstrap, ibootstrap=None):
Marius Isken's avatar
Marius Isken committed
209
210
        if self._bootstrap_weights is None:
            self._bootstrap_weights = self.make_bootstrap_weights(
Sebastian Heimann's avatar
wip    
Sebastian Heimann committed
211
                nbootstrap, type='bayesian')
Marius Isken's avatar
Marius Isken committed
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263

        if ibootstrap is None:
            return self._bootstrap_weights
        else:
            return self._bootstrap_weights[ibootstrap, :]

    def get_target_weights(self):
        if self._target_weights is None:
            self._target_weights = num.array(
                [target.get_combined_weight(
                    apply_balancing_weights=self.apply_balancing_weights)
                 for target in self.targets], dtype=num.float)

        return self._target_weights

    def inter_group_weights(self, ns):
        exp, root = self.get_norm_functions()

        group, ngroups = self.get_group_mask()

        ws = num.zeros(self.ntargets)
        for igroup in xrange(ngroups):
            mask = group == igroup
            ws[mask] = 1.0 / root(num.nansum(exp(ns[mask])))

        return ws

    def inter_group_weights2(self, ns):
        exp, root = self.get_norm_functions()

        group, ngroups = self.get_group_mask()
        ws = num.zeros(ns.shape)
        for igroup in xrange(ngroups):
            mask = group == igroup
            ws[:, mask] = (1.0 / root(
                num.nansum(exp(ns[:, mask]), axis=1)))[:, num.newaxis]

        return ws

    def get_xref(self):
        return self.pack(self.base_source)

    def get_parameter_bounds(self):
        out = []
        for p in self.problem_parameters:
            r = self.ranges[p.name]
            out.append((r.start, r.stop))

        for target in self.targets:
            for p in target.target_parameters:
                r = target.target_ranges[p.name]
                out.append((r.start, r.stop))
Sebastian Heimann's avatar
wip    
Sebastian Heimann committed
264
265

        return num.array(out, dtype=num.float)
Marius Isken's avatar
Marius Isken committed
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347

    def get_dependant_bounds(self):
        return []

    def raise_invalid_norm_exponent(self):
        raise GrondError('invalid norm exponent' % self.norm_exponent)

    def get_norm_functions(self):
        if self.norm_exponent == 2:
            def sqr(x):
                return x**2

            return sqr, num.sqrt

        elif self.norm_exponent == 1:
            def noop(x):
                return x

            return noop, num.abs

        else:
            self.raise_invalid_norm_exponent()

    def bootstrap_misfit(self, ms, ns, ibootstrap=None):
        exp, root = self.get_norm_functions()

        w = self.get_bootstrap_weights(ibootstrap) * \
            self.get_target_weights() * self.inter_group_weights(ns)

        if ibootstrap is None:
            return root(
                num.nansum(exp(w*ms[num.newaxis, :]), axis=1) /
                num.nansum(exp(w*ns[num.newaxis, :]), axis=1))
        else:
            return root(num.nansum(exp(w*ms)) / num.nansum(exp(w*ns)))

    def bootstrap_misfits(self, misfits, ibootstrap):
        exp, root = self.get_norm_functions()

        w = self.get_bootstrap_weights(ibootstrap)[num.newaxis, :] * \
            self.get_target_weights()[num.newaxis, :] * \
            self.inter_group_weights2(misfits[:, :, 1])

        bms = root(num.nansum(exp(w*misfits[:, :, 0]), axis=1) /
                   num.nansum(exp(w*misfits[:, :, 1]), axis=1))

        # From Henriette
        # w = self.get_target_weights()[num.newaxis, :] * \
        #     self.inter_group_weights2(misfits[:, :, 1])
        # #w = self.get_bootstrap_weights(ibootstrap)[num.newaxis, :] * \
        # #    self.get_target_weights()[num.newaxis, :] * \
        # #    self.inter_group_weights2(misfits[:, :, 1])

        # bms = num.sqrt(num.nansum((w*misfits[:, :, 0])**2, axis=1))
        return bms

    def global_misfit(self, ms, ns):
        exp, root = self.get_norm_functions()

        ws = self.get_target_weights() * self.inter_group_weights(ns)
        m = root(num.nansum(exp(ws*ms)) / num.nansum(exp(ws*ns)))
        return m

    def global_misfits(self, misfits):
        exp, root = self.get_norm_functions()
        ws = self.get_target_weights()[num.newaxis, :] * \
            self.inter_group_weights2(misfits[:, :, 1])
        gms = root(num.nansum(exp(ws*misfits[:, :, 0]), axis=1) /
                   num.nansum(exp(ws*misfits[:, :, 1]), axis=1))
        return gms

    def global_contributions(self, misfits):
        exp, root = self.get_norm_functions()
        ws = self.get_target_weights()[num.newaxis, :] * \
            self.inter_group_weights2(misfits[:, :, 1])

        gcms = exp(ws*misfits[:, :, 0]) / \
            num.nansum(exp(ws*misfits[:, :, 1]), axis=1)[:, num.newaxis]

        return gcms

    def make_group_mask(self):
348
349
        family_names = set()
        families = num.zeros(len(self.targets), dtype=num.int)
Marius Isken's avatar
Marius Isken committed
350

351
352
353
        for itarget, target in enumerate(self.targets):
            family_names.add(target.normalisation_family)
            families[itarget] = len(family_names) - 1
Marius Isken's avatar
Marius Isken committed
354

355
        return families, len(family_names)
Marius Isken's avatar
Marius Isken committed
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382

    def get_group_mask(self):
        if self._group_mask is None:
            self._group_mask = self.make_group_mask()

        return self._group_mask

    def evaluate(self, x, dump_data=False):
        raise NotImplementedError()

    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