utils.hpp

Go to the documentation of this file.

#ifndef UTILS_HPP
#define UTILS_HPP 1

#include <vector>
#include <algorithm>
#include "universal_error.hpp"
#include <cassert>
#include "boost/container/flat_map.hpp"
#include <iostream>

using std::vector;

bool is_nan(double x);

vector<double> linspace(double xl, double xh, int n);

vector<double> arange(double x_min, double x_max, double dx);

template <class T> vector<T> operator+
(vector<T> const& v1, vector<T> const& v2)
{
        assert(v1.size() == v2.size() && "Vectors must have the same length");
        vector<T> res(v1.size());
        for (size_t i = 0; i < v1.size(); ++i)
                res[i] = v1[i] + v2[i];
        return res;
}

template<class T> vector<T> operator+
(const vector<T>& v, const T& t)
{
        vector<T> res(v.size());
        for (size_t i = 0, endp = v.size(); i < endp; ++i)
                res[i] = v[i] + t;
        return res;
}

template<class T> vector<T> operator+
(const T& t, const vector<T>& v)
{
        return v + t;
}

template<class T> vector<T>& operator+=
(vector<T>& v, const T& t)
{
        for (size_t i = 0, endp = v.size(); i < endp; ++i)
                v[i] += t;
        return v;
}

template <class T> vector<T> operator-
(vector<T> const& v1, vector<T> const& v2)
{
        if (v1.size() != v2.size())
                throw UniversalError("Vectors must have the same length");

        vector<T> res(v1.size());
        for (size_t i = 0; i < v1.size(); ++i)
                res[i] = v1[i] - v2[i];
        return res;
}

template <class T> vector<T> operator*
(double d, vector<T> const& v)
{
        vector<T> res(v.size());
        for (size_t i = 0; i < v.size(); ++i)
                res[i] = d*v[i];
        return res;
}

template <typename T>
T BiLinearInterpolation(const vector<T>& x, const vector<T>& y, const std::vector<std::vector<T> >& z, T xi, T yi);

template <typename T>
T BiLinearInterpolation(const vector<T>& x, const vector<T>& y, const std::vector<std::vector<T> >& z, T xi, T yi)
{
        typename vector<T>::const_iterator itx = upper_bound(x.begin(), x.end(), xi);
        if (itx == x.end())
        {
                if (xi > x.back())
                {
                        std::cout << "X too large in BiLinearInterpolation, x_i " << xi << " max X " << x.back() << std::endl;
                        throw;
                }
                else
                        --itx;
        }
        if (itx == x.begin())
        {
                if (*itx < x.at(0))
                {
                        std::cout << "X too small in BiLinearInterpolation, x_i " << xi << " min X " << x.at(0) << std::endl;
                        throw;
                }
        }
        typename vector<T>::const_iterator ity = upper_bound(y.begin(), y.end(), yi);
        if (ity == y.end())
        {
                if (yi > y.back())
                {
                        std::cout << "Y too large in LinearInterpolation, y_i " << yi << " max Y " << y.back() << std::endl;
                        throw;
                }
                else
                        --ity;
        }
        if (ity == y.begin())
        {
                if (*ity < y.at(0))
                {
                        std::cout << "Y too small in LinearInterpolation, y_i " << yi << " min Y " << y.at(0) << std::endl;
                        throw;
                }
        }

        T delta = 1.0 / ((*itx - *(itx - 1)) * (*ity - *(ity - 1)));
        size_t idx = static_cast<std::size_t>(itx - x.begin());
        size_t idy = static_cast<std::size_t>(ity - y.begin());
        return (z[idx - 1][idy - 1] * (*itx - xi) * (*ity - yi) + z[idx][idy - 1] * (xi - *(itx - 1)) * (*ity - yi)
                + z[idx - 1][idy] * (*itx - xi) * (yi - *(ity - 1)) + z[idx][idy] * (xi - *(itx - 1))
                * (yi - *(ity - 1))) * delta;
}

template <typename T>
T LinearInterpolation(const vector<T> &x, const vector<T> &y, T xi);

