Multidimensional matrices, Iverson notation inspired.

APL.py

# -*- coding: utf-8 -*-
import operator
import itertools

def main():
    shape = [3, 3, 2]
    values = [1, 3, 3, 3, 1, 3, 3, 3, 1, 2, 2, 2, 4, 2, 2, 4, 4, 2]

    a = numeric([1,2,3,4,5])
    print a
    # ⟮1 2 3 4 5⟯
    b = numeric([a*-2, a*-1, a*0, a*1, a*2])
    print b
    # ╭-2 -4 -6 -8 -10╮
    # │-1 -2 -3 -4  -5│
    # │ 0  0  0  0   0│
    # │ 1  2  3  4   5│
    # ╰ 2  4  6  8  10╯
    print fold_dimension(operator.add, b, 0)
    # ⟮-30 -15 0 15 30⟯
    print fold_dimension(operator.add, b, 1)
    # ⟮0 0 0 0 0⟯
    print diagonal(b, 0, 1)
    # ⟮-2 -2 0 4 10⟯
    print index_dimension(b, 0, 3)
    # ⟮-8 -4 0 4 8⟯
    print index_dimension(b, 1, 3)
    # ⟮1 2 3 4 5⟯
    print reverse_dimension(b, 0)
    # ╭-10 -8 -6 -4 -2╮
    # │ -5 -4 -3 -2 -1│
    # │  0  0  0  0  0│
    # │  5  4  3  2  1│
    # ╰ 10  8  6  4  2╯
    print reverse_dimension(b, 1)
    # ╭ 2  4  6  8  10╮
    # │ 1  2  3  4   5│
    # │ 0  0  0  0   0│
    # │-1 -2 -3 -4  -5│
    # ╰-2 -4 -6 -8 -10╯
    print roll_dimension(b, 0, 1)
    # ╭-4 -6 -8 -10 -2╮
    # │-2 -3 -4  -5 -1│
    # │ 0  0  0   0  0│
    # │ 2  3  4   5  1│
    # ╰ 4  6  8  10  2╯
    print roll_dimension(b, 0, 3)
    # ╭-8 -10 -2 -4 -6╮
    # │-4  -5 -1 -2 -3│
    # │ 0   0  0  0  0│
    # │ 4   5  1  2  3│
    # ╰ 8  10  2  4  6╯
    print roll_dimension(b, 1, 1)
    # ╭-1 -2 -3 -4  -5╮
    # │ 0  0  0  0   0│
    # │ 1  2  3  4   5│
    # │ 2  4  6  8  10│
    # ╰-2 -4 -6 -8 -10╯
    print roll_dimension(b, 1, 3)
    # ╭ 1  2  3  4   5╮
    # │ 2  4  6  8  10│
    # │-2 -4 -6 -8 -10│
    # │-1 -2 -3 -4  -5│
    # ╰ 0  0  0  0   0╯
    print transpose(b)
    # ╭ -2 -1 0 1  2╮
    # │ -4 -2 0 2  4│
    # │ -6 -3 0 3  6│
    # │ -8 -4 0 4  8│
    # ╰-10 -5 0 5 10╯
    print b
    # ╭-2 -4 -6 -8 -10╮
    # │-1 -2 -3 -4  -5│
    # │ 0  0  0  0   0│
    # │ 1  2  3  4   5│
    # ╰ 2  4  6  8  10╯
    print slice_dimension(slice_dimension(b, 0, 1, 4), 1, 1, 4)
    # ╭-2 -3 -4╮
    # │ 0  0  0│
    # ╰ 2  3  4╯

class Numeric(object):
    def __init__(self, shape, values):
        self.shape = shape
        self.values = values

    def __pos__(self):
        return self

    def __neg__(self):
        return unary_operation(operator.neg, self)

    def __add__(self, other):
        return binary_operation(operator.add, self, other)

    def __radd__(self, other):
        return binary_operation(operator.add, other, self)

    def __sub__(self, other):
        return binary_operation(operator.sub, self, other)

    def __rsub__(self, other):
        return binary_operation(operator.sub, other, self)

    def __mul__(self, other):
        return binary_operation(operator.mul, self, other)

    def __rmul__(self, other):
        return binary_operation(operator.mul, other, self)

    def __repr__(self):
        return "".join(generate_table(self.shape, self.values))

def generate_table(shape, values):
    indexer = indexer_vector(shape, shape)
    even = indexer[0::2]
    odds = indexer[1::2]
    rows = []
    stride = 1
    for count, _ in even:
        stride *= count
        rows.append(stride)
    row_stride = stride
    cols = []
    for count, _ in odds:
        stride *= count
        cols.append(stride)
    col_stride = stride
    while len(rows) > len(cols):
        cols.append(col_stride)

    col_size = [1 for col in range(row_stride)]
    for j, index in enumerate(take_indices(even + odds)):
        k = j % len(col_size)
        assert index < len(values), ("FAIL", j, index, len(values))
        col_size[k] = max(col_size[k], len(repr(values[index])))

    view1 = [u"\u256D", u"\u2502", u"\u2570", u"\u256E", u"\u2502", u"\u256F"]
    view2 = [u"\u239B", u"\u239C", u"\u239D", u"\u239E", u"\u239F", u"\u23A0"]
    view = view1

    sp = ""
    for j, i in enumerate(take_indices(even + odds)):
        k = j % len(col_size)
        yield sp
        for z, n in enumerate(reversed(rows)):
            ## z == 2, j < 32
            if j % n == 0:
                p = (j - j % row_stride)
                uc = (p % cols[-1-z] == 0)
                lc = ((p + row_stride) % cols[-1-z] == 0)
                if uc and lc:
                    yield u"\u27EE".encode('utf-8')
                elif uc:
                    yield view[0].encode('utf-8')
                elif lc:
                    yield view[2].encode('utf-8')
                else:
                    yield view[1].encode('utf-8')
        yield repr(values[i]).rjust(col_size[k])
        sp = " "
        for z, n in enumerate(rows):
            if (j+1) % n == 0:
                p = (j - j % row_stride)
                uc = ((p + row_stride) % cols[z] == 0)
                lc = (p % cols[z] == 0)
                if uc and lc:
                    yield u"\u27EF".encode('utf-8')
                elif uc:
                    yield view[5].encode('utf-8')
                elif lc:
                    yield view[3].encode('utf-8')
                else:
                    yield view[4].encode('utf-8')
        if (j+1) % len(col_size) == 0:
            sp = "\n"

def unary_operation(op, self):
    return Numeric(self.shape, [op(value) for value in self.values])

def binary_operation(op, self, other):
    self = to_numeric(self)
    other = to_numeric(other)
    shape = []
    for count1, count2 in zip_shapes(self.shape, other.shape):
        k = max(count1, count2)
        shape.append(k)
        assert count1 == count2 or count1 == 1 or count2 == 1, "shape error"
    indexer1 = indexer_vector(self.shape,  shape)
    indexer2 = indexer_vector(other.shape, shape)
    indices1 = take_indices(indexer1, 0)
    indices2 = take_indices(indexer2, 0)
    values = [op(self.values[i], other.values[j])
        for i, j in itertools.izip(indices1, indices2)]
    return Numeric(shape, values)

def fold_dimension(op, self, rank=0):
    if not (0 <= rank < len(self.shape)):
        return self
    indexer = indexer_vector(self.shape, self.shape)
    count, stride = indexer.pop(rank)
    if stride == 0:
        return self
    width = count*stride
    shape  = indexer_shape(indexer)
    values = [reduce(op,
        (self.values[j] for j in range(i, i+width, stride)))
        for i in take_indices(indexer, 0)]
    if len(shape) == 0 or all(i == 1 for i in shape):
        return values[0]
    return Numeric(shape, values)

def index_dimension(self, rank, index):
    if not (0 <= rank < len(self.shape)):
        assert index == 0, "index range error"
        return self
    indexer = indexer_vector(self.shape, self.shape)
    count, stride = indexer.pop(rank)
    assert 0 <= index < count, "index range error"
    shape  = indexer_shape(indexer)
    values = [self.values[i] for i in take_indices(indexer, index*stride)]
    if len(shape) == 0 or all(i == 1 for i in shape):
        return values[0]
    return Numeric(shape, values)

def diagonal(self, rank1, rank2):
    rank1, rank2 = min(rank1, rank2), max(rank1, rank2)
    assert rank_count(self.shape, rank1) == rank_count(self.shape, rank2), "shape error"
    if not (0 <= rank1 < rank2 < len(self.shape)):
        return self
    indexer = indexer_vector(self.shape, self.shape)
    indexer[rank1] = (indexer[rank1][0],
        indexer[rank1][1] + indexer.pop(rank2)[1])
    shape = self.shape[:rank2] + self.shape[rank2+1:]
    values = [self.values[i] for i in take_indices(indexer, 0)]
    return Numeric(shape, values)

