Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
stieltjes_example.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 <iostream>
45#include <iomanip>
46
47#include "Stokhos.hpp"
52
54
55struct pce_quad_func {
59 pce(pce_), basis(basis_), vec(2) {}
60
61 double operator() (const double& a, const double& b) const {
62 vec[0] = a;
63 vec[1] = b;
64 return pce.evaluate(vec);
65 }
68 mutable Teuchos::Array<double> vec;
69};
70
71double rel_err(double a, double b) {
72 return std::abs(a-b)/std::abs(b);
73}
74
75int main(int argc, char **argv)
76{
77 try {
78
79 const unsigned int d = 2;
80 const unsigned int pmin = 1;
81 const unsigned int pmax = 15;
82 const unsigned int np = pmax-pmin+1;
83 bool use_pce_quad_points = false;
84 bool normalize = false;
85 bool project_integrals = false;
86 bool lanczos = false;
87 bool sparse_grid = true;
88#ifndef HAVE_STOKHOS_DAKOTA
89 sparse_grid = false;
90#endif
91 Teuchos::Array<double> mean(np), mean_st(np), std_dev(np), std_dev_st(np);
92 Teuchos::Array<double> pt(np), pt_st(np);
93
94 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
95 Teuchos::Array<double> eval_pt(d, 0.5);
96 double pt_true;
97
98 // Loop over orders
99 unsigned int n = 0;
100 for (unsigned int p=pmin; p<=pmax; p++) {
101
102 std::cout << "p = " << p << std::endl;
103
104 // Create product basis
105 for (unsigned int i=0; i<d; i++)
106 bases[i] = Teuchos::rcp(new basis_type(p));
107 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
108 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));
109
110 // Create approximation
111 Stokhos::OrthogPolyApprox<int,double> x(basis), u(basis), v(basis),
112 w(basis), w2(basis);
113 for (unsigned int i=0; i<d; i++) {
114 x.term(i, 1) = 1.0;
115 }
116
117 double x_pt = x.evaluate(eval_pt);
118 pt_true = std::sin(x_pt)/(1.0 + exp(x_pt));
119
120 // Quadrature
121 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad;
122#ifdef HAVE_STOKHOS_DAKOTA
123 if (sparse_grid)
124 quad =
125 Teuchos::rcp(new Stokhos::SparseGridQuadrature<int,double>(basis, p));
126#endif
127 if (!sparse_grid)
128 quad =
129 Teuchos::rcp(new Stokhos::TensorProductQuadrature<int,double>(basis));
130
131 // Triple product tensor
132 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
133 basis->computeTripleProductTensor();
134
135 // Quadrature expansion
136 Stokhos::QuadOrthogPolyExpansion<int,double> quad_exp(basis, Cijk, quad);
137
138 // Compute PCE via quadrature expansion
139 quad_exp.sin(u,x);
140 quad_exp.exp(v,x);
141 quad_exp.plusEqual(v, 1.0);
142 quad_exp.divide(w,u,v);
143
144 // Compute Stieltjes basis
145 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > st_bases(2);
146 Teuchos::RCP< Stokhos::LanczosProjPCEBasis<int,double> > stp_basis_u, stp_basis_v;
147 Teuchos::RCP< Stokhos::LanczosPCEBasis<int,double> > st_basis_u, st_basis_v;
148 if (lanczos) {
149 if (project_integrals) {
150 stp_basis_u =
152 p, Teuchos::rcp(&u,false), Cijk, normalize, true));
153 stp_basis_v =
155 p, Teuchos::rcp(&v,false), Cijk, normalize, true));
156 st_bases[0] = stp_basis_u;
157 st_bases[1] = stp_basis_v;
158 }
159 else {
160 st_basis_u =
162 p, Teuchos::rcp(&u,false), quad, normalize, true));
163 st_basis_v =
165 p, Teuchos::rcp(&v,false), quad, normalize, true));
166 st_bases[0] = st_basis_u;
167 st_bases[1] = st_basis_v;
168 }
169 }
170 else {
171 st_bases[0] =
173 p, Teuchos::rcp(&u,false), quad, use_pce_quad_points,
174 normalize, project_integrals, Cijk));
175 st_bases[1] =
177 p, Teuchos::rcp(&v,false), quad, use_pce_quad_points,
178 normalize, project_integrals, Cijk));
179 }
180
181 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> >
182 st_basis =
183 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(st_bases));
184 //std::cout << *st_basis << std::endl;
185
186 Stokhos::OrthogPolyApprox<int,double> u_st(st_basis), v_st(st_basis),
187 w_st(st_basis);
188 if (lanczos) {
189 if (project_integrals) {
190 u_st.term(0, 0) = stp_basis_u->getNewCoeffs(0);
191 u_st.term(0, 1) = stp_basis_u->getNewCoeffs(1);
192 v_st.term(0, 0) = stp_basis_v->getNewCoeffs(0);
193 v_st.term(1, 1) = stp_basis_v->getNewCoeffs(1);
194 }
195 else {
196 u_st.term(0, 0) = st_basis_u->getNewCoeffs(0);
197 u_st.term(0, 1) = st_basis_u->getNewCoeffs(1);
198 v_st.term(0, 0) = st_basis_v->getNewCoeffs(0);
199 v_st.term(1, 1) = st_basis_v->getNewCoeffs(1);
200 }
201 }
202 else {
203 u_st.term(0, 0) = u.mean();
204 u_st.term(0, 1) = 1.0;
205 v_st.term(0, 0) = v.mean();
206 v_st.term(1, 1) = 1.0;
207 }
208
209 // Triple product tensor
210 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > st_Cijk =
211 st_basis->computeTripleProductTensor();
212
213 // Tensor product quadrature
214 Teuchos::RCP<const Stokhos::Quadrature<int,double> > st_quad;
215 if (!use_pce_quad_points) {
216#ifdef HAVE_STOKHOS_DAKOTA
217 if (sparse_grid)
218 st_quad = Teuchos::rcp(new Stokhos::SparseGridQuadrature<int,double>(st_basis, p));
219#endif
220 if (!sparse_grid)
221 st_quad = Teuchos::rcp(new Stokhos::TensorProductQuadrature<int,double>(st_basis));
222 }
223 else {
224 Teuchos::Array<double> st_points_0;
225 Teuchos::Array<double> st_weights_0;
226 Teuchos::Array< Teuchos::Array<double> > st_values_0;
227 st_bases[0]->getQuadPoints(p+1, st_points_0, st_weights_0, st_values_0);
228 Teuchos::Array<double> st_points_1;
229 Teuchos::Array<double> st_weights_1;
230 Teuchos::Array< Teuchos::Array<double> > st_values_1;
231 st_bases[1]->getQuadPoints(p+1, st_points_1, st_weights_1, st_values_1);
232 Teuchos::RCP< Teuchos::Array< Teuchos::Array<double> > > st_points =
233 Teuchos::rcp(new Teuchos::Array< Teuchos::Array<double> >(st_points_0.size()));
234 for (int i=0; i<st_points_0.size(); i++) {
235 (*st_points)[i].resize(2);
236 (*st_points)[i][0] = st_points_0[i];
237 (*st_points)[i][1] = st_points_1[i];
238 }
239 Teuchos::RCP< Teuchos::Array<double> > st_weights =
240 Teuchos::rcp(new Teuchos::Array<double>(st_weights_0));
241 Teuchos::RCP< const Stokhos::OrthogPolyBasis<int,double> > st_b =
242 st_basis;
243 st_quad =
244 Teuchos::rcp(new Stokhos::UserDefinedQuadrature<int,double>(st_b,
245 st_points,
246 st_weights));
247 }
248
249 // Quadrature expansion
251 st_Cijk,
252 st_quad);
253
254 // Compute w_st = u_st*v_st in Stieltjes basis
255 st_quad_exp.divide(w_st, u_st, v_st);
256
257 // Project w_st back to original basis
258 pce_quad_func st_func(w_st, *st_basis);
259 quad_exp.binary_op(st_func, w2, u, v);
260
261 // std::cout.precision(12);
262 // std::cout << w;
263 // std::cout << w2;
264 // std::cout << w_st;
265 mean[n] = w.mean();
266 mean_st[n] = w2.mean();
267 std_dev[n] = w.standard_deviation();
268 std_dev_st[n] = w2.standard_deviation();
269 pt[n] = w.evaluate(eval_pt);
270 pt_st[n] = w2.evaluate(eval_pt);
271 n++;
272 }
273
274 n = 0;
275 int wi=10;
276 std::cout << "Statistical error:" << std::endl;
277 std::cout << "p "
278 << std::setw(wi) << "mean" << " "
279 << std::setw(wi) << "mean_st" << " "
280 << std::setw(wi) << "std_dev" << " "
281 << std::setw(wi) << "std_dev_st" << " "
282 << std::setw(wi) << "point" << " "
283 << std::setw(wi) << "point_st" << std::endl;
284 for (unsigned int p=pmin; p<pmax; p++) {
285 std::cout.precision(3);
286 std::cout.setf(std::ios::scientific);
287 std::cout << p << " "
288 << std::setw(wi) << rel_err(mean[n], mean[np-1]) << " "
289 << std::setw(wi) << rel_err(mean_st[n], mean[np-1]) << " "
290 << std::setw(wi) << rel_err(std_dev[n], std_dev[np-1]) << " "
291 << std::setw(wi) << rel_err(std_dev_st[n], std_dev[np-1])
292 << " "
293 << std::setw(wi) << rel_err(pt[n], pt_true) << " "
294 << std::setw(wi) << rel_err(pt_st[n], pt_true)
295 << std::endl;
296 n++;
297 }
298
299 }
300 catch (std::exception& e) {
301 std::cout << e.what() << std::endl;
302 }
303}
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
Legendre polynomial basis.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
value_type mean() const
Compute mean of expansion.
value_type standard_deviation() const
Compute standard deviation of expansion.
value_type evaluate(const Teuchos::Array< value_type > &point) const
Evaluate polynomial approximation at a point.
Abstract base class for multivariate orthogonal polynomials.
void plusEqual(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const value_type &x)
Orthogonal polynomial expansions based on numerical quadrature.
void exp(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void divide(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)
void binary_op(const FuncT &func, 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)
Nonlinear binary function.
void sin(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
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.
int main(int argc, char **argv)
double rel_err(double a, double b)
Stokhos::LegendreBasis< int, double > basis_type
const Stokhos::OrthogPolyApprox< int, double > & pce
double operator()(const double &a, const double &b) const
const Stokhos::OrthogPolyBasis< int, double > & basis
pce_quad_func(const Stokhos::OrthogPolyApprox< int, double > &pce_, const Stokhos::OrthogPolyBasis< int, double > &basis_)
Teuchos::Array< double > vec