
function demo2D_cubvsQMC(domain_number)

%--------------------------------------------------------------------------
% Demo: demo2D_cubvsQMC
%--------------------------------------------------------------------------
% Object:
% Demo on algebraic polynomial cubature on bivariate domains and compares 
% results with those of a QMC rule.
% 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 QMC rule of (XQMC,wc) and an approximation IQMC of IR;
% 6. computes some statistics;
% 7. plots domain and pointsets.
%--------------------------------------------------------------------------
% Dates:
% Written on 02/06/2020: A. Sommariva.
% Modified on 18/06/2020: A. Sommariva.
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% Function to study. 
% See "gallery_2D.m" for details. 
%--------------------------------------------------------------------------

function_example=2;


%--------------------------------------------------------------------------
% 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


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

card_QMC=10^5;


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

[g,gstr]=gallery_2D(function_example);
domain_example=domstr2dom(domain_number);
domain_struct=define_domain(domain_example);

% ........ reference value ........
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;

for k=1:length(nV)
    
    n=nV(k);
    
    % ........ algebraic type 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);
    
    % ........ QMC rule ........

    tic; 
    [XQMC,domain_struct]=cub_QMC(domain_struct,card_QMC);
    cpu(k,2)=toc;
    gXQMC=feval(g,XQMC(:,1),XQMC(:,2)); wQMC=XQMC(:,3);
    IQMC(k)=wQMC'*gXQMC;
    AEQMC(k)=abs(IQMC(k)-IR); REQMC(k)=AEQMC(k)./abs(IR);
    cardXQMC(k)=size(XQMC,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 |  n  | card X | cardXQMC |   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),cardXQMC(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 |  n  |          If           |          IQMC           |   REA   |   REQMC   |  ');
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),IQMC(k),RE(k),REQMC(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 \n');
fprintf('\n \t .............................................................................');
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 IQMC    : 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,XQMC,XW)

















