Answer a question

I'd like to be able to do matrix operations using a Python class as the elements—in this case, a simple Galois field implementation. It implements the necessary __add__, __mul__, __sub__ etc.

At first, I thought this should be possible with numpy arrays, using the dtype parameter, but from the dtype documentation, it seems that dtype can't be an arbitrary Python class. For example, I have a class Galois which does operations modulo 2:

>>> from galois import Galois
>>> Galois(1) + Galois(0)
Galois(1)
>>> Galois(1) + Galois(1)
Galois(0)

I can try to use this in numpy:

>>> import numpy as np
>>> a = np.identity(4, Galois)
>>> a
array([[1, 0, 0, 0],
       [0, 1, 0, 0],
       [0, 0, 1, 0],
       [0, 0, 0, 1]], dtype=object)

But if I do operations on the matrices, the elements aren't following the methods of my class:

>>> b = np.identity(4, Galois)
>>> a+b
array([[2, 0, 0, 0],
       [0, 2, 0, 0],
       [0, 0, 2, 0],
       [0, 0, 0, 2]], dtype=object)

Is there any way to make this work with numpy?

Is there any other Python matrix library that can do matrix operations (including inversion) on an arbitrary number-like class?

Update

Thanks for the answers so far. But I'm still not able to really use it as I hoped. Adds and multiplies seem good, but not matrix inversion. For example, let's try to get the AES inverse S-box affine transform matrix from the forward S-box affine transform matrix.

class Galois(object):
    MODULO = 2

    def __init__(self, val):
        self.val = int(val) % self.MODULO

    def __add__(self, val):
        return self.__class__((self.val + int(val)) % self.MODULO)
    def __sub__(self, val):
        return self.__class__((self.val - int(val)) % self.MODULO)
    def __mul__(self, val):
        return self.__class__((self.val * int(val)) % self.MODULO)
    def __int__(self):
        return self.val
    def __repr__(self):
        return "%s(%d)" % (self.__class__.__name__, self.val)
    def __float__(self):
        return float(self.val)

if __name__ == "__main__":
    import numpy as np

    Gv = np.vectorize(Galois)

    a = Gv(np.identity(8)) + Gv(np.eye(8,8,-1)) + Gv(np.eye(8,8,-2)) + Gv(np.eye(8,8,-3)) + Gv(np.eye(8,8,-4)) + Gv(np.eye(8,8,4)) + Gv(np.eye(8,8,5)) + Gv(np.eye(8,8,6)) + Gv(np.eye(8,8,7))
    print np.matrix(a)
    print np.matrix(a).I

The result:

[[Galois(1) Galois(0) Galois(0) Galois(0) Galois(1) Galois(1) Galois(1)
  Galois(1)]
 [Galois(1) Galois(1) Galois(0) Galois(0) Galois(0) Galois(1) Galois(1)
  Galois(1)]
 [Galois(1) Galois(1) Galois(1) Galois(0) Galois(0) Galois(0) Galois(1)
  Galois(1)]
 [Galois(1) Galois(1) Galois(1) Galois(1) Galois(0) Galois(0) Galois(0)
  Galois(1)]
 [Galois(1) Galois(1) Galois(1) Galois(1) Galois(1) Galois(0) Galois(0)
  Galois(0)]
 [Galois(0) Galois(1) Galois(1) Galois(1) Galois(1) Galois(1) Galois(0)
  Galois(0)]
 [Galois(0) Galois(0) Galois(1) Galois(1) Galois(1) Galois(1) Galois(1)
  Galois(0)]
 [Galois(0) Galois(0) Galois(0) Galois(1) Galois(1) Galois(1) Galois(1)
  Galois(1)]]
[[ 0.4  0.4 -0.6  0.4  0.4 -0.6  0.4 -0.6]
 [-0.6  0.4  0.4 -0.6  0.4  0.4 -0.6  0.4]
 [ 0.4 -0.6  0.4  0.4 -0.6  0.4  0.4 -0.6]
 [-0.6  0.4 -0.6  0.4  0.4 -0.6  0.4  0.4]
 [ 0.4 -0.6  0.4 -0.6  0.4  0.4 -0.6  0.4]
 [ 0.4  0.4 -0.6  0.4 -0.6  0.4  0.4 -0.6]
 [-0.6  0.4  0.4 -0.6  0.4 -0.6  0.4  0.4]
 [ 0.4 -0.6  0.4  0.4 -0.6  0.4 -0.6  0.4]]

Not the result I hoped for. It seems that for the matrix inversion, numpy just converts the matrix to floats, then does the inversion with plain real numbers.

Answers

You can use object as the dtype, which will allow arbitrary Python objects. I don't think there's any way of specializing a numpy array to accept only one particular class of Python object.