Commit df3132b0 authored by Andreas Steinberg's avatar Andreas Steinberg

added stationarity test

parent 0b2f53ec
......@@ -130,6 +130,9 @@ def help_and_die(parser, message):
def command_init(args):
rundir_template = grond.rundir(
rundir_template='runs/aquila1bounds.run')
dataset_config = grond.DatasetConfig(
stations_path='stations.txt',
events_path='events.txt',
......@@ -143,7 +146,9 @@ def command_init(args):
store_id='global_2s',
inner_misfit_config=grond.InnerMisfitConfig(
fmin=0.01,
fmax=0.1))]
fmax=0.1,
tmin='P-5',
tmax='P+5'))]
s2 = math.sqrt(2.0)
......@@ -170,6 +175,9 @@ def command_init(args):
dataset_config=dataset_config,
target_configs=target_configs,
problem_config=problem_config)
print config
......
This diff is collapsed.
grond @ 2af71c5f
Subproject commit 2af71c5f5e005f17af624fa50535aa467d2f27c7
......@@ -18,32 +18,61 @@ as_km = dict(scale_factor=km, scale_unit='km')
class CMTProblem(core.Problem):
parameters = [
core.Parameter('time', 's', label='Time'),
core.Parameter('north_shift', 'm', label='Northing', **as_km),
core.Parameter('east_shift', 'm', label='Easting', **as_km),
core.Parameter('depth', 'm', label='Depth', **as_km),
core.Parameter('magnitude', label='Magnitude'),
core.Parameter('rmnn', label='$m_{nn} / M_0$'),
core.Parameter('rmee', label='$m_{ee} / M_0$'),
core.Parameter('rmdd', label='$m_{dd} / M_0$'),
core.Parameter('rmne', label='$m_{ne} / M_0$'),
core.Parameter('rmnd', label='$m_{nd} / M_0$'),
core.Parameter('rmed', label='$m_{ed} / M_0$'),
core.Parameter('duration', 's', label='Duration')]
dependants = [
core.Parameter('rel_moment_iso', label='$M_{0}^{ISO}/M_{0}$'),
core.Parameter('rel_moment_clvd', label='$M_{0}^{CLVD}/M_{0}$')]
base_source = gf.Source.T()
targets = List.T(core.MisfitTarget.T())
ranges = Dict.T(String.T(), gf.Range.T())
distance_min = Float.T(default=0.0)
nbootstrap = Int.T(default=10)
mt_type = StringChoice.T(default='full', choices=['full', 'deviatoric'])
sourcetype ="point-source"
#config = core.MisfitTarget.misfit_config
if sourcetype == "rectangular_source":
parameters = [
core.Parameter('time', 's', label='Time'),
core.Parameter('north_shift', 'm', label='Northing', **as_km),
core.Parameter('east_shift', 'm', label='Easting', **as_km),
core.Parameter('depth', 'm', label='Depth', **as_km),
core.Parameter('length'),
core.Parameter('width'),
core.Parameter('strike1'),
core.Parameter('dip1'),
core.Parameter('rake1'),
core.Parameter('magnitude'),
core.Parameter('nucleation_x'),
core.Parameter('nucleation_y')]
# dependants = []
#core.Parameter('rel_moment_iso', label='$M_{0}^{ISO}/M_{0}$'),
# core.Parameter('rel_moment_clvd', label='$M_{0}^{CLVD}/M_{0}$')]
base_source = gf.RectangularSource.T()
targets = List.T(core.MisfitTarget.T())
ranges = Dict.T(String.T(), gf.Range.T())
distance_min = Float.T(default=0.0)
nbootstrap = Int.T(default=10)
mt_type = StringChoice.T(default='full', choices=['full', 'deviatoric'])
else:
parameters = [
core.Parameter('time', 's', label='Time'),
core.Parameter('north_shift', 'm', label='Northing', **as_km),
core.Parameter('east_shift', 'm', label='Easting', **as_km),
core.Parameter('depth', 'm', label='Depth', **as_km),
core.Parameter('magnitude', label='Magnitude'),
core.Parameter('rmnn', label='$m_{nn} / M_0$'),
core.Parameter('rmee', label='$m_{ee} / M_0$'),
core.Parameter('rmdd', label='$m_{dd} / M_0$'),
core.Parameter('rmne', label='$m_{ne} / M_0$'),
core.Parameter('rmnd', label='$m_{nd} / M_0$'),
core.Parameter('rmed', label='$m_{ed} / M_0$'),
core.Parameter('duration', 's', label='Duration')]
dependants = []
core.Parameter('rel_moment_iso', label='$M_{0}^{ISO}/M_{0}$'),
core.Parameter('rel_moment_clvd', label='$M_{0}^{CLVD}/M_{0}$')]
base_source = gf.Source.T()
targets = List.T(core.MisfitTarget.T())
ranges = Dict.T(String.T(), gf.Range.T())
distance_min = Float.T(default=0.0)
nbootstrap = Int.T(default=10)
mt_type = StringChoice.T(default='full', choices=['full', 'deviatoric'])
def unpack(self, x):
d = self.parameter_dict(x)
......
This diff is collapsed.
......@@ -6,6 +6,7 @@ import time
import copy
import shutil
import os.path as op
import random
import numpy as num
......@@ -350,7 +351,7 @@ class MisfitTarget(gf.Target):
return tmin_fit, tmax_fit, tfade
def post_process(self, engine, source, tr_syn):
tr_syn = tr_syn.pyrocko_trace()
nslc = self.codes
......@@ -417,6 +418,20 @@ class MisfitTarget(gf.Target):
cache=True,
backazimuth=backazimuth)
noisesamples = 100
pertubate_noise = True
if pertubate_noise is not None:
for c in range (noisesamples):
l= random.uniform(tmin_obs+5, tmin_fit-20)
extracted_obs = tr_obs.chop(l,l+5, inplace=False)
ltrace= random.uniform(tmin_obs+5, tmin_fit-20)
extracted_obs.shift(ltrace)
tr_obs.add(extracted_obs, interpolate=True)
if tobs_shift != 0.0:
tr_obs = tr_obs.copy()
tr_obs.shift(-tobs_shift)
......@@ -1044,23 +1059,12 @@ def get_mean_x(xs):
return num.mean(xs, axis=0)
def get_mean_x_and_gm(problem, xs, misfits):
gms = problem.global_misfits(misfits)
return num.mean(xs, axis=0), num.mean(gms)
def get_best_x(problem, xs, misfits):
gms = problem.global_misfits(misfits)
ibest = num.argmin(gms)
return xs[ibest, :]
def get_best_x_and_gm(problem, xs, misfits):
gms = problem.global_misfits(misfits)
ibest = num.argmin(gms)
return xs[ibest, :], gms[ibest]
def get_mean_source(problem, xs):
x_mean = get_mean_x(xs)
source = problem.unpack(x_mean)
......@@ -1850,10 +1854,9 @@ def export(what, rundirs, type=None, pnames=None, filename=None):
if shortform:
print >>out, '#', ' '.join('%16s' % x for x in pnames)
def dump(x, gm, indices):
def dump(x, indices):
if type == 'vector':
print >>out, ' ', ' '.join('%16.7g' % v for v in x[indices]), \
'%16.7g' % gm
print >>out, ' ', ' '.join('%16.7g' % v for v in x[indices])
elif type == 'source':
source = problem.unpack(x)
......@@ -1878,14 +1881,7 @@ def export(what, rundirs, type=None, pnames=None, filename=None):
indices = num.array(
[problem.name_to_index(pname) for pname in pnames_take])
if type == 'vector' and what in ('best', 'mean', 'ensemble'):
extra = ['global_misfit']
else:
extra = []
new_header = '# ' + ' '.join(
'%16s' % x for x in pnames_take + extra)
new_header = '# ' + ' '.join('%16s' % x for x in pnames_take)
if type == 'vector' and header != new_header:
print >>out, new_header
......@@ -1894,18 +1890,16 @@ def export(what, rundirs, type=None, pnames=None, filename=None):
indices = None
if what == 'best':
x_best, gm_best = get_best_x_and_gm(problem, xs, misfits)
dump(x_best, gm_best, indices)
dump(get_best_x(problem, xs, misfits), indices)
elif what == 'mean':
x_mean, gm_mean = get_mean_x_and_gm(problem, xs, misfits)
dump(x_mean, gm_mean, indices)
dump(get_mean_x(xs), indices)
elif what == 'ensemble':
gms = problem.global_misfits(misfits)
isort = num.argsort(gms)
for i in isort:
dump(xs[i], gms[i], indices)
dump(xs[i], indices)
elif what == 'stats':
rs = make_stats(problem, xs, misfits, pnames_clean)
......
This diff is collapsed.
......@@ -363,8 +363,8 @@ class Dataset(object):
def get_waveform_restituted(
self,
obj, quantity='displacement',
tmin=None, tmax=None, tpad=0.,
tfade=0., freqlimits=None, deltat=None,
tmin=None, tmax=None, tpad=0., restitute=None,
tfade=0., freqlimits=None, deltat=None,
toffset_noise_extract=0.):
assert quantity == 'displacement' # others not yet implemented
......@@ -373,11 +373,18 @@ class Dataset(object):
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)
if restitute is not None:
resp = self.get_response(tr)
else:
resp = None
return tr.transfer(tfade=tfade, freqlimits=freqlimits,
transfer_function=resp, invert=True)
......@@ -386,7 +393,7 @@ class Dataset(object):
obj, quantity='displacement',
tmin=None, tmax=None, tpad=0.,
tfade=0., freqlimits=None, deltat=None, cache=None,
backazimuth=None,
backazimuth=None,
source=None,
target=None):
......@@ -421,6 +428,9 @@ class Dataset(object):
else:
return obj
syn_test = self.synthetic_test
toffset_noise_extract = 0.0
if syn_test:
......
This diff is collapsed.
import math
import re
import random
import logging
import os
import os.path as op
import numpy as num
from scipy import signal
......@@ -329,7 +327,7 @@ def draw_jointpar_figures(
exclude=None, include=None):
color = 'misfit'
# exclude = ['duration']
# exclude = ['magnitude','rel_moment_iso', 'rel_moment_clvd', 'rel_moment']
# include = ['magnitude', 'rel_moment_iso', 'rel_moment_clvd', 'depth']
neach = 6
figsize = (8, 8)
......@@ -347,10 +345,12 @@ def draw_jointpar_figures(
for ipar in xrange(problem.ncombined):
par = problem.combined[ipar]
lo, hi = bounds[ipar]
# exclude = ['rel_moment_iso']
if lo == hi:
if exclude is None:
exclude = []
if lo is None:
exclude = [ipar]
exclude.append(par.name)
xref = problem.pack(problem.base_source)
......@@ -1365,17 +1365,14 @@ def draw_hudson_figure(model, plt):
mt = mean_source.pyrocko_moment_tensor()
u, v = hudson.project(mt)
try:
beachball.plot_beachball_mpl(
mt, axes,
beachball_type='dc',
position=(u, v),
size=beachballsize,
color_t=color,
zorder=2,
linewidth=0.5)
except beachball.BeachballError, e:
logger.warn(str(e))
beachball.plot_beachball_mpl(
mt, axes,
beachball_type='dc',
position=(u, v),
size=beachballsize,
color_t=color,
zorder=2,
linewidth=0.5)
mt = best_source.pyrocko_moment_tensor()
u, v = hudson.project(mt)
......@@ -1391,17 +1388,14 @@ def draw_hudson_figure(model, plt):
mt = problem.base_source.pyrocko_moment_tensor()
u, v = hudson.project(mt)
try:
beachball.plot_beachball_mpl(
mt, axes,
beachball_type='dc',
position=(u, v),
size=beachballsize,
color_t='red',
zorder=2,
linewidth=0.5)
except beachball.BeachballError, e:
logger.warn(str(e))
beachball.plot_beachball_mpl(
mt, axes,
beachball_type='dc',
position=(u, v),
size=beachballsize,
color_t='red',
zorder=2,
linewidth=0.5)
return [fig]
......@@ -1426,19 +1420,6 @@ plot_dispatch = {
def save_figs(figs, plot_dirname, plotname, formats, dpi):
for fmt in formats:
if fmt not in ['pdf', 'png']:
raise core.GrondError('unavailable output format: %s' % fmt)
assert re.match(r'^[a-z_]+$', plotname)
# remove files from previous runs
pat = re.compile(r'^%s-[0-9]+\.(%s)$' % (plotname, '|'.join(formats)))
if op.exists(plot_dirname):
for entry in os.listdir(plot_dirname):
if pat.match(entry):
os.unlink(op.join(plot_dirname, entry))
fns = []
for ifig, fig in enumerate(figs):
for format in formats:
......
This diff is collapsed.
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment