Repository: dwavesystems/dwave-system
Pull request #562 (2025-03-24 22:19:08)

Add multiparameter pseudolikelihood estimation, include background susceptibility.
Replaces this pull request: https://github.com/dwavesystems/dwave-system/pull/552 which can be cancelled.

1) Add support for multiparameter pseudolikelihood estimation.
2) Add support for background susceptibility estimationr
3) Add type hinting
4) Correct PEP8 compliance issues (black linting)


Patch:
diff --git a/docs/conf.py b/docs/conf.py
index 163dcd8..4853867 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -50,7 +50,8 @@ exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store', 'sdk_index.rst']
 linkcheck_retries = 2
 linkcheck_anchors = False
 linkcheck_ignore = [r'https://cloud.dwavesys.com/leap',  # redirects, many checks
-                    r'https://www.jstor.org/stable']
+                    r'https://www.jstor.org/stable',
+                    r'https://doi.org/']
 
 pygments_style = 'sphinx'
 
@@ -81,4 +82,4 @@ html_sidebars = {"**": ["search-field", "sidebar-nav-bs"]}  # remove ads
 intersphinx_mapping = {'python': ('https://docs.python.org/3', None),
     'networkx': ('https://networkx.github.io/documentation/stable/', None),
     'dwave': ('https://docs.dwavequantum.com/en/latest/', None),
-    }
\ No newline at end of file
+    }
diff --git a/docs/utilities.rst b/docs/utilities.rst
index 5945513..5f97e40 100644
--- a/docs/utilities.rst
+++ b/docs/utilities.rst
@@ -37,4 +37,7 @@ Temperature and Unit-Conversion Utilities
    freezeout_effective_temperature
    h_to_fluxbias
    Ip_in_units_of_B
+   maximum_pseudolikelihood
    maximum_pseudolikelihood_temperature
+   background_susceptibility_ising
+   background_susceptibility_bqm
diff --git a/dwave/system/temperatures.py b/dwave/system/temperatures.py
index 1afab72..3a53dec 100644
--- a/dwave/system/temperatures.py
+++ b/dwave/system/temperatures.py
@@ -12,90 +12,126 @@
 #    See the License for the specific language governing permissions and
 #    limitations under the License.
 
-r"""The following effective temperature and bias estimators are provided:
+r"""The following parameter estimation methods are provided:
 
-- Maximum pseudo-likelihood is an efficient estimator for the temperature 
-  describing a classical Boltzmann distribution P(x) = \exp(-H(x)/T)/Z(T) 
-  given samples from that distribution, where H(x) is the classical energy 
+
+
+- Maximum pseudo-likelihood is an efficient estimator for the temperature
+  describing a classical Boltzmann distribution :math:`P(s) = \exp(-H(s)/T)/Z(T)`
+  given samples from that distribution, where :math:`H(x)` is the classical energy
   function. The following links describe features of the estimator in
-  application to equilibrium distribution drawn from binary quadratic models and
-  non-equilibrium distributions generated by annealing:
+  application to equilibrium distribution drawn from binary quadratic models
+  and non-equilibrium distributions generated by annealing:
   https://www.jstor.org/stable/25464568
-  https://doi.org/10.3389/fict.2016.00023 
+  https://doi.org/10.3389/fict.2016.00023
 
-- An effective temperature can be inferred assuming freeze-out during the 
-  anneal at s=t/t_a, an annealing schedule, and a device physical temperature.
+- An effective temperature can also be inferred assuming freeze-out during the
+  anneal at :math:`s=t/t_a`, an annealing schedule, and a device physical temperature.
   Necessary device-specific properties are published for online solvers in the
   :ref:`qpu_solver_properties_specific` section.
 
+- Maximum pseudo-likelihood can be used to infer multiple parameters of
+  an exponential (Boltzmann) distribution given fair samples. Code supports
+  inference for the probability distributions structured as
+  :math:`P(s) = \exp(-\text{sum}_i x_i H_i(s)) /Z(x)`, a generalization of a Boltzmann
+  distribution parameterized only by the temperature (:math:`T=1/x` for one :math:`H`).
+  :math:`H_i` are binary quadratic models defined on a common set of variables.
+  E.g.: To infer h1, h2 and J12 for a 2 qubit spin system at temperature 1 consider:
+  ``H_1 = dimod.BinaryQuadraticModel.from_ising({0:1, 0: 0}, {})``
+  ``H_2 = dimod.BinaryQuadraticModel.from_ising({0:0, 0: 1}, {})``
+  ``H_3 = dimod.BinaryQuadraticModel.from_ising({}, {(0,1): 1})``
+  This code is not optimized for inferring many local terms, but rather
+  parameters controlling the energy of many degrees of freedom (spins).
+  A practically important example is the case that :math:`H_1` is the
+  programmed Hamiltonian (:math:`x_1` is inverse temperature, beta) and
+  :math:`H_2` is a correction accounting for imbalance in h/J asymmetry or
+  background susceptibility.
+
 - The biases (h) equivalent to application of flux bias, or vice-versa,
-  can be inferred as a function of the anneal progress s=t/t_a by
+  can be inferred as a function of the anneal progress :math:`s=t/t_a` by
   device-specific unit conversion. The necessary parameters for estimation
   [Mafm, B(s)] are published for online solvers in the
   :ref:`qpu_solver_properties_specific` section
-""" 
+"""
 
 import warnings
 import numpy as np
+import networkx as nx
+
 import dimod
 from scipy import optimize
-from typing import Tuple, Union, Optional, Literal
-
-__all__ = ['effective_field', 'maximum_pseudolikelihood_temperature',
-           'freezeout_effective_temperature', 'fast_effective_temperature',
-           'Ip_in_units_of_B', 'h_to_fluxbias', 'fluxbias_to_h']
-
-def effective_field(bqm,
-                    samples=None,
-                    current_state_energy=False) -> (np.ndarray,list):
-    r'''Returns the effective field for all variables and all samples.
+from typing import Tuple, Union, Optional, Literal, List
+from collections import defaultdict
+
+__all__ = [
+    "effective_field",
+    "maximum_pseudolikelihood",
+    "maximum_pseudolikelihood_temperature",
+    "freezeout_effective_temperature",
+    "fast_effective_temperature",
+    "Ip_in_units_of_B",
+    "h_to_fluxbias",
+    "fluxbias_to_h",
+    "background_susceptibility_ising",
+    "background_susceptibility_bqm",
+]
+
+
+def effective_field(
+    bqm: dimod.BinaryQuadraticModel,
+    samples: Union[None, dimod.SampleSet, Tuple[np.ndarray, list]] = None,
+    current_state_energy: bool = False,
+) -> Tuple[np.ndarray, list]:
+    r"""Returns the effective field for all variables and all samples.
 
     The effective field with ``current_state_energy = False`` is the energy
-    attributable to setting a variable to value 1, conditioned on fixed values 
-    for all neighboring variables (relative to exclusion of the variable, and 
-    associated energy terms, from the problem). 
+    attributable to setting a variable to value 1, conditioned on fixed values
+    for all neighboring variables (relative to exclusion of the variable, and
+    associated energy terms, from the problem).
 
     The effective field with ``current_state_energy = True`` is the energy gained
-    by flipping the variable state against its current value (from say -1 to 1 
-    in the Ising case, or 0 to 1 in the QUBO case). A positive value indicates 
-    that the energy can be decreased by flipping the variable, hence the 
+    by flipping the variable state against its current value (from say -1 to 1
+    in the Ising case, or 0 to 1 in the QUBO case). A positive value indicates
+    that the energy can be decreased by flipping the variable, hence the
     variable is in a locally excited state.
-    If all values are negative  (positive) within a sample, that sample is a 
-    local minima (maxima).  
+    If all values are negative (positive) within a sample, that sample is a
+    local minima (maxima).
 
     Any BQM can be converted to an Ising model with
 
     .. math::
-        H(s) = Constant + \sum_i h_i s_i + 0.5 \sum_{i,j} J_{i,j} s_i s_j
+        H(s) = Constant + \sum_i h_i s_i + 0.5 \sum_{i, j} J_{i, j} s_i s_j
 
     with unique values of :math:`J` (symmetric) and :math:`h`. The sample
     dependent effect field on variable i, :math:`f_i(s)`, is then defined
 
     if current_state_energy == False:
         .. math::
-            f_i(s) = h_i + \sum_j J_{i,j} s_j
+            f_i(s) = h_i + \sum_j J_{i, j} s_j
     else:
-        .. math:: 
-            f_i(s) = 2 s_i [h_i + \sum_j J_{i,j}  s_j]
+        .. math::
+            f_i(s) = 2 s_i [h_i + \sum_j J_{i, j}  s_j]
 
     Args:
-        bqm (:obj:`dimod.BinaryQuadraticModel`): 
-            Binary quadratic model.
-        samples (samples_like or :obj:`~dimod.SampleSet`,optional):
-            A collection of raw samples. `samples_like` is an extension of
-            NumPy's array like structure. See :func:`dimod.sampleset.as_samples`.
-            By default, a single sample with all +1 assignments is used. 
-        current_state_energy (bool, optional, default=False): 
-            By default, returns the effective field (the energy 
-            contribution associated to a state assignment of 1). When set to 
-            True, returns the energy lost in flipping the value of each 
-            variable. Note current_state_energy is typically negative for 
+        bqm:
+            A dimod binary quadratic model.
+        samples:
+            Either a dimod SampleSet or a `samples_like` object, the later is
+            an extension of NumPy's array like structure.
+            See :func:`dimod.sampleset.as_samples`.
+            By default, a single sample with all +1 assignments is used.
+        current_state_energy:
+            By default, returns the effective field (the energy
+            contribution associated to a state assignment of 1). When set to
+            True, returns the energy lost in flipping the value of each
+            variable. Note current_state_energy is typically negative for
             positive temperature samples, meaning energy is not decreased
             by flipping the spin against its current assignment.
     Returns:
-        samples_like: 
+        samples_like:
             A Tuple of the effective_fields, and the variable labels.
-            Effective fields are returned as a :obj:`numpy.ndarray`.
+            Effective fields are returned as a :obj:`numpy.ndarray`,
+            variable labels are returned as a list.
             Rows index samples, and columns index variables in the order
             returned by variable labels.
 
@@ -109,137 +145,246 @@ def effective_field(bqm,
        >>> import numpy as np
        >>> from dwave.system.temperatures import effective_field
        >>> N = 5
-       >>> bqm = dimod.BinaryQuadraticModel.from_ising({}, {(i,i+1) : -0.5 for i in range(N-1)})
+       >>> bqm = dimod.BinaryQuadraticModel.from_ising({}, {(i, i+1) : -0.5 for i in range(N-1)})
        >>> var_labels = list(range(N))
        >>> samples = (np.ones(shape=(1,N)), var_labels)
-       >>> E = effective_field(bqm,samples,current_state_energy=True)
-       >>> print('Cost to flip spin against current assignment', E)     
+       >>> E = effective_field(bqm, samples, current_state_energy=True)
+       >>> print('Cost to flip spin against current assignment', E)
        Cost to flip spin against current assignment (array([[-1., -2., -2., -2., -1.]]), [0, 1, 2, 3, 4])
 
-    '''
+    """
     if samples is None:
-        samples = np.ones(shape=(1,bqm.num_variables))
+        samples = np.ones(shape=(1, bqm.num_variables))
         labels = bqm.variables
     else:
         samples, labels = dimod.sampleset.as_samples(samples)
 
     if bqm.vartype is dimod.BINARY:
-        bqm = bqm.change_vartype('SPIN', inplace=False)
-        samples = 2*samples - 1
+        bqm = bqm.change_vartype(dimod.SPIN, inplace=False)
+        samples = 2 * samples - 1
 
-    h, (irow, icol, qdata), offset = bqm.to_numpy_vectors(
-        variable_order=labels)
+    h, (irow, icol, qdata), offset = bqm.to_numpy_vectors(variable_order=labels)
     # eff_field = h + J*s OR diag(Q) + (Q-diag(Q))*b
-    effective_fields = np.tile(h[np.newaxis,:],(samples.shape[0],1))
+    effective_fields = np.tile(h[np.newaxis, :], (samples.shape[0], 1))
     for sI in range(samples.shape[0]):
-        np.add.at(effective_fields[sI,:],irow,qdata*samples[sI,icol])
-        np.add.at(effective_fields[sI,:],icol,qdata*samples[sI,irow])
+        np.add.at(effective_fields[sI, :], irow, qdata * samples[sI, icol])
+        np.add.at(effective_fields[sI, :], icol, qdata * samples[sI, irow])
 
     if current_state_energy is True:
-        #Ising: eff_field = 2*s*(h + J*s)
-        effective_fields = 2*samples*effective_fields
-
-    return (effective_fields,labels)
-
-def maximum_pseudolikelihood_temperature(bqm=None,
-                                         sampleset=None,
-                                         site_energy=None,
-                                         num_bootstrap_samples=0,
-                                         seed=None,
-                                         T_guess=None,
-                                         optimize_method='bisect',
-                                         T_bracket=(1e-3,1000)) -> Tuple[float,np.ndarray]:
-    r'''Returns a sampling-based temperature estimate.
-
-    The temperature T parameterizes the Boltzmann distribution as 
-    :math:`P(x) = \exp(-H(x)/T)/Z(T)`, where :math:`P(x)` is a probability over a state space, 
-    :math:`H(x)` is the energy function (BQM) and :math:`Z(T)` is a normalization. 
-    Given a sample set (:math:`S`), a temperature estimate establishes the 
+        effective_fields = 2 * samples * effective_fields
+
+    return (effective_fields, labels)
+
+
+def background_susceptibility_ising(
+    h: Union[np.ndarray, dict], J: Union[np.ndarray, dict]
+) -> Tuple:
+    """Create the field and couplings for the background susceptibility correction.
+
+    Background susceptibility is a significant source of systematic error
+    in annealing processors, it can be treated as a perturbation of the
+    programmed Hamiltonian.
+
+    The perturbed Hamiltonian can be most concisely expressed for an Ising
+    model in matrix vector notation.
+    Assuming J is a real symmetric (0 on diagonal) matrix, h and s are column
+    vectors, and chi is a small real value, then the Ising model is
+    is defined: :math:`H = k + (h + \chi J h)' s +  1/2 s' (J + \chi J^2) s`
+    ' denotes transpose and matrix multiplication applies.
+    The constant (k, irrelevant to sampled distributions) is defined
+    :math:`k = - \chi/2 Trace[J^2]`
+    This function returns the perturbative part
+    :math:`H(s) = k + h' J s +  1/2 s' (J + \chi J^2) s`
+
+    Args:
+       h: linear terms (fields) in the Hamiltonian.
+       J: quadratic terms (couplings) in the Hamiltonians
+
+    Returns:
+       A Tuple of fields, couplings and scalar constant. If `h` and `J` are of type
+       np.ndarray so are returned fields and couplings. Otherwise a tuple
+       of dictionaries is returned.
+
+    .. _Background susceptibility correction:
+    https://docs.dwavesys.com/docs/latest/c_qpu_ice.html#qpu-ice-background-susceptibility (accessed October 21, 2024)
+    """
+    if isinstance(h, np.ndarray) and isinstance(J, np.ndarray):
+        dh = J @ h  # matched dimensions assumed
+        dJ = J @ J  # J = J.T, a square matrix is assumed.
+        k = np.sum(np.diagonal(dJ)) / 2
+    else:
+        # Assume h and J have dictionary attributes
+        if len(h):
+            dh = {n: 0 for n in h.keys()}
+            for ij, Jval in J.items():
+                i, j = ij
+                dh[j] += h[i] * Jval
+                dh[i] += h[j] * Jval
+        else:
+            dh = {}
+
+        G = nx.from_edgelist(J.keys())
+        J = {
+            frozenset(e): Jval for e, Jval in J.items()
+        }  # Convert for unambiguous indexing
+        dJ = defaultdict(float)
+        for n in G.nodes():
+            neighs = list(G.neighbors(n))
+            for idx, n1 in enumerate(neighs):
+                for n2 in neighs[idx + 1 :]:
+                    dJ[tuple(frozenset((n1, n2)))] += (
+                        J[frozenset((n, n1))] * J[frozenset((n, n2))]
+                    )
+        dJ = dict(dJ)
+        k = 0
+
+    return dh, dJ, k
+
+
+def background_susceptibility_bqm(bqm: dimod.BinaryQuadraticModel, chi: Optional[float] = None):
+    """Create the binary quadratic model for the background susceptibility correction.
+
+    Background susceptibility is a perturbative correction to the programmed Hamiltonian:
+    bqm = bqm (no background susceptibility) + chi * dbqm (corrections)
+    If chi is None, dbqm is returned, otherwise the perturbed Hamiltonian is
+    returned.
+
+    Args:
+       bqm: A dimod binary quadratic model, describing the programmed Hamiltonian
+       chi: Scale of background susceptibility correction
+
+    Returns:
+       dimod.BinaryQuadraticModel: A dimod binary quadratic model decribing the
+       background susceptibility correction.
+
+    .. _Background susceptibility correction:
+    https://docs.dwavesys.com/docs/latest/c_qpu_ice.html#qpu-ice-background-susceptibility (accessed October 21, 2024)
+
+    """
+    source_type = bqm.vartype
+    bqm = bqm.change_vartype(dimod.SPIN)
+    dh, dJ, _ = background_susceptibility_ising(bqm.linear, bqm.quadratic)
+    dbqm = dimod.BinaryQuadraticModel(source_type).from_ising(dh, dJ)
+    if chi is not None:
+        dbqm = bqm.change_vartype(source_type) + chi * dbqm
+    return dbqm
+
+
+def maximum_pseudolikelihood_temperature(
+    bqm: Union[None, dimod.BinaryQuadraticModel] = None,
+    sampleset: Union[None, dimod.SampleSet, Tuple[np.ndarray, List]] = None,
+    en1: Optional[np.ndarray] = None,
+    num_bootstrap_samples: int = 0,
+    seed: Optional[int] = None,
+    T_guess: Optional[float] = None,
+    optimize_method: Optional[str] = None,
+    T_bracket: Tuple[float, float] = (1e-3, 1000),
+    sample_weights: Optional[np.ndarray] = None,
+) -> Tuple[float, np.ndarray]:
+    r"""Returns a sampling-based temperature estimate.
+
+    The temperature T parameterizes the Boltzmann distribution as
+    :math:`P(x) = \exp(-H(x)/T)/Z(T)`, where :math:`P(x)` is a probability over a state space,
+    :math:`H(x)` is the energy function (BQM) and :math:`Z(T)` the partition
+    function (a normalization constant).
+    Given a sample set (:math:`S`), a temperature estimate establishes the
     temperature that is most likely to have produced the sample set.
     An effective temperature can be derived from a sample set by considering the
