リンドブラッド方程式のモンテカルロソルバ(マルチスレッド版)

pysrc/qulacs/__init__.pyi

@@ -99,6 +99,10 @@ class GeneralQuantumOperator():
         """
     @typing.overload
     def add_operator(self, pauli_operator: PauliOperator) -> None: ...
+    def apply_to_state(self, work_state: QuantumStateBase, state_to_be_multiplied: QuantumStateBase, dst_state: QuantumStateBase) -> None: 
+        """
+        Apply observable to `state_to_be_multiplied`. The result is stored into `dst_state`.
+        """
     def copy(self) -> GeneralQuantumOperator: 
         """
         Create copied instance of General Quantum operator class
@@ -701,6 +705,10 @@ class QuantumState(QuantumStateBase):
         """
         Create copied instance
         """
+    def get_amplitude(self, index: int) -> complex: 
+        """
+        Get state vector
+        """
     def get_classical_value(self, index: int) -> int: 
         """
         Get classical 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"
@@ -355,6 +356,11 @@ QuantumGateBase* Measurement(
     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();
     char* sbuf;

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,253 @@
+
+#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();
+        auto prev_norm = prev_state->get_squared_norm();
+        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;
+            // evolve 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();
+
+            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) {
+                throw std::runtime_error(
+                    "Failed to find the exact jump time. Try with "
+                    "smaller dt.");
+            }
+        }
+        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(state, k1);
+
+        // k2
+        buffer->load(state);
+        buffer->add_state_with_coef(-1.i * dt / 2., k1);
+        _effective_hamiltonian->apply_to_state(buffer, k2);
+
+        // k3
+        buffer->load(state);
+        buffer->add_state_with_coef(-1.i * dt / 2., k2);
+        _effective_hamiltonian->apply_to_state(buffer, k3);
+
+        // k4
+        buffer->load(state);
+        buffer->add_state_with_coef(-1.i * dt, k3);
+        _effective_hamiltonian->apply_to_state(buffer, k4);
+
+        // store the previous state in buffer
+        buffer->load(state);
+
+        // add them together
+        state->add_state_with_coef(-1.i * dt / 6., k1);
+        state->add_state_with_coef(-1.i * dt / 3., k2);
+        state->add_state_with_coef(-1.i * dt / 3., k3);
+        state->add_state_with_coef(-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::stringstream error_message_stream;
+        error_message_stream
+            << "* Warning : Gate-matrix of noisy evolution cannot be "
+               "defined. Nothing has been done.";
+        throw std::invalid_argument(error_message_stream.str());
+    }
+
+    /**
+     * \~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();
+                if (norm <= r) {  // jump occured
+                    // evolve the state to the time such that norm=r
+                    double dt_target_norm;
+                    try {
+                        dt_target_norm = _find_collapse(
+                            k1, k2, k3, k4, buffer, state, r, dt);
+                    } catch (std::runtime_error& e) {
+                        throw std::runtime_error(
+                            "_find_collapse failed. Result is "
+                            "unreliable.");
+                        return;
+                    }
+
+                    // get cumulative distribution
+                    prob_sum = 0.;
+                    for (size_t k = 0; k < _c_ops.size(); k++) {
+                        _c_ops[k]->apply_to_state(state, buffer);
+                        cumulative_dist[k] =
+                            buffer->get_squared_norm() + 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(state, buffer);
+                    buffer->normalize(buffer->get_squared_norm());
+                    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(state->get_squared_norm());
+        delete k1;
+        delete k2;
+        delete k3;
+        delete k4;
+        delete buffer;
+    }
+};

src/cppsim/general_quantum_operator.cpp

@@ -429,7 +429,7 @@ GeneralQuantumOperator& GeneralQuantumOperator::operator+=(
     const GeneralQuantumOperator& target) {
     ITYPE i, j;
     auto terms = target.get_terms();
-#pragma omp parallel for
+    // #pragma omp parallel for
     for (i = 0; i < _operator_list.size(); i++) {
         auto pauli_operator = _operator_list[i];
         for (j = 0; j < terms.size(); j++) {
@@ -482,7 +482,7 @@ GeneralQuantumOperator& GeneralQuantumOperator::operator+=(
     const PauliOperator& target) {
     bool flag = true;
     ITYPE i;
-#pragma omp parallel for
+    // #pragma omp parallel for
     for (i = 0; i < _operator_list.size(); i++) {
         auto pauli_operator = _operator_list[i];
         auto pauli_x = pauli_operator->get_x_bits();
@@ -526,7 +526,7 @@ GeneralQuantumOperator& GeneralQuantumOperator::operator-=(
     const GeneralQuantumOperator& target) {
     ITYPE i, j;
     auto terms = target.get_terms();
-#pragma omp parallel for
+    // #pragma omp parallel for
     for (i = 0; i < _operator_list.size(); i++) {
         auto pauli_operator = _operator_list[i];
         for (j = 0; j < terms.size(); j++) {
@@ -634,7 +634,7 @@ GeneralQuantumOperator& GeneralQuantumOperator::operator*=(
     auto terms = copy->get_terms();
     auto target_terms = target.get_terms();
     ITYPE i, j;
-#pragma omp parallel for
+    //#pragma omp parallel for
     for (i = 0; i < terms.size(); i++) {
         auto pauli_operator = terms[i];
         for (j = 0; j < target_terms.size(); j++) {
@@ -653,7 +653,7 @@ GeneralQuantumOperator& GeneralQuantumOperator::operator*=(
     _operator_list.clear();
     ITYPE i;
     auto terms = copy->get_terms();
-#pragma omp parallel for
+    //#pragma omp parallel for
     for (i = 0; i < terms.size(); i++) {
         auto pauli_operator = terms[i];
         PauliOperator* product = new PauliOperator;
@@ -665,7 +665,7 @@ GeneralQuantumOperator& GeneralQuantumOperator::operator*=(
 
 GeneralQuantumOperator& GeneralQuantumOperator::operator*=(CPPCTYPE target) {
     ITYPE i;
-#pragma omp parallel for
+    //#pragma omp parallel for
     for (i = 0; i < _operator_list.size(); i++) {
         *_operator_list[i] *= target;
     }