adaptive_basisgen_mcmds_changed.m

function adaptive_basisgen_mcmds(step, loadfilename)
%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
%
%   Options for the experiments (like 2-D, 3-D case , etc) can be fixed
%   below.
%
% mu_names  = {'cone_weight','vx_weight','vy_weight'};
% mu_ranges = [0,1]^3
%
% steps:
%
% step 1: detailed simulation of single trajectory, plot
%         POD and save reduced basis
%         reduced simulation and error
% step 2: conduct a basis generation using the RB-method defined in the
%         options. 
%         detailed_data and model are saved in a mat-file
% step 3: the previously saved mat-file is loaded and an animation/a plot
%         of the basisgeneration and p-partition is made
% step 31: detailed and reduced simulation of previously computed basis
%
%     Experiments for the MCMDS paper:
% step 4: Basisgeneration using NO p-partition and a fixed-greedy
%         algorithm.
%         M-Train for greedy can be fixed, but will probably be a 3^p grid
%         detailed_data and model saved in mat file
%
% step 41: mat file from previous step is loaded. 
%          plot of grid?
%
% step 5: Basisgeneration using a fixed p-partition and a fixed-greedy
%         algorithm using the same M-Train grid as the previous step
%         No adaption of p-grid, so no N_max is given
%         detailed_data and model are saved
% step 51: plot of p-grid of previous step
%           ???other plots useful????
%
% step 6: Basisgeneration using an adaptive p-partition given an N_max of
%         basis vectors per partition and a fixed-greedy algorithm using
%         the same M-Train grid as in the previous two steps.
%         detailed_data and the model are saved
%
% step 61: plot of the p-grid... other plots???
%
%
%
% Bernard Haasdonk 4.2.2010
% Markus Dihlmann 15.02.2010

if nargin < 1
  step = 1; 
end;


% Options:
% ---------------------------

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';      

% Set Parameter-space dimensions:
mu_dimensions ='2D-case';
%mu_dimensions ='3D-case';
params.mu_dimensions=mu_dimensions;
random_seed = 1234;



% Set RB-construction method:
%RB_method = 'POD_of_single_trajectory';
%RB_method = 'POD_of_several_trajectories';
%RB_method = 'greedy_random_fixed';
%RB_method = 'greedy_uniform_fixed';
%RB_method = 'greedy_refined_uniform';
%RB_method = 'greedy_refined_adaptive'; 
%RB_method = 'p_part_model_uniform';
RB_method = 'p_part_model_adaptive';

