| 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(auto it = nonzero_histories.cbegin(); it != nonzero_histories.cend();
>> 132 ++it)
134 { 133 {
135 xi = nonzero_historie.second; << 134 xi = it->second;
136 sum_x2 += xi * xi; 135 sum_x2 += xi * xi;
137 var += (xi - mean) * (xi - mean); 136 var += (xi - mean) * (xi - mean);
138 shift += (xi - mean) * (xi - mean) * (xi - 137 shift += (xi - mean) * (xi - mean) * (xi - mean);
139 vov += (xi - mean) * (xi - mean) * (xi - m 138 vov += (xi - mean) * (xi - mean) * (xi - mean) * (xi - mean);
140 } 139 }
141 140
142 var += (n - nonzero_histories.size()) * mean 141 var += (n - nonzero_histories.size()) * mean * mean;
143 shift += (n - nonzero_histories.size()) * me 142 shift += (n - nonzero_histories.size()) * mean * mean * mean * (-1);
144 vov += (n - nonzero_histories.size()) * mean 143 vov += (n - nonzero_histories.size()) * mean * mean * mean * mean;
145 144
146 if(var != 0.0) 145 if(var != 0.0)
147 { 146 {
148 vov = vov / (var * var) - 1.0 / n; 147 vov = vov / (var * var) - 1.0 / n;
149 148
150 var = var / (n - 1); 149 var = var / (n - 1);
151 150
152 sd = std::sqrt(var); 151 sd = std::sqrt(var);
153 152
154 r = sd / mean / std::sqrt(G4double(n)); 153 r = sd / mean / std::sqrt(G4double(n));
155 154
156 r2eff = (1 - efficiency) / (efficiency * n 155 r2eff = (1 - efficiency) / (efficiency * n);
157 r2int = sum_x2 / (sum * sum) - 1 / (effici 156 r2int = sum_x2 / (sum * sum) - 1 / (efficiency * n);
158 157
159 shift = shift / (2 * var * n); 158 shift = shift / (2 * var * n);
160 159
161 fom = 1 / (r * r) / cpu_time.back(); 160 fom = 1 / (r * r) / cpu_time.back();
162 } 161 }
163 162
164 // Find Largest History 163 // Find Largest History
165 // G4double largest = 0.0; 164 // G4double largest = 0.0;
166 largest = 0.0; 165 largest = 0.0;
167 largest_score_happened = 0; 166 largest_score_happened = 0;
168 G4double spend_time_of_largest = 0.0; 167 G4double spend_time_of_largest = 0.0;
169 for(const auto& nonzero_historie : nonzero_h << 168 for(auto it = nonzero_histories.cbegin(); it != nonzero_histories.cend();
>> 169 ++it)
170 { 170 {
171 if(std::abs(nonzero_historie.second) > lar << 171 if(std::abs(it->second) > largest)
172 { 172 {
173 largest = nonzero_histori << 173 largest = it->second;
174 largest_score_happened = nonzero_histori << 174 largest_score_happened = it->first;
175 spend_time_of_largest = << 175 spend_time_of_largest = cpu_time[it->first + 1] - cpu_time[it->first];
176 cpu_time[nonzero_historie.first + 1] - <<
177 } 176 }
178 } 177 }
179 178
180 mean_1 = 0.0; 179 mean_1 = 0.0;
181 var_1 = 0.0; 180 var_1 = 0.0;
182 shift_1 = 0.0; 181 shift_1 = 0.0;
183 vov_1 = 0.0; 182 vov_1 = 0.0;
184 sd_1 = 0.0; 183 sd_1 = 0.0;
185 r_1 = 0.0; 184 r_1 = 0.0;
186 vov_1 = 0.0; 185 vov_1 = 0.0;
187 186
188 mean_1 = (sum + largest) / (n + 1); 187 mean_1 = (sum + largest) / (n + 1);
189 188
190 for(const auto& nonzero_historie : nonzero_h << 189 for(auto it = nonzero_histories.cbegin(); it != nonzero_histories.cend();
>> 190 ++it)
191 { 191 {
192 xi = nonzero_historie.second; << 192 xi = it->second;
193 var_1 += (xi - mean_1) * (xi - mean_1); 193 var_1 += (xi - mean_1) * (xi - mean_1);
194 shift_1 += (xi - mean_1) * (xi - mean_1) * 194 shift_1 += (xi - mean_1) * (xi - mean_1) * (xi - mean_1);
195 vov_1 += (xi - mean_1) * (xi - mean_1) * ( 195 vov_1 += (xi - mean_1) * (xi - mean_1) * (xi - mean_1) * (xi - mean_1);
196 } 196 }
197 xi = largest; 197 xi = largest;
198 var_1 += (xi - mean_1) * (xi - mean_1); 198 var_1 += (xi - mean_1) * (xi - mean_1);
199 shift_1 += (xi - mean_1) * (xi - mean_1) * ( 199 shift_1 += (xi - mean_1) * (xi - mean_1) * (xi - mean_1);
200 vov_1 += (xi - mean_1) * (xi - mean_1) * (xi 200 vov_1 += (xi - mean_1) * (xi - mean_1) * (xi - mean_1) * (xi - mean_1);
201 201
202 var_1 += (n - nonzero_histories.size()) * me 202 var_1 += (n - nonzero_histories.size()) * mean_1 * mean_1;
203 203
204 if(var_1 != 0.0) 204 if(var_1 != 0.0)
205 { 205 {
206 shift_1 += (n - nonzero_histories.size()) 206 shift_1 += (n - nonzero_histories.size()) * mean_1 * mean_1 * mean_1 * (-1);
207 vov_1 += (n - nonzero_histories.size()) * 207 vov_1 += (n - nonzero_histories.size()) * mean_1 * mean_1 * mean_1 * mean_1;
208 208
209 vov_1 = vov_1 / (var_1 * var_1) - 1.0 / (n 209 vov_1 = vov_1 / (var_1 * var_1) - 1.0 / (n + 1);
210 210
211 var_1 = var_1 / n; 211 var_1 = var_1 / n;
212 212
213 sd_1 = std::sqrt(var_1); 213 sd_1 = std::sqrt(var_1);
214 214
215 r_1 = sd_1 / mean_1 / std::sqrt(G4double(n 215 r_1 = sd_1 / mean_1 / std::sqrt(G4double(n + 1));
216 216
217 shift_1 = shift_1 / (2 * var_1 * (n + 1)); 217 shift_1 = shift_1 / (2 * var_1 * (n + 1));
218 218
219 fom_1 = 1 / (r * r) / (cpu_time.back() + s 219 fom_1 = 1 / (r * r) / (cpu_time.back() + spend_time_of_largest);
220 } 220 }
221 221
222 if(nonzero_histories.size() < 500) 222 if(nonzero_histories.size() < 500)
223 { 223 {
224 calcSLOPE = false; 224 calcSLOPE = false;
225 } 225 }
226 else 226 else
227 { 227 {
228 auto i = G4int(nonzero_histories.size()); << 228 G4int i = G4int(nonzero_histories.size());
229 229
230 // 5% criterion 230 // 5% criterion
231 auto j = G4int(i * 0.05); << 231 G4int j = G4int(i * 0.05);
232 while(G4int(largest_scores.size()) > j) 232 while(G4int(largest_scores.size()) > j)
233 { 233 {
234 largest_scores.pop_back(); 234 largest_scores.pop_back();
235 } 235 }
236 calc_slope_fit(largest_scores); 236 calc_slope_fit(largest_scores);
237 } 237 }
238 238
239 calc_grid_point_of_history(); 239 calc_grid_point_of_history();
240 calc_stat_history(); 240 calc_stat_history();
241 241
242 // statistics have been calculated and this 242 // statistics have been calculated and this function does not need
243 // to be called again until data has been ad 243 // to be called again until data has been added
244 statsAreUpdated = true; 244 statsAreUpdated = true;
245 } 245 }
246 246
247 void G4ConvergenceTester::calc_grid_point_of_h 247 void G4ConvergenceTester::calc_grid_point_of_history()
248 { 248 {
249 // histroy_grid [ 0,,,15 ] 249 // histroy_grid [ 0,,,15 ]
250 // history_grid [0] 1/16 ,,, history_grid [1 250 // history_grid [0] 1/16 ,,, history_grid [15] 16/16
251 // if number of event is x then history_grid 251 // if number of event is x then history_grid [15] become x-1.
252 // 16 -> noBinOfHisotry 252 // 16 -> noBinOfHisotry
253 253
254 for(G4int i = 1; i <= noBinOfHistory; ++i) 254 for(G4int i = 1; i <= noBinOfHistory; ++i)
255 { 255 {
256 history_grid[i - 1] = G4int(n / (G4double( 256 history_grid[i - 1] = G4int(n / (G4double(noBinOfHistory)) * i - 0.1);
257 } 257 }
258 } 258 }
259 259
260 void G4ConvergenceTester::calc_stat_history() 260 void G4ConvergenceTester::calc_stat_history()
261 { 261 {
262 if(history_grid[0] == 0) 262 if(history_grid[0] == 0)
263 { 263 {
264 showHistory = false; 264 showHistory = false;
265 return; 265 return;
266 } 266 }
267 267
268 for(G4int i = 0; i < noBinOfHistory; ++i) 268 for(G4int i = 0; i < noBinOfHistory; ++i)
269 { 269 {
270 G4int ith = history_grid[i]; 270 G4int ith = history_grid[i];
271 271
272 G4int nonzero_till_ith = 0; 272 G4int nonzero_till_ith = 0;
273 G4double xi; 273 G4double xi;
274 G4double mean_till_ith = 0.0; 274 G4double mean_till_ith = 0.0;
275 275
276 for(const auto& itr : nonzero_histories) 276 for(const auto& itr : nonzero_histories)
277 { 277 {
278 if(itr.first <= ith) 278 if(itr.first <= ith)
279 { 279 {
280 xi = itr.second; 280 xi = itr.second;
281 mean_till_ith += xi; 281 mean_till_ith += xi;
282 ++nonzero_till_ith; 282 ++nonzero_till_ith;
283 } 283 }
284 } 284 }
285 285
286 if(nonzero_till_ith == 0) 286 if(nonzero_till_ith == 0)
287 { <<
288 continue; 287 continue;
289 } <<
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 { <<
320 continue; 317 continue;
321 } <<
322 if(!(std::fabs(mean_till_ith) > 0.0)) 318 if(!(std::fabs(mean_till_ith) > 0.0))
323 { <<
324 continue; 319 continue;
325 } <<
326 if(!(std::fabs(sum_till_ith) > 0.0)) 320 if(!(std::fabs(sum_till_ith) > 0.0))
327 { <<
328 continue; 321 continue;
329 } <<
330 322
331 vov_till_ith = vov_till_ith / std::pow(var 323 vov_till_ith = vov_till_ith / std::pow(var_till_ith, 2.0) - 1.0 / (ith + 1);
332 vov_history[i] = vov_till_ith; 324 vov_history[i] = vov_till_ith;
333 325
334 var_till_ith = var_till_ith / (ith + 1 - 326 var_till_ith = var_till_ith / (ith + 1 - 1);
335 var_history[i] = var_till_ith; 327 var_history[i] = var_till_ith;
336 sd_history[i] = std::sqrt(var_till_ith); 328 sd_history[i] = std::sqrt(var_till_ith);
337 r_history[i] = 329 r_history[i] =
338 std::sqrt(var_till_ith) / mean_till_ith 330 std::sqrt(var_till_ith) / mean_till_ith / std::sqrt(1.0 * (ith + 1));
339 331
340 if(std::fabs(cpu_time[ith]) > 0.0 && std:: 332 if(std::fabs(cpu_time[ith]) > 0.0 && std::fabs(r_history[i]) > 0.0)
341 { 333 {
342 fom_history[i] = 1.0 / std::pow(r_histor 334 fom_history[i] = 1.0 / std::pow(r_history[i], 2.0) / cpu_time[ith];
343 } 335 }
344 else 336 else
345 { 337 {
346 fom_history[i] = 0.0; 338 fom_history[i] = 0.0;
347 } 339 }
348 340
349 shift_till_ith += 341 shift_till_ith +=
350 ((ith + 1) - nonzero_till_ith) * std::po 342 ((ith + 1) - nonzero_till_ith) * std::pow(mean_till_ith, 3.0) * (-1.0);
351 shift_till_ith = shift_till_ith / (2 * v 343 shift_till_ith = shift_till_ith / (2 * var_till_ith * (ith + 1));
352 shift_history[i] = shift_till_ith; 344 shift_history[i] = shift_till_ith;
353 345
354 e_history[i] = 1.0 * nonzero_till_ith / (i 346 e_history[i] = 1.0 * nonzero_till_ith / (ith + 1);
355 if(std::fabs(e_history[i]) > 0.0) 347 if(std::fabs(e_history[i]) > 0.0)
356 { 348 {
357 r2eff_history[i] = (1 - e_history[i]) / 349 r2eff_history[i] = (1 - e_history[i]) / (e_history[i] * (ith + 1));
358 350
359 r2int_history[i] = (sum_x2_till_ith) / s 351 r2int_history[i] = (sum_x2_till_ith) / std::pow(sum_till_ith, 2.0) -
360 1 / (e_history[i] * ( 352 1 / (e_history[i] * (ith + 1));
361 } 353 }
362 } 354 }
363 } 355 }
364 356
365 void G4ConvergenceTester::ShowResult(std::ostr 357 void G4ConvergenceTester::ShowResult(std::ostream& out)
366 { 358 {
367 // if data has been added since the last com 359 // if data has been added since the last computation of the statistical values
368 // (not statsAreUpdated) call calStat to rec 360 // (not statsAreUpdated) call calStat to recompute the statistical values
369 if(!statsAreUpdated) 361 if(!statsAreUpdated)
370 { 362 {
371 calStat(); 363 calStat();
372 } 364 }
373 365
374 out << std::setprecision(6); 366 out << std::setprecision(6);
375 367
376 out << G4endl; 368 out << G4endl;
377 out << "G4ConvergenceTester Output Result of 369 out << "G4ConvergenceTester Output Result of " << name << G4endl;
378 out << std::setw(20) << "EFFICIENCY = " << s 370 out << std::setw(20) << "EFFICIENCY = " << std::setw(13) << efficiency
379 << G4endl; 371 << G4endl;
380 out << std::setw(20) << "MEAN = " << std::se 372 out << std::setw(20) << "MEAN = " << std::setw(13) << mean << G4endl;
381 out << std::setw(20) << "VAR = " << std::set 373 out << std::setw(20) << "VAR = " << std::setw(13) << var << G4endl;
382 out << std::setw(20) << "SD = " << std::setw 374 out << std::setw(20) << "SD = " << std::setw(13) << sd << G4endl;
383 out << std::setw(20) << "R = " << std::setw( 375 out << std::setw(20) << "R = " << std::setw(13) << r << G4endl;
384 out << std::setw(20) << "SHIFT = " << std::s 376 out << std::setw(20) << "SHIFT = " << std::setw(13) << shift << G4endl;
385 out << std::setw(20) << "VOV = " << std::set 377 out << std::setw(20) << "VOV = " << std::setw(13) << vov << G4endl;
386 out << std::setw(20) << "FOM = " << std::set 378 out << std::setw(20) << "FOM = " << std::setw(13) << fom << G4endl;
387 379
388 out << std::setw(20) << "THE LARGEST SCORE = 380 out << std::setw(20) << "THE LARGEST SCORE = " << std::setw(13) << largest
389 << " and it happened at " << largest_sco 381 << " and it happened at " << largest_score_happened << "th event"
390 << G4endl; 382 << G4endl;
391 if(mean != 0) 383 if(mean != 0)
392 { 384 {
393 out << std::setw(20) << "Affected Mean = " 385 out << std::setw(20) << "Affected Mean = " << std::setw(13) << mean_1
394 << " and its ratio to original is " << 386 << " and its ratio to original is " << mean_1 / mean << G4endl;
395 } 387 }
396 else 388 else
397 { 389 {
398 out << std::setw(20) << "Affected Mean = " 390 out << std::setw(20) << "Affected Mean = " << std::setw(13) << mean_1
399 << G4endl; 391 << G4endl;
400 } 392 }
401 if(var != 0) 393 if(var != 0)
402 { 394 {
403 out << std::setw(20) << "Affected VAR = " 395 out << std::setw(20) << "Affected VAR = " << std::setw(13) << var_1
404 << " and its ratio to original is " << 396 << " and its ratio to original is " << var_1 / var << G4endl;
405 } 397 }
406 else 398 else
407 { 399 {
408 out << std::setw(20) << "Affected VAR = " 400 out << std::setw(20) << "Affected VAR = " << std::setw(13) << var_1
409 << G4endl; 401 << G4endl;
410 } 402 }
411 if(r != 0) 403 if(r != 0)
412 { 404 {
413 out << std::setw(20) << "Affected R = " << 405 out << std::setw(20) << "Affected R = " << std::setw(13) << r_1
414 << " and its ratio to original is " << 406 << " and its ratio to original is " << r_1 / r << G4endl;
415 } 407 }
416 else 408 else
417 { 409 {
418 out << std::setw(20) << "Affected R = " << 410 out << std::setw(20) << "Affected R = " << std::setw(13) << r_1 << G4endl;
419 } 411 }
420 if(shift != 0) 412 if(shift != 0)
421 { 413 {
422 out << std::setw(20) << "Affected SHIFT = 414 out << std::setw(20) << "Affected SHIFT = " << std::setw(13) << shift_1
423 << " and its ratio to original is " << 415 << " and its ratio to original is " << shift_1 / shift << G4endl;
424 } 416 }
425 else 417 else
426 { 418 {
427 out << std::setw(20) << "Affected SHIFT = 419 out << std::setw(20) << "Affected SHIFT = " << std::setw(13) << shift_1
428 << G4endl; 420 << G4endl;
429 } 421 }
430 if(fom != 0) 422 if(fom != 0)
431 { 423 {
432 out << std::setw(20) << "Affected FOM = " 424 out << std::setw(20) << "Affected FOM = " << std::setw(13) << fom_1
433 << " and its ratio to original is " << 425 << " and its ratio to original is " << fom_1 / fom << G4endl;
434 } 426 }
435 else 427 else
436 { 428 {
437 out << std::setw(20) << "Affected FOM = " 429 out << std::setw(20) << "Affected FOM = " << std::setw(13) << fom_1
438 << G4endl; 430 << G4endl;
439 } 431 }
440 432
441 if(!showHistory) 433 if(!showHistory)
442 { 434 {
443 out << "Number of events of this run is to 435 out << "Number of events of this run is too small to do convergence tests."
444 << G4endl; 436 << G4endl;
445 return; 437 return;
446 } 438 }
447 439
448 check_stat_history(out); 440 check_stat_history(out);
449 441
450 // check SLOPE and output result 442 // check SLOPE and output result
451 if(calcSLOPE) 443 if(calcSLOPE)
452 { 444 {
453 if(slope >= 3) 445 if(slope >= 3)
454 { 446 {
455 noPass++; 447 noPass++;
456 out << "SLOPE is large enough" << G4endl 448 out << "SLOPE is large enough" << G4endl;
457 } 449 }
458 else 450 else
459 { 451 {
460 out << "SLOPE is not large enough" << G4 452 out << "SLOPE is not large enough" << G4endl;
461 } 453 }
462 } 454 }
463 else 455 else
464 { 456 {
465 out << "Number of non zero history too sma 457 out << "Number of non zero history too small to calculate SLOPE" << G4endl;
466 } 458 }
467 459
468 out << "This result passes " << noPass << " 460 out << "This result passes " << noPass << " / " << noTotal
469 << " Convergence Test." << G4endl; 461 << " Convergence Test." << G4endl;
470 out << G4endl; 462 out << G4endl;
471 } 463 }
472 464
473 void G4ConvergenceTester::ShowHistory(std::ost 465 void G4ConvergenceTester::ShowHistory(std::ostream& out)
474 { 466 {
475 if(!showHistory) 467 if(!showHistory)
476 { 468 {
477 out << "Number of events of this run is to 469 out << "Number of events of this run is too small to show history."
478 << G4endl; 470 << G4endl;
479 return; 471 return;
480 } 472 }
481 473
482 out << std::setprecision(6); 474 out << std::setprecision(6);
483 475
484 out << G4endl; 476 out << G4endl;
485 out << "G4ConvergenceTester Output History o 477 out << "G4ConvergenceTester Output History of " << name << G4endl;
486 out << "i/" << noBinOfHistory << " till_ith 478 out << "i/" << noBinOfHistory << " till_ith mean" << std::setw(13)
487 << "var" << std::setw(13) << "sd" << std 479 << "var" << std::setw(13) << "sd" << std::setw(13) << "r" << std::setw(13)
488 << "vov" << std::setw(13) << "fom" << st 480 << "vov" << std::setw(13) << "fom" << std::setw(13) << "shift"
489 << std::setw(13) << "e" << std::setw(13) 481 << std::setw(13) << "e" << std::setw(13) << "r2eff" << std::setw(13)
490 << "r2int" << G4endl; 482 << "r2int" << G4endl;
491 for(G4int i = 1; i <= noBinOfHistory; i++) 483 for(G4int i = 1; i <= noBinOfHistory; i++)
492 { 484 {
493 out << std::setw(4) << i << " " << std::se 485 out << std::setw(4) << i << " " << std::setw(5) << history_grid[i - 1]
494 << std::setw(13) << mean_history[i - 1 486 << std::setw(13) << mean_history[i - 1] << std::setw(13)
495 << var_history[i - 1] << std::setw(13) 487 << var_history[i - 1] << std::setw(13) << sd_history[i - 1]
496 << std::setw(13) << r_history[i - 1] < 488 << std::setw(13) << r_history[i - 1] << std::setw(13)
497 << vov_history[i - 1] << std::setw(13) 489 << vov_history[i - 1] << std::setw(13) << fom_history[i - 1]
498 << std::setw(13) << shift_history[i - 490 << std::setw(13) << shift_history[i - 1] << std::setw(13)
499 << e_history[i - 1] << std::setw(13) < 491 << e_history[i - 1] << std::setw(13) << r2eff_history[i - 1]
500 << std::setw(13) << r2int_history[i - 492 << std::setw(13) << r2int_history[i - 1] << G4endl;
501 } 493 }
502 } 494 }
503 495
504 void G4ConvergenceTester::check_stat_history(s 496 void G4ConvergenceTester::check_stat_history(std::ostream& out)
505 { 497 {
506 // 1 sigma rejection for null hypothesis 498 // 1 sigma rejection for null hypothesis
507 499
508 std::vector<G4double> first_ally; 500 std::vector<G4double> first_ally;
509 std::vector<G4double> second_ally; 501 std::vector<G4double> second_ally;
510 502
511 // use 2nd half of hisories 503 // use 2nd half of hisories
512 std::size_t N = mean_history.size() / 2; << 504 G4int N = mean_history.size() / 2;
513 std::size_t i; << 505 G4int i;
514 506
515 G4double pearson_r; 507 G4double pearson_r;
516 G4double t; 508 G4double t;
517 509
518 first_ally.resize(N); 510 first_ally.resize(N);
519 second_ally.resize(N); 511 second_ally.resize(N);
520 512
521 G4double sum_of_var = 513 G4double sum_of_var =
522 std::accumulate(var_history.begin(), var_h 514 std::accumulate(var_history.begin(), var_history.end(), 0.0);
523 if(sum_of_var == 0.0) 515 if(sum_of_var == 0.0)
524 { 516 {
525 out << "Variances in all historical grids 517 out << "Variances in all historical grids are zero." << G4endl;
526 out << "Terminating checking behavior of s 518 out << "Terminating checking behavior of statistics numbers." << G4endl;
527 return; 519 return;
528 } 520 }
529 521
530 // Mean 522 // Mean
531 523
532 for(i = 0; i < N; ++i) 524 for(i = 0; i < N; ++i)
533 { 525 {
534 first_ally[i] = history_grid[N + i]; 526 first_ally[i] = history_grid[N + i];
535 second_ally[i] = mean_history[N + i]; 527 second_ally[i] = mean_history[N + i];
536 } 528 }
537 529
538 pearson_r = calc_Pearson_r((G4int)N, first_a << 530 pearson_r = calc_Pearson_r(N, first_ally, second_ally);
539 t = pearson_r * std::sqrt((N - 2) / 531 t = pearson_r * std::sqrt((N - 2) / (1 - pearson_r * pearson_r));
540 532
541 if(t < 0.429318) // Student t of (Degree of 533 if(t < 0.429318) // Student t of (Degree of freedom = N-2 )
542 { 534 {
543 out << "MEAN distribution is RANDOM" << G 535 out << "MEAN distribution is RANDOM" << G4endl;
544 noPass++; 536 noPass++;
545 } 537 }
546 else 538 else
547 { 539 {
548 out << "MEAN distribution is not RANDOM" < 540 out << "MEAN distribution is not RANDOM" << G4endl;
549 } 541 }
550 542
551 // R 543 // R
552 544
553 for(i = 0; i < N; ++i) 545 for(i = 0; i < N; ++i)
554 { 546 {
555 first_ally[i] = 1.0 / std::sqrt(G4double( 547 first_ally[i] = 1.0 / std::sqrt(G4double(history_grid[N + i]));
556 second_ally[i] = r_history[N + i]; 548 second_ally[i] = r_history[N + i];
557 } 549 }
558 550
559 pearson_r = calc_Pearson_r(G4int(N), first_a << 551 pearson_r = calc_Pearson_r(N, first_ally, second_ally);
560 t = pearson_r * std::sqrt((N - 2) / 552 t = pearson_r * std::sqrt((N - 2) / (1 - pearson_r * pearson_r));
561 553
562 if(t > 1.090546) 554 if(t > 1.090546)
563 { 555 {
564 out << "r follows 1/std::sqrt(N)" << G4end 556 out << "r follows 1/std::sqrt(N)" << G4endl;
565 noPass++; 557 noPass++;
566 } 558 }
567 else 559 else
568 { 560 {
569 out << "r does not follow 1/std::sqrt(N)" 561 out << "r does not follow 1/std::sqrt(N)" << G4endl;
570 } 562 }
571 563
572 if(is_monotonically_decrease(second_ally)) << 564 if(is_monotonically_decrease(second_ally) == true)
573 { 565 {
574 out << "r is monotonically decrease " << G 566 out << "r is monotonically decrease " << G4endl;
575 } 567 }
576 else 568 else
577 { 569 {
578 out << "r is NOT monotonically decrease " 570 out << "r is NOT monotonically decrease " << G4endl;
579 } 571 }
580 572
581 if(r_history.back() < 0.1) 573 if(r_history.back() < 0.1)
582 { 574 {
583 out << "r is less than 0.1. r = " << r_his 575 out << "r is less than 0.1. r = " << r_history.back() << G4endl;
584 noPass++; 576 noPass++;
585 } 577 }
586 else 578 else
587 { 579 {
588 out << "r is NOT less than 0.1. r = " << r 580 out << "r is NOT less than 0.1. r = " << r_history.back() << G4endl;
589 } 581 }
590 582
591 // VOV 583 // VOV
592 for(i = 0; i < N; ++i) 584 for(i = 0; i < N; ++i)
593 { 585 {
594 first_ally[i] = 1.0 / history_grid[N + i] 586 first_ally[i] = 1.0 / history_grid[N + i];
595 second_ally[i] = vov_history[N + i]; 587 second_ally[i] = vov_history[N + i];
596 } 588 }
597 589
598 pearson_r = calc_Pearson_r(G4int(N), first_a << 590 pearson_r = calc_Pearson_r(N, first_ally, second_ally);
599 t = pearson_r * std::sqrt((N - 2) / 591 t = pearson_r * std::sqrt((N - 2) / (1 - pearson_r * pearson_r));
600 592
601 if(t > 1.090546) 593 if(t > 1.090546)
602 { 594 {
603 out << "VOV follows 1/std::sqrt(N)" << G4e 595 out << "VOV follows 1/std::sqrt(N)" << G4endl;
604 noPass++; 596 noPass++;
605 } 597 }
606 else 598 else
607 { 599 {
608 out << "VOV does not follow 1/std::sqrt(N) 600 out << "VOV does not follow 1/std::sqrt(N)" << G4endl;
609 } 601 }
610 602
611 if(is_monotonically_decrease(second_ally)) << 603 if(is_monotonically_decrease(second_ally) == true)
612 { 604 {
613 out << "VOV is monotonically decrease " << 605 out << "VOV is monotonically decrease " << G4endl;
614 } 606 }
615 else 607 else
616 { 608 {
617 out << "VOV is NOT monotonically decrease 609 out << "VOV is NOT monotonically decrease " << G4endl;
618 } 610 }
619 611
620 // FOM 612 // FOM
621 613
622 for(i = 0; i < N; ++i) 614 for(i = 0; i < N; ++i)
623 { 615 {
624 first_ally[i] = history_grid[N + i]; 616 first_ally[i] = history_grid[N + i];
625 second_ally[i] = fom_history[N + i]; 617 second_ally[i] = fom_history[N + i];
626 } 618 }
627 619
628 pearson_r = calc_Pearson_r(G4int(N), std::mo << 620 pearson_r = calc_Pearson_r(N, first_ally, second_ally);
629 t = pearson_r * std::sqrt((N - 2) / 621 t = pearson_r * std::sqrt((N - 2) / (1 - pearson_r * pearson_r));
630 622
631 if(t < 0.429318) 623 if(t < 0.429318)
632 { 624 {
633 out << "FOM distribution is RANDOM" << G4e 625 out << "FOM distribution is RANDOM" << G4endl;
634 noPass++; 626 noPass++;
635 } 627 }
636 else 628 else
637 { 629 {
638 out << "FOM distribution is not RANDOM" << 630 out << "FOM distribution is not RANDOM" << G4endl;
639 } 631 }
640 } 632 }
641 633
642 G4double G4ConvergenceTester::calc_Pearson_r(G 634 G4double G4ConvergenceTester::calc_Pearson_r(G4int N,
643 s 635 std::vector<G4double> first_ally,
644 s 636 std::vector<G4double> second_ally)
645 { 637 {
646 G4double first_mean = 0.0; 638 G4double first_mean = 0.0;
647 G4double second_mean = 0.0; 639 G4double second_mean = 0.0;
648 640
649 G4int i; 641 G4int i;
650 for(i = 0; i < N; i++) 642 for(i = 0; i < N; i++)
651 { 643 {
652 first_mean += first_ally[i]; 644 first_mean += first_ally[i];
653 second_mean += second_ally[i]; 645 second_mean += second_ally[i];
654 } 646 }
655 first_mean = first_mean / N; 647 first_mean = first_mean / N;
656 second_mean = second_mean / N; 648 second_mean = second_mean / N;
657 649
658 G4double a = 0.0; 650 G4double a = 0.0;
659 for(i = 0; i < N; ++i) 651 for(i = 0; i < N; ++i)
660 { 652 {
661 a += (first_ally[i] - first_mean) * (secon 653 a += (first_ally[i] - first_mean) * (second_ally[i] - second_mean);
662 } 654 }
663 655
664 G4double b1 = 0.0; 656 G4double b1 = 0.0;
665 G4double b2 = 0.0; 657 G4double b2 = 0.0;
666 for(i = 0; i < N; ++i) 658 for(i = 0; i < N; ++i)
667 { 659 {
668 b1 += (first_ally[i] - first_mean) * (firs 660 b1 += (first_ally[i] - first_mean) * (first_ally[i] - first_mean);
669 b2 += (second_ally[i] - second_mean) * (se 661 b2 += (second_ally[i] - second_mean) * (second_ally[i] - second_mean);
670 } 662 }
671 663
672 G4double rds = a / std::sqrt(b1 * b2); 664 G4double rds = a / std::sqrt(b1 * b2);
673 665
674 return rds; 666 return rds;
675 } 667 }
676 668
677 G4bool G4ConvergenceTester::is_monotonically_d 669 G4bool G4ConvergenceTester::is_monotonically_decrease(
678 const std::vector<G4double>& ally) << 670 std::vector<G4double> ally)
679 { 671 {
680 for(auto it = ally.cbegin(); it != ally.cend 672 for(auto it = ally.cbegin(); it != ally.cend() - 1; ++it)
681 { 673 {
682 if(*it < *(it + 1)) 674 if(*it < *(it + 1))
683 { <<
684 return FALSE; 675 return FALSE;
685 } <<
686 } 676 }
687 677
688 ++noPass; 678 ++noPass;
689 return TRUE; 679 return TRUE;
690 } 680 }
691 681
692 void G4ConvergenceTester::calc_slope_fit(const << 682 void G4ConvergenceTester::calc_slope_fit(std::vector<G4double>)
693 { 683 {
694 // create PDF bins 684 // create PDF bins
695 G4double max = largest_scores.front(); 685 G4double max = largest_scores.front();
696 auto last = G4int(largest_scores.size()); << 686 G4int last = G4int(largest_scores.size());
697 G4double min = 0.0; 687 G4double min = 0.0;
698 if(largest_scores.back() != 0) 688 if(largest_scores.back() != 0)
699 { 689 {
700 min = largest_scores.back(); 690 min = largest_scores.back();
701 } 691 }
702 else 692 else
703 { 693 {
704 min = largest_scores[last - 1]; 694 min = largest_scores[last - 1];
705 last = last - 1; 695 last = last - 1;
706 } 696 }
707 697
708 if(max * 0.99 < min) 698 if(max * 0.99 < min)
709 { 699 {
710 // upper limit is assumed to have been rea 700 // upper limit is assumed to have been reached
711 slope = 10.0; 701 slope = 10.0;
712 return; 702 return;
713 } 703 }
714 704
715 std::vector<G4double> pdf_grid; 705 std::vector<G4double> pdf_grid;
716 706
717 pdf_grid.resize(noBinOfPDF + 1); // no grid 707 pdf_grid.resize(noBinOfPDF + 1); // no grid = no bins + 1
718 pdf_grid[0] = max; 708 pdf_grid[0] = max;
719 pdf_grid[noBinOfPDF] = min; 709 pdf_grid[noBinOfPDF] = min;
720 G4double log10_max = std::log10(max); 710 G4double log10_max = std::log10(max);
721 G4double log10_min = std::log10(min); 711 G4double log10_min = std::log10(min);
722 G4double log10_delta = log10_max - log10_min 712 G4double log10_delta = log10_max - log10_min;
723 for(G4int i = 1; i < noBinOfPDF; ++i) 713 for(G4int i = 1; i < noBinOfPDF; ++i)
724 { 714 {
725 pdf_grid[i] = std::pow(10.0, log10_max - l 715 pdf_grid[i] = std::pow(10.0, log10_max - log10_delta / 10.0 * (i));
726 } 716 }
727 717
728 std::vector<G4double> pdf; 718 std::vector<G4double> pdf;
729 pdf.resize(noBinOfPDF); 719 pdf.resize(noBinOfPDF);
730 720
731 for(G4int j = 0; j < last; ++j) 721 for(G4int j = 0; j < last; ++j)
732 { 722 {
733 for(G4int i = 0; i < 11; ++i) 723 for(G4int i = 0; i < 11; ++i)
734 { 724 {
735 if(largest_scores[j] >= pdf_grid[i + 1]) 725 if(largest_scores[j] >= pdf_grid[i + 1])
736 { 726 {
737 pdf[i] += 1.0 / (pdf_grid[i] - pdf_gri 727 pdf[i] += 1.0 / (pdf_grid[i] - pdf_grid[i + 1]) / n;
738 break; 728 break;
739 } 729 }
740 } 730 }
741 } 731 }
742 732
743 f_xi.resize(noBinOfPDF); 733 f_xi.resize(noBinOfPDF);
744 f_yi.resize(noBinOfPDF); 734 f_yi.resize(noBinOfPDF);
745 for(G4int i = 0; i < noBinOfPDF; ++i) 735 for(G4int i = 0; i < noBinOfPDF; ++i)
746 { 736 {
747 f_xi[i] = (pdf_grid[i] + pdf_grid[i + 1]) 737 f_xi[i] = (pdf_grid[i] + pdf_grid[i + 1]) / 2;
748 f_yi[i] = pdf[i]; 738 f_yi[i] = pdf[i];
749 } 739 }
750 740
751 // number of variables ( a and k ) 741 // number of variables ( a and k )
752 minimizer = new G4SimplexDownhill<G4Converge 742 minimizer = new G4SimplexDownhill<G4ConvergenceTester>(this, 2);
753 // G4double minimum = minimizer->GetMinimum 743 // G4double minimum = minimizer->GetMinimum();
754 std::vector<G4double> mp = minimizer->GetMin 744 std::vector<G4double> mp = minimizer->GetMinimumPoint();
755 G4double k = mp[1]; 745 G4double k = mp[1];
756 746
757 // G4cout << "SLOPE " << 1/mp[1]+1 << G4endl 747 // G4cout << "SLOPE " << 1/mp[1]+1 << G4endl;
758 // G4cout << "SLOPE a " << mp[0] << G4endl; 748 // G4cout << "SLOPE a " << mp[0] << G4endl;
759 // G4cout << "SLOPE k " << mp[1] << G4endl; 749 // G4cout << "SLOPE k " << mp[1] << G4endl;
760 // G4cout << "SLOPE minimum " << minimizer- 750 // G4cout << "SLOPE minimum " << minimizer->GetMinimum() << G4endl;
761 751
762 slope = 1 / mp[1] + 1; 752 slope = 1 / mp[1] + 1;
763 if(k < 1.0 / 9) // Please look Pareto distr 753 if(k < 1.0 / 9) // Please look Pareto distribution with "sigma=a" and "k"
764 { 754 {
765 slope = 10; 755 slope = 10;
766 } 756 }
767 if(slope > 10) 757 if(slope > 10)
768 { 758 {
769 slope = 10; 759 slope = 10;
770 } 760 }
771 } 761 }
772 762
773 G4double G4ConvergenceTester::slope_fitting_fu 763 G4double G4ConvergenceTester::slope_fitting_function(std::vector<G4double> x)
774 { 764 {
775 G4double a = x[0]; 765 G4double a = x[0];
776 G4double k = x[1]; 766 G4double k = x[1];
777 767
778 if(a <= 0) 768 if(a <= 0)
779 { 769 {
780 return 3.402823466e+38; // FLOAT_MAX 770 return 3.402823466e+38; // FLOAT_MAX
781 } 771 }
782 if(k == 0) 772 if(k == 0)
783 { 773 {
784 return 3.402823466e+38; // FLOAT_MAX 774 return 3.402823466e+38; // FLOAT_MAX
785 } 775 }
786 776
787 // f_xi and f_yi is filled at "calc_slope_fi 777 // f_xi and f_yi is filled at "calc_slope_fit"
788 778
789 G4double y = 0.0; 779 G4double y = 0.0;
790 for(G4int i = 0; i < G4int(f_yi.size()); ++i 780 for(G4int i = 0; i < G4int(f_yi.size()); ++i)
791 { 781 {
792 // if ( 1/a * ( 1 + k * f_xi [ i ] / a ) < 782 // if ( 1/a * ( 1 + k * f_xi [ i ] / a ) < 0 )
793 if((1 + k * f_xi[i] / a) < 0) 783 if((1 + k * f_xi[i] / a) < 0)
794 { 784 {
795 y += 3.402823466e+38; // FLOAT_MAX 785 y += 3.402823466e+38; // FLOAT_MAX
796 } 786 }
797 else 787 else
798 { 788 {
799 y += (f_yi[i] - 1 / a * std::pow(1 + k * 789 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 * 790 (f_yi[i] - 1 / a * std::pow(1 + k * f_xi[i] / a, -1 / k - 1));
801 } 791 }
802 } 792 }
803 793
804 return y; 794 return y;
805 } 795 }
806 796