benchmark_optimized_true.ipynb

{
 "cells": [
  {
   "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 10 days\n",
      "--------------------------------------------\n",
      "\n",
      "Academic license - for non-commercial use only\n"
     ]
    }
   ],
   "source": [
    "# Create a gurobi model without the annoying academic license message\n",
    "gurobi_env = Gurobi.Env()\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",
    "    #set_optimizer_attribute(model, \"Threads\", 4)\n",
    "    return model\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_____\n",
    "# Dual Holistic Regression | Benchmark"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Primal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "compute_primal (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function compute_primal(X, y, k, γ, t_α, σ_X, log_path)\n",
    "    \n",
    "    n, p = size(X)\n",
    "    \n",
    "    model = create_gurobi_model(LogFile=log_path)\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] >= t_α*s[i])\n",
    "    @constraint(model, [i=1:p], -β[i]/σ_X[i] + M2*(1-b[i]) >= t_α*s[i])\n",
    "\n",
    "    @objective(model, Min, 0.5*sum((y[i] - X[i,:]'β)^2 for i=1:n) + (0.5/γ)* sum(β[j]^2 for j=1:p))\n",
    "    JuMP.optimize!(model)\n",
    "    \n",
    "    return objective_value(model), value.(β)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "function g(s, X, y, D, Z, γ, t_α, σ_X, model; compute_β=false)\n",
    "    \n",
    "    zero_indexes = findall(x->x == 0, s)\n",
    "    \n",
    "    # Create DZ once\n",
    "    DZ = D*Z\n",
    "    \n",
    "    # Compute norm\n",
    "    function compute_DZ_square_norm(in_norm)\n",
    "        return in_norm' * DZ * in_norm\n",
    "    end\n",
    "    \n",
    "    λ = model[:λ]\n",
    "    obj_1 = 0.5*y'y\n",
    "    obj_2 = t_α*λ'*(s.*σ_X)\n",
    "    obj_3 = - 0.5 * compute_DZ_square_norm(X'y + λ)\n",
    "    @objective(model, Max, obj_2 + obj_3)\n",
    "    \n",
    "    JuMP.optimize!(model)\n",
    "    λ_value = value.(λ)  \n",
    "    obj_value = objective_value(model)\n",
    "\n",
    "    \n",
    "    β_pred = nothing\n",
    "    \n",
    "    if compute_β\n",
    "        \n",
    "        # Get size\n",
    "        n, p = size(X)\n",
    "        \n",
    "        # Compute β_s\n",
    "        sparsity_indexes = findall(x->x>0, s)\n",
    "        X_s = X[:, sparsity_indexes]\n",
    "        λ_s = λ_value[sparsity_indexes]\n",
    "        β_s = inv(I/γ + X_s'X_s)*(X_s'y + λ_s)\n",
    "\n",
    "        # Compute β\n",
    "        β_pred = zeros(p)\n",
    "        β_pred[sparsity_indexes] = β_s\n",
    "\n",
    "    end\n",
    "    \n",
    "    return β_pred, λ_value, obj_1 + obj_value\n",
    "    \n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_____\n",
    "## Dual Optimized"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "g_s (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function g_s(supp_ids, obj1, Xty_s, D_s, σ_X_s, M_s, γ, t_α)\n",
    "    \n",
    "    l = length(supp_ids)\n",
    "    \n",
    "    if l==0\n",
    "        return zeros(0), 0.0\n",
    "    end\n",
    "    \n",
    "    lower = [-Inf for _ in 1:l]\n",
    "    upper = zeros(l)\n",
    "    \n",
    "    #initial_λ = zeros(l) .- 1.0\n",
    "    initial_λ = min.(-1.0, -t_α*(I/γ + M_s)*σ_X_s + Xty_s)\n",
    "    #println(initial_λ)\n",
    "    \n",
    "    \n",
    "    D_sXty_s = D_s*Xty_s\n",
    "    \n",
    "    f(λ) = -t_α*λ'σ_X_s + 0.5(Xty_s+λ)'*(D_sXty_s+D_s*λ)\n",
    "    \n",
    "    ∇f = function(λ)\n",
    "        return -t_α*σ_X_s + D_sXty_s + D_s*λ\n",
    "    end\n",
    "    \n",
    "    #println(\"Init \", initial_λ, f(initial_λ), ∇f(initial_λ))\n",
    "    \n",
    "    inner_optimizer = GradientDescent()\n",
    "    res = Optim.optimize(f, ∇f, lower, upper, initial_λ, Fminbox(inner_optimizer), inplace=false)\n",
    "\n",
    "    output =  - res.minimizer, obj1 - res.minimum\n",
    "    \n",
    "    #println(\"output\", output)\n",
    "    \n",
    "    return output\n",
    "    \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "∇g_s (generic function with 2 methods)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function ∇g_s(supp_ids, supp_c_ids, Xty_s, Xty_s_c, D_s, σ_X_s, M_s_c, λ_s, γ, p)\n",
    "    \n",
    "    μ_s = Xty_s + λ_s\n",
    "    μ_s_c = Xty_s_c\n",
    "    \n",
    "    D_stμ_s = D_s'μ_s\n",
    "    \n",
    "    grad = zeros(2p)    \n",
    "    grad[supp_ids] = λ_s .* σ_X_s - 0.5(D_stμ_s).^2\n",
    "    grad[supp_c_ids] = - 0.5*γ^2*(M_s_c'D_stμ_s + μ_s_c).^2\n",
    "    \n",
    "    return grad\n",
    "    \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "get_support_ids (generic function with 1 method)"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function get_support_ids(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] == zero(eltype(s))\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([3, 4, 6], [1, 2, 5])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_support_ids([0, 0, 4, 4, 0, 4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "project (generic function with 1 method)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function project(s, X, Xty, σ_X, γ)\n",
    "    supp_ids, supp_c_ids = get_support_ids(s)\n",
    "    X_s = X[:, supp_ids]\n",
    "    X_s_c = X[:, supp_c_ids]\n",
    "    M_s = X_s'X_s\n",
    "    M_s_c = X_s'X_s_c\n",
    "    D_s = inv(I/γ + M_s)\n",
    "    Xty_s = Xty[supp_ids]\n",
    "    Xty_s_c = Xty[supp_c_ids]\n",
    "    σ_X_s = σ_X[supp_ids]\n",
    "    return supp_ids, supp_c_ids, D_s, Xty_s, Xty_s_c, σ_X_s, M_s, M_s_c\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_____\n",
    "## Dual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "compute_dual (generic function with 1 method)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function compute_dual(X_p, y, k, γ, t_α, σ_X_p, log_path)\n",
    "    n, p = size(X_p)\n",
    "    \n",
    "    # Extended Matrices\n",
    "    X = hcat(X_p, -X_p)\n",
    "    σ_X = [σ_X_p; σ_X_p]\n",
    "    Xty = X'y \n",
    "    obj1 = 0.5*y'y\n",
    "    \n",
    "    # Outer problem\n",
    "    miop = create_gurobi_model(LogFile=log_path)\n",
    "    \n",
    "    @variable(miop, s[1:2*p], Bin)\n",
    "    @variable(miop, t >= 0)\n",
    "    \n",
    "    @constraint(miop, sum(s) <= k)\n",
    "    @constraint(miop, [i=1:p], s[i]+s[p+i]<=1)\n",
    "    \n",
    "    # --- Cutting plane --- #\n",
    "    \n",
    "#    inner_model = create_gurobi_model()\n",
    "#    @variable(inner_model, λ[j=1:2k] >= 0)\n",
    "    \n",
    "    # Initial solution\n",
    "    s_init = zeros(2*p) #TODO: change this\n",
    "    s_init[1:k] .= 1\n",
    "    \n",
    "    supp_ids, supp_c_ids, D_s, Xty_s, Xty_s_c, σ_X_s, M_s, M_s_c = project(s_init, X, Xty, σ_X, γ)\n",
    "    \n",
    "    λ_s_init, g_init = g_s(supp_ids, obj1, Xty_s, D_s, σ_X_s, M_s, γ, t_α)\n",
    "    ∇g_init = ∇g_s(supp_ids, supp_c_ids, Xty_s, Xty_s_c, D_s, σ_X_s, M_s_c, λ_s_init, γ, p)\n",
    "    \n",
    "    @constraint(miop, t >= g_init + dot(∇g_init, s - s_init))\n",
    "    @objective(miop, Min, t)\n",
    "    \n",
    "    \n",
    "    function outer_approximation(cb_data)\n",
    "        \n",
    "        s_val = [callback_value(cb_data, s[i]) for i=1:2p]\n",
    "\n",
    "        supp_ids, supp_c_ids, D_s, Xty_s, Xty_s_c, σ_X_s, M_s, M_s_c = project(s_val, X, Xty, σ_X, γ)\n",
    "        \n",
    "        λ_s_val, g_val = g_s(supp_ids, obj1, Xty_s, D_s, σ_X_s, M_s, γ, t_α)\n",
    "        ∇g_val = ∇g_s(supp_ids, supp_c_ids, Xty_s, Xty_s_c, D_s, σ_X_s, M_s_c, λ_s_val, γ, p)\n",
    "        \n",
    "        con = @build_constraint(t >= g_val + dot(∇g_val, s - s_val))\n",
    "        MOI.submit(miop, MOI.LazyConstraint(cb_data), con)\n",
    "        \n",
    "    end\n",
    "    \n",
    "    MOI.set(miop, MOI.LazyConstraintCallback(), outer_approximation)\n",
    "    JuMP.optimize!(miop)\n",
    "    \n",
    "    s_opt = JuMP.value.(s)\n",
    "\n",
    "    supp_ids, supp_c_ids, D_s, Xty_s, Xty_s_c, σ_X_s, M_s, M_s_c = project(s_opt, X, Xty, σ_X, γ)\n",
    "    λ_s_opt, _ = g_s(supp_ids, obj1, Xty_s, D_s, σ_X_s, M_s, γ, t_α)\n",
    "    \n",
    "    β_opt = zeros(2p)\n",
    "    β_opt[supp_ids] = D_s*(Xty_s + λ_s_opt)\n",
    "    β_pred = β_opt[1:p] .- β_opt[p+1:end]\n",
    "    \n",
    "    return objective_value(miop), β_pred \n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "____\n",
    "## Synthetic Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "generate_synthetic_data (generic function with 1 method)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function generate_synthetic_data(n, p, k, NR)\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "get_t_α_and_σ_X (generic function with 1 method)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function get_t_α_and_σ_X(X, y, α, γ)\n",
    "    n, p = size(X)\n",
    "    \n",
    "    # Student law\n",
    "    t_α = quantile(TDist(n-p), 1 - α/2) # Beware: n-p-1 if we add intercept\n",
    "    \n",
    "    # Estimator σ\n",
    "    M = 1/γ*I + X'X\n",
    "    M_inv = inv(M)\n",
    "    σ_tilde = sqrt((y'*(I - X*M_inv*X')*y)/(n-p))\n",
    "    σ_X = σ_tilde * sqrt.(diag(M_inv))\n",
    "    \n",
    "    return t_α, σ_X\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "get_OR2 (generic function with 1 method)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function get_insample_R2(y_pred, y_true, p)\n",
    "    TSE = sum((y_pred[i]-y_true[i])^2 for i=1:p)\n",
    "    baseline_E = sum((sum(y_true)/length(y_true)-y_true[i])^2 for i=1:p)\n",
    "    return 1 - TSE/baseline_E\n",
    "end\n",
    "\n",
    "function get_OR2(y_pred, y_true, y_train, p)\n",
    "    TSE = sum((y_pred[i]-y_true[i])^2 for i=1:p)\n",
    "    baseline_E = sum((sum(y_train)/length(y_train)-y_true[i])^2 for i=1:p)\n",
    "    return 1 - TSE/baseline_E\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "____\n",
    "## Experiences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "list2string (generic function with 1 method)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function r(x,dig)\n",
    "    return round(x,digits = dig)\n",
    "end\n",
    "\n",
    "function list2string(list)\n",
    "    s=\"\"\n",
    "    for x in list \n",
    "        s*=x\n",
    "    end\n",
    "    s*=\"\\n\"\n",
    "    return s\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "write_to_file (generic function with 1 method)"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function write_to_file(file_path, str)\n",
    "    open(file_path, \"a+\") do io\n",
    "        write(io, str)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "NR = 0.001\n",
    "α = 0.05\n",
    "n_train = 10000\n",
    "n = 2*n_train\n",
    "p = 300\n",
    "k = 10\n",
    "γ = 1\n",
    "\n",
    "Random.seed!(2021)\n",
    "t_synthetic = @elapsed (X_p, y), (X_val, y_val), (X_test, y_test), β_true = generate_synthetic_data(n, p, k, NR);\n",
    "t_thresh_var = @elapsed t_α, σ_X_p = get_t_α_and_σ_X(X_p, y, α, γ)\n",
    "\n",
    "X = hcat(X_p, -X_p)\n",
    "M = X'X\n",
    "Xty = X'y\n",
    "σ_X = [σ_X_p; -σ_X_p]\n",
    "obj1 = 0.5*y'y\n",
    "\n",
    "model_inner = create_gurobi_model()\n",
    "@variable(model_inner, λ[1:2*p] >= 0)\n",
    "\n",
    "# Initial solution\n",
    "\n",
    "total_t_g_s = 0\n",
    "total_t_g = 0\n",
    "\n",
    "for _ in tqdm(1:100)\n",
    "\n",
    "    \n",
    "    s_init = zeros(2p)\n",
    "    idx = sample(1:2p, k)\n",
    "    s_init[idx] .= 1;\n",
    "    \n",
    "    #s_init = zeros(2*p) #TODO: change this\n",
    "    #s_init[1:k] .= 1\n",
    "\n",
    "    # Dual \n",
    "    total_t_g_s += @elapsed begin\n",
    "        supp_ids, supp_c_ids, D_s, Xty_s, Xty_s_c, σ_X_s, M_s_c = project(s_init, X, Xty, σ_X, γ)\n",
    "        λ_s_init, g_init = g_s(supp_ids, obj1, Xty_s, D_s, σ_X_s, γ, t_α)\n",
    "    end\n",
    "    \n",
    "    total_t_g += @elapsed begin\n",
    "        Z_init = Diagonal(s_init)\n",
    "        D_init = γ*inv(I + γ * Z_init * M)\n",
    "        _ , λ_init, g_init = g(s_init, X, y, D_init, Z_init, γ, t_α, σ_X, model_inner)\n",
    "    end\n",
    "    \n",
    "end\n",
    "\n",
    "println(total_t_g_s)\n",
    "println(total_t_g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "NR = 0.001\n",
    "α = 0.05\n",
    "n_train = 10000\n",
    "n = 2*n_train\n",
    "\n",
    "file_path = \"results/results_arie_9.csv\"\n",
    "\n",
    "write_to_file(file_path, \"algo,seed,n,p,k,γ,NR,α,t_algo,R2,OR2,t_synthetic,t_thresh_var\\n\")\n",
    "\n",
    "for seed in [2021, 42, 25, 1998]\n",
    "    for γ in [1.0, 5.0, 10.0, 100.0]\n",
    "        for p in [10, 50, 70, 100, 150, 200, 250, 500]\n",
    "            \n",
    "            # Set k\n",
    "            k = Int(p/10)\n",
    "\n",
    "            # Generate data\n",
    "            Random.seed!(seed)\n",
    "            t_synthetic = @elapsed (X_train, y_train), (X_val, y_val), (X_test, y_test), β_true = generate_synthetic_data(n, p, k, NR);\n",
    "            t_thresh_var = @elapsed t_α, σ_X = get_t_α_and_σ_X(X_train, y_train, α, γ)\n",
    "\n",
    "            # Solving Dual\n",
    "            #log_path_dual = \"logs/dual_n=$(n)_p=$(p)_k=$(k)_gamma=$(γ)_seed=$(seed).txt\"\n",
    "            log_path_dual = \"debug.txt\"\n",
    "            t_dual = @elapsed obj_value, β_dual = compute_dual(X_train, y_train, k, γ, t_α, σ_X, log_path_dual);\n",
    "            R2_dual = get_insample_R2(X_train*β_dual, y_train, p)\n",
    "            OR2_dual = get_OR2(X_test*β_dual, y_test, y_train, p)\n",
    "            list_dual = [x*\",\" for x in string.([seed, n_train, p,k,γ,NR,α,\n",
    "                            r(t_dual,3),r(R2_dual,3),r(OR2_dual,3),\n",
    "                            r(t_synthetic,3),r(t_thresh_var,3)])]\n",
    "            s = \"dual,\"*list2string(list_dual) \n",
    "            write_to_file(file_path, s)\n",
    "\n",
    "            # Solving Primal\n",
    "            #log_path_primal = \"logs/primal_n=$(n)_p=$(p)_k=$(k)_gamma=$(γ)_seed=$(seed).txt\"\n",
    "            log_path_primal = \"debug.txt\"\n",
    "            t_primal = @elapsed obj_value, β_primal = compute_primal(X_train, y_train, k, γ, t_α, σ_X, log_path_primal)\n",
    "            R2_primal = get_insample_R2(X_train*β_primal, y_train, p)\n",
    "            OR2_primal = get_OR2(X_test*β_primal, y_test, y_train, p)\n",
    "            list_primal = [x*\",\" for x in string.([seed, n_train,p,k,γ,NR,α,\n",
    "                                r(t_primal,3),r(R2_primal,3),r(OR2_primal,3),\n",
    "                            r(t_synthetic,3),r(t_thresh_var,3)])]\n",
    "            s = \"primal,\"*list2string(list_primal) \n",
    "            write_to_file(file_path, s)\n",
    "        end\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "384px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}