Plotting determinant probabilities

Question

This is a fun exercise that tries to answer the following question: if the elements from a square matrix are randomly chosen from 0 to 9, is it more like that the determinant would be even or odd?

I'm pretty sure there are different ways of optimizing this piece of code. The first one that comes to my mind is to use only generators (maybe?) and not use the for loop to fill values.

import itertools
import math

import matplotlib.pyplot as plt
import numpy as np

n = 4
lst = list(itertools.product([0, 1], repeat=n))

algo = list(itertools.product(lst, repeat=n))
determinants = np.abs(np.linalg.det(algo))

determinants = np.where(determinants % 2 == 0, 0, 1)

square = math.ceil((math.sqrt(len(determinants)))) ** 2

# fill with generic value if not a perfect square
for num in range(square - len(determinants)):
    determinants = np.append(determinants, -1) 


matrix = np.reshape(determinants, (int(math.sqrt(square)), int(math.sqrt(square))))


fig, ax = plt.subplots()  # or
ax.imshow(matrix, cmap="Greens")
plt.show()

All suggestions are welcome :D

Answer 1

The n you care about isn't really 4, it's 4**2 = 16.

Don't use lists, itertools or for loops in a Numpy context.

Don't use np.where in this case; you can replace that entire operation with % 2 alone.

Don't use math.* methods when there are Numpy equivalents.

Don't concatenate-append. Instead, just use slices on a pre-allocated array.

When you reshape to your matrix, leave one of the dimensions implied as -1.

Suggested

The fill step could look like

side = round(np.ceil(np.sqrt(len(determinants))))
filled = np.empty(side**2)
filled[:len(determinants)] = determinants
filled[len(determinants):] = -1

but your operation will always be a perfect square, so just delete that, leaving this code which is equivalent to your original:

import matplotlib.pyplot as plt
import numpy as np

n = 16

lst = np.arange(2**n)
as_bytes = np.stack(
    (np.right_shift(lst, n//2), lst),
    axis=-1,
).astype(np.uint8)
algo = np.unpackbits(
    as_bytes, axis=1,
).reshape((-1, 4, n//4))

determinants = np.abs(np.linalg.det(algo)) % 2
side = round(np.ceil(np.sqrt(len(determinants))))
matrix = determinants.reshape((side, -1))

fig, ax = plt.subplots()
ax.imshow(matrix, cmap="Greens")
plt.show()

Output