Commit bf3a7f7f authored by Sebastian Heimann's avatar Sebastian Heimann

reimplement highscore movie (wip)

parent e352e9dc
......@@ -11,6 +11,7 @@ from pyrocko import util, marker
from pyrocko.gf import Range
import grond
import grond.toy
logger = logging.getLogger('main')
km = 1e3
......@@ -31,6 +32,7 @@ subcommand_descriptions = {
'forward': 'run forward modelling',
'harvest': 'manually run harvesting',
'plot': 'plot optimization result',
'movie': 'visualize optimizer evolution',
'baraddur': 'start Barad-dur plotting, watching this directory',
'export': 'export results',
'qc-polarization': 'check sensor orientations with polarization analysis',
......@@ -46,6 +48,7 @@ subcommand_usages = {
'forward <configfile> <eventnames> ... [options]'),
'harvest': 'harvest <rundir> [options]',
'plot': 'plot <plotnames> <rundir> [options]',
'movie': 'movie <rundir> [options]',
'baraddur': 'plot-server',
'export': 'export (best|mean|ensemble|stats) <rundirs> ... [options]',
'qc-polarization': 'qc-polarization <configfile> <eventname> '
......@@ -70,6 +73,7 @@ Subcommands:
forward %(forward)s
harvest %(harvest)s
plot %(plot)s
movie %(movie)s
export %(export)s
qc-polarization %(qc_polarization)s
baraddur %(baraddur)s
......@@ -548,6 +552,27 @@ selected by specifying a comma-separated list.''' % (
die(str(e))
def command_movie(args):
def setup(parser):
pass
parser, options, args = cl_parse('movie', args, setup)
if len(args) != 3:
help_and_die(parser, 'one argument required')
run_path, xpar_name, ypar_name = args
from grond import plot
try:
plot.make_movie(run_path, xpar_name, ypar_name)
except grond.GrondError as e:
die(str(e))
def command_baraddur(args):
from grond.baraddur import Baraddur, BaraddurConfig
import signal
......
......@@ -439,6 +439,11 @@ class HighScoreOptimizer(Optimizer):
def niterations(self):
return sum([ph.niterations for ph in self.sampler_phases])
def get_movie_maker(self, problem, history, xpar_name, ypar_name):
from . import plot
return plot.HighScoreOptimizerPlot(
self, problem, history, xpar_name, ypar_name)
class HighScoreOptimizerConfig(OptimizerConfig):
......
import math
import numpy as num
from matplotlib import pyplot as plt
from pyrocko.plot import mpl_init, mpl_margins
from grond import plot
class HighScoreOptimizerPlot(object):
def __init__(self, optimizer, problem, history, xpar_name, ypar_name):
self.optimizer = optimizer
self.problem = problem
self.history = history
self.xpar_name = xpar_name
self.ypar_name = ypar_name
self.fontsize = 10.
self.movie_filename = None
self.show = True
self.iiter = None
def start(self):
self.iiter = 0
nfx = 1
nfy = 1
problem = self.problem
ixpar = problem.name_to_index(self.xpar_name)
iypar = problem.name_to_index(self.ypar_name)
mpl_init(fontsize=self.fontsize)
fig = plt.figure(figsize=(9.6, 5.4))
labelpos = mpl_margins(fig, nw=nfx, nh=nfy, w=7., h=5., wspace=7.,
hspace=2., units=self.fontsize)
xpar = problem.parameters[ixpar]
ypar = problem.parameters[iypar]
if xpar.unit == ypar.unit:
axes = fig.add_subplot(nfy, nfx, 1, aspect=1.0)
else:
axes = fig.add_subplot(nfy, nfx, 1)
labelpos(axes, 2.5, 2.0)
axes.set_xlabel(xpar.get_label())
axes.set_ylabel(ypar.get_label())
axes.get_xaxis().set_major_locator(plt.MaxNLocator(4))
axes.get_yaxis().set_major_locator(plt.MaxNLocator(4))
xref = problem.get_xref()
axes.axhline(xpar.scaled(xref[ixpar]), color='black', alpha=0.3)
axes.axvline(ypar.scaled(xref[iypar]), color='black', alpha=0.3)
self.fig = fig
self.problem = problem
self.xpar = xpar
self.ypar = ypar
self.axes = axes
self.ixpar = ixpar
self.iypar = iypar
from matplotlib import colors
n = self.optimizer.nbootstrap + 1
hsv = num.vstack((
num.random.uniform(0., 1., n),
num.random.uniform(0.5, 0.9, n),
num.repeat(0.7, n))).T
self.bcolors = colors.hsv_to_rgb(hsv[num.newaxis, :, :])[0, :, :]
self.bcolors[0, :] = [0., 0., 0.]
bounds = self.problem.get_combined_bounds()
self.xlim = plot.fixlim(*xpar.scaled(bounds[ixpar]))
self.ylim = plot.fixlim(*ypar.scaled(bounds[iypar]))
self.set_limits()
from matplotlib.colors import LinearSegmentedColormap
self.cmap = LinearSegmentedColormap.from_list('probability', [
(1.0, 1.0, 1.0),
(0.5, 0.9, 0.6)])
self.writer = None
if self.movie_filename:
from matplotlib.animation import FFMpegWriter
metadata = dict(title=problem.name, artist='Grond')
self.writer = FFMpegWriter(
fps=30,
metadata=metadata,
codec='libx264',
bitrate=200000)
self.writer.setup(self.fig, self.movie_filename, dpi=200)
#if self.show:
#plt.ion()
#plt.show()
def set_limits(self):
self.axes.set_xlim(*self.xlim)
self.axes.set_ylim(*self.ylim)
def draw_frame(self, iiter):
msize = 15.
self.axes.cla()
fx = self.problem.extract(self.history.models[:iiter], self.ixpar)
fy = self.problem.extract(self.history.models[:iiter], self.iypar)
self.axes.scatter(
self.xpar.scaled(fx),
self.ypar.scaled(fy),
color='black',
s=msize * 0.15, alpha=0.2, edgecolors='none')
def update(self, xhist, chains_i, ibase, jchoice, local_sxs, factor, phase,
compensate_excentricity):
msize = 15.
iiter_frame = len(xhist)
self.axes.cla()
if jchoice is not None and local_sxs is not None:
nx = 100
ny = 100
sx = local_sxs[jchoice][self.ixpar] * factor
sy = local_sxs[jchoice][self.iypar] * factor
p = num.zeros((ny, nx))
for j in [jchoice]: # range(self.problem.nbootstrap+1):
if compensate_excentricity:
ps = core.excentricity_compensated_probabilities(
xhist[chains_i[j, :], :], local_sxs[jchoice], 2.)
else:
ps = num.ones(chains_i.shape[1])
bounds = self.get_combined_bounds()
x = num.linspace(
bounds[self.ixpar][0], bounds[self.ixpar][1], nx)
y = num.linspace(
bounds[self.iypar][0], bounds[self.iypar][1], ny)
for ichoice in range(chains_i.shape[1]):
iiter = chains_i[j, ichoice]
vx = xhist[iiter, self.ixpar]
vy = xhist[iiter, self.iypar]
pdfx = 1.0 / math.sqrt(2.0 * sx**2 * math.pi) * num.exp(
-(x - vx)**2 / (2.0 * sx**2))
pdfy = 1.0 / math.sqrt(2.0 * sy**2 * math.pi) * num.exp(
-(y - vy)**2 / (2.0 * sy**2))
p += ps[ichoice] * pdfx[num.newaxis, :] * \
pdfy[:, num.newaxis]
self.axes.pcolormesh(x, y, p, cmap=self.cmap)
fx = self.problem.extract(xhist, self.ixpar)
fy = self.problem.extract(xhist, self.iypar)
self.axes.scatter(
self.xpar.scaled(fx),
self.ypar.scaled(fy),
color='black',
s=msize * 0.15, alpha=0.2, edgecolors='none')
for ibootstrap in range(self.optimizer.nbootstrap + 1):
iiters = chains_i[ibootstrap, :]
fx = self.problem.extract(xhist[iiters, :], self.ixpar)
fy = self.problem.extract(xhist[iiters, :], self.iypar)
nfade = 20
factors = 1.0 + 5.0 * (num.maximum(
0.0, iiters - (xhist.shape[0] - nfade)) / nfade)**2
msizes = msize * factors
self.axes.scatter(
self.xpar.scaled(fx),
self.ypar.scaled(fy),
color=self.bcolors[ibootstrap],
s=msizes, alpha=0.5, edgecolors='none')
# if ibase is not None:
# fx = self.problem.extract(
# xhist[(ibase, -1), :], self.ixpar)
# fy = self.problem.extract(
# xhist[(ibase, -1), :], self.iypar)
# self.axes.plot(
# self.xpar.scaled(fx),
# self.ypar.scaled(fy),
# color='black')
fx = self.problem.extract(xhist[-1:, :], self.ixpar)
fy = self.problem.extract(xhist[-1:, :], self.iypar)
self.axes.scatter(
self.xpar.scaled(fx),
self.ypar.scaled(fy),
s=msize * 5.0,
color='none',
edgecolors=self.bcolors[ibootstrap])
self.axes.annotate(
'%i (%s)' % (iiter_frame, phase),
xy=(0., 1.),
xycoords='axes fraction',
xytext=(self.fontsize/2., -self.fontsize/2.),
textcoords='offset points',
ha='left',
va='top',
fontsize=self.fontsize,
fontstyle='normal')
self.set_limits()
self.post_update()
if self.writer:
self.writer.grab_frame()
if self.show:
plt.draw()
def post_update(self):
pass
def finish(self):
if self.writer:
self.writer.finish()
if self.show:
plt.show()
#plt.ioff()
def render(self):
self.start()
self.draw_frame(100)
self.finish()
......@@ -2118,232 +2118,12 @@ def plot_result(dirname, plotnames_want,
plt.show()
class SolverPlot(object):
def __init__(
self, plt, xpar_name, ypar_name, show=False, update_every=1,
movie_filename=None):
self.plt = plt
self.xpar_name = xpar_name
self.ypar_name = ypar_name
self.movie_filename = movie_filename
self.show = show
self.update_every = update_every
self.fontsize = 10
def want_to_update(self, iiter):
return iiter % self.update_every == 0
def start(self, problem):
nfx = 1
nfy = 1
ixpar = problem.name_to_index(self.xpar_name)
iypar = problem.name_to_index(self.ypar_name)
mpl_init(fontsize=self.fontsize)
fig = plt.figure(figsize=(9.6, 5.4))
labelpos = mpl_margins(fig, nw=nfx, nh=nfy, w=7., h=5., wspace=7.,
hspace=2., units=self.fontsize)
xpar = problem.parameters[ixpar]
ypar = problem.parameters[iypar]
if xpar.unit == ypar.unit:
axes = fig.add_subplot(nfy, nfx, 1, aspect=1.0)
else:
axes = fig.add_subplot(nfy, nfx, 1)
labelpos(axes, 2.5, 2.0)
axes.set_xlabel(xpar.get_label())
axes.set_ylabel(ypar.get_label())
axes.get_xaxis().set_major_locator(plt.MaxNLocator(4))
axes.get_yaxis().set_major_locator(plt.MaxNLocator(4))
xref = problem.get_xref()
axes.axhline(xpar.scaled(xref[ixpar]), color='black', alpha=0.3)
axes.axvline(ypar.scaled(xref[iypar]), color='black', alpha=0.3)
self.fig = fig
self.problem = problem
self.xpar = xpar
self.ypar = ypar
self.axes = axes
self.ixpar = ixpar
self.iypar = iypar
from matplotlib import colors
n = problem.nbootstrap + 1
hsv = num.vstack((
num.random.uniform(0., 1., n),
num.random.uniform(0.5, 0.9, n),
num.repeat(0.7, n))).T
self.bcolors = colors.hsv_to_rgb(hsv[num.newaxis, :, :])[0, :, :]
self.bcolors[0, :] = [0., 0., 0.]
bounds = self.problem.get_combined_bounds()
self.xlim = fixlim(*xpar.scaled(bounds[ixpar]))
self.ylim = fixlim(*ypar.scaled(bounds[iypar]))
self.set_limits()
from matplotlib.colors import LinearSegmentedColormap
self.cmap = LinearSegmentedColormap.from_list('probability', [
(1.0, 1.0, 1.0),
(0.5, 0.9, 0.6)])
self.writer = None
if self.movie_filename:
from matplotlib.animation import FFMpegWriter
metadata = dict(title=problem.name, artist='Grond')
self.writer = FFMpegWriter(
fps=30,
metadata=metadata,
codec='libx264',
bitrate=200000)
self.writer.setup(self.fig, self.movie_filename, dpi=200)
if self.show:
plt.ion()
plt.show()
def set_limits(self):
self.axes.set_xlim(*self.xlim)
self.axes.set_ylim(*self.ylim)
def update(self, xhist, chains_i, ibase, jchoice, local_sxs, factor, phase,
compensate_excentricity):
msize = 15.
iiter_frame = len(xhist)
self.axes.cla()
if jchoice is not None and local_sxs is not None:
nx = 100
ny = 100
sx = local_sxs[jchoice][self.ixpar] * factor
sy = local_sxs[jchoice][self.iypar] * factor
p = num.zeros((ny, nx))
for j in [jchoice]: # range(self.problem.nbootstrap+1):
if compensate_excentricity:
ps = core.excentricity_compensated_probabilities(
xhist[chains_i[j, :], :], local_sxs[jchoice], 2.)
else:
ps = num.ones(chains_i.shape[1])
bounds = self.get_combined_bounds()
x = num.linspace(
bounds[self.ixpar][0], bounds[self.ixpar][1], nx)
y = num.linspace(
bounds[self.iypar][0], bounds[self.iypar][1], ny)
for ichoice in range(chains_i.shape[1]):
iiter = chains_i[j, ichoice]
vx = xhist[iiter, self.ixpar]
vy = xhist[iiter, self.iypar]
pdfx = 1.0 / math.sqrt(2.0 * sx**2 * math.pi) * num.exp(
-(x - vx)**2 / (2.0 * sx**2))
pdfy = 1.0 / math.sqrt(2.0 * sy**2 * math.pi) * num.exp(
-(y - vy)**2 / (2.0 * sy**2))
p += ps[ichoice] * pdfx[num.newaxis, :] * \
pdfy[:, num.newaxis]
self.axes.pcolormesh(x, y, p, cmap=self.cmap)
fx = self.problem.extract(xhist, self.ixpar)
fy = self.problem.extract(xhist, self.iypar)
self.axes.scatter(
self.xpar.scaled(fx),
self.ypar.scaled(fy),
color='black',
s=msize * 0.15, alpha=0.2, edgecolors='none')
for ibootstrap in range(self.problem.nbootstrap + 1):
iiters = chains_i[ibootstrap, :]
fx = self.problem.extract(xhist[iiters, :], self.ixpar)
fy = self.problem.extract(xhist[iiters, :], self.iypar)
nfade = 20
factors = 1.0 + 5.0 * (num.maximum(
0.0, iiters - (xhist.shape[0] - nfade)) / nfade)**2
msizes = msize * factors
self.axes.scatter(
self.xpar.scaled(fx),
self.ypar.scaled(fy),
color=self.bcolors[ibootstrap],
s=msizes, alpha=0.5, edgecolors='none')
# if ibase is not None:
# fx = self.problem.extract(
# xhist[(ibase, -1), :], self.ixpar)
# fy = self.problem.extract(
# xhist[(ibase, -1), :], self.iypar)
# self.axes.plot(
# self.xpar.scaled(fx),
# self.ypar.scaled(fy),
# color='black')
fx = self.problem.extract(xhist[-1:, :], self.ixpar)
fy = self.problem.extract(xhist[-1:, :], self.iypar)
self.axes.scatter(
self.xpar.scaled(fx),
self.ypar.scaled(fy),
s=msize * 5.0,
color='none',
edgecolors=self.bcolors[ibootstrap])
self.axes.annotate(
'%i (%s)' % (iiter_frame, phase),
xy=(0., 1.),
xycoords='axes fraction',
xytext=(self.fontsize/2., -self.fontsize/2.),
textcoords='offset points',
ha='left',
va='top',
fontsize=self.fontsize,
fontstyle='normal')
self.set_limits()
self.post_update()
if self.writer:
self.writer.grab_frame()
if self.show:
plt.draw()
def post_update(self):
pass
def finish(self):
if self.writer:
self.writer.finish()
def make_movie(dirname, xpar_name, ypar_name):
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)
if self.show:
plt.ioff()
movie_maker.render()
from .waveform.target import * # noqa
from .satellite.target import * # noqa
from .base import * # noqa
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