Inverse problem: Biology-informed assimilation of pheromone sensors data in the pheromone propagation model

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
from pathlib import Path
import os
import time

import pherosensor

from pheromone_dispersion.geom import MeshRect2D
from pheromone_dispersion.diffusion_tensor import DiffusionTensor
from pheromone_dispersion.velocity import Velocity
from pheromone_dispersion.source_term import Source
from pheromone_dispersion.convection_diffusion_2D import DiffusionConvectionReaction2DEquation

from source_localization.obs import Obs
from source_localization.control import Control
from source_localization.cost import Cost
from source_localization.population_dynamique import PopulationDynamicModel
from source_localization.adjoint_convection_diffusion_2D import AdjointDiffusionConvectionReaction2DEquation

from utils.plot_colormap import plot_colormap
from utils.plot_obs import plot_obs
from utils.plot_ctrl import plot_ctrl
from utils.plot_cost import plot_cost

path = os.getcwd()
path_data = path + '/data'
path_output = path + '/output'
if not os.path.isdir(Path(path_output)):
    os.makedirs(Path(path_output))
path_plot = path_output + '/plot'
if not os.path.isdir(Path(path_plot)):
    os.makedirs(Path(path_plot))

This notebook aims to present the biology-informed data assimilation (BI-DA) problem solved to infer the quantity of pheromone emitted from pheromone sensors data, the solvers implemented in the module source_localization and the different submodules on which the solvers are based.

In a first section, the data assimilation problem is presented. Then, the different submodules are presented and illustrated over a simple but representative case, in line with the case used in the notebook dedicated to modelling pheromone propagation. Finally, the numerical solvers of the BI-DA problem are presented and used to infer the quantity of pheromone emitted in different scenario of biological information addition.

Pheromone sensors data assimilation for the inference of the quantity of emitted pheromone

The goal of data assimilation (DA) is to estimate the quantity of pheromone emitted by the insects, also referred to as the control variable, using time-series of pheromone sensors spatially distributed over the landscape of interest.

To do so, we solve an inverse problem corresponding to the following direct problem describing pheromone propagation over the landscape (see equation (1) in the notebook modelling pheromone propagation for details)

\[\frac{\partial c}{\partial t}-\nabla\cdot(\mathbf{K}\nabla c)+\nabla\cdot(\vec{u}c)+\tau_{loss}c=s~\forall (x,y,t)\in\Omega\times[0;T] \quad \quad \quad \text{(1)}\]

In this equation, the unknown is the pheromone density \(c\), and \(\mathbf{K}\), \(\vec{u}\) and \(\tau_{loss}\) are known parameters describing respectively the diffusion tensor, the advection by the wind and a deposition parameter. The term \(s\) is the emission rate of pheromone, which is the target to be inferred by the inverse problem.

The method consists in finding the quantity of pheromone emitted \(s\) that enables the prediction of the pheromone propagation model (1) that predicts the closest outputs to the data.

Therefore, the DA problem is defined as a minimization problem over a cost function \(j\):

\[\begin{split}\begin{cases} \text{find }s_a(x,y,t)\text{ such that} \\ s_a(x,y,t)=\text{arg}\min_{s(x,y,t)}j(s)\text{ with } j(s)=j_{obs}(s)+j_{reg}(s) \end{cases}\end{split}\]

The term \(j_{obs}\) aims to minimize the gap between model predictions and data. It reads

\[j_{obs}(s)=\|m\left(c(s)\right)-m_{obs}\|_{\mathbf{R}^{-1}}^2\]

where \(m_{obs}\) are the observations, \(c(s)\) is the pheromone concentration map predicted by the pheromone propagation model (1) from the given pheromone emission \(s\), \(m(c)\) is the observation operator that models the observed variable from a given pheromone concentration map and \(\mathbf{R}\) is the covariance operator of the observation error. The observation operator depends on the properties of the sensor. For instance, pheromone sensors commonly use a pre-concentrator that actively filters air during a time window \(\delta_t\) and releases the accumulated quantity into the sensor for the measure. In this case, the observation operator is \((m(c))_{i,k}= \int_{t_{obs,k}-\delta_t}^{t_{obs,k}} c(\tau,x_{obs, i},y_{obs, i}) d\tau\) where \((x_{obs, i},y_{obs, i})\) is the position of the \(i^{th}\) sensor and \(t_{obs,k}\) the \(k^{th}\) time of observation.

The regularization term \(j_{reg}\) is commonly added in order to ensure the well-posedness of the optimization problem. In the BI-DA framework, this regularization term is based on prior biological information describing biological knowledge of the emitter insects. A short introduction of this biology-informed regularization term can be found below and a more extensive one can be found in Malou et al (2024).

Biology-informed regularization terms

The BI-DA framework is based on regularization terms modeling biological information that enable to drive the DA method toward an optimum \(s_a\) consistent with the biological knowledge introduced in \(j_{reg}\). In general terms, the regularization term reads:

\[j_{reg} = \sum_{i \in \mathcal{L}} \alpha_{reg,i} j_{reg,i} \quad \quad \quad \quad \text{ (2)}\]

with \(\mathcal{L}\) the list of considered BI regularization terms \(j_{reg,i}\) and \(\alpha_{reg,i}\) the associated weight coefficient. This allows for composite regularization terms that take into account several biological informations.

Several biology-informed regularization terms are shortly presented here.

• The Tikhonov regularization \(j_{reg,T}(s)=\|s-s_b\|^2_{\mathbf{C}^{-1}}\) enables to take into account information on the insects life habits, such as their preferred habitat and insect life cycle, through \(s_b\) a prior value of \(s\) and \(\mathbf{C}\) the background error covariance operator.

• The LASSO regularization \(j_{reg,LASSO}(s)=\|s\|_{1}\) enables to select important local pheromones emissions at the beginning of the infestation, when only few insects are expected. The LASSO term enables to discard sporious emissions.

• The population dynamics informed regularization \(j_{reg,PD}(s)=\|\partial_ts-\frac{\partial_tq}{q}s+q\nabla \left(\sum_i \mathcal{F}_i(s)\right)-q\sum_i\mathcal{R}_i(s)\|^2_{L^2}\) (with \(q\) the quantity of pheromone emitted per insect) enables to take into account the spatio-temporal evolution of the density of insects \(p=s/q\) described by the flux \(F_i\) and reaction \(R_i\) terms that model the behavior of the insect of interest.

Solving the optimization problem

To solve the minimization problem, gradient-based methods are used. Thus, at each descent step, the gradient \(\nabla_sj\) must be computed and reads:

\[\nabla_s j(s) = \nabla_s j_{obs}(s) + \nabla_s j_{reg}(s)\]

Note that we denote \(\nabla_s\) the gradient with respect to the control \(s\) to avoid confusion with \(\nabla\) the classical gradient with respect to the space variables. Since \(s\) is defined in a high-dimensional space and \(c(s)\) involves the computation of a PDE model, the numerical computation of its gradient is challenging. We use a variational method and compute the gradient by using the adjoint model of the CTM that reads:

