def METHOD_NAME(H0, W0):
  """
  Return matrices mapping from latent statespace to input space and output
  space. These are the matrices H and W, respectively, used in the sensorization
  nowcast_fusion kernel. Since some outputs may be indeterminate, the indices of the
  fully determined rows are returned. This may be used, for example, to find
  the set of outputs which make up the rows of the returned W matrix.
  inputs:
    H0: map from full statespace to inputs (I x S)
    W0: map from full statespace to outputs (O x S)
  outputs:
    - the matrix H, mapping subspace to inputs (I x S')
    - the matrix W, mapping subspace to outputs (O' x S')
    - list of row indices of W0 that make up W (O')
  notes:
    - S' <= S and O' <= O
    - for numerical stability, inputs should be matrices of Fractions
  """
  # helper function to convert a float matrix into a fraction matrix
  fractions = lambda X: np.array([[Fraction(x) for x in row] for row in X])
  # Find a set of basis vectors that span the same subspace (of the full
  # statespace) that is spanned by the input vectors in H0. The result is a
  # minimal set of elements from which all inputs can be unambiguously
  # determined.
  B = eliminate(H0.copy())
  # the dimensions of full statespace (number of columns)
  size = B.shape[1]
  # the dimensions of the subspace (number of non-empty rows)
  rank = np.sum(np.sum(np.abs(B), axis=1) > 0)
  # B should be a square matrix with rows of zeros below rows of basis vectors
  num_rows = B.shape[0]
  if num_rows < size:
    Z = fractions(np.zeros((size - num_rows, size)))
    B = np.vstack((B, Z))
  elif num_rows > size:
    B = B[:size, :]
  # Attempt to build each input and output vector as a linear combination of
  # the subspace basis vectors. Since B may not be full rank, it may not be
  # invertible. Instead, solve by eliminating the augmented matrix of B
  # (transposed) with the identity matrix. After elimination, the (transposed)
  # inverse of B is contained within the augmented matrix.
  I = fractions(np.eye(size))
  BtI = np.hstack((B.T, I))
  IBit = eliminate(BtI)
  Bi = IBit[:, size:].T
  # possible, or "actual", solutions are in the leftmost columns
  # impossible, or "pseudo", solutions are in the rightmost columns
  Bi_actual, Bi_pseudo = Bi[:, :rank], Bi[:, rank:]
  # compute H, the map from statespace B to inputs
  # all inputs are within the span of statespace B
  H = matmul(H0, Bi_actual)
  # compute W, the map from statespace B to outputs
  # outputs not within the span of statespace B must be excluded
  W_actual = matmul(W0, Bi_actual)
  W_pseudo = matmul(W0, Bi_pseudo)
  # only keep rows where the coeficient of all pseudo basis vectors is zero
  actual_rows = np.flatnonzero(np.sum(np.abs(W_pseudo), axis=1) == 0)
  W = W_actual[actual_rows, :]
  # return H, W, and the indices of the rows of W0 that make up W
  return H, W, actual_rows