
function y=evaluate_scattered(x,t,ft,method,RBF_type,f,domain_structure)

%--------------------------------------------------------------------------
% OBJECT
%--------------------------------------------------------------------------
% Known the samples "(t,ft)", it reconstructs the function values "y" at
% "x".
%--------------------------------------------------------------------------
% INPUT
%--------------------------------------------------------------------------
% x: point in which the function "f" will be approximated;
% t, ft: available samples "(t,ft)" from which "y" that approximates "f(x)"
%    will be computed;
% method: string that decides the method; choose between
%    * 'RBF_SH' taken from G. Fasshauer code: shepard2D.m;
%    * 'DISC' taken from Dell'Accio, Di Tommaso, Siar, Vianello codes;
%    * 'RBF' it evaluates a RBF interpolant with parameter choosen by LOOCV;
%    * 'ScatteredInterpolant' uses the Matlab built-in function
%       "scatteredInterpolant".
%    * 'Multinode' uses the nultinode operator approach.
%    * 'moveint2D' is a modified version of DISC for function evaluations.
%    * 'PUM' is a method based on partition of unity.
%
% RBF_type: parameter that determines the possible RBF to use:
%
% 1: phi=@(r) (1+r.*r).^(1/2);             % Multiquadric
% 2: phi=@(r) exp(-r.*r);                  % Gaussian
% 3: phi=@(r) (1+r.*r).^(-1/2);            % Inverse Multiquadric
% 4: phi=@(r) (1+4*r).*(max(0,(1-r))).^4;  % Wendland 2
% 5: phi=@(r) r.*r.*log(r+eps*(1-sign(r))); % TPS
% 6: phi=@(r) r.^3;                        % polyharmonic spline
% 7: phi=@(r) r.^5;                        % polyharmonic spline
% 8: phi=@(r) r.^7;                        % polyharmonic spline
% 9: phi=@(r) (max(0,(1-r))).^2;           % Wendland W0
% 10: phi=@(r) (35*r.^2+18*r+3).*(max(0,(1-r))).^6;         % Wendland W4
% 11: phi=@(r) (32*r.^3+25*r.^2+8*r+1).*(max(0,(1-r))).^8;  % Wendland W6
%
% f: function to evaluate (not used for computing evaluation).
% domain_structure: structure of the domain containing many informations, 
%      including the bounding box.
%--------------------------------------------------------------------------
% OUTPUT
%--------------------------------------------------------------------------
% y: approximation of "f(x)"
%--------------------------------------------------------------------------
% DATES
%--------------------------------------------------------------------------
% Written on May 27, 2023. (Alvise Sommariva)
% Modified on June 6, 2023 (Alvise Sommariva)
%--------------------------------------------------------------------------
% PAPER
%--------------------------------------------------------------------------
% 1. "Numerical differentiation on scattered data through multivariate
%    polynomial interpolation";
% 2. "DISC: an adaptive numerical Differentiator by local polynomial
%    Interpolation on multivariate SCattered data".
%--------------------------------------------------------------------------
% COPYRIGHT
%--------------------------------------------------------------------------
% Copyright (C) 2023-
%
% Authors:
% Francesco Dell'Accio, Filomena Di Tommaso, Najoua Siar,
% Alvise Sommariva, Marco Vianello.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <https://www.gnu.org/licenses/>.
%--------------------------------------------------------------------------

if nargin < 4, method='DISC'; end
if nargin < 5, RBF_type=10; end

% addpath(genpath('../EXTERNAL_ROUTINES/'));

rbf_str='';

switch method

    case {'mn4','mn6','mn8','mn9'}

        deg_mn=str2num(extractAfter(method,2));

        y=evaluate_scattered_MULTINODE(x,t,ft,deg_mn);


    case 'movint2D'
        y=evaluate_scattered_MOVINT2D(x,t,ft);

    case 'RBF'
        % note: we set warning off (frequant conditioning issues)
        warning off
        [y,rbf_str]=evaluate_scattered_RBF(x,t,ft,RBF_type);
        warning on
        fprintf('\n \t * RBF: '); disp(rbf_str);

    case 'ScatteredInterpolant'
        local_method='natural'; % 'linear', 'natural'
        ExtrapolationMethod='nearest'; % 'nearest', 'linear', or 'none'

        F=scatteredInterpolant(t(:,1),t(:,2),ft,local_method,...
            ExtrapolationMethod);
        y=F(x(:,1),x(:,2));

    case 'RBF_SH'
        y=evaluate_scattered_RBF_SH(x,t,ft);

    case 'DISC'
        y=evaluate_scattered_DISC(x,t,ft);

    case 'PUM'

        % IMPORTANT: Working only in subdomains of [0,1] x [0,1]

        % PUM SETTINGS.

        dbox=domain_structure.dbox;

        % Mapping points to unit square [0,1] x [0,1]
        [tR,dbox]=map2unitsquare(t,dbox);
        [xR,dbox]=map2unitsquare(x,dbox);

        warning off;
        y = BLOOCV_PUM(xR,tR,ft,RBF_type);
        warning on;

