rb_tutorial.m

function res = rb_tutorial(step)
%function res = rb_tutorial(step)
%
% step : index of experiment to be performed
%        1  - 10: Thermal block, elliptic PDE
%        11 - 17: Advection-Diffusion, parabolic PDE
%
% This program allows to reproduce the thermalblock and advection-diffusion 
% experiments of the RB-Summerschool and the pictures are used in the 
% RB-Tutorial book chapter.
%
% B. Haasdonk: Reduced Basis Methods for Parametrized PDEs --  
% A Tutorial Introduction for Stationary and Instationary Problems.
% Chapter in P. Benner, A. Cohen, M. Ohlberger and K. Willcox (eds.): 
% "Model Reduction and Approximation: Theory and Algorithms", SIAM, Philadelphia, 2016
%
% A Preprint version of this chapter is available at:
%
% http://www.simtech.uni-stuttgart.de/publikationen/prints.php?ID=938

% B. Haasdonk, 7.9.2013

if nargin<1
  step = 1
  help rb_tutorial;
  disp('now running step 1');
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%PART 1: Stationary Thermalblock example
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

switch step
 case 1 
  disp('thermal block for different configurations and corresponding plots');
  params = [];
  params.B1 = 2;
  params.B2 = 2;
  params.mu_range = [0.1;10];
  params.numintervals_per_block = 5;
  model = thermalblock_model_struct(params);
  model_data = gen_model_data(model);
  plot_params = [];
  plot_params.axis_equal = 1;
  plot_params.axis_tight = 1;
  
  % plot 1: four identical values
  figure;
  model = set_mu(model,[1,1,1,1]);
  sim_data = detailed_simulation(model,model_data);
  plot_sim_data(model,model_data,sim_data,plot_params);

  % plot 2: lower and upper values identical
  figure;
  model = set_mu(model,[0.1,0.1,10,10]);
  sim_data = detailed_simulation(model,model_data);
  plot_sim_data(model,model_data,sim_data,plot_params);

  % plot 3: four different values
  figure;
  model = set_mu(model,[0.3,0.2,0.2,0.7]);
  sim_data = detailed_simulation(model,model_data);
  plot_sim_data(model,model_data,sim_data,plot_params);

  % plot 4: larger division
  %  mu = rand(params.B1*params.B2,1)+0.1;
  figure;
  mu = [ 0.8094    0.8547    0.3760    0.7797    0.7551    0.2626, ...
         0.2190    0.5984    1.0597    0.4404    0.6853    0.3238, ...   
         0.8513    0.3551    0.6060    0.7991    0.9909    1.0593 ...    
         0.6472    0.2386    0.2493    0.3575    0.9407    0.3543 ...
         0.9143    0.3435    1.0293    0.4500    0.2966    0.3511 ...
         0.7160    0.5733    0.4517    0.9308    0.6853    0.6497]';
  params = [];
  params.B1 = 6;
  params.B2 = 6;
  params.mu_range = [0.1;10];
  model = thermalblock_model_struct(params);
  model_data = gen_model_data(model);
  model = set_mu(model,mu);
  sim_data = detailed_simulation(model,model_data);
  plot_sim_data(model,model_data,sim_data,plot_params);
    
 case 2 
  disp('thermal block detailed gui');
  params = [];
  params.B1 = 2;
  params.B2 = 2;
  params.numintervals_per_block = 50;
  params.mu_range = [0.1;10];
  params.numintervals_per_block = 5;
  model = thermalblock_model_struct(params);
  model_data = gen_model_data(model);
  plot_params = [];
  plot_params.axis_equal = 1;
  plot_params.axis_tight = 1;
  plot_params.yscale_uicontrols = 0.7;
  demo_detailed_gui(model,model_data,[],plot_params);
  set(gcf,'Position',[584   245   553   691]);
  
 case 3 
  disp('equidistant samples for basis generation, offline versus online time');
  
  params = [];
  params.B1 = 2;
  params.B2 = 2;
  params.numintervals_per_block = 20;
  params.mu_range = [0.1;10];
  disp('initializing model')
  model = thermalblock_model_struct(params);
 
  disp('-----------------------------------------------------')
  disp('OFFLINE PHASE:');
  disp('generating model data: grid, fem inner product matrix, etc.');
  tic
  model_data = gen_model_data(model);
  model.RB_generation_mode = 'lagrangian';
  model.RB_mu_list = lin_equidist_samples(model,[4,2,2,2]);
  disp('generating detailed data: A_q,f_q,l_q components, reduced basis.');
  detailed_data = gen_detailed_data(model,model_data);
  disp('generating reduced data: A_Nq,f_Nq,l_Nq components');
  reduced_data = gen_reduced_data(model, detailed_data);
  t_offline = toc;
  disp('Offline time: ')
  t_offline
  save('thermalblock_2x2_data');

  disp('-----------------------------------------------------')
  disp('ONLINE PHASE:');
  % online phase for single simulation: 
  disp('reduced simulation: assembling system and solve');
  model = set_mu(model,[1,1,1,1]);
  tic
  sim_data = rb_simulation(model,reduced_data);
  t_online = toc;

  disp('Online time: ')
  t_online

  disp('starting reduced basis gui:')
  plot_params = [];
  plot_params.axis_equal = 1;
  plot_params.axis_tight = 1;
  plot_params.yscale_uicontrols = 0.7;
  demo_rb_gui(model,detailed_data,reduced_data,plot_params);
  
 case 4 
  disp('1-parameter example: error, estimator, effectivity bound plot.');
  
  % generate reduced model 
  disp('generating reduced model')
  params = [];
  params.B1 = 2;
  params.B2 = 2;
  params.numintervals_per_block = 20;
  params.mu_range = [0.1;10];
  model = thermalblock_model_struct(params);  
  model_data = gen_model_data(model);
  model.RB_generation_mode = 'lagrangian';
  c = 0.1;
  mu_list = {[0.1,c,c,c],...
             [0.5,c,c,c],...
             [0.9,c,c,c],...
             [1.3,c,c,c],...
             [1.7,c,c,c]};
  model.RB_mu_list = mu_list;
  detailed_data = gen_detailed_data(model,model_data);
  reduced_data = gen_reduced_data(model,detailed_data);

  disp('computing parameter sweep')
  mu1s = linspace(0.1,4,1000);
  % rapidly computable error landscape:
  Deltas = zeros(length(mu1s),1);
  for n = 1:length(mu1s)
    fprintf('.');
    mu = [mu1s(n),c,c,c];
    model = set_mu(model,mu);
    sim_data = rb_simulation(model,reduced_data);
    Deltas(n) = sim_data.Delta;
  end;
  fprintf('\n');
  plot(mu1s,Deltas,'Linewidth',2);
  xlabel('mu1','Fontsize',15);
  ylabel('error/estimator','Fontsize',15);
  title('error estimator by parameter sweep','Fontsize',15);
  set(gca,'Yscale','log');
  legend('estimator \Delta_u(\mu)');
  
  % comparison with expensive true error 
  
  disp('computing true error for some parameters')
  mu1s = linspace(0.1,4,40);
  errs = zeros(length(mu1s),1);
  for n = 1:length(mu1s)
    fprintf('.');
    mu = [mu1s(n),c,c,c];
    model = set_mu(model,mu);
    sim_data = detailed_simulation(model,model_data);
    rb_sim_data = rb_simulation(model,reduced_data);
    % linear combination of reduced vector with basis:
    rb_sim_data= rb_reconstruction(model,detailed_data,rb_sim_data);
    err = sim_data.uh;
    err.dofs = err.dofs - rb_sim_data.uh.dofs;
    errs(n) = fem_h10_norm(err);
  end;
  fprintf('\n');
  hold on;
  plot(mu1s,errs,'r*','Markersize',10);
  legend({'estimator \Delta_u(\mu)','true error |u-u_N|'},'Fontsize',15);

  %set(gcf,'Position', [379   492   434   328]);
  %saveas(gcf,'error_parameter_sweep.png');
  %saveas(gcf,'error_parameter_sweep.eps','epsc');

  save('rb_tutorial_step4');
  
 case 5  
  disp('Effectivity plot');
  load('rb_tutorial_step4');
  
  disp('computing parameter sweep')
  mu1s = linspace(0.1,4,40);
  % rapidly computable error landscape:
  Deltas = zeros(length(mu1s),1);
  for n = 1:length(mu1s)
    fprintf('.');
    mu = [mu1s(n),0.1,0.1,0.1];
    model = set_mu(model,mu);
    sim_data = rb_simulation(model,reduced_data);
    Deltas(n) = sim_data.Delta;
  end;
  fprintf('\n');

  etas = Deltas.*(errs.^(-1));
  i = find(errs<1e-7);
  %i = find(etas>20);
  etas(i)=NaN;
  plot(mu1s,etas,'b*','Markersize',10);
  xlabel('\mu1','Fontsize',15);
  ylabel('effectivity','Fontsize',15);
  title('effectivity over parameter','Fontsize',15);
