Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
stieltjes_example5.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"
49
51#include "Stokhos_SDMUtils.hpp"
52
53#include "Teuchos_CommandLineProcessor.hpp"
54
55template <typename ordinal_type, typename scalar_type>
56void
57f_func(const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
58 const Teuchos::SerialDenseVector<ordinal_type,scalar_type>& x,
59 double shift,
60 Teuchos::SerialDenseVector<ordinal_type,scalar_type>& f)
61{
62 ordinal_type n = A.numCols();
63 Teuchos::SerialDenseVector<ordinal_type,scalar_type> y(n), s(n);
64 s.putScalar(shift);
65 y.assign(x);
66 y -= s;
67 f.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, A, y, 0.0);
68}
69
70template <typename ordinal_type, typename scalar_type>
72g_func(const Teuchos::SerialDenseVector<ordinal_type,scalar_type>& f)
73{
74 return std::exp(-f.dot(f));
75 //return -f.dot(f);
76}
77
80typedef Teuchos::SerialDenseMatrix<int,pce_type> SDM;
81typedef Teuchos::SerialDenseVector<int,pce_type> SDV;
82typedef Teuchos::SerialDenseMatrix<int,double> SDM0;
83typedef Teuchos::SerialDenseVector<int,double> SDV0;
84
85// measure transformation approaches
87const int num_mt_method = 3;
90const char *mt_method_names[] = { "Stieltjes", "Lanczos", "Gram-Schmidt" };
91
92double rel_err(double a, double b) {
93 return std::abs(a-b)/std::abs(b);
94}
95
96int main(int argc, char **argv)
97{
98 try {
99
100 // Setup command line options
101 Teuchos::CommandLineProcessor CLP;
102 CLP.setDocString(
103 "This example runs a Stieltjes-based dimension reduction example.\n");
104
105 int n = 4;
106 CLP.setOption("n", &n, "Number of random variables");
107
108 int m = 2;
109 CLP.setOption("m", &m, "Number of intermediate variables");
110
111 int pmin = 1;
112 CLP.setOption("pmin", &pmin, "Starting expansion order");
113
114 int pmax = 7;
115 CLP.setOption("pmax", &pmax, "Final expansion order");
116
117 MT_METHOD mt_method = MT_LANCZOS;
118 CLP.setOption("mt_method", &mt_method,
120 "Measure transformation method");
121
122 bool normalize = false;
123 CLP.setOption("normalize", "unnormalize", &normalize,
124 "Normalize PC basis");
125
126 bool project_integrals = false;
127 CLP.setOption("project", "no_project", &project_integrals,
128 "Project integrals");
129
130#ifdef HAVE_STOKHOS_DAKOTA
131 bool sparse_grid = true;
132 CLP.setOption("sparse", "tensor", &sparse_grid,
133 "Use Sparse or Tensor grid");
134#else
135 bool sparse_grid = false;
136#endif
137
138 double rank_threshold = 1.0e-12;
139 CLP.setOption("rank_threshold", &rank_threshold, "Rank threshold");
140
141 double reduction_tolerance = 1.0e-12;
142 CLP.setOption("reduction_tolerance", &reduction_tolerance, "Quadrature reduction tolerance");
143
144 bool verbose = false;
145 CLP.setOption("verbose", "quient", &verbose,
146 "Verbose output");
147
148 CLP.parse( argc, argv );
149
150 std::cout << "Summary of command line options:" << std::endl
151 << "\tm = " << m << std::endl
152 << "\tn = " << n << std::endl
153 << "\tmt_method = " << mt_method_names[mt_method] << std::endl
154 << "\tnormalize = " << normalize << std::endl
155 << "\tproject = " << project_integrals << std::endl
156 << "\tsparse = " << sparse_grid << std::endl
157 << "\trank_threshold = " << rank_threshold << std::endl
158 << "\treduction_tolerance = " << reduction_tolerance << std::endl;
159
160 int np = pmax-pmin+1;
161
162 Teuchos::Array<double> mean(np), mean_st(np), std_dev(np), std_dev_st(np);
163 Teuchos::Array<double> pt(np), pt_st(np);
164
165 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(n);
166
167 Teuchos::Array<double> eval_pt(n, 0.56789);
168 double pt_true = 0.0;
169
170 SDM0 B0(n,n), Q0, R0;
171 B0.random();
172 Teuchos::Array<double> w(n, 1.0);
173 Stokhos::QR_MGS(n, B0, w, Q0, R0);
174 SDM0 A0(Teuchos::View, B0, m, n);
175
176 // Loop over orders
177 int idx = 0;
178 for (int p=pmin; p<=pmax; p++) {
179
180 std::cout << "p = " << p;
181
182 // Create product basis
183 for (int i=0; i<n; i++)
184 bases[i] = Teuchos::rcp(new basis_type(p, normalize));
185 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
186 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));
187 std::cout << ", basis sz = " << basis->size();
188
189 // Quadrature
190 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad;
191#ifdef HAVE_STOKHOS_DAKOTA
192 if (sparse_grid)
193 quad =
194 Teuchos::rcp(new Stokhos::SparseGridQuadrature<int,double>(basis, p));
195#endif
196 if (!sparse_grid)
197 quad =
198 Teuchos::rcp(new Stokhos::TensorProductQuadrature<int,double>(basis));
199 std::cout << ", quad sz = " << quad->size();
200
201 // Triple product tensor
202 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
203 basis->computeTripleProductTensor();
204
205 // Quadrature expansion
206 Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion<int,double> > quad_exp =
207 Teuchos::rcp(new Stokhos::QuadOrthogPolyExpansion<int,double>(basis, Cijk, quad));
208
209 // Create approximation
210 SDV x(n);
211 for (int i=0; i<n; i++) {
212 x[i].copyForWrite();
213 x[i].reset(quad_exp);
214 if (normalize)
215 x[i].term(i, 1) = 1.0/std::sqrt(3);
216 else
217 x[i].term(i, 1) = 1.0;
218 }
219
220 // Create matrix
221 SDM A(m,n);
222 for (int i=0; i<m; i++)
223 for (int j=0; j<n; j++)
224 A(i,j) = A0(i,j);
225
226 // Compute f = A*(x-1)
227 SDV f(m);
228 f_func(A, x, 1.0, f);
229
230 //std::cout << "f = " << f << std::endl;
231
232 // Compute g = exp(-f^T*f)
233 pce_type g = g_func(f);
234
235 // compute true point
236 SDV0 x0(n), f0(m);
237 for (int i=0; i<n; i++)
238 x0[i] = x[i].evaluate(eval_pt);
239 f_func(A0, x0, 1.0, f0);
240 pt_true = g_func(f0);
241
242 // Compute reduced basis
243 Teuchos::ParameterList params;
244 params.set("Verbose", verbose);
245 if (mt_method == MT_GRAM_SCHMIDT)
246 params.set("Reduced Basis Method", "Monomial Proj Gram-Schmidt");
247 else if (mt_method == MT_LANCZOS)
248 params.set("Reduced Basis Method", "Product Lanczos");
249 else if (mt_method == MT_STIELTJES) {
250 params.set("Reduced Basis Method", "Product Lanczos");
251 params.set("Use Old Stieltjes Method", true);
252 }
253 params.set("Project", project_integrals);
254 params.set("Normalize", normalize);
255 params.set("Rank Threshold", rank_threshold);
256 Teuchos::ParameterList& red_quad_params =
257 params.sublist("Reduced Quadrature");
258 red_quad_params.set("Reduction Tolerance", reduction_tolerance);
259 red_quad_params.set("Verbose", verbose);
260 Teuchos::Array< Stokhos::OrthogPolyApprox<int,double> > pces(m);
261 for (int i=0; i<m; i++)
262 pces[i] = f[i].getOrthogPolyApprox();
264 Teuchos::RCP< Stokhos::ReducedPCEBasis<int,double> > gs_basis =
265 factory.createReducedBasis(p, pces, quad, Cijk);
266 Teuchos::RCP<const Stokhos::Quadrature<int,double> > gs_quad =
267 gs_basis->getReducedQuadrature();
268 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > gs_Cijk =
269 Teuchos::null;
270 Teuchos::RCP< Teuchos::ParameterList > gs_exp_params =
271 Teuchos::rcp(new Teuchos::ParameterList);
272 gs_exp_params->set("Use Quadrature for Times", true);
273 Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion<int,double> > st_quad_exp =
275 gs_basis,
276 gs_Cijk,
277 gs_quad,
278 gs_exp_params));
279 std::cout << ", red. basis sz = " <<gs_basis->size();
280 std::cout << ", red. quad sz = " << gs_quad->size();
281
282 SDV f_st(m);
283 for (int i=0; i<m; i++) {
284 f_st(i).copyForWrite();
285 f_st(i).reset(st_quad_exp);
286 gs_basis->transformFromOriginalBasis(f(i).coeff(), f_st(i).coeff());
287 }
288
289 // Compute g_st = exp(-f_st^T*f_st)
290 pce_type g_st = g_func(f_st);
291
292 // Project g_st back to original basis
293 pce_type g2(quad_exp);
294 gs_basis->transformToOriginalBasis(g_st.coeff(), g2.coeff());
295
296 // std::cout.precision(12);
297 // std::cout << g;
298 // std::cout << g2;
299 // std::cout << g_st;
300 mean[idx] = g.mean();
301 mean_st[idx] = g2.mean();
302 std_dev[idx] = g.standard_deviation();
303 std_dev_st[idx] = g2.standard_deviation();
304 pt[idx] = g.evaluate(eval_pt);
305 pt_st[idx] = g2.evaluate(eval_pt);
306 idx++;
307
308 std::cout << std::endl;
309 }
310
311 idx = 0;
312 int wi=10;
313 std::cout << "Statistical error:" << std::endl;
314 std::cout << "p "
315 << std::setw(wi) << "mean" << " "
316 << std::setw(wi) << "mean_st" << " "
317 << std::setw(wi) << "std_dev" << " "
318 << std::setw(wi) << "std_dev_st" << " "
319 << std::setw(wi) << "point" << " "
320 << std::setw(wi) << "point_st" << std::endl;
321 for (int p=pmin; p<pmax; p++) {
322 std::cout.precision(3);
323 std::cout.setf(std::ios::scientific);
324 std::cout << p << " "
325 << std::setw(wi) << rel_err(mean[idx], mean[np-1]) << " "
326 << std::setw(wi) << rel_err(mean_st[idx], mean[np-1]) << " "
327 << std::setw(wi) << rel_err(std_dev[idx], std_dev[np-1]) << " "
328 << std::setw(wi) << rel_err(std_dev_st[idx], std_dev[np-1])
329 << " "
330 << std::setw(wi) << rel_err(pt[idx], pt_true) << " "
331 << std::setw(wi) << rel_err(pt_st[idx], pt_true)
332 << std::endl;
333 idx++;
334 }
335
336 }
337 catch (std::exception& e) {
338 std::cout << e.what() << std::endl;
339 }
340}
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
Orthogonal polynomial expansions based on numerical quadrature.
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.
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
void QR_MGS(ordinal_type k, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, const Teuchos::Array< scalar_type > &w, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R)
Compute thin QR using modified Gram-Schmidt.
const char * mt_method_names[]
scalar_type g_func(const Teuchos::SerialDenseVector< ordinal_type, scalar_type > &f)
const MT_METHOD mt_method_values[]
Teuchos::SerialDenseMatrix< int, double > SDM0
int main(int argc, char **argv)
Teuchos::SerialDenseMatrix< int, pce_type > SDM
double rel_err(double a, double b)
Sacado::PCE::OrthogPoly< double, Stokhos::StandardStorage< int, double > > pce_type
@ MT_LANCZOS
@ MT_STIELTJES
@ MT_GRAM_SCHMIDT
Teuchos::SerialDenseVector< int, pce_type > SDV
Stokhos::LegendreBasis< int, double > basis_type
const int num_mt_method
Teuchos::SerialDenseVector< int, double > SDV0
void f_func(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, const Teuchos::SerialDenseVector< ordinal_type, scalar_type > &x, double shift, Teuchos::SerialDenseVector< ordinal_type, scalar_type > &f)