From 08b95e51c8bd19ed811bb47e0349ff0d970eb50a Mon Sep 17 00:00:00 2001
From: 40% <Zeroto521@gmail.com>
Date: Mon, 24 Nov 2025 16:20:29 +0800
Subject: [PATCH 01/20] Add test for matrix sum return type

Adds a test to ensure that summing a matrix variable along an axis returns a MatrixExpr type instead of MatrixVariable, addressing issue #1117.
---
 tests/test_matrix_variable.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/tests/test_matrix_variable.py b/tests/test_matrix_variable.py
index 27f549000..5119ffa18 100644
--- a/tests/test_matrix_variable.py
+++ b/tests/test_matrix_variable.py
@@ -521,6 +521,14 @@ def test_matrix_matmul_return_type():
     assert type(y @ z) is MatrixExpr
 
 
+def test_matrix_sum_return_type():
+    # test #1117, require returning type is MatrixExpr not MatrixVariable
+    m = Model()
+
+    x = m.addMatrixVar((3, 2))
+    assert type(x.sum(axis=1)) is MatrixExpr
+
+
 def test_broadcast():
     # test #1065
     m = Model()

From 0d26107451a59d1ed55323a77e399c6579fe3a92 Mon Sep 17 00:00:00 2001
From: 40% <Zeroto521@gmail.com>
Date: Mon, 24 Nov 2025 16:20:46 +0800
Subject: [PATCH 02/20] return `MatrixExpr` type

---
 src/pyscipopt/matrix.pxi | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pyscipopt/matrix.pxi b/src/pyscipopt/matrix.pxi
index 8353ed767..b34ce2325 100644
--- a/src/pyscipopt/matrix.pxi
+++ b/src/pyscipopt/matrix.pxi
@@ -53,7 +53,7 @@ class MatrixExpr(np.ndarray):
         if kwargs.get("axis") is None:
             # Speed up `.sum()` #1070
             return quicksum(self.flat)
-        return super().sum(**kwargs)
+        return super().sum(**kwargs).view(MatrixExpr)
 
     def __le__(self, other: Union[float, int, "Expr", np.ndarray, "MatrixExpr"]) -> MatrixExprCons:
         return _matrixexpr_richcmp(self, other, 1)

From 160435a416fffeed602cfdccb8918351e488c953 Mon Sep 17 00:00:00 2001
From: 40% <Zeroto521@gmail.com>
Date: Mon, 24 Nov 2025 16:22:07 +0800
Subject: [PATCH 03/20] Update CHANGELOG.md

---
 CHANGELOG.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index cddd24531..28e6c067d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,7 @@
 ## Unreleased
 ### Added
 ### Fixed
+- Fixed the type of MatrixExpr.sum(axis=...) result from MatrixVariable to MatrixExpr.
 ### Changed
 ### Removed
 

From ede6a42e2a62175bc452919bac7a9106160c4b47 Mon Sep 17 00:00:00 2001
From: 40% <Zeroto521@gmail.com>
Date: Wed, 17 Dec 2025 23:16:44 +0800
Subject: [PATCH 04/20] Speed up `MatrixExpr.sum(axis=...)`

Enhanced the MatrixExpr.sum method to accept axis as int or tuple, handle keepdims, and provide better error checking for axis bounds. This improves compatibility with numpy's sum behavior and allows more flexible summation over matrix expressions.
---
 src/pyscipopt/matrix.pxi | 63 ++++++++++++++++++++++++++++++++++++----
 1 file changed, 57 insertions(+), 6 deletions(-)

