Dual remote sensing

In this notebook we consider an application to dual remote sensing: a univariate field is remotely-sensed, with each measurement having a white noise error of a certain level. We therefore jointly estimate the noise level for each sensor in parallel with the parameters of the covariance model of the sensed field. We then proceed to "map" the field using both measurements.

Imports¶

In [1]:

import matplotlib.pyplot as plt

import matplotlib.pyplot as plt

In [2]:

from debiased_spatial_whittle.models import DualRemoteSensing
from debiased_spatial_whittle.models.univariate import SquaredExponentialModel
from debiased_spatial_whittle.grids import RectangularGrid
from debiased_spatial_whittle.sampling import MultivariateSamplerOnRectangularGrid

from debiased_spatial_whittle.models import DualRemoteSensing from debiased_spatial_whittle.models.univariate import SquaredExponentialModel from debiased_spatial_whittle.grids import RectangularGrid from debiased_spatial_whittle.sampling import MultivariateSamplerOnRectangularGrid

In [3]:

from debiased_spatial_whittle.inference import Periodogram, ExpectedPeriodogram, MultivariateDebiasedWhittle, Estimator

from debiased_spatial_whittle.inference import Periodogram, ExpectedPeriodogram, MultivariateDebiasedWhittle, Estimator

Grid and model specification¶

In [4]:

grid = RectangularGrid((64, 64), nvars=2)
latent_model = SquaredExponentialModel(rho=5.0)
model = DualRemoteSensing(latent_model, sigma1=0.1, sigma2=0.2)

grid = RectangularGrid((64, 64), nvars=2) latent_model = SquaredExponentialModel(rho=5.0) model = DualRemoteSensing(latent_model, sigma1=0.1, sigma2=0.2)

Sample a realization¶

In [5]:

sampler = MultivariateSamplerOnRectangularGrid(model, grid, p=2)
s = sampler()

sampler = MultivariateSamplerOnRectangularGrid(model, grid, p=2) s = sampler()

In [6]:

plt.figure()
plt.subplot(1, 2, 1)
plt.imshow(s[..., 0], cmap="seismic")
plt.subplot(1, 2, 2)
plt.imshow(s[..., 1], cmap="seismic")
plt.show()

plt.figure() plt.subplot(1, 2, 1) plt.imshow(s[..., 0], cmap="seismic") plt.subplot(1, 2, 2) plt.imshow(s[..., 1], cmap="seismic") plt.show()

Parameter inference¶

In [7]:

latent_model = SquaredExponentialModel(rho=1.0)
model = DualRemoteSensing(latent_model, sigma1=0.01, sigma2=0.01)
model.set_param_bounds(dict(sigma1=(0.001, 1.0), sigma2=(0.001, 1.0)))

latent_model = SquaredExponentialModel(rho=1.0) model = DualRemoteSensing(latent_model, sigma1=0.01, sigma2=0.01) model.set_param_bounds(dict(sigma1=(0.001, 1.0), sigma2=(0.001, 1.0)))

In [8]:

periodogram = Periodogram()
dbw = MultivariateDebiasedWhittle(periodogram, ExpectedPeriodogram(grid, periodogram))
estimator = Estimator(dbw)
estimator(model, s, opt_callback=lambda *args, **kwargs: print(*args, **kwargs))

periodogram = Periodogram() dbw = MultivariateDebiasedWhittle(periodogram, ExpectedPeriodogram(grid, periodogram)) estimator = Estimator(dbw) estimator(model, s, opt_callback=lambda *args, **kwargs: print(*args, **kwargs))

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[8], line 4
      2 dbw = MultivariateDebiasedWhittle(periodogram, ExpectedPeriodogram(grid, periodogram))
      3 estimator = Estimator(dbw)
----> 4 estimator(model, s, opt_callback=lambda *args, **kwargs: print(*args, **kwargs))

File ~/checkouts/readthedocs.org/user_builds/debiased-spatial-whittle/checkouts/latest/src/debiased_spatial_whittle/inference/likelihood.py:682, in Estimator.__call__(self, model, sample, opt_callback)
    672         opt_result = minimize(
    673             opt_func,
    674             x0,
   (...)
    679             options=self.optim_options,
    680         )
    681     else:
--> 682         opt_result = minimize(
    683             opt_func,
    684             x0,
    685             method=self.method,
    686             bounds=bounds,
    687             callback=opt_callback,
    688             options=self.optim_options,
    689         )
    690 model.update_free_parameters(opt_result.x)
    691 self.opt_result = opt_result

File ~/checkouts/readthedocs.org/user_builds/debiased-spatial-whittle/envs/latest/lib/python3.10/site-packages/scipy/optimize/_minimize.py:713, in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
    710     res = _minimize_newtoncg(fun, x0, args, jac, hess, hessp, callback,
    711                              **options)
    712 elif meth == 'l-bfgs-b':
--> 713     res = _minimize_lbfgsb(fun, x0, args, jac, bounds,
    714                            callback=callback, **options)
    715 elif meth == 'tnc':
    716     res = _minimize_tnc(fun, x0, args, jac, bounds, callback=callback,
    717                         **options)

