66 mutable Teuchos::Array<double>
vec;
86 const unsigned int d = 2;
87 const unsigned int pmin = 2;
88 const unsigned int pmax = 10;
89 const unsigned int np = pmax-pmin+1;
90 bool use_pce_quad_points =
false;
91 bool normalize =
true;
92 bool sparse_grid =
true;
93 #ifndef HAVE_STOKHOS_DAKOTA 96 Teuchos::Array<double> mean(np), mean_st(np), std_dev(np), std_dev_st(np);
97 Teuchos::Array<double> pt(np), pt_st(np);
99 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
100 Teuchos::Array<double> eval_pt(d, 0.5);
105 for (
unsigned int p=pmin; p<=pmax; p++) {
107 std::cout <<
"p = " << p << std::endl;
110 for (
unsigned int i=0; i<d; i++)
112 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
118 for (
unsigned int i=0; i<d; i++) {
122 double x_pt =
x.evaluate(eval_pt);
126 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad;
127 #ifdef HAVE_STOKHOS_DAKOTA 130 Teuchos::rcp(
new Stokhos::SparseGridQuadrature<int,double>(basis, p));
137 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
138 basis->computeTripleProductTensor();
149 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > st_bases(1);
150 Teuchos::RCP<const Stokhos::LanczosProjPCEBasis<int,double> > st_1d_basis
152 p, Teuchos::rcp(&u,
false), Cijk, normalize));
153 st_bases[0] = st_1d_basis;
155 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> >
161 u_st.term(0, 0) = st_1d_basis->getNewCoeffs(0);
162 u_st.term(0, 1) = st_1d_basis->getNewCoeffs(1);
165 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > st_Cijk =
166 st_basis->computeTripleProductTensor();
169 Teuchos::RCP<const Stokhos::Quadrature<int,double> > st_quad;
170 if (!use_pce_quad_points) {
171 #ifdef HAVE_STOKHOS_DAKOTA 173 st_quad = Teuchos::rcp(
new Stokhos::SparseGridQuadrature<int,double>(st_basis, p));
179 Teuchos::Array<double> st_points_0;
180 Teuchos::Array<double> st_weights_0;
181 Teuchos::Array< Teuchos::Array<double> > st_values_0;
182 st_bases[0]->getQuadPoints(p+1, st_points_0, st_weights_0, st_values_0);
183 Teuchos::Array<double> st_points_1;
184 Teuchos::Array<double> st_weights_1;
185 Teuchos::Array< Teuchos::Array<double> > st_values_1;
186 st_bases[1]->getQuadPoints(p+1, st_points_1, st_weights_1, st_values_1);
187 Teuchos::RCP< Teuchos::Array< Teuchos::Array<double> > > st_points =
188 Teuchos::rcp(
new Teuchos::Array< Teuchos::Array<double> >(st_points_0.size()));
189 for (
int i=0; i<st_points_0.size(); i++) {
190 (*st_points)[i].resize(2);
191 (*st_points)[i][0] = st_points_0[i];
192 (*st_points)[i][1] = st_points_1[i];
194 Teuchos::RCP< Teuchos::Array<double> > st_weights =
195 Teuchos::rcp(
new Teuchos::Array<double>(st_weights_0));
196 Teuchos::RCP< const Stokhos::OrthogPolyBasis<int,double> > st_b =
210 st_quad_exp.
exp(w_st, u_st);
221 mean_st[n] = w2.
mean();
222 std_dev[n] = w.standard_deviation();
224 pt[n] = w.evaluate(eval_pt);
231 std::cout <<
"Statistical error:" << std::endl;
233 << std::setw(wi) <<
"mean" <<
" " 234 << std::setw(wi) <<
"mean_st" <<
" " 235 << std::setw(wi) <<
"std_dev" <<
" " 236 << std::setw(wi) <<
"std_dev_st" <<
" " 237 << std::setw(wi) <<
"point" <<
" " 238 << std::setw(wi) <<
"point_st" << std::endl;
239 for (
unsigned int p=pmin; p<pmax; p++) {
240 std::cout.precision(3);
241 std::cout.setf(std::ios::scientific);
242 std::cout << p <<
" " 243 << std::setw(wi) <<
rel_err(mean[n], mean[np-1]) <<
" " 244 << std::setw(wi) <<
rel_err(mean_st[n], mean[np-1]) <<
" " 245 << std::setw(wi) <<
rel_err(std_dev[n], std_dev[np-1]) <<
" " 246 << std::setw(wi) <<
rel_err(std_dev_st[n], std_dev[np-1])
248 << std::setw(wi) <<
rel_err(pt[n], pt_true) <<
" " 249 << std::setw(wi) <<
rel_err(pt_st[n], pt_true)
255 catch (std::exception& e) {
256 std::cout << e.what() << std::endl;
void sin(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
int main(int argc, char **argv)
void exp(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
const Stokhos::OrthogPolyApprox< int, double > & pce
sin_func(const Stokhos::OrthogPolyApprox< int, double > &pce_)
pce_quad_func(const Stokhos::OrthogPolyApprox< int, double > &pce_, const Stokhos::OrthogPolyBasis< int, double > &basis_)
Teuchos::Array< double > vec
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
value_type standard_deviation() const
Compute standard deviation of expansion.
const Stokhos::OrthogPolyBasis< int, double > & basis
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
const Stokhos::OrthogPolyApprox< int, double > & pce
Stokhos::ClenshawCurtisLegendreBasis< int, double > basis_type
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
double rel_err(double a, double b)
KOKKOS_INLINE_FUNCTION PCE< Storage > exp(const PCE< Storage > &a)
double operator()(const Teuchos::Array< double > &x) const
value_type mean() const
Compute mean of expansion.
Legendre polynomial basis using Clenshaw-Curtis quadrature points.
KOKKOS_INLINE_FUNCTION PCE< Storage > sin(const PCE< Storage > &a)
value_type evaluate(const Teuchos::Array< value_type > &point) const
Evaluate polynomial approximation at a point.
Orthogonal polynomial expansions based on numerical quadrature.
double operator()(const double &a, const double &b) const
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules...
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.