
function demo3D_cub

%--------------------------------------------------------------------------
% Object:
% Demo on using "dPOLYFIT", "dPOLYVAL", "dLEB" for algebraic polynomial
% interpolation/least squares approximation on spherical domains.
% The determination of the X pointset varies with the choice of the
% parameter "pts_type".
%--------------------------------------------------------------------------
% Dates:
% Written on 29/10/2020: M. Vianello;
%
% Modified on:
% 13/11/2020: A. Sommariva.
%--------------------------------------------------------------------------

clear;

%--------------------------------------------------------------------------
% Function to study. The variable "example" can be:
%
%  case 1 f=@(x,y,z) 1./(1+25*(x+2*y+z).^2);
%  case 2 f=@(x,y,z) sin(x+2*y+1.5*z);
%  case 3 f=@(x,y,z) exp(x+2*y+1.5*z).*cos(10*(x+2*y+1.5*z)).*...
%             tanh(4*(x+2*y+1.5*z));
%  case 4 f=@(x,y,z) exp(x+2*y+1.5*z);
%  case 5 f=@(x,y,z) exp(-(x+2*y+1.5*z));
%  case 6 f=@(x,y,z) (1.2*(x+2*y+1.5*z).^5+3*(x+2*y+1.5*z)+...
%             0.02*(x+2*y+1.5*z).^4-2).^4;
%  case 7 f=@(x,y,z) (x+2*y+1.5*z).^20+1;
%  case 8 f=@(x,y,z) (1.2*(x+2*y+1.5*z).^3+0.02*(x+2*y+1.5*z).^4+...
%             3*(x+2*y+1.5*z)-2).^2;
%--------------------------------------------------------------------------

example=2;


%--------------------------------------------------------------------------
% Domain to be considered. The variable "domain_example" can be:
%
% domain_example='pyramid';
% domain_example='pyramid-boundary';
% domain_example='3D_rect';
% domain_example='rotation-domain';
%--------------------------------------------------------------------------

% domain_example='pyramid';
% domain_example='pyramid-boundary';
% domain_example='3D_rect';
% domain_example='rotation-domain';
domain_example='rotation-surface';


%--------------------------------------------------------------------------
% "pts_type" can be one of the following strings:
%
% 'tchakaloff-near-optimal', 'afek', 'tchakaloff-compression',
% 'standard-least-squares'
%--------------------------------------------------------------------------
% pts_type='tchakaloff-near-optimal';
% pts_type='afek';
pts_type='tchakaloff-compression';
% pts_type='standard-least-squares';

% Degrees in numerical experiments: can be a vector.
nV=8;

% plot final domain and pointsets
plot_domain=1;

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

[g,gs]=define_function(example);
domain_struct=define_domain(domain_example);

fprintf('\n \t ');

for k=1:length(nV)
    n=nV(k);
    fprintf('%3.0f :',n);
    [X,dbox]=define_wam_pointset(domain_struct,n);
    g_X=feval(g,X(:,1),X(:,2),X(:,3));

    [pts,weights,jvec]=define_source_pointset(pts_type,n,X,dbox,...
        domain_struct);
    g_pts=feval(g,pts(:,1),pts(:,2),pts(:,3));

    [coeff,R,jvec] = dPOLYFIT(n,pts,weights,g_pts,jvec,dbox,...
        domain_struct);
    p_X=dPOLYVAL(n,coeff,X,R,jvec,dbox,domain_struct);

    err(k)=norm(p_X-g_X)/norm(g_X);

    leb(k)=dLEB(n,pts,weights,X,jvec,dbox,domain_struct);

    card(k)=size(pts,1);
    cardX(k)=size(X,1);
end

% ..... statistics .....

