Unexpected time while using the second mode for basins_of_attraction

[awage]

The computation of basins with the second mode can take a long time depending on the radius parameter ε.

res = 600
ds = Systems.henon(zeros(2); a = 1.4, b = 0.3)
xg = yg = range(-2.,2.,length = res)
grid = (xg,yg)
@time basins, att = basins_of_attraction(grid, ds; show_progress = false)
@time basins, att = basins_of_attraction(grid, ds; attractors = att, show_progress = false )
@time basins, att = basins_of_attraction(grid, ds; attractors = att, ε = 1e-6, show_progress = false )

Results are:

@time basins, att = basins_of_attraction(grid, ds; show_progress = false)
 7.703854 seconds (43.50 M allocations: 1.376 GiB, 4.25% gc time, 62.09% compilation time)

@time basins, att = basins_of_attraction(grid, ds; attractors = att, show_progress = false )
123.987298 seconds (286.80 M allocations: 19.574 GiB, 2.83% gc time, 1.35% compilation time)

@time basins, att = basins_of_attraction(grid, ds; attractors = att, ε = 1e-6, show_progress = false )
  1.055172 seconds (4.02 M allocations: 180.549 MiB, 4.70% gc time, 54.13% compilation time)

Default radius is 1e-3. My guess is that the nearest neighbor search finds a lot of points for this radius and there are too much memory allocation in basins_low_level.jl :

        for (k, t) in bsn_nfo.search_trees # this is a `Dict`
            idxs = isearch(t, u_full_state, WithinRange(ε))
            if !isempty(idxs)
                nxt_clr = 2*k + 1
                break
            end
        end

Is there a way to just stop after finding one index instead of having all of them?

[Datseris]

Yes. We need to change this codeblock entirely. We should make it always find the first nearest neighbor with NeighborNumber(1) instead of WithinRange(ε)), and then compare the neighbor distance to ε instead of doing isempty(idxs).

[awage]

This code does the trick:

        for (k, t) in bsn_nfo.search_trees # this is a `Dict`
            idx = isearch(t, u_full_state, NeighborNumber(1))
            if norm(t.data[idx][1] - u_full_state) < ε
                nxt_clr = 2*k + 1
                break
            end
        end

Although t.data[idx][1] looks ugly it works.

Results:

julia> @time basins, att = basins_of_attraction(grid, ds; show_progress = false)
  2.984120 seconds (29.77 M allocations: 636.008 MiB, 4.04% gc time, 3.01% compilation time)

julia> @time basins, att = basins_of_attraction(grid, ds; attractors = att, show_progress = false )
  1.943639 seconds (7.90 M allocations: 658.578 MiB, 6.44% gc time, 33.87% compilation time)

julia> @time basins, att = basins_of_attraction(grid, ds; attractors = att, ε = 1e-6, show_progress = false )
  1.887366 seconds (8.33 M allocations: 673.577 MiB, 6.53% gc time, 32.44% compilation time)

Should I open a PR?

[Datseris]

I think you can do:

# initialize
dist = [Inf]
neighborindex = [0]
searchstate = [u_full_state]

# in the actual code / loop:
Neighborhood.NearestNeighbors.knn_point!(tree, u_full_state, false, dist, neighborindex, Neighborhood.alwaysfalse)

and dist is updated in place to contain the distance to the nearest neighbor. This is the most efficient possible vesrtion ever with 0 allocations.

(please test it I may have misstyped something or so)

yes for the PR

[awage]

Seems that NearestNeighbors is not a member of Neighborhood in the scope (at least in the REPL).

It is not exported in Neighborhood. I can add a PR to this other repository.

[Datseris]

Neighborhood.NearestNeighbors doesn't exist? How can this be, it is clearly imported in the module: https://github.com/JuliaNeighbors/Neighborhood.jl/blob/master/src/kdtree.jl#L1 . It is not exported, I know, which is why I preface it with Neighborhood.

[awage]

I don't know, package status: [645ca80c] Neighborhood v0.2.3

[Datseris]

Well, the name is usable nevertheless:

julia> using Neighborhood

julia> Neighborhood.NearestNeighbors
NearestNeighbors

[awage]

Yes it does call the function! I did not know that

I tried this solution and it still takes a very long time. The first solution you gave was lightning quick!

        dist = [Inf]; neighborindex = [0];
        for (k, t) in bsn_nfo.search_trees # this is a `Dict`
            Neighborhood.NearestNeighbors.knn_point!(t, u_full_state, false, dist, neighborindex, Neighborhood.alwaysfalse)
            if dist[1] < ε
                nxt_clr = 2*k + 1
                break
            end
        end

[awage]

There is something strange happening with the default distance 1e-3:

@time basins, att2 = basins_of_attraction(grid, ds; attractors = att, ε = 1e-4, show_progress = true )
  1.474840 seconds (5.96 M allocations: 431.481 MiB, 5.87% gc time)

julia> @time basins, att2 = basins_of_attraction(grid, ds; attractors = att, ε = 1e-2, show_progress = true )
  2.295727 seconds (8.15 M allocations: 631.706 MiB, 5.39% gc time)

@time basins, att2 = basins_of_attraction(grid, ds; attractors = att, ε = 1e-3, show_progress = true )
184.270209 seconds (451.65 M allocations: 40.269 GiB, 4.60% gc time)

But for 1e-3 it takes 3 min! It's confusing.

[Datseris]

you're still re-creating dist = [Inf]; neighborindex = [0]; all the time though at each Δt step, they should only be created once.

[awage]

I get it. I included these variables into the structure BasinsInfo. But still get the same result:

        for (k, t) in bsn_nfo.search_trees # this is a `Dict`
            Neighborhood.NearestNeighbors.knn_point!(t, u_full_state, false, bsn_nfo.dist, bsn_nfo.neighborindex, Neighborhood.alwaysfalse)
            #@show bsn_nfo.dist
            if bsn_nfo.dist[1] < ε
                nxt_clr = 2*k + 1
                break
            end
        end

@time basins, att2 = basins_of_attraction(grid, ds; attractors = att, ε = 1e-3, show_progress = false )
140.511056 seconds (3.37 M allocations: 147.891 MiB, 0.02% gc time, 0.48% compilation time)

Allocations are way down! But the search is not efficient for the radius 1e-3

[Datseris]

its still much faster than current approach so the PR should still be opened. I don't know what the deal with 1e-3.