function [rew,raw]=mcdcov(x,varargin)

%MCDCOV computes the MCD estimator of a multivariate data set.  This
% estimator is given by the subset of h observations with smallest covariance
% determinant.  The MCD location estimate is then the mean of those h points,
% and the MCD scatter estimate is their covariance matrix.  The default value
% of h is roughly 0.75n (where n is the total number of observations), but the
% user may choose each value between n/2 and n. Based on the raw estimates,
% weights are assigned to the observations such that outliers get zero weight.
% The reweighted MCD estimator is then given by the mean and covariance matrix
% of the cases with non-zero weight. To compute the MCD estimator,
% the FASTMCD algorithm is used.
%
% The MCD method is intended for continuous variables, and assumes that
% the number of observations n is at least 5 times the number of variables p.
% If p is too large relative to n, it would be better to first reduce
% p by variable selection or robust principal components (see the functions
% robpca.m and rapca.m).
%
% The MCD method was introduced in:
%
%   Rousseeuw, P.J. (1984), "Least Median of Squares Regression,"
%   Journal of the American Statistical Association, Vol. 79, pp. 871-881.
%
% The MCD is a robust method in the sense that the estimates are not unduly
% influenced by outliers in the data, even if there are many outliers.
% Due to the MCD's robustness, we can detect outliers by their large
% robust distances. The latter are defined like the usual Mahalanobis
% distance, but based on the MCD location estimate and scatter matrix
% (instead of the nonrobust sample mean and covariance matrix).
%
% The FASTMCD algorithm uses several time-saving techniques which
% make it available as a routine tool to analyze data sets with large n,
% and to detect deviating substructures in them. A full description of the
% algorithm can be found in:
%
%   Rousseeuw, P.J. and Van Driessen, K. (1999), "A Fast Algorithm for the
%   Minimum Covariance Determinant Estimator," Technometrics, 41, pp. 212-223.
%
% An important feature of the FASTMCD algorithm is that it allows for exact
% fit situations, i.e. when more than h observations lie on a (hyper)plane.
% Then the program still yields the MCD location and scatter matrix, the latter
% being singular (as it should be), as well as the equation of the hyperplane.
%
%
% Required input argument:
%    x : a vector or matrix whose columns represent variables, and rows represent observations.
%        Missing values (NaN's) and infinite values (Inf's) are allowed, since observations (rows)
%        with missing or infinite values will automatically be excluded from the computations.
%
% Optional input arguments:
%       cor : If non-zero, the robust correlation matrix will be
%             returned. The default value is 0.
%         h : The quantile of observations whose covariance determinant will
%             be minimized.  Any value between n/2 and n may be specified.
%             The default value is 0.75*n.
%     alpha : (1-alpha) measures the fraction of outliers the algorithm should
%             resist. Any value between 0.5 and 1 may be specified. (default = 0.75)
%    ntrial : The number of random trial subsamples that are drawn for
%             large datasets. The default is 500.
%     plots : If equal to one, a menu is shown which allows to draw several plots,
%             such as a distance-distance plot. (default)
%             If 'plots' is equal to zero, all plots are suppressed.
%             See also makeplot.m
%   classic : If equal to one, the classical mean and covariance matrix are computed as well.
%             (default = 0)
%
% Input arguments for advanced users:
%     Hsets : Instead of random trial h-subsets (default, Hsets = []), Hsets makes it possible to give certain
%             h-subsets as input. Hsets is a matrix that contains the indices of the observations of one
%             h-subset as a row.
%    factor : If not equal to 0 (default), the consistency factor is adapted. Only usefull in case of the
%             kmax approach.
%
% I/O: result=mcdcov(x,'alpha',0.75,'h',h,'ntrial',500)
%  If only one output argument is listed, only the final result ('result')
%  is returned.
%  The user should only give the input arguments that have to change their default value.
%  The name of the input arguments needs to be followed by their value.
%  The order of the input arguments is of no importance.
%
% Examples: [rew,raw]=mcdcov(x);
%           result=mcdcov(x,'h',20,'plots',0);
%           [rew,raw]=mcdcov(x,'alpha',0.8,'cor',0)
%
% The output structure 'raw' contains intermediate results, with the following
% fields :
%
%     raw.center : The raw MCD location of the data.
%        raw.cov : The raw MCD covariance matrix (multiplied by a consistency factor).
%        raw.cor : The raw MCD correlation matrix, if input argument 'cor' was non-zero.
%         raw.wt : Weights based on the estimated raw covariance matrix 'raw.cov' and
%                  the estimated raw location 'raw.center' of the data. These weights determine
%                  which observations are used to compute the final MCD estimates.
%  raw.objective : The determinant of the raw MCD covariance matrix.
%
% The output structure 'rew' contains the final results, namely:
%
%       rew.center : The robust location of the data, obtained after reweighting, if
%                    the raw MCD is not singular.  Otherwise the raw MCD center is
%                    given here.
%          rew.cov : The robust covariance matrix, obtained after reweighting, if the raw MCD
%                    is not singular.  Otherwise the raw MCD covariance matrix is given here.
%          rew.cor : The robust correlation matrix, obtained after reweighting, if
%                    options.cor was non-zero.
%            rew.h : The number of observations that have determined the MCD estimator,
%                    i.e. the value of h.
%    rew. Hsubsets : A structure that contains Hopt and Hfreq:
%                    Hopt  : The subset of h points whose covariance matrix has minimal determinant,
%                            ordered following increasing robust distances.
%                    Hfreq : The subset of h points which are the most frequently selected during the whole
%                            algorithm.
%        rew.alpha : (1-alpha) measures the fraction of outliers the algorithm should
%                    resist.
%           rew.md : The distance of each observation from the classical
%                    center of the data, relative to the classical shape
%                    of the data. Often, outlying points fail to have a
%                    large Mahalanobis distance because of the masking
%                    effect.
%           rew.rd : The distance of each observation to the final,
%                    reweighted MCD center of the data, relative to the
%                    reweighted MCD scatter of the data.  These distances allow
%                    us to easily identify the outliers. If the reweighted MCD
%                    is singular, raw.rd is given here.
%       rew.cutoff : Cutoff values for the robust and mahalanobis distances
%         rew.flag : Flags based on the reweighted covariance matrix and the
%                    reweighted location of the data.  These flags determine which
%                    observations can be considered as outliers. If the reweighted
%                    MCD is singular, raw.wt is given here.
%       rew.method : A character string containing information about the method and
%                    about singular subsamples (if any).
%        rew.plane : In case of an exact fit, rew.plane contains the coefficients
%                    of a (hyper)plane a_1(x_i1-m_1)+...+a_p(x_ip-m_p)=0
%                    containing at least h observations, where (m_1,...,m_p)
%                    is the MCD location of these observations.
%      rew.classic : If the input argument 'classic' is equal to one, this structure
%                    contains results of the classical analysis: center (sample mean),
%                    cov (sample covariance matrix), md (Mahalanobis distances), class ('COV').
%        rew.class : 'MCDCOV'
%            rew.X : If x is bivariate, same as the x in the call to mcdcov,
%                    without rows containing missing or infinite values.
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust.html
%
% Written by Katrien Van Driessen and Bjorn Rombouts
% Revisions by Sanne Engelen, Sabine Verboven
% Last Update: 09/04/2004, 01/08/2007

