[med-svn] [libhpptools] 02/03: Imported Upstream version 1.1.1
Andreas Tille
tille at debian.org
Thu Sep 15 08:02:24 UTC 2016
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository libhpptools.
commit f76ab1819e7240812e054a3711557b9adbe89320
Author: Andreas Tille <tille at debian.org>
Date: Thu Sep 15 08:52:32 2016 +0200
Imported Upstream version 1.1.1
---
README.org | 16 ++--
examples/benchmark-logsum.make | 6 +-
examples/sample-logdiff.make | 14 +++
examples/test-alg.make | 4 +-
examples/test-zstr.make | 61 +++++++-------
include/alg.hpp | 34 +++++---
include/logdiff.hpp | 187 +++++++++++++++++++++++++++++++++++++++++
include/logger.hpp | 11 +++
include/logsum.hpp | 47 +++++------
include/logsumset.hpp | 2 +-
10 files changed, 303 insertions(+), 79 deletions(-)
diff --git a/README.org b/README.org
index 4dd8d4c..0a55064 100644
--- a/README.org
+++ b/README.org
@@ -39,20 +39,26 @@ Collection of new and extended SL algorithms. Contents:
- =for_each[_it][_advance]=: Apply functor to elements in range. Options: take iterator pair or range as input; pass iterator instead of element to functor; store next iterator prior to applying functor (e.g. use: remove elements from a list in one pass.)
-- =(min|max|minmax)[_value]_of=: Find min and/or max in range. Options: take iterator pair or range as input; optionally use functor to extract key; return iterator or value.
+- =(min|max|minmax)[_value]_of=: Find min and/or max in range. Options: take iterator pair or range as input; optionally use functor to extract key; return iterator(s) or value(s).
- =mean_stdv_of=: Find mean and sample stdv of elements in range. Options: take iterator pair or range as input; optionally use functor to extract key.
-- =(equal|all|any)_of=: Check range of elements are equal, or that all are true, or that at least one is true. Options: take iterator pair or range as input; optionally use functor to extract key.
+- =(equal|all|any)_of=: Check that all elements in a range are equal, or that all are true, or that at least one is true. Options: take iterator pair or range as input; optionally use functor to extract key.
- =os_join=: Use =operator <<= overloads to print range or to convert it to string using a custom separator. Options: take iterator pair or range as input; optionally use functor to extract key.
#+BEGIN_EXAMPLE
#include "alg.hpp"
...
-std::vector< int > v{3, 6, 18};
-bool equal_mod_3 = alg::equal_of(v, [] (int i) { return i%3; }); // true
-std::cout << "v: " << alg::os_join(v, ", ") << std::endl; // prints "v: 3, 6, 18"
+std::list<int> l{3, 6, 10, 18};
+// erase elements == 0 mod 5 => l == {3, 6, 18}
+alg::for_each_it_advance(l, [&l] (std::list<int>::iterator it) { if (*it%5 == 0) l.erase(it); });
+// iterator to minimum value mod 5 => iterator to 6
+auto it_min_val_mod_5 = alg::min_of(l, [] (int i) { return i%5; });
+// all equal mod 3 => true
+bool equal_mod_3 = alg::equal_of(l, [] (int i) { return i%3; });
+// to ostream => "l: 3, 6, 18"
+std::cout << "l: " << alg::os_join(l, ", ") << std::endl;
#+END_EXAMPLE
***** logsum
diff --git a/examples/benchmark-logsum.make b/examples/benchmark-logsum.make
index 2ae5d23..c8dd484 100644
--- a/examples/benchmark-logsum.make
+++ b/examples/benchmark-logsum.make
@@ -1,13 +1,11 @@
SHELL := /bin/bash
-DIR =
-
.PHONY: all sample clean
all: benchmark-logsum
-benchmark-logsum: ../include/logsum.hpp ${DIR}/logsum.h ${DIR}/logsum.cpp
- g++ -std=c++11 -O2 -Wall -Wextra -pedantic -I ${DIR} -DBENCHMARK_p7LOGSUM -x c++ ../include/logsum.hpp ${DIR}/logsum.cpp -o $@
+benchmark-logsum: ../include/logsum.hpp
+ g++ -std=c++11 -O2 -Wall -Wextra -pedantic -DBENCHMARK_p7LOGSUM -x c++ ../include/logsum.hpp -o $@
run: benchmark-logsum
./benchmark-logsum 42 0
diff --git a/examples/sample-logdiff.make b/examples/sample-logdiff.make
new file mode 100644
index 0000000..47cb16b
--- /dev/null
+++ b/examples/sample-logdiff.make
@@ -0,0 +1,14 @@
+SHELL := /bin/bash
+
+.PHONY: all sample clean
+
+all: sample-logdiff
+
+sample-logdiff: ../include/logdiff.hpp
+ g++ -std=c++11 -O0 -g3 -ggdb -fno-eliminate-unused-debug-types -Wall -Wextra -pedantic -DSAMPLE_LOGDIFF -x c++ $< -o $@
+
+sample: sample-logdiff
+ ./sample-logdiff .5 .2
+
+clean:
+ rm -rf sample-logdiff
diff --git a/examples/test-alg.make b/examples/test-alg.make
index 3b77385..f4b514d 100644
--- a/examples/test-alg.make
+++ b/examples/test-alg.make
@@ -5,10 +5,10 @@ SHELL := /bin/bash
all: test-alg
%: %.cpp
- g++ -std=c++11 -O0 -g3 -ggdb -fno-eliminate-unused-debug-types -Wall -Wextra -pedantic -I../include -o $@ $^ -lz
+ ${DOCKER_CMD} ${CXX} -std=c++11 -O0 -g3 -ggdb -fno-eliminate-unused-debug-types -Wall -Wextra -pedantic -I../include -o $@ $^ -lz
test: test-alg
- ./test-alg
+ ${DOCKER_CMD} ./test-alg
clean:
rm -rf test-alg
diff --git a/examples/test-zstr.make b/examples/test-zstr.make
index 7cf165b..c995dc3 100644
--- a/examples/test-zstr.make
+++ b/examples/test-zstr.make
@@ -5,38 +5,39 @@ SHELL := /bin/bash
all: test-strict_fstream ztxtpipe zpipe zc
%: %.cpp
- g++ -std=c++11 -O0 -g3 -ggdb -fno-eliminate-unused-debug-types -Wall -Wextra -pedantic -I../include -o $@ $^ -lz
+ ${DOCKER_CMD} ${CXX} -std=c++11 -O0 -g3 -ggdb -fno-eliminate-unused-debug-types -Wall -Wextra -pedantic -I../include -o $@ $^ -lz
test: ztxtpipe zpipe zc
- diff -q ztxtpipe.cpp <(./ztxtpipe <ztxtpipe.cpp)
- diff -q ztxtpipe.cpp <(gzip <ztxtpipe.cpp | ./ztxtpipe)
- diff -q ztxtpipe.cpp <(gzip -9 <ztxtpipe.cpp | ./ztxtpipe)
- diff -q <(cat ztxtpipe.cpp ztxtpipe.cpp) <(cat ztxtpipe.cpp ztxtpipe.cpp | gzip | ./ztxtpipe)
- diff -q <(cat ztxtpipe.cpp ztxtpipe.cpp) <({ gzip <ztxtpipe.cpp; gzip <ztxtpipe.cpp; } | ./ztxtpipe)
- diff -q zpipe.cpp <(./zpipe <zpipe.cpp)
- diff -q zpipe.cpp <(gzip <zpipe.cpp | ./zpipe)
- diff -q zpipe.cpp <(gzip -9 <zpipe.cpp | ./zpipe)
- diff -q <(cat zpipe.cpp zpipe.cpp) <(cat zpipe.cpp zpipe.cpp | gzip | ./zpipe)
- diff -q <(cat zpipe.cpp zpipe.cpp) <({ gzip <zpipe.cpp; gzip <zpipe.cpp; } | ./zpipe)
- diff -q <(<zpipe.cpp | gzip) <(<zpipe.cpp | gzip | gzip | ./zpipe)
- diff -q zc.cpp <(./zc <zc.cpp)
- diff -q zc.cpp <(./zc - <zc.cpp)
- diff -q zc.cpp <(./zc - - <zc.cpp)
- diff -q zc.cpp <(./zc zc.cpp)
- diff -q zc.cpp <(<zc.cpp gzip | ./zc)
- diff -q zc.cpp <(<zc.cpp gzip | ./zc -)
- diff -q zc.cpp <(<zc.cpp gzip | ./zc - -)
- diff -q zc.cpp <(./zc <(<zc.cpp gzip))
- diff -q zc.cpp <(<zc.cpp ./zc -c | zcat)
- diff -q zc.cpp <(<zc.cpp ./zc -c - | zcat)
- diff -q zc.cpp <(<zc.cpp ./zc -c - - | zcat)
- diff -q zc.cpp <(./zc -c zc.cpp | zcat)
- diff -q <(cat zc.cpp zc.cpp) <(./zc -c zc.cpp - <zc.cpp | zcat)
- diff -q <(cat zc.cpp zc.cpp) <(./zc -c - zc.cpp <zc.cpp | zcat)
- diff -q <(cat zc.cpp zc.cpp) <(./zc -c zc.cpp zc.cpp | zcat)
- diff -q <(cat zc.cpp zc.cpp) <({ ./zc -c zc.cpp; ./zc -c zc.cpp; } | zcat)
- diff -q <(cat zc.cpp zc.cpp) <({ gzip <zc.cpp; ./zc -c zc.cpp; } | zcat)
- diff -q <(cat zc.cpp zc.cpp) <({ ./zc -c zc.cpp; gzip <zc.cpp; } | zcat)
+ cat ztxtpipe.cpp | ${DOCKER_CMD} ./ztxtpipe | diff -q - ztxtpipe.cpp
+ cat ztxtpipe.cpp | gzip | ${DOCKER_CMD} ./ztxtpipe | diff -q - ztxtpipe.cpp
+ cat ztxtpipe.cpp ztxtpipe.cpp | gzip | ${DOCKER_CMD} ./ztxtpipe | diff -q - <(cat ztxtpipe.cpp ztxtpipe.cpp)
+ { gzip <ztxtpipe.cpp; gzip <ztxtpipe.cpp; } | ${DOCKER_CMD} ./ztxtpipe | diff -q - <(cat ztxtpipe.cpp ztxtpipe.cpp)
+
+ cat zpipe.cpp | ${DOCKER_CMD} ./zpipe | diff -q - zpipe.cpp
+ cat zpipe.cpp | gzip | ${DOCKER_CMD} ./zpipe | diff -q - zpipe.cpp
+ cat zpipe.cpp zpipe.cpp | gzip | ${DOCKER_CMD} ./zpipe | diff -q - <(cat zpipe.cpp zpipe.cpp)
+ { gzip <zpipe.cpp; gzip <zpipe.cpp; } | ${DOCKER_CMD} ./zpipe | diff -q - <(cat zpipe.cpp zpipe.cpp)
+ cat zpipe.cpp | gzip | gzip | ${DOCKER_CMD} ./zpipe | diff -q - <(cat zpipe.cpp | gzip)
+
+ cat zc.cpp | ${DOCKER_CMD} ./zc | diff -q - zc.cpp
+ cat zc.cpp | ${DOCKER_CMD} ./zc - | diff -q - zc.cpp
+ cat zc.cpp | ${DOCKER_CMD} ./zc - - | diff -q - zc.cpp
+ ${DOCKER_CMD} ./zc zc.cpp | diff -q - zc.cpp
+ cat zc.cpp | gzip | ${DOCKER_CMD} ./zc | diff -q - zc.cpp
+ cat zc.cpp | gzip | ${DOCKER_CMD} ./zc - | diff -q - zc.cpp
+ cat zc.cpp | gzip | ${DOCKER_CMD} ./zc - - | diff -q - zc.cpp
+
+ cat zc.cpp | ${DOCKER_CMD} ./zc -c | zcat | diff -q - zc.cpp
+ cat zc.cpp | ${DOCKER_CMD} ./zc -c - | zcat | diff -q - zc.cpp
+ cat zc.cpp | ${DOCKER_CMD} ./zc -c - - | zcat | diff -q - zc.cpp
+ ${DOCKER_CMD} ./zc -c zc.cpp | zcat | diff -q - zc.cpp
+
+ cat zc.cpp | ${DOCKER_CMD} ./zc -c zc.cpp - | zcat | diff -q - <(cat zc.cpp zc.cpp)
+ cat zc.cpp | ${DOCKER_CMD} ./zc -c - zc.cpp | zcat | diff -q - <(cat zc.cpp zc.cpp)
+ ${DOCKER_CMD} ./zc -c zc.cpp zc.cpp | zcat | diff -q - <(cat zc.cpp zc.cpp)
+ { ${DOCKER_CMD} ./zc -c zc.cpp; ${DOCKER_CMD} ./zc -c zc.cpp; } | zcat | diff -q - <(cat zc.cpp zc.cpp)
+ { ${DOCKER_CMD} ./zc -c zc.cpp; gzip <zc.cpp; } | zcat | diff -q - <(cat zc.cpp zc.cpp)
+ { gzip <zc.cpp; ${DOCKER_CMD} ./zc -c zc.cpp; } | zcat | diff -q - <(cat zc.cpp zc.cpp)
@echo "all passed"
clean:
diff --git a/include/alg.hpp b/include/alg.hpp
index 313efa0..776c8f3 100644
--- a/include/alg.hpp
+++ b/include/alg.hpp
@@ -68,7 +68,10 @@
#define __ALG_HPP
#include <algorithm>
+#include <cmath>
#include <iostream>
+#include <iterator>
+#include <numeric>
#include <sstream>
#include <type_traits>
@@ -290,13 +293,13 @@ min_of(Forward_Range&& rg, Unary_Function&& fn = Unary_Function())
* @param fn Functor used to obtain key from value.
*/
template < class Forward_Iterator, class Unary_Function = detail::Identity >
-typename std::result_of< Unary_Function(typename Forward_Iterator::value_type) >::type
+typename std::result_of< Unary_Function(typename std::iterator_traits< Forward_Iterator >::value_type) >::type
min_value_of(Forward_Iterator first, Forward_Iterator last, Unary_Function&& fn = Unary_Function())
{
auto it = min_of(first, last, std::forward< Unary_Function >(fn));
return it != last
? fn(*it)
- : typename std::result_of< Unary_Function(typename Forward_Iterator::value_type) >::type();
+ : typename std::result_of< Unary_Function(typename std::iterator_traits< Forward_Iterator >::value_type) >::type();
}
/**
@@ -307,7 +310,7 @@ min_value_of(Forward_Iterator first, Forward_Iterator last, Unary_Function&& fn
template < class Forward_Range, class Unary_Function = detail::Identity >
auto
min_value_of(Forward_Range&& rg, Unary_Function&& fn = Unary_Function())
- -> typename std::result_of< Unary_Function(typename decltype(rg.begin())::value_type) >::type
+ -> typename std::result_of< Unary_Function(typename std::iterator_traits< decltype(rg.begin()) >::value_type) >::type
{
return min_value_of(rg.begin(), rg.end(), std::forward< Unary_Function >(fn));
}
@@ -357,13 +360,13 @@ max_of(Forward_Range&& rg, Unary_Function&& fn = Unary_Function())
* @param fn Functor used to obtain key from value.
*/
template < class Forward_Iterator, class Unary_Function = detail::Identity >
-typename std::result_of< Unary_Function(typename Forward_Iterator::value_type) >::type
+typename std::result_of< Unary_Function(typename std::iterator_traits< Forward_Iterator >::value_type) >::type
max_value_of(Forward_Iterator first, Forward_Iterator last, Unary_Function&& fn = Unary_Function())
{
auto it = max_of(first, last, std::forward< Unary_Function >(fn));
return it != last
? fn(*it)
- : typename std::result_of< Unary_Function(typename Forward_Iterator::value_type) >::type();
+ : typename std::result_of< Unary_Function(typename std::iterator_traits< Forward_Iterator >::value_type) >::type();
}
/**
@@ -374,7 +377,7 @@ max_value_of(Forward_Iterator first, Forward_Iterator last, Unary_Function&& fn
template < class Forward_Range, class Unary_Function = detail::Identity >
auto
max_value_of(Forward_Range&& rg, Unary_Function&& fn = Unary_Function())
- -> typename std::result_of< Unary_Function(typename decltype(rg.begin())::value_type) >::type
+ -> typename std::result_of< Unary_Function(typename std::iterator_traits< decltype(rg.begin()) >::value_type) >::type
{
return max_value_of(rg.begin(), rg.end(), std::forward< Unary_Function >(fn));
}
@@ -428,15 +431,15 @@ minmax_of(Forward_Range&& rg, Unary_Function&& fn = Unary_Function())
* @param fn Functor used to obtain key from value.
*/
template < class Forward_Iterator, class Unary_Function = detail::Identity >
-std::pair< typename std::result_of< Unary_Function(typename Forward_Iterator::value_type) >::type,
- typename std::result_of< Unary_Function(typename Forward_Iterator::value_type) >::type >
+std::pair< typename std::result_of< Unary_Function(typename std::iterator_traits< Forward_Iterator >::value_type) >::type,
+ typename std::result_of< Unary_Function(typename std::iterator_traits< Forward_Iterator >::value_type) >::type >
minmax_value_of(Forward_Iterator first, Forward_Iterator last, Unary_Function&& fn = Unary_Function())
{
auto p = minmax_of(first, last, std::forward< Unary_Function >(fn));
return p.first != last
? std::make_pair(fn(*p.first), fn(*p.second))
- : std::make_pair(typename std::result_of< Unary_Function(typename Forward_Iterator::value_type) >::type(),
- typename std::result_of< Unary_Function(typename Forward_Iterator::value_type) >::type());
+ : std::make_pair(typename std::result_of< Unary_Function(typename std::iterator_traits< Forward_Iterator >::value_type) >::type(),
+ typename std::result_of< Unary_Function(typename std::iterator_traits< Forward_Iterator >::value_type) >::type());
}
/**
@@ -447,8 +450,8 @@ minmax_value_of(Forward_Iterator first, Forward_Iterator last, Unary_Function&&
template < class Forward_Range, class Unary_Function = detail::Identity >
auto
minmax_value_of(Forward_Range&& rg, Unary_Function&& fn = Unary_Function())
- -> std::pair< typename std::result_of< Unary_Function(typename decltype(rg.begin())::value_type) >::type,
- typename std::result_of< Unary_Function(typename decltype(rg.begin())::value_type) >::type >
+ -> std::pair< typename std::result_of< Unary_Function(typename std::iterator_traits< decltype(rg.begin()) >::value_type) >::type,
+ typename std::result_of< Unary_Function(typename std::iterator_traits< decltype(rg.begin()) >::value_type) >::type >
{
return minmax_value_of(rg.begin(), rg.end(), std::forward< Unary_Function >(fn));
}
@@ -624,13 +627,15 @@ g++ -std=c++11 -D SAMPLE_ALG -x c++ alg.hpp -o sample-alg
*/
#include <vector>
+#include <array>
using namespace std;
using namespace alg;
int main()
{
- vector< unsigned > v{ 42, 1, 15 };
+ array< array< unsigned, 3 >, 3 > v2 = {{ {{ 42, 1, 15 }}, {{1, 2, 3}}, {{4, 5, 6}} }};
+ array< unsigned, 3 >& v = v2[0];
cout << "v: " << os_join(v.begin(), v.end(), ", ") << endl;
cout << "v[0]: " << os_join(v.begin(), v.begin() + 1, ", ") << endl;
cout << "v+1: " << os_join(v.begin(), v.end(), ", ", [] (unsigned i) { return i + 1; }) << endl;
@@ -639,6 +644,9 @@ int main()
cout << "u: " << os_join(vector< int >{ 5, 17, 6 }, ";") << endl;
string s = os_join(vector< char >{ 'a', 'b', 'c' }, "-");
cout << "s: " << s << endl;
+ cout << "min: " << min_value_of(v) << endl;
+ cout << "max: " << max_value_of(v) << endl;
+ cout << "minmax: " << minmax_value_of(v).first << "," << minmax_value_of(v).second << endl;
}
#endif
diff --git a/include/logdiff.hpp b/include/logdiff.hpp
new file mode 100644
index 0000000..7867ad5
--- /dev/null
+++ b/include/logdiff.hpp
@@ -0,0 +1,187 @@
+//
+// logdiff -- inspired by Sean Eddy's fast table-driven log sum
+//
+// The MIT License (MIT)
+//
+// (C) Matei David, Ontario Institute for Cancer Research, 2015
+//
+
+#ifndef __LOGDIFF_HPP
+#define __LOGDIFF_HPP
+
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+
+/*
+ * LOGDIFF_SCALE defines the precision of the calculation; the
+ * default of 1000.0 means rounding differences to the nearest 0.001
+ * nat. LOGDIFF_TBL defines the size of the lookup table; the
+ * default of 16000 means entries are calculated for differences of 0
+ * to 16.000 nats (when LOGDIFF_SCALE is 1000.0).
+ */
+#define LOGDIFF_TBL 16000
+#define LOGDIFF_SCALE 1000.0f
+
+namespace logdiff
+{
+
+namespace detail
+{
+
+struct LogDiff_Helper
+{
+ /* Function: table()
+ * Synposis: Holds the main lookup table used in logspace diff computations.
+ *
+ * Purpose: Encapsulate static array inside a function to avoid the need for a separate definition.
+ */
+ static inline float* table(void)
+ {
+ static float _table[LOGDIFF_TBL]; /* LOGDIFF_TBL=16000: (A-B) = 0..16 nats, steps of 0.001 */
+ return _table;
+ }
+
+ /* Function: init()
+ * Synopsis: Initialize the logdiff table.
+ *
+ * Purpose: Initialize the logdiff lookup table for logdiff().
+ * This function must be called once before any
+ * call to logdiff().
+ *
+ * Note: The precision of the lookup table is determined
+ * by the compile-time <LOGDIFF_TBL> constant.
+ *
+ * Returns: true on success.
+ */
+ static bool init()
+ {
+ static bool first_time = true;
+ if (not first_time)
+ {
+ return true;
+ }
+ first_time = false;
+
+ for (unsigned i = 0; i < LOGDIFF_TBL; i++)
+ {
+ auto x = -((double) i / LOGDIFF_SCALE);
+ auto e_x = std::exp(x);
+ table()[i] = std::log(1.0 - e_x);
+ }
+
+ return true;
+ }
+
+ /* Function: logdiff()
+ * Synopsis: Approximate $\log(e^a - e^b)$.
+ *
+ * Purpose: Returns a fast table-driven approximation to
+ * $\log(e^a - e^b)$.
+ *
+ * Either <a> or <b> (or both) may be $-\infty$,
+ * but neither may be $+\infty$ or <NaN>.
+ *
+ * Note: This function is a critical optimization target, because
+ * it's in the inner loop of generic Forward() algorithms.
+ */
+ static inline float logdiff(float a, float b)
+ {
+ if (a < b) std::swap(a, b);
+ return logdiff_nocomp(a, b);
+ }
+ static inline float logdiff_nocomp(float a, float b)
+ {
+ // static object whose constructor initializes the lookup table
+ static Table_Initializer _init;
+ (void)_init;
+
+ return (b == -INFINITY or a - b >= (float)(LOGDIFF_TBL - 1) / LOGDIFF_SCALE)
+ ? a
+ : a + table()[(int)((a - b) * LOGDIFF_SCALE)];
+ }
+
+ /* Function: logdiff_error()
+ * Synopsis: Compute absolute error in probability from logdiff.
+ *
+ * Purpose: Compute the absolute error in probability space
+ * resulting from logdiff()'s table lookup:
+ * approximation result - exact result.
+ */
+ static float logdiff_error(float a, float b)
+ {
+ if (a < b) std::swap(a, b);
+ return logdiff_error_nocomp(a, b);
+ }
+ static float logdiff_error_nocomp(float a, float b)
+ {
+ assert(a >= b);
+ if (a == b) return 0.0;
+ float approx = logdiff_nocomp(a, b);
+ float exact = std::log(std::exp(a) - std::exp(b));
+ return std::exp(approx) - std::exp(exact);
+ }
+
+ /* Struct: Table_Initializer
+ * Purpose: Initialize the lookup table on construction.
+ */
+ struct Table_Initializer
+ {
+ Table_Initializer()
+ {
+ init();
+ }
+ }; // struct Table_Initializer
+}; // struct LogDiff_Helper
+
+} // namespace detail
+
+/*
+ * Publish main methods outside of the helper struct.
+ */
+inline float LogDiff(float a, float b) { return detail::LogDiff_Helper::logdiff(a, b); }
+inline float LogDiffError(float a, float b) { return detail::LogDiff_Helper::logdiff_error(a, b); }
+
+} // namespace logdiff
+
+#endif
+
+#ifdef SAMPLE_LOGDIFF
+
+/*
+
+g++ -std=c++11 -DSAMPLE_LOGDIFF -x c++ logdiff.hpp -o sample-logdiff
+
+*/
+
+#include <iostream>
+#include <sstream>
+
+using namespace std;
+using namespace logdiff;
+
+int main(int argc, char* argv[])
+{
+ if (argc != 3)
+ {
+ cerr << "use: " << argv[0] << " <a> <b>" << endl;
+ return EXIT_FAILURE;
+ }
+
+ float a;
+ float b;
+ float result;
+
+ istringstream(argv[1]) >> a;
+ istringstream(argv[2]) >> b;
+
+ result = LogDiff(a, b);
+ cout << "logdiff(" << a << ", " << b << ") = " << result << endl;
+
+ result = log(exp(a) - exp(b));
+ cout << "log(exp(" << a << ") - exp(" << b << ")) = " << result << endl;
+
+ cout << "Absolute error in probability: " << LogDiffError(a, b) << endl;
+}
+
+#endif
diff --git a/include/logger.hpp b/include/logger.hpp
index 3b433af..47b54d4 100644
--- a/include/logger.hpp
+++ b/include/logger.hpp
@@ -253,11 +253,22 @@ private:
// we need 2-level indirection in order to trigger expansion after token pasting
// http://stackoverflow.com/questions/1597007/creating-c-macro-with-and-line-token-concatenation-with-positioning-macr
// http://stackoverflow.com/a/11763196/717706
+#ifdef WIN32
+#define __EXPAND(...) __VA_ARGS__
+#define __LOG_aux2(N, ...) __EXPAND(__LOG_ ## N (__VA_ARGS__))
+#else
#define __LOG_aux2(N, ...) __LOG_ ## N (__VA_ARGS__)
+#endif
+
#define __LOG_aux1(N, ...) __LOG_aux2(N, __VA_ARGS__)
#define __NARGS_AUX(_0, _1, _2, _3, _4, _5, _6, _7, _8, _9, ...) _9
+
+#ifdef WIN32
+#define __NARGS(...) __EXPAND(__NARGS_AUX(__VA_ARGS__, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0, 0))
+#else
#define __NARGS(...) __NARGS_AUX(__VA_ARGS__, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0, 0)
+#endif
#ifndef LOG_FACILITY
#define LOG_FACILITY "main"
diff --git a/include/logsum.hpp b/include/logsum.hpp
index ee22f2b..9bf565c 100644
--- a/include/logsum.hpp
+++ b/include/logsum.hpp
@@ -143,13 +143,14 @@ struct p7_FLogsum_Helper
{
// static object whose constructor initializes the lookup table
static Table_Initializer _init;
+ (void)_init;
const float max = ESL_MAX(a, b);
const float min = ESL_MIN(a, b);
- return (min == -eslINFINITY || (max-min) >= 15.7f)
+ return (min == -eslINFINITY or max - min >= (float)(p7_LOGSUM_TBL - 1) / p7_LOGSUM_SCALE)
? max
- : max + flogsum_lookup()[(int)((max-min)*p7_LOGSUM_SCALE)];
+ : max + flogsum_lookup()[(int)((max - min) * p7_LOGSUM_SCALE)];
}
/* Function: p7_FLogsumError()
@@ -259,20 +260,21 @@ Run:
#include <chrono>
#include <iostream>
#include <random>
-#include <set>
+#include <deque>
#include <sstream>
#include <vector>
-#include "logsum.h"
-
using namespace std;
+using namespace logsum;
int main(int argc, char* argv[])
{
if (argc < 3)
{
cerr << "use: " << argv[0] << " <seed> <version>" << endl
- << "where <version> is 0 for old, 1 for new" << endl;
+ << "where <version> means:" << endl
+ << " 0: use exp&log" << endl
+ << " 1: use table lookup" << endl;
return EXIT_FAILURE;
}
size_t seed = 0;
@@ -285,47 +287,44 @@ int main(int argc, char* argv[])
}
clog << "seed: " << seed << endl;
clog << "version: " << version << endl;
- const unsigned n = 1000000;
- set< float > s;
+ const unsigned n = 100000000;
+ std::deque< std::pair< float, float > > s;
mt19937 rg(seed);
uniform_real_distribution< float > unif;
for (unsigned i = 0; i < n; ++i)
{
- s.insert(unif(rg));
+ s.emplace_back(unif(rg), unif(rg));
}
- p7_FLogsumInit();
+ float res = 0.0;
// start timer
auto start_time = chrono::high_resolution_clock::now();
if (version == 0)
{
- while (s.size() > 1)
+ for (unsigned i = 0; i < s.size(); ++i)
{
- float a = *s.begin();
- s.erase(s.begin());
- float b = *s.begin();
- s.erase(s.begin());
- s.insert(::p7_FLogsum(a, b));
+ float& a = s[i].first;
+ float& b = s[i].second;
+ res = std::log(std::exp(a) + std::exp(b));
+ (void)res;
}
}
else if (version == 1)
{
- while (s.size() > 1)
+ for (unsigned i = 0; i < s.size(); ++i)
{
- float a = *s.begin();
- s.erase(s.begin());
- float b = *s.begin();
- s.erase(s.begin());
- s.insert(logsum::p7_FLogsum(a, b));
+ float& a = s[i].first;
+ float& b = s[i].second;
+ res = logsum::p7_FLogsum(a, b);
+ (void)res;
}
}
// end timer
auto end_time = chrono::high_resolution_clock::now();
- cout << "time: " << chrono::duration_cast< chrono::milliseconds >(end_time - start_time).count() << endl
- << "result: " << *s.begin() << endl;
+ cout << "time: " << chrono::duration_cast< chrono::milliseconds >(end_time - start_time).count() << endl;
}
#endif
diff --git a/include/logsumset.hpp b/include/logsumset.hpp
index 7e7aed5..23bd2ab 100644
--- a/include/logsumset.hpp
+++ b/include/logsumset.hpp
@@ -67,7 +67,7 @@ public:
assert(not std::isnan(b));
_val_set.erase(_val_set.begin());
#ifdef LOG
- if (b - a > 15.7 and b > -80)
+ if (not std::isinf(a) and b - a > 15.7 and b > -80)
{
LOG("logsumset", warning)
<< "precision loss: a=" << a << " b=" << b << std::endl;
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/libhpptools.git
More information about the debian-med-commit
mailing list