From d4808392cc17b2219252a6371f4120b66447cf97 Mon Sep 17 00:00:00 2001 From: Miguel Raggi Date: Fri, 5 Jan 2018 19:52:00 -0600 Subject: [PATCH] Multisets is now a random access container! --- Benchmarks/benchmarker.hpp | 14 ++++ Benchmarks/benchtable.hpp | 10 +++ Benchmarks/combinations_benchmark.hpp | 18 ---- Benchmarks/combinations_tree_benchmark.hpp | 11 --- Benchmarks/main.cpp | 98 +++++++++++++--------- Benchmarks/permutations_benchmark.hpp | 11 --- README.md | 3 + examples/multisets.cpp | 9 +- include/Discreture/Combinations.hpp | 4 +- include/Discreture/Misc.hpp | 4 + include/Discreture/Multisets.hpp | 76 ++++++++++++++++- include/Discreture/NaturalNumber.hpp | 12 ++- include/Discreture/NumberRange.hpp | 4 +- tests/multiset_tests.cpp | 49 ++++++++++- 14 files changed, 232 insertions(+), 91 deletions(-) diff --git a/Benchmarks/benchmarker.hpp b/Benchmarks/benchmarker.hpp index 3aacf63..c047566 100644 --- a/Benchmarks/benchmarker.hpp +++ b/Benchmarks/benchmarker.hpp @@ -57,3 +57,17 @@ double ForEachBenchmark(const Container& A) }); }); } + +template +double ConstructionBenchmark(const Container& A, int numtimes) +{ + return Benchmark([&A,numtimes]() + { + for (int i = 0; i < numtimes; ++i) + { + auto t = dscr::random::random_int(0,A.size()); + DoNotOptimize(A[t]); + } + }); +} + diff --git a/Benchmarks/benchtable.hpp b/Benchmarks/benchtable.hpp index ea9b635..d82ce70 100644 --- a/Benchmarks/benchtable.hpp +++ b/Benchmarks/benchtable.hpp @@ -144,3 +144,13 @@ BenchRow ProduceRowForEach(std::string name, const Container& A) return BenchRow(name, t, A.size()); } + +template +BenchRow ProduceRowConstruct(std::string name, const Container& A, int numtimes) +{ + double t = ConstructionBenchmark(A,numtimes); + + name += " Construct"; + + return BenchRow(name, t, numtimes); +} diff --git a/Benchmarks/combinations_benchmark.hpp b/Benchmarks/combinations_benchmark.hpp index 1730c5a..cdfe804 100644 --- a/Benchmarks/combinations_benchmark.hpp +++ b/Benchmarks/combinations_benchmark.hpp @@ -4,17 +4,6 @@ #include "do_not_optimize.hpp" #include "Probability.hpp" -inline void BM_CombinationsConstruct(int n, int k, int numtimes) -{ - dscr::combinations_fast X(n, k); - - for (int i = 0; i < numtimes; ++i) - { - auto t = dscr::random::random_int(0,X.size()); - DoNotOptimize(X[t]); - } -} - inline void BM_CombinationsNAP(int n, int k) { dscr::combinations::combination comb(k); @@ -41,13 +30,6 @@ inline void BM_combinationsIterator(int n, int k, long size) } } -inline void BM_combinations(int n, int k, long size) -{ - for (auto& t : dscr::combinations(n,k)) - { - DoNotOptimize(t); - } -} class algoT_iterator : public boost::iterator_facade< algoT_iterator, diff --git a/Benchmarks/combinations_tree_benchmark.hpp b/Benchmarks/combinations_tree_benchmark.hpp index abbe52c..7eeb030 100644 --- a/Benchmarks/combinations_tree_benchmark.hpp +++ b/Benchmarks/combinations_tree_benchmark.hpp @@ -37,17 +37,6 @@ inline long BM_CombinationsTreeFindAll(int n, int k) return size; } -inline void BM_CombinationsTreeConstruct(int n, int k, int numtimes) -{ - dscr::combinations_tree_fast X(n, k); - - for (int i = 0; i < numtimes; ++i) - { - auto t = dscr::random::random_int(0,X.size()); - DoNotOptimize(X[t]); - } -} - #include "external/euler314_combination_iterator.hpp" inline void BM_CombinationsTreeEuler314(int n,int k) diff --git a/Benchmarks/main.cpp b/Benchmarks/main.cpp index 3c53b74..1cd6a51 100644 --- a/Benchmarks/main.cpp +++ b/Benchmarks/main.cpp @@ -15,7 +15,7 @@ int main() using dscr::binomial; std::ios_base::sync_with_stdio(false); - dscr::Chronometer C; + dscr::Chronometer chrono; cout << "|============================== Starting Speed Tests =============================|" << endl; @@ -23,7 +23,7 @@ int main() const int n = 40; const int k = 10; const auto combs_size = binomial(n,k); - const int construct = 100000; + const int construct = 1000000; const int nperm = 12; @@ -33,65 +33,85 @@ int main() const int nmotzkin = 20; const int nmultiset = 19; - //fast tests -// const int n = 25; -// const int k = 12; -// const int combconstruct = 1000; -// -// const int nperm = 9; -// -// const int npart = 20; -// const int nsetpart = 7; -// const int ndyck = 10; -// const int nmotzkin = 10; - + dscr::combinations C(n,k); + dscr::combinations_fast CF(n,k); + dscr::combinations_tree CT(n,k); + dscr::combinations_tree_fast CTF(n,k); + + dscr::permutations P(nperm); + dscr::permutations_fast PF(nperm); + + dscr::dyck_paths DP(ndyck); + dscr::dyck_paths_fast DPF(ndyck); + dscr::motzkin_paths MP(ndyck); + dscr::motzkin_paths_fast MPF(ndyck); + + dscr::partitions PT(npart); + dscr::partitions_fast PTF(npart); + dscr::set_partitions SPT(nsetpart); + + dscr::multisets MS({4,2,3,1,0,1,5,0,5,4,0,1,1,5,2,0,2}); + dscr::multisets_fast MSF({4,2,3,1,0,1,5,0,5,4,0,1,1,5,2,0,2}); + + BenchRow::print_header(cout); BenchRow::print_line(cout); - cout << ProduceRowForEach("Combinations", dscr::combinations(n,k)); - cout << ProduceRowForEach("Combinations Stack", dscr::combinations_fast(n,k)); - cout << BenchRow("Combinations (No iterator)", Benchmark([](){BM_CombinationsNAP(n,k);}), combs_size); - cout << ProduceRowForward("Combinations", dscr::combinations(n,k)); - cout << ProduceRowForward("Combinations Stack", dscr::combinations_fast(n,k)); - cout << ProduceRowReverse("Combinations", dscr::combinations(n,k)); - cout << BenchRow("Combinations Construct", Benchmark([](){BM_CombinationsConstruct(n,k,construct);}), construct); + cout << ProduceRowForEach("Combinations", C); + cout << ProduceRowForEach("Combinations Stack", CF); + cout << ProduceRowForward("Combinations", C); + cout << ProduceRowForward("Combinations Stack", CF); + cout << ProduceRowReverse("Combinations", C); + cout << ProduceRowReverse("Combinations Stack", CF); + cout << ProduceRowConstruct("Combinations", C, construct); + cout << ProduceRowConstruct("Combinations Stack", CF, construct); BenchRow::print_line(cout); - cout << ProduceRowForEach("Combinations Tree", dscr::combinations_tree(n,k)); - cout << ProduceRowForEach("Combinations Tree Stack", dscr::combinations_tree_fast(n,k)); - cout << BenchRow("Combinations Tree (No iterator)", Benchmark([](){BM_CombinationsTreeNAP(n,k);}), combs_size); - cout << ProduceRowForward("Combinations Tree", dscr::combinations_tree(n,k)); - cout << ProduceRowForward("Combinations Tree Stack", dscr::combinations_tree_fast(n,k)); - cout << ProduceRowReverse("Combinations Tree", dscr::combinations_tree(n,k)); + cout << ProduceRowForEach("Combinations Tree", CT); + cout << ProduceRowForEach("Combinations Tree Stack", CTF); + cout << ProduceRowForward("Combinations Tree", CT); + cout << ProduceRowForward("Combinations Tree Stack", CTF); + cout << ProduceRowReverse("Combinations Tree", CT); + cout << ProduceRowReverse("Combinations Tree Stack", CTF); #ifdef TEST_GSL_COMBINATIONS cout << BenchRow("Combinations Tree GSL", Benchmark([](){BM_CombinationsTreeGSL(n,k);}), combs_size); #endif - cout << BenchRow("Combinations Tree Construct", Benchmark([](){BM_CombinationsTreeConstruct(n,k,construct);}), construct); + cout << ProduceRowConstruct("Combinations Tree", CT, construct); + cout << ProduceRowConstruct("Combinations Tree Stack", CTF, construct); BenchRow::print_line(cout); - cout << ProduceRowForward("Permutations", dscr::permutations_fast(nperm)); - cout << ProduceRowReverse("Permutations", dscr::permutations_fast(nperm)); - cout << BenchRow("Permutations Construct", Benchmark([](){BM_PermutationsConstruct(nperm,construct);}), construct); + cout << ProduceRowForward("Permutations", P); + cout << ProduceRowForward("Permutations Stack", PF); + cout << ProduceRowReverse("Permutations", P); + cout << ProduceRowReverse("Permutations Stack", PF); + cout << ProduceRowConstruct("Permutations", P, construct); + cout << ProduceRowConstruct("Permutations Stack", PF, construct); BenchRow::print_line(cout); - cout << ProduceRowForward("Multisets", dscr::multisets_fast({4,2,3,1,0,1,5,0,5,4,0,1,1,5,2,0,2})); - cout << ProduceRowReverse("Multisets", dscr::multisets_fast({4,2,3,1,0,1,5,0,5,4,0,1,1,5,2,0,2})); - + cout << ProduceRowForward("Multisets", MS); + cout << ProduceRowForward("Multisets Stack", MSF); + cout << ProduceRowReverse("Multisets", MS); + cout << ProduceRowReverse("Multisets Stack", MSF); + cout << ProduceRowConstruct("Multisets", MS, construct); + cout << ProduceRowConstruct("Multisets", MSF, construct); BenchRow::print_line(cout); - cout << ProduceRowForward("Dyck Paths", dscr::dyck_paths_fast(ndyck)); + cout << ProduceRowForward("Dyck Paths", DP); + cout << ProduceRowForward("Dyck Paths Stack", DPF); BenchRow::print_line(cout); - cout << ProduceRowForward("Motzkin Paths", dscr::motzkin_paths_fast(nmotzkin)); + cout << ProduceRowForward("Motzkin Paths", MP); + cout << ProduceRowForward("Motzkin Paths Stack", MPF); BenchRow::print_line(cout); - cout << ProduceRowForward("Partitions", dscr::partitions_fast(npart)); + cout << ProduceRowForward("Partitions", PT); + cout << ProduceRowForward("Partitions Stack", PTF); BenchRow::print_line(cout); - cout << ProduceRowForward("Set Partitions", dscr::set_partitions(nsetpart)); + cout << ProduceRowForward("Set Partitions", SPT); BenchRow::print_line(cout); - cout << "\nTotal Time taken = " << C.Reset() << "s" << endl; + cout << "\nTotal Time taken = " << chrono.Reset() << "s" << endl; return 0; } diff --git a/Benchmarks/permutations_benchmark.hpp b/Benchmarks/permutations_benchmark.hpp index f73b6f9..94a499b 100644 --- a/Benchmarks/permutations_benchmark.hpp +++ b/Benchmarks/permutations_benchmark.hpp @@ -4,17 +4,6 @@ #include "do_not_optimize.hpp" #include "Probability.hpp" -inline void BM_PermutationsConstruct(int n, int numtimes) -{ - dscr::permutations X(n); - - for (int i = 0; i < numtimes; ++i) - { - auto t = dscr::random::random_int(0,X.size()); - DoNotOptimize(X[t]); - } -} - inline void BM_PermutationsRandom(int n, int numtimes) { dscr::permutations X(n); diff --git a/README.md b/README.md index 13ee8c7..d288943 100644 --- a/README.md +++ b/README.md @@ -388,3 +388,6 @@ The important columns are Speed and Speed (w/o _fast). Higher is better. They me # Contributing Please help us testing, debugging, benchmarking, packaging for the various distros, etc. Also, if you use discreture for your own purposes, let us know! + +Here is what we hope to accomplish: + diff --git a/examples/multisets.cpp b/examples/multisets.cpp index 16846b2..268b904 100644 --- a/examples/multisets.cpp +++ b/examples/multisets.cpp @@ -19,12 +19,17 @@ int main() cout << "These are all the submultisets of " << total << endl; auto X = multisets(total); - + int i = 0; for (auto& x : X) { cout << x << endl; +// cout << x << " <---" << i << " == " << X.get_index(x) << endl; +// multisets::multiset u(total.size()); +// multisets::construct_multiset(u,total,i); +// cout << "Constructed: " << u << endl; +// ++i; } - +// return 0; cout << "And here they are in reverse order: " << endl; for (auto it = X.rbegin(); it != X.rend(); ++it) { diff --git a/include/Discreture/Combinations.hpp b/include/Discreture/Combinations.hpp index a646b21..8ffed07 100644 --- a/include/Discreture/Combinations.hpp +++ b/include/Discreture/Combinations.hpp @@ -9,6 +9,7 @@ #include "CombinationsTree.hpp" #include "CombinationsTreePrunned.hpp" #include "combinations_bf.hpp" //Horrible. Do NOT read. Please. But I can't find another way. Sorry about that. If you think you can do better, please, tell me about it. +#include "NaturalNumber.hpp" namespace dscr { @@ -183,14 +184,13 @@ class basic_combinations { IntType r = k - i; IntType first = r; - + while (binomial(first, r) <= m) { ++first; } --first; - data[r - 1] = first; m -= binomial(first, r); } diff --git a/include/Discreture/Misc.hpp b/include/Discreture/Misc.hpp index 22e63a0..6965c05 100644 --- a/include/Discreture/Misc.hpp +++ b/include/Discreture/Misc.hpp @@ -1,11 +1,14 @@ #pragma once + #include #include #include + #include #include #include + namespace dscr { @@ -42,4 +45,5 @@ inline T pow(T a, unsigned long n) return r; } + } diff --git a/include/Discreture/Multisets.hpp b/include/Discreture/Multisets.hpp index eeb1e1d..a48ef0e 100644 --- a/include/Discreture/Multisets.hpp +++ b/include/Discreture/Multisets.hpp @@ -1,6 +1,7 @@ #pragma once #include "VectorHelpers.hpp" #include "Misc.hpp" +#include "NaturalNumber.hpp" #include namespace dscr @@ -14,9 +15,9 @@ class basic_multisets using value_type = RAContainerInt; using multiset = value_type; - static bool can_increment(size_t index, const multiset& sub, const multiset& total) + static void next_multiset(multiset& sub, const multiset& total) { - return sub[index] + 1 <= total[index]; + next_multiset(sub,total,total.size()); } static void next_multiset(multiset& sub, const multiset& total, size_t n) @@ -51,6 +52,47 @@ class basic_multisets } } + static void construct_multiset(multiset& sub, const multiset& total, size_type m) + { + assert(sub.size() == total.size()); + size_type n = total.size(); + if (n == 0) + return; + for (auto& s : sub) s = 0; + std::vector coeffs(n); + coeffs[0] = 1; + for (size_type i = 1; i < n; ++i) + { + coeffs[i] = coeffs[i-1]*(total[i-1]+1); + } + + for (long i = n-1; i >= 0; --i) + { + size_type w = coeffs[i]; + auto t = biggest_upto_satisfying_predicate(total[i], [m,w](size_type a) + { + return a*w <= m; + }); + sub[i] = t; + m -= w*t; + if (m <= 0) + break; + } + } + + size_type get_index(const multiset& sub) const + { + assert(sub.size() == m_total.size()); + size_type coeff = 1; + size_type result = 0; + for (size_t i = 0; i < m_total.size(); ++i) + { + result += coeff*sub[i]; + coeff *= (m_total[i]+1); + } + return result; + } + //////////////////////////////////////////////////////////// /// \brief class of all submultisets of a given multiset, expressed as incidence vectors with multiplicities /// \param IntType can be an int, short, etc. @@ -103,7 +145,7 @@ class basic_multisets class iterator : public boost::iterator_facade< iterator, const multiset&, - boost::bidirectional_traversal_tag + boost::random_access_traversal_tag > { @@ -115,6 +157,11 @@ class basic_multisets } + size_type ID() const + { + return m_ID; + } + private: iterator(const multiset& total, size_type id) : m_ID(id), m_total(total) {} @@ -143,6 +190,17 @@ class basic_multisets return m_ID == it.m_ID; } + void advance(difference_type m) + { + m_ID += m; + construct_multiset(m_submulti,m_total,m_ID); + } + + difference_type distance_to(const iterator& it) const + { + return it.m_ID - m_ID; + } + private: size_type m_ID{0}; size_type m_n{0}; @@ -163,6 +221,13 @@ class basic_multisets return iterator(m_total,size()); } + multiset operator[](size_type t) const + { + auto w = begin(); + w += t; + return *w; + } + class reverse_iterator : public boost::iterator_facade< reverse_iterator, @@ -228,6 +293,11 @@ class basic_multisets private: RAContainerInt m_total; size_type m_size; + + static bool can_increment(size_t index, const multiset& sub, const multiset& total) + { + return sub[index] < total[index]; + } }; using multisets = basic_multisets; using multisets_fast = basic_multisets>; diff --git a/include/Discreture/NaturalNumber.hpp b/include/Discreture/NaturalNumber.hpp index 3b7a974..96edbd4 100644 --- a/include/Discreture/NaturalNumber.hpp +++ b/include/Discreture/NaturalNumber.hpp @@ -45,7 +45,7 @@ class basic_natural_number { public: - iterator(IntType t) : m_ID(t) { } + explicit iterator(IntType t=0) : m_ID(t) { } private: void increment() @@ -111,4 +111,14 @@ class basic_natural_number using natural_number = basic_natural_number; +template +static IntType biggest_upto_satisfying_predicate(IntType up_to, Pred pred) +{ + basic_natural_number M(++up_to); + + auto t = std::partition_point(M.begin(), M.end(), pred); + + return *t-1; +} + } // end namespace dscr; diff --git a/include/Discreture/NumberRange.hpp b/include/Discreture/NumberRange.hpp index 91d8e75..9089b03 100644 --- a/include/Discreture/NumberRange.hpp +++ b/include/Discreture/NumberRange.hpp @@ -70,13 +70,13 @@ class basic_number_range { public: - iterator(size_type t_from, size_type t_step = 1) : m_step(t_step), m_ID(t_from) { } + explicit iterator(size_type t_from = 0, size_type t_step = 1) : m_step(t_step), m_ID(t_from) { } size_type step() const { return m_step; } - + private: void increment() { diff --git a/tests/multiset_tests.cpp b/tests/multiset_tests.cpp index 2e63eea..6cfed20 100644 --- a/tests/multiset_tests.cpp +++ b/tests/multiset_tests.cpp @@ -31,6 +31,23 @@ std::vector get_random_multiset(int n) return result; } +template +void dumb_advance(Iter& it, int m) +{ + while (m > 0) + { + ++it; + --m; + } + + while (m < 0) + { + --it; + ++m; + } + +} + TEST(Multisets,ForwardIteration) { for (int n = 0; n < 12; ++n) @@ -40,9 +57,15 @@ TEST(Multisets,ForwardIteration) set S(X.begin(),X.end()); ASSERT_EQ(X.size(), S.size()); //check if all multisets are different + long i = 0; for (const auto& x : X) { check_multiset(x, total); + ASSERT_EQ(X.get_index(x),i); +// cout << "all good: " << x << endl; +// cout << "Lets see: " << X[i] << endl; + ASSERT_EQ(X[i],x); + ++i; } } } @@ -63,7 +86,29 @@ TEST(Multisets,ReverseIteration) } } -TEST(Multisets,BidirectionalIteration) +TEST(Multisets,Bidirectional) +{ + for (int n = 0; n < 8; ++n) + { + auto total = get_random_multiset(n); + multisets X(total); + if (X.size() < 2) + continue; + auto it = X.begin(); + int start = dscr::random::random_int(0,X.size()/2); + dumb_advance(it,start); + auto x = *it; + ASSERT_EQ(*it,X[start]); + int adv = random::random_int(0,X.size()/2); + dumb_advance(it,adv); + ASSERT_EQ(*it,X[start+adv]); + dumb_advance(it,-adv); + + ASSERT_EQ(*it,x); + } +} + +TEST(Multisets,RandomAccess) { for (int n = 0; n < 8; ++n) { @@ -77,7 +122,7 @@ TEST(Multisets,BidirectionalIteration) auto x = *it; int adv = random::random_int(0,X.size()/2); std::advance(it,adv); - + ASSERT_EQ(*it,X[start+adv]); std::advance(it,-adv); ASSERT_EQ(*it,x);