def transpose(self, depth=2):
    indexer = indexer_vector(self.shape, self.shape)
    indexer.extend((1,0) for i in range(len(indexer), depth))
    indexer[:depth] = reversed(indexer[:depth])
    shape  = indexer_shape(indexer)
    values = [self.values[i] for i in take_indices(indexer, 0)]
    return Numeric(shape, values)

def reverse_dimension(self, rank):
    if not (0 <= rank < len(self.shape)):
        return self
    offset = 0
    indexer = indexer_vector(self.shape, self.shape)
    shuffles = []
    for index, (count, stride) in enumerate(indexer):
        if index != rank:
            shuffle = (count, stride, stride-stride*count, count-1)
        else:
            offset += stride*(count-1)
            shuffle = (count, -stride, -stride+stride*count, count-1)
        shuffles.append(shuffle)
    values = [self.values[i] for i in take_shuffle(shuffles, offset)]
    return Numeric(self.shape, values)

def roll_dimension(self, rank, amount):
    if not (0 <= rank < len(self.shape)):
        return self
    offset = 0
    indexer = indexer_vector(self.shape, self.shape)
    shuffles = []
    for index, (count, stride) in enumerate(indexer):
        if index != rank:
            shuffle = (count, stride, stride-stride*count, count-1)
        else:
            amount %= count
            offset += stride*amount
            shuffle = (count, stride, stride-stride*count, count-1-amount)
        shuffles.append(shuffle)
    values = [self.values[i] for i in take_shuffle(shuffles, offset)]
    return Numeric(self.shape, values)

def slice_dimension(self, rank, start, stop):
    if not (0 <= rank < len(self.shape)):
        return self
    offset = 0
    indexer = indexer_vector(self.shape, self.shape)
    shuffles = []
    for index, (count, stride) in enumerate(indexer):
        if index != rank:
            shuffle = (count, stride, stride-stride*count, count-1)
        else:
            assert 0 <= start < stop <= count, "slice error"
            offset += stride*start
            count = stop-start
            shuffle = (count, stride, stride-stride*count, count-1)
        shuffles.append(shuffle)
    shape  = indexer_shape(shuffles)
    values = [self.values[i] for i in take_shuffle(shuffles, offset)]
    return Numeric(shape, values)

def numeric(objects):
    objects = [to_numeric(n) for n in objects]
    shape = []
    for num in objects:
        for i, count in enumerate(num.shape):
            if i >= len(shape):
                shape.append(count)
                continue
            assert shape[i] == count or count == 1, "shape error"
            shape[i] = max(shape[i], count)
    values = [num.values[i]
        for num in objects
        for i in take_indices(indexer_vector(num.shape, shape), 0)]
    shape.append(len(objects))
    return Numeric(shape, values)
    
def indexer_vector(shape, newshape, stride=1):
    vector = []
    for count, newcount in zip_shapes(shape, newshape):
        assert count == newcount or count == 1
        vector.append((newcount, (stride if count > 1 else 0)))
        stride *= count
    return vector

def indexer_shape(vector):
    shape = [e[0] for e in vector]
    while len(shape) > 0 and shape[-1] == (1,0):
        shape.pop()
    return shape

def take_indices(vector, offset=0):
    state = [0] * len(vector)
    while True:
        yield offset
        for i in range(len(vector)):
            offset   += vector[i][1]
            state[i] += 1
            if state[i] < vector[i][0]:
                break
            else:
                offset  -= state[i]*vector[i][1]
                state[i] = 0
        else:
            break

def take_shuffle(shuffles, offset=0):
    state = [0] * len(shuffles)
    while True:
        yield offset
        for i in range(len(shuffles)):
            count, step, skip, skipi = shuffles[i]
            offset += (step if state[i] != skipi else skip)
            state[i] += 1
            if state[i] < count:
                break
            else:
                state[i] = 0
        else:
            break

def rank_count(shape, rank):
    return shape[rank] if rank < len(shape) else 1

def zip_shapes(shape1, shape2):
    return itertools.izip_longest(shape1, shape2, fillvalue=1)

def to_numeric(value):
    if isinstance(value, (int, float, bool)):
        return Numeric([], [value])
    assert isinstance(value, Numeric), "requires support for this value"
    return value

if __name__=="__main__":
    main()