Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

Merge pull request #19 from kubkon/master

Port nonlinear least-squares fitting (gsl/gsl_multifit_nlin.h).
  • Loading branch information...
commit 16d7914cd5da46739b064251e685feb4ef98c45d 2 parents 9251fcb + 19b2af0
@twiecki authored
View
1  cython_gsl/__init__.pxd
@@ -119,3 +119,4 @@ from cython_gsl.gsl_roots cimport *
from cython_gsl.gsl_min cimport *
from cython_gsl.gsl_fit cimport *
from cython_gsl.gsl_multimin cimport *
+from cython_gsl.gsl_multifit_nlin cimport *
View
113 cython_gsl/gsl_multifit_nlin.pxd
@@ -0,0 +1,113 @@
+from cython_gsl cimport *
+
+cdef extern from "gsl/gsl_multifit_nlin.h":
+
+ int gsl_multifit_gradient (const gsl_matrix * J,
+ const gsl_vector * f,
+ gsl_vector * g) nogil
+
+ int gsl_multifit_covar (const gsl_matrix * J,
+ double epsrel,
+ gsl_matrix * covar) nogil
+
+ # Definition of vector-valued functions with parameters based on
+ # gsl_vector
+ ctypedef struct gsl_multifit_function:
+ int (* f) (const gsl_vector * x, void * params, gsl_vector * f) nogil
+ size_t n
+ size_t p
+ void * params
+
+ ctypedef struct gsl_multifit_fsolver_type:
+ const char * name
+ size_t size
+ int (* alloc) (void * state, size_t n, size_t p) nogil
+ int (* set) (void * state, gsl_multifit_function * function,
+ gsl_vector * x, gsl_vector * f, gsl_vector * dx) nogil
+ int (* iterate) (void * state, gsl_multifit_function * function,
+ gsl_vector * x, gsl_vector * f, gsl_vector * dx) nogil
+ void (* free) (void * state)
+
+ ctypedef struct gsl_multifit_fsolver:
+ const gsl_multifit_fsolver_type * type
+ gsl_multifit_function * function
+ gsl_vector * x
+ gsl_vector * f
+ gsl_vector * dx
+ void * state
+
+ gsl_multifit_fsolver * gsl_multifit_fsolver_alloc (const gsl_multifit_fsolver_type * T,
+ size_t n, size_t p) nogil
+
+ void gsl_multifit_fsolver_free (gsl_multifit_fsolver * s) nogil
+
+ int gsl_multifit_fsolver_set (gsl_multifit_fsolver * s, gsl_multifit_function * f,
+ const gsl_vector * x) nogil
+
+ int gsl_multifit_fsolver_iterate (gsl_multifit_fsolver * s) nogil
+
+ int gsl_multifit_fsolver_driver (gsl_multifit_fsolver * s, const size_t maxiter,
+ const double epsabs, const double epsrel) nogil
+
+ # Definition of vector-valued functions and gradient with parameters
+ # based on gsl_vector
+ ctypedef struct gsl_multifit_function_fdf:
+ int (* f) (const gsl_vector * x, void * params, gsl_vector * f) nogil
+ int (* df) (const gsl_vector * x, void * params, gsl_matrix * df) nogil
+ int (* fdf) (const gsl_vector * x, void * params, gsl_vector * f,
+ gsl_matrix * df) nogil
+ size_t n
+ size_t p
+ void * params
+
+ ctypedef struct gsl_multifit_fdfsolver_type:
+ const char * name
+ size_t size
+ int (* alloc) (void * state, size_t n, size_t p) nogil
+ int (* set) (void * state, gsl_multifit_function_fdf * fdf, gsl_vector * x,
+ gsl_vector * f, gsl_matrix * J, gsl_vector * dx) nogil
+ int (* iterate) (void * state, gsl_multifit_function_fdf * fdf,
+ gsl_vector * x, gsl_vector * f, gsl_matrix * J,
+ gsl_vector * dx) nogil
+ void (* free) (void * state)
+
+ ctypedef struct gsl_multifit_fdfsolver:
+ const gsl_multifit_fdfsolver_type * type
+ gsl_multifit_function_fdf * fdf
+ gsl_vector * x
+ gsl_vector * f
+ gsl_matrix * J
+ gsl_vector * dx
+ void * state
+
+ gsl_multifit_fdfsolver * gsl_multifit_fdfsolver_alloc (const gsl_multifit_fdfsolver_type * T,
+ size_t n, size_t p) nogil
+
+ int gsl_multifit_fdfsolver_set (gsl_multifit_fdfsolver * s,
+ gsl_multifit_function_fdf * fdf,
+ const gsl_vector * x) nogil
+
+ int gsl_multifit_fdfsolver_iterate (gsl_multifit_fdfsolver * s) nogil
+
+ int gsl_multifit_fdfsolver_driver (gsl_multifit_fdfsolver * s, const size_t maxiter,
+ const double epsabs, const double epsrel) nogil
+
+ void gsl_multifit_fdfsolver_free (gsl_multifit_fdfsolver * s) nogil
+
+ const char * gsl_multifit_fdfsolver_name (const gsl_multifit_fdfsolver * s) nogil
+
+ gsl_vector * gsl_multifit_fdfsolver_position (const gsl_multifit_fdfsolver * s) nogil
+
+ int gsl_multifit_test_delta (const gsl_vector * dx, const gsl_vector * x,
+ double epsabs, double epsrel) nogil
+
+ int gsl_multifit_test_gradient (const gsl_vector * g, double epsabs) nogil
+
+ int gsl_multifit_fdfsolver_dif_df (const gsl_vector * x, gsl_multifit_function_fdf * fdf,
+ const gsl_vector * f, gsl_matrix * J) nogil
+
+ int gsl_multifit_fdfsolver_dif_fdf (const gsl_vector * x, gsl_multifit_function_fdf * fdf,
+ gsl_vector * f, gsl_matrix * J) nogil
+
+ const gsl_multifit_fdfsolver_type * gsl_multifit_fdfsolver_lmder
+ const gsl_multifit_fdfsolver_type * gsl_multifit_fdfsolver_lmsder
View
125 cython_gsl/test/multifit_nlin.pyx
@@ -0,0 +1,125 @@
+from cython_gsl cimport *
+
+from libc.math cimport exp, sqrt
+
+ctypedef struct Data:
+ size_t n
+ double * y
+ double * sigma
+
+cdef int expb_f (const gsl_vector * x, void * data, gsl_vector * f) nogil:
+ cdef Data * d = <Data *> data
+ cdef size_t n = d.n
+ cdef double * y = d.y
+ cdef double * sigma = d.sigma
+
+ cdef double A = gsl_vector_get (x, 0)
+ cdef double lambd = gsl_vector_get (x, 1)
+ cdef double b = gsl_vector_get (x, 2)
+
+ cdef size_t i
+ cdef double t, Yi
+
+ for i from 0 <= i < n:
+ t = i
+ Yi = A * exp (-lambd * t) + b
+ gsl_vector_set (f, i, (Yi - y[i])/sigma[i])
+
+ return GSL_SUCCESS
+
+cdef int expb_df (const gsl_vector * x, void * data, gsl_matrix * J) nogil:
+ cdef Data * d = <Data *> data
+ cdef size_t n = d.n
+ cdef double * y = d.y
+ cdef double * sigma = d.sigma
+
+ cdef double A = gsl_vector_get (x, 0)
+ cdef double lambd = gsl_vector_get (x, 1)
+
+ cdef size_t i
+ cdef double t, s, e
+
+ for i from 0 <= i < n:
+ t = i
+ s = sigma[i]
+ e = exp (-lambd * t)
+ gsl_matrix_set (J, i, 0, e/s)
+ gsl_matrix_set (J, i, 1, -t * A * e/s)
+ gsl_matrix_set (J, i, 2, 1/s)
+
+ return GSL_SUCCESS
+
+cdef int expb_fdf (const gsl_vector * x, void * data, gsl_vector * f,
+ gsl_matrix * J) nogil:
+
+ expb_f (x, data, f)
+ expb_df (x, data, J)
+
+ return GSL_SUCCESS
+
+
+def t_gsl_multifit_nlin_example():
+ cdef const gsl_multifit_fdfsolver_type * T
+ cdef gsl_multifit_fdfsolver * s
+ cdef int status = GSL_CONTINUE
+ cdef unsigned int i
+ cdef unsigned int iter = 0
+ cdef size_t n = 40
+ cdef size_t p = 3
+
+ cdef gsl_matrix * covar = gsl_matrix_alloc (p, p)
+ cdef double y[40]
+ cdef double sigma[40]
+ cdef Data d
+ d.n = n
+ d.y = y
+ d.sigma = sigma
+
+ cdef gsl_multifit_function_fdf f
+ cdef double * x_init = [1.0, 0.0, 0.0]
+ cdef gsl_vector_view x = gsl_vector_view_array (x_init, p)
+ cdef const gsl_rng_type * type
+ cdef gsl_rng * r
+
+ gsl_rng_env_setup ()
+
+ type = gsl_rng_default
+ r = gsl_rng_alloc (type)
+
+ f.f = &expb_f
+ f.df = &expb_df
+ f.fdf = &expb_fdf
+ f.n = n
+ f.p = p
+ f.params = &d
+
+ cdef double t
+
+ for i from 0 <= i < n:
+ t = i;
+ y[i] = 1.0 + 5 * exp (-0.1 * t) + gsl_ran_gaussian (r, 0.1)
+ sigma[i] = 0.1
+
+ T = gsl_multifit_fdfsolver_lmsder
+ s = gsl_multifit_fdfsolver_alloc (T, n, p)
+ gsl_multifit_fdfsolver_set (s, &f, &x.vector)
+
+ while (status == GSL_CONTINUE and iter < 500):
+
+ iter += 1
+ status = gsl_multifit_fdfsolver_iterate (s)
+
+ if status:
+ break
+
+ status = gsl_multifit_test_delta (s.dx, s.x, 1e-4, 1e-4)
+
+ A = gsl_vector_get (s.x, 0)
+ lambd = gsl_vector_get (s.x, 1)
+ b = gsl_vector_get (s.x, 2)
+
+ gsl_multifit_fdfsolver_free (s)
+ gsl_matrix_free (covar)
+ gsl_rng_free (r)
+
+ return (A, lambd, b)
View
21 cython_gsl/test/test_multifit_nlin.py
@@ -0,0 +1,21 @@
+import unittest, multifit_nlin
+
+
+class MultifitNlinTest(unittest.TestCase):
+
+ def test_multifit_nlin_example(self):
+ """Tests implementation of nonlinear least-squares fitting
+ against an example provided by GSL docs.
+
+ The docs is available at:
+ http://www.gnu.org/software/gsl/manual/html_node/Example-programs-for-Nonlinear-Least_002dSquares-Fitting.html#Example-programs-for-Nonlinear-Least_002dSquares-Fitting
+ """
+
+ A, lambd, b = multifit_nlin.t_gsl_multifit_nlin_example()
+ self.assertAlmostEqual(A, 5.04536, 5)
+ self.assertAlmostEqual(lambd, 0.10405, 5)
+ self.assertAlmostEqual(b, 1.01925, 5)
+
+
+if __name__ == '__main__':
+ unittest.main()
View
164 examples/multifit_nlin.pyx
@@ -0,0 +1,164 @@
+from cython_gsl cimport *
+
+from libc.math cimport exp, sqrt, pow
+from libc.stdio cimport printf
+
+ctypedef struct Data:
+ size_t n
+ double * y
+ double * sigma
+
+cdef int expb_f (const gsl_vector * x, void * data, gsl_vector * f) nogil:
+ cdef Data * d = <Data *> data
+ cdef size_t n = d.n
+ cdef double * y = d.y
+ cdef double * sigma = d.sigma
+
+ cdef double A = gsl_vector_get (x, 0)
+ cdef double lambd = gsl_vector_get (x, 1)
+ cdef double b = gsl_vector_get (x, 2)
+
+ cdef size_t i
+ cdef double t, Yi
+
+ for i from 0 <= i < n:
+ # Model Yi = A * exp(-lambd * i) + b
+ t = i
+ Yi = A * exp (-lambd * t) + b
+ gsl_vector_set (f, i, (Yi - y[i])/sigma[i])
+
+ return GSL_SUCCESS
+
+cdef int expb_df (const gsl_vector * x, void * data, gsl_matrix * J) nogil:
+ cdef Data * d = <Data *> data
+ cdef size_t n = d.n
+ cdef double * y = d.y
+ cdef double * sigma = d.sigma
+
+ cdef double A = gsl_vector_get (x, 0)
+ cdef double lambd = gsl_vector_get (x, 1)
+
+ cdef size_t i
+ cdef double t, s, e
+
+ for i from 0 <= i < n:
+ # Jacobian matrix J(i,j) = dfi / dxj,
+ # where fi = (Yi - yi)/sigma[i],
+ # Yi = A * exp(-lambda * i) + b
+ # and the xj are the parameters (A, lambda, b)
+ t = i
+ s = sigma[i]
+ e = exp (-lambd * t)
+ gsl_matrix_set (J, i, 0, e/s)
+ gsl_matrix_set (J, i, 1, -t * A * e/s)
+ gsl_matrix_set (J, i, 2, 1/s)
+
+ return GSL_SUCCESS
+
+cdef int expb_fdf (const gsl_vector * x, void * data, gsl_vector * f,
+ gsl_matrix * J) nogil:
+
+ expb_f (x, data, f)
+ expb_df (x, data, J)
+
+ return GSL_SUCCESS
+
+cdef void print_state (size_t iter, gsl_multifit_fdfsolver * s) nogil:
+ printf ("iter: %3u x = % 15.8f % 15.8f % 15.8f |f(x)| = %g\n",
+ iter,
+ gsl_vector_get (s.x, 0),
+ gsl_vector_get (s.x, 1),
+ gsl_vector_get (s.x, 2),
+ gsl_blas_dnrm2 (s.f))
+
+def main():
+ cdef const gsl_multifit_fdfsolver_type * T
+ cdef gsl_multifit_fdfsolver * s
+ cdef int status = GSL_CONTINUE
+ cdef unsigned int i
+ cdef unsigned int iter = 0
+ cdef size_t n = 40
+ cdef size_t p = 3
+
+ cdef gsl_matrix * covar = gsl_matrix_alloc (p, p)
+ cdef double y[40]
+ cdef double sigma[40]
+ cdef Data d
+ d.n = n
+ d.y = y
+ d.sigma = sigma
+
+ cdef gsl_multifit_function_fdf f
+ cdef double * x_init = [1.0, 0.0, 0.0]
+ cdef gsl_vector_view x = gsl_vector_view_array (x_init, p)
+ cdef const gsl_rng_type * type
+ cdef gsl_rng * r
+
+ gsl_rng_env_setup ()
+
+ type = gsl_rng_default
+ r = gsl_rng_alloc (type)
+
+ f.f = &expb_f
+ f.df = &expb_df
+ f.fdf = &expb_fdf
+ f.n = n
+ f.p = p
+ f.params = &d
+
+ # This is the data to be fitted
+ cdef double t
+
+ for i from 0 <= i < n:
+ t = i;
+ y[i] = 1.0 + 5 * exp (-0.1 * t) + gsl_ran_gaussian (r, 0.1)
+ sigma[i] = 0.1
+ printf ("data: %u %g %g\n", i, y[i], sigma[i])
+
+ T = gsl_multifit_fdfsolver_lmsder
+ s = gsl_multifit_fdfsolver_alloc (T, n, p)
+ gsl_multifit_fdfsolver_set (s, &f, &x.vector)
+
+ print_state (iter, s)
+
+ while (status == GSL_CONTINUE and iter < 500):
+
+ iter += 1
+ status = gsl_multifit_fdfsolver_iterate (s)
+
+ printf ("status = %d\n", status)
+
+ print_state (iter, s)
+
+ if status:
+ break
+
+ status = gsl_multifit_test_delta (s.dx, s.x, 1e-4, 1e-4)
+
+ gsl_multifit_covar (s.J, 0.0, covar)
+
+ cdef double chi = gsl_blas_dnrm2 (s.f)
+ cdef double dof = n - p
+ cdef double c = GSL_MAX_DBL(1, chi / sqrt (dof))
+
+ printf("chisq/dof = %g\n", pow (chi, 2.0) / dof)
+
+ cdef double fit, err
+
+ fit = gsl_vector_get (s.x, 0)
+ err = sqrt (gsl_matrix_get (covar, 0, 0))
+ printf ("A = %.5f +/- %.5f\n", fit, c*err)
+
+ fit = gsl_vector_get (s.x, 1)
+ err = sqrt (gsl_matrix_get (covar, 1, 1))
+ printf ("lambda = %.5f +/- %.5f\n", fit, c*err)
+
+ fit = gsl_vector_get (s.x, 2)
+ err = sqrt (gsl_matrix_get (covar, 2, 2))
+ printf ("b = %.5f +/- %.5f\n", fit, c*err)
+
+ printf ("status = %d\n", status)
+
+ gsl_multifit_fdfsolver_free (s)
+ gsl_matrix_free (covar)
+ gsl_rng_free (r)
Please sign in to comment.
Something went wrong with that request. Please try again.