grond harvest: allow exporting of target results with --export-fits=...

src/apps/grond.py

@@ -807,17 +807,35 @@ def command_harvest(args):
                  'average misfit of all NEACH best in all chains, '
                  '3: harvesting is done on the global chain only, bootstrap '
                  'chains are excluded')
+        parser.add_option(
+            '--export-fits', dest='export_fits', default='',
+            help='additionally export details about the fit of individual '
+                 'targets. "best" - export fits of best model, "mean" - '
+                 'export fits of ensemble mean model, "ensemble" - export '
+                 'fits of all models in harvest ensemble.')
 
     parser, options, args = cl_parse('harvest', args, setup)
-    if len(args) != 1:
+    if len(args) < 1:
         help_and_die(parser, 'no rundir')
 
-    run_path, = args
-    grond.harvest(
-        run_path,
-        force=options.force,
-        nbest=options.neach,
-        weed=options.weed)
+    export_fits = []
+    if options.export_fits.strip():
+        export_fits = [x.strip() for x in options.export_fits.split(',')]
+
+    for run_path in args:
+        try:
+            grond.harvest(
+                run_path,
+                force=options.force,
+                nbest=options.neach,
+                weed=options.weed,
+                export_fits=export_fits)
+
+        except grond.DirectoryAlreadyExists as e:
+            die(str(e) + '\n    Use --force to overwrite.')
+
+        except grond.GrondError as e:
+            die(str(e))
 
 
 def command_cluster(args):

src/core.py

@@ -21,6 +21,8 @@ from .problems.base import Problem, load_problem_info_and_data, \
 
 from .optimisers.base import BadProblem
 from .targets.waveform.target import WaveformMisfitResult
+from .targets.base import dump_misfit_result_collection, \
+    MisfitResultCollection, MisfitResult, MisfitResultError
 from .meta import expand_template, GrondError, selected
 from .environment import Environment
 from .monitor import GrondMonitor
@@ -65,7 +67,7 @@ def lock_rundir(rundir):
         os.remove(statefn)
 
 
-class DirectoryAlreadyExists(Exception):
+class DirectoryAlreadyExists(GrondError):
     pass
 
 
@@ -146,7 +148,9 @@ def forward(env, show='filtered'):
     trace.snuffle(all_trs, markers=markers, stations=list(stations.values()))
 
 
-def harvest(rundir, problem=None, nbest=10, force=False, weed=0):
+def harvest(
+        rundir, problem=None, nbest=10, force=False, weed=0,
+        export_fits=[]):
 
     env = Environment([rundir])
     optimiser = env.get_optimiser()
@@ -166,7 +170,8 @@ def harvest(rundir, problem=None, nbest=10, force=False, weed=0):
         if force:
             shutil.rmtree(dumpdir)
         else:
-            raise DirectoryAlreadyExists(dumpdir)
+            raise DirectoryAlreadyExists(
+                'Harvest directory already exists: %s' % dumpdir)
 
     util.ensuredir(dumpdir)
 
@@ -205,6 +210,36 @@ def harvest(rundir, problem=None, nbest=10, force=False, weed=0):
     for i in ibests:
         problem.dump_problem_data(dumpdir, xs[i], misfits[i, :, :])
 
+    if export_fits:
+        env.setup_modelling()
+        problem = env.get_problem()
+        history = env.get_history(subset='harvest')
+
+        for what in export_fits:
+            if what == 'best':
+                models = [history.get_best_model()]
+            elif what == 'mean':
+                models = [history.get_mean_model()]
+            elif what == 'ensemble':
+                models = history.models
+            else:
+                raise GrondError(
+                    'Invalid option for harvest\'s export_fits argument: %s'
+                    % what)
+
+            results = []
+            for x in models:
+                results.append([
+                    (result if isinstance(result, MisfitResult)
+                     else MisfitResultError(message=str(result))) for
+                    result in problem.evaluate(x)])
+
+            emr = MisfitResultCollection(results=results)
+
+            dump_misfit_result_collection(
+                emr,
+                op.join(dumpdir, 'fits-%s.yaml' % what))
+
     logger.info('Done harvesting problem "%s".' % problem.name)
 
 
