Support multiple input/output states.

src/evaluate_circuit.h

@@ -50,6 +50,15 @@ namespace qflex {
 std::vector<std::vector<std::size_t>> read_grid_layout_from_stream(
     std::istream* grid_data, std::size_t I, std::size_t J);
 
+struct QflexFinalQubits {
+  std::vector<std::vector<std::size_t>> qubits;
+  std::vector<std::size_t> output_pos_map;
+  std::vector<std::vector<std::size_t>> output_values_map;
+};
+
+QflexFinalQubits get_final_qubits(
+    const QflexGrid& grid, const std::list<ContractionOperation>& ordering);
+
 /**
  * Determines the final qubit positions and output states for a given ordering.
  * @param input QflexInput generated from the command line.
@@ -58,12 +67,20 @@ std::vector<std::vector<std::size_t>> read_grid_layout_from_stream(
  * be populated by this method.
  * @param output_states vector of output states for the given contraction
  * ordering, to be populated by this method.
- * @return the final state vector, with 'x' for cut locations.
  */
-std::string get_output_states(
-    const QflexInput* input, const std::list<ContractionOperation>& ordering,
-    std::vector<std::vector<std::size_t>>* final_qubits,
-    std::vector<std::string>* output_states);
+std::vector<std::string> get_output_states(
+    const std::string& base_state, const QflexFinalQubits& final_qubits);
+
+void apply_terminal_cuts(
+    const QflexGrid& grid, const QflexFinalQubits& final_qubits,
+    std::vector<std::vector<std::vector<Tensor>>>* tensor_grid_3D_ptr);
+
+void apply_delta_output(
+    const QflexGrid& grid, const std::string& final_state,
+    const QflexFinalQubits& final_qubits,
+    const std::vector<std::vector<std::vector<Tensor>>>& tensor_grid_3D,
+    std::vector<std::vector<Tensor>>* tensor_grid_prt,
+    std::vector<s_type>* scratch_2D_ptr);
 
 /**
  * Evaluates a circuit and returns the final amplitudes of each state resulting
@@ -78,8 +95,8 @@ std::string get_output_states(
  * States for qubits with terminal cuts are listed at the end of the state
  * bitstring, in the order of their terminal cuts.
  */
-std::vector<std::pair<std::string, std::complex<double>>> EvaluateCircuit(
-    QflexInput* input);
+std::vector<std::tuple<std::string, std::string, std::complex<double>>>
+EvaluateCircuit(const QflexInput& input);
 
 }  // namespace qflex
 

src/input.h

@@ -12,8 +12,8 @@ struct QflexInput {
   QflexOrdering ordering;
   QflexCircuit circuit;
   QflexGrid grid;
-  std::string initial_state;
-  std::string final_state;
+  std::vector<std::string> initial_states;
+  std::vector<std::string> final_states;
 };
 
 }  // namespace qflex

src/main.cpp

@@ -32,8 +32,8 @@ tensor network, CPU-based simulator of large quantum circuits.
     qflex::utils::readable_memory_string(qflex::global::memory_limit), R"(].
     -t,--track-memory=<seconds>            If <verbosity_level> > 0, track memory usage [default: )",
     qflex::global::track_memory_seconds, R"(].
-    --initial-conf=<initial_conf>          Initial configuration.
-    --final-conf=<final_conf>              Final configuration.
+    --initial-conf=<initial_conf>          Initial configuration [default: 00...00].
+    --final-conf=<final_conf>              Final configuration [default: 00...00].
     --version                              Show version.
 )");
 
