2021_03_06_dual_approach.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A Dual Approach to Holistic Regression\n",
    "## March 6th, 2021"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Random, Distributions\n",
    "using LinearAlgebra\n",
    "using Gurobi, JuMP\n",
    "using DataFrames\n",
    "using CSV\n",
    "using StatsBase\n",
    "using Plots\n",
    "using ProgressBars\n",
    "using Optim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "--------------------------------------------\n",
      "Warning: your license will expire in 3 days\n",
      "--------------------------------------------\n",
      "\n",
      "Academic license - for non-commercial use only\n"
     ]
    }
   ],
   "source": [
    "gurobi_env = Gurobi.Env()\n",
    "\n",
    "function create_gurobi_model(; TimeLimit=-1, LogFile=nothing)\n",
    "    model = Model(optimizer_with_attributes(() -> Gurobi.Optimizer(gurobi_env)));\n",
    "    if TimeLimit >= 0\n",
    "        println(\"Set Gurobi TimeLimit.\")\n",
    "        set_optimizer_attribute(model, \"TimeLimit\", TimeLimit)\n",
    "    end\n",
    "    if LogFile != nothing\n",
    "        println(\"LogFile: $(LogFile).\")\n",
    "        set_optimizer_attribute(model, \"LogFile\", LogFile)\n",
    "    else\n",
    "        set_optimizer_attribute(model, \"OutputFlag\", 0)\n",
    "    end\n",
    "    set_optimizer_attribute(model, \"NumericFocus\", 3)\n",
    "    return model\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "____\n",
    "## 0. Utils"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "write_to_file (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function write_list(file_path, l)\n",
    "\n",
    "    if length(l) == 0\n",
    "        return\n",
    "    end\n",
    "    \n",
    "    str = \"\"\n",
    "    for e in l\n",
    "        try e = round(e, digits=3) catch _ end\n",
    "        str *= \"$(e),\"\n",
    "    end\n",
    "    open(file_path, \"a+\") do io\n",
    "        write(io, chop(str)*\"\\n\")\n",
    "    end\n",
    "    \n",
    "    return\n",
    "end\n",
    "\n",
    "function write_to_file(file_path, str)\n",
    "    open(file_path, \"a+\") do io\n",
    "        write(io, str)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([3, 4], [1, 2, 5])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function get_support(s)\n",
    "    supp = similar(s, Int)\n",
    "    count_supp = 1\n",
    "    \n",
    "    supp_c = similar(s, Int)\n",
    "    count_supp_c = 1\n",
    "    \n",
    "    @inbounds for i in eachindex(s)\n",
    "        supp[count_supp] = i\n",
    "        supp_c[count_supp_c] = i\n",
    "        is_zero = s[i] < 0.5\n",
    "        count_supp += !is_zero\n",
    "        count_supp_c += is_zero\n",
    "    end\n",
    "    return resize!(supp, count_supp-1), resize!(supp_c, count_supp_c-1)\n",
    "end\n",
    "\n",
    "get_support([0, 0, 1, 1, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "____\n",
    "## 1. Generate Synthetic Data  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "function generate_synthetic_data(n, p, k, NR; seed=-1)\n",
    "    \"\"\"\n",
    "        n = num. of samples\n",
    "        p = num. of features\n",
    "        k = num. of non zero coefficients\n",
    "        NR = noise ratio ~ σ_noise = NR * σ_y_true\n",
    "    \"\"\"\n",
    "    if seed >= 0\n",
    "        Random.seed!(seed)\n",
    "    end\n",
    "    \n",
    "    # Generate PD matrix\n",
    "    A = randn(p, p)\n",
    "    A = A'*A\n",
    "    Σ = (A' + A)/2\n",
    "    \n",
    "    # Generate data X\n",
    "    d = MvNormal(Σ)\n",
    "    X = rand(d, n)'I\n",
    "    \n",
    "    # Split data\n",
    "    index_train = 1:floor(Int, 0.5*n)\n",
    "    index_val = floor(Int, 0.5*n)+1:floor(Int, 0.75*n)\n",
    "    index_test = floor(Int, 0.75*n)+1:n\n",
    "    \n",
    "    X_train = X[index_train,:]\n",
    "    X_val = X[index_val,:]\n",
    "    X_test = X[index_test,:]\n",
    "    \n",
    "    # Center\n",
    "    μ_train = [mean(X_train[:, j]) for j=1:p]\n",
    "    for j=1:p\n",
    "         X_train[:,j] = X_train[:,j] .- μ_train[j]\n",
    "         X_val[:,j] = X_val[:,j] .- μ_train[j]\n",
    "         X_test[:,j] = X_test[:,j] .- μ_train[j]\n",
    "    end\n",
    "    \n",
    "    # Scale\n",
    "    σ_train = [norm(X_train[:, j]) for j=1:p]\n",
    "    for j=1:p\n",
    "         X_train[:,j] = X_train[:,j]/σ_train[j]\n",
    "         X_val[:,j] = X_val[:,j] ./ σ_train[j]\n",
    "         X_test[:,j] = X_test[:,j] ./ σ_train[j]\n",
    "    end\n",
    "    \n",
    "    # Generate β\n",
    "    β = zeros(p)\n",
    "    for j=1:k\n",
    "        β[floor(Int, j*p/k)] = 1.0*rand([-1, 1])\n",
    "    end\n",
    "    \n",
    "    # Noise\n",
    "    ϵ = rand(Normal(0, std(X*β)*NR), n)\n",
    "    \n",
    "    # Target\n",
    "    y_train = X_train*β + ϵ[index_train]\n",
    "    y_val = X_val*β + ϵ[index_val]\n",
    "    y_test = X_test*β + ϵ[index_test]\n",
    "            \n",
    "    return  (X_train, y_train), (X_val, y_val), (X_test, y_test), β\n",
    "end\n",
    "\n",
    "function get_t_α(n, p, α)\n",
    "    return quantile(TDist(n-p), 1 - α/2)\n",
    "end\n",
    "\n",
    "function get_σ_X(X, y, γ)\n",
    "    n, p = size(X)\n",
    "    \n",
    "    # Estimator σ\n",
    "    M_inv = inv(I/γ + X'X)\n",
    "    σ_tilde = sqrt((y'*(I - X*M_inv*X')*y)/(n-p))\n",
    "    σ_X = σ_tilde * sqrt.(diag(M_inv))\n",
    "    \n",
    "    return σ_X\n",
    "end\n",
    "\n",
    "function get_R2(y_pred, y_true, y_train)\n",
    "    SS_res = norm(y_true .- y_pred)\n",
    "    SS_tot = norm(y_true .- mean(y_train))\n",
    "    return 1 - (SS_res/SS_tot)^2\n",
    "end\n",
    "\n",
    ";"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Compute inner problems and gradients"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### a. Compute g_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "g_s (generic function with 1 method)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function g_s(D_s, b_s, σ_X_s; GD=true)\n",
    "    \n",
    "    # Get length of support of s\n",
    "    l = length(b_s)\n",
    "    \n",
    "    # Case s == 0\n",
    "    if l==0\n",
    "        return zeros(0), 0.0\n",
    "    end\n",
    "    \n",
    "    # Initial solution\n",
    "    λ_s0 = zeros(l) .+ 1.0\n",
    "\n",
    "    # Compute objective and gradient at the same time\n",
    "    function fg!(F, G, λ_s)\n",
    "        \n",
    "        μ_s = λ_s .+ b_s\n",
    "        β_s = D_s*μ_s\n",
    "        \n",
    "        if G != nothing\n",
    "            G .= β_s .- σ_X_s\n",
    "        end\n",
    "        \n",
    "        if F != nothing\n",
    "            return -λ_s'σ_X_s + 0.5*μ_s'β_s\n",
    "        end\n",
    "    end\n",
    "    \n",
    "    # Lagrangian multipliers constraint\n",
    "    lower = zeros(l)\n",
    "    upper = [Inf for _ in 1:l]\n",
    "\n",
    "    res = Optim.optimize(Optim.only_fg!(fg!), lower, upper, λ_s0, \n",
    "        Fminbox(GD ? GradientDescent() : LBFGS()))\n",
    "\n",
    "    return Optim.minimizer(res), - Optim.minimum(res)\n",
    "    \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "function g_s_gurobi(D_s, b_s, σ_X_s, model)\n",
    "    \n",
    "    # Get length of support of s\n",
    "    l = length(b_s)\n",
    "    if l==0\n",
    "        return zeros(0), 0.0\n",
    "    end\n",
    "\n",
    "    λ_s = model[:λ][1:l]\n",
    "    μ_s = λ_s .+ b_s\n",
    "    β_s = D_s*μ_s\n",
    "    \n",
    "    @objective(model, Max, λ_s'σ_X_s - 0.5*μ_s'β_s)\n",
    "    \n",
    "    optimize!(model)\n",
    "    \n",
    "    value.(λ_s), objective_value(model)\n",
    "end;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "function g_gurobi(supp, Z, D, b, σ_X, model)\n",
    "\n",
    "    # Create DZ once\n",
    "    DZ = D*Z\n",
    "    \n",
    "    λ = model[:λ]\n",
    "    μ = b + λ\n",
    "\n",
    "    @objective(model, Max, λ'*Z*σ_X - 0.5μ'*DZ*μ)\n",
    "    \n",
    "    optimize!(model)\n",
    "    \n",
    "    value.(λ)[supp], objective_value(model)\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### b. Compute ∇g_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "∇g_s (generic function with 1 method)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function ∇g_s(supp, supp_c, b, M, λ_s, D_s, σ_X_s, γ)\n",
    "    \n",
    "    β_s = D_s*(b[supp] .+ λ_s)\n",
    "  \n",
    "    grad = zeros(length(b))\n",
    "    grad[supp] = λ_s .* σ_X_s - (β_s .^ 2)/(2γ)\n",
    "    grad[supp_c] = - 0.5*γ*(b[supp_c] - M[supp_c, supp]*β_s).^2\n",
    "    \n",
    "    return grad\n",
    "    \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "∇g (generic function with 1 method)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function ∇g(supp, D, b, λ_s, σ_X, γ)\n",
    "    \n",
    "    λ = zeros(length(b))\n",
    "    λ[supp] = λ_s\n",
    "    \n",
    "    grad = λ .* σ_X - ((D'*(b + λ)).^ 2)/(2γ)\n",
    "\n",
    "    return grad\n",
    "    \n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_____\n",
    "## 3. Compare speed "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### /!\\ t_α is already in σ_X /!\\"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "n_train = 10000\n",
    "n = 2*n_train\n",
    "p = 100\n",
    "k = 10\n",
    "NR = 0.001\n",
    "α = 0.05\n",
    "t_α = get_t_α(n_train, p, α)\n",
    "γ = 1/sqrt(n_train)\n",
    "\n",
    "# Generate data\n",
    "(X_p, y), _, _, β_true = generate_synthetic_data(n, p, k, NR, seed=42);\n",
    "σ_X_p = t_α * get_σ_X(X_p, y, γ); #t_α is already in σ_X\n",
    "\n",
    "# Compute data in p dimensions\n",
    "M_p = X_p'X_p\n",
    "b_p = X_p'y\n",
    "\n",
    "# Compute data in 2p dimensions\n",
    "M = [M_p -M_p; -M_p  M_p]\n",
    "b = [b_p; -b_p];\n",
    "σ_X = [σ_X_p; σ_X_p] ;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create s\n",
    "s_true = vcat(β_true .> 0, β_true .< 0) .* 1 \n",
    "supp, supp_c = get_support(s_true)\n",
    "\n",
    "# Get projected variables\n",
    "b_s = b[supp];\n",
    "σ_X_s = σ_X[supp];\n",
    "D_s = inv(I/γ + M[supp, supp]);\n",
    "\n",
    "# Create model for g_s_gurobi\n",
    "model_inner_g_s = create_gurobi_model();\n",
    "@variable(model_inner_g_s, λ[1:k] >= 0)\n",
    "\n",
    "# Create model for g_gurobi\n",
    "model_inner_g = create_gurobi_model();\n",
    "@variable(model_inner_g, λ[1:2p] >= 0);\n",
    "\n",
    "Z = Diagonal(s_true);\n",
    "D = inv(I/γ + Z*M);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare g_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GD: -0.04579702702244497\n",
      "LBFGS: -0.04579702702244497\n",
      "Gurobi (g_s): -0.04579702842816385\n",
      "Gurobi (g): -0.04579702842816385\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "10×4 Array{Float64,2}:\n",
       " 1.67444e-18  6.68598e-18  3.15464e-8  3.15464e-8\n",
       " 0.112048     0.112048     0.112048    0.112048  \n",
       " 1.02156e-20  6.86706e-17  4.4832e-8   4.4832e-8 \n",
       " 0.0750941    0.0750941    0.0750945   0.0750945 \n",
       " 0.0427485    0.0427485    0.0427494   0.0427494 \n",
       " 0.042549     0.042549     0.04255     0.04255   \n",
       " 0.0776394    0.0776394    0.0776398   0.0776398 \n",
       " 1.6996e-17   3.3594e-17   4.57974e-6  4.57974e-6\n",
       " 4.70998e-20  1.13528e-17  4.41367e-7  4.41367e-7\n",
       " 1.63789e-18  3.9433e-17   1.26127e-8  1.26127e-8"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "λ_s_GD, g_s_GD = g_s(D_s, b_s, σ_X_s; GD=true)\n",
    "λ_s_LBFGS, g_s_LBFGS = g_s(D_s, b_s, σ_X_s; GD=false)\n",
    "λ_s_guro_s, g_s_guro_s = g_s_gurobi(D_s, b_s, σ_X_s, model_inner_g_s)\n",
    "λ_s_guro, g_s_guro = g_gurobi(supp, Z, D, b, σ_X, model_inner_g)\n",
    "\n",
    "# Compare objective values\n",
    "println(\"GD: \",g_s_GD)\n",
    "println(\"LBFGS: \",g_s_LBFGS)\n",
    "println(\"Gurobi (g_s): \", g_s_guro_s)\n",
    "println(\"Gurobi (g): \", g_s_guro)\n",
    "\n",
    "# Compare λ\n",
    "hcat(λ_s_GD, λ_s_LBFGS, λ_s_guro_s, λ_s_guro)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compute D_s: 0.023847533000000008\n",
      "Compute D: 3.6927635900000046\n"
     ]
    }
   ],
   "source": [
    "# Computational time for D_s compared to D (given s)\n",
    "\n",
    "total_time_D_s = 0\n",
    "total_time_D = 0\n",
    "for _ in 1:1000\n",
    "    \n",
    "    total_time_D_s += @elapsed begin \n",
    "        supp, supp_c = get_support(s_true)\n",
    "        b_s = b[supp];\n",
    "        σ_X_s = σ_X[supp];\n",
    "        D_s = inv(I/γ + M[supp, supp]);\n",
    "    end\n",
    "    \n",
    "    total_time_D += @elapsed begin\n",
    "        Z = Diagonal(s_true);\n",
    "        D = inv(I/γ + Z*M);\n",
    "    end\n",
    "end\n",
    "\n",
    "println(\"Compute D_s: \", total_time_D_s)\n",
    "println(\"Compute D: \", total_time_D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optim.jl + GD: 4.810573553999998\n",
      "Optim.jl + LBFGS: 1.779046216000001\n",
      "Gurobi g_s (model created outside): 1.3252974939999995\n",
      "Gurobi g (model created outside): 41.289308443000074\n"
     ]
    }
   ],
   "source": [
    "# Computational time of g_s and g for given D_s or D\n",
    "\n",
    "total_time_GD = 0\n",
    "total_time_LBFGS = 0\n",
    "total_time_Gurobi_g_s = 0\n",
    "total_time_Gurobi_g = 0\n",
    "\n",
    "for _ in 1:1000\n",
    "    total_time_GD += @elapsed g_s(D_s, b_s, σ_X_s; GD=true)\n",
    "    total_time_LBFGS += @elapsed g_s(D_s, b_s, σ_X_s; GD=false)\n",
    "    total_time_Gurobi_g_s += @elapsed g_s_gurobi(D_s, b_s, σ_X_s, model_inner_g_s)\n",
    "    total_time_Gurobi_g += @elapsed g_gurobi(supp, Z, D, b, σ_X, model_inner_g)\n",
    "end\n",
    "\n",
    "println(\"Optim.jl + GD: \",total_time_GD)\n",
    "println(\"Optim.jl + LBFGS: \", total_time_LBFGS)\n",
    "println(\"Gurobi g_s (model created outside): \", total_time_Gurobi_g_s)\n",
    "println(\"Gurobi g (model created outside): \", total_time_Gurobi_g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare ∇g_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "|| ∇g_s_guro - ∇g_s_GD || =  7.28425921268618e-10\n",
      "|| ∇g_s_guro - ∇g_guro || =  6.151124447574778e-18\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "10×3 Array{Float64,2}:\n",
       " -0.000204746  -0.000204746  -0.000204746\n",
       " -0.000297864  -0.000297864  -0.000297864\n",
       " -5.99079e-6   -5.99079e-6   -5.99079e-6 \n",
       " -0.000144345  -0.000144345  -0.000144345\n",
       " -2.23807e-5   -2.23807e-5   -2.23807e-5 \n",
       " -0.000190836  -0.000190836  -0.000190836\n",
       " -2.47693e-6   -2.47693e-6   -2.47693e-6 \n",
       " -3.31235e-6   -3.31235e-6   -3.31235e-6 \n",
       " -4.72655e-5   -4.72655e-5   -4.72655e-5 \n",
       " -0.00671411   -0.00671411   -0.00671411 "
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Compare the gradient for different optimal lambdas\n",
    "\n",
    "∇g_s_guro = ∇g_s(supp, supp_c, b, M, λ_s_guro, D_s, σ_X_s, γ)\n",
    "∇g_s_GD = ∇g_s(supp, supp_c, b, M, λ_s_GD, D_s, σ_X_s, γ)\n",
    "∇g_guro = ∇g(supp, D, b, λ_s_guro, σ_X, γ)\n",
    "\n",
    "println(\"|| ∇g_s_guro - ∇g_s_GD || =  \", norm(∇g_s_guro - ∇g_s_GD))\n",
    "println(\"|| ∇g_s_guro - ∇g_guro || =  \", norm(∇g_s_guro - ∇g_guro))\n",
    "hcat(∇g_s_GD, ∇g_s_guro, ∇g_guro)[1:10, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compute ∇g_s: 0.023095553000000008\n",
      "Compute ∇g: 0.02642997599999998\n"
     ]
    }
   ],
   "source": [
    "# Computational time for D_s compared to D (given s)\n",
    "\n",
    "total_time_∇g_s = 0\n",
    "total_time_∇g = 0\n",
    "\n",
    "for _ in 1:1000\n",
    "    \n",
    "    total_time_∇g_s += @elapsed begin \n",
    "        ∇g_s(supp, supp_c, b, M, λ_s_GD, D_s, σ_X_s, γ)\n",
    "    end\n",
    "    \n",
    "    total_time_∇g += @elapsed begin\n",
    "        ∇g(supp, D, b, λ_s_GD, σ_X, γ)\n",
    "    end\n",
    "end\n",
    "\n",
    "println(\"Compute ∇g_s: \", total_time_∇g_s)\n",
    "println(\"Compute ∇g: \", total_time_∇g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_____\n",
    "## 4. Compute Cutting plane algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "compute_primal (generic function with 1 method)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function compute_primal(X, y, k, γ, σ_X; TimeLimit=-1, LogFile=nothing)\n",
    "    \n",
    "    n, p = size(X)\n",
    "    \n",
    "    model = create_gurobi_model(;LogFile=LogFile, TimeLimit=TimeLimit)\n",
    "\n",
    "    # TODO: change big-M values\n",
    "    M1 = 1000\n",
    "    M2 = 1000\n",
    "\n",
    "    @variable(model, β[i=1:p])\n",
    "    @variable(model, s[i=1:p], Bin)\n",
    "    @variable(model, b[i=1:p], Bin)\n",
    "\n",
    "    @constraint(model, sum(s) <= k)\n",
    "    \n",
    "    @constraint(model, [i=1:p], β[i] <= M1*s[i])\n",
    "    @constraint(model, [i=1:p], β[i] >= -M1*s[i])\n",
    "\n",
    "    @constraint(model, [i=1:p], β[i]/σ_X[i] + M2*b[i] >= s[i])\n",
    "    @constraint(model, [i=1:p], -β[i]/σ_X[i] + M2*(1-b[i]) >= s[i])\n",
    "\n",
    "    yty = y'y\n",
    "    XtX = X'X\n",
    "    ytX = y'X\n",
    "    \n",
    "    @objective(model, Min, 0.5*(yty - 2*ytX*β + β'*XtX*β + (1/γ)*sum(β[j]^2 for j=1:p)))\n",
    "        \n",
    "    JuMP.optimize!(model)\n",
    "    \n",
    "    return objective_value(model), value.(β)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "compute_warm_start_primal (generic function with 1 method)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function compute_warm_start_primal(M_p, b_p, k, γ, σ_X_p, time_limit; LogFile=nothing)\n",
    "    \n",
    "    p = length(b_p)\n",
    "    \n",
    "    model = create_gurobi_model(;TimeLimit=time_limit, LogFile=LogFile)\n",
    "\n",
    "    # TODO: change big-M values\n",
    "    M1 = 1000\n",
    "    M2 = 1000\n",
    "\n",
    "    @variable(model, β[i=1:p])\n",
    "    @variable(model, s[i=1:p], Bin)\n",
    "    @variable(model, b[i=1:p], Bin)\n",
    "\n",
    "    @constraint(model, sum(s) <= k)\n",
    "    \n",
    "    @constraint(model, [i=1:p], β[i] <= M1*s[i])\n",
    "    @constraint(model, [i=1:p], β[i] >= -M1*s[i])\n",
    "\n",
    "    @constraint(model, [i=1:p], β[i]/σ_X_p[i] + M2*b[i] >= s[i])\n",
    "    @constraint(model, [i=1:p], -β[i]/σ_X_p[i] + M2*(1-b[i]) >= s[i])\n",
    "    \n",
    "    @objective(model, Min, 0.5*(- 2*b_p'β + β'*M_p*β + (1/γ)*sum(β[j]^2 for j=1:p)))\n",
    "    \n",
    "    optimize!(model)\n",
    "    \n",
    "    s_val = Int.(value.(s))\n",
    "    b_val = Int.(value.(b))\n",
    "    \n",
    "    return vcat(s_val .* (b_val .== 0), s_val .* (b_val .== 1))\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "compute_dual"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "WarmStart ∈ { :None, :RidgeStart, :PrimalStart }\n",
    "\"\"\"\n",
    "function compute_dual(X_p, y, k, γ, σ_X_p; LogFile=nothing, WarmStart=:None, TimeLimit=-1)\n",
    "    \n",
    "    # Get dimensions\n",
    "    n, p = size(X_p)\n",
    "    \n",
    "    # Constant\n",
    "    C = 0.5*y'y\n",
    " \n",
    "    # Compute data in p dimensions\n",
    "    M_p = X_p'X_p\n",
    "    b_p = X_p'y\n",
    "\n",
    "    # Compute data in 2p dimensions\n",
    "    M = [M_p -M_p; -M_p  M_p]\n",
    "    b = [b_p; -b_p];\n",
    "    σ_X = [σ_X_p; σ_X_p]\n",
    "    \n",
    "    # Outer problem\n",
    "    miop = create_gurobi_model(;LogFile=LogFile, TimeLimit=TimeLimit)\n",
    "    @variable(miop, s[1:2p], Bin)\n",
    "    @variable(miop, t >= -C)\n",
    "    @constraint(miop, sum(s) <= k)\n",
    "    @constraint(miop, [i=1:p], s[i]+s[p+i]<=1)\n",
    "    \n",
    "    # Initial solution\n",
    "    if (WarmStart == :RidgeStart)\n",
    "        s0 = zeros(2p)\n",
    "        β_ridge = inv(I/γ + M_p)*b_p\n",
    "        s0[findall(x -> x>0, β_ridge)] .= 1.0\n",
    "        s0[findall(x -> x<0, β_ridge) .+ p] .= 1.0\n",
    "    \n",
    "    elseif (WarmStart == :PrimalStart)\n",
    "        s0 = compute_warm_start_primal(M_p, b_p, k, γ, σ_X_p, 20; LogFile=LogFile) # TODO: change time limit\n",
    "    else \n",
    "        s0 = zeros(2p)\n",
    "        s0[1:k] .= 1\n",
    "    end\n",
    "\n",
    "    # Initial cut\n",
    "    supp, supp_c = get_support(s0)\n",
    "    D_s, b_s, σ_X_s = inv(I/γ + M[supp, supp]), b[supp], σ_X[supp]\n",
    "    λ_s0, g_s0 = g_s(D_s, b_s, σ_X_s; GD=true)\n",
    "    ∇g_s0 = ∇g_s(supp, supp_c, b, M, λ_s0, D_s, σ_X_s, γ)\n",
    "    offset = g_s0 - sum(∇g_s0[supp])\n",
    "    \n",
    "    @constraint(miop, t >= ∇g_s0's + offset)\n",
    "    @objective(miop, Min, t + C)\n",
    "    \n",
    "    # Cutting planes    \n",
    "    function outer_approximation(cb_data)\n",
    "        \n",
    "        # Get feasible solution\n",
    "        s_val = [callback_value(cb_data, s[i]) for i=1:2p]\n",
    "        \n",
    "        # Generate cut\n",
    "        supp, supp_c = get_support(s_val)\n",
    "        D_s, b_s, σ_X_s = inv(I/γ + M[supp, supp]), b[supp], σ_X[supp]\n",
    "        λ_s_val, g_s_val = g_s(D_s, b_s, σ_X_s; GD=true)\n",
    "        ∇g_s_val = ∇g_s(supp, supp_c, b, M, λ_s_val, D_s, σ_X_s, γ)\n",
    "        offset = g_s_val - sum(∇g_s_val[supp])\n",
    "        \n",
    "        con = @build_constraint(t >= ∇g_s_val's + offset)\n",
    "        MOI.submit(miop, MOI.LazyConstraint(cb_data), con)\n",
    "        \n",
    "    end\n",
    "    \n",
    "    MOI.set(miop, MOI.LazyConstraintCallback(), outer_approximation)\n",
    "    optimize!(miop)\n",
    "    \n",
    "    s_opt = value.(miop[:s])\n",
    "    \n",
    "    supp, supp_c = get_support(s_opt)\n",
    "    D_s, b_s, σ_X_s = inv(I/γ + M[supp, supp]), b[supp], σ_X[supp]\n",
    "    λ_s_opt, g_s_opt = g_s(D_s, b_s, σ_X_s; GD=true)\n",
    "    β = zeros(2p)\n",
    "    β[supp] = D_s*(λ_s_opt + b_s)\n",
    "    \n",
    "    return objective_value(miop), β[1:p] - β[p+1:end]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Set Gurobi TimeLimit.\n",
      "LogFile: debug.txt.\n",
      "Gurobi Optimizer version 9.0.3 build v9.0.3rc0 (mac64)\n",
      "Optimize a model with 401 rows, 300 columns and 1100 nonzeros\n",
      "Model fingerprint: 0x526306df\n",
      "Model has 5050 quadratic objective terms\n",
      "Variable types: 100 continuous, 200 integer (200 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 1e+03]\n",
      "  Objective range  [1e-02, 1e+00]\n",
      "  QObjective range [5e-06, 1e+02]\n",
      "  Bounds range     [0e+00, 0e+00]\n",
      "  RHS range        [1e+01, 1e+03]\n",
      "Found heuristic solution: objective 9.3237205\n",
      "Presolve time: 0.00s\n",
      "Presolved: 401 rows, 300 columns, 1100 nonzeros\n",
      "Presolved model has 5050 quadratic objective terms\n",
      "Variable types: 100 continuous, 200 integer (200 binary)\n",
      "\n",
      "Root relaxation: objective 9.233698e+00, 486 iterations, 0.03 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0    9.23370    0  159    9.32372    9.23370  0.97%     -    0s\n",
      "H    0     0                       9.2891843    9.23370  0.60%     -    0s\n",
      "H    0     0                       9.2767353    9.23370  0.46%     -    0s\n",
      "H    0     0                       9.2767088    9.23370  0.46%     -    0s\n",
      "     0     2    9.23370    0  159    9.27671    9.23370  0.46%     -    0s\n",
      "  4352  1739    9.25589   27  142    9.27671    9.24045  0.39%   7.8    5s\n",
      " 10706  6021    9.27656   70   71    9.27671    9.24342  0.36%   6.2   10s\n",
      " 16173  9505    9.26145   29  135    9.27671    9.24563  0.33%   5.8   15s\n",
      " 23137 13316     cutoff   46         9.27671    9.24657  0.32%   5.7   20s\n",
      " 30342 17904    9.25878   29  137    9.27671    9.24820  0.31%   5.5   25s\n",
      " 37418 21729    9.26809   38  120    9.27671    9.24887  0.30%   5.5   30s\n",
      " 43300 25398    9.25220   23  145    9.27671    9.24925  0.30%   5.5   35s\n",
      " 49923 28649     cutoff   26         9.27671    9.24974  0.29%   5.5   40s\n",
      " 56366 32058    9.25327   26  142    9.27671    9.25011  0.29%   5.6   45s\n",
      " 63168 35828    9.26586   37  122    9.27671    9.25072  0.28%   5.6   50s\n",
      " 70309 39563     cutoff   50         9.27671    9.25134  0.27%   5.6   55s\n",
      " 77511 43487    9.26526   39  118    9.27671    9.25167  0.27%   5.6   60s\n",
      " 84385 47064    9.27477   50  103    9.27671    9.25201  0.27%   5.6   65s\n",
      " 89428 49288    9.27333   45  111    9.27671    9.25219  0.26%   5.6   70s\n",
      "\n",
      "Explored 93052 nodes (524227 simplex iterations) in 72.89 seconds\n",
      "Thread count was 4 (of 4 available processors)\n",
      "\n",
      "Solution count 4: 9.27671 9.27674 9.28918 9.32372 \n",
      "\n",
      "Solve interrupted\n",
      "Best objective 9.276708834666e+00, best bound 9.252337720840e+00, gap 0.2627%\n",
      "\n",
      "User-callback calls 193009, time in user-callback 0.06 sec\n",
      " 73.451961 seconds (1.39 M allocations: 66.923 MiB, 0.04% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(9.276708834665603, [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0115632  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.00842449])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time primal_obj, β_primal = compute_primal(X_p, y, k, γ, σ_X_p, TimeLimit=120, LogFile=\"debug.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Set Gurobi TimeLimit.\n",
      "LogFile: debug.txt.\n",
      "Gurobi Optimizer version 9.0.3 build v9.0.3rc0 (mac64)\n",
      "Optimize a model with 102 rows, 201 columns and 601 nonzeros\n",
      "Model fingerprint: 0xa89f65ed\n",
      "Variable types: 1 continuous, 200 integer (200 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [4e-07, 1e+00]\n",
      "  Objective range  [1e+00, 1e+00]\n",
      "  Bounds range     [9e+00, 9e+00]\n",
      "  RHS range        [4e-04, 1e+01]\n",
      "Presolve time: 0.00s\n",
      "Presolved: 102 rows, 201 columns, 601 nonzeros\n",
      "Variable types: 1 continuous, 200 integer (200 binary)\n",
      "\n",
      "Root relaxation: objective 9.275754e+00, 22 iterations, 0.00 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0    9.27612    0    2          -    9.27612      -     -    0s\n",
      "     0     2    9.27612    0    2          -    9.27612      -     -    0s\n",
      "H    6     7                       9.2897739    9.27667  0.14%   2.7    0s\n",
      "*   15     5               7       9.2767353    9.27670  0.00%   2.6    0s\n",
      "\n",
      "Cutting planes:\n",
      "  Lazy constraints: 10\n",
      "\n",
      "Explored 19 nodes (71 simplex iterations) in 0.30 seconds\n",
      "Thread count was 4 (of 4 available processors)\n",
      "\n",
      "Solution count 2: 9.27674 9.28977 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Best objective 9.276735332741e+00, best bound 9.276704480677e+00, gap 0.0003%\n",
      "\n",
      "User-callback calls 116, time in user-callback 0.26 sec\n",
      "  1.364038 seconds (2.62 M allocations: 125.274 MiB, 6.32% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(9.276735332741374, [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.011561  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.00842449])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time dual_obj, β_dual = compute_dual(X_p, y, k, γ, σ_X_p, TimeLimit=120, LogFile=\"debug.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100×3 Array{Float64,2}:\n",
       " 0.0  0.0         0.0       \n",
       " 0.0  0.0         0.0       \n",
       " 0.0  0.0         0.0       \n",
       " 0.0  0.0         0.0       \n",
       " 0.0  0.0         0.0       \n",
       " 0.0  0.0         0.0       \n",
       " 0.0  0.0         0.0       \n",
       " 0.0  0.0         0.0       \n",
       " 0.0  0.0         0.0       \n",
       " 1.0  0.0115632   0.011561  \n",
       " 0.0  0.0         0.0       \n",
       " 0.0  0.0         0.0       \n",
       " 0.0  0.0         0.0       \n",
       " ⋮                          \n",
       " 0.0  0.0         0.0       \n",
       " 1.0  0.00842449  0.00842449\n",
       " 0.0  0.0         0.0       \n",
       " 0.0  0.0         0.0       \n",
       " 0.0  0.0         0.0       \n",
       " 0.0  0.0         0.0       \n",
       " 0.0  0.0         0.0       \n",
       " 0.0  0.0         0.0       \n",
       " 0.0  0.0         0.0       \n",
       " 0.0  0.0         0.0       \n",
       " 0.0  0.0         0.0       \n",
       " 1.0  0.00842449  0.00842449"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hcat(β_true, β_primal, β_dual)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_____\n",
    "## 5. Experiences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "300"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Parameters\n",
    "n_train = 10000\n",
    "n = 2*n_train\n",
    "NR = 0.05\n",
    "α = 0.05;\n",
    "t_max = 5*60 # 30 min max per solving"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "csv_path = \"results/2021_03_09-19h00.csv\"\n",
    "logfile = \"results/2021_03_09-19h00.log\"\n",
    "\n",
    "write_list(csv_path, [\"Algo\", \"Seed\", \"n\", \"p\", \"k_true\", \"k\", \"γ\", \n",
    "                      \"NR\", \"α\", \"obj\", \"R2\", \"OR2\", \"t_algo\", \"t_data\", \"t_variance\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Set Gurobi TimeLimit.\n",
      "LogFile: results/2021_03_09-19h00.log.\n",
      "Gurobi Optimizer version 9.0.3 build v9.0.3rc0 (mac64)\n",
      "Optimize a model with 102 rows, 201 columns and 601 nonzeros\n",
      "Model fingerprint: 0x1b31c87a\n",
      "Variable types: 1 continuous, 200 integer (200 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e-04, 3e+01]\n",
      "  Objective range  [1e+00, 1e+00]\n",
      "  Bounds range     [1e+04, 1e+04]\n",
      "  RHS range        [1e+00, 3e+01]\n",
      "Presolve time: 0.00s\n",
      "Presolved: 102 rows, 201 columns, 601 nonzeros\n",
      "Variable types: 1 continuous, 200 integer (200 binary)\n",
      "\n",
      "Root relaxation: objective 1.459672e+04, 20 iterations, 0.00 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0 14663.8707    0    2          - 14663.8707      -     -    0s\n",
      "H    0     0                    14850.270097 14663.8707  1.26%     -    0s\n",
      "     0     2 14696.6736    0   10 14850.2701 14696.6736  1.03%     -    0s\n",
      "* 1030   817              32    14815.428062 14714.0690  0.68%   5.4    3s\n",
      "  1302  1051 14782.2829   64    - 14815.4281 14714.6009  0.68%   5.8    5s\n",
      "* 1344  1011              49    14812.908955 14714.6009  0.66%   5.6    5s\n",
      "* 1807  1057              39    14768.950108 14715.2133  0.36%   5.5    6s\n",
      "* 2174  1224              71    14768.560182 14715.3574  0.36%   5.6    6s\n",
      "* 2893  1830             192    14767.378841 14715.8405  0.35%   5.5    8s\n",
      "* 2896  1832             189    14767.355501 14715.8405  0.35%   5.5    8s\n",
      "* 2897  1832             188    14767.205496 14715.8405  0.35%   5.5    8s\n",
      "* 3153  2043             194    14767.044737 14715.9484  0.35%   5.3    9s\n",
      "* 3155  2043             190    14766.946613 14715.9484  0.35%   5.3    9s\n",
      "  3693  2543 14756.8119  195    2 14766.9466 14716.0598  0.34%   5.2   10s\n",
      "* 4402  3122             201    14766.935204 14716.2641  0.34%   5.1   11s\n",
      "* 5298  3763             203    14760.257247 14716.4355  0.30%   5.0   14s\n",
      "  5671  4099     cutoff   47      14760.2572 14716.5042  0.30%   5.2   15s\n",
      "  7071  5323     cutoff  123      14760.2572 14716.8864  0.29%   5.1   20s\n",
      "  8973  6932     cutoff  144      14760.2572 14717.1705  0.29%   5.1   26s\n",
      "  9970  7836 14723.9027   46   11 14760.2572 14717.3016  0.29%   5.1   31s\n",
      " 11359  8943 14739.4287   95    3 14760.2572 14717.5261  0.29%   5.2   36s\n",
      " 12725 10169 14732.8527  116    8 14760.2572 14717.7353  0.29%   5.2   40s\n",
      " 14179 11381     cutoff   64      14760.2572 14717.9536  0.29%   5.2   45s\n",
      " 15645 12490     cutoff   99      14760.2572 14718.0456  0.29%   5.1   50s\n",
      " 18448 14474 14732.6459   47   24 14760.2572 14718.3791  0.28%   5.1   55s\n",
      " 20293 16367     cutoff   85      14760.2572 14718.6235  0.28%   5.1   60s\n",
      " 21566 17432 14731.7166   51    5 14760.2572 14718.7513  0.28%   5.1   65s\n",
      " 24331 19730     cutoff   68      14760.2572 14719.0092  0.28%   5.1   70s\n",
      " 27347 22269     cutoff  160      14760.2572 14719.1705  0.28%   5.1   75s\n",
      " 29520 24137     cutoff   35      14760.2572 14719.3563  0.28%   5.1   81s\n",
      " 30657 24975 14742.4434  139    5 14760.2572 14719.4496  0.28%   5.0   85s\n",
      " 31611 25400     cutoff  114      14760.2572 14719.5082  0.28%   5.0   92s\n",
      " 31854 25801     cutoff   44      14760.2572 14719.5438  0.28%   5.1   98s\n",
      " 32198 26187     cutoff  107      14760.2572 14719.5801  0.28%   5.1  101s\n",
      "*33075 26512              42    14759.117955 14719.6316  0.27%   5.1  104s\n",
      " 33319 26993     cutoff  110      14759.1180 14719.6815  0.27%   5.1  106s\n",
      " 34797 28040     cutoff   82      14759.1180 14719.7688  0.27%   5.1  111s\n",
      " 36001 29207 14730.4435   49   19 14759.1180 14719.8327  0.27%   5.1  116s\n",
      " 37497 30283     cutoff  173      14759.1180 14719.9590  0.27%   5.1  120s\n",
      " 38802 31581 14721.1709   31   29 14759.1180 14720.0338  0.26%   5.1  125s\n",
      "*39992 31766              37    14757.583916 14720.1182  0.25%   5.1  129s\n",
      " 40211 32059     cutoff  155      14757.5839 14720.1743  0.25%   5.1  130s\n",
      " 41720 33486 14730.6085   70   13 14757.5839 14720.2983  0.25%   5.1  135s\n",
      " 43503 34803 14756.7138  201    1 14757.5839 14720.4310  0.25%   5.1  141s\n",
      " 45090 36232 14734.2341   67    6 14757.5839 14720.4908  0.25%   5.1  146s\n",
      " 46070 36833     cutoff   63      14757.5839 14720.5609  0.25%   5.1  150s\n",
      " 47307 37740     cutoff  113      14757.5839 14720.6665  0.25%   5.1  156s\n",
      " 48212 38676     cutoff  107      14757.5839 14720.7242  0.25%   5.2  161s\n",
      " 49309 39440     cutoff  109      14757.5839 14720.7941  0.25%   5.2  165s\n",
      " 50908 40644 14731.0383   55    2 14757.5839 14720.8776  0.25%   5.2  172s\n",
      " 51882 41615 14736.9369  100   10 14757.5839 14720.9245  0.25%   5.2  176s\n",
      " 52958 42501 14730.9696   41   10 14757.5839 14720.9683  0.25%   5.2  181s\n",
      " 53807 43178 14732.4730   66   17 14757.5839 14720.9885  0.25%   5.2  186s\n",
      " 54722 43954     cutoff  148      14757.5839 14721.0646  0.25%   5.2  191s\n",
      " 55669 44811     cutoff  109      14757.5839 14721.0924  0.25%   5.2  196s\n",
      " 57259 46094     cutoff  108      14757.5839 14721.1691  0.25%   5.2  201s\n",
      " 58876 47322 14750.5577  110    - 14757.5839 14721.2292  0.25%   5.2  206s\n",
      " 60299 48533 14729.2934   72   15 14757.5839 14721.2811  0.25%   5.2  211s\n",
      " 61752 49740     cutoff  165      14757.5839 14721.3279  0.25%   5.2  216s\n",
      " 62731 50723 14739.0100   98   12 14757.5839 14721.3523  0.25%   5.2  221s\n",
      " 63953 51556     cutoff  137      14757.5839 14721.3994  0.25%   5.2  227s\n",
      " 64761 52384 14740.0351   76    7 14757.5839 14721.4306  0.24%   5.2  232s\n",
      " 65394 52873     cutoff  152      14757.5839 14721.4493  0.24%   5.2  235s\n",
      " 66565 53676     cutoff  121      14757.5839 14721.4918  0.24%   5.2  241s\n",
      " 67452 54555     cutoff   64      14757.5839 14721.5241  0.24%   5.2  248s\n",
      " 67974 55071 14742.1474   96    5 14757.5839 14721.5407  0.24%   5.2  250s\n",
      " 69022 55867     cutoff  124      14757.5839 14721.5896  0.24%   5.2  256s\n",
      " 69845 56460 14733.9167   85    7 14757.5839 14721.6042  0.24%   5.3  263s\n",
      " 70298 56931     cutoff  129      14757.5839 14721.6212  0.24%   5.3  266s\n",
      " 70866 57397     cutoff  139      14757.5839 14721.6302  0.24%   5.3  270s\n",
      " 71721 58086     cutoff   68      14757.5839 14721.6758  0.24%   5.3  277s\n",
      " 72236 58430 14739.2095   74    4 14757.5839 14721.6900  0.24%   5.3  280s\n",
      " 73168 59083     cutoff  139      14757.5839 14721.7407  0.24%   5.3  286s\n",
      " 73974 59868     cutoff  146      14757.5839 14721.7656  0.24%   5.3  291s\n",
      " 74429 60337     cutoff  148      14757.5839 14721.7965  0.24%   5.3  296s\n",
      " 75053 60767     cutoff  118      14757.5839 14721.8241  0.24%   5.3  300s\n",
      "\n",
      "Cutting planes:\n",
      "  Lazy constraints: 1103\n",
      "\n",
      "Explored 75528 nodes (402979 simplex iterations) in 300.05 seconds\n",
      "Thread count was 4 (of 4 available processors)\n",
      "\n",
      "Solution count 10: 14757.6 14759.1 14760.3 ... 14768.6\n",
      "\n",
      "Time limit reached\n",
      "Best objective 1.475758391576e+04, best bound 1.472183380666e+04, gap 0.2422%\n",
      "\n",
      "User-callback calls 160979, time in user-callback 6.28 sec\n",
      "Set Gurobi TimeLimit.\n",
      "LogFile: results/2021_03_09-19h00.log.\n"
     ]
    }
   ],
   "source": [
    "for seed ∈ [35]\n",
    "    for p ∈ [100, 500, 1000, 2500, 5000]\n",
    "\n",
    "        # True sparsity\n",
    "        k_true = 10\n",
    "\n",
    "        # Generate data\n",
    "        t_data = @elapsed (X_train, y_train), _, (X_test, y_test), β_true = generate_synthetic_data(n, p, k_true, NR, seed=seed);\n",
    "\n",
    "        # Significance\n",
    "        t_variance = @elapsed t_α = get_t_α(n_train, p, α)\n",
    "        \n",
    "        # Robustness\n",
    "        for γ ∈ [p/sqrt(n_train)]\n",
    "            \n",
    "            # Variance\n",
    "            t_variance += @elapsed σ_X = t_α * get_σ_X(X_train, y_train, γ)\n",
    "        \n",
    "            # Estimated sparsity\n",
    "            Random.seed!(seed)\n",
    "            k = 12\n",
    "            \n",
    "            # Dual\n",
    "            t_dual = @elapsed obj_dual, β_dual = compute_dual(X_train, y_train, k, γ, σ_X; LogFile=logfile, TimeLimit=t_max);\n",
    "            R2_dual = get_R2(X_train*β_dual, y_train, y_train)\n",
    "            OR2_dual = get_R2(X_test*β_dual, y_test, y_train)\n",
    "            write_list(csv_path, [\"dual\", seed, n_train, p, k_true, k, γ, \n",
    "                                  NR, α, obj_dual, R2_dual, OR2_dual, t_dual, t_data, t_variance])\n",
    "            \n",
    "            # Primal\n",
    "            t_primal = @elapsed obj_primal, β_primal = compute_primal(X_train, y_train, k, γ, σ_X; LogFile=logfile, TimeLimit=t_max);\n",
    "            R2_primal = get_R2(X_train*β_primal, y_train, y_train)\n",
    "            OR2_primal = get_R2(X_test*β_primal, y_test, y_train)\n",
    "            write_list(csv_path, [\"primal\", seed, n_train, p, k_true, k, γ, \n",
    "                                  NR, α, obj_primal, R2_primal, OR2_primal, t_primal, t_data, t_variance])\n",
    "            \n",
    "        end\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.0.5",
   "language": "julia",
   "name": "julia-1.0"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.0.5"
  },
  "toc-showtags": false
 },
 "nbformat": 4,
 "nbformat_minor": 4
}