@@ -799,6 +834,7 @@ def export(
 
 
 __all__ = '''
+    DirectoryAlreadyExists
     forward
     harvest
     cluster

src/targets/base.py

@@ -4,10 +4,10 @@ import numpy as num
 
 from pyrocko import gf
 from pyrocko.guts_array import Array
-from pyrocko.guts import Object, Float, Dict
+from pyrocko.guts import Object, Float, String, Dict, List, Choice, load, dump
 
 from grond.analysers.base import AnalyserResult
-from grond.meta import has_get_plot_classes
+from grond.meta import has_get_plot_classes, GrondError
 
 
 guts_prefix = 'grond'
@@ -188,8 +188,41 @@ class MisfitTarget(Object):
         return modelling_results[0]
 
 
+class MisfitResultError(Object):
+    message = String.T()
+
+
+class MisfitResultCollection(Object):
+    results = List.T(List.T(
+        Choice.T([MisfitResult.T(), MisfitResultError.T()])))
+
+
+def dump_misfit_result_collection(misfit_result_collection, path):
+    dump(misfit_result_collection, filename=path)
+
+
+def load_misfit_result_collection(path):
+    try:
+        obj = load(filename=path)
+
+    except OSError as e:
+        raise GrondError(
+            'Failed to read ensemble misfit results from file "%s" (%s)' % (
+                path, e))
+
+    if not isinstance(obj, MisfitResultCollection):
+        raise GrondError(
+            'File "%s" does not contain any misfit result collection.' % path)
+
+    return obj
+
+
 __all__ = '''
     TargetGroup
     MisfitTarget
     MisfitResult
+    MisfitResultError
+    dump_misfit_result_collection
+    load_misfit_result_collection
+    MisfitResultCollection
 '''.split()

src/targets/satellite/plot.py

@@ -347,7 +347,7 @@ Surface displacements derived from satellite data.
 data and (right) the model residual.
 ''')
 
-            stat_obs = result.statics_obs
+            stat_obs = result.statics_obs['displacement.los']
             stat_syn = result.statics_syn['displacement.los']
             res = stat_obs - stat_syn
 

src/targets/satellite/target.py

@@ -4,6 +4,7 @@ from scipy import linalg as splinalg
 
 from pyrocko import gf
 from pyrocko.guts import String, Bool, Dict, List
+from pyrocko.guts_array import Array
 
 import os
 
@@ -107,9 +108,13 @@ class SatelliteTargetGroup(TargetGroup):
 class SatelliteMisfitResult(gf.Result, MisfitResult):
     """Carries the observations for a target and corresponding synthetics."""
     statics_syn = Dict.T(
+        String.T(),
+        Array.T(dtype=num.float, shape=(None,), serialize_as='base64'),
         optional=True,
         help='Predicted static displacements for a target (synthetics).')
     statics_obs = Dict.T(
+        String.T(),
+        Array.T(dtype=num.float, shape=(None,), serialize_as='base64'),
         optional=True,
         help='Observed static displacement for a target.')
 
@@ -218,7 +223,7 @@ class SatelliteMisfitTarget(gf.SatelliteTarget, MisfitTarget):
 
         if self._result_mode == 'full':
             result.statics_syn = statics
-            result.statics_obs = quadtree.leaf_medians
+            result.statics_obs = {'displacement.los': obs}
 
         return result
 

src/targets/waveform/target.py

@@ -110,10 +110,6 @@ class DomainChoice(StringChoice):
         'cc_max_norm']
 
 
-class Trace(Object):
-    pass
-
-
 class WaveformMisfitConfig(MisfitConfig):
     quantity = gf.QuantityType.T(default='displacement')
     fmin = Float.T(default=0.0, help='minimum frequency of bandpass filter')
@@ -359,10 +355,10 @@ class WaveformMisfitResult(gf.Result, MisfitResult):
 
     A number of different waveform  or phase representations are possible.
     '''
-    processed_obs = Trace.T(optional=True)
-    processed_syn = Trace.T(optional=True)
-    filtered_obs = Trace.T(optional=True)
-    filtered_syn = Trace.T(optional=True)
+    processed_obs = trace.Trace.T(optional=True)
+    processed_syn = trace.Trace.T(optional=True)
+    filtered_obs = trace.Trace.T(optional=True)
+    filtered_syn = trace.Trace.T(optional=True)
     spectrum_obs = TraceSpectrum.T(optional=True)
     spectrum_syn = TraceSpectrum.T(optional=True)
 
@@ -370,7 +366,7 @@ class WaveformMisfitResult(gf.Result, MisfitResult):
     tobs_shift = Float.T(optional=True)
     tsyn_pick = Timestamp.T(optional=True)
     tshift = Float.T(optional=True)
-    cc = Trace.T(optional=True)
+    cc = trace.Trace.T(optional=True)
 
     piggyback_subresults = List.T(WaveformPiggybackSubresult.T())
 

test/test_examples.py

@@ -87,5 +87,6 @@ def run_example(project_name, config_path, quick_config_path, event_name):
         mod_conf.set_basepath(conf.get_basepath())
         config.write_config(mod_conf, quick_config_path)
         grond('go', quick_config_path, event_name)
+        grond('harvest', '--force', '--export-fits=best,mean', rundir_path)
         grond('report', rundir_path)
         # assert os.path.isdir(os.path.join('report', event_name))