Skip to content

Instantly share code, notes, and snippets.

Per Vognsen pervognsen

Block or report user

Report or block pervognsen

Hide content and notifications from this user.

Learn more about blocking users

Contact Support about this user’s behavior.

Learn more about reporting abuse

Report abuse
View GitHub Profile
View random.py
def random_sphere_point(n):
# Any isotropic distribution would do, but isotropic multivariate normals are just a product of univariate normals.
x = np.random.randn(n)
return x / np.linalg.norm(x)
def random_ball_point(n):
# The volume cdf scales as r^n, hence the inverse cdf scales as r^(1/n).
return random_sphere_point(n) * np.random.random()**(1/n)
View nice.py
A = blocks(A, (w, h))
B = blocks(B, (h, w))
C = blocks(out, (w, w))
for j in range(B.shape[0]):
for k in range(B.shape[1]):
np.copyto(column_panel, B[j, k])
for i in range(A.shape[0]):
np.copyto(row_panel, A[i, j])
microkernel(C[i, k], row_panel, column_panel)
View blocked_matrix_multiply.py
# Version 0: Use slice views to access a matrix like a block matrix.
def blocked_matrix_multiply0(A, B, block_size, panel_size):
m, p, n = A.shape[0], A.shape[1], B.shape[1]
w, h = block_size, panel_size
C = np.zeros((m, n))
for i in range(0, m, w):
for j in range(0, p, h):
row_panel = A[i:i+w, j:j+h]
for k in range(0, n, w):
column_panel = B[j:j+h, k:k+w]
View reshape.py
A = np.arange(8 * 4).reshape(8, 4)
print(A)
# [[ 0 1 2 3]
# [ 4 5 6 7]
# [ 8 9 10 11]
# [12 13 14 15]
# [16 17 18 19]
# [20 21 22 23]
# [24 25 26 27]
# [28 29 30 31]]
View linalg.py
def plu_inplace(A):
P = np.arange(len(A))
for i in range(len(A)):
j = i + np.argmax(np.abs(A[i:, i]))
P[[i, j]], A[[i, j]] = P[[j, i]], A[[j, i]]
A[i+1:, i] /= A[i, i]
A[i+1:, i+1:] -= np.outer(A[i+1:, i], A[i, i+1:])
return P, A, A
def l1_solve(L, y):
View tricks.py
# Tensor product contractions with Einstein's summation notation. Examples:
# Matrix-matrix multiply is tensor(A, 'ij', B, 'jk')
# Matrix-vector multiply is tensor(A, 'ij', v, 'j')
# Matrix trace is tensor(A, 'ii')
# Matrix transpose is tensor(A, 'ji')
# Matrix diagonal is tensor('i', A, 'ii')
# Inner product is tensor(v, 'i', w, 'i')
# Outer product is tensor(v, 'i', w, 'j')
def tensor(*args, **kwargs):
result = ''
View multidimensional_array_views.md

Multi-dimensional array views for systems programmers

As C programmers, most of us think of pointer arithmetic for multi-dimensional arrays in a nested way:

The address for a 1-dimensional array is base + x. The address for a 2-dimensional array is base + x + y*x_size for row-major layout and base + y + x*y_size for column-major layout. The address for a 3-dimensional array is base + x + (y + z*y_size)*x_size for row-column-major layout. And so on.

View key_zip.py
def cmp(x, y):
return (x > y) - (x < y)
nil = object()
def sorted_zip(xs, ys, key=lambda x: x, cmp=cmp, nil=nil):
xs = iter(xs)
ys = iter(ys)
x = next(xs, nil)
y = next(ys, nil)
View low_variance_sampling.py
class DiscreteDistribution:
def __init__(self):
self.total_weight = 0
self.points = {}
def add(self, point, weight):
assert point not in self.points
self.points[point] = weight
self.total_weight += weight
View rad.py
# Reverse-mode automatic differentiation
import math
# d(-x) = -dx
def func_neg(x):
return -x, [-1]
# d(x + y) = dx + dy
def func_add(x, y):
You can’t perform that action at this time.