%-------------------------------------------------------------------------%
%         This software is licensed by Creative Commons BY-NC-SA:         %
%      http://creativecommons.org/licenses/by-nc-sa/3.0/it/legalcode      %
%-------------------------------------------------------------------------%
%
% File: BLOOCV_PU(M,dsites,neval,npu,rbf,wf,f,rhs,r_min,h,P1,ep);
%
% Goal: script that performs partition of unity with variable patches
%       and shape parameters
%
% Inputs:  M:            space dimension
%          dsites:       NXM matrix representing a set of N data sites
%          neval:        number of evaluation points in one direction
%          npu:          number of PU subdomains in one direction
%          rbf:          radial basis function
%          wf:           weight function
%          f:            the test function
%          rhs:          the function values
%          r_min:        minimum number of points in each PU subdomain
%          h:            upper bound for the radius
%          P1:           parameter for the partition of the PU radius
%          ep:           the partition for the shape parameter
%
% Outputs:  epoints: the evaluation points
%           Pf:      the interpolant computed at the evaluation points
%
% Calls on: IntegerBased_MD_Structure, IntegerBased_MD_Neighbourhood,
%           IntegerBased_MD_RangeSearch, IntegerBased_MD_ContainingQuery, 
%           DistanceMatrix, MakeSDGrid.
%
% Remark:   DistanceMatrix, MakeSDGrid come from the book:
%           [G. E. Fasshauer, Meshfree approximation methods with Matlab, 
%           World Scientific, Singapore, 2007]. 
%
%-------------------------------------------------------------------------%
function [epoints Pf] = BLOOCV_PU(M,dsites,neval,npu,rbf,wf,f,rhs,...
    r_min,h,P1,ep)
% Create npu^M equally spaced PU centres
puctrs = MakeSDGrid(M,npu);
% Create neval^M equally spaced evaluation points
epoints = MakeSDGrid(M,neval);
puradius = 1./npu;  % Define the initial PU radius
wep = 1./puradius;  % Parameter for weight function
npu_M = size(puctrs,1); neval_M = size(epoints,1); % Initialize;
rbfctrs = dsites; % Define the RBF centres
Pf = zeros(neval_M,1);  % Initialize
% Compute Shepard evaluation matrix
DM_eval = DistanceMatrix(epoints,puctrs);
SEM = wf(wep,DM_eval);
SEM = spdiags(1./(SEM*ones(npu_M,1)),0,neval_M,neval_M)*SEM;
% Parameter for integer-based partitioning structure
q = 1./puradius;
% Build the partitioning structure for data sites and evaluation points
idx_ds = IntegerBased_MD_Structure(dsites,q,puradius,M);
idx_ep = IntegerBased_MD_Structure(epoints,q,puradius,M);
for j = 1:npu_M
    puradius = 1/npu; old_puradius = puradius; % Initialize
    % Find the block containing the j-th subdomain centre
    index1 = IntegerBased_MD_ContainingQuery(puctrs(j,:),q,puradius,M);
    % Find data sites located in the j-th subdomain
    [dxx dx] = IntegerBased_MD_Neighbourhood(dsites,idx_ds,index1,q,M,1);
    idx = IntegerBased_MD_RangeSearch(puctrs(j,:),puradius,dxx,dx);
    % Find a range for the PU radius
    if (length(idx) < r_min)
         t = 1; % Initialize
        while (length(idx) < r_min)
            puradius = puradius + 1/8*old_puradius;
            if puradius > t*old_puradius
                t = t + 1;
                [dxx dx] = IntegerBased_MD_Neighbourhood(dsites,idx_ds,...
                    index1,q,M,t);
            end
            % Find data sites located in the j-th subdomain
            idx = IntegerBased_MD_RangeSearch(puctrs(j,:),puradius,dxx,dx);
        end
    end
    % The BLOOCV scheme
    E = []; % Initialize
    for P = 1:length(ep)
        puradius1 = linspace(puradius,h*puradius,P1);
        for Q = 1:length(puradius1)
            n = 1;
            while n < 99999
                if puradius1(Q) > n*old_puradius
                    n = n + 1;
                else
                    t = n;
                    [dxx dx] = IntegerBased_MD_Neighbourhood(dsites,...
                        idx_ds,index1,q,M,t);
                    break
                end
            end
            % Find data sites located in the j-th subdomain
            idx = IntegerBased_MD_RangeSearch(puctrs(j,:),puradius1(Q),...
                dxx,dx);
            % Build local interpolation matrix for the j-th subdomain
            DM_data = DistanceMatrix(dsites(idx,:),rbfctrs(idx,:));
            % Perform BLOOCV scheme
            Aj = rbf(ep(P),DM_data);
            invAj = inv(Aj);
            % Estimate the error
            e = (invAj*rhs(idx))./diag(invAj);
            E(P,Q) = norm(e(:),inf);
        end
    end
    % Select the radius and the shape parameter
    [e1 idx_e1] = min(E); [e2 idx_e2] = min(e1); 
    ep1 = ep(idx_e1(idx_e2)); puradius2 = puradius1(idx_e2);
    % Find the points lying in the j-th patch of radius puradius2
    n = 1;
    while n < 99999
        if puradius2 > n*old_puradius
            n = n + 1;
        else
            t = n;
            [dxx dx] = IntegerBased_MD_Neighbourhood(dsites,idx_ds,...
                index1,q,M,t);
            break
        end
    end
    [edxx edx] = IntegerBased_MD_Neighbourhood(epoints,idx_ep,index1,q,...
        M,t);
    % Find the points located in the j-th subdomain
    idx = IntegerBased_MD_RangeSearch(puctrs(j,:),puradius2,dxx,dx);
    % Compute the distance matrix
    DM_data = DistanceMatrix(dsites(idx,:),rbfctrs(idx,:));
    % Compute the interpolation matrix
    IM = rbf(ep1,DM_data);
    % Find the evaluation points located in the j-th subdomain
    eidx = IntegerBased_MD_RangeSearch(puctrs(j,:),puradius2,edxx,edx);
    if  (~isempty(eidx))
        % Compute local evaluation matrix and the local RBF interpolant
        DM_eval = DistanceMatrix(epoints(eidx,:),rbfctrs(idx,:));
        EM = rbf(ep1,DM_eval); localfit = EM * (IM\rhs(idx));
        % Accumulate global fit
        Pf(eidx) = Pf(eidx) + localfit.*SEM(eidx,j);
    end
end
% Compute exact solution
exact = f(epoints);
% Compute errors on evaluation grid
maxerr = norm(Pf - exact,inf);
rms_err = norm(Pf - exact)/sqrt(neval_M);
fprintf('RMS error:       %e\n', rms_err);
fprintf('Maximum error:   %e\n', maxerr);