Reducing M in row_echelon() begins to produce NAN values after iterations.
Wait, I thought you were on course 2? 
Well, there was logic in the template code to handle the case that the pivot was zero, right? Did you rearrange any of that? Maybe it’s worth getting a clean copy and comparing that to your existing code.
To help in knowing what to look for, here are some experiments with Inf and NaN to show examples of the types of operations that can generate these values:
v = 42. * np.ones((1,4), dtype = 'float64')
print(f"type(v) = {type(v)}")
print(f"v = {v}")
w = np.zeros((1,4), dtype = 'float64')
print(f"type(w) = {type(w)}")
print(f"w = {w}")
z = v / w
print(f"z = {z}")
a = -1. * z
print(f"a = {a}")
b = z + 42.
print(f"b = {b}")
c = z - 42.
print(f"c = {c}")
d = z + z
print(f"d = {d}")
e = z - z
print(f"e = {e}")
f = z / z
print(f"f = {f}")
g = np.log(w)
print(f"g = {g}")
h = np.log(-1 * v)
print(f"h = {h}")
Running the above yields this output:
type(v) = <class 'numpy.ndarray'>
v = [[42. 42. 42. 42.]]
type(w) = <class 'numpy.ndarray'>
w = [[0. 0. 0. 0.]]
z = [[inf inf inf inf]]
a = [[-inf -inf -inf -inf]]
b = [[inf inf inf inf]]
c = [[inf inf inf inf]]
d = [[inf inf inf inf]]
e = [[nan nan nan nan]]
f = [[nan nan nan nan]]
g = [[-inf -inf -inf -inf]]
h = [[nan nan nan nan]]
Actually it occurred to me later that I left out the case that is perhaps the most likely in this instance:
j = w / w
print(f"j = {j}")
Which gives this result:
j = [[nan nan nan nan]]