package scientific.stats.continuous

import std.math.*
import std.unittest.*
import std.unittest.testmacro.*

import scientific.numbers.*
import scientific.stats.random.*

/*
 * Log of Probability density function
 */
public func exponpowLogPDF(x: Float64, k: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale

    if (y < 0.0) {
        throw IllegalArgumentException("exponpowLogPDF: input value out of bound.")
    }

    if (k <= 0.0) {
        throw IllegalArgumentException("exponpowLogPDF: shape parameter out of bound.")
    }

    let temp = -exp(pow(y, k))
    return log(k) + (k - 1.0) * log(y) + 1.0 + pow(y, k) + temp - log(scale)
}

/*
 * Probability density function
 */
public func exponpowPDF(x: Float64, k: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale

    if (y < 0.0) {
        throw IllegalArgumentException("exponpowPDF: input value out of bound.")
    }

    if (k <= 0.0) {
        throw IllegalArgumentException("exponpowPDF: shape parameter out of bound.")
    }

    return exp(exponpowLogPDF(x, k, loc: loc, scale: scale))
}

/*
 * Cumulative probability density function
 */
public func exponpowCDF(x: Float64, k: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale

    if (y < 0.0) {
        throw IllegalArgumentException("exponpowCDF: input value x out of bound.")
    }

    if (k <= 0.0) {
        throw IllegalArgumentException("exponpowCDF: shape parameter out of bound.")
    }

    let temp = -exp(pow(y, k)) + 1.0

    return -exp(temp) + 1.0
}

/*
 * Log of Cumulative probability density function
 */
public func exponpowLogCDF(x: Float64, k: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale

    if (y < 0.0) {
        throw IllegalArgumentException("exponpowLogCDF: input value x out of bound.")
    }

    if (k <= 0.0) {
        throw IllegalArgumentException("exponpowLogCDF: shape parameter out of bound.")
    }

    let temp = exponpowCDF(x, k, loc: loc, scale: scale)
    if (temp < 0.000001) {
        throw IllegalArgumentException("exponpowLogCDF: return-value too small.")
    }

    return log(temp) 
}


/*
 * ppf
 */
public func exponpowPPF(q: Float64, k: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    if (q <= 0.0 || q >= 1.0) {
        throw IllegalArgumentException("exponpowPPF: quantile out of bound.")
    }

    if (k <= 0.0) {
        throw IllegalArgumentException("exponpowPPF: shape parameter out of bound.")
    }

    let temp = -log(1.0 - q) + 1.0

    return pow(log(temp), 1.0 / k) * scale + loc
}


@Test
public class TestExponpow {
    @TestCase
    func testExponpow(): Unit {
        @Assert(approxEqual(exponpowLogPDF(2.0, 2.0, loc: 1.0, scale: 2.0), -0.7271725972476867, atol:1e-13))
        @Assert(approxEqual(exponpowLogCDF(2.0, 2.0, loc: 1.0, scale: 2.0), -1.3973452463038125, atol:1e-13))
        @Assert(approxEqual(exponpowPPF(0.2, 2.0, loc: 1.0, scale: 2.0),     1.8976062080849678, atol:1e-13))
    }
}