Added a test for circuit conversion to two qubit gates layers and back. Found and fixed the issue with some circuits (containing classically controlled gates) being converted like that for mps initial qubits layout optimization in the network execution.

Author: aromanro

Circuit/Circuit.h

@@ -2192,6 +2192,7 @@ class Circuit : public IOperation<Time> {
     layers.emplace_back(std::make_shared<Circuits::Circuit<Time>>());
 
     std::unordered_map<Types::qubit_t, Types::qubit_t> qubitsUsed;
+    std::unordered_map<size_t, size_t> classicalBitLayer;
 
     for (const auto &op : GetOperations()) {
       // check the instruction, see if a new layer is needed
@@ -2213,10 +2214,30 @@ class Circuit : public IOperation<Time> {
           }
         }
 
+        if (op->IsConditional()) {
+          const auto bits = op->AffectedBits();
+          for (const auto bit : bits)
+            maxLevel = std::max(maxLevel, classicalBitLayer[bit]);
+        }
+
         // now set all the qubits in the instruction to the max level
         for (Types::qubit_t qbit : qubits) qubitsUsed[qbit] = maxLevel;
 
-        layers[maxLevel > 0 ? maxLevel - 1 : 0]->AddOperation(op->Clone());
+        const size_t layerIdx = maxLevel > 0 ? maxLevel - 1 : 0;
+
+        // ensure enough layers exist for the computed level
+        while (layers.size() <= layerIdx)
+          layers.push_back(std::make_shared<Circuits::Circuit<Time>>());
+
+        layers[layerIdx]->AddOperation(op->Clone());
+
+        const auto writtenBits = op->AffectedBits();
+        if (!writtenBits.empty() && !op->IsConditional()) {
+          const size_t writtenLevel = maxLevel > 0 ? maxLevel : 1;
+          for (const auto bit : writtenBits)
+            classicalBitLayer[bit] =
+                std::max(classicalBitLayer[bit], writtenLevel);
+        }
       } else
         // add the instruction to the last layer
         layers.back()->AddOperation(op->Clone());
@@ -2242,6 +2263,11 @@ class Circuit : public IOperation<Time> {
 
     std::unordered_map<Types::qubit_t, Types::qubit_t> qubitsUsed;
 
+    // Track which layer (1-based, like qubitsUsed) each classical bit was last
+    // written to, so that conditional operations reading those bits are placed
+    // in the same layer or later.
+    std::unordered_map<size_t, size_t> classicalBitLayer;
+
     for (const auto &op : GetOperations()) {
       // check the instruction, see if a new layer is needed
 
@@ -2262,10 +2288,34 @@ class Circuit : public IOperation<Time> {
           }
         }
 
+        // For conditional operations, ensure this op is placed no earlier than
+        // the layer where the classical bits it depends on were written.
+        if (op->IsConditional()) {
+          const auto bits = op->AffectedBits();
+          for (const auto bit : bits)
+            maxLevel = std::max(maxLevel, classicalBitLayer[bit]);
+        }
+
         // now set all the qubits in the instruction to the max level
         for (Types::qubit_t qbit : qubits) qubitsUsed[qbit] = maxLevel;
 
-        layers[maxLevel > 0 ? maxLevel - 1 : 0]->AddOperation(op);
+        const size_t layerIdx = maxLevel > 0 ? maxLevel - 1 : 0;
+
+        // ensure enough layers exist for the computed level
+        while (layers.size() <= layerIdx)
+          layers.push_back(std::make_shared<Circuits::Circuit<Time>>());
+
+        layers[layerIdx]->AddOperation(op);
+
+        // Record the layer for classical bits written by measurements / resets
+        // so that later conditional ops respect the dependency.
+        const auto writtenBits = op->AffectedBits();
+        if (!writtenBits.empty() && !op->IsConditional()) {
+          const size_t writtenLevel = maxLevel > 0 ? maxLevel : 1;
+          for (const auto bit : writtenBits)
+            classicalBitLayer[bit] =
+                std::max(classicalBitLayer[bit], writtenLevel);
+        }
       } else
         // add the instruction to the last layer
         layers.back()->AddOperation(op);

Network/NetworkJob.h

@@ -18,6 +18,8 @@
 
 #include "../Simulators/MPSDummySimulator.h"
 
+#include "Network.h"
+
 namespace Network {
 
 template <typename Time = Types::time_type>
@@ -330,6 +332,8 @@ class ExecuteJob {
   std::string maxBondDim;
   std::string singularValueThreshold;
   std::string mpsSample;
+
+  std::shared_ptr<Network::INetwork<Time>> network;
 };
 
 }  // namespace Network

Network/SimpleDisconnectedNetwork.h

