function [c,p,rss,r2] = clsap(xj, f, n, varargin)
%%  CLSAP   Classical least squares approximation problem
%   
%   [C] = CLSAP(XJ, F, N) computes the coefficients from the polynomial of
%   least squares of degree N in Gram basis given the equidistant data
%   F(XJ) in the interval [-1+1/NUP, 1-1/NUP], where NUP is a large number,
%   NUP>>N.
% 
%   The polynomial is:
%      C(1)*G_N(x) + C(2)*G_{N-1}(x) +...+ C(N)*G_1(x) + C(N+1)*G_0(x),
%      where G_K is the orthonormal Gram polynomial of degree K and
%      parameter NUP.
%   
%   CLSAP(XJ, F, N, EXPLICIT_F) allows the user furnishes the string
%   EXPLICIT_F with the function F of a variable x. It is highly
%   recommended when F is explicitly known.
%
%   [C] = CLSAP(XJ, F, N, 'MON') also computes the coefficients from the
%   polynomial in monomial basis in the matrix C, where the first row
%   contains the coefficients in Gram basis, and the second one contains
%   the coefficients in monomial basis, both in descending order of degree. 
%
%   CLSAP(XJ, F, N, EXPLICIT_F, 'MON') or CLSAP(XJ, F, N, 'MON', EXPLICIT_F)
%   produces:
%           CLSAP(XJ, F, N, EXPLICIT_F) and CLSAP(XJ, F, N, 'MON').
%
%   [C,P] = CLSAP(XJ, F, N) returns also a vector P with the polynomial
%   evaluated at the points XJ.
%
%   [C,P,RSS] = CLSAP(XJ, F, N) returns also the residual sum of squares
%   RSS.
%
%   [C,P,RSS,R2] = CLSAP(XJ, F, N) returns also the coefficient of
%   determination R2.
%
%       
%   Calls GRAM, REC.
%
% 
%   Copyright 2020 by Dimitar K. Dimitrov and Lourenco L. Peixoto
%   All rights reserved.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   CLSAP implemented by Lourenco L. Peixoto, 2020 - see [1-3].
%
%   Basic references:
%   
%       [1] I. Area, D. K. Dimitrov, E. Godoy, and V. Paschoa, "Approximate
%   calculation of sums I: Bounds for the zeros of Gram polynomials", SIAM
%   J. Numer. Anal., 2014.
%
%       [2] I. Area, D. K. Dimitrov, E. Godoy, and V. Paschoa, "Approximate
%   calculation of sums II: Gaussian type quadrature", SIAM J. Numer.
%   Anal., 2016.
%   
%       [3] D. K. Dimitrov and L. L. Peixoto, "An efficient algorithm for
%   the classical least squares approximation", submitted, 2020.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% standard parameters:
base = 'gram';
explicit_f = [];

if (nargin > 3)
    if (nargin > 4)
        if ~strcmpi(varargin{1},'mon')
            % clsap(xj, f, n, 'explicit_f', 'mon')
            base = 'mon';
            explicit_f = varargin{1};
        else
            % clsap(xj, f, n, 'mon', 'explicit_f')
            base = 'mon';
            explicit_f = varargin{2};
        end
    else
        if ~strcmpi(varargin{1},'mon')
            % clsap(xj, f, n, 'explicit_f')
            explicit_f = varargin{1};
        else
            % clsap(xj, f, n, 'mon')
            base = 'mon';
        end
    end
end

% number of data:
nup = length(xj);

% check the mesh-points:
if ( xj(1)~=-1+1/nup || xj(end)~=1-1/nup )
    error('The mesh-points must be in interval [-1+1/NUP, 1-1/NUP], NUP>>N.')
end

m0 = floor(2.5*sqrt(nup));

% number of points for the quadrature:
m = min(100,m0);

% initialise:
rm = inf;
m_max = min(700,m0);
rm_max = 1e-5;

% increase m if necessary: 
while (rm > rm_max && m <= m_max )

    % check the accuracy:
    % computes F in the zeros of the Gram polynomial of degree m-5:
    [g,w] = gram(m-5,nup);
    [fg] = eval_f(f, g, nup, explicit_f); 

    % sum(f^2)/nup with m-5 points:
    q0 = w*transpose(fg.^2);

    % computes F in the zeros of the Gram polynomial of degree m:
    [g,w] = gram(m,nup);
    [fg] = eval_f(f, g, nup, explicit_f); 

    % sum(f^2)/nup with m points:
    q1 = w*transpose(fg.^2);

    rm = abs(q1-q0)/q1;
    
    % increase m:
    m = m + 50;
end

epsilon = (5*(rm + 2*eps))^2;

        
% evaluates the Gram polynomials of degree 0,...,N with parameter NUP at the
% zeros of the Gram polynomial of degree m:
[~,MG] = rec(n,nup,transpose(g)); % MG is a (N+1)-by-M matrix

% the vector with the inner products <F,G_k>, k=0,...,n, ie, the vector
% with the Fourier coefficients:
a = w.*fg*transpose(MG);

% find the null Fourier coefficients:
if (rm <= rm_max)
    lima = a.^2/q1;
    v = find(lima<=epsilon);
    a(v) = zeros(1,numel(v));
end

% coefficients in descending order of degree:
c = flip(a); % c(1)*G_n(x) + c(2)*G_{n-1}(x) +...+ c(n)*G_1(x) + c(n+1)*G_0(x)

