function [CM, ER, ERp] = test_ljmNL3 (cltr, clts, indx, ris, scale, Thres)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% [CM,ER] = test_ljmNL3 (cl, indx, ris, scale)
% Date 29 Dec. 2003
% Description: test using averaged local jet features with 3rd non
%   linearity. The classification is done taking the class with maximum
%   likelihood (hypothesis of equally probable classes)
% Use:  cl classes to consider, for ex. cl=1:12
%       indx data structure with images to be tested, for ex.
%           indx{1}=500:545; etc.
%       ris models of training
%       scale nr. of scale in wich invariants are computed
%       CM confusion matrix
%       ER error rate
% Author: Marc'Aurelio Ranzato  
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

DEBUG = 1;
addpath('netlab/')
pat = {'BACT' 'RBC' 'SPRM'};
% resolution statistics
resolution = ris.resolution;
NDS = 10;            % nr. of standard deviations 
% that is, don't do the test on resolution

ncl = length(clts);
%Thres = .9;
ncltr = length(cltr);

car = {'/' '-' '\' '|'};    

CM = zeros(ncl,ncl+1);
mog = ris.mog;
mediatot = ris.mediatot;
FLD = ris.FLD;


% MAIN LOOP: OUT -> class, IN -> index
fprintf(1,'TEST \n')
for jj=1:ncl

    fprintf(1,'cl %s: ',pat{clts(jj)})
    load(['DB/features/data_collection_avljNL3_cl' int2str(clts(jj)) '.mat'])
    load(['DB/features/resol_collection_avljNL3_cl' int2str(clts(jj)) '.mat'])
    
    for ii=1:length(indx{jj})
        
        fprintf(1,'\b');
        fprintf(1,'%c',car{mod(ii,4)+1});
        
        val = store(indx{jj}(ii)+1,:);
        risol = resol(indx{jj}(ii)+1);
        
        val = val - mediatot;
        val = (FLD' * val')';
        
        for jjj= 1:ncltr
            
            dprobab(jjj) = gmmprob( mog(jjj).mix, val );
            
        end
        
        if sum(dprobab) == 0
            Mp = 0;
            estimate = 1;
        else
            prob = dprobab/sum(dprobab);
            Mp = max(prob);
            estimate = find(prob == Mp);
            estimate = estimate(1);
        end
        
        if (risol >= resolution(estimate,1) + NDS*resolution(estimate,2) | risol <= resolution(estimate,1)-NDS*resolution(estimate,2)) & Mp <= Thres 
            
            CM(jj,1) = CM(jj,1) + 1;
            
        else
            
            CM(jj,estimate+1) = CM(jj,estimate+1) +1;   % update CM
            
        end
            
        
    end
    
    if DEBUG == 1
        if ncl ~= ncltr & jj == ncl
             fprintf(1,' %u on %u', CM(jj,1), length(indx{jj}) );
        else
            fprintf(1,' %u on %u', CM(jj,jj+1), length(indx{jj}) );
       end
    end

    fprintf(1,' \n');
  
end

% Computation of ER
Tot = sum(sum(CM));
ER = 1 - sum(diag(CM(:,2:end)))/Tot;

ERp = 1;
if ncl == ncltr
    
    for iii = 1 : size(CM,1)
        NN(iii) = sum(CM(iii,:));
        ERp = ERp - (CM(iii,iii+1)/NN(iii))/size(CM,1);
    end
    %%%ERp = ERp - CM(ncl,1)/(sum(CM(ncl,:))*ncl); % unknown class is OK for FA class
else
    for iii = 1 : size(CM,1)-1
        NN(iii) = sum(CM(iii,:));
        ERp = ERp - (CM(iii,iii+1)/NN(iii))/size(CM,1);
    end
    ERp = ERp - CM(ncl,1)/(sum(CM(ncl,:))*ncl);
    
end

rmpath('netlab/')