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

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

    if (y < 0.0 || y > (1.0/k)) {
        throw IllegalArgumentException("genhalflogisticLogPDF: input value out of bound.")
    }

    let t = 1.0 - k * y
    return log(2.0) + (1.0/k - 1.0) * log(t) - 2.0 * log(1.0 + pow(t, 1.0/k)) - log(scale)
}

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

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

    if (y < 0.0 || y > (1.0/k)) {
        throw IllegalArgumentException("genhalflogisticPDF: input value out of bound.")
    }

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

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

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

    if (y < 0.0 || y > (1.0/k)) {
        throw IllegalArgumentException("genhalflogisticCDF: input value out of bound.")
    }

    let val1 = 1.0 / k
    let val2 = 1.0 - k * y
    let temp = pow(val2, val1)
    return (1.0 - temp) / (1.0 + temp)
}

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

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

    if (y < 0.0 || y > (1.0/k)) {
        throw IllegalArgumentException("genhalflogisticLogCDF: input value out of bound.")
    }

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

    return log(temp) 
}


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

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

    let t = (1.0 - q) / (1.0 + q)
    let res = 1.0 / k * (1.0 - pow(t, k))
    return res * scale + loc
}


@Test
public class TestGenhalflogistic {
    @TestCase
    func testGenhalflogistic(): Unit {
        @Assert(approxEqual(genhalflogisticLogPDF(1.8, 2.0, loc: 1.0, scale: 2.0),  0.06543885741205324, atol:1e-13))
        @Assert(approxEqual(genhalflogisticLogCDF(1.8, 2.0, loc: 1.0, scale: 2.0), -0.9624236501192067,  atol:1e-13)) 
        @Assert(approxEqual(genhalflogisticPPF(0.2, 2.0, loc: 1.0, scale: 2.0),     1.5555555555555554,  atol:1e-13)) 
    }
}