Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

pFq Not compatible with TaylorSeries package #335

Closed
Li-shiyue opened this issue Sep 21, 2023 · 17 comments
Closed

pFq Not compatible with TaylorSeries package #335

Li-shiyue opened this issue Sep 21, 2023 · 17 comments

Comments

@Li-shiyue
Copy link

I tried to use HypergeometricFunctions and TaylorSeries at the same time as follows.

using HypergeometricFunctions
using Plots
using TaylorSeries
x = set_variables("x", order=100)[1];
Base.isless(y::Number, x::TaylorN) = isless(y, x.coeffs[1].coeffs[1])
Base.log1p(x::TaylorN) = log(1+x)
Base.eps(::Type{TaylorN{T}}) where T = eps(T)
Base.isless(x::TaylorN, y::Number) = isless(x.coeffs[1].coeffs[1], y)
HypergeometricFunctions.pFq((2.0,),(5.0,),x)

However, I got the error:
ERROR: MethodError: no method matching AbstractFloat(::TaylorN{Float64})

Closest candidates are:
(::Type{T})(::T) where T<:Number
@ Core boot.jl:792
(::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number}
@ Base char.jl:50
(::Type{T})(::Base.TwicePrecision) where T<:Number
@ Base twiceprecision.jl:266
..

@lbenet
Copy link
Member

lbenet commented Sep 22, 2023

Quoting this comment from #209:

There is author issue:

x = set_variables("x", order=100)[1]
Base.AbstractFloat(x::TaylorN{T}) where {T} = x.coeffs[1].coeffs[1]
_₂F₁(1.0,2.0,3.0,(x.-1))
0.6137056388801092

However,I can only get the first coeff, How to fix it? Can anybody help me

@lbenet
Copy link
Member

lbenet commented Sep 22, 2023

I am not sure why this is happening... Yet, I do notice that you are overloading isless; isn't that already done within the package? See #323. Aside from that, I would need more time to check what's going on...

@Li-shiyue
Copy link
Author

I am not sure why this is happening... Yet, I do notice that you are overloading isless; isn't that already done within the package? See #323. Aside from that, I would need more time to check what's going on...

Thanks for your reply, the "isless"problem is solved.
However,I try to add Base.AbstractFloat(x::TaylorN{T}) where {T} = x.coeffs[1].coeffs[1] to fix the new problem but I still don't understand how to get other coeffs

@lbenet
Copy link
Member

lbenet commented Sep 22, 2023

I think I begin to understand the problem. Somehow you need to obtain the coefficients of the series. To be concrete, let us build a simple example:

julia> x, y = set_variables("x y", order=8)
2-element Vector{TaylorN{Float64}}:
  1.0 x + 𝒪(‖x‖⁹)
  1.0 y + 𝒪(‖x‖⁹)

julia> ee = exp(x-y^3)
 1.0 + 1.0 x + 0.5+ 0.16666666666666666- 1.0+ 0.041666666666666664 x⁴ - 1.0 x y³ + 0.008333333333333333 x⁵ - 0.5 x² y³ + 0.001388888888888889 x⁶ - 0.16666666666666666 x³ y³ + 0.5 y⁶ + 0.0001984126984126984 x⁷ - 0.041666666666666664 x⁴ y³ + 0.5 x y⁶ + 2.48015873015873e-5 x⁸ - 0.008333333333333333 x⁵ y³ + 0.25 x² y⁶ + 𝒪(‖x‖⁹)

Let me assume that you require the coeff of the monomial x³ y³ and x⁵ y³. Then:

julia> getcoeff(ee, (3,3))
-0.16666666666666666

julia> getcoeff(ee, (5,3))
-0.008333333333333333

Another alternative is to get the homogeneous polynomial of specific degree, and then from the way they are ordered, the specific monomial. In order to know the actual order, use show_monomials(order) where order is the degree of the homogeneous polynomial. For the first example above:

julia> show_monomials(6)
 1  -->  x⁶
 2  -->  x⁵ y
 3  -->  x⁴ y²
 4  -->  x³ y³
 5  -->  x² y⁴
 6  -->  x y⁵
 7  -->  y⁶

julia> ee[6] # returns the homogeneous polynomial of degree 6
 0.001388888888888889 x⁶ - 0.16666666666666666 x³ y³ + 0.5 y⁶

julia> ee[6][4] # 4th coef, of homogeneous polynomial of order 6
-0.16666666666666666

The last line is equivalent to ee.coeffs[7].coeffs[4]; 7 = 6+1 (order of homogeneous polynomial +1).