% The FASTMCD algorithm works as follows:
%
%       The dataset contains n cases and p variables.
%       When n < 2*nmini (see below), the algorithm analyzes the dataset as a whole.
%       When n >= 2*nmini (see below), the algorithm uses several subdatasets.
%
%       When the dataset is analyzed as a whole, a trial subsample of p+1 cases
%       is taken, of which the mean and covariance matrix are calculated.
%       The h cases with smallest relative distances are used to calculate
%       the next mean and covariance matrix, and this cycle is repeated csteps1
%       times. For small n we consider all subsets of p+1 out of n, otherwise
%       the algorithm draws 500 random subsets by default.
%       Afterwards, the 10 best solutions (means and corresponding covariance
%       matrices) are used as starting values for the final iterations.
%       These iterations stop when two subsequent determinants become equal.
%       (At most csteps3 iteration steps are taken.) The solution with smallest
%       determinant is retained.
%
%       When the dataset contains more than 2*nmini cases, the algorithm does part
%       of the calculations on (at most) maxgroup nonoverlapping subdatasets, of
%       (roughly) maxobs cases.
%
%       Stage 1: For each trial subsample in each subdataset, csteps1 (see below) iterations are
%       carried out in that subdataset. For each subdataset, the 10 best solutions are
%       stored.
%
%       Stage 2 considers the union of the subdatasets, called the merged set.
%       (If n is large, the merged set is a proper subset of the entire dataset.)
%       In this merged set, each of the 'best solutions' of stage 1 are used as starting
%       values for csteps2 (sse below) iterations. Also here, the 10 best solutions are stored.
%
%       Stage 3 depends on n, the total number of cases in the dataset.
%       If n <= 5000, all 10 preliminary solutions are iterated.
%       If n > 5000, only the best preliminary solution is iterated.
%       The number of iterations decreases to 1 according to n*p (If n*p <= 100,000 we
%       iterate csteps3 (sse below) times, whereas for n*p > 1,000,000 we take only one iteration step).
%

if rem(nargin-1,2)~=0
    error('The number of input arguments should be odd!');
end
% Assigning some input parameters
data = x;
raw.cor = [];
rew.cor = [];
rew.plane = [];
% The maximum value for n (= number of observations) is:
% nmax=50000;
nmax=realmax;
% The maximum value for p (= number of variables) is:
pmax=50;
% To change the number of subdatasets and their size, the values of
% maxgroup and nmini can be changed.
maxgroup=5;
nmini=300;
% The number of iteration steps in stages 1,2 and 3 can be changed
% by adapting the parameters csteps1, csteps2, and csteps3.
csteps1=2;
csteps2=2;
csteps3=100;
% dtrial : number of subsamples if not all (p+1)-subsets will be considered.
dtrial=500;

if size(data,1)==1
    data=data';
end

% Observations with missing or infinite values are ommitted.
ok=all(isfinite(data),2);
data=data(ok,:);
xx=data;
[n,p]=size(data);
% Some checks are now performed.
if n==0
    error('All observations have missing or infinite values.')
end
if n > nmax
    error(['The program allows for at most ' int2str(nmax) ' observations.'])
end
if p > pmax
    error(['The program allows for at most ' int2str(pmax) ' variables.'])
end
if n < p
    error('Need at least (number of variables) observations.')
end

%internal variables
hmin=quanf(0.5,n,p);
%Assiging default values
h=quanf(0.75,n,p);
default=struct('alpha',0.75,'h',h,'plots',1,'ntrial',dtrial,'cor',0,'seed',0,'classic',0,'Hsets',[],'factor',0);
list=fieldnames(default);
options=default;
IN=length(list);
i=1;
counter=1;

%Reading optional inputarguments
if nargin>2
    %
    % placing inputfields in array of strings
    %
    for j=1:nargin-1
        if rem(j,2)~=0
            chklist{i}=varargin{j};
            i=i+1;
        end
    end
    dummy=sum(strcmp(chklist,'h')+2*strcmp(chklist,'alpha'));
    switch dummy
        case 0 % defaultvalues should be taken
            alfa=options.alpha;
            h=options.h;
        case 3
            error('Both input arguments alpha and h are provided. Only one is required.')
    end
    %
    % Checking which default parameters have to be changed
    % and keep them in the structure 'options'.
    %
    while counter<=IN
        index=strmatch(list(counter,:),chklist,'exact');
        if ~isempty(index) % in case of similarity
            for j=1:nargin-2 % searching the index of the accompanying field
                if rem(j,2)~=0 % fieldnames are placed on odd index
                    if strcmp(chklist{index},varargin{j})
                        I=j;
                    end
                end
            end
            options=setfield(options,chklist{index},varargin{I+1});
            index=[];
        end
        counter=counter+1;
    end
    if dummy==1% checking inputvariable h
        % hmin is the minimum number of observations whose covariance determinant
        % will be minimized.
        if isempty(options.Hsets)
            if options.h < hmin
                disp(['Warning: The MCD must cover at least ' int2str(hmin) ' observations.'])
                disp(['The value of h is set equal to ' int2str(hmin)])
                options.h = hmin;
            elseif options.h > n
                error('h is greater than the number of non-missings and non-infinites.')
            elseif options.h < p
                error(['h should be larger than the dimension ' int2str(p) '.'])
            end
        end
        options.alpha=options.h/n;
    elseif dummy==2
        if options.alpha < 0.5
            options.alpha=0.5;
            mess=sprintf(['Attention (mcdcov.m): Alpha should be larger than 0.5. \n',...
                'It is set to 0.5.']);
            disp(mess)
        end
        if options.alpha > 1
            options.alpha=0.75;
            mess=sprintf(['Attention (mcdcov.m): Alpha should be smaller than 1.\n',...
                'It is set to 0.75.']);
            disp(mess)
        end
        options.h=quanf(options.alpha,n,p);
    end
end

h=options.h;  %number of regular datapoints on which estimates are based. h=[alpha*n]
plots=options.plots; %relevant plots are plotted
alfa=options.alpha; %percentage of regular observations
ntrial=options.ntrial; %number of subsets to be taken in the first step
cor=options.cor; %correlation matrix
seed=options.seed; %seed of the random generator
cutoff.rd=sqrt(chi2inv(0.975,p)); %cutoff value for the robust distance
cutoff.md=cutoff.rd; %cutoff value for the mahalanobis distance
Hsets = options.Hsets;
if ~isempty(Hsets)
    Hsets_ind = 1;
else
    Hsets_ind = 0;
end
factor = options.factor;
if factor == 0
    factor_ind = 0;
else
    factor_ind = 1;
end

% Some initializations.
rew.flag=repmat(NaN,1,length(ok));
raw.wt=repmat(NaN,1,length(ok));
raw.rd=repmat(NaN,1,length(ok));
rew.rd=repmat(NaN,1,length(ok));
rew.mahalanobis=repmat(NaN,1,length(ok));
rew.method=sprintf('\nMinimum Covariance Determinant Estimator.');
correl=NaN;

%    z    : if at least h observations lie on a hyperplane, then z contains the
%           coefficients of that plane.
% weights : weights of the observations that are not excluded from the computations.
%           These are the observations that don't contain missing or infinite values.
% bestobj : best objective value found.
z(1:p)=0;
weights=zeros(1,n);
bestobj=inf;

% The classical estimates are computed.
[Pcl,Tcl,Lcl,rcl,centerXcl,cXcl] = classSVD(data);
clasmean=cXcl;
clascov=Pcl*diag(Lcl)*Pcl';

