Skip to content

Commit cf27073

Browse files
jishnubjishnub
authored
Move SpecialFunctions to an extension (#564)
1 parent f7385d5 commit cf27073

File tree

4 files changed

+176
-89
lines changed

4 files changed

+176
-89
lines changed

‎Project.toml

Lines changed: 9 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
name = "ApproxFunBase"
22
uuid = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"
3-
version = "0.9.9"
3+
version = "0.9.10"
44

55
[deps]
66
AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
@@ -24,6 +24,12 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
2424
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
2525
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
2626

27+
[weakdeps]
28+
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
29+
30+
[extensions]
31+
ApproxFunBaseSpecialFunctionsExt = "SpecialFunctions"
32+
2733
[compat]
2834
AbstractFFTs = "0.5, 1"
2935
Aqua = "0.6"
@@ -50,7 +56,8 @@ julia = "1.9"
5056
Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
5157
Infinities = "e1ba4f0e-776d-440f-acd9-e1d2e9742647"
5258
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
59+
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
5360
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
5461

5562
[targets]
56-
test = ["Aqua", "Random", "Infinities", "Test"]
63+
test = ["Aqua", "Random", "Infinities", "Test", "SpecialFunctions"]

‎ext/ApproxFunBaseSpecialFunctionsExt.jl

Lines changed: 152 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,152 @@
1+
module ApproxFunBaseSpecialFunctionsExt
2+
3+
using ApproxFunBase
4+
using SpecialFunctions
5+
import Calculus
6+
7+
# we need to import all special functions to use Calculus.symbolic_derivatives_1arg
8+
# we can't do importall Base as we replace some Base definitions
9+
import SpecialFunctions: airy, besselh,
10+
lfact, beta, lbeta,
11+
eta, zeta, polygamma, logabsgamma, loggamma,
12+
besselj, bessely, besseli, besselk, besselkx,
13+
hankelh1, hankelh2, hankelh1x, hankelh2x,
14+
# functions from Calculus.symbolic_derivatives_1arg
15+
erf, erfinv, erfc, erfcinv, erfi, gamma, lgamma,
16+
digamma, invdigamma, trigamma,
17+
airyai, airybi, airyaiprime, airybiprime,
18+
besselj0, besselj1, bessely0, bessely1,
19+
erfcx, dawson
20+
21+
besselh(ν,k::Integer,f::Fun) = k == 1 ? hankelh1(ν,f) : k == 2 ? hankelh2(ν,f) : throw(Base.Math.AmosException(1))
22+
23+
for jy in (:j, :y)
24+
bjy = Symbol(:bessel, jy)
25+
for ν in (0,1)
26+
bjynu = Symbol(bjy, ν)
27+
@eval SpecialFunctions.$bjynu(f::Fun) = $bjy($ν,f)
28+
end
29+
end
30+
31+
# Other special functions
32+
for f in [:logabsgamma]
33+
@eval function $f(z::Fun{<:ConstantSpace, <:Real})
34+
t = $f(Number(z))
35+
Fun(t[1], space(z)), t[2]
36+
end
37+
end
38+
function loggamma(z::Fun{<:ConstantSpace})
39+
t = loggamma(Number(z))
40+
Fun(t, space(z))
41+
end
42+
for f in [:gamma, :loggamma]
43+
@eval begin
44+
function $f(a, z::Fun{<:ConstantSpace})
45+
t = $f(a, Number(z))
46+
Fun(t, space(z))
47+
end
48+
end
49+
end
50+
51+
for f in [:besselj, :besselk, :besselkx, :bessely, :besseli,
52+
:hankelh1x, :hankelh2x, :hankelh1, :hankelh2]
53+
@eval $f(nu, x::Fun{<:ConstantSpace}) = Fun($f(nu, Number(x)), space(x))
54+
end
55+
56+
57+
for (funsym, exp) in Calculus.symbolic_derivatives_1arg()
58+
if isdefined(SpecialFunctions, funsym)
59+
@eval begin
60+
$(funsym)(z::Fun{<:ConstantSpace,<:Real}) =
61+
Fun($(funsym)(Number(z)),space(z))
62+
$(funsym)(z::Fun{<:ConstantSpace,<:Complex}) =
63+
Fun($(funsym)(Number(z)),space(z))
64+
$(funsym)(z::Fun{<:ConstantSpace}) =
65+
Fun($(funsym)(Number(z)),space(z))
66+
end
67+
end
68+
end
69+
70+
for op in (:(lfact),:(erfinv),:(erfcinv),:(beta),:(lbeta),
71+
:(eta),:(zeta),:(gamma),:(lgamma),
72+
:(polygamma),:(invdigamma),:(digamma),:(trigamma))
73+
@eval begin
74+
$op(f::Fun) = Fun($op ∘ f,domain(f))
75+
end
76+
end
77+
78+
79+
function airy(k::Number,f::Fun)
80+
if k == 0
81+
airyai(f)
82+
elseif k == 1
83+
airyaiprime(f)
84+
elseif k == 2
85+
airybi(f)
86+
elseif k == 3
87+
airybiprime(f)
88+
else
89+
error("invalid argument")
90+
end
91+
end
92+
93+
erfc(f::Fun) = 1-erf(f)
94+
95+
for (op,ODE,RHS,growth) in ((:(erf),"f'*D^2+(2f*f'^2-f'')*D","0",:(imag)),
96+
(:(erfi),"f'*D^2-(2f*f'^2+f'')*D","0",:(real)),
97+
(:(airyai),"f'*D^2-f''*D-f*f'^3","0",:(imag)),
98+
(:(airybi),"f'*D^2-f''*D-f*f'^3","0",:(imag)),
99+
(:(airyaiprime),"f'*D^2-f''*D-f*f'^3","airyai(f)*f'^3",:(imag)),
100+
(:(airybiprime),"f'*D^2-f''*D-f*f'^3","airybi(f)*f'^3",:(imag)))
101+
L,R = Meta.parse(ODE),Meta.parse(RHS)
102+
@eval begin
103+
function $op(fin::Fun)
104+
T = ApproxFunBase.cfstype(fin)
105+
f= ApproxFunBase.setcanonicaldomain(fin)
106+
xmin, xmax, opfxmin, opfxmax, opmax =
107+
ApproxFunBase.specialfunctionnormalizationpoint2($op, $growth, f, T)
108+
S = space(f)
109+
B=[Evaluation(S,xmin), Evaluation(S,xmax)]
110+
D=Derivative(S)
111+
u=\([B;$L], [opfxmin;opfxmax;$R]; tolerance=10eps(T)*opmax)
112+
113+
setdomain(u, domain(fin))
114+
end
115+
end
116+
end
117+
118+
# ODE gives the first order ODE a special function op satisfies,
119+
# RHS is the right hand side
120+
# growth says what to use to choose a good point to impose an initial condition
121+
for (op, ODE, RHS, growth) in (
122+
(:(erfcx), "D-2f*f'", "-2f'/Base.sqrt(π)", :(real)),
123+
(:(dawson), "D+2f*f'", "f'", :(real))
124+
)
125+
L,R = Meta.parse(ODE), Meta.parse(RHS)
126+
@eval begin
127+
# depice before doing op
128+
$op(f::Fun{<:PiecewiseSpace}) = Fun(map($op, components(f)),PiecewiseSpace)
129+
130+
# We remove the MappedSpace
131+
# function $op{MS<:MappedSpace}(f::Fun{MS})
132+
# g=exp(Fun(f.coefficients,space(f).space))
133+
# Fun(g.coefficients,MappedSpace(domain(f),space(g)))
134+
# end
135+
function $op(fin::Fun)
136+
f=ApproxFunBase.setcanonicaldomain(fin) # removes possible issues with roots
137+
138+
xmax, opfxmax, opmax =
139+
ApproxFunBase.specialfunctionnormalizationpoint($op,$growth,f)
140+
# we will assume the result should be smooth on the domain
141+
# even if f is not
142+
# This supports Line/Rays
143+
D=Derivative(domain(f))
144+
B=Evaluation(domainspace(D),xmax)
145+
u=\([B, $L], [opfxmax, $R]; tolerance=eps(ApproxFunBase.cfstype(fin))*opmax)
146+
147+
setdomain(u, domain(fin))
148+
end
149+
end
150+
end
151+
152+
end

‎src/ApproxFunBase.jl

Lines changed: 1 addition & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@ using IntervalSets
1616
using LinearAlgebra
1717
using LowRankMatrices
1818
using SparseArrays
19-
using SpecialFunctions
19+
# using SpecialFunctions
2020
using StaticArrays: SVector, @SArray, SArray
2121
import Statistics: mean
2222

@@ -60,22 +60,6 @@ import LinearAlgebra: BlasInt, BlasFloat, norm, ldiv!, mul!, det, cross,
6060

6161
import SparseArrays: blockdiag
6262

63-
# import Arpack: eigs
64-
65-
# we need to import all special functions to use Calculus.symbolic_derivatives_1arg
66-
# we can't do importall Base as we replace some Base definitions
67-
import SpecialFunctions: airy, besselh,
68-
lfact, beta, lbeta,
69-
eta, zeta, polygamma, logabsgamma, loggamma,
70-
besselj, bessely, besseli, besselk, besselkx,
71-
hankelh1, hankelh2, hankelh1x, hankelh2x,
72-
# functions from Calculus.symbolic_derivatives_1arg
73-
erf, erfinv, erfc, erfcinv, erfi, gamma, lgamma,
74-
digamma, invdigamma, trigamma,
75-
airyai, airybi, airyaiprime, airybiprime,
76-
besselj0, besselj1, bessely0, bessely1,
77-
erfcx, dawson
78-
7963
import BandedMatrices: bandrange, inbands_setindex!, bandwidth,
8064
colstart, colstop, colrange, rowstart, rowstop, rowrange,
8165
bandwidths, _BandedMatrix, BandedMatrix, isbanded

‎src/specialfunctions.jl

Lines changed: 14 additions & 70 deletions
Original file line numberDiff line numberDiff line change
@@ -164,12 +164,12 @@ end
164164
# ODE gives the first order ODE a special function op satisfies,
165165
# RHS is the right hand side
166166
# growth says what to use to choose a good point to impose an initial condition
167-
for (op, ODE, RHS, growth) in ((:(exp), "D-f'", "0", :(real)),
167+
for (op, ODE, RHS, growth) in (
168+
(:(exp), "D-f'", "0", :(real)),
168169
(:(asinh), "sqrt(f^2+1)*D", "f'", :(real)),
169170
(:(acosh), "sqrt(f^2-1)*D", "f'", :(real)),
170171
(:(atanh), "(1-f^2)*D", "f'", :(real)),
171-
(:(erfcx), "D-2f*f'", "-2f'/sqrt(π)", :(real)),
172-
(:(dawson), "D+2f*f'", "f'", :(real)))
172+
)
173173
L,R = Meta.parse(ODE), Meta.parse(RHS)
174174
@eval begin
175175
# depice before doing op
@@ -208,16 +208,12 @@ function specialfunctionnormalizationpoint2(op, growth, f, T = cfstype(f))
208208
xmin, xmax, opfxmin, opfxmax, opmax
209209
end
210210

211-
for (op,ODE,RHS,growth) in ((:(erf),"f'*D^2+(2f*f'^2-f'')*D","0",:(imag)),
212-
(:(erfi),"f'*D^2-(2f*f'^2+f'')*D","0",:(real)),
211+
for (op,ODE,RHS,growth) in (
213212
(:(sin),"f'*D^2-f''*D+f'^3","0",:(imag)),
214213
(:(cos),"f'*D^2-f''*D+f'^3","0",:(imag)),
215214
(:(sinh),"f'*D^2-f''*D-f'^3","0",:(real)),
216215
(:(cosh),"f'*D^2-f''*D-f'^3","0",:(real)),
217-
(:(airyai),"f'*D^2-f''*D-f*f'^3","0",:(imag)),
218-
(:(airybi),"f'*D^2-f''*D-f*f'^3","0",:(imag)),
219-
(:(airyaiprime),"f'*D^2-f''*D-f*f'^3","airyai(f)*f'^3",:(imag)),
220-
(:(airybiprime),"f'*D^2-f''*D-f*f'^3","airybi(f)*f'^3",:(imag)))
216+
)
221217
L,R = Meta.parse(ODE),Meta.parse(RHS)
222218
@eval begin
223219
function $op(fin::Fun)
@@ -234,8 +230,6 @@ for (op,ODE,RHS,growth) in ((:(erf),"f'*D^2+(2f*f'^2-f'')*D","0",:(imag)),
234230
end
235231
end
236232

237-
erfc(f::Fun) = 1-erf(f)
238-
239233

240234
exp2(f::Fun) = exp(log(2)*f)
241235
exp10(f::Fun) = exp(log(10)*f)
@@ -277,35 +271,8 @@ end
277271
sinpi(f::Fun) = sin(π*f)
278272
cospi(f::Fun) = cos(π*f)
279273

280-
function airy(k::Number,f::Fun)
281-
if k == 0
282-
airyai(f)
283-
elseif k == 1
284-
airyaiprime(f)
285-
elseif k == 2
286-
airybi(f)
287-
elseif k == 3
288-
airybiprime(f)
289-
else
290-
error("invalid argument")
291-
end
292-
end
293-
294-
besselh(ν,k::Integer,f::Fun) = k == 1 ? hankelh1(ν,f) : k == 2 ? hankelh2(ν,f) : throw(Base.Math.AmosException(1))
295-
296-
for jy in (:j, :y)
297-
bjy = Symbol(:bessel, jy)
298-
for ν in (0,1)
299-
bjynu = Symbol(bjy, ν)
300-
@eval SpecialFunctions.$bjynu(f::Fun) = $bjy($ν,f)
301-
end
302-
end
303-
304274
## Miscellaneous
305-
for op in (:(expm1),:(log1p),:(lfact),:(sinc),:(cosc),
306-
:(erfinv),:(erfcinv),:(beta),:(lbeta),
307-
:(eta),:(zeta),:(gamma),:(lgamma),
308-
:(polygamma),:(invdigamma),:(digamma),:(trigamma))
275+
for op in (:(expm1),:(log1p),:(sinc),:(cosc))
309276
@eval begin
310277
$op(f::Fun) = Fun($op ∘ f,domain(f))
311278
end
@@ -457,41 +424,18 @@ for (funsym, exp) in Calculus.symbolic_derivatives_1arg()
457424
funsym == :sign && continue
458425
funsym == :exp && continue
459426
funsym == :sqrt && continue
460-
@eval begin
461-
$(funsym)(z::Fun{<:ConstantSpace,<:Real}) =
462-
Fun($(funsym)(Number(z)),space(z))
463-
$(funsym)(z::Fun{<:ConstantSpace,<:Complex}) =
464-
Fun($(funsym)(Number(z)),space(z))
465-
$(funsym)(z::Fun{<:ConstantSpace}) =
466-
Fun($(funsym)(Number(z)),space(z))
467-
end
468-
end
469-
470-
# Other special functions
471-
for f in [:logabsgamma]
472-
@eval function $f(z::Fun{<:ConstantSpace, <:Real})
473-
t = $f(Number(z))
474-
Fun(t[1], space(z)), t[2]
475-
end
476-
end
477-
function loggamma(z::Fun{<:ConstantSpace})
478-
t = loggamma(Number(z))
479-
Fun(t, space(z))
480-
end
481-
for f in [:gamma, :loggamma]
482-
@eval begin
483-
function $f(a, z::Fun{<:ConstantSpace})
484-
t = $f(a, Number(z))
485-
Fun(t, space(z))
427+
if isdefined(Base, funsym)
428+
@eval begin
429+
$(funsym)(z::Fun{<:ConstantSpace,<:Real}) =
430+
Fun($(funsym)(Number(z)),space(z))
431+
$(funsym)(z::Fun{<:ConstantSpace,<:Complex}) =
432+
Fun($(funsym)(Number(z)),space(z))
433+
$(funsym)(z::Fun{<:ConstantSpace}) =
434+
Fun($(funsym)(Number(z)),space(z))
486435
end
487436
end
488437
end
489438

490-
for f in [:besselj, :besselk, :besselkx, :bessely, :besseli,
491-
:hankelh1x, :hankelh2x, :hankelh1, :hankelh2]
492-
@eval $f(nu, x::Fun{<:ConstantSpace}) = Fun($f(nu, Number(x)), space(x))
493-
end
494-
495439
# Roots
496440

497441
for op in (:(argmax),:(argmin))

0 commit comments

Comments
 (0)