[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