
function demo2D_cub_adaptive(domain_number)

%--------------------------------------------------------------------------
% Object:
%--------------------------------------------------------------------------
% 1. Demo on algebraic polynomial cubature on some bivariate domains.
% 2. At the moment it works only over polygonal domains.
% 3. It compares the results with some algebraic rules.
%--------------------------------------------------------------------------
% Dates:
%--------------------------------------------------------------------------
% Written on 29/10/2020: M. Vianello;
%
% Modified on:
% 10/07/2023: A. Sommariva.
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% 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='renka-brown'; function_example=3;


%--------------------------------------------------------------------------
% Domain to be considered. The variable "domain_number" can be:
%
%     case 1, domain_example='polygon';
%    * case 2, domain_example='disk';
%    * case 3, domain_example='lune';
%    * case 4, domain_example='circular-annular-sector';
%    * case 5, domain_example='sector';
%    * case 6, domain_example='asymmetric-circular-sector';
%    * case 7, domain_example='asymmetric-annulus';
%    * case 8, domain_example='vertical-circular-zone';
%    * case 9, domain_example='horizontal-circular-zone';
%    * case 10, domain_example='circular-segment';
%    * case 11, domain_example='symmetric-lens';
%    * case 12, domain_example='butterfly';
%    * case 13, domain_example='candy';
%    * case 14, domain_example='NURBS';
%    * case 15, domain_example='union-disks';
%    * case 16, domain_example='asymmetric-circular-sector';
%     case 17, domain_example='square';
%     case 18, domain_example='rectangle';
%     case 19, domain_example='triangle';
%    * case 20, domain_example='polygcirc';
%     case 21, domain_example='unit-simplex';
%     case 22, domain_example='unit-square[0,1]x[0,1]';
%
% NOTE: The domains with "*" have not adaptive cubature routines.
%--------------------------------------------------------------------------

if nargin < 1, domain_number=1; end

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




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

switch domain_number
    case {1,17,18,19,21,22}
        [g,gstr]=gallery_2D(function_example,gallery_type);
    otherwise
        warning('Domain not available, proposing a polygon.');
        domain_number=1;
        [g,gstr]=gallery_2D(function_example,gallery_type);
end

domain_str=domstr2dom(domain_number);
domain_struct=define_domain(domain_str);


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

% ......... comparisons ...........
for k=1:length(nV)

    n=nV(k);

    fprintf('\n \t \t -> ADE: %3.0f',n);

    % ........ full rule ........
    tic; [XW,dbox]=define_cub_rule(domain_struct,n); cpu(k,1)=toc;
    gX=feval(g,XW(:,1),XW(:,2)); w=XW(:,3);
    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:2); u=w;
    tic; [XC,wc,momerrL,dbox]=dCATCH(n,X,w); cpu(k,2)=toc;
    gXC=feval(g,XC(:,1),XC(:,2));
    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_2D(domain_struct,XW,XC)
hold on;
title_str=['Cubature pointset, ADE: ',num2str(n)];
title(title_str);
axis equal;
hold off;











