% SPDX-License-Identifier: MIT
% Copyright © 2021 Weiwei "William" Kong

function [summary_tables, comp_models] = run_CCM_benchmark(base_cc_model, fw_arr, solver_arr, fw_hparams_arr, name_arr)
% A wrapper that runs several solvers on a user-provided ConstrCompModel object.
%
% Arguments:
%  
%   base_cc_model (ConstrCompModel): A ConstrCompModel object for running the frameworks on
%
%   fw_arr ({function}): An array of frameworks to benchmark
%
%   solver_arr ({function}): An array of solvers that correspond to the frameworks
%
%   fw_hparams_arr ({struct}): An array of framework hyperparameters (can be left empty)
% 
%   name_arr ({string}): An array of framework names (can be left empty)
% 
% Returns:
%
%   Two structs. The first contains some summary tables related to the execution of the method. The second contains some elements
%   related to the solutions generated by the algorithms
% 

  % Check lengths.
  if (length(fw_arr) ~= length(solver_arr))
    error('The number of frameworks must equal the number of solvers!');
  end
  if (~isempty(fw_hparams_arr) && length(fw_arr) ~= length(fw_hparams_arr))
    error('The number of frameworks must equal the number of hyperparameters when the latter is nonempty!');
  end
  if (~isempty(name_arr) && length(solver_arr) ~= length(name_arr))
    error('The number of solvers must equal the number of names when the latter is nonempty!');
  end

  % Initialize.
  n_frameworks = length(fw_arr);
  summary_tables.runtime = table();
  summary_tables.iter = table();
  summary_tables.fval = table();
  summary_tables.pdata = table(); % Problem data.
  summary_tables.mdata = table(); % Model data.
  
  % Run an optimization routine for each algorithm.
  for i=1:n_frameworks
    % Prepare.
    constr_comp_model = copy(base_cc_model);
    framework = fw_arr{i};
    solver = solver_arr{i};
    if ~isempty(name_arr)
      framework_name = name_arr{i};
    else
      framework_name = [func2str(framework), '_', func2str(solver)];
    end
    constr_comp_model.solver = solver;
    constr_comp_model.framework = framework;
    if (~isempty(fw_hparams_arr))
      constr_comp_model.solver_hparams = fw_hparams_arr{i};
    end
    % Optimize.
    constr_comp_model.optimize;
    % Record the results.
    comp_models.(framework_name) = constr_comp_model;
    summary_tables.runtime = add_column(['t_', framework_name], constr_comp_model.runtime, summary_tables.runtime);
    summary_tables.iter = add_column(['iter_', framework_name], constr_comp_model.iter, summary_tables.iter);
    summary_tables.fval = add_column(['fval_', framework_name], constr_comp_model.f_at_x, summary_tables.fval);
    % Auxiliary results.
    if (isfield(constr_comp_model.history, 'c0'))
      summary_tables.mdata = add_column(['c0_', framework_name], constr_comp_model.history.c0, summary_tables.mdata);
    end
    if (isfield(constr_comp_model.history, 'c'))
      summary_tables.mdata = add_column(['c_', framework_name], constr_comp_model.history.c, summary_tables.mdata);
    end
    if (isfield(constr_comp_model.history, 'wavg_c'))
      summary_tables.mdata = add_column(['wavg_c_', framework_name], constr_comp_model.history.wavg_c, summary_tables.mdata);
    end
    if (isfield(constr_comp_model.history, 'stage'))
      summary_tables.mdata = add_column(['stage_', framework_name], constr_comp_model.history.stage, summary_tables.mdata);
    end
    if (isfield(constr_comp_model.history, 'outer_iter'))
      summary_tables.mdata = add_column(['oiter_', framework_name], constr_comp_model.history.outer_iter, summary_tables.mdata);
    end
    if (isfield(constr_comp_model.history, 'k0_avg'))
      summary_tables.mdata = add_column(['k0_avg_', framework_name], constr_comp_model.history.k0_avg, summary_tables.mdata);
    end
    if (isfield(constr_comp_model.history, 'k0_max'))
      summary_tables.mdata = add_column(['k0_max_', framework_name], constr_comp_model.history.k0_avg, summary_tables.mdata);
    end
  end
  
  % Add solver independent data.
  summary_tables.pdata = add_column('m', constr_comp_model.m, summary_tables.pdata);
  summary_tables.pdata = add_column('M', constr_comp_model.M, summary_tables.pdata);
  summary_tables.pdata = add_column('K_constr', constr_comp_model.K_constr, summary_tables.pdata);
  summary_tables.pdata = add_column('L_constr', constr_comp_model.L_constr, summary_tables.pdata);
  summary_tables.mdata = add_column('feas_tol', constr_comp_model.feas_tol, summary_tables.mdata);
  summary_tables.mdata = add_column('feas_type', convertCharsToStrings(constr_comp_model.feas_type), summary_tables.mdata);
  summary_tables.mdata = add_column('opt_tol', constr_comp_model.opt_tol, summary_tables.mdata);
  summary_tables.mdata = add_column('opt_type', convertCharsToStrings(constr_comp_model.opt_type), summary_tables.mdata);

  % Merge summary tables for easy reading
  summary_tables.all = ...
    [summary_tables.pdata, summary_tables.fval, summary_tables.iter, summary_tables.runtime, summary_tables.mdata];

end

% Utility function to add a column to a table if the value is nonempty
function out_tbl = add_column(col_name, col_val, in_tbl)
  if ~isempty(col_val)
    tmp_tbl = table(col_val);
    tmp_tbl.Properties.VariableNames = {col_name};
    out_tbl = [in_tbl, tmp_tbl];
  end
end