initial version

.gitignore

+build
+*.pyc

apps/grond

+#!/usr/bin/env python
+
+import math
+import sys
+import logging
+from optparse import OptionParser
+
+from pyrocko import util
+from pyrocko.gf import Range
+
+import grond
+from grond import plot
+
+logger = logging.getLogger('main')
+
+km = 1000.
+
+
+def d2u(d):
+    if isinstance(d, dict):
+        return dict((k.replace('-', '_'), v) for (k, v) in d.iteritems())
+    else:
+        return d.replace('-', '_')
+
+
+subcommand_descriptions = {
+    'init': 'print example configuration',
+    'go': 'run Grond optimization',
+    'forward': 'run forward modelling',
+    'harvest': 'manually run harvesting',
+    'map-geometry': 'make station map',
+    'plot': 'plot optimization result',
+    'export': 'export results',
+}
+
+subcommand_usages = {
+    'init': 'init [options]',
+    'go': 'go <configfile> [options]',
+    'forward': 'forward <rundir> [options]',
+    'harvest': 'harvest <rundir> [options]',
+    'map-geometry': 'map-geometry <configfile> [options]',
+    'plot': 'plot <plotnames> <rundir> [options]',
+    'export': 'export <rundir> <outputfile> [options]',
+}
+
+subcommands = subcommand_descriptions.keys()
+
+program_name = 'grond'
+
+usage_tdata = d2u(subcommand_descriptions)
+usage_tdata['program_name'] = program_name
+
+usage = '''%(program_name)s <subcommand> [options] [--] <arguments> ...
+
+Subcommands:
+
+    init          %(init)s
+    go            %(go)s
+    forward       %(forward)s
+    harvest       %(harvest)s
+    map-geometry  %(map_geometry)s
+    plot          %(plot)s
+    export        %(export)s
+
+To get further help and a list of available options for any subcommand run:
+
+    %(program_name)s <subcommand> --help
+
+''' % usage_tdata
+
+
+def add_common_options(parser):
+    parser.add_option(
+        '--loglevel',
+        action='store',
+        dest='loglevel',
+        type='choice',
+        choices=('critical', 'error', 'warning', 'info', 'debug'),
+        default='info',
+        help ='set logger level to '
+              '"critical", "error", "warning", "info", or "debug". '
+              'Default is "%default".')
+
+
+def process_common_options(options):
+    util.setup_logging(program_name, options.loglevel)
+
+
+def cl_parse(command, args, setup=None):
+    usage = subcommand_usages[command]
+    descr = subcommand_descriptions[command]
+
+    if isinstance(usage, basestring):
+        usage = [usage]
+
+    susage = '%s %s' % (program_name, usage[0])
+    for s in usage[1:]:
+        susage += '\n%s%s %s' % (' '*7, program_name, s)
+
+    parser = OptionParser(
+        usage=susage,
+        description=descr[0].upper() + descr[1:] + '.')
+
+    if setup:
+        setup(parser)
+
+    add_common_options(parser)
+    (options, args) = parser.parse_args(args)
+    process_common_options(options)
+    return parser, options, args
+
+
+def die(message, err=''):
+    if err:
+        sys.exit('%s: error: %s \n %s' % (program_name, message, err))
+    else:
+        sys.exit('%s: error: %s' % (program_name, message))
+
+
+def help_and_die(parser, message):
+    parser.print_help(sys.stderr)
+    sys.stderr.write('\n')
+    die(message)
+
+
+def command_init(args):
+
+    dataset_config = grond.DatasetConfig(
+        stations_path='stations.txt',
+        events_path='events.txt',
+        waveform_paths=['data'])
+
+    target_configs = [grond.TargetConfig(
+        distance_min=10*km,
+        distance_max=1000*km,
+        channels=['Z', 'R', 'T'],
+        interpolation='multilinear',
+        store_id='global_2s',
+        inner_misfit_config=grond.InnerMisfitConfig(
+            fmin=0.01,
+            fmax=0.1))]
+
+    s2 = math.sqrt(2.0)
+
+    problem_config = grond.CMTProblemConfig(
+        name_template='cmt_%(event_name)s',
+        distance_min=2.*km,
+        nbootstrap=100,
+        mt_type='deviatoric',
+        ranges=dict(
+            time=Range(0, 10.0, relative='add'),
+            north_shift=Range(-16*km, 16*km),
+            east_shift=Range(-16*km, 16*km),
+            depth=Range(1*km, 11*km),
+            magnitude=Range(4.0, 6.0),
+            rmnn=Range(-s2, s2),
+            rmee=Range(-s2, s2),
+            rmdd=Range(-s2, s2),
+            rmne=Range(-1.0, 1.0),
+            rmnd=Range(-1.0, 1.0),
+            rmed=Range(-1.0, 1.0),
+            duration=Range(1.0, 15.0)))
+
+    config = grond.Config(
+        dataset_config=dataset_config,
+        target_configs=target_configs,
+        problem_config=problem_config)
+
+    print config
+
+
+def command_go(args):
+    def setup(parser):
+        parser.add_option(
+            '--force', dest='force', action='store_true',
+            help='overwrite existing run directory')
+        parser.add_option(
+            '--status', dest='status', default='state',
+            help='status output selection (choices: state, matrix)')
+
+    parser, options, args = cl_parse('go', args, setup)
+    if len(args) != 1:
+        help_and_die(parser, 'no config file given')
+
+    config_path = args[0]
+    config = grond.read_config(config_path)
+    if options.status == 'quiet':
+        status = ()
+    else:
+        status = tuple(options.status.split(','))
+
+    grond.go(config, force=options.force, status=status)
+
+
+def command_forward(args):
+    parser, options, args = cl_parse('forward', args)
+    if len(args) != 1:
+        help_and_die(parser, 'incorrect number of arguments')
+
+    run_path, = args
+    grond.forward(run_path)
+
+
+def command_harvest(args):
+    def setup(parser):
+        parser.add_option(
+            '--force', dest='force', action='store_true',
+            help='overwrite existing harvest directory')
+        parser.add_option(
+            '--neach', dest='neach', type=int, default=10,
+            help='take NEACH best samples from each chain (default: 10)')
+
+    parser, options, args = cl_parse('harvest', args, setup)
+    if len(args) != 1:
+        help_and_die(parser, 'no rundir')
+
+    run_path, = args
+    grond.harvest(run_path, force=options.force, nbest=options.neach)
+
+
+def command_map_geometry(args):
+    parser, options, args = cl_parse('map-geometry', args)
+    if len(args) != 2:
+        help_and_die(parser, 'two arguments required')
+
+    config_path = args[0]
+    output_path = args[1]
+    config = grond.read_config(config_path)
+    plot.map_geometry(config, output_path)
+
+
+def command_plot(args):
+    parser, options, args = cl_parse('plot', args)
+    if len(args) != 2:
+        help_and_die(parser, 'two arguments required')
+
+    plotnames = args[0].split(',')
+    dirname = args[1]
+    plot.plot_result(dirname, plotnames)
+
+
+def command_export(args):
+    parser, options, args = cl_parse('export', args)
+    if len(args) not in (1, 2):
+        help_and_die(parser, 'one or two arguments required')
+
+    dirname = args[0]
+    if len(args) == 2:
+        filename = args[1]
+    else:
+        filename = None
+
+    grond.export(dirname, filename)
+
+
+if __name__ == '__main__':
+
+    if len(sys.argv) < 2:
+        sys.exit('Usage: %s' % usage)
+
+    args = list(sys.argv)
+    args.pop(0)
+    command = args.pop(0)
+
+    if command in subcommands:
+        globals()['command_' + d2u(command)](args)
+
+    elif command in ('--help', '-h', 'help'):
+        if command == 'help' and args:
+            acommand = args[0]
+            if acommand in subcommands:
+                globals()['command_' + acommand](['--help'])
+
+        sys.exit('Usage: %s' % usage)
+
+    else:
+        die('no such subcommand: %s' % command)