diff --git a/src/pyscipopt/matrix.pxi b/src/pyscipopt/matrix.pxi
index 8353ed767..c1762f5a9 100644
--- a/src/pyscipopt/matrix.pxi
+++ b/src/pyscipopt/matrix.pxi
@@ -4,7 +4,7 @@
 """
 
 import numpy as np
-from typing import Union
+from typing import Optional, Union
 
 
 def _is_number(e):
@@ -44,16 +44,67 @@ def _matrixexpr_richcmp(self, other, op):
 
 
 class MatrixExpr(np.ndarray):
-    def sum(self, **kwargs):
+
+    def sum(
+        self,
+        axis: Optional[tuple[int]] = None,
+        keepdims: bool = False,
+        **kwargs,
+    ) -> Union[Expr, MatrixExpr]:
         """
         Based on `numpy.ndarray.sum`, but returns a scalar if `axis=None`.
         This is useful for matrix expressions to compare with a matrix or a scalar.
         """
 
-        if kwargs.get("axis") is None:
-            # Speed up `.sum()` #1070
-            return quicksum(self.flat)
-        return super().sum(**kwargs)
+        if axis is None:
+            axis = tuple(range(self.ndim))
+
+        elif isinstance(axis, int):
+            if axis < -self.ndim or axis >= self.ndim:
+                raise np.exceptions.AxisError(
+                    f"axis {axis} is out of bounds for array of dimension {self.ndim}"
+                )
+            axis = (axis,)
+
+        elif isinstance(axis, tuple) and all(isinstance(i, int) for i in axis):
+            for i in axis:
+                if i < -self.ndim or i >= self.ndim:
+                    raise np.exceptions.AxisError(
+                        f"axis {i} is out of bounds for array of dimension {self.ndim}"
+                    )
+
+        else:
+            raise TypeError("'axis' must be an int or a tuple of ints")
+
+        if len(axis := tuple(i + self.ndim if i < 0 else i for i in axis)) == self.ndim:
+            res = quicksum(self.flat)
+            if keepdims:
+                return (
+                    np.array([res], dtype=object)
+                    .reshape([1] * self.ndim)
+                    .view(MatrixExpr)
+                )
+            return res
+
+        keep_axes = tuple(i for i in range(self.ndim) if i not in axis)
+        shape = (
+            tuple(1 if i in axis else self.shape[i] for i in range(self.ndim))
+            if keepdims
+            else tuple(self.shape[i] for i in keep_axes)
+        )
+        return (
+            np.fromiter(
+                map(
+                    quicksum,
+                    self.transpose(keep_axes + axis).reshape(
+                        -1, np.prod([self.shape[i] for i in axis])
+                    ),
+                ),
+                dtype=object,
+            )
+            .reshape(shape)
+            .view(MatrixExpr)
+        )
 
     def __le__(self, other: Union[float, int, "Expr", np.ndarray, "MatrixExpr"]) -> MatrixExprCons:
         return _matrixexpr_richcmp(self, other, 1)

From 3a757bb330b6bf9e2629b47a709a50adfef6e1a8 Mon Sep 17 00:00:00 2001
From: 40% <Zeroto521@gmail.com>
Date: Thu, 18 Dec 2025 09:08:00 +0800
Subject: [PATCH 05/20] Test `MatrixExpr.sum(axis=...)`

Removed n=200 from the sum performance test to limit test size. Added additional performance assertions comparing np.ndarray.sum and the optimized sum method for matrix variables.
---
 tests/test_matrix_variable.py | 14 +++++++++++++-
 1 file changed, 13 insertions(+), 1 deletion(-)

diff --git a/tests/test_matrix_variable.py b/tests/test_matrix_variable.py
index 5119ffa18..5c35d2f4e 100644
--- a/tests/test_matrix_variable.py
+++ b/tests/test_matrix_variable.py
@@ -211,7 +211,7 @@ def test_matrix_sum_argument():
     assert (m.getVal(y) == np.full((2, 4), 3)).all().all()
 
 
-@pytest.mark.parametrize("n", [50, 100, 200])
+@pytest.mark.parametrize("n", [50, 100])
 def test_sum_performance(n):
     model = Model()
     x = model.addMatrixVar((n, n))
@@ -228,6 +228,18 @@ def test_sum_performance(n):
 
     assert model.isGT(end_orig - start_orig, end_matrix - start_matrix)
 
+    # Original sum via `np.ndarray.sum`, `np.sum` will call subclass method
+    start_orig = time()
+    np.ndarray.sum(x, axis=0)
+    end_orig = time()
+
+    # Optimized sum via `quicksum`
+    start_matrix = time()
+    x.sum(axis=0)
+    end_matrix = time()
+
+    assert model.isGT(end_orig - start_orig, end_matrix - start_matrix)
+
 
 def test_add_cons_matrixVar():
     m = Model()

From a13b7ec3dc3df9f6985404598d0f13304e3acba2 Mon Sep 17 00:00:00 2001
From: 40% <Zeroto521@gmail.com>
Date: Thu, 18 Dec 2025 09:14:29 +0800
Subject: [PATCH 06/20] Test `MatrixExpr.sum(axis=tuple(range(ndim))`

Renamed test_matrix_sum_argument to test_matrix_sum_axis and updated tests to use explicit axis arguments in sum operations. This clarifies the behavior when summing over all axes and improves test coverage for axis handling.
---
 tests/test_matrix_variable.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/tests/test_matrix_variable.py b/tests/test_matrix_variable.py
index 5c35d2f4e..2b1e99e09 100644
--- a/tests/test_matrix_variable.py
+++ b/tests/test_matrix_variable.py
@@ -181,7 +181,8 @@ def test_expr_from_matrix_vars():
         for term, coeff in expr_list:
             assert len(term) == 3
 
-def test_matrix_sum_argument():
+
+def test_matrix_sum_axis():
     m = Model()
 
     # Return a array when axis isn't None
@@ -190,7 +191,8 @@ def test_matrix_sum_argument():
 
     # compare the result of summing 2d array to a scalar with a scalar
     x = m.addMatrixVar((2, 3), "x", "I", ub=4)
-    m.addMatrixCons(x.sum() == 24)
+    # `axis=tuple(range(x.ndim))` is `axis=None`
+    m.addMatrixCons(x.sum(axis=tuple(range(x.ndim))) == 24)
 
     # compare the result of summing 2d array to 1d array
     y = m.addMatrixVar((2, 4), "y", "I", ub=4)
@@ -204,7 +206,7 @@ def test_matrix_sum_argument():
     # to fix the element values
     m.addMatrixCons(z == np.ones((2, 3, 4)))
 
-    m.setObjective(x.sum() + y.sum() + z.sum(), "maximize")
+    m.setObjective(x.sum() + y.sum() + z.sum(tuple(range(z.ndim))), "maximize")
     m.optimize()
 
     assert (m.getVal(x) == np.full((2, 3), 4)).all().all()

From f9b4a1671b1498cb950cf8a8092444b7d01cfb7e Mon Sep 17 00:00:00 2001
From: 40% <Zeroto521@gmail.com>
Date: Thu, 18 Dec 2025 09:55:06 +0800
Subject: [PATCH 07/20] Refactor MatrixExpr.sum axis handling with numpy
 utility

Replaces manual axis validation and normalization in MatrixExpr.sum with numpy's normalize_axis_tuple for improved reliability and code clarity. Updates type hints and simplifies logic for summing across all axes.
---
 src/pyscipopt/matrix.pxi | 43 ++++++++++++----------------------------
 1 file changed, 13 insertions(+), 30 deletions(-)

diff --git a/src/pyscipopt/matrix.pxi b/src/pyscipopt/matrix.pxi
index c1762f5a9..9d377780c 100644
--- a/src/pyscipopt/matrix.pxi
+++ b/src/pyscipopt/matrix.pxi
@@ -3,9 +3,11 @@
 # TODO Add tests
 """
 
