Fix bit ordering and probabilities of samples in optimization_algorithm.py (backport #97)

.pylintdict

@@ -104,6 +104,7 @@ qp
 quadratically
 quadraticconstraint
 quadraticprogram
+quantuminstance
 qubit
 qubits
 qubo

qiskit_optimization/algorithms/grover_optimizer.py

@@ -59,10 +59,10 @@ def __init__(self, num_value_qubits: int, num_iterations: int = 3,
             TypeError: When there one of converters is an invalid type.
         """
         self._num_value_qubits = num_value_qubits
-        self._num_key_qubits = None
+        self._num_key_qubits = 0
         self._n_iterations = num_iterations
-        self._quantum_instance = None
-        self._circuit_results = {}  # type: ignore
+        self._quantum_instance = None  # type: Optional[QuantumInstance]
+        self._circuit_results = {}  # type: dict
 
         if quantum_instance is not None:
             self.quantum_instance = quantum_instance
@@ -168,7 +168,7 @@ def solve(self, problem: QuadraticProgram) -> OptimizationResult:
                 problem_.objective.linear[i] = -val
             for (i, j), val in problem_.objective.quadratic.to_dict().items():
                 problem_.objective.quadratic[i, j] = -val
-        self._num_key_qubits = len(problem_.objective.linear.to_array())  # type: ignore
+        self._num_key_qubits = len(problem_.objective.linear.to_array())
 
         # Variables for tracking the optimum.
         optimum_found = False
@@ -237,13 +237,14 @@ def solve(self, problem: QuadraticProgram) -> OptimizationResult:
                     threshold = optimum_value
 
                     # trace out work qubits and store samples
-                    if self._quantum_instance.is_statevector:  # type: ignore
+                    if self._quantum_instance.is_statevector:
                         indices = list(range(n_key, len(outcome)))
                         rho = partial_trace(self._circuit_results, indices)
                         self._circuit_results = np.diag(rho.data) ** 0.5
                     else:
-                        self._circuit_results = {i[0:n_key]: v for i,
-                                                 v in self._circuit_results.items()}
+                        self._circuit_results = {
+                            i[-1 * n_key:]: v for i, v in self._circuit_results.items()
+                        }
 
                     raw_samples = self._eigenvector_to_solutions(self._circuit_results,
                                                                  problem_init)
@@ -313,7 +314,7 @@ def _get_probs(self, qc: QuantumCircuit) -> Dict[str, float]:
             state = result.get_counts(qc)
             shots = self.quantum_instance.run_config.shots
             hist = {key[::-1]: val / shots for key, val in state.items() if val > 0}
-            self._circuit_results = {b[::-1]: np.sqrt(v / shots) for (b, v) in state.items()}
+            self._circuit_results = {b: (v / shots) ** 0.5 for (b, v) in state.items()}
         return hist
 
     @staticmethod

qiskit_optimization/algorithms/optimization_algorithm.py

@@ -503,22 +503,23 @@ def _eigenvector_to_solutions(eigenvector: Union[dict, np.ndarray, StateFn],
             TypeError: If the type of eigenvector is not supported.
         """
         if isinstance(eigenvector, DictStateFn):
-            eigenvector = {bitstr: val ** 2 for (bitstr, val) in eigenvector.primitive.items()}
+            eigenvector = eigenvector.primitive
         elif isinstance(eigenvector, StateFn):
             eigenvector = eigenvector.to_matrix()
 
         def generate_solution(bitstr, qubo, probability):
-            x = np.fromiter(list(bitstr), dtype=int)
+            x = np.fromiter(list(bitstr[::-1]), dtype=int)
             fval = qubo.objective.evaluate(x)
             return SolutionSample(x=x, fval=fval, probability=probability,
                                   status=OptimizationResultStatus.SUCCESS)
 
         solutions = []
         if isinstance(eigenvector, dict):
-            all_counts = sum(eigenvector.values())
+            # When eigenvector is a dict, square the values since the values are normalized.
+            # See https://github.com/Qiskit/qiskit-terra/pull/5496 for more details.
+            probabilities = {bitstr: val ** 2 for (bitstr, val) in eigenvector.items()}
             # iterate over all samples
-            for bitstr, count in eigenvector.items():
-                sampling_probability = count / all_counts
+            for bitstr, sampling_probability in probabilities.items():
                 # add the bitstring, if the sampling probability exceeds the threshold
                 if sampling_probability >= min_probability:
                     solutions.append(generate_solution(bitstr, qubo, sampling_probability))
@@ -531,7 +532,7 @@ def generate_solution(bitstr, qubo, probability):
             for i, sampling_probability in enumerate(probabilities):
                 # add the i-th state if the sampling probability exceeds the threshold
                 if sampling_probability >= min_probability:
-                    bitstr = '{:b}'.format(i).rjust(num_qubits, '0')[::-1]
+                    bitstr = "{:b}".format(i).rjust(num_qubits, "0")
                     solutions.append(generate_solution(bitstr, qubo, sampling_probability))
 
         else:

qiskit_optimization/problems/quadratic_program.py

@@ -1153,7 +1153,7 @@ def substitute_variables(
         Returns:
             An optimization problem by substituting variables with constants or other variables.
             If the substitution is valid, `QuadraticProgram.status` is still
-            `QuadraticProgram.Status.VALIAD`.
+            `QuadraticProgram.Status.VALID`.
             Otherwise, it gets `QuadraticProgram.Status.INFEASIBLE`.
 
         Raises:
@@ -1166,6 +1166,10 @@ def substitute_variables(
     def to_ising(self) -> Tuple[OperatorBase, float]:
         """Return the Ising Hamiltonian of this problem.
 
+        Variables are mapped to qubits in the same order, i.e.,
+        i-th variable is mapped to i-th qubit.
+        See https://github.com/Qiskit/qiskit-terra/issues/1148 for details.
+
         Returns:
             qubit_op: The qubit operator for the problem
             offset: The constant value in the Ising Hamiltonian.
@@ -1265,6 +1269,10 @@ def from_ising(self,
                    offset: float = 0.0, linear: bool = False) -> None:
         r"""Create a quadratic program from a qubit operator and a shift value.
 
+        Variables are mapped to qubits in the same order, i.e.,
+        i-th variable is mapped to i-th qubit.
+        See https://github.com/Qiskit/qiskit-terra/issues/1148 for details.
+
         Args:
             qubit_op: The qubit operator of the problem.
             offset: The constant value in the Ising Hamiltonian.
@@ -1464,7 +1472,7 @@ def substitute_variables(
         Returns:
             An optimization problem by substituting variables with constants or other variables.
             If the substitution is valid, `QuadraticProgram.status` is still
-            `QuadraticProgram.Status.VALIAD`.
+            `QuadraticProgram.Status.VALID`.
             Otherwise, it gets `QuadraticProgram.Status.INFEASIBLE`.
 
         Raises:

releasenotes/notes/fix-bit-ordering-e807ec9f4b206ec3.yaml

@@ -0,0 +1,8 @@
+---
+fixes:
+  - |
+    Fix bit ordering in :class`qiskit_optimization.algorithms.minimum_eigen_optimizer.MinimumEigenOptimizer`
+    with qasm_simulator.
+  - |
+    Fix probabilities of solution samples with qasm_simulator.
+    See https://github.com/Qiskit/qiskit-optimization/pull/97 for details.

test/algorithms/test_grover_optimizer.py

@@ -48,7 +48,6 @@ def validate_results(self, problem, results):
         # Get expected value.
         solver = MinimumEigenOptimizer(NumPyMinimumEigensolver())
         comp_result = solver.solve(problem)
-
         # Validate results.
         np.testing.assert_array_almost_equal(comp_result.x, results.x)
         self.assertEqual(comp_result.fval, results.fval)
@@ -137,6 +136,7 @@ def test_qubo_gas_int_paper_example(self, simulator):
     def test_converter_list(self):
         """Test converters list"""
         # Input.
+
         model = Model()
         x_0 = model.binary_var(name='x0')
         x_1 = model.binary_var(name='x1')
@@ -152,6 +152,7 @@ def test_converter_list(self):
                               converters=qp2qubo)
         results = gmf.solve(op)
         self.validate_results(op, results)
+
         # a list of converters
         ineq2eq = InequalityToEquality()
         int2bin = IntegerToBinary()
@@ -168,31 +169,62 @@ def test_converter_list(self):
                             quantum_instance=self.sv_simulator,
                             converters=invalid)
 
-    def test_samples_and_raw_samples(self):
+    @data("sv", "qasm")
+    def test_samples_and_raw_samples(self, simulator):
         """Test samples and raw_samples"""
         op = QuadraticProgram()
-        op.integer_var(0, 3, 'x')
-        op.binary_var('y')
-        op.minimize(linear={'x': 1, 'y': 2})
-        op.linear_constraint(linear={'x': 1, 'y': 1}, sense='>=', rhs=1, name='xy')
+        op.integer_var(0, 3, "x")
+        op.binary_var("y")
+        op.minimize(linear={"x": 1, "y": 2})
+        op.linear_constraint(linear={"x": 1, "y": 1}, sense=">=", rhs=1, name="xy")
+        q_instance = self.sv_simulator if simulator == "sv" else self.qasm_simulator
+        grover_optimizer = GroverOptimizer(8, num_iterations=10, quantum_instance=q_instance)
         opt_sol = 1
         success = OptimizationResultStatus.SUCCESS
-        algorithm_globals.random_seed = 1
-        grover_optimizer = GroverOptimizer(
-            8, num_iterations=5, quantum_instance=self.qasm_simulator)
-        result = grover_optimizer.solve(op)
-        self.assertEqual(len(result.samples), 8)
-        self.assertEqual(len(result.raw_samples), 32)
-        self.assertAlmostEqual(sum(s.probability for s in result.samples), 1)
-        self.assertAlmostEqual(sum(s.probability for s in result.raw_samples), 1)
-        self.assertAlmostEqual(min(s.fval for s in result.samples), 0)
-        self.assertAlmostEqual(min(s.fval for s in result.samples if s.status == success), opt_sol)
-        self.assertAlmostEqual(min(s.fval for s in result.raw_samples), opt_sol)
-        for sample in result.raw_samples:
+        results = grover_optimizer.solve(op)
+        self.assertEqual(len(results.samples), 8)
+        self.assertEqual(len(results.raw_samples), 32)
+        self.assertAlmostEqual(sum(s.probability for s in results.samples), 1)
+        self.assertAlmostEqual(sum(s.probability for s in results.raw_samples), 1)
+        self.assertAlmostEqual(min(s.fval for s in results.samples), 0)
+        self.assertAlmostEqual(min(s.fval for s in results.samples if s.status == success), opt_sol)
+        self.assertAlmostEqual(min(s.fval for s in results.raw_samples), opt_sol)
+        for sample in results.raw_samples:
+            self.assertEqual(sample.status, success)
+        np.testing.assert_array_almost_equal(results.x, results.samples[0].x)
+        self.assertAlmostEqual(results.fval, results.samples[0].fval)
+        self.assertEqual(results.status, results.samples[0].status)
+        self.assertAlmostEqual(results.fval, results.raw_samples[0].fval)
+        self.assertEqual(results.status, results.raw_samples[0].status)
+        np.testing.assert_array_almost_equal([1, 0, 0, 0, 0], results.raw_samples[0].x)
+
+    @data("sv", "qasm")
+    def test_bit_ordering(self, simulator):
+        """Test bit ordering"""
+        # test minimize
+        algorithm_globals.random_seed = 2
+        q_instance = self.sv_simulator if simulator == "sv" else self.qasm_simulator
+        mdl = Model("docplex model")
+        x = mdl.binary_var("x")
+        y = mdl.binary_var("y")
+        mdl.minimize(x - 2 * y)
+        op = QuadraticProgram()
+        op.from_docplex(mdl)
+        opt_sol = -2
+        success = OptimizationResultStatus.SUCCESS
+        grover_optimizer = GroverOptimizer(3, num_iterations=10, quantum_instance=q_instance)
+        results = grover_optimizer.solve(op)
+        self.assertEqual(results.fval, opt_sol)
+        np.testing.assert_array_almost_equal(results.x, [0, 1])
+        self.assertEqual(results.status, success)
+        results.raw_samples.sort(key=lambda x: x.probability, reverse=True)
+        self.assertAlmostEqual(sum(s.probability for s in results.samples), 1, delta=1e-5)
+        self.assertAlmostEqual(sum(s.probability for s in results.raw_samples), 1, delta=1e-5)
+        self.assertAlmostEqual(min(s.fval for s in results.samples), -2)
+        self.assertAlmostEqual(min(s.fval for s in results.samples if s.status == success), opt_sol)
+        self.assertAlmostEqual(min(s.fval for s in results.raw_samples), opt_sol)
+        for sample in results.raw_samples:
             self.assertEqual(sample.status, success)
-        np.testing.assert_array_almost_equal(result.x, result.raw_samples[0].x[0:2])
-        self.assertAlmostEqual(result.fval, result.raw_samples[0].fval)
-        self.assertEqual(result.status, result.raw_samples[0].status)
 
 
 if __name__ == '__main__':