Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
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 
44 #include "Stokhos_ProductBasis.hpp"
45 
46 template <typename ordinal_type, typename value_type, typename storage_type>
47 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
48 OrthogPolyApprox(
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 
65 template <typename ordinal_type, typename value_type, typename storage_type>
66 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
67 OrthogPolyApprox(const Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>& x) :
68  basis_(x.basis_),
69  coeff_(x.coeff_)
70 {
71 }
72 
73 template <typename ordinal_type, typename value_type, typename storage_type>
74 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
75 ~OrthogPolyApprox()
76 {
77 }
78 
79 template <typename ordinal_type, typename value_type, typename storage_type>
80 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>&
81 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
82 operator=(const Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>& x)
83 {
84  if (this != &x) {
85  basis_ = x.basis_;
86  coeff_ = x.coeff_;
87  }
88 
89  return *this;
90 }
91 
92 template <typename ordinal_type, typename value_type, typename storage_type>
93 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>&
94 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
95 operator=(const value_type& v)
96 {
97  coeff_.init(value_type(0));
98  coeff_.init(&v, 1);
99 
100  return *this;
101 }
102 
103 template <typename ordinal_type, typename value_type, typename storage_type>
104 void
105 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
106 init(const value_type& v)
107 {
108  coeff_.init(v);
109 }
110 
111 template <typename ordinal_type, typename value_type, typename storage_type>
112 void
113 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
114 init(const value_type* v)
115 {
116  coeff_.init(v);
117 }
118 
119 template <typename ordinal_type, typename value_type, typename storage_type>
120 void
121 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
122 load(value_type* v)
123 {
124  coeff_.load(v);
125 }
126 
127 template <typename ordinal_type, typename value_type, typename storage_type>
128 Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type, value_type> >
129 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
130 basis() const
131 {
132  return basis_;
133 }
134 
135 template <typename ordinal_type, typename value_type, typename storage_type>
136 void
137 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
138 reset(const Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type, value_type> >& new_basis, ordinal_type sz)
139 {
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 }
148 
149 template <typename ordinal_type, typename value_type, typename storage_type>
150 void
151 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
152 resize(ordinal_type sz)
153 {
154  coeff_.resize(sz);
155 }
156 
157 template <typename ordinal_type, typename value_type, typename storage_type>
158 ordinal_type
159 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
160 size() const
161 {
162  return coeff_.size();
163 }
164 
165 template <typename ordinal_type, typename value_type, typename storage_type>
166 typename Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::pointer
167 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
168 coeff()
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();
176 }
177 
178 template <typename ordinal_type, typename value_type, typename storage_type>
179 typename Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::const_pointer
180 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
181 coeff() 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 
191 template <typename ordinal_type, typename value_type, typename storage_type>
192 typename Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::reference
193 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
194 operator[](ordinal_type i)
195 {
196  return coeff_[i];
197 }
198 
199 template <typename ordinal_type, typename value_type, typename storage_type>
200 typename Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::const_reference
201 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
202 operator[](ordinal_type i) const
203 {
204  return coeff_[i];
205 }
206 
207 template <typename ordinal_type, typename value_type, typename storage_type>
208 typename Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::reference
209 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
210 term(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 
221 template <typename ordinal_type, typename value_type, typename storage_type>
222 typename Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::const_reference
223 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
224 term(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 
235 template <typename ordinal_type, typename value_type, typename storage_type>
236 const Stokhos::MultiIndex<ordinal_type>&
237 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
238 order(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 
245 template <typename ordinal_type, typename value_type, typename storage_type>
246 value_type
247 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
248 evaluate(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 
255 template <typename ordinal_type, typename value_type, typename storage_type>
256 value_type
257 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
258 evaluate(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 
268 template <typename ordinal_type, typename value_type, typename storage_type>
269 value_type
270 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
271 mean() const
272 {
273  return coeff_[0];
274 }
275 
276 template <typename ordinal_type, typename value_type, typename storage_type>
277 value_type
278 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
279 standard_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 
289 template <typename ordinal_type, typename value_type, typename storage_type>
290 value_type
291 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
292 two_norm() const
293 {
294  return std::sqrt(this->two_norm_squared());
295 }
296 
297 template <typename ordinal_type, typename value_type, typename storage_type>
298 value_type
299 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
300 two_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 
316 template <typename ordinal_type, typename value_type, typename storage_type>
317 value_type
318 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
319 inner_product(const Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>& b) const
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 
343 template <typename ordinal_type, typename value_type, typename storage_type>
344 std::ostream&
345 Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>::
346 print(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 
376 template <typename ordinal_type, typename value_type, typename storage_type>
377 std::ostream&
378 Stokhos::operator << (
379  std::ostream& os,
380  const Stokhos::OrthogPolyApprox<ordinal_type, value_type, storage_type>& a)
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 }
Stokhos::OrthogPolyApprox::two_norm_squared
value_type two_norm_squared() const
Compute the squared two-norm of expansion.
Definition: Stokhos_OrthogPolyApproxImp.hpp:300
Stokhos::StandardStorage::resize
void resize(const ordinal_type &sz)
Resize to new size (values are preserved)
Definition: Stokhos_StandardStorage.hpp:94
Stokhos::OrthogPolyApprox::const_reference
storage_type::const_reference const_reference
Definition: Stokhos_OrthogPolyApprox.hpp:67
Stokhos::OrthogPolyApprox::resize
void resize(ordinal_type sz)
Resize coefficient array (coefficients are preserved)
Definition: Stokhos_OrthogPolyApproxImp.hpp:152
Stokhos::StandardStorage::size
ordinal_type size() const
Return size.
Definition: Stokhos_StandardStorage.hpp:97
Stokhos::OrthogPolyApprox::size
ordinal_type size() const
Return size.
Definition: Stokhos_OrthogPolyApproxImp.hpp:160
Stokhos::OrthogPolyApprox::order
const MultiIndex< ordinal_type > & order(ordinal_type term) const
Get orders for a given term.
Definition: Stokhos_OrthogPolyApproxImp.hpp:238
Stokhos::MultiIndex::size
ordinal_type size() const
Size.
Definition: Stokhos_ProductBasisUtils.hpp:99
Stokhos::OrthogPolyApprox::init
void init(const value_type &v)
Initialize coefficients to value.
Definition: Stokhos_OrthogPolyApproxImp.hpp:106
Stokhos::OrthogPolyApprox::load
void load(value_type *v)
Load coefficients to an array of values.
Definition: Stokhos_OrthogPolyApproxImp.hpp:122
TotalOrderBasisUnitTest::value_type
double value_type
Definition: Stokhos_LexicographicTreeBasisUnitTest.cpp:70
Stokhos::OrthogPolyApprox::const_pointer
storage_type::const_pointer const_pointer
Definition: Stokhos_OrthogPolyApprox.hpp:69
Stokhos::OrthogPolyApprox::basis
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > basis() const
Return basis.
Definition: Stokhos_OrthogPolyApproxImp.hpp:130
Stokhos::OrthogPolyApprox::standard_deviation
value_type standard_deviation() const
Compute standard deviation of expansion.
Definition: Stokhos_OrthogPolyApproxImp.hpp:279
Stokhos::operator<<
std::ostream & operator<<(std::ostream &os, const ProductContainer< coeff_type > &vec)
Definition: Stokhos_ProductContainer.hpp:184
Stokhos::OrthogPolyApprox::mean
value_type mean() const
Compute mean of expansion.
Definition: Stokhos_OrthogPolyApproxImp.hpp:271
Stokhos::OrthogPolyApprox::coeff_
storage_type coeff_
OrthogPolyApprox coefficients.
Definition: Stokhos_OrthogPolyApprox.hpp:184
Stokhos::MultiIndex< ordinal_type >
val
expr val()
Stokhos::OrthogPolyApprox::operator[]
reference operator[](ordinal_type i)
Array access.
Definition: Stokhos_OrthogPolyApproxImp.hpp:194
TotalOrderBasisUnitTest::ordinal_type
int ordinal_type
Definition: Stokhos_LexicographicTreeBasisUnitTest.cpp:69
Stokhos::OrthogPolyApprox::basis_
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > basis_
Basis expansion is relative to.
Definition: Stokhos_OrthogPolyApprox.hpp:181
Stokhos::OrthogPolyApprox::evaluate
value_type evaluate(const Teuchos::Array< value_type > &point) const
Evaluate polynomial approximation at a point.
Definition: Stokhos_OrthogPolyApproxImp.hpp:248
Stokhos::OrthogPolyApprox::coeff
pointer coeff()
Return coefficient array.
Definition: Stokhos_OrthogPolyApproxImp.hpp:168
Sacado::UQ::sqrt
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:878
Stokhos_ProductBasis.hpp
j
j
Definition: Sacado_Fad_Exp_MP_Vector.hpp:527
Stokhos::OrthogPolyApprox::reset
void reset(const Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > &new_basis, ordinal_type sz=0)
Reset to a new basis.
Definition: Stokhos_OrthogPolyApproxImp.hpp:138
Stokhos::OrthogPolyApprox::reference
storage_type::reference reference
Definition: Stokhos_OrthogPolyApprox.hpp:66
Stokhos::OrthogPolyApprox::print
std::ostream & print(std::ostream &os) const
Print approximation in basis.
Definition: Stokhos_OrthogPolyApproxImp.hpp:346
Stokhos::OrthogPolyApprox::OrthogPolyApprox
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.
Definition: Stokhos_OrthogPolyApproxImp.hpp:48
Stokhos::OrthogPolyApprox::~OrthogPolyApprox
~OrthogPolyApprox()
Destructor.
Definition: Stokhos_OrthogPolyApproxImp.hpp:75
cusp::detail::device::x
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:265
Stokhos::OrthogPolyApprox
Class to store coefficients of a projection onto an orthogonal polynomial basis.
Definition: Stokhos_OrthogPolyApprox.hpp:63
Stokhos::OrthogPolyApprox::pointer
storage_type::pointer pointer
Definition: Stokhos_OrthogPolyApprox.hpp:68
Stokhos::OrthogPolyApprox::inner_product
value_type inner_product(const OrthogPolyApprox &b) const
Compute the L2 inner product of 2 PCEs.
Definition: Stokhos_OrthogPolyApproxImp.hpp:319
Stokhos::StandardStorage::init
void init(const_reference v)
Initialize values to a constant value.
Definition: Stokhos_StandardStorage.hpp:76
Stokhos::OrthogPolyApprox::two_norm
value_type two_norm() const
Compute the two-norm of expansion.
Definition: Stokhos_OrthogPolyApproxImp.hpp:292
Stokhos::OrthogPolyApprox::term
reference term(ordinal_type dimension, ordinal_type order)
Get coefficient term for given dimension and order.
Definition: Stokhos_OrthogPolyApproxImp.hpp:210
Stokhos::OrthogPolyApprox::operator=
OrthogPolyApprox & operator=(const OrthogPolyApprox &x)
Assignment operator (deep copy)
Stokhos::OrthogPolyBasis
Abstract base class for multivariate orthogonal polynomials.
Definition: Stokhos_OrthogPolyBasis.hpp:74
Generated by doxygen 1.8.20