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

function [summary_tables, comp_models] = run_CM_benchmark(base_c_model, solver_arr, solver_hparams_arr, name_arr)
% A wrapper that runs several solvers on a user-provided CompModel object.
%
% Arguments:
%  
%   base_cc_model (CompModel): A CompModel object for running the solvers on
%
%   solver_arr ({function}): An array of solvers that correspond to the frameworks
%
%   solver_hparams_arr ({struct}): An array of solver 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
% 

  % Initialize
  if (~isempty(solver_hparams_arr) && length(solver_arr) ~= length(solver_hparams_arr))
    error('The number of solvers must equal the number of hyperparameters arrays when the latter is nonempty!');
  end
  n_solvers = length(solver_arr);
  summary_tables.runtime = table();
  summary_tables.iter = table();
  summary_tables.fn_iter = table();
  summary_tables.grad_iter = table();
  summary_tables.fval = table();
  summary_tables.normV = table();
  summary_tables.pdata = table(); % Problem data.
  summary_tables.mdata = table(); % Model data.
  
  % Run an optimization routine for each algorithm.
  for i=1:n_solvers
    % Prepare.
    if iscell(base_c_model)
      comp_model = base_c_model{i};
    else
      comp_model = copy(base_c_model);
    end
    solver = solver_arr{i};
    if ~isempty(name_arr)
      solver_name = name_arr{i};
    else
      solver_name = func2str(solver);
    end
    comp_model.solver = solver;
    if (~isempty(solver_hparams_arr))
      comp_model.solver_hparams = solver_hparams_arr{i};
    end
    % Optimize.
    comp_model.optimize;
    % Record the results
    comp_models.(solver_name) = comp_model;
    summary_tables.runtime = add_column(['t_', solver_name], comp_model.runtime, summary_tables.runtime);
    summary_tables.iter = add_column(['iter_', solver_name], comp_model.iter, summary_tables.iter);
    summary_tables.fn_iter = add_column(['fn_iter_', solver_name], comp_model.fn_iter, summary_tables.fn_iter);
    summary_tables.grad_iter = add_column(['grad_iter_', solver_name], comp_model.grad_iter, summary_tables.grad_iter);
    summary_tables.fval = add_column(['fval_', solver_name], comp_model.f_at_x, summary_tables.fval);
  end
  
  % Add solver independent data.
  summary_tables.pdata = add_column('m', comp_model.m, summary_tables.pdata);
  summary_tables.pdata = add_column('M', comp_model.M, summary_tables.pdata);
  summary_tables.mdata = add_column('opt_tol', comp_model.opt_tol, summary_tables.mdata);
  summary_tables.mdata = add_column('opt_type', convertCharsToStrings(comp_model.opt_type), summary_tables.mdata);

  % Merge summary tables for easy reading
  summary_tables.all = ...
    [summary_tables.pdata, summary_tables.fval, summary_tables.normV, summary_tables.iter, summary_tables.fn_iter, ...
     summary_tables.grad_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