import numpy as np
import scipy.sparse as sps
[docs]
class Control:
r"""
Class containing the control variable :math:`s(x,y,t)`
that is the quantity of pheromone emitted by the insects
and all its prior known features.
Attributes
----------
value : ~numpy.ndarray
The current estimation of the control variable :math:`s(x,y,t)`
raveled in a (:attr:`msh.t_array.size` * :attr:`msh.y.size` * :attr:`msh.x.size`,)-shape array. Initialized by zeros.
background_value : ~numpy.ndarray
The background/prior estimation of the control variable :math:`s_b(x,y,t)`
raveled in a (:attr:`msh.t_array.size` * :attr:`msh.y.size` * :attr:`msh.x.size`,)-shape array.
C_inv : ~numpy.ndarray
The inverse of the covariance matrix of the background error
(error between the prior estatimation and the true value) :math:`\mathbf{C}^{-1}`.
Initialized and by default the identity.
population_dynamic_model : subclass of ~scipy.sparse.linalg.LinearOperator or ~numpy.ndarray
The population dynamique model :math:`\partial_ts+\mathcal{M}(s)=0` that
the control variable satisfies.
Examples of population dynamique models can be found
in the submodule :mod:`~source_localization.population_dynamique`.
The user can easily provide its own population dynamic model as long as it contains
a method :meth:`_matvec` that estimate the residual of the model :math:`s\mapsto\partial_ts+\mathcal{M}(s)` and
a method :meth:`_rmatvec` that estimate the adjoint of the derivative of the model with respect to the control variable
:math:`\delta s\mapsto(\partial_t+d_s\mathcal{M}(s))^*\cdot\delta s`.
If no population dynamic models are a priori known, set to `None`.
exclusion_domain : ~numpy.ndarray of boolean or list
The exclusion domain :math:`\Omega_{exc}` corresponding to the time and places
where it is a prior known that no insects will emit pheromones
for the log-barrier regularization term.
Is `True` where the pheromone emissions should be excluded.
If no exclusion zones are a priori known, set to an empty list.
log_barrier_threshold : float
The threshold of the log-barrier regularization :math:`\epsilon`.
"""
[docs]
def __init__(self, background_source, msh, exclusion_domain=None, log_barrier_threshold=1e-10, population_dynamique_model=None):
r"""
Constructor method
Parameters
----------
background_source : ~pheromone_dispersion.source_term.Source
The background value of the source term :math:`s_b(x,y,t)`.
msh : ~pheromone_dispersion.geom.MeshRect2D
The geometry of the domain.
population_dynamic_model : subclass of ~scipy.sparse.linalg.LinearOperator or ~numpy.ndarray or None, optional, default: None
The population dynamique model :math:`\partial_ts+\mathcal{M}(s)=0` that
the control variable satisfies. `None` if no population dynamic models are known.
exclusion_domain : ~numpy.ndarray of boolean or list, optional, default: None
The exclusion domain :math:`\Omega_{exc}`. `None` if no exclusion domain are known.
log_barrier_threshold : float, optional, default: 1e-10
The threshold of the log-barrier regularization :math:`\epsilon`.
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.MeshRect2Dcalc_dt_implicit_solver`.
"""
# To do:
# add tests to check the provided population_dynamique_model
# add tests and exceptions for the log-barrier reg
if msh.t_array is None:
raise ValueError(
"The array containing the time at which the direct model is solved, has not been correctly initialized."
+ " This is probably because the time step has not been computed so the CFL condition is satisfied."
+ " The method calc_CFL of the class MeshRect2D should then be executed first."
)
background_value = np.array([])
for it, t in enumerate(msh.t_array):
background_source.at_current_time(t)
value_prov = np.copy(background_source.value) # NECESSARY???
background_value = np.append(background_value, value_prov.ravel())
self.background_value = np.copy(background_value)
self.value = np.zeros(background_value.shape)
self.C_inv = sps.identity(np.size(background_value), dtype=int)
self.population_dynamic_model = population_dynamique_model
if exclusion_domain is None:
self.exclusion_domain = []
else:
self.exclusion_domain = exclusion_domain
self.log_barrier_threshold = log_barrier_threshold
[docs]
def apply_control(self, direct_model):
"""
apply the current value of the control variable as source term :math:`s(x,y,t)` of the given pheromone propagation model
Parameters
----------
direct_model: ~pheromone_dispersion.convection_diffusion_2D.DiffusionConvectionReaction2DEquation
The pheromone propagation model
Notes
-----
Update the attribute :attr:`~pheromone_dispersion.convection_diffusion_2D.DiffusionConvectionReaction2DEquation.S`
of the given object of the class :class:`~pheromone_dispersion.convection_diffusion_2D.DiffusionConvectionReaction2DEquation`
with the current value of the control variable contained in the attribute :attr:`value`.
"""
value_prov = np.copy(self.value) # NECESSARY???
direct_model.set_source(
value_prov.reshape((direct_model.msh.t_array.shape[0], direct_model.msh.y.shape[0], direct_model.msh.x.shape[0])),
t=direct_model.msh.t_array,
)