@@ -647,6 +647,8 @@ class SimpleDisconnectedNetwork : public INetwork<Time> {
         job->mpsSample = mpsSample;
         job->singularValueThreshold = singularValueThreshold;
 
+        job->network = BaseClass::getptr();
+
         if (optSim) {
           job->optSim = optSim->Clone();
           job->executedGates = executed;
@@ -669,6 +671,8 @@ class SimpleDisconnectedNetwork : public INetwork<Time> {
       job->mpsSample = mpsSample;
       job->singularValueThreshold = singularValueThreshold;
 
+      job->network = BaseClass::getptr();
+
       if (optSim) {
         optSim->SetMultithreading(true);
         job->optSim = optSim;
@@ -829,6 +833,8 @@ class SimpleDisconnectedNetwork : public INetwork<Time> {
         job->mpsSample = mpsSample;
         job->singularValueThreshold = singularValueThreshold;
 
+        job->network = BaseClass::getptr();
+
         if (optSim) {
           job->optSim = optSim->Clone();
           job->executedGates = executed;
@@ -851,6 +857,8 @@ class SimpleDisconnectedNetwork : public INetwork<Time> {
       job->mpsSample = mpsSample;
       job->singularValueThreshold = singularValueThreshold;
 
+      job->network = BaseClass::getptr();
+
       if (optSim) {
         optSim->SetMultithreading(true);
         job->optSim = optSim;
@@ -2142,9 +2150,8 @@ class SimpleDisconnectedNetwork : public INetwork<Time> {
 
           if (optimizeInitialQubitsMap) sim->SetInitialQubitsMap(optimalMap);
 
-          // TODO: this breaks scheduling, needs to be done before scheduling, but at that point we don't know it mps swaps optimization is going to be used or not
-          //auto optCirc = Circuits::Circuit<Time>::LayersToCircuit(layers);
-          //dcirc->SetOperations(optCirc->GetOperations());
+          auto optCirc = Circuits::Circuit<Time>::LayersToCircuit(layers);
+          dcirc->SetOperations(optCirc->GetOperations());
         }
       }
     }

tests/mpsswapsopttests.cpp

@@ -39,6 +39,83 @@ constexpr std::array numGates{10,  20,  30,  50,  80,  100, 150,
                               200, 270, 340, 500,     800, 1000};
 
 
+
+BOOST_DATA_TEST_CASE(ConvertForCuttingLayersRoundtripMatchesStatevector,
+                     numGates, nrGates) {
+  constexpr int nrQubits = 10;
+  const size_t nrStates = 1ULL << nrQubits;
+
+  // build a random circuit with 1- and 2-qubit gates
+  auto randomCirc = std::make_shared<Circuits::Circuit<>>();
+
+  std::mt19937 g(std::random_device{}());
+  std::uniform_real_distribution<double> dblDist(-2. * M_PI, 2. * M_PI);
+  std::uniform_int_distribution<int> gateDist(
+      0, static_cast<int>(Circuits::QuantumGateType::kCUGateType));
+
+  for (int gateNr = 0; gateNr < nrGates; ++gateNr) {
+    Types::qubits_vector qubits(nrQubits);
+    std::iota(qubits.begin(), qubits.end(), 0);
+    std::shuffle(qubits.begin(), qubits.end(), g);
+    auto q1 = qubits[0];
+    auto q2 = qubits[1];
+
+    const double param1 = dblDist(g);
+    const double param2 = dblDist(g);
+    const double param3 = dblDist(g);
+    const double param4 = dblDist(g);
+
+    Circuits::QuantumGateType gateType =
+        static_cast<Circuits::QuantumGateType>(gateDist(g));
+
+    auto theGate = Circuits::CircuitFactory<>::CreateGate(
+        gateType, q1, q2, 0, param1, param2, param3, param4);
+    randomCirc->AddOperation(theGate);
+  }
+
+  // --- Reference: execute the original circuit on a statevector simulator ---
+  auto svSim = Simulators::SimulatorsFactory::CreateSimulator(
+      Simulators::SimulatorType::kQCSim,
+      Simulators::SimulationType::kStatevector);
+  svSim->AllocateQubits(nrQubits);
+  svSim->Initialize();
+
+  Circuits::OperationState stateSV;
+  stateSV.AllocateBits(nrQubits);
+  randomCirc->Execute(svSim, stateSV);
+
+  // --- Convert the circuit the same way OptimizeMPSInitialQubitsMap does ---
+  auto convertedCirc = std::make_shared<Circuits::Circuit<>>(*randomCirc);
+  convertedCirc->ConvertForCutting();
+  auto layers = convertedCirc->ToMultipleQubitsLayersNoClone();
+  auto finalCirc = Circuits::Circuit<>::LayersToCircuit(layers);
+
+  // BOOST_TEST_MESSAGE("Original gate count: " << randomCirc->size());
+  // BOOST_TEST_MESSAGE("Converted gate count: " << finalCirc->size());
+  // BOOST_TEST_MESSAGE("Number of layers: " << layers.size());
+
+  // --- Execute the converted circuit on an MPS simulator (no bond dim cap) ---
+  auto mpsSim = Simulators::SimulatorsFactory::CreateSimulator(
+      Simulators::SimulatorType::kQCSim,
+      Simulators::SimulationType::kMatrixProductState);
+  mpsSim->AllocateQubits(nrQubits);
+  mpsSim->Initialize();
+
+  Circuits::OperationState stateMPS;
+  stateMPS.AllocateBits(nrQubits);
+  finalCirc->Execute(mpsSim, stateMPS);
+
+  // --- Compare amplitudes ---
+  for (size_t s = 0; s < nrStates; ++s) {
+    const auto aSV = svSim->Amplitude(s);
+    const auto aMPS = mpsSim->Amplitude(s);
+    BOOST_CHECK_PREDICATE(
+        checkClose,
+        (aSV)(aMPS)(1e-6));  // no bond dim limit, expect close match
+  }
+}
+
+
 BOOST_DATA_TEST_CASE(OptimalQubitsMapZeroSwapCost, numGates, nrGates) {
   constexpr int nrQubits = 10;