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