rb_tutorial.m

function res = rb_tutorial(step)

% script for RB Summerschool in Munich, 16.9.2013-19.9.2013
% script used for generating pictures for course

% B. Haasdonk, 7.9.2013

if nargin<1
  step = 1
end;

switch step
 case 1 % 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 % 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 % 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 % 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  % 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(etas>100);
  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 % 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 % 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',15);
    set(gca,'Xtick',[]);
    set(gca,'Ytick',[]);
  end;
  set(gcf,'Position',[32          49        1569         906]);
%  saveas(gcf,'orthonormal_basis.png');
%  saveas(gcf,'orthonormal_basis.eps','epsc');
  
 case 8 % 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 % 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 % 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 % 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');

 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');
Ltransp = chol(KN);
RB = Utot / Ltransp;