Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_OrthogPolyApproxImp.hpp
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
45
46template <typename ordinal_type, typename value_type, typename storage_type>
49 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type, value_type> >& theBasis,
50 ordinal_type sz,
51 const value_type* vals) :
52 basis_(theBasis),
53 coeff_(0)
54{
55 if (sz != 0)
56 coeff_.resize(sz);
57 else if (theBasis != Teuchos::null)
58 coeff_.resize(theBasis->size());
59 else
60 coeff_.resize(ordinal_type(1));
61 if (vals != NULL)
62 coeff_.init(vals);
63}
64
65template <typename ordinal_type, typename value_type, typename storage_type>
68 basis_(x.basis_),
69 coeff_(x.coeff_)
70{
71}
72
73template <typename ordinal_type, typename value_type, typename storage_type>
76{
77}
79template <typename ordinal_type, typename value_type, typename storage_type>
83{
84 if (this != &x) {
85 basis_ = x.basis_;
86 coeff_ = x.coeff_;
87 }
88
89 return *this;
90}
91
92template <typename ordinal_type, typename value_type, typename storage_type>
95operator=(const value_type& v)
97 coeff_.init(value_type(0));
98 coeff_.init(&v, 1);
100 return *this;
101}
102
103template <typename ordinal_type, typename value_type, typename storage_type>
104void
106init(const value_type& v)
107{
108 coeff_.init(v);
109}
110
111template <typename ordinal_type, typename value_type, typename storage_type>
112void
114init(const value_type* v)
115{
116 coeff_.init(v);
117}
119template <typename ordinal_type, typename value_type, typename storage_type>
120void
122load(value_type* v)
123{
124 coeff_.load(v);
125}
126
127template <typename ordinal_type, typename value_type, typename storage_type>
128Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type, value_type> >
130basis() const
131{
132 return basis_;
134
135template <typename ordinal_type, typename value_type, typename storage_type>
136void
138reset(const Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type, value_type> >& new_basis, ordinal_type sz)
140 basis_ = new_basis;
141 if (sz != 0)
142 resize(sz);
143 else if (basis_ != Teuchos::null)
144 resize(basis_->size());
145 else
146 resize(1);
147}
149template <typename ordinal_type, typename value_type, typename storage_type>
150void
152resize(ordinal_type sz)
153{
154 coeff_.resize(sz);
155}
156
157template <typename ordinal_type, typename value_type, typename storage_type>
158ordinal_type
160size() const
162 return coeff_.size();
163}
165template <typename ordinal_type, typename value_type, typename storage_type>
168coeff()
169{
170#ifdef STOKHOS_DEBUG
171 TEUCHOS_TEST_FOR_EXCEPTION(coeff_.size() == 0, std::logic_error,
172 "Stokhos::OrthogPolyApprox::coeff(): " <<
173 "Coefficient array is empty!");
174#endif
175 return coeff_.coeff();
177
178template <typename ordinal_type, typename value_type, typename storage_type>
181coeff() const
182{
183#ifdef STOKHOS_DEBUG
184 TEUCHOS_TEST_FOR_EXCEPTION(coeff_.size() == 0, std::logic_error,
185 "Stokhos::OrthogPolyApprox::coeff(): " <<
186 "Coefficient array is empty!");
187#endif
188 return coeff_.coeff();
189}
190
191template <typename ordinal_type, typename value_type, typename storage_type>
194operator[](ordinal_type i)
195{
196 return coeff_[i];
197}
198
199template <typename ordinal_type, typename value_type, typename storage_type>
202operator[](ordinal_type i) const
203{
204 return coeff_[i];
205}
206
207template <typename ordinal_type, typename value_type, typename storage_type>
210term(ordinal_type dimension, ordinal_type theOrder)
211{
212 Teuchos::RCP< const Stokhos::ProductBasis<ordinal_type, value_type> >
213 product_basis = Teuchos::rcp_dynamic_cast< const Stokhos::ProductBasis<ordinal_type, value_type> >(basis_, true);
214 ordinal_type d = basis_->dimension();
215 MultiIndex<ordinal_type> theTerm(d);
216 theTerm[dimension] = theOrder;
217 ordinal_type index = product_basis->index(theTerm);
218 return coeff_[index];
219}
220
221template <typename ordinal_type, typename value_type, typename storage_type>
224term(ordinal_type dimension, ordinal_type theOrder) const
225{
226 Teuchos::RCP< const Stokhos::ProductBasis<ordinal_type, value_type> >
227 product_basis = Teuchos::rcp_dynamic_cast< const Stokhos::ProductBasis<ordinal_type, value_type> >(basis_, true);
228 ordinal_type d = basis_->dimension();
229 MultiIndex<ordinal_type> theTerm(d);
230 theTerm[dimension] = theOrder;
231 ordinal_type index = product_basis->index(theTerm);
232 return coeff_[index];
233}
234
235template <typename ordinal_type, typename value_type, typename storage_type>
238order(ordinal_type theTerm) const
239{
240 Teuchos::RCP< const Stokhos::ProductBasis<ordinal_type, value_type> >
241 product_basis = Teuchos::rcp_dynamic_cast< const Stokhos::ProductBasis<ordinal_type, value_type> >(basis_, true);
242 return product_basis->term(theTerm);
243}
244
245template <typename ordinal_type, typename value_type, typename storage_type>
246value_type
248evaluate(const Teuchos::Array<value_type>& point) const
249{
250 Teuchos::Array<value_type> basis_vals(basis_->size());
251 basis_->evaluateBases(point, basis_vals);
252 return evaluate(point, basis_vals);
253}
254
255template <typename ordinal_type, typename value_type, typename storage_type>
256value_type
258evaluate(const Teuchos::Array<value_type>& point,
259 const Teuchos::Array<value_type>& basis_vals) const
260{
261 value_type val = value_type(0.0);
262 for (ordinal_type i=0; i<static_cast<ordinal_type>(coeff_.size()); i++)
263 val += coeff_[i]*basis_vals[i];
264
265 return val;
266}
267
268template <typename ordinal_type, typename value_type, typename storage_type>
269value_type
271mean() const
272{
273 return coeff_[0];
274}
275
276template <typename ordinal_type, typename value_type, typename storage_type>
277value_type
279standard_deviation() const
280{
281 value_type std_dev = 0.0;
282 for (ordinal_type i=1; i<static_cast<ordinal_type>(coeff_.size()); i++) {
283 std_dev += coeff_[i]*coeff_[i]*basis_->norm_squared(i);
284 }
285 std_dev = std::sqrt(std_dev);
286 return std_dev;
287}
288
289template <typename ordinal_type, typename value_type, typename storage_type>
290value_type
292two_norm() const
293{
294 return std::sqrt(this->two_norm_squared());
295}
296
297template <typename ordinal_type, typename value_type, typename storage_type>
298value_type
300two_norm_squared() const
301{
302 value_type nrm = 0.0;
303 if (basis_ == Teuchos::null) { // Check for special case of constants
304 TEUCHOS_TEST_FOR_EXCEPTION(
305 coeff_.size() != 1, std::logic_error,
306 "basis_ == null && coeff_.size() > 1");
307 nrm = coeff_[0]*coeff_[0];
308 }
309 else {
310 for (ordinal_type i=0; i<static_cast<ordinal_type>(coeff_.size()); i++)
311 nrm += coeff_[i]*coeff_[i]*basis_->norm_squared(i);
312 }
313 return nrm;
314}
315
316template <typename ordinal_type, typename value_type, typename storage_type>
317value_type
320{
321 // Check a and b are compatible
322 TEUCHOS_TEST_FOR_EXCEPTION(
323 basis_ == Teuchos::null && coeff_.size() != 1, std::logic_error,
324 "basis_ == null && coeff_.size() > 1");
325 TEUCHOS_TEST_FOR_EXCEPTION(
326 b.basis_ == Teuchos::null && b.coeff_.size() != 1, std::logic_error,
327 "b.basis_ == null && b.coeff_.size() > 1");
328 TEUCHOS_TEST_FOR_EXCEPTION(
329 coeff_.size() != b.coeff_.size() &&
330 coeff_.size() != 1 && b.coeff_.size() != 1, std::logic_error,
331 "Coefficient array sizes do not match");
332
333 value_type v = 0.0;
334 if (coeff_.size() == 1 || b.coeff_.size() == 1)
335 v = coeff_[0]*b.coeff_[0];
336 else
337 for (ordinal_type i=0; i<static_cast<ordinal_type>(coeff_.size()); i++)
338 v += coeff_[i]*b.coeff_[i]*basis_->norm_squared(i);
339
340 return v;
341}
342
343template <typename ordinal_type, typename value_type, typename storage_type>
344std::ostream&
346print(std::ostream& os) const
347{
348 os << "Stokhos::OrthogPolyApprox of size " << coeff_.size() << " in basis "
349 << "\n\t" << basis_->getName() << ":" << std::endl;
350
351 Teuchos::RCP< const Stokhos::ProductBasis<ordinal_type, value_type> >
352 product_basis = Teuchos::rcp_dynamic_cast< const Stokhos::ProductBasis<ordinal_type, value_type> >(basis_);
353
354 if (product_basis != Teuchos::null) {
355 for (ordinal_type i=0; i<static_cast<ordinal_type>(coeff_.size()); i++) {
356 const Stokhos::MultiIndex<ordinal_type>& trm = product_basis->term(i);
357 os << "\t\t(";
358 for (ordinal_type j=0; j<static_cast<ordinal_type>(trm.size())-1; j++)
359 os << trm[j] << ", ";
360 os << trm[trm.size()-1] << ") = " << coeff_[i] << std::endl;
361 }
362 }
363 else {
364 os << "[ ";
365
366 for (ordinal_type i=0; i<static_cast<ordinal_type>(coeff_.size()); i++) {
367 os << coeff_[i] << " ";
368 }
369
370 os << "]\n";
371 }
372
373 return os;
374}
375
376template <typename ordinal_type, typename value_type, typename storage_type>
377std::ostream&
379 std::ostream& os,
381{
382 os << "[ ";
383
384 for (ordinal_type i=0; i<static_cast<ordinal_type>(a.size()); i++) {
385 os << a[i] << " ";
386 }
387
388 os << "]\n";
389 return os;
390}
expr val()
A multidimensional index.
ordinal_type size() const
Size.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
storage_type::reference reference
void reset(const Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > &new_basis, ordinal_type sz=0)
Reset to a new basis.
value_type inner_product(const OrthogPolyApprox &b) const
Compute the L2 inner product of 2 PCEs.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > basis_
Basis expansion is relative to.
value_type mean() const
Compute mean of expansion.
storage_type::const_reference const_reference
std::ostream & print(std::ostream &os) const
Print approximation in basis.
void init(const value_type &v)
Initialize coefficients to value.
value_type two_norm() const
Compute the two-norm of expansion.
OrthogPolyApprox & operator=(const OrthogPolyApprox &x)
Assignment operator (deep copy)
value_type standard_deviation() const
Compute standard deviation of expansion.
storage_type coeff_
OrthogPolyApprox coefficients.
void resize(ordinal_type sz)
Resize coefficient array (coefficients are preserved)
const MultiIndex< ordinal_type > & order(ordinal_type term) const
Get orders for a given term.
void load(value_type *v)
Load coefficients to an array of values.
OrthogPolyApprox(const Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > &basis=Teuchos::null, ordinal_type sz=0, const value_type *vals=NULL)
Constructor with supplied size sz.
reference operator[](ordinal_type i)
Array access.
pointer coeff()
Return coefficient array.
value_type two_norm_squared() const
Compute the squared two-norm of expansion.
storage_type::const_pointer const_pointer
reference term(ordinal_type dimension, ordinal_type order)
Get coefficient term for given dimension and order.
ordinal_type size() const
Return size.
value_type evaluate(const Teuchos::Array< value_type > &point) const
Evaluate polynomial approximation at a point.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > basis() const
Return basis.
Abstract base class for multivariate orthogonal polynomials.
ordinal_type size() const
Return size.
void init(const_reference v)
Initialize values to a constant value.
void resize(const ordinal_type &sz)
Resize to new size (values are preserved)
std::ostream & operator<<(std::ostream &os, const ProductContainer< coeff_type > &vec)