Is this what you are after?

@lbenet
Copy link
Member

lbenet commented Sep 22, 2023

Two comments in your code here:

  • log1p is also defined, so you do not need to overload it
  • It seems you only care about one-independent variable: x = set_variables("x", order=100)[1]. If this is the case, I think it is better to use Taylor1: x = Taylor1(100).

@Li-shiyue
Copy link
Author

Two comments in your code here:

  • log1p is also defined, so you do not need to overload it
  • It seems you only care about one-independent variable: x = set_variables("x", order=100)[1]. If this is the case, I think it is better to use Taylor1: x = Taylor1(100).

My problem is I tried to use TaylorSeries on hyperfunctin.here is my code
x = set_variables("x",order=10)
_₂F₁(1.0, 2.0, 3.0, (x-1))
ERROR: MethodError: no method matching AbstractFloat(::TaylorN{Float64})
Closest candidates are:
(::Type{T})(::T) where T<:Number
@ Core boot.jl:792
(::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number}
@ Base char.jl:50
(::Type{T})(::Base.TwicePrecision) where T<:Number
@ Base twiceprecision.jl:266
...
So I add
x = set_variables("x",order=10)
Base.AbstractFloat(x::TaylorN) = x.coeffs[1].coeffs[1]
_₂F₁(1.0, 2.0, 3.0, (x-1))
then I got the first coeff,but I need rest of them

@Li-shiyue
Copy link
Author

I think I begin to understand the problem. Somehow you need to obtain the coefficients of the series. To be concrete, let us build a simple example:

julia> x, y = set_variables("x y", order=8)
2-element Vector{TaylorN{Float64}}:
  1.0 x + 𝒪(‖x‖⁹)
  1.0 y + 𝒪(‖x‖⁹)

julia> ee = exp(x-y^3)
 1.0 + 1.0 x + 0.5+ 0.16666666666666666- 1.0+ 0.041666666666666664 x⁴ - 1.0 x y³ + 0.008333333333333333 x⁵ - 0.5 x² y³ + 0.001388888888888889 x⁶ - 0.16666666666666666 x³ y³ + 0.5 y⁶ + 0.0001984126984126984 x⁷ - 0.041666666666666664 x⁴ y³ + 0.5 x y⁶ + 2.48015873015873e-5 x⁸ - 0.008333333333333333 x⁵ y³ + 0.25 x² y⁶ + 𝒪(‖x‖⁹)

Let me assume that you require the coeff of the monomial x³ y³ and x⁵ y³. Then:

julia> getcoeff(ee, (3,3))
-0.16666666666666666

julia> getcoeff(ee, (5,3))
-0.008333333333333333

Another alternative is to get the homogeneous polynomial of specific degree, and then from the way they are ordered, the specific monomial. In order to know the actual order, use show_monomials(order) where order is the degree of the homogeneous polynomial. For the first example above:

julia> show_monomials(6)
 1  -->  x⁶
 2  -->  x⁵ y
 3  -->  x⁴ y²
 4  -->  x³ y³
 5  -->  x² y⁴
 6  -->  x y⁵
 7  -->  y⁶

julia> ee[6] # returns the homogeneous polynomial of degree 6
 0.001388888888888889 x⁶ - 0.16666666666666666 x³ y³ + 0.5 y⁶

julia> ee[6][4] # 4th coef, of homogeneous polynomial of order 6
-0.16666666666666666

The last line is equivalent to ee.coeffs[7].coeffs[4]; 7 = 6+1 (order of homogeneous polynomial +1).

Is this what you are after?

however getcoeff can not be used on hyperfunction

@lbenet
Copy link
Member

lbenet commented Sep 23, 2023

I do not know HypergeometricFunctions, nor your actual application. It would be helpful if you provide details of what is your goal. My guess is that you want to produce the series itself, right?

@lbenet
Copy link
Member

lbenet commented Sep 23, 2023

The following seems to work... though I still do not know if this is what you are after:

julia> using HypergeometricFunctions

julia> using TaylorSeries

julia> t = Taylor1(5);

julia> pFq((0.0,), (0.0,), 0.1)
1.1051709180756477

julia> Base.float(t::Taylor1) = Taylor1(float(t.coeffs))

julia> pFq((0.0,), (0.0,), 0.1+t)
 1.1051709180756477 + 1.1051709180756477 t + 0.5525854590378239+ 0.18419515301260794+ 0.046048788253151986 t⁴ + 0.009209757650630397 t⁵ + 𝒪(t⁶)

@Li-shiyue
Copy link
Author

