Commit 7609fc91 authored by Sebastian Heimann's avatar Sebastian Heimann

highscoreoptimizer movies

parent 9f9aee18
......@@ -559,15 +559,19 @@ def command_movie(args):
parser, options, args = cl_parse('movie', args, setup)
if len(args) != 3:
if len(args) != 4:
help_and_die(parser, 'one argument required')
run_path, xpar_name, ypar_name = args
run_path, xpar_name, ypar_name, movie_filename_template = args
from grond import plot
movie_filename = movie_filename_template % {
'xpar': xpar_name,
'ypar': ypar_name}
try:
plot.make_movie(run_path, xpar_name, ypar_name)
plot.make_movie(run_path, xpar_name, ypar_name, movie_filename)
except grond.GrondError as e:
die(str(e))
......
from .base import * # noqa
from .highscore.optimizer import * # noqa
from .highscore.plot import * # noqa
......@@ -28,7 +28,7 @@ def excentricity_compensated_probabilities(xs, sbx, factor):
((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)))
# print(num.sort(num.sum(distances_sqr_all < 1.0, axis=1)))
probabilities /= num.sum(probabilities)
return probabilities
......@@ -261,49 +261,62 @@ class Chains(object):
(self.nchains, nlinks_cap), num.int)
self.nlinks = 0
self.accept_sum = num.zeros(self.nchains, dtype=num.int)
self.extend(0, history.nmodels, history.models)
self.nread = 0
history.add_listener(self)
self.bootstrap_weights = num.vstack((
num.ones((1, self.problem.ntargets)),
bootstrap_weights))
def append(self, iiter, model, misfits):
gbms = self.problem.combine_misfits(misfits, self.bootstrap_weights)
def goto(self, n=None):
if n is None:
n = self.history.nmodels
self.chains_m[:, self.nlinks] = gbms
self.chains_i[:, self.nlinks] = iiter
n = min(self.history.nmodels, n)
self.nlinks += 1
chains_m = self.chains_m
chains_i = self.chains_i
assert self.nread <= n
for ichain in range(chains_m.shape[0]):
isort = num.argsort(chains_m[ichain, :self.nlinks])
chains_m[ichain, :self.nlinks] = chains_m[ichain, isort]
chains_i[ichain, :self.nlinks] = chains_i[ichain, isort]
while self.nread < n:
misfits = self.history.misfits[self.nread, :]
gbms = self.problem.combine_misfits(
misfits, self.bootstrap_weights)
if self.nlinks == self.nlinks_cap:
accept = (chains_i[:, self.nlinks_cap-1] != iiter).astype(num.int)
self.nlinks -= 1
else:
accept = num.ones(self.nchains, dtype=num.int)
self.chains_m[:, self.nlinks] = gbms
self.chains_i[:, self.nlinks] = n-1
self.nlinks += 1
chains_m = self.chains_m
chains_i = self.chains_i
for ichain in range(chains_m.shape[0]):
isort = num.argsort(chains_m[ichain, :self.nlinks])
chains_m[ichain, :self.nlinks] = chains_m[ichain, isort]
chains_i[ichain, :self.nlinks] = chains_i[ichain, isort]
if self.nlinks == self.nlinks_cap:
accept = (
chains_i[:, self.nlinks_cap-1] != n-1).astype(num.int)
self.nlinks -= 1
else:
accept = num.ones(self.nchains, dtype=num.int)
self.accept_sum += accept
self.accept_sum += accept
self.nread += 1
return accept
def append(self, iiter, model, misfits):
self.goto(iiter)
def extend(self, ioffset, n, models, misfits):
for i in range(ioffset, ioffset+n):
self.append(i, models[i-ioffset, :], misfits[i-ioffset, :, :])
self.goto(ioffset + n)
def models(self, ichain=None):
def indices(self, ichain):
if ichain is not None:
return self.history.models[
self.chains_i[ichain, :self.nlinks], :]
return self.chains_i[ichain, :self.nlinks]
else:
return self.history.models[
self.chains_i[:, :self.nlinks].ravel(), :]
return self.chains_i[:, :self.nlinks].ravel()
def models(self, ichain=None):
return self.history.models[self.indices(ichain), :]
def model(self, ichain, ilink):
return self.history.models[self.chains_i[ichain, ilink], :]
......@@ -365,7 +378,7 @@ class HighScoreOptimizer(Optimizer):
return bms
def setup_chains(self, problem, history):
def chains(self, problem, history):
bootstrap_weights = self.get_bootstrap_weights(problem)
nlinks_cap = int(round(
......@@ -375,30 +388,28 @@ class HighScoreOptimizer(Optimizer):
problem, history, self.nbootstrap+1, nlinks_cap,
bootstrap_weights=bootstrap_weights)
def get_sampler_phase(self, iiter):
niter = 0
for phase in self.sampler_phases:
if iiter < niter + phase.niterations:
return phase, iiter - niter
niter += phase.niterations
assert False, 'out of bounds'
def optimize(self, problem, rundir=None):
if rundir is not None:
self.dump(filename=op.join(rundir, 'optimizer.yaml'))
history = ModelHistory(problem, path=rundir, mode='w')
chains = self.chains(problem, history)
chains = self.setup_chains(problem, history)
phases = list(self.sampler_phases)
phase = phases.pop(0)
niter = sum(phase.niterations for phase in self.sampler_phases)
niter = self.niterations
isbad_mask = None
iiter_phase_start = 0
iiter = 0
while iiter < niter:
if iiter - iiter_phase_start == phase.niterations:
phase = phases.pop(0)
iiter_phase_start = iiter
iiter_phase = iiter - iiter_phase_start
for iiter in range(niter):
phase, iiter_phase = self.get_sampler_phase(iiter)
x = phase.get_sample(problem, iiter_phase, chains)
if isbad_mask is not None and num.any(isbad_mask):
......@@ -434,16 +445,16 @@ class HighScoreOptimizer(Optimizer):
history.append(x, misfits)
iiter += 1
@property
def niterations(self):
return sum([ph.niterations for ph in self.sampler_phases])
def get_movie_maker(self, problem, history, xpar_name, ypar_name):
def get_movie_maker(
self, problem, history, xpar_name, ypar_name, movie_filename):
from . import plot
return plot.HighScoreOptimizerPlot(
self, problem, history, xpar_name, ypar_name)
self, problem, history, xpar_name, ypar_name, movie_filename)
class HighScoreOptimizerConfig(OptimizerConfig):
......
This diff is collapsed.
......@@ -2118,12 +2118,12 @@ def plot_result(dirname, plotnames_want,
plt.show()
def make_movie(dirname, xpar_name, ypar_name):
def make_movie(dirname, xpar_name, ypar_name, movie_filename):
optimizer_fn = op.join(dirname, 'optimizer.yaml')
optimizer = guts.load(filename=optimizer_fn)
problem = load_problem_info(dirname)
history = ModelHistory(problem, path=dirname)
movie_maker = optimizer.get_movie_maker(
problem, history, xpar_name, ypar_name)
problem, history, xpar_name, ypar_name, movie_filename)
movie_maker.render()
......@@ -3,11 +3,53 @@ import numpy as num
from pyrocko.guts import Object, Float, Dict, List, String
from pyrocko import gf
from grond import Problem, Parameter, MisfitTarget
from grond import Problem, Parameter, MisfitTarget, HighScoreOptimizerPlot
guts_prefix = 'grond.toy'
class ToyOptimizerPlot(HighScoreOptimizerPlot):
def set_source(self, source):
self._source = source
def set_targets(self, targets):
self._targets = targets
def set_contour_data(self, contour_data):
self._contour_data = contour_data
def set_limits(self):
self.axes.set_xlim(-10., 10.)
self.axes.set_ylim(-10., 10.)
def start(self):
HighScoreOptimizerPlot.start(self)
x = [getattr(t, self.xpar_name) for t in self._targets]
y = [getattr(t, self.ypar_name) for t in self._targets]
self.axes.plot(x, y, '^', color='black')
for ibootstrap, (xc, yc, zc) in enumerate(
self._contour_data['east', 'depth']):
zmin = num.min(zc)
if self.optimizer.nbootstrap < 5:
alpha = 1.0
else:
alpha = 0.5
self.axes.contour(
xc, yc, zc, [zmin + 0.01],
colors=[self.bcolors[ibootstrap]], alpha=alpha)
self.axes.plot(
getattr(self._source, self.xpar_name),
getattr(self._source, self.ypar_name),
'*', color='black')
class ToyTarget(MisfitTarget):
north = Float.T()
east = Float.T()
......
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