Merge pull request #4 from ScaleWeather/parallelism

Parallelism

.github/workflows/test.yml

@@ -14,4 +14,4 @@ jobs:
     - uses: actions/checkout@v3
     - name: Build and test with cargo
       run: |
-        cargo test
+        cargo test --all-features

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "mappers"
-version = "0.5.2"
+version = "0.6.0"
 edition = "2021"
 description = "Pure Rust geographical projections library"
 repository = "https://github.com/ScaleWeather/mappers"
@@ -18,8 +18,16 @@ float-cmp = { version = "0.9", default-features = false, features = ["std"] }
 thiserror = { version = "1.0", default-features = false }
 geographiclib-rs = { version = "^0.2.3", default-features = false }
 dyn-clone = { version = "^1.0.10", default-features = false }
-const_soft_float = { version = "^0.1.2", default-features = false }   
+const_soft_float = { version = "^0.1.2", default-features = false }
+rayon = { version = "^1.8", default-features = false, optional = true}
 
 [dev-dependencies]
 float-cmp = { version = "0.9", default-features = false, features = ["std"] }
-proj = { version = "0.27", default-features = false }
+# REQUIRES: libproj-dev, clang, libtiff-dev and sqlite (binary!)
+proj = { version = "0.27", default-features = false, features = ["pkg_config"] }
+rand = { version = "0.8"}
+
+[features]
+default = ["multithreading"]
+multithreading = ["rayon"]
+

src/lib.rs

@@ -55,10 +55,35 @@
 //!# Ok(())
 //!# }
 //!```
+//! 
 //! Some projections are mathematically exactly inversible, and technically
 //! geographical coordinates projected and inverse projected should be identical.
 //! However, in practice limitations of floating-point arithmetics will
 //! introduce some errors along the way, as shown in the example above.
+//! 
+//! ## Multithreading
+//! 
+//! For projecting multiple coordinates at once, the crate provides `_parallel`
+//! functions that are available in a (default) `multithreading` feature. These functions
+//! use `rayon` crate to parallelize the projection process. They are provided 
+//! mainly for convenience, as they are not much different than calling 
+//! `.par_iter()` on a slice of coordinates and mapping the projection function over it.
+//! 
+//!```
+//!# use mappers::{Ellipsoid, projections::LambertConformalConic, Projection, ProjectionError};
+//!#
+//!# fn main() -> Result<(), ProjectionError> {
+//! let lcc = LambertConformalConic::new(2.0, 0.0, 30.0, 60.0, Ellipsoid::WGS84)?;
+//!
+//! // Parallel functions use slices of tuples as input and output
+//! let geographic_coordinates = vec![(6.8651, 45.8326); 10];
+//!
+//! let map_coordinates = lcc.project_parallel(&geographic_coordinates)?;
+//! let inversed_coordinates = lcc.inverse_project_parallel(&map_coordinates)?;
+//!
+//!# Ok(())
+//!# }
+//!```
 
 use dyn_clone::DynClone;
 use std::fmt::Debug;
@@ -116,6 +141,63 @@ pub trait Projection: Debug + DynClone + Send + Sync {
 
     /// Same as [`Projection::inverse_project()`] but does not check the result.
     fn inverse_project_unchecked(&self, x: f64, y: f64) -> (f64, f64);
+
+    /// Function analogous to [`Projection::project()`] for projecting
+    /// multiple coordinates at once using multithreading.
+    #[cfg(feature = "multithreading")]
+    fn project_parallel(&self, coords: &[(f64, f64)]) -> Result<Vec<(f64, f64)>, ProjectionError> {
+        use rayon::iter::{IntoParallelRefIterator, ParallelIterator};
+
+        let xy_points = coords
+            .par_iter()
+            .map(|(lon, lat)| self.project(*lon, *lat))
+            .collect::<Result<Vec<_>, _>>()?;
+
+        Ok(xy_points)
+    }
+
+    /// Function analogous to [`Projection::inverse_project()`] for projecting
+    /// multiple coordinates at once using multithreading.
+    #[cfg(feature = "multithreading")]
+    fn inverse_project_parallel(
+        &self,
+        xy_points: &[(f64, f64)],
+    ) -> Result<Vec<(f64, f64)>, ProjectionError> {
+        use rayon::iter::{IntoParallelRefIterator, ParallelIterator};
+
+        let coords = xy_points
+            .par_iter()
+            .map(|(x, y)| self.inverse_project(*x, *y))
+            .collect::<Result<Vec<_>, _>>()?;
+
+        Ok(coords)
+    }
+
+    /// Same as [`Projection::project_parallel()`] but does not check the result.
+    #[cfg(feature = "multithreading")]
+    fn project_parallel_unchecked(&self, coords: &[(f64, f64)]) -> Vec<(f64, f64)> {
+        use rayon::iter::{IntoParallelRefIterator, ParallelIterator};
+
+        let xy_points = coords
+            .par_iter()
+            .map(|(lon, lat)| self.project_unchecked(*lon, *lat))
+            .collect::<Vec<_>>();
+
+        xy_points
+    }
+
+    /// Same as [`Projection::inverse_project_parallel()`] but does not check the result.
+    #[cfg(feature = "multithreading")]
+    fn inverse_project_parallel_unchecked(&self, xy_points: &[(f64, f64)]) -> Vec<(f64, f64)> {
+        use rayon::iter::{IntoParallelRefIterator, ParallelIterator};
+
+        let coords = xy_points
+            .par_iter()
+            .map(|(x, y)| self.inverse_project_unchecked(*x, *y))
+            .collect::<Vec<_>>();
+
+        coords
+    }
 }
 
 dyn_clone::clone_trait_object!(Projection);

