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