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

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

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

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

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

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

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

    return exp(weibull_minLogPDF(x, c, loc: loc, scale: scale))
}

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

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

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

    let temp = - pow(y, c)
    return 1.0 - exp(temp)
}

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

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

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

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


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

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

    let temp = - log(1.0 - q)
    return pow(temp, 1.0 / c) * scale + loc
}


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

    let r1 = r.nextFloat64()
    let temp = log(-log(1.0 - r1)) / c

    return pow(Float64.getE(), temp) * scale + loc
}

@Test
public class TestWeibullMin {
    @TestCase
    func testWeibull_min(): Unit {
        @Assert(approxEqual(weibull_minLogPDF(2.0, 2.0, loc: 1.0, scale: 2.0), -0.9431471805599453, atol:1e-13))
        @Assert(approxEqual(weibull_minLogCDF(2.0, 2.0, loc: 1.0, scale: 2.0), -1.508691549446032,  atol:1e-13))
        @Assert(approxEqual(weibull_minPPF(0.2, 2.0, loc: 1.0, scale: 2.0),     1.9447614541548777, atol:1e-13))
    }
}