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

    if (k > 0.0) {
        if (x > 1.0) {
            throw IllegalArgumentException("genextremeLogPDF: input value out of bound.")
        }
    } else if (k < 0.0) {
        if (x < 1.0 / k) {
            throw IllegalArgumentException("genextremeLogPDF: input value out of bound.")
        }
    } else {
        return gumbel_rLogPDF(x, loc: loc, scale: scale)
    }

    let temp = 1.0 - k * y
    return -pow(temp, 1.0/k) + (1.0/k - 1.0) * log(temp) - log(scale)
}

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

    if (k > 0.0) {
        if (x > 1.0) {
            throw IllegalArgumentException("genextremePDF: input value out of bound.")
        }
    } else if (k < 0.0) {
        if (x < 1.0 / k) {
            throw IllegalArgumentException("genextremePDF: input value out of bound.")
        }
    } else {
        return gumbel_rPDF(x, loc: loc, scale: scale)
    }

    return exp(genextremeLogPDF(x, k, loc: loc, scale: scale))

}

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

    if (k > 0.0) {
        if (x > 1.0) {
            throw IllegalArgumentException("genextremeCDF: input value out of bound.")
        }
    } else if (k < 0.0) {
        if (x < 1.0 / k) {
            throw IllegalArgumentException("genextremeCDF: input value out of bound.")
        }
    } else {
        return gumbel_rCDF(x, loc: loc, scale: scale)
    }

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

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

    if (k > 0.0) {
        if (x > 1.0) {
            throw IllegalArgumentException("genextremeLogCDF: input value out of bound.")
        }
    } else if (k < 0.0) {
        if (x < 1.0 / k) {
            throw IllegalArgumentException("genextremeLogCDF: input value out of bound.")
        }
    } else {
        return gumbel_rLogCDF(x, loc: loc, scale: scale)
    }

    let temp = log(1.0 - k * y) / k
    return -exp(temp)
}


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

    if (k == 0.0) {
        return gumbel_rPPF(q, loc: loc, scale: scale)
    }

    let t = -log(-log(q))
    let res = (1.0 - exp(-k * t)) / k
    return res * scale + loc
}


@Test
public class TestGenextreme {
    @TestCase
    func testGenextreme(): Unit {
        @Assert(approxEqual(genextremeLogPDF(0.0, 2.0, loc: 1.0, scale: 2.0), -2.453934333213013,  atol:1e-13))
        @Assert(approxEqual(genextremeLogCDF(0.0, 2.0, loc: 1.0, scale: 2.0), -1.414213562373095,  atol:1e-13)) 
        @Assert(approxEqual(genextremePPF(0.2, 2.0, loc: 1.0, scale: 2.0),    -0.5902903939802346, atol:1e-13)) 
    }
}