end






%--------------------------------------------------------------------------
%                          MULTINODE
%--------------------------------------------------------------------------

function y=evaluate_scattered_MULTINODE(x,t,ft,deg_mn)

%--------------------------------------------------------------------------
% OBJECT
%--------------------------------------------------------------------------
% Known the samples "(t,ft)", it reconstructs the function values "y" at "x"
% (shape parameter equal to "ep"), via multinode operator.
%--------------------------------------------------------------------------

% Multinode evaluation.

if nargin < 4, deg_mn=8; end
y=multinode_operator(x,t,ft,deg_mn);






%--------------------------------------------------------------------------
%                            MOVINT2D
%--------------------------------------------------------------------------

function [y,error_est]=evaluate_scattered_MOVINT2D(x,t,ft)

%--------------------------------------------------------------------------
% OBJECT
%--------------------------------------------------------------------------
% Known the samples "(t,ft)", it reconstructs the function values "y" at "x"
% (shape parameter equal to "ep").
%
% The code is taken from "movint2D" routine by F. Dell'Accio, F. Di Tommaso,
% N. Siar.
%--------------------------------------------------------------------------

% Multinode evaluation.
card_x=size(x,1);

% Compute the values at evaluation points "x".
for k=1:card_x

    xL=x(k,:);
    [yL,error_est]=movint2D(xL,t,ft);

    if isempty(yL)
        warning('movint2D failure: no evaluation, few points');
        fprintf('\n \n');
        yV=NaN;
    end

    y(k,1)=yL(end);

end









%--------------------------------------------------------------------------
%                            RBF + LOOCV
%--------------------------------------------------------------------------

function [y,rbf_str]=evaluate_scattered_RBF(x,t,ft,RBF_type)

%--------------------------------------------------------------------------
% OBJECT
%--------------------------------------------------------------------------
% Known the samples "(t,ft)", it reconstructs the function values "y" at "x".
% Using RBF interpolation.
%--------------------------------------------------------------------------

if nargin < 4, RBF_type=10; end

[rbf,rbf_str] = RBF(RBF_type);

ep=RBF_optimal_scale(t,ft,RBF_type); % scale parameter

N=size(t,1);

% Compute distance matrix between data sites and centers
DM_data = DistanceMatrix(t,t);

% Build collocation matrix
CM = rbf(ep*DM_data);

% Compute RBF coeffs.
coeffs=(CM\ft);

% Compute distance matrix between evaluation points and centers
DM_eval = DistanceMatrix(x,t);
EM = rbf(ep*DM_eval);

% Compute RBF least squares approximation
y = EM * coeffs;




function RBF_scale=RBF_optimal_scale(centers,f_centers,RBF_type,...
    minep,maxep)

