
function [XYZW_C,moments,moments_C,err]=rule_compression(XYZW,ade,bbox,...
    compr_type)

%--------------------------------------------------------------------------
% Object:
%
% Compression of trivariate cubature formulas by Tchakaloff points.
% It is useful in node reduction of algebraic cubature formulas.
% For further informations, see the web page: 
%         http://www.math.unipd.it/~marcov/Tchakaloff.html
%
%--------------------------------------------------------------------------
% Input:
%
% XYZW: rule to compress of total degree "ade". It is a matrix of dimension 
%  m x 4. The first 3 columns are the coordinates of the nodes,
%  while the 4th column represent the corresponding weights.
%
% ade: polynomial (total) degree of exactness.
%
% bbox: to determine the required compression, we make use of shifted 
%  tensorial Chebyshev basis is used, where "bbox" determines the shift. 
%  In general "bbox" contains the cubature domain, and it
%  represents the smaller hyperrectangle with sides parallel to the axes
%  containing the domain. 
%  If no information is given, "bbox" represents the smaller hyperrectangle
%  with sides parallel to the axes containing the cubature nodes.
%
% compr_type: determines the algorithm to use for compression
%  0: M. Dessole, F. Marcuzzi, M. Vianello algorithm, 
%  1: Matlab built-in lsqnonneg, 
%  2: M. Slawski implementation of Lawson-Hanson algorithm, 
%  3: QR extraction (possible negative weights).
%
%--------------------------------------------------------------------------
% Output:
%
% XYZW_C: compressed of total degree "ade". It is a matrix of dimension 
%  m x 4. The first 3 columns are the coordinates of the nodes,
%  while the 4th column represent the corresponding weights. The number
%  of nodes is at most (ade+1)*(ade+2)*(ade+3)/6. Furtermore, weights are
%  positive if Lawson-Hanson algorithm is applied (i.e. "compr_type" is set
%  equal to "0", "1" or "2".
%
% moments: shifted tensorial Chebyshev (relatively to "box") computed by
%  the starting rule (defined by XYZW).
%
% moments_C: shifted tensorial Chebyshev (relatively to "box") computed by
%  the compressed rule (defined by XYZW_C).
%
% err: compression error in norm 2, i.e. "norm(moments-moments_C,2)".
%
%--------------------------------------------------------------------------
% Subroutines called by "rule_compression":
%
% 1. chebvand (and its subroutines).
% 2. LHDM
% 3. lsqnonneg
% 4. lawsonhanson
%
%--------------------------------------------------------------------------
% Authors:
% Alvise Sommariva and Marco Vianello,
% University of Padova, September 7, 2020.
%--------------------------------------------------------------------------


% .......................... troubleshooting ..............................

if nargin < 3
    X=XYZW(:,1); Y=XYZW(:,2); Z=XYZW(:,3);
    bbox=[min(X) max(X) min(Y) max(Y) min(Z) max(Z)];
end

if nargin < 4
    compr_type=0;
end

% ....................... compression startup .............................

% Computing moments of shifted tensorial Chebyshev basis (w.r.t. "bbox").
V=chebvand(XYZW(:,1:3),ade,bbox); 
moments=(XYZW(:,4))'*V; moments=moments';

[Q,R]=qr(V,0); % change of basis for stability ("V" may be ill-condit.).
moments1=R'\moments;

% ............................. compression ...............................


switch compr_type
    case 0
        % M. Dessole, F. Marcuzzi, M. Vianello algorithm
        [w,resnorm,exitflag, outeriter] = LHDM(Q', moments1);
        ind1=find(w>0);
    case 1
        % Matlab built-in "lsqnonneg"
        [w,RESNORM,RESIDUAL,EXITFLAG,OUTPUT]=lsqnonneg(Q',moments1);
        ind1=find(w>0);
    case 2
        % Lawson-Hanson algorithm by M.Slawski
        w = lawsonhanson(Q', moments1);
        ind1=find(w>0);
    otherwise
        % Extraction by QR method (possible negative weights!)
        w=Q'\moments1;
        ind1=find(w ~= 0);
end

XYZW_C=[XYZW(ind1,1:3) w(ind1)];

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

if nargout >= 3
    % Vandermonde matrix restricted to extracted nodes.
    VC=V(ind1,:); 
    % Computing moments of compressed rule
    moments_C=(XYZW_C(:,4))'*VC; moments_C=moments_C';
    % Making comparison with exact (non compressed) moments.
    err=norm(moments-moments_C,2);
end









%--------------------------------------------------------------------------
% chebvand
%--------------------------------------------------------------------------

function V = chebvand(points,ade,bbox)

%--------------------------------------------------------------------------
% Object:
% Vandermonde matrix at "points" points, relative to Chebyshev basis of 
% degree "ade" on rectangle "bbox". The basis is of tensorial type, 
% of total degree.
%--------------------------------------------------------------------------
% Input:
% ade: polynomial degree.
% points: 3-column array of point coordinates.
% bbox: 6-component vector such that the cube:
% [bbox(1),bbox(2)] x [bbox(3),bbox(4)] x [bbox(5),bbox(6)] contains X.
%--------------------------------------------------------------------------
% Output:
% V: Chebyshev-Vandermonde matrix at "points", graded lexic. order.
%    size is M x N, where M is the number of nodes and N is
%    (ade+1)*(ade+2)*(ade+3)/6.
%--------------------------------------------------------------------------
% Authors:
% Alvise Sommariva and Marco Vianello
% University of Padova, September 1, 2020
%--------------------------------------------------------------------------

if nargin < 3 
    bbox=[];
end

if isempty(bbox)
    bbox=[min(points(:,1)) max(points(:,1)) ...
        min(points(:,2)) max(points(:,2)) ...
        min(points(:,3)) max(points(:,3))];
end

degrees=[];
for ii=0:ade
    for jj=0:(ade-ii)
        for kk=0: ade-jj-ii % computing contribution
            degrees=[degrees; ii jj kk];
        end
    end
end

% mapping the mesh in the cube [-1,1]^3
map=[(2*points(:,1)-bbox(2)-bbox(1))/(bbox(2)-bbox(1)) ...
    (2*points(:,2)-bbox(4)-bbox(3))/(bbox(4)-bbox(3)) ...
    (2*points(:,3)-bbox(6)-bbox(5))/(bbox(6)-bbox(5))];

% Chebyshev-Vandermonde matrix on the mesh
T1=chebpolys(ade,map(:,1));
T2=chebpolys(ade,map(:,2));
T3=chebpolys(ade,map(:,3));
V=T1(:,degrees(:,1)+1).*T2(:,degrees(:,2)+1).*T3(:,degrees(:,3)+1);








%--------------------------------------------------------------------------
% chebpolys
%--------------------------------------------------------------------------

function T=chebpolys(deg,x)

%--------------------------------------------------------------------------
% Object:
% This routine computes the Chebyshev-Vandermonde matrix on the real line
% by recurrence.
%--------------------------------------------------------------------------
% Input:
% deg: maximum polynomial degree
% x: 1-column array of abscissas
%--------------------------------------------------------------------------
% Output:
% T: Chebyshev-Vandermonde matrix at x, T(i,j+1)=T_j(x_i), j=0,...,deg.
%    Its size is M x N where M is the number of nodes and N the dimension 
%    of the polynomial space (i.e. "deg+1").
%--------------------------------------------------------------------------
% Authors:
% Alvise Sommariva and Marco Vianello
% University of Padova, September 1, 2020.
%--------------------------------------------------------------------------

T=zeros(length(x),deg+1);
t0=ones(length(x),1);
T(:,1)=t0;
t1=x;
T(:,2)=t1;

for j=2:deg
    t2=2*x.*t1-t0;
    T(:,j+1)=t2;
    t0=t1;
    t1=t2;
end
