From 75c3c220699884b542b5dc11babfd7577e559813 Mon Sep 17 00:00:00 2001 From: Vanya Belyaev Date: Fri, 13 Dec 2024 16:07:30 +0100 Subject: [PATCH] 1. add conversion from Taylor(polynomial) expansing to Pade approximation --- ReleaseNotes/release_notes.md | 3 +- source/include/Ostap/Rational.h | 43 +++++-- source/src/GSL_helpers.cpp | 23 +++- source/src/GSL_helpers.h | 32 +++-- source/src/Rational.cpp | 217 +++++++++++++++++++++++++++----- 5 files changed, 269 insertions(+), 49 deletions(-) diff --git a/ReleaseNotes/release_notes.md b/ReleaseNotes/release_notes.md index 59eacaa9..b2edaa1c 100644 --- a/ReleaseNotes/release_notes.md +++ b/ReleaseNotes/release_notes.md @@ -2,7 +2,8 @@ 1. add several functions, like `sec` , `csc` , `cot` , `cas` 1. add Bring radical - 1. add Pade intrpolation + 1. add Pade interpolation + 1. add conversion from Taylor(polynomial) expansing to Pade approximation ## Backward incompatible diff --git a/source/include/Ostap/Rational.h b/source/include/Ostap/Rational.h index 667faea0..75b35320 100644 --- a/source/include/Ostap/Rational.h +++ b/source/include/Ostap/Rational.h @@ -457,6 +457,9 @@ namespace Ostap // ====================================================================== } ; // ======================================================================== + /// forward declaration + class Polynomial ; // forward declatation + // ======================================================================== /** @class Pade * Pade-like rational function with the optional setting * of constituent shape-fixing zeroes and poles @@ -587,6 +590,36 @@ namespace Ostap const std::vector >& czeroes = std::vector > () , const std::vector >& cpoles = std::vector > () ) ; // ===================================================================== + /** constructor from the polynomial/Taylor expansion + * @param p polynomial expansion + * @param n degree of P(x) + * @param m degree of Q(x) + */ + Pade + ( const Ostap::Math::Polynomial& p , + const unsigned short n , + const unsigned short m ) ; + // ===================================================================== + /** constructor from the polynomial/Taylor expansion + * @param p polynomial expansion + * @param n degree of P(x) + */ + Pade + ( const Ostap::Math::Polynomial& p , + const unsigned short n ) ; + // ===================================================================== + /** constructor from the polynomial/Taylor expansion + * @param n degree of P(x) + * @param m degree of Q(x) + * @param p polynomial expansion + */ + Pade + ( const unsigned short n , + const unsigned short m , + const std::vector& p , + const double xmin = -1 , + const double xmax = 1 ) ; + // ===================================================================== public: // ====================================================================== /// get x-min @@ -734,14 +767,10 @@ namespace Ostap // ====================================================================== } // ======================================================================== - - // - void pade_intepolation - ( const std::vector& x , - const std::vector& y , - const unsigned short n ) ; - // + // Ostap::Math::Pade + void pade_ququ ( const std::vector& f , + const unsigned short n ) ; } // The end of namespace Ostap::Math diff --git a/source/src/GSL_helpers.cpp b/source/src/GSL_helpers.cpp index 8929bc52..b18f0acc 100644 --- a/source/src/GSL_helpers.cpp +++ b/source/src/GSL_helpers.cpp @@ -1,6 +1,10 @@ // ============================================================================ // Include files // ============================================================================ +// Ostap +// ============================================================================ +#include "Ostap/GSL_utils.h" +// ============================================================================ // GSL // ============================================================================ #include "gsl/gsl_linalg.h" @@ -161,9 +165,22 @@ Ostap::GSL_Permutation::~GSL_Permutation () m_permutation = nullptr ; } } - - - +// ============================================================================ +/// print operator +std::ostream& operator<<( std::ostream& s , + const Ostap::GSL_Matrix& m ) +{ + if ( m.matrix() == nullptr ) { s << "GSL_Matrix{nullptr}" ; return s ; } + return Ostap::Utils::toStream ( *m.matrix() , s ) ; +} +// ============================================================================ +/// print operator +std::ostream& operator<<( std::ostream& s , + const Ostap::GSL_Vector& v ) +{ + if ( v.vector() == nullptr ) { s << "GSL_Vector{nullptr}" ; return s ; } + return Ostap::Utils::toStream ( *v.vector() , s ) ; +} // ============================================================================ // The END // ============================================================================ diff --git a/source/src/GSL_helpers.h b/source/src/GSL_helpers.h index 2f039d59..b4ce4ecd 100644 --- a/source/src/GSL_helpers.h +++ b/source/src/GSL_helpers.h @@ -4,6 +4,10 @@ // ============================================================================ // Include files // ============================================================================ +// STD/STD +// ============================================================================ +#include +// ============================================================================ // GSL // ============================================================================ #include "gsl/gsl_linalg.h" @@ -68,13 +72,15 @@ namespace Ostap inline const gsl_matrix* matrix () const { return m_matrix ; } // ======================================================================== /// get the mattix element - double get ( const unsigned short n1 , - const unsigned short n2 ) const + double get + ( const unsigned short n1 , + const unsigned short n2 ) const { return gsl_matrix_get ( m_matrix , n1 , n2 ) ; } /// set the matrix element - void set ( const unsigned short n1 , - const unsigned short n2 , - const double value ) + void set + ( const unsigned short n1 , + const unsigned short n2 , + const double value ) { gsl_matrix_set ( m_matrix , n1 , n2 , value ) ; } // ======================================================================== private: @@ -125,11 +131,13 @@ namespace Ostap inline const gsl_vector* vector () const { return m_vector ; } // ======================================================================== /// get the vector element - double get ( const unsigned short n ) const + double get + ( const unsigned short n ) const { return gsl_vector_get ( m_vector , n ) ; } /// set the vector element - void set ( const unsigned short n , - const double value ) + void set + ( const unsigned short n , + const double value ) { gsl_vector_set ( m_vector , n , value ) ; } // ======================================================================== private: @@ -173,7 +181,13 @@ namespace Ostap // ======================================================================== }; // ========================================================================== -} // The end of namespace Ostap +} // The end of namespace Ostap +// ============================================================================ +/// print operator +std::ostream& operator<<( std::ostream& s , const Ostap::GSL_Matrix& m ) ; +// ============================================================================ +/// print operator +std::ostream& operator<<( std::ostream& s , const Ostap::GSL_Vector& v ) ; // ============================================================================ // The END // ============================================================================ diff --git a/source/src/Rational.cpp b/source/src/Rational.cpp index 149b4542..eaacba28 100644 --- a/source/src/Rational.cpp +++ b/source/src/Rational.cpp @@ -5,6 +5,7 @@ // Ostap // ============================================================================ #include "Ostap/Interpolation.h" +#include "Ostap/Polynomials.h" #include "Ostap/Rational.h" #include "Ostap/Clenshaw.h" #include "Ostap/GSL_utils.h" @@ -910,8 +911,6 @@ void Ostap::Math::Pade::swap ( Ostap::Math::Pade& right ) // Ostap::Math::swap ( m_workspace , right.m_workspace ) ; } -// ========================================================================== - // ========================================================================== /* Interpolatory constructor * @param table interpolation table @@ -953,7 +952,7 @@ Ostap::Math::Pade::Pade "Ostap::Math::Pade" ) ; // // ===================================================================== - Ostap::GSL_Matrix A { N , N } ; + Ostap::GSL_Matrix A { N , N , Ostap::GSL_Matrix::Zero() } ; Ostap::GSL_Vector b { N } ; // free column // (2) fill the matrix A and free column b for ( unsigned short j = 0 ; j < N ; ++j ) @@ -1010,35 +1009,195 @@ Ostap::Math::Pade::Pade // (4) Feed Pade with calculated parameters for ( unsigned short k = 0 ; k < N ; ++k ) { setPar ( k , x.get ( k ) ) ; } } -// ========================================================================== -/* create Pade function that interpolates the data - * @param table interpolation table - * @param n degree of P(x) - * @param zeros list of real constituent zeroes - * @param poles list of real constituent poles - * @param czeros (half) list of complex constituent zeroes - * @param cpoles (half) list of complex constituent poles + + +// =========================================================================== +/** constructor from the polynomial expansion + * @param n degree of P(x) + * @param m degree of Q(x) + * @param f polynomial expansion */ -// ========================================================================== -Ostap::Math::Pade -Ostap::Math::Interpolation::pade -( const Ostap::Math::Interpolation::Table& table , - const unsigned short n , - const std::vector& zeroes , - const std::vector& poles , - const std::vector >& czeroes , - const std::vector >& cpoles ) +// =========================================================================== +Ostap::Math::Pade::Pade +( const unsigned short n , + const unsigned short m , + const std::vector& f , + const double xmin , + const double xmax ) + : Pade ( std::vector ( m + n + 1 , 0.0 ) , + n , {} , {} , {} , {} , xmin , xmax ) { - Ostap::Assert ( n + 1 <= table.size() , - "Data table is too short!" , - "Ostap::Math::Interpolation::pade" ) ; - return Ostap::Math::Pade ( table , - n , - zeroes , - poles , - czeroes , - cpoles ) ; + Ostap::Assert ( n + m + 1 <= f.size () , + "Invalid Polynomial->Pade setting!" , + "Ostap::Math::Pade" ) ; + // + const unsigned short N = this->npars () ; + // + typedef std::vector::const_iterator CI ; + typedef std::vector::const_reverse_iterator CRI ; + // + Ostap::GSL_Matrix A { N , N , Ostap::GSL_Matrix::Zero() } ; + Ostap::GSL_Vector b { N } ; // free column + // + for ( unsigned short j = 0 ; j < N ; ++j ) + { + if ( j < n + 1 ) { A.set ( j , j , 1 ) ; } + // + if ( 1 <= j && j < m + 1 ) + { + CI i1 = f.begin() ; + CRI i2 { f.begin() + j } ; + for ( unsigned short k = 0 ; k < j ; ++k ) + { + A.set ( j , n + 1 + k , -1 * ( *i2++ ) ) ; + } + } + if ( m + 1 <= j ) + { + CRI i1 { f.begin() + ( j - m ) } ; + CRI i2 { f.begin() + j } ; + for ( unsigned short k = 0 ; k < m ; ++k ) + { + A.set ( j , n + 1 + k , -1 * ( *i2++ ) ) ; + } + } + b.set( j , f[j] ) ; + } + // =============================================================================== + // (3) solve the system Ax=b using LU decomposition with pivoting + + // (3.1) make LU decomposition with pivoting + Ostap::GSL_Permutation P { N } ; + int signum ; + int ierror = gsl_linalg_LU_decomp ( A.matrix() , P.permutation() , &signum ) ; + if ( ierror ) { gsl_error ( "Failure in LU-decomposition" , __FILE__ , __LINE__ , ierror ) ; } + Ostap::Assert ( !ierror , + "Failure in LU-decomposition!" , + "Ostap::Math::Pade" , 1100 + ierror ) ; + + // (3.2) solve the system Ax=b + Ostap::GSL_Vector x { N } ; // solution + ierror = gsl_linalg_LU_solve ( A.matrix(), P.permutation() , b.vector() , x.vector() ); + if ( ierror ) { gsl_error ( "Failure in LU-solve" , __FILE__ , __LINE__ , ierror ) ; } + Ostap::Assert ( !ierror , + "Failure in LU-solve!" , + "Ostap::Math::Pade" , 1200 + ierror ) ; + + // (4) Feed Pade with calculated parameters + for ( unsigned short k = 0 ; k < N ; ++k ) { setPar ( k , x.get ( k ) ) ; } + // ========================================================================== } +// ============================================================================ +/* constructor from polynomi expansion + * @param p polynomial expansion + * @param n degree of P(x) + */ +// ============================================================================= +Ostap::Math::Pade::Pade +( const Ostap::Math::Polynomial& p , + const unsigned short n , + const unsigned short m ) + : Pade ( n , m , p.pars() , p.xmin() , p.xmax() ) +{} +// =============================================================================== +/* constructor from the polynomial expansion + * @param p polynomial expansion + * @param n degree of P(x) + */ +// =============================================================================== +Ostap::Math::Pade::Pade +( const Ostap::Math::Polynomial& p , + const unsigned short n ) + : Pade ( p , n , n + 1 <= p.npars() ? p.npars() - n - 1 : 0 ) +{ + Ostap::Assert ( n + 1 <= p.npars () , + "Invalid Poliynomial->Pade setting!" , + "Ostap::Math::Pade" ) ; +} + + + +/** + +#include "Ostap/GSL_utils.h" +#include "GSL_helpers.h" +#include + +// Ostap::Math::Pade +void Ostap::Math::pade_ququ +( const std::vector& f , + const unsigned short n ) +{ + Ostap::Assert ( n + 1 <= f.size() , + "Insufficient array of coefficients!" + "Ostap::Math::pade" ) ; + + // + const unsigned short N = f.size() ; + const unsigned short m = f.size() - n - 1 ; + // + typedef std::vector::const_iterator CI ; + typedef std::vector::const_reverse_iterator CRI; + // + Ostap::GSL_Matrix A { N , N , Ostap::GSL_Matrix::Zero() } ; + Ostap::GSL_Vector b { N } ; // free column + // + for ( unsigned short j = 0 ; j < N ; ++j ) + { + if ( j < n + 1 ) { A.set ( j , j , 1 ) ; } + // + if ( 1 <= j && j < m + 1 ) + { + CI i1 = f.begin() ; + CRI i2 { f.begin() + j } ; + for ( unsigned short k = 0 ; k < j ; ++k ) + { + A.set ( j , n + 1 + k , -1 * ( *i2++ ) ) ; + } + } + if ( m + 1 <= j ) + { + CRI i1 { f.begin() + ( j - m ) } ; + CRI i2 { f.begin() + j } ; + for ( unsigned short k = 0 ; k < m ; ++k ) + { + A.set ( j , n + 1 + k , -1 * ( *i2++ ) ) ; + } + } + b.set( j , f[j] ) ; + } + // solve Ax=b + + + // (3) solve the system Ax=b using LU decomposition with pivoting + + // (3.1) make LU decomposition with pivoting + Ostap::GSL_Permutation P { N } ; + int signum ; + int ierror = gsl_linalg_LU_decomp ( A.matrix() , P.permutation() , &signum ) ; + if ( ierror ) { gsl_error ( "Failure in LU-decomposition" , __FILE__ , __LINE__ , ierror ) ; } + Ostap::Assert ( !ierror , + "Failure in LU-decomposition!" , + "Ostap::Math::Pade" , 1100 + ierror ) ; + + // (3.2) solve the system Ax=b + Ostap::GSL_Vector x { N } ; // solution + ierror = gsl_linalg_LU_solve ( A.matrix(), P.permutation() , b.vector() , x.vector() ); + if ( ierror ) { gsl_error ( "Failure in LU-solve" , __FILE__ , __LINE__ , ierror ) ; } + Ostap::Assert ( !ierror , + "Failure in LU-solve!" , + "Ostap::Math::Pade" , 1200 + ierror ) ; + + // (4) Feed Pade with calculated parameters + // for ( unsigned short k = 0 ; k < N ; ++k ) { setPar ( k , x.get ( k ) ) ; } + + + std::cout << (*A.matrix()) << std::endl ; + std::cout << (*x.vector()) << std::endl ; +} + +*/ + // ============================================================================ // The END // ============================================================================