Merge branch 'refactor' into solver_plot

Conflicts:
	src/core.py
	src/plot.py
	src/problems.py
	src/targets.py

apps/grond

@@ -218,7 +218,8 @@ def command_init(args):
             group='all',
             interpolation='multilinear',
             store_id='gf_store',
-            inner_satellite_misfit_config=grond.InnerSatelliteMisfitConfig(
+            kite_scenes=['*all'],
+            inner_misfit_config=grond.InnerSatelliteMisfitConfig(
                 use_weight_focal=False))]
 
         problem_config = grond.RectangularProblemConfig(

setup.py

 #!/usr/bin/env python
-
+import os
 from distutils.core import setup
 
+
+def grond_completion():
+    if os.access('/etc/bash_completion.d/', os.W_OK):
+        return [('/etc/bash_completion.d', ['extras/grond'])]
+    return []
+
+
 setup(
     name='grond',
     description='What do you want to bust today?!',
@@ -11,8 +18,7 @@ setup(
     packages=['grond', 'grond.baraddur'],
     scripts=['apps/grond'],
     package_dir={'grond': 'src'},
-    package_data={'grond': [],
-                  'grond': ['baraddur/templates/*.html',
+    package_data={'grond': ['baraddur/templates/*.html',
                             'baraddur/res/*']},
-    data_files=[('/etc/bash_completion.d', ['extras/grond'])],
+    data_files=[] + grond_completion(),
     )

src/baraddur/meta.py

+import numpy as num
+
+
+class ColorCycler(list):
+
+    def __init__(self, *args, **kwargs):
+        list.__init__(self, *args, **kwargs)
+        self.index = -1
+
+    def next(self):
+        self.index += 1
+        if self.index >= len(self):
+            self.index = 0
+        return self[self.index]
+
+
+def makeColorGradient(misfits, fr=1., fg=.5, fb=1.,
+                      pr=0, pg=2.5, pb=4):
+    misfits /= misfits.max()
+    r = num.sin(fr * misfits + pr) * 127 + 128
+    g = num.sin(fg * misfits + pg) * 127 + 128
+    b = num.sin(fb * misfits + pb) * 127 + 128
+    return ['#%02x%02x%02x' % (r[i], g[i], b[i]) for i in xrange(misfits.size)]

src/baraddur/server.py

@@ -7,6 +7,8 @@ import logging
 import numpy as num
 import socket
 
+from scipy import ndimage
+
 from collections import OrderedDict
 from datetime import datetime
 
@@ -24,52 +26,16 @@ from bokeh.models import ColumnDataSource
 from bokeh import layouts
 from bokeh.plotting import figure
 
+from .meta import ColorCycler
+
 logging.basicConfig(level=logging.INFO)
 logger = logging.getLogger('grond.baraddur')
 
 
-def makeColorGradient(misfits, fr=1., fg=.5, fb=1.,
-                      pr=0, pg=2.5, pb=4):
-    misfits /= misfits.max()
-    r = num.sin(fr * misfits + pr) * 127 + 128
-    g = num.sin(fg * misfits + pg) * 127 + 128
-    b = num.sin(fb * misfits + pb) * 127 + 128
-    return ['#%02x%02x%02x' % (r[i], g[i], b[i]) for i in xrange(misfits.size)]
-
-
-class BaraddurRequestHandler(RequestHandler):
-    def initialize(self, config):
-        self.config = config
-        self.model = config.model
-
-
-class BaraddurBokehHandler(BokehHandler):
-    def __init__(self, config, *args, **kwargs):
-        BokehHandler.__init__(self, *args, **kwargs)
-        self.config = config
-        self.model = config.model
-        self.nmodels = 0
-
-        self.source = ColumnDataSource()
-
-    def update_source(self, func):
-        new_nmodels, new_data = func(
-            skip_nmodels=self.nmodels,
-            keys=self.source.data.keys())
-        self.source.stream(new_data)
-        self.nmodels += new_nmodels
-
-    @gen.coroutine
-    def update_models(self):
-        return self.update_source(self.model.get_models)
-
-    @gen.coroutine
-    def update_misfits(self):
-        return self.update_source(self.model.get_misfits)
-
-
 class BaraddurModel(object):
     def __init__(self, rundir):
+        if rundir is None:
+            return
         logger.debug('Loading problem from %s' % rundir)
         self.rundir = op.abspath(rundir)
         self.problem = grond.core.load_problem_info(self.rundir)
@@ -184,6 +150,77 @@ class BaraddurModel(object):
         return nmodels, mf_dict
 
 
+# This should be realized in JS!
+class SmoothColumnDataSource(ColumnDataSource):
+    '''Gauss smoothed ColumnDataSource '''
+    def __init__(self, *args, **kwargs):
+        self._gaussian_kw = {
+            'sigma': 3,
+        }
+        for kw in ['sigma', 'order', 'output', 'mode', 'cval', 'truncate']:
+            if kw in kwargs.keys():
+                self._gaussian_kw[kw] = kwargs.pop(kw)
+        ColumnDataSource.__init__(self, *args, **kwargs)
+
+    def stream(self, data, **kwargs):
+        sigma = self._gaussian_kw['sigma']
+        if data.keys() == self.data.keys():
+            raise AttributeError('streaming data must represent'
+                                 ' existing column')
+        for key, arr in data.iteritems():
+            if arr.ndim > 1 or not isinstance(arr, num.ndarray):
+                raise AttributeError('data is not numpy.ndarray of 1d')
+
+            nvalues = arr.size
+            if nvalues < sigma:
+                self.data[key][-(sigma-nvalues):]
+                d = num.empty(sigma)
+                d[:(sigma-nvalues)] = self.data[key][-(sigma-nvalues):]
+                d[(sigma-nvalues):] = arr[:]
+                arr = ndimage.gaussian_filter1d(d, **self._gaussian_kw)
+                arr = arr[-nvalues:]
+            else:
+                arr = ndimage.gaussian_filter1d(arr, **self._gaussian_kw)
+            data[key] = arr
+
+        ColumnDataSource.stream(self, data, **kwargs)
+
+
+class BaraddurRequestHandler(RequestHandler):
+    def initialize(self, config):
+        self.config = config
+        self.model = config.model
+
+
+class BaraddurBokehHandler(BokehHandler):
+    def __init__(self, config, *args, **kwargs):
+        BokehHandler.__init__(self, *args, **kwargs)
+        self.config = config
+        self.model = config.model
+
+    @staticmethod
+    def update_source(doc, func):
+        new_nmodels, new_data = func(
+            skip_nmodels=doc.nmodels,
+            keys=doc.source.data.keys())
+        doc.source.stream(new_data)
+        doc.nmodels += new_nmodels
+
+    def update_models(self, doc):
+        return self.update_source(doc, self.model.get_models)
+
+    def update_misfits(self, doc):
+        return self.update_source(doc, self.model.get_misfits)
+
+    @staticmethod
+    def decorate_document(modify_document):
+        def decorator(handler, doc):
+            doc.nmodels = 0
+            doc.source = ColumnDataSource()
+            return modify_document(handler, doc)
+        return decorator
+
+
 class Rundirs(BaraddurRequestHandler):
 
     def get(self):
@@ -215,12 +252,13 @@ class Misfit(BaraddurRequestHandler):
 
     class MisfitsPlot(BaraddurBokehHandler):
 
+        @BaraddurBokehHandler.decorate_document
         def modify_document(self, doc):
-            self.source.add(num.ndarray(0, dtype=num.float32),
-                            'misfit_mean')
-            self.source.add(num.ndarray(0, dtype=num.float32),
-                            'niter')
-            self.update_misfits()
+            doc.source.add(num.ndarray(0, dtype=num.float32),
+                           'misfit_mean')
+            doc.source.add(num.ndarray(0, dtype=num.float32),
+                           'niter')
+            self.update_misfits(doc)
 
             plot = figure(webgl=True,
                           responsive=True,
@@ -228,10 +266,10 @@ class Misfit(BaraddurRequestHandler):
                           x_axis_label='Iteration #',
                           y_axis_label='Misfit')
             plot.scatter('niter', 'misfit_mean',
-                         source=self.source, alpha=.4)
+                         source=doc.source, alpha=.4)
 
             doc.add_root(plot)
-            doc.add_periodic_callback(self.update_misfits,
+            doc.add_periodic_callback(lambda: self.update_misfits(doc),
                                       self.config.plot_update_small)
 
     bokeh_handlers = {'misfit': MisfitsPlot}
@@ -253,27 +291,27 @@ class Parameters(BaraddurRequestHandler):
 
         ncols = 4
 
+        @BaraddurBokehHandler.decorate_document
         def modify_document(self, doc):
             problem = self.model.problem
+            self.model = BaraddurModel(self.config.rundir)
 
             for param in ['niter'] + [p.name for p in problem.parameters]:
-                self.source.add(
+                doc.source.add(
                     num.ndarray(0, dtype=num.float32),
                     param)
-            self.model = BaraddurModel(self.config.rundir)
-            self.update_models()
+            self.update_models(doc)
 
             plots = []
             for par in problem.parameters:
                 fig = figure(webgl=True,
-                             responsive=True,
                              x_axis_label='Iteration #',
-                             y_axis_label='%s [%s]' % (par.label, par.unit))
+                             y_axis_label='%s [%s]' % (par.name, par.unit))
                 fig.scatter('niter', par.name,
-                            source=self.source, alpha=.4)
+                            source=doc.source, alpha=.4)
                 plots.append(fig)
-            if (len(plots) % self.ncols) != 0:
-                plots += [None] * (self.ncols - (len(plots) % self.ncols))
+            # if (len(plots) % self.ncols) != 0:
+            #     plots += [None] * (self.ncols - (len(plots) % self.ncols))
 
             grid = layouts.gridplot(
                 plots,
@@ -281,7 +319,7 @@ class Parameters(BaraddurRequestHandler):
                 ncols=self.ncols)
 
             doc.add_root(grid)
-            doc.add_periodic_callback(self.update_models,
+            doc.add_periodic_callback(lambda: self.update_models(doc),
                                       self.config.plot_update_small)
 
     bokeh_handlers = {'parameters': ParameterPlots}
@@ -302,6 +340,7 @@ class Targets(BaraddurRequestHandler):
 
     class TargetContributionPlot(BaraddurBokehHandler):
 
+        @BaraddurBokehHandler.decorate_document
         def modify_document(self, doc):
             plot = figure(webgl=True,
                           responsive=True,
@@ -311,15 +350,16 @@ class Targets(BaraddurRequestHandler):
 
             target_ids = [t.id for t in self.model.targets]
             for param in ['niter'] + target_ids:
-                self.source.add(
+                doc.source.add(
                     num.ndarray(0, dtype=num.float32),
                     param)
-            self.update_misfits()
-            colors = ['blue', 'red']
+            self.update_misfits(doc)
+
+            colors = ColorCycler(['blue', 'red', 'green'])
             for it, target in enumerate(self.model.targets):
                 plot.scatter('niter', target.id,
-                             source=self.source,
-                             color=colors[it])
+                             source=doc.source,
+                             color=colors.next())
 
             # plot.patches(
             #     ['niter'] * self.model.ntargets,
@@ -327,7 +367,7 @@ class Targets(BaraddurRequestHandler):
             #     source=self.source)
 
             doc.add_root(plot)
-            doc.add_periodic_callback(self.update_misfits,
+            doc.add_periodic_callback(lambda: self.update_misfits(doc),
                                       self.config.plot_update_small)
 
     bokeh_handlers = {'targets': TargetContributionPlot}

src/core.py

@@ -300,14 +300,15 @@ def analyse(problem, niter=1000, show_progress=False):
         return
 
     wtargets = []
-    for target in problem.targets:
+    if not problem.has_waveforms:
+        return
+
+    for target in problem.waveform_targets:
         wtarget = copy.copy(target)
         wtarget.flip_norm = True
         wtarget.weight = 1.0
         wtargets.append(wtarget)
 
-    groups, ngroups = problem.get_group_mask()
-
     wproblem = problem.copy()
     wproblem.targets = wtargets
 
@@ -339,7 +340,6 @@ def analyse(problem, niter=1000, show_progress=False):
             isok_mask = num.logical_not(isbad_mask)
         else:
             isok_mask = None
-
         _, ms = wproblem.evaluate(x, mask=isok_mask)
         mss[iiter, :] = ms
 
@@ -351,15 +351,16 @@ def analyse(problem, niter=1000, show_progress=False):
     if show_progress:
         pbar.finish()
 
-    # mean_ms = num.mean(mss, axis=0)
-    # weights = 1.0 / mean_ms
-    weights = num.ones(wproblem.ntargets)
+    mean_ms = num.mean(mss, axis=0)
+    weights = 1. / mean_ms
+    groups, ngroups = wproblem.get_group_mask()
+
     for igroup in xrange(ngroups):
         weights[groups == igroup] /= (
             num.nansum(weights[groups == igroup]) /
             num.nansum(num.isfinite(weights[groups == igroup])))
 
-    for weight, target in zip(weights, problem.targets):
+    for weight, target in zip(weights, problem.waveform_targets):
         target.analysis_result = TargetAnalysisResult(
             balancing_weight=float(weight))
 
@@ -422,7 +423,7 @@ def solve(problem,
           plot=None):
 
     xbounds = num.array(problem.get_parameter_bounds(), dtype=num.float)
-    npar = xbounds.shape[0]
+    npar = problem.nparameters
 
     nlinks_cap = int(round(chain_length_factor * npar + 1))
     chains_m = num.zeros((1 + problem.nbootstrap, nlinks_cap), num.float)
@@ -447,7 +448,7 @@ def solve(problem,
     isbad_mask = None
     accept_sum = num.zeros(1 + problem.nbootstrap, dtype=num.int)
     accept_hist = num.zeros(niter, dtype=num.int)
-    pnames = [p.name for p in problem.parameters]
+    pnames = problem.parameter_names
 
     if plot:
         plot.start(problem)

src/dataset.py

@@ -345,7 +345,7 @@ class Dataset(object):
             return self.kite_scenes[0]
         else:
             for scene in self.kite_scenes:
-                if scene.meta.scene_id is scene_id:
+                if scene.meta.scene_id == scene_id:
                     return scene
         raise NotFound('No kite scene with id %s defined' % scene_id)
 
@@ -843,9 +843,13 @@ class DatasetConfig(HasPaths):
                 ds.add_waveforms(paths=fp(self.waveform_paths))
 
             if self.kite_scene_paths:
+                logger.info('Loading kite scenes...')
                 for path in self.kite_scene_paths:
                     for fn in glob.glob(xjoin(fp(path), '*.npz')):
                         ds.add_kite_scene(filename=fn)
+                if not ds.kite_scenes:
+                    logger.warning('Could not find any kite scenes at %s' %
+                                   self.kite_scene_paths)
 
             if self.clippings_path:
                 ds.add_clippings(markers_filename=fp(self.clippings_path))

src/meta.py

 import os.path as op
 from string import Template
-from pyrocko.guts import Object, String
+from pyrocko.guts import Object, String, Float
 
 
 def xjoin(basepath, path):
@@ -42,6 +42,89 @@ def expand_template(template, d):
             'malformed placeholder in template: "%s"' % template)
 
 
+class ADict(dict):
+    def __getattr__(self, k):
+        return self[k]
+
+    def __setattr__(self, k, v):
+        self[k] = v
+
+
+class Parameter(Object):
+    name__ = String.T()
+    unit = String.T(optional=True)
+    scale_factor = Float.T(default=1., optional=True)
+    scale_unit = String.T(optional=True)
+    label = String.T(optional=True)
+
+    def __init__(self, *args, **kwargs):
+        if len(args) >= 1:
+            kwargs['name'] = args[0]
+        if len(args) >= 2:
+            kwargs['unit'] = args[1]
+
+        self.groups = [None]
+        self._name = None
+
+        Object.__init__(self, **kwargs)
+
+    def get_label(self, with_unit=True):
+        l = [self.label or self.name]
+        if with_unit:
+            unit = self.get_unit_label()
+            if unit:
+                l.append('[%s]' % unit)
+
+        return ' '.join(l)
+
+    def set_groups(self, groups):
+        if not isinstance(groups, list):
+            raise AttributeError('Groups must be a list of strings.')
+        self.groups = groups
+
+    def _get_name(self):
+        if None not in self.groups:
+            return '%s:%s' % (':'.join(self.groups), self._name)
+        return self._name
+
+    def _set_name(self, value):
+        self._name = value
+
+    name = property(_get_name, _set_name)
+
+    @property
+    def name_nogroups(self):
+        return self._name
+
+    def get_value_label(self, value, format='%(value)g%(unit)s'):
+        value = self.scaled(value)
+        unit = self.get_unit_suffix()
+        return format % dict(value=value, unit=unit)
+
+    def get_unit_label(self):
+        if self.scale_unit is not None:
+            return self.scale_unit
+        elif self.unit:
+            return self.unit
+        else:
+            return None
+
+    def get_unit_suffix(self):
+        unit = self.get_unit_label()
+        if not unit:
+            return ''
+        else:
+            return ' %s' % unit
+
+    def scaled(self, x):
+        if isinstance(x, tuple):
+            return tuple(v/self.scale_factor for v in x)
+        if isinstance(x, list):
+            return list(v/self.scale_factor for v in x)
+        else:
+            return x/self.scale_factor
+
+
 class HasPaths(Object):
     path_prefix = Path.T(optional=True)
 

src/plot.py

@@ -10,7 +10,7 @@ from pyrocko import beachball, guts, trace, util, gf
 from pyrocko import hudson
 from grond import core
 from matplotlib import pyplot as plt
-from matplotlib import cm, patches
+from matplotlib import cm, patches, gridspec
 from pyrocko.cake_plot import colors, \
     str_to_mpl_color as scolor, light
 
@@ -80,6 +80,18 @@ def str_duration(t):
         return s + '%.1f d' % (t / (24.*3600.))
 
 
+def scale_axes(ax, scale):
+    from matplotlib.ticker import ScalarFormatter
+
+    class FormatScaled(ScalarFormatter):
+        @staticmethod
+        def __call__(value, pos):
+            return u'%d' % (value * scale)
+
+    ax.get_xaxis().set_major_formatter(FormatScaled())
+    ax.get_yaxis().set_major_formatter(FormatScaled())
+
+
 def eigh_sorted(mat):
     evals, evecs = num.linalg.eigh(mat)
     iorder = num.argsort(evals)
@@ -929,17 +941,137 @@ def plot_dtrace_vline(axes, t, space, **kwargs):
     axes.plot([t, t], [-1.0 - space, -1.0], **kwargs)
 
 
+def draw_fits_figures_statics(ds, model, plt):
+    from pyrocko.orthodrome import latlon_to_ne_numpy
+    problem = model.problem
+
+    for target in problem.targets:
+        target.set_dataset(ds)
+
+    gms = problem.global_misfits(model.misfits)
+    isort = num.argsort(gms)
+    gms = gms[isort]
+    xs = model.xs[isort, :]
+    xbest = xs[0, :]
+
+    source = problem.get_source(xbest)
+
+    ms, ns, results = problem.evaluate(xbest, result_mode='full')
+
+    figures = []
+
+    def decorateAxes(ax, title):
+        ax.set_title(title)
+        ax.set_xlabel('[km]')
+        scale_axes(ax, 1./km)
+        ax.set_aspect('equal')
+
+    def drawRectangularOutline(ax):
+        source.regularize()
+        fn, fe = source.outline(cs='xy').T
+        offset_n, offset_e = latlon_to_ne_numpy(
+            sat_target.lats[0], sat_target.lons[0],
+            source.lat, source.lon)
+        fn += offset_n
+        fe += offset_e
+        ax.plot(offset_e, offset_n, marker='o')
+        ax.plot(fe, fn, marker='o')
+        # ax.fill(fe, fn, color=(0.5, 0.5, 0.5), alpha=0.5)
+        # ax.plot(fe[:2], fn[:2], linewidth=2., color='black', alpha=0.5)
+
+    def mapDisplacementGrid(displacements, scene):
+        qt = scene.quadtree
+        array = num.empty_like(scene.displacement)
+        array.fill(num.nan)
+        for syn_v, l in zip(displacements, qt.leaves):
+            array[l._slice_rows, l._slice_cols] = syn_v
+
+        array[scene.displacement_mask] = num.nan
+        return array
+
+    def drawTiles(ax, scene):
+        rect = scene.quadtree.getMPLRectangles()
+        for r in rect:
+            r.set_edgecolor((.4, .4, .4))
+            r.set_linewidth(.5)
+            r.set_facecolor('none')
+        map(ax.add_artist, rect)
+
+        ax.scatter(scene.quadtree.leaf_coordinates[:, 0],
+                   scene.quadtree.leaf_coordinates[:, 1],
+                   s=.25, c='black', alpha=.1)
+
+    for sat_target, result in zip(problem.satellite_targets, results):
+        fig = plt.figure()
+        fig.set_size_inches(16, 4)
+        axes = []
+        gs = gridspec.GridSpec(1, 3,
+                               hspace=.0001, left=.06, bottom=.1,
+                               right=.9)
+        axes.append(plt.subplot(gs[0, 0]))
+        axes.append(plt.subplot(gs[0, 1]))
+        axes.append(plt.subplot(gs[0, 2]))
+
+        scene = ds.get_kite_scene(sat_target.scene_id)
+
+        stat_obs = result.statics_obs
+        cmw = cm.ScalarMappable(cmap='coolwarm')