Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_StieltjesUnitTest.cpp
Go to the documentation of this file.
1// $Id$
2// $Source$
3// @HEADER
4// ***********************************************************************
5//
6// Stokhos Package
7// Copyright (2009) Sandia Corporation
8//
9// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10// license for use of this work by or on behalf of the U.S. Government.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
40//
41// ***********************************************************************
42// @HEADER
43
44#include "Teuchos_UnitTestHarness.hpp"
45#include "Teuchos_TestingHelpers.hpp"
46#include "Teuchos_UnitTestRepository.hpp"
47#include "Teuchos_GlobalMPISession.hpp"
48
49#include "Stokhos.hpp"
52
53// Quadrature functor to be passed into quadrature expansion for mapping
54// from Stieltjes basis back to original PCE
59
60 double operator() (const double& a) const {
61 vec[0] = a;
62 return pce.evaluate(vec);
63 }
66 mutable Teuchos::Array<double> vec;
67};
68
69// Class encapsulating setup of the Stieltjes basis for a given PCE
70// u = Func(x), where Func is specified by a template-parameter
71template <typename Func>
73 typedef typename Func::OrdinalType OrdinalType;
74 typedef typename Func::ValueType ValueType;
76 Func func;
79 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<OrdinalType,ValueType> > basis;
80 Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion<OrdinalType,ValueType> > exp;
81 Teuchos::RCP<const Stokhos::StieltjesPCEBasis<OrdinalType,ValueType> > st_1d_basis;
82 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<OrdinalType,ValueType> > st_basis;
83 Teuchos::RCP<const Stokhos::Quadrature<OrdinalType,ValueType> > st_quad;
85
86 Stieltjes_PCE_Setup(bool use_pce_quad_points_) :
87 func(), use_pce_quad_points(use_pce_quad_points_)
88 {
89 rtol = 1e-8;
90 atol = 1e-12;
91 const OrdinalType d = 3;
92 const OrdinalType p = 5;
93
94 // Create product basis
95 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<OrdinalType,ValueType> > > bases(d);
96 for (OrdinalType i=0; i<d; i++)
97 bases[i] =
99 basis =
101
102 // Create approximation
103 sz = basis->size();
105 for (OrdinalType i=0; i<d; i++)
106 x.term(i, 1) = 1.0;
107
108 // Tensor product quadrature
109 Teuchos::RCP<const Stokhos::Quadrature<OrdinalType,ValueType> > quad =
111
112 // Triple product tensor
113 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
114 basis->computeTripleProductTensor();
115
116 // Quadrature expansion
117 exp = Teuchos::rcp(new Stokhos::QuadOrthogPolyExpansion<OrdinalType,ValueType>(basis, Cijk, quad));
118
119 // Compute PCE via quadrature expansion
120 u.reset(basis);
121 v.reset(basis);
122 func.eval(*exp, x, u);
123 exp->times(v,u,u);
124
125 // Compute Stieltjes basis
126 st_1d_basis =
128 p, Teuchos::rcp(&u,false), quad, use_pce_quad_points));
129 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<OrdinalType,ValueType> > > st_bases(1);
130 st_bases[0] = st_1d_basis;
131 st_basis =
132 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<OrdinalType,ValueType>(st_bases, 1e-15));
133 st_sz = st_basis->size();
136 u_st[0] = u.mean();
137 u_st[1] = 1.0;
138
139 // Tensor product quadrature
140 st_quad =
142
143 // Triple product tensor
144 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > st_Cijk =
145 st_basis->computeTripleProductTensor();
146
147 // Quadrature expansion
149 st_Cijk,
150 st_quad);
151
152 st_exp.times(v_st, u_st, u_st);
153 }
154
155};
156
157//
158// Stieltjes tests based on expansion of u = cos(x) where x is a U([-1,1])
159// random variable
160//
162
163 template <typename Ordinal_Type, typename Value_Type>
165 typedef Ordinal_Type OrdinalType;
166 typedef Value_Type ValueType;
167 static const bool is_even = true;
168 void
174 };
175 Stieltjes_PCE_Setup< Stieltjes_Cos_Func<int,double> > setup(true);
176
177 // Tests mapping from Stieltjes basis to original is correct
178 TEUCHOS_UNIT_TEST( Stokhos_StieltjesPCEBasis, CosMap ) {
180 setup.st_1d_basis->transformCoeffsFromStieltjes(setup.u_st.coeff(),
181 u2.coeff());
182 success = Stokhos::comparePCEs(setup.u, "u", u2, "u2", setup.rtol, setup.atol, out);
183 }
184
185 // Tests Stieltjes basis is orthogonal
186 TEUCHOS_UNIT_TEST( Stokhos_StieltjesPCEBasis, CosOrthog ) {
187 const Teuchos::Array<double>& norms = setup.st_1d_basis->norm_squared();
188 const Teuchos::Array<double>& weights = setup.st_quad->getQuadWeights();
189 const Teuchos::Array< Teuchos::Array<double> >& values =
190 setup.st_quad->getBasisAtQuadPoints();
191 Teuchos::SerialDenseMatrix<int,double> mat(setup.st_sz,
192 setup.st_sz);
193 for (int i=0; i<setup.st_sz; i++) {
194 for (int j=0; j<setup.st_sz; j++) {
195 for (unsigned int k=0; k<weights.size(); k++)
196 mat(i,j) += weights[k]*values[k][i]*values[k][j];
197 mat(i,j) /= std::sqrt(norms[i]*norms[j]);
198 }
199 mat(i,i) -= 1.0;
200 }
201 success = mat.normInf() < setup.atol;
202 if (!success) {
203 out << "\n Error, mat.normInf() < atol = " << mat.normInf()
204 << " < " << setup.atol << ": failed!\n";
205 out << "mat = " << printMat(mat) << std::endl;
206 }
207 }
208
209 // Tests PCE computed from Stieltjes basis is same as original
210 TEUCHOS_UNIT_TEST( Stokhos_StieltjesPCEBasis, CosPCE ) {
212 stieltjes_pce_quad_func quad_func(setup.v_st, *setup.st_basis);
213 setup.exp->unary_op(quad_func, v2, setup.u);
214 success = comparePCEs(setup.v, "v", v2, "v2", setup.rtol, setup.atol, out);
215 }
216
217 // Tests mean computed from Stieltjes basis is same as original
218 TEUCHOS_UNIT_TEST( Stokhos_StieltjesPCEBasis, CosMean ) {
219 success = Teuchos::testRelErr("v.mean()", setup.v.mean(),
220 "v_st.mean()", setup.v_st.mean(),
221 "rtol", setup.rtol,
222 "rtol", setup.rtol,
223 Teuchos::Ptr<std::ostream>(out.getOStream().get()));
224
225 }
226
227 // Tests mean standard deviation from Stieltjes basis is same as original
228 TEUCHOS_UNIT_TEST( Stokhos_StieltjesPCEBasis, CosStandardDeviation ) {
229 success = Teuchos::testRelErr("v.standard_deviation()",
230 setup.v.standard_deviation(),
231 "v_st.standard_devaition()",
232 setup.v_st.standard_deviation(),
233 "rtol", 1e-1,
234 "rtol", 1e-1,
235 Teuchos::Ptr<std::ostream>(out.getOStream().get()));
236 }
237
238}
239
240//
241// Stieltjes tests based on expansion of u = sin(x) where x is a U([-1,1])
242// random variable
243//
245
246 template <typename Ordinal_Type, typename Value_Type>
248 typedef Ordinal_Type OrdinalType;
249 typedef Value_Type ValueType;
250 static const bool is_even = false;
251 void
257 };
258 Stieltjes_PCE_Setup< Stieltjes_Sin_Func<int,double> > setup(true);
259
260 // Tests mapping from Stieltjes basis to original is correct
261 TEUCHOS_UNIT_TEST( Stokhos_StieltjesPCEBasis, SinMap ) {
263 setup.st_1d_basis->transformCoeffsFromStieltjes(setup.u_st.coeff(),
264 u2.coeff());
265 success = Stokhos::comparePCEs(setup.u, "u", u2, "u2", setup.rtol, setup.atol, out);
266 }
267
268 // Tests Stieltjes basis is orthogonal
269 TEUCHOS_UNIT_TEST( Stokhos_StieltjesPCEBasis, SinOrthog ) {
270 const Teuchos::Array<double>& norms = setup.st_1d_basis->norm_squared();
271 const Teuchos::Array<double>& weights = setup.st_quad->getQuadWeights();
272 const Teuchos::Array< Teuchos::Array<double> >& values =
273 setup.st_quad->getBasisAtQuadPoints();
274 Teuchos::SerialDenseMatrix<int,double> mat(setup.st_sz,
275 setup.st_sz);
276 for (int i=0; i<setup.st_sz; i++) {
277 for (int j=0; j<setup.st_sz; j++) {
278 for (unsigned int k=0; k<weights.size(); k++)
279 mat(i,j) += weights[k]*values[k][i]*values[k][j];
280 mat(i,j) /= std::sqrt(norms[i]*norms[j]);
281 }
282 mat(i,i) -= 1.0;
283 }
284 success = mat.normInf() < setup.atol;
285 if (!success) {
286 out << "\n Error, mat.normInf() < atol = " << mat.normInf()
287 << " < " << setup.atol << ": failed!\n";
288 out << "mat = " << printMat(mat) << std::endl;
289 }
290 }
291
292 // Tests PCE computed from Stieltjes basis is same as original
293 TEUCHOS_UNIT_TEST( Stokhos_StieltjesPCEBasis, SinPCE ) {
295 stieltjes_pce_quad_func quad_func(setup.v_st, *setup.st_basis);
296 setup.exp->unary_op(quad_func, v2, setup.u);
297 success = comparePCEs(setup.v, "v", v2, "v2", setup.rtol, setup.atol, out);
298 }
299
300 // Tests mean computed from Stieltjes basis is same as original
301 TEUCHOS_UNIT_TEST( Stokhos_StieltjesPCEBasis, SinMean ) {
302 success = Teuchos::testRelErr("v.mean()", setup.v.mean(),
303 "v_st.mean()", setup.v_st.mean(),
304 "rtol", setup.rtol,
305 "rtol", setup.rtol,
306 Teuchos::Ptr<std::ostream>(out.getOStream().get()));
307 }
308
309 // Tests mean standard deviation from Stieltjes basis is same as original
310 TEUCHOS_UNIT_TEST( Stokhos_StieltjesPCEBasis, SinStandardDeviation ) {
311 success = Teuchos::testRelErr("v.standard_deviation()",
312 setup.v.standard_deviation(),
313 "v_st.standard_devaition()",
314 setup.v_st.standard_deviation(),
315 "rtol", 1e-1,
316 "rtol", 1e-1,
317 Teuchos::Ptr<std::ostream>(out.getOStream().get()));
318 }
319
320}
321
322//
323// Stieltjes tests based on expansion of u = exp(x) where x is a U([-1,1])
324// random variable. For this test we don't use the PCE quad points and
325// instead use those generated for the Stieltjes basis
326//
328
329 template <typename Ordinal_Type, typename Value_Type>
331 typedef Ordinal_Type OrdinalType;
332 typedef Value_Type ValueType;
333 static const bool is_even = false;
334 void
340 };
341 Stieltjes_PCE_Setup< Stieltjes_Exp_Func<int,double> > setup(false);
342
343 // Tests mapping from Stieltjes basis to original is correct
344 TEUCHOS_UNIT_TEST( Stokhos_StieltjesPCEBasis, ExpMap ) {
346 setup.st_1d_basis->transformCoeffsFromStieltjes(setup.u_st.coeff(),
347 u2.coeff());
348 success = Stokhos::comparePCEs(setup.u, "u", u2, "u2", setup.rtol, setup.atol, out);
349 }
350
351 // Tests Stieltjes basis is orthogonal
352 TEUCHOS_UNIT_TEST( Stokhos_StieltjesPCEBasis, ExpOrthog ) {
353 const Teuchos::Array<double>& norms = setup.st_1d_basis->norm_squared();
354 const Teuchos::Array<double>& weights = setup.st_quad->getQuadWeights();
355 const Teuchos::Array< Teuchos::Array<double> >& values =
356 setup.st_quad->getBasisAtQuadPoints();
357 Teuchos::SerialDenseMatrix<int,double> mat(setup.st_sz,
358 setup.st_sz);
359 for (int i=0; i<setup.st_sz; i++) {
360 for (int j=0; j<setup.st_sz; j++) {
361 for (unsigned int k=0; k<weights.size(); k++)
362 mat(i,j) += weights[k]*values[k][i]*values[k][j];
363 mat(i,j) /= std::sqrt(norms[i]*norms[j]);
364 }
365 mat(i,i) -= 1.0;
366 }
367 success = mat.normInf() < setup.atol;
368 if (!success) {
369 out << "\n Error, mat.normInf() < atol = " << mat.normInf()
370 << " < " << setup.atol << ": failed!\n";
371 out << "mat = " << printMat(mat) << std::endl;
372 }
373 }
374
375 // Tests PCE computed from Stieltjes basis is same as original
376 TEUCHOS_UNIT_TEST( Stokhos_StieltjesPCEBasis, ExpPCE ) {
378 stieltjes_pce_quad_func quad_func(setup.v_st, *setup.st_basis);
379 setup.exp->unary_op(quad_func, v2, setup.u);
380 success = comparePCEs(setup.v, "v", v2, "v2", setup.rtol, setup.atol, out);
381 }
382
383 // Tests mean computed from Stieltjes basis is same as original
384 TEUCHOS_UNIT_TEST( Stokhos_StieltjesPCEBasis, ExpMean ) {
385 success = Teuchos::testRelErr("v.mean()", setup.v.mean(),
386 "v_st.mean()", setup.v_st.mean(),
387 "rtol", setup.rtol,
388 "rtol", setup.rtol,
389 Teuchos::Ptr<std::ostream>(out.getOStream().get()));
390 }
391
392 // Tests mean standard deviation from Stieltjes basis is same as original
393 TEUCHOS_UNIT_TEST( Stokhos_StieltjesPCEBasis, ExpStandardDeviation ) {
394 success = Teuchos::testRelErr("v.standard_deviation()",
395 setup.v.standard_deviation(),
396 "v_st.standard_devaition()",
397 setup.v_st.standard_deviation(),
398 "rtol", 1e-1,
399 "rtol", 1e-1,
400 Teuchos::Ptr<std::ostream>(out.getOStream().get()));
401 }
402
403}
404
405int main( int argc, char* argv[] ) {
406 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
407 return Teuchos::UnitTestRepository::runUnitTestsFromMain(argc, argv);
408}
int main(int argc, char *argv[])
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
void reset(const Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > &new_basis, ordinal_type sz=0)
Reset to a new basis.
value_type mean() const
Compute mean of expansion.
pointer coeff()
Return coefficient array.
Abstract base class for multivariate orthogonal polynomials.
Orthogonal polynomial expansions based on numerical quadrature.
Generates three-term recurrence using the Discretized Stieltjes procedure applied to a polynomial cha...
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
TEUCHOS_UNIT_TEST(Stokhos_StieltjesPCEBasis, CosMap)
Stieltjes_PCE_Setup< Stieltjes_Cos_Func< int, double > > setup(true)
TEUCHOS_UNIT_TEST(Stokhos_StieltjesPCEBasis, ExpMap)
Stieltjes_PCE_Setup< Stieltjes_Exp_Func< int, double > > setup(false)
Stieltjes_PCE_Setup< Stieltjes_Sin_Func< int, double > > setup(true)
TEUCHOS_UNIT_TEST(Stokhos_StieltjesPCEBasis, SinMap)
bool comparePCEs(const PCEType &a1, const std::string &a1_name, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &a2, const std::string &a2_name, const ValueType &rel_tol, const ValueType &abs_tol, Teuchos::FancyOStream &out)
void eval(Stokhos::QuadOrthogPolyExpansion< OrdinalType, ValueType > &exp, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &x, Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &u)
void eval(Stokhos::QuadOrthogPolyExpansion< OrdinalType, ValueType > &exp, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &x, Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &u)
void eval(Stokhos::QuadOrthogPolyExpansion< OrdinalType, ValueType > &exp, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &x, Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &u)
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > u
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > v
Teuchos::RCP< const Stokhos::Quadrature< OrdinalType, ValueType > > st_quad
Teuchos::RCP< const Stokhos::CompletePolynomialBasis< OrdinalType, ValueType > > st_basis
Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion< OrdinalType, ValueType > > exp
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > v_st
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > u_st
Teuchos::RCP< const Stokhos::CompletePolynomialBasis< OrdinalType, ValueType > > basis
Stieltjes_PCE_Setup(bool use_pce_quad_points_)
Teuchos::RCP< const Stokhos::StieltjesPCEBasis< OrdinalType, ValueType > > st_1d_basis
const Stokhos::OrthogPolyApprox< int, double > & pce
stieltjes_pce_quad_func(const Stokhos::OrthogPolyApprox< int, double > &pce_, const Stokhos::OrthogPolyBasis< int, double > &basis_)
double operator()(const double &a) const
const Stokhos::OrthogPolyBasis< int, double > & basis
void printMat(const char *name, Epetra_IntSerialDenseMatrix &matrix)