44 #include "Teuchos_UnitTestHarness.hpp" 45 #include "Teuchos_TestingHelpers.hpp" 46 #include "Teuchos_UnitTestRepository.hpp" 47 #include "Teuchos_GlobalMPISession.hpp" 48 #include "Teuchos_ArrayView.hpp" 58 template <
typename OrdinalType,
typename ValueType>
77 template <
typename ordinal_type>
86 template <
typename scalar_type>
89 for (
int i=0; i<
x.size(); i++)
95 template <
typename scalar_type>
98 for (
int i=0; i<
x.size(); i++)
112 out <<
"For n =" << n <<
", k = " << k <<
", n_choose_k = " << v1
113 <<
" != " << v2 << std::endl;
121 out <<
"For n =" << n <<
", k = " << k <<
", n_choose_k = " << v1
122 <<
" != " << v2 << std::endl;
133 typedef index_set_type::multiindex_type multiindex_type;
135 typedef std::set<multiindex_type, less_type> multiindex_set;
136 typedef multiindex_set::iterator iterator;
138 multiindex_set sortedIndexSet(indexSet.begin(), indexSet.end());
141 std::ostream_iterator<multiindex_type> out_iterator(out,
"\n");
142 out << std::endl <<
"Sorted total order index set (dimension = " <<
setup.d
143 <<
", order = " <<
setup.p <<
"):" << std::endl;
144 std::copy(sortedIndexSet.begin(), sortedIndexSet.end(), out_iterator);
147 iterator prev = sortedIndexSet.begin();
148 iterator curr = prev; ++curr;
149 while (curr != sortedIndexSet.end()) {
153 while (i <
setup.d && order_prev == order_curr) {
154 order_prev -= (*prev)[i];
155 order_curr -= (*curr)[i];
158 if (order_prev >= order_curr) {
159 out <<
"previous index " << *prev <<
" and current index " 160 << *curr <<
" are out of order" << std::endl;
173 typedef index_set_type::multiindex_type multiindex_type;
175 typedef std::set<multiindex_type, less_type> multiindex_set;
176 typedef multiindex_set::iterator iterator;
178 multiindex_set sortedIndexSet(indexSet.begin(), indexSet.end());
181 std::ostream_iterator<multiindex_type> out_iterator(out,
"\n");
182 out << std::endl <<
"Sorted lexicographic index set (dimension = " 183 <<
setup.d <<
", order = " <<
setup.p <<
"):" << std::endl;
184 std::copy(sortedIndexSet.begin(), sortedIndexSet.end(), out_iterator);
187 iterator prev = sortedIndexSet.begin();
188 iterator curr = prev; ++curr;
189 while (curr != sortedIndexSet.end()) {
191 while (i <
setup.d && (*prev)[i] == (*curr)[i]) ++i;
192 if (i ==
setup.d || (*prev)[i] >= (*curr)[i]) {
193 out <<
"previous index " << *prev <<
" and current index " 194 << *curr <<
" are out of order" << std::endl;
208 TEUCHOS_TEST_EQUALITY(less(-0.774597,-0.774597),
false, out, success);
209 TEUCHOS_TEST_EQUALITY(less(-0.774597+tol/2.0,-0.774597),
false, out, success);
210 TEUCHOS_TEST_EQUALITY(less(-0.774597-tol/2.0,-0.774597),
false, out, success);
211 TEUCHOS_TEST_EQUALITY(less(-0.774597,-0.774597+tol/2.0),
false, out, success);
212 TEUCHOS_TEST_EQUALITY(less(-0.774597,-0.774597-tol/2.0),
false, out, success);
213 TEUCHOS_TEST_EQUALITY(less(-0.774597,0.0),
true, out, success);
214 TEUCHOS_TEST_EQUALITY(less(0.0,-0.774597),
false, out, success);
222 term_type a(2), b(2);
224 a[0] = -0.774597; a[1] = -0.774597;
225 b[0] = -0.774597; b[1] = 0.0;
227 TEUCHOS_TEST_EQUALITY(less(a,b),
true, out, success);
228 TEUCHOS_TEST_EQUALITY(less(b,a),
false, out, success);
236 typedef index_set_type::multiindex_type multiindex_type;
237 typedef index_set_type::iterator iterator;
241 out << std::endl <<
"Total order index set (dimension = " <<
setup.d
242 <<
", order = " <<
setup.p <<
"):" << std::endl;
243 std::ostream_iterator<multiindex_type> out_iterator(out,
"\n");
244 std::copy(indexSet.begin(), indexSet.end(), out_iterator);
247 for (iterator i=indexSet.begin(); i!=indexSet.end(); ++i) {
248 if (i->order() < 0 || i->order() >
setup.p) {
249 out <<
"index " << *i <<
" does not lie in total order set! " 258 typedef std::set<multiindex_type, less_type> multiindex_set;
259 multiindex_set sortedIndexSet(indexSet.begin(), indexSet.end());
261 out <<
"sorted index set size = " << sortedIndexSet.size() << std::endl;
262 out <<
"expected index set size = " 264 if (static_cast<ordinal_type>(sortedIndexSet.size()) !=
270 AnisotropicTotalOrderIndexSet ) {
275 typedef index_set_type::multiindex_type multiindex_type;
276 typedef index_set_type::iterator iterator;
277 multiindex_type upper(
setup.d);
280 index_set_type indexSet(
setup.p, upper);
283 out << std::endl <<
"Anisotropic total order index set (dimension = " 284 <<
setup.d <<
", order = " <<
setup.p <<
", and component orders = " 285 << upper <<
"):" << std::endl;
286 std::ostream_iterator<multiindex_type> out_iterator(out,
"\n");
287 std::copy(indexSet.begin(), indexSet.end(), out_iterator);
290 for (iterator i=indexSet.begin(); i!=indexSet.end(); ++i) {
291 if (i->order() < 0 || i->order() >
setup.p || !i->termWiseLEQ(upper)) {
292 out <<
"index " << *i <<
" does not lie in total order set! " 325 typedef std::map<coeff_type, ordinal_type, less_type> basis_set_type;
326 typedef Teuchos::Array<coeff_type> basis_map_type;
327 basis_set_type basis_set;
328 basis_map_type basis_map;
330 indexSet, basis_set, basis_map);
334 Teuchos::Array< Stokhos::MultiIndex<ordinal_type> > terms;
335 Teuchos::Array<ordinal_type> num_terms;
340 TEUCHOS_TEST_EQUALITY(static_cast<ordinal_type>(basis_set.size()),
341 static_cast<ordinal_type>(basis_map.size()),
343 TEUCHOS_TEST_EQUALITY(static_cast<ordinal_type>(basis_set.size()),
344 static_cast<ordinal_type>(terms.size()),
346 TEUCHOS_TEST_EQUALITY(sz, static_cast<ordinal_type>(terms.size()),
349 std::ostream_iterator<ordinal_type> out_iterator(out,
" ");
353 out <<
"term " << basis_map[i] <<
" == " << terms[i] <<
" : ";
354 bool is_equal =
true;
356 is_equal = is_equal && terms[i][
j] == basis_map[i][
j];
358 out <<
"passed" << std::endl;
360 out <<
"failed" << std::endl;
365 TEUCHOS_TEST_EQUALITY(basis_set[basis_map[i]], i, out, success);
380 typedef std::map<coeff_type, ordinal_type, less_type> basis_set_type;
381 typedef Teuchos::Array<coeff_type> basis_map_type;
382 basis_set_type basis_set;
383 basis_map_type basis_map;
385 indexSet, basis_set, basis_map);
393 TEUCHOS_TEST_EQUALITY(static_cast<ordinal_type>(basis_set.size()),
394 static_cast<ordinal_type>(basis_map.size()),
396 TEUCHOS_TEST_EQUALITY(sz, static_cast<ordinal_type>(basis_set.size()),
399 std::ostream_iterator<ordinal_type> out_iterator(out,
" ");
403 out <<
"term " << basis_map[i] <<
" <= " <<
setup.p <<
" : ";
406 is_less = is_less && basis_map[i][
j] <=
setup.p;
408 out <<
"passed" << std::endl;
410 out <<
"failed" << std::endl;
415 TEUCHOS_TEST_EQUALITY(basis_set[basis_map[i]], i, out, success);
420 template <
typename ordinal_type>
427 template <
typename term_type>
437 template <
typename basis_set_type>
444 template <
typename term_type>
459 index_set_type indexSet(dim, 0, order);
464 typedef std::map<coeff_type, ordinal_type, less_type> basis_set_type;
465 typedef Teuchos::Array<coeff_type> basis_map_type;
466 basis_set_type basis_set;
467 basis_map_type basis_map;
469 indexSet, basis_set, basis_map);
472 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> > > bases(dim);
480 Teuchos::RCP<Cijk_type> Cijk =
481 Stokhos::ProductBasisUtils::computeTripleProductTensor<ordinal_type,value_type>(bases, basis_set, basis_map, pred, pred);
485 Teuchos::RCP<Cijk_type> Cijk2 =
486 basis->computeTripleProductTensor();
489 TEUCHOS_TEST_EQUALITY(Cijk->num_k(), Cijk2->num_k(), out, success);
490 TEUCHOS_TEST_EQUALITY(Cijk->num_entries(), Cijk2->num_entries(), out, success);
493 for (Cijk_type::k_iterator k_it=Cijk2->k_begin();
494 k_it!=Cijk2->k_end(); ++k_it) {
495 int k = Stokhos::index(k_it);
496 for (Cijk_type::kj_iterator j_it = Cijk2->j_begin(k_it);
497 j_it != Cijk2->j_end(k_it); ++j_it) {
498 int j = Stokhos::index(j_it);
499 for (Cijk_type::kji_iterator i_it = Cijk2->i_begin(j_it);
500 i_it != Cijk2->i_end(j_it); ++i_it) {
501 int i = Stokhos::index(i_it);
502 double c = Cijk->getValue(i,
j,k);
503 double c2 = Stokhos::value(i_it);
509 <<
"Check: rel_err( C(" << i <<
"," <<
j <<
"," << k <<
") )" 510 <<
" = " <<
"rel_err( " << c <<
", " << c2 <<
" ) = " << err
511 <<
" <= " << tol <<
" : ";
512 if (s) out <<
"Passed.";
513 else out <<
"Failed!";
516 success = success && s;
531 index_set_type indexSet(dim, 0, order);
536 typedef std::map<coeff_type, ordinal_type, less_type> basis_set_type;
537 typedef Teuchos::Array<coeff_type> basis_map_type;
538 basis_set_type basis_set;
539 basis_map_type basis_map;
541 indexSet, basis_set, basis_map);
544 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> > > bases(dim);
552 Teuchos::RCP<Cijk_type> Cijk =
553 Stokhos::ProductBasisUtils::computeTripleProductTensorNew<ordinal_type,value_type>(bases, basis_set, basis_map, pred, pred,
false);
557 Teuchos::RCP<Cijk_type> Cijk2 =
558 basis->computeTripleProductTensor();
561 TEUCHOS_TEST_EQUALITY(Cijk->num_k(), Cijk2->num_k(), out, success);
562 TEUCHOS_TEST_EQUALITY(Cijk->num_entries(), Cijk2->num_entries(), out, success);
565 for (Cijk_type::k_iterator k_it=Cijk2->k_begin();
566 k_it!=Cijk2->k_end(); ++k_it) {
567 int k = Stokhos::index(k_it);
568 for (Cijk_type::kj_iterator j_it = Cijk2->j_begin(k_it);
569 j_it != Cijk2->j_end(k_it); ++j_it) {
570 int j = Stokhos::index(j_it);
571 for (Cijk_type::kji_iterator i_it = Cijk2->i_begin(j_it);
572 i_it != Cijk2->i_end(j_it); ++i_it) {
573 int i = Stokhos::index(i_it);
574 double c = Cijk->getValue(i,
j,k);
575 double c2 = Stokhos::value(i_it);
581 <<
"Check: rel_err( C(" << i <<
"," <<
j <<
"," << k <<
") )" 582 <<
" = " <<
"rel_err( " << c <<
", " << c2 <<
" ) = " << err
583 <<
" <= " << tol <<
" : ";
584 if (s) out <<
"Passed.";
585 else out <<
"Failed!";
588 success = success && s;
602 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> > > bases(dim);
614 Teuchos::RCP<Cijk_type> Cijk =
618 Teuchos::RCP<Cijk_type> Cijk2 =
622 TEUCHOS_TEST_EQUALITY(Cijk->num_k(), Cijk2->num_k(), out, success);
623 TEUCHOS_TEST_EQUALITY(Cijk->num_entries(), Cijk2->num_entries(), out, success);
626 for (Cijk_type::k_iterator k_it=Cijk2->k_begin();
627 k_it!=Cijk2->k_end(); ++k_it) {
628 int k = Stokhos::index(k_it);
629 for (Cijk_type::kj_iterator j_it = Cijk2->j_begin(k_it);
630 j_it != Cijk2->j_end(k_it); ++j_it) {
631 int j = Stokhos::index(j_it);
632 for (Cijk_type::kji_iterator i_it = Cijk2->i_begin(j_it);
633 i_it != Cijk2->i_end(j_it); ++i_it) {
634 int i = Stokhos::index(i_it);
635 double c = Cijk->getValue(i,
j,k);
636 double c2 = Stokhos::value(i_it);
642 <<
"Check: rel_err( C(" << i <<
"," <<
j <<
"," << k <<
") )" 643 <<
" = " <<
"rel_err( " << c <<
", " << c2 <<
" ) = " << err
644 <<
" <= " << tol <<
" : ";
645 if (s) out <<
"Passed.";
646 else out <<
"Failed!";
649 success = success && s;
660 typedef index_set_type::multiindex_type multiindex_type;
662 typedef std::set<multiindex_type, less_type> multiindex_set;
663 typedef multiindex_set::iterator iterator;
665 multiindex_set sortedIndexSet(indexSet.begin(), indexSet.end());
669 iterator i = sortedIndexSet.begin();
671 while (i != sortedIndexSet.end()) {
673 out << *i <<
": index = " << idx <<
" mapped index = " << idx_mapping
675 if (idx == idx_mapping)
692 typedef index_set_type::multiindex_type multiindex_type;
694 typedef std::set<multiindex_type, less_type> multiindex_set;
695 typedef multiindex_set::iterator iterator;
697 multiindex_set sortedIndexSet(indexSet.begin(), indexSet.end());
701 iterator i = sortedIndexSet.begin();
703 while (i != sortedIndexSet.end()) {
705 out << *i <<
": index = " << idx <<
" mapped index = " << idx_mapping
707 if (idx == idx_mapping)
726 typedef index_set_type::multiindex_type multiindex_type;
728 typedef std::set<multiindex_type, less_type> multiindex_set;
729 index_set_type indexSet(d, 0, p);
730 multiindex_set sortedIndexSet(indexSet.begin(), indexSet.end());
734 multiindex_type index(d), num_terms(d), orders(d,p);
743 while (dim > 0 && index[dim] > orders[dim]) {
749 orders[i] = orders[i-1] - index[i-1];
751 if (index[dim] > orders[dim])
765 out << index <<
": index = " << idx <<
" mapped index = " << I
782 Teuchos::GlobalMPISession mpiSession(&argc, &
argv);
783 int res = Teuchos::UnitTestRepository::runUnitTestsFromMain(argc,
argv);
784 Teuchos::TimeMonitor::summarize(std::cout);
A functor for comparing floating-point numbers to some tolerance.
ordinal_type factorial(const ordinal_type &n)
const basis_set_type & basis_set
ordinal_type n_choose_k(const ordinal_type &n, const ordinal_type &k)
Compute bionomial coefficient (n ; k) = n!/( k! (n-k)! )
Teuchos::RCP< Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensorLTO(const TotalOrderBasis< ordinal_type, value_type, LexographicLess< MultiIndex< ordinal_type > > > &product_basis, bool symmetric=false)
Multivariate orthogonal polynomial basis generated from a total order tensor product of univariate po...
Container storing a term in a generalized tensor product.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
bool operator()(const term_type &term) const
TEUCHOS_UNIT_TEST(Stokhos_ProductBasisUtils, NChooseK)
static void buildProductBasis(const index_set_type &index_set, const growth_rule_type &growth_rule, basis_set_type &basis_set, basis_map_type &basis_map)
Generate a product basis from an index set.
std::enable_if< Kokkos::is_view_uq_pce< Kokkos::View< RD, RP... > >::value &&Kokkos::is_view_uq_pce< Kokkos::View< XD, XP... > >::value >::type sum(const Kokkos::View< RD, RP... > &r, const Kokkos::View< XD, XP... > &x)
int main(int argc, char *argv[])
A comparison functor implementing a strict weak ordering based total-order ordering, recursive on the dimension.
static ordinal_type compute_terms(ordinal_type p, ordinal_type d, ordinal_type &sz, Teuchos::Array< MultiIndex< ordinal_type > > &terms, Teuchos::Array< ordinal_type > &num_terms)
Compute the 2-D array of basis terms which maps a basis index into the orders for each basis dimensio...
UnitTestSetup< ordinal_type, value_type > setup
total_order_predicate(ordinal_type dim_, ordinal_type order_)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
scalar_type quad_func2(const Teuchos::Array< scalar_type > &x)
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
general_predicate(const basis_set_type &basis_set_)
An anisotropic total order index set.
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
KOKKOS_INLINE_FUNCTION PCE< Storage > exp(const PCE< Storage > &a)
Legendre polynomial basis.
Stokhos::Sparse3Tensor< int, double > Cijk_type
scalar_type quad_func1(const Teuchos::Array< scalar_type > &x)
bool operator()(const term_type &term) const
An isotropic total order index set.
A comparison functor implementing a strict weak ordering based lexographic ordering.
KOKKOS_INLINE_FUNCTION PCE< Storage > sin(const PCE< Storage > &a)
ordinal_type lexicographicMapping(const Stokhos::MultiIndex< ordinal_type > &index, ordinal_type max_order)
A tensor product index set.
ordinal_type totalOrderMapping(const Stokhos::MultiIndex< ordinal_type > &index)