Bounding box in higher dimensions

[rmcaixeta]

To deal with issue #1148

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 87.70%. Comparing base (86359bb) to head (57c0ffc).

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1149      +/-   ##
==========================================
- Coverage   87.72%   87.70%   -0.03%     
==========================================
  Files         191      191              
  Lines        6013     6000      -13     
==========================================
- Hits         5275     5262      -13     
  Misses        738      738              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[eliascarv]

Code generated by the _pboxes function for some dimensions:

julia> gen(1)
quote
    #= REPL[10]:24 =#
    p = first(points)
    #= REPL[10]:25 =#
    ℒ = lentype(p)
    #= REPL[10]:26 =#
    cmin1 = typemax(ℒ)
    #= REPL[10]:27 =#
    cmax1 = typemin(ℒ)
    #= REPL[10]:29 =#
    for p = points
        #= REPL[10]:30 =#
        c = CoordRefSystems.values(convert(Cartesian, coords(p)))
        #= REPL[10]:31 =#
        cmin1 = min(c[1], cmin1)
        #= REPL[10]:32 =#
        cmax1 = max(c[1], cmax1)
        #= REPL[10]:33 =#
    end
    #= REPL[10]:35 =#
    Box(withcrs(p, (cmin1,)), withcrs(p, (cmax1,)))
end

julia> gen(2)
quote
    #= REPL[10]:24 =#
    p = first(points)
    #= REPL[10]:25 =#
    ℒ = lentype(p)
    #= REPL[10]:26 =#
    cmin1 = typemax(ℒ)
    cmin2 = typemax(ℒ)
    #= REPL[10]:27 =#
    cmax1 = typemin(ℒ)
    cmax2 = typemin(ℒ)
    #= REPL[10]:29 =#
    for p = points
        #= REPL[10]:30 =#
        c = CoordRefSystems.values(convert(Cartesian, coords(p)))
        #= REPL[10]:31 =#
        cmin1 = min(c[1], cmin1)
        cmin2 = min(c[2], cmin2)
        #= REPL[10]:32 =#
        cmax1 = max(c[1], cmax1)
        cmax2 = max(c[2], cmax2)
        #= REPL[10]:33 =#
    end
    #= REPL[10]:35 =#
    Box(withcrs(p, (cmin1, cmin2)), withcrs(p, (cmax1, cmax2)))
end

julia> gen(3)
quote
    #= REPL[10]:24 =#
    p = first(points)
    #= REPL[10]:25 =#
    ℒ = lentype(p)
    #= REPL[10]:26 =#
    cmin1 = typemax(ℒ)
    cmin2 = typemax(ℒ)
    cmin3 = typemax(ℒ)
    #= REPL[10]:27 =#
    cmax1 = typemin(ℒ)
    cmax2 = typemin(ℒ)
    cmax3 = typemin(ℒ)
    #= REPL[10]:29 =#
    for p = points
        #= REPL[10]:30 =#
        c = CoordRefSystems.values(convert(Cartesian, coords(p)))
        #= REPL[10]:31 =#
        cmin1 = min(c[1], cmin1)
        cmin2 = min(c[2], cmin2)
        cmin3 = min(c[3], cmin3)
        #= REPL[10]:32 =#
        cmax1 = max(c[1], cmax1)
        cmax2 = max(c[2], cmax2)
        cmax3 = max(c[3], cmax3)
        #= REPL[10]:33 =#
    end
    #= REPL[10]:35 =#
    Box(withcrs(p, (cmin1, cmin2, cmin3)), withcrs(p, (cmax1, cmax2, cmax3)))
end

julia> gen(4)
quote
    #= REPL[10]:24 =#
    p = first(points)
    #= REPL[10]:25 =#
    ℒ = lentype(p)
    #= REPL[10]:26 =#
    cmin1 = typemax(ℒ)
    cmin2 = typemax(ℒ)
    cmin3 = typemax(ℒ)
    cmin4 = typemax(ℒ)
    #= REPL[10]:27 =#
    cmax1 = typemin(ℒ)
    cmax2 = typemin(ℒ)
    cmax3 = typemin(ℒ)
    cmax4 = typemin(ℒ)
    #= REPL[10]:29 =#
    for p = points
        #= REPL[10]:30 =#
        c = CoordRefSystems.values(convert(Cartesian, coords(p)))
        #= REPL[10]:31 =#
        cmin1 = min(c[1], cmin1)
        cmin2 = min(c[2], cmin2)
        cmin3 = min(c[3], cmin3)
        cmin4 = min(c[4], cmin4)
        #= REPL[10]:32 =#
        cmax1 = max(c[1], cmax1)
        cmax2 = max(c[2], cmax2)
        cmax3 = max(c[3], cmax3)
        cmax4 = max(c[4], cmax4)
        #= REPL[10]:33 =#
    end
    #= REPL[10]:35 =#
    Box(withcrs(p, (cmin1, cmin2, cmin3, cmin4)), withcrs(p, (cmax1, cmax2, cmax3, cmax4)))
end

[eliascarv]

Allocation in heterogeneous collections:
setup

b1 = Box((-3, -1), (0.5, 0.5))
b2 = b1 |> Translate(2, 2)
s = Sphere((0, 0), 2)
m1 = Multi([b1, s])
m2 = Multi([b1, s, b2])

master

julia> @allocated(boundingbox(m1))
2544

julia> @allocated(boundingbox(m2))
3504

PR

julia> @allocated(boundingbox(m1))
2672

julia> @allocated(boundingbox(m2))
3696

[juliohm]

@rmcaixeta , we explored this idea with generated functions and noticed that the compiler allocates more memory when the Multi geometries or GeometrySet are heterogeneous. The type instability leads to dynamic execution and more bytes per geometry.

Since boundingbox is a pretty fundamental function, which is used in various algorithms for optimization, we are a bit reluctant to support arbitrary N > 3 dimensions.

Could you please elaborate on your use case? If you are trying to simulate an abstract Gaussian process over N > 3 non-physical dimensions, you can either add this method locally to see if the rest of GeoStatsProcesses.jl works, or choose a more general package for Gaussian process simulation that is not developed with geospatial coordinates in mind.

[rmcaixeta]

Nice piece of code @eliascarv
Hi @juliohm . The simulation is made in higher dimension (to "remove" some nonstationary features) and mapped back to original space. But that's fine, it apparently works well locally if I just add the initial suggested function. without changing anything in Meshes.jl. I opened the issue because it'd change just a small piece of code, I didn't imagine it'd cause problems. You can cancel the PR then. Have a good conference in Belem soon, cheers

[juliohm]

Thank you! Closing the PR.