function [ evoldet, dx , proba, iter, msgerr]=tornsey(der, reseau, epscal, maxiter , ...
						      nsteps, istep, hbar, obj, evnt)
%tornsey Algorithme de Tornsey
% INPUT
%   der: matrice des derivées du modèle (nombre de points du réseau, nombre
%     de parametres du modèle)
%   reseau: les points de calcul (des temps)
%   nsteps: nombre d'étapes de calcul
%   istep: indice de l'étape avant exécution
%   hbar: barre de progression
%   obj: handle de l'appelant
%    evnt: 1 quand on fait passer des tests
%OUTPUT
% evoldet: matrice d'autant de maxiter lignes et de 2 colonnes
% Seules les lignes de 1:nombre d'iterations effectuées sont utiles.
% 1iere colonne: les déterminants des matrices d'information
% 2ieme colonne: le log de ces determinants ou 0 s'ils ne sont pas >0
% dx: les valeurs de la courbe a tracer (autant que de points de calcul)
% proba: les proba (autant que de points du réseau)
% iter: indice de la dernière itération valide
% msgerr: message d'erreur ou ''
%% =========================================================================
iter=0;
msgerr=''; % message d'erreur
nx=length(reseau); % nombre de temps de calcul
px=size(der,2); % nombre de parametres
proba=zeros(1,nx);
proba(:)=1/nx;
evoldet=zeros(maxiter,2);
dx=zeros(1,nx);
dummy=1/(nx+(0.05*nx));
if epscal > dummy
    epscal=dummy;
end

Minf = zeros(px,px); %matrice d'information
% NOTE l'inversion de Minf ne marche pas toujours (tmax=504) même avec
% l'utilisation de mldivide

%% Boucle des itérations
for iter=1:maxiter
    istep=istep+1; % l'étape de calcul pour la barre de progression
    waitbar(istep/nsteps,hbar, ['Run Tornsey algorithm, iteration ', ...
        num2str(iter)])
    % Check for Cancel button press
    if isempty(evnt)
    if getappdata(hbar,'canceling') ||  get(obj,'userdata')==0
        iter=iter-1;
        return
    end
    end
    
    for i=1:nx
        if proba(i)>0
            Minf=Minf+ transpose(der(i,:))*der(i, :)*proba(i);
        end
    end
    deter=det(Minf); % determinant
    evoldet(iter,1)=deter;
    if deter >0
        evoldet(iter,2)=log(deter);
    end
    
    % Calcul de l'inverse par utilisation de svd car
    % il arrive que la matrice à inverser soit singulière et alors inv ne marche pas
    if any(any(isnan(Minf)))
        iter=iter-1; % l'iteration en cours donne de mauvais resultats
        msgerr=['Information matrix non inversible. Out at iteration ', num2str(iter)];
        
        break; % sortie boucle sur les iterations
    end
    invm=svdInverse(Minf);
    
    if any(any(isnan(invm)))
        iter=iter-1; % l'iteration en cours donne de mauvais resultats
        msgerr=['Information matrix non inversible. Out at iteration ', num2str(iter)];
        
        break; % sortie boucle sur les iterations
    end
    
    for i=1:nx
        if proba(i)>0
            dx(i)=der(i,:)*invm*transpose(der(i,:));
            % On peut faire aussi, quand ça marche:      dx(i)=der(i,:)* (Minf\transpose(der(i,:)));
            % On peut faire aussi, quand ça marche:   dx(i)=(der(i,:)/Minf)*transpose(der(i,:));
            proba(i)=proba(i)*dx(i)/px;
        end
    end
    
    proba(proba < epscal)=0;
    if ~any(proba >0 & abs(dx-px)>epscal)
        iter=iter-1; % l'iteration en cours donne de mauvais resultats
        %         msgerr=['Sortie anticipée à l''itération ', num2str(iter)];
        msgerr=['Null proba or EPSILON not reached. Out at iteration ', num2str(iter)];
        break; % sortie boucle sur les iterations
    end
    Minf(:,:)=0;
    
end % fin boucle sur iter

end % fin tornsey

    
        
            
            
            
    
    

