'''Chemical Engineering Design Library (ChEDL). Utilities for process modeling. Copyright (C) 2024 Caleb Bell Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ''' from math import cos, erf, exp, isnan, log, pi, sin, sqrt 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) from fluids.numerics import ( array_as_tridiagonals, assert_close, assert_close1d, assert_close2d, solve_tridiagonal, subset_matrix, tridiagonals_as_array, argsort1d, sort_paired_lists, ) from fluids.numerics import numpy as np assert_allclose = np.testing.assert_allclose def get_rtol(matrix): """Set tolerance based on condition number""" cond = np.linalg.cond(matrix) machine_eps = np.finfo(float).eps # ≈ 2.2e-16 return min(10 * cond * machine_eps,100*cond * machine_eps if cond > 1e8 else 1e-9) def check_inv(matrix, rtol=None): just_return = False try: # This will fail for bad matrix (inconsistent size) inputs cond = np.linalg.cond(matrix) except: just_return = True py_fail = False numpy_fail = False try: result = inv(matrix) except: py_fail = True try: expected = np.linalg.inv(matrix) except: numpy_fail = True if py_fail and not numpy_fail: if not just_return and cond > 1e14: # Let ill conditioned matrices pass return raise ValueError(f"Inconsistent failure states: Python Fail: {py_fail}, Numpy Fail: {numpy_fail}") if py_fail and numpy_fail: return if not py_fail and numpy_fail: # We'll allow our inv to work with numbers closer to return if just_return: return # Convert result to numpy array if it isn't already result = np.array(result) # Compute infinity norm of input matrix matrix_norm = np.max(np.sum(np.abs(matrix), axis=1)) thresh = matrix_norm * np.finfo(float).eps # We also need to check against the values we get in the inverse; it is helpful # to zero out anything too close to "zero" relative to the values used in the matrix # This is very necessary, and was needed when testing on different CPU architectures inv_norm = np.max(np.sum(np.abs(result), axis=1)) if cond < 1e10: zero_thresh = thresh elif cond < 1e14: zero_thresh = 10*thresh else: zero_thresh = 100*thresh trivial_relative_to_norm_result = (np.abs(result)/inv_norm < zero_thresh) trivial_relative_to_norm_expected = (np.abs(expected)/inv_norm < zero_thresh) # Zero out in both matrices where either condition is met combined_relative_mask = np.logical_or( trivial_relative_to_norm_result, trivial_relative_to_norm_expected ) result[combined_relative_mask] = 0.0 expected[combined_relative_mask] = 0.0 # Check both directions numpy_zeros = (expected == 0.0) our_zeros = (result == 0.0) mask_exact_zeros = numpy_zeros | our_zeros # Where numpy has zeros but we don't; no cases require it but it makes sense to do result[mask_exact_zeros] = np.where(np.abs(result[mask_exact_zeros]) < thresh, 0.0, result[mask_exact_zeros]) # Where we have zeros but numpy doesn't - this is the one we discovered. Apply the check only to `numpy_zeros` expected[numpy_zeros] = np.where(np.abs(expected[numpy_zeros]) < thresh, 0.0, expected[numpy_zeros]) if rtol is None: rtol = get_rtol(matrix) # For each element, use absolute tolerance if the expected value is near zero # In the near zero for some element cases but where others aren't, the relative differences can be brutal relative # to the other numbers in the matrix so we have to treat them differently mask = np.abs(expected) < 1e-14 if mask.any(): assert_allclose(result[mask], expected[mask], atol=thresh) assert_allclose(result[~mask], expected[~mask], rtol=rtol) else: assert_allclose(result, expected, rtol=rtol) def format_matrix_error(matrix): """Format a detailed error message for matrix comparison failure""" def matrix_info(matrix): """Get diagnostic information about a matrix""" arr = np.array(matrix) return { 'condition_number': np.linalg.cond(arr), 'determinant': np.linalg.det(arr), 'shape': arr.shape } info = matrix_info(matrix) return ( f"\nMatrix properties:" f"\n Shape: {info['shape']}" f"\n Condition number: {info['condition_number']:.2e}" f"\n Determinant: {info['determinant']:.2e}" f"\nInput matrix:" f"\n{np.array2string(np.array(matrix), precision=6, suppress_small=True)}" ) # 1x1 matrices matrices_1x1 = [ [[2.0]], [[0.5]], [[-3.0]], [[1e-10]], [[1e10]], ] # 2x2 matrices - regular cases matrices_2x2 = [ [[1.0, 0.0], [0.0, 1.0]], # Identity matrix [[2.0, 1.0], [1.0, 2.0]], # Symmetric matrix [[1.0, 2.0], [3.0, 4.0]], # General case [[1e-5, 1.0], [1.0, 1e5]], # Poorly conditioned [[0.1, 0.2], [0.3, 0.4]], # Decimal values # All ones matrices [[1.0, 1.0], [1.0, 1.0]], [[1.0, 1.0], [1.0, -1.0]], [[1.0, -1.0], [-1.0, 1.0]], # Upper triangular [[2.0, 3.0], [0.0, 4.0]], [[1.0, 10.0], [0.0, 2.0]], [[5.0, -3.0], [0.0, 1.0]], # Lower triangular [[2.0, 0.0], [3.0, 4.0]], [[1.0, 0.0], [10.0, 2.0]], [[5.0, 0.0], [-3.0, 1.0]], # Rotation matrices (θ = 30°, 45°, 60°) [[0.866, -0.5], [0.5, 0.866]], [[0.707, -0.707], [0.707, 0.707]], [[0.5, -0.866], [0.866, 0.5]], # Reflection matrices [[1.0, 0.0], [0.0, -1.0]], [[-1.0, 0.0], [0.0, 1.0]], [[0.0, 1.0], [1.0, 0.0]], # Scaling matrices [[10.0, 0.0], [0.0, 0.1]], [[100.0, 0.0], [0.0, 0.01]], [[1000.0, 0.0], [0.0, 0.001]], # Nearly zero determinant (different from existing near-singular) [[1.0, 2.0], [0.5, 1.0 + 1e-12]], [[2.0, 4.0], [1.0, 2.0 + 1e-13]], [[3.0, 6.0], [1.5, 3.0 + 1e-11]], # Mixed scale [[1e6, 1e-6], [1e-6, 1e6]], [[1e8, 1e-4], [1e-4, 1e8]], [[1e10, 1e-2], [1e-2, 1e10]], # Nilpotent matrices [[0.0, 1.0], [0.0, 0.0]], [[0.0, 2.0], [0.0, 0.0]], [[0.0, 0.5], [0.0, 0.0]], # Hadamard matrices (normalized) [[1/sqrt(2), 1/sqrt(2)], [1/sqrt(2), -1/sqrt(2)]], [[1/sqrt(2), 1/sqrt(2)], [-1/sqrt(2), 1/sqrt(2)]], [[-1/sqrt(2), 1/sqrt(2)], [1/sqrt(2), 1/sqrt(2)]] ] # 2x2 matrices - nearly singular cases matrices_2x2_near_singular = [ [[1.0, 2.0], [1.0 + 1e-10, 2.0 + 1e-10]], # Almost linearly dependent rows [[1.0, 1.0], [1.0, 1.0 + 1e-10]], # Almost zero determinant [[1e5, 1e5], [1e5, 1e5 + 1.0]], # Scaled nearly singular [[1e-10, 1.0], [1.0, 1.0]], # One very small pivot [[1.0, -1e-10], [1e-10, 1.0]], # Almost identity with perturbation # Precision Loss in Subtraction [[1, 1 + 1e-12], [1, 1]], # Case 1 [[1e8, 1e8 + 1], [1e8 + 2, 1e8]], # Case 3 [[1, 1 + 1e-10], [1 + 2e-10, 1]], # Case 4 [[1e20, 1e20 + 10], [1e20 + 20, 1e20]], # Case 5 [[1, 1 + 1e-14], [1 + 1e-14, 1]], # Case 6 # Numerical Instability in Small Matrices [[1e-16, 2e-16], [2e-16, 1e-16]], # Case 1 [[1e-12, 1e-12], [1e-12, 1e-12 + 1e-14]], # Case 2 [[1e-8, 1e-8 + 1e-15], [1e-8 + 1e-15, 1e-8]], # Case 3 [[1, 1 + 1e-13], [1 + 1e-13, 1]], # Case 4 [[1, 1 - 1e-14], [1, 1]], # Case 5 [[1e-15, 1e-15 + 1e-16], [1e-15 + 1e-16, 1e-15]], # Case 6 # # Overflow and Underflow Risks - not a target # [[1e308, 1e-308], [1, 1e-308]], # Case 1 # # [[1e-308, 1e308], [1e-308, 1e-308]], # Case 2 # [[1e308, 1], [1, 1e-308]], # Case 3 # [[1e308, 1e-100], [1e-100, 1e-308]], # Case 4 # [[1e10, 1e-308], [1e-308, 1e10]], # Case 5 # [[1e-308, 1e308], [1e308, 1e-308]], # Case 6 # LU Decomposition Stability [[1e-15, 1], [1, 1]], # Case 1 [[1e-20, 1], [1, 1e-10]], # Case 2 [[1, 1], [1, 1 + 1e-15]], # Case 3 [[1, 1 + 1e-12], [1 + 1e-12, 1]], # Case 4 [[1, 1], [1, 1 + 1e-16]], # Case 5 [[1e-10, 1], [1, 1e-10 + 1e-15]] # Case 6 ] # 3x3 matrices - regular cases matrices_3x3 = [ [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]], # Identity matrix [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]], # Nearly singular [[1.0, 0.5, 0.3], [0.5, 2.0, 0.7], [0.3, 0.7, 3.0]], # Symmetric positive definite [[1e-3, 1.0, 1e3], [1.0, 1.0, 1.0], [1e3, 1.0, 1e-3]], # Poorly conditioned # All ones matrices [[1.0, 1.0, 1.0], [1.0, 1.0, 1.0], [1.0, 1.0, 1.0]], [[1.0, 1.0, -1.0], [1.0, -1.0, 1.0], [-1.0, 1.0, 1.0]], [[1.0, -1.0, 1.0], [-1.0, 1.0, -1.0], [1.0, -1.0, 1.0]], # Upper triangular [[2.0, 3.0, 4.0], [0.0, 5.0, 6.0], [0.0, 0.0, 7.0]], [[1.0, -2.0, 3.0], [0.0, 4.0, -5.0], [0.0, 0.0, 6.0]], [[10.0, 20.0, 30.0], [0.0, 40.0, 50.0], [0.0, 0.0, 60.0]], # Lower triangular [[2.0, 0.0, 0.0], [3.0, 4.0, 0.0], [5.0, 6.0, 7.0]], [[1.0, 0.0, 0.0], [-2.0, 3.0, 0.0], [4.0, -5.0, 6.0]], [[10.0, 0.0, 0.0], [20.0, 30.0, 0.0], [40.0, 50.0, 60.0]], # 3D Rotation matrices (around x, y, and z axes, 45 degrees) [[1.0, 0.0, 0.0], [0.0, 0.707, -0.707], [0.0, 0.707, 0.707]], [[0.707, 0.0, 0.707], [0.0, 1.0, 0.0], [-0.707, 0.0, 0.707]], [[0.707, -0.707, 0.0], [0.707, 0.707, 0.0], [0.0, 0.0, 1.0]], # Permutation matrices [[0.0, 1.0, 0.0], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0]], [[0.0, 0.0, 1.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]], [[1.0, 0.0, 0.0], [0.0, 0.0, 1.0], [0.0, 1.0, 0.0]], # Rank deficient (rank 2) [[1.0, 0.0, 1.0], [0.0, 1.0, 0.0], [2.0, 0.0, 2.0]], [[1.0, 1.0, 2.0], [2.0, 2.0, 4.0], [3.0, 3.0, 6.0]], [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]], # Skew-symmetric matrices [[0.0, 1.0, -2.0], [-1.0, 0.0, 3.0], [2.0, -3.0, 0.0]], [[0.0, 2.0, -1.0], [-2.0, 0.0, 4.0], [1.0, -4.0, 0.0]], [[0.0, 5.0, -3.0], [-5.0, 0.0, 1.0], [3.0, -1.0, 0.0]], # Toeplitz matrices [[1.0, 2.0, 3.0], [2.0, 1.0, 2.0], [3.0, 2.0, 1.0]], [[4.0, -1.0, 2.0], [-1.0, 4.0, -1.0], [2.0, -1.0, 4.0]], [[2.0, 3.0, 4.0], [3.0, 2.0, 3.0], [4.0, 3.0, 2.0]], # Circulant matrices [[1.0, 2.0, 3.0], [3.0, 1.0, 2.0], [2.0, 3.0, 1.0]], [[4.0, 1.0, 2.0], [2.0, 4.0, 1.0], [1.0, 2.0, 4.0]], [[2.0, 3.0, 1.0], [1.0, 2.0, 3.0], [3.0, 1.0, 2.0]], # # Mixed scale with near dependencies [[1e6, 1e-3, 1.0], [1e-3, 1e6, 1.0], [1.0, 1.0, 1e-6]], [[1e9, 1e-6, 1.0], [1e-6, 1e9, 1.0], [1.0, 1.0, 1e-9]], [[1e12, 1e-9, 1.0], [1e-9, 1e12, 1.0], [1.0, 1.0, 1e-12]], # Still challenging but more reasonable condition numbers [[1e3, 1e-2, 1.0], [1e-2, 1e3, 1.0], [1.0, 1.0, 1e-3]], # Condition number ~10^6 [[1e4, 1e-3, 1.0], [1e-3, 1e4, 1.0], [1.0, 1.0, 1e-4]], # Condition number ~10^8 [[1e5, 1e-4, 2.0], [1e-4, 1e5, 2.0], [2.0, 2.0, 1e-5]] # Condition number ~10^10 ] # 3x3 matrices - nearly singular cases matrices_3x3_near_singular = [ [[1.0, 2.0, 3.0], [2.0, 4.0, 6.0], [3.0, 6.0, 9.0 + 1e-10]], # Almost linearly dependent rows [[1e5, 1e5, 1e5], [1e5, 1e5, 1e5], [1e5, 1e5, 1e5 + 1.0]], # Almost zero determinant with scaling [[1.0, 0.0, 1.0], [0.0, 1.0, 1e-10], [0.0, 0.0, 1e-10]], # Nearly dependent columns [[1.0, 0.0, 1e-10], [0.0, 1.0, 1e-10], [1e-10, 1e-10, 1.0]], # Almost rank 2 [[1e8, 1e-8, 1.0], [1e-8, 1e8, 1.0], [1.0, 1.0, 1e-10 + 1.0]], # Scaled with small perturbations # Precision Loss in Subtraction [[1, 1 + 1e-12, 1], [1, 1, 1], [1, 1, 1]], # Case 1 [[1e8, 1e8 + 1, 1e8], [1e8, 1e8, 1e8 + 1e-10], [1e8 + 1e-10, 1e8, 1e8]], # Case 2 [[1e10, 1e10 + 1e-5, 1e10], [1e10, 1e10 + 1e-6, 1e10], [1e10 + 1e-4, 1e10, 1e10]], # Case 3 [[1, 1 + 1e-10, 1], [1, 1 + 1e-11, 1], [1, 1 + 1e-12, 1]], # Case 4 [[1e20, 1e20 + 10, 1e20], [1e20, 1e20, 1e20 + 20], [1e20, 1e20 + 30, 1e20]], # Case 5 [[1, 1 + 1e-14, 1], [1, 1 + 1e-13, 1], [1, 1 + 1e-12, 1]], # Case 6 # Numerical Instability in Small Matrices [[1e-16, 1e-16, 1e-16], [1e-16, 1e-16, 1e-16], [1e-16, 1e-16, 1e-16 + 1e-17]], # Case 1 [[1e-12, 1e-12, 1e-12], [1e-12, 1e-12 + 1e-14, 1e-12], [1e-12, 1e-12, 1e-12]], # Case 2 [[1e-8, 1e-8, 1e-8 + 1e-15], [1e-8, 1e-8 + 1e-15, 1e-8], [1e-8 + 1e-15, 1e-8, 1e-8]], # Case 3 [[1, 1 + 1e-13, 1], [1 + 1e-13, 1, 1], [1, 1 + 1e-13, 1]], # Case 4 [[1, 1 - 1e-14, 1], [1, 1, 1 - 1e-14], [1, 1, 1]], # Case 5 [[1e-15, 1e-15 + 1e-16, 1e-15], [1e-15, 1e-15, 1e-15 + 1e-16], [1e-15, 1e-15, 1e-15]], # Case 6 # # Overflow and Underflow Risks # [[1e308, 1e-308, 1e308], [1, 1e-308, 1], [1e308, 1, 1e-308]], # Case 1 # [[1e-308, 1e308, 1e-308], [1e-308, 1e-308, 1e308], [1e308, 1e-308, 1e-308]], # Case 2 # [[1e308, 1e-100, 1], [1, 1e308, 1e-308], [1e-308, 1, 1e308]], # Case 3 # [[1e308, 1e-308, 1], [1, 1e308, 1e-308], [1e-308, 1, 1e308]], # Case 4 # [[1e10, 1e-308, 1], [1e-308, 1e10, 1e-308], [1, 1e-308, 1e10]], # Case 5 # [[1e-308, 1e308, 1], [1e308, 1e-308, 1], [1, 1, 1e308]], # Case 6 # LU Decomposition Stability [[1e-15, 1, 1], [1, 1, 1], [1, 1, 1]], # Case 1 [[1e-20, 1, 1], [1, 1e-10, 1], [1, 1, 1e-10]], # Case 2 [[1, 1, 1], [1, 1, 1 + 1e-15], [1, 1, 1]], # Case 3 [[1, 1 + 1e-12, 1], [1 + 1e-12, 1, 1], [1, 1, 1]], # Case 4 [[1, 1, 1], [1, 1, 1 + 1e-16], [1, 1, 1]], # Case 5 [[1e-10, 1, 1], [1, 1e-10, 1], [1, 1, 1e-10 + 1e-15]], # Case 6 ] # 4x4 matrices - regular cases matrices_4x4 = [ [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]], # Identity matrix [[1.0, 2.0, 3.0, 4.0], [5.0, 6.0, 7.0, 8.0], [9.0, 10.0, 11.0, 12.0], [13.0, 14.0, 15.0, 16.0]], # Singular [[1.0, 0.1, 0.1, 0.1], [0.1, 2.0, 0.2, 0.2], [0.1, 0.2, 3.0, 0.3], [0.1, 0.2, 0.3, 4.0]], # Diagonally dominant # All ones matrices [[1.0, 1.0, 1.0, 1.0], [1.0, 1.0, 1.0, 1.0], [1.0, 1.0, 1.0, 1.0], [1.0, 1.0, 1.0, 1.0]], [[1.0, 1.0, -1.0, -1.0], [1.0, -1.0, 1.0, -1.0], [-1.0, 1.0, 1.0, -1.0], [-1.0, -1.0, -1.0, 1.0]], [[1.0, -1.0, 1.0, -1.0], [-1.0, 1.0, -1.0, 1.0], [1.0, -1.0, 1.0, -1.0], [-1.0, 1.0, -1.0, 1.0]], # Upper triangular [[1.0, 2.0, 3.0, 4.0], [0.0, 5.0, 6.0, 7.0], [0.0, 0.0, 8.0, 9.0], [0.0, 0.0, 0.0, 10.0]], [[2.0, -1.0, 3.0, -2.0], [0.0, 4.0, -5.0, 6.0], [0.0, 0.0, 7.0, -8.0], [0.0, 0.0, 0.0, 9.0]], [[10.0, 20.0, 30.0, 40.0], [0.0, 50.0, 60.0, 70.0], [0.0, 0.0, 80.0, 90.0], [0.0, 0.0, 0.0, 100.0]], # Lower triangular [[1.0, 0.0, 0.0, 0.0], [2.0, 3.0, 0.0, 0.0], [4.0, 5.0, 6.0, 0.0], [7.0, 8.0, 9.0, 10.0]], [[2.0, 0.0, 0.0, 0.0], [-1.0, 3.0, 0.0, 0.0], [4.0, -5.0, 6.0, 0.0], [-7.0, 8.0, -9.0, 10.0]], [[10.0, 0.0, 0.0, 0.0], [20.0, 30.0, 0.0, 0.0], [40.0, 50.0, 60.0, 0.0], [70.0, 80.0, 90.0, 100.0]], # Block diagonal (2x2 blocks) [[2.0, 1.0, 0.0, 0.0], [1.0, 2.0, 0.0, 0.0], [0.0, 0.0, 3.0, 1.0], [0.0, 0.0, 1.0, 3.0]], [[4.0, -1.0, 0.0, 0.0], [-1.0, 4.0, 0.0, 0.0], [0.0, 0.0, 5.0, -1.0], [0.0, 0.0, -1.0, 5.0]], [[1.0, 0.5, 0.0, 0.0], [0.5, 1.0, 0.0, 0.0], [0.0, 0.0, 2.0, 0.5], [0.0, 0.0, 0.5, 2.0]], # Block tridiagonal [[2.0, 1.0, 0.0, 0.0], [1.0, 2.0, 1.0, 0.0], [0.0, 1.0, 2.0, 1.0], [0.0, 0.0, 1.0, 2.0]], [[4.0, -1.0, 0.0, 0.0], [-1.0, 4.0, -1.0, 0.0], [0.0, -1.0, 4.0, -1.0], [0.0, 0.0, -1.0, 4.0]], [[3.0, 1.0, 0.0, 0.0], [1.0, 3.0, 1.0, 0.0], [0.0, 1.0, 3.0, 1.0], [0.0, 0.0, 1.0, 3.0]], # Sparse matrix patterns [[1.0, 0.0, 2.0, 0.0], [0.0, 3.0, 0.0, 4.0], [2.0, 0.0, 5.0, 0.0], [0.0, 4.0, 0.0, 6.0]], [[2.0, 0.0, 0.0, 1.0], [0.0, 3.0, 1.0, 0.0], [0.0, 1.0, 4.0, 0.0], [1.0, 0.0, 0.0, 5.0]], [[1.0, 1.0, 0.0, 0.0], [1.0, 2.0, 1.0, 0.0], [0.0, 1.0, 3.0, 1.0], [0.0, 0.0, 1.0, 4.0]], # Vandermonde matrices # [[1.0, 1.0, 1.0, 1.0], # failing on other CPUs in test_lu_4x4 # [1.0, 2.0, 4.0, 8.0], # [1.0, 3.0, 9.0, 27.0], # [1.0, 4.0, 16.0, 64.0]], [[1.0, 1.0, 1.0, 1.0], [1.0, -1.0, 1.0, -1.0], [1.0, -2.0, 4.0, -8.0], [1.0, -3.0, 9.0, -27.0]], [[1.0, 1.0, 1.0, 1.0], [1.0, 0.5, 0.25, 0.125], [1.0, 2.0, 4.0, 8.0], [1.0, 3.0, 9.0, 27.0]], # Hilbert matrix segments [[1.0, 1/2, 1/3, 1/4], [1/2, 1/3, 1/4, 1/5], [1/3, 1/4, 1/5, 1/6], [1/4, 1/5, 1/6, 1/7]], [[1.0, 1/3, 1/5, 1/7], [1/3, 1/5, 1/7, 1/9], [1/5, 1/7, 1/9, 1/11], [1/7, 1/9, 1/11, 1/13]], [[1/2, 1/3, 1/4, 1/5], [1/3, 1/4, 1/5, 1/6], [1/4, 1/5, 1/6, 1/7], [1/5, 1/6, 1/7, 1/8]], # Rank deficient (rank 3) [[1.0, 0.0, 0.0, 1.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [2.0, 0.0, 0.0, 2.0]], [[1.0, 1.0, 1.0, 3.0], [2.0, 2.0, 2.0, 6.0], [3.0, 3.0, 3.0, 9.0], [4.0, 4.0, 4.0, 12.0]], [[1.0, 2.0, 3.0, 6.0], [4.0, 5.0, 6.0, 15.0], [7.0, 8.0, 9.0, 24.0], [10.0, 11.0, 12.0, 33.0]], # Mixed scale with multiple near dependencies [[1e6, 1e-3, 1.0, 1e-6], [1e-3, 1e6, 1e-6, 1.0], [1.0, 1e-6, 1e6, 1e-3], [1e-6, 1.0, 1e-3, 1e6]], [[1e9, 1e-6, 1.0, 1e-9], [1e-6, 1e9, 1e-9, 1.0], [1.0, 1e-9, 1e9, 1e-6], [1e-9, 1.0, 1e-6, 1e9]], [[1e12, 1e-9, 1.0, 1e-12], [1e-9, 1e12, 1e-12, 1.0], [1.0, 1e-12, 1e12, 1e-9], [1e-12, 1.0, 1e-9, 1e12]] ] # 4x4 matrices - nearly singular cases matrices_4x4_near_singular = [ [[1.0, 2.0, 3.0, 4.0], [2.0, 4.0, 6.0, 8.0], [3.0, 6.0, 9.0, 12.0], [4.0, 8.0, 12.0, 16.0 + 1e-10]], # Almost linearly dependent rows [[1e3, 1e3, 1e3, 1e3], [1e3, 1e3, 1e3, 1e3], [1e3, 1e3, 1e3, 1e3], [1e3, 1e3, 1e3, 1e3 + 1.0]], # Almost zero determinant with scaling # [[1.0, 0.0, 0.0, 1.0], # [0.0, 1.0, 0.0, 1e-10], # [0.0, 0.0, 1.0, 1e-10], # [0.0, 0.0, 0.0, 1e-10]], # Nearly dependent columns, too hard on some CPUs # [[1e5, 1e-5, 1.0, 1.0], # [1e-5, 1e5, 1.0, 1.0], # [1.0, 1.0, 1e-10, 1.0], # [1.0, 1.0, 1.0, 1e-10 + 1.0]], # Mixed scaling with near dependencies # Precision Loss in Subtraction [[1, 1 + 1e-12, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1]], # Case 1 [[1e8, 1e8 + 1, 1e8, 1e8], [1e8, 1e8, 1e8 + 1e-10, 1e8], [1e8, 1e8, 1e8, 1e8], [1e8 + 1e-10, 1e8, 1e8, 1e8]], # Case 2 [[1e10, 1e10 + 1e-5, 1e10, 1e10], [1e10, 1e10, 1e10 + 1e-6, 1e10], [1e10 + 1e-4, 1e10, 1e10, 1e10], [1e10, 1e10, 1e10, 1e10 + 1e-3]], # Case 3 # [[1, 1 + 1e-10, 1, 1], [1, 1 + 1e-11, 1, 1], [1, 1, 1 + 1e-12, 1], [1, 1, 1, 1 + 1e-13]], # Case 4 [[1e20, 1e20 + 10, 1e20, 1e20], [1e20, 1e20, 1e20 + 20, 1e20], [1e20, 1e20 + 30, 1e20, 1e20], [1e20, 1e20, 1e20 + 40, 1e20]], # Case 5 [[1, 1 + 1e-14, 1, 1], [1, 1 + 1e-13, 1, 1], [1, 1, 1 + 1e-12, 1], [1, 1, 1, 1 + 1e-11]], # Case 6 # Numerical Instability in Small Matrices [[1e-16, 1e-16, 1e-16, 1e-16], [1e-16, 1e-16, 1e-16, 1e-16], [1e-16, 1e-16, 1e-16, 1e-16], [1e-16, 1e-16, 1e-16, 1e-16 + 1e-17]], # Case 1 [[1e-12, 1e-12, 1e-12, 1e-12], [1e-12, 1e-12 + 1e-14, 1e-12, 1e-12], [1e-12, 1e-12, 1e-12, 1e-12], [1e-12, 1e-12, 1e-12, 1e-12]], # Case 2 # [[1e-8, 1e-8, 1e-8 + 1e-15, 1e-8], [1e-8, 1e-8, 1e-8, 1e-8 + 1e-15], [1e-8, 1e-8, 1e-8, 1e-8], [1e-8 + 1e-15, 1e-8, 1e-8, 1e-8]], # Case 3 [[1, 1 + 1e-13, 1, 1], [1, 1, 1 + 1e-13, 1], [1, 1, 1, 1 + 1e-13], [1 + 1e-13, 1, 1, 1]], # Case 4 [[1, 1 - 1e-14, 1, 1], [1, 1, 1 - 1e-14, 1], [1, 1, 1, 1 - 1e-14], [1, 1, 1, 1]], # Case 5 [[1e-15, 1e-15 + 1e-16, 1e-15, 1e-15], [1e-15, 1e-15, 1e-15 + 1e-16, 1e-15], [1e-15, 1e-15, 1e-15, 1e-15], [1e-15, 1e-15, 1e-15, 1e-15]], # Case 6 # # Overflow and Underflow Risks - not a target # [[1e308, 1e-308, 1e308, 1e-308], [1e-308, 1e308, 1, 1e-308], [1e308, 1, 1e-308, 1], [1, 1e-308, 1e308, 1e-308]], # Case 1 # [[1e-308, 1e308, 1e-308, 1], [1e308, 1e-308, 1, 1e308], [1, 1e-308, 1e308, 1e-308], [1e-308, 1, 1, 1e308]], # Case 2 # [[1e308, 1e-100, 1, 1], [1, 1e308, 1e-308, 1], [1e-308, 1, 1e308, 1], [1, 1, 1e-308, 1e308]], # Case 3 # [[1e308, 1e-308, 1, 1], [1, 1e308, 1e-308, 1], [1e-308, 1, 1e308, 1], [1, 1e-308, 1, 1e308]], # Case 4 # [[1e10, 1e-308, 1, 1e-308], [1e-308, 1e10, 1e-308, 1], [1, 1e-308, 1e10, 1e-308], [1e-308, 1, 1e-308, 1e10]], # Case 5 # [[1e-308, 1e308, 1, 1], [1e308, 1e-308, 1, 1e308], [1, 1, 1e308, 1e-308], [1e-308, 1e308, 1, 1]], # Case 6 # LU Decomposition Stability [[1e-15, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1]], # Case 1 [[1e-20, 1, 1, 1], [1, 1e-10, 1, 1], [1, 1, 1, 1e-10], [1, 1, 1, 1]], # Case 2 [[1, 1, 1, 1], [1, 1, 1, 1 + 1e-15], [1, 1, 1, 1], [1, 1, 1, 1]], # Case 3 [[1, 1 + 1e-12, 1, 1], [1 + 1e-12, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1]], # Case 4 [[1, 1, 1, 1], [1, 1, 1, 1 + 1e-16], [1, 1, 1, 1], [1, 1, 1, 1]], # Case 5 [[1e-10, 1, 1, 1], [1, 1e-10, 1, 1], [1, 1, 1e-10, 1], [1, 1, 1, 1e-10 + 1e-15]], # Case 6 ] # Singular matrices that should raise exceptions matrices_singular = [ [[0.0]], # Singular 1x1 [[1.0, 0.0], [0.0, 0.0]], # Singular 2x2 [[1.0, 0.0], [0.0]] # Irregular matrix ] @pytest.mark.parametrize("matrix", matrices_1x1) def test_inv_1x1(matrix): try: check_inv(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) @pytest.mark.parametrize("matrix", matrices_2x2) def test_inv_2x2(matrix): try: check_inv(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) @pytest.mark.parametrize("matrix", matrices_2x2_near_singular) def test_inv_2x2_near_singular(matrix): try: check_inv(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) @pytest.mark.parametrize("matrix", matrices_3x3) def test_inv_3x3(matrix): try: check_inv(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) @pytest.mark.parametrize("matrix", matrices_3x3_near_singular) def test_inv_3x3_near_singular(matrix): try: check_inv(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) @pytest.mark.parametrize("matrix", matrices_4x4) def test_inv_4x4(matrix): try: check_inv(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) @pytest.mark.parametrize("matrix", matrices_4x4_near_singular) def test_inv_4x4_near_singular(matrix): try: check_inv(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) @pytest.mark.parametrize("matrix", matrices_singular) def test_inv_singular_matrices(matrix): try: check_inv(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) def check_solve(matrix, b=None): """Set tolerance based on condition number and check solution""" if b is None: # Create a right-hand side vector that's compatible with the matrix size b = [1.0] * len(matrix) just_return = False try: # This will fail for bad matrix (inconsistent size) inputs cond = np.linalg.cond(matrix) except: just_return = True py_fail = False numpy_fail = False try: result = solve(matrix, b) except: py_fail = True try: expected = np.linalg.solve(matrix, b) except: numpy_fail = True if py_fail and not numpy_fail: if not just_return and cond > 1e14: # Let ill conditioned matrices pass return raise ValueError(f"Inconsistent failure states: Python Fail: {py_fail}, Numpy Fail: {numpy_fail}") if py_fail and numpy_fail: return if not py_fail and numpy_fail: return if just_return: return # Convert result to numpy array if it isn't already result = np.array(result) expected = np.array(expected) # Compute infinity norm of input matrix matrix_norm = np.max(np.sum(np.abs(matrix), axis=1)) thresh = matrix_norm * np.finfo(float).eps # Get solution norms sol_norm = np.max(np.abs(result)) # Adjust tolerance based on condition number if cond < 1e10: zero_thresh = thresh rtol = 10 * cond * np.finfo(float).eps elif cond < 1e14: zero_thresh = 10*thresh rtol = 10 * cond * np.finfo(float).eps else: zero_thresh = 100*thresh rtol = 100 * cond * np.finfo(float).eps # Zero out small values relative to solution norm trivial_relative_to_norm_result = (np.abs(result)/sol_norm < zero_thresh) trivial_relative_to_norm_expected = (np.abs(expected)/sol_norm < zero_thresh) # Zero out in both solutions where either condition is met combined_relative_mask = np.logical_or( trivial_relative_to_norm_result, trivial_relative_to_norm_expected ) result[combined_relative_mask] = 0.0 expected[trivial_relative_to_norm_expected] = 0.0 np.testing.assert_allclose(result, expected, rtol=rtol) @pytest.mark.parametrize("matrix", matrices_1x1) def test_solve_1x1(matrix): try: check_solve(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) @pytest.mark.parametrize("matrix", matrices_2x2) def test_solve_2x2(matrix): try: check_solve(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) @pytest.mark.parametrize("matrix", matrices_2x2_near_singular) def test_solve_2x2_near_singular(matrix): try: check_solve(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) @pytest.mark.parametrize("matrix", matrices_3x3) def test_solve_3x3(matrix): try: check_solve(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) @pytest.mark.parametrize("matrix", matrices_3x3_near_singular) def test_solve_3x3_near_singular(matrix): try: check_solve(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) @pytest.mark.parametrize("matrix", matrices_4x4) def test_solve_4x4(matrix): try: check_solve(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) @pytest.mark.parametrize("matrix", matrices_4x4_near_singular) def test_solve_4x4_near_singular(matrix): try: check_solve(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) @pytest.mark.parametrize("matrix", matrices_singular) def test_solve_singular_matrices(matrix): try: check_solve(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) specific_rhs_cases = [ ([[2.0, 1.0], [1.0, 3.0]], [1.0, 1.0]), ([[2.0, 1.0], [1.0, 3.0]], [1.0, -1.0]), ([[2.0, 1.0], [1.0, 3.0]], [1e6, 1e-6]), ([[2.0, 1.0], [1.0, 3.0]], [0.0, 1.0]), ([[2.0, 1.0], [1.0, 3.0]], [1e-15, 1e-15]), ([[1.0, 0.0], [0.0, 1.0]], [-1.0, 1.0]), ([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]], [1.0, 2.0, 3.0]), ([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]], [1e3, 1e0, 1e-3]), ([[1e5, 0.0], [0.0, 1e-5]], [1.0, 1.0]), ([[1.0, 1.0], [1.0, 1.0 + 1e-10]], [1.0, 1.0]), ([[2.0, 0.0, 0.0, 0.0], [0.0, 2.0, 0.0, 0.0], [0.0, 0.0, 2.0, 0.0], [0.0, 0.0, 0.0, 2.0]], [1.0, -1.0, 1.0, -1.0]), ([[1.0, 0.1, 0.1, 0.1], [0.1, 1.0, 0.1, 0.1], [0.1, 0.1, 1.0, 0.1], [0.1, 0.1, 0.1, 1.0]], [1e-8, 1e-8, 1e8, 1e8]) ] @pytest.mark.parametrize("matrix,rhs", specific_rhs_cases) def test_solve_specific_rhs(matrix, rhs): try: check_solve(matrix, rhs) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}\nRHS: {rhs}" raise Exception(new_message) def test_py_solve_bad_cases(): j = [[-3.8789618086360855, -3.8439678951838587, -1.1398039850146757e-07], [1.878915113936518, 1.8439217680605073, 1.139794740950828e-07], [-1.0, -1.0, 0.0]] nv = [-1.4181331207951953e-07, 1.418121622354107e-07, 2.220446049250313e-16] import fluids.numerics calc = fluids.numerics.py_solve(j, nv) import numpy as np expect = np.linalg.solve(j, nv) fluids.numerics.assert_close1d(calc, expect, rtol=1e-4) specific_solution_cases = [ # Case 1 ([ [0.8660254037844387, -0.49999999999999994, 0.0], [0.49999999999999994, 0.8660254037844387, 0.0], [0.0, 0.0, 1.0]], [1, 2, 3], [1.8660254037844384, 1.2320508075688774, 3.0]), # Case 2 ([ [4, -1, 0, 0], [-1, 4, -1, 0], [0, -1, 4, -1], [0, 0, -1, 4]], [1, -1, 1, -1], [0.21052631578947367, -0.15789473684210525, 0.15789473684210528, -0.21052631578947367]), # Case 3 ([ [2, 1, 1], [0, 3, -1], [0, 0, 4]], [1, 2, 3], [-0.3333333333333333, 0.9166666666666666, 0.75]), # Case 4 ([ [3, 1, -2], [2, -3, 1], [-1, 2, 4]], [7, -1, 3], [1.981132075471698, 1.7735849056603772, 0.3584905660377357]), # Case 5 ([[3.0]], [6.0], [2.0]), # Case 6 ([ [0.7071067811865476, -0.7071067811865475], [0.7071067811865475, 0.7071067811865476]], [1.0, 2.0], [2.1213203435596424, 0.7071067811865478]), # Case 7 ([ [2, 1, 1], [0, 3, -1], [0, 0, 4]], [1, 2, 3], [-0.3333333333333333, 0.9166666666666666, 0.75]), # Case 8 ([ [4, -1, 0, 0, 0], [-1, 4, -1, 0, 0], [0, -1, 4, -1, 0], [0, 0, -1, 4, -1], [0, 0, 0, -1, 4]], [1, -1, 1, -1, 1], [0.21153846153846154, -0.15384615384615383, 0.17307692307692307, -0.15384615384615383, 0.21153846153846154]), # Case 9 ([ [2.0, 1.0, 0.0, 0.0, 0.0, 0.0], [1.0, 2.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 3.0, -1.0, 0.0, 0.0], [0.0, 0.0, -1.0, 3.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 4.0, -1.0], [0.0, 0.0, 0.0, 0.0, -1.0, 4.0]], [1, -1, 2, -2, 3, -3], [1.0, -1.0, 0.5, -0.5000000000000001, 0.6, -0.6]), # Case 10 ([ [2, 1, 0, 0, 0, 0, 0], [-1, 3, -1, 0, 0, 0, 0], [0, 1, 2, 1, 0, 0, 0], [0, 0, -1, 4, -1, 0, 0], [0, 0, 0, 1, 3, 1, 0], [0, 0, 0, 0, -1, 2, -1], [0, 0, 0, 0, 0, 1, 3]], [1, -1, 1, -1, 1, -1, 1], [0.4983480176211454, 0.0033039647577092373, 0.5115638766519823, -0.02643171806167401, 0.3827092511013216, -0.12169603524229072, 0.37389867841409696]), ] @pytest.mark.parametrize("matrix,rhs,expected", specific_solution_cases) def test_solve_specific_solutions(matrix, rhs, expected): result = solve(matrix, rhs) assert_allclose(result, expected, rtol=1e-15) def check_lu(matrix): """Compare our LU decomposition against SciPy's""" import numpy as np from scipy import linalg just_return = False try: # This will fail for bad matrix (inconsistent size) inputs cond = np.linalg.cond(matrix) except: just_return = True py_fail = False scipy_fail = False try: P, L, U = lu(matrix) except: py_fail = True try: p, l, u = linalg.lu(matrix) 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 # Convert results to numpy arrays P, L, U = np.array(P), np.array(L), np.array(U) # Compute infinity norm of input matrix matrix_norm = np.max(np.sum(np.abs(matrix), axis=1)) thresh = matrix_norm * np.finfo(float).eps # Verify L is lower triangular with unit diagonal tril_mask = np.tril(np.ones_like(L, dtype=bool)) assert_allclose(L[~tril_mask], 0, atol=thresh) assert_allclose(np.diag(L), 1, rtol=thresh) # Verify U is upper triangular triu_mask = np.triu(np.ones_like(U, dtype=bool)) assert_allclose(U[~triu_mask], 0, atol=thresh) # Check that P is a permutation matrix P_sum_rows = np.sum(P, axis=1) P_sum_cols = np.sum(P, axis=0) assert_allclose(P_sum_rows, np.ones(len(matrix)), rtol=thresh) assert_allclose(P_sum_cols, np.ones(len(matrix)), rtol=thresh) # Most importantly: verify that PA = LU PA = P @ matrix LU = L @ U assert_allclose(PA, LU, rtol=1e-13, atol=10*thresh) # Compare with SciPy's results: # Since pivot choices might differ, we compare # The upper triangular factor (which should be unique up to sign changes) # Normalize each row to handle sign differences U_normalized = U / (np.max(np.abs(U), axis=1, keepdims=True) + np.finfo(float).eps) u_normalized = u / (np.max(np.abs(u), axis=1, keepdims=True) + np.finfo(float).eps) if cond < 1e7: np.testing.assert_allclose(np.abs(U_normalized), np.abs(u_normalized), rtol=1e-13) specific_lu_cases = [ # Case 1 ( [[0.8660254037844387, -0.49999999999999994, 0.0], [0.49999999999999994, 0.8660254037844387, 0.0], [0.0, 0.0, 1.0]], [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]], [[1.0, 0.0, 0.0], [0.5773502691896256, 1.0, 0.0], [0.0, 0.0, 1.0]], [[0.8660254037844387, -0.49999999999999994, 0.0], [0.0, 1.1547005383792515, 0.0], [0.0, 0.0, 1.0]] ), # Case 2 ( [[4.0, -1.0, 0.0, 0.0], [-1.0, 4.0, -1.0, 0.0], [0.0, -1.0, 4.0, -1.0], [0.0, 0.0, -1.0, 4.0]], [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]], [[1.0, 0.0, 0.0, 0.0], [-0.25, 1.0, 0.0, 0.0], [0.0, -0.26666666666666666, 1.0, 0.0], [0.0, 0.0, -0.26785714285714285, 1.0]], [[4.0, -1.0, 0.0, 0.0], [0.0, 3.75, -1.0, 0.0], [0.0, 0.0, 3.7333333333333334, -1.0], [0.0, 0.0, 0.0, 3.732142857142857]] ), # Case 3 ( [[2.0, 1.0, 1.0], [0.0, 3.0, -1.0], [0.0, 0.0, 4.0]], [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]], [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]], [[2.0, 1.0, 1.0], [0.0, 3.0, -1.0], [0.0, 0.0, 4.0]] ), # Case 4 ( [[3.0, 1.0, -2.0], [2.0, -3.0, 1.0], [-1.0, 2.0, 4.0]], [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]], [[1.0, 0.0, 0.0], [0.6666666666666666, 1.0, 0.0], [-0.3333333333333333, -0.6363636363636365, 1.0]], [[3.0, 1.0, -2.0], [0.0, -3.6666666666666665, 2.333333333333333], [0.0, 0.0, 4.818181818181818]] ), # Case 5 ( [[3.0]], [[1.0]], [[1.0]], [[3.0]] ), # Case 6 ( [[0.7071067811865476, -0.7071067811865475], [0.7071067811865475, 0.7071067811865476]], [[1.0, 0.0], [0.0, 1.0]], [[1.0, 0.0], [0.9999999999999998, 1.0]], [[0.7071067811865476, -0.7071067811865475], [0.0, 1.414213562373095]] ), # Case 7 ( [[4.0, -1.0, 0.0, 0.0, 0.0], [-1.0, 4.0, -1.0, 0.0, 0.0], [0.0, -1.0, 4.0, -1.0, 0.0], [0.0, 0.0, -1.0, 4.0, -1.0], [0.0, 0.0, 0.0, -1.0, 4.0]], [[1.0, 0.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 0.0, 1.0]], [[1.0, 0.0, 0.0, 0.0, 0.0], [-0.25, 1.0, 0.0, 0.0, 0.0], [0.0, -0.26666666666666666, 1.0, 0.0, 0.0], [0.0, 0.0, -0.26785714285714285, 1.0, 0.0], [0.0, 0.0, 0.0, -0.2679425837320574, 1.0]], [[4.0, -1.0, 0.0, 0.0, 0.0], [0.0, 3.75, -1.0, 0.0, 0.0], [0.0, 0.0, 3.7333333333333334, -1.0, 0.0], [0.0, 0.0, 0.0, 3.732142857142857, -1.0], [0.0, 0.0, 0.0, 0.0, 3.7320574162679425]] ), # Case 8 ( [[2.0, 1.0, 0.0, 0.0, 0.0, 0.0], [1.0, 2.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 3.0, -1.0, 0.0, 0.0], [0.0, 0.0, -1.0, 3.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 4.0, -1.0], [0.0, 0.0, 0.0, 0.0, -1.0, 4.0]], [[1.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 1.0]], [[1.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.5, 1.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0, 0.0, 0.0], [0.0, 0.0, -0.3333333333333333, 1.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 0.0, -0.25, 1.0]], [[2.0, 1.0, 0.0, 0.0, 0.0, 0.0], [0.0, 1.5, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 3.0, -1.0, 0.0, 0.0], [0.0, 0.0, 0.0, 2.6666666666666665, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 4.0, -1.0], [0.0, 0.0, 0.0, 0.0, 0.0, 3.75]] ), ] @pytest.mark.parametrize("matrix,p_expected,l_expected,u_expected", specific_lu_cases) def test_lu_specific_cases(matrix, p_expected, l_expected, u_expected): p, l, u = lu(matrix) assert_allclose(p, p_expected, rtol=1e-15) assert_allclose(l, l_expected, rtol=1e-15) assert_allclose(u, u_expected, rtol=1e-15) @pytest.mark.parametrize("matrix", matrices_1x1) def test_lu_1x1(matrix): try: check_lu(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) @pytest.mark.parametrize("matrix", matrices_2x2) def test_lu_2x2(matrix): try: check_lu(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) @pytest.mark.parametrize("matrix", matrices_2x2_near_singular) def test_lu_2x2_near_singular(matrix): try: check_lu(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) @pytest.mark.parametrize("matrix", matrices_3x3) def test_lu_3x3(matrix): try: check_lu(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) @pytest.mark.parametrize("matrix", matrices_3x3_near_singular) def test_lu_3x3_near_singular(matrix): try: check_lu(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) @pytest.mark.parametrize("matrix", matrices_4x4) def test_lu_4x4(matrix): try: check_lu(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) @pytest.mark.parametrize("matrix", matrices_4x4_near_singular) def test_lu_4x4_near_singular(matrix): try: check_lu(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) @pytest.mark.parametrize("matrix", matrices_singular) def test_lu_singular_matrices(matrix): try: check_lu(matrix) except Exception as e: new_message = f"Original error: {str(e)}\nAdditional context: {format_matrix_error(matrix)}" raise Exception(new_message) def test_gelsd_basic(): """Test basic functionality with simple well-conditioned problems""" # Simple 2x2 system A = [[1.0, 2.0], [3.0, 4.0]] b = [5.0, 6.0] x, residuals, rank, s = gelsd(A, b) # Compare with numpy's lstsq x_numpy = np.linalg.lstsq(np.array(A), np.array(b), rcond=None)[0] assert_allclose(x, x_numpy, rtol=1e-14) assert rank == 2 assert len(s) == 2 def test_gelsd_overdetermined(): """Test overdetermined system (more equations than unknowns)""" A = [[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]] b = [7.0, 8.0, 9.0] x, residuals, rank, s = gelsd(A, b) # Verify dimensions assert len(x) == 2 assert rank == 2 assert len(s) == 2 # Check residuals are positive for overdetermined system assert residuals >= 0 def test_gelsd_underdetermined(): """Test underdetermined system (fewer equations than unknowns)""" A = [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]] b = [7.0, 8.0] x, residuals, rank, s = gelsd(A, b) # Verify dimensions assert len(x) == 3 assert rank == 2 assert len(s) == 2 assert residuals == 0.0 # Should be exactly solvable def test_gelsd_ill_conditioned(): """Test behavior with ill-conditioned matrix""" A = [[1.0, 1.0], [1.0, 1.0 + 1e-15]] b = [2.0, 2.0] x, residuals, rank, s = gelsd(A, b) # Matrix should be detected as rank deficient assert rank == 1 assert s[0]/s[1] > 1e14 # Check condition number def test_gelsd_zero_matrix(): """Test with zero matrix""" A = [[0.0, 0.0], [0.0, 0.0]] b = [1.0, 1.0] x, residuals, rank, s = gelsd(A, b) assert rank == 0 assert all(sv == 0 for sv in s) def test_gelsd_against_lapack(): """Compare results against LAPACK's dgelsd""" from scipy.linalg import lapack A = [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0], [10.0, 11.0, 12.0]] b = [13.0, 14.0, 15.0, 16.0] # Our implementation x1, residuals1, rank1, s1 = gelsd(A, b) # LAPACK implementation m, n = np.array(A).shape minmn = min(m, n) maxmn = max(m, n) x2, s2, rank2, info = lapack.dgelsd(A, b, lwork=10000, size_iwork=10000) x2 = x2.ravel() # Compare results assert_allclose(x1, x2[:n], rtol=1e-12, atol=1e-12) assert_allclose(s1, s2[:minmn], rtol=1e-12, atol=1e-12) assert rank1 == rank2 @pytest.mark.parametrize("A, b, name", [ # Standard square matrix ([[1.0, 2.0], [3.0, 4.0]], [5.0, 6.0], "2x2 well-conditioned"), # Original test case ([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0], [10.0, 11.0, 12.0]], [13.0, 14.0, 15.0, 16.0], "4x3 overdetermined"), # Overdetermined system ([[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]], [7.0, 8.0, 9.0], "3x2 overdetermined"), # # Nearly singular system # ([[1.0, 1.0], # [1.0, 1.0 + 1e-6]], # 1e-10 broke on some CPUs 1e-6 didn't help # [2.0, 2.0], # "2x2 nearly singular"), # Zero matrix ([[0.0, 0.0], [0.0, 0.0]], [1.0, 1.0], "2x2 zero matrix"), # Underdetermined system ([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]], [7.0, 8.0], "2x3 underdetermined"), # Ill-conditioned matrix ([[1e-10, 1.0], [1.0, 1.0]], [1.0, 2.0], "2x2 ill-conditioned"), # Larger system ([[1.0, 2.0, 3.0, 4.0], [5.0, 6.0, 7.0, 8.0], [9.0, 10.0, 11.0, 12.0], [13.0, 14.0, 15.0, 16.0], [17.0, 18.0, 19.0, 20.0]], [21.0, 22.0, 23.0, 24.0, 25.0], "5x4 larger system") ]) def test_gelsd_against_lapack2(A, b, name): """Compare GELSD results against LAPACK's dgelsd for various test cases""" from scipy.linalg import lapack try: # Our implementation x1, residuals1, rank1, s1 = gelsd(A, b) # LAPACK implementation m, n = np.array(A).shape minmn = min(m, n) maxmn = max(m, n) if len(b) < maxmn: b_padded = np.zeros(maxmn, dtype=np.float64) b_padded[:len(b)] = b b_arr = b_padded else: b_arr = np.array(b) x2, s2, rank2, info = lapack.dgelsd(A, b_arr, lwork=10000, size_iwork=10000) x2 = x2.ravel() # Compare results assert_allclose(x1, x2[:n], rtol=1e-12, atol=1e-12) assert_allclose(s1, s2[:minmn], rtol=1e-12, atol=1e-12) assert rank1 == rank2 except Exception as e: raise AssertionError(f"Failed for case: {name}\nError: {str(e)}") def test_gelsd_rcond(): A = [[0., 1., 0., 1., 2., 0.], [0., 2., 0., 0., 1., 0.], [1., 0., 1., 0., 0., 4.], [0., 0., 0., 2., 3., 0.]] A = np.array(A).T.tolist() b = [1, 0, 0, 0, 0, 0] # With rcond=-1, should give full rank x1, residuals1, rank1, s1 = gelsd(A, b, rcond=-1) assert rank1 == 4 # With default rcond, should detect rank deficiency x2, residuals2, rank2, s2 = gelsd(A, b) assert rank2 == 3 @pytest.mark.parametrize("m,n,n_rhs", [ (4, 2, 1), # Overdetermined, single RHS (4, 0, 1), # Empty columns (4, 2, 1), # Standard overdetermined (2, 4, 1), # Underdetermined ]) def test_gelsd_empty_and_shapes(m, n, n_rhs): """Test various matrix shapes including empty matrices""" # Create test matrices if m * n > 0: A = np.arange(m * n).reshape(m, n).tolist() else: A = np.zeros((m, n)).tolist() if m > 0: b = np.ones(m).tolist() else: b = np.ones(0).tolist() x, residuals, rank, s = gelsd(A, b) # Check dimensions assert len(x) == n assert len(s) == min(m, n) # Check rank assert rank == min(m, n) # For zero-sized matrices, solution should be zero if m == 0: assert_allclose(x, np.zeros(n)) # For overdetermined systems, check residuals if m > n and n > 0: r = np.array(b) - np.dot(A, x) expected_residuals = float(np.sum(r * r)) assert_allclose(residuals, expected_residuals, atol=1e-28) def test_gelsd_incompatible_dims(): """Test error handling for incompatible dimensions""" A = [[1.0, 2.0], [3.0, 4.0]] b = [1.0, 2.0, 3.0] # Wrong dimension with pytest.raises(ValueError): gelsd(A, b) def test_array_as_tridiagonals(): A = [[10.0, 2.0, 0.0, 0.0], [3.0, 10.0, 4.0, 0.0], [0.0, 1.0, 7.0, 5.0], [0.0, 0.0, 3.0, 4.0]] tridiagonals = array_as_tridiagonals(A) expect_diags = [[3.0, 1.0, 3.0], [10.0, 10.0, 7.0, 4.0], [2.0, 4.0, 5.0]] assert_allclose(tridiagonals[0], expect_diags[0], rtol=0, atol=0) assert_allclose(tridiagonals[1], expect_diags[1], rtol=0, atol=0) assert_allclose(tridiagonals[2], expect_diags[2], rtol=0, atol=0) A = np.array(A) tridiagonals = array_as_tridiagonals(A) assert_allclose(tridiagonals[0], expect_diags[0], rtol=0, atol=0) assert_allclose(tridiagonals[1], expect_diags[1], rtol=0, atol=0) assert_allclose(tridiagonals[2], expect_diags[2], rtol=0, atol=0) a, b, c = [3.0, 1.0, 3.0], [10.0, 10.0, 7.0, 4.0], [2.0, 4.0, 5.0] expect_mat = tridiagonals_as_array(a, b, c) assert_allclose(expect_mat, A, rtol=0, atol=0) d = [3.0, 4.0, 5.0, 6.0] solved_expect = [0.1487758945386064, 0.756120527306968, -1.001883239171375, 2.2514124293785316] assert_allclose(solve_tridiagonal(a, b, c, d), solved_expect, rtol=1e-12) def test_subset_matrix(): kijs = [[0, 0.00076, 0.00171], [0.00076, 0, 0.00061], [0.00171, 0.00061, 0]] expect = [[0, 0.00061], [0.00061, 0]] got = subset_matrix(kijs, [1,2]) assert_allclose(expect, got, atol=0, rtol=0) got = subset_matrix(kijs, slice(1, 3, 1)) assert_allclose(expect, got, atol=0, rtol=0) expect = [[0, 0.00171], [0.00171, 0]] got = subset_matrix(kijs, [0,2]) assert_allclose(expect, got, atol=0, rtol=0) got = subset_matrix(kijs, slice(0, 3, 2)) assert_allclose(expect, got, atol=0, rtol=0) expect = [[0, 0.00076], [0.00076, 0]] got = subset_matrix(kijs, [0,1]) assert_allclose(expect, got, atol=0, rtol=0) got = subset_matrix(kijs, slice(0, 2, 1)) assert_allclose(expect, got, atol=0, rtol=0) got = subset_matrix(kijs, [0,1, 2]) assert_allclose(kijs, got, atol=0, rtol=0) got = subset_matrix(kijs, slice(0, 3, 1)) assert_allclose(kijs, got, atol=0, rtol=0) def test_argsort1d(): def check_argsort1d(input_list, expected, error_message): numpy_argsort1d = lambda x: list(np.argsort(x)) assert argsort1d(input_list) == expected, error_message assert argsort1d(input_list) == numpy_argsort1d(input_list), error_message check_argsort1d([3, 1, 2], [1, 2, 0], "Failed on simple test case") check_argsort1d([-1, -3, -2], [1, 2, 0], "Failed with negative numbers") check_argsort1d([], [], "Failed on empty list") check_argsort1d([42], [0], "Failed with single element list") check_argsort1d([99, 21, 31, 80, 70], [1, 2, 4, 3, 0], "Mismatch with expected output") check_argsort1d([2, 3, 1, 5, 4], [2, 0, 1, 4, 3], "Mismatch with expected output") check_argsort1d([3.5, 1, 2.2], [1, 2, 0], "Failed with mixed floats and ints") check_argsort1d([0.1, 0.2, 0.3], [0, 1, 2], "Failed with floats") check_argsort1d([True, False, True], [1, 0, 2], "Failed with boolean values") check_argsort1d(['apple', 'banana', 'cherry'], [0, 1, 2], "Failed with strings") check_argsort1d([2, 3, 2, 3, 3], [0, 2, 1, 3, 4], "Failed with duplicate numbers") check_argsort1d([-3, -1, 0, 1, 3], [0, 1, 2, 3, 4], "Failed with negative and positive numbers") # infinities and nan behavior does not match # check_argsort1d([-np.inf, np.inf, np.nan, 0, -1], [0, 4, 3, 2, 1], "Failed with infinities and NaN") def test_eye(): # Test basic functionality assert eye(1) == [[1.0]] assert eye(2) == [[1.0, 0.0], [0.0, 1.0]] assert eye(3) == [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]] # Test with different dtypes assert eye(2, dtype=int) == [[1, 0], [0, 1]] assert eye(2, dtype=float) == [[1.0, 0.0], [0.0, 1.0]] # Test error cases with pytest.raises(ValueError): eye(0) # Zero size with pytest.raises(ValueError): eye(-1) # Negative size with pytest.raises(TypeError): eye(2.5) # Non-integer size # Test matrix properties def check_matrix_properties(matrix): N = len(matrix) # Check dimensions assert all(len(row) == N for row in matrix), "Matrix rows have inconsistent lengths" # Check diagonal elements assert all(matrix[i][i] == 1 for i in range(N)), "Diagonal elements are not 1" # Check off-diagonal elements assert all(matrix[i][j] == 0 for i in range(N) for j in range(N) if i != j), "Off-diagonal elements are not 0" # Test matrix properties for various sizes for size in [1, 2, 3, 4, 5, 10]: check_matrix_properties(eye(size)) # Test type consistency def check_type_consistency(matrix, expected_type): assert all(isinstance(x, expected_type) for row in matrix for x in row), f"Not all elements are of type {expected_type}" # Check type consistency for different dtypes check_type_consistency(eye(3, dtype=float), float) check_type_consistency(eye(3, dtype=int), int) def test_dot_product(): assert dot_product([1, 2, 3], [4, 5, 6]) == 32.0 assert dot_product([1, 0], [0, 1]) == 0.0 # Orthogonal vectors assert dot_product([1, 1], [1, 1]) == 2.0 # Parallel vectors assert_close(dot_product([0.1, 0.2], [0.3, 0.4]), 0.11) assert_close(dot_product([-1, -2], [3, 4]), -11.0) # Test properties of dot product def test_commutative(a, b): """Test if a·b = b·a""" assert_close(dot_product(a, b), dot_product(b, a), rtol=1e-14) def test_distributive(a, b, c): """Test if a·(b + c) = a·b + a·c""" # Create vector sum b + c vec_sum = [bi + ci for bi, ci in zip(b, c)] left = dot_product(a, vec_sum) right = dot_product(a, b) + dot_product(a, c) return assert_close(left, right, rtol=1e-14) # Test mathematical properties a, b, c = [1, 2], [3, 4], [5, 6] test_commutative(a, b) test_distributive(a, b, c) # Test error cases with pytest.raises(ValueError): dot_product([1, 2], [1, 2, 3]) # Different lengths def test_matrix_vector_dot(): """Test the matrix-vector dot product function""" # Basic multiplication matrix = [[1, 2], [3, 4]] vector = [1, 2] result = matrix_vector_dot(matrix, vector) assert_close1d(result, [5, 11]) # Identity matrix matrix = [[1, 0], [0, 1]] assert_close1d(matrix_vector_dot(matrix, [2, 3]), [2, 3]) # Zero matrix matrix = [[0, 0], [0, 0]] assert_close1d(matrix_vector_dot(matrix, [1, 1]), [0, 0]) # Rectangular matrix (more rows than columns) matrix = [[1, 2], [3, 4], [5, 6]] vector = [1, 2] assert_close1d(matrix_vector_dot(matrix, vector), [5, 11, 17]) # Error cases with pytest.raises(ValueError): matrix_vector_dot([[1, 2], [3, 4]], [1, 2, 3]) # Incompatible dimensions def test_transpose(): # Empty matrix and empty rows assert transpose([]) == [] assert transpose([[]]) == [] # 1x1 matrix assert transpose([[1]]) == [[1]] # 2x2 matrix assert transpose([[1, 2], [3, 4]]) == [[1, 3], [2, 4]] # 3x3 matrix assert transpose([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) == [[1, 4, 7], [2, 5, 8], [3, 6, 9]] # Rectangular matrices assert transpose([[1, 2, 3], [4, 5, 6]]) == [[1, 4], [2, 5], [3, 6]] # Single row/column assert transpose([[1, 2, 3]]) == [[1], [2], [3]] assert transpose([[1], [2], [3]]) == [[1, 2, 3]] # Mixed types result = transpose([[1, 2.5], [3, 4.2]]) assert result[0][0] == 1 assert abs(result[1][1] - 4.2) < 1e-10 # Float comparison with tolerance def test_matrix_multiply(): # 2x2 matrices A = [[1, 2], [3, 4]] B = [[5, 6], [7, 8]] result = matrix_multiply(A, B) expect = [[19, 22], [43, 50]] assert_close2d(result, expect) # Identity matrix I = [[1, 0], [0, 1]] assert_close2d(matrix_multiply(A, I), A) assert_close2d(matrix_multiply(I, A), A) # Zero matrix Z = [[0, 0], [0, 0]] assert_close2d(matrix_multiply(A, Z), Z) # Different shapes A = [[1, 2, 3], [4, 5, 6]] # 2x3 B = [[7, 8], [9, 10], [11, 12]] # 3x2 result = matrix_multiply(A, B) expect = [[58, 64], [139, 154]] assert_close2d(result, expect) # Very small numbers A = [[1e-15, 1e-15], [1e-15, 1e-15]] result = matrix_multiply(A, A) expect = [[2e-30, 2e-30], [2e-30, 2e-30]] assert_close2d(result, expect) # Very large numbers A = [[1e15, 1e15], [1e15, 1e15]] result = matrix_multiply(A, A) expect = [[2e30, 2e30], [2e30, 2e30]] assert_close2d(result, expect) # Mixed scales A = [[1e10, 1e-10], [1e-10, 1e10]] result = matrix_multiply(A, A) expect = [[1e20 + 1e-20, 2], [2, 1e20 + 1e-20]] assert_close2d(result, expect) # Single-element matrices A = [[2]] B = [[3]] C = [[6.0]] assert matrix_multiply(A, B) == C # Large matrices (should not raise error) A = [[i for i in range(10)] for _ in range(10)] B = [[i for i in range(10)] for _ in range(10)] C = matrix_multiply(A, B) assert len(C) == 10 and len(C[0]) == 10 # Empty matrices with pytest.raises(ValueError): matrix_multiply([], []) with pytest.raises(ValueError): matrix_multiply([[]], [[]]) with pytest.raises(ValueError): matrix_multiply([], [[1]]) with pytest.raises(ValueError): matrix_multiply([[1]], []) # Incompatible dimensions with pytest.raises(ValueError): matrix_multiply([[1, 2]], [[1], [2], [3]]) # Irregular matrices with pytest.raises(ValueError): matrix_multiply([[1, 2], [1]], [[1, 2]]) with pytest.raises(ValueError): matrix_multiply([[1, 2]], [[1], [1, 2]]) # Non-numeric values with pytest.raises(TypeError): A = [[1, 2, 'a']] B = [[4, 5], [6, 7], [8, 9]] matrix_multiply(A, B) def test_sum_matrix_rows(): """Test row-wise matrix summation""" # Basic functionality assert_close1d(sum_matrix_rows([[1, 2, 3], [4, 5, 6]]), [6.0, 15.0]) assert_close1d(sum_matrix_rows([[1], [2]]), [1.0, 2.0]) # Handle zeros assert_close1d(sum_matrix_rows([[0, 0], [0, 0]]), [0.0, 0.0]) # Mixed positive and negative assert_close1d(sum_matrix_rows([[-1, 2], [3, -4]]), [1.0, -1.0]) # Large numbers assert_close1d(sum_matrix_rows([[1e15, 1e15], [1e15, 1e15]]), [2e15, 2e15]) # Small numbers assert_close1d(sum_matrix_rows([[1e-15, 1e-15], [1e-15, 1e-15]]), [2e-15, 2e-15]) # Test with single row assert sum_matrix_cols([[1, 2, 3]]) == [1.0, 2.0, 3.0] # Test with single column assert sum_matrix_cols([[1], [2], [3]]) == [6.0] # Error cases with pytest.raises(ValueError): sum_matrix_rows([]) # Empty matrix with pytest.raises(ValueError): sum_matrix_rows([[]]) # Empty rows with pytest.raises(ValueError): sum_matrix_rows([[1, 2], [1]]) # Irregular rows # Test non-numeric values with pytest.raises(TypeError): sum_matrix_cols([[1, 'a'], [2, 3]]) def test_sum_matrix_cols(): """Test column-wise matrix summation""" # Basic functionality assert_close1d(sum_matrix_cols([[1, 2, 3], [4, 5, 6]]), [5.0, 7.0, 9.0]) assert_close1d(sum_matrix_cols([[1], [2]]), [3.0]) # Test with single row assert sum_matrix_rows([[1, 2, 3]]) == [6.0] # Test with single column assert sum_matrix_rows([[1], [2], [3]]) == [1.0, 2.0, 3.0] # Handle zeros assert_close1d(sum_matrix_cols([[0, 0], [0, 0]]), [0.0, 0.0]) # Mixed positive and negative assert_close1d(sum_matrix_cols([[-1, 2], [3, -4]]), [2.0, -2.0]) # Large numbers assert_close1d(sum_matrix_cols([[1e15, 1e15], [1e15, 1e15]]), [2e15, 2e15]) # Small numbers assert_close1d(sum_matrix_cols([[1e-15, 1e-15], [1e-15, 1e-15]]), [2e-15, 2e-15]) # Error cases with pytest.raises(ValueError): sum_matrix_cols([]) # Empty matrix with pytest.raises(ValueError): sum_matrix_cols([[]]) # Empty rows with pytest.raises(ValueError): sum_matrix_cols([[1, 2], [1]]) # Irregular rows with pytest.raises(TypeError): sum_matrix_rows([[1, 'a'], [2, 3]]) def test_scalar_add_matrices(): """Test matrix addition functionality""" # Basic functionality assert_close2d(scalar_add_matrices([[1, 2], [3, 4]], [[5, 6], [7, 8]]), [[6.0, 8.0], [10.0, 12.0]]) # Single element matrices assert_close2d(scalar_add_matrices([[1]], [[2]]), [[3.0]]) # Test with zeros assert_close2d(scalar_add_matrices([[0, 0], [0, 0]], [[0, 0], [0, 0]]), [[0.0, 0.0], [0.0, 0.0]]) # Mixed positive and negative assert_close2d(scalar_add_matrices([[-1, 2], [3, -4]], [[1, -2], [-3, 4]]), [[0.0, 0.0], [0.0, 0.0]], atol=1e-14) # Large numbers assert_close2d(scalar_add_matrices([[1e15, 1e15], [1e15, 1e15]], [[1e15, 1e15], [1e15, 1e15]]), [[2e15, 2e15], [2e15, 2e15]]) # Small numbers assert_close2d(scalar_add_matrices([[1e-15, 1e-15], [1e-15, 1e-15]], [[1e-15, 1e-15], [1e-15, 1e-15]]), [[2e-15, 2e-15], [2e-15, 2e-15]]) # Different shapes of matrices rect1 = [[1, 2, 3], [4, 5, 6]] rect2 = [[7, 8, 9], [10, 11, 12]] assert_close2d(scalar_add_matrices(rect1, rect2), [[8.0, 10.0, 12.0], [14.0, 16.0, 18.0]]) # Error cases with pytest.raises(ValueError): scalar_add_matrices([], []) # Empty matrices with pytest.raises(ValueError): scalar_add_matrices([[]], [[]]) # Empty rows with pytest.raises(ValueError): scalar_add_matrices([[1, 2], [1]], [[1, 2], [3, 4]]) # Irregular rows A with pytest.raises(ValueError): scalar_add_matrices([[1, 2], [3, 4]], [[1, 2], [3]]) # Irregular rows B with pytest.raises(ValueError): scalar_add_matrices([[1, 2]], [[1, 2, 3]]) # Incompatible shapes with pytest.raises(TypeError): scalar_add_matrices([[1, 'a']], [[1, 2]]) # Invalid type def test_scalar_subtract_matrices(): """Test matrix subtraction functionality""" # Basic functionality assert_close2d(scalar_subtract_matrices([[1, 2], [3, 4]], [[5, 6], [7, 8]]), [[-4.0, -4.0], [-4.0, -4.0]]) # Single element matrices assert_close2d(scalar_subtract_matrices([[1]], [[2]]), [[-1.0]]) # Test with zeros assert_close2d(scalar_subtract_matrices([[0, 0], [0, 0]], [[0, 0], [0, 0]]), [[0.0, 0.0], [0.0, 0.0]]) # Mixed positive and negative assert_close2d(scalar_subtract_matrices([[-1, 2], [3, -4]], [[1, -2], [-3, 4]]), [[-2.0, 4.0], [6.0, -8.0]]) # Large numbers assert_close2d(scalar_subtract_matrices([[1e15, 1e15], [1e15, 1e15]], [[1e15, 1e15], [1e15, 1e15]]), [[0.0, 0.0], [0.0, 0.0]]) # Small numbers assert_close2d(scalar_subtract_matrices([[1e-15, 1e-15], [1e-15, 1e-15]], [[1e-15, 1e-15], [1e-15, 1e-15]]), [[0.0, 0.0], [0.0, 0.0]]) # Different shapes of matrices rect1 = [[1, 2, 3], [4, 5, 6]] rect2 = [[7, 8, 9], [10, 11, 12]] assert_close2d(scalar_subtract_matrices(rect1, rect2), [[-6.0, -6.0, -6.0], [-6.0, -6.0, -6.0]]) # Error cases with pytest.raises(ValueError): scalar_subtract_matrices([], []) # Empty matrices with pytest.raises(ValueError): scalar_subtract_matrices([[]], [[]]) # Empty rows with pytest.raises(ValueError): scalar_subtract_matrices([[1, 2], [1]], [[1, 2], [3, 4]]) # Irregular rows A with pytest.raises(ValueError): scalar_subtract_matrices([[1, 2], [3, 4]], [[1, 2], [3]]) # Irregular rows B with pytest.raises(ValueError): scalar_subtract_matrices([[1, 2]], [[1, 2, 3]]) # Incompatible shapes with pytest.raises(TypeError): scalar_subtract_matrices([[1, 'a']], [[1, 2]]) # Invalid type def test_scalar_multiply_matrix(): """Test matrix scalar multiplication functionality""" # Basic functionality assert_close2d(scalar_multiply_matrix(2.0, [[1, 2], [3, 4]]), [[2.0, 4.0], [6.0, 8.0]]) # Single element matrices assert_close2d(scalar_multiply_matrix(3.0, [[2]]), [[6.0]]) # Test with zeros assert_close2d(scalar_multiply_matrix(0.0, [[1, 2], [3, 4]]), [[0.0, 0.0], [0.0, 0.0]]) # Test with negative scalar assert_close2d(scalar_multiply_matrix(-1.0, [[1, 2], [3, 4]]), [[-1.0, -2.0], [-3.0, -4.0]]) # Large numbers assert_close2d(scalar_multiply_matrix(1e15, [[1, 2], [3, 4]]), [[1e15, 2e15], [3e15, 4e15]]) # Small numbers assert_close2d(scalar_multiply_matrix(1e-15, [[1, 2], [3, 4]]), [[1e-15, 2e-15], [3e-15, 4e-15]]) # Rectangle matrix rect = [[1, 2, 3], [4, 5, 6]] assert_close2d(scalar_multiply_matrix(2.0, rect), [[2.0, 4.0, 6.0], [8.0, 10.0, 12.0]]) # Error cases with pytest.raises(ValueError): scalar_multiply_matrix(2.0, []) # Empty matrix with pytest.raises(ValueError): scalar_multiply_matrix(2.0, [[]]) # Empty rows with pytest.raises(TypeError): scalar_multiply_matrix(2.0, [[1, 'a']]) # Invalid type def test_scalar_divide_matrix(): """Test matrix scalar division functionality""" # Basic functionality assert_close2d(scalar_divide_matrix(2.0, [[2, 4], [6, 8]]), [[1.0, 2.0], [3.0, 4.0]]) # Single element matrices assert_close2d(scalar_divide_matrix(2.0, [[4]]), [[2.0]]) # Test with ones (identity case) assert_close2d(scalar_divide_matrix(1.0, [[1, 2], [3, 4]]), [[1.0, 2.0], [3.0, 4.0]]) # Test with negative scalar assert_close2d(scalar_divide_matrix(-2.0, [[2, 4], [6, 8]]), [[-1.0, -2.0], [-3.0, -4.0]]) # Large numbers assert_close2d(scalar_divide_matrix(1e15, [[1e15, 2e15], [3e15, 4e15]]), [[1.0, 2.0], [3.0, 4.0]]) # Small numbers assert_close2d(scalar_divide_matrix(1e-15, [[1e-15, 2e-15], [3e-15, 4e-15]]), [[1.0, 2.0], [3.0, 4.0]]) # Rectangle matrix rect = [[2, 4, 6], [8, 10, 12]] assert_close2d(scalar_divide_matrix(2.0, rect), [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]) # Error cases with pytest.raises(ValueError): scalar_divide_matrix(2.0, []) # Empty matrix with pytest.raises(ValueError): scalar_divide_matrix(2.0, [[]]) # Empty rows with pytest.raises(TypeError): scalar_divide_matrix(2.0, [[1, 'a']]) # Invalid type with pytest.raises(TypeError): scalar_divide_matrix('2', [[1, 2]]) # Invalid scalar type with pytest.raises(ZeroDivisionError): scalar_divide_matrix(0.0, [[1, 2]]) # Division by zero def test_sort_paired_lists(): assert sort_paired_lists([3, 1, 2], ['c', 'a', 'b']) == ([1, 2, 3], ['a', 'b', 'c']) assert sort_paired_lists([], []) == ([], []) assert sort_paired_lists([2, 2, 1], ['a', 'b', 'c']) == ([1, 2, 2], ['c', 'a', 'b']) assert sort_paired_lists([-3, -1, -2], ['c', 'a', 'b']) == ([-3, -2, -1], ['c', 'b', 'a']) temps = [300.5, 100.1, 200.7] props = ['hot', 'cold', 'warm'] assert sort_paired_lists(temps, props) == ([100.1, 200.7, 300.5], ['cold', 'warm', 'hot']) with pytest.raises(ValueError): # Test 6: Unequal length lists sort_paired_lists([1, 2], [1])