function [Weight, ProjectMatrix] = LocalFeatureExtraction(patterns, targets, KN, method)

%Final Code: Apr 10, 2006
%===============================================================
%FUNCTION: FeatureWeighting implementing the local feature extraction (LFE)
%algorithms, (diagonal matrix and full matrix). 
%Currently, it only supporting the binary problem. The extension to
%multi-class problem is easy.
%
%A New Interpretation of RELIEF and its Extension to Feature Extarction, NIPS paper and Pattern Recognition Paper
%
%==========================================================================
%INPUT:
%     patterns:  training data: Each column is an observation: X =[x1,x2...]
%     targets:   class label (-1,+1)
%     KN :       # of nearest neighbor
%     method:    knn: search for kNN; or weight: search for KNN and weight based on distance; 
%OUTPUT:
%    Weight: weight of features
%    ProjectMatrix: project matrix
%--------------------------------------------------------------------------
%University of Florida, by Yijun Sun
%LUD: September 1, 2005/Sep 25, 2005/ Feb 9/ Apr 10, 2006
%==========================================================================
N = length(targets);    % # of patterns
dim = size(patterns,1); % Data dimenionality

index0 = find(targets==-1);
index1 = find(targets==1);
N0 = length(index0); % # of patterns in class 1
N1 = N-N0;           % # of patterns in class 2

patterns0 = patterns(:,index0); %patterns in class 1
patterns1 = patterns(:,index1); %patterns in class 2
%--------------------------------------------------------------------------

NM = zeros(dim,KN); 
NH = zeros(dim,KN); 
Sigma_NM_NH = zeros(dim,dim);

switch method
    
    case {'weight'}
        for i = 1:N,
            Temp0             = (patterns0 - patterns(:,i)*ones(1,N0)).^2;
            dist0             = sum(Temp0,1);
            
            Temp1             = (patterns1 - patterns(:,i)*ones(1,N1)).^2;
            dist1             = sum(Temp1,1);
            
            [m0, indices0]    = sort(dist0);
            [m1, indices1]    = sort(dist1);
            
            if targets(i)==-1; 
                weight0 = (1-(m0(2:1+KN)/max(m0(2:1+KN))).^(1.5)).^3;
                weight0 = weight0/sum(weight0);
                weight1 = (1-(m1(1:KN)/max(m1(1:KN))).^(1.5)).^3;
                weight1 = weight1/sum(weight1);
                
                for k=1:KN
                    NH(:,k) = sqrt(weight0(k))*(patterns(:,i)-patterns0(:,indices0(1+k)));
                    NM(:,k) = sqrt(weight1(k))*(patterns(:,i)-patterns1(:,indices1(k)));
                end
            end
            
            if targets(i)==1;
                weight0 = (1-(m0(1:KN)/max(m0(1:KN))).^(1.5)).^3;
                weight0 = weight0/sum(weight0);
                weight1 = (1-(m1(2:KN+1)/max(m1(2:KN+1))).^(1.5)).^3;
                weight1 = weight1/sum(weight1);
                
                for k=1:KN
                    NM(:,k) = sqrt(weight0(k))*(patterns(:,i)-patterns0(:,indices0(k)));
                    NH(:,k) = sqrt(weight1(k))*(patterns(:,i)-patterns1(:,indices1(1+k)));
                end
            end 
            Sigma_NM_NH = Sigma_NM_NH + (NM*(NM')-NH*(NH'));
        end
        
    case{'knn'}
        
        for i = 1:N,
            Temp0             = (patterns0 - patterns(:,i)*ones(1,N0)).^2;
            dist0             = sum(Temp0,1);
            
            Temp1             = (patterns1 - patterns(:,i)*ones(1,N1)).^2;
            dist1             = sum(Temp1,1);
            
            [m0, indices0]    = sort(dist0);
            [m1, indices1]    = sort(dist1);

            if targets(i)==-1; 
                for k=1:KN
                    NH(:,k) = (patterns(:,i)-patterns0(:,indices0(1+k)));
                    NM(:,k) = (patterns(:,i)-patterns1(:,indices1(k)));
                end
            end
            
            if targets(i)==1;
                for k=1:KN
                    NM(:,k) =(patterns(:,i)-patterns0(:,indices0(k)));
                    NH(:,k) =(patterns(:,i)-patterns1(:,indices1(1+k)));
                end
            end 
            Sigma_NM_NH = Sigma_NM_NH + (NM*(NM')-NH*(NH'));
        end 
end

[V,D] = eig(Sigma_NM_NH);

