rbasis/basisgen/+Greedy/+Plugin/Default.m

classdef Default < Greedy.Plugin.Interface
  % Default implementation of a Greedy.Plugin.Interface interface class

  properties (Transient)
    % temporary handle to the last object of type Greedy.User.ReducedData
    % computed by prepare_reduced_data().
    reduced_data = [];
  end

  properties (Dependent)

    % @copybrief Greedy.Plugin.Interface.needs_preparation
    %
    % @copydetails Greedy.Plugin.Interface.needs_preparation
    %
    needs_preparation;
  end

  properties

    % @copybrief .Greedy.Plugin.Interface.indicator_mode
    %
    % @copydetails .Greedy.Plugin.Interface.indicator_mode
    indicator_mode = 'error';

    % boolean flag indicating whether the `L^2(\Omega)`-norm is used by
    % compute_error().
    %
    % Otherwise the `L^{\infty}(\Omega)`-norm is applied.
    use_l2_error   = true;

    % @copybrief Greedy.Plugin.Interface.relative_error
    %
    % @copydetails Greedy.Plugin.Interface.relative_error
    relative_error = false;
  end

  methods

    function ged = Default(generator)
      % function ged = Default(generator)
      % constructor for a greedy extension object
      %
      % Parameters:
      %   generator: object of type SnapshotsGenerator.Cached generating the
      %              (high dimensional) basis functions
      ged = ged@Greedy.Plugin.Interface(generator);
    end

    function need_prepare = get.needs_preparation(this)
      need_prepare = isequal(this.indicator_mode, 'error');
    end

    function [max_errs, max_err_sequence, max_mu_index] = ...
        error_indicators(this, rmodel, detailed_data, parameter_set, reuse_reduced_data)
      % function [max_errs, max_err_sequence, max_mu_index] = error_indicators(this, rmodel, detailed_data, parameter_set, reuse_reduced_data);
      % @copybrief Greedy.Plugin.Interface.error_indicators()
      %
      % @copydetails Greedy.Plugin.Interface.error_indicators()
      %
      % Parameters:
      %   rmodel: an object of type .Greedy.User.IReducedModel
      %   parameter_set : an object of type ParameterSampling.Interface
      %   detailed_data: of type .Greedy.User.IDetailedData
      %
      % @note This method is controlled by the indicator_mode variable: 
      %
      % If 'indicator_mode' equals
      %  - 'estimator', the method error_estimators() is called
      %  - 'error', compute_error() is used to compute "true" errors.

      if nargin < 5 || isempty(reuse_reduced_data)
        reuse_reduced_data = false;
      end

      if ~reuse_reduced_data
        prepare_reduced_data(this, rmodel, detailed_data);
      end

      if isequal(this.indicator_mode, 'error')
        nmus   = size(parameter_set, 1);

        tmperrs   = cell(1, nmus);
        tmpmaxerr = zeros(1, nmus);

        rd = this.reduced_data;
        %parfor i = 1:nmus
        for i = 1:nmus
          fprintf('.');
          trmodel = rmodel;
          tthis   = this;
          trd     = rd;
          set_mu(trmodel, parameter_set(i,:));

          tmperrs{i} = compute_error(tthis, trmodel, trd, detailed_data);
          tmpmaxerr(i) = max(tmperrs{i});
        end

        max_errs = tmpmaxerr;

        [dummy, max_mu_index] = max(tmpmaxerr);
        max_err_sequence = tmperrs{max_mu_index};

      else
        [max_errs, max_err_sequence, max_mu_index] = error_estimators(this, rmodel, detailed_data, parameter_set);
      end
    end

    function errs = compute_error(this, rmodel, reduced_data, detailed_data)
      % function errs = compute_error(this, rmodel, detailed_data)
      % @copybrief Greedy.Plugin.Interface.compute_error()
      %
      % @copydetails Greedy.Plugin.Interface.compute_error()
      %
      % Parameters:
      %   rmodel: of type .Greedy.User.IReducedModel
      %   detailed_data: of type .Greedy.User.IDetailedData
      %   reduced_data: of type .Greedy.User.ReducedData

      mu = get_mu(rmodel);
      dmodel = rmodel.detailed_model;
      set_mu(dmodel, mu);

      model_data = detailed_data.model_data;

      U       = generate(this.generator, dmodel, model_data);
      Uapprox = generate_reduced(this, rmodel, reduced_data, detailed_data, U);

      if isfield(model_data, 'grid')
        grid = model_data.grid;
      else
        grid = [];
      end

      dmodel = rmodel.detailed_model;

      if this.use_l2_error
        errs = dmodel.l2_error_sequence_algorithm(U, Uapprox, model_data);
        if this.relative_error
          errs = errs./dmodel.l2_error_sequence_algorithm(U, 0, model_data);
        end
      else
        errs = dmodel.linfty_error_sequence_algorithm(U, Uapprox, []);
        if this.relative_error
          errs = errs./dmodel.linfty_error_sequence_algorithm(U, 0, []);
        end
      end
    end

  end

  methods(Abstract)

    % function [max_errs, max_err_sequence, max_mu_index] = error_estimators(this, rmodel, detailed_data, M_train);
    % computes a posteriori error estimators for the reduced simulations for
    % every parameter vector from a given parameter set `{\cal
    % M}_{\text{train}}`.
    %
    % An a posteriori error estimator estimates an error
    %   ``\eta^k(\mu) \geq \| v_h(t^k;\mu) - v_{\text{red}}(t^k;\mu) \|``
    % for every time step `0\leq t^0, \ldots, t^K = T` and every parameter `\mu
    % \in {\cal M}`. The norm is a problem specific norm determined by the
    % options in 'rmodel'.
    %
    % This function's main use is to find a parameter vector
    % ``\mu_{\max} = \arg \sup_{\mu \in {\cal M}_{\text{train}} } \max_{k=0,\ldots,K} \eta^k(\mu).``
    %
    % @note The estimator must only depend on low dimensional data as computed
    %       by the prepare_reduced_data() method such that it is efficiently
    %       computable.
    %
    % Parameters:
    %   rmodel: an object of type .Greedy.User.IReducedModel specifying the
    %           basis generation process.
    %   detailed_data: of type .Greedy.User.IDetailedData
    %   M_train: a set of parameter vectors `{\cal M}_{\text{train}}` as
    %            returned by ParameterSampling.Interface.space .
    %
    % Return values:
    %   max_errs: a matrix of size 'n_parameters x model.nt+1' storing the
    %             error indicator `\eta^k(\mu)` for every `k=0,\ldots,K` and
    %             `\mu \in {\cal M}_{\text{train}}`.
    %   max_err_sequence: a sequence of error indicators `\eta^k(\mu_{\max})`
    %                     for every `k=0,\ldots,K`.
    %   max_mu_index: the index of the parameter vector `\mu_{\max}` in the
    %                 'parameter_set.set' matrix.
    %
    [max_errs, max_err_sequence, max_mu_index] = error_estimators(this, rmodel, detailed_data, M_train);

  end
end