-import numpy as np
 from typing import Optional, Union
 
+import numpy as np
+from numpy.lib.array_utils import normalize_axis_tuple
+
 
 def _is_number(e):
     try:
@@ -47,7 +49,7 @@ class MatrixExpr(np.ndarray):
 
     def sum(
         self,
-        axis: Optional[tuple[int]] = None,
+        axis: Optional[Union[int, tuple[int, ...]]] = None,
         keepdims: bool = False,
         **kwargs,
     ) -> Union[Expr, MatrixExpr]:
@@ -56,35 +58,16 @@ class MatrixExpr(np.ndarray):
         This is useful for matrix expressions to compare with a matrix or a scalar.
         """
 
-        if axis is None:
-            axis = tuple(range(self.ndim))
-
-        elif isinstance(axis, int):
-            if axis < -self.ndim or axis >= self.ndim:
-                raise np.exceptions.AxisError(
-                    f"axis {axis} is out of bounds for array of dimension {self.ndim}"
-                )
-            axis = (axis,)
-
-        elif isinstance(axis, tuple) and all(isinstance(i, int) for i in axis):
-            for i in axis:
-                if i < -self.ndim or i >= self.ndim:
-                    raise np.exceptions.AxisError(
-                        f"axis {i} is out of bounds for array of dimension {self.ndim}"
-                    )
-
-        else:
-            raise TypeError("'axis' must be an int or a tuple of ints")
-
-        if len(axis := tuple(i + self.ndim if i < 0 else i for i in axis)) == self.ndim:
+        axis = normalize_axis_tuple(
+            range(self.ndim) if axis is None else axis, self.ndim
+        )
+        if len(axis) == self.ndim:
             res = quicksum(self.flat)
-            if keepdims:
-                return (
-                    np.array([res], dtype=object)
-                    .reshape([1] * self.ndim)
-                    .view(MatrixExpr)
-                )
-            return res
+            return (
+                np.array([res], dtype=object).reshape([1] * self.ndim).view(MatrixExpr)
+                if keepdims
+                else res
+            )
 
         keep_axes = tuple(i for i in range(self.ndim) if i not in axis)
         shape = (

From fe874d91802514086a70adbe6e6d519ac4744d53 Mon Sep 17 00:00:00 2001
From: 40% <Zeroto521@gmail.com>
Date: Thu, 18 Dec 2025 09:55:38 +0800
Subject: [PATCH 08/20] Add tests for matrix sum error

Added tests to verify error handling in matrix variable sum operations for invalid axis types, out-of-range values, and duplicate axes.
---
 tests/test_matrix_variable.py | 21 +++++++++++++++++++++
 1 file changed, 21 insertions(+)

diff --git a/tests/test_matrix_variable.py b/tests/test_matrix_variable.py
index 2b1e99e09..da42ee1fb 100644
--- a/tests/test_matrix_variable.py
+++ b/tests/test_matrix_variable.py
@@ -182,6 +182,27 @@ def test_expr_from_matrix_vars():
             assert len(term) == 3
 
 
+def test_matrix_sum_error():
+    m = Model()
+    x = m.addMatrixVar((2, 3), "x", "I", ub=4)
+
+    # test axis type
+    with pytest.raises(TypeError):
+        x.sum("0")
+
+    # test axis value (out of range)
+    with pytest.raises(ValueError):
+        x.sum(2)
+
+    # test axis value (out of range)
+    with pytest.raises(ValueError):
+        x.sum((-3,))
+
+    # test axis value (duplicate)
+    with pytest.raises(ValueError):
+        x.sum((0, 0))
+
+
 def test_matrix_sum_axis():
     m = Model()
 

From fa7c88f375a8381e00ab1c92310c72f22255f6c4 Mon Sep 17 00:00:00 2001
From: 40% <Zeroto521@gmail.com>
Date: Thu, 18 Dec 2025 09:55:47 +0800
Subject: [PATCH 09/20] Add tests for matrix sum with keepdims parameter

Introduces test cases to verify the shape of matrix variable sums when using the keepdims argument, ensuring correct behavior for both full and axis-specific summation.
---
 tests/test_matrix_variable.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/tests/test_matrix_variable.py b/tests/test_matrix_variable.py
index da42ee1fb..ebebd9411 100644
--- a/tests/test_matrix_variable.py
+++ b/tests/test_matrix_variable.py
@@ -234,6 +234,14 @@ def test_matrix_sum_axis():
     assert (m.getVal(y) == np.full((2, 4), 3)).all().all()
 
 
+def test_matrix_sum_keepdims():
+    m = Model()
+    x = m.addMatrixVar((1, 2, 3), "x", "I", ub=4)
+
+    assert x.sum(keepdims=True).shape == (1, 1, 1)
+    assert x.sum(axis=1, keepdims=True).shape == (1, 1, 3)
+
+
 @pytest.mark.parametrize("n", [50, 100])
 def test_sum_performance(n):
     model = Model()

From a58feae9ecd4612224f1b568ecae4c9e69da2f6f Mon Sep 17 00:00:00 2001
From: 40% <Zeroto521@gmail.com>
Date: Thu, 18 Dec 2025 10:04:46 +0800
Subject: [PATCH 10/20] call `.sum` via positional argument

Replaces the use of the 'axis' keyword argument with a positional argument in the z.sum() method call to align with the expected function signature.
---
 tests/test_matrix_variable.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_matrix_variable.py b/tests/test_matrix_variable.py
index ebebd9411..35c703561 100644
--- a/tests/test_matrix_variable.py
+++ b/tests/test_matrix_variable.py
@@ -221,7 +221,7 @@ def test_matrix_sum_axis():
 
     # compare the result of summing 3d array to a 2d array with a 2d array
     z = m.addMatrixVar((2, 3, 4), "z", "I", ub=4)
-    m.addMatrixCons(z.sum(axis=2) == x)
+    m.addMatrixCons(z.sum(2) == x)
     m.addMatrixCons(z.sum(axis=1) == y)
 
     # to fix the element values

From fad458c7fa95e45f3cdd2c291e494291d7209605 Mon Sep 17 00:00:00 2001
From: 40% <Zeroto521@gmail.com>
Date: Thu, 18 Dec 2025 10:56:09 +0800
Subject: [PATCH 11/20] Refactor sum method to use np.apply_along_axis

Replaces np.fromiter with np.apply_along_axis for summing along specified axes in MatrixExpr. This simplifies the code and improves readability.
---
 src/pyscipopt/matrix.pxi | 12 +++++-------
 1 file changed, 5 insertions(+), 7 deletions(-)

diff --git a/src/pyscipopt/matrix.pxi b/src/pyscipopt/matrix.pxi
index 9d377780c..e5ba90920 100644
--- a/src/pyscipopt/matrix.pxi
+++ b/src/pyscipopt/matrix.pxi
@@ -76,14 +76,12 @@ class MatrixExpr(np.ndarray):
             else tuple(self.shape[i] for i in keep_axes)
         )
         return (
-            np.fromiter(
-                map(
-                    quicksum,
-                    self.transpose(keep_axes + axis).reshape(
-                        -1, np.prod([self.shape[i] for i in axis])
-                    ),
+            np.apply_along_axis(
+                quicksum,
+                -1,
+                self.transpose(keep_axes + axis).reshape(
+                    -1, np.prod([self.shape[i] for i in axis])
                 ),
-                dtype=object,
             )
             .reshape(shape)
             .view(MatrixExpr)

From 363c8b1ba5e11977b8fb3294fff3a849c1cb5f21 Mon Sep 17 00:00:00 2001
From: 40% <Zeroto521@gmail.com>
Date: Thu, 18 Dec 2025 11:07:55 +0800
Subject: [PATCH 12/20] Directly test the `.sum` result

Added a comment to clarify the purpose of the test_matrix_sum_axis function, indicating it compares the result of summing a matrix variable after optimization.
---
 tests/test_matrix_variable.py | 29 ++++++++++++++++++++++-------
 1 file changed, 22 insertions(+), 7 deletions(-)

diff --git a/tests/test_matrix_variable.py b/tests/test_matrix_variable.py
index 35c703561..3527328fd 100644
--- a/tests/test_matrix_variable.py
+++ b/tests/test_matrix_variable.py
@@ -19,7 +19,7 @@
     sin,
     sqrt,
 )
-from pyscipopt.scip import GenExpr
+from pyscipopt.scip import CONST, GenExpr
 
 
 def test_catching_errors():
@@ -204,6 +204,7 @@ def test_matrix_sum_error():
 
 
 def test_matrix_sum_axis():
+    # compare the result of summing matrix variable after optimization
     m = Model()
 
     # Return a array when axis isn't None
@@ -234,12 +235,26 @@ def test_matrix_sum_axis():
     assert (m.getVal(y) == np.full((2, 4), 3)).all().all()
 
 
-def test_matrix_sum_keepdims():
-    m = Model()
-    x = m.addMatrixVar((1, 2, 3), "x", "I", ub=4)
-
-    assert x.sum(keepdims=True).shape == (1, 1, 1)
-    assert x.sum(axis=1, keepdims=True).shape == (1, 1, 3)
+@pytest.mark.parametrize(
+    "axis, keepdims",
+    [
+        (0, False),
+        (0, True),
+        (1, False),
+        (1, True),
+        ((0, 2), False),
+        ((0, 2), True),
+    ],
+)
+def test_matrix_sum_result(axis, keepdims):
+    # directly compare the result of np.sum and MatrixExpr.sum
+    _getVal = np.vectorize(lambda e: e.terms[CONST])
+    a = np.arange(6).reshape((1, 2, 3))
+
+    np_res = a.sum(axis, keepdims=keepdims)
+    scip_res = MatrixExpr.sum(a, axis, keepdims=keepdims)
+    assert (np_res == _getVal(scip_res)).all()
+    assert np_res.shape == _getVal(scip_res).shape
 
 
 @pytest.mark.parametrize("n", [50, 100])

From 0e03a3b715ef095f71bf26a60925f49b27ee5887 Mon Sep 17 00:00:00 2001
From: 40% <Zeroto521@gmail.com>
Date: Thu, 18 Dec 2025 11:09:53 +0800
Subject: [PATCH 13/20] Update CHANGELOG.md

---
 CHANGELOG.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index db876adaa..78b1b40e3 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Unreleased
 ### Added
+- Speed up MatrixExpr.sum(axis=...) via quicksum
 ### Fixed
 - all fundamental callbacks now raise an error if not implemented
 - Fixed the type of MatrixExpr.sum(axis=...) result from MatrixVariable to MatrixExpr.

From 4c0bb703a0f8c0ba7ae75d768e288bdec3b5c865 Mon Sep 17 00:00:00 2001
From: 40% <Zeroto521@gmail.com>
Date: Thu, 18 Dec 2025 11:23:00 +0800
Subject: [PATCH 14/20] Expand docstring for MatrixExpr.sum method

The docstring for the MatrixExpr.sum method was updated to provide detailed information about parameters, return values, and behavior, improving clarity and alignment with numpy conventions.
---
 src/pyscipopt/matrix.pxi | 32 +++++++++++++++++++++++++++++---
 1 file changed, 29 insertions(+), 3 deletions(-)

diff --git a/src/pyscipopt/matrix.pxi b/src/pyscipopt/matrix.pxi
index e5ba90920..b8e8753c7 100644
--- a/src/pyscipopt/matrix.pxi
+++ b/src/pyscipopt/matrix.pxi
@@ -54,10 +54,36 @@ class MatrixExpr(np.ndarray):
         **kwargs,
     ) -> Union[Expr, MatrixExpr]:
         """
