RB_greedy_extension.m

Go to the documentation of this file.

function [detailed_data, reduced_data] = RB_greedy_extension(model,detailed_data)
%function [detailed_data,reduced_data] = ...
%             RB_greedy_extension(model,detailed_data);
%  function performing a greedy search loop for reduced basis
%  generation.
%
%  In each iteration the parameter vector with maximum posterior estimator or
%  true error is selected for "refinement". An orthonormal basis is generated.
%
%  The detailed_data is expected to be reasonably filled, that rb_simulations
%  are possible with the given arguments, i.e.  gen_reduced_data(model,
%  detailed_data). That means, that it also can contain an existing basis,
%  which is then extended in this loop.  the offline_data for after the last
%  extension is optionally returned.
%
%  The stopping of the extention is performed, if
%
%    - 'model.RB_stop_epsilon' is reached: basis vectors are added until
%             the maximum posterior estimator is smaller than this
%             value. M_val can be ommited in this case.
%    - 'model.RB_stop_max_val_train_ratio': basis vectors are added until the
%             ratio r of the max estimator on a validation set M_val with the
%             estimator on the training set M is larger than r_max. M_val is
%             required in this case
%    - 'model.RB_stop_timeout': basis vectors are added until toc gives a
%             time exceeding a maximum level given
%    - 'model.RB_stop_Nmax': If given number of basis vectors in full
%             basis is obtained.
%    - 'model.p_part_early_refinement' is set to one and the estimated
%             number of basis vectors necessary to reach etol is higher
%             than 'RB_stop_Nmax'.
%
% The following fields are created/extended in 'detailed_data.RB_info':
%    - 'M_first_errs' : a vector with the initial error-estimator for each
%               M_train. Is determined to check, whether any
%               extension has to be performed or not.
%    - 'stopped_on_epsilon' : flag indicating whether epsilon-stop-crit is reached
%    - 'stopped_on_timeout': flag indicating whether timeout is reached
%    - 'stopped_on_Nmax': flag indicating whether Nmax is reached
%    - 'stopped_on_max_val_train_ratio': flag indicating whether
%               epsilon-stop-crit is reached
%    - 'stopped_on_empty_extension' : flag indicating, whether the last
%               extension failed due to accuracy problems, i.e. new
%               basis-vector already in span of existing.
%    - 'stopped_on_Nlimit_estimation': flag indicating ifthe extension
%               process was stopped, because it was estimated that more
%               than Nmax basis vectoras are needed to reach etol.
%               Partition of the parameter domain can be initiated.
%    - 'max_err_sequence' : the decreasing maximum posterior estimators or true
%               errors in each iteration, which are present AFTER the
%               respective basis extension. So in case of 'stopped_on_epsilon',
%               the last entry of this field should be smaller than the desired
%               epsilon. The first entry is the maximal initial error/estimator.
%    - 'mu_sequence' : sequence of mu_values chosen for
%               basis extension i.e. list of worst-estimator mu_values as the
%               columns of the matrix mu_values. The last entry corresponds
%               to the errors/estimator values AFTER the last basis
%               extension.
%    - 'N_sequence' : sequence of basis sizes AFTER each extension. First
%               entry is basis size before first extension.
%    - 'r_value_sequence' : ratios of val/train error-estimator in case of
%               validation set providing, AFTER the respective basis-extension.
%    - 'toc_value_sequence' : toc_values after each basis extension
%    - 'M_last_errs' : the posterior estimators for all mu in M_train at final
%               time AFTER the last respective basis extension
%
%  Required fields in 'detailed_data.RB_info':
%     - 'M_train' : matrix with parameter-vectors as columns, over which greedy
%                   search is to be performed.
%     - 'M_val'   : matrix with parameter-vectors as columns, which are used
%                   for validation. Only required, if the field
%                   'model.stop_on_max_val_train_ratio' is set)
%
%  The following fields of detailed_data are created/extended:
%     - RB      : is extended by the number of basis-vectors, overall being the
%                 orthonormal basis
%     - RB_info : the fields of this info-structure are extended, as described
%                 above
%
% required fields of model:
%    RB_extension_algorithm: function pointer to either
%                         RB_extension_max_error_snapshot() or
%                         RB_extension_PCA_fixspace()
%    RB_error_indicator: either 'estimator' or 'error' indicating,
%                  whether an aposteriori indicator or true error
%                  is to be used
%    RB_detailed_val_savepath : path to the detailed data of the
%                  validation set in validation mode. Is generated
%                  if not available.
%    RB_detailed_train_savepath : path to the detailed data of the
%                  training set. Is generated if not available. Is
%                  only required in case of 'error' mode instead of
%                  'estimator'
%
% additional fields as required by the selected
% extension algorithm and the reduced-basis-simulation must be set in model.

