function runAL(PAROIAL, PARAPF, PARALG5, Nvec, Mvec, mvec, seedvec, dimMvec, dimNvec, options,filename)
%RUNAL Compare OIAL, APFIAL, and Algorithm 5 PAL on generated problems.
%
% Supported problem classes:
%   "quadratic" : test_fn_lin_constr_01
%   "logreg"    : test_fn_lin_constr_logreg_l1_01
%   "ML1"       : test_fn_lin_constr_ML1
%   "ML2"       : test_fn_lin_constr_ML2
%   "ML3"       : test_fn_lin_constr_ML3
%   "quantSDP4" : test_fn_lin_quant_SDP4

numofProblems = length(Nvec);

TAPF = zeros(numofProblems, 1);
PFeasAPF = zeros(numofProblems,1);
DFeasAPF = zeros(numofProblems,1);

TOIAL = zeros(numofProblems, 1);
PFeasOIAL = zeros(numofProblems,1);
DFeasOIAL = zeros(numofProblems,1);

TALG5 = zeros(numofProblems, 1);
PFeasALG5 = zeros(numofProblems,1);
DFeasALG5 = zeros(numofProblems,1);

for j = 1:numofProblems
    N = Nvec(j);
    M = Mvec(j);
    m = mvec(j);
    seed = seedvec(j);
    dimM = dimMvec(j);
    dimN = dimNvec(j);

    fprintf('\n============================================================\n')
    fprintf('Problem %d/%d [%s]: dimM=%d, dimN=%d, seed=%d\n', ...
    j, numofProblems, char(options.problemclass), dimM, dimN, seed);
    fprintf('============================================================\n')

    if options.problemclass=="quadratic"
        [Finitial,PARinitial] = test_fn_lin_constr_01( ...
            N, M, m, seed, dimM, dimN);
    elseif options.problemclass=="logreg"
        [Finitial,PARinitial] = test_fn_lin_constr_logreg_l1_01( ...
            N, M, m, seed, dimM, dimN);
    elseif options.problemclass=="ML1"
        [Finitial,PARinitial] = test_fn_lin_constr_ML1( ...
            N, M, m, seed, dimM, dimN);
    elseif options.problemclass=="ML2"
        [Finitial,PARinitial] = test_fn_lin_constr_ML2( ...
            N, M, m, seed, dimM, dimN);
    elseif options.problemclass=="ML3"
        [Finitial,PARinitial] = test_fn_lin_constr_ML3( ...
            N, M, m, seed, dimM, dimN);
    elseif options.problemclass=="quantSDP4"
        [Finitial,PARinitial] = test_fn_lin_quant_SDP4( ...
            N, M, m, seed, dimM, dimN);
    else
        error('runAL:UnknownProblemClass', ...
            'Unknown options.problemclass. Use quadratic, logreg, ML1, ML2, ML3, or quantSDP4.');
    end

    x0 = PARinitial.x0;

    F = struct();
    F.prox = Finitial.prox;
    F.Residual = @(x) Finitial.Residual(x);
    F.Adjoint = @(x) Finitial.Adjoint(x);
    F.finit = @(x) Finitial.fs(x);
    F.gradinit = @(x) Finitial.grad(x);

    p0 = rand(dimM,1);

    PARAPFj = copyProblemParams(PARAPF,PARinitial);
    PAROIALj = copyProblemParams(PAROIAL,PARinitial);
    PARALG5j = copyProblemParams(PARALG5,PARinitial);

    PARAPFj.Mbar0 = PARAPFj.LipEst/4;
    PARAPFj.Munderbar = PARAPFj.LipEst/4;

    PAROIALj.Mbar0 = PAROIALj.LipEst/4;
    PAROIALj.Munderbar = PAROIALj.LipEst/4;

    PARALG5j.mu0 = PARALG5j.LowerCurv;
    PARALG5j.rho0 = max(PARALG5j.rho0place, ...
        PARALG5j.mu0+sqrt(PARALG5j.mu0^2+4));

    [time,NormResidualz,NormResidualR] = APFIAL(F,PARAPFj,x0,p0);
    TAPF(j,1) = time;
    PFeasAPF(j,1) = NormResidualz;
    DFeasAPF(j,1) = NormResidualR;

    % OIAL
    PAROIALj.epsilon = min( ...
        (PAROIALj.c*(PAROIALj.epsilonhat)^2)/(4*PAROIALj.D), ...
        PAROIALj.rhohat);
    [time,NormResidualz,NormResidualR] = OIAL(F,PAROIALj,x0,p0);
    TOIAL(j,1) = time;
    PFeasOIAL(j,1) = NormResidualz;
    DFeasOIAL(j,1) = NormResidualR;

    % Algorithm 5, specialized to Ax=b.
    [time,NormResidualz,NormResidualR] = Alg5PAL(F,PARALG5j,x0,p0);
    TALG5(j,1) = time;
    PFeasALG5(j,1) = NormResidualz;
    DFeasALG5(j,1) = NormResidualR;
