From b88e8d6126e7490cfe905e1be17381f8290ca340 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20BRIOL?= <fbriol@gmail.com>
Date: Sun, 8 Sep 2024 10:16:53 +0200
Subject: [PATCH] =?UTF-8?q?=F0=9F=8E=A8=20Apply=20pre-commit=20hooks.?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 CODE_OF_CONDUCT.md                            |  2 +-
 scripts/stubgen.py                            |  3 +-
 src/pyinterp/backends/xarray.py               |  2 +-
 src/pyinterp/core/__init__.pyi                |  2 +-
 src/pyinterp/core/array.pyi                   |  8 +--
 src/pyinterp/core/include/pyinterp/axis.hpp   |  4 +-
 .../core/include/pyinterp/bivariate.hpp       |  8 +--
 .../core/include/pyinterp/dateutils.hpp       |  9 ++--
 .../core/include/pyinterp/detail/axis.hpp     | 14 +++--
 .../pyinterp/detail/axis/container.hpp        | 28 +++++-----
 .../pyinterp/detail/geometry/crossover.hpp    |  4 +-
 .../pyinterp/detail/geometry/rtree.hpp        | 34 ++++++-------
 .../pyinterp/detail/interpolation/akima.hpp   |  8 +--
 .../pyinterp/detail/interpolation/bicubic.hpp |  4 +-
 .../detail/interpolation/bilinear.hpp         |  4 +-
 .../pyinterp/detail/interpolation/cspline.hpp |  8 +--
 .../detail/interpolation/cspline_base.hpp     |  6 +--
 .../detail/interpolation/cspline_periodic.hpp |  4 +-
 .../detail/interpolation/interpolator_1d.hpp  | 12 ++---
 .../detail/interpolation/interpolator_2d.hpp  |  4 +-
 .../detail/interpolation/polynomial.hpp       | 13 +++--
 .../pyinterp/detail/interpolation/steffen.hpp |  8 +--
 .../include/pyinterp/detail/isviewstream.hpp  |  4 +-
 .../core/include/pyinterp/detail/math.hpp     |  8 +--
 .../include/pyinterp/detail/math/bicubic.hpp  |  4 +-
 .../detail/math/radial_basis_functions.hpp    | 16 +++---
 .../include/pyinterp/detail/math/spline2d.hpp |  4 +-
 .../pyinterp/detail/math/trivariate.hpp       | 12 ++---
 .../pyinterp/detail/math/window_functions.hpp |  4 +-
 src/pyinterp/core/include/pyinterp/fill.hpp   | 16 +++---
 src/pyinterp/core/include/pyinterp/frame.hpp  |  4 +-
 .../include/pyinterp/geodetic/algorithm.hpp   | 18 +++----
 .../include/pyinterp/geodetic/line_string.hpp | 14 +++--
 .../core/include/pyinterp/geodetic/point.hpp  | 12 ++---
 .../core/include/pyinterp/geohash/int64.hpp   | 12 ++---
 .../core/include/pyinterp/geohash/string.hpp  | 35 ++++++-------
 src/pyinterp/core/include/pyinterp/period.hpp | 14 +++--
 .../core/include/pyinterp/quadrivariate.hpp   |  4 +-
 src/pyinterp/core/include/pyinterp/rtree.hpp  | 51 +++++++++----------
 .../include/pyinterp/streaming_histogram.hpp  |  4 +-
 .../core/include/pyinterp/temporal_axis.hpp   | 13 +++--
 .../core/include/pyinterp/trivariate.hpp      |  4 +-
 .../core/lib/geodetic/line_string.cpp         | 12 ++---
 src/pyinterp/core/lib/geodetic/rtree.cpp      | 20 ++++----
 src/pyinterp/core/lib/geohash/int64.cpp       | 16 +++---
 src/pyinterp/core/lib/geohash/string.cpp      | 31 +++++------
 src/pyinterp/core/module/bicubic.cpp          |  8 +--
 src/pyinterp/orbit.py                         |  4 +-
 .../statistics/descriptive_descriptive.py     |  3 +-
 .../statistics/streaming_histogram.py         |  3 +-
 src/pyinterp/tests/__init__.py                |  2 +-
 src/pyinterp/tests/core/test_dateutils.py     |  2 +-
 src/pyinterp/tests/core/test_geodetic.py      |  2 +-
 .../tests/test_orbit_interpolation.py         |  4 +-
 54 files changed, 261 insertions(+), 288 deletions(-)

diff --git a/CODE_OF_CONDUCT.md b/CODE_OF_CONDUCT.md
index d267279e..653ea211 100644
--- a/CODE_OF_CONDUCT.md
+++ b/CODE_OF_CONDUCT.md
@@ -6,7 +6,7 @@ In the interest of fostering an open and welcoming environment, we as
 contributors and maintainers pledge to making participation in our project and
 our community a harassment-free experience for everyone, regardless of age, body
 size, disability, ethnicity, sex characteristics, gender identity and expression,
-level of experience, education, socio-economic status, nationality, personal
+level of experience, education, socioeconomic status, nationality, personal
 appearance, race, religion, or sexual identity and orientation.
 
 ## Our Standards
diff --git a/scripts/stubgen.py b/scripts/stubgen.py
index dbd0b49b..281c4958 100644
--- a/scripts/stubgen.py
+++ b/scripts/stubgen.py
@@ -2,7 +2,6 @@
 Generates stubs for the core modules.
 =====================================
 """
-from typing import List
 import os
 import pathlib
 import re
@@ -62,7 +61,7 @@ def fix_core_geodetic(src: pathlib.Path):
         stream.writelines(lines)
 
 
-def fix_core_fill(src: pathlib.Path, grids: List[str]):
+def fix_core_fill(src: pathlib.Path, grids: list[str]):
     core = src / 'pyinterp' / 'core' / 'fill.pyi'
     with core.open('r') as stream:
         lines = stream.readlines()
diff --git a/src/pyinterp/backends/xarray.py b/src/pyinterp/backends/xarray.py
index f9a44022..9f1b4b6a 100644
--- a/src/pyinterp/backends/xarray.py
+++ b/src/pyinterp/backends/xarray.py
@@ -10,7 +10,7 @@
 """
 from __future__ import annotations
 
-from typing import Hashable
+from collections.abc import Hashable
 import pickle
 
 import numpy
diff --git a/src/pyinterp/core/__init__.pyi b/src/pyinterp/core/__init__.pyi
index dc7ea957..7034dceb 100644
--- a/src/pyinterp/core/__init__.pyi
+++ b/src/pyinterp/core/__init__.pyi
@@ -1224,7 +1224,7 @@ class Period:
 
 
 Array1DPeriod = numpy.ndarray[
-    Literal['N'],
+    tuple[int],
     numpy.dtype[Period]  # type: ignore[type-var]
 ]
 
diff --git a/src/pyinterp/core/array.pyi b/src/pyinterp/core/array.pyi
index 48a6235d..3b1bd80b 100644
--- a/src/pyinterp/core/array.pyi
+++ b/src/pyinterp/core/array.pyi
@@ -2,10 +2,10 @@ from typing import Any, Literal
 
 import numpy
 
-_1D = Literal['N']
-_2D = Literal['N', 'M']
-_3D = Literal['N', 'M', 'O']
-_4D = Literal['N', 'M', 'O', 'P']
+_1D = tuple[int]
+_2D = tuple[int, int]
+_3D = tuple[int, int, int]
+_4D = tuple[int, int, int, int]
 
 Array1DBool = numpy.ndarray[_1D, numpy.dtype[numpy.bool_]]
 Array1DFloat32 = numpy.ndarray[_1D, numpy.dtype[numpy.float32]]
diff --git a/src/pyinterp/core/include/pyinterp/axis.hpp b/src/pyinterp/core/include/pyinterp/axis.hpp
index b846c036..c887bd30 100644
--- a/src/pyinterp/core/include/pyinterp/axis.hpp
+++ b/src/pyinterp/core/include/pyinterp/axis.hpp
@@ -123,8 +123,8 @@ class Axis : public detail::Axis<T>,
   ///   points on the axis or the value -1 if the *bounded* parameter is set
   ///   to false and the index looked for is located outside the limits of the
   ///   axis.
