utilities.hpp

Go to the documentation of this file.

// Copyright 2015-2018 The ALMA Project Developers
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//   http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
// implied. See the License for the specific language governing
// permissions and limitations under the License.

#pragma once


#include <tuple>
#include <vector>
#include <string>
#include <memory>
#include <limits>
#include <algorithm>
#include <iostream>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include <boost/math/common_factor.hpp>
#include <boost/algorithm/string.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/tokenizer.hpp>
#include <boost/mpi.hpp>
#include <constants.hpp>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <unsupported/Eigen/Splines>

namespace alma {
template <typename T, typename... Args>
std::unique_ptr<T> make_unique(Args&&... args) {
    return std::unique_ptr<T>(new T(std::forward<Args>(args)...));
}


template <typename T> class Container_comparator {
public:
    inline bool operator()(const T& op1, const T& op2) const {
        return std::lexicographical_compare(
            std::begin(op1), std::end(op1), std::begin(op2), std::end(op2));
    }
};

template <typename T>
std::tuple<std::vector<typename T::key_type>,
           std::vector<typename T::mapped_type>>
split_keys_and_values(const T& input) {
    auto size = input.size();

    std::vector<typename T::key_type> keys;
    keys.reserve(size);
    std::vector<typename T::mapped_type> values;
    values.reserve(size);

    for (auto& i : input) {
        keys.emplace_back(i.first);
        values.emplace_back(i.second);
    }
    return std::make_tuple(keys, values);
}


inline void flush_istream(std::istream& f) {
    if (f.peek() == '\n') {
        f.ignore(1, '\n');
    }
}


template <typename T>
inline std::vector<T> tokenize_homogeneous_line(std::string& line,
                                                const std::string sep = " \t") {
    std::vector<T> nruter;
    boost::char_separator<char> separator(sep.c_str());
    boost::tokenizer<boost::char_separator<char>> tokenizer(line, separator);

    for (auto i : tokenizer)
        nruter.emplace_back(boost::lexical_cast<T>(i));
    return nruter;
}


inline bool starts_with_character(const std::string& line,
                                  const std::string& chars) {
    if (line.length() == 0)
        return false;

    for (auto c : chars) {
        if (line[0] == c)
            return true;
    }
    return false;
}


inline bool almost_equal(const double a,
                         const double b,
                         const double abstol = 1e-8,
                         const double reltol = 1e-5) {
    // Method and default values inspired by
    // numpy.allclose().
    const auto tol = std::fabs(abstol) + std::fabs(reltol * b);

    return std::fabs(a - b) < tol;
}


template <typename T> class Min_keeper {
public:
    Min_keeper(const double _abstol = 1e-8, const double _reltol = 1e-5)
        : abstol(_abstol), reltol(_reltol) {
    }


    double get_value() const {
        return this->minimum;
    }


    std::vector<T> get_vector() const {
        return this->objects;
    }


    void update(const T& obj, double value) {
        if (almost_equal(value, this->minimum, this->abstol, this->reltol)) {
            this->objects.emplace_back(obj);
        }
        else if (value < this->minimum) {
            this->minimum = value;
            this->objects.clear();
            this->objects.emplace_back(obj);
        }
    }


private:
    double minimum = std::numeric_limits<double>::infinity();
    std::vector<T> objects;
    const double abstol;
    const double reltol;
};

inline double ssqrt(double x) {
    return std::copysign(std::sqrt(std::fabs(x)), x);
}


template <typename T> inline int signum(const T& arg) {
    T zero(0);

    return (zero < arg) - (arg < zero);
}


inline std::string get_timestamp() {
    // A much more elegant solution using put_time() exists.
    // We chose the traditional implementation because put time
    // is not implemented in GCC < 5.
    auto t = std::time(nullptr);
    auto localtime = std::localtime(&t);
    std::size_t size = 0;

    std::vector<char> buffer;

    do {
        size += 64;
        buffer.reserve(size);
        // 64 bytesshould be enough for everybody, but we
        // check the return code of std::strftime anyway to see
        // if we need to allocate more space.
    } while (std::strftime(
                 buffer.data(), size, "%Y-%m-%dT%H:%M:%S", localtime) == 0u);
    return std::string(buffer.data());
}


inline std::array<std::size_t, 2> my_jobs(std::size_t njobs,
                                          std::size_t nprocs,
                                          std::size_t my_id) {
    // If nprocs divides njobs, each process gets exactly its
    // fair share of jobs. If there is a remainder R , an extra
    // job is assigned to the first R processes. If njobs > nprocs,
    // some processes do not get any job.
    auto quot = njobs / nprocs;
    auto rem = njobs % nprocs;

    std::array<std::size_t, 2> nruter;
    nruter[0] = quot * my_id + std::min(rem, my_id);
    nruter[1] = quot * (my_id + 1) + std::min(rem, my_id + 1);
    return nruter;
}


template <typename T, typename U> inline U python_mod(const T& n, const U& d) {
    if (d < 0)
        return python_mod(-n, -d);
    else {
        auto nruter = n % d;

        if (nruter < 0)
            nruter += d;
        return nruter;
    }
}


inline double bose_einstein(double omega, double T) {
    constexpr double prefactor = 1e12 * constants::hbar / constants::kB;

    return 1. / (std::exp(prefactor * omega / T) - 1.);
}


inline double bose_einstein_kernel(double omega, double T) {
    constexpr double prefactor = 0.5 * 1e12 * constants::hbar / constants::kB;
    auto x = prefactor * omega / T;

    if (x < 1e-320) {
        return 1.;
    }
    else if (x > 300.) {
        return 0.;
    }
    else {
        auto s = x / std::sinh(x);
        return s * s;
    }
}


inline Eigen::VectorXi reduce_integers(
    const Eigen::Ref<const Eigen::VectorXi>& vector) {
    std::size_t d = vector.size();
    Eigen::VectorXi::Scalar gcd;

    if (d == 0)
        gcd = 1;
    else {
        gcd = std::abs(vector(0));

        for (std::size_t i = 1; i < d; ++i)
            gcd = boost::math::gcd(gcd, vector(i));
    }
    return Eigen::VectorXi(vector / gcd);
}


inline std::string engineer_format(double x, bool insert_space = false) {
    std::vector<char> prefix(
        {'f', 'p', 'n', 'u', 'm', ' ', 'k', 'M', 'G', 'T'});
    std::vector<double> scaling(
        {1e-15, 1e-12, 1e-9, 1e-6, 1e-3, 1.0, 1e3, 1e6, 1e9, 1e12, 1e15});

    int idx =
        static_cast<int>(std::floor((std::log10(x) + 15.0) / (3.0 - 1e-12)));

    if (idx < 0) {
        idx = 0;
    }

    if (idx > 9) {
        idx = 9;
    }
    std::stringstream result_builder;
    result_builder << x / scaling.at(idx);

    if (idx != 5) { // avoid introducing blank space
        if (insert_space) {
            result_builder << " ";
        }
        result_builder << prefix.at(idx);
    }
    return result_builder.str();
}


inline Eigen::VectorXd logSpace(double min, double max, int N) {
    Eigen::VectorXd logresult(N);

    logresult.setLinSpaced(N, std::log10(min), std::log10(max));
    return Eigen::pow(10.0, logresult.array());
}


inline Eigen::VectorXd splineInterpolation(
    const Eigen::Ref<const Eigen::VectorXd>& X,
    const Eigen::Ref<const Eigen::VectorXd>& Y,
    const Eigen::Ref<const Eigen::VectorXd>& Xtarget) {
    Eigen::VectorXd result(Xtarget.size());

    // map X values to 1D chord lengths ranging from 0 to 1
    Eigen::VectorXd X_norm =
        (X.array() - X.minCoeff()) / (X.maxCoeff() - X.minCoeff());

    // construct cubic spline over the base points
    typedef Eigen::Spline<double, 1> Spline1D;
    Spline1D S(Eigen::SplineFitting<Spline1D>::Interpolate(
        Y.transpose(), 3, X_norm.transpose()));

    // convert Xtarget values to their equivalent chord length
    Eigen::VectorXd Xtarget_norm =
        (Xtarget.array() - X.minCoeff()) / (X.maxCoeff() - X.minCoeff());

    // obtain spline value at each of these chord lengths

    for (int nx = 0; nx < Xtarget.size(); nx++) {
        result(nx) = S(Xtarget_norm(nx))(0);
    }

    return result;
}


inline Eigen::VectorXd linearInterpolation(
    const Eigen::Ref<const Eigen::VectorXd>& X,
    const Eigen::Ref<const Eigen::VectorXd>& Y,
    const Eigen::Ref<const Eigen::VectorXd>& Xtarget) {
    Eigen::VectorXd result(Xtarget.size());
    int N = X.size();

    // perform linear interpolation for each of the specified targets

    for (int nx = 0; nx < Xtarget.size(); nx++) {
        if (Xtarget(nx) < X(0)) {
            std::cout << "alma::linearInterpolation > " << std::endl;
            std::cout << "WARNING: encountered target element out of bounds."
                      << std::endl;
            std::cout << "Interpolated result will be clipped." << std::endl;

            result(nx) = X(0);
        }

        else if (Xtarget(nx) > X(N - 1)) {
            std::cout << "alma::linearInterpolation > " << std::endl;
            std::cout << "WARNING: encountered target element out of bounds."
                      << std::endl;
            std::cout << "Interpolated result will be clipped." << std::endl;

            result(nx) = X(N - 1);
        }

        else {
            // use bisection algorithm to determine base array bin the target
            // sits
            // in.

            int leftindex = 0;
            int rightindex = N - 1;

            while (rightindex - leftindex > 1) {
                int midindex = (leftindex + rightindex) / 2;

                if (X(midindex) > Xtarget(nx)) {
                    rightindex = midindex;
                }
                else {
                    leftindex = midindex;
                }
            }

            double yleft = Y(leftindex);
            double yright = Y(rightindex);
            double xleft = X(leftindex);
            double xright = X(rightindex);

            result(nx) = yleft + (Xtarget(nx) - xleft) * (yright - yleft) /
                                     (xright - xleft);
        }
    }

    return result;
}
} // namespace alma

namespace boost {
namespace serialization {
template <class Archive, typename S>
void save(Archive& ar, const Eigen::Triplet<S>& t, const unsigned int version) {
    ar& t.row();
    ar& t.col();
    ar& t.value();
}


template <class Archive, typename S>
void load(Archive& ar, Eigen::Triplet<S>& t, const unsigned int version) {
    int row;
    ar& row;
    int col;
    ar& col;
    S value;
    ar& value;

    t = Eigen::Triplet<S>(row, col, value);
}


template <class Archive, typename S>
void serialize(Archive& ar, Eigen::Triplet<S>& t, const unsigned int version) {
    split_free(ar, t, version);
}
} // namespace serialization
} // namespace boost