A toolkit for working with phylogenetic data.
v0.24.0
common.hpp
Go to the documentation of this file.
1 #ifndef GENESIS_UTILS_MATH_COMMON_H_
2 #define GENESIS_UTILS_MATH_COMMON_H_
3 
4 /*
5  Genesis - A toolkit for working with phylogenetic data.
6  Copyright (C) 2014-2020 Lucas Czech and HITS gGmbH
7 
8  This program is free software: you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation, either version 3 of the License, or
11  (at your option) any later version.
12 
13  This program is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with this program. If not, see <http://www.gnu.org/licenses/>.
20 
21  Contact:
22  Lucas Czech <lucas.czech@h-its.org>
23  Exelixis Lab, Heidelberg Institute for Theoretical Studies
24  Schloss-Wolfsbrunnenweg 35, D-69118 Heidelberg, Germany
25 */
26 
34 #include <algorithm>
35 #include <cassert>
36 #include <cmath>
37 #include <functional>
38 #include <limits>
39 #include <stdexcept>
40 #include <type_traits>
41 #include <utility>
42 #include <vector>
43 
44 namespace genesis {
45 namespace utils {
46 
47 // =================================================================================================
48 // Constants
49 // =================================================================================================
50 
54 constexpr double PI = 3.141592653589793238463;
55 
56 inline double circumference( double radius )
57 {
58  return 2 * PI * radius;
59 }
60 
61 // =================================================================================================
62 // Number Handling
63 // =================================================================================================
64 
71 template< typename T > inline constexpr
72 T abs_diff( T const& lhs, T const& rhs )
73 {
74  return ((lhs > rhs) ? (lhs - rhs) : (rhs - lhs));
75 }
76 
80 template <typename T> inline constexpr
81 int signum(T x, std::false_type )
82 {
83  // The tag type `std::false_type` is unnamed in order to avoid compiler warnings about an
84  // unused parameter. As this function is `constexpr`, we cannot shut this warning down by
85  // using `(void) param`, so making it unnamed is a reasonable solution in this case.
86  return T(0) < x;
87 }
88 
92 template <typename T> inline constexpr
93 int signum(T x, std::true_type )
94 {
95  // The tag type `std::false_type` is unnamed in order to avoid compiler warnings about an
96  // unused parameter. As this function is `constexpr`, we cannot shut this warning down by
97  // using `(void) param`, so making it unnamed is a reasonable solution in this case.
98  return (T(0) < x) - (x < T(0));
99 }
100 
109 template <typename T> inline constexpr
110 int signum(T x)
111 {
112  return signum( x, std::is_signed<T>() );
113 }
114 
119  double lhs,
120  double rhs,
121  double max_rel_diff = std::numeric_limits<double>::epsilon()
122 ) {
123  // Calculate the difference.
124  auto diff = std::abs( lhs - rhs );
125 
126  // Find the larger number.
127  auto largest = std::max( std::abs( lhs ), std::abs( rhs ));
128 
129  // Do the comparison.
130  return ( diff <= largest * max_rel_diff );
131 }
132 
136 inline double round_to( double x, size_t accuracy_order )
137 {
138  double factor = std::pow( 10, accuracy_order );
139  return std::round( x * factor ) / factor;
140 }
141 
152 inline size_t int_pow( size_t base, size_t exp )
153 {
154  // Using Exponentiation by squaring, see
155  // http://stackoverflow.com/a/101613/4184258
156  size_t result = 1;
157  while( exp ) {
158  if( exp & 1 ) {
159  result *= base;
160  }
161  exp >>= 1;
162  base *= base;
163  }
164  return result;
165 }
166 
172 inline bool is_valid_int_pow( size_t base, size_t exp )
173 {
174  return std::pow( base, exp ) < static_cast<double>( std::numeric_limits<size_t>::max() );
175 }
176 
177 // =================================================================================================
178 // Helper Functions
179 // =================================================================================================
180 
188 template <class ForwardIteratorA, class ForwardIteratorB>
189 std::pair<std::vector<double>, std::vector<double>> finite_pairs(
190  ForwardIteratorA first_a, ForwardIteratorA last_a,
191  ForwardIteratorB first_b, ForwardIteratorB last_b
192 ) {
193  // Prepare result.
194  std::vector<double> vec_a;
195  std::vector<double> vec_b;
196 
197  // Iterate in parallel.
198  while( first_a != last_a && first_b != last_b ) {
199  if( std::isfinite( *first_a ) && std::isfinite( *first_b ) ) {
200  vec_a.push_back( *first_a );
201  vec_b.push_back( *first_b );
202  }
203  ++first_a;
204  ++first_b;
205  }
206  if( first_a != last_a || first_b != last_b ) {
207  throw std::runtime_error(
208  "Ranges need to have same length."
209  );
210  }
211 
212  assert( vec_a.size() == vec_b.size() );
213  return { vec_a, vec_b };
214 }
215 
222 template <class ForwardIteratorA, class ForwardIteratorB>
224  ForwardIteratorA first_a, ForwardIteratorA last_a,
225  ForwardIteratorB first_b, ForwardIteratorB last_b,
226  std::function<void( double value_a, double value_b )> execute
227 ) {
228  // Iterate in parallel.
229  while( first_a != last_a && first_b != last_b ) {
230  if( std::isfinite( *first_a ) && std::isfinite( *first_b ) ) {
231  execute( *first_a, *first_b );
232  }
233  ++first_a;
234  ++first_b;
235  }
236 
237  // Check if the ranges have the same length.
238  if( first_a != last_a || first_b != last_b ) {
239  throw std::runtime_error( "Ranges need to have same length." );
240  }
241 }
242 
243 } // namespace utils
244 } // namespace genesis
245 
246 #endif // include guard
bool is_valid_int_pow(size_t base, size_t exp)
Return whether the given power can be stored within a size_t.
Definition: common.hpp:172
size_t int_pow(size_t base, size_t exp)
Calculate the power base^exp for positive integer values.
Definition: common.hpp:152
std::pair< std::vector< double >, std::vector< double > > finite_pairs(ForwardIteratorA first_a, ForwardIteratorA last_a, ForwardIteratorB first_b, ForwardIteratorB last_b)
Helper function that cleans two ranges of double of the same length from non-finite values...
Definition: common.hpp:189
constexpr int signum(T x, std::false_type)
Implementation of signum(T x) for unsigned types. See there for details.
Definition: common.hpp:81
Container namespace for all symbols of genesis in order to keep them separate when used as a library...
double round_to(double x, size_t accuracy_order)
Retun the value of x, rounded to the decimal digit given by accuracy_order.
Definition: common.hpp:136
constexpr double PI
Make the world go round.
Definition: common.hpp:54
void for_each_finite_pair(ForwardIteratorA first_a, ForwardIteratorA last_a, ForwardIteratorB first_b, ForwardIteratorB last_b, std::function< void(double value_a, double value_b)> execute)
Iterate two ranges of double values in parallel, and execute a function for each pair of values from ...
Definition: common.hpp:223
double circumference(double radius)
Definition: common.hpp:56
constexpr T abs_diff(T const &lhs, T const &rhs)
Calculate the absolute differenence between two values.
Definition: common.hpp:72
bool almost_equal_relative(double lhs, double rhs, double max_rel_diff=std::numeric_limits< double >::epsilon())
Check whether two doubles are almost equal, using a relative epsilon to compare them.
Definition: common.hpp:118