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

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

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

    return log(a) + log(c) + (a  - 1.0) * log(1.0 - exp(- pow(y, c))) - pow(y, c) + (c - 1.0) * log(y) - log(scale)
}

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

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

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

    return exp(exponweibLogPDF(x, a, c, loc: loc, scale: scale))
}

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

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

    if (a <= 0.0 || c <= 0.0) {
        throw IllegalArgumentException("exponweibCDF: shape parameter out of bound.")
    }
    
    return exp(exponweibLogCDF(x, a, c, loc: loc, scale: scale))
}

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

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

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

    return a * log(1.0 - exp(- pow(y, c))) 
}


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

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

    let temp = - log(1.0 - pow(q, 1.0 / a))

    return pow(temp, 1.0 / c) * scale + loc
}

@Test
public class TestExponweib {
    @TestCase
    func testExponweib(): Unit {
        @Assert(approxEqual(exponweibLogPDF(2.0, 2.0, 2.0, loc: 1.0, scale: 2.0), -1.758691549446032,  atol:1e-13))
        @Assert(approxEqual(exponweibLogCDF(2.0, 2.0, 2.0, loc: 1.0, scale: 2.0), -3.017383098892064,  atol:1e-13))
        @Assert(approxEqual(exponweibPPF(0.2, 2.0, 2.0, loc: 1.0, scale: 2.0),     2.5398488246794986, atol:1e-13))
    }
}