Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // G4SimplexDownhill inline methods implementation 27 // 28 // Author: Tatsumi Koi (SLAC/SCCS), 2007 29 // -------------------------------------------------------------------------- 30 31 #include <cfloat> 32 #include <iostream> 33 #include <numeric> 34 35 template <class T> 36 void G4SimplexDownhill<T>::init() 37 { 38 alpha = 2.0; // refrection coefficient: 0 < alpha 39 beta = 0.5; // contraction coefficient: 0 < beta < 1 40 gamma = 2.0; // expantion coefficient: 1 < gamma 41 42 maximum_no_trial = 10000; 43 max_se = FLT_MIN; 44 // max_ratio = FLT_EPSILON/1; 45 max_ratio = DBL_EPSILON / 1; 46 minimized = false; 47 } 48 49 /* 50 51 void G4SimplexDownhill<class T>:: 52 SetFunction( G4int n , G4double( *afunc )( std::vector < G4double > ) ) 53 { 54 numberOfVariable = n; 55 theFunction = afunc; 56 minimized = false; 57 } 58 59 */ 60 61 template <class T> 62 G4double G4SimplexDownhill<T>::GetMinimum() 63 { 64 initialize(); 65 66 // First Tryal; 67 68 // G4cout << "Begin First Trials" << G4endl; 69 doDownhill(); 70 // G4cout << "End First Trials" << G4endl; 71 72 auto it_minh = std::min_element(currentHeights.cbegin(), currentHeights.cend()); 73 G4int imin = 0; 74 G4int i = 0; 75 for(auto it = currentHeights.cbegin(); it != currentHeights.cend(); ++it) 76 { 77 if(it == it_minh) 78 { 79 imin = i; 80 } 81 ++i; 82 } 83 minimumPoint = currentSimplex[imin]; 84 85 // Second Trial 86 87 // std::vector< G4double > minimumPoint = currentSimplex[ 0 ]; 88 initialize(); 89 90 currentSimplex[numberOfVariable] = minimumPoint; 91 92 // G4cout << "Begin Second Trials" << G4endl; 93 doDownhill(); 94 // G4cout << "End Second Trials" << G4endl; 95 96 G4double sum = 97 std::accumulate(currentHeights.begin(), currentHeights.end(), 0.0); 98 G4double average = sum / (numberOfVariable + 1); 99 G4double minimum = average; 100 101 minimized = true; 102 103 return minimum; 104 } 105 106 template <class T> 107 void G4SimplexDownhill<T>::initialize() 108 { 109 currentSimplex.resize(numberOfVariable + 1); 110 currentHeights.resize(numberOfVariable + 1); 111 112 for(G4int i = 0; i < numberOfVariable; ++i) 113 { 114 std::vector<G4double> avec(numberOfVariable, 0.0); 115 avec[i] = 1.0; 116 currentSimplex[i] = std::move(avec); 117 } 118 119 // std::vector< G4double > avec ( numberOfVariable , 0.0 ); 120 std::vector<G4double> avec(numberOfVariable, 1); 121 currentSimplex[numberOfVariable] = std::move(avec); 122 } 123 124 template <class T> 125 void G4SimplexDownhill<T>::calHeights() 126 { 127 for(G4int i = 0; i <= numberOfVariable; ++i) 128 { 129 currentHeights[i] = getValue(currentSimplex[i]); 130 } 131 } 132 133 template <class T> 134 std::vector<G4double> G4SimplexDownhill<T>::calCentroid(G4int ih) 135 { 136 std::vector<G4double> centroid(numberOfVariable, 0.0); 137 138 G4int i = 0; 139 for(const auto & it : currentSimplex) 140 { 141 if(i != ih) 142 { 143 for(G4int j = 0; j < numberOfVariable; ++j) 144 { 145 centroid[j] += it[j] / numberOfVariable; 146 } 147 } 148 ++i; 149 } 150 151 return centroid; 152 } 153 154 template <class T> 155 std::vector<G4double> G4SimplexDownhill<T>::getReflectionPoint( 156 std::vector<G4double> p, std::vector<G4double> centroid) 157 { 158 // G4cout << "Reflection" << G4endl; 159 160 std::vector<G4double> reflectionP(numberOfVariable, 0.0); 161 162 for(G4int i = 0; i < numberOfVariable; ++i) 163 { 164 reflectionP[i] = (1 + alpha) * centroid[i] - alpha * p[i]; 165 } 166 167 return reflectionP; 168 } 169 170 template <class T> 171 std::vector<G4double> G4SimplexDownhill<T>::getExpansionPoint( 172 std::vector<G4double> p, std::vector<G4double> centroid) 173 { 174 // G4cout << "Expantion" << G4endl; 175 176 std::vector<G4double> expansionP(numberOfVariable, 0.0); 177 178 for(G4int i = 0; i < numberOfVariable; ++i) 179 { 180 expansionP[i] = (1 - gamma) * centroid[i] + gamma * p[i]; 181 } 182 183 return expansionP; 184 } 185 186 template <class T> 187 std::vector<G4double> G4SimplexDownhill<T>::getContractionPoint( 188 std::vector<G4double> p, std::vector<G4double> centroid) 189 { 190 std::vector<G4double> contractionP(numberOfVariable, 0.0); 191 192 for(G4int i = 0; i < numberOfVariable; ++i) 193 { 194 contractionP[i] = (1 - beta) * centroid[i] + beta * p[i]; 195 } 196 197 return contractionP; 198 } 199 200 template <class T> 201 G4bool G4SimplexDownhill<T>::isItGoodEnough() 202 { 203 G4double sum = 204 std::accumulate(currentHeights.begin(), currentHeights.end(), 0.0); 205 G4double average = sum / (numberOfVariable + 1); 206 207 G4double delta = 0.0; 208 for(G4int i = 0; i <= numberOfVariable; ++i) 209 { 210 delta += std::abs(currentHeights[i] - average); 211 } 212 213 G4bool result = false; 214 if (average > 0.0) 215 { 216 result = ((delta / (numberOfVariable + 1) / average) < max_ratio); 217 } 218 return result; 219 } 220 221 template <class T> 222 void G4SimplexDownhill<T>::doDownhill() 223 { 224 G4int nth_trial = 0; 225 226 while(nth_trial < maximum_no_trial) 227 { 228 calHeights(); 229 230 if(isItGoodEnough()) 231 { 232 break; 233 } 234 235 auto it_maxh = std::max_element(currentHeights.cbegin(), currentHeights.cend()); 236 auto it_minh = std::min_element(currentHeights.cbegin(), currentHeights.cend()); 237 238 G4double h_H = *it_maxh; 239 G4double h_L = *it_minh; 240 241 G4int ih = 0; 242 G4int il = 0; 243 G4double h_H2 = 0.0; 244 G4int i = 0; 245 for(auto it = currentHeights.cbegin(); it != currentHeights.cend(); ++it) 246 { 247 if(it == it_maxh) 248 { 249 ih = i; 250 } 251 else 252 { 253 h_H2 = std::max(h_H2, *it); 254 } 255 256 if(it == it_minh) 257 { 258 il = i; 259 } 260 ++i; 261 } 262 263 std::vector<G4double> centroidPoint = calCentroid(ih); 264 265 // REFLECTION 266 std::vector<G4double> reflectionPoint = 267 getReflectionPoint(currentSimplex[ih], centroidPoint); 268 269 G4double h = getValue(reflectionPoint); 270 271 if(h <= h_L) 272 { 273 // EXPANSION 274 std::vector<G4double> expansionPoint = 275 getExpansionPoint(reflectionPoint, std::move(centroidPoint)); 276 G4double hh = getValue(expansionPoint); 277 278 if(hh <= h_L) 279 { 280 // Replace 281 currentSimplex[ih] = std::move(expansionPoint); 282 // G4cout << "A" << G4endl; 283 } 284 else 285 { 286 // Replace 287 currentSimplex[ih] = std::move(reflectionPoint); 288 // G4cout << "B1" << G4endl; 289 } 290 } 291 else 292 { 293 if(h <= h_H2) 294 { 295 // Replace 296 currentSimplex[ih] = std::move(reflectionPoint); 297 // G4cout << "B2" << G4endl; 298 } 299 else 300 { 301 if(h <= h_H) 302 { 303 // Replace 304 currentSimplex[ih] = std::move(reflectionPoint); 305 // G4cout << "BC" << G4endl; 306 } 307 // CONTRACTION 308 std::vector<G4double> contractionPoint = 309 getContractionPoint(currentSimplex[ih], std::move(centroidPoint)); 310 G4double hh = getValue(contractionPoint); 311 if(hh <= h_H) 312 { 313 // Replace 314 currentSimplex[ih] = std::move(contractionPoint); 315 // G4cout << "C" << G4endl; 316 } 317 else 318 { 319 // Replace 320 for(G4int j = 0; j <= numberOfVariable; ++j) 321 { 322 std::vector<G4double> vec(numberOfVariable, 0.0); 323 for(G4int k = 0; k < numberOfVariable; ++k) 324 { 325 vec[k] = (currentSimplex[j][k] + currentSimplex[il][k]) / 2.0; 326 } 327 currentSimplex[j] = std::move(vec); 328 } 329 } 330 } 331 } 332 333 ++nth_trial; 334 } 335 } 336 337 template <class T> 338 std::vector<G4double> G4SimplexDownhill<T>::GetMinimumPoint() 339 { 340 if(!minimized) 341 { 342 GetMinimum(); 343 } 344 345 auto it_minh = std::min_element(currentHeights.cbegin(), currentHeights.cend()); 346 347 G4int imin = 0; 348 G4int i = 0; 349 for(auto it = currentHeights.cbegin(); it != currentHeights.cend(); ++it) 350 { 351 if(it == it_minh) 352 { 353 imin = i; 354 } 355 ++i; 356 } 357 minimumPoint = currentSimplex[imin]; 358 359 return minimumPoint; 360 } 361