function ris = training_ljmNL3 (cl, indx, N_GAUSS, comp, scale, tipo)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% ris = training_ljmNL3 (cl, indx, N_GAUSS, comp, scale, tipo)
%
% Date: 29 Dec 2003
% Description: training using averaged local jet features with third non
%   linearity. For each pixel a set of invariants are extracted and 
%   each of these invariants is splitted in positive negative and absolute
%   value (with respect to the background mean). The averaged values over
%   all image pixel are the feature vector. After dimension reduction
%   obtained thanks to FLD a MoG is estimated for all the data of each
%   class.
% USE:  cl classes to consider, ex. cl=1:12;
%       indx indexes of images to consider (the
%           same for all classes), i.e. it=1; indx{1}=0:499; etc.
%       N_GAUSS nr. of gaussians in MoG
%       comp nr. of dimensions in the final feature space
%       scale nr. of scale in wich invariants are computed; it is assumed
%           scale=3;
%       tipo define covariance matrix structure in MoG, ex. tipo = 'full'
%       ris variable storing the training models
% Author: Marc'Aurelio Ranzato
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf(1,'TRAINING \n')
    
addpath('netlab/')
n_tot_samp = 0;
for cc = 1:length(cl)
    n_samples(cc) = length(indx{cc});
    n_tot_samp = n_tot_samp + n_samples(cc);
end
    
storet = zeros(n_tot_samp,(scale + 1)*9*3);

MAX_N_ITER = 1000;  % see gmmem
% Initial condition for randn and rand
%%%load('performance_SEPT/BestInitCondSEPTSim.mat')
%%%s = stateo;

lc = length(cl);

index2 = 0;

for cc = 1:lc   % Load data
    %keyboard
    load(['DB/features/data_collection_avljNL3_cl' int2str(cl(cc)) '.mat'])
    storet(index2 + 1: index2 + n_samples(cc),:) = store(indx{cc}+1,:);  
    load(['DB/features/resol_collection_avljNL3_cl' int2str(cl(cc)) '.mat'])
    resolution(cc,1) = mean(resol(indx{cc}+1));
    resolution(cc,2) = std(resol(indx{cc}+1));
    
    index2 = index2 + n_samples( cc );
    
end
                

%%% F L D %%%
fprintf(1,'Computing FLD \n')
mediatot = mean(storet);
storet = storet - repmat(mediatot,size(storet,1),1);
FLD = fisher_l_d(storet',comp,n_samples);
storet = (FLD'*storet')';   % on the 1st row the projection along FLD of 1st image features

fprintf(1,'Computing MoG \n');
% stima MoG
index2 = 0;
for jj = 1:lc
    
    st = storet(index2 + 1: index2 + n_samples(jj),:);
    
    fprintf(1,'computing MoG-%u \n', cl(jj));
    randn('state',0);
    %%randn('state',sum(100*clock))
    %%s{1}{jj} = randn('state');
    %%%randn('state',s{1}{jj});
    mix = gmm(size(st,2), N_GAUSS{jj}, tipo{jj});
    option = zeros(1,14);
    option(1) = -1;
    option(3) = 1;
    option(5) = 1;
    option(14) = MAX_N_ITER;
    rand('state',0);
    randn('state',0);
    %%rand('state',sum(100*clock))
    %%s{2}{jj} = rand('state');
    %%randn('state',sum(100*clock))
    %%s{3}{jj} = randn('state');
    %%%rand('state',s{2}{jj});
    %%%randn('state',s{3}{jj});
    mix = gmminit(mix, st, option);
    mix = gmmem(mix, st, option);
    mog(jj).mix = mix;
    index2 = index2 + n_samples(jj);
    
end
    
%%ris.state =  s;
ris.mog = mog;
ris.FLD = FLD;
ris.mediatot = mediatot;
ris.resolution = resolution;

%save(['Exp1avLJNL3_s' int2str(scale) '_ng' int2str(N_GAUSS)  tipo '_d' int2str(comp) '_tots' int2str(n_tot_samp)],'ris');
rmpath('netlab/')
