code: def _fill_diagonals(self, m):
          
          assert m.shape == (self.nsites, N_CODON, N_CODON)
          for r in range(self.nsites):
              scipy.fill_diagonal(m[r], 0)
              m[r][self._diag_indices] -= scipy.sum(m[r], axis=1)
docs: Fills diagonals of `nsites` matrices in `m` so rows sum to 0.