
function demoS2_cub_ade

%--------------------------------------------------------------------------
% Object:
% Demo on cubature on spherical domains.
% It checks the algebraic degree of exactness of a rule.
%--------------------------------------------------------------------------
% Dates:
% Written on 15/06/2023: A. Sommariva.
%--------------------------------------------------------------------------

clear;

%--------------------------------------------------------------------------
% Function to study. 
% The variable "gallery_str" can be 'easy' or 'hard'.
%--------------------------------------------------------------------------

gallery_str='easy'; example=1;

%--------------------------------------------------------------------------
% Domain to be considered. The variable "domain_type" can be:
%
% 1: 'sphere';
% 2: 'spherical-rectangle';
% 3: 'spherical-triangle'.
% 4: 'spherical-polygon'.
%
% In case domain_type is 2, 3, 4, special examples are considered (for
% details, see the subroutine "define_domain" defined below).
%--------------------------------------------------------------------------

domain_type=4; subexample=1;

%--------------------------------------------------------------------------
% Degrees in numerical experiments: can be a vector.
%--------------------------------------------------------------------------

nV=2:2:10;


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



% ....... Apply settings to define domain, pointsets, and functions .......

% Domain
switch domain_type
    case 1, domain_str='sphere';
    case 2, domain_str='spherical-rectangle';
    case 3, domain_str='spherical-triangle';
    otherwise, domain_str='spherical-polygon';   
end

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

domain_struct=gallery_domains_S2(domain_str,subexample);
[g,gstr]=gallery_3D(gallery_str,example);

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

fprintf('\n \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 \n \t * Computing reference rule \n');

degR=max(nV)+10;
tic; 
[XWR,dboxR]=define_cub_rule(domain_struct,degR); 
gXR=feval(g,XWR(:,1),XWR(:,2),XWR(:,3)); wR=XWR(:,4);
IR=wR'*gXR;
cpuR=toc;


for k=1:length(nV)
    
    n=nV(k);
    dim_poly=(n+1)^2;
    
    fprintf('\n \t * Experiments for degree of precision n= %3.0f',n);

    % ... reference moments and determine those not too large or small ...
    X=XWR(:,1:3); u=XWR(:,4);

    V=vandermonde_matrix(n,X,dboxR',domain_struct,dim_poly); 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:3); wf=XW(:,4);
    Vf=vandermonde_matrix(n,Xf,dboxR',domain_struct,dim_poly); 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,[],[],(n+1)^2); cpu(k,2)=toc;
    Vc=vandermonde_matrix(n,XC,dboxR',domain_struct,dim_poly); 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 \n');

switch domain_type
    case 1
        fprintf('\n \t Domain: sphere');
    case 2
        fprintf('\n \t Domain: spherical-rectangle');
    case 3
        fprintf('\n \t Domain: spherical-triangle');
end

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

fprintf('\n \t * Plotting figure 1: compressed rule at ade %3.0f',...
    nV(end));
clear_figure(1)
figure(1);
title_str='Cubature pointset (compressed rule)';
plot_S2(domain_struct,[],XC,title_str);

fprintf('\n \t * Plotting figure 2: full rule at ade %3.0f',nV(end));
clear_figure(2)
figure(2);
title_str='Cubature pointset (full rule)';
plot_S2(domain_struct,Xf,[],title_str);

fprintf('\n \n')