%--------------------------------------------------------------------------
% Object:
%--------------------------------------------------------------------------
% As explained in [1], we determine an almost optimal parameter "RBF_scale"
% for RBF interpolation by means of a LOOCV type procedure.
%--------------------------------------------------------------------------
% Input:
%--------------------------------------------------------------------------
%
% centers: RBF centers (N x 2 matrix, where N is the number of centers);
%
% f_centers: values of the function at centers;
%
% RBF_type:
%      1: phi=@(r) (1+r.*r).^(1/2);             % Multiquadric
%      2: phi=@(r) exp(-r.*r);                  % Gaussian
%      3: phi=@(r) (1+r.*r).^(-1/2);            % Inverse Multiquadric
%      4: phi=@(r) (1+4*r).*(max(0,(1-r))).^4;  % Wendland 2
%      5: phi=@(r) r.*r.*log(r+eps*(1-sign(r))); % TPS
%      6: phi=@(r) r.^3;                        % polyharmonic spline
%      7: phi=@(r) r.^5;                        % polyharmonic spline
%      8: phi=@(r) r.^7;                        % polyharmonic spline
%      9: phi=@(r) (max(0,(1-r))).^2;             % Wendland W0
%      10: phi=@(r) (35*r.^2+18*r+3).*(max(0,(1-r))).^6;         % Wendland W4
%      11: phi=@(r) (32*r.^3+25*r.^2+8*r+1).*(max(0,(1-r))).^8;  % Wendland W6
%      12: phi=@(r) (sqrt(2)/(3*sqrt(pi))*(3*r^2*log(r/(1+sqrt(1-r.^2)))+...
%            (2*r^2+1).*sqrt(1-r^2)))*max(0,(1-r));  % Missing Wendland
%      13: phi=@(r) exp(-r);                     % Matern beta_1=3/2.
%      14: phi=@(r) (1+r).*exp(-r);              % Matern beta_2=5/2.
%
%      Note: Optimal scaling "RBF_scale" is sought for RBF_type not in
%            [5,8], otherwise we set "RBF_scale" to 1 (those RBF are scale
%            independent).
%
% minep,maxep: define range for searching a (near) optimal RBF shape
%       parameter, e.g. epsilon in [minep,maxep]=[0.5,15]
%--------------------------------------------------------------------------
% Output:
%--------------------------------------------------------------------------
%
% RBF_scale: almost optimal RBF shape parameter.
%
%--------------------------------------------------------------------------
%% Copyright (C) 2021- R. Cavoretto, A. De Rossi, A. Sommariva, M. Vianello
%%
%% This program is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
%% the Free Software Foundation; either version 2 of the License, or
%% (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; if not, write to the Free Software
%% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
%%
%% Authors:
%%          Roberto Cavoretto <roberto.cavoretto@unito.it>
%%          Alessandra De Rossi   <?alessandra.derossi@unito.it>
%%          Alvise Sommariva <alvise@euler.math.unipd.it>
%%          Marco Vianello   <marcov@euler.math.unipd.it>
%%
%% Last Update: July 9, 2021.
%--------------------------------------------------------------------------

% defaults
if nargin < 3, RBF_type=10; end % Wendland 4.
if nargin < 4, minep=.5; end
if nargin < 5, maxep=15; end

Npts=size(centers,1);

if (RBF_type >= 5 & RBF_type <= 8)
    RBF_scale=1;
else
    % .... define loocv ....
    switch RBF_type
        case 1, RBF_order=1;
        case 5, RBF_order=2;
        case 6, RBF_order=2;
        case 7, RBF_order=3;
        case 8, RBF_order=4;
        otherwise, RBF_order=0;
    end

    % ... adding polynomial components (if required) ...
    switch RBF_order
        case 1
            border_RBF_cubmat=ones(Npts,1);
        case 2
            border_RBF_cubmat=[ones(Npts,1)  centers];
        case 3
            c1=(centers(:,1)); c2=(centers(:,2));
            border_RBF_cubmat=[ones(Npts,1) centers c1.^2 c1.*c2 c2.^2];
        case 4
            c1=(centers(:,1)); c2=(centers(:,2));
            border_RBF_cubmat=[ones(Npts,1) centers c1.^2 c1.*c2 c2.^2 ...
                c1.^3 (c1.^2).*c2 c1.*(c2.^2) c2.^3];
        otherwise
            border_RBF_cubmat=[];
    end

    DM = DistanceMatrix(centers,centers);
    phi = RBF(RBF_type);

    RBF_scale = fminbnd(@(epx) CostEpsilonLOOCV(epx,...
        DM,phi,f_centers,border_RBF_cubmat,RBF_order),...
        minep,maxep);
end




function DM = DistanceMatrix(dsites,ctrs)

%--------------------------------------------------------------------------
% Object:
%--------------------------------------------------------------------------
% DM = DistanceMatrix(dsites,ctrs)
% Forms the distance matrix of two sets of points in R^s,
% i.e., DM(i,j) = || datasite_i - center_j ||_2.
%--------------------------------------------------------------------------
% Input:
%--------------------------------------------------------------------------
%   dsites: Mxs matrix representing a set of M data sites in R^s
%              (i.e., each row contains one s-dimensional point)
%   ctrs:   Nxs matrix representing a set of N centers in R^s
%              (one center per row)
%--------------------------------------------------------------------------
% Output:
%--------------------------------------------------------------------------
%   DM:     MxN matrix whose i,j position contains the Euclidean
%              distance between the i-th data site and j-th center
%--------------------------------------------------------------------------
% Last Update: July 9, 2021.
%--------------------------------------------------------------------------

[M,~] = size(dsites); [N,s] = size(ctrs);
DM = zeros(M,N);

%--------------------------------------------------------------------------
% Accumulate sum of squares of coordinate differences
% The ndgrid command produces two MxN matrices:
%
%   dr, consisting of N identical columns (each containing
%       the d-th coordinate of the M data sites);
%
%   cc, consisting of M identical rows (each containing
%       the d-th coordinate of the N centers).
%--------------------------------------------------------------------------

for d=1:s
    [dr,cc] = ndgrid(dsites(:,d),ctrs(:,d));
    DM = DM + (dr-cc).^2;
end
DM = sqrt(DM);




function [phi,phi_str] =  RBF (RBF_type)

%--------------------------------------------------------------------------
% Input:
%--------------------------------------------------------------------------
% RBF_type: the choosen RBF is the "RBF_type"-th in the list of functions
%     described in "RBF.m".
%
% 1: phi=@(r) (1+r.*r).^(1/2);             % Multiquadric
% 2: phi=@(r) exp(-r.*r);                  % Gaussian
% 3: phi=@(r) (1+r.*r).^(-1/2);            % Inverse Multiquadric
% 4: phi=@(r) (1+4*r).*(max(0,(1-r))).^4;  % Wendland 2
% 5: phi=@(r) r.*r.*log(r+eps*(1-sign(r))); % TPS
% 6: phi=@(r) r.^3;                        % polyharmonic spline
% 7: phi=@(r) r.^5;                        % polyharmonic spline
% 8: phi=@(r) r.^7;                        % polyharmonic spline
% 9: phi=@(r) (max(0,(1-r))).^2;             % Wendland W0
% 10: phi=@(r) (35*r.^2+18*r+3).*(max(0,(1-r))).^6;         % Wendland W4
% 11: phi=@(r) (32*r.^3+25*r.^2+8*r+1).*(max(0,(1-r))).^8;  % Wendland W6
% 12: phi=@(r) (sqrt(2)/(3*sqrt(pi))*(3*r^2*log(r/(1+sqrt(1-r.^2)))+...
%            (2*r^2+1).*sqrt(1-r^2)))*max(0,(1-r));  % Missing Wendland
% 13 phi=@(r) exp(-r);                     % Matern beta_1=3/2.
% 14 phi=@(r) (1+r).*exp(-r);              % Matern beta_2=5/2.
% 15: phi=@(r) (max(0,(1-r)))^10*(429*r^4 + 450*r^3 + 210*r^2 + 50*r + 5)
%--------------------------------------------------------------------------
% Output:
%--------------------------------------------------------------------------
% phi: RBF function.
% phi_str: string with the RBF "name".
%--------------------------------------------------------------------------
% Copyrights.
%--------------------------------------------------------------------------
%% Copyright (C) 2007-2019 Alvise Sommariva, Marco Vianello.
%%
%% This program is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
%% the Free Software Foundation; either version 2 of the License, or
%% (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program; if not, write to the Free Software
%% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
%%
%% Author:  Alvise Sommariva <alvise@euler.math.unipd.it>
%%          Marco Vianello   <marcov@euler.math.unipd.it>
%%
%% Date: February 14, 2020.
%--------------------------------------------------------------------------

switch RBF_type
    case 1
        phi_str='Multiquadric';
        phi=@(r) (1+r.*r).^(1/2);            % Multiquadric
    case 2
        phi_str='Gaussian';
        phi=@(r) exp(-r.*r);                  % Gaussian
    case 3
        phi_str='Inverse Multiquadric';
        phi=@(r) (1+r.*r).^(-1/2);            % Inverse Multiquadric
    case 4
        phi_str='Wendland 2';
        phi=@(r) (1+4*r).*(max(0,(1-r))).^4;  % Wendland 2
    case 5
        phi_str='TPS';
        pert=@(r) eps*(1-abs(sign(r)));
        phi=@(r) r.*r.*log(r+pert(r)); % TPS
    case 6
        phi_str='PH: r^3';
        phi=@(r) r.^3;                        % polyharmonic spline
    case 7
        phi_str='PH: r^5';
        phi=@(r) r.^5;                        % polyharmonic spline
    case 8
        phi_str='PH: r^7';
        phi=@(r) r.^7;                        % polyharmonic spline
    case 9
        phi_str='Wendland W0';
        phi=@(r) (max(0,(1-r))).^2;             % Wendland W0
    case 10
        phi_str='Wendland W4';
        phi=@(r) (35*r.^2+18*r+3).*(max(0,(1-r))).^6;  % Wendland W4
    case 11
        phi_str='Wendland W6';
        phi=@(r) (32*r.^3+25*r.^2+8*r+1).*(max(0,(1-r))).^8;  % Wendland W6
    case 12                                    % Missing Wendland
        phi_str='Missing Wendland';
        pert=@(r) eps*(1-abs(sign(r)));
        phi=@(r) (sqrt(2)/(3*sqrt(pi))*...
            (3*r.^2.*log( r./(1+sqrt(1-r.^2)) + pert(r) )+...
            (2*r.^2+1).*sqrt(1-r.^2))).*max(0,(1-r));
    case 13                             % Matern beta_1=(d+1)/2, where d=2.
        phi_str='Matern beta_1';
        phi=@(r) exp(-r);
    case 14                             % Matern beta_2=(d+3)/2, where d=2.
        phi_str='Matern beta_2';
        phi=@(r) (1+r).*exp(-r);
    case 15
        phi_str='Wendland W8';
        phi=@(r) (max(0,(1-r))).^10.*(429*r.^4 + 450*r.^3 + 210*r.^2 + ...
            50*r + 5);

    otherwise
        error('RBF type not implemented')
end




function ceps = CostEpsilonLOOCV(ep,r,phi,rhs,P,RBF_order)

%--------------------------------------------------------------------------
% Object:
%--------------------------------------------------------------------------
% Implements cost function for optimization of shape parameter via Rippa's
% LOOCV algorithm.
%--------------------------------------------------------------------------
% Last Update: July 9, 2021.
%--------------------------------------------------------------------------

A = phi(r.*ep);
if RBF_order >= 1
    zeromat = zeros(size(P,2));
    A = [A P; ...
        P' zeromat];
    rhs = [rhs; zeros(size(zeromat,1),1)];
end
invA = pinv(A);
EF = (invA*rhs)./diag(invA);
EF = EF(1:end-RBF_order);
ceps = norm(EF(:),inf);









%--------------------------------------------------------------------------
%                        RBF SHEPARD
%--------------------------------------------------------------------------

function y=evaluate_scattered_RBF_SH(x,t,ft,ep)

%--------------------------------------------------------------------------
% OBJECT
%--------------------------------------------------------------------------
% Known the samples "(t,ft)", it reconstructs the function values "y" at "x"
% (shape parameter equal to "ep").
%
% The code is taken from "shepard2D.m", by Greg Fasshauer.
%--------------------------------------------------------------------------

if nargin <= 3, ep=100; end

rbf = @(e,r) exp(-(e*r).^2);

N=size(t,1);
DM_eval = DistanceMatrix(x,t);
EM = rbf(ep,DM_eval);
EM = EM./repmat(EM*ones(N,1),1,N);
y = EM*ft;









%--------------------------------------------------------------------------
%                              DISC
%--------------------------------------------------------------------------

function y=evaluate_scattered_DISC(x,t,ft)

%--------------------------------------------------------------------------
% OBJECT
%--------------------------------------------------------------------------
% Known the samples "(t,ft)", it reconstructs the function values "y" at "x"
% (shape parameter equal to "ep").
%
% The code is taken from DISC routine by F. Dell'Accio, F. Di Tommaso,
% N. Siar, M. Vianello.
%--------------------------------------------------------------------------


% Number of points in which evaluate the function.
card_x=size(x,1);

% Compute the values at evaluation points "x".
for k=1:card_x

    xL=x(k,:);
    [yV,error]=derivative_nu_f_xb(t,ft,xL);

    if isempty(yV)
        warning('evaluation not performed');
        yV=NaN;
    end

    y(k,1)=yV(length(yV));
    error(error==0)=eps;

end






function [ptsR,dbox]=map2unitsquare(pts,dbox)

%--------------------------------------------------------------------------
% OBJECT
%--------------------------------------------------------------------------
% This routine maps the points "pts" containted in the rectangle defined by
% dbox=[xLimit; yLimit], to "ptsR" belonging to the reference square
% [0,1]^2.
%--------------------------------------------------------------------------

X=pts(:,1); Y=pts(:,2);

if nargin < 2
    xLimit=[min(X) max(X)]; yLimit=[min(Y) max(Y)]; dbox=[xLimit; yLimit];
end

xLimit=dbox(1,:);
a=xLimit(1); b=xLimit(2);  h=(b-a);
XR=(X-a)/h;

yLimit=dbox(2,:);
a=yLimit(1); b=yLimit(2);  h=(b-a);
YR=(Y-a)/h;

ptsR=[XR YR];








