function [Moyenne,ET,ICInf,ICSup, nmoda, indice_panel]= ...
    filtre_particulaire_sivip(fetat,hobs,paramobs,yobs,nombre_particules,tobs,...
    vobs, vobs_or,  ...
    typebruitparam,pourcentbruitetat,fenetrenoyau,...
    fenetrepertubation, icbinf, icbsup, indice_panel, ...
    nomparam, nomparamS, paramS, modelepop, sivipargu, evnt)
%filtre_particulaire Filtre particulaire avec convolution
% Mise en oeuvre du filtre particulaire avec convolution
% INPUT
%       fetat(fonction) : Xt+1=fetat(Xt,t, etc)
%       hobsbruit(fonction)  : Yt+1=hobs(Xt,bruit,paramobs)
%       paramobs : parametres lies aux observations, dim :(nbrepet,nbparam,nbobs)
%       yobs : variable d'observation
%       nombre_particules : nb de particules
%       tobs :  temps d'observations (nbobs)
%       vobs: temps de calcul (nt)
%       vobs_or: same in the genuine time unit 
%       typebruitparam : type de bruit de la pertubation (vecteur :
%                        nblignes de Xt=nbre de parametres)
%                        0 : pas de bruit
%                        3 : binomial
%                        1 : poisson
%                        2 : gaussien
%       pourcentbruitetat : pourcentage de bruit applique a l'etat ds [0,100]
%       fenetrenoyau : largeur de la fenetre du noyau >1
%       fenetrepertubation : largeur de la fenetre des pertubations >1
%        icbinf, icbsup: bornes de l'IC
%       paramS: structure contenant les parametres fixes
%       modelepop
%       sivipargu, les arguments pour le calcul des sivip
%        evnt=1 si passage de test, vide sinon
%
% OUTPUT
%        Moyenne : moyenne des etats a tout les instants (nbligne= nblignes de Xt, nbcol= nt)
%        ET : ecarts types des etats a tout les instants (nbligne=nblignes de Xt, nbcol= nt)
%        ICInf ICSup : I-C a 90%,95% ou 99% (selon icbinf, icbsup) des etats a tout les instants (nbligne=nblignes de Xt, nbcol= nt)
%        nmoda: les mod des parametres (1, nparam)
%======================================================

%% Initialisation %%

data=guidata(gcf);
% pour stocker ce qui est necessaire aux graphes pour restaurer
data.PIST={}; % les IST en % pour chaque temps ou ils sont calculés
data.Q2cumY1={};% les Q2cum pour chaque temps ou ils sont calculés
data.explY1c={}; % taux d'explication de Y, pour chaque temps ou ils sont calculés
data.tpsopt=[]; % temps du prochain temps optimal (un pour chaque graphe)
data.yretenu=[]; % valeur du prochain temps optimal (un pour chaque graphe)
data.iretenu=[]; % indice du parametre correspondant au  prochain temps optimal

nmod=0; % un seul modele ici
%% 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;


%% 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);

tist=0; %  le 1ier temps optima , initialisé au temps de la 1iere observation
tistcalcul=0; % le temps où calculer les prochains IST
pas_temps=vobs_or(2)-vobs_or(1); % pas de temps entre 2 temps de calcul

%% 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, evnt est non null
    if isempty(evnt)
        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, paramS);
    

    % 21/01/2014: rajout du schema d'Euler
    Xt(end, :) = Xt(end, :) + pas_temps .* 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 in the model calculation at observation ', ...
            num2str(i), ' and time ', num2str(ti)]);
        Idle(1);
        Moyenne=[];
        return;
    end

    % Ramener a des valeurs valides
    Xt(end, Xt(end, :) < 0) = 0.0;
    
    %% Calcul des SIVIP
    if ti_or==tistcalcul(end)
        %% Calcul des SIVIP
 
        [tist, tistcalcul, IST, PIST, Q2cumY1, explY1c, indice_panel, data ] =calculer_sivip(i, ti_or, tist, tistcalcul, ...
            fetat, Xt, vobs_or, sivipargu,     indice_panel, ...
            pourcentbruitetat,  nombre_particules, ...
            nomparam, nomparamS, paramS, data, evnt);
        data.PIST{end+1}=PIST;
        data.Q2cumY1{end+1}=Q2cumY1;
        data.explY1c{end+1}=explY1c;
 
 
    end
    
    % On regarde si l'utilisateur n'a pas demand� une interruption
    if isempty(evnt)
        if ~get(data.estimation_tag,'userdata')
            return
        end
    end
    
    % 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 isempty(evnt)
            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
    %  Maj et affichage barre de progression
    set(data.hPb,'Value',((i+3)/length(vobs))*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);
        
        %% Calcul des poids
        % differences entre les estimes (Ytheorique) et les observations (qxn)
        diff=(Ytheorique-repmat(yobs(:,indiceobs),1,nombre_particules));
        
        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
        %     if strcmp(get(data.nlaplacien, '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);
        
        % par defaut, redistribution est vrai et en pratique
        % on ne laisse pas le choix a l'utilisateur
        %    if strcmp(get(data.redistribution, '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))];
        
        % afficher le nbre de particules retenues
        set(data.com_tab_tag,'string',chaine_affiche)
        
        %% Stockage des resultats au temps
        Moyenne(:,i)=Xt*(poids(:));
        ET(:,i)=sqrt(var(Xt,poids,2));
        %% Calcul des IC a 95%
        [TAB,Indice]=sort(Xt,2);
        
        for j=1:size(Xt,1)
            % On regarde si l'utilisateur n'a pas demand� une interruption
            if isempty(evnt)
                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
        
        %% 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 isempty(evnt)
                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
       
     
    
    
    %%  Sauvegarde des temps optimaux et de calcul
    data.optima=tist;
    data.tistcalcul=tistcalcul;
    handles=guihandles(1);
    ndata=catstruct(data,handles);
    guidata(1,ndata);
    
end %% fin filtre_particulaire_sivip
%% ==============================================================

