diff --git a/CHANGELOG.md b/CHANGELOG.md index 428bb21..9a62923 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,18 @@ # Changelog +## 0.9.1 2025-12-31 + +### Changed + +* Rust + * Remove explicit const-unrolling for shallow loops (eliminates nested const-unrolling) + * Improves compile time + * Minimal effect on performance in Rust benchmarks + * About 20% improvement for single-sample throughput from Python +* Python + * Lock python package version to rust version + * Improve handling of `__version__` when package is not present + ## 0.9.0 2025-12-27 ### Added diff --git a/Cargo.lock b/Cargo.lock index ceff62d..461a64f 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -221,7 +221,7 @@ checksum = "b248f5224d1d606005e02c97f5aa4e88eeb230488bcc03bc9ca4d7991399f2b5" [[package]] name = "interpn" -version = "0.9.0" +version = "0.9.1" dependencies = [ "criterion", "crunchy", diff --git a/Cargo.toml b/Cargo.toml index 5ec06e1..a9ff868 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "interpn" -version = "0.9.0" +version = "0.9.1" edition = "2024" authors = ["James Logan "] license = "MIT OR Apache-2.0" diff --git a/docs/1d_quality_of_fit_Rectilinear.html b/docs/1d_quality_of_fit_Rectilinear.html index cf01bf4..b5392f5 100644 --- a/docs/1d_quality_of_fit_Rectilinear.html +++ b/docs/1d_quality_of_fit_Rectilinear.html @@ -1,2 +1,2 @@
-
\ No newline at end of file +
\ No newline at end of file diff --git a/docs/1d_quality_of_fit_Rectilinear.svg b/docs/1d_quality_of_fit_Rectilinear.svg index b258c2f..25f0bf9 100644 --- a/docs/1d_quality_of_fit_Rectilinear.svg +++ b/docs/1d_quality_of_fit_Rectilinear.svg @@ -1 +1 @@ -05−1−0.500.51−1−0.500.51−202−200μ−100μ0−202−0.200.2−202−2−1.5−1−0.50DataInterpNScipyInterpN ErrorScipy ErrorComparison — InterpN vs. Scipy w/ Cubic InterpolantRectilinear Gridxxxf(x)ErrorQuadraticSineStepError, QuadraticError, SineError, Step \ No newline at end of file +05−1−0.500.51−1−0.500.51−202−200μ−100μ0−202−0.200.2−202−2−1.5−1−0.50DataInterpNScipyInterpN ErrorScipy ErrorComparison — InterpN vs. Scipy w/ Cubic InterpolantRectilinear Gridxxxf(x)ErrorQuadraticSineStepError, QuadraticError, SineError, Step \ No newline at end of file diff --git a/docs/1d_quality_of_fit_Regular.html b/docs/1d_quality_of_fit_Regular.html index b4dd8e4..096796f 100644 --- a/docs/1d_quality_of_fit_Regular.html +++ b/docs/1d_quality_of_fit_Regular.html @@ -1,2 +1,2 @@
-
\ No newline at end of file +
\ No newline at end of file diff --git a/docs/1d_quality_of_fit_Regular.svg b/docs/1d_quality_of_fit_Regular.svg index fd3b512..475657f 100644 --- a/docs/1d_quality_of_fit_Regular.svg +++ b/docs/1d_quality_of_fit_Regular.svg @@ -1 +1 @@ -05−1−0.500.51012−202−10f−5f0−202−0.200.2−202012DataInterpNScipyInterpN ErrorScipy ErrorComparison — InterpN vs. Scipy w/ Cubic InterpolantRegular Gridxxxf(x)ErrorQuadraticSineStepError, QuadraticError, SineError, Step \ No newline at end of file +05−1−0.500.51012−202−10f−5f0−202−0.200.2−202012DataInterpNScipyInterpN ErrorScipy ErrorComparison — InterpN vs. Scipy w/ Cubic InterpolantRegular Gridxxxf(x)ErrorQuadraticSineStepError, QuadraticError, SineError, Step \ No newline at end of file diff --git a/docs/2d_quality_of_fit_Rectilinear.html b/docs/2d_quality_of_fit_Rectilinear.html index b218428..b7e5fe4 100644 --- a/docs/2d_quality_of_fit_Rectilinear.html +++ b/docs/2d_quality_of_fit_Rectilinear.html @@ -1,2 +1,2 @@
-
\ No newline at end of file +
\ No newline at end of file diff --git a/docs/2d_quality_of_fit_Rectilinear.svg b/docs/2d_quality_of_fit_Rectilinear.svg index 39d73b8..0b7ec17 100644 --- a/docs/2d_quality_of_fit_Rectilinear.svg +++ b/docs/2d_quality_of_fit_Rectilinear.svg @@ -1 +1 @@ -Sampled data1020304050−0.00200.002Quadratic Test Function w/ Cubic InterpolantRectilinear GridTruthInterpNScipyError, InterpNError, Scipy \ No newline at end of file +Sampled data1020304050−0.00200.002Quadratic Test Function w/ Cubic InterpolantRectilinear GridTruthInterpNScipyError, InterpNError, Scipy \ No newline at end of file diff --git a/docs/2d_quality_of_fit_Regular.html b/docs/2d_quality_of_fit_Regular.html index 1e4bb85..891a69a 100644 --- a/docs/2d_quality_of_fit_Regular.html +++ b/docs/2d_quality_of_fit_Regular.html @@ -1,2 +1,2 @@
-
\ No newline at end of file +
\ No newline at end of file diff --git a/docs/2d_quality_of_fit_Regular.svg b/docs/2d_quality_of_fit_Regular.svg index dc91922..b92b02b 100644 --- a/docs/2d_quality_of_fit_Regular.svg +++ b/docs/2d_quality_of_fit_Regular.svg @@ -1 +1 @@ -Sampled data1020304050−500f0500fQuadratic Test Function w/ Cubic InterpolantRegular GridTruthInterpNScipyError, InterpNError, Scipy \ No newline at end of file +Sampled data1020304050−500f0500fQuadratic Test Function w/ Cubic InterpolantRegular GridTruthInterpNScipyError, InterpNError, Scipy \ No newline at end of file diff --git a/docs/3d_throughput_vs_nobs_cubic.html b/docs/3d_throughput_vs_nobs_cubic.html index cbb7076..b2bb5b2 100644 --- a/docs/3d_throughput_vs_nobs_cubic.html +++ b/docs/3d_throughput_vs_nobs_cubic.html @@ -1,2 +1,2 @@
-
\ No newline at end of file +
\ No newline at end of file diff --git a/docs/3d_throughput_vs_nobs_cubic.svg b/docs/3d_throughput_vs_nobs_cubic.svg index b169d33..25ba095 100644 --- a/docs/3d_throughput_vs_nobs_cubic.svg +++ b/docs/3d_throughput_vs_nobs_cubic.svg @@ -1 +1 @@ -12510251002510002510k2468ScipyInterpNInterpolation on 20x20x20 GridWithout Preallocated Output, CubicNumber of Observation PointsNormalized Throughput \ No newline at end of file +12510251002510002510k2468ScipyInterpNInterpolation on 20x20x20 GridWithout Preallocated Output, CubicNumber of Observation PointsNormalized Throughput \ No newline at end of file diff --git a/docs/3d_throughput_vs_nobs_linear.html b/docs/3d_throughput_vs_nobs_linear.html index 0f6ef9d..388f324 100644 --- a/docs/3d_throughput_vs_nobs_linear.html +++ b/docs/3d_throughput_vs_nobs_linear.html @@ -1,2 +1,2 @@
-
\ No newline at end of file +
\ No newline at end of file diff --git a/docs/3d_throughput_vs_nobs_linear.svg b/docs/3d_throughput_vs_nobs_linear.svg index f8d0941..bc4a521 100644 --- a/docs/3d_throughput_vs_nobs_linear.svg +++ b/docs/3d_throughput_vs_nobs_linear.svg @@ -1 +1 @@ -12510251002510002510k0510152025ScipyInterpNInterpolation on 20x20x20 GridWithout Preallocated Output, LinearNumber of Observation PointsNormalized Throughput \ No newline at end of file +12510251002510002510k051015202530ScipyInterpNInterpolation on 20x20x20 GridWithout Preallocated Output, LinearNumber of Observation PointsNormalized Throughput \ No newline at end of file diff --git a/docs/3d_throughput_vs_nobs_prealloc_cubic.html b/docs/3d_throughput_vs_nobs_prealloc_cubic.html index e4e23c4..29e0750 100644 --- a/docs/3d_throughput_vs_nobs_prealloc_cubic.html +++ b/docs/3d_throughput_vs_nobs_prealloc_cubic.html @@ -1,2 +1,2 @@
-
\ No newline at end of file +
\ No newline at end of file diff --git a/docs/3d_throughput_vs_nobs_prealloc_cubic.svg b/docs/3d_throughput_vs_nobs_prealloc_cubic.svg index 54dc2db..1c7fb58 100644 --- a/docs/3d_throughput_vs_nobs_prealloc_cubic.svg +++ b/docs/3d_throughput_vs_nobs_prealloc_cubic.svg @@ -1 +1 @@ -12510251002510002510k51015ScipyInterpNInterpolation on 20x20x20 GridWith Preallocated Output, CubicNumber of Observation PointsNormalized Throughput \ No newline at end of file +12510251002510002510k05101520ScipyInterpNInterpolation on 20x20x20 GridWith Preallocated Output, CubicNumber of Observation PointsNormalized Throughput \ No newline at end of file diff --git a/docs/3d_throughput_vs_nobs_prealloc_linear.html b/docs/3d_throughput_vs_nobs_prealloc_linear.html index b9787b9..005952c 100644 --- a/docs/3d_throughput_vs_nobs_prealloc_linear.html +++ b/docs/3d_throughput_vs_nobs_prealloc_linear.html @@ -1,2 +1,2 @@
-
\ No newline at end of file +
\ No newline at end of file diff --git a/docs/3d_throughput_vs_nobs_prealloc_linear.svg b/docs/3d_throughput_vs_nobs_prealloc_linear.svg index 96250b6..7f2e854 100644 --- a/docs/3d_throughput_vs_nobs_prealloc_linear.svg +++ b/docs/3d_throughput_vs_nobs_prealloc_linear.svg @@ -1 +1 @@ -12510251002510002510k0102030405060ScipyInterpNInterpolation on 20x20x20 GridWith Preallocated Output, LinearNumber of Observation PointsNormalized Throughput \ No newline at end of file +12510251002510002510k0204060ScipyInterpNInterpolation on 20x20x20 GridWith Preallocated Output, LinearNumber of Observation PointsNormalized Throughput \ No newline at end of file diff --git a/docs/4d_throughput_vs_nobs_cubic.html b/docs/4d_throughput_vs_nobs_cubic.html index 35a035b..5820a75 100644 --- a/docs/4d_throughput_vs_nobs_cubic.html +++ b/docs/4d_throughput_vs_nobs_cubic.html @@ -1,2 +1,2 @@
-
\ No newline at end of file +
\ No newline at end of file diff --git a/docs/4d_throughput_vs_nobs_cubic.svg b/docs/4d_throughput_vs_nobs_cubic.svg index af21a7b..6c0f69c 100644 --- a/docs/4d_throughput_vs_nobs_cubic.svg +++ b/docs/4d_throughput_vs_nobs_cubic.svg @@ -1 +1 @@ -12510251002510002510k2468ScipyInterpNInterpolation on 20x...x20 4D GridWithout Preallocated Output, CubicNumber of Observation PointsNormalized Throughput \ No newline at end of file +12510251002510002510k2468ScipyInterpNInterpolation on 20x...x20 4D GridWithout Preallocated Output, CubicNumber of Observation PointsNormalized Throughput \ No newline at end of file diff --git a/docs/4d_throughput_vs_nobs_linear.html b/docs/4d_throughput_vs_nobs_linear.html index 4d74de7..1215987 100644 --- a/docs/4d_throughput_vs_nobs_linear.html +++ b/docs/4d_throughput_vs_nobs_linear.html @@ -1,2 +1,2 @@
-
\ No newline at end of file +
\ No newline at end of file diff --git a/docs/4d_throughput_vs_nobs_linear.svg b/docs/4d_throughput_vs_nobs_linear.svg index e3b2ea8..d81b614 100644 --- a/docs/4d_throughput_vs_nobs_linear.svg +++ b/docs/4d_throughput_vs_nobs_linear.svg @@ -1 +1 @@ -12510251002510002510k010203040ScipyInterpNInterpolation on 20x...x20 4D GridWithout Preallocated Output, LinearNumber of Observation PointsNormalized Throughput \ No newline at end of file +12510251002510002510k010203040ScipyInterpNInterpolation on 20x...x20 4D GridWithout Preallocated Output, LinearNumber of Observation PointsNormalized Throughput \ No newline at end of file diff --git a/docs/4d_throughput_vs_nobs_prealloc_cubic.html b/docs/4d_throughput_vs_nobs_prealloc_cubic.html index 4bc15b9..3a7af65 100644 --- a/docs/4d_throughput_vs_nobs_prealloc_cubic.html +++ b/docs/4d_throughput_vs_nobs_prealloc_cubic.html @@ -1,2 +1,2 @@
-
\ No newline at end of file +
\ No newline at end of file diff --git a/docs/4d_throughput_vs_nobs_prealloc_cubic.svg b/docs/4d_throughput_vs_nobs_prealloc_cubic.svg index abc8a41..2568752 100644 --- a/docs/4d_throughput_vs_nobs_prealloc_cubic.svg +++ b/docs/4d_throughput_vs_nobs_prealloc_cubic.svg @@ -1 +1 @@ -12510251002510002510k510ScipyInterpNInterpolation on 20x...x20 4D GridWith Preallocated Output, CubicNumber of Observation PointsNormalized Throughput \ No newline at end of file +12510251002510002510k51015ScipyInterpNInterpolation on 20x...x20 4D GridWith Preallocated Output, CubicNumber of Observation PointsNormalized Throughput \ No newline at end of file diff --git a/docs/4d_throughput_vs_nobs_prealloc_linear.html b/docs/4d_throughput_vs_nobs_prealloc_linear.html index fc5767d..cfec02b 100644 --- a/docs/4d_throughput_vs_nobs_prealloc_linear.html +++ b/docs/4d_throughput_vs_nobs_prealloc_linear.html @@ -1,2 +1,2 @@
-
\ No newline at end of file +
\ No newline at end of file diff --git a/docs/4d_throughput_vs_nobs_prealloc_linear.svg b/docs/4d_throughput_vs_nobs_prealloc_linear.svg index 2593256..d171e91 100644 --- a/docs/4d_throughput_vs_nobs_prealloc_linear.svg +++ b/docs/4d_throughput_vs_nobs_prealloc_linear.svg @@ -1 +1 @@ -12510251002510002510k020406080ScipyInterpNInterpolation on 20x...x20 4D GridWith Preallocated Output, LinearNumber of Observation PointsNormalized Throughput \ No newline at end of file +12510251002510002510k020406080ScipyInterpNInterpolation on 20x...x20 4D GridWith Preallocated Output, LinearNumber of Observation PointsNormalized Throughput \ No newline at end of file diff --git a/docs/nearest_quality_of_fit.html b/docs/nearest_quality_of_fit.html index 4a8ce70..96c98d1 100644 --- a/docs/nearest_quality_of_fit.html +++ b/docs/nearest_quality_of_fit.html @@ -1,2 +1,2 @@
-
\ No newline at end of file +
\ No newline at end of file diff --git a/docs/nearest_quality_of_fit.svg b/docs/nearest_quality_of_fit.svg index 503b8a5..eec8f87 100644 --- a/docs/nearest_quality_of_fit.svg +++ b/docs/nearest_quality_of_fit.svg @@ -1 +1 @@ -Grid samples−1012−101Nearest-Neighbor Quality of Fit — InterpN vs. Scipy griddata (nearest)TruthInterpNScipyError: InterpNError: ScipyScipy - InterpN \ No newline at end of file +Grid samples−1012−101Nearest-Neighbor Quality of Fit — InterpN vs. Scipy griddata (nearest)TruthInterpNScipyError: InterpNError: ScipyScipy - InterpN \ No newline at end of file diff --git a/docs/speedup_vs_dims_1000_obs_cubic.html b/docs/speedup_vs_dims_1000_obs_cubic.html index 407e9e8..f1762c8 100644 --- a/docs/speedup_vs_dims_1000_obs_cubic.html +++ b/docs/speedup_vs_dims_1000_obs_cubic.html @@ -1,2 +1,2 @@
-
\ No newline at end of file +
\ No newline at end of file diff --git a/docs/speedup_vs_dims_1000_obs_cubic.svg b/docs/speedup_vs_dims_1000_obs_cubic.svg index fa88b00..9d5a269 100644 --- a/docs/speedup_vs_dims_1000_obs_cubic.svg +++ b/docs/speedup_vs_dims_1000_obs_cubic.svg @@ -1 +1 @@ -12345624681012InterpN Speedup vs. ScipyCubic, 1000 Observation PointsNumber of DimensionsSpeedup vs. Scipy \ No newline at end of file +12345624681012InterpN Speedup vs. ScipyCubic, 1000 Observation PointsNumber of DimensionsSpeedup vs. Scipy \ No newline at end of file diff --git a/docs/speedup_vs_dims_1000_obs_linear.html b/docs/speedup_vs_dims_1000_obs_linear.html index 3b54ab0..9a83b47 100644 --- a/docs/speedup_vs_dims_1000_obs_linear.html +++ b/docs/speedup_vs_dims_1000_obs_linear.html @@ -1,2 +1,2 @@
-
\ No newline at end of file +
\ No newline at end of file diff --git a/docs/speedup_vs_dims_1000_obs_linear.svg b/docs/speedup_vs_dims_1000_obs_linear.svg index 6be7e2f..b806e7e 100644 --- a/docs/speedup_vs_dims_1000_obs_linear.svg +++ b/docs/speedup_vs_dims_1000_obs_linear.svg @@ -1 +1 @@ -123456246810InterpN Speedup vs. ScipyLinear, 1000 Observation PointsNumber of DimensionsSpeedup vs. Scipy \ No newline at end of file +1234562468101214InterpN Speedup vs. ScipyLinear, 1000 Observation PointsNumber of DimensionsSpeedup vs. Scipy \ No newline at end of file diff --git a/docs/speedup_vs_dims_1_obs_cubic.html b/docs/speedup_vs_dims_1_obs_cubic.html index 5970a21..6976bf7 100644 --- a/docs/speedup_vs_dims_1_obs_cubic.html +++ b/docs/speedup_vs_dims_1_obs_cubic.html @@ -1,2 +1,2 @@
-
\ No newline at end of file +
\ No newline at end of file diff --git a/docs/speedup_vs_dims_1_obs_cubic.svg b/docs/speedup_vs_dims_1_obs_cubic.svg index fa80536..e6af1cd 100644 --- a/docs/speedup_vs_dims_1_obs_cubic.svg +++ b/docs/speedup_vs_dims_1_obs_cubic.svg @@ -1 +1 @@ -12345651015InterpN Speedup vs. ScipyCubic, 1 Observation PointNumber of DimensionsSpeedup vs. Scipy \ No newline at end of file +12345651015InterpN Speedup vs. ScipyCubic, 1 Observation PointNumber of DimensionsSpeedup vs. Scipy \ No newline at end of file diff --git a/docs/speedup_vs_dims_1_obs_linear.html b/docs/speedup_vs_dims_1_obs_linear.html index e647826..f72ae9d 100644 --- a/docs/speedup_vs_dims_1_obs_linear.html +++ b/docs/speedup_vs_dims_1_obs_linear.html @@ -1,2 +1,2 @@
-
\ No newline at end of file +
\ No newline at end of file diff --git a/docs/speedup_vs_dims_1_obs_linear.svg b/docs/speedup_vs_dims_1_obs_linear.svg index 19791ce..f413929 100644 --- a/docs/speedup_vs_dims_1_obs_linear.svg +++ b/docs/speedup_vs_dims_1_obs_linear.svg @@ -1 +1 @@ -123456050100150200InterpN Speedup vs. ScipyLinear, 1 Observation PointNumber of DimensionsSpeedup vs. Scipy \ No newline at end of file +123456050100150200250InterpN Speedup vs. ScipyLinear, 1 Observation PointNumber of DimensionsSpeedup vs. Scipy \ No newline at end of file diff --git a/docs/throughput_vs_dims_1000_obs.html b/docs/throughput_vs_dims_1000_obs.html index af9d9ae..b0c0854 100644 --- a/docs/throughput_vs_dims_1000_obs.html +++ b/docs/throughput_vs_dims_1000_obs.html @@ -1,2 +1,2 @@
-
\ No newline at end of file +
\ No newline at end of file diff --git a/docs/throughput_vs_dims_1000_obs.svg b/docs/throughput_vs_dims_1000_obs.svg index fe778ad..e7c43f1 100644 --- a/docs/throughput_vs_dims_1000_obs.svg +++ b/docs/throughput_vs_dims_1000_obs.svg @@ -1 +1 @@ -246250.01250.12512460.0010.010.11LinearScipy RegularGridInterpolator LinearInterpN MultilinearRegularInterpN MultilinearRectilinearInterpN NearestRegularInterpN NearestRectilinearNumpy InterpCubicScipy RegularGridInterpolator CubicInterpN MulticubicRegularInterpN MulticubicRectilinearScipy RectBivariateSpline CubicInterpolation on 4x...x4 N-Dimensional Grid1000 Observation PointsNumber of DimensionsNormalized ThroughputLinearCubic \ No newline at end of file +246250.01250.12512460.0010.010.11LinearScipy RegularGridInterpolator LinearInterpN MultilinearRegularInterpN MultilinearRectilinearInterpN NearestRegularInterpN NearestRectilinearNumpy InterpCubicScipy RegularGridInterpolator CubicInterpN MulticubicRegularInterpN MulticubicRectilinearScipy RectBivariateSpline CubicInterpolation on 4x...x4 N-Dimensional Grid1000 Observation PointsNumber of DimensionsNormalized ThroughputLinearCubic \ No newline at end of file diff --git a/docs/throughput_vs_dims_1_obs.html b/docs/throughput_vs_dims_1_obs.html index 6283bd3..4541cd9 100644 --- a/docs/throughput_vs_dims_1_obs.html +++ b/docs/throughput_vs_dims_1_obs.html @@ -1,2 +1,2 @@
-
\ No newline at end of file +
\ No newline at end of file diff --git a/docs/throughput_vs_dims_1_obs.svg b/docs/throughput_vs_dims_1_obs.svg index c8364d7..f970e92 100644 --- a/docs/throughput_vs_dims_1_obs.svg +++ b/docs/throughput_vs_dims_1_obs.svg @@ -1 +1 @@ -246250.01250.12512460.01250.1251LinearScipy RegularGridInterpolator LinearInterpN MultilinearRegularInterpN MultilinearRectilinearInterpN NearestRegularInterpN NearestRectilinearNumpy InterpCubicScipy RegularGridInterpolator CubicInterpN MulticubicRegularInterpN MulticubicRectilinearScipy RectBivariateSpline CubicInterpolation on 4x...x4 N-Dimensional Grid1 Observation PointNumber of DimensionsNormalized ThroughputLinearCubic \ No newline at end of file +2460.001250.01250.125122460.01250.1251LinearScipy RegularGridInterpolator LinearInterpN MultilinearRegularInterpN MultilinearRectilinearInterpN NearestRegularInterpN NearestRectilinearNumpy InterpCubicScipy RegularGridInterpolator CubicInterpN MulticubicRegularInterpN MulticubicRectilinearScipy RectBivariateSpline CubicInterpolation on 4x...x4 N-Dimensional Grid1 Observation PointNumber of DimensionsNormalized ThroughputLinearCubic \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 0803f94..3b90d16 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "maturin" [project] name = "interpn" -version = "0.9.0" +dynamic = ["version"] description = "N-dimensional interpolation/extrapolation methods" authors = [{ name = "James Logan", email = "jlogan03@gmail.com" }] readme = "README.md" @@ -74,12 +74,7 @@ module-name = "interpn.interpn" python-packages = ["interpn"] python-source = "src" profile = "release" - -#rustc-args = [ -# "-Cprofile-use=./scripts/pgo-profiles/interpn.profdata", -# "-Cllvm-args=-pgo-warn-missing-function" -#] - +version = { source = "cargo" } include = [ { path = "all", format = ["sdist", "wheel"] }, ] diff --git a/src/interpn/__init__.py b/src/interpn/__init__.py index cefa6c4..fd2db38 100644 --- a/src/interpn/__init__.py +++ b/src/interpn/__init__.py @@ -7,7 +7,7 @@ from typing import Literal from collections.abc import Sequence -from importlib.metadata import version +from importlib.metadata import PackageNotFoundError, version from importlib.util import find_spec import numpy as np @@ -26,7 +26,10 @@ from .nearest_regular import NearestRegular from .nearest_rectilinear import NearestRectilinear -__version__ = version("interpn") +try: + __version__ = version("interpn") +except PackageNotFoundError: + __version__ = "0.0.0" __all__ = [ "__version__", diff --git a/src/lib.rs b/src/lib.rs index 5ca2735..88d8cab 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -88,8 +88,6 @@ // expanded code that is entirely in const. #![allow(clippy::absurd_extreme_comparisons)] -use crunchy::unroll; - pub mod multilinear; pub use multilinear::{MultilinearRectilinear, MultilinearRegular}; @@ -134,10 +132,11 @@ pub(crate) fn index_arr_fixed_dims( ) -> T { let mut i = 0; - unroll! { - for j < 7 in 0..N { - i += loc[j] * dimprod[j]; - } + // unroll! { + // for j < 7 in 0..N { + for j in 0..N { + i += loc[j] * dimprod[j]; + // } } data[i] diff --git a/src/multicubic/rectilinear.rs b/src/multicubic/rectilinear.rs index a92bbdb..3ab6d2d 100644 --- a/src/multicubic/rectilinear.rs +++ b/src/multicubic/rectilinear.rs @@ -293,21 +293,19 @@ impl<'a, T: Float, const N: usize> MulticubicRectilinear<'a, T, N> { let mut store = [[T::zero(); FP]; N]; let mut acc = 1; - unroll! { - for i < 5 in 0..N { - // Populate cumulative product of higher dimensions for indexing. - // - // Each entry is the cumulative product of the size of dimensions - // higher than this one, which is the stride between blocks - // relating to a given index along each dimension. - if const { i > 0 } { - acc *= self.dims[N - i]; - } - dimprod[N - i - 1] = acc; - - // Populate lower corner and saturation flag for each dimension - (origin[i], sat[i]) = self.get_loc(x[i], i)?; + for i in 0..N { + // Populate cumulative product of higher dimensions for indexing. + // + // Each entry is the cumulative product of the size of dimensions + // higher than this one, which is the stride between blocks + // relating to a given index along each dimension. + if i > 0 { + acc *= self.dims[N - i]; } + dimprod[N - i - 1] = acc; + + // Populate lower corner and saturation flag for each dimension + (origin[i], sat[i]) = self.get_loc(x[i], i)?; } // Recursive interpolation of one dependency tree at a time @@ -317,45 +315,40 @@ impl<'a, T: Float, const N: usize> MulticubicRectilinear<'a, T, N> { macro_rules! unroll_vertices_body { ($i:ident) => { // Index, interpolate, or pass on each level of the tree - unroll! { - for j < 5 in 0..N { // const loop - - // Most of these iterations will get optimized out - if const{j == 0} { // const branch - // At leaves, index values - - unroll!{ - for k < 5 in 0..N { // const loop - // Bit pattern in an integer matches C-ordered array indexing - // so we can just use the vertex index to index into the array - // by selecting the appropriate bit from the index. - const OFFSET: usize = const{($i & (3 << (2*k))) >> (2*k)}; - loc[k] = origin[k] + OFFSET; - } - } - const STORE_IND: usize = $i % FP; - store[0][STORE_IND] = index_arr_fixed_dims(loc, dimprod, self.vals); + for j in 0..N { + // Most of these iterations will get optimized out + if j == 0 { + // const branch + // At leaves, index values + for k in 0..N { + // Bit pattern in an integer matches C-ordered array indexing + // so we can just use the vertex index to index into the array + // by selecting the appropriate bit from the index. + let offset: usize = ($i & (3 << (2 * k))) >> (2 * k); + loc[k] = origin[k] + offset; } - else { // const branch - // For other nodes, interpolate on child values - - const Q: usize = const{FP.pow(j as u32)}; - const LEVEL: bool = const {($i + 1).is_multiple_of(Q)}; - const P: usize = const{(($i + 1) / Q).saturating_sub(1) % FP}; - const IND: usize = const{j.saturating_sub(1)}; - - if LEVEL { // const branch - let grid_cell = &self.grids[IND][origin[IND]..origin[IND] + 4]; - let interped = interp_inner::( - store[IND], - grid_cell.try_into().unwrap(), - x[IND], - sat[IND], - self.linearize_extrapolation - ); - - store[j][P] = interped; - } + const STORE_IND: usize = $i % FP; + store[0][STORE_IND] = index_arr_fixed_dims(loc, dimprod, self.vals); + } else { + // const branch + // For other nodes, interpolate on child values + let q: usize = FP.pow(j as u32); + let level: bool = ($i + 1).is_multiple_of(q); + let p: usize = (($i + 1) / q).saturating_sub(1) % FP; + let ind: usize = j.saturating_sub(1); + + if level { + // const branch + let grid_cell = &self.grids[ind][origin[ind]..origin[ind] + 4]; + let interped = interp_inner::( + store[ind], + grid_cell.try_into().unwrap(), + x[ind], + sat[ind], + self.linearize_extrapolation, + ); + + store[j][p] = interped; } } } diff --git a/src/multicubic/regular.rs b/src/multicubic/regular.rs index 463474f..78d0c10 100644 --- a/src/multicubic/regular.rs +++ b/src/multicubic/regular.rs @@ -357,30 +357,28 @@ impl<'a, T: Float, const N: usize> MulticubicRegular<'a, T, N> { let mut store = [[T::zero(); FP]; N]; let mut acc = 1; - unroll! { - for i < 5 in 0..N { - // Populate cumulative product of higher dimensions for indexing. - // - // Each entry is the cumulative product of the size of dimensions - // higher than this one, which is the stride between blocks - // relating to a given index along each dimension. - if const { i > 0 } { - acc *= self.dims[N - i]; - } - dimprod[N - i - 1] = acc; - - // Populate lower corner and saturation flag for each dimension - (origin[i], sat[i]) = self.get_loc(x[i], i)?; - - // Calculate normalized delta locations - // For the cubic method, the normalized coordinate `t` is always relative - // to cube index 1 (out of 0-3) - let index_one_loc = self.starts[i] - + self.steps[i] - * ::from(origin[i] + 1) - .ok_or("Unrepresentable coordinate value")?; - dts[i] = (x[i] - index_one_loc) / self.steps[i]; + for i in 0..N { + // Populate cumulative product of higher dimensions for indexing. + // + // Each entry is the cumulative product of the size of dimensions + // higher than this one, which is the stride between blocks + // relating to a given index along each dimension. + if i > 0 { + acc *= self.dims[N - i]; } + dimprod[N - i - 1] = acc; + + // Populate lower corner and saturation flag for each dimension + (origin[i], sat[i]) = self.get_loc(x[i], i)?; + + // Calculate normalized delta locations + // For the cubic method, the normalized coordinate `t` is always relative + // to cube index 1 (out of 0-3) + let index_one_loc = self.starts[i] + + self.steps[i] + * ::from(origin[i] + 1) + .ok_or("Unrepresentable coordinate value")?; + dts[i] = (x[i] - index_one_loc) / self.steps[i]; } // Recursive interpolation of one dependency tree at a time @@ -390,43 +388,38 @@ impl<'a, T: Float, const N: usize> MulticubicRegular<'a, T, N> { macro_rules! unroll_vertices_body { ($i:ident) => { // Index, interpolate, or pass on each level of the tree - unroll! { - for j < 5 in 0..N { // const loop - - // Most of these iterations will get optimized out - if const{j == 0} { // const branch - // At leaves, index values - - unroll!{ - for k < 5 in 0..N { // const loop - // Bit pattern in an integer matches C-ordered array indexing - // so we can just use the vertex index to index into the array - // by selecting the appropriate bit from the index. - const OFFSET: usize = const{($i & (3 << (2*k))) >> (2*k)}; - loc[k] = origin[k] + OFFSET; - } - } - const STORE_IND: usize = $i % FP; - store[0][STORE_IND] = index_arr_fixed_dims(loc, dimprod, self.vals); + for j in 0..N { + // Most of these iterations will get optimized out + if j == 0 { + // const branch + // At leaves, index values + for k in 0..N { + // Bit pattern in an integer matches C-ordered array indexing + // so we can just use the vertex index to index into the array + // by selecting the appropriate bit from the index. + let offset: usize = ($i & (3 << (2 * k))) >> (2 * k); + loc[k] = origin[k] + offset; } - else { // const branch - // For other nodes, interpolate on child values - - const Q: usize = const{FP.pow(j as u32)}; - const LEVEL: bool = const {($i + 1).is_multiple_of(Q)}; - const P: usize = const{(($i + 1) / Q).saturating_sub(1) % FP}; - const IND: usize = const{j.saturating_sub(1)}; - - if LEVEL { // const branch - let interped = interp_inner::( - &store[IND], - dts[IND], - sat[IND], - self.linearize_extrapolation - ); - - store[j][P] = interped; - } + const STORE_IND: usize = $i % FP; + store[0][STORE_IND] = index_arr_fixed_dims(loc, dimprod, self.vals); + } else { + // const branch + // For other nodes, interpolate on child values + let q: usize = FP.pow(j as u32); + let level: bool = ($i + 1).is_multiple_of(q); + let p: usize = (($i + 1) / q).saturating_sub(1) % FP; + let ind: usize = j.saturating_sub(1); + + if level { + // const branch + let interped = interp_inner::( + &store[ind], + dts[ind], + sat[ind], + self.linearize_extrapolation, + ); + + store[j][p] = interped; } } } diff --git a/src/multilinear/rectilinear.rs b/src/multilinear/rectilinear.rs index 00c2a84..a9181b1 100644 --- a/src/multilinear/rectilinear.rs +++ b/src/multilinear/rectilinear.rs @@ -255,21 +255,19 @@ impl<'a, T: Float, const N: usize> MultilinearRectilinear<'a, T, N> { let mut store = [[T::zero(); FP]; N]; let mut acc = 1; - unroll! { - for i < 7 in 0..N { - // Populate cumulative product of higher dimensions for indexing. - // - // Each entry is the cumulative product of the size of dimensions - // higher than this one, which is the stride between blocks - // relating to a given index along each dimension. - if const { i > 0 } { - acc *= self.dims[N - i]; - } - dimprod[N - i - 1] = acc; - - // Populate lower corner and saturation flag for each dimension - origin[i] = self.get_loc(x[i], i)?; + for i in 0..N { + // Populate cumulative product of higher dimensions for indexing. + // + // Each entry is the cumulative product of the size of dimensions + // higher than this one, which is the stride between blocks + // relating to a given index along each dimension. + if i > 0 { + acc *= self.dims[N - i]; } + dimprod[N - i - 1] = acc; + + // Populate lower corner and saturation flag for each dimension + origin[i] = self.get_loc(x[i], i)?; } // Recursive interpolation of one dependency tree at a time @@ -279,21 +277,17 @@ impl<'a, T: Float, const N: usize> MultilinearRectilinear<'a, T, N> { unroll! { for i < 64 in 0..nverts { // const loop // Index, interpolate, or pass on each level of the tree - unroll!{ - for j < 7 in 0..N { // const loop + for j in 0..N { // Most of these iterations will get optimized out - if const{j == 0} { // const branch + if j == 0 { // const branch // At leaves, index values - - unroll!{ - for k < 7 in 0..N { // const loop + for k in 0..N { // Bit pattern in an integer matches C-ordered array indexing // so we can just use the vertex index to index into the array // by selecting the appropriate bit from the index. - const OFFSET: usize = const{(i & (1 << k)) >> k}; - loc[k] = origin[k] + OFFSET; - } + let offset: usize = (i & (1 << k)) >> k; + loc[k] = origin[k] + offset; } const STORE_IND: usize = i % FP; store[0][STORE_IND] = index_arr_fixed_dims(loc, dimprod, self.vals); @@ -301,29 +295,29 @@ impl<'a, T: Float, const N: usize> MultilinearRectilinear<'a, T, N> { else { // const branch // For other nodes, interpolate on child values - const Q: usize = const{FP.pow(j as u32)}; - const LEVEL: bool = const {(i + 1).is_multiple_of(Q)}; - const P: usize = const{((i + 1) / Q).saturating_sub(1) % FP}; - const IND: usize = const{j.saturating_sub(1)}; + let q: usize = FP.pow(j as u32); + let level: bool = (i + 1).is_multiple_of(q); + let p: usize = ((i + 1) / q).saturating_sub(1) % FP; + let ind: usize = j.saturating_sub(1); - if LEVEL { // const branch - let x0 = self.grids[IND][origin[IND]]; - let x1 = self.grids[IND][origin[IND] + 1]; + if level { // const branch + let x0 = self.grids[ind][origin[ind]]; + let x1 = self.grids[ind][origin[ind] + 1]; let step = x1 - x0; - let t = (x[IND] - x0) / step; + let t = (x[ind] - x0) / step; - let y0 = store[IND][0]; - let dy = store[IND][1] - y0; + let y0 = store[ind][0]; + let dy = store[ind][1] - y0; #[cfg(not(feature = "fma"))] let interped = y0 + t * dy; #[cfg(feature = "fma")] let interped = t.mul_add(dy, y0); - store[j][P] = interped; + store[j][p] = interped; } } - } + } } } diff --git a/src/multilinear/regular.rs b/src/multilinear/regular.rs index c83b1af..12aeaf5 100644 --- a/src/multilinear/regular.rs +++ b/src/multilinear/regular.rs @@ -313,31 +313,29 @@ impl<'a, T: Float, const N: usize> MultilinearRegular<'a, T, N> { let mut store = [[T::zero(); FP]; N]; let mut acc = 1; - unroll! { - for i < 7 in 0..N { - // Populate cumulative product of higher dimensions for indexing. - // - // Each entry is the cumulative product of the size of dimensions - // higher than this one, which is the stride between blocks - // relating to a given index along each dimension. - if const { i > 0 } { - acc *= self.dims[N - i]; - } - dimprod[N - i - 1] = acc; - - // Populate lower corner and saturation flag for each dimension - origin[i] = self.get_loc(x[i], i)?; - let origin_f = ::from(origin[i]).ok_or("Unrepresentable coordinate value")?; + for i in 0..N { + // Populate cumulative product of higher dimensions for indexing. + // + // Each entry is the cumulative product of the size of dimensions + // higher than this one, which is the stride between blocks + // relating to a given index along each dimension. + if i > 0 { + acc *= self.dims[N - i]; + } + dimprod[N - i - 1] = acc; + // Populate lower corner and saturation flag for each dimension + origin[i] = self.get_loc(x[i], i)?; + let origin_f = + ::from(origin[i]).ok_or("Unrepresentable coordinate value")?; - // Calculate normalized delta locations - #[cfg(not(feature = "fma"))] - let index_zero_loc = self.starts[i] + self.steps[i] * origin_f; - #[cfg(feature = "fma")] - let index_zero_loc = self.steps[i].mul_add(origin_f, self.starts[i]); + // Calculate normalized delta locations + #[cfg(not(feature = "fma"))] + let index_zero_loc = self.starts[i] + self.steps[i] * origin_f; + #[cfg(feature = "fma")] + let index_zero_loc = self.steps[i].mul_add(origin_f, self.starts[i]); - dts[i] = (x[i] - index_zero_loc) / self.steps[i]; - } + dts[i] = (x[i] - index_zero_loc) / self.steps[i]; } // Recursive interpolation of one dependency tree at a time @@ -347,21 +345,17 @@ impl<'a, T: Float, const N: usize> MultilinearRegular<'a, T, N> { unroll! { for i < 64 in 0..nverts { // const loop // Index, interpolate, or pass on each level of the tree - unroll!{ - for j < 7 in 0..N { // const loop + for j in 0..N { // Most of these iterations will get optimized out - if const{j == 0} { // const branch + if j == 0 { // const branch // At leaves, index values - - unroll!{ - for k < 7 in 0..N { // const loop + for k in 0..N { // Bit pattern in an integer matches C-ordered array indexing // so we can just use the vertex index to index into the array // by selecting the appropriate bit from the index. - const OFFSET: usize = const{(i & (1 << k)) >> k}; - loc[k] = origin[k] + OFFSET; - } + let offset: usize = (i & (1 << k)) >> k; + loc[k] = origin[k] + offset; } const STORE_IND: usize = i % FP; store[0][STORE_IND] = index_arr_fixed_dims(loc, dimprod, self.vals); @@ -369,25 +363,24 @@ impl<'a, T: Float, const N: usize> MultilinearRegular<'a, T, N> { else { // const branch // For other nodes, interpolate on child values - const Q: usize = const{FP.pow(j as u32)}; - const LEVEL: bool = const {(i + 1).is_multiple_of(Q)}; - const P: usize = const{((i + 1) / Q).saturating_sub(1) % FP}; - const IND: usize = const{j.saturating_sub(1)}; + let q: usize = FP.pow(j as u32); + let level: bool = (i + 1).is_multiple_of(q); + let p: usize = ((i + 1) / q).saturating_sub(1) % FP; + let ind: usize = j.saturating_sub(1); - if LEVEL { // const branch - let y0 = store[IND][0]; - let dy = store[IND][1] - y0; - let t = dts[IND]; + if level { // const branch + let y0 = store[ind][0]; + let dy = store[ind][1] - y0; + let t = dts[ind]; #[cfg(not(feature = "fma"))] let interped = y0 + t * dy; #[cfg(feature = "fma")] let interped = t.mul_add(dy, y0); - store[j][P] = interped; + store[j][p] = interped; } } - } } } } diff --git a/src/nearest/rectilinear.rs b/src/nearest/rectilinear.rs index c3e2fc1..d7c9c44 100644 --- a/src/nearest/rectilinear.rs +++ b/src/nearest/rectilinear.rs @@ -25,7 +25,6 @@ //! rectilinear::interpn_alloc(grids, &z, &obs).unwrap(); //! ``` use crate::index_arr_fixed_dims; -use crunchy::unroll; use num_traits::Float; /// Evaluate multicubic interpolation on a regular grid in up to 6 dimensions. @@ -204,36 +203,31 @@ impl<'a, T: Float, const N: usize> NearestRectilinear<'a, T, N> { let half = T::one() / two; let mut acc = 1; - unroll! { - for i < 7 in 0..N { - // Populate cumulative product of higher dimensions for indexing. - // - // Each entry is the cumulative product of the size of dimensions - // higher than this one, which is the stride between blocks - // relating to a given index along each dimension. - if const { i > 0 } { - acc *= self.dims[N - i]; - } - dimprod[N - i - 1] = acc; - - // Populate lower corner and saturation flag for each dimension - let origin = self.get_loc(x[i], i)?; - - // Determine normalized delta - let x0 = self.grids[i][origin]; - let x1 = self.grids[i][origin + 1]; - let step = x1 - x0; - let dt = (x[i] - x0) / step; - - // Determine nearest index for this dimension based on distance - let offset = if dt <= half { - 0 - } else { - 1 - }; - - loc[i] = origin + offset; + for i in 0..N { + // NOTE: it is not necessary to explicitly unroll this + // Populate cumulative product of higher dimensions for indexing. + // + // Each entry is the cumulative product of the size of dimensions + // higher than this one, which is the stride between blocks + // relating to a given index along each dimension. + if i > 0 { + acc *= self.dims[N - i]; } + dimprod[N - i - 1] = acc; + + // Populate lower corner and saturation flag for each dimension + let origin = self.get_loc(x[i], i)?; + + // Determine normalized delta + let x0 = self.grids[i][origin]; + let x1 = self.grids[i][origin + 1]; + let step = x1 - x0; + let dt = (x[i] - x0) / step; + + // Determine nearest index for this dimension based on distance + let offset = if dt <= half { 0 } else { 1 }; + + loc[i] = origin + offset; } let interped = index_arr_fixed_dims(loc, dimprod, self.vals); diff --git a/src/nearest/regular.rs b/src/nearest/regular.rs index a8fdc21..a572301 100644 --- a/src/nearest/regular.rs +++ b/src/nearest/regular.rs @@ -27,7 +27,6 @@ //! regular::interpn_alloc(&dims, &starts, &steps, &z, &obs).unwrap(); //! ``` use crate::index_arr_fixed_dims; -use crunchy::unroll; use num_traits::{Float, NumCast}; /// Evaluate nearest-neighbor interpolation on a regular grid in up to 6 dimensions. @@ -252,42 +251,38 @@ impl<'a, T: Float, const N: usize> NearestRegular<'a, T, N> { let half = T::one() / two; let mut acc = 1; - unroll! { - for i < 7 in 0..N { - // Populate cumulative product of higher dimensions for indexing. - // - // Each entry is the cumulative product of the size of dimensions - // higher than this one, which is the stride between blocks - // relating to a given index along each dimension. - if const { i > 0 } { - acc *= self.dims[N - i]; - } - dimprod[N - i - 1] = acc; - - // Populate lower corner and saturation flag for each dimension - let origin = self.get_loc(x[i], i)?; - let origin_f = ::from(origin).ok_or("Unrepresentable coordinate value")?; - - // Calculate normalized delta locations - #[cfg(not(feature = "fma"))] - let index_zero_loc = self.starts[i] + self.steps[i] * origin_f; - #[cfg(feature = "fma")] - let index_zero_loc = self.steps[i].mul_add(origin_f, self.starts[i]); - - let dt = (x[i] - index_zero_loc) / self.steps[i]; - - // Determine nearest index for this dimension based on distance. - // NOTE: This method, despite including a division operation, - // is about 10-20% faster than just checking the left and right - // distances against each other directly. - let offset = if dt <= half { - 0 - } else { - 1 - }; - - loc[i] = origin + offset; + for i in 0..N { + // NOTE: it is not necessary to explicitly unroll this + // Populate cumulative product of higher dimensions for indexing. + // + // Each entry is the cumulative product of the size of dimensions + // higher than this one, which is the stride between blocks + // relating to a given index along each dimension. + if i > 0 { + acc *= self.dims[N - i]; } + dimprod[N - i - 1] = acc; + + // Populate lower corner and saturation flag for each dimension + let origin = self.get_loc(x[i], i)?; + let origin_f = + ::from(origin).ok_or("Unrepresentable coordinate value")?; + + // Calculate normalized delta locations + #[cfg(not(feature = "fma"))] + let index_zero_loc = self.starts[i] + self.steps[i] * origin_f; + #[cfg(feature = "fma")] + let index_zero_loc = self.steps[i].mul_add(origin_f, self.starts[i]); + + let dt = (x[i] - index_zero_loc) / self.steps[i]; + + // Determine nearest index for this dimension based on distance. + // NOTE: This method, despite including a division operation, + // is about 10-20% faster than just checking the left and right + // distances against each other directly. + let offset = if dt <= half { 0 } else { 1 }; + + loc[i] = origin + offset; } let interped = index_arr_fixed_dims(loc, dimprod, self.vals); diff --git a/uv.lock b/uv.lock index 62940f6..b0fe85f 100644 --- a/uv.lock +++ b/uv.lock @@ -292,7 +292,6 @@ wheels = [ [[package]] name = "interpn" -version = "0.9.0" source = { editable = "." } dependencies = [ { name = "numpy", version = "2.2.6", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" },