Skip to content
Snippets Groups Projects
Commit bcbc64a3 authored by MALOU THIBAULT's avatar MALOU THIBAULT
Browse files

update the documentation of the obs module to the Docstring NumpyDoc format

parent 4b399595
No related branches found
No related tags found
No related merge requests found
Pipeline #221518 passed
Hide whitespace changes
Inline Side-by-side
src/source_localization/obs.py
+ 175
− 117
View file @ bcbc64a3
import numpy as np
import scipy.sparse as sps
"""
Module that contains all the features related to the observations and the sensors
"""
class Obs:
"""
Class containing:
- the observations, its value, the time and space location of the observations,
- the estimation of the observed variable,
- the observations error covariance matrix,
- the observation operator and the adjoint operator of its derivative with respect of the state variable.
In this package, the observations are the accumulation of pheromone over a given time window.
Therefore, the observation operator is the integral of the pheromone concentration over this time window.
r"""
Class containing all the features related to the sensors and the observations.
This includes:
- the observations :math:`m^{obs}`,
the observation times :math:`t^{obs}` and
the observation positions :math:`X^{obs}=(x^{obs}, y^{obs})`,
- the estimation of the observed variable :math:`m(c(s))`,
- the concentration of pheromone :math:`c(s)` at the time and positions required to compute the observed variable,
- the observation operator :math:`c\mapsto m(c)`
- the adjoint operator of its derivative with respect of the state variable :math:`\phi\mapsto (d_cm(c))^*\cdot\phi`,
- the covariance matrix of the observation error :math:`\mathbf{R}`.
In this class, the observations are the accumulation of pheromone over a given time window.
Therefore, the observation operator is
the integral of the pheromone concentration over this time window:
.. math::
m(c)_{i}=\int_{t^{obs}_i-\delta t}^{t^{obs}_i}c(x^{obs}_i, y^{obs}_i, \tau)d\tau
Attributes
----------
msh : ~pheromone_dispersion.geom.MeshRect2D
The geometry of the domain.
t_obs : ~numpy.ndarray
The observation times :math:`t^{obs}`.
X_obs : ~numpy.ndarray
The locations of the observations :math:`X^{obs}=(x^{obs}, y^{obs})`.
d_obs : ~numpy.ndarray
The observations data :math:`m^{obs}`.
N_obs : int
The number of observations.
nb_sensors : int
The number of sensors.
X_sensors : ~numpy.ndarray
The positions of the sensors.
dt_obs_operator : float
Length of the time window of accumulation :math:`\delta t`.
If the observation is instantaneous, set to `0`.
nb_dt_in_obs_time_window : int
Number of time steps covering the time window of accumulation :math:`[t^{obs}_i-\delta t; t^{obs}_i]`.
index_obs_to_index_time_est : ~numpy.ndarray
The indexes of the times
at which the state variable is required to compute the observation variable given the index of an observation.
index_time_est_to_index_obs : dict
Dictionary where each key is a time index, and the value is an array of indexes of the observations that
require the estimate of the state variable at this time to estimate the observed variable.
c_est : ~numpy.ndarray
The estimation the concentration of pheromone :math:`c(s;x,y,t)`
at the times :math:`t` and positions :math:`(x,y)`
required to compute the observed variable :math:`m(c(s))`.
d_est : ~numpy.ndarray
The estimation of the accumulation of pheromone (observed variable)
computed using the pheromone propagation model :math:`m(c(s))`.
R_inv : :mod:`~scipy.sparse` matrix, default: :func:`~scipy.sparse.identity` matrix
The inverse of the covariance matrix of the observation error :math:`\mathbf{R}^{-1}`.
"""
def __init__(self, t_obs, X_obs, d_obs, msh, index_sensor=None, dt_obs_operator=0.0):
"""
Instanciation of the class.
- TO DO:
* update the documentation and comments for the new observation operator
* add an exception to make sure that the accumulation time window does not overlapp for a given sensor
* add an exception to make sure that no observations are made after
- input:
* t_obs:
| array of the shape (N_obs,),
| contains the observations times
* X_obs:
| array of the shape (N_obs,2),
| contains the locations of the observations
* d_obs:
| array of the shape (N_obs,),
| contains the values of the observations
* msh:
| object of the class MeshRect2D
- attributes:
* msh:
| object of the class MeshRect2D
* t_obs:
| array of the shape (N_obs,),
| contains the observations times
* X_obs:
| array of the shape (N_obs,2),
| contains the locations of the observations
* d_obs:
| array of the shape (N_obs,),
| contains the values of the observations
* N_obs:
| integer,
| contains the number of observations
* nb_sensors:
| integer,
| contains the number of sensors
* X_sensors:
| array of shape (nb_sensors, 2),
| contains the position of the sensors
* dt_obs_operator:
| float,
| contains the length of the time window of accumulation
* nb_dt_in_obs_time_window:
| integer,
| contains the number of time steps covering the time window of accumulation
* index_obs_to_index_time_est:
| array of shape (N_obs, nb_dt_in_obs_time_window),
| contains the indexes of the times at which the state variable is required to compute the observation variable
| given the index of an observation
* index_time_est_to_index_obs:
| dictionnary of array of integers,
| for a time as key, contains the array of index of the observations
| that require the estimate of the state variable at this time to estimate the observed variable
* c_est:
| array of the shape (N_obs, nb_dt_in_obs_time_window),
| contains the estimation of the state variable (concentration) computed by the direct model
| at the time and space location required to estimate the observed variable
* d_est:
| array of the shape (N_obs,),
| contains the estimation of the observed variable (accumulation of pheromone) computed by the direct model
| at the time and space location of the observations
* R_inv:
| sparse array of shape (N_obs, N_obs), the identity matrix by default
| contains the inverse of the observations error covariance matrix
def __init__(self, t_obs, X_obs, d_obs, msh, dt_obs_operator=0.0):
r"""
Constructor method.
Parameters
----------
t_obs : ~numpy.ndarray
The observation times :math:`t^{obs}`.
X_obs : ~numpy.ndarray
The locations of the observations :math:`X^{obs}=(x^{obs}, y^{obs})`.
d_obs : ~numpy.ndarray
The observations data :math:`m^{obs}`.
msh : ~pheromone_dispersion.geom.MeshRect2D
The geometry of the domain.
dt_obs_operator : float, default: 0
Length of the time window of accumulation :math:`\delta t`.
If set to `0`, the observation is instantaneous.
Raises
------
ValueError
if the time array has not been initialized using
the methods :func:`~pheromone_dispersion.geom.MeshRect2D.calc_dt_explicit_solver`
or :func:`~pheromone_dispersion.geom.MeshRect2D.calc_dt_implicit_solver`.
ValueError
if the sizes of
the arrays containing the observation times :math:`t^{obs}`,
the observation positions :math:`X^{obs}` and
the observations data :math:`m^{obs}`
do not coincide.
ValueError
if one observation time :math:`t^{obs}`
or one observation position :math:`X^{obs}`
is out of the time window :math:`[0;T]` or domain :math:`\Omega`.
"""
# To do:
# update the documentation and comments for the new observation operator
# add an exception to make sure that the accumulation time window does not overlapp for a given sensor
# add an exception to make sure that no observations are made after
self.msh = msh
self.t_obs = np.copy(t_obs)
self.X_obs = np.copy(X_obs)
......@@ -170,14 +195,31 @@ class Obs:
self.index_time_est_to_index_obs[time_idx] = [increment]
def obs_operator(self):
"""
the observation operator
i.e. integration of the state variable of the direct model over the time window [t_obs-dt_obs_operator; t_obs]
assuming that the integration time window for different data of a same sensor are not overlapping
r"""
The observation operator :math:`c\mapsto m(c)`.
In this class, the observation operator is,
for the :math:`i^{th}` observation time and position,
:math:`m(c)_{i}=\int_{t^{obs}_i-\delta t}^{t^{obs}_i}c(x^{obs}_i, y^{obs}_i, \tau)d\tau`.
Moreover, it is assumed that the integration time window for different data of a same sensor are not overlapping
- do:
* compute an estimation of the observed variable given an estimation of the state variable of the direct model
* store the estimation of the observed variable in the attribute d_est
Notes
-----
This method computes the observed variable and
store the values in the attribute :attr:`d_est`
using the attribute :attr:`c_est`
that contains the estimation the concentration of pheromone :math:`c(s;x,y,t)`
computed by the pheromone propagation model
at the times :math:`t` and positions :math:`(x,y)` required to compute the observed variable.
To get these :math:`c(s;x,y,t)`,
the attribute :attr:`c_est` should be computed in a first time
using the method :meth:`~pheromone_dispersion.convection_diffusion_2D.DiffusionConvectionReaction2DEquation.solver_est_at_obs_times`
of the class :meth:`~pheromone_dispersion.convection_diffusion_2D.DiffusionConvectionReaction2DEquation`.
Raises
------
ValueError
if the attribute :attr:`c_est` has not been computed in a first time.
"""
if np.array([np.isnan(self.c_est[i]) for i in range(self.N_obs)]).any():
......@@ -187,24 +229,33 @@ class Obs:
)
self.d_est = np.sum(self.c_est, axis=1) * self.msh.dt
def onesensor_adjoint_derivative_obs_operator(self, t, delta_c, index_sensor):
"""
compute the spatial map of the adjoint of the derivative of the one-sensor observation operator with respect to the state variable
at a given time
- input:
* t:
| float,
| the current time
* delta_c:
| numpy array of shape (N_obs, ),
| element of the space of the concentration (state variable)
* index_sensor:
| integer,
| index of the sensor to consider
- output:
| numpy array of shape (msh.x.size*msh.y.size, ),
| evaluation of the adjoint of the derivative of the one-sensor observation operator in delta_c for the index_sensor-th sensor
def onesensor_adjoint_derivative_obs_operator(self, t, phi, index_sensor):
r"""
compute the spatial map
of the adjoint of the derivative of the one-sensor observation operator
with respect to the state variable
:math:`\phi\mapsto (d_cm_i(c))^*\cdot\phi`,
at a given time :math:`t`,
for the :math:`i^{th}` sensor and
for an element of the observation space :math:`\phi`.
Parameters
----------
t : float or integer
The current time :math:`t`.
phi : ~numpy.ndarray
The element of the observation space :math:`\phi` at current time :math:`t`.
index_sensor: int
The index of the sensor :math:`i`.
Returns
-------
~numpy.ndarray
The spatial map of the image
of the adjoint of the derivative of the one-sensor observation operator
:math:`(d_cm_i(c))^*\cdot\phi(x,y)`
at the given time
raveled into a vector.
"""
index_t = np.argmin(np.abs(t - self.msh.t_array))
out = np.zeros((self.msh.y.size, self.msh.x.size))
......@@ -217,24 +268,31 @@ class Obs:
index_obs = index_obs_current_t[index]
index_x_est = np.argmin(np.abs(self.msh.x - X_sensor[0]))
index_y_est = np.argmin(np.abs(self.msh.y - X_sensor[1]))
out[index_y_est, index_x_est] = delta_c[index_obs] * self.msh.dt
out[index_y_est, index_x_est] = phi[index_obs] * self.msh.dt
return out.reshape((self.msh.y.size * self.msh.x.size,))
def adjoint_derivative_obs_operator(self, t, delta_c):
"""
compute the spatial map of the adjoint of the derivative of the observation operator with respect to the state variable
at a given time
- input:
* t:
| float,
| the current time
* delta_c:
| numpy array of shape (N_obs, ),
| element of the space of the concentration (state variable)
- output:
| numpy array of shape (msh.x.size*msh.y.size, ),
| evaluation of the adjoint of the derivative of the observation operator in delta_c
def adjoint_derivative_obs_operator(self, t, phi):
r"""
compute the spatial map of
the adjoint of the derivative of the observation operator with respect to the state variable
:math:`\phi\mapsto (d_cm(c))^*\cdot\phi`,
at a given time :math:`t`
and for an element of the observation space :math:`\phi`.
Parameters
----------
t : float or integer
The current time :math:`t`.
phi : ~numpy.ndarray
The element of the observation space :math:`\phi` at current time :math:`t`.
Returns
-------
~numpy.ndarray
The spatial map of the image of the adjoint of the derivative of the observation operator
:math:`(d_cm(c))^*\cdot\phi(x,y)`
at the given time
raveled into a vector.
"""
index_t = np.argmin(np.abs(t - self.msh.t_array))
out = np.zeros((self.msh.y.size, self.msh.x.size))
......@@ -253,5 +311,5 @@ class Obs:
# for each observation points, store the associated input at the correponding estimate point
for i_x, i_y, i_obs in zip(index_x_est, index_y_est, index_obs):
out[i_y, i_x] = delta_c[i_obs] * self.msh.dt
out[i_y, i_x] = phi[i_obs] * self.msh.dt
return out.reshape((self.msh.y.size * self.msh.x.size,))
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment