Python for Bioinformatics: Cramer's rule calculation

Friday, June 24, 2011

Cramer's rule calculation

We have algebra homework that involves using Cramer's rule to solve not only 2 x 2 but also 3 x 3 systems. It seems kind of silly since this method is overkill for 2 x 2, and would never be used for 4 x 4 or larger.

(Note on the wikipedia article, start about halfway down, where it says "Explicit formulas for small systems".)

Also, and this gets closer to the point, drilling by solving 3 x 3 matrices is not really about the rule, which is pretty simple. It's about making an easy problem harder by stuffing a lot of arithmetic into it. And to me, that is symptomatic of a big difficulty with math education as I'm encountering it through my son. Computers are much better at computing sums than humans. It's just silly to drill students on arithmetic. If you want to do something complicated, why not derive Cramer's rule?

So, I decided to write a solver for 3 x 3 systems in Python. I wouldn't say it's thoroughly debugged yet, so let me know if you run into a problem. With the example shown, I did get the same answer as this online calculator.

The first code segment contains the equations explicitly entered as an array. I'm sure you know how to modify it to read input from a file.

test.py

import numpy as np
import Cramer

def test_Cramer():
    L = [2, 3, 0, 5,
         1, 1, 1, 3,
         2,-1, 3, 7]   
    A = np.array(L)
    A.shape = (3,4)
    result = Cramer.solve(A)
    if result:
        x,y,z = result
        print 'solution'
        print 'x =', x
        print 'y =', y
        print 'z =', z, '\n'
        Cramer.check(A,x,y,z)

test_Cramer()

The output looks like this:

> python test.py 
solve
[[ 2  3  0  5]
 [ 1  1  1  3]
 [ 2 -1  3  7]] 

compute 3 x 3 det of
[[ 2  3  0]
 [ 1  1  1]
 [ 2 -1  3]]
D =  5 

compute 3 x 3 det of
[[ 5  3  0]
 [ 3  1  1]
 [ 7 -1  3]]
Dx =  14 

compute 3 x 3 det of
[[2 5 0]
 [1 3 1]
 [2 7 3]]
Dy =  -1 

compute 3 x 3 det of
[[ 2  3  5]
 [ 1  1  3]
 [ 2 -1  7]]
Dz =  2 

solution
x = 2.8
y = -0.2
z = 0.4 

check
row 0 = [2 3 0 5]
2.0*2.8 + 3.0*-0.2 + 0.0*0.4  = 5.0 

row 1 = [1 1 1 3]
1.0*2.8 + 1.0*-0.2 + 1.0*0.4  = 3.0 

row 2 = [ 2 -1  3  7]
2.0*2.8 + -1.0*-0.2 + 3.0*0.4  = 7.0

Cramer.py

import numpy as np

def det2x2(A, v=False):
    if v:  print 'compute 2 x 2 det of'
    if v:  print A
    assert A.shape == (2,2)
    return A[0][0]*A[1][1] - A[0][1]*A[1][0]

def det3x3(A):
    print 'compute 3 x 3 det of'
    print A
    assert A.shape == (3,3)
    a,b,c = A[0]
    c1 = a * det2x2(A[1:3,[1,2]])
    c2 = b * det2x2(A[1:3,[0,2]])
    c3 = c * det2x2(A[1:3,[0,1]])
    return c1 - c2 + c3
         
def solve(A):
    print 'solve'
    print A, '\n'
    assert A.shape == (3,4)
    D = det3x3(A[:,:3])
    print 'D = ', D, '\n'
    if D == 0:
        print 'no solution'
        return
    Dx = det3x3(A[:,[3,1,2]])
    print 'Dx = ', Dx, '\n'
    Dy = det3x3(A[:,[0,3,2]])
    print 'Dy = ', Dy, '\n'
    Dz = det3x3(A[:,[0,1,3]])
    print 'Dz = ', Dz, '\n'
    return Dx*1.0/D, Dy*1.0/D, Dz*1.0/D
    
def check(A,x,y,z):
    print 'check'
    for i,r in enumerate(A):
        print 'row', i, '=', r
        pL = list()
        for coeff,var in zip(r[:3],(x,y,z)):
            c = str(round(coeff,2))
            v = str(round(var,2))
            pL.append(c + '*' + v)
        print ' + '.join(pL),
        print ' =', r[0]*x + r[1]*y + r[2]*z, '\n'