function [idx,dm,mm,Ss] = kur_rce(x,mode)

%
% [idx,dm,mean,Cov] = kur_rce(x,mode)
%
% Outlier identification based on the study of univariate
% projections onto the directions maximizing and minimizing
% the kurtosis coefficient
%
% Inputs:   observations, x, (a matrix having a row for each observation)
%           mode, 0   use p maxizing and p minimizing directions
%                     iterate until no outliers are detected (default)
%                 >0  use p maxizing and p minimizing directions
%                     iterate "mode" times
%                 <0  use only maximizing directions
%                     iterate "mode" times
% Outputs:  idx, zero/one vector with ones in the positions of the
%                suspected outliers
%           dm,  robust Mahalanobis distances
%           mean, robust mean estimator
%           Ss,  robust covariance matrix estimator
%

% Daniel Pena / Francisco J Prieto 23/5/00

[n,p] = size(x);

%Quite esta restricción por obsoleta
% if (n > 750)|(p > 50),
%   disp('Data set too large for this code');
%   return
% end

if nargin < 2,
  mode = 0;
end
maxit = Inf;
if mode < 0,
  maxit = abs(mode);
  mode = -1;
elseif mode > 0,
  maxit = mode;
  mode = 1;
end

% Initialization of parameters

%% Cutoff points for projections

ctf1 = 3.05 + 0.324*p;
ctf2 = exp(1.522-0.268*log(p));
ctf3 = 3.5;
ctf4 = 3.0;

%% minimum numbers of observations

lmn1 = max(floor(0.5*(n+p+1)),2*p);
lmn2 = min(2*p,max(p+2,floor(0.25*n)));

% Removing suspect outliers

xx0 = [x (1:n)'];
xx = xx0;
nn = n;
it = 1;
while (it <= maxit),
  xx1 = xx(:,1:p);
  V = kur_nw(xx1,mode);
  pr = xx1*V;
  md = median(pr);
  prn = abs(pr - ones(nn,1)*md);
%   mad = median(prn);

  %Problema con el estimador robusto de varianza (mad) al tener muchos valores en 0 o haber distribuciones con poca varianza. 
  %Solución: Si alguna dimension es muy pequeña, tomo la media de entre los percentiles 50 y 75 de prn en esa
  %dimensión como estimacion del desvio std.
  prn_percentiles = prctile(prn, [50 95]);
  
  if( size(prn,2) == 1 )
      prn_percentiles = prn_percentiles';
  end
  
  mad = prn_percentiles(1,:);
  
  disp_ratio = prn_percentiles(2,:) ./ prn_percentiles(1,:);
  
  if( any(disp_ratio > 5) )
      % Para una distribucion normal, en 3 sigma encontramos el 99% de los
      % puntos, por este motivos pusimos el umbral en 5 para ver cuando la
      % estimacion es muy mala.
    dim_with_low_var = find(mad == 0);
    mad(dim_with_low_var) = mean( prctile( prn(:,dim_with_low_var), [50 75] ) );

    if( any(mad(dim_with_low_var) == 0) )
        
%         disp(['No se pudo estimar la varianza en las dimensiones ' num2str(dim_with_low_var) '. Posiblemente esos features estén mal calculados.' ]);
        
        %No se indica ningun outlier en este caso.
        idx = false(n,1);
        
        return
        
    end
    
  end
  
  tt = prn./(ones(nn,1)*mad);

  if p > 2,
    ctfo = ctf1 - (0:(p-2))*0.75*(ctf1-ctf2)/(p-2);
  else
    ctfo = ctf1;
  end

  if mode >= 0,
    tctf = ones(nn,1)*[ ctfo ctf2 ctf3*ones(1,p-1) ctf4 ];
  else
    if( p > 1)
        tctf = ones(nn,1)*[ ctfo ctf2 ];
    else
        tctf = ones(nn,1)*[ ctfo ];
    end
  end
  
  taux = tt./tctf;
  t = max(taux,[],2);

%   in = find(t > 1);
% Meti una ampliación del nivel de rechazo para outliers ya que parecía muy
% sensible.
  in = find(t > 5);
  nn = length(in);
  if nn == 0,
    break
  end
%   inn = find(t <= 1);
% Meti una ampliación del nivel de rechazo para outliers ya que parecía muy
% sensible.
  inn = find(t <= 5);
  nu = length(inn);

  if nu < lmn2,
    [v,ix] = sort(t);
    ixx = ix(1:lmn2);
    xx = xx(ixx,:);
    break
  end
  xx = xx(inn,:);
  [nn,pp] = size(xx);
  if nn <= lmn2,
    break
  end
  it = it + 1;
end

%% Extracting the indices of the outliers

idx = ones(n,1);
[nn,pp] = size(xx);
idx(xx(:,pp)) = zeros(nn,1);

% recheck observations and relabel them

%% Cutoffs for Mahalanobis distances, based on the Hotelling-t2
%% Fixed at 0.99 . To modifay change the value of  ctflvl

ctflvl = 0.99;
ctf = (p*(n-1)/(n-p))*finv(ctflvl,p,n-p);

%% Mahalanobis distances using center and scale estimators
%% based on good observations

sidx = sum(idx);
s1 = find(idx);
s2 = find(idx == 0);
if sidx > 0,
  xx1 = xx0(s1,:);
  xx1r = xx1(:,1:p);
end
xx2r = xx0(s2,1:p);

mm = mean(xx2r);
Ss = cov(xx2r);
if sidx > 0,
  aux1 = xx1r - ones(sidx,1)*mm;
  %Quito un poco del ruido que genera esta inversion de matriz
  warning off
  dd = Ss\aux1';
  warning on
  dms = sum((aux1.*dd')');
end

%% Ensure that at least lmn1 observations are considered good

ado = sidx + lmn1 - n;
if ado > 0,
  [dmss,idms] = sort(dms);
  idm1 = idms(1:ado);
  ido = xx0(s1,pp);
  s3 = ido(idm1);
  idx(s3) = zeros(ado,1);

  sidx = sum(idx);
  s1 = find(idx);
  s2 = find(idx == 0);
  if sidx > 0,
    xx1 = xx0(s1,:);
    xx1r = xx1(:,1:p);
  end
  xx2r = xx0(s2,1:p);

  mm = mean(xx2r);
  Ss = cov(xx2r);
  if sidx > 0,
    aux1 = xx1r - ones(sidx,1)*mm;
    dd = Ss\aux1';
    dms = sum((aux1.*dd')');
  end
end

%% Check remaining observations and relabel if appropriate

while sidx > 0,

  s1 = find(dms <= ctf);
  s2 = length(s1);
  if s2 == 0,
    break
  end
  s3 = xx1(s1,pp);
  idx(s3) = zeros(s2,1);
  sidx = sum(idx);

  s1 = find(idx);
  s2 = find(idx == 0);
  if sidx > 0,
    xx1 = xx0(s1,:);
    xx1r = xx1(:,1:p);
  end
  xx2r = xx0(s2,1:p);

  mm = mean(xx2r);
  Ss = cov(xx2r);
  if sidx > 0,
    aux1 = xx1r - ones(sidx,1)*mm;
    dd = Ss\aux1';
    dms = sum((aux1.*dd')');
  end

end

% Values to be returned

idx = idx == 1;

aux1 = x - ones(n,1)*mm;
%Quito un poco del ruido que genera esta inversion de matriz
warning off
dd = Ss\aux1';
warning on
dms = sum((aux1.*dd')');

dm = sqrt(dms);