#!/usr/bin/env python
# Author : Pierre Schnizer 
"""
This module describes routines for computing Chebyshev
approximations to univariate functions.  A Chebyshev approximation is a
truncation of the series f(x) = \sum c_n T_n(x), where the Chebyshev
polynomials T_n(x) = \cos(n \arccos x) provide an orthogonal basis of
polynomials on the interval [-1,1] with the weight function 1 /
\sqrt{1-x^2}.  The first few Chebyshev polynomials are, T_0(x) = 1,
T_1(x) = x, T_2(x) = 2 x^2 - 1.


def f(x, p):
    if x < 0.5:
        return 0.25
    else:
        return 0.75

     
i = 10000;
n = i

cs = cheb_series(40)
F = gsl_function(f, None)
cs.init(F, 0.0, 1.0)

for i in range(100):
    x = i * n / 100.
    r10 = cs.eval_n(10, x)
    r40 = cs.eval(x)
    print "%g %g %g %g" % (x, f(x, None), r10, r40)

"""
import _callback

from  _generic_solver import _workspace
from gsl_function import gsl_function

class cheb_series(_workspace):
    """
    This class manages all internal detail. It provides the space for a
    Chebyshev series of order N.

    """
    _alloc    = _callback.gsl_cheb_alloc
    _free     = _callback.gsl_cheb_free
    _init     = _callback.gsl_cheb_init
    _eval     = _callback.gsl_cheb_eval
    _eval_err = _callback.gsl_cheb_eval_err
    _eval_n   = _callback.gsl_cheb_eval_n
    _eval_n_err = _callback.gsl_cheb_eval_n_err
    #_eval_mode   = _callback.gsl_cheb_eval_mode
    #_eval_mode_e   = _callback.gsl_cheb_eval_mode_e
    _calc_deriv   = _callback.gsl_cheb_calc_deriv
    _calc_integ   = _callback.gsl_cheb_calc_integ
    _get_coeff = _callback.pygsl_cheb_get_coefficients
    _set_coeff = _callback.pygsl_cheb_set_coefficients
    _get_a = _callback.pygsl_cheb_get_a
    _set_a = _callback.pygsl_cheb_set_a
    _get_b = _callback.pygsl_cheb_get_b
    _set_b = _callback.pygsl_cheb_set_b
    _get_f = _callback.pygsl_cheb_get_f
    _set_f = _callback.pygsl_cheb_set_f
    _get_order_sp = _callback.pygsl_cheb_get_order_sp
    _set_order_sp = _callback.pygsl_cheb_set_order_sp

    def __init__(self, size):
        """
        input : n
            n ... number of coefficients
        """
        self._size = size
        _workspace.__init__(self, size)
        
    def init(self, f, a, b):
        """
        This function computes the Chebyshev approximation for the
        function F over the range (a,b) to the previously specified order.
        The computation of the Chebyshev approximation is an O(n^2)
        process, and requires n function evaluations.

        input : f, a, b
            f ... a gsl_function
            a ... lower limit
            b ... upper limit
        """
        return self._init(self._ptr, f.get_ptr(), a, b)

    def eval(self, x):
        """
        This function evaluates the Chebyshev series CS at a given point X

        input : x
             x ... value where the series shall be evaluated.
        """
        return self._eval(self._ptr, x)

    def eval_err(self, x):
        """
         This function computes the Chebyshev series  at a given point X,
         estimating both the series RESULT and its absolute error ABSERR.
         The error estimate is made from the first neglected term in the
         series.
         
         input : x
             x ... value where the error shall be evaluated.
        """
        return self._eval_err(self._ptr, x)

    def eval_n(self, order, x):
        """
         This function evaluates the Chebyshev series CS at a given point
         N, to (at most) the given order ORDER.

        input : n, x
             n ... number of cooefficients
             x ... value where the series shall be evaluated.
         
        """
        return self._eval_n(self._ptr, order, x)

    def eval_n_err(self, order, x):
        """
        This function evaluates a Chebyshev series CS at a given point X,
        estimating both the series RESULT and its absolute error ABSERR,
        to (at most) the given order ORDER.  The error estimate is made
        from the first neglected term in the series.

        input : n, x
             n ... number of cooefficients
             x ... value where the error shall be evaluated.

        """
        return self._eval_n_err(self._ptr, order, x)

#    def eval_mode(self, x, mode):
#        """
#        
#        """
#        return self._eval(self._ptr, x, mode)
#
#    def eval_mode_e(self, x, mode):
#        return self._eval(self._ptr, x, mode)

    def calc_deriv(self):
        """
        This method computes the derivative of the series CS. It returns
        a new instance of the cheb_series class.
        """
        tmp = cheb_series(self._size)
        self._calc_deriv(tmp._ptr, self._ptr)
        return tmp

    def calc_integ(self):
        """
        This method computes the integral of the series CS. It returns
        a new instance of the cheb_series class.

        """
        tmp = cheb_series(self._size)
        self._calc_integ(tmp._ptr, self._ptr)
        return tmp

    def get_coefficients(self):
        """
        Get the chebyshev coefficients. 
        """
        return self._get_coeff(self._ptr)


    def set_coefficients(self, coefs):
        """
        Sets the chebyshev coefficients. 
        """
        return self._set_coeff(self._ptr, coefs)

    def get_a(self):
        """
        Get the lower boundary of the current representation
        """
        return self._get_a(self._ptr)

    def set_a(self, a):
        """
        Set the lower boundary of the current representation
        """
        return self._set_a(self._ptr, a)

    def get_b(self):
        """
        Get the upper boundary of the current representation
        """
        return self._get_b(self._ptr)

    def set_b(self, a):
        """
        Set the upper boundary of the current representation
        """
        return self._set_b(self._ptr, a)

    def get_f(self):
        """
        Get the value f (what is it ?) The documentation does not tell anything
        about it.
        """
        return self._get_f(self._ptr)

    def set_f(self, a):
        """
        Set the value f (what is it ?)
        """
        return self._set_f(self._ptr, a)

    def get_order_sp(self):
        """
        Get the value f (what is it ?) The documentation does not tell anything
        about it.
        """
        return self._get_order_sp(self._ptr)

    def set_order_sp(self, a):
        """
        Set the value f (what is it ?)
        """
        return self._set_order_sp(self._ptr, a)



if __name__ == '__main__':
    def f(x, p):
        if x < 0.5:
            return 0.25
        else:
            return 0.75
        
     

    cs = cheb_series(40)
    F = gsl_function(f, None)
    cs.init(F, 0.0, 1.0)
    i = 10000;
    n = i

    for i in xrange(n):
        x = i / float(n)
        r10 = cs.eval_n(10, x)
        r40 = cs.eval(x)
        print "%g %g %g %g" % (x, f(x, None), r10, r40)