TrainingSetAdaptation.m

classdef TrainingSetAdaptation < Greedy.MetaInterface
  % greedy basis generation extension which adaptively refines the trainings
  % parameter set.
  %
  % @todo add reference to adaptive paper.

  properties (Constant)

    % cell array of field names to be copied to the detailed data node instance
    % during init_basis() call
    info_fields = { 'max_refinement_steps', 'max_refinement_level', ...
                    'refinement_mode', 'refinement_theta', 'element_indicator_mode', ...
                    'element_indicator_s_max', 'force_stop_at_max_refinement_level' };
  end

  properties

    % maximum number of refinement steps before the algorithm stops.
    %
    % If this level is reached, the algorithm terminates and returns the
    % generated reduced basis spaces.
    max_refinement_steps    = inf;

    % maximum number of refinement levels in the refineable trainings parameter
    % set
    %
    % If this level is reached, the algorithm terminates and returns the
    % generated reduced basis spaces.
    max_refinement_level    = 5;

    % string specifying the method type for the adaptation of the parameter
    % sampling.
    %
    % Possible values are:
    %  - 'uniform': the parameter sampling is refined uniformly
    %  - 'adaptive': the parameter sampling is refined adaptively, such that
    %                only regions of the parameter space with high estimated
    %                errors are refined.
    refinement_mode         = 'uniform'; % 'uniform' or 'adaptive'

    % ratio of elements which are refined during an adaptive refinement step
    %
    % This value `\theta` is multiplied by the number of leaf elements in the
    % current parameter grid in order to compute the number grid cells to be
    % refined.
    refinement_theta        = 0.5;

    % string specifying the indicator for the training set nodes.
    %
    % Possible values:
    %   - 'nodes_cogs'
    %   - 'nodes'
    %   - 'nodes_skippedrefs'
    %   - 'nodes_cogs_skippedrefs'
    element_indicator_mode  = 'nodes_cogs';

    % specifies maximum number of skipped refinements in case of
    % '*_skippedrefs' indicator modes.
    element_indicator_s_max = 100;

    % boolean flag indicating whether at the maximum refinement level the
    % algorithm shall stop.
    %
    % If this is set to false, no further refinement is done at the last level,
    % but the reduced basis spaces are generated, anyways.
    force_stop_at_max_refinement_level = true;
  end

  methods
    function psa = TrainingSetAdaptation(child)
      % function psa = TrainingSetAdaptation(child)
      % constructor extending a Greedy algorithm object
      %
      % Parameters:
      %   child: the greedy algorithm of type Greedy.Interface to be extended
      %          by training set adaptation.
      assert(isa(child, 'Greedy.Interface'));
      psa = psa@Greedy.MetaInterface(child);
    end

    function detailed_data = init_basis(this, rmodel, model_data)
      % function detailed_data = init_basis(this, rmodel, model_data)
      % @copybrief Greedy.Interface.init_basis()
      %
      % @copydoc Greedy.Interface.init_basis()
      %
      % Parameters:
      %   rmodel: of type Greedy.User.IReducedModel
      detailed_data = init_basis(this.child, rmodel, model_data);

      if isstruct(model_data)
        detailed_data = Greedy.DataTree.Detailed.InfoNode(detailed_data, 'paramspace_adapt');
        set_fields(detailed_data, this, Greedy.TrainingSetAdaptation.info_fields);
      end
    end

    function detailed_data = basis_extension(this, rmodel, detailed_data, checkpoint)
      % function detailed_data_tree = basis_extension(this, rmodel, detailed_data_tree, checkpoint)
      % @copybrief GreedyInterfacebasis_extension()
      %
      % @copydoc GreedyInterfacebasis_extension()
      %
      % Parameters:
      %   rmodel: of type Greedy.User.IReducedModel
      %   detailed_data: of type Greedy.User.IDetailedData
      %
      % copied from basisgen_refined written by
      % \author Bernard Haasdonk

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

      basetoc = get(checkpoint, 'parameter_toc', 0);

      alltimer = tic;

      while true
        base_detailed_data = get(detailed_data, 1);
        base_detailed_data = basis_extension(this.child, rmodel, base_detailed_data, checkpoint.child(1));

        dd_leaf = get_active_leaf(base_detailed_data, rmodel);
        refinement_step = get_field(dd_leaf, 'refinement_step', 0);

        if stopped_on_any_child(base_detailed_data, 'stopped_on_timeout')
          set_field(detailed_data, 'elapsed_time', toc(alltimer) + basetoc);
          return;
        end

        if ~stopped_on_active_child(base_detailed_data, 'stopped_on_max_val_train_ratio', rmodel)
          detailed_data.stop_flags = union(detailed_data.stop_flags, base_detailed_data.stop_flags);
          break;
        end

        if refinement_step >= get_field(detailed_data, 'max_refinement_steps')
          set_stop_flag(dd_leaf, 'stopped_on_max_refinement_level');
          break;
        end

        M_train = get_field(dd_leaf, 'M_train');
        if isempty(M_train) || ~isa(M_train, 'ParameterSampling.IRefineable')
          break;
        end

        levels = get_refinement_levels(M_train);
        maxlev = max(levels);
        if maxlev == get_field(detailed_data, 'max_refinement_level') - 1 ...
            && ~get_field(detailed_data, 'force_stop_at_max_refinement_level')

          set_field(dd_leaf, 'stop_max_val_train_ratio', 1e6);

        end

        % if refinement is still possible and required val-train ratio is not assured
        % refine the parameter space regardless of stopped_on_epsilon flag.
        r_values = get_field(dd_leaf, 'r_value_sequence');
        if maxlev < get_field(detailed_data, 'max_refinement_level') ...
            && r_values(end) > get_field(dd_leaf, 'stop_max_val_train_ratio')

          set_stop_flag(dd_leaf, 'stopped_on_epsilon', false);
        end

        if maxlev >= get_field(detailed_data, 'max_refinement_level')
          set_stop_flag(dd_leaf, 'stopped_on_max_refinement_level');
          break;
        end

        disp('detected grid refinement necessary!!');
        disp(['max val train ratio was: ', num2str(r_values(end))]);
        nleafelements = size(M_train);

        snapshot(dd_leaf, 'ParamSpaceAdapt', 'Parameter space adaptation', ...
          ['Refinement step no. ', num2str(refinement_step)]);

        switch get_field(detailed_data, 'refinement_mode')
          case 'uniform'
            set_field(dd_leaf, 'M_train', refine(M_train));
          case 'adaptive'

            eta = refinement_element_indicators( this, rmodel, detailed_data, ...
                                                 dd_leaf );

            [dummy, ind] = sort(-eta);
            nmax = ceil(nleafelements * get_field(detailed_data, 'refinement_theta'));

            set_field(dd_leaf, 'M_train',  refine(M_train, ind(1:nmax)));
        end

        set_stop_flag(dd_leaf, 'stopped_on_max_val_train_ratio', false);

        refinement_step = refinement_step + 1;
        set_field(dd_leaf, 'refinement_step', refinement_step);
        checkpoint = checkpoint.store(rmodel, detailed_data, 'refinement', ...
                                       struct('parameter_toc', toc(alltimer) + basetoc));
      end
      set_field(detailed_data, 'elapsed_time', toc(alltimer) + basetoc);
    end

    function eta = refinement_element_indicators( this, rmodel, detailed_data, dd_leaf )
      % function eta = refinement_element_indicators( this, rmodel, detailed_data, dd_leaf )
      % evaluates the refinement indicators
      %
      % Parameters:
      %   rmodel: of type Greedy.User.IReducedModel
      %   detailed_data: of type Greedy.User.IDetailedData
      %   dd_leaf: Detailed data node of type Greedy.DataTree.Detailed.ILeafNode
      %            which is currently extended leaf in the detailed data tree.
      Delta_train        = get_field(dd_leaf, 'M_last_errors');
      MMesh              = get_field(dd_leaf, 'M_train');
      base_detailed_data = get(detailed_data, 1);

      if isempty(Delta_train)
        Delta_train = error_indicators(this.child, rmodel, base_detailed_data, MMesh.sample);
      end

      if isequal(get_field(detailed_data, 'element_indicator_mode'), 'nodes_cogs') ...
          || isequal(get_field(detailed_data, 'element_indicator_mode'), 'nodes_cogs_skippedrefs')

        % determine all cog vertex estimators
        M_cog     = cogs(MMesh);
        Delta_cog = error_indicators(this.child, rmodel, base_detailed_data, M_cog, true);
        set_field(dd_leaf, 'Delta_cog', Delta_cog);
      end

      Max_Delta_vertex = elementwise_maximum_of_vertex_values(MMesh, Delta_train);

      switch get_field(detailed_data, 'element_indicator_mode')
        case 'nodes_cogs'
          %          eta(i) =  max_{mu in (V(e) and cog(e))}  Delta (Phi_t, mu)
          eta = max( Delta_cog(:), Max_Delta_vertex(:) );
        case 'nodes'
          %         eta(i) =  max_{mu in (V(e))}   Delta (Phi_t, mu)
          eta = Max_Delta_vertex(:);
        case 'nodes_skippedrefs'
          %          eta(i) =  max_{mu in (V(e))}   Delta (Phi_t, mu) +
          %                                     s(i)/s_max * epsilon
          % determine number of skipped refinements of all leaf elements
          s = skipped_refinements(MMesh);
          eta = Max_Delta_vertex(:) + ...
            max(Max_Delta_vertex) / ...
            get_field(detailed_data, 'element_indicator_s_max') * ...
            s(:);
        case 'nodes_cogs_skippedrefs'
          %          eta(i) =  max_{mu in (V(e) and cogs)}   Delta (Phi_t, mu) +
          %                                     s(i)/s_max * epsilon
          % determine number of skipped refinements of all leaf elements
          s = skipped_refinements(MMesh);
          eta = max(Delta_cog(:),Max_Delta_vertex(:)) + ...
            max(Max_Delta_vertex) / ...
            get_field(detailed_data, 'element_indicator_s_max') * ...
            s(:);
        otherwise
          error('element_indicator_mode in grid-refinement unknown');
      end;
    end

  end

end