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 kappa4LogPDF(x: Float64, a: Float64, c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale

    var res = 0.0
    if (a == 0.0 && c == 0.0) {
        res = - y - exp(-y)
    } else if (a != 0.0 && c == 0.0) {
        res = - y + (1.0 / a - 1.0) * log(1.0 - a * exp(-y))
    } else if (a == 0.0 && c != 0.0) {
        let k = 1.0 / c
        res = (k - 1.0) * log(1.0 - c * y) - pow(1.0 - c * y, k)
    } else {
        let h = 1.0 / a
        let k = 1.0 / c
        res = (k - 1.0) * log(1.0 - c * y) + (h - 1.0) * log(1.0 - a * pow(1.0 - c * y, k))
    }
    return res - log(scale)
}

/*
 * Probability density function
 */
public func kappa4PDF(x: Float64, a: Float64, c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let res = exp(kappa4LogPDF(x, a, c, loc: loc, scale: scale))
    return res
}

/*
 * Cumulative probability density function
 */
public func kappa4CDF(x: Float64, a: Float64, c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let res = exp(kappa4LogCDF(x, a, c, loc: loc, scale: scale))
    return res
}

/*
 * Log of Cumulative probability density function
 */
public func kappa4LogCDF(x: Float64, a: Float64, c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale
    var res = 0.0

    if (a == 0.0 && c == 0.0) {
        res = - exp(-y)
    } else if (a != 0.0 && c == 0.0) {
        res = (1.0 / a) * log(1.0 - a * exp(-y))
    } else if (a == 0.0 && c != 0.0) {
        let k = 1.0 / c
        res = - pow(1.0 - c * y, k)
    } else {
        let h = 1.0 / a
        let k = 1.0 / c
        res = h * log(1.0 - a * pow(1.0 - c * y, k))
    }

    return res 
}

/*
 * PPF
 */
public func kappa4PPF(q: Float64, a: Float64, c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    if (q <= 0.0 || q >= 1.0) {
        throw IllegalArgumentException("kappa4PPF: quantile out of bound.")
    }
    var res = 0.0

    if (a == 0.0 && c == 0.0) {
        res = -log(-log(q))
    } else if (a != 0.0 && c == 0.0) {
        res = -log(1.0 - pow(q, a)) + log(a)
    } else if (a == 0.0 && c != 0.0) {
        let k = 1.0 / c
        res = k * (1.0 - pow(-log(q), c))
    } else {
        let h = 1.0 / a
        let k = 1.0 / c
        let t1 = (1.0 - pow(q, a)) / a
        let t2 = 1.0 - pow(t1, c)
        res = k * t2
    }

    return res * scale + loc
}

@Test
public class TestKappa4 {
    @TestCase
    func testKappa4(): Unit {
        @Assert(approxEqual(kappa4LogPDF(2.0, 0.1, 0.0, loc: 1.0, scale: 2.0), -1.7562807463841041,  atol:1e-13))
        @Assert(approxEqual(kappa4LogCDF(2.0, 0.1, 0.0, loc: 1.0, scale: 2.0), -0.6257039620268432,  atol:1e-13)) 
        @Assert(approxEqual(kappa4PPF(0.2, 0.1, 0.0, loc: 1.0, scale: 2.0),     0.20701569101413364, atol:1e-6))
    }
}