Program Listing for File gate_noisy_evolution.hpp¶
↰ Return to documentation for file (/home/docs/checkouts/readthedocs.org/user_builds/qulacs-rtd/checkouts/v0.6.7/src/cppsim/gate_noisy_evolution.hpp
)
#pragma once
#include <Eigen/Dense>
#include <csim/stat_ops.hpp>
#include <csim/update_ops.hpp>
#include <csim/update_ops_dm.hpp>
#include <cstring>
#include <fstream>
#include <numeric>
#include "exception.hpp"
#include "gate.hpp"
#include "gate_merge.hpp"
#include "gate_named_pauli.hpp"
#include "general_quantum_operator.hpp"
#include "observable.hpp"
#include "pauli_operator.hpp"
#include "state.hpp"
#include "type.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-6; // accuracy in solving <psi|psi>=r
int _find_collapse_max_steps =
200; // maximum number of steps for while loop in _find_collapse
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);
virtual double _find_collapse_original(QuantumStateBase* k1,
QuantumStateBase* k2, QuantumStateBase* k3, QuantumStateBase* k4,
QuantumStateBase* prev_state, QuantumStateBase* now_state,
double target_norm, double dt);
virtual void _evolve_one_step(QuantumStateBase* k1, QuantumStateBase* k2,
QuantumStateBase* k3, QuantumStateBase* k4, QuantumStateBase* buffer,
QuantumStateBase* state, double dt);
public:
ClsNoisyEvolution(Observable* hamiltonian,
std::vector<GeneralQuantumOperator*> c_ops, double time,
double dt = 1e-6);
~ClsNoisyEvolution();
virtual void set_matrix(ComplexMatrix&) const override {
throw NotImplementedException(
"Error: "
"ClsNoisyEvolution::set_matrix(ComplexMatrix&): Gate-matrix of "
"noisy evolution cannot be defined.");
}
virtual void set_seed(int seed) override { _random.set_seed(seed); };
virtual ClsNoisyEvolution* copy() const override {
return new ClsNoisyEvolution(_hamiltonian, _c_ops, _time, _dt);
}
virtual GeneralQuantumOperator* get_effective_hamiltonian() const {
return _effective_hamiltonian->copy();
}
virtual void set_find_collapse_max_steps(int n) {
this->_find_collapse_max_steps = n;
}
virtual void update_quantum_state(QuantumStateBase* state) override;
virtual boost::property_tree::ptree to_ptree() const override;
};
/*
元のバージョンで、時間発展をルンゲクッタ法を使っていたところを、 行列を抜き出して対角化で求めます。
性質上、
・ hamiltonianやc_opsに出現するビットの種類数が多いと遅い、 多分3~4ビットあたりが境目?
・ なので、操作するビットがたくさんある場合、それが独立なビット集合に分けられる場合、 集合ごとにこのクラスを使ってください
・ ゲートを作るときに重い操作を行うので、同じゲートに対して何回も演算すると速い
*/
class ClsNoisyEvolution_fast : 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 _norm_tol = 1e-4; // accuracy in solving <psi|psi>=r
// rが乱数なら精度不要なので1e-4にした
int _find_collapse_max_steps = 100;
ComplexVector eigenvalue_mto;
QuantumGateMatrix* eigenMatrixGate = nullptr;
QuantumGateMatrix* eigenMatrixRevGate = nullptr;
bool gate_none = false;
virtual double _find_collapse(QuantumStateBase* prev_state,
QuantumStateBase* now_state, double target_norm, double t_step);
virtual void _evolve_one_step(QuantumStateBase* state, double t_step);
public:
ClsNoisyEvolution_fast(Observable* hamiltonian,
std::vector<GeneralQuantumOperator*> c_ops, double time);
~ClsNoisyEvolution_fast() {
delete _hamiltonian;
delete _effective_hamiltonian;
for (size_t k = 0; k < _c_ops.size(); k++) {
delete _c_ops[k];
delete _c_ops_dagger[k];
}
delete eigenMatrixGate;
delete eigenMatrixRevGate;
};
virtual void set_matrix(ComplexMatrix&) const override {
throw NotImplementedException(
"Error: "
"ClsNoisyEvolution::set_matrix(ComplexMatrix&): Gate-matrix of "
"noisy evolution cannot be defined.");
}
virtual void set_seed(int seed) override { _random.set_seed(seed); };
virtual ClsNoisyEvolution_fast* copy() const override {
return new ClsNoisyEvolution_fast(_hamiltonian, _c_ops, _time);
}
virtual GeneralQuantumOperator* get_effective_hamiltonian() const {
return _effective_hamiltonian->copy();
}
virtual void set_find_collapse_max_steps(int n) {
this->_find_collapse_max_steps = n;
}
virtual void change_time(double time) { _time = time; }
virtual void update_quantum_state(QuantumStateBase* state) override;
virtual boost::property_tree::ptree to_ptree() const override;
};
// noisyEvolution_auto
// c_opsなどの情報から、いろんなゲートが混在していても内部で自動的にfastの組に分けられ、最適化される
class ClsNoisyEvolution_auto : public QuantumGateBase {
private:
std::vector<ClsNoisyEvolution_fast*> gates;
std::vector<int> bit_number;
void bit_merge(int a, int b) {
if (a == b) {
return;
}
int n = bit_number.size();
for (int i = 0; i < n; i++) {
if (bit_number[i] == b) {
bit_number[i] = a;
}
}
}
public:
ClsNoisyEvolution_auto() {}
ClsNoisyEvolution_auto(Observable* hamiltonian,
std::vector<GeneralQuantumOperator*> c_ops, double time) {
// hamiltonianをunion-findしつつ、使われているかを確認
// 量少ないし、ufじゃなくて普通に実装します
// bit_numberを決めて、実際に独立なゲートに振り分ける
UINT n_qubit = hamiltonian->get_qubit_count();
bit_number = std::vector<int>(n_qubit, -1);
for (auto pauli : hamiltonian->get_terms()) {
for (int bit : pauli->get_index_list()) {
bit_number[bit] = bit;
}
}
for (size_t k = 0; k < c_ops.size(); k++) {
for (auto pauli : c_ops[k]->get_terms()) {
for (int bit : pauli->get_index_list()) {
bit_number[bit] = bit;
}
}
}
// 使ってないビットは-1にした
for (auto pauli : hamiltonian->get_terms()) {
int uf_a = -1;
for (int bit : pauli->get_index_list()) {
if (uf_a == -1) {
uf_a = bit;
} else {
bit_merge(uf_a, bit);
}
}
}
for (size_t k = 0; k < c_ops.size(); k++) {
int uf_a = -1;
for (auto pauli : c_ops[k]->get_terms()) {
for (int bit : pauli->get_index_list()) {
if (uf_a == -1) {
uf_a = bit;
} else {
bit_merge(uf_a, bit);
}
}
}
}
std::vector<bool> aru(n_qubit);
for (UINT i = 0; i < n_qubit; i++) {
if (bit_number[i] == -1) {
continue;
}
aru[bit_number[i]] = true;
}
std::vector<UINT> taiou(n_qubit);
UINT kaz = 0;
for (UINT i = 0; i < n_qubit; i++) {
if (aru[i]) {
taiou[i] = kaz;
kaz++;
}
}
for (UINT i = 0; i < n_qubit; i++) {
bit_number[i] = taiou[bit_number[i]];
}
std::vector<Observable*> hamiltonians(kaz);
for (UINT i = 0; i < kaz; i++) {
hamiltonians[i] = new Observable(n_qubit);
}
std::vector<std::vector<GeneralQuantumOperator*>> c_opss(kaz);
for (auto pauli : hamiltonian->get_terms()) {
int uf_a = -1;
for (int bit : pauli->get_index_list()) {
if (uf_a == -1) {
uf_a = bit_number[bit];
} else if (uf_a != bit_number[bit]) {
throw std::runtime_error("Error: bug");
}
}
(*hamiltonians[uf_a]) += (*pauli);
}
for (size_t k = 0; k < c_ops.size(); k++) {
int uf_a = -1;
for (auto pauli : c_ops[k]->get_terms()) {
for (int bit : pauli->get_index_list()) {
if (uf_a == -1) {
uf_a = bit_number[bit];
} else if (uf_a != bit_number[bit]) {
throw std::runtime_error("Error: bug");
}
}
}
c_opss[uf_a].push_back(c_ops[k]);
}
gates.resize(kaz);
for (UINT i = 0; i < kaz; i++) {
gates[i] =
new ClsNoisyEvolution_fast(hamiltonians[i], c_opss[i], time);
}
}
~ClsNoisyEvolution_auto() {
for (auto it : gates) {
delete it;
}
}
virtual void set_matrix(ComplexMatrix&) const override {
throw NotImplementedException(
"Error: "
"ClsNoisyEvolution::set_matrix(ComplexMatrix&): Gate-matrix of "
"noisy evolution cannot be defined.");
}
virtual void set_seed(int seed) override {
for (auto it : gates) {
it->set_seed(seed);
}
};
virtual void update_quantum_state(QuantumStateBase* state) override {
for (auto gate : gates) {
gate->update_quantum_state(state);
}
}
virtual QuantumGateBase* copy() const override {
auto tar = new ClsNoisyEvolution_auto();
tar->gates = std::vector<ClsNoisyEvolution_fast*>(this->gates.size());
for (UINT i = 0; i < this->gates.size(); i++) {
tar->gates[i] = this->gates[i]->copy();
}
tar->bit_number = this->bit_number;
return tar;
}
};