/* ----------------------------------------------------------------------- */
/**
 *
 * @file kolmogorov.cpp
 *
 * @brief Kolmogorov Probability distribution function
 *
 *
 */
/* ----------------------------------------------------------------------- */

#include <dbconnector/dbconnector.hpp>

#include "kolmogorov.hpp"

namespace madlib {
namespace modules {
namespace prob {

/**
 * @brief Komogorov cumulative distribution function: In-database interface
 */
AnyType
kolmogorov_cdf::run(AnyType &args) {
    return prob::cdf(kolmogorov(), args[0].getAs<double>());
}

double KolmogorovProb(double z);

/**
 * @brief  Calculates the Kolmogorov distribution function
 * @param  x Value to compare with
 * @return   Output probability that Kolmogorov's test statistic will exceed
 *            the value x assuming the null hypothesis.
 *
 *    The implementation is based on the method described in
 *    http://www.jstatsoft.org/v08/i18/paper - Section 3
 *    The limiting forms in the paper correspond to P(K <= x); since we need
 *    the P(K > x) we return 1 - P(K <= x).
 */
double KolmogorovProb(double x) {
   double u = std::fabs(x);
   double p = 0;
   if (u < 0.1) {
        // probability too close to 1
        p = 1;
    } else if (u <= 1) {
        // for small u we use the theta function inversion formula
        // p = 1 - (sqrt(2*pi)/u) * sum e^( -((2i-1) pi)^2 / (8u^2) )
        const double v = 1./(8*u*u);
        const double k1 = -v * pow(M_PI, 2);
        const double k2 = 9 * k1;
        const double k3 = 25 * k1;
        const double k4 = 49 * k1;
        const double w = sqrt(2 * M_PI) / u;
        p = 1 - w * (std::exp(k1) + std::exp(k2) + std::exp(k3) + std::exp(k4));
    } else if (u < 5) {
        // p = 2 * sum -1^(i-1) * e^( -2 * i^2 * u^2) )
        double prob_sum = 0;
        double sign = 1;
        double v = u*u;
        uint16_t maxi = std::max(uint16_t(1), uint16_t(4./u + 0.5));
        // since u is > 1, 4/u will be < 4
        for (uint16_t i=1; i <= maxi && i <= 4; i++) {
            prob_sum += sign * std::exp(-2 * (i*i) * v);
            sign = -sign;  //
        }
        p = 2 * prob_sum;
    } else {
        p = 0;
    }
   return p;
}

} // namespace prob
} // namespace modules
} // namespace regress