Geant4  10.02.p02
G4PolynomialPDF.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // -------------------------------------------------------------------
28 // GEANT4 Class file
29 //
30 //
31 // File name: G4PolynomialPDF
32 //
33 // Author: Jason Detwiler (jasondet@gmail.com)
34 //
35 // Creation date: Aug 2012
36 // -------------------------------------------------------------------
37 
38 #include "G4PolynomialPDF.hh"
39 #include "Randomize.hh"
40 //#include <time.h>
41 
42 using namespace std;
43 
45  G4double x1, G4double x2) :
46 fX1(x1), fX2(x2), fChanged(true), fTolerance(1.e-8)
47 {
48  if(coeffs != NULL) SetCoefficients(n, coeffs);
49  else if(n > 0) SetNCoefficients(n);
50 }
51 
53 {
54  while(i >= fCoefficients.size()) fCoefficients.push_back(0); /* Loop checking, 30-Oct-2015, G.Folger */
55  fCoefficients[i] = value;
56  fChanged = true;
57 }
58 
59 void G4PolynomialPDF::SetCoefficients(size_t nCoeffs,
60  const G4double* coefficients)
61 {
62  SetNCoefficients(nCoeffs);
63  for(size_t i=0; i<GetNCoefficients(); ++i) {
64  SetCoefficient(i, coefficients[i]);
65  }
66  fChanged = true;
67 }
68 
70 {
71  if(x2 <= x1) {
72  G4cout << "G4PolynomialPDF::SetDomain(): Invalide domain! "
73  << "(x1 = " << x1 << ", x2 = " << x2 << ")." << G4endl;
74  return;
75  }
76  fX1 = x1;
77  fX2 = x2;
78  fChanged = true;
79 }
80 
82 {
85  while(fCoefficients.size()) { /* Loop checking, 30-Oct-2015, G.Folger */
86  if(fCoefficients[fCoefficients.size()-1] == 0.0) fCoefficients.pop_back();
87  else break;
88  }
89 
90  G4double x1N = fX1, x2N = fX2;
91  G4double sum = 0;
92  for(size_t i=0; i<GetNCoefficients(); ++i) {
93  sum += GetCoefficient(i)*(x2N - x1N)/G4double(i+1);
94  x1N*=fX1;
95  x2N*=fX2;
96  }
97  if(sum <= 0) {
98  G4cout << "G4PolynomialPDF::Normalize() PDF has non-positive area: "
99  << sum << G4endl;
100  Dump();
101  return;
102  }
103 
104  for(size_t i=0; i<GetNCoefficients(); ++i) {
105  SetCoefficient(i, GetCoefficient(i)/sum);
106  }
107 }
108 
110 {
116  if(ddxPower < -1 || ddxPower > 2) {
117  G4cout << "G4PolynomialPDF::GetX(): ddxPower " << ddxPower
118  << " not implemented" << G4endl;
119  return 0;
120  }
121 
122  double f = 0; // return value
123  double xN = 1; // x to the power N
124  double x1N = 1; // endpoint x1 to the power N; only used by CDF
125  for(size_t i=0; i<=GetNCoefficients(); ++i) {
126  if(ddxPower == -1) { // CDF
127  if(i>0) f += GetCoefficient(i-1)*(xN - x1N)/i;
128  x1N *= fX1;
129  }
130  else if(ddxPower == 0 && i<GetNCoefficients()) f += GetCoefficient(i)*xN; // PDF
131  else if(ddxPower == 1) { // (d/dx) PDF
132  if(i<GetNCoefficients()-1) f += GetCoefficient(i+1)*xN*(i+1);
133  }
134  else if(ddxPower == 2) { // (d2/dx2) PDF
135  if(i<GetNCoefficients()-2) f += GetCoefficient(i+2)*xN*((i+2)*(i+1));
136  }
137  xN *= x;
138  }
139  return f;
140 }
141 
143 {
144  // ax2 + bx + c = 0
145  // p': 2ax + b = 0 -> = 0 at min: x_extreme = -b/2a
146 
147  if(x1 < fX1 || x2 > fX2 || x2 < x1) {
148  G4cout << "G4PolynomialPDF::HasNegativeMinimum(): Invalid range "
149  << x1 << " - " << x2 << G4endl;
150  return false;
151  }
152 
153  // If flat, then check anywhere.
154  if(GetNCoefficients() == 1) return (Evaluate(x1) < -fTolerance);
155 
156  // If linear, or if quadratic with negative second derivative,
157  // just check the endpoints
158  if(GetNCoefficients() == 2 ||
159  (GetNCoefficients() == 3 && GetCoefficient(2) < 0)) {
160  return (Evaluate(x1) < -fTolerance) || (Evaluate(x2) < -fTolerance);
161  }
162 
163  // If quadratic and second dervative is positive, check at the mininum
164  if(GetNCoefficients() == 3) {
165  G4double xMin = -GetCoefficient(1)*0.5/GetCoefficient(2);
166  if(xMin < x1) xMin = x1;
167  if(xMin > x2) xMin = x2;
168  return Evaluate(xMin) < -fTolerance;
169  }
170 
171  // Higher-order polynomials: consider any extremum between x1 and x2. If none
172  // are found, check the endpoints.
173  G4double extremum = GetX(0, x1, x2, 1);
174  if(Evaluate(extremum) < -fTolerance) return true;
175  else if(extremum <= x1+(x2-x1)*fTolerance ||
176  extremum >= x2-(x2-x1)*fTolerance) return false;
177  else return
178  HasNegativeMinimum(x1, extremum) || HasNegativeMinimum(extremum, x2);
179 }
180 
182 {
183  if(fChanged) {
184  Normalize();
185  if(HasNegativeMinimum(fX1, fX2)) {
186  G4cout << "G4PolynomialPDF::GetRandomX(): PDF has negative values, returning 0..." << G4endl;
187  return 0;
188  }
189  fChanged = false;
190  }
191  return EvalInverseCDF(G4UniformRand());
192 }
193 
195  G4int ddxPower, G4double guess, G4bool bisect)
196 {
203 
204  // input range checking
205  if(GetNCoefficients() == 0) {
206  G4cout << "G4PolynomialPDF::GetX(): no PDF defined!" << G4endl;
207  return x2;
208  }
209  if(ddxPower < -1 || ddxPower > 1) {
210  G4cout << "G4PolynomialPDF::GetX(): ddxPower " << ddxPower
211  << " not implemented" << G4endl;
212  return x2;
213  }
214  if(ddxPower == -1 && (p<0 || p>1)) {
215  G4cout << "G4PolynomialPDF::GetX(): p is out of range" << G4endl;
216  return fX2;
217  }
218 
219  // check limits
220  if(x2 <= x1 || x1 < fX1 || x2 > fX2) {
221  G4cout << "G4PolynomialPDF::GetX(): domain must have fX1 <= x1 < x2 <= fX2. "
222  << "You sent x1 = " << x1 << ", x2 = " << x2 << "." << G4endl;
223  return x2;
224  }
225 
226  // Return x2 for flat lines
227  if((ddxPower == 0 && GetNCoefficients() == 1) ||
228  (ddxPower == 1 && GetNCoefficients() == 2)) return x2;
229 
230  // Solve p = mx + b -> x = (p-b)/m for linear functions
231  if((ddxPower == -1 && GetNCoefficients() == 1) ||
232  (ddxPower == 0 && GetNCoefficients() == 2) ||
233  (ddxPower == 1 && GetNCoefficients() == 3)) {
234  G4double b = (ddxPower > -1) ? GetCoefficient(ddxPower) : -GetCoefficient(0)*fX1;
235  G4double slope = GetCoefficient(ddxPower+1);
236  if(ddxPower == 1) slope *= 2.;
237  G4double value = (p-b)/slope;
238  if(value < x1) {
239  return x1;
240  }
241  else if(value > x2) {
242  return x2;
243  }
244  else {
245  return value;
246  }
247  }
248 
249  // Solve quadratic equation for f-p=0 when f is quadratic
250  if((ddxPower == -1 && GetNCoefficients() == 2) ||
251  (ddxPower == 0 && GetNCoefficients() == 3) ||
252  (ddxPower == 1 && GetNCoefficients() == 4)) {
253  G4double c = -p + ((ddxPower > -1) ? GetCoefficient(ddxPower) : 0);
254  if(ddxPower == -1) c -= (GetCoefficient(0) + GetCoefficient(1)/2.*fX1)*fX1;
255  G4double b = GetCoefficient(ddxPower+1);
256  if(ddxPower == 1) b *= 2.;
257  G4double a = GetCoefficient(ddxPower+2);
258  if(ddxPower == 1) a *= 3;
259  else if(ddxPower == -1) a *= 0.5;
260  double sqrtFactor = b*b - 4.*a*c;
261  if(sqrtFactor < 0) return x2; // quadratic equation has no solution (p not in range of f)
262  sqrtFactor = sqrt(sqrtFactor)/2./fabs(a);
263  G4double valueMinus = -b/2./a - sqrtFactor;
264  if(valueMinus >= x1 && valueMinus <= x2) return valueMinus;
265  else if(valueMinus > x2) return x2;
266  G4double valuePlus = -b/2./a + sqrtFactor;
267  if(valuePlus >= x1 && valuePlus <= x2) return valuePlus;
268  else if(valuePlus < x1) return x2;
269  return (x1-valueMinus <= valuePlus-x2) ? x1 : x2;
270  }
271 
272  // f is non-trivial, so use Newton-Raphson
273  // start in the middle if no good guess is provided
274  if(guess < x1 || guess > x2) guess = (x2+x1)*0.5;
275  G4double lastChange = 1;
276  size_t iterations = 0;
277  while(fabs(lastChange) > fTolerance) { /* Loop checking, 02.11.2015, A.Ribon */
278  // calculate f and f' simultaneously
279  G4double f = -p;
280  G4double dfdx = 0;
281  G4double xN = 1;
282  G4double x1N = 1; // only used by CDF
283  for(size_t i=0; i<=GetNCoefficients(); ++i) {
284  if(ddxPower == -1) { // CDF
285  if(i>0) f += GetCoefficient(i-1)*(xN - x1N)/G4double(i);
286  if(i<GetNCoefficients()) dfdx += GetCoefficient(i)*xN;
287  x1N *= fX1;
288  }
289  else if(ddxPower == 0) { // PDF
290  if(i<GetNCoefficients()) f += GetCoefficient(i)*xN;
291  if(i+1<GetNCoefficients()) dfdx += GetCoefficient(i+1)*xN*G4double(i+1);
292  }
293  else { // ddxPower == 1: (d/dx) PDF
294  if(i+1<GetNCoefficients()) f += GetCoefficient(i+1)*xN*G4double(i+1);
295  if(i+2<GetNCoefficients()) dfdx += GetCoefficient(i+2)*xN*G4double(i+2);
296  }
297  xN *= guess;
298  }
299  if(f == 0) return guess;
300  if(dfdx == 0) {
301  G4cout << "G4PolynomialPDF::GetX(): got f != 0 but slope = 0 for ddxPower = "
302  << ddxPower << G4endl;
303  return x2;
304  }
305  lastChange = - f/dfdx;
306 
307  if(guess + lastChange < x1) {
308  lastChange = x1 - guess;
309  } else if(guess + lastChange > x2) {
310  lastChange = x2 - guess;
311  }
312 
313  guess += lastChange;
314  lastChange /= (fX2-fX1);
315 
316  ++iterations;
317  if(iterations > 50) {
318  if(p!=0) {
319  G4cout << "G4PolynomialPDF::GetX(): got stuck searching for " << p
320  << " between " << x1 << " and " << x2 << " with ddxPower = "
321  << ddxPower
322  << ". Last guess was " << guess << "." << G4endl;
323  }
324  if(ddxPower==-1 && bisect) {
325  G4cout << "Bisceting and trying again..." << G4endl;
326  return Bisect(p, x1, x2);
327  }
328  else return guess;
329  }
330  }
331  return guess;
332 }
333 
335  // Bisect to get 1% precision, then use Newton-Raphson
336  G4double z = (x2 + x1)/2.0; // [x1 z x2]
337  if((x2 - x1)/(fX2 - fX1) < 0.01) return GetX(p, fX1, fX2, -1, z, false);
338  G4double fz = Evaluate(z, -1) - p;
339  if(fz < 0) return Bisect(p, z, x2); // [z x2]
340  return Bisect(p, x1, z); // [x1 z]
341 }
342 
344 {
345  G4cout << "G4PolynomialPDF::Dump() - PDF(x) = ";
346  for(size_t i=0; i<GetNCoefficients(); i++) {
347  if(i>0) G4cout << " + ";
348  G4cout << GetCoefficient(i);
349  if(i>0) G4cout << "*x";
350  if(i>1) G4cout << "^" << i;
351  }
352  G4cout << G4endl;
353  G4cout << "G4PolynomialPDF::Dump() - Interval: " << fX1 << " <= x < "
354  << fX2 << G4endl;
355 }
356 
357 
std::vector< G4double > fCoefficients
void SetCoefficients(const std::vector< G4double > &v)
G4double z
Definition: TRTMaterials.hh:39
G4double GetRandomX()
void SetNCoefficients(size_t n)
G4PolynomialPDF(size_t n=0, const double *coeffs=nullptr, G4double x1=0, G4double x2=1)
G4double a
Definition: TRTMaterials.hh:39
G4double GetX(G4double p, G4double x1, G4double x2, G4int ddxPower=0, G4double guess=1.e99, G4bool bisect=true)
int G4int
Definition: G4Types.hh:78
G4double Bisect(G4double p, G4double x1, G4double x2)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4double GetCoefficient(size_t i) const
void SetDomain(G4double x1, G4double x2)
bool G4bool
Definition: G4Types.hh:79
const G4int n
void SetCoefficient(size_t i, G4double value)
G4bool HasNegativeMinimum(G4double x1, G4double x2)
const G4double x[NPOINTSGL]
size_t GetNCoefficients() const
#define G4endl
Definition: G4ios.hh:61
G4double EvalInverseCDF(G4double p)
G4double Evaluate(G4double x, G4int ddxPower=0)
double G4double
Definition: G4Types.hh:76