engines.h

Go to the documentation of this file.

/*
 * This file is part of Quantum++.
 *
 * MIT License
 *
 * Copyright (c) 2013 - 2020 Vlad Gheorghiu.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#ifndef CLASSES_CIRCUITS_ENGINES_H_
#define CLASSES_CIRCUITS_ENGINES_H_

namespace qpp {
class QEngine : public IDisplay, public IJSON {
  protected:
    const QCircuit* qc_; 

    struct state_ {
        const QCircuit* qc_;          
        ket psi_{};                   
        std::vector<double> probs_{}; 
        std::vector<idx> dits_{};     
        std::vector<idx> subsys_{}; 

        explicit state_(const QCircuit* qc) : qc_{qc} {
            // EXECEPTION CHECKS

            if (qc->get_nq() == 0)
                throw exception::ZeroSize("qpp::QEngine::state_::reset()");
            // END EXCEPTION CHECKS
            reset();
        }

        // silence -Weffc++ class has pointer data members
        state_(const state_&) = default;

        // silence -Weffc++ class has pointer data members
        state_& operator=(const state_&) = default;

        void reset() {
            psi_ = States::get_instance().zero(qc_->get_nq(), qc_->get_d());
            probs_ = std::vector<double>(qc_->get_nc(), 0);
            dits_ = std::vector<idx>(qc_->get_nc(), 0);
            subsys_ = std::vector<idx>(qc_->get_nq(), 0);
            std::iota(std::begin(subsys_), std::end(subsys_), 0);
        }
    } st_; 
    std::map<std::string, idx, internal::EqualSameSizeStringDits>
        stats_; 

    void set_measured_(idx i) {
        if (get_measured(i))
            throw exception::QuditAlreadyMeasured(
                "qpp::QEngine::set_measured_()");
        st_.subsys_[i] = static_cast<idx>(-1); // set qudit i to measured state
        for (idx m = i; m < qc_->get_nq(); ++m) {
            if (!get_measured(m)) {
                --st_.subsys_[m];
            }
        }
    }

    // giving a vector of non-measured qudits, get their relative position wrt
    // the measured qudits
    std::vector<idx> get_relative_pos_(std::vector<idx> v) {
        idx vsize = v.size();
        for (idx i = 0; i < vsize; ++i) {
            if (get_measured(v[i]))
                throw exception::QuditAlreadyMeasured(
                    "qpp::QEngine::get_relative_pos_()");
            v[i] = st_.subsys_[v[i]];
        }
        return v;
    }

  public:
    explicit QEngine(const QCircuit& qc)
        : qc_{std::addressof(qc)}, st_{qc_}, stats_{} {}

    // silence -Weffc++ class has pointer data members
    QEngine(const QEngine&) = default;

    // silence -Weffc++ class has pointer data members
    QEngine& operator=(const QEngine&) = default;

    QEngine(QCircuit&&) = delete;

    virtual ~QEngine() = default;

    // getters
    ket get_psi() const { return st_.psi_; }

    std::vector<idx> get_dits() const { return st_.dits_; }

    idx get_dit(idx i) const {
        if (i >= qc_->get_nc())
            throw exception::OutOfRange("qpp::QEngine::get_dit()");

        return st_.dits_[i];
    }

    std::vector<double> get_probs() const { return st_.probs_; }

    bool get_measured(idx i) const {
        return st_.subsys_[i] == static_cast<idx>(-1);
    }

    std::vector<idx> get_measured() const {
        std::vector<idx> result;
        for (idx i = 0; i < qc_->get_nq(); ++i)
            if (get_measured(i))
                result.emplace_back(i);

        return result;
    }

    std::vector<idx> get_non_measured() const {
        std::vector<idx> result;
        for (idx i = 0; i < qc_->get_nq(); ++i)
            if (!get_measured(i))
                result.emplace_back(i);

        return result;
    }

    const QCircuit& get_circuit() const& noexcept { return *qc_; }

    QCircuit get_circuit() const&& noexcept { return *qc_; }

    std::map<std::string, idx, internal::EqualSameSizeStringDits>
    get_stats() const {
        return stats_;
    }
    // end getters

    // setters
    QEngine& set_dit(idx i, idx value) {
        if (i >= qc_->get_nc())
            throw exception::OutOfRange("qpp::QEngine::set_dit()");
        st_.dits_[i] = value;

        return *this;
    }

    QEngine& set_psi(const ket& psi) {
        // EXCEPTION CHECKS

        idx n = get_non_measured().size();
        idx D = static_cast<idx>(std::llround(std::pow(qc_->get_d(), n)));
        if (static_cast<idx>(psi.rows()) != D)
            throw exception::DimsNotEqual("qpp::QEngine::set_psi()");
        // END EXCEPTION CHECKS

        st_.psi_ = psi;

        return *this;
    }
    // end setters

    QEngine& reset_stats() {
        stats_ = {};

        return *this;
    }

    QEngine& reset() {
        st_.reset();

        return *this;
    }

    virtual QEngine& execute(const QCircuit::iterator::value_type& elem) {
        // EXCEPTION CHECKS

        // iterator must point to the same quantum circuit description
        if (elem.value_type_qc_ != qc_)
            throw exception::InvalidIterator("qpp::QEngine::execute()");
        // the rest of exceptions are caught by the iterator::operator*()
        // END EXCEPTION CHECKS

        auto h_tbl = qc_->get_cmat_hash_tbl_();
        idx d = qc_->get_d();

        // gate step
        if (elem.type_ == QCircuit::StepType::GATE) {
            auto gates = qc_->get_gates_();

            idx q_ip =
                std::distance(std::begin(qc_->get_gates_()), elem.gates_ip_);

            std::vector<idx> ctrl_rel_pos;
            std::vector<idx> target_rel_pos =
                get_relative_pos_(gates[q_ip].target_);

            switch (gates[q_ip].gate_type_) {
                case QCircuit::GateType::NONE:
                    break;
                case QCircuit::GateType::SINGLE:
                case QCircuit::GateType::TWO:
                case QCircuit::GateType::THREE:
                case QCircuit::GateType::CUSTOM:
                    st_.psi_ = apply(st_.psi_, h_tbl[gates[q_ip].gate_hash_],
                                     target_rel_pos, d);
                    break;
                case QCircuit::GateType::FAN:
                    for (idx m = 0; m < gates[q_ip].target_.size(); ++m)
                        st_.psi_ =
                            apply(st_.psi_, h_tbl[gates[q_ip].gate_hash_],
                                  {target_rel_pos[m]}, d);
                    break;
                case QCircuit::GateType::SINGLE_CTRL_SINGLE_TARGET:
                case QCircuit::GateType::SINGLE_CTRL_MULTIPLE_TARGET:
                case QCircuit::GateType::MULTIPLE_CTRL_SINGLE_TARGET:
                case QCircuit::GateType::MULTIPLE_CTRL_MULTIPLE_TARGET:
                case QCircuit::GateType::CUSTOM_CTRL:
                    ctrl_rel_pos = get_relative_pos_(gates[q_ip].ctrl_);
                    st_.psi_ = applyCTRL(
                        st_.psi_, h_tbl[gates[q_ip].gate_hash_], ctrl_rel_pos,
                        target_rel_pos, d, gates[q_ip].shift_);
                    break;
                case QCircuit::GateType::SINGLE_cCTRL_SINGLE_TARGET:
                case QCircuit::GateType::SINGLE_cCTRL_MULTIPLE_TARGET:
                case QCircuit::GateType::MULTIPLE_cCTRL_SINGLE_TARGET:
                case QCircuit::GateType::MULTIPLE_cCTRL_MULTIPLE_TARGET:
                case QCircuit::GateType::CUSTOM_cCTRL:
                    if (st_.dits_.empty()) {
                        st_.psi_ =
                            apply(st_.psi_, h_tbl[gates[q_ip].gate_hash_],
                                  target_rel_pos, d);
                    } else {
                        bool should_apply = true;
                        idx first_dit;
                        // we have a shift
                        if (!gates[q_ip].shift_.empty()) {
                            first_dit = (st_.dits_[(gates[q_ip].ctrl_)[0]] +
                                         gates[q_ip].shift_[0]) %
                                        d;
                            for (idx m = 1; m < gates[q_ip].ctrl_.size(); ++m) {
                                if ((st_.dits_[(gates[q_ip].ctrl_)[m]] +
                                     gates[q_ip].shift_[m]) %
                                        d !=
                                    first_dit) {
                                    should_apply = false;
                                    break;
                                }
                            }
                        }
                        // no shift
                        else {
                            first_dit = st_.dits_[(gates[q_ip].ctrl_)[0]];
                            for (idx m = 1; m < gates[q_ip].ctrl_.size(); ++m) {
                                if (st_.dits_[(gates[q_ip].ctrl_)[m]] !=
                                    first_dit) {
                                    should_apply = false;
                                    break;
                                }
                            }
                        }
                        if (should_apply) {
                            st_.psi_ = apply(
                                st_.psi_,
                                powm(h_tbl[gates[q_ip].gate_hash_], first_dit),
                                target_rel_pos, d);
                        }
                    }
                    break;
            } // end switch on gate type
        }     // end if gate step

        // measurement step
        else if (elem.type_ == QCircuit::StepType::MEASUREMENT) {
            auto measurements = qc_->get_measurements_();
            idx m_ip = std::distance(std::begin(qc_->get_measurements_()),
                                     elem.measurements_ip_);

            std::vector<idx> target_rel_pos =
                get_relative_pos_(measurements[m_ip].target_);

            std::vector<idx> resZ;
            double probZ;

            idx mres = 0;
            std::vector<double> probs;
            std::vector<cmat> states;

            switch (measurements[m_ip].measurement_type_) {
                case QCircuit::MeasureType::NONE:
                    break;
                case QCircuit::MeasureType::MEASURE_Z:
                    std::tie(resZ, probZ, st_.psi_) =
                        measure_seq(st_.psi_, target_rel_pos, d);
                    st_.dits_[measurements[m_ip].c_reg_] = resZ[0];
                    st_.probs_[measurements[m_ip].c_reg_] = probZ;
                    set_measured_(measurements[m_ip].target_[0]);
                    break;
                case QCircuit::MeasureType::MEASURE_Z_MANY:
                    std::tie(resZ, probZ, st_.psi_) =
                        measure_seq(st_.psi_, target_rel_pos, d);
                    st_.dits_[measurements[m_ip].c_reg_] = multiidx2n(
                        resZ, std::vector<idx>(target_rel_pos.size(), d));
                    st_.probs_[measurements[m_ip].c_reg_] = probZ;
                    for (auto&& elem : measurements[m_ip].target_)
                        set_measured_(elem);
                    break;
                case QCircuit::MeasureType::MEASURE_V:
                    std::tie(mres, probs, states) = measure(
                        st_.psi_, h_tbl[measurements[m_ip].mats_hash_[0]],
                        target_rel_pos, d);
                    st_.psi_ = states[mres];
                    st_.dits_[measurements[m_ip].c_reg_] = mres;
                    st_.probs_[measurements[m_ip].c_reg_] = probs[mres];
                    set_measured_(measurements[m_ip].target_[0]);
                    break;
                case QCircuit::MeasureType::MEASURE_V_MANY:
                    std::tie(mres, probs, states) = measure(
                        st_.psi_, h_tbl[measurements[m_ip].mats_hash_[0]],
                        target_rel_pos, d);
                    st_.psi_ = states[mres];
                    st_.dits_[measurements[m_ip].c_reg_] = mres;
                    st_.probs_[measurements[m_ip].c_reg_] = probs[mres];
                    for (auto&& elem : measurements[m_ip].target_)
                        set_measured_(elem);
                    break;
                case QCircuit::MeasureType::MEASURE_Z_ND:
                    std::tie(resZ, probZ, st_.psi_) =
                        measure_seq(st_.psi_, target_rel_pos, d, false);
                    st_.dits_[measurements[m_ip].c_reg_] = resZ[0];
                    st_.probs_[measurements[m_ip].c_reg_] = probZ;
                    break;
                case QCircuit::MeasureType::MEASURE_Z_MANY_ND:
                    std::tie(resZ, probZ, st_.psi_) =
                        measure_seq(st_.psi_, target_rel_pos, d, false);
                    st_.dits_[measurements[m_ip].c_reg_] = multiidx2n(
                        resZ, std::vector<idx>(target_rel_pos.size(), d));
                    st_.probs_[measurements[m_ip].c_reg_] = probZ;
                    break;
                case QCircuit::MeasureType::MEASURE_V_ND:
                case QCircuit::MeasureType::MEASURE_V_MANY_ND:
                    std::tie(mres, probs, states) = measure(
                        st_.psi_, h_tbl[measurements[m_ip].mats_hash_[0]],
                        target_rel_pos, d, false);
                    st_.psi_ = states[mres];
                    st_.dits_[measurements[m_ip].c_reg_] = mres;
                    st_.probs_[measurements[m_ip].c_reg_] = probs[mres];
                    break;
                case QCircuit::MeasureType::RESET:
                case QCircuit::MeasureType::RESET_MANY:
                    st_.psi_ = qpp::reset(st_.psi_, target_rel_pos, d);
                    break;
                case QCircuit::MeasureType::DISCARD:
                    std::tie(std::ignore, std::ignore, st_.psi_) =
                        measure_seq(st_.psi_, target_rel_pos, d);
                    set_measured_(measurements[m_ip].target_[0]);
                    break;
                case QCircuit::MeasureType::DISCARD_MANY:
                    std::tie(std::ignore, std::ignore, st_.psi_) =
                        measure_seq(st_.psi_, target_rel_pos, d);
                    for (auto&& elem : measurements[m_ip].target_)
                        set_measured_(elem);
                    break;
            } // end switch on measurement type
        }     // end else if measurement step

        // no-op
        else if (elem.type_ == QCircuit::StepType::NOP) {
        }

        // otherwise
        else {
        }

        return *this;
    }

    QEngine& execute(const QCircuit::iterator& it) { return execute(*it); }

    QEngine& execute(idx reps = 1, bool clear_stats = true) {
        auto initial_engine_state = st_; // saves the engine entry state

        if (clear_stats)
            reset_stats();

        for (idx i = 0; i < reps; ++i) {
            // sets the state of the engine to the entry state
            st_ = initial_engine_state;

            for (auto&& elem : *qc_)
                (void) execute(elem);

            // we measured at least one qudit
            if (qc_->get_measurement_count() > 0) {
                std::vector<idx> m_res = get_dits();

                std::stringstream ss;
                ss << disp(m_res, " ", "", "");
                ++stats_[ss.str()];
            }
        }

        return *this;
    }

    std::string to_JSON(bool enclosed_in_curly_brackets = true) const override {
        std::string result;

        if (enclosed_in_curly_brackets)
            result += "{";

        std::ostringstream ss;
        ss << disp(get_measured(), ", ");
        result += "\"measured/discarded (destructive)\" : " + ss.str() + ", ";

        ss.str("");
        ss.clear();
        ss << disp(get_non_measured(), ", ");
        result += "\"non-measured/non-discarded\" : " + ss.str();

        ss.str("");
        ss.clear();
        result += ", \"last probs\": ";
        ss << disp(get_probs(), ", ");
        result += ss.str();

        ss.str("");
        ss.clear();
        result += ", \"last dits\": ";
        ss << disp(get_dits(), ", ");
        result += ss.str();

        ss.str("");
        ss.clear();

        // compute the statistics
        if (!stats_.empty()) {
            result += ", \"stats\": {";
            idx reps = 0;
            for (auto&& elem : get_stats())
                reps += elem.second;
            result += "\"reps\": " + std::to_string(reps) + ", ";
            result += "\"outcomes\": " + std::to_string(stats_.size()) + ", ";

            std::vector<idx> dits_dims(qc_->get_nc(), qc_->get_d());
            bool is_first = true;
            for (auto&& elem : get_stats()) {
                if (is_first)
                    is_first = false;
                else
                    ss << ", ";
                ss << "\""
                   << "[" << elem.first << "]"
                   << "\" : " << elem.second;
            }
            ss << '}';
            result += ss.str();
        }

        if (enclosed_in_curly_brackets)
            result += "}";

        return result;
    }

  private:
    std::ostream& display(std::ostream& os) const override {
        os << "measured/discarded (destructive): " << disp(get_measured(), ", ")
           << '\n';
        os << "non-measured/non-discarded: " << disp(get_non_measured(), ", ")
           << '\n';
        os << "last probs: " << disp(get_probs(), ", ") << '\n';
        os << "last dits: " << disp(get_dits(), ", ");

        // compute the statistics
        if (!stats_.empty()) {
            idx reps = 0;
            for (auto&& elem : get_stats())
                reps += elem.second;
            os << "\nstats:\n";
            os << '\t' << "reps: " << reps << '\n';
            os << '\t' << "outcomes: " << stats_.size() << '\n';
            std::vector<idx> dits_dims(qc_->get_nc(), qc_->get_d());
            bool is_first = true;
            for (auto&& elem : get_stats()) {
                if (is_first)
                    is_first = false;
                else
                    os << '\n';
                os << '\t' << "[" << elem.first << "]"
                   << ": " << elem.second;
            }
        }

        return os;
    }
}; /* class QEngine */

