what if we stopped using multi-dimensionnal arrays to store views on linear arrays ?

What of my first lesson of C was about matrices in C. Why ? I was studying physics and we were taught what we needed.
Since I left school, this first lesson still seems weired nowadays. Take for instance a bitmap, either the one from a picture or the raw video RAM you use, the matrice is a view.

You want the point at column 10 row 21 on a 1980 x 1200 screen ? Just get the value at the offset of 21 x 1980 + 10 and you have your pixel in 2D notation.

So that's exactly what I coded : a matrix that is a view on a linear table hence kicking a lot of possible optimization :
  • first think : by allocating a linear array you are less likely to spread your allocation on faraway memory pages (less page fault)
  • for rank by rank operation (logical, addition, substraction, research of occurences) you can probably use SIMD optimization (requires a better backend than a « list » (ex : array))
In my actual life I never saw any usage of any matrices of more than 2 dimensions, so as a lazy person I did not generalized for nth dimension. During the coding, I decided to have both
matrix[x][y]
and
matrix[(x,y)]
notation being equivalent. However in the process, I discovered that for addressing it may change nothing, but for slice it is another story :
matrix[(x,y):(dx,dy):(sx, sy)]
return a sub matrix of dimension dx, dy starting from coordinate x, y by step of sx, sy. It does not look phenomenal, but, how do you express a subset of n nth dimensional array with the current [::] slice notation ?
Because, of course, reading the code I see a door for extrapolating a sane consistent behaviour with 2 arbitrary dimensions to arbitrary sized array.

After all, « natural » Nth dimension array makes sense like a ... video frambuffer wich is a matrix view in terms of pixel, but remember a pixel is an array of values (R,G,B, alpha).

With this notation getting the Blue followed by the alpha from a frambuffer would look like :
screen = matrix(map(int,open("/dev/fb0", "rb").read()), dim = (1980,1020, 4))

# get blue (rgba[3]) and alpha (4) at position 243, 743 would be written : 
blue, alpha = screen[(243,743, 3):(243, 743, 4) ]