% Set, if plot should be colored or not:
options=[];
  options.simple_grid=1; %if =1 plot is not colored
  


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);
  switch mu_dimensions
      case '2D-case'
          mu=[0.5,1];
      case '3D-case'
          mu = [0.5,1,1];
  end

  % 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(rbmatlabtemp,[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(rbmatlabtemp,[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(rbmatlabtemp,rb_fname),'RB');
  
  ds_model.RB_basis_filename = fullfile(rbmatlabtemp,rb_fname);
  ds_model.RB_generation_mode = 'file_load'; 
    
  detailed_data = gen_detailed_data(ds_model,ds_model_data);

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

  ds_model = model;
  
  %perform single reduced simulation for given parameter
  online_data = gen_online_data(ds_model,detailed_data);
  ds_model.N = online_data.N;
  tic,
%  ds_model.nt = ds_model.nt/4;
%  disp('set nt to quarter!!');
  rb_sim_data = rb_simulation(ds_model,online_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;
  model.get_rb_from_detailed_data = @(detailed_data) detailed_data.V; %only improvised help-method for method rb_test_indicator.m
  
  switch mu_dimensions
      case '2D-case'
          mu = [0.5,1];
          model.RB_numintervals = [1,1];
      case '3D-case'
          mu = [0.5,1,1];
          model.RB_numintervals = [1,1,1];
  end

  % set mu in model and base_model!!!
  model = model.set_mu(model,mu);
  model.RB_stop_Nmax = 150;

  
  switch RB_method
      
      case 'POD_of_single_trajectory'
        model.RB_generation_mode = 'PCA_trajectory'; 
        
      case 'POD_of_several_trajectories'
          %%% basis generation. POD of several trajectories
          model.RB_generation_mode = 'PCA_trajectories'; 
          switch mu_dimensions
              case '2D-case'
                  model.RB_mu_list = {[0,0],[1,1]};
              case '3D-case'
                  model.RB_mu_list = {[0,0,0],[1,1,1]};
          end
          
      case 'greedy_random_fixed'
          %%% 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;
          
      case 'greedy_uniform_fixed'
          %%% basis generation: greedy_uniform_fixed
          model.RB_generation_mode = 'greedy_uniform_fixed'; 
          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;
      
      case 'greedy_refined_uniform'
          %%% basis generation: greedy_refined with uniform refinement
          model.RB_generation_mode = 'greedy_refined'; 
          model.RB_stop_max_val_train_ratio = 0.95;
          model.RB_refinement_mode = 'uniform'; 
          model.RB_val_rand_seed = 543;
          model.RB_M_val_size = 10;
          model.RB_max_refinement_level = 5;
          model.RB_error_indicator = 'estimator';
          model.RB_extension_algorithm = @RB_extension_PCA_fixspace;
          model.RB_stop_epsilon = 1e-8;
          model.RB_stop_timeout = 60*5; % 5 minutes
          model.RB_stop_Nmax = 4;
          %model.mu_ranges = {[0,1.0],[0,1.0],[0,1.0]};%,[0.4,0.6]};

      case 'greedy_refined_adaptive'
          % 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 = 0.7;
          model.RB_refinement_theta = 0.2;
          model.RB_element_indicator_mode='nodes_cogs_skippedrefs';
          model.RB_element_indicator_s_max = 5;
          model.RB_val_rand_seed = 543;
          model.RB_M_val_size = 10;
          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 = 45*60; % ... minutes
          model.RB_stop_Nmax = 40;
  
      case 'p_part_model_uniform'
          % convert model into a p_part model
          %!!usual model attributes must be fixed 
          model.RB_generation_mode = 'greedy_refined'; 
          model.RB_stop_max_val_train_ratio = 0.95;
          model.RB_refinement_mode = 'uniform'; 
          model.RB_val_rand_seed = 543;
          model.RB_M_val_size = 10;
          model.RB_max_refinement_level = 5;
          model.RB_error_indicator = 'estimator';
          model.RB_extension_algorithm = @RB_extension_PCA_fixspace;
          model.RB_stop_epsilon = 1e-8;
          model.RB_stop_timeout = 60*5; % 5 minutes
          model.RB_stop_Nmax = 4;
          params.p_part_generation_mode = 'uniform';
          switch mu_dimensions
              case '2D-case'
                  params.p_part_numintervals = [2,2];
                  model = p_part_model(model,params);
              case '3D-case'
                  params.p_part_numintervals = [2,2,2];
                  model = p_part_model(model,params);
          end
          
      case 'p_part_model_adaptive'
          model.verbose=8;
          model.RB_generation_mode = 'greedy_refined'; 
          model.RB_refinement_mode = 'adaptive'; 
          model.RB_stop_max_val_train_ratio = 2;
          model.RB_refinement_theta = 0.2;
          model.RB_element_indicator_mode='nodes_cogs_skippedrefs';
          model.RB_element_indicator_s_max = 5;
          model.RB_val_rand_seed = 543;
          model.RB_M_val_size = 10;
          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 = 15*60*60; % .
          model.RB_stop_Nmax = 2;       %maximal number of RB-vectors per p_partition
          model.p_part_max_refinement_level=1;  %%Maximale refinement level for p_partitions
          % convert model into a p_part model
          params.p_part_generation_mode = 'adaptive';
          switch mu_dimensions
              case '2D-case'
                  params.p_part_numintervals =[1,1];
              case '3D-case'
                  params.p_part_numintervals = [1,1,1];
          end
          model = p_part_model(model,params);
  
          
          
  end
  

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

  %model.verbose=8;
  model_data = gen_model_data(model);
  detailed_data = gen_detailed_data(model,model_data);
  
  save(fullfile(rbmatlabtemp,'p_part_basisgen_data'),'model','detailed_data');
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  
  
  
  
case 3 %load detailed data generated in previous step and perform animation and reduced simulations
  
  load(fullfile(rbmatlabtemp,'p_part_basisgen_data'));      
        
  % plot an animation of mesh refinement during basis generation
   
  if ~(strcmp(RB_method,'POD_of_single_trajectory')||strcmp(RB_method,'POD_of_several_trajectories')||strcmp(RB_method,'p_part_model_uniform')||strcmp(RB_method,'p_part_model_adaptive'))
    animate_basisgen(detailed_data,model,options);
  end
  
  %plot partitions of parameter space
  if (strcmp(RB_method,'p_part_model_uniform')||strcmp(RB_method,'p_part_model_adaptive'))
    plot_options=[];
    %nleaf_elements=get(detailed_data.pgrid,'nleafelements');
    %RB_N_str=[];
    %for i=1:nleaf_elements
    %   gid = lid2gid(detailed_data.pgrid,i);
    %   RB_N_str{i}=num2str(size(detailed_data.parts_detailed_data{gid}.RB_info.mu_ind_seq,2));
    %end
    
    figure(1)
    % plot_options.text=RB_N_str;
    plot_grid(detailed_data.pgrid,plot_options);
    
    %plot grid from greedy adaptive into the existant p-part grid
    if isfield(model,'base_model')
        if (isfield(model.base_model,'RB_refinement_mode')&&(strcmp(model.base_model.RB_refinement_mode,'adaptive')))
            options.figureNr=1;
            options.color='r';
            options.text=plot_options.text;
            animate_basisgen(detailed_data,model,options);
        end
    end
    
  end
  
  save(fullfile(rbmatlabtemp,detailed_data_fname),...
       'model','detailed_data')
 
   %model = model;
   
  % compare detailed and reduced simulation
  load(fullfile(rbmatlabtemp,detailed_data_fname),...
       'model','detailed_data')
  load(fullfile(rbmatlabtemp,[trajectory_fnbase,'1']),...
       'ds_sim_data');
  %model = model;
  
  %perform single reduced simulation for given parameter
  online_data = gen_online_data(model,detailed_data);
 if ~(strcmp(RB_method,'p_part_model_uniform')||strcmp(RB_method,'p_part_model_adaptive'))
  model.N = online_data.N;
 end
  tic,
%  model.nt = model.nt/4;
%  disp('set nt to quarter!!');
  rb_sim_data = rb_simulation(model,online_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
  if ~(strcmp(RB_method,'p_part_model_uniform')||strcmp(RB_method,'p_part_model_adaptive'))
    lin_evol_model = model.base_model;
  else
    lin_evol_model = model.base_model.base_model;
  end
  
  
  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 31 %load detailed data generated in previous step and perform animation and reduced simulations

  detailed_data_fname = 'save_case_4'    
      
  load(fullfile(rbmatlabtemp,detailed_data_fname));      
        
  % plot an animation of mesh refinement during basis generation
   
%  if ~(strcmp(RB_method,'POD_of_single_trajectory')||strcmp(RB_method,'POD_of_several_trajectories')||strcmp(RB_method,'p_part_model_uniform')||strcmp(RB_method,'p_part_model_adaptive'))
%    animate_basisgen(detailed_data,model,options);
% end
  
%  %plot partitions of parameter space
%  if (strcmp(RB_method,'p_part_model_uniform')||strcmp(RB_method,'p_part_model_adaptive'))
%    plot_options=[];
%    %nleaf_elements=get(detailed_data.pgrid,'nleafelements');
%    %RB_N_str=[];
%    %for i=1:nleaf_elements
%    %   gid = lid2gid(detailed_data.pgrid,i);
%    %   RB_N_str{i}=num2str(size(detailed_data.parts_detailed_data{gid}.RB_info.mu_ind_seq,2));
%    %end
    
%    figure(1)
    % plot_options.text=RB_N_str;
%    plot_grid(detailed_data.pgrid,plot_options);
    
    %plot grid from greedy adaptive into the existant p-part grid
%    if isfield(model,'base_model')
%        if (isfield(model.base_model,'RB_refinement_mode')&&(strcmp(model.base_model.RB_refinement_mode,'adaptive')))
%            options.figureNr=1;
%            options.color='r';
%            options.text=plot_options.text;
%            animate_basisgen(detailed_data,model,options);
%        end
%    end
    
%  end
  
%  save(fullfile(rbmatlabtemp,detailed_data_fname),...
%       'model','detailed_data')
 
   %model = model;
   
  % compare detailed and reduced simulation
  load(fullfile(rbmatlabtemp,detailed_data_fname),...
       'model','detailed_data')
   
  tr_fn =  fullfile(rbmatlabtemp,[trajectory_fnbase,'1']);
  if exist(tr_fn)
    load(tr_fn,...
         'ds_sim_data');
  else
    model_data = gen_model_data(model);
    ds_sim_data = detailed_simulation(model,model_data);
  end;
  %model = model;
  model.N = model.get_rb_size(model,detailed_data);
  
  %perform single reduced simulation for given parameter
  online_data = gen_online_data(model,detailed_data);
  RB_method = model.RB_generation_mode;
  if ~(strcmp(RB_method,'p_part_model_uniform')||strcmp(RB_method,'p_part_model_adaptive'))
    model.N = online_data.N;
  end
  tic,
  %  model.nt = model.nt/4;
  %  disp('set nt to quarter!!');
  rb_sim_data = rb_simulation(model,online_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
  if ~(strcmp(RB_method,'p_part_model_uniform')||strcmp(RB_method,'p_part_model_adaptive'))
    lin_evol_model = model.base_model;
  else
    lin_evol_model = model.base_model.base_model;
  end
  
  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 4
% step 4: Basisgeneration using NO p-partition and a fixed-greedy
%         algorithm.
%         M-Train for greedy can be fixed, but will probably be a 3^p grid
%         detailed_data and model saved in mat file
%

            model = fast_model(params);
            model.save_time_indices = 0:model.nt;
            model.theta = 0;
            model.get_rb_from_detailed_data = @(detailed_data) detailed_data.V; %only improvised help-method for method rb_test_indicator.m

            switch mu_dimensions
              case '2D-case'
                  mu = [0.5,1];
                  model.RB_numintervals = [3,3];
              case '3D-case'
                  mu = [0.5,1,1];
                  model.RB_numintervals = [2,2,2];
            end

            % set mu in model and base_model!!!
            model = model.set_mu(model,mu);

            model.RB_generation_mode = 'greedy_uniform_fixed'; 
            model.RB_error_indicator = 'estimator';
            model.RB_extension_algorithm = @RB_extension_PCA_fixspace;
            model.RB_stop_epsilon = 1e-2;
            model.RB_stop_timeout = 12*60*60; % 
            model.RB_stop_Nmax = 1000;
            model.verbose=8;
            
            
            model.filecache_ignore_fields_in_model=[model.filecache_ignore_fields_in_model,'mu_ranges']; %for making possible to read data from cache in p-partitions
                     
            model_data = gen_model_data(model);
            overall_generation_time = tic;
            detailed_data = gen_detailed_data(model,model_data);
            t_gen = toc(overall_generation_time);    %time necessary for generating the basis
            
            %construct name of file to save:
            %case_epsilon_N_parametergreedybase_coarsefactor
            file_name = ['save_case_4_',num2str(model.RB_stop_epsilon),'_',num2str(params.coarse_factor),'_',num2str(model.RB_numintervals(1)),'.mat'];
            save(fullfile(rbmatlabtemp,file_name),'model','detailed_data','t_gen');
            
            

    case 41
        
% step 41: mat file from previous step is loaded. 
%          plot of grid?
%

            load(fullfile(rbmatlabtemp,loadfilename));
            options.simple_grid=1; %if =1 plot is not colored
            model.random_seed=random_seed;
            
            
            
            online_data = gen_online_data(model,detailed_data);
            model.N = online_data.N; %Unschoen aber nötig...
            
            %reduced_timer = tic;
            %model=set_mu(model,[0.5 0.5]);
            %  model.nt = model.nt/4;
            %  disp('set nt to quarter!!');
            %rb_sim_data = rb_simulation(model,online_data);
            %t_red = toc(reduced_timer);
            
            model.random_seed=1234;%debug ...raus machen!!!
            [test_error,test_error_max,av_rb_sim_time] = calculate_test_error(model, online_data);
            
            %rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
            disp(['time for reduced simulation: ',num2str(av_rb_sim_time)]);
            disp('Results of basis generation without p_partition:');
            disp(['intervalls:',num2str(model.RB_numintervals(1))]);
            disp(['Test Error Average(estimated error): ',num2str(test_error)]);
            disp(['Maximal Test Error (estimated error):',num2str(test_error_max)]);
            disp(['Time for generation of detailed_data: ',num2str(t_gen)]);
            disp(['Epsilon criteria to stop was: ',num2str(model.RB_stop_epsilon)]);
            disp(['Nmax was: ',num2str(model.RB_stop_Nmax)]);
            disp(['Number of basis vectors: ',num2str(length(detailed_data.RB_info.mu_sequence))]);
            if(detailed_data.RB_info.stopped_on_epsilon)
                disp('Stopped basis generation due to reached epsilon criteria');
            end
            if(detailed_data.RB_info.stopped_on_timeout)
                disp('Stopped basis generation due to reached time limit');
            end
            if(detailed_data.RB_info.stopped_on_Nmax)
                disp('Stopped basis generation due to reached maximum number of basis vectors ');
            end
            
            if strcmp(mu_dimensions,'2D-case')
                disp(['vy_weight was: ',num2str(model.base_model.vy_weight)]);
            end
            
            %Write results to a protocol file
            settings=[];
            settings.setting = 'Without p-partitioning';
            settings.file_name = loadfilename;
            settings.Greedy_intervals = model.RB_numintervals(1);
            settings.stop_on_epsilon = model.RB_stop_epsilon;
            settings.stop_on_N_max = model.RB_stop_Nmax;
            if strcmp(mu_dimensions,'2D-case')
            settings.vy_weight = model.base_model.vy_weight;
            end
            results = [];
            results.Number_of_basisvectors_used = length(detailed_data.RB_info.mu_sequence);
            results.average_test_error = test_error;
            results.maximal_test_error = test_error_max;
            results.time_for_reduced_simulation = av_rb_sim_time;
            results.time_for_generation = t_gen;
            
              
            % plot basisgeneration plot
            animate_basisgen(detailed_data,model,options);
            
            save_result_protocol('p_part_test protokoll',settings, results);
            
    case 5
        
% step 5: Basisgeneration using a fixed p-partition and a fixed-greedy
%         algorithm using the same M-Train grid as the previous step
%         No adaption of p-grid, so no N_max for refinement is given
%         detailed_data and model are saved

             model = fast_model(params);
            model.save_time_indices = 0:model.nt;
            model.theta = 0;
            model.get_rb_from_detailed_data = @(detailed_data) detailed_data.V; %only improvised help-method for method rb_test_indicator.m

            switch mu_dimensions
              case '2D-case'
                  mu = [0.5,1];
                  model.RB_numintervals = [2,2];
              case '3D-case'
                  mu = [0.5,1,1];
                  model.RB_numintervals = [2,2,2];
            end

            % set mu in model and base_model!!!
            model = model.set_mu(model,mu);
            % convert model into a p_part model
            %!!usual model attributes must be fixed 
            model.verbose=8;
            model.RB_generation_mode = 'greedy_uniform_fixed'; 
            model.RB_error_indicator = 'estimator';
            model.RB_extension_algorithm = @RB_extension_PCA_fixspace;
            model.RB_stop_epsilon = 1e-4;
            model.RB_stop_timeout = 60*60*12; 
            model.RB_stop_Nmax = 1000;
            model.filecache_ignore_fields_in_model=[model.filecache_ignore_fields_in_model,'mu_ranges']; %for making possible to read data from cache in p-partitions
            params.p_part_generation_mode = 'uniform';
            switch mu_dimensions
              case '2D-case'
                  params.p_part_numintervals = [2,2];
                  model = p_part_model(model,params);
              case '3D-case'
                  params.p_part_numintervals = [2,2,2];
                  model = p_part_model(model,params);
            end
            
            model_data = gen_model_data(model);
            overall_generation_time = tic;
            detailed_data = gen_detailed_data(model,model_data);
            t_gen = toc(overall_generation_time);    %time necessary for generating the basis
            
            file_name = ['save_case_5_',num2str(model.base_model.RB_stop_epsilon),'_',num2str(params.coarse_factor),'_',num2str(model.base_model.RB_numintervals(1)),'.mat'];
            save(fullfile(rbmatlabtemp,file_name),'model','detailed_data','t_gen');
            
            

    case 51
% step 51: plot of p-grid of previous step
%           ???other plots useful????
%
            load(fullfile(rbmatlabtemp,loadfilename));
            options.simple_grid=1; %if =1 plot is not colored
            model.random_seed=random_seed;
            
            
            
            online_data = gen_online_data(model,detailed_data);
            %model.N = online_data.N; %Unschoen aber nötig...
            
            nleafs = get(detailed_data.pgrid, 'nleafelements');
            leafcogs = get(detailed_data.pgrid, 'leafcogs');
            %reduced_timer = tic;
            %for nl=1:nleafs
            %    model=set_mu(model,leafcogs(nl,:));
            %    rb_sim_data = rb_simulation(model,online_data);
                %rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
            %end
            %t_red=toc(reduced_timer)/nleafs;
            
            %debug
            t_gen =0;
            %end debug ...entfernen sobald aktuelle simulationen vorhanden 
            
            [test_error,test_error_max,av_rb_sim_time] = calculate_test_error(model, online_data);
            
            disp(['average time for reduced simulation: ',num2str(av_rb_sim_time)]);
            
            disp('Results of basis generation with fixed p_partition:');
            disp(['Test Error (estimated error): ',num2str(test_error)]);
            disp(['Maximal Test Error (estimated error):',num2str(test_error_max)]);
            disp(['Time for generation of detailed_data: ',num2str(t_gen)]);
            disp(['Epsilon criteria to stop was: ',num2str(model.base_model.RB_stop_epsilon)]);
            disp(['Nmax was: ',num2str(model.base_model.RB_stop_Nmax)]);
            disp(['Number of partitions: ',num2str(nleafs)]);
            for p=1:length(detailed_data.parts_detailed_data)
                disp(['Partition ',num2str(p),' :']);
                disp(['Number of basis vectors in part ',num2str(p),' : ',num2str(length(detailed_data.parts_detailed_data{p}.RB_info.mu_sequence))]);
                if(detailed_data.parts_detailed_data{p}.RB_info.stopped_on_epsilon)
                    disp('Stopped basis generation due to reached epsilon criteria');
                end
                if(detailed_data.parts_detailed_data{p}.RB_info.stopped_on_timeout)
                    disp('Stopped basis generation due to reached time limit');
                end
                if(detailed_data.parts_detailed_data{p}.RB_info.stopped_on_Nmax)
                    disp('Stopped basis generation due to reached maximum number of basis vectors ');
                end
            end
            
            if strcmp(mu_dimensions,'2D-case')
                disp(['vy_weight was: ',num2str(model.base_model.base_model.vy_weight)]);
            end
            
            %Write results to a protocol file
            settings=[];
            settings.setting = 'With 2X2 fixed p-partitioning';
            settings.file_name = loadfilename;
            settings.Greedy_intervals = model.base_model.RB_numintervals(1);
            settings.stop_on_epsilon = model.base_model.RB_stop_epsilon;
            settings.stop_on_N_max = model.base_model.RB_stop_Nmax;
            if strcmp(mu_dimensions,'2D-case')
            settings.vy_weight = model.base_model.base_model.vy_weight;
            end
            results = [];
            results.average_test_error = test_error;
            results.maximal_test_error = test_error_max;
            results.time_for_reduced_simulation = av_rb_sim_time;
            results.time_for_generation = t_gen;
            for p=1:length(detailed_data.parts_detailed_data)
                results.(['Number_of_basisvectors_used_in',num2str(p)]) = length(detailed_data.parts_detailed_data{p}.RB_info.mu_sequence);
            end
            
            animate_basisgen(detailed_data,model,options);
            
            save_result_protocol('p_part_test protokoll',settings, results);
            
    case 6
% step 6: Basisgeneration using an adaptive p-partition given an N_max of
%         basis vectors per partition and a fixed-greedy algorithm using
%         the same M-Train grid as in the previous two steps.
%         detailed_data and the model are saved
%

        model = fast_model(params);
        model.save_time_indices = 0:model.nt;
        model.theta = 0;
        model.get_rb_from_detailed_data = @(detailed_data) detailed_data.V; %only improvised help-method for method rb_test_indicator.m

        switch mu_dimensions
          case '2D-case'
              mu = [0.5,1];
              model.RB_numintervals = [2,2];
          case '3D-case'
              mu = [0.5,1,1];
              model.RB_numintervals = [2,2,2];
        end

        % set mu in model and base_model!!!
        model = model.set_mu(model,mu);
        model.verbose=8;
        
        model.RB_generation_mode = 'greedy_uniform_fixed'; 
        model.RB_error_indicator = 'estimator';
        model.RB_extension_algorithm = @RB_extension_PCA_fixspace;
        model.RB_stop_epsilon = 1e-4;
        model.RB_stop_timeout = 24*60*60; % .
        model.RB_stop_Nmax = 80;       %maximal number of RB-vectors per p_partition
        model.p_part_max_refinement_level=3;  %%Maximale refinement level for p_partitions
        model.filecache_ignore_fields_in_model=[model.filecache_ignore_fields_in_model,'mu_ranges']; %for making possible to read data from cache in p-partitions
        % convert model into a p_part model
        params.p_part_generation_mode = 'adaptive';
        switch mu_dimensions
          case '2D-case'
              params.p_part_numintervals =[1,1];
          case '3D-case'
              params.p_part_numintervals = [1,1,1];
        end
        model = p_part_model(model,params);
        
        model_data = gen_model_data(model);
        overall_generation_time = tic;
        detailed_data = gen_detailed_data(model,model_data);
        detailes_data = clean_up_part_detailed_data(detailed_data);
        t_gen = toc(overall_generation_time);    %time necessary for generating the basis
        
        file_name = ['save_case_6_',num2str(model.base_model.RB_stop_epsilon),'_',num2str(model.base_model.RB_stop_Nmax),'_',num2str(params.coarse_factor),'_',num2str(model.base_model.RB_numintervals(1)),'_',num2str(model.base_model.p_part_max_refinement_level),'.mat'];
        save(fullfile(rbmatlabtemp,file_name),'model','detailed_data','t_gen');

    case 61
% step 61: plot of the p-grid... other plots???
%

        load(fullfile(rbmatlabtemp,loadfilename));
            options.simple_grid=1; %if =1 plot is not colored
            model.random_seed = random_seed;
            
            
            
            
            online_data = gen_online_data(model,detailed_data);
            %model.N = online_data.N; %Unschoen aber nötig...
            
            nleafs = get(detailed_data.pgrid, 'nleafelements');
            leafcogs = get(detailed_data.pgrid, 'leafcogs');
            %reduced_timer=tic;
            %for nl=1:nleafs
            %    model=set_mu(model,leafcogs(nl,:));
            %    rb_sim_data = rb_simulation(model,online_data);
                %rb_sim_data = rb_reconstruction(model,detailed_data,rb_sim_data);
            %end
            %t_red=toc(reduced_timer)/nleafs;
            
            [test_error,test_error_max,av_rb_sim_time] = calculate_test_error(model, online_data);
            
            disp(['average time for reduced simulation: ',num2str(av_rb_sim_time)]);
            
            disp('Results of basis generation with adaptive p_partition:');
            disp(['Test Error (estimated error): ',num2str(test_error)]);
            disp(['Maximal Test Error (estimated error):',num2str(test_error_max)]);
            disp(['Time for generation of detailed_data: ',num2str(t_gen)]);
            disp(['Epsilon criteria to stop was: ',num2str(model.base_model.RB_stop_epsilon)]);
            disp(['Nmax was: ',num2str(model.base_model.RB_stop_Nmax)]);
            disp(['Number of partitions: ',num2str(nleafs)]); 
            for p=1:length(detailed_data.parts_detailed_data)
                disp(['Partition ',num2str(p),' :']);
                disp(['Number of basis vectors in part ',num2str(p),' : ',num2str(length(detailed_data.parts_detailed_data{p}.RB_info.mu_sequence))]);
                if(detailed_data.parts_detailed_data{p}.RB_info.stopped_on_epsilon)
                    disp('Stopped basis generation due to reached epsilon criteria');
                end
                if(detailed_data.parts_detailed_data{p}.RB_info.stopped_on_timeout)
                    disp('Stopped basis generation due to reached time limit');
                end
                if(detailed_data.parts_detailed_data{p}.RB_info.stopped_on_Nmax)
                    disp('Stopped basis generation due to reached maximum number of basis vectors ');
                end
            end

            if strcmp(mu_dimensions,'2D-case')
                disp(['vy_weight was: ',num2str(model.base_model.base_model.vy_weight)]);
            end
            
            %Write results to a protocol file
            settings=[];
            settings.setting = 'With adaptive p-partitioning';
            settings.file_name = loadfilename;
            settings.Greedy_intervals = model.base_model.RB_numintervals(1);
            settings.stop_on_epsilon = model.base_model.RB_stop_epsilon;
            settings.stop_on_N_max = model.base_model.RB_stop_Nmax;
            if strcmp(mu_dimensions,'2D-case')
            settings.vy_weight = model.base_model.base_model.vy_weight;
            end
            results = [];
            results.average_test_error = test_error;
            results.maximal_test_error = test_error_max;
            results.time_for_reduced_simulation = av_rb_sim_time;
            results.time_for_generation = t_gen;
            results.Number_of_partitions = nleafs;
            for p=1:length(detailed_data.parts_detailed_data)
                results.(['Number_of_basisvectors_used_in',num2str(p)]) = length(detailed_data.parts_detailed_data{p}.RB_info.mu_sequence);
            end
            
            animate_basisgen(detailed_data,model,options);
            
            save_result_protocol('p_part_test protokoll',settings, results);
            
            
            
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% step = 7 % time measurements
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 case 7
  
%  clear all;
  nreps = 10;
  t_detailed = zeros(1,nreps);
  
  % initialize detailed model and simulate
  
  load(fullfile(rbmatlabtemp,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(rbmatlabtemp,rb_fname);
  detailed_data = gen_detailed_data(model,model_data);
  online_data = gen_online_data(model,detailed_data);

  Ns = [150:-10:10]; %,30,20,10];
  Ns = Ns(find(online_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,online_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(rbmatlabtemp,'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(rbmatlabtemp,[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(rbmatlabtemp,[detailed_data_fname,'_greedy']);
    load(dfname);
    disp('loaded basis from step 9')
  else
    load(fullfile(rbmatlabtemp,detailed_data_fname),...
         'model','detailed_data') 
    
%    model = fast_model(params);
    model_data = gen_model_data(model);
    %  dfname = fullfile(rbmatlabtemp,'advection_fv_output_detailed_data_146');
%    model.RB_basis_filename = ...
%       fullfile(rbmatlabtemp,'advection_fv_output_basis');
%    model.RB_generation_mode = 'file_load';
%    detailed_data = gen_detailed_data(model,model_data);  
    %  dfname = fullfile(rbmatlabtemp,'advection_fv_output_detailed');
    % load(dfname);
    tmp = load(fullfile(rbmatlabtemp,[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;
  online_data = gen_online_data(model,detailed_data);
  rb_sim_data = rb_simulation(model,online_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

  online_data = gen_online_data(model,detailed_data);
  
  Ns = 150:-10:10;
  Ns = Ns(find(online_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, online_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(rbmatlabtemp,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(rbmatlabtemp,'advection_fv_output_detailed_data_146');
%  model.RB_basis_filename = ...
%      fullfile(rbmatlabtemp,'advection_fv_output_basis');
%  model.RB_generation_mode = 'file_load';
%  detailed_data = gen_detailed_data(model,model_data);  
%  save(fullfile(rbmatlabtemp,...
%               'advection_fv_output_detailed_2trajectories.mat'),...
%       'model','detailed_data')
%  keyboard;
  %  dfname = fullfile(rbmatlabtemp,'advection_fv_output_detailed');
  % load(dfname);
  tmp = load(fullfile(rbmatlabtemp,[trajectory_fnbase,'2']));
  ds_sim_data = tmp.ds_sim_data;
  disp('loaded basis from step 4')
  
  %  factors = 0:0.05:1;
  online_data = gen_online_data(model,detailed_data);
  Ns = [1,8,16,32,64];
  Ns = Ns(find(online_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, online_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(rbmatlabtemp,'parameter_sweep.eps'),'epsc');
  keyboard;
 otherwise
  disp('step number unknown');
end;

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;
switch params.mu_dimensions
    case '2D-case'
        lin_evol_model.mu_names = {'cone_amplitude','vx_weight'};
        lin_evol_model.mu_ranges = {[0,1],[0,1]};
    case '3D-case'
        lin_evol_model.mu_names = {'cone_amplitude','vx_weight','vy_weight'};
        lin_evol_model.mu_ranges = {[0,1],[0,1],[0,1]};
end
%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;


end


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

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

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;
end

function Bcomp = B_components(model,model_data) 
Bcomp = model_data.B_components;
end

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);
end

function [epsilon_test,epsilon_max,rb_sim_av_time] = calculate_test_error(model, online_data)
    display('Calculating estimated test error');
    N_test=500;%amount of Tests
    rand('state',model.random_seed);
    
    random_mu = rand_uniform(N_test, model.mu_ranges);
    
    error_max=0;
    rb_sim_time=0;
    est_error=zeros(1,N_test);
    for i=1:N_test
       %Choose a mu from range by random choice
       mu = random_mu(:,i);
       model = set_mu(model,mu);
       %perform a reduced simulation and take the estiamated error of last
       %time step
       rb_sim_timer = tic;
       rb_sim_data = rb_simulation(model, online_data);
       rb_sim_time = rb_sim_time + toc(rb_sim_timer);
       %add error to total error
       est_error(i) = rb_sim_data.DeltaX(end);
       
       if~(isfield(model,'verbose'))
           model.verbose = 20;
       end
       
       
       if(model.verbose>5)
           if (mod(i,10)==0)
            disp('...')
           end
       end
    end
   
    
    %calculate average estimated error
   epsilon_test = sum(est_error)/N_test;
   epsilon_max=max(est_error);
   rb_sim_av_time = rb_sim_time/N_test;
   
   %Visualisation of error map on parameter domain
   figure;
   scatter(random_mu(1,:),random_mu(2,:),40,est_error,'filled');
   colorbar;
end

function save_result_protocol(filename, settings, results)
    
    fid = fopen(fullfile(rbmatlabhome,filename),'a');
    fprintf(fid,'----------------------------------------------------------\n');
    fprintf(fid,'----------------------------------------------------------\n');
    fprintf(fid,'%s\n',settings.setting);
    fprintf(fid,'%s\n',settings.file_name);
    fprintf(fid,'-----------------------------------------------------------\n');
    
    settingnames=fieldnames(settings);
    resultnames=fieldnames(results);
    
    fprintf(fid,'Settings:\n');
    for i=1:length(settingnames)
        if isscalar(getfield(settings,settingnames{i}))
            text=num2str(getfield(settings,settingnames{i}));
        else
            text = num2str(getfield(settings,settingnames{i}));
        end
        fprintf(fid,'%s:\t %s\n',settingnames{i},text);
    end
    fprintf(fid,'\nResults:\n');
    for i=1:length(resultnames)
        if isscalar(getfield(results,resultnames{i}))
            text=num2str(getfield(results,resultnames{i}));
        else
            text = num2str(getfield(results,resultnames{i}));
        end
        fprintf(fid,'%s:\t %s\n',resultnames{i},text);
    end
    
    %save figures
    saveas(figure(1),fullfile(rbmatlabhome,'protocol_figures',[settings.file_name,'_error_map.fig']),'fig');
    saveas(figure(1),fullfile(rbmatlabhome,'protocol_figures',[settings.file_name,'_error_map.jpg']),'jpg');
    saveas(figure(1),fullfile(rbmatlabhome,'protocol_figures',[settings.file_name,'_error_map.eps']),'eps');
    clf;
    close(figure(1));
    saveas(figure(2),fullfile(rbmatlabhome,'protocol_figures',[settings.file_name,'.fig']),'fig');
    saveas(figure(2),fullfile(rbmatlabhome,'protocol_figures',[settings.file_name,'.jpg']),'jpg');
    saveas(figure(2),fullfile(rbmatlabhome,'protocol_figures',[settings.file_name,'.eps']),'eps');
    clf;
    close(figure(2));
    fclose(fid);
end