Repository: dwavesystems/dwave-system
Pull request #530 (2024-08-15 21:14:13)

add LinearAncillaComposite
In scenarios where the h biases are not available or sufficient in range, one can use auxiliary qubits polarized with a large flux bias and coupled to the data qubit. This method:

-  Requires that each qubit is coupled to a non-programmed qubit
-  But, ancilla qubits can be shared between data qubits
-  The magnitude of the bias is limited by the magnitude of the coupler, but in the case of asymetric coupling range, one can use the strongest coupling available and control the bias sign with the flux bias sign

This PR adds a LinearAncillaComposite that implements this technique.

For the reviewers, please take into account future extensions of this method:
- One can use multiple ancilla qubits to extend the h-range and avoid shrinking energy scales of the problem. This is relevant for QUBO --> Ising translation.
- While the current implementation may work on embedded problems as is, we need an embedding-aware version that relaxes the condition of each qubit being coupled to an ancilla to each chain being coupled to at least one ancilla.

Please educate me on how to use the correct Sphynx syntax


Patch:
diff --git a/.circleci/config.yml b/.circleci/config.yml
index 9a4aba0..271e28d 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -17,13 +17,11 @@ jobs:
     parameters:
       python-version:
         type: string
-      dependency-specifiers:
-        type: string
       integration-test-python-version:
         type: string
 
     docker:
-      - image: python:<< parameters.python-version >>
+      - image: cimg/python:<< parameters.python-version >>
 
     steps:
       - checkout
@@ -35,6 +33,7 @@ jobs:
             python -m venv env
             . env/bin/activate
             pip install -U pip setuptools
+            pip install .     # install dwave-system to break cycle created by dwave-drivers in requirements
             pip install -r requirements.txt -r tests/requirements.txt
 
       - run: &run-python-tests
@@ -56,16 +55,7 @@ jobs:
             . env/bin/activate
             pip install -U pip setuptools
             pip install -r tests/requirements.txt
-            pip install .
-
-      - when:
-          condition: << parameters.dependency-specifiers >>
-          steps:
-            - run:
-                name: sweep dependencies
-                command: |
-                  . env/bin/activate
-                  pip install << parameters.dependency-specifiers >>
+            pip install .[drivers] -I --no-cache-dir --extra-index-url https://pypi.dwavesys.com/simple
 
       - run: *run-python-tests
 
@@ -74,7 +64,7 @@ jobs:
 
   test-doctest:
     docker:
-      - image: python:3.12
+      - image: cimg/python:3.9
 
     steps:
       - checkout
@@ -84,6 +74,7 @@ jobs:
           command: |
             python -m venv env
             . env/bin/activate
+            pip install .     # install dwave-system to break cycle created by dwave-drivers in requirements
             pip install -r requirements.txt -r docs/requirements.txt
 
       - run:
@@ -106,7 +97,9 @@ jobs:
         type: string
       xcode:
         type: string
-        default: "16.0.0"
+        default: "14.3.0"
+
+    resource_class: macos.m1.medium.gen1
 
     macos:
       xcode: << parameters.xcode >>
@@ -144,6 +137,7 @@ jobs:
             python -m venv env
             . env/bin/activate
             pip install -U pip
+            pip install .     # install dwave-system to break cycle created by dwave-drivers in requirements
             pip install -r requirements.txt -r tests/requirements.txt
 
       - run: *run-python-tests
@@ -179,6 +173,7 @@ jobs:
           command: |
             env\Scripts\activate.ps1
             pip install -U pip
+            pip install .     # install dwave-system to break cycle created by dwave-drivers in requirements
             pip install -r requirements.txt -r tests\requirements.txt
 
       - run:
@@ -189,7 +184,7 @@ jobs:
 
   deploy:
     docker:
-      - image: python:3.12
+      - image: cimg/python:3.9
 
     steps:
       - checkout
@@ -199,8 +194,9 @@ jobs:
           command: |
             python -m virtualenv env
             . env/bin/activate
+            pip install .     # install dwave-system to break cycle created by dwave-drivers in requirements
             pip install -r requirements.txt
-            pip install twine wheel setuptools
+            pip install twine wheel
 
       - run:
           name: verify version matches tag
@@ -224,18 +220,11 @@ workflows:
   test:
     jobs:
       - test-linux:
-          name: test-linux-<< matrix.python-version >> | << matrix.dependency-specifiers >>
+          name: test-linux-<< matrix.python-version >>
           matrix:
             parameters:
-              python-version: &python-versions ["3.9", "3.10", "3.11", "3.12", &latest-python "3.13"]
-              dependency-specifiers:
-                - "dwave-cloud-client~=0.12.0"
-                - "dwave-cloud-client~=0.13.0"
+              python-version: &python-versions ["3.8", "3.9", "3.10", "3.11", &latest-python "3.12"]
               integration-test-python-version: &integration-python-versions [*latest-python]
-            exclude:
-              - python-version: "3.13"
-                dependency-specifiers: "dwave-cloud-client~=0.12.0"
-                integration-test-python-version: "3.13"
 
       - test-osx:
           name: test-osx-<< matrix.python-version >>
diff --git a/.gitignore b/.gitignore
index 3f155ac..2057891 100644
--- a/.gitignore
+++ b/.gitignore
@@ -106,4 +106,4 @@ ENV/
 *.sublime-project
 *.sublime-workspace
 
-generated/
\ No newline at end of file
+generated/
diff --git a/README.rst b/README.rst
index d2f482a..63400b9 100644
--- a/README.rst
+++ b/README.rst
@@ -40,10 +40,23 @@ Installation
 
     pip install dwave-system
 
+**Installation from PyPI with drivers:**
+
 .. note::
