Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_SmolyakBasisImp.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
42#include "Teuchos_TimeMonitor.hpp"
43#include "Teuchos_TestForException.hpp"
44
45template <typename ordinal_type, typename value_type, typename ordering_type>
46template <typename index_set_type>
49 const Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type, value_type> > >& bases_,
50 const index_set_type& index_set,
51 const value_type& sparse_tol_,
52 const ordering_type& coeff_compare) :
53 p(0),
54 d(bases_.size()),
55 sz(0),
56 bases(bases_),
57 sparse_tol(sparse_tol_),
58 max_orders(d),
59 basis_set(coeff_compare),
60 norms()
61{
62
63 // Generate index set for the final Smolyak coefficients
64 //
65 // The Smolyak operator is given by the formula
66 //
67 // A = \sum_{k\in\K} \bigotimes_{i=1}^d \Delta^i_{k_i}
68 //
69 // where \Delta^i_0 = 0, \Delta^i_{k_i} = L^i_{k_i} - L^i_{k_i-1},
70 // and K is the supplied index set. This becomes
71 //
72 // A = \sum_{k\in\tilde{K}} c(k) \bigotimes_{i=1}^d L^i_{k_i}
73 //
74 // for some new index set \tilde{K} and coefficient c(k). Using the
75 // formula (cf. G W Wasilkowski and H Wozniakowski, "Explicit cost bounds
76 // of algorithms for multivariate tensor product problems,"
77 // Journal of Complexity (11), 1995)
78 //
79 // \bigotimes_{i=1}^d \Delta^i_{k_i} =
80 // \sum_{\alpha\in\Alpha} (-1)^{|\alpha|}
81 // \bigotimes_{i=1}^d L^i_{k_i-\alpha_i}
82 //
83 // where \Alpha = {0,1}^d and |\alpha| = \alpha_1 + ... + \alpha_d, we
84 // iterate over K and \Alpha, compute k-\alpha and the corresponding
85 // coefficient contribution (-1)^{|\alpha|} and store these in a map.
86 // The keys of of this map with non-zero coefficients define
87 // \tilde{K} and c(k).
90 typedef std::map<multiindex_type,ordinal_type,index_compare> index_map_type;
91 ordinal_type dim = index_set.dimension();
92 alpha_set_type alpha_set(dim, 1);
93 typename alpha_set_type::iterator alpha_begin = alpha_set.begin();
94 typename alpha_set_type::iterator alpha_end = alpha_set.end();
95 typename index_set_type::iterator index_iterator = index_set.begin();
96 typename index_set_type::iterator index_end = index_set.end();
97 multiindex_type diff(dim);
98 index_map_type index_map;
99 for (; index_iterator != index_end; ++index_iterator) {
100 for (typename alpha_set_type::iterator alpha = alpha_begin;
101 alpha != alpha_end; ++alpha) {
102 bool valid_index = true;
103 for (ordinal_type i=0; i<dim; ++i) {
104 diff[i] = (*index_iterator)[i] - (*alpha)[i];
105 if (diff[i] < 0) {
106 valid_index = false;
107 break;
108 }
109 }
110 if (valid_index) {
111 ordinal_type alpha_order = alpha->order();
112 ordinal_type val;
113 if (alpha_order % 2 == 0)
114 val = 1;
115 else
116 val = -1;
117 typename index_map_type::iterator index_map_iterator =
118 index_map.find(diff);
119 if (index_map_iterator == index_map.end())
120 index_map[diff] = val;
121 else
122 index_map_iterator->second += val;
123 }
124 }
125 }
126
127 // Generate tensor product bases
128 typename index_map_type::iterator index_map_iterator = index_map.begin();
129 typename index_map_type::iterator index_map_end = index_map.end();
130 for (; index_map_iterator != index_map_end; ++index_map_iterator) {
131
132 // Skip indices with zero coefficient
133 if (index_map_iterator->second == 0)
134 continue;
135
136 // Apply growth rule to cofficient multi-index
137 multiindex_type coeff_growth_index(dim);
138 for (ordinal_type i=0; i<dim; ++i) {
139 coeff_growth_index[i] =
140 bases[i]->coefficientGrowth(index_map_iterator->first[i]);
141 }
142
143 // Build tensor product basis for given index
144 Teuchos::RCP<tensor_product_basis_type> tp =
145 Teuchos::rcp(new tensor_product_basis_type(
146 bases, sparse_tol, coeff_growth_index));
147
148 // Include coefficients in union over all sets
149 for (ordinal_type i=0; i<tp->size(); ++i)
150 basis_set[tp->term(i)] = ordinal_type(0);
151
152 tp_bases.push_back(tp);
153 sm_pred.tp_preds.push_back(
154 TensorProductPredicate<ordinal_type>(coeff_growth_index) );
155 smolyak_coeffs.push_back(index_map_iterator->second);
156 }
157 sz = basis_set.size();
158
159 // Generate linear odering of coefficients
160 ordinal_type idx = 0;
161 basis_map.resize(sz);
162 for (typename coeff_set_type::iterator i = basis_set.begin();
163 i != basis_set.end();
164 ++i) {
165 i->second = idx;
166 basis_map[idx] = i->first;
167 ++idx;
168 }
169
170 // Compute max coefficient orders
171 for (ordinal_type i=0; i<sz; ++i) {
172 for (ordinal_type j=0; j<dim; ++j)
173 if (basis_map[i][j] > max_orders[j])
174 max_orders[j] = basis_map[i][j];
175 }
176
177 // Resize bases to make sure they are high enough order
178 for (ordinal_type i=0; i<dim; i++)
179 if (bases[i]->order() < max_orders[i])
180 bases[i] = bases[i]->cloneWithOrder(max_orders[i]);
181
182 // Compute largest order
183 p = 0;
184 for (ordinal_type i=0; i<d; i++) {
185 if (max_orders[i] > p)
186 p = max_orders[i];
187 }
188
189 // Compute norms
190 norms.resize(sz);
191 value_type nrm;
192 for (ordinal_type k=0; k<sz; k++) {
193 nrm = value_type(1.0);
194 for (ordinal_type i=0; i<d; i++)
195 nrm = nrm * bases[i]->norm_squared(basis_map[k][i]);
196 norms[k] = nrm;
197 }
198
199 // Create name
200 name = "Smolyak basis (";
201 for (ordinal_type i=0; i<d-1; i++)
202 name += bases[i]->getName() + ", ";
203 name += bases[d-1]->getName() + ")";
204
205 // Allocate array for basis evaluation
206 basis_eval_tmp.resize(d);
207 for (ordinal_type j=0; j<d; j++)
208 basis_eval_tmp[j].resize(max_orders[j]+1);
209}
210
211template <typename ordinal_type, typename value_type, typename ordering_type>
214{
215}
216
217template <typename ordinal_type, typename value_type, typename ordering_type>
218ordinal_type
220order() const
221{
222 return p;
223}
224
225template <typename ordinal_type, typename value_type, typename ordering_type>
226ordinal_type
228dimension() const
229{
230 return d;
231}
232
233template <typename ordinal_type, typename value_type, typename ordering_type>
234ordinal_type
236size() const
237{
238 return sz;
239}
240
241template <typename ordinal_type, typename value_type, typename ordering_type>
242const Teuchos::Array<value_type>&
244norm_squared() const
245{
246 return norms;
247}
248
249template <typename ordinal_type, typename value_type, typename ordering_type>
250const value_type&
252norm_squared(ordinal_type i) const
253{
254 return norms[i];
255}
256
257template <typename ordinal_type, typename value_type, typename ordering_type>
258Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
261{
262#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
263 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Total Triple-Product Tensor Fill Time");
264#endif
265
267 bases, basis_set, basis_map, sm_pred, sm_pred, sparse_tol);
268}
269
270template <typename ordinal_type, typename value_type, typename ordering_type>
271Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
274{
275#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
276 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Total Triple-Product Tensor Fill Time");
277#endif
278
280 for (ordinal_type i=0; i<sm_pred.tp_preds.size(); ++i) {
281 k_pred.tp_preds.push_back(
282 TotalOrderPredicate<ordinal_type>(1, sm_pred.tp_preds[i].orders) );
283 }
284
286 bases, basis_set, basis_map, sm_pred, k_pred, sparse_tol);
287}
288
289template <typename ordinal_type, typename value_type, typename ordering_type>
290value_type
292evaluateZero(ordinal_type i) const
293{
294 // z = psi_{i_1}(0) * ... * psi_{i_d}(0) where i_1,...,i_d are the basis
295 // terms for coefficient i
296 value_type z = value_type(1.0);
297 for (ordinal_type j=0; j<d; j++)
298 z = z * bases[j]->evaluate(value_type(0.0), basis_map[i][j]);
299
300 return z;
301}
302
303template <typename ordinal_type, typename value_type, typename ordering_type>
304void
306evaluateBases(const Teuchos::ArrayView<const value_type>& point,
307 Teuchos::Array<value_type>& basis_vals) const
308{
309 for (ordinal_type j=0; j<d; j++)
310 bases[j]->evaluateBases(point[j], basis_eval_tmp[j]);
311
312 // Only evaluate basis upto number of terms included in basis_pts
313 for (ordinal_type i=0; i<sz; i++) {
314 value_type t = value_type(1.0);
315 for (ordinal_type j=0; j<d; j++)
316 t *= basis_eval_tmp[j][basis_map[i][j]];
317 basis_vals[i] = t;
318 }
319}
320
321template <typename ordinal_type, typename value_type, typename ordering_type>
322void
324print(std::ostream& os) const
325{
326 os << "Smolyak basis of order " << p << ", dimension " << d
327 << ", and size " << sz << ". Component bases:\n";
328 for (ordinal_type i=0; i<d; i++)
329 os << *bases[i];
330 os << "Basis vector norms (squared):\n\t";
331 for (ordinal_type i=0; i<static_cast<ordinal_type>(norms.size()); i++)
332 os << norms[i] << " ";
333 os << "\n";
334}
335
336template <typename ordinal_type, typename value_type, typename ordering_type>
339term(ordinal_type i) const
340{
341 return basis_map[i];
342}
343
344template <typename ordinal_type, typename value_type, typename ordering_type>
345ordinal_type
348{
349 typename coeff_set_type::const_iterator it = basis_set.find(term);
350 TEUCHOS_TEST_FOR_EXCEPTION(it == basis_set.end(), std::logic_error,
351 "Invalid term " << term);
352 return it->second;
353}
354
355template <typename ordinal_type, typename value_type, typename ordering_type>
356const std::string&
358getName() const
359{
360 return name;
361}
362
363template <typename ordinal_type, typename value_type, typename ordering_type>
364Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type, value_type> > >
366getCoordinateBases() const
367{
368 return bases;
369}
370
371template <typename ordinal_type, typename value_type, typename ordering_type>
374getMaxOrders() const
375{
376 return max_orders;
377}
expr val()
A comparison functor implementing a strict weak ordering based lexographic ordering.
Abstract base class for 1-D orthogonal polynomials.
static Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor(const Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > &bases, const basis_set_type &basis_set, const basis_map_type &basis_map, const coeff_predicate_type &coeff_pred, const k_coeff_predicate_type &k_coeff_pred, const value_type sparse_tol=1.0e-12)
Teuchos::Array< Teuchos::RCP< tensor_product_basis_type > > tp_bases
Tensor product bases comprising Smolyak set.
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 const Teuchos::Array< value_type > & norm_squared() const
Return array storing norm-squared of each basis polynomial.
ordinal_type dimension() const
Return dimension of basis.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeLinearTripleProductTensor() const
Compute linear triple product tensor where k = 0,1,..,d.
virtual ordinal_type index(const MultiIndex< ordinal_type > &term) const
Get index of the multivariate polynomial given orders of each coordinate.
Teuchos::Array< Teuchos::Array< value_type > > basis_eval_tmp
Temporary array used in basis evaluation.
Teuchos::Array< ordinal_type > smolyak_coeffs
Smolyak coefficients.
TensorProductBasis< ordinal_type, value_type, LexographicLess< coeff_type > > tensor_product_basis_type
std::string name
Name of basis.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
ordinal_type d
Total dimension of basis.
ordinal_type order() const
Return order of basis.
coeff_map_type basis_map
Basis map.
virtual ordinal_type size() const
Return total size of basis.
coeff_type max_orders
Maximum orders for each dimension.
Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > bases
Array of bases.
ordinal_type sz
Total size of basis.
value_type sparse_tol
Tolerance for computing sparse Cijk.
coeff_set_type basis_set
Basis set.
virtual const std::string & getName() const
Return string name of basis.
Teuchos::Array< value_type > norms
Norms.
virtual MultiIndex< ordinal_type > getMaxOrders() const
Return maximum order allowable for each coordinate basis.
SmolyakPredicate< TensorProductPredicate< ordinal_type > > sm_pred
Predicate for building sparse triple products.
Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > getCoordinateBases() const
Return coordinate bases.
SmolyakBasis(const Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > &bases, const index_set_type &index_set, const value_type &sparse_tol=1.0e-12, const coeff_compare_type &coeff_compare=coeff_compare_type())
Constructor.
virtual const MultiIndex< ordinal_type > & term(ordinal_type i) const
Get orders of each coordinate polynomial given an index i.
virtual void print(std::ostream &os) const
Print basis to stream os.
virtual ~SmolyakBasis()
Destructor.
ordinal_type p
Total order of basis.
virtual value_type evaluateZero(ordinal_type i) const
Evaluate basis polynomial i at zero.
Predicate functor for building sparse triple products.
Teuchos::Array< tp_predicate_type > tp_preds
Predicate functor for building sparse triple products.
Predicate functor for building sparse triple products based on total order.