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

    let y = (x - loc) / scale

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

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

    return (3.0 * log(k) - 0.5 * log(2.0 * Float64.getPI()) -
            log(helpSF(k)) + log(y) + 0.5 * log(1.0 - y * y) + (- 0.5 * k * k * (1.0 - y * y)) - log(scale))
}



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

    let y = (x - loc) / scale

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

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

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


/*
 * Survival function
 */
public func argusSF(x: Float64, k: Float64): Float64 {

    if (x <= 0.0 || x >= 1.0) {
        throw IllegalArgumentException("argusLogPDF: input value x out of bound.")
    }

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

    let temp = k * sqrt(1.0 - x * x)
    return helpSF(temp) / helpSF(k)
}


public func helpSF(k: Float64): Float64 {
    return normalCDF(k) - k * normalPDF(k) - 0.5
}


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

    let y = (x - loc) / scale

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

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

    return 1.0 - argusSF(y, k)
}


/*
 * Log of probability density function
 */
public func argusLogCDF(x: Float64, k: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {

    let y = (x - loc) / scale

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

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

    let temp = argusCDF(x, k, loc: loc, scale: scale)

    if (temp < 0.000001) {
        throw IllegalArgumentException("argusLogCDF: return-value too small.")
    }

    return log(temp)
}


@Test
public class TestArgus {
    @TestCase
    func testArgus(): Unit {
        @Assert(approxEqual(argusPDF(0.2, 1.0),                           0.4867897904846194,  atol:1e-13))
        @Assert(approxEqual(argusCDF(0.2, 1.0),                           0.04869236463755611, atol:1e-13))
        @Assert(approxEqual(argusLogCDF(0.2, 1.0),                       -3.022233044808508,   atol:1e-13))
        @Assert(approxEqual(argusLogPDF(1.2, 1.0, loc: 1.0, scale: 2.0), -2.105831422585523,   atol:1e-13))
    }
}