-  auto find_index(const pybind11::array_t<T> &coordinates, bool bounded) const
-      -> pybind11::array_t<int64_t> {
+  auto find_index(const pybind11::array_t<T> &coordinates,
+                  bool bounded) const -> pybind11::array_t<int64_t> {
     detail::check_array_ndim("coordinates", 1, coordinates);
 
     auto size = coordinates.size();
diff --git a/src/pyinterp/core/include/pyinterp/bivariate.hpp b/src/pyinterp/core/include/pyinterp/bivariate.hpp
index bc36a5a5..3cae9bdc 100644
--- a/src/pyinterp/core/include/pyinterp/bivariate.hpp
+++ b/src/pyinterp/core/include/pyinterp/bivariate.hpp
@@ -112,8 +112,8 @@ template <template <class> class Point, typename Coordinate, typename Type>
 auto bivariate(const Grid2D<Type> &grid, const pybind11::array_t<Coordinate> &x,
                const pybind11::array_t<Coordinate> &y,
                const BivariateInterpolator<Point, Coordinate> *interpolator,
-               const bool bounds_error, const size_t num_threads)
-    -> pybind11::array_t<Coordinate> {
+               const bool bounds_error,
+               const size_t num_threads) -> pybind11::array_t<Coordinate> {
   pyinterp::detail::check_array_ndim("x", 1, x, "y", 1, y);
   pyinterp::detail::check_ndarray_shape("x", x, "y", y);
 
@@ -167,8 +167,8 @@ void implement_bivariate_interpolator(pybind11::module &m,
     using CoordinateSystem::CoordinateSystem;
 
     auto evaluate(const Point<T> &p, const Point<T> &p0, const Point<T> &p1,
-                  const T &q00, const T &q01, const T &q10, const T &q11) const
-        -> T override {
+                  const T &q00, const T &q01, const T &q10,
+                  const T &q11) const -> T override {
       PYBIND11_OVERLOAD_PURE(T, CoordinateSystem, "evaluate", p, p0,  // NOLINT
                              p1,                                      // NOLINT
                              q00, q01, q10, q11);                     // NOLINT
diff --git a/src/pyinterp/core/include/pyinterp/dateutils.hpp b/src/pyinterp/core/include/pyinterp/dateutils.hpp
index d8d171c2..f190dbb6 100644
--- a/src/pyinterp/core/include/pyinterp/dateutils.hpp
+++ b/src/pyinterp/core/include/pyinterp/dateutils.hpp
@@ -310,9 +310,8 @@ class FractionalSeconds {
     return static_cast<int>(std::log10(order_of_magnitude_));
   }
 
-  [[nodiscard]] constexpr auto cast(const int64_t frac,
-                                    const int64_t scale) const noexcept
-      -> int64_t {
+  [[nodiscard]] constexpr auto cast(
+      const int64_t frac, const int64_t scale) const noexcept -> int64_t {
     return order_of_magnitude_ <= scale ? (scale / order_of_magnitude_) * frac
                                         : frac / (order_of_magnitude_ / scale);
   }
@@ -473,8 +472,8 @@ constexpr auto datetime_from_minutes(const int64_t minutes)
 }
 
 /// Convert a datetime64 to a string
-inline auto datetime64_to_string(const int64_t value, const DType &dtype)
-    -> std::string {
+inline auto datetime64_to_string(const int64_t value,
+                                 const DType &dtype) -> std::string {
   auto date = Date{};
   auto time = Time{};
   auto ss = std::stringstream{};
diff --git a/src/pyinterp/core/include/pyinterp/detail/axis.hpp b/src/pyinterp/core/include/pyinterp/detail/axis.hpp
index 7e5e0a2a..64c77c46 100644
--- a/src/pyinterp/core/include/pyinterp/detail/axis.hpp
+++ b/src/pyinterp/core/include/pyinterp/detail/axis.hpp
@@ -144,8 +144,8 @@ class Axis {
   /// @param start index of the first element to include in the slice
   /// @param count number of elements to include in the slice
   /// @return a slice of the axis
-  [[nodiscard]] constexpr auto slice(int64_t start, int64_t count) const
-      -> Vector<T> {
+  [[nodiscard]] constexpr auto slice(int64_t start,
+                                     int64_t count) const -> Vector<T> {
     if (start < 0 || start + count > size()) {
       throw std::out_of_range("axis index out of range");
     }
@@ -248,9 +248,8 @@ class Axis {
 
   /// Returns the normalized value of the coordinate with the respect to the
   /// axis definition.
-  [[nodiscard]] constexpr auto normalize_coordinate(const T coordinate,
-                                                    const T min) const noexcept
-      -> T {
+  [[nodiscard]] constexpr auto normalize_coordinate(
+      const T coordinate, const T min) const noexcept -> T {
     if (is_angle() && (coordinate >= min + circle_ || coordinate < min)) {
       return math::normalize_angle(coordinate, min, circle_);
     }
@@ -290,9 +289,8 @@ class Axis {
   /// is located before, or the value of the last element of this container if
   /// the requested value is located after.
   /// @return index of the grid point containing it or -1 if outside grid area
-  [[nodiscard]] constexpr auto find_nearest_index(T coordinate,
-                                                  const bool bounded) const
-      -> int64_t {
+  [[nodiscard]] constexpr auto find_nearest_index(
+      T coordinate, const bool bounded) const -> int64_t {
     coordinate = normalize_coordinate(coordinate);
     auto result = axis_->find_index(coordinate, bounded);
     if (result == -1 && is_circle_) {
diff --git a/src/pyinterp/core/include/pyinterp/detail/axis/container.hpp b/src/pyinterp/core/include/pyinterp/detail/axis/container.hpp
index 9a7301b1..73f1fc5a 100644
--- a/src/pyinterp/core/include/pyinterp/detail/axis/container.hpp
+++ b/src/pyinterp/core/include/pyinterp/detail/axis/container.hpp
@@ -70,8 +70,8 @@ class Abstract {
   /// @param start index of the first element to include in the slice
   /// @param count number of elements to include in the slice
   /// @return a slice of the axis
-  [[nodiscard]] virtual auto slice(int64_t start, int64_t count) const
-      -> Vector<T> = 0;
+  [[nodiscard]] virtual auto slice(int64_t start,
+                                   int64_t count) const -> Vector<T> = 0;
 
   /// Get the minimum coordinate value.
   ///
@@ -107,8 +107,8 @@ class Abstract {
   /// the requested value is located after.
   /// @return index of the requested value it or -1 if outside this coordinate
   /// system area.
-  [[nodiscard]] virtual auto find_index(T coordinate, bool bounded) const
-      -> int64_t = 0;
+  [[nodiscard]] virtual auto find_index(T coordinate,
+                                        bool bounded) const -> int64_t = 0;
 
   /// compare two variables instances
   ///
@@ -281,9 +281,8 @@ class Irregular : public Abstract<T> {
   }
 
   /// @copydoc Abstract::slice(const int64_t, const int64_t) const
-  [[nodiscard]] constexpr auto slice(const int64_t start,
-                                     const int64_t count) const noexcept
-      -> Vector<T> override {
+  [[nodiscard]] constexpr auto slice(const int64_t start, const int64_t count)
+      const noexcept -> Vector<T> override {
     return points_.segment(start, count);
   }
 
@@ -313,9 +312,8 @@ class Irregular : public Abstract<T> {
   }
 
   /// @copydoc Abstract::find_index(double,bool) const
-  [[nodiscard]] constexpr auto find_index(const T coordinate,
-                                          const bool bounded) const
-      -> int64_t override {
+  [[nodiscard]] constexpr auto find_index(
+      const T coordinate, const bool bounded) const -> int64_t override {
     if (this->is_ascending_) {
       return this->find_index(coordinate, bounded, size(), std::less<T>());
     }
@@ -435,9 +433,8 @@ class AbstractRegular : public Abstract<T> {
   }
 
   /// @copydoc Abstract::slice(const int64_t, const int64_t) const
-  [[nodiscard]] constexpr auto slice(const int64_t start,
-                                     const int64_t count) const noexcept
-      -> Vector<T> override {
+  [[nodiscard]] constexpr auto slice(const int64_t start, const int64_t count)
+      const noexcept -> Vector<T> override {
     auto result = Vector<T>(count);
     for (int64_t ix = 0; ix < count; ++ix) {
       result[ix] = coordinate_value(start + ix);
@@ -551,9 +548,8 @@ class Regular<T, typename std::enable_if<std::is_integral_v<T>>::type>
       : AbstractRegular<T>(start, stop, num), step_2_(this->step_ >> 1) {}
 
   /// @copydoc Abstract::find_index(double,bool) const
-  [[nodiscard]] constexpr auto find_index(T coordinate,
-                                          bool bounded) const noexcept
-      -> int64_t override {
+  [[nodiscard]] constexpr auto find_index(
+      T coordinate, bool bounded) const noexcept -> int64_t override {
     auto index =
         static_cast<int64_t>(round((coordinate - this->start_), this->step_));
 
diff --git a/src/pyinterp/core/include/pyinterp/detail/geometry/crossover.hpp b/src/pyinterp/core/include/pyinterp/detail/geometry/crossover.hpp
index 0fdb2bfe..1c9a7c97 100644
--- a/src/pyinterp/core/include/pyinterp/detail/geometry/crossover.hpp
+++ b/src/pyinterp/core/include/pyinterp/detail/geometry/crossover.hpp
@@ -36,8 +36,8 @@ constexpr auto is_retrograde(const T& x1, const T& x2, const T& x3) -> bool {
 /// @param point the point to search.
 /// @return the index of the nearest point.
 template <typename T>
-auto nearest_point(const LineString<T>& line, const Point2D<T>& point)
-    -> size_t {
+auto nearest_point(const LineString<T>& line,
+                   const Point2D<T>& point) -> size_t {
   auto min_distance = std::numeric_limits<T>::max();
   auto index = size_t(0);
   for (auto it = line.begin(); it != line.end(); ++it) {
diff --git a/src/pyinterp/core/include/pyinterp/detail/geometry/rtree.hpp b/src/pyinterp/core/include/pyinterp/detail/geometry/rtree.hpp
index 98dd7a6c..132a1356 100644
--- a/src/pyinterp/core/include/pyinterp/detail/geometry/rtree.hpp
+++ b/src/pyinterp/core/include/pyinterp/detail/geometry/rtree.hpp
@@ -131,8 +131,8 @@ class RTree {
   /// @overload query(const Point &, const Strategy &, const uint32_t) const
   ///
   /// Overload of the query method with the default strategy.
-  auto query(const Point &point, const uint32_t k) const
-      -> std::vector<result_t> {
+  auto query(const Point &point,
+             const uint32_t k) const -> std::vector<result_t> {
     return query(point, boost::geometry::default_strategy(), k);
   }
 
@@ -162,8 +162,8 @@ class RTree {
   /// @overload query_ball(const Point &, const Strategy &, const coordinate_t)
   ///
   /// Overload of the query_ball method with the default strategy.
-  auto query_ball(const Point &point, const coordinate_t radius) const
-      -> std::vector<result_t> {
+  auto query_ball(const Point &point,
+                  const coordinate_t radius) const -> std::vector<result_t> {
     return query_ball(point, boost::geometry::default_strategy(), radius);
   }
 
@@ -203,8 +203,8 @@ class RTree {
   /// @overload query_within(const Point &, const Strategy &, const uint32_t)
   ///
   /// Overload of the query_within method with the default strategy.
-  auto query_within(const Point &point, const uint32_t k) const
-      -> std::vector<result_t> {
+  auto query_within(const Point &point,
+                    const uint32_t k) const -> std::vector<result_t> {
     return query_within(point, boost::geometry::default_strategy(), k);
   }
 
@@ -259,8 +259,8 @@ class RTree {
   ///
   /// Overload of the value method with the default strategy.
   auto value(const Point &point, const std::optional<coordinate_t> &radius,
-             const uint32_t k, const bool within) const
-      -> std::vector<value_t> {
+             const uint32_t k,
+             const bool within) const -> std::vector<value_t> {
     return value(point, boost::geometry::default_strategy(), radius, k, within);
   }
 
@@ -516,11 +516,10 @@ class RTree {
   /// math::RBF<promotion_t> &, const coordinate_t, const uint32_t, const bool)
   ///
   /// Overload of the radial_basis_function method with the default strategy.
-  auto radial_basis_function(const Point &point,
-                             const math::RBF<promotion_t> &rbf,
-                             const coordinate_t radius, const uint32_t k,
-                             const bool within) const
-      -> std::pair<promotion_t, uint32_t> {
+  auto radial_basis_function(
+      const Point &point, const math::RBF<promotion_t> &rbf,
+      const coordinate_t radius, const uint32_t k,
+      const bool within) const -> std::pair<promotion_t, uint32_t> {
     return radial_basis_function(point, boost::geometry::default_strategy(),
                                  rbf, radius, k, within);
   }
@@ -582,11 +581,10 @@ class RTree {
   /// coordinate_t, const uint32_t, const bool)
   ///
   /// Overload of the window_function method with the default strategy.
-  auto window_function(const Point &point,
-                       const math::WindowFunction<coordinate_t> &wf,
-                       const coordinate_t arg, const coordinate_t radius,
-                       const uint32_t k, const bool within) const
-      -> std::pair<coordinate_t, uint32_t> {
+  auto window_function(
+      const Point &point, const math::WindowFunction<coordinate_t> &wf,
+      const coordinate_t arg, const coordinate_t radius, const uint32_t k,
+      const bool within) const -> std::pair<coordinate_t, uint32_t> {
     return window_function(point, boost::geometry::default_strategy(), wf, arg,
                            radius, k, within);
   }
diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/akima.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/akima.hpp
index bca47a61..b44d5940 100644
--- a/src/pyinterp/core/include/pyinterp/detail/interpolation/akima.hpp
+++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/akima.hpp
@@ -33,8 +33,8 @@ class Akima : public Interpolator1D<T> {
   /// @brief Compute the coefficients of the interpolation
   /// @param xa X-coordinates of the data points.
   /// @param ya Y-coordinates of the data points.
-  auto compute_coefficients(const Vector<T>& xa, const Vector<T>& ya)
-      -> void override;
+  auto compute_coefficients(const Vector<T>& xa,
+                            const Vector<T>& ya) -> void override;
 
   /// Interpolation
   /// @param xa X-coordinates of the data points.
@@ -55,8 +55,8 @@ class Akima : public Interpolator1D<T> {
 };
 
 template <typename T>
-auto Akima<T>::compute_coefficients(const Vector<T>& xa, const Vector<T>& ya)
-    -> void {
+auto Akima<T>::compute_coefficients(const Vector<T>& xa,
+                                    const Vector<T>& ya) -> void {
   Interpolator1D<T>::compute_coefficients(xa, ya);
   auto size = xa.size();
   if (m_.size() < size + 4) {
diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/bicubic.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/bicubic.hpp
index 12f7ef95..8c93f9f8 100644
--- a/src/pyinterp/core/include/pyinterp/detail/interpolation/bicubic.hpp
+++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/bicubic.hpp
@@ -29,8 +29,8 @@ class Bicubic : public Interpolator2D<T> {
   /// @param y The point where the interpolation must be calculated.
   /// @return The interpolated value at the coordinates x, y.
   constexpr auto interpolate_(const Vector<T> &xa, const Vector<T> &ya,
-                              const Matrix<T> &za, const T &x, const T &y) const
-      -> T override;
+                              const Matrix<T> &za, const T &x,
+                              const T &y) const -> T override;
 
   /// Compute the coefficients of the bicubic interpolation
   /// @param xa X-coordinates of the data points.
diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/bilinear.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/bilinear.hpp
index 80fdc437..ec08f494 100644
--- a/src/pyinterp/core/include/pyinterp/detail/interpolation/bilinear.hpp
+++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/bilinear.hpp
@@ -26,8 +26,8 @@ class Bilinear : public Interpolator2D<T> {
   /// @param x The point where the interpolation must be calculated.
   /// @param y The point where the interpolation must be calculated.
   constexpr auto interpolate_(const Vector<T> &xa, const Vector<T> &ya,
-                              const Matrix<T> &za, const T &x, const T &y) const
-      -> T override;
+                              const Matrix<T> &za, const T &x,
+                              const T &y) const -> T override;
 };
 
 template <typename T>
diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/cspline.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/cspline.hpp
index d5c79c6d..c3c81fe6 100644
--- a/src/pyinterp/core/include/pyinterp/detail/interpolation/cspline.hpp
+++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/cspline.hpp
@@ -24,8 +24,8 @@ class CSpline : public CSplineBase<T> {
   /// @brief Compute the coefficients of the interpolation
   /// @param xa X-coordinates of the data points.
   /// @param ya Y-coordinates of the data points.
-  inline auto compute_coefficients(const Vector<T> &xa, const Vector<T> &ya)
-      -> void override;
+  inline auto compute_coefficients(const Vector<T> &xa,
+                                   const Vector<T> &ya) -> void override;
 
   /// @brief Solve a symmetric tridiagonal system
   /// @param x The solution of the system
@@ -63,8 +63,8 @@ constexpr auto CSpline<T>::solve_symmetric_tridiagonal(T *x) -> void {
 }
 
 template <typename T>
-auto CSpline<T>::compute_coefficients(const Vector<T> &xa, const Vector<T> &ya)
-    -> void {
+auto CSpline<T>::compute_coefficients(const Vector<T> &xa,
+                                      const Vector<T> &ya) -> void {
   Interpolator1D<T>::compute_coefficients(xa, ya);
   const auto size = xa.size();
   const auto size_m2 = size - 2;
diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/cspline_base.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/cspline_base.hpp
index 29c0b458..fa22dd46 100644
--- a/src/pyinterp/core/include/pyinterp/detail/interpolation/cspline_base.hpp
+++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/cspline_base.hpp
@@ -16,14 +16,14 @@ class CSplineCoefficients {
   /// Constructor
   /// @param c0 The first derivative at the first point
   /// @param c1 The first derivative at the last point
-  constexpr CSplineCoefficients(const T &c0, const T &c1) : c0_(c0), c1_(c1){};
+  constexpr CSplineCoefficients(const T &c0, const T &c1) : c0_(c0), c1_(c1) {};
 
   /// Compute the coefficients of the cubic spline interpolation
   /// @param dx The distance between the two points
   /// @param dy The difference between the two points
   /// @return The coefficients of the cubic spline interpolation
-  constexpr auto operator()(const T &dx, const T &dy) const
-      -> std::tuple<T, T, T> {
+  constexpr auto operator()(const T &dx,
+                            const T &dy) const -> std::tuple<T, T, T> {
     constexpr auto third = T(1) / T(3);
     return {(dy / dx) - dx * (c1_ + 2 * c0_) * third, c0_,
             (c1_ - c0_) / (3 * dx)};
diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/cspline_periodic.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/cspline_periodic.hpp
index 28479616..db4b2466 100644
--- a/src/pyinterp/core/include/pyinterp/detail/interpolation/cspline_periodic.hpp
+++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/cspline_periodic.hpp
@@ -22,8 +22,8 @@ class CSplinePeriodic : public CSplineBase<T> {
   /// Compute the coefficients of the interpolation
   /// @param xa X-coordinates of the data points.
   /// @param ya Y-coordinates of the data points.
-  auto compute_coefficients(const Vector<T> &xa, const Vector<T> &ya)
-      -> void override;
+  auto compute_coefficients(const Vector<T> &xa,
+                            const Vector<T> &ya) -> void override;
 
   /// Solve a symmetric cyclic tridiagonal system
   /// @param x The solution of the system
diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/interpolator_1d.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/interpolator_1d.hpp
index 3847470b..5cb25a0a 100644
--- a/src/pyinterp/core/include/pyinterp/detail/interpolation/interpolator_1d.hpp
+++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/interpolator_1d.hpp
@@ -30,8 +30,8 @@ class Interpolator1D : public Interpolator<T> {
   /// @param ya Y-coordinates of the data points.
   /// @param x The points where the interpolation must be calculated.
   /// @return The interpolated values at the points x.
-  auto operator()(const Vector<T> &xa, const Vector<T> &ya, const Vector<T> &x)
-      -> Vector<T> {
+  auto operator()(const Vector<T> &xa, const Vector<T> &ya,
+                  const Vector<T> &x) -> Vector<T> {
     compute_coefficients(xa, ya);
     auto y = Vector<T>(x.size());
     for (Eigen::Index i = 0; i < x.size(); ++i) {
@@ -55,8 +55,8 @@ class Interpolator1D : public Interpolator<T> {
   /// @param ya Y-coordinates of the data points.
   /// @param x The points where the derivative must be calculated.
   /// @return The derivatives of the interpolation function at the points x.
-  auto derivative(const Vector<T> &xa, const Vector<T> &ya, const Vector<T> &x)
-      -> Vector<T> {
+  auto derivative(const Vector<T> &xa, const Vector<T> &ya,
+                  const Vector<T> &x) -> Vector<T> {
     compute_coefficients(xa, ya);
     auto y = Vector<T>(x.size());
     for (Eigen::Index i = 0; i < x.size(); ++i) {
@@ -75,8 +75,8 @@ class Interpolator1D : public Interpolator<T> {
                            const T &x) const -> T = 0;
 
   /// Check if the arrays are valid.
-  virtual auto compute_coefficients(const Vector<T> &xa, const Vector<T> &ya)
-      -> void {
+  virtual auto compute_coefficients(const Vector<T> &xa,
+                                    const Vector<T> &ya) -> void {
     if (xa.size() != ya.size()) {
       throw std::invalid_argument("xa and ya must have the same size");
     }
diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/interpolator_2d.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/interpolator_2d.hpp
index ddfe2e85..0fa77c25 100644
--- a/src/pyinterp/core/include/pyinterp/detail/interpolation/interpolator_2d.hpp
+++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/interpolator_2d.hpp
@@ -49,8 +49,8 @@ class Interpolator2D : public Interpolator<T> {
  protected:
   /// Interpolate the value of y at x using the index of the last search.
   virtual auto interpolate_(const Vector<T> &xa, const Vector<T> &ya,
-                            const Matrix<T> &za, const T &x, const T &y) const
-      -> T = 0;
+                            const Matrix<T> &za, const T &x,
+                            const T &y) const -> T = 0;
 
   /// Check if the arrays are valid.
   virtual auto compute_coefficients(const Vector<T> &xa, const Vector<T> &ya,
diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/polynomial.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/polynomial.hpp
index 4fea5eaf..a40ff8ec 100644
--- a/src/pyinterp/core/include/pyinterp/detail/interpolation/polynomial.hpp
+++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/polynomial.hpp
@@ -25,8 +25,8 @@ class Polynomial : public Interpolator1D<T> {
   /// Compute the coefficients of the interpolation
   /// @param xa X-coordinates of the data points.
   /// @param ya Y-coordinates of the data points.
-  constexpr auto compute_coefficients(const Vector<T> &xa, const Vector<T> &ya)
-      -> void override;
+  constexpr auto compute_coefficients(const Vector<T> &xa,
+                                      const Vector<T> &ya) -> void override;
 
   /// Compute the coefficients of the interpolation
   /// @param xa X-coordinates of the data points.
@@ -47,14 +47,13 @@ class Polynomial : public Interpolator1D<T> {
   /// @param ya Y-coordinates of the data points.
   /// @param x The point where the derivative must be calculated.
   /// @return The derivative of the interpolation function at the point x.
-  auto derivative_(const Vector<T> &xa, const Vector<T> &ya, const T &x) const
-      -> T override;
+  auto derivative_(const Vector<T> &xa, const Vector<T> &ya,
+                   const T &x) const -> T override;
 };
 
 template <typename T>
-constexpr auto Polynomial<T>::compute_coefficients(const Vector<T> &xa,
-                                                   const Vector<T> &ya)
-    -> void {
+constexpr auto Polynomial<T>::compute_coefficients(
+    const Vector<T> &xa, const Vector<T> &ya) -> void {
   Interpolator1D<T>::compute_coefficients(xa, ya);
   auto size = xa.size();
   if (work_.size() < size) {
diff --git a/src/pyinterp/core/include/pyinterp/detail/interpolation/steffen.hpp b/src/pyinterp/core/include/pyinterp/detail/interpolation/steffen.hpp
index d4c7137c..d3219fc2 100644
--- a/src/pyinterp/core/include/pyinterp/detail/interpolation/steffen.hpp
+++ b/src/pyinterp/core/include/pyinterp/detail/interpolation/steffen.hpp
@@ -39,8 +39,8 @@ class Steffen : public Interpolator1D<T> {
   /// Compute the coefficients of the interpolation
   /// @param xa X-coordinates of the data points.
   /// @param ya Y-coordinates of the data points.
-  constexpr auto compute_coefficients(const Vector<T> &xa, const Vector<T> &ya)
-      -> void override;
+  constexpr auto compute_coefficients(const Vector<T> &xa,
+                                      const Vector<T> &ya) -> void override;
 
   /// Return the sign of x multiplied by the sign of y
   static constexpr auto copysign(const T &x, const T &y) -> T {
@@ -91,8 +91,8 @@ constexpr auto Steffen<T>::compute_coefficients(const Vector<T> &xa,
 
 template <typename T>
 constexpr auto Steffen<T>::interpolate_(const Vector<T> &xa,
-                                        const Vector<T> &ya, const T &x) const
-    -> T {
+                                        const Vector<T> &ya,
+                                        const T &x) const -> T {
   auto where = this->search(xa, x);
   if (!where) {
     return std::numeric_limits<T>::quiet_NaN();
diff --git a/src/pyinterp/core/include/pyinterp/detail/isviewstream.hpp b/src/pyinterp/core/include/pyinterp/detail/isviewstream.hpp
index 9c689bdd..74b2c4e1 100644
--- a/src/pyinterp/core/include/pyinterp/detail/isviewstream.hpp
+++ b/src/pyinterp/core/include/pyinterp/detail/isviewstream.hpp
@@ -34,8 +34,8 @@ class sviewbuf : public std::streambuf {
     return gptr() - eback();
   }
 
-  auto seekpos(pos_type sp, std::ios_base::openmode which)
-      -> pos_type override {
+  auto seekpos(pos_type sp,
+               std::ios_base::openmode which) -> pos_type override {
     return seekoff(sp - pos_type(static_cast<off_type>(0)), std::ios_base::beg,
                    which);
   }
diff --git a/src/pyinterp/core/include/pyinterp/detail/math.hpp b/src/pyinterp/core/include/pyinterp/detail/math.hpp
index a233c4b7..cd0cc1ba 100644
--- a/src/pyinterp/core/include/pyinterp/detail/math.hpp
+++ b/src/pyinterp/core/include/pyinterp/detail/math.hpp
@@ -269,16 +269,16 @@ constexpr auto sinc(const T &x) noexcept -> T {
 /// True if a and b are two values identical to an epsilon.
 template <typename T, typename std::enable_if<std::is_integral<T>::value,
                                               T>::type * = nullptr>
-constexpr auto is_same(const T &a, const T &b, const T &epsilon) noexcept
-    -> bool {
+constexpr auto is_same(const T &a, const T &b,
+                       const T &epsilon) noexcept -> bool {
   return std::abs(a - b) <= epsilon;
 }
 
 /// True if a and b are two values identical to an epsilon.
 template <typename T, typename std::enable_if<std::is_floating_point<T>::value,
                                               T>::type * = nullptr>
-constexpr auto is_same(const T &a, const T &b, const T &epsilon) noexcept
-    -> bool {
+constexpr auto is_same(const T &a, const T &b,
+                       const T &epsilon) noexcept -> bool {
   auto diff = std::fabs(a - b);
   if (diff <= epsilon) {
     return true;
diff --git a/src/pyinterp/core/include/pyinterp/detail/math/bicubic.hpp b/src/pyinterp/core/include/pyinterp/detail/math/bicubic.hpp
index 6022cd30..8642cdca 100644
--- a/src/pyinterp/core/include/pyinterp/detail/math/bicubic.hpp
+++ b/src/pyinterp/core/include/pyinterp/detail/math/bicubic.hpp
@@ -22,8 +22,8 @@ class Bicubic {
       : interpolator_(interpolation::factory_2d<double>(kind)) {}
 
   /// Return the interpolated value of y for a given point x
-  auto interpolate(const double x, const double y, const Frame2D &xr)
-      -> double {
+  auto interpolate(const double x, const double y,
+                   const Frame2D &xr) -> double {
     return (*interpolator_)(*(xr.x()), *(xr.y()), *(xr.q()), x, y);
   }
 
diff --git a/src/pyinterp/core/include/pyinterp/detail/math/radial_basis_functions.hpp b/src/pyinterp/core/include/pyinterp/detail/math/radial_basis_functions.hpp
index d4e09f11..6c533c91 100644
--- a/src/pyinterp/core/include/pyinterp/detail/math/radial_basis_functions.hpp
+++ b/src/pyinterp/core/include/pyinterp/detail/math/radial_basis_functions.hpp
@@ -143,20 +143,20 @@ class RBF {
   }
 
   /// Gauss
-  static auto gaussian(const Eigen::Ref<const Matrix<T>> &r, const T epsilon)
-      -> Matrix<T> {
+  static auto gaussian(const Eigen::Ref<const Matrix<T>> &r,
+                       const T epsilon) -> Matrix<T> {
     return (-(epsilon * r).array().pow(2)).exp();
   }
 
   /// Linear spline
-  static auto linear(const Eigen::Ref<const Matrix<T>> &r, const T /*epsilon*/)
-      -> Matrix<T> {
+  static auto linear(const Eigen::Ref<const Matrix<T>> &r,
+                     const T /*epsilon*/) -> Matrix<T> {
     return r;
   }
 
   /// Cubic spline
-  static auto cubic(const Eigen::Ref<const Matrix<T>> &r, const T /*epsilon*/)
-      -> Matrix<T> {
+  static auto cubic(const Eigen::Ref<const Matrix<T>> &r,
+                    const T /*epsilon*/) -> Matrix<T> {
     return r.array().pow(3);
   }
 
@@ -183,8 +183,8 @@ class RBF {
   }
 
   /// Resolution of the linear system
-  static auto solve_linear_system(const Matrix<T> &A, const Vector<T> &di)
-      -> Vector<T> {
+  static auto solve_linear_system(const Matrix<T> &A,
+                                  const Vector<T> &di) -> Vector<T> {
     Eigen::FullPivLU<Matrix<T>> lu(A);
     return lu.solve(di);
   }
diff --git a/src/pyinterp/core/include/pyinterp/detail/math/spline2d.hpp b/src/pyinterp/core/include/pyinterp/detail/math/spline2d.hpp
index 563d05a0..90af18de 100644
--- a/src/pyinterp/core/include/pyinterp/detail/math/spline2d.hpp
+++ b/src/pyinterp/core/include/pyinterp/detail/math/spline2d.hpp
@@ -24,8 +24,8 @@ class Spline2D {
         y_interpolator_(interpolation::factory_1d<double>(kind)) {}
 
   /// Return the interpolated value of y for a given point x
-  auto interpolate(const double x, const double y, const Frame2D &xr)
-      -> double {
+  auto interpolate(const double x, const double y,
+                   const Frame2D &xr) -> double {
     // Spline interpolation as function of X-coordinate
     for (Eigen::Index ix = 0; ix < xr.y()->size(); ++ix) {
       column_(ix) = (*x_interpolator_)(*(xr.x()), xr.q()->col(ix), x);
diff --git a/src/pyinterp/core/include/pyinterp/detail/math/trivariate.hpp b/src/pyinterp/core/include/pyinterp/detail/math/trivariate.hpp
index 53b61a3d..9e3265cc 100644
--- a/src/pyinterp/core/include/pyinterp/detail/math/trivariate.hpp
+++ b/src/pyinterp/core/include/pyinterp/detail/math/trivariate.hpp
@@ -109,13 +109,11 @@ constexpr auto trivariate(
 /// @param q111 Point value for the coordinate (x1, y1, z1)
 /// @return interpolated value at coordinate (x, y, z)
 template <template <class> class Point, typename T>
-constexpr auto trivariate(const Point<T> &p, const Point<T> &p0,
-                          const Point<T> &p1, const T &q000, const T &q010,
-                          const T &q100, const T &q110, const T &q001,
-                          const T &q011, const T &q101, const T &q111,
-                          const Bivariate<Point, T> *bivariate,
-                          const z_method_t<T, T> &interpolator = &linear<T, T>)
-    -> T {
+constexpr auto trivariate(
+    const Point<T> &p, const Point<T> &p0, const Point<T> &p1, const T &q000,
+    const T &q010, const T &q100, const T &q110, const T &q001, const T &q011,
+    const T &q101, const T &q111, const Bivariate<Point, T> *bivariate,
+    const z_method_t<T, T> &interpolator = &linear<T, T>) -> T {
   auto z0 = bivariate->evaluate(p, p0, p1, q000, q010, q100, q110);
   auto z1 = bivariate->evaluate(p, p0, p1, q001, q011, q101, q111);
   return interpolator(boost::geometry::get<2>(p), boost::geometry::get<2>(p0),
diff --git a/src/pyinterp/core/include/pyinterp/detail/math/window_functions.hpp b/src/pyinterp/core/include/pyinterp/detail/math/window_functions.hpp
index 2b9ec221..3cdc1745 100644
--- a/src/pyinterp/core/include/pyinterp/detail/math/window_functions.hpp
+++ b/src/pyinterp/core/include/pyinterp/detail/math/window_functions.hpp
@@ -74,8 +74,8 @@ constexpr auto nuttall(const T &d, const T &r, const T & /*unused*/) -> T {
 
 /// kBlackman-Harris window function.
 template <typename T>
-constexpr auto blackman_harris(const T &d, const T &r, const T & /*unused*/)
-    -> T {
+constexpr auto blackman_harris(const T &d, const T &r,
+                               const T & /*unused*/) -> T {
   if (d <= r) {
     auto ratio = (d + r) / r;
     return static_cast<T>(0.35875 - 0.48829 * std::cos(pi<T>() * ratio) +
diff --git a/src/pyinterp/core/include/pyinterp/fill.hpp b/src/pyinterp/core/include/pyinterp/fill.hpp
index d7c13054..88536cb0 100644
--- a/src/pyinterp/core/include/pyinterp/fill.hpp
+++ b/src/pyinterp/core/include/pyinterp/fill.hpp
@@ -297,8 +297,8 @@ template <typename Type>
 auto gauss_seidel(pybind11::EigenDRef<Matrix<Type>> &grid,
                   const FirstGuess first_guess, const bool is_circle,
                   const size_t max_iterations, const Type epsilon,
-                  const Type relaxation, size_t num_threads)
-    -> std::tuple<size_t, Type> {
+                  const Type relaxation,
+                  size_t num_threads) -> std::tuple<size_t, Type> {
   /// If the grid doesn't have an undefined value, this routine has nothing more
   /// to do.
   if (!grid.hasNaN()) {
@@ -343,8 +343,8 @@ auto gauss_seidel(pybind11::EigenDRef<Matrix<Type>> &grid,
 
 // Get the indexes that frame a given index.
 inline auto frame_index(const int64_t index, const int64_t size,
-                        const bool is_angle, std::vector<int64_t> &frame)
-    -> void {
+                        const bool is_angle,
+                        std::vector<int64_t> &frame) -> void {
   // Index in the center of the window
   auto center = static_cast<int64_t>(frame.size() / 2);
 
@@ -371,8 +371,8 @@ inline auto frame_index(const int64_t index, const int64_t size,
 }
 
 /// Checking the size of the filter window.
-constexpr auto check_windows_size(const std::string &name1, const uint32_t size)
-    -> void {
+constexpr auto check_windows_size(const std::string &name1,
+                                  const uint32_t size) -> void {
   if (size < 1) {
     throw std::invalid_argument(name1 + " must be >= 1");
   }
@@ -409,8 +409,8 @@ enum ValueType {
 /// @return The grid will have all the NaN filled with extrapolated values.
 template <typename Type>
 auto loess(const Grid2D<Type> &grid, const uint32_t nx, const uint32_t ny,
-           const ValueType value_type, const size_t num_threads)
-    -> pybind11::array_t<Type> {
+           const ValueType value_type,
+           const size_t num_threads) -> pybind11::array_t<Type> {
   check_windows_size("nx", nx, "ny", ny);
   auto result = pybind11::array_t<Type>(
       pybind11::array::ShapeContainer{grid.x()->size(), grid.y()->size()});
diff --git a/src/pyinterp/core/include/pyinterp/frame.hpp b/src/pyinterp/core/include/pyinterp/frame.hpp
index ec7ba048..3f3386c0 100644
--- a/src/pyinterp/core/include/pyinterp/frame.hpp
+++ b/src/pyinterp/core/include/pyinterp/frame.hpp
@@ -67,8 +67,8 @@ auto load_frame(const Grid2D<DataType> &grid, const double x, const double y,
 template <typename DataType, typename AxisType>
 auto load_frame(const Grid3D<DataType, AxisType> &grid, const double x,
                 const double y, const AxisType z, const axis::Boundary boundary,
-                const bool bounds_error, detail::math::Frame3D<AxisType> &frame)
-    -> bool {
+                const bool bounds_error,
+                detail::math::Frame3D<AxisType> &frame) -> bool {
   const auto &x_axis = *grid.x();
   const auto &y_axis = *grid.y();
   const auto &z_axis = *grid.z();
diff --git a/src/pyinterp/core/include/pyinterp/geodetic/algorithm.hpp b/src/pyinterp/core/include/pyinterp/geodetic/algorithm.hpp
index efbc60b1..d5f05b63 100644
--- a/src/pyinterp/core/include/pyinterp/geodetic/algorithm.hpp
+++ b/src/pyinterp/core/include/pyinterp/geodetic/algorithm.hpp
@@ -54,9 +54,8 @@ template <typename Geometry>
 
 /// Calculate the perimeter
 template <typename Geometry>
-[[nodiscard]] inline auto perimeter(const Geometry &geometry,
-                                    const std::optional<Spheroid> &wgs)
-    -> double {
+[[nodiscard]] inline auto perimeter(
+    const Geometry &geometry, const std::optional<Spheroid> &wgs) -> double {
   auto spheroid =
       wgs.has_value()
           ? static_cast<boost::geometry::srs::spheroid<double>>(*wgs)
@@ -70,8 +69,8 @@ template <typename Geometry>
 template <typename Geometry1, typename Geometry2>
 [[nodiscard]] inline auto covered_by(
     const Geometry2 &geometry2, const Eigen::Ref<const Eigen::VectorXd> &lon,
-    const Eigen::Ref<const Eigen::VectorXd> &lat, const size_t num_threads)
-    -> pybind11::array_t<bool> {
+    const Eigen::Ref<const Eigen::VectorXd> &lat,
+    const size_t num_threads) -> pybind11::array_t<bool> {
   detail::check_eigen_shape("lon", lon, "lat", lat);
   auto size = lon.size();
   auto result =
@@ -108,11 +107,10 @@ template <typename Geometry1, typename Geometry2>
 
 /// Calculate the distance between two geometries.
 template <typename Geometry1, typename Geometry2>
-[[nodiscard]] inline auto distance(const Geometry1 &geometry1,
-                                   const Geometry2 &geometry2,
-                                   const DistanceStrategy strategy,
-                                   const std::optional<Spheroid> &wgs)
-    -> double {
+[[nodiscard]] inline auto distance(
+    const Geometry1 &geometry1, const Geometry2 &geometry2,
+    const DistanceStrategy strategy,
+    const std::optional<Spheroid> &wgs) -> double {
   auto spheroid =
       wgs.has_value()
           ? static_cast<boost::geometry::srs::spheroid<double>>(*wgs)
diff --git a/src/pyinterp/core/include/pyinterp/geodetic/line_string.hpp b/src/pyinterp/core/include/pyinterp/geodetic/line_string.hpp
index 486fba6e..70ab3305 100644
--- a/src/pyinterp/core/include/pyinterp/geodetic/line_string.hpp
+++ b/src/pyinterp/core/include/pyinterp/geodetic/line_string.hpp
@@ -57,9 +57,8 @@ class LineString : public boost::geometry::model::linestring<Point> {
   ///
   /// @param rhs the linestring to test.
   /// @return true if the linestrings intersect.
-  [[nodiscard]] auto intersects(const LineString& rhs,
-                                const std::optional<Spheroid>& wgs) const
-      -> bool;
+  [[nodiscard]] auto intersects(
+      const LineString& rhs, const std::optional<Spheroid>& wgs) const -> bool;
 
   /// Compute the intersection of this linestring with another linestring.
   ///
@@ -134,8 +133,8 @@ class LineString : public boost::geometry::model::linestring<Point> {
 
   /// Returns the curvilinear distance along the linestring.
   [[nodiscard]] auto curvilinear_distance(
-      DistanceStrategy strategy, const std::optional<Spheroid>& wgs) const
-      -> Vector<double>;
+      DistanceStrategy strategy,
+      const std::optional<Spheroid>& wgs) const -> Vector<double>;
 
   /// Simplify the linestring using the Douglas-Peucker algorithm.
   ///
@@ -152,9 +151,8 @@ class LineString : public boost::geometry::model::linestring<Point> {
   /// @param point the point to project.
   /// @param wgs the spheroid to use.
   /// @return the closest point on this linestring to the given point.
-  [[nodiscard]] auto closest_point(const Point& point,
-                                   const std::optional<Spheroid>& wgs) const
-      -> Point;
+  [[nodiscard]] auto closest_point(
+      const Point& point, const std::optional<Spheroid>& wgs) const -> Point;
 
   /// Compute the closest point on this linestring to a given points.
   ///
diff --git a/src/pyinterp/core/include/pyinterp/geodetic/point.hpp b/src/pyinterp/core/include/pyinterp/geodetic/point.hpp
index c6f94cdf..3ff9be3d 100644
--- a/src/pyinterp/core/include/pyinterp/geodetic/point.hpp
+++ b/src/pyinterp/core/include/pyinterp/geodetic/point.hpp
@@ -66,17 +66,15 @@ class Point : public detail::geometry::GeographicPoint2D<double> {
   inline auto lat(double const v) -> void { this->set<1>(v); }
 
   /// Calculate the distance between the two points
-  [[nodiscard]] auto distance(const Point &other,
-                              const DistanceStrategy strategy,
-                              const std::optional<Spheroid> &wgs) const
-      -> double {
+  [[nodiscard]] auto distance(
+      const Point &other, const DistanceStrategy strategy,
+      const std::optional<Spheroid> &wgs) const -> double {
     return geodetic::distance(*this, other, strategy, wgs);
   }
 
   /// Returns the azimuth between the two points
-  [[nodiscard]] auto azimuth(const Point &other,
-                             const std::optional<Spheroid> &wgs) const
-      -> double;
+  [[nodiscard]] auto azimuth(
+      const Point &other, const std::optional<Spheroid> &wgs) const -> double;
 
   /// Get a tuple that fully encodes the state of this instance
   [[nodiscard]] auto getstate() const -> pybind11::tuple {
diff --git a/src/pyinterp/core/include/pyinterp/geohash/int64.hpp b/src/pyinterp/core/include/pyinterp/geohash/int64.hpp
index 1f266d51..dcdff094 100644
--- a/src/pyinterp/core/include/pyinterp/geohash/int64.hpp
+++ b/src/pyinterp/core/include/pyinterp/geohash/int64.hpp
@@ -36,8 +36,8 @@ namespace pyinterp::geohash::int64 {
 }
 
 // Encode a point into geohash with the given precision
-[[nodiscard]] auto encode(const geodetic::Point &point, uint32_t precision)
-    -> uint64_t;
+[[nodiscard]] auto encode(const geodetic::Point &point,
+                          uint32_t precision) -> uint64_t;
 
 // Encode points into geohash with the given precision
 [[nodiscard]] inline auto encode(const Eigen::Ref<const Eigen::VectorXd> &lon,
@@ -57,8 +57,8 @@ namespace pyinterp::geohash::int64 {
 
 // Returns the region encoded by the integer geohash with the specified
 // precision.
-[[nodiscard]] auto bounding_box(uint64_t hash, uint32_t precision)
-    -> geodetic::Box;
+[[nodiscard]] auto bounding_box(uint64_t hash,
+                                uint32_t precision) -> geodetic::Box;
 
 // Decode a hash into a geographic point with the given precision.
 // If round is true, the coordinates of the points will be rounded to the
@@ -107,7 +107,7 @@ namespace pyinterp::geohash::int64 {
 }
 
 // Returns all the GeoHash codes within the box.
-[[nodiscard]] auto bounding_boxes(const geodetic::Box &box, uint32_t precision)
-    -> Vector<uint64_t>;
+[[nodiscard]] auto bounding_boxes(const geodetic::Box &box,
+                                  uint32_t precision) -> Vector<uint64_t>;
 
 }  // namespace pyinterp::geohash::int64
diff --git a/src/pyinterp/core/include/pyinterp/geohash/string.hpp b/src/pyinterp/core/include/pyinterp/geohash/string.hpp
index 58855edc..f16fc6fd 100644
--- a/src/pyinterp/core/include/pyinterp/geohash/string.hpp
+++ b/src/pyinterp/core/include/pyinterp/geohash/string.hpp
@@ -51,8 +51,8 @@ class Array {
                            capsule_);
   }
 
-  static auto get_info(const pybind11::array &hash, pybind11::ssize_t ndim)
-      -> pybind11::buffer_info;
+  static auto get_info(const pybind11::array &hash,
+                       pybind11::ssize_t ndim) -> pybind11::buffer_info;
 
  private:
   std::vector<char> *array_;
@@ -65,8 +65,8 @@ class Array {
 auto allocate_array(size_t size, uint32_t precision) -> Array;
 
 /// Encode a point into geohash with the given bit depth
-auto encode(const geodetic::Point &point, char *buffer, uint32_t precision)
-    -> void;
+auto encode(const geodetic::Point &point, char *buffer,
+            uint32_t precision) -> void;
 
 /// Encode points into geohash with the given bit depth
 [[nodiscard]] auto encode(const Eigen::Ref<const Eigen::VectorXd> &lon,
@@ -74,14 +74,14 @@ auto encode(const geodetic::Point &point, char *buffer, uint32_t precision)
                           uint32_t precision) -> pybind11::array;
 
 /// Returns the region encoded
-[[nodiscard]] auto bounding_box(const char *hash, size_t count)
-    -> geodetic::Box;
+[[nodiscard]] auto bounding_box(const char *hash,
+                                size_t count) -> geodetic::Box;
 
 /// Decode a hash into a spherical equatorial point. If round is true, the
 /// coordinates of the points will be rounded to the accuracy defined by the
 /// GeoHash.
-[[nodiscard]] auto decode(const char *hash, size_t count, bool round)
-    -> geodetic::Point;
+[[nodiscard]] auto decode(const char *hash, size_t count,
+                          bool round) -> geodetic::Point;
 
 /// Decode hashes into a spherical equatorial points. If round is true, the
 /// coordinates of the points will be rounded to the accuracy defined by the
@@ -113,22 +113,23 @@ auto encode(const geodetic::Point &point, char *buffer, uint32_t precision)
                                   uint32_t precision) -> pybind11::array;
 
 [[nodiscard]] auto bounding_boxes(const geodetic::Polygon &polygon,
-                                  uint32_t precision, size_t num_threads)
-    -> pybind11::array;
+                                  uint32_t precision,
+                                  size_t num_threads) -> pybind11::array;
 
 [[nodiscard]] auto bounding_boxes(const geodetic::MultiPolygon &polygons,
-                                  uint32_t precision, size_t num_threads)
-    -> pybind11::array;
+                                  uint32_t precision,
+                                  size_t num_threads) -> pybind11::array;
 
 /// Returns the start and end indexes of the different GeoHash boxes.
-[[nodiscard]] auto where(const pybind11::array &hash) -> std::unordered_map<
-    std::string,
-    std::tuple<std::tuple<int64_t, int64_t>, std::tuple<int64_t, int64_t>>>;
+[[nodiscard]] auto where(const pybind11::array &hash)
+    -> std::unordered_map<
+        std::string,
+        std::tuple<std::tuple<int64_t, int64_t>, std::tuple<int64_t, int64_t>>>;
 
 /// Transforms the given codes from one precision to another. If the given
 /// precision is higher than the precision of the given codes, the result
 /// contains a zoom in, otherwise it contains a zoom out.
-[[nodiscard]] auto transform(const pybind11::array &hash, uint32_t precision)
-    -> pybind11::array;
+[[nodiscard]] auto transform(const pybind11::array &hash,
+                             uint32_t precision) -> pybind11::array;
 
 }  // namespace pyinterp::geohash::string
diff --git a/src/pyinterp/core/include/pyinterp/period.hpp b/src/pyinterp/core/include/pyinterp/period.hpp
index 65e388f0..0398d4c0 100644
--- a/src/pyinterp/core/include/pyinterp/period.hpp
+++ b/src/pyinterp/core/include/pyinterp/period.hpp
@@ -154,8 +154,8 @@ struct Period {
   /// @param tolerance The maximum distance from the period to consider the date
   /// close.
   /// @returns true if the date is close to the period.
-  constexpr auto is_close(const int64_t date, const int64_t tolerance) const
-      -> bool {
+  constexpr auto is_close(const int64_t date,
+                          const int64_t tolerance) const -> bool {
     return (date >= begin - tolerance) && (date <= last + tolerance);
   }
 
@@ -245,9 +245,8 @@ class PeriodList {
   /// True if the date belongs to the period.
   ///
   /// The vectors provided should be sorted in chronological order.
-  auto cross_a_period(
-      const Eigen::Ref<const Eigen::Matrix<int64_t, -1, 1>> &dates) const
-      -> Eigen::Matrix<bool, -1, 1> {
+  auto cross_a_period(const Eigen::Ref<const Eigen::Matrix<int64_t, -1, 1>>
+                          &dates) const -> Eigen::Matrix<bool, -1, 1> {
     // The index of the closest period of the current date processed.
     auto first_index = int64_t(0);
 
@@ -309,9 +308,8 @@ class PeriodList {
   /// if the date belongs to the period.
   ///
   /// The vectors provided should be sorted in chronological order.
-  auto belong_to_a_period(
-      const Eigen::Ref<const Eigen::Matrix<int64_t, -1, 1>> &dates) const
-      -> Eigen::Matrix<bool, -1, 1> {
+  auto belong_to_a_period(const Eigen::Ref<const Eigen::Matrix<int64_t, -1, 1>>
+                              &dates) const -> Eigen::Matrix<bool, -1, 1> {
     // Flag equal to 1 if the date belongs to a period, 0 otherwise.
     auto flags = Eigen::Matrix<bool, -1, 1>(dates.size());
 
diff --git a/src/pyinterp/core/include/pyinterp/quadrivariate.hpp b/src/pyinterp/core/include/pyinterp/quadrivariate.hpp
index fa3b6940..f72b9a00 100644
--- a/src/pyinterp/core/include/pyinterp/quadrivariate.hpp
+++ b/src/pyinterp/core/include/pyinterp/quadrivariate.hpp
@@ -123,8 +123,8 @@ auto quadrivariate(const Grid4D<Type, AxisType> &grid,
                    const Bivariate4D<Point, Coordinate> *interpolator,
                    const std::optional<std::string> &z_method,
                    const std::optional<std::string> &u_method,
-                   const bool bounds_error, const size_t num_threads)
-    -> pybind11::array_t<Coordinate> {
+                   const bool bounds_error,
+                   const size_t num_threads) -> pybind11::array_t<Coordinate> {
   pyinterp::detail::check_array_ndim("x", 1, x, "y", 1, y, "z", 1, z, "u", 1,
                                      u);
   pyinterp::detail::check_ndarray_shape("x", x, "y", y, "z", z, "u", u);
diff --git a/src/pyinterp/core/include/pyinterp/rtree.hpp b/src/pyinterp/core/include/pyinterp/rtree.hpp
index 0d9d649f..b43c2e28 100644
--- a/src/pyinterp/core/include/pyinterp/rtree.hpp
+++ b/src/pyinterp/core/include/pyinterp/rtree.hpp
@@ -189,8 +189,8 @@ class RTree : public detail::geometry::RTree<Point, Type> {
   /// Search for the nearest K nearest neighbors of a given coordinates.
   auto value(const array_t &coordinates,
              const std::optional<coordinate_t> &radius, const uint32_t k,
-             const bool within, const size_t num_threads) const
-      -> pybind11::tuple {
+             const bool within,
+             const size_t num_threads) const -> pybind11::tuple {
     detail::check_array_ndim("coordinates", 2, coordinates);
     switch (coordinates.shape(1)) {
       case dimension_t::value - 1:
@@ -207,12 +207,10 @@ class RTree : public detail::geometry::RTree<Point, Type> {
   }
 
   /// TODO
-  auto inverse_distance_weighting(const array_t &coordinates,
-                                  const std::optional<coordinate_t> &radius,
-                                  const uint32_t k, const uint32_t p,
-                                  const bool within,
-                                  const size_t num_threads) const
-      -> pybind11::tuple {
+  auto inverse_distance_weighting(
+      const array_t &coordinates, const std::optional<coordinate_t> &radius,
+      const uint32_t k, const uint32_t p, const bool within,
+      const size_t num_threads) const -> pybind11::tuple {
     detail::check_array_ndim("coordinates", 2, coordinates);
 
     switch (coordinates.shape(1)) {
@@ -232,13 +230,11 @@ class RTree : public detail::geometry::RTree<Point, Type> {
   }
 
   /// TODO
-  auto radial_basis_function(const array_t &coordinates,
-                             const std::optional<coordinate_t> &radius,
-                             const uint32_t k, const RadialBasisFunction rbf,
-                             const std::optional<promotion_t> &epsilon,
-                             const promotion_t smooth, const bool within,
-                             const size_t num_threads) const
-      -> pybind11::tuple {
+  auto radial_basis_function(
+      const array_t &coordinates, const std::optional<coordinate_t> &radius,
+      const uint32_t k, const RadialBasisFunction rbf,
+      const std::optional<promotion_t> &epsilon, const promotion_t smooth,
+      const bool within, const size_t num_threads) const -> pybind11::tuple {
     detail::check_array_ndim("coordinates", 2, coordinates);
     switch (coordinates.shape(1)) {
       case dimension_t::value - 1:
@@ -263,8 +259,8 @@ class RTree : public detail::geometry::RTree<Point, Type> {
                        const std::optional<coordinate_t> &radius,
                        const uint32_t k, const WindowFunction wf,
                        const std::optional<coordinate_t> &arg,
-                       const bool within, const size_t num_threads) const
-      -> pybind11::tuple {
+                       const bool within,
+                       const size_t num_threads) const -> pybind11::tuple {
     detail::check_array_ndim("coordinates", 2, coordinates);
     switch (coordinates.shape(1)) {
       case dimension_t::value - 1:
@@ -286,8 +282,8 @@ class RTree : public detail::geometry::RTree<Point, Type> {
                          const std::optional<coordinate_t> &radius,
                          const uint32_t k, const CovarianceFunction cov,
                          const coordinate_t sigma, const coordinate_t lambda,
-                         const bool within, const size_t num_threads) const
-      -> pybind11::tuple {
+                         const bool within,
+                         const size_t num_threads) const -> pybind11::tuple {
     detail::check_array_ndim("coordinates", 2, coordinates);
 
     switch (coordinates.shape(1)) {
@@ -369,9 +365,8 @@ class RTree : public detail::geometry::RTree<Point, Type> {
   /// Create the cartesian point for the given coordinates: longitude and
   /// latitude in degrees, altitude in meters, then the other coordinates
   /// defined in a Euclidean space.
-  auto from_lon_lat_alt(
-      const Eigen::Map<const Vector<coordinate_t>> &coordinates) const
-      -> Point {
+  auto from_lon_lat_alt(const Eigen::Map<const Vector<coordinate_t>>
+                            &coordinates) const -> Point {
     auto ecef = coordinates_.lla_to_ecef(
         detail::geometry::EquatorialPoint3D<coordinate_t>{
             coordinates(0), coordinates(1), coordinates(2)});
@@ -553,8 +548,8 @@ class RTree : public detail::geometry::RTree<Point, Type> {
   auto _value(Converter converter,
               const pybind11::array_t<coordinate_t> &coordinates,
               const std::optional<coordinate_t> &radius, const uint32_t k,
-              const bool within, const size_t num_threads) const
-      -> pybind11::tuple {
+              const bool within,
+              const size_t num_threads) const -> pybind11::tuple {
     auto _coordinates = coordinates.template unchecked<2>();
     auto size = coordinates.shape(0);
 
@@ -685,8 +680,8 @@ class RTree : public detail::geometry::RTree<Point, Type> {
                           const coordinate_t radius, const uint32_t k,
                           const CovarianceFunction &cov,
                           const coordinate_t sigma, const coordinate_t lambda,
-                          const bool within, const size_t num_threads) const
-      -> pybind11::tuple {
+                          const bool within,
+                          const size_t num_threads) const -> pybind11::tuple {
     auto _coordinates = coordinates.template unchecked<2>();
     auto size = coordinates.shape(0);
 
@@ -802,8 +797,8 @@ class RTree : public detail::geometry::RTree<Point, Type> {
                         const pybind11::array_t<coordinate_t> &coordinates,
                         const coordinate_t radius, const uint32_t k,
                         const WindowFunction wf, const coordinate_t arg,
-                        const bool within, const size_t num_threads) const
-      -> pybind11::tuple {
+                        const bool within,
+                        const size_t num_threads) const -> pybind11::tuple {
     auto _coordinates = coordinates.template unchecked<2>();
     auto size = coordinates.shape(0);
 
diff --git a/src/pyinterp/core/include/pyinterp/streaming_histogram.hpp b/src/pyinterp/core/include/pyinterp/streaming_histogram.hpp
index 3da2f47f..ebf00d58 100644
--- a/src/pyinterp/core/include/pyinterp/streaming_histogram.hpp
+++ b/src/pyinterp/core/include/pyinterp/streaming_histogram.hpp
@@ -209,8 +209,8 @@ class StreamingHistogram {
 
   /// Allocates the accumulators needed for the result statistics.
   [[nodiscard]] static auto allocate_accumulators(
-      const size_t bin_count, const std::vector<pybind11::ssize_t> &shape)
-      -> Vector<Accumulators> {
+      const size_t bin_count,
+      const std::vector<pybind11::ssize_t> &shape) -> Vector<Accumulators> {
     auto result = Vector<Accumulators>(StreamingHistogram::shape_size(shape));
     std::for_each(result.data(), result.data() + result.size(),
                   [bin_count](auto &acc) { acc.resize(bin_count); });
diff --git a/src/pyinterp/core/include/pyinterp/temporal_axis.hpp b/src/pyinterp/core/include/pyinterp/temporal_axis.hpp
index bfaa87ca..ffd8fed2 100644
--- a/src/pyinterp/core/include/pyinterp/temporal_axis.hpp
+++ b/src/pyinterp/core/include/pyinterp/temporal_axis.hpp
@@ -21,8 +21,8 @@ namespace detail {
 
 /// Check the properties of the numpy array and returns the type of dates
 /// encoded in this parameter.
-inline auto dtype(const std::string &name, const pybind11::array &array)
-    -> dateutils::DType {
+inline auto dtype(const std::string &name,
+                  const pybind11::array &array) -> dateutils::DType {
   detail::check_array_ndim(name, 1, array);
   if ((array.flags() & pybind11::array::c_style) == 0) {
     throw pybind11::type_error(name + " must be a C-style contiguous array");
@@ -129,8 +129,8 @@ class TemporalAxis : public Axis<int64_t> {
   }
 
   /// @copydoc Axis::find_index(const pybind11::array_t<int64_t>, bool) const
-  auto find_index(const pybind11::array &coordinates, bool bounded) const
-      -> pybind11::array_t<int64_t> {
+  auto find_index(const pybind11::array &coordinates,
+                  bool bounded) const -> pybind11::array_t<int64_t> {
     const auto values = safe_cast("coordinates", coordinates);
     const auto size = values.size();
     auto result = pybind11::array_t<int64_t>(size);
@@ -248,9 +248,8 @@ class TemporalAxis : public Axis<int64_t> {
 
   /// Returns the result of a method wrapped in a pybind11::array.
   template <typename Getter>
-  inline auto as_datetype(const Getter &getter,
-                          const dateutils::DType &dtype) const
-      -> pybind11::array {
+  inline auto as_datetype(const Getter &getter, const dateutils::DType &dtype)
+      const -> pybind11::array {
     auto value = getter();
     return {pybind11::dtype(static_cast<std::string>(dtype)), {}, {}, &value};
   }
diff --git a/src/pyinterp/core/include/pyinterp/trivariate.hpp b/src/pyinterp/core/include/pyinterp/trivariate.hpp
index 8d53c30e..cbc04be1 100644
--- a/src/pyinterp/core/include/pyinterp/trivariate.hpp
+++ b/src/pyinterp/core/include/pyinterp/trivariate.hpp
@@ -85,8 +85,8 @@ auto trivariate(const Grid3D<Type, AxisType> &grid,
                 const pybind11::array_t<AxisType> &z,
                 const Bivariate3D<Point, Coordinate> *interpolator,
                 const std::optional<std::string> &z_method,
-                const bool bounds_error, const size_t num_threads)
-    -> pybind11::array_t<Coordinate> {
+                const bool bounds_error,
+                const size_t num_threads) -> pybind11::array_t<Coordinate> {
   pyinterp::detail::check_array_ndim("x", 1, x, "y", 1, y, "z", 1, z);
   pyinterp::detail::check_ndarray_shape("x", x, "y", y, "z", z);
   auto z_interpolation_method =
diff --git a/src/pyinterp/core/lib/geodetic/line_string.cpp b/src/pyinterp/core/lib/geodetic/line_string.cpp
index 2340cc50..031251e8 100644
--- a/src/pyinterp/core/lib/geodetic/line_string.cpp
+++ b/src/pyinterp/core/lib/geodetic/line_string.cpp
@@ -154,10 +154,9 @@ auto LineString::curvilinear_distance(DistanceStrategy strategy,
   return result;
 }
 
-auto LineString::simplify(const double tolerance,
-                          const DistanceStrategy strategy,
-                          const std::optional<Spheroid>& wgs) const
-    -> LineString {
+auto LineString::simplify(
+    const double tolerance, const DistanceStrategy strategy,
+    const std::optional<Spheroid>& wgs) const -> LineString {
   auto spheroid =
       wgs.has_value()
           ? static_cast<boost::geometry::srs::spheroid<double>>(*wgs)
@@ -210,9 +209,8 @@ inline auto closest_point(const LineString& ls, const Point& point,
 
 }  // namespace impl
 
-auto LineString::closest_point(const Point& point,
-                               const std::optional<Spheroid>& wgs) const
-    -> Point {
+auto LineString::closest_point(
+    const Point& point, const std::optional<Spheroid>& wgs) const -> Point {
   return wgs.has_value() ? impl::closest_point(*this, point, *wgs)
                          : impl::closest_point(*this, point);
 }
diff --git a/src/pyinterp/core/lib/geodetic/rtree.cpp b/src/pyinterp/core/lib/geodetic/rtree.cpp
index 0a6aedbb..3947d996 100644
--- a/src/pyinterp/core/lib/geodetic/rtree.cpp
+++ b/src/pyinterp/core/lib/geodetic/rtree.cpp
@@ -85,8 +85,8 @@ auto RTree::insert(const Eigen::Ref<const Vector<double>> &lon,
 
 auto RTree::query(const Eigen::Ref<const Vector<double>> &lon,
                   const Eigen::Ref<const Vector<double>> &lat, const uint32_t k,
-                  const bool within, const size_t num_threads) const
-    -> pybind11::tuple {
+                  const bool within,
+                  const size_t num_threads) const -> pybind11::tuple {
   detail::check_container_size("lon", lon, "lat", lat);
   auto size = lon.size();
 
@@ -193,15 +193,13 @@ auto RTree::inverse_distance_weighting(
   return pybind11::make_tuple(data, neighbors);
 }
 
-auto RTree::radial_basis_function(const Eigen::Ref<const Vector<double>> &lon,
-                                  const Eigen::Ref<const Vector<double>> &lat,
-                                  const std::optional<double> &radius,
-                                  const uint32_t k,
-                                  const detail::math::RadialBasisFunction rbf,
-                                  const std::optional<double> &epsilon,
-                                  const double smooth, const bool within,
-                                  const size_t num_threads) const
-    -> pybind11::tuple {
+auto RTree::radial_basis_function(
+    const Eigen::Ref<const Vector<double>> &lon,
+    const Eigen::Ref<const Vector<double>> &lat,
+    const std::optional<double> &radius, const uint32_t k,
+    const detail::math::RadialBasisFunction rbf,
+    const std::optional<double> &epsilon, const double smooth,
+    const bool within, const size_t num_threads) const -> pybind11::tuple {
   detail::check_container_size("lon", lon, "lat", lat);
   auto _radius = radius.value_or(std::numeric_limits<double>::max());
   auto size = lon.size();
diff --git a/src/pyinterp/core/lib/geohash/int64.cpp b/src/pyinterp/core/lib/geohash/int64.cpp
index 749c8a9f..b6b1ef6c 100644
--- a/src/pyinterp/core/lib/geohash/int64.cpp
+++ b/src/pyinterp/core/lib/geohash/int64.cpp
@@ -243,8 +243,8 @@ auto allocate_array(const size_t size) -> Vector<uint64_t> {
 }
 
 // ---------------------------------------------------------------------------
-auto encode(const geodetic::Point &point, const uint32_t precision)
-    -> uint64_t {
+auto encode(const geodetic::Point &point,
+            const uint32_t precision) -> uint64_t {
   auto result = encoder(point.lat(), point.lon());
   if (precision != 64) {
     result >>= (64 - precision);
@@ -253,8 +253,8 @@ auto encode(const geodetic::Point &point, const uint32_t precision)
 }
 
 // ---------------------------------------------------------------------------
-auto bounding_box(const uint64_t hash, const uint32_t precision)
-    -> geodetic::Box {
+auto bounding_box(const uint64_t hash,
+                  const uint32_t precision) -> geodetic::Box {
   auto full_hash = hash << (64U - precision);
   auto lat_lng_int = deinterleaver(full_hash);
   auto lat = codec::decode_range(std::get<0>(lat_lng_int), 90);
@@ -268,8 +268,8 @@ auto bounding_box(const uint64_t hash, const uint32_t precision)
 }
 
 // ---------------------------------------------------------------------------
-auto neighbors(const uint64_t hash, const uint32_t precision)
-    -> Eigen::Matrix<uint64_t, 8, 1> {
+auto neighbors(const uint64_t hash,
+               const uint32_t precision) -> Eigen::Matrix<uint64_t, 8, 1> {
   auto box = bounding_box(hash, precision);
   auto center = box.centroid();
   auto [lon_delta, lat_delta] = box.delta(false);
@@ -325,8 +325,8 @@ auto grid_properties(const geodetic::Box &box, const uint32_t precision)
 }
 
 // ---------------------------------------------------------------------------
-auto bounding_boxes(const geodetic::Box &box, const uint32_t precision)
-    -> Vector<uint64_t> {
+auto bounding_boxes(const geodetic::Box &box,
+                    const uint32_t precision) -> Vector<uint64_t> {
   // Grid resolution in degrees
   const auto [lng_err, lat_err] = error_with_precision(precision);
 
diff --git a/src/pyinterp/core/lib/geohash/string.cpp b/src/pyinterp/core/lib/geohash/string.cpp
index eb86a5b0..830edb13 100644
--- a/src/pyinterp/core/lib/geohash/string.cpp
+++ b/src/pyinterp/core/lib/geohash/string.cpp
@@ -17,8 +17,8 @@ namespace pyinterp::geohash::string {
 static const auto base32 = Base32();
 
 // ---------------------------------------------------------------------------
-auto Array::get_info(const pybind11::array &hash, const pybind11::ssize_t ndim)
-    -> pybind11::buffer_info {
+auto Array::get_info(const pybind11::array &hash,
+                     const pybind11::ssize_t ndim) -> pybind11::buffer_info {
   auto info = hash.request();
   auto dtype = hash.dtype();
   if ((hash.flags() & pybind11::array::c_style) == 0) {
@@ -114,15 +114,15 @@ auto bounding_box(const char *const hash, const size_t count) -> geodetic::Box {
 }
 
 // ---------------------------------------------------------------------------
-auto decode(const char *const hash, const size_t count, const bool round)
-    -> geodetic::Point {
+auto decode(const char *const hash, const size_t count,
+            const bool round) -> geodetic::Point {
   auto bbox = bounding_box(hash, count);
   return round ? bbox.round() : bbox.centroid();
 }
 
 // ---------------------------------------------------------------------------
-auto decode(const pybind11::array &hash, const bool round)
-    -> std::tuple<Eigen::VectorXd, Eigen::VectorXd> {
+auto decode(const pybind11::array &hash,
+            const bool round) -> std::tuple<Eigen::VectorXd, Eigen::VectorXd> {
   auto info = Array::get_info(hash, 1);
   auto count = info.strides[0];
   auto lon = Eigen::VectorXd(info.shape[0]);
@@ -264,8 +264,8 @@ auto mask_cell(const geodetic::Box &box, const Geometry &geometry,
 static auto select_cell(const double lng_err, const double lat_err,
                         const geodetic::Point &point_sw, const size_t lon_step,
                         const size_t lat_step, const uint32_t bits,
-                        const uint32_t precision, const Matrix<bool> &mask)
-    -> pybind11::array {
+                        const uint32_t precision,
+                        const Matrix<bool> &mask) -> pybind11::array {
   // Count the number of cells that are enclosed by the polygon
   auto size = std::count(mask.data(), mask.data() + mask.size(), true);
 
@@ -330,16 +330,17 @@ auto bounding_boxes(const geodetic::Polygon &polygon, const uint32_t precision,
 
 // ---------------------------------------------------------------------------
 auto bounding_boxes(const geodetic::MultiPolygon &polygons,
-                    const uint32_t precision, const size_t num_threads)
-    -> pybind11::array {
+                    const uint32_t precision,
+                    const size_t num_threads) -> pybind11::array {
   return bounding_boxes<geodetic::MultiPolygon>(polygons, precision,
                                                 num_threads);
 }
 
 // ---------------------------------------------------------------------------
-auto where(const pybind11::array &hash) -> std::unordered_map<
-    std::string,
-    std::tuple<std::tuple<int64_t, int64_t>, std::tuple<int64_t, int64_t>>> {
+auto where(const pybind11::array &hash)
+    -> std::unordered_map<std::string,
+                          std::tuple<std::tuple<int64_t, int64_t>,
+                                     std::tuple<int64_t, int64_t>>> {
   // Index shifts of neighboring pixels
   static const auto shift_row =
       std::array<int64_t, 8>{-1, -1, -1, 0, 1, 0, 1, 1};
@@ -456,8 +457,8 @@ static auto zoom_out(char *ptr, pybind11::ssize_t size, uint32_t from_precision,
 }
 
 // ---------------------------------------------------------------------------
-auto transform(const pybind11::array &hash, uint32_t precision)
-    -> pybind11::array {
+auto transform(const pybind11::array &hash,
+               uint32_t precision) -> pybind11::array {
   // Decode the information in the provided table.
   auto info = Array::get_info(hash, 1);
   auto size = info.shape[0];
diff --git a/src/pyinterp/core/module/bicubic.cpp b/src/pyinterp/core/module/bicubic.cpp
index a0974426..6ce1925e 100644
--- a/src/pyinterp/core/module/bicubic.cpp
+++ b/src/pyinterp/core/module/bicubic.cpp
@@ -41,8 +41,8 @@ template <typename DataType, typename Interpolator>
 auto bicubic(const Grid2D<DataType> &grid, const py::array_t<double> &x,
              const py::array_t<double> &y, Eigen::Index nx, Eigen::Index ny,
              const std::string &fitting_model, const std::string &boundary,
-             const bool bounds_error, size_t num_threads)
-    -> py::array_t<double> {
+             const bool bounds_error,
+             size_t num_threads) -> py::array_t<double> {
   detail::check_array_ndim("x", 1, x, "y", 1, y);
   detail::check_ndarray_shape("x", x, "y", y);
 
@@ -166,8 +166,8 @@ auto bicubic_4d(const Grid4D<DataType, AxisType> &grid,
                 const py::array_t<AxisType> &z, const py::array_t<double> &u,
                 Eigen::Index nx, Eigen::Index ny,
                 const std::string &fitting_model, const std::string &boundary,
-                const bool bounds_error, size_t num_threads)
-    -> py::array_t<double> {
+                const bool bounds_error,
+                size_t num_threads) -> py::array_t<double> {
   detail::check_array_ndim("x", 1, x, "y", 1, y, "z", 1, z, "u", 1, u);
   detail::check_ndarray_shape("x", x, "y", y, "z", z, "u", u);
   auto boundary_type = parse_axis_boundary(boundary);
diff --git a/src/pyinterp/orbit.py b/src/pyinterp/orbit.py
index 9ddb9bd8..d76c85bc 100644
--- a/src/pyinterp/orbit.py
+++ b/src/pyinterp/orbit.py
@@ -8,7 +8,7 @@
 """
 from __future__ import annotations
 
-from typing import Iterator
+from collections.abc import Iterator
 import dataclasses
 
 import numpy
@@ -284,7 +284,7 @@ def iterate(
 
             # Update of the number of the next pass to be generated
             absolute_pass_number += 1
-        return StopIteration
+        return StopIteration  # type: ignore[return-value]
 
 
 @dataclasses.dataclass(frozen=True)
diff --git a/src/pyinterp/statistics/descriptive_descriptive.py b/src/pyinterp/statistics/descriptive_descriptive.py
index 9eef08f1..0b234c3c 100644
--- a/src/pyinterp/statistics/descriptive_descriptive.py
+++ b/src/pyinterp/statistics/descriptive_descriptive.py
@@ -8,7 +8,8 @@
 """
 from __future__ import annotations
 
-from typing import Any, Iterable
+from typing import Any
+from collections.abc import Iterable
 import copy
 
 import dask.array.core
diff --git a/src/pyinterp/statistics/streaming_histogram.py b/src/pyinterp/statistics/streaming_histogram.py
index e5040b05..087e16d7 100644
--- a/src/pyinterp/statistics/streaming_histogram.py
+++ b/src/pyinterp/statistics/streaming_histogram.py
@@ -8,7 +8,8 @@
 """
 from __future__ import annotations
 
-from typing import Any, Iterable
+from typing import Any
+from collections.abc import Iterable
 
 import dask.array.core
 import numpy
diff --git a/src/pyinterp/tests/__init__.py b/src/pyinterp/tests/__init__.py
index 5edb9144..041248dc 100644
--- a/src/pyinterp/tests/__init__.py
+++ b/src/pyinterp/tests/__init__.py
@@ -162,7 +162,7 @@ def parallel_lines():
     return ROOT.joinpath('parallel_lines.json')
 
 
-def run(pattern: Optional[str] = None) -> None:
+def run(pattern: str | None = None) -> None:
     """Run tests.
 
     Args:
diff --git a/src/pyinterp/tests/core/test_dateutils.py b/src/pyinterp/tests/core/test_dateutils.py
index c8f1ee33..cd34292d 100644
--- a/src/pyinterp/tests/core/test_dateutils.py
+++ b/src/pyinterp/tests/core/test_dateutils.py
@@ -13,7 +13,7 @@
 
 
 def make_date(samples=10000,
-              resolution='us') -> Tuple[List[datetime.datetime], np.ndarray]:
+              resolution='us') -> tuple[list[datetime.datetime], np.ndarray]:
     epoch = datetime.datetime(1970, 1, 1)
     delta = datetime.datetime.now() - datetime.datetime(1970, 1, 1)
 
diff --git a/src/pyinterp/tests/core/test_geodetic.py b/src/pyinterp/tests/core/test_geodetic.py
index 95a4ceac..8f16aab1 100644
--- a/src/pyinterp/tests/core/test_geodetic.py
+++ b/src/pyinterp/tests/core/test_geodetic.py
@@ -702,7 +702,7 @@ def test_linestring():
     assert a == b
 
 
-def get_coordinates() -> Tuple[np.ndarray, np.ndarray]:
+def get_coordinates() -> tuple[np.ndarray, np.ndarray]:
     """Get coordinates for a line string."""
     lon = np.array([
         -9.72270435, -9.58849868, -9.45429341, -9.32009001, -9.18588996,
diff --git a/src/pyinterp/tests/test_orbit_interpolation.py b/src/pyinterp/tests/test_orbit_interpolation.py
index f8cd55f0..abbe4f83 100644
--- a/src/pyinterp/tests/test_orbit_interpolation.py
+++ b/src/pyinterp/tests/test_orbit_interpolation.py
@@ -14,7 +14,7 @@
 
 def load_test_ephemeris(
     filename: os.PathLike
-) -> Tuple[float, NDArray, NDArray, NDArray, np.timedelta64]:
+) -> tuple[float, NDArray, NDArray, NDArray, np.timedelta64]:
     """Loads the ephemeris from a text file.
 
     Args:
@@ -27,7 +27,7 @@ def load_test_ephemeris(
     with open(filename) as stream:
         lines = stream.readlines()
 
-    def to_dict(comments) -> Dict[str, float]:
+    def to_dict(comments) -> dict[str, float]:
         """Returns a dictionary describing the parameters of the orbit."""
         result = {}
         for item in comments: