package scientific.stats.continuous

import std.math.*
import std.unittest.*
import std.unittest.testmacro.*

import scientific.numbers.*
import scientific.stats.random.*

public let EULER = 0.577215664901532860606512090082402431042

/*
 * Log of Probability density function
 */
public func gumbel_rLogPDF(x: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale
    return -y - exp(-y) - log(scale)
}

/*
 * Probability density function
 */
public func gumbel_rPDF(x: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale
    return exp(gumbel_rLogPDF(x, loc: loc, scale: scale))
}

/*
 * Cumulative probability density function
 */
public func gumbel_rCDF(x: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale
    return exp(gumbel_rLogCDF(x, loc: loc, scale: scale))
}

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

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

    let temp = -log(q)
    return -log(temp) * scale + loc
}

/*
 * Mean
 */
public func gumbel_rMean(loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    return EULER * scale + loc
}

/*
 * Var
 */
public func gumbel_rVar(loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    return Float64.getPI() * Float64.getPI() / 6.0 * scale * scale
}

/*
 * Std
 */
public func gumbel_rStd(loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let temp = gumbel_rVar(loc: loc, scale: scale)

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

    return sqrt(temp)
}

@Test
public class TestGumbelR {
    @TestCase
    func testGumbel_r(): Unit {
        @Assert(approxEqual(gumbel_rLogPDF(2.0, loc: 1.0, scale: 2.0), -1.7996778402725786,   atol:1e-13))
        @Assert(approxEqual(gumbel_rLogCDF(2.0, loc: 1.0, scale: 2.0), -0.6065306597126334,   atol:1e-13))
        @Assert(approxEqual(gumbel_rPPF(0.2, loc: 1.0, scale: 2.0),     0.048230009345778924, atol:1e-13))
        @Assert(approxEqual(gumbel_rMean(loc: 1.0, scale: 2.0),         2.1544313298030655,   atol:1e-13))
        @Assert(approxEqual(gumbel_rVar(loc: 1.0, scale: 2.0),          6.579736267392906,    atol:1e-13))
    }
}