CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

flatToGaussian.cc
Go to the documentation of this file.
1 // $Id: flatToGaussian.cc,v 1.4 2003/08/13 20:00:12 garren Exp $
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- flatToGaussian ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 
10 // Contains the methods that depend on the 30K-footpring gaussTables.cdat.
11 //
12 // flatToGaussian (double x)
13 // inverseErf (double x)
14 // erf (double x)
15 
16 // =======================================================================
17 // M Fischler - Created 1/25/00.
18 //
19 // =======================================================================
20 
21 #include "CLHEP/Random/Stat.h"
22 #include "CLHEP/Random/defs.h"
23 #include "CLHEP/Units/PhysicalConstants.h"
24 #include <iostream>
25 #include <cmath>
26 
27 namespace CLHEP {
28 
29 double transformSmall (double r);
30 
31 
32 //
33 // Table of errInts, for use with transform(r) and quickTransform(r)
34 //
35 
36 #ifdef Traces
37 #define Trace1
38 #define Trace2
39 #define Trace3
40 #endif
41 
42 // Since all these are this is static to this compilation unit only, the
43 // info is establised a priori and not at each invocation.
44 
45 // The main data is of course the gaussTables table; the rest is all
46 // bookkeeping to know what the tables mean.
47 
48 #define Table0size 200
49 #define Table1size 250
50 #define Table2size 200
51 #define Table3size 250
52 #define Table4size 1000
53 #define TableSize (Table0size+Table1size+Table2size+Table3size+Table4size)
54 
55 static const int Tsizes[5] = { Table0size,
56  Table1size,
57  Table2size,
58  Table3size,
59  Table4size };
60 
61 #define Table0step (2.0E-13)
62 #define Table1step (4.0E-11)
63 #define Table2step (1.0E-8)
64 #define Table3step (2.0E-6)
65 #define Table4step (5.0E-4)
66 
67 static const double Tsteps[5] = { Table0step,
68  Table1step,
69  Table2step,
70  Table3step,
71  Table4step };
72 
73 #define Table0offset 0
74 #define Table1offset (2*(Table0size))
75 #define Table2offset (2*(Table0size + Table1size))
76 #define Table3offset (2*(Table0size + Table1size + Table2size))
77 #define Table4offset (2*(Table0size + Table1size + Table2size + Table3size))
78 
79 static const int Toffsets[5] = { Table0offset,
83  Table4offset };
84 
85  // Here comes the big (30K bytes) table, kept in a file ---
86 
87 static const double gaussTables [2*TableSize] = {
88 #include "gaussTables.cdat"
89 };
90 
91 
92 double HepStat::flatToGaussian (double r) {
93 
94  double sign = +1.0; // We always compute a negative number of
95  // sigmas. For r > 0 we will multiply by
96  // sign = -1 to return a positive number.
97 #ifdef Trace1
98 std::cout << "r = " << r << "\n";
99 #endif
100 
101  if ( r > .5 ) {
102  r = 1-r;
103  sign = -1.0;
104 #ifdef Trace1
105 std::cout << "r = " << r << "sign negative \n";
106 #endif
107  } else if ( r == .5 ) {
108  return 0.0;
109  }
110 
111  // Find a pointer to the proper table entries, along with the fraction
112  // of the way in the relevant bin dx and the bin size h.
113 
114  // Optimize for the case of table 4 by testing for that first.
115  // By removing that case from the for loop below, we save not only
116  // several table lookups, but also several index calculations that
117  // now become (compile-time) constants.
118  //
119  // Past the case of table 4, we need not be as concerned about speed since
120  // this will happen only .1% of the time.
121 
122  const double* tptr = 0;
123  double dx = 0;
124  double h = 0;
125 
126  // The following big if block will locate tptr, h, and dx from whichever
127  // table is applicable:
128 
129  register int index;
130 
131  if ( r >= Table4step ) {
132 
133  index = int((Table4size<<1) * r); // 1 to Table4size-1
134  if (index <= 0) index = 1; // in case of rounding problem
135  if (index >= Table4size) index = Table4size-1;
136  dx = (Table4size<<1) * r - index; // fraction of way to next bin
137  h = Table4step;
138 #ifdef Trace2
139 std::cout << "index = " << index << " dx = " << dx << " h = " << h << "\n";
140 #endif
141  index = (index<<1) + (Table4offset-2);
142  // at r = table4step+eps, index refers to the start of table 4
143  // and at r = .5 - eps, index refers to the next-to-last entry:
144  // remember, the table has two numbers per actual entry.
145 #ifdef Trace2
146 std::cout << "offset index = " << index << "\n";
147 #endif
148 
149  tptr = &gaussTables [index];
150 
151  } else if (r < Tsteps[0]) {
152 
153  // If r is so small none of the tables apply, use the asymptotic formula.
154  return (sign * transformSmall (r));
155 
156  } else {
157 
158  for ( int tableN = 3; tableN >= 0; tableN-- ) {
159  if ( r < Tsteps[tableN] ) {continue;} // can't happen when tableN==0
160 #ifdef Trace2
161 std::cout << "Using table " << tableN << "\n";
162 #endif
163  double step = Tsteps[tableN];
164  index = int(r/step); // 1 to TableNsize-1
165  // The following two tests should probably never be true, but
166  // roundoff might cause index to be outside its proper range.
167  // In such a case, the interpolation still makes sense, but we
168  // need to take care that tptr points to valid data in the right table.
169  if (index == 0) index = 1;
170  if (index >= Tsizes[tableN]) index = Tsizes[tableN] - 1;
171  dx = r/step - index; // fraction of way to next bin
172  h = step;
173 #ifdef Trace2
174 std::cout << "index = " << index << " dx = " << dx << " h = " << h << "\n";
175 #endif
176  index = (index<<1) + Toffsets[tableN] - 2;
177  tptr = &gaussTables [index];
178  break;
179  } // end of the for loop which finds tptr, dx and h in one of the tables
180 
181  // The code can only get to here by a break statement, having set dx etc.
182  // It not get to here without going through one of the breaks,
183  // because before the for loop we test for the case of r < Tsteps[0].
184 
185  } // End of the big if block.
186 
187  // At the end of this if block, we have tptr, dx and h -- and if r is less
188  // than the smallest step, we have already returned the proper answer.
189 
190  double y0 = *tptr++;
191  double d0 = *tptr++;
192  double y1 = *tptr++;
193  double d1 = *tptr;
194 
195 #ifdef Trace3
196 std::cout << "y0: " << y0 << " y1: " << y1 << " d0: " << d0 << " d1: " << d1 << "\n";
197 #endif
198 
199  double x2 = dx * dx;
200  double oneMinusX = 1 - dx;
201  double oneMinusX2 = oneMinusX * oneMinusX;
202 
203  double f0 = (2. * dx + 1.) * oneMinusX2;
204  double f1 = (3. - 2. * dx) * x2;
205  double g0 = h * dx * oneMinusX2;
206  double g1 = - h * oneMinusX * x2;
207 
208 #ifdef Trace3
209 std::cout << "f0: " << f0 << " f1: " << f1 << " g0: " << g0 << " g1: " << g1 << "\n";
210 #endif
211 
212  double answer = f0 * y0 + f1 * y1 + g0 * d0 + g1 * d1;
213 
214 #ifdef Trace1
215 std::cout << "variate is: " << sign*answer << "\n";
216 #endif
217 
218  return sign * answer;
219 
220 } // flatToGaussian
221 
222 
223 double transformSmall (double r) {
224 
225  // Solve for -v in the asymtotic formula
226  //
227  // errInt (-v) = std::exp(-v*v/2) 1 1*3 1*3*5
228  // ------------ * (1 - ---- + ---- - ----- + ... )
229  // v*std::sqrt(2*pi) v**2 v**4 v**6
230 
231  // The value of r (=errInt(-v)) supplied is going to less than 2.0E-13,
232  // which is such that v < -7.25. Since the value of r is meaningful only
233  // to an absolute error of 1E-16 (double precision accuracy for a number
234  // which on the high side could be of the form 1-epsilon), computing
235  // v to more than 3-4 digits of accuracy is suspect; however, to ensure
236  // smoothness with the table generator (which uses quite a few terms) we
237  // also use terms up to 1*3*5* ... *13/v**14, and insist on accuracy of
238  // solution at the level of 1.0e-7.
239 
240  // This routine is called less than one time in a trillion firings, so
241  // speed is of no concern. As a matter of technique, we terminate the
242  // iterations in case they would be infinite, but this should not happen.
243 
244  double eps = 1.0e-7;
245  double guess = 7.5;
246  double v;
247 
248  for ( int i = 1; i < 50; i++ ) {
249  double vn2 = 1.0/(guess*guess);
250  double s1 = -13*11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2*vn2;
251  s1 += 11*9*7*5*3 * vn2*vn2*vn2*vn2*vn2*vn2;
252  s1 += -9*7*5*3 * vn2*vn2*vn2*vn2*vn2;
253  s1 += 7*5*3 * vn2*vn2*vn2*vn2;
254  s1 += -5*3 * vn2*vn2*vn2;
255  s1 += 3 * vn2*vn2 - vn2 + 1.0;
256  v = std::sqrt ( 2.0 * std::log ( s1 / (r*guess*std::sqrt(CLHEP::twopi)) ) );
257  if ( std::abs(v-guess) < eps ) break;
258  guess = v;
259  }
260 
261  return -v;
262 
263 } // transformSmall()
264 
265 
266 double HepStat::inverseErf (double t) {
267 
268  // This uses erf(x) = 2*gaussCDF(std::sqrt(2)*x) - 1
269 
270  return std::sqrt(0.5) * flatToGaussian(.5*(t+1));
271 
272 }
273 
274 
275 double HepStat::erf (double x) {
276 
277 // Math:
278 //
279 // For any given x we can "quickly" find t0 = erfQ (x) = erf(x) + epsilon.
280 //
281 // Then we can find x1 = inverseErf (t0) = inverseErf (erf(x)+epsilon)
282 //
283 // Expanding in the small epsion,
284 //
285 // x1 = x + epsilon * [deriv(inverseErf(u),u) at u = t0] + O(epsilon**2)
286 //
287 // so epsilon is (x1-x) / [deriv(inverseErf(u),u) at u = t0] + O(epsilon**2)
288 //
289 // But the derivative of an inverse function is one over the derivative of the
290 // function, so
291 // epsilon = (x1-x) * [deriv(erf(v),v) at v = x] + O(epsilon**2)
292 //
293 // And the definition of the erf integral makes that derivative trivial.
294 // Ultimately,
295 //
296 // erf(x) = erfQ(x) - (inverseErf(erfQ(x))-x) * std::exp(-x**2) * 2/std::sqrt(CLHEP::pi)
297 //
298 
299  double t0 = erfQ(x);
300  double deriv = std::exp(-x*x) * (2.0 / std::sqrt(CLHEP::pi));
301 
302  return t0 - (inverseErf (t0) - x) * deriv;
303 
304 }
305 
306 
307 } // namespace CLHEP
308 
CLHEP::HepStat::erf
static double erf(double x)
Definition: flatToGaussian.cc:275
Table1size
#define Table1size
Definition: flatToGaussian.cc:49
Table2step
#define Table2step
Definition: flatToGaussian.cc:63
Table0offset
#define Table0offset
Definition: flatToGaussian.cc:73
CLHEP::HepStat::erfQ
static double erfQ(double x)
Definition: erfQ.cc:25
CLHEP::transformSmall
double transformSmall(double r)
Definition: flatToGaussian.cc:223
CLHEP::HepStat::flatToGaussian
static double flatToGaussian(double r)
Definition: flatToGaussian.cc:92
Table1step
#define Table1step
Definition: flatToGaussian.cc:62
TableSize
#define TableSize
Definition: flatToGaussian.cc:53
Table0step
#define Table0step
Definition: flatToGaussian.cc:61
Table0size
#define Table0size
Definition: flatToGaussian.cc:48
Table3size
#define Table3size
Definition: flatToGaussian.cc:51
f1
void(* f1)()
Definition: testCategories.cc:151
Table3offset
#define Table3offset
Definition: flatToGaussian.cc:76
CLHEP
Definition: ClhepVersion.h:13
Table2size
#define Table2size
Definition: flatToGaussian.cc:50
v
they are gone ZOOM Features Discontinued The following features of the ZOOM package were felt to be extreme overkill These have been after checking that no existing user code was utilizing as in SpaceVector v
Definition: keyMergeIssues.doc:324
Table3step
#define Table3step
Definition: flatToGaussian.cc:64
Table4size
#define Table4size
Definition: flatToGaussian.cc:52
CLHEP::HepStat::inverseErf
static double inverseErf(double t)
Definition: flatToGaussian.cc:266
Table4offset
#define Table4offset
Definition: flatToGaussian.cc:77
i
long i
Definition: JamesRandomSeeding.txt:27
x
any side effects of that construction would occur twice The semantics of throw x
Definition: whyZMthrowRethrows.txt:37
Table2offset
#define Table2offset
Definition: flatToGaussian.cc:75
Table4step
#define Table4step
Definition: flatToGaussian.cc:65
Table1offset
#define Table1offset
Definition: flatToGaussian.cc:74