@@ -72,17 +72,6 @@ int main(int argc, char** argv) {
       qflex::global::track_memory_seconds =
           args["<track_memory_seconds>"].asLong();
 
-    // Get initial/final configurations
-    if (static_cast<bool>(args["--initial-conf"]))
-      input.initial_state = args["--initial-conf"].asString();
-    else if (static_cast<bool>(args["<initial_conf>"]))
-      input.initial_state = args["<initial_conf>"].asString();
-
-    if (static_cast<bool>(args["--final-conf"]))
-      input.final_state = args["--final-conf"].asString();
-    else if (static_cast<bool>(args["<final_conf>"]))
-      input.final_state = args["<final_conf>"].asString();
-
     // Getting filenames
     std::string circuit_filename = static_cast<bool>(args["--circuit"])
                                        ? args["--circuit"].asString()
@@ -95,6 +84,63 @@ int main(int argc, char** argv) {
                                     ? args["--grid"].asString()
                                     : args["<grid_filename>"].asString();
 
+    // Load circuit
+    input.circuit.load(std::ifstream(circuit_filename));
+
+    // Load ordering
+    input.ordering.load(std::ifstream(ordering_filename));
+
+    // Load grid
+    input.grid.load(grid_filename);
+
+    // Get initial/final configurations
+    for (const auto& arg : {"--initial-conf", "<initial_conf>"})
+      if (static_cast<bool>(args[arg])) {
+        std::stringstream ss(args[arg].asString());
+        std::string val;
+        while (std::getline(ss, val, ',')) {
+          if (val == "00...00")
+            input.initial_states.push_back(
+                std::string(input.circuit.num_active_qubits, '0'));
+          else if (val.find_first_not_of("01") != std::string::npos)
+            throw ERROR_MSG(
+                "Initial configurations can only have 0 or 1 characters.");
+          else if (std::size(val) != input.circuit.num_active_qubits)
+            throw ERROR_MSG(
+                "Initial configurations must have the number of active "
+                "qubits.");
+          else
+            input.initial_states.push_back(val);
+        }
+      }
+
+    for (const auto& arg : {"--final-conf", "<final_conf>"})
+      if (static_cast<bool>(args[arg])) {
+        std::stringstream ss(args[arg].asString());
+        std::string val;
+        while (std::getline(ss, val, ',')) {
+          if (val == "00...00")
+            input.final_states.push_back(
+                std::string(input.circuit.num_active_qubits, '0'));
+          else if (val.find_first_not_of("01") != std::string::npos)
+            throw ERROR_MSG(
+                "Final configurations can only have 0 or 1 characters.");
+          else if (std::size(val) != input.circuit.num_active_qubits)
+            throw ERROR_MSG(
+                "Final configurations must have the number of active qubits.");
+          else
+            input.final_states.push_back(val);
+        }
+      }
+
+    // Delete duplicate initial/final configuration
+    input.initial_states.erase(std::unique(std::begin(input.initial_states),
+                                           std::end(input.initial_states)),
+                               std::end(input.initial_states));
+    input.final_states.erase(std::unique(std::begin(input.final_states),
+                                         std::end(input.final_states)),
+                             std::end(input.final_states));
+
     // Print OMP_NUM_THREADS and MKL_NUM_THREADS
     if (qflex::global::verbose > 0)
       for (const char* var : {"OMP_NUM_THREADS", "MKL_NUM_THREADS"})
@@ -114,19 +160,11 @@ int main(int argc, char** argv) {
       alarm(qflex::global::track_memory_seconds);
     }
 
-    // Load circuit
-    input.circuit.load(std::ifstream(circuit_filename));
-
-    // Load ordering
-    input.ordering.load(std::ifstream(ordering_filename));
-
-    // Load grid
-    input.grid.load(grid_filename);
-
     // Evaluating circuit.
-    std::vector<std::pair<std::string, std::complex<double>>> amplitudes;
+    std::vector<std::tuple<std::string, std::string, std::complex<double>>>
+        amplitudes;
     try {
-      amplitudes = qflex::EvaluateCircuit(&input);
+      amplitudes = qflex::EvaluateCircuit(input);
     } catch (const std::string& err_msg) {
       throw ERROR_MSG("Failed to call EvaluateCircuit(). Error:\n\t[", err_msg,
                       "]");
@@ -137,13 +175,9 @@ int main(int argc, char** argv) {
       qflex::memory::print_memory_usage();
 
     // Printing output.
-    for (std::size_t c = 0; c < amplitudes.size(); ++c) {
-      const auto& state = amplitudes[c].first;
-      const auto& amplitude = amplitudes[c].second;
-      std::cout << input.initial_state << " --> " << state << ": "
-                << std::real(amplitude) << " " << std::imag(amplitude)
+    for (const auto& [initial_state, final_state, amplitude] : amplitudes)
+      std::cout << initial_state << " --> " << final_state << ": " << amplitude
                 << std::endl;
-    }
 
   } catch (const std::exception& ex) {
     std::cerr << ex.what() << std::endl;

src/pybind_main.cpp

@@ -67,8 +67,7 @@ std::vector<std::string> LoadStates(const py::dict &options,
       states.push_back(in_state.cast<std::string>());
     } else
       throw ERROR_MSG("states must be strings.");
-  } else
-    states.push_back("");
+  }
   return states;
 }
 
@@ -84,27 +83,39 @@ std::vector<std::pair<std::string, std::complex<double>>> simulate(
     LoadData(options, "grid", input.grid);
     LoadData(options, "circuit", input.circuit);
     LoadData(options, "ordering", input.ordering);
-    const auto initial_states = LoadStates(options, "initial");
-    const auto final_states = LoadStates(options, "final");
 
-    // Define container for amplitudes
-    std::vector<std::pair<
-        std::string, std::vector<std::pair<std::string, std::complex<double>>>>>
-        amplitudes;
+    // Load initial/final states
+    input.initial_states = [&options, &input]() {
+      std::vector<std::string> states = LoadStates(options, "initial");
+      if (std::empty(states))
+        states.push_back(std::string(input.circuit.num_active_qubits, '0'));
+      return states;
+    }();
+    input.final_states = [&options, &input]() {
+      std::vector<std::string> states = LoadStates(options, "final");
+      if (std::empty(states))
+        states.push_back(std::string(input.circuit.num_active_qubits, '0'));
+      return states;
+    }();
 
-    // TODO: Pybind11 does not allow multiple types as output
-    if (std::size(initial_states) != 1 or std::size(final_states) != 1)
+    // TODO: Not yet supported in Cirq
+    if (std::size(input.initial_states) != 1 ||
+        std::size(input.final_states) != 1)
       throw ERROR_MSG("Not yet supported");
 
-    for (const auto &is : initial_states) {
-      for (const auto &fs : final_states) {
-        input.initial_state = is;
-        input.final_state = fs;
-        amplitudes.push_back({is, EvaluateCircuit(&input)});
-      }
-    }
+    // Remove duplicate initial/final states
+    input.initial_states.erase(std::unique(std::begin(input.initial_states),
+                                           std::end(input.initial_states)),
+                               std::end(input.initial_states));
+    input.final_states.erase(std::unique(std::begin(input.final_states),
+                                         std::end(input.final_states)),
+                             std::end(input.final_states));
+
+    // Define container for amplitudes
+    std::vector<std::tuple<std::string, std::string, std::complex<double>>>
+        amplitudes = EvaluateCircuit(input);
 
-    return std::get<1>(amplitudes[0]);
+    return {{std::get<1>(amplitudes[0]), std::get<2>(amplitudes[0])}};
 
     // Gently return an error msg if exception is known. Otherwise, rethrow
     // exception.