%%% Produces six LaTeX table files: quadratic, logistic regression,
%%% elastic-net regression (ML1), group-sparse Huberized SVM (ML2),
%%% Bayesian D-optimal experimental design (ML3), and Huberized quantum
%%% state tomography (quantSDP4).

%% Initialize
clc
clear
close all

addpath(genpath('../..'));

%% Parameters for 3 Algorithms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% OIAL / APF-IAL RFISTA inputs
PAROIAL.chi = 0.001;
PAROIAL.beta = 1.25;

%%% Additional RFISTA inputs
PAROIAL.mu = inf;
PAROIAL.Mbar0 = inf;
PAROIAL.Munderbar = inf;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

PAROIAL.rhohat = 1e-5;
PAROIAL.epsilonhat = 1e-5;
PAROIAL.verbose = true;
PAROIAL.inner_verbose = false;
PAROIAL.print_every = 100;

PARAPF = PAROIAL;

PAROIAL.c = 10;

PARAPF.c = 10;
PARAPF.tildepsilon = 10;
PARAPF.alpha = 1.1;
PARAPF.omega = 1.5;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Algorithm 5 PAL inputs (arXiv:2206.01209), specialized to A*x=b
PARALG5.verbose = true;
PARALG5.inner_verbose = false;
PARALG5.print_every = 100;
PARALG5.epsilonhat = PAROIAL.epsilonhat;

% Paper parameters: rho_k = rho0*zeta^k, eta_k = eta0*sigma^k.
PARALG5.mu0 = Inf;
PARALG5.rho0 = 10;
PARALG5.rho0place = 10;
PARALG5.eta0 = 10;
PARALG5.zeta = 1.1;
PARALG5.sigma = 1 / (PARALG5.zeta + 0.4);
PARALG5.delta = 0.8;
PARALG5.alpha0 = 1;
PARALG5.M_alg2 = 10;
PARALG5.maxiter = 1e6;

%% 1. Quadratic vector problem
Nvec = [1000, 1000, 1000, 1000, 1000, ...
    500, 500, 500, 500, 500, 200, 200, 200, 200, 200];
Mvec = Nvec;
mvec = [10, 10, 10, 10, 10, 100, 100, 100, 100, 100, ...
    5, 5, 5, 5, 5];
seedvec = (1:length(Nvec)) * 50;
dimMvec = [10, 50, 100, 150, 200, 250, 300, 350, 400, 450, ...
    500, 550, 600, 650, 700];
dimNvec = [50, 125, 200, 275, 350, 425, 500, 575, 650, 725, ...
    800, 875, 950, 1025, 1100];

filename = 'table_QuadraticLC.tex';
options.problemclass = "quadratic";

PAROIAL.D = sqrt(2);
PAROIAL.epsilon = min( ...
    (PAROIAL.c * (PAROIAL.epsilonhat)^2) / (4 * PAROIAL.D), ...
    PAROIAL.rhohat);


runAL(PAROIAL, PARAPF, PARALG5, Nvec, Mvec, mvec, ...
    seedvec, dimMvec, dimNvec, options, filename);
% 
%% 2. Logistic regression over an L1 ball
Nvec = [50, 75, 100, 150, 175, 200, 250, 300, 350, 400, ...
    450, 500, 550, 600, 650];
Mvec = [30, 50, 80, 120, 140, 200, 250, 300, 350, 400, ...
    450, 500, 600, 800, 1000];
mvec = [0, 0, 0, 0, 0, 0, 1e-4, ...
    1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 1e-1, 5e0];
seedvec = (1:length(Nvec)) * 50;
dimMvec = [50, 100, 125, 150, 175, 200, 250, 300, 350, 400, ...
    450, 500, 600, 800, 1000];
dimNvec = [100, 150, 200, 300, 400, 500, 750, 1000, 1200, 1400, ...
    1500, 1750, 2000, 2200, 2500];

filename = 'table_LogReg.tex';
options.problemclass = "logreg";

PAROIAL.D = 2; % diameter of the unit L1 ball
PAROIAL.epsilon = min( ...
    (PAROIAL.c * (PAROIAL.epsilonhat)^2) / (4 * PAROIAL.D), ...
    PAROIAL.rhohat);

runAL(PAROIAL, PARAPF, PARALG5, Nvec, Mvec, mvec, ...
    seedvec, dimMvec, dimNvec, options, filename);

%% 3. ML1: elastic-net least-squares regression
% A sample-rich regression grid. The smaller instances are deliberately
% more ill-conditioned, while the larger dense instances are better
% conditioned to keep the expected runtime moderate.
Nvec = [300, 360, 430, 500, 580, 670, 760, 860, 970, 1090, ...
    1220, 1360, 1400, 1480, 1560];
Mvec = [15, 22, 32, 45, 62, 85, 110, 140, 180, 225, ...
    280, 345, 380, 410, 480];
mvec = [0, 0, 0, 0, 0, 0, 1e-4, ...
    5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 1e-1, 5e-1, 5e0];
seedvec = 1103 + (0:(length(Nvec)-1)) * 79;
dimMvec = [50, 100, 140, 180, 200, 250, 350, 400, 450, 500, ...
    550, 650, 750, 850, 1000];
dimNvec = [100, 125, 200, 250, 280, 300, 400, 450, 520, 600, ...
    700, 850, 1000, 1200, 1500];

