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