function [nwam,nweights] = wamcomprex(deg,wam,kind)

% extracts a low cardinality WAM (Weakly Admissible Mesh) from a given one
% works in 2d and 3d
% works also on algebraic varieties with the appropriate polynomial space 
% see the web page: http://www.math.unipd.it/~marcov/Tchakaloff.html

% by Federico Piazzon, Alvise Sommariva and Marco Vianello
% version 2.0, University of Padova, May 16 - 2016

% DISCLAIMER: These programs are free software; you can redistribuite them
% and/or modify them under the terms of the GNU/General Public License as
% published by the Free Software Foundation; either version 2 of the 
% License, or (at your opinion) any later version

% input

% deg: polynomial degree
% wam: 2-column or 3-column array of WAM points  
% kind: selection mode 
% 1: Tchakaloff points, 2: approximate Fekete points 
 
% output
% nwam: 2-column or 3-column array of the extracted WAM points 
% nweights: 1-column array of weights 

% warning: could be slow at moderate degrees, especially in 3d 


% function body

% unit weights at the wam points
w=ones(length(wam(:,1)),1);

% dimension of the polynomial space of degree 2*deg
% number of variables
dim=length(wam(1,:));
if dim == 2
N=(2*deg+1)*(2*deg+2)/2;
end
if dim == 3
N=(2*deg+1)*(2*deg+2)*(2*deg+3)/6;
end;

if ~(length(wam(:,1))>N)
% no compression possible
nwam=wam;
nweights=w;
else
% Chebyshev-Vandermonde matrix at wam 
a=min(wam(:,1)); b=max(wam(:,1)); c=min(wam(:,2)); d=max(wam(:,2));
if dim == 2
rect=[a b c d];
end
if dim == 3
e=min(wam(:,3)); f=max(wam(:,3));
rect=[a b c d e f];
end

V=chebvand(2*deg,wam,rect);
% polynomial basis orthogonalization, column pivoting
[Q,R,E]=qr(V,0);
r=rank(R);
% restriction to a polyn. basis on wam (polyn. space dimension = r)
Q=real(Q(:,1:r));
% possible re-orthogonalization
% [Q,R]=qr(Q,0);

% discrete moments of the orthogonal basis
mom=Q'*w;

% extracts the new nodes and weights  
if kind == 1
% NonNegative Least Squares: Tchakaloff points 
ww=lsqnonneg(Q',mom);
else 
% QR with column pivoting: approximate Fekete points 
ww=Q'\mom;
end;
% indexes of nonvanishing weights
ind=find(abs(ww)>0);
% compression 
nwam=wam(ind,:);
nweights=ww(ind);

end


