function [Xrls, yrls, s]=RLTSreg(y, X, max_fits, max_points)
% % RLTSreg: Calculates the Reweighted Least Trimmed Squares (RTLS) regression 
% % 
% % Syntax;
% % 
% % [Xrls,yrls,s]=RLTSreg(y, X, max_fits, max_points);
% %
% % **********************************************************************
% % 
% % Description
% % 
% % Calculates the Reweighted Least Trimmed Squares (RTLS) regression 
% % data points and the scale parameter from the LMTS regression.
% % 
% % This program is a modification of RLSreg.  It has been modified to 
% % trim the input data sets and trim the number of combinations of 
% % line fits that are processed.  The trimming allows the program to 
% % accomodate large data sets.
% % 
% % Ths program performs the Least Median Trimmed Squares Robust
% % Regression for simple or multiple columns of data and outputs the
% % regression parameters. 
% % 
% % **********************************************************************
% % 
% % Input Variable Description
% % 
% % y is the vector of the dependent variable.
% % 
% % X is the data matrix of the independent variable.
% %
% % max_fits is the number of best fit pairs of data. 
% %     The maximum value is 10000.
% %     The default value is 1000 or the largest value allowed.  
% %
% % max_points is the number of data points for curve fitting.
% %     The maximum value is 100000.  
% %     The default value is 100000 or the largest value allowed.
% % 
% % **********************************************************************
% % 
% % Output Variable Description
% % 
% % Xrls is the X values matrix to be taken into account for RTLS.
% % 
% % yrls is the y values vector to be taken into account for RTLS.
% % 
% % s is RLS scale parameter.
% % 
% % 
% % **********************************************************************
% % 
% % Reference:
% % Rousseeuw PJ, Leroy AM (1987): Robust regression and outlier detection. Wiley.
% % 
% % **********************************************************************
% % 
% % This program is originally the work of
% % 
% % Alexandros Leontitsis
% % Institute of Mathematics and Statistics
% % University of Kent at Canterbury
% % Canterbury
% % Kent, CT2 7NF
% % U.K.
% % 
% % University e-mail: al10@ukc.ac.uk (until December 2002)
% % Lifetime e-mail: leoaleq@yahoo.com
% % Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
% % 
% % Sep 3, 2001.
% % 
% % **********************************************************************
% %
% % This program was modified by Edward L. Zechmann
% %
% %     date  1 February    2008    Updated comments.
% %                                 Added trimmed programs.
% % 
% % modified  14 February   2008    Updated comments.
% %                                 modified error handling.
% % 
% % modified  2 December    2008    Fixed a bug in trimming the input data 
% %                                 arrays.  This fix improves accuracy 
% %                                 for data sets with less than 1000 
% %                                 points.  
% %
% %
% % 
% % **********************************************************************
% %
% % Feel free to modify this code.
% % 
% % See also: RLTSregor, RMTVE, MTVE, LMSpolor, LMTSpol, LMSpol, LMTSreg, LMSreg, LMTSregor, LMSregor
% % 

flag=0;

if nargin < 1 || isempty(y) || ~isnumeric(y)
   warning('Not enough input arguments.  Return empty array.');
   y=3;
   flag=1;
end

if flag==1
    Xrls=[];
    yrls=[];
    s=[];
else

    % Estimate the LMS values
    if nargin < 2 || isempty(X) || ~isnumeric(X)
        X=(1:length(y))';
    end

    pp=size(X,2);

    % If not input, set the maximum number of fits
    if nargin < 3  || isempty(max_fits) || ~isnumeric(max_fits)
        % default value of max_fits is 1000
        max_fits=min([1000, nchoosek(n, pp+1)]);
    end

    % make sure that max_fits does not exceed 10000
    max_fits=min( [max_fits, nchoosek(n, pp+1), 10000]);

    % If max_points is not an input, set the maximum number of points
    % for the input arrays X and y to a reasonable value.
    if nargin < 4 || isempty(max_points) || ~isnumeric(max_points)
        max_points=max([min([n, 100000]), max_fits*(pp+1)]);
    end

    if max_points < max_fits
        max_points=max_fits;
    end
    
    % Compute the robust regression 
    LMSout=LMTSreg(y, X, max_fits, max_points);

    % p is the number of parameters to be estimated
    p=size(X,2)+1;

    % Calculate the residuals
    r=y-LMSout;

    % Estimate the preliminary scale parameter
    s=LMSsca(r,0,p);

    % Take into account a data point, if its residual is relatively small
    w=find(abs(r/s)<=2.5);
    Xrls=X(w,:);
    yrls=y(w);

end