{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def L(x, l):\n",
    "    return x**2 + 1 + l*(x-2)*(x-4)\n",
    "\n",
    "ls = np.arange(0, 10, 1)\n",
    "x = np.arange(-1, 5, .1)\n",
    "plt.figure(figsize=(7,7))\n",
    "plt.grid()\n",
    "for l in ls:\n",
    "    y = L(x, l)\n",
    "    plt.plot(x, y, label = 'lambda = {}'.format(l), alpha = .2)\n",
    "    \n",
    "    plt.plot(x, 5*np.ones(x.shape), 'k', lw = 2)    \n",
    "    ind_min = np.argmin(y)\n",
    "    plt.plot(x[ind_min], y[ind_min], 'o')\n",
    "\n",
    "plt.ylim([0,10])\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('f(x)')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dual(l):\n",
    "    return -9*l**2/(1+l) + 8*l + 1\n",
    "\n",
    "l = np.arange(0, 10, .1)\n",
    "y = dual(l)\n",
    "plt.figure(figsize=(7,7))\n",
    "plt.grid()\n",
    "plt.plot(l, y)\n",
    "plt.xlabel('lambda')\n",
    "plt.ylabel('l(lambda)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 2: Log Berrier Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradDescentBacktraking(x, f, grad_f):\n",
    "    alpha = 1\n",
    "    fx = f(x)\n",
    "    while True:\n",
    "        grad_fx = grad_f(x)\n",
    "        delta = -grad_fx/np.linalg.norm(grad_fx)\n",
    "        fx_next = f(x+alpha*delta)\n",
    "        while fx_next > fx + 0.5*grad_fx.T@(alpha*delta) or np.isnan(fx_next):\n",
    "            alpha *= 0.5\n",
    "            fx_next = f(x+alpha*delta)\n",
    "            \n",
    "        x += alpha*delta\n",
    "        fx = fx_next\n",
    "        \n",
    "        alpha = min(1.2*alpha, np.inf)\n",
    "        if np.linalg.norm(alpha*delta) < 1e-10:\n",
    "            break\n",
    "            \n",
    "    return x, fx\n",
    "\n",
    "def F(x, mu):\n",
    "    f = x.sum()\n",
    "    g1 = x.T@x-1\n",
    "    g2 = -x[0]\n",
    "    return f - mu*np.log(-g1) - mu*np.log(-g2)\n",
    "\n",
    "def dF(x, mu):\n",
    "    g1 = x.T@x-1\n",
    "    g2 = -x[0]\n",
    "    \n",
    "    df = np.array([1,1])\n",
    "    dg1 = 2*x\n",
    "    dg2 = np.array([-1,0])\n",
    "    return df - mu/g1*dg1 - mu/g2*dg2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import partial\n",
    "\n",
    "x = np.array([0.5, 0.5])\n",
    "g = np.array([x.T@x-1, -x[0]])\n",
    "mu = 1.\n",
    "l = -mu/g\n",
    "xs, ls = [x.copy()], [l.copy()]\n",
    "for i in range(10):\n",
    "    F_mu = partial(F, mu=mu)\n",
    "    dF_mu = partial(dF, mu=mu)\n",
    "    \n",
    "    x, fx = gradDescentBacktraking(x, F_mu, dF_mu) # inner loop\n",
    "    g = np.array([x.T@x-1, -x[0]])\n",
    "    l = -mu/g\n",
    "    xs.append(x.copy())\n",
    "    ls.append(l.copy())\n",
    "    mu *= 0.5\n",
    "    \n",
    "xs = np.stack(xs)\n",
    "ls = np.stack(ls)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(7,7))\n",
    "plt.title('central path')\n",
    "plt.axis('equal')\n",
    "plt.plot(xs[0,0], xs[0,1], '*')\n",
    "plt.plot(xs[:,0], xs[:,1], '*-')\n",
    "plt.grid()\n",
    "plt.xlabel('x1')\n",
    "plt.ylabel('x2')\n",
    "plt.show()\n",
    "\n",
    "plt.figure(figsize=(7,7))\n",
    "plt.title('stationarity (force balance)')\n",
    "plt.axis('equal')\n",
    "plt.axis([-3, 3,-4, 2])\n",
    "\n",
    "xstar = xs[-1,:]\n",
    "plt.plot(xstar[0], xstar[1], '*')\n",
    "\n",
    "labmda1_dg1 = ls[-1,0]*2*xs[-1,:]\n",
    "labmda2_dg2 = ls[-1,1]*np.array([-1,0])\n",
    "df = np.array([1,1])\n",
    "fg1 = plt.arrow(xs[-1,0], xs[-1,1], -labmda1_dg1[0], -labmda1_dg1[1], color='b', width = .1, label = 'force from g1')\n",
    "fg2 = plt.arrow(xs[-1,0], xs[-1,1], -labmda2_dg2[0], -labmda2_dg2[1], color='r', width = .1, label = 'force from g2')\n",
    "ff = plt.arrow(xs[-1,0], xs[-1,1], -df[0], -df[1], color='k', width = .1, label = 'force from f')\n",
    "plt.legend([fg1, fg2, ff], ['force from g1', 'force from g2', 'force from f'])\n",
    "plt.grid()\n",
    "plt.xlabel('x1')\n",
    "plt.ylabel('x2')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}