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

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

    let temp = abs(y)
    return  (a - 1.0) * log(temp) - log(2.0) - gammaLog(a) - temp - log(scale)
}

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

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

    return exp(dgammaLogPDF(x, a, loc: loc, scale: scale))
}

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

    let r1 = r.nextFloat64()
    let r2 = gamSample(r, a)

    var res = 0.0
    if (r1 >= 0.5) {
        res = r2
    } else {
        res = -r2
    }
    return res * scale + loc
}


/*
 * Mean
 */
public func dgammaMean(a: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    if (a <= 0.0) {
        throw IllegalArgumentException("dgammaMean: shape parameter out of bound.")
    }

    return loc
}


/*
 * Var
 */
public func dgammaVar(a: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    if (a <= 0.0) {
        throw IllegalArgumentException("dgammaVar: shape parameter out of bound.")
    }

    let temp = a * (a + 1.0)
    return temp * scale * scale
}


@Test
public class TestDgamma {
    @TestCase
    func testDgamma(): Unit {
        @Assert(approxEqual(dgammaLogPDF(2.0, 2.0, loc: 1.0, scale: 2.0), -2.579441541679836, atol:1e-13))
        @Assert(approxEqual(dgammaMean(2.0, loc: 1.0, scale: 2.0),         1.0,               atol:1e-13)) 
        @Assert(approxEqual(dgammaVar(2.0, loc: 1.0, scale: 2.0),          24.0,              atol:1e-13)) 
    }
}