/*
   Several people have asked about non-linear fitting routines.
   Here is my own rough cut at such a routine.  See also lmfit.i
   by Eric Thiebaut.  Obviously I have some reservations about
   whether mine actually works; if it really is broken and you
   fix it, let me know.  I thought it worthwhile to include both
   my very stripped down interface and Eric's very complicated one
   here; probably something intermediate (perhaps with hook functions
   to optionally provide some of the more recondite bells and whistles)
   would be optimal.
   -- Dave Munro
 */

func levmar(x, y, y_and_dyda, a, weight=)
/* DOCUMENT a= levmar(y, x, y_and_dyda, a0)

***WARNING****
I think this function is broken.  It works in the au_refl_coeffs
routine because it bails out after zero or one iterations, basically
leaving the (very good) initial guess unchanged.
***WARNING****

     Perform a non-linear least squares fit to data vectors (X,Y)
     using the Levenberg-Marquardt method.
     Y_AND_DYDA must be defined as:
          func Y_AND_DYDA (a, x, &dyda)
     returning both the function value y(x) and derivative dyda(x).
     The input vector a is the list of parameters to be minimized;
     A0 is the initial guess for these parameters.  If a vector x
     is passed to this function, it must return an equal length y,
     and a numberof(a)-by-numberof(x) array of dyda.

     If the weight= keyword is present, it should be a list of
     sigma values for the Y data.
 */
{
  local dyda;
  dkl= array(0., numberof(a), numberof(a));
  dkl(1:numberof(dkl):numberof(a)+1)= 1.0;
  sig2= 1.0/(is_void(weight)? array(1.,numberof(y)) : weight*weight);
  lambda= 0.001;

  for (i=0 ; i<1000 ; ++i) {
    dya= y - y_and_dyda(a, x, dyda);
    chi2= sum(sig2*dya*dya);
    akl= sig2(-,)*dyda;
    bk= akl(,+)*dya(+);
    akl= (akl(,+)*dyda(,+));

    for (;;) {
      da= LUsolve(akl*(1.+lambda*dkl),bk);
      dya= y - y_and_dyda(a+da, x, dyda);
      chi2t= sum(sig2*dya*dya);
      if (chi2t < chi2) break;
      /* step failed, try again */
      lambda*= 10.;
    }

    /* step succeeded, take it */
    a+= da;

    /* check for convergence */
    if (abs(chi2t-chi2)<0.1 || abs(chi2t-chi2)<1.e-3*chi2t) break;
    lambda*= 0.1;
  }
  if (i>=1000) error, "failed to converge after 1000 iterations";

  return a;
}