リンドブラッド方程式のモンテカルロソルバ

pysrc/qulacs/__init__.pyi

@@ -111,6 +111,10 @@ class GeneralQuantumOperator():
         """
         Get expectation value
         """
+    def get_expectation_value_single_thread(self, state: QuantumStateBase) -> complex: 
+        """
+        Get expectation value
+        """
     def get_qubit_count(self) -> int: 
         """
         Get qubit count
@@ -189,6 +193,10 @@ class Observable(GeneralQuantumOperator):
         """
         Get expectation value
         """
+    def get_expectation_value_single_thread(self, state: QuantumStateBase) -> float: 
+        """
+        Get expectation value
+        """
     def get_qubit_count(self) -> int: 
         """
         Get qubit count
@@ -455,6 +463,10 @@ class PauliOperator():
         """
         Get expectation value
         """
+    def get_expectation_value_single_thread(self, state: QuantumStateBase) -> complex: 
+        """
+        Get expectation value
+        """
     def get_index_list(self) -> typing.List[int]: 
         """
         Get list of target qubit indices

pysrc/qulacs/gate.pyi

@@ -28,6 +28,7 @@ def Identity(index: int) -> qulacs_core.QuantumGateBase: ...
 def IndependentXZNoise(index: int, prob: float) -> qulacs_core.QuantumGateBase: ...
 def Instrument(kraus_list: List[qulacs_core.QuantumGateBase], register: int) -> qulacs_core.QuantumGateBase: ...
 def Measurement(index: int, register: int) -> qulacs_core.QuantumGateBase: ...
+def NoisyEvolution(hamiltonian: qulacs_core.Observable, c_ops: List[qulacs_core.GeneralQuantumOperator], time: float, dt: float) -> qulacs_core.QuantumGateBase: ...
 def P0(index: int) -> qulacs_core.QuantumGateBase: ...
 def P1(index: int) -> qulacs_core.QuantumGateBase: ...
 def ParametricPauliRotation(index_list: List[int], pauli_ids: List[int], angle: float) -> qulacs_core.QuantumGate_SingleParameter: ...

pysrc/qulacs/gate/__init__.pyi

@@ -24,6 +24,7 @@ __all__ = [
     "IndependentXZNoise",
     "Instrument",
     "Measurement",
+    "NoisyEvolution",
     "P0",
     "P1",
     "ParametricPauliRotation",
@@ -144,6 +145,10 @@ def Measurement(index: int, register: int) -> qulacs_core.QuantumGateBase:
     """
     Create measurement gate
     """
+def NoisyEvolution(hamiltonian: qulacs_core.Observable, c_ops: typing.List[qulacs_core.GeneralQuantumOperator], time: float, dt: float) -> qulacs_core.QuantumGateBase:
+    """
+    Create noisy evolution
+    """
 def P0(index: int) -> qulacs_core.QuantumGateBase:
     """
     Create projection gate to |0> subspace

python/cppsim_wrapper.cpp

@@ -49,6 +49,7 @@ PYBIND11_MODULE(qulacs_core, m) {
         .def("get_coef", &PauliOperator::get_coef, "Get coefficient of Pauli term")
         .def("add_single_Pauli", &PauliOperator::add_single_Pauli, "Add Pauli operator to this term", py::arg("index"), py::arg("pauli_string"))
         .def("get_expectation_value", &PauliOperator::get_expectation_value, "Get expectation value", py::arg("state"))
+        .def("get_expectation_value_single_thread", &PauliOperator::get_expectation_value_single_thread, "Get expectation value", py::arg("state"))
         .def("get_transition_amplitude", &PauliOperator::get_transition_amplitude, "Get transition amplitude", py::arg("state_bra"), py::arg("state_ket"))
         .def("copy", &PauliOperator::copy, pybind11::return_value_policy::take_ownership, "Create copied instance of Pauli operator class")
         .def("get_pauli_string", &PauliOperator::get_pauli_string, "get pauli string")
@@ -78,6 +79,7 @@ PYBIND11_MODULE(qulacs_core, m) {
             return quantum_operator.get_term(index)->copy();
         }, pybind11::return_value_policy::take_ownership, "Get Pauli term", py::arg("index"))
         .def("get_expectation_value", &GeneralQuantumOperator::get_expectation_value, "Get expectation value", py::arg("state"))
