diff --git a/src/cppsim/gate_factory.cpp b/src/cppsim/gate_factory.cpp index ccf99b0a..7b5ee4f1 100644 --- a/src/cppsim/gate_factory.cpp +++ b/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 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; diff --git a/src/cppsim/gate_factory.hpp b/src/cppsim/gate_factory.hpp index fd8056f6..b9affb8d 100644 --- a/src/cppsim/gate_factory.hpp +++ b/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 c_ops, double time, double dt = 1e-6); + } // namespace gate diff --git a/src/cppsim/gate_noisy_evolution.hpp b/src/cppsim/gate_noisy_evolution.hpp new file mode 100644 index 00000000..5df03725 --- /dev/null +++ b/src/cppsim/gate_noisy_evolution.hpp @@ -0,0 +1,340 @@ + +#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 _c_ops; // collapse operator + std::vector + _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 =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 + * が行われる。 + * + * 割線法を利用する場合、normは時間に対して広義単調減少であることが必要。 + */ + virtual double _find_collapse(QuantumStateBase* k1, QuantumStateBase* k2, + QuantumStateBase* k3, QuantumStateBase* k4, + QuantumStateBase* prev_state, QuantumStateBase* now_state, + double target_norm, double dt, bool use_secant_method = true) { + if (!use_secant_method) { + return _find_collapse_original( + k1, k2, k3, k4, prev_state, now_state, target_norm, dt); + } + auto mae_norm = prev_state->get_squared_norm_single_thread(); + auto now_norm = now_state->get_squared_norm_single_thread(); + double t_mae = 0; + double t_now = dt; + + // mae now で挟み撃ちする + int search_count = 0; + + if (std::abs(mae_norm - target_norm) < _norm_tol) { + now_state->load(prev_state); + return 0; + } + if (std::abs(now_norm - target_norm) < _norm_tol) { + return dt; + } + if (mae_norm < target_norm) { + throw std::runtime_error( + "must be prev_state.norm() >= target_norm. "); + } + if (now_norm > target_norm) { + throw std::runtime_error( + "must be now_state.norm() <= target_norm. "); + } + + QuantumStateBase* mae_state = prev_state->copy(); + double target_norm_log = std::log(target_norm); + double mae_norm_log = std::log(mae_norm); + double now_norm_log = std::log(now_norm); + QuantumStateBase* buf_state = prev_state->copy(); + QuantumStateBase* bufB_state = prev_state->copy(); + while (true) { + // we expect norm to reduce as Ae^-a*dt, so use log. + + double t_guess = 0; + if (search_count <= 20) { + // use secant method + t_guess = t_mae + (t_now - t_mae) * + (mae_norm_log - target_norm_log) / + (mae_norm_log - now_norm_log); + } else { + // use bisection method + t_guess = (t_mae + t_now) / 2; + } + + // evolve by time t_guess + buf_state->load(prev_state); + _evolve_one_step(k1, k2, k3, k4, bufB_state, buf_state, t_guess); + + double buf_norm = buf_state->get_squared_norm_single_thread(); + if (std::abs(buf_norm - target_norm) < _norm_tol) { + now_state->load(buf_state); + delete mae_state; + delete buf_state; + delete bufB_state; + + return t_guess; + } else if (buf_norm < target_norm) { + now_state->load(buf_state); + t_now = t_guess; + now_norm = now_state->get_squared_norm_single_thread(); + now_norm_log = std::log(now_norm); + } else { + mae_state->load(buf_state); + t_mae = t_guess; + mae_norm = mae_state->get_squared_norm_single_thread(); + mae_norm_log = std::log(mae_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."); + } + } + //ここには来ない + throw std::runtime_error( + "unexpectedly come to end of _find_collapse function."); + } + + /** + * \~japanese-en collapse が起こるタイミング (norm = rになるタイミング) + * を見つける。 この関数内で norm=rになるタイミングまでの evolution + * が行われる。 + * + * 割線法を利用する場合、normは時間に対して広義単調減少であることが必要。 + */ + virtual double _find_collapse_original(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; + // 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_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) { + 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_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 c_ops, double time, + double dt = 1e-6) { + _hamiltonian = dynamic_cast(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 initial_squared_norm = state->get_squared_norm(); + double r = _random.uniform(); + std::vector 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; + } + _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; + dt_target_norm = + _find_collapse(k1, k2, k3, k4, buffer, state, r, dt); + + // 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(buffer->get_squared_norm_single_thread()); + state->load(buffer); + + // update dt to be consistent with the step size + t += dt_target_norm; + + // update random variable + r = _random.uniform(); + } else { // if jump did not occur, update t to the next time + t += dt; + } + } + + // normalize the state and finish + state->normalize_single_thread( + state->get_squared_norm_single_thread() / initial_squared_norm); + delete k1; + delete k2; + delete k3; + delete k4; + delete buffer; + } +}; diff --git a/test/cppsim/test_noisyevolution.cpp b/test/cppsim/test_noisyevolution.cpp new file mode 100644 index 00000000..b550ee9c --- /dev/null +++ b/test/cppsim/test_noisyevolution.cpp @@ -0,0 +1,263 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../util/util.hpp" + +TEST(NoisyEvolutionTest, simple_check) { + // just check runtime error + UINT n = 4; + QuantumCircuit circuit(n); + Observable hamiltonian(n); + hamiltonian.add_operator(1., "Z 0 Z 1"); + GeneralQuantumOperator op(n), op2(n); + std::vector c_ops; + op.add_operator(1., "Z 0"); + op2.add_operator(1., "Z 1"); + c_ops.push_back(&op); + c_ops.push_back(&op2); + double time = 1.; + double dt = .01; + auto gate = gate::NoisyEvolution(&hamiltonian, c_ops, time, dt); + circuit.add_gate(gate); + QuantumState state(n); + circuit.update_quantum_state(&state); +} + +TEST(NoisyEvolutionTest, unitary_evolution) { + // check unitary evolution under ZZ hamiltonian + UINT n = 2; + QuantumCircuit circuit(n), circuit_ref(n); + Observable observable(n); + observable.add_operator(1, "X 0"); + + // create hamiltonian and collapse operator + Observable hamiltonian(n); + hamiltonian.add_operator(1., "Z 0 Z 1"); + GeneralQuantumOperator op(n); + std::vector c_ops; + op.add_operator(0., "Z 0"); + c_ops.push_back(&op); + int n_steps = 10; + double time = 3.14 / n_steps; // evolve by exp(-i*pi*ZZ/n_steps) + double dt = .001; + auto gate = gate::NoisyEvolution(&hamiltonian, c_ops, time, dt); + circuit.add_gate(gate); + + // reference circuit + std::vector target_index_list{0, 1}; + std::vector pauli_id_list{3, 3}; + circuit_ref.add_multi_Pauli_rotation_gate( + target_index_list, pauli_id_list, -time * 2); + + QuantumState state(n); + QuantumState state_ref(n); + gate::H(0)->update_quantum_state(&state); + gate::H(0)->update_quantum_state(&state_ref); + for (int k = 0; k < n_steps; k++) { + circuit.update_quantum_state(&state); + circuit_ref.update_quantum_state(&state_ref); + auto exp = observable.get_expectation_value(&state); + auto exp_ref = observable.get_expectation_value(&state_ref); + ASSERT_NEAR(exp.real(), exp_ref.real(), 1e-6); + } +} + +TEST(NoisyEvolutionTest, error_scaling) { + // check runge kutta error for unitary evolution under ZZ hamiltonian + int n_dt = 4; + double dt_start = 0.1; // dt will be decreased by 1/dt_scale + double dt_scale = 1.189207115; // this choice should make the error 1/2 for + // each iteration + double time = 1.; + UINT n = 2; + Observable observable(n); + observable.add_operator(1, "X 0"); + + // create hamiltonian and collapse operator + Observable hamiltonian(n); + hamiltonian.add_operator(1., "Z 0 Z 1"); + GeneralQuantumOperator op(n); + std::vector c_ops; + op.add_operator(0., "Z 0"); + c_ops.push_back(&op); + // reference circuit + QuantumCircuit circuit_ref(n); + std::vector target_index_list{0, 1}; + std::vector pauli_id_list{3, 3}; + circuit_ref.add_multi_Pauli_rotation_gate( + target_index_list, pauli_id_list, -time * 2); + + double prev_error = 1; + auto dt = dt_start; + QuantumState state(n); + QuantumState state_ref(n); + for (int k = 0; k < n_dt; k++) { + state.set_zero_state(); + state_ref.set_zero_state(); + gate::H(0)->update_quantum_state(&state); + gate::H(0)->update_quantum_state(&state_ref); + QuantumCircuit circuit(n); + auto gate = gate::NoisyEvolution(&hamiltonian, c_ops, time, dt); + circuit.add_gate(gate); + circuit.update_quantum_state(&state); + circuit_ref.update_quantum_state(&state_ref); + state.add_state_with_coef(-1., &state_ref); + if (k != 0) { + ASSERT_NEAR(std::sqrt(state.get_squared_norm()) / prev_error, + 1. / std::pow(dt_scale, 4.), .1); + } + prev_error = std::sqrt(state.get_squared_norm()); + dt /= dt_scale; + } +} + +TEST(NoisyEvolutionTest, EffectiveHamiltonian) { + // 2 qubit dephasing dynamics with ZZ interaction + double time = 1.; + double dt = 0.01; + UINT n = 2; + // create hamiltonian and collapse operator + Observable hamiltonian(n); + hamiltonian.add_operator(1., "Z 0 Z 1"); + GeneralQuantumOperator op(n); + std::vector c_ops; + op.add_operator(1., "Z 0"); + c_ops.push_back(&op); + GeneralQuantumOperator op2(n); + op2.add_operator(1., "Z 1"); + c_ops.push_back(&op2); + + // noisy evolution gate + auto gate = dynamic_cast( + gate::NoisyEvolution(&hamiltonian, c_ops, time, dt)); + ASSERT_EQ(gate->get_effective_hamiltonian()->to_string(), + "(1,0) Z 0 Z 1 + (0,-1) I"); +} + +TEST(NoisyEvolutionTest, dephasing) { + // 2 qubit dephasing dynamics with ZZ interaction + double time = 2.; + double dt = 0.1; + double decay_rate = 0.2; + double hamiltonian_energy = 0.2; + double ref = 0.5936940289967207; // generated by qutip + double ref_withoutnoise = 0.696706573268861; // generated by qutip (needed + // for variance calculation) + UINT n = 2; + UINT n_samples = 1000; + Observable observable(n); + observable.add_operator(1, "X 0"); + // create hamiltonian and collapse operator + Observable hamiltonian(n); + hamiltonian.add_operator(hamiltonian_energy, "Z 0 Z 1"); + GeneralQuantumOperator op(n); + std::vector c_ops; + op.add_operator(decay_rate, "Z 0"); + c_ops.push_back(&op); + GeneralQuantumOperator op2(n); + op2.add_operator(decay_rate, "Z 1"); + c_ops.push_back(&op2); + + QuantumState state(n), init_state(n); + QuantumCircuit circuit(n); + gate::H(0)->update_quantum_state(&init_state); + gate::H(1)->update_quantum_state(&init_state); + circuit.add_gate(gate::NoisyEvolution(&hamiltonian, c_ops, time, dt)); + double exp = 0.; + for (int k = 0; k < n_samples; k++) { + state.load(&init_state); + circuit.update_quantum_state(&state); + exp += observable.get_expectation_value(&state).real() / n_samples; + } + // intrinsic variance is (

