Skip to content

Commit

Permalink
Improve decay heuristic
Browse files Browse the repository at this point in the history
  • Loading branch information
@OlivierHnt
OlivierHnt committed Oct 11, 2024
1 parent 41cbc4d commit 7224d61
Show file tree
Hide file tree
Showing 2 changed files with 19 additions and 54 deletions.
49 changes: 8 additions & 41 deletions src/sequence_spaces/arithmetic/elementary.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@ function Base.:^(a::Sequence{<:SequenceSpace}, p::Real)
A = fft(a, fft_size(space_c))
C = A .^ p
c = _call_ifft!(C, space_c, eltype(a))
ν̄ = (_geometric_rate(space(a), coefficients(a))[1] .+ 1) ./ 2
ν̄ = 0.5 .* (max.(_geometric_rate(space(a), coefficients(a))[1], _geometric_rate(space(c), coefficients(c))[1]) .- 1.0) .+ 1.0
_resolve_saturation!(x -> ^(x, p), c, a, ν̄)
return c
end
Expand All @@ -93,7 +93,7 @@ for f ∈ (:exp, :cos, :sin, :cosh, :sinh)
A = fft(a, fft_size(space_c))
C = $f.(A)
c = _call_ifft!(C, space_c, eltype(a))
ν̄ = (_geometric_rate(space(a), coefficients(a))[1] .+ 1) ./ 2
ν̄ = 0.5 .* (max.(_geometric_rate(space(a), coefficients(a))[1], _geometric_rate(space(c), coefficients(c))[1]) .- 1.0) .+ 1.0
_resolve_saturation!($f, c, a, ν̄)
return c
end
Expand Down Expand Up @@ -121,12 +121,14 @@ function _resolve_saturation!(f, c, a, ν)
ν⁻¹ = inv(ν)
C = max(_contour(f, a, ν), _contour(f, a, ν⁻¹))
CoefType = eltype(c)
min_ord = order(c)
for k ∈ indices(space(c))
if mag(c[k]) > mag(C / ν ^ abs(k))
min_ord = min(min_ord, abs(k))
c[k] = zero(CoefType)
end
end
return c
return c, min_ord
end

function _resolve_saturation!(f, c, a, ν::NTuple{N}) where {N}
Expand All @@ -135,12 +137,14 @@ function _resolve_saturation!(f, c, a, ν::NTuple{N}) where {N}
_mix_ = Iterators.product(ntuple(i -> getindex.(_tuple_, i), Val(N))...)
C = maximum(μ -> _contour(f, a, μ), _mix_)
CoefType = eltype(c)
min_ord = order(c)
for k ∈ indices(space(c))
if mag(c[k]) > mag(C / prod(ν .^ abs.(k)))
min_ord = min.(min_ord, abs.(k))
c[k] = zero(CoefType)
end
end
return c
return c, min_ord
end


Expand Down Expand Up @@ -206,40 +210,3 @@ function _boxes!(C, ν, ::Val{D}) where {D}
end
return C
end





# function _find_decay(f, a, approx, ν::NTuple{N}) where {N}
# T = IntervalArithmetic.numtype(eltype(approx))
# N_v = order(approx)
# N_fft = fft_size(space(approx))
# ν̄_max = _rate(geometricweight(a)) # _rate(geometricweight(approx))

# #

# mid_a = mid.(a) # unnecessary allocation

# #

# ν̄_opt = ν̄_max

# C_opt, m_opt = _error(f, mid_a, approx, ν̄_opt, ν, N_fft, N_v, T)
# error_opt = C_opt * m_opt

# for i ∈ 1:N
# for ν̄ᵢ_ ∈ LinRange(ν̄_max[i], ν[i], 10)
# ν̄_ = ntuple(j -> ifelse(j == i, ν̄ᵢ_, ν̄_opt[j]), Val(N))
# C_, m_ = _error(f, mid_a, approx, ν̄_, ν, N_fft, N_v, T)
# error_ = C_ * m_

# error_ > error_opt && break

# ν̄_opt = ν̄_
# error_opt = error_
# end
# end

# return ν̄_opt
# end
24 changes: 11 additions & 13 deletions src/sequence_spaces/validated_sequence.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -367,15 +367,15 @@ function Base.:^(a::ValidatedSequence, p::Real)
@assert isthinzero(sequence_error(a)) # TODO: lift restriction
seq_a = sequence(a)
_isconstant(seq_a) && return ValidatedSequence(_at_value(x -> ^(x, p), seq_a), banachspace(a))
ν_ = rate.(weight(banachspace(a)))
ν_ = rate(banachspace(a))