You could also imagine with a 3D variant do « CPU intensive and slow » blitting (which does beat the purpose of a blitter (video instruction used to manipulate fast memory zone, including sprite handling).
Juste imagine we could implement an elliptic safe behaviour such as we could write :
screen[(dx, dy):(dx+width, dy+height)] = sprite([ [1, 2, ... ]) # sprite being width x height size
All the code being a view, we could actually memory map (mmap) my inner property __val__ to an actual raw framebuffer. which could be usefull for my hello world done by manipulating the raw framebuffer to actually print each letters :D

By writing about something I coded on a wit, without finality, I discover it may actually be usefull in stopping to waste gazillions of Joule dereferencing references of data scattered among memo provoking huge « page fault » storms. The actual solution is not perfect, because it eliminate only one level of dereferentiation, however by using the venerable array (to try), I expect significant gains without even thinking of importing numpy.

For now I have the following test passing regarding ACTUAL usage of the lib regarding addressing memory consistently and extracting patches of memory.
ex :
matrix([[1,2], [3,4]])[(0,0):(0,1)] == [[1,2]] # return first row
matrix([[1,2], [3,4]])[(0,0):(1,0)] == [[1],[3]] # return first column

m= matrix([[1,2],[3,4], [5,6], [7,8],[9,10]])

assert m[1][0] == m[(1,0)]
assert m[3][1] == m[(3,1)]
assert m[0] == [1,2]
assert m[(0,0)]==1
assert m[(0,1)]==2
assert m[(1,0)]==3
assert m[(1,1)]==4
assert m[(0,0):(0,1)]==[m[0],]
# cannot be expressed with regular slice 
assert m[(0,0):(1,0)] ==  [[1], [3]]
assert m[(0,0):(1,1)]==[[1,2], [3,4]]

assert m[(0,0):(2,1):(2,1)] == [[1, 2], [], [5, 6]]
assert m[(0,0):(3,0)] ==  [[1], [3], [5], [7]]
assert m[(0,1):(3,1)] == [[2], [4], [6], [8]]
assert m[(0,0):(3,0):(2,1)] == [[1], [], [5]]
assert m[(0,0):(3,0):(3,1)] == [[1], [], [], [7]]

This kind of notation for writing matrices makes writing matrix multiplication a little more readable :
def mul(m,n):
    mx,my=m.dim
    nx,ny=n.dim
    v = mdim([0,] * (mx * ny),  dim = (mx, ny))
    stride=my-1
    for x in range(mx):
        for y in range(ny):
            zz  = sum(map(lambda e: e[0]*e[1],zip(
                m[(x,0):(x,0+stride)][0],  # xth line from m
                mdim(n[(0,y):(0+stride,y)]).__val__, # yth row from n
            )))
            v[(x,y)]=zz
    return v
It is useless, but I'm pretty proud because, it really is useful, was coded in a day and passed most tests I could think of.

Annex

The code :
class DimensionError(Exception):
  pass

class matrix(object):
    def __eq__(self, tocomp):
        other=tocomp
        if type(tocomp) != type(self):
            other=matrix(tocomp)
        return (self.__val__, self.dim) == (other.__val__, other.dim)

    def __init__(self,a,**kw):
        if kw.get("dim"):
            self.__val__ = a
            self.dim = kw["dim"]
        else:
            self.__val__,self.dim = list(self.flattener(a))

    def flattener(self, v):
        dimy=len(v[0])
        dimx=len(v)
        ret = [ j for i in v for j in i ]
        assert len(ret) == dimx*dimy
        return  ret,(dimx, dimy)

    def unflattener(self):
        return self[(0,0): (self.dim[0]-1, self.dim[1]-1)]

    def __setitem__(self, k, v):
        if type(k)==int:
            # set nth line with slice
            raise NotImplementedError()
        for e, p in  enumerate(zip(k,self.dim)):
            i,j=p
            try:
                assert  i < = j, f"ooupla dim:{e+1} {i}>{j} "
            except:
                raise IndexError
        x,y=k
        dy=self.dim[1]
        offby=x*(dy)+y
        self.__val__[offby]=v

    def __repr__(self):
        return repr(self.unflattener())

    def __getitem__(self, k):
        x,y = 0,0
        dx, dy=self.dim
        step=(1,1)
        if type(k)==int:
			return self[(k,0):(k,dy-1)][0]

        if not isinstance(k,slice):
            x, y=k
            for e, p in  enumerate(zip(k,self.dim)):
                i,j=p
                try:
                    assert  i <= j, f"ooupla dim:{e+1} {i}<={j} k={k}"
                except:
                    raise IndexError
            off=x*dy+y
            return self.__val__[off]

        start, stop, step=k.start,k.stop,k.step
        step = step or (1,1)
        x0, xn, xd=start[1], stop[1]+1, step[1]
        if x0 > xn or xn < 0:
            raise NotImplementedError()
        y0, yn, yd=start[0], stop[0]+1, step[0]
        if y0 > yn or yn < 0:
            raise NotImplementedError()

        ret = []
        
        for y in range(y0,yn, yd):
            ret+=[[],]
            for x in range(x0, xn, xd):
                # row padding
                # I HATE PYTHON for [[],] * n => shallow copies !
                # if not for this bug, padding would be as simple as 
                # if y-y0>=len(ret):
                #    ret+=[[],] * (y-y0-len(ret)+1)
                for k in range(len(ret)-1, y-y0):
                    ret+=[[],]
                ret[y-y0]+=[self[(y,x)]]
        return ret

def transpose(md):
    res = []
    for x in range(md.dim[1]):
        res+=[[0] * md.dim[0],]
        for y in range(md.dim[0]):
            res[x][y]=md[(y,x)]
    return matrix(res)

def mul(m,n):
    mx,my=m.dim
    nx,ny=n.dim
    v = matrix([0,] * (mx * ny),  dim = (mx, ny))
    stride=my-1
    for x in range(mx):
        for y in range(ny):
            #z=sum([ x*y for x,y in zip(matrix(n[(0,y):(0+stride,y)]).__val__,matrix(m[(x,0):(x,0+stride)]).__val__)])
            zz  = sum(map(lambda e: e[0]*e[1],zip(
                m[(x,0):(x,0+stride)][0],  # xth line from m
                matrix(n[(0,y):(0+stride,y)]).__val__, # yth row from n
            )))
            v[(x,y)]=zz
    return v

m= matrix([[1,2],[3,4], [5,6], [7,8],[9,10]])

assert m[1][0] == m[(1,0)]
assert m[2][0] == m[(2,0)]
assert m[1][1] == m[(1,1)]
assert m[3][1] == m[(3,1)]
assert m[(-1,-1)]==8
assert m[0] == [1,2]
assert m[1] == [3,4]
assert m[2] == [5,6]
assert m[3] == [7,8]
assert m[(0,0)]==1
assert m[(0,1)]==2
assert m[(1,0)]==3
assert m[(1,1)]==4
assert m[(2,0)]==5
assert m[(2,1)]==6
assert m[(0,0):(0,1)]==[m[0],]
assert m[(0,0):(1,0)] ==  [[1], [3]]
assert m[(0,0):(1,1)]==[[1,2], [3,4]]

assert m[(0,0):(2,1):(2,1)] == [[1, 2], [], [5, 6]]
assert m[(0,0):(3,0)] ==  [[1], [3], [5], [7]]
assert m[(0,1):(3,1)] == [[2], [4], [6], [8]]
assert m[(0,0):(3,0):(2,1)] == [[1], [], [5]]
assert m[(0,0):(3,0):(3,1)] == [[1], [], [], [7]]
assert m[(0,0)] == 1
assert m[(1,0)]==3

try:
    print(m[(6,1)])
except IndexError:
    pass

from time import time
m= matrix([[1,2],[3,4], [5,6], [7,8],[9,10]])
n= matrix([[11,12],[13,14], [15,16], [17,18],[19,20]])
start=time()
for i in range(100_000):
    c=matrix(list(map(lambda x:x[0]+x[1], zip(m.__val__, n.__val__))), dim=m.dim)
print("addition with zip")
print(time() - start)
start=time()
for i in range(100_000):
    dx,dy=m.dim
    c=[]
    for x in range(dx):
        c+=[[0] * (dy), ]
        for y in range(dy):
            c[x][y]=m[x][y]+n[x][y]
print("addition with full loops")
print(time() - start)

m = matrix([[1,2,0],[4,3,-1]])
n = matrix([[5,1],[2,3],[3,4]])
assert mul(m,n)==[[9, 7], [23, 9]]

n,m = m,n
assert mul(m,n) == [[9, 13, -1], [14, 13, -3], [19, 18, -4]]
m=matrix([[1],[-1],[1], [1]])
n=matrix([[-10,2,3, 4],])
assert mul(m,n)==[[-10, 2, 3, 4], [10, -2, -3, -4], [-10, 2, 3, 4], [-10, 2, 3, 4]]
assert mul(n,m)== [[-5]]
m=matrix([[1],[-1],[1]])
n=matrix([[-7,2,3],[1,-2,3],])
assert mul(m,n)==[[-7, 2, 3], [7, -2, -3], [-7, 2, 3]]
assert mul(n,m) == [[-6], [6]]
m= matrix([[1,2],[3,4], [5,6], [7,8],[9,10]])
assert m[1] == [3,4]
assert mul(transpose(m),m)==[[165, 190], [190, 220]]
assert transpose(m) == [[1, 3, 5, 7, 9], [2, 4, 6, 8, 10]]
assert m[1] == [3,4]

No comments: