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

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

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

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

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

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


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

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

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

    return log(temp)
}

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

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

    var t = 0.5
    if (x < 0.0) {
        t = -t
    } else if (x == 0.0) {
        t = 0.0
    }

    let temp = 0.5 + t - t * (1.0 - igam(1.0 / k, pow(abs(y), k)))
    return temp
}


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

    let r1 = gamSample(r, 1.0 / k)
    let t = pow(r1, 1.0 / k)
    var r2 = r.nextFloat64()
    if (r2 < 0.5) {
        r2 = -r2
    }
    return r2 * scale + loc
}


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

    return loc
}


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

    let t1 = gammaLog(1.0 / k)
    let t3 = gammaLog(3.0 / k)
    return exp(t3 - t1) * scale * scale
}

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

    let temp = gennormVar(k ,loc: loc, scale: scale)
    return sqrt(temp)
}


@Test
public class TestGennorm {
    @TestCase
    func testGennorm(): Unit {
        @Assert(approxEqual(gennormLogPDF(2.0, 2.0, loc: 1.0, scale: 2.0), -1.5155121234846454, atol:1e-13))
        @Assert(approxEqual(gennormLogCDF(2.0, 2.0, loc: 1.0, scale: 2.0), -0.2741080327843856, atol:1e-13))
        @Assert(approxEqual(gennormMean(2.0, loc: 1.0, scale: 2.0),         1.0,                atol:1e-13)) 
        @Assert(approxEqual(gennormVar(2.0, loc: 1.0, scale: 2.0),          2.0,                atol:1e-13))  
    }
}