% Bernard Haasdonk 26.5.2007

if ~isfield(model,'RB_stop_timeout')
    model.RB_stop_timeout             = inf;
end;

if ~isfield(model,'RB_stop_Nmax')
    model.RB_stop_Nmax                = inf;
end;

if ~isfield(model,'RB_stop_epsilon')
    model.RB_stop_epsilon             = 0;
end;

if ~isfield(model,'RB_stop_max_val_train_ratio')
    model.RB_stop_max_val_train_ratio = inf;
end;

if ~isfield(model,'RB_save_temp_detailed_data')
    model.RB_save_temp_detailed_data = 0;
end;

% check, whether val-error is wanted
check_val_error = 0;
if model.RB_stop_max_val_train_ratio < inf
    check_val_error                   = 1;
end;

if isfield(model,'p_part_early_refinement')
    enable_p_partition                = 1;
else
    enable_p_partition                = 0;
end

%hp-greedy:
if isfield(model,'h_part')
    enable_h_partition                = 1;
else
    enable_h_partition                = 0;
end



% make sure, that detailed-simulations for validation set are available
% in case of active validation criterion and wanted true error
if check_val_error
    if isempty(detailed_data.RB_info.M_val)
        error ('M_val must be provided in case of max_val_train_ratio mode!');
    end;
    if isequal(model.RB_error_indicator,'error')
        detailed_val_savepath = model.RB_detailed_val_savepath;
        save_detailed_simulations_ptr = model.save_detailed_simulations;
        save_detailed_simulations_ptr(model,detailed_data, ...
            detailed_data.RB_info.M_val, ...
            detailed_val_savepath);
    else
        detailed_val_savepath = [];
    end;
end;

% make sure, that detailed-simulations for the training set are available in
% case of wanted true error
if isequal(model.RB_error_indicator,'error') | isequal(model.RB_error_indicator,'projection-error')
    save_detailed_simulations_ptr = model.save_detailed_simulations;
    save_detailed_simulations_ptr(model,detailed_data,...
        detailed_data.RB_info.M_train, ...
        model.RB_detailed_train_savepath);
    detailed_train_savepath = model.RB_detailed_train_savepath;
else
    detailed_train_savepath = 'dummy-string';
end;

% initialize detailed_data.RB_info
if ~isfield(detailed_data,'RB_info')
    detailed_data.RB_info                    = [];
end;

N = model.get_rb_size(model,detailed_data);

if ~isfield(detailed_data.RB_info,'max_err_sequence')
    detailed_data.RB_info.max_err_sequence   = nan(1,N);
end
if ~isfield(detailed_data.RB_info,'mu_sequence')
    detailed_data.RB_info.mu_sequence        = ...
        nan(length(model.mu_names), N);
end
if ~isfield(detailed_data.RB_info,'mu_ind_seq')
    detailed_data.RB_info.mu_ind_seq         = ...
        nan(1, N);
end

if check_val_error && ~isfield(detailed_data.RB_info,'r_value_sequence')
    detailed_data.RB_info.r_value_sequence   = nan(1,N);
end

if ~isfield(detailed_data.RB_info,'toc_value_sequence')
    detailed_data.RB_info.toc_value_sequence = nan(1,N);
end

if ~isfield(detailed_data.RB_info, 'N_sequence')
    detailed_data.RB_info.N_sequence = [];
end

% (partially) retain info of previous runs
RB_info = detailed_data.RB_info;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initialize loop-break conditions %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
RB_info.stopped_on_epsilon             = 0;
RB_info.stopped_on_timeout             = 0;
RB_info.stopped_on_Nmax                = 0;
RB_info.stopped_on_max_val_train_ratio = 0;
RB_info.stopped_on_empty_extension     = 0;
RB_info.stopped_on_Nlimit_estimation   = 0;
N_limit_est_av=[]; % average over all estimations (pessimistic)
continue_loop = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initialize output quantities %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
max_pe     = [];
mu_values  = [];
mu_ind_seq = [];
toc_values = [];
if check_val_error
    r_values = [];
end;
N_limit_est = [];% estimated amount of N needed
post_errs_sequence = [];
post_est_sequence = [];

