Problem with Broadcasting

[ssfrr]

I'm having some trouble minimizing a function that includes complex numbers and broadcasting. I'm not sure which is causing the problem. On slack @dpsanders mentioned being a little concerned by complex numbers, but looking at the stack trace I'm actually wondering if its a broadcasting problem (it seems to be trying to assign an interval into the result array from a broadcasting operation, and for some reason the result array has float elements.

I'm coming from a DSP background, and the application is delay estimation between two signals where the delay is not an integer number of samples. The main insight here is that a delay in the time domain is the same as multiplying by a complex exponential in the frequency domain. The frequency of the exponential determines the delay, so we can estimate the delay by finding the frequency that minimizes the distance between the two spectra. Here's a MWE:

using Plots
using IntervalOptimisation
using IntervalArithmetic
using FFTW

# first we generate some test data
D = 0.68 # just picking a random target delay (in samples)
x1 = randn(1000)
N = length(x1) ÷ 2 + 1
ω = range(0, π, length=N)
# create a target signal by delaying our original signal and adding some noise
x2 = irfft(rfft(x1) .* exp.(-im .* ω .* D), length(x1)) .+ randn.() .* 0.1
plot([x1[1:40] x2[1:40]])

Which gives

We're trying to estimate D. Here's our error function (which we see has a minimum at 0.68):

X1 = rfft(x1)
X2 = rfft(x2)
err(D) = sum(@.(abs2(X2 - X1 * exp(-im*ω*D))))
x = range(-5,5,length=100)
plot(x, err.(x))

Trying to minimise it gives:

globmin, minxs = IntervalOptimisation.minimise(err, -2..2)

MethodError: no method matching Float64(::Interval{Float64})
Closest candidates are:
  Float64(::Real, !Matched::RoundingMode) where T<:AbstractFloat at rounding.jl:194
  Float64(::T<:Number) where T<:Number at boot.jl:741
  Float64(!Matched::Int8) at float.jl:60
  ...

Stacktrace:
 [1] convert(::Type{Float64}, ::Interval{Float64}) at ./number.jl:7
 [2] setindex!(::Array{Float64,1}, ::Interval{Float64}, ::Int64) at ./array.jl:767
 [3] macro expansion at ./broadcast.jl:843 [inlined]
 [4] macro expansion at ./simdloop.jl:73 [inlined]
 [5] copyto! at ./broadcast.jl:842 [inlined]
 [6] copyto! at ./broadcast.jl:797 [inlined]
 [7] copy at ./broadcast.jl:773 [inlined]
 [8] materialize(::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Nothing,typeof(abs2),Tuple{Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Nothing,typeof(-),Tuple{Array{Complex{Float64},1},Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Nothing,typeof(*),Tuple{Array{Complex{Float64},1},Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Nothing,typeof(exp),Tuple{Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1},Nothing,typeof(*),Tuple{Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0},Nothing,typeof(-),Tuple{Complex{Bool}}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Interval{Float64}}}}}}}}}}}) at ./broadcast.jl:753
 [9] err(::Interval{Float64}) at ./In[54]:3
 [10] #minimise#1(::Type, ::Float64, ::Function, ::typeof(err), ::Interval{Float64}) at /home/sfr/.julia/packages/IntervalOptimisation/A3NrJ/src/optimise.jl:34
 [11] minimise(::Function, ::Interval{Float64}) at /home/sfr/.julia/packages/IntervalOptimisation/A3NrJ/src/optimise.jl:19
 [12] top-level scope at In[55]:1

[ssfrr]

update - rewriting the broadcast in terms of map seems to work:

function err(D)
    errs = map(X1, X2, ω) do x1, x2, om
        abs2(x2 - x1 * exp(-im*om*D))
    end
    sum(errs)
end

So it seems the problem is with the broadcasting

[dpsanders]

(Oops, wrong button sorry.)

It works for me:

julia> err(D) = sum(@.(abs2(X2 - X1 * exp(-im*ω*D))))
err (generic function with 1 method)

julia> @time minimise(err, -5..5)
  0.034151 seconds (3.65 k allocations: 1.380 MiB)
([5366.36, 5405.64], Interval{Float64}[[0.679316, 0.679927], [0.679926, 0.680537], [0.678715, 0.679317], [0.680536, 0.681157], [0.678105, 0.678716], [0.677504, 0.678106], [0.681156, 0.681748], [0.681747, 0.682349], [0.676904, 0.677505], [0.682348, 0.682949], [0.676313, 0.676905], [0.682948, 0.683559], [0.675702, 0.676314], [0.683558, 0.68416], [0.675102, 0.675703]])

What version of IntervalArithmetic.jl et al are you using?

[ssfrr]

I just re-tried this with IntervalOptimisation#master and IntervalArithmetic#master and got the same error. I tried on both Julia 1.1 and Julia 1.2.

[ssfrr]

Also, there's a much more minimal MWE:

using IntervalArithmetic

randn(10) .+ (-2..2)

So this is actually an IntervalArithmetic problem, not IntervalOptimization. I can move the issue over there.

[ssfrr]

Oh, actually this is a dup of https://github.com/JuliaIntervals/IntervalArithmetic.jl/issues/255

[dpsanders]

Closing since needs to be fixed in IntervalArithmetic.jl.

[dpsanders]

Reopening since it's not fixed yet...