vi_rb.m

function res = vi_rb(step)
%
% experiments for SINUM ARticle 2012, VI-RB paper from Julien script
% this means, that this just reprograms the results, while results
% of paper have been obtained by non-RBmatlab code.

% Bernard Haasdonk

if nargin == 0 || isempty(step)
  step = 10; % == demo_vi
end;

%model = elastic_rope_model_old;
model = elastic_rope_model;
model_data = gen_model_data(model);

switch step
 case 0 % check detailed discretization by grid refinement
  
  %  model = set_mu(model,[21.7157/200, 0.1111])
  model = set_mu(model,[21.7157/200, 0.1111])
  % detailed simulation and plot
  disp('begin detailed simulation');
  sim_data = detailed_simulation(model,model_data);
  disp('plotting ...');
  plot_params.plot_title = 'detailed solution & obstacle';
  p = plot_sim_data(model,model_data,sim_data,plot_params);
  title('xnumintervals=200');
  
  params.xnumintervals = 400;
  model = elastic_rope_model(params);
  model_data = gen_model_data(model);
  model = set_mu(model,[21.7157/200, 0.1111])
  % detailed simulation and plot
  disp('begin detailed simulation');
  sim_data = detailed_simulation(model,model_data);
  disp('plotting ...');
  plot_params.plot_title = 'detailed solution & obstacle';
  figure;
  p = plot_sim_data(model,model_data,sim_data,plot_params);
  title('xnumintervals = 400');
  keyboard;
  
 case 1 % demonstration of model methods