end

writeLatexTable(filename, dimMvec, dimNvec, ...
    TOIAL, TAPF, TALG5, ...
    PFeasOIAL, PFeasAPF, PFeasALG5, ...
    DFeasOIAL, DFeasAPF, DFeasALG5);

system(['cat ', filename]);
end


function PARout = copyProblemParams(PARin,PARinitial)
PARout = PARin;
fields = {'constr_fn','grad_constr_fn','set_projector','K_constr', ...
    'prod_fn','norm_fn','M','m','x0','UpperCurv','LowerCurv', ...
    'CondNumber','w0','muEst','LipEst'};
for i = 1:length(fields)
    fld = fields{i};
    if isfield(PARinitial,fld)
        PARout.(fld) = PARinitial.(fld);
    end
end
end


function writeLatexTable(filename, dimMvec, dimNvec, ...
    TOIAL, TAPF, TALG5, ...
    PFeasOIAL, PFeasAPF, PFeasALG5, ...
    DFeasOIAL, DFeasAPF, DFeasALG5)

numofProblems = length(dimMvec);
fid = fopen(filename, 'w');
if fid < 0
    error('runAL:FileOpen','Could not open %s for writing.', filename);
end

fprintf(fid, '%s\n', '\begin{tabular}{|cc|ccc|ccc|ccc|} \hline');
fprintf(fid, '%s', [ ...
    '\multicolumn{2}{|c||}{Dimensions} & ', ...
    '\multicolumn{3}{c|}{Time (s)} & ', ...
    '\multicolumn{3}{c|}{Primal Feasibility Residual} & ', ...
    '\multicolumn{3}{c|}{Dual Stationarity Residual}']);
fprintf(fid, '%s\n', ...
    ' \\ \cline{1-2}\cline{3-5}\cline{6-8}\cline{9-11}');
fprintf(fid, '%s\n', [ ...
    '$m$ & $n$ & O-IAL & APF-IAL & ProxALM \cite{LuM} & ', ...
    'O-IAL & APF-IAL & ProxALM \cite{LuM} & ', ...
    'O-IAL & APF-IAL & ProxALM \cite{LuM} \\ \hline']);

for k = 1:numofProblems
    fprintf(fid, '%7.0f & %7.0f & ', dimMvec(k), dimNvec(k));
    fprintf(fid, '%7.2f & %7.2f & %7.2f & ', ...
        TOIAL(k), TAPF(k), TALG5(k));
    fprintf(fid, '%7.4e & %7.4e & %7.4e & ', ...
        PFeasOIAL(k), PFeasAPF(k), PFeasALG5(k));
    fprintf(fid, '%7.4e & %7.4e & %7.4e ', ...
        DFeasOIAL(k), DFeasAPF(k), DFeasALG5(k));
    fprintf(fid, '%s\n', '\\ \hline');
end
fprintf(fid, '%s\n', '\end{tabular}');
fclose(fid);
end
