Walking_in_the_shadow_code.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import rc\n",
    "from scipy.optimize import minimize_scalar, minimize\n",
    "import datetime\n",
    "import os\n",
    "import networkx as nx\n",
    "import pandas as pd \n",
    "from gurobipy import*\n",
    "from collections import Counter, defaultdict\n",
    "from scipy.io import loadmat\n",
    "from time import process_time\n",
    "from matplotlib import ticker\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "from autograd import grad as grad_a\n",
    "from sklearn.preprocessing import StandardScaler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Oracles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def line_search(x, d, gamma_max,func):\n",
    "    \n",
    "    #line-search using Brent's rule in scipy\n",
    "    \n",
    "    '''\n",
    "    Minimizes f over [x, y], i.e., f(x+gamma*d) as a function of scalar gamma in [0,gamma_max]\n",
    "    '''\n",
    "\n",
    "    def fun(gamma):\n",
    "        ls = x + gamma*d\n",
    "        return func(ls)\n",
    "\n",
    "    res = minimize_scalar(fun, bounds=(0, gamma_max), method='bounded')\n",
    "\n",
    "    gamma = res.x\n",
    "    ls = x + gamma*d        \n",
    "    return ls, gamma\n",
    "\n",
    "\n",
    "def segment_search(f, grad_f, x, y, tol=1e-6, stepsize=True):\n",
    "    \n",
    "    #line-search using golden-section rule coded from scratch\n",
    "    \n",
    "    '''\n",
    "    Minimizes f over [x, y], i.e., f(x+gamma*(y-x)) as a function of scalar gamma in [0,1]\n",
    "    '''\n",
    "    \n",
    "    # restrict segment of search to [x, y]\n",
    "    d = (y-x).copy()\n",
    "    left, right = x.copy(), y.copy()\n",
    "    \n",
    "    # if the minimum is at an endpoint\n",
    "    if np.dot(d, grad_f(x))*np.dot(d, grad_f(y)) >= 0:\n",
    "        if f(y) <= f(x):\n",
    "            return y, 1\n",
    "        else:\n",
    "            return x, 0\n",
    "    \n",
    "    # apply golden-section method to segment\n",
    "    gold = (1+np.sqrt(5))/2\n",
    "    improv = np.inf\n",
    "    while improv > tol:\n",
    "        old_left, old_right = left, right\n",
    "        new = left+(right-left)/(1+gold)\n",
    "        probe = new+(right-new)/2\n",
    "        if f(probe) <= f(new):\n",
    "            left, right = new, right\n",
    "        else:\n",
    "            left, right = left, probe\n",
    "        improv = np.linalg.norm(f(right)-f(old_right))+np.linalg.norm(f(left)-f(old_left))\n",
    "    x_min = (left+right)/2\n",
    "    \n",
    "    # compute step size gamma\n",
    "    gamma = 0\n",
    "    if stepsize == True:\n",
    "        for i in range(len(d)):\n",
    "            if d[i] != 0:\n",
    "                gamma = (x_min[i]-x[i])/d[i]\n",
    "                break\n",
    "                \n",
    "    return x_min, gamma\n",
    "\n",
    "\n",
    "\n",
    "def active_constraints_oracle(A,b,x):\n",
    "    \n",
    "    '''\n",
    "    finding active constraints for polytope of the form Ax <= b\n",
    "    '''\n",
    "    b_prime = []\n",
    "    A_prime = []\n",
    "    A_not = []\n",
    "    b_not = []\n",
    "    for i in range(len(A)):\n",
    "        if np.abs(np.round(np.dot(A[i],x),4) - b[i]) < 0.01:\n",
    "            A_prime.append(list(A[i]))\n",
    "            b_prime.append(b[i])\n",
    "        else:\n",
    "            A_not.append(list(A[i]))\n",
    "            b_not.append(b[i])\n",
    "            \n",
    "    return np.array(A_prime),np.array(b_prime),np.array(A_not),np.array(b_not)\n",
    "\n",
    "\n",
    "def max_step_size(A,b,x, d):\n",
    "    \n",
    "    \n",
    "    '''\n",
    "    finding maximum step-size: argmax{\\delta | x + \\delta d \\in P} where P is given by Ax <= b and d is a \n",
    "    feasible direction at x. Here we use Gurobi as a black box solver.\n",
    "    '''\n",
    "\n",
    "    m = Model(\"opt\")\n",
    "    n = len(A.T)\n",
    "\n",
    "    lam = m.addVar(lb=-GRB.INFINITY, name='lam')\n",
    "\n",
    "    m.update()              \n",
    "\n",
    "    objExp = lam\n",
    "    m.setObjective(objExp, GRB.MAXIMIZE)\n",
    "    \n",
    "    m.update()\n",
    "\n",
    "    #feasibility constraints\n",
    "    \n",
    "    for i in range(len(A)):\n",
    "        m.addConstr(np.dot(np.array(x),A[i]) + lam* np.dot(np.array(d),A[i]),'<=', b[i])\n",
    "        \n",
    "    m.update()\n",
    "\n",
    "    #optimize\n",
    "    m.setParam( 'OutputFlag', False )\n",
    "    m.optimize()\n",
    "    \n",
    "    return lam.x\n",
    "    \n",
    "\n",
    "def max_step_size_search(A,b,x, d):\n",
    "    \n",
    "    '''\n",
    "    finding maximum step-size: argmax{\\delta | x + \\delta d \\in P} where P is given by Ax <= b and d is a \n",
    "    feasible direction at x. Here we search over set of inactive constraints that might get violated.\n",
    "    '''\n",
    "\n",
    "    A_J,b_J = active_constraints_oracle(A,b,x)[2:]\n",
    "    slack = b_J - np.dot(A_J,x)\n",
    "    \n",
    "    try:\n",
    "        if len(d) > 1:\n",
    "            denom = np.dot(A_J,d_pi)\n",
    "            excess = [slack[i]/denom[i] for i in range(len(slack)) if denom[i] > 0]\n",
    "            gamma_max = min(excess)\n",
    "            return gamma_max\n",
    "    except:\n",
    "        return min(slack)\n",
    "\n",
    "\n",
    "def shadow_oracle(A,b,x, d):\n",
    "     \n",
    "    '''\n",
    "    Projecting d on set of feasible directions for polytope given by Ax <= b\n",
    "    '''\n",
    "    \n",
    "    #find matrix of active constraints\n",
    "    A_I = active_constraints_oracle(A,b,x)[0]\n",
    "\n",
    "    m = Model(\"opt\")\n",
    "    n = len(A_I.T)\n",
    "\n",
    "    z = []\n",
    "    for i in range(n):\n",
    "        z.append(m.addVar(lb=-GRB.INFINITY, name='z_{}'.format(i)))\n",
    "\n",
    "    m.update()              \n",
    "\n",
    "    objExp = quicksum(np.square(np.array([d[k] -z[k] for k in range(n)])))\n",
    "   \n",
    "    m.setObjective(objExp, GRB.MINIMIZE)\n",
    "    m.update()\n",
    "\n",
    "    #feasibility constraints\n",
    "    for i in range(len(A_I)):\n",
    "        m.addConstr(np.dot(np.array(z),A_I[i]),'<=', 0)\n",
    "        \n",
    "    m.update()\n",
    "\n",
    "    #optimize\n",
    "    m.setParam( 'OutputFlag', False )\n",
    "    m.optimize()\n",
    "\n",
    "    der = [i.x for i in z]\n",
    "\n",
    "    return np.array(der)\n",
    "\n",
    "\n",
    "def shadow_oracle_other(A,b,x, grad,M,vert_rep,vertices):\n",
    "     \n",
    "    '''\n",
    "    Computing projection if we know an upper bound M on the value at which normal cone is not changing\n",
    "    '''\n",
    "    \n",
    "    epsilon = M/2\n",
    "    \n",
    "    if vert_rep == True:\n",
    "        g = projection_oracle_vertices(vertices,x - epsilon*grad)\n",
    "        der = (g -x)/epsilon\n",
    "        \n",
    "    else:\n",
    "        projection_oracle(x - epsilon*grad,A,b)\n",
    "        der = (g -x)/epsilon\n",
    "\n",
    "    return der\n",
    "\n",
    "    return np.array(der)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Algorithms - FW Variants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FW"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def FW(x, lmo, epsilon,func,grad_f, f_tol, time_tol):\n",
    "    \n",
    "    #record primal gap, function value, and time every iteration\n",
    "    now=datetime.datetime.now()\n",
    "    primal_gap = []\n",
    "    function_value=[func(x)]\n",
    "    time = [0]\n",
    "    f_improv = np.inf\n",
    "\n",
    "    #initialize iteration count\n",
    "    t = 0    \n",
    "\n",
    "    while f_improv > f_tol and time[-1] < time_tol:\n",
    "        \n",
    "        start = process_time()\n",
    "        \n",
    "        #compute gradient\n",
    "        grad = grad_f(x)\n",
    "\n",
    "        #solve linear subproblem and compute FW direction\n",
    "        v = lmo(grad)\n",
    "        d_FW = v-x\n",
    "\n",
    "        #If primal gap is small enough - terminate\n",
    "        if np.dot(-grad,d_FW) <= epsilon:\n",
    "            break\n",
    "        else:\n",
    "            #update convergence data\n",
    "            primal_gap.append(np.dot(-grad,d_FW))\n",
    "\n",
    "        #Update next iterate by doing a feasible line-search\n",
    "        x, gamma = segment_search(func, grad_f, x, v)\n",
    "        \n",
    "        end = process_time()\n",
    "        \n",
    "        time.append(time[t] + end - start)\n",
    "        f_improv = function_value[-1] - func(x)\n",
    "        function_value.append(func(x))\n",
    "        \n",
    "        t+=1\n",
    "        \n",
    "    return x, function_value, time,t,primal_gap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# AFW "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Function to compute away vertex\n",
    "def away_step(grad, S):\n",
    "    \n",
    "    '''\n",
    "    Compute away vertex by searching over current active set\n",
    "    '''\n",
    "    \n",
    "    costs = {}\n",
    "    \n",
    "    for k,v in S.items():\n",
    "        cost = np.dot(k,grad)\n",
    "        costs[cost] = [k,v]\n",
    "    vertex, alpha = costs[max(costs.keys())]  \n",
    "    return vertex,alpha\n",
    "\n",
    "#Function to update active set\n",
    "def update_S(S,gamma, Away, vertex):\n",
    "    \n",
    "    '''\n",
    "    Update convex decompistion of active step after every iteration \n",
    "    '''\n",
    "    \n",
    "    S = S.copy()\n",
    "    vertex = tuple(vertex)\n",
    "    \n",
    "    if not Away:\n",
    "        if vertex not in S.keys():\n",
    "            S[vertex] = gamma\n",
    "        else:\n",
    "            S[vertex] *= (1-gamma)\n",
    "            S[vertex] += gamma\n",
    "            \n",
    "        for k in S.keys():\n",
    "            if k != vertex:\n",
    "                S[k] *= (1-gamma)\n",
    "    else:\n",
    "        for k in S.keys():\n",
    "            if k != vertex:\n",
    "                S[k] *= (1+gamma)\n",
    "            else:\n",
    "                S[k] *= (1+gamma)\n",
    "                S[k] -= gamma\n",
    "    return {k:v for k,v in S.items() if np.round(v,3) > 0}\n",
    "\n",
    "\n",
    "#AFW Algorithm\n",
    "def AFW(x, lmo, epsilon,func,grad_f, f_tol, time_tol):\n",
    "    \n",
    "    #record primal gap, function value, and time every iteration\n",
    "    now=datetime.datetime.now()\n",
    "    primal_gap = []\n",
    "    function_value=[func(x)]\n",
    "    time = [0]\n",
    "    f_improv = np.inf\n",
    "\n",
    "    #initialize starting point and active set\n",
    "    t = 0    \n",
    "    S={tuple(x): 1}\n",
    "\n",
    "    while f_improv > f_tol and time[-1] < time_tol:\n",
    "        \n",
    "        start = process_time()\n",
    "        \n",
    "        #compute gradient\n",
    "        grad = grad_f(x)\n",
    "\n",
    "        #solve linear subproblem and compute FW direction\n",
    "        v = lmo(grad)\n",
    "        d_FW = v-x\n",
    "\n",
    "        #If primal gap is small enough - terminate\n",
    "        if np.dot(-grad,d_FW) <= epsilon:\n",
    "            break\n",
    "        else:\n",
    "            #update convergence data\n",
    "            primal_gap.append(np.dot(-grad,d_FW))\n",
    "\n",
    "        #Compute away vertex and direction\n",
    "        a,alpha_a = away_step(grad, S)\n",
    "        d_A = x - a\n",
    "\n",
    "        #Check if FW gap is greater than away gap\n",
    "        if np.dot(-grad,d_FW) >= np.dot(-grad,d_A):\n",
    "            #choose FW direction\n",
    "            d = d_FW\n",
    "            vertex = v\n",
    "            gamma_max = 1\n",
    "            Away = False\n",
    "        else:\n",
    "            #choose Away direction\n",
    "            d = d_A\n",
    "            vertex = a\n",
    "            gamma_max = alpha_a/(1-alpha_a)\n",
    "            Away = True\n",
    "\n",
    "        #Update next iterate by doing a feasible line-search\n",
    "        x, gamma = segment_search(func, grad_f, x, x + gamma_max *d)\n",
    "\n",
    "        #update active set based on direction chosen\n",
    "        S = update_S(S,gamma, Away, vertex)\n",
    "        \n",
    "        end = process_time()\n",
    "        time.append(time[t] + end - start)\n",
    "        \n",
    "        f_improv = function_value[-1] - func(x)\n",
    "        function_value.append(func(x))\n",
    "        \n",
    "        t+=1\n",
    "        \n",
    "    return x, function_value, time,t,primal_gap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pairwise FW"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_S_PW(S,gamma, FW_vertex, Away_vertex):\n",
    "    \n",
    "    '''\n",
    "    Update convex decompistion of active step after pairwise direction is chosen \n",
    "    '''\n",
    "    \n",
    "    vertex1 = tuple(FW_vertex)\n",
    "    vertex2 = tuple(Away_vertex)\n",
    "    S = S.copy()\n",
    "    \n",
    "    if vertex1 not in S.keys():\n",
    "        S[vertex1] = gamma\n",
    "    else:\n",
    "        S[vertex1] += gamma\n",
    "    \n",
    "    S[vertex2] -= gamma\n",
    "    return {k:v for k,v in S.items() if np.round(v,4) > 0}\n",
    "\n",
    "\n",
    "#PFW Algorithm\n",
    "def PFW(x, lmo, epsilon,func,grad_f, f_tol, time_tol):\n",
    "    \n",
    "    #record primal gap, function value, and time every iteration\n",
    "    now=datetime.datetime.now()\n",
    "    primal_gap = []\n",
    "    function_value=[func(x)]\n",
    "    time = [0]\n",
    "    f_improv = np.inf\n",
    "\n",
    "    #initialize starting point and active set\n",
    "    t = 0    \n",
    "    S={tuple(x): 1}\n",
    "\n",
    "    while f_improv > f_tol and time[-1] < time_tol:\n",
    "        \n",
    "        start = process_time()\n",
    "        \n",
    "        #compute gradient\n",
    "        grad = grad_f(x)\n",
    "\n",
    "        #solve linear subproblem and compute FW direction\n",
    "        v = lmo(grad)\n",
    "        d_FW = v-x\n",
    "\n",
    "        #If primal gap is small enough - terminate\n",
    "        if np.dot(-grad,d_FW) <= epsilon:\n",
    "            break\n",
    "        else:\n",
    "            #update convergence data\n",
    "            primal_gap.append(np.dot(-grad,d_FW))\n",
    "\n",
    "        #Compute away vertex and direction\n",
    "        a,alpha_a = away_step(grad, S)\n",
    "        d_A = x - a\n",
    "\n",
    "        #Pairwise step\n",
    "        d =  d_FW + d_A\n",
    "        gamma_max = alpha_a\n",
    "\n",
    "        #Update next iterate by doing a feasible line-search\n",
    "        x, gamma = line_search(x, d, gamma_max,func)\n",
    "\n",
    "        #update active set\n",
    "        S = update_S_PW(S,gamma, v, a)\n",
    "        \n",
    "        end = process_time()\n",
    "        time.append(time[t] + end - start)\n",
    "        f_improv = function_value[-1] - func(x)\n",
    "        function_value.append(func(x))\n",
    "        \n",
    "        t+=1\n",
    "        \n",
    "    return x, function_value, time,t,primal_gap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Decoposition Invaraint FW (DICG)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#DICG algorithm\n",
    "def DICG(x, lmo, feasibility_oracle, epsilon,func,grad_f, f_tol, time_tol):\n",
    "    \n",
    "    #record primal gap, function value, and time every iteration\n",
    "    now=datetime.datetime.now()\n",
    "    primal_gap = []\n",
    "    function_value=[func(x)]\n",
    "    time = [0]\n",
    "    f_improv = np.inf\n",
    "\n",
    "    #initialize starting point and active set\n",
    "    t = 0    \n",
    "\n",
    "    while f_improv > f_tol and time[-1] < time_tol:\n",
    "        \n",
    "        start = process_time()\n",
    "        \n",
    "        #compute gradient\n",
    "        grad = grad_f(x)\n",
    "\n",
    "        #solve linear subproblem and compute FW direction\n",
    "        v = lmo(grad)\n",
    "        d_FW = v-x\n",
    "\n",
    "        #If primal gap is small enough - terminate\n",
    "        if np.dot(-grad,d_FW) <= epsilon:\n",
    "            break\n",
    "        else:\n",
    "            #update convergence data\n",
    "            primal_gap.append(np.dot(-grad,d_FW))\n",
    "\n",
    "        #create new gradient to find best away vertex\n",
    "        g = np.array([grad[i] if x[i] > 0 else -9e9 for i in range(len(x))])\n",
    "        a = lmo(-g)\n",
    "        d_A = x - a\n",
    "\n",
    "        #Pairwise step\n",
    "        d =  d_FW + d_A\n",
    "        gamma_max = feasibility_oracle(x,d)\n",
    "\n",
    "        #Update next iterate by doing a feasible line-search\n",
    "        x, gamma = segment_search(func, grad_f, x, x+gamma_max*d)\n",
    "        \n",
    "        end = process_time()\n",
    "        time.append(time[t] + end - start)\n",
    "        f_improv = function_value[-1] - func(x)\n",
    "        function_value.append(func(x))\n",
    "        \n",
    "        t+=1\n",
    "        \n",
    "    return x, function_value, time,t,primal_gap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Projected Gradient variants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Shadow Walk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trace_PW_curve(x,grad,shadow,feasibility_oracle,func,tol):\n",
    "    \n",
    "    '''\n",
    "    trace the piecewise linear projection curve\n",
    "    '''\n",
    "    \n",
    "    count = 0\n",
    "    t = 0\n",
    "    \n",
    "    while True:\n",
    "        #compute shadow of gradient\n",
    "        d_pi = shadow(x,grad)\n",
    "        t1 = process_time()\n",
    "        \n",
    "        #find maximum step size in which one can move along shadow\n",
    "        gamma_max = feasibility_oracle(x, d_pi)\n",
    "        \n",
    "        #do line search alog shadow\n",
    "        x, gamma = segment_search(func, grad_f, x, x + gamma_max *d_pi)\n",
    "        \n",
    "        #if optimal step is not maximal or we reached endpoint of curve - terminate\n",
    "        if abs(gamma - 1) > 0.001 or np.dot(d_pi,d_pi)**0.5 < tol:\n",
    "            break\n",
    "        else:\n",
    "            count +=1 #record breakpoints\n",
    "        t2 = process_time()\n",
    "        t += (t2 - t1)\n",
    "\n",
    "    return x,count,t\n",
    "\n",
    "\n",
    "#shadow descent algorithm:\n",
    "def SD(x, lmo, shadow,feasibility_oracle, epsilon,func,grad_f, f_tol, time_tol,tol):\n",
    "    \n",
    "    #record primal gap, function value, and time every iteration\n",
    "    now=datetime.datetime.now()\n",
    "    primal_gap = []\n",
    "    function_value=[func(x)]\n",
    "    time1 = [0]\n",
    "    time2 = [0]\n",
    "    f_improv = np.inf\n",
    "    counts = []\n",
    "\n",
    "    #initialize starting point and active set\n",
    "    t = 0    \n",
    "\n",
    "    while f_improv > f_tol and time1[-1] < time_tol:\n",
    "        \n",
    "        start = process_time()\n",
    "        \n",
    "        #compute gradient\n",
    "        grad = grad_f(x)\n",
    "\n",
    "        #solve linear subproblem and compute FW direction\n",
    "        s = lmo(grad)\n",
    "        d_FW = s-x\n",
    "\n",
    "        #If primal gap is small enough - terminate\n",
    "        if np.dot(-grad,d_FW) <= epsilon:\n",
    "            break\n",
    "        else:\n",
    "            #update convergence data\n",
    "            primal_gap.append(np.dot(-grad,d_FW))\n",
    "        \n",
    "        #Compute directional derivative\n",
    "        t1 = process_time()\n",
    "        d_pi = shadow(x,grad)\n",
    "        t2 = process_time()\n",
    "        \n",
    "        gamma_max = feasibility_oracle(x, d_pi)\n",
    "\n",
    "        #Update next iterate by doing a feasible line-search\n",
    "        y, gamma = segment_search(func, grad_f, x, x + gamma_max *d_pi)\n",
    "        \n",
    "        #check for boundary case\n",
    "        if abs(gamma - 1) < 0.001:\n",
    "            \n",
    "            trace = True\n",
    "            \n",
    "            #trace PW curve\n",
    "            t3 = process_time()\n",
    "            x, c, t_pw = trace_PW_curve(x,grad,shadow,feasibility_oracle,func,tol)\n",
    "            counts.append(c)\n",
    "            t4 = process_time()\n",
    "            \n",
    "        else:\n",
    "            trace = False\n",
    "            x = y\n",
    "            counts.append(0)\n",
    "\n",
    "        end = process_time()\n",
    "        time1.append(time1[t] + end - start)\n",
    "\n",
    "        if trace:\n",
    "            time2.append( time2[t] + end - start - (t2 -t1) -(t4 - t3) + t_pw)\n",
    "        else:\n",
    "            time2.append( time2[t] + end - start - (t2 -t1))\n",
    "            \n",
    "        f_improv = function_value[-1] - func(x)\n",
    "        function_value.append(func(x))\n",
    "        t+=1\n",
    "        \n",
    "    return x, function_value, time1,time2,t,primal_gap,counts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Shadow CG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#shadow conditional gradients\n",
    "def SCG(x, lmo, shadow,feasibility_oracle, epsilon,func,grad_f, f_tol, time_tol,tol):\n",
    "    \n",
    "    #record primal gap, function value, and time every iteration\n",
    "    now=datetime.datetime.now()\n",
    "    primal_gap = []\n",
    "    function_value=[func(x)]\n",
    "    time1 = [0]\n",
    "    time2 = [0]\n",
    "    f_improv = np.inf\n",
    "    counts = []\n",
    "    FW = []\n",
    "\n",
    "    #initialize starting point and active set\n",
    "    t = 0    \n",
    "\n",
    "    while f_improv > f_tol and time1[-1] < time_tol:\n",
    "        \n",
    "        start = process_time()\n",
    "        #compute gradient\n",
    "        grad = grad_f(x)\n",
    "\n",
    "        #solve linear subproblem and compute FW direction\n",
    "        v = lmo(grad)\n",
    "        d_FW = v-x\n",
    "\n",
    "        #If primal gap is small enough - terminate\n",
    "        if np.dot(-grad,d_FW) <= epsilon:\n",
    "            break\n",
    "        else:\n",
    "            #update convergence data\n",
    "            primal_gap.append(np.dot(-grad,d_FW))\n",
    "            \n",
    "        #Compute directional derivative\n",
    "        t1 = process_time()\n",
    "        d_pi = shadow(x,grad)\n",
    "        t2 = process_time()\n",
    "        gamma_max = feasibility_oracle(x, d_pi)\n",
    "        \n",
    "        #Check if FW gap direction is better than normalized shadow\n",
    "        if np.dot(-grad,d_FW) >= np.dot(-grad,d_pi/(np.dot(d_pi,d_pi)**0.5)):\n",
    "            \n",
    "            #wrap around using FW direction\n",
    "            d = d_FW\n",
    "            vertex = v\n",
    "            gamma_max = 1\n",
    "            \n",
    "            #record which direction is chosen \n",
    "            shadow_dir = False\n",
    "            FW.append(1) \n",
    "            \n",
    "        else:\n",
    "            \n",
    "            #choose shadow direction\n",
    "            d = d_pi\n",
    "            \n",
    "            #record which direction is chosen \n",
    "            shadow_dir = True\n",
    "            FW.append(0) \n",
    "\n",
    "        #Update next iterate by doing a feasible line-search\n",
    "        y, gamma = segment_search(func, grad_f, x, x + gamma_max*d)\n",
    "        \n",
    "        #check for boundary case in shadow steps\n",
    "        if shadow_dir and abs(gamma - 1) < 0.001:\n",
    "        \n",
    "            #trace PW curve\n",
    "            trace = True\n",
    "            t3 = process_time()\n",
    "            x, c, t_pw = trace_PW_curve(x,grad,shadow,feasibility_oracle,func,tol)\n",
    "            t4 = process_time()\n",
    "            counts.append(c)\n",
    "            \n",
    "        else:\n",
    "            trace = False\n",
    "            x= y\n",
    "            counts.append(0)\n",
    "            \n",
    "        end = process_time()\n",
    "\n",
    "        time1.append(time1[t] + end - start)\n",
    "\n",
    "        if trace:\n",
    "            time2.append( time2[t] + end - start - (t2 -t1) -(t4 - t3) + t_pw)\n",
    "        else:\n",
    "            time2.append( time2[t] + end - start - (t2 -t1))\n",
    "            \n",
    "        f_improv = function_value[-1] - func(x)\n",
    "        function_value.append(func(x))\n",
    "        t+=1\n",
    "        \n",
    "    return x, function_value, time1,time2,t,primal_gap,counts,FW"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PGD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PDG  algorithm:\n",
    "def PGD(x, lmo,proj,L, epsilon,func,grad_f, f_tol, time_tol):\n",
    "    \n",
    "    #record primal gap, function value, and time every iteration\n",
    "    now=datetime.datetime.now()\n",
    "    primal_gap = []\n",
    "    function_value=[func(x)]\n",
    "    time = [0]\n",
    "    f_improv = np.inf\n",
    "    counts = []\n",
    "\n",
    "    #initialize starting point and active set\n",
    "    t = 0    \n",
    "\n",
    "    while f_improv > f_tol and time[-1] < time_tol:\n",
    "        #compute gradient\n",
    "        grad = grad_f(x)\n",
    "\n",
    "        #solve linear subproblem and compute FW direction\n",
    "        s = lmo(grad)\n",
    "        d_FW = s-x\n",
    "\n",
    "        #If primal gap is small enough - terminate\n",
    "        if np.dot(-grad,d_FW) <= epsilon:\n",
    "            break\n",
    "        else:\n",
    "            #update convergence data\n",
    "            primal_gap.append(np.dot(-grad,d_FW))\n",
    "\n",
    "        # take a 1/L step-size in negative gradient direction and then project back\n",
    "        x = proj(x,grad,L)\n",
    "        \n",
    "        later=datetime.datetime.now()\n",
    "        time.append((later - now).total_seconds())\n",
    "        f_improv = function_value[-1] - func(x)\n",
    "        function_value.append(func(x))\n",
    "        t+=1\n",
    "        \n",
    "    return x, function_value, time,t,primal_gap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "styles = {'AFW':'-', 'PFW':'--', 'DICG':'-', 'Shadow':'-',\n",
    "          'PGD':'--', 'Shadow CG':'-'}\n",
    "colors = {'AFW':'tab:blue', 'PFW':'tab:blue', 'DICG':'tab:red', 'Shadow':'tab:green',\n",
    "          'PGD':'tab:orange', 'Shadow CG':'tab:orange'}\n",
    "labels = {'AFW':'AFW', 'PFW':'PFW', 'DICG':'DICG', 'PGD':'PGD',\n",
    "          'Shadow':'Shadow Walk', 'Shadow CG':'Shadow CG'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotter(res, styles, colors, labels,oracle = False,fstar = None, log = True, sci = False, outfilename=None):\n",
    "    \n",
    "    '''\n",
    "    plot primal gap specifying whether we want to assume oracle access for shadow or not\n",
    "    '''\n",
    "    \n",
    "    fig = plt.figure()\n",
    "    fig.set_figheight(3)\n",
    "    fig.set_figwidth(10)\n",
    "\n",
    "    plt.subplot(1, 2, 1)\n",
    "    \n",
    "    if fstar:\n",
    "        for alg in res.keys():\n",
    "            plt.plot(np.arange(len(res[alg][1])), res[alg][1] - fstar, linestyle=styles[alg], color=colors[alg],\n",
    "                     label=labels[alg])\n",
    "    else:\n",
    "        for alg in res.keys():\n",
    "            plt.plot(np.arange(len(res[alg][1])), res[alg][1], linestyle=styles[alg], color=colors[alg],\n",
    "                     label=labels[alg])\n",
    "    plt.legend()\n",
    "    plt.ylabel(r'$\\log (f(\\mathbf{x}_t) - f(\\mathbf{x}^*))$')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.xlim(0,60)\n",
    "    plt.tick_params(which='minor', left=False)\n",
    "    if log:\n",
    "        plt.yscale('log')\n",
    "    if sci:\n",
    "        plt.ticklabel_format(axis = 'y',style='sci',scilimits = (0,0))\n",
    "    plt.grid(linestyle=':')\n",
    "\n",
    "    plt.subplot(1, 2, 2)\n",
    "    if fstar:\n",
    "        for alg in res.keys():\n",
    "            if oracle and alg in ['Shadow', 'Shadow CG']:\n",
    "                plt.plot(res[alg][3], res[alg][1] - fstar, linestyle=styles[alg], color=colors[alg])\n",
    "                plt.xlabel('Wall-clock time assuming shadow oracle (s)')\n",
    "            else:\n",
    "                plt.xlabel('Wall-clock time (s)')\n",
    "                plt.plot(res[alg][2], res[alg][1] - fstar, linestyle=styles[alg], color=colors[alg])\n",
    "    else:\n",
    "        for alg in res.keys():\n",
    "            \n",
    "            if oracle and alg in ['Shadow', 'Shadow CG']:\n",
    "                plt.plot(res[alg][3], res[alg][1], linestyle=styles[alg], color=colors[alg])\n",
    "                plt.xlabel('Wall-clock time assuming shadow oracle (s)')\n",
    "            else:\n",
    "                plt.plot(res[alg][2], res[alg][1], linestyle=styles[alg], color=colors[alg])\n",
    "                plt.xlabel('Wall-clock time (s)')\n",
    "                \n",
    "    plt.ylabel(r'$\\log (f(\\mathbf{x}_t) - f(\\mathbf{x}^*))$')\n",
    "    plt.tick_params(which='both', left=False, labelleft=False)\n",
    "    plt.grid(linestyle=':')\n",
    "    if log:\n",
    "        plt.yscale('log')\n",
    "    if sci:\n",
    "        plt.ticklabel_format(axis = 'y',style='sci',scilimits = (0,0))\n",
    "\n",
    "    if outfilename is not None:\n",
    "        plt.savefig(outfilename, dpi=200, bbox_inches='tight')\n",
    "        \n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def plotter_gaps(res, styles, colors, labels,oracle = False, log = True, sci = False, outfilename=None):\n",
    "    \n",
    "    \n",
    "    '''\n",
    "    plot duality gap specifying whether we want to assume oracle access for shadow or not\n",
    "    '''\n",
    "    \n",
    "    fig = plt.figure()\n",
    "    fig.set_figheight(3)\n",
    "    fig.set_figwidth(10)\n",
    "\n",
    "    plt.subplot(1, 2, 1)\n",
    "    for alg in res.keys():\n",
    "        if alg in ['Shadow', 'Shadow CG']:\n",
    "            plt.plot(np.arange(len(res[alg][5])), res[alg][5], linestyle=styles[alg], color=colors[alg],\n",
    "                     label=labels[alg])\n",
    "        else:\n",
    "            plt.plot(np.arange(len(res[alg][4])), res[alg][4], linestyle=styles[alg], color=colors[alg],\n",
    "                     label=labels[alg])\n",
    "    plt.legend()\n",
    "    plt.ylabel(r'$\\log(\\mathrm{Duality}$' + ' gap)')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.tick_params(which='minor', left=False)\n",
    "    if log:\n",
    "        plt.yscale('log')\n",
    "    if sci:\n",
    "        plt.ticklabel_format(axis = 'y',style='sci',scilimits = (0,0))\n",
    "    plt.grid(linestyle=':')\n",
    "\n",
    "    plt.subplot(1, 2, 2)\n",
    "    for alg in res.keys():\n",
    "        if oracle and alg in ['Shadow', 'Shadow CG']:\n",
    "            plt.plot(res[alg][3][1:], res[alg][5], linestyle=styles[alg], color=colors[alg])\n",
    "            plt.xlabel('Wall-clock time assuming shadow oracle (s)')\n",
    "        elif alg in ['Shadow', 'Shadow CG']:\n",
    "            plt.plot(res[alg][2][1:], res[alg][5], linestyle=styles[alg], color=colors[alg])\n",
    "            plt.xlabel('Wall-clock time (s)')\n",
    "        else:\n",
    "            plt.plot(res[alg][2][1:], res[alg][4], linestyle=styles[alg], color=colors[alg])\n",
    "            plt.xlabel('Wall-clock time (s)')\n",
    "    if log:\n",
    "        plt.yscale('log')\n",
    "    if sci:\n",
    "        plt.ticklabel_format(axis = 'y',style='sci',scilimits = (0,0))\n",
    "    plt.ylabel(r'$\\log(\\mathrm{Duality}$' + ' gap)')\n",
    "    plt.tick_params(which='both', left=False, labelleft=False)\n",
    "    plt.grid(linestyle=':')\n",
    "\n",
    "    if outfilename is not None:\n",
    "        plt.savefig(outfilename, dpi=200, bbox_inches='tight')\n",
    "        \n",
    "    plt.show()\n",
    "\n",
    "def plotter_break_point(res, styles, colors, labels, m = False,outfilename=None):\n",
    "    \n",
    "    '''\n",
    "    plot number of oracle calls made\n",
    "    '''\n",
    "    \n",
    "    fig = plt.figure()\n",
    "    fig.set_figheight(4)\n",
    "    fig.set_figwidth(6)\n",
    "    \n",
    "    li = ['Shadow', 'Shadow CG']\n",
    "    \n",
    "    for alg in li:\n",
    "        plt.plot(np.arange(len(res[alg][6])), np.array(res[alg][6]) + 1, linestyle=styles[alg], color=colors[alg],\n",
    "                 label=labels[alg])\n",
    "    if m:\n",
    "        plt.axhline(m,  color='black', linestyle='--', label = r'$m$')\n",
    "        \n",
    "\n",
    "    plt.legend()\n",
    "    #plt.ylabel('Number of Iterations spent in '+r'$\\mathrm{Trace}(\\mathbf{x}_t, \\nabla f(\\mathbf{x}_t))$')\n",
    "    plt.ylabel('Number of calls to shadow oracle')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.tick_params(which='minor', left=False)\n",
    "    plt.grid(linestyle=':')\n",
    "\n",
    "    if outfilename is not None:\n",
    "        plt.savefig(outfilename, dpi=200, bbox_inches='tight')\n",
    "        \n",
    "    plt.show()\n",
    "\n",
    "    \n",
    "def plotter_cum(res, styles, colors, labels, m = False,outfilename=None):\n",
    "    \n",
    "    '''\n",
    "    plot cumulative number of shadow steps taken\n",
    "    '''\n",
    "    \n",
    "    fig = plt.figure()\n",
    "    fig.set_figheight(4)\n",
    "    fig.set_figwidth(6)\n",
    "    \n",
    "    li = ['Shadow', 'Shadow CG']\n",
    "    \n",
    "    for alg in li:\n",
    "        data = []\n",
    "        for i in res[alg][6]:\n",
    "            if i >= 1:\n",
    "                data.append(i+1)\n",
    "            else:\n",
    "                data.append(i)\n",
    "        if alg == 'Shadow':       \n",
    "            plt.plot(np.arange(len(data)), np.cumsum(np.array(data)+1), linestyle=styles[alg], color=colors[alg],\n",
    "                 label=labels[alg])\n",
    "        else:      \n",
    "            plt.plot(np.arange(len(data)), np.cumsum(np.array(data)), linestyle=styles[alg], color=colors[alg],\n",
    "                 label=labels[alg])\n",
    "\n",
    "    plt.legend()\n",
    "    plt.ylabel('Cumulative number of shadow steps')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.tick_params(which='minor', left=False)\n",
    "    plt.grid(linestyle=':')\n",
    "\n",
    "    if outfilename is not None:\n",
    "        plt.savefig(outfilename, dpi=200, bbox_inches='tight')\n",
    "        \n",
    "    plt.show()\n",
    "\n",
    "    \n",
    "    \n",
    "def plotter_wout_PGD(res, styles, colors, labels,oracle = False,fstar = None, log = True, sci = False, outfilename=None):\n",
    "    \n",
    "    \n",
    "    '''\n",
    "    plot primal gap without PGD -  specifying whether we want to assume oracle access for shadow or not\n",
    "    '''\n",
    "    \n",
    "    \n",
    "    fig = plt.figure()\n",
    "    fig.set_figheight(3)\n",
    "    fig.set_figwidth(10)\n",
    "\n",
    "    plt.subplot(1, 2, 1)\n",
    "    \n",
    "    if fstar:\n",
    "        for alg in res.keys():\n",
    "            if alg != 'PGD':\n",
    "                plt.plot(np.arange(len(res[alg][1])), res[alg][1] - fstar, linestyle=styles[alg], color=colors[alg],\n",
    "                         label=labels[alg])\n",
    "    else:\n",
    "        for alg in res.keys():\n",
    "            if alg != 'PGD':\n",
    "                plt.plot(np.arange(len(res[alg][1])), res[alg][1], linestyle=styles[alg], color=colors[alg],\n",
    "                     label=labels[alg])\n",
    "    plt.legend()\n",
    "    plt.ylabel(r'$\\log (f(\\mathbf{x}_t) - f(\\mathbf{x}^*))$')\n",
    "    plt.xlabel('Iteration')\n",
    "    #plt.xlim(0,60)\n",
    "    plt.tick_params(which='minor', left=False)\n",
    "    if log:\n",
    "        plt.yscale('log')\n",
    "    if sci:\n",
    "        plt.ticklabel_format(axis = 'y',style='sci',scilimits = (0,0))\n",
    "    plt.grid(linestyle=':')\n",
    "\n",
    "    plt.subplot(1, 2, 2)\n",
    "    if fstar:\n",
    "        for alg in res.keys():\n",
    "            if alg != 'PGD':\n",
    "                if oracle and alg in ['Shadow', 'Shadow CG']:\n",
    "                    plt.plot(res[alg][3], res[alg][1] - fstar, linestyle=styles[alg], color=colors[alg])\n",
    "                    plt.xlabel('Wall-clock time assuming shadow oracle (s)')\n",
    "                else:\n",
    "                    plt.xlabel('Wall-clock time (s)')\n",
    "                    plt.plot(res[alg][2], res[alg][1] - fstar, linestyle=styles[alg], color=colors[alg])\n",
    "    else:\n",
    "        for alg in res.keys():\n",
    "            if alg != 'PGD':\n",
    "                if oracle and alg in ['Shadow', 'Shadow CG']:\n",
    "                    plt.plot(res[alg][3], res[alg][1], linestyle=styles[alg], color=colors[alg])\n",
    "                    plt.xlabel('Wall-clock time assuming shadow oracle (s)')\n",
    "                else:\n",
    "                    plt.xlabel('Wall-clock time (s)')\n",
    "                    plt.plot(res[alg][2], res[alg][1], linestyle=styles[alg], color=colors[alg])\n",
    "\n",
    "    plt.ylabel(r'$\\log (f(\\mathbf{x}_t) - f(\\mathbf{x}^*))$')\n",
    "    plt.tick_params(which='both', left=False, labelleft=False)\n",
    "    plt.grid(linestyle=':')\n",
    "    if log:\n",
    "        plt.yscale('log')\n",
    "    if sci:\n",
    "        plt.ticklabel_format(axis = 'y',style='sci',scilimits = (0,0))\n",
    "\n",
    "    if outfilename is not None:\n",
    "        plt.savefig(outfilename, dpi=200, bbox_inches='tight')\n",
    "        \n",
    "    plt.show()\n",
    "\n",
    "\n",
    "\n",
    "def plotter_gaps_wout_PGD(res, styles, colors, labels,oracle = False, log = True, sci = False, outfilename=None):\n",
    "    \n",
    "        \n",
    "    '''\n",
    "    plot duality gap without PGD -  specifying whether we want to assume oracle access for shadow or not\n",
    "    '''\n",
    "    \n",
    "    \n",
    "    fig = plt.figure()\n",
    "    fig.set_figheight(3)\n",
    "    fig.set_figwidth(10)\n",
    "\n",
    "    plt.subplot(1, 2, 1)\n",
    "    for alg in res.keys():\n",
    "        if alg != 'PGD':\n",
    "            if alg in ['Shadow', 'Shadow CG']:\n",
    "                plt.plot(np.arange(len(res[alg][5])), res[alg][5], linestyle=styles[alg], color=colors[alg],\n",
    "                         label=labels[alg])\n",
    "            else:\n",
    "                plt.plot(np.arange(len(res[alg][4])), res[alg][4], linestyle=styles[alg], color=colors[alg],\n",
    "                         label=labels[alg])\n",
    "    plt.legend()\n",
    "    plt.ylabel(r'$\\log(\\mathrm{Duality}$' + ' gap)')\n",
    "    plt.xlabel('Iteration')\n",
    "    plt.tick_params(which='minor', left=False)\n",
    "    if log:\n",
    "        plt.yscale('log')\n",
    "    if sci:\n",
    "        plt.ticklabel_format(axis = 'y',style='sci',scilimits = (0,0))\n",
    "    plt.grid(linestyle=':')\n",
    "\n",
    "    plt.subplot(1, 2, 2)\n",
    "    for alg in res.keys():\n",
    "        if alg != 'PGD':\n",
    "            if oracle and alg in ['Shadow', 'Shadow CG']:\n",
    "                plt.plot(res[alg][3][1:], res[alg][5], linestyle=styles[alg], color=colors[alg])\n",
    "                plt.xlabel('Wall-clock time assuming shadow oracle (s)')\n",
    "            elif alg in ['Shadow', 'Shadow CG']:\n",
    "                plt.xlabel('Wall-clock time (s)')\n",
    "                plt.plot(res[alg][2][1:], res[alg][5], linestyle=styles[alg], color=colors[alg])\n",
    "            else:\n",
    "                plt.xlabel('Wall-clock time (s)')\n",
    "                plt.plot(res[alg][2][1:], res[alg][4], linestyle=styles[alg], color=colors[alg])\n",
    "    if log:\n",
    "        plt.yscale('log')\n",
    "    if sci:\n",
    "        plt.ticklabel_format(axis = 'y',style='sci',scilimits = (0,0))\n",
    "    plt.ylabel(r'$\\log(\\mathrm{Duality}$' + ' gap)')\n",
    "    plt.tick_params(which='both', left=False, labelleft=False)\n",
    "    plt.grid(linestyle=':')\n",
    "\n",
    "    if outfilename is not None:\n",
    "        plt.savefig(outfilename, dpi=200, bbox_inches='tight')\n",
    "        \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Shadow and Projection oracles using FW"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#projection oracle using FW\n",
    "def proj_FW(x_0,y,lmo,f_tol,time_tol):\n",
    "    \n",
    "    '''\n",
    "    use the FW algorithm to project y onto polytope P using x_0 as a starting point\n",
    "    '''\n",
    "\n",
    "    def f_proj(z):\n",
    "        return np.dot(z - y,z - y)\n",
    "\n",
    "    def grad_f_proj(z):\n",
    "        return 2*(z - y)\n",
    "\n",
    "\n",
    "    return FW(x_0, lmo, 0,f_proj,grad_f_proj, f_tol, time_tol)[0]\n",
    "\n",
    "\n",
    "def shadow_oracle_FW_binary(x_0,grad_0,lmo,start,tol,f_tol,time_tol):\n",
    "    \n",
    "    \n",
    "    '''\n",
    "    use the FW algorithm to compute shadow of gradinet by doing binary search to find the breakpoint\n",
    "    '''\n",
    "\n",
    "    def f_proj(z):\n",
    "        return np.dot(z - x_0 + epsilon*grad_0,z - x_0 + epsilon*grad_0)\n",
    "\n",
    "    def grad_f_proj(z):\n",
    "        return 2*(z - x_0 + epsilon*grad_0)\n",
    "\n",
    "    epsilon = start\n",
    "\n",
    "    while True:\n",
    "\n",
    "        z1 = FW(x_0, lmo, 0,f_proj,grad_f_proj, f_tol, time_tol)[0]\n",
    "        d_pi = (z1 - x_0)/epsilon\n",
    "\n",
    "        epsilon = epsilon / 2\n",
    "        z2 = FW(x_0, lmo, 0,f_proj,grad_f_proj, f_tol, time_tol)[0]\n",
    "\n",
    "        if np.dot( np.round(x_0 + epsilon*d_pi -z2,5), np.round(x_0 + epsilon*d_pi - z2,5))**0.5  <= tol:\n",
    "            break\n",
    "        else:\n",
    "            epsilon = epsilon/2\n",
    "\n",
    "    return d_pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computational Experiments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figures 2, 5, 6 and 7:  Video Colocalization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load data\n",
    "data = loadmat('aeroplane_data_small.mat')\n",
    "A, b = data['A'], np.array([data['b'][i][0] for i in range(data['b'].shape[0])])\n",
    "\n",
    "\n",
    "#define objective function and gradient\n",
    "f = lambda x: np.dot(x, np.dot(A, x))/2+np.dot(b, x)\n",
    "grad_f = lambda x: np.dot(A, x)+b\n",
    "\n",
    "\n",
    "#create lmo oracle to find shortest path in DAG\n",
    "I = np.identity(20)\n",
    "x = np.concatenate([np.array([1]+[0]*19) for i in range(33)], axis=0)\n",
    "x = np.concatenate([np.array([0]*9+[1] + [0]*10) for i in range(33)], axis=0)\n",
    "\n",
    "lmo = lambda g: np.concatenate([I[np.argmin(g[20*i:20*i+20])] for i in range(33)], axis=0)\n",
    "\n",
    "#number of frames\n",
    "frames = 33\n",
    "\n",
    "#number of boxes:\n",
    "boxes = 20\n",
    "\n",
    "#number of edges\n",
    "n = frames*boxes\n",
    "\n",
    "#generate constraints for polytope\n",
    "B = -np.eye(n)\n",
    "for i in range(frames-1,-1,-1):\n",
    "    c = np.zeros(n)\n",
    "    for j in range(boxes*i, boxes*i + boxes):\n",
    "        c[j] = 1\n",
    "    B = np.vstack((c,B))\n",
    "\n",
    "#polytope is now given by Bx <= e\n",
    "B = np.vstack((-B[:frames],B))\n",
    "e = [-1]*frames + [1]*frames +  [0]*n\n",
    "\n",
    "#compute Lipschitz constant of objective\n",
    "L = max(np.linalg.eig(A)[0])\n",
    "\n",
    "#choose desired feasibility oracle \n",
    "def feasibility_oracle (x,d): \n",
    "    return max_step_size(B,e,x, d)\n",
    "\n",
    "# choose projection oracle\n",
    "def proj (x,grad,L):\n",
    "    return proj_FW(x,x-grad/L,lmo,1e-7,np.inf)\n",
    "\n",
    "#choose desired shadow oracle - for example we can oracle that uses FW\n",
    "def shadow(x,d):\n",
    "    return shadow_oracle_FW_binary(x,d,lmo,10,0.01,1e-7,np.inf)\n",
    "\n",
    "f_tol,time_tol,epsilon,tol = -1, np.inf, 0.0001, 0.0001"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "res = {}\n",
    "print('AFW')\n",
    "res['AFW'] = AFW(x, lmo, epsilon,f,grad_f, f_tol, time_tol)\n",
    "print('PFW')\n",
    "res['PFW'] = PFW(x, lmo, epsilon,f,grad_f, f_tol, time_tol=np.inf)\n",
    "print('DICG')\n",
    "res['DICG'] = DICG(x, lmo, feasibility_oracle, epsilon,f,grad_f, f_tol, time_tol)\n",
    "print('Shadow')\n",
    "res['Shadow'] = SD(x, lmo, shadow,feasibility_oracle, epsilon,f,grad_f, f_tol, time_tol,tol)\n",
    "print('Shadow CG')\n",
    "res['Shadow CG'] = SCG(x, lmo, shadow,feasibility_oracle, epsilon,f,grad_f, f_tol, time_tol,tol)\n",
    "print('PGD')\n",
    "res['PGD'] = PGD(x, lmo,proj,L, epsilon,f,grad_f, f_tol, time_tol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_star = f(res['Shadow'][0])\n",
    "\n",
    "#plot primal gaps and duality gaps without assuming oracle access to shadow\n",
    "plotter(res, styles, colors, labels, False, f_star, log = True)\n",
    "plotter_gaps(res, styles, colors, labels, False,log = True)\n",
    "\n",
    "#plot primal gaps and duality gaps assuming oracle access to shadow and leaving PGD out for a better comparison\n",
    "plotter_wout_PGD(res, styles, colors, labels, True, f_star, log = True)\n",
    "plotter_gaps_wout_PGD(res, styles, colors, labels, True,log = True)\n",
    "\n",
    "#plot number of shadow oracle calls made and access cumulative number of shadow steps taken by Shadow-CG and Shadow-Walk\n",
    "plotter_break_point(res, styles, colors, labels)\n",
    "plotter_cum(res, styles, colors, labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figures 11, 12 and 13: Lasso Regression (small instance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#specify problem dimension and noise parameter\n",
    "m, n, sigma = 40, 60, 0.1\n",
    "np.random.seed(4)\n",
    "\n",
    "#generate problem data\n",
    "x_star = np.zeros(n)\n",
    "rand = np.random.randint(1,n,int(0.25*n))\n",
    "for i in rand[:int(len(rand)*0.5)]:\n",
    "    x_star[i] = 1\n",
    "    \n",
    "for i in rand[:int(len(rand)*0.5)]:\n",
    "    x_star[i] = -1\n",
    "\n",
    "A = np.random.normal(0.5,0.25,(m, n))\n",
    "y = np.dot(A, x_star)+sigma*np.random.randn(m)\n",
    "\n",
    "\n",
    "#define objective function and gradient\n",
    "f = lambda z: np.linalg.norm(y-np.dot(A, z[:n]-z[n:]), 2)**2\n",
    "\n",
    "def grad_f(x):\n",
    "    return np.concatenate([2*np.matmul(A.T,np.dot(A, x[:n]-x[n:]) - y),-2*np.matmul(A.T,np.dot(A, x[:n]-x[n:]) - y)])\n",
    "\n",
    "\n",
    "#define linear optimization oracle\n",
    "tau = np.linalg.norm(x_star, 1)\n",
    "V = tau*np.identity(2*n)\n",
    "x = V[np.random.randint(len(V))].copy()\n",
    "lmo = lambda g: V[np.argmin(g)]\n",
    "    \n",
    "#generate constraint matrix:\n",
    "B = np.vstack((np.ones(2*n), -np.ones(2*n), -np.eye(2*n)))\n",
    "e = [tau] + [-tau] + [0]*(2*n)\n",
    "#polytope is now given by Bx <= e\n",
    "\n",
    "#choose desired feasibility oracle \n",
    "#since our polytope is a scaled simplex, the feasibility oracle could simply be the following\n",
    "def feasibility_oracle (x,d):\n",
    "    eta = np.min([-x[i]/d[i] if d[i] < 0 else np.inf for i in range(len(x))])\n",
    "    eta = min(eta, 1)\n",
    "    return eta\n",
    "\n",
    "#choose desired shadow oracle - here for example we can gurobi to compute shadow\n",
    "def shadow(x,d):\n",
    "    return shadow_oracle(B,e,x, -d)\n",
    "\n",
    "f_tol,time_tol,epsilon,tol = 0.00001, np.inf, 0.1, 0.001"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "res_lasso_sparse_small = {}\n",
    "print('AFW')\n",
    "res_lasso_sparse_small['AFW'] = AFW(x, lmo, epsilon,f,grad_f, f_tol, time_tol)\n",
    "print('PFW')\n",
    "res_lasso_sparse_small['PFW'] = PFW(x, lmo, epsilon,f,grad_f, f_tol, time_tol=np.inf)\n",
    "print('DICG')\n",
    "res_lasso_sparse_small['DICG'] = DICG(x, lmo, feasibility_oracle, epsilon,f,grad_f, f_tol, time_tol)\n",
    "print('Shadow')\n",
    "res_lasso_sparse_small['Shadow'] = SD(x, lmo, shadow,feasibility_oracle, epsilon,f,grad_f, f_tol, time_tol,tol)\n",
    "print('Shadow CG')\n",
    "res_lasso_sparse_small['Shadow CG'] = SCG(x, lmo, shadow,feasibility_oracle, epsilon,f,grad_f, f_tol, time_tol,tol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_star = f(res_lasso_sparse_small['AFW'][0])\n",
    "\n",
    "#plot primal gaps and duality gaps without assuming oracle access to shadow\n",
    "plotter_wout_PGD(res_lasso_sparse_small, styles, colors, labels, False, f_star, log = True)\n",
    "plotter_gaps_wout_PGD(res_lasso_sparse_small, styles, colors, labels, False,log = True)\n",
    "\n",
    "#plot primal gaps and duality gaps assuming oracle access to shadow and leaving PGD out for a better comparison\n",
    "plotter_wout_PGD(res_lasso_sparse_small, styles, colors, labels, True, f_star, log = True)\n",
    "plotter_gaps_wout_PGD(res_lasso_sparse_small, styles, colors, labels, True,log = True)\n",
    "\n",
    "#plot number of shadow oracle calls made and access cumulative number of shadow steps taken by Shadow-CG and Shadow-Walk\n",
    "plotter_break_point(res_lasso_sparse_small, styles, colors, labels)\n",
    "plotter_cum(res_lasso_sparse_small, styles, colors, labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figures 8, 9 and 10: Lasso Regression (large instance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#specify problem dimension and noise parameter\n",
    "m, n, sigma = 50, 100, 0.1\n",
    "np.random.seed(11)\n",
    "\n",
    "#generate problem data\n",
    "x_star = np.zeros(n)\n",
    "rand = np.random.randint(1,n,int(0.25*n))\n",
    "for i in rand[:int(len(rand)*0.5)]:\n",
    "    x_star[i] = 1\n",
    "    \n",
    "for i in rand[:int(len(rand)*0.5)]:\n",
    "    x_star[i] = -1\n",
    "\n",
    "A = np.random.normal(0.5,0.25,(m, n))\n",
    "y = np.dot(A, x_star)+sigma*np.random.randn(m)\n",
    "\n",
    "\n",
    "#define objective function and gradient\n",
    "f = lambda z: np.linalg.norm(y-np.dot(A, z[:n]-z[n:]), 2)**2\n",
    "\n",
    "def grad_f(x):\n",
    "    return np.concatenate([2*np.matmul(A.T,np.dot(A, x[:n]-x[n:]) - y),-2*np.matmul(A.T,np.dot(A, x[:n]-x[n:]) - y)])\n",
    "\n",
    "\n",
    "#define linear optimization oracle\n",
    "tau = np.linalg.norm(x_star, 1)\n",
    "V = tau*np.identity(2*n)\n",
    "x = V[np.random.randint(len(V))].copy()\n",
    "lmo = lambda g: V[np.argmin(g)]\n",
    "    \n",
    "#generate constraint matrix:\n",
    "B = np.vstack((np.ones(2*n), -np.ones(2*n), -np.eye(2*n)))\n",
    "e = [tau] + [-tau] + [0]*(2*n)\n",
    "#polytope is now given by Bx <= e\n",
    "\n",
    "#choose desired feasibility oracle \n",
    "#since our polytope is a scaled simplex, the feasibility oracle could simply be the following\n",
    "def feasibility_oracle (x,d):\n",
    "    eta = np.min([-x[i]/d[i] if d[i] < 0 else np.inf for i in range(len(x))])\n",
    "    eta = min(eta, 1)\n",
    "    return eta\n",
    "\n",
    "#choose desired shadow oracle - here for example we can gurobi to compute shadow\n",
    "def shadow(x,d):\n",
    "    return shadow_oracle(B,e,x, -d)\n",
    "\n",
    "f_tol,time_tol,epsilon,tol = 0.00001, np.inf, 0.1, 0.001"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "res_lasso_sparse_large = {}\n",
    "print('AFW')\n",
    "res_lasso_sparse_large['AFW'] = AFW(x, lmo, epsilon,f,grad_f, f_tol, time_tol)\n",
    "print('PFW')\n",
    "res_lasso_sparse_large['PFW'] = PFW(x, lmo, epsilon,f,grad_f, f_tol, time_tol=np.inf)\n",
    "print('DICG')\n",
    "res_lasso_sparse_large['DICG'] = DICG(x, lmo, feasibility_oracle, epsilon,f,grad_f, f_tol, time_tol)\n",
    "print('Shadow')\n",
    "res_lasso_sparse_large['Shadow'] = SD(x, lmo, shadow,feasibility_oracle, epsilon,f,grad_f, f_tol, time_tol, tol)\n",
    "print('Shadow CG')\n",
    "res_lasso_sparse_large['Shadow CG'] = SCG(x, lmo, shadow,feasibility_oracle, epsilon,f,grad_f, f_tol, time_tol, tol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_star = f(res_lasso_sparse_large['Shadow'][0])\n",
    "\n",
    "#plot primal gaps and duality gaps without assuming oracle access to shadow\n",
    "plotter_wout_PGD(res_lasso_sparse_large, styles, colors, labels, False, f_star, log = True)\n",
    "plotter_gaps_wout_PGD(res_lasso_sparse_large, styles, colors, labels, False,log = True)\n",
    "\n",
    "#plot primal gaps and duality gaps assuming oracle access to shadow and leaving PGD out for a better comparison\n",
    "plotter_wout_PGD(res_lasso_sparse_large, styles, colors, labels, True, f_star, log = True)\n",
    "plotter_gaps_wout_PGD(res_lasso_sparse_large, styles, colors, labels, True,log = True)\n",
    "\n",
    "#plot number of shadow oracle calls made and access cumulative number of shadow steps taken by Shadow-CG and Shadow-Walk\n",
    "plotter_break_point(res_lasso_sparse_large, styles, colors, labels)\n",
    "plotter_cum(res_lasso_sparse_large, styles, colors, labels)"
   ]
  }
 ],
 "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}