Retrieval 4: Using parameter functions
About this retrieval example
This example shows how to use the rt1 python package together with scipy optimize to setup a retrieval procedure to:
obtain static and dynamic parameters from a series of incidence-angle dependent \(\sigma^0\) measurements.
use functions to represent model parameters
Show code cell source
%matplotlib widget
from rt1_model import RT1, surface, volume, set_loglevel
from scipy.optimize import least_squares
import matplotlib.pyplot as plt
import numpy as np
rand = np.random.RandomState(123456) # initialize a reproducible random state
set_loglevel("info")
Specify simulation and fit parameters
Set parameter values that are used to simulate the data
dB, sig0 = False, True
num = 100 # Number of measurements
incs = 30 # Available incidence angles per measurement
noise_sigma = 0.5 if dB is True else 1e-3 # Noise-level (sigma of gaussian noise)
inc = rand.normal(45, 10, (num, incs)).clip(20, 70) # Incidence angles
N = rand.normal(0.1, 0.1, (num, 1)).clip(0.01, 0.25) # NormBRDF values
phi = np.linspace(0, 2.* np.pi, num)[:,np.newaxis]
dphi = rand.randint(0, 200)/100
sim_params = dict(k=2.5, dphi=dphi, omega=0.2, N=N) # Simulation parameter values
const_params = dict(t_s=0.4, phi=phi, PI=np.pi) # Constant parameters (assumed to be known)
Set start values and boundaries for the fit
start_vals = dict(omega=[0.1], k=[1], N=[0.1]*num, dphi=[1])
bnd_vals = dict(omega=[(0.01, 0.5)], k=[(0.5, 5.)], N=[(0.01, 0.5)]*num, dphi=[(0., np.pi)])
Setup RT1 and create a simulated dataset
V = volume.Rayleigh()
SRF = surface.HG_nadirnorm(t="t_s", ncoefs=10)
R = RT1(V=V, SRF=SRF, int_Q=True, dB=dB, sig0=sig0)
R.set_monostatic(p_0=0)
R.NormBRDF = "N/3 + 5*N**2"
R.tau = "k * (phi + dphi) / (6*PI) * sin(phi + dphi)**(2)"
R.set_geometry(t_0=np.deg2rad(inc))
tot = R.calc(**sim_params, **const_params)[0]
tot += rand.normal(0, noise_sigma, tot.shape) # Add some random noise
Show code cell output
09:13:59.114 INFO: Evaluating coefficients for interaction-term...
09:13:59.253 INFO: Coefficients extracted, it took 0.02107 sec.
Show code cell source
f, (ax, ax2) = plt.subplots(2, figsize=(10, 3.5), sharex=True, gridspec_kw=dict(hspace=0.15))
f.canvas.header_visible = False
f.text(0.5, .99, "Resulting input timeseries for 'NormBRDF' and 'tau' parameters functions:", ha="center", va="top")
f.text(.5, .92, f"NormBRDF = {R._NormBRDF} tau = {R._tau} || dphi={dphi}", ha="center", va="top", fontsize=8)
ax.plot(R.NormBRDF, marker=".", lw=0.5, label="NormBRDF")
ax.set_ylabel("NormBRDF")
ax.legend(loc="upper left")
ax_t = ax.twinx()
ax_t.plot(N, marker=".", lw=0.5, label="N", c="C1")
ax_t.set_xlabel("# measurement")
ax_t.set_ylabel("N")
ax_t.legend(loc="upper left", bbox_to_anchor=(0, .8))
ax2.plot(R.tau, marker=".", lw=0.5, label="tau")
ax2.set_ylabel("tau")
ax2.legend(loc="upper left")
ax2_t = ax2.twinx()
ax2_t.plot(phi, marker=".", lw=0.5, label="phi", c="C1")
ax2_t.set_ylabel("phi")
ax2_t.legend(loc="upper left", bbox_to_anchor=(0, .8))
f.subplots_adjust(left=0.075, right=0.9)
Setup scipy optimize to fit RT1 model to the data
def parse_params(x):
"""Map 1D parameter array to dict {parameter_name: value(s)}."""
return dict(omega=x[0], k=x[1], dphi=x[2], N=x[3:][:, np.newaxis])
def fun(x, **kwargs):
"""Calculate residuals."""
R.update_params(**parse_params(x), **kwargs)
res = (R.calc()[0] - tot).ravel() # Ravel output because scipy requires 1D arrays
return res
def jac(x, **kwargs):
"""Calculate jacobian."""
R.update_params(**parse_params(x), **kwargs)
jac = R.jacobian(param_list=["omega", "k", "dphi", "N"], format="scipy_least_squares")
return jac
# Unpack start-values and boundaries as required by scipy optimize
x0 = [*start_vals["omega"], *start_vals["k"], *start_vals["dphi"], *start_vals["N"]]
bounds = list(zip(*[*bnd_vals["omega"], *bnd_vals["k"], *bnd_vals["dphi"], *bnd_vals["N"]]))
res = least_squares(
fun=fun,
x0=x0,
bounds=bounds,
jac=jac,
#x_scale="jac",
ftol=1e-5,
gtol=1e-5,
xtol=1e-5,
verbose=2,
kwargs=const_params, # pass constant parameters to "fun" and "jac"
)
# Unpack found parameters
found_params = parse_params(res.x)
# Calcuate total backscatter based on found parameters
found_tot = R.calc(**found_params, **const_params)[0]
Show code cell output
Iteration Total nfev Cost Cost reduction Step norm Optimality
0 1 1.8725e+01 9.91e+00
1 2 9.8003e-01 1.77e+01 7.67e-01 1.85e+00
2 4 6.9148e-01 2.89e-01 1.49e-01 2.47e+00
3 5 6.5098e-01 4.05e-02 3.15e-01 2.52e+00
4 6 3.5180e-01 2.99e-01 6.00e-02 1.77e+00
5 7 1.8636e-01 1.65e-01 1.34e-01 9.68e-01
6 8 8.1120e-02 1.05e-01 8.61e-02 3.66e-01
7 9 3.0622e-02 5.05e-02 9.21e-02 1.61e-01
8 11 2.9938e-02 6.84e-04 1.12e-01 2.34e-01
9 12 2.0934e-02 9.00e-03 1.98e-02 4.64e-02
10 14 2.0335e-02 5.99e-04 1.66e-02 5.60e-02
11 15 1.9267e-02 1.07e-03 2.94e-02 8.44e-02
12 16 1.7490e-02 1.78e-03 2.73e-02 7.31e-02
13 17 1.7447e-02 4.28e-05 6.57e-02 1.32e-01
14 18 1.4393e-02 3.05e-03 1.05e-02 4.57e-02
15 19 1.3480e-02 9.13e-04 3.17e-02 6.68e-02
16 20 1.2474e-02 1.01e-03 3.18e-02 6.64e-02
17 21 1.1467e-02 1.01e-03 3.23e-02 6.56e-02
18 22 1.0558e-02 9.09e-04 3.29e-02 6.46e-02
19 23 9.7179e-03 8.40e-04 3.32e-02 6.31e-02
20 24 8.9266e-03 7.91e-04 3.33e-02 6.09e-02
21 25 8.2019e-03 7.25e-04 3.33e-02 5.85e-02
22 26 7.5243e-03 6.78e-04 3.32e-02 5.91e-02
23 27 6.9017e-03 6.23e-04 3.29e-02 6.04e-02
24 28 6.3252e-03 5.76e-04 3.25e-02 6.14e-02
25 29 5.7961e-03 5.29e-04 3.20e-02 6.19e-02
26 30 5.3100e-03 4.86e-04 3.14e-02 6.22e-02
27 31 4.8662e-03 4.44e-04 3.08e-02 6.21e-02
28 32 4.4605e-03 4.06e-04 3.00e-02 6.18e-02
29 33 4.0928e-03 3.68e-04 2.92e-02 6.11e-02
30 34 3.7596e-03 3.33e-04 2.83e-02 6.02e-02
31 35 3.4590e-03 3.01e-04 2.73e-02 5.90e-02
32 36 3.1891e-03 2.70e-04 2.63e-02 5.75e-02
33 37 2.9500e-03 2.39e-04 2.52e-02 5.59e-02
34 38 2.7313e-03 2.19e-04 2.40e-02 5.38e-02
35 39 2.5437e-03 1.88e-04 2.28e-02 5.15e-02
36 40 2.3750e-03 1.69e-04 2.16e-02 4.91e-02
37 41 2.2291e-03 1.46e-04 2.04e-02 4.65e-02
38 42 2.1005e-03 1.29e-04 1.91e-02 4.39e-02
39 43 1.9906e-03 1.10e-04 1.79e-02 4.13e-02
40 44 1.8930e-03 9.76e-05 1.66e-02 3.84e-02
41 45 1.8115e-03 8.16e-05 1.53e-02 3.55e-02
42 46 1.7413e-03 7.02e-05 1.40e-02 3.26e-02
43 47 1.6825e-03 5.88e-05 1.28e-02 2.98e-02
44 48 1.6328e-03 4.97e-05 1.17e-02 2.70e-02
45 49 1.5917e-03 4.12e-05 1.06e-02 2.44e-02
46 50 1.5575e-03 3.42e-05 9.54e-03 2.19e-02
47 51 1.5299e-03 2.76e-05 8.62e-03 1.97e-02
48 52 1.5066e-03 2.33e-05 7.69e-03 1.74e-02
49 53 1.4878e-03 1.87e-05 6.78e-03 1.52e-02
50 54 1.4738e-03 1.41e-05 6.03e-03 1.33e-02
51 55 1.4614e-03 1.24e-05 5.28e-03 1.15e-02
52 56 1.4528e-03 8.67e-06 4.59e-03 9.80e-03
53 57 1.4453e-03 7.47e-06 3.91e-03 8.01e-03
54 58 1.4403e-03 4.96e-06 3.33e-03 6.50e-03
55 59 1.4362e-03 4.13e-06 2.73e-03 4.87e-03
56 60 1.4358e-03 4.25e-07 1.65e-03 6.14e-03
57 61 1.4329e-03 2.90e-06 5.92e-04 6.63e-04
58 62 1.4326e-03 2.82e-07 3.59e-04 8.07e-04
59 63 1.4322e-03 3.70e-07 7.25e-04 1.19e-03
60 64 1.4319e-03 3.51e-07 8.36e-04 1.57e-03
61 65 1.4315e-03 4.07e-07 6.42e-04 1.35e-03
62 66 1.4312e-03 2.29e-07 3.37e-04 1.24e-03
63 67 1.4310e-03 1.77e-07 2.69e-04 1.11e-03
64 68 1.4309e-03 1.35e-07 2.29e-04 1.02e-03
65 69 1.4308e-03 1.14e-07 2.03e-04 9.47e-04
66 70 1.4307e-03 9.95e-08 1.84e-04 8.82e-04
67 71 1.4306e-03 8.96e-08 1.67e-04 8.19e-04
68 72 1.4305e-03 8.00e-08 1.55e-04 7.63e-04
69 73 1.4305e-03 7.22e-08 1.43e-04 7.10e-04
70 74 1.4304e-03 6.38e-08 1.30e-04 6.63e-04
71 75 1.4303e-03 5.73e-08 1.21e-04 6.20e-04
72 76 1.4303e-03 5.14e-08 1.12e-04 5.81e-04
73 77 1.4302e-03 4.64e-08 1.04e-04 5.46e-04
74 78 1.4302e-03 4.21e-08 9.72e-05 5.13e-04
75 79 1.4302e-03 3.83e-08 9.10e-05 4.83e-04
76 80 1.4301e-03 3.49e-08 8.53e-05 4.55e-04
77 81 1.4301e-03 3.19e-08 8.02e-05 4.29e-04
78 82 1.4301e-03 2.92e-08 7.55e-05 4.05e-04
79 83 1.4300e-03 2.68e-08 7.08e-05 3.83e-04
80 84 1.4300e-03 2.47e-08 6.69e-05 3.63e-04
81 85 1.4300e-03 2.28e-08 6.33e-05 3.43e-04
82 86 1.4300e-03 2.11e-08 5.97e-05 3.25e-04
83 87 1.4299e-03 1.96e-08 5.67e-05 3.09e-04
84 88 1.4299e-03 1.82e-08 5.39e-05 2.93e-04
85 89 1.4299e-03 1.70e-08 5.13e-05 2.78e-04
86 90 1.4299e-03 1.59e-08 4.89e-05 2.64e-04
87 91 1.4299e-03 1.52e-08 4.74e-05 2.51e-04
88 92 1.4299e-03 1.40e-08 4.52e-05 2.39e-04
`ftol` termination condition is satisfied.
Function evaluations 92, initial cost 1.8725e+01, final cost 1.4299e-03, first-order optimality 2.39e-04.
| Parameter | Target value | Start value | Retrieved value | (Target - Retrieved) |
|---|---|---|---|---|
| omega | 0.200 | 0.100 | 0.200 | -0.000 |
| k | 2.500 | 1.000 | 2.505 | 0.005 |
| dphi | 1.740 | 1.000 | 1.740 | 0.000 |
| N (mean) | 0.108 | 0.100 | 0.108 | 0.000 |
Visualize Results
Plot timeseries
Show code cell source
f, (ax, ax2) = plt.subplots(2, figsize=(10, 4), sharex=True)
f.canvas.header_visible = False
# Plot retrieved parameter timeseries
ax.set_ylabel("N / tau")
ax.plot(sim_params["N"], marker=".", lw=0.25, label="target N", c="C0")
ax.plot(found_params["N"], marker="o", lw=0.25, markerfacecolor="none", label="retrieved N", c="C0")
R.update_params(**sim_params, **const_params)
ax.plot(R.tau, marker=".", lw=0.25, label="target tau", c="C1")
R.update_params(**found_params, **const_params)
ax.plot(R.tau, marker="o", lw=0.25, markerfacecolor="none", label="retrieved tau", c="C1")
# Plot backscatter timeseries
ax2.set_ylabel(r"$\sigma_0$ [dB]")
ax2.set_xlabel("# measurement")
ax2.plot(tot, lw=0, marker=".", c="C0", ms=3)
ax2.plot(found_tot, lw=0, marker="o", markerfacecolor="none", c="C1", ms=3)
ax.legend(loc="upper center", ncols=2, bbox_to_anchor=(0.5, 1.5))
f.tight_layout()
Initialize analyzer widget and overlay results
Show code cell source
analyze_params = {key: (*np.mean(np.atleast_2d(bnd_vals[key]), axis=0), found_params[key].mean()) for key in found_params}
# pick the measurement for which you want to analyze the components
R.update_params(phi=0)
ana = R.analyze(**analyze_params, range_parameter="k")
# Plot fit-data on top
ana.ax.scatter(inc, tot, c="k", s=3, zorder=0)
ana.ax.scatter(inc, found_tot, c="C0", s=1, zorder=0)
# Indicate fit-results in slider-axes
for key, s in ana.sliders.items():
if key in ["omega", "k", "dphi"]:
s.ax.plot(sim_params[key], np.mean(s.ax.get_ylim()), marker="o")
# Add text for static parameters
t = ana.f.text(
0.6,
0.95,
"\n".join(
[
f"{key:>8} = {found_params[key]:.3f} ({sim_params[key]:.2f}) "
rf"| $\Delta$ = {found_params[key] - sim_params[key]: .3f}"
for key in ["omega", "k", "dphi"]
]
),
va="top",
fontdict=dict(family="monospace", size=8),
)
09:14:04.672 INFO: Evaluating coefficients for interaction-term...
09:14:04.743 INFO: Coefficients extracted, it took 0.02186 sec.