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

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")
@@ -71,11 +72,14 @@ PYBIND11_MODULE(qulacs_core, m) {
         .def("get_qubit_count", &GeneralQuantumOperator::get_qubit_count, "Get qubit count")
         .def("get_state_dim", &GeneralQuantumOperator::get_state_dim, "Get state dimension")
         .def("get_term_count", &GeneralQuantumOperator::get_term_count, "Get count of Pauli terms")
+        .def("apply_to_state", py::overload_cast<QuantumStateBase*, const QuantumStateBase&, QuantumStateBase*>(&GeneralQuantumOperator::apply_to_state, py::const_), "Apply observable to `state_to_be_multiplied`. The result is stored into `dst_state`.",
+            py::arg("work_state"), py::arg("state_to_be_multiplied"), py::arg("dst_state"))
         //.def("get_term", &GeneralQuantumOperator::get_term, pybind11::return_value_policy::take_ownership)
         .def("get_term",[](const GeneralQuantumOperator& quantum_operator, const unsigned int index) {
             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));
@@ -117,6 +121,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"))
@@ -126,7 +133,7 @@ PYBIND11_MODULE(qulacs_core, m) {
             "Compute ground state eigenvalue by power method", py::arg("state"), py::arg("iter_count"), py::arg("mu") = 0.0)
         .def("solve_ground_state_eigenvalue_by_lanczos_method", &HermitianQuantumOperator::solve_ground_state_eigenvalue_by_lanczos_method,
             "Compute ground state eigenvalue by lanczos method", py::arg("state"), py::arg("iter_count"), py::arg("mu") = 0.0)
-        .def("apply_to_state", &HermitianQuantumOperator::apply_to_state, "Apply observable to `state_to_be_multiplied`. The result is stored into `dst_state`.",
+        .def("apply_to_state", py::overload_cast<QuantumStateBase*, const QuantumStateBase&, QuantumStateBase*>(&HermitianQuantumOperator::apply_to_state, py::const_), "Apply observable to `state_to_be_multiplied`. The result is stored into `dst_state`.",
             py::arg("work_state"), py::arg("state_to_be_multiplied"), py::arg("dst_state"))
         .def("__str__", &HermitianQuantumOperator::to_string, "to string")
         ;
@@ -164,11 +171,13 @@ 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;
         }, "Get state vector")
+        .def("get_amplitude", [](const QuantumState& state, const UINT index) -> CPPCTYPE {
+            return state.data_cpp()[index];
+        }, "Get state vector", py::arg("index"))
         .def("get_qubit_count", [](const QuantumState& state) -> unsigned int {return (unsigned int) state.qubit_count; }, "Get qubit count")
         .def("__repr__", [](const QuantumState &p) {return p.to_string();});
         ;
@@ -506,6 +515,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/gate_factory.cpp

@@ -19,6 +19,7 @@
 #include "gate_named_one.hpp"
 #include "gate_named_pauli.hpp"
 #include "gate_named_two.hpp"
+#include "gate_noisy_evolution.hpp"
 #include "gate_reflect.hpp"
 #include "gate_reversible.hpp"
 #include "type.hpp"
@@ -328,6 +329,10 @@ QuantumGateBase* Measurement(
     delete gate1;
     return new_gate;
 }