-    rate of excitations only. A maximum-pseudo-likelihood (MPL) estimator 
-    considers local excitations only, which are sufficient to establish a 
+    rate of excitations only. A maximum-pseudo-likelihood (MPL) estimator
+    considers local excitations only, which are sufficient to establish a
     temperature efficiently (in compute time and number of samples). If the BQM
-    consists of uncoupled variables then the estimator is equivalent to a 
+    consists of uncoupled variables then the estimator is equivalent to a
     maximum likelihood estimator.
 
-    The effective MPL temperature is defined by the solution T to 
+    The effective MPL temperature is defined by the solution T to
 
     .. math::
-       0 = \sum_i \sum_{s \in S} f_i(s) \exp(f_i(s)/T), 
+       0 = \sum_i \sum_{s \in S} f_i(s) \exp(f_i(s)/T),
 
-    where f is the energy lost in flipping spin i against its current 
+    where f is the energy lost in flipping spin i against its current
     assignment (the effective field).
 
-    The problem is a convex root solving problem, and is solved with SciPy 
+    The problem is a convex root solving problem, and is solved with SciPy
     optimize.
+    In the case that all samples define local minima or maxima (all energies
+    en1 are same sign) SciPy is bypassed and the maximizing value (T=0)
+    is returned. This limit might be realized in low temperature samplesets.
 
     If the distribution is not Boltzmann with respect to the BQM provided, as
     may be the case for heuristic samplers (such as annealers), the temperature
-    estimate can be interpreted as characterizing only a rate of local 
-    excitations. In the case of sample sets obtained from D-Wave annealing 
-    quantum computers the temperature can be identified with a physical 
+    estimate can be interpreted as characterizing only a rate of local
+    excitations. In the case of sample sets obtained from D-Wave annealing
+    quantum computers the temperature can be identified with a physical
     temperature via a late-anneal freeze-out phenomena.
 
     Args:
         bqm (:obj:`dimod.BinaryQuadraticModel`, optional):
             Binary quadratic model describing sample distribution.
-            If ``bqm`` and ``site_energy`` are both None, then by default 
-            100 samples are drawn using :class:`~dwave.system.samplers.DWaveSampler`, 
+            If ``bqm`` and ``site_energy`` are both None, then by default
+            100 samples are drawn using :class:`~dwave.system.samplers.DWaveSampler`,
             with ``bqm`` defaulted as described.
-        sampleset (:class:`~dimod.SampleSet`, optional):
+        sampleset (samples_like or :class:`~dimod.SampleSet`, optional):
             A set of samples, assumed to be fairly sampled from
             a Boltzmann distribution characterized by ``bqm``.
-        site_energy (samples_like, optional):
-            A Tuple of effective fields and site labels.
+        en1 (nd.array, optional):
+            Effective fields as an np.ndarray (site labels not required).
             Derived from the ``bqm`` and ``sampleset`` if not provided.
+            First dimension indexes samples, second dimension indexes sites.
+            Ordering doesn't matter, but should be consistent with sample_weights.
         num_bootstrap_samples (int, optional, default=0):
-            Number of bootstrap estimators to calculate.
+            Number of bootstrap estimators to calculate. For now, the sampleset
+            must have uniform sample_weights to deploy this option. An aggregated
+            or weighted sampleset must be disaggregated with repetitions.
         seed (int, optional)
             Seeds the bootstrap method (if provided) allowing reproducibility
             of the estimators.
         T_guess (float, optional):
-            User approximation to the effective temperature, must be 
-            a positive scalar value.
+            User approximation to the effective temperature, must be
+            a positive (non-zero) scalar value.
             Seeding the root-search method can enable faster convergence.
             By default, T_guess is ignored if it falls outside the range
             of ``T_bracket``.
