Program Listing for File Utilities.hpp

Return to documentation for file (include/jet/Utilities.hpp)

#pragma once

#include <algorithm>
#include <complex>
#include <iostream>
#include <numeric>
#include <string>
#include <vector>

#include "Abort.hpp"

namespace Jet {
namespace Utilities {

constexpr inline bool is_pow_2(size_t value)
{
    return static_cast<bool>(value && !(value & (value - 1)));
}

constexpr inline size_t fast_log2(size_t value)
{
    return static_cast<size_t>(std::numeric_limits<size_t>::digits -
                               __builtin_clzll((value)) - 1ULL);
}

template <class T1, class T2>
inline std::ostream &operator<<(std::ostream &os, const std::pair<T1, T2> &p)
{
    return os << '{' << p.first << ',' << p.second << '}';
}

template <class T>
inline std::ostream &operator<<(std::ostream &os, const std::vector<T> &v)
{
    os << '{';
    for (size_t i = 0; i < v.size(); i++) {
        if (i != 0) {
            os << "  ";
        }
        os << v[i];
    }
    os << '}';
    return os;
}

inline std::string GenerateStringIndex(size_t id)
{
    static const std::vector<std::string> alphabet = {
        "a", "b", "c", "d", "e", "f", "g", "h", "i", "j", "k", "l", "m",
        "n", "o", "p", "q", "r", "s", "t", "u", "v", "w", "x", "y", "z",
        "A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M",
        "N", "O", "P", "Q", "R", "S", "T", "U", "V", "W", "X", "Y", "Z"};
    const size_t div_id = id / alphabet.size();
    const std::string prefix = alphabet[id % alphabet.size()];
    const std::string suffix = (div_id == 0) ? "" : std::to_string(div_id - 1);
    return prefix + suffix;
}

template <class scalar_type_t>
inline size_t Order(const std::vector<std::complex<scalar_type_t>> &mat)
{
    // If sqrt() returns a value just under the true square root, increment n.
    size_t n = static_cast<size_t>(sqrt(mat.size()));
    n += n * n != mat.size();
    return n;
}

template <class scalar_type_t>
inline std::vector<std::complex<scalar_type_t>> Eye(size_t n)
{
    std::vector<std::complex<scalar_type_t>> eye(n * n, 0);
    for (size_t i = 0; i < n; i++) {
        eye[i * n + i] = 1;
    }
    return eye;
}

template <typename scalar_type_t>
inline std::vector<std::complex<scalar_type_t>>
MultiplySquareMatrices(const std::vector<std::complex<scalar_type_t>> &m1,
                       const std::vector<std::complex<scalar_type_t>> &m2,
                       size_t n)
{
    JET_ABORT_IF_NOT(m1.size() == n * n, "LHS matrix has the wrong order");
    JET_ABORT_IF_NOT(m2.size() == n * n, "RHS matrix has the wrong order");

    std::vector<std::complex<scalar_type_t>> product(n * n);
    for (size_t i = 0; i < n; i++) {
        for (size_t j = 0; j < n; j++) {
            for (size_t k = 0; k < n; k++) {
                product[i * n + j] += m1[i * n + k] * m2[k * n + j];
            }
        }
    }
    return product;
}

template <typename scalar_type_t>
inline std::vector<std::complex<scalar_type_t>>
Pow(const std::vector<std::complex<scalar_type_t>> &mat, size_t k)
{
    if (k == 1) {
        return mat;
    }

    const auto n = Order(mat);
    if (k == 0) {
        return Eye<scalar_type_t>(n);
    }

    std::vector<std::complex<scalar_type_t>> power = mat;
    for (size_t i = 2; i <= k; i++) {
        power = MultiplySquareMatrices(power, mat, n);
    }
    return power;
}

template <typename scalar_type_t>
inline std::vector<std::complex<scalar_type_t>>
operator+(const std::vector<std::complex<scalar_type_t>> &m1,
          const std::vector<std::complex<scalar_type_t>> &m2)
{
    JET_ABORT_IF_NOT(m1.size() == m2.size(), "Matrices have different sizes");

    std::vector<std::complex<scalar_type_t>> sum(m1.size());
    for (size_t i = 0; i < m1.size(); i++) {
        sum[i] = m1[i] + m2[i];
    }
    return sum;
}

template <typename scalar_type_t>
inline std::vector<std::complex<scalar_type_t>>
operator-(const std::vector<std::complex<scalar_type_t>> &m1,
          const std::vector<std::complex<scalar_type_t>> &m2)
{
    JET_ABORT_IF_NOT(m1.size() == m2.size(), "Matrices have different sizes");

    std::vector<std::complex<scalar_type_t>> diff(m1.size());
    for (size_t i = 0; i < m1.size(); i++) {
        diff[i] = m1[i] - m2[i];
    }
    return diff;
}

template <typename scalar_type_t>
inline std::vector<std::complex<scalar_type_t>>
operator*(const std::vector<std::complex<scalar_type_t>> &mat,
          std::complex<scalar_type_t> c)
{
    std::vector<std::complex<scalar_type_t>> product = mat;
    for (size_t i = 0; i < product.size(); i++) {
        product[i] *= c;
    }
    return product;
}

template <typename scalar_type_t>
inline std::vector<std::complex<scalar_type_t>>
DiagExp(const std::vector<std::complex<scalar_type_t>> &mat)
{
    const auto n = Order(mat);
    std::vector<std::complex<scalar_type_t>> diag(mat.size(), 0);
    for (size_t i = 0; i < n; i++) {
        diag[i * n + i] = std::exp(mat[i * n + i]);
    }
    return diag;
}

template <typename Tensor, typename T>
inline Tensor DiagMatrix(const std::vector<T> &vec)
{
    const size_t n = vec.size();
    Tensor tens({n, n});
    for (size_t i = 0; i < n; i++) {
        for (size_t j = 0; j < n; j++) {
            tens.SetValue({i, j}, (i == j) ? vec[i] : 0.0);
        }
    }
    return tens;
}

template <class T> inline bool InVector(const T &e, const std::vector<T> &v)
{
    return std::find(v.cbegin(), v.cend(), e) != v.cend();
}

template <typename T>
inline std::vector<T> VectorIntersection(const std::vector<T> &v1,
                                         const std::vector<T> &v2)
{
    std::vector<T> result;
    for (const auto &value : v1) {
        if (InVector(value, v2)) {
            result.emplace_back(value);
        }
    }
    return result;
}

template <typename T>
inline std::vector<T> VectorUnion(const std::vector<T> &v1,
                                  const std::vector<T> &v2)
{
    std::vector<T> result = v1;
    for (const auto &value : v2) {
        if (!InVector(value, v1)) {
            result.emplace_back(value);
        }
    }
    return result;
}

template <typename T>
inline std::vector<T> VectorSubtraction(const std::vector<T> &v1,
                                        const std::vector<T> &v2)
{
    std::vector<T> result;
    for (const auto &value : v1) {
        if (!InVector(value, v2)) {
            result.emplace_back(value);
        }
    }
    return result;
}

template <typename T>
inline std::vector<T> VectorDisjunctiveUnion(const std::vector<T> &v1,
                                             const std::vector<T> &v2)
{
    return VectorSubtraction(VectorUnion(v1, v2), VectorIntersection(v1, v2));
}

inline std::string JoinStringVector(const std::vector<std::string> &v)
{
    return std::accumulate(v.begin(), v.end(), std::string(""));
}

template <typename T>
inline std::vector<T> VectorConcatenation(const std::vector<T> &v1,
                                          const std::vector<T> &v2)
{
    std::vector<T> concat = v1;
    concat.insert(concat.end(), v2.begin(), v2.end());
    return concat;
}

inline size_t Factorial(size_t n)
{
    size_t prod = 1;
    for (size_t i = 2; i <= n; i++) {
        prod *= i;
    }
    return prod;
}

inline size_t ShapeToSize(const std::vector<size_t> &shape)
{
    size_t size = 1;
    for (const auto &dim : shape) {
        size *= dim;
    }
    return size;
}

inline std::vector<size_t> UnravelIndex(unsigned long long index,
                                        const std::vector<size_t> &shape)
{
    const size_t size = ShapeToSize(shape);
    JET_ABORT_IF(size <= index, "Linear index does not fit in the shape.");

    std::vector<size_t> multi_index(shape.size());
    for (int i = multi_index.size() - 1; i >= 0; i--) {
        multi_index[i] = index % shape[i];
        index /= shape[i];
    }
    return multi_index;
}

inline unsigned long long RavelIndex(const std::vector<size_t> &index,
                                     const std::vector<size_t> &shape)
{
    JET_ABORT_IF_NOT(index.size() == shape.size(),
                     "Number of index and shape dimensions must match.");

    size_t multiplier = 1;

    unsigned long long linear_index = 0;
    for (int i = index.size() - 1; i >= 0; i--) {
        JET_ABORT_IF(index[i] >= shape[i], "Index does not fit in the shape.");
        linear_index += index[i] * multiplier;
        multiplier *= shape[i];
    }
    return linear_index;
}

inline void SplitStringOnDelimiter(std::string &s, const std::string &delimiter,
                                   std::vector<std::string> &tokens)
{
    const size_t pos = s.find(delimiter);
    if (pos == std::string::npos) {
        tokens.emplace_back(s);
        return;
    }

    const auto token = s.substr(0, pos);
    tokens.emplace_back(token);

    s.erase(0, pos + delimiter.length());
    tokens.emplace_back(s);
}

inline std::vector<std::string>
SplitStringOnMultipleDelimiters(std::string s,
                                const std::vector<std::string> &delimiters)
{
    // Remove spaces.
    s.erase(std::remove_if(s.begin(), s.end(), isspace), s.end());

    std::vector<std::string> tokens = {s};
    for (std::size_t i = 0; i < delimiters.size(); i++) {
        tokens.pop_back();
        SplitStringOnDelimiter(s, delimiters[i], tokens);
    }

    // Remove empty tokens.
    const auto empty = [](const std::string &token) { return token.empty(); };
    tokens.erase(std::remove_if(tokens.begin(), tokens.end(), empty),
                 tokens.end());

    return tokens;
}

inline void SplitStringOnDelimiterRecursively(const std::string &s,
                                              const std::string &delimiter,
                                              std::vector<std::string> &tokens)
{
    const size_t pos = s.find(delimiter);
    if (pos == std::string::npos) {
        tokens.emplace_back(s);
    }
    else {
        tokens.emplace_back(s.begin(), s.begin() + pos);
        const auto remaining = s.substr(pos + delimiter.length());
        SplitStringOnDelimiterRecursively(remaining, delimiter, tokens);
    }
}

inline void ReplaceAllInString(std::string &s, const std::string &from,
                               const std::string &to)
{
    JET_ABORT_IF(from.empty(), "Cannot replace occurrences of an empty string");

    size_t pos = s.find(from, 0);
    while (pos != std::string::npos) {
        s.replace(pos, from.length(), to);
        // Skip over `to` in case part of it matches `from`.
        pos = s.find(from, pos + to.length());
    }
}

inline int GetTotalMemory()
{
    FILE *meminfo = fopen("/proc/meminfo", "r");
    if (meminfo == NULL) {
        return -1;
    }

    char line[256];
    while (fgets(line, sizeof(line), meminfo)) {
        int memTotal;
        if (sscanf(line, "MemTotal: %d kB", &memTotal) == 1) {
            fclose(meminfo);
            return memTotal;
        }
    }

    // Getting here means we were not able to find what we were looking for
    fclose(meminfo);
    return -1;
}

inline int GetAvailableMemory()
{
    /* Same function as above but it parses the meminfo file
       in order to obtain the current amount of physical memory available
    */
    FILE *meminfo = fopen("/proc/meminfo", "r");
    if (meminfo == NULL) {
        return -1;
    }

    char line[256];
    while (fgets(line, sizeof(line), meminfo)) {
        int memAvail;
        if (sscanf(line, "MemAvailable: %d kB", &memAvail) == 1) {
            fclose(meminfo);
            return memAvail;
        }
    }

    fclose(meminfo);
    return -1;
}

template <typename T> void FastCopy(const std::vector<T> &a, std::vector<T> &b)
{
#ifdef _OPENMP
    size_t max_right_dim = 1024;
    size_t size = a.size();
    if (b.size() != size)
        b.resize(size);
#pragma omp parallel for schedule(static, max_right_dim)
    for (std::size_t p = 0; p < size; ++p) {
        b[p] = a[p];
    }
#else
    b = a;
#endif
}

template <typename T>
bool VectorInVector(const std::vector<T> &v, const std::vector<T> &w)
{
    for (std::size_t i = 0; i < v.size(); ++i) {
        if (!InVector(v[i], w))
            return false;
    }
    return true;
}

}; // namespace Utilities
}; // namespace Jet