EIPOD.m

classdef EIPOD < Greedy.Plugin.EI
  % 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.

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

    function Uapprox = generate_reduced(this, rmodel, reduced_data, detailed_data, U)
      % function Uapprox = generate_reduced(this, [], [], detailed_data, U)
      % @copybrief Greedy.Plugin.Interface.generate_reduced()
      %
      % @copydetails Greedy.Plugin.Interface.generate_reduced()
      %
      % Parameters:
      %   rmodel: of type .Greedy.User.IReducedModel
      %   detailed_data: of type .Greedy.User.IDetailedData
      %   reduced_data: of type .Greedy.User.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 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 Greedy.Plugin.Interface.basis_extension()
      %
      % @copydetails Greedy.Plugin.Interface.basis_extension()
      %
      % Parameters:
      %   rmodel: of type .Greedy.User.IReducedModel
      %   detailed_data: of type .Greedy.DataTree.Detailed.EILeafNode
      %
      % Return values:
      %   detailed_data: of type .Greedy.DataTree.Detailed.EILeafNode
      %
      dmodel = rmodel.detailed_model;
      set_mu(dmodel, mu);
      LU = generate(this.generator, dmodel, detailed_data.model_data);
      LU_red = generate_reduced(this, rmodel, [], detailed_data, LU);

      if size(LU,1) == detailed_data.model_data.grid.nelements
        A = detailed_data.model_data.W;
      else
        A = detailed_data.model_data.diamondW;
      end
      QM = detailed_data.QM;
      zeta_M = PCA_fixspace(LU - LU_red, QM, A, 1, 'qr',eps);

      % residual
      r_M = zeta_M;

      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 and project it onto QM^orthogonal
        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));
      t_M   = find(abs(r_M) >= max_r - 1e-12, 1);

      if ~isempty(t_M)
        if ~isempty(find(t_M(1)==TM, 1))
          error('empirical interpolation selected magic point twice!!');
        end

        q_M = PCA_fixspace(r_M, QM, A, 1, 'qr', eps);
%        q_M = orth_r_M / orth_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');
      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