\[\partial_tc^* + \nabla\cdot(\mathbf{K}^T\nabla c^*)+\nabla(\vec{u}c^*)-(\nabla.\vec{u})c^*-\tau_{loss}c^* = \left(\frac{dm}{dc}(c(s))\right)^*\cdot2\mathbf{R}^{-1}\left(m(c(s))-m_{obs}\right) \quad \quad \quad \quad \text{ (3)}\]

with \(c^*\) the adjoint state.

The adjoint model is closed by the following final and boundary conditions:

• a null final condition \(c^*(x,y,t=T)=0\) \(\forall (x,y)\in\Omega\),

• a null diffusive flux \(\mathbf{K}^T\nabla c^*\cdot\vec{n}=0\) \(\forall (x,y)\in\partial\Omega\),

• a null outgoing convective flux \(\vec{u}c^*\cdot\vec{n}=0\) \(\forall (x,y)\in\partial\Omega\cap \{(x,y)|\vec{u}(x,y,t)\cdot\vec{n}>0\}~\forall t\in]0;T]\) with \(\vec{n}\) the outgoing normal vector.

Once the adjoint state is computed, the gradient \(\nabla_s j\) is computed using the expression:

\[\nabla_s j(s) = \nabla_s j_{reg}(s)-c^*(s)\]

The adjoint model is solved using an implicit downwind finite volume scheme on the same Cartesian grid as the direct problem. More details can be found in the notebooks located in the folder ./test_numerical_scheme/ that introduce and test the numerical scheme, including convergence, and the solvers.

Construction of the solver of the data assimilation problem

Specifying the geometry using the geom submodule

As mentionned in the notebook on modelling pheromone propagation, for now, the numerical solver of the pheromone propagation model has been implemented for a rectangular domain \(\Omega=[0;L_x]\times[0;L_y]\). Moreover, the mesh is supposed to be cartesian (rectangular cells) with constant space steps \(\Delta x\) and \(\Delta y\).

Such geometry can be generated as an object of the MeshRect2Dclass of the geom submodule.

Lx = 30 # the length along the x-axis
Ly = 40 # the length along the y-axis
Delta_x = 0.5 # the space step along the x-axis
Delta_y = 0.5 # the space step along the y-axis
T_final = 20 # the final time of the simulation
msh = MeshRect2D(Lx, Ly, Delta_x, Delta_y, T_final)
msh.calc_dt_implicit_solver(.1)

Specifying the pheromone propagation model and its solver using the convection_diffusion_2D submodule

The first step is specifying the pheromone propagation model using the convection_diffusion_2D submodule of the pheromone_propagation module. This requires to specify the wind velocity field, the diffusion tensor, the loss coefficient and to initialize a source term (potentially by filling it with 0 or a prior estimate).

The reader should refer to the dedicated tutorial notebook for detailed information on that topic.

U_vi = np.load(Path(path_data) / f'U_vi.npy')
U_hi = np.load(Path(path_data) / f'U_hi.npy')
U = Velocity(msh,U_vi, U_hi)

K = DiffusionTensor(U, 10., 10.)

deposition_coeff = 0.01*np.ones((msh.y.size, msh.x.size))

amp = 0.5
xc = 6.
yc = 6.
r = 1.5
S_value = np.zeros((msh.t_array.size, msh.y.size, msh.x.size))
xx, yy= np.meshgrid(msh.x, msh.y)
S_value[:, (xx-xc)**2 + (yy-yc)**2 < r**2] = amp

S_target = Source(msh, S_value, t=msh.t_array)

solver = 'implicit with stationnary matrix inversion'
tol = 1e-10
direct_model = DiffusionConvectionReaction2DEquation(
    U,
    K,
    deposition_coeff,
    S_target,
    msh,
    tol_inversion=tol,
    time_discretization=solver
)
fname = 'inv_matrix_implicit_solver'
direct_model.init_inverse_matrix(path_to_matrix=path_data, matrix_file_name=fname)

=== Load of the inverse of the matrix of the implicit part of the implicit with stationnary matrix inversion scheme ===

Specifying the pheromone sensors data and properties using the obs submodule

The sensors data and the observation operator are contained in an object of the Obs class of the obs submodule.

The user should provide the data array, that is of shape (nb data, 1), and arrays, resp. of shape (nb data, 1) and (nb data, 2), that contains the time and location of each data. The user can also specify the length of the time window over which the pheromone accumulates in the pre-concentration chamber of the sensor before measuring. If this time window length is not specified, it is assumed that the observation is instantaneous .

The Obs class also contains useful attributes and methods for pheromone sensors data and properties. For instance, the attribute c_est contains the estimate of the concentration in pheromone at the positions and times required to compute the estimate of the observed variable that is containes in the attribute d_est. Initially, these attributes are filled with NaN. The attribute c_est can be updated using the solver_est_at_obs_times method of the DiffusionConvectionReaction2DEquation class and the attribute d_est can be updated from the attribute c_est using the method obs_operator.

In the test case used here, the data are generated as following at random sensors location using the pheromone propagation model and the target source term introduced above.

