# Copyright 2022 D-Wave
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
r"""Some notes on the parameter-estimation methods provided below:
- :func:`fluxbias_to_h` and :func:`h_to_fluxbias`
The biases (:math:`h`) equivalent to application of :term:`flux bias` or
vice-versa can be inferred as a function of the normalized anneal fraction,
:math:`s=t/t_a`, by device-specific unit conversion. The needed parameters
for estimation (:math:`\textrm M_{\textrm{AFM}}`, :math:`B(s)`) are
published for online solvers in the :ref:`qpu_solver_properties_specific`
section.
- :func:`.freezeout_effective_temperature`
An effective temperature can be inferred from assuming freeze-out
during the anneal, at :math:`s=t/t_a`, an annealing schedule, and a physical
temperature for the device. The needed device-specific properties are
published for online solvers in the
:ref:`qpu_solver_properties_specific` section.
- :func:`.maximum_pseudolikelihood_temperature`
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. [Chat2007]_ and
[Ray2016]_ describe features of the estimator in application to equilibrium
distribution drawn from binary quadratic models and non-equilibrium
distributions generated by annealing.
- :func:`maximum_pseudolikelihood`
Maximum pseudo-likelihood can be used to infer multiple parameters of
an exponential (Boltzmann) distribution given fair samples. The method
supports inference for probability distributions structured as
:math:`P(s) = \exp(-\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 :math:`h_1`, :math:`h_2` and :math:`J_{12}` 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 method is not optimized for inferring many local terms, but rather
parameters controlling the energy of many degrees of freedom (spins).
An example of practical importance is the case where :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 :math:`h/J`
asymmetry or background susceptibility.
"""
import warnings
import numpy as np
import networkx as nx
import dimod
from scipy import optimize
from typing import Tuple, Union, Optional, Literal, List
from collections import defaultdict
__all__ = [
"background_susceptibility_bqm",
"background_susceptibility_ising",
"effective_field",
"fast_effective_temperature",
"fluxbias_to_h",
"freezeout_effective_temperature",
"h_to_fluxbias",
"Ip_in_units_of_B",
"maximum_pseudolikelihood",
"maximum_pseudolikelihood_temperature",
]
[docs]
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.
Depending on the value you set for the ``current_state_energy`` parameter,
the effective field is one of the following:
* ``False``: 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).
* ``True``: The energy difference from flipping the variable state from
its current value (e.g., from -1 to 1, for an :term:`Ising` model, or 0
to 1, for a :term:`QUBO`). 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) for a
sample, that sample is a local minima (maxima).
As shown in the :ref:`qpu_qubo_ising_transformations` section, any binary
quadratic model (BQM) can be converted to an Ising model using
.. math::
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 effective field on variable :math:`i`, :math:`f_i(s)`, is
then defined as
* ``current_state_energy == False``:
.. math::
f_i(s) = h_i + \sum_j J_{i, j} s_j
* ``current_state_energy == True``:
.. math::
f_i(s) = 2 s_i [h_i + \sum_j J_{i, j} s_j]
Args:
bqm:
A dimod binary quadratic model.
samples:
Either a :class:`~dimod.SampleSet` or a ``samples_like`` object (an
extension of the :std:doc:`NumPy <numpy:index>`
:term:`array_like <numpy:array_like>` structure); see
:func:`~dimod.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 with a state assignment of 1). When set to ``True``,
returns the energy lost in flipping the value of each variable. Note
that ``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:
A tuple of the effective fields and the variable labels. Effective
fields are returned as a :obj:`numpy.ndarray` object and variable
labels are returned as a list. Rows index samples and columns index
variables in the order returned by variable labels.
Examples:
For a ferromagnetic Ising chain :math:`H = - 0.5 \sum_i s_i s_{i+1}`
and for a ground state sample (all +1 values), the energy lost when
flipping any spin is equal to the number of couplers frustrated: -2 for
most the chain (variables 1,2,..,N-2), and -1 at the ends (variables 0
and N-1).
>>> import dimod
>>> 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)})
>>> 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)
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))
labels = bqm.variables
else:
samples, labels = dimod.sampleset.as_samples(samples)
if bqm.vartype is dimod.BINARY:
bqm = bqm.change_vartype(dimod.SPIN, inplace=False)
samples = 2 * samples - 1
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))
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])
if current_state_energy is True:
effective_fields = 2 * samples * effective_fields
return (effective_fields, labels)
[docs]
def background_susceptibility_ising(
h: Union[np.ndarray, dict], J: Union[np.ndarray, dict]
) -> Tuple:
r"""Calculate the biases and couplings due to background susceptibility.
:ref:`Background susceptibility <qpu_ice_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
:term:`Ising` model in matrix vector notation. Assuming :math:`J` is a real,
symmetric (0 on the diagonal) matrix, :math:`h` and :math:`s` are column
vectors, and :math:`\chi` is a small real value, then the Ising model is
is defined as,
.. math::
H = k + (h + \chi J h)' s + 1/2 s' (J + \chi J^2) s,
where the apostrophe (') denotes transpose and matrix multiplication
applies.
The constant (:math:`k`, irrelevant to sampled distributions) is defined as
:math:`k = - \frac{\chi}{2} \textrm{Trace}[J^2]`.
This function returns the perturbative part
:math:`H(s) = k + h' J s + \frac{1}{2} s' (J + \chi J^2) s`.
Args:
h: Linear terms (fields/biases) in the Hamiltonian.
J: Quadratic terms (couplings) in the Hamiltonian.
Returns:
Couplings and scalar constant as a tuple. If :math:`h` and :math:`J` are
of type :obj:`numpy.ndarray`, returned fields and couplings are too;
otherwise returns a tuple of dictionaries.
Examples:
This example calculates the background-susceptibility corrections, shown
symbolically in the
:ref:`Background Susceptibility <qpu_ice_background_susceptibility>`
section, for a three-qubit BQM with the following values:
* Biases :math:`h_1=0.3, h_2=0.8, h_3=-0.25`
* Couplings :math:`J_{1,2}=1.2, J_{2,3}=-0.4`
* :math:`\chi=-0.01`
>>> import dimod
>>> from dwave.system.temperatures import background_susceptibility_ising
...
>>> h1 = 0.3; h2 = 0.8; h3 = -0.25
>>> J12 = 1.2; J23 = -0.4
>>> h, J, _ = background_susceptibility_ising(
... {'1': h1, '2': h2, '3': h3},
... {('1', '2'): J12, ('2', '3'): J23})
>>> print(round(h['2'], 2))
0.46
The formula shown in the
:ref:`Background Susceptibility <qpu_ice_background_susceptibility>`
section for qubit 2, is
:math:`h_2 + h_1 \chi J_{1,2} + h3 \chi J_{2,3}`. The perturbation part
can be rearranged as
:math:`h_1 \chi J_{1,2} + h3 \chi J_{2,3} = \chi (h_1 J_{1,2} + h3 J_{2,3})`
and using the values of this example, setting :math:`\chi = 1` for
convenience, gives
:math:`(h_1 J_{1,2} + h3 J_{2,3}) = 0.3*1.2 - 0.25*(-0.4) = 0.459999`.
"""
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
[docs]
def background_susceptibility_bqm(bqm: dimod.BinaryQuadraticModel, chi: Optional[float] = None):
r"""Create the binary quadratic model for the background susceptibility correction.
:ref:`Background susceptibility <qpu_ice_background_susceptibility>` can be
treated as a perturbative correction to the programmed Hamiltonian:
.. math::
bqm = bqm \textrm{(no background susceptibility)} +
\chi \times dbqm \textrm{(corrections)}
Args:
bqm: A :ref:`dimod <index_dimod>` :term:`binary quadratic model` (BQM)
describing the programmed Hamiltonian.
chi: Scale of background susceptibility correction, :math:`\chi`. If
``chi`` is ``None``, :math:`dbqm` is returned; otherwise the
perturbed Hamiltonian is returned.
Returns:
:class:`~dimod.binary.binary_quadratic_model.BinaryQuadraticModel`:
The given BQM with values corrected for the given background
susceptibility.
Examples:
This example calculates the background-susceptibility corrections, shown
symbolically in the
:ref:`Background Susceptibility <qpu_ice_background_susceptibility>`
section, for a three-qubit BQM with the following values:
* Biases :math:`h_1=0.3, h_2=0.8, h_3=-0.25`
* Couplings :math:`J_{1,2}=1.2, J_{2,3}=-0.4`
* :math:`\chi=-0.01`
>>> import dimod
>>> from dwave.system.temperatures import background_susceptibility_bqm
...
>>> h1 = 0.3; h2 = 0.8; h3 = -0.25
>>> J12 = 1.2; J23 = -0.4
>>> chi = -0.01
>>> bqm = dimod.BQM.from_ising(
... {'1': h1, '2': h2, '3': h3},
... {('1', '2'): J12, ('2', '3'): J23})
>>> bqm_chi = background_susceptibility_bqm(bqm, chi)
>>> print(round(bqm_chi.linear['1'], 2))
0.29
The formula shown in the
:ref:`Background Susceptibility <qpu_ice_background_susceptibility>`
section for qubit 1, is :math:`h_1 + h_2 \chi J_{1,2}`; using the values
of this example gives
:math:`h_1 + h_2 \chi J_{1,2} = 0.3 + 0.8*(-0.01)*1.2 = 0.2904`.
"""
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
[docs]
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"""Return 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 (:term:`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
temperature efficiently (in compute time and number of samples). If the BQM
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
.. math::
0 = \sum_i \sum_{s \in S} f_i(s) \exp(f_i(s)/T),
where :math:`f` is the energy lost in flipping spin :math:`i` against its
current assignment (the effective field).
The problem is a convex root-solving problem, and is solved with SciPy's
``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
temperature via a late-anneal freeze-out phenomena.
Args:
bqm (:class:`~dimod.binary.BinaryQuadraticModel`, optional):
Binary quadratic model describing sample distribution.
If ``sampleset`` and ``en1`` are not provided, then by default
100 samples are drawn using
:class:`~dwave.system.samplers.DWaveSampler`.
sampleset (:class:`~dimod.SampleSet` or ``samples_like`` as described in the :func:`~dimod.as_samples` function, optional):
A set of samples, assumed to be fairly sampled from a Boltzmann
distribution characterized by ``bqm``.
en1 (:obj:`numpy.ndarray` object, optional):
Effective fields as an :obj:`numpy.ndarray` object (site labels not
required). If not provided, derived from the ``bqm`` and
``sampleset`` parameters. First dimension indexes samples and second
dimension indexes sites. Ordering does not matter but should be
consistent with the ``sample_weights`` parameter.
num_bootstrap_samples (int, optional, default=0):
Number of bootstrap estimators to calculate. Currently supported for
samplesets only with uniform ``sample_weights``. 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
(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=None):
Optimize method used by the SciPy ``root_scalar`` algorithm. The
default method (type of solver) works well under default operation;
the 'bisect' method 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)):
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.
sample_weights (:class:`numpy.ndarray`, optional):
A set of weights for the samples. If sampleset is of
type :obj:`~dimod.SampleSet`, this defaults to
``sampleset.record.num_occurrences``; otherwise uniform weighting is
the default.
Returns:
Tuple: The optimal parameters and a list of bootstrapped estimators
(``T_estimate``, ``T_bootstrap_estimates``):
* ``T_estimate``: Temperature estimate
* ``T_bootstrap_estimates``: List of bootstrap estimates
Examples:
This example drawa samples from a D-Wave quantum computer for a large
spin-glass problem (random couplers :math:`J`, zero external field
:math:`h`) and establishes 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
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
ground states).
>>> import dimod
>>> from dwave.system.temperatures import maximum_pseudolikelihood_temperature
>>> from dwave.system import DWaveSampler
>>> from random import random
...
>>> 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)
>>> print('Effective temperature ',T) # doctest: +SKIP
Effective temperature 0.24066488780293813
See also:
[Chat2007]_ and [Ray2016]_.
"""
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):
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:
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
[docs]
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:
r"""Maximimum pseudolikelihood estimator for exponential models.
Uses the SciPy ``optimize`` submodule to solve the maximum pseudolikelihood
problem of weight estimation for an exponential model with the exponent
defined by a weighted sum of binary quadratic models (:term:`BQM`).
Assuming the probability of states :math:`s` is given by an exponential
model, :math:`P(s)=\frac{1}{Z} \textrm{exp}(- x_i H_i(s))`, for a given list
of functions (BQMs), estimate :math:`x`. The exponential model is also
called a Boltzmann distribution.
This function assumes the BQMs are dense (a function of all, or most of,
the variables), although operations on sparse BQMs such as
:math:`H(s)=h_i s_i` or :math:`H(s)=J_{ij}s_i s_j` are technically allowed.
If all provided samples define local minima or maxima (all energies for the
``en1`` parameter are same sign), the limit of :math:`\pm \textrm{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 quantum annealing, so that this
pathological case can be realized in applications of interest.
Common reasons for failures to infer parameters include:
- Too few samples (insufficient to resolve parameters).
- A sampleset that is singular or insensitive with respect to some
parameter (e.g. local minima or plateus only).
- BQMs that are too closely related, or weakly dependent on, the sampleset
variability (e.g. nearly collinear).
Args:
en1: Effective fields as an :class:`numpy.ndarray` (site labels not
required). If not provided, derived from the ``bqms`` and
``sampleset`` parameters. First dimension indexes samples and second
dimension indexes sites. Ordering does not matter but should be
consistent with the ``sample_weights`` parameter.
bqms: List of binary quadratic models :math:`[H_i]` describing the
sample distribution as
:math:`P(s) \approx \textrm{exp}(sum_i x_i H_i(s))` for unknown
model parameters :math:`x`.
sampleset: Set of samples, as a :class:`~dimod.SampleSet` or a
``samples_like`` object (an extension of the
:std:doc:`NumPy <numpy:index>` :term:`array_like <numpy:array_like>`
structure); see :func:`~dimod.as_samples`.
num_bootstrap_samples: Number of bootstrap estimators to calculate.
Bootstrapped estimates can be used to reliably estimate variance and
bias if samples are uncorrelated. Currently supported for samplesets
only with uniform ``sample_weights``; 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 the SciPy ``root_scalar`` algorithm. The
default method (type of solver) works well under default operation;
the 'bisect' method can be numerically more stable for the scalar
case (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 your ``sampleset``
parameter is of type :obj:`~dimod.SampleSet`, defaults to
``sampleset.record.num_occurrences``; otherwise uniform weighting is
the default.
return_optimize_object: When ``True``, and if the SciPy ``optimize``
submodule 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 (as the
default) for single-parameter esimation, ``len(bqms)=1``; for
multi-parameter estimation only value ``False`` is supported.
use_jacobian: By default (``use_jacobian=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 this parameter to False.
Returns:
Tuple: Optimal parameters and a list of bootstrapped estimates
(``x_estimate``, ``x_bootstrap_estimates``):
* ``x_estimate``: parameter estimates
* ``x_bootstrap_estimates``: NumPy array of bootstrap estimators
Examples:
This example builds upon the
:func:`.maximum_pseudolikelihood_temperature` example.
Draw samples from a D-Wave quantum computer for a large spin-glass
problem (random couplers :math:`J`, zero external field :math:`h`).
Establish a temperature and estimate of the background susceptibility
:math:`\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
:math:`\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:
:func:`.maximum_pseudolikelihood_temperature`, [Chat2007]_ and [Ray2016]_.
"""
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(
"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(
"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:
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:
# 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:
if len(np.unique(sample_weights)) != 1:
raise ValueError(
"Bootstraps require uniform sample_weights (num_occurrences)"
)
prng = np.random.RandomState(seed)
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
[docs]
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,
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
: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.
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`.
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.
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.
``MAFM`` is ignored when ``Ip`` is specified.
units_Ip:
Units in which the persistent current, ``Ip``, is specified.
Allowed values are ``'uA'`` (microamps) and ``'A'`` (amps)
units_B:
Units in which the schedule ``B`` is specified. Allowed values
are ``'GHz'`` (gigahertz) and ``'J'`` (Joules).
units_MAFM:
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)`.
Examples:
This example calculates the persistent current of an |adv2| QPU at its
quantum critical point, as provided in the
:ref:`qpu_solver_properties_specific` section.
>>> from dwave.system.temperatures import Ip_in_units_of_B
>>> ip = Ip_in_units_of_B(B=2.308, units_B="GHz", MAFM=2.113, units_MAFM="pH")
"""
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":
B_multiplier = 1
else:
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)
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
[docs]
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
: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
freeze-out (quasi-static) hypothesis.
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`.
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.
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.
``MAFM`` is ignored when ``Ip`` is specified.
units_Ip:
Units in which the persistent current, ``Ip``, is specified.
Allowed values are ``'uA'`` (microamps) and ``'A'`` (amps)
units_B:
Units in which the schedule ``B`` is specified. Allowed values
are ``'GHz'`` (gigahertz) and ``'J'`` (Joules).
units_MAFM:
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
``h`` values.
Examples:
See the :ref:`qpu_config_emulate_with_fbo` section.
"""
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
[docs]
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.
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
freeze-out (quasi-static) hypothesis.
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:
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`.
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.
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.
``MAFM`` is ignored when ``Ip`` is specified.
units_Ip:
Units in which the persistent current, ``Ip``, is specified.
Allowed values are ``'uA'`` (microamps) and ``'A'`` (amps)
units_B:
Units in which the schedule ``B`` is specified. Allowed values
are ``'GHz'`` (gigahertz) and ``'J'`` (Joules).
units_MAFM:
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.
Examples:
See the :ref:`qpu_config_emulate_with_fbo` section.
"""
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
[docs]
def freezeout_effective_temperature(
freezeout_B: float, temperature: float, units_B: str = "GHz", units_T: str = "mK"
) -> float:
r"""Calculate an effective temperature for a freezeout point.
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
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, and :math:`s` is the normalized anneal time :math:`s = t/t_a`
(in :math:`[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
.. math::
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)`, and dynamics stop
abruptly at :math:`s=s^*`, the distribution of returned samples is well
described by a Boltzmann distribution,
.. math::
P(x) = \exp(- B(s^*) R E_{Ising}(x) / 2 k_B T),
where :math:`T` is the physical temperature and :math:`k_B` is the Boltzmann
constant. :math:`R` is a Hamiltonian rescaling factor:
* :math:`R=1` if the :term:`QPU` is operated with the
:ref:`parameter_qpu_auto_scale` set to ``False``.
* :math:`R \neq 1` if ``auto_scale=True`` (default) and
*additional scaling factors must be accounted for*.
See the :ref:`qpu_annealing` section for operation details of D-Wave
annealing quantum computers.
This 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)` and
:math:`B(s)`, and single-qubit freeze-out (:math:`s^*`, for simple uncoupled
Hamltonians), are properties reported for each device in the
:ref:`qpu_solver_properties_specific` section. These values (typically
specified in mK and GHz) allow 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 :math:`s` or piecemeal. Large problems may have
slow dynamics at small values of :math:`s`, in which cases :math:`A(s)`
cannot be ignored as a contributing factor to the distribution.
Args:
freezeout_B (float):
:math:`B(s^*)`, the problem Hamiltonian energy scale at freeze-out.
temperature (float):
:math:`T`, the physical temperature of the quantum computer.
units_B (string, 'GHz'):
Units in which ``freezeout_B`` is specified. Allowed values:
'GHz' (Giga-Hertz) and 'J' (Joules).
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.
Examples:
This example uses the
:ref:`published parameters <qpu_solver_properties_specific>`
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
>>> T = freezeout_effective_temperature(freezeout_B = 3.91, temperature = 15.4)
>>> print('Effective temperature at single qubit freeze-out is', T) # doctest: +ELLIPSIS
Effective temperature at single qubit freeze-out is 0.164...
See also:
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
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 *= h * 1e9
elif units_B == "J":
pass
else:
raise ValueError("Units must be 'J' (Joules) or 'mK' (milli-Kelvin)")
if units_T == "mK":
temperature = temperature * 1e-3
elif units_T == "K":
pass
else:
raise ValueError("Units must be 'K' (Kelvin) or 'mK' (milli-Kelvin)")
kB = 1.3806503e-23 # Joules/Kelvin
return 2 * temperature * kB / freezeout_B
[docs]
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"""Estimate 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.
Use the :func:`.maximum_pseudolikelihood_temperature` function for greater
control of the problem Hamiltonian or more-general samplers.
Args:
sampler (:class:`dimod.Sampler`, optional, default=\ :class:`~dwave.system.samplers.DWaveSampler`):
A dimod sampler.
num_reads (int, optional):
Number of reads to use. Defaults to 100 when not specified in the
``sampler_params`` parameter.
seed (int, optional):
Seeds the problem-generation process, allowing reproducibility from
pseudo-random samplers.
h_range (float, default = :math:`[-1/6.1,1/6.1]`):
Determines the range of external fields probed for temperature
inference. Default is based on D-Wave's Advantage processor, where
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
and other nonidealities such as precision limitations.
nodelist (list, optional):
Variables included in the estimator. Defaults to the node list
(``DWaveSampler.nodelist``) of the structured sampler (:term:`QPU`)
being used.
sampler_params (dict, optional):
Any additional non-default sampler parameters. If ``num_reads`` is
provided, must be compatible with ``num_reads`` parameter.
optimize_method (str, optional):
Optimization method used by the SciPy ``root_scalar`` algorithm. The
default method works well under default operation; the 'bisect'
method can be numerically more stable when operated without
defaults.
num_bootstrap_samples (int, optional, default=0):
Number of bootstrap samples to use for estimation of the standard
error. By default no bootstrapping is performed and the standard
error defaults to 0.
Returns:
Tuple[float, float]:
Effective temperature describing single-qubit problems in an
external field, and a standard error (:math:`\pm 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 ``num_reads``
specified in the ``sampler_params`` parameter.
See also:
* The :func:`.freezeout_effective_temperature` function may be used,
in combination with published device values, to estimate
single-qubit freeze-out in approximate agreement with empirical
estimates of this function.
* [Chat2007]_ and [Ray2016]_.
Examples:
This example draws samples from a quantum computer (using the
:class:`~dwave.system.samplers.DWaveSampler` class) and establishs its
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
Effective temperature at freeze-out is 0.21685104745347336
"""
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[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(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:
if sampler_params is None:
sampler_params0 = {}
else:
sampler_params0 = sampler_params.copy()
if num_reads is None:
# 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."
)
else:
sampler_params0["auto_scale"] = False
if num_bootstrap_samples is None:
num_bootstrap_samples = sampler_params0["num_reads"]
sampleset = sampler.sample(bqm, **sampler_params0)
T, Tboot = maximum_pseudolikelihood_temperature(
bqm,
sampleset,
optimize_method=optimize_method,
num_bootstrap_samples=num_bootstrap_samples,
)
if num_bootstrap_samples == 0:
return T, np.float64(0.0)
else:
return T, np.sqrt(np.var(Tboot))