-    As of ``dwave-system`` 1.28.0, support for ``dwave-drivers`` is removed (it
-    was used for calibration of qubits in chains via ``VirtualGraphComposite``,
-    but it's no longer required due to improved calibration of newer QPUs).
+    Prior to v0.3.0, running ``pip install dwave-system`` installed a driver dependency called ``dwave-drivers``
+    (previously also called ``dwave-system-tuning``). This dependency has a restricted license and has been made optional
+    as of v0.3.0, but is highly recommended. To view the license details:
+
+    .. code-block:: python
+
+        from dwave.drivers import __license__
+        print(__license__)
+
+To install with optional dependencies:
+
+.. code-block:: bash
+
+    pip install dwave-system[drivers] --extra-index-url https://pypi.dwavesys.com/simple
 
 **Installation from source:**
 
@@ -52,6 +65,12 @@ Installation
     pip install -r requirements.txt
     python setup.py install
 
+Note that installing from source installs ``dwave-drivers``. To uninstall the proprietary components:
+
+.. code-block:: bash
+
+    pip uninstall dwave-drivers
+
 .. installation-end-marker
 
 
diff --git a/docs/reference/composites.rst b/docs/reference/composites.rst
index a14d837..ae3d161 100644
--- a/docs/reference/composites.rst
+++ b/docs/reference/composites.rst
@@ -268,6 +268,41 @@ Methods
    VirtualGraphComposite.sample_ising
    VirtualGraphComposite.sample_qubo
 
+
+
+Linear Bias
+===========
+
+Composite for using auxiliary qubits to bias problem qubits.
+
+
+LinearAncillaComposite
+-----------------------
+
+.. autoclass:: LinearAncillaComposite
+
+Properties
+~~~~~~~~~~
+
+.. autosummary::
+   :toctree: generated/
+
+   LinearAncillaComposite.child
+   LinearAncillaComposite.children
+   LinearAncillaComposite.parameters
+   LinearAncillaComposite.properties
+
+Methods
+~~~~~~~
+
+.. autosummary::
+   :toctree: generated/
+
+   LinearAncillaComposite.sample
+   LinearAncillaComposite.sample_ising
+   LinearAncillaComposite.sample_qubo
+
+
 Reverse Anneal
 ==============
 
@@ -326,3 +361,4 @@ Methods
    ReverseAdvanceComposite.sample
    ReverseAdvanceComposite.sample_ising
    ReverseAdvanceComposite.sample_qubo
+
diff --git a/dwave/embedding/transforms.py b/dwave/embedding/transforms.py
index 3f31b0e..235cfd8 100644
--- a/dwave/embedding/transforms.py
+++ b/dwave/embedding/transforms.py
@@ -356,8 +356,8 @@ def embed_bqm(source_bqm, embedding=None, target_adjacency=None,
         >>> embedding = {'a': {0}, 'b': {1}, 'c': {2, 3}}
         >>> # Embed and show the chain strength
         >>> target_bqm = dwave.embedding.embed_bqm(bqm, embedding, target)
-        >>> print(target_bqm.quadratic[(2, 3)])
-        -1.999...
+        >>> target_bqm.quadratic[(2, 3)]
+        -1.9996979771955565
         >>> print(target_bqm.quadratic)  # doctest: +SKIP
         {(0, 1): 1.0, (0, 3): 1.0, (1, 2): 1.0, (2, 3): -1.9996979771955565}
 
diff --git a/dwave/system/cache/__init__.py b/dwave/system/cache/__init__.py
new file mode 100644
index 0000000..3a63283
--- /dev/null
+++ b/dwave/system/cache/__init__.py
@@ -0,0 +1,16 @@
+# Copyright 2018 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.
+
+from dwave.system.cache.cache_manager import *
+from dwave.system.cache.database_manager import *
diff --git a/dwave/system/cache/cache_manager.py b/dwave/system/cache/cache_manager.py
new file mode 100644
index 0000000..27095c5
--- /dev/null
+++ b/dwave/system/cache/cache_manager.py
@@ -0,0 +1,59 @@
+# Copyright 2018 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.
+
+import os
+
+import homebase
+
+from dwave.system.package_info import __version__
+
+__all__ = ['cache_file']
+
+APPNAME = 'dwave-virtual-graph-cache'
+"""The application name is used to determine the cache location."""
+
+APPAUTHOR = 'dwave-systems'
+"""The application author is used to determine the cache location."""
+
+DATABASENAME = 'virtual-graph_v%s.db' % __version__
+"""The name for the sqlite database itself. Based on the version of the package."""
+
+
+def cache_file(app_name=APPNAME, app_author=APPAUTHOR, filename=DATABASENAME):
+    """Returns the filename (including path) for the data cache.
+
+    The path will depend on the operating system, certain environmental
+    variables and whether it is being run inside a virtual environment.
+    See `homebase <https://github.com/dwavesystems/homebase>`_.
+
+    Args:
+        app_name (str, optional): The application name.
+            Default is given by :obj:`.APPNAME`.
+        app_author (str, optional): The application author. Default
+            is given by :obj:`.APPAUTHOR`.
+        filename (str, optional): The name of the database file.
+            Default is given by :obj:`DATABASENAME`.
+
+    Returns:
+        str: The full path to the file that can be used as a cache.
+
+    Notes:
+        Creates the directory if it does not already exist.
+
+        If run inside of a virtual environment, the cache will be stored
+        in `/path/to/virtualenv/data/app_name`
+
+    """
+    user_data_dir = homebase.user_data_dir(app_name=app_name, app_author=app_author, create=True)
+    return os.path.join(user_data_dir, filename)
diff --git a/dwave/system/cache/database_manager.py b/dwave/system/cache/database_manager.py
new file mode 100644
index 0000000..321ab4e
--- /dev/null
+++ b/dwave/system/cache/database_manager.py
@@ -0,0 +1,604 @@
+# Copyright 2018 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.
+
+"""Utilities for accessing the sqlite cache.
+
+Note:
+    All methods assume that nodes are integer-labeled.
+
+"""
+import sqlite3
+import json
+import os
+import struct
+import base64
+import datetime
+
+from dwave.system.cache.cache_manager import cache_file
+from dwave.system.cache.schema import schema
+from dwave.system.exceptions import MissingFluxBias
+
+
+def cache_connect(database=None):
+    """Returns a connection object to a sqlite database.
+
+    Args:
+        database (str, optional): The path to the database the user wishes
+            to connect to. If not specified, a default is chosen using
+            :func:`.cache_file`. If the special database name ':memory:'
+            is given, then a temporary database is created in memory.
+
+    Returns:
+        :class:`sqlite3.Connection`
+
+    """
+    if database is None:
+        database = cache_file()
+
+    if os.path.isfile(database):
+        # just connect to the database as-is
+        conn = sqlite3.connect(database)
+    else:
+        # we need to populate the database
+        conn = sqlite3.connect(database)
+        conn.executescript(schema)
+
+    with conn as cur:
+        # turn on foreign keys, allows deletes to cascade.
+        cur.execute("PRAGMA foreign_keys = ON;")
+
+    conn.row_factory = sqlite3.Row
+
+    return conn
+
+
+def insert_chain(cur, chain, encoded_data=None):
+    """Insert a chain into the cache.
+
+    Args:
+        cur (:class:`sqlite3.Cursor`):
+            An sqlite3 cursor. This function is meant to be run within a :obj:`with` statement.
+
+        chain (iterable):
+            A collection of nodes. Chains in embedding act as one node.
+
+        encoded_data (dict, optional):
+            If a dictionary is provided, it will be populated with the serialized data. This is
+            useful for preventing encoding the same information many times.
+
+    Notes:
+        This function assumes that the nodes in chain are index-labeled.
+
+    """
+    if encoded_data is None:
+        encoded_data = {}
+
+    if 'nodes' not in encoded_data:
+        encoded_data['nodes'] = json.dumps(sorted(chain), separators=(',', ':'))
+    if 'chain_length' not in encoded_data:
+        encoded_data['chain_length'] = len(chain)
+
+    insert = "INSERT OR IGNORE INTO chain(chain_length, nodes) VALUES (:chain_length, :nodes);"
+
+    cur.execute(insert, encoded_data)
+
+
+def iter_chain(cur):
+    """Iterate over all of the chains in the database.
+
+    Args:
+        cur (:class:`sqlite3.Cursor`):
+            An sqlite3 cursor. This function is meant to be run within a :obj:`with` statement.
+
+    Yields:
+        list: The chain.
+
+    """
+    select = "SELECT nodes FROM chain"
+    for nodes, in cur.execute(select):
+        yield json.loads(nodes)
+
+
+def insert_system(cur, system_name, encoded_data=None):
+    """Insert a system name into the cache.
+
+    Args:
+        cur (:class:`sqlite3.Cursor`):
+            An sqlite3 cursor. This function is meant to be run within a :obj:`with` statement.
+
+        system_name (str):
+            The unique name of a system
+
+        encoded_data (dict, optional):
+            If a dictionary is provided, it will be populated with the serialized data. This is
+            useful for preventing encoding the same information many times.
+
+    """
+    if encoded_data is None:
+        encoded_data = {}
+
+    if 'system_name' not in encoded_data:
+        encoded_data['system_name'] = system_name
+
+    insert = "INSERT OR IGNORE INTO system(system_name) VALUES (:system_name);"
+    cur.execute(insert, encoded_data)
+
+
+def iter_system(cur):
+    """Iterate over all system names in the cache.
+
+    Args:
+        cur (:class:`sqlite3.Cursor`):
+            An sqlite3 cursor. This function is meant to be run within a :obj:`with` statement.
+
+    Yields:
+        str: system name
+
+    """
+    select = "SELECT system_name FROM system"
+    for system_name, in cur.execute(select):
+        yield system_name
+
+
+def insert_flux_bias(cur, chain, system, flux_bias, chain_strength, encoded_data=None):
+    """Insert a flux bias offset into the cache.
+
+    Args:
+        cur (:class:`sqlite3.Cursor`):
+            An sqlite3 cursor. This function is meant to be run within a :obj:`with` statement.
+
+        chain (iterable):
+            A collection of nodes. Chains in embedding act as one node.
+
+        system (str):
+            The unique name of a system.
+
+        flux_bias (float):
+            The flux bias offset associated with the given chain.
+
+        chain_strength (float):
+            The magnitude of the negative quadratic bias that induces the given chain in an Ising
+            problem.
+
+        encoded_data (dict, optional):
+            If a dictionary is provided, it will be populated with the serialized data. This is
+            useful for preventing encoding the same information many times.
+
+    """
+    if encoded_data is None:
+        encoded_data = {}
+
+    insert_chain(cur, chain, encoded_data)
+    insert_system(cur, system, encoded_data)
+
+    if 'flux_bias' not in encoded_data:
+        encoded_data['flux_bias'] = _encode_real(flux_bias)
+    if 'chain_strength' not in encoded_data:
+        encoded_data['chain_strength'] = _encode_real(chain_strength)
+    if 'insert_time' not in encoded_data:
+        encoded_data['insert_time'] = datetime.datetime.now()
+
+    insert = \
+        """
+        INSERT OR REPLACE INTO flux_bias(chain_id, system_id, insert_time, flux_bias, chain_strength)
+        SELECT
+            chain.id,
+            system.id,
+            :insert_time,
+            :flux_bias,
+            :chain_strength
+        FROM chain, system
+        WHERE
+            chain.chain_length = :chain_length AND
+            chain.nodes = :nodes AND
+            system.system_name = :system_name;
+        """
+
+    cur.execute(insert, encoded_data)
+
+
+def iter_flux_bias(cur):
+    """Iterate over all flux biases in the cache.
+
+    Args:
+        cur (:class:`sqlite3.Cursor`):
+            An sqlite3 cursor. This function is meant to be run within a :obj:`with` statement.
+
+    Yields:
+        tuple: A 4-tuple:
+
+            list: The chain.
+
+            str: The system name.
+
+            float: The flux bias associated with the chain.
+
+            float: The chain strength associated with the chain.
+
+    """
+    select = \
+        """
+        SELECT nodes, system_name, flux_bias, chain_strength FROM flux_bias_view;
+        """
+
+    for nodes, system, flux_bias, chain_strength in cur.execute(select):
+        yield json.loads(nodes), system, _decode_real(flux_bias), _decode_real(chain_strength)
+
+
+def _encode_real(v):
+    """Encode real numbers as base 64 encoded little endian 8 byte floats."""
+    bytes_ = struct.pack('<d', v)
+    return base64.b64encode(bytes_)
+
+
+def _decode_real(blob):
+    """Inverse of _encode_real."""
+    bytes_ = base64.b64decode(blob)
+    return struct.unpack('<d', bytes_)[0]
+
+
+def get_flux_biases_from_cache(cur, chains, system_name, chain_strength, max_age=3600):
+    """Determine the flux biases for all of the the given chains, system and chain strength.
+
+    Args:
+        cur (:class:`sqlite3.Cursor`):
+            An sqlite3 cursor. This function is meant to be run within a :obj:`with` statement.
+
+        chains (iterable):
+            An iterable of chains. Each chain is a collection of nodes. Chains in embedding act as
+            one node.
+
+        system_name (str):
+            The unique name of a system.
+
+        chain_strength (float):
+            The magnitude of the negative quadratic bias that induces the given chain in an Ising
+            problem.
+
+        max_age (int, optional, default=3600):
+            The maximum age (in seconds) for the flux_bias offsets.
+
+    Returns:
+        dict: A dict where the keys are the nodes in the chains and the values are the flux biases.
+
+    """
+
+    select = \
+        """
+        SELECT
+            flux_bias
+        FROM flux_bias_view WHERE
+            chain_length = :chain_length AND
+            nodes = :nodes AND
+            chain_strength = :chain_strength AND
+            system_name = :system_name AND
+            insert_time >= :time_limit;
+        """
+
+    encoded_data = {'chain_strength': _encode_real(chain_strength),
+                    'system_name': system_name,
+                    'time_limit': datetime.datetime.now() + datetime.timedelta(seconds=-max_age)}
+
+    flux_biases = {}
+    for chain in chains:
+        encoded_data['chain_length'] = len(chain)
+        encoded_data['nodes'] = json.dumps(sorted(chain), separators=(',', ':'))
+
+        row = cur.execute(select, encoded_data).fetchone()
+        if row is None:
+            raise MissingFluxBias
+        flux_bias = _decode_real(*row)
+
+        if flux_bias == 0:
+            continue
+
+        flux_biases.update({v: flux_bias for v in chain})
+
+    return flux_biases
+
+
+def insert_graph(cur, nodelist, edgelist, encoded_data=None):
+    """Insert a graph into the cache.
+
+    A graph is stored by number of nodes, number of edges and a
+    json-encoded list of edges.
+
+    Args:
+        cur (:class:`sqlite3.Cursor`): An sqlite3 cursor. This function
+            is meant to be run within a :obj:`with` statement.
+        nodelist (list): The nodes in the graph.
+        edgelist (list): The edges in the graph.
+        encoded_data (dict, optional): If a dictionary is provided, it
+            will be populated with the serialized data. This is useful for
+            preventing encoding the same information many times.
+
+    Notes:
+        This function assumes that the nodes are index-labeled and range
+        from 0 to num_nodes - 1.
+
+        In order to minimize the total size of the cache, it is a good
+        idea to sort the nodelist and edgelist before inserting.
+
+    Examples:
+        >>> nodelist = [0, 1, 2]
+        >>> edgelist = [(0, 1), (1, 2)]
+        >>> with pmc.cache_connect(':memory:') as cur:
+        ...     pmc.insert_graph(cur, nodelist, edgelist)
+
+        >>> nodelist = [0, 1, 2]
+        >>> edgelist = [(0, 1), (1, 2)]
+        >>> encoded_data = {}
+        >>> with pmc.cache_connect(':memory:') as cur:
+        ...     pmc.insert_graph(cur, nodelist, edgelist, encoded_data)
+        >>> encoded_data['num_nodes']
+        3
+        >>> encoded_data['num_edges']
+        2
+        >>> encoded_data['edges']
+        '[[0,1],[1,2]]'
+
+    """
+    if encoded_data is None:
+        encoded_data = {}
+
+    if 'num_nodes' not in encoded_data:
+        encoded_data['num_nodes'] = len(nodelist)
+    if 'num_edges' not in encoded_data:
+        encoded_data['num_edges'] = len(edgelist)
+    if 'edges' not in encoded_data:
+        encoded_data['edges'] = json.dumps(edgelist, separators=(',', ':'))
+
+    insert = \
+        """
+        INSERT OR IGNORE INTO graph(num_nodes, num_edges, edges)
+        VALUES (:num_nodes, :num_edges, :edges);
+        """
+
+    cur.execute(insert, encoded_data)
+
+
+def iter_graph(cur):
+    """Iterate over all graphs in the cache.
+
+    Args:
+        cur (:class:`sqlite3.Cursor`): An sqlite3 cursor. This function
+            is meant to be run within a :obj:`with` statement.
+
+    Yields:
+        tuple: A 2-tuple containing:
+
+            list: The nodelist for a graph in the cache.
+
+            list: the edgelist for a graph in the cache.
+
+    Examples:
+        >>> nodelist = [0, 1, 2]
+        >>> edgelist = [(0, 1), (1, 2)]
+        >>> with pmc.cache_connect(':memory:') as cur:
+        ...     pmc.insert_graph(cur, nodelist, edgelist)
+        ...     list(pmc.iter_graph(cur))
+        [([0, 1, 2], [[0, 1], [1, 2]])]
+
+    """
+    select = """SELECT num_nodes, num_edges, edges from graph;"""
+    for num_nodes, num_edges, edges in cur.execute(select):
+        yield list(range(num_nodes)), json.loads(edges)
+
+
+def insert_embedding(cur, source_nodelist, source_edgelist, target_nodelist, target_edgelist,
+                     embedding, embedding_tag):
+    """Insert an embedding into the cache.
+
+    Args:
+        cur (:class:`sqlite3.Cursor`):
+            An sqlite3 cursor. This function is meant to be run within a :obj:`with` statement.
+
+        source_nodelist (list):
+            The nodes in the source graph. Should be integer valued.
+
+        source_edgelist (list):
+            The edges in the source graph.
+
+        target_nodelist (list):
+            The nodes in the target graph. Should be integer valued.
+
+        target_edgelist (list):
+            The edges in the target graph.
+
+        embedding (dict):
+            The mapping from the source graph to the target graph.
+            Should be of the form {v: {s, ...}, ...} where v is a variable in the
+            source model and s is a variable in the target model.
+
+        embedding_tag (str):
+            A string tag to associate with the embedding.
+
+    """
+    encoded_data = {}
+
+    # first we need to encode the graphs and create the embedding id
+
+    source_data = {}
+    insert_graph(cur, source_nodelist, source_edgelist, source_data)
+    encoded_data['source_edges'] = source_data['edges']
+    encoded_data['source_num_nodes'] = source_data['num_nodes']
+    encoded_data['source_num_edges'] = source_data['num_edges']
+
+    target_data = {}
+    insert_graph(cur, target_nodelist, target_edgelist, target_data)
+    encoded_data['target_edges'] = target_data['edges']
+    encoded_data['target_num_nodes'] = target_data['num_nodes']
+    encoded_data['target_num_edges'] = target_data['num_edges']
+
+    encoded_data['tag'] = embedding_tag
+
+    insert_embedding = \
+        """
+        INSERT OR REPLACE INTO embedding(
+            source_id,
+            target_id,
+            tag)
+        SELECT
+            source_graph.id,
+            target_graph.id,
+            :tag
+        FROM
+            graph 'source_graph',
+            graph 'target_graph'
+        WHERE
+            source_graph.edges = :source_edges AND
+            source_graph.num_nodes = :source_num_nodes AND
+            source_graph.num_edges = :source_num_edges AND
+            target_graph.edges = :target_edges AND
+            target_graph.num_nodes = :target_num_nodes AND
+            target_graph.num_edges = :target_num_edges
+        """
+
+    cur.execute(insert_embedding, encoded_data)
+
+    # now each chain needs to be inserted
+
+    insert_embedding_component = \
+        """
+        INSERT OR REPLACE INTO embedding_component(
+            source_node,
+            chain_id,
+            embedding_id)
+        SELECT
+            :source_node,
+            chain.id,
+            embedding.id
+        FROM
+            graph 'source_graph',
+            graph 'target_graph',
+            chain,
+            embedding
+        WHERE
+            source_graph.edges = :source_edges AND
+            source_graph.num_nodes = :source_num_nodes AND
+            target_graph.edges = :target_edges AND
+            target_graph.num_nodes = :target_num_nodes AND
+            embedding.source_id = source_graph.id AND
+            embedding.target_id = target_graph.id AND
+            embedding.tag = :tag AND
+            chain.nodes = :nodes AND
+            chain.chain_length = :chain_length
+        """
+
+    for v, chain in embedding.items():
+        chain_data = {'source_node': v}
+        insert_chain(cur, chain, chain_data)
+
+        encoded_data.update(chain_data)
+
+        cur.execute(insert_embedding_component, encoded_data)
+
+
+def select_embedding_from_tag(cur, embedding_tag, target_nodelist, target_edgelist):
+    """Select an embedding from the given tag and target graph.
+
+    Args:
+        cur (:class:`sqlite3.Cursor`):
+            An sqlite3 cursor. This function is meant to be run within a :obj:`with` statement.
+
+        source_nodelist (list):
+            The nodes in the source graph. Should be integer valued.
+
+        source_edgelist (list):
+            The edges in the source graph.
+
+        target_nodelist (list):
+            The nodes in the target graph. Should be integer valued.
+
+        target_edgelist (list):
+            The edges in the target graph.
+
+    Returns:
+        dict: The mapping from the source graph to the target graph.
+        In the form {v: {s, ...}, ...} where v is a variable in the
+        source model and s is a variable in the target model.
+
+    """
+    encoded_data = {'num_nodes': len(target_nodelist),
+                    'num_edges': len(target_edgelist),
+                    'edges': json.dumps(target_edgelist, separators=(',', ':')),
+                    'tag': embedding_tag}
+
+    select = \
+        """
+        SELECT
+            source_node,
+            chain
+        FROM
+            embedding_component_view
+        WHERE
+            embedding_tag = :tag AND
+            target_edges = :edges AND
+            target_num_nodes = :num_nodes AND
+            target_num_edges = :num_edges
+        """
+
+    embedding = {v: json.loads(chain) for v, chain in cur.execute(select, encoded_data)}
+    return embedding
+
+
+def select_embedding_from_source(cur, source_nodelist, source_edgelist,
+                                 target_nodelist, target_edgelist):
+    """Select an embedding from the source graph and target graph.
+
+    Args:
+        cur (:class:`sqlite3.Cursor`):
+            An sqlite3 cursor. This function is meant to be run within a :obj:`with` statement.
+
+        target_nodelist (list):
+            The nodes in the target graph. Should be integer valued.
+
+        target_edgelist (list):
+            The edges in the target graph.
+
+        embedding_tag (str):
+            A string tag to associate with the embedding.
+
+    Returns:
+        dict: The mapping from the source graph to the target graph.
+        In the form {v: {s, ...}, ...} where v is a variable in the
+        source model and s is a variable in the target model.
+
+    """
+    encoded_data = {'target_num_nodes': len(target_nodelist),
+                    'target_num_edges': len(target_edgelist),
+                    'target_edges': json.dumps(target_edgelist, separators=(',', ':')),
+                    'source_num_nodes': len(source_nodelist),
+                    'source_num_edges': len(source_edgelist),
+                    'source_edges': json.dumps(source_edgelist, separators=(',', ':'))}
+
+    select = \
+        """
+        SELECT
+            source_node,
+            chain
+        FROM
+            embedding_component_view
+        WHERE
+            source_num_edges = :source_num_edges AND
+            source_edges = :source_edges AND
+            source_num_nodes = :source_num_nodes AND
+
+            target_num_edges = :target_num_edges AND
+            target_edges = :target_edges AND
+            target_num_nodes = :target_num_nodes
+        """
+
+    embedding = {v: json.loads(chain) for v, chain in cur.execute(select, encoded_data)}
+    return embedding
diff --git a/dwave/system/cache/schema.py b/dwave/system/cache/schema.py
new file mode 100644
index 0000000..b60de41
--- /dev/null
+++ b/dwave/system/cache/schema.py
@@ -0,0 +1,107 @@
+# Copyright 2018 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 schema used by the sqlite database for storing the penalty models."""
+
+schema = \
+    """
+    CREATE TABLE IF NOT EXISTS chain(
+        nodes TEXT UNIQUE,  -- json list of nodes in the chain
+        chain_length INTEGER NOT NULL,  -- number of nodes in the chain
+        id INTEGER PRIMARY KEY);
+
+    CREATE TABLE IF NOT EXISTS system(
+        system_name TEXT UNIQUE,  -- the name of the system
+        id INTEGER PRIMARY KEY);
+
+    CREATE TABLE IF NOT EXISTS flux_bias(
+        chain_id INTEGER NOT NULL,
+        system_id INTEGER NOT NULL,
+        insert_time DATE NOT NULL,  -- flux biases are only valid for relatively short time
+        flux_bias BLOB NOT NULL,
+        chain_strength BLOB NOT NULL,  -- magnitude of the chain strength
+        FOREIGN KEY (chain_id) REFERENCES chain(id) ON DELETE CASCADE,
+        FOREIGN KEY (system_id) REFERENCES system(id) ON DELETE CASCADE,
+        PRIMARY KEY(chain_id, system_id, chain_strength));
+
+    CREATE VIEW IF NOT EXISTS flux_bias_view AS
+        SELECT
+            chain.chain_length,
+            chain.nodes,
+            system.system_name,
+            flux_bias.flux_bias,
+            flux_bias.insert_time,
+            flux_bias.chain_strength
+        FROM flux_bias, chain, system WHERE
+            flux_bias.system_id = system.id AND
+            flux_bias.chain_id = chain.id;
+
+    CREATE TABLE IF NOT EXISTS graph(
+        num_nodes INTEGER NOT NULL,  -- for integer-labeled graphs, num_nodes encodes all of the nodes
+        num_edges INTEGER NOT NULL,  -- redundant, allows for faster selects
+        edges TEXT NOT NULL,  -- json list of lists, should be sorted (with each edge sorted)
+        id INTEGER PRIMARY KEY,
+        CONSTRAINT graph UNIQUE (
+            num_nodes,
+            edges));
+
+    CREATE TABLE IF NOT EXISTS embedding(
+        source_id INTEGER NOT NULL,
+        target_id INTEGER NOT NULL,
+        tag TEXT,  -- user specified
+        id INTEGER PRIMARY KEY,
+        protected INTEGER DEFAULT 0,
+        CONSTRAINT embedding_constraint UNIQUE (target_id, tag),
+        CONSTRAINT embedding_constraint UNIQUE (target_id, source_id),
+        FOREIGN KEY (source_id) REFERENCES graph(id) ON DELETE CASCADE,
+        FOREIGN KEY (target_id) REFERENCES graph(id) ON DELETE CASCADE);
+
+    CREATE TABLE IF NOT EXISTS embedding_component(
+        source_node INTEGER NOT NULL,
+        chain_id INTEGER NOT NULL,
+        embedding_id INTEGER NOT NULL,
+        FOREIGN KEY (chain_id) REFERENCES chain(id) ON DELETE CASCADE,
+        FOREIGN KEY (embedding_id) REFERENCES embedding(id) ON DELETE CASCADE,
+        PRIMARY KEY(embedding_id, source_node));
+
+    CREATE VIEW IF NOT EXISTS embedding_component_view AS
+        SELECT
+            source_graph.num_nodes 'source_num_nodes',
+            source_graph.num_edges 'source_num_edges',
+            source_graph.edges 'source_edges',
+
+            target_graph.num_nodes 'target_num_nodes',
+            target_graph.num_edges 'target_num_edges',
+            target_graph.edges 'target_edges',
+
+            embedding_component.source_node,
+
+            chain.chain_length,
+            chain.nodes 'chain',
+
+            embedding.tag 'embedding_tag'
+        FROM
+            graph 'source_graph',
+            graph 'target_graph',
+            embedding,
+            embedding_component,
+            chain
+
+        WHERE
+            embedding.source_id = source_graph.id AND
+            embedding.target_id = target_graph.id AND
+            embedding_component.chain_id = chain.id AND
+            embedding.id = embedding_component.embedding_id;
+    """
diff --git a/dwave/system/composites/__init__.py b/dwave/system/composites/__init__.py
index b4e1902..6a847af 100644
--- a/dwave/system/composites/__init__.py
+++ b/dwave/system/composites/__init__.py
@@ -1,4 +1,4 @@
-# Copyright 2018 D-Wave Systems Inc.
+# Copyright 2024 D-Wave Inc.
 #
 #    Licensed under the Apache License, Version 2.0 (the "License");
 #    you may not use this file except in compliance with the License.
@@ -14,6 +14,7 @@
 
 from dwave.system.composites.cutoffcomposite import *
 from dwave.system.composites.embedding import *
+from dwave.system.composites.linear_ancilla import *
 from dwave.system.composites.tiling import *
 from dwave.system.composites.virtual_graph import *
 from dwave.system.composites.reversecomposite import *
diff --git a/dwave/system/composites/cutoffcomposite.py b/dwave/system/composites/cutoffcomposite.py
index 201aebf..54d7c51 100644
--- a/dwave/system/composites/cutoffcomposite.py
+++ b/dwave/system/composites/cutoffcomposite.py
@@ -67,9 +67,8 @@ class CutOffComposite(dimod.ComposedSampler):
         ...                            {'ab': 0.8, 'ac': 0.7, 'bc': -1},
         ...                            0,
         ...                            dimod.SPIN)
-        >>> samples = CutOffComposite(
-        ...     AutoEmbeddingComposite(sampler), 0.75).sample(bqm, num_reads=1000)
-        >>> print(samples.first.energy)
+        >>> CutOffComposite(AutoEmbeddingComposite(sampler), 0.75).sample(bqm,
+        ...                 num_reads=1000).first.energy
         -5.5
 
     """
@@ -239,8 +238,7 @@ class PolyCutOffComposite(dimod.ComposedPolySampler):
         >>> import dimod
         >>> sampler = dimod.HigherOrderComposite(dimod.ExactSolver())
         >>> poly = dimod.BinaryPolynomial({'a': 3, 'abc':-4, 'ac': 0.2}, dimod.SPIN)
-        >>> samples = PolyCutOffComposite(sampler, 1).sample_poly(poly)
-        >>> print(samples.first.sample['a'])
+        >>> PolyCutOffComposite(sampler, 1).sample_poly(poly).first.sample['a']
         -1
 
     """
diff --git a/dwave/system/composites/embedding.py b/dwave/system/composites/embedding.py
index 619e669..e29f95d 100644
--- a/dwave/system/composites/embedding.py
+++ b/dwave/system/composites/embedding.py
@@ -30,6 +30,7 @@ from warnings import warn
 
 import dimod
 import minorminer
+import functools
 
 from dwave.embedding import (target_to_source, unembed_sampleset, embed_bqm,
                              chain_to_quadratic, EmbeddedStructure)
@@ -81,7 +82,7 @@ class EmbeddingComposite(dimod.ComposedSampler):
        >>> h = {'a': -1., 'b': 2}
        >>> J = {('a', 'b'): 1.5}
        >>> sampleset = sampler.sample_ising(h, J, num_reads=100)
-       >>> print(sampleset.first.energy)
+       >>> sampleset.first.energy
        -4.5
 
 
diff --git a/dwave/system/composites/linear_ancilla.py b/dwave/system/composites/linear_ancilla.py
new file mode 100644
index 0000000..2f00342
--- /dev/null
+++ b/dwave/system/composites/linear_ancilla.py
@@ -0,0 +1,232 @@
+# Copyright 2024 D-Wave 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.
+
+"""Composite to implement linear coefficients with ancilla qubits biased with flux-bias offsets.
+"""
+
+import numbers
+
+from collections import defaultdict
+from typing import Sequence, Mapping, Any
+
+import dimod
+import numpy as np
+
+from dimod.decorators import nonblocking_sample_method
+
+
+__all__ = ["LinearAncillaComposite"]
+
+
+class LinearAncillaComposite(dimod.ComposedSampler, dimod.Structured):
+    """Implements linear biases as ancilla qubits polarized with strong flux biases.
+
+    Linear bias :math:`h_i` of qubit :math:`i` is implemented through a coupling 
+    :math:`J_{ij}` between the qubit and a neighboring qubit :math:`j` that has a 
+    large flux-bias offset.
+
+    Args:
+    child_sampler (:class:`dimod.Sampler`):
+        A dimod sampler, such as a :class:`~dwave.system.samplers.DWaveSampler()`, 
+        that has flux bias controls.
+
+    Examples:
+        This example submits a two-qubit problem consisting of linear biases with opposed signs 
+        and anti-ferromagnetic coupling. A D-Wave quantum computer solves it with the fast-anneal 
+        protocol using ancilla qubits to represent the linear biases.
+
+        >>> from dwave.system import DWaveSampler, EmbeddingComposite, LinearAncillaComposite
+        ...
+        >>> sampler = EmbeddingComposite(LinearAncillaComposite(DWaveSampler()))
+        >>> sampleset = sampler.sample_ising({0:1, 1:-1}, {(0, 1): 1}, fast_anneal=True) # doctest: +SKIP
+        >>> sampleset.first.energy # doctest: +SKIP
+        -3
+
+    """
+
+    def __init__(
+        self,
+        child_sampler: dimod.Sampler,
+    ):
+        self.children = [child_sampler]
+        self.parameters = child_sampler.parameters.copy()
+        self.properties = dict(child_properties=child_sampler.properties.copy())
+        self.nodelist = child_sampler.nodelist
+        self.edgelist = child_sampler.edgelist
+
+    def nodelist(self):
+        pass  # overwritten by init
+
+    def edgelist(self):
+        pass  # overwritten by init
+
+    children = None  # overwritten by init
+    """list [child_sampler]: List containing the structured sampler."""
+
+    parameters = None  # overwritten by init
+    """dict[str, list]: Parameters in the form of a dict.
+
+    For an instantiated composed sampler, keys are the keyword parameters
+    accepted by the child sampler and parameters added by the composite.
+    """
+
+    properties = None  # overwritten by init
+    """dict: Properties in the form of a dict.
+
+    Contains the properties of the child sampler.
+    """
+
+    @nonblocking_sample_method
+    def sample(
+        self,
+        bqm: dimod.BinaryQuadraticModel,
+        *,
+        h_tolerance: numbers.Number = 0,
+        default_flux_bias_range: tuple[float, float] = (-0.005, 0.005),
+        **parameters,
+    ):
+        """Sample from the provided binary quadratic model.
+
+        .. note::
+            This composite does not support the :ref:`param_autoscale` parameter; use the
+            :class:`~dwave.preprocessing.composites.ScaleComposite` for scaling.
+
+        Args:
+            bqm (:class:`~dimod.binary.BinaryQuadraticModel`):
+                Binary quadratic model to be sampled from.
+
+            h_tolerance (:class:`numbers.Number`):
+                Magnitude of the linear bias to be set directly on problem qubits; above this the bias
+                is emulated by the flux-bias offset to an ancilla qubit. Assumed to be positive. 
+                Defaults to zero.
+
+            default_flux_bias_range (:class:`tuple`):
+                Flux-bias range, as a two-tuple, supported by the QPU. The values must be large enough to
+                ensure qubits remain polarized throughout the annealing process.
+
+            **parameters:
+                Parameters for the sampling method, specified by the child
+                sampler.
+
+        Returns:
+            :class:`~dimod.SampleSet`.
+
+        """
+        if h_tolerance < 0:
+            raise ValueError("h_tolerance needs to be positive or zero")
+
+        child = self.child
+        qpu_properties = _innermost_child_properties(child)
+        target_graph = child.to_networkx_graph()
+        source_graph = dimod.to_networkx_graph(bqm)
+        extended_j_range = qpu_properties["extended_j_range"]
+        # flux_bias_range is not supported at the moment
+        flux_bias_range = qpu_properties.get("flux_bias_range", default_flux_bias_range)
+
+        # Positive couplings tend to have smaller control error,
+        # we default to them if they have the same magnitude than negative couplings
+        # https://docs.dwavesys.com/docs/latest/c_qpu_ice.html#overview-of-ice
+        largest_j = max(extended_j_range[::-1], key=abs)
+        largest_j_sign = np.sign(largest_j)
+
+        # To implement the bias sign through flux bias sign,
+        # we pick a range (magnitude) that we can sign-flip
+        fb_magnitude = min(abs(b) for b in flux_bias_range)
+        flux_biases = [0] * qpu_properties["num_qubits"]
+
+        _bqm = bqm.copy()
+        used_ancillas = defaultdict(list)
+        for variable, bias in bqm.iter_linear():
+            if abs(bias) <= h_tolerance:
+                continue
+            if abs(bias) - h_tolerance > abs(largest_j):
+                return ValueError(
+                    "linear biases larger than the strongest coupling are not supported"
+                )  # TODO: implement larger biases through multiple ancillas
+
+            available_ancillas = set(target_graph.adj[variable]) - source_graph.nodes()
+            if not len(available_ancillas):
+                raise ValueError(f"variable {variable} has no ancillas available")
+            unused_ancillas = available_ancillas - used_ancillas.keys()
+            if len(unused_ancillas):
+                ancilla = unused_ancillas.pop()
+                # bias sign is handled by the flux bias
+                flux_biases[ancilla] = np.sign(bias) * largest_j_sign * fb_magnitude
+                _bqm.add_interaction(
+                    variable, ancilla, (abs(bias) - h_tolerance) * largest_j_sign
+                )
+            else:
+                if qpu_properties["j_range"][0] <= bias <= qpu_properties["j_range"][1]:
+                    # If j can be sign-flipped, select the least used ancilla regardless of the flux bias sign
+                    ancilla = sorted(
+                        list(available_ancillas), key=lambda x: len(used_ancillas[x])
+                    )[0]
+                    _bqm.add_interaction(
+                        variable,
+                        ancilla,
+                        (bias - h_tolerance * np.sign(bias))
+                        * np.sign([flux_biases[ancilla]]),
+                    )
+                else:
+                    # Ancilla sharing is limited to flux biases with appropiate sign
+                    signed_ancillas = [
+                        ancilla
+                        for ancilla in available_ancillas
+                        if largest_j_sign
+                        == np.sign(flux_biases[ancilla] * bias)
+                    ]
+                    if not len(signed_ancillas):
+                        return ValueError(
+                            f"variable {variable} has no ancillas available"
+                        )
+                    else:
+                        ancilla = sorted(
+                            list(signed_ancillas), key=lambda x: len(used_ancillas[x])
+                        )[0]
+                        _bqm.add_interaction(
+                            variable,
+                            ancilla,
+                            largest_j_sign * (abs(bias) - h_tolerance),
+                        )
+
+            used_ancillas[ancilla].append(variable)
+            _bqm.set_linear(variable, h_tolerance * np.sign(bias))
+
+        sampleset = self.child.sample(_bqm, flux_biases=flux_biases, **parameters)
+        yield
+        yield dimod.SampleSet.from_samples_bqm(
+            [
+                {k: v for k, v in sample.items() if k not in used_ancillas}
+                for sample in sampleset.samples()
+            ],
+            bqm=bqm,
+            info=sampleset.info.update(used_ancillas),
+        )
+
+
+def _innermost_child_properties(sampler: dimod.Sampler) -> Mapping[str, Any]:
+    """Returns the properties of the inner-most child sampler in a composite.
+
+    Args:
+        sampler: A dimod sampler
+
+    Returns:
+        properties (dict): The properties of the inner-most sampler
+
+    """
+
+    try:
+        return _innermost_child_properties(sampler.child)
+    except AttributeError:
+        return sampler.properties
diff --git a/dwave/system/composites/tiling.py b/dwave/system/composites/tiling.py
index 8808a63..0752ef3 100644
--- a/dwave/system/composites/tiling.py
+++ b/dwave/system/composites/tiling.py
@@ -25,10 +25,14 @@ fully-yielded Advantage system.
 
 See `Ocean Glossary <https://docs.ocean.dwavesys.com/en/stable/concepts/index.html>`_
 for explanations of technical terms in descriptions of Ocean tools.
+
 """
 
+from math import sqrt, ceil
+
 import dimod
 import dwave_networkx as dnx
+import numpy as np
 
 import dwave.embedding
 
diff --git a/dwave/system/composites/virtual_graph.py b/dwave/system/composites/virtual_graph.py
index 622ef03..91991a0 100644
--- a/dwave/system/composites/virtual_graph.py
+++ b/dwave/system/composites/virtual_graph.py
@@ -12,79 +12,202 @@
 #    See the License for the specific language governing permissions and
 #    limitations under the License.
 
-"""Deprecated and functionality removed. Behaves identically to
-:class:`~dwave.system.composites.embedding.FixedEmbeddingComposite`, and will
-be completely removed in dwave-system 2.0.
-
-Virtual graphs are deprecated due to improved calibration of newer QPUs; to
-calibrate chains for residual biases, follow the instructions in the
-`shimming tutorial <https://github.com/dwavesystems/shimming-tutorial>`_.
-
-A :std:doc:`dimod composite <oceandocs:docs_dimod/reference/samplers>` that
-uses the D-Wave virtual graph feature for improved
+"""
+A :std:doc:`dimod composite <oceandocs:docs_dimod/reference/samplers>` that 
+uses the D-Wave virtual graph feature for improved 
 :std:doc:`minor-embedding <oceandocs:docs_system/intro>`.
 
-D-Wave *virtual graphs* simplify the process of minor-embedding by enabling you
-to more easily create, optimize, use, and reuse an embedding for a given working
-graph. When you submit an embedding and specify a chain strength using these
-tools, they automatically calibrate the qubits in a chain to compensate for the
+D-Wave *virtual graphs* simplify the process of minor-embedding by enabling you 
+to more easily create, optimize, use, and reuse an embedding for a given working 
+graph. When you submit an embedding and specify a chain strength using these 
+tools, they automatically calibrate the qubits in a chain to compensate for the 
 effects of biases that may be introduced as a result of strong couplings.
 
 See `Ocean Glossary <https://docs.ocean.dwavesys.com/en/stable/concepts/index.html>`_
 for explanations of technical terms in descriptions of Ocean tools.
 """
 
-import warnings
-
 import dimod
 
 from dwave.system.composites.embedding import FixedEmbeddingComposite
+from dwave.system.flux_bias_offsets import get_flux_biases
+
+
+FLUX_BIAS_KWARG = 'flux_biases'
 
 __all__ = ['VirtualGraphComposite']
 
 
 class VirtualGraphComposite(FixedEmbeddingComposite):
-    """Removed. Used to provide access to the D-Wave virtual graph feature for
-    minor-embedding, but now is just a thin wrapper around the
-    :class:`~dwave.system.composites.embedding.FixedEmbeddingComposite`.
+    """Composite to use the D-Wave virtual graph feature for minor-embedding.
+
+    Calibrates qubits in chains to compensate for the effects of biases and enables easy
+    creation, optimization, use, and reuse of an embedding for a given working graph.
+
+    Inherits from :class:`~dimod.ComposedSampler` and :class:`~dimod.Structured`.
+
+    Args:
+        sampler (:class:`DWaveSampler`):
+            A dimod :class:`~dimod.Sampler`. Typically a :obj:`DWaveSampler` or
+            derived composite sampler; other samplers may not work or make sense 
+            with this composite layer.
+
+        embedding (dict[hashable, iterable]):
+            Mapping from a source graph to the specified sampler's graph (the 
+            target graph).
+
+        chain_strength (float, optional, default=None):
+            Desired chain coupling strength. This is the magnitude of couplings 
+            between qubits in a chain. If None, uses the maximum available as 
+            returned by a SAPI query to the D-Wave solver.
+
+        flux_biases (list/False/None, optional, default=None):
+            Per-qubit flux bias offsets in the form of a list of lists, where 
+            each sublist is of length 2 and specifies a variable and the flux 
+            bias offset associated with that variable. Qubits in a chain with 
+            strong negative J values experience a J-induced bias; this parameter 
+            compensates by recalibrating to remove that bias. If False, no flux 
+            bias is applied or calculated. If None, flux biases are pulled from 
+            the database or calculated empirically.
+
+        flux_bias_num_reads (int, optional, default=1000):
+            Number of samples to collect per flux bias value to calculate 
+            calibration information.
+
+        flux_bias_max_age (int, optional, default=3600):
+            Maximum age (in seconds) allowed for a previously calculated flux 
+            bias offset to be considered valid.
+
+    .. attention::
+        D-Wave's *virtual graphs* feature can require many seconds of D-Wave 
+        quantum computer time to calibrate qubits to compensate for the effects 
+        of biases. If your account has limited D-Wave system access, consider 
+        using :class:`FixedEmbeddingComposite` instead.
+
+    Examples:
+        This example uses :class:`.VirtualGraphComposite` to instantiate a 
+        composed sampler that submits an Ising problem to a D-Wave solver. 
+        This simple three-variable problem is manually minor-embedded such that
+        variables ``a`` and ``b`` are represented by single qubits while variable 
+        ``c`` is represented by a four-qubit chain. The chain strength is set to 
+        the maximum allowed found from querying the solver's extended J range. 
+        The minor embedding shown below was for an execution of this example on a 
+        particular Advantage system; select a suitable embedding for the QPU you 
+        use.
+
+        >>> from dwave.system import DWaveSampler, VirtualGraphComposite
+        ...
+        >>> h = {'a': 1, 'b': -1}
+        >>> J = {('b', 'c'): -1, ('a', 'c'): -1}
+        ...
+        >>> qpu = DWaveSampler()
+        >>> embedding = {'a': [2656], 'c': [2641, 2642, 2643, 2644], 'b': [2659]}
+        >>> qpu.properties['extended_j_range']
+        [-2.0, 1.0]
+        >>> # Sample using VirtualGraphComposite
+        >>> sampler = VirtualGraphComposite(qpu, embedding, chain_strength=2) # doctest: +SKIP
+        >>> sampleset = sampler.sample_ising(h, J, num_reads=100) # doctest: +SKIP
+        >>> print(sampleset)    # doctest: +SKIP
+            a  b  c energy num_oc. chain_.
+        0 +1 +1 +1   -2.0      21     0.0
+        1 -1 +1 +1   -2.0      66     0.0
+        2 -1 -1 -1   -2.0       8     0.0
+        3 -1 +1 -1   -2.0       5     0.0
+        ['SPIN', 4 rows, 100 samples, 3 variables]
+
+    See `Ocean Glossary <https://docs.ocean.dwavesys.com/en/stable/concepts/index.html>`_
+    for explanations of technical terms in descriptions of Ocean tools.
 
