-- LMfit Version 1.1,  14-Dec-2001
~
Procedure LMfit((x,y,guess,weights) fits a user-given function to a 
                                     set of data using the Levenberg-Marquard
                                     algorithm as in Bevington.
Parameters:       x:        vector for independent variables' values
                  y:        vector for dependent variables' values
                  weights:  vector with weights or scalar 1 for non-weighted fit.
                  guess:    vector of initial guesses for coeffs

Globals to set:   LM_f(coeffs,x)   the function of x with coeffs to be fitted
                  epsilon:  accuracy criterion (1e-5 is good, 1e-8 is ok).
                  Maxiter:  Max. no. of iterations
                  LM_lambda:lambda parmeter, defaults to 0.001
                  LM_fast:  set to 1 for faster, but less accurate, derivatives

Globals to read:  coeffs:   vector of final fitted coefficients
                  LM_deltas(guess):   vector or scalar of initial stepsize
                  LM_iter: returns the no. of iterations performed
                  LM_chisq[:Maxiter]: Vector of squared errors
                  LM_lambdas[:Maxiter]: Vector of lambda values

If the weights are derived from 1-sigma errors, LMfit returns the Chi-squared/deg-of-freedom
of the fit.

Needs "matrix ops" to work: I have renamed the "invert" function to "minvert" to 
avoid a conflict with SVD. I have not had good results with the invert function of SVD.
It is used in LM_New.

H.-Ulrich (Uli) Wienands, Stanford Linear Accelerator Center, uli@slac.stanford.edu

~
LMfit(xx,yy,guess,weights) = LM_Init(xx,yy,guess,weights),
(LM_alphas(xx,yy,coeffs,weights),
 LM_New(xx,yy,coeffs,weights),
 progressbar(LM_iter/Maxiter*100),
 (LM_lambda:=LM_lambda*10,
  LM_New(xx,yy,coeffs,weights)) while (LM_err1>=LM_err0),
 LM_lambda:=LM_lambda/10,
 coeffs:=LM_B,
 LM_iter:=LM_iter+1,
 LM_chisq[LM_iter]:=LM_err1,
 LM_lambdas[LM_iter]:=LM_lambda) while (abs(LM_err1-LM_err0)/LM_err0>epsilon and LM_iter<Maxiter),
LM_err1/(count(xx)-count(coeffs))

--LM_Init initializes various arrays
LM_Init(xx,yy,guess,weights) = LM_iter:=0,
LM_nterms:=count(guess),
LM_numdata:=count(xx),
LM_chisq:=0[:Maxiter],
LM_lambdas:=0[:Maxiter],
LM_k:=1,
coeffs:=guess,
LM_deltas:=Munit(LM_nterms,LM_nterms)*LM_delta*guess,
dum:=LMDelta_min,
LM_err0:=LM_Err(xx,yy,guess,weights),
LM_err1:=LM_err0*2

-- LM_alphas calculates alpha & beta matrices

LM_alphas(xx,yy,coeffs,weights) = LM_l:=1,
LM_beta:=0[:LM_nterms],
LM_alpha:=Mzero(LM_nterms,LM_nterms),
LM_err0:=LM_Err(xx,yy,coeffs,weights),
((LM_D0:=LM_f(coeffs,xx[LM_l])) when LM_fast==1,,
 LM_Derivs:=LMDeriv(coeffs,xx,LM_l),
 LM_j:=1,
 (LM_k:=1,
  (LM_alpha[LM_j,LM_k]:=LM_alpha[LM_j,LM_k]+LM_Derivs[LM_j]*LM_Derivs[LM_k]*weights[LM_l],
   LM_k:=LM_k+1) while LM_k<=LM_j,
  LM_beta[LM_j]:=LM_beta[LM_j]+weights[LM_l]*(yy[LM_l]-LM_f(coeffs,xx[LM_l]))*LM_Derivs[LM_j],
  LM_j:=LM_j+1) while LM_j<=LM_nterms,
 LM_l:=LM_l+1) while LM_l<=LM_numdata,
LM_j:=1,
(LM_k:=1,
 (LM_alpha[LM_k,LM_j]:=LM_alpha[LM_j,LM_k],
  LM_k:=LM_k+1) while LM_k<=LM_j,
 LM_j:=LM_j+1) while LM_j<=LM_nterms

-- LM_New finds new coefficients
LM_New(xx,yy,coeffs,weights) = LM_alphap:=LM_alpha,LM_B:=coeffs,LM_j:=1,
(LM_alphap[LM_j,LM_j]:=LM_alphap[LM_j,LM_j]+LM_lambda,
 LM_j:=LM_j+1) while LM_j<=LM_nterms,
LM_epsp:=minvert(LM_alphap),
LM_j:=1,
(LM_B[LM_j]:=LM_B[LM_j]+sum(LM_beta[LM_kk]*LM_epsp[LM_j,LM_kk],LM_kk,1,LM_nterms),
 LM_j:=LM_j+1) while LM_j<=LM_nterms,
LM_err1:=LM_Err(xx,yy,LM_B,weights)

-- LMDelta_f sets up step-size deltas for derivatives

LMDelta_min = LM_i:=1,
              (LM_deltas[LM_i,LM_i]:=LM_delta when LM_deltas[LM_i,LM_i]==0,, 
               LM_i:=LM_i+1) while LM_i <= LM_nterms

-- LMDeriv calculates derivatives wrt coeffs
LMDeriv(coeffs,xx,LM_numpt)[i] = (LM_f(coeffs+LM_deltas[i],xx[LM_numpt])-LM_D0)/LM_deltas[i,i] when LM_fast==1,
(LM_f(coeffs+LM_deltas[i],xx[LM_numpt])-LM_f(coeffs-LM_deltas[i],xx[LM_numpt]))/2/LM_deltas[i,i] otherwise dim[LM_nterms]

-- LM_Err calculates the sum of squared errors
LM_Err(xx,yy,a,weights) = sum(((LM_f(a,xx[j])-yy[j])^2)*weights[j],j,1,count(xx))

-- Mzero returns a zero matrix
Mzero(rows,cols)[i,j] = 0 dim[rows,cols]

-- Munit returns a unit matrix
Munit(rows,cols)[i,j] = 1 when i==j, 0 otherwise dim[rows,cols]

-- Tunable constants:
LM_lambda=0.001
LM_delta=0.001
LM_fast=0