Inconsistent solutions when computing the stationary distribution for small Markov chain

Question

I have a small discrete Markov chain with the following transition matrix:

$$ P = \begin{bmatrix} 0 & 0 & 1-r & r \\ 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 \\ r & 1-r & 0 & 0 \end{bmatrix} $$

This should be ergodic, since the whole Markov chain is a single recurrent class. However, when I try to solve for the stationary distribution $\pi$ in two ways, I get inconsistent solutions.

Since $\pi^{\top} = \pi^{\top} P$ by definition, one way to brute force approximate the stationary distribution is just to raise $P$ to a high power.

r = 1/7
P = np.array([
    [0, 0  , 1-r, r],
    [0, 0  , 0, 1],
    [1, 0  , 0, 0],
    [r, 1-r, 0, 0]
])
assert np.allclose(P.sum(axis=1), 1)

Pn = np.linalg.matrix_power(P, 1000)
pi = Pn[0, :]
assert np.isclose(pi.sum(), 1)
print(pi.squeeze())
print(np.dot(pi, P).squeeze())

This prints something that is clearly not the stationary distribution, since $\pi^{\top} \neq \pi^{\top} P$:

[0.53846154 0.46153846 0.         0.        ]
[0.         0.         0.46153846 0.53846154]

However, if I solve the system of linear equations, I get a different (also non-stationary) answer:

A = np.array([
    [-1 , 0  , 1-r, r  ],
    [0  , -1 , 0  , 1  ],
    [1  , 0  , -1 , 0  ],
    [r  , 1-r, 0  , -1 ],
    [1  , 1  , 1  , 1  ]
])
b  = np.array([0, 0, 0, 0, 1])[:, None]
pi = np.linalg.lstsq(A, b, rcond=None)[0]
pi = pi.reshape(1, -1)
assert np.isclose(pi.sum(), 1)

print(pi.squeeze())
print(np.dot(pi, P).squeeze())

Here, this prints

[0.25 0.25 0.25 0.25]
[0.28571429 0.21428571 0.21428571 0.28571429]

What am I doing wrong when trying to solve for $\pi$?