-        Based on `numpy.ndarray.sum`, but returns a scalar if `axis=None`.
-        This is useful for matrix expressions to compare with a matrix or a scalar.
+        Return the sum of the array elements over the given axis.
+
+        Parameters
+        ----------
+        axis : None or int or tuple of ints, optional
+            Axis or axes along which a sum is performed. The default, axis=None, will
+            sum all of the elements of the input array. If axis is negative it counts
+            from the last to the first axis. If axis is a tuple of ints, a sum is
+            performed on all of the axes specified in the tuple instead of a single axis
+            or all the axes as before.
+
+        keepdims : bool, optional
+            If this is set to True, the axes which are reduced are left in the result as
+            dimensions with size one. With this option, the result will broadcast
+            correctly against the input array.
+
+        **kwargs : ignored
+            Additional keyword arguments are ignored. They exist for compatibility
+            with `numpy.ndarray.sum`.
+
+        Returns
+        -------
+        Expr
+            If the sum is performed over all axes, return an Expr.
+
+        MatrixExpr
+            - If the sum is performed over a subset of axes return a MatrixExpr.
+            - If `keepdims` is True, the returned MatrixExpr will have the same number 
+              of dimensions as the original array, with the reduced axes having size one.
         """
-
         axis = normalize_axis_tuple(
             range(self.ndim) if axis is None else axis, self.ndim
         )

From 8e0e372d010cc5fc66113f2bff662e58e5cb6cfb Mon Sep 17 00:00:00 2001
From: 40% <Zeroto521@gmail.com>
Date: Thu, 18 Dec 2025 11:30:25 +0800
Subject: [PATCH 15/20] Clarify MatrixExpr.sum docstring and note quicksum
 usage

Updated the docstring for MatrixExpr.sum to specify that it uses quicksum for speed optimization instead of numpy.ndarray.sum. Added a detailed note explaining the difference between quicksum (using __iadd__) and numpy's sum (using __add__).
---
 src/pyscipopt/matrix.pxi | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/src/pyscipopt/matrix.pxi b/src/pyscipopt/matrix.pxi
index b8e8753c7..e163a4789 100644
--- a/src/pyscipopt/matrix.pxi
+++ b/src/pyscipopt/matrix.pxi
@@ -54,7 +54,8 @@ class MatrixExpr(np.ndarray):
         **kwargs,
     ) -> Union[Expr, MatrixExpr]:
         """