#

space_c = _image_trunc(^, space(a))
A = fft(seq_a, fft_size(space_c))
C = A .^ p
c = _call_ifft!(C, space_c, eltype(a))
ν̄_ = (_geometric_rate(space(a), coefficients(a))[1] .+ 1) ./ 2
ν̄_ = interval.(0.5 .* (max.(_geometric_rate(space(a), coefficients(a))[1], _geometric_rate(space(c), coefficients(c))[1]) .- 1.0) .+ 1.0)
_resolve_saturation!(x -> ^(x, p), c, a, ν̄_)

#
Expand Down Expand Up @@ -417,20 +417,20 @@ for f ∈ (:exp, :cos, :sin, :cosh, :sinh)
function Base.$f(a::ValidatedSequence)
seq_a = sequence(a)
_isconstant(seq_a) & isthinzero(sequence_error(a)) && return ValidatedSequence(_at_value($f, seq_a), banachspace(a))
ν_ = rate.(weight(banachspace(a)))
ν_ = rate(banachspace(a))

#

space_c = _image_trunc($f, space(a))
A = fft(seq_a, fft_size(space_c))
C = $f.(A)
c = _call_ifft!(C, space_c, eltype(a))
ν̄_ = (_geometric_rate(space(a), coefficients(a))[1] .+ 1) ./ 2
_resolve_saturation!($f, c, a, ν̄_)
ν̄_ = interval.(0.5 .* (max.(_geometric_rate(space(a), coefficients(a))[1], _geometric_rate(space(c), coefficients(c))[1]) .- 1.0) .+ 1.0)
_, N_v = _resolve_saturation!($f, c, a, ν̄_)

#

error = prod(_error($f, seq_a, c, ν_, ν̄_))
error = prod(_error($f, seq_a, c, ν_, ν̄_, N_v))

#

Expand Down Expand Up @@ -458,25 +458,23 @@ end

#

function _error(f, a, approx, ν::Interval, ν̄)
function _error(f, a, approx, ν::Interval, ν̄, N_v)
ν̄⁻¹ = inv(ν̄)

C = max(_contour(f, a, ν), _contour(f, a, ν⁻¹))
C = max(_contour(f, a, ν̄), _contour(f, a, ν̄⁻¹))

N_v = order(approx) # TAG: should be the largest |n| such that the coefficients are zero beyond
q = sum(k -> (ν ^ ExactReal(k) + ν⁻¹ ^ ExactReal(k)) * ν ^ ExactReal(abs(k)), -N_v:N_v)
q = sum(k -> (ν̄ ^ ExactReal(k) + ν̄⁻¹ ^ ExactReal(k)) * ν ^ ExactReal(abs(k)), -N_v:N_v)

return C, q / prod(ν̄ ^ ExactReal( fft_size(space(approx)) ) - ExactReal(1)) + ExactReal(2) * ν̄ / (ν̄ - ν) * (ν * ν̄⁻¹) ^ ExactReal(N_v + 1)
end

function _error(f, a, approx, ν::NTuple{N,Interval}, ν̄) where {N}
function _error(f, a, approx, ν::NTuple{N,Interval}, ν̄, N_v) where {N}
ν̄⁻¹ = inv.(ν̄)
_tuple_ = tuple(ν̄, ν̄⁻¹)
_mix_ = Iterators.product(ntuple(i -> getindex.(_tuple_, i), Val(N))...)
C = maximum(μ -> _contour(f, a, μ), _mix_)

N_v = order(approx) # TAG: error, should be the largest |n| such that the coefficients are zero beyond
q = sum(k -> sum(μ -> prod(μ .^ ExactReal.(k)), _mix_) * prod(ν̄ .^ ExactReal.(abs.(k))), TensorIndices(ntuple(i -> -N_v[i]:N_v[i], Val(N))))
q = sum(k -> sum(μ -> prod(μ .^ ExactReal.(k)), _mix_) * prod(ν .^ ExactReal.(abs.(k))), TensorIndices(ntuple(i -> -N_v[i]:N_v[i], Val(N))))

return C, q / prod(ν̄ .^ ExactReal.( fft_size(space(approx)) ) .- ExactReal(1)) + ExactReal(2^N) * prod(ν̄ ./ (ν̄ .- ν) .* (ν .* ν̄⁻¹) .^ ExactReal.(N_v .+ 1))
end
Expand Down

0 comments on commit 7224d61

Please sign in to comment.