1#!/usr/bin/env python
2# Author : Pierre Schnizer
3# $Id$
4
5from __future__ import print_function
6import sys
7import unittest
8import pygsl._numobj as Numeric
9import random
10import pygsl
11from pygsl import Float
12from pygsl import multifit_nlin
13import copy
14
15exp = Numeric.exp
16_eps = 1e-7
17
18def testfunc(t, A = 1., _lambda = .1, b=.5):
19        return  A * exp(- _lambda * t) + b
20
21def exp_f(x, params):
22    A       = x[0]
23    lambda_ = x[1]
24    b       = x[2]
25
26    t       = params[0]
27    yi      = params[1]
28    sigma   = params[2]
29
30    Yi = testfunc(t, A, lambda_, b)
31    f = (Yi -yi) / sigma
32    return f
33
34def exp_df(x, params):
35    A       = x[0]
36    lambda_ = x[1]
37    b       = x[2]
38
39    t       = params[0]
40    yi      = params[1]
41    sigma   = params[2]
42
43    e = exp(-lambda_ * t)
44    e_s = e/sigma
45    df = Numeric.array((e_s, -t * A * e_s, 1/sigma))
46    df = Numeric.transpose(df)
47    return df
48
49def exp_fdf(x, params):
50    f = exp_f(x, params)
51    df = exp_df(x, params)
52    return f, df
53
54
55class DefaultCase(unittest.TestCase):
56    A       = 1.
57    lambda_ = .1
58    b       = .5
59
60    def _getn(self):
61        return 40
62
63    def _getp(self):
64        return 3
65
66
67    def setUp(self):
68        t = Numeric.arange(self._getn());
69        y = testfunc(t, self.A, self.lambda_, self.b)
70        sigma = Numeric.ones(self._getn()) * 0.1
71        self.data = Numeric.array((t,y,sigma), Float)
72        self.sys = multifit_nlin.gsl_multifit_function_fdf(exp_f, exp_df,
73                                                           exp_fdf, self.data,
74                                                           self._getn(),
75                                                           self._getp())
76
77    def _run(self, solver):
78        x = Numeric.array((1.0, 0.0, 0.0))
79        solver.set(x)
80        for iter in range(20):
81            status = solver.iterate()
82            assert(status == 0 or status == -2)
83            x  = solver.getx()
84            dx = solver.getdx()
85            f  = solver.getf()
86            J  = solver.getJ()
87            tdx = multifit_nlin.gradient(J, f)
88            status = multifit_nlin.test_delta(dx, x, 1e-8, 1e-8)
89            fn = Numeric.sqrt(Numeric.sum(f*f))
90            if status == 0:
91                break
92        else:
93            raise ValueError("Number of Iterations exceeded!")
94        assert(Numeric.absolute(x[0] - self.A) < _eps)
95        assert(Numeric.absolute(x[1] - self.lambda_) < _eps)
96        assert(Numeric.absolute(x[2] - self.b) < _eps)
97        covar =  multifit_nlin.covar(solver.getJ(), 0.0)
98
99    def test_lmsder(self):
100        solver = multifit_nlin.lmsder(self.sys, self._getn(), self._getp())
101        self._run(solver)
102
103    def test_lmder(self):
104        solver = multifit_nlin.lmder(self.sys, self._getn(), self._getp())
105        self._run(solver)
106
107    def test_lmniel(self):
108        try:
109             _type = multifit_nlin.lmniel
110        except AttributeError:
111             print("Solver lmniel not available for your GSL version")
112
113        solver = _type(self.sys, self._getn(), self._getp())
114        self._run(solver)
115
116
117if __name__ == '__main__':
118    unittest.main()
119