fprintf('\n \t ........................................................');
fprintf('\n \t Domain: '); disp(domain_example);
fprintf('\n \t Function: '); disp(gs);
fprintf('\n \t ........................................................');
fprintf('\n \t |  n  |   err   |   leb   |   card    |   mesh0   |');
fprintf('\n \t ........................................................');
for k=1:length(nV)
    fprintf('\n \t | %3.0f | %1.1e | %1.1e |   %5.0f   |   %5.0f   |',...
        nV(k),err(k),leb(k),card(k),cardX(k));
end
fprintf('\n \t ........................................................');
fprintf('\n \n');

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


clf;

figure(1);
plot(nV,leb,'-o'); hold on; title('Lebesgue constant'); hold off;

figure(2);
semilogy(nV,err); hold on; title('Approximation errors'); hold off;

if plot_domain
    figure(3)
    switch domain_struct.domain
        case 'pyramid'
            structure_basis=domain_struct.parms{1};
            V=domain_struct.parms{2};
            elev=domain_struct.parms{3};
            plot_pyramid(structure_basis,V,elev,X,pts);

        case '3D_rect'
            dbox=domain_struct.parms{1};
            plot_3D_rect(dbox,X);

        case {'rotation-domain','rotation-surface'}
            basis_structure=domain_struct.parms{1};
            arc_intv=domain_struct.parms{2};
            dir_axis=domain_struct.parms{3};
            x0=domain_struct.parms{4};
            plot_rotation_domain(basis_structure,arc_intv,dir_axis,...
                x0,pts);

        case 'pyramid-boundary'
            structure_basis.domain='pyramid';
            structure_basis=domain_struct.parms{1};
            V=domain_struct.parms{2};
            elev=domain_struct.parms{3};
            plot_pyramid(structure_basis,V,elev,X,pts);

        otherwise
            plot3(X(:,1),X(:,2),X(:,3),'*','Color','red');
            hold on;
            plot3(pts(:,1),pts(:,2),pts(:,3),'*','Color','green');
    end
end








    function [pts,weights,jvec,dbox]=define_source_pointset(pts_type,n,X,...
            dbox,domain_struct)

        if nargin < 3, dbox=[]; end

        switch pts_type
            case 'afek'
                %% APPROXIMATE FEKETE INTERPOLATION POINTS
                [pts,jvec,dbox] = dAPPROXFEK(n,X,0,dbox,domain_struct);
                weights=ones(length(pts(:,1)),1);

            case 'leja'
                %% DISCRETE LEJA INTERPOLATION POINTS
                [pts,jvec,dbox] = dAPPROXFEK(n,X,1,dbox,domain_struct);
                weights=ones(length(pts(:,1)),1);

            case 'tchakaloff-near-optimal'
                %% TCHAKALOFF NEAR-OPTIMAL REGRESSION POINTS
                [pts,weights,~,~,dbox]=dNORD(n,X,domain_struct);
                jvec=[];

            case 'tchakaloff-compression'
                %% TCHKALOFF COMPRESSION OF STANDARD LEAST SQUARES
                [pts,weights,~,dbox]=dCATCH(2*n,X,domain_struct);
                jvec=[];

            case 'standard-least-squares'
                %% STANDARD LEAST SQUARES
                pts=X;
                weights=ones(length(X(:,1)),1);
                jvec=[];
                dbox=[];


            otherwise
                warning('Pointset not implemented, using Approx. Fekete points');
                %% APPROXIMATE FEKETE INTERPOLATION POINTS
                [pts,jvec,dbox] = dAPPROXFEK(n,X,0,dbox,domain_struct);
                weights=ones(length(fek(:,1)),1);

        end

        if isempty(dbox)
            dbox=[];
            L=size(X,2);
            for k=1:L
                mx=min(X(:,k)); Mx=max(X(:,k)); dboxL=[mx; Mx]; dbox=[dbox dboxL];
            end
        end

    end % end: "define_source_pointset"
end








function [f,fs]=define_function(example)

