levmar(fcn, a0, ia, options, varargin)

function [a, chisq, ok, output] = levmar(fcn, a0, ia, options, varargin)
% mrqmin minimisation de levenberg marquardt 
%   [a, chisq, vals, ok, iter] = mrqmin(a0, fcn, ia, sig, varargin)
% la routine effectue l'optimisation, de sum((fcn(a, varargin)./sig).^2)
%
% paramtres d'entre
% -> fcn       fonction critre  minimiser vals = fcn(a, varargin)
%              cette fonction renvoie un vecteur colonne vals de  
%              dimension v. On cherche  minimiser 
%              sum(vals.^2)
% -> a0        valeur initiale du vecteur de paramtres  optimiser
% -> ia        logical les a(ia) doivent tre optimiss, pas les autres.
%              ( a == [] <=> a == ones(size(a)))
% -> varargin  donnes du problme passes  la fonction critre fcn
%
% paramtres de sortie
% <- a         solution nouvelle
% <- chisq     valeur du chi2  la solution 
% <- ok        0 (chec) ou 1 (succs) 
% <- output    structure : 
%    .iterations  nombre d'itrations effectues
%    .funcCount   nombre d'valuation de la fonction
%    .algorithm   nom de l'algorithme utilis

%gestion des options emprunte  fminsearch.
defaultopt = optimset('maxiter','200*numberOfVariables',...
   'TolCon',sqrt(eps),...
   'TolX',sqrt(eps));
% If just 'defaults' passed in, return the default options in a
if nargin==1 
   if nargout <= 1 & isequal(fcn,'defaults')
      a= defaultopt;
      return   
   else
      error('pas assez d''arguments pour levmar');
   end;
end;
if nargin < 3
   ia = ones(size(a0));
end;
if isempty(ia)
   ia = logical(ones(size(a0)));
else
   ia = logical(ia);
end;

if nargin<4, options = []; end
n = prod(size(a0(ia)));
numberOfVariables = n;

options = optimset(defaultopt,options);
tola = optimget(options,'TolCon');
tolr = optimget(options,'TolX');
maxiter = optimget(options,'maxiter');
% In case the defaults were gathered from calling: optimset('levmar'):
if ischar(maxiter)
   maxiter = eval(maxiter);
end

ok = 0; iter = 0; echelle = 0.001; fac = 10; dernum = 0;
try
   [vals dyda] = feval(fcn, a0, varargin{:});
   nb = 1;
catch
   vals = feval(fcn, a0, varargin{:});
   dyda = numjac(vals, a0, fcn, varargin{:});
   dernum=1;
   nb = prod(size(a0))+1;
end;
dyda(:,~ia)=[]; 
chisq = vals.'*vals;
ochisq = inf; atry = a0; a = a0; count = nb;
while 1
   iter = iter+1;
   if iter > maxiter
      ok = 0;
      break;
   end;
   m = dyda.'*dyda;
 
   m = m+diag(diag(m))*echelle;;
   da = -m\(dyda.'*vals);
   atry(ia) = a(ia)+da';
   if dernum
      vals = feval(fcn, atry, varargin{:});
      dyda = numjac(vals, atry, fcn, varargin{:});
   else
      [vals, dyda] = feval(fcn, atry, varargin{:}); 
   end;
   count = count+nb;
   dyda(:,~ia)=[];
   chisq = vals.'*vals;
   if chisq <= ochisq + eps
      dchi = abs(ochisq - chisq); 
      if (dchi < tola | dchi < tolr*chisq) 
         a = atry;
         ok = ok+1;
         if ok == 4
            ok = 1;
            break;
         end;
      end;
      ochisq = chisq;
      echelle = echelle/fac;
      a = atry;
   else
      echelle = fac*echelle;
      chisq = ochisq;
      ok = 0;
   end;
end;   
if nargout == 4
   output.iterations = iter;
   output.funcCount = count;
   output.algorithm = 'Levenberg-Marquardt';
   if dernum 
      output.algorithm = [ output.algorithm ' avec drives numriques'];
   else
      output.algorithm = [ output.algorithm ' avec drives analytiques'];
   end;
end;