docs: def run_tpm(system, steps, blackbox):
      """Iterate the TPM for the given number of timesteps.
  
      Returns:
          np.ndarray: tpm * (noise_tpm^(t-1))
      """
      # Generate noised TPM
      # Noise the connections from every output element to elements in other
      # boxes.
      node_tpms = []
      for node in system.nodes:
          node_tpm = node.tpm_on
          for input_node in node.inputs:
              if not blackbox.in_same_box(node.index, input_node):
                  if input_node in blackbox.output_indices:
                      node_tpm = marginalize_out([input_node], node_tpm)
  
          node_tpms.append(node_tpm)
  
      noised_tpm = rebuild_system_tpm(node_tpms)
      noised_tpm = convert.state_by_node2state_by_state(noised_tpm)
  
      tpm = convert.state_by_node2state_by_state(system.tpm)
  
      # Muliply by noise
      tpm = np.dot(tpm, np.linalg.matrix_power(noised_tpm, steps - 1))
  
      return convert.state_by_state2state_by_node(tpm)