Algorithm.m

classdef Algorithm < Greedy.Interface
  % default greedy basis generation class
  %
  % This implements an abstract version of the default greedy algorithm for
  % reduced basis and empirical interpolation basis generation. The actual
  % implementation of the extension and checking and error evaluation routines
  % is delegated to an underlying object implementing a Greedy.Plugin.Interface. See
  % \ref basisgen for example usage information.
  %
  % The entry point for the greedy algorithm is the method basis_extension().
  % See there for more details on the algorithm

  properties(Constant)
    % cell array of fields to be copied to the detailed data leaf node instance
    % during init_basis().
    info_fields = { 'M_train', 'M_validation', ...
                    'stop_epsilon', 'stop_max_val_train_ratio', ...
                    'val_train_ratio_seq_for_break', ...
                    'stop_timeout' };
  end

  properties(Dependent)
    % an object of type SnapshotsGenerator.Cached
    %
    % This object is inherited from the #detailed_extension.
    detailed_gen;
  end

  properties(SetAccess=private, Dependent)
    % @copybrief Greedy.Interface.generated_basis_type
    %
    % @copydoc Greedy.Interface.generated_basis_type
    generated_basis_type;

    % @copybrief Greedy.Interface.generated_basis_type
    %
    % @copydoc Greedy.Interface.id
    id;
  end

  properties(SetAccess = private)
    % object of type Greedy.Plugin.Interface
    detailed_extension;
  end

  properties

    % double specifying the maximum error indicator for which the basis
    % generation shall be stopped.
    stop_epsilon             = 1e-2;

    % integer specifying the number of seconds after which the basis generation
    % shall be stopped.
    stop_timeout             = inf;

    % positive double value specifying the maximum ratio between the maximum
    % error indicator over the validation paramter set and the maximum error
    % over the trainining parameter set, for which the basis generation is
    % stopped.
    %
    % Note, that #M_validation needs to be non-empty, if this value is set to
    % something different than 'inf'.
    stop_max_val_train_ratio = inf;

    % number of subsequent cases where the max_err_val_train_ratio is too high
    % that is needed in order to break with a
    % Greedy.Info.Base.stop_max_val_train_ratio condition.
    %
    % The default is '1', but especially for basis generation with
    % Greedy.Plugin.PODEI extension algorithm, it is a good idea to set it to
    % at least '2', because this algorithms can produce very strange basis
    % spaces in between.
    %
    val_train_ratio_seq_for_break = 1;

    % object of type ParameterSampling.Interface defining the training parameter set
    % `M_{\text{train}} \subset \cal M`
    M_train                  = ParameterSampling.Uniform(1);

    % object of type ParameterSampling.Interface defining the validation parameter set
    % `M_{\text{val}} \subset \cal M`
    M_validation             = [];
  end

  methods
    function bgd = Algorithm(detailed_extension, M_train, initial_val_param_set)
    % constructor for a greedy basis generation algorithm
    %
    % This constructs a greedy basis generation algorithm using the following
    %
    % Parameters:
    %   detailed_extension: the underlying object of type Greedy.Plugin.Interface.
    %   M_train:   an object of type ParameterSampling.Interface specifying an initial
    %              training parameter set.
    %   initial_val_param_set: an object of type ParameterSampling.Interface specifying an initial
    %              validation parameter set.
      if nargin >=2
        bgd.M_train = M_train;
      end
      if nargin >= 3
        bgd.M_validation = initial_val_param_set;
      end
      bgd.detailed_extension = detailed_extension;
    end

    function enabled = enable_validation(this, detailed_data)
    % function enabled = enable_validation(this, detailed_data)
    % indicates whether the validation routines shall be executed
    %
    % This method returns true when a validation parameter set is given and the
    % property #stop_max_val_train_ratio is not set to 'inf'
    %
    % Parameters:
    %   detailed_data: of type .Greedy.DataTree.Detailed.ILeafNode
    %
    % Return values:
    %   enabled: boolean value
        enabled = ~isempty(detailed_data) && ~isempty(get_field(detailed_data, 'M_validation')) ...
                  && ~isinf(get_field(detailed_data, 'stop_max_val_train_ratio'));
    end

    function btype = get.generated_basis_type(this)
      btype = this.detailed_extension.generated_basis_type;
    end

    function detailed_gen = get.detailed_gen(this)
      detailed_gen = this.detailed_extension.generator;
    end

    function id = get.id(this)
      id = this.detailed_extension.id;
    end

    function detailed_data_tree = basis_extension(this, rmodel, detailed_data_tree, checkpoint)
    % function detailed_data_tree = basis_extension(this, rmodel, detailed_data_tree, checkpoint)
    % greedy basis generation extension
    %
    % @pre If required by the @ref Greedy.Plugin.Interface.error_indicators()
    % "error_indicators()" routine, the method prepare() should have been
    % called beforehand
    % @pre The argument 'detailed_data_tree' should be initialized via by the
    % init_basis() method.
    %
    % During each extension step, the methods from the underlying #detailed_extension
    % are called in this order: 
    %  -# The maximum error indicators are gather by calling
    %     @ref Greedy.Plugin.Interface.error_indicators() "error_indicators()"
    %     for the training set parameters, and if existent also vor the
    %     validation set parameters.
    %  -# Before the actual basis extension, we check whether any termination
    %     condition has been reached by computing some default checks
    %     documented below and extension specific checks in
    %     @ref Greedy.Plugin.Interface.pre_check_for_end() "pre_check_for_end()".
    %  -# For the worst approximated parameter a basis extension call is
    %     initiated by calling the method
    %     @ref Greedy.Plugin.Interface.basis_extension() "basis_extension()".
    %
    % Default checks:
    %   By default this methods checks
    %     - whether maximum value for the error indicator dropped below the
    %       value of 'get_field(detailed_data, ''''stop_epsilon'''')',
    %     - whether the timeout barrier #stop_timeout has been reached, or
    %     - whether the ratio between the maximum error indicators over the
    %       training set and the validation set exceeds the barrier of
    %       'get_field(detailed_data, ''''stop_max_val_train_ratio'''')'
    %
    % Parameters:
    %  detailed_data_tree: either a struct containing high dimensional model data
    %                      needed to execute detailed simulations or an object
    %                      of type Greedy.DataTree.Detailed.INode in case of resume from a
    %                      checkpoint.
    %  checkpoint: object of type Greedy.Checkpoint specifying a given point in the
    %              algorithm where it basis generation can be resumed.
    %
    % Return values:
    %  detailed_data_tree: object of type Greedy.DataTree.Detailed.INode storing the
    %                      reduced basis information in the leaf nodes and
    %                      information on the reduced basis generation in every
    %                      node.

      if nargin == 3 || isempty(checkpoint)
        checkpoint = Greedy.Checkpoint;
      end

      basetoc = get(checkpoint, 'toc', 0);
      % prepare(this, rmodel, detailed_data_tree);

      alltimer = tic;

      dd_leaf = get_active_leaf(detailed_data_tree, rmodel);
      cur_M_train      = get_field(dd_leaf, 'M_train');
      cur_M_validation = get_field(dd_leaf, 'M_validation');

      % ATTENTION: We assume, that prepare has been called beforehand,
      % if necessary.
      % prepare(this, rmodel, detailed_data_tree);

      while true

        steptic = tic;

        rmodel.M = get_ei_size(detailed_data_tree);
        rmodel.N = get_rb_size(detailed_data_tree);

        fprintf('\n========================================================\n');
        disp(['Extended reduced bases to size N = ',...
              num2str(rmodel.N),...
              ', M = ', num2str(rmodel.M), '.']);

        disp(['Computing ', num2str(size(cur_M_train.sample,1)),...
          ' error indicators for basis extension:']);
        [errs, max_err_seq, muind] = error_indicators(this,...
                                         rmodel,...
                                         detailed_data_tree,...
                                         cur_M_train.sample);
        fprintf('\n');
        this.pretty_print_errs('Error sequence: ', errs);

        if enable_validation(this, dd_leaf)

          disp(['Computing ', num2str(size(cur_M_validation.sample,1)),...
            ' error indicators for validating basis extension:']);
          val_errs    = error_indicators(this,...
                                         rmodel,...
                                         detailed_data_tree,...
                                         cur_M_validation.sample,...
                                         true);

          r_value = max(val_errs) / max(errs);
          disp(['Max validation/training ratio: ', num2str(r_value)]);
          append_field(dd_leaf, 'r_value_sequence', r_value);
        end

        Greedy.Algorithm.push_back_extension_info(dd_leaf, errs, muind, max_err_seq);
        if pre_check_for_end_meta(this, rmodel, detailed_data_tree, toc(alltimer)+basetoc);
          set_field(dd_leaf, 'M_last_errors', errs);
          break;
        end

        mu = cur_M_train.sample(muind,:);
        disp(['Choosing mu = ', num2str(mu(:)'), ...
              ' for basis extension (error = ', num2str(max(errs)), ')']);

        detailed_data_tree = basis_extension( this.detailed_extension,...
                                              rmodel,...
                                              detailed_data_tree,...
                                              max_err_seq, mu );

        % re-read dd_leaf, because basis extension algorithm is allowed to change the pointer
        dd_leaf = get_active_leaf(detailed_data_tree, rmodel);

        rmodel.N = get_rb_size(detailed_data_tree);
        rmodel.M = get_ei_size(detailed_data_tree);

        %! \todo check whether we should add another line
        % errs = error_indicators(this, rmodel, this.detailed_gen, detailed_data);
        %  and do the validation afterwards!

        append_field(dd_leaf, 'toc_value_sequence', toc(steptic));

        checkpoint = checkpoint.store(rmodel, detailed_data_tree,...
                                       'greedy_extension', ...
                                       struct('toc', toc(alltimer)+basetoc));

      end

      detailed_data_tree = finalize(this.detailed_extension, rmodel, detailed_data_tree);
    end

    function [max_errs, max_err_sequence, max_mu_index] = error_indicators(this, rmodel, detailed_data_tree, parameter_set, reuse_reduced_data)
    % function [max_errs, max_err_sequence, max_mu_index] = error_indicators(this, rmodel, detailed_data_tree, parameter_set, reuse_reduced_data)
    % routine computing the error indicators
    %
    % Parameters:
    %   detailed_data_tree: an object of type GreedyDataTreeDetailed.INode storing the
    %                       reduced basis functions in the leaf nodes.
    %   parameter_set:      A matrix of dimension
    %                       `\text{npar} \times \dim(\cal M)` containing in
    %                       each row a parameter vector for which the error
    %                       indicator should be computed.
    %   reuse_reduced_data: boolean flat indicating whether it is necessary to
    %                       recompute the reduced data or whether it should
    %                       still be valid since the last call to
    %                       error_indicators().
    %
    % Return values:
    %   max_errs: For transient problems this returns the error sequence over
    %             time for the parameter 'max_mu_index' with the worst error.
    %   max_err_sequence: a vector of the same length as the number of training
    %                     set parameters, contining the maximum error in the
    %                     error sequence over time for each of those
    %                     parameters.
    %   max_mu_index:     the parameter index of the parameter `\mu_{\max}` in
    %                     the training set.
      if nargin < 5
        reuse_reduced_data = false;
      end
      if nargin < 4
        parameter_set = [];
      end
      [max_errs, max_err_sequence, max_mu_index] = error_indicators(this.detailed_extension, rmodel, detailed_data_tree, parameter_set, reuse_reduced_data);
    end

    % initialization routine for basis extension
    %
    % This method is run by the gen_detailed_data() method before the execution
    % of the init_basis() methods and should is used for preparation purposes of
    %  - the training and validation parameter sample and
    %  - caching detailed simulations if necessary for error_indicators()
    %    computations.
    function prepare(this, rmodel, model_data)
      dmodel = rmodel.detailed_model;
      if init_required(this.M_train)
        init_sample(this.M_train, rmodel);
      end

      if isa(model_data, 'Greedy.DataTree.Detailed.INode')
        dd_leaf = get_active_leaf(model_data, rmodel);
      else
        dd_leaf = [];
      end

      validation_enabled = false;
      if enable_validation(this, dd_leaf)
        Msamples = get_field(dd_leaf, 'M_validation');
        validation_enabled = true;
      elseif isempty(dd_leaf) && ~isempty(this.M_validation) && ~isinf(this.stop_max_val_train_ratio)
        Msamples = this.M_validation;
        validation_enabled = true;
      end
      if validation_enabled && init_required(Msamples)
        init_sample(Msamples, rmodel);
      end

      if this.detailed_extension.needs_preparation
        if isa(model_data, 'Greedy.DataTree.Detailed.INode')
          detailed_data = model_data;
          model_data    = detailed_data.model_data;
        end

        prepare(this.detailed_gen, dmodel, model_data, this.M_train.sample);

        % prepare validation set if existent
        if validation_enabled
          prepare(this.detailed_gen, dmodel, model_data, Msamples.sample);
        end
      end
    end

    function detailed_data_tree = init_basis(this, rmodel, model_data)
      % function detailed_data_tree = init_basis(this, model_data)
      % @copybrief GreedyInterface.init_basis()
      %
      % @copydoc GreedyInterface.init_basis()
      %
      % Parameters:
      %   rmodel: of type Greedy.User.IReducedModel
      if isa(model_data, 'Greedy.DataTree.Detailed.INode')
        this.M_train = get_field_on_active_child(model_data, 'M_train', rmodel);
      end
      M = this.M_train.sample;
      detailed_data_tree = init_basis(this.detailed_extension, rmodel, model_data, M);
      leaf_descrs = get_active_leaf_description(detailed_data_tree, rmodel);
      for leaf_descr = leaf_descrs
        detailed_data_leaf = get(detailed_data_tree, leaf_descr.basepath);

        if isempty(get_field(detailed_data_leaf, Greedy.Algorithm.info_fields{1}, []))
          set_fields(detailed_data_leaf, this, Greedy.Algorithm.info_fields);
        end
      end
    end

  end

  methods (Access=private)
    function breakcondition = pre_check_for_end_meta(this, rmodel, detailed_data_tree, time)

      dd_leaf = get_active_leaf(detailed_data_tree, rmodel);

      breakcondition = false;
      max_err_sequence = get_field(dd_leaf, 'max_err_sequence');
      if max_err_sequence(end) < get_field(dd_leaf, 'stop_epsilon')
        set_stop_flag(dd_leaf, 'stopped_on_epsilon');
        breakcondition          = true;
      elseif time > this.stop_timeout
        set_stop_flag(dd_leaf, 'stopped_on_timeout');
        breakcondition          = true;
      end
      if enable_validation(this, dd_leaf)
        val_seq_for_break = get_field(dd_leaf, 'val_train_ratio_seq_for_break');
        r_values          = get_field(dd_leaf, 'r_value_sequence');
        minindex = min(val_seq_for_break, length(r_values))-1;
        if all(r_values(end-minindex:end) > get_field(dd_leaf, 'stop_max_val_train_ratio'))
          set_stop_flag(dd_leaf, 'stopped_on_max_val_train_ratio');
          breakcondition                          = true;
        end
      end
      if ~breakcondition
        breakcondition = pre_check_for_end(this.detailed_extension, rmodel, detailed_data_tree);
      end
      set_field(detailed_data_tree, 'elapsed_time', time);

    end
  end

  methods (Access = private, Static)
    function pretty_print_errs(title, errs)
      nerrs = length(errs);
      rows  = floor(nerrs/6);
      resherrs = reshape(errs(1:rows*6), rows, 6);
      resterrs = errs(rows*6+1:end);
      format   = '%1.4e  ';
      disp(char(title, num2str(resherrs, format), num2str(resterrs(:)', format)));
    end

    function push_back_extension_info(detailed_data, errs, muind, max_err_seq)

      append_field(detailed_data, 'max_err_sequence', max(errs));
      append_field(detailed_data, 'errs_sequence', {errs});
      append_field(detailed_data, 'mu_ind_sequence', muind);
      cur_M_train = get_field(detailed_data, 'M_train');
      append_field(detailed_data, 'mu_sequence', cur_M_train.sample(muind,:)');
      set_field(detailed_data, 'M_last_errors', errs);

      [dummy, maxt]   = max(max_err_seq);
      append_field(detailed_data, 'max_time_index_sequence', maxt);
    end
  end


end