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