wip

apps/grond

@@ -174,20 +174,31 @@ def command_init(args):
 
 def command_check(args):
     def setup(parser):
-        pass
+        parser.add_option(
+            '--event', dest='event_name', metavar='EVENTNAME',
+            help='process only event EVENTNAME')
 
     parser, options, args = cl_parse('check', args, setup)
     if len(args) != 1:
         help_and_die(parser, 'no config file given')
 
+    event_names = None
+    if options.event_name:
+        event_names = [options.event_name]
+
     config_path = args[0]
     config = grond.read_config(config_path)
 
-    grond.check(config)
+    grond.check(
+        config,
+        event_names=event_names)
 
 
 def command_go(args):
     def setup(parser):
+        parser.add_option(
+            '--event', dest='event_name', metavar='EVENTNAME',
+            help='process only event EVENTNAME')
         parser.add_option(
             '--force', dest='force', action='store_true',
             help='overwrite existing run directory')
@@ -209,20 +220,36 @@ def command_go(args):
     else:
         status = tuple(options.status.split(','))
 
+    event_names = None
+    if options.event_name:
+        event_names = [options.event_name]
+
     grond.go(
         config,
+        event_names=event_names,
         force=options.force,
         status=status,
         nparallel=options.nparallel)
 
 
 def command_forward(args):
-    parser, options, args = cl_parse('forward', args)
+    def setup(parser):
+        parser.add_option(
+            '--event', dest='event_name', metavar='EVENTNAME',
+            help='process only event EVENTNAME')
+
+    parser, options, args = cl_parse('forward', args, setup)
     if len(args) != 1:
         help_and_die(parser, 'incorrect number of arguments')
 
+    event_names = None
+    if options.event_name:
+        event_names = [options.event_name]
+
     run_path, = args
-    grond.forward(run_path)
+    grond.forward(
+        run_path,
+        event_names=event_names)
 
 
 def command_harvest(args):

src/cmt.py

