CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: testVectorDists.cc,v 1.2 2003/08/13 20:00:13 garren Exp $ 00002 // -*- C++ -*- 00003 // ---------------------------------------------------------------------- 00004 00005 // ---------------------------------------------------------------------- 00006 // 00007 // testVectorDists -- tests of the correctness of vector random distributions 00008 // 00009 // Currently tested: 00010 // RandMultiGauss 00011 // 00012 // ---------------------------------------------------------------------- 00013 00014 #include "CLHEP/RandomObjects/defs.h" 00015 #include "CLHEP/Random/Randomize.h" 00016 #include "CLHEP/RandomObjects/RandMultiGauss.h" 00017 #include "CLHEP/Matrix/SymMatrix.h" 00018 #include "CLHEP/Matrix/Matrix.h" 00019 #include "CLHEP/Matrix/Vector.h" 00020 #include <iostream> 00021 00022 using namespace std; 00023 using namespace CLHEP; 00024 00025 static const int MultiGaussBAD = 1 << 0; 00026 00027 00028 static HepMatrix modifiedOutput(const HepMatrix& D) { 00029 HepMatrix DD (D); 00030 int n = DD.num_row(); 00031 int m = DD.num_col(); 00032 int i; 00033 int j; 00034 for ( i = 1; i <= n; ++i ) { 00035 for ( j = 1; j <= m; ++j ) { 00036 if ( DD(i,j)*DD(i,j) < 1.0e-24 * DD(i,i) * DD(j,j) ) DD (i,j) = 0; 00037 } 00038 } 00039 return DD; 00040 } 00041 00042 00043 // -------------- 00044 // RandMultiGauss 00045 // -------------- 00046 00047 int testRandMultiGauss( ) { 00048 00049 cout << "\n--------------------------------------------\n"; 00050 cout << "Test of RandMultiGauss distribution \n\n"; 00051 00052 long seed; 00053 cout << "Please enter an integer seed: "; 00054 cin >> seed; 00055 00056 int nvectors; 00057 cout << "How many vectors should we generate: "; 00058 cin >> nvectors; 00059 double rootn = sqrt((double)nvectors); 00060 00061 int nMu; 00062 int nS; 00063 cout << "Enter the dimensions of mu, then S (normally the same): "; 00064 cin >> nMu >> nS; 00065 if ( nMu != nS ) { 00066 cout << "Usually mu and S will be of equal dimensions.\n"; 00067 cout << "You may be testing the behavior when that is not the case.\n"; 00068 cout << "Please verify by re-entering the correct dimensions: "; 00069 cin >> nMu >> nS; 00070 } 00071 int dim = (nMu >= nS) ? nMu : nS; 00072 00073 HepVector mu(nMu); 00074 HepSymMatrix S(nS); 00075 00076 cout << "Enter mu, one component at a time: \n"; 00077 int imu; 00078 double muElement; 00079 for (imu = 1; imu <= nMu; imu++) { 00080 cout << imu << ": "; 00081 cin >> muElement; 00082 mu(imu) = muElement; 00083 } 00084 00085 cout << "Enter S, one white-space-separated row at a time. \n"; 00086 cout << "The length needed for each row is given in {braces}.\n"; 00087 cout << 00088 "The diagonal elements of S will be the first numbers on each line:\n"; 00089 int row, col; 00090 double sij; 00091 for (row = 1; row <= nS; row++) { 00092 cout << row << " {" << nS - row + 1 << " numbers}: "; 00093 for (col = row; col <= nS; col++) { 00094 cin >> sij; 00095 S(row, col) = sij; 00096 } 00097 } 00098 00099 cout << "mu is: \n"; 00100 cout << mu; 00101 cout << endl; 00102 00103 cout << "S is: \n"; 00104 cout << S << endl; 00105 00106 HepSymMatrix tempS ( S ); // Since diagonalize does not take a const s 00107 // we have to copy S. 00108 HepMatrix U = diagonalize ( &tempS ); 00109 HepSymMatrix D = S.similarityT(U); 00110 cout << "D = Diagonalized S is " << modifiedOutput(D) << endl; 00111 bool pd = true; 00112 for ( int npd = 1; npd <= dim; npd++) { 00113 if ( D(npd,npd) < 0 ) { 00114 pd = false; 00115 } 00116 } 00117 if (!pd) { 00118 cout << "S is not positive definite.\n" << 00119 "The negative elements of D will have been raised to zero.\n" << 00120 "The second moment matrix at the end will not match S.\n"; 00121 } 00122 00123 cout << "\nInstantiating distribution utilizing TripleRand engine...\n"; 00124 TripleRand eng (seed); 00125 RandMultiGauss dist (eng, mu, S); 00126 00127 HepVector x(dim); 00128 00129 cout << "\n Sample fire(): \n"; 00130 00131 x = dist.fire(); 00132 cout << x; 00133 00134 int ntrials; 00135 cout << "Normal operation will try a single group of " << nvectors 00136 << " random vectors.\n" 00137 << "Enter 1 for a single trial with " << nvectors 00138 << " random vectors.\n" 00139 << "Alternatively some number of groups of " << nvectors 00140 << " vectors can be produced to accumulate deviation statistics.\n" 00141 << "Enter " << 5000/(dim*(dim+1))+1 00142 << " or some other number of trials to do this: "; 00143 cin >> ntrials; 00144 if (ntrials < 1) return 0; 00145 00146 cout << "\n Testing fire() ... \n"; 00147 00148 // I'm going to accumulate correlation matrix by equation (28.9) RPP 00149 // and compare to the specified matrix S. That is, E(x-<x>)(y-<y>) should 00150 // be Sxy. 00151 // 00152 // For off-diagonal terms, Sxy = <xy> - <x><y>. 00153 // 00154 // For diagonal terms, Sxx = <x^2> - <x>^2. 00155 00156 HepSymMatrix Sumxy(nS); 00157 HepSymMatrix Dprime(dim); 00158 HepSymMatrix VarD(dim); 00159 HepSymMatrix Delta(dim); 00160 00161 int ipr = nvectors / 10; if (ipr < 1) ipr = 1; 00162 int in1 = 0; 00163 int in2 = 0; 00164 int in3 = 0; 00165 int nentries = 0; 00166 float binno; 00167 int nbin; 00168 int bins[30]; 00169 int ix, iy; 00170 // double root2 = sqrt(2.0); 00171 double sumDelta = 0.0; 00172 double sumDelta2 = 0.0; 00173 int nunder = 0; 00174 int nover = 0; 00175 double worstDeviation=0; 00176 00177 int k; 00178 for(k=0; k<30; ++k) { 00179 bins[k] = 0; 00180 } 00181 for(k=0; k<ntrials; ++k ) { 00182 HepVector sumx(dim,0); 00183 HepMatrix sumxy(dim,dim,0); 00184 for ( int ifire = 1; ifire <= nvectors; ifire++) { 00185 x = dist.fire(); 00186 for (ix = 1; ix <= dim; ix++) { 00187 sumx(ix) += x(ix); 00188 for (iy = 1; iy <= dim; iy++) { 00189 sumxy(ix,iy) += x(ix)*x(iy); 00190 } 00191 } 00192 if ( (ifire % ipr) == 0 && k == 0 ) { 00193 cout << ifire << ", "; 00194 } 00195 } 00196 if( k == 0 ) 00197 cout << "Statistics for the first (or only) trial of " << nvectors 00198 << " vectors:\n\n"; 00199 00200 sumx = sumx / nvectors; 00201 if( k == 0 ) cout << "\nAverage (should match mu): \n" << sumx << endl; 00202 for (ix = 1; ix <= dim; ix++) { 00203 for (iy = 1; iy <= dim; iy++) { 00204 sumxy(ix,iy) = sumxy(ix,iy) / nvectors - sumx(ix)*sumx(iy); 00205 } 00206 } 00207 if (pd) { 00208 if( k == 0 ) cout << "Second Moments (should match S)\n" << sumxy << endl; 00209 } else { 00210 if( k == 0 ) cout << "Second Moments \n" << sumxy << endl; 00211 } 00212 00213 // Now transform sumxy with the same matrix that diagonalized S. Call the 00214 // result Dprime. There is a bit of fooling around here because sumxy is a 00215 // general matrix and similarityT() acts on a symmetric matrix. 00216 00217 Sumxy.assign(sumxy); 00218 Dprime = Sumxy.similarityT(U); 00219 00220 if( k == 0 ) cout << "\nDprime = Second moment matrix transformed by the same matrix that diagonalized S\n" << Dprime << endl; 00221 00222 for( ix=1; ix<=dim; ++ix ) { 00223 for( iy=ix; iy<=dim; ++iy ) { 00224 if( ix == iy ) { 00225 VarD(ix,iy) = 2.0*Dprime(ix,iy)*Dprime(ix,iy)/rootn; 00226 } else { 00227 VarD(ix,iy) = Dprime(ix,ix)*Dprime(iy,iy)/rootn; 00228 } 00229 } 00230 } 00231 if( k == 0 ) cout << "\nThe expected variance for Dprime elements is \n" 00232 << VarD << endl; 00233 00234 for( ix=1; ix<=dim; ++ix ) { 00235 for( iy=ix; iy<=dim; ++iy ) { 00236 Delta(ix,iy) = sqrt(rootn)*(D(ix,iy)-Dprime(ix,iy))/sqrt(VarD(ix,iy)); 00237 if (k==0) { 00238 if (abs(Delta(ix,iy)) > abs(worstDeviation)) { 00239 worstDeviation = Delta(ix,iy); 00240 } 00241 } 00242 } 00243 } 00244 00245 if( k == 0 ) { 00246 cout 00247 << "\nDifferences between each element and its normative value,\n" 00248 << "scaled by the expected deviation (sqrt(variance)) are: \n" 00249 << Delta << endl; 00250 } 00251 00252 if( k == 0 ) { 00253 cout << 00254 "About 5% of the above values should have absolute value more than 2.\n" 00255 << "Deviations of more than 4 sigma would probably indicate a problem.\n"; 00256 } 00257 00258 // Do a little counting 00259 00260 for( ix=1; ix<=dim; ++ix ) { 00261 for( iy=ix; iy<=dim; ++iy ) { 00262 if( Delta(ix,iy) >= -1.0 && Delta(ix,iy) <= 1.0 ) in1++; 00263 if( Delta(ix,iy) >= -2.0 && Delta(ix,iy) <= 2.0 ) in2++; 00264 if( Delta(ix,iy) >= -3.0 && Delta(ix,iy) <= 3.0 ) in3++; 00265 sumDelta += Delta(ix,iy); 00266 sumDelta2 += Delta(ix,iy)*Delta(ix,iy); 00267 binno = 5.0*(Delta(ix,iy)+3.0); 00268 if( binno < 0.0 ) ++nunder; 00269 if( binno > 30.0 ) ++nover; 00270 if( binno >= 0.0 && binno < 30.0 ) { 00271 nbin = (int)binno; 00272 ++nentries; 00273 ++bins[nbin]; 00274 } 00275 } 00276 } 00277 } 00278 00279 if (ntrials == 1) { 00280 cout << "The worst deviation of any element of D in this trial was " 00281 << worstDeviation << "\n"; 00282 if (abs(worstDeviation) > 4) { 00283 cout << "\nREJECT this distribution: \n" 00284 << "This value indicates a PROBLEM!!!!\n\n"; 00285 return MultiGaussBAD; 00286 } else { 00287 return 0; 00288 } 00289 } 00290 00291 float ndf = ntrials*dim*(dim+1)/2.0; 00292 cout << "\nOut of a total of " << ndf << " entries" << endl; 00293 cout << "There are " << in1 << " within 1 sigma or " 00294 << 100.0*(float)in1/ndf << "%" << endl; 00295 cout << "There are " << in2 << " within 2 sigma or " 00296 << 100.0*(float)in2/ndf << "%" << endl; 00297 cout << "There are " << in3 << " within 3 sigma or " 00298 << 100.0*(float)in3/ndf << "%" << endl; 00299 double aveDelta = sumDelta/(double)ndf; 00300 double rmsDelta = sumDelta2/(double)ndf - aveDelta*aveDelta; 00301 cout << "\nFor dim = " << dim << " Average(Delta) = " << aveDelta << " RMS(Delta) = " << rmsDelta << endl; 00302 cout << "\nPoor man's histogram of deviations in 30 bins from -3.0 to 3.0" << endl; 00303 cout << "This should be a standard unit Gaussian.\n" << endl; 00304 for(k=0; k<30; ++k ) { 00305 cout << setw(3) << k+1 << " " << setw(4) << bins[k] << endl; 00306 } 00307 cout << "\nTotal number of entries in this histogram is " << nentries << endl; 00308 cout << "\twith " << nunder << " underflow(s) and " << nover << " overflow(s)" << endl; 00309 00310 int status; 00311 00312 cout << "The mean squared delta should be between .85 and 1.15; it is " 00313 << rmsDelta << "\n"; 00314 00315 if( abs(rmsDelta-1.0) <= 0.15 ) { 00316 status = false; 00317 } else { 00318 cout << "REJECT this distribution based on improper spread of delta!\n"; 00319 status = MultiGaussBAD; 00320 } 00321 if (abs(worstDeviation)>4) { 00322 cout << "REJECT this distribution: Bad correlation in first trial!\n"; 00323 status = MultiGaussBAD; 00324 } 00325 00326 return status; 00327 00328 } // testRandMultiGauss() 00329 00330 int main() { 00331 00332 int mask = 0; 00333 00334 mask |= testRandMultiGauss(); 00335 00336 return mask; 00337 00338 } 00339 00340