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

testVectorDists.cc
Go to the documentation of this file.
1 // $Id: testVectorDists.cc,v 1.2 2003/08/13 20:00:13 garren Exp $
2 // -*- C++ -*-
3 // ----------------------------------------------------------------------
4 
5 // ----------------------------------------------------------------------
6 //
7 // testVectorDists -- tests of the correctness of vector random distributions
8 //
9 // Currently tested:
10 // RandMultiGauss
11 //
12 // ----------------------------------------------------------------------
13 
15 #include "CLHEP/Random/Randomize.h"
17 #include "CLHEP/Matrix/SymMatrix.h"
18 #include "CLHEP/Matrix/Matrix.h"
19 #include "CLHEP/Matrix/Vector.h"
20 #include <iostream>
21 
22 using namespace std;
23 using namespace CLHEP;
24 
25 static const int MultiGaussBAD = 1 << 0;
26 
27 
28 static HepMatrix modifiedOutput(const HepMatrix& D) {
29  HepMatrix DD (D);
30  int n = DD.num_row();
31  int m = DD.num_col();
32  int i;
33  int j;
34  for ( i = 1; i <= n; ++i ) {
35  for ( j = 1; j <= m; ++j ) {
36  if ( DD(i,j)*DD(i,j) < 1.0e-24 * DD(i,i) * DD(j,j) ) DD (i,j) = 0;
37  }
38  }
39  return DD;
40 }
41 
42 
43 // --------------
44 // RandMultiGauss
45 // --------------
46 
48 
49  cout << "\n--------------------------------------------\n";
50  cout << "Test of RandMultiGauss distribution \n\n";
51 
52  long seed;
53  cout << "Please enter an integer seed: ";
54  cin >> seed;
55 
56  int nvectors;
57  cout << "How many vectors should we generate: ";
58  cin >> nvectors;
59  double rootn = sqrt((double)nvectors);
60 
61  int nMu;
62  int nS;
63  cout << "Enter the dimensions of mu, then S (normally the same): ";
64  cin >> nMu >> nS;
65  if ( nMu != nS ) {
66  cout << "Usually mu and S will be of equal dimensions.\n";
67  cout << "You may be testing the behavior when that is not the case.\n";
68  cout << "Please verify by re-entering the correct dimensions: ";
69  cin >> nMu >> nS;
70  }
71  int dim = (nMu >= nS) ? nMu : nS;
72 
73  HepVector mu(nMu);
74  HepSymMatrix S(nS);
75 
76  cout << "Enter mu, one component at a time: \n";
77  int imu;
78  double muElement;
79  for (imu = 1; imu <= nMu; imu++) {
80  cout << imu << ": ";
81  cin >> muElement;
82  mu(imu) = muElement;
83  }
84 
85  cout << "Enter S, one white-space-separated row at a time. \n";
86  cout << "The length needed for each row is given in {braces}.\n";
87  cout <<
88  "The diagonal elements of S will be the first numbers on each line:\n";
89  int row, col;
90  double sij;
91  for (row = 1; row <= nS; row++) {
92  cout << row << " {" << nS - row + 1 << " numbers}: ";
93  for (col = row; col <= nS; col++) {
94  cin >> sij;
95  S(row, col) = sij;
96  }
97  }
98 
99  cout << "mu is: \n";
100  cout << mu;
101  cout << endl;
102 
103  cout << "S is: \n";
104  cout << S << endl;
105 
106  HepSymMatrix tempS ( S ); // Since diagonalize does not take a const s
107  // we have to copy S.
108  HepMatrix U = diagonalize ( &tempS );
109  HepSymMatrix D = S.similarityT(U);
110  cout << "D = Diagonalized S is " << modifiedOutput(D) << endl;
111  bool pd = true;
112  for ( int npd = 1; npd <= dim; npd++) {
113  if ( D(npd,npd) < 0 ) {
114  pd = false;
115  }
116  }
117  if (!pd) {
118  cout << "S is not positive definite.\n" <<
119  "The negative elements of D will have been raised to zero.\n" <<
120  "The second moment matrix at the end will not match S.\n";
121  }
122 
123  cout << "\nInstantiating distribution utilizing TripleRand engine...\n";
124  TripleRand eng (seed);
125  RandMultiGauss dist (eng, mu, S);
126 
127  HepVector x(dim);
128 
129  cout << "\n Sample fire(): \n";
130 
131  x = dist.fire();
132  cout << x;
133 
134  int ntrials;
135  cout << "Normal operation will try a single group of " << nvectors
136  << " random vectors.\n"
137  << "Enter 1 for a single trial with " << nvectors
138  << " random vectors.\n"
139  << "Alternatively some number of groups of " << nvectors
140  << " vectors can be produced to accumulate deviation statistics.\n"
141  << "Enter " << 5000/(dim*(dim+1))+1
142  << " or some other number of trials to do this: ";
143  cin >> ntrials;
144  if (ntrials < 1) return 0;
145 
146  cout << "\n Testing fire() ... \n";
147 
148 // I'm going to accumulate correlation matrix by equation (28.9) RPP
149 // and compare to the specified matrix S. That is, E(x-<x>)(y-<y>) should
150 // be Sxy.
151 //
152 // For off-diagonal terms, Sxy = <xy> - <x><y>.
153 //
154 // For diagonal terms, Sxx = <x^2> - <x>^2.
155 
156  HepSymMatrix Sumxy(nS);
157  HepSymMatrix Dprime(dim);
158  HepSymMatrix VarD(dim);
159  HepSymMatrix Delta(dim);
160 
161  int ipr = nvectors / 10; if (ipr < 1) ipr = 1;
162  int in1 = 0;
163  int in2 = 0;
164  int in3 = 0;
165  int nentries = 0;
166  float binno;
167  int nbin;
168  int bins[30];
169  int ix, iy;
170 // double root2 = sqrt(2.0);
171  double sumDelta = 0.0;
172  double sumDelta2 = 0.0;
173  int nunder = 0;
174  int nover = 0;
175  double worstDeviation=0;
176 
177  int k;
178  for(k=0; k<30; ++k) {
179  bins[k] = 0;
180  }
181  for(k=0; k<ntrials; ++k ) {
182  HepVector sumx(dim,0);
183  HepMatrix sumxy(dim,dim,0);
184  for ( int ifire = 1; ifire <= nvectors; ifire++) {
185  x = dist.fire();
186  for (ix = 1; ix <= dim; ix++) {
187  sumx(ix) += x(ix);
188  for (iy = 1; iy <= dim; iy++) {
189  sumxy(ix,iy) += x(ix)*x(iy);
190  }
191  }
192  if ( (ifire % ipr) == 0 && k == 0 ) {
193  cout << ifire << ", ";
194  }
195  }
196  if( k == 0 )
197  cout << "Statistics for the first (or only) trial of " << nvectors
198  << " vectors:\n\n";
199 
200  sumx = sumx / nvectors;
201  if( k == 0 ) cout << "\nAverage (should match mu): \n" << sumx << endl;
202  for (ix = 1; ix <= dim; ix++) {
203  for (iy = 1; iy <= dim; iy++) {
204  sumxy(ix,iy) = sumxy(ix,iy) / nvectors - sumx(ix)*sumx(iy);
205  }
206  }
207  if (pd) {
208  if( k == 0 ) cout << "Second Moments (should match S)\n" << sumxy << endl;
209  } else {
210  if( k == 0 ) cout << "Second Moments \n" << sumxy << endl;
211  }
212 
213 // Now transform sumxy with the same matrix that diagonalized S. Call the
214 // result Dprime. There is a bit of fooling around here because sumxy is a
215 // general matrix and similarityT() acts on a symmetric matrix.
216 
217  Sumxy.assign(sumxy);
218  Dprime = Sumxy.similarityT(U);
219 
220  if( k == 0 ) cout << "\nDprime = Second moment matrix transformed by the same matrix that diagonalized S\n" << Dprime << endl;
221 
222  for( ix=1; ix<=dim; ++ix ) {
223  for( iy=ix; iy<=dim; ++iy ) {
224  if( ix == iy ) {
225  VarD(ix,iy) = 2.0*Dprime(ix,iy)*Dprime(ix,iy)/rootn;
226  } else {
227  VarD(ix,iy) = Dprime(ix,ix)*Dprime(iy,iy)/rootn;
228  }
229  }
230  }
231  if( k == 0 ) cout << "\nThe expected variance for Dprime elements is \n"
232  << VarD << endl;
233 
234  for( ix=1; ix<=dim; ++ix ) {
235  for( iy=ix; iy<=dim; ++iy ) {
236  Delta(ix,iy) = sqrt(rootn)*(D(ix,iy)-Dprime(ix,iy))/sqrt(VarD(ix,iy));
237  if (k==0) {
238  if (abs(Delta(ix,iy)) > abs(worstDeviation)) {
239  worstDeviation = Delta(ix,iy);
240  }
241  }
242  }
243  }
244 
245  if( k == 0 ) {
246  cout
247  << "\nDifferences between each element and its normative value,\n"
248  << "scaled by the expected deviation (sqrt(variance)) are: \n"
249  << Delta << endl;
250  }
251 
252  if( k == 0 ) {
253  cout <<
254  "About 5% of the above values should have absolute value more than 2.\n"
255  << "Deviations of more than 4 sigma would probably indicate a problem.\n";
256  }
257 
258 // Do a little counting
259 
260  for( ix=1; ix<=dim; ++ix ) {
261  for( iy=ix; iy<=dim; ++iy ) {
262  if( Delta(ix,iy) >= -1.0 && Delta(ix,iy) <= 1.0 ) in1++;
263  if( Delta(ix,iy) >= -2.0 && Delta(ix,iy) <= 2.0 ) in2++;
264  if( Delta(ix,iy) >= -3.0 && Delta(ix,iy) <= 3.0 ) in3++;
265  sumDelta += Delta(ix,iy);
266  sumDelta2 += Delta(ix,iy)*Delta(ix,iy);
267  binno = 5.0*(Delta(ix,iy)+3.0);
268  if( binno < 0.0 ) ++nunder;
269  if( binno > 30.0 ) ++nover;
270  if( binno >= 0.0 && binno < 30.0 ) {
271  nbin = (int)binno;
272  ++nentries;
273  ++bins[nbin];
274  }
275  }
276  }
277  }
278 
279  if (ntrials == 1) {
280  cout << "The worst deviation of any element of D in this trial was "
281  << worstDeviation << "\n";
282  if (abs(worstDeviation) > 4) {
283  cout << "\nREJECT this distribution: \n"
284  << "This value indicates a PROBLEM!!!!\n\n";
285  return MultiGaussBAD;
286  } else {
287  return 0;
288  }
289  }
290 
291  float ndf = ntrials*dim*(dim+1)/2.0;
292  cout << "\nOut of a total of " << ndf << " entries" << endl;
293  cout << "There are " << in1 << " within 1 sigma or "
294  << 100.0*(float)in1/ndf << "%" << endl;
295  cout << "There are " << in2 << " within 2 sigma or "
296  << 100.0*(float)in2/ndf << "%" << endl;
297  cout << "There are " << in3 << " within 3 sigma or "
298  << 100.0*(float)in3/ndf << "%" << endl;
299  double aveDelta = sumDelta/(double)ndf;
300  double rmsDelta = sumDelta2/(double)ndf - aveDelta*aveDelta;
301  cout << "\nFor dim = " << dim << " Average(Delta) = " << aveDelta << " RMS(Delta) = " << rmsDelta << endl;
302  cout << "\nPoor man's histogram of deviations in 30 bins from -3.0 to 3.0" << endl;
303  cout << "This should be a standard unit Gaussian.\n" << endl;
304  for(k=0; k<30; ++k ) {
305  cout << setw(3) << k+1 << " " << setw(4) << bins[k] << endl;
306  }
307  cout << "\nTotal number of entries in this histogram is " << nentries << endl;
308  cout << "\twith " << nunder << " underflow(s) and " << nover << " overflow(s)" << endl;
309 
310  int status;
311 
312  cout << "The mean squared delta should be between .85 and 1.15; it is "
313  << rmsDelta << "\n";
314 
315  if( abs(rmsDelta-1.0) <= 0.15 ) {
316  status = false;
317  } else {
318  cout << "REJECT this distribution based on improper spread of delta!\n";
319  status = MultiGaussBAD;
320  }
321  if (abs(worstDeviation)>4) {
322  cout << "REJECT this distribution: Bad correlation in first trial!\n";
323  status = MultiGaussBAD;
324  }
325 
326  return status;
327 
328 } // testRandMultiGauss()
329 
330 int main() {
331 
332  int mask = 0;
333 
334  mask |= testRandMultiGauss();
335 
336  return mask;
337 
338 }
339 
340 
#define double(obj)
Definition: excDblThrow.cc:32
int main()
STL namespace.
HepSymMatrix similarityT(const HepMatrix &hm1) const
Definition: SymMatrix.cc:816
Definition: excDblThrow.cc:17
int testRandMultiGauss()
void assign(const HepMatrix &hm2)
Definition: SymMatrix.cc:718
HepMatrix diagonalize(HepSymMatrix *s)