%-------------------------------------------------------------------------%
%
% File: WSVD_PU(dsites,rbf,wf,ep,epw,testfunction,epoints,puctrs,tol)
%
% Goal: script that performs partition of unity interpolation using
%       stable WSVD bases
%
% Inputs:  dsites:       NX2 matrix representing a set of N data sites
%          rbf:          radial basis function. 
%          wf:           weight function
%          ep:           shape parameter of RBF
%          epw:          shape parameter for local weight
%          testfunction: test function
%          epoints:      evaluation points
%          puctrs:       centres of PU subdomains 
%          tol:          tolerance used during the Lanczos procedure
%  
% Outputs:  epoints: the evaluation points
%           Pf:      the interpolant computed at the evaluation points
%
% Calls on: BlockBased2D_Structure, BlockBased2D_ContainingQuery,
%           BlockBased2D_RangeSearch, BlockBased2D_DistanceMatrix,
%           distanceMatrix, by G.E. Fasshauer,
%
% Author:  G.E. Fasshauer. The routine has been modified for a fast
%          computation of the WSVD-PU approximant.
%
%-------------------------------------------------------------------------%
function [epoints Pf] = WSVD_PU(dsites,rbf,wf,ep,epw,testfunction,...
    epoints,puctrs,tol)
if ep<0
    disp('ep must be grater than 0')
    return
end
npu = size(puctrs,1); neval = size(epoints,1);
% Define the radius of PU subdomains 
puradius = sqrt(2/npu);
wep = 1./puradius;  % Parameter for weight function  
% Parameter for q-by-q cell-based partitioning structure
q = ceil(1./puradius);
% Parameter for q-by-q cell-based partitioning structure
% Compute the testfunction at data sites
rhs = testfunction(dsites);
rbfctrs=dsites;
Pf = zeros(neval,1);  % Initialize
% Build cell-based partitioning structure for data sites and
% evaluation points
if q>2
    idx_ds = BlockBased2D_Structure(dsites,q,puradius,0);
    idx_ep = BlockBased2D_Structure(epoints,q,puradius,0);
end
% Compute Shepard evaluation matrix
DM_eval = BlockBased2D_DistanceMatrix(epoints,puctrs,wep);
SEM = wf(epw,DM_eval);
SEM = spdiags(1./(SEM*ones(npu,1)),0,neval,neval)*SEM;
for j = 1:npu
    if q>2
        % Find the square cell containing the j-th subdomain centre
        index = BlockBased2D_ContainingQuery(puctrs(j,:),q,puradius,...
            0);
        % Find data sites located in the j-th subdomain
        idx = BlockBased2D_RangeSearch(puctrs(j,:),puradius,...
            dsites(idx_ds{index},:),idx_ds{index});
    else
        idx = BlockBased2D_RangeSearch(puctrs(j,:),puradius,dsites);
    end   
    if (~isempty(idx))
        % Build local interpolation matrix for the j-th subdomain
        DM_data = distanceMatrix(dsites(idx,:),rbfctrs(idx,:));
        IM = rbf(ep,DM_data); 
        % Find the evaluation points located in the j-th subdomain
        if q>2
            eidx = BlockBased2D_RangeSearch(puctrs(j,:),puradius,....
                epoints(idx_ep{index},:),idx_ep{index});
        else
            eidx = BlockBased2D_RangeSearch(puctrs(j,:),puradius,epoints);            
        end
        if  (~isempty(eidx))
            % Compute local evaluation matrix and the local RBF interpolant
            DM_eval = distanceMatrix(epoints(eidx,:),rbfctrs(idx,:));
            % Compute the stable local approximant
            IM_WSVD = WSVD(IM,rhs(idx),tol,length(idx),rbf,ep);            
            EM = rbf(ep,DM_eval); localfit = EM*IM_WSVD;
            % Accumulate global fit
            Pf(eidx) = Pf(eidx) + localfit.*SEM(eidx,j);
        end
    end
end
% Compute exact solution
exact = testfunction(epoints);
% Compute errors on evaluation grid
maxerr = norm(Pf-exact,inf);
rms_err = norm(Pf-exact)/sqrt(neval);
fprintf('RMS error:       %e\n', rms_err);
fprintf('Maximum error:   %e\n', maxerr);