if p < 5
    eps=1e-12;
elseif p <= 8
    eps=1e-14;
else
    eps=1e-16;
end

% The standardization of the data will now be performed.
med=median(data);
mad=sort(abs(data-repmat(med,n,1)));
mad=mad(h,:);
ii=min(find(mad < eps));
if length(ii)
    % The h-th order statistic is zero for the ii-th variable. The array plane contains
    % all the observations which have the same value for the ii-th variable.
    plane=find(abs(data(:,ii)-med(ii)) < eps)';
    meanplane=mean(data(plane,:));
    weights(plane)=1;
    if p==1
        rew.flag=weights;
        raw.wt=weights;
        [raw.center,rew.center]=deal(meanplane);
        [raw.cov,rew.cov,raw.objective]=deal(0);
        if plots
            rew.method=sprintf('\nUnivariate location and scale estimation.');
            rew.method=strvcat(rew.method,sprintf('%g of the %g observations are identical.',length(plane),n));
        end
    else
        z(ii)=1;
        rew.plane=z;
        covplane=cov(data(plane,:));
        [raw.center,raw.cov,rew.center,rew.cov,raw.objective,raw.wt,rew.flag,...
            rew.method]=displ(3,length(plane),weights,n,p,meanplane,covplane,rew.method,z,ok,...
            raw.wt,rew.flag,0,NaN,h,ii);
    end
    rew.Hsubsets.Hopt = plane;
    rew.Hsubsets.Hfreq = plane;
    %classical analysis?
    if options.classic==1
        classic.cov=clascov;
        classic.center=clasmean;
        classic.class='COV';
    else
        classic=0;
    end
    %assigning the output
    rewo=rew;rawo=raw;
    rew=struct('center',{rewo.center},'cov',{rewo.cov},'cor',{cor},'h',{h},'Hsubsets',{rewo.Hsubsets},'alpha',{alfa},...
        'flag',{rewo.flag},'plane', {rewo.plane},'method',{rewo.method},'class',{'MCDCOV'},'classic',{classic},'X',{xx});
    raw=struct('center',{rawo.center},'cov',{rawo.cov},'cor',{rawo.cor},'objective',{rawo.objective},...
        'wt',{rawo.wt},'class',{'MCDCOV'},'classic',{classic},'X',{x});
    if size(data,2)~=2
        rew=rmfield(rew,'X');
        raw=rmfield(raw,'X');
    end
    return
end
data=(data-repmat(med,n,1))./repmat(mad,n,1);

% The standardized classical estimates are now computed.
clmean=mean(data);
clcov=cov(data);

