{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 138,
   "id": "c03a657e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a462ed55",
   "metadata": {},
   "source": [
    "Méthode Newton $\\mathbb{R}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "id": "e90827f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f1(x):\n",
    "    f = 2*x-np.log(x)\n",
    "    g = 2 - 1/x\n",
    "    h = 1/(x**2)\n",
    "    return f, g, h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "id": "ca26160d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f2(x):\n",
    "    f = x**3 - 6*x**2 + 4*x + 2\n",
    "    g = 3*x**2 - 12*x + 4\n",
    "    h = 6*x - 12\n",
    "    return f, g, h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "id": "1f6bfe41",
   "metadata": {},
   "outputs": [],
   "source": [
    "def NewtonR(fct, x0, maxiter, tol=1e-9):\n",
    "    x = np.zeros(maxiter)\n",
    "    x[0] = x0\n",
    "\n",
    "    for k in range(maxiter-1):\n",
    "        f, g, h = fct(x[k])\n",
    "        x[k+1] = x[k] - g/h\n",
    "\n",
    "        if abs(g) < tol:\n",
    "            return x[:k+1]\n",
    "\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "id": "c012277d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1eec176bb10>]"
      ]
     },
     "execution_count": 142,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 20708 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Ici ça marche super bien\n",
    "\n",
    "absc = np.linspace(0.01, 1.5, 1000)\n",
    "oord = f1(absc)[0]\n",
    "\n",
    "x0 = 0.1\n",
    "\n",
    "x = NewtonR(f1, x0, 1000)\n",
    "for k in range(len(x)):\n",
    "    f, g, h = f1(x[k])\n",
    "    plt.plot(x[k], f, 'ob')\n",
    "\n",
    "plt.plot(absc, oord, 'g')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "id": "087255d2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1eec408c050>]"
      ]
     },
     "execution_count": 143,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 24932 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Ici ça dépend du point initial\n",
    "\n",
    "absc = np.linspace(-2, 5, 1000)\n",
    "oord = f2(absc)[0]\n",
    "\n",
    "x0 = [3,0]\n",
    "c = ['b', 'r']\n",
    "\n",
    "for i, x0 in enumerate(x0):\n",
    "    x = NewtonR(f2, x0, 1000)\n",
    "    for k in range(len(x)):\n",
    "        f, g, h = f2(x[k])\n",
    "        plt.plot(x[k], f, 'o' + c[i])\n",
    "\n",
    "plt.plot(absc, oord, 'g')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f68a7951",
   "metadata": {},
   "source": [
    "Examen d'implémentation | Gradient et Newton"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06964cdd",
   "metadata": {},
   "source": [
    "$f(x_1,x_2)=2x_1^2-1.05x_1^4+x_1^6/6+x_1x_2+x_2^2$\n",
    "\n",
    "$\\nabla f(x_1,x_2)=\n",
    "\\begin{pmatrix}\n",
    "4x_1 - 4.2x_1^3 + x_1^5 + x_2 & x_1 + 2x_2 \\\\\n",
    "\\end{pmatrix}$\n",
    "\n",
    "$\\nabla^2 f(x_1,x_2)=\n",
    "\\begin{pmatrix}\n",
    "4 - 12.6x_1^2 + 5x_1^4 & 1 \\\\\n",
    "1 & 2 \\\\\n",
    "\\end{pmatrix}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "id": "7ba5e5c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fexamen(x, nargout=2):\n",
    "    x1 = x[0]\n",
    "    x2 = x[1]\n",
    "    f = 2*x1**2 - 1.05*x1**4 + x1**6/6 + x1*x2 + x2**2 \n",
    "\n",
    "    g = np.array([4*x1 - 4*1.05*x1**3 + x1**5 + x2, x1 + 2*x2])\n",
    "\n",
    "    if nargout==3:\n",
    "        H = np.array([[4 - 12*1.05*x1**2 + 5*x1**4, 1],\n",
    "                      [              1            , 2]])\n",
    "        return f, g, H\n",
    "\n",
    "    return f, g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "id": "7fb63347",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(np.float64(1.0), array([1., 2.]))\n"
     ]
    }
   ],
   "source": [
    "# Vérification de fexamen: OK\n",
    "x = np.array([0,1])\n",
    "print(fexamen(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "id": "449e09c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def w1(f,x,d,alpha,beta1):\n",
    "    f0, g0 = f(x) # à l'itéré courant\n",
    "    f1, g1 = f(x+alpha*d) # à l'itéré potentiel\n",
    "    return f1 <= f0 + alpha*beta1 * np.dot(d, g0)\n",
    "\n",
    "def w2(f,x,d,alpha,beta2):\n",
    "    f0, g0 = f(x) # à l'itéré courant\n",
    "    f1, g1 = f(x+alpha*d) # à l'itéré potentiel   \n",
    "    return np.dot(d, g1) >= beta2*np.dot(d, g0)\n",
    "\n",
    "def wolfebissection(f,x,d,alpha,beta1=1e-3,beta2=0.9):\n",
    "    aleft = 0\n",
    "    aright = np.inf\n",
    "\n",
    "    while True:\n",
    "        w_1 = w1(f,x,d,alpha,beta1)\n",
    "        w_2 = w2(f,x,d,alpha,beta2)\n",
    "\n",
    "        if w_1 and w_2:\n",
    "            break\n",
    "\n",
    "        elif not w_1:\n",
    "            aright = alpha\n",
    "            alpha = (aleft+aright)/2\n",
    "\n",
    "        elif not w_2:\n",
    "            aleft = alpha\n",
    "            if aright<np.inf:\n",
    "                alpha = (aleft+aright)/2\n",
    "            else:\n",
    "                alpha *= 2 # 2 est arbitraire\n",
    "\n",
    "    return alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "id": "3d5c52c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def methgradient(f,x0,maxiter=300,tol=1e-6):\n",
    "    x = x0.copy()\n",
    "    fobj = np.zeros(maxiter)\n",
    "    ngrad = np.zeros(maxiter)\n",
    "\n",
    "    for k in range(maxiter):\n",
    "        fx, g = f(x)\n",
    "        d = -g\n",
    "        alpha = wolfebissection(f,x,d,1)\n",
    "\n",
    "        x = x+alpha*d # calcul nouvel itéré\n",
    "\n",
    "        # sauvegarde des f et g\n",
    "        fx, g = f(x)\n",
    "        fobj[k] = fx\n",
    "        ngrad[k] = np.linalg.norm(g)\n",
    "\n",
    "        # critère d'arrêt\n",
    "        if np.linalg.norm(g) <= tol:\n",
    "            break\n",
    "\n",
    "    return x, fobj[:k+1], ngrad[:k+1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "id": "c2a657c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def methNewton(f,x0,maxiter=300,tol=1e-6):\n",
    "    x = x0.copy()\n",
    "    fobj = np.zeros(maxiter)\n",
    "    ngrad = np.zeros(maxiter)\n",
    "\n",
    "    for k in range(maxiter):\n",
    "        fx, g, H = f(x, nargout=3) # besoin de H\n",
    "        z = np.linalg.solve(H,g)\n",
    "\n",
    "        x = x-z # calcul nouvel itéré\n",
    "\n",
    "        # sauvegarde des f et g\n",
    "        fx, g = f(x)\n",
    "        fobj[k] = fx\n",
    "        ngrad[k] = np.linalg.norm(g)\n",
    "\n",
    "        # critère d'arrêt\n",
    "        if np.linalg.norm(g) <= tol:\n",
    "            break\n",
    "\n",
    "    return x, fobj[:k+1], ngrad[:k+1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "id": "274c9c44",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 31704 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 31296 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x0 = np.array([0.1, 0.1])\n",
    "\n",
    "xg, fg, ngg = methgradient(fexamen,x0)\n",
    "xn, fn, ngn = methNewton(fexamen,x0)\n",
    "\n",
    "plt.semilogy(fg, label=\"gradient\")\n",
    "plt.semilogy(fn, label=\"Newton\")\n",
    "plt.title(\"Evolution de la fonction objectif\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "plt.semilogy(ngg, label=\"gradient\")\n",
    "plt.semilogy(ngn, label=\"Newton\")\n",
    "plt.title(\"Evolution de la norme du gradient\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "id": "1016c4ec",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Avec la méthode du gradient : x=[-9.13230123e-08  9.55910628e-08]\n",
      "fobj=1.7087772621488852e-14 | grad=[-2.69700986e-07  9.98591133e-08] | H=[[4. 1.]\n",
      " [1. 2.]]\n",
      "EigResult(eigenvalues=array([4.41421356, 1.58578644]), eigenvectors=array([[ 0.92387953, -0.38268343],\n",
      "       [ 0.38268343,  0.92387953]]))\n",
      "Avec la méthode de Newton : x=[ 3.64930450e-08 -1.82465225e-08]\n",
      "fobj=2.3305490890415854e-15 | grad=[1.27725658e-07 4.33680869e-19] | H=[[4. 1.]\n",
      " [1. 2.]]\n",
      "EigResult(eigenvalues=array([4.41421356, 1.58578644]), eigenvectors=array([[ 0.92387953, -0.38268343],\n",
      "       [ 0.38268343,  0.92387953]]))\n"
     ]
    }
   ],
   "source": [
    "print(f\"Avec la méthode du gradient : x={xg}\")\n",
    "f, g, H = fexamen(xg, nargout=3)\n",
    "print(f\"fobj={f} | grad={g} | H={H}\")\n",
    "print(np.linalg.eig(H))\n",
    "# On a un point ou le gradient vaut 0 et la matrice hessienne est definie positive\n",
    "# => on a un minimum local\n",
    "\n",
    "print(f\"Avec la méthode de Newton : x={xn}\")\n",
    "f, g, H = fexamen(xn, nargout=3)\n",
    "print(f\"fobj={f} | grad={g} | H={H}\")\n",
    "print(np.linalg.eig(H))\n",
    "# On a un point ou le gradient vaut 0 et la matrice hessienne est indéfinie\n",
    "# => on a un point selle"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf54a678",
   "metadata": {},
   "source": [
    "$f(x)=(x_1-2x_2)^4 + (x_2-3x_3)^4 + (x_3-4x_4)^4 + \\ldots + (x_n-(x_1+1))^4$\n",
    "\n",
    "$x_1 = 2x_2 = 2*3*4(x_1+1) => -23x_1=24 => x_1=-24/23$\n",
    "\n",
    "$x_2 = 3x_3 = 3*4(x_1+1) => x_2 = \\ldots$\n",
    "\n",
    "$\\ldots$\n",
    "\n",
    "$x_i = \\frac{-n!}{i(n!-1)}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "id": "1a83735f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fsize(x,nargout=2):\n",
    "  n = len(x)\n",
    "\n",
    "  #Calcul de f\n",
    "  v      = 1.0*x\n",
    "  v[:-1] = v[:-1]-x[1:]*np.arange(2,1+len(x))\n",
    "  v[-1]  = x[-1]-x[0]-1\n",
    "  f      = np.sum(v**4)\n",
    "\n",
    "  #Calcul du gradient uniquement s'il est demandé\n",
    "  g     = 4*v**3\n",
    "  g[1:] = g[1:]-4*(v[:-1]**3)*np.arange(2,1+len(x))\n",
    "  g[0]  = g[0]-4*v[-1]**3\n",
    "\n",
    "  #Calcul de la hessienne uniquement si elle est demandée\n",
    "  if nargout==3:\n",
    "    H = np.zeros((n,n))\n",
    "    for i in range(1,n):\n",
    "      H[i-1:i+1,i-1:i+1] = H[i-1:i+1,i-1:i+1]+ 12*v[i-1]**2*np.array([[1,-(i+1)],[-(i+1),(i+1)**2]])\n",
    "    H[0,0]   += 12*v[-1]**2\n",
    "    H[0,-1]  -= 12*v[-1]**2\n",
    "    H[-1,0]  -= 12*v[-1]**2\n",
    "    H[-1,-1] += 12*v[-1]**2\n",
    "  if nargout==3:\n",
    "    return f, g, H\n",
    "  else:\n",
    "    return f, g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "id": "1f5737dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def methgradientTime(f,x0,maxtime,maxiter=300):\n",
    "    start_time = time.time()\n",
    "    x = x0.copy()\n",
    "    fobj = np.zeros(maxiter)\n",
    "    temps = np.zeros(maxiter)\n",
    "\n",
    "    for k in range(maxiter):\n",
    "        fx, g = f(x)\n",
    "        d = -g\n",
    "        alpha = wolfebissection(f,x,d,1)\n",
    "\n",
    "        x = x+alpha*d # calcul nouvel itéré\n",
    "\n",
    "        # sauvegarde des f et g\n",
    "        fx, g = f(x)\n",
    "        fobj[k] = fx\n",
    "        temps[k] = time.time() - start_time\n",
    "\n",
    "        if temps[k] > maxtime:\n",
    "            break\n",
    "\n",
    "        if np.linalg.norm(g) <= 1e-6:\n",
    "            break\n",
    "\n",
    "    return x, fobj[:k+1], temps[:k+1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "id": "0ee9c336",
   "metadata": {},
   "outputs": [],
   "source": [
    "def methNewtonTime(f,x0,maxtime,maxiter=300):\n",
    "    start_time = time.time()\n",
    "    x = x0.copy()\n",
    "    fobj = np.zeros(maxiter)\n",
    "    temps = np.zeros(maxiter)\n",
    "\n",
    "    for k in range(maxiter):\n",
    "        fx, g, H = f(x, nargout=3) # besoin de H\n",
    "        z = np.linalg.solve(H,g)\n",
    "\n",
    "        x = x-z # calcul nouvel itéré\n",
    "\n",
    "        # sauvegarde des f et g\n",
    "        fx, g = f(x)\n",
    "        fobj[k] = fx\n",
    "        temps[k] = time.time() - start_time\n",
    "\n",
    "        if temps[k] > maxtime:\n",
    "            break\n",
    "\n",
    "        # critère d'arrêt\n",
    "        if np.linalg.norm(g) <= 1e-6:\n",
    "            break\n",
    "\n",
    "    return x, fobj[:k+1], temps[:k+1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4b1b445",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 31496 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34076 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# La méthode de Newton a le temps de calculer moins d'itérés que le gradient\n",
    "# à cause du coût de calcul plus élevé. En terme de temps, le gradient est meilleur.\n",
    "\n",
    "n = 5000\n",
    "x0 = np.ones(n)\n",
    "\n",
    "maxtime = 15\n",
    "xg, fg, tg = methgradientTime(fsize, x0, maxtime)\n",
    "xn, fn, tn = methNewtonTime(fsize, x0, maxtime)\n",
    "\n",
    "plt.semilogy(fg, label=\"gradient\")\n",
    "plt.semilogy(fn, label=\"Newton\")\n",
    "plt.title(\"Evolution de la fonction objectif en fonction des itérations\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "plt.semilogy(tg, fg, label=\"gradient\")\n",
    "plt.semilogy(tn, fn, label=\"Newton\")\n",
    "plt.title(\"Evolution de la norme du gradient, en fonction du temps\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "optinonlin",
   "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.13.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