-        Return the sum of the array elements over the given axis.
+        Return the sum of the array elements over the given axis. Speed optimized for
+        matrix expressions via `quicksum` instead of `numpy.ndarray.sum`.
 
         Parameters
         ----------
@@ -83,6 +84,12 @@ class MatrixExpr(np.ndarray):
             - If the sum is performed over a subset of axes return a MatrixExpr.
             - If `keepdims` is True, the returned MatrixExpr will have the same number 
               of dimensions as the original array, with the reduced axes having size one.
+
+        Notes
+        -----
+        `quicksum` uses `__iadd__` to accumulate the sum, which modifies an existing
+        object in place. `numpy.ndarray.sum` uses `__add__`, which creates a new object
+        for each addition.
         """
         axis = normalize_axis_tuple(
             range(self.ndim) if axis is None else axis, self.ndim

From 4931f06854e71931b2d4940bf1ce74b7a80f3a3e Mon Sep 17 00:00:00 2001
From: 40% <Zeroto521@gmail.com>
Date: Thu, 18 Dec 2025 13:03:16 +0800
Subject: [PATCH 16/20] Split up two costing time test cases

Renamed the existing performance test to clarify it tests the case where axis is None. Added a new test to measure performance when summing along a specific axis.
---
 tests/test_matrix_variable.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/tests/test_matrix_variable.py b/tests/test_matrix_variable.py
index 3527328fd..e4758f077 100644
--- a/tests/test_matrix_variable.py
+++ b/tests/test_matrix_variable.py
@@ -258,7 +258,7 @@ def test_matrix_sum_result(axis, keepdims):
 
 
 @pytest.mark.parametrize("n", [50, 100])
-def test_sum_performance(n):
+def test_matrix_sum_axis_is_none_performance(n):
     model = Model()
     x = model.addMatrixVar((n, n))
 
@@ -274,6 +274,12 @@ def test_sum_performance(n):
 
     assert model.isGT(end_orig - start_orig, end_matrix - start_matrix)
 
+
+@pytest.mark.parametrize("n", [50, 100])
+def test_matrix_sum_axis_not_none_performance(n):
+    model = Model()
+    x = model.addMatrixVar((n, n))
+
     # Original sum via `np.ndarray.sum`, `np.sum` will call subclass method
     start_orig = time()
     np.ndarray.sum(x, axis=0)

From b89a111ca50d346c8697065048f689466693ff6b Mon Sep 17 00:00:00 2001
From: 40% <Zeroto521@gmail.com>
Date: Thu, 18 Dec 2025 18:11:03 +0800
Subject: [PATCH 17/20] Supports numpy 1.x

Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
---
 src/pyscipopt/matrix.pxi | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/src/pyscipopt/matrix.pxi b/src/pyscipopt/matrix.pxi
index e163a4789..0c989b04d 100644
--- a/src/pyscipopt/matrix.pxi
+++ b/src/pyscipopt/matrix.pxi
@@ -6,7 +6,12 @@
 from typing import Optional, Union
 
 import numpy as np
-from numpy.lib.array_utils import normalize_axis_tuple
+try:
+    # NumPy 2.x location
+    from numpy.lib.array_utils import normalize_axis_tuple
+except ImportError:
+    # Fallback for NumPy 1.x
+    from numpy.core.numeric import normalize_axis_tuple
 
 
 def _is_number(e):

From 5a4c01fc996eb4d6a719b7f02ff33ee1547c035b Mon Sep 17 00:00:00 2001
From: 40% <Zeroto521@gmail.com>
Date: Thu, 18 Dec 2025 18:12:43 +0800
Subject: [PATCH 18/20] Supports Python 3.8

---
 src/pyscipopt/matrix.pxi | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/pyscipopt/matrix.pxi b/src/pyscipopt/matrix.pxi
index e163a4789..921824595 100644
--- a/src/pyscipopt/matrix.pxi
+++ b/src/pyscipopt/matrix.pxi
@@ -3,7 +3,7 @@
 # TODO Add tests
 """
 
