{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cutting plane and bundle methods\n",
    "The second part of the tutorial focuses on cutting plane and bundle\n",
    "methods. We aim at resolving the following LASSO problem:\n",
    "$$\n",
    "\\min_r \\;f(r) =  \\| X r + b \\|_2^2 + \\lambda \\|r \\|_1\n",
    "$$\n",
    "with $\\lambda$ a given regularization parameter, $X$ and $b$ input data.\n",
    "\n",
    "## Settings\n",
    "We import the usual packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Printf, Random\n",
    "using JuMP, OSQP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fix seed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Random.seed!(2713);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "const LB = -1e20\n",
    "const UB =  1e20\n",
    "const EPS = 1e-8\n",
    "const SOLVER = OSQP.Optimizer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first generate artificial data to study the algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "include(\"data.jl\")\n",
    "\n",
    "nVariables = 50;\n",
    "nCassures = 50;\n",
    "xMin, xMax = build_bounds(nVariables)\n",
    "X = build_X(nVariables, nCassures, xMin, xMax);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cutting plane\n",
    "The cutting plane method builds a proxy $\\underline{f}_k$ for the original\n",
    "function $f$, such that $\\underline{f}_k$ is polyhedral and is a lower approximation:\n",
    "$$\n",
    "\\underline{f}_k \\leq f \\quad \\forall k\n",
    "$$\n",
    "If we have at disposition a collection of point $x_1, \\cdots, x_k$,\n",
    "with associated subgradients $g_1, \\cdots, g_k$, the function\n",
    "$\\underline{f}_k$ writes out\n",
    "$$\n",
    "\\begin{aligned}\n",
    "f_k(x) = min_x \\;& \\theta  \\\\\n",
    "         s.t. \\quad & \\theta \\geq g_k^\\top (x - x_k) + f(x_k)\n",
    "\\end{aligned}\n",
    "$$\n",
    "Using the JuMP modeler, the master problem at iteration 0 writes as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function init_master_cutting_plane(model, X, xMin, xMax)\n",
    "    nVariables, nCassures = size(X)\n",
    "    x_ = @variable(model, xMin[i] <= x_[i in 1:nVariables] <= xMax[i])\n",
    "    α_ = @variable(model, α_ >= LB)\n",
    "    @objective(model, Min, α_)\n",
    "    return x_, α_\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At each iteration $k$, we should add a new cut inside the problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function build_cut(X_, x_)\n",
    "    nVariables, nCassures = size(X_)\n",
    "    sub_gradient_ = Base.zeros(nVariables)\n",
    "    rhs_ = 0.0\n",
    "    x_tmp = JuMP.value.(x_)\n",
    "\n",
    "    for i in 1:nVariables\n",
    "        rhs_ += x_tmp[i] / nVariables\n",
    "        sub_gradient_[i] = 1.0 / nVariables\n",
    "        for j in 1:nCassures\n",
    "            if x_tmp[i] < X_[i,j]\n",
    "                rhs_+= j * 2 * abs(x_tmp[i] - X[i, j]) / nCassures / nVariables\n",
    "                sub_gradient_[i] += -j * 2.0/nCassures/nVariables\n",
    "            elseif x_tmp[i] > X_[i,j]\n",
    "                rhs_+= j * 3 * abs(x_tmp[i] - X[i, j]) / nCassures / nVariables\n",
    "                sub_gradient_[i] += j * 3.0 /nCassures/nVariables\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    return x_tmp, rhs_, sub_gradient_\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With these two ingredients, we could define the cutting plane algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function launch_cutting_plane(X, xMin, xMax, maxit=10000)\n",
    "    master = Model(SOLVER)\n",
    "    JuMP.set_silent(master)\n",
    "    # master = Model(solver=ClpSolver(LogLevel=0));\n",
    "    x, α = init_master_cutting_plane(master, X, xMin, xMax)\n",
    "    stop = false\n",
    "    n_iter = 0\n",
    "    lb, ub = LB, UB\n",
    "\n",
    "    best_ub = ub\n",
    "\n",
    "    for n_iter in 1:maxit\n",
    "        JuMP.optimize!(master)\n",
    "        lb = JuMP.value(α)\n",
    "        x_value, rhs, sub_gradient = build_cut(X, x)\n",
    "        # TODO: code cutting plane algorithm here \n",
    "    end\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bundle algorithm\n",
    "Comparing to the cutting plane method, the bundle algorithm adds a\n",
    "a quadratic penalization.\n",
    "$$\n",
    "\\begin{aligned}\n",
    "f_k(x) = min_x \\;& \\theta + \\frac 12 \\| x - x_k \\|_2^2 \\\\\n",
    "         s.t. \\quad & \\theta \\geq g_k^\\top (x - x_k) + f(x_k)\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function update_center(var_proximal, nVariables, center)\n",
    "    JuMP.fix.(var_proximal, center)\n",
    "end\n",
    "\n",
    "function launch_proximal(X, xMin, xMax, maxit=10000)\n",
    "    master = Model(SOLVER)\n",
    "    JuMP.set_silent(master)\n",
    "    nVariables, nCassures = size(X)\n",
    "    x, α = init_master_cutting_plane(master, X, xMin, xMax)\n",
    "    # Add a proximal variable\n",
    "    var_proximal = @variable(master, var_proximal[1:nVariables])\n",
    "\n",
    "    stop = false\n",
    "    lb = -Inf\n",
    "    ub = Inf\n",
    "    # Number of iteration\n",
    "    n_iter = 0\n",
    "    # Best upper bound\n",
    "    best_ub = ub\n",
    "    prediction = 0.0\n",
    "    # Number of serious step\n",
    "    nb_ss = 0\n",
    "    # Number of null step\n",
    "    nb_ns = 0\n",
    "    # Maximum number of update\n",
    "    nb_update = 3\n",
    "    step = \"NONE\"\n",
    "    weight = 1.0\n",
    "    tol = 1e-1\n",
    "\n",
    "    # The objective writes out as a QP\n",
    "    @objective(master, Min, α + sum(weight * var_proximal[i]^2 for i in 1:nVariables))\n",
    "\n",
    "    center = zeros(nVariables)\n",
    "    update_center(var_proximal, nVariables, center)\n",
    "\n",
    "    for n_iter in 1:maxit\n",
    "        JuMP.optimize!(master)\n",
    "        lb = JuMP.value(α)\n",
    "        x_value, ub, sub_gradient = build_cut(X, x)\n",
    "        # TODO: code bundle algorithm here\n",
    "    end\n",
    "end"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.2.0",
   "language": "julia",
   "name": "julia-1.2"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}