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