ei_detailed.m

Go to the documentation of this file.

function ei_detailed_data = ei_detailed(model, model_data, LU_fnames, params)
%function ei_detailed_data = ei_detailed(model, model_data, LU_fnames, params)
% constructs a collateral reduced basis and interpolation points for given
% operator evaluations.
%
% This function prepares the empirical interpolation for an arbitrary
% sequence of DOF-vectors. The DOF-vectors are expected to be
% stored in the file(s) given in 'LU_fnames', which is a cell array of strings
% or a string array.  The files are expected to contain a single matrix 'LU',
% which is a set of columnwise DOF-vectors. Usually the DOF-vectors correspond
% to operator evaluations, which are e.g. computed by
% ei_operator_collect_files().
%
% \attention Every field in the returned structure 'ei_detailed_data' is a
%   cell array, and this function only manipulates entries indexed by
%   '(params.index, params.time_index)', e.g.
%   'ei_detailed_data.QM{params.index, params.time_index}' holds the basis
%   functions for the empirical interpolation of the operator with the index
%   'params.index' for arguments lying in the time slice indexed by
%   'params.time_index'. As this function extends data given in 'model_data',
%   the empirical interpolation can such be computed for different operators
%   and time slices.
%
% optional fields of params:
%     ei_Mmax                 :  integer indicating, when stop is wanted after
%                                reaching the maximum number of collateral RB
%                                vectors (default = inf)
%     ei_target_error         :  string which is one of
%           - 'approx' for L2-projection error,
%           - 'interpol' for L2-interpolation error or
%           - 'linfty-interpol' for `L^{\infty}` interpolation error
%           .
%                                (default = 'approx')
%     ei_stop_epsilon         :  interpolation error limit stopping the basis
%                                extension when reached. (Default: 0)
%     time_index              :  index of time slice for which the empirical
%                                interpolation information is generated. See
%                                time splitting for details. (default = 1)
%     index                   :  an operator index that can be associated with
%                                the generated empirical interpolation data.
%     compute_lebesgue        :  compute the Lebesgue constant for the
%                                interpolation (default = false)
%     skip_search             :  boolean flag indicating whether the search for
%                                new basis vectors shall be skipped after a
%                                time domain split. If set to false, the basis
%                                vectors from the previous CRB step chosen for
%                                this time domain are used.  (default = false)
%
% optional fields of model_data:
%     time_split_map          :  a map between a time step index `t \in \{1,
%                                \ldots, model.nt+1\}` and a time slice index
%                                `s \in \{1, \ldots, S_{\max}\}`, where
%                                `S_{\max}` is the number of time slices.
%                                For more details look at documentation of
%                                'params.LU_fnames'.  (default = `t\mapsto 1
%                                ` for all `t`).
%
% parameters:
%     LU_fnames               :  a cell array of strings or an array of strings
%                                of file names relative to rbmatlabtemp(),
%                                where time series of DOF vectors are stored
%                                for which the collateral reduced basis space
%                                shall be trained. Usually, 'LU_fnames' is
%                                generated by ei_operator_collect_files(). The
%                                files should contain an array named 'LU' of
%                                dimension 'dofs x model.nt+1' or
%                                'dofs x model.nt+1+\#intermediate newton steps'
%                                In the first case, only the vectors 'LU(:,s)'
%                                with 'model_data.time_split_map(s)==params.time_index'
%                                are used for empirical interpolation.
%     params                  :  a struct containing control parameters for the
%                                empirical interpolation.
%
% return values:
%     ei_detailed_data        :  a struct containing the vectors and matrices
%                                storing all high dimensional data needed for
%                                empirical interpolation, plus some information
%                                on the basis generation. Note, that this might
%                                be an extension the parameter 'model_data'.
%
% generated fields of ei_detailed_data:
%     QM: a matrix of columnwise DOF-vectors `q_j` for interpolation
%         to be used as the colateral basis
%     TM: a vector for the set of 'magic points' by containing the index
%         numbers of the corresponding DOF-nodes (In case of piecewise constant
%         or linear basis-functions, the maxima `x_i` always can be found in a
%         cell centroid (deg=0) or a node (deg=1).
%     BM: the corresponding interpolation-matrix of dimension 'Mmax x Mmax',
%         i.e. interpolation can be done by solving the equation system `B
%         \sigma = \left(\zeta(x_1), \ldots, \zeta(x_M)\right)` then `Q
%         \sigma` is the empirical interpolation. Note, that this
%         reconstruction may not be done in the online phase!
%     ei_info: structure giving special details on creation:
%           - 'max_err_sequence': a vector of errors, AFTER the
%                corresponding new basis vector extension,
%           - 'extension_filenum': a vector with the filenumbers, from which
%                the snapshots are selected
%           - 'extension_filepos': a vector with the relative position of the
%                snapshot within its file extension_filenum. By the latter two
%                fields, a unique identification of the selected snapshots is
%                possible, i.e. a later translation into parameter vectors,
%                time, etc.
%           - 'stopped_on_Mmax': flag indicating, whether stop on maximum
%               'params.ei_Mmax' is a reason for termination
%           - 'stopped_on_monotonicity': flag indicating, that there is no
%                error decrease anymore, which incicates an accuracy problem,
%                so termination is performed
%           - 'stopped_on_duplicate_point': flag indicating, that the
%                empirical interpolation selected the same interpolation point
%                twice for interpolation, which indicates an accuracy limit
%           - 'stopped_on_epsilon': flag indicating, that the empirical
%                interpolation stopped because the interpolation error dropped
%                below the given limit 'params.ei_stop_epsilon'
%
%

% Bernard Haasdonk 15.8.2007

%ei_target_error = 'interpol';


grid = model_data.grid;

if ~isfield(model_data, 'time_split_map')
  model_data.time_split_map = [1:model.nt+1; ones(size(1:model.nt+1)) ]';
end

if ~isfield(params,'time_index')
  params.time_index = 1;
end

if ~isfield(params,'ei_target_error')
  params.ei_target_error = 'approx';
end;

if ~isfield(params,'ei_stop_epsilon')
  params.ei_stop_epsilon = 0;
end;

if ~isfield(params,'skip_search')
  params.skip_search = false;
end

%if ~isfield(params,'ei_rel_abs')
%  params.ei_rel_abs = 'absolute';
%end;

if ~isfield(params,'ei_Mmax')
  params.ei_Mmax = inf;
end;

if ischar(LU_fnames)
  LU_fnames = {LU_fnames};
end;

if ~isfield(params,'index')
  params.index = 1;
end

if ~isfield(params,'compute_lebesgue')
  params.compute_lebesgue = false;
end

pari = params.index;
ti   = params.time_index;

% get weight matrix for inner-product computation
W = model.get_inner_product_matrix(model_data);

% determine initial size of vectors:
Q = zeros(size(W,1),0);
T = zeros(0,1);
B = zeros(0,0);

% initialize RB_info data
max_errs          = [];
extension_filenum = [];
extension_filepos = [];

m = 0;
% difference to ei from literature:
% instead of simply first vector and normalization => search worst l2-error
% vector
if isfield(model_data, 'QM') ...
    && all( [pari,ti] <= size(model_data.QM) ) ...
    && size(model_data.QM{pari,ti},2) > 0
  m = size(model_data.QM{pari,ti},2);
  if(model.verbose > 0)
    disp(['Extending collateral reduced basis space of size ', num2str(m), '.']);
  end
  if params.skip_search
    Q = model_data.QM{pari,ti};
    T = model_data.TM{pari,ti};
    B = model_data.BM{pari,ti};
  end
  tei_info          = model_data.ei_info{pari,ti};
  max_errs          = [ Inf; tei_info.max_err_sequence ];
  extension_filenum = tei_info.extension_filenum;
  extension_filepos = tei_info.extension_filepos;
end

continue_loop              = 1;
stopped_on_Mmax            = 0;
stopped_on_monotonicity    = 0;
stopped_on_epsilon         = 0;
stopped_on_duplicate_point = 0;
if params.ei_Mmax < 1
  continue_loop = 0;
end;

use_l2_error = 0;
if ismember(params.ei_target_error,{'approx','interpol'})
  disp('using l2-error for ei offline')
  use_l2_error = 1;
else
  disp('using linfty-error for ei offline')
end;


while continue_loop
  m = m+1;

  disp(['determining ei-basis & magic point ',num2str(m),':']);
  % for each discrete function in LU
  % compute best-l2-approximation in current space spanned by Q
  %
  % if <u_H,v_H> = u' * W * v    (u,v coefficient vectors)
  % then min_p (u_H - sum p(i) q_i) is minimized by
  %
  % p = (Q' W Q)^(-1) Q^T W u
  %
  % So compute whole array of p-vectors and resulting approximation Q*p

  if params.skip_search && m <= length(extension_filenum)
    tmp = load(fullfile(rbmatlabtemp,LU_fnames{extension_filenum(m)}));
    if size(tmp.LU, 2) == tmodel.nt+1
      tslice = tmodel_data.time_split_map(:,2) == tparams.time_index;
      tslice = tmodel_data.time_split_map(tslice,1);
      LU = tmp.LU(:,tslice);
    else
      LU = tmp.LU;
    end
    zeta_M = tmp.LU(:,extension_filepos(m));
  else
    if strcmp(params.ei_target_error, 'approx') == 1
      A = (Q' * W * Q)^(-1) * Q' * W;
    else
      A = [];
    end

    tmp_max_is   = zeros(length(LU_fnames),1);
    tmp_max_errs = zeros(length(LU_fnames),1);
    tmp_zeta_M   = zeros(grid.nelements,length(LU_fnames));
    parfor i=1:length(LU_fnames)
      tmp         = load(fullfile(rbmatlabtemp,LU_fnames{i}));
      tmodel_data = model_data;
      tmodel      = model;
      tparams     = params;

      if size(tmp.LU, 2) == tmodel.nt+1
        tslice = tmodel_data.time_split_map(:,2) == tparams.time_index;
        tslice = tmodel_data.time_split_map(tslice,1);
        LU = tmp.LU(:,tslice);
      else
        LU = tmp.LU;
      end
      LUapprox=[];
      switch tparams.ei_target_error
       case 'approx' % approximation, i.e. L2-projection error
        P = A * LU;
        LUapprox = Q*P;
       % interpolation, i.e. empirical interpolation error
       case {'interpol','linfty-interpol'}
        coeff = B \ LU(T,:);
        LUapprox = Q * coeff;
       otherwise
        error('ei_target_error unknown!')
      end;
      % search maximum l2-error
      if use_l2_error
        errors = tmodel.l2_error_sequence_algorithm(LU,LUapprox,model_data.W);
      else % use infty
        errors = tmodel.linfty_error_sequence_algorithm(LU,LUapprox,[]);
      end;
      % choose maximum error-snapshot for CRB extension
      [tmp_max_err, tmp_max_i] = max(errors);
      tmp_max_is(i)            = tmp_max_i(1);
      tmp_max_errs(i)          = tmp_max_err(1);
      tmp_zeta_M(:,i)          = LU(:,tmp_max_i(1));
    end;


    % get largest from filewise extremes
    [max_err, max_in] = max(tmp_max_errs);
    disp(['detected maximum error for extension: ',num2str(max_err(1))]);
    max_errs          = [max_errs; max_err(1)];
    zeta_M            = tmp_zeta_M(:,max_in(1));
    extension_filenum = [extension_filenum; max_in(1)];
    extension_filepos = [extension_filepos; tmp_max_is(max_in(1))];

  end

  r_M = zeta_M;

  if (m > 1) % compute interpolant
    zeta_M_part = zeta_M(T); % get target values in interpolation points
    %sigma = bicgstab(B,zeta_M_part,[],1000); % get interpolation
    %coefficients
    sigma = B \ zeta_M_part;
    r_M   = zeta_M - Q * sigma; % compute residual
    if model.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, t_M] = max(abs(r_M));
  t_M          = find(abs(r_M) >= max_r - 1e-12,1);
  if ~isempty(find(t_M(1)==T, 1))
    disp('empirical interpolation selected magic point twice!!');
    stopped_on_duplicate_point = 1;
    continue_loop = 0;

  elseif (length(max_errs)>1) && (max_errs(end)>max_errs(end-1)) ...
        && isequal(params.ei_target_error,'approx')
    % check that l2-error sequence is decreasing in case of
    % approximation error (interpol-error can be
    % non-monotonic...). If not: accuracy
    % problem and interpolation should be stopped
    % stop empirical interpolation due to accuracy limit
    disp('ei detected accuracy limit, loop is stopped.')
    continue_loop           = 0;
    stopped_on_monotonicity = 1;
  elseif (length(max_errs) > 1) && (max_errs(end) < params.ei_stop_epsilon)
    disp('ei detected accuracy limit, loop is stopped.')
    continue_loop      = 0;
    stopped_on_epsilon = 1;
  else
    q_M = r_M/r_M(t_M(1)); % then sign is respected short version
    % (after verifying zero-1 structure...)
    %  B = [B, zeros(m-1,1) ; ...
    %       Q(t_M,:), 1 ];     % interpolation matrix mit B(i,j) = q_j(t_i)
    %  T = [T; t_M];  % new vector of interpolation DOFs
    %  Q = [Q, q_M];  % new matrix of basis functions
    B = [B,        q_M(T) ; ...
         Q(t_M,:), q_M(t_M) ];     % interpolation matrix mit B(i,j) = q_j(t_i)
    T = [T; t_M];  % new vector of interpolation DOFs
    Q = [Q, q_M]; % new matrix of basis functions
    if model.verbose > 1
      disp(['positive minus negative contributions of new basis function (should be zero:)', num2str(sum(q_M))]);
    end
  end;

  if params.ei_Mmax <= length(T)
    disp('ei reached prescribed number of CRB-vectors, loop is stopped')
    continue_loop   = 0;
    stopped_on_Mmax = 1;
  end;
end; % loop

A = (Q' * W * Q)^(-1) * Q' * W;

% determine final error
tmp_max_errs = zeros(1, length(LU_fnames));
parfor i=1:length(LU_fnames)
  tmodel      = model;
  tmp         = load(fullfile(rbmatlabtemp,LU_fnames{i}));
  tmodel_data = model_data;
  tparams     = params;

  tslice   = tmodel_data.time_split_map(:,2) == tparams.time_index;
  tslice   = tmodel_data.time_split_map(tslice,1);
  LU       = tmp.LU(:,tslice);
  LUapprox = [];
  switch params.ei_target_error
   case 'approx' % approximation, i.e. L2-projection error
    P        = A * LU;
    LUapprox = Q * P;
   case {'interpol','linfty-interpol'} % L2-/Linfty-interpolation error
    coeff    = B \ LU(T,:);
    LUapprox = Q * coeff;
   otherwise
    error('ei_target_error unknown!')
  end;

  % search maximum l2- or linfty error
  if use_l2_error
    tmp_errors = tmodel.l2_error_sequence_algorithm(LU,LUapprox,model_data.W);
    % choose maximum l2error-snapshot for CRB extension
  else % use linfty_error
    tmp_errors = tmodel.linfty_error_sequence_algorithm(LU,LUapprox,[]);
  end;
  [tmp_max_err, tmp_max_i] = max(tmp_errors);
  tmp_max_errs(i)          = tmp_max_err;
end;

max_err = max(tmp_max_errs);

%params.title = 'Interpolation basis functions q_m';
%plot_element_data_sequence(grid,Q,params);
%disp(['Interpolation matrix B should have lower-triangular structure,' ...
%     'one-diagonal:']);
%disp(B);
%u = zeros(grid.nelements,1);
%u(T) = 1;
%plot_element_data(grid,u,params);
%title('Interpolation points/DOFS');
%keyboard;

% return results
ei_detailed_data = model_data;
if ~isfield(ei_detailed_data,'QM')
  ei_detailed_data.QM      = cell(1,1);
  ei_detailed_data.TM      = cell(1,1);
  ei_detailed_data.BM      = cell(1,1);
  ei_detailed_data.ei_info = cell(1,1);
end

ei_detailed_data.QM{pari,ti}                                    = Q;
ei_detailed_data.TM{pari,ti}                                    = T;
ei_detailed_data.BM{pari,ti}                                    = B;
ei_detailed_data.ei_info{pari,ti}                               = [];
ei_detailed_data.ei_info{pari,ti}.max_err_sequence              = [max_errs(2:end); max_err];
ei_detailed_data.ei_info{pari,ti}.extension_filenum             = extension_filenum;
ei_detailed_data.ei_info{pari,ti}.extension_filepos             = extension_filepos;
ei_detailed_data.ei_info{pari,ti}.stopped_on_monotonicity       = stopped_on_monotonicity;
ei_detailed_data.ei_info{pari,ti}.stopped_on_epsilon            = stopped_on_epsilon;
ei_detailed_data.ei_info{pari,ti}.stopped_on_Mmax               = stopped_on_Mmax;
ei_detailed_data.ei_info{pari,ti}.ei_stopped_on_duplicate_point = ...
    stopped_on_duplicate_point;
if params.compute_lebesgue
  XI                                             = Q / B;
  ei_detailed_data.ei_info{pari,ti}.lebesgue     = max(sum(abs(XI),2));
  ei_detailed_data.ei_info{pari,ti}.max_lebesgue = 2^length(T) - 1;
end