Does any one have a working example of the Levenberg-Marquardt code that is available to see how to use this optimizer is supposed to function?

You need to include the following NUMRECs:
mrqmin.o mrqcof.o nrutil.o covsrt.o gaussj.o fgauss.o gasdev.o ran1.o

I use a Makefile, that's why I have .o-files (object-files) instead of the .c-files.

Regards
xxx


PS: My post is too late for, but might help sombody else.



/* Driver for routine mrqmin */

#include <stdio.h>
#include <math.h>
#define NRANSI
#include "nr.h"
#include "nrutil.h"

#define NPT 100
#define MA 6
#define SPREAD 0.001

int main(void)
{
long idum=(-911);
int i,*ia,iter,itst,j,k,mfit=MA;
float alamda,chisq,ochisq,*x,*y,*sig,**covar,**alpha;
static float a[MA+1]=
{0.0,5.0,2.0,3.0,2.0,5.0,3.0};
static float gues[MA+1]=
{0.0,4.5,2.2,2.8,2.5,4.9,2.8};

ia=ivector(1,MA);
x=vector(1,NPT);
y=vector(1,NPT);
sig=vector(1,NPT);
covar=matrix(1,MA,1,MA);
alpha=matrix(1,MA,1,MA);
/* First try a sum of two Gaussians */
for (i=1;i<=NPT;i++) {
x[i]=0.1*i;
y[i]=0.0;
for (j=1;j<=MA;j+=3) {
y[i] += a[j]*exp(-SQR((x[i]-a[j+1])/a[j+2]));
}
y[i] *= (1.0+SPREAD*gasdev(&idum));
sig[i]=SPREAD*y[i];
}
for (i=1;i<=mfit;i++) ia[i]=1;
for (i=1;i<=MA;i++) a[i]=gues[i];
for (iter=1;iter<=2;iter++) {
alamda = -1;
mrqmin(x,y,sig,NPT,a,ia,MA,covar,alpha,&chisq,fgauss,&alamda);
k=1;
itst=0;
for (;;) {
printf("\n%s %2d %17s %10.4f %10s %9.2e\n","Iteration #",k,
"chi-squared:",chisq,"alamda:",alamda);
printf("%8s %8s %8s %8s %8s %8s\n",
"a[1]","a[2]","a[3]","a[4]","a[5]","a[6]");
for (i=1;i<=6;i++) printf("%9.4f",a[i]);
printf("\n");
k++;
ochisq=chisq;
mrqmin(x,y,sig,NPT,a,ia,MA,covar,alpha,&chisq,fgauss,&alamda);
if (chisq > ochisq)
itst=0;
else if (fabs(ochisq-chisq) < 0.1)
itst++;
if (itst < 4) continue;
alamda=0.0;
mrqmin(x,y,sig,NPT,a,ia,MA,covar,alpha,&chisq,fgauss,&alamda);
printf("\nUncertainties:\n");
for (i=1;i<=6;i++) printf("%9.4f",sqrt(covar[i][i]));
printf("\n");
printf("\nExpected results:\n");
printf(" %7.2f %8.2f %8.2f %8.2f %8.2f %8.2f\n",
5.0,2.0,3.0,2.0,5.0,3.0);
break;
}
if (iter == 1) {
printf("press return to continue with constraint\n");
(void) getchar();
printf("holding a[2] and a[5] constant\n");
for (j=1;j<=MA;j++) a[j] += 0.1;
a[2]=2.0;
ia[2]=0;
a[5]=5.0;
ia[5]=0;
}
}
free_matrix(alpha,1,MA,1,MA);
free_matrix(covar,1,MA,1,MA);
free_vector(sig,1,NPT);
free_vector(y,1,NPT);
free_vector(x,1,NPT);
free_ivector(ia,1,MA);
return 0;
}
#undef NRANSI

I am new to programming and am trying to use LM algorithm. When writing nr.h and nrutil.h, do we use them as is? For nr.h, we write all the "#define" statements and "void" statements. And for nrutil.h, do we keep in the astericks? Thanks

xxxsemoi


Thanks for a "real example". It's nice to have something to play with to understand the interaction between all these complicated routines.

FYI, Here is a makefile that I used to compile and immediately run your example:

In this case, I call the makefile "Make_demo" and the execitable program "demo". Your driver code is in a file called "demo.c" and the makefile shown below is in the same directory as all the following .c and .h files:

covsrt.c
demo.c
fgauss.c
gasdev.c
gaussj.c
mrqcof.c
mrqmin.c
nrutil.c
ran1.c
nr.h
nrutil.h

To run the make utility and delete old objects files and executables, type

"make clean -f Make_demo"

To just compile and run the file, type

"make -f Make_demo"

