import numpy as np
from scipy.sparse.linalg import LinearOperator as LinOp
[docs]
class StationnaryPopulationDynamicModel(LinOp):
r"""
Class containing the operator of a stationary population dynamic model :math:`\partial p(x,y,t)=0`
with :math:`p` the density of insect
and a non-zeros constant emission rate :math:`q`.
As the quantity of pheromone emitted by the insects is :math:`s=pq`, this leads to :math:`\partial_ts(x,y,t)=0`.
The implementation of this linear operator is a subclass of the :py:class:`~scipy.sparse.linalg.LinearOperator` class.
Attributes
-----------
msh : ~pheromone_dispersion.geom.MeshRect2D
The geometry of the domain.
shape : tuple of int
Shape of the matrix of the linear operator.
The shape is (:attr:`msh.t_array.size-1` * :attr:`msh.y.size` * :attr:`msh.x.size`,
:attr:`msh.t_array.size` * :attr:`msh.y.size` * :attr:`msh.x.size`).
dtype : ~numpy.dtype
Data type of the elements of the matrix of the linear operator.
The type is `float64`
"""
[docs]
def __init__(self, msh):
r"""
Constructor method
Parameters
----------
msh : ~pheromone_dispersion.geom.MeshRect2D
The geometry of the domain.
"""
self.msh = msh
self.shape = (
self.msh.y.size * self.msh.x.size * (self.msh.t_array.size - 1),
self.msh.y.size * self.msh.x.size * self.msh.t_array.size,
)
self.dtype = np.dtype(float)
[docs]
def _matvec(self, x_out):
r"""
Compute the image (matrix-vector product) of the stationnary population dynamic model operator
:math:`\mathcal{M}:s\mapsto\partial_ts`
for a given quantity of pheromone emitted :math:`s(x,y,t)`.
Parameters
----------
x_out : ~numpy.ndarray
The quantity of pheromone emitted :math:`s(x,y,t)` raveled into a vector.
Returns
-------
~numpy.ndarray
Array containing the image :math:`s\mapsto\mathcal{M}(s)`.
"""
S = x_out.reshape((self.msh.t_array.shape[0], self.msh.y.shape[0], self.msh.x.shape[0]))
# time derivative of the S
d_tS = np.zeros((self.msh.t_array.shape[0] - 1, self.msh.y.shape[0], self.msh.x.shape[0]))
d_tS = S[1:, :, :] - S[:-1, :, :]
return d_tS.reshape((self.msh.y.size * self.msh.x.size * (self.msh.t_array.size - 1),)) / self.msh.dt
[docs]
def _rmatvec(self, x_out):
r"""
Compute the image (matrix-vector product) of the stationnary population dynamic model adjoint operator
:math:`\mathcal{M}^*:\phi\mapsto-\partial_t\phi`
with :math:`\phi(t=0)=\phi(t=T)=0`
for a given element :math:`\phi(x,y,t)`.
Parameters
----------
x_out : ~numpy.ndarray
The element :math:`\phi(x,y,t)` raveled into a vector.
Returns
-------
~numpy.ndarray
Array containing the image :math:`\phi\mapsto\mathcal{M}^*(\phi)`.
"""
S = x_out.reshape((self.msh.t_array.shape[0] - 1, self.msh.y.shape[0], self.msh.x.shape[0]))
# time derivative of the S
d_tS_adjoint = np.zeros((self.msh.t_array.shape[0], self.msh.y.shape[0], self.msh.x.shape[0]))
d_tS_adjoint[0, :, :] = -S[0, :, :]
d_tS_adjoint[1:-1, :, :] = -S[1:, :, :] + S[:-1, :, :]
d_tS_adjoint[-1, :, :] = S[-1, :, :]
return d_tS_adjoint.reshape((self.msh.y.size * self.msh.x.size * self.msh.t_array.size,)) / self.msh.dt
[docs]
class PopulationDynamicModel(LinOp):
r"""
Class containing the operator of a population dynamic model :math:`\partial p(x,y,t)=-\gamma(x,y,t) p(x,y,t)`
with :math:`p` the density of insect,
:math:`\gamma` the constant death rate,
and a non-zeros constant emission rate :math:`q`.
As the quantity of pheromone emitted by the insects is :math:`s=pq`, this leads to :math:`\partial s(x,y,t)=-\gamma s(x,y,t)`.
The implementation of this linear operator is a subclass of the :py:class:`~scipy.sparse.linalg.LinearOperator` class.
Attributes
-----------
msh : ~pheromone_dispersion.geom.MeshRect2D
The geometry of the domain.
death_rate: ~numpy.dtype
The death rate of the population dynamic model :math:`\gamma(x,y,t)`.
shape : tuple of int
Shape of the matrix of the linear operator.
The shape is (:attr:`msh.t_array.size-1` * :attr:`msh.y.size` * :attr:`msh.x.size`,
:attr:`msh.t_array.size` * :attr:`msh.y.size` * :attr:`msh.x.size`).
dtype : ~numpy.dtype
Data type of the elements of the matrix of the linear operator.
The type is `float64`.
"""
[docs]
def __init__(self, msh, death_rate=0.1):
r"""
Constructor method
Parameters
----------
msh : ~pheromone_dispersion.geom.MeshRect2D
The geometry of the domain.
death_rate: ~numpy.dtype, optional, float or None, default: 0.1
The death rate of the population dynamic model :math:`\gamma(x,y,t)`.
If `None`, set to `0`.
If `float`, set to constante.
"""
self.msh = msh
if isinstance(death_rate, float) or isinstance(death_rate, int):
self.death_rate = death_rate * np.ones((self.msh.t_array.shape[0], self.msh.y.shape[0], self.msh.x.shape[0]))
elif isinstance(death_rate, np.ndarray):
self.death_rate = death_rate
elif death_rate is None:
self.death_rate = np.zeros((self.msh.t_array.shape[0], self.msh.y.shape[0], self.msh.x.shape[0]))
self.shape = (
self.msh.y.size * self.msh.x.size * (self.msh.t_array.size - 1),
self.msh.y.size * self.msh.x.size * self.msh.t_array.size,
)
self.dtype = np.dtype(float)
[docs]
def _matvec(self, x_out):
r"""
Compute the image (matrix-vector product) of the population dynamic model operator :math:`\mathcal{M}:s\mapsto\partial_ts+\gamma s`
for a given quantity of pheromone emitted :math:`s(x,y,t)`.
Parameters
----------
x_out : ~numpy.ndarray
The quantity of pheromone emitted :math:`s(x,y,t)` raveled into a vector.
Returns
-------
~numpy.ndarray
Array containing the image :math:`s\mapsto\mathcal{M}(s)`.
"""
S = x_out.reshape((self.msh.t_array.shape[0], self.msh.y.shape[0], self.msh.x.shape[0]))
res = (S[1:, :, :] - S[:-1, :, :]) / self.msh.dt
res += np.multiply(self.death_rate[1:, :, :], S[1:, :, :])
return res.reshape((self.msh.y.size * self.msh.x.size * (self.msh.t_array.size - 1),))
[docs]
def _rmatvec(self, x_out):
r"""
Compute the image (matrix-vector product) of the population dynamic model adjoint operator
:math:`\mathcal{M}^*:\phi\mapsto-\partial_t\phi+\gamma\phi`
with :math:`\phi(t=0)=\phi(t=T)=0`
for a given element :math:`\phi(x,y,t)`.
Parameters
----------
x_out : ~numpy.ndarray
The element :math:`\phi(x,y,t)` raveled into a vector.
Returns
-------
~numpy.ndarray
Array containing the image :math:`\phi\mapsto\mathcal{M}^*(\phi)`.
"""
S = x_out.reshape((self.msh.t_array.shape[0] - 1, self.msh.y.shape[0], self.msh.x.shape[0]))
res = np.zeros((self.msh.t_array.shape[0], self.msh.y.shape[0], self.msh.x.shape[0]))
res[0, :, :] = -S[0, :, :] / self.msh.dt
res[1:-1, :, :] = (-S[1:, :, :] + S[:-1, :, :]) / self.msh.dt
res[-1, :, :] = S[-1, :, :] / self.msh.dt
res[1:, :, :] += np.multiply(self.death_rate[1:, :, :], S)
return res.reshape((self.msh.y.size * self.msh.x.size * self.msh.t_array.size,))