Ask Your Question

Revision history [back]

click to hide/show revision 1
initial version

Issue porting solve utility from sympy to sagemath

I'm trying to port a utility from sympy to sagemath.

Here is my sympy code:

from sympy import symbols, Eq, solve, Matrix, pprint
x, y, z = symbols('x y z')

def sympy_solve(augmented_matrix):
    # Get the number of variables
    num_variables = augmented_matrix.shape[1] - 1

    # Generate symbols for variables
    variables = symbols('x:' + str(num_variables))

    # Extract coefficients and constants from the augmented matrix
    coefficients = augmented_matrix[:, :-1]
    constants = augmented_matrix[:, -1]

    # Create equations from the coefficients and constants
    equations = []
    for i in range(len(constants)):
        equation = Eq(sum(coefficients[i, j] * variables[j] for j in range(num_variables)), constants[i])
        equations.append(equation)

    # Solve the equations
    solution = solve(equations, variables, dict=True)
    return solution

When I run it, I get the result I expect:

A = Matrix([
    [1,1,1,-1,1],
    [0,1,-1,1,-1],
    [3,0,6,-6,6],
    [0,-1,1,-1,1]
])
print(sympy_solve(A))
# [{
#   x0: -2*x2 + 2*x3 + 2, 
#   x1: x2 - x3 - 1
# }]

My current solution in sagemath looks like this:

def sagemath_solve(augmented_matrix):

    A = augmented_matrix[:, :-1]
    Y = augmented_matrix[:, -1]

    m, n = A.dimensions()
    p, q = Y.dimensions()

    if m!=p:
        raise RuntimeError("The matrices have different numbers of rows")
    X = vector([var("x_{}".format(i)) for i in [1..n]])

    sols = []
    for j in range(q):
        system = [A[i]*X==Y[i,j] for i in range(m)]
        sols += solve(system, *X)
    return sols

However, it is returning rationals:

A = matrix([
    [1,1,1,-1,1],
    [0,1,-1,1,-1],
    [3,0,6,-6,6],
    [0,-1,1,-1,1]
])
print(sagemath_solve(A))

# [[
#   x_1 == 2*r1 - 2*r2 + 2, 
#   x_2 == -r1 + r2 - 1, 
#   x_3 == r2, 
#   x_4 == r1
# ]]

Any thoughts on how to fix this?