Various improvements (#100)

* various improvements

* fix test

* fix test again

* add test

Author: Korijn

pylinalg/matrix.py

@@ -737,8 +737,7 @@ def _mat_inv(m) -> np.ndarray:
     det = n11 * t11 + n21 * t12 + n31 * t13 + n41 * t14
 
     if det == 0:
-        out.fill(0)
-        return out  # singular matrix
+        raise np.linalg.LinAlgError("Singular matrix")
 
     det_inv = 1 / det
 
@@ -853,11 +852,11 @@ def _mat_inv(m) -> np.ndarray:
 if int(np.__version__.split(".")[0]) >= 2:
     _default_mat_inv_method = "numpy"
 else:
-    _default_mat_inv_method = "manual"
+    _default_mat_inv_method = "python"
 
 
 def mat_inverse(
-    matrix, /, *, method=_default_mat_inv_method, dtype=None, out=None
+    matrix, /, *, method=_default_mat_inv_method, raise_err=False, dtype=None, out=None
 ) -> np.ndarray:
     """
     Compute the inverse of a matrix.
@@ -868,7 +867,9 @@ def mat_inverse(
         The matrix to invert.
     method : str, optional
         The method to use for inversion. The default is "numpy" when
-        numpy version is 2.0.0 or newer, otherwise "manual".
+        numpy version is 2.0.0 or newer, otherwise "python".
+    raise_err : bool, optional
+        Raise a ValueError if the matrix is singular. Default is False.
     dtype : data-type, optional
         Overrides the data type of the result.
     out : ndarray, optional
@@ -886,11 +887,11 @@ def mat_inverse(
     -----
     The default method is "numpy" when numpy version >= 2.0.0,
     which uses the `numpy.linalg.inv` function.
-    The alternative method is "manual", which uses a manual implementation of
-    the inversion algorithm. The manual method is used to avoid a performance
+    The alternative method is "python", which uses a pure python implementation of
+    the inversion algorithm. The python method is used to avoid a performance
     issue with `numpy.linalg.inv` on some platforms when numpy version < 2.0.0.
     See: https://github.com/pygfx/pygfx/issues/763
-    The manual method is slower than the numpy method, but it is guaranteed to work.
+    The python method is slower than the numpy method, but it is guaranteed to work.
 
     When the matrix is singular, it will return a matrix filled with zeros,
     This is a common behavior in real-time graphics applications.
@@ -899,18 +900,18 @@ def mat_inverse(
 
     fn = {
         "numpy": np.linalg.inv,
-        "manual": _mat_inv,
+        "python": _mat_inv,
     }[method]
 
     matrix = np.asarray(matrix)
     try:
         inverse = fn(matrix)
-    except np.linalg.LinAlgError:
+    except np.linalg.LinAlgError as err:
+        if raise_err:
+            raise ValueError("The provided matrix is not invertible.") from err
         inverse = np.zeros_like(matrix, dtype=dtype)
     if out is None:
-        if dtype is not None:
-            return inverse.astype(dtype, copy=False)
-        return inverse
+        return np.asarray(inverse, dtype=dtype)
     else:
         out[:] = inverse
     return out

pylinalg/misc.py

@@ -131,7 +131,7 @@ def quat_to_axis_angle(quaternion, /, *, out=None, dtype=None) -> np.ndarray:
     quaternion = np.asarray(quaternion)
 
     if out is None:
-        quaternion = quaternion.astype(dtype)
+        quaternion = quaternion.astype(dtype, copy=False)
         out = (
             quaternion[..., :3] / np.sqrt(1 - quaternion[..., 3] ** 2),
             2 * np.arccos(quaternion[..., 3]),

pylinalg/vector.py

@@ -115,7 +115,9 @@ def vec_transform(vectors, matrix, /, *, w=1, out=None, dtype=None) -> np.ndarra
     return out
 
 
-def vec_unproject(vector, matrix, /, *, depth=0, out=None, dtype=None) -> np.ndarray:
+def vec_unproject(
+    vector, matrix, /, *, matrix_is_inv=False, depth=0, out=None, dtype=None
+) -> np.ndarray:
     """
     Un-project a vector from 2D space to 3D space.
 
@@ -130,6 +132,9 @@ def vec_unproject(vector, matrix, /, *, depth=0, out=None, dtype=None) -> np.nda
         The vector to be un-projected.
     matrix: ndarray, [4, 4]
         The camera's intrinsic matrix.
+    matrix_is_inv: bool, optional
+        Default is False. If True, the provided matrix is assumed to be the
+        inverse of the camera's intrinsic matrix.
     depth : number, optional
         The distance of the unprojected vector from the camera.
     out : ndarray, optional
@@ -162,10 +167,12 @@ def vec_unproject(vector, matrix, /, *, depth=0, out=None, dtype=None) -> np.nda
     if out is None:
         out = np.empty((*result_shape, 3), dtype=dtype)
 
-    try:
-        inverse_projection = np.linalg.inv(matrix)
-    except np.linalg.LinAlgError as err:
-        raise ValueError("The provided matrix is not invertible.") from err
+    if matrix_is_inv:
+        inverse_projection = matrix
+    else:
+        from .matrix import mat_inverse
+
+        inverse_projection = mat_inverse(matrix, raise_err=True)
 
     vector_hom = np.empty((*result_shape, 4), dtype=dtype)
     vector_hom[..., 2] = depth
@@ -295,7 +302,7 @@ def vec_dist(vector_a, vector_b, /, *, out=None, dtype=None) -> np.ndarray:
 
     shape = vector_a.shape[:-1]
     if out is None:
-        out = np.linalg.norm(vector_a - vector_b, axis=-1).astype(dtype)
+        out = np.linalg.norm(vector_a - vector_b, axis=-1).astype(dtype, copy=False)
     elif len(shape) >= 0:
         out[:] = np.linalg.norm(vector_a - vector_b, axis=-1)
     else:
@@ -346,7 +353,7 @@ def vec_angle(vector_a, vector_b, /, *, out=None, dtype=None) -> np.ndarray:
     )
 
     if out is None:
