Repository: dwavesystems/dwave-system
Pull request #513 (2024-03-18 22:44:25)

Conversion from flux_biases (units of Phi0) to h (unitless) and vice-versa
Both flux_biases and h apply time dependent longitudinal fields (fields conjugate to sigma^z in the Hamiltonian description, see documentation).
At given point in the anneal s=t/t_a, h and flux_biases can apply equal longitudinal fields subject to rescaling. If a freeze-out phenomena is assumed at a given s, as relevant to quantum Boltzmann or Boltzmann sampling,  flux_biases can be replaced by h or vice-versa.
'h' is typically unitless, whereas flux has units of Phi0, the conversion is a function of published QPU properties B(s) and M_AFM. This pull request adds code for this conversion. 


jackraymond: @JoelPasvolsky suggested as a reviewer.
JoelPasvolsky: @jackraymond, I started reviewing now but is it possible you have a switched sign for your ``temperatures.h_to_fluxbias()`` function?

```
>>> fb = [0]*qpu.properties['num_qubits']
>>> fb[30] = temperatures.h_to_fluxbias(0.5)

>>> # What I expect is an equivalent to this (most spins down):

>>> sampleset = qpu.sample_ising({30: 0.5}, {}, num_reads=5000, auto_scale=False)
>>> print(sampleset)
  30 energy num_oc.
0 -1   -0.5    4923
1 +1    0.5      77

 >>> # What I get is the opposite (most spins up):

>>> sampleset = qpu.sample_ising({30: 0}, {}, num_reads=5000, auto_scale=False, flux_biases=fb)
>>> print(sampleset)
  30 energy num_oc.
0 -1    0.0     277
1 +1    0.0    4723
```
jackraymond: > @jackraymond, I started reviewing now but is it possible you have a switched sign for your `temperatures.h_to_fluxbias()` function?
> 
> ```
> >>> fb = [0]*qpu.properties['num_qubits']
> >>> fb[30] = temperatures.h_to_fluxbias(0.5)
> 
> >>> # What I expect is an equivalent to this (most spins down):
> 
> >>> sampleset = qpu.sample_ising({30: 0.5}, {}, num_reads=5000, auto_scale=False)
> >>> print(sampleset)
>   30 energy num_oc.
> 0 -1   -0.5    4923
> 1 +1    0.5      77
> 
>  >>> # What I get is the opposite (most spins up):
> 
> >>> sampleset = qpu.sample_ising({30: 0}, {}, num_reads=5000, auto_scale=False, flux_biases=fb)
> >>> print(sampleset)
>   30 energy num_oc.
> 0 -1    0.0     277
> 1 +1    0.0    4723
> ```

