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 | } |
---|