
function demoS2_cub_adaptive

%--------------------------------------------------------------------------
% Object:
%--------------------------------------------------------------------------
% Demo on adaptive cubature on spherical polygons.
%--------------------------------------------------------------------------
% Dates:
%--------------------------------------------------------------------------
% Written on 29/10/2020: M. Vianello;
%
% Modified on:
% 16/05/2023: A. Sommariva.
%--------------------------------------------------------------------------

clear;

%--------------------------------------------------------------------------
% gallery_type: string that chooses the function and can be
%
%         'easy' some rather well behaved functions
%                (set as "example" a value in 1,...,23);
%
%         'renka-brown' that are taken from the paper 1., mentioned below;
%                (set as "example" a value in 1,...,10);
%
%         'polynomials' polynomials of degree "deg";
%                (set as "example" a value in 1,...,2);
%            f_type=1: fixed polynomial of degree "deg";
%            f_type=2: random polynomial of degree "deg";
%
% some examples:
%
%          1. Franke function has the setting:
%            gallery_type='easy'; f_type=18;
%
%          2. 'exp(-((x-x0).^2+(y-y0).^2))' has the setting:
%            gallery_type='easy'; f_type=16;
%
%          3 '(pi/16+exp(-3.2)*x+sin(pi/9)*y).^deg' has the setting:
%             gallery_type='polynomials'; f_type=1;
%
% For details, see the file gallery_2D.m in "../EXTERNAL_ROUTINES/DCUB".
%--------------------------------------------------------------------------

gallery_type='hard'; function_example=3;


%--------------------------------------------------------------------------
% Domain to be considered. The variable "domain_str" can be:
%
% case 1, 'spherical-polygon'
% ...
%
% For the definition of these examples in this demo, see the routine
% "define_domain" at the bottom of this file.
% This function is fundamental to understand how to use our codes in
% specific domains.
%--------------------------------------------------------------------------

domain_str='spherical-polygon'; subexample=1;

% ADE in numerical experiments: can be a vector.
nV=2:2:10;




% ........................ Main code below ................................

[g,gstr]=gallery_3D(gallery_type,function_example);
domain_struct=gallery_domains_S2(domain_str,subexample);

% ........ reference value ........
degR=NaN; tol=10^(-14);
tic; [IR,IL,dbox,flag,iters]=cub_adaptive(domain_struct,g,tol); cpuR=toc;

for k=1:length(nV)

    n=nV(k);

    % ........ full rule ........
    tic; [XW,dbox]=define_cub_rule(domain_struct,n); cpu(k,1)=toc;
    gX=feval(g,XW(:,1),XW(:,2),XW(:,3)); w=XW(:,4);
    I(k)=w'*gX;
    AE(k)=abs(I(k)-IR); RE(k)=AE(k)./abs(IR);
    cardX(k)=size(XW,1);

    % ........ compressed rule ........
    X=XW(:,1:3); u=w;
    tic; [XC,wc,momerrL,dbox]=dCATCH(n,X,w); cpu(k,2)=toc;
    gXC=feval(g,XC(:,1),XC(:,2),XC(:,3));
    IC(k)=wc'*gXC;
    AEC(k)=abs(IC(k)-IR); REC(k)=AE(k)./abs(IR);
    cardXC(k)=size(XC,1);

end

% ..... statistics ....
fprintf('\n \t ..............................................................................................');
fprintf('\n \t Domain: '); disp(domain_str);
fprintf('\n \t Function: '); disp(gstr);
fprintf('\n \t ..............................................................................................');
fprintf('\n \t |  n  | card X | cardXC |   cpuf   |   cpuc   |  ');
fprintf('\n \t ..............................................................................................');
for k=1:length(nV)
    fprintf('\n \t | %3.0f | %6.0f | %6.0f | %1.2e | %1.2e |',...
        nV(k),cardX(k),cardXC(k),cpu(k,1),cpu(k,2));
end
fprintf('\n \t ..............................................................................................');
fprintf('\n \t | %3.0f | %6.0f | %6.0f | %1.2e | %1.2e |',...
    NaN,NaN,0,cpuR,0);
fprintf('\n \t ..............................................................................................');
fprintf('\n \n');


fprintf('\n \t ..............................................................................................');
fprintf('\n \t                                RESULTS BY ALGEBRAIC RULES ');
fprintf('\n \t ..............................................................................................');
fprintf('\n \t |  n  |          If           |          Ic           |   REf   |   REc   |  ');
fprintf('\n \t ..............................................................................................');
for k=1:length(nV)
    fprintf('\n \t | %3.0f | %1.15e | %1.15e | %1.1e | %1.1e |',...
        nV(k),I(k),IC(k),REC(k),REC(k));
end
fprintf('\n \t ..............................................................................................');
fprintf('\n \n')

fprintf('\n \t ..............................................................................................');
fprintf('\n \t                               RESULTS BY ADAPTIVE ROUTINE');
fprintf('\n \t ..............................................................................................');
fprintf('\n \t | its |          IH           |          IL           |   AE    |    RE   | flag | ');
fprintf('\n \t ..............................................................................................');
fprintf('\n \t | %3.0f | %1.15e | %1.15e | %1.1e | %1.1e |   %1.0g  |',...
    iters,IR,IL,abs(IR-IL),abs((IR-IL)/(IR+(IR == 0))),flag);
fprintf('\n \t ..............................................................................................');
fprintf('\n \n')
fprintf('\n \t ..............................................................................................');
fprintf('\n \t                                     ADDITIONAL NOTES');
fprintf('\n \t ..............................................................................................');
if flag == 0
fprintf('\n \t * the adaptive routine terminated successfully');
else
    fprintf('\n \t * the adaptive routine did not terminated successfully');
end
fprintf('\n \t * the cubature result provided by adaptive rule is: %1.15e',IR);
fprintf('\n \t * the absolute error estimate is                  : %1.2e',abs((IR-IL)));
fprintf('\n \t * the relative error estimate is                  : %1.2e',abs((IR-IL)/IR));
fprintf('\n \t * the cputime required by the adaptive rule is    : %1.2e',cpuR)
fprintf('\n \t ..............................................................................................');
fprintf('\n \n');




% % ..... plots .....

clear_figure(1)
figure(1);
plot_S2(domain_struct,XW,XC)
hold on;
title_str=['Cubature pointset, ADE: ',num2str(n)];
title(title_str);
axis equal;
hold off;











