00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include <iostream>
00033 #include <numeric>
00034 #include <cfloat>
00035
00036 template<class T> void G4SimplexDownhill<T>::init()
00037 {
00038 alpha = 2.0;
00039 beta = 0.5;
00040 gamma = 2.0;
00041
00042 maximum_no_trial = 10000;
00043 max_se = FLT_MIN;
00044
00045 max_ratio = DBL_EPSILON/1;
00046 minimized = false;
00047 }
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063 template<class T>
00064 G4double G4SimplexDownhill<T>::GetMinimum()
00065 {
00066
00067 initialize();
00068
00069
00070
00071
00072 doDownhill();
00073
00074
00075 std::vector< G4double >::iterator it_minh =
00076 std::min_element( currentHeights.begin() , currentHeights.end() );
00077 G4int imin = -1;
00078 G4int i = 0;
00079 for ( std::vector< G4double >::iterator it = currentHeights.begin();
00080 it != currentHeights.end(); it++ )
00081 {
00082 if ( it == it_minh )
00083 {
00084 imin = i;
00085 }
00086 i++;
00087 }
00088 minimumPoint = currentSimplex[ imin ];
00089
00090
00091
00092
00093 initialize();
00094
00095 currentSimplex[ numberOfVariable ] = minimumPoint;
00096
00097
00098 doDownhill();
00099
00100
00101 G4double sum = std::accumulate( currentHeights.begin() ,
00102 currentHeights.end() , 0.0 );
00103 G4double average = sum/(numberOfVariable+1);
00104 G4double minimum = average;
00105
00106 minimized = true;
00107
00108 return minimum;
00109
00110 }
00111
00112
00113
00114 template<class T>
00115 void G4SimplexDownhill<T>::initialize()
00116 {
00117
00118 currentSimplex.resize( numberOfVariable+1 );
00119 currentHeights.resize( numberOfVariable+1 );
00120
00121 for ( G4int i = 0 ; i < numberOfVariable ; i++ )
00122 {
00123 std::vector< G4double > avec ( numberOfVariable , 0.0 );
00124 avec[ i ] = 1.0;
00125 currentSimplex[ i ] = avec;
00126 }
00127
00128
00129 std::vector< G4double > avec ( numberOfVariable , 1 );
00130 currentSimplex[ numberOfVariable ] = avec;
00131
00132 }
00133
00134
00135
00136 template<class T>
00137 void G4SimplexDownhill<T>::calHeights()
00138 {
00139
00140 for ( G4int i = 0 ; i <= numberOfVariable ; i++ )
00141 {
00142 currentHeights[i] = getValue ( currentSimplex[i] );
00143 }
00144
00145 }
00146
00147
00148
00149 template<class T>
00150 std::vector< G4double > G4SimplexDownhill<T>::calCentroid( G4int ih )
00151 {
00152
00153 std::vector< G4double > centroid ( numberOfVariable , 0.0 );
00154
00155 G4int i = 0;
00156 for ( std::vector< std::vector< G4double > >::iterator
00157 it = currentSimplex.begin(); it != currentSimplex.end() ; it++ )
00158 {
00159 if ( i != ih )
00160 {
00161 for ( G4int j = 0 ; j < numberOfVariable ; j++ )
00162 {
00163 centroid[j] += (*it)[j]/numberOfVariable;
00164 }
00165 }
00166 i++;
00167 }
00168
00169 return centroid;
00170 }
00171
00172
00173
00174 template<class T>
00175 std::vector< G4double > G4SimplexDownhill<T>::
00176 getReflectionPoint( std::vector< G4double > p ,
00177 std::vector< G4double > centroid )
00178 {
00179
00180
00181 std::vector< G4double > reflectionP ( numberOfVariable , 0.0 );
00182
00183 for ( G4int i = 0 ; i < numberOfVariable ; i++ )
00184 {
00185 reflectionP[ i ] = ( 1 + alpha ) * centroid[ i ] - alpha * p[ i ];
00186 }
00187
00188 return reflectionP;
00189 }
00190
00191
00192
00193 template<class T>
00194 std::vector< G4double > G4SimplexDownhill<T>::
00195 getExpansionPoint( std::vector< G4double > p ,
00196 std::vector< G4double > centroid )
00197 {
00198
00199
00200 std::vector< G4double > expansionP ( numberOfVariable , 0.0 );
00201
00202 for ( G4int i = 0 ; i < numberOfVariable ; i++ )
00203 {
00204 expansionP[i] = ( 1 - gamma ) * centroid[i] + gamma * p[i];
00205 }
00206
00207 return expansionP;
00208 }
00209
00210 template<class T>
00211 std::vector< G4double > G4SimplexDownhill<T>::
00212 getContractionPoint( std::vector< G4double > p ,
00213 std::vector< G4double > centroid )
00214 {
00215
00216
00217 std::vector< G4double > contractionP ( numberOfVariable , 0.0 );
00218
00219 for ( G4int i = 0 ; i < numberOfVariable ; i++ )
00220 {
00221 contractionP[i] = ( 1 - beta ) * centroid[i] + beta * p[i];
00222 }
00223
00224 return contractionP;
00225 }
00226
00227
00228
00229 template<class T>
00230 G4bool G4SimplexDownhill<T>::isItGoodEnough()
00231 {
00232 G4bool result = false;
00233
00234 G4double sum = std::accumulate( currentHeights.begin() ,
00235 currentHeights.end() , 0.0 );
00236 G4double average = sum/(numberOfVariable+1);
00237
00238
00239 G4double delta = 0.0;
00240 for ( G4int i = 0 ; i <= numberOfVariable ; i++ )
00241 {
00242 delta += std::abs ( currentHeights[ i ] - average );
00243 }
00244
00245
00246
00247 if ( delta/(numberOfVariable+1)/average < max_ratio )
00248 {
00249 result = true;
00250 }
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269 return result;
00270 }
00271
00272
00273
00274 template<class T>
00275 void G4SimplexDownhill<T>::doDownhill()
00276 {
00277
00278 G4int nth_trial = 0;
00279
00280 while ( nth_trial < maximum_no_trial )
00281 {
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297 calHeights();
00298
00299 if ( isItGoodEnough() )
00300 {
00301 break;
00302 }
00303
00304 std::vector< G4double >::iterator it_maxh =
00305 std::max_element( currentHeights.begin() , currentHeights.end() );
00306 std::vector< G4double >::iterator it_minh =
00307 std::min_element( currentHeights.begin() , currentHeights.end() );;
00308
00309 G4double h_H = *it_maxh;
00310 G4double h_L = *it_minh;
00311
00312 G4int ih = 0;;
00313 G4int il = 0;
00314 G4double h_H2 =0.0;
00315 G4int i = 0;
00316 for ( std::vector< G4double >::iterator
00317 it = currentHeights.begin(); it != currentHeights.end(); it++ )
00318 {
00319 if ( it == it_maxh )
00320 {
00321 ih = i;
00322 }
00323 else
00324 {
00325 h_H2 = std::max( h_H2 , *it );
00326 }
00327
00328 if ( it == it_minh )
00329 {
00330 il = i;
00331 }
00332 i++;
00333 }
00334
00335
00336
00337
00338
00339 std::vector< G4double > centroidPoint = calCentroid ( ih );
00340
00341
00342 std::vector< G4double > reflectionPoint =
00343 getReflectionPoint( currentSimplex[ ih ] , centroidPoint );
00344
00345 G4double h = getValue( reflectionPoint );
00346
00347 if ( h <= h_L )
00348 {
00349
00350 std::vector< G4double > expansionPoint =
00351 getExpansionPoint( reflectionPoint , centroidPoint );
00352 G4double hh = getValue( expansionPoint );
00353
00354 if ( hh <= h_L )
00355 {
00356
00357 currentSimplex[ ih ] = expansionPoint;
00358
00359 }
00360 else
00361 {
00362
00363 currentSimplex[ ih ] = reflectionPoint;
00364
00365 }
00366 }
00367 else
00368 {
00369 if ( h <= h_H2 )
00370 {
00371
00372 currentSimplex[ ih ] = reflectionPoint;
00373
00374 }
00375 else
00376 {
00377 if ( h <= h_H )
00378 {
00379
00380 currentSimplex[ ih ] = reflectionPoint;
00381
00382 }
00383
00384 std::vector< G4double > contractionPoint =
00385 getContractionPoint( currentSimplex[ ih ] , centroidPoint );
00386 G4double hh = getValue( contractionPoint );
00387 if ( hh <= h_H )
00388 {
00389
00390 currentSimplex[ ih ] = contractionPoint;
00391
00392 }
00393 else
00394 {
00395
00396 for ( G4int j = 0 ; j <= numberOfVariable ; j++ )
00397 {
00398 std::vector< G4double > vec ( numberOfVariable , 0.0 );
00399 for ( G4int k = 0 ; k < numberOfVariable ; k++ )
00400 {
00401 vec[ k ] = ( currentSimplex[ j ][ k ]
00402 + currentSimplex[ il ][ k ] ) / 2.0;
00403 }
00404 currentSimplex[ j ] = vec;
00405 }
00406
00407 }
00408 }
00409
00410 }
00411
00412 nth_trial++;
00413 }
00414 }
00415
00416
00417
00418 template<class T>
00419 std::vector< G4double > G4SimplexDownhill<T>::GetMinimumPoint()
00420 {
00421 if ( minimized != true )
00422 {
00423 GetMinimum();
00424 }
00425
00426 std::vector< G4double >::iterator it_minh =
00427 std::min_element( currentHeights.begin() , currentHeights.end() );;
00428 G4int imin = -1;
00429 G4int i = 0;
00430 for ( std::vector< G4double >::iterator
00431 it = currentHeights.begin(); it != currentHeights.end(); it++ )
00432 {
00433 if ( it == it_minh )
00434 {
00435 imin = i;
00436 }
00437 i++;
00438 }
00439 minimumPoint = currentSimplex[ imin ];
00440
00441 return minimumPoint;
00442 }