-    .. deprecated:: 1.25.0
-        This class is deprecated due to improved calibration of newer QPUs and
-        will be removed in 1.27.0; to calibrate chains for residual biases, 
-        follow the instructions in the 
-        `shimming tutorial <https://github.com/dwavesystems/shimming-tutorial>`_.
+    """
 
-    .. versionremoved:: 1.28.0
-        This class is now only a pass-through wrapper around the
-        :class:`~dwave.system.composites.embedding.FixedEmbeddingComposite`.
+    def __init__(self, sampler, embedding,
+                 chain_strength=None,
+                 flux_biases=None,
+                 flux_bias_num_reads=1000,
+                 flux_bias_max_age=3600):
 
-        It will be completely removed in dwave-system 2.0.
+        super(VirtualGraphComposite, self).__init__(sampler, embedding)
+        self.parameters.update(apply_flux_bias_offsets=[])
 
-        For removal reasons and alternatives, see the deprecation note above.
-    """
+        # Validate the chain strength, or obtain it from J-range if chain strength is not provided.
+        self.chain_strength = _validate_chain_strength(sampler, chain_strength)
 
-    def __init__(self, sampler, embedding, chain_strength=None, **kwargs):
-        super(VirtualGraphComposite, self).__init__(sampler, embedding)
+        if flux_biases is False:  # use 'is' because bool(None) is False
+            # in this case we are done
+            self.flux_biases = None
+            return
 
