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

Last change on this file since 1036 was 850, checked in by garnier, 17 years ago

geant4.8.2 beta
File size: 11.4 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: HEAD $
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.