fix: validate indices and sizes instead of relying on asserts

README.md

@@ -124,6 +124,12 @@ f.d("y")   # df/dy
 >>> 1.0
 ```
 
+## Validation
+
+Arguments coming from Python are validated. Out-of-range Hessian indices raise `IndexError`; inconsistent sizes raise `RuntimeError` — `eval` with the wrong number of values, `set_hm` with a mismatched shape, a negative size, a variable index outside the gradient, or padding below the current size.
+
+In C++ the element accessors `g(i)`, `h(i)` and `h(i, j)` treat their indices as a precondition: like `std::vector::operator[]` they are only checked via `assert` in debug builds, so they stay free of branches in hot loops. Everything that creates, resizes or evaluates a scalar validates its arguments and throws `std::runtime_error`. Define `HYPERJET_NO_EXCEPTIONS` to fall back to `assert` throughout.
+
 ## NumPy Integration
 
 HyperJet scalars work with NumPy for vector and matrix operations.

include/hyperjet/hyperjet.h

@@ -95,6 +95,28 @@ class DDScalar {
     }
   }
 
+  HYPERJET_INLINE static void check_valid_index(const index i,
+                                                const index size) {
+    if constexpr (throw_exceptions()) {
+      if (i < 0 || size <= i) {
+        throw std::runtime_error("Invalid index");
+      }
+    } else {
+      assert(0 <= i && i < size && "Invalid index");
+    }
+  }
+
+  HYPERJET_INLINE static void check_pad_size(const index size,
+                                             const index new_size) {
+    if constexpr (throw_exceptions()) {
+      if (new_size < size) {
+        throw std::runtime_error("Invalid new size");
+      }
+    } else {
+      assert(new_size >= size && "Invalid new size");
+    }
+  }
+
   HYPERJET_INLINE static void check_equal_size(const index size_a,
                                                const index size_b) {
     if constexpr (TSize == Dynamic) {
@@ -340,6 +362,7 @@ class DDScalar {
 
   void resize(const index size) {
     static_assert(is_dynamic());
+    check_valid_size(size);
     m_size = size;
     const index n = data_length_from_size(size);
     m_data.resize(n);
@@ -348,6 +371,8 @@ class DDScalar {
   Type pad_left(const index new_size) const {
     static_assert(is_dynamic());
 
+    check_pad_size(size(), new_size);
+
     Type result = empty(new_size);
 
     const index head = new_size - size();
@@ -387,6 +412,8 @@ class DDScalar {
   Type pad_right(const index new_size) const {
     static_assert(is_dynamic());
 
+    check_pad_size(size(), new_size);
+
     Type result = empty(new_size);
 
     const index tail = new_size - size();
@@ -454,6 +481,7 @@ class DDScalar {
 
   static Type empty(const index size) {
     if constexpr (is_dynamic()) {
+      check_valid_size(size);
       const index n = data_length_from_size(size);
       const Data data(n);
       Type result(data, size);
@@ -479,6 +507,7 @@ class DDScalar {
 
   static Type zero(const index size) {
     if constexpr (is_dynamic()) {
+      check_valid_size(size);
       const Data data(data_length_from_size(size), 0);
       Type result(data, size);
       return result;
@@ -507,6 +536,7 @@ class DDScalar {
 
   static Type variable(const index i, const Scalar f) {
     static_assert(!is_dynamic());
+    check_valid_index(i, TSize);
     Type result = zero();
     result.f() = f;
     result.g(i) = 1;
@@ -516,6 +546,7 @@ class DDScalar {
   static Type variable(const index i, const Scalar f, const index size) {
     if constexpr (is_dynamic()) {
       Type result = zero(size);
+      check_valid_index(i, size);
       result.f() = f;
       result.g(i) = 1;
       return result;
@@ -635,6 +666,9 @@ class DDScalar {
   void hm(const std::string mode, Eigen::Ref<Matrix> out) const
     requires(order() == 2)
   {
+    check_equal_size(size(), out.rows());
+    check_equal_size(size(), out.cols());
+
     index it = 0;
 
     for (index i = 0; i < size(); i++) {
@@ -663,6 +697,9 @@ class DDScalar {
   void set_hm(const Eigen::Ref<const Matrix> &value)
     requires(order() == 2)
   {
+    check_equal_size(size(), value.rows());
+    check_equal_size(size(), value.cols());
+
     index it = 0;
 
     for (index i = 0; i < size(); i++) {
@@ -709,6 +746,8 @@ class DDScalar {
 
   Scalar eval(typename std::conditional < TSize == Dynamic, std::vector<Scalar>,
               std::array<Scalar, TSize<0 ? 0 : TSize>>::type d) const {
+    check_equal_size(size(), length(d));
+
     Scalar result = f();
 
     for (index i = 0; i < size(); i++) {

python/src/common.h

@@ -17,6 +17,18 @@ namespace hj = hyperjet;
 namespace py = pybind11;
 using namespace py::literals;
 
+// The accessors in the C++ header treat their indices as a precondition and
+// only check them via assert, which is compiled out in release builds. Python
+// callers are untrusted, so the indices have to be validated here.
+template <typename T>
+void check_index(const T &self, const hj::index row, const hj::index col) {
+  if (row < 0 || self.size() <= row || col < 0 || self.size() <= col) {
+    throw py::index_error("index (" + std::to_string(row) + ", " +
+                          std::to_string(col) + ") is out of range for size " +
+                          std::to_string(self.size()));
+  }
+}
+
 template <typename T> auto bind(py::module &m, const std::string &name) {
   py::class_<T> cls(m, name.c_str());
 
@@ -99,12 +111,21 @@ template <typename T> auto bind(py::module &m, const std::string &name) {
     cls.def(
            "h",
            [](const T &self, hj::index row, hj::index col) ->
-           typename T::Scalar { return self.h(row, col); },
+           typename T::Scalar {
+             check_index(self, row, col);
+
+             return self.h(row, col);
+           },
            "row"_a, "col"_a)
-        .def("set_h",
-             py::overload_cast<hj::index, hj::index, typename T::Scalar>(
-                 &T::set_h),
-             "row"_a, "col"_a, "value"_a)
+        .def(
+            "set_h",
+            [](T &self, hj::index row, hj::index col,
+               typename T::Scalar value) {
+              check_index(self, row, col);
+
+              self.set_h(row, col, value);
+            },
+            "row"_a, "col"_a, "value"_a)
         .def("hm", py::overload_cast<std::string>(&T::hm, py::const_),
              "mode"_a = "full")
         .def("set_hm", &T::set_hm, "value"_a);

python/tests/test_DDScalar.py

@@ -565,6 +565,79 @@ def test_hessian_access(ctx):
         u.set_hm([[1, 0], [0, 1]])
 
 
+@pytest.mark.parametrize("ctx", **test_data)
+def test_variable_with_index_out_of_range(ctx):
+    for i in [2, -1]:
+        if ctx.dtype.is_dynamic:
+            with pytest.raises(RuntimeError):
+                ctx.dtype.variable(i=i, f=1, size=2)
+        else:
+            with pytest.raises(RuntimeError):
+                ctx.dtype.variable(i=i, f=1)
+
+
+@pytest.mark.parametrize("ctx", **test_data)
+def test_negative_size(ctx):
+    if not ctx.dtype.is_dynamic:
+        return
+
+    with pytest.raises(RuntimeError):
+        ctx.dtype.empty(size=-1)
+
+    with pytest.raises(RuntimeError):
+        ctx.dtype.zero(size=-1)
+
+    with pytest.raises(RuntimeError):
+        ctx.dtype.constant(f=1, size=-1)
+
+    with pytest.raises(RuntimeError):
+        ctx.dtype.constant(f=1, size=2).resize(-1)
+
+
+@pytest.mark.parametrize("ctx", **test_data)
+def test_pad_below_current_size(ctx):
+    if not ctx.dtype.is_dynamic:
+        return
+
+    u = ctx.dtype.constant(f=1, size=3)
+
+    with pytest.raises(RuntimeError):
+        u.pad_right(new_size=2)
+
+    with pytest.raises(RuntimeError):
+        u.pad_left(new_size=2)
+
+
+@pytest.mark.parametrize("ctx", **test_data)
+def test_hessian_index_out_of_range(ctx):
+    if ctx.dtype.order == 1:
+        return
+
+    u = ctx.from_data([1, 2, 3, 4, 5, 6])
+
+    for row, col in [(2, 0), (0, 2), (-1, 0), (0, -1)]:
+        with pytest.raises(IndexError):
+            u.h(row, col)
+
+        with pytest.raises(IndexError):
+            u.set_h(row, col, 1)
+
+
+def test_set_hm_with_wrong_shape():
+    u = hj.DDScalar([1, 2, 3, 4, 5, 6])
+
+    with pytest.raises(RuntimeError):
+        u.set_hm(np.ones((3, 3)))
+
+    with pytest.raises(RuntimeError):
+        u.set_hm(np.ones((1, 1)))
+
+
+def test_init_by_array_with_inconsistent_shape():
+    with pytest.raises(RuntimeError):
+        hj.DDScalar(f=1, g=[1, 2, 3], hm=np.ones((2, 2)))
+
+
 def test_is_dynamic():
     assert_equal(static_set_2.u1.is_dynamic, False)
     assert_equal(dynamic_set_2.u1.is_dynamic, True)
@@ -596,6 +669,24 @@ def test_eval():
     assert_equal(u.eval([11, 12, 13]), 5031.5)
 
 
+def test_eval_with_wrong_number_of_values():
+    u = hj.DDScalar([1, 2, 3, 4, 5, 6])
+
+    with pytest.raises(RuntimeError):
+        u.eval([])
+
+    with pytest.raises(RuntimeError):
+        u.eval([1])
+
+    with pytest.raises(RuntimeError):
+        u.eval([1, 2, 3])
+
+    v = hj.DScalar([1, 2, 3])
+
+    with pytest.raises(RuntimeError):
+        v.eval([1])
+
+
 @pytest.mark.parametrize("ctx", **test_data)
 def test_pad_right(ctx):
     u = ctx.u9

test/src/test.cpp

@@ -520,6 +520,32 @@ TEST_CASE("Hessian access is restricted to second order") {
   CHECK_FALSE(HasHessianMatrix<DDScalar<1, double, Dynamic>>);
 }
 
+// For dynamic scalars the sizes of the arguments are only known at runtime, so
+// they have to be validated instead of relying on the type system.
+
+TEST_CASE("Dynamic size checks") {
+  using D = DDScalar<2, double, Dynamic>;
+
+  const auto a = D::constant(1.0, 2);
+
+  CHECK_THROWS_AS(a.eval(std::vector<double>{1.0}), std::runtime_error);
+  CHECK_THROWS_AS(a.eval(std::vector<double>{1.0, 2.0, 3.0}),
+                  std::runtime_error);
+  CHECK_NOTHROW(a.eval(std::vector<double>{1.0, 2.0}));
+
+  auto b = D::constant(1.0, 2);
+
+  Eigen::MatrixXd too_large = Eigen::MatrixXd::Ones(3, 3);
+  Eigen::MatrixXd matching = Eigen::MatrixXd::Ones(2, 2);
+
+  CHECK_THROWS_AS(b.set_hm(too_large), std::runtime_error);
+  CHECK_NOTHROW(b.set_hm(matching));
+
+  Eigen::MatrixXd out_too_small = Eigen::MatrixXd::Zero(1, 1);
+
+  CHECK_THROWS_AS(b.hm("full", out_too_small), std::runtime_error);
+}
+
 const SScalar<double> s1(3.0, {{"x", 1.0}, {"y", 6.0}, {"z", 4.0}});
 const SScalar<double> s2(4.0, {{"x", 7.0}, {"y", 1.0}});
 const SScalar<double> s3(0.3, {{"x", 0.1}, {"y", 0.8}, {"z", 0.2}});