Skip to content
GitLab
Commit a9446199 authored by Emmanuel Bertin's avatar Emmanuel Bertin
Browse files

Fixed model-fitting issue with empty images in dual-image mode (thanks to V.... 

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.
parent 78801391
Hide whitespace changes
Inline Side-by-side
src/levmar/lm_core.c
View file @ a9446199
...@@ -37,6 +37,7 @@ ...@@ -37,6 +37,7 @@
#define AX_EQ_B_QR LM_ADD_PREFIX(Ax_eq_b_QR) #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_QRLS LM_ADD_PREFIX(Ax_eq_b_QRLS)
#define AX_EQ_B_SVD LM_ADD_PREFIX(Ax_eq_b_SVD) #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 #else
#define AX_EQ_B_LU LM_ADD_PREFIX(Ax_eq_b_LU_noLapack) #define AX_EQ_B_LU LM_ADD_PREFIX(Ax_eq_b_LU_noLapack)
#endif /* HAVE_LAPACK */ #endif /* HAVE_LAPACK */
...@@ -293,11 +294,12 @@ if(!(k%100)){ ...@@ -293,11 +294,12 @@ if(!(k%100)){
/* solve augmented equations */ /* solve augmented equations */
#ifdef HAVE_LAPACK #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; * 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 * 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_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_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; //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; ...@@ -516,7 +518,7 @@ int (*linsolver)(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)=NULL;
if(!work){ if(!work){
worksz=LM_DIF_WORKSZ(m, n); //4*n+4*m + n*m + m*m; 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){ if(!work){
fprintf(stderr, LCAT(LEVMAR_DIF, "(): memory allocation request failed\n")); fprintf(stderr, LCAT(LEVMAR_DIF, "(): memory allocation request failed\n"));
return LM_ERROR; return LM_ERROR;
...@@ -683,11 +685,12 @@ if(!(k%100)){ ...@@ -683,11 +685,12 @@ if(!(k%100)){
/* solve augmented equations */ /* solve augmented equations */
#ifdef HAVE_LAPACK #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; * 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 * 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_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_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; //issolved=AX_EQ_B_QR(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_QR;
...@@ -837,3 +840,4 @@ if(!(k%100)){ ...@@ -837,3 +840,4 @@ if(!(k%100)){
#undef AX_EQ_B_QR #undef AX_EQ_B_QR
#undef AX_EQ_B_QRLS #undef AX_EQ_B_QRLS
#undef AX_EQ_B_SVD #undef AX_EQ_B_SVD
#undef AX_EQ_B_BK
src/levmar/lmbc.c
View file @ a9446199
...@@ -28,7 +28,7 @@ ...@@ -28,7 +28,7 @@
#include <math.h> #include <math.h>
#include <float.h> #include <float.h>
#include "lm.h" #include "levmar.h"
#include "compiler.h" #include "compiler.h"
#include "misc.h" #include "misc.h"
......
src/levmar/lmbc_core.c
View file @ a9446199
...@@ -44,6 +44,7 @@ ...@@ -44,6 +44,7 @@
#define AX_EQ_B_QR LM_ADD_PREFIX(Ax_eq_b_QR) #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_QRLS LM_ADD_PREFIX(Ax_eq_b_QRLS)
#define AX_EQ_B_SVD LM_ADD_PREFIX(Ax_eq_b_SVD) #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 #else
#define AX_EQ_B_LU LM_ADD_PREFIX(Ax_eq_b_LU_noLapack) #define AX_EQ_B_LU LM_ADD_PREFIX(Ax_eq_b_LU_noLapack)
#endif /* HAVE_LAPACK */ #endif /* HAVE_LAPACK */
...@@ -551,11 +552,12 @@ if(!(k%100)){ ...@@ -551,11 +552,12 @@ if(!(k%100)){
/* solve augmented equations */ /* solve augmented equations */
#ifdef HAVE_LAPACK #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; * 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 * 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_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_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; //issolved=AX_EQ_B_QR(jacTjac, jacTe, Dp, m); ++nlss; linsolver=AX_EQ_B_QR;
...@@ -943,3 +945,4 @@ int ret; ...@@ -943,3 +945,4 @@ int ret;
#undef AX_EQ_B_QR #undef AX_EQ_B_QR
#undef AX_EQ_B_QRLS #undef AX_EQ_B_QRLS
#undef AX_EQ_B_SVD #undef AX_EQ_B_SVD
#undef AX_EQ_B_BK
src/levmar/lmblec.c
View file @ a9446199
...@@ -28,17 +28,15 @@ ...@@ -28,17 +28,15 @@
#include <math.h> #include <math.h>
#include <float.h> #include <float.h>
#include "lm.h" #include "levmar.h"
#include "misc.h" #include "misc.h"
#ifndef HAVE_LAPACK #ifndef HAVE_LAPACK
#ifdef _MSC_VER #ifdef _MSC_VER
#pragma message("Combined box and linearly constrained optimization requires LAPACK and was not compiled!") #pragma message("Combined box and linearly constrained optimization requires LAPACK and was not compiled!")
/*
#else #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 #endif // _MSC_VER
#else // LAPACK present #else // LAPACK present
......
src/levmar/lmbleic.c 0 → 100644
View file @ a9446199
/////////////////////////////////////////////////////////////////////////////////
//
// 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 <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#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 */
src/levmar/lmbleic_core.c 0 → 100644
View file @ a9446199
/////////////////////////////////////////////////////////////////////////////////
//
// 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<n; ++i){
jacextimm=jacext+i*mm;
jacim=jac+i*m;
for(j=0; j<m; ++j)
jacextimm[j]=jacim[j]; //jacext[i*mm+j]=jac[i*m+j];
for(j=m; j<mm; ++j)
jacextimm[j]=0.0; //jacext[i*mm+j]=0.0;
}
}
/*
* This function is similar to LEVMAR_DER except that the minimization is
* performed subject to the box constraints lb[i]<=p[i]<=ub[i], the linear
* equation constraints A*p=b, A being k1xm, b k1x1, and the linear inequality
* constraints C*p>=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<m-k1){
fprintf(stderr, LCAT(LEVMAR_BLEIC_DER, "(): cannot solve a problem with fewer measurements + equality constraints [%d + %d] than unknowns [%d]\n"), n, k1, m);
return LM_ERROR;
}
k12=k1+k2;
ptr=(LM_REAL *)malloc((3*mm + k12*mm + k12 + n*m + (covar? mm*mm : 0))*sizeof(LM_REAL));
if(!ptr){
fprintf(stderr, LCAT(LEVMAR_BLEIC_DER, "(): memory allocation request failed\n"));
return LM_ERROR;
}
pext=ptr;
lbext=pext+mm;
ubext=lbext+mm;
Aext=ubext+mm;
bext=Aext+k12*mm;
data.jac=bext+k12;
covext=covar? data.jac+n*m : NULL;
data.nineqcnstr=k2;
data.func=func;
data.jacf=jacf;
data.adata=adata;
/* compute y s.t. C*p - y=d, i.e. y=C*p-d.
* y is stored in the last k2 elements of pext
*/
for(i=0; i<k2; ++i){
for(j=0, tmp=0.0; j<m; ++j)
tmp+=C[i*m+j]*p[j];
pext[j=i+m]=tmp-d[i];
/* surplus variables must be >=0 */
lbext[j]=0.0;
ubext[j]=LM_REAL_MAX;
}
/* set the first m elements of pext equal to p */
for(i=0; i<m; ++i){
pext[i]=p[i];
lbext[i]=lb? lb[i] : LM_REAL_MIN;
ubext[i]=ub? ub[i] : LM_REAL_MAX;
}
/* setup the constraints matrix */
/* original linear equation constraints */
for(i=0; i<k1; ++i){
for(j=0; j<m; ++j)
Aext[i*mm+j]=A[i*m+j];
for(j=m; j<mm; ++j)
Aext[i*mm+j]=0.0;
bext[i]=b[i];
}
/* linear equation constraints resulting from surplus variables */
for(i=0, ii=k1; i<k2; ++i, ++ii){
for(j=0; j<m; ++j)
Aext[ii*mm+j]=C[i*m+j];
for(j=m; j<mm; ++j)
Aext[ii*mm+j]=0.0;
Aext[ii*mm+m+i]=-1.0;
bext[ii]=d[i];
}
if(!info) info=locinfo; /* make sure that LEVMAR_BLEC_DER() is called with non-null info */
/* note that the default weights for the penalty terms are being used below */
ret=LEVMAR_BLEC_DER(LMBLEIC_FUNC, LMBLEIC_JACF, pext, x, mm, n, lbext, ubext, Aext, bext, k12, NULL, itmax, opts, info, work, covext, (void *)&data);
/* copy back the minimizer */
for(i=0; i<m; ++i)
p[i]=pext[i];
#if 0
printf("Surplus variables for the minimizer:\n");
for(i=m; i<mm; ++i)
printf("%g ", pext[i]);
printf("\n\n");
#endif
if(covar){
for(i=0; i<m; ++i){
for(j=0; j<m; ++j)
covar[i*m+j]=covext[i*mm+j];
}
}
free(ptr);
return ret;
}
/* Similar to the LEVMAR_BLEIC_DER() function above, except that the Jacobian is approximated
* with the aid of finite differences (forward or central, see the comment for the opts argument)
*/
int LEVMAR_BLEIC_DIF(
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 */
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[5], /* I: opts[0-3] = minim. options [\mu, \epsilon1, \epsilon2, \epsilon3, \delta]. Respectively the
* scale factor for initial \mu, stopping thresholds for ||J^T e||_inf, ||Dp||_2 and ||e||_2 and
* the step used in difference approximation to the Jacobian. Set to NULL for defaults to be used.
* If \delta<0, the Jacobian is approximated with central differences which are more accurate
* (but slower!) compared to the forward differences employed by default.
*/
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_DIF_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.
* 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(!C || !d){
fprintf(stderr, RCAT(LCAT(LEVMAR_BLEIC_DIF, "(): missing inequality constraints, use "), LEVMAR_BLEC_DIF) "() in this case!\n");
return LM_ERROR;
}
if(!A || !b) k1=0; // sanity check
mm=m+k2;
if(n<m-k1){
fprintf(stderr, LCAT(LEVMAR_BLEIC_DIF, "(): cannot solve a problem with fewer measurements + equality constraints [%d + %d] than unknowns [%d]\n"), n, k1, m);
return LM_ERROR;
}
k12=k1+k2;
ptr=(LM_REAL *)malloc((3*mm + k12*mm + k12 + (covar? mm*mm : 0))*sizeof(LM_REAL));
if(!ptr){
fprintf(stderr, LCAT(LEVMAR_BLEIC_DIF, "(): memory allocation request failed\n"));
return LM_ERROR;
}
pext=ptr;
lbext=pext+mm;
ubext=lbext+mm;
Aext=ubext+mm;
bext=Aext+k12*mm;
data.jac=NULL;
covext=covar? bext+k12 : NULL;
data.nineqcnstr=k2;
data.func=func;
data.jacf=NULL;
data.adata=adata;
/* compute y s.t. C*p - y=d, i.e. y=C*p-d.
* y is stored in the last k2 elements of pext
*/
for(i=0; i<k2; ++i){
for(j=0, tmp=0.0; j<m; ++j)
tmp+=C[i*m+j]*p[j];
pext[j=i+m]=tmp-d[i];
/* surplus variables must be >=0 */
lbext[j]=0.0;
ubext[j]=LM_REAL_MAX;
}
/* set the first m elements of pext equal to p */
for(i=0; i<m; ++i){
pext[i]=p[i];
lbext[i]=lb? lb[i] : LM_REAL_MIN;
ubext[i]=ub? ub[i] : LM_REAL_MAX;
}
/* setup the constraints matrix */
/* original linear equation constraints */
for(i=0; i<k1; ++i){
for(j=0; j<m; ++j)
Aext[i*mm+j]=A[i*m+j];
for(j=m; j<mm; ++j)
Aext[i*mm+j]=0.0;
bext[i]=b[i];
}
/* linear equation constraints resulting from surplus variables */
for(i=0, ii=k1; i<k2; ++i, ++ii){
for(j=0; j<m; ++j)
Aext[ii*mm+j]=C[i*m+j];
for(j=m; j<mm; ++j)
Aext[ii*mm+j]=0.0;
Aext[ii*mm+m+i]=-1.0;
bext[ii]=d[i];
}
if(!info) info=locinfo; /* make sure that LEVMAR_BLEC_DIF() is called with non-null info */
/* note that the default weights for the penalty terms are being used below */
ret=LEVMAR_BLEC_DIF(LMBLEIC_FUNC, pext, x, mm, n, lbext, ubext, Aext, bext, k12, NULL, itmax, opts, info, work, covext, (void *)&data);
/* copy back the minimizer */
for(i=0; i<m; ++i)
p[i]=pext[i];
#if 0
printf("Surplus variables for the minimizer:\n");
for(i=m; i<mm; ++i)
printf("%g ", pext[i]);
printf("\n\n");
#endif
if(covar){
for(i=0; i<m; ++i){
for(j=0; j<m; ++j)
covar[i*m+j]=covext[i*mm+j];
}
}
free(ptr);
return ret;
}
/* convenience wrappers to LEVMAR_BLEIC_DER/LEVMAR_BLEIC_DIF */
/* box & linear inequality constraints */
int LEVMAR_BLIC_DER(
void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata),
void (*jacf)(LM_REAL *p, LM_REAL *j, int m, int n, void *adata),
LM_REAL *p, LM_REAL *x, int m, int n,
LM_REAL *lb, LM_REAL *ub,
LM_REAL *C, LM_REAL *d, int k2,
int itmax, LM_REAL opts[4], LM_REAL info[LM_INFO_SZ], LM_REAL *work, LM_REAL *covar, void *adata)
{
return LEVMAR_BLEIC_DER(func, jacf, p, x, m, n, lb, ub, NULL, NULL, 0, C, d, k2, itmax, opts, info, work, covar, adata);
}
int LEVMAR_BLIC_DIF(
void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata),
LM_REAL *p, LM_REAL *x, int m, int n,
LM_REAL *lb, LM_REAL *ub,
LM_REAL *C, LM_REAL *d, int k2,
int itmax, LM_REAL opts[5], LM_REAL info[LM_INFO_SZ], LM_REAL *work, LM_REAL *covar, void *adata)
{
return LEVMAR_BLEIC_DIF(func, p, x, m, n, lb, ub, NULL, NULL, 0, C, d, k2, itmax, opts, info, work, covar, adata);
}
/* linear equation & inequality constraints */
int LEVMAR_LEIC_DER(
void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata),
void (*jacf)(LM_REAL *p, LM_REAL *j, int m, int n, void *adata),
LM_REAL *p, LM_REAL *x, int m, int n,
LM_REAL *A, LM_REAL *b, int k1,
LM_REAL *C, LM_REAL *d, int k2,
int itmax, LM_REAL opts[4], LM_REAL info[LM_INFO_SZ], LM_REAL *work, LM_REAL *covar, void *adata)
{
return LEVMAR_BLEIC_DER(func, jacf, p, x, m, n, NULL, NULL, A, b, k1, C, d, k2, itmax, opts, info, work, covar, adata);
}
int LEVMAR_LEIC_DIF(
void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata),
LM_REAL *p, LM_REAL *x, int m, int n,
LM_REAL *A, LM_REAL *b, int k1,
LM_REAL *C, LM_REAL *d, int k2,
int itmax, LM_REAL opts[5], LM_REAL info[LM_INFO_SZ], LM_REAL *work, LM_REAL *covar, void *adata)
{
return LEVMAR_BLEIC_DIF(func, p, x, m, n, NULL, NULL, A, b, k1, C, d, k2, itmax, opts, info, work, covar, adata);
}
/* linear inequality constraints */
int LEVMAR_LIC_DER(
void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata),
void (*jacf)(LM_REAL *p, LM_REAL *j, int m, int n, void *adata),
LM_REAL *p, LM_REAL *x, int m, int n,
LM_REAL *C, LM_REAL *d, int k2,
int itmax, LM_REAL opts[4], LM_REAL info[LM_INFO_SZ], LM_REAL *work, LM_REAL *covar, void *adata)
{
return LEVMAR_BLEIC_DER(func, jacf, p, x, m, n, NULL, NULL, NULL, NULL, 0, C, d, k2, itmax, opts, info, work, covar, adata);
}
int LEVMAR_LIC_DIF(
void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata),
LM_REAL *p, LM_REAL *x, int m, int n,
LM_REAL *C, LM_REAL *d, int k2,
int itmax, LM_REAL opts[5], LM_REAL info[LM_INFO_SZ], LM_REAL *work, LM_REAL *covar, void *adata)
{
return LEVMAR_BLEIC_DIF(func, p, x, m, n, NULL, NULL, NULL, NULL, 0, C, d, k2, itmax, opts, info, work, covar, adata);
}
/* undefine all. THIS MUST REMAIN AT THE END OF THE FILE */
#undef LMBLEIC_DATA
#undef LMBLEIC_ELIM
#undef LMBLEIC_FUNC
#undef LMBLEIC_JACF
#undef LEVMAR_FDIF_FORW_JAC_APPROX
#undef LEVMAR_COVAR
#undef LEVMAR_TRANS_MAT_MAT_MULT
#undef LEVMAR_BLEIC_DER
#undef LEVMAR_BLEIC_DIF
#undef LEVMAR_BLIC_DER
#undef LEVMAR_BLIC_DIF
#undef LEVMAR_LEIC_DER
#undef LEVMAR_LEIC_DIF
#undef LEVMAR_LIC_DER
#undef LEVMAR_LIC_DIF
#undef LEVMAR_BLEC_DER
#undef LEVMAR_BLEC_DIF
src/levmar/lmdemo.c
View file @ a9446199
...@@ -28,7 +28,7 @@ ...@@ -28,7 +28,7 @@
#include <math.h> #include <math.h>
#include <float.h> #include <float.h>
#include "lm.h" #include "levmar.h"
#ifndef LM_DBL_PREC #ifndef LM_DBL_PREC
#error Demo program assumes that levmar has been compiled with double precision, see 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; ...@@ -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<n; ++i){
t=10*i;
x[i]=p[0] + p[1]*exp(-p[3]*t) + p[2]*exp(-p[4]*t);
}
}
void jacosborne(double *p, double *jac, int m, int n, void *data)
{
register int i, j;
double t, tmp1, tmp2;
for(i=j=0; i<n; ++i){
t=10*i;
tmp1=exp(-p[3]*t);
tmp2=exp(-p[4]*t);
jac[j++]=1.0;
jac[j++]=tmp1;
jac[j++]=tmp2;
jac[j++]=-p[1]*t*tmp1;
jac[j++]=-p[2]*t*tmp2;
}
}
/* helical valley function, minimum at (1.0, 0.0, 0.0) */ /* helical valley function, minimum at (1.0, 0.0, 0.0) */
#ifndef M_PI #ifndef M_PI
#define M_PI 3.14159265358979323846 /* pi */ #define M_PI 3.14159265358979323846 /* pi */
...@@ -471,7 +501,7 @@ register int j=0; ...@@ -471,7 +501,7 @@ register int j=0;
jac[j++]=1.0; jac[j++]=1.0;
} }
/* Hock - Schittkowski (modified) problem 52 (box/linearly constrained), minimum at (-0.09, 0.03, 0.25, -0.19, 0.03) /* Hock - Schittkowski (modified #1) problem 52 (box/linearly constrained), minimum at (-0.09, 0.03, 0.25, -0.19, 0.03)
* constr1: p[0] + 3*p[1] = 0; * constr1: p[0] + 3*p[1] = 0;
* constr2: p[2] + p[3] - 2*p[4] = 0; * constr2: p[2] + p[3] - 2*p[4] = 0;
* constr3: p[1] - p[4] = 0; * constr3: p[1] - p[4] = 0;
...@@ -484,7 +514,7 @@ register int j=0; ...@@ -484,7 +514,7 @@ register int j=0;
* constr8: 0.0 <= p[4] <= 0.3; * constr8: 0.0 <= p[4] <= 0.3;
* *
*/ */
void modhs52(double *p, double *x, int m, int n, void *data) void mod1hs52(double *p, double *x, int m, int n, void *data)
{ {
x[0]=4.0*p[0]-p[1]; x[0]=4.0*p[0]-p[1];
x[1]=p[1]+p[2]-2.0; x[1]=p[1]+p[2]-2.0;
...@@ -492,7 +522,7 @@ void modhs52(double *p, double *x, int m, int n, void *data) ...@@ -492,7 +522,7 @@ void modhs52(double *p, double *x, int m, int n, void *data)
x[3]=p[4]-1.0; x[3]=p[4]-1.0;
} }
void jacmodhs52(double *p, double *jac, int m, int n, void *data) void jacmod1hs52(double *p, double *jac, int m, int n, void *data)
{ {
register int j=0; register int j=0;
...@@ -521,6 +551,60 @@ register int j=0; ...@@ -521,6 +551,60 @@ register int j=0;
jac[j++]=1.0; jac[j++]=1.0;
} }
/* Hock - Schittkowski (modified #2) problem 52 (linear inequality constrained), minimum at (0.5, 2.0, 0.0, 1.0, 1.0)
* A fifth term [(p[0]-0.5)^2] is added to the objective function and
* the equality contraints are replaced by the following inequalities:
* constr1: p[0] + 3*p[1] >= -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) /* Schittkowski (modified) problem 235 (box/linearly constrained), minimum at (-1.725, 2.9, 0.725)
* constr1: p[0] + p[2] = -1.0; * constr1: p[0] + p[2] = -1.0;
* *
...@@ -653,13 +737,55 @@ register int j=0; ...@@ -653,13 +737,55 @@ register int j=0;
jac[j+3]=R9*p[1]+2*p[3]; 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() int main()
{ {
register int i, j; register int i, j;
int problem, ret; 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) x[16]; // 16 is max(2, 3, 5, 6, 16)
int m, n; int m, n;
double opts[LM_OPTS_SZ], info[LM_INFO_SZ]; double opts[LM_OPTS_SZ], info[LM_INFO_SZ];
...@@ -669,6 +795,7 @@ char *probname[]={ ...@@ -669,6 +795,7 @@ char *probname[]={
"Powell's function", "Powell's function",
"Wood's function", "Wood's function",
"Meyer's (reformulated) problem", "Meyer's (reformulated) problem",
"Osborne's problem",
"helical valley function", "helical valley function",
"Boggs & Tolle's problem #3", "Boggs & Tolle's problem #3",
"Hock - Schittkowski problem #28", "Hock - Schittkowski problem #28",
...@@ -679,13 +806,15 @@ char *probname[]={ ...@@ -679,13 +806,15 @@ char *probname[]={
"hatfldb problem", "hatfldb problem",
"hatfldc problem", "hatfldc problem",
"equilibrium combustion problem", "equilibrium combustion problem",
"Hock - Schittkowski modified problem #52", "Hock - Schittkowski modified #1 problem #52",
"Schittkowski modified problem #235", "Schittkowski modified problem #235",
"Boggs & Tolle modified problem #7", "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[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! //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 */ /* uncomment the appropriate line below to select a minimization problem */
...@@ -695,12 +824,13 @@ char *probname[]={ ...@@ -695,12 +824,13 @@ char *probname[]={
//2; // Powell's function //2; // Powell's function
//3; // Wood's function //3; // Wood's function
4; // Meyer's (reformulated) problem 4; // Meyer's (reformulated) problem
//5; // helical valley function //5; // Osborne's problem
//6; // helical valley function
#ifdef HAVE_LAPACK #ifdef HAVE_LAPACK
//6; // Boggs & Tolle's problem 3 //7; // Boggs & Tolle's problem 3
//7; // Hock - Schittkowski problem 28 //8; // Hock - Schittkowski problem 28
//8; // Hock - Schittkowski problem 48 //9; // Hock - Schittkowski problem 48
//9; // Hock - Schittkowski problem 51 //10; // Hock - Schittkowski problem 51
#else // no LAPACK #else // no LAPACK
#ifdef _MSC_VER #ifdef _MSC_VER
#pragma message("LAPACK not available, some test problems cannot be used") #pragma message("LAPACK not available, some test problems cannot be used")
...@@ -709,21 +839,24 @@ char *probname[]={ ...@@ -709,21 +839,24 @@ char *probname[]={
#endif // _MSC_VER #endif // _MSC_VER
#endif /* HAVE_LAPACK */ #endif /* HAVE_LAPACK */
//10; // Hock - Schittkowski problem 01 //11; // Hock - Schittkowski problem 01
//11; // Hock - Schittkowski modified problem 21 //12; // Hock - Schittkowski modified problem 21
//12; // hatfldb problem //13; // hatfldb problem
//13; // hatfldc problem //14; // hatfldc problem
//14; // equilibrium combustion problem //15; // equilibrium combustion problem
#ifdef HAVE_LAPACK #ifdef HAVE_LAPACK
//15; // Hock - Schittkowski modified problem 52 //16; // Hock - Schittkowski modified #1 problem 52
//16; // Schittkowski modified problem 235 //17; // Schittkowski modified problem 235
//17; // Boggs & Tolle modified problem #7 //18; // Boggs & Tolle modified problem #7
//19; // Hock - Schittkowski modified #2 problem 52
//20; // Hock - Schittkowski modified problem #76"
#endif /* HAVE_LAPACK */ #endif /* HAVE_LAPACK */
switch(problem){ switch(problem){
default: fprintf(stderr, "unknown problem specified (#%d)! Note that some minimization problems require LAPACK.\n", problem); default: fprintf(stderr, "unknown problem specified (#%d)! Note that some minimization problems require LAPACK.\n", problem);
exit(1); exit(1);
break; break;
case 0: case 0:
/* Rosenbrock function */ /* Rosenbrock function */
m=2; n=2; m=2; n=2;
...@@ -798,10 +931,27 @@ char *probname[]={ ...@@ -798,10 +931,27 @@ char *probname[]={
for(i=0; i<n; ++i) printf("gradient %d, err %g\n", i, err[i]); for(i=0; i<n; ++i) printf("gradient %d, err %g\n", i, err[i]);
} }
*/ */
break; break;
case 5: case 5:
/* Osborne's data fitting problem */
{
double x33[]={
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};
m=5; n=33;
p[0]=0.5; p[1]=1.5; p[2]=-1.0; p[3]=1.0E-2; p[4]=2.0E-2;
ret=dlevmar_der(osborne, jacosborne, p, x33, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
//ret=dlevmar_dif(osborne, p, x33, m, n, 1000, opts, info, NULL, NULL, NULL); // no Jacobian
}
break;
case 6:
/* helical valley function */ /* helical valley function */
m=3; n=3; m=3; n=3;
p[0]=-1.0; p[1]=0.0; p[2]=0.0; p[0]=-1.0; p[1]=0.0; p[2]=0.0;
...@@ -811,7 +961,7 @@ char *probname[]={ ...@@ -811,7 +961,7 @@ char *probname[]={
break; break;
#ifdef HAVE_LAPACK #ifdef HAVE_LAPACK
case 6: case 7:
/* Boggs-Tolle problem 3 */ /* Boggs-Tolle problem 3 */
m=5; n=5; m=5; n=5;
p[0]=2.0; p[1]=2.0; p[2]=2.0; p[0]=2.0; p[1]=2.0; p[2]=2.0;
...@@ -826,7 +976,8 @@ char *probname[]={ ...@@ -826,7 +976,8 @@ char *probname[]={
//ret=dlevmar_lec_dif(bt3, p, x, m, n, A, b, 3, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, no Jacobian //ret=dlevmar_lec_dif(bt3, p, x, m, n, A, b, 3, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, no Jacobian
} }
break; break;
case 7:
case 8:
/* Hock - Schittkowski problem 28 */ /* Hock - Schittkowski problem 28 */
m=3; n=3; m=3; n=3;
p[0]=-4.0; p[1]=1.0; p[2]=1.0; p[0]=-4.0; p[1]=1.0; p[2]=1.0;
...@@ -840,7 +991,8 @@ char *probname[]={ ...@@ -840,7 +991,8 @@ char *probname[]={
//ret=dlevmar_lec_dif(hs28, p, x, m, n, A, b, 1, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, no Jacobian //ret=dlevmar_lec_dif(hs28, p, x, m, n, A, b, 1, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, no Jacobian
} }
break; break;
case 8:
case 9:
/* Hock - Schittkowski problem 48 */ /* Hock - Schittkowski problem 48 */
m=5; n=5; m=5; n=5;
p[0]=3.0; p[1]=5.0; p[2]=-3.0; p[0]=3.0; p[1]=5.0; p[2]=-3.0;
...@@ -855,7 +1007,8 @@ char *probname[]={ ...@@ -855,7 +1007,8 @@ char *probname[]={
//ret=dlevmar_lec_dif(hs48, p, x, m, n, A, b, 2, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, no Jacobian //ret=dlevmar_lec_dif(hs48, p, x, m, n, A, b, 2, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, no Jacobian
} }
break; break;
case 9:
case 10:
/* Hock - Schittkowski problem 51 */ /* Hock - Schittkowski problem 51 */
m=5; n=5; m=5; n=5;
p[0]=2.5; p[1]=0.5; p[2]=2.0; p[0]=2.5; p[1]=0.5; p[2]=2.0;
...@@ -870,9 +1023,10 @@ char *probname[]={ ...@@ -870,9 +1023,10 @@ char *probname[]={
//ret=dlevmar_lec_dif(hs51, p, x, m, n, A, b, 3, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, no Jacobian //ret=dlevmar_lec_dif(hs51, p, x, m, n, A, b, 3, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, no Jacobian
} }
break; break;
#endif /* HAVE_LAPACK */ #endif /* HAVE_LAPACK */
case 10: case 11:
/* Hock - Schittkowski problem 01 */ /* Hock - Schittkowski problem 01 */
m=2; n=2; m=2; n=2;
p[0]=-2.0; p[1]=1.0; p[0]=-2.0; p[1]=1.0;
...@@ -887,7 +1041,8 @@ char *probname[]={ ...@@ -887,7 +1041,8 @@ char *probname[]={
ret=dlevmar_bc_der(hs01, jachs01, p, x, m, n, lb, ub, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian ret=dlevmar_bc_der(hs01, jachs01, p, x, m, n, lb, ub, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
} }
break; break;
case 11:
case 12:
/* Hock - Schittkowski (modified) problem 21 */ /* Hock - Schittkowski (modified) problem 21 */
m=2; n=2; m=2; n=2;
p[0]=-1.0; p[1]=-1.0; p[0]=-1.0; p[1]=-1.0;
...@@ -902,7 +1057,8 @@ char *probname[]={ ...@@ -902,7 +1057,8 @@ char *probname[]={
ret=dlevmar_bc_der(hs21, jachs21, p, x, m, n, lb, ub, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian ret=dlevmar_bc_der(hs21, jachs21, p, x, m, n, lb, ub, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
} }
break; break;
case 12:
case 13:
/* hatfldb problem */ /* hatfldb problem */
m=4; n=4; m=4; n=4;
p[0]=p[1]=p[2]=p[3]=0.1; p[0]=p[1]=p[2]=p[3]=0.1;
...@@ -919,7 +1075,8 @@ char *probname[]={ ...@@ -919,7 +1075,8 @@ char *probname[]={
ret=dlevmar_bc_der(hatfldb, jachatfldb, p, x, m, n, lb, ub, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian ret=dlevmar_bc_der(hatfldb, jachatfldb, p, x, m, n, lb, ub, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
} }
break; break;
case 13:
case 14:
/* hatfldc problem */ /* hatfldc problem */
m=4; n=4; m=4; n=4;
p[0]=p[1]=p[2]=p[3]=0.9; p[0]=p[1]=p[2]=p[3]=0.9;
...@@ -935,7 +1092,8 @@ char *probname[]={ ...@@ -935,7 +1092,8 @@ char *probname[]={
ret=dlevmar_bc_der(hatfldc, jachatfldc, p, x, m, n, lb, ub, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian ret=dlevmar_bc_der(hatfldc, jachatfldc, p, x, m, n, lb, ub, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
} }
break; break;
case 14:
case 15:
/* equilibrium combustion problem */ /* equilibrium combustion problem */
m=5; n=5; m=5; n=5;
p[0]=p[1]=p[2]=p[3]=p[4]=0.0001; p[0]=p[1]=p[2]=p[3]=p[4]=0.0001;
...@@ -951,9 +1109,10 @@ char *probname[]={ ...@@ -951,9 +1109,10 @@ char *probname[]={
ret=dlevmar_bc_der(combust, jaccombust, p, x, m, n, lb, ub, 5000, opts, info, NULL, NULL, NULL); // with analytic Jacobian ret=dlevmar_bc_der(combust, jaccombust, p, x, m, n, lb, ub, 5000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
} }
break; break;
#ifdef HAVE_LAPACK #ifdef HAVE_LAPACK
case 15: case 16:
/* Hock - Schittkowski modified problem 52 */ /* Hock - Schittkowski modified #1 problem 52 */
m=5; n=4; m=5; n=4;
p[0]=2.0; p[1]=2.0; p[2]=2.0; p[0]=2.0; p[1]=2.0; p[2]=2.0;
p[3]=2.0; p[4]=2.0; p[3]=2.0; p[4]=2.0;
...@@ -970,11 +1129,12 @@ char *probname[]={ ...@@ -970,11 +1129,12 @@ char *probname[]={
lb[0]=-0.09; lb[1]=0.0; lb[2]=-DBL_MAX; lb[3]=-0.2; lb[4]=0.0; lb[0]=-0.09; lb[1]=0.0; lb[2]=-DBL_MAX; lb[3]=-0.2; lb[4]=0.0;
ub[0]=DBL_MAX; ub[1]=0.3; ub[2]=0.25; ub[3]=0.3; ub[4]=0.3; ub[0]=DBL_MAX; ub[1]=0.3; ub[2]=0.25; ub[3]=0.3; ub[4]=0.3;
ret=dlevmar_blec_der(modhs52, jacmodhs52, p, x, m, n, lb, ub, A, b, 3, weights, 1000, opts, info, NULL, NULL, NULL); // box & lin. constraints, analytic Jacobian ret=dlevmar_blec_der(mod1hs52, jacmod1hs52, p, x, m, n, lb, ub, A, b, 3, weights, 1000, opts, info, NULL, NULL, NULL); // box & lin. constraints, analytic Jacobian
//ret=dlevmar_blec_dif(modhs52, p, x, m, n, lb, ub, A, b, 3, weights, 1000, opts, info, NULL, NULL, NULL); // box & lin. constraints, no Jacobian //ret=dlevmar_blec_dif(mod1hs52, p, x, m, n, lb, ub, A, b, 3, weights, 1000, opts, info, NULL, NULL, NULL); // box & lin. constraints, no Jacobian
} }
break; break;
case 16:
case 17:
/* Schittkowski modified problem 235 */ /* Schittkowski modified problem 235 */
m=3; n=2; m=3; n=2;
p[0]=-2.0; p[1]=3.0; p[2]=1.0; p[0]=-2.0; p[1]=3.0; p[2]=1.0;
...@@ -993,7 +1153,8 @@ char *probname[]={ ...@@ -993,7 +1153,8 @@ char *probname[]={
//ret=dlevmar_blec_dif(mods235, p, x, m, n, lb, ub, A, b, 2, NULL, 1000, opts, info, NULL, NULL, NULL); // box & lin. constraints, no Jacobian //ret=dlevmar_blec_dif(mods235, p, x, m, n, lb, ub, A, b, 2, NULL, 1000, opts, info, NULL, NULL, NULL); // box & lin. constraints, no Jacobian
} }
break; break;
case 17:
case 18:
/* Boggs & Tolle modified problem 7 */ /* Boggs & Tolle modified problem 7 */
m=5; n=5; m=5; n=5;
p[0]=-2.0; p[1]=1.0; p[2]=1.0; p[3]=1.0; p[4]=1.0; p[0]=-2.0; p[1]=1.0; p[2]=1.0; p[3]=1.0; p[4]=1.0;
...@@ -1012,6 +1173,47 @@ char *probname[]={ ...@@ -1012,6 +1173,47 @@ char *probname[]={
//ret=dlevmar_blec_dif(modbt7, p, x, m, n, lb, ub, A, b, 3, NULL, 10000, opts, info, NULL, NULL, NULL); // box & lin. constraints, no Jacobian //ret=dlevmar_blec_dif(modbt7, p, x, m, n, lb, ub, A, b, 3, NULL, 10000, opts, info, NULL, NULL, NULL); // box & lin. constraints, no Jacobian
} }
break; break;
case 19:
/* Hock - Schittkowski modified #2 problem 52 */
m=5; n=5;
p[0]=2.0; p[1]=2.0; p[2]=2.0;
p[3]=2.0; p[4]=2.0;
for(i=0; i<n; i++) x[i]=0.0;
{
double C[3*5]={1.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, -2.0, 0.0, -1.0, 0.0, 0.0, 1.0},
d[3]={-1.0, -2.0, -7.0};
ret=dlevmar_bleic_der(mod2hs52, jacmod2hs52, p, x, m, n, NULL, NULL, NULL, NULL, 0, C, d, 3, 1000, opts, info, NULL, NULL, NULL); // lin. ineq. constraints, analytic Jacobian
//ret=dlevmar_bleic_dif(mod2hs52, p, x, m, n, NULL, NULL, NULL, NULL, 0, C, d, 3, 1000, opts, info, NULL, NULL, NULL); // lin. ineq. constraints, no Jacobian
}
break;
case 20:
/* Hock - Schittkowski modified problem 76 */
m=4; n=4;
p[0]=0.5; p[1]=0.5; p[2]=0.5; p[3]=0.5;
for(i=0; i<n; i++) x[i]=0.0;
{
double A[1*4]={0.0, 1.0, 4.0, 0.0},
b[1]={1.5};
double C[2*4]={-1.0, -2.0, -1.0, -1.0, -3.0, -1.0, -2.0, 1.0},
d[2]={-5.0, -0.4};
double lb[4]={0.0, 0.0, 0.0, 0.0};
ret=dlevmar_bleic_der(modhs76, jacmodhs76, p, x, m, n, lb, NULL, A, b, 1, C, d, 2, 1000, opts, info, NULL, NULL, NULL); // lin. ineq. constraints, analytic Jacobian
//ret=dlevmar_bleic_dif(modhs76, p, x, m, n, lb, NULL, A, b, 1, C, d, 2, 1000, opts, info, NULL, NULL, NULL); // lin. ineq. constraints, no Jacobian
/* variations:
* if no lb is used, the minimizer is (-0.1135922 0.1330097 0.3417476 0.07572816)
* if the rhs of constr2 is 4.0, the minimizer is (0.0, 0.166667, 0.333333, 0.0)
*/
}
break;
#endif /* HAVE_LAPACK */ #endif /* HAVE_LAPACK */
} /* switch */ } /* switch */
......
src/levmar/lmlec.c
View file @ a9446199
...@@ -27,7 +27,7 @@ ...@@ -27,7 +27,7 @@
#include <stdlib.h> #include <stdlib.h>
#include <math.h> #include <math.h>
#include "lm.h" #include "levmar.h"
#include "misc.h" #include "misc.h"
...@@ -35,10 +35,8 @@ ...@@ -35,10 +35,8 @@
#ifdef _MSC_VER #ifdef _MSC_VER
#pragma message("Linearly constrained optimization requires LAPACK and was not compiled!") #pragma message("Linearly constrained optimization requires LAPACK and was not compiled!")
/*
#else #else
#warning Linearly constrained optimization requires LAPACK and was not compiled! //#warning Linearly constrained optimization requires LAPACK and was not compiled!
*/
#endif // _MSC_VER #endif // _MSC_VER
#else // LAPACK present #else // LAPACK present
......
src/levmar/matlab/CMakeLists.txt 0 → 100644
View file @ a9446199
# 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)
src/levmar/matlab/README.txt
View file @ a9446199
This directory contains a matlab MEX interface to levmar. This interface 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. 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 FILES
The following files are included: The following files are included:
......
src/levmar/matlab/jacbt3.m
View file @ a9446199
function jac = bt3_jac(p, adata) function jac = jacbt3(p, adata)
n=5; n=5;
m=5; m=5;
......
src/levmar/matlab/jacexpfit.m
View file @ a9446199
function jac = expfit_jac(p, data) function jac = jacexpfit(p, data)
n=data; n=data;
m=max(size(p)); m=max(size(p));
......
src/levmar/matlab/jachs01.m
View file @ a9446199
function jac = hs01_jac(p) function jac = jachs01(p)
m=2; m=2;
jac(1, 1:m)=[-20.0*p(1), 10.0]; jac(1, 1:m)=[-20.0*p(1), 10.0];
......
src/levmar/matlab/jacmeyer.m
View file @ a9446199
function jac = meyer_jac(p, data1, data2) function jac = jacmeyer(p, data1, data2)
n=16; n=16;
m=3; m=3;
......
src/levmar/matlab/jacmodhs52.m
View file @ a9446199
function jac = modhs52_jac(p) function jac = jacmodhs52(p)
m=5; m=5;
jac(1, 1:m)=[4.0, -1.0, 0.0, 0.0, 0.0]; jac(1, 1:m)=[4.0, -1.0, 0.0, 0.0, 0.0];
......
src/levmar/matlab/jacmodhs76.m 0 → 100644
View file @ a9446199
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)];
src/levmar/matlab/jacosborne.m 0 → 100644
View file @ a9446199
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
src/levmar/matlab/levmar.c
View file @ a9446199
...@@ -24,7 +24,7 @@ ...@@ -24,7 +24,7 @@
#include <string.h> #include <string.h>
#include <ctype.h> #include <ctype.h>
#include <lm.h> #include <levmar.h>
#include <mex.h> #include <mex.h>
...@@ -46,6 +46,10 @@ ...@@ -46,6 +46,10 @@
#define MIN_CONSTRAINED_BC 1 #define MIN_CONSTRAINED_BC 1
#define MIN_CONSTRAINED_LEC 2 #define MIN_CONSTRAINED_LEC 2
#define MIN_CONSTRAINED_BLEC 3 #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 { struct mexdata {
/* matlab names of the fitting function & its Jacobian */ /* matlab names of the fitting function & its Jacobian */
...@@ -68,9 +72,9 @@ static void matlabFmtdErrMsgTxt(char *fmt, ...) ...@@ -68,9 +72,9 @@ static void matlabFmtdErrMsgTxt(char *fmt, ...)
char buf[256]; char buf[256];
va_list args; va_list args;
va_start(args, fmt); va_start(args, fmt);
vsprintf(buf, fmt, args); vsprintf(buf, fmt, args);
va_end(args); va_end(args);
mexErrMsgTxt(buf); mexErrMsgTxt(buf);
} }
...@@ -81,9 +85,9 @@ static void matlabFmtdWarnMsgTxt(char *fmt, ...) ...@@ -81,9 +85,9 @@ static void matlabFmtdWarnMsgTxt(char *fmt, ...)
char buf[256]; char buf[256];
va_list args; va_list args;
va_start(args, fmt); va_start(args, fmt);
vsprintf(buf, fmt, args); vsprintf(buf, fmt, args);
va_end(args); va_end(args);
mexWarnMsgTxt(buf); mexWarnMsgTxt(buf);
} }
...@@ -183,7 +187,7 @@ double *mp; ...@@ -183,7 +187,7 @@ double *mp;
else if(!mxIsDouble(lhs[0]) || mxIsComplex(lhs[0]) || !(mxGetM(lhs[0])==1 || mxGetN(lhs[0])==1) || else if(!mxIsDouble(lhs[0]) || mxIsComplex(lhs[0]) || !(mxGetM(lhs[0])==1 || mxGetN(lhs[0])==1) ||
__MAX__(mxGetM(lhs[0]), mxGetN(lhs[0]))!=n){ __MAX__(mxGetM(lhs[0]), mxGetN(lhs[0]))!=n){
fprintf(stderr, "levmar: '%s' should produce a real vector with %d elements (got %d).\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; ret=1;
} }
/* delete the matrix created by matlab */ /* delete the matrix created by matlab */
...@@ -220,20 +224,26 @@ double *mp; ...@@ -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_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_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_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[]) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *Prhs[])
{ {
register int i; register int i;
register double *pdbl; register double *pdbl;
mxArray **prhs=(mxArray **)&Prhs[0], *At; mxArray **prhs=(mxArray **)&Prhs[0], *At, *Ct;
struct mexdata mdata; struct mexdata mdata;
int len, status; int len, status;
double *p, *p0, *ret, *x; 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 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 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 */ /* parse input args; start by checking their number */
if((nrhs<5)) if((nrhs<5))
...@@ -362,6 +372,10 @@ if(!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) || !(mxGetM(prhs[1])==1 && mxGet ...@@ -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, "bc", 2)) mintype=MIN_CONSTRAINED_BC;
else if(!strncmp(minhowto, "lec", 3)) mintype=MIN_CONSTRAINED_LEC; else if(!strncmp(minhowto, "lec", 3)) mintype=MIN_CONSTRAINED_LEC;
else if(!strncmp(minhowto, "blec", 4)) mintype=MIN_CONSTRAINED_BLEC; 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); else matlabFmtdErrMsgTxt("levmar: unknown minimization type '%s'.", minhowto);
mxFree(minhowto); mxFree(minhowto);
...@@ -378,7 +392,7 @@ if(!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) || !(mxGetM(prhs[1])==1 && mxGet ...@@ -378,7 +392,7 @@ if(!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) || !(mxGetM(prhs[1])==1 && mxGet
* upon the minimization type determined above * upon the minimization type determined above
*/ */
/** lb, ub **/ /** 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 */ /* 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[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(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 ...@@ -397,7 +411,7 @@ if(!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) || !(mxGetM(prhs[1])==1 && mxGet
} }
/** A, b **/ /** 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 */ /* 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[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(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 ...@@ -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 /* arguments below this point are assumed to be extra arguments passed
* to every invocation of the fitting function and its Jacobian * to every invocation of the fitting function and its Jacobian
*/ */
...@@ -471,8 +506,8 @@ extraargs: ...@@ -471,8 +506,8 @@ extraargs:
covar=mxMalloc(m*m*sizeof(double)); covar=mxMalloc(m*m*sizeof(double));
/* invoke levmar */ /* invoke levmar */
if(!lb && !ub){ switch(mintype){
if(!A && !b){ /* no constraints */ case MIN_UNCONSTRAINED: /* no constraints */
if(havejac) if(havejac)
status=dlevmar_der(func, jacfunc, p, x, m, n, itmax, opts, info, NULL, covar, (void *)&mdata); status=dlevmar_der(func, jacfunc, p, x, m, n, itmax, opts, info, NULL, covar, (void *)&mdata);
else else
...@@ -481,8 +516,18 @@ extraargs: ...@@ -481,8 +516,18 @@ extraargs:
fflush(stderr); fflush(stderr);
fprintf(stderr, "LEVMAR: calling dlevmar_der()/dlevmar_dif()\n"); fprintf(stderr, "LEVMAR: calling dlevmar_der()/dlevmar_dif()\n");
#endif /* DEBUG */ #endif /* DEBUG */
} break;
else{ /* linear constraints */ 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 #ifdef HAVE_LAPACK
if(havejac) if(havejac)
status=dlevmar_lec_der(func, jacfunc, p, x, m, n, A, b, Arows, itmax, opts, info, NULL, covar, (void *)&mdata); 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: ...@@ -496,35 +541,86 @@ extraargs:
fflush(stderr); fflush(stderr);
fprintf(stderr, "LEVMAR: calling dlevmar_lec_der()/dlevmar_lec_dif()\n"); fprintf(stderr, "LEVMAR: calling dlevmar_lec_der()/dlevmar_lec_dif()\n");
#endif /* DEBUG */ #endif /* DEBUG */
} break;
} case MIN_CONSTRAINED_BLEC: /* box & linear equation constraints */
else{ #ifdef HAVE_LAPACK
if(!A && !b){ /* box constraints */
if(havejac) 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 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 #ifdef DEBUG
fflush(stderr); 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 */ #endif /* DEBUG */
} break;
else{ /* box & linear constraints */ case MIN_CONSTRAINED_BLEIC: /* box, linear equation & inequalities constraints */
#ifdef HAVE_LAPACK #ifdef HAVE_LAPACK
if(havejac) 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 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 #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 */ #endif /* HAVE_LAPACK */
#ifdef DEBUG #ifdef DEBUG
fflush(stderr); 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 */ #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 #ifdef DEBUG
fflush(stderr); fflush(stderr);
printf("LEVMAR: minimization returned %d in %g iter, reason %g\n\tSolution: ", status, info[5], info[6]); printf("LEVMAR: minimization returned %d in %g iter, reason %g\n\tSolution: ", status, info[5], info[6]);
...@@ -568,6 +664,7 @@ cleanup: ...@@ -568,6 +664,7 @@ cleanup:
/* cleanup */ /* cleanup */
mxDestroyArray(mdata.rhs[0]); mxDestroyArray(mdata.rhs[0]);
if(A) mxFree(A); if(A) mxFree(A);
if(C) mxFree(C);
mxFree(mdata.fname); mxFree(mdata.fname);
if(havejac) mxFree(mdata.jacname); if(havejac) mxFree(mdata.jacname);
......
src/levmar/matlab/levmar.m
View file @ a9446199
...@@ -2,11 +2,17 @@ function [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, typ ...@@ -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 % LEVMAR matlab MEX interface to the levmar non-linear least squares minimization
% library available from http://www.ics.forth.gr/~lourakis/levmar/ % 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, '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, '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, '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, '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 % 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. % 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 ...@@ -44,14 +50,21 @@ function [ret, popt, info, covar]=levmar(fname, jacname, p0, x, itmax, opts, typ
% 'bc' specifies minimization subject to box constraints. % 'bc' specifies minimization subject to box constraints.
% 'lec' specifies minimization subject to linear equation constraints. % 'lec' specifies minimization subject to linear equation constraints.
% 'blec' specifies minimization subject to box and 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 % 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 % - 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, b: k x m matrix and k vector specifying linear equation constraints for p, i.e. A*p=b
% A should have full rank. % 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 % - 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 % 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. % minimization is to be carried out, default weights are used.
......
src/levmar/matlab/lmdemo.m
View file @ a9446199
% Demo program for levmar's MEX-file interface
% Performs minimization of several test problems
format long;
% Unconstrained minimization % 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) % 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'); ...@@ -38,6 +43,24 @@ disp('Meyer''s (reformulated) problem');
popt 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 % Linear constraints
% Boggs-Tolle problem 3 % Boggs-Tolle problem 3
...@@ -72,7 +95,7 @@ popt ...@@ -72,7 +95,7 @@ popt
% Box and linear constraints % 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]; p0=[2.0, 2.0, 2.0, 2.0, 2.0];
x=[0.0, 0.0, 0.0, 0.0]; x=[0.0, 0.0, 0.0, 0.0];
lb=[-0.09, 0.0, -realmax, -0.2, 0.0]; lb=[-0.09, 0.0, -realmax, -0.2, 0.0];
...@@ -84,10 +107,10 @@ b=[0.0, 0.0, 0.0]'; ...@@ -84,10 +107,10 @@ b=[0.0, 0.0, 0.0]';
options=[1E-03, 1E-15, 1E-15, 1E-20]; 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); [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 popt
% Hock-Schittkowski modified problem 235 % Schittkowski modified problem 235
p0=[-2.0, 3.0, 1.0]; p0=[-2.0, 3.0, 1.0];
x=[0.0, 0.0]; x=[0.0, 0.0];
lb=[-realmax, 0.1, 0.7]; lb=[-realmax, 0.1, 0.7];
...@@ -101,3 +124,17 @@ options=[1E-03, 1E-15, 1E-15, 1E-20]; ...@@ -101,3 +124,17 @@ options=[1E-03, 1E-15, 1E-15, 1E-20];
disp('Hock-Schittkowski modified problem 235'); disp('Hock-Schittkowski modified problem 235');
popt 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
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment