Native Simulator: change prereq check in state injection to return false instead of throw

src/Simulation/Native/src/simulator/capi.cpp

@@ -45,7 +45,7 @@ extern "C"
         return Microsoft::Quantum::Simulator::get(id)->JointEnsembleProbability(bv, qv);
     }
 
-    MICROSOFT_QUANTUM_DECL void InjectState(
+    MICROSOFT_QUANTUM_DECL bool InjectState(
         _In_ unsigned sid,
         _In_ unsigned n,
         _In_reads_(n) unsigned* q,
@@ -61,7 +61,7 @@ extern "C"
         }
         std::vector<unsigned> qubits(q, q + n);
 
-        Microsoft::Quantum::Simulator::get(sid)->InjectState(qubits, amplitudes);
+        return Microsoft::Quantum::Simulator::get(sid)->InjectState(qubits, amplitudes);
     }
 
     MICROSOFT_QUANTUM_DECL void allocateQubit(_In_ unsigned id, _In_ unsigned q)

src/Simulation/Native/src/simulator/capi.hpp

@@ -35,7 +35,7 @@ extern "C"
         _In_reads_(n) int* b,
         _In_reads_(n) unsigned* q);
 
-    MICROSOFT_QUANTUM_DECL void InjectState(
+    MICROSOFT_QUANTUM_DECL bool InjectState(
         _In_ unsigned sid,
         _In_ unsigned n,
         _In_reads_(n) unsigned* q, // The listed qubits must be unentangled and in state |0>

src/Simulation/Native/src/simulator/local_test.cpp

@@ -390,7 +390,7 @@ TEST_CASE("permute_basis", "[local_test]")
     const double amp = 1.0 / std::sqrt(5);
     std::vector<ComplexType> amplitudes = {{amp, 0.0}, {amp, 0.0}, {amp, 0.0}, {0.0, 0.0},
                                            {0.0, 0.0}, {0.0, 0.0}, {amp, 0.0}, {amp, 0.0}};
-    psi.inject_state({q0, q1, q2}, amplitudes);
+    REQUIRE(psi.inject_state({q0, q1, q2}, amplitudes));
 
     SECTION("identity permutation")
     {
@@ -453,7 +453,7 @@ TEST_CASE("permute_basis depends on the order of logical qubits (2)", "[local_te
     // Inject state, which would allow us to easily check permutations. It's not a normalized state but for this
     // test it doesn't matter.
     std::vector<ComplexType> amplitudes = {{0.0, 0.0}, {1.0, 0.0}, {2.0, 0.0}, {3.0, 0.0}};
-    psi.inject_state({q0, q1}, amplitudes);
+    REQUIRE(psi.inject_state({q0, q1}, amplitudes));
     // after the state injection, positions of the qubits are q0:0 and q1:1
 
     SECTION("q0-q1 order (matches the current positions of the qubits in the standard basis)")
@@ -496,7 +496,7 @@ TEST_CASE("permute_basis depends on the order of logical qubits (3)", "[local_te
     // test it doesn't matter.
     std::vector<ComplexType> amplitudes = {{0.0, 0.0}, {1.0, 0.0}, {2.0, 0.0}, {3.0, 0.0},
                                            {4.0, 0.0}, {5.0, 0.0}, {6.0, 0.0}, {7.0, 0.0}};
-    psi.inject_state({q0, q1, q2}, amplitudes);
+    REQUIRE(psi.inject_state({q0, q1, q2}, amplitudes));
 
     SECTION("q0-q1-q2 order")
     {
@@ -544,7 +544,7 @@ TEST_CASE("Inject total cat state", "[local_test]")
     std::vector<ComplexType> amplitudes = {{amp, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {amp, 0.0}};
     REQUIRE(amplitudes.size() == N);
 
-    sim.InjectState(qs, amplitudes);
+    REQUIRE(sim.InjectState(qs, amplitudes));
 
     // undo the injected state back to |00>
     sim.CX({qs[0]}, qs[1]);
@@ -571,13 +571,13 @@ TEST_CASE("Should fail to inject state if qubits aren't all |0>", "[local_test]"
 
     // unentangled but not |0>
     sim.H(qs[1]);
-    REQUIRE_THROWS(sim.InjectState(qs, amplitudes));
-    REQUIRE_THROWS(sim.InjectState({qs[0], qs[1]}, amplitudes_sub));
+    REQUIRE_FALSE(sim.InjectState(qs, amplitudes));
+    REQUIRE_FALSE(sim.InjectState({qs[0], qs[1]}, amplitudes_sub));
 
     // entanglement doesn't make things any better
     sim.CX({qs[1]}, qs[2]);
-    REQUIRE_THROWS(sim.InjectState(qs, amplitudes));
-    REQUIRE_THROWS(sim.InjectState({qs[0], qs[1]}, amplitudes_sub));
+    REQUIRE_FALSE(sim.InjectState(qs, amplitudes));
+    REQUIRE_FALSE(sim.InjectState({qs[0], qs[1]}, amplitudes_sub));
 }
 
 TEST_CASE("Inject total state on reordered qubits", "[local_test]")
@@ -599,7 +599,7 @@ TEST_CASE("Inject total state on reordered qubits", "[local_test]")
 
     // Notice, that we are listing the qubits in order that doesn't match their allocation order. We are saying here,
     // that InjectState should create Bell pair from qs[1] and qs[2]!
-    sim.InjectState({qs[1], qs[2], qs[0]}, amplitudes);
+    REQUIRE(sim.InjectState({qs[1], qs[2], qs[0]}, amplitudes));
     REQUIRE((sim.isclassical(qs[0]) && !sim.M(qs[0])));
 
     // undo the state change and check that the whole system is back to |000>
@@ -652,7 +652,7 @@ TEST_CASE("Inject state on two qubits out of three", "[local_test]")
         sim.H(q0);
     }
 
-    sim.InjectState({x, y}, amplitudes);
+    REQUIRE(sim.InjectState({x, y}, amplitudes));
 
     // undo the state injection with quantum op and check that the qubits we injected state for are back to |0>
     sim.H(x);
@@ -697,7 +697,7 @@ TEST_CASE("Perf of injecting equal superposition state", "[skip]") // local micr
         std::vector<ComplexType> amplitudes(N, {amp, 0.0});
 
         auto start = high_resolution_clock::now();
-        sim.InjectState(qs, amplitudes);
+        REQUIRE(sim.InjectState(qs, amplitudes));
         sim.M(qs[0]); // to have the same overhead compared to preparation test case
         std::cout << "    Total state injection:\t";
         std::cout << duration_cast<microseconds>(high_resolution_clock::now() - start).count();
@@ -712,7 +712,7 @@ TEST_CASE("Perf of injecting equal superposition state", "[skip]") // local micr
         std::vector<ComplexType> amplitudes(N, {amp, 0.0});
 
         auto start = std::chrono::high_resolution_clock::now();
-        sim.InjectState(qs, amplitudes);
+        REQUIRE(sim.InjectState(qs, amplitudes));
         sim.H(q_last);
         sim.M(qs[0]); // to have the same overhead compared to preparation test case
         std::cout << "  Partial state injection:\t";
@@ -761,7 +761,7 @@ TEST_CASE("Perf of injecting cat state", "[skip]") // local micro_benchmark
         amplitudes[N - 1] = {amp, 0.0};
 
         auto start = std::chrono::high_resolution_clock::now();
-        sim.InjectState(qs, amplitudes);
+        REQUIRE(sim.InjectState(qs, amplitudes));
         sim.M(qs[0]); // to have the same overhead compared to preparation test case
         std::cout << "    Total cat state injection:\t";
         std::cout << duration_cast<microseconds>(high_resolution_clock::now() - start).count();
@@ -778,7 +778,7 @@ TEST_CASE("Perf of injecting cat state", "[skip]") // local micro_benchmark
         amplitudes[N - 1] = {amp, 0.0};
 
         auto start = std::chrono::high_resolution_clock::now();
-        sim.InjectState(qs, amplitudes);
+        REQUIRE(sim.InjectState(qs, amplitudes));
         sim.CX({qs[0]}, q_last);
         sim.M(qs[0]); // to have the same overhead compared to preparation test case
         std::cout << "  Partial cat state injection:\t";

src/Simulation/Native/src/simulator/simulator.hpp

@@ -59,10 +59,10 @@ class Simulator : public Microsoft::Quantum::Simulator::SimulatorInterface
         return p;
     }
 
-    void InjectState(const std::vector<logical_qubit_id>& qubits, const std::vector<ComplexType>& amplitudes)
+    bool InjectState(const std::vector<logical_qubit_id>& qubits, const std::vector<ComplexType>& amplitudes)
     {
         recursive_lock_type l(mutex());
-        psi.inject_state(qubits, amplitudes);
+        return psi.inject_state(qubits, amplitudes);
     }
 
     bool isclassical(logical_qubit_id q)

src/Simulation/Native/src/simulator/simulatorinterface.hpp

@@ -30,7 +30,7 @@ class SimulatorInterface
 
     virtual double JointEnsembleProbability(std::vector<Gates::Basis> bs, std::vector<unsigned> qs) = 0;
 
-    virtual void InjectState(
+    virtual bool InjectState(
         const std::vector<logical_qubit_id>& qubits,
         const std::vector<ComplexType>& amplitudes) = 0;
 

src/Simulation/Native/src/simulator/wavefunction.hpp

@@ -565,27 +565,27 @@ class Wavefunction
         return kernels::jointprobability(wfn_, bs, get_qubit_positions(qs));
     }
 
-    /// \pre: Each qubit, listed in `q`, must be unentangled and in state |0>.
+    /// \pre: Each qubit, listed in `q`, must be unentangled and in state |0>. If the prerequisite isn't satisfied,
+    /// the method returns `false` and leaves the state of the system unchanged.
     /// Place qubits, listed in `q` into superposition of basis vectors with provided `amplitudes`, where the order of
-    /// qubits in array `q` defines the standard computational basis in little endian order.
-    void inject_state(const std::vector<logical_qubit_id>& qubits, const std::vector<ComplexType>& amplitudes)
+    /// qubits in array `q` defines the standard computational basis in little endian order. Returns `true` if the state
+    /// is successfuly injected.
+    bool inject_state(const std::vector<logical_qubit_id>& qubits, const std::vector<ComplexType>& amplitudes)
     {
         assert((static_cast<size_t>(1) << qubits.size()) == amplitudes.size());
 
         flush();
 
-        // Check prerequisites.
-        std::vector<positional_qubit_id> positions = get_qubit_positions(qubits);
-        for (positional_qubit_id p : positions)
+        if (qubits.size() == num_qubits_)
         {
-            if (!kernels::isclassical(wfn_, p) || kernels::getvalue(wfn_, p) != 0)
+            // Check prerequisites. In the case of total state injection the wave function must consist of a single
+            // term |0...0> (so we can avoid checking each qubit individually).
+            double eps = 100. * std::numeric_limits<double>::epsilon();
+            if (std::norm(wfn_[0]) < 1.0 - eps)
             {
-                throw std::runtime_error("Cannot prepare state of entangled qubits or if they are not in state |0>");
+                return false;
             }
-        }
 
-        if (qubits.size() == num_qubits_)
-        {
             // For full state injection we can copy the user's wave function into our store and reorder the
             // positions map without doing any math.
             for (unsigned i = 0; i < qubits.size(); i++)
@@ -596,6 +596,16 @@ class Wavefunction
         }
         else
         {
+            // Check prerequisites.
+            std::vector<positional_qubit_id> positions = get_qubit_positions(qubits);
+            for (positional_qubit_id p : positions)
+            {
+                if (!kernels::isclassical(wfn_, p) || kernels::getvalue(wfn_, p) != 0)
+                {
+                    return false;
+                }
+            }
+
             // The current state can be thought of as Sum(a_i*|i>|0...0>), after the state injection it will become
             // Sum(a_i*|i>Sum(b_j*|j>)) = Sum(a_i*b_j|i>|j>). Thus, to compute amplitude of a term |k> after the state
             // injection we need to find the corresponding |i> vector from the original wave function and |j> vector
@@ -619,6 +629,8 @@ class Wavefunction
             }
             std::swap(wfn_, wfn_new);
         }
+
+        return true;
     }
 
     /// measure a qubit