index =find(diag(D)>0);
Weight = 0;
for n = 1:length(index);
    Weight = Weight+D(index(n),index(n))*V(:,index(n))*(V(:,index(n))');
end

[D,index] =sort(diag(D));
D = flipud(D); index = flipud(index);

D(D<0)=0;
ProjectMatrix = zeros(dim,dim);
for n=1:dim
    ProjectMatrix(n,:) = sqrt(D(n))*(V(:,index(n))');%
end

return

%=====================================END OF CODE =========================
%==========================================================================






%case{'ArFace_knn'}
    
    for i = 1:N,
        for n= 1:N0
            dist0(n)             = sum((patterns0(:,n) - patterns(:,i)).^2);
        end
        
        for n=1:N1
            n
            dist1(n)             = sum((patterns1(:,n) - patterns(:,i)).^2);
        end            
        [m0, indices0]    = sort(dist0);
        [m1, indices1]    = sort(dist1);
        
        keyboard
        if targets(i)==-1; 
            for k=1:KN
                NH(:,k) = (patterns(:,i)-patterns0(:,indices0(1+k)));
                NM(:,k) = (patterns(:,i)-patterns1(:,indices1(k)));
            end
        end
        
        if targets(i)==1;
            for k=1:KN
                NM(:,k) =(patterns(:,i)-patterns0(:,indices0(k)));
                NH(:,k) =(patterns(:,i)-patterns1(:,indices1(1+k)));
            end
        end 
        Sigma_NM_NH = Sigma_NM_NH + NM*(NM');
        Sigma_NM_NH = Sigma_NM_NH - NH*(NH');%/(N);
    end 


%     figure;imagesc(ProjectMatrix'*ProjectMatrix); 
%     title('Reconstructed Distance Matrix')

% ratio =0;
% DiscriminantPower = sum(D(D>0));
% N=0;
% while ratio<R
%     N=N+1;
%     ratio = sum(D(1:N))/DiscriminantPower;
% end

% S = 0;
% for n=1:N;
%     S=S+D(n)*V(:,index(n))*(V(:,index(n))');
% end



















%==========Feature extraction from the principal components=============
[D,index] =sort(diag(D));
D = flipud(D); index = flipud(index);
ratio =0;
DiscriminantPower = sum(D(D>0));
N=0;
while ratio<0.8
    N=N+1;
    ratio = sum(D(1:N))/DiscriminantPower;
end

S = 0;
for n=1:N;
    S=S+D(n)*V(:,index(n))*(V(:,index(n))');
end
figure;imagesc(S); title('Reconstructed Distance Matrix from PC')
s2 = diag(S);

Ratio = sum(D(1:N))/sum(D(D>0));
figure; plot(s1,'-o');
hold on;plot(s2,'-r*');legend('1','2')
title(['diagonal terms: Ration =' num2str(Ratio*100) ' %'])
keyboard


return


%===============================================================
%FUNCTION: FeatureWeighting implementing the feature weighting algorithms
%(diagonal matrix and full matrix).
%INPUT:
%     patterns:  training data: Each column is an observation
%     targets:   class label (-1,+1)
%    test_patterns:
%    test_targets:
%     KN :  # of nearest neighbor
%     it: # of iteration
%   distance: 'Euclidean' or 'Block'
%OUTPUT:
%    Weight: weight of features
%--------------------------------------------------------------------------
%University of Florida, by Yijun Sun
%September 1, 2005
%==========================================================================
N = length(targets);% # of patterns
dim = size(patterns,1); % Data dimenionality

index0 = find(targets==-1);
index1 = find(targets==1);
N0 = length(index0);
N1 = N-N0;

patterns0 = patterns(:,index0);
patterns1 = patterns(:,index1);
ProjectMatrix = eye(dim,dim);
Old_Weight = eye(dim,dim);
%--------------------------------------------------------------------------
for t=1:it
    
    NM = zeros(dim,N,KN); 
    NH = zeros(dim,N,KN); 
    Sigma_NM = zeros(dim,dim);
    Sigma_NH = zeros(dim,dim);
    
    New_patterns = ProjectMatrix*patterns;
    New_patterns0 = New_patterns(:,index0);
    New_patterns1 = New_patterns(:,index1);
    
    for i = 1:N,
        %         Temp0             = (patterns0 - patterns(:,i)*ones(1,N0)).^2;
        %         dist0             = sum(Temp0,1);
        %         
        %         Temp1             = (patterns1 - patterns(:,i)*ones(1,N1)).^2;
        %         dist1             = sum(Temp1,1);
        
        Temp0             = (New_patterns0(1:2,:) - New_patterns(1:2,i)*ones(1,N0)).^2;
        dist0             = sum(Temp0,1);
        
        Temp1             = (New_patterns1(1:2,:) - New_patterns(1:2,i)*ones(1,N1)).^2;
        dist1             = sum(Temp1,1);
        
        [m0, indices0]    = sort(dist0);
        [m1, indices1]    = sort(dist1);
        
        if targets(i)==-1; 
            for k=1:KN
                NH(:,i,k) = (patterns(:,i)-patterns0(:,indices0(1+k)));
                NM(:,i,k) = (patterns(:,i)-patterns1(:,indices1(k)));
            end
        else
            for k=1:KN
                NM(:,i,k) = (patterns(:,i)-patterns0(:,indices0(k)));
                NH(:,i,k) = (patterns(:,i)-patterns1(:,indices1(1+k)));
            end
        end    
    end
    
    for k=1:KN
        Sigma_NM = Sigma_NM + NM(:,:,k)*(NM(:,:,k)')/(N*KN);
        Sigma_NH = Sigma_NH + NH(:,:,k)*(NH(:,:,k)')/(N*KN);
    end
    figure;imagesc(Sigma_NM-Sigma_NH);title('Sigma(m-h)')
    [V,D] = eig((Sigma_NM-Sigma_NH));
    
    index =find(diag(D)>0);
    Weight = 0;
    for n = 1:length(index);
        Weight = Weight+D(index(n),index(n))*V(:,index(n))*(V(:,index(n))');
    end
    diff(t) = norm(Weight-Old_Weight,'fro');
    Old_Weight =Weight;
    
    %     figure;imagesc(Weight); 
    %     title('Reconstructed Distance Matrix')
    w(:,t) = (diag(Weight)); 
    
    %      D = inv(D);
    %      D(D<0)=0;
    %      ProjectMatrix = sqrt(D)*(V');
    
    Tem = diag(Sigma_NM-Sigma_NH); Tem(Tem<0)=0;
    ProjectMatrix = sqrt(diag(Tem)); keyboard
end

% S = (ProjectMatrix')*ProjectMatrix;
% figure;imagesc(S)
figure;
plot(diff,'-o')

figure;
for n=1:10
    subplot(2,5,n)
    plot(w(:,n),'-o')
    title(num2str(n))
end
keyboard