%  set(gca,'Yscale','log');
  
  hold on;
  gammas = mu1s;
  alpha = 0.1;
  effectivities = gammas/alpha;
  plot(mu1s,effectivities,'r-','Linewidth',2);
  
  legend({'effectivity \Delta_u/|u-u_N|',...
          'upper bound \gamma(\mu)/\alpha(\mu)'},...
         'Fontsize',15,...
         'Location','NorthWest');
  
%  set(gcf,'Position', [379   492   434   328]);
%  saveas(gcf,'effectivity_parameter_sweep.png');
%  saveas(gcf,'effectivity_parameter_sweep.eps','epsc');
  
 case 6 
  disp('Convergence plot for test error on equidistant samples');

    params = [];
    params.B1 = 3;
    params.B2 = 3;
    params.numintervals_per_block = 20;
    params.mu_range = [0.5;2]; 
    model = thermalblock_model_struct(params);  
    model_data = gen_model_data(model);
    model.RB_generation_mode = 'model_RB_basisgen';
    model.RB_basisgen = @lagrangian_orth_basisgen;
    Ns = 2:8; % linear dependent for higher values...
    
    maxDeltas = zeros(length(Ns),1);
    maxerrs = zeros(length(Ns),1);
    
    disp('Precomputing test snapshots');
    ntest = 100;
    [U,test_mus] = filecache_function(@precompute_test_snapshots,...
                                      params,model,model_data,ntest);
       
    disp('Generating basis and determining max error, estimators');
    tic;
    for Ni = 1:length(Ns)
      N = Ns(Ni)
      train_mu1s = linspace(params.mu_range(1),params.mu_range(2),N);
      train_mus = [train_mu1s(:), ones(N,params.B1*params.B2-1)*1];      
      model.RB_mu_list = mat2cell(train_mus,ones(N,1),params.B1*params.B2);
      detailed_data = gen_detailed_data(model,model_data);
      reduced_data = gen_reduced_data(model,detailed_data);    
      for i = 1:ntest
        model = set_mu(model,test_mus(i,:));
        rb_sim_data = rb_simulation(model, reduced_data);
        rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
        err = rb_sim_data.uh;
        err.dofs = err.dofs - U(:,i);
        errnorm = fem_h10_norm(err);
        if maxerrs(Ni)<errnorm
          maxerrs(Ni)=errnorm;
        end;
        if maxDeltas(Ni)<rb_sim_data.Delta
          maxDeltas(Ni)=rb_sim_data.Delta;
        end;
      end; % i
    end; % Ni
    
    time_step_6 = toc  
    save('rb_tutorial_step6');
    % case 6.1 % plot of error convergenc 
    %  load('rb_tutorial_step6');
    plot(Ns,maxDeltas,'b-','Linewidth',2);
    hold on
    plot(Ns,maxerrs,'r:','Linewidth',2);
    set(gca,'Yscale','log');
    legend({'estimator \Delta_u(\mu)','error |u(\mu)-u_N(\mu)|'},'Fontsize',15);
    title('test error decay for growing sample size','Fontsize',15);
    xlabel('sample/basis Size N','Fontsize',15);
    ylabel('error/estimator','Fontsize',15);
    
    %  set(gcf,'Position', [379   492   434   328]);
    %  saveas(gcf,'lagrangian_error_convergence.png');
    %  saveas(gcf,'lagrangian_error_convergence.eps','epsc');

 case 7 
  disp('plot basis');
  load('rb_tutorial_step6');
  plot_params = [];
  plot_params.axis_equal = 1;
  plot_params.axis_tight = 1;
  plot_params.plot = @plot_vertex_data;
  plot_params.title='Orthonormal Reduced Basis \Phi_N';
  % plot as sequence
  plot_sequence(detailed_data.RB,detailed_data.grid,plot_params);
    
  % or in single plot
  plot_params.show_colorbar = 0;
  plot_params.no_lines = 1;
  figure;
  for n = 1:8
    subplot(2,4,n);
    plot_vertex_data(detailed_data.grid,detailed_data.RB(:,n), ...
                     plot_params);
    title(['n=',num2str(n)],'Fontsize',10);
    set(gca,'Xtick',[]);
    set(gca,'Ytick',[]);
  end;
  set(gcf,'Position',[  32         254        1554         701 ]);  
  saveasshown(gcf,'orthonormal_basis.png');
  saveasshown(gcf,'orthonormal_basis.eps','epsc');
  
 case 8 
  disp('Greedy error convergence');
  
  params = [];
  params.B1 = 2;
  params.B2 = 2;
  params.numintervals_per_block = 20;
  params.mu_range = [0.5;2]; 
  model = thermalblock_model_struct(params);  
  model_data = gen_model_data(model);
  model.RB_stop_epsilon = 1e-6;
  model.RB_stop_Nmax = 50;
  model.RB_train_size = 1000;
  
  disp('Generating basis and determining max error, estimators');
  tic;
  detailed_data = gen_detailed_data(model,model_data);
  reduced_data = gen_reduced_data(model,detailed_data);
  
  disp('Precomputing test snapshots');
  ntest = 100;
  [U,test_mus] = filecache_function(@precompute_test_snapshots,...
                                    params,model,model_data,ntest);
  
  Ns = 1:size(detailed_data.RB,2);
  maxDeltas = zeros(length(Ns),1);
  maxerrs = zeros(length(Ns),1);
  for Ni = 1:length(Ns)
    N = Ns(Ni)
    model.N = N;
    reduced_data_subset = model.reduced_data_subset(model,reduced_data);    
    for i = 1:ntest
      model = set_mu(model,test_mus(i,:));
      rb_sim_data = rb_simulation(model, reduced_data_subset);
      rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
      err = rb_sim_data.uh;
      err.dofs = err.dofs - U(:,i);
      errnorm = fem_h10_norm(err);
      if maxerrs(Ni)<errnorm
        maxerrs(Ni)=errnorm;
      end;
      if maxDeltas(Ni)<rb_sim_data.Delta
        maxDeltas(Ni)=rb_sim_data.Delta;
      end;
    end; % i
  end; % Ni
  time_step_6 = toc  
  save('rb_tutorial_step8');
 
 case 8.1 
  disp('plot of greedy errors');
  load('rb_tutorial_step8');

  plot(Ns,maxDeltas,'b-','Linewidth',2);
  hold on
  plot(Ns,maxerrs,'r:','Linewidth',2);
  set(gca,'Yscale','log');
  % if wanted, training error can be plotted
  %  plot(Ns(1:end-2),detailed_data.RB_info.max_err_est_sequence(2:end),...
%       'g-.','Linewidth',2);
  set(gca,'Yscale','log');
  legend({'test estimator \Delta_u(\mu)','test error |u(\mu)-u_N(\mu)|'},'Fontsize',15);
  title('Error decay for Greedy basis','Fontsize',15);
  xlabel('sample/basis Size N','Fontsize',15);
  ylabel('error/estimator','Fontsize',15);

