Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
43 #include "Stokhos_BasisFactory.hpp"
47 
48 template <typename ordinal_type, typename value_type>
51  ordinal_type max_p,
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
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 =
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") {
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
103  if (2*max_p > pce_basis->order()) {
104 
105  // Basis
107  ordinal_type dim = prod_basis->dimension();
109  for (ordinal_type i=0; i<dim; i++)
110  bases[i] = prod_basis->getCoordinateBases()[i]->cloneWithOrder(2*max_p);
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();
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 
178 template <typename ordinal_type, typename value_type>
182  ordinal_type max_p,
183  value_type threshold,
187  Teuchos::Array<ordinal_type>& num_terms_,
188  SDM& Qp_, SDM& A_, SDM& F_)
189 {
191  ordinal_type pce_sz_ = pce_basis_->size();
192 
193  // Get quadrature data
194  const Teuchos::Array<value_type>& weights = quad->getQuadWeights();
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);
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 
286 template <typename ordinal_type, typename value_type>
289 {
290 }
291 
292 template <typename ordinal_type, typename value_type>
295 order() const
296 {
297  return p;
298 }
299 
300 template <typename ordinal_type, typename value_type>
303 dimension() const
304 {
305  return d;
306 }
307 
308 template <typename ordinal_type, typename value_type>
311 size() const
312 {
313  return sz;
314 }
315 
316 template <typename ordinal_type, typename value_type>
320 {
321  return norms;
322 }
323 
324 template <typename ordinal_type, typename value_type>
325 const value_type&
328 {
329  return norms[i];
330 }
331 
332 template <typename ordinal_type, typename value_type>
336 
337 {
338  return Teuchos::null;
339 }
340 
341 template <typename ordinal_type, typename value_type>
345 
346 {
347  return Teuchos::null;
348 }
349 
350 template <typename ordinal_type, typename value_type>
354 {
355  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented!");
356 }
357 
358 template <typename ordinal_type, typename value_type>
359 void
362  Teuchos::Array<value_type>& basis_vals) const
363 {
364  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented!");
365 }
366 
367 template <typename ordinal_type, typename value_type>
368 const std::string&
370 getName() const
371 {
372  return name;
373 }
374 
375 template <typename ordinal_type, typename value_type>
376 void
378 print(std::ostream& os) const
379 {
380  os << "Gram-Schmidt basis of order " << p << ", dimension " << d
381  << ", and size " << sz << ". Matrix coefficients:\n";
382  os << 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 
389 template <typename ordinal_type, typename value_type>
390 void
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 
411 template <typename ordinal_type, typename value_type>
412 void
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 
433 template <typename ordinal_type, typename value_type>
437 {
438  return reduced_quad;
439 }
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
virtual void evaluateBases(const Teuchos::ArrayView< const value_type > &point, Teuchos::Array< value_type > &basis_vals) const
Evaluate basis polynomials at given point point.
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.
virtual const std::string & getName() const
Return string name of basis.
Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > reduced_quad
Reduced quadrature object.
virtual Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > getReducedQuadrature() const
Get reduced quadrature object.
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.
KOKKOS_INLINE_FUNCTION PCE< Storage > pow(const PCE< Storage > &a, const PCE< Storage > &b)
T & get(ParameterList &l, const std::string &name)
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
RCP< ParameterList > sublist(const RCP< ParameterList > &paramList, const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
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.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
int quad2(const Epetra_Map &map, bool verbose)
bool isParameter(const std::string &name) const
SDM Q
Values of transformed basis at quadrature points.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual value_type evaluateZero(ordinal_type i) const
Evaluate basis polynomial i at zero.
Abstract base class for quadrature methods.
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.
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
void resize(size_type new_size, const value_type &x=value_type())
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Generate a basis from a given set of PCE expansions that is orthogonal with respect to the product me...
void push_back(const value_type &x)
int reshape(OrdinalType numRows, OrdinalType numCols)
virtual const Teuchos::Array< value_type > & norm_squared() const
Return array storing norm-squared of each basis polynomial.
size_type size() const
ordinal_type dimension() const
Return dimension of basis.
SDM Qp
Coefficients of transformed basis in original basis.
virtual ordinal_type size() const
Return total size of basis.
Teuchos::Array< Stokhos::MultiIndex< ordinal_type > > terms
2-D array of basis terms
#define TEUCHOS_ASSERT(assertion_test)
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 Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeLinearTripleProductTensor() const
Compute linear triple product tensor where k = 0,1,..,d.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > pce_basis
Original pce basis.
virtual void print(std::ostream &os) const
Print basis to stream os.
Teuchos::Array< ordinal_type > num_terms
Number of terms up to each order.
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules...
OrdinalType numRows() const
Encapsulate various orthogonalization (ie QR) methods.