% The univariate non-classical case is now handled.
if p==1 & h~=n
    [rew.center, rewsca, weights,raw.center,raw.cov,rawdist,raw.wt,Hopt]=unimcd(data,h);
    rew.Hsubsets.Hopt = Hopt';
    rew.Hsubsets.Hfreq = Hopt';
    raw.rd=sqrt(rawdist');
    rew.cov=rewsca^2;
    raw.objective=raw.cov*prod(mad)^2;
    mah=(data-rew.center).^2/rew.cov;
    rew.rd=sqrt(mah');
    rew.flag= rew.rd <= cutoff.rd;
    [raw.cov,raw.center]=trafo(raw.cov,raw.center,med,mad,p);
    [rew.cov,rew.center]=trafo(rew.cov,rew.center,med,mad,p);
    rew.mahalanobis=abs(data'-clmean)/sqrt(clcov);
    spec.ask=1;
    %classical analysis?
    if options.classic==1
        classic.cov=clascov;
        classic.center=clasmean;
        classic.md=rew.mahalanobis;
        classic.class='COV';
    else
        classic=0;
    end
    %assigning the output
    rewo=rew;rawo=raw;
    rew=struct('center',{rewo.center},'cov',{rewo.cov},'cor',{rewo.cor},'h',{h},'Hsubsets',{rewo.Hsubsets},...
        'alpha',{alfa},'rd',{rewo.rd},'cutoff',{cutoff},'flag',{rewo.flag}, 'plane',{[]},'method',{rewo.method},...
        'class',{'MCDCOV'},'md',{rewo.mahalanobis},'classic',{classic},'X',{xx});
    raw=struct('center',{rawo.center},'cov',{rawo.cov},'cor',{rawo.cor},'objective',{rawo.objective},...
        'rd',{rawo.rd},'cutoff',{cutoff},'wt',{rawo.wt}, 'class',{'MCDCOV'},'classic',{classic},'X',{x});
    if size(data,2)~=2
        rew=rmfield(rew,'X');
        raw=rmfield(raw,'X');
    end
    try
        if plots & options.classic
            makeplot(rew,'classic',1)
        elseif plots
            makeplot(rew)
        end
    catch %output must be given even if plots are interrupted
        %> delete(gcf) to get rid of the menu
        return
    end
    return
end
%exact fit situation
if rcl < p
    % all observations lie on a hyperplane.
    z = Pcl(:,1);
    rew.plane=z;
    weights(1:n)=1;
    if cor
        correl=clcov./(sqrt(diag(clcov))*sqrt(diag(clcov))');
    end
    [clcov,clmean]=trafo(clcov,clmean,med,mad,p);
    [raw.center,raw.cov,rew.center,rew.cov,raw.objective,raw.wt,rew.flag,...
        rew.method]=displ(1,n,weights,n,p,clmean,clcov,rew.method,z./mad',ok,...
        raw.wt,rew.flag,cor,correl);
    if cor
        [rew.cor,raw.cor]=deal(correl);
    end
    %classical analysis?
    if options.classic==1
        classic.cov=clascov;
        classic.center=clasmean;
        classic.class='COV';
    else
        classic=0;
    end
    rew.Hsubsets.Hopt=1:n;
    rew.Hsubsets.Hfreq=1:n;
    %assigning the output
    rewo=rew;rawo=raw;
    rew=struct('center',{rewo.center},'cov',{rewo.cov},'cor',{rewo.cor},'h',{h},'Hsubsets',{rewo.Hsubsets},'alpha',{alfa},...
        'rd',{rewo.rd},'cutoff',{cutoff},'flag',{rewo.flag},'plane',{rewo.plane},'method',{rewo.method},...
        'class',{'MCDCOV'},'classic',{classic},'X',{xx});
    raw=struct('center',{rawo.center},'cov',{rawo.cov},'cor',{rawo.cor},'objective',{rawo.objective},...
        'cutoff',{cutoff},'wt',{rawo.wt}, 'class',{'MCDCOV'},'classic',{classic},'X',{x});
    if size(data,2)~=2
        rew=rmfield(rew,'X');
        raw=rmfield(raw,'X');
    end
    return
end

% The classical case is now handled.
if h==n
    if plots
        msg=sprintf('The MCD estimates based on %g observations are equal to the classical estimates.\n',h);
        rew.method=strvcat(rew.method,msg);
    end
    raw.center=clmean;
    raw.cov=clcov;
    raw.objective=det(clcov);
    mah=libra_mahalanobis(data,clmean,'cov',clcov);
    rew.mahalanobis=sqrt(mah);
    raw.rd=rew.mahalanobis;
    weights= mah <= cutoff.rd^2;
    raw.wt=weights;
    [rew.center,rew.cov]=weightmecov(data,weights);
    if cor
        raw.cor=raw.cov./(sqrt(diag(raw.cov))*sqrt(diag(raw.cov))');
        rew.cor=rew.cov./(sqrt(diag(rew.cov))*sqrt(diag(rew.cov))');
    else
        raw.cor=0;
        rew.cor=0;
    end
    if det(rew.cov) < exp(-50*p)
        [center,covar,z,correl,plane,count]=fit(data,NaN,med,mad,p,z,cor,rew.center,rew.cov,n);
        rew.plane=z;
        if cor
            correl=covar./(sqrt(diag(cov))*sqrt(diag(covar))');
        end
        rew.method=displrw(count,n,p,center,covar,rew.method,z,cor,correl);
        [raw.cov,raw.center]=trafo(raw.cov,raw.center,med,mad,p);
        [rew.cov,rew.center]=trafo(rew.cov,rew.center,med,mad,p);
        rew.rd=raw.rd;
    else
        mah=libra_mahalanobis(data,rew.center,'cov',rew.cov);
        weights = mah <= cutoff.md^2;
        [raw.cov,raw.center]=trafo(raw.cov,raw.center,med,mad,p);
        [rew.cov,rew.center]=trafo(rew.cov,rew.center,med,mad,p);
        rew.rd=sqrt(mah);
    end
    raw.objective=raw.objective*prod(mad)^2;
    rew.flag=weights;
    %classical analysis?
    if options.classic==1
        classic.cov=clascov;
        classic.center=clasmean;
        classic.md=rew.mahalanobis;
        classic.class='COV';
    else
        classic=0;
    end
    %assigning Hsubsets:
    rew.Hsubsets.Hopt = 1:n;
    rew.Hsubsets.Hfreq = 1:n;
    %assigning the output
    rewo=rew;rawo=raw;
    rew=struct('center',{rewo.center},'cov',{rewo.cov},'cor',{rewo.cor},'h',{h},'Hsubsets',{rewo.Hsubsets},'alpha',{alfa},...
        'rd',{rewo.rd},'cutoff',{cutoff},'flag',{rewo.flag},'plane',{rewo.plane},...
        'method',{rewo.method},'class',{'MCDCOV'},'md',{rewo.mahalanobis},'classic',{classic},'X',{xx});
    raw=struct('center',{rawo.center},'cov',{rawo.cov},'cor',{rawo.cor},'objective',{rawo.objective},...
        'rd',{rawo.rd},'cutoff',{cutoff},'wt',{rawo.wt}, 'class',{'MCDCOV'},'classic',{classic},'X',{x});
    if size(data,2)~=2
        rew=rmfield(rew,'X');
        raw=rmfield(raw,'X');
    end
    try
        if plots & options.classic
            makeplot(rew,'classic',1)
        elseif plots
            makeplot(rew)
        end
    catch %output must be given even if plots are interrupted
        %> delete(gcf) to get rid of the menu
        return
    end
    return
end
percent=h/n;
teller = zeros(1,n+1);

if Hsets_ind
    csteps = csteps1;
    inplane = NaN;
    fine = 0;
    part = 0;
    final = 1;
    tottimes = 0;
    nsamp = size(Hsets,1);
    obsingroup = n;
else

    %  If n >= 2*nmini the dataset will be divided into subdatasets.  For n < 2*nmini the set
    %  will be treated as a whole.

    if n >= 2*nmini
        maxobs=maxgroup*nmini;
        if n >= maxobs
            ngroup=maxgroup;
            group(1:maxgroup)=nmini;
        else
            ngroup=floor(n/nmini);
            minquan=floor(n/ngroup);
            group(1)=minquan;
            for s=2:ngroup
                group(s)=minquan+double(rem(n,ngroup)>=s-1);
            end
        end
        part=1;
        adjh=floor(group(1)*percent);
        nsamp=floor(ntrial/ngroup);
        minigr=sum(group);
        obsingroup=fillgroup(n,group,ngroup,seed);
        % obsingroup : i-th row contains the observations of the i-th group.
        % The last row (ngroup+1-th) contains the observations for the 2nd stage
        % of the algorithm.
    else
        [part,group,ngroup,adjh,minigr,obsingroup]=deal(0,n,1,h,n,n);
        replow=[50,22,17,15,14,zeros(1,45)];
        if n < replow(p)
            % All (p+1)-subsets will be considered.
            al=1;
            perm=[1:p,p];
            nsamp=nchoosek(n,p+1);
        else
            al=0;
            nsamp=ntrial;
        end
    end
    % some further initialisations.

    csteps=csteps1;
    inplane=NaN;
    % tottimes : the total number of iteration steps.
    % fine     : becomes 1 when the subdatasets are merged.
    % final    : becomes 1 for the final stage of the algorithm.
    [tottimes,fine,final,prevdet]=deal(0);

    if part
        % bmean1 : contains, for the first stage of the algorithm, the means of the ngroup*10
        %          best estimates.
        % bcov1  : analogous to bmean1, but now for the covariance matrices.
        % bobj1  : analogous to bmean1, but now for the objective values.
        % coeff1 : if in the k-th subdataset there are at least adjh observations that lie on
        %          a hyperplane then the coefficients of this plane will be stored in the
        %          k-th column of coeff1.
        coeff1=repmat(NaN,p,ngroup);
        bobj1=repmat(inf,ngroup,10);
        bmean1=cell(ngroup,10);
        bP1 = cell(ngroup,10);
        bL1 = cell(ngroup,10);
        [bmean1{:}]=deal(NaN);
        [bL1{:}]=deal(NaN);
        [bP1{:}]=deal(NaN);
    end

    % bmean : contains the means of the ten best estimates obtained in the second stage of the
    %         algorithm.
    % bcov  : analogous to bmean, but now for the covariance matrices.
    % bobj  : analogous to bmean, but now for the objective values.
    % coeff : analogous to coeff1, but now for the merged subdataset.
    % If the data is not split up, the 10 best estimates obtained after csteps1 iterations
    % will be stored in bmean, bcov and bobj.
    coeff=repmat(NaN,p,1);
    bobj=repmat(inf,1,10);
    bmean=cell(1,10);
    bP = cell(1,10);
    bL = cell(1,10);
    [bmean{:}]=deal(NaN);
    [bP{:}]=deal(NaN);
    [bL{:}] = deal(NaN);
end

seed=0;
tellerwhilelus = 0;

while final~=2
    if fine | (~part & final)
        if ~Hsets_ind
            nsamp=10;
        end
        if final
            adjh=h;
            ngroup=1;
            if n*p <= 1e+5
                csteps=csteps3;
            elseif n*p <=1e+6
                csteps=10-(ceil(n*p/1e+5)-2);
            else
                csteps=1;
            end
            if n > 5000

                nsamp=1;
            end
        else
            adjh=floor(minigr*percent);
            csteps=csteps2;
        end
    end

    % found : becomes 1 if we have a singular intermediate MCD estimate.
    found=0;

    for k=1:ngroup
        if ~fine
            found=0;
        end
        for i=1:nsamp
            tottimes=tottimes+1;
            % ns becomes 1 if we have a singular trial subsample and if there are at
            % least adjh observations in the subdataset that lie on the concerning hyperplane.
            % In that case we don't have to take C-steps. The determinant is zero which is
            % already the lowest possible value. If ns=1, no C-steps will be taken and we
            % start with the next sample. If we, for the considered subdataset, haven't
            % already found a singular MCD estimate, then the results must be first stored in
            % bmean, bcov, bobj or in bmean1, bcov1 and bobj1.  If we, however, already found
            % a singular result for that subdataset, then the results won't be stored
            % (the hyperplane we just found is probably the same as the one we found earlier.
            % We then let adj be zero. This will guarantee us that the results won't be
            % stored) and we start immediately with the next sample.
            adj=1;
            ns=0;

            % For the second and final stage of the algorithm the array sortdist(1:adjh)
            % contains the indices of the observations corresponding to the adjh observations
            % with minimal relative distances with respect to the best estimates of the
            % previous stage. An exception to this, is when the estimate of the previous
            % stage is singular.  For the second stage we then distinguish two cases :
            %
            % 1. There aren't adjh observations in the merged set that lie on the hyperplane.
            %    The observations on the hyperplane are then extended to adjh observations by
            %    adding the observations of the merged set with smallest orthogonal distances
            %    to that hyperplane.
            % 2. There are adjh or more observations in the merged set that lie on the
            %    hyperplane. We distinguish two cases. We haven't or have already found such
            %    a hyperplane. In the first case we start with a new sample.  But first, we
            %    store the results in bmean1, bcov1 and bobj1. In the second case we
            %    immediately start with a new sample.
            %
            % For the final stage we do the same as 1. above (if we had h or more observations
            % on the hyperplane we would already have found it).

            if ~Hsets_ind

                if final
                    if ~isinf(bobj(i))
                        meanvct=bmean{i};
                        P=bP{i};
                        L = bL{i};
                        if bobj(i)==0
                            [dis,sortdist]=sort(abs(sum((data-repmat(meanvct,n,1))'.*repmat(coeff,1,n))));
                        else
                            [dis,sortdist]=mahal2((data-repmat(meanvct,n,1))*P,sqrt(L),part,fine,final,k,obsingroup);
                        end
                    else
                        break
                    end
                elseif fine
                    if ~isinf(bobj1(k,i))
                        meanvct=bmean1{k,i};
                        P=bP1{k,i};
                        L=bL1{k,i};
                        if bobj1(k,i)==0
                            [dis,ind]=sort(abs(sum((data(obsingroup{end},:)-repmat(meanvct,minigr,1))'.*repmat(coeff1(:,k),1,minigr))));
                            sortdist=obsingroup{end}(ind);
                            if dis(adjh) < 1e-8
                                if found==0
                                    obj=0;
                                    coeff=coeff1(:,k);
                                    found=1;
                                else
                                    adj=0;
                                end
                                ns=1;
                            end
                        else
                            [dis,sortdist]=mahal2((data(obsingroup{end},:)-repmat(meanvct,minigr,1))*P,sqrt(L),part,fine,final,k,obsingroup);
                        end
                    else
                        break;
                    end
                else
                    if ~part
                        if al
                            k=p+1;
                            perm(k)=perm(k)+1;
                            while ~(k==1 |perm(k) <=(n-(p+1-k)))
                                k=k-1;
                                perm(k)=perm(k)+1;
                                for j=(k+1):p+1
                                    perm(j)=perm(j-1)+1;
                                end
                            end
                            index=perm;
                        else
                            [index,seed]=randomset(n,p+1,seed);
                        end
                    else
                        [index,seed]=randomset(group(k),p+1,seed);
                        index=obsingroup{k}(index);
                    end

                    [P,T,L,r,centerX,meanvct] = classSVD(data(index,:));
                    if r==0
                        ns = 1;
                        rew.center=meanvct;
                        rew.cov=zeros(p,p);
                        rew.flag = zeros(1,n);
                        rew.flag(index) = 1;
                        rew.Hsubsets.Hopt=1:n;
                        rew.Hsubsets.Hfreq=1:n;
                    elseif (r < p) & (r~=0)
                        % The trial subsample is singular.
                        % We distinguish two cases :
                        %
                        % 1. There are adjh or more observations in the subdataset that lie
                        %    on the hyperplane. If the data is not split up, we have adjh=h and thus
                        %    an exact fit. If the data is split up we distinguish two cases.
                        %    We haven't or have already found such a hyperplane.  In the first case
                        %    we check if there are more than h observations in the entire set
                        %    that lie on the hyperplane. If so, we have an exact fit situation.
                        %    If not, we start with a new trial subsample.  But first, the
                        %    results must be stored bmean1, bcov1 and bobj1.  In the second case
                        %    we immediately start with a new trial subsample.
                        %
                        % 2. There aren't adjh observations in the subdataset that lie on the
                        %    hyperplane. We then extend the trial subsample until it isn't singular
                        %    anymore.


                        % eigvct : contains the coefficients of the hyperplane.

                        eigvct = P(:,1);

                        if ~part
                            dist=abs(sum((data-repmat(meanvct,n,1))'.*repmat(eigvct,1,n)));
                        else
                            dist=abs(sum((data(obsingroup{k},:)-repmat(meanvct,group(k),1))'.*repmat(eigvct,1,group(k))));
                        end

                        obsinplane=find(dist < 1e-8);
                        % count : number of observations that lie on the hyperplane.
                        count=length(obsinplane);

                        if count >= adjh
                            if ~part
                                [center,covar,eigvct,correl]=fit(data,obsinplane,med,mad,p,eigvct,cor);
                                rew.plane=eigvct;
                                weights(obsinplane)=1;
                                [raw.center,raw.cov,rew.center,rew.cov,raw.objective,...
                                    raw.wt,rew.flag,rew.method]=displ(2,count,weights,n,p,center,covar,...
                                    rew.method,eigvct,ok,raw.wt,rew.flag,cor,correl);
                                if cor
                                    [rew.cor,raw.cor]=deal(correl);
                                end
                                rew.Hsubsets.Hopt=obsinplane;
                                rew.Hsubsets.Hfreq=obsinplane;
                                return
                            elseif found==0
                                if count>=h
                                    [center,covar,eigvct,correl]=fit(data,obsinplane,med,mad,p,eigvct,cor);
                                    rew.plane=eigvct;
                                    weights(obsinplane)=1;
                                    [raw.center,raw.cov,rew.center,rew.cov,raw.objective,...
                                        raw.wt,rew.flag,rew.method,varargout]=displ(2,count2,weights,n,p,center,covar,...
                                        rew.method,eigvct,ok,raw.wt,rew.flag,cor,correl);
                                    if cor
                                        [rew.cor,raw.cor]=deal(correl);
                                    end
                                    rew.Hsubsets.Hopt=obsinplane;
                                    rew.Hsubsets.Hfreq=obsinplane;
                                    return
                                end
                                obj=0;
                                inplane(k)=count;
                                coeff1(:,k)=eigvct;
                                found=1;
                                ns=1;
                            else
                                ns=1;
                                adj=0;
                            end
                        else
                            covmat = cov(data(index,:));
                            meanvct = mean(data(index,:));
                            while det(covmat) < exp(-50*p)
                                [index1,seed]=addobs(index,n,seed);
                                [covmat,meanvct] = updatecov(data(index,:),covmat,meanvct,data(setdiff(index1,index),:),[],1);
                                index = index1;
                            end
                        end
                    end

                    if ~ns
                        if ~part
                            [dis,sortdist] = mahal2((data - repmat(meanvct,n,1))*P,sqrt(L),part,fine,final,k,obsingroup);
                        else
                            [dis,sortdist] = mahal2((data(obsingroup{k},:) - repmat(meanvct,group(k),1))*P,sqrt(L),part,fine,final,k,obsingroup);
                        end
                    end
                end
            end

            if ~ns
                for j=1:csteps
                    tottimes=tottimes+1;
                    if j == 1
                        if Hsets_ind
                            obs_in_set = Hsets(i,:);
                        else
                            obs_in_set = sort(sortdist(1:adjh));
                            teller(obs_in_set) = teller(obs_in_set) + 1;
                            teller(end) = teller(end) + 1;
                        end
                    else
                        % The observations correponding to the adjh smallest mahalanobis
                        % distances determine the subset for the next iteration.
                        if ~part
                            [dis2,sortdist] = mahal2((data - repmat(meanvct,n,1))*P,sqrt(L),part,fine,final,k,obsingroup);
                        else
                            if final
                                [dis2,sortdist] = mahal2((data - repmat(meanvct,n,1))*P,sqrt(L),part,fine,final,k,obsingroup);
                            elseif fine
                                [dis2,sortdist] = mahal2((data(obsingroup{end},:) - repmat(meanvct,minigr,1))*P,sqrt(L),part,fine,final,k,obsingroup);
                            else
                                [dis2,sortdist] = mahal2((data(obsingroup{k},:) - repmat(meanvct,group(k),1))*P,sqrt(L),part,fine,final,k,obsingroup);
                            end
                        end
                        % Creation of a H-subset.
                        obs_in_set=sort(sortdist(1:adjh));
                        teller(obs_in_set) = teller(obs_in_set) + 1;
                        teller(end) = teller(end) + 1;
                    end
                    [P,T,L,r,centerX,meanvct] = classSVD(data(obs_in_set,:));
                    if r==0
                        rew.center=meanvct;
                        rew.cov=zeros(p,p);
                        rew.flag=(data==data(obs_in_set(1),:));
                        rew.Hsubset.Hopt=obs_in_set(1);
                        rew.Hsubset.Hfreq=obs_in_set(1);
                    else

                        obj=prod(L);

                        if obj < exp(-50*p)
                            % The adjh-subset is singular. If adjh=h we have an exact fit situation.
                            % If adjh < h we distinguish two cases :
                            %
                            % 1. We haven't found earlier a singular adjh-subset. We first check if
                            %    in the entire set there are h observations that lie on the hyperplane.
                            %    If so, we have an exact fit situation. If not, we stop taking C-steps
                            %    (the determinant is zero which is the lowest possible value) and
                            %    store the results in the appropriate arrays.  We then begin with
                            %    the next trial subsample.
                            %
                            % 2. We have, for the concerning subdataset, already found a singular
                            %    adjh-subset. We then immediately begin with the next trial subsample.

                            if ~part | final | (fine & n==minigr)
                                covmat=cov(data(obs_in_set,:));
                                [center,covar,z,correl,obsinplane,count]=fit(data,NaN,med,mad,p,NaN,...
                                    cor,meanvct,covmat,n);
                                rew.plane=z;
                                weights(obsinplane)=1;
                                [raw.center,raw.cov,rew.center,rew.cov,raw.objective,...
                                    raw.wt,rew.flag,rew.method]=displ(2,count,weights,n,p,center,covar,...
                                    rew.method,z,ok,raw.wt,rew.flag,cor,correl);
                                if cor
                                    [rew.cor,raw.cor]=deal(correl);
                                end
                                rew.Hsubset.Hopt=obsinplane;
                                rew.Hsubset.Hfreq=obsinplane;
                                return
                            elseif found==0
                                eigvct = V(:,1);
                                dist=abs(sum((data-repmat(meanvct,n,1))'.*repmat(eigvct,1,n)));
                                obsinplane=find(dist<1e-8);
                                count=length(obsinplane);
                                if count >= h
                                    [center,covar,eigvct,correl]=fit(data,obsinplane,med,mad,p,eigvct,cor);
                                    rew.plane=eigvct;
                                    weights(obsinplane)=1;
                                    [raw.center,raw.cov,rew.center,rew.cov,raw.objective,...
                                        raw.wt,rew.flag,rew.method]=displ(2,count,weights,n,p,center,covar,...
                                        rew.method,eigvct,ok,raw.wt,rew.flag,cor,correl);
                                    if cor
                                        [rew.cor,raw.cor]=deal(correl);
                                    end
                                    rew.Hsubset.Hopt=obsinplane;
                                    rew.Hsubset.Hfreq=obsinplane;
                                    return
                                end
                                obj=0;
                                found=1;
                                if ~fine
                                    coeff1(:,k)=eigvct;
                                    dist=abs(sum((data(obsingroup{k},:)-repmat(meanvct,group(k),1))'.*repmat(eigvct,1,group(k))));
                                    inplane(k)=length(dist(dist<1e-8));
                                else
                                    coeff=eigvct;
                                    dist=abs(sum((data(obsingroup{end},:)-repmat(meanvct,minigr,1))'.*repmat(eigvct,1,minigr)));
                                    inplane=length(dist(dist<1e-8));
                                end
                                break;
                            else
                                adj=0;
                                break;
                            end
                        end
                    end
                    % We stop taking C-steps when two subsequent determinants become equal.
                    % We have then reached convergence.
                    if j >= 2 & obj == prevdet
                        break;
                    end
                    prevdet=obj;

                end % C-steps

            end

            % After each iteration, it has to be checked whether the new solution
            % is better than some previous one.  A distinction is made between the
            % different stages of the algorithm:
            %
            %  - Let us first consider the first (second) stage of the algorithm.
            %    We distinguish two cases if the objective value is lower than the largest
            %    value in bobj1 (bobj) :
            %
            %      1. The new objective value did not yet occur in bobj1 (bobj).  We then store
            %         this value, the corresponding mean and covariance matrix at the right
            %         place in resp. bobj1 (bobj), bmean1 (bmean) and bcov1 (bcov).
            %         The objective value is inserted by shifting the greater determinants
            %         upwards. We perform the same shifting in bmean1 (bmean) and bcov1 (bcov).
            %
            %      2. The new objective value already occurs in bobj1 (bobj). A comparison is
            %         made between the new mean vector and covariance matrix and those
            %         estimates with the same determinant. When for an equal determinant,
            %         the mean vector or covariance matrix do not correspond, the new results
            %         will be stored in bobj1 (bobj), bmean1 (bmean) and bcov1 (bcov).
            %
            %    If the objective value is not lower than the largest value in bobj1 (bobj),
            %    nothing happens.
            %
            %  - For the final stage of the algorithm, only the best solution has to be kept.
            %    We then check if the objective value is lower than the till then lowest value.
            %    If so, we have a new best solution. If not, nothing happens.


            if ~final & adj
                if fine | ~part
                    if obj < max(bobj) & ~ns
                        [bmean,bP,bL,bobj]=insertion(bmean,bP,bL,bobj,meanvct,P,L,obj,1,eps);
                    end
                else
                    if obj < max(bobj1(k,:)) & ~ns
                        [bmean1,bP1,bL1,bobj1]=insertion(bmean1,bP1,bL1,bobj1,meanvct,P,L,obj,k,eps);
                    end
                end
            end

            if final & obj< bestobj
                % bestset           : the best subset for the whole data.
                % bestobj           : objective value for this set.
                % initmean, initcov : resp. the mean and covariance matrix of this set.
                bestset=obs_in_set;
                bestobj=obj;
                initmean=meanvct;
                initcov=cov(data(bestset,:));
                raw.initcov=cov(data(bestset,:));
            end

        end % nsamp
    end % ngroup


    if part & ~fine
        fine=1;
    elseif (part & fine & ~final) | (~part & ~final)
        final=1;
    else
        final=2;
    end

end % while loop

[P,T,L,r,centerX,cX] = classSVD(data(bestset,:));
mah=libra_mahalanobis((data - repmat(cX,n,1))*P,zeros(size(P,2),1),'cov',L);
sortmah=sort(mah);

[sortset,indbestset] = sort(mah(bestset));
sortbestset = bestset(indbestset);
rew.Hsubsets.Hopt = sortbestset;

if ~factor_ind
    factor = sortmah(h)/chi2inv(h/n,p);
else
    factor = sortmah(h)/chi2inv(h/n,p/2);
end

raw.cov=factor*initcov;
% We express the results in the original units.
[raw.cov,raw.center]=trafo(raw.cov,initmean,med,mad,p);
raw.objective=bestobj*prod(mad)^2;

if cor
    raw.cor=raw.cov./(sqrt(diag(raw.cov))*sqrt(diag(raw.cov))');
end

% the mahalanobis distances are computed without the factor, therefore we
% have to correct for it now.
mah=mah/factor;
raw.rd=sqrt(mah);
weights=mah<=cutoff.md^2;
raw.wt=weights;
[rew.center,rew.cov]=weightmecov(data,weights);
[trcov,trcenter]=trafo(rew.cov,rew.center,med,mad,p);

% determination of Hfreq:
[telobs,indobs] = greatsort(teller(1:(end - 1)));
rew.Hsubsets.Hfreq = indobs(1:(h));
if size(rew.Hsubsets.Hfreq,2) == 1
    rew.Hsubsets.Hfreq = rew.Hsubsets.Hfreq';
end

if cor
    rew.cor=rew.cov./(sqrt(diag(rew.cov))*sqrt(diag(rew.cov))');
end

if prod(sqrt(L)) < exp(-50*p)
    [center,covar,z,correl,plane,count]=fit(data,NaN,med,mad,p,z,cor,rew.center,rew.cov,n);
    rew.plane=z;
    if cor
        correl=covar./(sqrt(diag(covar))*sqrt(diag(covar))');
        [rew.cor,raw.cor] = deal(correl);
    end
    rew.method=displrw(count,n,p,center,covar,rew.method,z,cor,correl);
    rew.flag=weights;
    rew.rd=raw.rd;
else
    mah=libra_mahalanobis(data,rew.center,'cov',rew.cov);
    rew.flag=(mah <= cutoff.md^2);
    rew.rd=sqrt(mah);
end

rew.mahalanobis=sqrt(libra_mahalanobis(data,clmean,'cov',clcov));
rawo=raw;
reso=rew;
if options.classic==1
    classic.cov=clascov;
    classic.center=clasmean;
    classic.md=rew.mahalanobis;
    classic.flag = (classic.md <= cutoff.md);
    if options.cor==1
        classic.cor=clascov./(sqrt(diag(clascov))*sqrt(diag(clascov))');
    end
    classic.class='COV';
else
    classic=0;
end
%assigning the output
rew=struct('center',{trcenter},'cov',{trcov},'cor',{reso.cor},'h',{h},'Hsubsets',{reso.Hsubsets},'alpha',{alfa},...
    'rd',{reso.rd},'flag',{reso.flag},'md',{reso.mahalanobis},'cutoff',{cutoff},...
    'plane',{reso.plane},'method',{reso.method},'class',{'MCDCOV'},'classic',{classic},'X',{xx});
raw=struct('center',{rawo.center},'cov',{rawo.cov},'cor',{rawo.cor},'objective',{rawo.objective},...
    'rd',{rawo.rd},'cutoff',{cutoff},...
    'wt',{rawo.wt},'class',{'MCDCOV'},'classic',{classic},'X',{x});

if size(data,2)~=2
    rew=rmfield(rew,'X');
    raw=rmfield(raw,'X');
end

try
    if plots & options.classic
        makeplot(rew,'classic',1)
    elseif plots
        makeplot(rew)
    end
catch %output must be given even if plots are interrupted
    %> delete(gcf) to get rid of the menu
end

%-----------------------------------------------------------------------------------------
function [raw_center,raw_cov,center,covar,raw_objective,raw_wt,mcd_wt,method]=displ(exactfit,...
    count,weights,n,p,center,covar,method,z,ok,raw_wt,mcd_wt,cor,correl,varargin)

% Determines some fields of the output argument REW for the exact fit situation.  It also
% displays and writes the messages concerning the exact fit situation.  If the raw MCD
% covariance matrix is not singular but the reweighted is, then the function displrw is
% called instead of this function.

[raw_center,center]=deal(center);
[raw_cov,cover]=deal(covar);
raw_objective=0;
mcd_wt=weights;
raw_wt=weights;

switch exactfit
    case 1
        msg='The covariance matrix of the data is singular.';
    case 2
        msg='The covariance matrix has become singular during the iterations of the MCD algorithm.';
    case 3
        msg=sprintf('The %g-th order statistic of the absolute deviation of variable %g is zero. ',varargin{1},varargin{2});
end

msg=sprintf([msg '\nThere are %g observations in the entire dataset of %g observations that lie on the \n'],count,n);
switch p
    case 2
        msg=sprintf([msg 'line with equation %g(x_i1-m_1)%+g(x_i2-m_2)=0 \n'],z);
        msg=sprintf([msg 'where the mean (m_1,m_2) of these observations is the MCD location']);
    case 3
        msg=sprintf([msg 'plane with equation %g(x_i1-m_1)%+g(x_i2-m_2)%+g(x_i3-m_3)=0 \n'],z);
        msg=sprintf([msg 'where the mean (m_1,m_2,m_3) of these observations is the MCD location']);
    otherwise
        msg=sprintf([msg 'hyperplane with equation a_1 (x_i1-m_1) + ... + a_p (x_ip-m_p) = 0 \n']);
        msg=sprintf([msg 'with coefficients a_i equal to : \n\n']);
        msg=sprintf([msg sprintf('%g  ',z)]);
        msg=sprintf([msg '\n\nand where the mean (m_1,...,m_p) of these observations is the MCD location']);
end

method=strvcat(method,[msg '.']);
disp(method)

%-----------------------------------------------------------------------------------------
function method=displrw(count,n,p,center,covar,method,z,cor,correl)

% Displays and writes messages in the case the reweighted robust covariance matrix
% is singular.

msg=sprintf('The reweighted MCD scatter matrix is singular. \n');
msg=sprintf([msg 'There are %g observations in the entire dataset of %g observations that lie on the\n'],count,n);

switch p
    case 2
        msg=sprintf([msg 'line with equation %g(x_i1-m_1)%+g(x_i2-m_2)=0 \n\n'],z);
        msg=sprintf([msg 'where the mean (m_1,m_2) of these observations is : \n\n']);
    case 3
        msg=sprintf([msg 'plane with equation %g(x_i1-m_1)%+g(x_i2-m_2)%+g(x_i3-m_3)=0 \n\n'],z);
        msg=sprintf([msg 'where the mean (m_1,m_2,m_3) of these observations is : \n\n']);
    otherwise
        msg=sprintf([msg 'hyperplane with equation a_1 (x_i1-m_1) + ... + a_p (x_ip-m_p) = 0 \n']);
        msg=sprintf([msg 'with coefficients a_i equal to : \n\n']);
        msg=sprintf([msg sprintf('%g  ',z)]);
        msg=sprintf([msg '\n\nand where the mean (m_1,...,m_p) of these observations is : \n\n']);
end

msg=sprintf([msg sprintf('%g  ',center)]);
msg=sprintf([msg '\n\nTheir covariance matrix equals : \n\n']);
msg=sprintf([msg sprintf([repmat('% 13.5g ',1,p) '\n'],covar)]);
if cor
    msg=sprintf([msg '\n\nand their correlation matrix equals : \n\n']);
    msg=sprintf([msg sprintf([repmat('% 13.5g ',1,p) '\n'],correl)]);
end

method=strvcat(method,msg);

%-----------------------------------------------------------------------------------------
function [initmean,initcov,z,correl,varargout]=fit(dat,plane,med,mad,p,z,cor,varargin)

% This function is called in the case of an exact fit. It computes the correlation
% matrix and transforms the coefficients of the hyperplane, the mean, the covariance
% and the correlation matrix to the original units.

if isnan(plane)
    [meanvct,covmat,n]=deal(varargin{:});
    [z, eigvl]=eigs(covmat,1,0,struct('disp',0));
    dist=abs(sum((dat-repmat(meanvct,n,1))'.*repmat(z,1,n)));
    plane=find(dist < 1e-8);
    varargout{1}=plane;
    varargout{2}=length(plane);
end

z=z./mad';
[initcov,initmean]=trafo(cov(dat(plane,:)),mean(dat(plane,:)),med,mad,p);
if cor
    correl=initcov./(sqrt(diag(initcov))*sqrt(diag(initcov))');
else
    correl=NaN;
end

%------------------------------------------------------------------------------------------
function obsingroup=fillgroup(n,group,ngroup,seed)

% Creates the subdatasets.

obsingroup=cell(1,ngroup+1);

jndex=0;
for k=1:ngroup
    for m=1:group(k)
        [random,seed]=uniran(seed);
        ran=floor(random*(n-jndex)+1);
        jndex=jndex+1;
        if jndex==1
            index(1,jndex)=ran;
            index(2,jndex)=k;
        else
            index(1,jndex)=ran+jndex-1;
            index(2,jndex)=k;
            ii=min(find(index(1,1:jndex-1) > ran-1+[1:jndex-1]));
            if length(ii)
                index(1,jndex:-1:ii+1)=index(1,jndex-1:-1:ii);
                index(2,jndex:-1:ii+1)=index(2,jndex-1:-1:ii);
                index(1,ii)=ran+ii-1;
                index(2,ii)=k;
            end
        end
    end
    obsingroup{k}=index(1,index(2,:)==k);
    obsingroup{ngroup+1}=[obsingroup{ngroup+1},obsingroup{k}];
end

%-----------------------------------------------------------------------------------------
function [ranset,seed]=randomset(tot,nel,seed)

% This function is called if not all (p+1)-subsets out of n will be considered.
% It randomly draws a subsample of nel cases out of tot.

for j=1:nel
    [random,seed]=uniran(seed);
    num=floor(random*tot)+1;
    if j > 1
        while any(ranset==num)
            [random,seed]=uniran(seed);
            num=floor(random*tot)+1;
        end
    end
    ranset(j)=num;
end

%-----------------------------------------------------------------------------------------
function [index,seed]=addobs(index,n,seed)

% Extends a trial subsample with one observation.

jndex=length(index);
[random,seed]=uniran(seed);
ran=floor(random*(n-jndex)+1);
jndex=jndex+1;
index(jndex)=ran+jndex-1;
ii=min(find(index(1:jndex-1) > ran-1+[1:jndex-1]));
if length(ii)~=0
    index(jndex:-1:ii+1)=index(jndex-1:-1:ii);
    index(ii)=ran+ii-1;
end

%-----------------------------------------------------------------------------------------

function mahsort=mahal(dat,meanvct,covmat,part,fine,final,k,obsingroup,group,minigr,n,nvar)

% Orders the observations according to the mahalanobis distances.

if ~part | final
    [dis,ind]=sort(libra_mahalanobis(dat,meanvct,'cov',covmat));
    mahsort=ind;
elseif fine
    [dis,ind]=sort(libra_mahalanobis(dat(obsingroup{end},:),meanvct,'cov',covmat));
    mahsort=obsingroup{end}(ind);
else
    [dis,ind]=sort(libra_mahalanobis(dat(obsingroup{k},:),meanvct,'cov',covmat));
    mahsort=obsingroup{k}(ind);
end

%-----------------------------------------------------------------------------------------

function [dis,mahsort]=mahal2(score,sca,part,fine,final,k,obsingroup)

% Orders the observations according to the mahalanobis distances for a diagonal
% covariance matrix and zero mean. sca contains the squareroot of the diagonal elements.

if ~part | final
    [dis,ind]=sort(libra_mahalanobis(score,zeros(size(score,2),1),'cov',sca.^2));
    mahsort=ind;
elseif fine
    [dis,ind]=sort(libra_mahalanobis(score,zeros(size(score,2),1),'cov',sca.^2));
    mahsort=obsingroup{end}(ind);
else
    [dis,ind]=sort(libra_mahalanobis(score,zeros(size(score,2),1),'cov',sca.^2));
    mahsort=obsingroup{k}(ind);
end

%------------------------------------------------------------------------------------------

function [covmat,meanvct]=trafo(covmat,meanvct,med,mad,nvar)

% Transforms a mean vector and a covariance matrix to the original units.

covmat=covmat.*repmat(mad,nvar,1).*repmat(mad',1,nvar);
meanvct=meanvct.*mad+med;

%-----------------------------------------------------------------------------------------
function [bestmean,bestP,bestL,bobj]=insertion(bestmean,bestP,bestL,bobj,meanvct,P,L,obj,row,eps)

% Stores, for the first and second stage of the algorithm, the results in the appropriate
% arrays if it belongs to the 10 best results.

insert=1;

equ=find(obj==bobj(row,:));

for j=equ
    if (meanvct==bestmean{row,j})
        if all(P==bestP{row,j})
            if all(L==bestL{row,j})
                insert=0;
            end
        end
    end
end

if insert
    ins=min(find(obj < bobj(row,:)));

    if ins==10
        bestmean{row,ins}=meanvct;
        bestP{row,ins} = P;
        bestL{row,ins} = L;
        bobj(row,ins)=obj;
    else
        [bestmean{row,ins+1:10}]=deal(bestmean{row,ins:9});
        bestmean{row,ins}=meanvct;
        [bestP{row,ins+1:10}] = deal(bestP{row,ins:9});
        bestP{row,ins} = P;
        [bestL{row,ins+1:10}] = deal(bestL{row,ins:9});
        bestL{row,ins} = L;
        bobj(row,ins+1:10)=bobj(row,ins:9);
        bobj(row,ins)=obj;
    end

end

%-----------------------------------------------------------------------------------------
function quan=quanf(alfa,n,rk)
quan=floor(2*floor((n+rk+1)/2)-n+2*(n-floor((n+rk+1)/2))*alfa);
%--------------------------------------------------------------------------