Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
gram_schmidt_example3.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"
48#include "Stokhos_Sacado.hpp"
50
51#include "Teuchos_CommandLineProcessor.hpp"
52#include "Teuchos_TabularOutputter.hpp"
53
54// quadrature methods
56static const int num_quadrature_method = 2;
58 TENSOR, SPARSE };
59static const char *quadrature_method_names[] = {
60 "Tensor", "Sparse" };
61
62// reduced basis methods
64static const int num_reduced_basis_method = 5;
67static const char *reduced_basis_method_names[] = {
68 "Lanczos", "Monomial-Proj-GS", "Monomial-Proj-GS2", "Monomial-GS", "Lanczos-GS" };
69
70// orthogonalization methods
72static const int num_orthogonalization_method = 8;
75static const char *orthogonalization_method_names[] = {
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" };
77
78// quadrature reduction methods
80static const int num_quad_reduction_method = 4;
83static const char *quad_reduction_method_names[] = {
84 "None", "Q Squared", "Q Squared2", "Q2" };
85
86// solver methods
88static const int num_solver_method = 7;
91static const char *solver_method_names[] = {
92 "TRSM", "GLPK", "Clp", "Clp-IP", "qpOASES", "Basis Pursuit", "Orthogonal Matching Pursuit" };
93
96
97template <class ScalarType>
98inline ScalarType f(const Teuchos::Array<ScalarType>& x,
99 double a, double b) {
100 ScalarType y = a;
101 int n = x.size();
102 for (int i=0; i<n; i++)
103 y += x[i] / (i+1.0);
104 y = b + 1.0/y;
105 return y;
106}
107
108template <class ScalarType>
109inline ScalarType g(const Teuchos::Array<ScalarType>& x,
110 const ScalarType& y) {
111 ScalarType z = y;
112 //ScalarType z = 0.0;
113 int n = x.size();
114 for (int i=0; i<n; i++)
115 z += x[i];
116 z = std::exp(z);
117 //z = y*z;
118 return z;
119}
120
121double coeff_error(const pce_type& z, const pce_type& z2) {
122 double err_z = 0.0;
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;
127 }
128 return err_z;
129}
130
133 const Teuchos::Array<double>& weights = quad.getQuadWeights();
134 const Teuchos::Array< Teuchos::Array<double> >& vals =
136 int nqp = quad.size();
137 int npc = basis.size();
138 double max_err = 0.0;
139
140 // Loop over all basis function combinations
141 for (int i=0; i<npc; ++i) {
142 for (int j=0; j<npc; ++j) {
143
144 // Compute inner product of ith and jth basis function
145 double err = 0.0;
146 for (int k=0; k<nqp; ++k)
147 err += weights[k]*vals[k][i]*vals[k][j];
148
149 // Subtract off what it should be
150 if (i == j)
151 err -= basis.norm_squared(i);
152
153 // Accumulate max error
154 if (std::abs(err) > max_err)
155 max_err = std::abs(err);
156 }
157 }
158 return max_err;
159}
160
161struct MyOptions {
162 int d;
163 int d2;
165 int p_end;
167 double pole;
168 double shift;
174 int level;
181 bool use_Q;
186};
187
189public:
203};
204
205void compute_pces(bool compute_z_red, int p, const MyOptions& options,
206 MyResults& results)
207{
208 // Create product basis
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 =
213 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));
214
215 results.basis_size = basis->size();
216
217 // Quadrature
218 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad;
219 if (options.quad_method == TENSOR)
220 quad =
221 Teuchos::rcp(new Stokhos::TensorProductQuadrature<int,double>(basis));
222 else if (options.quad_method == SPARSE) {
223#ifdef HAVE_STOKHOS_DAKOTA
224 if (options.level == -1)
225 quad =
226 Teuchos::rcp(new Stokhos::SparseGridQuadrature<int,double>(basis));
227 else
228 quad =
229 Teuchos::rcp(new Stokhos::SparseGridQuadrature<int,double>(
230 basis, p, 1e-12, Pecos::SLOW_RESTRICTED_GROWTH));
231#else
232 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Sparse grid quadrature only supported when compiled with Dakota!");
233#endif
234 }
235
236 results.quad_size = quad->size();
237
238 // Triple product tensor
239 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
240 basis->computeTripleProductTensor();
241
242 // Quadrature expansion
243 Teuchos::RCP<Stokhos::QuadOrthogPolyExpansion<int,double> > quad_exp =
244 Teuchos::rcp(new Stokhos::QuadOrthogPolyExpansion<int,double>(basis,
245 Cijk,
246 quad));
247
248 // Create approximation
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++) {
254 x[i].copyForWrite();
255 x[i].reset(quad_exp);
256 x[i].term(i,1) = 1.0 / basis_vals[i+1];
257 }
258 Teuchos::Array<pce_type> x2(options.d2);
259 for (int i=0; i<options.d2; i++) {
260 x2[i] = x[i];
261 }
262
263 // Compute PCE via quadrature expansion
264 pce_type y = f(x, options.pole, options.shift);
265 results.z = g(x2, y);
266
267 if (!compute_z_red)
268 return;
269
270 // Create new basis from (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");
278 else if (options.reduced_basis_method == MONOMIAL_PROJ_GS2)
279 params.set("Reduced Basis Method", "Monomial Proj Gram-Schmidt2");
280 else if (options.reduced_basis_method == MONOMIAL_GS)
281 params.set("Reduced Basis Method", "Monomial Gram-Schmidt");
282 else if (options.reduced_basis_method == LANCZOS)
283 params.set("Reduced Basis Method", "Product Lanczos");
284 else if (options.reduced_basis_method == LANCZOS_GS)
285 params.set("Reduced Basis Method", "Product Lanczos Gram-Schmidt");
286 params.set("Verbose", options.verbose);
287 params.set("Project", options.project);
288 //params.set("Normalize", false);
289 params.set("Use Old Stieltjes Method", options.use_stieltjes);
290 params.set("Orthogonalization Method",
292 params.set("Rank Threshold", options.rank_threshold);
293 Teuchos::ParameterList& red_quad_params =
294 params.sublist("Reduced Quadrature");
295 red_quad_params.set(
296 "Reduced Quadrature Method",
298 red_quad_params.set(
299 "Solver Method", solver_method_names[options.solver_method]);
300 red_quad_params.set(
301 "Eliminate Dependent Rows", options.eliminate_dependent_rows);
302 red_quad_params.set("Write MPS File", false);
303 red_quad_params.set("Reduction Tolerance", options.reduction_tolerance);
304 red_quad_params.set("Verbose", options.verbose);
305 red_quad_params.set("Objective Value", options.objective_value);
306 red_quad_params.set("Q2 Rank Threshold", options.rank_threshold2);
307 red_quad_params.set(
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 =
315 factory.createReducedBasis(p, pces, quad, Cijk);
316 Teuchos::RCP<const Stokhos::Quadrature<int,double> > gs_quad =
317 gs_basis->getReducedQuadrature();
318
319 results.reduced_basis_size = gs_basis->size();
320 results.reduced_quad_size = gs_quad->size();
321
322 // Triple product tensor & expansion
323 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > gs_Cijk;
324 Teuchos::RCP< Teuchos::ParameterList > gs_exp_params =
325 Teuchos::rcp(new Teuchos::ParameterList);
326 if (options.reduced_basis_method == LANCZOS)
327 gs_Cijk = gs_basis->computeTripleProductTensor();
328 else {
329 gs_Cijk = Teuchos::null;
330 gs_exp_params->set("Use Quadrature for Times", true);
331 }
332 Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion<int,double> > gs_quad_exp =
334 gs_basis, gs_Cijk, gs_quad, gs_exp_params));
335
336 // Create new expansions
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());
342 }
343 pce_type y_gs(gs_quad_exp);
344 gs_basis->transformFromOriginalBasis(y.coeff(), y_gs.coeff());
345
346 // Compute z_gs = g(x2_gs, y_gs) in Gram-Schmidt basis
347 pce_type z_gs = g(x2_gs, y_gs);
348
349 // Project z_gs back to original basis
350 pce_type z2(quad_exp);
351 gs_basis->transformToOriginalBasis(z_gs.coeff(), z2.coeff());
352 results.z_red = z2;
353
354 // Compute discrete orthogonality error
355 results.disc_orthog_error = disc_orthog_error(*gs_basis, *gs_quad);
356}
357
358int main(int argc, char **argv)
359{
360 try {
361
362 // Setup command line options
363 Teuchos::CommandLineProcessor CLP;
364 CLP.setDocString(
365 "This example runs a Gram-Schmidt-based dimension reduction example.\n");
366 MyOptions options;
367
368 options.d = 4;
369 CLP.setOption("d", &options.d, "Stochastic dimension");
370
371 options.d2 = 1;
372 CLP.setOption("d2", &options.d2, "Intermediate stochastic dimension");
373
374 options.p_begin = 1;
375 CLP.setOption("p_begin", &options.p_begin, "Starting polynomial order");
376
377 options.p_end = 5;
378 CLP.setOption("p_end", &options.p_end, "Ending polynomial order");
379
380 options.p_truth = 7;
381 CLP.setOption("p_truth", &options.p_truth, "Truth polynomial order");
382
383 options.pole = 10.0;
384 CLP.setOption("pole", &options.pole, "Pole location");
385
386 options.shift = 0.0;
387 CLP.setOption("shift", &options.shift, "Shift location");
388
389 options.rank_threshold = 1.0e-120;
390 CLP.setOption("rank_threshold", &options.rank_threshold, "Rank threshold");
391
392 options.rank_threshold2 = 1.0e-120;
393 CLP.setOption("rank_threshold2", &options.rank_threshold2, "Rank threshold for Q2");
394
395 options.reduction_tolerance = 1.0e-12;
396 CLP.setOption("reduction_tolerance", &options.reduction_tolerance, "Quadrature reduction tolerance");
397
398 options.verbose = false;
399 CLP.setOption("verbose", "quiet", &options.verbose, "Verbose output");
400
401 options.quad_method = TENSOR;
402 CLP.setOption("quadrature_method", &options.quad_method,
404 quadrature_method_names, "Quadrature method");
405
406 options.level = -1;
407 CLP.setOption("level", &options.level,
408 "Sparse grid level (set to -1 to use default)");
409
411 CLP.setOption("reduced_basis_method", &options.reduced_basis_method,
413 reduced_basis_method_names, "Reduced basis method");
414
416 CLP.setOption("orthogonalization_method",
417 &options.orthogonalization_method,
421 "Orthogonalization method");
422
424 CLP.setOption("reduced_quadrature_method", &options.quad_reduction_method,
426 quad_reduction_method_names, "Reduced quadrature method");
427
428 options.solver_method = TRSM;
429 CLP.setOption("solver_method", &options.solver_method, num_solver_method,
431 "Reduced quadrature solver method");
432
434 CLP.setOption("quad_orthogonalization_method",
439 "Quadrature Orthogonalization method");
440
441 options.eliminate_dependent_rows = true;
442 CLP.setOption("cpqr", "no-cpqr", &options.eliminate_dependent_rows,
443 "Eliminate dependent rows in quadrature constraint matrix");
444
445 options.use_Q = true;
446 CLP.setOption("use-Q", "no-use-Q", &options.use_Q, "Use Q in LP");
447
448 options.restrict_r = false;
449 CLP.setOption("restrict-rank", "no-restrict-rank", &options.restrict_r,
450 "Restrict rank in LP");
451
452 options.objective_value = 0.0;
453 CLP.setOption("objective_value", &options.objective_value,
454 "Value for LP objective function");
455
456 options.project = true;
457 CLP.setOption("project", "no-project", &options.project,
458 "Use Projected Lanczos Method");
459
460 options.use_stieltjes = false;
461 CLP.setOption("stieltjes", "no-stieltjes", &options.use_stieltjes,
462 "Use Old Stieltjes Method");
463
464 CLP.parse( argc, argv );
465
466 std::cout << "Summary of command line options:" << std::endl
467 << "\tquadrature_method = "
469 << std::endl
470 << "\tlevel = " << options.level
471 << std::endl
472 << "\treduced_basis_method = "
474 << std::endl
475 << "\torthogonalization_method = "
477 << std::endl
478 << "\tquadrature_reduction_method = "
480 << std::endl
481 << "\tsolver_method = "
482 << solver_method_names[options.solver_method] << std::endl
483 << "\tquad_orthogonalization_method = "
485 << std::endl
486 << "\tcpqr = " << options.eliminate_dependent_rows
487 << std::endl
488 << "\tuse-Q = " << options.use_Q << std::endl
489 << "\trestrict-rank = " << options.restrict_r << std::endl
490 << "\tobjective_value = " << options.objective_value << std::endl
491 << "\tproject = " << options.project << std::endl
492 << "\tstieljtes = " << options.use_stieltjes << 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
500 << "\trank_threshold = " << options.rank_threshold << std::endl
501 << "\trank_threshold2 = " << options.rank_threshold2 << std::endl
502 << "\treduction_tolerance = " << options.reduction_tolerance
503 << std::endl
504 << "\tverbose = " << options.verbose << std::endl
505 << std::endl << std::endl;
506
507 std::stringstream ss;
508 typedef Teuchos::TabularOutputter TO;
509 TO out(ss);
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);
519 out.outputHeader();
520
521 MyResults results_truth;
522 compute_pces(false, options.p_truth, options, results_truth);
523
524 int n = options.p_end - options.p_begin + 1;
525 Teuchos::Array<MyResults> results(n);
526 for (int i=0; i<n; ++i) {
527 int p = options.p_begin + i;
528 compute_pces(true, p, options, results[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());
534 results[i].coeff_error = coeff_error(results_truth.z,
535 results[i].z);
536
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,
543 results[i].z_red);
544
545 out.outputField(p);
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);
550 out.outputField(results[i].coeff_error);
551 out.outputField(results[i].reduced_coeff_error);
552 out.outputField(results[i].disc_orthog_error);
553 out.nextRow();
554 }
555 std::cout << std::endl << ss.str() << std::endl;
556
557 }
558 catch (std::exception& e) {
559 std::cout << e.what() << std::endl;
560 }
561}
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[]
Orthogonalization_Method
@ HOUSEHOLDER_NP
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
Reduced_Basis_Method
@ MONOMIAL_PROJ_GS2
@ MONOMIAL_PROJ_GS
Quadrature_Reduction_Method
static const int num_orthogonalization_method
Quadrature_Method quad_method
Solver_Method solver_method
Orthogonalization_Method quad_orthogonalization_method
Quadrature_Reduction_Method quad_reduction_method
Reduced_Basis_Method reduced_basis_method
Orthogonalization_Method orthogonalization_method