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

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

    return log(k) - gammaLog(1.0 / k) - pow(abs(y), k) - log(scale)
}

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

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

    let temp = halfgennormLogPDF(x, k, loc:loc, scale: scale)
    return exp(temp)
}

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

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

    return igam(1.0 / k, pow(y, k))
}

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

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

    let temp = halfgennormCDF(x, k, loc: loc, scale: scale)

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

    return log(temp) 
}


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

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


@Test
public class TsetHalfgennorm {
    @TestCase
    func testHalfgennorm(): Unit {
        @Assert(approxEqual(halfgennormLogPDF(2.0, 2.0, loc: 1.0, scale: 2.0), -0.8223649429247,    atol:1e-13))
        @Assert(approxEqual(halfgennormLogCDF(2.0, 2.0, loc: 1.0, scale: 2.0), -0.652965625676331,  atol:1e-13))
        @Assert(approxEqual(halfgennormPPF(0.2, 2.0, loc: 1.0, scale: 2.0),     1.3582869092425833, atol:1e-13))
    }
}