switch example
    case 1
        f=@(x,y,z) 1./(1+25*(x+2*y+z).^2);
        fs='1./(1+25*(x+2*y+z).^2)';
    case 2
        f=@(x,y,z) sin(x+2*y+1.5*z);
        fs='sin(x+2*y+1.5*z)';
    case 3
        f=@(x,y,z) exp(x+2*y+1.5*z).*cos(10*(x+2*y+1.5*z)).*...
            tanh(4*(x+2*y+1.5*z));
        fs='exp(x+2*y+1.5*z)*cos(10*(x+2*y+1.5*z))*tanh(4*(x+2*y+1.5*z));';
    case 4
        f=@(x,y,z) exp(x+2*y+1.5*z);
        fs='exp(x+2*y+1.5*z)';
    case 5
        f=@(x,y,z) exp(-(x+2*y+1.5*z));
        fs='exp(-(x+2*y+1.5*z))';
    case 6
        f=@(x,y,z) (1.2*(x+2*y+1.5*z).^5+3*(x+2*y+1.5*z)+0.02*(x+2*y+1.5*z).^4-2).^4;
        fs='(1.2*(x+2*y+1.5*z).^5+3*(x+2*y+1.5*z)+0.02*(x+2*y+1.5*z).^4-2).^4';
    case 7
        f=@(x,y,z) (x+2*y+1.5*z).^20+1;
        fs='(x+2*y+1.5*z).^20+1';
    case 8
        f=@(x,y,z) (1.2*(x+2*y+1.5*z).^3+0.02*(x+2*y+1.5*z).^4+3*(x+2*y+1.5*z)-2).^2;
        fs='(1.2*(x+2*y+1.5*z).^3+0.02*(x+2*y+1.5*z).^4+3*(x+2*y+1.5*z)-2).^2';
end

end % end: "define_function"









function [domain_struct]=define_domain(example)

domain_struct.dbox=[];
domain_struct.domain=example;