% monomial basis:
if strcmp(base,'mon')
    cm = monomial(n,nup,a); % monomial basis
    c = [c; cm]; % first row contains Gram basis, the second one contains the monomial ones. Both in descending order of degree.
end
    

% Optionally evaluate the polynomial and compute RSS and R2:
if nargout>1 % polynomial at mesh-points:
    v = find(a~=0,1,'last'); % use the precise degree
    p = clenshaw(v-1,nup,a(1:v),xj); % Clenshaw's algorithm
    if nargout>2 % the RSS:
        rss = sum((f-p).^2);
        if nargout>3 % the coefficient of determination R2:
            r2 = 1 - rss / ( sum(f.^2) - (sum(f))^2/nup );
        end
    end
end


end

function [c] = monomial(n,nup,a)
% Computes the coefficients in monomial basis.

% the matrix with the coefficients from all polynomials of degree 0,...,N:
B = coefmatrix_full(n,nup);

% the coefficients in monomial basis:
c = sum(transpose(a).*B,1);

c = flip(c);

end

function [fg] = eval_f(f, g, nup, explicit_f)
%     EVAL_F evaluate F in the points G. It uses linear interpolant when
%     the function is unknown. 

% evaluate f(g):
if ischar(explicit_f)
    % eval_f(f, g, nup, explicit_f)
    x=g;
    fg = transpose(eval(explicit_f));
else
    % eval_f(f, g, nup, []) 
    fg = linear_interp(f,g,nup);
end
    

end
        

function [fg] = linear_interp(f,g,nup)
% approximate the data F(G) using linear interpolant between the mesh
% points, where x_t < G < x_{t+1}, t=0,...,NUP-2.

% identify every zero G between the mesh points:
t = transpose(floor(.5*nup*(g+1-1/nup)));

x0 = -1+(1+2*t)/nup; x1 = x0+2/nup; % abscissas
y0 = f(t+1); y1 = f(t+2);           % F(x0), F(x1)
a1 = (y1-y0)./(x1-x0);              % linear coefficient
a0 = y0-a1.*x0;                     % independent coefficient
fg = a0+a1.*transpose(g);           % linear interpolant
end




function B = coefmatrix_full(max,nup)
% Compute all coefficients from Gram polynomials of degree = 0,...,MAX

% B is the square matrix of order MAX+1 with the coefficients of the
% orthonormal Gram polynomials G_0, G_1, ..., G_MAX with parameter NUP.
%
%       | b_{0,0}     0         0       0   ...    0|
%       | b_{1,0}  b_{1,1}      0       0   ...    0|
%       | b_{2,0}  b_{2,1}   b_{2,2}    0   ...    0|
%   B = |    .                                     .| 
%       |    .                                     .| 
%       |    .                                     .| 
%       | b_{n,0}  b_{n,1}   b_{n,2}   ...   b_{n,n}|
% 

B = zeros(max+1);

B(1,1) = 1;

if (max>0)
    k=1:max;
    alp = sqrt(nup^2*(k.^2-.25)./(k.^2.*(nup^2-k.^2))); % alp(k) = alpha_{k-1}, 
    % where alpha is a term in the REC of the orthonormal polynomials.
    
    B(2,2) = 2*alp(1);
    
    if (max>1)
        B(3,3) = 2*alp(2)*B(2,2); B(3,1) = -alp(2)/alp(1);
        for i=4:max+1
            im1 = i-1;
            im2 = im1-1;
            B(i,i) = 2*alp(im1)*B(im1,im1);
            for j=i-2:-2:2
                B(i,j) = 2*alp(im1)*B(im1,j-1) - alp(im1)/alp(im2)*B(im2,j);
            end
            B(i,1) = -alp(im1)/alp(im2)*B(im2,1);
        end
    end
end

end



function s = clenshaw(n,nup,d,x)
% Clenshaw's algorithm to the Gram orthonormal polynomial.

% S = CLENSHAW(N,NUP,D,X) returns the summation:
%
%   D(1)*G_0(X) + D(2)*G_1(X) + ... + D(n+1)*G_n(X),
% 
% where G obeys the relation: G_k - a_k*G_{k-1} - b_k*G_{k-2} = 0, b_k~=0,
% k=2,...,N.
% 
% G_k is the orthonormal Gram polynomial of degree k and parameter NUP. 
%
%
% See (Deuflhard, 1976) and (Clenshaw, 1955).

% computes alpha_0,...,alpha_{n-1}:
k=1:n;
alp = sqrt(nup^2*(k.^2-.25)./(k.^2.*(nup^2-k.^2))); 

nt = length(x);
G0 = ones(1,nt);

if ( n > 1 )
    % initialise:
    u(n,1:nt) = d(n+1);
    u(n-1,:) =  2*alp(n)*x.*u(n,:) + d(n); 

    for k=n-1:-1:2
        u(k-1,:) = 2*alp(k)*x.*u(k,:) - alp(k+1)/alp(k) *u(k+1,:) + d(k);  
    end
    
    G1 = 2*alp(1)*x;
    s = (d(1) - alp(2)/alp(1) * u(2,:)).*G0 + u(1,:).*G1;
    
elseif ( n > 0 ) % n=1
   
    G1 = 2*alp(1)*x;
    s = d(1)*G0 + d(2).*G1;
    
else % n=0
    s = d(1)*G0;
    
end

end