Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
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
58  pce(pce_), basis(basis_), vec(1) {}
59 
60  double operator() (const double& a) const {
61  vec[0] = a;
62  return pce.evaluate(vec);
63  }
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
71 template <typename Func>
73  typedef typename Func::OrdinalType OrdinalType;
74  typedef typename Func::ValueType ValueType;
76  Func func;
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
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
111 
112  // Triple product tensor
115 
116  // Quadrature expansion
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));
130  st_bases[0] = st_1d_basis;
131  st_basis =
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
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 //
161 namespace StieltjesCosTest {
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
172  exp.cos(u,x);
173  }
174  };
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();
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 = " << 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 //
244 namespace StieltjesSinTest {
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
255  exp.sin(u,x);
256  }
257  };
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();
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 = " << 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 //
327 namespace StieltjesExpTest {
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
338  exp.exp(u,x);
339  }
340  };
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();
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 = " << 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 
405 int main( int argc, char* argv[] ) {
406  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
408 }
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
Teuchos::RCP< const Stokhos::Quadrature< OrdinalType, ValueType > > st_quad
void sin(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void eval(Stokhos::QuadOrthogPolyExpansion< OrdinalType, ValueType > &exp, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &x, Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &u)
void times(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
value_type evaluate(const Teuchos::Array< value_type > &point) const
Evaluate polynomial approximation at a point.
TEUCHOS_UNIT_TEST(Stokhos_StieltjesPCEBasis, ExpMap)
void exp(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Stieltjes_PCE_Setup< Stieltjes_Sin_Func< int, double > > setup(true)
void cos(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
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)
pointer coeff()
Return coefficient array.
void eval(Stokhos::QuadOrthogPolyExpansion< OrdinalType, ValueType > &exp, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &x, Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &u)
Stieltjes_PCE_Setup< Stieltjes_Cos_Func< int, double > > setup(true)
static int runUnitTestsFromMain(int argc, char *argv[])
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void reset(const Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > &new_basis, ordinal_type sz=0)
Reset to a new basis.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
Stieltjes_PCE_Setup(bool use_pce_quad_points_)
double operator()(const double &a) const
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > u_st
TEUCHOS_UNIT_TEST(Stokhos_StieltjesPCEBasis, CosMap)
Teuchos::RCP< const Stokhos::CompletePolynomialBasis< OrdinalType, ValueType > > st_basis
value_type mean() const
Compute mean of expansion.
KOKKOS_INLINE_FUNCTION PCE< Storage > exp(const PCE< Storage > &a)
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > u
int main(int argc, char **argv)
ScalarTraits< ScalarType >::magnitudeType normInf() const
Teuchos::Array< double > vec
Teuchos::RCP< const Stokhos::StieltjesPCEBasis< OrdinalType, ValueType > > st_1d_basis
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > v
Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion< OrdinalType, ValueType > > exp
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > v_st
const Stokhos::OrthogPolyBasis< int, double > & basis
size_type size() const
stieltjes_pce_quad_func(const Stokhos::OrthogPolyApprox< int, double > &pce_, const Stokhos::OrthogPolyBasis< int, double > &basis_)
virtual ordinal_type size() const
Return total size of basis.
TEUCHOS_UNIT_TEST(Stokhos_StieltjesPCEBasis, SinMap)
const Stokhos::OrthogPolyApprox< int, double > & pce
bool testRelErr(const std::string &v1_name, const Scalar &v1, const std::string &v2_name, const Scalar &v2, const std::string &maxRelErr_error_name, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType &maxRelErr_error, const std::string &maxRelErr_warning_name, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType &maxRelErr_warning, const Ptr< std::ostream > &out)
Teuchos::RCP< const Stokhos::CompletePolynomialBasis< OrdinalType, ValueType > > basis
Stieltjes_PCE_Setup< Stieltjes_Exp_Func< int, double > > setup(false)
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules...
void eval(Stokhos::QuadOrthogPolyExpansion< OrdinalType, ValueType > &exp, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &x, Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &u)
reference term(ordinal_type dimension, ordinal_type order)
Get coefficient term for given dimension and order.