import os
from logging import info
import itertools
import numpy as np
import pandas as pd
[docs]
def check_adjtltest(self, period_order, mapper, transform_pipe):
"""Run the adjoint / tangent-linear dot-product test on each transform.
For every ``(date, transform)`` pair processed in the forward
direction, computes the two sides of the dot-product identity used to
validate the adjoint of a transform::
<dx | H^T(H(dx))> vs <H(dx) | H(dx)>
``dx_in`` (left-hand side) is accumulated from the transform's inputs,
by summing the element-wise product of the tangent-linear increment
(``incr``) and the adjoint sensitivity (``delta_adj_out`` if present,
else ``adj_out``) of each precursor. ``dx_out`` (right-hand side) is
accumulated the same way from the transform's outputs and successors.
The two quantities should be equal, up to floating-point round-off, if
the transform's adjoint is correctly implemented; the relative
difference is reported in units of machine epsilon.
Results for all checked transforms are logged and written to
``{self.workdir}/check_transforms.log``.
Args:
self: obsoperator instance exposing ``data_tl`` and ``data_adj``
(datastores keyed by transform name, then by sub-simulation
date/id ``ddi``, holding the tangent-linear increments and
adjoint outputs for each input/output tracer ``trid``) as well
as ``workdir``.
period_order (list[tuple]): ordered ``(ddi, transform, direction)``
triplets describing the transforms scheduled to run, as built
by the transform pipeline. Only ``direction == "forward"``
entries are checked here (the tangent-linear/adjoint pair for
a transform is compared once, not once per direction).
mapper (dict): per-transform metadata (``inputs``, ``outputs``,
``precursors``, ``successors``, ``subsimus``) describing how
sub-simulations and tracers (``trid``) are linked across
transforms.
transform_pipe: container exposing each transform instance as an
attribute (``getattr(transform_pipe, transform)``); used here
only to retrieve the plugin name of transforms and their
precursors/successors.
Returns:
None. Results are emitted via the ``info`` logger and written to
``check_transforms.log`` in ``self.workdir``.
"""
df = None
# Check the adjoint for each transformations
for ddi, transform, direction in period_order:
if transform not in self.data_adj \
or transform not in self.data_tl \
or direction != "forward":
continue
data_tl = self.data_tl[transform]
data_adj = self.data_adj[transform]
dx_in = 0
dx_out = 0
if ddi not in data_tl:
continue
# dx_in = <dx | H^T(H(dx))>, accumulated over the transform's
# inputs and their precursors
list_trids = data_tl.get(ddi, {"inputs": {}})["inputs"]
for trid in list_trids:
# Skip if data not in adj or tl datastore
if ddi not in data_adj:
continue
if data_tl[ddi]["inputs"].get(trid, {}) == {} \
and data_adj[ddi]["inputs"].get(trid, {}) == {}:
continue
# Force using the transform input itself if no precursor
auto_precursor = False
if trid not in mapper[transform]["precursors"]:
precursors = [transform]
auto_precursor = False
list_di = mapper[transform]["subsimus"][ddi]["outputs"][trid]
else:
precursors = mapper[transform]["precursors"][trid]
list_di = mapper[transform]["subsimus"][ddi]["inputs"][trid]
for di, precursor in itertools.product(list_di, precursors):
# Get value from dump2input transform if force_dump
precursor_name = getattr(
transform_pipe, precursor).plugin.name
if precursor_name == "dump2inputs":
if precursor in self.data_tl \
and precursor in self.data_adj:
tl = self.data_tl[precursor][di][
"outputs"][trid][di][transform]
adj = self.data_adj[precursor][di][
"outputs"][trid][di][transform]
else:
continue
else:
tl = data_tl[ddi]["inputs"][trid][di].get(
precursor, {ddi: data_tl[ddi]["inputs"][trid][di]})
adj = data_adj[ddi]["inputs"][trid][di].get(
precursor, {ddi: data_adj[ddi]["inputs"][trid][di]})
# Loop over precursor sub-simulations
subsimus_in = mapper[precursor]["subsimus"]
for precursor_di in subsimus_in:
# Skip if subsimu not in precursor subsimulations
if di not in subsimus_in[precursor_di]['outputs'].get(trid, []):
continue
tl_precursor = tl[precursor_di]
adj_precursor = adj[precursor_di]
# Deal with sparse data
if isinstance(tl_precursor, pd.DataFrame) \
and isinstance(adj_precursor, pd.DataFrame):
tl_precursor = tl_precursor["maindata"]
adj_precursor = adj_precursor["maindata"]
if "incr" not in tl_precursor or "adj_out" not in adj_precursor:
continue
data_tl_in = tl_precursor["incr"]
data_adj_in = adj_precursor.get(
"delta_adj_out", adj_precursor["adj_out"])
dx_in += np.float64((data_tl_in * data_adj_in).sum())
# dx_out = <H(dx) | H(dx)>, accumulated over the transform's
# outputs and their successors
for trid in mapper[transform]["outputs"]:
if ddi not in data_tl or ddi not in data_adj:
continue
if trid not in data_tl[ddi]["outputs"] \
or trid not in data_adj[ddi]["outputs"]:
continue
for di in mapper[transform]["subsimus"][ddi]["outputs"][trid]:
for successor in mapper[transform]["successors"][trid]:
# Get value from loadfromoutputs transform if force_loadin
successor_name = getattr(
transform_pipe, successor).plugin.name
if successor_name == "loadfromoutputs" \
and successor in self.data_tl \
and successor in self.data_adj:
if di not in self.data_tl[successor] \
or di not in self.data_adj[successor]:
continue
tl = self.data_tl[successor][di]["inputs"][trid][di][
transform]
adj = self.data_adj[successor][di]["inputs"][trid][di][
transform]
else:
tl = data_tl[ddi]["outputs"][trid][di][successor]
adj = data_adj[ddi]["outputs"][trid][di][successor]
# Loop over successor sub-simulations
subsimus_in = mapper[successor]["subsimus"]
for successor_di in subsimus_in:
# Skip if subsimu not in successor subsimulations
if di not in subsimus_in[successor_di]['inputs'].get(trid, []):
continue
tl_successor = tl[successor_di]
adj_successor = adj[successor_di]
# Deal with sparse data
if isinstance(tl_successor, pd.DataFrame) \
and isinstance(adj_successor, pd.DataFrame):
tl_successor = tl_successor["maindata"]
adj_successor = adj_successor["maindata"]
if "incr" not in tl_successor or "adj_out" not in adj_successor:
continue
data_tl_out = tl_successor["incr"]
data_adj_out = adj_successor["adj_out"]
dx_out += np.float64(
(data_tl_out * data_adj_out).sum())
# Duplicate input / output for toobsvect and fromcontrol
transf_name = getattr(transform_pipe, transform).plugin.name
if transf_name == "toobsvect":
dx_out = dx_in
# elif transf_name == "fromcontrol":
# print(__file__)
# import code
# code.interact(local=dict(locals(), **globals()))
# dx_in = dx_out
# Skip transforms for which neither side of the test could be
# computed for this sub-simulation (nothing to report)
if dx_in == dx_out == 0:
continue
# Relative difference between the two sides of the dot-product
# test, expressed in units of machine epsilon
diff = (
(dx_in / dx_out - 1) / np.finfo(np.float64).eps
if dx_out != 0 else np.nan
)
df_tr = pd.DataFrame(
data=[[transform, transf_name, ddi, dx_in, dx_out, diff]],
columns=["id", "name", "date", "dx_in", "dx_out", "diff"]
)
df = df_tr if df is None else pd.concat([df, df_tr], ignore_index=True)
# Format the results table and write it to the log / logfile
df["date"] = df["date"].map(lambda x: x.isoformat(sep=" "))
df["dx_in"] = df["dx_in"].map(lambda x: f"{x:.15f}")
df["dx_out"] = df["dx_out"].map(lambda x: f"{x:.15f}")
df["diff"] = df["diff"].map(lambda x: f"{x:.1E}")
df.columns = ["transform_ID", "transform_name", "sub-simulation",
"<dx|H*(H(dx))>", "<H(dx)|H(dx)>", "difference"]
df_str = df.to_string(index=False)
# Write to logfile
info("check transforms results:\n" + df_str)
adjtl_log = os.path.join(self.workdir, "check_transforms.log")
with open(adjtl_log, "w") as f:
f.write(df_str)