%  set(gcf,'Position', [379   492   434   328]);
%  saveas(gcf,'greedy_test_error.png');
%  saveas(gcf,'greedy_test_error.eps','epsc');

  % possible further cases: 
  % plot of sample points: not spectacular...
  %  load('rb_tutorial_step8');
  %  mus = detailed_data.RB_info.max_mu_sequence;
  %  mumat = cell2mat(mus)
  %  plot(mumat(1,:),mumat(2,:),'*');

 case 9 
  disp('experiment with basis from 4');
  
  load('rb_tutorial_step4');
  
  mus= {[1,0.2,0.2,0.2],[1,0.1,1,0.1]};
  for i = 1:length(mus)
    disp('-------------------------------------------')
    disp(['Experiment for mu = (',num2str(mus{i}),')']);
    model = set_mu(model,mus{i});
    sim_data = detailed_simulation(model,model_data);
    rb_sim_data = rb_simulation(model,reduced_data);
    rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
    disp(['output s(mu) = ', num2str(sim_data.s)]);
    disp(['output s_N(mu) = ', num2str(rb_sim_data.s)]);
    disp(['output error s(mu)-s_N(mu) = ', ...
          num2str(sim_data.s-rb_sim_data.s)]);
    disp(['output error bound Delta_s(\mu) = ', num2str(rb_sim_data.Delta_s)]);
    err = sim_data.uh;
    err.dofs = err.dofs - rb_sim_data.uh.dofs;
    disp(['field error |u-u_N| = ',num2str(fem_h10_norm(err))]);
    disp(['field error bound Delta_u(\mu) = ', num2str(rb_sim_data.Delta)]);
  end;
  
 case 10 
  disp('convergence curves, training error for different block numbers');

  B1S = [2,3,4];
  B2S = [2,3,4];
  
  max_err_est_sequences = {};
  
  for bi = 1:length(B1S)
    params = [];
    params.B1 = B1S(bi);
    params.B2 = B2S(bi);
    disp('--------------------------------------------------');
    disp(['generating basis for (B1,B2) = (',...
          num2str(params.B1),',',num2str(params.B2),')']);
    params.numintervals_per_block = 20;
    params.mu_range = [0.5;2]; 
    model = thermalblock_model_struct(params);  
    model_data = gen_model_data(model);
    model.RB_stop_epsilon = 1e-6;
    model.RB_stop_Nmax = 50;
    model.RB_train_size = 1000;
    detailed_data = gen_detailed_data(model,model_data);    
    max_err_est_sequences = [max_err_est_sequences,...
                    {detailed_data.RB_info.max_err_est_sequence}]
  end;
  