Base.float(t::Taylor1) = Taylor1(float(t.coeffs))

Thanks,that helps me!

@Li-shiyue
Copy link
Author

I do not know HypergeometricFunctions, nor your actual application. It would be helpful if you provide details of what is your goal. My guess is that you want to produce the series itself, right?

function bbb(params,t)
a,b,c = params
s = a + b
return pFq((a,),(s,),c*(t-1))
end
t = Taylor1(5)
Base.float(t::Taylor1) = Taylor1(float(t.coeffs))
aaa([1,2,4],t)

ERROR: MethodError: no method matching eps(::Type{Taylor1{Float64}})

Closest candidates are:
eps()
@ Base float.jl:957
eps(::FixedPointNumbers.FixedPoint)
@ FixedPointNumbers ~/.julia/packages/FixedPointNumbers/HAGk2/src/FixedPointNumbers.jl:222
eps(::DualNumbers.Dual)
@ DualNumbers ~/.julia/packages/DualNumbers/5knFX/src/dual.jl:57
...
How can I solved eps error? My goal is to get the Taylor series

@edwardcao3026
Copy link

I think he is interested in finding series coefficients to a customized function involving Hypergeometric functions.

I do not know HypergeometricFunctions, nor your actual application. It would be helpful if you provide details of what is your goal. My guess is that you want to produce the series itself, right?

function bbb(params,t) a,b,c = params s = a + b return pFq((a,),(s,),c*(t-1)) end t = Taylor1(5) Base.float(t::Taylor1) = Taylor1(float(t.coeffs)) aaa([1,2,4],t)

ERROR: MethodError: no method matching eps(::Type{Taylor1{Float64}})

Closest candidates are: eps() @ Base float.jl:957 eps(::FixedPointNumbers.FixedPoint) @ FixedPointNumbers ~/.julia/packages/FixedPointNumbers/HAGk2/src/FixedPointNumbers.jl:222 eps(::DualNumbers.Dual) @ DualNumbers ~/.julia/packages/DualNumbers/5knFX/src/dual.jl:57 ... How can I solved eps error? My goal is to get the Taylor series

@Li-shiyue
Copy link
Author

I do not know HypergeometricFunctions, nor your actual application. It would be helpful if you provide details of what is your goal. My guess is that you want to produce the series itself, right?

Yes, I tried to get Taylor coefficients series with hypergeometric functions
I define a function but get the eps error

t = Taylor1(5)
Base.float(t::Taylor1) = Taylor1(float(t.coeffs))
pFq((1.0,),(2.0,),t-1)
function f(a,b,t)
return pFq((a,),(b,),(t-1))*pFq((b,),(a,),(t+1))
end
f(2,3,t)
ERROR: MethodError: no method matching eps(::Type{Taylor1{Float64}})

Closest candidates are:
eps()
@ Base float.jl:957
eps(!Matched::FixedPointNumbers.FixedPoint)
@ FixedPointNumbers ~/.julia/packages/FixedPointNumbers/HAGk2/src/FixedPointNumbers.jl:222
eps(!Matched::DualNumbers.Dual)
@ DualNumbers ~/.julia/packages/DualNumbers/5knFX/src/dual.jl:57
...
I don't know how to solved it

@lbenet
Copy link
Member

lbenet commented Sep 23, 2023

I will have a look later... I need to see how is it used in HyperGeometricFunctions. I am not sure if by using a similar trick as I did before is enough or not, or if the result is supposed to be a scalar.

@lbenet
Copy link
Member

lbenet commented Sep 25, 2023

First, let me indicate where the problem appears, which is not included in your message above; I essentially use the same code you have:

julia> using HypergeometricFunctions

julia> using TaylorSeries

julia> Base.float(t::Taylor1) = Taylor1(float(t.coeffs))

julia> t = Taylor1(5)
 1.0 t + 𝒪(t⁶)

julia> f1 = pFq((1.0,), (2.0,), t-1) # seemingly ok
 0.6321205588285577 + 0.26424111765711533 t + 0.08030139707139416+ 0.018988156876153774+ 0.003659846827343676 t⁴ + 0.0005941848175816562 t⁵ + 𝒪(t⁶)

julia> f2 = pFq((3.0,), (2.0,), t-1) # throws error related to eps
ERROR: MethodError: no method matching eps(::Type{Taylor1{Float64}})

Closest candidates are:
  eps()
   @ Base float.jl:957
  eps(::AbstractFloat)
   @ Base float.jl:953
  eps(::T) where T<:Dates.TimeType
   @ Dates /opt/julia/julia-1.9.2/share/julia/stdlib/v1.9/Dates/src/types.jl:435
  ...

