docs: def simulate_indices(self, ts_length, init=None, num_reps=None,
                           random_state=None):
          """
          Simulate time series of state transitions, where state indices
          are returned.
  
          Parameters
          ----------
          ts_length : scalar(int)
              Length of each simulation.
  
          init : int or array_like(int, ndim=1), optional
              Initial state(s). If None, the initial state is randomly
              drawn.
  
          num_reps : scalar(int), optional(default=None)
              Number of repetitions of simulation.
  
          random_state : int or np.random.RandomState, optional
              Random seed (integer) or np.random.RandomState instance to
              set the initial state of the random number generator for
              reproducibility. If None, a randomly initialized RandomState
              is used.
  
          Returns
          -------
          X : ndarray(ndim=1 or 2)
              Array containing the state values of the sample path(s). See
              the `simulate` method for more information.
  
          """
          random_state = check_random_state(random_state)
          dim = 1  # Dimension of the returned array: 1 or 2
  
          msg_out_of_range = 'index {init} is out of the state space'
  
          try:
              k = len(init)  # init is an array
              dim = 2
              init_states = np.asarray(init, dtype=int)
              # Check init_states are in the state space
              if (init_states >= self.n).any() or (init_states < -self.n).any():
                  idx = np.where(
                      (init_states >= self.n) + (init_states < -self.n)
                  )[0][0]
                  raise ValueError(msg_out_of_range.format(init=idx))
              if num_reps is not None:
                  k *= num_reps
                  init_states = np.tile(init_states, num_reps)
          except TypeError:  # init is a scalar(int) or None
              k = 1
              if num_reps is not None:
                  dim = 2
                  k = num_reps
              if init is None:
                  init_states = random_state.randint(self.n, size=k)
              elif isinstance(init, numbers.Integral):
                  # Check init is in the state space
                  if init >= self.n or init < -self.n:
                      raise ValueError(msg_out_of_range.format(init=init))
                  init_states = np.ones(k, dtype=int) * init
              else:
                  raise ValueError(
                      'init must be int, array_like of ints, or None'
                  )
  
          # === set up array to store output === #
          X = np.empty((k, ts_length), dtype=int)
  
          # Random values, uniformly sampled from [0, 1)
          random_values = random_state.random_sample(size=(k, ts_length-1))
  
          # Generate sample paths and store in X
          if not self.is_sparse:  # Dense
              _generate_sample_paths(
                  self.cdfs, init_states, random_values, out=X
              )
          else:  # Sparse
              _generate_sample_paths_sparse(
                  self.cdfs1d, self.P.indices, self.P.indptr, init_states,
                  random_values, out=X
              )
  
          if dim == 1:
              return X[0]
          else:
              return X