Yep - convention for flux_biases is opposite to the convention for h. Minus sign matters!
jackraymond: > Also need to be added to [this file](https://docs.ocean.dwavesys.com/en/stable/docs_system/reference/utilities.html) I'm not sure that having these function under **Temperature Utilities** is optimal.

This repo contains conversions from GHz to kB T , and some other related functions. Not sure there is a better spot. 
The temperature estimators can use this function directly under the freeze-out hypothesis to estimate effective temperature, so the code will interface directly with other functions in the module. 

However, happy to consider other options. 
JoelPasvolsky: > This repo contains conversions from GHz to kB T , and some other related functions.

I wasn't concerned with the source file, I mean that users looking to use flux bias as a substitution for an ``h`` bias might not look at this temperature-related [**doc** section](https://docs.ocean.dwavesys.com/en/stable/docs_system/reference/utilities.html#module-dwave.system.temperatures).
For the docs, should these be bumped out to general [Utilities page](https://docs.ocean.dwavesys.com/en/stable/docs_system/reference/utilities.html#module-dwave.system.utilities) functions?
jackraymond: > > This repo contains conversions from GHz to kB T , and some other related functions.
> 
> I wasn't concerned with the source file, I mean that users looking to use flux bias as a substitution for an `h` bias might not look at this temperature-related [**doc** section](https://docs.ocean.dwavesys.com/en/stable/docs_system/reference/utilities.html#module-dwave.system.temperatures). For the docs, should these be bumped out to general [Utilities page](https://docs.ocean.dwavesys.com/en/stable/docs_system/reference/utilities.html#module-dwave.system.utilities) functions?

I'll create a 'Unit conversion' subsection, unless you have a better suggestion
randomir: @JoelPasvolsky, test failures here are caused by (private) changes in cloud-client 0.11.3. The fix for mock tests that tweak client internals is in https://github.com/dwavesystems/dwave-system/pull/515.

Patch:
diff --git a/docs/reference/utilities.rst b/docs/reference/utilities.rst
index 3006e51..0eedf54 100644
--- a/docs/reference/utilities.rst
+++ b/docs/reference/utilities.rst
@@ -20,8 +20,9 @@ Utilities
 
    coupling_groups
 
-Temperature Utilities
----------------------
+
+Temperature and Unit-Conversion Utilities
+-----------------------------------------
 
 .. automodule:: dwave.system.temperatures
 
@@ -31,6 +32,9 @@ Temperature Utilities
    :toctree: generated/
 
    effective_field
-   maximum_pseudolikelihood_temperature
-   freezeout_effective_temperature
    fast_effective_temperature
+   fluxbias_to_h
+   freezeout_effective_temperature
+   h_to_fluxbias
+   Ip_in_units_of_B
+   maximum_pseudolikelihood_temperature
diff --git a/dwave/system/temperatures.py b/dwave/system/temperatures.py
index f18744c..2495552 100644
--- a/dwave/system/temperatures.py
+++ b/dwave/system/temperatures.py
@@ -12,7 +12,7 @@
 #    See the License for the specific language governing permissions and
 #    limitations under the License.
 
-"""The following effective temperature estimators are provided:
+"""The following effective temperature and bias 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) 
@@ -28,16 +28,22 @@
   Necessary device-specific properties are published for online solvers:
   https://docs.dwavesys.com/docs/latest/doc_physical_properties.html
 
+- 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
+  device-specific unit conversion. The necessary parameters for estimation
+  [Mafm, B(s)] 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
+from typing import Tuple, Union, Optional, Literal
 
 __all__ = ['effective_field', 'maximum_pseudolikelihood_temperature',
-           'freezeout_effective_temperature', 'fast_effective_temperature']
+           'freezeout_effective_temperature', 'fast_effective_temperature',
+           'Ip_in_units_of_B', 'h_to_fluxbias', 'fluxbias_to_h']
 
 def effective_field(bqm,
                     samples=None,
@@ -355,6 +361,219 @@ def maximum_pseudolikelihood_temperature(bqm=None,
     
     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]:
+    """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 :ref:`documented form <sysdocs_gettingstarted:doc_qpu>` 
+    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:`system documentation <sysdocs_gettingstarted:doc_qpu_characteristics>`. 
+            This parameter is ignored when ``Ip`` is specified.
+        
+        MAFM:
+            Mutual inductance, :math:`M_{AFM}`, specified for each quantum 
+            computer in the 
+            :ref:`system documentation <sysdocs_gettingstarted:doc_qpu_characteristics>`. 
+            ``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)`.
+    """
+    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
+
+
+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]:
+    """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:`system documentation <sysdocs_gettingstarted:doc_qpu_characteristics>`. 
+            This parameter is ignored when ``Ip`` is specified.
+        
+        MAFM:
+            Mutual inductance, :math:`M_{AFM}`, specified for each quantum 
+            computer in the 
+            :ref:`system documentation <sysdocs_gettingstarted:doc_qpu_characteristics>`. 
+            ``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.
+    """
+    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]:
+    """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:`system documentation <sysdocs_gettingstarted:doc_qpu_characteristics>`. 
+            This parameter is ignored when ``Ip`` is specified.
+        
+        MAFM:
+            Mutual inductance, :math:`M_{AFM}`, specified for each quantum 
+            computer in the 
+            :ref:`system documentation <sysdocs_gettingstarted:doc_qpu_characteristics>`. 
+            ``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.
+    """
+    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:
     '''Provides an effective temperature as a function of freezeout information.
     
@@ -441,7 +660,7 @@ def freezeout_effective_temperature(freezeout_B, temperature, units_B = 'GHz', u
     
     #Convert units_B to Joules
     if units_B == 'GHz':
-        h = 6.626068e-34 #J/Hz
+        h = 6.62607e-34 #J/Hz
         freezeout_B = freezeout_B *h
         freezeout_B *= 1e9
     elif units_B == 'J':


Test patch:
diff --git a/tests/test_temperatures.py b/tests/test_temperatures.py
index df98e07..1e683ef 100644
--- a/tests/test_temperatures.py
+++ b/tests/test_temperatures.py
@@ -15,16 +15,38 @@
 import unittest
 import numpy as np
 import dimod
+from itertools import product
 
 from dwave.system.temperatures import (maximum_pseudolikelihood_temperature,
                                        effective_field,
                                        freezeout_effective_temperature,
-                                       fast_effective_temperature)
+                                       fast_effective_temperature,
+                                       Ip_in_units_of_B,
+                                       h_to_fluxbias,
+                                       fluxbias_to_h)
 
 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']
+        for uIp, uB, uMAFM in product(uIps, uBs, uMAFMs):
+            _ = 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)
+        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))
+
     def test_effective_field(self):
         # For a simple model of independent spins H = sum_i s_i
         # The effective field is 1 (setting a spin to 1, costs 1 unit of energy,