mauro3 / WhereTheWaterFlows.jl

Hydrolocial water flow routing on digital elevation models
MIT License
16 stars 3 forks source link

Don't use recursion #13

Open mauro3 opened 3 years ago

mauro3 commented 3 years ago

There are several recursive algorithms in this package. This means that they could potentially run into Stackoverflows if the catchments and/or streams get too big. I haven't run into this yet (with max stream-length of 4k cells and max catchment size of 5e6 cells). Recursion can be re-written in terms of loops, see e.g. https://raganwald.com/2018/05/27/tail.html.

mauro3 commented 3 years ago

Here pathological example which hits a StackOverflow

ulia> using WhereTheWaterFlows, PyPlot

julia> const WWF = WhereTheWaterFlows
WhereTheWaterFlows

julia> function ramp(l,w)
       y=-w:w
       x=1:l
       dem = (1 .+ y.^2) .* x'.*0.00001
       return x,y,dem
       end
ramp (generic function with 1 method)

julia> x,y,dem = ramp(3*10^4,3); @time area, slen, dir, nout, nin, pits, c, bnds = waterflows(dem, drain_pits=false);
ERROR: TaskFailedException:
StackOverflowError:
Stacktrace:
 [1] _flowrouting_catchments!(::Array{Float64,2}, ::Array{Int64,2}, ::Array{Int64,2}, ::Array{Int8,2}, ::Array{Float64,2}, ::Int64, ::CartesianIndex{2}) at /home/mauro/.julia/dev/WhereTheWaterFlows/src/WhereTheWaterFlows.jl:259

This starts to happen from a stream length of about 13000 cells (or about 25000 cells without Threads.@threads in _flowrouting_catchments) for my computer. Note that the determining factor is maximum(slen) and not the catchment-size, e.g.

julia> x,y,dem = ramp(1.2*10^4,3); @time area, slen, dir, nout, nin, pits, c, bnds = waterflows(dem, drain_pits=false);
  0.025947 seconds (48.09 k allocations: 8.585 MiB)

runs just as well as the 1000x larger catchment

julia> x,y,dem = ramp(1.2*10^4,3000); @time area, slen, dir, nout, nin, pits, c, bnds = waterflows(dem, drain_pits=false); maximum(slen)
 17.063569 seconds (123.87 k allocations: 2.699 GiB, 1.92% gc time)
11999

This is because the recursion is depth-first, I think.

mauro3 commented 4 months ago

Note that on Linux ulimit -s 32768 would quadruple the usual stack size. However, that only seems to apply to the main thread, thus multi-threaded needs to be turned off for this to work.

mauro3 commented 4 months ago

This is semi-addressed with https://github.com/mauro3/WhereTheWaterFlows.jl/commit/2f3776b2e101c475a0e869627f0ed4bd55a9a691