Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
stieltjes_example4.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
55template <int Num>
57 static const int N = Num;
58 double a;
59
60 s_quad_func(double a_) : a(a_) {}
61
62 double operator() (double x[Num]) const {
63 double prod = 1;
64 for (int i=0; i<Num; i++)
65 prod *= x[i];
66 return 1.0 / (prod - a);
67 }
68};
69
70struct pce_quad_func {
74 pce(pce_), basis(basis_), vec(1) {}
75
76 double operator() (const double& a) const {
77 vec[0] = a;
78 return pce.evaluate(vec);
79 }
82 mutable Teuchos::Array<double> vec;
83};
84
85double rel_err(double a, double b) {
86 return std::abs(a-b)/std::abs(b);
87}
88
89int main(int argc, char **argv)
90{
91 try {
92
93 const int d = 5;
94 const int pmin = 1;
95 const int pmax = 10;
96 const int np = pmax-pmin+1;
97 const double a = 1.5;
98 bool use_pce_quad_points = false;
99 bool normalize = false;
100 bool project_integrals = false;
101 bool lanczos = true;
102 bool sparse_grid = true;
103#ifndef HAVE_STOKHOS_DAKOTA
104 sparse_grid = false;
105#endif
106 Teuchos::Array<double> mean(np), mean_st(np), std_dev(np), std_dev_st(np);
107 Teuchos::Array<double> pt(np), pt_st(np);
108
109 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
110
111 Teuchos::Array<double> eval_pt(d, 0.56789);
112 double pt_true;
113
114 // Loop over orders
115 int n = 0;
116 for (int p=pmin; p<=pmax; p++) {
117
118 std::cout << "p = " << p << std::endl;
119
120 // Create product basis
121 for (int i=0; i<d; i++)
122 bases[i] = Teuchos::rcp(new basis_type(p));
123 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
124 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));
125
126 // Create approximation
127 Stokhos::OrthogPolyApprox<int,double> s(basis), t1(basis), t2(basis);
128 Teuchos::Array< Stokhos::OrthogPolyApprox<int,double>* > xi(d);
129 for (int i=0; i<d; i++) {
130 xi[i] = new Stokhos::OrthogPolyApprox<int,double>(basis);
131 xi[i]->term(i, 1) = 1.0;
132 }
133
134 // Quadrature
135 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad;
136#ifdef HAVE_STOKHOS_DAKOTA
137 if (sparse_grid)
138 quad =
139 Teuchos::rcp(new Stokhos::SparseGridQuadrature<int,double>(basis, p));
140#endif
141 if (!sparse_grid)
142 quad =
143 Teuchos::rcp(new Stokhos::TensorProductQuadrature<int,double>(basis));
144
145 // Triple product tensor
146 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
147 basis->computeTripleProductTensor();
148
149 // Quadrature expansion
150 Stokhos::QuadOrthogPolyExpansion<int,double> quad_exp(basis, Cijk, quad);
151
152 // Compute PCE via quadrature expansion
153 s_quad_func<d> s_func(a);
156 for (int i=0; i<d; i++)
157 xip[i] = xi[i];
158 quad_exp.nary_op(s_func,s,xip);
159 quad_exp.divide(t1,1.0,s);
160 delete [] xip;
161
162 // compute true point
163 Teuchos::Array<double> xx(d);
164 for (int i=0; i<d; i++)
165 xx[i] = xi[i]->evaluate(eval_pt);
166 pt_true = s_func(&(xx[0]));
167 pt_true = 1.0/pt_true;
168
169 // Compute Stieltjes basis
170 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > st_bases(1);
171 Teuchos::RCP< Stokhos::LanczosProjPCEBasis<int,double> > stp_basis_s;
172 Teuchos::RCP< Stokhos::LanczosPCEBasis<int,double> > st_basis_s;
173 if (lanczos) {
174 if (project_integrals) {
175 stp_basis_s =
177 p, Teuchos::rcp(&s,false), Cijk, normalize, true));
178 st_bases[0] = stp_basis_s;
179 }
180 else {
181 st_basis_s =
183 p, Teuchos::rcp(&s,false), quad, normalize, true));
184 st_bases[0] = st_basis_s;
185 }
186 }
187 else {
188 st_bases[0] =
190 p, Teuchos::rcp(&s,false), quad, use_pce_quad_points,
191 normalize, project_integrals, Cijk));
192 }
193
194 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> >
195 st_basis =
196 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(st_bases));
197 //std::cout << *st_basis << std::endl;
198
199 Stokhos::OrthogPolyApprox<int,double> s_st(st_basis), t_st(st_basis);
200 if (lanczos) {
201 if (project_integrals) {
202 s_st.term(0, 0) = stp_basis_s->getNewCoeffs(0);
203 s_st.term(0, 1) = stp_basis_s->getNewCoeffs(1);
204 }
205 else {
206 s_st.term(0, 0) = st_basis_s->getNewCoeffs(0);
207 s_st.term(0, 1) = st_basis_s->getNewCoeffs(1);
208 }
209 }
210 else {
211 s_st.term(0, 0) = s.mean();
212 s_st.term(0, 1) = 1.0;
213 }
214
215 // Triple product tensor
216 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > st_Cijk =
217 st_basis->computeTripleProductTensor();
218
219 // Tensor product quadrature
220 Teuchos::RCP<const Stokhos::Quadrature<int,double> > st_quad;
221#ifdef HAVE_STOKHOS_DAKOTA
222 if (sparse_grid)
223 st_quad = Teuchos::rcp(new Stokhos::SparseGridQuadrature<int,double>(st_basis, p));
224#endif
225 if (!sparse_grid)
226 st_quad = Teuchos::rcp(new Stokhos::TensorProductQuadrature<int,double>(st_basis));
227
228 // Quadrature expansion
230 st_Cijk,
231 st_quad);
232
233 // Compute t_st = 1/s_st in Stieltjes basis
234 st_quad_exp.divide(t_st, 1.0, s_st);
235
236 // Project t_st back to original basis
237 pce_quad_func st_func(t_st, *st_basis);
238 quad_exp.unary_op(st_func, t2, s);
239
240 // std::cout.precision(12);
241 // std::cout << w;
242 // std::cout << w2;
243 // std::cout << w_st;
244 mean[n] = t1.mean();
245 mean_st[n] = t2.mean();
246 std_dev[n] = t1.standard_deviation();
247 std_dev_st[n] = t2.standard_deviation();
248 pt[n] = t1.evaluate(eval_pt);
249 pt_st[n] = t2.evaluate(eval_pt);
250 n++;
251
252 for (int i=0; i<d; i++)
253 delete xi[i];
254 }
255
256 n = 0;
257 int wi=10;
258 std::cout << "Statistical error:" << std::endl;
259 std::cout << "p "
260 << std::setw(wi) << "mean" << " "
261 << std::setw(wi) << "mean_st" << " "
262 << std::setw(wi) << "std_dev" << " "
263 << std::setw(wi) << "std_dev_st" << " "
264 << std::setw(wi) << "point" << " "
265 << std::setw(wi) << "point_st" << std::endl;
266 for (int p=pmin; p<pmax; p++) {
267 std::cout.precision(3);
268 std::cout.setf(std::ios::scientific);
269 std::cout << p << " "
270 << std::setw(wi) << rel_err(mean[n], mean[np-1]) << " "
271 << std::setw(wi) << rel_err(mean_st[n], mean[np-1]) << " "
272 << std::setw(wi) << rel_err(std_dev[n], std_dev[np-1]) << " "
273 << std::setw(wi) << rel_err(std_dev_st[n], std_dev[np-1])
274 << " "
275 << std::setw(wi) << rel_err(pt[n], pt_true) << " "
276 << std::setw(wi) << rel_err(pt_st[n], pt_true)
277 << std::endl;
278 n++;
279 }
280
281 }
282 catch (std::exception& e) {
283 std::cout << e.what() << std::endl;
284 }
285}
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.
Orthogonal polynomial expansions based on numerical quadrature.
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 nary_op(const FuncT &func, OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > **a)
void unary_op(const FuncT &func, OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Nonlinear unary function.
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
s_quad_func(double a_)
double operator()(double x[Num]) const
static const int N