@@ -105,7 +105,7 @@ class CMTProblem(core.Problem):
         rm6 = m6 / mt.scalar_moment()
 
         x = num.array([
-            source.time,
+            source.time - self.base_source.time,
             source.north_shift,
             source.east_shift,
             source.depth,
@@ -267,7 +267,7 @@ class CMTProblemConfig(core.ProblemConfig):
             event.depth = 0.
 
         base_source = gf.MTSource.from_pyrocko_event(event)
-        base_source.stf = gf.HalfSinusoidSTF(duration=1.0)
+        base_source.stf = gf.HalfSinusoidSTF(duration=event.duration or 0.0)
 
         subs = dict(
             event_name=event.name,

src/core.py

@@ -12,7 +12,7 @@ import numpy as num
 from pyrocko.guts import load, Object, String, Float, Int, Bool, List, \
     StringChoice, Dict, Timestamp
 from pyrocko import orthodrome as od, gf, trace, guts, util, weeding
-from pyrocko import parimap, model
+from pyrocko import parimap, model, gui_util
 
 from grond import dataset
 
@@ -356,15 +356,17 @@ class MisfitTarget(gf.Target):
                 mv, mn, tr_syn_proc, tr_obs_proc = tr_syn.misfit(
                     tr_obs, ms, nocache=True, debug=True)
 
-            result = MisfitResult(misfit_value=mv, misfit_norm=mn)
+            result = MisfitResult(
+                misfit_value=float(mv), misfit_norm=float(mn))
+
             if self._return_traces:
                 result.filtered_obs = tr_obs
                 result.filtered_syn = tr_syn
                 result.processed_obs = tr_obs_proc
                 result.processed_syn = tr_syn_proc
                 result.taper = ms.taper
-                result.tobs_shift = tobs_shift
-                result.tsyn_pick = tsyn
+                result.tobs_shift = float(tobs_shift)
+                result.tsyn_pick = float(tsyn)
 
             return result
 
@@ -460,23 +462,45 @@ class HasPaths(Object):
                 for p in path]
 
 
+class RandomResponse(trace.FrequencyResponse):
+
+    scale = Float.T(default=0.0)
+
+    def set_random_state(self, rstate):
+        self._rstate = rstate
+
+    def evaluate(self, freqs):
+        n = freqs.size
+        return 1.0 + freqs*(self._rstate.normal(scale=self.scale, size=n) + 
+            0.0J * self._rstate.normal(scale=self.scale, size=n))
+
+
+class SyntheticWaveformNotAvailable(Exception):
+    pass
+
+
 class SyntheticTest(Object):
     random_seed = Int.T(default=0)
     inject_solution = Bool.T(default=False)
+    ignore_data_availability = Bool.T(default=False)
+    add_real_noise = Bool.T(default=False)
+    toffset_real_noise = Float.T(default=-3600.)
     x = Dict.T(String.T(), Float.T())
 
     def __init__(self, **kwargs):
         Object.__init__(self, **kwargs)
         self._synthetics = None
+        self._rstate = num.random.RandomState(self.random_seed)
 
     def set_config(self, config):
         self._config = config
 
     def get_problem(self):
+        raise Exception('TODO: fixme (event_names)')
+
         ds = self._config.get_dataset()
         events = ds.get_events()
-        rstate = num.random.RandomState(self.random_seed)
-        event = events[rstate.randint(0, len(events))]
+        event = events[0]
         return self._config.get_problem(event)
 
     def get_x_random(self):
@@ -484,11 +508,10 @@ class SyntheticTest(Object):
         xbounds = num.array(problem.bounds(), dtype=num.float)
         npar = xbounds.shape[0]
 
-        rstate = num.random.RandomState(self.random_seed)
         x = num.zeros(npar, dtype=num.float)
         while True:
             for i in xrange(npar):
-                x[i] = rstate.uniform(xbounds[i, 0], xbounds[i, 1])
+                x[i] = self._rstate.uniform(xbounds[i, 0], xbounds[i, 1])
 
             try:
                 x = problem.preconstrain(x)
@@ -506,7 +529,13 @@ class SyntheticTest(Object):
                 problem.parameter_array(self.x))
 
         else:
-            x = self.get_x_random()
+            print problem.base_source
+            x = problem.preconstrain(
+                problem.pack(
+                    problem.base_source))
+
+            print x
+            print problem.unpack(x)
 
         return x
 
@@ -527,9 +556,22 @@ class SyntheticTest(Object):
             if result.trace.codes == nslc:
                 tr = result.trace.pyrocko_trace()
                 tr.extend(tmin - tfade * 2.0, tmax + tfade * 2.0)
+                tr2 = tr.copy()
+
+                randomresponse = RandomResponse(scale=10.)
+                randomresponse.set_random_state(self._rstate)
+
                 tr = tr.transfer(tfade=tfade, freqlimits=freqlimits)
+                tr2 = tr2.transfer(tfade=tfade, freqlimits=freqlimits,
+                    transfer_function=randomresponse)
+
                 tr.chop(tmin, tmax)
-                return tr
+                tr2.chop(tmin, tmax)
+
+                #trace.snuffle([tr, tr2])
+                return tr2
+
+        return None
 
 
 class DatasetConfig(HasPaths):
@@ -838,6 +880,9 @@ def stations_mean_latlondist(stations):
 
 def read_config(path):
     config = load(filename=path)
+    if not isinstance(config, Config):
+        raise GrondError('invalid Grond configuration in file "%s"' % path)
+
     config.set_basepath(op.dirname(path) or '.')
     return config
 
@@ -1142,27 +1187,45 @@ def bootstrap_outliers(problem, misfits, std_factor=1.0):
     return num.where(gms > m+s)[0]
 
 
-def forward(rundir):
+def forward(rundir_or_config_path, event_names=None):
 
-    # config = guts.load(filename=op.join('.', 'grond_td.conf'))
-    # config.set_basepath('.')
+    if os.path.isdir(rundir_or_config_path):
+        config = guts.load(filename=op.join(rundir_or_config_path, 'config.yaml'))
+        config.set_basepath(rundir)
+        problem, xs, misfits = load_problem_info_and_data(rundir, subset='harvest')
 
-    config = guts.load(filename=op.join(rundir, 'config.yaml'))
-    config.set_basepath(rundir)
-    ds = config.get_dataset()
+        gms = problem.global_misfits(misfits)
+        ibest = num.argmin(gms)
+        xbest = xs[ibest, :]
 
-    problem, xs, misfits = load_problem_info_and_data(rundir, subset='harvest')
-    for target in problem.targets:
-        target.set_dataset(ds)
+        ds = config.get_dataset()
+        problem.set_engine(config.engine_config.get_engine())
 
-    gms = problem.global_misfits(misfits)
-    isort = num.argsort(gms)
-    gms = gms[isort]
-    xs = xs[isort, :]
+        for target in problem.targets:
+            target.set_dataset(ds)
+
+        payload = [(problem, xbest)]
+
+    else:
+        config = read_config(rundir_or_config_path)
+        ds = config.get_dataset()
+        events = ds.get_events(event_names=event_names)
+
+        payload = []
+        for event in events:
+            problem = config.get_problem(event)
+            xref = problem.preconstrain(
+                problem.pack(problem.base_source))
+            payload.append((problem, xref))
 
     all_trs = []
-    for xbest in xs[:1, :]:
-        ms, ns, results = problem.evaluate(xbest, return_traces=True)
+    events = []
+    for (problem, x) in payload:
+        ds.empty_cache()
+        ms, ns, results = problem.evaluate(x, return_traces=True)
+
+        event = problem.unpack(x).pyrocko_event()
+        events.append(event)
 
         for result in results:
             if result:
@@ -1171,7 +1234,11 @@ def forward(rundir):
                 all_trs.append(result.filtered_obs)
                 all_trs.append(result.filtered_syn)
 
-    trace.snuffle(all_trs)
+    markers = []
+    for ev in events:
+        markers.append(gui_util.EventMarker(ev))
+
+    trace.snuffle(all_trs, markers=markers, stations=ds.get_stations())
 
 
 def harvest(rundir, problem=None, nbest=10, force=False, weed=0):
@@ -1236,9 +1303,9 @@ def check_problem(problem):
 g_state = {}
 
 
-def check(config):
+def check(config, event_names=None):
     ds = config.get_dataset()
-    events = ds.get_events()
+    events = ds.get_events(event_names=event_names)
     nevents = len(events)
 
     from matplotlib import pyplot as plt
@@ -1257,6 +1324,7 @@ def check(config):
             results_list = []
             for i in xrange(10):
                 x = problem.random_uniform(xbounds)
+                print x
                 ms, ns, results = problem.evaluate(x, return_traces=True)
                 results_list.append(results)
 
@@ -1304,7 +1372,10 @@ def check(config):
 
                         axes.plot(xdata, ydata*0.5 + 1.5, color=color)
                         if result.tsyn_pick:
-                            axes.axvline(result.tsyn_pick, color=(0.7, 0.7, 0.7), zorder=2)
+                            axes.axvline(
+                                result.tsyn_pick,
+                                color=(0.7, 0.7, 0.7),
+                                zorder=2)
 
                         t = result.processed_syn.get_xdata()
                         taper = result.taper
@@ -1326,18 +1397,19 @@ def check(config):
                 str(e)))
 
 
