import numpy as np
import copy
from logging import info
def _apply_finite_difference_tl(data_tl_ref, data_tl_prime):
"""Turn two saved forward runs into a finite-difference TL increment.
For every transform / sub-simulation / trid / date recorded in
``data_tl_prime`` (the perturbed run ``H(xb + dx)``), computes
``incr = spec_prime - spec_ref`` against the matching entry in
``data_tl_ref`` (the reference run ``H(xb)``) and stores it back into
``data_tl_prime`` in place, for both its ``"outputs"`` and
``"inputs"`` sides. Entries without a ``"spec"`` field, or whose
``"spec"`` is an object dtype, are left untouched.
Args:
data_tl_ref: ``obsoper.data_tl`` saved from the reference run.
data_tl_prime: ``obsoper.data_tl`` from the perturbed run, mutated
in place with the computed ``"incr"`` entries.
"""
for transform in data_tl_prime:
ref0 = data_tl_ref[transform]
prime0 = data_tl_prime[transform]
for ddi in prime0:
ref1 = ref0[ddi]
prime1 = prime0[ddi]
for inout in ("outputs", "inputs"):
for trid in prime1[inout]:
ref2 = ref1[inout][trid]
prime2 = prime1[inout][trid]
for ddi_trid in prime2:
ref3 = ref2[ddi_trid]
prime3 = prime2[ddi_trid]
for transf_trid in prime3:
if "spec" not in ref3[transf_trid]:
continue
if ref3[transf_trid]["spec"].dtype != "O":
prime3[transf_trid]["incr"] = (
prime3[transf_trid]["spec"]
- ref3[transf_trid]["spec"]
)
[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:
_apply_finite_difference_tl(data_tl_save, obsoper.data_tl)
# 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(f'Machine accuracy: {accuracy}')
info(f'< H.dx, H.dx > = {scaleprod1:.17E}')
info(f'< dx, H*(H.dx) > = {scaleprod2:.17E}')
info(f'The difference is {np.abs(scaleprod2 / scaleprod1 - 1) / accuracy:.1E} times the machine accuracy')
# Return test
return np.floor(np.abs(scaleprod2 / scaleprod1 - 1) / accuracy)
