{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "63d68a42",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import pherosensor\n",
    "\n",
    "from pheromone_dispersion.geom import MeshRect2D\n",
    "from pheromone_dispersion.velocity import Velocity\n",
    "from pheromone_dispersion.diffusion_tensor import DiffusionTensor\n",
    "from pheromone_dispersion.diffusion_operator import Diffusion\n",
    "from pheromone_dispersion.advection_operator import Advection, AdvectionAdjoint\n",
    "from pheromone_dispersion.reaction_operator import Reaction\n",
    "from source_localization.population_dynamique import PopulationDynamicModel\n",
    "from source_localization.population_dynamique import StationnaryPopulationDynamicModel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "7ff10728",
   "metadata": {},
   "outputs": [],
   "source": [
    "Lx = 50 # the length along the x-axis\n",
    "Ly = 60 # the length along the y-axis\n",
    "Delta_x = 0.4 # 0.4 the space step along the x-axis\n",
    "Delta_y = 0.4 # 0.4 the space step along the y-axis\n",
    "T_final = 5 # the final time of the simulation\n",
    "msh = MeshRect2D(Lx, Ly, Delta_x, Delta_y, T_final)\n",
    "msh.calc_dt_implicit_solver(.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a80eaeed",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "v1 = np.random.normal(0, 1, size=msh.x.size*msh.y.size)\n",
    "v2 = np.random.normal(0, 1, size=msh.x.size*msh.y.size)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76fea742",
   "metadata": {},
   "source": [
    "# Test of the adjoint operators \n",
    "\n",
    "This notebook aims to check numerically that the discretized adjoint operators are indeed the adjoint of the discretized direct operators. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "864e983f-17b8-4780-a550-435e21cf6894",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "x1 = np.random.normal(0, 1, size=msh.x.size*msh.y.size*(msh.t_array.size-1))\n",
    "x2 = np.random.normal(0, 1, size=msh.x.size*msh.y.size*msh.t_array.size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "00c8d8af",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "test of the adjoint operator of the population dynamic model operator:\n",
      "<v_1,Av_2> =  12503.694354331254\n",
      "<A*v_1,v_2> =  12503.694354331188\n",
      "|<v_1,Av_2>-<A*v_1,v_2>| =  6.548361852765083e-11\n"
     ]
    }
   ],
   "source": [
    "pop_dyn_model = PopulationDynamicModel(msh)\n",
    "print(\"test of the adjoint operator of the population dynamic model operator:\")\n",
    "print(\"<v_1,Av_2> = \",np.dot(x1,pop_dyn_model.matvec(x2)))\n",
    "print(\"<A*v_1,v_2> = \",np.dot(pop_dyn_model.rmatvec(x1), x2))\n",
    "print(\"|<v_1,Av_2>-<A*v_1,v_2>| = \", np.abs(np.dot(x1,pop_dyn_model.matvec(x2))-np.dot(pop_dyn_model.rmatvec(x1),x2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "cd692889-4695-483e-bada-20c4b26ab801",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "test of the adjoint operator of the stationnary population dynamic model operator:\n",
      "<v_1,Av_2> =  12455.477432781594\n",
      "<A*v_1,v_2> =  12455.47743278158\n",
      "|<v_1,Av_2>-<A*v_1,v_2>| =  1.4551915228366852e-11\n"
     ]
    }
   ],
   "source": [
    "stat_pop_dyn_model = StationnaryPopulationDynamicModel(msh)\n",
    "print(\"test of the adjoint operator of the stationnary population dynamic model operator:\")\n",
    "print(\"<v_1,Av_2> = \",np.dot(x1,stat_pop_dyn_model.matvec(x2)))\n",
    "print(\"<A*v_1,v_2> = \",np.dot(stat_pop_dyn_model.rmatvec(x1), x2))\n",
    "print(\"|<v_1,Av_2>-<A*v_1,v_2>| = \", np.abs(np.dot(x1,stat_pop_dyn_model.matvec(x2))-np.dot(stat_pop_dyn_model.rmatvec(x1),x2)))   "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b1fd7e3",
   "metadata": {},
   "source": [
    "## Test of the adjoint of the advection operator\n",
    "\n",
    "Recall that the advection operator is: $A:c(x,y)\\mapsto \\nabla\\cdot(U(x,y) c(x,y))~\\forall (x,y)\\in\\Omega$ such that $U(x,y)c(x,y)\\cdot \\vec{n} = 0~\\forall(x,y)\\in\\partial\\Omega_i$ with $\\partial\\Omega_i = \\{ (x,y)\\in\\partial\\Omega, U(x,y)\\cdot\\vec{n}<0\\}$.\\\n",
    "On the other hand, the adjoint operator of the advection operator is $A^*:c(x,y)\\mapsto -U(x,y)\\cdot\\nabla c(x,y)=-\\nabla\\cdot(U(x,y) c(x,y))+c(x,y)\\nabla\\cdot U(x,y)~\\forall (x,y)\\in\\Omega$ such that $U(x,y)c(x,y)\\cdot \\vec{n} = 0~\\forall(x,y)\\in\\partial\\Omega_o$ with $\\partial\\Omega_o = \\{ (x,y)\\in\\partial\\Omega, U(x,y)\\cdot\\vec{n}>0\\}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "ef95c7ff",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "U_hi = np.random.normal(0, 1, size=(msh.y_horizontal_interface.size,msh.x.size,2))\n",
    "U_vi = np.random.normal(0, 1, size=(msh.y.size,msh.x_vertical_interface.size,2))\n",
    "\n",
    "U = Velocity(msh, U_vi, U_hi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "cf57659f",
   "metadata": {},
   "outputs": [],
   "source": [
    "A_flux_term = Advection(U, msh)\n",
    "A_divU_term = Reaction(U.div, msh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "321294b1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "test of the adjoint operator of the advection operator:\n",
      "<v_1,Av_2> =  -676.6346071435416\n",
      "<A*v_1,v_2> =  -676.6346071435415\n",
      "|<v_1,Av_2>-<A*v_1,v_2>| =  1.1368683772161603e-13\n"
     ]
    }
   ],
   "source": [
    "print(\"test of the adjoint operator of the advection operator:\")\n",
    "print(\"<v_1,Av_2> = \",np.dot(v1,A_flux_term.matvec(v2)))\n",
    "print(\"<A*v_1,v_2> = \",np.dot(A_flux_term.rmatvec(v1)+A_divU_term.matvec(v1),v2))\n",
    "print(\"|<v_1,Av_2>-<A*v_1,v_2>| = \", np.abs(np.dot(v1,A_flux_term.matvec(v2))-np.dot(A_flux_term.rmatvec(v1)+A_divU_term.matvec(v1),v2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f3c694f",
   "metadata": {},
   "source": [
    "## Test of the adjoint of the diffusion operator\n",
    "\n",
    "Recall that the diffusion operator is: $D:c\\mapsto \\nabla\\cdot(\\mathbf{K}(x,y)\\nabla c(x,y))~\\forall (x,y)\\in\\Omega$ such that $K\\nabla c\\cdot \\vec{n} = 0~\\forall(x,y)\\in\\partial\\Omega$.\\\n",
    "This operator is: $D^*:c\\mapsto \\nabla\\cdot(\\mathbf{K}^T(x,y)\\nabla c(x,y))~\\forall (x,y)\\in\\Omega$ such that $K^T\\nabla c\\cdot \\vec{n} = 0~\\forall(x,y)\\in\\partial\\Omega$.\\\n",
    "Let us note that the diffusion operator is auto-adjoint ($D^*=D$) for a symmetric diffusion tensor ($\\mathbf{K}=\\mathbf{K}^T$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "c1e688de",
   "metadata": {},
   "outputs": [],
   "source": [
    "U_hi = np.ones((msh.y_horizontal_interface.size,msh.x.size,2))\n",
    "U_hi[:,:,1]=0\n",
    "U_vi = np.ones((msh.y.size,msh.x_vertical_interface.size,2))\n",
    "U_vi[:,:,1]=0\n",
    "U = Velocity(msh, U_vi, U_hi)\n",
    "\n",
    "K_u_t = 0.01\n",
    "K_u = 0.01\n",
    "K = DiffusionTensor(U, K_u, K_u_t)\n",
    "\n",
    "# change for a "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "6aaabc1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "D = Diffusion(K, msh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "3249df5b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "test of the adjoint operator of the diffusion operator\n",
      "<v_1,Dv_2> =  13.741038678354187\n",
      "<D*v_1,v_2> =  13.741038678354208\n",
      "|<v_1,Dv_2>-<D*v_1,v_2>| =  2.1316282072803006e-14\n"
     ]
    }
   ],
   "source": [
    "print(\"test of the adjoint operator of the diffusion operator\")\n",
    "print(\"<v_1,Dv_2> = \",np.dot(v1,D.matvec(v2)))\n",
    "print(\"<D*v_1,v_2> = \",np.dot(D.matvec(v1),v2))\n",
    "print(\"|<v_1,Dv_2>-<D*v_1,v_2>| = \", np.abs(np.dot(v1,D.matvec(v2))-np.dot(D.matvec(v1),v2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fafa200",
   "metadata": {},
   "source": [
    "## Test of the adjoint of the reaction operator\n",
    "\n",
    "Recall that the reaction operator is: $D:c\\mapsto\\tau_{loss}c$.\\\n",
    "This operator is auto-adjoint: $R^*=R$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "670c9549",
   "metadata": {},
   "outputs": [],
   "source": [
    "tau_loss = np.random.normal(0, 1, size=(msh.y.size,msh.x.size))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "0cecaaeb",
   "metadata": {},
   "outputs": [],
   "source": [
    "R = Reaction(tau_loss, msh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "b026834d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "test of the adjoint operator of the diffusion operator\n",
      "<v_1,Dv_2> =  151.0434490129216\n",
      "<D*v_1,v_2> =  151.0434490129216\n",
      "|<v_1,Dv_2>-<D*v_1,v_2>| =  0.0\n"
     ]
    }
   ],
   "source": [
    "print(\"test of the adjoint operator of the diffusion operator\")\n",
    "print(\"<v_1,Dv_2> = \",np.dot(v1,R.matvec(v2)))\n",
    "print(\"<D*v_1,v_2> = \",np.dot(R.matvec(v1),v2))\n",
    "print(\"|<v_1,Dv_2>-<D*v_1,v_2>| = \", np.abs(np.dot(v1,R.matvec(v2))-np.dot(R.matvec(v1),v2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0696c32e",
   "metadata": {},
   "source": [
    "## Test of the adjoint operators of the population dynamic model operators"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41f35d06-b48a-4d36-935d-6da5a56e8ce9",
   "metadata": {},
   "source": [
    "Recall that the population dynamic model considered in the toolbox until now is $\\partial_t s = \\gamma s$, and when $\\gamma=0$, we have a stationary population dynamic model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "3a58e1f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "v1 = np.random.normal(0, 1, size=msh.x.size*msh.y.size*(msh.t_array.size-1))\n",
    "v2 = np.random.normal(0, 1, size=msh.x.size*msh.y.size*msh.t_array.size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "0f07c242",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "test of the adjoint operator of the population dynamique\n",
      "<v_1,Dv_2> =  12503.694354331254\n",
      "<D*v_1,v_2> =  12503.694354331188\n",
      "|<v_1,Dv_2>-<D*v_1,v_2>| =  6.548361852765083e-11\n"
     ]
    }
   ],
   "source": [
    "pdm = PopulationDynamicModel(msh)\n",
    "\n",
    "print(\"test of the adjoint operator of the population dynamique\")\n",
    "print(\"<v_1,Dv_2> = \",np.dot(v1,pdm.matvec(v2)))\n",
    "print(\"<D*v_1,v_2> = \",np.dot(pdm.rmatvec(v1),v2))\n",
    "print(\"|<v_1,Dv_2>-<D*v_1,v_2>| = \", np.abs(np.dot(v1,pdm.matvec(v2))-np.dot(pdm.rmatvec(v1),v2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "2d4d495c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "test of the adjoint operator of the stationnary population dynamique\n",
      "<v_1,Dv_2> =  12455.477432781594\n",
      "<D*v_1,v_2> =  12455.47743278158\n",
      "|<v_1,Dv_2>-<D*v_1,v_2>| =  1.4551915228366852e-11\n"
     ]
    }
   ],
   "source": [
    "spdm = StationnaryPopulationDynamicModel(msh)\n",
    "\n",
    "print(\"test of the adjoint operator of the stationnary population dynamique\")\n",
    "print(\"<v_1,Dv_2> = \",np.dot(v1,spdm.matvec(v2)))\n",
    "print(\"<D*v_1,v_2> = \",np.dot(spdm.rmatvec(v1),v2))\n",
    "print(\"|<v_1,Dv_2>-<D*v_1,v_2>| = \", np.abs(np.dot(v1,spdm.matvec(v2))-np.dot(spdm.rmatvec(v1),v2)))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pherosensor-new",
   "language": "python",
   "name": "pherosensor-new"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}