rb_tutorial_standalone.m

function res = rb_tutorial_standalone(step)
%function res = rb_tutorial_standalone(step)
% 
% step : index of experiment to be performed
%        1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 100
%
% This program allows to reproduce the thermalblock experiments of
% the RB-Tutorial use in the 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
%
% This is a modified version of the original rb_tutorial.m in
% RBmatlab, which is available for download at www.morepas.org.
% The modifications are such that this experiments can be run
% standalone, i.e. without the RBmatlab classes. But only the stationary
% experiments have been transferred, i.e. the instationary experiments can
% be found in rb_tutorial.m
% For making this file independent of RBmatlab, 
% the complete FEM-discretization has been precomputed in
% a file data_rb_tutorial_standalone.mat, which is used here. 
% (This data can be (re-)generated in step 100 if RBmatlab is installed.)

% B. Haasdonk 10.4.2015

if nargin < 1
  help rb_tutorial_standalone
  disp('Now running step 1:')
  step = 1
end;

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 = standalone_thermalblock_model(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 = standalone_thermalblock_model(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('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 = standalone_thermalblock_model(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 + error estimator components, reduced basis.');
  detailed_data = gen_detailed_data(model,model_data);
  disp('generating reduced data: A_Nq,f_Nq,l_Nq + error estimator 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
  rb_sim_data = rb_simulation(model,reduced_data);
  t_online = toc;

  disp('Online time: ')
  t_online
      
 case 3 
  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 = standalone_thermalblock_model(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 - rb_sim_data.uh;
    errs(n) = sqrt(err'*detailed_data.K*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_standalone_step3');
  
 case 4  
  disp('Effectivity plot');
  load('rb_tutorial_standalone_step3');
  
  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);
  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 5 
    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 = standalone_thermalblock_model(params);  
    model_data = gen_model_data(model);
%    model.RB_generation_mode = 'greedy_rand_uniform_with_orthogonalization';
    model.RB_generation_mode = 'lagrangian_with_orthonormalization';
    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] = 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 = train_mus;
      %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 - U(:,i);
        errnorm = sqrt(err' * detailed_data.K * 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_5 = toc  
    save('rb_tutorial_standalone_step5');
    % case 6.1 % plot of error convergenc 
    %  load('rb_tutorial_standalone_step5');
    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 6 
  disp('plot basis');
  
  load('rb_tutorial_standalone_step5');
  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';
    
  % or in single plot
  plot_params.show_colorbar = 0;
  plot_params.no_lines = 1;
  figure;
  for n = 1:8
    subplot(2,4,n);
    phi_ext = zeros(length(detailed_data.grid.X),1);
    % expand vector to full size including dirichlet values
    phi_ext(detailed_data.grid.non_dirichlet_indices) = detailed_data.RB(:,n);
    plot_vertex_data(detailed_data.grid,phi_ext, ...
                     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 7   
  disp('Greedy error convergence');
  
  params = [];
  params.B1 = 2;
  params.B2 = 2;
  params.numintervals_per_block = 20;
  params.mu_range = [0.5;2]; 
  model = standalone_thermalblock_model(params);  
  model_data = gen_model_data(model);
  model.RB_generation_mode = 'greedy_rand_uniform_with_orthogonalization';
  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);
  [U,test_mus] = 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 = 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 - U(:,i);
      errnorm = sqrt(err'*detailed_data.K*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_standalone_step7');
 
 case 8 
  disp('plot of greedy errors');
  load('rb_tutorial_standalone_step7');

  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_standalone_step8');
  %  mus = detailed_data.RB_info.max_mu_sequence;
  %  mumat = cell2mat(mus)
  %  plot(mumat(1,:),mumat(2,:),'*');

 case 9 
  disp('experiment with basis from step 3');
  
  load('rb_tutorial_standalone_step3');
  
  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 - rb_sim_data.uh;
    errnorm = sqrt(err'*detailed_data.K*err);
    disp(['field error |u-u_N| = ',num2str(errnorm)]);
    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 = standalone_thermalblock_model(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_standalone_step10');
% case 10.1 % plot training error curves
%  load('rb_tutorial_standalone_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');

 case 100 
  disp('generate data file data_rb_tutorial_standalone.mat (with RBmatlab!)');

  if ~exist('rbmatlabhome','file')
    error('This step can only be executed with RBmatlab installed')
    error(['This step is not required, if the data_rb_tutorial_standalone.mat file is', ...
          'provided.']);
  end;
  
  save_varlist = {};  
  params.mu_range = [0.1;10];
  B1_list = [2,6,2,3,4];
  B2_list = [2,6,2,3,4];
  numintervals_per_block_list = [5,5,20,20,20];
  for i = 1:length(B1_list);
    params.B1 = B1_list(i);
    params.B2 = B2_list(i);
    params.numintervals_per_block = numintervals_per_block_list(i);
    model = thermalblock_model_struct(params);
    %    model_data = lin_stat_gen_model_data(model);
    model_data = [];
    md.grid = construct_grid(model);
    md.df_info = feminfo(model,md.grid);
    model.decomp_mode = 1;
    [A_comp,f_comp] = model.operators(model,md);
    % only Neumann in our case:
    f = f_comp{3}; 
    % only Diffusion components:
    A_comp = A_comp(1:(model.B1*model.B2)); 
    l_comp = model.operators_output(model,md);
    % cut down matrices to inner dofs:
    dirichlet_indices = md.df_info.dirichlet_gids;
    dirichlet_flag = zeros(md.grid.nvertices,1);
    dirichlet_flag(dirichlet_indices)=1;
    non_dirichlet_indices = find(dirichlet_flag==0);
    for i = 1:length(A_comp);
      A_comp_i = A_comp{i};
      A_comp_i = A_comp_i(non_dirichlet_indices,non_dirichlet_indices);  
      A_comp{i} = A_comp_i;
    end;
    f = f(non_dirichlet_indices);  
    l = l_comp{1};
    l = l(non_dirichlet_indices);  
    % export only structures, not classes:
    model_data.grid.X = md.grid.X;
    model_data.grid.Y = md.grid.Y;
    model_data.grid.VI = md.grid.VI;
    model_data.grid.dirichlet_indices = md.df_info.dirichlet_gids;
    model_data.grid.non_dirichlet_indices = non_dirichlet_indices;
    model_data.A_comp = A_comp;
    model_data.f = f;
    model_data.l = l;
    model_data.K = md.df_info.regularized_h10_inner_product_matrix(...
        non_dirichlet_indices,non_dirichlet_indices);

    model_data_varname = ['model_data_',num2str(params.B1),'_',num2str(params.B2),'_',...
                num2str(params.numintervals_per_block)];
    
    eval([model_data_varname, ' = model_data;']);
    save_varlist = [save_varlist, {model_data_varname}];
  end;  
  save('data_rb_tutorial_standalone.mat',save_varlist{:});  
  
 otherwise
  error('step unknown!');
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% auxiliary functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% model generation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function model = standalone_thermalblock_model(params)
% Thermal Block example 
%
% i.e. 
% - div ( a(x) grad u(x)) = q(x)    on Omega
%                   u(x)) = g_D(x)  on Gamma_D
%     a(x) (grad u(x)) n) = g_N(x)  on Gamma_N
%                       s = l(u) linear output functional
%
% where Omega = [0,1]^2 is divided into B1 * B2 blocks 
% QA := B1*B2. The heat conductivities are given by mu_i:
%
%   -------------------------
%   | ...   | ..... | mu_QA |
%   -------------------------
%   | ....  | ...   | ..... |
%   -------------------------
%   | mu_1  | ...   | mu_B1 |
%   -------------------------
%
% `Gamma_D =` upper edge
% `Gamma_N = boundary(Omega)\Gamma_D`
%
% `a(x) = mu_i(x) if x\in B_i`
% `q(x) = 0`
% `g_D  = 0` on top
% `g_N  = 1` on lower edge, 0 otherwise
% `l(u)` = average over lower edge
%
% The discretization information is fully cut out here and loaded
% from disc in gen_model_data

model = [];
if nargin == 0
    params.B1=3;
    params.B2=3;
end
model.B1 = params.B1;
model.B2 = params.B2;
model.number_of_blocks = params.B1*params.B2;

if ~isfield(params,'numintervals_per_block')
  params.numintervals_per_block = 5;
end;
model.numintervals_per_block = params.numintervals_per_block;

mu_names = {};
mu_ranges = {};
if ~isfield(params,'mu_range')
  params.mu_range = [0.1,10];
end;
mu_range = params.mu_range;
for p = 1:model.number_of_blocks
  mu_names = [mu_names,{['mu',num2str(p)]}];
  mu_ranges = [mu_ranges,{mu_range}];
end;
model.mu_names = mu_names;
model.mu_ranges = mu_ranges;

%default values 1 everywhere
model.mus = ones(model.number_of_blocks,1);

% basis generation settings
%model.RB_train_rand_seed = 100;
model.RB_train_size = 1000;
model.RB_stop_epsilon = 1e-3;
model.RB_stop_Nmax = 100;
model.RB_generation_mode = 'greedy_rand_uniform_with_orthogonalization';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%% high level functions %%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function model = set_mu(model,mu)
if length(mu)~=model.number_of_blocks
  error('length of mu does not fit to number of blocks!');
end;
model.mus = [mu(:)];


function model_data = gen_model_data(model)
% simply load precomputed matrices and grid information
model_data_varname = ['model_data_',num2str(model.B1),'_',num2str(model.B2),'_',...
                    num2str(model.numintervals_per_block)];
load('data_rb_tutorial_standalone',model_data_varname);
eval(['model_data =', model_data_varname,';']);


function detailed_data = gen_detailed_data(model,model_data)
detailed_data = model_data;
detailed_data = rb_basis_generation(model,detailed_data);


function detailed_data = rb_basis_generation(model,detailed_data)
switch model.RB_generation_mode 
 
 case 'lagrangian'   
   Utot = [];
   for m = 1:size(model.RB_mu_list,1)
     model = set_mu(model,model.RB_mu_list(m,:));
     sim_data = detailed_simulation(model,detailed_data);
     if isempty(Utot)
       Utot = sim_data.uh;
     else
       Utot = [Utot sim_data.uh];
     end;
   end;
   detailed_data.RB = Utot;
   
 case 'lagrangian_with_orthonormalization'   
  
  Utot = [];
  for m = 1:size(model.RB_mu_list,1)
    model = set_mu(model,model.RB_mu_list(m,:));
    sim_data = detailed_simulation(model,detailed_data);
    if isempty(Utot)
      Utot = sim_data.uh;
    else
      Utot = [Utot, sim_data.uh];
    end;
  end;
  % Gram-Schmidt orthonormalization
  KN = Utot'* detailed_data.K *Utot; 
  KN = 0.5*(KN + KN');

  % in case of singular KN produce less RB vectors: 
  % Ltransp = chol(KN);   
  % RB = Utot / Ltransp;   
  
  %[Ltransp,p] = chol(KN);
  %if p==0
  %  r = size(Utot,2);
  %else
  %  r = p-1;
  %end;
  %RB = Utot(:,1:r) / Ltransp;
  %detailed_data.RB = RB;
  
  RB = orthonormalize_qr(Utot,detailed_data.K,1e-12);
  detailed_data.RB = RB;
  
 case 'greedy_rand_uniform_with_orthogonalization'

  sim_data = detailed_simulation(model, detailed_data);
  n = sqrt(sim_data.uh'*detailed_data.K*sim_data.uh);
  detailed_data.RB = sim_data.uh / n;
  reduced_data = gen_reduced_data(model, detailed_data);
  
  mus = rand_uniform(model.RB_train_size,model.mu_ranges);
  max_err_est = 1e10;
  max_err_est_sequence = [];
  max_mu_sequence = {};
  
  while( max_err_est> model.RB_stop_epsilon) && ...
        (size(detailed_data.RB,2)< model.RB_stop_Nmax)
    max_err_est = 0;
    mu_max = [];
    for i = 1:size(mus,2);
      model = set_mu(model, mus(:,i));
      rb_sim_data = rb_simulation(model,reduced_data);
      %    disp(['rb sim: Delta = ',num2str(rb_sim_data.Delta)]);
      if rb_sim_data.Delta > max_err_est
        mu_max = mus(:,i);
        max_err_est = rb_sim_data.Delta; 
      end;
    end;
    disp(['max error estimator: ',num2str(max_err_est),...
          ' for mu = ',num2str(mu_max')]);
    max_err_est_sequence = [max_err_est_sequence,max_err_est];
    max_mu_sequence = [max_mu_sequence,{mu_max}];
    model = set_mu(model,mu_max);
    sim_data = detailed_simulation(model,detailed_data);
    pr = sim_data.uh - ...
         detailed_data.RB * (detailed_data.RB' * detailed_data.K * ...
                             sim_data.uh);
    n = sqrt(max(pr'*detailed_data.K*pr,0));
    detailed_data.RB = [detailed_data.RB,pr/n];  
    reduced_data = gen_reduced_data(model,detailed_data);
    disp(['new basis size: ',num2str(size(detailed_data.RB,2))]);
  end;
  
  RB_info.max_mu_sequence = max_mu_sequence;
  RB_info.max_err_est_sequence = max_err_est_sequence;
  detailed_data.RB_info = RB_info;
 
 otherwise
  error('RB_generation_mode unknown!!!');

end;


function sim_data = detailed_simulation(model,model_data)
sim_data =[];
A_coeff = model.mus;
Q_A = length(A_coeff);
A = model_data.A_comp{1} * A_coeff(1);
for q=2:Q_A
  if A_coeff(q)~=0
    A = A + A_coeff(q)*model_data.A_comp{q};
  end;
end;
uh = A\model_data.f;    
sim_data.uh = uh;
sim_data.s = (model_data.l(:)') * sim_data.uh;


function p = plot_sim_data(model,detailed_data,sim_data, ...
                                    plot_params)
if nargin<4
  plot_params = [];
end;

% expand vector to full size including dirichlet values
uh_ext = zeros(length(detailed_data.grid.X),1);
uh_ext(detailed_data.grid.non_dirichlet_indices) = sim_data.uh(:);

p = plot_vertex_data(detailed_data.grid,uh_ext,plot_params);
hold on;

% plot coarse mesh
if ~isfield(plot_params,'plot_blocks')
  plot_params.plot_blocks = 1;
end;
if plot_params.plot_blocks
  X = [0:1/model.B1:1];
  Y = [0:1/model.B2:1];
  l1 = line([X;X],...
            [zeros(1,model.B1+1);...
             ones(1,model.B1+1)]);
  set(l1,'color',[0,0,0],'linestyle','-.');
  %keyboard;
  l2 = line([zeros(1,model.B2+1);...
             ones(1,model.B2+1)],...
            [Y;Y]);
  set(l2,'color',[0,0,0],'linestyle','-.');
  p = [p(:);l1(:);l2(:)];
end;


function p = plot_vertex_data(grid,data,plot_params);
if nargin<3
  plot_params = [];
end;
if ~isfield(plot_params,'title');
  plot_params.title = '';
end;
if ~isfield(plot_params,'axis_equal');
  plot_params.axis_equal = 0;
end;
if ~isfield(plot_params,'no_lines');
  plot_params.no_lines = 1;
end;
if ~isfield(plot_params,'show_colorbar');
  plot_params.show_colorbar = 1;
end;
if ~isfield(plot_params,'colorbar_location');
  plot_params.colorbar_location = 'EastOutside';
end;
% expand vector to full size including dirichlet values
XX = grid.X(grid.VI(:));        
XX = reshape(XX,size(grid.VI)); % nelements*nneigh matrix
YY = grid.Y(grid.VI(:));       
YY = reshape(YY,size(grid.VI)); % nelements*nneigh matrix
CC = data(grid.VI);
p = patch(XX',YY',0*CC', CC');
c = jet(256); 
colormap(c);
if plot_params.axis_equal
  axis equal;
  axis tight;
end;
if plot_params.no_lines
  set(p,'linestyle','none');
end;
if plot_params.show_colorbar
  if isfield(plot_params,'clim')
    set(gca,'Clim',plot_params.clim)
  end;
  colorbar(plot_params.colorbar_location);  
end;
title(plot_params.title);


function reduced_data = gen_reduced_data(model,detailed_data)
reduced_data = [];
N = size(detailed_data.RB,2);
reduced_data.N = N;
A_comp = detailed_data.A_comp;
f = detailed_data.f;
Q_A = length(A_comp);
reduced_data.Q_A = Q_A;
reduced_data.AN_comp = zeros(N^2,Q_A);
for q = 1:Q_A
  AN_comp_q = detailed_data.RB'*A_comp{q}*detailed_data.RB;
  reduced_data.AN_comp(:,q) = AN_comp_q(:); 
end;
% f,l non parametric, hence simple projection
Q_f = 1;
reduced_data.fN = detailed_data.RB' * detailed_data.f;
reduced_data.lN = detailed_data.RB' * detailed_data.l; 

% plus error estimation quantities
% G = (v_r^q, v_r^q) = v_r^q' * K * v_r^q
% with {v_r^q}_q = (v_f^q, v_a^q'n)_{q'n}
% G = [Gff, Gfa; Gfa', Gaa];
%
% (v_f^q,v_f^q') = v_f^q' * K * v_f^q
% matrices of coefficient vectors of Riesz-representers:
% K * v_f^q = f^q      (coefficient vector equations)
K = detailed_data.K;
K_v_f = detailed_data.f;
v_f = K \ K_v_f;
Q_r = Q_f + N * Q_A;
K_v_r = zeros(length(detailed_data.f),Q_r);
K_v_r = [detailed_data.f];
for q = 1:Q_A
  K_v_r = [K_v_r, detailed_data.A_comp{q}*detailed_data.RB];
end;
v_r = K \ K_v_r;
G = v_r' * K_v_r; 
reduced_data.G = G;


function reduced_data_ss = reduced_data_subset(model,reduced_data)
reduced_data_ss = reduced_data;

if reduced_data.N ~= model.N
  % extract correct N-sized submatrices from reduced_data
  if model.N > reduced_data.N
    error('N too large for current size of reduced basis!');
  end;
  N = model.N;
  dummy = zeros(reduced_data.N);
  dummy(1:N,1:N)= 1;
  ind = find(dummy);
  reduced_data_ss.AN_comp   =  reduced_data.AN_comp(ind,:);
  reduced_data_ss.fN = reduced_data.fN(1:N);  
  reduced_data_ss.lN = reduced_data.lN(1:N);  
  % inner product matrix of riesz-representers of residual components:
  Q_f = 1;
  Q_a = size(reduced_data_ss.AN_comp,2);

  ind = reshape(1:reduced_data.N*Q_a,reduced_data.N,Q_a)'; 
  ind = ind(:,1:N)';
  ind = [1; ind(:) + 1];  
  reduced_data_ss.G = reduced_data.G(ind,ind);
  
  reduced_data_ss.N = N;

end


function rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data)
rb_sim_data.uh = ...
    detailed_data.RB(:,1:length(rb_sim_data.uN)) * rb_sim_data.uN;


function rb_sim_data = rb_simulation(model,reduced_data)
rb_sim_data =[];
A_coeff = model.mus(:);
N = reduced_data.N;
AN = reshape(reduced_data.AN_comp * A_coeff, N , N);
uN = AN\reduced_data.fN;
rb_sim_data.uN = uN;
rb_sim_data.s = reduced_data.lN(:)' * rb_sim_data.uN;
% for elliptic compliant case, X-norm (=H10-norm) error estimator:
Q_r = size(reduced_data.G,1);
neg_auN_coeff = -uN * A_coeff'; 
f_coeff = 1;
res_coeff = [f_coeff; neg_auN_coeff(:)];
res_norm_sqr = res_coeff' * reduced_data.G * res_coeff;
% prevent possibly negative numerical disturbances:
res_norm_sqr = max(res_norm_sqr,0);    
res_norm = sqrt(res_norm_sqr);
alpha_LB = min(model.mus);
rb_sim_data.Delta = ...
    res_norm/alpha_LB; 
rb_sim_data.Delta_s = ...
            res_norm_sqr/alpha_LB; 


function mu_list = lin_equidist_samples(model,numintervals)
mu_list = [];
for i = 1:length(model.mus);
  mu_range = model.mu_ranges{i};
  n = numintervals(i)+1;
  mi = linspace(mu_range(1),mu_range(2),n);
  if i==1
    mmi = mi;
  else
    mmi = ones(size(mu_list,1),1) * mi(:)';
  end;
  mmi_vec = mmi(:);
  mu_list = [repmat(mu_list,n,1), mmi_vec];
end;


function [U,test_mus] = precompute_test_snapshots(params,model,model_data,ntest)
% function precomputing test snapshots for step 5
test_mu1s = rand_uniform(ntest,{params.mu_range});
test_mus = [test_mu1s(:), ones(ntest,params.B1*params.B2-1)*1];    
n = length(model_data.grid.non_dirichlet_indices);
U = zeros(n,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(:)]; 
end;
fprintf('\n');

function Xon = orthonormalize_qr(X,A,epsilon)
% function Xon = orthonormalize_qr(X[,A,epsilon])
%
% orthonormalization of vectors in columns of 
% matrix X to Xon by qr-decomposition of vectors and 
% the inner product matrix A can be specified additionally, i.e.
% <x,x> = x' * A * x. Then the Cholesky factoriziation of A is 
% performed A = R_M' * R_M
%
% Then a QR-decomposition of R_M * X is performed which yields the
% desired new vectors
if nargin < 3
  epsilon = 1e-7; % => for nonlinear evolution this value is required.
end;
if isempty(X)
  Xon = zeros(size(X));
  return;
end;
if nargin<2
  A = 1;
end;
% check on identity of vectors (can happen, that numerics does not detect
% this afterwards due to rounding errors!!)
for i=1:(size(X,2)-1)
  for j=(i+1):size(X,2)
    if isequal(X(:,i),X(:,j))
      X(:,j) = 0;
    end;
  end;
end;
Xon = X;
% incomplete cholesky of inner-product matrix:
if matlab_version > 7.9
  opts.shape = 'upper';
    opts.type = 'ict';
  R_M = ichol(sparse(A), opts);
else
     R_M = cholinc(sparse(A),'inf');
end;
% qr decomposition of R_M * X with permutation indices E
[Q, R, E]= qr(R_M * Xon,0 );
% sort such that linear independent first and Q * R = R_M * Xon
Xon = Xon(:,E);
% search nonvanishing diagonal entries of R
ind = find(abs(diag(R))>epsilon);
Xon = Xon(:,ind);
Rind = R(ind,ind); 
Xon = Xon(:,ind) / Rind;
% resort vectors such that original order is recovered most:
K = abs(X'* A * Xon);
% search in each row maximum entry
ind = [];
for i = 1:size(Xon,2)
  [max_val, max_i] = max(K(i,:)); 
  if length(max_i)>1
    error('found multiple matching vectors!');
  end;
  if max_val > 0
    ind = [ind, max_i];
  end;
  K(:,max_i) = 0; % "delete" from further search
end;
all_ind = ones(1,size(Xon,2));
all_ind(ind) = 0;
rest_ind = find(all_ind);
% resort:
if ~isempty(rest_ind)
  Xon = [Xon(:,ind),Xon(:,rest_ind)];
else
  Xon = Xon(:,ind);  
end;

function M = rand_uniform(N,intervals)
%function M = rand_uniform(N,intervals)
%
% function generating uniformly distributed random data in a hypercube
% N vectors are generated as columns of the matrix M, intervals is a 
% cell array indicating the borders of the intervals
%  
% example: rand_uniform(100,{[0,1],[100,1000]})
  mu_dim = length(intervals);
  M = rand(mu_dim,N);
  for i=1:mu_dim
    mu_range = intervals{i};
    M(i,:) = M(i,:)*(mu_range(2)-mu_range(1))+mu_range(1);
  end;  

function res = matlab_version
%function res = matlab_version                                                                                                                                      
%                                                                                                                                                                   
% return numerical value of MATLAB version. Currently used to                                                                                                      
% switch between ichol and cholinc.                                                                                                                                                                     
v = ver('MATLAB');
res = str2num(v.Version);