Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
stieltjes_example2.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
42#include <iostream>
43#include <iomanip>
44
45#include "Stokhos.hpp"
48
49//typedef Stokhos::LegendreBasis<int,double> basis_type;
51
52struct pce_quad_func {
56 pce(pce_), basis(basis_), vec(2) {}
57
58 double operator() (const double& a, const double& b) const {
59 vec[0] = a;
60 vec[1] = b;
61 return pce.evaluate(vec);
62 }
65 mutable Teuchos::Array<double> vec;
66};
67
68struct sin_func {
70 pce(pce_) {}
71 double operator() (const Teuchos::Array<double>& x) const {
72 return std::sin(pce.evaluate(x));
73 }
75};
76
77struct exp_func {
79 pce(pce_) {}
80 double operator() (const Teuchos::Array<double>& x) const {
81 return 1.0 + std::exp(pce.evaluate(x));
82 }
84};
85
86double rel_err(double a, double b) {
87 return std::abs(a-b)/std::abs(b);
88}
89
90int main(int argc, char **argv)
91{
92 try {
93
94 const unsigned int d = 2;
95 const unsigned int pmin = 1;
96 const unsigned int pmax = 10;
97 const unsigned int np = pmax-pmin+1;
98 bool use_pce_quad_points = false;
99 bool normalize = false;
100 bool sparse_grid = true;
101#ifndef HAVE_STOKHOS_DAKOTA
102 sparse_grid = false;
103#endif
104 Teuchos::Array<double> mean(np), mean_st(np), std_dev(np), std_dev_st(np);
105 Teuchos::Array<double> pt(np), pt_st(np);
106
107 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
108 Teuchos::Array<double> eval_pt(d, 0.5);
109 double pt_true;
110
111 // Loop over orders
112 unsigned int n = 0;
113 for (unsigned int p=pmin; p<=pmax; p++) {
114
115 std::cout << "p = " << p << std::endl;
116
117 // Create product basis
118 for (unsigned int i=0; i<d; i++)
119 bases[i] = Teuchos::rcp(new basis_type(p));
120 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
121 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));
122
123 // Create approximation
124 Stokhos::OrthogPolyApprox<int,double> x(basis), u(basis), v(basis),
125 w(basis), w2(basis);
126 for (unsigned int i=0; i<d; i++) {
127 x.term(i, 1) = 1.0;
128 }
129
130 double x_pt = x.evaluate(eval_pt);
131 pt_true = std::sin(x_pt)/(1.0 + exp(x_pt));
132
133 // Quadrature
134 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad;
135#ifdef HAVE_STOKHOS_DAKOTA
136 if (sparse_grid)
137 quad =
138 Teuchos::rcp(new Stokhos::SparseGridQuadrature<int,double>(basis, p));
139#endif
140 if (!sparse_grid)
141 quad =
142 Teuchos::rcp(new Stokhos::TensorProductQuadrature<int,double>(basis));
143
144 // Triple product tensor
145 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
146 basis->computeTripleProductTensor();
147
148 // Quadrature expansion
149 Stokhos::QuadOrthogPolyExpansion<int,double> quad_exp(basis, Cijk, quad);
150
151 // Compute PCE via quadrature expansion
152 quad_exp.sin(u,x);
153 quad_exp.exp(v,x);
154 quad_exp.plusEqual(v, 1.0);
155 quad_exp.divide(w,u,v);
156
157 // Compute Stieltjes basis
158 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > st_bases(2);
159 Teuchos::RCP<sin_func> sf = Teuchos::rcp(new sin_func(x));
160 Teuchos::RCP<exp_func> ef = Teuchos::rcp(new exp_func(x));
161 st_bases[0] =
163 p, sf, quad, use_pce_quad_points, normalize));
164 st_bases[1] =
166 p, ef, quad, use_pce_quad_points, normalize));
167 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> >
168 st_basis =
169 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(st_bases));
170 //std::cout << *st_basis << std::endl;
171
172 Stokhos::OrthogPolyApprox<int,double> u_st(st_basis), v_st(st_basis),
173 w_st(st_basis);
174 u_st.term(0, 0) = u.mean();
175 u_st.term(0, 1) = 1.0;
176 v_st.term(0, 0) = v.mean();
177 v_st.term(1, 1) = 1.0;
178
179 // Triple product tensor
180 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > st_Cijk =
181 st_basis->computeTripleProductTensor();
182
183 // Tensor product quadrature
184 Teuchos::RCP<const Stokhos::Quadrature<int,double> > st_quad;
185 if (!use_pce_quad_points) {
186#ifdef HAVE_STOKHOS_DAKOTA
187 if (sparse_grid)
188 st_quad = Teuchos::rcp(new Stokhos::SparseGridQuadrature<int,double>(st_basis, p));
189#endif
190 if (!sparse_grid)
191 st_quad = Teuchos::rcp(new Stokhos::TensorProductQuadrature<int,double>(st_basis));
192 }
193 else {
194 Teuchos::Array<double> st_points_0;
195 Teuchos::Array<double> st_weights_0;
196 Teuchos::Array< Teuchos::Array<double> > st_values_0;
197 st_bases[0]->getQuadPoints(p+1, st_points_0, st_weights_0, st_values_0);
198 Teuchos::Array<double> st_points_1;
199 Teuchos::Array<double> st_weights_1;
200 Teuchos::Array< Teuchos::Array<double> > st_values_1;
201 st_bases[1]->getQuadPoints(p+1, st_points_1, st_weights_1, st_values_1);
202 Teuchos::RCP< Teuchos::Array< Teuchos::Array<double> > > st_points =
203 Teuchos::rcp(new Teuchos::Array< Teuchos::Array<double> >(st_points_0.size()));
204 for (int i=0; i<st_points_0.size(); i++) {
205 (*st_points)[i].resize(2);
206 (*st_points)[i][0] = st_points_0[i];
207 (*st_points)[i][1] = st_points_1[i];
208 }
209 Teuchos::RCP< Teuchos::Array<double> > st_weights =
210 Teuchos::rcp(new Teuchos::Array<double>(st_weights_0));
211 Teuchos::RCP< const Stokhos::OrthogPolyBasis<int,double> > st_b =
212 st_basis;
213 st_quad =
214 Teuchos::rcp(new Stokhos::UserDefinedQuadrature<int,double>(st_b,
215 st_points,
216 st_weights));
217 }
218
219 // Quadrature expansion
221 st_Cijk,
222 st_quad);
223
224 // Compute w_st = u_st*v_st in Stieltjes basis
225 st_quad_exp.divide(w_st, u_st, v_st);
226
227 // Project w_st back to original basis
228 pce_quad_func st_func(w_st, *st_basis);
229 quad_exp.binary_op(st_func, w2, u, v);
230
231 // std::cout.precision(12);
232 // std::cout << w;
233 // std::cout << w2;
234 // std::cout << w_st;
235 mean[n] = w.mean();
236 mean_st[n] = w2.mean();
237 std_dev[n] = w.standard_deviation();
238 std_dev_st[n] = w2.standard_deviation();
239 pt[n] = w.evaluate(eval_pt);
240 pt_st[n] = w2.evaluate(eval_pt);
241 n++;
242 }
243
244 n = 0;
245 int wi=10;
246 std::cout << "Statistical error:" << std::endl;
247 std::cout << "p "
248 << std::setw(wi) << "mean" << " "
249 << std::setw(wi) << "mean_st" << " "
250 << std::setw(wi) << "std_dev" << " "
251 << std::setw(wi) << "std_dev_st" << " "
252 << std::setw(wi) << "point" << " "
253 << std::setw(wi) << "point_st" << std::endl;
254 for (unsigned int p=pmin; p<pmax; p++) {
255 std::cout.precision(3);
256 std::cout.setf(std::ios::scientific);
257 std::cout << p << " "
258 << std::setw(wi) << rel_err(mean[n], mean[np-1]) << " "
259 << std::setw(wi) << rel_err(mean_st[n], mean[np-1]) << " "
260 << std::setw(wi) << rel_err(std_dev[n], std_dev[np-1]) << " "
261 << std::setw(wi) << rel_err(std_dev_st[n], std_dev[np-1])
262 << " "
263 << std::setw(wi) << rel_err(pt[n], pt_true) << " "
264 << std::setw(wi) << rel_err(pt_st[n], pt_true)
265 << std::endl;
266 n++;
267 }
268
269 }
270 catch (std::exception& e) {
271 std::cout << e.what() << std::endl;
272 }
273}
Legendre polynomial basis using Clenshaw-Curtis quadrature points.
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
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 functional map...
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::ClenshawCurtisLegendreBasis< int, double > basis_type
const Stokhos::OrthogPolyApprox< int, double > & pce
double operator()(const Teuchos::Array< double > &x) const
exp_func(const Stokhos::OrthogPolyApprox< int, double > &pce_)
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
sin_func(const Stokhos::OrthogPolyApprox< int, double > &pce_)
double operator()(const Teuchos::Array< double > &x) const
const Stokhos::OrthogPolyApprox< int, double > & pce