From 8edb481749b155e2a69f1585a826d16177e2f334 Mon Sep 17 00:00:00 2001 From: Quba1 <22771850+Quba1@users.noreply.github.com> Date: Wed, 24 Jan 2024 19:33:34 +0100 Subject: [PATCH 1/2] add parallel functions with tests --- .github/workflows/test.yml | 2 +- Cargo.toml | 14 ++++- src/lib.rs | 51 ++++++++++++++++ tests/multithreading.rs | 118 +++++++++++++++++++++++++++++++++++++ 4 files changed, 181 insertions(+), 4 deletions(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 9a156c8..469bec3 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -14,4 +14,4 @@ jobs: - uses: actions/checkout@v3 - name: Build and test with cargo run: | - cargo test \ No newline at end of file + cargo test --all-features \ No newline at end of file diff --git a/Cargo.toml b/Cargo.toml index 94aad04..faed1f2 100644 --- a/Cargo.toml +++ b/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"] + diff --git a/src/lib.rs b/src/lib.rs index 18ad10f..3d1239b 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -116,6 +116,57 @@ 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); + + #[cfg(feature = "multithreading")] + fn project_parallel(&self, coords: &[(f64, f64)]) -> Result, ProjectionError> { + use rayon::iter::{IntoParallelRefIterator, ParallelIterator}; + + let xy_points = coords + .par_iter() + .map(|(lon, lat)| self.project(*lon, *lat)) + .collect::, _>>()?; + + Ok(xy_points) + } + + #[cfg(feature = "multithreading")] + fn inverse_project_parallel( + &self, + xy_points: &[(f64, f64)], + ) -> Result, ProjectionError> { + use rayon::iter::{IntoParallelRefIterator, ParallelIterator}; + + let coords = xy_points + .par_iter() + .map(|(x, y)| self.inverse_project(*x, *y)) + .collect::, _>>()?; + + Ok(coords) + } + + #[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::>(); + + xy_points + } + + #[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::>(); + + coords + } } dyn_clone::clone_trait_object!(Projection); diff --git a/tests/multithreading.rs b/tests/multithreading.rs index 08068c6..94ccd80 100644 --- a/tests/multithreading.rs +++ b/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::>(); + let lons = std::iter::repeat(lons) + .take(100_000) + .flatten() + .collect::>(); + + let lats = (0..1000).map(|_| rng.sample(&lat_range)).collect::>(); + let lats = std::iter::repeat(lats) + .take(100_000) + .flatten() + .collect::>(); + + 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); + } +} From 15cafd045d7941921c36fb706d3d5d759d9f369d Mon Sep 17 00:00:00 2001 From: Quba1 <22771850+Quba1@users.noreply.github.com> Date: Wed, 24 Jan 2024 19:49:27 +0100 Subject: [PATCH 2/2] add docs for multithreading --- src/lib.rs | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/src/lib.rs b/src/lib.rs index 3d1239b..0a58455 100644 --- a/src/lib.rs +++ b/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; @@ -117,6 +142,8 @@ 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, ProjectionError> { use rayon::iter::{IntoParallelRefIterator, ParallelIterator}; @@ -129,6 +156,8 @@ pub trait Projection: Debug + DynClone + Send + Sync { 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, @@ -144,6 +173,7 @@ pub trait Projection: Debug + DynClone + Send + Sync { 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}; @@ -156,6 +186,7 @@ pub trait Projection: Debug + DynClone + Send + Sync { 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};