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 // G4SimplexDownhill inline methods implementa 26 // G4SimplexDownhill inline methods implementation 27 // 27 // 28 // Author: Tatsumi Koi (SLAC/SCCS), 2007 28 // Author: Tatsumi Koi (SLAC/SCCS), 2007 29 // ------------------------------------------- 29 // -------------------------------------------------------------------------- 30 30 31 #include <cfloat> 31 #include <cfloat> 32 #include <iostream> 32 #include <iostream> 33 #include <numeric> 33 #include <numeric> 34 34 35 template <class T> 35 template <class T> 36 void G4SimplexDownhill<T>::init() 36 void G4SimplexDownhill<T>::init() 37 { 37 { 38 alpha = 2.0; // refrection coefficient: 0 38 alpha = 2.0; // refrection coefficient: 0 < alpha 39 beta = 0.5; // contraction coefficient: 39 beta = 0.5; // contraction coefficient: 0 < beta < 1 40 gamma = 2.0; // expantion coefficient: 1 < 40 gamma = 2.0; // expantion coefficient: 1 < gamma 41 41 42 maximum_no_trial = 10000; 42 maximum_no_trial = 10000; 43 max_se = FLT_MIN; 43 max_se = FLT_MIN; 44 // max_ratio = FLT_EPSILON/1; 44 // max_ratio = FLT_EPSILON/1; 45 max_ratio = DBL_EPSILON / 1; 45 max_ratio = DBL_EPSILON / 1; 46 minimized = false; 46 minimized = false; 47 } 47 } 48 48 49 /* 49 /* 50 50 51 void G4SimplexDownhill<class T>:: 51 void G4SimplexDownhill<class T>:: 52 SetFunction( G4int n , G4double( *afunc )( std 52 SetFunction( G4int n , G4double( *afunc )( std::vector < G4double > ) ) 53 { 53 { 54 numberOfVariable = n; 54 numberOfVariable = n; 55 theFunction = afunc; 55 theFunction = afunc; 56 minimized = false; 56 minimized = false; 57 } 57 } 58 58 59 */ 59 */ 60 60 61 template <class T> 61 template <class T> 62 G4double G4SimplexDownhill<T>::GetMinimum() 62 G4double G4SimplexDownhill<T>::GetMinimum() 63 { 63 { 64 initialize(); 64 initialize(); 65 65 66 // First Tryal; 66 // First Tryal; 67 67 68 // G4cout << "Begin First Trials" << G4endl; 68 // G4cout << "Begin First Trials" << G4endl; 69 doDownhill(); 69 doDownhill(); 70 // G4cout << "End First Trials" << G4endl; 70 // G4cout << "End First Trials" << G4endl; 71 71 72 auto it_minh = std::min_element(currentHeigh << 72 std::vector<G4double>::const_iterator it_minh = >> 73 std::min_element(currentHeights.cbegin(), currentHeights.cend()); 73 G4int imin = 0; 74 G4int imin = 0; 74 G4int i = 0; 75 G4int i = 0; 75 for(auto it = currentHeights.cbegin(); it != 76 for(auto it = currentHeights.cbegin(); it != currentHeights.cend(); ++it) 76 { 77 { 77 if(it == it_minh) 78 if(it == it_minh) 78 { 79 { 79 imin = i; 80 imin = i; 80 } 81 } 81 ++i; 82 ++i; 82 } 83 } 83 minimumPoint = currentSimplex[imin]; 84 minimumPoint = currentSimplex[imin]; 84 85 85 // Second Trial 86 // Second Trial 86 87 87 // std::vector< G4double > minimumPoint = cu 88 // std::vector< G4double > minimumPoint = currentSimplex[ 0 ]; 88 initialize(); 89 initialize(); 89 90 90 currentSimplex[numberOfVariable] = minimumPo 91 currentSimplex[numberOfVariable] = minimumPoint; 91 92 92 // G4cout << "Begin Second Trials" << G4endl 93 // G4cout << "Begin Second Trials" << G4endl; 93 doDownhill(); 94 doDownhill(); 94 // G4cout << "End Second Trials" << G4endl; 95 // G4cout << "End Second Trials" << G4endl; 95 96 96 G4double sum = 97 G4double sum = 97 std::accumulate(currentHeights.begin(), cu 98 std::accumulate(currentHeights.begin(), currentHeights.end(), 0.0); 98 G4double average = sum / (numberOfVariable + 99 G4double average = sum / (numberOfVariable + 1); 99 G4double minimum = average; 100 G4double minimum = average; 100 101 101 minimized = true; 102 minimized = true; 102 103 103 return minimum; 104 return minimum; 104 } 105 } 105 106 106 template <class T> 107 template <class T> 107 void G4SimplexDownhill<T>::initialize() 108 void G4SimplexDownhill<T>::initialize() 108 { 109 { 109 currentSimplex.resize(numberOfVariable + 1); 110 currentSimplex.resize(numberOfVariable + 1); 110 currentHeights.resize(numberOfVariable + 1); 111 currentHeights.resize(numberOfVariable + 1); 111 112 112 for(G4int i = 0; i < numberOfVariable; ++i) 113 for(G4int i = 0; i < numberOfVariable; ++i) 113 { 114 { 114 std::vector<G4double> avec(numberOfVariabl 115 std::vector<G4double> avec(numberOfVariable, 0.0); 115 avec[i] = 1.0; 116 avec[i] = 1.0; 116 currentSimplex[i] = std::move(avec); << 117 currentSimplex[i] = avec; 117 } 118 } 118 119 119 // std::vector< G4double > avec ( numberOfVa 120 // std::vector< G4double > avec ( numberOfVariable , 0.0 ); 120 std::vector<G4double> avec(numberOfVariable, 121 std::vector<G4double> avec(numberOfVariable, 1); 121 currentSimplex[numberOfVariable] = std::move << 122 currentSimplex[numberOfVariable] = avec; 122 } 123 } 123 124 124 template <class T> 125 template <class T> 125 void G4SimplexDownhill<T>::calHeights() 126 void G4SimplexDownhill<T>::calHeights() 126 { 127 { 127 for(G4int i = 0; i <= numberOfVariable; ++i) 128 for(G4int i = 0; i <= numberOfVariable; ++i) 128 { 129 { 129 currentHeights[i] = getValue(currentSimple 130 currentHeights[i] = getValue(currentSimplex[i]); 130 } 131 } 131 } 132 } 132 133 133 template <class T> 134 template <class T> 134 std::vector<G4double> G4SimplexDownhill<T>::ca 135 std::vector<G4double> G4SimplexDownhill<T>::calCentroid(G4int ih) 135 { 136 { 136 std::vector<G4double> centroid(numberOfVaria 137 std::vector<G4double> centroid(numberOfVariable, 0.0); 137 138 138 G4int i = 0; 139 G4int i = 0; 139 for(const auto & it : currentSimplex) << 140 for(auto it = currentSimplex.cbegin(); it != currentSimplex.cend(); ++it) 140 { 141 { 141 if(i != ih) 142 if(i != ih) 142 { 143 { 143 for(G4int j = 0; j < numberOfVariable; + 144 for(G4int j = 0; j < numberOfVariable; ++j) 144 { 145 { 145 centroid[j] += it[j] / numberOfVariabl << 146 centroid[j] += (*it)[j] / numberOfVariable; 146 } 147 } 147 } 148 } 148 ++i; 149 ++i; 149 } 150 } 150 151 151 return centroid; 152 return centroid; 152 } 153 } 153 154 154 template <class T> 155 template <class T> 155 std::vector<G4double> G4SimplexDownhill<T>::ge 156 std::vector<G4double> G4SimplexDownhill<T>::getReflectionPoint( 156 std::vector<G4double> p, std::vector<G4doubl 157 std::vector<G4double> p, std::vector<G4double> centroid) 157 { 158 { 158 // G4cout << "Reflection" << G4endl; 159 // G4cout << "Reflection" << G4endl; 159 160 160 std::vector<G4double> reflectionP(numberOfVa 161 std::vector<G4double> reflectionP(numberOfVariable, 0.0); 161 162 162 for(G4int i = 0; i < numberOfVariable; ++i) 163 for(G4int i = 0; i < numberOfVariable; ++i) 163 { 164 { 164 reflectionP[i] = (1 + alpha) * centroid[i] 165 reflectionP[i] = (1 + alpha) * centroid[i] - alpha * p[i]; 165 } 166 } 166 167 167 return reflectionP; 168 return reflectionP; 168 } 169 } 169 170 170 template <class T> 171 template <class T> 171 std::vector<G4double> G4SimplexDownhill<T>::ge 172 std::vector<G4double> G4SimplexDownhill<T>::getExpansionPoint( 172 std::vector<G4double> p, std::vector<G4doubl 173 std::vector<G4double> p, std::vector<G4double> centroid) 173 { 174 { 174 // G4cout << "Expantion" << G4endl; 175 // G4cout << "Expantion" << G4endl; 175 176 176 std::vector<G4double> expansionP(numberOfVar 177 std::vector<G4double> expansionP(numberOfVariable, 0.0); 177 178 178 for(G4int i = 0; i < numberOfVariable; ++i) 179 for(G4int i = 0; i < numberOfVariable; ++i) 179 { 180 { 180 expansionP[i] = (1 - gamma) * centroid[i] 181 expansionP[i] = (1 - gamma) * centroid[i] + gamma * p[i]; 181 } 182 } 182 183 183 return expansionP; 184 return expansionP; 184 } 185 } 185 186 186 template <class T> 187 template <class T> 187 std::vector<G4double> G4SimplexDownhill<T>::ge 188 std::vector<G4double> G4SimplexDownhill<T>::getContractionPoint( 188 std::vector<G4double> p, std::vector<G4doubl 189 std::vector<G4double> p, std::vector<G4double> centroid) 189 { 190 { 190 std::vector<G4double> contractionP(numberOfV 191 std::vector<G4double> contractionP(numberOfVariable, 0.0); 191 192 192 for(G4int i = 0; i < numberOfVariable; ++i) 193 for(G4int i = 0; i < numberOfVariable; ++i) 193 { 194 { 194 contractionP[i] = (1 - beta) * centroid[i] 195 contractionP[i] = (1 - beta) * centroid[i] + beta * p[i]; 195 } 196 } 196 197 197 return contractionP; 198 return contractionP; 198 } 199 } 199 200 200 template <class T> 201 template <class T> 201 G4bool G4SimplexDownhill<T>::isItGoodEnough() 202 G4bool G4SimplexDownhill<T>::isItGoodEnough() 202 { 203 { >> 204 G4bool result = false; >> 205 203 G4double sum = 206 G4double sum = 204 std::accumulate(currentHeights.begin(), cu 207 std::accumulate(currentHeights.begin(), currentHeights.end(), 0.0); 205 G4double average = sum / (numberOfVariable + 208 G4double average = sum / (numberOfVariable + 1); 206 209 207 G4double delta = 0.0; 210 G4double delta = 0.0; 208 for(G4int i = 0; i <= numberOfVariable; ++i) 211 for(G4int i = 0; i <= numberOfVariable; ++i) 209 { 212 { 210 delta += std::abs(currentHeights[i] - aver 213 delta += std::abs(currentHeights[i] - average); 211 } 214 } 212 215 213 G4bool result = false; << 216 if(delta / (numberOfVariable + 1) / average < max_ratio) 214 if (average > 0.0) << 215 { 217 { 216 result = ((delta / (numberOfVariable + 1) << 218 result = true; 217 } 219 } >> 220 218 return result; 221 return result; 219 } 222 } 220 223 221 template <class T> 224 template <class T> 222 void G4SimplexDownhill<T>::doDownhill() 225 void G4SimplexDownhill<T>::doDownhill() 223 { 226 { 224 G4int nth_trial = 0; 227 G4int nth_trial = 0; 225 228 226 while(nth_trial < maximum_no_trial) 229 while(nth_trial < maximum_no_trial) 227 { 230 { 228 calHeights(); 231 calHeights(); 229 232 230 if(isItGoodEnough()) 233 if(isItGoodEnough()) 231 { 234 { 232 break; 235 break; 233 } 236 } 234 237 235 auto it_maxh = std::max_element(currentHei << 238 std::vector<G4double>::const_iterator it_maxh = 236 auto it_minh = std::min_element(currentHei << 239 std::max_element(currentHeights.cbegin(), currentHeights.cend()); >> 240 std::vector<G4double>::const_iterator it_minh = >> 241 std::min_element(currentHeights.cbegin(), currentHeights.cend()); 237 242 238 G4double h_H = *it_maxh; 243 G4double h_H = *it_maxh; 239 G4double h_L = *it_minh; 244 G4double h_L = *it_minh; 240 245 241 G4int ih = 0; 246 G4int ih = 0; 242 G4int il = 0; 247 G4int il = 0; 243 G4double h_H2 = 0.0; 248 G4double h_H2 = 0.0; 244 G4int i = 0; 249 G4int i = 0; 245 for(auto it = currentHeights.cbegin(); it 250 for(auto it = currentHeights.cbegin(); it != currentHeights.cend(); ++it) 246 { 251 { 247 if(it == it_maxh) 252 if(it == it_maxh) 248 { 253 { 249 ih = i; 254 ih = i; 250 } 255 } 251 else 256 else 252 { 257 { 253 h_H2 = std::max(h_H2, *it); 258 h_H2 = std::max(h_H2, *it); 254 } 259 } 255 260 256 if(it == it_minh) 261 if(it == it_minh) 257 { 262 { 258 il = i; 263 il = i; 259 } 264 } 260 ++i; 265 ++i; 261 } 266 } 262 267 263 std::vector<G4double> centroidPoint = calC 268 std::vector<G4double> centroidPoint = calCentroid(ih); 264 269 265 // REFLECTION 270 // REFLECTION 266 std::vector<G4double> reflectionPoint = 271 std::vector<G4double> reflectionPoint = 267 getReflectionPoint(currentSimplex[ih], c 272 getReflectionPoint(currentSimplex[ih], centroidPoint); 268 273 269 G4double h = getValue(reflectionPoint); 274 G4double h = getValue(reflectionPoint); 270 275 271 if(h <= h_L) 276 if(h <= h_L) 272 { 277 { 273 // EXPANSION 278 // EXPANSION 274 std::vector<G4double> expansionPoint = 279 std::vector<G4double> expansionPoint = 275 getExpansionPoint(reflectionPoint, std << 280 getExpansionPoint(reflectionPoint, centroidPoint); 276 G4double hh = getValue(expansionPoint); 281 G4double hh = getValue(expansionPoint); 277 282 278 if(hh <= h_L) 283 if(hh <= h_L) 279 { 284 { 280 // Replace 285 // Replace 281 currentSimplex[ih] = std::move(expansi << 286 currentSimplex[ih] = expansionPoint; 282 // G4cout << "A" << G4endl; 287 // G4cout << "A" << G4endl; 283 } 288 } 284 else 289 else 285 { 290 { 286 // Replace 291 // Replace 287 currentSimplex[ih] = std::move(reflect << 292 currentSimplex[ih] = reflectionPoint; 288 // G4cout << "B1" << G4endl; 293 // G4cout << "B1" << G4endl; 289 } 294 } 290 } 295 } 291 else 296 else 292 { 297 { 293 if(h <= h_H2) 298 if(h <= h_H2) 294 { 299 { 295 // Replace 300 // Replace 296 currentSimplex[ih] = std::move(reflect << 301 currentSimplex[ih] = reflectionPoint; 297 // G4cout << "B2" << G4endl; 302 // G4cout << "B2" << G4endl; 298 } 303 } 299 else 304 else 300 { 305 { 301 if(h <= h_H) 306 if(h <= h_H) 302 { 307 { 303 // Replace 308 // Replace 304 currentSimplex[ih] = std::move(refle << 309 currentSimplex[ih] = reflectionPoint; 305 // G4cout << "BC" << G4endl; 310 // G4cout << "BC" << G4endl; 306 } 311 } 307 // CONTRACTION 312 // CONTRACTION 308 std::vector<G4double> contractionPoint 313 std::vector<G4double> contractionPoint = 309 getContractionPoint(currentSimplex[i << 314 getContractionPoint(currentSimplex[ih], centroidPoint); 310 G4double hh = getValue(contractionPoin 315 G4double hh = getValue(contractionPoint); 311 if(hh <= h_H) 316 if(hh <= h_H) 312 { 317 { 313 // Replace 318 // Replace 314 currentSimplex[ih] = std::move(contr << 319 currentSimplex[ih] = contractionPoint; 315 // G4cout << "C" << G4endl; 320 // G4cout << "C" << G4endl; 316 } 321 } 317 else 322 else 318 { 323 { 319 // Replace 324 // Replace 320 for(G4int j = 0; j <= numberOfVariab 325 for(G4int j = 0; j <= numberOfVariable; ++j) 321 { 326 { 322 std::vector<G4double> vec(numberOf 327 std::vector<G4double> vec(numberOfVariable, 0.0); 323 for(G4int k = 0; k < numberOfVaria 328 for(G4int k = 0; k < numberOfVariable; ++k) 324 { 329 { 325 vec[k] = (currentSimplex[j][k] + 330 vec[k] = (currentSimplex[j][k] + currentSimplex[il][k]) / 2.0; 326 } 331 } 327 currentSimplex[j] = std::move(vec) << 332 currentSimplex[j] = vec; 328 } 333 } 329 } 334 } 330 } 335 } 331 } 336 } 332 337 333 ++nth_trial; 338 ++nth_trial; 334 } 339 } 335 } 340 } 336 341 337 template <class T> 342 template <class T> 338 std::vector<G4double> G4SimplexDownhill<T>::Ge 343 std::vector<G4double> G4SimplexDownhill<T>::GetMinimumPoint() 339 { 344 { 340 if(!minimized) << 345 if(minimized != true) 341 { 346 { 342 GetMinimum(); 347 GetMinimum(); 343 } 348 } 344 349 345 auto it_minh = std::min_element(currentHeigh << 350 std::vector<G4double>::const_iterator it_minh = >> 351 std::min_element(currentHeights.cbegin(), currentHeights.cend()); 346 352 347 G4int imin = 0; 353 G4int imin = 0; 348 G4int i = 0; 354 G4int i = 0; 349 for(auto it = currentHeights.cbegin(); it != 355 for(auto it = currentHeights.cbegin(); it != currentHeights.cend(); ++it) 350 { 356 { 351 if(it == it_minh) 357 if(it == it_minh) 352 { 358 { 353 imin = i; 359 imin = i; 354 } 360 } 355 ++i; 361 ++i; 356 } 362 } 357 minimumPoint = currentSimplex[imin]; 363 minimumPoint = currentSimplex[imin]; 358 364 359 return minimumPoint; 365 return minimumPoint; 360 } 366 } 361 367