Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_GramSchmidtUnitTest.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"
54
55// Quadrature functor to be passed into quadrature expansion for mapping
56// from Gram-Schmidt basis back to original PCE
62
63 double operator() (const double& a) const {
64 vec[0] = a;
65 return pce.evaluate(vec);
66 }
69 mutable Teuchos::Array<double> vec;
70};
76
77 double operator() (const double& a, const double& b) const {
78 vec[0] = a;
79 vec[1] = b;
80 return pce.evaluate(vec);
81 }
84 mutable Teuchos::Array<double> vec;
85};
86
87// Class encapsulating setup of the Gram-Schmidt basis for a given PCE
88template <typename OrdinalType, typename ValueType>
90 ValueType rtol, atol;
91 OrdinalType sz, gs_sz;
92 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<OrdinalType,ValueType> > basis;
93 Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion<OrdinalType,ValueType> > exp;
94 Teuchos::RCP<const Stokhos::GramSchmidtBasis<OrdinalType,ValueType> > gs_basis;
95 Teuchos::RCP<const Stokhos::Quadrature<OrdinalType,ValueType> > gs_quad;
97
99 {
100 rtol = 1e-7;
101 atol = 1e-12;
102 const OrdinalType d = 3;
103 const OrdinalType p = 5;
104
105 // Create product basis
106 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<OrdinalType,ValueType> > > bases(d);
107 for (OrdinalType i=0; i<d; i++)
108 bases[i] =
110 basis =
111 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<OrdinalType,ValueType>(bases, 1e-15));
112
113 // Create approximation
114 sz = basis->size();
116 for (OrdinalType i=0; i<d; i++)
117 x.term(i, 1) = 1.0;
118
119 // Tensor product quadrature
120 Teuchos::RCP<const Stokhos::Quadrature<OrdinalType,ValueType> > quad =
122
123 // Triple product tensor
124 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
125 basis->computeTripleProductTensor();
126
127 // Quadrature expansion
128 Teuchos::RCP<Teuchos::ParameterList> exp_params =
129 Teuchos::rcp(new Teuchos::ParameterList);
130 exp_params->set("Use Quadrature for Times", true);
131 exp = Teuchos::rcp(new Stokhos::QuadOrthogPolyExpansion<OrdinalType,ValueType>(basis, Cijk, quad, exp_params));
132
133 // Compute PCE via quadrature expansion
134 u.reset(basis);
135 v.reset(basis);
136 w.reset(basis);
137 exp->sin(u,x);
138 exp->exp(v,x);
139 exp->times(w,u,v);
140
141 // Compute Stieltjes basis
142 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<OrdinalType,ValueType> > > st_bases(2);
143 st_bases[0] = Teuchos::rcp(new Stokhos::StieltjesPCEBasis<OrdinalType,ValueType>(p, Teuchos::rcp(&u,false), quad, true));
144 st_bases[1] = Teuchos::rcp(new Stokhos::StieltjesPCEBasis<OrdinalType,ValueType>(p, Teuchos::rcp(&v,false), quad, true));
145 Teuchos::RCP<const Stokhos::OrthogPolyBasis<OrdinalType,ValueType> > st_basis =
146 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<OrdinalType,ValueType>(st_bases, 1e-15));
147 Stokhos::OrthogPolyApprox<OrdinalType,ValueType> u_st(st_basis), v_st(st_basis);
148 u_st.term(0, 0) = u.mean();
149 u_st.term(0, 1) = 1.0;
150 v_st.term(0, 0) = v.mean();
151 v_st.term(1, 1) = 1.0;
152
153 // Use Gram-Schmidt to orthogonalize Stieltjes basis of u and v
154 Teuchos::Array<ValueType> st_points_0;
155 Teuchos::Array<ValueType> st_weights_0;
156 Teuchos::Array< Teuchos::Array<ValueType> > st_values_0;
157 st_bases[0]->getQuadPoints(p+1, st_points_0, st_weights_0, st_values_0);
158 Teuchos::Array<ValueType> st_points_1;
159 Teuchos::Array<ValueType> st_weights_1;
160 Teuchos::Array< Teuchos::Array<ValueType> > st_values_1;
161 st_bases[1]->getQuadPoints(p+1, st_points_1, st_weights_1, st_values_1);
162 Teuchos::RCP< Teuchos::Array< Teuchos::Array<ValueType> > > st_points =
163 Teuchos::rcp(new Teuchos::Array< Teuchos::Array<ValueType> >(st_points_0.size()));
164 for (int i=0; i<st_points_0.size(); i++) {
165 (*st_points)[i].resize(2);
166 (*st_points)[i][0] = st_points_0[i];
167 (*st_points)[i][1] = st_points_1[i];
168 }
169 Teuchos::RCP< Teuchos::Array<ValueType> > st_weights =
170 Teuchos::rcp(new Teuchos::Array<ValueType>);
171 *st_weights = st_weights_0;
172
173 gs_basis = Teuchos::rcp(new Stokhos::GramSchmidtBasis<OrdinalType,ValueType>(st_basis, *st_points, *st_weights, 1e-15));
174 gs_sz = gs_basis->size();
175
176 // Triple product tensor
177 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > gs_Cijk =
178 gs_basis->computeTripleProductTensor();
179
180 // Create quadrature for Gram-Schmidt basis using quad points and
181 // and weights from original basis mapped to Stieljtes basis
182 Teuchos::RCP< const Teuchos::Array< Teuchos::Array<ValueType> > > points =
183 st_points;
184 Teuchos::RCP< const Teuchos::Array<ValueType> > weights = st_weights;
185 gs_quad = Teuchos::rcp(new Stokhos::UserDefinedQuadrature<OrdinalType,ValueType>(gs_basis, points, weights));
186
187 // Gram-Schmidt quadrature expansion
189 gs_Cijk,
190 gs_quad,
191 exp_params);
192
196
197 // Map expansion in Stieltjes basis to Gram-Schmidt basis
198 gs_basis->transformCoeffs(u_st.coeff(), u_gs.coeff());
199 gs_basis->transformCoeffs(v_st.coeff(), v_gs.coeff());
200
201 // Compute w_gs = u_gs*v_gs in Gram-Schmidt basis
202 gs_exp.times(w_gs, u_gs, v_gs);
203 }
204
205};
206
208 GramSchmidt_PCE_Setup<int,double> setup;
209
210 // Tests mapping from Gram-Schmidt basis to original is correct
211 TEUCHOS_UNIT_TEST( Stokhos_GramSchmidtBasis, Map ) {
213 gram_schmidt_pce_binary_quad_func quad_func(setup.u_gs, *setup.gs_basis);
214 setup.exp->binary_op(quad_func, u2, setup.u, setup.v);
215 success = comparePCEs(setup.u, "u", u2, "u2", setup.rtol, setup.atol, out);
216 }
217
218 // Tests Gram-Schmidt basis is orthogonal
219 TEUCHOS_UNIT_TEST( Stokhos_GramSchmidtBasis, Orthog ) {
220 const Teuchos::Array<double>& norms = setup.gs_basis->norm_squared();
221 const Teuchos::Array<double>& weights = setup.gs_quad->getQuadWeights();
222 const Teuchos::Array< Teuchos::Array<double> >& values =
223 setup.gs_quad->getBasisAtQuadPoints();
224 Teuchos::SerialDenseMatrix<int,double> mat(setup.gs_sz, setup.gs_sz);
225 for (int i=0; i<setup.gs_sz; i++) {
226 for (int j=0; j<setup.gs_sz; j++) {
227 for (int k=0; k<weights.size(); k++)
228 mat(i,j) += weights[k]*values[k][i]*values[k][j];
229 mat(i,j) /= std::sqrt(norms[i]*norms[j]);
230 }
231 mat(i,i) -= 1.0;
232 }
233 double tol = 1e-4;
234 success = mat.normInf() < tol;
235 if (!success) {
236 out << "\n Error, mat.normInf() < tol = " << mat.normInf()
237 << " < " << tol << ": failed!\n";
238 out << "mat = " << printMat(mat) << std::endl;
239 }
240 }
241
242 // Tests PCE computed from Gram-Schmidt basis is same as original
243 TEUCHOS_UNIT_TEST( Stokhos_GramSchmidtBasis, PCE ) {
245 gram_schmidt_pce_binary_quad_func quad_func(setup.w_gs, *setup.gs_basis);
246 setup.exp->binary_op(quad_func, w2, setup.u, setup.v);
247 success = comparePCEs(setup.w, "w", w2, "w2", setup.rtol, setup.atol, out);
248 }
249
250 // Tests mean computed from Gram-Schmidt basis is same as original
251 TEUCHOS_UNIT_TEST( Stokhos_GramSchmidtBasis, Mean ) {
252 success = Teuchos::testRelErr("w.mean()", setup.w.mean(),
253 "w_gs.mean()", setup.w_gs.mean(),
254 "rtol", setup.rtol,
255 "rtol", setup.rtol,
256 Teuchos::Ptr<std::ostream>(out.getOStream().get()));
257 }
258
259 // Tests mean standard deviation from Gram-Schmidt basis is same as original
260 TEUCHOS_UNIT_TEST( Stokhos_GramSchmidtBasis, StandardDeviation ) {
261 success = Teuchos::testRelErr("w.standard_deviation()",
262 setup.w.standard_deviation(),
263 "w_gs.standard_devaition()",
264 setup.w_gs.standard_deviation(),
265 "rtol", 1e-3,
266 "rtol", 1e-3,
267 Teuchos::Ptr<std::ostream>(out.getOStream().get()));
268 }
269
270}
271
272int main( int argc, char* argv[] ) {
273 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
274 return Teuchos::UnitTestRepository::runUnitTestsFromMain(argc, argv);
275}
int main(int argc, char *argv[])
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Transforms a non-orthogonal multivariate basis to an orthogonal one using the Gram-Schmit procedure.
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.
reference term(ordinal_type dimension, ordinal_type order)
Get coefficient term for given dimension and order.
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_GramSchmidtBasis, Map)
GramSchmidt_PCE_Setup< int, double > setup
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > v_gs
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > u_gs
Teuchos::RCP< const Stokhos::CompletePolynomialBasis< OrdinalType, ValueType > > basis
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > w_gs
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > u
Teuchos::RCP< const Stokhos::Quadrature< OrdinalType, ValueType > > gs_quad
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > w
Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion< OrdinalType, ValueType > > exp
Teuchos::RCP< const Stokhos::GramSchmidtBasis< OrdinalType, ValueType > > gs_basis
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > v
gram_schmidt_pce_binary_quad_func(const Stokhos::OrthogPolyApprox< int, double > &pce_, const Stokhos::OrthogPolyBasis< int, double > &basis_)
const Stokhos::OrthogPolyApprox< int, double > & pce
double operator()(const double &a, const double &b) const
const Stokhos::OrthogPolyBasis< int, double > & basis
gram_schmidt_pce_unary_quad_func(const Stokhos::OrthogPolyApprox< int, double > &pce_, const Stokhos::OrthogPolyBasis< int, double > &basis_)
const Stokhos::OrthogPolyApprox< int, double > & pce
const Stokhos::OrthogPolyBasis< int, double > & basis
void printMat(const char *name, Epetra_IntSerialDenseMatrix &matrix)