function [Moyenne,ET,ICInf,ICSup, nmoda, indice_panel, VraisM, structHisto]= ...
    filtre_particulaire(fetat,hobs,paramobs,yobs,nombre_particules,tobs,...
    vobs,  vobs_or, deltaT, ...
    typebruitparam,pourcentbruitetat,fenetrenoyau,...
    fenetrepertubation, icbinf, icbsup, indice_panel, ...
    nbrep, ...
    nmod, nomparam, nomparamS, donnees, modelepop, ...
    affichage, nbhisto_tot,  nbhisto_deb,nbhisto_fin)
%filtre_particulaire Manage the particles filtering process using a convolution kernel
% (Mise en oeuvre du filtre particulaire avec convolution)
% INPUT
%       fetat(function) : Xt+1=fetat(Xt,t, etc)+ something
%       hobsbruit(function)  : Yt+1=hobs(Xt,bruit,paramobs)
%       paramobs : parameters realtive to the observations, size :(nbrepet,nbparam,nbobs)
%       yobs : observations variable
%       nombre_particules : number of particles
%       tobs : observations times, size: (nbobs)
%       vobs: computational times, size: (nt)
%       vobs_or: same in the genuine time unit
%       deltaT: computational time step
%       typebruitparam : pertubation noise type, vector of length >= number of parameters
%                        0 : no noise
%                        3 : binomial
%                        1 : poisson
%                        2 : gaussian
%       pourcentbruitetat : percentage of noise applied to the state
%       ('applique a l'etat') (inside [0,100])
%       fenetrenoyau : width of the kernel window  (>1)
%       fenetrepertubation : width of the pertubations window (>1)
%       icbinf, icbsup: IC bounds
%       affichage : 1 if histograms should be plotted as calculation goes,
%                   0 otherwise.
%       nbrep : number of the current dynamics ('repetitions')
%       nmod:  0 if one model, only
%              1 or 2: number of the current model (1 or 2) in case of comparison of two models
%       donnees: structure that contains the fixed parameters and the number of dynamics
%       nbhisto_tot, nbhisto_deb,nbhisto_fin: if affichage, total number of
%       possible histograms, number of histograms at the beginning, number
%       of histograms at the end.
%
% OUTPUT
%        Moyenne : mean of the states (etats) at all the times, size: (nrow(Xt), nt)
%        ET : standard deviation of the states (etats) at all the times, size: (nrow(Xt), nt)
%        ICInf ICSup : IC at level 90%,95% or 99% (according to icbinf and icbsup) of the states
%           at all the times, size: (nrow(Xt), nt)
%        VraisM : likelihood of the model (used in case of models
%        comparison only, i.e when nmod>0), size: (nbobs)
%        nmoda: histograms mod, whatever they are plotted or not, size: (1,nbparam)
%        structHisto: structure that contains all is necessary to save/restore the histograms
%======================================================

%% Initialisation %%
VraisM={};
structHisto={};
data=guidata(gcf);
nbobs = length(tobs); % nbre d'observations

% Initialisation des particules en utilisant les lois a priori
[Xt,Xmin,Xmax]=initialisation_particules(nombre_particules, nmod);
nbparam=size(Xt,1)-1;
Moyenne=zeros(size(Xt,1),length(vobs));
ET=Moyenne;
ICInf=Moyenne;
ICSup=Moyenne;
epsilonw=1e-40;

if nmod >0 % cas de la comparaison de modèles
    strnmod = num2str(nmod);
    VraisM=zeros(1,nbobs);
else
    strnmod='';
end


% pour l'ajout du bruit on cherche les indices correspondants aux
% differents types de bruit (3 types)
indicebruit=cell(1,3);
for ib=1:(2)
    indicebruit{ib}=find(typebruitparam==ib);
end
indicebruit{3}=find(typebruitparam>=3);

nbhisto_cur=0; % nbre courant d'histo (trac�s ou non)


%% Boucle sur le temps
for i=1:length(vobs)
    
    % On regarde si l'utilisateur n'a pas demand� une interruption
    % Quand on fait exécuter un fichier test, data.estimation_tag est un
    % cell array
    if ~iscell(get(data.estimation_tag,'userdata'))
        if ~get(data.estimation_tag,'userdata')
            return;
        end
    end
    ti=vobs(i);
    ti_or= vobs_or(i); % time in genuine time unit
    % Evolution des particules
    [Xt, dXt]=fetat(Xt,ti_or,pourcentbruitetat,  ...
        nomparam, nomparamS, donnees);
    
 % 21/01/2014: rajout du schema d'Euler
    Xt(end, :) = Xt(end, :) + deltaT .* dXt;
  
    % 5/12/2014: rajout du test isinf
    if any(isnan(Xt(end, :))) || any(isinf(Xt(end, :)))
        % Il y a eu une erreur de calcul dans le modele
        % on positionne le handle de la routine appelante
        h=gcbo;
        set(h, 'userdata', 0);
        errordlg(['An error occurred during model calculation at the observation ', ...
            num2str(i), ' at the time ', num2str(ti)]);
        Idle(1);
        Moyenne=[];
        return;
    end
    
    
    
     
    % Ramener a des valeurs valides
    Xt(end, Xt(end, :) < 0) = 0.0;
    
    % Calcul de Moyenne, ET, IC
    poids=1/size(Xt,2).*ones(1,size(Xt,2));
    Moyenne(:,i)=Xt(:,end);
    ET(:,i)=sqrt(var(Xt,poids,2));
    [TAB,Indice]=sort(Xt,2);
    
    

    
    for j=1:size(Xt,1)
        % On regarde si l'utilisateur n'a pas demand� une interruption
        if ~iscell(get(data.estimation_tag,'userdata'))
            if ~get(data.estimation_tag,'userdata')
                return
            end
        end
        
        w=poids(Indice(j,:));
        w=cumsum(w);
        II=find(w>=icbinf | w<= icbsup);
        if isempty(II)
            II=[1 n];
        end
        ICInf(j,i)=TAB(j,II(1));
        ICSup(j,i)=TAB(j,II(end));
    end
    refresh
    
    set(data.hPb,'Value',(i+length(vobs)*(nbrep-1))/(length(vobs)*data.nombre_dynamique)*100);
    refresh
    drawnow
    
    
    % Trouver l'observation qui correspond au temps de calcul
    
    indiceobs=find(tobs==ti);
    if ~isempty(indiceobs)
        % Un temps d'observation correspond au temps de calcul
        chaine=[' Observation #' num2str(indiceobs)] ;
        
        % calcul de l'observation theorique (nbre de répétitions xnombre de particules)
        Ytheorique=hobs(Xt,paramobs(:,:,indiceobs), modelepop);
        
        % On regarde si l'utilisateur n'a pas demand� une interruption
        if ~iscell(get(data.estimation_tag,'userdata'))
            if ~get(data.estimation_tag,'userdata')
                return
            end
        end
        
        %% Calcul des poids
        % differences entre les estimes (Ytheorique) et les observations (qxn)
        diff=(Ytheorique-repmat(yobs(:,indiceobs),1,nombre_particules));
        
        
        % Calcul de la vraisemblance du modèle dans le cas de la
        % comparaison de modèles
        if nmod >0
            VraisM(indiceobs) = calcul_vraisemblance( Ytheorique, diff,nombre_particules);
        end
        clear Ytheorique
        % ecart-type de ces differences (qx1)
        covdiff=std(diff,0,2);
        
        % Largeur de la fenetre du noyau (qx1)
        hn=covdiff./(nombre_particules^(1/fenetrenoyau));
        clear covdiff
        % on l'utilise pour normaliser diff
        diffnormalise=diag(1./hn)*diff;
        clear diff hn
        
        % Poids
        % par defaut, nlaplacien faux et en pratique
        % on ne laisse pas le choix a l'utilisateur
        %        nomvar= ['nlaplacien' strnmod];
        %         if strcmp(get(data.(nomvar), 'Checked'),'on')
        %             poids=exp(-0.5*sum(abs(diffnormalise),1));
        %         else
        poids=exp(-0.5*sum(diffnormalise.^2,1));
        %         end
        % normalisation des poids
        
        poids=poids/sum(poids);
        
        % nomvar=['redistribution' strnmod];
        
        % par defaut, redistribution est vrai et en pratique
        % on ne laisse pas le choix a l'utilisateur
        %        if strcmp(get(data.(nomvar), 'Checked'),'on')
        Neff=1/sum(poids.^2);
        if Neff<nombre_particules/(log10(nombre_particules)^i)
            poids=1/nombre_particules.*ones(1,nombre_particules);
        end
        %        end
        
        
        
        %% Echantillonnage sur les particules de poids > eps
        ind=find(poids>epsilonw);
        chaine_affiche=[chaine, ' => number of retained particles = ', ...
            num2str(length(ind))];
        % ecrire plutôt le no de la répétition si on a demandé une étude
        % de stabilité
        if data.nombre_dynamique >1
            chaine_affiche=[ 'Dynamics #' num2str(nbrep)];
        end
        
        % On regarde si l'utilisateur n'a pas demand� une interruption
        if ~iscell(get(data.estimation_tag,'userdata'))
            if ~get(data.estimation_tag,'userdata')
                return
            end
        end
        
        %% afficher le nbre de particules retenues ou le modèle courant
        switch nmod
            case 0
                set(data.com_tab_tag,'string',chaine_affiche)
            case 1
                set(data.com_tab_tag,'string','Numerator model identification in progress ... ')
            case 2
                set(data.com_tab_tag,'string', 'Denominator model identification in progress ... ');
        end
        
        %% Stockage des resultats au temps
        Moyenne(:,i)=Xt*(poids(:));
        ET(:,i)=sqrt(var(Xt,poids,2));
        %% Calcul des IC a 90% ou 95% ou 99%
        [TAB,Indice]=sort(Xt,2);
        
        for j=1:size(Xt,1)
            % On regarde si l'utilisateur n'a pas demand� une interruption
            if ~iscell(get(data.estimation_tag,'userdata'))
                if ~get(data.estimation_tag,'userdata')
                    return
                end
            end
            
            w=poids(Indice(j,:));
            w=cumsum(w);
            II=find(w>=icbinf | w<= icbsup);
            if isempty(II)
                II=[1 n];
            end
            ICInf(j,i)=TAB(j,II(1));
            ICSup(j,i)=TAB(j,II(end));
        end
        
        clear TAB
        
        
        %% Histogrammes
        if (affichage)
            
            nbhisto_cur = nbhisto_cur+1;
            
            if nbhisto_cur <= nbhisto_deb ||...
                    nbhisto_cur> (nbhisto_tot - nbhisto_fin)
                % sauvegarde des infos necessaires aux histos
                structHisto{end+1}= {indiceobs, Xt};
                
                
            end
            
        end %fin affichage
        
        
        %% Calcul des mod
        for j=1:nbparam
            [nm, nxout] = hist(Xt(j,:), 20);
            indicem = find(nm==max(nm));
            nmoda(1,j) = nxout(indicem(1)); % On prend l'abscisse du 1ier max
        end
        
        %% Tirage multinomial
        ind=randsample(ind,nombre_particules,true,poids(ind));
        Xt=Xt(:,ind);
        Xtempo=Xt; % sauvegarde des particules
        
        %% Rajout d'un bruit apres re-echantillonage (etape de regularisation)
        % Calcul de l'ecart type  de la perturbation
        ETPerturbation=std(Xtempo,0,2);
        ETPerturbation=ETPerturbation/(nombre_particules^(1/fenetrepertubation));
        particulesachanger=1:nombre_particules;
        % pour toutes les particules
        % tant qu'elles ne sont pas dans [Xmin,Xmax]
        while(~isempty(particulesachanger))
            % On regarde si l'utilisateur n'a pas demand� une interruption
            if ~iscell(get(data.estimation_tag,'userdata'))
                if ~get(data.estimation_tag,'userdata')
                    return
                end
            end
            
            % variance dependante du nombre de particule et de l'etat
            % Diversification des particules par regularisation
            % poisson
            indp=1; % indice de Poisson dans la liste des types de bruits
            if (~isempty(indicebruit{indp}))
                ETPert=ETPerturbation(indicebruit{indp}).^2;
                ETPert=repmat(ETPert,1,length(particulesachanger));
                bruit=poissrnd(ETPert)-round(ETPert);
                Xt(indicebruit{indp},particulesachanger)=Xtempo(indicebruit{indp},particulesachanger)+bruit;
                clear ETPert bruit
            end
            % gaussien
            indp=2; % indice de gaussien dans la liste des types de bruits
            if (~isempty(indicebruit{indp}))
                bruit = diag(ETPerturbation(indicebruit{indp}))*randn(length(indicebruit{indp}),length(particulesachanger));
                Xt(indicebruit{indp},particulesachanger)=Xtempo(indicebruit{indp},particulesachanger)+bruit;
                clear bruit
            end
            % binomial
            indp=3; % indice de binomial  dans la liste des types de bruits
            if (~isempty(indicebruit{indp}))
                bruit=round(diag(ETPerturbation(indicebruit{indp}))*randn(length(indicebruit{indp}),length(particulesachanger)));
                Xt(indicebruit{indp},particulesachanger)=Xtempo(indicebruit{indp},particulesachanger)+bruit;
                clear bruit
            end
            % on cherche les particules qui ne sont pas dans le domaine de recherche
            particulesachanger=find(any(Xt(1:nbparam,:)<Xmin,1)+any(Xmax<Xt(1:nbparam,:),1));
        end
    end % fin un temps d'observation=le temps de calcul
    
end % fin i (boucle sur vobs)



handles=guihandles(1);
ndata=catstruct(data,handles);
guidata(1,ndata);

end %% fin filtre_particulaire
%% ==============================================================
