Skip to content

Commit

Permalink
added fitting method "mpfit"
Browse files Browse the repository at this point in the history
this method is another implementation of Lev-Mar, from MPFIT library,
which is based on MINPACK-1. MPFIT is not committed into SVN, but will
be included in tarball.
  • Loading branch information
wojdyr committed Jan 1, 2012
1 parent 7ecec52 commit c6a2854
Show file tree
Hide file tree
Showing 17 changed files with 341 additions and 56 deletions.
8 changes: 7 additions & 1 deletion INSTALL
Expand Up @@ -17,7 +17,9 @@ $ ./configure (optionally followed by options, see: Configuration Options)
$ make
# make install

Uninstallation of the program installed from source is also easy:
Uninstallation
==============

either do "make uninstall" or manually find and remove all files, e.g.:
find /usr/local/ -name \*fityk\* -exec rm -r \{} \;

Expand Down Expand Up @@ -51,6 +53,10 @@ To install the latest version from Git
======================================

Install automake, autoconf, libtool.

Download http://www.physics.wisc.edu/~craigm/idl/down/cmpfit-1.2.tar.gz
unpack it and copy cmpfit-1.2/mpfit.* to fityk/src/cmpfit/

Before the ./configure step, run:
$ autoreconf -iv

Expand Down
4 changes: 3 additions & 1 deletion NEWS
@@ -1,4 +1,6 @@
User-visible changes in version 1.1.2 (2011-1 - ):
User-visible changes in version 1.1.2 (2012-01- ):
* added fitting method "mpfit" -- wrapper around MINPACK-1 based MPFIT library
(http://www.physics.wisc.edu/~craigm/idl/cmpfit.html)

User-visible changes in version 1.1.1 (2011-09-28):
* option exit_on_warning=0/1 was replaced with on_error=stop/exit; new option
Expand Down
24 changes: 3 additions & 21 deletions TODO
@@ -1,6 +1,4 @@

* fitting (Lev-Mar): a better way of handling singular matrix error

* FFT data transformation using FFTW library,
fft[-inv]-xx, where xx is re, im, amp, e.g. @n = fft-re(@3)

Expand Down Expand Up @@ -76,25 +74,9 @@
http://dx.doi.org/10.1002/cem.1005
http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2631518/

* fitting: consider using 3rd-party library for nonlinear fitting/optimization.
The list of libraries in this category:
NLopt: http://ab-initio.mit.edu/wiki/index.php/NLopt
ALGLIB: http://www.alglib.net/ (lsfitcreatewfg)
MINUIT
COOOL: http://coool.mines.edu/
GAUL: http://gaul.sourceforge.net/,
WNLIB: http://www.willnaylor.com/wnlib.html
GSL and 3rd-party extensions to GSL,
netlib: http://www.netlib.org/
OOL: http://OOL.sf.net,
scitbx::lbfgs http://cctbx.sourceforge.net/
http://sal.jyu.fi/B/3/index.shtml,
http://coin-or.org
http://www-unix.mcs.anl.gov/otc/Guide/faq/nonlinear-programming-faq.html
levmar: http://www.ics.forth.gr/~lourakis/levmar/
OpenOpt: http://openopt.org/
http://www.applied-mathematics.net/optimization/rosenbrock.html
MPFIT: http://www.physics.wisc.edu/~craigm/idl/cmpfit.html
* fitting: consider using NLopt for alternative optimization methods
(http://ab-initio.mit.edu/wiki/index.php/NLopt) and a library for
Genetic Algorithms.

* fitting: in addition to the least squares fit, add so-called robust fit
(it is called robust fit in the Numerical Recipes book), i.e. optimize
Expand Down
5 changes: 5 additions & 0 deletions autogen.sh
Expand Up @@ -2,6 +2,11 @@

set -x
(cd doc && make)
# download and unpack cmpfit-1.2.tar.gz
# from http://www.physics.wisc.edu/~craigm/idl/cmpfit.html
echo "copy content of cmpfit-1.2 (or later) to src/cmpfit/"
#rm -rf src/cmpfit/
#cp -r cmpfit-1.2/ src/cmpfit/
autoreconf --install --verbose \
&& ./configure "$@"

15 changes: 12 additions & 3 deletions doc/fit.rst
Expand Up @@ -91,8 +91,8 @@ The fitting method can be set using the set command::

set fitting_method = method

where method is one of: ``levenberg_marquardt``, ``nelder_mead_simplex``,
``genetic_algorithms``.
where method is one of: ``levenberg_marquardt``, ``mpfit``,
``nelder_mead_simplex``, ``genetic_algorithms``.

All non-linear fitting methods are iterative, and there are three common
stopping criteria:
Expand Down Expand Up @@ -220,7 +220,7 @@ Levenberg-Marquardt

This is a standard nonlinear least-squares routine, and involves
computing the first derivatives of functions. For a description
of the L-M method see *Numerical Recipes*, chapter 15.5
of the algorithm see *Numerical Recipes*, chapter 15.5
or Siegmund Brandt, *Data Analysis*, chapter 10.15.
Essentially, it combines an inverse-Hessian method with a steepest
descent method by introducing a |lambda| factor. When |lambda| is equal
Expand Down Expand Up @@ -249,6 +249,15 @@ criteria.

.. |lambda| replace:: *λ*

Since version 1.1.2 it is possible to use another implementations of
the Levenberg-Marquardt method, from MPFIT_ library.
To switch between the two implementation use command::

set fitting_method = mpfit # switch to MPFIT
set fitting_method = levenberg_marquardt # switch back to default one

.. _MPFIT: http://www.physics.wisc.edu/~craigm/idl/cmpfit.html

.. _nelder:

Nelder-Mead Downhill Simplex
Expand Down
7 changes: 2 additions & 5 deletions doc/index.rst
Expand Up @@ -110,11 +110,8 @@ after the end of subscription.

<div class="smallfont">

Students and home users may donate 20% of the normal price to
``wojdyr@gmail.com`` using
`PayPal <https://www.paypal.com/cgi-bin/webscr?cmd=_donations&business=E98FRTPDBQ3L6&lc=US&currency_code=USD&item_name=Fityk>`_,
`MoneyBookers <https://www.moneybookers.com/app/payment.pl?pay_to_email=wojdyr@gmail.com&language=EN&detail1_text=The+amount+can+be+changed+at+the+end+of+the+URL&detail1_description=Fityk&currency=USD&amount=53>`_
or `Flattr <https://flattr.com/donation/give/to/wojdyr>`_.
There is a discount for home users and students,
email wojdyr@gmail.com for details.

.. raw:: html

Expand Down
5 changes: 5 additions & 0 deletions doc/intro.rst
Expand Up @@ -63,3 +63,8 @@ and send me either the modified version or a patch.
The following people have contributed to this manual (in chronological order):
Marcin Wojdyr (maintainer), Stan Gierlotka, Jaap Folmer, Michael Richardson.

One of the fitting methods uses MPFIT_ library (MINPACK-1 Least Squares
Fitting Library in C), which includes software developed by
the University of Chicago, as Operator of Argonne National Laboratory.

.. _MPFIT: http://www.physics.wisc.edu/~craigm/idl/cmpfit.html
3 changes: 2 additions & 1 deletion src/LMfit.cpp
Expand Up @@ -19,7 +19,8 @@ LMfit::LMfit(Ftk* F)

LMfit::~LMfit() {}

// WSSR is also called chi2
// note: WSSR is also called chi2

void LMfit::init()
{
alpha.resize(na_*na_);
Expand Down
13 changes: 6 additions & 7 deletions src/LMfit.h
@@ -1,15 +1,15 @@
// This file is part of fityk program. Copyright (C) Marcin Wojdyr
// Licence: GNU General Public License ver. 2+

#ifndef FITYK__LMFIT__H__
#define FITYK__LMFIT__H__
#include "common.h"
/// Simple implementation of the Levenberg-Marquardt method,
/// uses Jordan elimination with partial pivoting.

#ifndef FITYK_LMFIT_H_
#define FITYK_LMFIT_H_
#include <vector>
#include <map>
#include <string>
#include "common.h"
#include "fit.h"

/// Levenberg-Marquardt method
class LMfit : public Fit
{
public:
Expand All @@ -28,4 +28,3 @@ class LMfit : public Fit
};

#endif

165 changes: 165 additions & 0 deletions src/MPfit.cpp
@@ -0,0 +1,165 @@
// This file is part of fityk program. Copyright (C) Marcin Wojdyr
// Licence: GNU General Public License ver. 2+

#include "MPfit.h"
#include "logic.h"
#include "data.h"

using namespace std;


MPfit::MPfit(Ftk* F)
: Fit(F, "mpfit")
{
// 0 value means default
mp_conf_.ftol = 0.;
mp_conf_.xtol = 0.;
mp_conf_.gtol = 0.;
mp_conf_.epsfcn = 0.;
mp_conf_.stepfactor = 0.;
mp_conf_.covtol = 0.;
mp_conf_.maxiter = 0;
mp_conf_.maxfev = 0;
mp_conf_.nprint = 0;
mp_conf_.douserscale = 0;
mp_conf_.nofinitecheck = 0;
mp_conf_.iterproc = NULL;
}

MPfit::~MPfit()
{
}

void MPfit::init()
{
}

int calculate_for_mpfit(int m, int npar, double *par, double *deviates,
double **derivs, void *mpfit)
{
return static_cast<MPfit*>(mpfit)->
calculate(m, npar, par, deviates, derivs);
}

/*
void on_mpfit_iteration(void *mpfit)
{
static_cast<MPfit*>(mpfit)->on_iteration();
}
*/

int MPfit::on_iteration()
{
return (int) common_termination_criteria(iter_nr_-start_iter_);
}

int MPfit::calculate(int /*m*/, int npar, double *par, double *deviates,
double **derivs)
{
int stop = on_iteration();
if (stop)
return -1; // error code reserved for user function

vector<realt> A(par, par+npar);
//printf("wssr=%g, p0=%g, p1=%g\n", compute_wssr(A, dmdm_), par[0], par[1]);
if (!derivs)
compute_deviates(A, deviates);
else {
++iter_nr_;
compute_derivatives_mp(A, dmdm_, derivs, deviates);
}
return 0;
}

static
const char* mpstatus_to_string(int n)
{
switch (n) {
case MP_ERR_INPUT: return "General input parameter error";
case MP_ERR_NAN: return "User function produced non-finite values";
case MP_ERR_FUNC: return "No user function was supplied";
case MP_ERR_NPOINTS: return "No user data points were supplied";
case MP_ERR_NFREE: return "No free parameters";
case MP_ERR_MEMORY: return "Memory allocation error";
case MP_ERR_INITBOUNDS:
return "Initial values inconsistent w constraints";
case MP_ERR_BOUNDS: return "Initial constraints inconsistent";
case MP_ERR_PARAM: return "General input parameter error";
case MP_ERR_DOF: return "Not enough degrees of freedom";

// Potential success status codes
case MP_OK_CHI: return "Convergence in chi-square value";
case MP_OK_PAR: return "Convergence in parameter value";
case MP_OK_BOTH: return "Convergence in chi2 and parameter value";
case MP_OK_DIR: return "Convergence in orthogonality";
case MP_MAXITER: return "Maximum number of iterations reached";
case MP_FTOL: return "ftol is too small; no further improvement";
case MP_XTOL: return "xtol is too small; no further improvement";
case MP_GTOL: return "gtol is too small; no further improvement";

// user-defined codes
case -1: return "One of user-defined criteria stopped fitting.";
default: return "unexpected status code";
}
}

void MPfit::autoiter()
{
start_iter_ = iter_nr_;
wssr_before_ = compute_wssr(a_orig_, dmdm_);

int m = 0;
v_foreach (DataAndModel*, i, dmdm_)
m += (*i)->data()->get_n();

// this is also handled in Fit::common_termination_criteria()
mp_conf_.maxiter = max_iterations_;
mp_conf_.maxfev = F_->get_settings()->max_wssr_evaluations;

//mp_conf_.ftol = F_->get_settings()->lm_stop_rel_change;

double *perror = new double[na_];
double *a = new double[na_];
mp_par *pars = new mp_par[na_];
for (int i = 0; i < na_; ++i) {
a[i] = a_orig_[i];
mp_par& p = pars[i];
p.fixed = 0;
p.limited[0] = 0; // no lower limit
p.limited[1] = 0; // no upper limit
p.limits[0] = 0.;
p.limits[1] = 0.;
p.parname = NULL;
p.step = 0.; // step size for finite difference
p.relstep = 0.; // relative step size for finite difference
p.side = 3; // Sidedness of finite difference derivative
// 0 - one-sided derivative computed automatically
// 1 - one-sided derivative (f(x+h) - f(x) )/h
// -1 - one-sided derivative (f(x) - f(x-h))/h
// 2 - two-sided derivative (f(x+h) - f(x-h))/(2*h)
// 3 - user-computed analytical derivatives
p.deriv_debug = 0;
p.deriv_reltol = 0.;
p.deriv_abstol = 0.;
}

// zero result_
result_.bestnorm = result_.orignorm = 0.;
result_.niter = result_.nfev = result_.status = 0;
result_.npar = result_.nfree = result_.npegged = result_.nfunc = 0;
result_.resid = result_.covar = NULL;
result_.xerror = perror;

int status = mpfit(calculate_for_mpfit, m, na_, a, pars, &mp_conf_, this,
&result_);
//printf("%d :: %d\n", iter_nr_, result_.niter);
//soft_assert(iter_nr_ + 1 == result_.niter);
//soft_assert(result_.nfev + 1 == evaluations_);
soft_assert(status == result_.status);
F_->msg("mpfit status: " + S(mpstatus_to_string(status)));
post_fit(vector<realt>(a, a+na_), result_.bestnorm);
delete [] pars;
delete [] a;
delete [] perror;
}

36 changes: 36 additions & 0 deletions src/MPfit.h
@@ -0,0 +1,36 @@
// This file is part of fityk program. Copyright (C) Marcin Wojdyr
// Licence: GNU General Public License ver. 2+

/// wrapper around MPFIT (cmpfit) library,
/// http://www.physics.wisc.edu/~craigm/idl/cmpfit.html
/// which is Levenberg-Marquardt implementation based on MINPACK-1

#ifndef FITYK_MPFIT_H_
#define FITYK_MPFIT_H_
#include "fit.h"
#include "cmpfit/mpfit.h"


struct mp_config_struct;

/// Wrapper around CMPFIT
class MPfit : public Fit
{
public:
MPfit(Ftk* F);
~MPfit();
virtual void init(); // called before do_iteration()/autoiter()
void autoiter();

// implementation (must be public to be called inside callback function)
int calculate(int m, int npar, double *par, double *deviates,
double **derivs);
int on_iteration();
private:
mp_config_struct mp_conf_;
mp_result result_;
int start_iter_;
};

#endif

4 changes: 2 additions & 2 deletions src/Makefile.am
Expand Up @@ -30,11 +30,11 @@ libfityk_la_SOURCES = logic.cpp view.cpp lexer.cpp eparser.cpp cparser.cpp \
vm.h transform.h \
settings.h ui.h GAfit.h LMfit.h guess.h NMfit.h \
model.h fit.h voigt.h numfuncs.h \
swig/fityk_lua.cpp swig/luarun.h
swig/fityk_lua.cpp swig/luarun.h \
MPfit.cpp MPfit.h cmpfit/mpfit.c cmpfit/mpfit.h

include_HEADERS = fityk.h


# --- directory cli/ ---
cli_cfityk_SOURCES = cli/gnuplot.cpp cli/main.cpp cli/gnuplot.h
cli_cfityk_LDADD = libfityk.la $(READLINE_LIBS)
Expand Down

0 comments on commit c6a2854

Please sign in to comment.