Commit aacba13a authored by Sebastian Heimann's avatar Sebastian Heimann
Browse files

corrected harvesting, improved plots, hudson plot added

parent 0922e8a0
...@@ -20,22 +20,22 @@ as_km = dict(scale_factor=km, scale_unit='km') ...@@ -20,22 +20,22 @@ as_km = dict(scale_factor=km, scale_unit='km')
class CMTProblem(core.Problem): class CMTProblem(core.Problem):
parameters = [ parameters = [
core.Parameter('time', 's'), core.Parameter('time', 's', label='Time'),
core.Parameter('north_shift', 'm', **as_km), core.Parameter('north_shift', 'm', label='Northing', **as_km),
core.Parameter('east_shift', 'm', **as_km), core.Parameter('east_shift', 'm', label='Easting', **as_km),
core.Parameter('depth', 'm', **as_km), core.Parameter('depth', 'm', label='Depth', **as_km),
core.Parameter('magnitude'), core.Parameter('magnitude', label='Magnitude'),
core.Parameter('rmnn'), core.Parameter('rmnn', label='$m_{nn} / M_0$'),
core.Parameter('rmee'), core.Parameter('rmee', label='$m_{ee} / M_0$'),
core.Parameter('rmdd'), core.Parameter('rmdd', label='$m_{dd} / M_0$'),
core.Parameter('rmne'), core.Parameter('rmne', label='$m_{ne} / M_0$'),
core.Parameter('rmnd'), core.Parameter('rmnd', label='$m_{nd} / M_0$'),
core.Parameter('rmed'), core.Parameter('rmed', label='$m_{ed} / M_0$'),
core.Parameter('duration')] core.Parameter('duration', 's', label='Duration')]
dependants = [ dependants = [
core.Parameter('rel_moment_iso'), core.Parameter('rel_moment_iso', label='$M_{0}^{ISO}/M_{0}$'),
core.Parameter('rel_moment_clvd')] core.Parameter('rel_moment_clvd', label='$M_{0}^{CLVD}/M_{0}$')]
base_source = gf.Source.T() base_source = gf.Source.T()
targets = List.T(core.MisfitTarget.T()) targets = List.T(core.MisfitTarget.T())
...@@ -43,7 +43,7 @@ class CMTProblem(core.Problem): ...@@ -43,7 +43,7 @@ class CMTProblem(core.Problem):
ranges = Dict.T(String.T(), gf.Range.T()) ranges = Dict.T(String.T(), gf.Range.T())
distance_min = Float.T(default=0.0) distance_min = Float.T(default=0.0)
nbootstrap = Int.T(default=10) nbootstrap = Int.T(default=10)
mt_type = StringChoice.T(choices=['full', 'deviatoric']) mt_type = StringChoice.T(default='full', choices=['full', 'deviatoric'])
def unpack(self, x): def unpack(self, x):
d = self.parameter_dict(x) d = self.parameter_dict(x)
...@@ -114,6 +114,15 @@ class CMTProblem(core.Problem): ...@@ -114,6 +114,15 @@ class CMTProblem(core.Problem):
return x return x
def random_uniform(self, xbounds):
x = num.zeros(self.nparameters)
for i in [0, 1, 2, 3, 4, 11]:
x[i] = num.random.uniform(xbounds[i, 0], xbounds[i, 1])
x[5:11] = mtm.random_m6()
return x.tolist()
def preconstrain(self, x): def preconstrain(self, x):
d = self.parameter_dict(x) d = self.parameter_dict(x)
......
...@@ -33,6 +33,7 @@ class Parameter(Object): ...@@ -33,6 +33,7 @@ class Parameter(Object):
unit = String.T(optional=True) unit = String.T(optional=True)
scale_factor = Float.T(default=1., optional=True) scale_factor = Float.T(default=1., optional=True)
scale_unit = String.T(optional=True) scale_unit = String.T(optional=True)
label = String.T(optional=True)
def __init__(self, *args, **kwargs): def __init__(self, *args, **kwargs):
if len(args) >= 1: if len(args) >= 1:
...@@ -42,15 +43,35 @@ class Parameter(Object): ...@@ -42,15 +43,35 @@ class Parameter(Object):
Object.__init__(self, **kwargs) Object.__init__(self, **kwargs)
def get_label(self): def get_label(self, with_unit=True):
l = [self.name] l = [self.label or self.name]
if self.scale_unit is not None: if with_unit:
l.append('[%s]' % self.scale_unit) unit = self.get_unit_label()
elif self.unit is not None: if unit:
l.append('[%s]' % self.unit) l.append('[%s]' % unit)
return ' '.join(l) return ' '.join(l)
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): def scaled(self, x):
if isinstance(x, tuple): if isinstance(x, tuple):
return tuple(v/self.scale_factor for v in x) return tuple(v/self.scale_factor for v in x)
...@@ -72,7 +93,7 @@ class Problem(Object): ...@@ -72,7 +93,7 @@ class Problem(Object):
name = String.T() name = String.T()
parameters = List.T(Parameter.T()) parameters = List.T(Parameter.T())
dependants = List.T(Parameter.T()) dependants = List.T(Parameter.T())
apply_balancing_weights = Bool.T() apply_balancing_weights = Bool.T(default=True)
def __init__(self, **kwargs): def __init__(self, **kwargs):
Object.__init__(self, **kwargs) Object.__init__(self, **kwargs)
...@@ -822,6 +843,7 @@ def solve(problem, ...@@ -822,6 +843,7 @@ def solve(problem,
xhist = num.zeros((niter, npar)) xhist = num.zeros((niter, npar))
isbad_mask = None isbad_mask = None
accept_sum = num.zeros(1 + problem.nbootstrap, dtype=num.int) accept_sum = num.zeros(1 + problem.nbootstrap, dtype=num.int)
accept_hist = num.zeros(niter, dtype=num.int)
while iiter < niter: while iiter < niter:
...@@ -843,10 +865,7 @@ def solve(problem, ...@@ -843,10 +865,7 @@ def solve(problem,
ntries_preconstrain += 1 ntries_preconstrain += 1
if mbx is None or iiter < niter_inject + niter_uniform: if mbx is None or iiter < niter_inject + niter_uniform:
x = [] x = problem.random_uniform(xbounds)
for i in xrange(npar):
v = num.random.uniform(xbounds[i, 0], xbounds[i, 1])
x.append(v)
else: else:
# jchoice = num.random.randint(0, 1 + problem.nbootstrap) # jchoice = num.random.randint(0, 1 + problem.nbootstrap)
jchoice = num.argmin(accept_sum) jchoice = num.argmin(accept_sum)
...@@ -932,6 +951,7 @@ def solve(problem, ...@@ -932,6 +951,7 @@ def solve(problem,
accept = num.ones(1 + problem.nbootstrap, dtype=num.int) accept = num.ones(1 + problem.nbootstrap, dtype=num.int)
accept_sum += accept accept_sum += accept
accept_hist[iiter] = num.sum(accept)
lines = [] lines = []
if 'state' in status: if 'state' in status:
...@@ -1105,12 +1125,11 @@ def harvest(rundir, problem=None, nbest=10, force=False, weed=0): ...@@ -1105,12 +1125,11 @@ def harvest(rundir, problem=None, nbest=10, force=False, weed=0):
if gms[ibest] > mean_gm_best + std_gm_best: if gms[ibest] > mean_gm_best + std_gm_best:
ibad.add(ibootstrap) ibad.add(ibootstrap)
print ibad
ibests_list = [ ibests_list = [
ibests_ for (ibootstrap, ibests_) in enumerate(ibests_list) ibests_ for (ibootstrap, ibests_) in enumerate(ibests_list)
if ibootstrap not in ibad] if ibootstrap not in ibad]
ibests = num.vstack(ibests_list) ibests = num.concatenate(ibests_list)
if weed == 2: if weed == 2:
ibests = ibests[gms[ibests] < mean_gm_best] ibests = ibests[gms[ibests] < mean_gm_best]
...@@ -1143,7 +1162,7 @@ def go(config, force=False, nparallel=1, status=('state',)): ...@@ -1143,7 +1162,7 @@ def go(config, force=False, nparallel=1, status=('state',)):
g_state[id(g_data)] = g_data g_state[id(g_data)] = g_data
for x in parimap.parimap( for x in parimap.parimap(
process_event, process_event,
xrange(nevents), xrange(nevents),
[id(g_data)] * nevents, [id(g_data)] * nevents,
nprocs=nparallel): nprocs=nparallel):
...@@ -1171,6 +1190,12 @@ def process_event(ievent, g_data_id): ...@@ -1171,6 +1190,12 @@ def process_event(ievent, g_data_id):
problem = config.get_problem(event) problem = config.get_problem(event)
# FIXME
synt = ds.synthetic_test
if synt and synt.inject_solution:
problem.base_source = problem.unpack(synt.get_x())
check(problem) check(problem)
rundir = config.rundir_template % dict( rundir = config.rundir_template % dict(
......
import math import math
import random import random
import logging
import os.path as op import os.path as op
import numpy as num import numpy as num
from scipy import signal from scipy import signal
from pyrocko import automap, moment_tensor as mtm, beachball, guts, trace, util from pyrocko import automap, moment_tensor as mtm, beachball, guts, trace, util
from pyrocko import hudson
from grond import core from grond import core
from matplotlib import pyplot as plt from matplotlib import pyplot as plt
from matplotlib import cm, patches from matplotlib import cm, patches
from pyrocko.cake_plot import mpl_init, labelspace, colors, \ from pyrocko.cake_plot import mpl_init, labelspace, colors, \
str_to_mpl_color as scolor, light str_to_mpl_color as scolor, light
logger = logging.getLogger('grond.plot')
km = 1000. km = 1000.
...@@ -24,6 +28,21 @@ def nextpow2(i): ...@@ -24,6 +28,21 @@ def nextpow2(i):
return 2**int(math.ceil(math.log(i)/math.log(2.))) return 2**int(math.ceil(math.log(i)/math.log(2.)))
def get_mean_source(problem, xs):
x_mean = num.mean(xs, axis=0)
source = problem.unpack(x_mean)
return source
def get_best_source(problem, xs, misfits):
gms = problem.global_misfits(misfits)
print gms.min(), gms.mean()
ibest = num.argmin(gms)
x_best = xs[ibest, :]
source = problem.unpack(x_best)
return source
def fixlim(lo, hi): def fixlim(lo, hi):
if lo == hi: if lo == hi:
return lo - 1.0, hi + 1.0 return lo - 1.0, hi + 1.0
...@@ -311,7 +330,16 @@ def draw_sequence_figures(model, plt, misfit_cutoff=None): ...@@ -311,7 +330,16 @@ def draw_sequence_figures(model, plt, misfit_cutoff=None):
def draw_jointpar_figures( def draw_jointpar_figures(
model, plt, misfit_cutoff=None, ibootstrap=None, color=None): model, plt, misfit_cutoff=None, ibootstrap=None, color=None,
exclude=None):
color = 'magnitude'
#exclude = ['duration']
neach = 6
figsize = (8, 8)
#cmap = cm.YlOrRd
#cmap = cm.jet
cmap = cm.coolwarm
problem = model.problem problem = model.problem
if not problem: if not problem:
...@@ -340,8 +368,6 @@ def draw_jointpar_figures( ...@@ -340,8 +368,6 @@ def draw_jointpar_figures(
nmodels = xs.shape[0] nmodels = xs.shape[0]
color = 'misfit'
if color == 'dist': if color == 'dist':
mx = num.mean(xs, axis=0) mx = num.mean(xs, axis=0)
cov = num.cov(xs.T) cov = num.cov(xs.T)
...@@ -353,38 +379,47 @@ def draw_jointpar_figures( ...@@ -353,38 +379,47 @@ def draw_jointpar_figures(
color = iorder color = iorder
elif color in problem.parameter_names: elif color in problem.parameter_names:
print 'xx'
ind = problem.name_to_index(color) ind = problem.name_to_index(color)
color = ordersort(problem.extract(xs, ind)) color = ordersort(problem.extract(xs, ind))
# npar = problem.nparameters smap = {}
ncomb = problem.ncombined iselected = 0
for ipar in xrange(problem.ncombined):
par = problem.combined[ipar]
if exclude and par.name in exclude:
continue
neach = 8 smap[iselected] = ipar
cmap = cm.YlOrRd iselected += 1
cmap = cm.jet
#cmap = cm.coolwarm nselected = iselected
nfig = (ncomb-1) / neach + 1 if nselected == 0:
return
nfig = (nselected-2) / neach + 1
figs = [] figs = []
for ifig in xrange(nfig): for ifig in xrange(nfig):
figs_row = [] figs_row = []
for jfig in xrange(nfig): for jfig in xrange(nfig):
if ifig >= jfig: if ifig >= jfig:
figs_row.append(plt.figure()) figs_row.append(plt.figure(figsize=figsize))
else: else:
figs_row.append(None) figs_row.append(None)
figs.append(figs_row) figs.append(figs_row)
for ipar in xrange(ncomb): for iselected in xrange(nselected):
ipar = smap[iselected]
ypar = problem.combined[ipar] ypar = problem.combined[ipar]
for jpar in xrange(ipar): for jselected in xrange(iselected):
jpar = smap[jselected]
xpar = problem.combined[jpar] xpar = problem.combined[jpar]
ixg = (ipar - 1) ixg = (iselected - 1)
iyg = jpar iyg = jselected
ix = ixg % neach ix = ixg % neach
iy = iyg % neach iy = iyg % neach
...@@ -398,15 +433,48 @@ def draw_jointpar_figures( ...@@ -398,15 +433,48 @@ def draw_jointpar_figures(
axes = fig.add_subplot(*aind) axes = fig.add_subplot(*aind)
axes.set_xlim(*fixlim(*xpar.scaled(bounds[jpar]))) axes.axvline(0., color=scolor('aluminium3'), lw=0.5)
axes.set_ylim(*fixlim(*ypar.scaled(bounds[ipar]))) axes.axhline(0., color=scolor('aluminium3'), lw=0.5)
for spine in axes.spines.values():
spine.set_edgecolor(scolor('aluminium5'))
spine.set_linewidth(0.5)
xmin, xmax = fixlim(*xpar.scaled(bounds[jpar]))
ymin, ymax = fixlim(*ypar.scaled(bounds[ipar]))
if ix == 0 or jselected + 1 == iselected:
for (xpos, xoff, x) in [(0.0, 10., xmin), (1.0, -10., xmax)]:
axes.annotate(
'%.2g%s' % (x, xpar.get_unit_suffix()),
xy=(xpos, 1.05),
xycoords='axes fraction',
xytext=(xoff, 5.),
textcoords='offset points',
verticalalignment='bottom',
horizontalalignment='left',
rotation=45.)
if iy == neach - 1 or jselected + 1 == iselected:
for (ypos, yoff, y) in [(0., 10., ymin), (1.0, -10., ymax)]:
axes.annotate(
'%.2g%s' % (y, ypar.get_unit_suffix()),
xy=(1.0, ypos),
xycoords='axes fraction',
xytext=(5., yoff),
textcoords='offset points',
verticalalignment='bottom',
horizontalalignment='left',
rotation=45.)
axes.set_xlim(xmin, xmax)
axes.set_ylim(ymin, ymax)
axes.get_xaxis().set_ticks([]) axes.get_xaxis().set_ticks([])
axes.get_yaxis().set_ticks([]) axes.get_yaxis().set_ticks([])
if ipar == ncomb - 1 or ix == neach - 1: if iselected == nselected - 1 or ix == neach - 1:
axes.annotate( axes.annotate(
xpar.get_label(), xpar.get_label(with_unit=False),
xy=(0.5, -0.05), xy=(0.5, -0.05),
xycoords='axes fraction', xycoords='axes fraction',
verticalalignment='top', verticalalignment='top',
...@@ -415,21 +483,21 @@ def draw_jointpar_figures( ...@@ -415,21 +483,21 @@ def draw_jointpar_figures(
if iy == 0: if iy == 0:
axes.annotate( axes.annotate(
ypar.get_label(), ypar.get_label(with_unit=False),
xy=(-0.05, 0.5), xy=(-0.05, 0.5),
xycoords='axes fraction', xycoords='axes fraction',
verticalalignment='center', verticalalignment='top',
horizontalalignment='right') horizontalalignment='right',
rotation=45.)
fx = problem.extract(xs, jpar) fx = problem.extract(xs, jpar)
fy = problem.extract(xs, ipar) fy = problem.extract(xs, ipar)
axes.scatter( axes.scatter(
xpar.scaled(fx), xpar.scaled(fx),
ypar.scaled(fy), ypar.scaled(fy),
c=color, c=color,
s=3, alpha=0.5, lw=0, cmap=cmap) s=3, alpha=0.5, cmap=cmap, edgecolors='none')
cov = num.cov((xpar.scaled(fx), ypar.scaled(fy))) cov = num.cov((xpar.scaled(fx), ypar.scaled(fy)))
evals, evecs = eigh_sorted(cov) evals, evecs = eigh_sorted(cov)
...@@ -445,8 +513,17 @@ def draw_jointpar_figures( ...@@ -445,8 +513,17 @@ def draw_jointpar_figures(
fx = problem.extract(xref, jpar) fx = problem.extract(xref, jpar)
fy = problem.extract(xref, ipar) fy = problem.extract(xref, ipar)
axes.axvline(xpar.scaled(fx), color='black', alpha=0.3)
axes.axhline(ypar.scaled(fy), color='black', alpha=0.3) ref_color = scolor('aluminium6')
ref_color_light = 'none'
axes.plot(
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))
def draw_solution_figure( def draw_solution_figure(
...@@ -640,7 +717,6 @@ def draw_bootstrap_figure(model, plt): ...@@ -640,7 +717,6 @@ def draw_bootstrap_figure(model, plt):
bms_softclip = num.where(bms > 1.0, 0.1 * num.log10(bms) + 1.0, bms) bms_softclip = num.where(bms > 1.0, 0.1 * num.log10(bms) + 1.0, bms)
axes.plot(imodels, bms_softclip[isort_bms], color='red', alpha=0.2) axes.plot(imodels, bms_softclip[isort_bms], color='red', alpha=0.2)
isort = num.argsort(gms)[::-1] isort = num.argsort(gms)[::-1]
iorder = num.empty(isort.size) iorder = num.empty(isort.size)
iorder[isort] = imodels iorder[isort] = imodels
...@@ -654,16 +730,12 @@ def draw_bootstrap_figure(model, plt): ...@@ -654,16 +730,12 @@ def draw_bootstrap_figure(model, plt):
axes.axhline(m, color='black') axes.axhline(m, color='black')
axes.axhline(m-s, color='black', alpha=0.5) axes.axhline(m-s, color='black', alpha=0.5)
#axes.plot(imodels[ibests], gms_softclip[isort[ibests]], 'x', color='black')
#axes.plot(imodels[ibests], gms_softclip[isort][ibests], '+', color='black')
#iii = isort[-1]
#axes.plot(imodels[iii], gms_softclip[isort][iii], 'o')
axes.plot(imodels, gms_softclip[isort], color='black') axes.plot(imodels, gms_softclip[isort], color='black')
axes.set_xlim(imodels[0], imodels[-1]) axes.set_xlim(imodels[0], imodels[-1])
axes.set_xlabel('Tested model, sorted descending by global misfit value') axes.set_xlabel('Tested model, sorted descending by global misfit value')
def gather(l, key, sort=None, filter=None): def gather(l, key, sort=None, filter=None):
d = {} d = {}
for x in l: for x in l:
...@@ -859,7 +931,7 @@ def draw_fits_figures(ds, model, plt): ...@@ -859,7 +931,7 @@ def draw_fits_figures(ds, model, plt):
ixx = ix/nxmax ixx = ix/nxmax
iyy = iy/nymax iyy = iy/nymax
if (iyy, ixx) not in figures: if (iyy, ixx) not in figures:
figures[iyy, ixx] = plt.figure() figures[iyy, ixx] = plt.figure(figsize=(16, 9))
fig = figures[iyy, ixx] fig = figures[iyy, ixx]
...@@ -992,10 +1064,93 @@ def draw_fits_figures(ds, model, plt): ...@@ -992,10 +1064,93 @@ def draw_fits_figures(ds, model, plt):
title += ' (%i/%i, %i/%i)' % (iyy+1, nyy, ixx+1, nxx) title += ' (%i/%i, %i/%i)' % (iyy+1, nyy, ixx+1, nxx)
fig.suptitle(title, fontsize=fontsize) fig.suptitle(title, fontsize=fontsize)
fig.savefig('fits-%s-%i-%i.pdf' % (
title, iyy, ixx))
plt.show() plt.show()
def draw_hudson_figure(model, plt):
color = 'black'
fontsize = 12.
markersize = fontsize * 1.5
markersize_small = markersize * 0.2
beachballsize = markersize
beachballsize_small = beachballsize * 0.5
width = 7.
figsize = (width, width / (4./3.))