template <typename T>
T LinearInterpolation(const vector<T> &x, const vector<T> &y, T xi)
{
        typename vector<T>::const_iterator it = upper_bound(x.begin(), x.end(), xi);
        if (it == x.end())
        {
                std::cout << "X too large in LinearInterpolation, x_i " << xi << " max X " << x.back() << std::endl;
                throw;
        }
        if (it == x.begin())
        {
                if (*it < x.at(0))
                {
                        std::cout << "X too small in LinearInterpolation, x_i " << xi << " min X " << x.at(0) << std::endl;
                        throw;
                }
        }
        if (*it == xi)
                return y[static_cast<std::size_t>(it - x.begin())];

        return y[static_cast<std::size_t>(it - x.begin())] + (xi - *it)*        (y[static_cast<std::size_t>(it - 1 - x.begin())] 
                - y[static_cast<std::size_t>(it - x.begin())]) / (*(it - 1) - *it);
}

template <typename T>
T LinearInterpolation(typename vector<T>::const_iterator itx_begin, typename vector<T>::const_iterator itx_end,
        typename vector<T>::const_iterator ity_begin, T xi);

template <typename T>
T LinearInterpolation(typename vector<T>::const_iterator itx_begin, typename vector<T>::const_iterator itx_end,
        typename vector<T>::const_iterator ity_begin, T xi)
{
        typename vector<T>::const_iterator it = upper_bound(itx_begin, itx_end, xi);
        if (it == itx_end)
        {
                std::cout << "X too large in LinearInterpolation, x_i " << xi << " max X " << *itx_end << std::endl;
                throw;
        }
        if (it == itx_begin)
        {
                if (*it < *itx_begin)
                {
                        std::cout << "X too small in LinearInterpolation, x_i " << xi << " min X " << *itx_begin << std::endl;
                        throw;
                }
        }
        const size_t toadd = static_cast<std::size_t>(it - itx_begin);
        if (*it == xi)
                return *(ity_begin + toadd);

        return *(ity_begin + toadd) + (xi - *it) * (*(ity_begin + toadd + 1) - *(ity_begin + toadd)) / (*(it - 1) - *it);
}

double min(vector<double> const& v);

double max(vector<double> const& v);

template <class T> void RemoveVector
(vector<T> &v, vector<int> &indeces)
{
        if (indeces.empty())
                return;
        sort(indeces.begin(), indeces.end());
        vector<T> result;
        result.reserve(v.size() - indeces.size());
        int counter = 0;
        for (size_t i = 0; i < static_cast<size_t>(indeces.back()); ++i)
        {
                if (size_t(indeces[size_t(counter)]) == i)
                        ++counter;
                else
                        result.push_back(v[i]);
        }
        for (size_t i = static_cast<size_t>(indeces.back()) + 1; i < v.size(); ++i)
                result.push_back(v[i]);
        v = result;
}

template <class T> void RemoveVector
(vector<T> &v, vector<size_t> &indeces)
{
        if (indeces.empty())
                return;
        sort(indeces.begin(), indeces.end());
        vector<T> result;
        result.reserve(v.size() - indeces.size());
        size_t counter = 0;
        for (size_t i = 0; i < indeces.back(); ++i)
        {
                if (indeces[counter] == i)
                        ++counter;
                else
                        result.push_back(v[i]);
        }
        for (size_t i = indeces.back() + 1; i < v.size(); ++i)
                result.push_back(v[i]);
        v = result;
}
template <class T> vector<T> VectorValues(vector<T> const&v, vector<int> const &index)
{
        if (index.empty() || v.empty())
                return vector<T>();

        vector<T> result(index.size());
        for (size_t i = 0; i < index.size(); ++i)
                result.at(i) = v.at(static_cast<size_t>(index.at(i)));
        return result;
}

template <class T> vector<T> VectorValues
(vector<T> const&v,
        vector<size_t> const &index)
{
        if (index.empty() || v.empty())
                return vector<T>();

        vector<T> result(index.size());
        for (size_t i = 0; i < index.size(); ++i)
                result.at(i) = v.at(index.at(i));
        return result;
}

template <class T> T VectorSum(vector<T> const&v)
{
        if (v.empty())
                return 0;
        T result = v[0];
        for (size_t i = 1; i < v.size(); ++i)
                result += v[i];
        return result;
}

template <class T> vector<T> unique(vector<T> const& v)
{
        size_t n = v.size();
        vector<T> res;
        res.reserve(n);
        if (n == 0)
                return res;
        res.push_back(v[0]);
        for (typename vector<T>::const_iterator it = v.begin() + 1; it != v.end(); ++it)
                if (*it == *(it - 1))
                        continue;
                else
                        res.push_back(*it);
        return res;
}

