function convergence_output = convergence_study(convergence_function,structure_function,iterative_params,fix_params) 
%% convergence_output = convergence_study(convergence_function,structure_function,iterative_params,fix_params) 
% Make the convergence study of the specified functions over an iterative
% parameter (number of grids points, number of time points...). The
% convergence functions must take as inputs the structures [L,LX]. These 
% structures are given by an initialisation function (could be fbt, rzp,
% fge...) that must be iterated over the convergence loop. [L,LX] change at
% every loop according to the change in the input parameters of the
% initialisation function. The inputs of the initialisation function could 
% be partially fixed (never change) and partially variable over the 
% convergence loop (for example, one can vary the time step dt or the grid 
% dimensions nr, nz). The output of the convergence functions are then 
% collected in a structure. The element of the structure are named as the
% convergence functions. If one of the functions is an anonymous handle 
% function, its name will f#, where # is the index were the function is
% placed in the cell array.
%
% Input
%   convergence_function: function handle or cell array of function
%     handles. These are the function over which the convergence study is
%     performed. The input of these functions must be [L,LX] obtained by
%     structure_function. The output must be an array.
%   structure_function: function handle to initialise the common variables 
%     for all the functions contained in convergence_function (if more than
%     one is used). Its input is a combination of iterative_params and
%     fix_params. Its output are [L,LX], to be given to
%     convergence_function.
%   iterative_params: cell array of parameters over which the convergence 
%     study is run. Could change more than one parameteter each time, but
%     the dimension of the different parameters must be the same (nstep). 
%   fix_params: cell array of input parameters for the structure_function
%     which don't change during the convergence loop.
%
% Output
%   convergence_output: struct with the outputs of each function. The
%     elements of the structure are arrays were the first dimension refers  
%     to the number of convergence steps (nsteps), while the other
%     dimensions are the same as the outputs of the convergence_function.
%
% An example of its use is the following:
%   
%   %Define convergence functions 
%   f1 = @(L,LX) sum(LX.Iy(:)); f2 = @(L,LX) sum(L.rrx(:));
%   convergence_function = {f1,f2};
% 
%   %Define structure_function
%   structure_function = @rzp;
% 
%   %Define iterative params (with the same number of elements):
%   nstep = 3;
%   nz = 2.^(log2(8)+(1:nstep)-1);
%   nr = round(nz.*28./32);
%   iterative_params = {repmat({'nr'},1,nstep), nr, repmat({'nz'},1,nstep), nz};
% 
%   %Define fix_paramas
%   fix_params = {'tcv',61400,0.5,'selx',''};
%   
%   %Call the convergence study for all the functions defined
%   C = convergence_study(convergence_function,structure_function,iterative_params,fix_params);
% 
%   %Plot the output of the convergence study
%   plot_convergence_study(1./nr,C.f1); plot_convergence_study(1./nr,C.f2);
%
%
% See also plot_convergence_study
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

assert(~any(diff(cellfun(@numel,iterative_params))),'Cell array iterative_params must have the same number of elements inside all cell');

nstep = numel(iterative_params{1});
nfunc = numel(convergence_function);

%% Concatentate the two cell arrays for fixed and iterative paramaters
structure_function_args = structure_convert(fix_params,iterative_params);
func_string = cell(1,nfunc);

%% Convergence loop
for ii = 1:nstep
    %Produce iterated [L,LX] structure with structure_function
    [L,LX] = structure_function(structure_function_args{ii}{:});
    
    %Produce the convergence output from the convergence_function
    for jj = 1:nfunc
      %Anonymous function handling 
      func_string{jj} = func2str(convergence_function{jj});
      if strcmp(func_string{jj}(1),'@')
        if strcmp(func_string{jj}(1:7),'@(L,LX)')
          func_string{jj} = sprintf('f%d',jj);
        else
          error('Anonymous handle functions must get as inputs @(L,LX)');
        end
      end
      Citer.(func_string{jj}){ii} = convergence_function{jj}(L,LX); 
    end
end

%% Reconvert the output such as the first dimension refers to the number of convergence steps (nsteps)
for jj = 1:nfunc
  nn = ndims(Citer.(func_string{jj}){1})+1;
  convergence_output.(func_string{jj}) = permute(cat(nn,Citer.(func_string{jj}){:}),[nn,1:nn-1]);
end
end

%% Cat the cell arrays for the fixed and iterative parameters
function S = structure_convert(fix_cell,iterative_cell)
  S = cell(1,numel(iterative_cell{1}));
  for ii = 1:numel(iterative_cell{1})
    S{ii} = fix_cell;
    for jj = 1:length(iterative_cell)
      S{ii} = [S{ii} iterative_cell{jj}(ii)];
    end
  end
end