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

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

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

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

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

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

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

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

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

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

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

    let res = 1.0 - igam(k, 1.0 / y)
    return res
}

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

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

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

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

    return log(temp) 
}


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

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

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

/*
 * Mean
 */
public func invgammaMean(k: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    var res = Float64.Inf
    if (k > 1.0) {
        res = 1.0 / (k - 1.0)
        return res * scale + loc
    }
    return res
}

/*
 * Var
 */
public func invgammaVar(k: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    var res = Float64.Inf
    if (k > 2.0) {
        res = 1.0 / ((k - 1.0) * (k - 1.0) * (k - 2.0))
        return res * scale * scale
    }
    return res
}

/*
 * Std
 */
public func invgammaStd(k: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    if (k <= 2.0) {
        return Float64.Inf
    }

    let temp = invgammaVar(k, loc: loc, scale: scale)

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

    return sqrt(temp)
}

@Test
public class TestInvgamma {
    @TestCase
    func testInvgamma(): Unit {
        @Assert(approxEqual(invgammaLogPDF(2.0, 2.0, loc: 1.0, scale: 2.0), -0.6137056388801095, atol:1e-13))
        @Assert(approxEqual(invgammaLogCDF(2.0, 2.0, loc: 1.0, scale: 2.0), -0.9013877113318907, atol:1e-13))
        @Assert(approxEqual(invgammaPPF(0.2, 2.0, loc: 1.0, scale: 2.0),     1.6679338826284824, atol:1e-13))
        @Assert(approxEqual(invgammaMean(3.0, loc: 1.0, scale: 2.0),         2.0,                atol:1e-13))
        @Assert(approxEqual(invgammaVar(3.0, loc: 1.0, scale: 2.0),          1.0,                atol:1e-13))
    }
}