See about adding null space (#72)

* See about adding null space

* .

* .

fluids/friction.py

@@ -355,10 +355,6 @@ def Colebrook(Re, eD, tol=None):
     for hundreds of thousand of points within the region 1E-12 < Re < 1E12
     and 0 < eD < 0.1.
 
-    The numerical solution attempts the secant method using `scipy`'s `newton`
-    solver, and in the event of nonconvergence, attempts the `fsolve` solver
-    as well. An initial guess is provided via the `Clamond` function.
-
     The numerical and analytical solution take similar amounts of time; the
     `mpmath` solution used when `tol=0` is approximately 45 times slower. This
     function takes approximately 8 us normally.

fluids/numerics/__init__.py

@@ -30,7 +30,8 @@
 from fluids.numerics.arrays import (solve as py_solve, inv, dot_product, norm2, matrix_vector_dot, eye,
                      array_as_tridiagonals, tridiagonals_as_array, transpose,
                      solve_tridiagonal, subset_matrix, argsort1d, shape, sort_paired_lists,
-                     stack_vectors, matrix_multiply, transpose, eye, inv, sum_matrix_rows, gelsd)
+                     stack_vectors, matrix_multiply, transpose, eye, inv, sum_matrix_rows, gelsd,
+                     null_space)
 
 from fluids.numerics.special import (py_hypot, py_cacos, py_catan, py_catanh, 
                                      trunc_exp, trunc_log, cbrt, factorial, comb)

fluids/numerics/arrays.py

@@ -48,7 +48,7 @@
            'argsort1d', 'lu', 'gelsd', 'matrix_vector_dot', 'matrix_multiply',
            'sum_matrix_rows', 'sum_matrix_cols', 'sort_paired_lists',
            'scalar_divide_matrix', 'scalar_multiply_matrix', 'scalar_subtract_matrices', 'scalar_add_matrices',
-           'stack_vectors']
+           'stack_vectors', 'null_space']
 primitive_containers = frozenset([list, tuple])
 
 def transpose(matrix):
@@ -1587,4 +1587,52 @@ def gelsd(a, b, rcond=None):
         # Compute residuals as |b - Ax|^2
         diff = [b[i] - Ax[i] for i in range(m)]
         residuals = dot_product(diff, diff)
-    return x, residuals, rank, s
+    return x, residuals, rank, s
+
+def null_space(a, rcond=None):
+    """
+    Construct an orthonormal basis for the null space of A using SVD.
+
+    Parameters
+    ----------
+    a : list[list[float]]
+        Input matrix A of shape (M, N)
+    rcond : float, optional
+        Relative condition number. Singular values ``s`` smaller than
+        ``rcond * max(s)`` are considered zero.
+        Default: floating point eps * max(M,N).
+
+    Returns
+    -------
+    Z : list[list[float]]
+        Orthonormal basis for the null space of A.
+        K = dimension of effective null space, as determined by rcond
+    """
+    import numpy as np
+
+    # Get dimensions and handle empty cases
+    m = len(a)
+    n = len(a[0]) if m > 0 else 0
+
+    if m == 0 or n == 0:
+        return []  # Empty matrix
+
+    # Use numpy only for SVD computation
+    U, s, Vt = np.linalg.svd(np.array(a, dtype=np.float64), full_matrices=True)
+
+    # Convert numpy arrays to Python lists
+    s = s.tolist()
+    Vt = Vt.tolist()
+
+    # Set default rcond
+    if rcond is None:
+        rcond = max(m, n) * 2.2e-16  # Approximate machine epsilon for float64
+
+    # Determine effective null space dimension using rcond
+    tol = max(s) * rcond if s else 0.0
+    num = sum(sv > tol for sv in s)
+    # Extract null space basis
+    V = transpose(Vt)  # V is transpose of Vt
+    Z = [row[num:] for row in V]  # Extract last N - num columns
+
+    return Z

tests/test_numerics_arrays.py

