
function demo2D_cub_ade(domain_number)

%--------------------------------------------------------------------------
% Object:
% Demo on algebraic polynomial cubature on bivariate domains.
% This routine tests the algebraic degree of exactness "ADE" of the rule.
% 1. defines the domain;
% 2. computes a reference rule;
% 3. -> for each degree of exactness "deg"
%     a. computes the moments "momsR" of a suitable polynomial basis of 
%        degree "deg";
%     b. computes a starting rule (X,w) with ADE equal to "deg";
%     c. checks the quality of the approximation;
%     d. computes a compressed rule (XC,wc) from (X,w);
%     e. checks the quality of the approximation;
% 4. computes and writes some statistics;
% 5. plots domain and the last pointsets.
%--------------------------------------------------------------------------
% Dates:
% Written on 29/10/2020: M. Vianello;
%
% Modified on:
% 17/06/2023: A. Sommariva.
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% 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=3; end

% Algebraic degree of exactness in numerical experiments: can be a vector.
nV=2:15;






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

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

% ........ reference value ........
degR=max(nV)+20;
tic; [XWR,dboxR]=define_cub_rule(domain_struct,degR); cpuR=toc;
X=XWR(:,1:2); u=XWR(:,3);

for k=1:length(nV)
    
    n=nV(k);
    dim_poly=(n+1)^2;
    
    % ... reference moments and determine those not too large or small ...
    V=vandermonde_matrix(n,X,dboxR,domain_struct); momsR=V'*u;
    tol=10^(-12); iok=find(momsR >= tol & momsR <= 1/tol);
    
    % ........ full rule ........
    tic; XW=define_cub_rule(domain_struct,n); cpu(k,1)=toc;
    Xf=XW(:,1:2); wf=XW(:,3);
    Vf=vandermonde_matrix(n,Xf,dboxR,domain_struct); momsf=Vf'*wf;
    AE(k)=norm(momsf-momsR);
    RE(k)=norm((momsf(iok)-momsR(iok))./(momsR(iok)));
    cardX(k)=size(XW,1);
    
    % ........ compressed rule ........
    tic; [XC,wc,momerr(k)]=dCATCH(n,Xf,wf,[],[],[]); cpu(k,2)=toc;
    Vc=vandermonde_matrix(n,XC,dboxR,domain_struct); momsc=Vc'*wc;
    AEC(k)=norm(momsc-momsR);
    REC(k)=norm((momsc(iok)-momsR(iok))./(momsR(iok)));
    cardXC(k)=size(XC,1);
    
end

% ..... statistics ....

fprintf('\n \t ..........................................................');
fprintf('\n \t Object: testing algebraic degree of exactness of a rule. \n');
fprintf('\n \t We compute all the moments at degree "n" by a higher degree');
fprintf('\n \t rule and then we test the moments by the rules in analysis');
fprintf('\n \n \t ..........................................................');

fprintf('\n \t Domain: '); disp(domain_struct.domain)


fprintf('\n \n \t ..........................................................');
fprintf('\n \t |  n  | card X | cardXC |   cpuf   |   cpuc   |  momerr  | ');
fprintf('\n \t ..........................................................');
for k=1:length(nV)
    fprintf('\n \t | %3.0f | %6.0f | %6.0f | %1.2e | %1.2e | %1.2e |',...
        nV(k),cardX(k),cardXC(k),cpu(k,1),cpu(k,2),momerr(k));
end
fprintf('\n \t ..........................................................');
fprintf('\n \t | %3.0f | %6.0f |        | %1.2e |          |          |',...
    degR,length(XWR(:,1)),cpuR);
fprintf('\n \t ..........................................................');
fprintf('\n \n');
fprintf('\n \t ..........................................................');
fprintf('\n \t |  n  |   AEf   |   REf   |   AEc   |   REc   |  ');
fprintf('\n \t ..........................................................');
for k=1:length(nV)
    fprintf('\n \t | %3.0f | %1.1e | %1.1e | %1.1e | %1.1e |',...
        nV(k),AE(k),RE(k),AEC(k),REC(k));
end
fprintf('\n \t .......................................................... \n \n');

fprintf('\n \t ...................... Legend ...........................');
fprintf('\n \t n: tested algebraic degree of exactness');
fprintf('\n \t card X  : cardinality of the full rule');
fprintf('\n \t card XC : cardinality of the compressed rule');
fprintf('\n \t cpuf    : cputime for determining full rule');
fprintf('\n \t cpuc    : cputime for determining compressed rule');
fprintf('\n \t momerr  : compression absolute moment error in norm 2');
fprintf('\n \t AEf     : max absolute moment error full rule');
fprintf('\n \t REf     : max relative (filtered) moment error full rule');
fprintf('\n \t AEc     : max absolute moment error compressed rule');
fprintf('\n \t REc     : max relative (filtered) moment error comp. rule');
fprintf('\n \t .........................................................');
fprintf('\n \n');

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

clear_figure(1)
figure(1);
plot_2D(domain_struct,XW,XC);
