#!/usr/bin/env python
# Author : Pierre Schnizer 
"""
Wrapper over the functions as described in Chapter 33 of the
reference manual.

Routines for approximating a data set with a non linear function

"""

import _callback
from gsl_function import gsl_multifit_function_fdf, gsl_multifit_function
from _generic_solver import _generic_solver

class _fsolver(_generic_solver):
    type = None
    _alloc     = _callback.gsl_multifit_fsolver_alloc  
    _free      = _callback.gsl_multifit_fsolver_free   
    _set       = _callback.gsl_multifit_fsolver_set    
    _name      = _callback.gsl_multifit_fsolver_name   
    _position  = _callback.gsl_multifit_fsolver_position
    _iterate   = _callback.gsl_multifit_fsolver_iterate
    _getx      = _callback.gsl_multifit_fsolver_getx
    _getdx      = _callback.gsl_multifit_fsolver_getdx
    _getf      = _callback.gsl_multifit_fsolver_getf

    def __init__(self, system, n, p):
        self._ptr = None
        assert(self._free != None)
        assert(self._alloc != None)
        assert(self._set != None)
        assert(self._name != None)
        assert(self._iterate != None)
        self._ptr = self._alloc(self.type, n, p)
        self._isset = 0
        self.system = system

    def set(self, x):
        f = self.system.get_ptr()
        self._set(self._ptr, f, x)
        self._isset = 1

    def position(self, x):
        f = self.system.get_ptr()
        self._position(self._ptr, f, x)
        self._isset = 1

    def getx(self):
        return self._getx(self._ptr)

    def getdx(self):
        return self._getdx(self._ptr)

    def getf(self):
        return self._getf(self._ptr)


        
class _fdfsolver(_fsolver):
    type = None
    _alloc     = _callback.gsl_multifit_fdfsolver_alloc  
    _free      = _callback.gsl_multifit_fdfsolver_free   
    _set       = _callback.gsl_multifit_fdfsolver_set    
    _name      = _callback.gsl_multifit_fdfsolver_name   
    _position  = _callback.gsl_multifit_fdfsolver_position
    _iterate   = _callback.gsl_multifit_fdfsolver_iterate
    _getx      = _callback.gsl_multifit_fdfsolver_getx
    _getdx     = _callback.gsl_multifit_fdfsolver_getdx
    _getf      = _callback.gsl_multifit_fdfsolver_getf
    _getJ      = _callback.gsl_multifit_fdfsolver_getJ

    def getJ(self):
        return self._getJ(self._ptr)

class lmder(_fdfsolver):
    """
    This is an unscaled version of the LMDER algorithm.  The elements
    of the diagonal scaling matrix D are set to 1.  This algorithm may
    be useful in circumstances where the scaled version of LMDER
    converges too slowly, or the function is already scaled
    appropriately.    
    """
    type = _callback.cvar.gsl_multifit_fdfsolver_lmder
    
class lmsder(_fdfsolver):
    """
         This is a robust and efficient version of the Levenberg-Marquardt
     algorithm as implemented in the scaled LMDER routine in MINPACK.
     Minpack was written by Jorge J. More', Burton S. Garbow and
     Kenneth E. Hillstrom.

     The algorithm uses a generalized trust region to keep each step
     under control.  In order to be accepted a proposed new position x'
     must satisfy the condition |D (x' - x)| < \delta, where D is a
     diagonal scaling matrix and \delta is the size of the trust
     region.  The components of D are computed internally, using the
     column norms of the Jacobian to estimate the sensitivity of the
     residual to each component of x.  This improves the behavior of the
     algorithm for badly scaled functions.

     On each iteration the algorithm attempts to minimize the linear
     system |F + J p| subject to the constraint |D p| < \Delta.  The
     solution to this constrained linear system is found using the
     Levenberg-Marquardt method.

     The proposed step is now tested by evaluating the function at the
     resulting point, x'.  If the step reduces the norm of the function
     sufficiently, and follows the predicted behavior of the function
     within the trust region. then it is accepted and size of the trust
     region is increased.  If the proposed step fails to improve the
     solution, or differs significantly from the expected behavior
     within the trust region, then the size of the trust region is
     decreased and another trial step is computed.

     The algorithm also monitors the progress of the solution and
     returns an error if the changes in the solution are smaller than
     the machine precision.  The possible error codes are,

    `GSL_ETOLF'
          the decrease in the function falls below machine precision

    `GSL_ETOLX'
          the change in the position vector falls below machine

    `GSL_ETOLG'
          the norm of the gradient, relative to the norm of the
          function, falls below machine precision

     These error codes indicate that further iterations will be
     unlikely to change the solution from its current value.

    """
    type = _callback.cvar.gsl_multifit_fdfsolver_lmsder 

def test_delta(dx, x, epsabs, epsrel):
    """
    This function tests for the convergence of the sequence by
     comparing the last step DX with the absolute error EPSABS and
     relative error EPSREL to the current position X.  The test returns
     `GSL_SUCCESS' if the following condition is achieved,

          |dx_i| < epsabs + epsrel |x_i|

     for each component of X and returns `GSL_CONTINUE' otherwise.

    """
    return _callback.gsl_multifit_test_delta(dx, x, epsabs, epsrel)

def test_gradient(dx, epsabs):
    """
     This function tests the residual gradient G against the absolute
     error bound EPSABS.  Mathematically, the gradient should be
     exactly zero at the minimum. The test returns `GSL_SUCCESS' if the
     following condition is achieved,

          \sum_i |g_i| < epsabs

     and returns `GSL_CONTINUE' otherwise.  This criterion is suitable
     for situations where the precise location of the minimum, x, is
     unimportant provided a value can be found where the gradient is
     small enough.

    """
    return _callback.gsl_multifit_test_gradient(dx, epsabs)

def gradient(J, f):
    """
      This function computes the gradient G of \Phi(x) = (1/2)
     ||F(x)||^2 from the Jacobian matrix J and the function values F,
     using the formula g = J^T f.

    """
    return _callback.gsl_multifit_gradient(J, f)

def covar(J, epsrel):
    """
    This function uses the Jacobian matrix J to compute the covariance
    matrix of the best-fit parameters, COVAR.  The parameter EPSREL is
    used to remove linear-dependent columns when J is rank deficient.
    
    The covariance matrix is given by,
    
    covar = (J^T J)^{-1}
    
    and is computed by QR decomposition of J with column-pivoting.  Any
    columns of R which satisfy
    
          |R_{kk}| <= epsrel |R_{11}|

    are considered linearly-dependent and are excluded from the
    covariance matrix (the corresponding rows and columns of the
    covariance matrix are set to zero).

    input : J, epsrel
       J       ... Jakobian matrix
       epsrel  
    """
    return _callback.gsl_multifit_covar(J, epsrel)