Commit e3d2bb38d8151d216c632cd7eaacf6b21c917e4e - libhpptools

Imported Upstream version 1.1.1 Andreas Tille 7 years ago

10 changed file(s) with 304 addition(s) and 80 deletion(s). Raw diff Collapse all Expand all

+11

-5

README.org less more

38	38
39	39	- =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.)
40	40
41		- =(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.
	41	- =(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).
42	42
43	43	- =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.
44	44
45		- =(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.
	45	- =(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.
46	46
47	47	- =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.
48	48
49	49	#+BEGIN_EXAMPLE
50	50	#include "alg.hpp"
51	51	...
52		std::vector< int > v{3, 6, 18};
53		bool equal_mod_3 = alg::equal_of(v, [] (int i) { return i%3; }); // true
54		std::cout << "v: " << alg::os_join(v, ", ") << std::endl; // prints "v: 3, 6, 18"
	52	std::list<int> l{3, 6, 10, 18};
	53	// erase elements == 0 mod 5 => l == {3, 6, 18}
	54	alg::for_each_it_advance(l, [&l] (std::list<int>::iterator it) { if (*it%5 == 0) l.erase(it); });
	55	// iterator to minimum value mod 5 => iterator to 6
	56	auto it_min_val_mod_5 = alg::min_of(l, [] (int i) { return i%5; });
	57	// all equal mod 3 => true
	58	bool equal_mod_3 = alg::equal_of(l, [] (int i) { return i%3; });
	59	// to ostream => "l: 3, 6, 18"
	60	std::cout << "l: " << alg::os_join(l, ", ") << std::endl;
55	61	#+END_EXAMPLE
56	62
57	63	***** logsum

-4

examples/benchmark-logsum.make less more

0	0	SHELL := /bin/bash
1
2		DIR =
3	1
4	2	.PHONY: all sample clean
5	3
6	4	all: benchmark-logsum
7	5
8		benchmark-logsum: ../include/logsum.hpp ${DIR}/logsum.h ${DIR}/logsum.cpp
9		g++ -std=c++11 -O2 -Wall -Wextra -pedantic -I ${DIR} -DBENCHMARK_p7LOGSUM -x c++ ../include/logsum.hpp ${DIR}/logsum.cpp -o $@
	6	benchmark-logsum: ../include/logsum.hpp
	7	g++ -std=c++11 -O2 -Wall -Wextra -pedantic -DBENCHMARK_p7LOGSUM -x c++ ../include/logsum.hpp -o $@
10	8
11	9	run: benchmark-logsum
12	10	./benchmark-logsum 42 0

+14

-0

examples/sample-logdiff.make less more

	0	SHELL := /bin/bash
	1
	2	.PHONY: all sample clean
	3
	4	all: sample-logdiff
	5
	6	sample-logdiff: ../include/logdiff.hpp
	7	g++ -std=c++11 -O0 -g3 -ggdb -fno-eliminate-unused-debug-types -Wall -Wextra -pedantic -DSAMPLE_LOGDIFF -x c++ $< -o $@
	8
	9	sample: sample-logdiff
	10	./sample-logdiff .5 .2
	11
	12	clean:
	13	rm -rf sample-logdiff

-2

examples/test-alg.make less more

4	4	all: test-alg
5	5
6	6	%: %.cpp
7		g++ -std=c++11 -O0 -g3 -ggdb -fno-eliminate-unused-debug-types -Wall -Wextra -pedantic -I../include -o $@ $^ -lz
	7	${DOCKER_CMD} ${CXX} -std=c++11 -O0 -g3 -ggdb -fno-eliminate-unused-debug-types -Wall -Wextra -pedantic -I../include -o $@ $^ -lz
8	8
9	9	test: test-alg
10		./test-alg
	10	${DOCKER_CMD} ./test-alg
11	11
12	12	clean:
13	13	rm -rf test-alg

+31

-30

examples/test-zstr.make less more

4	4	all: test-strict_fstream ztxtpipe zpipe zc
5	5
6	6	%: %.cpp
7		g++ -std=c++11 -O0 -g3 -ggdb -fno-eliminate-unused-debug-types -Wall -Wextra -pedantic -I../include -o $@ $^ -lz
	7	${DOCKER_CMD} ${CXX} -std=c++11 -O0 -g3 -ggdb -fno-eliminate-unused-debug-types -Wall -Wextra -pedantic -I../include -o $@ $^ -lz
8	8
9	9	test: ztxtpipe zpipe zc
10		diff -q ztxtpipe.cpp <(./ztxtpipe <ztxtpipe.cpp)
11		diff -q ztxtpipe.cpp <(gzip <ztxtpipe.cpp \| ./ztxtpipe)
12		diff -q ztxtpipe.cpp <(gzip -9 <ztxtpipe.cpp \| ./ztxtpipe)
13		diff -q <(cat ztxtpipe.cpp ztxtpipe.cpp) <(cat ztxtpipe.cpp ztxtpipe.cpp \| gzip \| ./ztxtpipe)
14		diff -q <(cat ztxtpipe.cpp ztxtpipe.cpp) <({ gzip <ztxtpipe.cpp; gzip <ztxtpipe.cpp; } \| ./ztxtpipe)
15		diff -q zpipe.cpp <(./zpipe <zpipe.cpp)
16		diff -q zpipe.cpp <(gzip <zpipe.cpp \| ./zpipe)
17		diff -q zpipe.cpp <(gzip -9 <zpipe.cpp \| ./zpipe)
18		diff -q <(cat zpipe.cpp zpipe.cpp) <(cat zpipe.cpp zpipe.cpp \| gzip \| ./zpipe)
19		diff -q <(cat zpipe.cpp zpipe.cpp) <({ gzip <zpipe.cpp; gzip <zpipe.cpp; } \| ./zpipe)
20		diff -q <(<zpipe.cpp \| gzip) <(<zpipe.cpp \| gzip \| gzip \| ./zpipe)
21		diff -q zc.cpp <(./zc <zc.cpp)
22		diff -q zc.cpp <(./zc - <zc.cpp)
23		diff -q zc.cpp <(./zc - - <zc.cpp)
24		diff -q zc.cpp <(./zc zc.cpp)
25		diff -q zc.cpp <(<zc.cpp gzip \| ./zc)
26		diff -q zc.cpp <(<zc.cpp gzip \| ./zc -)
27		diff -q zc.cpp <(<zc.cpp gzip \| ./zc - -)
28		diff -q zc.cpp <(./zc <(<zc.cpp gzip))
29		diff -q zc.cpp <(<zc.cpp ./zc -c \| zcat)
30		diff -q zc.cpp <(<zc.cpp ./zc -c - \| zcat)
31		diff -q zc.cpp <(<zc.cpp ./zc -c - - \| zcat)
32		diff -q zc.cpp <(./zc -c zc.cpp \| zcat)
33		diff -q <(cat zc.cpp zc.cpp) <(./zc -c zc.cpp - <zc.cpp \| zcat)
34		diff -q <(cat zc.cpp zc.cpp) <(./zc -c - zc.cpp <zc.cpp \| zcat)
35		diff -q <(cat zc.cpp zc.cpp) <(./zc -c zc.cpp zc.cpp \| zcat)
36		diff -q <(cat zc.cpp zc.cpp) <({ ./zc -c zc.cpp; ./zc -c zc.cpp; } \| zcat)
37		diff -q <(cat zc.cpp zc.cpp) <({ gzip <zc.cpp; ./zc -c zc.cpp; } \| zcat)
38		diff -q <(cat zc.cpp zc.cpp) <({ ./zc -c zc.cpp; gzip <zc.cpp; } \| zcat)
	10	cat ztxtpipe.cpp \| ${DOCKER_CMD} ./ztxtpipe \| diff -q - ztxtpipe.cpp
	11	cat ztxtpipe.cpp \| gzip \| ${DOCKER_CMD} ./ztxtpipe \| diff -q - ztxtpipe.cpp
	12	cat ztxtpipe.cpp ztxtpipe.cpp \| gzip \| ${DOCKER_CMD} ./ztxtpipe \| diff -q - <(cat ztxtpipe.cpp ztxtpipe.cpp)
	13	{ gzip <ztxtpipe.cpp; gzip <ztxtpipe.cpp; } \| ${DOCKER_CMD} ./ztxtpipe \| diff -q - <(cat ztxtpipe.cpp ztxtpipe.cpp)
	14
	15	cat zpipe.cpp \| ${DOCKER_CMD} ./zpipe \| diff -q - zpipe.cpp
	16	cat zpipe.cpp \| gzip \| ${DOCKER_CMD} ./zpipe \| diff -q - zpipe.cpp
	17	cat zpipe.cpp zpipe.cpp \| gzip \| ${DOCKER_CMD} ./zpipe \| diff -q - <(cat zpipe.cpp zpipe.cpp)
	18	{ gzip <zpipe.cpp; gzip <zpipe.cpp; } \| ${DOCKER_CMD} ./zpipe \| diff -q - <(cat zpipe.cpp zpipe.cpp)
	19	cat zpipe.cpp \| gzip \| gzip \| ${DOCKER_CMD} ./zpipe \| diff -q - <(cat zpipe.cpp \| gzip)
	20
	21	cat zc.cpp \| ${DOCKER_CMD} ./zc \| diff -q - zc.cpp
	22	cat zc.cpp \| ${DOCKER_CMD} ./zc - \| diff -q - zc.cpp
	23	cat zc.cpp \| ${DOCKER_CMD} ./zc - - \| diff -q - zc.cpp
	24	${DOCKER_CMD} ./zc zc.cpp \| diff -q - zc.cpp
	25	cat zc.cpp \| gzip \| ${DOCKER_CMD} ./zc \| diff -q - zc.cpp
	26	cat zc.cpp \| gzip \| ${DOCKER_CMD} ./zc - \| diff -q - zc.cpp
	27	cat zc.cpp \| gzip \| ${DOCKER_CMD} ./zc - - \| diff -q - zc.cpp
	28
	29	cat zc.cpp \| ${DOCKER_CMD} ./zc -c \| zcat \| diff -q - zc.cpp
	30	cat zc.cpp \| ${DOCKER_CMD} ./zc -c - \| zcat \| diff -q - zc.cpp
	31	cat zc.cpp \| ${DOCKER_CMD} ./zc -c - - \| zcat \| diff -q - zc.cpp
	32	${DOCKER_CMD} ./zc -c zc.cpp \| zcat \| diff -q - zc.cpp
	33
	34	cat zc.cpp \| ${DOCKER_CMD} ./zc -c zc.cpp - \| zcat \| diff -q - <(cat zc.cpp zc.cpp)
	35	cat zc.cpp \| ${DOCKER_CMD} ./zc -c - zc.cpp \| zcat \| diff -q - <(cat zc.cpp zc.cpp)
	36	${DOCKER_CMD} ./zc -c zc.cpp zc.cpp \| zcat \| diff -q - <(cat zc.cpp zc.cpp)
	37	{ ${DOCKER_CMD} ./zc -c zc.cpp; ${DOCKER_CMD} ./zc -c zc.cpp; } \| zcat \| diff -q - <(cat zc.cpp zc.cpp)
	38	{ ${DOCKER_CMD} ./zc -c zc.cpp; gzip <zc.cpp; } \| zcat \| diff -q - <(cat zc.cpp zc.cpp)
	39	{ gzip <zc.cpp; ${DOCKER_CMD} ./zc -c zc.cpp; } \| zcat \| diff -q - <(cat zc.cpp zc.cpp)
39	40	@echo "all passed"
40	41
41	42	clean:

+21

-13

include/alg.hpp less more

67	67	#define __ALG_HPP
68	68
69	69	#include <algorithm>
	70	#include <cmath>
70	71	#include <iostream>
	72	#include <iterator>
	73	#include <numeric>
71	74	#include <sstream>
72	75	#include <type_traits>
73	76

289	292	* @param fn Functor used to obtain key from value.
290	293	*/
291	294	template < class Forward_Iterator, class Unary_Function = detail::Identity >
292		typename std::result_of< Unary_Function(typename Forward_Iterator::value_type) >::type
	295	typename std::result_of< Unary_Function(typename std::iterator_traits< Forward_Iterator >::value_type) >::type
293	296	min_value_of(Forward_Iterator first, Forward_Iterator last, Unary_Function&& fn = Unary_Function())
294	297	{
295	298	auto it = min_of(first, last, std::forward< Unary_Function >(fn));
296	299	return it != last
297	300	? fn(*it)
298		: typename std::result_of< Unary_Function(typename Forward_Iterator::value_type) >::type();
	301	: typename std::result_of< Unary_Function(typename std::iterator_traits< Forward_Iterator >::value_type) >::type();
299	302	}
300	303
301	304	/**

306	309	template < class Forward_Range, class Unary_Function = detail::Identity >
307	310	auto
308	311	min_value_of(Forward_Range&& rg, Unary_Function&& fn = Unary_Function())
309		-> typename std::result_of< Unary_Function(typename decltype(rg.begin())::value_type) >::type
	312	-> typename std::result_of< Unary_Function(typename std::iterator_traits< decltype(rg.begin()) >::value_type) >::type
310	313	{
311	314	return min_value_of(rg.begin(), rg.end(), std::forward< Unary_Function >(fn));
312	315	}

356	359	* @param fn Functor used to obtain key from value.
357	360	*/
358	361	template < class Forward_Iterator, class Unary_Function = detail::Identity >
359		typename std::result_of< Unary_Function(typename Forward_Iterator::value_type) >::type
	362	typename std::result_of< Unary_Function(typename std::iterator_traits< Forward_Iterator >::value_type) >::type
360	363	max_value_of(Forward_Iterator first, Forward_Iterator last, Unary_Function&& fn = Unary_Function())
361	364	{
362	365	auto it = max_of(first, last, std::forward< Unary_Function >(fn));
363	366	return it != last
364	367	? fn(*it)
365		: typename std::result_of< Unary_Function(typename Forward_Iterator::value_type) >::type();
	368	: typename std::result_of< Unary_Function(typename std::iterator_traits< Forward_Iterator >::value_type) >::type();
366	369	}
367	370
368	371	/**

373	376	template < class Forward_Range, class Unary_Function = detail::Identity >
374	377	auto
375	378	max_value_of(Forward_Range&& rg, Unary_Function&& fn = Unary_Function())
376		-> typename std::result_of< Unary_Function(typename decltype(rg.begin())::value_type) >::type
	379	-> typename std::result_of< Unary_Function(typename std::iterator_traits< decltype(rg.begin()) >::value_type) >::type
377	380	{
378	381	return max_value_of(rg.begin(), rg.end(), std::forward< Unary_Function >(fn));
379	382	}

427	430	* @param fn Functor used to obtain key from value.
428	431	*/
429	432	template < class Forward_Iterator, class Unary_Function = detail::Identity >
430		std::pair< typename std::result_of< Unary_Function(typename Forward_Iterator::value_type) >::type,
431		typename std::result_of< Unary_Function(typename Forward_Iterator::value_type) >::type >
	433	std::pair< typename std::result_of< Unary_Function(typename std::iterator_traits< Forward_Iterator >::value_type) >::type,
	434	typename std::result_of< Unary_Function(typename std::iterator_traits< Forward_Iterator >::value_type) >::type >
432	435	minmax_value_of(Forward_Iterator first, Forward_Iterator last, Unary_Function&& fn = Unary_Function())
433	436	{
434	437	auto p = minmax_of(first, last, std::forward< Unary_Function >(fn));
435	438	return p.first != last
436	439	? std::make_pair(fn(p.first), fn(p.second))
437		: std::make_pair(typename std::result_of< Unary_Function(typename Forward_Iterator::value_type) >::type(),
438		typename std::result_of< Unary_Function(typename Forward_Iterator::value_type) >::type());
	440	: std::make_pair(typename std::result_of< Unary_Function(typename std::iterator_traits< Forward_Iterator >::value_type) >::type(),
	441	typename std::result_of< Unary_Function(typename std::iterator_traits< Forward_Iterator >::value_type) >::type());
439	442	}
440	443
441	444	/**

446	449	template < class Forward_Range, class Unary_Function = detail::Identity >
447	450	auto
448	451	minmax_value_of(Forward_Range&& rg, Unary_Function&& fn = Unary_Function())
449		-> std::pair< typename std::result_of< Unary_Function(typename decltype(rg.begin())::value_type) >::type,
450		typename std::result_of< Unary_Function(typename decltype(rg.begin())::value_type) >::type >
	452	-> std::pair< typename std::result_of< Unary_Function(typename std::iterator_traits< decltype(rg.begin()) >::value_type) >::type,
	453	typename std::result_of< Unary_Function(typename std::iterator_traits< decltype(rg.begin()) >::value_type) >::type >
451	454	{
452	455	return minmax_value_of(rg.begin(), rg.end(), std::forward< Unary_Function >(fn));
453	456	}

623	626	*/
624	627
625	628	#include <vector>
	629	#include <array>
626	630
627	631	using namespace std;
628	632	using namespace alg;
629	633
630	634	int main()
631	635	{
632		vector< unsigned > v{ 42, 1, 15 };
	636	array< array< unsigned, 3 >, 3 > v2 = {{ {{ 42, 1, 15 }}, {{1, 2, 3}}, {{4, 5, 6}} }};
	637	array< unsigned, 3 >& v = v2[0];
633	638	cout << "v: " << os_join(v.begin(), v.end(), ", ") << endl;
634	639	cout << "v[0]: " << os_join(v.begin(), v.begin() + 1, ", ") << endl;
635	640	cout << "v+1: " << os_join(v.begin(), v.end(), ", ", [] (unsigned i) { return i + 1; }) << endl;

638	643	cout << "u: " << os_join(vector< int >{ 5, 17, 6 }, ";") << endl;
639	644	string s = os_join(vector< char >{ 'a', 'b', 'c' }, "-");
640	645	cout << "s: " << s << endl;
	646	cout << "min: " << min_value_of(v) << endl;
	647	cout << "max: " << max_value_of(v) << endl;
	648	cout << "minmax: " << minmax_value_of(v).first << "," << minmax_value_of(v).second << endl;
641	649	}
642	650
643	651	#endif

+187

-0

include/logdiff.hpp less more

	0	//
	1	// logdiff -- inspired by Sean Eddy's fast table-driven log sum
	2	//
	3	// The MIT License (MIT)
	4	//
	5	// (C) Matei David, Ontario Institute for Cancer Research, 2015
	6	//
	7
	8	#ifndef __LOGDIFF_HPP
	9	#define __LOGDIFF_HPP
	10
	11	#include <algorithm>
	12	#include <cassert>
	13	#include <cmath>
	14
	15	/*
	16	* LOGDIFF_SCALE defines the precision of the calculation; the
	17	* default of 1000.0 means rounding differences to the nearest 0.001
	18	* nat. LOGDIFF_TBL defines the size of the lookup table; the
	19	* default of 16000 means entries are calculated for differences of 0
	20	* to 16.000 nats (when LOGDIFF_SCALE is 1000.0).
	21	*/
	22	#define LOGDIFF_TBL 16000
	23	#define LOGDIFF_SCALE 1000.0f
	24
	25	namespace logdiff
	26	{
	27
	28	namespace detail
	29	{
	30
	31	struct LogDiff_Helper
	32	{
	33	/* Function: table()
	34	* Synposis: Holds the main lookup table used in logspace diff computations.
	35	*
	36	* Purpose: Encapsulate static array inside a function to avoid the need for a separate definition.
	37	*/
	38	static inline float* table(void)
	39	{
	40	static float _table[LOGDIFF_TBL]; /* LOGDIFF_TBL=16000: (A-B) = 0..16 nats, steps of 0.001 */
	41	return _table;
	42	}
	43
	44	/* Function: init()
	45	* Synopsis: Initialize the logdiff table.
	46	*
	47	* Purpose: Initialize the logdiff lookup table for logdiff().
	48	* This function must be called once before any
	49	* call to logdiff().
	50	*
	51	* Note: The precision of the lookup table is determined
	52	* by the compile-time <LOGDIFF_TBL> constant.
	53	*
	54	* Returns: true on success.
	55	*/
	56	static bool init()
	57	{
	58	static bool first_time = true;
	59	if (not first_time)
	60	{
	61	return true;
	62	}
	63	first_time = false;
	64
	65	for (unsigned i = 0; i < LOGDIFF_TBL; i++)
	66	{
	67	auto x = -((double) i / LOGDIFF_SCALE);
	68	auto e_x = std::exp(x);
	69	table()[i] = std::log(1.0 - e_x);
	70	}
	71
	72	return true;
	73	}
	74
	75	/* Function: logdiff()
	76	* Synopsis: Approximate $\log(e^a - e^b)$.
	77	*
	78	* Purpose: Returns a fast table-driven approximation to
	79	* $\log(e^a - e^b)$.
	80	*
	81	* Either <a> or <b> (or both) may be $-\infty$,
	82	* but neither may be $+\infty$ or <NaN>.
	83	*
	84	* Note: This function is a critical optimization target, because
	85	* it's in the inner loop of generic Forward() algorithms.
	86	*/
	87	static inline float logdiff(float a, float b)
	88	{
	89	if (a < b) std::swap(a, b);
	90	return logdiff_nocomp(a, b);
	91	}
	92	static inline float logdiff_nocomp(float a, float b)
	93	{
	94	// static object whose constructor initializes the lookup table
	95	static Table_Initializer _init;
	96	(void)_init;
	97
	98	return (b == -INFINITY or a - b >= (float)(LOGDIFF_TBL - 1) / LOGDIFF_SCALE)
	99	? a
	100	: a + table()[(int)((a - b) * LOGDIFF_SCALE)];
	101	}
	102
	103	/* Function: logdiff_error()
	104	* Synopsis: Compute absolute error in probability from logdiff.
	105	*
	106	* Purpose: Compute the absolute error in probability space
	107	* resulting from logdiff()'s table lookup:
	108	* approximation result - exact result.
	109	*/
	110	static float logdiff_error(float a, float b)
	111	{
	112	if (a < b) std::swap(a, b);
	113	return logdiff_error_nocomp(a, b);
	114	}
	115	static float logdiff_error_nocomp(float a, float b)
	116	{
	117	assert(a >= b);
	118	if (a == b) return 0.0;
	119	float approx = logdiff_nocomp(a, b);
	120	float exact = std::log(std::exp(a) - std::exp(b));
	121	return std::exp(approx) - std::exp(exact);
	122	}
	123
	124	/* Struct: Table_Initializer
	125	* Purpose: Initialize the lookup table on construction.
	126	*/
	127	struct Table_Initializer
	128	{
	129	Table_Initializer()
	130	{
	131	init();
	132	}
	133	}; // struct Table_Initializer
	134	}; // struct LogDiff_Helper
	135
	136	} // namespace detail
	137
	138	/*
	139	* Publish main methods outside of the helper struct.
	140	*/
	141	inline float LogDiff(float a, float b) { return detail::LogDiff_Helper::logdiff(a, b); }
	142	inline float LogDiffError(float a, float b) { return detail::LogDiff_Helper::logdiff_error(a, b); }
	143
	144	} // namespace logdiff
	145
	146	#endif
	147
	148	#ifdef SAMPLE_LOGDIFF
	149
	150	/*
	151
	152	g++ -std=c++11 -DSAMPLE_LOGDIFF -x c++ logdiff.hpp -o sample-logdiff
	153
	154	*/
	155
	156	#include <iostream>
	157	#include <sstream>
	158
	159	using namespace std;
	160	using namespace logdiff;
	161
	162	int main(int argc, char* argv[])
	163	{
	164	if (argc != 3)
	165	{
	166	cerr << "use: " << argv[0] << " <a> <b>" << endl;
	167	return EXIT_FAILURE;
	168	}
	169
	170	float a;
	171	float b;
	172	float result;
	173
	174	istringstream(argv[1]) >> a;
	175	istringstream(argv[2]) >> b;
	176
	177	result = LogDiff(a, b);
	178	cout << "logdiff(" << a << ", " << b << ") = " << result << endl;
	179
	180	result = log(exp(a) - exp(b));
	181	cout << "log(exp(" << a << ") - exp(" << b << ")) = " << result << endl;
	182
	183	cout << "Absolute error in probability: " << LogDiffError(a, b) << endl;
	184	}
	185
	186	#endif

+11

-0

include/logger.hpp less more

252	252	// we need 2-level indirection in order to trigger expansion after token pasting
253	253	// http://stackoverflow.com/questions/1597007/creating-c-macro-with-and-line-token-concatenation-with-positioning-macr
254	254	// http://stackoverflow.com/a/11763196/717706
	255	#ifdef WIN32
	256	#define __EXPAND(...) __VA_ARGS__
	257	#define __LOG_aux2(N, ...) __EXPAND(__LOG_ ## N (__VA_ARGS__))
	258	#else
255	259	#define __LOG_aux2(N, ...) __LOG_ ## N (__VA_ARGS__)
	260	#endif
	261
256	262	#define __LOG_aux1(N, ...) __LOG_aux2(N, __VA_ARGS__)
257	263
258	264	#define __NARGS_AUX(_0, _1, _2, _3, _4, _5, _6, _7, _8, _9, ...) _9
	265
	266	#ifdef WIN32
	267	#define __NARGS(...) __EXPAND(__NARGS_AUX(__VA_ARGS__, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0, 0))
	268	#else
259	269	#define __NARGS(...) __NARGS_AUX(__VA_ARGS__, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0, 0)
	270	#endif
260	271
261	272	#ifndef LOG_FACILITY
262	273	#define LOG_FACILITY "main"

+24

-25

include/logsum.hpp less more

142	142	{
143	143	// static object whose constructor initializes the lookup table
144	144	static Table_Initializer _init;
	145	(void)_init;
145	146
146	147	const float max = ESL_MAX(a, b);
147	148	const float min = ESL_MIN(a, b);
148	149
149		return (min == -eslINFINITY \|\| (max-min) >= 15.7f)
	150	return (min == -eslINFINITY or max - min >= (float)(p7_LOGSUM_TBL - 1) / p7_LOGSUM_SCALE)
150	151	? max
151		: max + flogsum_lookup()[(int)((max-min)*p7_LOGSUM_SCALE)];
	152	: max + flogsum_lookup()[(int)((max - min) * p7_LOGSUM_SCALE)];
152	153	}
153	154
154	155	/* Function: p7_FLogsumError()

258	259	#include <chrono>
259	260	#include <iostream>
260	261	#include <random>
261		#include <set>
	262	#include <deque>
262	263	#include <sstream>
263	264	#include <vector>
264	265
265		#include "logsum.h"
266
267	266	using namespace std;
	267	using namespace logsum;
268	268
269	269	int main(int argc, char* argv[])
270	270	{
271	271	if (argc < 3)
272	272	{
273	273	cerr << "use: " << argv[0] << " <seed> <version>" << endl
274		<< "where <version> is 0 for old, 1 for new" << endl;
	274	<< "where <version> means:" << endl
	275	<< " 0: use exp&log" << endl
	276	<< " 1: use table lookup" << endl;
275	277	return EXIT_FAILURE;
276	278	}
277	279	size_t seed = 0;

284	286	}
285	287	clog << "seed: " << seed << endl;
286	288	clog << "version: " << version << endl;
287		const unsigned n = 1000000;
288		set< float > s;
	289	const unsigned n = 100000000;
	290	std::deque< std::pair< float, float > > s;
289	291	mt19937 rg(seed);
290	292	uniform_real_distribution< float > unif;
291	293	for (unsigned i = 0; i < n; ++i)
292	294	{
293		s.insert(unif(rg));
294		}
295		p7_FLogsumInit();
	295	s.emplace_back(unif(rg), unif(rg));
	296	}
	297	float res = 0.0;
296	298
297	299	// start timer
298	300	auto start_time = chrono::high_resolution_clock::now();
299	301
300	302	if (version == 0)
301	303	{
302		while (s.size() > 1)
	304	for (unsigned i = 0; i < s.size(); ++i)
303	305	{
304		float a = *s.begin();
305		s.erase(s.begin());
306		float b = *s.begin();
307		s.erase(s.begin());
308		s.insert(::p7_FLogsum(a, b));
	306	float& a = s[i].first;
	307	float& b = s[i].second;
	308	res = std::log(std::exp(a) + std::exp(b));
	309	(void)res;
309	310	}
310	311	}
311	312	else if (version == 1)
312	313	{
313		while (s.size() > 1)
	314	for (unsigned i = 0; i < s.size(); ++i)
314	315	{
315		float a = *s.begin();
316		s.erase(s.begin());
317		float b = *s.begin();
318		s.erase(s.begin());
319		s.insert(logsum::p7_FLogsum(a, b));
	316	float& a = s[i].first;
	317	float& b = s[i].second;
	318	res = logsum::p7_FLogsum(a, b);
	319	(void)res;
320	320	}
321	321	}
322	322
323	323	// end timer
324	324	auto end_time = chrono::high_resolution_clock::now();
325	325
326		cout << "time: " << chrono::duration_cast< chrono::milliseconds >(end_time - start_time).count() << endl
327		<< "result: " << *s.begin() << endl;
	326	cout << "time: " << chrono::duration_cast< chrono::milliseconds >(end_time - start_time).count() << endl;
328	327	}
329	328
330	329	#endif

-1

include/logsumset.hpp less more

66	66	assert(not std::isnan(b));
67	67	_val_set.erase(_val_set.begin());
68	68	#ifdef LOG
69		if (b - a > 15.7 and b > -80)
	69	if (not std::isinf(a) and b - a > 15.7 and b > -80)
70	70	{
71	71	LOG("logsumset", warning)
72	72	<< "precision loss: a=" << a << " b=" << b << std::endl;