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

    if (k <= 0.0) {
        throw IllegalArgumentException("rdistLogPDF: shape parameter out of bound.")
    }
    if (y < -1.0 || y > 1.0) {
        throw IllegalArgumentException("rdistLogPDF: input value out of bound.")
    }

    let res = (0.5 * k - 1.0) * log(1.0 - y * y) - betaLog(0.5, 0.5 * k)
    return res - log(scale)
}

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

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

    if (y < -1.0 || y > 1.0) {
        throw IllegalArgumentException("rdistPDF: input value out of bound.")
    }

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

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

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

    if (y < -1.0 || y > 1.0) {
        throw IllegalArgumentException("rdistPDF: input value out of bound.")
    }

    let res = betaCDF(0.5 * (y + 1.0), 0.5 * k, 0.5 * k)
    return res
}

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

    if (k <= 0.0) {
        throw IllegalArgumentException("rdistLogCDF: shape parameter out of bound.")
    }
    if (y < -1.0 || y > 1.0) {
        throw IllegalArgumentException("rdistLogCDF: input value out of bound.")
    }

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

    return log(temp)
}


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

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

    let res = 2.0 * betaPPF(q, 0.5 * k, 0.5 * k) - 1.0
    return res * scale + loc
}

@Test
public class TestRdist {
    @TestCase
    func testRdist(): Unit {
        @Assert(approxEqual(rdistLogPDF(2.0, 2.0, loc: 1.0, scale: 2.0), -1.3862943611198906, atol:1e-13))
        @Assert(approxEqual(rdistLogCDF(2.0, 2.0, loc: 1.0, scale: 2.0), -0.2876820724517809, atol:1e-13))
        @Assert(approxEqual(rdistPPF(0.2, 2.0, loc: 1.0, scale: 2.0),    -0.2,                atol:1e-6))
    }
}