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: BasisPolynomial
10 //- Description: Class implementation of base class for basis polynomials
11 //-
12 //- Owner: Mike Eldred, Sandia National Laboratories
13
14 #include "BasisPolynomial.hpp"
15 #include "HermiteOrthogPolynomial.hpp"
16 #include "LegendreOrthogPolynomial.hpp"
17 #include "JacobiOrthogPolynomial.hpp"
18 #include "LaguerreOrthogPolynomial.hpp"
19 #include "GenLaguerreOrthogPolynomial.hpp"
20 #include "ChebyshevOrthogPolynomial.hpp"
21 #include "NumericGenOrthogPolynomial.hpp"
22 #include "LagrangeInterpPolynomial.hpp"
23 #include "HermiteInterpPolynomial.hpp"
24 #include "PiecewiseInterpPolynomial.hpp"
25 #include "KrawtchoukOrthogPolynomial.hpp"
26 #include "MeixnerOrthogPolynomial.hpp"
27 #include "CharlierOrthogPolynomial.hpp"
28 #include "HahnOrthogPolynomial.hpp"
29
30
31 namespace Pecos {
32
33 /** This constructor is the one which must build the base class data
34 for all derived classes. get_polynomial() instantiates a derived
35 class letter and the derived constructor selects this base class
36 constructor in its initialization list (to avoid recursion in the
37 base class constructor calling get_polynomial() again). Since the
38 letter IS the representation, its rep pointer is set to NULL. */
BasisPolynomial(BaseConstructor)39 BasisPolynomial::BasisPolynomial(BaseConstructor):// basisPolyType(-1),
40 wtFactor(1.), ptFactor(1.)
41 { /* empty ctor */ }
42
43
44 /** The default constructor is used in Array<BasisPolynomial>
45 instantiations and by the alternate envelope constructor. polyRep
46 is NULL in this case (problem_db is needed to build a meaningful
47 instance). */
BasisPolynomial()48 BasisPolynomial::BasisPolynomial()
49 { /* empty ctor */ }
50
51
52 /** Envelope constructor which does not require access to problem_db.
53 This constructor executes get_polynomial(type), which invokes the
54 default constructor of the derived letter class, which in turn
55 invokes the BaseConstructor of the base class. */
BasisPolynomial(short poly_type,short rule)56 BasisPolynomial::BasisPolynomial(short poly_type, short rule):
57 // Set the rep pointer to the appropriate derived type
58 polyRep(get_polynomial(poly_type, rule))
59 {
60 if (poly_type && !polyRep) // bad type, insufficient memory
61 abort_handler(-1);
62 }
63
64
65 /** Used only by the envelope constructor to initialize polyRep to the
66 appropriate derived type. */
67 std::shared_ptr<BasisPolynomial>
get_polynomial(short poly_type,short rule)68 BasisPolynomial::get_polynomial(short poly_type, short rule)
69 {
70 std::shared_ptr<BasisPolynomial> polynomial;
71 // In orthogonal polynomial and global interpolation polynomial cases,
72 // basisPolyType is not available at construct time, but is thereafter.
73 switch (poly_type) {
74 case NO_POLY:
75 break;
76 case HERMITE_ORTHOG: // var_type == "normal"
77 polynomial = (rule) ? std::make_shared<HermiteOrthogPolynomial>(rule) :
78 std::make_shared<HermiteOrthogPolynomial>();
79 if (polynomial) polynomial->basisPolyType = poly_type; break;
80 case LEGENDRE_ORTHOG: // var_type == "uniform"
81 polynomial = (rule) ? std::make_shared<LegendreOrthogPolynomial>(rule) :
82 std::make_shared<LegendreOrthogPolynomial>();
83 if (polynomial) polynomial->basisPolyType = poly_type; break;
84 case LAGUERRE_ORTHOG: // var_type == "exponential"
85 polynomial = std::make_shared<LaguerreOrthogPolynomial>();
86 if (polynomial) polynomial->basisPolyType = poly_type; break;
87 case JACOBI_ORTHOG: // var_type == "beta"
88 polynomial = std::make_shared<JacobiOrthogPolynomial>();
89 if (polynomial) polynomial->basisPolyType = poly_type; break;
90 case GEN_LAGUERRE_ORTHOG: // var_type == "gamma"
91 polynomial = std::make_shared<GenLaguerreOrthogPolynomial>();
92 if (polynomial) polynomial->basisPolyType = poly_type; break;
93 case CHEBYSHEV_ORTHOG: // for Clenshaw-Curtis and Fejer
94 polynomial = (rule) ? std::make_shared<ChebyshevOrthogPolynomial>(rule) :
95 std::make_shared<ChebyshevOrthogPolynomial>();
96 if (polynomial) polynomial->basisPolyType = poly_type; break;
97 case NUM_GEN_ORTHOG:
98 polynomial = std::make_shared<NumericGenOrthogPolynomial>();
99 if (polynomial) polynomial->basisPolyType = poly_type; break;
100 case LAGRANGE_INTERP:
101 polynomial = std::make_shared<LagrangeInterpPolynomial>();
102 if (polynomial) polynomial->basisPolyType = poly_type; break;
103 case HERMITE_INTERP:
104 polynomial = (rule) ? std::make_shared<HermiteInterpPolynomial>(rule) :
105 std::make_shared<HermiteInterpPolynomial>();
106 if (polynomial) polynomial->basisPolyType = poly_type; break;
107 // PIECEWISE options include poly order, point type, and point data order:
108 // LINEAR/QUADRATIC/CUBIC covers poly order, rule covers EQUIDISTANT/GENERAL
109 // point type, and data order is inferred from poly order (grads for CUBIC).
110 case PIECEWISE_LINEAR_INTERP: case PIECEWISE_QUADRATIC_INTERP:
111 case PIECEWISE_CUBIC_INTERP:
112 polynomial = (rule) ?
113 std::make_shared<PiecewiseInterpPolynomial>(poly_type, rule) :
114 std::make_shared<PiecewiseInterpPolynomial>(poly_type);
115 break;
116 // Some Discrete orthogonal polynomials
117 case KRAWTCHOUK_DISCRETE: // var_type == "binomial"
118 polynomial = std::make_shared<KrawtchoukOrthogPolynomial>();
119 if (polynomial) polynomial->basisPolyType = poly_type; break;
120 case MEIXNER_DISCRETE: // var_type == "negative binomial"
121 polynomial = std::make_shared<MeixnerOrthogPolynomial>();
122 if (polynomial) polynomial->basisPolyType = poly_type; break;
123 case CHARLIER_DISCRETE: // var_type == "poisson"
124 polynomial = std::make_shared<CharlierOrthogPolynomial>();
125 if (polynomial) polynomial->basisPolyType = poly_type; break;
126 case HAHN_DISCRETE: // var_type == "hygergeometric"
127 polynomial = std::make_shared<HahnOrthogPolynomial>();
128 if (polynomial) polynomial->basisPolyType = poly_type; break;
129 default:
130 PCerr << "Error: BasisPolynomial type " << poly_type << " not available."
131 << std::endl;
132 break;
133 }
134 return polynomial;
135 }
136
137
138 /** Copy constructor manages sharing of polyRep. */
BasisPolynomial(const BasisPolynomial & polynomial)139 BasisPolynomial::BasisPolynomial(const BasisPolynomial& polynomial):
140 polyRep(polynomial.polyRep)
141 { /* empty ctor */ }
142
143
operator =(const BasisPolynomial & polynomial)144 BasisPolynomial BasisPolynomial::operator=(const BasisPolynomial& polynomial)
145 {
146 polyRep = polynomial.polyRep;
147 return *this; // calls copy constructor since returned by value
148 }
149
150
~BasisPolynomial()151 BasisPolynomial::~BasisPolynomial()
152 { /* empty dtor */ }
153
154
type1_value(unsigned short n)155 Real BasisPolynomial::type1_value(unsigned short n)
156 {
157 if (!polyRep) {
158 PCerr << "Error: type1_value(unsigned short) not available for this basis "
159 << "polynomial type." << std::endl;
160 abort_handler(-1);
161 }
162 return polyRep->type1_value(n);
163 }
164
165
type1_value(Real x,unsigned short n)166 Real BasisPolynomial::type1_value(Real x, unsigned short n)
167 {
168 if (!polyRep) {
169 PCerr << "Error: type1_value(Real, unsigned short) not available for this "
170 << "basis polynomial type." << std::endl;
171 abort_handler(-1);
172 }
173 return polyRep->type1_value(x, n);
174 }
175
176
type2_value(Real x,unsigned short n)177 Real BasisPolynomial::type2_value(Real x, unsigned short n)
178 {
179 if (!polyRep) {
180 PCerr << "Error: type2_value(Real, unsigned short) not available for this "
181 << "basis polynomial type." << std::endl;
182 abort_handler(-1);
183 }
184 return polyRep->type2_value(x, n);
185 }
186
187
type1_gradient(unsigned short n)188 Real BasisPolynomial::type1_gradient(unsigned short n)
189 {
190 if (!polyRep) {
191 PCerr << "Error: type1_gradient(unsigned short) not available for this "
192 << "basis polynomial type." << std::endl;
193 abort_handler(-1);
194 }
195 return polyRep->type1_gradient(n);
196 }
197
198
type1_gradient(Real x,unsigned short n)199 Real BasisPolynomial::type1_gradient(Real x, unsigned short n)
200 {
201 if (!polyRep) {
202 PCerr << "Error: type1_gradient(Real, unsigned short) not available for "
203 << "this basis polynomial type." << std::endl;
204 abort_handler(-1);
205 }
206 return polyRep->type1_gradient(x, n);
207 }
208
209
type2_gradient(Real x,unsigned short n)210 Real BasisPolynomial::type2_gradient(Real x, unsigned short n)
211 {
212 if (!polyRep) {
213 PCerr << "Error: type2_gradient(Real, unsigned short) not available for "
214 << "this basis polynomial type." << std::endl;
215 abort_handler(-1);
216 }
217 return polyRep->type2_gradient(x, n);
218 }
219
220
type1_hessian(Real x,unsigned short n)221 Real BasisPolynomial::type1_hessian(Real x, unsigned short n)
222 {
223 if (!polyRep) {
224 PCerr << "Error: type1_hessian(Real, unsigned short) not available for "
225 << "this basis polynomial type." << std::endl;
226 abort_handler(-1);
227 }
228 return polyRep->type1_hessian(x, n);
229 }
230
231
norm_squared(unsigned short n)232 Real BasisPolynomial::norm_squared(unsigned short n)
233 {
234 if (!polyRep) {
235 PCerr << "Error: norm_squared(unsigned short) not available for this basis "
236 << "polynomial type." << std::endl;
237 abort_handler(-1);
238 }
239 return polyRep->norm_squared(n);
240 }
241
242
collocation_points(unsigned short n)243 const RealArray& BasisPolynomial::collocation_points(unsigned short n)
244 {
245 if (!polyRep) {
246 PCerr << "Error: collocation_points() not available for this basis "
247 << "polynomial type." << std::endl;
248 abort_handler(-1);
249 }
250 return polyRep->collocation_points(n);
251 }
252
253
type1_collocation_weights(unsigned short n)254 const RealArray& BasisPolynomial::type1_collocation_weights(unsigned short n)
255 {
256 if (!polyRep) {
257 PCerr << "Error: type1_collocation_weights() not available for this basis "
258 << "polynomial type." << std::endl;
259 abort_handler(-1);
260 }
261 return polyRep->type1_collocation_weights(n);
262 }
263
264
type2_collocation_weights(unsigned short n)265 const RealArray& BasisPolynomial::type2_collocation_weights(unsigned short n)
266 {
267 if (!polyRep) {
268 PCerr << "Error: type2_collocation_weights() not available for this basis "
269 << "polynomial type." << std::endl;
270 abort_handler(-1);
271 }
272 return polyRep->type2_collocation_weights(n);
273 }
274
275
set_new_point(Real x,short order)276 void BasisPolynomial::set_new_point(Real x, short order)
277 {
278 if (polyRep)
279 polyRep->set_new_point(x, order);
280 else {
281 PCerr << "Error: set_new_point(Real, short) not available for this basis "
282 << "polynomial type." << std::endl;
283 abort_handler(-1);
284 }
285 }
286
287
288 void BasisPolynomial::
set_new_point(Real x,short order,const UShortArray & delta_key)289 set_new_point(Real x, short order, const UShortArray& delta_key)
290 {
291 if (polyRep)
292 polyRep->set_new_point(x, order, delta_key);
293 else {
294 PCerr << "Error: set_new_point(Real, short, UShortArray) not available for "
295 << "this basis polynomial type." << std::endl;
296 abort_handler(-1);
297 }
298 }
299
300
exact_index() const301 size_t BasisPolynomial::exact_index() const
302 {
303 if (!polyRep) {
304 PCerr << "Error: exact_index() not available for this basis polynomial "
305 << "type." << std::endl;
306 abort_handler(-1);
307 }
308 return polyRep->exact_index();
309 }
310
311
exact_delta_index() const312 size_t BasisPolynomial::exact_delta_index() const
313 {
314 if (!polyRep) {
315 PCerr << "Error: exact_delta_index() not available for this basis "
316 << "polynomial type." << std::endl;
317 abort_handler(-1);
318 }
319 return polyRep->exact_delta_index();
320 }
321
322
barycentric_value_factors() const323 const RealVector& BasisPolynomial::barycentric_value_factors() const
324 {
325 if (!polyRep) {
326 PCerr << "Error: barycentric_value_factors() not available for this basis "
327 << "polynomial type." << std::endl;
328 abort_handler(-1);
329 }
330 return polyRep->barycentric_value_factors();
331 }
332
333
barycentric_value_factor(unsigned short i) const334 Real BasisPolynomial::barycentric_value_factor(unsigned short i) const
335 {
336 if (!polyRep) {
337 PCerr << "Error: barycentric_value_factor() not available for this basis "
338 << "polynomial type." << std::endl;
339 abort_handler(-1);
340 }
341 return polyRep->barycentric_value_factor(i);
342 }
343
344
barycentric_gradient_factors() const345 const RealVector& BasisPolynomial::barycentric_gradient_factors() const
346 {
347 if (!polyRep) {
348 PCerr << "Error: barycentric_gradient_factors() not available for this "
349 << "basis polynomial type." << std::endl;
350 abort_handler(-1);
351 }
352 return polyRep->barycentric_gradient_factors();
353 }
354
355
barycentric_gradient_factor(unsigned short i) const356 Real BasisPolynomial::barycentric_gradient_factor(unsigned short i) const
357 {
358 if (!polyRep) {
359 PCerr << "Error: barycentric_gradient_factor() not available for this "
360 << "basis polynomial type." << std::endl;
361 abort_handler(-1);
362 }
363 return polyRep->barycentric_gradient_factor(i);
364 }
365
366
barycentric_value_factor_sum() const367 Real BasisPolynomial::barycentric_value_factor_sum() const
368 {
369 if (!polyRep) {
370 PCerr << "Error: barycentric_value_factor_sum() not available for this "
371 << "basis polynomial type." << std::endl;
372 abort_handler(-1);
373 }
374 return polyRep->barycentric_value_factor_sum();
375 }
376
377
barycentric_difference_product() const378 Real BasisPolynomial::barycentric_difference_product() const
379 {
380 if (!polyRep) {
381 PCerr << "Error: barycentric_difference_product() not available for this "
382 << "basis polynomial type." << std::endl;
383 abort_handler(-1);
384 }
385 return polyRep->barycentric_difference_product();
386 }
387
388
reset_gauss()389 void BasisPolynomial::reset_gauss()
390 {
391 if (polyRep)
392 polyRep->reset_gauss();
393 else {
394 PCerr << "Error: reset_gauss() not available for this basis polynomial "
395 << "type." << std::endl;
396 abort_handler(-1);
397 }
398 }
399
400
parameter_update() const401 bool BasisPolynomial::parameter_update() const
402 {
403 if (polyRep)
404 return polyRep->parameter_update();
405 else // default for non-parametric polynomials (e.g., PiecewiseInterp)
406 return false;
407 }
408
409
points_defined(unsigned short order) const410 bool BasisPolynomial::points_defined(unsigned short order) const
411 {
412 if (polyRep)
413 return polyRep->points_defined(order);
414 else // default
415 return false;
416 }
417
418
type1_weights_defined(unsigned short order) const419 bool BasisPolynomial::type1_weights_defined(unsigned short order) const
420 {
421 if (polyRep)
422 return polyRep->type1_weights_defined(order);
423 else // default
424 return false;
425 }
426
427
type2_weights_defined(unsigned short order) const428 bool BasisPolynomial::type2_weights_defined(unsigned short order) const
429 {
430 if (polyRep)
431 return polyRep->type2_weights_defined(order);
432 else // default
433 return false;
434 }
435
436
point_factor()437 Real BasisPolynomial::point_factor()
438 {
439 if (polyRep)
440 return polyRep->point_factor();
441 else // default is used whenever ptFactor does not need to be updated
442 return ptFactor;
443 }
444
445
weight_factor()446 Real BasisPolynomial::weight_factor()
447 {
448 if (polyRep)
449 return polyRep->weight_factor();
450 else // default is used whenever wtFactor does not need to be updated
451 return wtFactor;
452 }
453
454
pull_parameter(short dist_param,Real & param) const455 void BasisPolynomial::pull_parameter(short dist_param, Real& param) const
456 {
457 if (polyRep)
458 polyRep->pull_parameter(dist_param, param);
459 else {
460 PCerr << "Error: pull_parameter(Real) not available for this basis "
461 << "polynomial type." << std::endl;
462 abort_handler(-1);
463 }
464 }
465
466
467 void BasisPolynomial::
pull_parameter(short dist_param,unsigned int & param) const468 pull_parameter(short dist_param, unsigned int& param) const
469 {
470 if (polyRep)
471 polyRep->pull_parameter(dist_param, param);
472 else {
473 PCerr << "Error: pull_parameter(unsigned int) not available for this basis "
474 << "polynomial type." << std::endl;
475 abort_handler(-1);
476 }
477 }
478
479
push_parameter(short dist_param,Real param)480 void BasisPolynomial::push_parameter(short dist_param, Real param)
481 {
482 if (polyRep)
483 polyRep->push_parameter(dist_param, param);
484 else {
485 PCerr << "Error: push_parameter(Real) not available for this basis "
486 << "polynomial type." << std::endl;
487 abort_handler(-1);
488 }
489 }
490
491
push_parameter(short dist_param,unsigned int param)492 void BasisPolynomial::push_parameter(short dist_param, unsigned int param)
493 {
494 if (polyRep)
495 polyRep->push_parameter(dist_param, param);
496 else {
497 PCerr << "Error: push_parameter(unsigned int) not available for this basis "
498 << "polynomial type." << std::endl;
499 abort_handler(-1);
500 }
501 }
502
503
parameterized() const504 bool BasisPolynomial::parameterized() const
505 {
506 if (polyRep)
507 return polyRep->parameterized();
508 else
509 return false; // default if not overridden
510 }
511
512
collocation_rule(short rule)513 void BasisPolynomial::collocation_rule(short rule)
514 {
515 if (polyRep)
516 polyRep->collocation_rule(rule);
517 else {
518 PCerr << "Error: collocation_rule(short) not available for this basis "
519 << "polynomial type." << std::endl;
520 abort_handler(-1);
521 }
522 }
523
524
collocation_rule() const525 short BasisPolynomial::collocation_rule() const
526 {
527 if (!polyRep) {
528 PCerr << "Error: collocation_rule() not available for this basis "
529 << "polynomial type." << std::endl;
530 abort_handler(-1);
531 }
532 return polyRep->collocation_rule();
533 }
534
535
interpolation_size() const536 size_t BasisPolynomial::interpolation_size() const
537 {
538 if (!polyRep) {
539 PCerr << "Error: interpolation_size() not available for this basis "
540 << "polynomial type." << std::endl;
541 abort_handler(-1);
542 }
543 return polyRep->interpolation_size();
544 }
545
546
interpolation_points(const RealArray & interpolation_pts)547 void BasisPolynomial::interpolation_points(const RealArray& interpolation_pts)
548 {
549 if (polyRep)
550 polyRep->interpolation_points(interpolation_pts);
551 else {
552 PCerr << "Error: interpolation_points() not available for this basis "
553 << "polynomial type." << std::endl;
554 abort_handler(-1);
555 }
556 }
557
558
interpolation_points() const559 const RealArray& BasisPolynomial::interpolation_points() const
560 {
561 if (!polyRep) {
562 PCerr << "Error: interpolation_points() not available for this basis "
563 << "polynomial type." << std::endl;
564 abort_handler(-1);
565 }
566 return polyRep->interpolation_points();
567 }
568
569
length_scale() const570 Real BasisPolynomial::length_scale() const
571 {
572 if (!polyRep) {
573 PCerr << "Error: length_scale() not available for this basis polynomial "
574 << "type." << std::endl;
575 abort_handler(-1);
576 }
577 return polyRep->length_scale();
578 }
579
580
precompute_rules(unsigned short order)581 void BasisPolynomial::precompute_rules(unsigned short order)
582 {
583 if (polyRep)
584 polyRep->precompute_rules(order); // else no-op
585 //else {
586 // PCerr << "Error: precompute_rules() not available for this basis "
587 // << "polynomial type." << std::endl;
588 // abort_handler(-1);
589 //}
590 }
591
592 } // namespace Pecos
593