cubegrid_interp_greedy.m

Go to the documentation of this file.

function [interp_func, varargout] = cubegrid_interp_greedy(func, params)
% function [interp_func, varargout] = cubegrid_interp_greedy(func, params)
% method determining interpolation points for rbf (radial basis function)
% interpolation of a given function via cubegrid refinements in
% greedy-fashion
%
% parameters:
% - func
% - params
%
% required fields of params:
% - ranges

% check for parameters
if ~isfield(params, 'tolerance')
    params.tolerance = 1e-3;
end

if ~isfield(params, 'nmax')
    params.nmax = 4*2^length(params.range);
end

if ~isfield(params, 'domain_corners')
    params.domain_corners = 1;
end

if ~isfield(params, 'midpoint_version')
    params.midpoint_version = 1;
elseif ~params.midpoint_version
    params.domain_corners = 1;
end

if ~isfield(params, 'verbose')
    params.verbose = 0;
end

if ~isfield(params, 'kernel_type')
    params.kernel_type = 'gaussian_loc';
end

if ~isfield(params, 'dilation')
    params.dilation = 0.3;
end

if ~isfield(params, 'penalize_factor')
    params.penalize_factor = 0.5;
end

if params.verbose
    disp('enter cubegrid interpolation greedy');
end

if nargout > 1
    varargout{1} = [];
end

start_algorithm = tic;

% create cubegrid
params.numintervals = ones(length(params.range), 1);
grid = cubegrid(params);

% construct interpolant function
interp_func = RbfInterpolant(zeros(grid.dimension, 0), [], params.kernel_type, params.dilation);

% input for adaptive algorithm
continue_loop = 1;
num_points = 0;
C = [];
C_added = [];
V = [];
penalize_count = 0;
new_ids = 1;
all_points = [];
all_values = [];
max_err_seq = [];
dilation_seq = {};
num_points_seq = [];

if params.verbose
    disp('start greedy algorithm');
end

