advection_fv_output_mcmds.m

Go to the documentation of this file.

function advection_fv_output_mcmds(step)
% small script implementing a simple advection example
% for producing matrices to be used in the RB-DS framework
% Discretization with FV Functions. Used for the mcmds lin-ds paper.
%
% mu_names  = {'cone_weight','vx_weight','vy_weight'};
% mu_ranges = [0,1]^3
%
% possible steps:
%
% step = 1; % initialization of model and plot of model data
% step = 2; % detailed simulation of lin_evol
% step = 2.5; % conversion to DS primal model, lin_ds detailed simulation 
% step = 3; % conversion to DS primal model, DS detailed simulation 
%             and comparison with lin_evol
% step = 4 % compute several trajectories.
% step = 4.5 % compute reduced basis based on 2 trajecories.
% step = 4.75 % comparison of detailed and reduced simulation
% step = 5 % figures of detailed simulations based on step 4
% step = 6 % figures of reduced simulations based on step 4
% step = 7 % time measurements based on step 4
% step = 8; % compare original and accelerated ds_model
% step = 9; % generate greedy basis on subset of parameter domain
% step = 10;  % reduced simulations and error plots
% step = 11; % test of matlab ode solver for resulting sparse system
% step = 12; % further plots for MCMDS paper

% Bernard Haasdonk 25.8.2009

% cartesian grid

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 = 1; 
%verbose(1);

switch step
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %step = 1; % initialization of model and plot of model data
 case 1
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
%  verbose(1,'initializaton of lin_evol_model and plot model_data');
  disp('initializaton of lin_evol_model and plot model_data');
  model = advection_fv_output_model(params);
  model_data = gen_model_data(model);
  params.axis_equal = 1;
  params.axis_tight = 1;
  plot(model_data.grid,params);
  %    inspect(model_data.grid)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %step = 2; % detailed simulation of lin_evol
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 case 2
%  verbose(1,'detailed simulation of lin_evol_model and plot');
  disp('detailed simulation of lin_evol_model and plot');
  model = advection_fv_output_model(params);
  model_data = gen_model_data(model);
  %    model.debug = 1;
  %    disp('model debugging turned on.')
  model.cone_weight = 0.0;
  model.vx_weight = 0.0;
  model.vy_weight = 0.0;
  model.vx_weight = 1;
  model.vy_weight = 1;
  model.cone_weight = 1.0;
%  model.vx_weight = 1;
%  model.vy_weight = 0;
  sim_data = detailed_simulation(model,model_data);
  plot_sim_data(model,model_data,sim_data,[]);

  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%step = 2.5; % conversion to DS primal model, DS detailed simulation 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 case 2.5
  disp('initialization of lin_ds model from lin_evol model');
  %  verbose(1,'initialization of lin_ds model from lin_evol model and ...
  %       plot');
  
  lin_evol_model = advection_fv_output_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;
  
  ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
  %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];
  % 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);

  params.no_lines = 1;
  params.clim = [0,1];
  params.plot_title='detailed implicit lin_ds';
  lin_evol_from_ds_sim_data.U = ds_sim_data.X;
  lin_evol_from_ds_sim_data.y = ds_sim_data.Y;
  plot_sim_data(lin_evol_model,lin_evol_model_data,lin_evol_from_ds_sim_data,params);
    
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%step = 3; % conversion to DS primal model, DS detailed simulation 
        % and comparison with lin_evol
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 case 3 % initialization of DS model and lin_ds detailed simulation
  disp('initialization of lin_ds model from lin_evol model and plot');