-from typing import Optional, Union
+from typing import Optional, Tuple, Union
 
 import numpy as np
 from numpy.lib.array_utils import normalize_axis_tuple
@@ -49,7 +49,7 @@ class MatrixExpr(np.ndarray):
 
     def sum(
         self,
-        axis: Optional[Union[int, tuple[int, ...]]] = None,
+        axis: Optional[Union[int, Tuple[int, ...]]] = None,
         keepdims: bool = False,
         **kwargs,
     ) -> Union[Expr, MatrixExpr]:

From 8ac4bf05101c0898dfe8d97014ba2ad44fe59454 Mon Sep 17 00:00:00 2001
From: 40% <Zeroto521@gmail.com>
Date: Thu, 18 Dec 2025 21:03:52 +0800
Subject: [PATCH 19/20] Simplify a bit

Refactored the sum method in MatrixExpr to clarify axis typing and simplify the application of np.apply_along_axis. The new implementation improves readability and maintains the intended behavior.
---
 src/pyscipopt/matrix.pxi | 16 ++++------------
 1 file changed, 4 insertions(+), 12 deletions(-)

diff --git a/src/pyscipopt/matrix.pxi b/src/pyscipopt/matrix.pxi
index e54df802f..90ca74a5e 100644
--- a/src/pyscipopt/matrix.pxi
+++ b/src/pyscipopt/matrix.pxi
@@ -96,7 +96,7 @@ class MatrixExpr(np.ndarray):
         object in place. `numpy.ndarray.sum` uses `__add__`, which creates a new object
         for each addition.
         """
-        axis = normalize_axis_tuple(
+        axis: Tuple[int, ...] = normalize_axis_tuple(
             range(self.ndim) if axis is None else axis, self.ndim
         )
         if len(axis) == self.ndim:
@@ -113,17 +113,9 @@ class MatrixExpr(np.ndarray):
             if keepdims
             else tuple(self.shape[i] for i in keep_axes)
         )
-        return (
-            np.apply_along_axis(
-                quicksum,
-                -1,
-                self.transpose(keep_axes + axis).reshape(
-                    -1, np.prod([self.shape[i] for i in axis])
-                ),
-            )
-            .reshape(shape)
-            .view(MatrixExpr)
-        )
+        return np.apply_along_axis(
+            quicksum, -1, self.transpose(keep_axes + axis).reshape(shape + (-1,))
+        ).view(MatrixExpr)
 
     def __le__(self, other: Union[float, int, "Expr", np.ndarray, "MatrixExpr"]) -> MatrixExprCons:
         return _matrixexpr_richcmp(self, other, 1)

From c11743d131c7543570048f1e0e8915aa4cf31df0 Mon Sep 17 00:00:00 2001
From: Joao-Dionisio <up201606210@up.pt>
Date: Thu, 25 Dec 2025 13:48:43 +0100
Subject: [PATCH 20/20] suggestion

---
 src/pyscipopt/matrix.pxi | 21 +++++----------------
 1 file changed, 5 insertions(+), 16 deletions(-)

diff --git a/src/pyscipopt/matrix.pxi b/src/pyscipopt/matrix.pxi
index 90ca74a5e..1a6a09cf3 100644
--- a/src/pyscipopt/matrix.pxi
+++ b/src/pyscipopt/matrix.pxi
@@ -4,7 +4,6 @@
 """
 
 from typing import Optional, Tuple, Union
