function [DRy elite] = DR(P,T,PT,yy,classifier)

if nargin < 5
    classifier = 'NNKlite';
end

[q, Q, NN, class] = nc(T); % Class information
[N na] = size(P);
n = size(PT,1);

% Create Validation set (V matrix and TV vector)

[ITR ITE TTR TTE] = KFCV(P,T,4);

th(4) = 0;
for k = 1:4
    
    P1 = ITR{k}; T1 = TTR{k};
    V = ITE{k}; TV = TTE{k};
    N1 = length(T1); m = length(TV);
    for i = 1:q
        temp = (1:m)'.*(TV==Q(i));
        classV{i} = temp(temp>0);
        LC(i) = length(classV{i});
    end

    DRY(m) = 0;
    [y DCy] = feval(classifier,P1,T1,P1);
    DCY = sum(DCy);
    AR = sum(y==T1'); % Self-assessment of classifier
    for i = 1:m
        X = V(i,:);
        CX = (yy(i)==Q);
        P2 = [P1; X];
        T2 = [T1; Q(CX)];
        [y DCy] = feval(classifier,P2,T2,P1);
        d = sum((rm(X,(N-m),1) - P1).^2,2);
        classi = classV{CX};
        otherclasses = 1:N1;
        otherclasses(classi) = [];
        z1 = sqrt(LC(CX) / (sum(1./(d(classi) + eps)) - eps));
        z2 = sqrt(( N1 - LC(CX) ) / sum(1./(d(otherclasses) + eps) - eps));
        Z = max( ((z2 - z1) / (z1 + z2)) , 0);
        DRY(i) = ((sum(y==T1')/AR)*sum(DCy)/DCY * Z).^(0.5);
    end


    TH = 1.05; M = 0;
    for i = 1:20
        TH = TH - 0.05;
        temp = (1:m).*(DRY>=TH);
        ELITE = temp(temp>0);
        if isempty(ELITE)
            Z(i) = NaN;
        else
            Z(i) = (sum(y(ELITE)==TV(ELITE)')) / length(ELITE);
            if Z(i) > M
                M = Z(i);
                ind2 = i;
            end
        end
    end

    [M ind1] = max(Z);
    th(k) = 1.05 - 0.025*(ind1 + ind2);
end % of k loop
    

th = median(th)
DRY(1:n) = 0;
[y DCy] = feval(classifier,P,T,P);
DCY = sum(DCy);
AR = sum(y==T'); % Self-assessment of classifier
for i = 1:n
    X = PT(i,:);
    CX = (yy(i)==Q);
    P2 = [P; X];
    T2 = [T; Q(CX)];
    [y DCy] = feval(classifier,P2,T2,P);
    d = sum((rm(X,N,1) - P).^2,2);
    classi = class{CX};
    otherclasses = 1:N;
    otherclasses(classi) = [];
    z1 = sqrt(NN(CX) / (sum(1./(d(classi) + eps)) - eps));
    z2 = sqrt(( N - NN(CX) ) / sum(1./(d(otherclasses) + eps) - eps));
    Z = max( ((z2 - z1) / (z1 + z2)) , 0);
    DRY(i) = ((sum(y==T')/AR)*sum(DCy)/DCY * Z).^(0.5);
end

f = (DRY >= th);
temp = (1:n).*f;
elite = temp(temp>0);
DRy = min(DRY,1);