-        warnings.warn(
-            "'VirtualGraphComposite' functionality is removed due to improved "
-            "calibration of newer QPUs and in future will raise an exception. "
-            "Currently it's equivalent to 'FixedEmbeddingComposite'. If needed, "
-            "follow the instructions in the shimming tutorial at "
-            "https://github.com/dwavesystems/shimming-tutorial instead.",
-            DeprecationWarning, stacklevel=2
-        )
-
-        if chain_strength is not None:
-            warnings.warn(
-                "'chain_strength' parameter is ignored since dwave-system 1.28.",
-                DeprecationWarning, stacklevel=2)
-
-        # for API backwards compatibility
-        self.parameters.update(apply_flux_bias_offsets=[])
-        self.chain_strength = chain_strength
-        self.flux_biases = None
+        if FLUX_BIAS_KWARG not in sampler.parameters:
+            raise ValueError("Given child sampler does not accept flux_biases.")
+
+        # come back as a dict
+        flux_biases = get_flux_biases(sampler, embedding,
+                                      num_reads=flux_bias_num_reads,
+                                      chain_strength=self.chain_strength,
+                                      max_age=flux_bias_max_age)
+
+        self.flux_biases = [flux_biases.get(v, 0.0) for v in range(sampler.properties['num_qubits'])]
+
+        return
 
     @dimod.bqm_structured
     def sample(self, bqm, apply_flux_bias_offsets=True, **kwargs):
+        """Sample from the given Ising model.
+
+        Args:
+
+            h (list/dict):
+                Linear biases of the Ising model. If a list, the list's indices
+                are used as variable labels.
+
+            J (dict of (int, int):float):
+                Quadratic biases of the Ising model.
+
+            apply_flux_bias_offsets (bool, optional):
+                If True, use the calculated flux_bias offsets (if available).
+
+            **kwargs:
+                Optional keyword arguments for the sampling method, specified per solver.
+
+        See `Ocean Glossary <https://docs.ocean.dwavesys.com/en/stable/concepts/index.html>`_
+        for explanations of technical terms in descriptions of Ocean tools.
+
+        """
+
+        if apply_flux_bias_offsets:
+            if self.flux_biases is not None:
+                kwargs[FLUX_BIAS_KWARG] = self.flux_biases
+
+        kwargs.setdefault('chain_strength', self.chain_strength)
+
         return super(VirtualGraphComposite, self).sample(bqm, **kwargs)
+
+
+def _validate_chain_strength(sampler, chain_strength):
+    """Validate the provided chain strength, checking J-ranges of the sampler's 
+    children.
+
+    Args:
+        chain_strength (float): The provided chain strength.  Use None to use 
+            J-range.
+
+    Returns (float):
+        A valid chain strength, either provided or based on available J-range.  
+        Positive finite float.
+
+    """
+    properties = sampler.properties
+
+    if 'extended_j_range' in properties:
+        max_chain_strength = - min(properties['extended_j_range'])
+    elif 'j_range' in properties:
+        max_chain_strength = - min(properties['j_range'])
+    else:
+        raise ValueError("input sampler should have 'j_range' and/or 'extended_j_range' property.")
+
+    if chain_strength is None:
+        chain_strength = max_chain_strength
+    elif chain_strength > max_chain_strength:
+        raise ValueError("Provided chain strength exceedds the allowed range.")
+
+    return chain_strength
diff --git a/dwave/system/exceptions.py b/dwave/system/exceptions.py
index ca04985..770d41f 100644
--- a/dwave/system/exceptions.py
+++ b/dwave/system/exceptions.py
@@ -13,6 +13,10 @@
 #    limitations under the License.
 
 
+class MissingFluxBias(Exception):
+    """"""
+
+
 class FailoverCondition(Exception):
     """QPU or Solver API call failed with an error that might be mitigated by
     retrying on a different solver.
diff --git a/dwave/system/flux_bias_offsets.py b/dwave/system/flux_bias_offsets.py
index b0f9d1e..af87d77 100644
--- a/dwave/system/flux_bias_offsets.py
+++ b/dwave/system/flux_bias_offsets.py
@@ -15,34 +15,74 @@
 
 import warnings
 
+import dimod
 
-def get_flux_biases(*args, **kwargs):
-    """Removed. Used to get the flux bias offsets for sampler and embedding,
-    but now just returns an empty dictionary.
+from dwave.system.cache.database_manager import get_flux_biases_from_cache, cache_connect, insert_flux_bias
+from dwave.system.exceptions import MissingFluxBias
+from dwave.system.samplers.dwave_sampler import DWaveSampler
 
-    .. versionremoved:: 1.28.0
-        Due to improved calibration of newer QPUs, balancing chains with
-        the flux biases, as it used to be implemented by this function (and the
-        obsoleted package ``dwave-drivers``) is no longer supported.
 
-        To calibrate chains for residual biases, follow the instructions in the
-        `shimming tutorial <https://github.com/dwavesystems/shimming-tutorial>`_.
+def get_flux_biases(sampler, embedding, chain_strength, num_reads=1000, max_age=3600):
+    """Get the flux bias offsets for sampler and embedding.
 
-        This function (and its submodule) will be completely removed in
-        dwave-system 2.0.
+    Args:
+        sampler (:obj:`.DWaveSampler`):
+            A D-Wave sampler.
+
+        embedding (dict[hashable, iterable]):
+            Mapping from a source graph to the specified sampler’s graph (the target graph). The
+            keys of embedding should be nodes in the source graph, the values should be an iterable
+            of nodes in the target graph.
+
+        chain_strength (number):
+            Desired chain coupling strength. This is the magnitude of couplings between qubits
+            in a chain.
+
+        num_reads (int, optional, default=1000):
+            The number of reads per system call if new flux biases need to be calculated.
+
+        max_age (int, optional, default=3600):
+            The maximum age (in seconds) allowed for previously calculated flux bias offsets.
 
     Returns:
-        dict:
-            An empty dict, since dwave-system 1.28.
-            Flux biases are not calculated/set for nodes/chains anymore.
+        dict: A dict where the keys are the nodes in the chains and the values are the flux biases.
+
     """
 
-    warnings.warn(
-        "'get_flux_biases' functionality is removed due to improved calibration "
-        "of newer QPUs and in future will raise an exception; if needed, "
-        "follow the instructions in the shimming tutorial at "
-        "https://github.com/dwavesystems/shimming-tutorial instead. ",
-        DeprecationWarning, stacklevel=2
-    )
+    if not isinstance(sampler, dimod.Sampler):
+        raise TypeError("input sampler should be DWaveSampler")
+
+    # try to read the chip_id, otherwise get the name
+    system_name = sampler.properties.get('chip_id', str(sampler.__class__))
+
+    try:
+        with cache_connect() as cur:
+            fbo = get_flux_biases_from_cache(cur, embedding.values(), system_name,
+                                             chain_strength=chain_strength,
+                                             max_age=max_age)
+        return fbo
+    except MissingFluxBias:
+        pass
+
+    # if dwave-drivers is not available, then we can't calculate the biases
+    try:
+        import dwave.drivers as drivers
+    except ImportError:
+        msg = ("dwave-drivers not found, cannot calculate flux biases. dwave-drivers can be "
+               "installed with "
+               "'pip install dwave-drivers --extra-index-url https://pypi.dwavesys.com/simple'. "
+               "See documentation for dwave-drivers license.")
+        raise RuntimeError(msg)
+
+    fbo = drivers.oneshot_flux_bias(sampler, embedding.values(), num_reads=num_reads,
+                                    chain_strength=chain_strength,
+                                    label='VirtualGraph flux bias measurements')
+
+    # store them in the cache
+    with cache_connect() as cur:
+        for chain in embedding.values():
+            v = next(iter(chain))
+            flux_bias = fbo.get(v, 0.0)
+            insert_flux_bias(cur, chain, system_name, flux_bias, chain_strength)
 
-    return {}
+    return fbo
diff --git a/dwave/system/package_info.py b/dwave/system/package_info.py
index f83e474..87c9752 100644
--- a/dwave/system/package_info.py
+++ b/dwave/system/package_info.py
@@ -14,7 +14,22 @@
 
 __all__ = ['__version__', '__author__', '__authoremail__', '__description__']
 
-__version__ = '1.27.0'
+__version__ = '1.25.0'
 __author__ = 'D-Wave Systems Inc.'
 __authoremail__ = 'tools@dwavesys.com'
 __description__ = 'All things D-Wave System.'
+
+
+# register the non-open-source packages required for dwave-drivers to work
+contrib = [{
+    'name': 'drivers',
+    'title': 'D-Wave Drivers',
+    'description': 'These drivers enable some automated performance-tuning features.',
+    'license': {
+        'name': 'D-Wave EULA',
+        'url': 'https://docs.ocean.dwavesys.com/eula',
+    },
+    'requirements': [
+        'dwave-drivers>=0.4.4',
+    ]
+}]
diff --git a/dwave/system/samplers/clique.py b/dwave/system/samplers/clique.py
index b39443c..ed0fd7a 100644
--- a/dwave/system/samplers/clique.py
+++ b/dwave/system/samplers/clique.py
@@ -17,6 +17,7 @@ from typing import Tuple
 
 import dimod
 import networkx as nx
+import dwave_networkx as dnx
 
 from minorminer.busclique import find_clique_embedding, busgraph_cache
 from dwave.preprocessing import ScaleComposite
diff --git a/dwave/system/samplers/dwave_sampler.py b/dwave/system/samplers/dwave_sampler.py
index bbcf64a..48886d7 100644
--- a/dwave/system/samplers/dwave_sampler.py
+++ b/dwave/system/samplers/dwave_sampler.py
@@ -26,7 +26,7 @@ from collections import defaultdict
 from typing import Optional, Dict
 
 import dimod