% start iteration
while continue_loop
    
    % determine new points for interpolant function
    new_vertices = zeros(0, grid.dimension);
    if ~params.midpoint_version || (params.domain_corners && grid.refine_steps == 0)
        
        % determine new vertices
        old_ids = grid.creation_step < grid.creation_step(new_ids(1));
        new_vertex_indices = ...
            setdiff(grid.vertexindex(new_ids, :), grid.vertexindex(old_ids, :));
        
        % detect midpoint vertex
        new_vertices = grid.vertex(new_vertex_indices, :);
        [~, dupl_id, ~] = detect_duplicate_rows([C; new_vertices], 1e-10);
        new_vertices(dupl_id - size(C, 1), :) = [];
    end
    
    % gather midpoints
    C_new = get_midpoints(grid, new_ids);
    C = [C; C_new];
    
    % compute new midpoint data
    tmp_values = cell(length(new_ids), 1);
    for i = 1:length(new_ids)
        tmp_values{i} = func(C_new(i, :)');
    end
    V = [V; cell2mat(tmp_values)];
    
    % compute errors
    errors = min_max_errors(interp_func, C, V, params);
    max_error = max(errors);
    max_err_seq = [max_err_seq, max_error];
    num_points_seq = [num_points_seq, interp_func.num_points];
    if strncmp(interp_func.kernel_type, 'gaussian', 8)
        dilation_seq = [dilation_seq, interp_func.dilation];
    end
    
    % stopping criteria
    continue_loop = max_error >= params.tolerance;
    
    if params.verbose > 1
        disp([num2str(interp_func.num_points), '(', ...
              num2str(num_points), ') points, max error = ', ...
              num2str(max_error)]);
    end
    
%     if num_points > 0 && exist('rbhdd_plot_infsup_results', 'file') == 2
%         infsup_constant = RbfInterpolant(interp_func.data_points, ...
%                                          interp_func.data_values, ...
%                                          interp_func.kernel_type, ...
%                                          interp_func.dilation);
%         rbhdd_plot_infsup_results;
%         if ~isfield(varargout{1}, 'max_test_errs')
%             varargout{1}.max_test_errs = [];
%         end
%         varargout{1}.max_test_errs = [varargout{1}.max_test_errs, act_val];  
%     end
    
    if num_points >= params.nmax
        warning('rbmatlab:cubegrid_interp_greedy', ...
                ['iteration maximum reached, maximal error: ', num2str(max_error)])
        continue_loop = 0;
    end
    
    if continue_loop
        
        % compute new vertex data
        tmp_values = cell(size(new_vertices, 1), 1);
        for i = 1:size(new_vertices, 1)
            tmp_values{i} = func(new_vertices(i, :)');
        end
        
        % add new vertex data
        interp_func.add_data(new_vertices', cell2mat(tmp_values), ...
            (interp_func.dilation(end)+(interp_func.num_points>0)*10*eps) * ...
            ones(size(new_vertices, 1)-(interp_func.num_points==0), 1));
        all_points = [all_points; new_vertices];
        all_values = [all_values; cell2mat(tmp_values)];
        num_points = num_points + size(new_vertices, 1);
        
        % again compute errors
        errors = min_max_errors(interp_func, C, V, params);
        max_error = max(errors);
        max_err_seq = [max_err_seq, max_error];
        num_points_seq = [num_points_seq, interp_func.num_points];
        if strncmp(interp_func.kernel_type, 'gaussian', 8)
            dilation_seq = [dilation_seq, interp_func.dilation];
        end
        
        % stopping criteria
        continue_loop = max_error >= params.tolerance;
        
        if params.verbose > 1
            disp([num2str(interp_func.num_points), '(', ...
                  num2str(num_points), ') points, max error = ', ...
                  num2str(max_error)]);
        end

%         if num_points > 0 && exist('rbhdd_plot_infsup_results', 'file') == 2
%             infsup_constant = RbfInterpolant(interp_func.data_points, ...
%                                              interp_func.data_values, ...
%                                              interp_func.kernel_type, ...
%                                              interp_func.dilation);
%             rbhdd_plot_infsup_results;
%             if ~isfield(varargout{1}, 'max_test_errs')
%                 varargout{1}.max_test_errs = [];
%             end
%             varargout{1}.max_test_errs = [varargout{1}.max_test_errs, act_val];  
%         end
        
        if num_points >= params.nmax
            warning('rbmatlab:cubegrid_interp_greedy', ...
                    ['iteration maximum reached, maximal error: ', num2str(max_error)])
            continue_loop = 0;
        end
    end
    
    if continue_loop
        
        all_points = [all_points; C_new];
        all_values = [all_values; V(end+1-(length(new_ids):-1:1))];
        num_points = num_points + length(new_ids);
        
        % penalize un-refined for long time
        widths = sqrt(get_volume(grid, find(grid.isleaf)));
        inf_ind = errors == inf;
        indicators = params.penalize_factor*max(widths)^2/min(widths)...
            *penalize_count/max([penalize_count; 1]) + errors/max_error;
        indicators(inf_ind) = inf;
        
        % determine new refinement id
        [~, max_id] = max(indicators);
        refine_id = grid.gid2lid(max_id);
        
        % refine
        nelements_old = grid.nelements;
        refine(grid, refine_id);
        new_ids = nelements_old+1:grid.nelements;
        
        % add data
        interp_func.add_data(C(max_id, :)', V(max_id), interp_func.dilation(end)+10*eps);
        C_added = [C_added; max_id];
        penalize = setdiff(1:size(C, 1), C_added);
        penalize_count(max_id) = 0;
        penalize_count(penalize) = penalize_count(penalize) + 1;
        penalize_count = [penalize_count; zeros(length(new_ids), 1)];
        
    end
    
end

if nargout > 1
    varargout{1}.max_err_seq = max_err_seq;
    varargout{1}.all_points = all_points;
    varargout{1}.all_values = all_values;
    varargout{1}.num_points_seq = num_points_seq;
    if strncmp(interp_func.kernel_type, 'gaussian', 8)
        varargout{1}.dilation_seq = dilation_seq;
    end
    varargout{1}.comp_time = toc(start_algorithm);
end

if params.verbose > 1
    disp(['terminate cubegrid interpolation greedy after ', num2str(toc(start_algorithm)), ...
        ' seconds with ', num2str(interp_func.num_points), ' points.']);
end

end %function

function C = get_midpoints(cube_grid, ids)
% computes midpoints of elements

if nargin < 2
    ids = 1:cube_grid.nelements;
end

n = length(ids);

C = zeros(n, cube_grid.dimension);

for i = 1:cube_grid.dimension
    
    X = reshape(cube_grid.vertex(cube_grid.vertexindex(ids, :), i), ...
        n, size(cube_grid.vertexindex, 2));
    C(:, i) = 0.5^cube_grid.dimension * sum(X, 2);
end

end

function errors = min_max_errors(interp_func, points, data, params)
% if gaussian kernel:
% determines good dilation parameter(s) via heuristic procedure
% always: gives back absolute interpolation errors

errors = test_errors(interp_func, points, data);

if strncmp(interp_func.kernel_type, 'gaussian', 8)
    
    old_max_error = inf;
    
    while max(errors) < old_max_error
        
        old_max_error = max(errors);
        dilation = mean(interp_func.dilation) * ...
            ones(size(interp_func.dilation));
        
        interp_func.dilation = 0.5*dilation;
        errors = test_errors(interp_func, points, data);
        
        if max(errors) > old_max_error
        
            interp_func.dilation = 2*dilation;
            errors = test_errors(interp_func, points, data);
            
            if max(errors) > old_max_error
                
                factors = 0.75:0.125:1.5;
                nf = length(factors);
                errors = zeros(size(points, 1), nf);
                
                for i = 1:nf
                    interp_func.dilation = factors(i)*dilation;
                    errors(:, i) = test_errors(interp_func, points, data);
                end
                
                [~, index] = min(max(errors, [], 1));
                
                interp_func.dilation = factors(index)*dilation;
                errors = errors(:, index);
                
                break
            end
        end
    end
    
end

if strcmp(interp_func.kernel_type, 'gaussian_loc')
    
    dilation = interp_func.dilation;
    
    for cur_ind = [1:interp_func.num_points, 1:interp_func.num_points]
        
        old_max_error = max(errors);
        
        dilation(cur_ind) = 0.5*dilation(cur_ind);
        interp_func.dilation = dilation;
        errors = test_errors(interp_func, points, data);
        
        if max(errors) > old_max_error
            
            dilation(cur_ind) = 4*dilation(cur_ind);
            interp_func.dilation = dilation;
            errors = test_errors(interp_func, points, data);
            
            if max(errors) > old_max_error
                
                dilation(cur_ind) = 0.5*dilation(cur_ind);
                interp_func.dilation = dilation;
                errors = test_errors(interp_func, points, data);
            end
        end
    end
    
end

end

function errors = test_errors(interp_func, points, data)
% computes absolute interpolation errors

errors = inf*ones(size(points, 1), 1);

if interp_func.num_points > 0
    
    for i = 1:length(errors)
        errors(i) = abs(interp_func(points(i, :)') - data(i)) ...
            / abs(data(i));
    end
end

end