06_Poisson_INCG/Poisson_INCG.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Coefficient field inversion in an elliptic partial differential equation\n",
    "\n",
    "We consider the estimation of a coefficient in an elliptic partial\n",
    "differential equation as a model problem. Depending on the\n",
    "interpretation of the unknowns and the type of measurements, this\n",
    "model problem arises, for instance, in inversion for groundwater flow\n",
    "or heat conductivity.  It can also be interpreted as finding a\n",
    "membrane with a certain spatially varying stiffness. Let\n",
    "$\\Omega\\subset\\mathbb{R}^n$, $n\\in\\{1,2,3\\}$ be an open, bounded\n",
    "domain and consider the following problem:\n",
    "\n",
    "$$\n",
    "\\min_{m} J(m):=\\frac{1}{2}\\int_\\Omega (u-u_d)^2\\, dx + \\frac{\\gamma}{2}\\int_\\Omega|\\nabla m|^2\\,dx,\n",
    "$$\n",
    "\n",
    "where $u$ is the solution of\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "\\quad -\\nabla\\cdot(\\exp(m)\\nabla u) &= f \\text{ in }\\Omega,\\\\\n",
    "u &= 0 \\text{ on }\\partial\\Omega.\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "Here $m\\in U_{ad}:=\\{m\\in H^1(\\Omega) \\bigcap L^{\\infty}(\\Omega)\\}$ the unknown coefficient field, $u_d$ denotes (possibly noisy) data, $f\\in H^{-1}(\\Omega)$ a given force, and $\\gamma\\ge 0$ the regularization parameter.\n",
    "\n",
    "### The variational (or weak) form of the state equation:\n",
    "\n",
    "Find $u\\in H_0^1(\\Omega)$ such that \n",
    "\n",
    "$$(\\exp(m)\\nabla u,\\nabla v) - (f,v) = 0, \\text{ for all } v\\in H_0^1(\\Omega),$$\n",
    "\n",
    "where $H_0^1(\\Omega)$ is the space of functions vanishing on $\\partial\\Omega$ with square integrable derivatives. Here, $(\\cdot\\,,\\cdot)$ denotes the $L^2$-inner product, i.e, for scalar functions $u,v \\in L^2(\\Omega)$  we denote \n",
    "\n",
    "$$(u,v) := \\int_\\Omega u(x) v(x) \\,dx.$$\n",
    "\n",
    "### Gradient evaluation:\n",
    "\n",
    "The Lagrangian functional $\\mathscr{L}:H_0^1(\\Omega)\\times H^1(\\Omega)\\times H_0^1(\\Omega)\\rightarrow \\mathbb{R}$ is given by\n",
    "\n",
    "$$\n",
    "\\mathscr{L}(u,m,p):= \\frac{1}{2}(u-u_d,u-u_d) +\n",
    "\\frac{\\gamma}{2}(\\nabla m, \\nabla m) +  (\\exp(m)\\nabla u,\\nabla p) - (f,p).\n",
    "$$\n",
    "\n",
    "Then the gradient of the cost functional $\\mathcal{J}(m)$ with respect to the parameter $m$ is\n",
    "\n",
    "$$\n",
    "    \\mathcal{G}(m)(\\tilde m) := \\gamma(\\nabla m, \\nabla \\tilde{m}) +\n",
    "     (\\tilde{m}\\exp(m)\\nabla u, \\nabla p) \\quad \\forall \\tilde{m} \\in H^1(\\Omega),\n",
    "$$\n",
    "\n",
    "where $u \\in H_0^1(\\Omega)$ is the solution of the forward problem,\n",
    "\n",
    "$$ \\mathscr{L}_p(u,m,p)(\\tilde{p})  := (\\exp(m)\\nabla u, \\nabla \\tilde{p}) - (f,\\tilde{p}) = 0\n",
    "\\quad \\forall \\tilde{p} \\in H_0^1(\\Omega), $$\n",
    "\n",
    "and $p \\in H_0^1(\\Omega)$ is the solution of the adjoint problem,\n",
    "\n",
    "$$ \\mathscr{L}_u(u,m,p)(\\tilde{u}) := (\\exp(m)\\nabla p, \\nabla \\tilde{u}) + (u-u_d,\\tilde{u}) = 0\n",
    "\\quad \\forall \\tilde{u} \\in H_0^1(\\Omega).$$\n",
    "\n",
    "### Hessian action:\n",
    "\n",
    "To evaluate the action $\\mathcal{H}(m)(\\hat{m})$ of the Hessian is a given direction $\\hat{m}$ , we consider variations of the meta-Lagrangian functional\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mathscr{L}^H(u,m,p; \\hat{u}, \\hat{m}, \\hat{p}) := & {} & {} \\\\\n",
    "{} & \\gamma(\\nabla m, \\nabla \\tilde{m}) + (\\tilde{m}\\exp(m)\\nabla u, \\nabla p) & \\text{gradient}\\\\\n",
    "{} & + (\\exp(m)\\nabla u, \\nabla \\hat{p}) - (f,\\hat{p}) & \\text{forward eq}\\\\\n",
    "{} & + (\\exp(m)\\nabla p, \\nabla \\hat{u}) + (u-u_d,\\hat{u}) & \\text{adjoint eq}.\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Then action of the Hessian is a given direction $\\hat{m}$ is\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "(\\tilde{m}, \\mathcal{H}(m)(\\hat{m}) ) & := \\mathscr{L}^H_m(u,m,p; \\hat{u}, \\hat{m}, \\hat{p})(\\tilde{m}) \\\\\n",
    "{} & =\n",
    "(\\tilde{m} \\exp(m) \\nabla \\hat{u}, \\nabla{p}) + \\gamma (\\nabla \\hat{m}, \\nabla \\tilde{m}) + (\\tilde{m} \\hat{m} \\exp(m)\\nabla u, \\nabla p) + (\\tilde{m}\\exp(m) \\nabla u, \\nabla \\hat{p}) \\quad \\forall \\tilde{m} \\in H^1(\\Omega),\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "where \n",
    "\n",
    "- $u\\in H^1_0(\\Omega)$ and $p \\in H^1_0(\\Omega)$ are the solution of the forward and adjoint problem, respectively;\n",
    "\n",
    "- $\\hat{u} \\in H^1_0(\\Omega)$ is the solution of the incremental forward problem,\n",
    "\n",
    "$$\n",
    "\\mathscr{L}^H_p(u,m,p; \\hat{u}, \\hat{m}, \\hat{p})(\\tilde{p}) := (\\exp(m) \\nabla \\hat{u}, \\nabla \\tilde{p}) + (\\hat{m} \\exp(m) \\nabla u, \\nabla \\tilde p) = 0 \\quad \\forall \\tilde{p} \\in H_0^1(\\Omega);\n",
    "$$\n",
    "\n",
    "\n",
    "- and $\\hat{p} \\in H^1_0(\\Omega)$ is the solution of the incremental adjoint problem,\n",
    "$$\n",
    "\\mathscr{L}^H_u(u,m,p; \\hat{u}, \\hat{m}, \\hat{p})(\\tilde{u}) := (\\hat{u}, \\tilde{u}) + (\\hat{m} \\exp(m)\\nabla p, \\nabla \\tilde{u}) + (\\exp(m) \\nabla \\tilde u, \\nabla \\hat{p}) = 0 \\quad \\forall \\tilde{u} \\in H_0^1(\\Omega).\n",
    "$$\n",
    "\n",
    "### Inexact Newton-CG:\n",
    "\n",
    "Written in abstract form, the Newton Method computes an update direction $\\hat{m}_k$ by solving the linear system \n",
    "\n",
    "$$\n",
    "(\\tilde{m}, \\mathcal{H}(m_k)(\\hat{m}_k) ) = -\\mathcal{G}(m_k)(\\tilde m) \\quad \\forall \\tilde{m} \\in H^1(\\Omega),\n",
    "$$\n",
    "\n",
    "where the evaluation of the gradient $\\mathcal{G}(m_k)$ involve the solution $u_k$ and $p_k$ of the forward and adjoint problem (respectively) for $m = m_k$.\n",
    "Similarly, the Hessian action $\\mathcal{H}(m_k)(\\hat{m}_k)$ requires to additional solve the incremental forward and adjoint problems.\n",
    "\n",
    "### Discrete Newton system:\n",
    "$\n",
    "\\def\\tu{\\tilde u}\n",
    "\\def\\tm{\\tilde m}\n",
    "\\def\\tp{\\tilde p}\n",
    "\\def\\hu{\\hat u}\n",
    "\\def\\hp{\\hat p}\n",
    "\\def\\hm{\\hat m}\n",
    "$\n",
    "\n",
    "$\n",
    "\\def\\bu{{\\bf u}}\n",
    "\\def\\bm{{\\bf m}}\n",
    "\\def\\bp{{\\bf p}}\n",
    "\\def\\btu{{\\bf \\tilde u}}\n",
    "\\def\\btm{{\\bf \\tilde m}}\n",
    "\\def\\btp{{\\bf \\tilde p}}\n",
    "\\def\\bhu{{\\bf \\hat u}}\n",
    "\\def\\bhm{{\\bf \\hat m}}\n",
    "\\def\\bhp{{\\bf \\hat p}}\n",
    "\\def\\bg{{\\bf g}}\n",
    "$\n",
    "\n",
    "$\n",
    "\\def\\bA{{\\bf A}}\n",
    "\\def\\bC{{\\bf C}}\n",
    "\\def\\bH{{\\bf H}}\n",
    "\\def\\bR{{\\bf R}}\n",
    "\\def\\bW{{\\bf W}}\n",
    "$\n",
    "\n",
    "Let us denote the vectors corresponding to the discretization of the functions $u_k, m_k, p_k$ by $\\bu_k, \\bm_k, \\bp_k$ and of the functions $\\hu_k, \\hm_k, \\hp_k$ by $\\bhu_k, \\bhm_k,\\bhp_k$.\n",
    "\n",
    "Then, the discretization of the above system is given by the following symmetric linear system:\n",
    "\n",
    "$$\n",
    "  \\bH_k \\, \\bhm_k = -\\bg_k.\n",
    "$$\n",
    "\n",
    "The gradient $\\bg_k$ is computed using the following three steps\n",
    "\n",
    "- Given $\\bm_k$ we solve the forward problem\n",
    "\n",
    "$$ \\bA_k \\bu_k = {\\bf f}, $$\n",
    "\n",
    "where $\\bA_k \\bu_k$ stems from the discretization $(\\exp(m_k)\\nabla u_k, \\nabla \\tilde{p})$, and ${\\bf f}$ stands for the discretization of the right hand side $f$.\n",
    "\n",
    "- Given $\\bm_k$ and $\\bu_k$ solve the adjoint problem\n",
    "\n",
    "$$ \\bA_k^T \\bp_k = - \\bW_{\\scriptsize\\mbox{uu}}\\,(\\bu_k-\\bu_d) $$\n",
    "\n",
    "where $\\bA_k^T \\bp_k$ stems from the discretization of $(\\exp(m_k)\\nabla \\tilde{u}, \\nabla p_k)$, $\\bW_{\\scriptsize\\mbox{uu}}$ is the mass matrix corresponding to the $L^2$ inner product in the state space, and $\\bu_d$ stems from the data.\n",
    "\n",
    "- Define the gradient \n",
    "\n",
    "$$ \\bg_k = \\bR \\bm_k + \\bC_k^T \\bp_k, $$\n",
    "\n",
    "where $\\bR$ is the matrix stemming from discretization of the regularization operator $\\gamma ( \\nabla \\hat{m}, \\nabla \\tilde{m})$, and $\\bC_k$ stems from discretization of the term $(\\tilde{m}\\exp(m_k)\\nabla u_k, \\nabla p_k)$.\n",
    "\n",
    "Similarly the action of the Hessian $\\bH_k \\, \\bhm_k$ in a direction $\\bhm_k$ (by using the CG algorithm we only need the action of $\\bH_k$ to solve the Newton step) is given by\n",
    "\n",
    "- Solve the incremental forward problem\n",
    "\n",
    "$$ \\bA_k \\bhu_k = -\\bC_k \\bhm_k, $$\n",
    "\n",
    "where $\\bC_k \\bm_k$ stems from discretization of $(\\hat{m} \\exp(m_k) \\nabla u_k, \\nabla \\tilde p)$.\n",
    "\n",
    "- Solve the incremental adjoint problem\n",
    "\n",
    "$$ \\bA_k^T \\bhp_k = -(\\bW_{\\scriptsize\\mbox{uu}} \\bhu_k + \\bW_{\\scriptsize\\mbox{um}}\\,\\bhm_k),$$\n",
    "\n",
    "where $\\bW_{\\scriptsize\\mbox{um}}\\,\\bhm_k$ stems for the discretization of $(\\hat{m}_k \\exp(m_k)\\nabla p_k, \\nabla \\tilde{u})$.\n",
    "\n",
    "- Define the Hessian action\n",
    "\n",
    "$$\n",
    "  \\bH_k \\, \\bhm = \\underbrace{(\\bR + \\bW_{\\scriptsize\\mbox{mm}})}_{\\text{Hessian of the regularization}} \\bhm +\n",
    "    \\underbrace{(\\bC_k^{T}\\bA_k^{-T} (\\bW_{\\scriptsize\\mbox{uu}}\n",
    "    \\bA_k^{-1} \\bC_k - \\bW_{\\scriptsize\\mbox{um}}) -\n",
    "    \\bW_{\\scriptsize\\mbox{mu}} \\bA_k^{-1}\n",
    "    \\bC_k)}_{\\text{Hessian of the data misfit}}\\;\\bhm.\n",
    "$$\n",
    "\n",
    "### Goals:\n",
    "\n",
    "By the end of this notebook, you should be able to:\n",
    "\n",
    "- solve the forward and adjoint Poisson equations\n",
    "- understand the inverse method framework\n",
    "- visualise and understand the results\n",
    "- modify the problem and code\n",
    "\n",
    "### Mathematical tools used:\n",
    "\n",
    "- Finite element method\n",
    "- Derivation of gradiant and Hessian via the adjoint method\n",
    "- inexact Newton-CG\n",
    "- Armijo line search\n",
    "\n",
    "### List of software used:\n",
    "\n",
    "- <a href=\"http://fenicsproject.org/\">FEniCS</a>, a parallel finite element element library for the discretization of partial differential equations\n",
    "- <a href=\"http://www.mcs.anl.gov/petsc/\">PETSc</a>, for scalable and efficient linear algebra operations and solvers\n",
    "- <a href=\"http://matplotlib.org/\">Matplotlib</a>, a python package used for plotting the results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up\n",
    "\n",
    "### Import dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from dolfin import *\n",
    "import numpy as np\n",
    "\n",
    "import sys\n",
    "sys.path.append( \"../hippylib\" )\n",
    "from hippylib import *\n",
    "\n",
    "import logging\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "sys.path.append( \"../hippylib/tutorial\" )\n",
    "import nb\n",
    "\n",
    "logging.getLogger('FFC').setLevel(logging.WARNING)\n",
    "logging.getLogger('UFL').setLevel(logging.WARNING)\n",
    "set_log_active(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model set up:\n",
    "\n",
    "As in the introduction, the first thing we need to do is set up the numerical model.  In this cell, we set the mesh, the finite element functions $u, p, g$ corresponding to state, adjoint and coefficient/gradient variables, and the corresponding test functions and the parameters for the optimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# create mesh and define function spaces\n",
    "nx = 64\n",
    "ny = 64\n",
    "mesh = UnitSquareMesh(nx, ny)\n",
    "Vm = FunctionSpace(mesh, 'Lagrange', 1)\n",
    "Vu = FunctionSpace(mesh, 'Lagrange', 2)\n",
    "\n",
    "# The true and inverted parameter\n",
    "mtrue = interpolate(Expression('log(2 + 7*(pow(pow(x[0] - 0.5,2) + pow(x[1] - 0.5,2),0.5) > 0.2))', degree=5),Vm)\n",
    "m = interpolate(Expression(\"log(2.0)\", degree=1),Vm)\n",
    "\n",
    "# define function for state and adjoint\n",
    "u = Function(Vu)\n",
    "p = Function(Vu)\n",
    "\n",
    "# define Trial and Test Functions\n",
    "u_trial, p_trial, m_trial = TrialFunction(Vu), TrialFunction(Vu), TrialFunction(Vm)\n",
    "u_test, p_test, m_test = TestFunction(Vu), TestFunction(Vu), TestFunction(Vm)\n",
    "\n",
    "# initialize input functions\n",
    "f = Constant(\"1.0\")\n",
    "u0 = Constant(\"0.0\")\n",
    "\n",
    "# plot\n",
    "plt.figure(figsize=(15,5))\n",
    "nb.plot(mesh,subplot_loc=121, mytitle=\"Mesh\", show_axis='on')\n",
    "nb.plot(mtrue,subplot_loc=122, mytitle=\"True parameter field\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# set up dirichlet boundary conditions\n",
    "def boundary(x,on_boundary):\n",
    "    return on_boundary\n",
    "\n",
    "bc_state = DirichletBC(Vu, u0, boundary)\n",
    "bc_adj = DirichletBC(Vu, Constant(0.), boundary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up synthetic observations:\n",
    "\n",
    "- Propose a coefficient field $m_{\\rm true}$ shown above\n",
    "- The weak form of the pde: \n",
    "    Find $u\\in H_0^1(\\Omega)$ such that $\\underbrace{(\\exp(m_{\\rm true})\\nabla u,\\nabla v)}_{\\; := \\; a_{pde}} - \\underbrace{(f,v)}_{\\; := \\;L_{pde}} = 0, \\text{ for all } v\\in H_0^1(\\Omega)$.\n",
    "\n",
    "- Perturb the solution: $u = u + \\eta$, where $\\eta \\sim \\mathcal{N}(0, \\sigma)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# noise level\n",
    "noise_level = 0.05\n",
    "\n",
    "# weak form for setting up the synthetic observations\n",
    "a_goal = inner(exp(mtrue) * nabla_grad(u_trial), nabla_grad(u_test)) * dx\n",
    "L_goal = f * u_test * dx\n",
    "\n",
    "# solve the forward/state problem to generate synthetic observations\n",
    "goal_A, goal_b = assemble_system(a_goal, L_goal, bc_state)\n",
    "\n",
    "utrue = Function(Vu)\n",
    "solve(goal_A, utrue.vector(), goal_b)\n",
    "\n",
    "ud = Function(Vu)\n",
    "ud.assign(utrue)\n",
    "\n",
    "# perturb state solution and create synthetic measurements ud\n",
    "# ud = u + ||u||/SNR * random.normal\n",
    "MAX = ud.vector().norm(\"linf\")\n",
    "noise = Vector()\n",
    "goal_A.init_vector(noise,1)\n",
    "noise.set_local( noise_level * MAX * np.random.normal(0, 1, len(ud.vector().array())) )\n",
    "bc_adj.apply(noise)\n",
    "\n",
    "ud.vector().axpy(1., noise)\n",
    "\n",
    "# plot\n",
    "nb.multi1_plot([utrue, ud], [\"State solution with mtrue\", \"Synthetic observations\"])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The cost function evaluation:\n",
    "\n",
    "$$\n",
    "J(m):=\\underbrace{\\frac{1}{2}\\int_\\Omega (u-u_d)^2\\, dx}_{\\text{misfit} } + \\underbrace{\\frac{\\gamma}{2}\\int_\\Omega|\\nabla m|^2\\,dx}_{\\text{reg}}\n",
    "$$\n",
    "\n",
    "In the code below, $\\bW$ and $\\bR$ are symmetric positive definite matrices that stem from finite element discretization of the misfit and regularization component of the cost functional, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# regularization parameter\n",
    "gamma = 1e-8\n",
    "\n",
    "# weak for for setting up the misfit and regularization compoment of the cost\n",
    "W_equ   = inner(u_trial, u_test) * dx\n",
    "R_equ   = gamma * inner(nabla_grad(m_trial), nabla_grad(m_test)) * dx\n",
    "\n",
    "W = assemble(W_equ)\n",
    "R = assemble(R_equ)\n",
    "\n",
    "# refine cost function\n",
    "def cost(u, ud, m, W, R):\n",
    "    diff = u.vector() - ud.vector()\n",
    "    reg = 0.5 * m.vector().inner(R*m.vector() ) \n",
    "    misfit = 0.5 * diff.inner(W * diff)\n",
    "    return [reg + misfit, misfit, reg]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting up the state equations, right hand side for the adjoint and the necessary matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# weak form for setting up the state equation\n",
    "a_state = inner(exp(m) * nabla_grad(u_trial), nabla_grad(u_test)) * dx\n",
    "L_state = f * u_test * dx\n",
    "\n",
    "# weak form for setting up the adjoint equation\n",
    "a_adj = inner(exp(m) * nabla_grad(p_trial), nabla_grad(p_test)) * dx\n",
    "L_adj = -inner(u - ud, p_test) * dx\n",
    "\n",
    "# weak form for setting up matrices\n",
    "Wum_equ = inner(exp(m) * m_trial * nabla_grad(p_test), nabla_grad(p)) * dx\n",
    "C_equ   = inner(exp(m) * m_trial * nabla_grad(u), nabla_grad(u_test)) * dx\n",
    "Wmm_equ = inner(exp(m) * m_trial * m_test *  nabla_grad(u),  nabla_grad(p)) * dx\n",
    "\n",
    "M_equ   = inner(m_trial, m_test) * dx\n",
    "\n",
    "# assemble matrix M\n",
    "M = assemble(M_equ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initial guess\n",
    "We solve the state equation and compute the cost functional for the initial guess of the parameter ``a_ini``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# solve state equation\n",
    "state_A, state_b = assemble_system (a_state, L_state, bc_state)\n",
    "solve (state_A, u.vector(), state_b)\n",
    "\n",
    "# evaluate cost\n",
    "[cost_old, misfit_old, reg_old] = cost(u, ud, m, W, R)\n",
    "\n",
    "# plot\n",
    "plt.figure(figsize=(15,5))\n",
    "nb.plot(m,subplot_loc=121, mytitle=\"m_ini\", vmin=mtrue.vector().min(), vmax=mtrue.vector().max())\n",
    "nb.plot(u,subplot_loc=122, mytitle=\"u(m_ini)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The reduced Hessian apply to a vector $\\bhm$:\n",
    "\n",
    "Here we describe how to apply the reduced Hessian operator to a vector $\\bhm$. For an opportune choice of the regularization, the reduced Hessian operator evaluated in a neighborhood of the solution is positive define, whereas far from the solution the reduced Hessian may be indefinite. On the constrary, the Gauss-Newton approximation of the Hessian is always positive defined.\n",
    "\n",
    "For this reason, it is beneficial to perform a few initial Gauss-Newton steps (5 in this particular example) to accelerate the convergence of the inexact Newton-CG algorithm.\n",
    "\n",
    "The Hessian apply reads:\n",
    "$$\n",
    "\\begin{align}\n",
    "\\bhu &= -\\bA^{-1} \\bC \\bhm\\, & \\text{linearized forward}\\\\\n",
    "\\bhp &= -\\bA^{-T} (\\bW_{\\scriptsize\\mbox{uu}} \\bhu +\n",
    "\\bW_{\\scriptsize\\mbox{um}}\\,\\bhm) & \\text{adjoint}\\\\\n",
    "\\bH \\bhm &= (\\bR + \\bW_{\\scriptsize\\mbox{mm}})\\bhm + \\bC^T \\bhp + \\bW_{\\scriptsize\\mbox{mu}} \\bhu.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The Gauss-Newton Hessian apply is obtained by dropping the second derivatives operators $\\bW_{\\scriptsize\\mbox{um}}\\,\\bhm$, $\\bW_{\\scriptsize\\mbox{mm}}\\bf \\bhm$, and $\\bW_{\\scriptsize\\mbox{mu}} \\bhu$:\n",
    "$$\n",
    "\\begin{align}\n",
    "\\bhu &= -\\bA^{-1} \\bC \\bf \\bhm\\, & \\text{linearized forward}\\\\\n",
    "\\bhp &= -\\bA^{-T} \\bW_{\\scriptsize\\mbox{uu}} \\bhu & \\text{adjoint}\\\\\n",
    "\\bH_{\\rm GN} \\bhm &= \\bR \\bhm + \\bC^T \\bhp.\n",
    "\\end{align}\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Class HessianOperator to perform Hessian apply to a vector\n",
    "class HessianOperator():\n",
    "    cgiter = 0\n",
    "    def __init__(self, R, Wmm, C, A, adj_A, W, Wum, use_gaussnewton=False):\n",
    "        self.R = R\n",
    "        self.Wmm = Wmm\n",
    "        self.C = C\n",
    "        self.A = A\n",
    "        self.adj_A = adj_A\n",
    "        self.W = W\n",
    "        self.Wum = Wum\n",
    "        self.use_gaussnewton = use_gaussnewton\n",
    "        \n",
    "        # incremental state\n",
    "        self.du = Vector()\n",
    "        self.A.init_vector(self.du,0)\n",
    "        \n",
    "        #incremental adjoint\n",
    "        self.dp = Vector()\n",
    "        self.adj_A.init_vector(self.dp,0)\n",
    "        \n",
    "        # auxiliary vectors\n",
    "        self.CT_dp = Vector()\n",
    "        self.C.init_vector(self.CT_dp, 1)\n",
    "        self.Wum_du = Vector()\n",
    "        self.Wum.init_vector(self.Wum_du, 1)\n",
    "        \n",
    "    def init_vector(self, v, dim):\n",
    "        self.R.init_vector(v,dim)\n",
    "\n",
    "    # Hessian performed on v, output as generic vector y\n",
    "    def mult(self, v, y):\n",
    "        self.cgiter += 1\n",
    "        y.zero()\n",
    "        if self.use_gaussnewton:\n",
    "            self.mult_GaussNewton(v,y)\n",
    "        else:\n",
    "            self.mult_Newton(v,y)\n",
    "            \n",
    "    # define (Gauss-Newton) Hessian apply H * v\n",
    "    def mult_GaussNewton(self, v, y):\n",
    "        \n",
    "        #incremental forward\n",
    "        rhs = -(self.C * v)\n",
    "        bc_adj.apply(rhs)\n",
    "        solve (self.A, self.du, rhs)\n",
    "        \n",
    "        #incremental adjoint\n",
    "        rhs = - (self.W * self.du)\n",
    "        bc_adj.apply(rhs)\n",
    "        solve (self.adj_A, self.dp, rhs)\n",
    "        \n",
    "        # Reg/Prior term\n",
    "        self.R.mult(v,y)\n",
    "        \n",
    "        # Misfit term\n",
    "        self.C.transpmult(self.dp, self.CT_dp)\n",
    "        y.axpy(1, self.CT_dp)\n",
    "        \n",
    "    # define (Newton) Hessian apply H * v\n",
    "    def mult_Newton(self, v, y):\n",
    "        \n",
    "        #incremental forward\n",
    "        rhs = -(self.C * v)\n",
    "        bc_adj.apply(rhs)\n",
    "        solve (self.A, self.du, rhs)\n",
    "        \n",
    "        #incremental adjoint\n",
    "        rhs = -(self.W * self.du) -  self.Wum * v\n",
    "        bc_adj.apply(rhs)\n",
    "        solve (self.adj_A, self.dp, rhs)\n",
    "        \n",
    "        #Reg/Prior term\n",
    "        self.R.mult(v,y)\n",
    "        y.axpy(1., self.Wmm*v)\n",
    "        \n",
    "        #Misfit term\n",
    "        self.C.transpmult(self.dp, self.CT_dp)\n",
    "        y.axpy(1., self.CT_dp)\n",
    "        self.Wum.transpmult(self.du, self.Wum_du)\n",
    "        y.axpy(1., self.Wum_du)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The inexact Newton-CG optimization with Armijo line search:\n",
    "\n",
    "We solve the constrained optimization problem using the inexact Newton-CG method with Armijo line search.\n",
    "\n",
    "The stopping criterion is based on a relative reduction of the norm of the gradient (i.e. $\\frac{\\|g_{n}\\|}{\\|g_{0}\\|} \\leq \\tau$).\n",
    "\n",
    "First, we compute the gradient by solving the state and adjoint equation for the current parameter $m$, and then substituing the current state $u$, parameter $m$ and adjoint $p$ variables in the weak form expression of the gradient:\n",
    "$$ (g, \\tilde{m}) = \\gamma(\\nabla m, \\nabla \\tilde{m}) +(\\tilde{m}\\nabla u, \\nabla p).$$\n",
    "\n",
    "Then, we compute the Newton direction $\\hat m$ by iteratively solving $\\mathcal{H} {\\hat m} = -g$.\n",
    "The Newton system is solved inexactly by early termination of conjugate gradient iterations via Eisenstat–Walker (to prevent oversolving) and Steihaug  (to avoid negative curvature) criteria.\n",
    "\n",
    "Finally, the Armijo line search uses backtracking to find $\\alpha$ such that a sufficient reduction in the cost functional is achieved.\n",
    "More specifically, we use backtracking to find $\\alpha$ such that:\n",
    "$$J( m + \\alpha \\hat m ) \\leq J(m) + \\alpha c_{\\rm armijo} (\\hat m,g). $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# define parameters for the optimization\n",
    "tol = 1e-8\n",
    "c = 1e-4\n",
    "maxiter = 12\n",
    "plot_on = False\n",
    "\n",
    "# initialize iter counters\n",
    "iter = 1\n",
    "total_cg_iter = 0\n",
    "converged = False\n",
    "\n",
    "# initializations\n",
    "g, m_delta = Vector(), Vector()\n",
    "R.init_vector(m_delta,0)\n",
    "R.init_vector(g,0)\n",
    "\n",
    "m_prev = Function(Vm)\n",
    "\n",
    "print \"Nit   CGit   cost          misfit        reg           sqrt(-G*D)    ||grad||       alpha  tolcg\"\n",
    "\n",
    "while iter <  maxiter and not converged:\n",
    "\n",
    "    # assemble matrix C\n",
    "    C =  assemble(C_equ)\n",
    "\n",
    "    # solve the adoint problem\n",
    "    adjoint_A, adjoint_RHS = assemble_system(a_adj, L_adj, bc_adj)\n",
    "    solve(adjoint_A, p.vector(), adjoint_RHS)\n",
    "\n",
    "    # assemble W_ua and R\n",
    "    Wum = assemble (Wum_equ)\n",
    "    Wmm = assemble (Wmm_equ)\n",
    "\n",
    "    # evaluate the  gradient\n",
    "    CT_p = Vector()\n",
    "    C.init_vector(CT_p,1)\n",
    "    C.transpmult(p.vector(), CT_p)\n",
    "    MG = CT_p + R * m.vector()\n",
    "    solve(M, g, MG)\n",
    "\n",
    "    # calculate the norm of the gradient\n",
    "    grad2 = g.inner(MG)\n",
    "    gradnorm = sqrt(grad2)\n",
    "\n",
    "    # set the CG tolerance (use Eisenstat–Walker termination criterion)\n",
    "    if iter == 1:\n",
    "        gradnorm_ini = gradnorm\n",
    "    tolcg = min(0.5, sqrt(gradnorm/gradnorm_ini))\n",
    "\n",
    "    # define the Hessian apply operator (with preconditioner)\n",
    "    Hess_Apply = HessianOperator(R, Wmm, C, state_A, adjoint_A, W, Wum, use_gaussnewton=(iter<6) )\n",
    "    P = R + gamma * M\n",
    "    Psolver = PETScKrylovSolver(\"cg\", amg_method())\n",
    "    Psolver.set_operator(P)\n",
    "    \n",
    "    solver = CGSolverSteihaug()\n",
    "    solver.set_operator(Hess_Apply)\n",
    "    solver.set_preconditioner(Psolver)\n",
    "    solver.parameters[\"rel_tolerance\"] = tolcg\n",
    "    solver.parameters[\"zero_initial_guess\"] = True\n",
    "    solver.parameters[\"print_level\"] = -1\n",
    "\n",
    "    # solve the Newton system H a_delta = - MG\n",
    "    solver.solve(m_delta, -MG)\n",
    "    total_cg_iter += Hess_Apply.cgiter\n",
    "    \n",
    "    # linesearch\n",
    "    alpha = 1\n",
    "    descent = 0\n",
    "    no_backtrack = 0\n",
    "    m_prev.assign(m)\n",
    "    while descent == 0 and no_backtrack < 10:\n",
    "        m.vector().axpy(alpha, m_delta )\n",
    "\n",
    "        # solve the state/forward problem\n",
    "        state_A, state_b = assemble_system(a_state, L_state, bc_state)\n",
    "        solve(state_A, u.vector(), state_b)\n",
    "\n",
    "        # evaluate cost\n",
    "        [cost_new, misfit_new, reg_new] = cost(u, ud, m, W, R)\n",
    "\n",
    "        # check if Armijo conditions are satisfied\n",
    "        if cost_new < cost_old + alpha * c * MG.inner(m_delta):\n",
    "            cost_old = cost_new\n",
    "            descent = 1\n",
    "        else:\n",
    "            no_backtrack += 1\n",
    "            alpha *= 0.5\n",
    "            m.assign(m_prev)  # reset a\n",
    "\n",
    "    # calculate sqrt(-G * D)\n",
    "    graddir = sqrt(- MG.inner(m_delta) )\n",
    "\n",
    "    sp = \"\"\n",
    "    print \"%2d %2s %2d %3s %8.5e %1s %8.5e %1s %8.5e %1s %8.5e %1s %8.5e %1s %5.2f %1s %5.3e\" % \\\n",
    "        (iter, sp, Hess_Apply.cgiter, sp, cost_new, sp, misfit_new, sp, reg_new, sp, \\\n",
    "         graddir, sp, gradnorm, sp, alpha, sp, tolcg)\n",
    "\n",
    "    if plot_on:\n",
    "        nb.multi1_plot([m,u,p], [\"m\",\"u\",\"p\"], same_colorbar=False)\n",
    "        plt.show()\n",
    "    \n",
    "    # check for convergence\n",
    "    if gradnorm < tol and iter > 1:\n",
    "        converged = True\n",
    "        print \"Newton's method converged in \",iter,\"  iterations\"\n",
    "        print \"Total number of CG iterations: \", total_cg_iter\n",
    "        \n",
    "    iter += 1\n",
    "    \n",
    "if not converged:\n",
    "    print \"Newton's method did not converge in \", maxiter, \" iterations\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nb.multi1_plot([mtrue, m], [\"mtrue\", \"m\"])\n",
    "nb.multi1_plot([u,p], [\"u\",\"p\"], same_colorbar=False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Copyright (c) 2016, The University of Texas at Austin & University of California, Merced.\n",
    "All Rights reserved.\n",
    "See file COPYRIGHT for details.\n",
    "\n",
    "This file is part of the hIPPYlib library. For more information and source code\n",
    "availability see https://hippylib.github.io.\n",
    "\n",
    "hIPPYlib is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License (as published by the Free Software Foundation) version 2.0 dated June 1991."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}