Skip to content

GitLab

Commit 579da771 authored by Sebastian Heimann's avatar Sebastian Heimann
Browse files

Merge branch 'master' of gitext:heimann/grond into refactor 

Conflicts:
	src/core.py
parents 8a618494 61ad6a4a
Hide whitespace changes
Inline Side-by-side
src/__init__.py
View file @ 579da771
......@@ -3,3 +3,4 @@ from .dataset import * # noqa
from .problems import * # noqa
from .targets import * # noqa
from .meta import * # noqa
from .synthetic_tests import * #noqa
src/core.py
View file @ 579da771
......@@ -51,6 +51,8 @@ class SolverConfig(Object):
default='multivariate_normal')
scatter_scale_transition = Float.T(default=2.0)
scatter_scale = Float.T(default=1.0)
chain_length_factor = Float.T(default=8.0)
compensate_excentricity = Bool.T(default=True)
def get_solver_kwargs(self):
return dict(
......@@ -60,7 +62,9 @@ class SolverConfig(Object):
niter_non_explorative=self.niter_non_explorative,
sampler_distribution=self.sampler_distribution,
scatter_scale_transition=self.scatter_scale_transition,
scatter_scale=self.scatter_scale)
scatter_scale=self.scatter_scale,
chain_length_factor=self.chain_length_factor,
compensate_excentricity=self.compensate_excentricity)
class EngineConfig(HasPaths):
......@@ -232,6 +236,14 @@ def read_config(path):
return config
def write_config(config, path):
basepath = config.get_basepath()
dirname = op.dirname(path) or '.'
config.change_basepath(dirname)
guts.dump(config, filename=path)
config.change_basepath(basepath)
def analyse(problem, niter=1000, show_progress=False):
if niter == 0:
return
......@@ -259,11 +271,10 @@ def analyse(problem, niter=1000, show_progress=False):
isbad_mask = None
for iiter in xrange(niter):
break
while True:
x = []
for i in xrange(npar):
v = rstate.uniform(xbounds[i, 0], xbounds[i, 1])
for ipar in xrange(npar):
v = rstate.uniform(xbounds[ipar, 0], xbounds[ipar, 1])
x.append(v)
try:
......@@ -302,6 +313,37 @@ def analyse(problem, niter=1000, show_progress=False):
balancing_weight=float(weight))
def excentricity_compensated_choice(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_all = math.sqrt(num.sum(
# ((xs[num.newaxis, :, :] - xs[:, num.newaxis, :]) *
# scale[num.newaxis, num.newaxis, :])**2, axis=2))
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)
r = num.random.random()
ichoice = num.searchsorted(num.cumsum(probabilities), r)
ichoice = min(ichoice, xs.shape[0]-1)
return xs[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 solve(problem,
rundir=None,
niter_uniform=1000,
......@@ -310,14 +352,16 @@ def solve(problem,
niter_non_explorative=0,
scatter_scale_transition=2.0,
scatter_scale=1.0,
chain_length_factor=8.0,
xs_inject=None,
sampler_distribution='multivariate_normal',
compensate_excentricity=True,
status=()):
xbounds = num.array(problem.bounds(), dtype=num.float)
npar = xbounds.shape[0]
nlinks_cap = 8 * npar + 1
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
......@@ -384,10 +428,16 @@ def solve(problem,
jchoice = num.argmin(accept_sum)
if phase in ('transition', 'explorative'):
ichoice = num.random.randint(0, nlinks)
xb = xhist[chains_i[jchoice, ichoice]]
if compensate_excentricity:
xchoice = excentricity_compensated_choice(
xhist[chains_i[jchoice, :], :], sbx, 3.)
else:
ichoice = num.random.randint(0, nlinks)
xchoice = xhist[chains_i[jchoice, ichoice]]
else:
xb = mxs[jchoice]
xchoice = mxs[jchoice]
if sampler_distribution == 'multivariate_normal':
ntries_sample = 0
......@@ -397,7 +447,7 @@ def solve(problem,
while True:
ntries_sample += 1
vs = num.random.multivariate_normal(
xb, factor**2 * covs[jchoice])
xchoice, factor**2 * covs[jchoice])
ok_mask = num.logical_and(
xbounds[:, 0] <= vs, vs <= xbounds[:, 1])
......@@ -420,28 +470,38 @@ def solve(problem,
x = vs.tolist()
if sampler_distribution == 'normal':
x = []
for i in xrange(npar):
ntry = 0
while True:
if sbx[i] > 0.:
v = num.random.normal(
xb[i], factor*sbx[i])
else:
v = xb[i]
if xbounds[i, 0] <= v and v <= xbounds[i, 1]:
break
if ntry > 1000:
raise GrondError(
'failed to produce a suitable '
'candidate sample from normal '
'distribution')
ntry += 1
x.append(v)
ncandidates = 1
xcandidates = num.zeros((ncandidates, npar))
for icandidate in xrange(ncandidates):
for ipar in xrange(npar):
ntry = 0
while True:
if sbx[ipar] > 0.:
v = num.random.normal(
xchoice[ipar], factor*sbx[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[jchoice, :], :],
sbx,
factor)
try:
x = problem.preconstrain(x)
......@@ -508,7 +568,7 @@ def solve(problem,
lines = []
if 'state' in status:
lines.append('%i' % iiter)
lines.append('%s, %i' % (problem.name, iiter))
lines.append(''.join('-X'[int(acc)] for acc in accept))
xhist[iiter, :] = x
......@@ -838,16 +898,20 @@ def check(
freqlimits[3] = 0.5/deltat
freqlimits = tuple(freqlimits)
trs_projected, trs_restituted, trs_raw = \
ds.get_waveform(
target.codes,
tmin=tmin+tobs_shift,
tmax=tmax+tobs_shift,
tfade=tfade,
freqlimits=freqlimits,
deltat=deltat,
backazimuth=target.get_backazimuth_for_waveform(),
debug=True)
try:
trs_projected, trs_restituted, trs_raw = \
ds.get_waveform(
target.codes,
tmin=tmin+tobs_shift,
tmax=tmax+tobs_shift,
tfade=tfade,
freqlimits=freqlimits,
deltat=deltat,
backazimuth=target.get_backazimuth_for_waveform(),
debug=True)
except dataset.NotFound, e:
logger.warn(str(e))
continue
trs_projected = copy.deepcopy(trs_projected)
trs_restituted = copy.deepcopy(trs_restituted)
......@@ -975,8 +1039,8 @@ def check(
event.name or util.time_to_str(event.time),
str(e)))
if show_waveforms:
trace.snuffle(trs_all, stations=ds.get_stations(), markers=markers)
if show_waveforms:
trace.snuffle(trs_all, stations=ds.get_stations(), markers=markers)
g_state = {}
......@@ -1246,6 +1310,7 @@ __all__ = '''
Config
load_problem_info_and_data
read_config
write_config
forward
harvest
go
......
src/plot.py
View file @ 579da771
......@@ -164,7 +164,7 @@ class GrondModel(object):
listener()
def draw_sequence_figures(model, plt, misfit_cutoff=None):
def draw_sequence_figures(model, plt, misfit_cutoff=None, sort_by='iteration'):
problem = model.problem
imodels = num.arange(model.nmodels)
......@@ -182,7 +182,13 @@ def draw_sequence_figures(model, plt, misfit_cutoff=None):
isort = num.argsort(gms)[::-1]
imodels = imodels[isort]
if sort_by == 'iteration':
imodels = imodels[isort]
elif sort_by == 'misfit':
imodels = num.arange(imodels.size)
else:
assert False
gms = gms[isort]
gms_softclip = gms_softclip[isort]
xs = xs[isort, :]
......@@ -313,7 +319,7 @@ def draw_sequence_figures(model, plt, misfit_cutoff=None):
def draw_jointpar_figures(
model, plt, misfit_cutoff=None, ibootstrap=None, color=None,
exclude=None, include=None):
exclude=None, include=None, draw_ellipses=False):
color = 'misfit'
# exclude = ['duration']
......@@ -493,17 +499,18 @@ def draw_jointpar_figures(
c=color,
s=msize, alpha=0.5, cmap=cmap, edgecolors='none')
cov = num.cov((xpar.scaled(fx), ypar.scaled(fy)))
evals, evecs = eigh_sorted(cov)
evals = num.sqrt(evals)
ell = patches.Ellipse(
xy=(num.mean(xpar.scaled(fx)), num.mean(ypar.scaled(fy))),
width=evals[0]*2,
height=evals[1]*2,
angle=num.rad2deg(num.arctan2(evecs[1][0], evecs[0][0])))
ell.set_facecolor('none')
axes.add_artist(ell)
if draw_ellipses:
cov = num.cov((xpar.scaled(fx), ypar.scaled(fy)))
evals, evecs = eigh_sorted(cov)
evals = num.sqrt(evals)
ell = patches.Ellipse(
xy=(num.mean(xpar.scaled(fx)), num.mean(ypar.scaled(fy))),
width=evals[0]*2,
height=evals[1]*2,
angle=num.rad2deg(num.arctan2(evecs[1][0], evecs[0][0])))
ell.set_facecolor('none')
axes.add_artist(ell)
fx = problem.extract(xref, jpar)
fy = problem.extract(xref, ipar)
......
src/synthetic_tests.py
View file @ 579da771
......@@ -4,6 +4,9 @@ from pyrocko import trace
from pyrocko.guts import (Object, Dict, String, Float, Bool, Int)
guts_prefix = 'grond'
class SyntheticWaveformNotAvailable(Exception):
pass
......@@ -13,6 +16,7 @@ class SyntheticTest(Object):
respect_data_availability = Bool.T(default=False)
real_noise_scale = Float.T(default=0.0)
white_noise_scale = Float.T(default=0.0)
relative_white_noise_scale = Float.T(default=0.0)
random_response_scale = Float.T(default=0.0)
real_noise_toffset = Float.T(default=-3600.)
random_seed = Int.T(optional=True)
......@@ -56,26 +60,33 @@ class SyntheticTest(Object):
for iresult, result in enumerate(results):
tr = result.trace.pyrocko_trace()
tfade = tr.tmax - tr.tmin
tr_orig = tr.copy()
tr.extend(tr.tmin - tfade, tr.tmax + tfade)
rstate = num.random.RandomState(
(self.random_seed or 0) + iresult)
if self.random_response_scale != 0:
tf = RandomResponse(scale=self.random_response_scale)
rstate = num.random.RandomState(
(self.random_seed or 0) + iresult)
tf.set_random_state(rstate)
tr = tr.transfer(
tfade=tfade,
transfer_function=tf)
if self.white_noise_scale != 0.0:
rstate = num.random.RandomState(
(self.random_seed or 0) + iresult)
u = rstate.normal(
scale=self.white_noise_scale,
size=tr.data_len())
tr.ydata += u
if self.relative_white_noise_scale != 0.0:
u = rstate.normal(
scale=self.relative_white_noise_scale * num.std(
tr_orig.ydata),
size=tr.data_len())
tr.ydata += u
synthetics[result.trace.codes] = tr
self._synthetics = synthetics
......@@ -111,3 +122,9 @@ class RandomResponse(trace.FrequencyResponse):
return 1.0 + (
self._rstate.normal(scale=self.scale, size=n) +
0.0J * self._rstate.normal(scale=self.scale, size=n))
__all__ = '''
SyntheticTest
SyntheticWaveformNotAvailable
'''.split()
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