------------makefile "Make_demo" ---------------------------------
CC = gcc
CFLAGS = -g -Wall
LIBS = -lm

program = demo
objects = demo.o nrutil.o covsrt.o gaussj.o fgauss.o gasdev.o ran1.o mrqcof.o mrqmin.o


all: $(program) check

$(program): $(objects)
@echo "--------------------------------------------------"
@echo "Building L-M 'demo' program"
@echo "--------------------------------------------------"
$(CC) $(CFLAGS) -o $@ $^ $(LIBS)

check:
@echo "--------------------------------------------------"
@echo "Running L-M 'demo' program"
@echo "--------------------------------------------------"
./$(program)

clean:
rm -f gmon.out *.o *~

distclean: clean
rm -f $(program)
----------------------------------------------------------------------


Thanks again,


Steve Licata



You need to include the following NUMRECs:
mrqmin.o mrqcof.o nrutil.o covsrt.o gaussj.o fgauss.o gasdev.o ran1.o

I use a Makefile, that's why I have .o-files (object-files) instead of the .c-files.

Regards
xxx


PS: My post is too late for, but might help sombody else.



/* Driver for routine mrqmin */

#include <stdio.h>
#include <math.h>
#define NRANSI
#include "nr.h"
#include "nrutil.h"

#define NPT 100
#define MA 6
#define SPREAD 0.001

int main(void)
{
long idum=(-911);
int i,*ia,iter,itst,j,k,mfit=MA;
float alamda,chisq,ochisq,*x,*y,*sig,**covar,**alpha;
static float a[MA+1]=
{0.0,5.0,2.0,3.0,2.0,5.0,3.0};
static float gues[MA+1]=
{0.0,4.5,2.2,2.8,2.5,4.9,2.8};

ia=ivector(1,MA);
x=vector(1,NPT);
y=vector(1,NPT);
sig=vector(1,NPT);
covar=matrix(1,MA,1,MA);
alpha=matrix(1,MA,1,MA);
/* First try a sum of two Gaussians */
for (i=1;i<=NPT;i++) {
x[i]=0.1*i;
y[i]=0.0;
for (j=1;j<=MA;j+=3) {
y[i] += a[j]*exp(-SQR((x[i]-a[j+1])/a[j+2]));
}
y[i] *= (1.0+SPREAD*gasdev(&idum));
sig[i]=SPREAD*y[i];
}
for (i=1;i<=mfit;i++) ia[i]=1;
for (i=1;i<=MA;i++) a[i]=gues[i];
for (iter=1;iter<=2;iter++) {
alamda = -1;
mrqmin(x,y,sig,NPT,a,ia,MA,covar,alpha,&chisq,fgauss,&alamda);
k=1;
itst=0;
for (;;) {
printf("\n%s %2d %17s %10.4f %10s %9.2e\n","Iteration #",k,
"chi-squared:",chisq,"alamda:",alamda);
printf("%8s %8s %8s %8s %8s %8s\n",
"a[1]","a[2]","a[3]","a[4]","a[5]","a[6]");
for (i=1;i<=6;i++) printf("%9.4f",a[i]);
printf("\n");
k++;
ochisq=chisq;
mrqmin(x,y,sig,NPT,a,ia,MA,covar,alpha,&chisq,fgauss,&alamda);
if (chisq > ochisq)
itst=0;
else if (fabs(ochisq-chisq) < 0.1)
itst++;
if (itst < 4) continue;
alamda=0.0;
mrqmin(x,y,sig,NPT,a,ia,MA,covar,alpha,&chisq,fgauss,&alamda);
printf("\nUncertainties:\n");
for (i=1;i<=6;i++) printf("%9.4f",sqrt(covar[i][i]));
printf("\n");
printf("\nExpected results:\n");
printf(" %7.2f %8.2f %8.2f %8.2f %8.2f %8.2f\n",
5.0,2.0,3.0,2.0,5.0,3.0);
break;
}
if (iter == 1) {
printf("press return to continue with constraint\n");
(void) getchar();
printf("holding a[2] and a[5] constant\n");
for (j=1;j<=MA;j++) a[j] += 0.1;
a[2]=2.0;
ia[2]=0;
a[5]=5.0;
ia[5]=0;
}
}
free_matrix(alpha,1,MA,1,MA);
free_matrix(covar,1,MA,1,MA);
free_vector(sig,1,NPT);
free_vector(y,1,NPT);
free_vector(x,1,NPT);
free_ivector(ia,1,MA);
return 0;
}
#undef NRANSI

In the orignal paper published by Marquardt, it is advised to scale the parameter vector and the matrix. In this code however, the parameter and the so-called Hessian matrix etc aren't scaled using the co-variance. Any idea why this is?