import numpy as np
import copy
from logging import info
[docs]
def testing_TL_linearity(self, accuracy):
"""Verify that the tangent-linear operator is linear: ``TL(λ·dx) == λ·TL(dx)``.
Runs the TL operator twice — once for ``dx`` and once for ``λ·dx`` —
and prints the ratio of the two scalar products. The ratio should equal 1
to machine precision for a truly linear operator.
Args:
self (Plugin): mode plugin with all configuration attributes.
accuracy (float): machine epsilon (``np.finfo(np.float64).eps``).
"""
# 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
# Scaling factor
lambdalin = getattr(self, 'lambdalin', 0.1)
# Running the tangent linear code of the model at x for incr
info("Testing the linearity of the tangent-linear")
obsvect = obsoper.obsoper(
controlvect, obsvect, 'tl',
datei=datei, datef=datef,
workdir=workdir,
reload_results=getattr(self, 'reload_results', False),
run_id=0)
# Computing TL_x(incr) * lambda
scaleprod1 = lambdalin * np.nansum(obsvect.dy)
# Running the tangent linear code of the model at x for lambda * incr
controlvect.dx *= lambdalin
obsvect = obsoper.obsoper(controlvect, obsvect, 'tl',
datei=datei, datef=datef,
workdir=workdir,
reload_results=getattr(self, 'reload_results',
False),
run_id=1)
# Computing TL_x(lambda * incr)
scaleprod2 = np.nansum(obsvect.dy)
# Final verbose
info('Machine accuracy: {}'.format(accuracy))
info('lambda TL_H(x)(dx) = {:.17E}'.format(scaleprod1))
info('TL_H(x)(lambda.dx), = {:.17E}'.format(scaleprod2))
info(
"Result of TL linearity test: {:.17E} = 1 for lambda = {:.17E} ?".format(
scaleprod1 / scaleprod2, lambdalin))
info('The difference is {:.1E} times the machine accuracy'
.format(np.abs(scaleprod1 - scaleprod2) / accuracy / scaleprod2))
[docs]
def testing_TL_Taylor(self, accuracy):
"""Run a Taylor-expansion test of the tangent-linear operator.
Checks the first-order convergence property:
``‖H(x + λ·dx) − H(x) − TL(λ·dx)‖ = O(λ²)``.
Iterates ``taylor_iter`` times, halving ``λ`` by factor ``lambdalin`` each
time, and prints a summary table of norms and ratios.
Args:
self (Plugin): mode plugin with all configuration attributes.
accuracy (float): machine epsilon (``np.finfo(np.float64).eps``).
"""
# 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
# Computes forward model with basic inputs
obsvect = obsoper.obsoper(controlvect, obsvect, 'fwd',
datei=datei, datef=datef,
workdir=workdir,
reload_results=getattr(self, 'reload_results',
False),
run_id=0)
scaleprod0 = copy.deepcopy(obsvect.ysim)
info(f'H(x) = {np.nanmean(scaleprod0):.17E} +/- {np.nanstd(scaleprod0):.17E}')
# number of Taylor's loops:
nloops = getattr(self, 'taylor_iter', 1)
# Scaling factor
lambdalin = getattr(self, 'lambdalin', 0.1)
# Loops on smaller and smaller incr
lambdataylor = lambdalin
x0 = copy.deepcopy(controlvect.x)
dx0 = copy.deepcopy(controlvect.dx)
ratios = []
differences = []
list_scaleprod1 = []
list_scaleprod2 = []
list_lambdataylor = []
for i in range(nloops):
list_lambdataylor.append(lambdataylor)
# Computes perturbed forward model at x + lambda * dx
controlvect.x = x0 + lambdataylor * dx0
obsvect = obsoper.obsoper(
controlvect, obsvect, 'fwd',
datei=datei, datef=datef,
workdir=workdir,
reload_results=getattr(self, 'reload_results', False),
run_id=i + 1)
scaleprod1 = copy.deepcopy(obsvect.ysim)
list_scaleprod1.append(scaleprod1)
info(f'H(x+lambda*incr) = {np.nanmean(scaleprod1):.17E} +/- {np.nanstd(scaleprod1):.17E}')
# Computes TL at x for lambda * dx
controlvect.x = x0
controlvect.dx = lambdataylor * dx0
obsvect = obsoper.obsoper(
controlvect, obsvect, 'tl',
datei=datei, datef=datef,
workdir=workdir,
reload_results=getattr(self, 'reload_results', False),
run_id=i + 100)
scaleprod2 = copy.deepcopy(obsvect.dy)
list_scaleprod2.append(scaleprod2)
info(f'TL_H(x)(lambda.dx) = {np.nanmean(scaleprod2):.17E} +/- {np.nanstd(scaleprod2):.17E}')
ratios.append((scaleprod1 - scaleprod0) / scaleprod2)
differences.append(np.abs((scaleprod1 - scaleprod0) - scaleprod2))
# info('Ratio = {:.17E} for lambda = {:.17E}'.format(ratios[i],
# lambdataylor))
lambdataylor *= lambdalin
info("Summary:")
info("H(x) = {:.17E} ".format(np.nansum(scaleprod0 ** 2) ** 0.5))
info("lambda , || TL_H(x)(lambda.dx) || , < H(x+lambda*incr) - H(x) | TL_H(x)(lambda.dx) > , (1 - X) / lambda ")
info("\n".join(
[f"{lambdataylor:.17E} , "
f"{np.nansum(scaleprod2 ** 2):.17E} "
f"{np.nansum((scaleprod1 - scaleprod0) * scaleprod2):.17E} "
f"{1 - (np.nansum((scaleprod1 - scaleprod0) * scaleprod2) / np.nansum(scaleprod2 ** 2)):.17E}"
for scaleprod1, scaleprod2, lambdataylor
in zip(list_scaleprod1, list_scaleprod2, list_lambdataylor)
]
)
)
print(__file__)
import code
code.interact(local=dict(locals(), **globals()))