+        .def("get_expectation_value_single_thread", &GeneralQuantumOperator::get_expectation_value, "Get expectation value", py::arg("state"))
         .def("get_transition_amplitude", &GeneralQuantumOperator::get_transition_amplitude, "Get transition amplitude", py::arg("state_bra"), py::arg("state_ket"))
         .def("__str__", &GeneralQuantumOperator::to_string, "to string")
         //.def_static("get_split_GeneralQuantumOperator", &(GeneralQuantumOperator::get_split_observable));
@@ -118,6 +120,9 @@ PYBIND11_MODULE(qulacs_core, m) {
         .def("get_expectation_value", [](const HermitianQuantumOperator& observable, const QuantumStateBase* state) {
                                           double res = observable.get_expectation_value(state).real();
                                           return res;}, "Get expectation value", py::arg("state"))
+        .def("get_expectation_value_single_thread", [](const HermitianQuantumOperator& observable, const QuantumStateBase* state) {
+                                          double res = observable.get_expectation_value_single_thread(state).real();
+                                          return res;}, "Get expectation value", py::arg("state"))
         .def("get_transition_amplitude", &HermitianQuantumOperator::get_transition_amplitude, "Get transition amplitude", py::arg("state_bra"), py::arg("state_ket"))
         //.def_static("get_split_Observable", &(HermitianQuantumOperator::get_split_observable));
         .def("add_random_operator", &HermitianQuantumOperator::add_random_operator, "Add random pauli operator", py::arg("operator_count"))
@@ -165,7 +170,6 @@ PYBIND11_MODULE(qulacs_core, m) {
         .def("to_string",&QuantumState::to_string, "Get string representation")
         .def("sampling", (std::vector<ITYPE> (QuantumState::*)(UINT))&QuantumState::sampling, "Sampling measurement results", py::arg("count"))
         .def("sampling", (std::vector<ITYPE>(QuantumState::*)(UINT, UINT))&QuantumState::sampling, "Sampling measurement results", py::arg("count"), py::arg("seed"))
-
         .def("get_vector", [](const QuantumState& state) {
         Eigen::VectorXcd vec = Eigen::Map<Eigen::VectorXcd>(state.data_cpp(), state.dim);
         return vec;
@@ -507,6 +511,7 @@ PYBIND11_MODULE(qulacs_core, m) {
     mgate.def("Instrument", &gate::Instrument, pybind11::return_value_policy::take_ownership, "Create instruments", py::arg("kraus_list"), py::arg("register"));
     mgate.def("Adaptive", py::overload_cast<QuantumGateBase *, std::function<bool(const std::vector<unsigned int>&)>> (&gate::Adaptive), pybind11::return_value_policy::take_ownership, "Create adaptive gate", py::arg("gate"), py::arg("condition"));
     mgate.def("Adaptive",  py::overload_cast<QuantumGateBase *, std::function<bool(const std::vector<unsigned int>&, unsigned int)>, unsigned int> (&gate::Adaptive), pybind11::return_value_policy::take_ownership, "Create adaptive gate", py::arg("gate"), py::arg("condition"),py::arg("id"));
+    mgate.def("NoisyEvolution", &gate::NoisyEvolution, pybind11::return_value_policy::take_ownership, "Create noisy evolution", py::arg("hamiltonian"), py::arg("c_ops"), py::arg("time"), py::arg("dt"));
 
     py::class_<QuantumGate_SingleParameter, QuantumGateBase>(m, "QuantumGate_SingleParameter")
         .def("get_parameter_value", &QuantumGate_SingleParameter::get_parameter_value, "Get parameter value")

src/cppsim/state_dm.hpp

@@ -342,6 +342,7 @@ class DensityMatrixCpu : public QuantumStateBase {
 
     /**
      * \~japanese-en 量子状態を足しこむ
+     * TODO: implement this in single_thread
      */
     virtual void add_state_with_coef_single_thread(
         CPPCTYPE coef, const QuantumStateBase* state) override {

test/cppsim/test_noisyevolution.cpp

@@ -196,6 +196,7 @@ TEST(NoisyEvolutionTest, T1T2) {
     double ref = -0.3135191750739427;  // generated by qutip
     UINT n = 2;
     UINT n_samples = 100;
+
     Observable observable(n);
     observable.add_operator(1, "X 0");
     // create hamiltonian and collapse operator
@@ -223,7 +224,6 @@ TEST(NoisyEvolutionTest, T1T2) {
         gate::H(0)->update_quantum_state(&state);
         gate::H(1)->update_quantum_state(&state);
         circuit.update_quantum_state(&state);
-        state.normalize(state.get_squared_norm());
         exp += observable.get_expectation_value(&state).real() / n_samples;
     }
     std::cout << "NoisyEvolution: " << exp << " ref: " << ref << std::endl;