%  model = set_mu(model,[21.7157/200, 0.1111])
%  model = set_mu(model,[21.7157/200, 0.1111])
  model = set_mu(model,[0.075, 0.4])
  % detailed simulation and plot
  disp('begin detailed simulation');
  sim_data = detailed_simulation(model,model_data);
  disp('plotting ...');
  plot_params.plot_title = 'detailed solution & obstacle';
  plot_params.solution_component = 'lambda';
  p = plot_sim_data(model,model_data,sim_data,plot_params);
    
  disp('press enter to continue');
  pause
  
  % basis generation, reduced simulation, error comp.

  % rb generation
  disp('compute snapshots');
  detailed_data = gen_detailed_data(model,model_data);
  %save('detailed_data_vi.mat','model','detailed_data');
  
  reduced_data = gen_reduced_data(model,detailed_data);
  
  % rb simulation and plot
  mus = {[30/200,0],detailed_data.RB_info.mu_train(1,:)};
  %  model = model.set_mu(model,[30,0]); % outside training-sample
  
  for mui = 1:2;
    if mui==1
      disp('-------------------------------------')
      disp('simulation for non-training parameter');
    else
      disp('-------------------------------------')
      disp('simulation for training parameter');
    end;
    model = model.set_mu(model,mus{mui});
    
    sim_data = detailed_simulation(model,model_data);
    disp('plotting ...');
    plot_params.plot_title = 'detailed solution & obstacle';
    figure;
    subplot(1,2,1);
    p = plot_sim_data(model,model_data,sim_data,plot_params);
    
    disp('begin reduced simulation');
    rb_sim_data = rb_simulation(model,reduced_data);
    rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
    disp('plotting ...');
    plot_params.plot_title = 'reduced solution & obstacle';
    subplot(1,2,2);
    p = plot_sim_data(model,model_data,rb_sim_data,plot_params);
    
    % error computation
    disp('compute errors and norms');
    K = model.get_inner_product_matrix(model,model_data);
    U_err = sqrt((rb_sim_data.U-sim_data.U)' * K * (rb_sim_data.U- ...
                                                    sim_data.U));
    L_R = K \ (rb_sim_data.L-sim_data.L);
    L_err = sqrt((rb_sim_data.L-sim_data.L)' * L_R);
    disp(['error (U-UN): ',num2str(U_err)]);
    disp(['error (L-LN): ',num2str(L_err)]);
    
    U_norm = sqrt(sim_data.U' * K * sim_data.U);
    L_R = K \ sim_data.L;
    L_norm = sqrt(sim_data.L' * L_R);
    disp(['norm U: ',num2str(U_norm)]);
    disp(['norm L: ',num2str(L_norm)]);
    
    % error-estimator:
    model.enable_error_estimator = 1;
    reduced_data = gen_reduced_data(model,detailed_data);
    rb_sim_data = rb_simulation(model,reduced_data);
    disp(['error estimator (U-UN) <= Delta_u = ',num2str(rb_sim_data.Delta_u)]);
    disp(['error estimator (L-LN) <= Delta_l = ',num2str(rb_sim_data.Delta_lambda)]);
    
  end;
  keyboard;

 case 1.5 % check accuracy of rb-scheme for large Lagrange basis: 289 vectors
  % IMPORTANT STEP HIGHLIGHTING ACCURACY PROBLEM
  
  model.RB_numintervals = [16,16];
  model.RB_generation_mode = 'pure_snapshots'; % no supr.
  model.enable_error_estimator = 1;
%  model.detailed_simulation = @my_new_detailed_simulation;
  disp('generating Lagrange basis');
  detailed_data = filecache_function(@gen_detailed_data,model,model_data);
  save('lagrange-basis-detailed-data','model','detailed_data');

  reduced_data = gen_reduced_data(model,detailed_data);

  if 0
  disp('searching maximum error estimator over train set')
  max_Delta_u = 0;
  max_Delta_lambda = 0;
  for m = 1:289
    model = set_mu(model,detailed_data.RB_info.mu_train(m,:));
    rb_sim_data = rb_simulation(model,reduced_data);
    if max_Delta_u <     rb_sim_data.Delta_u
      max_Delta_u = rb_sim_data.Delta_u
      disp(['found for m = ',num2str(m)]);
    end;
    if max_Delta_lambda <     rb_sim_data.Delta_lambda
      max_Delta_lambda = rb_sim_data.Delta_lambda
      disp(['found for m = ',num2str(m)]);
    end;
  end;
  disp('maximum error estimators:')
  max_Delta_u
  max_Delta_lambda
  %%%%%%%% RESULT:
  %%%% max_Delta_u = 0.0024 !!!!!
  %%%% max_Delta_lambda = 0.0721 !!!!!!
  % for parameter m = 256
  end % if 0
  
  disp('searching maximum rb error over train set')
  max_err_u = 0;
  max_err_lambda = 0;
  K = model.get_inner_product_matrix(model,model_data);
  for m = 1:289
    model = set_mu(model,detailed_data.RB_info.mu_train(m,:));
    rb_sim_data = rb_simulation(model,reduced_data);
    rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
    sim_data = detailed_simulation(model,model_data);
    err_u = sqrt((rb_sim_data.U-sim_data.U)' * K * (rb_sim_data.U- ...
                                                    sim_data.U));
    L_R = K \ (rb_sim_data.L-sim_data.L);
    err_lambda = sqrt((rb_sim_data.L-sim_data.L)' * L_R);
    if max_err_u < err_u
      max_err_u = err_u
      disp(['found for m = ',num2str(m)]);
    end;
    if max_err_lambda < err_lambda
      max_err_lambda = err_lambda
      disp(['found for m = ',num2str(m)]);
    end;
  end;
  disp('maximum rb errors:')
  max_err_u
  max_err_lambda
  %%%%%%%% RESULT with standard quadprog options:
  %%%% max_err_u =      0.0010 !!!!!
  %%%% max_err_lambda = 0.0313 !!!!!!
  % for parameter m = 256!!!
  % 
  % NEW RESULTS with MaxIter=1000 as quadprog options :-) !!!
  %max_err_u = 4.4909e-009
  %max_err_lambda = 1.3476e-007
  
 case 1.6 % check accuracy of rb-scheme for large Lagrange basis: 289 vectors
          % IMPORTANT STEP HIGHLIGHTING ACCURACY PROBLEM

  model.RB_numintervals = [16,16];
  model.RB_generation_mode = 'pure_snapshots'; % no supr.
  model.enable_error_estimator = 1;
%  model.detailed_simulation = @my_new_detailed_simulation;
  disp('generating Lagrange basis');
  detailed_data = filecache_function(@gen_detailed_data,model,model_data);
  save('lagrange-basis-detailed-data','model','detailed_data');
  
  reduced_data = gen_reduced_data(model,detailed_data);
  
  K = model.get_inner_product_matrix(model,model_data);
  m  = 256;
  model = set_mu(model,detailed_data.RB_info.mu_train(m,:));
  rb_sim_data = rb_simulation(model,reduced_data);
  rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
  sim_data = detailed_simulation(model,model_data);
  err_u = sqrt((rb_sim_data.U-sim_data.U)' * K * (rb_sim_data.U- ...
                                                  sim_data.U));
  L_R = K \ (rb_sim_data.L-sim_data.L);
  err_lambda = sqrt((rb_sim_data.L-sim_data.L)' * L_R);
  disp('rb errors for single nasty sample:')
  err_u
  err_lambda
  %%%%%%%% RESULT for standard QUADPROG parameter:
  %%%% max_err_u =      0.0010 !!!!!
  %%%% max_err_lambda = 0.0313 !!!!!!
  %%%%%%%% RESULT with maxiter = 1000:
  % err_u =  5.2329e-013
  % err_lambda =  1.0892e-011
  % Very nice, but 5 seconds "reduced" computation time.

 case 1.7 % generate coarse lagrange basis
  
  model.RB_numintervals = [4,4];
  model.RB_generation_mode = 'pure_snapshots'; % no supr.
  model.enable_error_estimator = 1;
  model.supremizer_enrichment = 0;
  %  model.detailed_simulation = @my_new_detailed_simulation;
  disp('generating Lagrange basis');
  detailed_data = filecache_function(@gen_detailed_data,model,model_data);
  save('lagrange-basis-5x5-no-supr-detailed-data','model','detailed_data');
  model.supremizer_enrichment = 1;
  save('lagrange-basis-5x5-with-supr-detailed-data','model','detailed_data');
 
 case 2 % greedy basis generation for paper!

  %  model.RB_numintervals = [4,4]; % works perfect: N=25 exact representation
  %  model.RB_numintervals = [8,8]; 
  model.RB_numintervals = [16,16];
  % weights of error components
  model.w_u = 1;
  model.w_lambda = 1;
  model.RB_stop_Nmax = 210;
  model.RB_stop_epsilon = 1e-8; 
  model.supremizer_enrichment = 1; % on/off
                                       
  % WN2: 
  % rb-estimator with W subbasis for dual reduced basis:
  model.RB_generation_mode = 'greedy-rb-estimator';
  model.supremizer_enrichment = 1; % on/off
  model.enable_error_estimator = 1;
  model.RB_stop_Nmax = 10; % 4 suffice here!
  %model.RB_numintervals = [16,16];
  model.dummy = 2;
  model.dual_lin_independence_mode = 1; % use W  subbasis 
  
  disp('--------------------------------------------------');
  disp(model.RB_generation_mode);
  detailed_data = filecache_function(@gen_detailed_data,model,model_data);
  save('greedy-rb-estimator-with-supr-detailed-data-lin-indep-1',...
       'model','detailed_data');
  max_err_sequence = detailed_data.RB_info.max_err_sequence;
  plot(max_err_sequence);
  set(gca,'Yscale','log');
  title('greedy rb-estimator, with supr.-enrichment, lin-indep=1');
  disp('press enter to continue');
  plot_detailed_data(model,detailed_data,[]);
  keyboard;
  
  % WN1: 
  model.RB_generation_mode = 'greedy-rb-estimator';
  model.dual_lin_independence_mode = 0; % use dual snapshots 
  model.supremizer_enrichment = 1; % on/off
  model.enable_error_estimator = 1;
%  model.RB_stop_Nmax = 300; 
  model.RB_stop_Nmax = 50; 
  disp('--------------------------------------------------');
  disp(model.RB_generation_mode);
  detailed_data = filecache_function(@gen_detailed_data,model,model_data);
  save('greedy-rb-estimator-with-supr-detailed-data-50',...
       'model','detailed_data');
  max_err_sequence = detailed_data.RB_info.max_err_sequence;
  plot(max_err_sequence);
  set(gca,'Yscale','log');
  title('greedy rb-estimator convergence, with supr.-enrichment');
  disp('press enter to continue');
  plot_detailed_data(model,detailed_data,[]);
  keyboard;
  
 case 2.1 % greedy error convergence for different error indicators
  
%  model.RB_numintervals = [4,4]; % works perfect: N=25 exact representation
%  model.RB_numintervals = [8,8]; 
  model.RB_numintervals = [16,16];
  % weights of error components
  model.w_u = 1;
  model.w_lambda = 1;
  model.RB_stop_Nmax = 210;
  model.RB_stop_epsilon = 1e-8; 
  model.supremizer_enrichment = 1; % on/off
  model.enable_error_estimator = 0;
  model.dual_lin_independence_mode = 0; % use dual snapshots
%  model.dual_lin_independence_mode = 1; % use W  basis vectors
 
  %  model.RB_generation_mode = 'greedy-true-projection-error';
  %  disp('--------------------------------------------------');
  %  disp(model.RB_generation_mode);
  %  detailed_data = filecache_function(@gen_detailed_data,model, ...
  %                                  model_data);
  %  save('greedy-true-projection-error-detailed-data','model','detailed_data');
  %  max_err_sequence1 = detailed_data.RB_info.max_err_sequence;
  %  plot(max_err_sequence1);
  %  set(gca,'Yscale','log');
  %  title('greedy true projection error convergence');
  %  %  plot_detailed_data(model,detailed_data,[]);
  %  keyboard;
  %  disp('press enter to continue');
  %%  pause;

  % rb-estimator with W subbasis for dual reduced basis:
  model.RB_generation_mode = 'greedy-rb-estimator';
  model.supremizer_enrichment = 1; % on/off
  model.enable_error_estimator = 1;
  model.RB_stop_Nmax = 20; % 4 suffice here!
  model.RB_numintervals = [128,128];
  model.dual_lin_independence_mode = 1; % use W  subbasis 
  
  disp('--------------------------------------------------');
  disp(model.RB_generation_mode);
  detailed_data = filecache_function(@gen_detailed_data,model,model_data);
  save('greedy-rb-estimator-with-supr-detailed-data-lin-indep-1-128sqr',...
       'model','detailed_data');
  max_err_sequence7 = detailed_data.RB_info.max_err_sequence;
  plot(max_err_sequence7);
  set(gca,'Yscale','log');
  title('greedy rb-estimator, with supr.-enrichment, lin-indep=1');
  disp('press enter to continue');
  plot_detailed_data(model,detailed_data,[]);
  keyboard;

  % rb-estimator with W subbasis for dual reduced basis:
  model.RB_generation_mode = 'greedy-rb-estimator';
  model.supremizer_enrichment = 1; % on/off
  model.enable_error_estimator = 1;
  model.RB_stop_Nmax = 10; % 4 suffice here!
  model.RB_numintervals = [16,16];
  model.dummy = 1;
  model.dual_lin_independence_mode = 1; % use W  subbasis 
  
  disp('--------------------------------------------------');
  disp(model.RB_generation_mode);
  detailed_data = filecache_function(@gen_detailed_data,model,model_data);
  save('greedy-rb-estimator-with-supr-detailed-data-lin-indep-1',...
       'model','detailed_data');
  max_err_sequence7 = detailed_data.RB_info.max_err_sequence;
  plot(max_err_sequence7);
  set(gca,'Yscale','log');
  title('greedy rb-estimator, with supr.-enrichment, lin-indep=1');
  disp('press enter to continue');
  plot_detailed_data(model,detailed_data,[]);
  keyboard;

    % rb-estimator with W subbasis for dual reduced basis:
  model.RB_generation_mode = 'greedy-rb-estimator';
  model.supremizer_enrichment = 0; % on/off
  model.enable_error_estimator = 1;
  model.RB_stop_Nmax = 100; % must suffice here!
  model.RB_numintervals = [16,16];
  model.dual_lin_independence_mode = 1; % use W  subbasis 
  
  disp('--------------------------------------------------');
  disp(model.RB_generation_mode);
  detailed_data = filecache_function(@gen_detailed_data,model,model_data);
  save('greedy-rb-estimator-no-supr-detailed-data-lin-indep-1',...
       'model','detailed_data');
  max_err_sequence7 = detailed_data.RB_info.max_err_sequence;
  plot(max_err_sequence7);
  set(gca,'Yscale','log');
  title('greedy rb-estimator, no supr.-enrichment, lin-indep=1');
  disp('press enter to continue');
  plot_detailed_data(model,detailed_data,[]);
  keyboard;

  
  
  
  model.RB_numintervals = [16,16];

  

  % rb-estimator with real large training data set
  model.RB_generation_mode = 'greedy-rb-estimator';
  model.supremizer_enrichment = 1; % on/off
  model.enable_error_estimator = 1;
  model.RB_stop_Nmax = 100; % 100 must suffice for now.
  model.RB_numintervals = [128,128];
  disp('--------------------------------------------------');
  disp(model.RB_generation_mode);
  detailed_data = filecache_function(@gen_detailed_data,model,model_data);
  save('greedy-rb-estimator-with-supr-detailed-data-128sqr_train',...
       'model','detailed_data');
  max_err_sequence7 = detailed_data.RB_info.max_err_sequence;
  plot(max_err_sequence7);
  set(gca,'Yscale','log');
  title('greedy rb-estimator convergence, with supr.-enrichment');
  disp('press enter to continue');
  plot_detailed_data(model,detailed_data,[]);
  keyboard;

  model.RB_numintervals = [16,16];
  
  model.RB_generation_mode = 'greedy-rb-estimator';
  model.supremizer_enrichment = 1; % on/off
  model.w_u = 1;
  model.w_lambda = 0.01;
  model.enable_error_estimator = 1;
  model.RB_stop_Nmax = 50; 
  disp('--------------------------------------------------');
  disp(model.RB_generation_mode);
  detailed_data = filecache_function(@gen_detailed_data,model,model_data);
  save('greedy-rb-estimator-with-supr-detailed-data-50-w0p01',...
       'model','detailed_data');
  max_err_sequence = detailed_data.RB_info.max_err_sequence;
  plot(max_err_sequence);
  set(gca,'Yscale','log');
  title('greedy rb-estimator convergence, with supr.-enrichment');
  disp('press enter to continue');
  plot_detailed_data(model,detailed_data,[]);
  keyboard;

  model.w_u = 1;
  model.w_lambda = 1;
 
  model.RB_generation_mode = 'greedy-rb-estimator';
  model.supremizer_enrichment = 1; % on/off
  model.enable_error_estimator = 1;
  model.RB_stop_Nmax = 300; 
  disp('--------------------------------------------------');
  disp(model.RB_generation_mode);
  detailed_data = filecache_function(@gen_detailed_data,model,model_data);
  save('greedy-rb-estimator-with-supr-detailed-data-300',...
       'model','detailed_data');
  max_err_sequence7 = detailed_data.RB_info.max_err_sequence;
  plot(max_err_sequence7);
  set(gca,'Yscale','log');
  title('greedy rb-estimator convergence, with supr.-enrichment');
  disp('press enter to continue');
  plot_detailed_data(model,detailed_data,[]);
  keyboard;
  
%  model.RB_generation_mode = 'greedy-rb-estimator';
%  model.supremizer_enrichment = 1; % on/off
%  model.enable_error_estimator = 1;
%  model.RB_stop_Nmax = 210; 
%  disp('--------------------------------------------------');
%  disp(model.RB_generation_mode);
%  detailed_data = filecache_function(@gen_detailed_data,model,model_data);
%  save('greedy-rb-estimator-with-supr-detailed-data','model','detailed_data');
%  max_err_sequence4 = detailed_data.RB_info.max_err_sequence;
%  plot(max_err_sequence4);
%  set(gca,'Yscale','log');
%  title('greedy rb-estimator convergence, supr.-enrichment');
%  disp('press enter to continue');
%  plot_detailed_data(model,detailed_data,[]);
%  keyboard;
  
%  model.RB_generation_mode = 'greedy-rb-estimator';
%  model.supremizer_enrichment = 0; % on/off
%  model.RB_stop_Nmax = 210; 
%  disp('--------------------------------------------------');
%  disp(model.RB_generation_mode);
%  detailed_data = filecache_function(@gen_detailed_data,model,model_data);
%  save('greedy-rb-estimator-no-supr-detailed-data','model','detailed_data');
%  max_err_sequence1 = detailed_data.RB_info.max_err_sequence;
%  plot(max_err_sequence1);
%  set(gca,'Yscale','log');
%  title('greedy rb-estimator convergence, no supr.-enrichment');
%  disp('press enter to continue');
%  plot_detailed_data(model,detailed_data,[]);
%  keyboard;
  
  model.RB_generation_mode = 'greedy-rb-estimator';
  model.supremizer_enrichment = 0; % on/off
  model.RB_stop_Nmax = 300; 
  disp('--------------------------------------------------');
  disp(model.RB_generation_mode);
  detailed_data = filecache_function(@gen_detailed_data,model,model_data);
  save('greedy-rb-estimator-no-supr-detailed-data-300',...
       'model','detailed_data');
  max_err_sequence5 = detailed_data.RB_info.max_err_sequence;
  plot(max_err_sequence5);
  set(gca,'Yscale','log');
  title('greedy rb-estimator convergence, no supr.-enrichment');
  disp('press enter to continue');
  plot_detailed_data(model,detailed_data,[]);
  keyboard;

  model.RB_generation_mode = 'greedy-true-rb-error';
  model.supremizer_enrichment = 0; % on/off
  model.RB_stop_Nmax = 210; % such that VN = V, WN = W;
  disp('--------------------------------------------------');
  disp(model.RB_generation_mode);
  detailed_data = filecache_function(@gen_detailed_data,model,model_data);
  save('greedy-true-rb-error-no-supr-detailed-data','model','detailed_data');
  max_err_sequence3 = detailed_data.RB_info.max_err_sequence;
  plot(max_err_sequence3);
  set(gca,'Yscale','log');
  title('greedy true rb-error convergence, no supremizers');
  %  plot_detailed_data(model,detailed_data,[]);
  disp('press enter to continue');
  keyboard;
%  pause;

  model.RB_generation_mode = 'greedy-true-rb-error';
  model.supremizer_enrichment = 1; % on/off
  model.RB_stop_Nmax = 210; % takes too long for larger values
  disp('--------------------------------------------------');
  disp(model.RB_generation_mode);
  detailed_data = filecache_function(@gen_detailed_data,model,model_data);
  save('greedy-true-rb-error-with-supr-detailed-data','model','detailed_data');
  max_err_sequence2 = detailed_data.RB_info.max_err_sequence;
  plot(max_err_sequence2);
  set(gca,'Yscale','log');
  title('greedy true rb-error convergence, supr.-enrichment');
  disp('press enter to continue');
  plot_detailed_data(model,detailed_data,[]);
  keyboard;
  pause;

  % plot all convergence curves: non comparable!!!!!
%  keyboard;
%  figure;
%  plot([max_err_sequence1;max_err_sequence3;max_err_sequence2;max_err_sequence4]')
%  legend('true-projection-error','true-rb-error-wo-supr',...
%        'true-rb-error-with-supr',...
%        'rb-estimator');
%  
%  plot([max_err_sequence3;max_err_sequence2;max_err_sequence4]')
%  legend('true-rb-error-wo-supr',...
%        'true-rb-error-with-supr',...
%        'rb-estimator-with-supr');
  
 case 2.2
  % separate greedy: first u then lambda based on projection error 
  
  disp('please implement monotonicity check on lamdba before running!');
  keyboard;
  
  %  model.RB_numintervals = [4,4]; % works perfect: N=25 exact representation
  %  model.RB_numintervals = [8,8]; 

  model.RB_numintervals = [16,16];
  % weights of error components
  model.w_u = 1;
  model.w_lambda = 1;
%  model.RB_stop_Nmax = 100;
  model.RB_stop_Nmax = 100;
  model.RB_stop_epsilon = 1e-13; % very accurate! 
  model.supremizer_enrichment = 0; % on/off
  model.enable_error_estimator = 0;
 
  model.RB_generation_mode = 'separate-greedy';
  disp('--------------------------------------------------');
  disp(model.RB_generation_mode);
  detailed_data = filecache_function(@gen_detailed_data,model, ...
                                     model_data);
  detailed_data = gen_detailed_data(model, model_data);
  save('separate-greedy-detailed-data','model','detailed_data');
  figure,plot(detailed_data.RB_info.max_u_err_sequence);
  set(gca,'Yscale','log');
  title('greedy u projection error convergence');
  %  plot_detailed_data(model,detailed_data,[]);
  figure,plot(detailed_data.RB_info.max_lambda_err_sequence);
  set(gca,'Yscale','log');
  title('greedy lambda projection error convergence');
  %  plot_detailed_data(model,detailed_data,[]);
  
  disp('press enter to continue');
  keyboard;
%  pause;

 case 3 % check all bases from 2!! interactive demo of mu distribution
  
%  fn =         'greedy-true-rb-error-with-supr-detailed-data'; % N_V = 200, N_S = 210
%  fn =         'greedy-true-rb-error-no-supr-detailed-data'; % N_V = 200, N_S = 210
  %fn = 'greedy-rb-estimator-no-supr-detailed-data-300.mat'; % N_V = 200, N_S = 289
%  fn = 'greedy-rb-estimator-with-supr-detailed-data'; % N_V = 200, N_S = 210
%  fn = 'greedy-rb-estimator-with-supr-detailed-data-300';
%  fn = 'separate-greedy-detailed-data'; % N_S = 100, N_V = 70;
%  fn = 'greedy-rb-estimator-with-supr-detailed-data-300';
%  fn = 'greedy-rb-estimator-with-supr-detailed-data-50-w0p01';

  fn = 'greedy-rb-estimator-with-supr-detailed-data-lin-indep-1';
    
  load(fn);
  
  plot_detailed_data(model,detailed_data);

  if isfield(detailed_data.RB_info,'mu_sequence');
    mu_sequence = detailed_data.RB_info.mu_sequence;
    inspect_mu_distribution(mu_sequence,model);
  end;
  if isfield(detailed_data.RB_info,'mu_u_sequence');
    mu_sequence = detailed_data.RB_info.mu_u_sequence;
    inspect_mu_distribution(mu_sequence,model);
  end;
  if isfield(detailed_data.RB_info,'mu_lambda_sequence');
    mu_sequence = detailed_data.RB_info.mu_lambda_sequence;
    inspect_mu_distribution(mu_sequence,model);
  end;
  
  keyboard;
  
 case 3.1 % investigate basis from 2
  
  load('greedy-true-rb-error-detailed-data','model','detailed_data');
  
  % detect duplicate lambda columns:
  % did not find...

  % compute gramian matrix of lambda snapshots
  
  K = model.get_inner_product_matrix(model,model_data);
  L = detailed_data.RB_L;
  GL = L' * (K\ L);
  GL = (GL+GL')/2;
  e = eig(GL);
  plot(e(end:-1:1));
  set(gca,'Yscale','log')
  
  % about 60 dominant eigenvalues => 66 dimensional space should suffice
  % but larger set to be expected.
 
 case 3.2 % investigate greedy-separate basis from 2.1

  % test error for maximum accuracy
  % without supremizers (supremizers not sufficient for stabilita
  % due to different snapshot sets.)
   
  % is error comparable to train projection error?

  % test rb-error: comparable to traing error?
  
 case 4 % determine statistics of the different bases from 2
  
  % FoR PAPER:
  
  %fn = 'greedy-rb-estimator-with-supr-detailed-data-lin-indep-1';
  %par.numintervals = [16,16]; % for training data
  %Ns = [1:10];

  fn = 'greedy-rb-estimator-with-supr-detailed-data-50';
  par.numintervals = [16,16]; % for training data
  Ns = [1:50];

  % OLD:
  
    %  fn='separate-greedy-detailed-data';
  %  par.numintervals = [16,16]; % for training data
  %  Ns = [1:50, 60,70,80];


%  fn = 'greedy-rb-estimator-with-supr-detailed-data-300';
%  par.numintervals = [16,16]; % for training data
%  Ns = [1:50, 60:10:280, 289, 290:10:300];
  
  %  fn = 'lagrange-basis-5x5-with-supr-detailed-data';
  %  %   'lagrange-basis-5x5-no-supr-detailed-data'};
  %    par.numintervals = [4,4]; % for training data
  
  % 'greedy-true-projection-error-detailed-data'
  %  fns ={...
  %,...
  %  fns = {'greedy-rb-estimator-with-supr-detailed-data'};
  %    'greedy-true-rb-error-no-supr-detailed-data',...
  % fns = { 'lagrange-basis-detailed-data'};
  %     'greedy-rb-estimator-no-supr-detailed-data',...
  %     'greedy-true-rb-error-with-supr-detailed-data',...
  %
  % Desired statistics for evaluation of basis:
  % maximum train rb-error for different basis sizes
  % mean train rb-error for different basis sizes + std
  % maximum test rb-error for different basis sizes
  % mean test rb-error for different basis sizes + std
  % maximum train rb-error-estimator for different basis sizes
  % mean train rb-error-estimator for different basis sizes + std
  % maximum test rb-error-estimator for different basis sizes
  % mean test rb-error-estimator for different basis sizes + std
  % average online runtime 
  % tilde_N_V, N_V, N_S
  % List of Test-parameters

  % generate train and test data
  
  % for rb-error on training-grid:
  %par.numintervals = [16,16];
  par.range = model.mu_ranges; 
  mu_mesh = cubegrid(par);
  mu.train = get(mu_mesh,'vertex')';
  % for rb-error on test-grid:
  par.numintervals = [32,32];
  mu_mesh = cubegrid(par);
  mu.test = get(mu_mesh,'vertex')';
  %  mu_test = rand_uniform(ntest,model.mu_ranges);
    
  modes = {'test','train'};
  error_indicators = {'estimator','error'};
  
  % initialize empty statistics fields of bases if not existing
  %for fni = 1:length(fns)
  %  fn = fns{fni};
    clear('model','detailed_data');
    load(fn);
    if ~isfield(model,'supremizer_enrichment')
      model.supremizer_enrichment = 1;
    end;
    if ~isfield(model,'dual_lin_independence_mode')
      model.dual_lin_independence_mode = 0;
    end;
    
    % computing rb sizes of ceil(Nfact*N_S), etc.
    %Nfactors = [0:0.05:1]; % 21 sampling points should be fine
    % compute Nfactors from desired N_S list:
    if ~isfield(detailed_data,'L_no_post_processing')
      N_S = size(detailed_data.RB_L,2);
      Ns = Ns(find(Ns<=N_S));
      if Ns(end)~=N_S
        Ns = [Ns, N_S];
      end;
    end;
%    Nfactors = Ns/N_S;

%    Nfactors = [0,0.5,1];
    
%    nfact = length(Nfactors);
    nfact = length(Ns);

    if ~isfield(detailed_data,'statistics')
%      detailed_data.statistics.Nfactors = Nfactors;
      detailed_data.statistics.Ns = Ns;
      detailed_data.statistics.train.error.absolute.max = NaN * Ns;
      detailed_data.statistics.train.error.absolute.mean = NaN * Ns;
      detailed_data.statistics.train.error.absolute.std = NaN * Ns;
      detailed_data.statistics.train.estimator.absolute.max = NaN * Ns;
      detailed_data.statistics.train.estimator.absolute.mean = NaN * Ns;
      detailed_data.statistics.train.estimator.absolute.std = NaN * Ns;
      detailed_data.statistics.test.error.absolute.max = NaN * Ns;
      detailed_data.statistics.test.error.absolute.mean = NaN * Ns;
      detailed_data.statistics.test.error.absolute.std = NaN * Ns;
      detailed_data.statistics.test.estimator.absolute.max = NaN * Ns;
      detailed_data.statistics.test.estimator.absolute.mean = NaN * Ns;
      detailed_data.statistics.test.estimator.absolute.std = NaN * Ns;
      detailed_data.statistics.train.error.relative.max = NaN * Ns;
      detailed_data.statistics.train.error.relative.mean = NaN * Ns;
      detailed_data.statistics.train.error.relative.std = NaN * Ns;
      detailed_data.statistics.train.estimator.relative.max = NaN * Ns;
      detailed_data.statistics.train.estimator.relative.mean = NaN * Ns;
      detailed_data.statistics.train.estimator.relative.std = NaN * Ns;
      detailed_data.statistics.test.error.relative.max = NaN * Ns;
      detailed_data.statistics.test.error.relative.mean = NaN * Ns;
      detailed_data.statistics.test.error.relative.std = NaN * Ns;
      detailed_data.statistics.test.estimator.relative.max = NaN * Ns;
      detailed_data.statistics.test.estimator.relative.mean = NaN * Ns;
      detailed_data.statistics.test.estimator.relative.std = NaN * Ns;
      detailed_data.statistics.runtime.std = NaN * Ns;
      detailed_data.statistics.runtime.mean = NaN * Ns;
      detailed_data.statistics.tildeNV = NaN * Ns;
      detailed_data.statistics.NV = NaN * Ns;
      detailed_data.statistics.NS = NaN * Ns;
      detailed_data.statistics.NW = NaN * Ns;
      detailed_data.statistics.test.mu_sequence = mu.test ;
      detailed_data.statistics.train.mu_sequence = mu.train ;
    end;
    save(fn,'model','detailed_data');
%  end;
  
  % for all nfactors in bisecting order:
  %  for nfi = 1:length(Nfactors);
  for nfi = 1:nfact;
    %  for nfi = [21,1,10,5,15,3,7,13,18,2,4,6,8,9,11,12,14,16,17,19,20];
    %  for nfi = 1:length(Nfactors)
    %    Nfactor = Nfactors(nfi);
    %    disp(['Starting evaluation for Nfactor ',num2str(Nfactor)]);
    NS = Ns(nfi);
    disp(['Starting evaluation for N ',num2str(NS)]);
%    for fni = 1:length(fns)
%      fn = fns{fni};
      disp(['  Starting evaluation for file ',fn]);
%      load(fn);
      K = model.get_inner_product_matrix(model,model_data); 
      
      for modei = 2:-1:1 % train / test
        mode = modes{modei};
        disp(['    Starting mode ',mode]);
        disp('       computing all snapshots');
        mu_sequence = getfield(mu,mode); % get train/test-sequence
        nsamples = size(mu_sequence,2);
        [Utest,Ltest] = filecache_function(...
            @compute_test_snapshots,model,model_data,mu_sequence);
        
        %       for erri = 1:2
        %         err_indicator = error_indicators{erri};
        model.enable_error_estimator = 1;
        %         disp(['       Starting error indicator: ', ...
        %               err_indicator]);        
        err_struct = getfield(...
            detailed_data.statistics,mode);
        
        if isnan(err_struct.error.absolute.mean(nfi)) % if result not already existing
          det_data = detailed_data;
          if isfield(detailed_data,'L_no_post_processing')
            tildeNVmax = size(detailed_data.RB_U_no_supr,2);
            NSmax = size(detailed_data.L_no_post_processing,2);
            Nfactor = NS/NSmax;
            tildeNV = ceil(tildeNVmax*Nfactor);     
%           NS = ceil(NSmax*Nfactor);
            U_no_supr = det_data.RB_U_no_supr(:,1:tildeNV);
          elseif isfield(detailed_data,'RB_U_no_supr')
            tildeNVmax = size(detailed_data.RB_U_no_supr,2);
            NSmax = size(detailed_data.RB_L,2);
            Nfactor = NS/NSmax;
            tildeNV = ceil(tildeNVmax*Nfactor);
%           NS = ceil(NSmax*Nfactor);
            U_no_supr = det_data.RB_U_no_supr(:,1:tildeNV);
          else
            tildeNVmax = size(detailed_data.RB_U,2);       
            NSmax = size(detailed_data.RB_L,2);
            Nfactor = NS/NSmax;
            tildeNV = ceil(tildeNVmax*Nfactor);
            %       NS = ceil(NSmax*Nfactor);
            U_no_supr =  det_data.RB_U(:,1:tildeNV);
          end;
          det_data.RB_L = det_data.RB_L(:,1:NS);
          if model.dual_lin_independence_mode == 1
            det_data.RB_L = det_data.L_no_post_processing(:,1:NS);
            dual_speye = speye(size(det_data.RB_L,1));
            L_dofs = find(max(abs(det_data.RB_L),[],2)>1e-11);
            det_data.RB_L = full(dual_speye(:,L_dofs));
          end;
          NW = size(det_data.RB_L,2);
          if model.supremizer_enrichment
            det_data.RB_U = improve_conditioning(...
                [U_no_supr,K\det_data.RB_L],K);
            NV = size(det_data.RB_U,2);
          else % no supr.-enrichment
            det_data.RB_U = improve_conditioning(U_no_supr,K);
            NV = tildeNV;
          end;
          reduced_data = gen_reduced_data(model,det_data);
          err_us = zeros(1,nsamples);
          err_lambdas = zeros(1,nsamples);
          delta_us = zeros(1,nsamples);
          delta_lambdas = zeros(1,nsamples);
          norm_us = zeros(1,nsamples);
          norm_lambdas = zeros(1,nsamples);
          runtimes = zeros(1,nsamples);
          for mui = 1:nsamples
            tic;
            fprintf('.');
            model = model.set_mu(model,mu_sequence(:,mui));
            rb_sim_data = rb_simulation(model,reduced_data);
            rb_sim_data = rb_reconstruction(model,det_data,rb_sim_data);
            err_us(mui) = sqrt((rb_sim_data.U-Utest(:,mui))'*K*...
                               (rb_sim_data.U-Utest(:,mui)));
            err_lambdas(mui) = sqrt((rb_sim_data.L-Ltest(:,mui))'*(K\ ...
                                                  (rb_sim_data.L- ...
                                                   Ltest(:,mui))));
            norm_us(mui) = sqrt(Utest(:,mui)'*K*Utest(:,mui));
            norm_lambdas(mui) = sqrt(Ltest(:,mui)'*(K\Ltest(:,mui)));
            delta_us(mui) = rb_sim_data.Delta_u;
            delta_lambdas(mui) = rb_sim_data.Delta_lambda;
%           if delta_lambdas(mui)<err_lambdas(mui)
%             disp('check accuracy!');
%             keyboard;
%           end;
%           if delta_us(mui)<err_us(mui)
%             disp('check accuracy!');
%             keyboard;
%           end;
            runtimes(mui) = toc;
          end; % end mu loop
          
          % compute relative errors
          rel_err_us = err_us;
          rel_delta_us = delta_us;
          i = find(norm_us~=0);
          if ~isempty(i)
            rel_err_us(i) = err_us(i).*(norm_us(i).^-1);
            rel_delta_us(i) = delta_us(i).*(norm_us(i).^-1);
          end;
          i = find(norm_us==0);
          if ~isempty(i)
            rel_err_us(i) = NaN;
            rel_delta_us(i) = NaN;
          end;
          
          rel_err_lambdas = err_lambdas;
          rel_delta_lambdas = delta_lambdas;
          i = find(norm_lambdas~=0);
          if ~isempty(i)
            rel_err_lambdas(i) = err_lambdas(i).*(norm_lambdas(i).^-1);
            rel_delta_lambdas(i) = delta_lambdas(i).*(norm_lambdas(i).^-1);
          end;
          i = find(norm_lambdas==0);
          if ~isempty(i)
            rel_err_lambdas(i) = NaN;
            rel_delta_lambdas(i) = NaN;
          end;
          
          err_struct.error.absolute.max(nfi) = max(err_us + err_lambdas);
          err_struct.error.absolute.mean(nfi)  = mean(err_us + err_lambdas);
          err_struct.error.absolute.std(nfi)  = std(err_us +  err_lambdas);
          err_struct.error.relative.max(nfi)  = max(rel_err_us + rel_err_lambdas);
          err_struct.error.relative.mean(nfi)  = max(rel_err_us + rel_err_lambdas);
          err_struct.error.relative.std(nfi)  = std(rel_err_us + rel_err_lambdas);

          err_struct.estimator.absolute.max(nfi)  = max(delta_us + delta_lambdas);
          err_struct.estimator.absolute.mean(nfi)  = mean(delta_us + delta_lambdas);
          err_struct.estimator.absolute.std(nfi)  = std(delta_us + delta_lambdas);
          err_struct.estimator.relative.max(nfi)  = max(rel_delta_us + rel_delta_lambdas);
          err_struct.estimator.relative.mean(nfi)  = mean(rel_delta_us + rel_delta_lambdas);
          err_struct.estimator.relative.std(nfi)  = std(rel_delta_us + rel_delta_lambdas);
                  
          detailed_data.statistics.tildeNV(nfi) = tildeNV; 
          detailed_data.statistics.NV(nfi) = NV; 
          detailed_data.statistics.NW(nfi) = NW; 
          detailed_data.statistics.NS(nfi) = NS;
          if (modei==1) 
            detailed_data.statistics.runtime.mean(nfi) = mean(runtimes); 
            detailed_data.statistics.runtime.std(nfi) = std(runtimes); 
          end;
          eval(['detailed_data.statistics.',mode,' = err_struct;']); 
          %         disp('check correct setting of detailed_data');
          %         keyboard;
          fprintf('\n');
          save(fn,'model','detailed_data');
        end; % if result not already existing
      end; % modei
    end;
%  end;
  
 case 5 % plots of different statistics from case 4

  fn = 'greedy-rb-estimator-with-supr-detailed-data-lin-indep-1';
  load(fn);
  stat = detailed_data.statistics
  p1 = plot(stat.NS,stat.train.error.absolute.max,'ko');
  hold on;
  p2 = plot(stat.NS,stat.train.estimator.absolute.max,'b*');
%  p3 = plot(detailed_data.RB_info.max_err_sequence);
  legend({'error','estimator'});
  xlabel('cone spanning elements');
  ylabel('accuracy');
  title('training set error/estimator convergence');
  set(gca,'Yscale','log');
    
  % test figure
  figure;
  fn = 'greedy-rb-estimator-with-supr-detailed-data-lin-indep-1';
  load(fn);
  stat = detailed_data.statistics
  p1 = plot(stat.NS,stat.test.error.absolute.max,'ko');
  hold on;
  p2 = plot(stat.NS,stat.test.estimator.absolute.max,'b*');
  legend({'error','estimator'});
  xlabel('cone spanning elements');
  ylabel('accuracy');
  title('test set error/estimator convergence');
  set(gca,'Yscale','log');
  keyboard;

    
  % separate greedy basis
  
    % train figure
  figure;
  fn = 'separate-greedy-detailed-data';
  load(fn);
  stat = detailed_data.statistics
  p1 = plot(stat.NS,stat.train.error.absolute.max,'ko');
  hold on;
  p2 = plot(stat.NS,stat.train.estimator.absolute.max,'b*');
%  p3 = plot(detailed_data.RB_info.max_err_sequence);
  legend({'error','estimator'});
  xlabel('cone spanning elements N_S');
  ylabel('accuracy');
  title('training set error/estimator convergence');
  set(gca,'Yscale','log');
    
  % test figure
  figure;
  fn = 'separate-greedy-detailed-data';
  load(fn);
  stat = detailed_data.statistics
  p1 = plot(stat.NS,stat.test.error.absolute.max,'ko');
  hold on;
  p2 = plot(stat.NS,stat.test.estimator.absolute.max,'b*');
  legend({'error','estimator'});
  xlabel('cone spanning elements N_S');
  ylabel('accuracy');
  title('test set error/estimator convergence');
  set(gca,'Yscale','log');
  keyboard;

 case 6
  
  % PLOTS FOR PAPER!!!

    % test figure
  figure;
  fn = 'greedy-rb-estimator-with-supr-detailed-data-50';
%  fn = 'greedy-rb-estimator-with-supr-detailed-data-300';
  load(fn);
  stat = detailed_data.statistics
  p1 = plot(stat.NS(1:50),stat.test.error.absolute.max(1:50),'-b+');
  set(p1,'Markersize',10);
  hold on;
  p2 = plot(stat.NS(1:50),stat.train.error.absolute.max(1:50),'-.g*');
  set(p2,'Color',[0,0.5,0],'Linewidth',2,'Markersize',10);
  set(gca,'Fontsize',15);
  legend({'test-error','train-error'},'Fontsize',20);
  xlabel('N_W','Fontsize',20);
  ylabel('accuracy','Fontsize',20);
  title('test/train error convergence, W_N^{(1)}','Fontsize',20);
  set(gca,'Yscale','log');
  set(gca,'Ylim',[0.02,2]);
%  set(gca,'Yticklabel',['1e-01';'1e-00']);  
%  set(gca,'Xticklabel',[' 0';'  ';'10';'  ';'20';'  ';'30';'  ';'40';'  ';'50']);
%  saveas(gcf,'greedy-test-train-error-WN1-new.eps','epsc');
  saveas(gcf,'greedy-test-train-error-WN1-newest.eps','epsc');
  
  % plot of better basis  
  figure;
  fn = 'greedy-rb-estimator-with-supr-detailed-data-lin-indep-1';
  load(fn);
  stat = detailed_data.statistics
  p1 = plot(stat.NW(1:4),stat.test.error.absolute.max(1:4),'-b+');
  set(p1,'Markersize',10);
  hold on;
  p2 = plot(stat.NW(1:4),stat.train.error.absolute.max(1:4),'-.g*');
  set(p2,'Color',[0,0.5,0],'Linewidth',2,'Markersize',10);
  set(gca,'Fontsize',15);
  legend({'test-error','train-error'},'Location','SouthWest',...
         'Fontsize',20);
  xlabel('N_W','Fontsize',20);
  ylabel('accuracy','Fontsize',20);
  title('test/train error convergence, W_N^{(2)}','Fontsize',20);
  set(gca,'Yscale','log');
  set(gca,'Xlim',[10,60]);
  set(gca,'Ylim',[1e-14,11]);
%  set(gca,'Xticklabel',...
%         ['10';'  ';'20';'  ';'30';'  ';'40';'  ';'50';'  ';'60']);
%  set(gca,'Yticklabel',...
%         ['1e-14';'1e-12';'1e-10';'1e-08';'1e-06';'1e-04';'1e-02';'1e-00']);
%  saveas(gcf,'greedy-test-train-error-WN2-new.eps','epsc');
  saveas(gcf,'greedy-test-train-error-WN2-newest.eps','epsc');

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
  % OLD 
  % OLD 
  % OLD 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
  
 case 9
  
  % plot of maximum training rb-error and estimator over dimension N 
  %            for rb-estimator-wo-supr-300 
  % plot of maximum test rb-error and estimator over online-runtime 
  %            for rb-estimator-wo-supr-300 
  
  disp('rb-estimator-greedy-with-supr-300: train and test error decay:')
  % train figure
  figure;
  fn = 'greedy-rb-estimator-with-supr-detailed-data-300';
  load(fn);
  stat = detailed_data.statistics
  p1 = plot(stat.NS,stat.train.error.absolute.max,'ko');
  hold on;
  p2 = plot(stat.NS,stat.train.estimator.absolute.max,'b*');
  p3 = plot(detailed_data.RB_info.max_err_sequence);
  legend({'error','estimator','estimator'});
  xlabel('cone spanning elements N_S');
  ylabel('accuracy');
  title('training set error/estimator convergence');
  set(gca,'Yscale','log');
    
  % test figure
  figure;
  fn = 'greedy-rb-estimator-with-supr-detailed-data-300';
  load(fn);
  stat = detailed_data.statistics
  p1 = plot(stat.NS,stat.test.error.absolute.max,'ko');
  hold on;
  p2 = plot(stat.NS,stat.test.estimator.absolute.max,'b*');
  legend({'error','estimator'});
  xlabel('cone spanning elements N_S');
  ylabel('accuracy');
  title('test set error/estimator convergence');
  set(gca,'Yscale','log');


  disp('enter dbcont for continuing');
  keyboard;
  
  % => bad news: slow convergence and non-monotonocity...
  % => good news: error and estimator quite close!!
  % => good message: train error comparable to test error

  disp('Lagrange basis 5x5: nice training error, bad test error decay')
  
  fn = 'lagrange-basis-5x5-with-supr-detailed-data'; 
  %  'lagrange-basis-5x5-no-supr-detailed-data'};
  figure;
%  fn = 'greedy-rb-estimator-no-supr-detailed-data-300';
  load(fn);
  stat = detailed_data.statistics
  p1 = plot(stat.NS,stat.train.error.max,'ko');
  hold on;
  p2 = plot(stat.NS,stat.train.estimator.max,'b*');
%  p3 = plot(detailed_data.RB_info.max_err_sequence);
  legend({'error','estimator'});
  xlabel('cone spanning elements N_S');
  ylabel('accuracy');
  title('lagrange basis 5x5 training set error/estimator convergence');
  set(gca,'Yscale','log');
    
  % test figure
  figure;
%  fn = 'greedy-rb-estimator-no-supr-detailed-data-300';
  load(fn);
  stat = detailed_data.statistics
  p1 = plot(stat.NS,stat.test.error.max,'ko');
  hold on;
  p2 = plot(stat.NS,stat.test.estimator.max,'b*');
  legend({'error','estimator'});
  xlabel('cone spanning elements N_S');
  ylabel('accuracy');
  title('lagrange basis 5x5 test set error/estimator convergence');
  set(gca,'Yscale','log');

  disp('enter dbcont for continuing');
  keyboard;
  
  % comparison of different lagrange basis sizes
  
  fn = 'lagrange-basis-5x5-with-supr-detailed-data'; 
  %  'lagrange-basis-5x5-no-supr-detailed-data'};
  figure;
  %  fn = 'greedy-rb-estimator-no-supr-detailed-data-300';
  load(fn);
  stat = detailed_data.statistics
  p1 = plot(stat.NS,stat.test.error.max,'ko');
  hold on;
  p2 = plot(stat.NS,stat.test.estimator.max,'k*');
  fn = 'lagrange-basis-detailed-data'; 
  %  fn = 'greedy-rb-estimator-no-supr-detailed-data-300';
  load(fn);
  stat = detailed_data.statistics
  p1 = plot(stat.NS,stat.test.error.max,'ko');
  hold on;
  p2 = plot(stat.NS,stat.test.estimator.max,'k*');
%  p3 = plot(detailed_data.RB_info.max_err_sequence);
  legend({'error, 5x5','estimator 5x5','error 17x17','estimator 17x17'});
  xlabel('cone spanning elements N_S');
  ylabel('accuracy');
  title('lagrange-bases: test set error/estimator convergence');
  set(gca,'Yscale','log');
  
  % plot of train rb-estimator with and without supremizers o. online-runtime 
  % plot of test rb-estimator with and without supremizers over online-runtime 
  
  % ????
  
  % plot of train rb-estimator for greedy-rb-estim and
  %          greedy-rb-error wo supr over online-runtime 
  % plot of test rb-estimator for greedy-rb-estim and
  %          greedy-rb-error wo supr over online-runtime 
  
  % ????

  % plot of training rb-error for separate greedy basis
  % plot of test rb-error for separate greedy basis

  % ????

%  plot(Ns,errs)
%  legend(fns)
%  set(gca,'Yscale','log')
%  title(['rb-error decay on n=',num2str(ntest),' gridpoints']);
%  keyboard;    
%
% good news: rb_error estimator yields quite good rb-error 
% even better than greedy-proj-error approach!!!
  
 case 10 % interactive demo of reduced simulation
  
  %  demo_vi

  fn = 'greedy-rb-estimator-with-supr-detailed-data-50'; 
  %  fn = 'greedy-rb-estimator-detailed-data';
  load(fn);
  %  model = elastic_rope_model;
  model.enable_error_estimator = 0;
%  model.enable_error
  params.plot_title = 'elastic rope';
  params.no_lines = 0;
  params.axis_equal = 0;
  params.axis_tight = 1;
  params.solution_component = 'u';
%  params.solution_component = 'lambda';
  demo_rb_gui(model,detailed_data,[],params);
  %  demo_detailed_gui(model,detailed_data,[],params);

 case 11 % interactive demo of detailed simulation
    
%  fn = 'greedy-rb-estimator-detailed-data';
%  load(fn);
  detailed_data = gen_model_data(model);
  detailed_data.RB_U = [];
  detailed_data.RB_lambda = [];
  params.plot_title = 'elastic rope';
  params.no_lines = 0;
  params.axis_equal = 0;
  params.axis_tight = 1;
  params.solution_component = 'u';
%  params.solution_component = 'lambda';
  demo_detailed_gui(model,detailed_data,[],params);
  
 case 12 % interactive demo for error in u or lambda:
  fn = 'greedy-rb-estimator-detailed-data';
  load(fn);
  params.plot_title = 'elastic rope';
  params.no_lines = 0;
  params.axis_equal = 0;
  params.axis_tight = 1;
  % demo_rb_gui(model,detailed_data,[],params);
  % demo_detailed_gui(model,detailed_data,[],params);

  %  params.solution_component = 'u';
  % model.set_dofs_in_sim_data = @my_set_u_dofs_in_sim_data;

  % select here, whether u or lambda is to be visualized
  params.solution_component = 'lambda';
  model.set_dofs_in_sim_data = @my_set_lambda_dofs_in_sim_data;
  model.get_dofs_from_sim_data = @my_get_lambda_dofs_from_sim_data;
  demo_rb_error_gui(model,detailed_data,[],params);
 
 case 13 % debugging of detailed solver

  model = elastic_rope_model;
  model.xnumintervals = 202;
  model = set_mu(model,[0.075, 0.4])
  % detailed simulation and plot
  model_data = gen_model_data(model);
  disp('begin detailed simulation');
  sim_data = detailed_simulation(model,model_data);
  disp('plotting ...');
  plot_params.plot_title = 'detailed solution & obstacle';
%  plot_params.solution_component = 'lambda';
  p = plot_sim_data(model,model_data,sim_data,plot_params);
  save('debug_rope.mat','sim_data','model');
  keyboard;
  
 case 14
  julien = [];
  load A;
  julien.A = A(2:200,2:200); 
  load B;
  julien.B = B(2:200,2:200); 
  load f;
  julien.f = f(2:200); 
  load g;
  julien.g = g(2:200); 
  load U;
  julien.u = U(2:200); 
  load Lambda;
  julien.lambda = Lambda(2:200); 

  load debug_rope.mat;

  keyboard;

 case 15 % CPOD Basis generation 
  
%  model.RB_numintervals = [16,16];
%  % weights of error components
%  model.w_u = 1;
%  model.w_lambda = 1;
%%  model.RB_stop_Nmax = 210;
%  model.RB_stop_epsilon = 1e-8; 
%  model.supremizer_enrichment = 1; % on/off
  
%  model.RB_generation_mode = 'greedy-true-rb-error';
%  model.supremizer_enrichment = 1; % on/off
%  model.enable_error_estimator = 0;
%  model.RB_stop_Nmax = 10; % 4 suffice here!
                           
%  %model.RB_numintervals = [16,16];
%  model.dummy = 2;
%  model.dual_lin_independence_mode = 1; % use W  subbasis, i.e. WN^(2)
  
%  disp('--------------------------------------------------');
%  disp(model.RB_generation_mode);
%  detailed_data = filecache_function(@gen_detailed_data,model,model_data);
%  save('greedy-rb-estimator-with-supr-detailed-data-lin-indep-1',...
%       'model','detailed_data');
%  max_err_sequence = detailed_data.RB_info.max_err_sequence;
%  plot(max_err_sequence);
%  set(gca,'Yscale','log');
%  title('greedy rb-error, with supr.-enrichment, lin-indep=1');
%  disp('press enter to continue');
%  plot_detailed_data(model,detailed_data,[]);
%  keyboard;

  model.RB_generation_mode = 'CPOD_space_projection_error';
  
  %  NWs = [1,5,10,20]  
  %  for n = 1:length(NWs);
  model.RB_numintervals = [16,16];
  NW = 2;
  model.RB_stop_Nmax = NW;
  model.supremizer_enrichment = 1; % on/off
  model.dummy = rand(1);
  
  disp('--------------------------------------------------');
  disp(model.RB_generation_mode);
  detailed_data = filecache_function(@gen_detailed_data,model,model_data);
  save('CPOD-space-projection-with-supr-detailed-data',...
       'model','detailed_data');
  max_err_sequence = detailed_data.RB_info.max_err_sequence;
  plot(max_err_sequence);
  set(gca,'Yscale','log');
  title('CPOD-space-projection, with supr.-enrichment');
  disp('press enter to continue');
  plot_detailed_data(model,detailed_data,[]);
  keyboard;
  
 otherwise
  error('step unknown')
end;

% later: check projection-greedy basis and determine maximum RB
% error => what factor?
function sim_data = my_set_u_dofs_in_sim_data(sim_data,dofs) 
    sim_data.U = dofs;
function sim_data = my_set_lambda_dofs_in_sim_data(sim_data,dofs) 
    sim_data.L = dofs;
function dofs = my_get_lambda_dofs_from_sim_data(sim_data) 
    dofs = sim_data.L;

function [Utest,Ltest] = compute_test_snapshots(model,model_data,mu_test);
ntest = size(mu_test,2);
Utest = zeros(length(model_data.grid.X)-2,0);
Ltest = zeros(length(model_data.grid.X)-2,0);
for mui = 1:ntest
  fprintf('.');
  model = model.set_mu(model,mu_test(:,mui));
  sim_data = detailed_simulation(model,model_data);
  Utest = [Utest,sim_data.U];
  Ltest = [Ltest,sim_data.L];
end;
fprintf('\n');