docs: def kron_with_controls(*matrices: np.ndarray) -> np.ndarray:
      """Computes the kronecker product of a sequence of matrices and controls.
  
      Use linalg.CONTROL_TAG to represent controls. Any entry of the output
      matrix corresponding to a situation where the control is not satisfied will
      be overwritten by identity matrix elements.
  
      The control logic works by imbuing NaN with the meaning "failed to meet one
      or more controls". The normal kronecker product then spreads the per-item
      NaNs to all the entries in the product that need to be replaced by identity
      matrix elements. This method rewrites those NaNs. Thus CONTROL_TAG can be
      the matrix [[NaN, 0], [0, 1]] or equivalently [[NaN, NaN], [NaN, 1]].
  
      Because this method re-interprets NaNs as control-failed elements, it won't
      propagate error-indicating NaNs from its input to its output in the way
      you'd otherwise expect.
  
      Args:
          *matrices: The matrices and controls to combine with the kronecker
              product.
  
      Returns:
          The resulting matrix.
      """
  
      product = kron(*matrices)
  
      # The NaN from CONTROL_TAG spreads to everywhere identity belongs.
      for i in range(product.shape[0]):
          for j in range(product.shape[1]):
              if np.isnan(product[i, j]):
                  product[i, j] = 1 if i == j else 0
  
      return product