File ~/checkouts/readthedocs.org/user_builds/debiased-spatial-whittle/envs/latest/lib/python3.10/site-packages/scipy/optimize/_lbfgsb_py.py:347, in _minimize_lbfgsb(fun, x0, args, jac, bounds, disp, maxcor, ftol, gtol, eps, maxfun, maxiter, iprint, callback, maxls, finite_diff_rel_step, **unknown_options)
    344         iprint = disp
    346 # _prepare_scalar_function can use bounds=None to represent no bounds
--> 347 sf = _prepare_scalar_function(fun, x0, jac=jac, args=args, epsilon=eps,
    348                               bounds=bounds,
    349                               finite_diff_rel_step=finite_diff_rel_step)
    351 func_and_grad = sf.fun_and_grad
    353 fortran_int = _lbfgsb.types.intvar.dtype

File ~/checkouts/readthedocs.org/user_builds/debiased-spatial-whittle/envs/latest/lib/python3.10/site-packages/scipy/optimize/_optimize.py:288, in _prepare_scalar_function(fun, x0, jac, args, bounds, epsilon, finite_diff_rel_step, hess)
    284     bounds = (-np.inf, np.inf)
    286 # ScalarFunction caches. Reuse of fun(x) during grad
    287 # calculation reduces overall function evaluations.
--> 288 sf = ScalarFunction(fun, x0, args, grad, hess,
    289                     finite_diff_rel_step, bounds, epsilon=epsilon)
    291 return sf

File ~/checkouts/readthedocs.org/user_builds/debiased-spatial-whittle/envs/latest/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:166, in ScalarFunction.__init__(self, fun, x0, args, grad, hess, finite_diff_rel_step, finite_diff_bounds, epsilon)
    163     self.f = fun_wrapped(self.x)
    165 self._update_fun_impl = update_fun
--> 166 self._update_fun()
    168 # Gradient evaluation
    169 if callable(grad):

File ~/checkouts/readthedocs.org/user_builds/debiased-spatial-whittle/envs/latest/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:262, in ScalarFunction._update_fun(self)
    260 def _update_fun(self):
    261     if not self.f_updated:
--> 262         self._update_fun_impl()
    263         self.f_updated = True

File ~/checkouts/readthedocs.org/user_builds/debiased-spatial-whittle/envs/latest/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:163, in ScalarFunction.__init__.<locals>.update_fun()
    162 def update_fun():
--> 163     self.f = fun_wrapped(self.x)

File ~/checkouts/readthedocs.org/user_builds/debiased-spatial-whittle/envs/latest/lib/python3.10/site-packages/scipy/optimize/_differentiable_functions.py:145, in ScalarFunction.__init__.<locals>.fun_wrapped(x)
    141 self.nfev += 1
    142 # Send a copy because the user may overwrite it.
    143 # Overwriting results in undefined behaviour because
    144 # fun(self.x) will change self.x, with the two no longer linked.
--> 145 fx = fun(np.copy(x), *args)
    146 # Make sure the function returns a true scalar
    147 if not np.isscalar(fx):

File ~/checkouts/readthedocs.org/user_builds/debiased-spatial-whittle/checkouts/latest/src/debiased_spatial_whittle/inference/likelihood.py:699, in Estimator._get_opt_func.<locals>.func(param_values)
    697 def func(param_values):
    698     model.update_free_parameters(param_values)
--> 699     return self.likelihood(z, model).item()

File ~/checkouts/readthedocs.org/user_builds/debiased-spatial-whittle/checkouts/latest/src/debiased_spatial_whittle/inference/likelihood.py:74, in MultivariateDebiasedWhittle.__call__(self, z, model, params_for_gradient)
     67 def __call__(
     68     self,
     69     z: xp.ndarray,
     70     model: CovarianceModel,
     71     params_for_gradient: list[ModelParameter] = None,
     72 ):
     73     """Computes the likelihood for these data"""
---> 74     p = self.periodogram([z[..., 0], z[..., 1]])
     75     ep = self.expected_periodogram(model)
     76     n_spatial_dim = p.ndim - 2

File ~/checkouts/readthedocs.org/user_builds/debiased-spatial-whittle/checkouts/latest/src/debiased_spatial_whittle/inference/periodogram.py:201, in Periodogram.__call__(self, sample)
    199         z_values = sample.values * self.taper(sample.grid.n)
    200 else:
--> 201     z_values = sample * self.taper(sample.shape)
    202 f = 1 / prod_list(z_values.shape) * xp.abs(fftn(z_values)) ** 2
    203 if isinstance(sample, SampleOnRectangularGrid):

AttributeError: 'list' object has no attribute 'shape'

Mapping¶

In [9]:

import numpy as np

import numpy as np

In [10]:

xs = np.stack(np.mgrid[:64, :64], -1).reshape(-1, 2)
ys = model.predict(
    (xs, xs), (s[..., 0].reshape((-1, 1)), s[..., 1].reshape((-1, 1))), xs
)
ys = ys.reshape((64, 64))
plt.figure()
plt.imshow(ys, cmap="seismic")
plt.show()

xs = np.stack(np.mgrid[:64, :64], -1).reshape(-1, 2) ys = model.predict( (xs, xs), (s[..., 0].reshape((-1, 1)), s[..., 1].reshape((-1, 1))), xs ) ys = ys.reshape((64, 64)) plt.figure() plt.imshow(ys, cmap="seismic") plt.show()