Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/global/HEPNumerics/include/G4SimplexDownhill.icc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /global/HEPNumerics/include/G4SimplexDownhill.icc (Version 11.3.0) and /global/HEPNumerics/include/G4SimplexDownhill.icc (Version 10.1)


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