-
 import numpy as np
 try:
     # NumPy 2.x location
@@ -59,8 +58,7 @@ class MatrixExpr(np.ndarray):
         **kwargs,
     ) -> Union[Expr, MatrixExpr]:
         """
-        Return the sum of the array elements over the given axis. Speed optimized for
-        matrix expressions via `quicksum` instead of `numpy.ndarray.sum`.
+        Return the sum of the array elements over the given axis.
 
         Parameters
         ----------
@@ -82,19 +80,10 @@ class MatrixExpr(np.ndarray):
 
         Returns
         -------
-        Expr
-            If the sum is performed over all axes, return an Expr.
-
-        MatrixExpr
-            - If the sum is performed over a subset of axes return a MatrixExpr.
-            - If `keepdims` is True, the returned MatrixExpr will have the same number 
-              of dimensions as the original array, with the reduced axes having size one.
-
-        Notes
-        -----
-        `quicksum` uses `__iadd__` to accumulate the sum, which modifies an existing
-        object in place. `numpy.ndarray.sum` uses `__add__`, which creates a new object
-        for each addition.
+        Expr or MatrixExpr
+            If the sum is performed over all axes, return an Expr, otherwise return
+            a MatrixExpr.
+
         """
         axis: Tuple[int, ...] = normalize_axis_tuple(
             range(self.ndim) if axis is None else axis, self.ndim