
function demo_cubature_cputime

% OBJECT:
% the procedure tests the cubature on curvilinear polygons

n=6; % degree of precision
tests=1; % number of tests
extraction_method=1; % 1: lsqnonneg, 2: fast LH, 0: backslash
domain_type=1; % non convex domain "domain_type=0", 
               % convex domain "domain_type=1"
               
               


for kk=1:tests
    [cpuL(kk),xyw,Z,Zin,res]=demosplinecub3(n,domain_type,extraction_method);
end

fprintf('\n \n \t ------------ STATS ------------');
fprintf('\n \n \t * degree: %4.0f',n);
fprintf('\n \t * cardinality polynomial space: %8.0f',(n+1)*(n+2)/2);
fprintf('\n \t * cardinality extracted nodes : %8.0f',size(xyw,1));
fprintf('\n \t * points analysed     : %8.0f',size(Z,1));
fprintf('\n \t * points inside domain: %8.0f \n',size(Zin,1));
fprintf('\n \t * cubature conditioning : %1.3e',...
    sum(abs(xyw(:,3)))/sum((xyw(:,3))));
fprintf('\n \t * residuals             : %1.3e \n',res);
fprintf('\n \t tests    : %3.0f',tests);
fprintf('\n \t aver cpu : %1.3e',mean(cpuL));
fprintf('\n \n')




%==========================================================================
% Attached functions
%==========================================================================


%==========================================================================
% demosplinecub3
%==========================================================================

function [cpuL,xyw,Z,Zin,res]=demosplinecub3(n,domain_type,extraction_method)

% ................... troubleshooting / settings ..........................

if nargin < 1
    n=25;
end

if nargin < 2
    domain_type=0;
end

if nargin < 3
    extraction_method=2;
end

do_stats=0; % stats: set "do_stats=1"
do_plot=1;  % plot: set "do_plot=1"
example=1;

% ..................... determining problem ...............................

switch example
    case 1
        vert=[-1 0; -2 -1; -1.5 -2; 0 -1.6; 0 -1];
end

% .................. computing spline boundary ............................

% ... straight line component ...

t=linspace(0,1,size(vert,1)-1);
L=length(vert(:,1))-1;

for i=1:L
    for j=1:length(t)
        x(j)=vert(i,1)+t(j)*(vert(i+1,1)-vert(i,1));
        y(j)=vert(i,2)+t(j)*(vert(i+1,2)-vert(i,2));
    end
    Sx(i)=spline(t,x);
    Sy(i)=spline(t,y);
end


% ... arc component ...
a=0; b=pi/2; m=101; theta=linspace(a,b,m);
switch domain_type
    case 0
        
        x=cos(3*pi/2-theta);
        y=sin(3*pi/2-theta);
        
    case 1
        x=-1+cos(theta);
        y=-1+sin(theta);
        
end

xV=[vert(:,1); x']; yV=[vert(:,2); y'];
for kk=1:length(xV)
fprintf('\n \t %1.15e %1.15e',xV(kk),yV(kk));
end
pause



Sx(L+1)=spline(theta,x);
Sy(L+1)=spline(theta,y);

for ii=1:L
    SxL=Sx(i);
    SxL.breaks;
end

% .................. determining cubature rule ............................

tic;
[xyw,res,Z,Zin]=splcub(n,Sx,Sy,extraction_method);
cpuL=toc;

% format long e
% xyw(:,1:2)

% .............................. stats ....................................

if do_stats == 1
    fprintf('\n \n \t * degree: %4.0f',n);
    fprintf('\n \t * cardinality nodes: %8.0f',size(xyw,1));
    fprintf('\n \t * points analysed: %8.0f inside: %8.0f',...
        size(Z,1),size(Zin,1));
    fprintf('\n \t * ratios: %1.3e',sum(abs(xyw(:,3)))/sum((xyw(:,3))));
end

% ............................ plotting sets ..............................

if do_plot == 1
    
    % ... plot boundary ...
    clf; hold on;
    
    % ... plot all points ...
    plot(Z(:,1),Z(:,2),'ro');
    
    % ... plot points inside the domain ...
    plot(Zin(:,1),Zin(:,2),'go');
    
    % ... plot cubature nodes ...
    plot(xyw(:,1),xyw(:,2),'g*','MarkerEdgeColor','k',...
        'MarkerFaceColor','g',...
        'MarkerSize',4)
    
    % ... plot boundary ...
    for k=1:length(Sx)
        SxL=Sx(k);
        SyL=Sy(k);
        t0=SxL.breaks(1);
        t1=SxL.breaks(end);
        tt=linspace(t0,t1,1000);
        xx=ppval(SxL,tt);
        yy=ppval(SyL,tt);
        plot(xx,yy,'k-')
    end
    
    axis equal;
    hold off;
end




