# Author : Pierre Schnizer 
import sys
from pygsl import  roots
import unittest
#import unittestgui
import pygsl._numobj as Numeric
import pygsl
sys.stdout = sys.stderr

_eps = 1e-3
def quadratic(x, params):
    #sys.stderr.write(str(params))
    a = params[0]
    b = params[1]
    c = params[2]
    tmp =  a * x ** 2 + b * x + c 
    #sys.stderr.write("\ta %f b %f c %f x %30.25f  y %30.25f\n" % (a, b, c, x, tmp))
    return tmp

def quadratic_deriv(x, params):
    a = params[0]
    b = params[1]
    c = params[2]
    dy= 2.0 * a * x + b
    #sys.stderr.write("\ta %f b %f c %f x %30.25f dy %30.25f\n" % (a, b, c, x, dy))
    return dy

def quadratic_fdf(x, params):
    y = quadratic(x, params)
    dy = quadratic_deriv(x, params)
    return y , dy

class DefaultTest(unittest.TestCase):

    def _run(self, solver):
        solver.set(0.0, 5.0)
        iter = 0
        r_expected = Numeric.sqrt(5.0)
        #print "Testing solver ", solver.name() 
        #print "%5s [%9s %9s]  %9s  %10s %9s" % ("iter", "upper", "lower", "root", "err", "err(est)")
        for iter in range(100):           
            status = solver.iterate()
            x_lo = solver.x_lower()
            x_up = solver.x_upper()
            status = roots.test_interval(x_lo, x_up, 0, 0.001)
            r = solver.root()
            if status == 0:
                break
        else:
            raise ValueError, "Exeeded maximum number of iterations!"
        assert(Numeric.absolute(r-r_expected) < _eps)
            #print "%5d [%.7f %.7f]  %.7f  % .6f % .6f" %(iter, x_lo, x_up, r, r -r_expected, x_up - x_lo)
        #print "Convereged :"    
        #print "%5d [%.7f %.7f]  %.7f  % .6f % .6f" %(iter, x_lo, x_up, r, r -r_expected, x_up - x_lo)

class FTest(DefaultTest):        
    def setUp(self):
        a = 1.0
        b = 0.0
        c = -5.0
        quadratic(0, (a,b,c))
        self.sys = roots.gsl_function(quadratic, (a,b,c))

    def test_brent(self):
        solver = roots.brent(self.sys)
        self._run(solver)

    def test_bisection(self):
        solver = roots.bisection(self.sys)
        self._run(solver)

    def test_falsepos(self):
        solver = roots.falsepos(self.sys)
        self._run(solver)


class FdFTest(DefaultTest):        
    def setUp(self):
        a = 1.0
        b = 0.0
        c = -5.0
        self.sys = roots.gsl_function_fdf(quadratic, quadratic_deriv, quadratic_fdf, Numeric.array((a,b,c)))

    def _run(self, solver):
        x = 5.0
        solver.set(x)
        r_expected = Numeric.sqrt(5.0)
        #print "Testing solver ", solver.name() 
        #print "%5s %9s  %10s %9s" % ("iter", "root", "err", "err(est)")
        ok = 1
        for iter in range(10):
            status = solver.iterate()
            x0 = x
            x = solver.root()
            status = roots.test_delta(x, x0, 0.0, 1e-3)
            r = solver.root()
            if status == 0:
                #print "Convereged :"
                ok = 0
            #print "%5d  %.7f  % .6f % .6f" %(iter, r, r -r_expected, x - x0)
            if ok == 0:
                break
        else:
            raise ValueError, "Exeeded maximum number of iterations!"
        assert(Numeric.absolute(r-r_expected) < _eps)
        
    def test_newton(self):
        solver = roots.newton(self.sys)
        self._run(solver)

    def test_secant(self):
        solver = roots.secant(self.sys)
        self._run(solver)

    def test_steffenson(self):
        solver = roots.steffenson(self.sys)
        self._run(solver)

def fn1(x, params):
    return Numeric.cos(x) + 1.0


def suite():
    suite = unittest.TestSuite()
    suite.addTest(FdFTest("test_secant"))
    suite.addTest(FdFTest("test_steffenson"))
    suite.addTest(FdFTest("test_newton"))
    suite.addTest(FTest("test_bisection"))
    suite.addTest(FTest("test_brent"))
    suite.addTest(FTest("test_falsepos"))
    return suite

def run():
    runner = unittest.TextTestRunner()
    runner.run(suite())
    
if __name__ == '__main__':
    run()