entropies.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 ENTROPY_H_
#define ENTROPY_H_

namespace qpp {
template <typename Derived>
double entropy(const Eigen::MatrixBase<Derived>& A) {
    const dyn_mat<typename Derived::Scalar>& rA = A.derived();

    // EXCEPTION CHECKS

    // check zero-size
    if (!internal::check_nonzero_size(rA))
        throw exception::ZeroSize("qpp::entropy()");

    // check square matrix
    if (!internal::check_square_mat(rA))
        throw exception::MatrixNotSquare("qpp::entropy()");
    // END EXCEPTION CHECKS

    dmat ev = svals(rA); // get the singular values
    double result = 0;
    for (idx i = 0; i < static_cast<idx>(ev.rows()); ++i)
        if (ev(i) != 0) // not identically zero
            result -= ev(i) * std::log2(ev(i));

    return result;
}

inline double entropy(const std::vector<double>& prob) {
    // EXCEPTION CHECKS

    // check zero-size
    if (!internal::check_nonzero_size(prob))
        throw exception::ZeroSize("qpp::entropy()");
    // END EXCEPTION CHECKS

    double result = 0;
    for (idx i = 0; i < prob.size(); ++i)
        if (std::abs(prob[i]) != 0) // not identically zero
            result -= std::abs(prob[i]) * std::log2(std::abs(prob[i]));

    return result;
}

template <typename Derived>
double renyi(const Eigen::MatrixBase<Derived>& A, double alpha) {
    const dyn_mat<typename Derived::Scalar>& rA = A.derived();

    // EXCEPTION CHECKS

    // check zero-size
    if (!internal::check_nonzero_size(rA))
        throw exception::ZeroSize("qpp::renyi()");

    // check square matrix
    if (!internal::check_square_mat(rA))
        throw exception::MatrixNotSquare("qpp::renyi()");

    if (alpha < 0)
        throw exception::OutOfRange("qpp::renyi()");
    // END EXCEPTION CHECKS

    if (alpha == 0) // H max
        return std::log2(rA.rows());

    if (alpha == 1) // Shannon/von-Neumann
        return entropy(rA);

    if (alpha == infty) // H min
        return -std::log2(svals(rA)[0]);

    dmat sv = svals(rA); // get the singular values
    double result = 0;
    for (idx i = 0; i < static_cast<idx>(sv.rows()); ++i)
        if (sv(i) != 0) // not identically zero
            result += std::pow(sv(i), alpha);

    return std::log2(result) / (1 - alpha);
}

inline double renyi(const std::vector<double>& prob, double alpha) {
    // EXCEPTION CHECKS

    // check zero-size
    if (!internal::check_nonzero_size(prob))
        throw exception::ZeroSize("qpp::renyi()");

    if (alpha < 0)
        throw exception::OutOfRange("qpp::renyi()");

    if (alpha == 0) // H max
        return std::log2(prob.size());

    if (alpha == 1) // Shannon/von-Neumann
        return entropy(prob);
    // END EXCEPTION CHECKS

    if (alpha == infty) // H min
    {
        double max = 0;
        for (idx i = 0; i < prob.size(); ++i)
            if (std::abs(prob[i]) > max)
                max = std::abs(prob[i]);

        return -std::log2(max);
    }

    double result = 0;
    for (idx i = 0; i < prob.size(); ++i)
        if (std::abs(prob[i]) != 0) // not identically zero
            result += std::pow(std::abs(prob[i]), alpha);

    return std::log2(result) / (1 - alpha);
}

template <typename Derived>
double tsallis(const Eigen::MatrixBase<Derived>& A, double q) {
    const dyn_mat<typename Derived::Scalar>& rA = A.derived();

    // EXCEPTION CHECKS

    // check zero-size
    if (!internal::check_nonzero_size(rA))
        throw exception::ZeroSize("qpp::tsallis()");

    // check square matrix
    if (!internal::check_square_mat(rA))
        throw exception::MatrixNotSquare("qpp::tsallis()");

    if (q < 0)
        throw exception::OutOfRange("qpp::tsallis()");
    // END EXCEPTION CHECKS

    if (q == 1) // Shannon/von-Neumann with base e logarithm
        return entropy(rA) * std::log(2);

    dmat ev = svals(rA); // get the singular values
    double result = 0;
    for (idx i = 0; i < static_cast<idx>(ev.rows()); ++i)
        if (ev(i) != 0) // not identically zero
            result += std::pow(ev(i), q);

    return (result - 1) / (1 - q);
}

inline double tsallis(const std::vector<double>& prob, double q) {
    // EXCEPTION CHECKS

    // check zero-size
    if (!internal::check_nonzero_size(prob))
        throw exception::ZeroSize("qpp::tsallis()");

    if (q < 0)
        throw exception::OutOfRange("qpp::tsallis()");
    // END EXCEPTION CHECKS

    if (q == 1) // Shannon/von-Neumann with base e logarithm
        return entropy(prob) * std::log(2);

    double result = 0;
    for (idx i = 0; i < prob.size(); ++i)
        if (std::abs(prob[i]) != 0) // not identically zero
            result += std::pow(std::abs(prob[i]), q);

    return (result - 1) / (1 - q);
}

template <typename Derived>
double qmutualinfo(const Eigen::MatrixBase<Derived>& A,
                   const std::vector<idx>& subsysA,
                   const std::vector<idx>& subsysB,
                   const std::vector<idx>& dims) {
    const dyn_mat<typename Derived::Scalar>& rA = A.derived();

    // EXCEPTION CHECKS

    // check zero-size
    if (!internal::check_nonzero_size(rA))
        throw exception::ZeroSize("qpp::qmutualinfo()");

    // check that dims is a valid dimension vector
    if (!internal::check_dims(dims))
        throw exception::DimsInvalid("qpp::qmutualinfo()");

    // check square matrix
    if (!internal::check_square_mat(rA))
        throw exception::MatrixNotSquare("qpp::qmutualinfo()");

    // check that dims match the dimension of A
    if (!internal::check_dims_match_mat(dims, rA))
        throw exception::DimsMismatchMatrix("qpp::qmutualinfo()");

    // check that subsys are valid
    if (!internal::check_subsys_match_dims(subsysA, dims) ||
        !internal::check_subsys_match_dims(subsysB, dims))
        throw exception::SubsysMismatchDims("qpp::qmutualinfo()");
    // END EXCEPTION CHECKS

    // The full system indexes {0,1,...,n-1}
    std::vector<idx> full_system(dims.size());
    std::iota(std::begin(full_system), std::end(full_system), 0);

    // Sorted input subsystems
    std::vector<idx> subsysAsorted{subsysA};
    std::vector<idx> subsysBsorted{subsysB};

    // sort the input subsystems (as needed by std::set_union)
    std::sort(std::begin(subsysAsorted), std::end(subsysAsorted));
    std::sort(std::begin(subsysBsorted), std::end(subsysBsorted));

    // construct the union of A and B
    std::vector<idx> subsysAB;
    std::set_union(std::begin(subsysAsorted), std::end(subsysAsorted),
                   std::begin(subsysBsorted), std::end(subsysBsorted),
                   std::back_inserter(subsysAB));

    // construct the complements
    std::vector<idx> subsysA_bar = complement(subsysA, dims.size());
    std::vector<idx> subsysB_bar = complement(subsysB, dims.size());
    std::vector<idx> subsysAB_bar = complement(subsysAB, dims.size());

    cmat rhoA = ptrace(rA, subsysA_bar, dims);
    cmat rhoB = ptrace(rA, subsysB_bar, dims);
    cmat rhoAB = ptrace(rA, subsysAB_bar, dims);

    return entropy(rhoA) + entropy(rhoB) - entropy(rhoAB);
}

template <typename Derived>
double qmutualinfo(const Eigen::MatrixBase<Derived>& A,
                   const std::vector<idx>& subsysA,
                   const std::vector<idx>& subsysB, idx d = 2) {
    const dyn_mat<typename Derived::Scalar>& rA = A.derived();

    // EXCEPTION CHECKS

    // check zero-size
    if (!internal::check_nonzero_size(rA))
        throw exception::ZeroSize("qpp::qmutualinfo()");

    // check valid dims
    if (d < 2)
        throw exception::DimsInvalid("qpp::qmutualinfo()");
    // END EXCEPTION CHECKS

    idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
    std::vector<idx> dims(n, d); // local dimensions vector

    return qmutualinfo(rA, subsysA, subsysB, dims);
}

} /* namespace qpp */

#endif /* ENTROPY_H_ */