1#  _______________________________________________________________________
2#
3#  PECOS: Parallel Environment for Creation Of Stochastics
4#  Copyright (c) 2011, Sandia National Laboratories.
5#  This software is distributed under the GNU Lesser General Public License.
6#  For more information, see the README file in the top Pecos directory.
7#  _______________________________________________________________________
8
9import unittest, numpy, os
10from PyDakota.univariate_polynomials import *
11
12def test_orthogonal_polynomial(samples, poly, true_values, true_derivatives,
13                               true_norms=None):
14    """
15    Test the following functions of a polynomial basis :
16        type1_value
17        type1_gradient
18        norm_squared (if true_norms is not None and poly is an
19                      OrthogonalPolynomial)
20
21    Parameters
22    ----------
23    true_values : np.ndarray (num_terms x num_samples)
24       Exact values of the polynomial for each degree
25
26    true_derivaties : np.ndarray (num_terms x num_samples)
27       Exact values of the 1st derivative of the polynomial for each degree
28
29    true_norms : np.ndarray (num_terms-1)
30        l2_norm of basis from degree 1 and up
31    """
32    assert samples.ndim==1
33    nsamples = samples.shape[0]
34    assert true_values.shape[1]==nsamples
35    nterms = true_values.shape[0]
36    assert true_derivatives.shape[1]==nsamples
37
38    for i in range(nterms):
39      for j in range(nsamples):
40        assert numpy.allclose(poly.type1_value(samples[j],i), true_values[i,j])
41
42    degree = nterms-1
43    assert numpy.allclose(poly.values(samples, degree), true_values.T)
44
45    gradient_nterms=true_derivatives.shape[0]
46    for i in range(gradient_nterms):
47      for j in range(nsamples):
48        assert numpy.allclose(
49          poly.type1_gradient(samples[j],i), true_derivatives[i,j])
50
51    if true_norms is not None:
52        for i in xrange( true_norms.shape[0] ):
53            assert numpy.allclose( true_norms[i], poly.norm_squared( i+1 ) )
54
55
56
57
58def test_orthogonality(poly, degree, eps=2*numpy.finfo( numpy.float ).eps):
59    x = poly.collocation_points(degree+1)
60    w = poly.type1_collocation_weights(degree+1)
61    V = poly.values(x,degree)
62    nterms = V.shape[1]
63    for i in xrange(nterms):
64        for j in xrange(nterms):
65            integral = numpy.sum(V[:,i]*V[:,j]*w)
66            if i!=j:
67                assert numpy.absolute(integral)<eps,(i,j,numpy.absolute(integral))
68            else:
69                l2_norm = poly.norm_squared(i)
70                msg = (j,numpy.absolute(l2_norm-integral))
71                assert numpy.absolute(l2_norm-integral)<eps, msg
72
73class TestUnivariatePolynomials(unittest.TestCase):
74    def setUp( self ):
75        pass
76
77    def test_legendre_polynomial(self):
78        """Test Legendre which are orthogonal to uniform random variables"""
79        poly_type=LEGENDRE_ORTHOG
80        poly = BasisPolynomial(poly_type)
81        x = numpy.linspace(-1., 1., 100)
82        true_values = numpy.array( [0.*x+1., x, ( 3. * x**2 - 1. ) / 2.,
83                                    ( 5. * x**3 - 3.*x ) / 2.,
84                                    ( 35. * x**4 - 30 * x**2 + 3. ) / 8.,
85                                    x * ( 63. * x**4 - 70. * x**2 + 15) / 8.])
86
87        true_derivatives = numpy.array( [0.*x, 0.*x+1., 3.*x] )
88        true_norms = numpy.array([1./3.,1./5.,1./7.,1./9.])
89        test_orthogonal_polynomial(
90            x, poly, true_values, true_derivatives, true_norms )
91
92        degree = 100
93        test_orthogonality(poly, degree, eps=1e-14)
94
95    def test_jacobi_polynomials( self ):
96        """ Jacobi( 0.5, 0.5 ) which are orthogonal to Beta(1.5,1.5)
97        random variables"""
98        poly_type=JACOBI_ORTHOG
99        poly = BasisPolynomial(poly_type)
100        poly.alpha_stat(1.5); poly.beta_stat(1.5)
101        x = numpy.linspace(-1., 1., 100)
102        true_values = numpy.array(
103            [0.*x+1.,
104            3. * x / 2, 5. * ( 4. * x**2 - 1. ) / 8.,
105            35. * x * (2. * x**2 - 1. ) / 16.,
106            63. * ( 16. * x**4 - 12 * x**2 + 1. )/128.,
107            231. * x * ( 16. * x**4 - 16.*x**2+3)/256.])
108        true_derivatives = numpy.array([0.*x,0.*x+3./2.,5.*x])
109
110        # answers found using wolfram alpha
111        #N[int[(JacobiP[5,1/2,1/2,2*x-1])^2*(x)^(1/2)*(1-x)^(1/2)/Gamma[1.5]^2*Gamma[3],{x,0,1}],16]
112        true_norms = numpy.array([0.5625,0.390625, 0.299072265625,
113                                  0.242248535156245, 0.2035560607910156])
114        test_orthogonal_polynomial(
115            x, poly, true_values, true_derivatives, true_norms )
116
117        degree = 98# probem with degree > 98
118        test_orthogonality(poly, degree, eps=1e-14)
119
120if __name__ == '__main__':
121    unittest.main()
122