import numpy as np
import copy
from logging import info
[docs]
def testing_AD_whole(self, accuracy, testspace, **kwargs):
"""Verify the adjoint identity ``<H·dx, H·dx> == <dx, H*·H·dx>``.
Runs either the tangent-linear operator or two finite-difference forward
runs (when ``use_forward`` is set) to obtain ``H·dx``, then runs the
adjoint to obtain ``H*·H·dx``. The two inner products are printed and
their relative discrepancy expressed as a multiple of machine epsilon is
returned.
Args:
self (Plugin): mode plugin with all configuration attributes.
accuracy (float): machine epsilon (``np.finfo(np.float64).eps``).
testspace (str): ``'control'`` to work in control space or ``'chi'``
to include the square-root-B mapping.
**kwargs: forwarded verbatim to the observation operator.
Returns:
float: ``floor(|<dx, H*(H·dx)> / <H·dx, H·dx> - 1| / epsilon)``,
the test ratio in multiples of machine accuracy. Values up to a
few tens are typically acceptable.
"""
# Working directory
workdir = self.workdir
# Control vector
controlvect = self.controlvect
# Observation operator
obsoper = self.obsoperator
# Obsvervation vector
obsvect = self.obsvect
# Simulation window
datei = self.datei
datef = self.datef
# Some verbose
info("Computing the test of the adjoint")
# Running the tangent linear code of the model: H(xb)(dx)
if not self.use_forward:
obsvect = obsoper.obsoper(
controlvect,
obsvect,
"tl",
datei=datei,
datef=datef,
workdir=workdir,
reload_results=self.reload_results,
check_transforms=self.check_transforms,
**kwargs
)
else:
# Reference simulation: H(xb)
obsvect_ref = obsoper.obsoper(
controlvect,
obsvect,
"fwd",
datei=datei,
datef=datef,
workdir=workdir,
reload_results=self.reload_results,
check_transforms=self.check_transforms,
run_id=1,
**kwargs
)
ysim_ref = copy.deepcopy(obsvect_ref.ysim)
# Save check transforms
if self.check_transforms:
data_tl_save = copy.deepcopy(obsoper.data_tl)
del obsoper.data_tl
# Perturbed simulation: H(xb + dx)
controlvect.x = copy.deepcopy(controlvect.xb + controlvect.dx_save)
obsvect_dy = obsoper.obsoper(
controlvect,
obsvect,
"fwd",
datei=datei,
datef=datef,
workdir=workdir,
reload_results=self.reload_results,
check_transforms=self.check_transforms,
run_id=2,
**kwargs
)
# Compute the TL of check_transform
if self.check_transforms:
for transform in obsoper.data_tl:
data_tl_ref0 = data_tl_save[transform]
data_tl_prime0 = obsoper.data_tl[transform]
for ddi in data_tl_prime0:
data_tl_ref1 = data_tl_ref0[ddi]
data_tl_prime1 = data_tl_prime0[ddi]
for trid in data_tl_prime1["outputs"]:
data_tl_ref2 = data_tl_ref1["outputs"][trid]
data_tl_prime2 = data_tl_prime1["outputs"][trid]
for ddi_trid in data_tl_prime2:
data_tl_ref3 = data_tl_ref2[ddi_trid]
data_tl_prime3 = data_tl_prime2[ddi_trid]
for transf_trid in data_tl_prime3:
if "spec" not in data_tl_ref3[transf_trid]:
continue
if data_tl_ref3[transf_trid]["spec"].dtype != "O":
data_tl_prime3[transf_trid]["incr"] = \
data_tl_prime3[transf_trid]["spec"] \
- data_tl_ref3[transf_trid]["spec"]
for trid in data_tl_prime1["inputs"]:
data_tl_ref2 = data_tl_ref1["inputs"][trid]
data_tl_prime2 = data_tl_prime1["inputs"][trid]
for ddi_trid in data_tl_prime2:
data_tl_ref3 = data_tl_ref2[ddi_trid]
data_tl_prime3 = data_tl_prime2[ddi_trid]
for transf_trid in data_tl_prime3:
if "spec" not in data_tl_ref3[transf_trid]:
continue
if data_tl_ref3[transf_trid]["spec"].dtype != "O":
data_tl_prime3[transf_trid]["incr"] = \
data_tl_prime3[transf_trid]["spec"] \
- data_tl_ref3[transf_trid]["spec"]
# Replace dy by difference between ref and perturbed control vector
obsvect.dy = obsvect_dy.ysim - ysim_ref
# Putting xb and x back to there original values
controlvect.x -= controlvect.dx_save
# Computing < H.dx, H.dx >
scaleprod1 = np.nansum(np.power(obsvect.dy, 2))
# Putting increments in the observation vector
obsvect.dy = np.where(np.isnan(obsvect.dy), 0, obsvect.dy)
# Running the observation operator
controlvect = obsoper.obsoper(
controlvect,
obsvect,
"adj",
datei=datei,
datef=datef,
workdir=workdir,
reload_results=self.reload_results,
check_transforms=self.check_transforms,
**kwargs
)
# Computing < dx, H*(H.dx) >
if testspace == "control":
scaleprod2 = np.nansum(controlvect.dx_save * controlvect.dx)
elif testspace == "chi":
scaleprod2 = np.nansum(
controlvect.sqrtbprod_ad(
controlvect.dx, **kwargs) * controlvect.chi
)
# Final verbose
info('Machine accuracy: {}'.format(accuracy))
info('< H.dx, H.dx > = {:.17E}'.format(scaleprod1))
info('< dx, H*(H.dx) > = {:.17E}'.format(scaleprod2))
info('The difference is {:.1E} times the machine accuracy'
.format(np.abs(scaleprod2 / scaleprod1 - 1) / accuracy))
# Return test
return np.floor(np.abs(scaleprod2 / scaleprod1 - 1) / accuracy)
