Merge branch 'master' of gitext:heimann/grond

Conflicts:
	src/solvers/highscore.py

apps/grond

@@ -171,15 +171,15 @@ def command_init(args):
             events_path='events.txt',
             waveform_paths=['data'])
 
-        target_configs = [grond.TargetConfig(
-            super_group='time_domain',
-            group='all',
+        target_groups = [grond.WaveformTargetGroup(
+            normalisation_family='time_domain',
+            path='all',
             distance_min=10*km,
             distance_max=1000*km,
             channels=['Z', 'R', 'T'],
             interpolation='multilinear',
             store_id='gf_store',
-            inner_misfit_config=grond.InnerMisfitConfig(
+            misfit_config=grond.WaveformMisfitConfig(
                 fmin=0.01,
                 fmax=0.1))]
 
@@ -213,14 +213,21 @@ def command_init(args):
             kite_scene_paths=['scenes'],
             )
 
-        target_configs = [grond.TargetConfig(
-            super_group='insar_target',
-            group='all',
+        target_groups = [grond.SatelliteTargetGroup(
+            normalisation_family='insar_target',
+            path='all',
             interpolation='multilinear',
             store_id='gf_store',
             kite_scenes=['*all'],
-            inner_misfit_config=grond.InnerSatelliteMisfitConfig(
-                use_weight_focal=False))]
+            misfit_config=grond.SatelliteMisfitConfig(
+                use_weight_focal=False,
+                optimize_orbital_ramp=True,
+                ranges={
+                    'offset': '-0.5 .. 0.5',
+                    'ramp_north': '-1e-4 .. 1e4',
+                    'ramp_east': '-1e-4 .. 1e4'
+                    }
+                ))]
 
         problem_config = grond.RectangularProblemConfig(
             name_template='rect_source',
@@ -242,7 +249,7 @@ def command_init(args):
     config = grond.Config(
         rundir_template=op.join(op.abspath(op.curdir), 'rundir'),
         dataset_config=dataset_config,
-        target_configs=target_configs,
+        target_groups=target_groups,
         problem_config=problem_config,
         engine_config=engine_config)
 
@@ -289,8 +296,9 @@ def command_events(args):
     config_path = args[0]
     config = grond.read_config(config_path)
 
+    print 'Available Events:'
     for event_name in grond.get_event_names(config):
-        print event_name
+        print '* %s' % event_name
 
 
 def command_check(args):
@@ -355,6 +363,18 @@ def command_go(args):
             help='start the Barad-dur plotting server')
 
     parser, options, args = cl_parse('go', args, setup)
+
+    if len(args) == 1:
+        config_path = args[0]
+        config = grond.read_config(config_path)
+
+        print 'Available Events:'
+        for event_name in grond.get_event_names(config):
+            print '* %s' % event_name
+        print '\n'
+
+        help_and_die(parser, 'missing arguments')
+
     if len(args) < 2:
         help_and_die(parser, 'missing arguments')
 

setup.py

@@ -16,7 +16,7 @@ setup(
     author='Sebastian Heimann',
     author_email='sebastian.heimann@gfz-potsdam.de',
     packages=['grond', 'grond.baraddur', 'grond.problems', 'grond.solvers',
-              'grond.analysers', 'grond.listeners'],
+              'grond.analysers', 'grond.listeners', 'grond.targets'],
     scripts=['apps/grond'],
     package_dir={'grond': 'src'},
     package_data={'grond': ['baraddur/templates/*.html',

src/core.py

@@ -14,9 +14,9 @@ from pyrocko import parimap, model, marker as pmarker
 from .dataset import DatasetConfig, NotFound
 from .problems.base import ProblemConfig, Problem
 from .solvers.base import SolverConfig
+from .targets.base import TargetGroup
 from .analysers.base import AnalyserConfig
 from .listeners import TerminalListener
-from .targets import TargetConfig
 from .meta import Path, HasPaths, expand_template, xjoin, GrondError, Notifier
 
 logger = logging.getLogger('grond.core')
@@ -77,7 +77,7 @@ class EngineConfig(HasPaths):
 class Config(HasPaths):
     rundir_template = Path.T()
     dataset_config = DatasetConfig.T()
-    target_configs = List.T(TargetConfig.T())
+    target_groups = List.T(TargetGroup.T())
     problem_config = ProblemConfig.T()
     analyser_config = AnalyserConfig.T(default=AnalyserConfig.D())
     solver_config = SolverConfig.T(default=SolverConfig.D())
@@ -96,9 +96,9 @@ class Config(HasPaths):
         ds = self.get_dataset(event.name)
 
         targets = []
-        for igroup, target_config in enumerate(self.target_configs):
-            targets.extend(target_config.get_targets(
-                ds, event, 'group_%i' % igroup))
+        for igroup, target_group in enumerate(self.target_groups):
+            targets.extend(target_group.get_targets(
+                ds, event, 'target.%i' % igroup))
 
         return targets
 

src/listeners/curses.py

@@ -119,7 +119,7 @@ class CursesListener(object):
     def set_screen(self, scr):
         self.scr = scr
 
-    def set_state(self, state):
+    def state(self, state):
         self.state = state
         self.parameter_pad.update(state)
         self.header_pad.update(state)

src/listeners/terminal.py

 from pyrocko import util
-import progressbar as pbar
 from .base import Listener
 
 
@@ -45,14 +44,18 @@ class TerminalListener(Listener):
         def l(t):
             lines.append(t)
 
+        def fmt(s):
+            return util.gform(s, significant_digits=(self.col_width-1-6)/2)
+
         out_ln = self.row_name +\
             ''.join([self.parameter_fmt] * len(state.parameter_sets))
         col_param_width = max([len(p) for p in state.parameter_names]) + 2
 
         l('Problem name: {s.problem_name}'
-          '\t({s.runtime:s} - remaining {s.runtime_remaining})'
+          '\t({s.runtime:s} - remaining {s.runtime_remaining}'
+          ' @ {s.iter_per_second:.1f} iter/s)'
             .format(s=state))
-        l('Iteration {s.iiter} / {s.niter}\t\t({s.iter_per_second:.1f} iter/s)'
+        l('Iteration {s.iiter} / {s.niter}'
           .format(s=state))
 
         l(out_ln.format(
@@ -61,9 +64,6 @@ class TerminalListener(Listener):
             col_width=self.col_width,
             type='s'))
 
-        def fmt(s):
-            return util.gform(s, significant_digits=(self.col_width-1-6)/2)
-
         for ip, parameter_name in enumerate(state.parameter_names):
             l(out_ln.format(
                 parameter_name,

src/plot.py

@@ -1151,17 +1151,21 @@ def draw_fits_figures(ds, model, plt):
 
         target_to_result[target] = result
 
-        dtrace.meta = dict(super_group=target.super_group, group=target.group)
+        dtrace.meta = dict(
+            normalisation_family=target.normalisation_family,
+            path=target.path)
         dtraces.append(dtrace)
 
         result.processed_syn.meta = dict(
-            super_group=target.super_group, group=target.group)
+            normalisation_family=target.normalisation_family,
+            path=target.path)
 
         all_syn_trs.append(result.processed_syn)
 
         if result.spectrum_syn:
             result.spectrum_syn.meta = dict(
-                super_group=target.super_group, group=target.group)
+                normalisation_family=target.normalisation_family,
+                path=target.path)
 
             all_syn_specs.append(result.spectrum_syn)
 
@@ -1170,7 +1174,7 @@ def draw_fits_figures(ds, model, plt):
         return []
 
     def skey(tr):
-        return tr.meta['super_group'], tr.meta['group']
+        return tr.meta['normalisation_family'], tr.meta['path']
 
     trace_minmaxs = trace.minmax(all_syn_trs, skey)
 
@@ -1185,7 +1189,7 @@ def draw_fits_figures(ds, model, plt):
 
     cg_to_targets = gather(
         problem.waveform_targets,
-        lambda t: (t.super_group, t.group, t.codes[3]),
+        lambda t: (t.normalisation_family, t.path, t.codes[3]),
         filter=lambda t: t in target_to_result)
 
     cgs = sorted(cg_to_targets.keys())
@@ -1283,7 +1287,8 @@ def draw_fits_figures(ds, model, plt):
 
                 target = frame_to_target[iy, ix]
 
-                amin, amax = trace_minmaxs[target.super_group, target.group]
+                amin, amax = trace_minmaxs[
+                    target.normalisation_family, target.path]
                 absmax = max(abs(amin), abs(amax))
 
                 ny_this = nymax  # min(ny, nymax)
@@ -1343,7 +1348,8 @@ def draw_fits_figures(ds, model, plt):
 
                 elif target.misfit_config.domain == 'frequency_domain':
 
-                    asmax = amp_spec_maxs[target.super_group, target.group]
+                    asmax = amp_spec_maxs[
+                        target.normalisation_family, target.path]
                     fmin, fmax = \
                         target.misfit_config.get_full_frequency_range()
 
@@ -1848,7 +1854,8 @@ class SolverPlot(object):
 
         self.bcolors = colors.hsv_to_rgb(hsv[num.newaxis, :, :])[0, :, :]
 
-        bounds = self.problem.get_parameter_bounds() + self.problem.get_dependant_bounds()
+        bounds = self.problem.get_parameter_bounds()\
+            + self.problem.get_dependant_bounds()
 
         self.xlim = fixlim(*xpar.scaled(bounds[ixpar]))
         self.ylim = fixlim(*ypar.scaled(bounds[iypar]))

src/problems/base.py

@@ -7,7 +7,7 @@ from pyrocko import gf, util, guts
 from pyrocko.guts import Object, String, Bool, List, Dict, Int
 
 from ..meta import ADict, Parameter, GrondError
-from ..targets import MisfitTarget, MisfitSatelliteTarget
+from ..targets import WaveformMisfitTarget, SatelliteMisfitTarget
 
 
 guts_prefix = 'grond'
@@ -152,12 +152,12 @@ class Problem(Object):
     @property
     def satellite_targets(self):
         return [t for t in self.targets
-                if isinstance(t, MisfitSatelliteTarget)]
+                if isinstance(t, SatelliteMisfitTarget)]
 
     @property
     def waveform_targets(self):
         return [t for t in self.targets
-                if isinstance(t, MisfitTarget)]
+                if isinstance(t, WaveformMisfitTarget)]
 
     @property
     def has_statics(self):
@@ -338,17 +338,14 @@ class Problem(Object):
         return gcms
 
     def make_group_mask(self):
-        super_group_names = set()
-        groups = num.zeros(len(self.targets), dtype=num.int)
-        ngroups = 0
-        for itarget, target in enumerate(self.targets):
-            if target.super_group not in super_group_names:
-                super_group_names.add(target.super_group)
-                ngroups += 1
+        family_names = set()
+        families = num.zeros(len(self.targets), dtype=num.int)
 
-            groups[itarget] = ngroups - 1
+        for itarget, target in enumerate(self.targets):
+            family_names.add(target.normalisation_family)
+            families[itarget] = len(family_names) - 1
 
-        return groups, ngroups
+        return families, len(family_names)
 
     def get_group_mask(self):
         if self._group_mask is None:

src/problems/cmt.py

 import numpy as num
 import math
+import logging
 
 from pyrocko import gf, util, moment_tensor as mtm
 from pyrocko.guts import String, Float, Dict, Int, StringChoice
 
-from .base import Problem, ProblemConfig, logger
+from .base import Problem, ProblemConfig
 from ..meta import Forbidden, expand_template, Parameter
 
 guts_prefix = 'grond'
-logger = logger.getChild('cmt')
+logger = logging.getLogger('grond.problems').getChild('cmt')
 km = 1e3
 as_km = dict(scale_factor=km, scale_unit='km')
 

src/problems/double_dc.py

 import numpy as num
+import logging
 
 from pyrocko import gf, util
 from pyrocko.guts import String, Float, Dict, Int
 
-from .base import Problem, ProblemConfig, logger
+from .base import Problem, ProblemConfig
 from ..meta import Forbidden, expand_template, Parameter
 
 
 guts_prefix = 'grond'
-logger = logger.getChild('grond.doubleDC')
+logger = logging.getLogger('grond.problems').getChild('double_dc')
 km = 1e3
 as_km = dict(scale_factor=km, scale_unit='km')
 

src/problems/rectangular.py

 import numpy as num
+import logging
 
 from pyrocko import gf
 from pyrocko.guts import String, Bool, Float, Dict, Int
 
-from .base import Problem, ProblemConfig, logger
+from .base import Problem, ProblemConfig
 from ..meta import expand_template, Parameter
 
 guts_prefix = 'grond'
-logger = logger.getChild('rectangular')
+logger = logging.getLogger('grond.problems').getChild('rectangular')
 km = 1e3
 as_km = dict(scale_factor=km, scale_unit='km')
 

src/solvers/highscore.py

+import math
+import logging
+
+import numpy as num
+
+from collections import OrderedDict
+
+from pyrocko.guts import StringChoice, Int, Float, Bool
+
+from ..meta import GrondError, Forbidden
+from .base import Solver, SolverConfig
+
+guts_prefix = 'grond'
+
+logger = logging.getLogger('grond.solver.highscore')
+
+
+def excentricity_compensated_probabilities(xs, sbx, factor):
+    inonflat = num.where(sbx != 0.0)[0]
+    scale = num.zeros_like(sbx)
+    scale[inonflat] = 1.0 / (sbx[inonflat] * (factor if factor != 0. else 1.0))
+    distances_sqr_all = num.sum(
+        ((xs[num.newaxis, :, :] - xs[:, num.newaxis, :]) *
+         scale[num.newaxis, num.newaxis, :])**2, axis=2)
+    probabilities = 1.0 / num.sum(distances_sqr_all < 1.0, axis=1)
+    print num.sort(num.sum(distances_sqr_all < 1.0, axis=1))
+    probabilities /= num.sum(probabilities)
+    return probabilities
+
+
+def excentricity_compensated_choice(xs, sbx, factor):
+    probabilities = excentricity_compensated_probabilities(
+        xs, sbx, factor)
+    r = num.random.random()
+    ichoice = num.searchsorted(num.cumsum(probabilities), r)
+    ichoice = min(ichoice, xs.shape[0]-1)
+    return ichoice
+
+
+def select_most_excentric(xcandidates, xs, sbx, factor):
+    inonflat = num.where(sbx != 0.0)[0]
+    scale = num.zeros_like(sbx)
+    scale[inonflat] = 1.0 / (sbx[inonflat] * (factor if factor != 0. else 1.0))
+    distances_sqr_all = num.sum(
+        ((xcandidates[num.newaxis, :, :] - xs[:, num.newaxis, :]) *
+         scale[num.newaxis, num.newaxis, :])**2, axis=2)
+    # print num.sort(num.sum(distances_sqr_all < 1.0, axis=0))
+    ichoice = num.argmin(num.sum(distances_sqr_all < 1.0, axis=0))
+    return xcandidates[ichoice]
+
+
+def local_std(xs):
+    ssbx = num.sort(xs, axis=0)
+    dssbx = num.diff(ssbx, axis=0)
+    mdssbx = num.median(dssbx, axis=0)
+    return mdssbx * dssbx.shape[0] / 2.6
+
+
+def solve(problem,
+          rundir=None,
+          niter_uniform=1000,
+          niter_transition=1000,
+          niter_explorative=10000,
+          niter_non_explorative=0,
+          scatter_scale_transition=2.0,
+          scatter_scale=1.0,
+          chain_length_factor=8.0,
+          sampler_distribution='multivariate_normal',
+          standard_deviation_estimator='median_density_single_chain',
+          compensate_excentricity=True,
+          xs_inject=None,
+          status=(),
+          plot=None,
+          notifier=None,
+          state=None):
+
+    xbounds = num.array(problem.get_parameter_bounds(), dtype=num.float)
+    npar = problem.nparameters
+
+    nlinks_cap = int(round(chain_length_factor * npar + 1))
+    chains_m = num.zeros((1 + problem.nbootstrap, nlinks_cap), num.float)
+    chains_i = num.zeros((1 + problem.nbootstrap, nlinks_cap), num.int)
+    nlinks = 0
+    mbx = None
+
+    if xs_inject is not None and xs_inject.size != 0:
+        niter_inject = xs_inject.shape[0]
+    else:
+        niter_inject = 0
+
+    niter = niter_inject + niter_uniform + niter_transition + \
+        niter_explorative + niter_non_explorative
+
+    iiter = 0
+    sbx = None
+    mxs = None
+    covs = None
+    local_sxs = None
+    xhist = num.zeros((niter, npar))
+    isbad_mask = None
+    accept_sum = num.zeros(1 + problem.nbootstrap, dtype=num.int)
+    accept_hist = num.zeros(niter, dtype=num.int)
+    pnames = problem.parameter_names
+
+    state.problem_name = problem.name
+    state.parameter_names = problem.parameter_names + ['Misfit']
+    npar = len(state.parameter_names)
+
+    state.parameter_sets = OrderedDict()
+
+    state.parameter_sets['BS mean'] = num.zeros(npar)
+    state.parameter_sets['BS std'] = num.zeros(npar)
+    state.parameter_sets['Global mean'] = num.zeros(npar)
+    state.parameter_sets['Global std'] = num.zeros(npar)
+    state.parameter_sets['Global best'] = num.zeros(npar)
+
+    for par in state.parameter_sets.values():
+        par.fill(num.nan)
+
+    state.niter = niter
+
+    if plot:
+        plot.start(problem)
+
+    while iiter < niter:
+        ibootstrap_choice = None
+        ichoice = None
+
+        if iiter < niter_inject:
+            phase = 'inject'
+        elif iiter < niter_inject + niter_uniform:
+            phase = 'uniform'
+        elif iiter < niter_inject + niter_uniform + niter_transition:
+            phase = 'transition'
+        elif iiter < niter_inject + niter_uniform + niter_transition + \
+                niter_explorative:
+            phase = 'explorative'
+        else:
+            phase = 'non_explorative'
+
+        factor = 0.0
+        if phase == 'transition':
+            T = float(niter_transition)
+            A = scatter_scale_transition
+            B = scatter_scale
+            tau = T/(math.log(A) - math.log(B))
+            t0 = math.log(A) * T / (math.log(A) - math.log(B))
+            t = float(iiter - niter_uniform - niter_inject)
+            factor = num.exp(-(t-t0) / tau)
+
+        elif phase in ('explorative', 'non_explorative'):
+            factor = scatter_scale
+
+        ntries_preconstrain = 0
+        ntries_sample = 0
+
+        if phase == 'inject':
+            x = xs_inject[iiter, :]
+        else:
+            while True:
+                ntries_preconstrain += 1
+
+                if mbx is None or phase == 'uniform':
+                    x = problem.random_uniform(xbounds)
+                else:
+                    # ibootstrap_choice = num.random.randint(
+                    #     0, 1 + problem.nbootstrap)
+                    ibootstrap_choice = num.argmin(accept_sum)
+
+                    if phase in ('transition', 'explorative'):
+
+                        if compensate_excentricity:
+                            ichoice = excentricity_compensated_choice(
+                                xhist[chains_i[ibootstrap_choice, :], :],
+                                local_sxs[ibootstrap_choice], 2.)
+
+                            xchoice = xhist[
+                                chains_i[ibootstrap_choice, ichoice], :]
+
+                        else:
+                            ichoice = num.random.randint(0, nlinks)
+                            xchoice = xhist[
+                                chains_i[ibootstrap_choice, ichoice], :]
+                    else:
+                        xchoice = mxs[ibootstrap_choice]
+
+                    if sampler_distribution == 'multivariate_normal':
+                        ntries_sample = 0
+
+                        ntry = 0
+                        ok_mask_sum = num.zeros(npar, dtype=num.int)
+                        while True:
+                            ntries_sample += 1
+                            vs = num.random.multivariate_normal(
+                                xchoice, factor**2 * covs[ibootstrap_choice])
+
+                            ok_mask = num.logical_and(
+                                xbounds[:, 0] <= vs, vs <= xbounds[:, 1])
+
+                            if num.all(ok_mask):
+                                break
+
+                            ok_mask_sum += ok_mask
+
+                            if ntry > 1000:
+                                raise GrondError(
+                                    'failed to produce a suitable candidate '
+                                    'sample from multivariate normal '
+                                    'distribution, (%s)' %
+                                    ', '.join('%s:%i' % xx for xx in
+                                              zip(pnames, ok_mask_sum)))
+
+                            ntry += 1
+
+                        x = vs.tolist()
+
+                    if sampler_distribution == 'normal':
+                        ncandidates = 1
+                        xcandidates = num.zeros((ncandidates, npar))
+                        for icandidate in xrange(ncandidates):
+                            for ipar in xrange(npar):
+                                ntry = 0
+                                while True:
+                                    if local_sxs[ibootstrap_choice][ipar] > 0.:
+                                        v = num.random.normal(
+                                            xchoice[ipar],
+                                            factor*local_sxs[
+                                                ibootstrap_choice][ipar])
+                                    else:
+                                        v = xchoice[ipar]
+
+                                    if xbounds[ipar, 0] <= v and \
+                                            v <= xbounds[ipar, 1]:
+
+                                        break
+
+                                    if ntry > 1000:
+                                        raise GrondError(
+                                            'failed to produce a suitable '
+                                            'candidate sample from normal '
+                                            'distribution')
+
+                                    ntry += 1
+
+                                xcandidates[icandidate, ipar] = v
+
+                        x = select_most_excentric(
+                            xcandidates,
+                            xhist[chains_i[ibootstrap_choice, :], :],
+                            local_sxs[ibootstrap_choice],
+                            factor)
+
+                try:
+                    x = problem.preconstrain(x)
+                    break
+
+                except Forbidden:
+                    pass
+
+        ibase = None
+        if ichoice is not None and ibootstrap_choice is not None:
+            ibase = chains_i[ibootstrap_choice, ichoice]