Stacktrace:
 [1] pFqweniger::Tuple{Float64}, β::Tuple{Float64}, z::Taylor1{Float64}; kmax::Int64)
   @ HypergeometricFunctions ~/.julia/packages/HypergeometricFunctions/LZbIV/src/weniger.jl:158
...

The reason I copy-pasted part of the error message is because it points out to the line where the error is thrown: in this case, it is this line of HypergeometricFunctions.jl, where eps is used for a Taylor1 of for the type Taylor1. If you check the code, it seems that some precision check is needed there, which is broken for Taylor1s.

So I extended eps in what seems to me the natural way:

julia> Base.eps(a::Taylor1) = eps(constant_term(a))

julia> Base.eps(::Type{Taylor1{T}}) where {T<:Number} = eps(T)

With these changes, it seems to work:

julia> f2 = pFq((3.0,),(2.0,),t-1)
 0.1839397205857213 + 0.36787944117141885 t + 0.27590958088092615+ 0.12262648021343794+ 0.038320783898288496 t⁴ + 0.009196673989352056 t⁵ + 𝒪(t⁶)

I am not sure if the result is correct/consistent; you should check that with known concrete values, because it could be that a different output of the eps functions may be needed to ensure convergence. (The repo includes some references, which I guess are worth reading to check the precision issues; I haven't done so.)

I hope this helps!

@Li-shiyue
Copy link
Author

First, let me indicate where the problem appears, which is not included in your message above; I essentially use the same code you have:

julia> using HypergeometricFunctions

julia> using TaylorSeries

julia> Base.float(t::Taylor1) = Taylor1(float(t.coeffs))

julia> t = Taylor1(5)
 1.0 t + 𝒪(t⁶)

julia> f1 = pFq((1.0,), (2.0,), t-1) # seemingly ok
 0.6321205588285577 + 0.26424111765711533 t + 0.08030139707139416+ 0.018988156876153774+ 0.003659846827343676 t⁴ + 0.0005941848175816562 t⁵ + 𝒪(t⁶)

julia> f2 = pFq((3.0,), (2.0,), t-1) # throws error related to eps
ERROR: MethodError: no method matching eps(::Type{Taylor1{Float64}})

Closest candidates are:
  eps()
   @ Base float.jl:957
  eps(::AbstractFloat)
   @ Base float.jl:953
  eps(::T) where T<:Dates.TimeType
   @ Dates /opt/julia/julia-1.9.2/share/julia/stdlib/v1.9/Dates/src/types.jl:435
  ...

Stacktrace:
 [1] pFqweniger::Tuple{Float64}, β::Tuple{Float64}, z::Taylor1{Float64}; kmax::Int64)
   @ HypergeometricFunctions ~/.julia/packages/HypergeometricFunctions/LZbIV/src/weniger.jl:158
...

The reason I copy-pasted part of the error message is because it points out to the line where the error is thrown: in this case, it is this line of HypergeometricFunctions.jl, where eps is used for a Taylor1 of for the type Taylor1. If you check the code, it seems that some precision check is needed there, which is broken for Taylor1s.

So I extended eps in what seems to me the natural way:

julia> Base.eps(a::Taylor1) = eps(constant_term(a))

julia> Base.eps(::Type{Taylor1{T}}) where {T<:Number} = eps(T)

With these changes, it seems to work:

julia> f2 = pFq((3.0,),(2.0,),t-1)
 0.1839397205857213 + 0.36787944117141885 t + 0.27590958088092615+ 0.12262648021343794+ 0.038320783898288496 t⁴ + 0.009196673989352056 t⁵ + 𝒪(t⁶)

I am not sure if the result is correct/consistent; you should check that with known concrete values, because it could be that a different output of the eps functions may be needed to ensure convergence. (The repo includes some references, which I guess are worth reading to check the precision issues; I haven't done so.)

I hope this helps!

This works for me!!! Thanks!
but, I am still curious about the difference between TaylorN and Taylor1.
so I tried set_variables("x y", order=8) to find the two variables coefficient
The supplement above didn't work.
function bbb(params)
x,y = params
return _₁F₁(2,3,x+0.5*y-1)
end
x, y = set_variables("x y", order=8)
bbb([x,y])

@lbenet
Copy link
Member

lbenet commented Oct 3, 2023

Closing this because the discussion is following in #336; feel free to reopen if you think there is more to discuss here...

@lbenet lbenet closed this as completed Oct 4, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants