
function demo2D_cub_battery

%--------------------------------------------------------------------------
% 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 ADE in a list
%     a. computes a full rule (X,w) at ADE (X,w);
%     b. computes a compressed rule (XC,wc) from (X,w);
%     c. tests the cubature approximation on several integrands ;
% 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

%--------------------------------------------------------------------------
% 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'; exampleV=1:8;


% ADE 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=30;
tic; [XWR,dboxR]=define_cub_rule(domain_struct,degR); cpuR=toc;
Xr=XWR(:,1:2); wr=XWR(:,3);

refV=[]; recV=[];

for k=1:length(nV)
    
    n=nV(k);

    % ........ full rule ........
    XW=define_cub_rule(domain_struct,n); Xf=XW(:,1:2); wf=XW(:,3);
    
    % ........ compressed rule ........
    [Xc,wc,momerr(k)]=dCATCH(n,Xf,wf,[],[],[]);

    for example=exampleV
        [f,fstr]=gallery_2D(example,gallery_type,n);
        Ir(example)=wr'*feval(f,Xr(:,1),Xr(:,2));
        If(example)=wf'*feval(f,Xf(:,1),Xf(:,2));
        Ic(example)=wc'*feval(f,Xc(:,1),Xc(:,2));
    end

    % ....... stats ........
    fprintf('\n \n')
    fprintf('\n \t * degree: %3.0f',n);
    fprintf('\n \t .....................................................');
    for k=1:length(exampleV)
        IrL=Ir(k); IfL=If(k); IcL=Ic(k);
        den=(1-sign(abs(IrL)))+abs(IrL);
        % relative errors
        ref(k)=abs(IrL-IfL)/den;
        rec(k)=abs(IrL-IcL)/den;
        fprintf('\n \t | %2.0f | %1.15e | %1.2e | %1.2e |',...
            exampleV(k),IrL,ref(k),rec(k));
    end
    fprintf('\n \t .....................................................');

    refV=[refV; ref]; recV=[recV; rec]; 

end

fprintf('\n \n');

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

clf;

figure(1);
plot_2D(domain_struct,XW,Xc);

for k=1:length(exampleV)
     legend_arg{k,1}=num2str(exampleV(k));
end

legend_arg=legend_arg(1:length(exampleV));


figure(2)
for k=1:length(exampleV)
    cL=rand(3,1); cL=cL/norm(cL);
    semilogy(nV',refV(:,k),'color',cL,'LineWidth',4);
    hold on;
end

legend(legend_arg);

figure(3)
for k=1:length(exampleV)
    cL=rand(3,1); cL=cL/norm(cL);
    semilogy(nV',recV(:,k),'color',cL,'LineWidth',4);
    hold on;
end

legend(legend_arg);











