[1350] | 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 | // Convergence Tests for Monte Carlo results. |
---|
| 29 | // |
---|
| 30 | // Reference |
---|
| 31 | // MCNP(TM) -A General Monte Carlo N-Particle Transport Code |
---|
| 32 | // Version 4B |
---|
| 33 | // Judith F. Briesmeister, Editor |
---|
| 34 | // LA-12625-M, Issued: March 1997, UC 705 and UC 700 |
---|
| 35 | // CHAPTER 2. GEOMETRY, DATA, PHYSICS, AND MATHEMATICS |
---|
| 36 | // VI. ESTIMATION OF THE MONTE CARLO PRECISION |
---|
| 37 | // |
---|
| 38 | // Positives numbers are assumed for inputs |
---|
| 39 | // |
---|
| 40 | // Koi, Tatsumi (SLAC/SCCS) |
---|
| 41 | // |
---|
| 42 | |
---|
| 43 | #include "G4ConvergenceTester.hh" |
---|
| 44 | |
---|
| 45 | G4ConvergenceTester::G4ConvergenceTester() |
---|
| 46 | : n(0), sum(0.), mean(0.), var(0.), sd(0.), r(0.), efficiency(0.), |
---|
| 47 | r2eff(0.), r2int(0.), shift(0.), vov(0.), fom(0.), largest(0.), |
---|
| 48 | largest_score_happened(0), mean_1(0.), var_1(0.), sd_1(0.), r_1(0.), |
---|
| 49 | shift_1(0.), vov_1(0.), fom_1(0.), noBinOfHistory(16), slope(0.), |
---|
| 50 | noBinOfPDF(10), minimizer(0), noPass(0), noTotal(8) |
---|
| 51 | { |
---|
| 52 | nonzero_histories.clear(); |
---|
| 53 | largest_scores.clear(); |
---|
| 54 | largest_scores.push_back( 0.0 ); |
---|
| 55 | |
---|
| 56 | history_grid.resize( noBinOfHistory , 0 ); |
---|
| 57 | mean_history.resize( noBinOfHistory , 0.0 ); |
---|
| 58 | var_history.resize( noBinOfHistory , 0.0 ); |
---|
| 59 | sd_history.resize( noBinOfHistory , 0.0 ); |
---|
| 60 | r_history.resize( noBinOfHistory , 0.0 ); |
---|
| 61 | vov_history.resize( noBinOfHistory , 0.0 ); |
---|
| 62 | fom_history.resize( noBinOfHistory , 0.0 ); |
---|
| 63 | shift_history.resize( noBinOfHistory , 0.0 ); |
---|
| 64 | e_history.resize( noBinOfHistory , 0.0 ); |
---|
| 65 | r2eff_history.resize( noBinOfHistory , 0.0 ); |
---|
| 66 | r2int_history.resize( noBinOfHistory , 0.0 ); |
---|
| 67 | |
---|
| 68 | timer = new G4Timer(); |
---|
| 69 | timer->Start(); |
---|
| 70 | cpu_time.clear(); |
---|
| 71 | cpu_time.push_back( 0.0 ); |
---|
| 72 | } |
---|
| 73 | |
---|
| 74 | |
---|
| 75 | |
---|
| 76 | G4ConvergenceTester::~G4ConvergenceTester() |
---|
| 77 | { |
---|
| 78 | delete timer; |
---|
| 79 | } |
---|
| 80 | |
---|
| 81 | |
---|
| 82 | |
---|
| 83 | void G4ConvergenceTester::AddScore( G4double x ) |
---|
| 84 | { |
---|
| 85 | |
---|
| 86 | //G4cout << x << G4endl; |
---|
| 87 | |
---|
| 88 | timer->Stop(); |
---|
| 89 | cpu_time.push_back( timer->GetSystemElapsed() + timer->GetUserElapsed() ); |
---|
| 90 | |
---|
| 91 | if ( x == 0.0 ) |
---|
| 92 | { |
---|
| 93 | } |
---|
| 94 | else |
---|
| 95 | { |
---|
| 96 | nonzero_histories.insert( std::pair< G4int , G4double > ( n , x ) ); |
---|
| 97 | if ( x > largest_scores.back() ) |
---|
| 98 | { |
---|
| 99 | // Following serch should become faster if begin from bottom. |
---|
| 100 | std::vector< G4double >::iterator it; |
---|
| 101 | for ( it = largest_scores.begin() ; it != largest_scores.end() ; it++ ) |
---|
| 102 | { |
---|
| 103 | if ( x > *it ) |
---|
| 104 | { |
---|
| 105 | largest_scores.insert( it , x ); |
---|
| 106 | break; |
---|
| 107 | } |
---|
| 108 | } |
---|
| 109 | |
---|
| 110 | if ( largest_scores.size() > 201 ) |
---|
| 111 | { |
---|
| 112 | largest_scores.pop_back(); |
---|
| 113 | } |
---|
| 114 | //G4cout << largest_scores.size() << " " << largest_scores.front() << " " << largest_scores.back() << G4endl; |
---|
| 115 | } |
---|
| 116 | sum += x; |
---|
| 117 | } |
---|
| 118 | |
---|
| 119 | n++; |
---|
| 120 | return; |
---|
| 121 | } |
---|
| 122 | |
---|
| 123 | |
---|
| 124 | |
---|
| 125 | void G4ConvergenceTester::calStat() |
---|
| 126 | { |
---|
| 127 | |
---|
| 128 | |
---|
| 129 | efficiency = double( nonzero_histories.size() ) / n; |
---|
| 130 | |
---|
| 131 | mean = sum / n; |
---|
| 132 | |
---|
| 133 | G4double sum_x2 = 0.0; |
---|
| 134 | var = 0.0; |
---|
| 135 | shift = 0.0; |
---|
| 136 | vov = 0.0; |
---|
| 137 | |
---|
| 138 | G4double xi; |
---|
| 139 | std::map< G4int , G4double >::iterator it; |
---|
| 140 | for ( it = nonzero_histories.begin() ; it != nonzero_histories.end() ; it++ ) |
---|
| 141 | { |
---|
| 142 | xi = it->second; |
---|
| 143 | sum_x2 += xi * xi; |
---|
| 144 | var += ( xi - mean ) * ( xi - mean ); |
---|
| 145 | shift += ( xi - mean ) * ( xi - mean ) * ( xi - mean ); |
---|
| 146 | vov += ( xi - mean ) * ( xi - mean ) * ( xi - mean ) * ( xi - mean ); |
---|
| 147 | } |
---|
| 148 | |
---|
| 149 | var += ( n - nonzero_histories.size() ) * mean * mean; |
---|
| 150 | shift += ( n - nonzero_histories.size() ) * mean * mean * mean * ( -1 ); |
---|
| 151 | vov += ( n - nonzero_histories.size() ) * mean * mean * mean * mean; |
---|
| 152 | |
---|
| 153 | vov = vov / ( var * var ) - 1.0 / n; |
---|
| 154 | |
---|
| 155 | var = var/(n-1); |
---|
| 156 | |
---|
| 157 | sd = std::sqrt ( var ); |
---|
| 158 | |
---|
| 159 | r = sd / mean / std::sqrt ( G4double(n) ); |
---|
| 160 | |
---|
| 161 | r2eff = ( 1 - efficiency ) / ( efficiency * n ); |
---|
| 162 | r2int = sum_x2 / ( sum * sum ) - 1 / ( efficiency * n ); |
---|
| 163 | |
---|
| 164 | shift = shift / ( 2 * var * n ); |
---|
| 165 | |
---|
| 166 | fom = 1 / (r*r) / cpu_time.back(); |
---|
| 167 | |
---|
| 168 | // Find Largest History |
---|
| 169 | //G4double largest = 0.0; |
---|
| 170 | largest = 0.0; |
---|
| 171 | largest_score_happened = 0; |
---|
| 172 | G4double spend_time_of_largest = 0.0; |
---|
| 173 | for ( it = nonzero_histories.begin() ; it != nonzero_histories.end() ; it++ ) |
---|
| 174 | { |
---|
| 175 | if ( std::abs ( it->second ) > largest ) |
---|
| 176 | { |
---|
| 177 | largest = it->second; |
---|
| 178 | largest_score_happened = it->first; |
---|
| 179 | spend_time_of_largest = cpu_time [ it->first+1 ] - cpu_time [ it->first ]; |
---|
| 180 | } |
---|
| 181 | } |
---|
| 182 | |
---|
| 183 | mean_1 = 0.0; |
---|
| 184 | var_1 = 0.0; |
---|
| 185 | shift_1 = 0.0; |
---|
| 186 | vov_1 = 0.0; |
---|
| 187 | sd_1 = 0.0; |
---|
| 188 | r_1 = 0.0; |
---|
| 189 | vov_1 = 0.0; |
---|
| 190 | |
---|
| 191 | // G4cout << "The largest history = " << largest << G4endl; |
---|
| 192 | |
---|
| 193 | mean_1 = ( sum + largest ) / ( n + 1 ); |
---|
| 194 | |
---|
| 195 | sum_x2 = 0.0; |
---|
| 196 | for ( it = nonzero_histories.begin() ; it != nonzero_histories.end() ; it++ ) |
---|
| 197 | { |
---|
| 198 | xi = it->second; |
---|
| 199 | sum_x2 += xi * xi; |
---|
| 200 | var_1 += ( xi - mean_1 ) * ( xi - mean_1 ); |
---|
| 201 | shift_1 += ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 ); |
---|
| 202 | vov_1 += ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 ); |
---|
| 203 | } |
---|
| 204 | xi = largest; |
---|
| 205 | sum_x2 += xi * xi; |
---|
| 206 | var_1 += ( xi - mean_1 ) * ( xi - mean_1 ); |
---|
| 207 | shift_1 += ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 ); |
---|
| 208 | vov_1 += ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 ); |
---|
| 209 | |
---|
| 210 | var_1 += ( n - nonzero_histories.size() ) * mean_1 * mean_1; |
---|
| 211 | shift_1 += ( n - nonzero_histories.size() ) * mean_1 * mean_1 * mean_1 * ( -1 ); |
---|
| 212 | vov_1 += ( n - nonzero_histories.size() ) * mean_1 * mean_1 * mean_1 * mean_1; |
---|
| 213 | |
---|
| 214 | vov_1 = vov_1 / ( var_1 * var_1 ) - 1.0 / ( n + 1 ); |
---|
| 215 | |
---|
| 216 | var_1 = var_1 / n ; |
---|
| 217 | |
---|
| 218 | sd_1 = std::sqrt ( var_1 ); |
---|
| 219 | |
---|
| 220 | r_1 = sd_1 / mean_1 / std::sqrt ( G4double(n + 1) ); |
---|
| 221 | |
---|
| 222 | shift_1 = shift_1 / ( 2 * var_1 * ( n + 1 ) ); |
---|
| 223 | |
---|
| 224 | fom_1 = 1 / ( r * r ) / ( cpu_time.back() + spend_time_of_largest ); |
---|
| 225 | |
---|
| 226 | if ( nonzero_histories.size() < 500 ) |
---|
| 227 | { |
---|
| 228 | G4cout << "Number of non zero history too small to calcuarte SLOPE" << G4endl; |
---|
| 229 | } |
---|
| 230 | else |
---|
| 231 | { |
---|
| 232 | G4int i = int ( nonzero_histories.size() ); |
---|
| 233 | |
---|
| 234 | // 5% criterion |
---|
| 235 | G4int j = int ( i * 0.05 ); |
---|
| 236 | while ( int( largest_scores.size() ) > j ) |
---|
| 237 | { |
---|
| 238 | largest_scores.pop_back(); |
---|
| 239 | } |
---|
| 240 | calc_slope_fit( largest_scores ); |
---|
| 241 | } |
---|
| 242 | |
---|
| 243 | calc_grid_point_of_history(); |
---|
| 244 | calc_stat_history(); |
---|
| 245 | } |
---|
| 246 | |
---|
| 247 | |
---|
| 248 | |
---|
| 249 | void G4ConvergenceTester::calc_grid_point_of_history() |
---|
| 250 | { |
---|
| 251 | |
---|
| 252 | // histroy_grid [ 0,,,15 ] |
---|
| 253 | // history_grid [0] 1/16 ,,, history_grid [15] 16/16 |
---|
| 254 | // if number of event is x then history_grid [15] become x-1. |
---|
| 255 | // 16 -> noBinOfHisotry |
---|
| 256 | |
---|
| 257 | G4int i; |
---|
| 258 | for ( i = 1 ; i <= noBinOfHistory ; i++ ) |
---|
| 259 | { |
---|
| 260 | history_grid [ i-1 ] = int ( n / ( double( noBinOfHistory ) ) * i - 0.1 ); |
---|
| 261 | //G4cout << "history_grid " << i-1 << " " << history_grid [ i-1 ] << G4endl; |
---|
| 262 | } |
---|
| 263 | |
---|
| 264 | } |
---|
| 265 | |
---|
| 266 | |
---|
| 267 | |
---|
| 268 | void G4ConvergenceTester::calc_stat_history() |
---|
| 269 | { |
---|
| 270 | // G4cout << "i/16 till_ith mean var sd r vov fom shift e r2eff r2int" << G4endl; |
---|
| 271 | |
---|
| 272 | G4int i; |
---|
| 273 | for ( i = 1 ; i <= noBinOfHistory ; i++ ) |
---|
| 274 | { |
---|
| 275 | |
---|
| 276 | G4int ith = history_grid [ i-1 ]; |
---|
| 277 | |
---|
| 278 | G4int nonzero_till_ith = 0; |
---|
| 279 | G4double xi; |
---|
| 280 | G4double mean_till_ith = 0.0; |
---|
| 281 | std::map< G4int , G4double >::iterator it; |
---|
| 282 | |
---|
| 283 | for ( it = nonzero_histories.begin() ; it !=nonzero_histories.end() ; it++ ) |
---|
| 284 | { |
---|
| 285 | if ( it->first <= ith ) |
---|
| 286 | { |
---|
| 287 | xi = it->second; |
---|
| 288 | mean_till_ith += xi; |
---|
| 289 | nonzero_till_ith++; |
---|
| 290 | } |
---|
| 291 | } |
---|
| 292 | |
---|
| 293 | mean_till_ith = mean_till_ith / ( ith+1 ); |
---|
| 294 | mean_history [ i-1 ] = mean_till_ith; |
---|
| 295 | |
---|
| 296 | G4double sum_x2_till_ith = 0.0; |
---|
| 297 | G4double var_till_ith = 0.0; |
---|
| 298 | G4double vov_till_ith = 0.0; |
---|
| 299 | G4double shift_till_ith = 0.0; |
---|
| 300 | |
---|
| 301 | for ( it = nonzero_histories.begin() ; it !=nonzero_histories.end() ; it++ ) |
---|
| 302 | { |
---|
| 303 | if ( it->first <= ith ) |
---|
| 304 | { |
---|
| 305 | xi = it->second; |
---|
| 306 | sum_x2_till_ith += xi * xi; |
---|
| 307 | var_till_ith += ( xi - mean_till_ith ) * ( xi - mean_till_ith ); |
---|
| 308 | shift_till_ith += ( xi - mean_till_ith ) * ( xi - mean_till_ith ) * ( xi - mean_till_ith ); |
---|
| 309 | vov_till_ith += ( xi - mean_till_ith ) * ( xi - mean_till_ith ) * ( xi - mean_till_ith ) * ( xi - mean_till_ith ); |
---|
| 310 | } |
---|
| 311 | } |
---|
| 312 | |
---|
| 313 | var_till_ith += ( (ith+1) - nonzero_till_ith ) * mean_till_ith * mean_till_ith; |
---|
| 314 | |
---|
| 315 | vov_till_ith += ( (ith+1) - nonzero_till_ith ) * mean_till_ith * mean_till_ith * mean_till_ith * mean_till_ith ; |
---|
| 316 | vov_till_ith = vov_till_ith / ( var_till_ith * var_till_ith ) - 1.0 / (ith+1); |
---|
| 317 | vov_history [ i-1 ] = vov_till_ith; |
---|
| 318 | |
---|
| 319 | var_till_ith = var_till_ith / ( ith+1 - 1 ); |
---|
| 320 | var_history [ i-1 ] = var_till_ith; |
---|
| 321 | sd_history [ i-1 ] = std::sqrt( var_till_ith ); |
---|
| 322 | r_history [ i-1 ] = std::sqrt( var_till_ith ) / mean_till_ith / std::sqrt ( 1.0*(ith+1) ); |
---|
| 323 | |
---|
| 324 | fom_history [ i-1 ] = 1 / ( r_history [ i-1 ] * r_history [ i-1 ] ) / cpu_time [ ith ]; |
---|
| 325 | |
---|
| 326 | shift_till_ith += ( (ith+1) - nonzero_till_ith ) * mean_till_ith * mean_till_ith * mean_till_ith * ( -1 ); |
---|
| 327 | shift_till_ith = shift_till_ith / ( 2 * var_till_ith * (ith+1) ); |
---|
| 328 | shift_history [ i-1 ] = shift_till_ith; |
---|
| 329 | |
---|
| 330 | e_history [ i-1 ] = 1.0*nonzero_till_ith / (ith+1); |
---|
| 331 | r2eff_history [ i-1 ] = ( 1 - e_history [ i-1 ] ) / ( e_history [ i-1 ] * (ith+1) ); |
---|
| 332 | |
---|
| 333 | G4double sum_till_ith = mean_till_ith * (ith+1); |
---|
| 334 | r2int_history [ i-1 ] = ( sum_x2_till_ith ) / ( sum_till_ith * sum_till_ith ) - 1 / ( e_history [ i-1 ] * (ith+1) ); |
---|
| 335 | |
---|
| 336 | } |
---|
| 337 | |
---|
| 338 | } |
---|
| 339 | |
---|
| 340 | |
---|
| 341 | |
---|
| 342 | void G4ConvergenceTester::ShowResult() |
---|
| 343 | { |
---|
| 344 | calStat(); |
---|
| 345 | |
---|
| 346 | G4cout << "EFFICIENCY = " << efficiency << G4endl; |
---|
| 347 | G4cout << "MEAN = " << mean << G4endl; |
---|
| 348 | G4cout << "VAR = " << var << G4endl; |
---|
| 349 | G4cout << "SD = " << sd << G4endl; |
---|
| 350 | G4cout << "R = "<< r << G4endl; |
---|
| 351 | G4cout << "SHIFT = "<< shift << G4endl; |
---|
| 352 | G4cout << "VOV = "<< vov << G4endl; |
---|
| 353 | G4cout << "FOM = "<< fom << G4endl; |
---|
| 354 | |
---|
| 355 | G4cout << "THE LARGEST SCORE = " << largest << " and it happend at " << largest_score_happened << "th event" << G4endl; |
---|
| 356 | G4cout << "Affected Mean = " << mean_1 << " and its ratio to orignal is " << mean_1/mean << G4endl; |
---|
| 357 | G4cout << "Affected VAR = " << var_1 << " and its ratio to orignal is " << var_1/var << G4endl; |
---|
| 358 | G4cout << "Affected R = " << r_1 << " and its ratio to orignal is " << r_1/r << G4endl; |
---|
| 359 | G4cout << "Affected SHIFT = " << shift_1 << " and its ratio to orignal is " << shift_1/shift << G4endl; |
---|
| 360 | G4cout << "Affected FOM = " << fom_1 << " and its ratio to orignal is " << fom_1/fom << G4endl; |
---|
| 361 | |
---|
| 362 | check_stat_history(); |
---|
| 363 | |
---|
| 364 | // check SLOPE and output result |
---|
| 365 | if ( slope >= 3 ) |
---|
| 366 | { |
---|
| 367 | noPass++; |
---|
| 368 | G4cout << "SLOPE is large enough" << G4endl; |
---|
| 369 | } |
---|
| 370 | else |
---|
| 371 | { |
---|
| 372 | G4cout << "SLOPE is not large enough" << G4endl; |
---|
| 373 | } |
---|
| 374 | |
---|
| 375 | G4cout << "This result passes " << noPass << " / "<< noTotal << " Convergence Test." << G4endl; |
---|
| 376 | G4cout << G4endl; |
---|
| 377 | |
---|
| 378 | } |
---|
| 379 | |
---|
| 380 | void G4ConvergenceTester::ShowHistory() |
---|
| 381 | { |
---|
| 382 | G4cout << "i/" << noBinOfHistory << " till_ith mean var sd r vov fom shift e r2eff r2int" << G4endl; |
---|
| 383 | for ( G4int i = 1 ; i <= noBinOfHistory ; i++ ) |
---|
| 384 | { |
---|
| 385 | G4cout << i << " " |
---|
| 386 | << history_grid [ i-1 ] << " " |
---|
| 387 | << mean_history [ i-1 ] << " " |
---|
| 388 | << var_history [ i-1 ] << " " |
---|
| 389 | << sd_history [ i-1 ] << " " |
---|
| 390 | << r_history [ i-1 ] << " " |
---|
| 391 | << vov_history [ i-1 ] << " " |
---|
| 392 | << fom_history [ i-1 ] << " " |
---|
| 393 | << shift_history [ i-1 ] << " " |
---|
| 394 | << e_history [ i-1 ] << " " |
---|
| 395 | << r2eff_history [ i-1 ] << " " |
---|
| 396 | << r2int_history [ i-1 ] << " " |
---|
| 397 | << G4endl; |
---|
| 398 | } |
---|
| 399 | } |
---|
| 400 | |
---|
| 401 | void G4ConvergenceTester::check_stat_history() |
---|
| 402 | { |
---|
| 403 | |
---|
| 404 | // 1 sigma rejection for null hypothesis |
---|
| 405 | |
---|
| 406 | std::vector<G4double> first_ally; |
---|
| 407 | std::vector<G4double> second_ally; |
---|
| 408 | |
---|
| 409 | // use 2nd half of hisories |
---|
| 410 | G4int N = mean_history.size() / 2; |
---|
| 411 | G4int i; |
---|
| 412 | |
---|
| 413 | G4double pearson_r; |
---|
| 414 | G4double t; |
---|
| 415 | |
---|
| 416 | first_ally.resize( N ); |
---|
| 417 | second_ally.resize( N ); |
---|
| 418 | |
---|
| 419 | // Mean |
---|
| 420 | |
---|
| 421 | for ( i = 0 ; i < N ; i++ ) |
---|
| 422 | { |
---|
| 423 | first_ally [ i ] = history_grid [ N + i ]; |
---|
| 424 | second_ally [ i ] = mean_history [ N + i ]; |
---|
| 425 | } |
---|
| 426 | |
---|
| 427 | pearson_r = calc_Pearson_r ( N , first_ally , second_ally ); |
---|
| 428 | t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) ); |
---|
| 429 | |
---|
| 430 | if ( t < 0.429318 ) // Student t of (Degree of freedom = N-2 ) |
---|
| 431 | { |
---|
| 432 | G4cout << "MEAN distribution is RANDOM" << G4endl; |
---|
| 433 | noPass++; |
---|
| 434 | } |
---|
| 435 | else |
---|
| 436 | { |
---|
| 437 | G4cout << "MEAN distribution is not RANDOM" << G4endl; |
---|
| 438 | } |
---|
| 439 | |
---|
| 440 | |
---|
| 441 | // R |
---|
| 442 | |
---|
| 443 | for ( i = 0 ; i < N ; i++ ) |
---|
| 444 | { |
---|
| 445 | first_ally [ i ] = 1.0 / std::sqrt ( G4double(history_grid [ N + i ]) ); |
---|
| 446 | second_ally [ i ] = r_history [ N + i ]; |
---|
| 447 | } |
---|
| 448 | |
---|
| 449 | pearson_r = calc_Pearson_r ( N , first_ally , second_ally ); |
---|
| 450 | t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) ); |
---|
| 451 | |
---|
| 452 | if ( t > 1.090546 ) |
---|
| 453 | { |
---|
| 454 | G4cout << "r follows 1/std::sqrt(N)" << G4endl; |
---|
| 455 | noPass++; |
---|
| 456 | } |
---|
| 457 | else |
---|
| 458 | { |
---|
| 459 | G4cout << "r does not follow 1/std::sqrt(N)" << G4endl; |
---|
| 460 | } |
---|
| 461 | |
---|
| 462 | if ( is_monotonically_decrease( second_ally ) == true ) |
---|
| 463 | { |
---|
| 464 | G4cout << "r is monotonically decrease " << G4endl; |
---|
| 465 | } |
---|
| 466 | else |
---|
| 467 | { |
---|
| 468 | G4cout << "r is NOT monotonically decrease " << G4endl; |
---|
| 469 | } |
---|
| 470 | |
---|
| 471 | if ( r_history.back() < 0.1 ) |
---|
| 472 | { |
---|
| 473 | G4cout << "r is less than 0.1. r = " << r_history.back() << G4endl; |
---|
| 474 | noPass++; |
---|
| 475 | } |
---|
| 476 | else |
---|
| 477 | { |
---|
| 478 | G4cout << "r is NOT less than 0.1. r = " << r_history.back() << G4endl; |
---|
| 479 | } |
---|
| 480 | |
---|
| 481 | |
---|
| 482 | // VOV |
---|
| 483 | for ( i = 0 ; i < N ; i++ ) |
---|
| 484 | { |
---|
| 485 | first_ally [ i ] = 1.0 / history_grid [ N + i ]; |
---|
| 486 | second_ally [ i ] = vov_history [ N + i ]; |
---|
| 487 | } |
---|
| 488 | |
---|
| 489 | pearson_r = calc_Pearson_r ( N , first_ally , second_ally ); |
---|
| 490 | t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) ); |
---|
| 491 | |
---|
| 492 | if ( t > 1.090546 ) |
---|
| 493 | { |
---|
| 494 | G4cout << "VOV follows 1/std::sqrt(N)" << G4endl; |
---|
| 495 | noPass++; |
---|
| 496 | } |
---|
| 497 | else |
---|
| 498 | { |
---|
| 499 | G4cout << "VOV does not follow 1/std::sqrt(N)" << G4endl; |
---|
| 500 | } |
---|
| 501 | |
---|
| 502 | if ( is_monotonically_decrease( second_ally ) == true ) |
---|
| 503 | { |
---|
| 504 | G4cout << "VOV is monotonically decrease " << G4endl; |
---|
| 505 | } |
---|
| 506 | else |
---|
| 507 | { |
---|
| 508 | G4cout << "VOV is NOT monotonically decrease " << G4endl; |
---|
| 509 | } |
---|
| 510 | |
---|
| 511 | // FOM |
---|
| 512 | |
---|
| 513 | for ( i = 0 ; i < N ; i++ ) |
---|
| 514 | { |
---|
| 515 | first_ally [ i ] = history_grid [ N + i ]; |
---|
| 516 | second_ally [ i ] = fom_history [ N + i ]; |
---|
| 517 | } |
---|
| 518 | |
---|
| 519 | pearson_r = calc_Pearson_r ( N , first_ally , second_ally ); |
---|
| 520 | t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) ); |
---|
| 521 | |
---|
| 522 | if ( t < 0.429318 ) |
---|
| 523 | { |
---|
| 524 | G4cout << "FOM distribution is RANDOM" << G4endl; |
---|
| 525 | noPass++; |
---|
| 526 | } |
---|
| 527 | else |
---|
| 528 | { |
---|
| 529 | G4cout << "FOM distribution is not RANDOM" << G4endl; |
---|
| 530 | } |
---|
| 531 | |
---|
| 532 | } |
---|
| 533 | |
---|
| 534 | |
---|
| 535 | |
---|
| 536 | G4double G4ConvergenceTester::calc_Pearson_r ( G4int N , std::vector<G4double> first_ally , std::vector<G4double> second_ally ) |
---|
| 537 | { |
---|
| 538 | G4double first_mean = 0.0; |
---|
| 539 | G4double second_mean = 0.0; |
---|
| 540 | |
---|
| 541 | G4int i; |
---|
| 542 | for ( i = 0 ; i < N ; i++ ) |
---|
| 543 | { |
---|
| 544 | first_mean += first_ally [ i ]; |
---|
| 545 | second_mean += second_ally [ i ]; |
---|
| 546 | } |
---|
| 547 | first_mean = first_mean / N; |
---|
| 548 | second_mean = second_mean / N; |
---|
| 549 | |
---|
| 550 | G4double a = 0.0; |
---|
| 551 | for ( i = 0 ; i < N ; i++ ) |
---|
| 552 | { |
---|
| 553 | a += ( first_ally [ i ] - first_mean ) * ( second_ally [ i ] - second_mean ); |
---|
| 554 | } |
---|
| 555 | |
---|
| 556 | G4double b1 = 0.0; |
---|
| 557 | G4double b2 = 0.0; |
---|
| 558 | for ( i = 0 ; i < N ; i++ ) |
---|
| 559 | { |
---|
| 560 | b1 += ( first_ally [ i ] - first_mean ) * ( first_ally [ i ] - first_mean ); |
---|
| 561 | b2 += ( second_ally [ i ] - second_mean ) * ( second_ally [ i ] - second_mean ); |
---|
| 562 | } |
---|
| 563 | |
---|
| 564 | G4double r = a / std::sqrt ( b1 * b2 ); |
---|
| 565 | |
---|
| 566 | return r; |
---|
| 567 | } |
---|
| 568 | |
---|
| 569 | |
---|
| 570 | |
---|
| 571 | G4bool G4ConvergenceTester::is_monotonically_decrease ( std::vector<G4double> ally ) |
---|
| 572 | { |
---|
| 573 | |
---|
| 574 | std::vector<G4double>::iterator it; |
---|
| 575 | for ( it = ally.begin() ; it != ally.end() - 1 ; it++ ) |
---|
| 576 | { |
---|
| 577 | if ( *it < *(it+1) ) return FALSE; |
---|
| 578 | } |
---|
| 579 | |
---|
| 580 | noPass++; |
---|
| 581 | return TRUE; |
---|
| 582 | } |
---|
| 583 | |
---|
| 584 | |
---|
| 585 | |
---|
| 586 | //void G4ConvergenceTester::calc_slope_fit ( std::vector<G4double> largest_socres ) |
---|
| 587 | void G4ConvergenceTester::calc_slope_fit ( std::vector<G4double> ) |
---|
| 588 | { |
---|
| 589 | |
---|
| 590 | // create PDF bins |
---|
| 591 | G4double max = largest_scores.front(); |
---|
| 592 | G4int last = int ( largest_scores.size() ); |
---|
| 593 | G4double min = 0.0; |
---|
| 594 | if ( largest_scores.back() != 0 ) |
---|
| 595 | { |
---|
| 596 | min = largest_scores.back(); |
---|
| 597 | } |
---|
| 598 | else |
---|
| 599 | { |
---|
| 600 | min = largest_scores[ last-1 ]; |
---|
| 601 | last = last - 1; |
---|
| 602 | } |
---|
| 603 | |
---|
| 604 | //G4cout << "largest " << max << G4endl; |
---|
| 605 | //G4cout << "last " << min << G4endl; |
---|
| 606 | |
---|
| 607 | if ( max*0.99 < min ) |
---|
| 608 | { |
---|
| 609 | // upper limit is assumed to have been reached |
---|
| 610 | slope = 10.0; |
---|
| 611 | return; |
---|
| 612 | } |
---|
| 613 | |
---|
| 614 | std::vector < G4double > pdf_grid; |
---|
| 615 | |
---|
| 616 | pdf_grid.resize( noBinOfPDF+1 ); // no grid = no bins + 1 |
---|
| 617 | pdf_grid[ 0 ] = max; |
---|
| 618 | pdf_grid[ noBinOfPDF ] = min; |
---|
| 619 | G4double log10_max = std::log10( max ); |
---|
| 620 | G4double log10_min = std::log10( min ); |
---|
| 621 | G4double log10_delta = log10_max - log10_min; |
---|
| 622 | for ( G4int i = 1 ; i < noBinOfPDF ; i++ ) |
---|
| 623 | { |
---|
| 624 | pdf_grid[i] = std::pow ( 10.0 , log10_max - log10_delta/10.0*(i) ); |
---|
| 625 | //G4cout << "pdf i " << i << " " << pdf_grid[i] << G4endl; |
---|
| 626 | } |
---|
| 627 | |
---|
| 628 | std::vector < G4double > pdf; |
---|
| 629 | pdf.resize( noBinOfPDF ); |
---|
| 630 | |
---|
| 631 | for ( G4int j=0 ; j < last ; j ++ ) |
---|
| 632 | { |
---|
| 633 | for ( G4int i = 0 ; i < 11 ; i++ ) |
---|
| 634 | { |
---|
| 635 | if ( largest_scores[j] >= pdf_grid[i+1] ) |
---|
| 636 | { |
---|
| 637 | pdf[i] += 1.0 / ( pdf_grid[i] - pdf_grid[i+1] ) / n; |
---|
| 638 | //G4cout << "pdf " << j << " " << i << " " << largest_scores[j] << " " << G4endl; |
---|
| 639 | break; |
---|
| 640 | } |
---|
| 641 | } |
---|
| 642 | } |
---|
| 643 | |
---|
| 644 | f_xi.resize( noBinOfPDF ); |
---|
| 645 | f_yi.resize( noBinOfPDF ); |
---|
| 646 | for ( G4int i = 0 ; i < noBinOfPDF ; i++ ) |
---|
| 647 | { |
---|
| 648 | //G4cout << "pdf i " << i << " " << (pdf_grid[i]+pdf_grid[i+1])/2 << " " << pdf[i] << G4endl; |
---|
| 649 | f_xi[i] = (pdf_grid[i]+pdf_grid[i+1])/2; |
---|
| 650 | f_yi[i] = pdf[i]; |
---|
| 651 | } |
---|
| 652 | |
---|
| 653 | // number of variables ( a and k ) |
---|
| 654 | minimizer = new G4SimplexDownhill<G4ConvergenceTester> ( this , 2 ); |
---|
| 655 | //G4double minimum = minimizer->GetMinimum(); |
---|
| 656 | std::vector<G4double> mp = minimizer->GetMinimumPoint(); |
---|
| 657 | G4double k = mp[1]; |
---|
| 658 | |
---|
| 659 | //G4cout << "SLOPE " << 1/mp[1]+1 << G4endl; |
---|
| 660 | //G4cout << "SLOPE a " << mp[0] << G4endl; |
---|
| 661 | //G4cout << "SLOPE k " << mp[1] << G4endl; |
---|
| 662 | //G4cout << "SLOPE minimum " << minimizer->GetMinimum() << G4endl; |
---|
| 663 | |
---|
| 664 | slope = 1/mp[1]+1; |
---|
| 665 | if ( k < 1.0/9 ) // Please look Pareto distribution with "sigma=a" and "k" |
---|
| 666 | { |
---|
| 667 | slope = 10; |
---|
| 668 | } |
---|
| 669 | if ( slope > 10 ) |
---|
| 670 | { |
---|
| 671 | slope = 10; |
---|
| 672 | } |
---|
| 673 | } |
---|
| 674 | |
---|
| 675 | |
---|
| 676 | |
---|
| 677 | G4double G4ConvergenceTester::slope_fitting_function ( std::vector< G4double > x ) |
---|
| 678 | { |
---|
| 679 | |
---|
| 680 | G4double a = x[0]; |
---|
| 681 | G4double k = x[1]; |
---|
| 682 | |
---|
| 683 | if ( a <= 0 ) |
---|
| 684 | { |
---|
| 685 | return 3.402823466e+38; // FLOAT_MAX |
---|
| 686 | } |
---|
| 687 | if ( k == 0 ) |
---|
| 688 | { |
---|
| 689 | return 3.402823466e+38; // FLOAT_MAX |
---|
| 690 | } |
---|
| 691 | |
---|
| 692 | // f_xi and f_yi is filled at "calc_slope_fit" |
---|
| 693 | |
---|
| 694 | G4double y = 0.0; |
---|
| 695 | G4int i; |
---|
| 696 | for ( i = 0 ; i < int ( f_yi.size() ) ; i++ ) |
---|
| 697 | { |
---|
| 698 | //if ( 1/a * ( 1 + k * f_xi [ i ] / a ) < 0 ) |
---|
| 699 | if ( ( 1 + k * f_xi [ i ] / a ) < 0 ) |
---|
| 700 | { |
---|
| 701 | y +=3.402823466e+38; // FLOAT_MAX |
---|
| 702 | } |
---|
| 703 | else |
---|
| 704 | { |
---|
| 705 | y += ( f_yi [ i ] - 1/a*std::pow ( 1 + k * f_xi [ i ] / a , - 1/k - 1 ) ) * ( f_yi [ i ] - 1/a*std::pow ( 1 + k * f_xi [ i ] / a , - 1/k - 1 ) ); |
---|
| 706 | } |
---|
| 707 | } |
---|
| 708 | // G4cout << "y = " << y << G4endl; |
---|
| 709 | |
---|
| 710 | return y; |
---|
| 711 | } |
---|