import numpy as np
from logging import info
import pickle
from .sampling import compute_sample
[docs]
def execute(self, **kwargs):
"""Performs a variational inversion given a minimizer method and a
simulator (i.e. a function to minimize and its gradient)
Args:
self (Plugin): definition of the mode set-up
"""
# Working directory
workdir = self.workdir
# Control vector
controlvect = self.controlvect
# Observation operator
obsoper = self.obsoperator
# Simulation window
datei = self.datei
datef = self.datef
# Minimizer
minimizer = self.minimizer
# Simulator
simulator = self.simulator
# Obsvervation vector
obsvect = self.obsvect
# Check that the observation vector is not empty
if obsvect.obsvect_mask.sum() == 0:
raise Exception("WARNING: trying to run an inversion with no "
"observation points! Please check your monitor file")
# If no Monte Carlo estimation of posterior uncertainties,
# Do simple inversion
if getattr(self, "montecarlo", None) is None:
# Some verbose
towrite = """
Running a variational inversion with the following modules:
Minimizer: {}
Simulator: {}
Size of the control vector: {}
Size of the observation vector: {}
""".format(
minimizer.plugin.name, simulator.plugin.name,
controlvect.dim, obsvect.dim
)
info(towrite)
# Initial run of the simulator as a starting point for the minimizer
costinit, gradinit = simulator.simul(
controlvect.chi, run_id=-1, **kwargs)
zgnorm_init = np.sqrt(np.sum(gradinit ** 2))
info("Nb of observations: " + str(len(obsoper.obsvect.yobs)))
info("Initial cost: " + str(costinit))
info("Initial gradient norm: " + str(zgnorm_init))
# Runs the minimizer
chiopt = minimizer.minimize(costinit, gradinit,
controlvect.chi, **kwargs)
# Save uncertainty reduction to the controlvect
if getattr(minimizer, "save_uncertainties", False):
# Convert the eigen vectors back to the control space
eigenvects = []
for pevec in minimizer.pevecs:
eigenvects.append(controlvect.sqrtbprod(pevec, **kwargs)
- controlvect.xb)
eigenvects = np.array(eigenvects)
cntrl_file = "{}/evectors.pickle".format(workdir)
with open(cntrl_file, "wb") as f:
pickle.dump(eigenvects, f, pickle.HIGHEST_PROTOCOL)
# Computes Pa from eigen vectors
# too big, do this as post-processing by block
# controlvect.pa = (controlvect.build_b(controlvect)
# - (eigenvects.T.dot(eigenvects)))
# Save results
controlvect.x = controlvect.sqrtbprod(chiopt, **kwargs)
controlvect.dump(
"{}/controlvect_final.pickle".format(workdir),
to_netcdf=controlvect.save_out_netcdf or self.save_out_netcdf,
dir_netcdf="{}/controlvect/".format(workdir),
)
# Last call to the simulator for final diagnostics
costend, gradend = simulator.simul(chiopt, run_id=-2, **kwargs)
zgnorm_final = np.sqrt(np.dot(gradend, gradend))
info("Original cost: " + str(costinit))
info("Final cost: " + str(costend))
info("Initial gradient norm: " + str(zgnorm_init))
info("Final gradient norm: " + str(zgnorm_final))
info("Ratio final/initial gradient norm: " +
str(zgnorm_final / zgnorm_init))
# Compute chi-square test
chi_square = 2 * costend / len(obsoper.obsvect.yobs)
with open("{}/chi_square_test.txt".format(workdir), "w") as f:
f.write(str(chi_square))
# Return values
return controlvect, obsvect
# If Monte Carlo, call individual inversions
montecarlo = self.montecarlo
ensemble_dir = "{}/ensemble/".format(workdir)
nsample = montecarlo.nsample
# Build ensemble of x
if getattr(montecarlo, "perturb_x", True):
if getattr(montecarlo, "seed_x", None) is not None:
np.random.seed(montecarlo.seed_x)
x_sample = np.random.normal(size=(nsample, controlvect.dim))
x_sample = controlvect.sqrtbprod(x_sample.T, ensemble=True).T
else:
x_sample = np.ones((nsample, controlvect.dim)) \
* controlvect.xb[np.newaxis]
# Build ensemble of y
if getattr(montecarlo, "perturb_y", True):
if getattr(montecarlo, "seed_y", None) is not None:
np.random.seed(montecarlo.seed_y)
y_sample = np.random.normal(size=(nsample, obsvect.dim))
for k, y in enumerate(y_sample):
y_sample[k] = obsvect.rinvprod(y, inverse=False)
else:
y_sample = np.ones((nsample, obsvect.dim)) \
* obsvect.yobs[np.newaxis]
# Compute simulation for reference control vector as well
x_sample = np.append(controlvect.xb[np.newaxis], x_sample, axis=0)
y_sample = np.append(obsvect.yobs[np.newaxis], y_sample, axis=0)
nvalid = compute_sample(self, controlvect, obsvect, x_sample, y_sample)
# Load posterior control vectors
h_dir = "{workdir}/ensemble/xa_emsemble/".format(workdir=workdir)
xa_sample = np.zeros((nvalid + 1, controlvect.dim))
for k in range(nvalid + 1):
xa_file = "{}/controlvect_final_{:04d}.pickle".format(h_dir, k)
controlvect.load(xa_file)
xa_sample[k] = controlvect.x
xa_mean = xa_sample[1:].mean(axis=0)
dx = xa_sample[1:] - xa_mean
info("Pa calculation")
if getattr(montecarlo, "fullPa_calc", True):
pa = dx.T.dot(dx) / (nvalid - 1)
else:
pa = np.sum(dx ** 2, axis=0) / (nvalid - 1)
info(dx.shape)
info(pa.shape)
controlvect.x[:] = xa_mean
controlvect.pa = pa
# Save results
controlvect.dump(
"{}/controlvect_final.pickle".format(workdir),
to_netcdf=controlvect.save_out_netcdf or self.save_out_netcdf,
dir_netcdf="{}/controlvect/".format(workdir),
)
# Stop here if not requested to re-run prior and posterior forwards
if not montecarlo.aggregate_results:
return controlvect, obsvect
# Run a prior forward
controlvect.chi = controlvect.sqrtbprod_ad(0 * controlvect.xb)
costinit, gradinit = simulator.simul(controlvect.chi, run_id=-1,
**kwargs)
zgnorm = np.sqrt(np.sum(gradinit ** 2))
# Re-run a forward with the posterior mean to compare with observations
controlvect.x[:] = xa_mean
controlvect.chi = controlvect.sqrtbprod_ad(
controlvect.x - controlvect.xb, inverse=True)
costfinal, gradfinal = simulator.simul(controlvect.chi, run_id=-2,
**kwargs)
zgnorm = np.sqrt(np.dot(gradfinal, gradfinal)) / zgnorm
info("Original cost: " + str(costinit))
info("Final cost: " + str(costfinal))
info("Ratio final/initial gradient norm: " + str(zgnorm))
# Compute chi-square test
chi_square = 2 * costfinal / len(obsoper.obsvect.yobs)
with open("{}/chi_square_test.txt".format(workdir), "w") as f:
f.write(str(chi_square))
info("The inversion exited with a chi-square test of: {}"
.format(chi_square))
return controlvect, obsvect
