package scientific.stats.continuous

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

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

foreign func psi(t: Float64): Float64

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

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

    return log(k) - y - (k + 1.0) * log(1.0 + exp(-y)) - log(scale)
}

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

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

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


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

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

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

    return log(temp)
}

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

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

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


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

    let ppf = -log(pow(q, -1.0 / k) - 1.0)

    return ppf * scale + loc
}


@Test
public class TestGenlogistic {
    @TestCase
    func testGenlogistic(): Unit {
        @Assert(approxEqual(genlogisticLogPDF(2.0, 2.0, loc: 1.0, scale: 2.0), -1.9222309525403203, atol:1e-13))
        @Assert(approxEqual(genlogisticLogCDF(2.0, 2.0, loc: 1.0, scale: 2.0), -0.9481539683602133, atol:1e-13))
        @Assert(approxEqual(genlogisticPPF(0.2, 2.0, loc: 1.0, scale: 2.0),     0.5761292889993161, atol:1e-13))
    }
}