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 truncweibull_minLogPDF(x: Float64, a: Float64, b: Float64, c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale
    var paramCheck = true
    if (0.0 <= a && 0.0 < c && a < b) {
        paramCheck = false
    }
    if (paramCheck) {
        throw IllegalArgumentException("truncweibull_minLogPDF: shape parameter out of bound.")
    }

    if (y <= a || y > b) {
        throw IllegalArgumentException("truncweibull_minLogPDF: input value out of bound.")
    }

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

    return log(temp)
}

/*
 * Probability density function
 */
public func truncweibull_minPDF(x: Float64, a: Float64, b: Float64, c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale
    var paramCheck = true
    if (0.0 <= a && 0.0 < c && a < b) {
        paramCheck = false
    }
    if (paramCheck) {
        throw IllegalArgumentException("truncweibull_minPDF: shape parameter out of bound.")
    }
    if (y <= a || y > b) {
        throw IllegalArgumentException("truncweibull_minPDF: input value out of bound.")
    }

    let t = exp(-pow(a, c)) - exp(-pow(b, c))
    let res = c * pow(y, c - 1.0) * exp(-pow(y, c)) / t
    return res / scale
}

/*
 * Cumulative probability density function
 */
public func truncweibull_minCDF(x: Float64, a: Float64, b: Float64, c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale
    var paramCheck = true
    if (0.0 <= a && 0.0 < c && a < b) {
        paramCheck = false
    }
    if (paramCheck) {
        throw IllegalArgumentException("truncweibull_minLogCDF: shape parameter out of bound.")
    }
    if (y <= a || y > b) {
        throw IllegalArgumentException("truncweibull_minLogCDF: input value out of bound.")
    }

    let t1 = exp(-pow(a, c)) - exp(-pow(y, c))
    let t2 = exp(-pow(a, c)) - exp(-pow(b, c))
    return t1 / t2
}

/*
 * Log of Cumulative probability density function
 */
public func truncweibull_minLogCDF(x: Float64, a: Float64, b: Float64, c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale
    var paramCheck = true
    if (0.0 <= a && 0.0 < c && a < b) {
        paramCheck = false
    }
    if (paramCheck) {
        throw IllegalArgumentException("truncweibull_minLogCDF: shape parameter out of bound.")
    }
    if (y <= a || y > b) {
        throw IllegalArgumentException("truncweibull_minLogCDF: input value out of bound.")
    }

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

    return log(temp)
}

/*
 * PPF
 */
public func truncweibull_minPPF(q: Float64, a: Float64, b: Float64, c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    var paramCheck = true
    if (0.0 <= a && 0.0 < c && a < b) {
        paramCheck = false
    }
    if (paramCheck) {
        throw IllegalArgumentException("truncweibull_minPPF: shape parameter out of bound.")
    }
    
    if (q <= 0.0 || q >= 1.0) {
        throw IllegalArgumentException("truncweibull_minPPF: quantile out of bound.")
    }

    let t = -log((1.0 - q) * exp(-pow(a, c)) + q * exp(-pow(b, c)))
    let res = pow(t, 1.0 / c)
    return res * scale + loc
}

@Test
public class TestTruncWeibull {
    @TestCase
    func testTruncweibull_min(): Unit {
        @Assert(approxEqual(truncweibull_minLogPDF(4.0, 1.0, 2.0, 1.0, loc: 1.0, scale: 2.0), -0.7344720351728634,  atol:1e-13))
        @Assert(approxEqual(truncweibull_minLogCDF(4.0, 1.0, 2.0, 1.0, loc: 1.0, scale: 2.0), -0.47407698418010646, atol:1e-13))
        @Assert(approxEqual(truncweibull_minPPF(0.2, 1.0, 2.0, 1.0, loc: 1.0, scale: 2.0),     3.270320549673619,   atol:1e-6))
    }
}