Math.h
Go to the documentation of this file.
1 /*!
2  *
3  *
4  * \brief Very basic math abstraction layer.
5  *
6  * \par
7  * This file serves as a minimal abstraction layer.
8  * Inclusion of this file makes some frequently used
9  * functions, constants, and header file inclusions
10  * OS-, compiler-, and version-independent.
11  *
12  *
13  *
14  *
15  * \author -
16  * \date -
17  *
18  *
19  * \par Copyright 1995-2015 Shark Development Team
20  *
21  * <BR><HR>
22  * This file is part of Shark.
23  * <http://image.diku.dk/shark/>
24  *
25  * Shark is free software: you can redistribute it and/or modify
26  * it under the terms of the GNU Lesser General Public License as published
27  * by the Free Software Foundation, either version 3 of the License, or
28  * (at your option) any later version.
29  *
30  * Shark is distributed in the hope that it will be useful,
31  * but WITHOUT ANY WARRANTY; without even the implied warranty of
32  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
33  * GNU Lesser General Public License for more details.
34  *
35  * You should have received a copy of the GNU Lesser General Public License
36  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
37  *
38  */
39 #ifndef SHARK_CORE_MATH_H
40 #define SHARK_CORE_MATH_H
41 
42 #include <boost/math/constants/constants.hpp>
43 #include <boost/math/special_functions/sign.hpp>
44 #include <boost/utility/enable_if.hpp>
45 #include <boost/type_traits/is_arithmetic.hpp>
46 #include <limits>
47 #include <cmath>
48 #include <shark/Core/Exception.h>
49 
50 namespace shark {
51 
52  // define doxygen group for shark global functions, var, etc. (should only appear once in all shark files)
53  /**
54  * \defgroup shark_globals shark_globals
55  *
56  * \brief Several mathematical, linear-algebra, or other functions within Shark are not part of any particular class. They are collected here in the doxygen group "shark_globals"
57  *
58  * @{
59  *
60  */
61 
62  /**
63  * \brief Constant for sqrt( 2 * pi ).
64  */
65  static const double SQRT_2_PI = boost::math::constants::root_two_pi<double>();
66 // static const double SQRT_PI = boost::math::constants::root_pi<double>();
67 
68  /// Maximum allowed input value for exp.
69  template<class T>
71  return boost::math::constants::ln_two<T>()*std::numeric_limits<T>::max_exponent;
72  }
73  /// Minimum value for exp(x) allowed so that it is not 0.
74  template<class T>
76  return boost::math::constants::ln_two<T>()*std::numeric_limits<T>::min_exponent;
77  }
78 
79  /// Calculates x^2.
80  template <class T>
81  inline typename boost::enable_if<boost::is_arithmetic<T>, T>::type sqr( const T & x) {
82  return x * x;
83  }
84 
85  /// Calculates x^3.
86  template <class T> inline T cube( const T & x) {
87  return x * x * x;
88  }
89 
90  ///\brief Logistic function/logistic function.
91  ///
92  ///Calculates the sigmoid function 1/(1+exp(-x)). The type must be arithmetic. For example
93  ///float,double,long double, int,... but no custom Type.
94  template<class T>
95  typename boost::enable_if<boost::is_arithmetic<T>, T>::type sigmoid(T x){
96  if(x < minExpInput<T>()) {
97  return 1;
98  }
99  if(x > maxExpInput<T>()) {
100  return 0;
101  }
102  return 1. / (1.+ std::exp(-x));
103  }
104 
105  ///\brief Thresholded exp function, over- and underflow safe.
106  ///
107  ///Replaces the value of exp(x) for numerical reasons by the a threshold value if it gets too large.
108  ///Use it only, if there is no other way to get the function stable!
109  ///
110  ///@param x the exponent
111  template<class T>
112  T safeExp(T x ){
113  if(x > maxExpInput<T>()){
114  //std::cout<<"warning, x too high"<<std::endl;
116  }
117  //Allow Koenig Lookup
118  using std::exp;
119  return exp(x);
120  }
121  ///\brief Thresholded log function, over- and underflow safe.
122  ///
123  ///Replaces the value of log(x) for numerical reasons by the a threshold value if it gets too low.
124  ///Use it only, if there is no other way to get the function stable!
125  ///@param x the exponent
126  template<class T>
127  T safeLog(T x)
128  {
129 
130  if(x> -std::numeric_limits<T>::epsilon() && x < std::numeric_limits<T>::epsilon()){
131  //std::cout<<"warning, x too low"<<std::endl;
132  return boost::math::sign(x)*std::numeric_limits<T>::min_exponent;
133  }
134 
135  //Allow Koenig Lookup
136  using std::log;
137  return log(x);
138  };
139 
140  ///\brief Numerically stable version of the function log(1+exp(x)).
141  ///
142  ///Numerically stable version of the function log(1+exp(x)).
143  ///This function is the integral of the famous sigmoid function.
144  ///The type must be arithmetic. For example
145  ///float,double,long double, int,... but no custom Type.
146  template<class T>
147  typename boost::enable_if<boost::is_arithmetic<T>, T>::type softPlus(T x){
148  if(x > maxExpInput<T>()){
149  return x;
150  }
151  if(x < minExpInput<T>()){
152  return 0;
153  }
154  return std::log(1+std::exp(x));
155  }
156 
157  ///\brief Numerically stable version of the function log(1+exp(x)). calculated with float precision to save some time
158  ///
159  ///Numerically stable version of the function log(1+exp(x)).
160  ///This function is the integral of the famous sigmoid function.
161  inline double softPlus(double x){
162  if(x > 15){
163  return x;
164  }
165  if(x < -17){
166  return 0;
167  }
168  return std::log(1.0f+std::exp(float(x)));
169  }
170 
171  ///brief lets x have the same sign as y.
172  ///
173  ///This is the famous well known copysign function from fortran.
174  template<class T>
175  T copySign(T x, T y){
176  using std::abs;
177  if(y > 0){
178  return abs(x);
179  }
180  return -abs(x);
181  }
182 
183  /// \brief Calculates the trigamma function -the derivative of the digamma function
184  ///
185  /// The trigamma function is defined as
186  /// \f[\Psi^{(1)}(x)= \frac{\partial^2}{\partial ^2} \log \gamma(x)\f]
187  ///
188  /// It has poles as x=-1,-2,-3,.. at which the function value is infinity.
189  /// For more information see http://mathworld.wolfram.com/TrigammaFunction.html
190  ///
191  /// We calculate it in two different ways. If the argument is small a taylor expansion
192  /// is used. otherwise the argument is transformed to a valu large enough such that
193  /// the asymptotic series holds which is then easy to calculate with high precision
194  inline double trigamma(double x){
195  double tiny = 1e-4;//threshold for taylor expansion validity around 0
196  double large = 15;//threshold above which the asymptotic formula is valid
197  //coefficients for taylor expansion - value at 0
198  double trigamma1 = 1.6449340668482264365; // pi^2/6 = Zeta(2)
199  double tetragamma1 = -2.404113806319188570799476; // -2 Zeta(3) */
200  //even Bernoulli coefficients B_2..B_10 for the asymptotic formula
201  double b2 = 1./6;
202  double b4 = -1./30;
203  double b6 = 1./42;
204  double b8 = -1./30;
205  double b10 = 5./66;
206  // singularities are at -1, -2, -3,...
207  if((x <= 0) && (std::floor(x) == x)) {
208  return -std::numeric_limits<double>::infinity();
209  }
210  //treating of general negative arguments
211  if(x < 0)
212  {
213  //use reflection relation: -psi^1(1-x)-psi^1(x)=pi* d/dz cot(pi*x)
214  //=> psi^1(x) = -psi^1(1-x)-pi* d/dz cot(pi*x)
215  double pi = boost::math::constants::pi<double>();
216  double s = std::sin(pi*x);
217  double ddz_cot_pi_x=-pi / (s * s);
218  return -trigamma(1-x) - pi * ddz_cot_pi_x;
219  }
220  // Use Taylor series if argument <= tiny
221  //for this make use of the entity trigamma(x)=trigamma(x+1)+1/x^2
222  //to move the origin to 1 and calculate the taylor expansion at 1
223  if(x <= tiny) {
224  return 1/(x*x) + trigamma1 + tetragamma1*x;
225  }
226 
227  //otherwise make use of the same entity
228  //for n steps until x+n > large, which is the threshold for the
229  // asymptotic formula to converge.
230  double result = 0;
231  while(x < large) {
232  result += 1/(x*x);
233  x++;
234  }
235  //now apply the asymptotic formula
236  //trigamma(x)=1/x+1/(2*x*x)+sum_k B_{2k}/z^{2k+1}
237  if(x >= large) {
238  double r = 1/(x*x);
239  double t = (b4 + r*(b6 + r*(b8 + r*b10)));
240  result += 0.5*r + (1 + r*(b2 + r*t))/x;
241  }
242  return result;
243  }
244 
245  /// \brief Calculates the tetragamma function - the derivative of the trigamma function
246  ///
247  /// The trigamma function is defined as
248  /// \f[\Psi^{(2)}(x)= \frac{\partial}{\partial x}\Psi^{(1)}(x)= \frac{\partial^3}{\partial ^3} \log \gamma(x)\f]
249  ///
250  /// The function is undefined for x=-1,-2,-3,... and an exception is thrown.
251  /// For more information see http://mathworld.wolfram.com/PolygammaFunction.html
252  ///
253  /// We calculate it in two different ways. If the argument is small a taylor expansion
254  /// is used. otherwise the argument is transformed to a valu large enough such that
255  /// the asymptotic series holds which is then easy to calculate with high precision
256  inline double tetragamma(double x)
257  {
258  double tiny = 1e-4;//threshold for use of taylor expansion if value is close to 0
259  double large = 18;//threshold above which the asymptotic formula is valid
260  //coefficients for taylor expansion at 1
261  double tetragamma1 = -2.404113806319188570799476; /* -2 Zeta(3) */
262  double pentagamma1 = 6.49393940226682914909602217; /* 6 Zeta(4) */
263  //even Bernoulli coefficients B_2..B_10 for the asymptotic formula
264  double b2 = 1./6;
265  double b4 = -1./30;
266  double b6 = 1./42;
267  double b8 = -1./30;
268  double b10 = 5./66;
269  // singularities are at -1, -2, -3,...
270  if((x <= 0) && (std::floor(x) == x)) {
271  throw SHARKEXCEPTION("[tetragamma] negative whole numbers ar enot allowed");
272  }
273  //treating of general negative arguments
274  if(x < 0)
275  {
276  //use derivative of reflection relation of trigamma:
277  //-psi^2(1-x)-psi^2(x)=pi* d^2/dz^2 cot(pi*x)
278  //=> psi^2(x) = -psi^2(1-x)-pi* d^2/dz^2 cot(pi*x)
279  //and d^2/dz^2 cot(pi*x)= -pi d/dz 1/sin(pi*x)^2
280  // = 2*pi^2 cos(pi*x)/sin(pi*x)^3
281  double pi = boost::math::constants::pi<double>();
282  double s_pi_x = std::sin(pi*x);
283  double c_pi_x = std::cos(pi*x);
284  double ddz2_cot_pi_x=2*sqr(pi)*c_pi_x/cube(s_pi_x);
285  return -tetragamma(1-x) - pi * ddz2_cot_pi_x;
286  }
287  // Use Taylor series if argument <= tiny
288  //for this make use of the entity tetragamma(x)=tetragamma(x+1)-2/x^3
289  //to move the origin to 1 and calculate the taylor expansion at 1
290  if(x <= tiny) {
291  return -2/cube(x) + tetragamma1 + pentagamma1*x;
292  }
293 
294  //otherwise make use of the same entity
295  //for n steps until x+n > large, which is the threshold for the
296  // asymptotic formula to converge.
297  double result = 0;
298  while(x < large) {
299  result -= 2/cube(x);
300  x++;
301  }
302  //now apply the asymptotic formula - the derivative of the asymptotic formula
303  // of the trigamma function
304  //tetragamma(x)=-1/x^2-1/(x^3)-sum_k (2k+1)*B_{2k}/z^{2(k+1)}
305 
306  if(x >= large) {
307  double r = 1/(x*x);
308  double t = (5*b4 + r*(7*b6 + r*(9*b8 + r*11*b10)));
309  result -= r/x + r*(1 + r*(3*b2 + r*t));
310  }
311  return result;
312 
313  }
314 
315 }
316 /** @}*/
317 
318 #endif