From a9446199f40b1b6af1ff85ddf5357873f56a7b70 Mon Sep 17 00:00:00 2001 From: Emmanuel Bertin Date: Thu, 26 Aug 2010 14:30:26 +0000 Subject: [PATCH] Fixed model-fitting issue with empty images in dual-image mode (thanks to V. de Lapparent for reporting). Updated LevMar to version 2.5. Updated URLs in various places. Pushed version number to 2.12.1. --- README | 7 +- configure | 20 +- configure.ac | 2 +- man/sex.1 | 8 +- man/sex.1.in | 6 +- src/ldactoasc.h | 2 +- src/levmar/Axb.c | 2 +- src/levmar/Axb_core.c | 452 +++++++++++++++++++++------ src/levmar/CMakeLists.txt | 14 +- src/levmar/Makefile.am | 13 +- src/levmar/Makefile.icc | 22 +- src/levmar/Makefile.in | 19 +- src/levmar/Makefile.vc | 19 +- src/levmar/README | 2 +- src/levmar/README.txt | 18 +- src/levmar/compiler.h | 4 + src/levmar/expfit.c | 2 +- src/levmar/levmar.h | 370 ++++++++++++++++++++++ src/levmar/lm.c | 2 +- src/levmar/lm.h | 120 +++++++- src/levmar/lm_core.c | 10 +- src/levmar/lmbc.c | 2 +- src/levmar/lmbc_core.c | 5 +- src/levmar/lmblec.c | 6 +- src/levmar/lmbleic.c | 89 ++++++ src/levmar/lmbleic_core.c | 506 +++++++++++++++++++++++++++++++ src/levmar/lmdemo.c | 274 ++++++++++++++--- src/levmar/lmlec.c | 6 +- src/levmar/matlab/CMakeLists.txt | 58 ++++ src/levmar/matlab/README.txt | 1 + src/levmar/matlab/jacbt3.m | 2 +- src/levmar/matlab/jacexpfit.m | 2 +- src/levmar/matlab/jachs01.m | 2 +- src/levmar/matlab/jacmeyer.m | 2 +- src/levmar/matlab/jacmodhs52.m | 2 +- src/levmar/matlab/jacmodhs76.m | 7 + src/levmar/matlab/jacosborne.m | 11 + src/levmar/matlab/levmar.c | 159 ++++++++-- src/levmar/matlab/levmar.m | 17 +- src/levmar/matlab/lmdemo.m | 43 ++- src/levmar/matlab/modhs76.m | 7 + src/levmar/matlab/osborne.m | 7 + src/levmar/misc.c | 2 +- src/levmar/misc_core.c | 47 ++- src/profit.c | 11 +- xsl/sextractor.xsl | 10 +- 46 files changed, 2095 insertions(+), 297 deletions(-) create mode 100644 src/levmar/levmar.h create mode 100644 src/levmar/lmbleic.c create mode 100644 src/levmar/lmbleic_core.c create mode 100644 src/levmar/matlab/CMakeLists.txt create mode 100644 src/levmar/matlab/jacmodhs76.m create mode 100644 src/levmar/matlab/jacosborne.m create mode 100644 src/levmar/matlab/modhs76.m create mode 100644 src/levmar/matlab/osborne.m diff --git a/README b/README index 6d7cb02..a29a3bb 100644 --- a/README +++ b/README @@ -9,9 +9,10 @@ input image is a MEF. configuration files like default.param are still needed, though. The SExtractor homepage is -http://terapix.iap.fr/soft/sextractor +http://astromatic.net/software/sextractor In case of problems, questions or suggestions related to the software, please -refer to the TERAPIX forum: -http://terapix.iap.fr/forum/ +refer to the SExtractor forum: +http://astromatic.net/forum/forumdisplay.php?fid=4 Emmanuel Bertin. + diff --git a/configure b/configure index 54decfa..ab8f499 100755 --- a/configure +++ b/configure @@ -1,6 +1,6 @@ #! /bin/sh # Guess values for system-dependent variables and create Makefiles. -# Generated by GNU Autoconf 2.63 for sextractor 2.12.0. +# Generated by GNU Autoconf 2.63 for sextractor 2.12.1. # # Report bugs to . # @@ -750,8 +750,8 @@ SHELL=${CONFIG_SHELL-/bin/sh} # Identity of this package. PACKAGE_NAME='sextractor' PACKAGE_TARNAME='sextractor' -PACKAGE_VERSION='2.12.0' -PACKAGE_STRING='sextractor 2.12.0' +PACKAGE_VERSION='2.12.1' +PACKAGE_STRING='sextractor 2.12.1' PACKAGE_BUGREPORT='bertin@iap.fr' ac_unique_file="src/makeit.c" @@ -1508,7 +1508,7 @@ if test "$ac_init_help" = "long"; then # Omit some internal or obsolete options to make the list less imposing. # This message is too long to be a string in the A/UX 3.1 sh. cat <<_ACEOF -\`configure' configures sextractor 2.12.0 to adapt to many kinds of systems. +\`configure' configures sextractor 2.12.1 to adapt to many kinds of systems. Usage: $0 [OPTION]... [VAR=VALUE]... @@ -1578,7 +1578,7 @@ fi if test -n "$ac_init_help"; then case $ac_init_help in - short | recursive ) echo "Configuration of sextractor 2.12.0:";; + short | recursive ) echo "Configuration of sextractor 2.12.1:";; esac cat <<\_ACEOF @@ -1711,7 +1711,7 @@ fi test -n "$ac_init_help" && exit $ac_status if $ac_init_version; then cat <<\_ACEOF -sextractor configure 2.12.0 +sextractor configure 2.12.1 generated by GNU Autoconf 2.63 Copyright (C) 1992, 1993, 1994, 1995, 1996, 1998, 1999, 2000, 2001, @@ -1725,7 +1725,7 @@ cat >config.log <<_ACEOF This file contains any messages produced by compilers while running configure, to aid debugging if configure makes a mistake. -It was created by sextractor $as_me 2.12.0, which was +It was created by sextractor $as_me 2.12.1, which was generated by GNU Autoconf 2.63. Invocation command line was $ $0 $@ @@ -2428,7 +2428,7 @@ fi # Define the identity of the package. PACKAGE='sextractor' - VERSION='2.12.0' + VERSION='2.12.1' cat >>confdefs.h <<_ACEOF @@ -28593,7 +28593,7 @@ exec 6>&1 # report actual input values of CONFIG_FILES etc. instead of their # values after options handling. ac_log=" -This file was extended by sextractor $as_me 2.12.0, which was +This file was extended by sextractor $as_me 2.12.1, which was generated by GNU Autoconf 2.63. Invocation command line was CONFIG_FILES = $CONFIG_FILES @@ -28656,7 +28656,7 @@ Report bugs to ." _ACEOF cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1 ac_cs_version="\\ -sextractor config.status 2.12.0 +sextractor config.status 2.12.1 configured by $0, generated by GNU Autoconf 2.63, with options \\"`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`\\" diff --git a/configure.ac b/configure.ac index d5fcfab..c929061 100644 --- a/configure.ac +++ b/configure.ac @@ -6,7 +6,7 @@ define([AC_CACHE_LOAD],) define([AC_CACHE_SAVE],) # This is your standard Bertin source code... -AC_INIT(sextractor, 2.12.0, [bertin@iap.fr]) +AC_INIT(sextractor, 2.12.1, [bertin@iap.fr]) AC_CONFIG_SRCDIR(src/makeit.c) AC_CONFIG_AUX_DIR(autoconf) AM_CONFIG_HEADER(config.h) diff --git a/man/sex.1 b/man/sex.1 index 0d97758..c8cc08a 100644 --- a/man/sex.1 +++ b/man/sex.1 @@ -1,4 +1,4 @@ -.TH SEXTRACTOR "1" "August 2010" "SWarp 2.12.0" "User Commands" +.TH SEXTRACTOR "1" "August 2010" "SWarp 2.12.1" "User Commands" .SH NAME sex \- extract a source catalog from an astronomical FITS image .SH SYNOPSIS @@ -18,7 +18,7 @@ SExtractor is a program that builds a catalogue of objects from an astronomical image. Although it is particularly oriented towards reduction of large scale galaxy-survey data, it performs rather well on moderately crowded star fields. .RE -See http://terapix.iap.fr/soft/sextractor for more details. +See http://astromatic.net/software/sextractor for more details. .SS "Operation modes:" .TP \fB\-h\fR, \fB\-\-help\fR @@ -43,5 +43,5 @@ Report bugs to . .PP The full documentation for .B SExtractor -is maintained as a Postscript manual available at -.B http://terapix.iap.fr/soft/sextractor +is maintained as a PDF manual available at +.B http://astromatic.net/software/sextractor diff --git a/man/sex.1.in b/man/sex.1.in index 06325dc..c63631b 100644 --- a/man/sex.1.in +++ b/man/sex.1.in @@ -18,7 +18,7 @@ SExtractor is a program that builds a catalogue of objects from an astronomical image. Although it is particularly oriented towards reduction of large scale galaxy-survey data, it performs rather well on moderately crowded star fields. .RE -See http://terapix.iap.fr/soft/sextractor for more details. +See http://astromatic.net/software/sextractor for more details. .SS "Operation modes:" .TP \fB\-h\fR, \fB\-\-help\fR @@ -43,5 +43,5 @@ Report bugs to . .PP The full documentation for .B SExtractor -is maintained as a Postscript manual available at -.B http://terapix.iap.fr/soft/sextractor +is maintained as a PDF manual available at +.B http://astromatic.net/software/sextractor diff --git a/src/ldactoasc.h b/src/ldactoasc.h index 9957be0..ffa0a56 100644 --- a/src/ldactoasc.h +++ b/src/ldactoasc.h @@ -30,7 +30,7 @@ #define MYVERSION VERSION #endif #define COPYRIGHT "Emmanuel BERTIN " -#define WEBSITE "http://astromatic.iap.fr/software/sextractor/" +#define WEBSITE "http://astromatic.net/software/sextractor/" #define INSTITUTE "IAP http://www.iap.fr" /*----------------------------- Physical constants --------------------------*/ diff --git a/src/levmar/Axb.c b/src/levmar/Axb.c index a2b43fb..7fce7c2 100644 --- a/src/levmar/Axb.c +++ b/src/levmar/Axb.c @@ -28,7 +28,7 @@ #include #include -#include "lm.h" +#include "levmar.h" #include "misc.h" #if !defined(LM_DBL_PREC) && !defined(LM_SNGL_PREC) diff --git a/src/levmar/Axb_core.c b/src/levmar/Axb_core.c index 9c7b88e..46b7051 100644 --- a/src/levmar/Axb_core.c +++ b/src/levmar/Axb_core.c @@ -26,6 +26,7 @@ #error This file should not be compiled directly! #endif + #ifdef LINSOLVERS_RETAIN_MEMORY #define __STATIC__ static #else @@ -36,7 +37,6 @@ /* prototypes of LAPACK routines */ - #define GEQRF LM_MK_LAPACK_NAME(geqrf) #define ORGQR LM_MK_LAPACK_NAME(orgqr) #define TRTRS LM_MK_LAPACK_NAME(trtrs) @@ -47,6 +47,8 @@ #define GETRS LM_MK_LAPACK_NAME(getrs) #define GESVD LM_MK_LAPACK_NAME(gesvd) #define GESDD LM_MK_LAPACK_NAME(gesdd) +#define SYTRF LM_MK_LAPACK_NAME(sytrf) +#define SYTRS LM_MK_LAPACK_NAME(sytrs) /* QR decomposition */ extern int GEQRF(int *m, int *n, LM_REAL *a, int *lda, LM_REAL *tau, LM_REAL *work, int *lwork, int *info); @@ -74,12 +76,17 @@ extern int GESVD(char *jobu, char *jobvt, int *m, int *n, LM_REAL *a, int *lda, extern int GESDD(char *jobz, int *m, int *n, LM_REAL *a, int *lda, LM_REAL *s, LM_REAL *u, int *ldu, LM_REAL *vt, int *ldvt, LM_REAL *work, int *lwork, int *iwork, int *info); +/* LDLt/UDUt factorization and systems solution */ +extern int SYTRF(char *uplo, int *n, LM_REAL *a, int *lda, int *ipiv, LM_REAL *work, int *lwork, int *info); +extern int SYTRS(char *uplo, int *n, int *nrhs, LM_REAL *a, int *lda, int *ipiv, LM_REAL *b, int *ldb, int *info); + /* precision-specific definitions */ #define AX_EQ_B_QR LM_ADD_PREFIX(Ax_eq_b_QR) #define AX_EQ_B_QRLS LM_ADD_PREFIX(Ax_eq_b_QRLS) #define AX_EQ_B_CHOL LM_ADD_PREFIX(Ax_eq_b_Chol) #define AX_EQ_B_LU LM_ADD_PREFIX(Ax_eq_b_LU) #define AX_EQ_B_SVD LM_ADD_PREFIX(Ax_eq_b_SVD) +#define AX_EQ_B_BK LM_ADD_PREFIX(Ax_eq_b_BK) /* * This function returns the solution of Ax = b @@ -105,8 +112,8 @@ __STATIC__ int buf_sz=0; static int nb=0; /* no __STATIC__ decl. here! */ -LM_REAL *a, *qtb, *tau, *r, *work; -int a_sz, qtb_sz, tau_sz, r_sz, tot_sz; +LM_REAL *a, *tau, *r, *work; +int a_sz, tau_sz, r_sz, tot_sz; register int i, j; int info, worksz, nrhs=1; register LM_REAL sum; @@ -115,8 +122,9 @@ register LM_REAL sum; #ifdef LINSOLVERS_RETAIN_MEMORY { if(buf) free(buf); - buf_sz=0; buf=NULL; + buf_sz=0; + return 1; } #else @@ -125,10 +133,8 @@ register LM_REAL sum; /* calculate required memory size */ a_sz=m*m; - qtb_sz=m; tau_sz=m; r_sz=m*m; /* only the upper triangular part really needed */ - if(!nb){ LM_REAL tmp; @@ -137,7 +143,7 @@ register LM_REAL sum; nb=((int)tmp)/m; // optimal worksize is m*nb } worksz=nb*m; - tot_sz=a_sz + qtb_sz + tau_sz + r_sz + worksz; + tot_sz=a_sz + tau_sz + r_sz + worksz; #ifdef LINSOLVERS_RETAIN_MEMORY if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */ @@ -160,8 +166,7 @@ register LM_REAL sum; #endif /* LINSOLVERS_RETAIN_MEMORY */ a=buf; - qtb=a+a_sz; - tau=qtb+qtb_sz; + tau=a+a_sz; r=tau+tau_sz; work=r+r_sz; @@ -209,15 +214,15 @@ register LM_REAL sum; } } - /* Q is now in a; compute Q^T b in qtb */ + /* Q is now in a; compute Q^T b in x */ for(i=0; ibuf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */ @@ -330,8 +331,7 @@ register LM_REAL sum; #endif /* LINSOLVERS_RETAIN_MEMORY */ a=buf; - atb=a+a_sz; - tau=atb+atb_sz; + tau=a+a_sz; r=tau+tau_sz; work=r+r_sz; @@ -340,11 +340,11 @@ register LM_REAL sum; for(j=0; jbuf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */ @@ -488,16 +484,15 @@ int info, nrhs=1; #endif /* LINSOLVERS_RETAIN_MEMORY */ a=buf; - b=a+a_sz; - /* store A into a anb B into b. A is assumed symmetric, + /* store A into a and B into x. A is assumed symmetric, * hence no transposition is needed */ for(i=0; ibuf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */ @@ -643,15 +633,14 @@ LM_REAL *a, *b; #endif /* LINSOLVERS_RETAIN_MEMORY */ a=buf; - b=a+a_sz; - ipiv=(int *)(b+b_sz); + ipiv=(int *)(a+a_sz); - /* store A (column major!) into a and B into b */ + /* store A (column major!) into a and B into x */ for(i=0; ibuf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */ + if(buf) free(buf); /* free previously allocated memory */ + + buf_sz=tot_sz; + buf=(LM_REAL *)malloc(buf_sz); + if(!buf){ + fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_BK) "() failed!\n"); + exit(1); + } + } +#else + buf_sz=tot_sz; + buf=(LM_REAL *)malloc(buf_sz); + if(!buf){ + fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_BK) "() failed!\n"); + exit(1); + } +#endif /* LINSOLVERS_RETAIN_MEMORY */ + + a=buf; + work=a+a_sz; + ipiv=(int *)(work+work_sz); + + /* store A into a and B into x; A is assumed to be symmetric, hence + * the column and row major order representations are the same + */ + for(i=0; ibuf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */ if(buf) free(buf); /* free previously allocated memory */ buf_sz=tot_sz; - buf=(void *)malloc(buf_sz); + buf=(LM_REAL *)malloc(buf_sz); if(!buf){ fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n"); exit(1); @@ -935,26 +1044,26 @@ LM_REAL *a, *b; } #else buf_sz=tot_sz; - buf=(void *)malloc(buf_sz); + buf=(LM_REAL *)malloc(buf_sz); if(!buf){ fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n"); exit(1); } #endif /* LINSOLVERS_RETAIN_MEMORY */ - ipiv=(int *)buf; - a=(LM_REAL *)(ipiv + ipiv_sz); - b=a+a_sz; + a=buf; + ipiv=(int *)(a+a_sz); - /* store A (column major!) into a and B into b */ + /* store A (column major!) into a and B into x */ for(i=0; ibuf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */ +// if(buf) free(buf); /* free previously allocated memory */ +// +// buf_sz=tot_sz; +// buf=(void *)malloc(tot_sz); +// if(!buf){ +// fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n"); +// exit(1); +// } +// } +//#else +// buf_sz=tot_sz; +// buf=(void *)malloc(tot_sz); +// if(!buf){ +// fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n"); +// exit(1); +// } +//#endif /* LINSOLVERS_RETAIN_MEMORY */ +// +// a=buf; +// work=a+a_sz; +// idx=(int *)(work+work_sz); +// +// /* avoid destroying A, B by copying them to a, x resp. */ +// for(i=0; imax) +// max=tmp; +// if(max==0.0){ +// fprintf(stderr, RCAT("Singular matrix A in ", AX_EQ_B_LU) "()!\n"); +//#ifndef LINSOLVERS_RETAIN_MEMORY +// free(buf); +//#endif +// +// return 0; +// } +// work[i]=LM_CNST(1.0)/max; +// } +// +// for(j=0; j=max){ +// max=tmp; +// maxi=i; +// } +// } +// if(j!=maxi){ +// for(k=0; k=0; --i){ +// sum=x[i]; +// for(j=i+1; j # http://www.ics.forth.gr/~lourakis/levmar/ noinst_LIBRARIES = liblevmar.a -liblevmar_a_SOURCES = Axb.c lmbc.c lm.c lmblec.c lmlec.c misc.c \ - compiler.h lm.h misc.h +liblevmar_a_SOURCES = Axb.c lmbc.c lm.c lmblec.c lmbleic.c lmlec.c misc.c \ + compiler.h levmar.h misc.h EXTRA_liblevmar_a_SOURCES = Axb_core.c lmbc_core.c lm_core.c \ - lmblec_core.c lmlec_core.c misc_core.c \ - LICENSE README README.txt lmdemo.c + lmblec_core.c lmbleic_core.c lmlec_core.c \ + misc_core.c \ + LICENSE README README.txt \ + Makefile.icc Makefile.vc levmar.vcproj lmdemo.vcproj \ + expfit.c lmdemo.c lm.h matlab diff --git a/src/levmar/Makefile.icc b/src/levmar/Makefile.icc index 85bb0d5..a2ab2e7 100644 --- a/src/levmar/Makefile.icc +++ b/src/levmar/Makefile.icc @@ -11,8 +11,8 @@ ARCHFLAGS=-march=pentium4 -mcpu=pentium4 CFLAGS=$(CONFIGFLAGS) $(ARCHFLAGS) -O3 -tpp7 -xW -ip -ipo -unroll #-g LAPACKLIBS_PATH=/usr/local/lib # WHEN USING LAPACK, CHANGE THIS TO WHERE YOUR COMPILED LIBS ARE! LDFLAGS=-L$(LAPACKLIBS_PATH) -L. -LIBOBJS=lm.o Axb.o misc.o lmlec.o lmbc.o lmblec.o -LIBSRCS=lm.c Axb.c misc.c lmlec.c lmbc.c lmblec.c +LIBOBJS=lm.o Axb.o misc.o lmlec.o lmbc.o lmblec.o lmbleic.o +LIBSRCS=lm.c Axb.c misc.c lmlec.c lmbc.c lmblec.c lmbleic.c DEMOBJS=lmdemo.o DEMOSRCS=lmdemo.c AR=xiar @@ -24,6 +24,9 @@ LAPACKLIBS=-llapack -lblas -lf2c # comment this line if you are not using LAPACK # The following works with the ATLAS updated lapack and Linux_P4SSE2 from http://www.netlib.org/atlas/archives/linux/ #LAPACKLIBS=-L/usr/local/atlas/lib -llapack -lcblas -lf77blas -latlas -lf2c +#LAPACKLIBS=-llapack -lgoto2 -lpthread -lf2c # This works with GotoBLAS + # from http://www.tacc.utexas.edu/research-development/tacc-projects/ + LIBS=$(LAPACKLIBS) all: liblevmar.a lmdemo @@ -35,14 +38,15 @@ liblevmar.a: $(LIBOBJS) lmdemo: $(DEMOBJS) liblevmar.a $(CC) $(ARCHFLAGS) $(LDFLAGS) $(DEMOBJS) -o lmdemo -llevmar $(LIBS) -lm -lm.o: lm.c lm_core.c lm.h misc.h compiler.h -Axb.o: Axb.c Axb_core.c lm.h misc.h -misc.o: misc.c misc_core.c lm.h misc.h -lmlec.o: lmlec.c lmlec_core.c lm.h misc.h -lmbc.o: lmbc.c lmbc_core.c lm.h misc.h compiler.h -lmblec.o: lmblec.c lmblec_core.c lm.h misc.h +lm.o: lm.c lm_core.c levmar.h misc.h compiler.h +Axb.o: Axb.c Axb_core.c levmar.h misc.h +misc.o: misc.c misc_core.c levmar.h misc.h +lmlec.o: lmlec.c lmlec_core.c levmar.h misc.h +lmbc.o: lmbc.c lmbc_core.c levmar.h misc.h compiler.h +lmblec.o: lmblec.c lmblec_core.c levmar.h misc.h +lmbleic.o: lmbleic.c lmbleic_core.c levmar.h misc.h -lmdemo.o: lm.h +lmdemo.o: levmar.h clean: @rm -f $(LIBOBJS) $(DEMOBJS) diff --git a/src/levmar/Makefile.in b/src/levmar/Makefile.in index 236a594..78aafd0 100644 --- a/src/levmar/Makefile.in +++ b/src/levmar/Makefile.in @@ -50,7 +50,8 @@ ARFLAGS = cru liblevmar_a_AR = $(AR) $(ARFLAGS) liblevmar_a_LIBADD = am_liblevmar_a_OBJECTS = Axb.$(OBJEXT) lmbc.$(OBJEXT) lm.$(OBJEXT) \ - lmblec.$(OBJEXT) lmlec.$(OBJEXT) misc.$(OBJEXT) + lmblec.$(OBJEXT) lmbleic.$(OBJEXT) lmlec.$(OBJEXT) \ + misc.$(OBJEXT) liblevmar_a_OBJECTS = $(am_liblevmar_a_OBJECTS) DEFAULT_INCLUDES = -I.@am__isrc@ -I$(top_builddir) depcomp = $(SHELL) $(top_srcdir)/autoconf/depcomp @@ -195,16 +196,19 @@ top_builddir = @top_builddir@ top_srcdir = @top_srcdir@ # Makefile.am for the levmar Levenberg-Marquardt fitting library -# Copyright (C) 2008 Emmanuel Bertin. +# Copyright (C) 2007-2010 Emmanuel Bertin. # levmar code by Manolis Lourakis # http://www.ics.forth.gr/~lourakis/levmar/ noinst_LIBRARIES = liblevmar.a -liblevmar_a_SOURCES = Axb.c lmbc.c lm.c lmblec.c lmlec.c misc.c \ - compiler.h lm.h misc.h +liblevmar_a_SOURCES = Axb.c lmbc.c lm.c lmblec.c lmbleic.c lmlec.c misc.c \ + compiler.h levmar.h misc.h EXTRA_liblevmar_a_SOURCES = Axb_core.c lmbc_core.c lm_core.c \ - lmblec_core.c lmlec_core.c misc_core.c \ - LICENSE README README.txt lmdemo.c + lmblec_core.c lmbleic_core.c lmlec_core.c \ + misc_core.c \ + LICENSE README README.txt \ + Makefile.icc Makefile.vc levmar.vcproj lmdemo.vcproj \ + expfit.c lmdemo.c lm.h matlab all: all-am @@ -255,12 +259,15 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/Axb.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/Axb_core.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/expfit.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lm.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lm_core.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lmbc.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lmbc_core.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lmblec.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lmblec_core.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lmbleic.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lmbleic_core.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lmdemo.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lmlec.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/lmlec_core.Po@am__quote@ diff --git a/src/levmar/Makefile.vc b/src/levmar/Makefile.vc index 52e3910..b8ae99a 100644 --- a/src/levmar/Makefile.vc +++ b/src/levmar/Makefile.vc @@ -21,8 +21,8 @@ CONFIGFLAGS=#/ULINSOLVERS_RETAIN_MEMORY CFLAGS=$(CONFIGFLAGS) /I. /MD /W3 /EHsc /O2 $(SPOPTFLAGS) # /Wall LAPACKLIBS_PATH=C:\src\lib # WHEN USING LAPACK, CHANGE THIS TO WHERE YOUR COMPILED LIBS ARE! LDFLAGS=/link /subsystem:console /opt:ref /libpath:$(LAPACKLIBS_PATH) /libpath:. -LIBOBJS=lm.obj Axb.obj misc.obj lmlec.obj lmbc.obj lmblec.obj -LIBSRCS=lm.c Axb.c misc.c lmlec.c lmbc.c lmblec.c +LIBOBJS=lm.obj Axb.obj misc.obj lmlec.obj lmbc.obj lmblec.obj lmbleic.obj +LIBSRCS=lm.c Axb.c misc.c lmlec.c lmbc.c lmblec.c lmbleic.c DEMOBJS=lmdemo.obj DEMOSRCS=lmdemo.c AR=lib /nologo @@ -40,14 +40,15 @@ levmar.lib: $(LIBOBJS) lmdemo.exe: $(DEMOBJS) levmar.lib $(CC) $(DEMOBJS) $(LDFLAGS) /out:lmdemo.exe $(LIBS) -lm.obj: lm.c lm_core.c lm.h misc.h compiler.h -Axb.obj: Axb.c Axb_core.c lm.h misc.h -misc.obj: misc.c misc_core.c lm.h misc.h -lmlec.obj: lmlec.c lmlec_core.c lm.h misc.h -lmbc.obj: lmbc.c lmbc_core.c lm.h misc.h compiler.h -lmblec.obj: lmblec.c lmblec_core.c lm.h misc.h +lm.obj: lm.c lm_core.c levmar.h misc.h compiler.h +Axb.obj: Axb.c Axb_core.c levmar.h misc.h +misc.obj: misc.c misc_core.c levmar.h misc.h +lmlec.obj: lmlec.c lmlec_core.c levmar.h misc.h +lmbc.obj: lmbc.c lmbc_core.c levmar.h misc.h compiler.h +lmblec.obj: lmblec.c lmblec_core.c levmar.h misc.h +lmbleic.obj: lmbleic.c lmbleic_core.c levmar.h misc.h -lmdemo.obj: lm.h +lmdemo.obj: levmar.h clean: -del $(LIBOBJS) $(DEMOBJS) diff --git a/src/levmar/README b/src/levmar/README index f5dbc85..8abf809 100644 --- a/src/levmar/README +++ b/src/levmar/README @@ -1,4 +1,4 @@ -The levmar v2.4 library has been included in this package untouched, except for +The levmar v2.5 library has been included in this package untouched, except for three warnings removed, LU decomposition replaced with calls to ATLAS-Lapack routines, Hessian matrix inversion done with SVD, and an AutoMakefile added. Emmanuel Bertin diff --git a/src/levmar/README.txt b/src/levmar/README.txt index 45d2bfa..44da929 100644 --- a/src/levmar/README.txt +++ b/src/levmar/README.txt @@ -1,6 +1,6 @@ ************************************************************** LEVMAR - version 2.4 + version 2.5 By Manolis Lourakis Institute of Computer Science @@ -12,15 +12,15 @@ GENERAL This is levmar, a copylefted C/C++ implementation of the Levenberg-Marquardt non-linear least squares algorithm. levmar includes double and single precision LM versions, both -with analytic and finite difference approximated jacobians. levmar also has some support -for constrained non-linear least squares, allowing linear equation and box constraints. -You have the following options regarding the solution of the underlying augmented normal -equations: +with analytic and finite difference approximated Jacobians. levmar also has some support +for constrained non-linear least squares, allowing linear equation, box and linear +inequality constraints. The following options regarding the solution of the underlying +augmented normal equations are offered: 1) Assuming that you have LAPACK (or an equivalent vendor library such as ESSL, MKL, NAG, ...) installed, you can use the included LAPACK-based solvers (default). -2) If you don't have LAPACK or decide not to use it, undefine HAVE_LAPACK in lm.h +2) If you don't have LAPACK or decide not to use it, undefine HAVE_LAPACK in levmar.h and a LAPACK-free, LU-based linear systems solver will by used. Also, the line setting the variable LAPACKLIBS in the Makefile should be commented out. @@ -38,12 +38,12 @@ LICENSE levmar is released under the GNU Public License (GPL), which can be found in the included LICENSE file. Note that under the terms of GPL, commercial use is allowed only if a software employing levmar is also published in source under the GPL. However, if you are interested -in using levmar in a proprietary commercial apprlication, a commercial license for levmar +in using levmar in a proprietary commercial application, a commercial license for levmar can be obtained by contacting the author using the email address at the end of this file. COMPILATION - You might first consider setting a few configuration options at the top of - lm.h. See the accompanying comments for more details. + levmar.h. See the accompanying comments for more details. - On a Linux/Unix system, typing "make" will build both levmar and the demo program using gcc. Alternatively, if Intel's C++ compiler is installed, it @@ -63,7 +63,7 @@ COMPILATION for details. MATLAB INTERFACE -Since version 2.2, the levmar distrubution includes a matlab interface. +Since version 2.2, the levmar distribution includes a matlab interface. See the 'matlab' subdirectory for more information and examples of use. Notice that *_core.c files are not to be compiled directly; For example, diff --git a/src/levmar/compiler.h b/src/levmar/compiler.h index c15e220..706924c 100644 --- a/src/levmar/compiler.h +++ b/src/levmar/compiler.h @@ -38,4 +38,8 @@ #define LM_FINITE finite // other than MSVC, ICC, GCC, let's hope this will work #endif +#ifdef _MSC_VER // avoid deprecation warnings in VS2005 +#define _CRT_SECURE_NO_WARNINGS +#endif + #endif /* _COMPILER_H_ */ diff --git a/src/levmar/expfit.c b/src/levmar/expfit.c index 43790eb..21d54f4 100644 --- a/src/levmar/expfit.c +++ b/src/levmar/expfit.c @@ -23,7 +23,7 @@ #include #include -#include +#include #ifndef LM_DBL_PREC #error Example program assumes that levmar has been compiled with double precision, see LM_DBL_PREC! diff --git a/src/levmar/levmar.h b/src/levmar/levmar.h new file mode 100644 index 0000000..34d718b --- /dev/null +++ b/src/levmar/levmar.h @@ -0,0 +1,370 @@ +/* +//////////////////////////////////////////////////////////////////////////////////// +// +// Prototypes and definitions for the Levenberg - Marquardt minimization algorithm +// Copyright (C) 2004 Manolis Lourakis (lourakis at ics forth gr) +// Institute of Computer Science, Foundation for Research & Technology - Hellas +// Heraklion, Crete, Greece. +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +//////////////////////////////////////////////////////////////////////////////////// +*/ + +#ifndef _LEVMAR_H_ +#define _LEVMAR_H_ + + +/************************************* Start of configuration options *************************************/ + +/* specify whether to use LAPACK or not. The first option is strongly recommended */ +/*#define HAVE_LAPACK*/ /* use LAPACK */ +#undef HAVE_LAPACK /* uncomment this to force not using LAPACK */ + +/* to avoid the overhead of repeated mallocs(), routines in Axb.c can be instructed to + * retain working memory between calls. Such a choice, however, renders these routines + * non-reentrant and is not safe in a shared memory multiprocessing environment. + * Bellow, this option is turned on only when not compiling with OpenMP. + */ +#if !defined(_OPENMP) +#define LINSOLVERS_RETAIN_MEMORY /* comment this if you don't want routines in Axb.c retain working memory between calls */ +#endif + +/* determine the precision variants to be build. Default settings build + * both the single and double precision routines + */ +#define LM_DBL_PREC /* comment this if you don't want the double precision routines to be compiled */ +#define LM_SNGL_PREC /* comment this if you don't want the single precision routines to be compiled */ + +/****************** End of configuration options, no changes necessary beyond this point ******************/ + + +#ifdef __cplusplus +extern "C" { +#endif + + +#define FABS(x) (((x)>=0.0)? (x) : -(x)) + +/* work arrays size for ?levmar_der and ?levmar_dif functions. + * should be multiplied by sizeof(double) or sizeof(float) to be converted to bytes + */ +#define LM_DER_WORKSZ(npar, nmeas) (2*(nmeas) + 4*(npar) + (nmeas)*(npar) + (npar)*(npar)) +#define LM_DIF_WORKSZ(npar, nmeas) (4*(nmeas) + 4*(npar) + (nmeas)*(npar) + (npar)*(npar)) + +/* work arrays size for ?levmar_bc_der and ?levmar_bc_dif functions. + * should be multiplied by sizeof(double) or sizeof(float) to be converted to bytes + */ +#define LM_BC_DER_WORKSZ(npar, nmeas) (2*(nmeas) + 4*(npar) + (nmeas)*(npar) + (npar)*(npar)) +#define LM_BC_DIF_WORKSZ(npar, nmeas) LM_BC_DER_WORKSZ((npar), (nmeas)) /* LEVMAR_BC_DIF currently implemented using LEVMAR_BC_DER()! */ + +/* work arrays size for ?levmar_lec_der and ?levmar_lec_dif functions. + * should be multiplied by sizeof(double) or sizeof(float) to be converted to bytes + */ +#define LM_LEC_DER_WORKSZ(npar, nmeas, nconstr) LM_DER_WORKSZ((npar)-(nconstr), (nmeas)) +#define LM_LEC_DIF_WORKSZ(npar, nmeas, nconstr) LM_DIF_WORKSZ((npar)-(nconstr), (nmeas)) + +/* work arrays size for ?levmar_blec_der and ?levmar_blec_dif functions. + * should be multiplied by sizeof(double) or sizeof(float) to be converted to bytes + */ +#define LM_BLEC_DER_WORKSZ(npar, nmeas, nconstr) LM_LEC_DER_WORKSZ((npar), (nmeas)+(npar), (nconstr)) +#define LM_BLEC_DIF_WORKSZ(npar, nmeas, nconstr) LM_LEC_DIF_WORKSZ((npar), (nmeas)+(npar), (nconstr)) + +/* work arrays size for ?levmar_bleic_der and ?levmar_bleic_dif functions. + * should be multiplied by sizeof(double) or sizeof(float) to be converted to bytes + */ +#define LM_BLEIC_DER_WORKSZ(npar, nmeas, nconstr1, nconstr2) LM_BLEC_DER_WORKSZ((npar)+(nconstr2), (nmeas)+(nconstr2), (nconstr1)+(nconstr2)) +#define LM_BLEIC_DIF_WORKSZ(npar, nmeas, nconstr1, nconstr2) LM_BLEC_DIF_WORKSZ((npar)+(nconstr2), (nmeas)+(nconstr2), (nconstr1)+(nconstr2)) + +#define LM_OPTS_SZ 5 /* max(4, 5) */ +#define LM_INFO_SZ 10 +#define LM_ERROR -1 +#define LM_INIT_MU 1E-03 +#define LM_STOP_THRESH 1E-17 +#define LM_DIFF_DELTA 1E-06 +#define LM_VERSION "2.5 (December 2009)" + +#ifdef LM_DBL_PREC +/* double precision LM, with & without Jacobian */ +/* unconstrained minimization */ +extern int dlevmar_der( + void (*func)(double *p, double *hx, int m, int n, void *adata), + void (*jacf)(double *p, double *j, int m, int n, void *adata), + double *p, double *x, int m, int n, int itmax, double *opts, + double *info, double *work, double *covar, void *adata); + +extern int dlevmar_dif( + void (*func)(double *p, double *hx, int m, int n, void *adata), + double *p, double *x, int m, int n, int itmax, double *opts, + double *info, double *work, double *covar, void *adata); + +/* box-constrained minimization */ +extern int dlevmar_bc_der( + void (*func)(double *p, double *hx, int m, int n, void *adata), + void (*jacf)(double *p, double *j, int m, int n, void *adata), + double *p, double *x, int m, int n, double *lb, double *ub, + int itmax, double *opts, double *info, double *work, double *covar, void *adata); + +extern int dlevmar_bc_dif( + void (*func)(double *p, double *hx, int m, int n, void *adata), + double *p, double *x, int m, int n, double *lb, double *ub, + int itmax, double *opts, double *info, double *work, double *covar, void *adata); + +#ifdef HAVE_LAPACK +/* linear equation constrained minimization */ +extern int dlevmar_lec_der( + void (*func)(double *p, double *hx, int m, int n, void *adata), + void (*jacf)(double *p, double *j, int m, int n, void *adata), + double *p, double *x, int m, int n, double *A, double *b, int k, + int itmax, double *opts, double *info, double *work, double *covar, void *adata); + +extern int dlevmar_lec_dif( + void (*func)(double *p, double *hx, int m, int n, void *adata), + double *p, double *x, int m, int n, double *A, double *b, int k, + int itmax, double *opts, double *info, double *work, double *covar, void *adata); + +/* box & linear equation constrained minimization */ +extern int dlevmar_blec_der( + void (*func)(double *p, double *hx, int m, int n, void *adata), + void (*jacf)(double *p, double *j, int m, int n, void *adata), + double *p, double *x, int m, int n, double *lb, double *ub, double *A, double *b, int k, double *wghts, + int itmax, double *opts, double *info, double *work, double *covar, void *adata); + +extern int dlevmar_blec_dif( + void (*func)(double *p, double *hx, int m, int n, void *adata), + double *p, double *x, int m, int n, double *lb, double *ub, double *A, double *b, int k, double *wghts, + int itmax, double *opts, double *info, double *work, double *covar, void *adata); + +/* box, linear equations & inequalities constrained minimization */ +extern int dlevmar_bleic_der( + void (*func)(double *p, double *hx, int m, int n, void *adata), + void (*jacf)(double *p, double *j, int m, int n, void *adata), + double *p, double *x, int m, int n, double *lb, double *ub, + double *A, double *b, int k1, double *C, double *d, int k2, + int itmax, double *opts, double *info, double *work, double *covar, void *adata); + +extern int dlevmar_bleic_dif( + void (*func)(double *p, double *hx, int m, int n, void *adata), + double *p, double *x, int m, int n, double *lb, double *ub, + double *A, double *b, int k1, double *C, double *d, int k2, + int itmax, double *opts, double *info, double *work, double *covar, void *adata); + +/* box & linear inequality constraints */ +extern int dlevmar_blic_der( + void (*func)(double *p, double *hx, int m, int n, void *adata), + void (*jacf)(double *p, double *j, int m, int n, void *adata), + double *p, double *x, int m, int n, double *lb, double *ub, double *C, double *d, int k2, + int itmax, double opts[4], double info[LM_INFO_SZ], double *work, double *covar, void *adata); + +extern int dlevmar_blic_dif( + void (*func)(double *p, double *hx, int m, int n, void *adata), + double *p, double *x, int m, int n, double *lb, double *ub, double *C, double *d, int k2, + int itmax, double opts[5], double info[LM_INFO_SZ], double *work, double *covar, void *adata); + +/* linear equation & inequality constraints */ +extern int dlevmar_leic_der( + void (*func)(double *p, double *hx, int m, int n, void *adata), + void (*jacf)(double *p, double *j, int m, int n, void *adata), + double *p, double *x, int m, int n, double *A, double *b, int k1, double *C, double *d, int k2, + int itmax, double opts[4], double info[LM_INFO_SZ], double *work, double *covar, void *adata); + +extern int dlevmar_leic_dif( + void (*func)(double *p, double *hx, int m, int n, void *adata), + double *p, double *x, int m, int n, double *A, double *b, int k1, double *C, double *d, int k2, + int itmax, double opts[5], double info[LM_INFO_SZ], double *work, double *covar, void *adata); + +/* linear inequality constraints */ +extern int dlevmar_lic_der( + void (*func)(double *p, double *hx, int m, int n, void *adata), + void (*jacf)(double *p, double *j, int m, int n, void *adata), + double *p, double *x, int m, int n, double *C, double *d, int k2, + int itmax, double opts[4], double info[LM_INFO_SZ], double *work, double *covar, void *adata); + +extern int dlevmar_lic_dif( + void (*func)(double *p, double *hx, int m, int n, void *adata), + double *p, double *x, int m, int n, double *C, double *d, int k2, + int itmax, double opts[5], double info[LM_INFO_SZ], double *work, double *covar, void *adata); +#endif /* HAVE_LAPACK */ + +#endif /* LM_DBL_PREC */ + + +#ifdef LM_SNGL_PREC +/* single precision LM, with & without Jacobian */ +/* unconstrained minimization */ +extern int slevmar_der( + void (*func)(float *p, float *hx, int m, int n, void *adata), + void (*jacf)(float *p, float *j, int m, int n, void *adata), + float *p, float *x, int m, int n, int itmax, float *opts, + float *info, float *work, float *covar, void *adata); + +extern int slevmar_dif( + void (*func)(float *p, float *hx, int m, int n, void *adata), + float *p, float *x, int m, int n, int itmax, float *opts, + float *info, float *work, float *covar, void *adata); + +/* box-constrained minimization */ +extern int slevmar_bc_der( + void (*func)(float *p, float *hx, int m, int n, void *adata), + void (*jacf)(float *p, float *j, int m, int n, void *adata), + float *p, float *x, int m, int n, float *lb, float *ub, + int itmax, float *opts, float *info, float *work, float *covar, void *adata); + +extern int slevmar_bc_dif( + void (*func)(float *p, float *hx, int m, int n, void *adata), + float *p, float *x, int m, int n, float *lb, float *ub, + int itmax, float *opts, float *info, float *work, float *covar, void *adata); + +#ifdef HAVE_LAPACK +/* linear equation constrained minimization */ +extern int slevmar_lec_der( + void (*func)(float *p, float *hx, int m, int n, void *adata), + void (*jacf)(float *p, float *j, int m, int n, void *adata), + float *p, float *x, int m, int n, float *A, float *b, int k, + int itmax, float *opts, float *info, float *work, float *covar, void *adata); + +extern int slevmar_lec_dif( + void (*func)(float *p, float *hx, int m, int n, void *adata), + float *p, float *x, int m, int n, float *A, float *b, int k, + int itmax, float *opts, float *info, float *work, float *covar, void *adata); + +/* box & linear equation constrained minimization */ +extern int slevmar_blec_der( + void (*func)(float *p, float *hx, int m, int n, void *adata), + void (*jacf)(float *p, float *j, int m, int n, void *adata), + float *p, float *x, int m, int n, float *lb, float *ub, float *A, float *b, int k, float *wghts, + int itmax, float *opts, float *info, float *work, float *covar, void *adata); + +extern int slevmar_blec_dif( + void (*func)(float *p, float *hx, int m, int n, void *adata), + float *p, float *x, int m, int n, float *lb, float *ub, float *A, float *b, int k, float *wghts, + int itmax, float *opts, float *info, float *work, float *covar, void *adata); + +/* box, linear equations & inequalities constrained minimization */ +extern int slevmar_bleic_der( + void (*func)(float *p, float *hx, int m, int n, void *adata), + void (*jacf)(float *p, float *j, int m, int n, void *adata), + float *p, float *x, int m, int n, float *lb, float *ub, + float *A, float *b, int k1, float *C, float *d, int k2, + int itmax, float *opts, float *info, float *work, float *covar, void *adata); + +extern int slevmar_bleic_dif( + void (*func)(float *p, float *hx, int m, int n, void *adata), + float *p, float *x, int m, int n, float *lb, float *ub, + float *A, float *b, int k1, float *C, float *d, int k2, + int itmax, float *opts, float *info, float *work, float *covar, void *adata); + +/* box & linear inequality constraints */ +extern int slevmar_blic_der( + void (*func)(float *p, float *hx, int m, int n, void *adata), + void (*jacf)(float *p, float *j, int m, int n, void *adata), + float *p, float *x, int m, int n, float *lb, float *ub, float *C, float *d, int k2, + int itmax, float opts[4], float info[LM_INFO_SZ], float *work, float *covar, void *adata); + +extern int slevmar_blic_dif( + void (*func)(float *p, float *hx, int m, int n, void *adata), + float *p, float *x, int m, int n, float *lb, float *ub, float *C, float *d, int k2, + int itmax, float opts[5], float info[LM_INFO_SZ], float *work, float *covar, void *adata); + +/* linear equality & inequality constraints */ +extern int slevmar_leic_der( + void (*func)(float *p, float *hx, int m, int n, void *adata), + void (*jacf)(float *p, float *j, int m, int n, void *adata), + float *p, float *x, int m, int n, float *A, float *b, int k1, float *C, float *d, int k2, + int itmax, float opts[4], float info[LM_INFO_SZ], float *work, float *covar, void *adata); + +extern int slevmar_leic_dif( + void (*func)(float *p, float *hx, int m, int n, void *adata), + float *p, float *x, int m, int n, float *A, float *b, int k1, float *C, float *d, int k2, + int itmax, float opts[5], float info[LM_INFO_SZ], float *work, float *covar, void *adata); + +/* linear inequality constraints */ +extern int slevmar_lic_der( + void (*func)(float *p, float *hx, int m, int n, void *adata), + void (*jacf)(float *p, float *j, int m, int n, void *adata), + float *p, float *x, int m, int n, float *C, float *d, int k2, + int itmax, float opts[4], float info[LM_INFO_SZ], float *work, float *covar, void *adata); + +extern int slevmar_lic_dif( + void (*func)(float *p, float *hx, int m, int n, void *adata), + float *p, float *x, int m, int n, float *C, float *d, int k2, + int itmax, float opts[5], float info[LM_INFO_SZ], float *work, float *covar, void *adata); +#endif /* HAVE_LAPACK */ + +#endif /* LM_SNGL_PREC */ + +/* linear system solvers */ +#ifdef HAVE_LAPACK + +#ifdef LM_DBL_PREC +extern int dAx_eq_b_QR(double *A, double *B, double *x, int m); +extern int dAx_eq_b_QRLS(double *A, double *B, double *x, int m, int n); +extern int dAx_eq_b_Chol(double *A, double *B, double *x, int m); +extern int dAx_eq_b_LU(double *A, double *B, double *x, int m); +extern int dAx_eq_b_SVD(double *A, double *B, double *x, int m); +extern int dAx_eq_b_BK(double *A, double *B, double *x, int m); +#endif /* LM_DBL_PREC */ + +#ifdef LM_SNGL_PREC +extern int sAx_eq_b_QR(float *A, float *B, float *x, int m); +extern int sAx_eq_b_QRLS(float *A, float *B, float *x, int m, int n); +extern int sAx_eq_b_Chol(float *A, float *B, float *x, int m); +extern int sAx_eq_b_LU(float *A, float *B, float *x, int m); +extern int sAx_eq_b_SVD(float *A, float *B, float *x, int m); +extern int sAx_eq_b_BK(float *A, float *B, float *x, int m); +#endif /* LM_SNGL_PREC */ + +#else /* no LAPACK */ + +#ifdef LM_DBL_PREC +extern int dAx_eq_b_LU_noLapack(double *A, double *B, double *x, int n); +#endif /* LM_DBL_PREC */ + +#ifdef LM_SNGL_PREC +extern int sAx_eq_b_LU_noLapack(float *A, float *B, float *x, int n); +#endif /* LM_SNGL_PREC */ + +#endif /* HAVE_LAPACK */ + +/* Jacobian verification, double & single precision */ +#ifdef LM_DBL_PREC +extern void dlevmar_chkjac( + void (*func)(double *p, double *hx, int m, int n, void *adata), + void (*jacf)(double *p, double *j, int m, int n, void *adata), + double *p, int m, int n, void *adata, double *err); +#endif /* LM_DBL_PREC */ + +#ifdef LM_SNGL_PREC +extern void slevmar_chkjac( + void (*func)(float *p, float *hx, int m, int n, void *adata), + void (*jacf)(float *p, float *j, int m, int n, void *adata), + float *p, int m, int n, void *adata, float *err); +#endif /* LM_SNGL_PREC */ + +/* standard deviation, coefficient of determination (R2) & Pearson's correlation coefficient for best-fit parameters */ +#ifdef LM_DBL_PREC +extern double dlevmar_stddev( double *covar, int m, int i); +extern double dlevmar_corcoef(double *covar, int m, int i, int j); +extern double dlevmar_R2(void (*func)(double *p, double *hx, int m, int n, void *adata), double *p, double *x, int m, int n, void *adata); + +#endif /* LM_DBL_PREC */ + +#ifdef LM_SNGL_PREC +extern float slevmar_stddev( float *covar, int m, int i); +extern float slevmar_corcoef(float *covar, int m, int i, int j); +extern float slevmar_R2(void (*func)(float *p, float *hx, int m, int n, void *adata), float *p, float *x, int m, int n, void *adata); +#endif /* LM_SNGL_PREC */ + +#ifdef __cplusplus +} +#endif + +#endif /* _LEVMAR_H_ */ diff --git a/src/levmar/lm.c b/src/levmar/lm.c index 48839a5..81d8acb 100644 --- a/src/levmar/lm.c +++ b/src/levmar/lm.c @@ -28,7 +28,7 @@ #include #include -#include "lm.h" +#include "levmar.h" #include "compiler.h" #include "misc.h" diff --git a/src/levmar/lm.h b/src/levmar/lm.h index 83a39f0..34d718b 100644 --- a/src/levmar/lm.h +++ b/src/levmar/lm.h @@ -19,15 +19,15 @@ //////////////////////////////////////////////////////////////////////////////////// */ -#ifndef _LM_H_ -#define _LM_H_ +#ifndef _LEVMAR_H_ +#define _LEVMAR_H_ /************************************* Start of configuration options *************************************/ /* specify whether to use LAPACK or not. The first option is strongly recommended */ -#define HAVE_LAPACK /* use LAPACK */ -#undef HAVE_LAPACK /* uncomment this to force not using LAPACK */ +/*#define HAVE_LAPACK*/ /* use LAPACK */ +#undef HAVE_LAPACK /* uncomment this to force not using LAPACK */ /* to avoid the overhead of repeated mallocs(), routines in Axb.c can be instructed to * retain working memory between calls. Such a choice, however, renders these routines @@ -78,13 +78,19 @@ extern "C" { #define LM_BLEC_DER_WORKSZ(npar, nmeas, nconstr) LM_LEC_DER_WORKSZ((npar), (nmeas)+(npar), (nconstr)) #define LM_BLEC_DIF_WORKSZ(npar, nmeas, nconstr) LM_LEC_DIF_WORKSZ((npar), (nmeas)+(npar), (nconstr)) +/* work arrays size for ?levmar_bleic_der and ?levmar_bleic_dif functions. + * should be multiplied by sizeof(double) or sizeof(float) to be converted to bytes + */ +#define LM_BLEIC_DER_WORKSZ(npar, nmeas, nconstr1, nconstr2) LM_BLEC_DER_WORKSZ((npar)+(nconstr2), (nmeas)+(nconstr2), (nconstr1)+(nconstr2)) +#define LM_BLEIC_DIF_WORKSZ(npar, nmeas, nconstr1, nconstr2) LM_BLEC_DIF_WORKSZ((npar)+(nconstr2), (nmeas)+(nconstr2), (nconstr1)+(nconstr2)) + #define LM_OPTS_SZ 5 /* max(4, 5) */ #define LM_INFO_SZ 10 #define LM_ERROR -1 #define LM_INIT_MU 1E-03 #define LM_STOP_THRESH 1E-17 #define LM_DIFF_DELTA 1E-06 -#define LM_VERSION "2.4 (April 2009)" +#define LM_VERSION "2.5 (December 2009)" #ifdef LM_DBL_PREC /* double precision LM, with & without Jacobian */ @@ -136,6 +142,56 @@ extern int dlevmar_blec_dif( void (*func)(double *p, double *hx, int m, int n, void *adata), double *p, double *x, int m, int n, double *lb, double *ub, double *A, double *b, int k, double *wghts, int itmax, double *opts, double *info, double *work, double *covar, void *adata); + +/* box, linear equations & inequalities constrained minimization */ +extern int dlevmar_bleic_der( + void (*func)(double *p, double *hx, int m, int n, void *adata), + void (*jacf)(double *p, double *j, int m, int n, void *adata), + double *p, double *x, int m, int n, double *lb, double *ub, + double *A, double *b, int k1, double *C, double *d, int k2, + int itmax, double *opts, double *info, double *work, double *covar, void *adata); + +extern int dlevmar_bleic_dif( + void (*func)(double *p, double *hx, int m, int n, void *adata), + double *p, double *x, int m, int n, double *lb, double *ub, + double *A, double *b, int k1, double *C, double *d, int k2, + int itmax, double *opts, double *info, double *work, double *covar, void *adata); + +/* box & linear inequality constraints */ +extern int dlevmar_blic_der( + void (*func)(double *p, double *hx, int m, int n, void *adata), + void (*jacf)(double *p, double *j, int m, int n, void *adata), + double *p, double *x, int m, int n, double *lb, double *ub, double *C, double *d, int k2, + int itmax, double opts[4], double info[LM_INFO_SZ], double *work, double *covar, void *adata); + +extern int dlevmar_blic_dif( + void (*func)(double *p, double *hx, int m, int n, void *adata), + double *p, double *x, int m, int n, double *lb, double *ub, double *C, double *d, int k2, + int itmax, double opts[5], double info[LM_INFO_SZ], double *work, double *covar, void *adata); + +/* linear equation & inequality constraints */ +extern int dlevmar_leic_der( + void (*func)(double *p, double *hx, int m, int n, void *adata), + void (*jacf)(double *p, double *j, int m, int n, void *adata), + double *p, double *x, int m, int n, double *A, double *b, int k1, double *C, double *d, int k2, + int itmax, double opts[4], double info[LM_INFO_SZ], double *work, double *covar, void *adata); + +extern int dlevmar_leic_dif( + void (*func)(double *p, double *hx, int m, int n, void *adata), + double *p, double *x, int m, int n, double *A, double *b, int k1, double *C, double *d, int k2, + int itmax, double opts[5], double info[LM_INFO_SZ], double *work, double *covar, void *adata); + +/* linear inequality constraints */ +extern int dlevmar_lic_der( + void (*func)(double *p, double *hx, int m, int n, void *adata), + void (*jacf)(double *p, double *j, int m, int n, void *adata), + double *p, double *x, int m, int n, double *C, double *d, int k2, + int itmax, double opts[4], double info[LM_INFO_SZ], double *work, double *covar, void *adata); + +extern int dlevmar_lic_dif( + void (*func)(double *p, double *hx, int m, int n, void *adata), + double *p, double *x, int m, int n, double *C, double *d, int k2, + int itmax, double opts[5], double info[LM_INFO_SZ], double *work, double *covar, void *adata); #endif /* HAVE_LAPACK */ #endif /* LM_DBL_PREC */ @@ -191,6 +247,56 @@ extern int slevmar_blec_dif( void (*func)(float *p, float *hx, int m, int n, void *adata), float *p, float *x, int m, int n, float *lb, float *ub, float *A, float *b, int k, float *wghts, int itmax, float *opts, float *info, float *work, float *covar, void *adata); + +/* box, linear equations & inequalities constrained minimization */ +extern int slevmar_bleic_der( + void (*func)(float *p, float *hx, int m, int n, void *adata), + void (*jacf)(float *p, float *j, int m, int n, void *adata), + float *p, float *x, int m, int n, float *lb, float *ub, + float *A, float *b, int k1, float *C, float *d, int k2, + int itmax, float *opts, float *info, float *work, float *covar, void *adata); + +extern int slevmar_bleic_dif( + void (*func)(float *p, float *hx, int m, int n, void *adata), + float *p, float *x, int m, int n, float *lb, float *ub, + float *A, float *b, int k1, float *C, float *d, int k2, + int itmax, float *opts, float *info, float *work, float *covar, void *adata); + +/* box & linear inequality constraints */ +extern int slevmar_blic_der( + void (*func)(float *p, float *hx, int m, int n, void *adata), + void (*jacf)(float *p, float *j, int m, int n, void *adata), + float *p, float *x, int m, int n, float *lb, float *ub, float *C, float *d, int k2, + int itmax, float opts[4], float info[LM_INFO_SZ], float *work, float *covar, void *adata); + +extern int slevmar_blic_dif( + void (*func)(float *p, float *hx, int m, int n, void *adata), + float *p, float *x, int m, int n, float *lb, float *ub, float *C, float *d, int k2, + int itmax, float opts[5], float info[LM_INFO_SZ], float *work, float *covar, void *adata); + +/* linear equality & inequality constraints */ +extern int slevmar_leic_der( + void (*func)(float *p, float *hx, int m, int n, void *adata), + void (*jacf)(float *p, float *j, int m, int n, void *adata), + float *p, float *x, int m, int n, float *A, float *b, int k1, float *C, float *d, int k2, + int itmax, float opts[4], float info[LM_INFO_SZ], float *work, float *covar, void *adata); + +extern int slevmar_leic_dif( + void (*func)(float *p, float *hx, int m, int n, void *adata), + float *p, float *x, int m, int n, float *A, float *b, int k1, float *C, float *d, int k2, + int itmax, float opts[5], float info[LM_INFO_SZ], float *work, float *covar, void *adata); + +/* linear inequality constraints */ +extern int slevmar_lic_der( + void (*func)(float *p, float *hx, int m, int n, void *adata), + void (*jacf)(float *p, float *j, int m, int n, void *adata), + float *p, float *x, int m, int n, float *C, float *d, int k2, + int itmax, float opts[4], float info[LM_INFO_SZ], float *work, float *covar, void *adata); + +extern int slevmar_lic_dif( + void (*func)(float *p, float *hx, int m, int n, void *adata), + float *p, float *x, int m, int n, float *C, float *d, int k2, + int itmax, float opts[5], float info[LM_INFO_SZ], float *work, float *covar, void *adata); #endif /* HAVE_LAPACK */ #endif /* LM_SNGL_PREC */ @@ -204,6 +310,7 @@ extern int dAx_eq_b_QRLS(double *A, double *B, double *x, int m, int n); extern int dAx_eq_b_Chol(double *A, double *B, double *x, int m); extern int dAx_eq_b_LU(double *A, double *B, double *x, int m); extern int dAx_eq_b_SVD(double *A, double *B, double *x, int m); +extern int dAx_eq_b_BK(double *A, double *B, double *x, int m); #endif /* LM_DBL_PREC */ #ifdef LM_SNGL_PREC @@ -212,6 +319,7 @@ extern int sAx_eq_b_QRLS(float *A, float *B, float *x, int m, int n); extern int sAx_eq_b_Chol(float *A, float *B, float *x, int m); extern int sAx_eq_b_LU(float *A, float *B, float *x, int m); extern int sAx_eq_b_SVD(float *A, float *B, float *x, int m); +extern int sAx_eq_b_BK(float *A, float *B, float *x, int m); #endif /* LM_SNGL_PREC */ #else /* no LAPACK */ @@ -259,4 +367,4 @@ extern float slevmar_R2(void (*func)(float *p, float *hx, int m, int n, void *ad } #endif -#endif /* _LM_H_ */ +#endif /* _LEVMAR_H_ */ diff --git a/src/levmar/lm_core.c b/src/levmar/lm_core.c index 2a47884..e86d3c9 100644 --- a/src/levmar/lm_core.c +++ b/src/levmar/lm_core.c @@ -37,6 +37,7 @@ #define AX_EQ_B_QR LM_ADD_PREFIX(Ax_eq_b_QR) #define AX_EQ_B_QRLS LM_ADD_PREFIX(Ax_eq_b_QRLS) #define AX_EQ_B_SVD LM_ADD_PREFIX(Ax_eq_b_SVD) +#define AX_EQ_B_BK LM_ADD_PREFIX(Ax_eq_b_BK) #else #define AX_EQ_B_LU LM_ADD_PREFIX(Ax_eq_b_LU_noLapack) #endif /* HAVE_LAPACK */ @@ -293,11 +294,12 @@ if(!(k%100)){ /* solve augmented equations */ #ifdef HAVE_LAPACK - /* 5 alternatives are available: LU, Cholesky, 2 variants of QR decomposition and SVD. + /* 6 alternatives are available: LU, Cholesky, 2 variants of QR decomposition, SVD and LDLt. * Cholesky is the fastest but might be inaccurate; QR is slower but more accurate; * SVD is the slowest but most accurate; LU offers a tradeoff between accuracy and speed */ + //issolved=AX_EQ_B_BK(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_BK; issolved=AX_EQ_B_LU(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_LU; //issolved=AX_EQ_B_CHOL(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_CHOL; //issolved=AX_EQ_B_QR(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_QR; @@ -516,7 +518,7 @@ int (*linsolver)(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)=NULL; if(!work){ worksz=LM_DIF_WORKSZ(m, n); //4*n+4*m + n*m + m*m; - work=(LM_REAL *)malloc(worksz*sizeof(LM_REAL)); /* allocate a big chunk in one step */ + work=(LM_REAL *)calloc(1,worksz*sizeof(LM_REAL)); /* allocate a big chunk in one step */ if(!work){ fprintf(stderr, LCAT(LEVMAR_DIF, "(): memory allocation request failed\n")); return LM_ERROR; @@ -683,11 +685,12 @@ if(!(k%100)){ /* solve augmented equations */ #ifdef HAVE_LAPACK - /* 5 alternatives are available: LU, Cholesky, 2 variants of QR decomposition and SVD. + /* 6 alternatives are available: LU, Cholesky, 2 variants of QR decomposition, SVD and LDLt. * Cholesky is the fastest but might be inaccurate; QR is slower but more accurate; * SVD is the slowest but most accurate; LU offers a tradeoff between accuracy and speed */ + //issolved=AX_EQ_B_BK(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_BK; issolved=AX_EQ_B_LU(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_LU; //issolved=AX_EQ_B_CHOL(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_CHOL; //issolved=AX_EQ_B_QR(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_QR; @@ -837,3 +840,4 @@ if(!(k%100)){ #undef AX_EQ_B_QR #undef AX_EQ_B_QRLS #undef AX_EQ_B_SVD +#undef AX_EQ_B_BK diff --git a/src/levmar/lmbc.c b/src/levmar/lmbc.c index 6495d34..17c4a24 100644 --- a/src/levmar/lmbc.c +++ b/src/levmar/lmbc.c @@ -28,7 +28,7 @@ #include #include -#include "lm.h" +#include "levmar.h" #include "compiler.h" #include "misc.h" diff --git a/src/levmar/lmbc_core.c b/src/levmar/lmbc_core.c index 20600dc..36ac57e 100644 --- a/src/levmar/lmbc_core.c +++ b/src/levmar/lmbc_core.c @@ -44,6 +44,7 @@ #define AX_EQ_B_QR LM_ADD_PREFIX(Ax_eq_b_QR) #define AX_EQ_B_QRLS LM_ADD_PREFIX(Ax_eq_b_QRLS) #define AX_EQ_B_SVD LM_ADD_PREFIX(Ax_eq_b_SVD) +#define AX_EQ_B_BK LM_ADD_PREFIX(Ax_eq_b_BK) #else #define AX_EQ_B_LU LM_ADD_PREFIX(Ax_eq_b_LU_noLapack) #endif /* HAVE_LAPACK */ @@ -551,11 +552,12 @@ if(!(k%100)){ /* solve augmented equations */ #ifdef HAVE_LAPACK - /* 5 alternatives are available: LU, Cholesky, 2 variants of QR decomposition and SVD. + /* 6 alternatives are available: LU, Cholesky, 2 variants of QR decomposition, SVD and LDLt. * Cholesky is the fastest but might be inaccurate; QR is slower but more accurate; * SVD is the slowest but most accurate; LU offers a tradeoff between accuracy and speed */ + //issolved=AX_EQ_B_BK(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_BK; issolved=AX_EQ_B_LU(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_LU; //issolved=AX_EQ_B_CHOL(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_CHOL; //issolved=AX_EQ_B_QR(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_QR; @@ -943,3 +945,4 @@ int ret; #undef AX_EQ_B_QR #undef AX_EQ_B_QRLS #undef AX_EQ_B_SVD +#undef AX_EQ_B_BK diff --git a/src/levmar/lmblec.c b/src/levmar/lmblec.c index 4d41309..2a92e29 100644 --- a/src/levmar/lmblec.c +++ b/src/levmar/lmblec.c @@ -28,17 +28,15 @@ #include #include -#include "lm.h" +#include "levmar.h" #include "misc.h" #ifndef HAVE_LAPACK #ifdef _MSC_VER #pragma message("Combined box and linearly constrained optimization requires LAPACK and was not compiled!") -/* #else -#warning Combined box and linearly constrained optimization requires LAPACK and was not compiled! -*/ +//#warning Combined box and linearly constrained optimization requires LAPACK and was not compiled! #endif // _MSC_VER #else // LAPACK present diff --git a/src/levmar/lmbleic.c b/src/levmar/lmbleic.c new file mode 100644 index 0000000..7deaf99 --- /dev/null +++ b/src/levmar/lmbleic.c @@ -0,0 +1,89 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Levenberg - Marquardt non-linear minimization algorithm +// Copyright (C) 2009 Manolis Lourakis (lourakis at ics forth gr) +// Institute of Computer Science, Foundation for Research & Technology - Hellas +// Heraklion, Crete, Greece. +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +/******************************************************************************* + * Wrappers for linear inequality constrained Levenberg-Marquardt minimization. + * The same core code is used with appropriate #defines to derive single and + * double precision versions, see also lmbleic_core.c + *******************************************************************************/ + +#include +#include +#include +#include + +#include "levmar.h" +#include "misc.h" + + +#ifndef HAVE_LAPACK + +#ifdef _MSC_VER +#pragma message("Linear inequalities constrained optimization requires LAPACK and was not compiled!") +#else +//#warning Linear inequalities constrained optimization requires LAPACK and was not compiled! +#endif // _MSC_VER + +#else // LAPACK present + +#if !defined(LM_DBL_PREC) && !defined(LM_SNGL_PREC) +#error At least one of LM_DBL_PREC, LM_SNGL_PREC should be defined! +#endif + + +#ifdef LM_SNGL_PREC +/* single precision (float) definitions */ +#define LM_REAL float +#define LM_PREFIX s + +#define LM_REAL_MAX FLT_MAX +#define LM_REAL_MIN -FLT_MAX +#define __SUBCNST(x) x##F +#define LM_CNST(x) __SUBCNST(x) // force substitution + +#include "lmbleic_core.c" // read in core code + +#undef LM_REAL +#undef LM_PREFIX +#undef LM_REAL_MAX +#undef LM_REAL_MIN +#undef __SUBCNST +#undef LM_CNST +#endif /* LM_SNGL_PREC */ + +#ifdef LM_DBL_PREC +/* double precision definitions */ +#define LM_REAL double +#define LM_PREFIX d + +#define LM_REAL_MAX DBL_MAX +#define LM_REAL_MIN -DBL_MAX +#define LM_CNST(x) (x) + +#include "lmbleic_core.c" // read in core code + +#undef LM_REAL +#undef LM_PREFIX +#undef LM_REAL_MAX +#undef LM_REAL_MIN +#undef LM_CNST +#endif /* LM_DBL_PREC */ + +#endif /* HAVE_LAPACK */ + diff --git a/src/levmar/lmbleic_core.c b/src/levmar/lmbleic_core.c new file mode 100644 index 0000000..7ff1d34 --- /dev/null +++ b/src/levmar/lmbleic_core.c @@ -0,0 +1,506 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Levenberg - Marquardt non-linear minimization algorithm +// Copyright (C) 2009 Manolis Lourakis (lourakis at ics forth gr) +// Institute of Computer Science, Foundation for Research & Technology - Hellas +// Heraklion, Crete, Greece. +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#ifndef LM_REAL // not included by lmbleic.c +#error This file should not be compiled directly! +#endif + + +/* precision-specific definitions */ +#define LMBLEIC_DATA LM_ADD_PREFIX(lmbleic_data) +#define LMBLEIC_ELIM LM_ADD_PREFIX(lmbleic_elim) +#define LMBLEIC_FUNC LM_ADD_PREFIX(lmbleic_func) +#define LMBLEIC_JACF LM_ADD_PREFIX(lmbleic_jacf) +#define LEVMAR_BLEIC_DER LM_ADD_PREFIX(levmar_bleic_der) +#define LEVMAR_BLEIC_DIF LM_ADD_PREFIX(levmar_bleic_dif) +#define LEVMAR_BLIC_DER LM_ADD_PREFIX(levmar_blic_der) +#define LEVMAR_BLIC_DIF LM_ADD_PREFIX(levmar_blic_dif) +#define LEVMAR_LEIC_DER LM_ADD_PREFIX(levmar_leic_der) +#define LEVMAR_LEIC_DIF LM_ADD_PREFIX(levmar_leic_dif) +#define LEVMAR_LIC_DER LM_ADD_PREFIX(levmar_lic_der) +#define LEVMAR_LIC_DIF LM_ADD_PREFIX(levmar_lic_dif) +#define LEVMAR_BLEC_DER LM_ADD_PREFIX(levmar_blec_der) +#define LEVMAR_BLEC_DIF LM_ADD_PREFIX(levmar_blec_dif) +#define LEVMAR_TRANS_MAT_MAT_MULT LM_ADD_PREFIX(levmar_trans_mat_mat_mult) +#define LEVMAR_COVAR LM_ADD_PREFIX(levmar_covar) +#define LEVMAR_FDIF_FORW_JAC_APPROX LM_ADD_PREFIX(levmar_fdif_forw_jac_approx) + +struct LMBLEIC_DATA{ + LM_REAL *jac; + int nineqcnstr; // #inequality constraints + void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata); + void (*jacf)(LM_REAL *p, LM_REAL *jac, int m, int n, void *adata); + void *adata; +}; + + +/* wrapper ensuring that the user-supplied function is called with the right number of variables (i.e. m) */ +static void LMBLEIC_FUNC(LM_REAL *pext, LM_REAL *hx, int mm, int n, void *adata) +{ +struct LMBLEIC_DATA *data=(struct LMBLEIC_DATA *)adata; +int m; + + m=mm-data->nineqcnstr; + (*(data->func))(pext, hx, m, n, data->adata); +} + + +/* wrapper for computing the Jacobian at pext. The Jacobian is nxmm */ +static void LMBLEIC_JACF(LM_REAL *pext, LM_REAL *jacext, int mm, int n, void *adata) +{ +struct LMBLEIC_DATA *data=(struct LMBLEIC_DATA *)adata; +int m; +register int i, j; +LM_REAL *jac, *jacim, *jacextimm; + + m=mm-data->nineqcnstr; + jac=data->jac; + + (*(data->jacf))(pext, jac, m, n, data->adata); + + for(i=0; i=d, C being k2xm, d k2x1. + * + * The inequalities are converted to equations by introducing surplus variables, + * i.e. c^T*p >= d becomes c^T*p - y = d, with y>=0. To transform all inequalities + * to equations, a total of k2 surplus variables are introduced; a problem with only + * box and linear constraints results then and is solved with LEVMAR_BLEC_DER() + * Note that opposite direction inequalities should be converted to the desired + * direction by negating, i.e. c^T*p <= d becomes -c^T*p >= -d + * + * This function requires an analytic Jacobian. In case the latter is unavailable, + * use LEVMAR_BLEIC_DIF() bellow + * + */ +int LEVMAR_BLEIC_DER( + void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata), /* functional relation describing measurements. A p \in R^m yields a \hat{x} \in R^n */ + void (*jacf)(LM_REAL *p, LM_REAL *j, int m, int n, void *adata), /* function to evaluate the Jacobian \part x / \part p */ + LM_REAL *p, /* I/O: initial parameter estimates. On output has the estimated solution */ + LM_REAL *x, /* I: measurement vector. NULL implies a zero vector */ + int m, /* I: parameter vector dimension (i.e. #unknowns) */ + int n, /* I: measurement vector dimension */ + LM_REAL *lb, /* I: vector of lower bounds. If NULL, no lower bounds apply */ + LM_REAL *ub, /* I: vector of upper bounds. If NULL, no upper bounds apply */ + LM_REAL *A, /* I: equality constraints matrix, k1xm. If NULL, no linear equation constraints apply */ + LM_REAL *b, /* I: right hand constraints vector, k1x1 */ + int k1, /* I: number of constraints (i.e. A's #rows) */ + LM_REAL *C, /* I: inequality constraints matrix, k2xm */ + LM_REAL *d, /* I: right hand constraints vector, k2x1 */ + int k2, /* I: number of inequality constraints (i.e. C's #rows) */ + int itmax, /* I: maximum number of iterations */ + LM_REAL opts[4], /* I: minim. options [\mu, \epsilon1, \epsilon2, \epsilon3]. Respectively the scale factor for initial \mu, + * stopping thresholds for ||J^T e||_inf, ||Dp||_2 and ||e||_2. Set to NULL for defaults to be used + */ + LM_REAL info[LM_INFO_SZ], + /* O: information regarding the minimization. Set to NULL if don't care + * info[0]= ||e||_2 at initial p. + * info[1-4]=[ ||e||_2, ||J^T e||_inf, ||Dp||_2, mu/max[J^T J]_ii ], all computed at estimated p. + * info[5]= # iterations, + * info[6]=reason for terminating: 1 - stopped by small gradient J^T e + * 2 - stopped by small Dp + * 3 - stopped by itmax + * 4 - singular matrix. Restart from current p with increased mu + * 5 - no further error reduction is possible. Restart with increased mu + * 6 - stopped by small ||e||_2 + * 7 - stopped by invalid (i.e. NaN or Inf) "func" values. This is a user error + * info[7]= # function evaluations + * info[8]= # Jacobian evaluations + * info[9]= # linear systems solved, i.e. # attempts for reducing error + */ + LM_REAL *work, /* working memory at least LM_BLEIC_DER_WORKSZ() reals large, allocated if NULL */ + LM_REAL *covar, /* O: Covariance matrix corresponding to LS solution; mxm. Set to NULL if not needed. */ + void *adata) /* pointer to possibly additional data, passed uninterpreted to func & jacf. + * Set to NULL if not needed + */ +{ + struct LMBLEIC_DATA data; + LM_REAL *ptr, *pext, *Aext, *bext, *covext; /* corresponding to p, A, b, covar for the full set of variables; + pext=[p, surplus], pext is mm, Aext is (k1+k2)xmm, bext (k1+k2), covext is mmxmm + */ + LM_REAL *lbext, *ubext; // corresponding to lb, ub for the full set of variables + int mm, ret, k12; + register int i, j, ii; + register LM_REAL tmp; + LM_REAL locinfo[LM_INFO_SZ]; + + if(!jacf){ + fprintf(stderr, RCAT("No function specified for computing the Jacobian in ", LEVMAR_BLEIC_DER) + RCAT("().\nIf no such function is available, use ", LEVMAR_BLEIC_DIF) RCAT("() rather than ", LEVMAR_BLEIC_DER) "()\n"); + return LM_ERROR; + } + + if(!C || !d){ + fprintf(stderr, RCAT(LCAT(LEVMAR_BLEIC_DER, "(): missing inequality constraints, use "), LEVMAR_BLEC_DER) "() in this case!\n"); + return LM_ERROR; + } + + if(!A || !b) k1=0; // sanity check + + mm=m+k2; + + if(n=0 */ + lbext[j]=0.0; + ubext[j]=LM_REAL_MAX; + } + /* set the first m elements of pext equal to p */ + for(i=0; i=0 */ + lbext[j]=0.0; + ubext[j]=LM_REAL_MAX; + } + /* set the first m elements of pext equal to p */ + for(i=0; i #include -#include "lm.h" +#include "levmar.h" #ifndef LM_DBL_PREC #error Demo program assumes that levmar has been compiled with double precision, see LM_DBL_PREC! @@ -162,6 +162,36 @@ double ui, tmp; } } +/* Osborne's problem, minimum at (0.3754, 1.9358, -1.4647, 0.0129, 0.0221) */ +void osborne(double *p, double *x, int m, int n, void *data) +{ +register int i; +double t; + + for(i=0; i= -1.0; + * constr2: p[2] + p[3] - 2*p[4] >= -2.0; + * constr3: p[1] - p[4] <= 7.0; + * + * + */ +void mod2hs52(double *p, double *x, int m, int n, void *data) +{ + x[0]=4.0*p[0]-p[1]; + x[1]=p[1]+p[2]-2.0; + x[2]=p[3]-1.0; + x[3]=p[4]-1.0; + x[4]=p[0]-0.5; +} + +void jacmod2hs52(double *p, double *jac, int m, int n, void *data) +{ +register int j=0; + + jac[j++]=4.0; + jac[j++]=-1.0; + jac[j++]=0.0; + jac[j++]=0.0; + jac[j++]=0.0; + + jac[j++]=0.0; + jac[j++]=1.0; + jac[j++]=1.0; + jac[j++]=0.0; + jac[j++]=0.0; + + jac[j++]=0.0; + jac[j++]=0.0; + jac[j++]=0.0; + jac[j++]=1.0; + jac[j++]=0.0; + + jac[j++]=0.0; + jac[j++]=0.0; + jac[j++]=0.0; + jac[j++]=0.0; + jac[j++]=1.0; + + jac[j++]=1.0; + jac[j++]=0.0; + jac[j++]=0.0; + jac[j++]=0.0; + jac[j++]=0.0; +} + /* Schittkowski (modified) problem 235 (box/linearly constrained), minimum at (-1.725, 2.9, 0.725) * constr1: p[0] + p[2] = -1.0; * @@ -653,13 +737,55 @@ register int j=0; jac[j+3]=R9*p[1]+2*p[3]; } +/* Hock - Schittkowski (modified) problem 76 (linear inequalities & equations constrained), minimum at (0.0, 0.00909091, 0.372727, 0.354545) + * The non-squared terms in the objective function have been removed, the rhs of constr2 has been changed to 0.4 (from 4) + * and constr3 has been changed to an equation. + * + * constr1: p[0] + 2*p[1] + p[2] + p[3] <= 5; + * constr2: 3*p[0] + p[1] + 2*p[2] - p[3] <= 0.4; + * constr3: p[1] + 4*p[2] = 1.5; + * + */ +void modhs76(double *p, double *x, int m, int n, void *data) +{ + x[0]=p[0]; + x[1]=sqrt(0.5)*p[1]; + x[2]=p[2]; + x[3]=sqrt(0.5)*p[3]; +} + +void jacmodhs76(double *p, double *jac, int m, int n, void *data) +{ +register int j=0; + + jac[j++]=1.0; + jac[j++]=0.0; + jac[j++]=0.0; + jac[j++]=0.0; + + jac[j++]=0.0; + jac[j++]=sqrt(0.5); + jac[j++]=0.0; + jac[j++]=0.0; + + jac[j++]=0.0; + jac[j++]=0.0; + jac[j++]=1.0; + jac[j++]=0.0; + + jac[j++]=0.0; + jac[j++]=0.0; + jac[j++]=0.0; + jac[j++]=sqrt(0.5); +} + int main() { register int i, j; int problem, ret; -double p[5], // 6 is max(2, 3, 5) +double p[5], // 5 is max(2, 3, 5) x[16]; // 16 is max(2, 3, 5, 6, 16) int m, n; double opts[LM_OPTS_SZ], info[LM_INFO_SZ]; @@ -669,6 +795,7 @@ char *probname[]={ "Powell's function", "Wood's function", "Meyer's (reformulated) problem", + "Osborne's problem", "helical valley function", "Boggs & Tolle's problem #3", "Hock - Schittkowski problem #28", @@ -679,13 +806,15 @@ char *probname[]={ "hatfldb problem", "hatfldc problem", "equilibrium combustion problem", - "Hock - Schittkowski modified problem #52", + "Hock - Schittkowski modified #1 problem #52", "Schittkowski modified problem #235", "Boggs & Tolle modified problem #7", + "Hock - Schittkowski modified #2 problem #52", + "Hock - Schittkowski modified problem #76", }; opts[0]=LM_INIT_MU; opts[1]=1E-15; opts[2]=1E-15; opts[3]=1E-20; - opts[4]=LM_DIFF_DELTA; // relevant only if the Jacobian is approximated using finite differences; specifies forward differencing + opts[4]= LM_DIFF_DELTA; // relevant only if the Jacobian is approximated using finite differences; specifies forward differencing //opts[4]=-LM_DIFF_DELTA; // specifies central differencing to approximate Jacobian; more accurate but more expensive to compute! /* uncomment the appropriate line below to select a minimization problem */ @@ -695,12 +824,13 @@ char *probname[]={ //2; // Powell's function //3; // Wood's function 4; // Meyer's (reformulated) problem - //5; // helical valley function + //5; // Osborne's problem + //6; // helical valley function #ifdef HAVE_LAPACK - //6; // Boggs & Tolle's problem 3 - //7; // Hock - Schittkowski problem 28 - //8; // Hock - Schittkowski problem 48 - //9; // Hock - Schittkowski problem 51 + //7; // Boggs & Tolle's problem 3 + //8; // Hock - Schittkowski problem 28 + //9; // Hock - Schittkowski problem 48 + //10; // Hock - Schittkowski problem 51 #else // no LAPACK #ifdef _MSC_VER #pragma message("LAPACK not available, some test problems cannot be used") @@ -709,21 +839,24 @@ char *probname[]={ #endif // _MSC_VER #endif /* HAVE_LAPACK */ - //10; // Hock - Schittkowski problem 01 - //11; // Hock - Schittkowski modified problem 21 - //12; // hatfldb problem - //13; // hatfldc problem - //14; // equilibrium combustion problem + //11; // Hock - Schittkowski problem 01 + //12; // Hock - Schittkowski modified problem 21 + //13; // hatfldb problem + //14; // hatfldc problem + //15; // equilibrium combustion problem #ifdef HAVE_LAPACK - //15; // Hock - Schittkowski modified problem 52 - //16; // Schittkowski modified problem 235 - //17; // Boggs & Tolle modified problem #7 + //16; // Hock - Schittkowski modified #1 problem 52 + //17; // Schittkowski modified problem 235 + //18; // Boggs & Tolle modified problem #7 + //19; // Hock - Schittkowski modified #2 problem 52 + //20; // Hock - Schittkowski modified problem #76" #endif /* HAVE_LAPACK */ switch(problem){ default: fprintf(stderr, "unknown problem specified (#%d)! Note that some minimization problems require LAPACK.\n", problem); exit(1); break; + case 0: /* Rosenbrock function */ m=2; n=2; @@ -798,10 +931,27 @@ char *probname[]={ for(i=0; i #include -#include "lm.h" +#include "levmar.h" #include "misc.h" @@ -35,10 +35,8 @@ #ifdef _MSC_VER #pragma message("Linearly constrained optimization requires LAPACK and was not compiled!") -/* #else -#warning Linearly constrained optimization requires LAPACK and was not compiled! -*/ +//#warning Linearly constrained optimization requires LAPACK and was not compiled! #endif // _MSC_VER #else // LAPACK present diff --git a/src/levmar/matlab/CMakeLists.txt b/src/levmar/matlab/CMakeLists.txt new file mode 100644 index 0000000..3f529f6 --- /dev/null +++ b/src/levmar/matlab/CMakeLists.txt @@ -0,0 +1,58 @@ +# CMake file for levmar's MEX-file; see http://www.cmake.org +# Requires FindMatlab.cmake included with cmake + +PROJECT(LEVMARMEX) +#CMAKE_MINIMUM_REQUIRED(VERSION 1.4) + +INCLUDE("C:/Program Files/CMake 2.4/share/cmake-2.4/Modules/FindMatlab.cmake") + +# f2c is sometimes equivalent to libF77 & libI77; in that case, set HAVE_F2C to 0 +SET(HAVE_F2C 1 CACHE BOOL "Do we have f2c or F77/I77?" ) + +# the directory where the lapack/blas/f2c libraries reside +SET(LAPACKBLAS_DIR /usr/lib CACHE PATH "Path to lapack/blas libraries") + +# the directory where levmar.h resides +SET(LM_H_DIR .. CACHE PATH "Path to levmar.h") +# the directory where the levmar library resides +SET(LEVMAR_DIR .. CACHE PATH "Path to levmar library") + +# actual names for the lapack/blas/f2c libraries +SET(LAPACK_LIB lapack CACHE STRING "The name of the lapack library") +SET(BLAS_LIB blas CACHE STRING "The name of the blas library") +IF(HAVE_F2C) + SET(F2C_LIB f2c CACHE STRING "The name of the f2c library") +ELSE(HAVE_F2C) + SET(F77_LIB libF77 CACHE STRING "The name of the F77 library") + SET(I77_LIB libI77 CACHE STRING "The name of the I77 library") +ENDIF(HAVE_F2C) + +########################## NO CHANGES BEYOND THIS POINT ########################## + +INCLUDE_DIRECTORIES(${LM_H_DIR}) +LINK_DIRECTORIES(${LAPACKBLAS_DIR} ${LEVMAR_DIR}) + +SET(SRC levmar.c) + +# naming conventions for the generated file's suffix +IF(WIN32) + SET(SUFFIX ".mexw32") +ELSE(WIN32) + SET(SUFFIX ".mexglx") +ENDIF(WIN32) + +SET(OUTNAME "levmar${SUFFIX}") + +ADD_LIBRARY(${OUTNAME} MODULE ${SRC}) + +IF(HAVE_F2C) + ADD_CUSTOM_COMMAND(OUTPUT ${OUTNAME} + DEPENDS ${SRC} + COMMAND mex + ARGS -O -I${LM_H_DIR} ${SRC} -I${MATLAB_INCLUDE_DIR} -L${LAPACKBLAS_DIR} -L${LEVMAR_DIR} -L${MATLAB_MEX_LIBRARY} -llevmar -l${LAPACK_LIB} -l${BLAS_LIB} -l${F2C_LIB} -output ${MATLAB_LIBRARIES} ${OUTNAME}) +ELSE(HAVE_F2C) + ADD_CUSTOM_COMMAND(OUTPUT ${OUTNAME} + DEPENDS ${SRC} + COMMAND mex + ARGS -O -I${LM_H_DIR} ${SRC} -I${MATLAB_INCLUDE_DIR} -L${LAPACKBLAS_DIR} -L${LEVMAR_DIR} -L${MATLAB_MEX_LIBRARY} -llevmar -l${LAPACK_LIB} -l${BLAS_LIB} -l${F77_LIB} -l${I77_LIB} ${MATLAB_LIBRARIES} -output ${OUTNAME}) +ENDIF(HAVE_F2C) diff --git a/src/levmar/matlab/README.txt b/src/levmar/matlab/README.txt index c35637e..cdc84c9 100644 --- a/src/levmar/matlab/README.txt +++ b/src/levmar/matlab/README.txt @@ -1,5 +1,6 @@ This directory contains a matlab MEX interface to levmar. This interface has been tested with Matlab v. 6.5 R13 under linux and v. 7.4 R2007 under Windows. +Users have also reported success with GNU Octave. FILES The following files are included: diff --git a/src/levmar/matlab/jacbt3.m b/src/levmar/matlab/jacbt3.m index becb3b5..64383f2 100644 --- a/src/levmar/matlab/jacbt3.m +++ b/src/levmar/matlab/jacbt3.m @@ -1,4 +1,4 @@ -function jac = bt3_jac(p, adata) +function jac = jacbt3(p, adata) n=5; m=5; diff --git a/src/levmar/matlab/jacexpfit.m b/src/levmar/matlab/jacexpfit.m index 13747fd..978c7b7 100644 --- a/src/levmar/matlab/jacexpfit.m +++ b/src/levmar/matlab/jacexpfit.m @@ -1,4 +1,4 @@ -function jac = expfit_jac(p, data) +function jac = jacexpfit(p, data) n=data; m=max(size(p)); diff --git a/src/levmar/matlab/jachs01.m b/src/levmar/matlab/jachs01.m index 49ee83a..e08e5e3 100644 --- a/src/levmar/matlab/jachs01.m +++ b/src/levmar/matlab/jachs01.m @@ -1,4 +1,4 @@ -function jac = hs01_jac(p) +function jac = jachs01(p) m=2; jac(1, 1:m)=[-20.0*p(1), 10.0]; diff --git a/src/levmar/matlab/jacmeyer.m b/src/levmar/matlab/jacmeyer.m index 9d8ff5d..475d361 100644 --- a/src/levmar/matlab/jacmeyer.m +++ b/src/levmar/matlab/jacmeyer.m @@ -1,4 +1,4 @@ -function jac = meyer_jac(p, data1, data2) +function jac = jacmeyer(p, data1, data2) n=16; m=3; diff --git a/src/levmar/matlab/jacmodhs52.m b/src/levmar/matlab/jacmodhs52.m index 30b524b..2be4924 100644 --- a/src/levmar/matlab/jacmodhs52.m +++ b/src/levmar/matlab/jacmodhs52.m @@ -1,4 +1,4 @@ -function jac = modhs52_jac(p) +function jac = jacmodhs52(p) m=5; jac(1, 1:m)=[4.0, -1.0, 0.0, 0.0, 0.0]; diff --git a/src/levmar/matlab/jacmodhs76.m b/src/levmar/matlab/jacmodhs76.m new file mode 100644 index 0000000..5ffb3dc --- /dev/null +++ b/src/levmar/matlab/jacmodhs76.m @@ -0,0 +1,7 @@ +function jac = jacmodhs76(p) + m=4; + + jac(1, 1:m)=[1.0, 0.0, 0.0, 0.0]; + jac(2, 1:m)=[0.0, sqrt(0.5), 0.0, 0.0]; + jac(3, 1:m)=[0.0, 0.0, 1.0, 0.0]; + jac(4, 1:m)=[0.0, 0.0, 0.0, sqrt(0.5)]; diff --git a/src/levmar/matlab/jacosborne.m b/src/levmar/matlab/jacosborne.m new file mode 100644 index 0000000..22082c6 --- /dev/null +++ b/src/levmar/matlab/jacosborne.m @@ -0,0 +1,11 @@ +function jac = jacosborne(p) + n=33; + m=5; + + for i=1:n + t=10*(i-1); + tmp1=exp(-p(4)*t); + tmp2=exp(-p(5)*t); + + jac(i, 1:m)=[1.0, tmp1, tmp2, -p(2)*t*tmp1, -p(3)*t*tmp2]; + end diff --git a/src/levmar/matlab/levmar.c b/src/levmar/matlab/levmar.c index 77129ef..b2fb6a2 100644 --- a/src/levmar/matlab/levmar.c +++ b/src/levmar/matlab/levmar.c @@ -24,7 +24,7 @@ #include #include -#include +#include #include @@ -46,6 +46,10 @@ #define MIN_CONSTRAINED_BC 1 #define MIN_CONSTRAINED_LEC 2 #define MIN_CONSTRAINED_BLEC 3 +#define MIN_CONSTRAINED_BLEIC 4 +#define MIN_CONSTRAINED_BLIC 5 +#define MIN_CONSTRAINED_LEIC 6 +#define MIN_CONSTRAINED_LIC 7 struct mexdata { /* matlab names of the fitting function & its Jacobian */ @@ -68,9 +72,9 @@ static void matlabFmtdErrMsgTxt(char *fmt, ...) char buf[256]; va_list args; - va_start(args, fmt); - vsprintf(buf, fmt, args); - va_end(args); + va_start(args, fmt); + vsprintf(buf, fmt, args); + va_end(args); mexErrMsgTxt(buf); } @@ -81,9 +85,9 @@ static void matlabFmtdWarnMsgTxt(char *fmt, ...) char buf[256]; va_list args; - va_start(args, fmt); - vsprintf(buf, fmt, args); - va_end(args); + va_start(args, fmt); + vsprintf(buf, fmt, args); + va_end(args); mexWarnMsgTxt(buf); } @@ -183,7 +187,7 @@ double *mp; else if(!mxIsDouble(lhs[0]) || mxIsComplex(lhs[0]) || !(mxGetM(lhs[0])==1 || mxGetN(lhs[0])==1) || __MAX__(mxGetM(lhs[0]), mxGetN(lhs[0]))!=n){ fprintf(stderr, "levmar: '%s' should produce a real vector with %d elements (got %d).\n", - dat->fname, m, __MAX__(mxGetM(lhs[0]), mxGetN(lhs[0]))); + dat->fname, n, __MAX__(mxGetM(lhs[0]), mxGetN(lhs[0]))); ret=1; } /* delete the matrix created by matlab */ @@ -220,20 +224,26 @@ double *mp; [ret, p, info, covar]=levmar_bc (f, j, p0, x, itmax, opts, 'bc', lb, ub, ...) [ret, p, info, covar]=levmar_lec (f, j, p0, x, itmax, opts, 'lec', A, b, ...) [ret, p, info, covar]=levmar_blec(f, j, p0, x, itmax, opts, 'blec', lb, ub, A, b, wghts, ...) + +[ret, p, info, covar]=levmar_bleic(f, j, p0, x, itmax, opts, 'bleic', lb, ub, A, b, C, d, ...) +[ret, p, info, covar]=levmar_blic (f, j, p0, x, itmax, opts, 'blic', lb, ub, C, d, ...) +[ret, p, info, covar]=levmar_leic (f, j, p0, x, itmax, opts, 'leic', A, b, C, d, ...) +[ret, p, info, covar]=levmar_lic (f, j, p0, x, itmax, opts, 'lic', C, d, ...) + */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *Prhs[]) { register int i; register double *pdbl; -mxArray **prhs=(mxArray **)&Prhs[0], *At; +mxArray **prhs=(mxArray **)&Prhs[0], *At, *Ct; struct mexdata mdata; int len, status; double *p, *p0, *ret, *x; -int m, n, havejac, Arows, itmax, nopts, mintype, nextra; +int m, n, havejac, Arows, Crows, itmax, nopts, mintype, nextra; double opts[LM_OPTS_SZ]={LM_INIT_MU, LM_STOP_THRESH, LM_STOP_THRESH, LM_STOP_THRESH, LM_DIFF_DELTA}; double info[LM_INFO_SZ]; -double *lb=NULL, *ub=NULL, *A=NULL, *b=NULL, *wghts=NULL, *covar=NULL; +double *lb=NULL, *ub=NULL, *A=NULL, *b=NULL, *wghts=NULL, *C=NULL, *d=NULL, *covar=NULL; /* parse input args; start by checking their number */ if((nrhs<5)) @@ -362,6 +372,10 @@ if(!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) || !(mxGetM(prhs[1])==1 && mxGet else if(!strncmp(minhowto, "bc", 2)) mintype=MIN_CONSTRAINED_BC; else if(!strncmp(minhowto, "lec", 3)) mintype=MIN_CONSTRAINED_LEC; else if(!strncmp(minhowto, "blec", 4)) mintype=MIN_CONSTRAINED_BLEC; + else if(!strncmp(minhowto, "bleic", 5)) mintype=MIN_CONSTRAINED_BLEIC; + else if(!strncmp(minhowto, "blic", 4)) mintype=MIN_CONSTRAINED_BLIC; + else if(!strncmp(minhowto, "leic", 4)) mintype=MIN_CONSTRAINED_LEIC; + else if(!strncmp(minhowto, "lic", 3)) mintype=MIN_CONSTRAINED_BLIC; else matlabFmtdErrMsgTxt("levmar: unknown minimization type '%s'.", minhowto); mxFree(minhowto); @@ -378,7 +392,7 @@ if(!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) || !(mxGetM(prhs[1])==1 && mxGet * upon the minimization type determined above */ /** lb, ub **/ - if(nrhs>=7 && (mintype==MIN_CONSTRAINED_BC || mintype==MIN_CONSTRAINED_BLEC)){ + if(nrhs>=7 && (mintype==MIN_CONSTRAINED_BC || mintype==MIN_CONSTRAINED_BLEC || mintype==MIN_CONSTRAINED_BLIC || mintype==MIN_CONSTRAINED_BLEIC)){ /* check if the next two arguments are real row or column vectors */ if(mxIsDouble(prhs[5]) && !mxIsComplex(prhs[5]) && (mxGetM(prhs[5])==1 || mxGetN(prhs[5])==1)){ if(mxIsDouble(prhs[6]) && !mxIsComplex(prhs[6]) && (mxGetM(prhs[6])==1 || mxGetN(prhs[6])==1)){ @@ -397,7 +411,7 @@ if(!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) || !(mxGetM(prhs[1])==1 && mxGet } /** A, b **/ - if(nrhs>=7 && (mintype==MIN_CONSTRAINED_LEC || mintype==MIN_CONSTRAINED_BLEC)){ + if(nrhs>=7 && (mintype==MIN_CONSTRAINED_LEC || mintype==MIN_CONSTRAINED_BLEC || mintype==MIN_CONSTRAINED_LEIC || mintype==MIN_CONSTRAINED_BLEIC)){ /* check if the next two arguments are a real matrix and a real row or column vector */ if(mxIsDouble(prhs[5]) && !mxIsComplex(prhs[5]) && mxGetM(prhs[5])>=1 && mxGetN(prhs[5])>=1){ if(mxIsDouble(prhs[6]) && !mxIsComplex(prhs[6]) && (mxGetM(prhs[6])==1 || mxGetN(prhs[6])==1)){ @@ -428,6 +442,27 @@ if(!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) || !(mxGetM(prhs[1])==1 && mxGet } } } + + /** C, d **/ + if(nrhs>=7 && (mintype==MIN_CONSTRAINED_BLEIC || mintype==MIN_CONSTRAINED_BLIC || mintype==MIN_CONSTRAINED_LEIC || mintype==MIN_CONSTRAINED_LIC)){ + /* check if the next two arguments are a real matrix and a real row or column vector */ + if(mxIsDouble(prhs[5]) && !mxIsComplex(prhs[5]) && mxGetM(prhs[5])>=1 && mxGetN(prhs[5])>=1){ + if(mxIsDouble(prhs[6]) && !mxIsComplex(prhs[6]) && (mxGetM(prhs[6])==1 || mxGetN(prhs[6])==1)){ + if((i=mxGetN(prhs[5]))!=m) + matlabFmtdErrMsgTxt("levmar: C must have %d columns, got %d.", m, i); + if((i=__MAX__(mxGetM(prhs[6]), mxGetN(prhs[6])))!=(Crows=mxGetM(prhs[5]))) + matlabFmtdErrMsgTxt("levmar: d must have %d elements, got %d.", Crows, i); + + Ct=prhs[5]; + d=mxGetPr(prhs[6]); + C=getTranspose(Ct); + + prhs+=2; + nrhs-=2; + } + } + } + /* arguments below this point are assumed to be extra arguments passed * to every invocation of the fitting function and its Jacobian */ @@ -471,8 +506,8 @@ extraargs: covar=mxMalloc(m*m*sizeof(double)); /* invoke levmar */ - if(!lb && !ub){ - if(!A && !b){ /* no constraints */ + switch(mintype){ + case MIN_UNCONSTRAINED: /* no constraints */ if(havejac) status=dlevmar_der(func, jacfunc, p, x, m, n, itmax, opts, info, NULL, covar, (void *)&mdata); else @@ -481,8 +516,18 @@ extraargs: fflush(stderr); fprintf(stderr, "LEVMAR: calling dlevmar_der()/dlevmar_dif()\n"); #endif /* DEBUG */ - } - else{ /* linear constraints */ + break; + case MIN_CONSTRAINED_BC: /* box constraints */ + if(havejac) + status=dlevmar_bc_der(func, jacfunc, p, x, m, n, lb, ub, itmax, opts, info, NULL, covar, (void *)&mdata); + else + status=dlevmar_bc_dif(func, p, x, m, n, lb, ub, itmax, opts, info, NULL, covar, (void *)&mdata); +#ifdef DEBUG + fflush(stderr); + fprintf(stderr, "LEVMAR: calling dlevmar_bc_der()/dlevmar_bc_dif()\n"); +#endif /* DEBUG */ + break; + case MIN_CONSTRAINED_LEC: /* linear equation constraints */ #ifdef HAVE_LAPACK if(havejac) status=dlevmar_lec_der(func, jacfunc, p, x, m, n, A, b, Arows, itmax, opts, info, NULL, covar, (void *)&mdata); @@ -496,35 +541,86 @@ extraargs: fflush(stderr); fprintf(stderr, "LEVMAR: calling dlevmar_lec_der()/dlevmar_lec_dif()\n"); #endif /* DEBUG */ - } - } - else{ - if(!A && !b){ /* box constraints */ + break; + case MIN_CONSTRAINED_BLEC: /* box & linear equation constraints */ +#ifdef HAVE_LAPACK if(havejac) - status=dlevmar_bc_der(func, jacfunc, p, x, m, n, lb, ub, itmax, opts, info, NULL, covar, (void *)&mdata); + status=dlevmar_blec_der(func, jacfunc, p, x, m, n, lb, ub, A, b, Arows, wghts, itmax, opts, info, NULL, covar, (void *)&mdata); else - status=dlevmar_bc_dif(func, p, x, m, n, lb, ub, itmax, opts, info, NULL, covar, (void *)&mdata); + status=dlevmar_blec_dif(func, p, x, m, n, lb, ub, A, b, Arows, wghts, itmax, opts, info, NULL, covar, (void *)&mdata); +#else + mexErrMsgTxt("levmar: no box & linear constraints support, HAVE_LAPACK was not defined during MEX-file compilation."); +#endif /* HAVE_LAPACK */ + #ifdef DEBUG fflush(stderr); - fprintf(stderr, "LEVMAR: calling dlevmar_bc_der()/dlevmar_bc_dif()\n"); + fprintf(stderr, "LEVMAR: calling dlevmar_blec_der()/dlevmar_blec_dif()\n"); #endif /* DEBUG */ - } - else{ /* box & linear constraints */ + break; + case MIN_CONSTRAINED_BLEIC: /* box, linear equation & inequalities constraints */ #ifdef HAVE_LAPACK if(havejac) - status=dlevmar_blec_der(func, jacfunc, p, x, m, n, lb, ub, A, b, Arows, wghts, itmax, opts, info, NULL, covar, (void *)&mdata); + status=dlevmar_bleic_der(func, jacfunc, p, x, m, n, lb, ub, A, b, Arows, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata); else - status=dlevmar_blec_dif(func, p, x, m, n, lb, ub, A, b, Arows, wghts, itmax, opts, info, NULL, covar, (void *)&mdata); + status=dlevmar_bleic_dif(func, p, x, m, n, lb, ub, A, b, Arows, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata); #else - mexErrMsgTxt("levmar: no box & linear constraints support, HAVE_LAPACK was not defined during MEX-file compilation."); + mexErrMsgTxt("levmar: no box, linear equation & inequality constraints support, HAVE_LAPACK was not defined during MEX-file compilation."); #endif /* HAVE_LAPACK */ #ifdef DEBUG fflush(stderr); - fprintf(stderr, "LEVMAR: calling dlevmar_blec_der()/dlevmar_blec_dif()\n"); + fprintf(stderr, "LEVMAR: calling dlevmar_bleic_der()/dlevmar_bleic_dif()\n"); #endif /* DEBUG */ - } + break; + case MIN_CONSTRAINED_BLIC: /* box, linear inequalities constraints */ +#ifdef HAVE_LAPACK + if(havejac) + status=dlevmar_bleic_der(func, jacfunc, p, x, m, n, lb, ub, NULL, NULL, 0, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata); + else + status=dlevmar_bleic_dif(func, p, x, m, n, lb, ub, NULL, NULL, 0, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata); +#else + mexErrMsgTxt("levmar: no box & linear inequality constraints support, HAVE_LAPACK was not defined during MEX-file compilation."); +#endif /* HAVE_LAPACK */ + +#ifdef DEBUG + fflush(stderr); + fprintf(stderr, "LEVMAR: calling dlevmar_blic_der()/dlevmar_blic_dif()\n"); +#endif /* DEBUG */ + break; + case MIN_CONSTRAINED_LEIC: /* linear equation & inequalities constraints */ +#ifdef HAVE_LAPACK + if(havejac) + status=dlevmar_bleic_der(func, jacfunc, p, x, m, n, NULL, NULL, A, b, Arows, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata); + else + status=dlevmar_bleic_dif(func, p, x, m, n, NULL, NULL, A, b, Arows, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata); +#else + mexErrMsgTxt("levmar: no linear equation & inequality constraints support, HAVE_LAPACK was not defined during MEX-file compilation."); +#endif /* HAVE_LAPACK */ + +#ifdef DEBUG + fflush(stderr); + fprintf(stderr, "LEVMAR: calling dlevmar_leic_der()/dlevmar_leic_dif()\n"); +#endif /* DEBUG */ + break; + case MIN_CONSTRAINED_LIC: /* linear inequalities constraints */ +#ifdef HAVE_LAPACK + if(havejac) + status=dlevmar_bleic_der(func, jacfunc, p, x, m, n, NULL, NULL, NULL, NULL, 0, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata); + else + status=dlevmar_bleic_dif(func, p, x, m, n, NULL, NULL, NULL, NULL, 0, C, d, Crows, itmax, opts, info, NULL, covar, (void *)&mdata); +#else + mexErrMsgTxt("levmar: no linear equation & inequality constraints support, HAVE_LAPACK was not defined during MEX-file compilation."); +#endif /* HAVE_LAPACK */ + +#ifdef DEBUG + fflush(stderr); + fprintf(stderr, "LEVMAR: calling dlevmar_lic_der()/dlevmar_lic_dif()\n"); +#endif /* DEBUG */ + break; + default: + mexErrMsgTxt("levmar: unexpected internal error."); } + #ifdef DEBUG fflush(stderr); printf("LEVMAR: minimization returned %d in %g iter, reason %g\n\tSolution: ", status, info[5], info[6]); @@ -568,6 +664,7 @@ cleanup: /* cleanup */ mxDestroyArray(mdata.rhs[0]); if(A) mxFree(A); + if(C) mxFree(C); mxFree(mdata.fname); if(havejac) mxFree(mdata.jacname); diff --git a/src/levmar/matlab/levmar.m b/src/levmar/matlab/levmar.m index 7adecaa..3dfaad0 100644 --- a/src/levmar/matlab/levmar.m +++ b/src/levmar/matlab/levmar.m @@ -2,11 +2,17 @@ function [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, typ % LEVMAR matlab MEX interface to the levmar non-linear least squares minimization % library available from http://www.ics.forth.gr/~lourakis/levmar/ % -% levmar can be used in any of the 4 following ways: +% levmar can be used in any of the 8 following ways: % [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'unc', ...) % [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'bc', lb, ub, ...) % [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'lec', A, b, ...) % [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'blec', lb, ub, A, b, wghts, ...) +% +% [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'bleic', lb, ub, A, b, C, d, ...) +% [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'blic', lb, ub, C, d, ...) +% [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'leic', A, b, C, d, ...) +% [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, 'lic', C, d, ...) +% % % The dots at the end denote any additional, problem specific data that are passed uniterpreted to % all invocations of fname and jacname, see below for details. @@ -44,14 +50,21 @@ function [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, typ % 'bc' specifies minimization subject to box constraints. % 'lec' specifies minimization subject to linear equation constraints. % 'blec' specifies minimization subject to box and linear equation constraints. +% 'bleic' specifies minimization subject to box, linear equation and inequality constraints. +% 'blic' specifies minimization subject to box and linear inequality constraints. +% 'leic' specifies minimization subject to linear equation and inequality constraints. +% 'lic' specifies minimization subject to linear inequality constraints. % If omitted, a default of 'unc' is assumed. Depending on the minimization type, the MEX -% interface will invoke one of dlevmar_XXX, dlevmar_bc_XXX, dlevmar_lec_XXX or dlevmar_blec_XXX +% interface will invoke one of dlevmar_XXX, dlevmar_bc_XXX, dlevmar_lec_XXX, dlevmar_blec_XXX or dlevmar_bleic_XXX % % - lb, ub: vectors of doubles specifying lower and upper bounds for p, respectively % % - A, b: k x m matrix and k vector specifying linear equation constraints for p, i.e. A*p=b % A should have full rank. % +% - C, d: k x m matrix and k vector specifying linear inequality constraints for p, i.e. C*p>=d +% A should have full rank. +% % - wghts: vector of doubles specifying the weights for the penalty terms corresponding to % the box constraints, see lmblec_core.c for more details. If omitted and a 'blec' type % minimization is to be carried out, default weights are used. diff --git a/src/levmar/matlab/lmdemo.m b/src/levmar/matlab/lmdemo.m index 4393011..d4cb96a 100644 --- a/src/levmar/matlab/lmdemo.m +++ b/src/levmar/matlab/lmdemo.m @@ -1,3 +1,8 @@ +% Demo program for levmar's MEX-file interface +% Performs minimization of several test problems + +format long; + % Unconstrained minimization % fitting the exponential model x_i=p(1)*exp(-p(2)*i)+p(3) of expfit.c to noisy measurements obtained with (5.0 0.1 1.0) @@ -38,6 +43,24 @@ disp('Meyer''s (reformulated) problem'); popt +% Osborne's problem +p0=[0.5, 1.5, -1.0, 1.0E-2, 2.0E-2]; + +x=[8.44E-1, 9.08E-1, 9.32E-1, 9.36E-1, 9.25E-1, 9.08E-1, 8.81E-1,... +8.5E-1, 8.18E-1, 7.84E-1, 7.51E-1, 7.18E-1, 6.85E-1, 6.58E-1,... +6.28E-1, 6.03E-1, 5.8E-1, 5.58E-1, 5.38E-1, 5.22E-1, 5.06E-1,... +4.9E-1, 4.78E-1, 4.67E-1, 4.57E-1, 4.48E-1, 4.38E-1, 4.31E-1,... +4.24E-1, 4.2E-1, 4.14E-1, 4.11E-1, 4.06E-1]; + + +options=[1E-03, 1E-15, 1E-15, 1E-20, 1E-06]; + +[ret, popt, info, covar]=levmar('osborne', 'jacosborne', p0, x, 200, options); +%[ret, popt, info, covar]=levmar('osborne', p0, x, 200, options, 'unc'); +disp('Osborne''s problem'); +popt + + % Linear constraints % Boggs-Tolle problem 3 @@ -72,7 +95,7 @@ popt % Box and linear constraints -% Hock-Schittkowski modified problem 52 +% Hock-Schittkowski modified problem 52 (#1) p0=[2.0, 2.0, 2.0, 2.0, 2.0]; x=[0.0, 0.0, 0.0, 0.0]; lb=[-0.09, 0.0, -realmax, -0.2, 0.0]; @@ -84,10 +107,10 @@ b=[0.0, 0.0, 0.0]'; options=[1E-03, 1E-15, 1E-15, 1E-20]; [ret, popt, info, covar]=levmar('modhs52', 'jacmodhs52', p0, x, 200, options, 'blec', lb, ub, A, b); -disp('Hock-Schittkowski modified problem 52'); +disp('Hock-Schittkowski modified problem 52 (#1)'); popt -% Hock-Schittkowski modified problem 235 +% Schittkowski modified problem 235 p0=[-2.0, 3.0, 1.0]; x=[0.0, 0.0]; lb=[-realmax, 0.1, 0.7]; @@ -101,3 +124,17 @@ options=[1E-03, 1E-15, 1E-15, 1E-20]; disp('Hock-Schittkowski modified problem 235'); popt +% Box, linear equation & inequality constraints +p0=[0.5, 0.5, 0.5, 0.5]; +x=[0.0, 0.0, 0.0, 0.0]; +lb=[0.0, 0.0, 0.0, 0.0]; +ub=[realmax, realmax, realmax, realmax]; +A=[0.0, 1.0, 4.0, 0.0]; +b=[1.5]'; +C=[-1.0, -2.0, -1.0, -1.0; + -3.0, -1.0, -2.0, 1.0]; +d=[-5.0, -0.4]'; + +[ret, popt, info, covar]=levmar('modhs76', 'jacmodhs76', p0, x, 200, options, 'bleic', lb, ub, A, b, C, d); +disp('Hock-Schittkowski modified problem 76'); +popt diff --git a/src/levmar/matlab/modhs76.m b/src/levmar/matlab/modhs76.m new file mode 100644 index 0000000..1e49583 --- /dev/null +++ b/src/levmar/matlab/modhs76.m @@ -0,0 +1,7 @@ +function x = modhs76(p) + + x(1)=p(1); + x(2)=sqrt(0.5)*p(2); + x(3)=p(3); + x(4)=sqrt(0.5)*p(4); + diff --git a/src/levmar/matlab/osborne.m b/src/levmar/matlab/osborne.m new file mode 100644 index 0000000..2d15373 --- /dev/null +++ b/src/levmar/matlab/osborne.m @@ -0,0 +1,7 @@ +function x = osborne(p) + n=33; + + for i=1:n + t=10*(i-1); + x(i)=p(1) + p(2)*exp(-p(4)*t) + p(3)*exp(-p(5)*t); + end diff --git a/src/levmar/misc.c b/src/levmar/misc.c index 25a72d5..1c6e996 100644 --- a/src/levmar/misc.c +++ b/src/levmar/misc.c @@ -28,7 +28,7 @@ #include #include -#include "lm.h" +#include "levmar.h" #include "misc.h" #if !defined(LM_DBL_PREC) && !defined(LM_SNGL_PREC) diff --git a/src/levmar/misc_core.c b/src/levmar/misc_core.c index 45fbc36..92722bc 100644 --- a/src/levmar/misc_core.c +++ b/src/levmar/misc_core.c @@ -579,8 +579,8 @@ int rnk; LM_REAL fact; #ifdef HAVE_LAPACK - rnk=LEVMAR_PSEUDOINVERSE(JtJ, C, m); - if(!rnk) return 0; + rnk=LEVMAR_PSEUDOINVERSE(JtJ, C, m); + if(!rnk) return 0; #else #ifdef _MSC_VER #pragma message("LAPACK not available, LU will be used for matrix inversion when computing the covariance; this might be unstable at times") @@ -592,10 +592,9 @@ LM_REAL fact; // rnk=LEVMAR_LUINVERSE(JtJ, C, m); - rnk = SVDINV(JtJ, C, m); + rnk = SVDINV(JtJ, C, m); - if (!rnk) - return 0; + if (!rnk) return 0; // rnk=m; /* assume full rank */ #endif /* HAVE_LAPACK */ @@ -689,12 +688,6 @@ register int i; #ifdef HAVE_LAPACK -/* compute the Cholesky decomposition of C in W, s.t. C=W^t W and W is upper triangular */ -int LEVMAR_CHOLESKY(LM_REAL *C, LM_REAL *W, int m) -{ -register int i, j; -int info; - /* compute the Cholesky decomposition of C in W, s.t. C=W^t W and W is upper triangular */ int LEVMAR_CHOLESKY(LM_REAL *C, LM_REAL *W, int m) { @@ -711,13 +704,13 @@ int info; POTF2("U", (int *)&m, W, (int *)&m, (int *)&info); /* error treatment */ if(info!=0){ - if(info<0){ + if(info<0){ fprintf(stderr, "LAPACK error: illegal value for argument %d of dpotf2 in %s\n", -info, LCAT(LEVMAR_CHOLESKY, "()")); - } - else{ - fprintf(stderr, "LAPACK error: the leading minor of order %d is not positive definite,\n%s()\n", info, - RCAT("and the Cholesky factorization could not be completed in ", LEVMAR_CHOLESKY)); - } + } + else{ + fprintf(stderr, "LAPACK error: the leading minor of order %d is not positive definite,\n%s()\n", info, + RCAT("and the Cholesky factorization could not be completed in ", LEVMAR_CHOLESKY)); + } return LM_ERROR; } @@ -783,12 +776,12 @@ register LM_REAL sum0=0.0, sum1=0.0, sum2=0.0, sum3=0.0; switch(n - i){ case 7 : e[i]=x[i]-y[i]; sum0+=e[i]*e[i]; ++i; - case 6 : e[i]=x[i]-y[i]; sum0+=e[i]*e[i]; ++i; - case 5 : e[i]=x[i]-y[i]; sum0+=e[i]*e[i]; ++i; - case 4 : e[i]=x[i]-y[i]; sum0+=e[i]*e[i]; ++i; + case 6 : e[i]=x[i]-y[i]; sum1+=e[i]*e[i]; ++i; + case 5 : e[i]=x[i]-y[i]; sum2+=e[i]*e[i]; ++i; + case 4 : e[i]=x[i]-y[i]; sum3+=e[i]*e[i]; ++i; case 3 : e[i]=x[i]-y[i]; sum0+=e[i]*e[i]; ++i; - case 2 : e[i]=x[i]-y[i]; sum0+=e[i]*e[i]; ++i; - case 1 : e[i]=x[i]-y[i]; sum0+=e[i]*e[i]; ++i; + case 2 : e[i]=x[i]-y[i]; sum1+=e[i]*e[i]; ++i; + case 1 : e[i]=x[i]-y[i]; sum2+=e[i]*e[i]; //++i; } } } @@ -818,12 +811,12 @@ register LM_REAL sum0=0.0, sum1=0.0, sum2=0.0, sum3=0.0; switch(n - i){ case 7 : e[i]=-y[i]; sum0+=e[i]*e[i]; ++i; - case 6 : e[i]=-y[i]; sum0+=e[i]*e[i]; ++i; - case 5 : e[i]=-y[i]; sum0+=e[i]*e[i]; ++i; - case 4 : e[i]=-y[i]; sum0+=e[i]*e[i]; ++i; + case 6 : e[i]=-y[i]; sum1+=e[i]*e[i]; ++i; + case 5 : e[i]=-y[i]; sum2+=e[i]*e[i]; ++i; + case 4 : e[i]=-y[i]; sum3+=e[i]*e[i]; ++i; case 3 : e[i]=-y[i]; sum0+=e[i]*e[i]; ++i; - case 2 : e[i]=-y[i]; sum0+=e[i]*e[i]; ++i; - case 1 : e[i]=-y[i]; sum0+=e[i]*e[i]; ++i; + case 2 : e[i]=-y[i]; sum1+=e[i]*e[i]; ++i; + case 1 : e[i]=-y[i]; sum2+=e[i]*e[i]; //++i; } } } diff --git a/src/profit.c b/src/profit.c index 35a08bb..adc84dd 100644 --- a/src/profit.c +++ b/src/profit.c @@ -9,7 +9,7 @@ * * Contents: Fit an arbitrary profile combination to a detection. * -* Last modify: 19/08/2010 +* Last modify: 26/08/2010 * *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ @@ -31,7 +31,7 @@ #include "globals.h" #include "prefs.h" #include "fits/fitscat.h" -#include "levmar/lm.h" +#include "levmar/levmar.h" #include "fft.h" #include "fitswcs.h" #include "check.h" @@ -199,14 +199,13 @@ OUTPUT Pointer to an allocated fit structure (containing details about the fit). NOTES It is a modified version of the lm_minimize() of lmfit. AUTHOR E. Bertin (IAP) -VERSION 08/03/2010 +VERSION 26/08/2010 ***/ void profit_fit(profitstruct *profit, picstruct *field, picstruct *wfield, objstruct *obj, obj2struct *obj2) { profitstruct pprofit; - profitstruct hdprofit; patternstruct *pattern; psfstruct *psf; checkstruct *check; @@ -2450,7 +2449,7 @@ INPUT Pointer to the profit structure, OUTPUT -. NOTES -. AUTHOR E. Bertin (IAP) -VERSION 02/07/2010 +VERSION 26/08/2010 ***/ void profit_resetparam(profitstruct *profit, paramenum paramtype) { @@ -2631,6 +2630,8 @@ void profit_resetparam(profitstruct *profit, paramenum paramtype) if (parammin!=parammax && (param<=parammin || param>=parammax)) param = (parammin+parammax)/2.0; + else if (parammin==0.0 && parammax==0.0) + parammax = 1.0; profit_setparam(profit, paramtype, param, parammin, parammax); return; diff --git a/xsl/sextractor.xsl b/xsl/sextractor.xsl index a6b7a1f..279b31f 100644 --- a/xsl/sextractor.xsl +++ b/xsl/sextractor.xsl @@ -16,8 +16,8 @@ - -