%  save('rb_tutorial_step10');
% case 10.1 % plot training error curves
%  load('rb_tutorial_step10');
  
  figure;
  hold on;
  legstr = {};
  styles = {'b-','r-','g-','k-'}
  for i=1:length(max_err_est_sequences)
    seq = max_err_est_sequences{i};
    plot(0:length(seq)-1,seq,styles{i});
    legstr = [legstr,{['(B1,B2)= (',...
                       num2str(B1S(i)),',',num2str(B2S(i)),')']}] 
  end;
  legend(legstr,'Fontsize',15);
  title('greedy training error decay','Fontsize',15);
  xlabel('basis size N','Fontsize',15);
  ylabel('training error','Fontsize',15);
  set(gca,'Yscale','log');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%PART 2: Instationary Advection-diffusion example
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 case 11 
  disp('full simulation of advection-diffusion problem.');
  
  params = [];
  params.coarse_factor = 4;
  model = advection_diffusion_fv_model(params);
  model_data = gen_model_data(model);
  plot_params = [];
  plot_params.axis_equal = 1;
  plot_params.axis_tight = 1;
  plot_params.no_lines = 1; 
  
  model = set_mu(model,[1,0.1]);
  tic; sim_data = detailed_simulation(model,model_data);
  disp('time for full simulation (seconds): ')
  t = toc
  plot_sim_data(model,model_data,sim_data,plot_params);
  
 case 12 
  disp('generate plots of varying mu and t:');
  
  params = [];
  params.coarse_factor = 4;
  model = advection_diffusion_fv_model(params);
  model_data = gen_model_data(model);
  plot_params = [];
  plot_params.axis_equal = 1;
  plot_params.axis_tight = 1;
  plot_params.no_lines = 1; 
  plot_params.show_colorbar = 0;   
  plot_params.clim = [0,1];
  
  mus = {[0,0],[1,0],[1,1]};
  ts = [1,129,257];
  snapshots = cell(3,3);
  for mi = 1:length(mus)  
    model = set_mu(model,mus{mi});
    sim_data = detailed_simulation(model,model_data);
    for ti = 1:length(ts);
      snapshots{mi,ti} = sim_data.U(:,ts(ti));
      plot_element_data(model_data.grid,...
                        sim_data.U(:,ts(ti)),plot_params);
      set(gca,'Ytick',[]);
      set(gca,'Xtick',[]);
      saveas(gcf,['advection_diffusion_snapshot-m',num2str(mi),...
                  '-t',num2str(ti),'.png']);
      saveas(gcf,['advection_diffusion_snapshot-m',num2str(mi),...
                  '-t',num2str(ti),'.eps'],'epsc');
    end;
  end;

  % single joint plot
  figure;
  for mi = 1:length(mus)  
    for ti = 1:length(ts);   
      subplot(3,3,(ti-1)*3 + mi);
      plot_element_data(model_data.grid,...
                        snapshots{mi,ti} ,plot_params);
      set(gca,'Ytick',[]);
      set(gca,'Xtick',[]);
      if ti == 1
        mu = mus{mi};
        title(['\mu = (',num2str(mu(1)),',',num2str(mu(2)),')^T'],...
              'Fontsize',15);
      end;
    end;
  end;
  pos = [354         335        1241         613]
  set(gcf,'Position',pos);
  disp('please do a manual => File => Export Setup => Export');
  disp('as files advection_diffusion_snapshots.png and *.eps');
