# Copyright 2022 D-Wave Systems Inc.
#
# 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.
"""The following effective temperature estimators 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
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:
https://www.jstor.org/stable/25464568
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.
Necessary device-specific properties are published for online solvers:
https://docs.dwavesys.com/docs/latest/doc_physical_properties.html
"""
import warnings
import numpy as np
import dimod
from scipy import optimize
from typing import Tuple
__all__ = ['effective_field', 'maximum_pseudolikelihood_temperature',
'freezeout_effective_temperature', 'fast_effective_temperature']
[docs]def effective_field(bqm,
samples=None,
current_state_energy=False) -> (np.ndarray,list):
'''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).
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
variable is in a locally excited state.
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
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
else:
.. 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:`.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
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`.
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), the energy lost when flipping
any spin is equal to the number of couplers frustrated: -2 in the center
of the chain (variables 1,2,..,N-2), and -1 at the end (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('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:
#Ising: eff_field = 2*s*(h + J*s)
effective_fields = 2*samples*effective_fields
return (effective_fields,labels)
[docs]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]:
'''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
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 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
optimize.
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 (: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`,
with ``bqm`` defaulted as described.
sampleset (: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.
Derived from the ``bqm`` and ``sampleset`` if not provided.
num_bootstrap_samples (int, optional, default=0):
Number of bootstrap estimators to calculate.
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.
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).
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).
Returns:
Tuple of float and NumPy array:
(T_estimate,T_bootstrap_estimates)
*T_estimate*: a temperature estimate
*T_bootstrap_estimates*: a numpy array of bootstrap estimators
Examples:
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.
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:
https://doi.org/10.3389/fict.2016.00023
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
else:
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
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:
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]
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]
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)
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
if num_bootstrap_samples > 0:
#By bootstrapping with respect to samples we
if root_results is not None:
x0 = root_results.root
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
[docs]def freezeout_effective_temperature(freezeout_B, temperature, units_B = 'GHz', units_T = 'mK') -> float:
'''Provides an effective temperature as a function of freezeout information.
See https://docs.dwavesys.com/docs/latest/c_qpu_annealing.html 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
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
.. 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 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,
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 properties: https://docs.dwavesys.com/docs/latest/doc_physical_properties.html
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.
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, optional, 'GHz'):
Units in which ``freezeout_B`` is specified. Allowed values:
'GHz' (Giga-Hertz) and 'J' (Joules).
units_T (string, optional, '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
published parameters <https://docs.dwavesys.com/docs/latest/doc_physical_properties.html>
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.626068e-34 #J/Hz
freezeout_B = freezeout_B *h
freezeout_B *= 1e9
elif units_B == 'J':
pass
else:
raise ValueException("Units must be 'J' (Joules) "
"or 'mK' (milli-Kelvin)")
if units_T == 'mK':
temperature = temperature * 1e-3
elif units_T == 'K':
pass
else:
raise ValueException("Units must be 'K' (Kelvin) "
"or 'mK' (milli-Kelvin)")
kB = 1.3806503e-23 # J/K
return 2*temperature*kB/freezeout_B
[docs]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]:
'''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.
Args:
sampler (:class:`dimod.Sampler`, optional, default=\ :class:`~dwave.system.samplers.DWaveSampler`):
A dimod sampler.
num_reads (int, optional):
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]):
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
(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.
sampler_params (dict, optional):
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
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 is defaulted to 0.
Returns:
Tuple[float, float]:
The effective temperature describing single qubit problems in an
external field, and a standard error (+/- 1 sigma).
By default the confidence interval is set as 0.
See also:
https://doi.org/10.3389/fict.2016.00023
https://www.jstor.org/stable/25464568
Examples:
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
See also:
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[1] > sampler.properties['h_range'][1]:
raise ValueError('h_range[1] exceeds programmable range')
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)}, {})
#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.")
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))