-def go(config, force=False, nparallel=1, status=('state',)):
+def go(config, event_names=None, force=False, nparallel=1, status=('state',)):
 
     status = tuple(status)
 
     ds = config.get_dataset()
-    events = ds.get_events()
+    events = ds.get_events(event_names=event_names)
+
     nevents = len(events)
 
     if nevents == 0:
         raise GrondError('no events found')
 
-    g_data = (config, force, status, nparallel)
+    g_data = (config, force, status, nparallel, event_names)
 
     g_state[id(g_data)] = g_data
 
@@ -1352,14 +1424,14 @@ def go(config, force=False, nparallel=1, status=('state',)):
 
 def process_event(ievent, g_data_id):
 
-    config, force, status, nparallel = g_state[g_data_id]
+    config, force, status, nparallel, event_names = g_state[g_data_id]
 
     if nparallel > 1:
         status = ()
 
     ds = config.get_dataset()
 
-    events = ds.get_events()
+    events = ds.get_events(event_names=event_names)
     nevents = len(events)
 
     event = events[ievent]

src/dataset.py

@@ -315,7 +315,13 @@ class Dataset(object):
 
         return trs
 
-    def get_waveform_raw(self, obj, tmin=None, tmax=None, tpad=0.):
+    def get_waveform_raw(
+            self, obj,
+            tmin=None,
+            tmax=None,
+            tpad=0.,
+            toffset_noise_extract=0.):
+
         net, sta, loc, cha = self.get_nslc(obj)
 
         if self.is_blacklisted((net, sta, loc, cha)):