template <class T> vector<int> unique_index(vector<T> const& v)
{
        if (v.empty())
                return vector<int>();

        vector<int> res;
        res.push_back(0);
        for (size_t i = 1; i < v.size(); ++i)
                if (v[i] != v[i - 1])
                        res.push_back(static_cast<int>(i));
        return res;
}

template <class T> vector<T> RemoveList(vector<T> const&v, vector<T> const&list)
{
        vector<T> res;
        for (size_t i = 0; i < v.size(); ++i)
                if (!binary_search(list.begin(), list.end(), v[i]))
                        res.push_back(v[i]);
        return res;
}

template <class T> void RemoveVal(vector<T> &vec, T val)
{
        for (size_t i = 0; i < vec.size(); ++i)
        {
                if (vec[i] == val)
                {
                        vec.erase(vec.begin() + static_cast<long>(i));
                        return;
                }
        }
}

template <class T> bool InVector(vector<T> const&vec, T val)
{
        int n = vec.size();
        for (int i = 0; i < n; ++i)
                if (vec[i] == val)
                        return true;
        return false;
}

template <class T> int IndexInVector(vector<T> const&vec, T val)
{
        int n = vec.size();
        for (int i = 0; i < n; ++i)
                if (vec[i] == val)
                        return i;
        throw UniversalError("Value not found in vector");
}
template <class T> void ReArrangeVector(vector<T> &v, vector<int> const& indeces)
{
        const vector<T> temp = v;
        for (size_t i = 0; i < v.size(); ++i)
                v[i] = temp[static_cast<size_t>(indeces[i])];
}
namespace
{
        template<class T> struct index_cmp
        {
                explicit index_cmp(const T _arr) : arr(_arr) {}
                bool operator()(const size_t a, const size_t b) const
                {
                        return arr[a] < arr[b];
                }
        private:
                const T arr;
        };
}
template<class T> void sort_index(const vector<T> & arr, vector<int>& res)
{
        res.resize(arr.size());
        for (size_t i = 0; i < res.size(); ++i)
                res[i] = static_cast<int>(i);
        sort(res.begin(), res.end(), index_cmp<vector<T> >(arr));
}

template<class T> void sort_index(const vector<T> & arr, vector<size_t>& res)
{
        res.resize(arr.size());
        for (size_t i = 0; i < res.size(); ++i)
                res[i] = i;
        sort(res.begin(), res.end(), index_cmp<vector<T> >(arr));
}


template<class T> vector<size_t> sort_index(const vector<T> & arr)
{
        vector<size_t> res(arr.size());
        for (size_t i = 0; i < res.size(); ++i)
                res[i] = i;
        sort(res.begin(), res.end(), index_cmp<vector<T> >(arr));
        return res;
}

namespace
{
        template <class RAIter, class Compare>
        class PairComp
        {
        public:
                Compare comp;
                explicit PairComp(Compare comp_) : comp(comp_) {}
                bool operator() (const std::pair<size_t, RAIter>& a,
                        const std::pair<size_t, RAIter>& b) const {
                        return comp(*a.second, *b.second);
                }
        };
}

template <class RAIter, class Compare>
void sort_index(RAIter iterBegin, RAIter iterEnd, Compare comp,
        std::vector<size_t>& indexes)
{

        std::vector< std::pair<size_t, RAIter> > pv;
        pv.reserve(iterEnd - iterBegin);

        RAIter iter;
        size_t k;
        for (iter = iterBegin, k = 0; iter != iterEnd; iter++, k++) {
                pv.push_back(std::pair<int, RAIter>(k, iter));
        }
        PairComp<RAIter, Compare> compy(comp);
        std::sort(pv.begin(), pv.end(), compy);

        indexes.resize(pv.size());
        for (size_t i = 0; i < pv.size(); ++i)
                indexes[i] = pv[i].first;
}

template<class T> vector<T> join(vector<T> const& v1,
        vector<T> const& v2)
{
        if (v1.empty() && v2.empty())
                return vector<T>();
        vector<T> res(v1.size() + v2.size());
        if (!v1.empty())
                copy(v1.begin(), v1.end(), res.begin());
        if (!v2.empty())
                copy(v2.begin(), v2.end(), res.begin() + static_cast<int>(v1.size()));
        return res;
}