% perform extension loop:
while continue_loop
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % determine mu-worst error-indicator for current training set %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    reduced_data = gen_reduced_data(model,detailed_data);
    
    model.N = reduced_data.N;
    post_errs = rb_test_indicator(model, ...
        detailed_data, reduced_data,...
        detailed_data.RB_info.M_train,...
        detailed_train_savepath);
    % store real error as well, if estimator and error are computed for
    % 'ei_estimator_test'
    if isequal(model.RB_error_indicator,'ei_estimator_test')
        post_errs_sequence = [post_errs_sequence, post_errs.error];
        post_est_sequence  = [post_est_sequence, post_errs.estimator];
        post_errs          = post_errs.estimator;
    end
    
    [max_errs, mu_inds] = max(post_errs);
    max_err             = max_errs(1);      % maximum error (estimation)
    mu_ind              = mu_inds(1);       % parameter index in RB_info.M_train
    
    if isnan(max_err)
        keyboard
    end;
    
    % hp-greedy: first chosen mu is the first training parameter
    if enable_h_partition && model.N == 0
        mu_ind = 1;
    end;
    
    mu = detailed_data.RB_info.M_train(:,mu_ind);
    
    max_pe     = [max_pe, max_err];
    mu_ind_seq = [mu_ind_seq, mu_ind];
    mu_values  = [mu_values, mu];
    
    % update mu-vector is stored in variable mu
    disp(['Detected maximum error prediction ',num2str(max_err),...
        ' for mu=[',num2str(mu'),']']);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % check loop-break conditions %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if max_err < model.RB_stop_epsilon
        continue_loop      = 0;
        RB_info.stopped_on_epsilon = 1;
    end;
    
    if model.get_rb_size(model,detailed_data) >= model.RB_stop_Nmax
        continue_loop   = 0;
        RB_info.stopped_on_Nmax = 1;
    end;
    
    if toc > model.RB_stop_timeout
        RB_info.stopped_on_timeout = 1;
        continue_loop      = 0;
    end;
    
    if check_val_error
        val_est = rb_test_indicator(model,...
            detailed_data,reduced_data,...
            detailed_data.RB_info.M_val, ...
            detailed_val_savepath);
        
        r_values = [r_values, max(val_est) / max_err];
        %debug
        disp(['val_train ratio:  ',num2str(r_values(end))]);
        %end debug
        if (r_values(end) > model.RB_stop_max_val_train_ratio)
            RB_info.stopped_on_max_val_train_ratio = 1;
            continue_loop = 0;
        end;
    end;
    
    %in models using adaptive p-partition early stopping:
    %an estimater extrapolates the curve N->epsilon and estimates if the
    %epsilon_tol can be reached with a basis size smaller than Nmax. If no,
    %the field RB_stopped_on_Nlimit_estimation is set to 1
    
    
    if enable_p_partition && model.p_part_early_refinement
        %estimate N total needed to reach etol:
        err_seq=max_pe;
        etol = model.RB_stop_epsilon;
        Nmax = model.RB_stop_Nmax;
        
        s=20; %stepsize backwards
        
        
        if(length(err_seq)>(s))
            
            %for i=(s+1):length(err_seq);
            i=length(err_seq);
            a=(log(err_seq(i-s)/err_seq(i)))/(s);
            N_limit_est=[N_limit_est,i-1/a*log(etol/err_seq(i))];
            %end;
            
            
            
            %for i=1:length(N_limit_est);
            %summe = summe+N_limit_est(i);
            if(length(N_limit_est_av>0))
                N_1=length(N_limit_est_av);
                newAv=N_1/(N_1+1)*N_limit_est_av(end)+(1/(N_1+1))*N_limit_est(end);
                N_limit_est_av=[N_limit_est_av,newAv];
                
            else
                N_limit_est_av=N_limit_est(end);
            end
            %end;
            disp(['actual estimation of N needed to reach etol: ',num2str(N_limit_est_av(end))]);
            if(N_limit_est_av(end)>Nmax)
                RB_info.stopped_on_Nlimit_estimation = 1;
                continue_loop                = 0;
            end
        end
        
        
    end %of P-Part_early_refinement
    
    % always store computational time
    toc_values = [toc_values, toc];
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % update information structure and possibly save intermediate state to file %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    RB_info.max_err_sequence = [RB_info.max_err_sequence, max_pe(end)];
    RB_info.mu_ind_seq       = [RB_info.mu_ind_seq, mu_ind_seq(end)];
    RB_info.mu_sequence      = [RB_info.mu_sequence, mu_values(:, end)];
    RB_info.N_sequence       = [RB_info.N_sequence, model.N];
    if isequal(model.RB_error_indicator,'ei_estimator_test')
        RB_info.post_errs_sequence = post_errs_sequence;
        RB_info.post_est_sequence  = post_est_sequence;
    end
    if check_val_error
        RB_info.r_value_sequence = [RB_info.r_value_sequence, r_values(end)];
    end;
    RB_info.toc_value_sequence = [RB_info.toc_value_sequence, toc_values(end)];
    
    if enable_p_partition && model.p_part_early_refinement
        RB_info.N_limit_est    = N_limit_est;
        RB_info.N_limit_est_av = N_limit_est_av;
    end
    
    if numel(max_pe) == 1
        RB_info.M_first_errs = post_errs;
    end
    
    detailed_data.RB_info = RB_info;
    
    if model.RB_save_temp_detailed_data
        save(fullfile(rbmatlabtemp,'temp_detailed_data.mat'),'detailed_data', ...
            'model')
    end;
    
    if continue_loop
        
        %%%%%%%%%%%%%%%%
        % extend basis %
        %%%%%%%%%%%%%%%%
        
        % Use original params with new mu such that filecaching works in the
        % extension routine!
        
        %oldmu        = model.get_mu(model);
        %model        = model.set_mu(model, mu);
        %[RBext, par] = model.RB_extension_algorithm(model, ...
        %                                            detailed_data);
        %model        = model.set_mu(model, oldmu);
        oldmu = get_mu(model);
        model = set_mu(model, mu);
        RBext = model.RB_extension_algorithm(model, detailed_data);
        model = set_mu(model, oldmu);
        
        % if empty: no further processing has to be performed
        if ~isempty(RBext)
            
            if model.verbose > 29
                disp(par);
            end;
            RB = [model.get_rb_from_detailed_data(detailed_data), RBext];
            % perform permutation that maintains original order if possible:
            %RBo = model_orthonormalize_qr(model,detailed_data,RB,eps);
            %K = model.inner_product(model,detailed_data,RBo,RB);
            %[m,i] = max(abs(K));
            %RB = RBo(:,i);
            %    if model.debug
            %  K  = model.inner_product(model,detailed_data,RB,RB);
            %  if max(max(abs(K-eye(size(K)))))>1e-10
            %           %error('error in orthonormal basis extension!');
            %  end;
            %   % end;
            detailed_data = model.set_rb_in_detailed_data(...
                detailed_data,...
                RB);
            %detailed_data.RB = [detailed_data.RB, RBext];
            N = model.get_rb_size(model,detailed_data);
            detailed_data.N = N;
            disp(['Extended RB to length ', num2str(N)]);
            
            %orthogonality check
            
            RB = model.get_rb_from_detailed_data(detailed_data);
            W = model.get_inner_product_matrix(detailed_data);
            disp(['checking orthogonality with tolerance ',num2str(1e-12),...
                '...',...
                num2str(max(max(RB'*W*RB-eye(size(RB,2)))))]);
            if(max(max(RB'*W*RB-eye(size(RB,2))))>1e-12)
                disp('orthogonalising again...')
                detailed_data = model.set_rb_in_detailed_data(detailed_data,...
                    model.orthonormalize(RB, W));
                RB = model.get_rb_from_detailed_data(detailed_data);
                disp(['new orthogonality value:  ',num2str(max(max(RB'*W*RB-eye(size(RB,2)))))]);
            end
            
        end % if ~isempty(RB)
        
        % model.N is not set yet: quick-fix of problem with "empty extension" after
        % orthogonalizing
        if model.N == model.get_rb_size(model, detailed_data)
            RB_info.stopped_on_empty_extension = 1;
            continue_loop              = 0;
        end
        
        %  keyboard;
        
        % check on accuracy problem (non-decreasing error)
        %  if length(max_pe)>1 & max_pe(end)>max_pe(end-1)
        %    stopped_on_error_increase = 1;
        %    continue_loop = 0;
        %    disp('stopped on error-increase! Please check!');
        %    keyboard;
        %  else
        %    stopped_on_error_increase = 0;
        %  end;
        
        %    model.Nmax = size(detailed_data.RB,2);
        %    model.N = model.Nmax;
        %    [mu_ind, post_errs] = detect_mu_worst(detailed_data,model);
        %    max_err = max(post_errs);
        
        
    end; % if continue_loop
    
end; % while continue_loop