%  verbose(1,'initialization of lin_ds model from lin_evol model and ...
%         plot');
  
  lin_evol_model = advection_fv_output_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;
  
  ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
%  ds_model.theta = 0.0;
  %ds_model.u_function_ptr = @(params) 0.1; % define new
  
  ds_model_data = gen_model_data(ds_model);
  mu = [0.5,1,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);
  lin_evol_model = lin_evol_model.set_mu(lin_evol_model,mu);
  lin_evol_sim_data = detailed_simulation(lin_evol_model,...
                                          lin_evol_model_data);
  err = lin_evol_sim_data.U - ds_sim_data.X;
  disp(['max error in DOFs of lin_evol and ds simulation: ',...
       num2str(max(abs(err(:))))]);

  params.plot_title='detailed lin_evol';
  params.no_lines = 1;
  plot_sim_data(lin_evol_model,lin_evol_model_data,lin_evol_sim_data,params);
  params.plot_title='detailed lin_ds';
  lin_evol_from_ds_sim_data.U = ds_sim_data.X;
  lin_evol_from_ds_sim_data.y = ds_sim_data.Y;
  plot_sim_data(lin_evol_model,lin_evol_model_data,lin_evol_from_ds_sim_data,params);
  params.plot_title='error explicit lin_ds and lin_evol';
  error_sim_data = lin_evol_sim_data;
  error_sim_data.U = error_sim_data.U-ds_sim_data.X;
  plot_sim_data(lin_evol_model,lin_evol_model_data,error_sim_data,params);
  keyboard;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %step = 4 % reduced Basis generation by single trajectories and
          %comparison of detailed and reduced simulation
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
 case 4 % generate some trajectories

  ds_model = fast_model(params); 
%  lin_evol_model = advection_fv_output_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;
%  
%  ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
  %ds_model.u_function_ptr = @(params) 0.1; % define new
  ds_model_data = gen_model_data(ds_model);
%  mu = [1,1,0.5];
%  mu = [0.5,0.5,0.5];

  mus = {[0,0,0],[1,1,1],[1,1,0.5],[0,0.5,1]};
  for m = 1:length(mus)
    disp(['computing trajectory number ',num2str(m)]);
    mu = mus{m};
    ds_model = ds_model.set_mu(ds_model,mu);
    %ds_model.theta = 0;
    %ds_model.theta = 1;
    %disp('set theta to 1!!!');
    tic;
    ds_sim_data = detailed_simulation(ds_model,ds_model_data);
    t_det = toc;
    disp(['time for full simulation: ',num2str(t_det)]);
    save(fullfile(rbmatlabresult,['trajectory',num2str(m)]),...
         'ds_model','ds_model_data','ds_sim_data');
  end;
  
 case 4.5    % compute reduced basis.
  % make sure that the previous trajectories are computed!!
  
  disp('PCA of snapshots takes a while...');
  load(fullfile(rbmatlabresult,'trajectory1'),...
       'ds_model','ds_model_data','ds_sim_data');
  epsilon = 1e-4;
  RB1 = PCA_fixspace(ds_sim_data.X,[],ds_model_data.G,1,[], ...
                     epsilon);
  %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;
  % select first vector only (stationary solution)
  RB = RB1(:,1);
  disp(['reduced to ',num2str(size(RB,2)),' basis vectors.']);    
%  keyboard;
  
  load(fullfile(rbmatlabresult,'trajectory2'),...
       'ds_model','ds_model_data','ds_sim_data');

  epsilon = 1e-12;
%  ds_sim_data.X = ds_sim_data.X-...
%      RB*(RB'*(ds_model_data.G*ds_sim_data.X));
  RB2 = PCA_fixspace(ds_sim_data.X,RB,ds_model_data.G,[],'gram-schmidt',epsilon);
  %RB = orthonormalize_qr(ds_sim_data.X,ds_model_data.G);
  % select correct orthogonal set
  K = RB2' * ds_model_data.G * RB2;
  i = find(abs(diag(K)-1)>1e-2);
  i = min(i);
  if isempty(i)
    i = size(RB2,2)+1;
  end;
  disp(['found ',num2str(size(RB2,2)),' basis vectors after traj. 2']);
  RB2 = RB2(:,1:(i-1));
  disp(['reduced to ',num2str(size(RB2,2)),' basis vectors after traj. 2']);
  RB = [RB, RB2];
%  RB = orthonormalize(RB);
  keyboard;

% check orthonormality
  K = RB' * ds_model_data.G * RB;
  i = find(abs(diag(K)-1)>1e-2);
  i = min(i);
  if isempty(i)
    i = size(RB,2)+1;
  end;
  % make consistent with preprint: 127 vectors overall
  if i>128
    i = 128;
  end;
  
  RB = RB(:,1:(i-1));
  disp(['reduced to ',num2str(size(RB,2)),' basis vectors overall']);
  keyboard;  
  
  V = RB;
  W = ds_model_data.G * V;  
  rbfname = fullfile(rbmatlabresult,'advection_fv_output_basis');       
  save(rbfname,'RB');
  
  ds_model.RB_basis_filename = rbfname;
  ds_model.RB_generation_mode = 'file_load'; 
    
  detailed_data = gen_detailed_data(ds_model,ds_model_data);

  model = ds_model;
  save(fullfile(rbmatlabresult,'advection_fv_output_detailed_data_2traj'),...
       'model','detailed_data')
  
 case 4.75 % compare reduced and detailed simulation
  
  load(fullfile(rbmatlabresult,'advection_fv_output_detailed_data_2traj'),...
       'model','detailed_data')
  load(fullfile(rbmatlabresult,'trajectory2'),...
       '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');
   
  keyboard;
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 case 5 % generate plots of detailed data for paper: 
 % step = 5 % figures of detailed simuations based on step 4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  disp('generate plots of detailed_data for morepas poster');

  load(fullfile(rbmatlabresult,'advection_fv_output_detailed_data_2traj'),...
       'model','detailed_data')
  
%  model = fast_model(params); 
%  rbfname = 'advection_fv_output_basis';       
%  load(fullfile(rbmatlabresult,rbfname),'ds_model','ds_model_data');

  % make sure, that these trajectories exist.
  trajectories = [3,4];

  for tr = 1:length(trajectories);
    trname = ['trajectory',num2str(trajectories(tr))];
    % plots of trajectory  
    load(fullfile(rbmatlabresult,trname));
    lin_evol_model = ds_model.base_model;
    lin_evol_model_data = gen_model_data(lin_evol_model);  
    
    params.no_lines = 1;
    params.show_colorbar = 0;
    params.axis_equal = 1;
    params.plot = lin_evol_model.plot;
    lin_evol_model.plot(ds_sim_data.X(:,1),lin_evol_model_data.grid, ...
                        params);
    set(gca,'Xtick',[])
    set(gca,'Ytick',[])
    saveas(gcf,fullfile(rbmatlabresult,[trname,'_t1.eps']),'epsc');
    close(gcf);
    
    lin_evol_model.plot(ds_sim_data.X(:,64),lin_evol_model_data.grid, ...
                        params);
    set(gca,'Xtick',[])
    set(gca,'Ytick',[])
    saveas(gcf,fullfile(rbmatlabresult,[trname,'_tmiddle.eps']),'epsc');
    close(gcf);
    
    lin_evol_model.plot(ds_sim_data.X(:,end),lin_evol_model_data.grid, ...
                        params);
    set(gca,'Xtick',[])
    set(gca,'Ytick',[])
    saveas(gcf,fullfile(rbmatlabresult,[trname,'_tend.eps']),'epsc');
    close(gcf);
    
    p = ds_model.plot_sim_data(ds_model, ds_model_data, ds_sim_data,params);
    % Ylim = get(gca,'Ylim');
    % Ylim = [Ylim(1),Ylim(2)*1.1];
    set(gca,'Ylim',[0,0.09]);
    set(p(2),'Color',[0,0.8,0]);
    saveas(gcf,fullfile(rbmatlabresult,[trname,'_output.eps']),'epsc');
    close(gcf);
    close(gcf);    
  end;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% step = 6 % figures of reduced simuations based on step 4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 case 6 % generate plots of reduced model for poster/paper 
  
  disp('generate plots of reduced simulations output estimation');

  load(fullfile(rbmatlabresult,'advection_fv_output_detailed_data_2traj'),...
       'model','detailed_data') 
  %  ds_model = model;
  
  % as reference
  %  load(fullfile(rbmatlabresult,'trajectory1'));
  load(fullfile(rbmatlabresult,'trajectory2'));
  ds_model.enable_error_estimator = 1;
  
%  ds_model.error_estimation = 1;
%  ds_model.RB_generation_mode = 'file_load'
%  ds_model.RB_basis_filename = fullfile(...
%      rbmatlabresult,'advection_fv_output_basis');
%  detailed_data = gen_detailed_data(ds_model,ds_model_data);
  reduced_data = gen_reduced_data(ds_model,detailed_data);
  
  Ns = [150:-10:10];
  
  Ns = Ns(find(reduced_data.N>=Ns));
  
  for nind = 1:length(Ns)
    % perform reduced simulation
    ds_model.N = Ns(nind);
    rb_sim_data = rb_simulation(ds_model,reduced_data);
    rb_sim_data = rb_reconstruction(ds_model,detailed_data,rb_sim_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; ...
                 ds_sim_data.Y]);
    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');
    
    l = legend(['y\^(t): reduced output'],...
           'y\^(t)+\Delta_y(t)',...
           'y\^(t)-\Delta_y(t)',...
           ['y(t): detailed output'],'Location','SouthEast');
    set(l,'Fontsize',15);
    title(['k = ',num2str(ds_model.N)],'Fontsize',15);
    saveas(gcf,...
           fullfile(rbmatlabresult,...
                    ['output_estimation',num2str(ds_model.N),'.eps']),'epsc');
    close(gcf);
  end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% step = 7 % time measurements
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 case 7
  disp('time measurements for morepas poster');
  
%  clear all;
  nreps = 10;
  t_detailed = zeros(1,nreps);
  
  % initialize detailed model and simulate
  
  load(fullfile(rbmatlabresult,'advection_fv_output_detailed_data_2traj'),...
       'model','detailed_data') 
  ds_model = model;
    
%  lin_evol_model = advection_fv_output_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;
%  
%  ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
%  % set accelerated data functions from below
%  ds_model.A_function_ptr = @(model,model_data) ...
%      eval_affine_decomp(@A_components,...
%                        @A_coefficients,...
%                        model,model_data);
%  
%  ds_model.B_function_ptr = @(model,model_data) ...
%      eval_affine_decomp(@B_components,...
%                        @B_coefficients,...
%                        model,model_data);
%  %ds_model.u_function_ptr = @(params) 0.1; % define new
  ds_model_data = gen_model_data(ds_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;
    ds_model = ds_model.set_mu(ds_model,mu);
    tic,
    ds_sim_data = detailed_simulation(ds_model,ds_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');
  ds_model.RB_generation_mode = 'file_load';
  ds_model.RB_basis_filename = ...
      fullfile(rbmatlabresult,'advection_fv_output_basis');
  detailed_data = gen_detailed_data(ds_model,ds_model_data);
  ds_model.enable_error_estimator = 1;
  reduced_data = gen_reduced_data(ds_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;
%    ds_model.error_estimation = err_est;
    ds_model.enable_error_estimator = err_est;
    disp(['err_est= ',num2str(err_est)]);
    for n_ind = 1:length(Ns)
      disp(['M = ',num2str(Ns(n_ind))]);
      ds_model.N = Ns(n_ind);
      for r = 1:nreps
        disp(['repetition = ',num2str(r)])
        mu = [1,1,1]-rand(1,3)*0.01;
        ds_model = ds_model.set_mu(ds_model,mu);
        tic;
        rb_sim_data = rb_simulation(ds_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_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;
  
  ds_model = lin_ds_from_lin_evol_model(lin_evol_model);
  %ds_model.u_function_ptr = @(params) 0.1; % define new
  ds_model_data = gen_model_data(ds_model);

  % set accelerated data functions
  ds_model_fast = lin_ds_from_lin_evol_model(lin_evol_model);
  ds_model_fast.A_function_ptr = @(model,model_data) ...
      eval_affine_decomp(@A_components,...
                         @A_coefficients,...
                         model,model_data);
  
  ds_model_fast.B_function_ptr = @(model,model_data) ...
      eval_affine_decomp(@B_components,...
                         @B_coefficients,...
                         model,model_data);
  
 %  ds_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;
  ds_model = ds_model.set_mu(ds_model,mu);
  ds_model_fast = ds_model_fast.set_mu(ds_model_fast,mu);
  tic,
  ds_sim_data = detailed_simulation(ds_model,ds_model_data);
  t_detailed = toc;
  tic,
  ds_sim_data_fast = detailed_simulation(ds_model_fast,ds_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;
  ds_model = fast_model(params);
  ds_model.cone_range = [0.4,0.6];
  ds_model.base_model.cone_range = ds_model.cone_range; 
  %  ds_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
  ds_model.mu_ranges = {[0,1],[0.4,0.6],[0.4,0.6]};
  %  ds_model.rb_basis_generation_ptr = ...
  %      @basisgen_POD_greedy_uniform_fixed;
  
  %  ds_model.RB_generation_mode = 'uniform_fixed';
  %  ds_model.RB_num_intervals = [2,2,2];
  
  ds_model_data = gen_model_data(ds_model);
    
  disp('generating basis from POD-Greedy');
  
  % for basis generation, set further fields:

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

  detailed_data = gen_detailed_data(ds_model,ds_model_data);
  
  model = ds_model;
  dfname = fullfile(rbmatlabresult,'advection_fv_output_lin_ds_detailed');
  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,'advection_fv_output_lin_ds_detailed');
    load(dfname);
    disp('loaded basis from step 9')
  else
    load(fullfile(rbmatlabresult,'advection_fv_output_detailed_data_2traj'),...
         '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,'trajectory2'));
    ds_sim_data = tmp.ds_sim_data;
    mu_from_file = get_mu(tmp.ds_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');
  
  l = legend(['y\^(t): reduced output'],...
         'y\^(t)+\Delta_y(t)',...
         'y\^(t)-\Delta_y(t)','Location','NorthWest');
  set(l,'Fontsize',15);
  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; % make 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,'advection_fv_output_detailed_data_2traj'),...
       'model','detailed_data') 
  model.enable_error_estimator = 1;
  % 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,'trajectory2'));
  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,4,8,16];
  Ns = Ns(find(reduced_data.N>=Ns));
  factors = 0:0.05:1;
  err_estim = zeros(length(factors),length(Ns));

%  if 0
  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;
  
%  end; % if 0
  
  
%  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','Fontsize',15);
  ylabel('\Delta_y(\mu)','Fontsize',15);  
  title('error estimator from parameter sweep','Fontsize',15)
  l = legend({'k=1','k=4','k=8','k=16'},'Location','NorthWest');
  set(l,'Fontsize',15);
  saveas(gcf,fullfile(rbmatlabresult,'parameter_sweep.eps'),'epsc');
  % plot of true errors
  
  model_data = gen_model_data(model);
  
  figure;
%  factors = [0,0.5,1];
  errs = zeros(length(factors),length(Ns));
  for i =1:length(factors)
    model = model.set_mu(model, factors(i)*[1,1,1]);
    ds_sim_data = filecache_function(@detailed_simulation,model,model_data);
    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));
      e = rb_sim_data.Y-ds_sim_data.Y;
      e = max(sqrt(max(sum(e.*e),0)));
      errs(i,ni) = e;
      %    effectivity(ni) = err_estim(ni)/err(ni);   
    end;
  end;
  %  keyboard;
  p = plot(factors,errs');
  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;  
  % same range as first plot
  set(gca,'Ylim',[0,10]);
  
  xlabel('parameter \mu_1=\mu_2=\mu_3','Fontsize',15);
  ylabel('|y(\mu)-y\^(\mu)|','Fontsize',15);  
  title('output errors from parameter sweep','Fontsize',15)
  l = legend({'k=1','k=4','k=8','k=16'},'Location','NorthWest');
  set(l,'Fontsize',15);
  saveas(gcf,fullfile(rbmatlabresult,'parameter_sweep_errs.eps'),'epsc');
  keyboard;
  % plot of true errors
  
 otherwise
  disp('step number unknown');
end;

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

function ds_model = fast_model(params);
  % fast_model
lin_evol_model = advection_fv_output_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;

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

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

%ds_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) 
  % A_components
Acomp = model_data.A_components;

function Acoeff = A_coefficients(model); 
  % A_coefficients

%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.vy_weight]*(1-model.t);
%  keyboard;
%end;
  
function Bcomp = B_components(model,model_data) 
  % B_components
Bcomp = model_data.B_components;

function Bcoeff = B_coefficients(model); 
  % B_coefficients

% 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.vy_weight]*(1-model.t);
Q_0 = model.base_model.cone_number;
dcoeff = zeros(1,Q_0);
max_pos = 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(:);
%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);