Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // G4ConvergenceTester class implementation 27 // 28 // Author: Koi, Tatsumi (SLAC/SCCS) 29 // ------------------------------------------- 30 31 #include "G4ConvergenceTester.hh" 32 #include "G4AutoLock.hh" 33 #include <iomanip> 34 35 namespace 36 { 37 G4Mutex aMutex = G4MUTEX_INITIALIZER; 38 } 39 40 G4ConvergenceTester::G4ConvergenceTester(const 41 : name(theName) 42 { 43 nonzero_histories.clear(); 44 largest_scores.clear(); 45 largest_scores.push_back(0.0); 46 47 history_grid.resize(noBinOfHistory, 0); 48 mean_history.resize(noBinOfHistory, 0.0); 49 var_history.resize(noBinOfHistory, 0.0); 50 sd_history.resize(noBinOfHistory, 0.0); 51 r_history.resize(noBinOfHistory, 0.0); 52 vov_history.resize(noBinOfHistory, 0.0); 53 fom_history.resize(noBinOfHistory, 0.0); 54 shift_history.resize(noBinOfHistory, 0.0); 55 e_history.resize(noBinOfHistory, 0.0); 56 r2eff_history.resize(noBinOfHistory, 0.0); 57 r2int_history.resize(noBinOfHistory, 0.0); 58 59 timer = new G4Timer(); 60 timer->Start(); 61 cpu_time.clear(); 62 cpu_time.push_back(0.0); 63 } 64 65 G4ConvergenceTester::~G4ConvergenceTester() 66 { 67 delete timer; 68 } 69 70 void G4ConvergenceTester::AddScore(G4double x) 71 { 72 G4AutoLock l(&aMutex); 73 74 timer->Stop(); 75 cpu_time.push_back(timer->GetSystemElapsed() 76 77 if(x < 0.0) 78 { 79 std::ostringstream message; 80 message << "Expecting zero or positive num 81 << "but received a negative number 82 G4Exception("G4ConvergenceTester::AddScore 83 JustWarning, message); 84 } 85 86 if(x == 0.0) 87 { 88 } 89 else 90 { 91 nonzero_histories.insert(std::pair<G4int, 92 if(x > largest_scores.back()) 93 { 94 // Following search should become faster 95 for(auto it = largest_scores.begin(); it 96 { 97 if(x > *it) 98 { 99 largest_scores.insert(it, x); 100 break; 101 } 102 } 103 104 if(largest_scores.size() > 201) 105 { 106 largest_scores.pop_back(); 107 } 108 } 109 sum += x; 110 } 111 112 // Data has been added so statistics have no 113 statsAreUpdated = false; 114 ++n; 115 l.unlock(); 116 return; 117 } 118 119 void G4ConvergenceTester::calStat() 120 { 121 if (n == 0) return; 122 123 efficiency = G4double(nonzero_histories.size 124 125 mean = sum / n; 126 127 G4double sum_x2 = 0.0; 128 var = 0.0; 129 shift = 0.0; 130 vov = 0.0; 131 132 G4double xi; 133 for(const auto& nonzero_historie : nonzero_h 134 { 135 xi = nonzero_historie.second; 136 sum_x2 += xi * xi; 137 var += (xi - mean) * (xi - mean); 138 shift += (xi - mean) * (xi - mean) * (xi - 139 vov += (xi - mean) * (xi - mean) * (xi - m 140 } 141 142 var += (n - nonzero_histories.size()) * mean 143 shift += (n - nonzero_histories.size()) * me 144 vov += (n - nonzero_histories.size()) * mean 145 146 if(var != 0.0) 147 { 148 vov = vov / (var * var) - 1.0 / n; 149 150 var = var / (n - 1); 151 152 sd = std::sqrt(var); 153 154 r = sd / mean / std::sqrt(G4double(n)); 155 156 r2eff = (1 - efficiency) / (efficiency * n 157 r2int = sum_x2 / (sum * sum) - 1 / (effici 158 159 shift = shift / (2 * var * n); 160 161 fom = 1 / (r * r) / cpu_time.back(); 162 } 163 164 // Find Largest History 165 // G4double largest = 0.0; 166 largest = 0.0; 167 largest_score_happened = 0; 168 G4double spend_time_of_largest = 0.0; 169 for(const auto& nonzero_historie : nonzero_h 170 { 171 if(std::abs(nonzero_historie.second) > lar 172 { 173 largest = nonzero_histori 174 largest_score_happened = nonzero_histori 175 spend_time_of_largest = 176 cpu_time[nonzero_historie.first + 1] - 177 } 178 } 179 180 mean_1 = 0.0; 181 var_1 = 0.0; 182 shift_1 = 0.0; 183 vov_1 = 0.0; 184 sd_1 = 0.0; 185 r_1 = 0.0; 186 vov_1 = 0.0; 187 188 mean_1 = (sum + largest) / (n + 1); 189 190 for(const auto& nonzero_historie : nonzero_h 191 { 192 xi = nonzero_historie.second; 193 var_1 += (xi - mean_1) * (xi - mean_1); 194 shift_1 += (xi - mean_1) * (xi - mean_1) * 195 vov_1 += (xi - mean_1) * (xi - mean_1) * ( 196 } 197 xi = largest; 198 var_1 += (xi - mean_1) * (xi - mean_1); 199 shift_1 += (xi - mean_1) * (xi - mean_1) * ( 200 vov_1 += (xi - mean_1) * (xi - mean_1) * (xi 201 202 var_1 += (n - nonzero_histories.size()) * me 203 204 if(var_1 != 0.0) 205 { 206 shift_1 += (n - nonzero_histories.size()) 207 vov_1 += (n - nonzero_histories.size()) * 208 209 vov_1 = vov_1 / (var_1 * var_1) - 1.0 / (n 210 211 var_1 = var_1 / n; 212 213 sd_1 = std::sqrt(var_1); 214 215 r_1 = sd_1 / mean_1 / std::sqrt(G4double(n 216 217 shift_1 = shift_1 / (2 * var_1 * (n + 1)); 218 219 fom_1 = 1 / (r * r) / (cpu_time.back() + s 220 } 221 222 if(nonzero_histories.size() < 500) 223 { 224 calcSLOPE = false; 225 } 226 else 227 { 228 auto i = G4int(nonzero_histories.size()); 229 230 // 5% criterion 231 auto j = G4int(i * 0.05); 232 while(G4int(largest_scores.size()) > j) 233 { 234 largest_scores.pop_back(); 235 } 236 calc_slope_fit(largest_scores); 237 } 238 239 calc_grid_point_of_history(); 240 calc_stat_history(); 241 242 // statistics have been calculated and this 243 // to be called again until data has been ad 244 statsAreUpdated = true; 245 } 246 247 void G4ConvergenceTester::calc_grid_point_of_h 248 { 249 // histroy_grid [ 0,,,15 ] 250 // history_grid [0] 1/16 ,,, history_grid [1 251 // if number of event is x then history_grid 252 // 16 -> noBinOfHisotry 253 254 for(G4int i = 1; i <= noBinOfHistory; ++i) 255 { 256 history_grid[i - 1] = G4int(n / (G4double( 257 } 258 } 259 260 void G4ConvergenceTester::calc_stat_history() 261 { 262 if(history_grid[0] == 0) 263 { 264 showHistory = false; 265 return; 266 } 267 268 for(G4int i = 0; i < noBinOfHistory; ++i) 269 { 270 G4int ith = history_grid[i]; 271 272 G4int nonzero_till_ith = 0; 273 G4double xi; 274 G4double mean_till_ith = 0.0; 275 276 for(const auto& itr : nonzero_histories) 277 { 278 if(itr.first <= ith) 279 { 280 xi = itr.second; 281 mean_till_ith += xi; 282 ++nonzero_till_ith; 283 } 284 } 285 286 if(nonzero_till_ith == 0) 287 { 288 continue; 289 } 290 291 mean_till_ith = mean_till_ith / (ith + 1 292 mean_history[i] = mean_till_ith; 293 294 G4double sum_x2_till_ith = 0.0; 295 G4double var_till_ith = 0.0; 296 G4double vov_till_ith = 0.0; 297 G4double shift_till_ith = 0.0; 298 299 for(const auto& itr : nonzero_histories) 300 { 301 if(itr.first <= ith) 302 { 303 xi = itr.second; 304 sum_x2_till_ith += std::pow(xi, 2.0); 305 var_till_ith += std::pow(xi - mean_til 306 shift_till_ith += std::pow(xi - mean_t 307 vov_till_ith += std::pow(xi - mean_til 308 } 309 } 310 311 var_till_ith += 312 ((ith + 1) - nonzero_till_ith) * std::po 313 vov_till_ith += 314 ((ith + 1) - nonzero_till_ith) * std::po 315 316 G4double sum_till_ith = mean_till_ith * (i 317 318 if(!(std::fabs(var_till_ith) > 0.0)) 319 { 320 continue; 321 } 322 if(!(std::fabs(mean_till_ith) > 0.0)) 323 { 324 continue; 325 } 326 if(!(std::fabs(sum_till_ith) > 0.0)) 327 { 328 continue; 329 } 330 331 vov_till_ith = vov_till_ith / std::pow(var 332 vov_history[i] = vov_till_ith; 333 334 var_till_ith = var_till_ith / (ith + 1 - 335 var_history[i] = var_till_ith; 336 sd_history[i] = std::sqrt(var_till_ith); 337 r_history[i] = 338 std::sqrt(var_till_ith) / mean_till_ith 339 340 if(std::fabs(cpu_time[ith]) > 0.0 && std:: 341 { 342 fom_history[i] = 1.0 / std::pow(r_histor 343 } 344 else 345 { 346 fom_history[i] = 0.0; 347 } 348 349 shift_till_ith += 350 ((ith + 1) - nonzero_till_ith) * std::po 351 shift_till_ith = shift_till_ith / (2 * v 352 shift_history[i] = shift_till_ith; 353 354 e_history[i] = 1.0 * nonzero_till_ith / (i 355 if(std::fabs(e_history[i]) > 0.0) 356 { 357 r2eff_history[i] = (1 - e_history[i]) / 358 359 r2int_history[i] = (sum_x2_till_ith) / s 360 1 / (e_history[i] * ( 361 } 362 } 363 } 364 365 void G4ConvergenceTester::ShowResult(std::ostr 366 { 367 // if data has been added since the last com 368 // (not statsAreUpdated) call calStat to rec 369 if(!statsAreUpdated) 370 { 371 calStat(); 372 } 373 374 out << std::setprecision(6); 375 376 out << G4endl; 377 out << "G4ConvergenceTester Output Result of 378 out << std::setw(20) << "EFFICIENCY = " << s 379 << G4endl; 380 out << std::setw(20) << "MEAN = " << std::se 381 out << std::setw(20) << "VAR = " << std::set 382 out << std::setw(20) << "SD = " << std::setw 383 out << std::setw(20) << "R = " << std::setw( 384 out << std::setw(20) << "SHIFT = " << std::s 385 out << std::setw(20) << "VOV = " << std::set 386 out << std::setw(20) << "FOM = " << std::set 387 388 out << std::setw(20) << "THE LARGEST SCORE = 389 << " and it happened at " << largest_sco 390 << G4endl; 391 if(mean != 0) 392 { 393 out << std::setw(20) << "Affected Mean = " 394 << " and its ratio to original is " << 395 } 396 else 397 { 398 out << std::setw(20) << "Affected Mean = " 399 << G4endl; 400 } 401 if(var != 0) 402 { 403 out << std::setw(20) << "Affected VAR = " 404 << " and its ratio to original is " << 405 } 406 else 407 { 408 out << std::setw(20) << "Affected VAR = " 409 << G4endl; 410 } 411 if(r != 0) 412 { 413 out << std::setw(20) << "Affected R = " << 414 << " and its ratio to original is " << 415 } 416 else 417 { 418 out << std::setw(20) << "Affected R = " << 419 } 420 if(shift != 0) 421 { 422 out << std::setw(20) << "Affected SHIFT = 423 << " and its ratio to original is " << 424 } 425 else 426 { 427 out << std::setw(20) << "Affected SHIFT = 428 << G4endl; 429 } 430 if(fom != 0) 431 { 432 out << std::setw(20) << "Affected FOM = " 433 << " and its ratio to original is " << 434 } 435 else 436 { 437 out << std::setw(20) << "Affected FOM = " 438 << G4endl; 439 } 440 441 if(!showHistory) 442 { 443 out << "Number of events of this run is to 444 << G4endl; 445 return; 446 } 447 448 check_stat_history(out); 449 450 // check SLOPE and output result 451 if(calcSLOPE) 452 { 453 if(slope >= 3) 454 { 455 noPass++; 456 out << "SLOPE is large enough" << G4endl 457 } 458 else 459 { 460 out << "SLOPE is not large enough" << G4 461 } 462 } 463 else 464 { 465 out << "Number of non zero history too sma 466 } 467 468 out << "This result passes " << noPass << " 469 << " Convergence Test." << G4endl; 470 out << G4endl; 471 } 472 473 void G4ConvergenceTester::ShowHistory(std::ost 474 { 475 if(!showHistory) 476 { 477 out << "Number of events of this run is to 478 << G4endl; 479 return; 480 } 481 482 out << std::setprecision(6); 483 484 out << G4endl; 485 out << "G4ConvergenceTester Output History o 486 out << "i/" << noBinOfHistory << " till_ith 487 << "var" << std::setw(13) << "sd" << std 488 << "vov" << std::setw(13) << "fom" << st 489 << std::setw(13) << "e" << std::setw(13) 490 << "r2int" << G4endl; 491 for(G4int i = 1; i <= noBinOfHistory; i++) 492 { 493 out << std::setw(4) << i << " " << std::se 494 << std::setw(13) << mean_history[i - 1 495 << var_history[i - 1] << std::setw(13) 496 << std::setw(13) << r_history[i - 1] < 497 << vov_history[i - 1] << std::setw(13) 498 << std::setw(13) << shift_history[i - 499 << e_history[i - 1] << std::setw(13) < 500 << std::setw(13) << r2int_history[i - 501 } 502 } 503 504 void G4ConvergenceTester::check_stat_history(s 505 { 506 // 1 sigma rejection for null hypothesis 507 508 std::vector<G4double> first_ally; 509 std::vector<G4double> second_ally; 510 511 // use 2nd half of hisories 512 std::size_t N = mean_history.size() / 2; 513 std::size_t i; 514 515 G4double pearson_r; 516 G4double t; 517 518 first_ally.resize(N); 519 second_ally.resize(N); 520 521 G4double sum_of_var = 522 std::accumulate(var_history.begin(), var_h 523 if(sum_of_var == 0.0) 524 { 525 out << "Variances in all historical grids 526 out << "Terminating checking behavior of s 527 return; 528 } 529 530 // Mean 531 532 for(i = 0; i < N; ++i) 533 { 534 first_ally[i] = history_grid[N + i]; 535 second_ally[i] = mean_history[N + i]; 536 } 537 538 pearson_r = calc_Pearson_r((G4int)N, first_a 539 t = pearson_r * std::sqrt((N - 2) / 540 541 if(t < 0.429318) // Student t of (Degree of 542 { 543 out << "MEAN distribution is RANDOM" << G 544 noPass++; 545 } 546 else 547 { 548 out << "MEAN distribution is not RANDOM" < 549 } 550 551 // R 552 553 for(i = 0; i < N; ++i) 554 { 555 first_ally[i] = 1.0 / std::sqrt(G4double( 556 second_ally[i] = r_history[N + i]; 557 } 558 559 pearson_r = calc_Pearson_r(G4int(N), first_a 560 t = pearson_r * std::sqrt((N - 2) / 561 562 if(t > 1.090546) 563 { 564 out << "r follows 1/std::sqrt(N)" << G4end 565 noPass++; 566 } 567 else 568 { 569 out << "r does not follow 1/std::sqrt(N)" 570 } 571 572 if(is_monotonically_decrease(second_ally)) 573 { 574 out << "r is monotonically decrease " << G 575 } 576 else 577 { 578 out << "r is NOT monotonically decrease " 579 } 580 581 if(r_history.back() < 0.1) 582 { 583 out << "r is less than 0.1. r = " << r_his 584 noPass++; 585 } 586 else 587 { 588 out << "r is NOT less than 0.1. r = " << r 589 } 590 591 // VOV 592 for(i = 0; i < N; ++i) 593 { 594 first_ally[i] = 1.0 / history_grid[N + i] 595 second_ally[i] = vov_history[N + i]; 596 } 597 598 pearson_r = calc_Pearson_r(G4int(N), first_a 599 t = pearson_r * std::sqrt((N - 2) / 600 601 if(t > 1.090546) 602 { 603 out << "VOV follows 1/std::sqrt(N)" << G4e 604 noPass++; 605 } 606 else 607 { 608 out << "VOV does not follow 1/std::sqrt(N) 609 } 610 611 if(is_monotonically_decrease(second_ally)) 612 { 613 out << "VOV is monotonically decrease " << 614 } 615 else 616 { 617 out << "VOV is NOT monotonically decrease 618 } 619 620 // FOM 621 622 for(i = 0; i < N; ++i) 623 { 624 first_ally[i] = history_grid[N + i]; 625 second_ally[i] = fom_history[N + i]; 626 } 627 628 pearson_r = calc_Pearson_r(G4int(N), std::mo 629 t = pearson_r * std::sqrt((N - 2) / 630 631 if(t < 0.429318) 632 { 633 out << "FOM distribution is RANDOM" << G4e 634 noPass++; 635 } 636 else 637 { 638 out << "FOM distribution is not RANDOM" << 639 } 640 } 641 642 G4double G4ConvergenceTester::calc_Pearson_r(G 643 s 644 s 645 { 646 G4double first_mean = 0.0; 647 G4double second_mean = 0.0; 648 649 G4int i; 650 for(i = 0; i < N; i++) 651 { 652 first_mean += first_ally[i]; 653 second_mean += second_ally[i]; 654 } 655 first_mean = first_mean / N; 656 second_mean = second_mean / N; 657 658 G4double a = 0.0; 659 for(i = 0; i < N; ++i) 660 { 661 a += (first_ally[i] - first_mean) * (secon 662 } 663 664 G4double b1 = 0.0; 665 G4double b2 = 0.0; 666 for(i = 0; i < N; ++i) 667 { 668 b1 += (first_ally[i] - first_mean) * (firs 669 b2 += (second_ally[i] - second_mean) * (se 670 } 671 672 G4double rds = a / std::sqrt(b1 * b2); 673 674 return rds; 675 } 676 677 G4bool G4ConvergenceTester::is_monotonically_d 678 const std::vector<G4double>& ally) 679 { 680 for(auto it = ally.cbegin(); it != ally.cend 681 { 682 if(*it < *(it + 1)) 683 { 684 return FALSE; 685 } 686 } 687 688 ++noPass; 689 return TRUE; 690 } 691 692 void G4ConvergenceTester::calc_slope_fit(const 693 { 694 // create PDF bins 695 G4double max = largest_scores.front(); 696 auto last = G4int(largest_scores.size()); 697 G4double min = 0.0; 698 if(largest_scores.back() != 0) 699 { 700 min = largest_scores.back(); 701 } 702 else 703 { 704 min = largest_scores[last - 1]; 705 last = last - 1; 706 } 707 708 if(max * 0.99 < min) 709 { 710 // upper limit is assumed to have been rea 711 slope = 10.0; 712 return; 713 } 714 715 std::vector<G4double> pdf_grid; 716 717 pdf_grid.resize(noBinOfPDF + 1); // no grid 718 pdf_grid[0] = max; 719 pdf_grid[noBinOfPDF] = min; 720 G4double log10_max = std::log10(max); 721 G4double log10_min = std::log10(min); 722 G4double log10_delta = log10_max - log10_min 723 for(G4int i = 1; i < noBinOfPDF; ++i) 724 { 725 pdf_grid[i] = std::pow(10.0, log10_max - l 726 } 727 728 std::vector<G4double> pdf; 729 pdf.resize(noBinOfPDF); 730 731 for(G4int j = 0; j < last; ++j) 732 { 733 for(G4int i = 0; i < 11; ++i) 734 { 735 if(largest_scores[j] >= pdf_grid[i + 1]) 736 { 737 pdf[i] += 1.0 / (pdf_grid[i] - pdf_gri 738 break; 739 } 740 } 741 } 742 743 f_xi.resize(noBinOfPDF); 744 f_yi.resize(noBinOfPDF); 745 for(G4int i = 0; i < noBinOfPDF; ++i) 746 { 747 f_xi[i] = (pdf_grid[i] + pdf_grid[i + 1]) 748 f_yi[i] = pdf[i]; 749 } 750 751 // number of variables ( a and k ) 752 minimizer = new G4SimplexDownhill<G4Converge 753 // G4double minimum = minimizer->GetMinimum 754 std::vector<G4double> mp = minimizer->GetMin 755 G4double k = mp[1]; 756 757 // G4cout << "SLOPE " << 1/mp[1]+1 << G4endl 758 // G4cout << "SLOPE a " << mp[0] << G4endl; 759 // G4cout << "SLOPE k " << mp[1] << G4endl; 760 // G4cout << "SLOPE minimum " << minimizer- 761 762 slope = 1 / mp[1] + 1; 763 if(k < 1.0 / 9) // Please look Pareto distr 764 { 765 slope = 10; 766 } 767 if(slope > 10) 768 { 769 slope = 10; 770 } 771 } 772 773 G4double G4ConvergenceTester::slope_fitting_fu 774 { 775 G4double a = x[0]; 776 G4double k = x[1]; 777 778 if(a <= 0) 779 { 780 return 3.402823466e+38; // FLOAT_MAX 781 } 782 if(k == 0) 783 { 784 return 3.402823466e+38; // FLOAT_MAX 785 } 786 787 // f_xi and f_yi is filled at "calc_slope_fi 788 789 G4double y = 0.0; 790 for(G4int i = 0; i < G4int(f_yi.size()); ++i 791 { 792 // if ( 1/a * ( 1 + k * f_xi [ i ] / a ) < 793 if((1 + k * f_xi[i] / a) < 0) 794 { 795 y += 3.402823466e+38; // FLOAT_MAX 796 } 797 else 798 { 799 y += (f_yi[i] - 1 / a * std::pow(1 + k * 800 (f_yi[i] - 1 / a * std::pow(1 + k * 801 } 802 } 803 804 return y; 805 } 806