The science/engineering application I'm working on has lots of linear algebra matrix multiplications, therefore I use Numpy matrices. However, there are many functions in python that interchangeably accept matrix or array types. Nice, no? Well, not really. Let me demonstrate the problem with an example:
from scipy.linalg import expm
from numpy import matrix
# Setup input variable as matrix
A = matrix([[ 0, -1.0, 0, 0],
[ 0, 0, 0, 1.0],
[ 0, 0, 0, 0],
[ 0, 0, 1.0, 0]])
# Do some computation with that input
B = expm(A)
b1 = B[0:2, 2:4]
b2 = B[2:4, 2:4].T
# Compute and Print the desired output
print "The innocent but wrong answer:"
print b2 * b1
print "The answer I should get:"
print matrix(b2) * matrix(b1)
When run you get:
The innocent but wrong answer:
[[-0.16666667 -0.5 ]
[ 0. 1. ]]
The answer I should get, since I expected everything to still be matrices:
[[ 0.33333333 0.5 ]
[ 0.5 1. ]]
Any tips or advice on how to avoid this sort of a mix up? Its really messy to keep wrapping variables in matrix() calls to ensure they still are matrices. It seems there is no standard in this regard, and so it can lead to bugs that are hard to detect.