Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_MonomialProjGramSchmidtPCEBasis2Imp.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
47
48template <typename ordinal_type, typename value_type>
51 ordinal_type max_p,
52 const Teuchos::Array< Stokhos::OrthogPolyApprox<ordinal_type, value_type> >& pce,
53 const Teuchos::RCP<const Stokhos::Quadrature<ordinal_type, value_type> >& quad,
54 const Teuchos::ParameterList& params_) :
55 name("Monomial Proj Gram Schmidt PCE Basis"),
56 params(params_),
57 pce_basis(pce[0].basis()),
58 pce_sz(pce_basis->size()),
59 p(max_p),
60 d(pce.size()),
61 verbose(params.get("Verbose", false)),
62 rank_threshold(params.get("Rank Threshold", 1.0e-12)),
63 orthogonalization_method(params.get("Orthogonalization Method",
64 "Householder"))
65{
66 // Check for pce's that are constant and don't represent true random
67 // dimensions
68 Teuchos::Array< Stokhos::OrthogPolyApprox<ordinal_type, value_type> > non_const_pce;
69 for (ordinal_type i=0; i<pce.size(); i++) {
70 if (pce[i].standard_deviation() > 1.0e-15)
71 non_const_pce.push_back(pce[i]);
72 }
73 d = non_const_pce.size();
74
75 // Build Q, Qp matrices
76 SDM A, F;
77 sz = buildQ(max_p, rank_threshold, non_const_pce, quad, terms, num_terms,
78 Qp, A, F);
79 Q.reshape(A.numRows(), sz);
80 ordinal_type ret =
81 Q.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, A, Qp, 0.0);
82 TEUCHOS_ASSERT(ret == 0);
83
84//print_matlab(std::cout << "Qp = ", Qp);
85
86 // Compute reduced quadrature rule
87 Teuchos::ParameterList quad_params = params.sublist("Reduced Quadrature");
89 quad_params);
90 SDM Q2;
91 if (quad_params.isParameter("Reduced Quadrature Method") &&
92 quad_params.get<std::string>("Reduced Quadrature Method") == "Q2") {
93 Teuchos::Array< Stokhos::MultiIndex<ordinal_type> > terms2;
94 Teuchos::Array<ordinal_type> num_terms2;
95 value_type rank_threshold2 = quad_params.get("Q2 Rank Threshold",
97 SDM Qp2, A2, F2;
98
99 // Build basis, quadrature of order 2*max_p
100 Teuchos::Array< Stokhos::OrthogPolyApprox<ordinal_type, value_type> > pce2(non_const_pce);
101 Teuchos::RCP<const Stokhos::Quadrature<ordinal_type, value_type> > quad2 = quad;
102 Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type,value_type> > basis2 = pce_basis;
103 if (2*max_p > pce_basis->order()) {
104
105 // Basis
106 Teuchos::RCP<const Stokhos::ProductBasis<ordinal_type,value_type> > prod_basis = Teuchos::rcp_dynamic_cast<const Stokhos::ProductBasis<ordinal_type,value_type> >(pce_basis);
107 ordinal_type dim = prod_basis->dimension();
108 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> > > bases(dim);
109 for (ordinal_type i=0; i<dim; i++)
110 bases[i] = prod_basis->getCoordinateBases()[i]->cloneWithOrder(2*max_p);
111 Teuchos::RCP< const Stokhos::CompletePolynomialBasis<ordinal_type,value_type> > cp_basis2 = Teuchos::rcp(new Stokhos::CompletePolynomialBasis<ordinal_type,value_type>(bases));
112 basis2 = cp_basis2;
113 //quad_params.sublist("Basis").set("Stochastic Galerkin Basis", basis2);
114 std::cout << "built new basis of dimension " << basis2->dimension()
115 << " and order " << basis2->order()
116 << " with total size " << basis2->size() << std::endl;
117
118 // Quadrature
119 // quad_params.sublist("Quadrature").set("Quadrature Order", 2*max_p);
120 // quad2 =
121 // Stokhos::QuadratureFactory<ordinal_type,value_type>::create(quad_params);
122
124 std::cout << "built new quadrature with total size " << quad2->size()
125 << std::endl;
126
127 // Project pce to new basis
128 for (ordinal_type i=0; i<d; i++) {
129 pce2[i].reset(basis2); // this keeps lower order coeffs and sets
130 // higher order ones to 0
131 }
132 }
133
134 // Build Q matrix of order 2*max_p
135 ordinal_type sz2 =
136 buildQ(2*max_p, rank_threshold2, pce2, quad2, terms2, num_terms2,
137 Qp2, A2, F2);
138 //print_matlab(std::cout << "Qp2 = ", Qp2);
139
140 // Get quadrature data
141 const Teuchos::Array<value_type>& weights = quad->getQuadWeights();
142 const Teuchos::Array< Teuchos::Array<value_type> >& points =
143 quad->getQuadPoints();
144 ordinal_type nqp = weights.size();
145
146 // Original basis at quadrature points -- needed to transform expansions
147 // in this basis back to original
148 ordinal_type pce_sz2 = basis2->size();
149 SDM AA(nqp, pce_sz2);
150 Teuchos::Array<value_type> basis_vals(pce_sz2);
151 for (ordinal_type i=0; i<nqp; i++) {
152 basis2->evaluateBases(points[i], basis_vals);
153 for (ordinal_type j=0; j<pce_sz2; j++)
154 AA(i,j) = basis_vals[j];
155 }
156 Q2.reshape(nqp, sz2);
157 ret = Q2.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, AA, Qp2, 0.0);
158 TEUCHOS_ASSERT(ret == 0);
159 reduced_quad = quad_factory.createReducedQuadrature(Q, Q2, F, weights);
160
161 // // Get quadrature data
162 // const Teuchos::Array<value_type>& weights2 = quad2->getQuadWeights();
163 // ordinal_type nqp2 = weights2.size();
164 // Q2.reshape(nqp2, sz2);
165 // ret = Q2.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, A2, Qp2, 0.0);
166 // TEUCHOS_ASSERT(ret == 0);
167 // reduced_quad = quad_factory.createReducedQuadrature(Q, Q2, F2, weights2);
168 }
169 else {
170 const Teuchos::Array<value_type>& weights = quad->getQuadWeights();
171 reduced_quad = quad_factory.createReducedQuadrature(Q, Q2, F, weights);
172 }
173
174 // Basis is orthonormal by construction
175 norms.resize(sz, 1.0);
176}
177
178template <typename ordinal_type, typename value_type>
179ordinal_type
181buildQ(
182 ordinal_type max_p,
183 value_type threshold,
184 const Teuchos::Array< Stokhos::OrthogPolyApprox<ordinal_type, value_type> >& pce,
185 const Teuchos::RCP<const Stokhos::Quadrature<ordinal_type, value_type> >& quad,
186 Teuchos::Array< Stokhos::MultiIndex<ordinal_type> >& terms_,
187 Teuchos::Array<ordinal_type>& num_terms_,
188 SDM& Qp_, SDM& A_, SDM& F_)
189{
190 Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type,value_type> > pce_basis_ = pce[0].basis();
191 ordinal_type pce_sz_ = pce_basis_->size();
192
193 // Get quadrature data
194 const Teuchos::Array<value_type>& weights = quad->getQuadWeights();
195 const Teuchos::Array< Teuchos::Array<value_type> >& points =
196 quad->getQuadPoints();
197 const Teuchos::Array< Teuchos::Array<value_type> >& basis_values =
198 quad->getBasisAtQuadPoints();
199 ordinal_type nqp = weights.size();
200
201 // Original basis at quadrature points -- needed to transform expansions
202 // in this basis back to original
203 A_.reshape(nqp, pce_sz_);
204 for (ordinal_type i=0; i<nqp; i++)
205 for (ordinal_type j=0; j<pce_sz_; j++)
206 A_(i,j) = basis_values[i][j];
207
208 // Compute norms of each pce for rescaling
209 Teuchos::Array<value_type> pce_norms(d, 0.0);
210 for (ordinal_type j=0; j<d; j++) {
211 for (ordinal_type i=0; i<pce_sz_; i++)
212 pce_norms[j] += (pce[j])[i]*(pce[j])[i]*pce_basis_->norm_squared(i);
213 pce_norms[j] = std::sqrt(pce_norms[j]);
214 }
215
216 // Compute F matrix -- PCEs evaluated at all quadrature points
217 // Since F is used in the reduced quadrature below as the quadrature points
218 // for this reduced basis, does scaling by the pce_norms mess up the points?
219 // No -- F essentially defines the random variables this basis is a function
220 // of, and thus they can be scaled in any way we want. Because we don't
221 // explicitly write the basis in terms of F, the scaling is implicit.
222 F_.reshape(nqp, d);
223 Teuchos::Array< Teuchos::Array<value_type> > values(nqp);
224 for (ordinal_type i=0; i<nqp; i++)
225 for (ordinal_type j=0; j<d; j++)
226 F_(i,j) = pce[j].evaluate(points[i], basis_values[i]);
227
228 // Build the reduced basis
229 // Compute basis terms -- 2-D array giving powers for each linear index
230 ordinal_type max_sz;
231 CPBUtils::compute_terms(max_p, d, max_sz, terms_, num_terms_);
232
233 // Compute B matrix -- monomials in F
234 // for i=0,...,nqp-1
235 // for j=0,...,sz-1
236 // B(i,j) = F(i,1)^terms_[j][1] * ... * F(i,d)^terms_[j][d]
237 // where sz is the total size of a basis up to order p and terms_[j]
238 // is an array of powers for each term in the total-order basis
239 SDM B(nqp, max_sz);
240 for (ordinal_type i=0; i<nqp; i++) {
241 for (ordinal_type j=0; j<max_sz; j++) {
242 B(i,j) = 1.0;
243 for (ordinal_type k=0; k<d; k++)
244 B(i,j) *= std::pow(F_(i,k), terms_[j][k]);
245 }
246 }
247
248 // Project B into original basis -- should use SPAM for this
249 SDM Bp(pce_sz_, max_sz);
250 const Teuchos::Array<value_type>& basis_norms =
251 pce_basis_->norm_squared();
252 for (ordinal_type i=0; i<pce_sz_; i++) {
253 for (ordinal_type j=0; j<max_sz; j++) {
254 Bp(i,j) = 0.0;
255 for (ordinal_type k=0; k<nqp; k++)
256 Bp(i,j) += weights[k]*B(k,j)*A_(k,i);
257 Bp(i,j) /= basis_norms[i];
258 }
259 }
260
261 // Rescale columns of Bp to have unit norm
262 for (ordinal_type j=0; j<max_sz; j++) {
263 value_type nrm = 0.0;
264 for (ordinal_type i=0; i<pce_sz_; i++)
265 nrm += Bp(i,j)*Bp(i,j)*basis_norms[i];
266 nrm = std::sqrt(nrm);
267 for (ordinal_type i=0; i<pce_sz_; i++)
268 Bp(i,j) /= nrm;
269 }
270
271 // Compute our new basis -- each column of Qp is the coefficients of the
272 // new basis in the original basis. Constraint pivoting so first d+1
273 // columns and included in Qp.
274 Teuchos::Array<value_type> w(pce_sz_, 1.0);
275 SDM R;
276 Teuchos::Array<ordinal_type> piv(max_sz);
277 for (int i=0; i<d+1; i++)
278 piv[i] = 1;
280 ordinal_type sz_ = SOF::createOrthogonalBasis(
281 orthogonalization_method, threshold, verbose, Bp, w, Qp_, R, piv);
282
283 return sz_;
284}
285
286template <typename ordinal_type, typename value_type>
289{
290}
291
292template <typename ordinal_type, typename value_type>
293ordinal_type
295order() const
296{
297 return p;
298}
299
300template <typename ordinal_type, typename value_type>
301ordinal_type
303dimension() const
304{
305 return d;
306}
307
308template <typename ordinal_type, typename value_type>
309ordinal_type
311size() const
312{
313 return sz;
314}
315
316template <typename ordinal_type, typename value_type>
317const Teuchos::Array<value_type>&
319norm_squared() const
320{
321 return norms;
322}
323
324template <typename ordinal_type, typename value_type>
325const value_type&
327norm_squared(ordinal_type i) const
328{
329 return norms[i];
330}
331
332template <typename ordinal_type, typename value_type>
333Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
336
337{
338 return Teuchos::null;
339}
340
341template <typename ordinal_type, typename value_type>
342Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
345
346{
347 return Teuchos::null;
348}
349
350template <typename ordinal_type, typename value_type>
351value_type
353evaluateZero(ordinal_type i) const
354{
355 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented!");
356}
357
358template <typename ordinal_type, typename value_type>
359void
361evaluateBases(const Teuchos::ArrayView<const value_type>& point,
362 Teuchos::Array<value_type>& basis_vals) const
363{
364 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented!");
365}
366
367template <typename ordinal_type, typename value_type>
368const std::string&
370getName() const
371{
372 return name;
373}
374
375template <typename ordinal_type, typename value_type>
376void
378print(std::ostream& os) const
379{
380 os << "Gram-Schmidt basis of order " << p << ", dimension " << d
381 << ", and size " << sz << ". Matrix coefficients:\n";
382 os << printMat(Qp) << std::endl;
383 os << "Basis vector norms (squared):\n\t";
384 for (ordinal_type i=0; i<sz; i++)
385 os << norms[i] << " ";
386 os << "\n";
387}
388
389template <typename ordinal_type, typename value_type>
390void
392transformToOriginalBasis(const value_type *in, value_type *out,
393 ordinal_type ncol, bool transpose) const
394{
395 if (transpose) {
396 SDM zbar(Teuchos::View, const_cast<value_type*>(in), ncol, ncol, sz);
397 SDM z(Teuchos::View, out, ncol, ncol, pce_sz);
398 ordinal_type ret =
399 z.multiply(Teuchos::NO_TRANS, Teuchos::TRANS, 1.0, zbar, Qp, 0.0);
400 TEUCHOS_ASSERT(ret == 0);
401 }
402 else {
403 SDM zbar(Teuchos::View, const_cast<value_type*>(in), sz, sz, ncol);
404 SDM z(Teuchos::View, out, pce_sz, pce_sz, ncol);
405 ordinal_type ret =
406 z.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Qp, zbar, 0.0);
407 TEUCHOS_ASSERT(ret == 0);
408 }
409}
410
411template <typename ordinal_type, typename value_type>
412void
414transformFromOriginalBasis(const value_type *in, value_type *out,
415 ordinal_type ncol, bool transpose) const
416{
417 if (transpose) {
418 SDM z(Teuchos::View, const_cast<value_type*>(in), ncol, ncol, pce_sz);
419 SDM zbar(Teuchos::View, out, ncol, ncol, sz);
420 ordinal_type ret =
421 zbar.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, z, Qp, 0.0);
422 TEUCHOS_ASSERT(ret == 0);
423 }
424 else {
425 SDM z(Teuchos::View, const_cast<value_type*>(in), pce_sz, pce_sz, ncol);
426 SDM zbar(Teuchos::View, out, sz, sz, ncol);
427 ordinal_type ret =
428 zbar.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Qp, z, 0.0);
429 TEUCHOS_ASSERT(ret == 0);
430 }
431}
432
433template <typename ordinal_type, typename value_type>
434Teuchos::RCP<const Stokhos::Quadrature<ordinal_type, value_type> >
437{
438 return reduced_quad;
439}
int quad2(const Epetra_Map &map, bool verbose)
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
SDM Q
Values of transformed basis at quadrature points.
Teuchos::Array< Stokhos::MultiIndex< ordinal_type > > terms
2-D array of basis terms
virtual void transformFromOriginalBasis(const value_type *in, value_type *out, ordinal_type ncol=1, bool transpose=false) const
Transform coefficients from original basis to this basis.
virtual const std::string & getName() const
Return string name of basis.
virtual void transformToOriginalBasis(const value_type *in, value_type *out, ordinal_type ncol=1, bool transpose=false) const
Transform coefficients to original basis from this basis.
virtual ordinal_type size() const
Return total size of basis.
virtual Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > getReducedQuadrature() const
Get reduced quadrature object.
Teuchos::Array< ordinal_type > num_terms
Number of terms up to each order.
ordinal_type buildQ(ordinal_type max_p, value_type threshold, const Teuchos::Array< Stokhos::OrthogPolyApprox< ordinal_type, value_type > > &pce, const Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > &quad, Teuchos::Array< Stokhos::MultiIndex< ordinal_type > > &terms_, Teuchos::Array< ordinal_type > &num_terms_, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Qp_, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &A_, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &F_)
Build the reduced basis, parameterized by total order max_p.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > SDM
virtual void print(std::ostream &os) const
Print basis to stream os.
Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > reduced_quad
Reduced quadrature object.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
virtual const Teuchos::Array< value_type > & norm_squared() const
Return array storing norm-squared of each basis polynomial.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > pce_basis
Original pce basis.
virtual value_type evaluateZero(ordinal_type i) const
Evaluate basis polynomial i at zero.
MonomialProjGramSchmidtPCEBasis2(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::ParameterList &params=Teuchos::ParameterList())
Constructor.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeLinearTripleProductTensor() const
Compute linear triple product tensor where k = 0,1,..,d.
virtual void evaluateBases(const Teuchos::ArrayView< const value_type > &point, Teuchos::Array< value_type > &basis_vals) const
Evaluate basis polynomials at given point point.
SDM Qp
Coefficients of transformed basis in original basis.
A multidimensional index.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
Encapsulate various orthogonalization (ie QR) methods.
Abstract base class for quadrature methods.
Generate a basis from a given set of PCE expansions that is orthogonal with respect to the product me...
virtual Teuchos::RCP< const Stokhos::UserDefinedQuadrature< ordinal_type, value_type > > createReducedQuadrature(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Q, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Q2, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &F, const Teuchos::Array< value_type > &weights) const
Get reduced quadrature object.
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
void printMat(const char *name, Epetra_IntSerialDenseMatrix &matrix)