/*
 * 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<double, double> getMean(std::list<std::pair<Type, Type>> &list) {
                std::pair<double, double> result(0, 0);
                std::pair<int, int> counter(0, 0);

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

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

                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<double>rankify(std::vector<Type>& X) {
                std::vector<double> Rank_X;
                if (X.empty()) {
                    return Rank_X;
                }
                Rank_X.resize(X.size());
                int i = 0;
                for (typename std::vector<Type>::iterator it1 = X.begin(); it1 != X.end(); ++it1) {
                    int r = 1, s = 1;

                    for (typename std::vector<Type>::iterator it2 = X.begin(); it2 != X.end(); ++it2) {
                        if (it1 == it2) {
                            continue;
                        }
                        if (*it2 < *it1) r++;
                        if (*it1 == *it2) s++;
                    }

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

                // Return Rank Vector
                return Rank_X;
            }

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

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

                if (list.empty()) {
                    return false;
                }

                std::vector<double> X, Y;
                for (auto elt : list) {
                    if (!isNAN(elt.first) && !isNAN(elt.second)) {
                        X.push_back((double)elt.first);
                        Y.push_back((double)elt.second);
                    }
                }

                if (X.empty()) {
                    return false;
                }

                std::vector<double> rank_X = rankify(X);
                std::vector<double> rank_Y = rankify(Y);
                std::list<std::pair<double, double>> rankList;
                
                for (int i = 0; i < X.size(); i++)
                    rankList.push_back(std::make_pair(rank_X[i], rank_Y[i]));
                
                return getPearson(rankList, result);
            }

            template <typename Type>
            bool getKendall(std::list<std::pair<Type, Type>> &list, double &result) {
                result = NAN;

                if (list.empty()) {
                    return false;
                }

                std::vector<double> X, Y;
                for (auto elt : list) {
                    if (!isNAN(elt.first) && !isNAN(elt.second)) {
                        X.push_back((double)elt.first);
                        Y.push_back((double)elt.second);
                    }
                }

                if (X.empty()) {
                    return false;
                }
                
                std::vector<double> rank_x = rankify(X);
                std::vector<double> rank_y = rankify(Y);

                long nc = 0;
                long nd = 0;
                for (int i = 0; i < rank_x.size()-1; i++){
                    for (int j = i+1; j<rank_x.size(); j++){
                        if (((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;
                        if (((rank_x[i] > rank_x[j]) && (rank_y[i] <  rank_y[j])) || 
                            ((rank_x[i] < rank_x[j]) && (rank_y[i] >  rank_y[j])))
                            ++nd;

                    }
                    
                }
                result = ((double)(nc-nd)/(0.5*rank_x.size()*(rank_x.size()-1)));
                return true;
            }

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