@@ -337,10 +343,16 @@ class Dataset(object):
                     'waveform clipped', (net, sta, loc, cha))
 
         trs = self.pile.all(
-            tmin=tmin, tmax=tmax, tpad=tpad,
+            tmin=tmin+toffset_noise_extract,
+            tmax=tmax+toffset_noise_extract,
+            tpad=tpad,
             trace_selector=lambda tr: tr.nslc_id == (net, sta, loc, cha),
             want_incomplete=False)
 
+        if toffset_noise_extract != 0.0:
+            for tr in trs:
+                tr.shift(-toffset_noise_extract)
+
         if len(trs) == 1:
             return trs[0]
 
@@ -352,15 +364,21 @@ class Dataset(object):
             self,
             obj, quantity='displacement',
             tmin=None, tmax=None, tpad=0.,
-            tfade=0., freqlimits=None, deltat=None):
+            tfade=0., freqlimits=None, deltat=None,
+            toffset_noise_extract=0.):
+
+        assert quantity == 'displacement'  # others not yet implemented
 
-        tr = self.get_waveform_raw(obj, tmin=tmin, tmax=tmax, tpad=tpad+tfade)
+        tr = self.get_waveform_raw(
+            obj, tmin=tmin, tmax=tmax, tpad=tpad+tfade,
+            toffset_noise_extract=toffset_noise_extract)
 
         if deltat is not None:
             tr.downsample_to(deltat, snap=True)
             tr.deltat = deltat
 
         resp = self.get_response(tr)
+        print resp
         return tr.transfer(tfade=tfade, freqlimits=freqlimits,
                            transfer_function=resp, invert=True)
 
@@ -373,6 +391,8 @@ class Dataset(object):
             source=None,
             target=None):
 
+        assert quantity == 'displacement'  # others not yet implemented
+
         if cache is True:
             cache = self._cache
 
@@ -394,14 +414,26 @@ class Dataset(object):
             else:
                 return obj
 
-        if self.synthetic_test:
-            tr = self.synthetic_test.get_waveform(
-                nslc, tmin, tmax,
-                tfade=tfade, freqlimits=freqlimits)
-            if cache is not None:
-                cache[tr.nslc_id, tmin, tmax] = tr
+        syn_test = self.synthetic_test
+        toffset_noise_extract = 0.0
+        if syn_test:
+            if syn_test.ignore_data_availability:
+                if syn_test.add_real_noise:
+                    raise DatasetError(
+                        'ignore_data_availability=True and '
+                        'add_real_noise=True cannot be combined.')
 
-            return tr
+                tr = syn_test.get_waveform(
+                    nslc, tmin, tmax,
+                    tfade=tfade, freqlimits=freqlimits)
+
+                if cache is not None:
+                    cache[tr.nslc_id, tmin, tmax] = tr
+
+                return tr
+
+            if syn_test.add_real_noise:
+                toffset_noise_extract = syn_test.toffset_real_noise
 
         abs_delays = []
         for ocha in 'ENZRT':
