/*
 * Toolbox.hh
 *
 *  Created on: Jun 23, 2018
 *      Author: benjamin
 */

#ifndef TOOLBOX_HH_
#define TOOLBOX_HH_

#include "ParamData.hh"

namespace AMDA {
    namespace Parameters {
        namespace StatisticFunctions {

            template <typename Type>
            Type min(Type a, Type b) {
                if (!std::isfinite(a) || isNAN(a))
                    return b;
                if (!std::isfinite(b) || isNAN(b))
                    return a;
                return std::min(a, b);
            }

            template <typename Type>
            std::vector<Type> min(std::vector<Type> a, std::vector<Type> b) {
                std::vector<Type> res = a;
                if (a.size() != b.size()) {
                    BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_OPER_NOT_ALLOWED) << AMDA::ex_msg("Vectors don't have the same size"));
                }
                for (unsigned int i = 0; i < a.size(); ++i) {
                    res[i] = StatisticFunctions::min(a[i], b[i]);
                }
                return res;
            }

            template <typename Type>
            Tab2DData<Type> min(Tab2DData<Type> a, Tab2DData<Type> b) {
                Tab2DData<Type> res(a);
                if (a.getDim1Size() != b.getDim1Size() || a.getDim2Size() != b.getDim2Size()) {
                    BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_OPER_NOT_ALLOWED) << AMDA::ex_msg("Tab2DData don't have the same dimensions"));
                }
                for (int i = 0; i < a.getDim1Size(); ++i) {
                    for (int j = 0; j < a.getDim2Size(); ++j) {
                        res[i][j] = StatisticFunctions::min(a[i][j], b[i][j]);
                    }
                }
                return res;
            }

            template <typename Type>
            Type max(Type a, Type b) {
                if (!std::isfinite(a) || isNAN(a))
                    return b;
                if (!std::isfinite(b) || isNAN(b))
                    return a;
                return std::max(a, b);
            }

            template <typename Type>
            std::vector<Type> max(std::vector<Type> a, std::vector<Type> b) {
                std::vector<Type> res = a;
                if (a.size() != b.size()) {
                    BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_OPER_NOT_ALLOWED) << AMDA::ex_msg("Vectors don't have the same size"));
                }
                for (unsigned int i = 0; i < a.size(); ++i) {
                    res[i] = StatisticFunctions::max(a[i], b[i]);
                }
                return res;
            }

            template <typename Type>
            Tab2DData<Type> max(Tab2DData<Type> a, Tab2DData<Type> b) {
                Tab2DData<Type> res(a);
                if (a.getDim1Size() != b.getDim1Size() || a.getDim2Size() != b.getDim2Size()) {
                    BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_OPER_NOT_ALLOWED) << AMDA::ex_msg("Tab2DData don't have the same dimensions"));
                }
                for (int i = 0; i < a.getDim1Size(); ++i) {
                    for (int j = 0; j < a.getDim2Size(); ++j) {
                        res[i][j] = StatisticFunctions::max(a[i][j], b[i][j]);
                    }
                }
                return res;
            }

            template <typename Type>
            Type pow(Type a, double exp) {
                return std::pow(a, exp);
            }

            template <typename Type>
            std::vector<Type> pow(std::vector<Type> a, double exp) {
                std::vector<Type> res = a;
                for (unsigned int i = 0; i < a.size(); ++i) {
                    res[i] = StatisticFunctions::pow(a[i], exp);
                }
                return res;
            }

            template <typename Type>
            Tab2DData<Type> pow(Tab2DData<Type> a, double exp) {
                Tab2DData<Type> res = a;
                for (int i = 0; i < a.getDim1Size(); ++i) {
                    for (int j = 0; j < a.getDim2Size(); ++j) {
                        res[i][j] = StatisticFunctions::pow(a[i][j], exp);
                    }
                }
                return res;
            }

            template <typename Type>
            Type square(Type a) {
                return StatisticFunctions::pow(a, 2);
            }

            template <typename Type>
            std::vector<Type> square(std::vector<Type> a) {
                return StatisticFunctions::pow(a, 2);
            }

            template <typename Type>
            Tab2DData<Type> square(Tab2DData<Type> a) {
                return StatisticFunctions::pow(a, 2);
            }

            template <typename Type>
            Type root_square(Type a) {
                if (isNAN(a) || a < 0) {
                    Type res;
                    res << NotANumber();
                    return res;
                }
                return std::sqrt(a);
            }

            template <typename Type>
            std::vector<Type> root_square(std::vector<Type> a) {
                std::vector<Type> res = a;
                for (unsigned int i = 0; i < a.size(); ++i) {
                    res[i] = StatisticFunctions::root_square(a[i]);
                }
                return res;
            }

            template <typename Type>
            Tab2DData<Type> root_square(Tab2DData<Type> a) {
                Tab2DData<Type> res = a;
                for (int i = 0; i < a.getDim1Size(); ++i) {
                    for (int j = 0; j < a.getDim2Size(); ++j) {
                        res[i][j] = StatisticFunctions::root_square(a[i][j]);
                    }
                }
                return res;
            }

            template <typename Type>
            Type div(Type a, Type b) {
                if (isNAN(a) || isNAN(b) || (b == 0)) {
                    Type res;
                    res << NotANumber();
                    return res;
                }
                return a / b;
            }

            template <typename Type>
            std::vector<Type> div(std::vector<Type> a, std::vector<Type> b) {
                std::vector<Type> res = a;
                if (a.size() != b.size()) {
                    BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_OPER_NOT_ALLOWED) << AMDA::ex_msg("Vectors don't have the same size"));
                }
                for (unsigned int i = 0; i < a.size(); ++i) {
                    res[i] = StatisticFunctions::div(a[i], b[i]);
                }
                return res;
            }

            template <typename Type>
            Tab2DData<Type> div(Tab2DData<Type> a, Tab2DData<Type> b) {
                Tab2DData<Type> res(a);
                if (a.getDim1Size() != b.getDim1Size() || a.getDim2Size() != b.getDim2Size()) {
                    BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_OPER_NOT_ALLOWED) << AMDA::ex_msg("Tab2DData don't have the same dimensions"));
                }
                for (int i = 0; i < a.getDim1Size(); ++i) {
                    for (int j = 0; j < a.getDim2Size(); ++j) {
                        res[i][j] = StatisticFunctions::div(a[i][j], b[i][j]);
                    }
                }
                return res;
            }

            template <typename Type>
            std::pair<Type, Type> getMean(std::list<std::pair<Type, Type>> &list) {
                std::pair<Type, Type> result(0, 0);
                std::pair<int, int> counter(0, 0);

                for (auto elt : list) {
                    if (!isNAN(elt.first)) {
                        result.first += elt.first;
                        counter.first += 1;
                    }
                    if (!isNAN(elt.second)) {
                        result.second += elt.second;
                        counter.second += 1;
                    }
                }
                if (counter.first != 0)
                    result.first /= counter.first;

                if (counter.second != 0)
                    result.second /= counter.second;

                return result;
            }

            template <typename Type>
            std::pair<Type, Type> getStd(std::list<std::pair<Type, Type>> &list) {
                std::pair<Type, Type> mean = getMean(list);
                std::pair<Type, Type> result(0, 0);
                int counter1 = 0;
                int counter2 = 0;
                for (auto elt : list) {
                    if (!isNAN(elt->first)) {
                        result.first += (elt->first - mean.first)*(elt->first - mean.first);
                        counter1 += 1;
                    }
                    if (!isNAN(elt->second)) {
                        result.second += (elt->second - mean.second)*(elt->second - mean.second);
                        counter2 += 1;
                    }
                }
                if (counter1 != 0) {
                    result.first /= counter1;
                } else {
                    result.first << NotANumber();
                }
                if (counter2 != 0) {
                    result.second /= counter2;
                } else {
                    result.second << NotANumber();
                }
                return std::sqrt(result);
            }

            template <typename Type>
            std::vector<Type>rankify(std::vector<Type>& X) {

                int N = X.size();

                // Rank Vector
                std::vector<Type> Rank_X(N);

                for (int i = 0; i < N; i++) {
                    int r = 1, s = 1;

                    // Count no of smaller elements
                    // in 0 to i-1
                    for (int j = 0; j < N; j++) {
                        if (X[j] < X[i]) r++;
                        if (X[j] == X[i] && i != j) s++;
                    }

                    // Count no of smaller elements
                    // in i+1 to N-1
                    /**
                    for (int j = i + 1; j < N; j++) {
                        if (X[j] < X[i]) r++;
                        if (X[j] == X[i]) s++;
                    }*/

                    // Use Fractional Rank formula
                    // fractional_rank = r + (n-1)/2
                    Rank_X[i] = (Type) r + (s - 1) * 0.5;
                }

                // Return Rank Vector
                return Rank_X;
            }

            template <typename Type>
            bool getCovariance(std::list<std::pair<Type, Type>> &list, Type &result) {
                if (list.empty()) {
                    return false;
                }
                std::pair<Type, Type> mean = getMean(list);
                result = 0;
                int counter = 0;
                for (auto elt : list) {
                    if (!isNAN(elt.first) && !isNAN(elt.second)) {
                        result += (elt.first - mean.first)*(elt.second - mean.second);
                        counter += 1;
                    }
                }
                if (counter != 0)
                    result /= counter;
                return true;
            }

            template <typename Type>
            bool getPearson(std::list<std::pair<Type, Type>> &list, Type &result) {
                if (list.empty()) {
                    return false;
                }
                result = 0;
                int counter = 0;
                Type sum1 = 0, sum2 = 0;
                Type sum1Sq = 0, sum2Sq = 0;
                Type sum12 = 0;
                for (auto elt : list) {
                    if (!isNAN(elt.first) && !isNAN(elt.second)) {
                        sum1 += elt.first;
                        sum2 += elt.second;
                        sum1Sq += elt.first * elt.first;
                        sum2Sq += elt.second * elt.second;
                        sum12 += elt.first * elt.second;
                        counter += 1;
                    }
                }
                if (counter > 0)
                    result = (counter * sum12 - sum1 * sum2) /
                    (std::sqrt((counter * sum1Sq - sum1 * sum1)*(counter * sum2Sq - sum2 * sum2)));
                return true;
            }

            template <typename Type>
            bool getSpearman(std::list<std::pair<Type, Type>> &list, Type &result) {
                if (list.empty()) {
                    return false;
                }

                std::vector<Type> X, Y;
                for (auto elt : list) {
                    if (!isNAN(elt.first) && !isNAN(elt.second)) {
                        X.push_back(elt.first);
                        Y.push_back(elt.second);
                    }
                }
                int n = X.size();
                if (n == 0)
                    return true;
                std::vector<Type> rank_X = rankify(X);
                std::vector<Type> rank_Y = rankify(Y);
                std::list<std::pair<Type, Type>> rankList;
                
                for (int i = 0; i < n; i++)
                    rankList.push_back(std::make_pair(rank_X[i], rank_Y[i]));
                
                result = 0;
                getPearson(rankList, result);
                return true;
            }

            template <typename Type>
            bool getKendall(std::list<std::pair<Type, Type>> &list, Type &result) {
                if (list.empty()) {
                    return false;
                }
                result = 0;
                std::vector<Type> X, Y;
                for (auto elt : list) {
                    if (!isNAN(elt.first) && !isNAN(elt.second)) {
                        X.push_back(elt.first);
                        Y.push_back(elt.second);
                    }
                }
                int n = X.size();
                if (n == 0)
                    return true;
                std::vector<Type> rank_x = rankify(X);
                std::vector<Type> rank_y = rankify(Y);
                Type nc = 0;
                for (int i = 0; i < n; i++){
                    for (int j=i;j<n;j++){
                        if((i != j )&& ((rank_x[i] >  rank_x[j] &&  rank_y[i] >  rank_y[j]) || 
                                (rank_x[i]  < rank_x[j] &&  rank_y[i] <  rank_y[j])))
                            nc += 1;
                    }
                    
                }
                result = (Type) (4*nc - n*(n-1))/(Type) (n * (n - 1));

                return true;
            }

        } /* namespace StatisticFunctions */
    } /* namespace Parameters */
} /* namespace AMDA */
#endif /* TOOLBOX_HH_ */