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

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

    if (a <= 0.0 || c == 0.0) {
        throw IllegalArgumentException("gengammaLogPDF: shape parameter out of bound.")
    }

    return log(abs(c)) + (a * c - 1.0) * log(y) - pow(y, c) - gammaLog(a) - log(scale)
}

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

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

    if (a <= 0.0 || c == 0.0) {
        throw IllegalArgumentException("gengammaPDF: shape parameter out of bound.")
    }

    return exp(gengammaLogPDF(x, a, c, loc: loc, scale: scale))
}

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

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

    if (a <= 0.0 || c == 0.0) {
        throw IllegalArgumentException("gengammaCDF: shape parameter out of bound.")
    }
    
    let t = pow(y, c)
    let val1 = igam(a, t)
    let val2 = 1.0 - igam(a, t)
    var res = 0.0
    if (c > 0.0) {
        res = val1
    } else {
        res = val2
    }
    return res
}

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

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

    if (a <= 0.0 || c == 0.0) {
        throw IllegalArgumentException("gengammaLogCDF: shape parameter out of bound.")
    }

    let temp = gengammaCDF(x, a, c, loc: loc, scale: scale)

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

    return log(temp)
}

/*
 * Sample
 */
public func gengammaSample(r: Random, a: Float64, c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    if (a <= 0.0 || c == 0.0) {
        throw IllegalArgumentException("gengammaSample: shape parameter out of bound.")
    }

    let r1 = gamSample(r, a)
    let r2 = pow(r1, 1.0 / c)
    return r2 * scale + loc
}

@Test
public class testGengamma {
    @TestCase
    func testGengamma(): Unit {
        @Assert(approxEqual(gengammaLogPDF(2.0, 2.0, 3.0, loc: 1.0, scale: 2.0), -3.1852707946915624, atol:1e-13))
        @Assert(approxEqual(gengammaLogCDF(2.0, 2.0, 3.0, loc: 1.0, scale: 2.0), -4.934927177496508,  atol:1e-13))
    }
}