tutorials/Simple quadratic problems.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple quadratic convex problems\n",
    "\n",
    "@Author: Ettore Biondi - ebiondi@caltech.edu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Instantiation of vectors and operator\n",
    "For testing the library we will be using a discretized version of the following operator:\n",
    "\\begin{align}\n",
    "y = \\frac{d^2f(x)}{dx^2},\n",
    "\\end{align}\n",
    "in which we simply compute the second-order derivative of a function $f(x)$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Importing necessary libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING! DATAPATH not found. The folder /tmp will be used to write binary files\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/nas/home/fpicetti/miniconda3/envs/occd/lib/python3.10/site-packages/dask_jobqueue/core.py:20: FutureWarning: tmpfile is deprecated and will be removed in a future release. Please use dask.utils.tmpfile instead.\n",
      "  from distributed.utils import tmpfile\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import occamypy as o\n",
    "\n",
    "# Plotting\n",
    "from matplotlib import rcParams\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "import matplotlib.pyplot as plt\n",
    "rcParams.update({\n",
    "    'image.cmap'     : 'gray',\n",
    "    'image.aspect'   : 'auto',\n",
    "    'image.interpolation': None,\n",
    "    'axes.grid'      : True,\n",
    "    'figure.figsize' : (6, 4),\n",
    "    'savefig.dpi'    : 300,\n",
    "    'axes.labelsize' : 14,\n",
    "    'axes.titlesize' : 16,\n",
    "    'font.size'      : 14,\n",
    "    'legend.fontsize': 14,\n",
    "    'xtick.labelsize': 14,\n",
    "    'ytick.labelsize': 14,\n",
    "    'text.usetex'    : True,\n",
    "    'font.family'    : 'serif',\n",
    "    'font.serif'     : 'Latin Modern Roman',\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 200  # Number of points of the function f(x)\n",
    "dx = 1/N  # Sampling of the function\n",
    "\n",
    "# The stencil used is simply: (f(ix-1)-2f(ix)+f(ix+1))/(dx*dx)\n",
    "D2 = np.zeros((N, N), dtype=np.float64)\n",
    "np.fill_diagonal(D2, -2 / (dx * dx))\n",
    "np.fill_diagonal(D2[1:], 1 / (dx * dx))\n",
    "np.fill_diagonal(D2[:,1:], 1 / (dx * dx))\n",
    "\n",
    "f = o.VectorNumpy((N, 1)).zero()\n",
    "y = f.clone()\n",
    "\n",
    "D2 = o.Matrix(o.VectorNumpy(D2), f, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we set any inversion problem, we study some of the properties of the constructed operator Deriv2Op."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dot-product tests of forward and adjoint operators\n",
      "--------------------------------------------------\n",
      "Applying forward operator add=False\n",
      " Runs in: 0.004538536071777344 seconds\n",
      "Applying adjoint operator add=False\n",
      " Runs in: 0.0005342960357666016 seconds\n",
      "Dot products add=False: domain=2.381674e+05 range=2.381674e+05 \n",
      "Absolute error: 5.820766e-11\n",
      "Relative error: 2.443982e-16 \n",
      "\n",
      "Applying forward operator add=True\n",
      " Runs in: 7.200241088867188e-05 seconds\n",
      "Applying adjoint operator add=True\n",
      " Runs in: 0.00022292137145996094 seconds\n",
      "Dot products add=True: domain=4.763347e+05 range=4.763347e+05 \n",
      "Absolute error: 1.164153e-10\n",
      "Relative error: 2.443982e-16 \n",
      "\n",
      "-------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "# Verifying operator adjointness through dot-product test\n",
    "D2.dotTest(verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computing maximum and minimum eigenvalues of the operator using the power iteration method and compare them against the ones computed using numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Maximum eigenvalue: -159990.22855475376 (Power method), -159990.2285552523 (NUMPY)\n",
      "Minimum eigenvalue: -9.771445230813697 (Power method), -9.77144474784679 (NUMPY)\n"
     ]
    }
   ],
   "source": [
    "egsOp = D2.powerMethod(verbose=False, eval_min=True, tol=1e-300)\n",
    "egsNp, _ = np.linalg.eig(D2.matrix.getNdArray())\n",
    "egsNp = egsNp[egsNp.argsort()[::-1]] #Sorting the eigenvalues\n",
    "print(\"\\nMaximum eigenvalue: %s (Power method), %s (NUMPY)\" % (egsOp[0], egsNp[-1]))\n",
    "print(  \"Minimum eigenvalue: %s (Power method), %s (NUMPY)\" % (egsOp[1], egsNp[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the matrix is negative definite. The small mismatch in the estimated eigenvalues is due to the dependence of the power method on the initial random eigenvector."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inversion tests\n",
    "We will now focus our attention on inverting a function knowning its second-order derivative. In this case we will assume that $y$ is constant and equal to $1$. Therefore, we expect to obtain a parabola with positive curvature. Given the chosen boundary conditions we know that the matrix is invertible since all eigenvalues have the same sign and are different then zero.\n",
    "We will solve the following objective functions using linear conjugate-gradient methods:\n",
    "\\begin{equation*}\n",
    "\\phi_1(\\mathbf{f}) = \\frac{1}{2}\\|D_2\\mathbf{f}-\\mathbf{y}\\|_2^2\n",
    "\\end{equation*}\n",
    "and\n",
    "\\begin{equation*}\n",
    "\\phi_2(\\mathbf{f}) = \\frac{1}{2}\\mathbf{f}^T D_2 \\mathbf{f} - \\mathbf{f}^{T} \\mathbf{y},\n",
    "\\end{equation*}\n",
    "where $D_2$ represents the discretized second-order derivative operator, while $\\mathbf{f}$ and $\\mathbf{y}$ are the discretized representations of $f$ and $y$, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "y.set(1.) # y = 1\n",
    "# Note that f = 0\n",
    "Phi1 = o.LeastSquares(f.clone().zero(), y, D2)\n",
    "Phi2 = o.LeastSquaresSymmetric(f.clone().zero(), y, D2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantiation of solver objects\n",
    "First, we create two different solver object for solving the two inversion problem stated above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create LCG solver\n",
    "LCGsolver = o.CG(o.BasicStopper(niter=2000))\n",
    "LCGsolver.setDefaults(save_obj=True)\n",
    "\n",
    "# Create LCG solver for symmetric systems\n",
    "SLCG = o.CGsym(o.BasicStopper(niter=2000))\n",
    "SLCG.setDefaults(save_obj=True, save_model=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Secondly, we run the solvers to minimize the objective functions previously defined."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "##########################################################################################\n",
      "            CG Solver\n",
      "    Restart folder: /tmp/restart_2022-04-22T01-50-53.913509/\n",
      "##########################################################################################\n",
      "\n",
      "iter = 0000, obj = 1.00000e+02, rnorm = 1.41e+01, gnorm = 5.66e+04, feval = 00002\n",
      "iter = 0001, obj = 9.98000e+01, rnorm = 1.41e+01, gnorm = 4.66e+04, feval = 00004\n",
      "iter = 0002, obj = 9.96483e+01, rnorm = 1.41e+01, gnorm = 4.20e+04, feval = 00006\n",
      "iter = 0003, obj = 9.95208e+01, rnorm = 1.41e+01, gnorm = 3.90e+04, feval = 00008\n",
      "iter = 0004, obj = 9.94086e+01, rnorm = 1.41e+01, gnorm = 3.69e+04, feval = 00010\n",
      "iter = 0005, obj = 9.93072e+01, rnorm = 1.41e+01, gnorm = 3.52e+04, feval = 00012\n",
      "iter = 0006, obj = 9.92139e+01, rnorm = 1.41e+01, gnorm = 3.39e+04, feval = 00014\n",
      "iter = 0007, obj = 9.91271e+01, rnorm = 1.41e+01, gnorm = 3.28e+04, feval = 00016\n",
      "iter = 0008, obj = 9.90455e+01, rnorm = 1.41e+01, gnorm = 3.18e+04, feval = 00018\n",
      "iter = 0009, obj = 9.89684e+01, rnorm = 1.41e+01, gnorm = 3.10e+04, feval = 00020\n",
      "iter = 0010, obj = 9.88950e+01, rnorm = 1.41e+01, gnorm = 3.03e+04, feval = 00022\n",
      "iter = 0011, obj = 9.88249e+01, rnorm = 1.41e+01, gnorm = 2.96e+04, feval = 00024\n",
      "iter = 0012, obj = 9.87576e+01, rnorm = 1.41e+01, gnorm = 2.91e+04, feval = 00026\n",
      "iter = 0013, obj = 9.86929e+01, rnorm = 1.40e+01, gnorm = 2.85e+04, feval = 00028\n",
      "iter = 0014, obj = 9.86304e+01, rnorm = 1.40e+01, gnorm = 2.80e+04, feval = 00030\n",
      "iter = 0015, obj = 9.85700e+01, rnorm = 1.40e+01, gnorm = 2.76e+04, feval = 00032\n",
      "iter = 0016, obj = 9.85114e+01, rnorm = 1.40e+01, gnorm = 2.72e+04, feval = 00034\n",
      "iter = 0017, obj = 9.84545e+01, rnorm = 1.40e+01, gnorm = 2.68e+04, feval = 00036\n",
      "iter = 0018, obj = 9.83992e+01, rnorm = 1.40e+01, gnorm = 2.64e+04, feval = 00038\n",
      "iter = 0019, obj = 9.83453e+01, rnorm = 1.40e+01, gnorm = 2.61e+04, feval = 00040\n",
      "iter = 0020, obj = 9.82927e+01, rnorm = 1.40e+01, gnorm = 2.58e+04, feval = 00042\n",
      "iter = 0021, obj = 9.82413e+01, rnorm = 1.40e+01, gnorm = 2.55e+04, feval = 00044\n",
      "iter = 0022, obj = 9.81911e+01, rnorm = 1.40e+01, gnorm = 2.52e+04, feval = 00046\n",
      "iter = 0023, obj = 9.81420e+01, rnorm = 1.40e+01, gnorm = 2.49e+04, feval = 00048\n",
      "iter = 0024, obj = 9.80939e+01, rnorm = 1.40e+01, gnorm = 2.47e+04, feval = 00050\n",
      "iter = 0025, obj = 9.80468e+01, rnorm = 1.40e+01, gnorm = 2.45e+04, feval = 00052\n",
      "iter = 0026, obj = 9.80005e+01, rnorm = 1.40e+01, gnorm = 2.42e+04, feval = 00054\n",
      "iter = 0027, obj = 9.79551e+01, rnorm = 1.40e+01, gnorm = 2.40e+04, feval = 00056\n",
      "iter = 0028, obj = 9.79104e+01, rnorm = 1.40e+01, gnorm = 2.38e+04, feval = 00058\n",
      "iter = 0029, obj = 9.78666e+01, rnorm = 1.40e+01, gnorm = 2.36e+04, feval = 00060\n",
      "iter = 0030, obj = 9.78234e+01, rnorm = 1.40e+01, gnorm = 2.34e+04, feval = 00062\n",
      "iter = 0031, obj = 9.77810e+01, rnorm = 1.40e+01, gnorm = 2.32e+04, feval = 00064\n",
      "iter = 0032, obj = 9.77392e+01, rnorm = 1.40e+01, gnorm = 2.30e+04, feval = 00066\n",
      "iter = 0033, obj = 9.76980e+01, rnorm = 1.40e+01, gnorm = 2.29e+04, feval = 00068\n",
      "iter = 0034, obj = 9.76575e+01, rnorm = 1.40e+01, gnorm = 2.27e+04, feval = 00070\n",
      "iter = 0035, obj = 9.76175e+01, rnorm = 1.40e+01, gnorm = 2.25e+04, feval = 00072\n",
      "iter = 0036, obj = 9.75780e+01, rnorm = 1.40e+01, gnorm = 2.24e+04, feval = 00074\n",
      "iter = 0037, obj = 9.75391e+01, rnorm = 1.40e+01, gnorm = 2.22e+04, feval = 00076\n",
      "iter = 0038, obj = 9.75007e+01, rnorm = 1.40e+01, gnorm = 2.21e+04, feval = 00078\n",
      "iter = 0039, obj = 9.74628e+01, rnorm = 1.40e+01, gnorm = 2.20e+04, feval = 00080\n",
      "iter = 0040, obj = 9.74253e+01, rnorm = 1.40e+01, gnorm = 2.18e+04, feval = 00082\n",
      "iter = 0041, obj = 9.73883e+01, rnorm = 1.40e+01, gnorm = 2.17e+04, feval = 00084\n",
      "iter = 0042, obj = 9.73517e+01, rnorm = 1.40e+01, gnorm = 2.16e+04, feval = 00086\n",
      "iter = 0043, obj = 9.73156e+01, rnorm = 1.40e+01, gnorm = 2.14e+04, feval = 00088\n",
      "iter = 0044, obj = 9.72799e+01, rnorm = 1.39e+01, gnorm = 2.13e+04, feval = 00090\n",
      "iter = 0045, obj = 9.72445e+01, rnorm = 1.39e+01, gnorm = 2.12e+04, feval = 00092\n",
      "iter = 0046, obj = 9.72096e+01, rnorm = 1.39e+01, gnorm = 2.11e+04, feval = 00094\n",
      "iter = 0047, obj = 9.71750e+01, rnorm = 1.39e+01, gnorm = 2.10e+04, feval = 00096\n",
      "iter = 0048, obj = 9.71407e+01, rnorm = 1.39e+01, gnorm = 2.09e+04, feval = 00098\n",
      "iter = 0049, obj = 9.71068e+01, rnorm = 1.39e+01, gnorm = 2.08e+04, feval = 00100\n",
      "iter = 0050, obj = 9.70733e+01, rnorm = 1.39e+01, gnorm = 2.07e+04, feval = 00102\n",
      "iter = 0051, obj = 9.70401e+01, rnorm = 1.39e+01, gnorm = 2.06e+04, feval = 00104\n",
      "iter = 0052, obj = 9.70071e+01, rnorm = 1.39e+01, gnorm = 2.05e+04, feval = 00106\n",
      "iter = 0053, obj = 9.69745e+01, rnorm = 1.39e+01, gnorm = 2.04e+04, feval = 00108\n",
      "iter = 0054, obj = 9.69422e+01, rnorm = 1.39e+01, gnorm = 2.03e+04, feval = 00110\n",
      "iter = 0055, obj = 9.69102e+01, rnorm = 1.39e+01, gnorm = 2.02e+04, feval = 00112\n",
      "iter = 0056, obj = 9.68785e+01, rnorm = 1.39e+01, gnorm = 2.01e+04, feval = 00114\n",
      "iter = 0057, obj = 9.68471e+01, rnorm = 1.39e+01, gnorm = 2.00e+04, feval = 00116\n",
      "iter = 0058, obj = 9.68159e+01, rnorm = 1.39e+01, gnorm = 1.99e+04, feval = 00118\n",
      "iter = 0059, obj = 9.67850e+01, rnorm = 1.39e+01, gnorm = 1.99e+04, feval = 00120\n",
      "iter = 0060, obj = 9.67543e+01, rnorm = 1.39e+01, gnorm = 1.98e+04, feval = 00122\n",
      "iter = 0061, obj = 9.67239e+01, rnorm = 1.39e+01, gnorm = 1.97e+04, feval = 00124\n",
      "iter = 0062, obj = 9.66937e+01, rnorm = 1.39e+01, gnorm = 1.96e+04, feval = 00126\n",
      "iter = 0063, obj = 9.66638e+01, rnorm = 1.39e+01, gnorm = 1.95e+04, feval = 00128\n",
      "iter = 0064, obj = 9.66341e+01, rnorm = 1.39e+01, gnorm = 1.95e+04, feval = 00130\n",
      "iter = 0065, obj = 9.66046e+01, rnorm = 1.39e+01, gnorm = 1.94e+04, feval = 00132\n",
      "iter = 0066, obj = 9.65753e+01, rnorm = 1.39e+01, gnorm = 1.93e+04, feval = 00134\n",
      "iter = 0067, obj = 9.65463e+01, rnorm = 1.39e+01, gnorm = 1.92e+04, feval = 00136\n",
      "iter = 0068, obj = 9.65175e+01, rnorm = 1.39e+01, gnorm = 1.92e+04, feval = 00138\n",
      "iter = 0069, obj = 9.64887e+01, rnorm = 1.39e+01, gnorm = 1.90e+04, feval = 00140\n",
      "iter = 0070, obj = 9.64595e+01, rnorm = 1.39e+01, gnorm = 1.86e+04, feval = 00142\n",
      "iter = 0071, obj = 9.64289e+01, rnorm = 1.39e+01, gnorm = 1.78e+04, feval = 00144\n",
      "iter = 0072, obj = 9.63958e+01, rnorm = 1.39e+01, gnorm = 1.69e+04, feval = 00146\n",
      "iter = 0073, obj = 9.63594e+01, rnorm = 1.39e+01, gnorm = 1.60e+04, feval = 00148\n",
      "iter = 0074, obj = 9.63193e+01, rnorm = 1.39e+01, gnorm = 1.51e+04, feval = 00150\n",
      "iter = 0075, obj = 9.62750e+01, rnorm = 1.39e+01, gnorm = 1.42e+04, feval = 00152\n",
      "iter = 0076, obj = 9.62261e+01, rnorm = 1.39e+01, gnorm = 1.33e+04, feval = 00154\n",
      "iter = 0077, obj = 9.61718e+01, rnorm = 1.39e+01, gnorm = 1.25e+04, feval = 00156\n",
      "iter = 0078, obj = 9.61117e+01, rnorm = 1.39e+01, gnorm = 1.17e+04, feval = 00158\n",
      "iter = 0079, obj = 9.60451e+01, rnorm = 1.39e+01, gnorm = 1.10e+04, feval = 00160\n",
      "iter = 0080, obj = 9.59710e+01, rnorm = 1.39e+01, gnorm = 1.02e+04, feval = 00162\n",
      "iter = 0081, obj = 9.58885e+01, rnorm = 1.38e+01, gnorm = 9.51e+03, feval = 00164\n",
      "iter = 0082, obj = 9.57994e+01, rnorm = 1.38e+01, gnorm = 1.48e+04, feval = 00166\n",
      "iter = 0083, obj = 9.57926e+01, rnorm = 1.38e+01, gnorm = 1.07e+04, feval = 00168\n",
      "iter = 0084, obj = 9.56911e+01, rnorm = 1.38e+01, gnorm = 9.04e+03, feval = 00170\n",
      "iter = 0085, obj = 9.56825e+01, rnorm = 1.38e+01, gnorm = 1.78e+04, feval = 00172\n",
      "iter = 0086, obj = 9.55757e+01, rnorm = 1.38e+01, gnorm = 1.31e+04, feval = 00174\n",
      "iter = 0087, obj = 9.55701e+01, rnorm = 1.38e+01, gnorm = 9.11e+03, feval = 00176\n",
      "iter = 0088, obj = 9.55005e+01, rnorm = 1.38e+01, gnorm = 3.82e+04, feval = 00178\n",
      "iter = 0089, obj = 9.54397e+01, rnorm = 1.38e+01, gnorm = 8.27e+03, feval = 00180\n",
      "iter = 0090, obj = 9.54344e+01, rnorm = 1.38e+01, gnorm = 1.24e+04, feval = 00182\n",
      "iter = 0091, obj = 9.53471e+01, rnorm = 1.38e+01, gnorm = 3.82e+04, feval = 00184\n",
      "iter = 0092, obj = 9.52903e+01, rnorm = 1.38e+01, gnorm = 9.52e+03, feval = 00186\n",
      "iter = 0093, obj = 9.52862e+01, rnorm = 1.38e+01, gnorm = 8.57e+03, feval = 00188\n",
      "iter = 0094, obj = 9.52539e+01, rnorm = 1.38e+01, gnorm = 3.14e+04, feval = 00190\n",
      "iter = 0095, obj = 9.51354e+01, rnorm = 1.38e+01, gnorm = 2.39e+04, feval = 00192\n",
      "iter = 0096, obj = 9.51176e+01, rnorm = 1.38e+01, gnorm = 7.11e+03, feval = 00194\n",
      "iter = 0097, obj = 9.51134e+01, rnorm = 1.38e+01, gnorm = 1.03e+04, feval = 00196\n",
      "iter = 0098, obj = 9.50593e+01, rnorm = 1.38e+01, gnorm = 3.70e+04, feval = 00198\n",
      "iter = 0099, obj = 9.49394e+01, rnorm = 1.38e+01, gnorm = 2.30e+04, feval = 00200\n",
      "iter = 0100, obj = 9.49226e+01, rnorm = 1.38e+01, gnorm = 7.51e+03, feval = 00202\n",
      "iter = 0101, obj = 9.49194e+01, rnorm = 1.38e+01, gnorm = 7.56e+03, feval = 00204\n",
      "iter = 0102, obj = 9.48992e+01, rnorm = 1.38e+01, gnorm = 2.46e+04, feval = 00206\n",
      "iter = 0103, obj = 9.47920e+01, rnorm = 1.38e+01, gnorm = 4.04e+04, feval = 00208\n",
      "iter = 0104, obj = 9.47042e+01, rnorm = 1.38e+01, gnorm = 1.53e+04, feval = 00210\n",
      "iter = 0105, obj = 9.46962e+01, rnorm = 1.38e+01, gnorm = 6.60e+03, feval = 00212\n",
      "iter = 0106, obj = 9.46933e+01, rnorm = 1.38e+01, gnorm = 7.43e+03, feval = 00214\n",
      "iter = 0107, obj = 9.46815e+01, rnorm = 1.38e+01, gnorm = 1.84e+04, feval = 00216\n",
      "iter = 0108, obj = 9.45974e+01, rnorm = 1.38e+01, gnorm = 4.02e+04, feval = 00218\n",
      "iter = 0109, obj = 9.44682e+01, rnorm = 1.37e+01, gnorm = 2.83e+04, feval = 00220\n",
      "iter = 0110, obj = 9.44356e+01, rnorm = 1.37e+01, gnorm = 9.94e+03, feval = 00222\n",
      "iter = 0111, obj = 9.44317e+01, rnorm = 1.37e+01, gnorm = 5.75e+03, feval = 00224\n",
      "iter = 0112, obj = 9.44286e+01, rnorm = 1.37e+01, gnorm = 7.69e+03, feval = 00226\n",
      "iter = 0113, obj = 9.44184e+01, rnorm = 1.37e+01, gnorm = 1.61e+04, feval = 00228\n",
      "iter = 0114, obj = 9.43520e+01, rnorm = 1.37e+01, gnorm = 3.91e+04, feval = 00230\n",
      "iter = 0115, obj = 9.41911e+01, rnorm = 1.37e+01, gnorm = 3.70e+04, feval = 00232\n",
      "iter = 0116, obj = 9.41331e+01, rnorm = 1.37e+01, gnorm = 1.62e+04, feval = 00234\n",
      "iter = 0117, obj = 9.41219e+01, rnorm = 1.37e+01, gnorm = 7.76e+03, feval = 00236\n",
      "iter = 0118, obj = 9.41190e+01, rnorm = 1.37e+01, gnorm = 5.04e+03, feval = 00238\n",
      "iter = 0119, obj = 9.41165e+01, rnorm = 1.37e+01, gnorm = 7.51e+03, feval = 00240\n",
      "iter = 0120, obj = 9.41064e+01, rnorm = 1.37e+01, gnorm = 1.61e+04, feval = 00242\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter = 0121, obj = 9.40543e+01, rnorm = 1.37e+01, gnorm = 3.20e+04, feval = 00244\n",
      "iter = 0122, obj = 9.39411e+01, rnorm = 1.37e+01, gnorm = 4.37e+04, feval = 00246\n",
      "iter = 0123, obj = 9.37967e+01, rnorm = 1.37e+01, gnorm = 3.13e+04, feval = 00248\n",
      "iter = 0124, obj = 9.37544e+01, rnorm = 1.37e+01, gnorm = 1.42e+04, feval = 00250\n",
      "iter = 0125, obj = 9.37464e+01, rnorm = 1.37e+01, gnorm = 7.09e+03, feval = 00252\n",
      "iter = 0126, obj = 9.37437e+01, rnorm = 1.37e+01, gnorm = 4.84e+03, feval = 00254\n",
      "iter = 0127, obj = 9.37419e+01, rnorm = 1.37e+01, gnorm = 5.55e+03, feval = 00256\n",
      "iter = 0128, obj = 9.37376e+01, rnorm = 1.37e+01, gnorm = 1.04e+04, feval = 00258\n",
      "iter = 0129, obj = 9.37140e+01, rnorm = 1.37e+01, gnorm = 2.07e+04, feval = 00260\n",
      "iter = 0130, obj = 9.36340e+01, rnorm = 1.37e+01, gnorm = 3.98e+04, feval = 00262\n",
      "iter = 0131, obj = 9.34869e+01, rnorm = 1.37e+01, gnorm = 4.41e+04, feval = 00264\n",
      "iter = 0132, obj = 9.33549e+01, rnorm = 1.37e+01, gnorm = 3.36e+04, feval = 00266\n",
      "iter = 0133, obj = 9.32981e+01, rnorm = 1.37e+01, gnorm = 1.74e+04, feval = 00268\n",
      "iter = 0134, obj = 9.32866e+01, rnorm = 1.37e+01, gnorm = 7.56e+03, feval = 00270\n",
      "iter = 0135, obj = 9.32833e+01, rnorm = 1.37e+01, gnorm = 5.67e+03, feval = 00272\n",
      "iter = 0136, obj = 9.32815e+01, rnorm = 1.37e+01, gnorm = 4.49e+03, feval = 00274\n",
      "iter = 0137, obj = 9.32798e+01, rnorm = 1.37e+01, gnorm = 5.42e+03, feval = 00276\n",
      "iter = 0138, obj = 9.32762e+01, rnorm = 1.37e+01, gnorm = 7.43e+03, feval = 00278\n",
      "iter = 0139, obj = 9.32671e+01, rnorm = 1.37e+01, gnorm = 1.45e+04, feval = 00280\n",
      "iter = 0140, obj = 9.32352e+01, rnorm = 1.37e+01, gnorm = 2.58e+04, feval = 00282\n",
      "iter = 0141, obj = 9.31132e+01, rnorm = 1.36e+01, gnorm = 4.32e+04, feval = 00284\n",
      "iter = 0142, obj = 9.28735e+01, rnorm = 1.36e+01, gnorm = 3.78e+04, feval = 00286\n",
      "iter = 0143, obj = 9.27669e+01, rnorm = 1.36e+01, gnorm = 2.95e+04, feval = 00288\n",
      "iter = 0144, obj = 9.27118e+01, rnorm = 1.36e+01, gnorm = 1.37e+04, feval = 00290\n",
      "iter = 0145, obj = 9.27027e+01, rnorm = 1.36e+01, gnorm = 7.51e+03, feval = 00292\n",
      "iter = 0146, obj = 9.26985e+01, rnorm = 1.36e+01, gnorm = 3.74e+03, feval = 00294\n",
      "iter = 0147, obj = 9.26961e+01, rnorm = 1.36e+01, gnorm = 3.29e+03, feval = 00296\n",
      "iter = 0148, obj = 9.26935e+01, rnorm = 1.36e+01, gnorm = 5.01e+03, feval = 00298\n",
      "iter = 0149, obj = 9.26860e+01, rnorm = 1.36e+01, gnorm = 9.65e+03, feval = 00300\n",
      "iter = 0150, obj = 9.26522e+01, rnorm = 1.36e+01, gnorm = 1.73e+04, feval = 00302\n",
      "iter = 0151, obj = 9.24926e+01, rnorm = 1.36e+01, gnorm = 3.42e+04, feval = 00304\n",
      "iter = 0152, obj = 9.22854e+01, rnorm = 1.36e+01, gnorm = 3.01e+04, feval = 00306\n",
      "iter = 0153, obj = 9.21362e+01, rnorm = 1.36e+01, gnorm = 3.20e+04, feval = 00308\n",
      "iter = 0154, obj = 9.19493e+01, rnorm = 1.36e+01, gnorm = 1.19e+04, feval = 00310\n",
      "iter = 0155, obj = 9.19361e+01, rnorm = 1.36e+01, gnorm = 1.23e+04, feval = 00312\n",
      "iter = 0156, obj = 9.19314e+01, rnorm = 1.36e+01, gnorm = 4.70e+03, feval = 00314\n",
      "iter = 0157, obj = 9.19266e+01, rnorm = 1.36e+01, gnorm = 2.37e+03, feval = 00316\n",
      "iter = 0158, obj = 9.19243e+01, rnorm = 1.36e+01, gnorm = 5.13e+03, feval = 00318\n",
      "iter = 0159, obj = 9.19233e+01, rnorm = 1.36e+01, gnorm = 3.65e+03, feval = 00320\n",
      "iter = 0160, obj = 9.19181e+01, rnorm = 1.36e+01, gnorm = 4.92e+03, feval = 00322\n",
      "iter = 0161, obj = 9.19036e+01, rnorm = 1.36e+01, gnorm = 2.24e+04, feval = 00324\n",
      "iter = 0162, obj = 9.18569e+01, rnorm = 1.36e+01, gnorm = 1.87e+04, feval = 00326\n",
      "iter = 0163, obj = 9.16903e+01, rnorm = 1.35e+01, gnorm = 5.44e+04, feval = 00328\n",
      "iter = 0164, obj = 9.15757e+01, rnorm = 1.35e+01, gnorm = 2.71e+04, feval = 00330\n",
      "iter = 0165, obj = 9.12069e+01, rnorm = 1.35e+01, gnorm = 7.56e+04, feval = 00332\n",
      "iter = 0166, obj = 9.10037e+01, rnorm = 1.35e+01, gnorm = 2.07e+04, feval = 00334\n",
      "iter = 0167, obj = 9.09239e+01, rnorm = 1.35e+01, gnorm = 3.74e+04, feval = 00336\n",
      "iter = 0168, obj = 9.08715e+01, rnorm = 1.35e+01, gnorm = 8.63e+03, feval = 00338\n",
      "iter = 0169, obj = 9.08660e+01, rnorm = 1.35e+01, gnorm = 1.04e+04, feval = 00340\n",
      "iter = 0170, obj = 9.08569e+01, rnorm = 1.35e+01, gnorm = 4.40e+03, feval = 00342\n",
      "iter = 0171, obj = 9.08562e+01, rnorm = 1.35e+01, gnorm = 3.63e+03, feval = 00344\n",
      "iter = 0172, obj = 9.08536e+01, rnorm = 1.35e+01, gnorm = 4.09e+03, feval = 00346\n",
      "iter = 0173, obj = 9.08531e+01, rnorm = 1.35e+01, gnorm = 1.89e+03, feval = 00348\n",
      "iter = 0174, obj = 9.08526e+01, rnorm = 1.35e+01, gnorm = 3.72e+03, feval = 00350\n",
      "iter = 0175, obj = 9.08506e+01, rnorm = 1.35e+01, gnorm = 3.40e+03, feval = 00352\n",
      "iter = 0176, obj = 9.08498e+01, rnorm = 1.35e+01, gnorm = 4.65e+03, feval = 00354\n",
      "iter = 0177, obj = 9.08413e+01, rnorm = 1.35e+01, gnorm = 1.20e+04, feval = 00356\n",
      "iter = 0178, obj = 9.08339e+01, rnorm = 1.35e+01, gnorm = 7.66e+03, feval = 00358\n",
      "iter = 0179, obj = 9.08272e+01, rnorm = 1.35e+01, gnorm = 1.37e+04, feval = 00360\n",
      "iter = 0180, obj = 9.07326e+01, rnorm = 1.35e+01, gnorm = 4.30e+04, feval = 00362\n",
      "iter = 0181, obj = 9.06355e+01, rnorm = 1.35e+01, gnorm = 2.36e+04, feval = 00364\n",
      "iter = 0182, obj = 9.05784e+01, rnorm = 1.35e+01, gnorm = 4.08e+04, feval = 00366\n",
      "iter = 0183, obj = 9.00352e+01, rnorm = 1.34e+01, gnorm = 8.76e+04, feval = 00368\n",
      "iter = 0184, obj = 8.97072e+01, rnorm = 1.34e+01, gnorm = 3.45e+04, feval = 00370\n",
      "iter = 0185, obj = 8.96426e+01, rnorm = 1.34e+01, gnorm = 3.18e+04, feval = 00372\n",
      "iter = 0186, obj = 8.94410e+01, rnorm = 1.34e+01, gnorm = 6.11e+04, feval = 00374\n",
      "iter = 0187, obj = 8.92568e+01, rnorm = 1.34e+01, gnorm = 1.88e+04, feval = 00376\n",
      "iter = 0188, obj = 8.92460e+01, rnorm = 1.34e+01, gnorm = 7.08e+03, feval = 00378\n",
      "iter = 0189, obj = 8.92411e+01, rnorm = 1.34e+01, gnorm = 1.07e+04, feval = 00380\n",
      "iter = 0190, obj = 8.92277e+01, rnorm = 1.34e+01, gnorm = 8.71e+03, feval = 00382\n",
      "iter = 0191, obj = 8.92245e+01, rnorm = 1.34e+01, gnorm = 3.72e+03, feval = 00384\n",
      "iter = 0192, obj = 8.92237e+01, rnorm = 1.34e+01, gnorm = 3.41e+03, feval = 00386\n",
      "iter = 0193, obj = 8.92222e+01, rnorm = 1.34e+01, gnorm = 4.01e+03, feval = 00388\n",
      "iter = 0194, obj = 8.92215e+01, rnorm = 1.34e+01, gnorm = 1.68e+03, feval = 00390\n",
      "iter = 0195, obj = 8.92214e+01, rnorm = 1.34e+01, gnorm = 1.43e+03, feval = 00392\n",
      "iter = 0196, obj = 8.92209e+01, rnorm = 1.34e+01, gnorm = 3.37e+03, feval = 00394\n",
      "iter = 0197, obj = 8.92193e+01, rnorm = 1.34e+01, gnorm = 4.12e+03, feval = 00396\n",
      "iter = 0198, obj = 8.92186e+01, rnorm = 1.34e+01, gnorm = 1.97e+03, feval = 00398\n",
      "iter = 0199, obj = 8.92182e+01, rnorm = 1.34e+01, gnorm = 3.12e+03, feval = 00400\n",
      "iter = 0200, obj = 8.92128e+01, rnorm = 1.34e+01, gnorm = 1.29e+04, feval = 00402\n",
      "iter = 0201, obj = 8.91858e+01, rnorm = 1.34e+01, gnorm = 1.81e+04, feval = 00404\n",
      "iter = 0202, obj = 8.91645e+01, rnorm = 1.34e+01, gnorm = 1.32e+04, feval = 00406\n",
      "iter = 0203, obj = 8.91540e+01, rnorm = 1.34e+01, gnorm = 1.09e+04, feval = 00408\n",
      "iter = 0204, obj = 8.91370e+01, rnorm = 1.34e+01, gnorm = 2.07e+04, feval = 00410\n",
      "iter = 0205, obj = 8.89646e+01, rnorm = 1.33e+01, gnorm = 6.17e+04, feval = 00412\n",
      "iter = 0206, obj = 8.84547e+01, rnorm = 1.33e+01, gnorm = 8.97e+04, feval = 00414\n",
      "iter = 0207, obj = 8.80282e+01, rnorm = 1.33e+01, gnorm = 4.73e+04, feval = 00416\n",
      "iter = 0208, obj = 8.79071e+01, rnorm = 1.33e+01, gnorm = 3.53e+04, feval = 00418\n",
      "iter = 0209, obj = 8.77829e+01, rnorm = 1.33e+01, gnorm = 5.18e+04, feval = 00420\n",
      "iter = 0210, obj = 8.73399e+01, rnorm = 1.32e+01, gnorm = 9.03e+04, feval = 00422\n",
      "iter = 0211, obj = 8.66031e+01, rnorm = 1.32e+01, gnorm = 6.73e+04, feval = 00424\n",
      "iter = 0212, obj = 8.63922e+01, rnorm = 1.31e+01, gnorm = 2.46e+04, feval = 00426\n",
      "iter = 0213, obj = 8.63691e+01, rnorm = 1.31e+01, gnorm = 9.95e+03, feval = 00428\n",
      "iter = 0214, obj = 8.63630e+01, rnorm = 1.31e+01, gnorm = 8.55e+03, feval = 00430\n",
      "iter = 0215, obj = 8.63474e+01, rnorm = 1.31e+01, gnorm = 1.51e+04, feval = 00432\n",
      "iter = 0216, obj = 8.63192e+01, rnorm = 1.31e+01, gnorm = 9.48e+03, feval = 00434\n",
      "iter = 0217, obj = 8.63154e+01, rnorm = 1.31e+01, gnorm = 2.92e+03, feval = 00436\n",
      "iter = 0218, obj = 8.63150e+01, rnorm = 1.31e+01, gnorm = 1.60e+03, feval = 00438\n",
      "iter = 0219, obj = 8.63148e+01, rnorm = 1.31e+01, gnorm = 1.65e+03, feval = 00440\n",
      "iter = 0220, obj = 8.63142e+01, rnorm = 1.31e+01, gnorm = 2.54e+03, feval = 00442\n",
      "iter = 0221, obj = 8.63130e+01, rnorm = 1.31e+01, gnorm = 2.45e+03, feval = 00444\n",
      "iter = 0222, obj = 8.63126e+01, rnorm = 1.31e+01, gnorm = 7.21e+02, feval = 00446\n",
      "iter = 0223, obj = 8.63126e+01, rnorm = 1.31e+01, gnorm = 8.32e+02, feval = 00448\n",
      "iter = 0224, obj = 8.63119e+01, rnorm = 1.31e+01, gnorm = 2.69e+03, feval = 00450\n",
      "iter = 0225, obj = 8.63099e+01, rnorm = 1.31e+01, gnorm = 3.75e+03, feval = 00452\n",
      "iter = 0226, obj = 8.63090e+01, rnorm = 1.31e+01, gnorm = 3.06e+03, feval = 00454\n",
      "iter = 0227, obj = 8.63087e+01, rnorm = 1.31e+01, gnorm = 1.72e+03, feval = 00456\n",
      "iter = 0228, obj = 8.63080e+01, rnorm = 1.31e+01, gnorm = 2.72e+03, feval = 00458\n",
      "iter = 0229, obj = 8.62963e+01, rnorm = 1.31e+01, gnorm = 1.15e+04, feval = 00460\n",
      "iter = 0230, obj = 8.62444e+01, rnorm = 1.31e+01, gnorm = 4.46e+04, feval = 00462\n",
      "iter = 0231, obj = 8.59920e+01, rnorm = 1.31e+01, gnorm = 4.85e+04, feval = 00464\n",
      "iter = 0232, obj = 8.56478e+01, rnorm = 1.31e+01, gnorm = 2.82e+04, feval = 00466\n",
      "iter = 0233, obj = 8.55580e+01, rnorm = 1.31e+01, gnorm = 1.86e+04, feval = 00468\n",
      "iter = 0234, obj = 8.53869e+01, rnorm = 1.31e+01, gnorm = 8.20e+04, feval = 00470\n",
      "iter = 0235, obj = 8.44416e+01, rnorm = 1.30e+01, gnorm = 8.75e+04, feval = 00472\n",
      "iter = 0236, obj = 8.26137e+01, rnorm = 1.29e+01, gnorm = 1.14e+05, feval = 00474\n",
      "iter = 0237, obj = 8.22005e+01, rnorm = 1.28e+01, gnorm = 6.06e+04, feval = 00476\n",
      "iter = 0238, obj = 8.18526e+01, rnorm = 1.28e+01, gnorm = 4.37e+04, feval = 00478\n",
      "iter = 0239, obj = 8.16727e+01, rnorm = 1.28e+01, gnorm = 5.17e+04, feval = 00480\n",
      "iter = 0240, obj = 8.15321e+01, rnorm = 1.28e+01, gnorm = 4.87e+04, feval = 00482\n",
      "iter = 0241, obj = 8.03312e+01, rnorm = 1.27e+01, gnorm = 9.98e+04, feval = 00484\n",
      "iter = 0242, obj = 7.97336e+01, rnorm = 1.26e+01, gnorm = 1.05e+05, feval = 00486\n",
      "iter = 0243, obj = 7.90122e+01, rnorm = 1.26e+01, gnorm = 2.74e+04, feval = 00488\n",
      "iter = 0244, obj = 7.89756e+01, rnorm = 1.26e+01, gnorm = 2.54e+04, feval = 00490\n",
      "iter = 0245, obj = 7.89336e+01, rnorm = 1.26e+01, gnorm = 6.04e+03, feval = 00492\n",
      "iter = 0246, obj = 7.89326e+01, rnorm = 1.26e+01, gnorm = 3.09e+03, feval = 00494\n",
      "iter = 0247, obj = 7.89316e+01, rnorm = 1.26e+01, gnorm = 2.29e+03, feval = 00496\n",
      "iter = 0248, obj = 7.89314e+01, rnorm = 1.26e+01, gnorm = 1.27e+03, feval = 00498\n",
      "iter = 0249, obj = 7.89308e+01, rnorm = 1.26e+01, gnorm = 4.13e+03, feval = 00500\n",
      "iter = 0250, obj = 7.89295e+01, rnorm = 1.26e+01, gnorm = 2.83e+03, feval = 00502\n",
      "iter = 0251, obj = 7.89281e+01, rnorm = 1.26e+01, gnorm = 4.32e+03, feval = 00504\n",
      "iter = 0252, obj = 7.89274e+01, rnorm = 1.26e+01, gnorm = 1.59e+03, feval = 00506\n",
      "iter = 0253, obj = 7.89272e+01, rnorm = 1.26e+01, gnorm = 1.86e+03, feval = 00508\n",
      "iter = 0254, obj = 7.89269e+01, rnorm = 1.26e+01, gnorm = 8.07e+02, feval = 00510\n",
      "iter = 0255, obj = 7.89269e+01, rnorm = 1.26e+01, gnorm = 4.22e+02, feval = 00512\n",
      "iter = 0256, obj = 7.89268e+01, rnorm = 1.26e+01, gnorm = 1.37e+03, feval = 00514\n",
      "iter = 0257, obj = 7.89267e+01, rnorm = 1.26e+01, gnorm = 1.03e+03, feval = 00516\n",
      "iter = 0258, obj = 7.89266e+01, rnorm = 1.26e+01, gnorm = 1.56e+03, feval = 00518\n",
      "iter = 0259, obj = 7.89264e+01, rnorm = 1.26e+01, gnorm = 6.97e+02, feval = 00520\n",
      "iter = 0260, obj = 7.89264e+01, rnorm = 1.26e+01, gnorm = 2.68e+02, feval = 00522\n",
      "iter = 0261, obj = 7.89264e+01, rnorm = 1.26e+01, gnorm = 4.18e+02, feval = 00524\n",
      "iter = 0262, obj = 7.89264e+01, rnorm = 1.26e+01, gnorm = 2.58e+02, feval = 00526\n",
      "iter = 0263, obj = 7.89264e+01, rnorm = 1.26e+01, gnorm = 3.54e+02, feval = 00528\n",
      "iter = 0264, obj = 7.89262e+01, rnorm = 1.26e+01, gnorm = 1.99e+03, feval = 00530\n",
      "iter = 0265, obj = 7.89257e+01, rnorm = 1.26e+01, gnorm = 1.94e+03, feval = 00532\n",
      "iter = 0266, obj = 7.89254e+01, rnorm = 1.26e+01, gnorm = 2.87e+03, feval = 00534\n",
      "iter = 0267, obj = 7.89233e+01, rnorm = 1.26e+01, gnorm = 5.85e+03, feval = 00536\n",
      "iter = 0268, obj = 7.89216e+01, rnorm = 1.26e+01, gnorm = 3.66e+03, feval = 00538\n",
      "iter = 0269, obj = 7.89199e+01, rnorm = 1.26e+01, gnorm = 5.88e+03, feval = 00540\n",
      "iter = 0270, obj = 7.89165e+01, rnorm = 1.26e+01, gnorm = 5.02e+03, feval = 00542\n",
      "iter = 0271, obj = 7.89156e+01, rnorm = 1.26e+01, gnorm = 1.96e+03, feval = 00544\n",
      "iter = 0272, obj = 7.89151e+01, rnorm = 1.26e+01, gnorm = 3.58e+03, feval = 00546\n",
      "iter = 0273, obj = 7.89135e+01, rnorm = 1.26e+01, gnorm = 3.53e+03, feval = 00548\n",
      "iter = 0274, obj = 7.89125e+01, rnorm = 1.26e+01, gnorm = 4.32e+03, feval = 00550\n",
      "iter = 0275, obj = 7.89084e+01, rnorm = 1.26e+01, gnorm = 1.02e+04, feval = 00552\n",
      "iter = 0276, obj = 7.88913e+01, rnorm = 1.26e+01, gnorm = 1.21e+04, feval = 00554\n",
      "iter = 0277, obj = 7.88823e+01, rnorm = 1.26e+01, gnorm = 1.04e+04, feval = 00556\n",
      "iter = 0278, obj = 7.88650e+01, rnorm = 1.26e+01, gnorm = 2.10e+04, feval = 00558\n",
      "iter = 0279, obj = 7.88125e+01, rnorm = 1.26e+01, gnorm = 2.25e+04, feval = 00560\n",
      "iter = 0280, obj = 7.87930e+01, rnorm = 1.26e+01, gnorm = 7.36e+03, feval = 00562\n",
      "iter = 0281, obj = 7.87882e+01, rnorm = 1.26e+01, gnorm = 1.05e+04, feval = 00564\n",
      "iter = 0282, obj = 7.87724e+01, rnorm = 1.26e+01, gnorm = 1.41e+04, feval = 00566\n",
      "iter = 0283, obj = 7.87614e+01, rnorm = 1.26e+01, gnorm = 7.84e+03, feval = 00568\n",
      "iter = 0284, obj = 7.87493e+01, rnorm = 1.25e+01, gnorm = 1.76e+04, feval = 00570\n",
      "iter = 0285, obj = 7.85940e+01, rnorm = 1.25e+01, gnorm = 4.43e+04, feval = 00572\n",
      "iter = 0286, obj = 7.84559e+01, rnorm = 1.25e+01, gnorm = 3.42e+04, feval = 00574\n",
      "iter = 0287, obj = 7.82204e+01, rnorm = 1.25e+01, gnorm = 7.84e+04, feval = 00576\n",
      "iter = 0288, obj = 7.63551e+01, rnorm = 1.24e+01, gnorm = 1.65e+05, feval = 00578\n",
      "iter = 0289, obj = 7.46982e+01, rnorm = 1.22e+01, gnorm = 8.06e+04, feval = 00580\n",
      "iter = 0290, obj = 7.38587e+01, rnorm = 1.22e+01, gnorm = 9.97e+04, feval = 00582\n",
      "iter = 0291, obj = 7.20386e+01, rnorm = 1.20e+01, gnorm = 9.76e+04, feval = 00584\n",
      "iter = 0292, obj = 7.14985e+01, rnorm = 1.20e+01, gnorm = 4.01e+04, feval = 00586\n",
      "iter = 0293, obj = 7.12858e+01, rnorm = 1.19e+01, gnorm = 6.56e+04, feval = 00588\n",
      "iter = 0294, obj = 7.05149e+01, rnorm = 1.19e+01, gnorm = 1.12e+05, feval = 00590\n",
      "iter = 0295, obj = 6.84456e+01, rnorm = 1.17e+01, gnorm = 1.05e+05, feval = 00592\n",
      "iter = 0296, obj = 6.73479e+01, rnorm = 1.16e+01, gnorm = 1.53e+05, feval = 00594\n",
      "iter = 0297, obj = 6.39557e+01, rnorm = 1.13e+01, gnorm = 3.13e+05, feval = 00596\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter = 0298, obj = 4.98865e+01, rnorm = 9.99e+00, gnorm = 4.18e+05, feval = 00598\n",
      "iter = 0299, obj = 3.35917e+01, rnorm = 8.20e+00, gnorm = 2.33e+05, feval = 00600\n",
      "iter = 0300, obj = 1.97988e+01, rnorm = 6.29e+00, gnorm = 3.08e+05, feval = 00602\n",
      "iter = 0301, obj = 9.29499e+00, rnorm = 4.31e+00, gnorm = 2.08e+05, feval = 00604\n",
      "iter = 0302, obj = 7.40993e+00, rnorm = 3.85e+00, gnorm = 1.18e+05, feval = 00606\n",
      "iter = 0303, obj = 6.99169e+00, rnorm = 3.74e+00, gnorm = 4.24e+04, feval = 00608\n",
      "iter = 0304, obj = 6.12529e+00, rnorm = 3.50e+00, gnorm = 6.53e+04, feval = 00610\n",
      "iter = 0305, obj = 5.84772e+00, rnorm = 3.42e+00, gnorm = 2.88e+04, feval = 00612\n",
      "iter = 0306, obj = 5.73638e+00, rnorm = 3.39e+00, gnorm = 5.93e+04, feval = 00614\n",
      "iter = 0307, obj = 5.32788e+00, rnorm = 3.26e+00, gnorm = 8.54e+04, feval = 00616\n",
      "iter = 0308, obj = 2.75342e-01, rnorm = 7.42e-01, gnorm = 5.77e+04, feval = 00618\n",
      "iter = 0309, obj = 1.46648e-01, rnorm = 5.42e-01, gnorm = 3.66e+04, feval = 00620\n",
      "iter = 0310, obj = 9.50978e-02, rnorm = 4.36e-01, gnorm = 1.23e+04, feval = 00622\n",
      "iter = 0311, obj = 7.03705e-02, rnorm = 3.75e-01, gnorm = 1.49e+04, feval = 00624\n",
      "iter = 0312, obj = 5.82868e-02, rnorm = 3.41e-01, gnorm = 1.59e+04, feval = 00626\n",
      "iter = 0313, obj = 3.13612e-02, rnorm = 2.50e-01, gnorm = 4.86e+03, feval = 00628\n",
      "iter = 0314, obj = 3.01569e-02, rnorm = 2.46e-01, gnorm = 3.67e+03, feval = 00630\n",
      "iter = 0315, obj = 2.96609e-02, rnorm = 2.44e-01, gnorm = 1.59e+03, feval = 00632\n",
      "iter = 0316, obj = 2.89121e-02, rnorm = 2.40e-01, gnorm = 3.36e+03, feval = 00634\n",
      "iter = 0317, obj = 2.80082e-02, rnorm = 2.37e-01, gnorm = 4.75e+03, feval = 00636\n",
      "iter = 0318, obj = 2.40827e-02, rnorm = 2.19e-01, gnorm = 3.89e+03, feval = 00638\n",
      "iter = 0319, obj = 2.29647e-02, rnorm = 2.14e-01, gnorm = 4.47e+03, feval = 00640\n",
      "iter = 0320, obj = 2.17708e-02, rnorm = 2.09e-01, gnorm = 3.47e+03, feval = 00642\n",
      "iter = 0321, obj = 1.47556e-02, rnorm = 1.72e-01, gnorm = 1.32e+04, feval = 00644\n",
      "iter = 0322, obj = 5.28578e-03, rnorm = 1.03e-01, gnorm = 6.27e+03, feval = 00646\n",
      "iter = 0323, obj = 1.93650e-03, rnorm = 6.22e-02, gnorm = 4.17e+03, feval = 00648\n",
      "iter = 0324, obj = 1.48571e-03, rnorm = 5.45e-02, gnorm = 1.05e+03, feval = 00650\n",
      "iter = 0325, obj = 1.32329e-03, rnorm = 5.14e-02, gnorm = 1.62e+03, feval = 00652\n",
      "iter = 0326, obj = 1.17840e-03, rnorm = 4.85e-02, gnorm = 1.34e+03, feval = 00654\n",
      "iter = 0327, obj = 9.29601e-04, rnorm = 4.31e-02, gnorm = 1.36e+03, feval = 00656\n",
      "iter = 0328, obj = 8.64856e-04, rnorm = 4.16e-02, gnorm = 5.27e+02, feval = 00658\n",
      "iter = 0329, obj = 8.50669e-04, rnorm = 4.12e-02, gnorm = 2.93e+02, feval = 00660\n",
      "iter = 0330, obj = 8.44711e-04, rnorm = 4.11e-02, gnorm = 3.59e+02, feval = 00662\n",
      "iter = 0331, obj = 8.04424e-04, rnorm = 4.01e-02, gnorm = 7.73e+02, feval = 00664\n",
      "iter = 0332, obj = 7.47683e-04, rnorm = 3.87e-02, gnorm = 1.12e+03, feval = 00666\n",
      "iter = 0333, obj = 4.56847e-04, rnorm = 3.02e-02, gnorm = 1.93e+03, feval = 00668\n",
      "iter = 0334, obj = 3.16368e-04, rnorm = 2.52e-02, gnorm = 8.38e+02, feval = 00670\n",
      "iter = 0335, obj = 2.49844e-04, rnorm = 2.24e-02, gnorm = 7.59e+02, feval = 00672\n",
      "iter = 0336, obj = 2.21775e-04, rnorm = 2.11e-02, gnorm = 5.80e+02, feval = 00674\n",
      "iter = 0337, obj = 1.64330e-04, rnorm = 1.81e-02, gnorm = 1.22e+03, feval = 00676\n",
      "iter = 0338, obj = 3.28861e-05, rnorm = 8.11e-03, gnorm = 5.60e+02, feval = 00678\n",
      "iter = 0339, obj = 2.14097e-05, rnorm = 6.54e-03, gnorm = 3.63e+02, feval = 00680\n",
      "iter = 0340, obj = 1.11507e-05, rnorm = 4.72e-03, gnorm = 1.29e+02, feval = 00682\n",
      "iter = 0341, obj = 1.06153e-05, rnorm = 4.61e-03, gnorm = 4.45e+01, feval = 00684\n",
      "iter = 0342, obj = 1.04430e-05, rnorm = 4.57e-03, gnorm = 5.59e+01, feval = 00686\n",
      "iter = 0343, obj = 1.01959e-05, rnorm = 4.52e-03, gnorm = 4.29e+01, feval = 00688\n",
      "iter = 0344, obj = 1.00086e-05, rnorm = 4.47e-03, gnorm = 6.74e+01, feval = 00690\n",
      "iter = 0345, obj = 9.19456e-06, rnorm = 4.29e-03, gnorm = 7.06e+01, feval = 00692\n",
      "iter = 0346, obj = 9.00833e-06, rnorm = 4.24e-03, gnorm = 2.88e+01, feval = 00694\n",
      "iter = 0347, obj = 8.87419e-06, rnorm = 4.21e-03, gnorm = 4.56e+01, feval = 00696\n",
      "iter = 0348, obj = 8.75230e-06, rnorm = 4.18e-03, gnorm = 3.58e+01, feval = 00698\n",
      "iter = 0349, obj = 8.08523e-06, rnorm = 4.02e-03, gnorm = 1.62e+02, feval = 00700\n",
      "iter = 0350, obj = 2.90767e-06, rnorm = 2.41e-03, gnorm = 1.49e+02, feval = 00702\n",
      "iter = 0351, obj = 2.13994e-06, rnorm = 2.07e-03, gnorm = 5.81e+01, feval = 00704\n",
      "iter = 0352, obj = 1.86481e-06, rnorm = 1.93e-03, gnorm = 6.56e+01, feval = 00706\n",
      "iter = 0353, obj = 1.38554e-06, rnorm = 1.66e-03, gnorm = 5.93e+01, feval = 00708\n",
      "iter = 0354, obj = 1.17868e-06, rnorm = 1.54e-03, gnorm = 5.68e+01, feval = 00710\n",
      "iter = 0355, obj = 5.12514e-07, rnorm = 1.01e-03, gnorm = 9.65e+01, feval = 00712\n",
      "iter = 0356, obj = 1.06203e-07, rnorm = 4.61e-04, gnorm = 2.62e+01, feval = 00714\n",
      "iter = 0357, obj = 5.53705e-08, rnorm = 3.33e-04, gnorm = 1.92e+01, feval = 00716\n",
      "iter = 0358, obj = 2.64433e-08, rnorm = 2.30e-04, gnorm = 5.93e+00, feval = 00718\n",
      "iter = 0359, obj = 2.44984e-08, rnorm = 2.21e-04, gnorm = 3.56e+00, feval = 00720\n",
      "iter = 0360, obj = 2.24334e-08, rnorm = 2.12e-04, gnorm = 6.35e+00, feval = 00722\n",
      "iter = 0361, obj = 1.93468e-08, rnorm = 1.97e-04, gnorm = 4.16e+00, feval = 00724\n",
      "iter = 0362, obj = 1.79202e-08, rnorm = 1.89e-04, gnorm = 4.39e+00, feval = 00726\n",
      "iter = 0363, obj = 1.43709e-08, rnorm = 1.70e-04, gnorm = 3.19e+00, feval = 00728\n",
      "iter = 0364, obj = 1.39977e-08, rnorm = 1.67e-04, gnorm = 1.12e+00, feval = 00730\n",
      "iter = 0365, obj = 1.39388e-08, rnorm = 1.67e-04, gnorm = 1.09e+00, feval = 00732\n",
      "iter = 0366, obj = 1.38076e-08, rnorm = 1.66e-04, gnorm = 8.38e-01, feval = 00734\n",
      "iter = 0367, obj = 1.37591e-08, rnorm = 1.66e-04, gnorm = 8.76e-01, feval = 00736\n",
      "iter = 0368, obj = 1.34702e-08, rnorm = 1.64e-04, gnorm = 2.65e+00, feval = 00738\n",
      "iter = 0369, obj = 1.07050e-08, rnorm = 1.46e-04, gnorm = 4.03e+00, feval = 00740\n",
      "iter = 0370, obj = 8.79042e-09, rnorm = 1.33e-04, gnorm = 5.87e+00, feval = 00742\n",
      "iter = 0371, obj = 3.70225e-09, rnorm = 8.60e-05, gnorm = 8.41e+00, feval = 00744\n",
      "iter = 0372, obj = 1.17343e-09, rnorm = 4.84e-05, gnorm = 2.29e+00, feval = 00746\n",
      "iter = 0373, obj = 1.02173e-09, rnorm = 4.52e-05, gnorm = 7.19e-01, feval = 00748\n",
      "iter = 0374, obj = 9.27027e-10, rnorm = 4.31e-05, gnorm = 1.01e+00, feval = 00750\n",
      "iter = 0375, obj = 7.81920e-10, rnorm = 3.95e-05, gnorm = 1.02e+00, feval = 00752\n",
      "iter = 0376, obj = 6.92414e-10, rnorm = 3.72e-05, gnorm = 1.21e+00, feval = 00754\n",
      "iter = 0377, obj = 5.22944e-10, rnorm = 3.23e-05, gnorm = 1.91e+00, feval = 00756\n",
      "iter = 0378, obj = 3.33336e-10, rnorm = 2.58e-05, gnorm = 5.44e-01, feval = 00758\n",
      "iter = 0379, obj = 3.16895e-10, rnorm = 2.52e-05, gnorm = 8.58e-02, feval = 00760\n",
      "iter = 0380, obj = 3.16093e-10, rnorm = 2.51e-05, gnorm = 9.12e-02, feval = 00762\n",
      "iter = 0381, obj = 3.14535e-10, rnorm = 2.51e-05, gnorm = 2.36e-01, feval = 00764\n",
      "iter = 0382, obj = 3.06968e-10, rnorm = 2.48e-05, gnorm = 2.04e-01, feval = 00766\n",
      "iter = 0383, obj = 3.01151e-10, rnorm = 2.45e-05, gnorm = 3.09e-01, feval = 00768\n",
      "iter = 0384, obj = 2.80177e-10, rnorm = 2.37e-05, gnorm = 8.89e-01, feval = 00770\n",
      "iter = 0385, obj = 9.24071e-11, rnorm = 1.36e-05, gnorm = 5.19e-01, feval = 00772\n",
      "iter = 0386, obj = 6.52657e-11, rnorm = 1.14e-05, gnorm = 7.56e-01, feval = 00774\n",
      "iter = 0387, obj = 3.22234e-11, rnorm = 8.03e-06, gnorm = 5.42e-01, feval = 00776\n",
      "iter = 0388, obj = 1.07378e-11, rnorm = 4.63e-06, gnorm = 2.72e-01, feval = 00778\n",
      "iter = 0389, obj = 6.21057e-12, rnorm = 3.52e-06, gnorm = 1.08e-01, feval = 00780\n",
      "iter = 0390, obj = 5.80519e-12, rnorm = 3.41e-06, gnorm = 8.17e-02, feval = 00782\n",
      "iter = 0391, obj = 3.70255e-12, rnorm = 2.72e-06, gnorm = 8.16e-02, feval = 00784\n",
      "iter = 0392, obj = 3.41827e-12, rnorm = 2.61e-06, gnorm = 7.24e-02, feval = 00786\n",
      "iter = 0393, obj = 2.61143e-12, rnorm = 2.29e-06, gnorm = 9.37e-02, feval = 00788\n",
      "iter = 0394, obj = 1.37797e-12, rnorm = 1.66e-06, gnorm = 1.23e-01, feval = 00790\n",
      "iter = 0395, obj = 8.87239e-13, rnorm = 1.33e-06, gnorm = 3.43e-02, feval = 00792\n",
      "iter = 0396, obj = 7.54208e-13, rnorm = 1.23e-06, gnorm = 4.13e-02, feval = 00794\n",
      "iter = 0397, obj = 6.58235e-13, rnorm = 1.15e-06, gnorm = 3.64e-02, feval = 00796\n",
      "iter = 0398, obj = 5.40312e-13, rnorm = 1.04e-06, gnorm = 1.08e-02, feval = 00798\n",
      "iter = 0399, obj = 5.33993e-13, rnorm = 1.03e-06, gnorm = 1.29e-02, feval = 00800\n",
      "iter = 0400, obj = 4.65828e-13, rnorm = 9.65e-07, gnorm = 1.89e-02, feval = 00802\n",
      "iter = 0401, obj = 4.50120e-13, rnorm = 9.49e-07, gnorm = 1.53e-02, feval = 00804\n",
      "iter = 0402, obj = 4.27915e-13, rnorm = 9.25e-07, gnorm = 1.39e-02, feval = 00806\n",
      "iter = 0403, obj = 4.05532e-13, rnorm = 9.01e-07, gnorm = 2.05e-02, feval = 00808\n",
      "iter = 0404, obj = 3.81808e-13, rnorm = 8.74e-07, gnorm = 1.10e-02, feval = 00810\n",
      "iter = 0405, obj = 3.75166e-13, rnorm = 8.66e-07, gnorm = 8.34e-03, feval = 00812\n",
      "iter = 0406, obj = 3.71874e-13, rnorm = 8.62e-07, gnorm = 3.05e-03, feval = 00814\n",
      "iter = 0407, obj = 3.71122e-13, rnorm = 8.62e-07, gnorm = 4.02e-03, feval = 00816\n",
      "iter = 0408, obj = 3.69605e-13, rnorm = 8.60e-07, gnorm = 3.26e-03, feval = 00818\n",
      "iter = 0409, obj = 3.68958e-13, rnorm = 8.59e-07, gnorm = 3.02e-03, feval = 00820\n",
      "iter = 0410, obj = 3.68370e-13, rnorm = 8.58e-07, gnorm = 1.79e-03, feval = 00822\n",
      "iter = 0411, obj = 3.67552e-13, rnorm = 8.57e-07, gnorm = 5.07e-03, feval = 00824\n",
      "iter = 0412, obj = 3.62704e-13, rnorm = 8.52e-07, gnorm = 8.26e-03, feval = 00826\n",
      "iter = 0413, obj = 3.53361e-13, rnorm = 8.41e-07, gnorm = 1.31e-02, feval = 00828\n",
      "iter = 0414, obj = 3.35062e-13, rnorm = 8.19e-07, gnorm = 8.69e-03, feval = 00830\n",
      "iter = 0415, obj = 3.29676e-13, rnorm = 8.12e-07, gnorm = 1.17e-02, feval = 00832\n",
      "iter = 0416, obj = 3.07496e-13, rnorm = 7.84e-07, gnorm = 1.22e-02, feval = 00834\n",
      "iter = 0417, obj = 2.79825e-13, rnorm = 7.48e-07, gnorm = 1.25e-02, feval = 00836\n",
      "iter = 0418, obj = 2.73747e-13, rnorm = 7.40e-07, gnorm = 6.29e-03, feval = 00838\n",
      "iter = 0419, obj = 2.68533e-13, rnorm = 7.33e-07, gnorm = 1.25e-02, feval = 00840\n",
      "iter = 0420, obj = 2.30644e-13, rnorm = 6.79e-07, gnorm = 2.49e-02, feval = 00842\n",
      "iter = 0421, obj = 1.78410e-13, rnorm = 5.97e-07, gnorm = 3.14e-02, feval = 00844\n",
      "iter = 0422, obj = 6.67242e-14, rnorm = 3.65e-07, gnorm = 2.46e-02, feval = 00846\n",
      "iter = 0423, obj = 4.36609e-14, rnorm = 2.96e-07, gnorm = 5.63e-03, feval = 00848\n",
      "iter = 0424, obj = 4.04943e-14, rnorm = 2.85e-07, gnorm = 6.23e-03, feval = 00850\n",
      "iter = 0425, obj = 3.68576e-14, rnorm = 2.72e-07, gnorm = 7.50e-03, feval = 00852\n",
      "iter = 0426, obj = 3.30861e-14, rnorm = 2.57e-07, gnorm = 3.44e-03, feval = 00854\n",
      "iter = 0427, obj = 3.25568e-14, rnorm = 2.55e-07, gnorm = 2.22e-03, feval = 00856\n",
      "iter = 0428, obj = 3.22894e-14, rnorm = 2.54e-07, gnorm = 7.05e-04, feval = 00858\n",
      "iter = 0429, obj = 3.22656e-14, rnorm = 2.54e-07, gnorm = 5.70e-04, feval = 00860\n",
      "iter = 0430, obj = 3.21875e-14, rnorm = 2.54e-07, gnorm = 8.79e-04, feval = 00862\n",
      "iter = 0431, obj = 3.19757e-14, rnorm = 2.53e-07, gnorm = 2.64e-03, feval = 00864\n",
      "iter = 0432, obj = 2.70359e-14, rnorm = 2.33e-07, gnorm = 6.31e-03, feval = 00866\n",
      "iter = 0433, obj = 2.53153e-14, rnorm = 2.25e-07, gnorm = 3.08e-03, feval = 00868\n",
      "iter = 0434, obj = 2.50119e-14, rnorm = 2.24e-07, gnorm = 1.63e-03, feval = 00870\n",
      "iter = 0435, obj = 2.45260e-14, rnorm = 2.21e-07, gnorm = 1.57e-03, feval = 00872\n",
      "iter = 0436, obj = 2.39502e-14, rnorm = 2.19e-07, gnorm = 3.77e-03, feval = 00874\n",
      "iter = 0437, obj = 1.99398e-14, rnorm = 2.00e-07, gnorm = 8.07e-03, feval = 00876\n",
      "iter = 0438, obj = 1.64989e-14, rnorm = 1.82e-07, gnorm = 3.45e-03, feval = 00878\n",
      "iter = 0439, obj = 1.50793e-14, rnorm = 1.74e-07, gnorm = 5.99e-03, feval = 00880\n",
      "iter = 0440, obj = 1.19652e-14, rnorm = 1.55e-07, gnorm = 6.92e-03, feval = 00882\n",
      "iter = 0441, obj = 7.75628e-15, rnorm = 1.25e-07, gnorm = 5.79e-03, feval = 00884\n",
      "iter = 0442, obj = 3.26849e-15, rnorm = 8.09e-08, gnorm = 4.98e-03, feval = 00886\n",
      "iter = 0443, obj = 1.96500e-15, rnorm = 6.27e-08, gnorm = 1.15e-03, feval = 00888\n",
      "iter = 0444, obj = 1.89977e-15, rnorm = 6.16e-08, gnorm = 7.11e-04, feval = 00890\n",
      "iter = 0445, obj = 1.79183e-15, rnorm = 5.99e-08, gnorm = 1.01e-03, feval = 00892\n",
      "iter = 0446, obj = 1.73697e-15, rnorm = 5.89e-08, gnorm = 1.02e-03, feval = 00894\n",
      "iter = 0447, obj = 1.49079e-15, rnorm = 5.46e-08, gnorm = 2.57e-03, feval = 00896\n",
      "iter = 0448, obj = 9.77349e-16, rnorm = 4.42e-08, gnorm = 1.03e-03, feval = 00898\n",
      "iter = 0449, obj = 9.26887e-16, rnorm = 4.31e-08, gnorm = 4.94e-04, feval = 00900\n",
      "iter = 0450, obj = 8.28870e-16, rnorm = 4.07e-08, gnorm = 1.03e-03, feval = 00902\n",
      "iter = 0451, obj = 7.57633e-16, rnorm = 3.89e-08, gnorm = 5.36e-04, feval = 00904\n",
      "iter = 0452, obj = 7.47690e-16, rnorm = 3.87e-08, gnorm = 3.51e-04, feval = 00906\n",
      "iter = 0453, obj = 7.36963e-16, rnorm = 3.84e-08, gnorm = 2.05e-04, feval = 00908\n",
      "iter = 0454, obj = 7.35612e-16, rnorm = 3.84e-08, gnorm = 5.11e-05, feval = 00910\n",
      "iter = 0455, obj = 7.35140e-16, rnorm = 3.83e-08, gnorm = 1.11e-04, feval = 00912\n",
      "iter = 0456, obj = 7.33973e-16, rnorm = 3.83e-08, gnorm = 1.09e-04, feval = 00914\n",
      "iter = 0457, obj = 7.33336e-16, rnorm = 3.83e-08, gnorm = 8.36e-05, feval = 00916\n",
      "iter = 0458, obj = 7.31843e-16, rnorm = 3.83e-08, gnorm = 9.49e-05, feval = 00918\n",
      "iter = 0459, obj = 7.31347e-16, rnorm = 3.82e-08, gnorm = 9.11e-05, feval = 00920\n",
      "iter = 0460, obj = 7.30758e-16, rnorm = 3.82e-08, gnorm = 9.93e-05, feval = 00922\n",
      "iter = 0461, obj = 6.92405e-16, rnorm = 3.72e-08, gnorm = 8.06e-04, feval = 00924\n",
      "iter = 0462, obj = 6.30754e-16, rnorm = 3.55e-08, gnorm = 1.16e-03, feval = 00926\n",
      "iter = 0463, obj = 5.39675e-16, rnorm = 3.29e-08, gnorm = 8.36e-04, feval = 00928\n",
      "iter = 0464, obj = 3.61770e-16, rnorm = 2.69e-08, gnorm = 9.75e-04, feval = 00930\n",
      "iter = 0465, obj = 3.15836e-16, rnorm = 2.51e-08, gnorm = 6.46e-04, feval = 00932\n",
      "iter = 0466, obj = 3.01237e-16, rnorm = 2.45e-08, gnorm = 2.42e-04, feval = 00934\n",
      "iter = 0467, obj = 2.92042e-16, rnorm = 2.42e-08, gnorm = 3.44e-04, feval = 00936\n",
      "iter = 0468, obj = 2.79823e-16, rnorm = 2.37e-08, gnorm = 5.05e-04, feval = 00938\n",
      "iter = 0469, obj = 2.62312e-16, rnorm = 2.29e-08, gnorm = 3.25e-04, feval = 00940\n",
      "iter = 0470, obj = 2.20013e-16, rnorm = 2.10e-08, gnorm = 6.94e-04, feval = 00942\n",
      "iter = 0471, obj = 2.04881e-16, rnorm = 2.02e-08, gnorm = 2.31e-04, feval = 00944\n",
      "iter = 0472, obj = 1.89051e-16, rnorm = 1.94e-08, gnorm = 7.43e-04, feval = 00946\n",
      "iter = 0473, obj = 1.18701e-16, rnorm = 1.54e-08, gnorm = 6.95e-04, feval = 00948\n",
      "iter = 0474, obj = 8.39335e-17, rnorm = 1.30e-08, gnorm = 5.76e-04, feval = 00950\n",
      "iter = 0475, obj = 4.70420e-17, rnorm = 9.70e-09, gnorm = 8.63e-04, feval = 00952\n",
      "iter = 0476, obj = 8.26902e-18, rnorm = 4.07e-09, gnorm = 1.39e-04, feval = 00954\n",
      "iter = 0477, obj = 7.59483e-18, rnorm = 3.90e-09, gnorm = 9.05e-05, feval = 00956\n",
      "iter = 0478, obj = 6.94317e-18, rnorm = 3.73e-09, gnorm = 5.33e-05, feval = 00958\n",
      "iter = 0479, obj = 6.62564e-18, rnorm = 3.64e-09, gnorm = 6.22e-05, feval = 00960\n",
      "iter = 0480, obj = 6.45186e-18, rnorm = 3.59e-09, gnorm = 3.86e-05, feval = 00962\n",
      "iter = 0481, obj = 6.19638e-18, rnorm = 3.52e-09, gnorm = 8.14e-05, feval = 00964\n",
      "iter = 0482, obj = 5.11900e-18, rnorm = 3.20e-09, gnorm = 1.24e-04, feval = 00966\n",
      "iter = 0483, obj = 3.21854e-18, rnorm = 2.54e-09, gnorm = 1.02e-04, feval = 00968\n",
      "iter = 0484, obj = 2.89830e-18, rnorm = 2.41e-09, gnorm = 2.92e-05, feval = 00970\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter = 0485, obj = 2.74799e-18, rnorm = 2.34e-09, gnorm = 6.39e-05, feval = 00972\n",
      "iter = 0486, obj = 2.24175e-18, rnorm = 2.12e-09, gnorm = 6.87e-05, feval = 00974\n",
      "iter = 0487, obj = 1.78427e-18, rnorm = 1.89e-09, gnorm = 4.58e-05, feval = 00976\n",
      "iter = 0488, obj = 1.72935e-18, rnorm = 1.86e-09, gnorm = 5.82e-06, feval = 00978\n",
      "iter = 0489, obj = 1.72792e-18, rnorm = 1.86e-09, gnorm = 3.80e-06, feval = 00980\n",
      "iter = 0490, obj = 1.72159e-18, rnorm = 1.86e-09, gnorm = 6.84e-06, feval = 00982\n",
      "iter = 0491, obj = 1.71999e-18, rnorm = 1.85e-09, gnorm = 4.21e-06, feval = 00984\n",
      "iter = 0492, obj = 1.71820e-18, rnorm = 1.85e-09, gnorm = 3.66e-06, feval = 00986\n",
      "iter = 0493, obj = 1.71538e-18, rnorm = 1.85e-09, gnorm = 3.85e-06, feval = 00988\n",
      "iter = 0494, obj = 1.71488e-18, rnorm = 1.85e-09, gnorm = 1.47e-06, feval = 00990\n",
      "iter = 0495, obj = 1.71463e-18, rnorm = 1.85e-09, gnorm = 1.97e-06, feval = 00992\n",
      "iter = 0496, obj = 1.71423e-18, rnorm = 1.85e-09, gnorm = 2.39e-06, feval = 00994\n",
      "iter = 0497, obj = 1.70901e-18, rnorm = 1.85e-09, gnorm = 1.18e-05, feval = 00996\n",
      "iter = 0498, obj = 1.68729e-18, rnorm = 1.84e-09, gnorm = 1.79e-05, feval = 00998\n",
      "iter = 0499, obj = 1.66690e-18, rnorm = 1.83e-09, gnorm = 1.72e-05, feval = 01000\n",
      "iter = 0500, obj = 1.64116e-18, rnorm = 1.81e-09, gnorm = 1.22e-05, feval = 01002\n",
      "iter = 0501, obj = 1.63390e-18, rnorm = 1.81e-09, gnorm = 1.26e-05, feval = 01004\n",
      "iter = 0502, obj = 1.59055e-18, rnorm = 1.78e-09, gnorm = 1.86e-05, feval = 01006\n",
      "iter = 0503, obj = 1.57066e-18, rnorm = 1.77e-09, gnorm = 1.47e-05, feval = 01008\n",
      "iter = 0504, obj = 1.53570e-18, rnorm = 1.75e-09, gnorm = 1.75e-05, feval = 01010\n",
      "iter = 0505, obj = 1.48959e-18, rnorm = 1.73e-09, gnorm = 3.45e-05, feval = 01012\n",
      "iter = 0506, obj = 1.31383e-18, rnorm = 1.62e-09, gnorm = 2.84e-05, feval = 01014\n",
      "iter = 0507, obj = 1.24966e-18, rnorm = 1.58e-09, gnorm = 2.85e-05, feval = 01016\n",
      "iter = 0508, obj = 1.22404e-18, rnorm = 1.56e-09, gnorm = 7.72e-06, feval = 01018\n",
      "iter = 0509, obj = 1.21043e-18, rnorm = 1.56e-09, gnorm = 1.32e-05, feval = 01020\n",
      "iter = 0510, obj = 1.18727e-18, rnorm = 1.54e-09, gnorm = 1.21e-05, feval = 01022\n",
      "iter = 0511, obj = 1.17578e-18, rnorm = 1.53e-09, gnorm = 1.38e-05, feval = 01024\n",
      "iter = 0512, obj = 1.14499e-18, rnorm = 1.51e-09, gnorm = 2.37e-05, feval = 01026\n",
      "iter = 0513, obj = 9.85113e-19, rnorm = 1.40e-09, gnorm = 4.52e-05, feval = 01028\n",
      "iter = 0514, obj = 6.99893e-19, rnorm = 1.18e-09, gnorm = 7.27e-05, feval = 01030\n",
      "iter = 0515, obj = 3.23181e-19, rnorm = 8.04e-10, gnorm = 3.45e-05, feval = 01032\n",
      "iter = 0516, obj = 2.29682e-19, rnorm = 6.78e-10, gnorm = 4.20e-05, feval = 01034\n",
      "iter = 0517, obj = 9.35155e-20, rnorm = 4.32e-10, gnorm = 4.19e-05, feval = 01036\n",
      "iter = 0518, obj = 1.71754e-20, rnorm = 1.85e-10, gnorm = 8.93e-06, feval = 01038\n",
      "iter = 0519, obj = 7.96273e-21, rnorm = 1.26e-10, gnorm = 7.70e-06, feval = 01040\n",
      "iter = 0520, obj = 3.77812e-21, rnorm = 8.69e-11, gnorm = 4.87e-06, feval = 01042\n",
      "iter = 0521, obj = 1.33310e-21, rnorm = 5.16e-11, gnorm = 5.07e-06, feval = 01044\n",
      "iter = 0522, obj = 4.24434e-22, rnorm = 2.91e-11, gnorm = 2.17e-06, feval = 01046\n",
      "iter = 0523, obj = 2.36707e-22, rnorm = 2.18e-11, gnorm = 8.52e-07, feval = 01048\n",
      "iter = 0524, obj = 1.91836e-22, rnorm = 1.96e-11, gnorm = 5.65e-07, feval = 01050\n",
      "iter = 0525, obj = 1.67103e-22, rnorm = 1.83e-11, gnorm = 4.95e-07, feval = 01052\n",
      "iter = 0526, obj = 1.47280e-22, rnorm = 1.72e-11, gnorm = 5.52e-07, feval = 01054\n",
      "iter = 0527, obj = 1.24152e-22, rnorm = 1.58e-11, gnorm = 4.04e-07, feval = 01056\n",
      "iter = 0528, obj = 1.18713e-22, rnorm = 1.54e-11, gnorm = 1.55e-07, feval = 01058\n",
      "iter = 0529, obj = 9.86537e-23, rnorm = 1.40e-11, gnorm = 6.48e-07, feval = 01060\n",
      "iter = 0530, obj = 7.97220e-23, rnorm = 1.26e-11, gnorm = 3.21e-07, feval = 01062\n",
      "iter = 0531, obj = 7.53422e-23, rnorm = 1.23e-11, gnorm = 1.42e-07, feval = 01064\n",
      "iter = 0532, obj = 7.10839e-23, rnorm = 1.19e-11, gnorm = 1.57e-07, feval = 01066\n",
      "iter = 0533, obj = 6.90866e-23, rnorm = 1.18e-11, gnorm = 1.75e-07, feval = 01068\n",
      "iter = 0534, obj = 6.69278e-23, rnorm = 1.16e-11, gnorm = 2.01e-07, feval = 01070\n",
      "iter = 0535, obj = 6.34731e-23, rnorm = 1.13e-11, gnorm = 9.26e-08, feval = 01072\n",
      "iter = 0536, obj = 5.61061e-23, rnorm = 1.06e-11, gnorm = 5.28e-07, feval = 01074\n",
      "iter = 0537, obj = 1.72560e-23, rnorm = 5.87e-12, gnorm = 4.33e-07, feval = 01076\n",
      "iter = 0538, obj = 1.12986e-23, rnorm = 4.75e-12, gnorm = 1.52e-07, feval = 01078\n",
      "iter = 0539, obj = 8.88374e-24, rnorm = 4.22e-12, gnorm = 5.97e-08, feval = 01080\n",
      "iter = 0540, obj = 8.78987e-24, rnorm = 4.19e-12, gnorm = 1.95e-08, feval = 01082\n",
      "Objective function didn't reduce, will terminate solver:\n",
      "\tobj_new = 8.78987e-24\tobj_cur = 8.78987e-24\n",
      "##########################################################################################\n",
      "            CG Solver end\n",
      "##########################################################################################\n",
      "\n"
     ]
    }
   ],
   "source": [
    "LCGsolver.run(Phi1, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "##########################################################################################\n",
      "            Symmetric CG Solver\n",
      "    Restart folder: /tmp/restart_2022-04-22T01-50-53.913671/\n",
      "##########################################################################################\n",
      "\n",
      "iter = 0000, obj = 0.00000e+00, rnorm = 1.41e+01, feval = 00001\n",
      "iter = 0001, obj = 2.50000e-01, rnorm = 1.41e+02, feval = 00003\n",
      "iter = 0002, obj = 4.95025e-01, rnorm = 1.39e+02, feval = 00005\n",
      "iter = 0003, obj = 7.35125e-01, rnorm = 1.38e+02, feval = 00007\n",
      "iter = 0004, obj = 9.70350e-01, rnorm = 1.36e+02, feval = 00009\n",
      "iter = 0005, obj = 1.20075e+00, rnorm = 1.35e+02, feval = 00011\n",
      "iter = 0006, obj = 1.42638e+00, rnorm = 1.34e+02, feval = 00013\n",
      "iter = 0007, obj = 1.64728e+00, rnorm = 1.32e+02, feval = 00015\n",
      "iter = 0008, obj = 1.86350e+00, rnorm = 1.31e+02, feval = 00017\n",
      "iter = 0009, obj = 2.07510e+00, rnorm = 1.29e+02, feval = 00019\n",
      "iter = 0010, obj = 2.28213e+00, rnorm = 1.28e+02, feval = 00021\n",
      "iter = 0011, obj = 2.48463e+00, rnorm = 1.27e+02, feval = 00023\n",
      "iter = 0012, obj = 2.68265e+00, rnorm = 1.25e+02, feval = 00025\n",
      "iter = 0013, obj = 2.87625e+00, rnorm = 1.24e+02, feval = 00027\n",
      "iter = 0014, obj = 3.06548e+00, rnorm = 1.22e+02, feval = 00029\n",
      "iter = 0015, obj = 3.25038e+00, rnorm = 1.21e+02, feval = 00031\n",
      "iter = 0016, obj = 3.43100e+00, rnorm = 1.19e+02, feval = 00033\n",
      "iter = 0017, obj = 3.60740e+00, rnorm = 1.18e+02, feval = 00035\n",
      "iter = 0018, obj = 3.77963e+00, rnorm = 1.17e+02, feval = 00037\n",
      "iter = 0019, obj = 3.94773e+00, rnorm = 1.15e+02, feval = 00039\n",
      "iter = 0020, obj = 4.11175e+00, rnorm = 1.14e+02, feval = 00041\n",
      "iter = 0021, obj = 4.27175e+00, rnorm = 1.12e+02, feval = 00043\n",
      "iter = 0022, obj = 4.42778e+00, rnorm = 1.11e+02, feval = 00045\n",
      "iter = 0023, obj = 4.57988e+00, rnorm = 1.10e+02, feval = 00047\n",
      "iter = 0024, obj = 4.72810e+00, rnorm = 1.08e+02, feval = 00049\n",
      "iter = 0025, obj = 4.87250e+00, rnorm = 1.07e+02, feval = 00051\n",
      "iter = 0026, obj = 5.01313e+00, rnorm = 1.05e+02, feval = 00053\n",
      "iter = 0027, obj = 5.15003e+00, rnorm = 1.04e+02, feval = 00055\n",
      "iter = 0028, obj = 5.28325e+00, rnorm = 1.03e+02, feval = 00057\n",
      "iter = 0029, obj = 5.41285e+00, rnorm = 1.01e+02, feval = 00059\n",
      "iter = 0030, obj = 5.53887e+00, rnorm = 9.97e+01, feval = 00061\n",
      "iter = 0031, obj = 5.66138e+00, rnorm = 9.83e+01, feval = 00063\n",
      "iter = 0032, obj = 5.78040e+00, rnorm = 9.69e+01, feval = 00065\n",
      "iter = 0033, obj = 5.89600e+00, rnorm = 9.55e+01, feval = 00067\n",
      "iter = 0034, obj = 6.00822e+00, rnorm = 9.40e+01, feval = 00069\n",
      "iter = 0035, obj = 6.11713e+00, rnorm = 9.26e+01, feval = 00071\n",
      "iter = 0036, obj = 6.22275e+00, rnorm = 9.12e+01, feval = 00073\n",
      "iter = 0037, obj = 6.32515e+00, rnorm = 8.98e+01, feval = 00075\n",
      "iter = 0038, obj = 6.42438e+00, rnorm = 8.84e+01, feval = 00077\n",
      "iter = 0039, obj = 6.52048e+00, rnorm = 8.70e+01, feval = 00079\n",
      "iter = 0040, obj = 6.61350e+00, rnorm = 8.56e+01, feval = 00081\n",
      "iter = 0041, obj = 6.70350e+00, rnorm = 8.41e+01, feval = 00083\n",
      "iter = 0042, obj = 6.79053e+00, rnorm = 8.27e+01, feval = 00085\n",
      "iter = 0043, obj = 6.87463e+00, rnorm = 8.13e+01, feval = 00087\n",
      "iter = 0044, obj = 6.95585e+00, rnorm = 7.99e+01, feval = 00089\n",
      "iter = 0045, obj = 7.03425e+00, rnorm = 7.85e+01, feval = 00091\n",
      "iter = 0046, obj = 7.10988e+00, rnorm = 7.71e+01, feval = 00093\n",
      "iter = 0047, obj = 7.18278e+00, rnorm = 7.57e+01, feval = 00095\n",
      "iter = 0048, obj = 7.25300e+00, rnorm = 7.42e+01, feval = 00097\n",
      "iter = 0049, obj = 7.32060e+00, rnorm = 7.28e+01, feval = 00099\n",
      "iter = 0050, obj = 7.38563e+00, rnorm = 7.14e+01, feval = 00101\n",
      "iter = 0051, obj = 7.44813e+00, rnorm = 7.00e+01, feval = 00103\n",
      "iter = 0052, obj = 7.50815e+00, rnorm = 6.86e+01, feval = 00105\n",
      "iter = 0053, obj = 7.56575e+00, rnorm = 6.72e+01, feval = 00107\n",
      "iter = 0054, obj = 7.62098e+00, rnorm = 6.58e+01, feval = 00109\n",
      "iter = 0055, obj = 7.67388e+00, rnorm = 6.43e+01, feval = 00111\n",
      "iter = 0056, obj = 7.72450e+00, rnorm = 6.29e+01, feval = 00113\n",
      "iter = 0057, obj = 7.77290e+00, rnorm = 6.15e+01, feval = 00115\n",
      "iter = 0058, obj = 7.81913e+00, rnorm = 6.01e+01, feval = 00117\n",
      "iter = 0059, obj = 7.86323e+00, rnorm = 5.87e+01, feval = 00119\n",
      "iter = 0060, obj = 7.90525e+00, rnorm = 5.73e+01, feval = 00121\n",
      "iter = 0061, obj = 7.94525e+00, rnorm = 5.59e+01, feval = 00123\n",
      "iter = 0062, obj = 7.98328e+00, rnorm = 5.44e+01, feval = 00125\n",
      "iter = 0063, obj = 8.01938e+00, rnorm = 5.30e+01, feval = 00127\n",
      "iter = 0064, obj = 8.05360e+00, rnorm = 5.16e+01, feval = 00129\n",
      "iter = 0065, obj = 8.08600e+00, rnorm = 5.02e+01, feval = 00131\n",
      "iter = 0066, obj = 8.11663e+00, rnorm = 4.88e+01, feval = 00133\n",
      "iter = 0067, obj = 8.14553e+00, rnorm = 4.74e+01, feval = 00135\n",
      "iter = 0068, obj = 8.17275e+00, rnorm = 4.60e+01, feval = 00137\n",
      "iter = 0069, obj = 8.19835e+00, rnorm = 4.45e+01, feval = 00139\n",
      "iter = 0070, obj = 8.22238e+00, rnorm = 4.31e+01, feval = 00141\n",
      "iter = 0071, obj = 8.24488e+00, rnorm = 4.17e+01, feval = 00143\n",
      "iter = 0072, obj = 8.26590e+00, rnorm = 4.03e+01, feval = 00145\n",
      "iter = 0073, obj = 8.28550e+00, rnorm = 3.89e+01, feval = 00147\n",
      "iter = 0074, obj = 8.30373e+00, rnorm = 3.75e+01, feval = 00149\n",
      "iter = 0075, obj = 8.32063e+00, rnorm = 3.61e+01, feval = 00151\n",
      "iter = 0076, obj = 8.33625e+00, rnorm = 3.46e+01, feval = 00153\n",
      "iter = 0077, obj = 8.35065e+00, rnorm = 3.32e+01, feval = 00155\n",
      "iter = 0078, obj = 8.36388e+00, rnorm = 3.18e+01, feval = 00157\n",
      "iter = 0079, obj = 8.37598e+00, rnorm = 3.04e+01, feval = 00159\n",
      "iter = 0080, obj = 8.38700e+00, rnorm = 2.90e+01, feval = 00161\n",
      "iter = 0081, obj = 8.39700e+00, rnorm = 2.76e+01, feval = 00163\n",
      "iter = 0082, obj = 8.40603e+00, rnorm = 2.62e+01, feval = 00165\n",
      "iter = 0083, obj = 8.41413e+00, rnorm = 2.47e+01, feval = 00167\n",
      "iter = 0084, obj = 8.42135e+00, rnorm = 2.33e+01, feval = 00169\n",
      "iter = 0085, obj = 8.42775e+00, rnorm = 2.19e+01, feval = 00171\n",
      "iter = 0086, obj = 8.43338e+00, rnorm = 2.05e+01, feval = 00173\n",
      "iter = 0087, obj = 8.43828e+00, rnorm = 1.91e+01, feval = 00175\n",
      "iter = 0088, obj = 8.44250e+00, rnorm = 1.77e+01, feval = 00177\n",
      "iter = 0089, obj = 8.44610e+00, rnorm = 1.62e+01, feval = 00179\n",
      "iter = 0090, obj = 8.44913e+00, rnorm = 1.48e+01, feval = 00181\n",
      "iter = 0091, obj = 8.45163e+00, rnorm = 1.34e+01, feval = 00183\n",
      "iter = 0092, obj = 8.45365e+00, rnorm = 1.20e+01, feval = 00185\n",
      "iter = 0093, obj = 8.45525e+00, rnorm = 1.06e+01, feval = 00187\n",
      "iter = 0094, obj = 8.45648e+00, rnorm = 9.17e+00, feval = 00189\n",
      "iter = 0095, obj = 8.45738e+00, rnorm = 7.75e+00, feval = 00191\n",
      "iter = 0096, obj = 8.45800e+00, rnorm = 6.32e+00, feval = 00193\n",
      "iter = 0097, obj = 8.45840e+00, rnorm = 4.90e+00, feval = 00195\n",
      "iter = 0098, obj = 8.45863e+00, rnorm = 3.46e+00, feval = 00197\n",
      "iter = 0099, obj = 8.45873e+00, rnorm = 2.00e+00, feval = 00199\n",
      "iter = 0100, obj = 8.45875e+00, rnorm = 4.71e-16, feval = 00201\n",
      "Objective function variation not monotonic, will terminate solver:\n",
      "\tobj_old=8.45873e+00\tobj_cur=8.45875e+00\tobj_new=8.45875e+00\n",
      "##########################################################################################\n",
      "            Symmetric CG Solver end\n",
      "##########################################################################################\n",
      "\n"
     ]
    }
   ],
   "source": [
    "SLCG.run(Phi2, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEoCAYAAAC6v50/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAjTUlEQVR4nO3de3hc9X3n8fd3dL+PZeviG7bHBoGdAJYNgeAUOZEpbRKWZG3zpJvrk0aQNGm6zdYu3fZptptdVt60uw1Ju1Zobs0mxTZJA6GlQYCSQLIYS+EOxvZgwFcZZNmWZVu33/4xR2I0Htm6zpkz83k96NGcy5zz/Y3M+cz5nZs55xARkewW8rsAERHxn8JAREQUBiIiojAQEREUBiIigsJAZBQza/S7BhE/KAxERgv7XYCIHxQGIiKiMBAZZmZhoNvnMkR8YboCWbKdd5wgAqwCHvZetzjnuv2sSySVtGcgWc3MmoGIc64FeNg5twNoAbb7W5lIauX6XYCIX8ysCQh7QTDCOddtZpVmFnHORb1564FvAlHgXi80Epe3ntheRbc3KgrUO+e2TGQeET8oDCSbNQNL4oa7416H42d0znWYGcDm4YCIZ2Zbie1ZbEkYty9ueDuxIBlzHhG/KAwkK5lZBGJ7Ad5wI7DLex0GKpNs9CNjBEEzsC/J3sJ2oCtunqfGmOe8ZYqkmsJAslWX9zMsHHfA+E5gc/zM8WGRMD4MbHLOWeI051zreOcR8ZvCQLKSd1wgGn9cAEb69M87jgCsI3amUaLVwMU26I2J83jhEvYGW3XmkvhNZxNJNtsArPcCYJ13QBnn3O1J5h21QTezcNwGvTvZwuNubdGdOI+3R3Abse6opO8XSSVdZyBZzztTqDvZ8YC4eY4752bFDW9yzm3xuoAecc6tSpi/KX7vwszak8xzHFiiMJB0oG4ikdiB4fNOFYWRA82bgS5vD6KSWJdRF4x0N20YPkBMbA8gWTfT8Dz7iB0wDhM7M6l7+psjMnHaM5Csl/gtXiQb6ZiByOizikSykvYMREREewYiIqIwEBERFAYiIoLCQEREUBiIiAgKAxERQWEgIiIoDEREBIWBiIigMBARERQGIiKCwkBERFAYiIgIAX64TTgcdsuWLfO7jJQ4ffo0JSUlfpeREmpr5smWdkL6t7W9vf1N51xVsmmBDYOamhp27drldxkp0dbWRkNDg99lpITamnmypZ2Q/m01s9fGmqZuIhERURiIiIjCQEREUBiIiAgKAxERQWEgIiIE+NTS3gHHIy8dxQwMw/sPMyPkjTNvHAnDZvGvYzOMmuaNH34PccMjr89b5oWXw8g8w+8/fz1j1Xqqz3H8dN/oto6jXq+i5PMOzyAiQoDDoLPX8envZsd1BgA8+vCMLPZC4fV2wCYPlfgAjp/mTUoSkOcvZ6QGb5Vnz5yl+KnHzqspJ2SEzEZ+h+zt9Y4M8/bw8O/83BAVRXnkhozcnNh7c0NGKDT6d47FXicG93BbzhuXEKaj548NlBflUlteSE15IRVFeVSVFVCYlzPJv5TIzApsGMwtCfGTP7gBBzjnvN8ADudgyI0e73B4/40Mu5Hh2HzEj0+YFlv228sfazkkjnfD7xpjPUmWQ8I6X9mzh2XLll14PQnDI9W68a1n1Pikn9nY60n8jBJrGvnMxqiXkffBkSNHqK4Jj/78HQw5x8CQY2jIMeStb+Rv7E0f8l4PevMMOcfx3iH2dJ5icNAx6ByDQ7HlDCb8DAzFfWgzJGSwrLqUP3//cn7rsqQXgYr4JrBhUJADVy0M+11GSrT17afhhiV+l5ESsSs4V/qy7iEvEOJjYSR4R43zfuNGDY9alnN09/Zz5ORZOk+e5eTZAQ509fLAs4f54j/9hp9vWjsDLRCZvMCGgch0C4WSHUeZ/LGV4vxc5oWLRo1rXF7DLV9/gm1PvUF23FlLgkJnE4mk0JULwly9MMx9HQf9LkVkFIWBSIp94Mq5vHT4JEdPD/ldisgIhYFIir338moAXuoa9LkSkbcpDERSbMmcEmaX5LO3W3sGkj4UBiIpZmasWjSLZ44N8MrRU36XIwIoDER88UeNl3FuEH7nb3/J5h3P8tDzRxhMwbUOImNJizAws7CZNZrZejNrNrOw3zWJzKTl88ppfk8RH1o5n58+e4g7vt/O/c/oDCPxT1qEAbARCDvndsQNi2S0WYUhvrrhKp740/cC8OapPp8rkmw2qTAws3oza08yPmJmm7xv+ZvG+w3fOdcSFwQRIDqZukSCqKwwD4BT5wZ8rkSy2YSvQDazRqALqE8yeatzbp03XxRoBm6fwLLDAM651onWJRJUOSGjtCCXU2f7/S5FstiEw2B4Q33+XRstkjBf1Mw24oWBmTUlWVw0YcN/p3Nuw0RrEgm6ssJcTp3VnoH4ZzrvTVQPdCeONLOIcy7qnGu50JvNbD1wl/e6UXsHkk1iYaA9A/HPdIZBZZJxXUD4Ym/0up6agW5vj2PzGPM1AU0AVVVVtLW1TbLUYOnp6VFbM1B8W4fOneH1w70Z2fZs/ZsGTVrctdTbC1g6jvlagBaAuro619DQMMOVpYfYbZ0b/C4jJbK1rd95dSdv9fTR0LDG36JmQLb+TYNmOk8tTbYXUEmSriMRGa2sME/dROKr6QyDDpJ0FTnndJqoyEXEzibSAWTxz7SFQeJG3zu7aNt0LV8kk1UU5XHybD/9g7p5nfhjwmEwfEGZ97rZO/g7bMPwRWfAeufcuK8xEMlmVy+soH/Q8ZvXu/0uRbLUZK8zaAW2JJkWjRuvU0NFxun6pXPICRltuzu5dkmyE/NEZla63JtIJKtVFOVxw7I5/Pg3BxlQV5H4QGEgkiY++q5LOHziLNvbD/hdimQhhYFImli3vIZrl1TyVw+8yN5OPfRGUkthIJImzIy7P7KScwOD/OTpQ36XI1kmLa5AFpGYmvJCrl4Y5kcdByktyOWKueWsWTaHUMgu/maRKdCegUia+fj1i+kbHOKuf32Zj39rJ//hnid56PkjfpclGU57BiJp5taV87l15XxO9Pazo+MA9/wyyh3fb+enX1jDO+ZX+F2eZCjtGYikqYriPD69ZgkP/dFvUZyfww92vu53SZLBFAYiaa6iKI+Vl4R54eAJv0uRDKYwEAmAS6vL2NPZw9CQ87sUyVAKA5EAqKsto7dvkIPdZ/wuRTKUwkAkAC6rKQXglaO6GE1mhsJAJACWVZcBsFthIDNEYSASABVFecytKGTP0R6/S5EMpTAQCYi62jJ2vtrF2f5Bv0uRDKQwEAmI318T4WD3GW7b+mv2v3na73IkwygMRAJizaVz+Mbv1bP/rV4+8e2dPPziUfbqdFOZJrodhUiAvP/Kucwuzecz393FZ763C4idaXTfZ99NWWGez9VJkGnPQCRgrovM5sn//D5+/Ll386V1l/HK0R7dyE6mTGEgEkDF+bmsvGQWn3/vMhbNLmb7Lj0dTaZGYSASYGbGx69fzM79XTyx902/y5EAUxiIBNxHrl3IkjklfO7/dtD+Wpff5UhAKQxEAq44P5dvf/IaSgty2Xzfczq7SCZFYSCSARbPKeH2GyPs7ezh8MmzfpcjAaQwEMkQCyuLAThyQmEgE6cwEMkQteWFgMJAJictwsDM6s1sn5m1ez/NftckEjRzK7wwUDeRTEK6XIFc6ZxbCrFgALr9LUckeCqK8ijIDXHkhB6AIxM3qT0D75t8e5LxETPbZGaN3u/weJbnnGuNG4w456KTqUskm5kZ88NF7O3Uba5l4iYcBmbW6L2sTzJ5q3Nui7dx3wFMqLvHzDY553ZMtCYRiVm3vIZf7HmTTnUVyQRNOAycc63OuY7E8WYWSZgvCmyMm96U5KcxYTHXTLQeEXnbh+sXMDjkeGx3p9+lSMBM5zGDpH39ZhZxzkWdcy0XerMXJuFprEck61xWU0p1WQH3P3OIutpyFs8uJlyc73dZEgDTGQaVScZ1MbENvI4ViEyBmXHzO2r53q9f44m9T2AGqy6ZRW1FITetqOWWq+b5XaKkKXNucpeum5lzzlnccBOwzjm3IW7cPmBDsm6lSa6zCWgCqKqqWrVt27bpWGza6+npobS01O8yUkJtnTrnHJ29joM9Q+zpHmJf9yCvnxzi7CDctaaIuaWpPaNcf9P0sXbt2nbn3Opk06ZzzyDZXkAl03iaqNfV1AJQV1fnGhoapmvRaa2trQ21NfOksq3PHzzBB+5+nNKFl9OQ4r0D/U2DYTq/InSQpKtIp4mK+G9ZdSkhgz067VTGMG1hkLjR9w4IZ0c/jkiaK8zLYdHsEl44eMLvUiRNTeo6AzPb5L1uTjg9dMPwRWfAeufc7dNVqIhMzU3La3h0d6dOO5WkJnudwRbnnDnnNsdfPeydQrpleJ7pLVVEpuJza5dRV1PGZ7/fzqFu3bJCRkuLG9WJyMyrKMrjnk/ETiT5zPd20T845HNFkk4UBiJZZMGsYr78wRW8cOgkLx466Xc5kkYUBiJZZs2lcwB49kC3v4VIWlEYiGSZ+eEiZpfk8+wBnVkkb1MYiGQZM2NpVSmvd/X6XYqkEYWBSBaqrSjUE9FkFIWBSBaaW1HI4RNnmey9ySTzKAxEslBtRSF9A0Mc7+33uxRJEwoDkSxUW14IwGE9L1k8CgORLHTJ7GIA9hzVjeskRmEgkoUury2nrDCXJ199y+9SJE0oDESyUE7IuHZxJU9Gu/wuRdKEwkAkS70rUkn0zdN06hRTQWEgkrWui8wGoOUXUc70DfpcjfhNYSCSpVbMq6DximruefxVbvn647S/pi6jbDadz0AWkQDJCRn3fOIa2nZ38ic7nuWj9+xk/aoFzCrOo6q8kA+vnE9JgTYR2UJ/aZEs11BXzYNfWMMf/KCDB549xIkz/TgHP999jJaPrSIUMr9LlBRQGIgI1eWFbL/j3QAMDTm+86v9/NVPX+QHO1/no9ct8rk6SQUdMxCRUUIh41M3LGbR7GIe3/Om3+VIiigMROQ8ZsY751fw3EE98yBbKAxEJKl3zK/gYPcZjp/u87sUSQGFgYgkdVlNKQD7jun+RdlAYSAiSUXmxMIg+uZpnyuRVFAYiEhSC2YVkZdjvKowyAoKAxFJKjcnxOLZJTz0/BFee0uBkOkUBiIypr/84Aq6Tvfxgbsf52cvHPG7HJlBvoSBmYXNrClhuNHM1ptZs5mF/ahLREZbc+kcfvqFNSyZU8Id32/n2QPdfpckM8SvPYMIcHvc8EYg7JzbETcsImlgYWUx3//9d1FZUsBf/+wVv8uRGXLR21GYWT3wTefcqoTxEWA90AHUAy3Oue7xrNQ512FmXXHDLXGTI8DD41mOiKRGeWEeH7l2IXc/upc7f/Qc1WUFzA8XcfM7aykvzPO7PJkGFwwDM2sEuoht7BNtdc6t8+aLAs2M/rY/YcPdQ8651qksR0Sm38euX8ST0S4eeOYQPecGAPjZi0e45xPX+FyZTIcLhsHwRtls9F0Lvb2C+PmiZrYRLwzijwfEiY5jI3+nc27DxYoWkdSrLitk2x3XAzAwOMTnf/Abdu7XMxAyxWTvWloPdCeONLOIcy6a0O0zLma2HrjLe92ovQOR9JWbE2L14lk89MIRuk73UVmS73dJMkWTPYBcmWRcFxAez5u97qfVXgAMDzcDj5hZ+yRrEpEUilSVABDV7Soygi/PM/C+9c9KGF56sfd53U9NAFVVVbS1tc1UiWmlp6dHbc1AQW/r8TNDAPzjw7vouXTsPYOgt3MigtzWyYZBsr2ASpJ0HU0nr/upBaCurs41NDTM5OrSRltbG2pr5smEtv7bsV38fH8Xf/PpG8kZ44lomdDO8QpyWyfbTdRBkq4i51x0auWISJC8/8q5dPf289Lhk36XIlM0qTBI3Oh7Zxdtm5aKRCQwrlkc+0745Ks6qyjoxnOdQb33uhl4OO4snw1mtgnvojPn3JSuMRCR4JkXLqK6rIAXD2nPIOjGc51BK7AlybRo3HidBiqSpZZVl+oBOBlAdy0VkSlZWhULA+ec36XIFCgMRGRKllaVcOrsAEdOnvW7FJkChYGITMm1S2YD8NjLx3yuRKZCYSAiU3LF3DIWzS7my/e/wCtHT/ldjkySwkBEpsTM+Mbv1TPkHN94bK+OHQSUwkBEpuwd8yvYsHohP3n6EHc/utfvcmQSFAYiMi2+fMtybl5Ry92P7qHrdJ/f5cgEKQxEZFoU5ObQdGOE/kHHE3vf9LscmSCFgYhMmyvnV1BWmEvbbp1ZFDQKAxGZNrk5IX57RS0PPX945NGYEgwKAxGZVh+7bhG9/YN85ru7OK1ACAyFgYhMq6sWhvlfG69m5/4ubv7bX/DY6/2c7R/0uyy5CF+edCYime3WlfMpzs+h+aGX+e6LZ/jhf/kZKxeGWb14FuGifIrycyj2fuaHi7m0ppTCvBy/y85qCgMRmRE3rajlfVfU8PX7HqGneB6/jr7F37XtI9k1afm5IdZdUcN7L6/md95ZS3G+Nk2ppk9cRGZMTsi4qiqXhoblAAwOOXr7BjjTN0hv3yCn+wbYd+w07fu7eODZwzz43GG+9+v9bLvjegpytaeQSgoDEUmZnJBRVphHWWHeyLgV8yq45ap5/MUHlvOTpw/xpe3P8K3H9/PZhqU+Vpp9dABZRNJCbk6If79qAZfXlvHkq2/5XU7WURiISFq5vLaMV47o7qeppjAQkbRSV1vOoRNnOXGm3+9SsorCQETSyspLwgA89nKnv4VkGR1AFpG0cu3iSiJzSvjL+1/g0Ikz3LyiltLCXArzciiPO/As00thICJpJRQytn5sFf/9X15iy0O72fLQ7pFpt61eyEevW8S8cCGzSwt8rDLzKAxEJO1cWlPGtz91Le2vHef1rtOcPjfIC4dO8sOdr3PvrjcA+MP3LuOPb6rzudLMoTAQkbS1atEsVi2aNTJ8x40Rdh85xYPPHeZrj+7l5NkB/ux3ryA/V4c/p0phICKBsWh2CYtml/C+K2ooK8zlO7/aT1F+Dptvvtzv0gJPcSoigZMTMr5y6zv5d1fP47u/2s+pszoNdaoUBiISWB+/fjG9fYO0vnTU71ICz5cwMLOwmTXFDdeb2T4za/d+mv2oS0SC5eqFYYrzc3jmjRN+lxJ4fh0ziAC3Ay3ecKVzbinEggHo9qkuEQmQnJCxfG45zx1UGEzVRfcMvG/t7UnGR8xsk5k1er/D412pc64D6Iobbo2bHHHORce7LBHJblcuCPP8wROc6dPT1KbigmFgZo3ey/okk7c657Z4G/IdwJS7dsxsk3Nux1SXIyLZo/GKas4NDPH3P9/H4RNn/C4nsC7YTTT8jd3MRo03s0jCfFEz20is64f44wFxogl7AMlcc7GCRUTivSsym+Vzy/naI3v42iN72HRzHZ9rWOZ3WYEz2WMGSfv1zSzinIs651rOf8uFeQETnmQ9IpKlckLGA19Yw0uHT3L3o3v465+9wodXLqC2otDv0gJlsmcTVSYZ18U4N+Ze99NqM1ufMEnHCkRkwnJCxjvmV/Bnv3sFQ85x89/+gm1PvcHzB0+w71gPb3T1cuB4L4e6zzAwOOR3uWnJXLKnUyfOZOaccxY33ASsc85tiBu3D9jgHRyeEd56mwCqqqpWbdu2baZWlVZ6enooLS31u4yUUFszT6rb+eyxAX68t59XTyTf6Jfnw4bL8nnPgum/A2q6/03Xrl3b7pxbnWzaZLuJku0FVDLDp4R63U8tAHV1da6hoWEmV5c22traUFszT7a0NdXtbAA+P+TY9dpxunv7ONM/SP+gY8g5+gaGuP/pQ/zD811cX38ljctrpnXdQf6bTjYMOkjSVaRTQkUkHYRCxrVLkvVmw23XLGTtV9v42qN7WHlJWLfC9kzqmEHiRt87+JsdfTYiEmh5OSG+dNNlPHfwBKv/WysfvPtxvvPEq4ynyzyTXXDPwDvQW++9bgYejjs9dIOZbSK2l1DvnLt9RisVEZkmH1q5gLqaclpfOsojL3fy5Qde5JGXO/mLDyznspoyv8vzxXiuM2gFtiSZFo0bf7HrB0RE0sryeeUsn1fO59cu45u/jPJ3bfu480fPcd9n3+13ab7Q8wxEJKuFQsbtNy7FAf/jX19m37Eellal7xlBM0W3sBYRAdavWkBhXoi/b9vndym+UBiIiABzSgv40MoF/Mtzhzk3kH03vVMYiIh4blpeQ2/fIE9Guy4+c4ZRGIiIeOoXzQJg95FTPleSegoDERFPeWEuRXk5HD151u9SUk5hICLiMTOqywvoPHXO71JSTmEgIhKnuqxAewYiItmuuqyQY9ozEBHJbtXlBRzsPsNjuzsZGsqe+xUpDERE4tx69XzKi/L41LefYtVXHuZL257JiusOdDsKEZE4Vy0M88Tm9/LQC0doe7mT+zoOMC9cyJduqvO7tBmlPQMRkQT5uSFuuWoef3Pb1aytq+KHO9/g+Ok+v8uaUQoDEZEL+OQNS+g6fY533fUIj7x01O9yZoy6iURELuDGy6p48A/fw3+892n+6J+e5ouNl1JelEdBboj8nBD5ud5PToh93YMsees088JF5OUE67u2wkBE5CKumFvOP3zyGjb+n1/zlQdfuuC8//X/tZETMi6tLuXTa5bwjvkVVJcVMKs4n1DIUlTxxCkMRETGYX64iF9uWsuJM/30nBugb3CIvgHvx3u96zdPM29JHa939fLPTx/kT3Y8O/L+soJcvrrxKn57Ra2PrRibwkBEZJxCIWNWST6zSvKTTu8/kEvD6oUAfPF9l7Kns4fosdN0njrLvU+9wX/a/gzXRWZTUZSXyrLHJVidWiIiAZGbE+KKueW8/8q5fOqGJXx1w1WcOjvA/c8c8ru0pBQGIiIpsGJeOdVlBbTvT89nJSgMRERSwMyov2QW7a8fx7n0u82FwkBEJEVurKvija4zPL73Tb9LOY/CQEQkRT5cP59ZxXn882/S77iBwkBEJEUKcnO4ckGYFw6d8LuU8ygMRERSaMW8cvZ29nC2P73uhKowEBFJocvnljMw5Nj/1mm/Sxkl5RedmVkYWA2EgWuAu7xJo8Y557pTXZuIyEybW1EIwNGT57g8jS5G9mPPYCMQds7tiBtONk5EJOPUlA2HQXo9Z/miewZmVg980zm3KmF8BFgPdAD1QMt4vs0751riBiPAw8651sRxFy9dRCR4qssLAOgMUhiYWSPQRWxjn2irc26dN18UaAZuH++Kve4i4oMg2TgRkUxSmJdDRVEeR0+e87uUUS4YBsMbZbPRt1319gri54ua2Ua8MDCzpiSLiyZs5O90zm1ImCfZOBGRjFJTXhC8bqIx1APdiSPNLOKciyZ0BZ3HzNbjHTg2s0bnXGuycZOsTUQkrV1SWcLeYz1+lzHKZA8gVyYZ10XsbKAL8rqemoFHzKx9rHEiIpmqflGY6LHTdKXRc5VTfmqp941/aZJJycaN4nU/NQFUVVXR1tY2vcWlqZ6eHrU1A2VLW7OlnTD+tuZ0xS44+96Dv+Dq6vR4rMxkq0i2F1BJkq6j6eR1P7UA1NXVuYaGhplcXdpoa2tDbc082dLWbGknjL+t7+ob5H/u+jf6yhfQ0HD5zBc2DpPtJuogSVeRcy46tXJERDJfUX4OK+aVs+u1436XMmJSYZC40ffOLto2LRWJiGSBlZfM4rkDJxgaSo9nG4znOoN673Uzoy8Q22Bmm/AuOnPOjfsaAxGRbLesupQz/YMcPnmW+eEiv8sZ13UGrcCWJNOiceN1GqiIyAQsrSoFYF9nT1qEge5aKiLig2XVsTDY25ke1xsoDEREfDCnNJ/ywlz2pcnFZwoDEREfmBlLq0sVBiIi2W5ZVSn7jqXHQ24UBiIiPllaXcqxU+fo7vX/thQKAxERn6xaNAuAJ/a+5XMlCgMREd+sXBgmXJzHjzoO4Jy/F58pDEREfJKbE6LptyI88nInf3rfc/QNDPlXi29rFhERPnvjUs70DXL3o3t5+egp/vz9V3DN4mRPCZhZ2jMQEfGRmfGlm+r437ddzYGuXj7xrZ08e6A75XVoz0BEJA3cunI+7146m1u/8QS3fP0JLq8t44Zlc/jgVfO4emF4xtevPQMRkTRRXV7Ijz53A5tvvpzKknz+8dev8eG/e4L217pmfN0KAxGRNFJbUchnG5byg89cxy83r2XIwZOvKgxERLJWTXkhVWUFvJqCq5QVBiIiaWzJ7BL2v6UwEBHJakvmlLC3s4fBGX4imsJARCSNveeyORzv7ecXrxyb0fUoDERE0thNy2uZW1HIn//z8/zxvU/z/METM7IehYGISBrLzw3xtY+spKIoj6de6+LEmf4ZWY8uOhMRSXPXLK7kX774nhldh/YMREREYSAiIgoDERFBYSAiIigMREQEhYGIiKAwEBERFAYiIgKYczN786OZYmangN1+15Eic4A3/S4iRdTWzJMt7YT0b+si51xVsglBvgJ5t3Nutd9FpIKZ7VJbM0+2tDVb2gnBbqu6iURERGEgIiLBDoMWvwtIIbU1M2VLW7OlnRDgtgb2ALKIiEyfIO8ZiIjINAnU2URmFgHWAx1APdDinOv2tagpMLN64JvOuVUJ48dsZ1A/A6+tjd7gNcBnxtOmILbXzIbbGSbW1nudcx3etIxqazwz2+qcuz1uOKPaambNwD5gG7AR6HLO7fCmBb+tzrnA/AAPx72OAFv9rmkKbWkk9g/DTaSdQfwMiG0Um+KG1wPtGdze40A4G9oaV+95/5Yzra1AM+C8v++mTGtrYLqJvHQd4ZyLEkvnQHLOtTrv22K8C7UzwJ/BamBz3HArUG9m4Qxt7yr39je/SqALMvZvOywCdA8PZGhbn3LOmXNulnNuy/DITGlrYMKA2DeP7sSRiR92BrhQOwP5GTjnWoENcaMi3vhuMrO90bjBDcS+UUIGthXAzNY7r7skTka2FUa6PONlRFuDFAaVScZ1EeuCyCQXamdgP4OEvaDbgOFvVhnZXjOLmNkmYLsXhpCBbfU2atEkkzKurUDEOx4UNbPmuFDIiLYGKQwkA5hZGKh3zm2+2LxB5u0dtACrzGy93/XMoPpk3Z2ZyDm3xeve7Qa2Att9LmlaBSkMkqVpJUl2wQLuQu3MhM+gmdFdRhnbXm+jsR3Y7oVgRrXV+5bcOsbkjGorjHyRAUbCfrirJyPaGqRTSztIssuV0D+bCcZsp5kx1rQU1DVlXrfJZudct5mFvY1lRrXX20A2u7dPF97l/a4kw9rq2ejVDhA2syZiAZFRbR3+uwKrkkzOiLYGJgziPlhgpK9ym38VzYwLtTPIn4HXVbIj7iybRm8409rbBdwbN7waiA7/z59JbY07FgKMXGfQEjccPy3QbSUW6ncNDwz/e4bM+X82ULejSLx4I/70rqDxvmnUE/u2sYXYucit3rQx2xnEz8CreV/C6Khzbmnc9Exq73re/ja4itiewnAYZFRbYaT7pIm3/y1v9TaCGdXWuAsnu4Gl8ce9MqGtgQoDERGZGUE6gCwiIjNEYSAiIgoDERFRGIiICAoDERFBYSAiIigMJAuZWXvcA2iybv0iyeg6A8k6ZnbezdXibo8xE+sbtexk6xfxm/YMJOskCYIIbz+Sc1olW7aCQNKRwkBk9FPYgrRskWmjMJCsYmb1Xp99kzfcSOxWxOvMrGn4NsVm1ugNr/cehD48bp83vsnMtseNr/fGbYpb13nLTlx/XE1NceuMJNYaN21rij4qyTKBuWupyHRwznWY2b1xw61mtoHYjQJ3wEjXzmbn3DpvuNLMmpxzLWbWSuwZx7eb2fDtqbd64zrMbKuZNXoPQTlv2cCo9Xvrah5elzeuPW55w+trGZ5mZpF0vAWyBJvCQOR864HuhGfdDt/HvhvvDqxxff+rvGc0RIjdrTQ8wXU9nDCuazhQgLcYfcfXtHxkogSfwkAkjrdBn03sFtvDG/vEA75dCcN3mtlbxO5vP+Y39jG+0c+eSr0i00XHDERGP4IwQuzhNPF7BcP9/+cZfi6F93zcKN5D0OPmT1x2oq2J6wIqEx8cIzLTtGcgWcXr+rmNWFdMq7cB3wps9p5I1eE9mKXZOxjcQWwD3xr3cJOImXV7G+xdxLqU1hPb8G8n9ozn+OMJI8tOtn5vXU3E9irqvfcn1rqLWJisBm43s80zdV2EZCdddCYiIuomEhERhYGIiKAwEBERFAYiIoLCQEREUBiIiAgKAxERQWEgIiLA/wcOdV6cK1uCnQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.semilogy((LCGsolver.obj / LCGsolver.obj[0]))\n",
    "ax.set_xlabel(\"iteration\")\n",
    "ax.set_xlim(0)\n",
    "plt.suptitle(r\"$\\phi_{LCG}$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEoCAYAAACuBsGbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAjo0lEQVR4nO3deXxU9b3/8dc3G0kIZoEISGQJSGSXJKLVFqENtS5dtCC1tdVaDe2tt+39VaFU2z6qtRbu7ertbcFudqMQqrZurURNtXYRElaBgAQ07IEQspFt5vv7Y05wHMOWmeTMmXk/H488mHMymfP5Dsk733zP93yPsdYiIiLelOB2ASIi0nsKcRERD1OIi4h4mEJcRMTDFOIiIh6mEBcR8TCFuEQtY0yJ2zWIRDuFuESzLLcLEIl2CnEREQ9TiEtUMsZkAQ0ulyES9Ywuu5do4oyD5wNFwBrn8XJrbYObdYlEK/XEJWoYY5YA+dba5cAaa+1qYDlQ5m5lItErye0CRACMMaVAlhPgJ1lrG4wxOcaYfGttTR8evxB4BKgBVjq/QEKfM5fAXwYNzq4aoNBau7Q3zxOJBIW4RIslwJig7Yagx1l9fXBrbZUxBmBRT78sjDHLCPx1sDRk366Q55UR+CVw2ueJRIpCXFxnjMmHQK/b2S4B1jmPs4CcvuyFB+mxt+8M8+zqoXdeBtSHPG/tKZ7XH/VLHFKISzSoJygMCQyrNDiPFwOLgp9sjFkIVDmb64AS4FJr7SIn9EuB1c5+eCtAC4O+rjCkt3zyF0fIsbKAhdZaE/o5a235uT5PJNJ0YlNc5wR2TXePvJsztvy2cXJnX40TjHOAYgLB3P21JQSGYrIIhHeR89x1wHxrbbmzPT+kjDkEZsOEKgbOJoRLQp9njCkxxsx1PrLO4jVEzplCXKLFPGCuE9JznBOdWGsXhDyv3Pl899hzOTAXWOZ8fj5Qbq2tIhDM3TNbioGVcHL4JnR4420hbIzJcnrnWZxivnrIsgANoc8L+mWRoymS0lc0T1yiijNLpOEUY9MlOMMgzvNKnMdl1tp5znMqrbVFPTxeBixzTmAuJBDiNU7YY4w5Zq3NDjrWQue1s4Dnu18n6POloTNpgo8XtO8YMEYhLn1FY+ISbfJ7mt7n6D7ZWQgUB41pr3R68Dm8fVw7+BdBTndgE+gx5wBVTq98EVAf9BpzcMbonSmO87pPWjpf+46pkI7u5+1yjp1FYLZLw1m2XeScqScuUaWnHu65fC1Qf5pfAiIxR2PiEm3qz/yUtzgnDxc6m3MU4BJv1BMXT3OGQ7o/yvtpPrlI1FCIi4h4mIZTREQ8TCEuIuJhCnEREQ9TiIuIeJhCXETEwxTiIiIedsYQN8YUGmMqe9ifb4xZ2H2xhVZpExHpf6edJ+4sOFQPVIauk2yMWWOtneM8ziewRkToinMiItKHzupiH2OMDQ5xJ7SXdYe4s+9tq8CJiEjf6+2YeCE9rLEcuqi/iIj0rd6GeE4P++rphxvaiojIW/plPXFnidBSgNTU1KKRI0f2x2Gjjt/vJyEhPicEqe1qezyKZPt37NhxxFqbG7q/tyHeU687h1PcxspZH3o5QEFBga2uru7lYb2toqKCWbNmuV2GK9T2WW6X4Yp4bjtEtv3GmDd62t/bXxFV9DCkomVARUT6V69CPDSsnROaqyJSkYiInLXTDqd035jWebwEWOPcwRsC9xNcSKBXXqg54iIi/e+0Ie4EdjmwtIfP1QTtLw/9vIiI9L34PW0sIhIDFOIiIh6mEBcR8TCFuIiIhynERUQ8TCEuIuJhCnEREQ9TiIuIeJhCXETEwxTiIiIephAXEfEwhbiIiIcpxEVEPEwhLiLiYQpxEREPU4iLiHiYQlxExMMU4iIiHqYQFxHxMIW4iIiHKcRFRDxMIS4i4mEKcRERD1OIi4h4mEJcRMTDFOIiIh6W5HYBIiKn0uXzc6LTR1unn7ZOH22dPtq7Ao87uvy0d/lp7/LR4bN0dvnp8Pnp9Pnp9Fm6fH66/JYun6XLH3jsC/rw27f+9fsJ/GvB2sA+Cye3bXdBtvsfi7WnKDpIXV0bK/dWAvCR6SO4etKwiL9HCnERiThrLa0dPo6c8LNl33GOn+jk+IlOGk900tjWSeOJLpraOmlq76KprYuWduejw0er8++JDh8dPn9E6jEGkhIMiQmGRGNICH1sDAkGjDEkJIDBYAwkGIMBMAT+JfAcZxfGnOKAjpYWP400A9DQ2hGRtoRSiIvIWfH7LfWtHRxubOdQUxt1Te0caW6nrqmdo80d1Ld0cKS5nfqWDhpaO98K4L/9/R2vlWBgUGoyGQOSGJSaxMABSWSmpzAiO5H0lCTSUwL/piUnkp6SSGpKIqlJCaQmJzIgKYEByYHtFOdjQFICKYmJJCUaUpISSE5IICnRkJRoSE5IICHhDGnbRyoqKpg166o+PYZCXEQAOH6ik73HWtl77AT7G05w4Hgb+xpOcPB4GwePt3G4qY1O3zvHEDIGJDE4I4WcgSnkZaczNS+T7IEpZKencKi2hsunTyEzLfnkx3lpyQxMSTzZo5XwKMRF4oS1lkON7ew+0sKeo4GP2vpW3jjaypv1rTS1db3t+QOSErggK43hmalcNiaHoZmpDB00gKHnpXL+eQPIzUgld9AA0lIST3nMiopaZvXBOLC8RSEuEmM6fX7eONrKzkNN7DzczK66wEdNXQutHb6Tz0tJTCAvJ42ROekUjszmwpw0LsxOZ0R2GiOy0sgZmKLesgeEFeLGmHygBKgH8oHl1tqGCNQlImehrqmdbQca2Xqgke0HGtl+sImaupa3nRDMy04jPzeDS0fnkD9kIKOHDGTMkIEMz0wj0aWxYomccHvic621S7s3jDFLgEVhvqaI9OBQYxsbaxvYsu84m/cd57X9jRxuaj/5+eGZqRQMG8RVBbmMP38Q44cOYuz5A0lP0R/csSzc/905wNKg7awwX09EgBMdPjbtbaDqzQY21B5jY+1xDja2AYGZHePOz+DdFw1h0gWZTBg+iInDzyMrPcXlqsUNYf+KNsasAeYBxUBZ2BWJxKEjze2s21PP2j3HWLenntf2N9LlD8wEGTU4ncvyc5iWl8W0CzOZMPw89a7lpHC/E+YBzwO7gYeCh1ZE5NSONrfzz5qj/KvmKP+qqef1w4ELQgYkJXDJhVmUzsynaFQ200dmkzNQPWw5NWPP5trRU32xMSXOw3xgGbCopyA3xpQCpQC5ublFq1at6vUxvay5uZmMjAy3y3BFvLc9OW0g1fU+XjvqY+tRP7VNgROPqYkwPjuRgpwECrITGZWZQHIMnWyM5/93iGz7Z8+eXWmtLQ7d3+sQd2amnDyx6WxXAmNON0OloKDAVldX9+qYXhe4emuW22W4It7abq1lV10LFdWHefxf1ew8Dh1dflKSEigelc2V44ZwxdjBTBmRSVJi7K5DF2//76Ei2X5jTI8hHs5wSglQ3r1hra0xxiwHcoCGMF5XxJM6uvy8urue8m2HeH77IWrrTwBwQYbhk5ePZub4XGaMzjntxTEi5yqcEC8H5gJVwTuttTVhVSTiIc3tXby4/TB/fe0gf6uuo6m9i9TkBK4cO4QFM8cyqyCX1ze+yqxZE90uVWJUr0Pc6XnXGGMWAjUEeuDLIlaZSJQ6fqKT8q2HeGbzAV5+/QgdXX4GD0zh2inDmTNxKFeOG/K23vbrLtYqsS+s2SnW2tWRKkQkmrW0d7Fm6yGe3Lifl3bW0emzjMhK45bLRvGBycMoGpWtqx/FFZpsKnIKnT4/L+2o4/H1+yjfdoi2Tj/DM1O59V2juW7qcC65MEtri4jrFOIiQay1vLa/kdWVe3ly436OtnSQnZ7M3KI8PjRtBMWjsl1bm1qkJwpxEeBYSwePr99HWeVeth1oJCUpgZIJ53Pj9DyuKsglOYanAYq3KcQlbllr+WfNUVa8Wstftxykw+dnal4mD3x4Eh+aNoLM9GS3SxQ5I4W4xJ3jrZ2UVdby+3+/Sc2RFs5LTeLjl41k/qUXMmH4eW6XJ3JOFOISN7bub+TRf+zhTxv30dbpp2hUNt977ziunTKc1GRdgCPepBCXmObzW9ZsPcgvX9nDv3fXk5qcwA3TR3DL5aOYdEGm2+WJhE0hLjGppb2LsnW1/OKVPbxZ38qIrDS+eu3FzC8eqbFuiSkKcYkpR5rb+dUre/jNv97g+IlOCkdmsfiai5kzcWhMLzQl8UshLjGhtr6VZS/tomzdXjp8fq6eOIw7nTW5RWKZQlw8bVddM//34i6e2LCPRGO4sXAEd87MZ2xu/K5hLfFFIS6e9PrhJn70/Os8uWk/A5ISuPVdoymdmc+wzFS3SxPpVwpx8ZRddc38sHwnT27aT1pyIgtmjuWO94xhSMYAt0sTcYVCXDyhtr6VHz6/k8eq9pLqhHfpzHzdf1LinkJcotqR5nYefn4nv3/1TYwxfPrKMXxu1lj1vEUcCnGJSs3tXfzs5RoeeamGti4/NxVfyBfeN47hmWlulyYSVRTiElW6fH5Wrqvl+2t2cKS5g2unDOPu9xeQr9kmIj1SiEvU+NuOOh58eis7DjVz6ehsHvlUMdNHap63yOkoxMV1NXXNPPDUVl6srmPU4HR+8olCPjB5mO6aI3IWFOLimsa2Th5+fie/fGUPqcmJLL7mYm67cjQDkrSioMjZUohLv7PW8sSGfTz49HaOtrRzU9GF3H11AbmDNONE5FwpxKVfVR9s4mtPbOHVPfVMy8vkF7cVMzUvy+2yRDxLIS79ot1n+c6z2/nZyzUMSk3iOzdO4abiC3XTYZEwKcSlz71YfZh7/36CIyd2cVNxHl+5ZoKutBSJEIW49Jmjze3c/9RW/rRhPxcMNKwsvZzL8ge7XZZITFGIS8R1n7i8/8mtNLd38aWSi5iUsE8BLtIHFOISUQePt3Hv45t5fvthCkdm8Z2PTmX80EFUVOx3uzSRmKQQl4iw1rK6ci/3P7WVTp+fr10/kduuGE2iTlyK9CmFuITtcFMbX31sC+XbDjFjdA5L505l9JCBbpclEhcU4hKWZzYf4N7HN9PS4eO+6yZw+5VjNG1QpB8pxKVXmto6+cafX+Oxqn1MzcvkezdNY9z5g9wuSyTuKMTlnK3dU89/rdzA/oYTfPF9F3HXe8eRnJjgdlkicUkhLmety+fn4Rde5+EXdpKXnU7ZZ6+gaJSWihVxU9ghboyZC+QANQDW2vJwX1Oiz76GE3zpD+tZu+cYNxaO4P4PTyZjgPoAIm4L66fQCfB8a+1SY0w+UAYURaQyiRrPvXaQe1Zvosvn5wfzL+Ej00e4XZKIOMLtSi2x1o4FsNbWoACPKR1dfpb8ZTs///tupozI5OGbp2vqoEiU6XWIG2MKgYbux9baqkgVJe6rrW/lrhXr2VjbwG1XjGbxtRfrZg0iUSicnng+UO8MqZQbYxYCVRoT974Xtx/mSys34LeWn95SyAcmD3e7JBE5BWOt7d0XGlMKLLPWGmc7C9htrX3HdAXnuaUAubm5RatWrep1wV7W3NxMRkb03rXdby1PvN7Jn3d1MnJQAndNH8D56ZGZOhjtbe9Lant8th0i2/7Zs2dXWmuLQ/eHE+JzgcXW2qKgfRYY64yP96igoMBWV1f36pheV1FRwaxZs9wuo0cNrR3854r1vLzzCPOK8njgI5NJTY7c8Ek0t72vqe2z3C7DNZFsvzGmxxAPZzilpzHwBqA+jNcUF2zd38iC367j0PF2HrpxCjfPGOl2SSJylnod4tbaGmNMvTEmy1rb4EwxrLHWNkSuPOlrf964n4WrN5KVlsLKBZczfaQu3hHxknCnGM4DFhtjdgFjgfeFX5L0B7/f8t/PVfOTil1cOjqbH3+ikPMHpbpdloico7BC3Ol1L4pMKdJfmto6+a+VGyjfdpibZ4zkmx+aREqS1j4R8SJdNx1n3jzaymceXUvNkRYe+PAkbrl8FMZo6VgRr1KIx5G1e+pZ8JtKfH7Lb26fwRXjhrhdkoiESSEeJ/5YuZfFj20mLzuNn992KWN0+bxITFCIxzhrLd9fs4MfvfA678ofzE9uKSQrPcXtskQkQhTiMayjy89X/riJx9bvY15RHg/eMEUnMEVijEI8Rh0/0cmC36zjXzX1fHnOeO567zidwBSJQQrxGHTg+Alu/cWr7D7SwvfnT+OG6XlulyQifUQhHmOqDzZx2y9fpamti0c/rRkoIrFOIR5DXt1dzx2PriU1OZFVC97FxAvOc7skEeljCvEYsWbrIe76fRV52Wk8evsM8rLT3S5JRPqBQjwGlK2r5SuPbWbyiEx+edul5AzUFEKReKEQ97hHXqrhwWe28Z6LhvDTW4oYqDvQi8QV/cR7VPBFPNdNGc735k/TPTBF4pBC3IP8fssDT2/ll6/sYX7xhXz7xikkJmgOuEg8Uoh7jM9vWfzYJlat28vtV47ha9dP0EU8InFMIe4hXT4/Xy7byJ827OcL77uI/yq5SAEuEucU4h7R0eXni39Yz7NbDnLP1QV8fvY4t0sSkSigEPeA9i4fn/9dFeXbDvO16yfymXePcbskEYkSCvEo197l43O/reKF7Yd54MOT+OS7RrtdkohEEYV4FGvr9PG531byYnUdD94wmU9cNsrtkkQkyijEo1Rbp4/P/raSiuo6vn3DFD5+2Ui3SxKRKKQQj0IdXX7+43dVVFTX8dCNU7h5hgJcRHqm27xEmU6fn7t+HxgDf/CGyQpwETkthXgU6fIFphE+t/UQ3/zQJI2Bi8gZKcSjhM9vubtsI89sPsh9103g1itGu12SiHiAQjwKWGu574nNPLFhP/dcXcAd78l3uyQR8QiFuMustTzw1DZWvFrL52eP1ZWYInJOFOIu+375Tn7xym4+feVo7n5/gdvliIjHKMRd9PO/7+ZHz+/kpuI8vn79RC1mJSLnTCHuktWVe3ngqa1cM3kYD904VQEuIr2iEHfBc68dZNEfN/Gei4bwg49dohs6iEivKcT72b9rjnLXivVMGZHJT28p0i3VRCQsEQtxY8yySL1WrKpt8nPHr9dxYXYav7ztUt3UWETCFpEQN8YUAqWReK1YVVvfynfXtTEwJYlff+YysgemuF2SiMSASHUF84GGCL1WzKlv6eDWX7xKp9+y6jMzGJGV5nZJIhIjwu6JG2PmWmtXR6KYWHSiw8cdj65lb8MJvliYyvihg9wuSURiSFghbozJB2oiVEvM8fktX/zDetbXNvDD+ZcwPlsnMUUksoy1tvdfHNQLN8Ycs9Zmn+J5pThj5rm5uUWrVq3q9TG9wlrLb7d18PybXXzi4hTmjE6mubmZjIwMt0tzhdqutsejSLZ/9uzZldba4nd8wlrbqw+gBMgK2j52Nl83fvx4Gw8eeWmXHbXoKfvAk6+d3Pfiiy+6V5DL1Pb4FM9ttzay7QfW2R4yNdwTmzcFXWmY5fS4y621cT3E8pctB3nwmW1cM3kYX712gtvliEgM63WIW2vLg7eNMcustcvDL8nbNtQ28KWV65mWl8X3519Cgq7GFJE+FInZKVnGmIXO4yXOyc64VFvfyh2PriV30AB+dmsxqck6kSkifSvseeLW2gZgqfMRt5raOrnj0XW0d/n5Q+nlDMkY4HZJIhIHdN13BPj8li+sWM/rdc08+ukZjDtfc8FFpH9oAawIePDpbbxYXcc3PzSJd180xO1yRCSOKMTDtOLVN/nFK7u5/cox3HK57k4vIv1LIR6Gf9cc5WtPbOGq8bnce52mEopI/1OI91JtfSuf+10VIwen86Obp+vGDiLiCoV4L7S0d3Hnr9fR5fPzs08Vk5mW7HZJIhKnNDvlHFlrubtsIzsONfGrT88gPzd+14UQEfepJ36O/q9iF89uOchXr53AzPG5bpcjInFOIX4OXth+iP95rpqPXHIBn3n3GLfLERFRiJ+tmrpmvrhiAxOHn8dDN04laOEvERHXKMTPQkt7Fwt+U0lyUgLLPllEWorWRBGR6KAQPwNrLQtXb2JXXTP/e/N08rLT3S5JROQkhfgZPPJyDU9vPsDCD1zMFeN0Sb2IRBeF+Gn8Y9cRvvPsdq6ZPIwFM+N2hV0RiWIK8VM4eLyNL6xYz5ghA/nvedN0IlNEopIu9ulBp8/Pf66ooqXdx4o7LydjgN4mEYlOSqceLP3LdtbuOcYPP3YJFw3V2uAiEr00nBLiL1sO8MjLu/nk5aP48CUj3C5HROS0FOJB3jzayj1lm5ial8l912tpWRGJfgpxR3uXj8//vgpj4McfL2RAki7oEZHopzFxx0PPbGfzvuMs+2QRF+bogh4R8Qb1xIFnNx/gV//Yw+1XjuHqScPcLkdE5KzFfYjX1rey8I+bmJaXyVeuudjtckREzklch3hgPvh6AP7344WkJMX12yEiHhTXY+LffW4HG2ob+PHHCzUOLiKeFLddz5d21PHTv+3i5hkjuW7qcLfLERHplbgM8bqmdv7fqg2MH5rB16+f6HY5IiK9FnfDKdZa7lm9kaa2Ln53x+W6wYOIeFrc9cR/9Y89VFTXcd91EygYpnVRRMTb4irEtx1o5KFntlMy4XxuuXyU2+WIiIQtbkK8rdPHF1asJzM9mSUf1Y2ORSQ2xM2Y+EPPbGPn4WZ+ffsMBmcMcLscEZGIiIueeEX1YR795xvcfuUYZo7PdbscEZGICasnbowpBEqczUuBO621DeEWFUn1LR3cs3oT44dmsPADBW6XIyISUb0OcWNMFlBsrV3qbM8FngeKIlNa+Ky1LH5sE8dbO3n00zNITdZ0QhGJLeEMpxQDi4K2y4FCJ9yjQlnlXv762iG+/P7xTLzgPLfLERGJuF6HuLW2HJgXtCvf2d8QZk0RUVvfyv1PbmXGmBzueE++2+WIiPSJsE5sWmurgjbnA0vDKycy/H7L3WUbsdby3XnTSEzQdEIRiU3GWhv+iwSGUMqstXNO8flSoBQgNze3aNWqVWEf83T+uqeTFds7uH1yCjPzkvv0WOeiubmZjIwMt8twhdqutsejSLZ/9uzZldba4tD9kQrxZcCisxlKKSgosNXV1WEf81R2Hmriuof/znvGDeFntxZH1UU9FRUVzJo1y+0yXKG2z3K7DFfEc9shsu03xvQY4mHPEzfGLMQJcLdPanb5/Hy5bCMDUxJ56KNToirARUT6Qlgh7kwrXB3UAy85zdP73LKXati09zgPfGQy5w9KdbMUEZF+Ec488XygzHncvbsGWB1+Wedu+8FGflC+g+umDOf6qRe4UYKISL/rdYhba2uAqBiv6PT5ubtsI+elJnP/hye5XY6ISL+JiQWwflKxiy37GvnJJwq1uJWIxBXPL4C1/WAjD7+wkw9Ou4BrpuhemSISXzwd4l0+PwtXb+K81GS++SENo4hI/PH0cMojL+9m097j/PjjheQMTHG7HBGRfufZnvjrh5v5fvkOPjBpGNdOGeZ2OSIirvBkiPv9lkV/3ER6SiL3f2SSLuoRkbjlyRD/zb/eoPKNY3z9+om6qEdE4prnQnxfwwmW/mU7V43P5YbpI9wuR0TEVZ4KcWst9z6+GQs8eMNkDaOISNzzVIj/acN+KqrruOfqAvKy090uR0TEdZ4J8fqWDr755GtMH5nFp9412u1yRESigmdC/FtPb6W5vYslH52qO/WIiDg8EeJ/33mEx6r2sWDmWMYPHeR2OSIiUSPqQ7yt08e9T2xm9OB07nrvOLfLERGJKlF/2f3DL+zkjaOt/P6Oy0hNTnS7HBGRqBLVPfHqg00s+1sNHy3M44pxQ9wuR0Qk6kRtiPv9lvue2Myg1CTuvW6C2+WIiESlqA3x1ZV7WbvnGIuvnaAVCkVETiEqQ7y+pYNvP7uNGaNzmFuY53Y5IiJRKypD/KFnttHc1sW3bphMguaEi4icUtSF+Ku76ymr3MudM/M1J1xE5AyiKsQ7fX7ue2IzI7LS+MJ7L3K7HBGRqBdVIf6rV/aw41Az3/jgRNJSNCdcRORMoibEDx5v4wflO3jvxeczZ+JQt8sREfGEqAnxbz29lU6/5RsfnKh1wkVEzlJUhPgrrx/hqU0H+I9ZYxk1eKDb5YiIeIbrId7R5efrf9rCyJx0PnvVWLfLERHxFNdD/NF/7GFXXQvf+OBELXAlInKOXA3xw41vncx83wSdzBQROVeuhvhDz26n02f5+vUT3SxDRMSzXAvxtXvqeXz9Pu6cOYbRQ3QyU0SkN1wJcZ/f8o0/vcbwzFQ+P1t36xER6a2w7uxjjMkH5gJVQCGw3FrbcKavW7m2lq0HGnn45umkp0T9zYVERKJWuAm6zFo7B8AYUwMsARac7gv8Fv7nuWpmjM7h+qnDwzy8iEh86/VwitMLP8laWwPcdKava2i3HGvt4Ou6MlNEJGzhjIkXAg2hO0PDPVRjh+Vjl45k8ojMMA4tIiIQXojn9LCvHsg60wHvfv/4MA4rIiLd+uWsojGmFCh1NtuHDErd0h/HjUJDgCNuF+EStT0+xXPbIbLtH9XTznBCvKdedw49DLFYa5cDywGMMeustcVhHNez1Ha1Pd7Ec9uhf9ofznBKFT0MqTgnOEVEpB/0OsRDw9o5obkq7IpEROSshTsmPs8YsxDnYh9r7WnniDuWh3lML1Pb45PaHr/6vP3GWtvXxxARkT7i+nriIiLSe/01xbBXa6x4lTGmEChxNi8F7uxubzy9F8aYZcFDbPHQdmPMXAIn/GsArLXlzv54aHs+ge/7eiCfoDbGWvudn/FHrLVFIftP2c4+ew+stX3+AawJepxPYM2Vfjl2f38QmHZZGrQ9F6iMt/fC+Sa18fR94PxfLwxqX1z9v3e3PWh7SSy2n8Avqnd8f5+pnX31HvT5cEpv11jxsGJgUdB2OVBojMmKs/cin6BrBuKk7UustUsh0D7r9NLipO0Ac0K2syD22m+tLbfWVoXuP107+/I96I8x8V6tseJVNvDn87ygXfnO/gbi5L0wxsy11q4O2R3TbXf+vG4IehwsptsezBizxumwlABlzu54af/p2tln70F/hHiv1ljxspDf0vOBpc7jmH8vnG/Kni74ivW25wP1zph4jTFmoRNkEPtt7zaPQFt3E5hyXO7sj5f2n66dffYe6I4MfcgYk0Xgmzn0z8xYVthDLzwe5AAl9q319ZcTCLNsV6vqX91DifnAMmMM3cNL0nf6oyd+1musxKAlvH1oJabfC6fnWX6KT8d02wm07+RfYM7wWfd5kFhv+8khA2e8eDkwFljsdGRivv2O07Wzz96D/uiJx+UaK86VrIustQ3GmCznhzoe3oubgm72keWsYFlO7Lf9HSe6eOuHN9bbDoEZGyd/gVtra5y/RnKIj/bDadrp/Ez0yXvQ5yEe1AAgPtZYccZFV9u35oCWONsx/V4EjYECJ+eJLw/aDv5crLW9xhhT3/0Lu/vcgPM90BDLbXeU89Yc6JO6QyoO2n/arOvLn/1+uew+dJJ7LI+TOW3dFbK7xlo7NujzMf1eOH9ClxIYTlpKYD5sTay33Wn3YgL//2OBh+wpLvSItbbDyc5L94ntHKA8KMRjpv3OsGEhb31/r7GnuKgruJ199R5o7RQREQ/T2ikiIh6mEBcR8TCFuIiIhynERUQ8TCEuIuJhCnEREQ9TiItnGGMqgxaVirvji/RE88TFM4wxhaHrOActadAXx3vba/d0fBG3qScuntFDgHffDizienptBbhEI4W4eNmiMz8lKl9bJGIU4uIJxphCZ0y61NkuIbBOxxxjTKmzbgnGmBJne64xZknQvl3O/lJjTFnQ/kJn38KgY73jtUOPH1RTadAx80NrDfrcsn56qyTO6KYQ4gnW2ipjzMqg7XJjzDwCiw+thpNDIIuCbsyQY4wptdYuN8aUA0XW2gXGmHXOyyxz9lUZY5YZY0qc9bDf8drA247vHGtJ8A0/jDGVQa/Xfbzl3Z8zxuTH4PKr4jKFuMSSuQSWfQ2+x2WR828DzuqSQWPbRUHLxuZwbrfKmgusCdlX3/2LADjK21ezjMXbkUkUUIhLTHCCeDCBZX+7Qzr0RGR9yPZiY8xRYDU93xf05Gv30IMeHE69IpGiMXHxsoagx/nASgLrPJ90qnnd3WtCW2uXOgGdBeQEPT/0tUMtCz0WkBN6YwyRvqaeuHiCM0Qyn8CQRffNBpYBi5w7plQ5N55Y4pykrCIQzOXO15YA+caYBido1xEYeplLILDLCNwPNXi8/ORr93R851ilBHrxhc7Xh9a6jsAvgWJggTFmUV/Na5f4pIt9REQ8TMMpIiIephAXEfEwhbiIiIcpxEVEPEwhLiLiYQpxEREPU4iLiHiYQlxExMP+P4kmWN6g6fcjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(SLCG.obj)\n",
    "ax.set_xlabel(\"iteration\")\n",
    "ax.set_xlim(0), ax.set_ylim(0,10)\n",
    "plt.suptitle(r\"$\\phi_{symLCG}$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also, let's compare the two inverted functions with the analytical solution.\n",
    "To find the solution for the continuous case we need three conditions:\n",
    "\\begin{equation}\n",
    "\\frac{d^2f(x)}{dx^2}=1,\\\\\n",
    "f(x=0)=0,\\\\\n",
    "f(x=x_f)=0.\n",
    "\\end{equation}\n",
    "$x = 0$ and $x = x_f$ are not sampled and lay outside of the interval $\\mathbf{x}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.linspace(dx, N*dx, N)\n",
    "alpha = 0.5\n",
    "beta  = -(X[-1] + dx) * 0.5\n",
    "gamma = 0.0\n",
    "f_an  = alpha * X * X + beta * X + gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaIAAAEoCAYAAAAaIPXRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABRFUlEQVR4nO3deVxU1fvA8c8BVNxHcCk1F3CpbLFBy2xTw8r2DLS9vn4TstU2yPb6VgYtvzZLsNWyUmgx23SmRMushMlyS5NxSVtcYFRUQOD8/pjFAUdgYOAOw/N+vXg5M/fMvc+9wjxzlnuO0lojhBBCGCXM6ACEEEI0b5KIhBBCGEoSkRBCCENJIhJCCGEoSURCCCEMJYlICCGEoSQRiaCilDIrpVKUUoVKqXzX4xij46ovpVSG65zi/XhPvOt9GUqphIaMTwgjKbmPSAQjpVQekKu1TjY6FgClVJLWOrOe+8gDUrXW1lqW10AnYAhAbd9XX77OVSmVBsRorRMbIwbRvEiNSASrAqMDqGJ0Yx5MKWUG7Fprh9ba2lhJyMXXuVqA2Y0Yg2hGIowOQIhgp5RKAoxoHnQ09gEPd66NnAhFMyOJSDQJrr6VNNfTiTg/LGOAaK11qqtMAjDF9Xqi1trq6l+yAHYgWWttV0qluJ7H4Kx1ZHvt3w5kcLBWYHE9jnG9D611uldch+yryjYHftbuXLWhZK9j2l0/7uaxWFe5NCAJZ3NfZm2ukdcxkryfe73/kHN1xVPp2D7OH9f2dNfrtYrFFYcdMLm2O+rbBCqaIK21/MhP0P3gTAAZVV6LB/KBeK/X8gFz1TJV3pfi9TgLSKhyHLPrcQKQ59qHGUjz2meejxir21cakOS1zQRo79hrOH+zj/PwdW6WKsepzTVKcZ+b13kn1HCuvuLJqnKcGMBS21hcx02q8v6kqseWn9D/kT4i0ZQU4PzW7d1M5K6NAJ4mpCjXt3g3B4CrdpSgvWotOD9Mk73KmbWzT8amq9QivFW3L6WUCWfy83yz11o7AFstz/NwfNWsHD7KHPYauWJLA6Z6bR+Pn02Prusb730crbUd57V3jwys8f8LSHTF5H5/rj9xiNAgTXOiqbFXee4Aoqq8lokzuSS7PhTnuF6PBxxVhlDHUvmDser+D6e6fcVjQP+Ol+qu0RCczV8O90Zdt5FwQ3wcx33s0YDV67nPWLSzSTQZKFRK2YDZ2qvZUzQfkohEKMrA2cSWTOVv5Cac/Tje39CrdsI7atq5qzZ02H25+qoaa9SfqSHLK6ViXDWV+h7XJ631aHftCucXByQZNT/SNCdCjuuDs8BHQrDhownK3TTkB3MN+/K5rYFUrQ3WxIaPJFLNNTAf5nUrvs8xBlhWm0DcAyZczaDpQBzOZkLRzEgiEsGqth+wpsO8ngHMwKvG46q95PqYpWBcDcfw7mOJAWzV7cuVCDO9R6a5PujN1cRbG5X6V1z7dNfOquPZ7oot2z0qzms/7mtwyLn62qHW2gZYvZsm3f1yVfrNDhsLYKo6eo/aN42KECJNcyKoeDXTxOD8oErB2ecTg2totlIqRTuHFafg7KuIUkoVVPkAzARivftCwNMUlKaUisJVW9IHhy6ncnDocra7SUo7h3xnuoZL57sHIRxuX65/k13TEyVwsLnPBkxRSjl0NffluK6B+1zTcPad2LTWDqVUapUh01acTVruId41XiOtdaIrbvd+orziPuRcq8ST4m4689qPOznGaq3jfJyDz1g4OIjEncxjcA71Fs2MTPEjhBDCUNI0J4QQwlCSiIQQQhhKEpEQQghDSSISQghhKElEQgghDCWJSAghhKEkEQkhhDCUJCIhhBCGkkQkhBDCUJKIhBBCGEoSkRBCCEMF5aSn7tUvcU4SaQYyq05eWZuy/uxHCCGEMYJy0lOllEVrPdr1OAZI1Von+1vWn/0IIYQwRtA1zXlNKQ941k/xuV5MdWX92Y8QQgjjBF0iwtmE5qj6YtXEUouy/uxHCCGEQYIxEflambMA36tQVlfWn/0IIYQwSFAOVmhMrqWKkwAiIyPjevXqZXBEwaGiooKwsGD8ntL45FocJNfiILkWB61bt26H1rpLXd8fjInIV60lCh/NbDWUrdV+XEskZwIMHDhQ5/66ksXrdnDBCUf6G3dIycnJYcSIEUaHERTkWhwk1+IguRZOX674mwtO6L6pPvsIxkRkw0ezmmuwQa3LKqXwYz8AFBRrTp36LUUlZcR2PYOjj+jgV+BCCNGcrN+2h5tn2eq9n6CrV1ZNFK7BBXO8nyulTDWVrWk/vpRVQM/lr9Fx1nje+GRBfU5DCCFC3pvzcrhr8+R67yfoEpFLolIqRSkVDyRUufcnjcrDsKsrW922Q3Rspdizdgm/bdnLyg+eYEdRScBOSAghQknh3lKWv/cYd3ywvt77CsamOXdtJt311FplW6IfZQ+7zZdW4dB39HX8sPZ5luWt4J1vf+Xui0+uyykIIURIe3fxavJyfwnIvoK1RmSYxAl3cuyRrSnYV8GvWc9QfKDc6JCEECKoHCivYOHrj7CtqIL+XVrVe39BWSMy0vkndGfexVfyaYvZVHRazWe/bGXcyTKku64qKirYsWMHDoeD8vKmmdQ7duzImjVrjA4jKDTlaxEeHo7JZKJz584y7LqevlzxN9tMx3HmoGV0GjCEPz6ZV6/9SSKqokV4GCePv5f2Od9whNrCO4s+JnHoHbhG4Qk/bdmyBaUUffr0oUWLFk3yOu7Zs4f27dsbHUZQaKrXQmvNgQMH+Pfff9myZQtyv2Ddaa154/sNbO0yHC4cztj4AcytZyKSrwU+jD8llvf1uQCcu/sjlqzfaXBETdfevXvp0aMHLVu2bJJJSIQGpRQtW7akR48e7N271+hwmrS8TYX8tmUXAC0jwrh6WP2TuiQiHzq2aUHJiddTeKAF9t9yee/9WUaH1KRJM4gIFvK7WH9vLshFVRwA4LLBPejcrv59RPK/chhXnHUiV3/biQmfFbPu0/9j/bYio0MSQghD/Vmwjz9evxX9+jhiN33Cf07vE5D9SiI6jL6d29Lr/DsAWLpyE2/OW2hwRMJIDoeD5GRZyko0b2999SM//rqeTYUH6NM3NmCzz0giqsb1V13JKf2iKC6D5e89RuHeUqNDEgZwOBzMmTOHjIwM0tPTa36DECGoqKSMZTMfoahUY+7dgWuu+2/A9i2JqBrDYqLoG38d8ScPYNvR43j/581GhyQMkpSUBEBKSsoh2xwORyNHI0Tjy/ppA333raRFGMSedSlnDajzZNuHkERUDaUUCTfezR8jn2dHx0G888NGSssqjA5LNDKTyURmZibZ2dlYrZUn6MjMzKz1fvypTaWmppKZmRkUNbCq5yyan7LyCjZ89wEvnV3B4tv7cuHEBwkLC9woWElENbjoxO50be8cFbJtTwlfrPjL4IhEY9uwYQN5eXkUFBRgsx2cadhmsxETE4PJZKrVfpKSkmqVWLKzswEYN24cMTGyoLAw3he//UVCyccAfGMaS+LJgf29lERUg5YRYVx3am/n4wO7mTX3K7TWBkclGtPcuXOJi4sjKSmpUtNcRkYG8fHxtd6PyWRi586dNTbl2e12hg4dislkIiEhoa5hB4w0PTZvWmt+/HYuJ4RtYIfuQPtTrqN1y/CAHkMSUS1cdUpvBmydx87XrmHLrPtYapcbXJsTm83GkCFDKr1mt9uJjY31e1/jx49nzpxqVyMJKjabjYKCAqPDEAZamr+TFbOe4PmlJbxVPIqrTx8Y8GPIFD+1ENW2JUNGXMD3H2bw0/oC3nr7bYb/7x6jwxKNIDMzk5ycHMaMGQOA2WwGnP0m7sdudrsdm83GsmXLSE5OxmazYbfbK9WizGYzU6dO9Qx+qMpms2GxWDCZTBQUFBATE0NycjKpqakAWCwWsrKysNls5ObmEhMTg91uJz4+npiYGGw2GxMnTiQ5OdmzLT8/n+TkZBwOB1arFZPJdNjje593VFQUBQUFxMfHk5qaypQpUwgPD+w3YRH83pr5Dl/88g/fRsCd4xOIDsANrFVJIqqlpEtG8fP0GL7JtfPHZy+y5uYkjjlSVnCtiz73fWHYsTc+fYFf5ceNG+czceTn5x/SLGe1Wj3lkpOTsVgsxMbGHjLSrrqmLrPZzOjRoysli/j4ePLy8sjIyGDIkCHY7XZSU1OxWCye98XFxZGXl4fZbPaUd78/Li6OoUOHkpCQgNlsplOnTtUmotTUVEaPHk18fDxWq7VSMnz22WcBZ6JKS0vz2Txpt9vJzs729J+5z9fdzFjTdhE8fv9nN+vnvgDA8BNjmHjxqAY5jjTN1VKfzm05OvE+olsrTjftYPY3PxodkmgEubm5DB48+JDXHQ4HUVGVV6IfN865XuOyZctITHQum5Wfn1/vGEwmE3FxcYAzUWVnZzN69OhKZaKiojyj26Kjoz3l3dtqO+jBXWvy1fdltx9c9DgqKoqEhIRDytlsNtLS0khJSfFst1gslWqS1W0XwSXz61x+WbuVMAUDE1Po07ltgxxHakR+mHhlIqMK5zC21c+8vvZtthSeTs9ObYwOSzQgm83GyJEjD3nd3XTmPWLO/Tg7O9uTgOx2e0BGvnknvZ07G66PMjc315NcvJNtdna2J7mCs4nQ+7nbxIkTycvLq/Sau5nQPTtF1eTs3i6Cy9+79vPZumKOSnqN7pvnc+OV42p+Ux1JIvLDiUeZ+LDHtYzd8TNXhH3DtJwVpF52itFhNTn+No8ZadmyZdx+++2HvB4dHV0pyaSnp2MymRgyZIgnIdntdp+JqLbDvQ8nOTn5kOmG3H059eXuV4KDSclut7Ns2TKysrLYs2ePZ1tGRkal92ZnZx8yqAMq96v5ilFqQ8Hp7SUbKavQ7GnTk47n3ctJvTo12LEkEfnpnNHn8f17gzg9fBXhtjdxnHsSpjYtjQ5LNBCbzeazaS4hIYHs7GzPB6vZbPb0dSQnJ3vuBara72Gz2Q5pVqu6ffbs2Z7mNHeTm91ux2QyeQYluG94dQ9QyMrKOuT97v4kd9JIS0tjzpw5OBwO0tPTfc4SERMTw/jx48nMzCQ/P99zTu79H4476Xo3CXpvcyc475GG7r6uqKioQ5KaMNae4gO8/9PBmWQmntnANVattfy4fgYMGKBrUlFRoe99Ml2nDG+p+0VH6Be/yK3xPU3RwoULA7Kf1atXB2Q/RigsLNTx8fF69+7dPrcnJSX5vc+UlBRdWFhYz8gaR1ZW1iGv7d69W1sslkrnnpWVpQsLC3VeXt4h18RisXge+9oeHx+v8/PzAxx59QL1Oxmov5Fg9KrlN31e/Gn62Dte16OeXajLyyuqLQ/k6np89gbdYAWlVIxSKkUpFe/611SXskops+u1FKVUVnX78TM+hp07nrn5YazfWca3GQ9QfKBpLoEtfLPb7SQmJpKbm+uzH8QtMTHRr+lvHA4H0dHR9W6aawwOh8Nnv82GDRs8tRer1Up6ejqzZ8/GZDJ5Rvylp6djtVqxWq2Vmurc2zMzM7FarZ5BF9I/FFxKyyqwTH+Ir61LUFm3M/GMmIBO5+NTfbJYQ/wAFq/HMUCGv2UBE5DktS0ByKvp2LWpEWmtdWlZub4s4TIN6L5RLfTMJetr9b6mpLnXiDIyMnRGRobWWh+2RuQuV9saTlpaWiBCaxTeNRlv1V2LpkJqRNXL+nmTHtgtUgP6kovP0/tLy2p8D6FUI1JKVfpqpLW2Az6HatRQdgiQ6rXZCpgDVStqER7GeZP+R8+OEewuLifr448pr5Bpf0JJUlJSjTd9usvVlq8+mWAViIEPounRWvNR5lTW/ltMl7ZhjLrpKSJbNPxNzME2WMEMOKq+qJSKcSWa2pa1KqW821RiALTWh5Svq6tPG8icK+9nU5tj+C2iPfNX/cP5xx8ZqN2LJqQpNLUJURs5a7dznSmPwWe3IrfVUK47a1CjHDeoakRAlI/XCnA2tflVVmtt83p9PBDQ+fTbtorgvEuv4kCL9gBkLMp3NwMKIUSTo7Vm/oLPGdN2LTefZmLIjc822ojgYKsRBZyrOc6stfY5ZlYplQQkAXTp0oWcnJxa77s/mogwKKuAX7fsYvrH33JMdGjMxVVUVOTXtTicjh07eu49aarKy8ub/DkESihci+Li4oD8bgfqbyRYrC0oZ9T29yAc3isfTf/Wexvt/BolEbk+7Kubqtiitbbiu/YThY8mOD/KpgGHHfqktc4EMgEGDhyoR4wYUU2Yh/pp3woWLJhPp7wMbC1vYtLjd/v1/mCVk5ODv9fClzVr1tC+ffv6B2SgPXv2NPlzCJRQuBaRkZGcdNJJ9d5PoP5GgkXW1Fe5dP8yitu2Yudx/2XSeSMa7diNkohcH/a1YcNHk5uP/qFalVVKpQCpWmuHUsoUyD4it4lnxPDzU69hzbVz5v7nWZF0I8f37BjowwghRINZsWUXq95/kj6/FzHh4kHcc9vQRj1+UPURVU04rpFxc7yfu0e+1aJsApDtlXwaZBhQ385tOXb8g4QrWLLqL6bPmt0QhxFCiAaTMWsO36/6i7IKCD/tJvo20OSmhxNUicgl0X2TKpCgtfaeVCuNysO5fZZ1JaUsIF8ppZVS2vXeBjHp6ss54/ielGvY9Okz/PFv025DF0I0H+u37WHNnKco13D6cT2ZdOXYRo8h6AYruGo67hFu1irbEmtT1vV6A98KfNAxR3bg6MRUru6dwhUn7OBZy1Ieveacxjq8EELU2Ws5dlp2jeUo0xYGXH6vIeusBWONqElKuu5qup14Bu0iKuj9+xts2rnX6JBEALmXMBAilGwp3Mfc5VtZc+KdhCVlk3T9tYbEIYkoQE7q1YkfjrwBgCvCv+Xdb3KNDUgEjMPhYM6cOWRkZJCeHtDb0YQw1IzFdspcs8IMi+1KXO+GW+qhOpKIAujic8/BUm5m8879bJ7zEH/v2m90SCJA3FP5+Jqmp7qlv4UIVtv3lPDxD6s8z28Z2c+wWCQRBdCwmCiyOYdjpu3lw3nfM33uYqNDEgFgMpnIzMwkOzv7kNm2MzNre2cCftWm3OsNBUMNzJ8ZxkXT8frXP/HvtKvo+9VNDO4azun9OhsWiySiAFJKccU1/yWub0f2lGp+fvN+dhSVGB2WqKcNGzaQl5dHQUEBNtvBmaNsNhsxMTG1nmsuKSmpVonFvajeuHHjZIkE0SB27T/AkhlT2L63gsJCB8nxJ6BUo43vOoQkogAbMaAL/c+7HoCffv6FGV/9ZHBEor7mzp1LXFwcSUlJlZrmMjIy/Jql2mQysXPnzhqb8ux2O0OHDsVkMh2ywqsRpOkx9Lyz8Ddylzo/mwacPZ5zjzN2wuagG77d1CmluCY5hX+WzqVFnzje/203N19wgI5tWhgdWnB51M/ZJx7dFdj3+8Fms/HQQw9Veq3qste1NX78eObMmePX8hFGstlsFBQUGB2GCKCikjI2fjWNjq00praRjL/l4YZf+K4GkogawHnHd+f5a99h7bYiOABv/7CRO+L7Gx2WqIPMzExycnIYM2YM4FxlFJz9Ju7Hbna7HZvNxrJly0hOTsZms2G32yvVosxmM1OnTj1sIrLZbFgsFkwmEwUFBcTExJCcnExqqnN5LYvFQlZWFjabjdzcXGJiYrDb7cTHxxMTE4PNZmPixIkkJyd7tuXn55OcnIzD4cBqtWIymWpMhJmZmURFRVFQUEB8fDypqalMmTKF8PDQmNS3OXtvyTruiVrIMze35c7iZC4+qafRIUkiaghhYYpbRsZy5+xfAXjrhw3894y+tGsll9ujHjWUgLy/lsaNG+czceTn5x/SLGe1Wj3lkpOTsVgsxMbGHjLSrrqmLvdy2t7JIj4+nry8PDIyMhgyZAh2u53U1FQsFovnfXFxceTl5WE2mz3l3e+Pi4tj6NChJCQkYDab6dSpU7WJKDU1ldGjRxMfH4/Vaq2UDJ999lkAT5IFPImwamL2V2ZmJmlpaT6bPO12O9nZ2Z4+Ofc1dDdd1rRdOO0tKaPguzforgr4XffitIRbaBFufA+N8RGEqItO6E6vqDYAtNvwDTMXrTY4IlEXubm5DB48+JDXHQ4HUVGV59wdN845+9SyZctITHROApKfn1/vGEwmE3FxcYAzUWVnZzN6dOVVTaKiojyj26Kjoz3l3dtqO+jBXWvy1fflTjzuMgkJCSQkJGAymeqdhNxxJiQkHHJsm81GWloaKSkpnu0Wi6VS7bS67eKg95b8wfUVnzgfR17B2LheBkfkJF/RG0hEeBiTRsTy5m2jsf68jugdv3D9WV/SVmpFTYrNZmPkyJGHvO5uOvMeMed+nJ2d7UlAdrs9ICPfvJPezp07672/w8nNzfUkAu9km52d7UmuBQUFzJ49m4SEBGJiYjwJ2F0rcdeOUlNTGT9+PLNnz/YkTpvN5kkQdru9Us3MYrF4juFt4sSJ5OXlVXrN3fTonvGiasJ3bxcH7S0p46v/uwNL/p9cf2ZvjvvPNUFRGwKpETWoy8096dD3RACWLf6WNxeuNDgi4a9ly5b5XLsmOjraU0MA5z1CmZmZ2Gw2T0Ky2+2VyrjVd2lxd/+TN3dfTn25+5XAmZTMZjN2u51ly5Z5kkZMTAzjx49n9OjRlQZsOBwOT4IGDjl3d01lyJAhnsfevJOgW3Z2NkOGDDkkTu/akK/zltrQod5ZtJpli7/lozVlvL0rLmhqQyCJqEG1jAhj7B1P079LK/7eXY7l5bspKikzOizhB5vN5rNpLiEhoVIyMJvNntpDcnIy2dnZ2Gw2n81MVZvVqm6fPXs2WVlZWK1WbDYbVquV2bNne5reYmJiPDe8Wq1W0tPTycrKOuT9NpuN7OxscnNzycjIwOFwkJmZicPhOOz9TO4kk5mZicVi8ZyHe//gTDgpKSnk5+eTmprKnDlzPNfAYrEQHx+PzWZj6NChmM3mStfBZDJhMpk826vjTuTezYze29z/eidDu91OYmKizAtYxd6SMua/fA9/7S6nf5dWXH77VFpGBNHHv9Zaflw/AwYM0IFWWlaux14xXg/rGa4/uu4I/Zp1ZcCP0RAWLlwYkP2sXr06IPsxQmFhoY6Pj9e7d+/2uT0pKcnvfaakpOjCwsJ6RtY4srKyDnlt8+bNlV7PyMjQeXl5nucJCQlaa63T0tI8r8fHx2uttc7Pz9cpKSlaa+d1yM/P95SxWCyVrmdWVpYuLCzUeXl5h1xni8Xieexre3x8vM7Pzz/seQXqdzJQfyON4bWc9Xr4wC4a0GPHj9MlB8oDun8gV9fjs1c6LBpYi/Awxk1+muP7reTY8D+Z+t3r7DktnfaRcl9RsHKPSktOTvbZZ+GWmJh42KYhXxwOB9HR0fVummsMDofDZx/LL7/8QkFBgWf2h5iYmErNYEOHDvXU4FJSUnA4HJX6h9yPo6OjPTUlu91ORkaGZ8CFewi8e5Tf6NGjSU9P9xzHu6nOvT0zM9PTZzR69GjpH/Kyr7SMzMV2dl7yBqM2fcLYOx4MrtoQSI3I+6chakRaO2tFDzz5pNaPdNB/P9xbT1uwokGOE0jNvUaUkZGhMzIytNb6sDUid7na1nDS0tICEVqj8K51eKvuWuTn5+v8/HxdWFjoqfkEo+ZWI5qes173Tv1c9079XA+f+k3Aa0NaS42oSWgRHsaJo69m9bwPOUZtZMvXL7D7jNfoILWioFXbmQ+SkpJqPQWOr5m7g1VdBj44HA4KCgo8w6mF8dy1IbdbRvYLvtoQMny70VxmPoqUzy4hb+ZUft/5Fp3PSuaui6vvrBVNQ1NoamsMMlIt+MxctJp2Wf+h/YmXcWDQZSTEGT+Lgi/BlxpDVER4GMMvnsj2kpZsK6pg4at3s2v/AaPDEkKEqH2lZXz9cgqLV2xly9eZ3HRm76CsDYEkokZ1qbkng85x3vz30/dLeH2+rYZ3CCFE3byTsxrbd98AcMLIixh/cl+DIzq8oEtESqkYpVSKUire9a+pvmWVUhkNFa8/IsLDuPL2xzn2yNb06daBTxb9zK59UisSQgTWnuIDfP/mQ/zpKCO2c0suC7b7hqoIxsgytNbpWmsrkA1U1+tZY1mllBkImjn3LznpKHrePJNtie+xpVUMb3x/6J33QghRH2/nrOH5gctZeH0bTr3kGsYP83/JksYUVIlIKVVp8L/W2g6Mq2fZGMARoBDrLTxMcfOFp3qev7lkI4V7Sw2MSAgRSgr2llK6NINuykFUr/5cNunRoK4NQZAlIsCMj6RRNenUtqxSKkFrnR3A+ALiwhO6069rOwCKiot5bcGvBkckqrN8+fJaLfHtD5vNRlxcHJmZmXXeR1xcnGfaH6NiaIj91fe8mru3vvmVCXwKwAftrueSIFhvqCbBloiifLxWAJj8LetKSEHZ7hUeprhr9AB6b51P5DtX8N3/JfHPrmKjwxKH8dZbb5GREdhuRve6Qf6oer/SjBkz6jXRaV1iaIj9Bfq8mrN/dhVjm3EnP64v5MfygZw1ZjzhBq++WhuhfB+RuTa1IaVUEq4+pC5dupCTk9PQcQHQWms6R2oW/1vM1sJfeeCVd7n2tOBZxbWoqCgg16Jjx47s2bOn/gEZqHfv3tjtdn777Tf69g3cyKN27dpRXFxcq+uzYcMGfv31Vy699FLPa/3796/3tfUnBoDy8vJqy/q7v4Y6r+oUFxcH5Hc7UH8jgfTe0nw+W/AD75dqxk06g5t2rCUnZ53RYdWoURKR68O+ut4yi2vAga/aTxS++3gOW1YpFQ/Uqm6vtc4EMgEGDhyoR4wYUZu3BUR4j+PYmvM+P60vYPNnz9A7+Rf6dm7baMevTk5ODoG4FmvWrKF9+/b1D8ggVquVyy67jO+//55Zs2YFdMaAVq1aERkZWavrM23aNEaPHh3wa+lPDAB79uyptqy/+2uo86pOZGSkz6U9/BWov5FA2bxzH39OmURRqebk2E7cMvk+zhzQxeiwaqVRmua01pla69RqftxJw4aPJjfXQISqaio7TimV5EqCJtfjoJoJ8Yz+nTkm4W7CFBxRsok3Ps8xOqRGo5Sq9NPY22vLZrPRt29fz9IOVbe5+0SsViuZmZmHLD/gnsQzMzOz2n6m7OxsYmNjPe93OBzExcWRnp6O1WrFbrdjsVg8yzj46o+x2Wye8lar1dPkVdsYqnIvG2G1WsnOzva8170v9zn7WnOp6vVxn2OnTp08y2f4c16+jleb69+cvPBxDt/l/g7A0ZdP5oz+nQ2OyA/1maiuIX6APK/HMTiHaHs/N9WmbJV96tocu6EmPa1O3qYCnXnrWVo/0kHPeeBCvWKLo9Fj8KWhJz0FKv009vbaysjI8Ez0CVRa8kBr53IG3ssQmM3mSksQxMTEeCZFTUpKqjSZaFpammdiVfexvCcL9V5uISkp6ZBlGbzfX1hYqM1mc6W43Nv8iaHq/r3PNysrS//666+eZR28z/lw+/OOQ2vnEg3e+6zpvPLz86s9Xk3X35dQnPR07T+7dZ/7Pten3ZCqzxt5iv55w85GPT71nPQ02AYrACS6b1IFErTW3l9x0qg8RLu6siilTEqpFNfjtGCrEQGYe3Vi4+A7OaDDGRv+HbM+X2B0SI2i6i9iY2+vDfeS32+99ZZnmYHZs2dXKhMdHV1p4baoqKhKne95eXmYTCbsdjsFBQXVTpCalJTkqXXZ7Xa/5m6bM2dOpeUR0tLSPBO3+hODt4SEBBITEz3LMCQkJDB37txDFvZzL9/gS3R0dK3PwZfs7Oxqj1fT9W8unluwFq1hS7czOOLalxjax9dYruAVdIMVtLNpzd1+YK2yLbG2ZV3bHa7tgR17G2DXXziK2a+M4Jrwbzht83R+so/ilJj6/QGL+rPb7aSlpXn6RaKiopg4caJf/URTp04lOjqahISEWq2Rk5CQQGZmJlFRUSQkJBw2Ln/W2/E3BreoqCjy8/M9q74mJibSs2dPWrZsWet9+MPXee3cubNBjhVKfv3TwfxV/3qe333OQAOjqZtgrBE1OwO6tSf/6FuYtUoz+bWFTJv2Yp2/xYvA8LUwnDsxeC8RXh1330xKSopn0Tb3MgmHk5ycTFpa2iEzens/99UnEx8fT25u7iHHr0sMblOnTgWcw7Ldyfc///nPIedfUFBQ7XBr7xpK1RhrOq/k5GS/j9fcvPT6mxy9ZgboMi44/kiO69HR6JD8V592vVD7MaKPyG3zzr363OGDNKDj+nTUllX/GBaL1s17YTyLxaLNZrOOj4/XhYWFnj6irKwsbTKZtNls1llZWTovL89TLi8vz7M9KSlJFxYW6sLCQp2QkKCzsrK0xWLxLIedl5dX6b1V+zR8LUGen5/v6U9xL7Fd9f0Wi0WnpaVpi8XiWWq7rjFo7eyrycrKqvTe3bt3a4vFojMyMjzHc7/X1/4KCws9/VIWi0XHx8frhIQET59Vbc+rpuP5uv6HE0p9REv+2K5PHdBZA/q8s8z6j3/3GBIH9ewjMvzDP5h+jExEWmv98Dtf646RSgN6/MSbdVl5hWGxNOdEVFV1q5I2hKqd98Gksa9FQwiVRFReXqGvuuUeDeg2LdD3vDLbsFjqm4ikaS6I3DZ2BKedOoSIMHDYl/PZr1uNDkk0kuTkZOx2O1ar9bB9Q0J4m/frVn6fNx2AM4cN5s6rLzI4oroLusEKzVnndq0YcduLbBuwiN9Nx/PcgnWcf/yRtIoINzo00cASExOx2Wx+DSYQzVdJWTnLvnqHOZfAo9+35djbXqS7qbXRYdWZJKIgkzxmCO+v2AP7DrClcD/vLt3EjWfIh1Ook8534Y9ZP+QzoXgmMVFhHHPJjZVm9G+KpGkuyHSIbMHtZx+cc+4l6xoc+2SZCCGE0679B/hnYQYxYf9grziCI0Ym07F1C6PDqhdJREHo6lN608fUkqOXP8u+V8fy0ieLjQ5JCBEkXv18KS1+fpPScs0bkddx9fDgXvSuNiQRBaGWEWGkjDmWbWvyyN9Ryg+vTmbzzn1Gh1VnFRUVRocgBND0fxf/cuxn8ct3kPKlg4s+ieDkMdeHRB+yJKIgNeaE7gxKmAxAzk+reObtOcYGVEdt27Zl69atlJaWOu8XEMIAWmtKS0vZunUrbdsGxwz3dfH8rHks/nE5ANEjb+SiE3sYG1CAyGCFIKWUYvLtd7Jp/uvkrfub1j++zC+bL+GkXp2MDs0vPXv2ZMeOHWzatImysjKjw6mT4uJiIiMjjQ4jKDTlaxEREUHHjh3p3LkJzUrtZfVfu1n+9v3sOwDDB3bm1skphDWBRe9qQxJREIvr3YmT/vs07/x5B307rCdl7jwG33ptvZY1aGxhYWF07dqVrl27Gh1KneXk5ARk/ZpQINfCGFprnvhiNduH3c2o8nRiL53M8NimmVB9kaa5IHfntZfzTdsxAFy67TUWrPrH4IiEEI3NumYbP+TvZFfbXmw8bxr3Jk8wOqSAkkQU5Pp0bstO8204dFuGh69m4WczKS1r2h2uQojaKy2r4Kkv13ieX3VyL/p3a7qrHvsiiagJSDonjmnll/Hk4hI+f/kR3s5ZZXRIQohG8tai1UR+fT+RpYW0j4zgztEDjA4p4CQRNQGd2raky8hJvPmbZtU/JcxLn8SOohKjwxJCNLDCvaV8/uwtfP3NUtp+fAt3nN2fqLYNsx6UkSQRNRE3nHUscRdeA8B3i3/guQ/nGxyREKKh/V/Wt3yf8x0Afc+8lOtO7WNsQA1EElET0TIijJvue4qTYzsRHqb4aeGXrP5rt9FhCSEayPpte8ibcSeOYk1cn47ceO9UWkaE5ke2DN9uQkYO7MrxNz5H0T8V2COP4PHPV/HBxGFNaji3EKJ23vh0PjPP+Id01YptZ97HuccdaXRIDSY002uIUkpx34QEits4fyF/tBcwX4ZzCxFyFq7dxug/XyY6soITR47hnlsnhfQXzqBLREqpGKVUilIq3vWvqa5llVIJSqkk1/aQmGe/X9f2XDust+f51Ll5FB8oNzAiIUQglZSVs+CTdxgVvpzdujX5g+5gUPeORofVoIIuEQEZWut0rbUVyAbS6lJWKZUAxGitMwF7DftpUibH96dL+H4G5NzNlv+7nBnWX40OSQgRIK9bfmP351PZua+C6SqRSU18raHaCKpEpJSqtAKc1toOjKtj2TStdbp7m9Y6LsDhGsbUpiW3nHMcG+0bse8s5etnbuLvXfuNDksIUU//7Crmy2eSeXVJIed+WE7Xs2+jc7tWRofV4IIqEQFmwFH1xapJp6aySinPNtfjkHPN6QMwX3ojAAu/W8bUN5rm7NxCiIOmvplNzvfLAIg9bwLXDO9ncESNI9gSUZSP1woAk59lY4ACV/Oc3d2PFKggg0FEeBiT73uM047pRuc2il6/vcIP63cYHZYQoo5+su9kzYePs+8AnHZMVyanPkZEeLB9RDeMUB2+HQXEa61HAyilMoENwCFrKCilkoAkgC5dupCTk9OIYdZf/ysf5eN9D9G11Tomv/86+04/hYgATA1fVFTU5K5FQ5FrcZBci4MCeS3KKzSPLi3m7zOe5Ly2z9F5+HhK/lxBzp8B2X3Qa5RE5Pqwr249W4trwIGv2k8UPprgaihbANjcL2qtHUopk1IqxtWXhNe2TCATYODAgXrEiBHVn0yQOcY8jPee/ZW7eJ/byt5mYdgl3DjimHrvNycnh6Z2LRqKXIuD5FocFMhrMXPpRv7cswoi2rDxtEd54+6z6G5qHZB9NwWNkohcH/a1YcNHk1vV5FFT2cOMt3fgTFAhpVuHSDqOmkz+twuJDfubWfOe5l9zBt06NM3Fy4RobnYWlfDcgnWe57eO6teskhAEWR9R1YTjGqQwx/u5+16h6sq6thW4y7q22bXWjgYM3zDXndGf6a0nMmHufp555V2eeiPb6JCEELU0deY8TNn/IWr3OnpHt+HGM/oaHVKjq1UiUkp1UEpNVErNUUrN9/qZrZS6USnVIYAxJXoNLkjQWid7bUuj8hDt6somAlNczYLJwNkBjDGotAgPY+y4G1hX0omiUk3ejLv5IV8GLggR7H6y7+TnV29n8cq/aJPzNI9eNIhWEeFGh9XoamyaU0rdi7MJbLbWeoaP7ScB9yuldmitn61vQK7aTLrrqbXKtkQ/yjqA1PrG01ScGhuNeeJz/JJ6LUvWbOPFqQ8xZPq0kJ0kUYimrrSsguefepDvVv1N6wgYfMOTjDy6q9FhGaLaTylXEsrUWk/RWi/3VUZr/YvW+j5ghlLqngaIUdTS/RMSGHnGycR2bsk/pS2Z8Z2vrjUhRDB4a9Eain+aBcCoM0/mgeSrDY7IONXWiLTWz9R2R1rrXUC9a0Si7o7oGMnFD7zJ71+tYVN4a1765g8uPOFIeke3NTo0IYSXPwv2UbroeT4ZC8/Fdsd02xsc2bF5DVDwJu02IWbCWcdwdK9uAJSUVfDgpyvRWhsclRDCTWvNtI/mk6Q+JUwptp/+IP8ZcazRYRnKr0SklLpXKfWH63FfpdSohglL1FVEeBhTLzsB9z2tW759m09tm40NSgjh8fWKvxmz6VlaqTKyys/kynFXNJsZFA7H3/uI7Frr/gBa6w1KKZRSY7XWHzdAbKKOju/ZkeuH92HxIxfy7S8b6bh3E6Pen0vHNi2MDk2IZq2opIwZT9xO4bKfePr8KNaNuJfEXodM+NLs+JuGba5h2zeCMxk1QEwiAO4+ZyAdYp3zvS6a/zVTZ31pcERCiGc+/Ibvv/yEnI3lTN1i5tYLQn+Jh9rwNxEl4ZwOp59Sar1Sahm+Jx8VBmvXKoLkh19maIwJR7Hmx5dvIXdjyE0sIUSTkbepkMX/N5HteysY3Ks91z2UIa0ULv4molyt9Uda6/u01v1wJqbcBohLBMB5x3cn7sanOaN3BC8OdzBnzruUlMlqrkI0tpKycu6b/TP7iopoFQ5DbnyKy+J6GR1W0KhL05z3TApxOJdcEEHqoUnXcdl1EzjxiHBuLXqZDOtKo0MSotl5dWE+fxSU8c/l7zBs4mM8cusEDjMnZrPkVyLSWm/QWr+utd7teqkTvmfGFkGiu6k1UaPvZlVFb3qFbaf9kqdY/dfumt8ohAiIdf/u4dWc9QAoFcYN11xLz05tDI4quBw2ESmlOtY0h5zW+hmt9bde7+kQ4HnnRABce1p/3u58DwcqFPq3eTz59JOUlVcYHZYQIa+8QvNwxhzKD5QAYO5l4tpT+xgbVBA6bCJyzZQwvrb3CimlzgbGedWWRJAIC1PcdOVYJiw/lv9+tp8f33+BjAW/Gh2WECHvrUWrWT7tViLfvYoj9qzh6ctPIDwAC1eGmmqb5lyTnBYqpaYrpe5RSo1SSvVx1Xz6uJ7fq5R6zVlcv944YQt/xXZpx7BJ04iJbol9Zynznrye/O1FRoclRMj6s2Afnzx+A2u3FbO3pIxx55zBgG7tjQ4rKNXYR+Sa1PQm4COcq6zehHPG65sAM5CttZ7k3UQnglPyOSdyylV3ooAflq3kgbe/pqJCpv8RItAqKjQPPTcN63fOQcXDEpK49YKhBkcVvGo9s4Lr5tVDloEQTUdEeBgP3nc/jvW57Og1ih93tmbm0o3ccFrzW4hLiIY0a+l6JpS8TfHREWxv2ZP7HnxClmSpRk3LQJwtSzuElmO7dyB+8jS2mU4EIO3rtWzcsdfgqIQIHRt37KVwfjpnttvCs5f34vxHP+QkmcanWjWlaBNew7OVUoMbMBbRSG4Z1Y8B3doBsP9AOfd8mEe5NNEJUW8VFZppH3zMJPURAC+1m8yt5w8xOKrgV1MiigKGuG5iHQzIFQ0BrSLCeX7cYCJUOcfY0vh96hgyFiw3OiwhmrzXv10Jn6VQVlbGzPJzuO6q65rl0t/+qs2ouTScN67eD6Qrpf5QSs32Sk6iCTquR0duPvVI1ub9yKq/9jH38av5/R8ZeS9EXeVvL+KTx67hWes/XJhVjuO0BzmuR0ejw2oSajNqboPrxtVxQCrOWlEmruTkmvzUMyO3aDpuuyCOU6+9hzAFlqVrePSJxyktkxtdhfBXeYXmkafTsSz5DQV0jb+Jm+KPNzqsJsPfKX5maK13aa2/cScn1+Sn44ANgUhGSqkYpVSKUire9a+pLmVd25KUUgk17ae5ahEexmP3T+Gc046nQysYuW0Wr8mNrkL47bWc9az/bBrlGuJPGchDDz4io+T8EJAr5ZqFITkQ+wIytNbpWmsrkI2zabAuZRO01pla62ytdTowJUDxhZR+Xdtz2SOz+Oimo0k+dh9dlv6PXzYXGh2WEE3Gr386eMH6B7sufJ5zTj2WC+6fybHdZaYzfwQyZafWd2YFpVSlmby11nacta26lB1d5S2m+sQWym4cdRxf9XuIEh3BVeHfMPv9Nykpk1F0QtRkX2kZk2cvp6xCs79VZ9pdO4NbLpAxXf4KWCIK0GqtZnzM5l016dS2rFLKopQyKaXigawAxBeSwsIUd1x1KdO0M4+Ptj/Dhz+tNzgqIYLf45/Y2OC6D69dqwheGD+YiHBpkvNXsF0xX6u9FuC7NlNT2URXmQ2A2dV8Jw7jqKg29Lr4Pp5c0Y0xb/zLmrfv48tftxodlhBBK+/vUhY/fAGxC26lVamDxy4exFFRsrxDXdR6ip8maAjOUX4xQIZSCldfUSVKqSScK83SpUsXcnJyGjXIYBKtNb8PupW2X07mp/UFvJp6FfvveITo1sH2faVxFRUVNevfC29yLZwcJRXkvJvG0nU76NhqJxdesJ2o3X+QkyMtCXXRKInI9WEfW00Ri6vG4qv2E4XvxfcOW9bVPGd2Jx6llBXIU0plaq0r7UtrnYlzODoDBw7UI0aMqPmEQph52OnsW/8jH3/wAVbrYjqduZm377u+WU9dn5OTQ3P/vXCTawFaayY+/jLfWJzzPI+6eCzTptxCxzYtDI6s6WqUROT6sK8NGz6a3FwDEWpd1pX4rFVey+TwSU24dGzdgilPvcKeDcvhyOP4bncXpi/K55aR/YwOTYig8NbitZy7dRqLO4TRvfsR3DN1uiShegqqNpeqCcdVs5nj/dx9P1ANZa1AfE37F74N6RPFsf99gXUDrgfgecs6GdItBLD8Twdh1kdI7P4P8yb25pwHZjE8trPRYTV5QZWIXBLdN6nivBfI+/6kNCoP0fZZ1pVw7K5tCa4aUkajnUEIuCimBUN6O2cMLq/Q3PZ+Hrv2HzA4KiGMs2v/Aea8+xo3hH/NAR3Oi+3v5J6EM40OKyQE3WAFVxJxDyqwVtmW6EfZ7IaKsTkID1O8cMVgxrz4HZ02LmDHO69yf9g0Xrn3PyjVfPuLRPOktebBae9xUl4aHA8vqKsZedIgmT0hQOQqisPq2akNaZefQPlP77Lmn2IWvXgbb3y7yuiwhGh07yxew6IXbuO/H+/ijp+7cdzYKXRpIx+fgSJXUlTr/OOPJP6+d+nRMYJVf+1j9pSx0l8kmpVVf+3igymJrNi6lyM7hNNh3AuMOaG70WGFFElEokb3Xx3PiBvupUUYHKM2kTXzFRz7So0OS4gGV1RSRuq02SxetppwBWdeO5kHrj7X6LBCTtD1EYng0yoinKkPP0h3/RfpnT5h94FpTJ11PE/+9xLCmvH9RSK0aa25N+tXVpb3IO7a++i03cZTjzxKZAtZ6C7QpEYkauWoqDaMTnqaL8tPpoPazzV/PswbOWuMDkuIBpO52M5XK/8BYEvX07nxf28S06WdwVGFJklEotZGDzqCNUOfYmNFN7qXbGTpM1fww/odRoclRMD9sH4HGVlfep7fMLwPlwzuYWBEoU0SkfDL7RfE8VLU/Zz59j7e+noFaXdexZbCfUaHJUTA/OXYT9pd17Au83YGrn2TIb07cf/5xxgdVkiTRCT80iI8jNQJVxB7+gUAfPWlhXsffoziA+UGRyZE/ZWUlXPv41P5+sv57D0ArcMOMO1qs9wv1MDk6gq/desQyUPPv8Xok/tzoAI2W95iyscr0FoW0xNN25PvW8h54ylKy+HsuBjuf/ZtunWINDqskCeJSNRJXO8ornl2LuedfSr/XjiNT37ZyltLNhodlhB19uGP+Yxe9xiPntmCuN7tuCJtLif39bXsmQg0SUSizq49/WhOnDSdsgjnYmBPfrmGpfk7DY5KCP/9aN9J6Rf3cUbEai41d+HcR7P476hBRofVbEgiEnWmlOKxSwYx+CgTAOXlB3g8PZ3NO2Xwgmg6Nu/ch+WdJ7kufAElOoKnOzzIA1fFy5yKjUgSkaiXVhHhTL8mji6tFT3nJWF9PZ17p6Swu1hm6hbBb0/xAe65P5U5019ifUEFT4VP4u7/XkublnKvf2OSRCTq7YiOkUy/4VTamqIpq4D5777KnWmvUVZeYXRoQhxWeYXmnudexzJzGiu3lTN5eV8uueFuuptaGx1asyOJSAREXO9OJL/0BacO6MyuYs1vb6Ty9MdLZCSdCFpPz87hmxcm4yjWDOsfxfXPfoq5Vyejw2qWJBGJgBk7pC9j0+ZxwXEmssZGMHrFPcz8/g+jwxLiELOXbebtvG30OKIz/bu04qL/fUTiyTFGh9VsSUOoCKi7LzmFXbs/JDI/mVPU72xecA8LO89g5DHdjA5NCAAWrt3G/Z+spDyiHRvPn855PcpITTzL6LCaNakRiYBSSvHAlfE8E/Uo+3VLEsMXYXnlDlZu3WV0aEKwcusubn4vl/IKZ5PxcT1MPJd8EeEyi7yhJBGJgItsEc59E67kfy0n8+qyUtJem82Dt13Dpp17jQ5NNGN/FuzjwduupfPHE2hRVkQPU2veumEo7VpJw5DRJBGJBtGlfStuuPF2vqo4hQoNn3/6OXfffQc7ikqMDk00Q459pdydcjdfzJ3H96v/5ehNWbwzYShdZfqeoBB0iUgpFaOUSlFKxbv+NdVQ3qyUyqvvfkTgDejWnodf+5hzhw+iXIPlwze56ZW57C0pMzo00YwUHyhn8uPP8fW7GZRVwDnDjubBtAz6dW1vdGjCJRjrpBla69EASik7kAYk+yqolIoHCgBzffYjGs7QPlFMnr6AshtOR8UMJ9fRlptn2Xj9+iG0CA+670EixBwor+CxmV9yy/5p/NxBEd21KzdPm88pMdFGhya8BNUngVKq0vhJrbUdGHe48lprq9baVt/9iIZ13vHduWHaQv6IvRKAReu2k/rRb1RUyD1GouFUVGie+GAhN226i6Gd9vDCDSdx9fNfcYm5l9GhiSqCKhHhrNk4qr5YNbE04n5EgFwzrDe3n93f8/z7b77gkZlfyQ2vokForXkqezFXrbud3mHb+K2iLyvOnE7yOYONDk34EGxNc77mXC8ATAbtRwTQnfH92b6nmJwvP2L9h0/w8bxWRLb6gvuvGCETTIqAeubj75lz1/l0OqGCsCF9mXvcyzx4YZz8ngWpYEtEjU4plQQkAXTp0oWcnBxjAwoSRUVFDXItRnfSbOnRhb9ahrH67/1kpV7M9h1vcOnxXQN+rEBpqGvRFDWFa/HV2kKsU29gxdZ9pO2NwHbGPVzdpYxFixYF9DhN4Vo0FY2SiFwf9rHVFLFora34rrVE4aOZrQa13o/WOhPIBBg4cKAeMWKEn4cKTTk5OTTUtTjjzLO4vU1r5j6VzK9/FnHm56msH/QFN446tkGOV18NeS2ammC/Fu8u3cgPL12BbdNuurYLY8w9L/FKyg0NMjAm2K9FU9Ioicj1YV8bNnw0q7kGG/gjUPsRDaBFeBgv3H0DB0r209H6AM8O28m3C2/mvVYzuOa0/jXvQAgfZi7dyMNzV9Hv+DF02/wu59z2BC+mJMnozCYgqP6HqiYK1+CCOd7Pa3M/UE37EcZrFRHOy/dNQo2bTiHtOTv8Fzp/PYn3lsgkqcJ/7/zgTEIA6/skcNrj85j28N1Etgg3ODJRG0GViFwS3TeiAglaa+97f9LwGobtvlnV9TjN9Z7a7EcEgdYtw3l0YiJPRj3JLt2Gsyp+ZsGU0WQuWG50aKIJefPblUx7Z5bn+Um9TLw9KZ72kS0MjEr4I+gGK7hqM+mup9Yq2xKrPLe6yqRTRXX7EcGjXasIHkm+msczFBtfn8RHKzfx+6SRFL9k4fYLhhgdnghyb367ktcnjST/zwKOviKVbkMu4J0JJ0sSamKCsUYkmpkOkS149KaraXPBw3RuE8Yyu4OZt57Nc3N/lPuMxGG99vUvzLhpBEvX7UCjiYluLUmoiZJEJIJCu1YRTH/8Hi68K52u7cLYs6+Et75dwTPz10oyEpVorXl+Xi5bpiXw4x876dRaccGtj/Hq/+6TJNREBV3TnGi+2rSM4LWHJ3NLi1Ys2x7O7tbdeTUnn8J9B/jfJYOIkNFPzV55hWZq9mIuXnkHx5u3s39XBwpOf4CXH75TklATJn/ZIqhEtghn2pRJDI072D/0wc+bufXVzyg+UG5gZMJoJWXlPDxzPlesvIkTwjawWXej5XXvMv2xuyUJNXGSiETQiWwRzvRr47jspB4AHLPsCT69P4Hkh9LZtf+AwdEJI+wtKePmJ17ht2cSOLJsK2sqjuLNAdP5338ulCHaIUCa5kRQahEexnOJJ9K5teLjD1fzb1E5c194gGLHv7yY9jRHdJQFzZqLf3YVc+eDj2B581l27qtgUsduHHvPDB65dBhhssR3SJAakQhaYWGKBy4+kaSMRZx2TFd2lWgWvvsSzz33JKv/2m10eKIRrPprF8l33cMnr6Wzc18FJ8d2YnDqp6ReJkkolEgiEkEv+dzBTH5nKecMjeWDy1vzTNhLWKbfg2XVP0aHJhrQt7//S+L0paztOIy+nSMZZe7DHW//yF2XDJNZtEOMJCLRJCQMjeH+d39ged8bqNCKO8Jms/fDCWR+s0qGd4egmUs3cuM7uewrLae4ZSfa3vAmj77/E1edPsDo0EQDkEQkmoyzBnblopuf5v5W91GkI+m8dREZE4Zx72sfUVImI+pCQUlZOVPemMcHj1xPeUUFAD1MrXn/zgs5Y2DwLhUi6kcSkWhSBnRrT8odd/FI9HPcMr8c2+YiPphyBTfeP5V/dxcbHZ6oh7937WdCyuO8e/dYvrYs5tjlz3HiUSY+veU0BnRrb3R4ogHJqDnR5ES1bclTN19JWMvW8OT1/LZlLx+/8BB/0oWnki5neL/ORoco/PSjfSf/d/+NzM3+lJJyOL5HWwaMuZGnJw6jdUsZnh3qJBGJJqlVRDjpyWOJ7RvDnPsupUXbjqwN68k1b/zEXaMHcPOIfjKqqgnQWvP2d+toabmfxzpb+CocTj+xL+PS5jLx7ONkUEIzIYlINFlKKW469yRO6JfHbbPyYF8ZFRqeXbCOZfZ/eeHKk+nUtqXRYYrDKNxbSvqH8xm36VFOCl9PSXQkN95yHeNu/R/DYqKNDk80IukjEk3e8NjOfDn5LE7u41yUt8PezSyZcjb/nXANP+TvMDg64csP63cw4YYr6P7p1ZwUtp6tOpr7Oz3DvY88J0moGZIakQgJ3TpE8v7EU3hmwVoWpT3K/B2l5M/JomDND5x930zuTTxLpoIJAgfKK3j2o8XMf+IaFq3Yylfh0K73cP4d+QxPXXIKrSLk/6g5khqRCBkR4WFMGXMMd73yGRedP4rWEbBoxVZK307k4Rems+ZvmY3BSOu37SHhlcW8e8cYFq3YSusIOPf8czhqUjb3XT5cklAzJolIhJxzjjuSGR98wVUPvsjl5mgeGXaAp4se4IdXbyLzm9WUlVcYHWKzcqC8gle+/YPzX/yeX//ey1Enn8PR3SK58oH/Y8Z7c4kfdKTRIQqDSSISIalbh0hmPHwbl760lOkkUIHiv+Ff0H/uxVx7270s/9NhdIjNwsqtuxj/2Os8u2Adpa4vAPbjkrnjvWXMePgOunaQyWuFJCIRwpRSXHNaf8bc9hJTTM+wvvwI0r/ayIevPc9dF5/II+/MZ0+xLCvREPaXlvPU7BwmX2Lm8yeT6bntewBO6NmRebedzk3xx8nweuERdIlIKRWjlEpRSsW7/jXVUN6slMo7zOsprp+smvYjQldsl3Y8dfsErGdko7ofT8twWLh8M6/efD6X3/8SX674W+arCxCtNZ8t38qEG6/jlaR4vsmzo4GO22w8cP4xfDxpOEcf0cHoMEWQCcZRcxla69EASik7kAYk+yqolIoHCgBzlddNwBCtdbrreQLwDRDXcGGLYNYiPIybzjme8+OW8PgrM7C99RBoWNtiADfPsnFy3ygevOAYTuhpMjrUJmvl1l3M/HguCdtf4ZqylXy4u5xjj2zN0Osf4tG7b6dP57ZGhyiCVFAlIqVUjPdzrbVdKTWOwyQirbXV9b6qm4YAqUCm67kVyFJKmbTWjkDGLJqWXtFtmPHw7cy7+HLSs3IoqHA2Cvy8oYCkx17k1EGx3H39WHp2amNwpE3Hv7uLeemjbznensnTYTmEhWm294vi6quHceEtTzPulFhphhPVCramOTPgqPpi1QRVE1eCSvR6Kcb1+iH7Fs2PUoqLT+rJFw+N5z+n9SEiTKF1Bfu+eobp917LxIuH8eScRewsKjE61KC2fU8Jj75n5YYLT+WVSRdy5J8Wygnj9fILeCfuY159/QOuOFWmWhI1C6oaERDl47UCwOTvjrTWNq+n44H0OsYkQlTH1i145KJBXHdqH9I++p6VHdpRvHUvX3+/gp/zRsHNV9Hy9NvoHynDvb0V7C0lY3E+3894gCVffcquYmf/2iu/d+Kr4S8y8bJzie3SzuAoRVMSbIko4Fz9RWZ3v5OP7UlAEkCXLl3IyclpvOCCWFFRUbO6FlcMbMXaJ2bRbf5n5H/1Jj0ji3mg3Wfs/sXKzPJzuH19AWfFRhHdOtgaERrPtn0VfLl+Hz9aLBSXQ+92vdhVrInr04He8ddxxjmXcEKXCP5clcufRgfbCJrb30hDUo0xWsj1YR9bTRGL1trqGlSQ7J00lFKFQJzW2l7N/rXW2mf9XymVAaTWpllu4MCBeu3atTUVaxZycnIYMWKE0WE0uooKzdcr/sL6+YckFGczPHw1AB+ugTc3HUn/S+7kqsTLGR4b3SxmhtZa8/OGAt7+MIsNX03n9+GPgTr4/fWU0jxuvDGZ+GO6Novr4a25/o34opTK01oPqev7G6VGpLXOrLkUADZ8NM9Vl4Sqo5RKwZWEZKCCqI2wMMX5J/ZgzAl38e3v15A6fy5jdsxk6qKl/PbvWiw/3cSy1++j5yX3cPHlVzE2ricdIlsYHXbAOfaVMm/5FubPeol8y3v89MdONBB/xBz+iL2KAd3aMTl+AOcNOl/6gES9BVXTnGuUnOe5a5DCnCrPC2qTUFy1q2yvsvFAdiDjFaFLKcXZx3Rj1NETmfZRLHHXLyHq6xn8uHoLuXYHel84j85bzdNf/87ZR3fjohO7M2JglyY9seqB8goWr9vO4h9/oov9Ey5R37Fo6WZ+/OMArcLh1ON6EdVnEO9MOJkz+3dudjUg0XCCKhG5JLpqMjacfTveQ7fTAAuuYdmu+4jMrsdpHGziiwGyXK+732tHEpHwk1KK4zqHc2vaw/xx153M/Oo7fvniLVaajgeg+EAF1rw1LH1gBD0GDmbgxbdwQfxITo2JbhJJqaSsnB/yd7IgZwlrvpjOYLWOx04qAFfol5q7sLV9N0685hH+e/Eotq7J46wBXYwNWoScoEtErmY49wg3a5VtiVWeW11l0qu8bgfk65oIqP7d2vO/G86n6Mpz+MS2hVk/beb3f/bQNz+L+Zv28Mum7/hiwXf8fkI0Iy9JYHfPkfQ77mRO7deZ2C7tgqIGobXGvmMvS9bv4Ps/drDupwU45k3ltz+L0MCKTop7T4zma30yv3QaQ///nM/7g3t6FhjcusbY+EVoCrpEJESwa9cqgmtP7cO1p/Zh3b97+HRJNO3aRPL3L9+Su347A9vs4e6wD+CvD/h7axTvvn4ES7a356hhF3Hc8HM44ahoBnXvSMc2Dd+3tGv/AX7/eze2jTtYueQr/lq5hF/63uDZ3rplD/K3FtEiHE7qG8WRx5/KtLgHuWzYIC7v1r7B4xMCJBEJUS8DurUnZeww9GUf8OuWXXz2/S/sXmVhdtkqRoYv50hVwJZVf/HZz6XwuYWOrRTDzh7N78ffTg9Ta2K6tKVvxzD6HtmZIztG0q1DJEd0jMTUuiWRLcJqrEXtLy1n594SCvaW8veuYv4s2Mf2bf9Qsi2fiJ1r6Fa8gXezF/DbxkJ2lThHyJ5485k42jvvEd/fqjMXTriNk84ey5hhJ3BSr06Ey+AD0cgkEQkRAEopBh9lYvCVI4GRbNq5ly/X/MOG321s7z2T04uW8cfmHfxbVE5F604AbHXsZ6tjP5u/TOaV9X/RtX0rTO1bc/7pg+jW8yj20obSFu1Zm7+V3UX7UeFhgKIg9gJ2tu7N/gPllJZVcPQPD+D4+0927dnH344SFv+nDSd0c3XyRMDrRUXsKtF07xBO3x6dacFuhh7bjdP7dea0fp2J7XJ+UDQbiuZLEpEQDaB3dFtuOD0WTo+lYkIC67btYdnGQlb98jNrilrRwqE4UO6soezZtZvC/ZrC/cWwrZgnhtuIj/jNuSMNo5fsxWov9+x71BXR7Ord3fN855aN5G7Y5Xn+e2E4LbocxR+6J+vpxUmJbTkj9hTMQ4dj7m2if9f2UusRQUUSkRANLCxMcfQRHZzLHwzrDUBpWQWbC/Zi376X/HPz+POPFezctIaiv/PJ6dqW5eUHaK330p79dOuzmrPa7kFr51RDJa27evbdMjyM7iOu5soz99K2W2+6DoijIOYYNnRux7Hd2jEmui0twpvvbBCiaZBEJIQBWkaE0a9re/p1bQ+DjoCR/Q8pU1JWTlFxGWeXVVBRoanQGq2hVYswIiPCiWwR7upHGmPAGQgROJKIhAhSrSLCadUu+O9FEqK+pM4uhBDCUJKIhBBCGEoSkRBCCENJIhJCCGEoSURCCCEMJYlICCGEoSQRCSGEMJQkIiGEEIaSRCSEEMJQkoiEEEIYShKREEIIQ0kiEkIIYShJREIIIQwVdIlIKRWjlEpRSsW7/jXVUN6slMqroUxGQIMUQggRMEGXiIAMrXW61toKZANphyuolIp3PTRXU8YMJAU2RCGEEIESVIlIKRXj/VxrbQfGHa681tqqtbbVsNsYwFH/6IQQQjSEoEpEOGs2jqovVk1QtaWUStBaZ9c3KCGEEA0n2BJRlI/XCgCTvztyJS97fQMSQgjRsEJ5qXBzbWpDSqkkDvYhlSilVjZsWE1GZ2CH0UEECbkWB8m1OEiuxUED6/PmRklErg/72GqKWFyDE3zVfqLws4/HNYjBWpuyWutMINP1vlyt9RB/jhWq5FocJNfiILkWB8m1OEgplVuf9zdKInJ92NeGDR/Nc65BC/4ap5RyPza5kqG1jvsSQgjRQIKqaU5rbfdKHu5+njlVnhdorR017KdSbUgpleFHMhRCCNGIgm2wAkCi+4ZWIEFrney1LQ2v4dzum15dj9O87itybzdV2V7T6DtJVgfJtThIrsVBci0OkmtxUL2uhdJaByoQIYQQwm/BWCMSQgjRjARVH1FjcDXPJeAcGGEGMg/X5+RP2abIz2thBtxNn0OBiaFyLer6/+zqe0yuqVxT4u+1UEol4BxgZIdD+2ebsjp8VsTjHPkbU13Zpsj19z9Dax1XQ7m6fWZqrZvVD86h4u7HMTjntqt32ab4U9vzwzmkPsnreQKQZ3T8RvxOeJUzO/98jI/fqGvh+j1I8SobMr8TdbgWKVWepxkdfwCvQ3xtf9/r+pnZrJrm/JnLzt9575oaP89vCJDq9dwKmGuaGb0pqMf/c8jNYViHa5GmtU53l9U1fFtuSupwLUZXeW4KdExG0bWb07Nen5nNKhHh31x2AZ33LgjV+vy0s7kl0eulGNfrh7y/CfL7/zmE5zCs9bVwNdU4vB6Hmrr8XlhcI3XjgawGjC1Y1fkzs7klIn/msgvYvHdByq/zq/KNaDyQ3gAxGcGv6xDicxj6cy1igAJXH5Hd65aLUOHv33+i6z0bcE4vFjJ9ZX6o82dmsxusIOrH1Rxn1lpXbYpoLmo1h2EzEAXEu38PlFKZOD+EOxkalXHczdcxQIZSCnezpahZc6sR+TOXXUDmvQtidT2/NCo30zV1tb4O/sxh2ET5+/fhqSW7mmlNIdR07c/vRQyuWpB2zuASC0wJhT5UP9X5M7O51Yj8mcsukPPeBSO/z881S0Wq1tqhlDKFSB+Rv9chlOcw9PfvoyoHzg+jUODPtaj0BUU7pyrLJLS+uNZGnT8zm1WNqOoF8TWXnftbTE1lmzp/roXreQKQ7ZV8QqI/wM/fCavWOtP943otM0SSUF3+Pgrcz919ZyHy5cTfvw8rPv4eQuX3ojqB+sxsdlP8VL3hyrsdVymVhXMcfGZNZUNBba+Fq1x+lbfbtdbVLe3RZPjzO+F6zYRzDas0nIM2MkLlQ8fPvw8TMAXn70YsMDVUEhH4fS0ScPYP2XHWCkKlluxukjZz8PfdvWxPwD4zm10iEkIIEVyaVdOcEEKI4COJSAghhKEkEQkhhDCUJCIhhBCGkkQkhBDCUJKIhBBCGEoSkRBCCENJIhJCCGEoSURCCCEMJYlICCGEoZrb7NtCNCnuJQaAoUCG63FMqM17KJo3qREJEdziXQvxLcM5uWo2kGxwTEIElCQiIYKbexr9oUAWQKjMei6Em8y+LUQToJTKdycgpVRMqCwxIARIH5EQQcu1Iq4DyHX96+4zcq97I0RIkBqREEHKtSCZCWfSGYJrGW5XP5EQIUMSkRBCCEPJYAUhhBCGkkQkhBDCUJKIhBBCGEoSkRBCCENJIhJCCGEoSURCCCEMJYlICCGEoSQRCSGEMNT/A5Fq0YT1Kn3PAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(X, Phi1.model.plot(), linewidth=3, label=\"$\\hat{f}(x)$ from $\\phi_{LCG}$\")\n",
    "ax.plot(X, Phi2.model.plot(), linewidth=2, dashes=[6, 2], label=\"$\\hat{f}(x)$ from $\\phi_{SymLCG}$\")\n",
    "ax.plot(X, f_an, 'k', linewidth=2, dashes=[2, 2], label=\"Analytical solution\")\n",
    "\n",
    "ax.legend()\n",
    "ax.set_xlabel(r\"$x$\"), ax.set_ylabel(\"$f(x)$\")\n",
    "ax.set_ylim(-0.14, 0), ax.set_xlim(0,1)\n",
    "\n",
    "plt.suptitle(\"Inverted functions\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's try to solve both inversions using the inverse of $D_2$ as a preconditioner"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# P = [D_2]^-1\n",
    "P = o.Matrix(o.VectorNumpy(np.linalg.inv(D2.matrix.getNdArray())), f, y)\n",
    "\n",
    "Phi1Prec = o.LeastSquares(f.clone(), y, D2, prec=P * P)\n",
    "Phi2Prec = o.LeastSquaresSymmetric(f.clone(), y, D2, prec=P)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "LCGsolver.setDefaults() # Re-setting default solver values\n",
    "SLCG.setDefaults() # Re-setting default solver values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "##########################################################################################\n",
      "            Preconditioned CG Solver\n",
      "    Restart folder: /tmp/restart_2022-04-22T01-50-53.913509/\n",
      "##########################################################################################\n",
      "\n",
      "iter = 0000, obj = 1.00000e+02, rnorm = 1.41e+01, gnorm = 5.66e+04, feval = 00002\n",
      "iter = 0001, obj = 1.26970e-22, rnorm = 1.59e-11, gnorm = 2.10e-06, feval = 00004\n",
      "iter = 0002, obj = 7.95245e-47, rnorm = 1.26e-23, gnorm = 3.61e-19, feval = 00006\n",
      "Objective function didn't reduce, will terminate solver:\n",
      "\tobj_new = 7.95245e-47\tobj_cur = 7.95245e-47\n",
      "##########################################################################################\n",
      "            Preconditioned CG Solver end\n",
      "##########################################################################################\n",
      "\n"
     ]
    }
   ],
   "source": [
    "LCGsolver.run(Phi1Prec, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "##########################################################################################\n",
      "            Preconditioned Symmetric CG Solver\n",
      "    Restart folder: /tmp/restart_2022-04-22T01-50-53.913671/\n",
      "##########################################################################################\n",
      "\n",
      "iter = 0000, obj = 0.00000e+00, rnorm = 1.41e+01, feval = 00001\n",
      "iter = 0001, obj = 8.45875e+00, rnorm = 1.69e-11, feval = 00003\n",
      "Objective function variation not monotonic, will terminate solver:\n",
      "\tobj_old=0.00000e+00\tobj_cur=8.45875e+00\tobj_new=8.45875e+00\n",
      "##########################################################################################\n",
      "            Preconditioned Symmetric CG Solver end\n",
      "##########################################################################################\n",
      "\n"
     ]
    }
   ],
   "source": [
    "SLCG.run(Phi2Prec, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, we converge to the global minimum in effectively one iteration."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}