adaptive_basisgen_mcmds.m

function adaptive_basisgen_mcmds(step)
%function adaptive_basisgen_mcmds(step)

% small script performing the experiments for the (second half) of
% the MCMDS paper. Experiments are based on a simple advection example
%
% mu_names  = {'cone_weight','vx_weight','vy_weight'};
% mu_ranges = [0,1]^3
%
% possible steps:
%
% step 1: detailed simulation of single trajectory, plot
%         POD and save reduced basis
%         reduced simulation and error

% Bernard Haasdonk 4.2.2010

if nargin < 1
  step = 1; 
end;

params = [];
model = [];
%params.coarse_factor = 2; % fact 2 => ndofs = 16384, nt = 1024, 15s detailed
%params.coarse_factor = 1; % fact 1 => ndofs = 32768, nt = 2048, 60s detailed
params.coarse_factor = 2; 
%verbose(1);
trajectory_fnbase = 'advection_trajectory';
detailed_data_fname = 'advection_detailed_data';
rb_fname = 'advection_rb';      

switch step
  
 case 1
  disp('initialization of lin_ds model from lin_evol model');
  %  verbose(1,'initialization of lin_ds model from lin_evol model and ...
  %       plot');
  
  ds_model = fast_model(params);
  ds_model.save_time_indices = 0:ds_model.nt;
  
  %ds_model.u_function_ptr = @(params) 0.1; % define new
%  ds_model.theta = 1;
  ds_model.theta = 0;
  ds_model_data = gen_model_data(ds_model);
%  mu = [0.5,1,1];
  mu = [0.5,1];
  % set mu in model and base_model!!!
  ds_model = ds_model.set_mu(ds_model,mu);
  ds_sim_data = detailed_simulation(ds_model,ds_model_data);

  % save trajectory
  save(fullfile(rbmatlabresult,[trajectory_fnbase,num2str(1)]),...
       'ds_model','ds_model_data','ds_sim_data');
  
  % plot trajectory
  params.no_lines = 1;
  params.clim = [0,1];
  params.plot_title='detailed implicit lin_ds';
  plot_sim_data(ds_model,ds_model_data,ds_sim_data,params);  
  
  disp('PCA of snapshots takes a while...');
  disp('Loading of trajectory...');
  load(fullfile(rbmatlabresult,[trajectory_fnbase,'1']),...
       'ds_model','ds_model_data','ds_sim_data');
  disp('...done.');
  % smaller epsilon does result in more than 1 vector for mu=0,0,0
  epsilon = 1e-6;
  disp('Calling PCA_Fixspace...');
  RB = PCA_fixspace(ds_sim_data.X,[],ds_model_data.G,150,[],epsilon);
  disp('...done.');
  %RB = orthonormalize_qr(ds_sim_data.X,ds_model_data.G);
  % select correct orthogonal set
  K = RB' * ds_model_data.G * RB;
  i = find(abs(diag(K)-1)>1e-2);
  i = min(i);
  disp(['found ',num2str(size(RB,2)),' basis vectors after traj. 1']);
  if isempty(i)
    i = size(RB,2)+1;
  end;
%  keyboard;
  RB = RB(:,1:(i-1));
  disp(['reduced to ',num2str(size(RB,2)),' basis vectors.']);    
%  keyboard;
  
  V = RB;
  W = ds_model_data.G * V;  
  save(fullfile(rbmatlabresult,rb_fname),'RB');
  
  ds_model.RB_basis_filename = fullfile(rbmatlabresult,rb_fname);
  ds_model.RB_generation_mode = 'file_load'; 
    
  detailed_data = gen_detailed_data(ds_model,ds_model_data);

  model = ds_model;
  save(fullfile(rbmatlabresult,detailed_data_fname),...
       'model','detailed_data')
 
  % compare detailed and reduced simulation
  load(fullfile(rbmatlabresult,detailed_data_fname),...
       'model','detailed_data')
  load(fullfile(rbmatlabresult,[trajectory_fnbase,'1']),...
       'ds_sim_data');

  ds_model = model;
  
  %perform single reduced simulation for given parameter
  reduced_data = gen_reduced_data(ds_model,detailed_data);
  ds_model.N = reduced_data.N;
  tic,
%  ds_model.nt = ds_model.nt/4;
%  disp('set nt to quarter!!');
  rb_sim_data = rb_simulation(ds_model,reduced_data);
  t_red = toc;
  rb_sim_data = rb_reconstruction(ds_model,detailed_data,rb_sim_data)
  disp(['time for reduced simulation: ',num2str(t_red)]);
  
  % plot of reduced trajectory
  lin_evol_model = ds_model.base_model;
  lin_evol_model_data = gen_model_data(lin_evol_model);  
  
  params.no_lines = 1;
  params.show_colorbar = 1;
  params.axis_equal = 1;
  params.plot = lin_evol_model.plot;
  % plot single snapshot
  %lin_evol_model.plot(rb_sim_data.X(:,1),lin_evol_model_data.grid, ...
  %                   params);
  
  % plot sequence:
  params.title = 'reduced solution';
  plot_sequence(rb_sim_data.X,lin_evol_model_data.grid,params)  

  % plot sequence:
  params.title = 'error';
  plot_sequence((rb_sim_data.X-ds_sim_data.X),...
                lin_evol_model_data.grid,params)  

  % plot of reduced output
  figure,plot(rb_sim_data.time,[rb_sim_data.Y;ds_sim_data.Y]); 
  legend('reduced output','detailed output');
 
 case 2 
  
  % generate more intelligent basis
  
  disp('initialization of lin_ds model from lin_evol model');
  %  verbose(1,'initialization of lin_ds model from lin_evol model and ...
  %       plot');
  
  model = fast_model(params);
  model.save_time_indices = 0:model.nt;
  model.theta = 0;
  
%  mu = [0.5,1,1];
  mu = [0.5,1];
  % set mu in model and base_model!!!
  model = model.set_mu(model,mu);
  model.RB_stop_Nmax = 150;

  %%%% basis generation: POD of single trajectory:
   model.RB_generation_mode = 'PCA_trajectory'; 

  %%% basis generation. POD of several trajectories
  %model.RB_generation_mode = 'PCA_trajectories'; 
  %model.RB_mu_list = {[0,0,0],[1,1,1]};
  
  %%% basis generation: greedy_random_fixed
  %model.RB_generation_mode = 'greedy_random_fixed'; 
  %model.RB_train_rand_seed = 12356;
  %model.RB_train_size = 5;
  %model.RB_error_indicator = 'estimator';
  %model.RB_stop_epsilon = 1e-2;
  %model.RB_stop_timeout = 1*60; % 1 minute
  %model.RB_stop_Nmax = 5;
  %model.RB_extension_algorithm = @RB_extension_PCA_fixspace;
  
  %%% basis generation: greedy_uniform_fixed
  %model.RB_generation_mode = 'greedy_uniform_fixed'; 
  %model.RB_numintervals = [1,1,1];
  %model.RB_error_indicator = 'estimator';
  %model.RB_extension_algorithm = @RB_extension_PCA_fixspace;
  %model.RB_stop_epsilon = 1e-2;
  %model.RB_stop_timeout = 1*60; % 1 minute
  %model.RB_stop_Nmax = 5;

  %%% basis generation: greedy_refined with uniform refinement
  %model.RB_generation_mode = 'greedy_refined'; 
  %model.RB_stop_max_val_train_ratio = 5;
  %model.RB_refinement_mode = 'uniform'; 
  %model.RB_val_rand_seed = 543;
  %model.RB_M_val_size = 10;
  %model.RB_numintervals = [1,1,1];
  %model.RB_max_refinement_level = 10;
  %model.RB_error_indicator = 'estimator';
  %model.RB_extension_algorithm = @RB_extension_PCA_fixspace;
  %model.RB_stop_epsilon = 1e-8;
  %model.RB_stop_timeout = 60*60; % 1 hour
  %model.RB_stop_max_val_train_ratio = 4; 
  %model.RB_stop_Nmax = 150;
  %model.mu_ranges = {[0.4,0.6],[0.4,0.6],[0.4,0.6]};
      
  %%% basis generation: greedy_refined with local refinement
  %% RB_detailed_val_savepath
  %model.RB_generation_mode = 'greedy_refined'; 
  %model.RB_refinement_mode = 'adaptive'; 
  %model.RB_stop_max_val_train_ratio = 5;
  %model.RB_refinement_theta = 0.2;
  %model.RB_val_rand_seed = 543;
  %model.RB_M_val_size = 10;
  %model.RB_numintervals = [1,1,1];
  %model.RB_max_refinement_level = 10;
  %model.RB_error_indicator = 'estimator';
  %model.RB_extension_algorithm = @RB_extension_PCA_fixspace;
  %model.RB_stop_epsilon = 1e-2;
  %model.RB_stop_timeout = 1*60; % 1 minute
  %model.RB_stop_Nmax = 5;

  % convert model into a p_part model
  %params.p_part_generation_mode = 'uniform';
  %params.p_part_numintervals = [2,1,1];
  %model = p_part_model(model,params);

  % convert model into a p_part model
  %params.p_part_generation_mode = 'adaptive';
  %params.p_part_numintervals = [2,1,1];
  %model = p_part_model(model,params);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  model_data = gen_model_data(model);
  detailed_data = gen_detailed_data(model,model_data);
  
  
  model = model;
  save(fullfile(rbmatlabresult,detailed_data_fname),...
       'model','detailed_data')
 
  % compare detailed and reduced simulation
  load(fullfile(rbmatlabresult,detailed_data_fname),...
       'model','detailed_data')
  load(fullfile(rbmatlabresult,[trajectory_fnbase,'1']),...
       'ds_sim_data');

  model = model;
  
  %perform single reduced simulation for given parameter
  reduced_data = gen_reduced_data(model,detailed_data);
  model.N = reduced_data.N;
  tic,
%  model.nt = model.nt/4;
%  disp('set nt to quarter!!');
  rb_sim_data = rb_simulation(model,reduced_data);
  t_red = toc;
  rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data)
  disp(['time for reduced simulation: ',num2str(t_red)]);
  
  % plot of reduced trajectory
  lin_evol_model = model.base_model;
  lin_evol_model_data = gen_model_data(lin_evol_model);  
  
  params.no_lines = 1;
  params.show_colorbar = 1;
  params.axis_equal = 1;
  params.plot = lin_evol_model.plot;
  % plot single snapshot
  %lin_evol_model.plot(rb_sim_data.X(:,1),lin_evol_model_data.grid, ...
  %                   params);
  
  % plot sequence:
  params.title = 'reduced solution';
  plot_sequence(rb_sim_data.X,lin_evol_model_data.grid,params)  

  % plot sequence:
  params.title = 'error';
  plot_sequence((rb_sim_data.X-ds_sim_data.X),...
                lin_evol_model_data.grid,params)  

  % plot of reduced output
  figure,plot(rb_sim_data.time,[rb_sim_data.Y;ds_sim_data.Y]); 
  legend('reduced output','detailed output');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% step = 7 % time measurements
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 case 7
  
%  clear all;
  nreps = 10;
  t_detailed = zeros(1,nreps);
  
  % initialize detailed model and simulate
  
  load(fullfile(rbmatlabresult,detailed_data_fname),...
       'model','detailed_data') 
  model = model;
    
%  lin_evol_model = advection_fv_output_vconst_model(params);
  %lin_evol_model.debug = 1;
  %lin_evol_model_data = gen_model_data(lin_evol_model);
  
  % set some additional fields in model which will be copied into
  % ds model depending on whether constants are known or not:
%  estimate_constants = 0;
%  if estimate_constants
%    %%% if constant estimaton has to be performed:
%    lin_evol_model.estimate_lin_ds_nmu = 4;
%    lin_evol_model.estimate_lin_ds_nX = 10;
%    lin_evol_model.estimate_lin_ds_nt = 4;
%  else
%    %%% otherwise:
%    % the following values are rough bounds, so could be 
%    % non-rigorous, then larger search must be performed.
%    lin_evol_model.state_bound_constant_C1 = 1;
%    lin_evol_model.output_bound_constant_C2 = 1.2916;
%  end;
%  
%  model = lin_ds_from_lin_evol_model(lin_evol_model);
%  % set accelerated data functions from below
%  model.A_function_ptr = @(model,model_data) ...
%      eval_affine_decomp(@A_components,...
%                        @A_coefficients,...
%                        model,model_data);
%  
%  model.B_function_ptr = @(model,model_data) ...
%      eval_affine_decomp(@B_components,...
%                        @B_coefficients,...
%                        model,model_data);
%  %model.u_function_ptr = @(params) 0.1; % define new
  model_data = gen_model_data(model);
  
  disp('skipping detailed simulations...')
%  if 0
    
  for r = 1:nreps
    disp(['repetition = ',num2str(r)])
    mu = [1,1,0.5]-rand(1,3)*0.01;
    model = model.set_mu(model,mu);
    tic,
    ds_sim_data = detailed_simulation(model,model_data);
    t_detailed(1,r) = toc;
  end;

%  end;
  
  disp('--------------------------------------------------')
  disp('detailed simulation')
  disp(t_detailed);
  disp(['mean t_detailed=',num2str(mean(t_detailed))]);
  disp('--------------------------------------------------')
  
  % reduced simulation with and without error estimation  
  
  clear('ds_sim_data');
  model.RB_generation_mode = 'file_load';
  model.RB_basis_filename = ...
      fullfile(rbmatlabresult,rb_fname);
  detailed_data = gen_detailed_data(model,model_data);
  reduced_data = gen_reduced_data(model,detailed_data);

  Ns = [150:-10:10]; %,30,20,10];
  Ns = Ns(find(reduced_data.N>=Ns));
    
  t_reduced_with_err_est = zeros(length(Ns),nreps);
  t_reduced_without_err_est = zeros(length(Ns),nreps);
  
  %   reduced simulation with and without error estimation
  
  for err_est = 0:1;
    model.error_estimation = err_est;
    disp(['err_est= ',num2str(err_est)]);
    for n_ind = 1:length(Ns)
      disp(['M = ',num2str(Ns(n_ind))]);
      model.N = Ns(n_ind);
      for r = 1:nreps
        disp(['repetition = ',num2str(r)])
        mu = [1,1,1]-rand(1,3)*0.01;
        model = model.set_mu(model,mu);
        tic;
        rb_sim_data = rb_simulation(model,reduced_data);
        if size(rb_sim_data.Xr,1)~=Ns(n_ind)
          disp('N dimension not set correctly!');
          keyboard;
        end;
        if err_est == 1
          t_reduced_with_err_est(n_ind,r) = toc;
        else
          t_reduced_without_err_est(n_ind,r) = toc;
        end;
      end;
    end;
  end;
    
  % output measurements

  disp('--------------------------------------------------')
  disp('detailed simulation')
  disp(t_detailed);
  disp(['mean t_detailed=',num2str(mean(t_detailed))]);

  disp('--------------------------------------------------')
  disp('reduced simulation with error estimation')
  disp(t_reduced_with_err_est);
  disp(['mean t_reduced_with_err_est=',...
        num2str(mean(t_reduced_with_err_est,2)')]);
  disp(['Ns =',num2str(Ns)]);
  
  disp('--------------------------------------------------')
  disp('reduced simulation without error estimation')
  disp(t_reduced_without_err_est);
  disp(['mean t_reduced_without_err_est=',...
        num2str(mean(t_reduced_without_err_est,2)')]);
  disp(['Ns =',num2str(Ns)]);
  disp('--------------------------------------------------')

  save(fullfile(rbmatlabresult,'time_measurements'));
  
%   dim_x = 65536, nt = 2048
% 
%mean t_detailed=64.3806
%--------------------------------------------------
%reduced simulation with error estimation
%    3.4542    3.5969    3.5659    3.4216    3.6170
%    3.3648    3.3816    3.4218    3.4391    3.4552
%    3.1643    3.1312    3.1244    3.1433    3.1720
%    3.0594    3.0496    3.0628    3.0486    3.0387
%    2.9993    2.9953    3.0078    2.9413    2.9603
%    2.9555    2.9093    2.9091    2.9510    2.9256%%%
%
%mean t_reduced_with_err_est=3.5311      3.4125       3.147      3.0518      2.9808      2.9301
%Ns =60  50  40  30  20  10
%--------------------------------------------------
%reduced simulation without error estimation
%    2.3975    2.4220    2.2924    2.4490    2.3302
%    2.3322    2.3198    2.3425    2.3056    2.3512
%    2.2227    2.2734    2.2210    2.2582    2.2504
%    2.2531    2.2180    2.2429    2.2212    2.2612
%    2.2283    2.1979    2.2581    2.2241    2.1957
%    2.2274    2.2014    2.2577    2.1856    2.2268
%
%mean t_reduced_without_err_est=2.3782      2.3303      2.2451      2.2393      2.2208      2.2198
%Ns =60  50  40  30  20  10
%--------------------------------------------------
   
 case 8 % compare original and accelerated detailed model
  
  disp('compare original and accelerated lin_ds models');
  
  %  clear all;
  % initialize detailed model and simulate  
  lin_evol_model = advection_fv_output_vconst_model(params);
  estimate_constants = 0;
  if estimate_constants
    %%% if constant estimaton has to be performed:
    lin_evol_model.estimate_lin_ds_nmu = 4;
    lin_evol_model.estimate_lin_ds_nX = 10;
    lin_evol_model.estimate_lin_ds_nt = 4;
  else
    %%% otherwise:
    % the following values are rough bounds, so could be 
    % non-rigorous, then larger search must be performed.
    lin_evol_model.state_bound_constant_C1 = 1;
    lin_evol_model.output_bound_constant_C2 = 1.2916;
  end;
  
  model = lin_ds_from_lin_evol_model(lin_evol_model);
  %model.u_function_ptr = @(params) 0.1; % define new
  model_data = gen_model_data(model);

  % set accelerated data functions
  model_fast = lin_ds_from_lin_evol_model(lin_evol_model);
  model_fast.A_function_ptr = @(model,model_data) ...
      eval_affine_decomp(@A_components,...
                         @A_coefficients,...
                         model,model_data);
  
  model_fast.B_function_ptr = @(model,model_data) ...
      eval_affine_decomp(@B_components,...
                         @B_coefficients,...
                         model,model_data);
  
 %  model_fast.C_function_ptr = @(model,model_data) ...
 %      eval_affine_decomp(@C_components,...
 %                       @C_coefficients,...
 %                       model,model_data);
  
  %  disp('skipping detailed simulations...')
  %  if 0
    
  mu = [1,1,0.5]-rand(1,3)*0.01;
  model = model.set_mu(model,mu);
  model_fast = model_fast.set_mu(model_fast,mu);
  tic,
  ds_sim_data = detailed_simulation(model,model_data);
  t_detailed = toc;
  tic,
  ds_sim_data_fast = detailed_simulation(model_fast,model_data);
  t_detailed_fast = toc;

  disp(['original t = ',num2str(t_detailed),...
        ', accelerated t = ',num2str(t_detailed_fast)]);
  if ~isequal(ds_sim_data,ds_sim_data_fast)
    disp('simulation results of original and accelerated model differ!')
  else
    disp('simulation results of original and accelerated model equal!')
  end;
  
  keyboard;
  
 case 9  % generate basis from POD-Greedy
  
  filecache_clear;
  model = fast_model(params);
  model.cone_range = [0.4,0.6];
  model.base_model.cone_range = model.cone_range; 
  %  model.
  % shrink range to 'single point'
%  r = [0.4,0.6]; % => 2 -> 0.8 error nach 10 basisvektoren
%  r = [0.0,1.0]; % => 4.6 -> 1.9 error nach 10 basisvektoren
%  r = [0.45,0.55]; % => 2.38 -> 0.51 error nach 10 basisvektoren
  model.mu_ranges = {[0,1],[0.4,0.6],[0.4,0.6]};
  %  model.rb_basis_generation_ptr = ...
  %      @basisgen_POD_greedy_uniform_fixed;
  
  %  model.RB_generation_mode = 'uniform_fixed';
  %  model.RB_num_intervals = [2,2,2];
  
  model_data = gen_model_data(model);
    
  disp('generating basis from POD-Greedy');
  
  % for basis generation, set further fields:

  model.verbose = 5; % => only dots in greedy loop
  model.RB_generation_mode = 'greedy_uniform_fixed';
  model.RB_numintervals = [1,1,1]; % => 0.5 in cone_pos obtained  
%  model.RB_numintervals = [4,4,4];  
  model.RB_stop_epsilon = 1e-5;
  model.RB_stop_timeout = 3600;
  model.RB_stop_Nmax =150;
%  model.RB_stop_timeout = 600;
%  model.RB_stop_Nmax =10;
  model.RB_extension_algorithm = @RB_extension_PCA_fixspace;
  model.RB_error_indicator = 'estimator';
  % default:
  % model.get_estimator_from_sim_data = @(sim_data) sim_data.DeltaX(end);

  detailed_data = gen_detailed_data(model,model_data);
  
  model = model;
  dfname = fullfile(rbmatlabresult,[detailed_data_fname,'_greedy']);
  save(dfname,'model','detailed_data');
  
  disp('successfully generated and saved detailed_data');

  params.plot_title = 'reduced basis';
  params.plot = model.base_model.plot;
  lin_evol_model_data = gen_model_data(model.base_model);
  plot_sequence(detailed_data.V,lin_evol_model_data.grid,params)  

  keyboard;

 case 10  % reduced simulations and error plots

  %basis_from_step = 9; % step 9
  basis_from_step = 4; % step 4

  if basis_from_step == 9
    dfname = fullfile(rbmatlabresult,[detailed_data_fname,'_greedy']);
    load(dfname);
    disp('loaded basis from step 9')
  else
    load(fullfile(rbmatlabresult,detailed_data_fname),...
         'model','detailed_data') 
    
%    model = fast_model(params);
    model_data = gen_model_data(model);
    %  dfname = fullfile(rbmatlabresult,'advection_fv_output_detailed_data_146');
%    model.RB_basis_filename = ...
%       fullfile(rbmatlabresult,'advection_fv_output_basis');
%    model.RB_generation_mode = 'file_load';
%    detailed_data = gen_detailed_data(model,model_data);  
    %  dfname = fullfile(rbmatlabresult,'advection_fv_output_detailed');
    % load(dfname);
    tmp = load(fullfile(rbmatlabresult,[trajectory_fnbase,'2']));
    ds_sim_data = tmp.ds_sim_data;
    mu_from_file = get_mu(tmp.model);
    disp('loaded basis from step 4')
    model = model.set_mu(model,mu_from_file);
  end;
  % trajectory2 above is made with this 
  %   model = model.set_mu(model,[1,1,1]);
  %  model = model.set_mu(model,[1,1,1]);
  % trajectory3 above is made with this: 
  %model = model.set_mu(model,[1,1,0.5]);
  %model = model.set_mu(model,[0.5,0.6,0.6]);
  %model = model.set_mu(model,[0.5,0.4,0.4]);
  %model = model.set_mu(model,[0.5,0.5,0.5]);
  model.N = size(detailed_data.V,2);
  %model.N = 150;
  reduced_data = gen_reduced_data(model,detailed_data);
  rb_sim_data = rb_simulation(model,reduced_data);
  
  params.no_lines = 1;
  params.show_colorbar = 1;
  params.axis_equal = 1;
  
  % plot of reduced output
  figure;
  p = plot(rb_sim_data.time,...
           [rb_sim_data.Y;...
            rb_sim_data.Y+rb_sim_data.DeltaY; ...
            rb_sim_data.Y-rb_sim_data.DeltaY]);
  colors = {[1,0,0],[1,0,0],[1,0,0],[0,1,0]};
  linestyles = {'-','-.','-.','-'};
  widths = [2,2,2,2];
  for i = 1:length(p)
    set(p(i),'Color',colors{i},'Linestyle',linestyles{i},...
             'Linewidth',widths(i));
  end;
  % set marker on rb_sim_data:
  Xdata = get(p(1),'XData');
  Ydata = get(p(1),'YData');
  Zdata = get(p(1),'ZData');
  set(p(1),'XData',Xdata(1:4:end));
  set(p(1),'YData',Ydata(1:4:end));
  set(p(1),'ZData',Zdata(1:4:end));
  set(p(1),'Marker','o');
  
  legend(['y\^(t): reduced output'],...
         'y\^(t)+\Delta_y(t)',...
         'y\^(t)-\Delta_y(t)','Location','NorthWest');
  title(['k = ',num2str(model.N)],'Fontsize',15);

  % table of error, error estimator and effectivity

  reduced_data = gen_reduced_data(model,detailed_data);
  
  Ns = 150:-10:10;
  Ns = Ns(find(reduced_data.N>=Ns));
  
  err = zeros(1,length(Ns));
  err_estim = zeros(1,length(Ns));
  effectivity = zeros(1,length(Ns));
  
  for ni = 1:length(Ns)
    model.N = Ns(ni);
    rb_sim_data = rb_simulation(model, reduced_data);
    err(ni) = abs(ds_sim_data.Y(end) - rb_sim_data.Y(end));
    err_estim(ni) = rb_sim_data.DeltaY(end);
    effectivity(ni) = err_estim(ni)/err(ni);
  end;
  Ns
  err
  err_estim
  effectivity
  disp('finished computing effectivities')
  keyboard;
    
 case 11 % test of matlab ode solver for resulting sparse system
 
  params.coarse_factor = 8; % nake model smaller
  model = fast_model(params);
  model.RB_generation_mode = 'from_detailed_data'; 
  
  % detailed simulation:
  model_data = gen_model_data(model);
  model_data.RB = zeros(size(model_data.G,1),0);
  detailed_data = gen_detailed_data(model,model_data);  
  
  model.decomp_mode = 0; % complete
  model = model.set_time(model,0)
  model = model.set_mu(model,[0.5,1,1]);
  
  A = @(model,model_data,t) ...
      model.A_function_ptr(model.set_time(model,t),model_data);
  
  B = @(model,model_data,t) ...
      model.B_function_ptr(model.set_time(model,t),model_data); 
  
  x0 = model.x0_function_ptr(model,model_data);
  tic;
  [t,x] = ode15s(@(t,x) A(model,model_data,t)*x + B(model,model_data,t), ...
                 [0 1], x0);
  time_ode15 = toc;
  disp(['ode15 takes ',num2str(time_ode15),' sec']);
  
  params.plot_title = 'ode15 solution';
  params.plot = model.base_model.plot;
  params.axis_equal = 1;
  lin_evol_model_data = gen_model_data(model.base_model);
  plot_sequence((x'),lin_evol_model_data.grid,params)  

  tic;
  sim_data = detailed_simulation(model,model_data);
  time_rbmatlab = toc;
  disp(['rbmatlab takes ',num2str(time_rbmatlab),' sec']);
  params.plot_title = 'rbmatlab solution';
  params.plot = model.base_model.plot;
  params.axis_equal = 1;
  lin_evol_model_data = gen_model_data(model.base_model);
  plot_sequence(sim_data.X,lin_evol_model_data.grid,params)  

  % further plots for MCMCS paper
 case 12
  
  load(fullfile(rbmatlabresult,detailed_data_fname),...
       'model','detailed_data') 
  
  % plot of error estimator over parameter variation
%  model = fast_model(params);
%  model_data = gen_model_data(model);
  %  dfname = fullfile(rbmatlabresult,'advection_fv_output_detailed_data_146');
%  model.RB_basis_filename = ...
%      fullfile(rbmatlabresult,'advection_fv_output_basis');
%  model.RB_generation_mode = 'file_load';
%  detailed_data = gen_detailed_data(model,model_data);  
%  save(fullfile(rbmatlabresult,...
%               'advection_fv_output_detailed_2trajectories.mat'),...
%       'model','detailed_data')
%  keyboard;
  %  dfname = fullfile(rbmatlabresult,'advection_fv_output_detailed');
  % load(dfname);
  tmp = load(fullfile(rbmatlabresult,[trajectory_fnbase,'2']));
  ds_sim_data = tmp.ds_sim_data;
  disp('loaded basis from step 4')
  
  %  factors = 0:0.05:1;
  reduced_data = gen_reduced_data(model,detailed_data);
  Ns = [1,8,16,32,64];
  Ns = Ns(find(reduced_data.N>=Ns));
  factors = 0:0.05:1;
  err_estim = zeros(length(factors),length(Ns));
  for ni = 1:length(Ns)
    model.N = Ns(ni);
    for i =1:length(factors)
      model = model.set_mu(model, factors(i)*[1,1,1]);
      rb_sim_data = rb_simulation(model, reduced_data);
      %    err(ni) = abs(ds_sim_data.Y(end) - rb_sim_data.Y(end));
      err_estim(i,ni) = rb_sim_data.DeltaY(end);
      %    effectivity(ni) = err_estim(ni)/err(ni);   
    end;
  end;
%  keyboard;
  p = plot(factors,err_estim');
  colors = {[1,0,0],[0,0.6,0],[0,0,1],[0,0,0]};
  linestyles = {'-','-.','-','-.'};
  for i = 1:4;
    set(p(i),'Linewidth',2,'Color',colors{i},'Linestyle',linestyles{i});
  end;  
  
  xlabel('parameter \mu_1=\mu_2=\mu_3');
  ylabel('\Delta_y(\mu)');  
  title('error estimator from parameter sweep')
  legend({'k=1','k=8','k=16','k=32','k=64'},'Location','NorthWest');
  saveas(gcf,fullfile(rbmatlabresult,'parameter_sweep.eps'),'epsc');
  keyboard;
 otherwise
  disp('step number unknown');
end;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% generation of model from lin_evol_model and replace 
% by fast vectors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function model = fast_model(params);
lin_evol_model = advection_fv_output_vconst_model(params);
lin_evol_model.vy_weight = 0.5;
lin_evol_model.mu_names = {'cone_amplitude','vx_weight'};
lin_evol_model.mu_ranges = {[0,1],[0,1]};
%lin_evol_model.debug = 1;
%lin_evol_model_data = gen_model_data(lin_evol_model);

% set some additional fields in model which will be copied into
% ds model depending on whether constants are known or not:
estimate_constants = 0;
if estimate_constants
  %%% if constant estimaton has to be performed:
  lin_evol_model.estimate_lin_ds_nmu = 4;
  lin_evol_model.estimate_lin_ds_nX = 10;
  lin_evol_model.estimate_lin_ds_nt = 4;
else
  %%% otherwise:
  % the following values are rough bounds, so could be 
  % non-rigorous, then larger search must be performed.
  lin_evol_model.state_bound_constant_C1 = 1;
  lin_evol_model.output_bound_constant_C2 = 1.2916;
end;

model = lin_ds_from_lin_evol_model(lin_evol_model);
% set accelerated data functions from below
model.A_function_ptr = @(model,model_data) ...
    eval_affine_decomp(@A_components,...
                       @A_coefficients,...
                       model,model_data);

model.B_function_ptr = @(model,model_data) ...
    eval_affine_decomp(@B_components,...
                       @B_coefficients,...
                       model,model_data);
%model.u_function_ptr = @(params) 0.1; % define new

%model   .theta = 1;
%disp('set theta to 1!!!');
%keyboard;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% auxialiary coefficient functions for acceleration of 
% advection model: explicit implementation of coefficients
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Acomp = A_components(model,model_data) 
Acomp = model_data.A_components;

function Acoeff = A_coefficients(model); 
%model.base_model.decomp_mode = 2;
%model.base_model.t = model.t;
%mu = get_mu(model);
%model.base_model = model.set_mu(model.base_model,mu);
%[L_I_coeff, L_E_coeff, b_coeff] = ...
%    model.base_model.operators_ptr(model.base_model, );
%% L_E = Id + Delta t * A 
%Acoeff = L_E_coeff(2:end)/model.base_model.dt;
%keyboard;
%Acoeff = -model.base_model.dt * ...
%        model.base_model.velocity_coefficients_ptr(model.base_model);
%if model.t>0
  Acoeff = - [model.vx_weight, model.base_model.vy_weight]*0.5; %*(1-model.t);
%  keyboard;
%end;
  
function Bcomp = B_components(model,model_data) 
Bcomp = model_data.B_components;

function Bcoeff = B_coefficients(model); 
% hmmm the coefficients are now called twice in A and B
% should somehow be cached later...
%model.base_model.decomp_mode = 2;
%model.base_model.t = model.t;
%mu = get_mu(model);
%model.base_model = model.set_mu(model.base_model,mu);
%[L_I_coeff, L_E_coeff, b_coeff] = ...
%    model.base_model.operators_ptr(model.base_model, );
%% L_E = Id + Delta t * A 
%Bcoeff2 = b_coeff/model.base_model.dt;
vcoeff = - [model.vx_weight, model.base_model.vy_weight]*0.5;%*(1-model.t);
Q_0 = model.base_model.cone_number;
dcoeff = zeros(1,Q_0);
max_pos = model.base_model.cone_weight * (Q_0-1)+1;
t = model.t;
for q = 1:Q_0  
  dcoeff(q) = (1-(max_pos-q))*(1-t) * ((max_pos>=q) && (max_pos < q+1)) ...
      + (1+(max_pos-q))*(1-t) * ((max_pos>=q-1) && (max_pos < q));
end;
v = -(vcoeff(:)*(dcoeff(:)'));
Bcoeff =  v(:)*model.cone_amplitude;
%if model.t>0
%keyboard;
%end;

%function Ccomp = C_components(model,model_data) 
%Ccomp = model_data.C_components;
%
%function Ccoeff = C_coefficients(model); 
%model.base_model.decomp_mode = 2;
%model.base_model.t = model.t;
%mu = get_mu(model);
%model.base_model = model.set_mu(model.base_model,mu);
%Ccoeff = ...
%    model.base_model.operators_output(model.base_model);


%| \docupdate