%-------------------------------------------------------------------------%
%         This software is licensed by Creative Commons BY-NC-SA:         %
%      http://creativecommons.org/licenses/by-nc-sa/3.0/it/legalcode      %
%-------------------------------------------------------------------------%
%
% File: PUM_3D_CSRBF(dsites,rbf,wf,ep,testfunction,neval,npu)
%
% Goal: script that performs partition of unity interpolation using
%       compactly supported RBFs
%
% Inputs:  dsites:       NX3 matrix representing a set of N data sites
%          rbf:          radial basis function. The CSRBF and the weight 
%                        function used with this code must be given 
%                        in shifted form rbf2(u) = rbf(r), u = 1-e*r. 
%                        For example, the Wendland C2
%                        rbf = @(e,r) max(1-e*r,0).^4.*(4*e*r+1);
%                        becomes rbf2 @(e,r) r.^4.*(5*spones(r)-4*r)
%          wf:           weight function
%          ep:           shape parameter of RBF
%          testfunction: test function
%          neval:        number of evaluation points in one direction
%          npu:          number of PU subdomains in one direction
%  
% Outputs:  epoints: the evaluation points
%           Pf:      the interpolant computed at the evaluation points
%
% Calls on: BlockBased3D_Structure, BlockBased3D_ContainingQuery,
%           BlockBased3D_RangeSearch, BlockBased3D_DistanceMatrix,
%           inhull: by  J. D'Errico, from MATLAB Central File Exchange
%
% Remark:  The original routine comes from the book: 
%          [G. E. Fasshauer, Meshfree approximation methods with Matlab, 
%          World Scientific, Singapore, 2007].
%
%-------------------------------------------------------------------------%
function [epoints Pf] = PUM_3D_CSRBF(dsites,rbf,wf,ep,...
    testfunction,neval,npu)
if ep<0
    disp('ep must be grater than 0')
    return
end
% Define the domain
max_clmns = max(dsites,[],1); min_clmns = min(dsites,[],1);
min_dsites = min(min_clmns); max_dsites = max(max_clmns);
% Create npu-by-npu-by-npu equally spaced centers of PU subdomains 
pugridx = linspace(min_clmns(1),max_clmns(1),npu);
pugridy = linspace(min_clmns(2),max_clmns(2),npu); 
pugridz = linspace(min_clmns(3),max_clmns(3),npu); 
[xpu,ypu,zpu] = meshgrid(pugridx,pugridy,pugridz);
puctrs = [xpu(:) ypu(:) zpu(:)];
% Create neval-by-neval-by-neval equally spaced evaluation points
epgridx = linspace(min_clmns(1),max_clmns(1),neval); 
epgridy = linspace(min_clmns(2),max_clmns(2),neval); 
epgridz = linspace(min_clmns(3),max_clmns(3),neval); 
[xe,ye,ze] = meshgrid(epgridx,epgridy,epgridz);
epoints = [xe(:) ye(:) ze(:)];
% Define the radius of PU subdomains 
puradius = (max_dsites-min_dsites).*sqrt(2)./npu;
wep = 1./puradius;  % Parameter for weight function  
% Parameter for q-by-q-by-q cell-based partitioning structure
q = ceil((max_dsites-min_dsites)./puradius);
K = convhulln(dsites); % Compute the convex hull 
% Reduce points in taking only those contained in the convex hull
in_epoints = inhull(epoints,dsites,K); in_puctrs = inhull(puctrs,dsites,K);
epoints = epoints(in_epoints,:); puctrs = puctrs(in_puctrs,:); 
npu_ConvexHull = size(puctrs,1); neval_ConvexHull = size(epoints,1);
rbfctrs = dsites; 
if npu_ConvexHull == 0
    disp('npu_ConvexHull=0; increase parameter npu')
    return
end
% Compute the testfunction at data sites
rhs = testfunction(dsites(:,1),dsites(:,2),dsites(:,3));
Pf = zeros(neval_ConvexHull,1);  % Initialize
% Compute Shepard evaluation matrix
DM_eval = BlockBased3D_DistanceMatrix(epoints,puctrs,wep);
SEM = wf(ep,DM_eval);
SEM = spdiags(1./(SEM*ones(npu_ConvexHull,1)),0,neval_ConvexHull,...
    neval_ConvexHull)*SEM;
% Build cell-based partitioning structure for data sites and 
% evaluation points
if q>2
idx_ds = BlockBased3D_Structure(dsites,q,puradius,min_dsites); 
idx_ep = BlockBased3D_Structure(epoints,q,puradius,min_dsites); 
end
for j = 1:npu_ConvexHull
    if q>2
        % Find the cell containing the j-th subdomain centre
        index = BlockBased3D_ContainingQuery(puctrs(j,:),q,puradius,...
            min_dsites);
        % Find data sites located in the j-th subdomain
        idx = BlockBased3D_RangeSearch(puctrs(j,:),puradius,...
            dsites(idx_ds{index},:),idx_ds{index});
    else
        idx = BlockBased3D_RangeSearch(puctrs(j,:),puradius,dsites);
    end   
    if (~isempty(idx))
        % Build local interpolation matrix for the j-th subdomain
        DM_data = BlockBased3D_DistanceMatrix(dsites(idx,:),...
            rbfctrs(idx,:),ep);
        IM = rbf(ep,DM_data);
        % Find the evaluation points located in the j-th subdomain
        if q>2
            eidx = BlockBased3D_RangeSearch(puctrs(j,:),puradius,....
                epoints(idx_ep{index},:),idx_ep{index});
        else
            eidx = BlockBased3D_RangeSearch(puctrs(j,:),puradius,epoints);            
        end
        if  (~isempty(eidx))
            % Compute local evaluation matrix and the local RBF interpolant
            DM_eval = BlockBased3D_DistanceMatrix(epoints(eidx,:),...
                rbfctrs(idx,:),ep);
            EM = rbf(ep,DM_eval); localfit = EM * (IM\rhs(idx));
            % Accumulate global fit
            Pf(eidx) = Pf(eidx) + localfit.*SEM(eidx,j);
        end
    end
end
% Compute exact solution
exact = testfunction(epoints(:,1),epoints(:,2),epoints(:,3));
% Compute errors on evaluation grid
maxerr = norm(Pf-exact,inf);
rms_err = norm(Pf-exact)/sqrt(neval_ConvexHull);
fprintf('RMS error:       %e\n', rms_err);
fprintf('Maximum error:   %e\n', maxerr);