@@ -23,7 +23,7 @@
 
 import pytest
 from fluids.numerics.arrays import (inv, solve, lu, gelsd, eye, dot_product, transpose, matrix_vector_dot, matrix_multiply, sum_matrix_rows, sum_matrix_cols,
-    scalar_divide_matrix, scalar_multiply_matrix, scalar_subtract_matrices, scalar_add_matrices)
+    scalar_divide_matrix, scalar_multiply_matrix, scalar_subtract_matrices, scalar_add_matrices, null_space)
 from fluids.numerics import (
     array_as_tridiagonals,
     assert_close,
@@ -2138,3 +2138,186 @@ def test_sort_paired_lists():
         # Test 6: Unequal length lists
         sort_paired_lists([1, 2], [1])
 
+
+def format_matrix_error_null_space(matrix):
+    """Format a detailed error message for matrix comparison failure"""
+    def matrix_info(matrix):
+        """Get diagnostic information about a matrix"""
+        arr = np.array(matrix)
+        rank = np.linalg.matrix_rank(arr)
+        shape = arr.shape
+        try:
+            cond = np.linalg.cond(arr)
+        except:
+            cond = float('inf')
+        # Only compute determinant for square matrices
+        det = np.linalg.det(arr) if shape[0] == shape[1] else None
+        return {
+            'rank': rank,
+            'condition_number': cond,
+            'shape': shape,
+            'null_space_dim': shape[1] - rank,
+            'determinant': det
+        }
+    info = matrix_info(matrix)
+
+    msg = (
+        f"\nMatrix properties:"
+        f"\n  Shape: {info['shape']}"
+        f"\n  Rank: {info['rank']}"
+        f"\n  Null space dimension: {info['null_space_dim']}"
+        f"\n  Condition number: {info['condition_number']:.2e}"
+    )
+    if info['determinant'] is not None:
+        msg += f"\n  Determinant: {info['determinant']:.2e}"
+    msg += (
+        f"\nInput matrix:"
+        f"\n{np.array2string(np.array(matrix), precision=6, suppress_small=True)}"
+    )
+    return msg
+
+def check_null_space(matrix, rtol=None):
+    """Check if null_space implementation matches scipy behavior"""
+    import scipy
+    just_return = False
+    try:
+        # This will fail for bad matrix inputs
+        cond = np.linalg.cond(np.array(matrix))
+    except:
+        just_return = True
+
+    py_fail = False
+    scipy_fail = False
+
+    try:
+        result = null_space(matrix)
+        if not result:  # Empty result is valid for some cases
+            return
+        result = np.array(result)
+    except:
+        py_fail = True
+
+    # Convert to numpy array if not already
+    matrix = np.array(matrix)
+    try:
+        expected = scipy.linalg.null_space(matrix, rcond=rtol)
+    except:
+        scipy_fail = True
+
+    if py_fail and not scipy_fail:
+        if not just_return and cond > 1e14:
+            # Let ill conditioned matrices pass
+            return 
+        raise ValueError(f"Inconsistent failure states: Python Fail: {py_fail}, SciPy Fail: {scipy_fail}")
+    if py_fail and scipy_fail:
+        return
+    if not py_fail and scipy_fail:
+        return
+    if just_return:
+        return
+
+
+    if rtol is None:
+        rtol = get_rtol(matrix)
+
+    # Compute matrix norm for threshold
+    matrix_norm = np.max(np.sum(np.abs(matrix), axis=1))
+    thresh = matrix_norm * np.finfo(float).eps
+
+
+    # We need to handle sign ambiguity in the basis vectors
+    # Both +v and -v are valid basis vectors
+    # Compare shapes first
+    assert result.shape == expected.shape, \
+           f"Shape mismatch: got {result.shape}, expected {expected.shape}"
+
+    if result.shape[1] > 0:  # Only if we have vectors
+        # For each column in result, check if it matches any expected column or its negative
+        used_expected_cols = set()
+
+        for i in range(result.shape[1]):
+            res_col = result[:, i].reshape(-1, 1)
+            found_match = False
+
+            # Try matching with each unused expected column
+            for j in range(expected.shape[1]):
+                if j in used_expected_cols:
+                    continue
+
+                exp_col = expected[:, j].reshape(-1, 1)
+
+                # Check both orientations with looser tolerance
+                matches_positive = np.allclose(res_col, exp_col, rtol=1e-7, atol=1e-10)
+                matches_negative = np.allclose(res_col, -exp_col, rtol=1e-7, atol=1e-10)
+
+                if matches_positive or matches_negative:
+                    used_expected_cols.add(j)
+                    found_match = True
+                    break
+
+            assert found_match, f"Column {i} doesn't match any expected column in either orientation"
+
+        # Verify orthonormality
+        gram = result.T @ result
+        assert_allclose(gram, np.eye(gram.shape[0]), rtol=1e-10, atol=100*thresh,
+                       err_msg="Basis vectors are not orthonormal")
+
+        # Verify it's actually a null space
+        product = matrix @ result
+        assert_allclose(product, np.zeros_like(product), atol=100*thresh,
+                       err_msg="Result vectors are not in the null space")
+
+# Test cases specific to null space
+matrices_full_rank = [
+    [[1.0]],  # 1x1
+    [[1.0, 0.0], [0.0, 1.0]],  # 2x2 identity
+    [[1.0, 2.0], [3.0, 4.0]],  # 2x2 full rank
+]
+
+matrices_rank_deficient = [
+    [[1.0, 1.0], [1.0, 1.0]],  # 2x2 rank 1
+    [[1.0, 2.0, 3.0], [2.0, 4.0, 6.0]],  # 2x3 rank 1
+    [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 0.0]],  # 3x3 rank 2
+]
+
+matrices_tall = [
+    [[1.0], [0.0], [0.0]],  # 3x1
+    [[1.0, 0.0], [0.0, 1.0], [1.0, 1.0]],  # 3x2
+]
+
+matrices_wide = [
+    [[1.0, 0.0, 0.0]],  # 1x3
+    [[1.0, 0.0, 1.0], [0.0, 1.0, 1.0]],  # 2x3
+]
+
+@pytest.mark.parametrize("matrix", matrices_full_rank)
+def test_null_space_full_rank(matrix):
+    try:
+        check_null_space(matrix)
+    except Exception as e:
+        new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error_null_space(matrix)}"
+        raise Exception(new_message)
+
+@pytest.mark.parametrize("matrix", matrices_rank_deficient)
+def test_null_space_rank_deficient(matrix):
+    try:
+        check_null_space(matrix)
+    except Exception as e:
+        new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error_null_space(matrix)}"
+        raise Exception(new_message)
+
+@pytest.mark.parametrize("matrix", matrices_tall)
+def test_null_space_tall(matrix):
+    try:
+        check_null_space(matrix)
+    except Exception as e:
+        new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error_null_space(matrix)}"
+        raise Exception(new_message)
+
+@pytest.mark.parametrize("matrix", matrices_wide)
+def test_null_space_wide(matrix):
+    try:
+        check_null_space(matrix)
+    except Exception as e:
+        new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error_null_space(matrix)}"
+        raise Exception(new_message)