if not (Path(path_data) / 'd_obs_noised.npy').exists():
    from random import choices, seed
    print("--- generating synthetic data from the target source term and random sensors location---")
    dt_sensor = 2.
    nb_sensor = 30
    seed(2)
    random_index_sensor_loc = choices(range(0, msh.x.size*msh.y.size), k=nb_sensor)
    X_sensor = np.zeros((nb_sensor, 2))
    index_sensor = 0
    for i, (xi,yi) in enumerate(zip(np.repeat(msh.x,msh.y.size),np.tile(msh.y,msh.x.size))):
        if i in random_index_sensor_loc:
            X_sensor[index_sensor,:]=[xi,yi]
            index_sensor += 1
    nb_obs_per_sensor = int(msh.T_final // dt_sensor)
    t_obs_per_sensor = np.array([dt_sensor*(i+1) for i in range(nb_obs_per_sensor)])
    X_obs = np.array([[x,y] for x,y in zip(np.tile(X_sensor[:,0],nb_obs_per_sensor),np.tile(X_sensor[:,1],nb_obs_per_sensor))])
    t_obs = np.repeat(t_obs_per_sensor, nb_sensor)
    d_obs = np.zeros_like(t_obs)
    obs_prov = Obs(t_obs, X_obs, d_obs, msh, dt_obs_operator=dt_sensor)
    direct_model.solver_est_at_obs_times(obs_prov)
    obs_prov.obs_operator()
    noise = np.random.normal(0., 0.01*np.max(obs_prov.d_est), size=obs_prov.d_est.size)
    noise_level = msh.mass_cell * msh.dt * (noise).dot(noise)
    np.save(Path(path_data) / 't_obs.npy', obs_prov.t_obs)
    np.save(Path(path_data) / 'X_obs.npy', obs_prov.X_obs)
    np.save(Path(path_data) / 'd_obs_no_noise.npy', obs_prov.d_est)
    np.save(Path(path_data) / 'd_obs_noised.npy', obs_prov.d_est+noise)
    np.save(Path(path_data) / 'noise_level.npy', noise_level)

--- generating synthetic data from the target source term and random sensors location---
t = 20.000 / 20.000 s

dt_sensor = 2.
t_obs = np.load(Path(path_data) / 't_obs.npy')
X_obs = np.load(Path(path_data) / 'X_obs.npy')
d_obs = np.load(Path(path_data) / 'd_obs_noised.npy')
obs = Obs(t_obs, X_obs, d_obs, msh, dt_obs_operator=dt_sensor)

noise_level = np.load(Path(path_data) / 'noise_level.npy')
print("||epsilon||_2^2 = ", noise_level)

||epsilon||_2^2 =  1.67083409330364e-05

Moreover, the data can be plotted at a given time using the plot_obs method of the utils submodule.

plot_obs(msh, t_obs[50], obs, figsize=(20, 15))

The Obs class also contains the method adjoint_derivative_obs_operator that computes the spatial map of \((d_cm)^*\cdot\delta c\), i.e. the adjoint of the derivative of the observation operator with respect to the state variable evaluated along \(\delta c\) at a given time. This method is needed to compute the source function of the adjoint model (see eq. (3) ).

It also contains the method onesensor_adjoint_derivative_obs_operator that is similar to the adjoint_derivative_obs_operator method but for the one-sensor observation operator, i.e. the observation operator obtained with only one sensor, and that can be used to compute the one-sensor-adjoint state.

The Obs class also have the attribute R_inv that contains \(\mathbf{R}^{-1}\), the inverse of the observation error covariance error. In the absence of dedicated data, this attribute is set to the identity matrix. However, it will be possible in the future to specify this covariance matrix.

Specifying the control variable and some prior knowledge about it using the control submodule

The control variable to be optimized (the quantity of emitted pheromones) and all the tools to model prior knowledge about it are contained in the control submodule.

The user should provide an object of the class Source that will be used as background/prior value of the control. When initialized, the control value is set to 0.

amp = 0.5
xc = 10.
yc = 10.
c = 10.
S_value_prior = np.zeros((msh.t_array.size, msh.y.size, msh.x.size))
xx, yy= np.meshgrid(msh.x, msh.y)
S_value_prior[:, np.maximum(np.abs(xx-xc), np.abs(yy-yc) ) < c] = amp

S_prior = Source(msh, S_value_prior, t=msh.t_array)

t_on = 1.
S_prior.at_current_time(t_on)
plot_colormap(msh,S_prior.value,'S_{prior}', 'g.m^{-2}.s^{-1}', cmap="jet",title=f't = {"{:.2f}".format(t_on)} s')

The user can include more prior knowledge of the control variable such as the population dynamic model. In such case, the user should specify the matrix or linear operator (by means of its matrix-vector product and of its transpose matrix-vector product) used to compute the residual of the population dynamic model.

Up to now, the submodule population_dynamic contains a stationnary population dynamique model \(\partial_t p = 0\) and a population dynamique model with linear reaction term \(\partial_t p = \gamma p\) modelling that a given proportion \(\gamma\) of the population stops to emit pheromone at each time step. The user can also easily specify its own population dynamic model using the different differential operators included in the toolbox (see the different notebooks to use differential operators).

In the following, for sake of simplicity, a population dynamic model with linear reaction term is used with a null reaction term, which is then equivalent to a stationary population dynamic.

death_rate = np.zeros((msh.t_array.size, msh.y.size, msh.x.size))
population_dynamic_model = PopulationDynamicModel(msh, death_rate=death_rate)

ctrl = Control(S_prior, msh, population_dynamique_model=population_dynamic_model)

For sake of illustration, the target source term that has been used to generate the data is stored in a dedicated object of the Controlclass.

ctrl_target = Control(S_target, msh)
ctrl_target.value = ctrl_target.background_value

Specifying the adjoint model and its solver using the adjoint_convection_diffusion_2D submodule

The adjoint model and its solver are here contained in the AdjointDiffusionConvectionReaction2DEquation class of the submodule adjoint_convection_diffusion_2D.

Similarly to the direct model, building an object of this class requires to specify the wind velocity field, the diffusion tensor and the loss coefficient.

solver = 'implicit with stationnary matrix inversion'
tol = 1e-10
adjoint_model =  AdjointDiffusionConvectionReaction2DEquation(U, K, deposition_coeff, msh, tol_inversion=tol, time_discretization=solver)
fname = 'inv_matrix_implicit_solver.npy'
adjoint_model.init_inverse_matrix(path_to_matrix=path_data, matrix_file_name=fname)

=== Load of the inverse of the matrix of the implicit part of the implicit with stationnary matrix inversion scheme ===

A special care has been put on the time discretization so that the discretized adjoint problem is effectively the adjoint (i.e. the matrix transposition) of the discretized direct problem.

The reader can find more details in the notebooks located in the folder ./test_numerical_scheme/ that introduces and test the numerical schemes and their convergence. These notebooks also include test of the adjoint operators and of the gradient.

Specifying the cost function using the cost submodule

The cost function, its different terms and its gradient are contained in the Cost class of the cost submodule. The user should provide an object of the class Obs that will be used to compute the observation term of the cost function, and an object of the class Control that will be used to compute the regularization terms based on biological prior knowledge on the control variable and updated during the optimization process

cost = Cost(msh, obs, ctrl)

The class Cost contains several useful attributes and methods that will be presented below. This class contains the attribute j_obs that is the observation term of the cost function, the attribute j_reg that is the dictionnary containing the different regularization terms, the attribute alpha that is the dictionnary of weight coefficient of the different regularization terms (see eq. (2)). The attribute implemented_regularization_types contains the list of keywords of the implemented regularization terms. The keywords of the attributes j_reg and alpha corresponds to keywords present in the attribute implemented_regularization_types.

print(cost.implemented_regularization_types)

['Tikhonov', 'Population dynamic', 'LASSO', 'time group LASSO', 'logarithm barrier']

The list of keyword of regularization terms to be considered and the associated dictionnary of weight coefficient can be provided as input of the instanciation of the class Cost or can be added and modified afterwards as in the following. If no regularization terms and its weigth coefficient are specified, than the inference will be made without regularization terms by default.

The class Cost also contains the methods obs_cost and reg_cost that computes resp. the observation and regularization terms and store the result in the associated attribute, as well as a method cost_and_gradient that brings together several methods and returns the evaluation of the cost function and its gradient from an array of shape (nb time point * nb cell in y axis * nb cell in x axis, 1) containing a quantity of pheromone emitted.

Inference of the quantity of emitted pheromones

The inference is performed using the minimize method of the class Cost. This method takes as inputs:

• the direct model that has been instanciated previously,

• the adjoint model that has been instanciated previously,

• the keyword of the optimization algorithm (see below).

The user can also provide the following optional intputs.

• The input options is the dictionnary of optional inputs of the optimization algorithm, such as tolerance of the termination criteria, the maximal number of iteration or the step size of the gradient descent.

• The input path_save is the string of the path to the folder in which the outputs will be saved. If not specified, the results of the optimization will not be save, just returned as outputs of the method.

• The input j_obs_threshold is the value of the threshold of the observation term of the cost function under which the optimization algorithm is stopped. This is particularly useful when some properties of the observation error are known and that we want to satisfy the Morozov’s principle. For instance, one do not want the observation term of the cost function to be less than the square of the norm of the noise.

• The input s_init enables to specify the initial value of the control variable. If not, the initial value is 0.

• The input restart_flag is a boolean flag that is False by default. It the flag is True, then the algorithm will load the outputs that has been previously saved in the folder given by the input path_save and start again where the algorithm stopped.

When the optimization algorithm reached the termination criteria, the method minimize returns:

• the direct model with the quantity of pheromone emitted resulting from the optimization algorithm as source term,

• the array containing the observation term of the cost function at each iteration of the optimization algorithm,

• the dictionnary containing the arrays containg the regularization terms of the cost function at each iteration of the optimization algorithm,

• the array containg the resulting quantity of pheromone emitted.

When the input path_save is specified, in addition to the outputs detailled above, the method saves:

• the weight coefficients of each regularization terms, in order to keep track of the test case,

• the array containg the resulting quantity of pheromone emitted in a file which name includes the number of the last iteration, to keep track of the different runs in case of restart,

• the array containing the step size at each iteration of the optimization algorithm, to keep track of the different runs in case of restart,

• the sum of the step size over the iterations of the optimization algorithm and the trajectory integrated quantity of pheromone emitted, to compute the trajectory-integrated one-sensor adjoint states. (This feature will be improved in a future update of the module).

Inference without regularization

In a first time, the inference is illustrated without regularization. In this situation, the user should use a differentiable optimization algorithm such as the gradient descent algorithm. This algorithm has been implemented in the submodule gradient_descent.

Warning: the values of hyperparameter used here, especially the values of the step size and the maximal number of iterations, are not ideal and more attention should be given in real cases but they allow fast computations and illustrative results in this notebook.

options = {
    "nit_max": 20,
    "step size": 5. * 1.0e+1
}

direct_model_no_reg, j_obs_no_reg, _, S_a_no_reg = cost.minimize(
    direct_model,
    adjoint_model,
    "gradient descent",
    options=options,
    path_save=path_output,
    j_obs_threshold = noise_level
)

initialization of the optimization process, computing the initial cost
Optimizing using the gradient descent algorithm
iteration 1: j = 1.167344e-02, j_obs = 1.167344e-02, ite comp time = 3.655e+00s -
iteration 2: j = 4.660842e-03, j_obs = 4.660842e-03, ite comp time = 6.494e+00s -
iteration 3: j = 3.549630e-03, j_obs = 3.549630e-03, ite comp time = 6.414e+00s -
iteration 4: j = 3.090900e-03, j_obs = 3.090900e-03, ite comp time = 6.297e+00s -
iteration 5: j = 2.786551e-03, j_obs = 2.786551e-03, ite comp time = 5.623e+00s -
iteration 6: j = 2.553989e-03, j_obs = 2.553989e-03, ite comp time = 5.277e+00s -
iteration 7: j = 2.365947e-03, j_obs = 2.365947e-03, ite comp time = 5.333e+00s -
iteration 8: j = 2.208759e-03, j_obs = 2.208759e-03, ite comp time = 5.149e+00s -
iteration 9: j = 2.074275e-03, j_obs = 2.074275e-03, ite comp time = 5.373e+00s -
iteration 10: j = 1.957173e-03, j_obs = 1.957173e-03, ite comp time = 5.246e+00s -
iteration 11: j = 1.853772e-03, j_obs = 1.853772e-03, ite comp time = 5.280e+00s -
iteration 12: j = 1.761407e-03, j_obs = 1.761407e-03, ite comp time = 5.952e+00s -
iteration 13: j = 1.678092e-03, j_obs = 1.678092e-03, ite comp time = 6.220e+00s -
iteration 14: j = 1.602304e-03, j_obs = 1.602304e-03, ite comp time = 6.819e+00s -
iteration 15: j = 1.532857e-03, j_obs = 1.532857e-03, ite comp time = 6.523e+00s -
iteration 16: j = 1.468809e-03, j_obs = 1.468809e-03, ite comp time = 6.525e+00s -
iteration 17: j = 1.409406e-03, j_obs = 1.409406e-03, ite comp time = 6.429e+00s -
iteration 18: j = 1.354034e-03, j_obs = 1.354034e-03, ite comp time = 6.441e+00s -
iteration 19: j = 1.302192e-03, j_obs = 1.302192e-03, ite comp time = 6.730e+00s -
iteration 20: j = 1.253466e-03, j_obs = 1.253466e-03, ite comp time = 6.252e+00s -

final iteration:
number of iteration = 20
j = 0.10737760430141632 % of the initial cost
j_obs = 0.10737760430141632 % of the initial cost

The resulting optimal value of the quantity of emitted pheromones can be plotted at the a given time using the plot_ctrl method of the utils module. The cost function and its different terms can also be plotted using the plot_cost method of the same module (the empty inputs correspond to the regularization terms, see below).

plot_ctrl(
    msh,
    msh.t_array[1],
    cost.ctrl,
    ctrl_target=ctrl_target.value,
    cmap="jet",
    label='s_o',
    obs=cost.obs,
    figsize=(20, 15)
)
plot_ctrl(
    msh,
    msh.t_array[-1],
    cost.ctrl,
    ctrl_target=ctrl_target.value,
    cmap="jet",
    label='s_o',
    obs=cost.obs,
    figsize=(20, 15)
)

/home/slabarthe/Documents/SVN_GIT/PHEROSENSOR/pherosensor-toolbox/src/utils/plot_ctrl.py:106: UserWarning: No contour levels were found within the data range.
  contour = ax.contour(msh.x, msh.y, ctrl_target_current_t, colors='k', linestyles='dashed', levels=0, linewidths=4)

plot_cost(
    j_obs_no_reg,
    {},
    noise_level,
    {},
    {},
    figsize=(20, 15)
)

Inference with a population dynamic-informed regularization term

Now, the inference is illustrated with the population dynamic-informed regularization.

One can add a regularization term to the cost function by using the add_reg method of the object of the class Cost. It takes as input the keyword of the regularization term and the value of the weight coefficient associated. If the regularization term is already in the cost function, this method will only update the value of the weight coefficient.

cost.add_reg("Population dynamic", 5. * 1e-2)

As the population dynamic-informed regularization term is differentiable, the gradient descent algorithm is used.

One may want to start from where the optimization stopped in the case without regularization (which is referred to as “hot start” in the litterature). To do so, the input ‘restart_flag` is set to True.

options = {
    "nit_max": 100,
    "step size": 5 * 1.0e-2
}

direct_model_pop_dyn_reg, j_obs_pop_dyn_reg, j_reg_pop_dyn_reg, S_a_pop_dyn_reg = cost.minimize(
    direct_model,
    adjoint_model,
    "gradient descent",
    options=options,
    path_save=path_output,
    restart_flag=True,
    j_obs_threshold = noise_level
)

restarting the optimization process from the data contained in the dir /home/slabarthe/Documents/SVN_GIT/PHEROSENSOR/pherosensor-toolbox/tutorials/output
Warning: the regularization Population dynamic term has not been found in the given directory. It is assumed that this regularization term was not used previously and thus is set to 0.
Optimizing using the gradient descent algorithm
iteration 21: j = 1.639804e-02, j_obs = 1.207513e-03, ite comp time = 6.609e+00s -
iteration 22: j = 9.358158e-03, j_obs = 1.208904e-03, ite comp time = 6.403e+00s -
iteration 23: j = 7.672829e-03, j_obs = 1.209167e-03, ite comp time = 6.625e+00s -
iteration 24: j = 6.808572e-03, j_obs = 1.209675e-03, ite comp time = 6.597e+00s -
iteration 25: j = 6.243985e-03, j_obs = 1.209887e-03, ite comp time = 6.350e+00s -
iteration 26: j = 5.829556e-03, j_obs = 1.210224e-03, ite comp time = 6.179e+00s -
iteration 27: j = 5.503943e-03, j_obs = 1.210407e-03, ite comp time = 5.965e+00s -
iteration 28: j = 5.236778e-03, j_obs = 1.210670e-03, ite comp time = 6.208e+00s -
iteration 29: j = 5.010742e-03, j_obs = 1.210833e-03, ite comp time = 6.005e+00s -
iteration 30: j = 4.815256e-03, j_obs = 1.211054e-03, ite comp time = 5.956e+00s -
iteration 31: j = 4.643237e-03, j_obs = 1.211202e-03, ite comp time = 5.847e+00s -
iteration 32: j = 4.489867e-03, j_obs = 1.211395e-03, ite comp time = 6.055e+00s -
iteration 33: j = 4.351593e-03, j_obs = 1.211532e-03, ite comp time = 6.156e+00s -
iteration 34: j = 4.225846e-03, j_obs = 1.211705e-03, ite comp time = 6.253e+00s -
iteration 35: j = 4.110601e-03, j_obs = 1.211832e-03, ite comp time = 6.062e+00s -
iteration 36: j = 4.004340e-03, j_obs = 1.211990e-03, ite comp time = 6.161e+00s -
iteration 37: j = 3.905805e-03, j_obs = 1.212109e-03, ite comp time = 6.047e+00s -
iteration 38: j = 3.814033e-03, j_obs = 1.212254e-03, ite comp time = 6.016e+00s -
iteration 39: j = 3.728191e-03, j_obs = 1.212366e-03, ite comp time = 6.040e+00s -
iteration 40: j = 3.647632e-03, j_obs = 1.212501e-03, ite comp time = 5.930e+00s -
iteration 41: j = 3.571777e-03, j_obs = 1.212607e-03, ite comp time = 5.921e+00s -
iteration 42: j = 3.500174e-03, j_obs = 1.212734e-03, ite comp time = 5.973e+00s -
iteration 43: j = 3.432403e-03, j_obs = 1.212834e-03, ite comp time = 6.018e+00s -
iteration 44: j = 3.368137e-03, j_obs = 1.212953e-03, ite comp time = 6.147e+00s -
iteration 45: j = 3.307059e-03, j_obs = 1.213049e-03, ite comp time = 5.965e+00s -
iteration 46: j = 3.248926e-03, j_obs = 1.213161e-03, ite comp time = 6.175e+00s -
iteration 47: j = 3.193494e-03, j_obs = 1.213252e-03, ite comp time = 6.408e+00s -
iteration 48: j = 3.140574e-03, j_obs = 1.213359e-03, ite comp time = 6.200e+00s -
iteration 49: j = 3.089975e-03, j_obs = 1.213446e-03, ite comp time = 6.268e+00s -
iteration 50: j = 3.041549e-03, j_obs = 1.213547e-03, ite comp time = 6.073e+00s -
iteration 51: j = 2.995140e-03, j_obs = 1.213630e-03, ite comp time = 6.012e+00s -
iteration 52: j = 2.950630e-03, j_obs = 1.213726e-03, ite comp time = 6.154e+00s -
iteration 53: j = 2.907892e-03, j_obs = 1.213805e-03, ite comp time = 6.078e+00s -
iteration 54: j = 2.866828e-03, j_obs = 1.213897e-03, ite comp time = 6.034e+00s -
iteration 55: j = 2.827331e-03, j_obs = 1.213973e-03, ite comp time = 6.107e+00s -
iteration 56: j = 2.789322e-03, j_obs = 1.214061e-03, ite comp time = 6.053e+00s -
iteration 57: j = 2.752711e-03, j_obs = 1.214134e-03, ite comp time = 1.217e+01s -
iteration 58: j = 2.717430e-03, j_obs = 1.214217e-03, ite comp time = 1.644e+01s -
iteration 59: j = 2.683402e-03, j_obs = 1.214287e-03, ite comp time = 1.548e+01s -
iteration 60: j = 2.650569e-03, j_obs = 1.214367e-03, ite comp time = 1.702e+01s -
iteration 61: j = 2.618866e-03, j_obs = 1.214434e-03, ite comp time = 1.623e+01s -
iteration 62: j = 2.588243e-03, j_obs = 1.214511e-03, ite comp time = 1.625e+01s -
iteration 63: j = 2.558642e-03, j_obs = 1.214575e-03, ite comp time = 1.736e+01s -
iteration 64: j = 2.530021e-03, j_obs = 1.214648e-03, ite comp time = 1.618e+01s -
iteration 65: j = 2.502329e-03, j_obs = 1.214710e-03, ite comp time = 1.597e+01s -
iteration 66: j = 2.475528e-03, j_obs = 1.214780e-03, ite comp time = 1.760e+01s -
iteration 67: j = 2.449575e-03, j_obs = 1.214839e-03, ite comp time = 1.871e+01s -
iteration 68: j = 2.424436e-03, j_obs = 1.214906e-03, ite comp time = 1.772e+01s -
iteration 69: j = 2.400072e-03, j_obs = 1.214963e-03, ite comp time = 1.663e+01s -
iteration 70: j = 2.376454e-03, j_obs = 1.215027e-03, ite comp time = 1.875e+01s -
iteration 71: j = 2.353546e-03, j_obs = 1.215081e-03, ite comp time = 1.816e+01s -
iteration 72: j = 2.331323e-03, j_obs = 1.215143e-03, ite comp time = 1.987e+01s -
iteration 73: j = 2.309754e-03, j_obs = 1.215195e-03, ite comp time = 1.843e+01s -
iteration 74: j = 2.288815e-03, j_obs = 1.215254e-03, ite comp time = 1.875e+01s -
iteration 75: j = 2.268479e-03, j_obs = 1.215304e-03, ite comp time = 1.750e+01s -
iteration 76: j = 2.248724e-03, j_obs = 1.215360e-03, ite comp time = 1.820e+01s -
iteration 77: j = 2.229525e-03, j_obs = 1.215408e-03, ite comp time = 1.797e+01s -
iteration 78: j = 2.210865e-03, j_obs = 1.215462e-03, ite comp time = 1.759e+01s -
iteration 79: j = 2.192718e-03, j_obs = 1.215507e-03, ite comp time = 1.860e+01s -
iteration 80: j = 2.175070e-03, j_obs = 1.215560e-03, ite comp time = 1.895e+01s -
iteration 81: j = 2.157899e-03, j_obs = 1.215603e-03, ite comp time = 1.905e+01s -
iteration 82: j = 2.141190e-03, j_obs = 1.215653e-03, ite comp time = 1.654e+01s -
iteration 83: j = 2.124924e-03, j_obs = 1.215694e-03, ite comp time = 1.599e+01s -
iteration 84: j = 2.109087e-03, j_obs = 1.215742e-03, ite comp time = 1.632e+01s -
iteration 85: j = 2.093662e-03, j_obs = 1.215781e-03, ite comp time = 1.585e+01s -
iteration 86: j = 2.078638e-03, j_obs = 1.215827e-03, ite comp time = 1.517e+01s -
iteration 87: j = 2.063996e-03, j_obs = 1.215865e-03, ite comp time = 1.691e+01s -
iteration 88: j = 2.049728e-03, j_obs = 1.215909e-03, ite comp time = 1.589e+01s -
iteration 89: j = 2.035817e-03, j_obs = 1.215944e-03, ite comp time = 1.634e+01s -
iteration 90: j = 2.022255e-03, j_obs = 1.215986e-03, ite comp time = 1.512e+01s -
iteration 91: j = 2.009026e-03, j_obs = 1.216020e-03, ite comp time = 1.650e+01s -
iteration 92: j = 1.996123e-03, j_obs = 1.216060e-03, ite comp time = 1.619e+01s -
iteration 93: j = 1.983532e-03, j_obs = 1.216093e-03, ite comp time = 1.619e+01s -
iteration 94: j = 1.971245e-03, j_obs = 1.216131e-03, ite comp time = 1.560e+01s -
iteration 95: j = 1.959251e-03, j_obs = 1.216162e-03, ite comp time = 1.680e+01s -
iteration 96: j = 1.947543e-03, j_obs = 1.216198e-03, ite comp time = 1.657e+01s -
iteration 97: j = 1.936108e-03, j_obs = 1.216227e-03, ite comp time = 1.686e+01s -
iteration 98: j = 1.924942e-03, j_obs = 1.216262e-03, ite comp time = 1.796e+01s -
iteration 99: j = 1.914032e-03, j_obs = 1.216290e-03, ite comp time = 1.752e+01s -
iteration 100: j = 1.903374e-03, j_obs = 1.216323e-03, ite comp time = 1.253e+03s -
iteration 101: j = 1.892958e-03, j_obs = 1.216349e-03, ite comp time = 1.628e+01s -
iteration 102: j = 1.882778e-03, j_obs = 1.216380e-03, ite comp time = 1.474e+01s -
iteration 103: j = 1.872825e-03, j_obs = 1.216405e-03, ite comp time = 8.805e+01s -
iteration 104: j = 1.863095e-03, j_obs = 1.216435e-03, ite comp time = 6.269e+00s -
iteration 105: j = 1.853578e-03, j_obs = 1.216458e-03, ite comp time = 5.585e+00s -
iteration 106: j = 1.844272e-03, j_obs = 1.216487e-03, ite comp time = 6.424e+00s -
iteration 107: j = 1.835166e-03, j_obs = 1.216509e-03, ite comp time = 6.427e+00s -
iteration 108: j = 1.826259e-03, j_obs = 1.216536e-03, ite comp time = 6.527e+00s -
iteration 109: j = 1.817541e-03, j_obs = 1.216556e-03, ite comp time = 6.285e+00s -
iteration 110: j = 1.809010e-03, j_obs = 1.216582e-03, ite comp time = 6.368e+00s -
iteration 111: j = 1.800658e-03, j_obs = 1.216601e-03, ite comp time = 6.228e+00s -
iteration 112: j = 1.792483e-03, j_obs = 1.216625e-03, ite comp time = 6.170e+00s -
iteration 113: j = 1.784477e-03, j_obs = 1.216643e-03, ite comp time = 6.180e+00s -
iteration 114: j = 1.776637e-03, j_obs = 1.216666e-03, ite comp time = 6.321e+00s -
iteration 115: j = 1.768958e-03, j_obs = 1.216683e-03, ite comp time = 6.363e+00s -
iteration 116: j = 1.761436e-03, j_obs = 1.216704e-03, ite comp time = 6.830e+00s -
iteration 117: j = 1.754066e-03, j_obs = 1.216720e-03, ite comp time = 6.875e+00s -
iteration 118: j = 1.746845e-03, j_obs = 1.216740e-03, ite comp time = 7.198e+00s -
iteration 119: j = 1.739767e-03, j_obs = 1.216754e-03, ite comp time = 7.369e+00s -
iteration 120: j = 1.732832e-03, j_obs = 1.216773e-03, ite comp time = 6.777e+00s -

final iteration:
number of iteration = 120
j = 0.1484422225954371 % of the initial cost
j_obs = 0.10423431620529669 % of the initial cost
j_Population dynamic = 0.04420790639014042 % of the initial cost

ctrl_pop_dyn_reg = Control(direct_model_pop_dyn_reg.S, msh)
ctrl_pop_dyn_reg.value = ctrl_pop_dyn_reg.background_value

plot_ctrl(
    msh,
    msh.t_array[1],
    ctrl_pop_dyn_reg,
    ctrl_target=ctrl_target.value,
    cmap="jet",
    label='s_o',
    obs=cost.obs,
    figsize=(20, 15)
)
plot_ctrl(
    msh,
    msh.t_array[-1],
    ctrl_pop_dyn_reg,
    ctrl_target=ctrl_target.value,
    cmap="jet",
    label='s_o',
    obs=cost.obs,
    figsize=(20, 15)
)

/home/slabarthe/Documents/SVN_GIT/PHEROSENSOR/pherosensor-toolbox/src/utils/plot_ctrl.py:106: UserWarning: No contour levels were found within the data range.
  contour = ax.contour(msh.x, msh.y, ctrl_target_current_t, colors='k', linestyles='dashed', levels=0, linewidths=4)

plot_cost(
    j_obs_pop_dyn_reg,
    j_reg_pop_dyn_reg,
    noise_level,
    {'Population dynamic': 'PD'},
    {'Population dynamic': 'g'},
    figsize=(20, 15)
)

Let us note that if the optimization is started again and no regularization terms are found in the provided folder, the regularization terms are set to 0.

Inference with the LASSO regularization term

Now, the inference is illustrated with both the population dynamic-informed and the LASSO regularizations.

cost.add_reg('LASSO',  10 * 1e-4)

In this case, the weight coefficient of the population dynamic-informed regularization term is modified using the modify_weight_reg method of the object of the class Cost. This method takes as input the keyword of the regularization term and the new value of the weight coefficient. If the regularization term was not already considered, then this method acts as the add_reg method. If the new weight coefficient is 0, then the regularization term will be removed.

cost.modify_weight_reg("Population dynamic", 2. * 1e-4)

As the cost function is no longer differentiable due to the LASSO regularization term, the user should use a convexe non-differentiable optimization algorithm such as the proximal gradient algorithm. This algorithm has been implemented in the submodule proximal_gradient. This type of algorithm should also be used when the group-LASSO regularization term is used.

options = {
    "nit_max": 40,
    "step size": 10 * 1.0e-1
}

direct_model_LASSO_reg, j_obs_LASSO_reg, j_reg_LASSO_reg, S_a_LASSO_reg = cost.minimize(
    direct_model,
    adjoint_model,
    "proximal gradient",
    options=options,
    path_save=path_output,
    restart_flag=True,
    j_obs_threshold = noise_level
)

restarting the optimization process from the data contained in the dir /home/slabarthe/Documents/SVN_GIT/PHEROSENSOR/pherosensor-toolbox/tutorials/output
Warning: the regularization LASSO term has not been found in the given directory. It is assumed that this regularization term was not used previously and thus is set to 0.
Optimizing using the proximal gradient algorithm
iteration 121: j = 9.183761e-02, j_obs = 1.216786e-03, ite comp time = 1.498e+01s -
iteration 122: j = 9.125290e-02, j_obs = 1.216103e-03, ite comp time = 1.456e+01s -
iteration 123: j = 9.069233e-02, j_obs = 1.215469e-03, ite comp time = 1.552e+01s -
iteration 124: j = 9.015351e-02, j_obs = 1.214884e-03, ite comp time = 1.513e+01s -
iteration 125: j = 8.963486e-02, j_obs = 1.214345e-03, ite comp time = 1.576e+01s -
iteration 126: j = 8.913538e-02, j_obs = 1.213851e-03, ite comp time = 1.620e+01s -
iteration 127: j = 8.865464e-02, j_obs = 1.213400e-03, ite comp time = 1.629e+01s -
iteration 128: j = 8.819243e-02, j_obs = 1.212991e-03, ite comp time = 1.346e+01s -
iteration 129: j = 8.774803e-02, j_obs = 1.212623e-03, ite comp time = 1.277e+01s -
iteration 130: j = 8.732052e-02, j_obs = 1.212291e-03, ite comp time = 1.256e+01s -
iteration 131: j = 8.690920e-02, j_obs = 1.211995e-03, ite comp time = 1.208e+01s -
iteration 132: j = 8.651345e-02, j_obs = 1.211733e-03, ite comp time = 1.258e+01s -
iteration 133: j = 8.613257e-02, j_obs = 1.211502e-03, ite comp time = 1.321e+01s -
iteration 134: j = 8.576591e-02, j_obs = 1.211301e-03, ite comp time = 1.217e+01s -
iteration 135: j = 8.541292e-02, j_obs = 1.211129e-03, ite comp time = 1.265e+01s -
iteration 136: j = 8.507305e-02, j_obs = 1.210984e-03, ite comp time = 1.329e+01s -
iteration 137: j = 8.474563e-02, j_obs = 1.210865e-03, ite comp time = 1.222e+01s -
iteration 138: j = 8.443013e-02, j_obs = 1.210771e-03, ite comp time = 1.289e+01s -
iteration 139: j = 8.412587e-02, j_obs = 1.210700e-03, ite comp time = 1.222e+01s -
iteration 140: j = 8.383223e-02, j_obs = 1.210651e-03, ite comp time = 1.213e+01s -
iteration 141: j = 8.354878e-02, j_obs = 1.210623e-03, ite comp time = 1.267e+01s -
iteration 142: j = 8.327504e-02, j_obs = 1.210616e-03, ite comp time = 1.237e+01s -
iteration 143: j = 8.301037e-02, j_obs = 1.210627e-03, ite comp time = 1.180e+01s -
iteration 144: j = 8.275438e-02, j_obs = 1.210656e-03, ite comp time = 1.299e+01s -
iteration 145: j = 8.250657e-02, j_obs = 1.210702e-03, ite comp time = 1.309e+01s -
iteration 146: j = 8.226660e-02, j_obs = 1.210764e-03, ite comp time = 1.270e+01s -
iteration 147: j = 8.203406e-02, j_obs = 1.210840e-03, ite comp time = 1.405e+01s -
iteration 148: j = 8.180871e-02, j_obs = 1.210929e-03, ite comp time = 1.319e+01s -
iteration 149: j = 8.159025e-02, j_obs = 1.211031e-03, ite comp time = 1.328e+01s -
iteration 150: j = 8.137838e-02, j_obs = 1.211144e-03, ite comp time = 1.293e+01s -
iteration 151: j = 8.117296e-02, j_obs = 1.211267e-03, ite comp time = 1.359e+01s -
iteration 152: j = 8.097372e-02, j_obs = 1.211399e-03, ite comp time = 1.548e+01s -
iteration 153: j = 8.078034e-02, j_obs = 1.211539e-03, ite comp time = 1.506e+01s -
iteration 154: j = 8.059263e-02, j_obs = 1.211687e-03, ite comp time = 1.676e+01s -
iteration 155: j = 8.041036e-02, j_obs = 1.211841e-03, ite comp time = 1.535e+01s -
iteration 156: j = 8.023325e-02, j_obs = 1.212000e-03, ite comp time = 1.546e+01s -
iteration 157: j = 8.006114e-02, j_obs = 1.212165e-03, ite comp time = 1.495e+01s -
iteration 158: j = 7.989385e-02, j_obs = 1.212333e-03, ite comp time = 1.602e+01s -
iteration 159: j = 7.973123e-02, j_obs = 1.212505e-03, ite comp time = 1.558e+01s -
iteration 160: j = 7.957308e-02, j_obs = 1.212679e-03, ite comp time = 1.667e+01s -

final iteration:
number of iteration = 160
j = 6.81659129779809 % of the initial cost
j_obs = 0.10388358493708338 % of the initial cost
j_Population dynamic = 0.0001659359187830108 % of the initial cost
j_LASSO = 6.712541776942223 % of the initial cost

ctrl_LASSO = Control(direct_model_LASSO_reg.S, msh)
ctrl_LASSO.value = ctrl_LASSO.background_value
plot_ctrl(
    msh,
    msh.t_array[1],
    ctrl_LASSO,
    ctrl_target=ctrl_target.value,
    cmap="jet",
    label='s_o',
    obs=cost.obs,
    figsize=(20, 15)
)
plot_ctrl(
    msh,
    msh.t_array[-1],
    ctrl_LASSO,
    ctrl_target=ctrl_target.value,
    cmap="jet",
    label='s_o',
    obs=cost.obs,
    figsize=(20, 15)
)

/home/slabarthe/Documents/SVN_GIT/PHEROSENSOR/pherosensor-toolbox/src/utils/plot_ctrl.py:106: UserWarning: No contour levels were found within the data range.
  contour = ax.contour(msh.x, msh.y, ctrl_target_current_t, colors='k', linestyles='dashed', levels=0, linewidths=4)

Inference with the Tikhonov regularization term

Now, the inference is illustrated with both the population dynamic-informed and the Tikhonov regularizations.

cost.add_reg('Tikhonov',5 * 1e-5)

One can remove a regularization term from the cost function using the remove_reg method of the object of the class Cost. This method just take the keyword of the regularization term to remove.

cost.remove_reg("LASSO")

As the LASSO regularization was removed and as the Tikhonov regularization term is differentiable, the cost function is differentiable and thus, the gradient descent algorithm will be used.

options = {
    "nit_max": 40,
    "step size": 15 * 1.0e-1
}

direct_model_Tikhonov_reg, j_obs_Tikhonov_reg, j_reg_Tikhonov_reg, S_a_Tikhonov_reg = cost.minimize(
    direct_model,
    adjoint_model,
    "gradient descent",
    options=options,
    path_save=path_output,
    restart_flag=True,
    j_obs_threshold = noise_level
)

restarting the optimization process from the data contained in the dir /home/slabarthe/Documents/SVN_GIT/PHEROSENSOR/pherosensor-toolbox/tutorials/output
Warning: the regularization Tikhonov term has not been found in the given directory. It is assumed that this regularization term was not used previously and thus is set to 0.
Optimizing using the gradient descent algorithm
iteration 161: j = 9.809640e-02, j_obs = 1.212855e-03, ite comp time = 7.005e+00s -
iteration 162: j = 9.805680e-02, j_obs = 1.206557e-03, ite comp time = 6.695e+00s -
iteration 163: j = 9.801895e-02, j_obs = 1.201276e-03, ite comp time = 6.908e+00s -
iteration 164: j = 9.798274e-02, j_obs = 1.196922e-03, ite comp time = 7.180e+00s -
iteration 165: j = 9.794805e-02, j_obs = 1.193410e-03, ite comp time = 6.803e+00s -
iteration 166: j = 9.791479e-02, j_obs = 1.190662e-03, ite comp time = 6.769e+00s -
iteration 167: j = 9.788286e-02, j_obs = 1.188607e-03, ite comp time = 6.745e+00s -
iteration 168: j = 9.785217e-02, j_obs = 1.187181e-03, ite comp time = 7.210e+00s -
iteration 169: j = 9.782265e-02, j_obs = 1.186322e-03, ite comp time = 7.362e+00s -
iteration 170: j = 9.779423e-02, j_obs = 1.185978e-03, ite comp time = 7.062e+00s -
iteration 171: j = 9.776682e-02, j_obs = 1.186096e-03, ite comp time = 6.797e+00s -
iteration 172: j = 9.774037e-02, j_obs = 1.186631e-03, ite comp time = 7.066e+00s -
iteration 173: j = 9.771482e-02, j_obs = 1.187542e-03, ite comp time = 7.111e+00s -
iteration 174: j = 9.769011e-02, j_obs = 1.188788e-03, ite comp time = 7.212e+00s -
iteration 175: j = 9.766619e-02, j_obs = 1.190336e-03, ite comp time = 7.219e+00s -
iteration 176: j = 9.764301e-02, j_obs = 1.192152e-03, ite comp time = 7.208e+00s -
iteration 177: j = 9.762053e-02, j_obs = 1.194207e-03, ite comp time = 6.509e+00s -
iteration 178: j = 9.759871e-02, j_obs = 1.196475e-03, ite comp time = 6.116e+00s -
iteration 179: j = 9.757750e-02, j_obs = 1.198930e-03, ite comp time = 6.041e+00s -
iteration 180: j = 9.755688e-02, j_obs = 1.201551e-03, ite comp time = 6.692e+00s -
iteration 181: j = 9.753680e-02, j_obs = 1.204316e-03, ite comp time = 6.005e+00s -
iteration 182: j = 9.751724e-02, j_obs = 1.207207e-03, ite comp time = 5.804e+00s -
iteration 183: j = 9.749817e-02, j_obs = 1.210208e-03, ite comp time = 6.274e+00s -
iteration 184: j = 9.747956e-02, j_obs = 1.213302e-03, ite comp time = 6.097e+00s -
iteration 185: j = 9.746138e-02, j_obs = 1.216476e-03, ite comp time = 6.232e+00s -
iteration 186: j = 9.744362e-02, j_obs = 1.219717e-03, ite comp time = 6.217e+00s -
iteration 187: j = 9.742624e-02, j_obs = 1.223013e-03, ite comp time = 5.758e+00s -
iteration 188: j = 9.740923e-02, j_obs = 1.226354e-03, ite comp time = 5.960e+00s -
iteration 189: j = 9.739258e-02, j_obs = 1.229731e-03, ite comp time = 5.683e+00s -
iteration 190: j = 9.737625e-02, j_obs = 1.233136e-03, ite comp time = 5.730e+00s -
iteration 191: j = 9.736025e-02, j_obs = 1.236559e-03, ite comp time = 5.984e+00s -
iteration 192: j = 9.734454e-02, j_obs = 1.239996e-03, ite comp time = 5.841e+00s -
iteration 193: j = 9.732911e-02, j_obs = 1.243438e-03, ite comp time = 5.850e+00s -
iteration 194: j = 9.731396e-02, j_obs = 1.246881e-03, ite comp time = 6.898e+00s -
iteration 195: j = 9.729906e-02, j_obs = 1.250320e-03, ite comp time = 5.560e+00s -
iteration 196: j = 9.728442e-02, j_obs = 1.253751e-03, ite comp time = 5.714e+00s -
iteration 197: j = 9.727000e-02, j_obs = 1.257168e-03, ite comp time = 5.829e+00s -
iteration 198: j = 9.725581e-02, j_obs = 1.260570e-03, ite comp time = 5.725e+00s -
iteration 199: j = 9.724184e-02, j_obs = 1.263952e-03, ite comp time = 5.738e+00s -
iteration 200: j = 9.722807e-02, j_obs = 1.267312e-03, ite comp time = 5.790e+00s -

final iteration:
number of iteration = 200
j = 8.328998044207413 % of the initial cost
j_obs = 0.10856372400338492 % of the initial cost
j_Population dynamic = 0.00017074318343414518 % of the initial cost
j_Tikhonov = 8.220263577020594 % of the initial cost

ctrl_Tikhonov = Control(direct_model_Tikhonov_reg.S, msh)
ctrl_Tikhonov.value = ctrl_Tikhonov.background_value
plot_ctrl(
    msh,
    msh.t_array[1],
    ctrl_Tikhonov,
    ctrl_target=ctrl_target.value,
    cmap="jet",
    label='s_o',
    obs=cost.obs,
    figsize=(20, 15)
)
plot_ctrl(
    msh,
    msh.t_array[-1],
    ctrl_Tikhonov,
    ctrl_target=ctrl_target.value,
    cmap="jet",
    label='s_o',
    obs=cost.obs,
    figsize=(20, 15)
)