MercuryDPM  Beta
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ExtendedMath.cc
Go to the documentation of this file.
1 //Copyright (c) 2013-2014, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 
26 #include "ExtendedMath.h"
27 
28 #include <cmath>
29 #include <assert.h>
31 //#include <limits>
32 //#include <iostream>
33 //#include <sys/stat.h>
34 #include <iomanip>
35 //#include <cmath>
36 //#include <sstream>
37 //#include <cstdlib>
38 //#include <limits>
39 //#include <string>
40 
48 {
49  const Mdouble ep = 1e-5;
50 
51  if (gamma_in > 1.0 + ep)
52  {
53  return ((gamma_in - 1) * gamma(gamma_in - 1));
54  }
55  else
56  {
57 
58  if ((gamma_in - ep < 1.0) && (gamma_in + ep > 1.0))
59  return 1.0;
60  else if ((gamma_in - ep < 0.5) && (gamma_in + ep > 0.5))
61  return constants::sqrt_pi;
62  else
63  return std::numeric_limits<Mdouble>::quiet_NaN();
64  }
65 } //end func gamma
66 
70 Mdouble mathsFunc::chi_squared(const Mdouble x, const unsigned int k)
71 {
72 
73  Mdouble prefactor = pow(2, k / 2.0) * gamma(k / 2.0);
74  Mdouble mainfactor = pow(x, k / 2.0 - 1) * exp(x / -2.0);
75 
76  return mainfactor / prefactor;
77 
78 }
79 
86 Mdouble mathsFunc::chi_squared_prob(const Mdouble x_max, const unsigned int k)
87 {
88 
89 //The current value was picked by tried were it stopped effect the 4 d.p.
90  const int num_steps_per_unit = 100;
91  Mdouble sum = 0;
92  Mdouble x = 0;
93  long int num_steps = static_cast<int>(num_steps_per_unit * x_max);
94 //Use trapezional rule, but ignoring the ends
95  for (int i = 0; i < num_steps; i++)
96  {
97  x = x_max / num_steps * (i + 0.5);
98  sum = sum + chi_squared(x, k);
99  }
100  return 1.0 - sum * x_max / num_steps;
101 
102 }
103 
104 Mdouble mathsFunc::goldenSectionSearch(Mdouble(*function)(const Mdouble), Mdouble min, Mdouble cur, Mdouble max, Mdouble endCondition, Mdouble curVal)
105 {
106  if (std::abs(max - min) < endCondition)
107  {
108  return 0.5 * (min + max);
109  }
110  std::cout << "Min=" << min << " Max=" << max << " diff=" << max - min << std::endl;
111  double resphi = 2 - 0.5 * (1 + std::sqrt(5));
112  double x;
113  if (max - cur > cur - min)
114  {
115  x = cur + resphi * (max - cur);
116  }
117  else
118  {
119  x = cur - resphi * (cur - min);
120  }
121  if(std::isnan(curVal))
122  curVal = function(cur);
123  double xVal = function(x);
124  if (xVal < curVal)
125  {
126  if (max - cur > cur - min)
127  {
128  return goldenSectionSearch(function, cur, x, max, endCondition, xVal);
129  }
130  else
131  {
132  return goldenSectionSearch(function, min, x, cur, endCondition, xVal);
133  }
134  }
135  else
136  {
137  if (max - cur > cur - min)
138  {
139  return goldenSectionSearch(function, min, cur, x, endCondition, curVal);
140  }
141  else
142  {
143  return goldenSectionSearch(function, x, cur, max, endCondition, curVal);
144  }
145  }
146 }
147 
148 bool mathsFunc::isEqual(Mdouble v1, Mdouble v2, double absError)
149 {
150  return std::abs(v1 - v2) <= absError;
151 }
152 
153 
154 bool mathsFunc::isEqual(Vec3D v1, Vec3D v2, double absError)
155 {
156  return isEqual(v1.X,v2.X,absError)&&isEqual(v1.Y,v2.Y,absError)&&isEqual(v1.Z,v2.Z,absError);
157 }
158 
160 {
161  const Mdouble *p = coef;
162  Mdouble b0 = *p++;
163  Mdouble b1 = 0, b2;
164  int i = N - 1;
165 
166  assert(i > 0);
167  do
168  {
169  b2 = b1;
170  b1 = b0;
171  b0 = x * b1 - b2 + *p++;
172  }
173  while ( --i);
174 
175  return (0.5 * (b0 - b2));
176 }
177 
179 {
180  // Cooefficients for [0..8]
181  const Mdouble A[] =
182  {
183  -4.415341646479339379501E-18,
184  3.330794518822238097831E-17,
185  -2.431279846547954693591E-16,
186  1.715391285555133030611E-15,
187  -1.168533287799345168081E-14,
188  7.676185498604935616881E-14,
189  -4.856446783111929460901E-13,
190  2.955052663129639834611E-12,
191  -1.726826291441555707231E-11,
192  9.675809035373236912241E-11,
193  -5.189795601635262906661E-10,
194  2.659823724682386650351E-9,
195  -1.300025009986248042121E-8,
196  6.046995022541918949321E-8,
197  -2.670793853940611733911E-7,
198  1.117387539120103718151E-6,
199  -4.416738358458750563591E-6,
200  1.644844807072889708931E-5,
201  -5.754195010082103703981E-5,
202  1.885028850958416557291E-4,
203  -5.763755745385823658851E-4,
204  1.639475616941335798421E-3,
205  -4.324309995050575944301E-3,
206  1.054646039459499831831E-2,
207  -2.373741480589946881561E-2,
208  4.930528423967070848781E-2,
209  -9.490109704804764442101E-2,
210  1.716209015222087753491E-1,
211  -3.046826723431983986831E-1,
212  6.767952744094760849951E-1
213  };
214 
215  // Cooefficients for [8..infinity]
216  const Mdouble B[] =
217  {
218  -7.233180487874753954561E-18,
219  -4.830504485944182071261E-18,
220  4.465621420296759999011E-17,
221  3.461222867697461093101E-17,
222  -2.827623980516583484941E-16,
223  -3.425485619677219134621E-16,
224  1.772560133056526383601E-15,
225  3.811680669352622420751E-15,
226  -9.554846698828307648701E-15,
227  -4.150569347287222086631E-14,
228  1.540086217521409826911E-14,
229  3.852778382742142701141E-13,
230  7.180124451383666233671E-13,
231  -1.794178531506806117781E-12,
232  -1.321581184044771311881E-11,
233  -3.149916527963241364541E-11,
234  1.188914710784643834241E-11,
235  4.940602388224969589101E-10,
236  3.396232025708386345151E-9,
237  2.266668990498178064591E-8,
238  2.048918589469063741831E-7,
239  2.891370520834756482971E-6,
240  6.889758346916823984261E-5,
241  3.369116478255694089901E-3,
242  8.044904110141088316081E-1
243  };
244 
245  if (x < 0)
246  x = -x;
247 
248  if (x <= 8.0)
249  {
250  Mdouble y = (x / 2.0) - 2.0;
251  return (chebyshev(y, A, 30));
252  }
253 
254  return (chebyshev(32.0 / x - 2.0, B, 25) / sqrt(x));
255 }
256 
258 {
259  if (x < 0)
260  x = -x;
261  return exp(x) * I0_exp(x);
262 }
Mdouble chebyshev(Mdouble x, const Mdouble coef[], int N)
Mdouble X
the vector components
Definition: Vector.h:52
Mdouble chi_squared_prob(const Mdouble x, const unsigned int k)
This is the function which actually gives the probability back using a chi squared test...
Definition: ExtendedMath.cc:86
Mdouble I0_exp(Mdouble x)
const Mdouble sqrt_pi
Definition: ExtendedMath.h:43
double Mdouble
bool isEqual(Mdouble v1, Mdouble v2, double absError)
Compares the difference of two Mdouble with an absolute error, useful in UnitTests.
Mdouble Y
Definition: Vector.h:52
double goldenSectionSearch(double(*function)(const double), double min, double cur, double max, double endCondition, double curVal=std::numeric_limits< Mdouble >::quiet_NaN())
This function performs a golden section search to find the location of the minimum of a function...
Mdouble gamma(Mdouble gamma_in)
This is the gamma function returns the true value for the half integer value.
Definition: ExtendedMath.cc:47
Mdouble chi_squared(const Mdouble x, const unsigned int k)
This is a chi_squared function return the value x and degrees of freedom k.
Definition: ExtendedMath.cc:70
Mdouble I0(Mdouble x)
Implementation of a 3D vector (by Vitaliy).
Definition: Vector.h:45
Mdouble Z
Definition: Vector.h:52