Repository: dwavesystems/dwave-system (language: python)

Issue:
Allow greedy substitute sampler in `MockDWaveSampler` to try harder
**Current Problem**
Some doctests (that use `MockDWaveSampler` implicitly) are failing (sporadically) because they were written with `TabuSampler` level of result quality in mind, yet we switched to greedy sampler in #443.

One example of a such failure is visible [here](https://app.circleci.com/pipelines/github/dwavesystems/dwave-system/340/workflows/4d665991-7fb7-422f-9270-bbe37f28e73a/jobs/6919), as part of #436.

**Proposed Solution**
Introduce "oversampling rate" parameter used by the greedy sampler in `MockDWaveSampler`. Oversampling would default to none, but could be tweaked by doctests init code. Mock sampler would take `num_reads` * `oversampling_rate` greedy samples, but would truncate to and return only the best `num_reads` samples.

Albeit this is only a probabilistic fix, setting oversampling rate to >10 is shown to regularly resolve all failures of current tests.

**Alternatives Considered**
1. Simplify problems used in (mock) tests and doctests to strictly convex, thus guaranteed to be solved by the greedy substitute sampler. This might hurt the pedagogical aspect of examples in docs.
2. Revert to `TabuSampler`. Goes against the rationale to switch to greedy in the first place.

/cc @jackraymond @JoelPasvolsky 


Discussion / hints:
An alternative possibility, that I like is for MockDWaveSampler() to swap in the ground state for the first sample when num_variables<=X, where X could be as large as 16. Determining the ground state by exhaustive enumeration is fast enough in this regime, and it would allow us to avoid probabilistic guarantees on our test up to a scale X. Any "test" bigger than X~10 that is relying on probabilistic guarantees is probably one that we should replace.  
@jackraymond, I initially didn't like the discontinuities this approach has -- only one ground state is inserted, only up to a certain problem size. But on second thought, the ground state inserted can conceivably be seen as generated by greedy descent starting from a "lucky" initial state. And it's great that's it's deterministic for the class of problems we should use in our (doc)tests anyway. So, let's go with your proposal instead! I'll open a PR.

Fix patch:
diff --git a/docs/conf.py b/docs/conf.py
index cc865c1f..4e18dae0 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -98,8 +98,6 @@
 modindex_common_prefix = ['dwave-system.']
 
 doctest_global_setup = """
-from __future__ import print_function, division
-
 import dimod
 from dwave.embedding import *
 
diff --git a/requirements.txt b/requirements.txt
index 5a56eafb..6105aa33 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,12 +1,12 @@
 --extra-index-url https://pypi.dwavesys.com/simple
 
-dimod==0.11.0
+dimod==0.11.6
 dwave-preprocessing==0.4.0
-dwave-cloud-client==0.10.0
-dwave-networkx==0.8.10
+dwave-cloud-client==0.10.2
+dwave-networkx==0.8.12
 dwave-drivers==0.4.4
-dwave-greedy
+dwave-greedy==0.2.5
 homebase==1.0.1
-minorminer==0.2.8.post1
-numpy==1.20.0; python_version < "3.10"
-numpy==1.21.4; python_version >= "3.10"
+minorminer==0.2.9
+numpy==1.21.6
+scipy==1.7.3    # fix py310/win curve_fit import error in 1.9.3
diff --git a/setup.py b/setup.py
index b924e28f..581f66e1 100644
--- a/setup.py
+++ b/setup.py
@@ -33,7 +33,7 @@
                     'homebase>=1.0.0,<2.0.0',
                     'minorminer>=0.2.8,<0.3.0',
                     'numpy>=1.17.3,<2.0.0',
-                    'dwave-tabu>=0.4.2',
+                    'dwave-greedy>=0.2.0',
                     ]
 
 # NOTE: dwave-drivers can also be installed with `dwave install drivers`,


Test patch:
diff --git a/dwave/system/testing.py b/dwave/system/testing.py
index f94b3736..d49adb1c 100644
--- a/dwave/system/testing.py
+++ b/dwave/system/testing.py
@@ -12,24 +12,21 @@
 #    See the License for the specific language governing permissions and
 #    limitations under the License.
 
+import concurrent.futures
 import warnings
-import collections
+from uuid import uuid4
+
+import numpy as np
 import dimod
 import dwave.cloud.computation
-import unittest.mock as mock
-import dwave_networkx as dnx
-
 from dwave.system import qpu_graph
-from uuid import uuid4
-
 
 try:
     #Simplest and efficient choice returning low energy states:
     from greedy import SteepestDescentSampler as SubstituteSampler
 except ImportError:
     from dimod import SimulatedAnnealingSampler as SubstituteSampler
-import concurrent.futures
-import numpy as np
+
 
 class MockDWaveSampler(dimod.Sampler, dimod.Structured):
     """Mock sampler modeled after DWaveSampler that can be used for tests.
@@ -56,7 +53,7 @@ class MockDWaveSampler(dimod.Sampler, dimod.Structured):
             List of nodes to exclude (along with associated edges), for
             emulation of unyielded qubits. This parameter is made redundant
             by the use of nodelist.
-        
+
         broken_edges (iterable of (int,int) tuples, optional):
             List of edges to exclude, for emulation of unyielded edges.
             This parameter is made redundant by the use of edgelist.
@@ -65,14 +62,14 @@ class MockDWaveSampler(dimod.Sampler, dimod.Structured):
             Supported options are 'chimera', 'pegasus' or 'zephyr'.
             The default is 'chimera' when the value is not specified as part 
             of the from the ``properties`` argument.
-            
+
         topology_shape (string, optional):
             A list of three numbers [m,n,t] for Chimera, defaulted as [4,4,4]. 
             A list of one number [m] for Pegasus, defaulted as [3].
             A list of two numbers [m,t] for Zephyr, defaulted as [2,4].
             The default above apply only when the value is not 
             specified as part of the from the ``properties`` argument.
-            
+
         parameter_warnings (bool, optional, default=True):
             The MockSampler is adaptive with respect to ``num_reads``,
             ``answer_mode`` and ``max_answers`` and ``label`` 
@@ -80,6 +77,16 @@ class MockDWaveSampler(dimod.Sampler, dimod.Structured):
             dwave-greedy is installed. All other parameters are ignored and a 
             warning will be raised 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 cut offs
+            should be used since solution time increases exponentially with
+            problem size.
+            Set ``exact_solver_cutoff`` to zero to disable exact ground state
+            calculation.
+
     Examples
         The first example creates a MockSampler without reference to a
         particular online system, with a fully-yielded Chimera C16 structure.
@@ -105,17 +112,24 @@ class MockDWaveSampler(dimod.Sampler, dimod.Structured):
 
     """
     # 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, **config):
+                 topology_type=None, topology_shape=None,
+                 parameter_warnings=True,
+                 exact_solver_cutoff=EXACT_SOLVER_CUTOFF_DEFAULT,
+                 **config):
         self.parameter_warnings = parameter_warnings
+        self.exact_solver_cutoff = exact_solver_cutoff
 
         #Parse or default topology dependent arguments:
         if properties is not None and 'topology' in properties:
@@ -134,9 +148,9 @@ def __init__(self,
         else:
             if topology_type is None:
                 topology_type = 'chimera'
-            shape_defaults = {'chimera' : [4,4,4],
-                              'pegasus' : [3],
-                              'zephyr' : [2,4]}
+            shape_defaults = {'chimera': [4,4,4],
+                              'pegasus': [3],
+                              'zephyr': [2,4]}
             if topology_type in shape_defaults:
                 if topology_shape is None:
                     topology_shape = shape_defaults[topology_type]
@@ -144,9 +158,10 @@ def __init__(self,
                 raise ValueError("Only 'chimera', 'pegasus' and 'zephyr' "
                                  "topologies are supported")
         self.properties = {
-            'chip_id' : 'MockDWaveSampler',
-            'topology' : {'type': topology_type, 'shape': topology_shape}}
-        
+            '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()
@@ -158,7 +173,7 @@ def __init__(self,
         solver_graph = qpu_graph(self.properties['topology']['type'],
                                  self.properties['topology']['shape'],
                                  self.nodelist, self.edgelist)
-        
+
         if broken_nodes is None and broken_edges is None:
             self.nodelist = sorted(solver_graph.nodes)
             self.edgelist = sorted(tuple(sorted(edge))
@@ -167,7 +182,7 @@ def __init__(self,
             if broken_nodes is None:
                 broken_nodes = []
             self.nodelist = sorted(set(solver_graph.nodes).difference(broken_nodes))
-            if broken_edges == None:
+            if broken_edges is None:
                 broken_edges = []
             self.edgelist = sorted(tuple(sorted((u, v))) for
                                    u, v in solver_graph.edges
@@ -177,20 +192,19 @@ def __init__(self,
                                    and (v, u) not in broken_edges)
         #Finalize yield-dependent properties:
         self.properties.update({
-            '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
-                                      else [0,0] for i in range(len(self.nodelist))]})
+            '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
+                                     else [0, 0] for i in range(len(self.nodelist))]})
         # Non-topology-dependent properties and parameters mocked from
         # Advantage_system4.1 accessed February 10th 2022
         # with simplified lists for large parameters, and modified
         # topology arguments per MockSolver initialization:
         # See also:
         # https://docs.dwavesys.com/docs/latest/c_solver_parameters.html
-        
+
         self.parameters = {
-            
             'anneal_offsets': ['parameters'],
             'anneal_schedule': ['parameters'],
             'annealing_time': ['parameters'],
@@ -208,9 +222,9 @@ def __init__(self,
             'reduce_intersample_correlation': ['parameters'],
             'reinitialize_state': ['parameters'],
             'warnings': [],
-            'label': []}
-        
-        
+            'label': [],
+        }
+
         self.properties.update({
             'h_range': [-4.0, 4.0],
             'j_range': [-1.0, 1.0],
@@ -283,21 +297,24 @@ def __init__(self,
             'readout_thermalization_range': [0.0, 10000.0],
             'tags': [],
             'category': 'qpu',
-            'quota_conversion_rate': 1})
-        
+            'quota_conversion_rate': 1,
+        })
+
         if properties is not None:
             # provided properties overwrite defaults, note that
             # particular care should be taken with respect to
             # topology-dependent arguments:
             self.properties.update(properties)
-        
+
     @classmethod
     def from_qpu_sampler(cls, sampler):
-        return cls(properties=sampler.properties,nodelist=sampler.nodelist,edgelist=sampler.edgelist)
+        return cls(properties=sampler.properties,
+                   nodelist=sampler.nodelist,
+                   edgelist=sampler.edgelist)
 
     @dimod.bqm_structured
     def sample(self, bqm, **kwargs):
-        
+
         # Check kwargs compatibility with parameters and substitute sampler:
         mocked_parameters={'answer_mode',
                            'max_answers',
@@ -315,7 +332,7 @@ def sample(self, bqm, **kwargs):
             # (2) We could also consider using 'large_sparse_opt' to
             # exploit fixed connectivity at large scale for current QPU
             # designs, but unlikely to be a significant inefficiency.
-            
+
             # The fallback sampler (SA) is an inferior choice, but
             # greedy may not always be installed, unlike dimod.
         for kw in kwargs:
@@ -349,13 +366,12 @@ def sample(self, bqm, **kwargs):
         flux_biases = kwargs.get('flux_biases')
         if flux_biases is not None:
             self.flux_biases_flag = True
-            
-            
+
         substitute_kwargs = {'num_reads' : kwargs.get('num_reads')}
         if substitute_kwargs['num_reads'] is None:
             substitute_kwargs['num_reads'] = 1
-        
-        if SubstituteSampler is not dimod.SimulatedAnnealingSampler:            
+
+        if SubstituteSampler is not dimod.SimulatedAnnealingSampler:
             initial_state = kwargs.get('initial_state')
             if initial_state is not None:
                 # Initial state format is a list of (qubit,values)
@@ -366,20 +382,26 @@ def sample(self, bqm, **kwargs):
                     np.array([pair[1] for pair in initial_state
                               if pair[1]!=3],dtype=float),
                     [pair[0] for pair in initial_state if pair[1]!=3])
-                
+
         ss = SubstituteSampler().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'
             ss = ss.aggregate()
-        
+
         max_answers = kwargs.get('max_answers')
         if max_answers is not None:
             # Truncate sampleset if requested. Do not reorder (per DWaveSampler())
             ss = ss.truncate(max_answers)
-        
+
         return ss
 
     def to_networkx_graph(self):
@@ -424,7 +446,8 @@ def sample_dqm(self, dqm, **kwargs):
     def min_time_limit(self, dqm):
         # not caring about the problem, just returning the min
         return self.properties['minimum_time_limit'][0][1]
-    
+
+
 class MockLeapHybridSolver:
 
     properties = {'supported_problem_types': ['bqm'],
diff --git a/tests/test_embedding_composite.py b/tests/test_embedding_composite.py
index 4be9bc43..d5a2d4fb 100644
--- a/tests/test_embedding_composite.py
+++ b/tests/test_embedding_composite.py
@@ -17,6 +17,7 @@
 import warnings
 
 from collections.abc import Mapping
+from unittest import mock
 
 import dimod
 import dwave_networkx as dnx
@@ -31,7 +32,7 @@
                                      AutoEmbeddingComposite,
                                      )
 
-from dwave.system.testing import MockDWaveSampler, mock
+from dwave.system.testing import MockDWaveSampler
 from dwave.embedding import chain_breaks
 from dwave.system.warnings import ChainStrengthWarning
 
diff --git a/tests/test_mock_sampler.py b/tests/test_mock_sampler.py
index c5665bc3..cea4049a 100644
--- a/tests/test_mock_sampler.py
+++ b/tests/test_mock_sampler.py
@@ -13,8 +13,12 @@
 #    limitations under the License.
 
 """who watches the watchmen?"""
-import unittest
 import os
+import unittest
+from unittest import mock
+
+import numpy as np
+import dimod
 import dimod.testing as dit
 
 from dwave.system import DWaveSampler
@@ -22,6 +26,7 @@
 from dwave.cloud.exceptions import ConfigFileError, SolverNotFoundError
 from dimod import DiscreteQuadraticModel, ExtendedVartype, SampleSet
 
+
 class TestMockDWaveSampler(unittest.TestCase):
     @unittest.skipIf(os.getenv('SKIP_INT_TESTS'), "Skipping integration test.")
     def test_properties_and_params(self):
@@ -39,48 +44,40 @@ def test_properties_and_params(self):
         mock = MockDWaveSampler.from_qpu_sampler(sampler)
         self.assertEqual(mock.nodelist.sort(),sampler.nodelist.sort())
         self.assertEqual(mock.edgelist.sort(),sampler.edgelist.sort())
-        
-        
+
     def test_sampler(self):
         sampler = MockDWaveSampler()
         dit.assert_sampler_api(sampler)
         dit.assert_structured_api(sampler)
-        
+
     def test_mock_parameters(self):
         # Single cell Chimera (DW2000Q) solver:
         sampler = MockDWaveSampler(topology_type='chimera', topology_shape = [1,1,4])
         num_reads = 10
-        ss = sampler.sample_ising({0 : -1}, {}, num_reads=num_reads)
+        ss = sampler.sample_ising({0: -1}, {}, num_reads=num_reads)
         self.assertEqual(sum(ss.record.num_occurrences),num_reads)
-        self.assertTrue(len(ss.record.num_occurrences)<=2) #Binary state histogram
-        ss = sampler.sample_ising({0 : -1}, {}, num_reads=num_reads,
+        self.assertTrue(len(ss.record.num_occurrences) <= 2) #Binary state histogram
+        ss = sampler.sample_ising({0: -1}, {}, num_reads=num_reads,
                                   answer_mode='raw')
-        self.assertEqual(len(ss.record.num_occurrences),num_reads)
+        self.assertEqual(len(ss.record.num_occurrences), num_reads)
         
-        ss = sampler.sample_ising({0 : -1}, {}, num_reads=num_reads,
+        ss = sampler.sample_ising({0: -1}, {}, num_reads=num_reads,
                                   answer_mode='raw', max_answers=2)
-        self.assertEqual(len(ss.record.num_occurrences),2)
-
-            
-        try:
-            from greedy import SteepestDescentSampler as SubstituteSampler
-            mock_fallback_substitute = False
-        except:
-            mock_fallback_substitute = True
-        if not mock_fallback_substitute:
-            # If the greedy is available, we can mock more parameters (than
-            # is possible with the fallback dimod.SimulatedAnnealingSampler()
-            # method
+        self.assertEqual(len(ss.record.num_occurrences), 2)
 
+        # note: greedy should be available; it's a package/test requirement
+        # disable exact ground state calc
+        with mock.patch.object(sampler, 'exact_solver_cutoff', 0):
             #QPU format initial states:
             initial_state = [(i,1) if i%4==0 else (i,3) for i in range(8)]
-            ss = sampler.sample_ising({0 : 1, 4 : 1}, {(0,4) : -2},
+            ss = sampler.sample_ising({0: 1, 4: 1}, {(0, 4): -2},
                                       num_reads=1,
                                       answer_mode='raw',
                                       initial_state=initial_state)
             #The initialized state is a local minima, and should
             #be returned under greedy descent mocking:
             self.assertEqual(ss.record.energy[0],0)
+
     def test_unmock_parameters(self):
         # Check valid DWaveSampler() parameter, that is not mocked throws a warning:
         sampler = MockDWaveSampler()
@@ -95,7 +92,39 @@ def test_unmock_parameters(self):
         with warnings.catch_warnings():
             warnings.simplefilter(action="error", category=UserWarning)
             ss = sampler.sample_ising({0 : -1}, {}, annealing_time=123)
-    
+
+    def test_ground_override(self):
+        # bqm with a local minimum at (-1, -1)
+        bqm = dimod.BQM.from_ising({'a': -1, 'b': -1}, {'ab': -2})
+        local_minimum = [-1, -1]    # energy: 0
+        ground_state = [1, 1]       # energy: -4
+        local_minimum_state = list(zip(bqm.variables, local_minimum))
+
+        # simple two connected variables sampler (path(2))
+        sampler = MockDWaveSampler(nodelist=['a', 'b'], edgelist=[('a', 'b')],
+                                   exact_solver_cutoff=len(bqm.variables))
+
+        # disable exact ground state calc and start from a local minimum -> greedy should stay stuck
+        with mock.patch.object(sampler, 'exact_solver_cutoff', 0):
+            ss = sampler.sample(bqm, initial_state=local_minimum_state)
+            np.testing.assert_array_equal(ss.record[0].sample, local_minimum)
+
+        # boundary on which exact ground state calc kicks in
+        with mock.patch.object(sampler, 'exact_solver_cutoff', len(bqm.variables)):
+            ss = sampler.sample(bqm, initial_state=local_minimum_state)
+            np.testing.assert_array_equal(ss.record[0].sample, ground_state)
+
+        # double-check the default
+        self.assertEqual(sampler.exact_solver_cutoff, len(bqm.variables))
+        ss = sampler.sample(bqm)
+        np.testing.assert_array_equal(ss.record[0].sample, ground_state)
+
+    def test_empty_bqm(self):
+        sampler = MockDWaveSampler()
+        bqm = dimod.BQM('SPIN')
+        ss = sampler.sample(bqm)
+        self.assertIs(ss.vartype, bqm.vartype)
+
     def test_chimera_topology(self):
         grid_parameter = 5
         tile_parameter = 2