From 69e23bddf6468efc9962b322871c8025e840d5cd Mon Sep 17 00:00:00 2001 From: Daksha Date: Sun, 4 Jan 2026 02:08:57 +0530 Subject: [PATCH] Fix cubic interpolation bounds handling --- src/grid/cubic.py | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/src/grid/cubic.py b/src/grid/cubic.py index 6d1a389b..0bfd7679 100644 --- a/src/grid/cubic.py +++ b/src/grid/cubic.py @@ -20,7 +20,7 @@ r"""Hyper Rectangular Grid In Either Two or Three Dimensions.""" import numpy as np -from scipy.interpolate import CubicSpline, RegularGridInterpolator +from scipy.interpolate import CubicSpline, RegularGridInterpolator, interpn from sympy import symbols from sympy.functions.combinatorial.numbers import bell @@ -165,6 +165,20 @@ def interpolate(self, points, values, use_log=False, nu_x=0, nu_y=0, nu_z=0, met interpolate = RegularGridInterpolator((x, y, z), values, method=method) return interpolate(points) + # Use scipy interpn for cubic interpolation when no derivatives are requested. + if method == "cubic" and nu_x == 0 and nu_y == 0 and nu_z == 0 and not use_log: + x, y, z = self.get_points_along_axes() + values_reshaped = values.reshape(self.shape) + + return interpn( + (x, y, z), + values_reshaped, + points, + method="cubic", + bounds_error=False, + fill_value=None, + ) + # Interpolate the Z-Axis. def z_spline(z, x_index, y_index, nu_z=nu_z): # x_index, y_index is assumed to be in the grid while z is not assumed. @@ -293,7 +307,8 @@ def index_to_coordinates(self, index): r"""Convert grid point index to its (i, j) or (i, j, k) integer coordinates in the grid. For 3D grid it has a shape of :math:`(N_x, N_y, N_z)` denoting the number of points in - :math:`x`\, :math:`y`\, and :math:`z` directions. So, each grid point has a :math:`(i, j, k)` + :math:`x`\, :math:`y`\, and :math:`z` directions. + So, each grid point has a :math:`(i, j, k)` integer coordinate where :math:`0 <= i <= N_x - 1`\, :math:`0 <= j <= N_y - 1`\, and :math:`0 <= k <= N_z - 1`\. Two-dimensional case similarly follows. Assumes the grid enumerates in the last coordinate first (with others fixed), following the @@ -708,7 +723,7 @@ def from_cube(cls, fname, weight="Trapezoid", return_data=False): if not fname.endswith(".cube"): raise ValueError("Argument fname should be a cube file with *.cube extension!") - with open(fname, "rt") as f: + with open(fname) as f: # skip the title and second line f.readline() f.readline()