ei_error_convergence.m

function l = ei_error_convergence(detailed_data,detailed_path, ...
                                        ei_operator_path, params)
%function l = ei_error_convergence(detailed_data,detailed_path, ...
%                                       ei_operator_path, params)
%  
% function generating a error-convergence graph over the set 
% of parameters given in the paths. The handle to the plot is returned.
% The function performs the detailed_ei_interpolated simulation for
% increasing number of basis vectors M until the value given in params.M.
% (1) the maximum approximation error is plotted, i.e. projection of all 
% ei-snapshots on the increasing ei-spaces, and computation of the
% maximum L^infty(L2)-error, (2) the interpolation error,
% i.e. interpolation of all snapshots, error computation and maximum
% search and (3) the error between detailed-ei-simulation and the real
% detailed simulations is computed. 
% the detailed_path must contain detailed simulations, e.g. generated by 
% save_detailed_simulations, the ei_operator_path must contain the
% corresponding space-discretization operator results, e.g. generated by
% save_ei_operator_trajectories
  
% Bernard Haasdonk 7.9.2007

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%%%% general settings to be used
  Msamples = params.Msamples; 
  %  Mmax = size(detailed_data.BM,2);
  Mmax = params.Mmax; % size(detailed_data.BM,2);
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
%  detailed_path = params.RB_detailed_train_savepath;
%  ei_operator_path = params.ei_operator_savepath;

  errnames = {'approximation','interpolation','ei-simulation'}
  
  Ms = round((0:(Msamples-1)) * ...
             (Mmax-1)/(Msamples-1)+1);

  settings = load(fullfile(rbmatlabtemp,...
                           detailed_path,'settings.mat'));
  
  % storage for maximum l2-errors (i.e. linfty([0,T],l2)-norm)
  approx_errs = zeros(Msamples,1);
  interpol_errs = zeros(Msamples,1);
  ei_simulation_errs = zeros(Msamples,1);
  
  for mu_ind = 1:size(settings.M,2);
%  for mu_ind = 1:1;
    disp(['processing parameter vector ',num2str(mu_ind),...
          '/',num2str(size(settings.M,2))]);
    params = set_mu(settings.M(:,mu_ind),params);
    if ismember('ei-simulation',errnames)
      sim_data = load_detailed_simulation(mu_ind,settings.savepath,params);
      U_H = sim_data.U;
    end;
    for Mind = 1:Msamples
      fprintf('.');
      M = Ms(Mind)
      params.M = M;
      
      % detailed simulation with ei of the operator
      if ismember('ei-simulation',errnames)
        U_ei_H = detailed_ei_nonlin_evol_simulation(detailed_data, ...
                                                    params);
        
        l2err = fv_l2_error(U_H,U_ei_H,detailed_data.W);
        max_l2err = max(l2err);
        if max_l2err > ei_simulation_errs(Mind)
          ei_simulation_errs(Mind) = max_l2err;
        end;
        clear('U_ei_H');
      end;
      
      fname = ['LU',num2str(mu_ind)];
      tmp = load(fullfile(rbmatlabtemp,ei_operator_path,fname));
      L_E_U = tmp.LU;
      clear('tmp');
      % approximation error
      if ismember('approximation',errnames)
        % 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
        W = feval(params.get_inner_product_matrix(detailed_data),...
                  detailed_data.grid,...
                  params);
        A = (detailed_data.QM(:,1:M)' * W * ...
             detailed_data.QM(:,1:M))^(-1) * ...
            detailed_data.QM(:,1:M)' * W;
        
        L_E_U_appr = detailed_data.QM(:,1:M) * A * L_E_U;
        
        l2err = fv_l2_error(L_E_U_appr,L_E_U,detailed_data.W);
        max_l2err = max(l2err);
        if max_l2err > approx_errs(Mind)
          approx_errs(Mind) = max_l2err;
        end;
        clear('L_E_U_appr');
      end;
      
      % interpolation error
      if ismember('interpolation',errnames)
        % for each discrete function in LU
        % compute interpolation in current space spanned by Q
        W = params.get_inner_product_matrix(detailed_data);
        
        QM = detailed_data.QM(:,1:M);
        
        % get interpolation coefficients
        coeff = detailed_data.BM(1:M,1:M) \ L_E_U(detailed_data.TM(1:M),:);
        L_E_U_interpol = QM * coeff;
        
        l2err = fv_l2_error(L_E_U_interpol,L_E_U,detailed_data.W);
        max_l2err = max(l2err);
        if max_l2err > interpol_errs(Mind)
          interpol_errs(Mind) = max_l2err;
        end;
        clear('L_E_U_interpol');
      end;
    end;
    fprintf('.');
    % save workspace
    disp('workspace is saved in ws.mat in case of break....');
    save('ws');  
  end; 
  disp('finished!!!');
  l = plot(Ms,[approx_errs,interpol_errs,ei_simulation_errs]);
  set(l,'Linewidth',2);
  set(l(1),'LineStyle','-');
  set(l(2),'LineStyle','-.');
  set(l(3),'LineStyle',':');  
  set(gca,'Yscale','log');
  legend({'max approx-error','max interpol-error','max-ei-simulation-error'});
  title('empirical interpolation error'); 
  xlabel('M');
  ylabel('L2-error');
  % save workspace
  disp('workspace is saved in ws.mat for modification of plots.');
  save('ws');  
% TO BE ADJUSTED TO NEW SYNTAX
%| \docupdate