LibGEOS.jl
LibGEOS.jl copied to clipboard
JuliaGeo
Crashes when using @threads with intersection functions
Currently, I have two arrays of type Array{Polygon,1}, called grid_1 and grid_2. I would like to find the intersecting area of every point of grid_1 with every point of grid_2 (i.e., yield a matrix with dimensions length(grid_1) by length(grid_2). I am wondering if there is anyway that may be faster that running a two for loops as follows?
polygon_intersection = Matrix{Float64}(undef,length(grid_1),length(grid_2));
for i = 1:length(grid_1)
for j = 1:length(grid_2)
polygon_intersection[i,j] = area(intersection(grid_1[i],grid_2[j]))
end
end
I was thinking of using something of the form area.(intersection.(grid_1[i],grid_2)) and then only have to loop over one dimension, but that does not appear to be implemented (and I am not sure these array calculations would even be faster than a for loop). Any suggestions are appreciated!
In julia there is generally no performance penalty on loops, so I wouldn't worry about getting rid of them.
From https://trac.osgeo.org/geos/wiki/RFC3
Function names in the new API will be updated with an _r, as is the familiar C standard for reentrant/thread safe versions.
GEOSIntersection_r is used here
https://github.com/JuliaGeo/LibGEOS.jl/blob/v0.6.8/src/geos_functions.jl#L423-L429
EDIT: Be careful with using threads.
I am not sure if this is worth opening another issue, but I am having some troubles with these functions being thread safe I believe. When I run the following MWE code (written for this question)
using LibGEOS
grid_1 = Dict();
# c = center, ll = lower left corner, ul = upper left corner
# ur = upper right corner, lr = lower right corner
grid_1["x_c"] = collect(0:0.5:100);
grid_1["y_c"] = collect(0:0.5:100);
grid_1["x_ll"] = grid_1["x_c"] .- (0.5/2);
grid_1["x_ul"] = grid_1["x_c"] .- (0.5/2);
grid_1["x_ur"] = grid_1["x_c"] .+ (0.5/2);
grid_1["x_lr"] = grid_1["x_c"] .+ (0.5/2);
grid_1["y_ll"] = grid_1["y_c"] .- (0.5/2);
grid_1["y_ul"] = grid_1["y_c"] .+ (0.5/2);
grid_1["y_ur"] = grid_1["y_c"] .+ (0.5/2);
grid_1["y_lr"] = grid_1["y_c"] .- (0.5/2);
grid_1["value"] = rand(length(grid_1["x_c"]));
grid_1["uncertainty"] = rand(length(grid_1["x_c"]));
grid_2 = Dict();
grid_2["x_c"] = collect(1:0.5:101);
grid_2["y_c"] = collect(1:0.5:101);
grid_2["x_ll"] = grid_2["x_c"] .- (0.5/2);
grid_2["x_ul"] = grid_2["x_c"] .- (0.5/2);
grid_2["x_ur"] = grid_2["x_c"] .+ (0.5/2);
grid_2["x_lr"] = grid_2["x_c"] .+ (0.5/2);
grid_2["y_ll"] = grid_2["y_c"] .- (0.5/2);
grid_2["y_ul"] = grid_2["y_c"] .+ (0.5/2);
grid_2["y_ur"] = grid_2["y_c"] .+ (0.5/2);
grid_2["y_lr"] = grid_2["y_c"] .- (0.5/2);
grid_2["value"] = rand(length(grid_2["x_c"]));
grid_2["uncertainty"] = rand(length(grid_2["x_c"]));
for grid in (grid_1,grid_2)
grid["polygon"] = Vector{Polygon}(undef,length(grid["x_c"]));
# makes polygons from the pixcorners for both grids
for idx = 1:length(grid["polygon"])
grid["polygon"][idx] = Polygon([[[grid["x_ll"][idx],grid["y_ll"][idx]],[grid["x_ul"][idx],grid["y_ul"][idx]],[grid["x_ur"][idx],grid["y_ur"][idx]],[grid["x_lr"][idx],grid["y_lr"][idx]],[grid["x_ll"][idx],grid["y_ll"][idx]]]]);
end
grid["polygon_area"] = area.(grid["polygon"]);
end
# averages grid_1 value to grid_2
function F_l2_to_l2_regrid(grid_1,grid_2)
grid_2["grid_1_value"] = Vector{Float64}(undef,length(grid_2["polygon"]));
Threads.@threads for j = 1:length(grid_2["polygon"])
A = Vector{Float64}(undef,length(grid_1["polygon"]));
B = Vector{Float64}(undef,length(grid_1["polygon"]));
p = 1;
polygon_intersection = Vector{Float64}(undef,length(grid_1["polygon"]));
for i = 1:length(grid_1["polygon"])
polygon_intersection[i] = area(intersection(grid_1["polygon"][i],grid_2["polygon"][j]));
A[i] = grid_1["value"][i]*polygon_intersection[i]/(grid_1["uncertainty"][i]^p)/grid_2["polygon_area"][j];
B[i] = polygon_intersection[i]/(grid_1["uncertainty"][i]^p)/grid_2["polygon_area"][j];
end
grid_2["grid_1_value"][j] = sum(A)/sum(B);
end
end
F_l2_to_l2_regrid(grid_1,grid_2)
I get a very long error that kills my Julia session, specifically:
*** Error in `julia': double free or corruption (fasttop): 0x00002b760c0cc660 ***
======= Backtrace: =========
/lib64/libc.so.6(+0x81329)[0x2b75f0fa0329]
/n/home06/nbalasus/.julia/artifacts/514b2d8cdf634f84bc9d53a25f76ec2826d964ce/lib/libgeos-3.9.0.so(_ZNK4geos4geom8Geometry19getEnvelopeInternalEv+0x48)[0x2b763b6cd988]
/n/home06/nbalasus/.julia/artifacts/514b2d8cdf634f84bc9d53a25f76ec2826d964ce/lib/libgeos-3.9.0.so(_ZNK4geos4geom7Polygon23computeEnvelopeInternalEv+0x19)[0x2b763b6ddd59]
/n/home06/nbalasus/.julia/artifacts/514b2d8cdf634f84bc9d53a25f76ec2826d964ce/lib/libgeos-3.9.0.so(_ZNK4geos4geom8Geometry19getEnvelopeInternalEv+0x2a)[0x2b763b6cd96a]
/n/home06/nbalasus/.julia/artifacts/514b2d8cdf634f84bc9d53a25f76ec2826d964ce/lib/libgeos-3.9.0.so(_ZN4geos9operation9overlayng11OverlayUtil13isEnvDisjointEPKNS_4geom8GeometryES6_PKNS3_14PrecisionModelE+0x92)[0x2b763b767892]
/n/home06/nbalasus/.julia/artifacts/514b2d8cdf634f84bc9d53a25f76ec2826d964ce/lib/libgeos-3.9.0.so(_ZN4geos9operation9overlayng9OverlayNG9getResultEv+0x46)[0x2b763b764b26]
/n/home06/nbalasus/.julia/artifacts/514b2d8cdf634f84bc9d53a25f76ec2826d964ce/lib/libgeos-3.9.0.so(_ZN4geos9operation9overlayng9OverlayNG7overlayEPKNS_4geom8GeometryES6_iPKNS3_14PrecisionModelE+0x65)[0x2b763b764cd5]
/n/home06/nbalasus/.julia/artifacts/514b2d8cdf634f84bc9d53a25f76ec2826d964ce/lib/libgeos-3.9.0.so(_ZN4geos9operation9overlayng15OverlayNGRobust7OverlayEPKNS_4geom8GeometryES6_i+0xa1)[0x2b763b765e91]
/n/home06/nbalasus/.julia/artifacts/514b2d8cdf634f84bc9d53a25f76ec2826d964ce/lib/libgeos-3.9.0.so(_ZN4geos4geom16HeuristicOverlayEPKNS0_8GeometryES3_i+0x235)[0x2b763b6d6a05]
/n/home06/nbalasus/.julia/artifacts/514b2d8cdf634f84bc9d53a25f76ec2826d964ce/lib/libgeos-3.9.0.so(_ZNK4geos4geom8Geometry12intersectionEPKS1_+0x64)[0x2b763b6cec24]
/n/home06/nbalasus/.julia/artifacts/514b2d8cdf634f84bc9d53a25f76ec2826d964ce/lib/libgeos_c.so(GEOSIntersection_r+0x27)[0x2b763b3c96f7]
[0x2b762c695582]
[0x2b762c695d0e]
[0x2b762c694767]
[0x2b762c69539f]
[0x2b762c6953bd]
/n/sw/helmod/apps/centos7/Comp/gcc/9.3.0-fasrc01/julia/1.5.1-fasrc01/bin/../lib/libjulia.so.1(+0xccc86)[0x2b75f0211c86]
======= Memory map: ========
00400000-00402000 r-xp 00000000 00:2d 4474665034 /n/sw/helmod/apps/centos7/Comp/gcc/9.3.0-fasrc01/julia/1.5.1-fasrc01/bin/julia
00602000-00603000 r--p 00002000 00:2d 4474665034 /n/sw/helmod/apps/centos7/Comp/gcc/9.3.0-fasrc01/julia/1.5.1-fasrc01/bin/julia
00603000-00604000 rw-p 00003000 00:2d 4474665034 /n/sw/helmod/apps/centos7/Comp/gcc/9.3.0-fasrc01/julia/1.5.1-fasrc01/bin/julia
01d89000-029a5000 rw-p 00000000 00:00 0 [heap]
2b75eff21000-2b75eff43000 r-xp 00000000 fd:00 33669052 /usr/lib64/ld-2.17.so
2b75eff43000-2b75eff45000 rw-p 00000000 00:00 0
2b75eff45000-2b75eff48000 r--p 00000000 00:00 0
2b75eff48000-2b75eff58000 rwxp 00000000 00:00 0
2b75eff58000-2b75eff5b000 rw-p 00000000 00:00 0
2b75eff60000-2b75f0139000 rw-p 00000000 00:00 0
2b75f0139000-2b75f0140000 r--s 00000000 fd:00 100816988 /usr/lib64/gconv/gconv-modules.cache
2b75f0142000-2b75f0143000 r--p 00021000 fd:00 33669052 /usr/lib64/ld-2.17.so
2b75f0143000-2b75f0144000 rw-p 00022000 fd:00 33669052 /usr/lib64/ld-2.17.so
2b75f0144000-2b75f0145000 rw-p 00000000 00:00 0
2b75f0145000-2b75f03c3000 r-xp 00000000 00:2d 4450091073 /n/sw/helmod/apps/centos7/Comp/gcc/9.3.0-fasrc01/julia/1.5.1-fasrc01/lib/libjulia.so.1.5
2b75f03c3000-2b75f05c2000 ---p 0027e000 00:2d 4450091073 /n/sw/helmod/apps/centos7/Comp/gcc/9.3.0-fasrc01/julia/1.5.1-fasrc01/lib/libjulia.so.1.5
2b75f05c2000-2b75f05c8000 r--p 0027d000 00:2d 4450091073 /n/sw/helmod/apps/centos7/Comp/gcc/9.3.0-fasrc01/julia/1.5.1-fasrc01/lib/libjulia.so.1.5
2b75f05c8000-2b75f061a000 rw-p 00283000 00:2d 4450091073 /n/sw/helmod/apps/centos7/Comp/gcc/9.3.0-fasrc01/julia/1.5.1-fasrc01/lib/libjulia.so.1.5
2b75f061a000-2b75f08f7000 rw-p 00000000 00:00 0
2b75f08f7000-2b75f08f9000 r-xp 00000000 fd:00 33689620 /usr/lib64/libdl-2.17.so
2b75f08f9000-2b75f0af9000 ---p 00002000 fd:00 33689620 /usr/lib64/libdl-2.17.so
2b75f0af9000-2b75f0afa000 r--p 00002000 fd:00 33689620 /usr/lib64/libdl-2.17.so
2b75f0afa000-2b75f0afb000 rw-p 00003000 fd:00 33689620 /usr/lib64/libdl-2.17.so
2b75f0afb000-2b75f0b02000 r-xp 00000000 fd:00 33872652 /usr/lib64/librt-2.17.so
2b75f0b02000-2b75f0d01000 ---p 00007000 fd:00 33872652 /usr/lib64/librt-2.17.so
2b75f0d01000-2b75f0d02000 r--p 00006000 fd:00 33872652 /usr/lib64/librt-2.17.so
2b75f0d02000-2b75f0d03000 rw-p 00007000 fd:00 33872652 /usr/lib64/librt-2.17.so
2b75f0d03000-2b75f0d1a000 r-xp 00000000 fd:00 33669082 /usr/lib64/libpthread-2.17.so
2b75f0d1a000-2b75f0f19000 ---p 00017000 fd:00 33669082 /usr/lib64/libpthread-2.17.so
2b75f0f19000-2b75f0f1a000 r--p 00016000 fd:00 33669082 /usr/lib64/libpthread-2.17.so
2b75f0f1a000-2b75f0f1b000 rw-p 00017000 fd:00 33669082 /usr/lib64/libpthread-2.17.so
2b75f0f1b000-2b75f0f1f000 rw-p 00000000 00:00 0
2b75f0f1f000-2b75f10e3000 r-xp 00000000 fd:00 33669060 /usr/lib64/libc-2.17.so
2b75f10e3000-2b75f12e2000 ---p 001c4000 fd:00 33669060 /usr/lib64/libc-2.17.so
2b75f12e2000-2b75f12e6000 r--p 001c3000 fd:00 33669060 /usr/lib64/libc-2.17.so
2b75f12e6000-2b75f12e8000 rw-p 001c7000 fd:00 33669060 /usr/lib64/libc-2.17.so
2b75f12e8000-2b75f12ed000 rw-p 00000000 00:00 0
2b75f12ed000-2b75f131b000 r-xp 00000000 00:2d 4474995838 /n/sw/helmod/apps/centos7/Comp/gcc/9.3.0-fasrc01/julia/1.5.1-fasrc01/lib/julia/libunwind.so.8.0.1
2b75f131b000-2b75f151a000 ---p 0002e000 00:2d 4474995838 /n/sw/helmod/apps/centos7/Comp/gcc/9.3.0-fasrc01/julia/1.5.1-fasrc01/lib/julia/libunwind.so.8.0.1
2b75f151a000-2b75f151b000 r--p 0002d000 00:2d 4474995838 /n/sw/helmod/apps/centos7/Comp/gcc/9.3.0-fasrc01/julia/1.5.1-fasrc01/lib/julia/libunwind.so.8.0.1
2b75f151b000-2b75f151c000 rw-p 0002e000 00:2d 4474995838 /n/sw/helmod/apps/centos7/Comp/gcc/9.3.0-fasrc01/julia/1.5.1-fasrc01/lib/julia/libunwind.so.8.0.1
2b75f151c000-2b75f1526000 rw-p 00000000 00:00 0
2b75f1526000-2b75f4746000 r-xp 00000000 00:2d 4431219297 /n/sw/helmod/apps/centos7/Comp/gcc/9.3.0-fasrc01/julia/1.5.1-fasrc01/lib/julia/libLLVM-9jl.so
2b75f4746000-2b75f4946000 ---p 03220000 00:2d 4431219297 /n/sw/helmod/apps/centos7/Comp/gcc/9.3.0-fasrc01/julia/1.5.1-fasrc01/lib/julia/libLLVM-9jl.so
2b75f4946000-2b75f4c1d000 r--p 03220000 00:2d 4431219297 /n/sw/helmod/apps/centos7/Comp/gcc/9.3.0-fasrc01/julia/1.5.1-fasrc01/lib/julia/libLLVM-9jl.so
2b75f4c1d000-2b75f4c44000 rw-p 034f7000 00:2d 4431219297 /n/sw/helmod/apps/centos7/Comp/gcc/9.3.0-fasrc01/julia/1.5.1-fasrc01/lib/julia/libLLVM-9jl.so
2b75f4c44000-2b75f4c89000 rw-p 00000000 00:00 0
2b75f4c89000-2b75f4e52000 r-xp 00000000 00:2d 4475100338 /n/sw/helmod/apps/centos7/Core/gcc/9.3.0-fasrc01/lib64/libstdc++.so.6.0.28
2b75f4e52000-2b75f5052000 ---p 001c9000 00:2d 4475100338 /n/sw/helmod/apps/centos7/Core/gcc/9.3.0-fasrc01/lib64/libstdc++.so.6.0.28
2b75f5052000-2b75f505d000 r--p 001c9000 00:2d 4475100338 /n/sw/helmod/apps/centos7/Core/gcc/9.3.0-fasrc01/lib64/libstdc++.so.6.0.28
2b75f505d000-2b75f5060000 rw-p 001d4000 00:2d 4475100338 /n/sw/helmod/apps/centos7/Core/gcc/9.3.0-fasrc01/lib64/libstdc++.so.6.0.28
2b75f5060000-2b75f5063000 rw-p 00000000 00:00 0
2b75f5063000-2b75f5164000 r-xp 00000000 fd:00 33689622 /usr/lib64/libm-2.17.so
2b75f5164000-2b75f5363000 ---p 00101000 fd:00 33689622 /usr/lib64/libm-2.17.so
2b75f5363000-2b75f5364000 r--p 00100000 fd:00 33689622 /usr/lib64/libm-2.17.so
2b75f5364000-2b75f5365000 rw-p 00101000 fd:00 33689622 /usr/lib64/libm-2.17.so
2b75f5365000-2b75f537b000 r-xp 00000000 00:2d 4475273610 /n/sw/helmod/apps/centos7/Core/gcc/9.3.0-fasrc01/lib64/libgcc_s.so.1
2b75f537b000-2b75f557b000 ---p 00016000 00:2d 4475273610 /n/sw/helmod/apps/centos7/Core/gcc/9.3.0-fasrc01/lib64/libgcc_s.so.1
2b75f557b000-2b75f557c000 r--p 00016000 00:2d 4475273610 /n/sw/helmod/apps/centos7/Core/gcc/9.3.0-fasrc01/lib64/libgcc_s.so.1
2b75f557c000-2b75f557d000 rw-p 00017000 00:2d 4475273610 /n/sw/helmod/apps/centos7/Core/gcc/9.3.0-fasrc01/lib64/libgcc_s.so.1
2b75f557d000-2b75fbabf000 r--p 00000000 fd:00 67214576 /usr/lib/locale/locale-archive
2b75fbabf000-2b75fbac0000 ---p 00000000 00:00 0
2b75fbac0000-2b75fbcc0000 rw-p 00000000 00:00 0
2b75fbcc0000-2b75fbcc1000 ---p 00000000 00:00 0
2b75fbcc1000-2b75fcbcf000 rw-p 00000000 00:00 0
2b75fcbcf000-2b75fd42f000 r-xp 00000000 00:2d 4474897380 /n/sw/helmod/apps/centos7/Comp/gcc/9.3.0-fasrc01/julia/1.5.1-fasrc01/lib/julia/sys.so
2b75fd42f000-2b75fd62e000 ---p 00860000 00:2d 4474897380 /n/sw/helmod/apps/centos7/Comp/gcc/9.3.0-fasrc01/julia/1.5.1-fasrc01/lib/julia/sys.so
2b75fd62e000-2b75fd651000 r--p 0085f000 00:2d 4474897380 /n/sw/helmod/apps/centos7/Comp/gcc/9.3.0-fasrc01/julia/1.5.1-fasrc01/lib/julia/sys.so
2b75fd651000-2b7604e8a000 rw-p 00882000 00:2d 4474897380 /n/sw/helmod/apps/centos7/Comp/gcc/9.3.0-fasrc01/julia/1.5.1-fasrc01/lib/julia/sys.so
2b7604e8a000-2b7604ea7000 rw-p 00000000 00:00 0
2b7604ea7000-2b7608eab000 rw-p 00000000 00:00 0
2b7608eab000-2b7609217000 rw-p 00000000 00:00 0
2b7609217000-2b760921f000 ---p 00000000 00:00 0
2b760921f000-2b760a01a000 rw-p 00000000 00:00 0
2b760a01a000-2b760a01b000 ---p 00000000 00:00 0
2b760a01b000-2b760a21b000 rw-p 00000000 00:00 0
2b760a21b000-2b760a21c000 ---p 00000000 00:00 0
2b760a21c000-2b760a41c000 rw-p 00000000 00:00 0
2b760a41c000-2b760a41d000 ---p 00000000 00:00 0
2b760a41d000-2b760a61d000 rw-p 00000000 00:00 0
2b760a61d000-2b760a61e000 ---p 00000000 00:00 0
2b760a61e000-2b760a81e000 rw-p 00000000 00:00 0
2b760a81e000-2b760a81f000 ---p 00000000 00:00 0
2b760a81f000-2b760aa1f000 rw-p 00000000 00:00 0
2b760aa1f000-2b760aa20000 ---p 00000000 00:00 0
2b760aa20000-2b760ac20000 rw-p 00000000 00:00 0
2b760ac20000-2b760ac21000 ---p 00000000 00:00 0
2b760ac21000-2b760ae21000 rw-p 00000000 00:00 0
2b760ae21000-2b760ae22000 ---p 00000000 00:00 0
2b760ae22000-2b760b622000 rw-p 00000000 00:00 0
2b760b622000-2b760b62a000 ---p 00000000 00:00 0
2b760b62a000-2b760ba22000 rw-p 00000000 00:00 0
2b760ba22000-2b760ba2a000 ---p 00000000 00:00 0
2b760ba2a000-2b760bf22000 rw-p 00000000 00:00 0
2b760c000000-2b760c155000 rw-p 00000000 00:00 0
2b760c155000-2b7610000000 ---p 00000000 00:00 0
signal (6): Aborted
in expression starting at REPL[30]:1
gsignal at /lib64/libc.so.6 (unknown line)
abort at /lib64/libc.so.6 (unknown line)
__libc_message at /lib64/libc.so.6 (unknown line)
_int_free at /lib64/libc.so.6 (unknown line)
_ZNK4geos4geom8Geometry19getEnvelopeInternalEv at /n/home06/nbalasus/.julia/artifacts/514b2d8cdf634f84bc9d53a25f76ec2826d964ce/lib/libgeos-3.9.0.so (unknown line)
_ZNK4geos4geom7Polygon23computeEnvelopeInternalEv at /n/home06/nbalasus/.julia/artifacts/514b2d8cdf634f84bc9d53a25f76ec2826d964ce/lib/libgeos-3.9.0.so (unknown line)
_ZNK4geos4geom8Geometry19getEnvelopeInternalEv at /n/home06/nbalasus/.julia/artifacts/514b2d8cdf634f84bc9d53a25f76ec2826d964ce/lib/libgeos-3.9.0.so (unknown line)
_ZN4geos9operation9overlayng11OverlayUtil13isEnvDisjointEPKNS_4geom8GeometryES6_PKNS3_14PrecisionModelE at /n/home06/nbalasus/.julia/artifacts/514b2d8cdf634f84bc9d53a25f76ec2826d964ce/lib/libgeos-3.9.0.so (unknown line)
_ZN4geos9operation9overlayng9OverlayNG9getResultEv at /n/home06/nbalasus/.julia/artifacts/514b2d8cdf634f84bc9d53a25f76ec2826d964ce/lib/libgeos-3.9.0.so (unknown line)
_ZN4geos9operation9overlayng9OverlayNG7overlayEPKNS_4geom8GeometryES6_iPKNS3_14PrecisionModelE at /n/home06/nbalasus/.julia/artifacts/514b2d8cdf634f84bc9d53a25f76ec2826d964ce/lib/libgeos-3.9.0.so (unknown line)
_ZN4geos9operation9overlayng15OverlayNGRobust7OverlayEPKNS_4geom8GeometryES6_i at /n/home06/nbalasus/.julia/artifacts/514b2d8cdf634f84bc9d53a25f76ec2826d964ce/lib/libgeos-3.9.0.so (unknown line)
_ZN4geos4geom16HeuristicOverlayEPKNS0_8GeometryES3_i at /n/home06/nbalasus/.julia/artifacts/514b2d8cdf634f84bc9d53a25f76ec2826d964ce/lib/libgeos-3.9.0.so (unknown line)
_ZNK4geos4geom8Geometry12intersectionEPKS1_ at /n/home06/nbalasus/.julia/artifacts/514b2d8cdf634f84bc9d53a25f76ec2826d964ce/lib/libgeos-3.9.0.so (unknown line)
GEOSIntersection_r at /n/home06/nbalasus/.julia/artifacts/514b2d8cdf634f84bc9d53a25f76ec2826d964ce/lib/libgeos_c.so (unknown line)
GEOSIntersection_r at /n/home06/nbalasus/.julia/packages/LibGEOS/JGKwE/src/libgeos_api.jl:353 [inlined]
intersection at /n/home06/nbalasus/.julia/packages/LibGEOS/JGKwE/src/geos_functions.jl:424 [inlined]
intersection at /n/home06/nbalasus/.julia/packages/LibGEOS/JGKwE/src/geos_functions.jl:424
intersection at /n/home06/nbalasus/.julia/packages/LibGEOS/JGKwE/src/geos_operations.jl:72
macro expansion at ./REPL[29]:12 [inlined]
#4#threadsfor_fun at ./threadingconstructs.jl:81
#4#threadsfor_fun at ./threadingconstructs.jl:48
unknown function (ip: 0x2b762c6953bc)
jl_apply at /scratch/pkrastev/lmod_build/julia-1.5.1/src/julia.h:1690 [inlined]
start_task at /scratch/pkrastev/lmod_build/julia-1.5.1/src/task.c:707
unknown function (ip: (nil))
Allocations: 1456768 (Pool: 1456085; Big: 683); GC: 2
Aborted (core dumped)
This error is avoided by removing the Threads.@threads. Additionally, it oddly works when I reduce the number of loops (e.g., changing it to Threads.@threads for j = 1:10). I have 4 cores available.
I see this crashing as well. Here is a bit simpler example that crashes every time I run it (Julia 1.7.1, GEOS_jll 3.10.0+0, LibGEOS.jl 0.6.8)
using LibGEOS
function findintersecting(g1, g2)
n1, n2 = length(g1), length(g2)
polygon_intersection = fill(0.0,n1,n2)
Threads.@threads for i = 1:n1
for j = 1:n2
if intersects(g1[i], g2[j])
polygon_intersection[i,j] = 1
end
end
end
sum(polygon_intersection)
end
pnts = 1.0*[[-1,-1],[+1,-1],[+1,+1],[-1,+1]] .+ [0.00*randn(2) for i=1:4]
geom = [[pnts; [pnts[1]]]]
g1 = [LibGEOS.Polygon(geom) for i=1:500]
findintersecting(g1,g1)
If I add some randomness
pnts = 1.0*[[-1,-1],[+1,-1],[+1,+1],[-1,+1]] .+ [0.01*randn(2) for i=1:4]
the crashing still happens, but much more seldomly.
Even simpler crasher
import LibGEOS
function cra(n)
p = [[[-1.,-1],[+1,-1],[+1,+1],[-1,+1],[-1,-1]]]
g1 = [LibGEOS.Polygon(p) for i=1:n]
r = fill(false, n, n)
Threads.@threads for i=1:n
for j=1:n
r[i,j] = LibGEOS.intersects(g1[i],g1[j])
end
end
return r
end
cra(1000)
@jaakkor2, about https://github.com/JuliaGeo/LibGEOS.jl/issues/91#issuecomment-1003350071
I haven't tried threading yet, but I thought the idea of the _r functions that take a GEOSContext is that you create a separate context per thread to avoid these issues?
@jaakkor2, about #91 (comment)
I haven't tried threading yet, but I thought the idea of the
_rfunctions that take aGEOSContextis that you create a separate context per thread to avoid these issues?
Sorry for confusion. It seems that using one global context and threading is asking for trouble.
https://lists.osgeo.org/pipermail/geos-devel/2017-November/008157.html discussion here says that the API is thread-safe, but the implementation is not. And likely Julia interface needs modification..
Ah good to know. Would still be interesting to try it out with thread local contexts, to see if any remaining thread safety issues cause problems for examples like this.
This still crashes
import LibGEOS
ctx = [LibGEOS.GEOSContext() for i=1:Threads.nthreads()]
function cra(n, ctx)
p = [[[-1.,-1],[+1,-1],[+1,+1],[-1,+1],[-1,-1]]]
g1 = [LibGEOS.Polygon(p) for i=1:n]
r = fill(false, n, n)
Threads.@threads for i=1:n
for j=1:n
r[i,j] = LibGEOS.intersects(g1[i].ptr, g1[j].ptr, ctx[Threads.threadid()])
end
end
return r
end
cra(500, ctx)
If Threads.@threads is in front of the inner for loop, then this seems to work 🤷
import LibGEOS
ctx = [LibGEOS.GEOSContext() for i=1:Threads.nthreads()]
function cra(n, ctx)
p = [[[-1.,-1],[+1,-1],[+1,+1],[-1,+1],[-1,-1]]]
g1 = [LibGEOS.Polygon(p) for i=1:n]
r = fill(false, n, n)
for j=1:n
Threads.@threads for i=1:n
r[i,j] = LibGEOS.intersects(g1[i].ptr, g1[j].ptr, ctx[Threads.threadid()])
end
end
return r
end
cra(500, ctx)
Even simpler crasher
import LibGEOS function cra(n) p = [[[-1.,-1],[+1,-1],[+1,+1],[-1,+1],[-1,-1]]] g1 = [LibGEOS.Polygon(p) for i=1:n] r = fill(false, n, n) Threads.@threads for i=1:n for j=1:n r[i,j] = LibGEOS.intersects(g1[i],g1[j]) end end return r end cra(1000)
@jaakkor2 The following seems to work using @spawn instead of @threads.
import LibGEOS
function cra(n)
p = [[[-1.,-1],[+1,-1],[+1,+1],[-1,+1],[-1,-1]]]
g1 = [LibGEOS.Polygon(p) for i=1:n]
r = fill(false, n, n)
task = Threads.@spawn for i=1:n
for j=1:n
r[i,j] = LibGEOS.intersects(g1[i],g1[j])
end
end
wait(task)
return r
end
cra(1000)
Sorry, the above is not correct and operates on a single core. My apologies, I am quite new to Julia!
From https://trac.osgeo.org/geos/wiki/RFC3
Function names in the new API will be updated with an _r, as is the familiar C standard for reentrant/thread safe versions.
GEOSIntersection_ris used here https://github.com/JuliaGeo/LibGEOS.jl/blob/v0.6.8/src/geos_functions.jl#L423-L429EDIT: Be careful with using threads.
Unfortunately the contexts don't seem to be published.
Hmm good point about locking access to share resources. And you are right that these evaled functions don't include the optional context argument, it would be good to add them. If you want to test the contexts for threading purposes, you can replace intersection(g1, g2) with intersection(g1.ptr, g2.ptr, context) however, as was done in the example above.
Unfortunately the contexts don't seem to be published.
I've responded in https://discourse.julialang.org/t/improving-performance-on-nested-for-loops-sparsearrays-libgeos/74038/11: I don’t think it was an intentional design for it to remain internal, just that we didn’t implement it all the way through back then.
I was playing around with this and it still crashes even if all of the writing to shared resources is removed. Here is an example that crashes:
function cra(n)
ctx = [LibGEOS.GEOSContext() for i=1:Threads.nthreads()]
p = [[[-1.,-1],[+1,-1],[+1,+1],[-1,+1],[-1,-1]]]
g1 = [LibGEOS.Polygon(p) for i=1:n]
Threads.@threads for i=1:n
for j=1:n
LibGEOS.intersects(g1[i].ptr, g1[j].ptr, ctx[Threads.threadid()])
end
end
end
cra(1000)
It is the same as the example above, except that it does not write the results to a shared array and creates the context for the threads within the loop.
When debugging with the VSC debugger I get the following error:
I have also gotten the error when I switched the j and the i:
Do we think that this is caused by the C backend? I tried setting a breakpoint in malloc_error_break in VSC as suggested but it isn't breaking...
@skygering I haven't been able to reproduce any of the crashes on LibGEOS v1.7.2, see below:
| | |_| | | | (_| | | Version 1.8.0 (2022-08-17)
_/ |\__'_|_|_|\__'_| | Official https://julialang.org/ release
|__/ |
...
(@v1.8) pkg> st
Status `~/.julia/environments/v1.8/Project.toml`
...
[cf35fbd7] GeoInterface v1.0.1
[a90b1aa1] LibGEOS v0.7.2
Can I doublecheck which version of Julia and LibGEOS are you reproducing the crashes on?
For https://github.com/JuliaGeo/LibGEOS.jl/issues/91#issuecomment-1004870593, I get:
| | |_| | | | (_| | | Version 1.8.0 (2022-08-17)
_/ |\__'_|_|_|\__'_| | Official https://julialang.org/ release
|__/ |
julia> using LibGEOS
julia> function findintersecting(g1, g2)
n1, n2 = length(g1), length(g2)
polygon_intersection = fill(0.0,n1,n2)
Threads.@threads for i = 1:n1
for j = 1:n2
if intersects(g1[i], g2[j])
polygon_intersection[i,j] = 1
end
end
end
sum(polygon_intersection)
end
findintersecting (generic function with 1 method)
julia> pnts = 1.0*[[-1,-1],[+1,-1],[+1,+1],[-1,+1]] .+ [0.00*randn(2) for i=1:4]
4-element Vector{Vector{Float64}}:
[-1.0, -1.0]
[1.0, -1.0]
[1.0, 1.0]
[-1.0, 1.0]
julia> geom = [[pnts; [pnts[1]]]]
1-element Vector{Vector{Vector{Float64}}}:
[[-1.0, -1.0], [1.0, -1.0], [1.0, 1.0], [-1.0, 1.0], [-1.0, -1.0]]
julia> g1 = [LibGEOS.Polygon(geom) for i=1:500]
500-element Vector{Polygon}:
Polygon(Ptr{Nothing} @0x00006000028030c0)
Polygon(Ptr{Nothing} @0x0000600002800230)
Polygon(Ptr{Nothing} @0x00006000028004b0)
Polygon(Ptr{Nothing} @0x0000600002802c60)
Polygon(Ptr{Nothing} @0x0000600002800a50)
Polygon(Ptr{Nothing} @0x00006000028000a0)
Polygon(Ptr{Nothing} @0x0000600002803d90)
Polygon(Ptr{Nothing} @0x0000600002800410)
Polygon(Ptr{Nothing} @0x0000600002801720)
Polygon(Ptr{Nothing} @0x0000600002801c70)
⋮
Polygon(Ptr{Nothing} @0x000060000281eb20)
Polygon(Ptr{Nothing} @0x000060000281d5e0)
Polygon(Ptr{Nothing} @0x000060000281eda0)
Polygon(Ptr{Nothing} @0x000060000281ee90)
Polygon(Ptr{Nothing} @0x000060000281e120)
Polygon(Ptr{Nothing} @0x000060000281caf0)
Polygon(Ptr{Nothing} @0x000060000281f250)
Polygon(Ptr{Nothing} @0x000060000281c870)
Polygon(Ptr{Nothing} @0x000060000281f610)
Polygon(Ptr{Nothing} @0x000060000281fe30)
julia> findintersecting(g1,g1)
250000.0
For https://github.com/JuliaGeo/LibGEOS.jl/issues/91#issuecomment-1004910824, I get:
| | |_| | | | (_| | | Version 1.8.0 (2022-08-17)
_/ |\__'_|_|_|\__'_| | Official https://julialang.org/ release
|__/ |
julia> import LibGEOS
julia> ctx = [LibGEOS.GEOSContext() for i=1:Threads.nthreads()]
1-element Vector{LibGEOS.GEOSContext}:
LibGEOS.GEOSContext(Ptr{Nothing} @0x000000012f809e00)
julia> function cra(n, ctx)
p = [[[-1.,-1],[+1,-1],[+1,+1],[-1,+1],[-1,-1]]]
g1 = [LibGEOS.Polygon(p) for i=1:n]
r = fill(false, n, n)
Threads.@threads for i=1:n
for j=1:n
r[i,j] = LibGEOS.intersects(g1[i].ptr, g1[j].ptr, ctx[Threads.threadid()])
end
end
return r
end
cra (generic function with 1 method)
julia> cra(500, ctx)
500×500 Matrix{Bool}:
1 1 1 1 1 1 1 1 1 1 1 1 1 … 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 … 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
⋮ ⋮ ⋮ ⋱ ⋮ ⋮
1 1 1 1 1 1 1 1 1 1 1 1 1 … 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 … 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
For https://github.com/JuliaGeo/LibGEOS.jl/issues/91#issuecomment-1247348423, I get:
| | |_| | | | (_| | | Version 1.8.0 (2022-08-17)
_/ |\__'_|_|_|\__'_| | Official https://julialang.org/ release
|__/ |
julia> import LibGEOS
julia> function cra(n)
ctx = [LibGEOS.GEOSContext() for i=1:Threads.nthreads()]
p = [[[-1.,-1],[+1,-1],[+1,+1],[-1,+1],[-1,-1]]]
g1 = [LibGEOS.Polygon(p) for i=1:n]
Threads.@threads for i=1:n
for j=1:n
LibGEOS.intersects(g1[i].ptr, g1[j].ptr, ctx[Threads.threadid()])
end
end
end
cra (generic function with 1 method)
julia> cra(1000)
julia>
Crashes for me with LibGEOS v0.7.2 and GeoInterface v1.0.1 if i start Julia with --threads 4, does not crash if started without threads.
I also have LibGEOS v0.7.2 and GeoInterface v1.0.1. I am also running with 4 threads like @jaakkor2 mentioned. With the last example that I provided, I am crashing every time.
It doesn't seem like v0.7.2 includes the new Polygon updates, which had memory issues, but I don't think it is that given that I did not have to change the code for making Polygons from vectors.
Looking at the source code of LibGEOS, we need to construct geometries unique to a context.
So this works:
function cra(n)
@info "Using $(Threads.nthreads()) threads"
out = zeros(n, n) # don't use falses, it's not thread safe ?!
ctx = [LibGEOS.GEOSContext() for i=1:Threads.nthreads()]
p = [[[-1.,-1],[+1,-1],[+1,+1],[-1,+1],[-1,-1]]]
g1 = [LibGEOS.Polygon(p) for i=1:n]
# Generate a unique per thread/context geometry (!)
g2 = [LibGEOS.cloneGeom.(getfield.(g1, :ptr), Ref(c)) for c in ctx]
Threads.@threads for i=1:n
for j=1:n
ti = Threads.threadid()
out[i, j] = LibGEOS.intersects(g2[ti][i], g2[ti][j], ctx[ti])
end
end
sum(out)
end
cra(1000)
# [ Info: Using 8 threads
# 1.0e6
This also works for me! Thank you so much @evetion!
However, it is kind of annoying that we have to make a geometry for each thread though... I guess we shouldn't have a huge number of threads but that doesn't seem the best for performance, especially if we are working with a huge number of polygons. I don't really see a way to create less though right off the bat. This almost seems to get rid of the benefit of threads over @distributed in Julia, since we now have to have a copy of every geometry for each thread. Rather than the threads safely sharing memory, we are just creating a new version of every object in memory to prevent data races.
In fact, this seems to work as well, which does not use contexts but rather just makes nthreads copies of every polygon. Thoughts?
function cra2(n)
nthreads= Threads.nthreads()
@info "Using $(nthreads) threads"
out = zeros(n, n) # don't use falses, it's not thread safe ?!
p = [[[-1.,-1],[+1,-1],[+1,+1],[-1,+1],[-1,-1]]]
g1 = [LibGEOS.Polygon(p) for i=1:n]
g2 = [LibGEOS.cloneGeom.(getfield.(g1, :ptr)) for i in 1:nthreads]
Threads.@threads for i=1:n
for j=1:n
ti = Threads.threadid()
out[i, j] = LibGEOS.intersects(g2[ti][i], g2[ti][j])
end
end
sum(out)
end
I am also wondering if we should expose the context parameter within constructors like Polygon, since in most cases they calls functions createPolygon (which does accept a context), in addition to functions that take in the geometries like intersects, intersection, etc.
I am also wondering if we should expose the context parameter within constructors like Polygon, since in most cases they calls functions createPolygon (which does accept a context), in addition to functions that take in the geometries like intersects, intersection, etc.
Yeah the current state wasn't by design, so I think that will be very welcome :)
I have started the PR to expose context within the functions and types :)
However, I am still not convinced on the use of contexts as used above. It does stop the errors, but as I showed above this doesn't actually have to do with contexts, and rather with having a separate copy of each geometry used within each thread. This seems really inefficent. This might be the design of the underlying library, but then it seems like it really isn't thread-safe. They have just eliminated all possible race conditions by creating "separate" memory for each thread. Thoughts?
Yeah, unfortunately it seems that way.. https://www.mail-archive.com/[email protected]/msg05011.html
To keep this conversation going, I am also having issues with distributed computing using Julia's @distributed and getting seg faults.
For example:
using Distributed
addprocs(3)
@everywhere using LibGEOS
coord = [[[0.,0.],[+1.,0.],[+1.,+1.],[0.,+1.],[0.,0.]]]
coord_lst = [coord for i=1:5]
@distributed (+) for c in coord_lst
rect = LibGEOS.Polygon(c)
LibGEOS.area(rect)
end
gives 5 as expected.
However,
coord = [[[0.,0.],[+1.,0.],[+1.,+1.],[0.,+1.],[0.,0.]]]
coord_lst = [coord for i=1:5]
poly_lst = [LibGEOS.Polygon(c) for c in coord_lst]
@distributed (+) for rect in poly_lst
LibGEOS.area(rect)
end
gives a seg fault. Here is the beginning of the error message. It is really long.
The stack trace at the bottom of the error message says that there is an unsafe read:
I am a bit surprised by this error given that area seems like it should be a read only function and each polygon is in its own process.
Luckily I think that making the polygons within each process actually works for what I need to do, but I thought I would bring it up.
In the example, coord_lst is sharing the memory of coord. With this it seems to work coord_lst = [copy(coord) for i=1:5].