EIC code displayed by LXR master/​include/​Geant4/​G4SimplexDownhill.icc	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2025-09-15 08:59:24

 //
 // ********************************************************************
 // * License and Disclaimer                                           *
 // *                                                                  *
 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
 // * conditions of the Geant4 Software License,  included in the file *
 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
 // * include a list of copyright holders.                             *
 // *                                                                  *
 // * Neither the authors of this software system, nor their employing *
 // * institutes,nor the agencies providing financial support for this *
 // * work  make  any representation or  warranty, express or implied, *
 // * regarding  this  software system or assume any liability for its *
 // * use.  Please see the license in the file  LICENSE  and URL above *
 // * for the full disclaimer and the limitation of liability.         *
 // *                                                                  *
 // * This  code  implementation is the result of  the  scientific and *
 // * technical work of the GEANT4 collaboration.                      *
 // * By using,  copying,  modifying or  distributing the software (or *
 // * any work based  on the software)  you  agree  to acknowledge its *
 // * use  in  resulting  scientific  publications,  and indicate your *
 // * acceptance of all terms of the Geant4 Software license.          *
 // ********************************************************************
 //
 // G4SimplexDownhill inline methods implementation
 //
 // Author: Tatsumi Koi (SLAC/SCCS), 2007
 // --------------------------------------------------------------------------
 
 #include <cfloat>
 #include <iostream>
 #include <numeric>
 
 template <class T>
 void G4SimplexDownhill<T>::init()
 {
   alpha = 2.0;  // refrection coefficient:  0 < alpha
   beta  = 0.5;  // contraction coefficient:   0 < beta < 1
   gamma = 2.0;  // expantion coefficient:  1 < gamma
 
   maximum_no_trial = 10000;
   max_se           = FLT_MIN;
   // max_ratio = FLT_EPSILON/1;
   max_ratio = DBL_EPSILON / 1;
   minimized = false;
 }
 
 /*
 
 void G4SimplexDownhill<class T>::
 SetFunction( G4int n , G4double( *afunc )( std::vector < G4double > ) )
 {
    numberOfVariable = n;
    theFunction = afunc;
    minimized = false;
 }
 
 */
 
 template <class T>
 G4double G4SimplexDownhill<T>::GetMinimum()
 {
   initialize();
 
   // First Tryal;
 
   // G4cout << "Begin First Trials" << G4endl;
   doDownhill();
   // G4cout << "End First Trials" << G4endl;
 
   auto it_minh = std::min_element(currentHeights.cbegin(), currentHeights.cend());
   G4int imin = 0;
   G4int i    = 0;
   for(auto it = currentHeights.cbegin(); it != currentHeights.cend(); ++it)
   {
     if(it == it_minh)
     {
       imin = i;
     }
     ++i;
   }
   minimumPoint = currentSimplex[imin];
 
   // Second Trial
 
   // std::vector< G4double > minimumPoint = currentSimplex[ 0 ];
   initialize();
 
   currentSimplex[numberOfVariable] = minimumPoint;
 
   // G4cout << "Begin Second Trials" << G4endl;
   doDownhill();
   // G4cout << "End Second Trials" << G4endl;
 
   G4double sum =
     std::accumulate(currentHeights.begin(), currentHeights.end(), 0.0);
   G4double average = sum / (numberOfVariable + 1);
   G4double minimum = average;
 
   minimized = true;
 
   return minimum;
 }
 
 template <class T>
 void G4SimplexDownhill<T>::initialize()
 {
   currentSimplex.resize(numberOfVariable + 1);
   currentHeights.resize(numberOfVariable + 1);
 
   for(G4int i = 0; i < numberOfVariable; ++i)
   {
     std::vector<G4double> avec(numberOfVariable, 0.0);
     avec[i]           = 1.0;
     currentSimplex[i] = std::move(avec);
   }
 
   // std::vector< G4double > avec ( numberOfVariable , 0.0 );
   std::vector<G4double> avec(numberOfVariable, 1);
   currentSimplex[numberOfVariable] = std::move(avec);
 }
 
 template <class T>
 void G4SimplexDownhill<T>::calHeights()
 {
   for(G4int i = 0; i <= numberOfVariable; ++i)
   {
     currentHeights[i] = getValue(currentSimplex[i]);
   }
 }
 
 template <class T>
 std::vector<G4double> G4SimplexDownhill<T>::calCentroid(G4int ih)
 {
   std::vector<G4double> centroid(numberOfVariable, 0.0);
 
   G4int i = 0;
   for(const auto & it : currentSimplex)
   {
     if(i != ih)
     {
       for(G4int j = 0; j < numberOfVariable; ++j)
       {
         centroid[j] += it[j] / numberOfVariable;
       }
     }
     ++i;
   }
 
   return centroid;
 }
 
 template <class T>
 std::vector<G4double> G4SimplexDownhill<T>::getReflectionPoint(
   std::vector<G4double> p, std::vector<G4double> centroid)
 {
   // G4cout << "Reflection" << G4endl;
 
   std::vector<G4double> reflectionP(numberOfVariable, 0.0);
 
   for(G4int i = 0; i < numberOfVariable; ++i)
   {
     reflectionP[i] = (1 + alpha) * centroid[i] - alpha * p[i];
   }
 
   return reflectionP;
 }
 
 template <class T>
 std::vector<G4double> G4SimplexDownhill<T>::getExpansionPoint(
   std::vector<G4double> p, std::vector<G4double> centroid)
 {
   // G4cout << "Expantion" << G4endl;
 
   std::vector<G4double> expansionP(numberOfVariable, 0.0);
 
   for(G4int i = 0; i < numberOfVariable; ++i)
   {
     expansionP[i] = (1 - gamma) * centroid[i] + gamma * p[i];
   }
 
   return expansionP;
 }
 
 template <class T>
 std::vector<G4double> G4SimplexDownhill<T>::getContractionPoint(
   std::vector<G4double> p, std::vector<G4double> centroid)
 {
   std::vector<G4double> contractionP(numberOfVariable, 0.0);
 
   for(G4int i = 0; i < numberOfVariable; ++i)
   {
     contractionP[i] = (1 - beta) * centroid[i] + beta * p[i];
   }
 
   return contractionP;
 }
 
 template <class T>
 G4bool G4SimplexDownhill<T>::isItGoodEnough()
 {
   G4double sum =
     std::accumulate(currentHeights.begin(), currentHeights.end(), 0.0);
   G4double average = sum / (numberOfVariable + 1);
 
   G4double delta = 0.0;
   for(G4int i = 0; i <= numberOfVariable; ++i)
   {
     delta += std::abs(currentHeights[i] - average);
   }
 
   G4bool result = false;
   if (average > 0.0)
   {
     result = ((delta / (numberOfVariable + 1) / average) < max_ratio);
   }
   return result;
 }
 
 template <class T>
 void G4SimplexDownhill<T>::doDownhill()
 {
   G4int nth_trial = 0;
 
   while(nth_trial < maximum_no_trial)
   {
     calHeights();
 
     if(isItGoodEnough())
     {
       break;
     }
 
     auto it_maxh = std::max_element(currentHeights.cbegin(), currentHeights.cend());
     auto it_minh = std::min_element(currentHeights.cbegin(), currentHeights.cend());
 
     G4double h_H = *it_maxh;
     G4double h_L = *it_minh;
 
     G4int ih      = 0;
     G4int il      = 0;
     G4double h_H2 = 0.0;
     G4int i       = 0;
     for(auto it = currentHeights.cbegin(); it != currentHeights.cend(); ++it)
     {
       if(it == it_maxh)
       {
         ih = i;
       }
       else
       {
         h_H2 = std::max(h_H2, *it);
       }
 
       if(it == it_minh)
       {
         il = i;
       }
       ++i;
     }
 
     std::vector<G4double> centroidPoint = calCentroid(ih);
 
     // REFLECTION
     std::vector<G4double> reflectionPoint =
       getReflectionPoint(currentSimplex[ih], centroidPoint);
 
     G4double h = getValue(reflectionPoint);
 
     if(h <= h_L)
     {
       // EXPANSION
       std::vector<G4double> expansionPoint =
         getExpansionPoint(reflectionPoint, std::move(centroidPoint));
       G4double hh = getValue(expansionPoint);
 
       if(hh <= h_L)
       {
         // Replace
         currentSimplex[ih] = std::move(expansionPoint);
         // G4cout << "A" << G4endl;
       }
       else
       {
         // Replace
         currentSimplex[ih] = std::move(reflectionPoint);
         // G4cout << "B1" << G4endl;
       }
     }
     else
     {
       if(h <= h_H2)
       {
         // Replace
         currentSimplex[ih] = std::move(reflectionPoint);
         // G4cout << "B2" << G4endl;
       }
       else
       {
         if(h <= h_H)
         {
           // Replace
           currentSimplex[ih] = std::move(reflectionPoint);
           // G4cout << "BC" << G4endl;
         }
         // CONTRACTION
         std::vector<G4double> contractionPoint =
           getContractionPoint(currentSimplex[ih], std::move(centroidPoint));
         G4double hh = getValue(contractionPoint);
         if(hh <= h_H)
         {
           // Replace
           currentSimplex[ih] = std::move(contractionPoint);
           // G4cout << "C" << G4endl;
         }
         else
         {
           // Replace
           for(G4int j = 0; j <= numberOfVariable; ++j)
           {
             std::vector<G4double> vec(numberOfVariable, 0.0);
             for(G4int k = 0; k < numberOfVariable; ++k)
             {
               vec[k] = (currentSimplex[j][k] + currentSimplex[il][k]) / 2.0;
             }
             currentSimplex[j] = std::move(vec);
           }
         }
       }
     }
 
     ++nth_trial;
   }
 }
 
 template <class T>
 std::vector<G4double> G4SimplexDownhill<T>::GetMinimumPoint()
 {
   if(!minimized)
   {
     GetMinimum();
   }
 
   auto it_minh = std::min_element(currentHeights.cbegin(), currentHeights.cend());
 
   G4int imin = 0;
   G4int i    = 0;
   for(auto it = currentHeights.cbegin(); it != currentHeights.cend(); ++it)
   {
     if(it == it_minh)
     {
       imin = i;
     }
     ++i;
   }
   minimumPoint = currentSimplex[imin];
 
   return minimumPoint;
 }

This page was automatically generated by the 2.3.7 LXR engine. The LXR team