{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Plotting-central-path---of-an-LP\" data-toc-modified-id=\"Plotting-central-path---of-an-LP-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Plotting central path - of an LP</a></span></li><li><span><a href=\"#Solving-by-IPM\" data-toc-modified-id=\"Solving-by-IPM-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Solving by IPM</a></span></li></ul></div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# uncomment and run once\n",
    "# Using Pkg\n",
    "#Pkg.add(\"Plots\")\n",
    "#Pkg.add(\"LinearAlgebra\")\n",
    "#Pkg.add(\"OPtim\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Pkg.add(\"PlotlyJS\")\n",
    "#Pkg.add(\"ORCA\")# if you want to try plotlyjs()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots, LinearAlgebra, Optim \n",
    "plotlyjs() #Plots has multiple plotting backend. plotlyjs() has more interactivity in notebook but sometimes buggy\n",
    "#pyplot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = x -> (1-x[1])^2 + 5*(x[2] - x[1]^2)^2\n",
    "df = x -> [2*(10*x[1]^3-10*x[1]*x[2]+x[1]-1), 10*(x[2]-x[1]^2)]\n",
    "hf = x -> [60x[1]^2-20x[2]+2 -20x[1]; -20x[1] 10]\n",
    "\n",
    "xdomain = -2:0.01:2\n",
    "ydomain = -2:0.01:2\n",
    "\n",
    "lev_line = contour(xdomain,ydomain,(x,y)->f([x,y]), levels = 200);\n",
    "lev_line = plot!(lev_line,[1],[1],seriestype = :scatter, color = :red,label=\"solution\");\n",
    "lev_line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = Float64[-0.5, 2]\n",
    "α0 = 0.02\n",
    "IT_MAX = 5000\n",
    "ε = 1e-6\n",
    "\n",
    "function alg_grad(x0,f,df;α0=α0)\n",
    "    x = x0\n",
    "    \n",
    "     pts = [x]\n",
    "    vals = [f(x)]\n",
    "    grads = [df(x)]\n",
    "\n",
    "    stop_test = false\n",
    "    it = 0 #iteration number\n",
    "\n",
    "    while !stop_test\n",
    "        g = - df(x)       # descent direction\n",
    "        α = α0            # step size\n",
    "        x += α*g          # step\n",
    "\n",
    "        # keeping track of results\n",
    "        push!(pts,x)\n",
    "        push!(vals,f(x))\n",
    "        push!(grads,g)\n",
    "        \n",
    "        # stopping test\n",
    "        it += 1\n",
    "        stop_test = (it > IT_MAX - 1) || (norm(g)<ε)\n",
    "    end\n",
    "    \n",
    "    println(\"gradient algorithm with fixed step α=\",α0, \" ended in \", it, \" iterations\")\n",
    "    return pts, vals, grads\n",
    "   \n",
    "end\n",
    "\n",
    "pts, vals, grads = alg_grad(x0,f,df,α0=0.02);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function plot_trajectory(pts;color=\"blue\",label=`AUTO`,plt=lev_line)\n",
    "    plt = deepcopy(plt)\n",
    "    plot!(plt,[p[1] for p in pts],[p[2] for p in pts],color=color,label=label,shape=:cross)\n",
    "    return plt\n",
    "end\n",
    "\n",
    "plt = plot_trajectory(pts,label=\"gradient\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a very simple linear search guaranteeing that value is indeed decreasing\n",
    "# if you set it_max = 1 you only keep the step = 1. \n",
    "# You can do this to explore why it is not always a good idea\n",
    "function reducing_step(f,x,d;reducing_factor=.9,it_max = 10, verbose = true)\n",
    "    f0 = f(x)\n",
    "    α = 1\n",
    "    it = 0\n",
    "    while (f(x+α*d) > f(x) + 1e-8) && (it <= it_max)\n",
    "        α = α*reducing_factor\n",
    "        it += 1\n",
    "        verbose && println(\"step reduced, it = \", it) \n",
    "    end\n",
    "    return α\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function alg_Newton(x0,f,df, hf;reduced_step=false)\n",
    "    x = x0\n",
    "    \n",
    "    pts = [x]\n",
    "    vals = [f(x)]\n",
    "    grads = [df(x)]\n",
    "\n",
    "    stop_test = false\n",
    "    it = 0 #iteration number\n",
    "\n",
    "    while !stop_test\n",
    "        g = df(x)         # gradient\n",
    "        H = hf(x)         # hessian\n",
    "        d = - H \\ g       # descent direction \n",
    "        reduced_step ? α=reducing_step(f,x,d) : α = 1   # step size\n",
    "        x += α*d          # step\n",
    "\n",
    "        # keeping track of results\n",
    "        push!(pts,x)\n",
    "        push!(vals,f(x))\n",
    "        push!(grads,g)\n",
    "        \n",
    "        # stopping test\n",
    "        it += 1\n",
    "        stop_test = (it > IT_MAX - 1) || (norm(g)<ε)\n",
    "    end\n",
    "    \n",
    "    println(\"Newton algorithm ended in \", it, \" iterations\")\n",
    "    return pts, vals, grads\n",
    "   \n",
    "end\n",
    "\n",
    "pts,vals,grads = alg_Newton(x0,f,df, hf)\n",
    "pts2,vals2,grads2 = alg_Newton(x0,f,df, hf,reduced_step=true)\n",
    "\n",
    "plt=plot_trajectory(pts,label=\"Newton\",color=\"green\",plt = plt)\n",
    "plt=plot_trajectory(pts2,label=\"Newton_reduced\",color=\"violet\",plt = plt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting central path - of an LP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Optim, LinearAlgebra, Plots\n",
    "\n",
    "# constructing the LP min {c'x | Ax <= b}\n",
    "A = [-1 0; 0 -1 ; 0.5 1; 1 -1]\n",
    "b = [0 ; 0 ; 1; 0.5]\n",
    "\n",
    "c = [0.1, 1]\n",
    "\n",
    "\n",
    "m,n = size(A)\n",
    "\n",
    "function ln(x) #extended ln to avoid undefinite value\n",
    "    x >= 10^(-6) ? (return log(x)) : (return -10 + 10^8*x)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Logarithmic penalty\n",
    "ϕ(x) =  - sum([ln(b[i]-A[i,:]'x) for i in 1:m])\n",
    "\n",
    "function gradient_phi(x)\n",
    "    grad = zeros(n)\n",
    "    # A REMPLIR \n",
    "    for i in 1:m\n",
    "        grad += 1/(b[i]-A[i,:]'x)*A[i,:]\n",
    "    end\n",
    "    return grad\n",
    "end\n",
    "gradient_phi(x0)\n",
    "\n",
    "function finite_diff(x,δ)\n",
    "  println((ϕ(x + [δ, 0]) - ϕ(x - [δ, 0]))/2δ)\n",
    "  println((ϕ(x + [0, δ]) - ϕ(x - [0, δ]))/2δ)\n",
    "end\n",
    "\n",
    "finite_diff(x0,0.001) # checking gradient by finite difference\n",
    "gradient_phi(x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function central_point(t,c)\n",
    "    f(x) = t*c'*x  + ϕ(x)\n",
    "    x0 = [0.4,0.4]\n",
    "    res = optimize(f, x0)\n",
    "    return Optim.minimizer(res) # using an external optimization library\n",
    "end\n",
    "\n",
    "function central_path(t_list,c)\n",
    "    path = []\n",
    "    for t in t_list\n",
    "        x = central_point(t,c)\n",
    "        push!(path,x)\n",
    "        #println(x[1],' ',x[2])\n",
    "    end\n",
    "    return path\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plotting tool\n",
    "function plot_trajectory(pts;color=\"blue\",label=\"\")\n",
    "    plt = plot(Shape([0, 0.5, 1,0],[0, 0., 0.5,1]),opacity=.5,label=\"\")\n",
    "    plot!(plt,[p[1] for p in pts],[p[2] for p in pts],color=color,label=label,shape=:cross)\n",
    "    return plt\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plotting the central path\n",
    "t0 = 1\n",
    "\n",
    "c_path = central_path([t0*1.1^t for t in 0:50],c)\n",
    "plot_trajectory(c_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving by IPM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function Newton_direction(x,t,c) # compute Newton direction from x for t*c'x+ϕ(x)\n",
    "    # A REMPLIR\n",
    "    d = [1/(b[i]-A[i,:]'x) for i in 1:m]\n",
    "    G = t*c + A'd \n",
    "    H = A'Diagonal(d)^2 *A\n",
    "    return - H \\ G\n",
    "end\n",
    "\n",
    "function Newton(x,t,c,n=5)\n",
    "    f(x) = t*c'x + ϕ(x)\n",
    "    for _ in 1:n\n",
    "        d = Newton_direction(x,t,c)\n",
    "        α = reducing_step(f,x,d)\n",
    "        x += α*Newton_direction(x,t,c) \n",
    "        #println(\"(\",x[1],\",\",x[2],\")\")\n",
    "    end\n",
    "    return x\n",
    "end\n",
    "\n",
    "# Checking that Newton's algorithm is converging\n",
    "Newton(x0,t0,c,10) - central_point(t0,c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function IPM(x0,t0,ρ,c; N_outer = 5, N_inner = 3,compute_central_points=false)\n",
    "    xout_list = [x0] # list of initial points of each outer iteration\n",
    "    x_list =[x0] # list of all the points\n",
    "    cp_list =[] # list of the central points associated to t for each outer iteration\n",
    "    x = x0\n",
    "    t = t0\n",
    "    \n",
    "    # A REMPLIR\n",
    "    for n_outer in 1:N_outer\n",
    "        f(x) = t*c'x + ϕ(x)\n",
    "        for n_inner in 1:N_inner\n",
    "            d = Newton_direction(x,t,c)\n",
    "            α = reducing_step(f,x,d)\n",
    "            x += α*Newton_direction(x,t,c)\n",
    "            push!(x_list,x)\n",
    "        end\n",
    "        push!(xout_list,x)\n",
    "        compute_central_points && push!(cp_list,central_point(t,c))\n",
    "        t = ρ*t\n",
    "         \n",
    "    end\n",
    "    \n",
    "    # if compute_central_points is true return 3 list otherwise 2 \n",
    "    compute_central_points ? (return x_list,xout_list,cp_list) : (return x_list,xout_list)\n",
    "    \n",
    "end\n",
    "\n",
    "x_list,xout_list = IPM(x0,t0,ρ,c; N_outer = 5, N_inner = 3)\n",
    "x_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = [0.1,1] # cost\n",
    "x0 = [0.5,0.5]\n",
    "\n",
    "#### Algorithm parameters\n",
    "ρ = 2. # multiplier of t between \n",
    "t0 = 1 # initial t\n",
    "\n",
    "N_outer = 10 # number of outer iterations\n",
    "N_inner = 3  # number of inner iterations\n",
    "\n",
    "# recompute central path if you change c = [0.1,1] \n",
    "#c_path = central_path([t0*1.1^t for t in 0:50],c)\n",
    "\n",
    "\n",
    "\n",
    "function plot_IPM(x0,t0,ρ,c; N_outer = 5, N_inner = 5,plot_central_points=true)\n",
    "    \n",
    "    #compute trajectory\n",
    "    if plot_central_points\n",
    "        x_list,xout_list, cp_list = IPM(x0,t0,ρ,c, N_outer = N_outer, N_inner = N_inner,compute_central_points=true)\n",
    "    else\n",
    "        x_list,xout_list = IPM(x0,t0,ρ,c, N_outer = N_outer, N_inner = N_inner)\n",
    "    end\n",
    "    \n",
    "    # plot polyhedron\n",
    "    plt = plot(Shape([0, 0.5, 1,0],[0, 0., 0.5,1]),opacity=.5,label=\"\",size=(600,400))\n",
    "    \n",
    "    # plot -c (scaled)\n",
    "    plot!(plt,[x0[1],x0[1]-0.1*c[1]],[x0[2],x0[2]-0.1*c[2]],arrow=2, label=\"-c\",lw=2)\n",
    "    \n",
    "    #plot central path\n",
    "    plot!(plt,[p[1] for p in c_path],[p[2] for p in c_path],color=\"blue\",label=\"central path\",lw=2)\n",
    "    \n",
    "    # plot IPM trajectory\n",
    "    plot!(plt,[p[1] for p in x_list],[p[2] for p in x_list],color=\"orange\",label=\"IPM\",shape=:cross,lw=2)\n",
    "    # overline outer points\n",
    "    plot!(plt,[p[1] for p in xout_list],[p[2] for p in xout_list],\n",
    "        seriestype = :scatter,color=\"red\", markersize=6,label=\"\",shape=:cross)\n",
    "    \n",
    "    #plot central points\n",
    "    if plot_central_points\n",
    "        plot!(plt,[p[1] for p in cp_list],[p[2] for p in cp_list],\n",
    "            seriestype = :scatter,color=\"violet\", markersize=4,label=\"central points\",shape=:diamond)\n",
    "    end\n",
    "    return plt\n",
    "end\n",
    "\n",
    "plt = plot_IPM(x0,t0,ρ,c, N_outer =N_outer, N_inner = N_inner);\n",
    "gui(plt)"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.4.0",
   "language": "julia",
   "name": "julia-1.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.4.0"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}