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

RandMultiGauss.cc
Go to the documentation of this file.
1 // $Id: RandMultiGauss.cc,v 1.3 2003/08/13 20:00:13 garren Exp $
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- RandMultiGauss ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 
10 // =======================================================================
11 // Mark Fischler - Created: 17th September 1998
12 // =======================================================================
13 
14 // Some theory about how to get the Multivariate Gaussian from a bunch
15 // of Gaussian deviates. For the purpose of this discussion, take mu = 0.
16 //
17 // We want an n-vector x with distribution (See table 28.1 of Review of PP)
18 //
19 // exp ( .5 * x.T() * S.inverse() * x )
20 // f(x;S) = ------------------------------------
21 // sqrt ( (2*pi)^n * S.det() )
22 //
23 // Suppose S = U * D * U.T() with U orthogonal ( U*U.T() = 1 ) and D diagonal.
24 // Consider a random n-vector y such that each element y(i)is distributed as a
25 // Gaussian with sigma = sqrt(D(i,i)). Then the distribution of y is the
26 // product of n Gaussains which can be written as
27 //
28 // exp ( .5 * y.T() * D.inverse() * y )
29 // f(y;D) = ------------------------------------
30 // sqrt ( (2*pi)^n * D.det() )
31 //
32 // Now take an n-vector x = U * y (or y = U.T() * x ). Then substituting,
33 //
34 // exp ( .5 * x * U * D.inverse() U.T() * x )
35 // f(x;D,U) = ------------------------------------------
36 // sqrt ( (2*pi)^n * D.det() )
37 //
38 // and this simplifies to the desired f(x;S) when we note that
39 // D.det() = S.det() and U * D.inverse() * U.T() = S.inverse()
40 //
41 // So the strategy is to diagonalize S (finding U and D), form y with each
42 // element a Gaussian random with sigma of sqrt(D(i,i)), and form x = U*y.
43 // (It turns out that the D information needed is the sigmas.)
44 // Since for moderate or large n the work needed to diagonalize S can be much
45 // greater than the work to generate n Gaussian variates, we save the U and
46 // sigmas for both the default S and the latest S value provided.
47 
50 #include <cmath> // for log()
51 
52 namespace CLHEP {
53 
54 // ------------
55 // Constructors
56 // ------------
57 
59  const HepVector & mu,
60  const HepSymMatrix & S )
61  : localEngine(&anEngine),
62  deleteEngine(false),
63  set(false),
64  nextGaussian(0.0)
65 {
66  if (S.num_row() != mu.num_row()) {
67  std::cerr << "In constructor of RandMultiGauss distribution: \n" <<
68  " Dimension of mu (" << mu.num_row() <<
69  ") does not match dimension of S (" << S.num_row() << ")\n";
70  std::cerr << "---Exiting to System\n";
71  exit(1);
72  }
73  defaultMu = mu;
74  defaultSigmas = HepVector(S.num_row());
75  prepareUsigmas (S, defaultU, defaultSigmas);
76 }
77 
79  const HepVector & mu,
80  const HepSymMatrix & S )
81  : localEngine(anEngine),
82  deleteEngine(true),
83  set(false),
84  nextGaussian(0.0)
85 {
86  if (S.num_row() != mu.num_row()) {
87  std::cerr << "In constructor of RandMultiGauss distribution: \n" <<
88  " Dimension of mu (" << mu.num_row() <<
89  ") does not match dimension of S (" << S.num_row() << ")\n";
90  std::cerr << "---Exiting to System\n";
91  exit(1);
92  }
93  defaultMu = mu;
94  defaultSigmas = HepVector(S.num_row());
95  prepareUsigmas (S, defaultU, defaultSigmas);
96 }
97 
99  : localEngine(&anEngine),
100  deleteEngine(false),
101  set(false),
102  nextGaussian(0.0)
103 {
104  defaultMu = HepVector(2,0);
105  defaultU = HepMatrix(2,1);
106  defaultSigmas = HepVector(2);
107  defaultSigmas(1) = 1.;
108  defaultSigmas(2) = 1.;
109 }
110 
112  : localEngine(anEngine),
113  deleteEngine(true),
114  set(false),
115  nextGaussian(0.0)
116 {
117  defaultMu = HepVector(2,0);
118  defaultU = HepMatrix(2,1);
119  defaultSigmas = HepVector(2);
120  defaultSigmas(1) = 1.;
121  defaultSigmas(2) = 1.;
122 }
123 
125  if ( deleteEngine ) delete localEngine;
126 }
127 
128 // ----------------------------
129 // prepareUsigmas()
130 // ----------------------------
131 
132 void RandMultiGauss::prepareUsigmas( const HepSymMatrix & S,
133  HepMatrix & U,
134  HepVector & sigmas ) {
135 
136  HepSymMatrix tempS ( S ); // Since diagonalize does not take a const s
137  // we have to copy S.
138 
139  U = diagonalize ( &tempS ); // S = U Sdiag U.T()
140  HepSymMatrix D = S.similarityT(U); // D = U.T() S U = Sdiag
141  for (int i = 1; i <= S.num_row(); i++) {
142  double s2 = D(i,i);
143  if ( s2 > 0 ) {
144  sigmas(i) = sqrt ( s2 );
145  } else {
146  std::cerr << "In RandMultiGauss distribution: \n" <<
147  " Matrix S is not positive definite. Eigenvalues are:\n";
148  for (int ixx = 1; ixx <= S.num_row(); ixx++) {
149  std::cerr << " " << D(ixx,ixx) << std::endl;
150  }
151  std::cerr << "---Exiting to System\n";
152  exit(1);
153  }
154  }
155 } // prepareUsigmas
156 
157 // -----------
158 // deviates()
159 // -----------
160 
161 HepVector RandMultiGauss::deviates ( const HepMatrix & U,
162  const HepVector & sigmas,
163  HepRandomEngine * engine,
164  bool& available,
165  double& next)
166 {
167  // Returns vector of gaussian randoms based on sigmas, rotated by U,
168  // with means of 0.
169 
170  int n = sigmas.num_row();
171  HepVector v(n); // The vector to be returned
172 
173  double r,v1,v2,fac;
174 
175  int i = 1;
176  if (available) {
177  v(1) = next;
178  i = 2;
179  available = false;
180  }
181 
182  while ( i <= n ) {
183  do {
184  v1 = 2.0 * engine->flat() - 1.0;
185  v2 = 2.0 * engine->flat() - 1.0;
186  r = v1*v1 + v2*v2;
187  } while ( r > 1.0 );
188  fac = sqrt(-2.0*log(r)/r);
189  v(i++) = v1*fac;
190  if ( i <= n ) {
191  v(i++) = v2*fac;
192  } else {
193  next = v2*fac;
194  available = true;
195  }
196  }
197 
198  for ( i = 1; i <= n; i++ ) {
199  v(i) *= sigmas(i);
200  }
201 
202  return U*v;
203 
204 } // deviates()
205 
206 // ---------------
207 // fire signatures
208 // ---------------
209 
211  // Returns a pair of unit normals, using the S and mu set in constructor,
212  // utilizing the engine belonging to this instance of RandMultiGauss.
213 
214  return defaultMu + deviates ( defaultU, defaultSigmas,
215  localEngine, set, nextGaussian );
216 
217 } // fire();
218 
219 
221 
222  HepMatrix U;
223  HepVector sigmas;
224 
225  if (mu.num_row() == S.num_row()) {
226  prepareUsigmas ( S, U, sigmas );
227  return mu + deviates ( U, sigmas, localEngine, set, nextGaussian );
228  } else {
229  std::cerr << "In firing RandMultiGauss distribution with explicit mu and S: \n"
230  << " Dimension of mu (" << mu.num_row() <<
231  ") does not match dimension of S (" << S.num_row() << ")\n";
232  std::cerr << "---Exiting to System\n";
233  exit(1);
234  }
235  return mu; // This line cannot be reached. But without returning
236  // some HepVector here, KCC 3.3 complains.
237 
238 } // fire(mu, S);
239 
240 
241 // --------------------
242 // fireArray signatures
243 // --------------------
244 
245 void RandMultiGauss::fireArray( const int size, HepVector* array ) {
246 
247  int i;
248  for (i = 0; i < size; ++i) {
249  array[i] = defaultMu + deviates ( defaultU, defaultSigmas,
250  localEngine, set, nextGaussian );
251  }
252 
253 } // fireArray ( size, vect )
254 
255 
256 void RandMultiGauss::fireArray( const int size, HepVector* array,
257  const HepVector& mu, const HepSymMatrix& S ) {
258 
259  // For efficiency, we diagonalize S once and generate all the vectors based
260  // on that U and sigmas.
261 
262  HepMatrix U;
263  HepVector sigmas;
264  HepVector mu_ (mu);
265 
266  if (mu.num_row() == S.num_row()) {
267  prepareUsigmas ( S, U, sigmas );
268  } else {
269  std::cerr <<
270  "In fireArray for RandMultiGauss distribution with explicit mu and S: \n"
271  << " Dimension of mu (" << mu.num_row() <<
272  ") does not match dimension of S (" << S.num_row() << ")\n";
273  std::cerr << "---Exiting to System\n";
274  exit(1);
275  }
276 
277  int i;
278  for (i=0; i<size; ++i) {
279  array[i] = mu_ + deviates(U, sigmas, localEngine, set, nextGaussian);
280  }
281 
282 } // fireArray ( size, vect, mu, S )
283 
284 // ----------
285 // operator()
286 // ----------
287 
289  return fire();
290 }
291 
292 HepVector RandMultiGauss::operator()
293  ( const HepVector& mu, const HepSymMatrix& S ) {
294  return fire(mu,S);
295 }
296 
297 
298 } // namespace CLHEP
virtual int num_row() const
Definition: Vector.cc:117
int num_row() const
HepSymMatrix similarityT(const HepMatrix &hm1) const
Definition: SymMatrix.cc:816
#define exit(x)
Definition: excDblThrow.cc:17
HepMatrix diagonalize(HepSymMatrix *s)
RandMultiGauss(HepRandomEngine &anEngine, const HepVector &mu, const HepSymMatrix &S)
void fireArray(const int size, HepVector *array)