Skip to content

Commit c72cbc0

Browse files
authored
improve inference in tolerance in exp (#152)
1 parent c557301 commit c72cbc0

File tree

2 files changed

+42
-24
lines changed

2 files changed

+42
-24
lines changed

src/specialfunctions.jl

Lines changed: 36 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -119,7 +119,7 @@ end
119119

120120

121121
# Default is just try solving ODE
122-
function ^(f::Fun{S,T},β) where {S,T}
122+
function ^(f::Fun, β)
123123
A=Derivative()-β*differentiate(f)/f
124124
B=Evaluation(leftendpoint(domain(f)))
125125
[B;A]\[first(f)^β;0]
@@ -149,34 +149,38 @@ atan(f::Fun)=cumsum(f'/(1+f^2))+atan(first(f))
149149
# condition in calculating secial functions
150150
function specialfunctionnormalizationpoint(op,growth,f)
151151
g=chop(growth(f),eps(cfstype(f)))
152-
xmin = isempty(g.coefficients) ? leftendpoint(domain(g)) : argmin(g)
153-
xmax = isempty(g.coefficients) ? rightendpoint(domain(g)) : argmax(g)
152+
d = domain(g)
153+
T = eltype(d)
154+
xmin = isempty(g.coefficients) ? leftendpoint(d) : T(argmin(g))::T
155+
xmax = isempty(g.coefficients) ? rightendpoint(d) : T(argmax(g))::T
154156
opfxmin,opfxmax = op(f(xmin)),op(f(xmax))
155157
opmax = maximum(abs,(opfxmin,opfxmax))
156-
if abs(opfxmin) == opmax xmax,opfxmax = xmin,opfxmin end
158+
if abs(opfxmin) == opmax
159+
xmax,opfxmax = xmin,opfxmin
160+
end
157161
xmax,opfxmax,opmax
158162
end
159163

160164
# ODE gives the first order ODE a special function op satisfies,
161165
# RHS is the right hand side
162166
# growth says what to use to choose a good point to impose an initial condition
163-
for (op,ODE,RHS,growth) in ((:(exp),"D-f'","0",:(real)),
164-
(:(asinh),"sqrt(f^2+1)*D","f'",:(real)),
165-
(:(acosh),"sqrt(f^2-1)*D","f'",:(real)),
166-
(:(atanh),"(1-f^2)*D","f'",:(real)),
167-
(:(erfcx),"D-2f*f'","-2f'/sqrt(π)",:(real)),
168-
(:(dawson),"D+2f*f'","f'",:(real)))
169-
L,R = Meta.parse(ODE),Meta.parse(RHS)
167+
for (op, ODE, RHS, growth) in ((:(exp), "D-f'", "0", :(real)),
168+
(:(asinh), "sqrt(f^2+1)*D", "f'", :(real)),
169+
(:(acosh), "sqrt(f^2-1)*D", "f'", :(real)),
170+
(:(atanh), "(1-f^2)*D", "f'", :(real)),
171+
(:(erfcx), "D-2f*f'", "-2f'/sqrt(π)", :(real)),
172+
(:(dawson), "D+2f*f'", "f'", :(real)))
173+
L,R = Meta.parse(ODE), Meta.parse(RHS)
170174
@eval begin
171175
# depice before doing op
172-
$op(f::Fun{<:PiecewiseSpace}) = Fun(map(f->$op(f),components(f)),PiecewiseSpace)
176+
$op(f::Fun{<:PiecewiseSpace}) = Fun(map($op, components(f)),PiecewiseSpace)
173177

174178
# We remove the MappedSpace
175179
# function $op{MS<:MappedSpace}(f::Fun{MS})
176180
# g=exp(Fun(f.coefficients,space(f).space))
177181
# Fun(g.coefficients,MappedSpace(domain(f),space(g)))
178182
# end
179-
function $op(fin::Fun{S,T}) where {S,T}
183+
function $op(fin::Fun)
180184
f=setcanonicaldomain(fin) # removes possible issues with roots
181185

182186
xmax,opfxmax,opmax=specialfunctionnormalizationpoint($op,$growth,f)
@@ -185,7 +189,7 @@ for (op,ODE,RHS,growth) in ((:(exp),"D-f'","0",:(real)),
185189
# This supports Line/Rays
186190
D=Derivative(domain(f))
187191
B=Evaluation(domainspace(D),xmax)
188-
u=\([B,eval($L)],Any[opfxmax,eval($R)];tolerance=eps(T)*opmax)
192+
u=\([B, $L], Any[opfxmax, $R]; tolerance=eps(cfstype(fin))*opmax)
189193

190194
setdomain(u,domain(fin))
191195
end
@@ -451,11 +455,11 @@ for (funsym, exp) in Calculus.symbolic_derivatives_1arg()
451455
funsym == :exp && continue
452456
funsym == :sqrt && continue
453457
@eval begin
454-
$(funsym)(z::Fun{CS,T}) where {CS<:ConstantSpace,T<:Real} =
458+
$(funsym)(z::Fun{<:ConstantSpace,<:Real}) =
455459
Fun($(funsym)(Number(z)),space(z))
456-
$(funsym)(z::Fun{CS,T}) where {CS<:ConstantSpace,T<:Complex} =
460+
$(funsym)(z::Fun{<:ConstantSpace,<:Complex}) =
457461
Fun($(funsym)(Number(z)),space(z))
458-
$(funsym)(z::Fun{CS}) where {CS<:ConstantSpace} =
462+
$(funsym)(z::Fun{<:ConstantSpace}) =
459463
Fun($(funsym)(Number(z)),space(z))
460464
end
461465
end
@@ -464,23 +468,31 @@ end
464468

465469
for op in (:(argmax),:(argmin))
466470
@eval begin
467-
function $op(f::Fun{S,T}) where {S<:RealSpace,T<:Real}
471+
function $op(f::Fun{<:RealSpace,<:Real})
468472
# need to check for zero as extremal_args is not defined otherwise
473+
d = domain(f)
474+
T = eltype(d)
469475
iszero(f) && return leftendpoint(domain(f))
470476
# the following avoids warning when differentiate(f)==0
471477
pts = extremal_args(f)
472478
# the extra real avoids issues with complex round-off
473-
pts[$op(real(f.(pts)))]
479+
v = map(realf, pts)::Vector
480+
x = pts[convert(Int, $op(v))::Int]
481+
convert(T, x)::T
474482
end
475483

476-
function $op(f::Fun{S,T}) where {S,T}
484+
function $op(f::Fun)
477485
# need to check for zero as extremal_args is not defined otherwise
486+
d = domain(f)
487+
T = eltype(d)
478488
iszero(f) && return leftendpoint(domain(f))
479489
# the following avoids warning when differentiate(f)==0
480490
pts = extremal_args(f)
481-
fp=f.(pts)
491+
fp = map(f, pts)
482492
@assert norm(imag(fp))<100eps()
483-
pts[$op(real(fp))]
493+
v = real(fp)::Vector
494+
x = pts[convert(Int, $op(v))::Int]
495+
convert(T, x)::T
484496
end
485497
end
486498
end
@@ -512,8 +524,8 @@ end
512524
for op in (:(maximum),:(minimum),:(extrema))
513525
@eval function $op(f::Fun{S,T}) where {S<:RealSpace,T<:Real}
514526
pts = iszero(f') ? [leftendpoint(domain(f))] : extremal_args(f)
515-
516-
$op(f.(pts))
527+
v = map(f, pts)
528+
$op(v)
517529
end
518530
end
519531

test/runtests.jl

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -215,5 +215,11 @@ end
215215
end
216216
end
217217

218+
@testset "misc" begin
219+
a = @inferred ApproxFunBase.specialfunctionnormalizationpoint(exp,real,Fun())
220+
@test a[1] == 1
221+
@test a[2] exp(1)
222+
end
223+
218224
@time include("ETDRK4Test.jl")
219225
include("show.jl")

0 commit comments

Comments
 (0)