EI.m

classdef EI < Greedy.Plugin.EICommon
  % Plugin for the Greedy.Algorithm class generating a collateral reduced basis
  % space plus interpolation DOFs and a local grid.
  %
  % This can be used for the empirical interpolation of parametrized functions
  % or operator evaluations.

  properties (Hidden=true)
    % temporary matrix generated by prepare_reduced_data()
    A = [];
  end

  properties (Constant)
    generated_basis_type = 'ei';

  end
  properties (Constant, Hidden=true)

    % cell array of field names that shall be copied to the generated
    % Greedy.DataTree.Detailed.EILeafNode instance.
    info_fields = { 'stop_Mmax', 'ei_target_error', 'noof_extensions', ...
      'minimum_residual', 'enable_regularisation' };
  end

  properties
    % maximum number of generated collateral reduced basis vectors.
    stop_Mmax                       = inf;

    % error mode for the empirical interpolation error
    %
    % This can be either 'interpol' or 'approx' for the true interpolation
    % error or the projection error of the operator interpolation error in the
    % collateral reduced basis space.
    ei_target_error                 = 'approx';

    % boolean flag indicating that the error indicator which is assumed to be
    % minimized is not equal to the ones returned by error_indicators()
    %
    % This is set to 'true' by the Greedy.Plugin.PODEI class.
    target_error_external           = false;

    % specifies how many collateral reduced basis functions shall be added in
    % each extension step.
    %
    % This makes only sense to be set greater than '1' for trajectories, i.e.
    % when more than one possible candidates for basis extension exist for one
    % parameter.
    noof_extensions                 = 1;

    % if the residual is below this barrier, the basis extension step is skipped.
    %
    % This option is only used in case of #target_error_external .
    minimum_residual                = 1e-6;

    % if this regularisation flag is set to true, the collateral reduced basis
    % spaces are forced to have zero mean, enforcing global conservation in
    % case of operator interpolation.
    enable_regularisation = false;
  end

  methods
    function eide = EI(generator)
      % constructor storing the snapshot generator
      %
      % Parameters:
      %   generator: snapshot generator of type SnapshotsGenerator.Cached.
      eide = eide@Greedy.Plugin.EICommon(generator);
    end

    function set.ei_target_error(this, value)
%      if isequal(value, 'linfty-interpol')
%        this.use_l2_error = 0;
%      else
%        this.use_l2_error = 1;
%      end
      this.ei_target_error = value;
    end

    function detailed_data = init_basis(this, rmodel, model_data, M_train)
      % function detailed_data = init_basis(this, rmodel, model_data, M_train)
      % @copybrief .Greedy.Plugin.Interface.init_basis()
      %
      % @copydetails Greedy.Plugin.Interface.init_basis()
      %
      % Parameters:
      %   rmodel: of type .Greedy.User.IReducedModel
      %
      % Return values:
      %   detailed_data: of type .Greedy.DataTree.Detailed.EILeafNode
      %
      % This method initializes the detailed data node.
      if isa(model_data, 'Greedy.DataTree.Detailed.INode')
        detailed_data = model_data;
        model_data    = detailed_data.model_data;
      else
        detailed_data = Greedy.DataTree.Detailed.EILeafNode(this, model_data);
        set_fields(detailed_data, this, Greedy.Plugin.EI.info_fields);
      end
      ophandle = this.generator.ophandle;
      W = ophandle.inner_product_matrix(model_data);

      % determine initial size of vectors:
      detailed_data.QM = zeros(size(W,1),0);
      detailed_data.TM = zeros(0,1);
      detailed_data.BM = zeros(0,0);

    end

    function prepare_reduced_data(this, rmodel, detailed_data)
      % function prepare_reduced_data(this, rmodel, detailed_data)
      % @copybrief GreedyPlugin::Interfaceprepare_reduced_data()
      %
      % @copydetails GreedyPlugin::Interfaceprepare_reduced_data()
      %
      % Parameters:
      %   rmodel: of type ::GreedyUser::IReducedModel
      %   detailed_data: of type ::GreedyDataTreeDetailed.EILeafNode
      %
      % This method pre-computes the normal solution matrix for
      % `L^2`-approximation in the collateral reduced space, in case
      % 'ei_target_error' is set to 'approx'.
      ophandle = this.generator.ophandle;
      if isequal(get_field(detailed_data, 'ei_target_error'), 'approx')
        W  = ophandle.inner_product_matrix(detailed_data.model_data);
        QM = detailed_data.QM;
        this.A = (QM' * W * QM)^(-1) * QM' * W;
      else
        this.A = [];
      end
    end

    function Uapprox = generate_reduced(this, rmodel, reduced_data, detailed_data, U)
      % function Uapprox = generate_reduced(this, [], [], detailed_data, U)
      % @copybrief GreedyPlugin::Interfacegenerate_reduced()
      %
      % @copydetails GreedyPlugin::Interfacegenerate_reduced()
      %
      % Parameters:
      %   rmodel: of type ::GreedyUser::IReducedModel
      %   detailed_data: of type ::GreedyUser::IDetailedData
      %   reduced_data: of type ::GreedyUser::IReducedData

      switch get_field(detailed_data, 'ei_target_error')
        case 'approx'
          P = this.A * U;
          Uapprox = detailed_data.QM*P;
        case {'interpol', 'linfty-interpol'}
          coeff = detailed_data.BM \ U(detailed_data.TM,:);
          Uapprox = detailed_data.QM * coeff;
        otherwise
          error('ei_target_error unknown!');
      end
    end

    function [breakloop, reason] = pre_check_for_end(this, rmodel, detailed_data)
      % function [breakloop, reason] = pre_check_for_end(this, rmodel, detailed_data)
      % @copybrief Greedy.Plugin.Interface.pre_check_for_end()
      %
      % @copydetails Greedy.Plugin.Interface.pre_check_for_end()
      %
      % Parameters:
      %   rmodel: of type .Greedy.User.IReducedModel
      %   detailed_data: of type .Greedy.DataTree.Detailed.EILeafNode
      breakloop = false;
      reason = '';
      dd_leaf = get_active_leaf(detailed_data, rmodel);

      if length(detailed_data.TM) >= get_field(dd_leaf, 'stop_Mmax')
        set_stop_flag(dd_leaf, 'stopped_on_max_basis_size');
        breakloop                          = true;
        reason = 'EI stopped on maximum basis size';
      elseif stopped_on_any_child(dd_leaf, 'stopped_on_empty_extension')
        breakloop = true;
        reason = 'EI stopped on empty extension';
      elseif stopped_on_all_leafs(dd_leaf, 'stopped_on_epsilon')
        breakloop = true;
        reason = 'EI stopped on epsilon';
      else
        max_errs = get_field(dd_leaf, 'max_err_sequence');
        % check for monotonicity, but not in case of external error estimation,
        % because this does not make any sense then.
        if length(max_errs) > 1 && (max_errs(end)>max_errs(end-1)) ...
           && isequal(get_field(dd_leaf, 'ei_target_error'), 'approx') ...
           && ~this.target_error_external
          set_stop_flag(dd_leaf, 'stopped_on_monotonicity');
          breakloop                        = true;
          reason = 'EI stopped on monotonicity';
        end
      end
    end

    function detailed_data = basis_extension(this, rmodel, detailed_data, max_err_seq, mu)
      % function detailed_data = basis_extension(this, rmodel, detailed_data, max_err_seq, mu)
      % @copybrief GreedyPlugin::Interfacebasis_extension()
      %
      % @copydetails GreedyPlugin::Interfacebasis_extension()
      %
      % Parameters:
      %   rmodel: of type ::GreedyUser::IReducedModel
      %   detailed_data: of type ::GreedyDataTreeDetailed.EILeafNode
      %
      % Return values:
      %   detailed_data: of type ::GreedyDataTreeDetailed.EILeafNode
      %
      dmodel = rmodel.detailed_model;
      set_mu(dmodel, mu);
      LU = generate(this.generator, dmodel, detailed_data.model_data);

      for ext_num = 1:get_field(detailed_data, 'noof_extensions')
        if ext_num >= 2
          % recompute error
          prepare_reduced_data(this, rmodel, detailed_data);
          max_err_seq = compute_error(this, rmodel, [], detailed_data);
        end

        [dummy, max_i] = max(max_err_seq);
        zeta_M         = LU(:,max_i(1));

        % residual
        r_M = zeta_M;

        QM = detailed_data.QM;
        TM = detailed_data.TM;
        BM = detailed_data.BM;

        m = size(QM, 2);
        %% compute interpolant
        if (m > 0)

          % get target values in interpolation points
          zeta_M_part = zeta_M(TM);

          % get interpolation coefficients
          %sigma = bicgstab(BM,zeta_M_part,[],1000);
          sigma = BM \ zeta_M_part;

          % compute residual
          r_M   = zeta_M - QM * sigma;
          if rmodel.verbose > 1
            disp(['positive minus negative contributions of residual (should be zero:)', num2str(sum(r_M))]);
            disp(['positive minus negative contributions of zeta_M   (should be zero:)', num2str(sum(zeta_M))]);
          end
%          if(sum(r_M) > 1e-4)
%            keyboard;
%          end
        end;
        max_r = max(abs(r_M));
        if (this.target_error_external || ext_num >=2)...
            && max_r < get_field(detailed_data, 'minimum_residual') * 10^(ext_num-1)

          disp(['skipping basis extension after ', num2str(ext_num), 'extension steps...']);
          ext_num = ext_num - 1;
          break;
        end
        t_M   = find(abs(r_M) >= max_r - 1e-12, 1);

        if ~isempty(t_M)
          if ~isempty(find(t_M(1)==TM, 1))
            % if target error is external, this error is not immediately fatal!
            if this.target_error_external
              disp('skipping EIDetailed basis extension because a magic point is already selected');
              disp('This might happen because of very bad parameter choice for ei extension.');
              disp(['Error at twice selected point is: ', num2str(max_r)]);
              disp('The error should be very low!');
              ext_num = ext_num - 1;
              break;
            end
            error('empirical interpolation selected magic point twice!!');
          end

          q_M = r_M / r_M(t_M(1)); % then sign is respected short version
          % (after verifying zero-1 structure...)
          %  BM = [BM, zeros(m-1,1) ; ...
          %       QM(t_M,:), 1 ];     % interpolation matrix mit BM(i,j) = q_j(t_i)
          %  TM = [TM; t_M];  % new vector of interpolation DOFs
          %  QM = [QM, q_M];  % new matrix of basis functions
          if rmodel.verbose > 1
            disp(['positive minus negative contributions of new basis function (should be zero:)', num2str(sum(q_M))]);
          end
          if abs(sum(q_M)) > 1e-12 && get_field(detailed_data, 'enable_regularisation')
            subset = abs(q_M) > 1e-8;
            if length(subset) > 10
              q_M(subset) = q_M(subset) - sum(q_M)/length(q_M(subset));
            end
            if rmodel.verbose > 1
              disp('WARNING: applied regularisation!');
              disp(['positive minus negative contributions of new basis function (should be zero:)', num2str(sum(q_M))]);
            end
          end


          BM = [BM,        q_M(TM) ; ...
               QM(t_M,:), q_M(t_M) ];   % interpolation matrix mit BM(i,j) = q_j(t_i)
          TM = [TM; t_M];                % new vector of interpolation DOFs
          QM = [QM, q_M];                % new matrix of basis functions

          detailed_data.BM = BM;
          detailed_data.TM = TM;
          detailed_data.QM = QM;


        else
          set_stop_flag(detailed_data, 'stopped_on_empty_extension');
          break;
        end
      end

      if ext_num > 0
        append_field(detailed_data, 'ei_ext_count', ext_num);
      end

    end

    function cop = copy(this)
      % function cop = copy(this)
      % deep copy of Greedy.Plugin.EI object
      %
      % Return values:
      %   cop: copied object of type Greedy.Plugin.EI
      cop = EI(this.generator);

      fnames = fieldnames(this);
      for i=1:length(fnames)
        cop.(fnames{i}) = this.(fnames{i});
      end
    end

  end
end