%  saveas(gcf,'advection_diffusion_snapshots.png');
%  saveas(gcf,'advection_diffusion_snapshots.eps','epsc');

 case 13 
  disp('POD-Greedy basis generation, caution takes +- 1 hour!');
  
  params = [];
  params.coarse_factor = 4;
  model = advection_diffusion_fv_model(params);
  model_data = gen_model_data(model);

  % POD-Greedy settings:
  
  % choose initial data components as initial basis:
  model.rb_init_data_basis = @RB_init_data_basis;
  % or choose empty initial basis:
  %model.rb_init_data_basis = @RB_init_basis_empty;
  
  model.RB_generation_mode = 'greedy_uniform_fixed';
  model.RB_extension_algorithm = @RB_extension_PCA_fixspace;
  model.RB_stop_epsilon = 1e-2; 
  model.RB_stop_timeout = 2*60*60; % 2 hours
  model.RB_stop_Nmax = 100;
  % decide, whether estimator or true error is error-indicator for greedy
  % params.RB_error_indicator = 'error'; % true error, i.e. strong POD-greedy
  model.RB_error_indicator = 'estimator'; % Delta from rb_simulation
                                          
  % choose a 10x10 uniform grid
  model.RB_numintervals = 9 * ones(size(model.mu_names));

  detailed_data = gen_detailed_data(model,model_data);
 
  save rb_tutorial_step_13.mat;
  
 case 14
  disp('plot first 16 snapshots and point distribution');
  
  load rb_tutorial_step_13.mat;

  plot_params = [];
  plot_params.axis_equal = 1;
  plot_params.axis_tight = 1;
  plot_params.no_lines = 1; 
  plot_params.show_colorbar = 0;   
  
  plot_detailed_data(model,detailed_data);
  
  disp('inspect diagrams and press enter to continue')
  pause;
    
  close(4); % runtime plot not too interesting
  title('maximal POD-Greedy training estimator','Fontsize',15);
  xlabel('basis vector number N','Fontsize',15);
  ylabel('estimator','Fontsize',15);
  saveas(gcf,'PODgreedy_train_estimator_advection_diffusion.png');
  saveas(gcf,'PODgreedy_train_estimator_advection_diffusion.eps','epsc');
  close(3); % training error decay only partially interesting
  view(2); 
  title('POD-Greedy parameter selection','Fontsize',15);
  xlabel('mu1','Fontsize',15);
  ylabel('mu2','Fontsize',15);
  set(gcf,'Color',[1,1,1]);
