G4SimplexDownhill.icc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 // Author: Tatsumi Koi (SLAC/SCCS), 2007
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;  // refrection coefficient:  0 < alpha
00039    beta = 0.5;   // contraction coefficient:   0 < beta < 1 
00040    gamma = 2.0;  // expantion coefficient:  1 < gamma
00041 
00042    maximum_no_trial = 10000;
00043    max_se = FLT_MIN;
00044    //max_ratio = FLT_EPSILON/1;
00045    max_ratio = DBL_EPSILON/1;
00046    minimized = false;
00047 }
00048 
00049 
00050 /*
00051 
00052 void G4SimplexDownhill<class T>::
00053 SetFunction( G4int n , G4double( *afunc )( std::vector < G4double > ) ) 
00054 {
00055    numberOfVariable = n; 
00056    theFunction = afunc;
00057    minimized = false; 
00058 }
00059 
00060 */
00061 
00062 
00063 template<class T> 
00064 G4double G4SimplexDownhill<T>::GetMinimum()
00065 {
00066 
00067    initialize();
00068 
00069 // First Tryal;
00070   
00071    //G4cout << "Begin First Trials" << G4endl;
00072    doDownhill();
00073    //G4cout << "End First Trials" << G4endl;
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 // Second Trial
00091 
00092    //std::vector< G4double > minimumPoint = currentSimplex[ 0 ];
00093    initialize();
00094 
00095    currentSimplex[ numberOfVariable ] = minimumPoint;
00096 
00097    //G4cout << "Begin Second Trials" << G4endl;
00098    doDownhill();
00099    //G4cout << "End Second Trials" << G4endl;
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    //std::vector< G4double > avec ( numberOfVariable , 0.0 ); 
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    //G4cout << "Reflection" << G4endl;
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    //G4cout << "Expantion" << G4endl;
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    //G4cout << "Contraction" << G4endl;
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    //G4cout << "average " << average << G4endl;
00238 
00239    G4double delta = 0.0; 
00240    for ( G4int i = 0 ; i <= numberOfVariable ; i++ )
00241    {
00242       delta += std::abs ( currentHeights[ i ] - average );  
00243    }
00244    //G4cout << "ratio of delta to average is "
00245    //       << delta / (numberOfVariable+1) / average << G4endl;
00246 
00247    if ( delta/(numberOfVariable+1)/average < max_ratio )
00248    {
00249       result = true; 
00250    }
00251   
00252 /*
00253    G4double sigma = 0.0; 
00254    G4cout << "average " << average << G4endl;
00255    for ( G4int i = 0 ; i <= numberOfVariable ; i++ )
00256    {
00257       sigma += ( currentHeights[ i ] - average )
00258               *( currentHeights[ i ] - average );  
00259    }
00260 
00261    G4cout << "standard error of hs "
00262           << std::sqrt ( sigma ) / (numberOfVariable+1) << G4endl;
00263    if ( std::sqrt ( sigma ) / (numberOfVariable+1) < max_se )
00264    {
00265       result = true; 
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       G4cout << "Begining " << nth_trial << "th trial " << G4endl;
00285       for ( G4int j = 0 ; j <= numberOfVariable ; j++ )
00286       {
00287          G4cout << "SimplexPoint " << j << ": "; 
00288          for ( G4int i = 0 ; i < numberOfVariable ; i++ )
00289          {
00290              G4cout << currentSimplex[j][i] 
00291                        << " ";
00292          }
00293          G4cout << G4endl; 
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       //G4cout << "max " << h_H << " " << ih << G4endl;
00336       //G4cout << "max-dash " << h_H2 << G4endl;
00337       //G4cout << "min " << h_L << " " << il << G4endl;
00338 
00339       std::vector< G4double > centroidPoint = calCentroid ( ih );
00340 
00341       // REFLECTION
00342       std::vector< G4double > reflectionPoint =
00343         getReflectionPoint( currentSimplex[ ih ] , centroidPoint );
00344 
00345       G4double h = getValue( reflectionPoint );
00346 
00347       if ( h <= h_L ) 
00348       {
00349       // EXPANSION
00350          std::vector< G4double > expansionPoint =
00351            getExpansionPoint( reflectionPoint , centroidPoint );
00352          G4double hh = getValue( expansionPoint );
00353           
00354          if ( hh <= h_L )
00355          {
00356          // Replace 
00357             currentSimplex[ ih ] = expansionPoint;
00358             //G4cout << "A" << G4endl;
00359          } 
00360          else 
00361          {
00362          // Replace 
00363             currentSimplex[ ih ] = reflectionPoint;
00364             //G4cout << "B1" << G4endl;
00365          }
00366       } 
00367       else 
00368       {
00369          if ( h <= h_H2 )
00370          {
00371          // Replace 
00372             currentSimplex[ ih ] = reflectionPoint;
00373             //G4cout << "B2" << G4endl;
00374          }
00375          else
00376          {
00377             if ( h <= h_H )
00378             {
00379             // Replace 
00380               currentSimplex[ ih ] = reflectionPoint;
00381               //G4cout << "BC" << G4endl;
00382             }
00383             // CONTRACTION
00384             std::vector< G4double > contractionPoint =
00385               getContractionPoint( currentSimplex[ ih ] , centroidPoint );
00386             G4double hh = getValue( contractionPoint );
00387             if ( hh <= h_H )
00388             {
00389             // Replace 
00390                currentSimplex[ ih ] = contractionPoint;
00391                //G4cout << "C" << G4endl;
00392             }
00393             else
00394             {
00395             // Replace
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                //G4cout << "D" << G4endl;
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 }

Generated on Mon May 27 17:49:50 2013 for Geant4 by  doxygen 1.4.7