Document that masked arrays require special care

Masked values propagate, which sometimes results in a counter-intuitive behavior of our return values. We should document this.

Original description

This minimal example works just fine:

# Create matrices
N = 10
H_0_AA = np.diag(np.arange(N, dtype=float))
H_0_BB = 10 * np.eye(N) + np.diag(np.arange(N))
h_p_AB = np.random.rand(N, N)
H_p_AB = {(1,): h_p_AB}
H_p_BA = {(1,): h_p_AB.T}
h_p_AA = np.random.rand(N, N)
H_p_AA = {(1,): 0.5 * (h_p_AA + h_p_AA.T)}
h_p_BB = np.random.rand(N, N)
H_p_BB = {(1,): 0.5 * (h_p_BB + h_p_BB.T)}

H = to_BlockSeries(H_0_AA, H_0_BB, H_p_AA, H_p_BB, H_p_AB)
H_tilde, U, U_adjoint = general(H)

H_AA = sum(H_tilde.evaluated[0, 0, :3])
H_AA

However, if we don't pass H_p_AA, the resulting array H_AA is empty. If we don't pass H_p_BB, H_BB = sum(H_tilde.evaluated[1, 1, :3]) is empty. Here's a minimal code the reproduces the issue:

# Create matrices
N = 10
H_0_AA = np.diag(np.arange(N, dtype=float))
H_0_BB = 10 * np.eye(N) + np.diag(np.arange(N))
h_p_AB = np.random.rand(N, N)
H_p_AB = {(1,): h_p_AB}
H_p_BA = {(1,): h_p_AB.T}
H_p_AA = {}
h_p_BB = np.random.rand(N, N)
H_p_BB = {(1,): 0.5 * (h_p_BB + h_p_BB.T)}

H = to_BlockSeries(H_0_AA, H_0_BB, H_p_AA, H_p_BB, H_p_AB)
H_tilde, U, U_adjoint = general(H)

H_AA = sum(H_tilde.evaluated[0, 0, :3])
H_AA