function leb = dLEB(deg,nodes,w,X,jvec,dbox)

%--------------------------------------------------------------------------
% Object:
% This routine computes the Lebesgue constant at "X" of d-variate  
% weighted least-squares polynomial fitting on "nodes".
%--------------------------------------------------------------------------
% Input:
% deg: polynomial degree;
% nodes: d-column array of points on which the Lebesgue constant is 
%    required, at degree "deg";
% w: 1-column array of positive weights or positive scalar in
%    case of equal weights;
% X: d-column array of control point coordinates, useful to estimate the
%    Lebesgue constant; 
% * jvec: vector of column indexes, selects a polynomial basis; 
% * dbox: variable that defines a hyperrectangle with sides parallel to the
%    axis, containing the domain. 
%    If "dbox" is not provided, it is the smaller "hyperrectangle", with 
%    sides parallel to the cartesian axes, containing the pointset "X".
%    It is a matrix with dimension "2 x d", where "d" is the dimension of
%    the space in which it is embedded the domain. 
%    For example, for a 2-sphere, it is "d=3", for a 2 dimensional polygon 
%    it is "d=2".
%    As example, the set "[-1,1] x [0,1]" is described as "[-1 0; 1 1]".
%
% Note: 
% The variables with an asterisk "*" are not mandatory and can be also set
% as empty matrix.
%--------------------------------------------------------------------------
% Output:
% leb: Lebesgue constant estimate based on evaluations at "X".
%--------------------------------------------------------------------------
% Dates:
% Written by M. Dessole, F. Marcuzzi, M. Vianello - July 2020
%
% Modified by M. Dessole, A. Sommariva, M. Vianello - November 2020
%--------------------------------------------------------------------------


% .........................  Function Body ................................

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

if nargin < 6, dbox=[]; end
if nargin < 5, jvec=[]; end
if nargin < 4, X=[]; end
if nargin < 3, w=[]; end
if isempty(X), X=nodes; end
if isempty(w), w=1; end





% ..... main code below .....

if isempty(dbox)
Z=[X; nodes];
a=min(Z); b=max(Z); dbox=[a;b];
end

[UX,~] = dCHEBVAND(deg,X,dbox);
[VY,jvec,~,R,~] = dORTHVAND(deg,nodes,w,jvec,[],dbox);

VX = UX(:,jvec)/R;   
VYVXt = VY*VX';    
if isscalar(w)
    V = w*VYVXt;
else
    V = zeros(size(VYVXt));
    for k=1:size(V,2)
        V(:,k)=VYVXt(:,k).*w;
    end
end

leb = norm(V,1);