switch example

    case 'pyramid'
        % P: column array of 2D mesh points (cone basis)
        % V: vertex (cone vertex)
        % elev: elevation of the cut in percentage

        % basis_struct.domain='disk';
        % basis_struct.domain='polygon';
        basis_struct.domain='polygonal-boundary';

        switch basis_struct.domain
            case 'disk' % cone basis (disk)
                x1=1; x2=1.5; r1=0.5; % center and radius
                basis_struct.parms=[x1 x2 r1];

            case {'polygon','polygonal-boundary'} % cone basis (polygon)
                basis_struct.parms=[...
                    0.1 0; 0.7 0.2; 1 0.5; 0.75 0.85; 0.5 1; 0 0.25; 0.1 0];

        end

        V=[0 0 1];                 % vertex
        elev=1;                    % elevation

        domain_struct.parms={basis_struct,V,elev};


    case 'pyramid-boundary'
        % P: column array of 2D mesh points (cone basis)
        % V: vertex (cone vertex)
        % elev: elevation of the cut in percentage

        % basis_struct.domain='circle-arc';
        basis_struct.domain='polygonal-boundary';
        % basis_struct.domain='asymmetric-circular-sector-boundary';
        % basis_struct.domain='asymmetric-annulus-boundary';
        % basis_struct.domain='vertical-circular-zone-boundary';
        % basis_struct.domain='horizontal-circular-zone-boundary';
        % basis_struct.domain='circular-segment-boundary';
        % basis_struct.domain='symmetric-lens-boundary';
        % basis_struct.domain='butterfly-boundary';
        % basis_struct.domain='candy-boundary';
        % basis_struct.domain='sector-boundary;

        switch basis_struct.domain
            case 'circle-arc' % cone basis (disk)
                x1=1; x2=1.5; r1=0.5; % center and radius
                arc_angles=[0 pi/4];
                basis_struct.parms=[x1 x2 r1 arc_angles];

            case 'polygonal-boundary' % cone basis (polygon)

                basis_struct.parms=[0.1 0; 0.7 0.2; 1 0.5; 0.75 0.85; ...
                    0.5 1; 0 0.25; 0.1 0];


            case 'asymmetric-circular-sector-boundary'
                % in this case "domain_struct.parms" is a vector [a b c d r]
                a=1; b=1; c=2; d=2; r=3;
                basis_struct.parms=[a b c d r];

            case 'asymmetric-annulus-boundary'
                % in this case "domain_struct.parms" is a vector [a b c d r1 r2]
                a=1; b=1; c=2; d=2; r1=3; r2=4;
                basis_struct.parms=[a b c d r1 r2];

            case 'vertical-circular-zone-boundary'
                % in this case "domain_struct.parms" is a vector [a b r alpha]
                a=1; b=1; r=4; alpha=-pi/4;
                basis_struct.parms=[a b r alpha];

            case 'horizontal-circular-zone-boundary'
                % in this case "domain_struct.parms" is a vector [a b r alpha]
                a=1; b=1; r=4; alpha=-pi/4;
                basis_struct.parms=[a b r alpha];

            case 'circular-segment-boundary'
                % in this case "domain_parms" is a vector [a b r alpha beta]
                a=1; b=1; r=4; alpha=-pi/4; beta=pi/6;
                basis_struct.parms=[a b r alpha beta];

            case 'symmetric-lens-boundary'
                % in this case "domain_parms" is a vector [a r] (with a < r)
                a=1; r=3;
                basis_struct.parms=[a r];

            case 'butterfly-boundary'
                % in this case "domain_parms" is a vector [a b r alpha beta]
                a=1; b=2; r=3; alpha=-pi/4; beta=pi/6;
                basis_struct.parms=[a b r alpha beta];

            case 'candy-boundary'
                % in this case "domain_parms" is a vector [a r alpha] with
                % "-alpha > acos(a/r)"
                a=0.5; r=1; alpha=-1.5;
                basis_struct.parms=[a r alpha];

            case 'sector-boundary'
                center=[1 1.5]; r=0.5; alpha=pi/4; beta=pi/2;
                basis_struct.parms=[center r alpha beta];

        end

        V=[0 0 1];                 % vertex
        elev=0.7;                    % elevation

        domain_struct.parms={basis_struct,V,elev};








    case '3D_rect'

        dBOX=[-1 0.5; -1 0.5; -1 0.5]';
        domain_struct.parms={dBOX};







    case 'rotation-domain'


        % basis_domain.domain='polygonal-boundary';
        % basis_domain.domain='disk';
        % basis_domain.domain='circle-arc';
        % basis_domain.domain='polygon';
        % basis_domain.domain='asymmetric-circular-sector-boundary';
        % basis_domain.domain='asymmetric-annulus-boundary';
        % basis_domain.domain='vertical-circular-zone-boundary';
        % basis_domain.domain='horizontal-circular-zone-boundary';
        % basis_domain.domain='circular-segment-boundary';
        % basis_domain.domain='symmetric-lens-boundary';
        % basis_domain.domain='butterfly-boundary';
        % basis_domain.domain='candy-boundary';
        % basis_domain.domain='sector-boundary;
        % basis_domain.domain='asymmetric-circular-sector';
        % basis_domain.domain='asymmetric-annulus';
        % basis_domain.domain='vertical-circular-zone';
        % basis_domain.domain='horizontal-circular-zone';
        % basis_domain.domain='circular-segment';
        % basis_domain.domain='symmetric-lens';
        % basis_domain.domain='butterfly';
        % basis_domain.domain='candy';
        % basis_domain.domain='sector;

        switch basis_domain.domain

            case 'asymmetric-circular-sector'
                % in this case "domain_struct.parms" is a vector [a b c d r]
                a=1; b=1; c=2; d=2; r=3;
                basis_domain.parms=[a b c d r];

            case 'asymmetric-annulus'
                % in this case "domain_struct.parms" is a vector [a b c d r1 r2]
                a=1; b=1; c=2; d=2; r1=3; r2=4;
                basis_domain.parms=[a b c d r1 r2];

            case 'vertical-circular-zone'
                % in this case "domain_struct.parms" is a vector [a b r alpha]
                a=1; b=1; r=4; alpha=-pi/4;
                basis_domain.parms=[a b r alpha];

            case 'horizontal-circular-zone'
                % in this case "domain_struct.parms" is a vector [a b r alpha]
                a=1; b=1; r=4; alpha=-pi/4;
                basis_domain.parms=[a b r alpha];

            case 'circular-segment'
                % in this case "domain_parms" is a vector [a b r alpha beta]
                a=1; b=1; r=4; alpha=-pi/4; beta=pi/6;
                basis_domain.parms=[a b r alpha beta];

            case 'symmetric-lens'
                % in this case "domain_parms" is a vector [a r] (with a < r)
                a=1; r=3;
                basis_domain.parms=[a r];

            case 'butterfly'
                % in this case "domain_parms" is a vector [a b r alpha beta]
                a=1; b=2; r=3; alpha=-pi/4; beta=pi/6;
                basis_domain.parms=[a b r alpha beta];

            case 'candy'
                % in this case "domain_parms" is a vector [a r alpha] with
                % "-alpha > acos(a/r)"
                a=0.5; r=1; alpha=-1.5;
                basis_domain.parms=[a r alpha];

            case 'sector'
                center=[1 1.5]; r=0.5; alpha=pi/4; beta=pi/2;
                basis_domain.parms=[center r alpha beta];


            case 'asymmetric-circular-sector-boundary'
                % in this case "domain_struct.parms" is a vector [a b c d r]
                a=1; b=1; c=2; d=2; r=3;
                basis_domain.parms=[a b c d r];

            case 'asymmetric-annulus-boundary'
                % in this case "domain_struct.parms" is a vector [a b c d r1 r2]
                a=1; b=1; c=2; d=2; r1=3; r2=4;
                basis_domain.parms=[a b c d r1 r2];

            case 'vertical-circular-zone-boundary'
                % in this case "domain_struct.parms" is a vector [a b r alpha]
                a=1; b=1; r=4; alpha=-pi/4;
                basis_domain.parms=[a b r alpha];

            case 'horizontal-circular-zone-boundary'
                % in this case "domain_struct.parms" is a vector [a b r alpha]
                a=1; b=1; r=4; alpha=-pi/4;
                basis_domain.parms=[a b r alpha];

            case 'circular-segment-boundary'
                % in this case "domain_parms" is a vector [a b r alpha beta]
                a=1; b=1; r=4; alpha=-pi/4; beta=pi/6;
                basis_domain.parms=[a b r alpha beta];

            case 'symmetric-lens-boundary'
                % in this case "domain_parms" is a vector [a r] (with a < r)
                a=1; r=3;
                basis_domain.parms=[a r];

            case 'butterfly-boundary'
                % in this case "domain_parms" is a vector [a b r alpha beta]
                a=1; b=2; r=3; alpha=-pi/4; beta=pi/6;
                basis_domain.parms=[a b r alpha beta];

            case 'candy-boundary'
                % in this case "domain_parms" is a vector [a r alpha] with
                % "-alpha > acos(a/r)"
                a=0.5; r=1; alpha=-1.5;
                basis_domain.parms=[a r alpha];

            case 'sector-boundary'
                center=[1 1.5]; r=0.5; alpha=pi/4; beta=pi/2;
                basis_domain.parms=[center r alpha beta];

            case 'circle-arc'
                x1=0; x2=0; r1=1; arc_angles=[-pi pi]; % center and radius
                basis_domain.parms=[x1 x2 r1 arc_angles];

            case 'disk'   % domain basis (disk)
                x1=0; x2=0; r1=1;  % center and radius
                basis_domain.parms=[x1 x2 r1];

            case {'polygon','polygonal-boundary'} % domain basis (polygon)

                basis_domain.parms=[0.1 0; 0.7 0.2; 1 0.5; 0.75 0.85; ...
                    0.5 1; 0 0.25; 0.1 0];

        end

        angle_intv=[0 pi/4]; % angle rotation
        dir_axis=[0 1 0]; % axis direction
        x0=[0 0 0];       % point in the axis

        domain_struct.domain='rotation-domain';
        domain_struct.parms={basis_domain,angle_intv,dir_axis,x0};








    case 'rotation-surface'


        % basis_domain.domain='polygonal-boundary';
        % basis_domain.domain='circle-arc';
        % basis_domain.domain='asymmetric-circular-sector-boundary';
        % basis_domain.domain='asymmetric-annulus-boundary';
        % basis_domain.domain='vertical-circular-zone-boundary';
        % basis_domain.domain='horizontal-circular-zone-boundary';
        % basis_domain.domain='circular-segment-boundary';
        % basis_domain.domain='symmetric-lens-boundary';
        % basis_domain.domain='butterfly-boundary';
        % basis_domain.domain='candy-boundary';
        basis_domain.domain='sector-boundary';

        switch basis_domain.domain

            case 'asymmetric-circular-sector-boundary'
                % in this case "domain_struct.parms" is a vector [a b c d r]
                a=1; b=1; c=2; d=2; r=3;
                basis_domain.parms=[a b c d r];

            case 'asymmetric-annulus-boundary'
                % in this case "domain_struct.parms" is a vector [a b c d r1 r2]
                a=1; b=1; c=2; d=2; r1=3; r2=4;
                basis_domain.parms=[a b c d r1 r2];

            case 'vertical-circular-zone-boundary'
                % in this case "domain_struct.parms" is a vector [a b r alpha]
                a=1; b=1; r=4; alpha=-pi/4;
                basis_domain.parms=[a b r alpha];

            case 'horizontal-circular-zone-boundary'
                % in this case "domain_struct.parms" is a vector [a b r alpha]
                a=1; b=1; r=4; alpha=-pi/4;
                basis_domain.parms=[a b r alpha];

            case 'circular-segment-boundary'
                % in this case "domain_parms" is a vector [a b r alpha beta]
                a=1; b=1; r=4; alpha=-pi/4; beta=pi/6;
                basis_domain.parms=[a b r alpha beta];

            case 'symmetric-lens-boundary'
                % in this case "domain_parms" is a vector [a r] (with a < r)
                a=1; r=3;
                basis_domain.parms=[a r];

            case 'butterfly-boundary'
                % in this case "domain_parms" is a vector [a b r alpha beta]
                a=1; b=2; r=3; alpha=-pi/4; beta=pi/6;
                basis_domain.parms=[a b r alpha beta];

            case 'candy-boundary'
                % in this case "domain_parms" is a vector [a r alpha] with
                % "-alpha > acos(a/r)"
                a=0.5; r=1; alpha=-1.5;
                basis_domain.parms=[a r alpha];

            case 'sector-boundary'
                center=[1 1.5]; r=0.5; alpha=pi/4; beta=pi/2;
                basis_domain.parms=[center r alpha beta];

            case 'circle-arc'
                x1=3; x2=0; r1=1; arc_angles=[-pi pi]; % center and radius
                basis_domain.parms=[x1 x2 r1 arc_angles];

            case {'polygonal-boundary'} % domain basis (polygon)
                basis_domain.parms=[0.1 0; 0.7 0.2; 1 0.5; 0.75 0.85; ...
                    0.5 1; 0 0.25; 0.1 0];

        end

        angle_intv=[-pi pi]; % angle rotation
        dir_axis=[0 1 0]; % axis direction
        x0=[0 0 0];       % point in the axis

        domain_struct.domain='rotation-surface';
        domain_struct.parms={basis_domain,angle_intv,dir_axis,x0};








    otherwise
        domain_struct.parms=[];

end

