51#include "Teuchos_CommandLineProcessor.hpp"
52#include "Teuchos_TabularOutputter.hpp"
68 "Lanczos",
"Monomial-Proj-GS",
"Monomial-Proj-GS2",
"Monomial-GS",
"Lanczos-GS" };
76 "Householder",
"Householder without Pivoting",
"Classical Gram-Schmidt",
"Modified Gram-Schmidt",
"Modified Gram-Schmidt with Reorthogonalization",
"Modified Gram-Schmidt without Pivoting",
"Modified Gram-Schmidt without Pivoting with Reorthogonalization",
"SVD" };
84 "None",
"Q Squared",
"Q Squared2",
"Q2" };
92 "TRSM",
"GLPK",
"Clp",
"Clp-IP",
"qpOASES",
"Basis Pursuit",
"Orthogonal Matching Pursuit" };
97template <
class ScalarType>
98inline ScalarType
f(
const Teuchos::Array<ScalarType>& x,
102 for (
int i=0; i<n; i++)
108template <
class ScalarType>
109inline ScalarType
g(
const Teuchos::Array<ScalarType>& x,
110 const ScalarType& y) {
114 for (
int i=0; i<n; i++)
123 int n = std::min(z.size(), z2.size());
124 for (
int i=0; i<n; i++) {
125 double ew = std::abs(z.coeff(i)-z2.coeff(i));
126 if (ew > err_z) err_z = ew;
134 const Teuchos::Array< Teuchos::Array<double> >& vals =
136 int nqp = quad.
size();
137 int npc = basis.
size();
138 double max_err = 0.0;
141 for (
int i=0; i<npc; ++i) {
142 for (
int j=0;
j<npc; ++
j) {
146 for (
int k=0; k<nqp; ++k)
147 err += weights[k]*vals[k][i]*vals[k][
j];
154 if (std::abs(err) > max_err)
155 max_err = std::abs(err);
209 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(options.
d);
210 for (
int i=0; i<options.
d; i++)
211 bases[i] = Teuchos::rcp(
new basis_type(p,
true));
212 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
218 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad;
223#ifdef HAVE_STOKHOS_DAKOTA
224 if (options.
level == -1)
226 Teuchos::rcp(
new Stokhos::SparseGridQuadrature<int,double>(basis));
229 Teuchos::rcp(
new Stokhos::SparseGridQuadrature<int,double>(
230 basis, p, 1e-12, Pecos::SLOW_RESTRICTED_GROWTH));
232 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Sparse grid quadrature only supported when compiled with Dakota!");
239 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
240 basis->computeTripleProductTensor();
243 Teuchos::RCP<Stokhos::QuadOrthogPolyExpansion<int,double> > quad_exp =
249 Teuchos::Array<double> point(options.
d, 1.0);
250 Teuchos::Array<double> basis_vals(basis->size());
251 basis->evaluateBases(point, basis_vals);
252 Teuchos::Array<pce_type> x(options.
d);
253 for (
int i=0; i<options.
d; i++) {
255 x[i].reset(quad_exp);
256 x[i].term(i,1) = 1.0 / basis_vals[i+1];
258 Teuchos::Array<pce_type> x2(options.
d2);
259 for (
int i=0; i<options.
d2; i++) {
265 results.
z =
g(x2, y);
271 Teuchos::Array< Stokhos::OrthogPolyApprox<int,double> > pces(options.
d2+1);
272 for (
int i=0; i<options.
d2; i++)
273 pces[i] = x2[i].getOrthogPolyApprox();
274 pces[options.
d2] = y.getOrthogPolyApprox();
275 Teuchos::ParameterList params;
277 params.set(
"Reduced Basis Method",
"Monomial Proj Gram-Schmidt");
279 params.set(
"Reduced Basis Method",
"Monomial Proj Gram-Schmidt2");
281 params.set(
"Reduced Basis Method",
"Monomial Gram-Schmidt");
283 params.set(
"Reduced Basis Method",
"Product Lanczos");
285 params.set(
"Reduced Basis Method",
"Product Lanczos Gram-Schmidt");
286 params.set(
"Verbose", options.
verbose);
287 params.set(
"Project", options.
project);
289 params.set(
"Use Old Stieltjes Method", options.
use_stieltjes);
290 params.set(
"Orthogonalization Method",
293 Teuchos::ParameterList& red_quad_params =
294 params.sublist(
"Reduced Quadrature");
296 "Reduced Quadrature Method",
302 red_quad_params.set(
"Write MPS File",
false);
304 red_quad_params.set(
"Verbose", options.
verbose);
308 "Orthogonalization Method",
310 red_quad_params.set(
"Use Q in LP", options.
use_Q);
311 red_quad_params.set(
"Restrict Rank", options.
restrict_r);
312 red_quad_params.set(
"Order Restriction", 2*p);
314 Teuchos::RCP< Stokhos::ReducedPCEBasis<int,double> > gs_basis =
316 Teuchos::RCP<const Stokhos::Quadrature<int,double> > gs_quad =
317 gs_basis->getReducedQuadrature();
323 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > gs_Cijk;
324 Teuchos::RCP< Teuchos::ParameterList > gs_exp_params =
325 Teuchos::rcp(
new Teuchos::ParameterList);
327 gs_Cijk = gs_basis->computeTripleProductTensor();
329 gs_Cijk = Teuchos::null;
330 gs_exp_params->set(
"Use Quadrature for Times",
true);
332 Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion<int,double> > gs_quad_exp =
334 gs_basis, gs_Cijk, gs_quad, gs_exp_params));
337 Teuchos::Array<pce_type> x2_gs(options.
d2);
338 for (
int i=0; i<options.
d2; i++) {
339 x2_gs[i].copyForWrite();
340 x2_gs[i].reset(gs_quad_exp);
341 gs_basis->transformFromOriginalBasis(x2[i].coeff(), x2_gs[i].coeff());
344 gs_basis->transformFromOriginalBasis(y.coeff(), y_gs.coeff());
351 gs_basis->transformToOriginalBasis(z_gs.coeff(), z2.coeff());
363 Teuchos::CommandLineProcessor
CLP;
365 "This example runs a Gram-Schmidt-based dimension reduction example.\n");
369 CLP.setOption(
"d", &options.
d,
"Stochastic dimension");
372 CLP.setOption(
"d2", &options.
d2,
"Intermediate stochastic dimension");
375 CLP.setOption(
"p_begin", &options.
p_begin,
"Starting polynomial order");
378 CLP.setOption(
"p_end", &options.
p_end,
"Ending polynomial order");
381 CLP.setOption(
"p_truth", &options.
p_truth,
"Truth polynomial order");
384 CLP.setOption(
"pole", &options.
pole,
"Pole location");
387 CLP.setOption(
"shift", &options.
shift,
"Shift location");
399 CLP.setOption(
"verbose",
"quiet", &options.
verbose,
"Verbose output");
407 CLP.setOption(
"level", &options.
level,
408 "Sparse grid level (set to -1 to use default)");
416 CLP.setOption(
"orthogonalization_method",
421 "Orthogonalization method");
431 "Reduced quadrature solver method");
434 CLP.setOption(
"quad_orthogonalization_method",
439 "Quadrature Orthogonalization method");
443 "Eliminate dependent rows in quadrature constraint matrix");
445 options.
use_Q =
true;
446 CLP.setOption(
"use-Q",
"no-use-Q", &options.
use_Q,
"Use Q in LP");
449 CLP.setOption(
"restrict-rank",
"no-restrict-rank", &options.
restrict_r,
450 "Restrict rank in LP");
454 "Value for LP objective function");
457 CLP.setOption(
"project",
"no-project", &options.
project,
458 "Use Projected Lanczos Method");
462 "Use Old Stieltjes Method");
466 std::cout <<
"Summary of command line options:" << std::endl
467 <<
"\tquadrature_method = "
470 <<
"\tlevel = " << options.
level
472 <<
"\treduced_basis_method = "
475 <<
"\torthogonalization_method = "
478 <<
"\tquadrature_reduction_method = "
481 <<
"\tsolver_method = "
483 <<
"\tquad_orthogonalization_method = "
488 <<
"\tuse-Q = " << options.
use_Q << std::endl
489 <<
"\trestrict-rank = " << options.
restrict_r << std::endl
491 <<
"\tproject = " << options.
project << std::endl
493 <<
"\tp_begin = " << options.
p_begin << std::endl
494 <<
"\tp_end = " << options.
p_end << std::endl
495 <<
"\tp_truth = " << options.
p_truth << std::endl
496 <<
"\td = " << options.
d << std::endl
497 <<
"\td2 = " << options.
d2 << std::endl
498 <<
"\tpole = " << options.
pole << std::endl
499 <<
"\tshift = " << options.
shift << std::endl
504 <<
"\tverbose = " << options.
verbose << std::endl
505 << std::endl << std::endl;
507 std::stringstream ss;
508 typedef Teuchos::TabularOutputter TO;
510 out.setFieldTypePrecision(TO::DOUBLE, 5);
511 out.pushFieldSpec(
"\"Order\"", TO::INT);
512 out.pushFieldSpec(
"\"Basis Size\"", TO::INT);
513 out.pushFieldSpec(
"\"Quad. Size\"", TO::INT);
514 out.pushFieldSpec(
"\"Red. Basis Size\"", TO::INT);
515 out.pushFieldSpec(
"\"Red. Quad. Size\"", TO::INT);
516 out.pushFieldSpec(
"\"Coeff. Error\"", TO::DOUBLE);
517 out.pushFieldSpec(
"\"Red. Coeff. Error\"", TO::DOUBLE);
518 out.pushFieldSpec(
"\"Disc. Orthog. Error\"", TO::DOUBLE);
525 Teuchos::Array<MyResults> results(n);
526 for (
int i=0; i<n; ++i) {
529 results[i].mean_error =
530 std::abs(results_truth.
z.mean()-results[i].z.mean()) /
531 std::abs(results_truth.
z.mean());
532 results[i].std_dev_error =
533 std::abs(results_truth.
z.standard_deviation()-results[i].z.standard_deviation()) / std::abs(results_truth.
z.standard_deviation());
537 results[i].reduced_mean_error =
538 std::abs(results_truth.
z.mean()-results[i].z_red.mean()) /
539 std::abs(results_truth.
z.mean());
540 results[i].reduced_std_dev_error =
541 std::abs(results_truth.
z.standard_deviation()-results[i].z_red.standard_deviation()) / std::abs(results_truth.
z.standard_deviation());
542 results[i].reduced_coeff_error =
coeff_error(results_truth.
z,
546 out.outputField(results[i].basis_size);
547 out.outputField(results[i].quad_size);
548 out.outputField(results[i].reduced_basis_size);
549 out.outputField(results[i].reduced_quad_size);
551 out.outputField(results[i].reduced_coeff_error);
555 std::cout << std::endl << ss.str() << std::endl;
558 catch (std::exception& e) {
559 std::cout << e.what() << std::endl;
double reduced_mean_error
double reduced_std_dev_error
double reduced_coeff_error
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
Abstract base class for multivariate orthogonal polynomials.
virtual const Teuchos::Array< value_type > & norm_squared() const =0
Return array storing norm-squared of each basis polynomial.
virtual ordinal_type size() const =0
Return total size of basis.
Orthogonal polynomial expansions based on numerical quadrature.
Abstract base class for quadrature methods.
virtual ordinal_type size() const =0
Get number of quadrature points.
virtual const Teuchos::Array< value_type > & getQuadWeights() const =0
Get quadrature weights.
virtual const Teuchos::Array< Teuchos::Array< value_type > > & getBasisAtQuadPoints() const =0
Get values of basis at quadrature points.
Generate a basis from a given set of PCE expansions that is orthogonal with respect to the product me...
virtual Teuchos::RCP< Stokhos::ReducedPCEBasis< ordinal_type, value_type > > createReducedBasis(ordinal_type p, const Teuchos::Array< Stokhos::OrthogPolyApprox< ordinal_type, value_type > > &pce, const Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > &quad, const Teuchos::RCP< const Stokhos::Sparse3Tensor< ordinal_type, value_type > > &Cijk) const
Get reduced quadrature object.
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
Sacado::ETPCE::OrthogPoly< double, Stokhos::StandardStorage< int, double > > pce_type
double disc_orthog_error(const Stokhos::OrthogPolyBasis< int, double > &basis, const Stokhos::Quadrature< int, double > &quad)
static const Solver_Method solver_method_values[]
static const char * solver_method_names[]
double coeff_error(const pce_type &z, const pce_type &z2)
static const Quadrature_Reduction_Method quad_reduction_method_values[]
static const int num_quadrature_method
int main(int argc, char **argv)
@ ORTHOGONAL_MATCHING_PURSUIT
static const char * reduced_basis_method_names[]
void compute_pces(bool compute_z_red, int p, const MyOptions &options, MyResults &results)
static const Orthogonalization_Method orthogonalization_method_values[]
static const char * orthogonalization_method_names[]
static const int num_solver_method
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
static const char * quadrature_method_names[]
static const Quadrature_Method quadrature_method_values[]
static const int num_quad_reduction_method
static const char * quad_reduction_method_names[]
Stokhos::LegendreBasis< int, double > basis_type
static const Reduced_Basis_Method reduced_basis_method_values[]
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
static const int num_reduced_basis_method
Quadrature_Reduction_Method
static const int num_orthogonalization_method
Quadrature_Method quad_method
Solver_Method solver_method
double reduction_tolerance
Orthogonalization_Method quad_orthogonalization_method
Quadrature_Reduction_Method quad_reduction_method
Reduced_Basis_Method reduced_basis_method
Orthogonalization_Method orthogonalization_method
bool eliminate_dependent_rows