| # Gaudi Framework, version v25r0 |

Home |  | Generated: Mon Feb 17 2014

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
NumericalDefiniteIntegral.cpp
Go to the documentation of this file.
1 // $Id: NumericalDefiniteIntegral.cpp,v 1.4 2007/11/20 13:00:17 marcocle Exp $
2 // ============================================================================
3 // Include
4 // ============================================================================
5 // STD & STL
6 // ============================================================================
7 #include <vector>
8 #include <algorithm>
9 // ============================================================================
10 // GSL
11 // ============================================================================
12 #include "gsl/gsl_errno.h"
13 #include "gsl/gsl_integration.h"
14 // ============================================================================
15 // GaudiKernel
16 // ============================================================================
17 #include "GaudiKernel/GaudiException.h"
18 // ============================================================================
19 // GaudiMath
20 // ============================================================================
21 #include "GaudiMath/NumericalDefiniteIntegral.h"
22 #include "GaudiMath/NumericalDerivative.h"
23 // ============================================================================
24 // Local
25 // ============================================================================
26 #include "Helpers.h"
27 // ============================================================================
28
29
30 #ifdef __ICC
31 // disable icc remark #1572: floating-point equality and inequality comparisons are unreliable
32 // The comparison are meant
33 #pragma warning(disable:1572)
34 #endif
35 #ifdef WIN32
36 // Disable the warning
37 // C4996: 'std::copy': Function call with parameters that may be unsafe
38 // The parameters are checked
39 #pragma warning(disable:4996)
40 #endif
41
42 namespace Genfun
43 {
44 namespace GaudiMathImplementation
45 {
46
47 struct NumericalDefiniteIntegral::_Workspace
48 { gsl_integration_workspace* ws ; };
49 struct NumericalDefiniteIntegral::_Function
50 { gsl_function* fn ; };
51
52 // ========================================================================
54 // ========================================================================
55 FUNCTION_OBJECT_IMP( NumericalDefiniteIntegral )
56 // ========================================================================
57
58
109 NumericalDefiniteIntegral::NumericalDefiniteIntegral
110 ( const AbsFunction& function ,
111 const size_t index ,
112 const double a ,
113 const double b ,
114 const GaudiMath::Integration::Type type ,
115 const GaudiMath::Integration::KronrodRule rule ,
116 const double epsabs ,
117 const double epsrel ,
118 const size_t size )
119 : AbsFunction ()
120 , m_function ( function.clone () )
121 , m_DIM ( 0 )
122 , m_index ( index )
123 , m_a ( a )
124 , m_b ( b )
125 , m_ia ( false )
126 , m_ib ( false )
127 , m_type ( type )
128 , m_category ( GaudiMath::Integration::Finite )
129 , m_rule ( rule )
130 , m_points ( )
131 , m_pdata ( 0 )
132 , m_epsabs ( epsabs )
133 , m_epsrel ( epsrel )
134 , m_result ( GSL_NEGINF )
135 , m_error ( GSL_POSINF )
136 , m_size ( size )
137 , m_ws ()
138 , m_argument ()
139 , m_argF ()
140 {
141 if ( GaudiMath::Integration::Fixed == m_rule )
142 { m_rule = GaudiMath::Integration::Default ; }
143 if ( function.dimensionality() < 2 )
144 { Exception("::constructor: invalid dimensionality ") ; }
145 if ( m_index >= function.dimensionality() )
146 { Exception("::constructor: invalid variable index") ; }
147
148 m_DIM = function.dimensionality() - 1 ;
149 m_argument = Argument( m_DIM ) ;
150 m_argF = Argument( m_DIM + 1 ) ;
151
152 }
153
154
164 NumericalDefiniteIntegral::NumericalDefiniteIntegral
165 ( const AbsFunction& function ,
166 const size_t index ,
167 const double a ,
168 const double b ,
169 const NumericalDefiniteIntegral::Points& points ,
170 const double epsabs ,
171 const double epsrel ,
172 const size_t size )
173 : AbsFunction ()
174 , m_function ( function.clone() )
175 , m_DIM ( 0 )
176 , m_index ( index )
177 , m_a ( a )
178 , m_b ( b )
179 , m_ia ( false )
180 , m_ib ( false )
181 , m_type ( GaudiMath::Integration:: Other )
182 , m_category ( GaudiMath::Integration:: Singular )
183 , m_rule ( GaudiMath::Integration:: Fixed )
184 , m_points ( points )
185 , m_pdata ( 0 )
186 , m_epsabs ( epsabs )
187 , m_epsrel ( epsrel )
188 //
189 , m_result ( GSL_NEGINF )
190 , m_error ( GSL_POSINF )
191 //
192 , m_size ( size )
193 , m_ws ( 0 )
194 , m_argument ()
195 , m_argF ()
196 {
197 if ( function.dimensionality() < 2 )
198 { Exception("::constructor: invalid dimensionality ") ; }
199 if ( m_index >= function.dimensionality() )
200 { Exception("::constructor: invalid variable index") ; }
201
202 m_DIM = function.dimensionality() - 1 ;
203 m_argument = Argument( m_DIM ) ;
204 m_argF = Argument( m_DIM + 1 ) ;
205
206 const double l1 = std::min ( a , b ) ;
207 const double l2 = std::max ( a , b ) ;
208 m_points.push_back ( l1 ) ;
209 m_points.push_back ( l2 ) ;
210 std::sort ( m_points.begin() , m_points.end() ) ;
211 m_points.erase( std::unique( m_points.begin () ,
212 m_points.end () ) ,
213 m_points.end() );
214
215 Points::iterator lower =
216 std::lower_bound ( m_points.begin () , m_points.end () , l1 ) ;
217 m_points.erase ( m_points.begin () , lower ) ;
218 Points::iterator upper =
219 std::upper_bound ( m_points.begin () , m_points.end () , l2 ) ;
220 m_points.erase ( upper , m_points.end () ) ;
221
222 m_pdata = new double[ m_points.size() ] ;
223 std::copy( m_points.begin() , m_points.end() , m_pdata );
224 }
225
251 NumericalDefiniteIntegral::NumericalDefiniteIntegral
252 ( const AbsFunction& function ,
253 const size_t index ,
254 const double a ,
255 const GaudiMath::Integration::Inf /* b */ ,
256 const double epsabs ,
257 const double epsrel ,
258 const size_t size )
259 : AbsFunction()
260 , m_function ( function.clone() )
261 , m_DIM ( 0 )
262 , m_index ( index )
263 , m_a ( a )
264 , m_b ( GSL_POSINF )
265 , m_ia ( false )
266 , m_ib ( true )
267 , m_type ( GaudiMath::Integration:: Other )
268 , m_category ( GaudiMath::Integration:: Infinite )
269 , m_rule ( GaudiMath::Integration:: Fixed )
270 , m_points ( )
271 , m_pdata ( 0 )
272 , m_epsabs ( epsabs )
273 , m_epsrel ( epsrel )
274 //
275 , m_result ( GSL_NEGINF )
276 , m_error ( GSL_POSINF )
277 //
278 , m_size ( size )
279 , m_ws ( 0 )
280 , m_argument ()
281 , m_argF ()
282 {
283 if ( function.dimensionality() < 2 )
284 { Exception("::constructor: invalid dimensionality ") ; }
285 if ( m_index >= function.dimensionality() )
286 { Exception("::constructor: invalid variable index") ; }
287
288 m_DIM = function.dimensionality() - 1 ;
289 m_argument = Argument( m_DIM ) ;
290 m_argF = Argument( m_DIM + 1 ) ;
291
292 }
293
294
320 NumericalDefiniteIntegral::NumericalDefiniteIntegral
321 ( const AbsFunction& function ,
322 const size_t index ,
323 const GaudiMath::Integration::Inf /* a */ ,
324 const double b ,
325 const double epsabs ,
326 const double epsrel ,
327 const size_t size )
328 : AbsFunction()
329 , m_function ( function.clone() )
330 , m_DIM ( 0 )
331 , m_index ( index )
332 , m_a ( GSL_NEGINF )
333 , m_b ( b )
334 , m_ia ( true )
335 , m_ib ( false )
336 , m_type ( GaudiMath::Integration:: Other )
337 , m_category ( GaudiMath::Integration:: Infinite )
338 , m_rule ( GaudiMath::Integration:: Fixed )
339 , m_points ( )
340 , m_pdata ( 0 )
341 , m_epsabs ( epsabs )
342 , m_epsrel ( epsrel )
343 //
344 , m_result ( GSL_NEGINF )
345 , m_error ( GSL_POSINF )
346 //
347 , m_size ( size )
348 , m_ws ( 0 )
349 , m_argument ()
350 , m_argF ()
351 {
352 if ( function.dimensionality() < 2 )
353 { Exception("::constructor: invalid dimensionality ") ; }
354 if ( m_index >= function.dimensionality() )
355 { Exception("::constructor: invalid variable index") ; }
356
357 m_DIM = function.dimensionality() - 1 ;
358 m_argument = Argument( m_DIM ) ;
359 m_argF = Argument( m_DIM + 1 ) ;
360 }
361
385 NumericalDefiniteIntegral::NumericalDefiniteIntegral
386 ( const AbsFunction& function ,
387 const size_t index ,
388 const float epsabs ,
389 const float epsrel ,
390 const size_t size )
391 : AbsFunction()
392 , m_function ( function.clone() )
393 , m_DIM ( 0 )
394 , m_index ( index )
395 , m_a ( GSL_NEGINF )
396 , m_b ( GSL_POSINF )
397 , m_ia ( true )
398 , m_ib ( true )
399 , m_type ( GaudiMath::Integration:: Other )
400 , m_category ( GaudiMath::Integration:: Infinite )
401 , m_rule ( GaudiMath::Integration:: Fixed )
402 , m_points ( )
403 , m_pdata ( 0 )
404 , m_epsabs ( epsabs )
405 , m_epsrel ( epsrel )
406 //
407 , m_result ( GSL_NEGINF )
408 , m_error ( GSL_POSINF )
409 //
410 , m_size ( size )
411 , m_ws ( 0 )
412 , m_argument ()
413 , m_argF ()
414 {
415 if ( function.dimensionality() < 2 )
416 { Exception("::constructor: invalid dimensionality ") ; }
417 if ( m_index >= function.dimensionality() )
418 { Exception("::constructor: invalid variable index") ; }
419
420 m_DIM = function.dimensionality() - 1 ;
421 m_argument = Argument( m_DIM ) ;
422 m_argF = Argument( m_DIM + 1 ) ;
423
424 }
425
427 NumericalDefiniteIntegral::NumericalDefiniteIntegral
428 ( const NumericalDefiniteIntegral& right )
429 : AbsFunction ()
430 , m_function ( right.m_function->clone() )
431 , m_DIM ( right.m_DIM )
432 , m_index ( right.m_index )
433 , m_a ( right.m_a )
434 , m_b ( right.m_b )
435 , m_ia ( right.m_ia )
436 , m_ib ( right.m_ib )
437 , m_type ( right.m_type )
438 , m_category ( right.m_category )
439 , m_rule ( right.m_rule )
440 , m_points ( right.m_points )
441 , m_pdata ( 0 )
442 , m_epsabs ( right.m_epsabs )
443 , m_epsrel ( right.m_epsrel )
444 , m_result ( GSL_NEGINF )
445 , m_error ( GSL_POSINF )
446 , m_size ( right.m_size )
447 , m_ws ( 0 )
448 , m_argument ( right.m_argument )
449 , m_argF ( right.m_argF )
450 {
451 if( 0 != right.m_pdata )
452 {
453 m_pdata = new double[m_points.size()] ;
454 std::copy( m_points.begin() , m_points.end() , m_pdata );
455 }
456 }
457
458 NumericalDefiniteIntegral::~NumericalDefiniteIntegral()
459 {
460 if( 0 != m_function ) { delete m_function ; m_function = 0 ; }
461 if( 0 != m_pdata ) { delete[] m_pdata ; m_pdata = 0 ; }
462 }
463
464 // ========================================================================
465 // throw the exception
466 // ========================================================================
467 StatusCode NumericalDefiniteIntegral::Exception
468 ( const std::string& message ,
469 const StatusCode& sc ) const
470 {
471 throw GaudiException( "NumericalDefiniteIntegral::" + message ,
472 "*GaudiMath*" , sc ) ;
473 return sc ;
474 }
475 // ========================================================================
476
477 // ========================================================================
479 // ========================================================================
480 double NumericalDefiniteIntegral::operator()
481 ( const double argument ) const
482 {
483 // reset the result and the error
484 m_result = GSL_NEGINF ;
485 m_error = GSL_POSINF ;
486
487 // check the argument
488 if( 1 != m_DIM ) { Exception ( "operator(): invalid argument size " ) ; };
489
490 m_argument[0] = argument ;
491 return (*this) ( m_argument );
492 }
493 // ========================================================================
494
495 // ========================================================================
497 // ========================================================================
498 Genfun::Derivative
499 NumericalDefiniteIntegral::partial ( unsigned int idx ) const
500 {
501 if ( idx >= m_DIM )
502 { Exception ( "::partial(i): invalid variable index " ) ; };
503 //
504 const AbsFunction& aux = NumericalDerivative( *this , idx ) ;
505 return Genfun::FunctionNoop( &aux ) ;
506 }
507 // ========================================================================
508
509
510 // ========================================================================
512 // ========================================================================
513 double NumericalDefiniteIntegral::operator()
514 ( const Argument& argument ) const
515 {
516 // reset the result and the error
517 m_result = GSL_NEGINF ;
518 m_error = GSL_POSINF ;
519
520 // check the argument
521 if( argument.dimension() != m_DIM )
522 { Exception ( "operator(): invalid argument size " ) ; };
523
524 // copy the argument
525
526 for( size_t i = 0 ; i < m_DIM ; ++i )
527 {
528 m_argument [i] = argument [i] ;
529 const size_t j = i < m_index ? i : i + 1 ;
530 m_argF [j] = argument [i] ;
531 };
532
533 // create the helper object
534 GSL_Helper helper( *m_function , m_argF , m_index );
535
536 // use GSL to evaluate the numerical derivative
537 gsl_function F ;
538 F.function = &GSL_Adaptor ;
539 F.params = &helper ;
540 _Function F1 ;
541 F1.fn = &F ;
542
543 if ( GaudiMath::Integration::Infinite == category () )
544 { return QAGI ( &F1 ) ; } // RETURN
545
546 if ( m_a == m_b )
547 {
548 m_result = 0 ; m_error = 0 ; // EXACT
549 return m_result ; // RETURN
550 }
551
552 if ( GaudiMath::Integration::Singular == category () )
553 { return QAGP ( &F1 ) ; } // RETURN
554 else if ( GaudiMath::Integration::Finite == category () )
555 if ( GaudiMath::Integration::NonAdaptive == type () )
556 { return QNG ( &F1 ) ; } // RETURN
557 else if ( GaudiMath::Integration::Adaptive == type () )
558 { return QAG ( &F1 ) ; } // RETURN
559 else if ( GaudiMath::Integration::AdaptiveSingular == type () )
560 { return QAGS ( &F1 ) ; } // RETURN
561 else
562 { Exception ( "::operator(): invalid type " ); }
563 else
564 { Exception ( "::operator(): invalid category " ); }
565
566 return 0 ;
567 }
568 // ========================================================================
569
570 // ========================================================================
572 // ========================================================================
573 NumericalDefiniteIntegral::_Workspace*
574 NumericalDefiniteIntegral::allocate() const
575 {
576 if ( 0 != m_ws ) { return m_ws; }
577 gsl_integration_workspace* aux =
578 gsl_integration_workspace_alloc( size () );
579 if ( 0 == aux ) { Exception ( "allocate()::invalid workspace" ) ; };
580 m_ws = new _Workspace() ;
581 m_ws->ws = aux ;
582 return m_ws ;
583 }
584 // ========================================================================
585
586 // ========================================================================
587 // adaptive integration on infinite interval
588 // ========================================================================
589 double NumericalDefiniteIntegral::QAGI ( _Function* F ) const
590 {
591 // check the argument
592 if ( 0 == F ) { Exception("::QAGI: invalid function"); }
593
594 // allocate workspace
595 if ( 0 == ws() ) { allocate() ; }
596
597 int ierror = 0 ;
598
599 if ( m_ia && m_ib )
600 {
601 ierror = gsl_integration_qagi ( F->fn ,
602 m_epsabs , m_epsrel ,
603 size () , ws()->ws ,
604 &m_result , &m_error ) ;
605 }
606 else if ( m_ia )
607 {
608 ierror = gsl_integration_qagil ( F->fn , m_b ,
609 m_epsabs , m_epsrel ,
610 size () , ws()->ws ,
611 &m_result , &m_error ) ;
612 }
613 else if ( m_ib )
614 {
615 ierror = gsl_integration_qagiu ( F->fn , m_a ,
616 m_epsabs , m_epsrel ,
617 size () , ws()->ws ,
618 &m_result , &m_error ) ;
619 }
620 else
621 { Exception ( "::QAGI: invalid mode" ) ; };
622
623 if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGI" ,
624 __FILE__ , __LINE__ , ierror ) ;}
625
626 return m_result ;
627 }
628 // ========================================================================
629
630 // ========================================================================
631 // adaptive integration with known singular points
632 // ========================================================================
633 double NumericalDefiniteIntegral::QAGP( _Function* F ) const
634 {
635 if( 0 == F ) { Exception("QAGP::invalid function") ; }
636
637 // no known singular points ?
638 if( points().empty() || 0 == m_pdata ) { return QAGS( F ) ; }
639
640 const size_t npts = points().size();
641
642 // use GSL
643 int ierror =
644 gsl_integration_qagp ( F->fn ,
645 m_pdata , npts ,
646 m_epsabs , m_epsrel ,
647 size () , ws()->ws ,
648 &m_result , &m_error ) ;
649
650 if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGI " ,
651 __FILE__ , __LINE__ , ierror ) ; }
652
653 // the sign
654 if ( m_a > m_b ) { m_result *= -1 ; }
655
656 return m_result ;
657 }
658 // ========================================================================
659
660 // ========================================================================
661 // non-adaptive integration
662 // ========================================================================
663 double NumericalDefiniteIntegral::QNG ( _Function* F ) const
664 {
665 if( 0 == F ) { Exception("QNG::invalid function") ; }
666
667 // integration limits
668 const double low = std::min ( m_a , m_b ) ;
669 const double high = std::max ( m_a , m_b ) ;
670
671 size_t neval = 0 ;
672 int ierror =
673 gsl_integration_qng ( F->fn ,
674 low , high ,
675 m_epsabs , m_epsrel ,
676 &m_result , &m_error , &neval ) ;
677
678 if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QNG " ,
679 __FILE__ , __LINE__ , ierror ) ; }
680
681 // sign
682 if ( m_a > m_b ) { m_result *= -1 ; }
683
684 return m_result ;
685 }
686 // ========================================================================
687
688
689 // ========================================================================
690 // simple adaptive integration
691 // ========================================================================
692 double NumericalDefiniteIntegral::QAG ( _Function* F ) const
693 {
694 if( 0 == F ) { Exception("QAG::invalid function") ; }
695
696 // allocate workspace
697 if( 0 == ws () ) { allocate () ; }
698
699 // integration limits
700 const double low = std::min ( m_a , m_b ) ;
701 const double high = std::max ( m_a , m_b ) ;
702
703 int ierror =
704 gsl_integration_qag ( F->fn ,
705 low , high ,
706 m_epsabs , m_epsrel ,
707 size () , (int) rule() , ws ()->ws ,
708 &m_result , &m_error );
709
710 if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAG " ,
711 __FILE__ , __LINE__ , ierror ) ; }
712
713 // sign
714 if ( m_a > m_b ) { m_result *= -1 ; }
715
716 return m_result ;
717 }
718 // ========================================================================
719
720 // ========================================================================
721 // adaptive integration with singularities
722 // ========================================================================
723 double NumericalDefiniteIntegral::QAGS ( _Function* F ) const
724 {
725 if( 0 == F ) { Exception("QAG::invalid function") ; }
726
727 if ( m_a == m_b )
728 {
729 m_result = 0 ;
730 m_error = 0 ; // EXACT !
731 return m_result ;
732 }
733
734 // allocate workspace
735 if( 0 == ws () ) { allocate () ; }
736
737 // integration limits
738 const double low = std::min ( m_a , m_b ) ;
739 const double high = std::max ( m_a , m_b ) ;
740
741 int ierror =
742 gsl_integration_qags ( F->fn ,
743 low , high ,
744 m_epsabs , m_epsrel ,
745 size () , ws()->ws ,
746 &m_result , &m_error );
747
748 if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGS " ,
749 __FILE__ , __LINE__ , ierror ) ; }
750
751 // sign
752 if ( m_a > m_b ) { m_result *= -1 ; }
753
754 return m_result ;
755 }
756 // ========================================================================
757
758
759 } // end of namespace GaudiMathImplementation
760 } // end of namespace Genfun
761
762
763
764 // ============================================================================
765 // The END
766 // ============================================================================