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 9.0.p2)


  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: geant4-09-00-patch-01 $
 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    std::vector< G4double > 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    std::vector< G4double > minimumPoint = currentSimplex[ imin ];
                                                   >> 441 
                                                   >> 442    return minimumPoint; 
360 }                                                 443 }
361                                                   444