00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043 #include "G4ConvergenceTester.hh"
00044
00045 G4ConvergenceTester::G4ConvergenceTester()
00046 : n(0), sum(0.), mean(0.), var(0.), sd(0.), r(0.), efficiency(0.),
00047 r2eff(0.), r2int(0.), shift(0.), vov(0.), fom(0.), largest(0.),
00048 largest_score_happened(0), mean_1(0.), var_1(0.), sd_1(0.), r_1(0.),
00049 shift_1(0.), vov_1(0.), fom_1(0.), noBinOfHistory(16), slope(0.),
00050 noBinOfPDF(10), minimizer(0), noPass(0), noTotal(8), statsAreUpdated(true)
00051 {
00052 nonzero_histories.clear();
00053 largest_scores.clear();
00054 largest_scores.push_back( 0.0 );
00055
00056 history_grid.resize( noBinOfHistory , 0 );
00057 mean_history.resize( noBinOfHistory , 0.0 );
00058 var_history.resize( noBinOfHistory , 0.0 );
00059 sd_history.resize( noBinOfHistory , 0.0 );
00060 r_history.resize( noBinOfHistory , 0.0 );
00061 vov_history.resize( noBinOfHistory , 0.0 );
00062 fom_history.resize( noBinOfHistory , 0.0 );
00063 shift_history.resize( noBinOfHistory , 0.0 );
00064 e_history.resize( noBinOfHistory , 0.0 );
00065 r2eff_history.resize( noBinOfHistory , 0.0 );
00066 r2int_history.resize( noBinOfHistory , 0.0 );
00067
00068 timer = new G4Timer();
00069 timer->Start();
00070 cpu_time.clear();
00071 cpu_time.push_back( 0.0 );
00072 }
00073
00074
00075
00076 G4ConvergenceTester::~G4ConvergenceTester()
00077 {
00078 delete timer;
00079 }
00080
00081
00082
00083 void G4ConvergenceTester::AddScore( G4double x )
00084 {
00085
00086
00087
00088 timer->Stop();
00089 cpu_time.push_back( timer->GetSystemElapsed() + timer->GetUserElapsed() );
00090
00091 if ( x == 0.0 )
00092 {
00093 }
00094 else
00095 {
00096 nonzero_histories.insert( std::pair< G4int , G4double > ( n , x ) );
00097 if ( x > largest_scores.back() )
00098 {
00099
00100 std::vector< G4double >::iterator it;
00101 for ( it = largest_scores.begin() ; it != largest_scores.end() ; it++ )
00102 {
00103 if ( x > *it )
00104 {
00105 largest_scores.insert( it , x );
00106 break;
00107 }
00108 }
00109
00110 if ( largest_scores.size() > 201 )
00111 {
00112 largest_scores.pop_back();
00113 }
00114
00115 }
00116 sum += x;
00117 }
00118
00119
00120 statsAreUpdated = false;
00121 n++;
00122 return;
00123 }
00124
00125
00126
00127 void G4ConvergenceTester::calStat()
00128 {
00129
00130
00131 efficiency = double( nonzero_histories.size() ) / n;
00132
00133 mean = sum / n;
00134
00135 G4double sum_x2 = 0.0;
00136 var = 0.0;
00137 shift = 0.0;
00138 vov = 0.0;
00139
00140 G4double xi;
00141 std::map< G4int , G4double >::iterator it;
00142 for ( it = nonzero_histories.begin() ; it != nonzero_histories.end() ; it++ )
00143 {
00144 xi = it->second;
00145 sum_x2 += xi * xi;
00146 var += ( xi - mean ) * ( xi - mean );
00147 shift += ( xi - mean ) * ( xi - mean ) * ( xi - mean );
00148 vov += ( xi - mean ) * ( xi - mean ) * ( xi - mean ) * ( xi - mean );
00149 }
00150
00151 var += ( n - nonzero_histories.size() ) * mean * mean;
00152 shift += ( n - nonzero_histories.size() ) * mean * mean * mean * ( -1 );
00153 vov += ( n - nonzero_histories.size() ) * mean * mean * mean * mean;
00154
00155 vov = vov / ( var * var ) - 1.0 / n;
00156
00157 var = var/(n-1);
00158
00159 sd = std::sqrt ( var );
00160
00161 r = sd / mean / std::sqrt ( G4double(n) );
00162
00163 r2eff = ( 1 - efficiency ) / ( efficiency * n );
00164 r2int = sum_x2 / ( sum * sum ) - 1 / ( efficiency * n );
00165
00166 shift = shift / ( 2 * var * n );
00167
00168 fom = 1 / (r*r) / cpu_time.back();
00169
00170
00171
00172 largest = 0.0;
00173 largest_score_happened = 0;
00174 G4double spend_time_of_largest = 0.0;
00175 for ( it = nonzero_histories.begin() ; it != nonzero_histories.end() ; it++ )
00176 {
00177 if ( std::abs ( it->second ) > largest )
00178 {
00179 largest = it->second;
00180 largest_score_happened = it->first;
00181 spend_time_of_largest = cpu_time [ it->first+1 ] - cpu_time [ it->first ];
00182 }
00183 }
00184
00185 mean_1 = 0.0;
00186 var_1 = 0.0;
00187 shift_1 = 0.0;
00188 vov_1 = 0.0;
00189 sd_1 = 0.0;
00190 r_1 = 0.0;
00191 vov_1 = 0.0;
00192
00193
00194
00195 mean_1 = ( sum + largest ) / ( n + 1 );
00196
00197 for ( it = nonzero_histories.begin() ; it != nonzero_histories.end() ; it++ )
00198 {
00199 xi = it->second;
00200 var_1 += ( xi - mean_1 ) * ( xi - mean_1 );
00201 shift_1 += ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 );
00202 vov_1 += ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 );
00203 }
00204 xi = largest;
00205 var_1 += ( xi - mean_1 ) * ( xi - mean_1 );
00206 shift_1 += ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 );
00207 vov_1 += ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 ) * ( xi - mean_1 );
00208
00209 var_1 += ( n - nonzero_histories.size() ) * mean_1 * mean_1;
00210 shift_1 += ( n - nonzero_histories.size() ) * mean_1 * mean_1 * mean_1 * ( -1 );
00211 vov_1 += ( n - nonzero_histories.size() ) * mean_1 * mean_1 * mean_1 * mean_1;
00212
00213 vov_1 = vov_1 / ( var_1 * var_1 ) - 1.0 / ( n + 1 );
00214
00215 var_1 = var_1 / n ;
00216
00217 sd_1 = std::sqrt ( var_1 );
00218
00219 r_1 = sd_1 / mean_1 / std::sqrt ( G4double(n + 1) );
00220
00221 shift_1 = shift_1 / ( 2 * var_1 * ( n + 1 ) );
00222
00223 fom_1 = 1 / ( r * r ) / ( cpu_time.back() + spend_time_of_largest );
00224
00225 if ( nonzero_histories.size() < 500 )
00226 {
00227 G4cout << "Number of non zero history too small to calculate SLOPE" << G4endl;
00228 }
00229 else
00230 {
00231 G4int i = int ( nonzero_histories.size() );
00232
00233
00234 G4int j = int ( i * 0.05 );
00235 while ( int( largest_scores.size() ) > j )
00236 {
00237 largest_scores.pop_back();
00238 }
00239 calc_slope_fit( largest_scores );
00240 }
00241
00242 calc_grid_point_of_history();
00243 calc_stat_history();
00244
00245
00246
00247 statsAreUpdated = true;
00248 }
00249
00250
00251
00252 void G4ConvergenceTester::calc_grid_point_of_history()
00253 {
00254
00255
00256
00257
00258
00259
00260 G4int i;
00261 for ( i = 1 ; i <= noBinOfHistory ; i++ )
00262 {
00263 history_grid [ i-1 ] = int ( n / ( double( noBinOfHistory ) ) * i - 0.1 );
00264
00265 }
00266
00267 }
00268
00269
00270
00271 void G4ConvergenceTester::calc_stat_history()
00272 {
00273
00274
00275 G4int i;
00276 for ( i = 1 ; i <= noBinOfHistory ; i++ )
00277 {
00278
00279 G4int ith = history_grid [ i-1 ];
00280
00281 G4int nonzero_till_ith = 0;
00282 G4double xi;
00283 G4double mean_till_ith = 0.0;
00284 std::map< G4int , G4double >::iterator it;
00285
00286 for ( it = nonzero_histories.begin() ; it !=nonzero_histories.end() ; it++ )
00287 {
00288 if ( it->first <= ith )
00289 {
00290 xi = it->second;
00291 mean_till_ith += xi;
00292 nonzero_till_ith++;
00293 }
00294 }
00295
00296 mean_till_ith = mean_till_ith / ( ith+1 );
00297 mean_history [ i-1 ] = mean_till_ith;
00298
00299 G4double sum_x2_till_ith = 0.0;
00300 G4double var_till_ith = 0.0;
00301 G4double vov_till_ith = 0.0;
00302 G4double shift_till_ith = 0.0;
00303
00304 for ( it = nonzero_histories.begin() ; it !=nonzero_histories.end() ; it++ )
00305 {
00306 if ( it->first <= ith )
00307 {
00308 xi = it->second;
00309 sum_x2_till_ith += xi * xi;
00310 var_till_ith += ( xi - mean_till_ith ) * ( xi - mean_till_ith );
00311 shift_till_ith += ( xi - mean_till_ith ) * ( xi - mean_till_ith ) * ( xi - mean_till_ith );
00312 vov_till_ith += ( xi - mean_till_ith ) * ( xi - mean_till_ith ) * ( xi - mean_till_ith ) * ( xi - mean_till_ith );
00313 }
00314 }
00315
00316 var_till_ith += ( (ith+1) - nonzero_till_ith ) * mean_till_ith * mean_till_ith;
00317
00318 vov_till_ith += ( (ith+1) - nonzero_till_ith ) * mean_till_ith * mean_till_ith * mean_till_ith * mean_till_ith ;
00319 vov_till_ith = vov_till_ith / ( var_till_ith * var_till_ith ) - 1.0 / (ith+1);
00320 vov_history [ i-1 ] = vov_till_ith;
00321
00322 var_till_ith = var_till_ith / ( ith+1 - 1 );
00323 var_history [ i-1 ] = var_till_ith;
00324 sd_history [ i-1 ] = std::sqrt( var_till_ith );
00325 r_history [ i-1 ] = std::sqrt( var_till_ith ) / mean_till_ith / std::sqrt ( 1.0*(ith+1) );
00326
00327 fom_history [ i-1 ] = 1 / ( r_history [ i-1 ] * r_history [ i-1 ] ) / cpu_time [ ith ];
00328
00329 shift_till_ith += ( (ith+1) - nonzero_till_ith ) * mean_till_ith * mean_till_ith * mean_till_ith * ( -1 );
00330 shift_till_ith = shift_till_ith / ( 2 * var_till_ith * (ith+1) );
00331 shift_history [ i-1 ] = shift_till_ith;
00332
00333 e_history [ i-1 ] = 1.0*nonzero_till_ith / (ith+1);
00334 r2eff_history [ i-1 ] = ( 1 - e_history [ i-1 ] ) / ( e_history [ i-1 ] * (ith+1) );
00335
00336 G4double sum_till_ith = mean_till_ith * (ith+1);
00337 r2int_history [ i-1 ] = ( sum_x2_till_ith ) / ( sum_till_ith * sum_till_ith ) - 1 / ( e_history [ i-1 ] * (ith+1) );
00338
00339 }
00340
00341 }
00342
00343
00344
00345 void G4ConvergenceTester::ShowResult(std::ostream& out)
00346 {
00347
00348
00349 if(!statsAreUpdated) { calStat(); }
00350
00351 out << "EFFICIENCY = " << efficiency << G4endl;
00352 out << "MEAN = " << mean << G4endl;
00353 out << "VAR = " << var << G4endl;
00354 out << "SD = " << sd << G4endl;
00355 out << "R = "<< r << G4endl;
00356 out << "SHIFT = "<< shift << G4endl;
00357 out << "VOV = "<< vov << G4endl;
00358 out << "FOM = "<< fom << G4endl;
00359
00360 out << "THE LARGEST SCORE = " << largest << " and it happend at " << largest_score_happened << "th event" << G4endl;
00361 out << "Affected Mean = " << mean_1 << " and its ratio to orignal is " << mean_1/mean << G4endl;
00362 out << "Affected VAR = " << var_1 << " and its ratio to orignal is " << var_1/var << G4endl;
00363 out << "Affected R = " << r_1 << " and its ratio to orignal is " << r_1/r << G4endl;
00364 out << "Affected SHIFT = " << shift_1 << " and its ratio to orignal is " << shift_1/shift << G4endl;
00365 out << "Affected FOM = " << fom_1 << " and its ratio to orignal is " << fom_1/fom << G4endl;
00366
00367 check_stat_history(out);
00368
00369
00370 if ( slope >= 3 )
00371 {
00372 noPass++;
00373 out << "SLOPE is large enough" << G4endl;
00374 }
00375 else
00376 {
00377 out << "SLOPE is not large enough" << G4endl;
00378 }
00379
00380 out << "This result passes " << noPass << " / "<< noTotal << " Convergence Test." << G4endl;
00381 out << G4endl;
00382
00383 }
00384
00385 void G4ConvergenceTester::ShowHistory(std::ostream& out)
00386 {
00387 out << "i/" << noBinOfHistory << " till_ith mean var sd r vov fom shift e r2eff r2int" << G4endl;
00388 for ( G4int i = 1 ; i <= noBinOfHistory ; i++ )
00389 {
00390 out << i << " "
00391 << history_grid [ i-1 ] << " "
00392 << mean_history [ i-1 ] << " "
00393 << var_history [ i-1 ] << " "
00394 << sd_history [ i-1 ] << " "
00395 << r_history [ i-1 ] << " "
00396 << vov_history [ i-1 ] << " "
00397 << fom_history [ i-1 ] << " "
00398 << shift_history [ i-1 ] << " "
00399 << e_history [ i-1 ] << " "
00400 << r2eff_history [ i-1 ] << " "
00401 << r2int_history [ i-1 ] << " "
00402 << G4endl;
00403 }
00404 }
00405
00406 void G4ConvergenceTester::check_stat_history(std::ostream& out)
00407 {
00408
00409
00410
00411 std::vector<G4double> first_ally;
00412 std::vector<G4double> second_ally;
00413
00414
00415 G4int N = mean_history.size() / 2;
00416 G4int i;
00417
00418 G4double pearson_r;
00419 G4double t;
00420
00421 first_ally.resize( N );
00422 second_ally.resize( N );
00423
00424
00425
00426 for ( i = 0 ; i < N ; i++ )
00427 {
00428 first_ally [ i ] = history_grid [ N + i ];
00429 second_ally [ i ] = mean_history [ N + i ];
00430 }
00431
00432 pearson_r = calc_Pearson_r ( N , first_ally , second_ally );
00433 t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) );
00434
00435 if ( t < 0.429318 )
00436 {
00437 out << "MEAN distribution is RANDOM" << G4endl;
00438 noPass++;
00439 }
00440 else
00441 {
00442 out << "MEAN distribution is not RANDOM" << G4endl;
00443 }
00444
00445
00446
00447
00448 for ( i = 0 ; i < N ; i++ )
00449 {
00450 first_ally [ i ] = 1.0 / std::sqrt ( G4double(history_grid [ N + i ]) );
00451 second_ally [ i ] = r_history [ N + i ];
00452 }
00453
00454 pearson_r = calc_Pearson_r ( N , first_ally , second_ally );
00455 t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) );
00456
00457 if ( t > 1.090546 )
00458 {
00459 out << "r follows 1/std::sqrt(N)" << G4endl;
00460 noPass++;
00461 }
00462 else
00463 {
00464 out << "r does not follow 1/std::sqrt(N)" << G4endl;
00465 }
00466
00467 if ( is_monotonically_decrease( second_ally ) == true )
00468 {
00469 out << "r is monotonically decrease " << G4endl;
00470 }
00471 else
00472 {
00473 out << "r is NOT monotonically decrease " << G4endl;
00474 }
00475
00476 if ( r_history.back() < 0.1 )
00477 {
00478 out << "r is less than 0.1. r = " << r_history.back() << G4endl;
00479 noPass++;
00480 }
00481 else
00482 {
00483 out << "r is NOT less than 0.1. r = " << r_history.back() << G4endl;
00484 }
00485
00486
00487
00488 for ( i = 0 ; i < N ; i++ )
00489 {
00490 first_ally [ i ] = 1.0 / history_grid [ N + i ];
00491 second_ally [ i ] = vov_history [ N + i ];
00492 }
00493
00494 pearson_r = calc_Pearson_r ( N , first_ally , second_ally );
00495 t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) );
00496
00497 if ( t > 1.090546 )
00498 {
00499 out << "VOV follows 1/std::sqrt(N)" << G4endl;
00500 noPass++;
00501 }
00502 else
00503 {
00504 out << "VOV does not follow 1/std::sqrt(N)" << G4endl;
00505 }
00506
00507 if ( is_monotonically_decrease( second_ally ) == true )
00508 {
00509 out << "VOV is monotonically decrease " << G4endl;
00510 }
00511 else
00512 {
00513 out << "VOV is NOT monotonically decrease " << G4endl;
00514 }
00515
00516
00517
00518 for ( i = 0 ; i < N ; i++ )
00519 {
00520 first_ally [ i ] = history_grid [ N + i ];
00521 second_ally [ i ] = fom_history [ N + i ];
00522 }
00523
00524 pearson_r = calc_Pearson_r ( N , first_ally , second_ally );
00525 t = pearson_r * std::sqrt ( ( N - 2 ) / ( 1 - pearson_r * pearson_r ) );
00526
00527 if ( t < 0.429318 )
00528 {
00529 out << "FOM distribution is RANDOM" << G4endl;
00530 noPass++;
00531 }
00532 else
00533 {
00534 out << "FOM distribution is not RANDOM" << G4endl;
00535 }
00536
00537 }
00538
00539
00540
00541 G4double G4ConvergenceTester::calc_Pearson_r ( G4int N , std::vector<G4double> first_ally , std::vector<G4double> second_ally )
00542 {
00543 G4double first_mean = 0.0;
00544 G4double second_mean = 0.0;
00545
00546 G4int i;
00547 for ( i = 0 ; i < N ; i++ )
00548 {
00549 first_mean += first_ally [ i ];
00550 second_mean += second_ally [ i ];
00551 }
00552 first_mean = first_mean / N;
00553 second_mean = second_mean / N;
00554
00555 G4double a = 0.0;
00556 for ( i = 0 ; i < N ; i++ )
00557 {
00558 a += ( first_ally [ i ] - first_mean ) * ( second_ally [ i ] - second_mean );
00559 }
00560
00561 G4double b1 = 0.0;
00562 G4double b2 = 0.0;
00563 for ( i = 0 ; i < N ; i++ )
00564 {
00565 b1 += ( first_ally [ i ] - first_mean ) * ( first_ally [ i ] - first_mean );
00566 b2 += ( second_ally [ i ] - second_mean ) * ( second_ally [ i ] - second_mean );
00567 }
00568
00569 G4double rds = a / std::sqrt ( b1 * b2 );
00570
00571 return rds;
00572 }
00573
00574
00575
00576 G4bool G4ConvergenceTester::is_monotonically_decrease ( std::vector<G4double> ally )
00577 {
00578
00579 std::vector<G4double>::iterator it;
00580 for ( it = ally.begin() ; it != ally.end() - 1 ; it++ )
00581 {
00582 if ( *it < *(it+1) ) return FALSE;
00583 }
00584
00585 noPass++;
00586 return TRUE;
00587 }
00588
00589
00590
00591
00592 void G4ConvergenceTester::calc_slope_fit ( std::vector<G4double> )
00593 {
00594
00595
00596 G4double max = largest_scores.front();
00597 G4int last = int ( largest_scores.size() );
00598 G4double min = 0.0;
00599 if ( largest_scores.back() != 0 )
00600 {
00601 min = largest_scores.back();
00602 }
00603 else
00604 {
00605 min = largest_scores[ last-1 ];
00606 last = last - 1;
00607 }
00608
00609
00610
00611
00612 if ( max*0.99 < min )
00613 {
00614
00615 slope = 10.0;
00616 return;
00617 }
00618
00619 std::vector < G4double > pdf_grid;
00620
00621 pdf_grid.resize( noBinOfPDF+1 );
00622 pdf_grid[ 0 ] = max;
00623 pdf_grid[ noBinOfPDF ] = min;
00624 G4double log10_max = std::log10( max );
00625 G4double log10_min = std::log10( min );
00626 G4double log10_delta = log10_max - log10_min;
00627 for ( G4int i = 1 ; i < noBinOfPDF ; i++ )
00628 {
00629 pdf_grid[i] = std::pow ( 10.0 , log10_max - log10_delta/10.0*(i) );
00630
00631 }
00632
00633 std::vector < G4double > pdf;
00634 pdf.resize( noBinOfPDF );
00635
00636 for ( G4int j=0 ; j < last ; j ++ )
00637 {
00638 for ( G4int i = 0 ; i < 11 ; i++ )
00639 {
00640 if ( largest_scores[j] >= pdf_grid[i+1] )
00641 {
00642 pdf[i] += 1.0 / ( pdf_grid[i] - pdf_grid[i+1] ) / n;
00643
00644 break;
00645 }
00646 }
00647 }
00648
00649 f_xi.resize( noBinOfPDF );
00650 f_yi.resize( noBinOfPDF );
00651 for ( G4int i = 0 ; i < noBinOfPDF ; i++ )
00652 {
00653
00654 f_xi[i] = (pdf_grid[i]+pdf_grid[i+1])/2;
00655 f_yi[i] = pdf[i];
00656 }
00657
00658
00659 minimizer = new G4SimplexDownhill<G4ConvergenceTester> ( this , 2 );
00660
00661 std::vector<G4double> mp = minimizer->GetMinimumPoint();
00662 G4double k = mp[1];
00663
00664
00665
00666
00667
00668
00669 slope = 1/mp[1]+1;
00670 if ( k < 1.0/9 )
00671 {
00672 slope = 10;
00673 }
00674 if ( slope > 10 )
00675 {
00676 slope = 10;
00677 }
00678 }
00679
00680
00681
00682 G4double G4ConvergenceTester::slope_fitting_function ( std::vector< G4double > x )
00683 {
00684
00685 G4double a = x[0];
00686 G4double k = x[1];
00687
00688 if ( a <= 0 )
00689 {
00690 return 3.402823466e+38;
00691 }
00692 if ( k == 0 )
00693 {
00694 return 3.402823466e+38;
00695 }
00696
00697
00698
00699 G4double y = 0.0;
00700 G4int i;
00701 for ( i = 0 ; i < int ( f_yi.size() ) ; i++ )
00702 {
00703
00704 if ( ( 1 + k * f_xi [ i ] / a ) < 0 )
00705 {
00706 y +=3.402823466e+38;
00707 }
00708 else
00709 {
00710 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 ) );
00711 }
00712 }
00713
00714
00715 return y;
00716 }