-        optimize_method (str,optional,default='bisect'):
-            SciPy method used for optimization. Options are 'bisect' and
-            None (the default SciPy optimize method). 
+        optimize_method (str, optional, default=None):
+            Optimize method used by SciPy ``root_scalar`` method. The default
+            method works well under default operation, 'bisect' can be
+            numerically more stable for the scalar case (temperature estimation
+            only.
         T_bracket (list or Tuple of 2 floats, optional, default=(0.001,1000)):
-            If excitations are absent, temperature is defined as zero, otherwise 
-            this defines the range of Temperatures over which to attempt a fit when
-            using the 'bisect' ``optimize_method`` (the default).
-
+            Relevant only if optimize_method='bisect'.
+            If excitations are absent, temperature is defined as zero, otherwise
+            this defines the range of Temperatures over which to attempt a fit when.
+        sample_weights (np.ndarray, optional):
+            A set of weights for the samples. If sampleset is of
+            type :obj:`~dimod.SampleSet` set this is default to
+            sampleset.record.num_occurrences, otherwise uniform weighting is
+            the default.
     Returns:
-        Tuple of float and NumPy array:
-            (T_estimate,T_bootstrap_estimates)
+        Tuple: The optimal parameters and a list of bootstrapped estimators (T_estimate, T_bootstrap_estimates):
 
-            *T_estimate*: a temperature estimate
-            *T_bootstrap_estimates*: a numpy array of bootstrap estimators
+        * *T_estimate*: a temperature estimate
+        * *T_bootstrap_estimates*: a list of bootstrap estimates
 
     Examples:
-       Draw samples from a D-Wave Quantum Computer for a large spin-glass 
+       Draw samples from a D-Wave Quantum Computer for a large spin-glass
        problem (random couplers J, zero external field h).
-       Establish a temperature estimate by maximum pseudo-likelihood. 
+       Establish a temperature estimate by maximum pseudo-likelihood.
 
        Note that due to the complicated freeze-out properties of hard models,
-       such as large scale spin-glasses, deviation from a classical Boltzmann 
+       such as large scale spin-glasses, deviation from a classical Boltzmann
        distribution is anticipated.
        Nevertheless, the T estimate can always be interpreted as an estimator
-       of local excitations rates. For example T will be 0 if only 
-       local minima are returned (even if some of the local minima are not 
+       of local excitations rates. For example T will be 0 if only
+       local minima are returned (even if some of the local minima are not
        ground states).
 
        >>> import dimod
        >>> from dwave.system.temperatures import maximum_pseudolikelihood_temperature
        >>> from dwave.system import DWaveSampler
        >>> from random import random
-       >>> sampler = DWaveSampler() 
+       >>> sampler = DWaveSampler()
        >>> bqm = dimod.BinaryQuadraticModel.from_ising({}, {e : 1-2*random() for e in sampler.edgelist})
        >>> sampleset = sampler.sample(bqm, num_reads=100, auto_scale=False)
-       >>> T,T_bootstrap =  maximum_pseudolikelihood_temperature(bqm,sampleset) 
+       >>> T,T_bootstrap =  maximum_pseudolikelihood_temperature(bqm, sampleset)
        >>> print('Effective temperature ',T)    # doctest: +SKIP
        Effective temperature  0.24066488780293813
 
@@ -249,164 +394,552 @@ def maximum_pseudolikelihood_temperature(bqm=None,
 
         https://www.jstor.org/stable/25464568
 
-    '''
-
-    T_estimate = 0
-    T_bootstrap_estimates = np.zeros(num_bootstrap_samples)
-
-    #Check for largest local excitation in every sample, and over all samples
-    if site_energy is None:
-        if bqm is None or sampleset is None:
-            raise ValueError('site_energy can only be derived if both'
-                             'bqm and sampleset are provided as arguments')
-        site_energy = effective_field(bqm,
-                                      sampleset,
-                                      current_state_energy = True)
-    max_excitation = np.max(site_energy[0],axis=1)
-    max_excitation_all =  np.max(max_excitation)
-    if max_excitation_all <= 0:
-        #There are no local excitations present in the sample set, therefore
-        #the temperature is estimated as 0. 
-        pass
+    """
+    if T_guess:
+        x0 = [-1 / T_guess]
     else:
-        
+        x0 = None
+
+    if optimize_method == "bisect":
+        bisect_bracket = [-1 / T_bracket[i] for i in range(2)]
+
+        if not 0 <= T_bracket[0] < T_bracket[1]:
+            raise ValueError("Bad T_bracket, must be positive ordered scalars.")
+
+        if x0 is not None and (x0[0] < bisect_bracket[0] or x0[0] > bisect_bracket[1]):
+            raise ValueError(f"x0 {x0[0]} and T_bracket {T_bracket} are inconsistent")
+
+        kwargs_opt = {"bracket": bisect_bracket}
+    else:
+        kwargs_opt = {}
+
+    x, x_bs = maximum_pseudolikelihood(
+        en1=en1,
+        bqms=[bqm],
+        sampleset=sampleset,
+        num_bootstrap_samples=num_bootstrap_samples,
+        seed=seed,
+        x0=x0,
+        optimize_method=optimize_method,
+        kwargs_opt=kwargs_opt,
+        sample_weights=sample_weights,
+    )
+    # exp(-1/Teff Hp) = exp(estimator[0] Hp)
+    return -1 / x, -1 / x_bs
+
+
+def _get_en1_for_pseudo_likelihood(en1, bqms, sampleset, sample_weights):
+    """Get effective fields for maximum pseudolikelihood."""
+    if en1 is None:
+        if bqms is None or sampleset is None:
+            raise ValueError(
+                "en1 can only be derived if both "
+                "bqms and sampleset are provided as arguments"
+            )
+
+        en1 = np.array(
+            [
+                effective_field(bqm, sampleset, current_state_energy=True)[0]
+                for bqm in bqms
+            ]
+        )
+
+    if en1.ndim == 2 or en1.shape[0] == 1:
+        # Scalar method 'inverse temperature only' for one Hamiltonian
+        en1 = en1.reshape(en1.shape[-2:])  # f_{i}(s) - column i, row s.
+
+    if sample_weights is None:
+        if isinstance(sampleset, dimod.sampleset.SampleSet):
+            sample_weights = sampleset.record.num_occurrences / np.sum(
+                sampleset.record.num_occurrences
+            )
+        else:
+            sample_weights = np.ones(en1.shape[-2])
+
+    if len(sample_weights) != en1.shape[-2]:
+        raise ValueError(
+            "The sample weights, if not defaulted, must match the sampleset shape, "
+            f"sample_weights.shape={sample_weights.shape}, en1.shape[-2]=={en1.shape[-2]}"
+        )
+
+    if any(sample_weights < 0) or np.max(sample_weights) == 0:
+        raise ValueError(
+            "sample weights must be non-negative with atleast one positive value"
+        )
+
+    if any(sample_weights == 0):
+        sample_weights = sample_weights[sample_weights > 0]
+        en1 = en1[:, sample_weights > 0, :]
+
+    return en1, sample_weights
+
+
+def _x_bootstrap_pseudo_likelihood(en1, max_excitation, num_bootstrap_samples):
+    """Get parameter estimates and bootstrap estimators for special case."""
+    # Only local minima (or maxima) observed
+    if max_excitation == 0:
+        warnings.warn(
+            "All local fields are zero, there is no gradient and the "
+            "parameter value is unconstrained. Zero is assigned."
+        )
+        x = 0
+    else:
+        x = np.sign(max_excitation) * float("Inf")
+
+    if en1.ndim == 3:
+        warnings.warn(
+            "An exponential model with ill-defined parameters fits the data; consider "
+            "taking more samples, modifying bqms or exploiting a perturbative approach."
+        )
+        x_bootstraps = x[np.newaxis, :] * np.ones((num_bootstrap_samples, 1))
+    else:
+        # Infinite results from all local minima, well defined limit
+        x_bootstraps = x * np.ones(num_bootstrap_samples)
+
+    return x, x_bootstraps
+
+
+def _create_d_mean_log_pseudo_likelihood(en1, degenerate_fields, sample_weights):
+    """Construct d_mean_log_pseudo_likelihood function for root finding routine."""
+    if en1.ndim == 2 or en1.shape[0] == 1:
+        # Derivative of mean (w.r.t samples) log pseudo liklihood amounts
+        # to local energy matching criteria
+        # O = sum_i \sum_s w(s) f_i(s) P(s, i) #s = sample, i = variable index
+        # f_i(s) is energy lost in flipping spin s_i against current assignment.
+        # P(s, i) = 1/(1 + exp[x f_i(s)]), probability to flip against current state.
+        if degenerate_fields:
+            # Histogram effective fields for faster processing
+            sw_replicated = np.tile(
+                sample_weights[:, np.newaxis], reps=(1, en1.shape[-1])
+            )
+            en1_u = np.unique(en1)
+            sw_u = np.histogram(
+                en1, bins=np.append(en1_u, float("Inf")), weights=sw_replicated
+            )[0]
+
+            def d_mean_log_pseudo_likelihood(x):
+                with warnings.catch_warnings():  # Overflow errors are safe
+                    warnings.simplefilter(action="ignore", category=RuntimeWarning)
+                    expFactor = np.exp(en1_u * x)
+                return np.sum(sw_u * en1_u / (1 + expFactor))
+
+        else:
+            # Faster if degeneracy of local fields is small
+            def d_mean_log_pseudo_likelihood(x):
+                with warnings.catch_warnings():  # Overflow errors are safe
+                    warnings.simplefilter(action="ignore", category=RuntimeWarning)
+                    expFactor = np.exp(en1 * x)
+                return np.sum(sample_weights * np.sum(en1 / (1 + expFactor), axis=1))
+
+    else:
+        if en1.ndim != 3:
+            raise ValueError("en1 must be 2d for a single bqm, or 3d for multiple bqms")
+
+        if degenerate_fields is not False:
+            raise ValueError(
+                "Exploiting degenerate multi-dimensional fields is not supported."
+            )
+
         def d_mean_log_pseudo_likelihood(x):
-            #Derivative of mean (w.r.t samples) log pseudo liklihood amounts
-            #to local energy matching criteria
-            #O = sum_i \sum_s f_i(s) P(s,i) #s = sample, i = variable index
-            #f_i(s) is energy lost in flipping spin s_i against current assignment.
-            #P(s,i) = 1/(1 + exp[x f_i(s)]), probability to flip against current state.
-            #x = -1/T
-            with warnings.catch_warnings():
-                #Overflow errors are safe:
-                warnings.simplefilter(action='ignore', category=RuntimeWarning)
-                expFactor = np.exp(site_energy[0]*x)
-            return np.sum(site_energy[0]/(1 + expFactor))
-        
-        #Ensures good gradient method, except pathological cases
-        if T_guess is None:
-            x0 = -1/max_excitation_all
+            with warnings.catch_warnings():  # Overflow errors are safe
+                warnings.simplefilter(action="ignore", category=RuntimeWarning)
+                expFactor = np.exp(np.sum(en1 * x[:, np.newaxis, np.newaxis], axis=0))
+            return np.sum(
+                sample_weights[np.newaxis, :]
+                * np.sum(en1 / (1 + expFactor[np.newaxis, :, :]), axis=-1),
+                axis=-1,
+            )
+
+    return d_mean_log_pseudo_likelihood
+
+
+def _create_dd_mean_log_pseudo_likelihood(
+    en1, degenerate_fields, sample_weights, use_jacobian
+):
+    """Construct dd_mean_log_pseudo_likelihood function for root finding routine."""
+    if use_jacobian is False:
+        dd_mean_log_pseudo_likelihood = None
+
+    elif en1.ndim == 2:
+        if degenerate_fields:
+            # Histogram effective fields for faster processing
+            sw_replicated = np.tile(
+                sample_weights[:, np.newaxis], reps=(1, en1.shape[-1])
+            )
+            en1_u = np.unique(en1)
+            sw_u = np.histogram(
+                en1, bins=np.append(en1_u, float("Inf")), weights=sw_replicated
+            )[0]
+
+            def dd_mean_log_pseudo_likelihood(x):  # second (scalar) derivative
+                with warnings.catch_warnings():  # expFactor=Inf causes an irrelevant warning
+                    warnings.simplefilter(action="ignore", category=RuntimeWarning)
+                    expFactor = np.exp(en1_u * x)
+                    norm = expFactor + 2 + 1 / expFactor
+                return -np.sum(sw_u * en1_u * en1_u / norm)
+
         else:
-            if T_guess < T_bracket[0]:
-                x0 = -1/T_bracket[0]
-            elif T_guess > T_bracket[1]:
-                x0 = -1/T_bracket[1]
-            else:
-                x0 = -1/T_guess
-        root_results = None
-        if optimize_method == 'bisect':
-            #Check T_bracket
-            if not 0 <= T_bracket[0] < T_bracket[1]:
-                raise ValueError('Bad T_bracket, must be positive ordered scalars.')
-            #Convert to bisection bracket for -1/T
-            bisect_bracket = [-1/T_bracket[i] for i in range(2)]
-            
-            if d_mean_log_pseudo_likelihood(bisect_bracket[0]) <0:
+
+            def dd_mean_log_pseudo_likelihood(x):
+                with warnings.catch_warnings():  # expFactor=Inf causes an irrelevant warning
+                    warnings.simplefilter(action="ignore", category=RuntimeWarning)
+                    expFactor = np.exp(en1 * x)
+                    norm = expFactor + 2 + 1 / expFactor
+                return -np.sum(sample_weights * np.sum(en1 * en1 / norm, axis=1))
+
+    else:
+
+        def dd_mean_log_pseudo_likelihood(x):  # jacobian
+            # [Only] If few Hamiltonians efficient, this is the expected use case
+            with warnings.catch_warnings():  # expFactor=Inf causes an irrelevant warning
+                warnings.simplefilter(action="ignore", category=RuntimeWarning)
+                expFactor = np.exp(np.sum(en1 * x[:, np.newaxis, np.newaxis], axis=0))
+                norm = expFactor + 2 + 1 / expFactor
+            return -np.sum(
+                [
+                    [
+                        sample_weights
+                        * np.sum(en1[i, :, :] * en1[j, :, :] / norm, axis=1)
+                        for i in range(en1.shape[0])
+                    ]
+                    for j in range(en1.shape[0])
+                ],
+                axis=-1,
+            )
+
+    return dd_mean_log_pseudo_likelihood
+
+
+def maximum_pseudolikelihood(
+    en1: Optional[np.ndarray] = None,
+    bqms: Union[None, List[dimod.BinaryQuadraticModel]] = None,
+    sampleset: Union[None, dimod.SampleSet, Tuple[np.ndarray, List]] = None,
+    num_bootstrap_samples: int = 0,
+    seed: Optional[int] = None,
+    x0: Union[None, List, np.ndarray] = None,
+    optimize_method: Optional[str] = None,
+    kwargs_opt: Optional[dict] = None,
+    sample_weights: Optional[np.ndarray] = None,
+    return_optimize_object: bool = False,
+    degenerate_fields: Optional[bool] = None,
+    use_jacobian: bool = True,
+) -> Tuple:
+    """Maximimum pseudolikelihood estimator for exponential models.
+
+    Uses the SciPy optimize method to solve the maximum pseudolikelihood problem
+    of weight estimation for an exponential model with exponent defined by a
+    weighted sum of bqms.
+    Assuming the Probability of states s is given by an exponential model
+    :math:`P(s)=1/Z exp(- x_i H_i(s))` for a given list of functions (bqms)
+    estimate x. The exponential model is also called a Boltzmann distribution.
+    Code is designed assuming the bqms are dense (a function of all or most of
+    the variables), although technically operation although sparse bqms
+    like :math:`H(s)=h_i s_i` or :math:`H(s)=J_{ij}s_i s_j` are technically allowed.
+
+    In the case that all samples define local minima or maxima (all energies
+    en1 are same sign) the limit of +/-Inf defines a maximum. In this special
+    case scipy is not called and the limit result is returned. Note that,
+    excitations can be rare in low-temperature samplesets such as those
+    derived from annealing, so that this pathological case can be realized in
+    applications of interest.
+
+
+    Note common reasons for parameter inference failure include:
+    - Too few samples, insufficient to resolve parameters
+    - The samplest is singular or insensitive with respect to some parameter
+      (e.g. local minima or plateus only).
+    - bqms are too closely related, or weakly dependent on the sampleset
+      variability (e.g. nearly collinear).
+
+    Args:
+        en1: Effective fields as an np.ndarray (site labels not required).
+            Derived from the ``bqms`` and ``sampleset`` if not provided.
+            First dimension indexes samples, second dimension indexes sites.
+            Ordering doesn't matter, but should be consistent with sample_weights.
+        bqms: A list of Binary quadratic models [H_i] describing the sample
+            distribution as :math:`P(s) ~ exp(sum_i x_i H_i(s))` for unknown
+            model parameters x.
+        sampleset: A set of samples, as a dimod Sampleset or samples-like object.
+        num_bootstrap_samples: Number of bootstrap estimators to calculate.
+            Bootstrapped estimates can be used to reliably estimate variance and
+            bias if samples are uncorrelated. For
+            now, the sampleset must have uniform sample_weights to deploy this
+            option -- an aggregated or weighted sampleset must be disaggregated
+            (raw format) with repetitions.
+        seed: Seeds the bootstrap method (if provided) allowing reproducibility
+            of the estimators.
+        x0: Initial guess for the fitting parameters. Should have the same
+            length as bqms when provided.
+        optimize_method (str, optional, default=None):
+            Optimize method used by SciPy ``root_scalar`` method. The default
+            method works well under default operation, 'bisect' can be
+            numerically more stable for the scalar case (inverse temperature
+            estimation only).
+        kwargs_opt: Arguments used by the SciPy optimization methods. If using the
+            'bisect' optimization method and for a single bqm, bounds for the fitting
+            parameter can be set using a ``tuple[float, float]`` for ``bracket``.
+        sample_weights: A set of weights for the samples. If sampleset is of
+            type :obj:`~dimod.SampleSet` set this is default to
+            sampleset.record.num_occurrences, otherwise uniform weighting is
+            the default.
+        return_optimize_object: When True, and if ``sciPy`` optimization is
+            invoked, the associated OptimizeResult is returned. Otherwise only
+            the optimal parameter values are returned as floats.
+        degenerate_fields: If effective fields are degenerate owing to sparse
+            connectivity, low precision and/or large number of samples, or low
+            entropy then histogramming (aggregating) of fields is used to
+            accelerate the search stage. A value True is supported (and default)
+            for single parameter esimation ``len(bqms)=1``; for multi-parameter
+            estimation only value False is supported.
+        use_jacobian: By default (True) the second derivative of the
+            pseudolikelihood is calculated and used by ScipY root finding
+            methods. The associated complexity of this non-essential
+            calculation is quadratic in len(bqms); use of the second derivative
+            is disabled by setting the value to False.
+    Returns:
+        Tuple: The optimal parameters and a list of bootstrapped estimates (x_estimate, x_bootstrap_estimates):
+
+        * *x_estimate*: parameter estimates
+        * *x_bootstrap_estimates*: a numpy array of bootstrap estimators
+
+    Examples:
+        This example builds upon the example for maximum_pseudolikelihood_temperature
+
+        Draw samples from a D-Wave Quantum Computer for a large spin-glass
+        problem (random couplers J, zero external field h).
+        Establish a temperature and estimate to the background susceptibility chi
+        which is expected to be a small negative value.
+        Since background susceptiblity is a perturbation, and the problem Hamiltonian
+        and correction Hamiltonian are correlated, a large number of samples can be
+        required for confidence. Other perturbative corrections such as flux noise
+        can also interfere with estimation of small parameters like chi.
+
+        >>> import dimod
+        >>> from dwave.system.temperatures import maximum_pseudolikelihood
+        >>> from dwave.system.temperatures import background_susceptibility_bqm
+        >>> from dwave.system import DWaveSampler
+        >>> from random import random
+        >>> sampler = DWaveSampler()
+        >>> bqm1 = dimod.BinaryQuadraticModel.from_ising({}, {e : 1-2*random() for e in sampler.edgelist})
+        >>> bqm2 = background_susceptibility_bqm(bqm1)
+        >>> sampleset = sampler.sample(bqm1, num_reads=1000, auto_scale=False)
+        >>> params, _ =  maximum_pseudolikelihood(bqms=[bqm1, bqm2], sampleset=sampleset)
+        >>> print('Effective temperature ', -1/params[0], 'Background susceptibliity', params[1]/params[0])    # doctest: +SKIP
+        Effective temperature  0.22298662677716122 Background susceptibliity -0.009343961890466117
+
+    See also:
+        The function :class:`~dwave.system.temperatures.maximum_pseudolikelihood_temperature`
+        https://doi.org/10.3389/fict.2016.00023
+
+        https://www.jstor.org/stable/25464568
+
+    """
+    if kwargs_opt is None:
+        kwargs_opt = {}
+
+    root_results = None
+
+    if degenerate_fields is None:
+        if en1 is None:
+            degenerate_fields = len(bqms) == 1
+        else:
+            degenerate_fields = en1.ndim == 2
+
+    en1, sample_weights = _get_en1_for_pseudo_likelihood(
+        en1, bqms, sampleset, sample_weights
+    )
+
+    max_excitation = np.max(en1, axis=(-1, -2))
+    min_excitation = np.min(en1, axis=(-1, -2))
+    prod_minmax = max_excitation * min_excitation
+
+    if (en1.ndim == 2 and prod_minmax >= 0) or en1.ndim == 3 and all(prod_minmax >= 0):
+        return _x_bootstrap_pseudo_likelihood(
+            en1, max_excitation, num_bootstrap_samples
+        )
+
+    else:
+        # To do: A simple example that shows an estimate of background
+        # susceptiblity separated from (inverse temperature) may be in order:
+        # Consider H(s) = - h s_1 - h s_3 + 2h s_2 - s_1 s_2 - s_2 s_3
+        # e.g. suppose the sampleset contains only + + and - - .
+        # Now we have the perturbative correction ..
+        # We can infer chi based on ground state manifold (only)! Assuming
+        # perturbative.
+        # There are no local excitations present in the sample set, therefore
+        # the temperature is estimated as 0.'
+        if en1.ndim == 2 or en1.shape[0] == 1:
+            # Scalar method 'inverse temperature' for one Hamiltonian
+            en1 = en1.reshape(en1.shape[-2:])  # f_{i}(s) - column i, row s.
+            if x0 is None:
+                x0 = -1 / max_excitation
+            elif hasattr(x0, "__iter__"):
+                x0 = x0[0]  # Scalar
+        else:
+            # Multi-parameter generalization
+            if x0 is None:
+                # Assume all bqms except the first are perturbative in absence
+                # of additional context.
+                x0 = np.zeros(en1.shape[0])
+                x0[0] = -1 / max_excitation[0]  # Smallest gap
+
+        d_mean_log_pseudo_likelihood = _create_d_mean_log_pseudo_likelihood(
+            en1, degenerate_fields, sample_weights
+        )
+
+        if optimize_method == "bisect" and en1.ndim == 2:
+            bisect_bracket = kwargs_opt.pop("bracket", (1e-3, 1000))
+            # bisect can be relatively robust, since we can have a problem of
+            # vanishing gradients given conservative bounds.
+            if d_mean_log_pseudo_likelihood(bisect_bracket[0]) < 0:
                 warnings.warn(
-                    'Temperature is less than T_bracket[0], or perhaps negative: ' 
-                    'rescaling the Hamiltonian, modification of T_bracket[0], or '
-                    'a change of the optimize_method (to None) '
-                    'can resolve this issue assuming a numerical precision issue '
-                    'induced by very large or small effective fields is not the cause. '
-                    'Automated precision requirements mean that this routine works '
-                    'best when excitations (effective fields) are O(1).')
-                T_estimate = T_bracket[0]
+                    "value should be positive at bisect_bracket[0]."
+                    "value returned matches the lower bound"
+                    "This might be resolved by a more sensible "
+                    "scaling of the bqm, or a change of the optimize_method."
+                )
+                x = bisect_bracket[0]
             elif d_mean_log_pseudo_likelihood(bisect_bracket[1]) > 0:
                 warnings.warn(
-                    'Temperature is greater than T_bracket[1], or perhaps negative:'
-                    'rescaling the Hamiltonian, modification of T_bracket[1] or '
-                    'a change of the optimize_method (to None) '
-                    'can resolve this issue assuming a numerical precision issue '
-                    'induced by very large or small effective fields is not the cause. '
-                    'Automated precision requirements mean that this routine works '
-                    'best when excitations (effective fields) are O(1).')
-                T_estimate = T_bracket[1]
+                    "value should be positive at bisect_bracket[0]. "
+                    "value returned matches the upper bound. "
+                    "This might be resolved by a more sensible "
+                    "scaling of the bqm, or a change of the optimize_method."
+                )
+                x = bisect_bracket[1]
             else:
-                #Bounds are good:
-                if x0 < bisect_bracket[0] or x0 > bisect_bracket[1]:
-                    x0 = (bisect_bracket[0] + bisect_bracket[1])/2 
-                root_results = optimize.root_scalar(f=d_mean_log_pseudo_likelihood, x0 = x0,
-                                                    method=optimize_method, bracket=bisect_bracket)
-                T_estimate = -1/(root_results.root)
+                root_results = optimize.root_scalar(
+                    f=d_mean_log_pseudo_likelihood,
+                    x0=x0,
+                    method=optimize_method,
+                    bracket=bisect_bracket,
+                    **kwargs_opt,
+                )
+                x = root_results.root
         else:
-            #For very large or small site energies this may be numerically
-            #unstable, therefore bisection search is the default.
-            def dd_mean_log_pseudo_likelihood(x):
-                
-                import warnings
-                with warnings.catch_warnings():
-                    warnings.simplefilter(action='ignore', category=RuntimeWarning)
-                    #overflow errors are harmless, +Inf is a safe value.
-                    expFactor = np.exp(site_energy[0]*x)
-                    #divide by zero (1/expFactor) and divide by +Inf errors are harmless
-                    return np.sum(-site_energy[0]*site_energy[0]/(expFactor + 2 + 1/expFactor))
-
-            root_results = optimize.root_scalar(f=d_mean_log_pseudo_likelihood, x0 = x0,
-                                                fprime=dd_mean_log_pseudo_likelihood)
-            T_estimate = -1/root_results.root
+            # Typically preferrable to bisect, but can be numerically unstable
+            # given large variance in effective fields or poor initial condition
+            # choices.
+            dd_mean_log_pseudo_likelihood = _create_dd_mean_log_pseudo_likelihood(
+                en1, degenerate_fields, sample_weights, use_jacobian
+            )
+
+            # Use of root finding routines are preferred to fsolve; best option
+            # is not obvious
+            if en1.ndim == 2:
+                root_results = optimize.root_scalar(
+                    f=d_mean_log_pseudo_likelihood,
+                    x0=x0,
+                    x1=x0 / 2,  # Naive choice should be fine if secant defaulted.
+                    fprime=dd_mean_log_pseudo_likelihood,
+                    **kwargs_opt,
+                )
+                x = root_results.root
+            else:
+                root_results = optimize.root(
+                    fun=d_mean_log_pseudo_likelihood,
+                    x0=x0,
+                    jac=dd_mean_log_pseudo_likelihood,
+                    **kwargs_opt,
+                )
+                x = root_results.x
+        x0 = x
         if num_bootstrap_samples > 0:
-            #By bootstrapping with respect to samples we 
-            if root_results is not None:
-                x0 = root_results.root
+            if len(np.unique(sample_weights)) != 1:
+                raise ValueError(
+                    "Bootstraps require uniform sample_weights (num_occurrences)"
+                )
             prng = np.random.RandomState(seed)
-            num_samples = site_energy[0].shape[0]
-            for bs in range(num_bootstrap_samples):
-                indices = np.random.choice(
-                    num_samples,
-                    num_bootstrap_samples,
-                    replace=True)
-                T_bootstrap_estimates[bs],_ = maximum_pseudolikelihood_temperature(
-                    site_energy = (site_energy[0][indices,:],site_energy[1]),
-                    num_bootstrap_samples = 0,
-                    T_guess = T_estimate)
-
-    return T_estimate, T_bootstrap_estimates
-
-def Ip_in_units_of_B(Ip: Union[None, float, np.ndarray]=None,
-                     B: Union[None, float, np.ndarray]=1.391,
-                     MAFM: Optional[float]=1.647,
-                     units_Ip: Optional[str]='uA',
-                     units_B: Literal['GHz', 'J'] = 'GHz',
-                     units_MAFM : Optional[str]='pH') -> Union[float, np.ndarray]:
+            num_samples = en1.shape[0]
+            x_bootstraps = []
+            for _ in range(num_bootstrap_samples):
+                indices = prng.choice(num_samples, num_bootstrap_samples, replace=True)
+                if en1.ndim == 2:
+                    x_bs, _ = maximum_pseudolikelihood(
+                        en1=en1[indices, :],
+                        num_bootstrap_samples=0,
+                        x0=x,
+                        optimize_method=optimize_method,
+                        kwargs_opt=kwargs_opt,
+                        return_optimize_object=return_optimize_object,
+                    )
+                else:
+                    x_bs, _ = maximum_pseudolikelihood(
+                        en1=en1[:, indices, :],
+                        num_bootstrap_samples=0,
+                        x0=x,
+                        optimize_method=optimize_method,
+                        kwargs_opt=kwargs_opt,
+                        return_optimize_object=return_optimize_object,
+                    )
+                x_bootstraps.append(x_bs)
+
+            if return_optimize_object is False:
+                x_bootstraps = np.array(x_bootstraps)
+
+        else:
+            if return_optimize_object is False:
+                if en1.ndim == 2:
+                    x_bootstraps = np.empty(0)
+                else:
+                    x_bootstraps = np.empty((0, len(x)))
+            else:
+                x_bootstraps = []
+
+    if return_optimize_object and root_results is not None:
+        # Return supplementary information on the optimization when available
+        return root_results, x_bootstraps
+
+    return x, x_bootstraps
+
+
+def Ip_in_units_of_B(
+    Ip: Union[None, float, np.ndarray] = None,
+    B: Union[None, float, np.ndarray] = 1.391,
+    MAFM: Optional[float] = 1.647,
+    units_Ip: Optional[str] = "uA",
+    units_B: Literal["GHz", "J"] = "GHz",
+    units_MAFM: Optional[str] = "pH",
+) -> Union[float, np.ndarray]:
     r"""Estimate qubit persistent current :math:`I_p(s)` in schedule units.
 
-    Under a simple, noiseless freeze-out model, you can substitute flux biases 
-    for programmed linear biases, ``h``, in the standard transverse-field Ising 
-    model as implemented on D-Wave quantum computers. Perturbations in ``h`` are 
-    not, however, equivalent to flux perturbations with respect to dynamics 
-    because of differences in the dependence on the anneal fraction, :math:`s`: 
-    :math:`I_p(s) \propto \sqrt(B(s))`. The physical origin of each term is different, 
+    Under a simple, noiseless freeze-out model, you can substitute flux biases
+    for programmed linear biases, ``h``, in the standard transverse-field Ising
+    model as implemented on D-Wave quantum computers. Perturbations in ``h`` are
+    not, however, equivalent to flux perturbations with respect to dynamics
+    because of differences in the dependence on the anneal fraction, :math:`s`:
+    :math:`I_p(s) \propto \sqrt(B(s))`. The physical origin of each term is different,
     and so precision and noise models also differ.
 
     Assume a Hamiltonian in the documented form,
     :math:numref:`qpu_equation_quantum_hamiltonian`, of the :ref:`qpu_annealing`
-    section with an additional flux-bias-dependent component 
+    section with an additional flux-bias-dependent component
     :math:`H(s) \Rightarrow H(s) - H_F(s) \sum_i \Phi_i \sigma^z_i`,
-    where :math:`\Phi_i` are flux biases (in units of :math:`\Phi_0`), 
-    :math:`\sigma^z_i` is the Pauli-z operator, and 
-    :math:`H_F(s) = Ip(s) \Phi_0`. Schedules for D-Wave quantum computers 
-    specify energy in units of Joule or GHz. 
+    where :math:`\Phi_i` are flux biases (in units of :math:`\Phi_0`),
+    :math:`\sigma^z_i` is the Pauli-z operator, and
+    :math:`H_F(s) = Ip(s) \Phi_0`. Schedules for D-Wave quantum computers
+    specify energy in units of Joule or GHz.
 
     Args:
         Ip:
-            Persistent current, :math:`I_p(s)`, in units of amps or 
-            microamps. When not provided, inferred from :math:`M_{AFM}` 
-            and and :math:`B(s)` based on the relation 
-            :math:`B(s) = 2 M_{AFM} I_p(s)^2`. 
+            Persistent current, :math:`I_p(s)`, in units of amps or
+            microamps. When not provided, inferred from :math:`M_{AFM}`
+            and and :math:`B(s)` based on the relation
+            :math:`B(s) = 2 M_{AFM} I_p(s)^2`.
 
         B:
-            Annealing schedule field, :math:`B(s)`, associated with the 
-            problem Hamiltonian. Schedules are provided for each quantum 
-            computer in the :ref:`qpu_solver_properties_specific` section. 
+            Annealing schedule field, :math:`B(s)`, associated with the
+            problem Hamiltonian. Schedules are provided for each quantum
+            computer in the :ref:`qpu_solver_properties_specific` section.
             This parameter is ignored when ``Ip`` is specified.
 
         MAFM:
-            Mutual inductance, :math:`M_{AFM}`, specified for each quantum 
-            computer in the :ref:`qpu_solver_properties_specific` section. 
+            Mutual inductance, :math:`M_{AFM}`, specified for each quantum
+            computer in the :ref:`qpu_solver_properties_specific` section.
             ``MAFM`` is ignored when ``Ip`` is specified.
 
         units_Ip:
-            Units in which the persistent current, ``Ip``, is specified. 
+            Units in which the persistent current, ``Ip``, is specified.
             Allowed values are ``'uA'`` (microamps) and ``'A'`` (amps)
 
         units_B:
@@ -414,80 +947,84 @@ def Ip_in_units_of_B(Ip: Union[None, float, np.ndarray]=None,
             are ``'GHz'`` (gigahertz) and ``'J'`` (Joules).
 
         units_MAFM:
-            Units in which the mutual inductance, ``MAFM``, is specified. Allowed 
+            Units in which the mutual inductance, ``MAFM``, is specified. Allowed
             values are ``'pH'`` (picohenry) and ``'H'`` (Henry).
 
     Returns:
         :math:`I_p(s)` with units matching the Hamiltonian :math:`B(s)`.
     """
-    h = 6.62607e-34  # Plank's constant for converting energy in Hertz to Joules 
+    h = 6.62607e-34  # Plank's constant for converting energy in Hertz to Joules
     Phi0 = 2.0678e-15  # superconducting magnetic flux quantum (h/2e); units: Weber=J/A
 
-    if units_B == 'GHz':
-        B_multiplier = 1e9*h  # D-Wave schedules use GHz by convention
-    elif units_B == 'J':
+    if units_B == "GHz":
+        B_multiplier = 1e9 * h  # D-Wave schedules use GHz by convention
+    elif units_B == "J":
         B_multiplier = 1
     else:
-        raise ValueError('Schedule B must be in units GHz or J, ' 
-                         f'but given {units_B}')
+        raise ValueError(
+            "Schedule B must be in units GHz or J, " f"but given {units_B}"
+        )
     if Ip is None:
-        B = B*B_multiplier # To Joules
-        if units_MAFM == 'pH':
-            MAFM = MAFM*1e-12  # conversion from picohenry to Henry
-        elif units_MAFM != 'H':
-            raise ValueError('MAFM must be in units pH or H, ' 
-                             f'but given {units_MAFM}')
-        Ip = np.sqrt(B/(2*MAFM))  # Units of A = C/s, O(1e-6) 
+        B = B * B_multiplier  # To Joules
+        if units_MAFM == "pH":
+            MAFM = MAFM * 1e-12  # conversion from picohenry to Henry
+        elif units_MAFM != "H":
+            raise ValueError(
+                "MAFM must be in units pH or H, " f"but given {units_MAFM}"
+            )
+        Ip = np.sqrt(B / (2 * MAFM))  # Units of A = C/s, O(1e-6)
     else:
-        if units_Ip == 'uA':
-            Ip = Ip*1e-6  # Conversion from microamps to amp
-        elif units_Ip != 'A':
-            raise ValueError('Ip must be in units uA or A, ' 
-                             f'but given {units_Ip}')
-
-    return Ip*Phi0/B_multiplier
-
-
-def h_to_fluxbias(h: Union[float, np.ndarray]=1,
-                  Ip: Optional[float]=None,
-                  B: float=1.391, MAFM: Optional[float]=1.647,
-                  units_Ip: Optional[str]='uA',
-                  units_B : str='GHz',
-                  units_MAFM : Optional[str]='pH') -> Union[float, np.ndarray]:
+        if units_Ip == "uA":
+            Ip = Ip * 1e-6  # Conversion from microamps to amp
+        elif units_Ip != "A":
+            raise ValueError("Ip must be in units uA or A, " f"but given {units_Ip}")
+
+    return Ip * Phi0 / B_multiplier
+
+
+def h_to_fluxbias(
+    h: Union[float, np.ndarray] = 1,
+    Ip: Optional[float] = None,
+    B: float = 1.391,
+    MAFM: Optional[float] = 1.647,
+    units_Ip: Optional[str] = "uA",
+    units_B: str = "GHz",
+    units_MAFM: Optional[str] = "pH",
+) -> Union[float, np.ndarray]:
     r"""Convert problem Hamiltonian bias ``h`` to equivalent flux bias.
 
-    Unitless bias ``h`` is converted to the equivalent flux bias in units 
+    Unitless bias ``h`` is converted to the equivalent flux bias in units
     :math:`\Phi_0`, the magnetic flux quantum.
 
     The dynamics of ``h`` and flux bias differ, as described in the
     :func:`Ip_in_units_of_B` function.
-    Equivalence at a specific point in the anneal is valid under a 
+    Equivalence at a specific point in the anneal is valid under a
     freeze-out (quasi-static) hypothesis.
 
-    Defaults are based on the published physical properties of 
-    `Leap <https://cloud.dwavesys.com/leap/>`_\ 's  
+    Defaults are based on the published physical properties of
+    `Leap <https://cloud.dwavesys.com/leap/>`_\ 's
     ``Advantage_system4.1`` solver at single-qubit freezeout (:math:`s=0.612`).
 
     Args:
         Ip:
-            Persistent current, :math:`I_p(s)`, in units of amps or 
-            microamps. When not provided, inferred from :math:`M_{AFM}` 
-            and and :math:`B(s)` based on the relation 
-            :math:`B(s) = 2 M_{AFM} I_p(s)^2`. 
+            Persistent current, :math:`I_p(s)`, in units of amps or
+            microamps. When not provided, inferred from :math:`M_{AFM}`
+            and and :math:`B(s)` based on the relation
+            :math:`B(s) = 2 M_{AFM} I_p(s)^2`.
 
         B:
-            Annealing schedule field, :math:`B(s)`, associated with the 
-            problem Hamiltonian. Schedules are provided for each quantum 
+            Annealing schedule field, :math:`B(s)`, associated with the
+            problem Hamiltonian. Schedules are provided for each quantum
             computer in the :ref:`qpu_solver_properties_specific` section.
             This parameter is ignored when ``Ip`` is specified.
 
         MAFM:
-            Mutual inductance, :math:`M_{AFM}`, specified for each quantum 
+            Mutual inductance, :math:`M_{AFM}`, specified for each quantum
             computer in the :ref:`qpu_solver_properties_specific` section.
             ``MAFM`` is ignored when ``Ip`` is specified.
 
         units_Ip:
-            Units in which the persistent current, ``Ip``, is specified. 
+            Units in which the persistent current, ``Ip``, is specified.
             Allowed values are ``'uA'`` (microamps) and ``'A'`` (amps)
 
         units_B:
@@ -495,60 +1032,66 @@ def h_to_fluxbias(h: Union[float, np.ndarray]=1,
             are ``'GHz'`` (gigahertz) and ``'J'`` (Joules).
 
         units_MAFM:
-            Units in which the mutual inductance, ``MAFM``, is specified. Allowed 
+            Units in which the mutual inductance, ``MAFM``, is specified. Allowed
             values are ``'pH'`` (picohenry) and ``'H'`` (Henry).
 
     Returns:
-        Flux-bias values producing equivalent longitudinal fields to the given 
+        Flux-bias values producing equivalent longitudinal fields to the given
         ``h`` values.
     """
-    Ip = Ip_in_units_of_B(Ip, B, MAFM,
-                          units_Ip, units_B, units_MAFM)  # Convert/Create Ip in units of B, scalar
-    # B(s)/2 h_i = Ip(s) phi_i 
-    return -B/2/Ip*h
-
-def fluxbias_to_h(fluxbias: Union[float, np.ndarray]=1,
-                  Ip: Optional[float]=None,
-                  B: float=1.391, MAFM: Optional[float]=1.647, 
-                  units_Ip: Optional[str]='uA',
-                  units_B : str='GHz', units_MAFM : Optional[str]='pH') -> Union[float, np.ndarray]:
+    Ip = Ip_in_units_of_B(
+        Ip, B, MAFM, units_Ip, units_B, units_MAFM
+    )  # Convert/Create Ip in units of B, scalar
+    # B(s)/2 h_i = Ip(s) phi_i
+    return -B / 2 / Ip * h
+
+
+def fluxbias_to_h(
+    fluxbias: Union[float, np.ndarray] = 1,
+    Ip: Optional[float] = None,
+    B: float = 1.391,
+    MAFM: Optional[float] = 1.647,
+    units_Ip: Optional[str] = "uA",
+    units_B: str = "GHz",
+    units_MAFM: Optional[str] = "pH",
+) -> Union[float, np.ndarray]:
     r"""Convert flux biases to equivalent problem Hamiltonian bias ``h``.
 
-    Converts flux biases in units of :math:`\Phi_0`, the magnetic flux quantum, 
-    to equivalent problem Hamiltonian biases ``h``, which are unitless. 
+    Converts flux biases in units of :math:`\Phi_0`, the magnetic flux quantum,
+    to equivalent problem Hamiltonian biases ``h``, which are unitless.
 
     The dynamics of ``h`` and flux bias differ, as described in the
     :func:`Ip_in_units_of_B` function.
-    Equivalence at a specific point in the anneal is valid under a 
+    Equivalence at a specific point in the anneal is valid under a
     freeze-out (quasi-static) hypothesis.
 
-    Defaults are based on the published physical properties of 
-    `Leap <https://cloud.dwavesys.com/leap/>`_\ 's  
+    Defaults are based on the published physical properties of
+    `Leap <https://cloud.dwavesys.com/leap/>`_\ 's
     ``Advantage_system4.1`` solver at single-qubit freezeout (:math:`s=0.612`).
 
     Args:
-        fluxbias: 
+        fluxbias:
              A flux bias in units of :math:`\Phi_0`.
 
         Ip:
-            Persistent current, :math:`I_p(s)`, in units of amps or 
-            microamps. When not provided, inferred from :math:`M_{AFM}` 
-            and and :math:`B(s)` based on the relation 
-            :math:`B(s) = 2 M_{AFM} I_p(s)^2`. 
-    
+            Persistent current, :math:`I_p(s)`, in units of amps or
+            microamps. When not provided, inferred from :math:`M_{AFM}`
+            and and :math:`B(s)` based on the relation
+            :math:`B(s) = 2 M_{AFM} I_p(s)^2`.
+
         B:
-            Annealing schedule field, :math:`B(s)`, associated with the 
-            problem Hamiltonian. Schedules are provided for each quantum 
+            Annealing schedule field, :math:`B(s)`, associated with the
+            problem Hamiltonian. Schedules are provided for each quantum
             computer in the :ref:`qpu_solver_properties_specific` section.
             This parameter is ignored when ``Ip`` is specified.
 
         MAFM:
-            Mutual inductance, :math:`M_{AFM}`, specified for each quantum 
+            Mutual inductance, :math:`M_{AFM}`, specified for each quantum
             computer in the :ref:`qpu_solver_properties_specific` section.
             ``MAFM`` is ignored when ``Ip`` is specified.
 
         units_Ip:
-            Units in which the persistent current, ``Ip``, is specified. 
+            Units in which the persistent current, ``Ip``, is specified.
             Allowed values are ``'uA'`` (microamps) and ``'A'`` (amps)
 
         units_B:
@@ -556,38 +1099,41 @@ def fluxbias_to_h(fluxbias: Union[float, np.ndarray]=1,
             are ``'GHz'`` (gigahertz) and ``'J'`` (Joules).
 
         units_MAFM:
-            Units in which the mutual inductance, ``MAFM``, is specified. Allowed 
+            Units in which the mutual inductance, ``MAFM``, is specified. Allowed
             values are ``'pH'`` (picohenry) and ``'H'`` (Henry).
 
     Returns:
         ``h`` values producing equivalent longitudinal fields to the flux biases.
     """
-    Ip = Ip_in_units_of_B(Ip, B, MAFM,
-                          units_Ip, units_B, units_MAFM)  # Convert/Create Ip in units of B, scalar
-    # B(s)/2 h_i = Ip(s) phi_i 
-    return -2*Ip/B*fluxbias
+    Ip = Ip_in_units_of_B(
+        Ip, B, MAFM, units_Ip, units_B, units_MAFM
+    )  # Convert/Create Ip in units of B, scalar
+    # B(s)/2 h_i = Ip(s) phi_i
+    return -2 * Ip / B * fluxbias
 
 
-def freezeout_effective_temperature(freezeout_B, temperature, units_B = 'GHz', units_T = 'mK') -> float:
-    r'''Provides an effective temperature as a function of freezeout information.
+def freezeout_effective_temperature(
+    freezeout_B: float, temperature: float, units_B: str = "GHz", units_T: str = "mK"
+) -> float:
+    r"""Provides an effective temperature as a function of freezeout information.
 
-    See the :ref:`qpu_annealing` section for a 
+    See the :ref:`qpu_annealing` section for a
     complete summary of D-Wave annealing quantum computer operation.
 
     A D-Wave annealing quantum computer is assumed to implement a Hamiltonian
-    :math:`H(s) = B(s)/2 H_P - A(s)/2 H_D`, where: :math:`H_P` is the unitless 
+    :math:`H(s) = B(s)/2 H_P - A(s)/2 H_D`, where: :math:`H_P` is the unitless
     diagonal problem Hamiltonian,
     :math:`H_D` is the unitless driver Hamiltonian, :math:`B(s)` is the problem energy scale; A(s)
     is the driver energy scale, amd :math:`s` is the normalized anneal
     time :math:`s = t/t_a` (in [0,1]).
     Diagonal elements of :math:`H_P`, indexed by the spin state :math:`x`, are equal to
-    the energy of a classical Ising spin system 
+    the energy of a classical Ising spin system
 
     .. math::
-        E_{Ising}(x) = \sum_i h_i x_i + \sum_{i>j} J_{i,j} x_i x_j
+        E_{Ising}(x) = \sum_i h_i x_i + \sum_{i>j} J_{i, j} x_i x_j
 
-    If annealing achieves a thermally equilibrated distribution 
-    over decohered states at large :math:`s` where :math:`A(s) \ll B(s)`, 
+    If annealing achieves a thermally equilibrated distribution
+    over decohered states at large :math:`s` where :math:`A(s) \ll B(s)`,
     and dynamics stop abruptly at :math:`s=s^*`, the distribution of returned
     samples is well described by a Boltzmann distribution:
 
@@ -595,22 +1141,22 @@ def freezeout_effective_temperature(freezeout_B, temperature, units_B = 'GHz', u
         P(x) = \exp(- B(s^*) R E_{Ising}(x) / 2 k_B T)
 
     where T is the physical temperature, and :math:`k_B` is the Boltzmann constant.
-    R is a Hamiltonain rescaling factor, if a QPU is operated with auto_scale=False, 
+    R is a Hamiltonain rescaling factor, if a QPU is operated with auto_scale=False,
     then R=1.
     The function calculates the unitless effective temperature as :math:`T_{eff} = 2 k_B T/B(s^*)`.
 
-    Device temperature :math:`T`, annealing schedules {:math:`A(s)`, :math:`B(s)`} and 
-    single-qubit freeze-out (:math:`s^*`, for simple uncoupled Hamltonians) are reported 
+    Device temperature :math:`T`, annealing schedules {:math:`A(s)`, :math:`B(s)`} and
+    single-qubit freeze-out (:math:`s^*`, for simple uncoupled Hamltonians) are reported
     device properties: see the :ref:`qpu_solver_properties_specific` section.
-    These values (typically specified in mK and GHz) allows the calculation of an effective 
-    temperature for simple Hamiltonians submitted to D-Wave quantum computers. Complicated 
-    problems exploiting embeddings, or with many coupled variables, may freeze out at 
-    different values of s or piecemeal). Large problems may have slow dynamics at small 
+    These values (typically specified in mK and GHz) allows the calculation of an effective
+    temperature for simple Hamiltonians submitted to D-Wave quantum computers. Complicated
+    problems exploiting embeddings, or with many coupled variables, may freeze out at
+    different values of s or piecemeal). Large problems may have slow dynamics at small
     values of s, so :math:`A(s)` cannot be ignored as a contributing factor to the distribution.
 
     Note that for QPU solvers this temperature estimate applies to problems
-    submitted with no additional scaling factors (sampling with ``auto_scale = False``). 
-    If ``auto_scale=True`` (default) additional scaling factors must be accounted for. 
+    submitted with no additional scaling factors (sampling with ``auto_scale = False``).
+    If ``auto_scale=True`` (default) additional scaling factors must be accounted for.
 
     Args:
         freezeout_B (float):
@@ -619,22 +1165,22 @@ def freezeout_effective_temperature(freezeout_B, temperature, units_B = 'GHz', u
         temperature (float):
             :math:`T`, the physical temperature of the quantum computer.
 
-        units_B (string, optional, 'GHz'):
+        units_B (string, 'GHz'):
             Units in which ``freezeout_B`` is specified. Allowed values:
             'GHz' (Giga-Hertz) and 'J' (Joules).
 
-        units_T (string, optional, 'mK'):
+        units_T (string, 'mK'):
             Units in which the ``temperature`` is specified. Allowed values:
             'mK' (milli-Kelvin) and 'K' (Kelvin).
 
     Returns:
-        float : The effective (unitless) temperature. 
+        float : The effective (unitless) temperature.
 
     Examples:
 
-       This example uses the 
+       This example uses the
        :ref:`published parameters <qpu_solver_properties_specific>`_
-       for the Advantage_system4.1 QPU solver as of November 22nd 2021: 
+       for the Advantage_system4.1 QPU solver as of November 22nd 2021:
        :math:`B(s=0.612) = 3.91` GHz , :math:`T = 15.4` mK.
 
        >>> from dwave.system.temperatures import freezeout_effective_temperature
@@ -644,72 +1190,81 @@ def freezeout_effective_temperature(freezeout_B, temperature, units_B = 'GHz', u
 
     See also:
 
-        The function :class:`~dwave.system.temperatures.fast_effective_temperature` 
+        The function :class:`~dwave.system.temperatures.fast_effective_temperature`
         estimates the temperature for single-qubit Hamiltonians, in approximate
-        agreement with estimates by this function at reported single-qubit 
+        agreement with estimates by this function at reported single-qubit
         freeze-out values :math:`s^*` and device physical parameters.
-    '''
-
-    #Convert units_B to Joules
-    if units_B == 'GHz':
-        h = 6.62607e-34 #J/Hz
-        freezeout_B = freezeout_B *h
-        freezeout_B *= 1e9
-    elif units_B == 'J':
+    """
+    # Convert units_B to Joules
+    if units_B == "GHz":
+        h = 6.62607e-34  # J/Hz
+        freezeout_B *= h * 1e9
+    elif units_B == "J":
         pass
     else:
-        raise ValueException("Units must be 'J' (Joules) "
-                             "or 'mK' (milli-Kelvin)")
+        raise ValueError("Units must be 'J' (Joules) or 'mK' (milli-Kelvin)")
 
-    if units_T == 'mK':
+    if units_T == "mK":
         temperature = temperature * 1e-3
-    elif units_T == 'K':
+    elif units_T == "K":
         pass
     else:
-        raise ValueException("Units must be 'K' (Kelvin) "
-                             "or 'mK' (milli-Kelvin)")
-    kB = 1.3806503e-23 # J/K
+        raise ValueError("Units must be 'K' (Kelvin) or 'mK' (milli-Kelvin)")
+    kB = 1.3806503e-23  # Joules/Kelvin
 
-    return 2*temperature*kB/freezeout_B
+    return 2 * temperature * kB / freezeout_B
 
-def fast_effective_temperature(sampler=None, num_reads=None, seed=None,
-                               h_range=(-1/6.1,1/6.1), sampler_params=None,
-                               optimize_method=None,
-                               num_bootstrap_samples=0) -> Tuple[np.float64,np.float64]:
-    r'''Provides an estimate to the effective temperature, :math:`T`, of a sampler.
 
-    This function submits a set of single-qubit problems to a sampler and 
+def fast_effective_temperature(
+    sampler: dimod.Sampler,
+    num_reads: Optional[int] = None,
+    seed: Union[None, int, np.random.RandomState] = None,
+    h_range: Tuple = (-1 / 6.1, 1 / 6.1),
+    nodelist: Optional[List] = None,
+    sampler_params: dict = None,
+    optimize_method: Optional[str] = "bisect",
+    num_bootstrap_samples: Optional[int] = 0,
+) -> Tuple[np.float64, np.float64]:
+    r"""Provides an estimate to the effective temperature, :math:`T`, of a sampler.
+
+    This function submits a set of single-qubit problems to a sampler and
     uses the rate of excitations to infer a maximum-likelihood estimate of temperature.
+    For greater control of the problem Hamiltonian or more general samplers,
+    `maximum_pseudolikelihood_temperature` should be preferred.
 
     Args:
         sampler (:class:`dimod.Sampler`, optional, default=\ :class:`~dwave.system.samplers.DWaveSampler`):
-            A dimod sampler. 
+            A dimod sampler.
 
         num_reads (int, optional):
-            Number of reads to use. Default is 100 if not specified in 
+            Number of reads to use. Default is 100 if not specified in
             ``sampler_params``.
 
         seed (int, optional):
             Seeds the problem generation process. Allowing reproducibility
             from pseudo-random samplers.
 
-        h_range (float, optional, default = [-1/6.1,1/6.1]):
+        h_range (float, default = [-1/6.1,1/6.1]):
             Determines the range of external fields probed for temperature
             inference. Default is based on a D-Wave Advantage processor, where
-            single-qubit freeze-out implies an effective temperature of 6.1 
+            single-qubit freeze-out implies an effective temperature of 6.1
             (see :class:`~dwave.system.temperatures.freezeout_effective_temperature`).
             The range should be chosen inversely proportional to the anticipated
-            temperature for statistical efficiency, and to accomodate precision 
+            temperature for statistical efficiency, and to accomodate precision
             and other nonidealities such as precision limitations.
 
+        nodelist (list, optional):
+            Variables included in the estimator. If not provided defaulted to
+            the nodelist of the structured sampler (QPU) specified.
+
         sampler_params (dict, optional):
-            Any additional non-defaulted sampler parameterization. If 
-            ``num_reads`` is a key, must be compatible with ``num_reads`` 
+            Any additional non-defaulted sampler parameterization. If
+            ``num_reads`` is a key, must be compatible with ``num_reads``
             argument.
 
         optimize_method (str, optional):
             Optimize method used by SciPy ``root_scalar`` method. The default
-            method works well under default operation, 'bisect' can be 
+            method works well under default operation, 'bisect' can be
             numerically more stable when operated without defaults.
 
         num_bootstrap_samples (int, optional, default=0):
@@ -720,9 +1275,13 @@ def fast_effective_temperature(sampler=None, num_reads=None, seed=None,
     Returns:
         Tuple[float, float]:
             The effective temperature describing single qubit problems in an
-            external field, and a standard error (+/- 1 sigma). 
+            external field, and a standard error (+/- 1 sigma).
             By default the confidence interval is set as 0.
 
+    Raises:
+        ValueError: If the sampler is not structured, and no nodelist is provided.
+        ValueError: If num_reads is inconsistent with sampler_params
+
     See also:
 
         https://doi.org/10.3389/fict.2016.00023
@@ -730,72 +1289,85 @@ def fast_effective_temperature(sampler=None, num_reads=None, seed=None,
         https://www.jstor.org/stable/25464568
 
     Examples:
-       Draw samples from a :class:`~dwave.system.samplers.DWaveSampler`, and establish the temperature
+       Draw samples from a :class:`~dwave.system.samplers.DWaveSampler`
+       and establish the temperature
 
        >>> from dwave.system.temperatures import fast_effective_temperature
        >>> from dwave.system import DWaveSampler
        >>> sampler = DWaveSampler()
        >>> T, _ = fast_effective_temperature(sampler)
        >>> print('Effective temperature at freeze-out is',T)    # doctest: +SKIP
-       0.21685104745347336
+       Effective temperature at freeze-out is  0.21685104745347336
 
     See also:
-        The function :class:`~dwave.system.temperatures.freezeout_effective_temperature` 
-        may be used in combination with published device values to estimate single-qubit 
+        The function :class:`~dwave.system.temperatures.freezeout_effective_temperature`
+        may be used in combination with published device values to estimate single-qubit
         freeze-out, in approximate agreement with empirical estimates of this function.
 
         https://doi.org/10.3389/fict.2016.00023
 
         https://www.jstor.org/stable/25464568
-    '''
+    """
 
     if sampler is None:
         from dwave.system import DWaveSampler
+
         sampler = DWaveSampler()
 
-    if 'h_range' in sampler.properties:
-        if h_range[0] < sampler.properties['h_range'][0]:
-            raise ValueError('h_range[0] exceeds programmable range')
+    if "h_range" in sampler.properties:
+        if h_range[0] < sampler.properties["h_range"][0]:
+            raise ValueError("h_range[0] exceeds programmable range")
 
-        if h_range[1] > sampler.properties['h_range'][1]:
-            raise ValueError('h_range[1] exceeds programmable range')
+        if h_range[1] > sampler.properties["h_range"][1]:
+            raise ValueError("h_range[1] exceeds programmable range")
+    if nodelist is None:
+        if hasattr(sampler, "nodelist"):
+            nodelist = sampler.nodelist
+        else:
+            raise ValueError(
+                "nodelist is not provided and cannot be inferred from the sampler."
+            )
+    num_vars = len(nodelist)
 
     prng = np.random.RandomState(seed)
-    h_values = h_range[0] + (h_range[1]-h_range[0])*prng.rand(len(sampler.nodelist))
-    bqm = dimod.BinaryQuadraticModel.from_ising({var: h_values[idx] for idx,var in enumerate(sampler.nodelist)}, {})
+    h_values = h_range[0] + (h_range[1] - h_range[0]) * prng.rand(num_vars)
+    bqm = dimod.BinaryQuadraticModel.from_ising(
+        {var: h_values[idx] for idx, var in enumerate(nodelist)}, {}
+    )
 
-    #Create local sampling_params copy - default necessary additional fields:
+    # Create local sampling_params copy - default necessary additional fields:
     if sampler_params is None:
         sampler_params0 = {}
     else:
         sampler_params0 = sampler_params.copy()
     if num_reads is None:
-        #Default is 100, makes efficient use of QPU access time:
-        if 'num_reads' not in sampler_params0:
-            sampler_params0['num_reads'] = 100 
-    elif ('num_reads' in sampler_params0
-        and sampler_params0['num_reads'] != num_reads):
-        raise ValueError("sampler_params['num_reads'] != num_reads, "
-                         "incompatible input arguments.")
+        # Default is 1000, makes efficient use of QPU access time:
+        if "num_reads" not in sampler_params0:
+            sampler_params0["num_reads"] = 1000
+    elif sampler_params0.get("num_reads") != num_reads:
+        raise ValueError(
+            "sampler_params['num_reads'] != num_reads, incompatible input arguments."
+        )
     else:
-        sampler_params0['num_reads'] = num_reads
-    if ('auto_scale' in sampler_params0
-        and sampler_params0['auto_scale'] is not False):
-        raise ValueError("sampler_params['auto_scale'] == False, "
-                         "is required by this method.")
+        sampler_params0["num_reads"] = num_reads
+    if "auto_scale" in sampler_params0 and sampler_params0["auto_scale"] is not False:
+        raise ValueError(
+            "sampler_params['auto_scale'] == False, is required by this method."
+        )
     else:
-        sampler_params0['auto_scale'] = False
+        sampler_params0["auto_scale"] = False
 
     if num_bootstrap_samples is None:
-        num_bootstrap_samples = sampler_params0['num_reads'] 
+        num_bootstrap_samples = sampler_params0["num_reads"]
 
     sampleset = sampler.sample(bqm, **sampler_params0)
 
-    T,Tboot = maximum_pseudolikelihood_temperature(
+    T, Tboot = maximum_pseudolikelihood_temperature(
         bqm,
         sampleset,
         optimize_method=optimize_method,
-        num_bootstrap_samples=num_bootstrap_samples)
+        num_bootstrap_samples=num_bootstrap_samples,
+    )
 
     if num_bootstrap_samples == 0:
         return T, np.float64(0.0)


Test patch:
diff --git a/tests/test_temperatures.py b/tests/test_temperatures.py
index 1e683ef..720e0cf 100644
--- a/tests/test_temperatures.py
+++ b/tests/test_temperatures.py
@@ -15,37 +15,42 @@
 import unittest
 import numpy as np
 import dimod
+import warnings
 from itertools import product
 
-from dwave.system.temperatures import (maximum_pseudolikelihood_temperature,
-                                       effective_field,
-                                       freezeout_effective_temperature,
-                                       fast_effective_temperature,
-                                       Ip_in_units_of_B,
-                                       h_to_fluxbias,
-                                       fluxbias_to_h)
+from dwave.system.temperatures import (
+    maximum_pseudolikelihood,
+    maximum_pseudolikelihood_temperature,
+    effective_field,
+    freezeout_effective_temperature,
+    fast_effective_temperature,
+    Ip_in_units_of_B,
+    h_to_fluxbias,
+    fluxbias_to_h,
+    background_susceptibility_ising,
+    background_susceptibility_bqm,
+)
 
 from dwave.system.testing import MockDWaveSampler
 
+
 class TestTemperatures(unittest.TestCase):
     def test_Ip_in_units_of_B(self):
-        uBs = ['J', 'GHz']
-        uIps = ['A', 'uA']
-        uMAFMs = ['H', 'pH']
+        uBs = ["J", "GHz"]
+        uIps = ["A", "uA"]
+        uMAFMs = ["H", "pH"]
         for uIp, uB, uMAFM in product(uIps, uBs, uMAFMs):
-            _ = Ip_in_units_of_B(units_Ip=uIp,
-                                 units_B=uB,
-                                 units_MAFM=uMAFM)
+            Ip_in_units_of_B(units_Ip=uIp, units_B=uB, units_MAFM=uMAFM)
 
     def test_fluxbias_h(self):
         phi = np.random.random()
         h = fluxbias_to_h(phi)
         phi2 = h_to_fluxbias(h)
-        self.assertLess(abs(phi-phi2), 1e-15)
+        self.assertLess(abs(phi - phi2), 1e-15)
         phi = np.random.random(10)
         h = fluxbias_to_h(phi)
         phi2 = h_to_fluxbias(h)
-        self.assertTrue(np.all(np.abs(phi-phi2) < 1e-15))
+        self.assertTrue(np.all(np.abs(phi - phi2) < 1e-15))
 
     def test_effective_field(self):
         # For a simple model of independent spins H = sum_i s_i
@@ -55,112 +60,299 @@ class TestTemperatures(unittest.TestCase):
         num_samples = 2
         var_labels = list(range(num_var))
         bqm = dimod.BinaryQuadraticModel.from_ising({var: 1 for var in var_labels}, {})
-        samples_like = (np.ones(shape=(num_samples,num_var)),var_labels)
-        E = effective_field(bqm,
-                            samples_like)
-        self.assertTrue(np.array_equal(np.ones(shape=(num_samples,num_var)), E[0]))
-        self.assertTrue(num_var==len(E[1]))
+        samples_like = (np.ones(shape=(num_samples, num_var)), var_labels)
+        E = effective_field(bqm, samples_like)
+        self.assertTrue(np.array_equal(np.ones(shape=(num_samples, num_var)), E[0]))
+        self.assertEqual(num_var, len(E[1]))
         # energy lost in flipping from sample value (1) to -1 is H(1) - H(-1) = +2.
-        E = effective_field(bqm,
-                            samples_like,
-                            current_state_energy=True)
-        self.assertTrue(np.array_equal(2*np.ones(shape=(num_samples,num_var)), E[0]))
+        E = effective_field(bqm, samples_like, current_state_energy=True)
+        self.assertTrue(np.array_equal(2 * np.ones(shape=(num_samples, num_var)), E[0]))
 
     def test_effective_field_vartype(self):
         # Check effective fields are identical whether using bqm or ising model
         num_var = 4
         var_labels = list(range(num_var))
-        bqm = dimod.BinaryQuadraticModel.from_ising({var: var for var in var_labels}, {(var1,var2) : var1%2 + var2%2 - 1 for var1 in var_labels for var2 in var_labels})
-        E_ising = effective_field(bqm,current_state_energy=True)
-        bqm.change_vartype('BINARY',inplace=True)
-        E_bqm = effective_field(bqm,current_state_energy=True)
-        self.assertTrue(bqm.vartype==dimod.BINARY) 
+        bqm = dimod.BinaryQuadraticModel.from_ising(
+            {var: var for var in var_labels},
+            {
+                (var1, var2): var1 % 2 + var2 % 2 - 1
+                for var1 in var_labels
+                for var2 in var_labels
+            },
+        )
+        E_ising = effective_field(bqm, current_state_energy=True)
+        bqm.change_vartype("BINARY", inplace=True)
+        E_bqm = effective_field(bqm, current_state_energy=True)
+        self.assertTrue(bqm.vartype == dimod.BINARY)
         self.assertTrue(np.array_equal(E_ising[0], E_bqm[0]))
 
+    def test_background_susceptibility(self):
+        # A Hamiltonian with + + + and - - - as ground states.
+        # Symmetry is broken
+        n = 3
+        dh = 1 / 4
+        h = np.array([dh, -2 * dh, dh])
+        J = np.array([[0, -1, 0], [-1, 0, -1], [0, -1, 0]])
+        dh, dJ, k = background_susceptibility_ising(h, J)
+        # Assert expected dh and dJ values.
+        # ([2+3], [1+3], [1+2])
+
+        Jd = {
+            (n1, n2): J[n1, n2] for n2 in range(n) for n1 in range(n2) if J[n1, n2] != 0
+        }
+        hd = {n: h[n] for n in range(n)}
+        bqm = dimod.BinaryQuadraticModel("SPIN").from_ising(hd, Jd)
+        dh, dJ, _ = background_susceptibility_ising(hd, Jd)
+        # Assert sparse and dense method match
+        dbqm = dimod.BinaryQuadraticModel("SPIN").from_ising(dh, dJ)
+
+        chi = -1 / 2**6
+        bqmPdbqm = bqm + chi * dbqm
+        self.assertEqual(bqmPdbqm, background_susceptibility_bqm(bqm, chi=chi))
+
+    def test_maximum_pseudolikelihood_bqms(self):
+        """Tests for parameters beyond those applicable to maximum_pseudolikelihood_temparature."""
+        # h1 s1 + h2 s2 + J12 s1 s2; coefficients to be inferred:
+        x = [0.5, -0.4, 0.3]
+        # bqms on space {-1,1}^2, defined without coefficients
+        bqms = [
+            dimod.BinaryQuadraticModel("SPIN").from_ising(
+                {j: (j == i) for j in range(2)}, {}
+            )
+            for i in range(2)
+        ] + [dimod.BinaryQuadraticModel("SPIN").from_ising({}, {(0, 1): 1})]
+        bqm = sum(
+            [-x * bqm for bqm, x in zip(bqms, x)]
+        )  # total bqm weighted by coefficients
+        ss = dimod.ExactSolver().sample(bqm)
+        # In practice, by sampling - for testing exact ratios
+        exact_unnormalized_weights = np.exp(
+            -ss.record.energy + np.min(ss.record.energy)
+        )
+        # A discretized version of weights, to check reproducibility at high precision
+        ss.record.num_occurrences = np.ceil(1000 * exact_unnormalized_weights)
+        # Correct outcome indendent of formatting / detailed root finding spec.
+        beta_by_sampling = None
+        x_by_sampling = None
+        for roo, df, uj, sw in product(
+            [True, False],
+            [True, False],
+            [True, False],
+            [None, exact_unnormalized_weights],
+        ):
+            # Note that, if sample weights are None, we use num_occurrences,
+            # otherwise sample weights are exactly matched to Boltzmann
+            # frequencies (limit of many fair samples).
+            res = maximum_pseudolikelihood(
+                bqms=[bqm],  # Recover minus Inverse temperature -1
+                sampleset=ss,
+                sample_weights=sw,
+                return_optimize_object=roo,
+                degenerate_fields=df,
+                use_jacobian=uj,
+            )
+            if roo:
+                self.assertTrue(res[0].converged)
+                x_ret = res[0].root
+            else:
+                x_ret = res[0]
+            if sw is None:
+                # Sampled distribution, consistent inference:
+                if beta_by_sampling is None:
+                    beta_by_sampling = x_ret
+                    self.assertLess(abs(beta_by_sampling + 1), 1e-2)  # Close
+                else:
+                    self.assertAlmostEqual(x_ret, beta_by_sampling)
+            else:
+                # Exact distribution, exact inference
+                self.assertAlmostEqual(x_ret, -1)
+            if df:
+                # Multidimensional histogramming not (currently) supported
+                with self.assertRaises(ValueError):
+                    maximum_pseudolikelihood(
+                        bqms=bqms,  # Recover h1, h2, J12
+                        sampleset=ss,
+                        sample_weights=exact_unnormalized_weights,
+                        return_optimize_object=roo,
+                        degenerate_fields=df,
+                        use_jacobian=uj,
+                    )
+                continue
+
+            res = maximum_pseudolikelihood(
+                bqms=bqms,  # Recover h1, h2, J12
+                sampleset=ss,
+                sample_weights=sw,
+                return_optimize_object=roo,
+                degenerate_fields=df,
+                use_jacobian=uj,
+            )
+            if roo:
+                self.assertTrue(res[0].success)
+                x_ret = res[0].x
+            else:
+                x_ret = res[0]
+            if sw is None:
+                # Sampled distribution, consistent inference:
+                if x_by_sampling is None:
+                    x_by_sampling = x_ret
+                    self.assertTrue(
+                        all(
+                            abs(a / b - 1) < 1e-2 for a, b in zip(x_by_sampling, x_ret)
+                        ),
+                        1e-2,
+                    )  # Close
+                else:
+                    self.assertTrue(
+                        all(abs(v1 - v2) < 1e-4 for v1, v2 in zip(x_ret, x_by_sampling))
+                    )
+            else:
+                # Exact weights, exact recovery:
+                self.assertTrue(all(abs(v1 - v2) < 1e-4 for v1, v2 in zip(x_ret, x)))
+
+    def test_maximum_pseudolikelihood_instructive_examples(self):
+        """2 degenerate ground states, symmetry is broken by background
+        susceptibility. beta=1 and beta chi=0.02 are recovered given exact
+        sample weights"""
+
+        chi = -0.02
+        dh = 1 / 4
+        bqm = dimod.BinaryQuadraticModel("SPIN").from_ising(
+            {0: dh, 1: -2 * dh, 2: dh}, {(i, i + 1): -1 for i in range(2)}
+        )
+
+        dbqm = background_susceptibility_bqm(bqm)
+        bqmPdbqm = bqm + chi * dbqm
+        ss = dimod.ExactSolver().sample(bqmPdbqm)
+        sample_weights = np.exp(-ss.record.energy + np.min(ss.record.energy))
+        xtup, _ = maximum_pseudolikelihood(
+            bqms=[bqm, dbqm], sampleset=ss, sample_weights=sample_weights
+        )
+        self.assertAlmostEqual(xtup[0], -1)  # unperturbed Hamiltonian
+        self.assertAlmostEqual(xtup[1], -chi)  # susceptibility
+
     def test_maximum_pseudolikelihood_temperature(self):
         # Single variable H = s_i problem with mean energy (-15 + 5)/20 = -0.5
         # 5 measured excitations out of 20.
         # This implies an effective temperature 1/atanh(0.5)
-        site_energy = np.array([2]*5 + [-2]*15)
-        site_names = ['a']
-        for optimize_method in [None,'bisect']:
+        en1 = np.array([2] * 5 + [-2] * 15)
+        for optimize_method in ["bisect", None]:
             T = maximum_pseudolikelihood_temperature(
-                site_energy = (site_energy[:,np.newaxis],site_names),
-                optimize_method = optimize_method)
-            self.assertTrue(type(T) is tuple and len(T)==2)
-            self.assertTrue(np.abs(T[0]-1/np.arctanh(0.5))<1e-8)
-        
+                en1=en1[:, np.newaxis], optimize_method=optimize_method
+            )
+            self.assertTrue(type(T) is tuple and len(T) == 2)
+            self.assertLess(
+                np.abs(T[0] - 1 / np.arctanh(0.5)),
+                1e-8,
+                f"T={1/np.arctanh(0.5)} expected, but T={T[0]}; "
+                f"optimize_method={optimize_method}",
+            )
+
         # Single variable H = s_i problem with mean energy (-5 + 5)/10 = 0
         # This implies an infinite temperature (up to numerical tolerance
         # threshold of scipy optimize.)
-        site_energy = np.array([1]*5 + [-1]*5)
-        T_bracket = [0.1,1]
-        with self.assertWarns(UserWarning) as w:
+        en1 = np.array([1] * 5 + [-1] * 5)
+        T_bracket = [0.1, 1]
+        with self.assertWarns(UserWarning):
             # Returned value should match upper bracket value and
             # throw a warning.
             # Temperature is infinite (excitations and relaxations)
             # are equally likely.
             T = maximum_pseudolikelihood_temperature(
-                site_energy = (site_energy[:,np.newaxis],[1]),
-                T_bracket=T_bracket)
-            self.assertTrue(type(T) is tuple and len(T)==2)
-            self.assertTrue(T[0]==T_bracket[1])
-        
-        # Single variable H = s_i problem with no sample excitations 
-        # This implies zero temperature 
-        # Any bounds on T_bracket should be ignored. 
-        site_energy = np.array([-1]*5)
-        import warnings
+                en1=en1[:, np.newaxis], optimize_method="bisect", T_bracket=T_bracket
+            )
+            self.assertTrue(type(T) is tuple and len(T) == 2)
+            self.assertTrue(T[0] == T_bracket[1])
+
+        # Single variable H = s_i problem with no sample excitations
+        # This implies zero temperature
+        # Any bounds on T_bracket should be ignored.
+        en1 = np.array([-1] * 5)
         with warnings.catch_warnings():
-            #Ignore expected 'out of T_bracket bound' warning:
-            warnings.simplefilter(action='ignore', category=UserWarning)
+            # Ignore expected 'out of T_bracket bound' warning:
+            warnings.simplefilter(action="ignore", category=UserWarning)
             T = maximum_pseudolikelihood_temperature(
-                site_energy = (site_energy[:,np.newaxis],[1]),
-                T_bracket=T_bracket)
-        self.assertTrue(type(T) is tuple and len(T)==2)
-        self.assertTrue(T[0]==0)
+                en1=en1[:, np.newaxis], T_bracket=T_bracket
+            )
+        self.assertTrue(type(T) is tuple and len(T) == 2)
+        self.assertTrue(T[0] == 0)
 
     def test_freezeout_effective_temperature(self):
         # 24mK and 1GHz line up conveniently for T=1.00
-        BsGHz=1
-        TmK=24
-        T = freezeout_effective_temperature(BsGHz,TmK)
-        self.assertTrue(np.round(T*100)==100)
-        
+        BsGHz = 1
+        TmK = 24
+        T = freezeout_effective_temperature(BsGHz, TmK)
+        self.assertTrue(np.round(T * 100) == 100)
+
         # https://docs.dwavesys.com/docs/latest/doc_physical_properties.html
         # Accessed November 12th, 2021
         # Advantage_system4.1 (Advantage): B(s=0.612) = 3.91GHz , T = 15.4mK
         # T_eff = 0.16
         # DW_2000Q_6 (DW2000Q-LN): B(s=0.719) = 6.54GHz , T = 13.5mK
-        # T_eff = 0.086 
-        BsGHz=3.91
-        TmK=15.4
-        T = freezeout_effective_temperature(BsGHz,TmK)
-        self.assertTrue(np.round(T*100)==16)
-        
-        BsGHz=6.54
-        TmK=13.5
-        T = freezeout_effective_temperature(BsGHz,TmK)
-        self.assertTrue(np.round(T*1000)==86)
-        
-        
+        # T_eff = 0.086
+        BsGHz = 3.91
+        TmK = 15.4
+        T = freezeout_effective_temperature(BsGHz, TmK)
+        self.assertTrue(np.round(T * 100) == 16)
+
+        BsGHz = 6.54
+        TmK = 13.5
+        T = freezeout_effective_temperature(BsGHz, TmK)
+        self.assertTrue(np.round(T * 1000) == 86)
+
     def test_fast_effective_temperature(self):
         # Initializing in a ground state, all effective
         # fields must be non-negative.
         sampler = MockDWaveSampler()
-        T, sigma = fast_effective_temperature(sampler=sampler)
-        # Simulated Annealer will tend to return only ground
-        # states, hence temperature 0. But exceptions are
-        # possible. Hence no assertion on value.
-        
+        with warnings.catch_warnings():  # Suppress MockDWaveSampler "no auto_scale" warning
+            warnings.simplefilter(action="ignore", category=UserWarning)
+            T, sigma = fast_effective_temperature(sampler=sampler)
+            # MockDWaveSampler() returns only local minima for high precision
+            # problems (ExactSolver or SteepestDescentSolver)
+            self.assertEqual(T, 0)
 
     def test_bootstrap_errors(self):
-        site_energy = np.array([2]*25 +  [-2]*75)
+        en1 = np.array([2] * 25 + [-2] * 75)
         num_bootstrap_samples = 100
-        
-        T,Tb = maximum_pseudolikelihood_temperature(site_energy = (site_energy[:,np.newaxis],[1]),num_bootstrap_samples = num_bootstrap_samples)
-        
+
+        T, Tb = maximum_pseudolikelihood_temperature(
+            en1=en1[:, np.newaxis], num_bootstrap_samples=num_bootstrap_samples
+        )
+
         # Add test to check bootstrap estimator implementation.
         # T = 1/np.arctanh(0.5). With high probability bootstrapped values
         # are finite and will throw no warnings.
         self.assertTrue(len(Tb) == num_bootstrap_samples)
+
+    def test_sample_weights(self):
+        n = 3
+        bqm = dimod.BinaryQuadraticModel("BINARY").from_qubo(
+            {(i, j): np.random.normal() for i in range(n) for j in range(i, n)}
+        )
+        ss = dimod.ExactSolver().sample(bqm)
+        for Texact in [1, 0.1 + 0.9 * np.random.random()]:
+            sample_weights = np.exp(
+                -1 / Texact * (ss.record.energy + np.min(ss.record.energy))
+            )
+            Ttup, _ = maximum_pseudolikelihood_temperature(
+                bqm=bqm, sampleset=ss, sample_weights=sample_weights
+            )
+            # Previously 1e-8, but less accurate in some environments
+            # tested by circleCI
+            self.assertLess(abs(Ttup - Texact), 1e-2)
+
+    def test_background_susceptibility_ising(self):
+        #
+        n = 3
+        h = np.random.normal(size=n)
+        J = np.random.normal(size=(n, n))
+        J = J + J.transpose() - 2 * np.diag(np.diag(J))
+        dh, dJ, k = background_susceptibility_ising(h, J)
+        self.assertEqual(type(dh), np.ndarray)
+        self.assertEqual(type(dJ), np.ndarray)
+        self.assertEqual(type(k), np.float64)
+        self.assertEqual(dh.shape, h.shape)
+        self.assertEqual(J.shape, dJ.shape)
+        h_dict = {i: h[i] for i in range(n)}
+        J_dict = {(i, j): J[i, j] for i in range(n) for j in range(i + 1, n)}
+        dh_dict, dJ_dict, k = background_susceptibility_ising(h_dict, J_dict)