%  export_fig('PODgreedy_points_advection_diffusion.png',gcf);
  saveas(gcf,'PODgreedy_points_advection_diffusion.png');
  saveas(gcf,'PODgreedy_points_advection_diffusion.eps','epsc');
  
  figure;
  for n = 1:16
    subplot(4,4,n);
    plot_element_data(detailed_data.grid,detailed_data.RB(:,n), ...
                      plot_params);
    title(['n=',num2str(n)],'Fontsize',15);
    set(gca,'Xtick',[]);
    set(gca,'Ytick',[]);
  end;
  set(gcf,'Position',[32          49        1569         906]);
  saveas(gcf,'basis_advection_diffusion.png');
  saveas(gcf,'basis_advection_diffusion.eps','epsc');
  
%  plot_params = [];
%  plot_params.axis_equal = 1;
%  plot_params.axis_tight = 1;
%  plot_params.no_lines = 1; 

 case 15 
  disp('reduced simulation');
  
  load rb_tutorial_step_13.mat;
  reduced_data = gen_reduced_data(model,detailed_data);
  tic

  model = set_mu(model,[1,0.1]);
  rb_sim_data = rb_simulation(model,reduced_data);
  disp('time for reduced simulation (in seconds): ')
  t = toc
  rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);

  plot_params = [];
  plot_params.axis_equal = 1;
  plot_params.axis_tight = 1;
  plot_params.no_lines = 1; 
  plot_params.show_colorbar = 0;   

  plot_sim_data(model,model_data,rb_sim_data,plot_params);
  disp('press enter for rb simulation gui')
  pause;
  
  demo_rb_gui(model,detailed_data,[],plot_params,'rb simulation');

  disp('press enter for rb error simulation gui')
  pause;

  plot_params.show_colorbar = 1;   
  demo_rb_error_gui(model,detailed_data,[],plot_params,'rb error simulation');

 case 16 
  disp('error, estimator by parameter sweep');
  
  % sweep along the diagonal of the parameter space, i.e. 
  % some training points, but ´mostly test-points
  
  load rb_tutorial_step_13.mat;
  reduced_data = gen_reduced_data(model,detailed_data);
  
  ntest = 101; % higher later, take care that 10x10 training grid is contained
              
  %  RandStream.setDefaultStream(RandStream('mt19937ar','seed', ...
  %                                      12345));
  %  mu_test = rand_uniform(ntest,model.mu_ranges);
  
  mu_test = [1;1]*(0:(ntest-1))/(ntest-1)
  
  errs_max_t = zeros(ntest,1);
  errs_final_t = zeros(ntest,1);
  Deltas_final_t = zeros(ntest,1);
  
  W = model.get_inner_product_matrix(detailed_data);
  for mi = 1:ntest
    disp(['processing parameter ',num2str(mi),'/',num2str(ntest)]);
    model = model.set_mu(model,mu_test(:,mi));
    sim_data = detailed_simulation(model,model_data);
    rb_sim_data = rb_simulation(model,reduced_data);
    rb_sim_data = rb_reconstruction(model,detailed_data, ...
                                    rb_sim_data);
    errs = rb_sim_data.U(:,model.save_time_indices+1) -sim_data.U;
    err_norms = sqrt(sum(errs.* (W * errs)));
    errs_final_t(mi) = err_norms(end);
    errm = max(err_norms);
    errs_max_t(mi) = errm(1);
    Deltas_final_t(mi) = rb_sim_data.Delta(end);
  end;
  p = plot(mu_test(2,:),[errs_max_t(:),Deltas_final_t(:)]);
  set(p,'Linewidth',2);
  set(p(2),'Linestyle','-.');
  title('error and error estimator at final time t=T','Fontsize',15); 
  legend({'error','estimator'},'Fontsize',15);
  set(gca,'Yscale','log');
  ylabel('error/estimator','Fontsize',15);
  xlabel('parameter scaling factor s','Fontsize',15);
  
  saveas(gcf,'advection-diffusion-error-parameter-sweep.png');
  saveas(gcf,'advection-diffusion-error-parameter-sweep.eps','epsc');
  saveas(gcf,'advection-diffusion-error-parameter-sweep.fig');
  
 case 17 
  disp('Effectivity computation for 200 random points');

  % plot: test effectivity over diffusivity parameter
  
  load rb_tutorial_step_13.mat;
  reduced_data = gen_reduced_data(model,detailed_data);

  if matlab_version> 7.9
    s = RandStream('mt19937ar','Seed',12345)
    RandStream.setGlobalStream(s);
  else
    RandStream.setDefaultStream(RandStream('mt19937ar','seed', ...
                                 12345));
  end;
  
  ntest = 200; % higher later
  mu_test = rand_uniform(ntest,model.mu_ranges);
  errs_max_t = zeros(ntest,1);
  errs_final_t = zeros(ntest,1);
  Deltas_final_t = zeros(ntest,1);

  for mi = 1:ntest
    disp(['processing parameter ',num2str(mi),'/',num2str(ntest)]);
    model = model.set_mu(model,mu_test(:,mi));
    sim_data = detailed_simulation(model,model_data);
    rb_sim_data = rb_simulation(model,reduced_data);
    rb_sim_data = rb_reconstruction(model,detailed_data, ...
                                    rb_sim_data);
    errs = rb_sim_data.U(:,model.save_time_indices+1) -sim_data.U;
    err_norms = sqrt(sum(errs.* (model.get_inner_product_matrix(detailed_data) * errs)));
    errs_final_t(mi) = err_norms(end);
    errm = max(err_norms);
    errs_max_t(mi) = errm(1);
    Deltas_final_t(mi) = rb_sim_data.Delta(end);
  end;
