Merge pull request #623 from RocketPy-Team/mnt/fix-warnings-and-numpy…

…-upgrade

MNT: Fix warnings in test suite and adds support for numpy 2.0

CHANGELOG.md

@@ -65,6 +65,7 @@ You can install this version by running `pip install rocketpy==1.3.0`
 
 ### Changed
 
+- MNT: Fix warnings in test suite and adds support for numpy 2.0 [#623](https://github.com/RocketPy-Team/RocketPy/pull/623)
 - REL: Bump versioning to RocketPy v1.3.0 [#614](https://github.com/RocketPy-Team/RocketPy/pull/614)
 - ENH: Adds StochasticModel.visualize_attributes() method [#612](https://github.com/RocketPy-Team/RocketPy/pull/612)
 - DOC: Monte carlo documentation updates [#607](https://github.com/RocketPy-Team/RocketPy/pull/607)

rocketpy/_encoders.py

@@ -29,7 +29,7 @@ def default(self, o):
             ),
         ):
             return int(o)
-        elif isinstance(o, (np.float_, np.float16, np.float32, np.float64)):
+        elif isinstance(o, (np.float16, np.float32, np.float64)):
             return float(o)
         elif isinstance(o, np.ndarray):
             return o.tolist()

rocketpy/environment/environment.py

@@ -570,7 +570,7 @@ def set_gravity_model(self, gravity=None):
         >>> g_0 = 9.80665
         >>> env_cte_g = Environment(gravity=g_0)
         >>> env_cte_g.gravity([0, 100, 1000])
-        [9.80665, 9.80665, 9.80665]
+        [np.float64(9.80665), np.float64(9.80665), np.float64(9.80665)]
 
         It's also possible to variate the gravity acceleration by defining
         its function of height:
@@ -667,7 +667,7 @@ def set_elevation(self, elevation="Open-Elevation"):
         None
         """
         if elevation != "Open-Elevation" and elevation != "SRTM":
-            self.elevation = elevation
+            self.elevation = float(elevation)
         # elif elevation == "SRTM" and self.latitude != None and self.longitude != None:
         #     # Trigger the authentication flow.
         #     #ee.Authenticate()
@@ -682,7 +682,7 @@ def set_elevation(self, elevation="Open-Elevation"):
         #     self.elevation = elev
 
         elif self.latitude is not None and self.longitude is not None:
-            self.elevation = self.__fetch_open_elevation()
+            self.elevation = float(self.__fetch_open_elevation())
             print("Elevation received: ", self.elevation)
         else:
             raise ValueError(
@@ -1722,7 +1722,7 @@ def process_windy_atmosphere(self, model="ECMWF"):
         self.max_expected_height = max(altitude_array[0], altitude_array[-1])
 
         # Get elevation data from file
-        self.elevation = response["header"]["elevation"]
+        self.elevation = float(response["header"]["elevation"])
 
         # Compute info data
         self.atmospheric_model_init_date = netCDF4.num2date(
@@ -2423,9 +2423,9 @@ def process_forecast_reanalysis(self, file, dictionary):
                 f_x_y2 = ((x2 - x) / (x2 - x1)) * f_x1_y2 + (
                     (x - x1) / (x2 - x1)
                 ) * f_x2_y2
-                self.elevation = ((y2 - y) / (y2 - y1)) * f_x_y1 + (
-                    (y - y1) / (y2 - y1)
-                ) * f_x_y2
+                self.elevation = float(
+                    ((y2 - y) / (y2 - y1)) * f_x_y1 + ((y - y1) / (y2 - y1)) * f_x_y2
+                )
             except:
                 raise ValueError(
                     "Unable to read surface elevation data. Check file and dictionary."
@@ -2793,9 +2793,9 @@ def process_ensemble(self, file, dictionary):
                 f_x_y2 = ((x2 - x) / (x2 - x1)) * f_x1_y2 + (
                     (x - x1) / (x2 - x1)
                 ) * f_x2_y2
-                self.elevation = ((y2 - y) / (y2 - y1)) * f_x_y1 + (
-                    (y - y1) / (y2 - y1)
-                ) * f_x_y2
+                self.elevation = float(
+                    ((y2 - y) / (y2 - y1)) * f_x_y1 + ((y - y1) / (y2 - y1)) * f_x_y2
+                )
             except:
                 raise ValueError(
                     "Unable to read surface elevation data. Check file and dictionary."
@@ -3090,15 +3090,15 @@ def calculate_density_profile(self):
 
         >>> env = Environment()
         >>> env.calculate_density_profile()
-        >>> env.density(0)
+        >>> float(env.density(0))
         1.225000018124288
 
         Creating an Environment object and calculating the density
         at 1000m above Sea Level:
 
         >>> env = Environment()
         >>> env.calculate_density_profile()
-        >>> env.density(1000)
+        >>> float(env.density(1000))
         1.1116193933422585
         """
         # Retrieve pressure P, gas constant R and temperature T

rocketpy/mathutils/function.py

@@ -16,7 +16,14 @@
 import numpy as np
 from scipy import integrate, linalg, optimize
 
-NUMERICAL_TYPES = (float, int, complex, np.ndarray, np.integer, np.floating)
+# Numpy 1.x compatibility,
+# TODO: remove these lines when all dependencies support numpy>=2.0.0
+if np.lib.NumpyVersion(np.__version__) >= "2.0.0b1":
+    from numpy import trapezoid  # pragma: no cover
+else:
+    from numpy import trapz as trapezoid  # pragma: no cover
+
+NUMERICAL_TYPES = (float, int, complex, np.integer, np.floating)
 INTERPOLATION_TYPES = {
     "linear": 0,
     "polynomial": 1,
@@ -549,9 +556,9 @@ def set_discrete(
             func.set_interpolation(interpolation)
             func.set_extrapolation(extrapolation)
         elif func.__dom_dim__ == 2:
-            lower = 2 * [lower] if isinstance(lower, (int, float)) else lower
-            upper = 2 * [upper] if isinstance(upper, (int, float)) else upper
-            sam = 2 * [samples] if isinstance(samples, (int, float)) else samples
+            lower = 2 * [lower] if isinstance(lower, NUMERICAL_TYPES) else lower
+            upper = 2 * [upper] if isinstance(upper, NUMERICAL_TYPES) else upper
+            sam = 2 * [samples] if isinstance(samples, NUMERICAL_TYPES) else samples
             # Create nodes to evaluate function
             xs = np.linspace(lower[0], upper[0], sam[0])
             ys = np.linspace(lower[1], upper[1], sam[1])
@@ -829,26 +836,26 @@ def get_value(self, *args):
         ...    [(0, 0), (1, 1), (1.5, 2.25), (2, 4), (2.5, 6.25), (3, 9), (4, 16)]
         ... )
         >>> f3.get_value(2)
-        4.0
+        np.float64(4.0)
         >>> f3.get_value(2.5)
-        6.25
+        np.float64(6.25)
         >>> f3.get_value([1, 2, 3])
-        [1.0, 4.0, 9.0]
+        [np.float64(1.0), np.float64(4.0), np.float64(9.0)]
         >>> f3.get_value([1, 2.5, 4.0])
-        [1.0, 6.25, 16.0]
+        [np.float64(1.0), np.float64(6.25), np.float64(16.0)]
 
         Testing with ndarray source (2 dimensions):
         >>> f4 = Function(
         ...    [(0, 0, 0), (1, 1, 1), (1, 2, 2), (2, 4, 8), (3, 9, 27)]
         ... )
         >>> f4.get_value(1, 1)
-        1.0
+        np.float64(1.0)
         >>> f4.get_value(2, 4)
-        8.0
+        np.float64(8.0)
         >>> abs(f4.get_value(1, 1.5) - 1.5) < 1e-2  # the interpolation is not perfect
-        True
+        np.True_
         >>> f4.get_value(3, 9)
-        27.0
+        np.float64(27.0)
         """
         if len(args) != self.__dom_dim__:
             raise ValueError(
@@ -860,7 +867,7 @@ def get_value(self, *args):
             # if the function is 1-D:
             if self.__dom_dim__ == 1:
                 # if the args is a simple number (int or float)
-                if isinstance(args[0], (int, float, complex)):
+                if isinstance(args[0], NUMERICAL_TYPES):
                     return self.source(args[0])
                 # if the arguments are iterable, we map and return a list
                 if isinstance(args[0], Iterable):
@@ -869,7 +876,7 @@ def get_value(self, *args):
             # if the function is n-D:
             else:
                 # if each arg is a simple number (int or float)
-                if all(isinstance(arg, (int, float, complex)) for arg in args):
+                if all(isinstance(arg, NUMERICAL_TYPES) for arg in args):
                     return self.source(*args)
                 # if each arg is iterable, we map and return a list
                 if all(isinstance(arg, Iterable) for arg in args):
@@ -880,7 +887,7 @@ def get_value(self, *args):
 
         # Returns value for other interpolation type
         else:  # interpolation is "polynomial", "spline", "akima" or "linear"
-            if isinstance(args[0], (int, float, complex, np.integer)):
+            if isinstance(args[0], NUMERICAL_TYPES):
                 args = [list(args)]
 
         x = list(args[0])
@@ -1330,19 +1337,19 @@ def plot_2d(
         if callable(self.source):
             # Determine boundaries
             lower = [0, 0] if lower is None else lower
-            lower = 2 * [lower] if isinstance(lower, (int, float)) else lower
+            lower = 2 * [lower] if isinstance(lower, NUMERICAL_TYPES) else lower
             upper = [10, 10] if upper is None else upper
-            upper = 2 * [upper] if isinstance(upper, (int, float)) else upper
+            upper = 2 * [upper] if isinstance(upper, NUMERICAL_TYPES) else upper
         else:
             # Determine boundaries
             x_data = self.x_array
             y_data = self.y_array
             x_min, x_max = x_data.min(), x_data.max()
             y_min, y_max = y_data.min(), y_data.max()
             lower = [x_min, y_min] if lower is None else lower
-            lower = 2 * [lower] if isinstance(lower, (int, float)) else lower
+            lower = 2 * [lower] if isinstance(lower, NUMERICAL_TYPES) else lower
             upper = [x_max, y_max] if upper is None else upper
-            upper = 2 * [upper] if isinstance(upper, (int, float)) else upper
+            upper = 2 * [upper] if isinstance(upper, NUMERICAL_TYPES) else upper
             # Plot data points if force_data = True
             if force_data:
                 axes.scatter(x_data, y_data, self.source[:, -1])
@@ -1513,7 +1520,8 @@ def compare_plots(
                     ax.scatter(points[0], points[1], marker="o")
 
         # Setup legend
-        ax.legend(loc="best", shadow=True)
+        if any([plot[1] for plot in plots]):
+            ax.legend(loc="best", shadow=True)
 
         # Turn on grid and set title and axis
         plt.grid(True)
@@ -1875,7 +1883,9 @@ def __add__(self, other):
                 return Function(lambda x: (self.get_value_opt(x) + other(x)))
         # If other is Float except...
         except AttributeError:
-            if isinstance(other, NUMERICAL_TYPES):
+            if isinstance(other, NUMERICAL_TYPES) or self.__is_single_element_array(
+                other
+            ):
                 # Check if Function object source is array or callable
                 if isinstance(self.source, np.ndarray):
                     # Operate on grid values
@@ -1997,7 +2007,9 @@ def __mul__(self, other):
             source = np.column_stack((self.x_array, self.y_array * other.y_array))
             outputs = f"({self.__outputs__[0]}*{other.__outputs__[0]})"
             return Function(source, inputs, outputs, interp, extrap)
-        elif isinstance(other, NUMERICAL_TYPES):
+        elif isinstance(other, NUMERICAL_TYPES) or self.__is_single_element_array(
+            other
+        ):
             if not self_source_is_array:
                 return Function(lambda x: (self.get_value_opt(x) * other), inputs)
             source = np.column_stack((self.x_array, np.multiply(self.y_array, other)))
@@ -2080,7 +2092,9 @@ def __truediv__(self, other):
                 return Function(lambda x: (self.get_value_opt(x) / other(x)))
         # If other is Float except...
         except AttributeError:
-            if isinstance(other, NUMERICAL_TYPES):
+            if isinstance(other, NUMERICAL_TYPES) or self.__is_single_element_array(
+                other
+            ):
                 # Check if Function object source is array or callable
                 if isinstance(self.source, np.ndarray):
                     # Operate on grid values
@@ -2119,7 +2133,7 @@ def __rtruediv__(self, other):
             A Function object which gives the result of other(x)/self(x).
         """
         # Check if Function object source is array and other is float
-        if isinstance(other, NUMERICAL_TYPES):
+        if isinstance(other, NUMERICAL_TYPES) or self.__is_single_element_array(other):
             if isinstance(self.source, np.ndarray):
                 # Operate on grid values
                 ys = other / self.y_array
@@ -2187,7 +2201,9 @@ def __pow__(self, other):
                 return Function(lambda x: (self.get_value_opt(x) ** other(x)))
         # If other is Float except...
         except AttributeError:
-            if isinstance(other, NUMERICAL_TYPES):
+            if isinstance(other, NUMERICAL_TYPES) or self.__is_single_element_array(
+                other
+            ):
                 # Check if Function object source is array or callable
                 if isinstance(self.source, np.ndarray):
                     # Operate on grid values
@@ -2226,7 +2242,7 @@ def __rpow__(self, other):
             A Function object which gives the result of other(x)**self(x).
         """
         # Check if Function object source is array and other is float
-        if isinstance(other, NUMERICAL_TYPES):
+        if isinstance(other, NUMERICAL_TYPES) or self.__is_single_element_array(other):
             if isinstance(self.source, np.ndarray):
                 # Operate on grid values
                 ys = other**self.y_array
@@ -2375,7 +2391,7 @@ def integral(self, a, b, numerical=False):
                     # self.__extrapolation__ = 'zero'
                     pass
         elif self.__interpolation__ == "linear" and numerical is False:
-            # Integrate from a to b using np.trapz
+            # Integrate from a to b using np.trapezoid
             x_data = self.x_array
             y_data = self.y_array
             # Get data in interval
@@ -2394,8 +2410,7 @@ def integral(self, a, b, numerical=False):
                 y_integration_data = np.concatenate(([self(a)], y_integration_data))
                 x_integration_data = np.concatenate((x_integration_data, [b]))
                 y_integration_data = np.concatenate((y_integration_data, [self(b)]))
-            # Integrate using np.trapz
-            ans = np.trapz(y_integration_data, x_integration_data)
+            ans = trapezoid(y_integration_data, x_integration_data)
         else:
             # Integrate numerically
             ans, _ = integrate.quad(self, a, b, epsabs=1e-4, epsrel=1e-3, limit=1000)
@@ -2606,16 +2621,16 @@ def is_strictly_bijective(self):
         Examples
         --------
         >>> f = Function([[0, 0], [1, 1], [2, 4]])
-        >>> f.isbijective()
-        True
-        >>> f.is_strictly_bijective()
+        >>> f.isbijective() == True
         True
+        >>> f.is_strictly_bijective() == True
+        np.True_
 
         >>> f = Function([[-1, 1], [0, 0], [1, 1], [2, 4]])
         >>> f.isbijective()
         False
         >>> f.is_strictly_bijective()
-        False
+        np.False_
 
         A Function which is not "strictly" bijective, but is bijective, can be
         constructed as x^2 defined at -1, 0 and 2.
@@ -2624,7 +2639,7 @@ def is_strictly_bijective(self):
         >>> f.isbijective()
         True
         >>> f.is_strictly_bijective()
-        False
+        np.False_
         """
         if isinstance(self.source, np.ndarray):
             # Assuming domain is sorted, range must also be
@@ -2900,6 +2915,10 @@ def savetxt(
             file.write(header_line + newline)
             np.savetxt(file, data_points, fmt=fmt, delimiter=delimiter, newline=newline)
 
+    @staticmethod
+    def __is_single_element_array(var):
+        return isinstance(var, np.ndarray) and var.size == 1
+
     # Input validators
     def __validate_source(self, source):
         """Used to validate the source parameter for creating a Function object.
@@ -2957,7 +2976,7 @@ def __validate_source(self, source):
                 )
             return source
 
-        if isinstance(source, (int, float)):
+        if isinstance(source, NUMERICAL_TYPES):
             # Convert number source into vectorized lambda function
             temp = 1 * source
 
@@ -3172,7 +3191,7 @@ def calc_output(func, inputs):
         output_data = []
         for key in sorted(source.keys()):
             i = np.linspace(key[0], key[1], datapoints)
-            i = i[~np.in1d(i, input_data)]
+            i = i[~np.isin(i, input_data)]
             input_data = np.concatenate((input_data, i))
 
             f = Function(source[key])

rocketpy/plots/flight_plots.py

@@ -66,6 +66,14 @@ def trajectory_3d(self):
         min_x = min(self.flight.x[:, 1])
         max_y = max(self.flight.y[:, 1])
         min_y = min(self.flight.y[:, 1])
+        min_xy = min(min_x, min_y)
+        max_xy = max(max_x, max_y)
+
+        # avoids errors when x_lim and y_lim are the same
+        if abs(min_z - max_z) < 1e-5:
+            max_z += 1
+        if abs(min_xy - max_xy) < 1e-5:
+            max_xy += 1
 
         _ = plt.figure(figsize=(9, 9))
         ax1 = plt.subplot(111, projection="3d")
@@ -109,8 +117,8 @@ def trajectory_3d(self):
         ax1.set_ylabel("Y - North (m)")
         ax1.set_zlabel("Z - Altitude Above Ground Level (m)")
         ax1.set_title("Flight Trajectory")
-        ax1.set_xlim(min_x, max_x)
-        ax1.set_ylim(min_y, max_y)
+        ax1.set_xlim(min_xy, max_xy)
+        ax1.set_ylim(min_xy, max_xy)
         ax1.set_zlim(min_z, max_z)
         ax1.view_init(15, 45)
         ax1.set_box_aspect(None, zoom=0.95)  # 95% for label adjustment