@@ -438,6 +470,7 @@ class Dataset(object):
                         trs.append(self.get_waveform_restituted(
                             station.nsl() + (cha,),
                             tmin=tmin, tmax=tmax, tpad=tpad+abs_delay_max,
+                            toffset_noise_extract=toffset_noise_extract,
                             tfade=tfade, freqlimits=freqlimits, deltat=deltat))
 
                     trs_projected.extend(
@@ -455,6 +488,22 @@ class Dataset(object):
                 if tmin is not None and tmax is not None:
                     tr.chop(tmin, tmax)
 
+            if syn_test:
+                trs_projected_synthetic = []
+                for tr in trs_projected:
+                    tr_syn = syn_test.get_waveform(
+                        tr.nslc_id, tmin, tmax,
+                        tfade=tfade, freqlimits=freqlimits)
+
+                    if tr_syn:
+                        if syn_test.add_real_noise:
+                            tr_syn = tr_syn.copy()
+                            tr_syn.add(tr)
+
+                        trs_projected_synthetic.append(tr_syn)
+
+                trs_projected = trs_projected_synthetic
+
             if cache is not None:
                 for tr in trs_projected:
                     cache[tr.nslc_id, tmin, tmax] = tr
@@ -469,10 +518,11 @@ class Dataset(object):
             cache[nslc, tmin, tmax] = e
             raise
 
-    def get_events(self, magmin=None):
+    def get_events(self, magmin=None, event_names=None):
         evs = []
         for ev in self.events:
-            if magmin is None or ev.magnitude >= magmin:
+            if ((magmin is None or ev.magnitude >= magmin) and
+                    (event_names is None or ev.name in event_names)):
                 evs.append(ev)
 
         return evs

src/plot.py

@@ -333,6 +333,14 @@ def draw_jointpar_figures(
     xs = model.xs
 
     bounds = problem.bounds() + problem.dependant_bounds()
+    for ipar in xrange(problem.ncombined):
+        par = problem.combined[ipar]
+        lo, hi = bounds[ipar]
+        if lo == hi:
+            if exclude is None:
+                exclude = []
+
+            exclude.append(par.name)
 
     xref = problem.pack(problem.base_source)
 
@@ -504,10 +512,10 @@ def draw_jointpar_figures(
                 xpar.scaled(fx), ypar.scaled(fy), 's',
                 mew=1.5, ms=5, color=ref_color_light, mec=ref_color)
 
-    for jfig, figs_row in enumerate(figs):
-        for ifig, fig in enumerate(figs_row):
-            if fig is not None:
-                fig.savefig('jointpar-%i-%i.pdf' % (jfig, ifig))
+    #for jfig, figs_row in enumerate(figs):
+    #    for ifig, fig in enumerate(figs_row):
+    #        if fig is not None:
+    #            fig.savefig('jointpar-%i-%i.pdf' % (jfig, ifig))
 
 
 def draw_solution_figure(
@@ -854,6 +862,8 @@ def draw_fits_figures(ds, model, plt):
     dtraces = []
     for target, result in zip(problem.targets, results):
         if result is None:
+            print 'xxx'
+            print target
             dtraces.append(None)
             continue
 
@@ -861,6 +871,7 @@ def draw_fits_figures(ds, model, plt):
         w = target.get_combined_weight(problem.apply_balancing_weights)
 
         if target.misfit_config.domain != 'time_domain':
+            dtraces.append(None)
             continue
 
         for tr in (
@@ -879,12 +890,12 @@ def draw_fits_figures(ds, model, plt):
                 (result.processed_syn.get_ydata() -
                  result.processed_obs.get_ydata())**2))
         dtraces.append(dtrace)
-
         all_syn_trs.append(result.processed_syn)
 
     amin, amax = trace.minmax(all_syn_trs, lambda tr: None)[None]
 
-    dmin, dmax = trace.minmax([x for x in dtraces if x], lambda tr: None)[None]
+    dmin, dmax = trace.minmax(
+        [x for x in dtraces if x is not None], lambda tr: None)[None]
 
     for tr in dtraces:
         if tr:
@@ -970,6 +981,7 @@ def draw_fits_figures(ds, model, plt):
                 if (iy, ix) not in frame_to_target:
                     continue
 
+
                 ixx = ix/nxmax
                 iyy = iy/nymax
                 if (iyy, ixx) not in figures:
@@ -978,6 +990,11 @@ def draw_fits_figures(ds, model, plt):
                 fig = figures[iyy, ixx]
 
                 target = frame_to_target[iy, ix]
+                print target
+                print target.misfit_config.domain
+
+                if target.misfit_config.domain != 'time_domain':
+                    continue
 
                 ny_this = min(ny, nymax)
                 nx_this = min(nx, nxmax)
@@ -1175,7 +1192,7 @@ def draw_hudson_figure(model, plt):
         position=(u, v),
         size=beachballsize,
         color_t=color,
-        zorder=-1,
+        zorder=2,
         linewidth=0.5)
 
     mt = best_source.pyrocko_moment_tensor()
@@ -1187,7 +1204,19 @@ def draw_hudson_figure(model, plt):
         mew=1,
         mec='black',
         color='none',
-        zorder=1)
+        zorder=-2)
+
+    mt = problem.base_source.pyrocko_moment_tensor()
+    u, v = hudson.project(mt)
+
+    beachball.plot_beachball_mpl(
+        mt, axes,
+        beachball_type='dc',
+        position=(u, v),
+        size=beachballsize,
+        color_t='red',
+        zorder=2,
+        linewidth=0.5)
 
     fig.savefig('hudson.pdf')
 
@@ -1246,6 +1275,7 @@ def plot_result(dirname, plotnames_want):
             if plotname in plotnames_want: