Merge pull request #29 from qiskit-community/mixture_collapse

Add mixture collapse option and fix C++ instability

Author: haggaila

lindbladmpo/LindbladMPOSolver.py

@@ -928,6 +928,7 @@ def verify_parameters(
                 or (key == "b_save_final_state")
                 or (key == "b_initial_rho_compression")
                 or (key == "b_apply_gate_compression")
+                or (key == "b_collapse_equal_mixture")
             ):
                 if not isinstance(parameters[key], bool):
                     check_msg += (

lindbladmpo/examples/simulation_building/LindbladMatrixSolver.py

@@ -572,6 +572,7 @@ def verify_parameters(
             "trotter_order",
             "b_initial_rho_compression",
             "b_apply_gate_compression",
+            "b_collapse_equal_mixture",
             "max_dim_rho",
             "cut_off",
             "cut_off_rho",

src/ModelParameters.h

@@ -80,6 +80,7 @@ class ModelParameters : public SimulationParameters
 
         operator[]("custom_observables") = "";
         operator[]("collapse") = "";
+        operator[]("b_collapse_equal_mixture") = "";
     }
 };
 

src/lindbladian.h

@@ -210,7 +210,7 @@ bool SetLindbladian(SpinHalfSystem &C, ModelParameters param, Lattice2d L)
         g_2 = vector<double>(N, g_2[0]);
     if (g_3_len == 1)
         g_3 = vector<double>(N, g_3[0]);
-    if (g_2_len == 1)
+    if (g_4_len == 1)
         g_4 = vector<double>(N, g_4[0]);
 
     for (int i = 1; i <= int(N); i++)

src/lindbladmpo.cc

@@ -29,7 +29,7 @@ using namespace std::chrono;
 
 stream2d cout2 = stream2d(&cerr, NULL);
 
-const string SOLVER_VERSION = "0.3.1";
+const string SOLVER_VERSION = "0.3.2";
 
 const double IMAGINARY_THRESHOLD = 1e-4;
 // Threshold for the imaginary value of a quantity that should be real, to issue a warning
@@ -870,10 +870,13 @@ int main(int argc, char *argv[])
           << "\n";
     cout2.flush();
 
+    const bool b_collapse_equal_mixture = param.boolval("b_collapse_equal_mixture");
     if (collapse.size())
     {
         auto t_collapse_start = steady_clock::now();
         cout2 << "\nStarting evaluation of " << collapse.size() << " collapse operators.\n";
+        if (b_collapse_equal_mixture)
+            cout2 << "Creating an equal mixture of collapsed states.\n";
 
         int i_coll = 0;
         MPS rho_0(C.rho); // Keep a copy of the uncollapsed state
@@ -902,6 +905,14 @@ int main(int argc, char *argv[])
                         exit(1);
                 i_op++;
             }
+            if (b_collapse_equal_mixture)
+            {
+                Cplx z = C.trace_rho();
+                if (std::abs(z) < _2_N)
+                    cout2 << "\t\tNote: " << "Tr{rho} = " << z << " encountered during collapse projectors, "
+                          << "this is smaller than 2^(-N)!";
+                C.rho /= z;
+            }
             if (i_coll == 0)
                 rho_c = MPS(C.rho);
             else

test/model/test_model.py

@@ -628,30 +628,49 @@ def test_apply_gates_2(self):
 
     def test_apply_gates_cx_sx_h(self):
         """Compare the application of gates at intermediate times."""
-        N = 8
+        N = 13
         t_final = 1
         solver_params1 = {
-            "tau": 0.02,
+            "tau": 0.1,
             "t_final": t_final,
             "N": N,
-            "h_z": 0.43,
+            "h_z": 0.0,
             "1q_components": ["x", "y", "z"],
             "2q_components": ["xx", "yy", "zz", "xy", "xz", "yz"],
             "b_quiet": True,
-            "init_product_state": ["+x", "+x", "+z", "+z", "+z", "+z", "-z", "-z"],
+            "init_product_state": [
+                "+x",
+                "+x",
+                "+z",
+                "+z",
+                "+z",
+                "+z",
+                "-z",
+                "-z",
+                "+z",
+                "+z",
+                "+z",
+                "+z",
+                "+z",
+            ],
             "apply_gates": [
-                (0.0, "cz", 0, 1),
-                (0.0, "h", 2),
-                (0.0, "h", 3),
-                (0.0, "cz", 2, 3),
-                (0.0, "h", 4),
-                (0.0, "h", 5),
-                (0.0, "h", 5),
-                (0.0, "cx", 4, 5),
-                (0.0, "h", 5),
-                (0.0, "sx", 6),
-                (0.0, "sx", 7),
-                (0.0, "cz", 6, 7),
+                (0.1, "cz", 0, 1),
+                (0.2, "h", 2),
+                (0.3, "h", 3),
+                (0.4, "cz", 2, 3),
+                (0.5, "h", 4),
+                (0.6, "h", 5),
+                (0.7, "h", 5),
+                (0.8, "cx", 4, 5),
+                (0.9, "h", 5),
+                (0.3, "sx", 6),
+                (0.4, "sx", 7),
+                (0.5, "cz", 6, 7),
+                (0.7, "h", 10),
+                (0.7, "cx", 10, 8),
+                (0.7, "h", 12),
+                (0.8, "cx", 12, 11),
+                (0.9, "cx", 8, 11),
             ],
         }
         s_files_prefix1 = s_output_path + "test_apply_gates_cx_sx_h"
@@ -668,8 +687,8 @@ def test_apply_gates_cx_sx_h(self):
             solver1.result["obs-2q"][("xz", (2, 3))][1][-1],
         )
         self.assertAlmostEqual(
-            solver1.result["obs-2q"][("xx", (0, 1))][1][-1],
-            solver1.result["obs-2q"][("xx", (2, 3))][1][-1],
+            solver1.result["obs-2q"][("yy", (0, 1))][1][-1],
+            solver1.result["obs-2q"][("yy", (2, 3))][1][-1],
         )
         self.assertAlmostEqual(
             solver1.result["obs-2q"][("xy", (0, 1))][1][-1],
@@ -680,13 +699,17 @@ def test_apply_gates_cx_sx_h(self):
             solver1.result["obs-2q"][("xz", (4, 5))][1][-1],
         )
         self.assertAlmostEqual(
-            solver1.result["obs-2q"][("xx", (0, 1))][1][-1],
-            solver1.result["obs-2q"][("xx", (4, 5))][1][-1],
+            solver1.result["obs-2q"][("yy", (0, 1))][1][-1],
+            solver1.result["obs-2q"][("yy", (4, 5))][1][-1],
         )
         self.assertAlmostEqual(
             solver1.result["obs-2q"][("xy", (0, 1))][1][-1],
             solver1.result["obs-2q"][("xy", (4, 5))][1][-1],
         )
+        self.assertAlmostEqual(
+            solver1.result["global"][("s_2", ())][1][0],
+            0.0,
+        )
 
         (_, c_data1), _ = prepare_concurrence_data(solver1.result, (0, 1))
         (_, c_data2), _ = prepare_concurrence_data(solver1.result, (6, 7))