function [fek,P,leb] = approxfek0(deg,mymesh,leja,controlmesh,rect);

% extracts approximate Fekete or Leja interpolation points from a 2d 
% or 3d mesh/cloud and estimates their Lebesgue constant

% by Alvise Sommariva and Marco Vianello, University of Padova 
% November 22, 2012

% see, e.g., the papers: 
%
% A. Sommariva and M. Vianello, Computing approximate Fekete points by QR 
% factorizations of Vandermonde matrices, Comput. Math. Appl. 57 (2009)
% preprint: http://www.math.unipd.it/~marcov/pdf/afek.pdf
% 
% L. Bos, S. De Marchi, A. Sommariva and M. Vianello, Computing 
% multivariate Fekete and Leja points by numerical linear algebra, 
% SIAM J. Numer. Anal. 48 (2010)
% preprint: http://www.math.unipd.it/~marcov/pdf/leja.pdf


% INPUT: 

% deg = interpolation degree

% mymesh = 2- or 3-column array of mesh points coordinates, from which the 
% approximate Fekete or Leja points will be extracted; if mymesh is 
% already a set of interpolation points, their Lebesgue constant will be 
% estimated on controlmesh 

% leja=0: computes approximate Fekete points
% leja=1: computes approximate Leja points

% controlmesh = 2- or 3-column array of mesh point coordinates
% (possibly finer mesh to evaluate the Lebesgue constant of the points) 
% if controlmesh == [ ], then it will be set controlmesh = mymesh

% rect = 4- or 6-component vector such that the rectangle 
% [rect(1),rect(2)] x [rect(3),rect(4)] in 2d
% or [rect(1),rect(2)] x [rect(3),rect(4)] x [rect(5),rect(6)] in 3d
% contains the meshes
% if rect == [ ], it will be automatically tailored to the meshes 


% OUTPUT: 

% fek = 2- or 3-column array of approximate Fekete or Leja 
% points coordinates

% P = transformation matrix for the total-degree product Chebyshev basis 
% of the rectangle, namely 
% [new basis] = [old basis] * P
% [new Vandermonde] = [old Vandermonde] * P

% leb = approximate Lebesgue constant
% a call without this output argument will directly avoid to compute it 


% FUNCTION BODY

% setting rectangle if not defined 
if size(rect) == [0 0]
allpts=[mymesh;controlmesh];
rect(1)=min(allpts(:,1));
rect(2)=max(allpts(:,1));
rect(3)=min(allpts(:,2));
rect(4)=max(allpts(:,2));
if length(mymesh(1,:)) == 3
rect(5)=min(allpts(:,3));
rect(6)=max(allpts(:,3));
end;
end;

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

% approximate Fekete or Leja points extraction from mymesh 

% computing the rectangular Chebyshev-Vandermonde matrix on the mesh 
W=chebvand(deg,mymesh,rect);

% change of basis via 2 successive QR factorizations of the 
% Vandermonde matrices; to keep the original basis set ort=0
ort=1; 
B=W;
if ort==1
[Q,R]=qr(B,0);
P1=inv(R);
B=B*P1;
[Q,R]=qr(B,0);
P2=inv(R);
B=B*P2;
P=P1*P2;
else
P=eye(dim);
end;

% P is the basis transformation matrix: 
% now B = W * P is the Chebyshev-Vandermonde matrix on mymesh 
% in the new basis and is nearly-orthogonal 
%(unless W is too severely ill conditioned, say cond(W)>1/eps) 

% extracting the interpolation points and computing the weights 

if leja == 0
% approximate Fekete points and corresponding weights
% by QR factorization with column pivoting of B'=(W*P)'
w=B'\ones(dim,1); 
ind=find(abs(w)>0);
fek=mymesh(ind,:);
else
% approximate Leja points and corresponding weights  
% by LU factorization with row pivoting of B=W*P
[L,U,perm]=lu(B,'vector');
fek=mymesh(perm(1:dim),:);
end;

% estimating the Lebesgue constant

% square Chebyshev-Vandermonde matrix at the interpolation points
V=chebvand(deg,fek,rect);

% setting controlmesh if not defined and computing the rectangular 
% Chebyshev-Vandermonde matrix on controlmesh
if size(controlmesh) == [0 0]
controlmesh=mymesh;
controlV=W;
else
controlV=chebvand(deg,controlmesh,rect);
end;

if nargout == 3
% numerical maximization of the Lebesgue function on controlmesh, 
% stabilizing by the basis transformation matrix 
leb=norm((V*P)'\(controlV*P)',1);
end;

end


function V = chebvand(deg,gmesh,rect);

% computes the Chebyshev-Vandermonde matrix on a 2d or 3d given mesh 
% in the total-degree product Chebyshev basis of a rectangle containing 
% the mesh points

% INPUT:

% deg = polynomial degree

% gmesh = 2- or 3-column array of mesh point coordinates

% rect = 4- or 6-component vector such that the rectangle 
% [rect(1),rect(2)] x [rect(3),rect(4)] in 2d
% or [rect(1),rect(2)] x [rect(3),rect(4)] x [rect(5),rect(6)] in 3d
% contains the mesh


% OUTPUT:

% V = Chebyshev-Vandermonde matrix 


% FUNCTION BODY

if length(gmesh(1,:)) == 2 

% couples with length less or equal to deg
j=(0:1:deg);
[j1,j2]=meshgrid(j);
good=find(j1(:)+j2(:)<deg+1);
couples=[j1(good) j2(good)];

% mapping the mesh in the square [-1,1]^2
a=rect(1);b=rect(2);c=rect(3);d=rect(4);
map=[(2*gmesh(:,1)-b-a)/(b-a) (2*gmesh(:,2)-d-c)/(d-c)];

% Chebyshev-Vandermonde matrix on the mesh
V1=cos(couples(:,1)*acos(map(:,1)'));
V2=cos(couples(:,2)*acos(map(:,2)'));
V=(V1.*V2)';

end;


if length(gmesh(1,:)) == 3 
    
%triples with length less or equal to deg
j=(0:1:deg);
[j1,j2,j3]=meshgrid(j,j,j);
good=find(j1(:)+j2(:)+j3(:)<deg+1);
triples=[j1(good) j2(good) j3(good)];

%mapping the mesh in the cube [-1,1]^3
a=rect(1);b=rect(2);
c=rect(3);d=rect(4);
e=rect(5);f=rect(6);
map=[(2*gmesh(:,1)-b-a)/(b-a) (2*gmesh(:,2)-d-c)/(d-c) ... 
(2*gmesh(:,3)-f-e)/(f-e)];

%Chebyshev-Vandermonde matrix on the mesh
V1=cos(triples(:,1)*acos(map(:,1)'));
V2=cos(triples(:,2)*acos(map(:,2)'));
V3=cos(triples(:,3)*acos(map(:,3)'));
V=((V1.*V2).*V3)';

end; 

end

