Fix 3D spline evaluation layout handling for multivariate data

.gitignore

@@ -7,6 +7,7 @@
 # These are backup files generated by rustfmt
 **/*.rs.bk
 
-# IDE
+# IDE & Agents
 .idea/
 .vscode/
+.codex

CHANGELOG.md

@@ -1,8 +1,15 @@
 # Changelog
 
-## v0.4.0 (2022)
+## v0.5.0 (01.05.2026)
+
+* Fix 3D spline evaluation layout handling for multivariate data.
+* Update dependencies for `ndarray` 0.17 and replace deprecated reshape calls with explicit layout-aware reshaping.
+
+
+## v0.4.0 (20.04.2022)
 * Update to ndarray 0.15 and sprs 0.11
 
+
 ## v0.3.0 (22.03.2020)
 
 * Add `Real` pub trait. Trait `Real` is only implemented for `f32` and `f64`, 

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "csaps"
-version = "0.4.0"
+version = "0.5.0"
 authors = ["Eugene Prilepin <esp.home@gmail.com>"]
 
 description = "Cubic spline approximation (smoothing)"
@@ -17,19 +17,18 @@ edition = "2021"
 
 
 [badges]
-travis-ci = { repository = "espdev/csaps-rs", branch = "master" }
 coveralls = { repository = "espdev/csaps-rs", branch = "master", service = "github" }
 
 
 [dependencies]
 num-traits = "0.2.14"
-ndarray = "0.16.1"
+ndarray = "0.17.2"
 sprs-ldl = "0.10.0"
-sprs = "0.11.2"
+sprs = "0.11.4"
 almost = "0.2.0"
-itertools = "0.13.0"
-thiserror = "1.0.30" 
+itertools = "0.14.0"
+thiserror = "2.0.18"
 
 [dev-dependencies]
 approx = "0.5.1"
-ndarray = {version = "0.16.1", features = ["approx"]}
+ndarray = {version = "0.17.2", features = ["approx"]}

src/ndarrayext.rs

@@ -1,12 +1,25 @@
 use itertools::Itertools;
-use ndarray::{prelude::*, IntoDimension, Slice};
+use ndarray::{prelude::*, IntoDimension, Order, Slice};
 
 use crate::{
     util::dim_from_vec,
     CsapsError::{ReshapeFrom2d, ReshapeTo2d},
     Real, Result,
 };
 
+pub(crate) fn reshape_order<T, D>(data: &ArrayView<'_, T, D>) -> Order
+where
+    D: Dimension,
+{
+    if data.is_standard_layout() {
+        Order::RowMajor
+    } else if data.ndim() > 1 && data.raw_view().reversed_axes().is_standard_layout() {
+        Order::ColumnMajor
+    } else {
+        Order::RowMajor
+    }
+}
+
 pub fn diff<'a, T: 'a, D, V>(data: V, axis: Option<Axis>) -> Array<T, D>
 where
     T: Real<T>,
@@ -43,7 +56,10 @@ where
     let axis_size = shape[axis.0];
     let new_shape = [numel / axis_size, axis_size];
 
-    match data_view.permuted_axes(axes).into_shape(new_shape) {
+    let data_view = data_view.permuted_axes(axes);
+    let order = reshape_order(&data_view);
+
+    match data_view.into_shape_with_order((new_shape, order)) {
         Ok(view_2d) => Ok(view_2d),
         Err(error) => Err(ReshapeTo2d {
             input_shape: shape,
@@ -62,7 +78,9 @@ where
     let shape = data.shape().to_vec();
     let new_shape = [shape[0..(ndim - 1)].iter().product(), shape[ndim - 1]];
 
-    match data.into_shape(new_shape) {
+    let order = reshape_order(&data);
+
+    match data.into_shape_with_order((new_shape, order)) {
         Ok(data_2d) => Ok(data_2d),
         Err(error) => Err(ReshapeTo2d {
             input_shape: shape,
@@ -93,7 +111,9 @@ where
     let new_shape: D = dim_from_vec(ndim, new_shape_vec.clone());
     let data_view = data.into();
 
-    match data_view.into_shape(new_shape) {
+    let order = reshape_order(&data_view);
+
+    match data_view.into_shape_with_order((new_shape, order)) {
         Ok(view_nd) => {
             let mut axes_tmp: Vec<usize> = (0..ndim).collect();
             let end_axis = axes_tmp.pop().unwrap();

src/ndg/evaluate.rs

@@ -1,6 +1,10 @@
 use ndarray::{Array, ArrayView1, Dimension};
 
-use crate::{ndarrayext::to_2d_simple, util::dim_from_vec, NdSpline, Real};
+use crate::{
+    ndarrayext::{reshape_order, to_2d_simple},
+    util::dim_from_vec,
+    NdSpline, Real,
+};
 
 use super::{util::permute_axes, GridCubicSmoothingSpline, NdGridSpline};
 
@@ -36,8 +40,10 @@ where
                 coeffs_shape[ndim_m1] = xi_ax.len();
                 let shape: D = dim_from_vec(self.ndim, coeffs_shape);
 
+                let order = reshape_order(&coeffs_2d.view());
+
                 coeffs_2d
-                    .into_shape(shape)
+                    .into_shape_with_order((shape, order))
                     .unwrap()
                     .permuted_axes(permuted_axes.clone())
                     .to_owned()

src/ndg/make.rs

@@ -1,7 +1,10 @@
 use ndarray::Dimension;
 
 use crate::util::dim_from_vec;
-use crate::{ndarrayext::to_2d_simple, CubicSmoothingSpline, Real, RealRef, Result};
+use crate::{
+    ndarrayext::{reshape_order, to_2d_simple},
+    CubicSmoothingSpline, Real, RealRef, Result,
+};
 
 use super::{util::permute_axes, GridCubicSmoothingSpline, NdGridSpline};
 
@@ -54,9 +57,11 @@ where
                 coeffs_shape[ndim_m1] = spline.pieces() * spline.order();
                 let new_shape: D = dim_from_vec(ndim, coeffs_shape);
 
-                spline
-                    .coeffs()
-                    .into_shape(new_shape)
+                let coeffs = spline.coeffs();
+                let order = reshape_order(&coeffs);
+
+                coeffs
+                    .into_shape_with_order((new_shape, order))
                     .unwrap()
                     .permuted_axes(permuted_axes.clone())
                     .to_owned()

src/umv/evaluate.rs

@@ -44,14 +44,15 @@ where
         // Returns NxM array of coeffs values for given 1xM indices array
         // where N is ndim and M is the size of xi
         let get_indexed_coeffs = |inds: &Array1<usize>| {
-            // Returns Nx1 2-d array of coeffs by given index
-            let coeffs_by_index = |&index| coeffs.slice(s![.., index]).insert_axis(Axis(1));
+            let mut indexed_coeffs = Array2::<T>::zeros((coeffs.nrows(), inds.len()));
 
-            // Get the M-sized vector of coeffs values Nx1 arrays
-            // for all dimensions for 1xM indices array
-            let indexed_coeffs: Vec<_> = inds.iter().map(coeffs_by_index).collect();
+            for (col, &index) in inds.iter().enumerate() {
+                indexed_coeffs
+                    .column_mut(col)
+                    .assign(&coeffs.slice(s![.., index]));
+            }
 
-            concatenate(Axis(1), &indexed_coeffs).unwrap()
+            indexed_coeffs
         };
 
         // Vectorized computing the spline pieces (polynoms) on the given data sites