#ifndef VTR_MATH_H
#define VTR_MATH_H

#include <map>
#include <cmath>

#include "vtr_assert.h"

namespace vtr {
    /*********************** Math operations *************************************/
    int ipow(int base, int exp);

    template<typename X, typename Y>
    Y linear_interpolate_or_extrapolate(const std::map<X,Y> *xy_map, X requested_x);

    constexpr int nint(float val) { return static_cast<int>(val + 0.5); }

    //Returns a 'safe' ratio which evaluates to zero if the denominator is zero
    template<typename T>
    T safe_ratio(T numerator, T denominator) {
        if (denominator == T(0)) {
            return 0;
        }
        return numerator / denominator;
    }

    template<typename InputIterator>
    double median(InputIterator first, InputIterator last) {
        auto len = std::distance(first, last);
        auto iter = first + len / 2;

        if (len % 2 == 0) {
            return (*iter + *(iter + 1)) / 2;        
        } else {
            return *iter;        
        }
    }

    template<typename Container>
    double median(Container c) {
        return median(std::begin(c), std::end(c));
    }


    template<typename InputIterator>
    double geomean(InputIterator first, InputIterator last, double init=1.) {
        //Compute the geometric mean of the elments in range [first, last)
        //
        //To avoid potential round-off issues we transform the standard formula:
        //
        //      geomean = ( v_1 * v_2 * ... * v_n) ^ (1/n)
        //
        //by taking the log:
        //
        //      geomean = exp( (1 / n) * (log(v_1) + log(v_2) + ... + log(v_n)))

        double log_sum = std::log(init);
        size_t n = 0;
        for(auto iter = first; iter != last; ++iter) {
            log_sum += std::log(*iter);
            n += 1;
        }

        VTR_ASSERT(n > 0.);

        return std::exp( (1. / n) * log_sum );
    }

    template<typename Container>
    double geomean(Container c) {
        return geomean(std::begin(c), std::end(c));
    }


    template<typename InputIterator>
    double arithmean(InputIterator first, InputIterator last, double init=0.) {
        double sum = init;
        size_t n = 0;
        for (auto iter = first; iter != last; ++iter) {
            sum += *iter;
            n += 1;
        }

        VTR_ASSERT(n > 0.);
        return sum / n;
    }

    template<typename Container>
    double arithmean(Container c) {
        return arithmean(std::begin(c), std::end(c));
    }


    //Return the greatest common divisor of x and y
    // Note that T should be an integral type
    template<typename T>
    static T gcd(T x, T y){
        static_assert(std::is_integral<T>::value, "T must be integral");
        //Euclidean algorithm
        if (y == 0){
            return x;
        }
        return gcd(y, x % y);
    }

    //Return the least common multiple of x and y
    // Note that T should be an integral type
    template<typename T>
    T lcm(T x, T y) {
        static_assert(std::is_integral<T>::value, "T must be integral");

        if (x == 0 && y == 0) {
            return 0;
        } else {
            return (x / gcd(x,y)) * y;
        }
    }

    constexpr double DEFAULT_REL_TOL = 1e-9;
    constexpr double DEFAULT_ABS_TOL = 0;

    template<class T>
    bool isclose(T a, T b, T rel_tol, T abs_tol) {

        if (std::isinf(a) && std::isinf(b)) return (std::signbit(a) == std::signbit(b));
        if (std::isnan(a) && std::isnan(b)) return false;

        T abs_largest = std::max(std::abs(a), std::abs(b));
        return std::abs(a - b) <= std::max(rel_tol * abs_largest, abs_tol);
    }

    template<class T>
    bool isclose(T a, T b) {
        return isclose<T>(a, b, DEFAULT_REL_TOL, DEFAULT_ABS_TOL);
    }

}

#endif