geos/par: warmin up prepped geom with a contains() call makes it thread-safe

.github/workflows/rust.yml

@@ -12,9 +12,17 @@ jobs:
             override: true
             components: rustfmt, clippy
             profile: minimal
+
+      - uses: actions/setup-python@v2
+        with:
+          python-version: 3.9
+          architecture: x64
+
       - name: Install deps
         run: |
           sudo apt-get -y update
           sudo apt-get -y install build-essential libssl-dev
+          pip install numpy
+
       - run: cargo build --features geos/static --verbose
       - run: cargo test --features geos/static --verbose

Cargo.toml

@@ -20,16 +20,21 @@ path = "src/devel/make_bitmap.rs"
 [dependencies]
 geos = { version = "8" }
 lazy_static = "1.4.0"
-numpy = "0.14"
+numpy = { version = "0.14", features = [ "rayon" ] }
 pyo3 = { version = "0.14" }
 roaring = "0.7.0"
 rust-embed = "5.9.0"
 xz2 = "0.1.6"
-openssl = { version = "0.10", features = ["vendored"] }
+openssl = { version = "0.10" }
+ndarray = { version = "0.15.3", features = [ "rayon" ] }
+
+[dev-dependencies]
+rayon = "1"
 
 [patch.crates-io]
 geos-sys = { version = "2", git = "https://github.com/gauteh/geos-sys", branch = "static-build" }
 geos = { version = "8", git = "https://github.com/gauteh/geos", branch = "static" }
+#geos-sys = { path = "../../dev/geos-sys" }
 
 [build-dependencies]
 reqwest = { version = "0.11", features = [ "blocking" ] }

pyproject.toml

@@ -3,4 +3,4 @@ requires = ["maturin>=0.11,<0.12"]
 build-backend = "maturin"
 
 [tool.maturin]
-cargo-extra-args = "--features extension-module,geos/static"
+cargo-extra-args = "--features extension-module,geos/static,openssl/vendored"

src/lib.rs

@@ -127,6 +127,20 @@ impl RoaringLandmask {
         )
         .to_owned()
     }
+
+    pub fn contains_many_par(
+        &self,
+        py: Python,
+        x: PyReadonlyArrayDyn<f64>,
+        y: PyReadonlyArrayDyn<f64>,
+    ) -> Py<PyArray<bool, numpy::IxDyn>> {
+        let x = x.as_array();
+        let y = y.as_array();
+
+        use ndarray::Zip;
+        let contains = Zip::from(&x).and(&y).par_map_collect(|x, y| self.contains(*x, *y));
+        PyArray::from_owned_array(py, contains).to_owned()
+    }
 }
 
 #[cfg(test)]
@@ -157,4 +171,72 @@ mod tests {
 
         b.iter(|| mask.contains(5., 65.6))
     }
+
+    #[bench]
+    fn test_contains_many(b: &mut Bencher) {
+        let mask = RoaringLandmask::new().unwrap();
+
+        let (x, y): (Vec<f64>, Vec<f64>) = (0..360 * 2)
+            .map(|v| v as f64 * 0.5 - 180.)
+            .map(|x| {
+                (0..180 * 2)
+                    .map(|y| y as f64 * 0.5 - 90.)
+                    .map(move |y| (x, y))
+            })
+            .flatten()
+            .unzip();
+
+        pyo3::prepare_freethreaded_python();
+        pyo3::Python::with_gil(|py| {
+
+            let x = PyArray::from_vec(py, x);
+            let y = PyArray::from_vec(py, y);
+
+            println!("testing {} points..", x.len());
+
+            b.iter(|| {
+                let len = x.len();
+
+                let x = x.to_dyn().readonly();
+                let y = y.to_dyn().readonly();
+
+                let onland = mask.contains_many(py, x, y);
+                assert!(onland.as_ref(py).len() == len);
+            })
+        })
+    }
+
+    #[bench]
+    fn test_contains_many_par(b: &mut Bencher) {
+        let mask = RoaringLandmask::new().unwrap();
+
+        let (x, y): (Vec<f64>, Vec<f64>) = (0..360 * 2)
+            .map(|v| v as f64 * 0.5 - 180.)
+            .map(|x| {
+                (0..180 * 2)
+                    .map(|y| y as f64 * 0.5 - 90.)
+                    .map(move |y| (x, y))
+            })
+            .flatten()
+            .unzip();
+
+        pyo3::prepare_freethreaded_python();
+        pyo3::Python::with_gil(|py| {
+
+            let x = PyArray::from_vec(py, x);
+            let y = PyArray::from_vec(py, y);
+
+            println!("testing {} points..", x.len());
+
+            b.iter(|| {
+                let len = x.len();
+
+                let x = x.to_dyn().readonly();
+                let y = y.to_dyn().readonly();
+
+                let onland = mask.contains_many_par(py, x, y);
+                assert!(onland.as_ref(py).len() == len);
+            })
+        })
+    }
 }

src/shapes.rs

@@ -23,14 +23,24 @@ impl Drop for Gshhg {
     }
 }
 
-// PreparedGeometry is Send, Geometry is Send. *mut Geometry is never modified.
+// PreparedGeometry is Send+Sync, Geometry is Send+Sync. *mut Geometry is never modified.
 unsafe impl Send for Gshhg {}