setup.py

+#!/usr/bin/env python
+
+from distutils.core import setup
+
+setup(
+    name='grond',
+    description='What do you want to bust today?!',
+    author='Sebastian Heimann',
+    author_email='sebastian.heimann@gfz-potsdam.de',
+    packages=['grond'],
+    package_dir={'grond': 'src'},
+    scripts=['apps/grond'],
+    package_data={'grond': []})

src/__init__.py

+from core import *  # noqa
+from cmt import *  # noqa
+from dataset import *  # noqa

src/cmt.py

+import math
+import logging
+
+import numpy as num
+
+from pyrocko import gf, moment_tensor as mtm, util
+from pyrocko.guts import Float, String, Dict, List, Int, StringChoice
+from grond import core
+
+guts_prefix = 'grond'
+
+
+logger = logging.getLogger('grond.cmt')
+
+km = 1000.
+
+as_km = dict(scale_factor=km, scale_unit='km')
+
+
+class CMTProblem(core.Problem):
+
+    parameters = [
+        core.Parameter('time', 's'),
+        core.Parameter('north_shift', 'm', **as_km),
+        core.Parameter('east_shift', 'm', **as_km),
+        core.Parameter('depth', 'm', **as_km),
+        core.Parameter('magnitude'),
+        core.Parameter('rmnn'),
+        core.Parameter('rmee'),
+        core.Parameter('rmdd'),
+        core.Parameter('rmne'),
+        core.Parameter('rmnd'),
+        core.Parameter('rmed'),
+        core.Parameter('duration')]
+
+    dependants = [
+        core.Parameter('rel_moment_iso'),
+        core.Parameter('rel_moment_clvd')]
+
+    base_source = gf.Source.T()
+    targets = List.T(core.MisfitTarget.T())
+
+    ranges = Dict.T(String.T(), gf.Range.T())
+    distance_min = Float.T(default=0.0)
+    nbootstrap = Int.T(default=10)
+    mt_type = StringChoice.T(choices=['full', 'deviatoric'])
+
+    def unpack(self, x):
+        d = self.parameter_dict(x)
+        rm6 = num.array([d.rmnn, d.rmee, d.rmdd, d.rmne, d.rmnd, d.rmed],
+                        dtype=num.float)
+
+        m0 = mtm.magnitude_to_moment(d.magnitude)
+        m6 = rm6 * m0
+
+        p = {}
+        for k in self.base_source.keys():
+            if k in d:
+                p[k] = float(
+                    self.ranges[k].make_relative(self.base_source[k], d[k]))
+
+        stf = gf.HalfSinusoidSTF(duration=float(d.duration))
+
+        source = self.base_source.clone(m6=m6, stf=stf, **p)
+        return source
+
+    def make_dependant(self, xs, pname):
+        if xs.ndim == 1:
+            return self.make_dependant(xs[num.newaxis, :], pname)[0]
+
+        if pname in ('rel_moment_iso', 'rel_moment_clvd'):
+            y = num.zeros(xs.shape[0])
+            for i, x in enumerate(xs):
+                source = self.unpack(x)
+                mt = source.pyrocko_moment_tensor()
+                res = mt.standard_decomposition()
+
+                if pname == 'rel_moment_iso':
+                    ratio_iso, m_iso = res[0][1:3]
+                    y[i] = ratio_iso * num.sign(m_iso[0, 0])
+                else:
+                    ratio_clvd, m_clvd = res[2][1:3]
+                    evals, evecs = mtm.eigh_check(m_clvd)
+                    ii = num.argmax(num.abs(evals))
+                    y[i] = ratio_clvd * num.sign(evals[ii])
+
+            return y
+
+        else:
+            raise KeyError(pname)
+
+    def extract(self, xs, i):
+        if xs.ndim == 1:
+            return self.extract(xs[num.newaxis, :], i)[0]
+
+        if i < self.nparameters:
+            return xs[:, i]
+        else:
+            return self.make_dependant(
+                xs, self.dependants[i-self.nparameters].name)
+
+    def pack(self, source):
+        m6 = source.m6
+        mt = source.pyrocko_moment_tensor()
+        rm6 = m6 / mt.scalar_moment()
+
+        x = num.array([
+            source.time,
+            source.north_shift,
+            source.east_shift,
+            source.depth,
+            mt.moment_magnitude(),
+            ] + rm6.tolist() + [source.stf.duration], dtype=num.float)
+
+        return x
+
+    def preconstrain(self, x):
+
+        d = self.parameter_dict(x)
+        m6 = num.array([d.rmnn, d.rmee, d.rmdd, d.rmne, d.rmnd, d.rmed],
+                       dtype=num.float)
+
+        m9 = mtm.symmat6(*m6)
+        if self.mt_type == 'deviatoric':
+            trace_m = num.trace(m9)
+            m_iso = num.diag([trace_m / 3., trace_m / 3., trace_m / 3.])
+            m9 -= m_iso
+
+        m0_unscaled = math.sqrt(num.sum(m9.A**2)) / math.sqrt(2.)
+
+        m9 /= m0_unscaled
+        m6 = mtm.to6(m9)
+        d.rmnn, d.rmee, d.rmdd, d.rmne, d.rmnd, d.rmed = m6
+        x = self.parameter_array(d)
+
+        source = self.unpack(x)
+        if any(self.distance_min > source.distance_to(t)
+               for t in self.targets):
+            raise core.Forbidden()
+
+        return x
+
+    def bounds(self):
+        out = []
+        for p in self.parameters:
+            r = self.ranges[p.name]
+            out.append((r.start, r.stop))
+
+        return out
+
+    def dependant_bounds(self):
+        out = [
+            (-1., 1.),
+            (-1., 1.)]
+
+        return out
+
+    def evaluate(self, x, return_traces=False):
+        source = self.unpack(x)
+        engine = gf.get_engine()
+        for target in self.targets:
+            target.set_return_traces(return_traces)
+
+        resp = engine.process(source, self.targets)
+        data = []
+        results = []
+        for source, target, result in resp.iter_results('results'):
+            if isinstance(result, gf.SeismosizerError):
+                logger.debug(
+                    '%s.%s.%s.%s: %s' % (target.codes + (str(result),)))
+
+                data.append((None, None))
+                if return_traces:
+                    results.append(None)
+            else:
+                data.append((result.misfit_value, result.misfit_norm))
+                if return_traces:
+                    results.append(result)
+
+        ms, ns = num.array(data, dtype=num.float).T
+        if return_traces:
+            return ms, ns, results
+        else:
+            return ms, ns
+
+    def forward(self, x):
+        source = self.unpack(x)
+        engine = gf.get_engine()
+        plain_targets = [target.get_plain_target() for target in self.targets]
+
+        resp = engine.process(source, plain_targets)
+        results = []
+        for source, target, result in resp.iter_results('results'):
+            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
+
+    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 bootstrap_misfit(self, ms, ns, ibootstrap=None):
+        w = self.get_bootstrap_weights(ibootstrap) * self.get_target_weights()
+        if ibootstrap is None:
+            return num.sqrt(
+                num.nansum((w*ms[num.newaxis, :])**2, axis=1) /
+                num.nansum((w*ns[num.newaxis, :])**2, axis=1))
+        else:
+            return num.sqrt(num.nansum((w*ms)**2) / num.nansum((w*ns)**2))
+
+    def bootstrap_misfits(self, misfits, ibootstrap):
+        w = self.get_bootstrap_weights(ibootstrap)[num.newaxis, :] * \
+            self.get_target_weights()[num.newaxis, :]
+
+        bms = num.sqrt(num.nansum((w*misfits[:, :, 0])**2, axis=1) /
+                       num.nansum((w*misfits[:, :, 1])**2, axis=1))
+        return bms
+
+    def global_misfit(self, ms, ns):
+        ws = self.get_target_weights()
+        m = num.sqrt(num.nansum((ws*ms)**2) / num.nansum((ws*ns)**2))
+        return m
+
+    def global_misfits(self, misfits):
+        ws = self.get_target_weights()[num.newaxis, :]
+        gms = num.sqrt(num.nansum((ws*misfits[:, :, 0])**2, axis=1) /
+                       num.nansum((ws*misfits[:, :, 1])**2, axis=1))
+        return gms
+
+    def global_contributions(self, misfits):
+        ws = self.get_target_weights()[num.newaxis, :]
+        gcms = (ws*misfits[:, :, 0])**2 / \
+            num.nansum((ws*misfits[:, :, 1])**2, axis=1)[:, num.newaxis]
+
+        return gcms
+
+
+class CMTProblemConfig(core.ProblemConfig):
+
+    ranges = Dict.T(String.T(), gf.Range.T())
+    distance_min = Float.T(default=0.0)
+    nbootstrap = Int.T(default=10)
+    mt_type = StringChoice.T(choices=['full', 'deviatoric'])
+
+    def get_problem(self, event, targets):
+        if event.depth is None:
+            event.depth = 0.
+
+        base_source = gf.MTSource.from_pyrocko_event(event)
+        base_source.stf = gf.HalfSinusoidSTF(duration=1.0)
+
+        subs = dict(
+            event_name=event.name,
+            event_time=util.time_to_str(event.time))
+
+        problem = CMTProblem(
+            name=self.name_template % subs,
+            apply_balancing_weights=self.apply_balancing_weights,
+            base_source=base_source,
+            targets=targets,
+            ranges=self.ranges,
+            distance_min=self.distance_min,
+            nbootstrap=self.nbootstrap,
+            mt_type=self.mt_type)
+
+        return problem
+
+
+__all__ = [
+    'CMTProblem',
+    'CMTProblemConfig',
+]

src/dataset.py

+
+import logging
+from collections import defaultdict
+import numpy as num
+from pyrocko import util, pile, model, config, trace, snuffling
+from pyrocko.fdsn import enhanced_sacpz, station as fs
+from pyrocko.guts import Object, Tuple, String, Float, dump_all, load_all
+
+
+logger = logging.getLogger('grond.dataset')
+
+
+class InvalidObject(Exception):
+    pass
+
+
+class NotFound(Exception):
+    pass
+
+
+class StationCorrection(Object):
+    codes = Tuple.T(4, String.T())
+    delay = Float.T()
+    factor = Float.T()
+
+
+def load_station_corrections(filename):
+    scs = load_all(filename=filename)
+    for sc in scs:
+        assert isinstance(sc, StationCorrection)
+