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):
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]
