gsl_multifit.h

Go to the documentation of this file.

/* $Id: gsl_multifit.h 14889 2013-07-10 20:48:10Z predoehl $ */
/* {{{=========================================================================== *
   |
   |  Copyright (c) 1994-2011 by Kobus Barnard (author)
   |
   |  Personal and educational use of this code is granted, provided that this
   |  header is kept intact, and that the authorship is not misrepresented, that
   |  its use is acknowledged in publications, and relevant papers are cited.
   |
   |  For other use contact the author (kobus AT cs DOT arizona DOT edu).
   |
   |  Please note that the code in this file has not necessarily been adequately
   |  tested. Naturally, there is no guarantee of performance, support, or fitness
   |  for any particular task. Nonetheless, I am interested in hearing about
   |  problems that you encounter.
   |
   |  Author:  Kyle Simek
 * =========================================================================== }}}*/

// vim: tabstop=4 shiftwidth=4 foldmethod=marker

#ifndef KJB_CPP_GSL_MULTIFIT_H
#define KJB_CPP_GSL_MULTIFIT_H

#include <l_cpp/l_exception.h>
#include <gsl_cpp/gsl_util.h>
#include <gsl_cpp/gsl_vector.h>
#include <boost/function/function1.hpp>

#ifdef KJB_HAVE_GSL
#include <gsl/gsl_multifit_nlin.h>
#else
#warning "Gnu Scientific Library is absent, yet essential to this program."

/*
 * Here is some fakery to make client programs compile even when GSL is absent.
 */
typedef void gsl_multifit_fdfsolver_type;
typedef void gsl_multifit_fdfsolver;
struct gsl_matrix {};
#define gsl_multifit_fdfsolver_lmsder ((void*)0)
#define gsl_multifit_fdfsolver_lmder ((void*)0)
#define gsl_multifit_fdfsolver_alloc(x,z) 0
#define gsl_set_error_handler_off() (0)
struct gsl_multifit_function_fdf {
    double (*f)( const gsl_vector*, void* , gsl_vector * );
    void (*df)( const gsl_vector*, void*, gsl_matrix* );
    void (*fdf)( const gsl_vector*, void*, gsl_vector*, gsl_matrix* );
    int n;
    size_t p;
    void *params;
};

typedef void gsl_multifit_fsolver_type;
typedef void gsl_multifit_fsolver;
// #define gsl_multifit_fsolver_nmsimplex2 ((void*)0)
#define gsl_multifit_fsolver_alloc(x,z) 0
//struct gsl_multifit_function {
//    size_t n; void *params;
//    double (*f)( const gsl_vector*, void* );
//};

#endif


namespace kjb {

class Gsl_multifit_fdf {

#ifdef KJB_HAVE_GSL

    struct basic_Gsl_multifit_fdf {
        gsl_multifit_fdfsolver* p;
        size_t dim;
        size_t num_obs;

        basic_Gsl_multifit_fdf( 
            const gsl_multifit_fdfsolver_type* type,
            size_t n,
            size_t p
        )
        :   p( (n>0 && p>0) ? gsl_multifit_fdfsolver_alloc(type, n, p) : 00 ),
            dim( p ),
            num_obs( n )
        {
            ETX_2( 0 == n, "basic_Gsl_multifit_fdf(): zero observations" );
            ETX_2( 0 == p, "basic_Gsl_multifit_fdf(): Zero dimensions" );
            ETX_2( 00 == p, "basic_Gsl_multifit_fdf(): Memory alloc failed" );
        }

        ~basic_Gsl_multifit_fdf()
        {
            gsl_multifit_fdfsolver_free( p );
        }

        void swap( basic_Gsl_multifit_fdf& that )
        {
            using std::swap;
            swap( p, that.p );
        }

    private:
        // this class is not copyable, nor assignable, nor has a default ctor.
        basic_Gsl_multifit_fdf( const basic_Gsl_multifit_fdf& );
        basic_Gsl_multifit_fdf();
        basic_Gsl_multifit_fdf& operator=( const basic_Gsl_multifit_fdf& );
    };

    basic_Gsl_multifit_fdf m_fit;
#endif

    bool m_verbose;

    boost::function1<double, const gsl_vector*> boost_f_;
    bool m_initialized;

    mutable Gsl_Vector* m_gradient;

public:
    Gsl_multifit_fdf(
        const gsl_multifit_fdfsolver_type* type,
        gsl_multifit_function_fdf* f,
        const gsl_vector* x0,
        size_t data_count,
        bool verbosity = true
    )
#ifdef KJB_HAVE_GSL
    :   m_fit( type, data_count, x0 ? x0 -> size : 0 ),
        m_verbose( verbosity ),
        m_initialized(false),
        m_gradient(NULL)
    {
        if(x0->size != f->p)
        {
            KJB_THROW(Dimension_mismatch);
        }

        initialize(f, x0);
    }
#else
    {
        KJB_THROW_2( Missing_dependency, "GNU GSL" );
    }
#endif

    Gsl_multifit_fdf(
            size_t data_count,
            size_t size,
            const gsl_multifit_fdfsolver_type* type = gsl_multifit_fdfsolver_lmsder)
#ifdef KJB_HAVE_GSL
    :   m_fit( type, data_count, size),
        m_verbose( false ),
        m_initialized(false),
        m_gradient(NULL)
    {}
#else
    {
        KJB_THROW_2( Missing_dependency, "GNU GSL" );
    }
#endif

    ~Gsl_multifit_fdf()
    {
        delete m_gradient;
    }

    void initialize(
        gsl_multifit_function_fdf* f,
        const gsl_vector* x0
    )
    {
#ifdef KJB_HAVE_GSL
        GSL_ETX( gsl_multifit_fdfsolver_set( m_fit.p, f, x0 ) );
        m_initialized = true;
#endif
    }

    size_t dim() const
    {
#ifdef KJB_HAVE_GSL
        return m_fit.dim;
#endif
    }

    size_t num_observations() const
    {
#ifdef KJB_HAVE_GSL
        return m_fit.num_obs;
#endif
    }

    int iterate()
    {
#ifdef KJB_HAVE_GSL
        if(!m_initialized)
        {
            KJB_THROW_2(Runtime_error, "Solver not initialized");
        }

        int rc = gsl_multifit_fdfsolver_iterate( m_fit.p );
        if (m_verbose) gsl_iterate_EPE(rc);
        return rc;
#endif
    }

    gsl_vector* position() const
    {
#ifdef KJB_HAVE_GSL
        return gsl_multifit_fdfsolver_position( m_fit.p );
#endif
    }

    gsl_matrix* jacobian() const
    {
#ifdef KJB_HAVE_GSL
        return m_fit.p->J;
#endif
    }


    const gsl_vector* min() const
    {
#ifdef KJB_HAVE_GSL
        return m_fit.p->f;
#endif
    }

    int test_delta( double epsabs, double epsrel) const
    {
#ifdef KJB_HAVE_GSL
        
        return gsl_multifit_test_delta( m_fit.p->dx, m_fit.p->x, epsabs, epsrel );
#endif
    }

    int test_gradient( double epsabs = 1e-6) const
    {
#ifdef KJB_HAVE_GSL
        if(m_gradient == NULL)
            m_gradient = new Gsl_Vector(dim());

        gsl_multifit_gradient(m_fit.p->J, m_fit.p->f, *m_gradient);
        return gsl_multifit_test_gradient(*m_gradient, epsabs);
#endif
    }

    void swap( Gsl_multifit_fdf& that )
    {
#ifdef KJB_HAVE_GSL
        if ( this == &that )
        {
            return;
        }
        m_fit.swap( that.m_fit );

        using std::swap;
        swap( m_verbose, that.m_verbose );
        swap( m_gradient, that.m_gradient );
#endif
    }

    bool verbosity( bool v )
    {
        bool oldv = m_verbose;
        m_verbose = v;
        return oldv;
    }

    bool verbosity( )
    {
        return m_verbose;
    }
};

} // namespace kjb

#endif