By default, use apply_unitary_effect_to_state in VariationalStudy eva…

…luations (#269)

* add LiH test

* add h2 test

* fix test function

* edit tests

* implement spatial low rank

* lint

* define LetterWithSubscripts and change ansatz params method

* fix tests and update notebook

* add LiH test

* add h2 test

* fix test function

* edit tests

* implement spatial low rank

* lint

* polish low rank trotter step

* define low rank ansatz

* add default initial params test

* add final swap gates

* add lih test

* fix factor of 2

* fix low rank

* port low rank changes

* remove non-unitary case

* update trotter algorithm

* start on ansatz

* update ansatz

* lint

* increase test tol: lih hamiltonian has off-diagonal one-body

* unused variable

* update low rank test and doc

* clean up split operator test

* remove debug statement

* allow objective circuit output to be numpy array

* more flexibility in VariationalBlackBox

* hide black_box_type and document it

* switch default black box to use unitary simulation

* use apply_unitary_effect_to_state in default initial params test

* add qubit permutation to low rank ansatz

* lint

* rename DEFAULT to UNITARY_SIMULATE

openfermioncirq/variational/ansatzes/default_initial_params_test.py

@@ -81,15 +81,12 @@ def test_trotter_ansatzes_default_initial_params_iterations_1(
                 occupied_orbitals=range(len(qubits) // 2))
     )
 
-    simulator = cirq.google.XmonSimulator()
-
     # Compute value using ansatz circuit and objective
-    result = simulator.simulate(
-            preparation_circuit + ansatz.circuit,
-            param_resolver=
-                ansatz.param_resolver(ansatz.default_initial_params()),
-            qubit_order=ansatz.qubit_permutation(qubits)
-    )
+    circuit = (preparation_circuit + ansatz.circuit
+              ).with_parameters_resolved_by(
+                      ansatz.param_resolver(ansatz.default_initial_params()))
+    result = circuit.apply_unitary_effect_to_state(
+            qubit_order=ansatz.qubit_permutation(qubits))
     obj_val = objective.value(result)
 
     # Compute value using study
@@ -120,8 +117,8 @@ def test_trotter_ansatzes_default_initial_params_iterations_1(
                 order=order,
                 algorithm=trotter_algorithm)
     )
-    result = simulator.simulate(preparation_circuit + simulation_circuit)
-    final_state = result.final_state
+    final_state = (preparation_circuit + simulation_circuit
+                  ).apply_unitary_effect_to_state()
     correct_val = openfermion.expectation(
             objective._hamiltonian_linear_op, final_state).real
 
@@ -165,15 +162,12 @@ def test_trotter_ansatzes_default_initial_params_iterations_2(
                 occupied_orbitals=range(len(qubits) // 2))
     )
 
-    simulator = cirq.google.XmonSimulator()
-
     # Compute value using ansatz circuit and objective
-    result = simulator.simulate(
-            preparation_circuit + ansatz.circuit,
-            param_resolver=
-                ansatz.param_resolver(ansatz.default_initial_params()),
-            qubit_order=ansatz.qubit_permutation(qubits)
-    )
+    circuit = (preparation_circuit + ansatz.circuit
+              ).with_parameters_resolved_by(
+                      ansatz.param_resolver(ansatz.default_initial_params()))
+    result = circuit.apply_unitary_effect_to_state(
+            qubit_order=ansatz.qubit_permutation(qubits))
     obj_val = objective.value(result)
 
     # Compute value using study
@@ -218,8 +212,8 @@ def test_trotter_ansatzes_default_initial_params_iterations_2(
                 order=order,
                 algorithm=trotter_algorithm)
     )
-    result = simulator.simulate(preparation_circuit + simulation_circuit)
-    final_state = result.final_state
+    final_state = (preparation_circuit + simulation_circuit
+                  ).apply_unitary_effect_to_state()
     correct_val = openfermion.expectation(
             objective._hamiltonian_linear_op, final_state).real
 

openfermioncirq/variational/ansatzes/low_rank.py

@@ -104,7 +104,6 @@ def __init__(self,
                  include_all_z: bool=False,
                  adiabatic_evolution_time: Optional[float]=None,
                  spin_basis: bool=True,
-                 omit_final_swaps: bool=False,
                  qubits: Optional[Sequence[cirq.QubitId]]=None
                  ) -> None:
         """
@@ -128,8 +127,6 @@ def __init__(self,
                 of the entries of the two-body tensor of the Hamiltonian.
             spin_basis: Whether the Hamiltonian is given in the spin orbital
                 (rather than spatial orbital) basis.
-            omit_final_swaps: If this is set to True, then SWAP gates at
-                the end of the circuit may be omitted.
             qubits: Qubits to be used by the ansatz circuit. If not specified,
                 then qubits will automatically be generated by the
                 `_generate_qubits` method.
@@ -139,7 +136,6 @@ def __init__(self,
         self.final_rank = final_rank
         self.include_all_cz = include_all_cz
         self.include_all_z = include_all_z
-        self.omit_final_swaps = omit_final_swaps
 
         if adiabatic_evolution_time is None:
             adiabatic_evolution_time = (
@@ -262,11 +258,14 @@ def two_body_interaction(p, q, a, b) -> cirq.OP_TREE:
             # Undo final basis transformation.
             yield bogoliubov_transform(qubits, prior_basis_matrix)
 
-        if not self.omit_final_swaps:
-            # If the number of swap networks was odd, swap the qubits back
-            if self.iterations & 1 and len(self.eigenvalues) & 1:
-                yield swap_network(qubits)
-
+    def qubit_permutation(self, qubits: Sequence[cirq.QubitId]
+                          ) -> Sequence[cirq.QubitId]:
+        """The qubit permutation induced by the ansatz circuit."""
+        # An odd number of swap networks reverses the qubit ordering
+        if self.iterations & 1 and len(self.eigenvalues) & 1:
+            return qubits[::-1]
+        else:
+            return qubits
 
     def default_initial_params(self) -> numpy.ndarray:
         """Approximate evolution by H(t) = T + (t/A)V.

openfermioncirq/variational/ansatzes/low_rank_test.py

@@ -54,8 +54,7 @@ def test_low_rank_trotter_ansatz_param_bounds():
 
 def test_low_rank_trotter_ansatz_circuit():
 
-    ansatz = LowRankTrotterAnsatz(
-            lih_hamiltonian, final_rank=2, omit_final_swaps=True)
+    ansatz = LowRankTrotterAnsatz(lih_hamiltonian, final_rank=2)
     circuit = ansatz.circuit
     cirq.DropNegligible().optimize_circuit(circuit)
     assert circuit.to_text_diagram(transpose=True).strip() == """

openfermioncirq/variational/study.py

@@ -77,7 +77,7 @@ def __init__(self,
                  target: Optional[float]=None,
                  black_box_type: Type[
                      variational_black_box.VariationalBlackBox]=
-                     variational_black_box.DEFAULT,
+                     variational_black_box.UNITARY_SIMULATE,
                  datadir: Optional[str]=None) -> None:
         """
         Args:

openfermioncirq/variational/study_test.py

@@ -26,7 +26,7 @@
 from openfermioncirq.variational.study import (
         VariationalStudy)
 from openfermioncirq.variational.variational_black_box import (
-        DefaultVariationalBlackBox)
+        XmonSimulateVariationalBlackBox)
 from openfermioncirq.testing import (
         ExampleAlgorithm,
         ExampleAnsatz,
@@ -41,14 +41,18 @@
 
 a, b = test_ansatz.qubits
 preparation_circuit = cirq.Circuit.from_ops(cirq.X(a))
-test_study = VariationalStudy('test_study',
-                              test_ansatz,
-                              test_objective,
-                              preparation_circuit=preparation_circuit)
-test_study_noisy = VariationalStudy('test_study_noisy',
-                                    test_ansatz,
-                                    test_objective_noisy,
-                                    preparation_circuit=preparation_circuit)
+test_study = VariationalStudy(
+        'test_study',
+        test_ansatz,
+        test_objective,
+        preparation_circuit=preparation_circuit,
+        black_box_type=variational_black_box.XMON_SIMULATE)
+test_study_noisy = VariationalStudy(
+        'test_study_noisy',
+        test_ansatz,
+        test_objective_noisy,
+        preparation_circuit=preparation_circuit,
+        black_box_type=variational_black_box.XMON_SIMULATE)
 
 
 def test_variational_study_circuit():
@@ -68,7 +72,7 @@ def test_variational_study_optimize_and_extend_and_summary():
 
     study = VariationalStudy(
             'study', test_ansatz, test_objective,
-            black_box_type=variational_black_box.DEFAULT_STATEFUL,
+            black_box_type=variational_black_box.XMON_SIMULATE_STATEFUL,
             target=-10.5)
     assert len(study.trial_results) == 0
     assert study.target == -10.5
@@ -111,7 +115,7 @@ def test_variational_study_optimize_and_extend_and_summary():
     assert result.repetitions == 1
     assert all(
             result.data_frame['optimal_parameters'].apply(
-                lambda x: DefaultVariationalBlackBox(
+                lambda x: XmonSimulateVariationalBlackBox(
                     test_ansatz, test_objective).evaluate(x))
             == result.data_frame['optimal_value'])
     assert isinstance(result.results[0].cost_spent, float)
@@ -140,7 +144,8 @@ def test_variational_study_save_load():
             study_name,
             test_ansatz,
             test_objective,
-            datadir=datadir)
+            datadir=datadir,
+            black_box_type=variational_black_box.XMON_SIMULATE_STATEFUL)
     study.optimize(
             OptimizationParams(
                 ScipyOptimizationAlgorithm(
@@ -179,36 +184,36 @@ def test_variational_study_save_load():
 
 
 def test_variational_black_box_dimension():
-    black_box = DefaultVariationalBlackBox(test_ansatz, test_objective)
+    black_box = XmonSimulateVariationalBlackBox(test_ansatz, test_objective)
     assert black_box.dimension == 2
 
 
 def test_variational_black_box_bounds():
-    black_box = DefaultVariationalBlackBox(test_ansatz, test_objective)
+    black_box = XmonSimulateVariationalBlackBox(test_ansatz, test_objective)
     assert black_box.bounds == test_study.ansatz.param_bounds()
 
 
 def test_variational_black_box_noise_bounds():
-    black_box = DefaultVariationalBlackBox(test_ansatz, test_objective)
+    black_box = XmonSimulateVariationalBlackBox(test_ansatz, test_objective)
     assert black_box.noise_bounds(100) == (-numpy.inf, numpy.inf)
 
 
 def test_variational_black_box_evaluate():
-    black_box = DefaultVariationalBlackBox(test_ansatz, test_objective)
+    black_box = XmonSimulateVariationalBlackBox(test_ansatz, test_objective)
     numpy.testing.assert_allclose(
             black_box.evaluate(test_ansatz.default_initial_params()), 0.0)
     numpy.testing.assert_allclose(
             black_box.evaluate(numpy.array([0.0, 0.5])), 1.0)
 
 
 def test_variational_black_box_evaluate_with_cost():
-    black_box = DefaultVariationalBlackBox(test_ansatz, test_objective)
+    black_box = XmonSimulateVariationalBlackBox(test_ansatz, test_objective)
     numpy.testing.assert_allclose(
             black_box.evaluate_with_cost(
                 test_ansatz.default_initial_params(), 2.0),
             0.0)
 
-    black_box_noisy = DefaultVariationalBlackBox(
+    black_box_noisy = XmonSimulateVariationalBlackBox(
             test_ansatz, test_objective_noisy)
     numpy.random.seed(33534)
     noisy_val = black_box_noisy.evaluate_with_cost(

openfermioncirq/variational/variational_black_box.py

@@ -83,7 +83,27 @@ def noise_bounds(self,
         return self.objective.noise_bounds(cost, confidence)
 
 
-class DefaultVariationalBlackBox(VariationalBlackBox):
+class UnitarySimulateVariationalBlackBox(VariationalBlackBox):
+
+    def evaluate_noiseless(self,
+                           x: numpy.ndarray) -> float:
+        """Evaluate parameters with a noiseless simulation."""
+        # Default: evaluate using apply_unitary_effect_to_state
+        circuit = (self.preparation_circuit + self.ansatz.circuit
+                  ).with_parameters_resolved_by(self.ansatz.param_resolver(x))
+        final_state = circuit.apply_unitary_effect_to_state(
+                qubit_order=self.ansatz.qubit_permutation(self.ansatz.qubits))
+        return self.objective.value(final_state)
+
+
+class UnitarySimulateVariationalStatefulBlackBox(
+        UnitarySimulateVariationalBlackBox,
+        StatefulBlackBox):
+    """A stateful black box encapsulating a variational objective function."""
+    pass
+
+
+class XmonSimulateVariationalBlackBox(VariationalBlackBox):
 
     def evaluate_noiseless(self,
                            x: numpy.ndarray) -> float:
@@ -97,11 +117,13 @@ def evaluate_noiseless(self,
         return self.objective.value(result)
 
 
-class DefaultVariationalStatefulBlackBox(DefaultVariationalBlackBox,
-                                         StatefulBlackBox):
+class XmonSimulateVariationalStatefulBlackBox(XmonSimulateVariationalBlackBox,
+                                              StatefulBlackBox):
     """A stateful black box encapsulating a variational objective function."""
     pass
 
 
-DEFAULT = DefaultVariationalBlackBox
-DEFAULT_STATEFUL = DefaultVariationalStatefulBlackBox
+UNITARY_SIMULATE = UnitarySimulateVariationalBlackBox
+UNITARY_SIMULATE_STATEFUL = UnitarySimulateVariationalStatefulBlackBox
+XMON_SIMULATE = XmonSimulateVariationalBlackBox
+XMON_SIMULATE_STATEFUL = XmonSimulateVariationalStatefulBlackBox