
function demo2D_cub(domain_number)

%--------------------------------------------------------------------------
% Object:
% Demo on algebraic polynomial cubature on bivariate domains.
% The demo:
% 1. defines the test function;
% 2. defines the domain;
% 3. computes a reference value IR of the integral;
% 4. computes a starting rule (X,w) and an approximation If of IR;
% 5. computes a compressed rule of (XC,wc) and an approximation Ic of IR;
% 6. computes some statistics;
% 7. plots domain and pointsets.
%--------------------------------------------------------------------------
% DATES:
%--------------------------------------------------------------------------
% Written on 29/10/2020: A. Sommariva and M. Vianello;
%
% Modified on:
% 17/06/2023: A. Sommariva.
%--------------------------------------------------------------------------
% COPYRIGHT
%--------------------------------------------------------------------------
% Copyright (C) 2023-
%
% Authors:
% 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/>.
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% 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]';
%
%--------------------------------------------------------------------------

if nargin < 1, domain_number=14; end

%--------------------------------------------------------------------------
% 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='easy'; function_example=16;


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


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


% Domain from a gallery
domain_example=domstr2dom(domain_number);
domain_struct=define_domain(domain_example);

% Integrand from a gallery
[g,gstr]=gallery_2D(function_example,gallery_type,nV);

% ........ reference value ........

fprintf('\n \t * computing reference rule');
degR=max(nV)+10;
tic; [XWR,dbox]=define_cub_rule(domain_struct,degR); cpuR=toc;
gXR=feval(g,XWR(:,1),XWR(:,2)); wR=XWR(:,3);
IR=wR'*gXR;

fprintf('\n \t * computing test rules');
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)); 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;
    if size(X,1) > (n+1)*(n+2)/2
        tic; [XC,wc,momerrL,dbox]=dCATCH(n,X,w); cpu(k,2)=toc;
    else
        XC=XW; wc=w; cpu(k,2)=0;
    end
    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_example);
fprintf('\n \t Function: '); disp(gstr);
fprintf('\n \t .............................................................................');
fprintf('\n \t | ade | 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 |',...
    degR,length(XWR(:,1)),0,cpuR,0);
fprintf('\n \t .............................................................................');
fprintf('\n \n');

fprintf('\n \t .............................................................................');
fprintf('\n \t | ade |          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 \t | %3.0f | %1.15e | %1.15e | %1.1e | %1.1e |',...
    degR,IR,IR,0,0);
fprintf('\n \t .............................................................................');
fprintf('\n \n');
fprintf('\n \t .............................................................................');
fprintf('\n \t ade   : algebraic degree of exactness');
fprintf('\n \t card X: cardinality full rule');
fprintf('\n \t cardXC: cardinality compressed rule');
fprintf('\n \t cpuf  : cputime full rule');
fprintf('\n \t cpuf  : cputime compressed rule');
fprintf('\n \t If    : integral computed by full rule');
fprintf('\n \t Ic    : integral computed by compressed rule');
fprintf('\n \t REf   : relative error of integral computed by full rule');
fprintf('\n \t REc   : relative error of integral computed by comp. rule');
fprintf('\n \t .............................................................................');
fprintf('\n \n');




% ..... plots .....
clear_figure(1)
figure(1);
plot_2D(domain_struct,XW,XC)

clear_figure(2)
figure(2)
plot_2D(domain_struct)





















