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

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

    let res = log(2.0) - y - 2.0 * log(1.0 + exp(-y))

    return res - log(scale)
}

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

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

    let temp = halflogisticLogPDF(x, loc: loc, scale: scale)

    return exp(temp)
}


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

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

    return tanh(0.5 * y)
}


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

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

    let temp = halflogisticCDF(x, loc: loc, scale: scale)

    if (temp < 0.000001) {
        throw IllegalArgumentException("halflogisticLogCDF: return-value too small.")
    }

    return log(temp)
}


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

    let res = 2.0 * atanh(q)
    return res * scale + loc
}


@Test
public class TestHalfLogistic {
    @TestCase
    func testHalflogistic(): Unit {
        @Assert(approxEqual(halflogisticLogPDF(3.0, loc: 2.0, scale: 1.0), -0.9333761944765004, atol:1e-13))
        @Assert(approxEqual(halflogisticLogCDF(3.0, loc: 2.0, scale: 1.0), -0.7719368329053048, atol:1e-13))
        @Assert(approxEqual(halflogisticPPF(0.7, loc: 2.0, scale: 1.0),     3.734601055388106,  atol:1e-13))
    }
}