template <typename NoiseModel>
class QNoisyEngine : public QEngine {
    NoiseModel noise_;                            
    std::vector<std::vector<idx>> noise_results_; 
  public:
    explicit QNoisyEngine(const QCircuit& qc, const NoiseModel& noise)
        : QEngine{qc}, noise_{noise}, noise_results_(qc.get_step_count()) {
        // EXCEPTION CHECKS

        // check noise has the correct dimensionality
        if (qc.get_d() != noise.get_d())
            throw exception::DimsNotEqual("qpp::QNoisyEngine::QNoisyEngine()");
        // END EXCEPTION CHECKS
    }

    using QEngine::execute;

    QNoisyEngine& execute(const QCircuit::iterator::value_type& elem) override {
        // get the relative position of the target
        std::vector<idx> target_rel_pos = get_relative_pos_(get_non_measured());
        if (elem.type_ != QCircuit::StepType::MEASUREMENT) {
            // apply the noise
            for (auto&& i : target_rel_pos) {
                st_.psi_ = noise_(st_.psi_, i);
                // record the Kraus operator that occurred
                noise_results_[elem.ip_].emplace_back(noise_.get_last_idx());
            }
        }
        // execute the circuit step
        (void) QEngine::execute(elem);

        return *this;
    }

    // getters
    std::vector<std::vector<idx>> get_noise_results() const {
        return noise_results_;
    }
    // end getters
}; /* class QNoisyEngine */

} /* namespace qpp */

#endif /* CLASSES_CIRCUITS_ENGINES_H_ */