{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computer Session 3\n",
    "\n",
    "\n",
    "\n",
    "## A first geometric optimisation problem\n",
    "\n",
    "We consider a convex polygon inscribed in the unit disk in $\\R^d$. In other\n",
    "words, we consider $0= \\theta_1\\leq \\theta_2\\leq \\dots\\leq \\theta_n< 2\\pi$, and\n",
    "the associated polygon is the convex hull of $\\{e^{i\\theta_k}\\}_{k=1,\\dots n}$.\n",
    "\n",
    "Find numerically the polygon that maximize the perimeter. Your function to maximize must\n",
    "take, as an argument, the geometric quantities of the problem: you should fix the number of sides or of angles. Then, computing the\n",
    "perimeter, minimise it under constraints, using (for instance) a penalised\n",
    "method (Exercise 3.2). Compute the gradient by finite differences.\n",
    "\n",
    "Plot your results graphically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 72508 characters]",
      "text/plain": [
       "<Figure size 700x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "def generate_thetas(n):\n",
    "    random_steps = np.random.random(n)\n",
    "    return np.concatenate(\n",
    "        ([0], np.cumsum(random_steps / np.sum(random_steps) * (2 * np.pi)))\n",
    "    )\n",
    "\n",
    "\n",
    "n = 4\n",
    "thetas = generate_thetas(n)\n",
    "thetas_inf = np.linspace(0, 2 * np.pi, 1000)\n",
    "\n",
    "plt.figure(figsize=(7, 7))\n",
    "plt.plot(np.cos(thetas), np.sin(thetas), label=\"polygon\")\n",
    "plt.scatter(np.cos(thetas), np.sin(thetas), color=\"red\", label=\"vertices\")\n",
    "plt.plot(np.cos(thetas_inf), np.sin(thetas_inf), \"k--\", label=\"unit circle\", alpha=0.5)\n",
    "plt.legend()\n",
    "plt.title(f\"Polygon with {n} sides\")\n",
    "plt.grid(True)\n",
    "plt.xlim(-1.1, 1.1)\n",
    "plt.ylim(-1.1, 1.1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 67416 characters]",
      "text/plain": [
       "<Figure size 700x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal angles (radians): [0.         2.0943487  4.1884815  6.28318531 6.28318531]\n",
      "Maximum perimeter: 3.463947248823501\n",
      "2 * pi = 6.283185307179586\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 58576 characters]",
      "text/plain": [
       "<Figure size 700x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal angles (radians): [0.         1.57079587 3.14158997 4.71238704 6.28318531 6.28318531]\n",
      "Maximum perimeter: 4.242639312441127\n",
      "2 * pi = 6.283185307179586\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 70772 characters]",
      "text/plain": [
       "<Figure size 700x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal angles (radians): [0.         1.25678882 2.51244268 3.76820788 5.02645179 6.28318531\n",
      " 6.28318531]\n",
      "Maximum perimeter: 4.702203347007718\n",
      "2 * pi = 6.283185307179586\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 69316 characters]",
      "text/plain": [
       "<Figure size 700x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal angles (radians): [0.         1.04746453 2.09466627 3.14153418 4.188399   5.2357538\n",
      " 6.28318531 6.28318531]\n",
      "Maximum perimeter: 4.999797344580875\n",
      "2 * pi = 6.283185307179586\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 71700 characters]",
      "text/plain": [
       "<Figure size 700x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal angles (radians): [0.         0.89680577 1.79461604 2.69230295 3.58961474 4.4879636\n",
      " 5.38591918 6.28318531 6.28318531]\n",
      "Maximum perimeter: 5.206903626576404\n",
      "2 * pi = 6.283185307179586\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 73900 characters]",
      "text/plain": [
       "<Figure size 700x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal angles (radians): [0.         0.78574785 1.57141937 2.35536796 3.14089114 3.92706882\n",
      " 4.71302233 5.49934121 6.28318531 6.28318531]\n",
      "Maximum perimeter: 5.359003436175853\n",
      "2 * pi = 6.283185307179586\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 71144 characters]",
      "text/plain": [
       "<Figure size 700x700 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal angles (radians): [0.         0.69852563 1.39640614 2.09468989 2.79342176 3.49154254\n",
      " 4.18880895 4.88742344 5.58668458 6.28318531 6.28318531]\n",
      "Maximum perimeter: 5.473854663395368\n",
      "2 * pi = 6.283185307179586\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "\n",
    "def polygon_perimeter(theta, n):\n",
    "    points = np.array([[np.cos(t), np.sin(t)] for t in theta])\n",
    "    perimeter = 0\n",
    "    for i in range(n - 1):\n",
    "        perimeter += np.sqrt(np.sum((points[i + 1] - points[i]) ** 2))\n",
    "    perimeter += np.sqrt(np.sum((points[0] - points[-1]) ** 2))\n",
    "    return -perimeter\n",
    "\n",
    "\n",
    "def constraint_increasing(theta):\n",
    "    return np.array([theta[i + 1] - theta[i] for i in range(len(theta) - 1)])\n",
    "\n",
    "\n",
    "def optimize_polygon(n):\n",
    "    theta0 = generate_thetas(n)\n",
    "\n",
    "    constraints = [\n",
    "        {\"type\": \"ineq\", \"fun\": constraint_increasing},\n",
    "        {\"type\": \"eq\", \"fun\": lambda x: x[0]},\n",
    "        {\"type\": \"ineq\", \"fun\": lambda x: 2 * np.pi - x[-1]},\n",
    "    ]\n",
    "\n",
    "    result = minimize(\n",
    "        lambda x: polygon_perimeter(x, n),\n",
    "        theta0,\n",
    "        constraints=constraints,\n",
    "        method=\"SLSQP\",\n",
    "    )\n",
    "\n",
    "    return result.x\n",
    "\n",
    "\n",
    "def plot_perimeter(n):\n",
    "    optimal_angles = optimize_polygon(n + 1)\n",
    "    plt.figure(figsize=(7, 7))\n",
    "    t = np.linspace(0, 2 * np.pi, 100)\n",
    "    plt.plot(np.cos(t), np.sin(t), \"k--\", alpha=0.5, label=\"unit circle\")\n",
    "\n",
    "    points = np.array([[np.cos(t), np.sin(t)] for t in optimal_angles])\n",
    "    points = np.vstack([points, points[0]])\n",
    "    plt.plot(points[:, 0], points[:, 1], \"b-\", linewidth=2, label=\"optimal polygon\")\n",
    "    plt.scatter(points[:-1, 0], points[:-1, 1], color=\"red\", label=\"vertices\")\n",
    "\n",
    "    plt.legend()\n",
    "    plt.axis(\"equal\")\n",
    "    plt.grid(True)\n",
    "    plt.title(f\"Optimal {n}-sided Polygon Inscribed in Unit Circle\")\n",
    "    plt.xlabel(\"x\")\n",
    "    plt.ylabel(\"y\")\n",
    "    plt.axis(\"equal\")\n",
    "    plt.grid(True)\n",
    "    plt.show()\n",
    "\n",
    "    print(f\"Optimal angles (radians): {optimal_angles}\")\n",
    "    print(f\"Maximum perimeter: {-polygon_perimeter(optimal_angles, n)}\")\n",
    "    print(f\"2 * pi = {2 * np.pi}\")\n",
    "\n",
    "\n",
    "for n in np.arange(3, 10, 1):\n",
    "    plot_perimeter(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Projection on convex sets $K=\\{x\\in\\R^d, Cx\\leq b\\}$\n",
    "\n",
    "Let compute the projection of a point $x_0$ on the convex set $K=\\{x\\in\\R^d,\n",
    "Cx\\leq b\\}$, where $C$ is a matrix and $b$ is a vector. The projection is\n",
    "defined as: \n",
    "$$\n",
    "\\Pi_K(x_0) = \\arg\\min_{x\\in K} \\|x-x_0\\|^2.$$\n",
    "\n",
    "1. Define the projection operator when $C=\\mathrm{Id}$ and $b=0$. We denote by $\\Pi_{\\mathbb{K}_-}$ this projection operator. \n",
    "   How can you compute it using Python?\n",
    "2. Write the KKT conditions for the quadratic optimization problem when $K=\\{x\\in\\R^d, Cx\\leq b\\}$: \n",
    "    $$\n",
    "    \\min_{x\\in K}\\left\\{ \\frac12 ||x||^2-\\langle x_0,x\\rangle\\right\\}.$$\n",
    "    How is it related to the projection operator?\n",
    "3. Show that, for all $\\tau>0$, we have \n",
    "    $$\\forall \\mu\\in\\mathrm{K}_-, \\quad\n",
    "    \\langle\\lambda^*-\\mu,\\lambda^*-(\\lambda^*+\\tau(b-Cx^*))\\rangle\\leq 0.$$\n",
    "    Express this condition in terms of the projection operator\n",
    "    $\\Pi_{\\mathbb{K}_-}$.\n",
    "4. Implement Uzawa’s algorithm to compute numerically the projection operator on $K=\\{x\\in\\R^d, Cx\\leq b\\}$:\n",
    "    $$\\lambda_{n+1}=\\Pi_{\\mathbb{K}_-}(\\lambda_n+\\tau(b-Cx_n)),\\text{ then }\n",
    "    x_{n+1}=C^T\\lambda_{n+1}+x_0.$$\n",
    "5. You can adapte the code of the first exercise to use the projection operator. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial point: [ 2 -1  3  0 -5]\n",
      "Projection of x0 onto K: [ 0 -1  0  0 -5]\n"
     ]
    }
   ],
   "source": [
    "x0 = np.array([2, -1, 3, 0, -5])\n",
    "\n",
    "\n",
    "def K(x):\n",
    "    return np.minimum(x, 0)\n",
    "\n",
    "\n",
    "print(f\"Initial point: {x0}\")\n",
    "print(f\"Projection of x0 onto K: {K(x0)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question 2**: KKT Conditions\n",
    "\n",
    "First, let's identify the components:\n",
    "\n",
    "Objective function: $f(x) = \\frac{1}{2}||x||^2 - \\langle x_0, x \\rangle$\n",
    "\n",
    "Inequality constraints: $g_i(x) = (Cx - b)_i \\leq 0$ for $i = 1, 2, ..., m$ where $m$ is the number of rows in $C$\n",
    "\n",
    "The Lagrangian is: $$L(x, \\lambda) = \\frac{1}{2}||x||^2 - \\langle x_0, x \\rangle + \\lambda^T(Cx - b)$$\n",
    "\n",
    "The KKT conditions are:\n",
    "\n",
    "- Stationarity: $\\nabla_x L(x^*, \\lambda^*) = 0$ $$x^* - x_0 + C^T\\lambda^* = 0$$ $$\\Rightarrow x^* = x_0 - C^T\\lambda^*$$\n",
    "- Primal feasibility: $Cx^* \\leq b$\n",
    "- Dual feasibility: $\\lambda^* \\geq 0$\n",
    "- Complementary slackness: $\\lambda_i^*(Cx^* - b)_i = 0$ for all $i$\n",
    "\n",
    "The original projection problem is: $$\\Pi_K(x_0) = \\arg \\min_{x\\in K} ||x-x_0||^2 $$\n",
    "\n",
    "This is equivalent to: $$\\arg \\min_{x\\in K} \\frac{1}{2}||x-x_0||^2 = \\arg \\min_{x \\in K} \\frac{1}{2}||x||^2 - \\langle x_0, x \\rangle + \\frac{1}{2}||x_0||^2$$\n",
    "\n",
    "Since $\\frac{1}{2}||x_0||^2$ is a constant that doesn't affect the minimization, this is equivalent to: $$\\arg\\min_{x\\in K} \\frac{1}{2}||x||^2 - \\langle x_0, x \\rangle$$\n",
    "\n",
    "Therefore, the solution $x^*$ to our KKT system is precisely the projection of $x_0$ onto $K$, i.e., $\\Pi_K(x_0) = x^*$.\n",
    "\n",
    "From the stationarity condition, we have: $$x^* = x_0 - C^T\\lambda^*$$\n",
    "\n",
    "This shows that the projection is obtained by shifting the original point $x_0$ by $C^T\\lambda^*$, where $\\lambda^*$ are the optimal Lagrange multipliers that satisfy all the KKT conditions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The polygonal isoperimetric inequality\n",
    "\n",
    "Given an perimeter constraint $L$ and a fixed number of sides $N$, what is the polygon with $N$ sides that maximises the enclosed area? We represent a polygon in $\\R^2$ by a collection of point $\\{(x_i,y_i)\\}_{i=1,\\dots,N}$, ordered in a counterclockwise fashion.\n",
    " We admit that the area of this polygon is given by \n",
    "$$A=\\frac12\\sum_{k=1}^N (x_ky_{k+1}-y_{k}x_{k+1}).$$\n",
    "\n",
    "\n",
    "You must enter as parameters the number of sides of the polygon you seek and use whichever algorithm you see fit (typically, a Uzawa method).\n",
    "\n",
    "Plot your results graphically."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The obstacle problem\n",
    "\n",
    "Let $g$ be a given continuous function on the interval $[0,1]$. We consider an \\textit{obstacle problem}: find a function $ u: [0,1] \\longrightarrow \\mathbb{R}$ such that:\n",
    "\n",
    "$$\n",
    "\\left\\{\n",
    "\\begin{array}{ll}\n",
    "-u''(x) \\geq 1 & x \\in (0,1) \\\\\n",
    "u(x) \\geq g(x) & x \\in (0,1) \\\\\n",
    "\\left( -u''(x) - 1 \\right)(u(x) - g(x)) = 0 & x \\in (0,1) \\\\\n",
    "u(0) = u(1) = 0 & \\\\\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "The first equation represents a minimum concavity for the function $u$, the second equation represents the obstacle: $u$ must remain above $g$. The third equation expresses the fact that we must satisfy at least one of the two previous equations with equality: either we solve $-u''(x) = 1$, or $u(x) = g(x)$, and we are on the obstacle."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Associated minimization problem\n",
    "\n",
    "We discretize this problem by introducing a uniform mesh: $x_j = jh$, where $h$ is the step size of the mesh, and $j \\in \\{0, \\dots, n+1\\}$, with $n \\geq 1$ an integer and $h = \\frac{1}{n+1}$. Let $ g_j = g(x_j) $ for $ j \\in \\{0, \\dots, n+1\\} $. We seek values $ u_j = u(x_j) $ for $ j \\in \\{0, \\dots, n+1\\} $ such that:\n",
    "\n",
    "$$\n",
    "\\left\\{\n",
    "\\begin{array}{ll}\n",
    "\\displaystyle -\\frac{u_{j-1} - 2u_j + u_{j+1}}{h^2} \\geq 1 & j \\in \\{0, \\dots, n+1\\} \\\\\n",
    "u_j \\geq g_j & j \\in \\{0, \\dots, n+1\\} \\\\\n",
    "\\displaystyle \\left( -\\frac{u_{j-1} - 2u_j + u_{j+1}}{h^2} \\right) (u_j - g_j) = 0 & j \\in \\{0, \\dots, n+1\\} \\\\\n",
    "u_0 = u_{n+1} = 0 & \\\\\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "Recall that $ -\\frac{u_{j-1} - 2u_j + u_{j+1}}{h^2} $ is the finite difference approximation of $ -u''(x_j) $.\n",
    "\n",
    "We introduce the matrix $ A \\in S_n(\\mathbb{R}) $, defined by:\n",
    "\n",
    "$$\n",
    "A = \\frac{1}{h^2} \\left( \\begin{array}{ccccc}\n",
    "2 & -1 & 0 & \\dots & 0 \\\\\n",
    "-1 & 2 & -1 & \\vdots \\\\\n",
    "0 & \\ddots & \\ddots & \\ddots & 0 \\\\\n",
    "\\vdots & & -1 & 2 & -1 \\\\\n",
    "0 & \\dots & 0 & -1 & 2\n",
    "\\end{array} \\right).\n",
    "$$\n",
    "\n",
    "We also define the column vectors $ b $ and $ g $ of $ \\mathbb{R}^n $ as follows:\n",
    "\n",
    "$$\n",
    "b = \\left( \\begin{array}{c}\n",
    "1 \\\\\n",
    "\\vdots \\\\\n",
    "1\n",
    "\\end{array} \\right), \\quad g = \\left( \\begin{array}{c}\n",
    "g_1 \\\\\n",
    "\\vdots \\\\\n",
    "g_n\n",
    "\\end{array} \\right).\n",
    "$$\n",
    "\n",
    "Recall that if $ u = \\left( \\begin{array}{c} u_1 \\\\ \\vdots \\\\ u_n \\end{array} \\right) $, then\n",
    "\n",
    "$$\n",
    "u \\text{ is a solution of (2)} \\Longleftrightarrow u \\text{ is a solution of } \\left\\{\n",
    "\\begin{array}{l}\n",
    "\\displaystyle \\min_{v \\in K} \\left\\{ \\frac{1}{2}(Av, v) - (b, v) \\right\\} \\\\\n",
    "K = \\{v \\in \\mathbb{R}^n : v \\geq g \\}\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "We define the function $J$ on $\\mathbb{R}^n$ by:\n",
    "\n",
    "$$\n",
    "J(v) = \\frac{1}{2} (Av, v) - (b, v).\n",
    "$$\n",
    "\n",
    "### Problem Resolution\n",
    "\n",
    "Consider the specific case where $g(x) = \\max \\left( 0, 1 - 100(x - 0.7)^2 \\right)$.\n",
    "\n",
    "We want to solve this problem using the projected gradient algorithm. Denote by $\\Pi_K =\\{v\\geq g\\}$ the projection onto the convex set $K$. Without proving it, we can use the fact that:\n",
    "\n",
    "$$\n",
    "\\Pi_K(v) = \\left( \\max(v_i, g_i) \\right)_{1 \\leq i \\leq n}.\n",
    "$$\n",
    "\n",
    "1. Write a gradient method with a fixed step size to find the minimum of $J$ over $\\mathbb{R}^n$. Test this method and verify that it converges to the desired result.\n",
    "2. Adapt the previous program to implement the projected gradient method with a constant step size. For instance, you can choose the optimal step size $\\rho_{\\textrm{opt}} = \\frac{2}{\\lambda_1(A) + \\lambda_n(A)}$ of the gradient method without constraints. Also, don't forget to set a maximum number of iterations and a relevant stopping criterion.\n",
    "3. Test the program for different values of $n$. Represent the graph of the solution and the graph of the obstacle on the same plot. Verify that if $u(x) \\neq g(x)$, then $-u''(x) = 1$ (using the finite difference approximation).\n",
    "4. Implement the Uzawa algorithm to solve this problem and compare the efficiency of the methods implemented. Test with different choices of obstacles.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "studies",
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