%  keyboard;
%  [mu_sort,ind ] = sort(mu_test(2,:)); 
%  plot(mu_test(2,ind),[Deltas_final_t(ind)./errs_max_t(ind)]);

  [mu_sort,ind ] = sort(mu_test(1,:)); 
  plot(mu_test(1,ind),[Deltas_final_t(ind)./errs_max_t(ind)],'b*');
  ylim = get(gca,'Ylim');
  ylim(1) = 0;
  set(gca,'Ylim',ylim);
  
  % plot: maximum test error over Basis size
  title('effectivities at final time t=T over test set','Fontsize',15); 
  ylabel('effectivity','Fontsize',15);
  xlabel('velocity weight mu1','Fontsize',15);
  
  saveas(gcf,'advection-diffusion-effectivity.png');
  saveas(gcf,'advection-diffusion-effectivity.eps','epsc');
  saveas(gcf,'advection-diffusion-effectivity.fig');
   
 otherwise
  disp('step number unknown');
end;


function mu_list = lin_equidist_samples(model,numintervals)
% generate cell array of parameter vectors 
% numintervals a vector of number of equal sized intervals per 
% parameter direction 
p = [];
p.range = model.mu_ranges;
p.numintervals = numintervals;
c = cubegrid(p);
mu_list = mat2cell(c.vertex,ones(size(c.vertex,1),1),...
                   size(c.vertex,2));  

