source: trunk/source/global/HEPNumerics/include/G4SimplexDownhill.icc @ 1350

Last change on this file since 1350 was 1337, checked in by garnier, 14 years ago

tag geant4.9.4 beta 1 + modifs locales
File size: 11.5 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26//
27// $Id: G4SimplexDownhill.icc,v 1.2 2007/05/11 13:05:53 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-01 $
29//
30// Author: Tatsumi Koi (SLAC/SCCS), 2007
31// --------------------------------------------------------------------------
32
33#include <iostream>
34#include <numeric>
35#include <cfloat>
36 
37template<class T> void G4SimplexDownhill<T>::init()
38{
39   alpha = 2.0;  // refrection coefficient:  0 < alpha
40   beta = 0.5;   // contraction coefficient:   0 < beta < 1
41   gamma = 2.0;  // expantion coefficient:  1 < gamma
42
43   maximum_no_trial = 10000;
44   max_se = FLT_MIN;
45   //max_ratio = FLT_EPSILON/1;
46   max_ratio = DBL_EPSILON/1;
47   minimized = false;
48}
49
50
51/*
52
53void G4SimplexDownhill<class T>::
54SetFunction( G4int n , G4double( *afunc )( std::vector < G4double > ) )
55{
56   numberOfVariable = n;
57   theFunction = afunc;
58   minimized = false;
59}
60
61*/
62
63
64template<class T>
65G4double G4SimplexDownhill<T>::GetMinimum()
66{
67
68   initialize();
69
70// First Tryal;
71 
72   //G4cout << "Begin First Trials" << G4endl;
73   doDownhill();
74   //G4cout << "End First Trials" << G4endl;
75
76   std::vector< G4double >::iterator it_minh =
77      std::min_element( currentHeights.begin() , currentHeights.end() );
78   G4int imin = -1;
79   G4int i = 0;
80   for ( std::vector< G4double >::iterator it = currentHeights.begin();
81         it != currentHeights.end(); it++ )
82   {
83      if ( it == it_minh )
84      {
85         imin = i;
86      }
87      i++;
88   }
89   std::vector< G4double > minimumPoint = currentSimplex[ imin ];
90
91// Second Trial
92
93   //std::vector< G4double > minimumPoint = currentSimplex[ 0 ];
94   initialize();
95
96   currentSimplex[ numberOfVariable ] = minimumPoint;
97
98   //G4cout << "Begin Second Trials" << G4endl;
99   doDownhill();
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;
106
107   minimized = true;
108
109   return minimum;
110
111}
112
113
114
115template<class T>
116void G4SimplexDownhill<T>::initialize()
117{
118
119   currentSimplex.resize( numberOfVariable+1 );
120   currentHeights.resize( numberOfVariable+1 );
121
122   for ( G4int i = 0 ; i < numberOfVariable ; i++ )
123   {
124      std::vector< G4double > avec ( numberOfVariable , 0.0 );
125      avec[ i ] = 1.0;
126      currentSimplex[ i ] = avec;
127   }
128
129   //std::vector< G4double > avec ( numberOfVariable , 0.0 );
130   std::vector< G4double > avec ( numberOfVariable , 1 );
131   currentSimplex[ numberOfVariable ] = avec;
132
133}
134
135
136
137template<class T>
138void G4SimplexDownhill<T>::calHeights()
139{
140
141   for ( G4int i = 0 ; i <= numberOfVariable ; i++ )
142   {
143      currentHeights[i] = getValue ( currentSimplex[i] );
144   }
145
146}
147
148
149
150template<class T>
151std::vector< G4double > G4SimplexDownhill<T>::calCentroid( G4int ih )
152{
153
154    std::vector< G4double > centroid ( numberOfVariable , 0.0 );
155
156    G4int i = 0;
157    for ( std::vector< std::vector< G4double > >::iterator
158          it = currentSimplex.begin(); it != currentSimplex.end() ; it++ )
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    }
169
170    return centroid;
171}
172
173
174 
175template<class T>
176std::vector< G4double > G4SimplexDownhill<T>::
177getReflectionPoint( std::vector< G4double > p ,
178                    std::vector< G4double > centroid )
179{
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}
191
192
193
194template<class T>
195std::vector< G4double > G4SimplexDownhill<T>::
196getExpansionPoint( 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;
209}
210
211template<class T>
212std::vector< G4double > G4SimplexDownhill<T>::
213getContractionPoint( std::vector< G4double > p ,
214                     std::vector< G4double > centroid )
215{
216   //G4cout << "Contraction" << G4endl;
217
218   std::vector< G4double > contractionP ( numberOfVariable , 0.0 );
219
220   for ( G4int i = 0 ; i < numberOfVariable ; i++ )
221   {
222      contractionP[i] = ( 1 - beta ) * centroid[i] + beta * p[i];
223   }
224   
225   return contractionP;
226}
227
228
229
230template<class T>
231G4bool G4SimplexDownhill<T>::isItGoodEnough()
232{
233   G4bool result = false;
234
235   G4double sum = std::accumulate( currentHeights.begin() ,
236                                   currentHeights.end() , 0.0 );
237   G4double average = sum/(numberOfVariable+1);
238   //G4cout << "average " << average << G4endl;
239
240   G4double delta = 0.0;
241   for ( G4int i = 0 ; i <= numberOfVariable ; i++ )
242   {
243      delta += std::abs ( currentHeights[ i ] - average ); 
244   }
245   //G4cout << "ratio of delta to average is "
246   //       << delta / (numberOfVariable+1) / average << G4endl;
247
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;
271}
272
273
274 
275template<class T>
276void G4SimplexDownhill<T>::doDownhill()
277{
278
279   G4int nth_trial = 0;
280
281   while ( nth_trial < maximum_no_trial )
282   {
283
284/*
285      G4cout << "Begining " << nth_trial << "th trial " << G4endl;
286      for ( G4int j = 0 ; j <= numberOfVariable ; j++ )
287      {
288         G4cout << "SimplexPoint " << j << ": ";
289         for ( G4int i = 0 ; i < numberOfVariable ; i++ )
290         {
291             G4cout << currentSimplex[j][i]
292                       << " ";
293         }
294         G4cout << G4endl;
295      }
296*/
297
298      calHeights();
299
300      if ( isItGoodEnough() )
301      {
302         break;
303      }
304
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++ )
319      {
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++;
334      }
335
336      //G4cout << "max " << h_H << " " << ih << G4endl;
337      //G4cout << "max-dash " << h_H2 << G4endl;
338      //G4cout << "min " << h_L << " " << il << G4endl;
339
340      std::vector< G4double > centroidPoint = calCentroid ( ih );
341
342      // REFLECTION
343      std::vector< G4double > reflectionPoint =
344        getReflectionPoint( currentSimplex[ ih ] , centroidPoint );
345
346      G4double h = getValue( reflectionPoint );
347
348      if ( h <= h_L )
349      {
350      // EXPANSION
351         std::vector< G4double > expansionPoint =
352           getExpansionPoint( reflectionPoint , centroidPoint );
353         G4double hh = getValue( expansionPoint );
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
369      {
370         if ( h <= h_H2 )
371         {
372         // Replace
373            currentSimplex[ ih ] = reflectionPoint;
374            //G4cout << "B2" << G4endl;
375         }
376         else
377         {
378            if ( h <= h_H )
379            {
380            // Replace
381              currentSimplex[ ih ] = reflectionPoint;
382              //G4cout << "BC" << G4endl;
383            }
384            // CONTRACTION
385            std::vector< G4double > contractionPoint =
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
411      }
412
413      nth_trial++;
414   }
415}
416
417
418 
419template<class T>
420std::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;
443}
Note: See TracBrowser for help on using the repository browser.