+QuantumGateBase* NoisyEvolution(Observable* hamiltonian,
+    std::vector<GeneralQuantumOperator*> c_ops, double time, double dt) {
+    return new ClsNoisyEvolution(hamiltonian, c_ops, time, dt);
+}
 
 QuantumGateBase* create_quantum_gate_from_string(std::string gate_string) {
     const char* gateString = gate_string.c_str();

src/cppsim/gate_factory.hpp

@@ -433,4 +433,15 @@ DllExport QuantumGateBase* AmplitudeDampingNoise(
 DllExport QuantumGateBase* Measurement(
     UINT target_index, UINT classical_register_address);
 
+/**
+ * Noisy Evolution
+ * TODO: do this comment
+ *
+ * @param[in] target_index ターゲットとなる量子ビットの添え字
+ * @param[in] classical_register_address 測定値を格納する古典レジスタの場所
+ * @return 作成されたゲートのインスタンス
+ */
+DllExport QuantumGateBase* NoisyEvolution(Observable* hamiltonian,
+    std::vector<GeneralQuantumOperator*> c_ops, double time, double dt = 1e-6);
+
 }  // namespace gate

src/cppsim/gate_noisy_evolution.hpp

@@ -0,0 +1,249 @@
+
+#pragma once
+
+#include "gate.hpp"
+#include "general_quantum_operator.hpp"
+#include "observable.hpp"
+#include "state.hpp"
+#include "utility.hpp"
+
+class ClsNoisyEvolution : public QuantumGateBase {
+private:
+    Random _random;
+    Observable* _hamiltonian;
+    GeneralQuantumOperator* _effective_hamiltonian;
+    std::vector<GeneralQuantumOperator*> _c_ops;  // collapse operator
+    std::vector<GeneralQuantumOperator*>
+        _c_ops_dagger;        // collapse operator dagger
+    double _time;             // evolution time
+    double _dt;               // step size for runge kutta evolution
+    double _norm_tol = 1e-8;  // accuracy in solving <psi|psi>=r
+    int _find_collapse_max_steps =
+        1000;  // maximum number of steps for while loop in _find_collapse
+
+    /**
+     * \~japanese-en collapse が起こるタイミング (norm = rになるタイミング)
+     * を見つける。 この関数内で norm=rになるタイミングまでの evolution
+     * が行われる。
+     */
+    virtual double _find_collapse(QuantumStateBase* k1, QuantumStateBase* k2,
+        QuantumStateBase* k3, QuantumStateBase* k4,
+        QuantumStateBase* prev_state, QuantumStateBase* now_state,
+        double target_norm, double dt) {
+        auto now_norm = now_state->get_squared_norm_single_thread();
+        auto prev_norm = prev_state->get_squared_norm_single_thread();
+        auto t_guess = dt;
+        int search_count = 0;
+        while (std::abs(now_norm - target_norm) > _norm_tol) {
+            // we expect norm to reduce as Ae^-a*dt, so we first calculate the
+            // coefficient a
+            auto a = std::log(now_norm / prev_norm) / t_guess;
+            // then guess the time to become target norm as (target_norm) =
+            // (prev_norm)e^-a*(t_guess) which means -a*t_guess =
+            // log(target_norm/prev_norm)
+            t_guess = std::log(target_norm / prev_norm) / a;
+            // evole by time t_guess
+            now_state->load(prev_state);
+            _evolve_one_step(k1, k2, k3, k4, prev_state, now_state, t_guess);
+            now_norm = now_state->get_squared_norm_single_thread();
+            search_count++;
+            // avoid infinite loop
+            // It sometimes fails to find t_guess to reach the target norm.
+            // More likely to happen when dt is not small enough compared to the
+            // relaxation times
+            if (search_count > _find_collapse_max_steps) {
+                std::cerr << "Failed to find the exact jump time. Try with "
+                             "smaller dt."
+                          << std::endl;
+                return -1.;
+            }
+        }
+        return t_guess;
+    }
+
+    /**
+     * \~japanese-en runge-kutta を 1 step 進める
+     * この関数が終わった時点で、時間発展前の状態は buffer に格納される。
+     */
+    virtual void _evolve_one_step(QuantumStateBase* k1, QuantumStateBase* k2,
+        QuantumStateBase* k3, QuantumStateBase* k4, QuantumStateBase* buffer,
+        QuantumStateBase* state, double dt) {
+        // Runge-Kutta evolution
+        // k1
+        _effective_hamiltonian->apply_to_state_single_thread(state, k1);
+
+        // k2
+        buffer->load(state);
+        buffer->add_state_with_coef_single_thread(-1.i * dt / 2., k1);
+        _effective_hamiltonian->apply_to_state_single_thread(buffer, k2);
+
+        // k3
+        buffer->load(state);
+        buffer->add_state_with_coef_single_thread(-1.i * dt / 2., k2);
+        _effective_hamiltonian->apply_to_state_single_thread(buffer, k3);
+
+        // k4
+        buffer->load(state);
+        buffer->add_state_with_coef_single_thread(-1.i * dt, k3);
+        _effective_hamiltonian->apply_to_state_single_thread(buffer, k4);
+
+        // store the previous state in buffer
+        buffer->load(state);
+
+        // add them together
+        state->add_state_with_coef_single_thread(-1.i * dt / 6., k1);
+        state->add_state_with_coef_single_thread(-1.i * dt / 3., k2);
+        state->add_state_with_coef_single_thread(-1.i * dt / 3., k3);
+        state->add_state_with_coef_single_thread(-1.i * dt / 6., k4);
+    }
+
+public:
+    ClsNoisyEvolution(Observable* hamiltonian,
+        std::vector<GeneralQuantumOperator*> c_ops, double time,
+        double dt = 1e-6) {
+        _hamiltonian = dynamic_cast<Observable*>(hamiltonian->copy());
+        for (auto const& op : c_ops) {
+            _c_ops.push_back(op->copy());
+            _c_ops_dagger.push_back(op->get_dagger());
+        }
+        _effective_hamiltonian = hamiltonian->copy();
+        for (size_t k = 0; k < _c_ops.size(); k++) {
+            auto cdagc = (*_c_ops_dagger[k]) * (*_c_ops[k]) * (-.5i);
+            *_effective_hamiltonian += cdagc;
+        }
+        _time = time;
+        _dt = dt;
+    };
+    ~ClsNoisyEvolution() {
+        delete _hamiltonian;
+        delete _effective_hamiltonian;
+        for (size_t k = 0; k < _c_ops.size(); k++) {
+            delete _c_ops[k];
+            delete _c_ops_dagger[k];
+        }
+    };
+
+    /**
+     * \~japanese-en 自身のゲート行列をセットする
+     *
+     * @param matrix 行列をセットする変数の参照
+     */
+    virtual void set_matrix(ComplexMatrix& matrix) const override {
+        std::cerr << "* Warning : Gate-matrix of noisy evolution cannot be "
+                     "defined. Nothing has been done."
+                  << std::endl;
+    }
+
+    /**
+     * \~japanese-en 乱数シードをセットする
+     *
+     * @param seed シード値
+     */
+    virtual void set_seed(int seed) override { _random.set_seed(seed); };
+
+    virtual QuantumGateBase* copy() const override {
+        return new ClsNoisyEvolution(_hamiltonian, _c_ops, _time, _dt);
+    }
+
+    /**
+     * \~japanese-en NoisyEvolution が使用する有効ハミルトニアンを得る
+     */
+    virtual GeneralQuantumOperator* get_effective_hamiltonian() const {
+        return _effective_hamiltonian->copy();
+    }
+
+    /**
+     * \~japanese-en collapse 時間を探すときに許す最大ループ数をセットする
+     *
+     * @param n ステップ数
+     */
+    virtual void set_find_collapse_max_steps(int n) {
+        this->_find_collapse_max_steps = n;
+    }
+
+    /**
+     * \~japanese-en 量子状態を更新する
+     *
+     * @param state 更新する量子状態
+     */
+    virtual void update_quantum_state(QuantumStateBase* state) {
+        double r = _random.uniform();
+        std::vector<double> cumulative_dist(_c_ops.size());
+        double prob_sum = 0;
+        auto k1 = state->copy();  // vector for Runge-Kutta k
+        auto k2 = state->copy();  // vector for Runge-Kutta k
+        auto k3 = state->copy();  // vector for Runge-Kutta k
+        auto k4 = state->copy();  // vector for Runge-Kutta k
+        auto buffer = state->copy();
+        double t = 0;
+
+        while (std::abs(t - _time) >
+               1e-10 * _time) {  // For machine precision error.
+            // For final time, we modify the step size to match the total
+            // execution time
+            auto dt = _dt;
+            if (t + _dt > _time) {
+                dt = _time - t;
+            }
+            while (std::abs(dt) >
+                   1e-10 * _dt) {  // dt is decreased by dt_target_norm if jump
+                                   // occurs and if jump did not occure dt will
+                                   // be set to zero.
+                // evolve the state by dt
+                _evolve_one_step(k1, k2, k3, k4, buffer, state, dt);
+                // check if the jump should occur or not
+                auto norm = state->get_squared_norm_single_thread();
+                if (norm <= r) {  // jump occured
+                    // evolve the state to the time such that norm=r
+                    auto dt_target_norm =
+                        _find_collapse(k1, k2, k3, k4, buffer, state, r, dt);
+                    if (dt_target_norm < 0) {
+                        std::cerr
+                            << "_find_collapse failed. Result is unreliable."
+                            << std::endl;
+                        return;
+                    }
+                    // get cumulative distribution
+                    prob_sum = 0.;
+                    for (size_t k = 0; k < _c_ops.size(); k++) {
+                        _c_ops[k]->apply_to_state_single_thread(state, buffer);
+                        cumulative_dist[k] =
+                            buffer->get_squared_norm_single_thread() + prob_sum;
+                        prob_sum = cumulative_dist[k];
+                    }
+
+                    // determine which collapse operator to be applied
+                    auto jump_r = _random.uniform() * prob_sum;
+                    auto ite = std::lower_bound(
+                        cumulative_dist.begin(), cumulative_dist.end(), jump_r);
+                    auto index = std::distance(cumulative_dist.begin(), ite);
+
+                    // apply the collapse operator and normalize the state
+                    _c_ops[index]->apply_to_state_single_thread(state, buffer);
+                    buffer->normalize_single_thread(
+                        buffer->get_squared_norm_single_thread());
+                    state->load(buffer);
+
+                    // update dt to be consistent with the step size
+                    t += dt_target_norm;
+                    dt -= dt_target_norm;
+
+                    // update random variable
+                    r = _random.uniform();
+                } else {  // if jump did not occur, update t to the next time
+                          // and break the loop
+                    t += dt;
+                    dt = 0.;
+                }
+            }
+        }
+
+        // normalize the state and finish
+        state->normalize_single_thread(state->get_squared_norm_single_thread());
+        delete k1;
+        delete k2;
+        delete k3;
+        delete k4;
+        delete buffer;
+    }
+};