function [U,test_mus] = precompute_test_snapshots(params,model,model_data,ntest)
% function precomputing test snapshots for step 6
test_mu1s = rand_uniform(ntest,{params.mu_range});
test_mus = [test_mu1s(:), ones(ntest,params.B1*params.B2-1)*1];    
U = zeros(model_data.df_info.ndofs,0);
for i = 1:ntest
  fprintf('.');
  model = set_mu(model,test_mus(i,:));
  sim_data = detailed_simulation(model, model_data);
  U = [U,sim_data.uh.dofs(:)]; 
end;
fprintf('\n');

function RB = lagrangian_orth_basisgen(model, ...
                                                  detailed_data)
Utot = [];
for m = 1:length(model.RB_mu_list)
  model = set_mu(model,model.RB_mu_list{m});
  sim_data = detailed_simulation(model,detailed_data);
  if isempty(Utot)
    Utot = model.get_dofs_from_sim_data(sim_data);
  else
    U = model.get_dofs_from_sim_data(sim_data);
    Utot = [Utot U];
  end;
end;
% Gram-Schmidt orthonormalization
%KN = Utot'*model.get_inner_product_matrix(detailed_data)*Utot; 
%KN = 0.5*(KN + KN');
% This does no longer work with recent matlab:
%Ltransp = chol(KN);
%RB = Utot / Ltransp;
RB = orthonormalize_qr(Utot,model.get_inner_product_matrix(detailed_data),1e-12);