template<class T> class BinaryOperation
{
public:

        virtual T operator()(T const& t1, T const& t2) const = 0;

        virtual ~BinaryOperation(void) {}
};

template<class T> vector<T> binary_unite(vector<T> const& v1,
        vector<T> const& v2,
        BinaryOperation<T> const& bin_op)
{
        assert(v1.size() == v2.size());

        vector<T> res(v1.size());
        for (size_t i = 0; i < v1.size(); ++i)
                res[i] = bin_op(v1[i], v2[i]);
        return res;
}

template<class T> class UnaryOperation
{
public:

        virtual T operator()(T const& t) const = 0;

        virtual ~UnaryOperation(void) {}
};

template<class T> vector<T> apply_to_each_term(vector<T> const& v,
        UnaryOperation<T> const& un_op)
{
        vector<T> res(v.size());
        for (size_t i = 0; i < v.size(); ++i)
                res[i] = un_op(v[i]);
        return res;
}

template<class T> T pair_member(const std::pair<T, T>& p, int index)
{
        assert((0 == index) || (1 == index));
        return 0 == index ? p.first : p.second;
}

template<class T> void set_pair_member(std::pair<T, T>& p, int index, const T& val)
{
        assert((0 == index) || (1 == index));
        if (0 == index)
                p.first = val;
        else
                p.second = val;
}

template<class T, class S> vector<T> list_static_cast(const vector<S>& source)
{
        vector<T> res(source.size());
        for (size_t i = 0; i < res.size(); ++i)
                res[i] = static_cast<T>(source[i]);
        return res;
}

template<class T> void insert_all_to_back(vector<T>& subject, const vector<T>& addendum)
{
        if (!addendum.empty())
                subject.insert(subject.end(), addendum.begin(), addendum.end());
}
template<class Iter, class T>
Iter binary_find(Iter begin, Iter end, T val)
{
        // Finds the lower bound in at most log(last - first) + 1 comparisons
        Iter i = std::lower_bound(begin, end, val);

        if (i != end && !(val < *i))
                return i; // found
        else
                return end; // not found
}

template<class S, class T> typename vector<T>::const_reference safe_retrieve
(vector<T> const& data, vector<S> const& keys, const S& key)
{
        assert(data.size() == keys.size());
        typename vector<S>::const_iterator it = binary_find(keys.begin(), keys.end(), key);
        assert(it != keys.end());
        size_t index = static_cast<size_t>(it - keys.begin());
        return data.at(index);
}

template<class S, class T> typename vector<T>::reference safe_retrieve
(vector<T> &data, vector<S> const& keys, const S& key)
{
        assert(data.size() == keys.size());
        typename vector<S>::const_iterator it = binary_find(keys.begin(), keys.end(), key);
        assert(it != keys.end());
        return data[static_cast<size_t>(it - keys.begin())];
}

template<class S, class T> const T& safe_retrieve
(const boost::container::flat_map<S, T>& m,
        const S& s)
{
        typename boost::container::flat_map<S, T>::const_iterator it =
                m.find(s);
        assert(it != m.end());
        return it->second;
}
template<class T> vector<vector<T> > CombineVectors(vector<vector<vector<T> > > const& data)
{
        size_t counter = 0;
        size_t out_size = data.size();
        for (size_t i = 0; i < out_size; ++i)
                counter += data[i].size();
        vector<vector<T> > res(counter);
        counter = 0;
        for (size_t i = 0; i < out_size; ++i)
        {
                for (size_t j = 0; j < data[i].size(); ++j)
                {
                        res[counter] = data[i][j];
                        ++counter;
                }
        }
        return res;
}

template<class T> vector<T>  CombineVectors(vector<vector<T> > const& data)
{
        vector<T> res;
        for (size_t i = 0; i < data.size(); ++i)
                res.insert(res.end(), data[i].begin(), data[i].end());
        return res;
}

template<class S, class T> T& safe_retrieve
(boost::container::flat_map<S, T>& m,
        const S& s)
{
        typename boost::container::flat_map<S, T>::iterator it =
                m.find(s);
        assert(it != m.end());
        return it->second;
}

double fastsqrt(double x);

#endif // UTILS_HPP