From facb9273c571b6bef430e71b31ff229abd6d1a46 Mon Sep 17 00:00:00 2001
From: Mark Abraham <mark.j.abraham@gmail.com>
Date: Mon, 11 Jun 2018 08:30:19 +0200
Subject: [PATCH] Bump lmfit support to version 7
The breaking API change means that distributions will not be able to
build GROMACS reliably with support for lmfit of arbitrary versions.
Nor does lmfit provide any support for querying the version. Tools
that call such fitting will now issue a fatal error.
Updated the FindLmfit and other cmake code to make better use of
modern cmake idioms for managing build targets. lmfit support is now
handled by a tri-state GMX_USE_LMFIT cmake variable, which defaults to
using the bundled version.
Reorganized aspects of gmx_lmcurve to better encapsulate the
dependency, now that we have to support the absence of lmfit.
Updated install guide and COPYING file.
Refs #2533
Change-Id: I6b14e08be216f803326b0ad9215b8d89bb0962c1
---
diff --git a/COPYING b/COPYING
index 4cc52a2..ba2710e 100644
--- a/COPYING
+++ b/COPYING
@@ -1249,7 +1249,7 @@
--
Copyright (c) 1980-1999 University of Chicago,
as operator of Argonne National Laboratory
- Copyright (c) 2004-2013 Joachim Wuttke, Forschungszentrum Juelich GmbH
+ Copyright (c) 2004-2015 Joachim Wuttke, Forschungszentrum Juelich GmbH
All rights reserved.
diff --git a/cmake/FindLmfit.cmake b/cmake/FindLmfit.cmake
index 5bf489a..9dd23f1 100644
--- a/cmake/FindLmfit.cmake
+++ b/cmake/FindLmfit.cmake
@@ -1,7 +1,7 @@
#
# This file is part of the GROMACS molecular simulation package.
#
-# Copyright (c) 2016, by the GROMACS development team, led by
+# Copyright (c) 2016,2018, by the GROMACS development team, led by
# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
# and including many others, as listed in the AUTHORS file in the
# top-level source directory and at http://www.gromacs.org.
@@ -32,19 +32,20 @@
# To help us fund GROMACS development, we humbly ask that you cite
# the research papers on the package. Check out http://www.gromacs.org.
-# This package tries to find an external lmfit library. It is intended
-# to work with pkg-config, because that is the mechanism supported in
-# lmfit. Upon exit, the following variables may be set:
+# This package tries to find an external lmfit library, version
+# 7. Note that pkg-config support was removed before version 7, so is
+# no longer supported here. Upon exit, the following variables may be
+# set:
#
# LMFIT_FOUND - lmfit was found
# LMFIT_INCLUDE_DIR - lmfit include directory
-# LMFIT_LIBRARIES - lmfit libraries
+# LMFIT_LIBRARY - lmfit library
# LMFIT_LINKS_OK - lmfit libraries link correctly
# LMFIT_VERSION - lmfit version string as "major.minor"
#
-# If you cannot use pkg-config for some reason, then setting
-# LMFIT_INCLUDE_DIRS and LMFIT_LIBRARY_DIRS on the cmake command line
-# to suitable values will work.
+# CMake will search the CMAKE_PREFIX_PATH in the usual way, but if you
+# need more control then setting LMFIT_INCLUDE_DIR and LMFIT_LIBRARY
+# on the cmake command line to suitable values will work.
include(CMakePushCheckState)
cmake_push_check_state()
@@ -55,12 +56,12 @@
# lmfit doesn't support CMake-based find_package version
# checking in 6.1, so this code does nothing.
if(LMFIT_FIND_VERSION_EXACT)
- pkg_check_modules(PC_LMFIT lmfit=${LMFIT_FIND_VERSION})
+ pkg_check_modules(PC_LMFIT QUIET lmfit=${LMFIT_FIND_VERSION})
else()
- pkg_check_modules(PC_LMFIT lmfit>=${LMFIT_FIND_VERSION})
+ pkg_check_modules(PC_LMFIT QUIET lmfit>=${LMFIT_FIND_VERSION})
endif()
else()
- pkg_check_modules(PC_LMFIT lmfit)
+ pkg_check_modules(PC_LMFIT QUIET lmfit)
if (PC_LMFIT_VERSION)
string(REGEX REPLACE "^([0-9]+):([0-9]+)" "\\1.\\2" LMFIT_VERSION "${PC_LMFIT_VERSION}")
endif()
@@ -68,20 +69,15 @@
endif()
# Try to find lmfit, perhaps with help from pkg-config
-find_path(LMFIT_INCLUDE_DIRS lmcurve.h HINTS "${PC_LMFIT_INCLUDE_DIRS}" PATH_SUFFIXES include)
-find_library(LMFIT_LIBRARY_DIRS NAMES lmfit HINTS "${PC_LMFIT_LIBRARY_DIRS}" PATH_SUFFIXES lib64 lib)
+find_path(LMFIT_INCLUDE_DIR lmcurve.h HINTS "${PC_LMFIT_INCLUDE_DIRS}" PATH_SUFFIXES include)
+find_library(LMFIT_LIBRARY NAMES lmfit HINTS "${PC_LMFIT_LIBRARY_DIRS}" PATH_SUFFIXES lib64 lib)
# Make sure we can also link, so that cross-compilation is properly supported
-if (LMFIT_INCLUDE_DIRS AND LMFIT_LIBRARY_DIRS)
+if (LMFIT_INCLUDE_DIR AND LMFIT_LIBRARY)
include(CheckCXXSourceCompiles)
- set(CMAKE_REQUIRED_INCLUDES ${LMFIT_INCLUDE_DIRS})
- set(CMAKE_REQUIRED_LIBRARIES ${LMFIT_LIBRARY_DIRS})
+ set(CMAKE_REQUIRED_INCLUDES ${LMFIT_INCLUDE_DIR})
+ set(CMAKE_REQUIRED_LIBRARIES ${LMFIT_LIBRARY})
check_cxx_source_compiles("#include <lmcurve.h>\nint main(){lmcurve(0,0,0,0,0,0,0,0);}" LMFIT_LINKS_OK)
-endif()
-
-if (LMFIT_LINKS_OK)
- set(LMFIT_INCLUDE_DIR ${LMFIT_INCLUDE_DIRS})
- set(LMFIT_LIBRARIES ${LMFIT_LIBRARY_DIRS})
endif()
include(FindPackageHandleStandardArgs)
@@ -90,11 +86,19 @@
LMFIT_FOUND
REQUIRED_VARS
LMFIT_INCLUDE_DIR
- LMFIT_LIBRARIES
+ LMFIT_LIBRARY
LMFIT_LINKS_OK
VERSION_VAR
LMFIT_VERSION)
-mark_as_advanced(LMFIT_INCLUDE_DIRS LMFIT_LIBRARY_DIRS)
+mark_as_advanced(LMFIT_INCLUDE_DIR LMFIT_LIBRARY)
+
+if (LMFIT_FOUND)
+ add_library(lmfit INTERFACE IMPORTED)
+ set_target_properties(lmfit PROPERTIES
+ INTERFACE_INCLUDE_DIRECTORIES "${LMFIT_INCLUDE_DIR}"
+ INTERFACE_LINK_LIBRARIES "${LMFIT_LIBRARY}"
+ )
+endif()
cmake_pop_check_state()
diff --git a/cmake/gmxManageLmfit.cmake b/cmake/gmxManageLmfit.cmake
index 182928f..5953a19 100644
--- a/cmake/gmxManageLmfit.cmake
+++ b/cmake/gmxManageLmfit.cmake
@@ -1,7 +1,7 @@
#
# This file is part of the GROMACS molecular simulation package.
#
-# Copyright (c) 2016, by the GROMACS development team, led by
+# Copyright (c) 2016,2018, by the GROMACS development team, led by
# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
# and including many others, as listed in the AUTHORS file in the
# top-level source directory and at http://www.gromacs.org.
@@ -32,21 +32,16 @@
# To help us fund GROMACS development, we humbly ask that you cite
# the research papers on the package. Check out http://www.gromacs.org.
-set(GMX_LMFIT_MINIMUM_REQUIRED_VERSION "6.1")
-set(GMX_BUNDLED_LMFIT_DIR "${CMAKE_SOURCE_DIR}/src/external/lmfit")
+# Note that lmfit does not have a stable API, so GROMACS only supports
+# the same version that it bundles.
+set(GMX_LMFIT_REQUIRED_VERSION "7.0")
-option(GMX_EXTERNAL_LMFIT "Use external lmfit instead of compiling the version bundled with GROMACS." OFF)
-mark_as_advanced(GMX_EXTERNAL_LMFIT)
+include(gmxOptionUtilities)
-macro(manage_lmfit)
- if(GMX_EXTERNAL_LMFIT)
- # Find an external lmfit library.
- find_package(Lmfit ${GMX_LMFIT_MINIMUM_REQUIRED_VERSION})
- if(NOT LMFIT_FOUND)
- message(FATAL_ERROR "External lmfit could not be found, please adjust your pkg-config path to include the lmfit.pc file")
- endif()
- endif()
-endmacro()
+gmx_option_multichoice(GMX_USE_LMFIT
+ "Use external lmfit instead of compiling the version bundled with GROMACS."
+ INTERNAL
+ INTERNAL EXTERNAL NONE)
macro(get_lmfit_properties LMFIT_SOURCES_VAR LMFIT_LIBRARIES_VAR LMFIT_INCLUDE_DIR_VAR LMFIT_INCLUDE_DIR_ORDER_VAR)
if (GMX_EXTERNAL_LMFIT)
@@ -62,5 +57,29 @@
endif()
endmacro()
-manage_lmfit()
+function(manage_lmfit)
+ if(GMX_USE_LMFIT STREQUAL "INTERNAL")
+ set(BUNDLED_LMFIT_DIR "${CMAKE_SOURCE_DIR}/src/external/lmfit")
+ file(GLOB LMFIT_SOURCES ${BUNDLED_LMFIT_DIR}/*.cpp)
+ # Create a fake library for lmfit for libgromacs to depend on
+ add_library(lmfit INTERFACE)
+ target_sources(lmfit INTERFACE ${LMFIT_SOURCES})
+ target_include_directories(lmfit INTERFACE
+ $<BUILD_INTERFACE:${BUNDLED_LMFIT_DIR}>)
+ target_link_libraries(libgromacs PRIVATE lmfit)
+
+ set(HAVE_LMFIT_VALUE TRUE)
+ elseif(GMX_USE_LMFIT STREQUAL "EXTERNAL")
+ # Find an external lmfit library.
+ find_package(Lmfit ${GMX_LMFIT_MINIMUM_REQUIRED_VERSION})
+ if(NOT LMFIT_FOUND)
+ message(FATAL_ERROR "External lmfit could not be found, please adjust your pkg-config path to include the lmfit.pc file")
+ endif()
+ target_link_libraries(libgromacs PRIVATE lmfit)
+ set(HAVE_LMFIT_VALUE TRUE)
+ else()
+ set(HAVE_LMFIT_VALUE FALSE)
+ endif()
+ set(HAVE_LMFIT ${HAVE_LMFIT_VALUE} CACHE BOOL "Whether lmfit library support is enabled")
+endfunction()
diff --git a/docs/CMakeLists.txt b/docs/CMakeLists.txt
index 9b38d75..3f050b5 100644
--- a/docs/CMakeLists.txt
+++ b/docs/CMakeLists.txt
@@ -366,7 +366,7 @@
REQUIRED_CUDA_COMPUTE_CAPABILITY REGRESSIONTEST_VERSION
SOURCE_MD5SUM REGRESSIONTEST_MD5SUM_STRING
GMX_TNG_MINIMUM_REQUIRED_VERSION
- GMX_LMFIT_MINIMUM_REQUIRED_VERSION
+ GMX_LMFIT_REQUIRED_VERSION
COMMENT "Configuring Sphinx configuration file")
gmx_add_sphinx_input_file(${SPHINX_CONFIG_VARS_FILE})
gmx_add_sphinx_source_files(FILES ${SPHINX_SOURCE_FILES})
diff --git a/docs/conf-vars.py.cmakein b/docs/conf-vars.py.cmakein
index f4c3e02..038014a 100644
--- a/docs/conf-vars.py.cmakein
+++ b/docs/conf-vars.py.cmakein
@@ -1,7 +1,7 @@
#
# This file is part of the GROMACS molecular simulation package.
#
-# Copyright (c) 2015,2016,2017, by the GROMACS development team, led by
+# Copyright (c) 2015,2016,2017,2018, by the GROMACS development team, led by
# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
# and including many others, as listed in the AUTHORS file in the
# top-level source directory and at http://www.gromacs.org.
@@ -48,5 +48,5 @@
('SOURCE_MD5SUM', '@SOURCE_MD5SUM@'),
('REGRESSIONTEST_MD5SUM', '@REGRESSIONTEST_MD5SUM_STRING@'),
('GMX_TNG_MINIMUM_REQUIRED_VERSION', '@GMX_TNG_MINIMUM_REQUIRED_VERSION@'),
- ('GMX_LMFIT_MINIMUM_REQUIRED_VERSION', '@GMX_LMFIT_MINIMUM_REQUIRED_VERSION@')
+ ('GMX_LMFIT_REQUIRED_VERSION', '@GMX_LMFIT_REQUIRED_VERSION@')
]
diff --git a/docs/install-guide/index.rst b/docs/install-guide/index.rst
index 48fe2c3..5324629 100644
--- a/docs/install-guide/index.rst
+++ b/docs/install-guide/index.rst
@@ -340,10 +340,14 @@
by setting ``-DGMX_EXTERNAL_TNG=yes``, but TNG
|GMX_TNG_MINIMUM_REQUIRED_VERSION| is bundled in the |Gromacs|
source already.
-* An external lmfit library for Levenberg-Marquardt curve fitting
- can be used by setting ``-DGMX_EXTERNAL_LMFIT=yes``, but lmfit
- |GMX_LMFIT_MINIMUM_REQUIRED_VERSION| is bundled in the |Gromacs|
- source already.
+* The lmfit library for Levenberg-Marquardt curve fitting is used in
+ |Gromacs|. Only lmfit |GMX_LMFIT_REQUIRED_VERSION| is supported. A
+ reduced version of that library is bundled in the |Gromacs|
+ distribution, and the default build uses it. That default may be
+ explicitly enabled with ``-DGMX_USE_LMFIT=internal``. To use an
+ external lmfit library, set ``-DGMX_USE_LMFIT=external``, and adjust
+ ``CMAKE_PREFIX_PATH`` as needed. lmfit support can be disabled with
+ ``-DGMX_USE_LMFIT=none``.
* zlib is used by TNG for compressing some kinds of trajectory data
* Building the |Gromacs| documentation is optional, and requires
ImageMagick, pdflatex, bibtex, doxygen, python 2.7, sphinx
diff --git a/src/config.h.cmakein b/src/config.h.cmakein
index b7e5964..d0c3c35 100644
--- a/src/config.h.cmakein
+++ b/src/config.h.cmakein
@@ -377,6 +377,9 @@
/* Define if we have feenableexcept */
#cmakedefine01 HAVE_FEENABLEEXCEPT
+/* Define if we have lmfit support */
+#cmakedefine01 HAVE_LMFIT
+
/*! \endcond */
#endif
diff --git a/src/external/lmfit/README b/src/external/lmfit/README
index 3c0557b..4b8d7c6 100644
--- a/src/external/lmfit/README
+++ b/src/external/lmfit/README
@@ -1,5 +1,5 @@
This directory contains only the barebones version of the
-lmfit package, version 6.1. Full version is available from
+lmfit package, version 7.0. Full version is available from
http://apps.jcns.fz-juelich.de/doku/sc/lmfit
The license under which lmfit is bundled may be found in the GROMACS
diff --git a/src/external/lmfit/lmmin.cpp b/src/external/lmfit/lmmin.cpp
index eb2d298..813fe5d 100644
--- a/src/external/lmfit/lmmin.cpp
+++ b/src/external/lmfit/lmmin.cpp
@@ -16,40 +16,41 @@
#include <assert.h>
#include <stdlib.h>
#include <stdio.h>
-#include <cmath>
+#include <math.h>
#include <float.h>
#include "lmmin.h"
-using namespace std;
-
-static double lm_enorm(int n, const double* x);
-
-#define MIN(a, b) (((a) <= (b)) ? (a) : (b))
-#define MAX(a, b) (((a) >= (b)) ? (a) : (b))
-#define SQR(x) (x) * (x)
+#define MIN(a,b) (((a)<=(b)) ? (a) : (b))
+#define MAX(a,b) (((a)>=(b)) ? (a) : (b))
+#define SQR(x) (x)*(x)
/* Declare functions that do the heavy numerics.
Implementions are in this source file, below lmmin.
- Dependences: lmmin calls lmpar, which calls qrfac and qrsolv. */
-static void lm_lmpar(const int n, double* r, const int ldr, const int* Pivot,
- const double* diag, const double* qtb, const double delta,
- double* par, double* x, double* Sdiag, double* aux, double* xdi);
-static void lm_qrfac(const int m, const int n, double* A, int* Pivot, double* Rdiag,
- double* Acnorm, double* W);
-static void lm_qrsolv(const int n, double* r, const int ldr, const int* Pivot,
- const double* diag, const double* qtb, double* x,
- double* Sdiag, double* W);
+ Dependences: lmmin calls qrfac and lmpar; lmpar calls qrsolv. */
+void lm_lmpar(
+ const int n, double *const r, const int ldr, int *const ipvt,
+ double *const diag, double *const qtb, double delta, double *const par,
+ double *const x,
+ double *const sdiag, double *const aux, double *const xdi);
+void lm_qrfac(
+ const int m, const int n, double *const a, int *const ipvt,
+ double *const rdiag, double *const acnorm, double *const wa );
+void lm_qrsolv(
+ const int n, double *const r, const int ldr, int *const ipvt,
+ double *const diag, double *const qtb, double *const x,
+ double *const sdiag, double *const wa);
-/******************************************************************************/
-/* Numeric constants */
-/******************************************************************************/
-/* Set machine-dependent constants to values from float.h. */
-#define LM_MACHEP DBL_EPSILON /* resolution of arithmetic */
-#define LM_DWARF DBL_MIN /* smallest nonzero number */
+/*****************************************************************************/
+/* Numeric constants */
+/*****************************************************************************/
+
+/* machine-dependent constants from float.h */
+#define LM_MACHEP DBL_EPSILON /* resolution of arithmetic */
+#define LM_DWARF DBL_MIN /* smallest nonzero number */
#define LM_SQRT_DWARF sqrt(DBL_MIN) /* square should not underflow */
#define LM_SQRT_GIANT sqrt(DBL_MAX) /* square should not overflow */
-#define LM_USERTOL 30 * LM_MACHEP /* users are recommended to require this */
+#define LM_USERTOL 30*LM_MACHEP /* users are recommended to require this */
/* If the above values do not work, the following seem good for an x86:
LM_MACHEP .555e-16
@@ -65,19 +66,19 @@
LM_USER_TOL 1.e-14
*/
-/* Predefined control parameter sets (msgfile=NULL means stdout). */
const lm_control_struct lm_control_double = {
- LM_USERTOL, LM_USERTOL, LM_USERTOL, LM_USERTOL,
- 100., 100, 1, NULL, 0, -1, -1};
+ LM_USERTOL, LM_USERTOL, LM_USERTOL, LM_USERTOL, 100., 100, 1,
+ NULL, 0, -1, -1 };
const lm_control_struct lm_control_float = {
- 1.e-7, 1.e-7, 1.e-7, 1.e-7,
- 100., 100, 1, NULL, 0, -1, -1};
+ 1.e-7, 1.e-7, 1.e-7, 1.e-7, 100., 100, 1,
+ NULL, 0, -1, -1 };
-/******************************************************************************/
-/* Message texts (indexed by status.info) */
-/******************************************************************************/
-const char* lm_infmsg[] = {
+/*****************************************************************************/
+/* Message texts (indexed by status.info) */
+/*****************************************************************************/
+
+const char *lm_infmsg[] = {
"found zero (sum of squares below underflow limit)",
"converged (the relative error in the sum of squares is at most tol)",
"converged (the relative error of the parameter vector is at most tol)",
@@ -90,66 +91,86 @@
"crashed (not enough memory)",
"exploded (fatal coding error: improper input parameters)",
"stopped (break requested within function evaluation)",
- "found nan (function value is not-a-number or infinite)"};
+ "found nan (function value is not-a-number or infinite)",
+ "won't fit (no free parameter)"
+};
-/******************************************************************************/
-/* Monitoring auxiliaries. */
-/******************************************************************************/
+const char *lm_shortmsg[] = {
+ "found zero", // 0
+ "converged (f)", // 1
+ "converged (p)", // 2
+ "converged (2)", // 3
+ "degenerate", // 4
+ "call limit", // 5
+ "failed (f)", // 6
+ "failed (p)", // 7
+ "failed (o)", // 8
+ "no memory", // 9
+ "invalid input", // 10
+ "user break", // 11
+ "found nan", // 12
+ "no free par" // 13
+};
-static void lm_print_pars(int nout, const double* par, FILE* fout)
+
+/*****************************************************************************/
+/* Monitoring auxiliaries. */
+/*****************************************************************************/
+
+static void lm_print_pars(const int nout, const double *par, FILE* fout)
{
- int i;
- for (i = 0; i < nout; ++i)
- fprintf(fout, " %16.9g", par[i]);
+ fprintf(fout, " pars:");
+ for (int i = 0; i < nout; ++i)
+ fprintf(fout, " %23.16g", par[i]);
fprintf(fout, "\n");
}
-/******************************************************************************/
-/* lmmin (main minimization routine) */
-/******************************************************************************/
-void lmmin(const int n, double* x, const int m, const void* data,
- void (*evaluate)(const double* par, const int m_dat,
- const void* data, double* fvec, int* userbreak),
- const lm_control_struct* C, lm_status_struct* S)
+/*****************************************************************************/
+/* lmmin (main minimization routine) */
+/*****************************************************************************/
+
+void lmmin(
+ const int n, double *const x, const int m, const double* y,
+ const void *const data,
+ void (*const evaluate)(
+ const double *const par, const int m_dat, const void *const data,
+ double *const fvec, int *const userbreak),
+ const lm_control_struct *const C, lm_status_struct *const S)
{
int j, i;
- double actred, dirder, fnorm, fnorm1, gnorm, pnorm, prered, ratio, step,
- sum, temp, temp1, temp2, temp3;
-
- /*** Initialize internal variables. ***/
+ double actred, dirder, fnorm, fnorm1, gnorm, pnorm,
+ prered, ratio, step, sum, temp, temp1, temp2, temp3;
+ static double p1 = 0.1, p0001 = 1.0e-4;
int maxfev = C->patience * (n+1);
- int inner_success; /* flag for loop control */
- double lmpar = 0; /* Levenberg-Marquardt parameter */
+ int inner_success; /* flag for loop control */
+ double lmpar = 0; /* Levenberg-Marquardt parameter */
double delta = 0;
double xnorm = 0;
double eps = sqrt(MAX(C->epsilon, LM_MACHEP)); /* for forward differences */
- int nout = C->n_maxpri == -1 ? n : MIN(C->n_maxpri, n);
+ int nout = C->n_maxpri==-1 ? n : MIN(C->n_maxpri, n);
- /* Reinterpret C->msgfile=NULL as stdout (which is unavailable for
- compile-time initialization of lm_control_double and similar). */
+ /* The workaround msgfile=NULL is needed for default initialization */
FILE* msgfile = C->msgfile ? C->msgfile : stdout;
- /*** Default status info; must be set before first return statement. ***/
-
- S->outcome = 0; /* status code */
+ /* Default status info; must be set ahead of first return statements */
+ S->outcome = 0; /* status code */
S->userbreak = 0;
- S->nfev = 0; /* function evaluation counter */
+ S->nfev = 0; /* function evaluation counter */
- /*** Check input parameters for errors. ***/
+/*** Check input parameters for errors. ***/
- if (n <= 0) {
+ if ( n < 0 ) {
fprintf(stderr, "lmmin: invalid number of parameters %i\n", n);
- S->outcome = 10;
+ S->outcome = 10; /* invalid parameter */
return;
}
if (m < n) {
fprintf(stderr, "lmmin: number of data points (%i) "
- "smaller than number of parameters (%i)\n",
- m, n);
+ "smaller than number of parameters (%i)\n", m, n);
S->outcome = 10;
return;
}
@@ -173,93 +194,99 @@
}
if (C->scale_diag != 0 && C->scale_diag != 1) {
fprintf(stderr, "lmmin: logical variable scale_diag=%i, "
- "should be 0 or 1\n",
- C->scale_diag);
+ "should be 0 or 1\n", C->scale_diag);
S->outcome = 10;
return;
}
- /*** Allocate work space. ***/
+/*** Allocate work space. ***/
/* Allocate total workspace with just one system call */
- char* ws;
- if ((ws = (char *)malloc((2*m + 5*n + m*n) * sizeof(double) +
- n * sizeof(int))) == NULL) {
+ char *ws;
+ if ( ( ws = static_cast<char *>(malloc(
+ (2*m+5*n+m*n)*sizeof(double) + n*sizeof(int)) ) ) == NULL ) {
S->outcome = 9;
return;
}
/* Assign workspace segments. */
- char* pws = ws;
- double* fvec = (double*)pws;
- pws += m * sizeof(double) / sizeof(char);
- double* diag = (double*)pws;
- pws += n * sizeof(double) / sizeof(char);
- double* qtf = (double*)pws;
- pws += n * sizeof(double) / sizeof(char);
- double* fjac = (double*)pws;
- pws += n * m * sizeof(double) / sizeof(char);
- double* wa1 = (double*)pws;
- pws += n * sizeof(double) / sizeof(char);
- double* wa2 = (double*)pws;
- pws += n * sizeof(double) / sizeof(char);
- double* wa3 = (double*)pws;
- pws += n * sizeof(double) / sizeof(char);
- double* wf = (double*)pws;
- pws += m * sizeof(double) / sizeof(char);
- int* Pivot = (int*)pws;
- //pws += n * sizeof(int) / sizeof(char);
+ char *pws = ws;
+ double *fvec = (double*) pws; pws += m * sizeof(double)/sizeof(char);
+ double *diag = (double*) pws; pws += n * sizeof(double)/sizeof(char);
+ double *qtf = (double*) pws; pws += n * sizeof(double)/sizeof(char);
+ double *fjac = (double*) pws; pws += n*m*sizeof(double)/sizeof(char);
+ double *wa1 = (double*) pws; pws += n * sizeof(double)/sizeof(char);
+ double *wa2 = (double*) pws; pws += n * sizeof(double)/sizeof(char);
+ double *wa3 = (double*) pws; pws += n * sizeof(double)/sizeof(char);
+ double *wf = (double*) pws; pws += m * sizeof(double)/sizeof(char);
+ int *ipvt = (int*) pws; pws += n * sizeof(int) /sizeof(char);
- /* Initialize diag. */
- if (!C->scale_diag)
+ /* Initialize diag */ // TODO: check whether this is still needed
+ if (!C->scale_diag) {
for (j = 0; j < n; j++)
- diag[j] = 1;
-
- /*** Evaluate function at starting point and calculate norm. ***/
-
- if (C->verbosity) {
- fprintf(msgfile, "lmmin start ");
- lm_print_pars(nout, x, msgfile);
+ diag[j] = 1.;
}
- (*evaluate)(x, m, data, fvec, &(S->userbreak));
- if (C->verbosity > 4)
- for (i = 0; i < m; ++i)
- fprintf(msgfile, " fvec[%4i] = %18.8g\n", i, fvec[i]);
- S->nfev = 1;
- if (S->userbreak)
- goto terminate;
- fnorm = lm_enorm(m, fvec);
- if (C->verbosity)
- fprintf(msgfile, " fnorm = %18.8g\n", fnorm);
- if (!std::isfinite(fnorm)) {
+/*** Evaluate function at starting point and calculate norm. ***/
+
+ if( C->verbosity&1 )
+ fprintf(msgfile, "lmmin start (ftol=%g gtol=%g xtol=%g)\n",
+ C->ftol, C->gtol, C->xtol);
+ if( C->verbosity&2 )
+ lm_print_pars(nout, x, msgfile);
+ (*evaluate)(x, m, data, fvec, &(S->userbreak));
+ if( C->verbosity&8 )
+ {
+ if (y) {
+ for( i=0; i<m; ++i )
+ fprintf(msgfile, " i, f, y-f: %4i %18.8g %18.8g\n",
+ i, fvec[i], y[i]-fvec[i]);
+ } else {
+ for( i=0; i<m; ++i )
+ fprintf(msgfile, " i, f: %4i %18.8g\n", i, fvec[i]);
+ }
+ }
+ S->nfev = 1;
+ if ( S->userbreak )
+ goto terminate;
+ if ( n == 0 ) {
+ S->outcome = 13; /* won't fit */
+ goto terminate;
+ }
+ fnorm = lm_fnorm(m, fvec, y);
+ if( C->verbosity&2 )
+ fprintf(msgfile, " fnorm = %24.16g\n", fnorm);
+ if( !isfinite(fnorm) ){
+ if( C->verbosity )
+ fprintf(msgfile, "nan case 1\n");
S->outcome = 12; /* nan */
goto terminate;
- } else if (fnorm <= LM_DWARF) {
+ } else if( fnorm <= LM_DWARF ){
S->outcome = 0; /* sum of squares almost zero, nothing to do */
goto terminate;
}
- /*** The outer loop: compute gradient, then descend. ***/
+/*** The outer loop: compute gradient, then descend. ***/
- for (int outer = 0;; ++outer) {
+ for( int outer=0; ; ++outer ) {
- /** Calculate the Jacobian. **/
+/*** [outer] Calculate the Jacobian. ***/
+
for (j = 0; j < n; j++) {
temp = x[j];
- step = MAX(eps * eps, eps * fabs(temp));
+ step = MAX(eps*eps, eps * fabs(temp));
x[j] += step; /* replace temporarily */
(*evaluate)(x, m, data, wf, &(S->userbreak));
++(S->nfev);
- if (S->userbreak)
+ if ( S->userbreak )
goto terminate;
for (i = 0; i < m; i++)
fjac[j*m+i] = (wf[i] - fvec[i]) / step;
x[j] = temp; /* restore */
}
- if (C->verbosity >= 10) {
+ if ( C->verbosity&16 ) {
/* print the entire matrix */
- printf("\nlmmin Jacobian\n");
+ printf("Jacobian\n");
for (i = 0; i < m; i++) {
printf(" ");
for (j = 0; j < n; j++)
@@ -268,34 +295,39 @@
}
}
- /** Compute the QR factorization of the Jacobian. **/
+/*** [outer] Compute the QR factorization of the Jacobian. ***/
- /* fjac is an m by n array. The upper n by n submatrix of fjac is made
- * to contain an upper triangular matrix R with diagonal elements of
- * nonincreasing magnitude such that
- *
- * P^T*(J^T*J)*P = R^T*R
- *
- * (NOTE: ^T stands for matrix transposition),
- *
- * where P is a permutation matrix and J is the final calculated
- * Jacobian. Column j of P is column Pivot(j) of the identity matrix.
- * The lower trapezoidal part of fjac contains information generated
- * during the computation of R.
- *
- * Pivot is an integer array of length n. It defines a permutation
- * matrix P such that jac*P = Q*R, where jac is the final calculated
- * Jacobian, Q is orthogonal (not stored), and R is upper triangular
- * with diagonal elements of nonincreasing magnitude. Column j of P
- * is column Pivot(j) of the identity matrix.
- */
+/* fjac is an m by n array. The upper n by n submatrix of fjac
+ * is made to contain an upper triangular matrix R with diagonal
+ * elements of nonincreasing magnitude such that
+ *
+ * P^T*(J^T*J)*P = R^T*R
+ *
+ * (NOTE: ^T stands for matrix transposition),
+ *
+ * where P is a permutation matrix and J is the final calculated
+ * Jacobian. Column j of P is column ipvt(j) of the identity matrix.
+ * The lower trapezoidal part of fjac contains information generated
+ * during the computation of R.
+ *
+ * ipvt is an integer array of length n. It defines a permutation
+ * matrix P such that jac*P = Q*R, where jac is the final calculated
+ * Jacobian, Q is orthogonal (not stored), and R is upper triangular
+ * with diagonal elements of nonincreasing magnitude. Column j of P
+ * is column ipvt(j) of the identity matrix.
+ */
- lm_qrfac(m, n, fjac, Pivot, wa1, wa2, wa3);
- /* return values are Pivot, wa1=rdiag, wa2=acnorm */
+ lm_qrfac(m, n, fjac, ipvt, wa1, wa2, wa3);
+ /* return values are ipvt, wa1=rdiag, wa2=acnorm */
- /** Form Q^T * fvec, and store first n components in qtf. **/
- for (i = 0; i < m; i++)
- wf[i] = fvec[i];
+/*** [outer] Form Q^T * fvec, and store first n components in qtf. ***/
+
+ if (y)
+ for (i = 0; i < m; i++)
+ wf[i] = fvec[i] - y[i];
+ else
+ for (i = 0; i < m; i++)
+ wf[i] = fvec[i];
for (j = 0; j < n; j++) {
temp3 = fjac[j*m+j];
@@ -311,15 +343,16 @@
qtf[j] = wf[j];
}
- /** Compute norm of scaled gradient and detect degeneracy. **/
+/*** [outer] Compute norm of scaled gradient and detect degeneracy. ***/
+
gnorm = 0;
for (j = 0; j < n; j++) {
- if (wa2[Pivot[j]] == 0)
+ if (wa2[ipvt[j]] == 0)
continue;
sum = 0;
for (i = 0; i <= j; i++)
sum += fjac[j*m+i] * qtf[i];
- gnorm = MAX(gnorm, fabs(sum / wa2[Pivot[j]] / fnorm));
+ gnorm = MAX(gnorm, fabs( sum / wa2[ipvt[j]] / fnorm ));
}
if (gnorm <= C->gtol) {
@@ -327,8 +360,10 @@
goto terminate;
}
- /** Initialize or update diag and delta. **/
- if (!outer) { /* first iteration only */
+/*** [outer] Initialize / update diag and delta. ***/
+
+ if ( !outer ) {
+ /* first iteration only */
if (C->scale_diag) {
/* diag := norms of the columns of the initial Jacobian */
for (j = 0; j < n; j++)
@@ -337,129 +372,140 @@
for (j = 0; j < n; j++)
wa3[j] = diag[j] * x[j];
xnorm = lm_enorm(n, wa3);
- if (C->verbosity >= 2) {
- fprintf(msgfile, "lmmin diag ");
- lm_print_pars(nout, x, msgfile); // xnorm
- fprintf(msgfile, " xnorm = %18.8g\n", xnorm);
- }
- /* Only now print the header for the loop table. */
- if (C->verbosity >= 3) {
- fprintf(msgfile, " o i lmpar prered"
- " ratio dirder delta"
- " pnorm fnorm");
- for (i = 0; i < nout; ++i)
- fprintf(msgfile, " p%i", i);
- fprintf(msgfile, "\n");
- }
} else {
xnorm = lm_enorm(n, x);
}
- if (!std::isfinite(xnorm)) {
+ if( !isfinite(xnorm) ){
+ if( C->verbosity )
+ fprintf(msgfile, "nan case 2\n");
S->outcome = 12; /* nan */
goto terminate;
}
- /* Initialize the step bound delta. */
- if (xnorm)
+ /* initialize the step bound delta. */
+ if ( xnorm )
delta = C->stepbound * xnorm;
else
delta = C->stepbound;
+ /* only now print the header for the loop table */
+ if( C->verbosity&2 ) {
+ fprintf(msgfile, " #o #i lmpar prered actred"
+ " ratio dirder delta"
+ " pnorm fnorm");
+ for (i = 0; i < nout; ++i)
+ fprintf(msgfile, " p%i", i);
+ fprintf(msgfile, "\n");
+ }
} else {
if (C->scale_diag) {
for (j = 0; j < n; j++)
- diag[j] = MAX(diag[j], wa2[j]);
+ diag[j] = MAX( diag[j], wa2[j] );
}
}
- /** The inner loop. **/
+/*** The inner loop. ***/
int inner = 0;
do {
- /** Determine the Levenberg-Marquardt parameter. **/
- lm_lmpar(n, fjac, m, Pivot, diag, qtf, delta, &lmpar,
+/*** [inner] Determine the Levenberg-Marquardt parameter. ***/
+
+ lm_lmpar(n, fjac, m, ipvt, diag, qtf, delta, &lmpar,
wa1, wa2, wf, wa3);
/* used return values are fjac (partly), lmpar, wa1=x, wa3=diag*x */
- /* Predict scaled reduction. */
+ /* predict scaled reduction */
pnorm = lm_enorm(n, wa3);
- if (!std::isfinite(pnorm)) {
+ if( !isfinite(pnorm) ){
+ if( C->verbosity )
+ fprintf(msgfile, "nan case 3\n");
S->outcome = 12; /* nan */
goto terminate;
}
- temp2 = lmpar * SQR(pnorm / fnorm);
+ temp2 = lmpar * SQR( pnorm / fnorm );
for (j = 0; j < n; j++) {
wa3[j] = 0;
for (i = 0; i <= j; i++)
- wa3[i] -= fjac[j*m+i] * wa1[Pivot[j]];
+ wa3[i] -= fjac[j*m+i] * wa1[ipvt[j]];
}
- temp1 = SQR(lm_enorm(n, wa3) / fnorm);
- if (!std::isfinite(temp1)) {
+ temp1 = SQR( lm_enorm(n, wa3) / fnorm );
+ if( !isfinite(temp1) ){
+ if( C->verbosity )
+ fprintf(msgfile, "nan case 4\n");
S->outcome = 12; /* nan */
goto terminate;
}
- prered = temp1 + 2*temp2;
+ prered = temp1 + 2 * temp2;
dirder = -temp1 + temp2; /* scaled directional derivative */
- /* At first call, adjust the initial step bound. */
- if (!outer && pnorm < delta)
+ /* at first call, adjust the initial step bound. */
+ if ( !outer && !inner && pnorm < delta )
delta = pnorm;
- /** Evaluate the function at x + p. **/
+/*** [inner] Evaluate the function at x + p. ***/
+
for (j = 0; j < n; j++)
wa2[j] = x[j] - wa1[j];
- (*evaluate)(wa2, m, data, wf, &(S->userbreak));
+
+ (*evaluate)( wa2, m, data, wf, &(S->userbreak) );
++(S->nfev);
- if (S->userbreak)
+ if ( S->userbreak )
goto terminate;
- fnorm1 = lm_enorm(m, wf);
- if (!std::isfinite(fnorm1)) {
- S->outcome = 12; /* nan */
- goto terminate;
+ fnorm1 = lm_fnorm(m, wf, y);
+ // exceptionally, for this norm we do not test for infinity
+ // because we can deal with it without terminating.
+
+/*** [inner] Evaluate the scaled reduction. ***/
+
+ /* actual scaled reduction (supports even the case fnorm1=infty) */
+ if (p1 * fnorm1 < fnorm)
+ actred = 1 - SQR(fnorm1 / fnorm);
+ else
+ actred = -1;
+
+ /* ratio of actual to predicted reduction */
+ ratio = prered ? actred/prered : 0;
+
+ if( C->verbosity&32 ) {
+ if (y) {
+ for( i=0; i<m; ++i )
+ fprintf(msgfile, " i, f, y-f: %4i %18.8g %18.8g\n",
+ i, fvec[i], y[i]-fvec[i]);
+ } else {
+ for( i=0; i<m; ++i )
+ fprintf(msgfile, " i, f, y-f: %4i %18.8g\n",
+ i, fvec[i]);
+ }
}
-
- /** Evaluate the scaled reduction. **/
-
- /* Actual scaled reduction. */
- actred = 1 - SQR(fnorm1 / fnorm);
-
- /* Ratio of actual to predicted reduction. */
- ratio = prered ? actred / prered : 0;
-
- if (C->verbosity == 2) {
- fprintf(msgfile, "lmmin (%i:%i) ", outer, inner);
- lm_print_pars(nout, wa2, msgfile); // fnorm1,
- } else if (C->verbosity >= 3) {
- printf("%3i %2i %9.2g %9.2g %14.6g"
+ if( C->verbosity&2 ) {
+ printf("%3i %2i %9.2g %9.2g %9.2g %14.6g"
" %9.2g %10.3e %10.3e %21.15e",
- outer, inner, lmpar, prered, ratio,
+ outer, inner, lmpar, prered, actred, ratio,
dirder, delta, pnorm, fnorm1);
for (i = 0; i < nout; ++i)
fprintf(msgfile, " %16.9g", wa2[i]);
fprintf(msgfile, "\n");
}
- /* Update the step bound. */
- if (ratio <= 0.25) {
- if (actred >= 0)
- temp = 0.5;
- else if (actred > -99) /* -99 = 1-1/0.1^2 */
- temp = MAX(dirder / (2*dirder + actred), 0.1);
- else
- temp = 0.1;
- delta = temp * MIN(delta, pnorm / 0.1);
- lmpar /= temp;
- } else if (ratio >= 0.75) {
- delta = 2 * pnorm;
- lmpar *= 0.5;
- } else if (!lmpar) {
- delta = 2 * pnorm;
- }
+ /* update the step bound */
+ if (ratio <= 0.25) {
+ if (actred >= 0)
+ temp = 0.5;
+ else
+ temp = 0.5 * dirder / (dirder + 0.5 * actred);
+ if (p1 * fnorm1 >= fnorm || temp < p1)
+ temp = p1;
+ delta = temp * MIN(delta, pnorm / p1);
+ lmpar /= temp;
+ } else if (lmpar == 0 || ratio >= 0.75) {
+ delta = 2 * pnorm;
+ lmpar *= 0.5;
+ }
- /** On success, update solution, and test for convergence. **/
+/*** [inner] On success, update solution, and test for convergence. ***/
- inner_success = ratio >= 1e-4;
- if (inner_success) {
+ inner_success = ratio >= p0001;
+ if ( inner_success ) {
- /* Update x, fvec, and their norms. */
+ /* update x, fvec, and their norms */
if (C->scale_diag) {
for (j = 0; j < n; j++) {
x[j] = wa2[j];
@@ -472,172 +518,193 @@
for (i = 0; i < m; i++)
fvec[i] = wf[i];
xnorm = lm_enorm(n, wa2);
- if (!std::isfinite(xnorm)) {
+ if( !isfinite(xnorm) ){
+ if( C->verbosity )
+ fprintf(msgfile, "nan case 6\n");
S->outcome = 12; /* nan */
goto terminate;
}
fnorm = fnorm1;
}
- /* Convergence tests. */
+ /* convergence tests */
S->outcome = 0;
- if (fnorm <= LM_DWARF)
- goto terminate; /* success: sum of squares almost zero */
- /* Test two criteria (both may be fulfilled). */
+ if( fnorm<=LM_DWARF )
+ goto terminate; /* success: sum of squares almost zero */
+ /* test two criteria (both may be fulfilled) */
if (fabs(actred) <= C->ftol && prered <= C->ftol && ratio <= 2)
- S->outcome = 1; /* success: x almost stable */
+ S->outcome = 1; /* success: x almost stable */
if (delta <= C->xtol * xnorm)
S->outcome += 2; /* success: sum of squares almost stable */
if (S->outcome != 0) {
goto terminate;
}
- /** Tests for termination and stringent tolerances. **/
- if (S->nfev >= maxfev) {
+/*** [inner] Tests for termination and stringent tolerances. ***/
+
+ if ( S->nfev >= maxfev ){
S->outcome = 5;
goto terminate;
}
- if (fabs(actred) <= LM_MACHEP && prered <= LM_MACHEP &&
- ratio <= 2) {
+ if ( fabs(actred) <= LM_MACHEP &&
+ prered <= LM_MACHEP && ratio <= 2 ){
S->outcome = 6;
goto terminate;
}
- if (delta <= LM_MACHEP * xnorm) {
+ if ( delta <= LM_MACHEP*xnorm ){
S->outcome = 7;
goto terminate;
}
- if (gnorm <= LM_MACHEP) {
+ if ( gnorm <= LM_MACHEP ){
S->outcome = 8;
goto terminate;
}
- /** End of the inner loop. Repeat if iteration unsuccessful. **/
- ++inner;
- } while (!inner_success);
+/*** [inner] End of the loop. Repeat if iteration unsuccessful. ***/
- }; /*** End of the outer loop. ***/
+ ++inner;
+ } while ( !inner_success );
+
+/*** [outer] End of the loop. ***/
+
+ };
terminate:
- S->fnorm = lm_enorm(m, fvec);
- if (C->verbosity >= 2)
- printf("lmmin outcome (%i) xnorm %g ftol %g xtol %g\n", S->outcome,
- xnorm, C->ftol, C->xtol);
- if (C->verbosity & 1) {
- fprintf(msgfile, "lmmin final ");
- lm_print_pars(nout, x, msgfile); // S->fnorm,
- fprintf(msgfile, " fnorm = %18.8g\n", S->fnorm);
+ S->fnorm = lm_fnorm(m, fvec, y);
+ if( C->verbosity&1 )
+ fprintf(msgfile, "lmmin terminates with outcome %i\n", S->outcome);
+ if( C->verbosity&2 )
+ lm_print_pars(nout, x, msgfile);
+ if( C->verbosity&8 ) {
+ if (y) {
+ for( i=0; i<m; ++i )
+ fprintf(msgfile, " i, f, y-f: %4i %18.8g %18.8g\n",
+ i, fvec[i], y[i]-fvec[i] );
+ } else {
+ for( i=0; i<m; ++i )
+ fprintf(msgfile, " i, f, y-f: %4i %18.8g\n", i, fvec[i]);
+ }
}
- if (S->userbreak) /* user-requested break */
+ if( C->verbosity&2 )
+ fprintf(msgfile, " fnorm=%24.16g xnorm=%24.16g\n", S->fnorm, xnorm);
+ if ( S->userbreak ) /* user-requested break */
S->outcome = 11;
- /*** Deallocate the workspace. ***/
+/*** Deallocate the workspace. ***/
free(ws);
} /*** lmmin. ***/
-/******************************************************************************/
-/* lm_lmpar (determine Levenberg-Marquardt parameter) */
-/******************************************************************************/
-static void lm_lmpar(const int n, double* r, const int ldr, const int* Pivot,
- const double* diag, const double* qtb, const double delta,
- double* par, double* x, double* Sdiag, double* aux, double* xdi)
+/*****************************************************************************/
+/* lm_lmpar (determine Levenberg-Marquardt parameter) */
+/*****************************************************************************/
+
+void lm_lmpar(
+ const int n, double *const r, const int ldr, int *const ipvt,
+ double *const diag, double *const qtb, double delta, double *const par,
+ double *const x, double *const sdiag, double *const aux, double *const xdi)
+{
/* Given an m by n matrix A, an n by n nonsingular diagonal matrix D,
* an m-vector b, and a positive number delta, the problem is to
* determine a parameter value par such that if x solves the system
*
* A*x = b and sqrt(par)*D*x = 0
*
- * in the least squares sense, and dxnorm is the Euclidean norm of D*x,
- * then either par=0 and (dxnorm-delta) < 0.1*delta, or par>0 and
- * abs(dxnorm-delta) < 0.1*delta.
+ * in the least squares sense, and dxnorm is the euclidean
+ * norm of D*x, then either par=0 and (dxnorm-delta) < 0.1*delta,
+ * or par>0 and abs(dxnorm-delta) < 0.1*delta.
*
* Using lm_qrsolv, this subroutine completes the solution of the
- * problem if it is provided with the necessary information from the
- * QR factorization, with column pivoting, of A. That is, if A*P = Q*R,
- * where P is a permutation matrix, Q has orthogonal columns, and R is
- * an upper triangular matrix with diagonal elements of nonincreasing
- * magnitude, then lmpar expects the full upper triangle of R, the
- * permutation matrix P, and the first n components of Q^T*b. On output
- * lmpar also provides an upper triangular matrix S such that
+ * problem if it is provided with the necessary information from
+ * the QR factorization, with column pivoting, of A. That is, if
+ * A*P = Q*R, where P is a permutation matrix, Q has orthogonal
+ * columns, and R is an upper triangular matrix with diagonal
+ * elements of nonincreasing magnitude, then lmpar expects the
+ * full upper triangle of R, the permutation matrix P, and the
+ * first n components of Q^T*b. On output lmpar also provides an
+ * upper triangular matrix S such that
*
* P^T*(A^T*A + par*D*D)*P = S^T*S.
*
* S is employed within lmpar and may be of separate interest.
*
- * Only a few iterations are generally needed for convergence of the
- * algorithm. If, however, the limit of 10 iterations is reached, then
- * the output par will contain the best value obtained so far.
+ * Only a few iterations are generally needed for convergence
+ * of the algorithm. If, however, the limit of 10 iterations
+ * is reached, then the output par will contain the best value
+ * obtained so far.
*
* Parameters:
*
* n is a positive integer INPUT variable set to the order of r.
*
- * r is an n by n array. On INPUT the full upper triangle must contain
- * the full upper triangle of the matrix R. On OUTPUT the full upper
- * triangle is unaltered, and the strict lower triangle contains the
- * strict upper triangle (transposed) of the upper triangular matrix S.
+ * r is an n by n array. On INPUT the full upper triangle
+ * must contain the full upper triangle of the matrix R.
+ * On OUTPUT the full upper triangle is unaltered, and the
+ * strict lower triangle contains the strict upper triangle
+ * (transposed) of the upper triangular matrix S.
*
- * ldr is a positive integer INPUT variable not less than n which
- * specifies the leading dimension of the array R.
+ * ldr is a positive integer INPUT variable not less than n
+ * which specifies the leading dimension of the array R.
*
- * Pivot is an integer INPUT array of length n which defines the
- * permutation matrix P such that A*P = Q*R. Column j of P is column
- * Pivot(j) of the identity matrix.
+ * ipvt is an integer INPUT array of length n which defines the
+ * permutation matrix P such that A*P = Q*R. Column j of P
+ * is column ipvt(j) of the identity matrix.
*
- * diag is an INPUT array of length n which must contain the diagonal
- * elements of the matrix D.
+ * diag is an INPUT array of length n which must contain the
+ * diagonal elements of the matrix D.
*
* qtb is an INPUT array of length n which must contain the first
* n elements of the vector Q^T*b.
*
- * delta is a positive INPUT variable which specifies an upper bound
- * on the Euclidean norm of D*x.
+ * delta is a positive INPUT variable which specifies an upper
+ * bound on the euclidean norm of D*x.
*
- * par is a nonnegative variable. On INPUT par contains an initial
- * estimate of the Levenberg-Marquardt parameter. On OUTPUT par
- * contains the final estimate.
+ * par is a nonnegative variable. On INPUT par contains an
+ * initial estimate of the Levenberg-Marquardt parameter.
+ * On OUTPUT par contains the final estimate.
*
- * x is an OUTPUT array of length n which contains the least-squares
- * solution of the system A*x = b, sqrt(par)*D*x = 0, for the output par.
+ * x is an OUTPUT array of length n which contains the least
+ * squares solution of the system A*x = b, sqrt(par)*D*x = 0,
+ * for the output par.
*
- * Sdiag is an array of length n needed as workspace; on OUTPUT it
- * contains the diagonal elements of the upper triangular matrix S.
+ * sdiag is an array of length n needed as workspace; on OUTPUT
+ * it contains the diagonal elements of the upper triangular
+ * matrix S.
*
* aux is a multi-purpose work array of length n.
*
* xdi is a work array of length n. On OUTPUT: diag[j] * x[j].
*
*/
-{
int i, iter, j, nsing;
double dxnorm, fp, fp_old, gnorm, parc, parl, paru;
double sum, temp;
static double p1 = 0.1;
- /*** Compute and store in x the Gauss-Newton direction. If the Jacobian
- is rank-deficient, obtain a least-squares solution. ***/
+/*** lmpar: compute and store in x the gauss-newton direction. if the
+ jacobian is rank-deficient, obtain a least squares solution. ***/
nsing = n;
for (j = 0; j < n; j++) {
aux[j] = qtb[j];
- if (r[j*ldr+j] == 0 && nsing == n)
+ if (r[j * ldr + j] == 0 && nsing == n)
nsing = j;
if (nsing < n)
aux[j] = 0;
}
- for (j = nsing-1; j >= 0; j--) {
- aux[j] = aux[j] / r[j+ldr*j];
+ for (j = nsing - 1; j >= 0; j--) {
+ aux[j] = aux[j] / r[j + ldr * j];
temp = aux[j];
for (i = 0; i < j; i++)
- aux[i] -= r[j*ldr+i] * temp;
+ aux[i] -= r[j * ldr + i] * temp;
}
for (j = 0; j < n; j++)
- x[Pivot[j]] = aux[j];
+ x[ipvt[j]] = aux[j];
- /*** Initialize the iteration counter, evaluate the function at the origin,
- and test for acceptance of the Gauss-Newton direction. ***/
+/*** lmpar: initialize the iteration counter, evaluate the function at the
+ origin, and test for acceptance of the gauss-newton direction. ***/
for (j = 0; j < n; j++)
xdi[j] = diag[j] * x[j];
@@ -645,65 +712,67 @@
fp = dxnorm - delta;
if (fp <= p1 * delta) {
#ifdef LMFIT_DEBUG_MESSAGES
- printf("debug lmpar nsing=%d, n=%d, terminate[fp<=p1*del]\n", nsing, n);
+ printf("debug lmpar nsing %d n %d, terminate (fp<p1*delta)\n",
+ nsing, n);
#endif
*par = 0;
return;
}
- /*** If the Jacobian is not rank deficient, the Newton step provides a
- lower bound, parl, for the zero of the function. Otherwise set this
- bound to zero. ***/
+/*** lmpar: if the jacobian is not rank deficient, the newton
+ step provides a lower bound, parl, for the zero of
+ the function. otherwise set this bound to zero. ***/
parl = 0;
if (nsing >= n) {
for (j = 0; j < n; j++)
- aux[j] = diag[Pivot[j]] * xdi[Pivot[j]] / dxnorm;
+ aux[j] = diag[ipvt[j]] * xdi[ipvt[j]] / dxnorm;
for (j = 0; j < n; j++) {
sum = 0;
for (i = 0; i < j; i++)
- sum += r[j*ldr+i] * aux[i];
- aux[j] = (aux[j] - sum) / r[j+ldr*j];
+ sum += r[j * ldr + i] * aux[i];
+ aux[j] = (aux[j] - sum) / r[j + ldr * j];
}
temp = lm_enorm(n, aux);
parl = fp / delta / temp / temp;
}
- /*** Calculate an upper bound, paru, for the zero of the function. ***/
+/*** lmpar: calculate an upper bound, paru, for the zero of the function. ***/
for (j = 0; j < n; j++) {
sum = 0;
for (i = 0; i <= j; i++)
- sum += r[j*ldr+i] * qtb[i];
- aux[j] = sum / diag[Pivot[j]];
+ sum += r[j * ldr + i] * qtb[i];
+ aux[j] = sum / diag[ipvt[j]];
}
gnorm = lm_enorm(n, aux);
paru = gnorm / delta;
if (paru == 0)
paru = LM_DWARF / MIN(delta, p1);
- /*** If the input par lies outside of the interval (parl,paru),
- set par to the closer endpoint. ***/
+/*** lmpar: if the input par lies outside of the interval (parl,paru),
+ set par to the closer endpoint. ***/
*par = MAX(*par, parl);
*par = MIN(*par, paru);
if (*par == 0)
*par = gnorm / dxnorm;
- /*** Iterate. ***/
+/*** lmpar: iterate. ***/
- for (iter = 0;; iter++) {
+ for (iter=0; ; iter++) {
- /** Evaluate the function at the current value of par. **/
+ /** evaluate the function at the current value of par. **/
+
if (*par == 0)
*par = MAX(LM_DWARF, 0.001 * paru);
temp = sqrt(*par);
for (j = 0; j < n; j++)
aux[j] = temp * diag[j];
- lm_qrsolv(n, r, ldr, Pivot, aux, qtb, x, Sdiag, xdi);
- /* return values are r, x, Sdiag */
+ lm_qrsolv(n, r, ldr, ipvt, aux, qtb, x, sdiag, xdi);
+ /* return values are r, x, sdiag */
for (j = 0; j < n; j++)
xdi[j] = diag[j] * x[j]; /* used as output */
@@ -711,49 +780,58 @@
fp_old = fp;
fp = dxnorm - delta;
- /** If the function is small enough, accept the current value
+ /** if the function is small enough, accept the current value
of par. Also test for the exceptional cases where parl
is zero or the number of iterations has reached 10. **/
- if (fabs(fp) <= p1 * delta ||
- (parl == 0 && fp <= fp_old && fp_old < 0) || iter == 10) {
+
+ if (fabs(fp) <= p1 * delta
+ || (parl == 0 && fp <= fp_old && fp_old < 0)
+ || iter == 10) {
#ifdef LMFIT_DEBUG_MESSAGES
- printf("debug lmpar nsing=%d, iter=%d, "
- "par=%.4e [%.4e %.4e], delta=%.4e, fp=%.4e\n",
+ printf("debug lmpar nsing %d iter %d "
+ "par %.4e [%.4e %.4e] delta %.4e fp %.4e\n",
nsing, iter, *par, parl, paru, delta, fp);
#endif
break; /* the only exit from the iteration. */
}
- /** Compute the Newton correction. **/
+ /** compute the Newton correction. **/
+
for (j = 0; j < n; j++)
- aux[j] = diag[Pivot[j]] * xdi[Pivot[j]] / dxnorm;
+ aux[j] = diag[ipvt[j]] * xdi[ipvt[j]] / dxnorm;
for (j = 0; j < n; j++) {
- aux[j] = aux[j] / Sdiag[j];
- for (i = j+1; i < n; i++)
- aux[i] -= r[j*ldr+i] * aux[j];
+ aux[j] = aux[j] / sdiag[j];
+ for (i = j + 1; i < n; i++)
+ aux[i] -= r[j * ldr + i] * aux[j];
}
temp = lm_enorm(n, aux);
parc = fp / delta / temp / temp;
- /** Depending on the sign of the function, update parl or paru. **/
+ /** depending on the sign of the function, update parl or paru. **/
+
if (fp > 0)
parl = MAX(parl, *par);
- else /* fp < 0 [the case fp==0 is precluded by the break condition] */
+ else if (fp < 0)
paru = MIN(paru, *par);
+ /* the case fp==0 is precluded by the break condition */
- /** Compute an improved estimate for par. **/
+ /** compute an improved estimate for par. **/
+
*par = MAX(parl, *par + parc);
+
}
} /*** lm_lmpar. ***/
-/******************************************************************************/
-/* lm_qrfac (QR factorization, from lapack) */
-/******************************************************************************/
+/*****************************************************************************/
+/* lm_qrfac (QR factorization, from lapack) */
+/*****************************************************************************/
-static void lm_qrfac(const int m, const int n, double* A, int* Pivot, double* Rdiag,
- double* Acnorm, double* W)
+void lm_qrfac(
+ const int m, const int n, double *const A, int *const Pivot,
+ double *const Rdiag, double *const Acnorm, double *const W)
+{
/*
* This subroutine uses Householder transformations with column pivoting
* to compute a QR factorization of the m by n matrix A. That is, qrfac
@@ -772,27 +850,27 @@
*
* n is an INPUT parameter set to the number of columns of A.
*
- * A is an m by n array. On INPUT, A contains the matrix for which the
- * QR factorization is to be computed. On OUTPUT the strict upper
- * trapezoidal part of A contains the strict upper trapezoidal part
- * of R, and the lower trapezoidal part of A contains a factored form
- * of Q (the non-trivial elements of the vectors w described above).
+ * A is an m by n array. On INPUT, A contains the matrix for
+ * which the QR factorization is to be computed. On OUTPUT
+ * the strict upper trapezoidal part of A contains the strict
+ * upper trapezoidal part of R, and the lower trapezoidal
+ * part of A contains a factored form of Q (the non-trivial
+ * elements of the vectors w described above).
*
* Pivot is an integer OUTPUT array of length n that describes the
- * permutation matrix P. Column j of P is column Pivot(j) of the
- * identity matrix.
+ * permutation matrix P:
+ * Column j of P is column ipvt(j) of the identity matrix.
*
- * Rdiag is an OUTPUT array of length n which contains the diagonal
- * elements of R.
+ * Rdiag is an OUTPUT array of length n which contains the
+ * diagonal elements of R.
*
- * Acnorm is an OUTPUT array of length n which contains the norms of
- * the corresponding columns of the input matrix A. If this information
- * is not needed, then Acnorm can share storage with Rdiag.
+ * Acnorm is an OUTPUT array of length n which contains the norms
+ * of the corresponding columns of the input matrix A. If this
+ * information is not needed, then Acnorm can share storage with Rdiag.
*
* W is a work array of length n.
*
*/
-{
int i, j, k, kmax;
double ajnorm, sum, temp;
@@ -802,16 +880,19 @@
/** Compute initial column norms;
initialize Pivot with identity permutation. ***/
+
for (j = 0; j < n; j++) {
W[j] = Rdiag[j] = Acnorm[j] = lm_enorm(m, &A[j*m]);
Pivot[j] = j;
}
/** Loop over columns of A. **/
- assert(n <= m);
+
+ assert( n <= m );
for (j = 0; j < n; j++) {
/** Bring the column of largest norm into the pivot position. **/
+
kmax = j;
for (k = j+1; k < n; k++)
if (Rdiag[k] > Rdiag[kmax])
@@ -834,6 +915,7 @@
/** Compute the Householder reflection vector w_j to reduce the
j-th column of A to a multiple of the j-th unit vector. **/
+
ajnorm = lm_enorm(m-j, &A[j*m+j]);
if (ajnorm == 0) {
Rdiag[j] = 0;
@@ -850,6 +932,7 @@
/** Apply the Householder transformation U_w := 1 - 2*w_j.w_j/|w_j|^2
to the remaining columns, and update the norms. **/
+
for (k = j+1; k < n; k++) {
/* Compute scalar product w_j * a_j. */
sum = 0;
@@ -866,12 +949,12 @@
/* No idea what happens here. */
if (Rdiag[k] != 0) {
temp = A[m*k+j] / Rdiag[k];
- if (fabs(temp) < 1) {
- Rdiag[k] *= sqrt(1 - SQR(temp));
+ if ( fabs(temp)<1 ) {
+ Rdiag[k] *= sqrt(1-SQR(temp));
temp = Rdiag[k] / W[k];
} else
temp = 0;
- if (temp == 0 || 0.05 * SQR(temp) <= LM_MACHEP) {
+ if ( temp == 0 || 0.05 * SQR(temp) <= LM_MACHEP ) {
Rdiag[k] = lm_enorm(m-j-1, &A[m*k+j+1]);
W[k] = Rdiag[k];
}
@@ -882,13 +965,16 @@
}
} /*** lm_qrfac. ***/
-/******************************************************************************/
-/* lm_qrsolv (linear least-squares) */
-/******************************************************************************/
-static void lm_qrsolv(const int n, double* r, const int ldr, const int* Pivot,
- const double* diag, const double* qtb, double* x,
- double* Sdiag, double* W)
+/*****************************************************************************/
+/* lm_qrsolv (linear least-squares) */
+/*****************************************************************************/
+
+void lm_qrsolv(
+ const int n, double *const r, const int ldr, int *const ipvt,
+ double *const diag, double *const qtb, double *const x,
+ double *const sdiag, double *const wa)
+{
/*
* Given an m by n matrix A, an n by n diagonal matrix D, and an
* m-vector b, the problem is to determine an x which solves the
@@ -921,145 +1007,152 @@
*
* n is a positive integer INPUT variable set to the order of R.
*
- * r is an n by n array. On INPUT the full upper triangle must contain
- * the full upper triangle of the matrix R. On OUTPUT the full upper
- * triangle is unaltered, and the strict lower triangle contains the
- * strict upper triangle (transposed) of the upper triangular matrix S.
+ * r is an n by n array. On INPUT the full upper triangle must
+ * contain the full upper triangle of the matrix R. On OUTPUT
+ * the full upper triangle is unaltered, and the strict lower
+ * triangle contains the strict upper triangle (transposed) of
+ * the upper triangular matrix S.
*
- * ldr is a positive integer INPUT variable not less than n which
- * specifies the leading dimension of the array R.
+ * ldr is a positive integer INPUT variable not less than n
+ * which specifies the leading dimension of the array R.
*
- * Pivot is an integer INPUT array of length n which defines the
- * permutation matrix P such that A*P = Q*R. Column j of P is column
- * Pivot(j) of the identity matrix.
+ * ipvt is an integer INPUT array of length n which defines the
+ * permutation matrix P such that A*P = Q*R. Column j of P
+ * is column ipvt(j) of the identity matrix.
*
- * diag is an INPUT array of length n which must contain the diagonal
- * elements of the matrix D.
+ * diag is an INPUT array of length n which must contain the
+ * diagonal elements of the matrix D.
*
* qtb is an INPUT array of length n which must contain the first
* n elements of the vector Q^T*b.
*
- * x is an OUTPUT array of length n which contains the least-squares
- * solution of the system A*x = b, D*x = 0.
+ * x is an OUTPUT array of length n which contains the least
+ * squares solution of the system A*x = b, D*x = 0.
*
- * Sdiag is an OUTPUT array of length n which contains the diagonal
- * elements of the upper triangular matrix S.
+ * sdiag is an OUTPUT array of length n which contains the
+ * diagonal elements of the upper triangular matrix S.
*
- * W is a work array of length n.
+ * wa is a work array of length n.
*
*/
-{
int i, kk, j, k, nsing;
double qtbpj, sum, temp;
double _sin, _cos, _tan, _cot; /* local variables, not functions */
- /*** Copy R and Q^T*b to preserve input and initialize S.
- In particular, save the diagonal elements of R in x. ***/
+/*** qrsolv: copy R and Q^T*b to preserve input and initialize S.
+ In particular, save the diagonal elements of R in x. ***/
for (j = 0; j < n; j++) {
for (i = j; i < n; i++)
- r[j*ldr+i] = r[i*ldr+j];
- x[j] = r[j*ldr+j];
- W[j] = qtb[j];
+ r[j * ldr + i] = r[i * ldr + j];
+ x[j] = r[j * ldr + j];
+ wa[j] = qtb[j];
}
- /*** Eliminate the diagonal matrix D using a Givens rotation. ***/
+/*** qrsolv: eliminate the diagonal matrix D using a Givens rotation. ***/
for (j = 0; j < n; j++) {
- /*** Prepare the row of D to be eliminated, locating the diagonal
- element using P from the QR factorization. ***/
+/*** qrsolv: prepare the row of D to be eliminated, locating the
+ diagonal element using P from the QR factorization. ***/
- if (diag[Pivot[j]] != 0) {
- for (k = j; k < n; k++)
- Sdiag[k] = 0;
- Sdiag[j] = diag[Pivot[j]];
+ if (diag[ipvt[j]] == 0)
+ goto L90;
+ for (k = j; k < n; k++)
+ sdiag[k] = 0;
+ sdiag[j] = diag[ipvt[j]];
- /*** The transformations to eliminate the row of D modify only
- a single element of Q^T*b beyond the first n, which is
- initially 0. ***/
+/*** qrsolv: the transformations to eliminate the row of D modify only
+ a single element of Q^T*b beyond the first n, which is initially 0. ***/
- qtbpj = 0;
- for (k = j; k < n; k++) {
+ qtbpj = 0;
+ for (k = j; k < n; k++) {
- /** Determine a Givens rotation which eliminates the
- appropriate element in the current row of D. **/
- if (Sdiag[k] == 0)
- continue;
- kk = k + ldr * k;
- if (fabs(r[kk]) < fabs(Sdiag[k])) {
- _cot = r[kk] / Sdiag[k];
- _sin = 1 / hypot(1, _cot);
- _cos = _sin * _cot;
- } else {
- _tan = Sdiag[k] / r[kk];
- _cos = 1 / hypot(1, _tan);
- _sin = _cos * _tan;
- }
+ /** determine a Givens rotation which eliminates the
+ appropriate element in the current row of D. **/
- /** Compute the modified diagonal element of R and
- the modified element of (Q^T*b,0). **/
- r[kk] = _cos * r[kk] + _sin * Sdiag[k];
- temp = _cos * W[k] + _sin * qtbpj;
- qtbpj = -_sin * W[k] + _cos * qtbpj;
- W[k] = temp;
+ if (sdiag[k] == 0)
+ continue;
+ kk = k + ldr * k;
+ if (fabs(r[kk]) < fabs(sdiag[k])) {
+ _cot = r[kk] / sdiag[k];
+ _sin = 1 / sqrt(1 + SQR(_cot));
+ _cos = _sin * _cot;
+ } else {
+ _tan = sdiag[k] / r[kk];
+ _cos = 1 / sqrt(1 + SQR(_tan));
+ _sin = _cos * _tan;
+ }
- /** Accumulate the tranformation in the row of S. **/
- for (i = k+1; i < n; i++) {
- temp = _cos * r[k*ldr+i] + _sin * Sdiag[i];
- Sdiag[i] = -_sin * r[k*ldr+i] + _cos * Sdiag[i];
- r[k*ldr+i] = temp;
- }
+ /** compute the modified diagonal element of R and
+ the modified element of (Q^T*b,0). **/
+
+ r[kk] = _cos * r[kk] + _sin * sdiag[k];
+ temp = _cos * wa[k] + _sin * qtbpj;
+ qtbpj = -_sin * wa[k] + _cos * qtbpj;
+ wa[k] = temp;
+
+ /** accumulate the tranformation in the row of S. **/
+
+ for (i = k + 1; i < n; i++) {
+ temp = _cos * r[k * ldr + i] + _sin * sdiag[i];
+ sdiag[i] = -_sin * r[k * ldr + i] + _cos * sdiag[i];
+ r[k * ldr + i] = temp;
}
}
- /** Store the diagonal element of S and restore
+ L90:
+ /** store the diagonal element of S and restore
the corresponding diagonal element of R. **/
- Sdiag[j] = r[j*ldr+j];
- r[j*ldr+j] = x[j];
+
+ sdiag[j] = r[j * ldr + j];
+ r[j * ldr + j] = x[j];
}
- /*** Solve the triangular system for z. If the system is singular, then
- obtain a least-squares solution. ***/
+/*** qrsolv: solve the triangular system for z. If the system is
+ singular, then obtain a least squares solution. ***/
nsing = n;
for (j = 0; j < n; j++) {
- if (Sdiag[j] == 0 && nsing == n)
+ if (sdiag[j] == 0 && nsing == n)
nsing = j;
if (nsing < n)
- W[j] = 0;
+ wa[j] = 0;
}
- for (j = nsing-1; j >= 0; j--) {
+ for (j = nsing - 1; j >= 0; j--) {
sum = 0;
- for (i = j+1; i < nsing; i++)
- sum += r[j*ldr+i] * W[i];
- W[j] = (W[j] - sum) / Sdiag[j];
+ for (i = j + 1; i < nsing; i++)
+ sum += r[j * ldr + i] * wa[i];
+ wa[j] = (wa[j] - sum) / sdiag[j];
}
- /*** Permute the components of z back to components of x. ***/
+/*** qrsolv: permute the components of z back to components of x. ***/
for (j = 0; j < n; j++)
- x[Pivot[j]] = W[j];
+ x[ipvt[j]] = wa[j];
} /*** lm_qrsolv. ***/
-/******************************************************************************/
-/* lm_enorm (Euclidean norm) */
-/******************************************************************************/
-static double lm_enorm(int n, const double* x)
+/*****************************************************************************/
+/* lm_enorm (Euclidean norm) */
+/*****************************************************************************/
+
+double lm_enorm(const int n, const double *const x)
+{
/* This function calculates the Euclidean norm of an n-vector x.
*
- * The Euclidean norm is computed by accumulating the sum of squares
- * in three different sums. The sums of squares for the small and large
- * components are scaled so that no overflows occur. Non-destructive
- * underflows are permitted. Underflows and overflows do not occur in
- * the computation of the unscaled sum of squares for the intermediate
- * components. The definitions of small, intermediate and large components
+ * The Euclidean norm is computed by accumulating the sum of
+ * squares in three different sums. The sums of squares for the
+ * small and large components are scaled so that no overflows
+ * occur. Non-destructive underflows are permitted. Underflows
+ * and overflows do not occur in the computation of the unscaled
+ * sum of squares for the intermediate components.
+ * The definitions of small, intermediate and large components
* depend on two constants, LM_SQRT_DWARF and LM_SQRT_GIANT. The main
- * restrictions on these constants are that LM_SQRT_DWARF**2 not underflow
- * and LM_SQRT_GIANT**2 not overflow.
+ * restrictions on these constants are that LM_SQRT_DWARF**2 not
+ * underflow and LM_SQRT_GIANT**2 not overflow.
*
* Parameters:
*
@@ -1067,9 +1160,8 @@
*
* x is an INPUT array of length n.
*/
-{
int i;
- double agiant, s1, s2, s3, xabs, x1max, x3max;
+ double agiant, s1, s2, s3, xabs, x1max, x3max, temp;
s1 = 0;
s2 = 0;
@@ -1078,27 +1170,33 @@
x3max = 0;
agiant = LM_SQRT_GIANT / n;
- /** Sum squares. **/
+ /** sum squares. **/
+
for (i = 0; i < n; i++) {
xabs = fabs(x[i]);
if (xabs > LM_SQRT_DWARF) {
- if (xabs < agiant) {
- s2 += SQR(xabs);
- } else if (xabs > x1max) {
- s1 = 1 + s1 * SQR(x1max / xabs);
+ if ( xabs < agiant ) {
+ s2 += xabs * xabs;
+ } else if ( xabs > x1max ) {
+ temp = x1max / xabs;
+ s1 = 1 + s1 * SQR(temp);
x1max = xabs;
} else {
- s1 += SQR(xabs / x1max);
+ temp = xabs / x1max;
+ s1 += SQR(temp);
}
- } else if (xabs > x3max) {
- s3 = 1 + s3 * SQR(x3max / xabs);
+ } else if ( xabs > x3max ) {
+ temp = x3max / xabs;
+ s3 = 1 + s3 * SQR(temp);
x3max = xabs;
} else if (xabs != 0) {
- s3 += SQR(xabs / x3max);
+ temp = xabs / x3max;
+ s3 += SQR(temp);
}
}
- /** Calculate the norm. **/
+ /** calculation of norm. **/
+
if (s1 != 0)
return x1max * sqrt(s1 + (s2 / x1max) / x1max);
else if (s2 != 0)
@@ -1110,3 +1208,80 @@
return x3max * sqrt(s3);
} /*** lm_enorm. ***/
+
+
+/*****************************************************************************/
+/* lm_fnorm (Euclidean norm of difference) */
+/*****************************************************************************/
+
+double lm_fnorm(const int n, const double *const x, const double *const y)
+{
+/* This function calculates the Euclidean norm of an n-vector x-y.
+ *
+ * The Euclidean norm is computed by accumulating the sum of
+ * squares in three different sums. The sums of squares for the
+ * small and large components are scaled so that no overflows
+ * occur. Non-destructive underflows are permitted. Underflows
+ * and overflows do not occur in the computation of the unscaled
+ * sum of squares for the intermediate components.
+ * The definitions of small, intermediate and large components
+ * depend on two constants, LM_SQRT_DWARF and LM_SQRT_GIANT. The main
+ * restrictions on these constants are that LM_SQRT_DWARF**2 not
+ * underflow and LM_SQRT_GIANT**2 not overflow.
+ *
+ * Parameters:
+ *
+ * n is a positive integer INPUT variable.
+ *
+ * x, y are INPUT arrays of length n.
+ */
+ if (!y)
+ return lm_enorm(n, x);
+ int i;
+ double agiant, s1, s2, s3, xabs, x1max, x3max, temp;
+
+ s1 = 0;
+ s2 = 0;
+ s3 = 0;
+ x1max = 0;
+ x3max = 0;
+ agiant = LM_SQRT_GIANT / n;
+
+ /** sum squares. **/
+
+ for (i = 0; i < n; i++) {
+ xabs = fabs(x[i]-y[i]);
+ if (xabs > LM_SQRT_DWARF) {
+ if ( xabs < agiant ) {
+ s2 += xabs * xabs;
+ } else if ( xabs > x1max ) {
+ temp = x1max / xabs;
+ s1 = 1 + s1 * SQR(temp);
+ x1max = xabs;
+ } else {
+ temp = xabs / x1max;
+ s1 += SQR(temp);
+ }
+ } else if ( xabs > x3max ) {
+ temp = x3max / xabs;
+ s3 = 1 + s3 * SQR(temp);
+ x3max = xabs;
+ } else if (xabs != 0) {
+ temp = xabs / x3max;
+ s3 += SQR(temp);
+ }
+ }
+
+ /** calculation of norm. **/
+
+ if (s1 != 0)
+ return x1max * sqrt(s1 + (s2 / x1max) / x1max);
+ else if (s2 != 0)
+ if (s2 >= x3max)
+ return sqrt(s2 * (1 + (x3max / s2) * (x3max * s3)));
+ else
+ return sqrt(x3max * ((s2 / x3max) + (x3max * s3)));
+ else
+ return x3max * sqrt(s3);
+
+} /*** lm_fnorm. ***/
diff --git a/src/external/lmfit/lmmin.h b/src/external/lmfit/lmmin.h
index 17b4880..d2ccfae 100644
--- a/src/external/lmfit/lmmin.h
+++ b/src/external/lmfit/lmmin.h
@@ -17,9 +17,15 @@
#include "lmstruct.h"
-/*! \brief
- * Levenberg-Marquardt minimization.
- *
+/* Levenberg-Marquardt minimization. */
+void lmmin(
+ const int n_par, double* par, const int m_dat, const double* y,
+ const void* data,
+ void (*evaluate)(
+ const double* par, const int m_dat, const void* data,
+ double* fvec, int* userbreak),
+ const lm_control_struct* control, lm_status_struct* status);
+/*
* This routine contains the core algorithm of our library.
*
* It minimizes the sum of the squares of m nonlinear functions
@@ -29,14 +35,16 @@
*
* Parameters:
*
- * \param[in] n_par The number of variables (INPUT, positive integer).
- * \param par The parameters to be fitted
- * x is the solution vector (INPUT/OUTPUT, array of length n).
+ * n_par is the number of variables (INPUT, positive integer).
+ *
+ * par is the solution vector (INPUT/OUTPUT, array of length n).
* On input it must be set to an estimated solution.
* On output it yields the final estimate of the solution.
*
- * m is the number of functions to be minimized (INPUT, positive integer).
+ * m_dat is the number of functions to be minimized (INPUT, positive integer).
* It must fulfill m>=n.
+ *
+ * y contains data to be fitted. Use a null pointer if there are no data.
*
* data is a pointer that is ignored by lmmin; it is however forwarded
* to the user-supplied functions evaluate and printout.
@@ -56,9 +64,9 @@
* status contains OUTPUT variables that inform about the fit result,
* as declared and explained in lmstruct.h
*/
-void lmmin( const int n_par, double *par, const int m_dat, const void *data,
- void (*evaluate) (const double *par, const int m_dat, const void *data,
- double *fvec, int *userbreak),
- const lm_control_struct *control, lm_status_struct *status );
+
+/* Refined calculation of Eucledian norm. */
+double lm_enorm(const int, const double*);
+double lm_fnorm(const int, const double*, const double*);
#endif /* LMMIN_H */
diff --git a/src/external/lmfit/lmstruct.h b/src/external/lmfit/lmstruct.h
index d545ff1..1f513b3 100644
--- a/src/external/lmfit/lmstruct.h
+++ b/src/external/lmfit/lmstruct.h
@@ -39,23 +39,23 @@
stepbound itself. In most cases stepbound should lie
in the interval (0.1,100.0). Generally, the value
100.0 is recommended. */
- int patience; /* Used to set the maximum number of function evaluations
+ int patience; /* Used to set the maximum number of function evaluations
to patience*(number_of_parameters+1). */
- int scale_diag; /* If 1, the variables will be rescaled internally.
+ int scale_diag; /* If 1, the variables will be rescaled internally.
Recommended value is 1. */
FILE* msgfile; /* Progress messages will be written to this file. */
- int verbosity; /* OR'ed: 1: print some messages; 2: print Jacobian. */
- int n_maxpri; /* -1, or max number of parameters to print. */
- int m_maxpri; /* -1, or max number of residuals to print. */
+ int verbosity; /* OR'ed: 1: print some messages; 2: print Jacobian. */
+ int n_maxpri; /* -1, or max number of parameters to print. */
+ int m_maxpri; /* -1, or max number of residuals to print. */
} lm_control_struct;
/* Collection of output parameters for status info. */
typedef struct {
- double fnorm; /* norm of the residue vector fvec. */
- int nfev; /* actual number of iterations. */
- int outcome; /* Status indicator. Nonnegative values are used as index
- for the message text lm_infmsg, set in lmmin.c. */
- int userbreak; /* Set when function evaluation requests termination. */
+ double fnorm; /* norm of the residue vector fvec. */
+ int nfev; /* actual number of iterations. */
+ int outcome; /* Status indicator. Nonnegative values are used as index
+ for the message text lm_infmsg, set in lmmin.c. */
+ int userbreak; /* Set when function evaluation requests termination. */
} lm_status_struct;
/* Preset (and recommended) control parameter settings. */
diff --git a/src/gromacs/CMakeLists.txt b/src/gromacs/CMakeLists.txt
index 6f1b3e9..7aab572 100644
--- a/src/gromacs/CMakeLists.txt
+++ b/src/gromacs/CMakeLists.txt
@@ -138,10 +138,6 @@
tmpi_get_source_list(THREAD_MPI_SOURCES ${CMAKE_SOURCE_DIR}/src/external/thread_mpi/src)
list(APPEND LIBGROMACS_SOURCES ${THREAD_MPI_SOURCES})
-get_lmfit_properties(LMFIT_SOURCES LMFIT_LIBRARIES_TO_LINK LMFIT_INCLUDE_DIRECTORY LMFIT_INCLUDE_DIR_ORDER)
-include_directories(${LMFIT_INCLUDE_DIR_ORDER} SYSTEM "${LMFIT_INCLUDE_DIRECTORY}")
-list(APPEND LIBGROMACS_SOURCES ${LMFIT_SOURCES})
-
configure_file(version.h.cmakein version.h)
gmx_install_headers(
analysisdata.h
@@ -260,6 +256,7 @@
VERSION ${LIBRARY_VERSION}
COMPILE_FLAGS "${OpenMP_C_FLAGS}")
+manage_lmfit()
gmx_write_installed_header_list()
# Only install the library in mdrun-only mode if it is actually necessary
diff --git a/src/gromacs/correlationfunctions/expfit.cpp b/src/gromacs/correlationfunctions/expfit.cpp
index 4fc9fd4..d2cbd4b 100644
--- a/src/gromacs/correlationfunctions/expfit.cpp
+++ b/src/gromacs/correlationfunctions/expfit.cpp
@@ -3,7 +3,7 @@
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
@@ -48,15 +48,10 @@
#include <string.h>
-#include <cmath>
-
#include <algorithm>
-
-#include <lmstruct.h>
#include "gromacs/correlationfunctions/integrate.h"
#include "gromacs/fileio/xvgr.h"
-#include "gromacs/math/functions.h"
#include "gromacs/math/vec.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/futil.h"
@@ -394,9 +389,6 @@
}
}
-/*! \brief function type for passing to fitting routine */
-typedef double (*t_lmcurve)(double x, const double *a);
-
/*! \brief array of fitting functions corresponding to the pre-defined types */
t_lmcurve lmcurves[effnNR+1] = {
lmc_exp_one_parm, lmc_exp_one_parm, lmc_exp_two_parm,
@@ -414,104 +406,6 @@
return 0.0;
}
return lmcurves[eFitFn](x, parm);
-}
-
-/*! \brief lmfit_exp supports fitting of different functions
- *
- * This routine calls the Levenberg-Marquardt non-linear fitting
- * routine for fitting a data set with errors to a target function.
- * Fitting routines included in gromacs in src/external/lmfit.
- */
-static gmx_bool lmfit_exp(int nfit,
- const double x[],
- const double y[],
- const double dy[],
- double parm[],
- gmx_bool bVerbose,
- int eFitFn,
- int nfix)
-{
- double chisq, ochisq;
- gmx_bool bCont;
- int j;
- int maxiter = 100;
- lm_control_struct control;
- lm_status_struct *status;
- int nparam = effnNparams(eFitFn);
- int p2;
- gmx_bool bSkipLast;
-
- if ((eFitFn < 0) || (eFitFn >= effnNR))
- {
- fprintf(stderr, "fitfn = %d, should be in the range 0..%d\n",
- eFitFn, effnNR-1);
- return FALSE;
- }
- /* Using default control structure for double precision fitting that
- * comes with the lmfit package (i.e. from the include file).
- */
- control = lm_control_double;
- control.verbosity = (bVerbose ? 1 : 0);
- control.n_maxpri = 0;
- control.m_maxpri = 0;
-
- snew(status, 1);
- /* Initial params */
- chisq = 1e12;
- j = 0;
- if (bVerbose)
- {
- printf("%4s %10s Parameters\n", "Step", "chi^2");
- }
- /* Check whether we have to skip some params */
- if (nfix > 0)
- {
- do
- {
- p2 = 1 << (nparam-1);
- bSkipLast = ((p2 & nfix) == p2);
- if (bSkipLast)
- {
- nparam--;
- nfix -= p2;
- }
- }
- while ((nparam > 0) && (bSkipLast));
- if (bVerbose)
- {
- printf("Using %d out of %d parameters\n", nparam, effnNparams(eFitFn));
- }
- }
- do
- {
- ochisq = chisq;
- gmx_lmcurve(nparam, parm, nfit, x, y, dy,
- lmcurves[eFitFn], &control, status);
- chisq = gmx::square(status->fnorm);
- if (bVerbose)
- {
- printf("status: fnorm = %g, nfev = %d, userbreak = %d\noutcome = %s\n",
- status->fnorm, status->nfev, status->userbreak,
- lm_infmsg[status->outcome]);
- }
- if (bVerbose)
- {
- int mmm;
- printf("%4d %8g", j, chisq);
- for (mmm = 0; (mmm < effnNparams(eFitFn)); mmm++)
- {
- printf(" %8g", parm[mmm]);
- }
- printf("\n");
- }
- j++;
- bCont = (fabs(ochisq - chisq) > fabs(control.ftol*chisq));
- }
- while (bCont && (j < maxiter));
-
- sfree(status);
-
- return TRUE;
}
/*! \brief Ensure the fitting parameters are well-behaved.
diff --git a/src/gromacs/correlationfunctions/gmx_lmcurve.cpp b/src/gromacs/correlationfunctions/gmx_lmcurve.cpp
index f3a103c..8b3870e 100644
--- a/src/gromacs/correlationfunctions/gmx_lmcurve.cpp
+++ b/src/gromacs/correlationfunctions/gmx_lmcurve.cpp
@@ -1,7 +1,7 @@
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2016, by the GROMACS development team, led by
+ * Copyright (c) 2016,2018, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
@@ -38,14 +38,30 @@
* Defines a driver routine for lmfit, and a callback for it to use.
*
* \author David van der Spoel <david.vanderspoel@icm.uu.se>
+ * \author Mark Abraham <mark.j.abraham@gmail.com>
* \ingroup module_correlationfunctions
*/
#include "gmxpre.h"
#include "gmx_lmcurve.h"
+#include "config.h"
+
+#include <cmath>
+
+#if HAVE_LMFIT
#include <lmmin.h>
#include <lmstruct.h>
+#endif
+
+#include "gromacs/correlationfunctions/expfit.h"
+#include "gromacs/math/functions.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
+
+#if HAVE_LMFIT
+
+extern t_lmcurve lmcurves[effnNR+1];
typedef struct {
const double* t;
@@ -73,7 +89,8 @@
*info = 0;
}
-void gmx_lmcurve(
+//! Calls lmmin with the given data, with callback function \c f.
+static void gmx_lmcurve(
const int n_par, double* par, const int m_dat,
const double* t, const double* y, const double *dy,
double (*f)(double t, const double* par),
@@ -81,6 +98,114 @@
{
lmcurve_data_struct data = { t, y, dy, f };
- lmmin(n_par, par, m_dat, (const void*)&data, lmcurve_evaluate,
+ lmmin(n_par, par, m_dat, nullptr, (const void*)&data, lmcurve_evaluate,
control, status);
}
+
+#endif
+
+bool lmfit_exp(int nfit,
+ const double x[],
+ const double y[],
+ const double dy[],
+ double parm[],
+ bool bVerbose,
+ int eFitFn,
+ int nfix)
+{
+ if ((eFitFn < 0) || (eFitFn >= effnNR))
+ {
+ fprintf(stderr, "fitfn = %d, should be in the range 0..%d\n",
+ eFitFn, effnNR-1);
+ return false;
+ }
+#if HAVE_LMFIT
+ double chisq, ochisq;
+ gmx_bool bCont;
+ int j;
+ int maxiter = 100;
+ lm_control_struct control;
+ lm_status_struct *status;
+ int nparam = effnNparams(eFitFn);
+ int p2;
+ gmx_bool bSkipLast;
+
+ /* Using default control structure for double precision fitting that
+ * comes with the lmfit package (i.e. from the include file).
+ */
+ control = lm_control_double;
+ control.verbosity = (bVerbose ? 1 : 0);
+ control.n_maxpri = 0;
+ control.m_maxpri = 0;
+
+ snew(status, 1);
+ /* Initial params */
+ chisq = 1e12;
+ j = 0;
+ if (bVerbose)
+ {
+ printf("%4s %10s Parameters\n", "Step", "chi^2");
+ }
+ /* Check whether we have to skip some params */
+ if (nfix > 0)
+ {
+ do
+ {
+ p2 = 1 << (nparam-1);
+ bSkipLast = ((p2 & nfix) == p2);
+ if (bSkipLast)
+ {
+ nparam--;
+ nfix -= p2;
+ }
+ }
+ while ((nparam > 0) && (bSkipLast));
+ if (bVerbose)
+ {
+ printf("Using %d out of %d parameters\n", nparam, effnNparams(eFitFn));
+ }
+ }
+ do
+ {
+ ochisq = chisq;
+ gmx_lmcurve(nparam, parm, nfit, x, y, dy,
+ lmcurves[eFitFn], &control, status);
+ chisq = gmx::square(status->fnorm);
+ if (bVerbose)
+ {
+ printf("status: fnorm = %g, nfev = %d, userbreak = %d\noutcome = %s\n",
+ status->fnorm, status->nfev, status->userbreak,
+ lm_infmsg[status->outcome]);
+ }
+ if (bVerbose)
+ {
+ int mmm;
+ printf("%4d %8g", j, chisq);
+ for (mmm = 0; (mmm < effnNparams(eFitFn)); mmm++)
+ {
+ printf(" %8g", parm[mmm]);
+ }
+ printf("\n");
+ }
+ j++;
+ bCont = (fabs(ochisq - chisq) > fabs(control.ftol*chisq));
+ }
+ while (bCont && (j < maxiter));
+
+ sfree(status);
+#else
+ gmx_fatal(FARGS, "This build of GROMACS was not configured with support "
+ "for lmfit, so the requested fitting cannot be performed. "
+ "See the install guide for instructions on how to build "
+ "GROMACS with lmfit supported.");
+ GMX_UNUSED_VALUE(nfit);
+ GMX_UNUSED_VALUE(x);
+ GMX_UNUSED_VALUE(y);
+ GMX_UNUSED_VALUE(dy);
+ GMX_UNUSED_VALUE(parm);
+ GMX_UNUSED_VALUE(bVerbose);
+ GMX_UNUSED_VALUE(eFitFn);
+ GMX_UNUSED_VALUE(nfix);
+#endif
+ return true;
+}
diff --git a/src/gromacs/correlationfunctions/gmx_lmcurve.h b/src/gromacs/correlationfunctions/gmx_lmcurve.h
index 0f6c957..fa77671 100644
--- a/src/gromacs/correlationfunctions/gmx_lmcurve.h
+++ b/src/gromacs/correlationfunctions/gmx_lmcurve.h
@@ -1,7 +1,7 @@
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2016, by the GROMACS development team, led by
+ * Copyright (c) 2016,2018, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
@@ -43,13 +43,22 @@
#ifndef GMX_CORRELATION_FUNCTIONS_GMX_LMCURVE_H
#define GMX_CORRELATION_FUNCTIONS_GMX_LMCURVE_H
-#include <lmstruct.h>
+/*! \brief function type for passing to fitting routine */
+typedef double (*t_lmcurve)(double x, const double *a);
-//! Calls lmmin with the given data, with callback function \c f.
-void gmx_lmcurve( const int n_par, double *par, const int m_dat,
- const double *t, const double *y, const double *dy,
- double (*f)(double t, const double *par ),
- const lm_control_struct *control,
- lm_status_struct *status );
+/*! \brief lmfit_exp supports fitting of different functions
+ *
+ * This routine calls the Levenberg-Marquardt non-linear fitting
+ * routine for fitting a data set with errors to a target function.
+ * Fitting routines included in gromacs in src/external/lmfit.
+ */
+bool lmfit_exp(int nfit,
+ const double x[],
+ const double y[],
+ const double dy[],
+ double parm[],
+ bool bVerbose,
+ int eFitFn,
+ int nfix);
#endif
diff --git a/src/gromacs/correlationfunctions/tests/expfit.cpp b/src/gromacs/correlationfunctions/tests/expfit.cpp
index eb93463..5d5bd96 100644
--- a/src/gromacs/correlationfunctions/tests/expfit.cpp
+++ b/src/gromacs/correlationfunctions/tests/expfit.cpp
@@ -1,7 +1,7 @@
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2014,2015,2017, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015,2017,2018, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
@@ -43,6 +43,8 @@
#include "gmxpre.h"
#include "gromacs/correlationfunctions/expfit.h"
+
+#include "config.h"
#include <cmath>
@@ -152,6 +154,10 @@
//static var
std::vector<ExpfitData> ExpfitTest::data_;
+// TODO calling test() leads to a fatal error, which we could in
+// principle test for.
+#if HAVE_LMFIT
+
TEST_F (ExpfitTest, EffnEXP1) {
double param[] = {25};
test(effnEXP1, param, 1e-5, 0);
@@ -202,6 +208,8 @@
test(effnPRES, param, 1e-4, 1);
}
+#endif
+
}
}