package scientific.stats.continuous

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

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


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

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

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

    let temp = -(y - 1.0) * (y - 1.0) / (2.0 * y * k * k)
    return log(1.0 + y) - log(2.0 * k) - 0.5 * log(2.0 * Float64.getPI()) - 1.5 * log(y) + temp - log(scale)
}

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

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

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

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

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

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

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

    let temp = 1.0 / k * (sqrt(y) - 1.0 / sqrt(y))

    return normalCDF(temp)
}

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

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

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

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

    return log(temp) 
}


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

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

    let temp = k * normalPPF(q)
    let t = temp + sqrt(temp * temp + 4.0)

    return (0.25 * t * t) * scale + loc
}

/*
 * Mean of the distribution.
 */
public func fatiguelifeMean(k: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    if (k <= 0.0) {
        throw IllegalArgumentException("fatiguelifeMean: shape parameter out of bound.")
    }

    return (0.5 * k * k + 1.0) * scale + loc
}  


/*
 * Var of the distribution.
 */
public func fatiguelifeVar(k: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    if (k <= 0.0) {
        throw IllegalArgumentException("fatiguelifeVar: shape parameter out of bound.")
    }

    let temp = 5.0 * k * k + 4.0

    return temp * k * k / 4.0 * scale * scale
}  


/*
 * Std of the distribution.
 */
public func fatiguelifeStd(k: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    if (k <= 0.0) {
        throw IllegalArgumentException("fatiguelifeStd: shape parameter out of bound.")
    }

    return sqrt(fatiguelifeVar(k, loc: loc, scale: scale))
}  

@Test
public class TestFatigueLife {
    @TestCase
    func testFatiguelife(): Unit {
        @Assert(approxEqual(fatiguelifeLogPDF(2.0, 2.0, loc: 1.0, scale: 2.0), -1.6156941959364262, atol:1e-13))
        @Assert(approxEqual(fatiguelifeLogCDF(2.0, 2.0, loc: 1.0, scale: 2.0), -1.0165619839535647, atol:1e-13)) 
        @Assert(approxEqual(fatiguelifePPF(0.2, 2.0, loc: 1.0, scale: 2.0),     1.4332098420421369, atol:1e-9)) 
        @Assert(approxEqual(fatiguelifeMean(2.0, loc: 1.0, scale: 2.0),         7.0,                atol:1e-13)) 
        @Assert(approxEqual(fatiguelifeVar(2.0, loc: 1.0, scale: 2.0),          96.0,               atol:1e-13)) 
    }
}