tests/multithreading.rs

@@ -30,3 +30,121 @@ fn arc_interop() {
         assert_approx_eq!(f64, geo_coords.1, 32.68850065409422, epsilon = 0.000_000_1);
     }
 }
+
+#[cfg(feature = "multithreading")]
+#[test]
+fn parallelism_checked() {
+    let proj = AzimuthalEquidistant::new(30.0, 30.0, Ellipsoid::WGS84).unwrap();
+
+    let coords = vec![(25.0, 45.0); 100];
+    let xy_points = vec![(200_000.0, 300_000.0); 100];
+
+    let map_coords = proj.project_parallel(&coords).unwrap();
+    let geo_coords = proj.inverse_project_parallel(&xy_points).unwrap();
+
+    for map_coord in map_coords {
+        assert_approx_eq!(f64, map_coord.0, -398563.2994422894, epsilon = 0.000_000_1);
+        assert_approx_eq!(f64, map_coord.1, 1674853.7525355904, epsilon = 0.000_000_1);
+    }
+
+    for geo_coord in geo_coords {
+        assert_approx_eq!(f64, geo_coord.0, 32.132000374279365, epsilon = 0.000_000_1);
+        assert_approx_eq!(f64, geo_coord.1, 32.68850065409422, epsilon = 0.000_000_1);
+    }
+}
+
+#[cfg(feature = "multithreading")]
+#[test]
+fn parallelism_unchecked() {
+    let proj = AzimuthalEquidistant::new(30.0, 30.0, Ellipsoid::WGS84).unwrap();
+
+    let coords = vec![(25.0, 45.0); 100];
+    let xy_points = vec![(200_000.0, 300_000.0); 100];
+
+    let map_coords = proj.project_parallel_unchecked(&coords);
+    let geo_coords = proj.inverse_project_parallel_unchecked(&xy_points);
+
+    for map_coord in map_coords {
+        assert_approx_eq!(f64, map_coord.0, -398563.2994422894, epsilon = 0.000_000_1);
+        assert_approx_eq!(f64, map_coord.1, 1674853.7525355904, epsilon = 0.000_000_1);
+    }
+
+    for geo_coord in geo_coords {
+        assert_approx_eq!(f64, geo_coord.0, 32.132000374279365, epsilon = 0.000_000_1);
+        assert_approx_eq!(f64, geo_coord.1, 32.68850065409422, epsilon = 0.000_000_1);
+    }
+}
+
+#[cfg(feature = "multithreading")]
+#[test]
+fn parallelism_arc_interop() {
+    let proj = AzimuthalEquidistant::new(30.0, 30.0, Ellipsoid::WGS84).unwrap();
+    let proj = Arc::new(proj);
+    let mut handles = vec![];
+
+    for _ in 0..10 {
+        let proj = Arc::clone(&proj);
+        let handle = thread::spawn(move || {
+            let coords = vec![(25.0, 45.0); 100];
+            let xy_points = vec![(200_000.0, 300_000.0); 100];
+
+            let map_coords = proj.project_parallel_unchecked(&coords);
+            let geo_coords = proj.inverse_project_parallel_unchecked(&xy_points);
+
+            (map_coords, geo_coords)
+        });
+        handles.push(handle);
+    }
+
+    for handle in handles {
+        let (map_coords, geo_coords) = handle.join().unwrap();
+
+        for map_coord in map_coords {
+            assert_approx_eq!(f64, map_coord.0, -398563.2994422894, epsilon = 0.000_000_1);
+            assert_approx_eq!(f64, map_coord.1, 1674853.7525355904, epsilon = 0.000_000_1);
+        }
+
+        for geo_coord in geo_coords {
+            assert_approx_eq!(f64, geo_coord.0, 32.132000374279365, epsilon = 0.000_000_1);
+            assert_approx_eq!(f64, geo_coord.1, 32.68850065409422, epsilon = 0.000_000_1);
+        }
+    }
+}
+
+#[cfg(feature = "multithreading")]
+#[ignore = "this test takes a long time to run"]
+#[test]
+fn parallelism_checked_long() {
+    use mappers::projections::LambertConformalConic;
+    use rand::{distributions::Uniform, Rng};
+
+    let proj = LambertConformalConic::new(0.0, 30.0, 30.0, 60.0, Ellipsoid::WGS84).unwrap();
+    let mut rng = rand::thread_rng();
+    let lon_range = Uniform::new(-15.0, 15.0);
+    let lat_range = Uniform::new(40.0, 50.0);
+
+    let lons = (0..1000).map(|_| rng.sample(&lon_range)).collect::<Vec<_>>();
+    let lons = std::iter::repeat(lons)
+        .take(100_000)
+        .flatten()
+        .collect::<Vec<_>>();
+
+    let lats = (0..1000).map(|_| rng.sample(&lat_range)).collect::<Vec<_>>();
+    let lats = std::iter::repeat(lats)
+        .take(100_000)
+        .flatten()
+        .collect::<Vec<_>>();
+
+    let ref_coords: Vec<(f64, f64)> = lons.iter().copied().zip(lats.iter().copied()).collect();
+
+    println!("Starting projection");
+
+    let map_coords = proj.project_parallel(&ref_coords).unwrap();
+
+    let geo_coords = proj.inverse_project_parallel(&map_coords).unwrap();
+
+    for (ref_coord, geo_coord) in ref_coords.iter().zip(geo_coords.iter()) {
+        assert_approx_eq!(f64, ref_coord.0, geo_coord.0, epsilon = 0.000_1);
+        assert_approx_eq!(f64, ref_coord.1, geo_coord.1, epsilon = 0.000_1);
+    }
+}