+unsafe impl Sync for Gshhg {}
+
+// `PreparededGeometry::contains` needs a call to `contains` before it is thread-safe:
+// https://github.com/georust/geos/issues/95
+fn warmup_prepped(prepped: &PreparedGeometry<'_>) {
+    let point = CoordSeq::new_from_vec(&[&[0.0, 0.0]]).unwrap();
+    let point = Geometry::create_point(point).unwrap();
+    prepped.contains(&point).unwrap();
+}
 
 impl Clone for Gshhg {
     fn clone(&self) -> Self {
         let gptr = self.geom.clone();
         debug_assert!(gptr != self.geom);
         let prepped = unsafe { (&*gptr).to_prepared_geom().unwrap() };
+        warmup_prepped(&prepped);
 
         Gshhg {
             geom: gptr,
@@ -45,6 +55,7 @@ impl Gshhg {
         let gptr = Box::into_raw(bxd);
         let prepped = unsafe { (&*gptr).to_prepared_geom() }
             .map_err(|_| io::Error::new(io::ErrorKind::InvalidData, "cannot prepare geomtry"))?;
+        warmup_prepped(&prepped);
 
         Ok(Gshhg {
             geom: gptr,
@@ -53,15 +64,19 @@ impl Gshhg {
     }
 
     pub fn from_compressed<P: AsRef<Path>>(path: P) -> io::Result<Gshhg> {
+        let g = Gshhg::get_geometry_from_compressed(path)?;
+
+        Gshhg::from_geom(g)
+    }
+
+    pub fn get_geometry_from_compressed<P: AsRef<Path>>(path: P) -> io::Result<Geometry<'static>> {
         let fd = File::open(path)?;
         let fd = io::BufReader::new(fd);
         let mut fd = xz2::bufread::XzDecoder::new(fd);
         let mut buf = Vec::new();
         fd.read_to_end(&mut buf)?;
 
-        let g = geos::Geometry::new_from_wkb(&buf).unwrap();
-
-        Gshhg::from_geom(g)
+        Ok(geos::Geometry::new_from_wkb(&buf).unwrap())
     }
 }
 

tests/test_geos_par_prepped.rs

@@ -0,0 +1,54 @@
+use roaring_landmask::Gshhg;
+use geos::{CoordSeq, Geom, Geometry};
+
+#[ignore]
+#[test]
+fn test_par_prepped_no_warmup() {
+    use rayon::prelude::*;
+
+    // let s = Gshhg::new().unwrap();
+    // let prepped = &s.prepped;
+
+    let g = Gshhg::get_geometry_from_compressed(
+        "gshhs/gshhs_f_-180.000000E-90.000000N180.000000E90.000000N.wkb.xz",
+    ).unwrap();
+    let prepped = g.to_prepared_geom().unwrap();
+
+    (0..10000).into_par_iter().for_each(|k| {
+
+        let x = k % 180;
+        let y = (k / 180) % 90;
+
+
+        let point = CoordSeq::new_from_vec(&[&[x as f64, y as f64]]).unwrap();
+        let point = Geometry::create_point(point).unwrap();
+        prepped.contains(&point).unwrap();
+    });
+}
+
+#[test]
+fn test_par_prepped_with_warmup() {
+    use rayon::prelude::*;
+
+    // let s = Gshhg::new().unwrap();
+    // let prepped = &s.prepped;
+
+    let g = Gshhg::get_geometry_from_compressed(
+        "gshhs/gshhs_f_-180.000000E-90.000000N180.000000E90.000000N.wkb.xz",
+    ).unwrap();
+    let prepped = g.to_prepared_geom().unwrap();
+
+    let point = CoordSeq::new_from_vec(&[&[10., 50.]]).unwrap();
+    let point = Geometry::create_point(point).unwrap();
+    prepped.contains(&point).unwrap();
+
+    (0..10000).into_par_iter().for_each(|k| {
+        let x = k % 180;
+        let y = (k / 180) % 90;
+
+
+        let point = CoordSeq::new_from_vec(&[&[x as f64, y as f64]]).unwrap();
+        let point = Geometry::create_point(point).unwrap();
+        prepped.contains(&point).unwrap();
+    });
+}

tests/test_roaring.py

@@ -21,3 +21,37 @@ def test_landmask_many(benchmark):
 
   print ("points:", len(xx.ravel()))
   benchmark(l.contains_many, xx.ravel(), yy.ravel())
+
+def test_landmask_many_few(benchmark):
+  l = RoaringLandmask.new()
+
+  x = np.linspace(-180, 180, 10)
+  y = np.linspace(-90, 90, 5)
+
+  xx, yy = np.meshgrid(x,y)
+
+  print ("points:", len(xx.ravel()))
+  benchmark(l.contains_many, xx.ravel(), yy.ravel())
+
+def test_landmask_many_par(benchmark):
+  l = RoaringLandmask.new()
+
+  x = np.arange(-180, 180, .5)
+  y = np.arange(-90, 90, .5)
+
+  xx, yy = np.meshgrid(x,y)
+
+  print ("points:", len(xx.ravel()))
+  benchmark(l.contains_many_par, xx.ravel(), yy.ravel())
+
+def test_landmask_many_par_few(benchmark):
+  l = RoaringLandmask.new()
+
+  x = np.linspace(-180, 180, 10)
+  y = np.linspace(-90, 90, 5)
+
+  xx, yy = np.meshgrid(x,y)
+
+  print ("points:", len(xx.ravel()))
+  benchmark(l.contains_many_par, xx.ravel(), yy.ravel())
+