-        out = np.arccos(the_cos).astype(dtype)
+        out = np.arccos(the_cos).astype(dtype, copy=False)
     elif len(shape) >= 0:
         out[:] = np.arccos(the_cos)
     else:

pyproject.toml

@@ -2,7 +2,7 @@
 
 [project]
 name = "pylinalg"
-version = "0.6.2"
+version = "0.6.3"
 description = "Linear algebra utilities for Python"
 readme = "README.md"
 license = { file = "LICENSE" }

tests/test_matrix.py

@@ -385,15 +385,15 @@ def test_mat_euler_vs_scipy():
 def test_mat_inverse():
     a = la.mat_from_translation([1, 2, 3])
     np_inv = la.mat_inverse(a, method="numpy")
-    manual_inv = la.mat_inverse(a, method="manual")
-    npt.assert_array_equal(np_inv, manual_inv)
+    python_inv = la.mat_inverse(a, method="python")
+    npt.assert_array_equal(np_inv, python_inv)
 
     # test for singular matrix
     b = la.mat_from_scale([0, 2, 3])
     np_inv = la.mat_inverse(b, method="numpy")
-    manual_inv = la.mat_inverse(b, method="manual")
+    python_inv = la.mat_inverse(b, method="python")
     npt.assert_array_equal(np_inv, np.zeros((4, 4)))
-    npt.assert_array_equal(manual_inv, np.zeros((4, 4)))
+    npt.assert_array_equal(python_inv, np.zeros((4, 4)))
 
 
 def test_mat_has_shear():

tests/test_pylinalg.py

@@ -58,9 +58,9 @@ def popattr(key):
         # confirm the signature of each callable
         sig = inspect.signature(getattr(la, key))
 
-        assert (
-            sig.return_annotation is not inspect.Signature.empty
-        ), f"Missing return annotation on {key}"
+        assert sig.return_annotation is not inspect.Signature.empty, (
+            f"Missing return annotation on {key}"
+        )
         if sig.return_annotation is bool:
             key_parts = key.split("_")
             assert key_parts[1] in ("is", "has")
@@ -77,9 +77,9 @@ def popattr(key):
                     assert param.KEYWORD_ONLY
                     has_out = True
             assert has_out, f"Function {key} does not have 'out' keyword-only argument"
-            assert (
-                has_dtype
-            ), f"Function {key} does not have 'dtype' keyword-only argument"
+            assert has_dtype, (
+                f"Function {key} does not have 'dtype' keyword-only argument"
+            )
 
     # assert that all callables are available in __all__
     assert set(__all__) == set(callables)

tests/test_vectors.py

@@ -303,6 +303,16 @@ def test_vec_unproject_exceptions():
         la.vec_unproject(vector, matrix)
 
 
+def test_vec_unproject_is_inverse():
+    a = la.mat_perspective(-1, 1, -1, 1, 1, 100)
+    a_inv = la.mat_inverse(a)
+    vecs = np.array([[1, 2], [4, 5], [7, 8]])
+
+    expected = la.vec_unproject(vecs, a)
+    actual = la.vec_unproject(vecs, a_inv, matrix_is_inv=True)
+    npt.assert_array_equal(expected, actual)
+
+
 def test_vector_apply_rotation_ordered():
     """Test that a positive pi/2 rotation about the z-axis and then the y-axis
     results in a different output then in standard rotation ordering."""