_withoutnoise^2-

_withnoise^2). + // take 5-sigma for assertion. Correct code should violate this assertion by + // probabilty of only 3e-4 %. + auto fivesigma = 5 * sqrt(ref_withoutnoise * ref_withoutnoise - ref * ref) / + sqrt(n_samples); + ASSERT_NEAR(exp, ref, fivesigma); +} + +TEST(NoisyEvolutionTest, T1T2) { + // 2 qubit dephasing dynamics with ZZ interaction + double time = 2.; + double dt = 0.1; + double decay_rate_z = 0.2; + double decay_rate_p = 0.6; + double decay_rate_m = 0.1; + double hamiltonian_energy = 1.; + 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 + Observable hamiltonian(n); + hamiltonian.add_operator(hamiltonian_energy, "Z 0 Z 1"); + std::vector c_ops; + for (int k = 0; k < 6; k++) c_ops.push_back(new GeneralQuantumOperator(n)); + c_ops[0]->add_operator(decay_rate_z, "Z 0"); + c_ops[1]->add_operator(decay_rate_z, "Z 1"); + c_ops[2]->add_operator(decay_rate_p / 2, "X 0"); + c_ops[2]->add_operator(decay_rate_p / 2 * 1.i, "Y 0"); + c_ops[3]->add_operator(decay_rate_p / 2, "X 1"); + c_ops[3]->add_operator(decay_rate_p / 2 * 1.i, "Y 1"); + c_ops[4]->add_operator(decay_rate_m / 2, "X 0"); + c_ops[4]->add_operator(-decay_rate_m / 2 * 1.i, "Y 0"); + c_ops[5]->add_operator(decay_rate_m / 2, "X 1"); + c_ops[5]->add_operator(-decay_rate_m / 2 * 1.i, "Y 1"); + + QuantumState state(n); + QuantumCircuit circuit(n); + circuit.add_gate(gate::NoisyEvolution(&hamiltonian, c_ops, time, dt)); + double exp = 0.; + for (int k = 0; k < n_samples; k++) { + state.set_zero_state(); + 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; + ASSERT_NEAR(exp, ref, .1); +} + +TEST(NoisyEvolutionTest, check_inf_occurence) { + // test case for checking inf + double time = 1.; + double dt = 0.1; + double decay_rate = 1.; + double hamiltonian_energy = 1.; + UINT n = 2; + UINT n_samples = 100; + Observable observable(n); + observable.add_operator(1, "X 0"); + // create hamiltonian and collapse operator + Observable hamiltonian(n); + hamiltonian.add_operator(hamiltonian_energy, "Z 0 Z 1"); + GeneralQuantumOperator op(n); + std::vector c_ops; + op.add_operator(decay_rate, "Z 0"); + c_ops.push_back(&op); + GeneralQuantumOperator op2(n); + op2.add_operator(decay_rate, "Z 1"); + c_ops.push_back(&op2); + + QuantumState state(n); + QuantumCircuit circuit(n); + circuit.add_gate(gate::NoisyEvolution(&hamiltonian, c_ops, time, dt)); + double exp = 0.; + for (int k = 0; k < n_samples; k++) { + circuit.update_quantum_state(&state); + ASSERT_FALSE( + std::isinf(observable.get_expectation_value(&state).real())); + } +}