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
9 //- Class: JacobiOrthogPolynomial
10 //- Description: Class for Jacobi Orthogonal Polynomial
11 //-
12 //- Owner: Mike Eldred, Sandia National Laboratories
13
14 #ifndef JACOBI_ORTHOG_POLYNOMIAL_HPP
15 #define JACOBI_ORTHOG_POLYNOMIAL_HPP
16
17 #include "OrthogonalPolynomial.hpp"
18 #include "BetaRandomVariable.hpp"
19
20
21 namespace Pecos {
22
23 /// Derived orthogonal polynomial class for Jacobi polynomials
24
25 /** The JacobiOrthogPolynomial class evaluates a univariate Jacobi
26 polynomial P^(alpha,beta)_n of a particular order. These
27 polynomials are orthogonal with respect to the weight function
28 (1-x)^alpha (1+x)^beta when integrated over the support range of
29 [-1,+1]. This corresponds to the probability density function
30 f(x) = (1-x)^alpha (1+x)^beta / (2^(alpha+beta+1) B(alpha+1,beta+1))
31 for the beta distribution for [L,U]=[-1,1], where common
32 statistical PDF notation conventions (see, e.g., the uncertain
33 variables section in the DAKOTA Reference Manual) and the
34 Abramowiz and Stegun orthogonal polynomial conventions are
35 inverted and require conversion in this case (alpha_poly =
36 beta_stat - 1; beta_poly = alpha_stat - 1 with the poly
37 definitions used in both cases above). It enables (mixed)
38 multidimensional orthogonal polynomial basis functions within
39 OrthogPolyApproximation. A special case is the
40 LegendreOrthogPolynomial (implemented separately), for which
41 alpha_poly = beta_poly = 0. */
42
43 class JacobiOrthogPolynomial: public OrthogonalPolynomial
44 {
45 public:
46
47 //
48 //- Heading: Constructor and destructor
49 //
50
51 /// default constructor
52 JacobiOrthogPolynomial();
53 /// standard constructor
54 JacobiOrthogPolynomial(Real alpha_stat, Real beta_stat);
55 /// destructor
56 ~JacobiOrthogPolynomial();
57
58 //
59 //- Heading: Virtual function redefinitions
60 //
61
62 /// calculate and return wtFactor based on alphaPoly and betaPoly
63 Real weight_factor();
64
65 protected:
66
67 //
68 //- Heading: Virtual function redefinitions
69 //
70
71 Real type1_value(Real x, unsigned short order);
72 Real type1_gradient(Real x, unsigned short order);
73 Real type1_hessian(Real x, unsigned short order);
74 Real norm_squared(unsigned short order);
75
76 const RealArray& collocation_points(unsigned short order);
77 const RealArray& type1_collocation_weights(unsigned short order);
78
79 void pull_parameter(short dist_param, Real& param) const;
80 void push_parameter(short dist_param, Real param);
81 bool parameterized() const;
82
83 Real length_scale() const;
84
85 private:
86
87 //
88 //- Heading: Data
89 //
90
91 /// the alpha parameter for the Jacobi polynomial as defined by
92 /// Abramowitz and Stegun (differs from statistical PDF notation)
93 Real alphaPoly;
94 /// the beta parameter for the Jacobi polynomial as defined by
95 /// Abramowitz and Stegun (differs from statistical PDF notation)
96 Real betaPoly;
97 };
98
99
JacobiOrthogPolynomial()100 inline JacobiOrthogPolynomial::JacobiOrthogPolynomial():
101 alphaPoly(0.), betaPoly(0.)
102 { collocRule = GAUSS_JACOBI; }
103
104
105 inline JacobiOrthogPolynomial::
JacobiOrthogPolynomial(Real alpha_stat,Real beta_stat)106 JacobiOrthogPolynomial(Real alpha_stat, Real beta_stat):
107 alphaPoly(beta_stat-1.), betaPoly(alpha_stat-1.) // inverted conventions
108 { collocRule = GAUSS_JACOBI; }
109
110
~JacobiOrthogPolynomial()111 inline JacobiOrthogPolynomial::~JacobiOrthogPolynomial()
112 { }
113
114
115 inline void JacobiOrthogPolynomial::
pull_parameter(short dist_param,Real & param) const116 pull_parameter(short dist_param, Real& param) const
117 {
118 // return BE_ALPHA,BETA (unconverted)
119 switch (dist_param) {
120 case JACOBI_ALPHA: param = alphaPoly; break;
121 case JACOBI_BETA: param = betaPoly; break;
122 case BE_ALPHA: param = betaPoly + 1.; break; // beta poly to alpha stat
123 case BE_BETA: param = alphaPoly + 1.; break; // alpha poly to beta stat
124 default:
125 PCerr << "Error: unsupported distribution parameter in JacobiOrthog"
126 << "Polynomial::parameter()." << std::endl;
127 abort_handler(-1); break;
128 }
129 }
130
131
push_parameter(short dist_param,Real param)132 inline void JacobiOrthogPolynomial::push_parameter(short dist_param, Real param)
133 {
134 // *_stat() routines are called for each approximation build from
135 // PolynomialApproximation::update_basis_distribution_parameters().
136 // Logic for first pass included for completeness, but should not be needed.
137 if (collocPointsMap.empty() || collocWeightsMap.empty()) { // first pass
138 switch (dist_param) {
139 case JACOBI_ALPHA: alphaPoly = param; break;
140 case JACOBI_BETA: betaPoly = param; break;
141 case BE_ALPHA: betaPoly = param - 1.; break; // alpha stat to beta poly
142 case BE_BETA: alphaPoly = param - 1.; break; // beta stat to alpha poly
143 }
144 }
145 else
146 switch (dist_param) {
147 case JACOBI_ALPHA:
148 if (!real_compare(alphaPoly, param))
149 { alphaPoly = param; reset_gauss(); }
150 break;
151 case JACOBI_BETA:
152 if (!real_compare(betaPoly, param))
153 { betaPoly = param; reset_gauss(); }
154 break;
155 case BE_ALPHA: {
156 Real bp = param - 1.; // alpha stat to beta poly
157 if (!real_compare(betaPoly, bp))
158 { betaPoly = bp; reset_gauss(); }
159 break;
160 }
161 case BE_BETA: {
162 Real ap = param - 1.; break; // beta stat to alpha poly
163 if (!real_compare(alphaPoly, ap))
164 { alphaPoly = ap; reset_gauss(); }
165 break;
166 }
167 }
168 }
169
170
parameterized() const171 inline bool JacobiOrthogPolynomial::parameterized() const
172 { return true; }
173
174
175 /** return max(mean, stdev) on [-1,1]. */
length_scale() const176 inline Real JacobiOrthogPolynomial::length_scale() const
177 {
178 Real mean, stdev;
179 // BetaRandomVariable uses alpha_stat, beta_stat:
180 BetaRandomVariable::
181 moments_from_params(betaPoly+1., alphaPoly+1., -1., 1., mean, stdev);
182 return std::max(mean, stdev);
183 }
184
185 } // namespace Pecos
186
187 #endif
188