-import dwave_networkx as dnx
+
 from dimod.exceptions import BinaryQuadraticModelStructureError
 from dwave.cloud.client import Client
 from dwave.cloud.exceptions import (
@@ -38,6 +38,8 @@ from dwave.cloud.exceptions import (
 from dwave.system.exceptions import FailoverCondition, RetryCondition
 from dwave.system.warnings import WarningHandler, WarningAction
 
+import dwave_networkx as dnx
+
 __all__ = ['DWaveSampler', 'qpu_graph']
 
 
@@ -155,7 +157,7 @@ class DWaveSampler(dimod.Sampler, dimod.Structured):
         >>> sampleset = sampler.sample_ising({qubit_a: -1, qubit_b: 1},
         ...                                  {},
         ...                                  num_reads=100)
-        >>> print(sampleset.first.sample[qubit_a] == 1 and sampleset.first.sample[qubit_b] == -1)
+        >>> sampleset.first.sample[qubit_a] == 1 and sampleset.first.sample[qubit_b] == -1
         True
 
     See `Ocean Glossary <https://docs.ocean.dwavesys.com/en/stable/concepts/index.html>`_
@@ -396,7 +398,7 @@ class DWaveSampler(dimod.Sampler, dimod.Structured):
             >>> sampleset = sampler.sample_ising({qubit_a: -1, qubit_b: 1},
             ...                                  {},
             ...                                  num_reads=100)
-            >>> print(sampleset.first.sample[qubit_a] == 1 and sampleset.first.sample[qubit_b] == -1)
+            >>> sampleset.first.sample[qubit_a] == 1 and sampleset.first.sample[qubit_b] == -1
             True
 
         See `Ocean Glossary <https://docs.ocean.dwavesys.com/en/stable/concepts/index.html>`_
@@ -591,3 +593,4 @@ class DWaveSampler(dimod.Sampler, dimod.Structured):
         return qpu_graph(self.properties['topology']['type'],
                          self.properties['topology']['shape'],
                          self.nodelist, self.edgelist)
+        
diff --git a/dwave/system/samplers/leap_hybrid_sampler.py b/dwave/system/samplers/leap_hybrid_sampler.py
index 96c03eb..629432a 100644
--- a/dwave/system/samplers/leap_hybrid_sampler.py
+++ b/dwave/system/samplers/leap_hybrid_sampler.py
@@ -25,7 +25,7 @@ from typing import Any, Dict, List, NamedTuple, Optional
 import dimod
 import dwave.optimization
 import numpy
-from dwave.cloud.client import Client
+from dwave.cloud import Client
 
 from dwave.system.utilities import classproperty, FeatureFlags
 
diff --git a/dwave/system/temperatures.py b/dwave/system/temperatures.py
index 30cdbd4..ebe3c54 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.
 
-r"""The following effective temperature and bias 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) 
@@ -48,13 +48,13 @@ __all__ = ['effective_field', 'maximum_pseudolikelihood_temperature',
 def effective_field(bqm,
                     samples=None,
                     current_state_energy=False) -> (np.ndarray,list):
-    r'''Returns the effective field for all variables and all samples.
-
+    '''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 
@@ -64,12 +64,11 @@ def effective_field(bqm,
     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
+    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::
@@ -77,7 +76,7 @@ def effective_field(bqm,
     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.
@@ -104,7 +103,7 @@ def effective_field(bqm,
        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
@@ -115,7 +114,7 @@ def effective_field(bqm,
        >>> 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))
@@ -126,7 +125,7 @@ def effective_field(bqm,
     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
@@ -138,7 +137,7 @@ def effective_field(bqm,
     if current_state_energy is True:
         #Ising: eff_field = 2*s*(h + J*s)
         effective_fields = 2*samples*effective_fields
-
+        
     return (effective_fields,labels)
 
 def maximum_pseudolikelihood_temperature(bqm=None,
@@ -149,8 +148,8 @@ def maximum_pseudolikelihood_temperature(bqm=None,
                                          T_guess=None,
                                          optimize_method='bisect',
                                          T_bracket=(1e-3,1000)) -> Tuple[float,np.ndarray]:
-    r'''Returns a sampling-based temperature estimate.
-
+    '''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. 
@@ -162,15 +161,15 @@ def maximum_pseudolikelihood_temperature(bqm=None,
     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.
 
@@ -180,7 +179,7 @@ def maximum_pseudolikelihood_temperature(bqm=None,
     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.
@@ -231,7 +230,7 @@ def maximum_pseudolikelihood_temperature(bqm=None,
        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
@@ -242,15 +241,15 @@ def maximum_pseudolikelihood_temperature(bqm=None,
        >>> 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)
 
@@ -340,7 +339,7 @@ def maximum_pseudolikelihood_temperature(bqm=None,
                     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
@@ -359,7 +358,7 @@ def maximum_pseudolikelihood_temperature(bqm=None,
                     site_energy = (site_energy[0][indices,:],site_energy[1]),
                     num_bootstrap_samples = 0,
                     T_guess = T_estimate)
-
+    
     return T_estimate, T_bootstrap_estimates
 
 def Ip_in_units_of_B(Ip: Union[None, float, np.ndarray]=None,
@@ -368,7 +367,7 @@ def Ip_in_units_of_B(Ip: Union[None, float, np.ndarray]=None,
                      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.
+    """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 
@@ -377,7 +376,7 @@ def Ip_in_units_of_B(Ip: Union[None, float, np.ndarray]=None,
     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`,
@@ -392,14 +391,14 @@ def Ip_in_units_of_B(Ip: Union[None, float, np.ndarray]=None,
             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 
@@ -418,12 +417,13 @@ def Ip_in_units_of_B(Ip: Union[None, float, np.ndarray]=None,
             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':
@@ -455,7 +455,7 @@ def h_to_fluxbias(h: Union[float, np.ndarray]=1,
                   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.
+    """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.
@@ -475,14 +475,14 @@ def h_to_fluxbias(h: Union[float, np.ndarray]=1,
             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 
@@ -501,6 +501,7 @@ def h_to_fluxbias(h: Union[float, np.ndarray]=1,
             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.
@@ -515,7 +516,7 @@ def fluxbias_to_h(fluxbias: Union[float, np.ndarray]=1,
                   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``.
+    """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. 
@@ -545,7 +546,7 @@ def fluxbias_to_h(fluxbias: Union[float, np.ndarray]=1,
             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 
@@ -563,7 +564,7 @@ def fluxbias_to_h(fluxbias: Union[float, np.ndarray]=1,
         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.
     """
@@ -574,10 +575,11 @@ def fluxbias_to_h(fluxbias: Union[float, np.ndarray]=1,
 
 
 def freezeout_effective_temperature(freezeout_B, temperature, units_B = 'GHz', units_T = 'mK') -> float:
-    r'''Provides an effective temperature as a function of freezeout information.
-
+    '''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 
@@ -587,7 +589,7 @@ def freezeout_effective_temperature(freezeout_B, temperature, units_B = 'GHz', u
     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
 
@@ -598,7 +600,7 @@ def freezeout_effective_temperature(freezeout_B, temperature, units_B = 'GHz', u
 
     .. 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.
@@ -623,7 +625,7 @@ def freezeout_effective_temperature(freezeout_B, temperature, units_B = 'GHz', u
 
         temperature (float):
             :math:`T`, the physical temperature of the quantum computer.
-
+        
         units_B (string, optional, 'GHz'):
             Units in which ``freezeout_B`` is specified. Allowed values:
             'GHz' (Giga-Hertz) and 'J' (Joules).
@@ -632,21 +634,22 @@ def freezeout_effective_temperature(freezeout_B, temperature, units_B = 'GHz', u
             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` 
@@ -654,7 +657,7 @@ def freezeout_effective_temperature(freezeout_B, temperature, units_B = 'GHz', u
         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
@@ -665,7 +668,7 @@ def freezeout_effective_temperature(freezeout_B, temperature, units_B = 'GHz', u
     else:
         raise ValueException("Units must be 'J' (Joules) "
                              "or 'mK' (milli-Kelvin)")
-
+    
     if units_T == 'mK':
         temperature = temperature * 1e-3
     elif units_T == 'K':
@@ -674,15 +677,15 @@ def freezeout_effective_temperature(freezeout_B, temperature, units_B = 'GHz', u
         raise ValueException("Units must be 'K' (Kelvin) "
                              "or 'mK' (milli-Kelvin)")
     kB = 1.3806503e-23 # J/K
-
+    
     return 2*temperature*kB/freezeout_B
 
 def fast_effective_temperature(sampler=None, num_reads=None, seed=None,
                                h_range=(-1/6.1,1/6.1), sampler_params=None,
                                optimize_method=None,
                                num_bootstrap_samples=0) -> Tuple[np.float64,np.float64]:
-    r'''Provides an estimate to the effective temperature, :math:`T`, of a sampler.
-
+    '''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.
 
@@ -706,7 +709,7 @@ def fast_effective_temperature(sampler=None, num_reads=None, seed=None,
             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`` 
@@ -727,23 +730,23 @@ def fast_effective_temperature(sampler=None, num_reads=None, seed=None,
             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 
@@ -752,23 +755,24 @@ def fast_effective_temperature(sampler=None, num_reads=None, seed=None,
         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 = {}
@@ -793,7 +797,7 @@ def fast_effective_temperature(sampler=None, num_reads=None, seed=None,
 
     if num_bootstrap_samples is None:
         num_bootstrap_samples = sampler_params0['num_reads'] 
-
+        
     sampleset = sampler.sample(bqm, **sampler_params0)
 
     T,Tboot = maximum_pseudolikelihood_temperature(
@@ -802,6 +806,7 @@ def fast_effective_temperature(sampler=None, num_reads=None, seed=None,
         optimize_method=optimize_method,
         num_bootstrap_samples=num_bootstrap_samples)
 
+    
     if num_bootstrap_samples == 0:
         return T, np.float64(0.0)
     else:
diff --git a/releasenotes/notes/add-linear-ancilla-composite-3281ed6733b0f0c7.yaml b/releasenotes/notes/add-linear-ancilla-composite-3281ed6733b0f0c7.yaml
new file mode 100644
index 0000000..ab7ecf1
--- /dev/null
+++ b/releasenotes/notes/add-linear-ancilla-composite-3281ed6733b0f0c7.yaml
@@ -0,0 +1,4 @@
+---
+features:
+  - |
+    Add `LinearAncillaComposite` for implementing linear coefficients through ancilla qubits polarized with strong flux biases 
diff --git a/requirements.txt b/requirements.txt
index bc0a3e7..c06f72b 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,15 +1,18 @@
-dimod==0.12.18
-dwave-optimization==0.4.0
-dwave-preprocessing==0.6.7
+--extra-index-url https://pypi.dwavesys.com/simple
+
+dimod==0.12.13
+dwave-optimization==0.1.0
+dwave-preprocessing==0.6.4
 dwave-cloud-client==0.12.0
-dwave-networkx==0.8.14
-dwave-samplers==1.4.0
+dwave-networkx==0.8.10
+dwave-drivers==0.4.4
+dwave-samplers==1.2.0
 homebase==1.0.1
-minorminer==0.2.16
-networkx~=3.0
+minorminer==0.2.14
+networkx==2.8.8
 
 numpy==1.23.5;python_version<"3.11"  # for compatibility with scipy 1.9.3
-numpy==2.1.0;python_version>="3.13"
+oldest-supported-numpy;python_version>="3.11"
 
 scipy==1.9.3;python_version<"3.11"
-scipy==1.14.1;python_version>="3.11"
+scipy==1.11.2;python_version>="3.11"
diff --git a/setup.py b/setup.py
index f701582..b5f1495 100644
--- a/setup.py
+++ b/setup.py
@@ -24,9 +24,9 @@ os.chdir(setup_folder_loc)
 exec(open(os.path.join(".", "dwave", "system", "package_info.py")).read())
 
 
-install_requires = ['dimod>=0.12.18,<0.14.0',
-                    'dwave-optimization>=0.1.0,<0.6',
-                    'dwave-cloud-client>=0.12.0,<0.14.0',
+install_requires = ['dimod>=0.12.7,<0.14.0',
+                    'dwave-optimization>=0.1.0,<0.3',
+                    'dwave-cloud-client>=0.12.0,<0.13.0',
                     'dwave-networkx>=0.8.10',
                     'dwave-preprocessing>=0.5.0',
                     'homebase>=1.0.0,<2.0.0',
@@ -36,11 +36,18 @@ install_requires = ['dimod>=0.12.18,<0.14.0',
                     'scipy>=1.7.3',
                     ]
 
-python_requires = '>=3.9'
+# NOTE: dwave-drivers can also be installed with `dwave install drivers`,
+# and the exact requirements for `drivers` contrib package are defined in
+# `dwave.system.package_info`, the `contrib` dict.
+extras_require = {'drivers': ['dwave-drivers>=0.4.0,<0.5.0'],
+                  }
+
+python_requires = '>=3.8'
 
 packages = ['dwave',
             'dwave.embedding',
             'dwave.system',
+            'dwave.system.cache',
             'dwave.system.composites',
             'dwave.system.samplers',
             ]
@@ -50,11 +57,11 @@ classifiers = [
     'Operating System :: OS Independent',
     'Development Status :: 3 - Alpha',
     'Programming Language :: Python :: 3 :: Only',
+    'Programming Language :: Python :: 3.8',
     'Programming Language :: Python :: 3.9',
     'Programming Language :: Python :: 3.10',
     'Programming Language :: Python :: 3.11',
     'Programming Language :: Python :: 3.12',
-    'Programming Language :: Python :: 3.13',
 ]
 
 setup(
@@ -67,7 +74,13 @@ setup(
     url='https://github.com/dwavesystems/dwave-system',
     license='Apache 2.0',
     packages=packages,
+    entry_points={
+        'dwave_contrib': [
+            'dwave-system = dwave.system.package_info:contrib'
+        ]
+    },
     install_requires=install_requires,
+    extras_require=extras_require,
     python_requires=python_requires,
     classifiers=classifiers,
     zip_safe=False


Test patch:
diff --git a/dwave/system/testing.py b/dwave/system/testing.py
index 810fa3a..9c5fec8 100644
--- a/dwave/system/testing.py
+++ b/dwave/system/testing.py
@@ -28,15 +28,9 @@ from dwave.system import qpu_graph
 class MockDWaveSampler(dimod.Sampler, dimod.Structured):
     """Mock sampler modeled after DWaveSampler that can be used for tests.
 
-    Properties and topology parameters are populated to qualitatively match
+    Properties and topology parameters are populated qualitatively matching
     online systems, and a placeholder sampler routine based on steepest descent
-    is instantiated by default.
-
-    The :attr:`.EXACT_SOLVER_CUTOFF_DEFAULT` defines the problem size threshold for using the exact solver.
-    For problems with fewer variables than this threshold, the exact ground state is computed 
-    using a brute-force solver. This provides a reproducible solution for small problem sizes.
-
-    For larger problems, the `SteepestDescentSampler` is used as a placeholder solver.
+    is instantiated.
 
     Args:
         nodelist (iterable of ints, optional):
@@ -79,36 +73,16 @@ class MockDWaveSampler(dimod.Sampler, dimod.Structured):
             parameters. By default ``initial_state`` can also be mocked, if
             dwave-greedy is installed. All other parameters are ignored and a 
             warning will be raised by default.
-        
-        substitute_sampler (:class:`~dimod.Sampler`, optional, default=SteepestDescentSampler()):
-            The sampler to be used as a substitute when executing the mock sampler. 
-            By default, :class:`~dwave.samplers.SteepestDescentSampler` is employed, which performs a 
-            deterministic steepest descent optimization on the BQM. Supported options are
-            any dimod-compatible sampler to customize the sampling behavior of 
-            `MockDWaveSampler()`.
-
-        substitute_kwargs (dict, optional, default={}):
-            A dictionary of keyword arguments to pass to the `substitute_sampler`'s 
-            `sample` method. This allows users to configure the substitute sampler 
-            with specific parameters like `num_reads`, `initial_state`, or other 
-            sampler-specific options. If not provided, an empty dictionary is used 
-            by default.
 
         exact_solver_cutoff (int, optional, default=:attr:`EXACT_SOLVER_CUTOFF_DEFAULT`):
             For problems smaller or equal in size to ``exact_solver_cutoff``, the
             first sample in any sampleset returned by the sampling routines
             is replaced by a reproducible ground state (determined exactly with
-            a brute-force :class:`~dimod.ExactSolver`). Only small cutoffs
+            a brute-force :class:`~dimod.ExactSolver`). Only small cut offs
             should be used since solution time increases exponentially with
             problem size.
-            
-            - When ``substitute_sampler`` is not provided, the default value is 
-            ``EXACT_SOLVER_CUTOFF_DEFAULT`` (e.g., 16).
-            - When ``substitute_sampler`` is provided, the default value is 
-            ``0``, disabling exact ground state calculation.
-            
             Set ``exact_solver_cutoff`` to zero to disable exact ground state
-            calculation explicitly.
+            calculation.
 
     Examples
         The first example creates a MockSampler without reference to a
@@ -134,43 +108,23 @@ class MockDWaveSampler(dimod.Sampler, dimod.Structured):
         -1
 
     """
+    # Feature suggestion - add seed as an optional input, to allow reproducibility.
+
     nodelist = None
     edgelist = None
     properties = None
     parameters = None
 
+    # by default, use ExactSolver for problems up to size (inclusive):
+    EXACT_SOLVER_CUTOFF_DEFAULT = 16
+
     def __init__(self,
                  nodelist=None, edgelist=None, properties=None,
                  broken_nodes=None, broken_edges=None,
                  topology_type=None, topology_shape=None,
                  parameter_warnings=True,
-                 substitute_sampler=None, 
-                 substitute_kwargs=None, 
-                 exact_solver_cutoff=None,
+                 exact_solver_cutoff=EXACT_SOLVER_CUTOFF_DEFAULT,
                  **config):
-        
-        self.mocked_parameters={'answer_mode',
-                                'max_answers',
-                                'num_reads',
-                                'label',
-                                'initial_state'}
-
-        EXACT_SOLVER_CUTOFF_DEFAULT = 16
-
-        if substitute_sampler is None:
-            substitute_sampler = SteepestDescentSampler()
-            if exact_solver_cutoff is None:
-                exact_solver_cutoff = EXACT_SOLVER_CUTOFF_DEFAULT
-        else:
-            if exact_solver_cutoff is None:
-                exact_solver_cutoff = 0
-
-        self.substitute_sampler = substitute_sampler
-
-        if substitute_kwargs is None:
-            substitute_kwargs = {} 
-        self.substitute_kwargs = substitute_kwargs
-
         self.parameter_warnings = parameter_warnings
         self.exact_solver_cutoff = exact_solver_cutoff
 
@@ -204,6 +158,7 @@ class MockDWaveSampler(dimod.Sampler, dimod.Structured):
             'chip_id': 'MockDWaveSampler',
             'topology': {'type': topology_type, 'shape': topology_shape}
         }
+
         #Create graph object, introduce defects per input arguments
         if nodelist is not None:
             self.nodelist = nodelist.copy()
@@ -216,11 +171,6 @@ class MockDWaveSampler(dimod.Sampler, dimod.Structured):
                                  self.properties['topology']['shape'],
                                  self.nodelist, self.edgelist)
 
-        if topology_type == 'pegasus':
-            m = self.properties['topology']['shape'][0]
-            num_qubits = m*(m-1)*24  # fabric_only=True technicality
-        else:
-            num_qubits = len(solver_graph)
         if broken_nodes is None and broken_edges is None:
             self.nodelist = sorted(solver_graph.nodes)
             self.edgelist = sorted(tuple(sorted(edge))
@@ -239,7 +189,7 @@ class MockDWaveSampler(dimod.Sampler, dimod.Structured):
                                    and (v, u) not in broken_edges)
         #Finalize yield-dependent properties:
         self.properties.update({
-            'num_qubits': num_qubits,
+            'num_qubits': len(solver_graph),
             'qubits': self.nodelist.copy(),
             'couplers': self.edgelist.copy(),
             'anneal_offset_ranges': [[-0.5, 0.5] if i in self.nodelist
@@ -330,7 +280,7 @@ class MockDWaveSampler(dimod.Sampler, dimod.Structured):
             'default_programming_thermalization': 1000.0,
             'default_readout_thermalization': 0.0,
             'extended_j_range': [-2.0, 1.0],
-            'fast_anneal_time_range': [0.005, 2000.0],
+            'fast_anneal_time_range': [0.005, 83000.0],
             'h_gain_schedule_range': [-3.0, 3.0],
             'max_anneal_schedule_points': 12,
             'max_h_gain_schedule_points': 20,
@@ -363,9 +313,15 @@ class MockDWaveSampler(dimod.Sampler, dimod.Structured):
     def sample(self, bqm, **kwargs):
 
         # Check kwargs compatibility with parameters and substitute sampler:
+        mocked_parameters={'answer_mode',
+                           'max_answers',
+                           'num_reads',
+                           'label',
+                           'initial_state',
+                           }
         for kw in kwargs:
             if kw in self.parameters:
-                if self.parameter_warnings and kw not in self.mocked_parameters:
+                if self.parameter_warnings and kw not in mocked_parameters:
                     warnings.warn(f'{kw!r} parameter is valid for DWaveSampler(), '
                                   'but not mocked in MockDWaveSampler().')
             else:
@@ -388,21 +344,19 @@ class MockDWaveSampler(dimod.Sampler, dimod.Structured):
         label = kwargs.get('label')
         if label is not None:
             info.update(problem_label=label)
-        
-        # Special handling of flux_biases, for compatibility with virtual graphs
+
+        #Special handling of flux_biases, for compatibility with virtual graphs
+
         flux_biases = kwargs.get('flux_biases')
         if flux_biases is not None:
             self.flux_biases_flag = True
 
-        # Create a local dictionary combining self.substitute_kwargs and relevant kwargs
-        substitute_kwargs = self.substitute_kwargs.copy()
+        substitute_kwargs = {'num_reads' : kwargs.get('num_reads')}
+        if substitute_kwargs['num_reads'] is None:
+            substitute_kwargs['num_reads'] = 1
 
-        # Handle 'num_reads', defaulting to 1 if not provided
-        num_reads = kwargs.get('num_reads', substitute_kwargs.get('num_reads', 1))
-        substitute_kwargs['num_reads'] = num_reads
-
-        if 'initial_state' in kwargs:
-            initial_state = kwargs['initial_state']
+        initial_state = kwargs.get('initial_state')
+        if initial_state is not None:
             # Initial state format is a list of (qubit,values)
             # value=3 denotes an unused variable (should be absent
             # from bqm). 
@@ -412,17 +366,15 @@ class MockDWaveSampler(dimod.Sampler, dimod.Structured):
                           if pair[1]!=3],dtype=float),
                 [pair[0] for pair in initial_state if pair[1]!=3])
 
-        sampler_kwargs = kwargs.copy()
-        sampler_kwargs.update(substitute_kwargs)
-
-        ss = self.substitute_sampler.sample(bqm, **sampler_kwargs)
+        ss = SteepestDescentSampler().sample(bqm, **substitute_kwargs)
         ss.info.update(info)
+
         # determine ground state exactly for small problems
         if 0 < len(bqm) <= self.exact_solver_cutoff and len(ss) >= 1:
             ground = dimod.ExactSolver().sample(bqm).truncate(1)
             ss.record[0].sample = ground.record[0].sample
             ss.record[0].energy = ground.record[0].energy
-        
+
         answer_mode = kwargs.get('answer_mode')
         if answer_mode is None or answer_mode == 'histogram':
             # Default for DWaveSampler() is 'histogram'
diff --git a/tests/qpu/test_dwavecliquesampler.py b/tests/qpu/test_dwavecliquesampler.py
index 9392a9a..183e466 100644
--- a/tests/qpu/test_dwavecliquesampler.py
+++ b/tests/qpu/test_dwavecliquesampler.py
@@ -15,6 +15,7 @@
 import itertools
 import os
 import unittest
+import warnings
 
 import dimod
 
diff --git a/tests/qpu/test_dwavesampler.py b/tests/qpu/test_dwavesampler.py
index 31f0e15..0a07747 100644
--- a/tests/qpu/test_dwavesampler.py
+++ b/tests/qpu/test_dwavesampler.py
@@ -21,6 +21,7 @@ import numpy
 import dimod
 from dwave.cloud.exceptions import ConfigFileError, SolverNotFoundError
 from dwave.cloud.client import Client
+from dwave_networkx.generators.pegasus import pegasus_graph
 
 from dwave.system.samplers import DWaveSampler
 
@@ -168,14 +169,14 @@ class TestMissingQubits(unittest.TestCase):
 class TestClientSelection(unittest.TestCase):
 
     def test_client_type(self):
-        with mock.patch('dwave.cloud.client.qpu.Client') as qpu:
+        with mock.patch('dwave.cloud.qpu.Client') as qpu:
             self.assertEqual(DWaveSampler().client, qpu())
             self.assertEqual(DWaveSampler(client='qpu').client, qpu())
 
-        with mock.patch('dwave.cloud.client.sw.Client') as sw:
+        with mock.patch('dwave.cloud.sw.Client') as sw:
             self.assertEqual(DWaveSampler(client='sw').client, sw())
 
-        with mock.patch('dwave.cloud.client.hybrid.Client') as hybrid:
+        with mock.patch('dwave.cloud.hybrid.Client') as hybrid:
             self.assertEqual(DWaveSampler(client='hybrid').client, hybrid())
 
     def test_base_client(self):
diff --git a/tests/qpu/test_leaphybrid_common.py b/tests/qpu/test_leaphybrid_common.py
index 93a9adc..0b42563 100644
--- a/tests/qpu/test_leaphybrid_common.py
+++ b/tests/qpu/test_leaphybrid_common.py
@@ -19,7 +19,6 @@ import unittest
 
 from parameterized import parameterized_class
 
-import dimod
 from dwave.cloud.exceptions import ConfigFileError, SolverNotFoundError
 from dwave.cloud.testing import isolated_environ
 from dwave.system import LeapHybridSampler, LeapHybridDQMSampler, LeapHybridCQMSampler
@@ -84,32 +83,3 @@ class TestLegacySolverSelection(unittest.TestCase):
             self.assertIn(self.problem_type, sampler.solver.supported_problem_types)
 
         sampler.client.close()
-
-
-@parameterized_class(
-    ("sampler_cls", "sample_meth", "problem_gen"), [
-        (LeapHybridSampler, "sample", lambda self: dimod.BQM.from_qubo({})),
-        (LeapHybridDQMSampler, "sample_dqm", lambda self: dimod.DQM.from_numpy_vectors([0], [0], ([], [], []))),
-        (LeapHybridCQMSampler, "sample_cqm", lambda self: dimod.CQM.from_bqm(dimod.BQM.from_qubo({'ab': 1}))),
-    ])
-@unittest.skipIf(os.getenv('SKIP_INT_TESTS'), "Skipping integration test.")
-class TestSamplesetInterface(unittest.TestCase):
-
-    def test_wait_id_availability(self):
-        # verify https://github.com/dwavesystems/dwave-system/issues/540 is fixed
-        try:
-            sampler = self.sampler_cls()
-        except (ValueError, ConfigFileError):
-            raise unittest.SkipTest(f"{self.sampler_cls} not available")
-
-        problem = self.problem_gen()
-        ss = getattr(sampler, self.sample_meth)(problem)
-
-        with self.subTest("sampleset.wait_id() exists"):
-            pid = ss.wait_id()
-            self.assertIsInstance(pid, str)
-
-        with self.subTest("sampleset.wait_id() exists post-resolve"):
-            ss.resolve()
-            pid_post = ss.wait_id()
-            self.assertEqual(pid, pid_post)
diff --git a/tests/qpu/test_leaphybriddqmsampler.py b/tests/qpu/test_leaphybriddqmsampler.py
index f0b2815..27414a9 100644
--- a/tests/qpu/test_leaphybriddqmsampler.py
+++ b/tests/qpu/test_leaphybriddqmsampler.py
@@ -64,6 +64,9 @@ class TestLeapHybridSampler(unittest.TestCase):
         np.testing.assert_array_almost_equal(dqm.energies(sampleset),
                                              sampleset.record.energy)
 
+    # NOTE: enable when problem labelling deployed to prod
+    @unittest.skipIf(sampler is not None and 'cloud' in sampler.client.endpoint,
+                     "labels not supported in production")
     def test_problem_labelling(self):
         dqm = dimod.DQM()
         u = dqm.add_variable(2, 'a')
diff --git a/tests/qpu/test_leaphybridsampler.py b/tests/qpu/test_leaphybridsampler.py
index 3bd1289..10d3112 100644
--- a/tests/qpu/test_leaphybridsampler.py
+++ b/tests/qpu/test_leaphybridsampler.py
@@ -35,6 +35,9 @@ class TestLeapHybridSampler(unittest.TestCase):
         sampleset = sampler.sample_ising({'a': -1}, {'ab': 1})
         sampleset.resolve()
 
+    # NOTE: enable when problem labelling deployed to prod
+    @unittest.skipIf(sampler is not None and 'cloud' in sampler.client.endpoint,
+                     "labels not supported in production")
     def test_problem_labelling(self):
         bqm = dimod.BQM.from_ising({'a': -1}, {'ab': 1})
         label = 'problem label'
diff --git a/tests/test_database_manager.py b/tests/test_database_manager.py
new file mode 100644
index 0000000..2e67d98
--- /dev/null
+++ b/tests/test_database_manager.py
@@ -0,0 +1,173 @@
+# Copyright 2018 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.
+
+import unittest
+import time
+import homebase
+
+import networkx as nx
+
+from dwave.system.exceptions import MissingFluxBias
+import dwave.system.cache as cache
+
+tmp_database_name = 'tmp_test_database_manager_{}.db'.format(time.time())
+# test_database_path = cache.cache_file(filename=tmp_database_name)
+# conn = cache.cache_connect(test_database_path)
+
+
+class TestCacheManager(unittest.TestCase):
+    def test_same_database(self):
+        """multiple calls to get_database_path"""
+        db1 = cache.cache_file()
+        db2 = cache.cache_file()
+
+        self.assertEqual(db1, db2)
+
+
+class TestDatabaseManager(unittest.TestCase):
+    def setUp(self):
+        # new connection for just this test
+        self.clean_connection = cache.cache_connect(':memory:')
+
+        # test_database_path = cache.cache_file(filename=tmp_database_name)
+        # test_connection = cache.cache_connect(test_database_path)
+
+    def tearDown(self):
+        # close the memory connection
+        self.clean_connection.close()
+
+    def test_insert_retrieve_chain(self):
+        """insert and retrieve some chains"""
+        conn = self.clean_connection
+
+        chain = [0, 1, 3]
+
+        # insert the chain twice, should only result in one net insert
+        with conn as cur:
+            cache.insert_chain(cur, chain)
+            cache.insert_chain(cur, chain)
+
+        with conn as cur:
+            all_chains = list(cache.iter_chain(cur))
+        self.assertEqual(len(all_chains), 1)
+        returned_chain, = all_chains
+        self.assertEqual(len(returned_chain), len(chain))
+        self.assertEqual(returned_chain, chain)
+
+        # insert a new chain
+        with conn as cur:
+            cache.insert_chain(cur, [0, 2, 3])
+        with conn as cur:
+            all_chains = list(cache.iter_chain(cur))
+        self.assertEqual(len(all_chains), 2)  # should now be two chains
+
+    def test_insert_flux_bias(self):
+        conn = self.clean_connection
+
+        with conn as cur:
+            cache.insert_flux_bias(cur, [0, 1, 2], 'test_system', .1, 1)
+
+        with conn as cur:
+            # try dumping
+            flux_biases = list(cache.iter_flux_bias(cur))
+        self.assertEqual(len(flux_biases), 1)
+        self.assertEqual(flux_biases[0], ([0, 1, 2], 'test_system', .1, 1))
+
+        # retrieve by name
+        with conn as cur:
+            biases = cache.get_flux_biases_from_cache(cur, [[0, 1, 2]], 'test_system', 1)
+
+        for v, fbo in biases.items():
+            self.assertIn(v, [0, 1, 2])
+            self.assertEqual(fbo, .1)
+
+        # now get something wrong out
+        with self.assertRaises(MissingFluxBias):
+            with conn as cur:
+                biases = cache.get_flux_biases_from_cache(cur, [[0, 1, 2]], 'another_system', 1)
+
+    def test_graph_insert_retrieve(self):
+        conn = self.clean_connection
+
+        graph = nx.barbell_graph(8, 8)
+        nodelist = sorted(graph)
+        edgelist = sorted(sorted(edge) for edge in graph.edges)
+
+        with conn as cur:
+            cache.insert_graph(cur, nodelist, edgelist)
+
+            # should only be one graph
+            graphs = list(cache.iter_graph(cur))
+            self.assertEqual(len(graphs), 1)
+            (nodelist_, edgelist_), = graphs
+            self.assertEqual(nodelist, nodelist_)
+            self.assertEqual(edgelist, edgelist_)
+
+        # trying to reinsert should still result in only one graph
+        with conn as cur:
+            cache.insert_graph(cur, nodelist, edgelist)
+            graphs = list(cache.iter_graph(cur))
+            self.assertEqual(len(graphs), 1)
+
+        # inserting with an empty dict as encoded_data should populate it
+        encoded_data = {}
+        with conn as cur:
+            cache.insert_graph(cur, nodelist, edgelist, encoded_data)
+        self.assertIn('num_nodes', encoded_data)
+        self.assertIn('num_edges', encoded_data)
+        self.assertIn('edges', encoded_data)
+
+        # now adding another graph should result in two items
+        graph = nx.complete_graph(4)
+        nodelist = sorted(graph)
+        edgelist = sorted(sorted(edge) for edge in graph.edges)
+        with conn as cur:
+            cache.insert_graph(cur, nodelist, edgelist)
+            graphs = list(cache.iter_graph(cur))
+            self.assertEqual(len(graphs), 2)
+
+    def test_insert_embedding(self):
+        test_database_path = cache.cache_file(filename=tmp_database_name)
+        conn = cache.cache_connect(test_database_path)
+
+        source_nodes = [0, 1, 2]
+        source_edges = [[0, 1], [0, 2], [1, 2]]
+        target_nodes = [0, 1, 2, 3]
+        target_edges = [[0, 1], [0, 3], [1, 2], [2, 3]]
+
+        embedding = {0: [0], 1: [1], 2: [2, 3]}
+
+        with conn as cur:
+            cache.insert_embedding(cur, source_nodes, source_edges, target_nodes, target_edges,
+                                     embedding, 'tag1')
+
+            embedding_ = cache.select_embedding_from_tag(cur, 'tag1', target_nodes, target_edges)
+
+            self.assertEqual(embedding, embedding_)
+
+        # now reinsert but with a different embedding
+        embedding = {0: [0, 1], 1: [2], 2: [3]}
+        with conn as cur:
+            cache.insert_embedding(cur, source_nodes, source_edges, target_nodes, target_edges,
+                                     embedding, 'tag1')
+
+            # get it back
+            embedding_ = cache.select_embedding_from_tag(cur, 'tag1', target_nodes, target_edges)
+
+            self.assertEqual(embedding, embedding_)
+
+            # get it back from source graph
+            embedding_ = cache.select_embedding_from_source(cur, source_nodes, source_edges,
+                                                              target_nodes, target_edges)
+            self.assertEqual(embedding, embedding_)
diff --git a/tests/test_linear_ancilla_composite.py b/tests/test_linear_ancilla_composite.py
new file mode 100644
index 0000000..91816ba
--- /dev/null
+++ b/tests/test_linear_ancilla_composite.py
@@ -0,0 +1,119 @@
+#  Copyright 2024 D-Wave 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.
+
+import collections
+import unittest
+
+import dimod
+import networkx as nx
+import numpy as np
+
+from dwave.system.testing import MockDWaveSampler
+from dwave.system import LinearAncillaComposite
+
+
+class TestLinearAncillaComposite(unittest.TestCase):
+    def setUp(self):
+        self.qpu = MockDWaveSampler(properties=dict(extended_j_range=[-2, 1]))
+        self.tracked_qpu = dimod.TrackingComposite(self.qpu)
+
+        self.sampler = LinearAncillaComposite(
+            dimod.StructureComposite(
+                self.tracked_qpu,
+                nodelist=self.qpu.nodelist,
+                edgelist=self.qpu.edgelist,
+            )
+        )
+
+        self.submask = nx.subgraph(
+            self.sampler.to_networkx_graph(),
+            list(self.sampler.nodelist)[::2],
+        )
+
+        # this problem should run
+        self.linear_problem = dimod.BinaryQuadraticModel.from_ising(
+            {i: (-1) ** i for i in self.submask.nodes()},
+            {},
+        )
+
+        # this problem shouldn't run
+        self.linear_problem_full_graph = dimod.BinaryQuadraticModel.from_ising(
+            {i: (-1) ** i for i in self.qpu.nodelist},
+            {},
+        )
+
+    def test_only_quadratic(self):
+        """if no linear biases, the bqm remains intact"""
+
+        bqm = dimod.generators.ran_r(1, self.submask, seed=1)
+        self.sampler.sample(bqm)
+        self.assertEqual(bqm, self.tracked_qpu.input["bqm"])
+
+    def test_h_tolerance_too_large(self):
+        """if h tolerance is larger than the linear biases,
+        the bqm remains intact
+        """
+
+        self.sampler.sample(self.linear_problem, h_tolerance=1.01)
+        self.assertEqual(self.linear_problem, self.tracked_qpu.input["bqm"])
+
+    def test_intermediate_h_tolerance(self):
+        """check the desired h-tolerance is left in the qubit bias"""
+
+        h_tolerance = 0.5
+        self.sampler.sample(self.linear_problem, h_tolerance=h_tolerance)
+        for variable, bias in self.tracked_qpu.input["bqm"].linear.items():
+            if variable in self.linear_problem.variables:  # skip the ancillas
+                self.assertEqual(
+                    bias,
+                    np.sign(self.linear_problem.get_linear(variable)) * h_tolerance,
+                )
+
+    def test_no_ancillas_available(self):
+        """send a problem that uses all the qubits, not leaving any ancillas available"""
+
+        with self.assertRaises(ValueError):
+            ss = self.sampler.sample(self.linear_problem_full_graph)
+
+    def test_ancillas_present(self):
+        """check the solver used ancillas"""
+
+        self.sampler.sample(self.linear_problem)
+        self.assertGreater(
+            len(self.tracked_qpu.input["bqm"].variables),
+            len(self.linear_problem.variables),
+        )
+
+    def test_ancilla_cleanup(self):
+        """check the problem returned has no additional variables"""
+
+        sampleset = self.sampler.sample(self.linear_problem)
+        self.assertEqual(
+            len(self.linear_problem.variables),
+            len(sampleset.variables),
+        )
+
+    def test_flux_biases_present(self):
+        """check flux biases are applied to non-data qubits"""
+
+        self.sampler.sample(self.linear_problem)
+        flux_biases = np.array(self.tracked_qpu.input["flux_biases"])
+
+        # flux biases are used
+        self.assertGreater(sum(flux_biases != 0), 0)
+
+        # the qubits with flux biases are not data qubits
+        for qubit, flux_bias in enumerate(flux_biases):
+            if flux_bias != 0:
+                self.assertNotIn(qubit, self.linear_problem.variables)
diff --git a/tests/test_mock_sampler.py b/tests/test_mock_sampler.py
index 87b5a31..4a052df 100644
--- a/tests/test_mock_sampler.py
+++ b/tests/test_mock_sampler.py
@@ -25,7 +25,6 @@ from dwave.system import DWaveSampler
 from dwave.system.testing import MockDWaveSampler, MockLeapHybridDQMSampler
 from dwave.cloud.exceptions import ConfigFileError, SolverNotFoundError
 from dimod import DiscreteQuadraticModel, ExtendedVartype, SampleSet
-from dwave.samplers import SteepestDescentSolver
 
 
 class TestMockDWaveSampler(unittest.TestCase):
@@ -217,94 +216,7 @@ class TestMockDWaveSampler(unittest.TestCase):
                                    broken_edges=delete_edges)
         self.assertTrue(len(sampler.nodelist)==4)
         self.assertTrue(len(sampler.edgelist)==1)
-
-    def test_custom_substitute_sampler(self):
-        """Test that MockDWaveSampler uses the provided custom substitute_sampler."""
-
-        # Define a sampler that always returns the a constant (excited) state
-        class SteepestAscentSolver(SteepestDescentSolver):
-            def sample(self, bqm, **kwargs):
-                # Return local (or global)  maxima instead of local minima
-                # NOTE: energy returned is not faithful to the original bqm (energy calculated as `-bqm`)
-                return super().sample(-bqm, **kwargs)
-
-        inverted_sampler = SteepestAscentSolver()
-
-        # Create a simple BQM
-        bqm = dimod.BQM({'a': 1, 'b': 1}, {}, 0.0, vartype="SPIN")
-
-        # Instantiate MockDWaveSampler with nodelist and edgelist including 'a' and 'b'
-        sampler = MockDWaveSampler(
-            substitute_sampler=inverted_sampler,
-            nodelist=['a', 'b'],
-            edgelist=[('a', 'b')]
-        )
-
-        # First Subtest: First sample does not use ExactSampler();
-        # Second sample does not use SteepestDescentSampler()
-        with self.subTest("Sampler without ExactSampler"):
-            ss = sampler.sample(bqm, num_reads=2)
-            self.assertEqual(sampler.exact_solver_cutoff, 0)
-            self.assertEqual(ss.record.sample.shape, (1,2), 'Unique sample expected')
-            self.assertTrue(np.all(ss.record.sample==1), 'Excited states expected')
-
-        sampler = MockDWaveSampler(
-            substitute_sampler=inverted_sampler,
-            nodelist=['a', 'b'],
-            edgelist=[('a', 'b')],
-            exact_solver_cutoff=2
-        )
-        # Second Subtest: First sample uses ExactSampler();
-        # Second sampler uses inverted sampler. Explicit exact_solver_cutoff overrides substitute_sampler.
-        with self.subTest("Sampler with ExactSampler and substitute_sampler"):
-            ss = sampler.sample(bqm, num_reads=2, answer_mode='raw')
-            self.assertEqual(sampler.exact_solver_cutoff, 2)
-            self.assertEqual(ss.record.sample.shape, (2,2), 'Non-unique samples expected')
-            self.assertTrue(np.all(ss.record.sample[0,:] == -1), 'Excited states expected')
-            self.assertTrue(np.all(ss.record.sample[1,:] == 1), 'Excited states expected')
-
-    def test_mocking_sampler_params(self):
-        """Test that substitute_kwargs are correctly passed to the substitute_sampler."""
-
-        # Define a constant sampler that checks for a custom parameter
-        class ConstantSampler(dimod.Sampler):
-            properties = {}
-            parameters = {'custom_param': [], 'num_reads': []}
-
-            def sample(self, bqm, **kwargs):
-                custom_param = kwargs.get('custom_param')
-                num_reads = kwargs.get('num_reads')
-                # Raise exception if parameters passed incorrectly
-                if custom_param != 'test_value':
-                    raise ValueError("custom_param not passed correctly")
-                if num_reads != 10:
-                    raise ValueError(f"num_reads not passed correctly, expected 10, got {num_reads}")
-                # Return a default sample
-                sample = {v: -1 for v in bqm.variables}
-                return dimod.SampleSet.from_samples_bqm(sample, bqm)
-
-        constant_sampler = ConstantSampler()
-
-        # Create a simple BQM
-        bqm = dimod.BQM({'a': 1, 'b': 1}, {('a', 'b'): 1}, 0.0, vartype="SPIN")
-
-        # Instantiate MockDWaveSampler with nodelist and edgelist including 'a' and 'b'
-        sampler = MockDWaveSampler(
-            substitute_sampler=constant_sampler,
-            substitute_kwargs={'custom_param': 'test_value'},
-            nodelist=['a', 'b'],
-            edgelist=[('a', 'b')]
-        )
-
-        # Sample using the MockDWaveSampler
-        ss = sampler.sample(bqm, num_reads=10)
-
-        # Check that the sample returned is as expected from the custom sampler
-        expected_sample = {'a': -1, 'b': -1}
-        self.assertEqual(ss.first.sample, expected_sample)
-        self.assertEqual(ss.first.energy, bqm.energy(expected_sample))
-
-
+        
 class TestMockLeapHybridDQMSampler(unittest.TestCase):
     def test_sampler(self):
         sampler = MockLeapHybridDQMSampler()
diff --git a/tests/test_virtual_graph_composite.py b/tests/test_virtual_graph_composite.py
index 0993928..5cfcefd 100644
--- a/tests/test_virtual_graph_composite.py
+++ b/tests/test_virtual_graph_composite.py
@@ -28,22 +28,20 @@ class TestVirtualGraphWithMockDWaveSampler(unittest.TestCase):
 
     def test_smoke(self):
         child_sampler = MockDWaveSampler()
-        with self.assertWarns(DeprecationWarning):
-            sampler = VirtualGraphComposite(child_sampler, {'a': [0]}, flux_bias_num_reads=1)
+        sampler = VirtualGraphComposite(child_sampler, {'a': [0]}, flux_bias_num_reads=1)
 
         # depending on how recenlty flux bias data was gathered, this may be true
         child_sampler.flux_biases_flag = False
-        child_sampler.mocked_parameters.add('flux_biases')  # Don't raise warning
+
         if sampler.flux_biases:
             sampler.sample_ising({'a': -1}, {})
             self.assertTrue(child_sampler.flux_biases_flag)  # true when some have been provided to sample_ising
 
     def test_structure_keyword_setting(self):
-        with self.assertWarns(DeprecationWarning):
-            sampler = VirtualGraphComposite(self.sampler, embedding={'a': set(range(8)),
-                                                                     'b': set(range(8, 16)),
-                                                                     'c': set(range(16, 24))},
-                                            flux_biases=False)
+        sampler = VirtualGraphComposite(self.sampler, embedding={'a': set(range(8)),
+                                                                 'b': set(range(8, 16)),
+                                                                 'c': set(range(16, 24))},
+                                        flux_biases=False)
 
         nodelist, edgelist, adj = sampler.structure
         self.assertEqual(nodelist, ['a', 'b', 'c'])
@@ -51,11 +49,10 @@ class TestVirtualGraphWithMockDWaveSampler(unittest.TestCase):
         self.assertEqual(adj, {'a': {'b'}, 'b': {'a', 'c'}, 'c': {'b'}})
 
         # unlike variable names
-        with self.assertWarns(DeprecationWarning):
-            sampler = VirtualGraphComposite(self.sampler, embedding={'a': set(range(8)),
-                                                                     1: set(range(8, 16)),
-                                                                     'c': set(range(16, 24))},
-                                            flux_biases=False)
+        sampler = VirtualGraphComposite(self.sampler, embedding={'a': set(range(8)),
+                                                                 1: set(range(8, 16)),
+                                                                 'c': set(range(16, 24))},
+                                        flux_biases=False)
         nodelist, edgelist, adj = sampler.structure
         self.assertEqual(set(nodelist), {'a', 1, 'c'})
         self.assertEqual(adj, {'a': {1}, 1: {'a', 'c'}, 'c': {1}})
@@ -78,19 +75,54 @@ class TestVirtualGraphWithMockDWaveSampler(unittest.TestCase):
         __, __, adj = sampler.structure
         embedding = {v: (v,) for v in adj}
 
-        with self.assertWarns(DeprecationWarning):
-            sampler = VirtualGraphComposite(sampler, embedding=embedding, flux_biases=False)
+        sampler = VirtualGraphComposite(sampler, embedding=embedding, flux_biases=False)
 
         self.assertEqual(sampler.embedding, embedding)
 
     def test_simple_complete_graph_sample_ising(self):
         """sample_ising on a K4."""
 
-        with self.assertWarns(DeprecationWarning):
-            K4 = VirtualGraphComposite(self.sampler, embedding={0: {0, 4},
-                                                                1: {1, 5},
-                                                                2: {2, 6},
-                                                                3: {3, 7}},
-                                       flux_biases=False)
+        K4 = VirtualGraphComposite(self.sampler, embedding={0: {0, 4},
+                                                            1: {1, 5},
+                                                            2: {2, 6},
+                                                            3: {3, 7}},
+                                   flux_biases=False)
 
         K4.sample_ising({0: .1, 1: .2}, {(0, 1): 1.5})
+
+
+class Test_ValidateChainStrength(unittest.TestCase):
+
+    def test_no_properties(self):
+        from dwave.system.composites.virtual_graph import _validate_chain_strength
+
+        Sampler = collections.namedtuple('Sampler', ['properties'])
+        sampler = Sampler({})
+
+        with self.assertRaises(ValueError):
+            _validate_chain_strength(sampler, None)
+
+        with self.assertRaises(ValueError):
+            _validate_chain_strength(sampler, 1.0)
+
+    def test_j_range(self):
+        from dwave.system.composites.virtual_graph import _validate_chain_strength
+
+        Sampler = collections.namedtuple('Sampler', ['properties'])
+        sampler = Sampler({'j_range': [-1.0, 1.0]})
+
+        self.assertEqual(_validate_chain_strength(sampler, None), 1.0)
+        self.assertEqual(_validate_chain_strength(sampler, .5), .5)
+
+        with self.assertRaises(ValueError):
+            _validate_chain_strength(sampler, 1.5)
+
+    def test_extended_j_range(self):
+        from dwave.system.composites.virtual_graph import _validate_chain_strength
+
+        Sampler = collections.namedtuple('Sampler', ['properties'])
+        sampler = Sampler({'j_range': [-1.0, 1.0], 'extended_j_range': [-2.0, 2.0]})
+
+        self.assertEqual(_validate_chain_strength(sampler, None), 2.0)
+        self.assertEqual(_validate_chain_strength(sampler, .5), .5)
+        self.assertEqual(_validate_chain_strength(sampler, 1.5), 1.5)