filename = 'table_ML1_ElasticNet.tex';
options.problemclass = "ML1";

PAROIAL.D = 1;
PAROIAL.epsilon = min( ...
    (PAROIAL.c * (PAROIAL.epsilonhat)^2) / (4 * PAROIAL.D), ...
    PAROIAL.rhohat);

runAL(PAROIAL, PARAPF, PARALG5, Nvec, Mvec, mvec, ...
    seedvec, dimMvec, dimNvec, options, filename);

% %% 4. ML2: group-sparse Huberized-hinge SVM
% Uses data-derived matched-cohort fairness constraints on average scores.
% A distinct classification grid with different sample sizes, dimensions,
% curvature pairs, constraint counts, and random seeds from both logreg and
% ML1. As above, the largest dense instances have milder conditioning.
Nvec = [260, 330, 410, 500, 600, 710, 830, 1500, 2500, 3200, ...
    4300, 4500, 5000, 5800, 7000];
Mvec = [20, 30, 45, 65, 90, 120, 155, 200, 255, 320, ...
    395, 480, 575, 680, 800];
mvec = [0, 0, 0, 0, 2e-3, 5e-3, 1e-2, ...
    2.5e-2, 5e-2, 1e-1, 2e-1, 4e-1, 8e-1, 1e0, 5e0];
seedvec = 2207 + (0:(length(Nvec)-1)) * 83;
dimMvec = [100, 125, 160, 220, 300, 340, 410, 500, 580, 630, ...
    720, 800, 1000, 1200, 2000];
dimNvec = [200, 300, 350, 420, 500, 650, 800, 1000, 1300, 1800, ...
    2200, 2500, 2800, 3500, 4000];

filename = 'table_ML2_HuberizedSVM.tex';
options.problemclass = "ML2";

PAROIAL.D = 1;
PAROIAL.epsilon = min( ...
    (PAROIAL.c * (PAROIAL.epsilonhat)^2) / (4 * PAROIAL.D), ...
    PAROIAL.rhohat);

runAL(PAROIAL, PARAPF, PARALG5, Nvec, Mvec, mvec, ...
    seedvec, dimMvec, dimNvec, options, filename);

%% 5. ML3: Bayesian D-optimal experimental design over a box
% Uses one total-budget constraint and sparse category-budget constraints.
% Here N is the dimension of the Bayesian information matrix, dimN is the
% number of candidate experiments/design weights, and dimM is the number
% of linear resource constraints. All six input vectors differ from those
% used for the other four problem classes.
Nvec = [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, ...
    13, 14, 15, 16, 17];
Mvec = [5, 8, 12, 15, 20, 25, 30, 40, 50, 60, ...
    70, 80, 85, 90, 100];
mvec = [0, 0, 0, 0, 0, 0, 0, ...
    0, 5e-5, 1e-5, 5e-4, 1e-4, 5e-3, 1e-3, 5e-2];
seedvec = 5101 + (0:(length(Nvec)-1)) * 101;
dimMvec = [5, 10, 20, 30, 50, 60, 75, 85, 100, 120, 130, 140, 160, ...
    180 200];
dimNvec = [25, 30, 50, 70, 80, 90, 115, 125, 175, 200, ...
    220, 230, 250, 275, 300];

filename = 'table_ML3_BayesianDOptimal.tex';
options.problemclass = "ML3";

PAROIAL.D = 2; % diameter of [0,2/sqrt(dimN)]^dimN
PAROIAL.epsilon = min( ...
    (PAROIAL.c * (PAROIAL.epsilonhat)^2) / (4 * PAROIAL.D), ...
    PAROIAL.rhohat);

runAL(PAROIAL, PARAPF, PARALG5, Nvec, Mvec, mvec, ...
    seedvec, dimMvec, dimNvec, options, filename);


%% 6. quantSDP4: Huberized quantum state tomography
% Here N is the number of measured observables, dimN=d^2 is the real
% dimension of the d-by-d Hermitian matrix variable, and dimM is the number
% of exact linear expectation constraints.
Nvec = [60, 75, 95, 120, 150, 180, 215, 250, 285, 345, ...
    425, 515, 660, 800, 970];
Mvec = [14, 21, 31, 45, 63, 86, 115, 150, 192, 242, ...
    300, 368, 445, 532, 630];
mvec = [0, 0, 0, 0, 0, 1e-4, 5e-4, 1e-3, 5e-3, 1e-2, ...
    5e-2, 1e-1, 5e-1, 1e0, 5e0];
seedvec = 7301 + (0:(length(Nvec)-1)) * 107;
dimMvec = [8, 10, 12, 15, 18, 21, 24, 27, 30, 35, ...
    42, 50, 61, 73, 86];
dimNvec = [144, 196, 256, 324, 400, 484, 576, 676, 784, 1024, ...
    1296, 1600, 2304, 3136, 4096];

filename = 'table_quantSDP4_QuantumTomography.tex';
options.problemclass = "quantSDP4";

PAROIAL.D = sqrt(2); % diameter of the density-matrix spectrahedron
PAROIAL.epsilon = min( ...
    (PAROIAL.c * (PAROIAL.epsilonhat)^2) / (4 * PAROIAL.D), ...
    PAROIAL.rhohat);

runAL(PAROIAL, PARAPF, PARALG5, Nvec, Mvec, mvec, ...
    seedvec, dimMvec, dimNvec, options, filename);

% profile report
