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;
87 bool sparse_grid =
true;
88#ifndef HAVE_STOKHOS_DAKOTA
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);
94 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
95 Teuchos::Array<double> eval_pt(d, 0.5);
100 for (
unsigned int p=pmin; p<=pmax; p++) {
102 std::cout <<
"p = " << p << std::endl;
105 for (
unsigned int i=0; i<d; i++)
107 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
113 for (
unsigned int i=0; i<d; i++) {
117 double x_pt = x.evaluate(eval_pt);
118 pt_true = std::sin(x_pt)/(1.0 +
exp(x_pt));
121 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad;
122#ifdef HAVE_STOKHOS_DAKOTA
125 Teuchos::rcp(
new Stokhos::SparseGridQuadrature<int,double>(basis, p));
132 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
133 basis->computeTripleProductTensor();
141 quad_exp.plusEqual(v, 1.0);
142 quad_exp.divide(w,u,v);
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;
149 if (project_integrals) {
152 p, Teuchos::rcp(&u,
false), Cijk, normalize,
true));
155 p, Teuchos::rcp(&v,
false), Cijk, normalize,
true));
156 st_bases[0] = stp_basis_u;
157 st_bases[1] = stp_basis_v;
162 p, Teuchos::rcp(&u,
false), quad, normalize,
true));
165 p, Teuchos::rcp(&v,
false), quad, normalize,
true));
166 st_bases[0] = st_basis_u;
167 st_bases[1] = st_basis_v;
173 p, Teuchos::rcp(&u,
false), quad, use_pce_quad_points,
174 normalize, project_integrals, Cijk));
177 p, Teuchos::rcp(&v,
false), quad, use_pce_quad_points,
178 normalize, project_integrals, Cijk));
181 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> >
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);
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);
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;
210 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > st_Cijk =
211 st_basis->computeTripleProductTensor();
214 Teuchos::RCP<const Stokhos::Quadrature<int,double> > st_quad;
215 if (!use_pce_quad_points) {
216#ifdef HAVE_STOKHOS_DAKOTA
218 st_quad = Teuchos::rcp(
new Stokhos::SparseGridQuadrature<int,double>(st_basis, p));
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];
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 =
255 st_quad_exp.divide(w_st, u_st, v_st);
259 quad_exp.binary_op(st_func, w2, u, v);
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);
276 std::cout <<
"Statistical error:" << std::endl;
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])
293 << std::setw(wi) <<
rel_err(pt[n], pt_true) <<
" "
294 << std::setw(wi) <<
rel_err(pt_st[n], pt_true)
300 catch (std::exception& e) {
301 std::cout << e.what() << std::endl;