docs: def sample(self, bqm, beta_range=None, num_reads=10, sweeps=1000,
                 beta_schedule_type="geometric", seed=None,
                 interrupt_function=None, initial_states=None):
          """Sample from a binary quadratic model using an implemented sample method.
  
          Args:
              bqm (:obj:`dimod.BinaryQuadraticModel`):
                  The binary quadratic model to be sampled.
  
              beta_range (tuple, optional):
                  A 2-tuple defining the beginning and end of the beta schedule, where beta is the
                  inverse temperature. The schedule is applied linearly in beta. Default range is set
                  based on the total bias associated with each node.
  
              num_reads (int, optional, default=10):
                  Each read is the result of a single run of the simulated annealing algorithm.
  
              sweeps (int, optional, default=1000):
                  Number of sweeps or steps.
  
              beta_schedule_type (string, optional, default='geometric'):
                  Beta schedule type, or how the beta values are interpolated between
                  the given 'beta_range'. Supported values are:
  
                  * linear
                  * geometric
  
              seed (int, optional):
                  Seed to use for the PRNG. Specifying a particular seed with a constant
                  set of parameters produces identical results. If not provided, a random seed
                  is chosen.
  
              initial_states (tuple(numpy.ndarray, dict), optional):
                  A tuple where the first value is a numpy array of initial states to seed the
                  simulated annealing runs, and the second is a dict defining a linear variable
                  labelling.
  
              interrupt_function (function, optional):
                  If provided, interrupt_function is called with no parameters between each sample of
                  simulated annealing. If the function returns True, then simulated annealing will
                  terminate and return with all of the samples and energies found so far.
  
          Returns:
              :obj:`dimod.Response`: A `dimod` :obj:`~dimod.Response` object.
  
          Examples:
              This example runs simulated annealing on a binary quadratic model with some
              different input parameters.
  
              >>> import dimod
              >>> import neal
              ...
              >>> sampler = neal.SimulatedAnnealingSampler()
              >>> bqm = dimod.BinaryQuadraticModel({'a': .5, 'b': -.5}, {('a', 'b'): -1}, 0.0, dimod.SPIN)
              >>> # Run with default parameters
              >>> response = sampler.sample(bqm)
              >>> # Run with specified parameters
              >>> response = sampler.sample(bqm, seed=1234, beta_range=[0.1, 4.2],
              ...                                num_reads=1, sweeps=20,
              ...                                beta_schedule_type='geometric')
              >>> # Reuse a seed
              >>> a1 = next((sampler.sample(bqm, seed=88)).samples())['a']
              >>> a2 = next((sampler.sample(bqm, seed=88)).samples())['a']
              >>> a1 == a2
              True
  
          """
  
          # if already index-labelled, just continue
          if all(v in bqm.linear for v in range(len(bqm))):
              _bqm = bqm
              use_label_map = False
          else:
              try:
                  inverse_mapping = dict(enumerate(sorted(bqm.linear)))
              except TypeError:
                  # in python3 unlike types cannot be sorted
                  inverse_mapping = dict(enumerate(bqm.linear))
              mapping = {v: i for i, v in iteritems(inverse_mapping)}
  
              _bqm = bqm.relabel_variables(mapping, inplace=False)
              use_label_map = True
  
          # beta_range, sweeps are handled by simulated_annealing
  
          if not isinstance(num_reads, Integral):
              raise TypeError("'samples' should be a positive integer")
          if num_reads < 1:
              raise ValueError("'samples' should be a positive integer")
  
          if not (seed is None or isinstance(seed, Integral)):
              raise TypeError("'seed' should be None or a positive integer")
          if isinstance(seed, Integral) and not (0 < seed < (2**64 - 1)):
              error_msg = "'seed' should be an integer between 0 and 2^64 - 1"
              raise ValueError(error_msg)
  
          if interrupt_function is None:
              def interrupt_function():
                  return False
  
          num_variables = len(_bqm)
  
          # get the Ising linear biases
          linear = _bqm.spin.linear
          h = [linear[v] for v in range(num_variables)]
  
          quadratic = _bqm.spin.quadratic
          if len(quadratic) > 0:
              couplers, coupler_weights = zip(*iteritems(quadratic))
              couplers = map(lambda c: (c[0], c[1]), couplers)
              coupler_starts, coupler_ends = zip(*couplers)
          else:
              coupler_starts, coupler_ends, coupler_weights = [], [], []
  
          if beta_range is None:
              beta_range = _default_ising_beta_range(linear, quadratic)
  
          sweeps_per_beta = max(1, sweeps // 1000.0)
          num_betas = int(math.ceil(sweeps / sweeps_per_beta))
          if beta_schedule_type == "linear":
              # interpolate a linear beta schedule
              beta_schedule = np.linspace(*beta_range, num=num_betas)
          elif beta_schedule_type == "geometric":
              # interpolate a geometric beta schedule
              beta_schedule = np.geomspace(*beta_range, num=num_betas)
          else:
              raise ValueError("Beta schedule type {} not implemented".format(beta_schedule_type))
  
          if seed is None:
              # pick a random seed
              seed = randint(0, (1 << 64 - 1))
  
          np_rand = np.random.RandomState(seed % 2**32)
  
          states_shape = (num_reads, num_variables)
  
          if initial_states is not None:
              initial_states_array, init_label_map = initial_states
              if not initial_states_array.shape == states_shape:
                  raise ValueError("`initial_states` must have shape "
                                   "{}".format(states_shape))
              if init_label_map is not None:
                  get_label = inverse_mapping.get if use_label_map else lambda i: i
                  initial_states_array = initial_states_array[:, [init_label_map[get_label(i)] for i in range(num_variables)]]
              numpy_initial_states = np.ascontiguousarray(initial_states_array, dtype=np.int8)
          else:
              numpy_initial_states = 2*np_rand.randint(2, size=(num_reads, num_variables)).astype(np.int8) - 1
  
          # run the simulated annealing algorithm
          samples, energies = sa.simulated_annealing(num_reads, h,
                                                     coupler_starts, coupler_ends,
                                                     coupler_weights,
                                                     sweeps_per_beta, beta_schedule,
                                                     seed,
                                                     numpy_initial_states,
                                                     interrupt_function)
          off = _bqm.spin.offset
          info = {
              "beta_range": beta_range,
              "beta_schedule_type": beta_schedule_type
          }
          response = dimod.SampleSet.from_samples(
              samples,
              energy=energies+off,
              info=info,
              vartype=dimod.SPIN
          )
  
          response.change_vartype(_bqm.vartype, inplace=True)
          if use_label_map:
              response.relabel_variables(inverse_mapping, inplace=True)
  
          return response