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 wrapcauchyLogPDF(x: Float64, k: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale
    if (k <= 0.0 || k >= 1.0) {
        throw IllegalArgumentException("wrapcauchyLogPDF: shape parameter out of bound.")
    }
    if (y < 0.0 || y > 2.0 * Float64.getPI()) {
        throw IllegalArgumentException("wrapcauchyLogPDF: input value out of bound.")
    }

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

    return log(temp)
}

/*
 * Probability density function
 */
public func wrapcauchyPDF(x: Float64, k: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale
    if (k <= 0.0 || k >= 1.0) {
        throw IllegalArgumentException("wrapcauchyPDF: shape parameter out of bound.")
    }
    if (y < 0.0 || y > 2.0 * Float64.getPI()) {
        throw IllegalArgumentException("wrapcauchyLogPDF: input value out of bound.")
    }

    let res = (1.0 - k * k) / (2.0 * Float64.getPI() * (1.0 + k * k - 2.0 * k * cos(y)))
    return res / scale
}

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

    if (k <= 0.0 || k >= 1.0) {
        throw IllegalArgumentException("wrapcauchyCDF: shape parameter out of bound.")
    }
    if (y < 0.0 || y > 2.0 * Float64.getPI()) {
        throw IllegalArgumentException("wrapcauchyLogPDF: input value out of bound.")
    }

    var res = 0.0
    let t = (1.0 + k) / (1.0 - k)
    if (y < Float64.getPI()) {
        res = 1.0 / Float64.getPI() * atan(t * tan(0.5 * y))
    } else {
        res = 1.0 - 1.0 / Float64.getPI() * atan(t * tan((2.0 * Float64.getPI() - y) / 2.0))
    }
    return res
}

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

    if (k <= 0.0 || k >= 1.0) {
        throw IllegalArgumentException("wrapcauchyLogCDF: shape parameter out of bound.")
    }
    if (y < 0.0 || y > 2.0 * Float64.getPI()) {
        throw IllegalArgumentException("wrapcauchyLogPDF: input value out of bound.")
    }

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

    return log(temp)
}

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

    var res = 0.0
    let t = (1.0 - k) / (1.0 + k)
    if (q < 0.5) {
        res = 2.0 * atan(t * tan(Float64.getPI() * q))
    } else {
        res = 2.0 * Float64.getPI() - 2.0 * atan(t * tan(Float64.getPI() * (1.0 - q)))
    }

    return res * scale + loc
}

@Test
public class TestWrapCauchy {
    @TestCase
    func testWrapcauchy(): Unit {
        @Assert(approxEqual(wrapcauchyLogPDF(2.0, 0.2, loc: 1.0, scale: 2.0), -2.1992843009410437, atol:1e-13))
        @Assert(approxEqual(wrapcauchyLogCDF(2.0, 0.2, loc: 1.0, scale: 2.0), -2.1504610604374745, atol:1e-13))
        @Assert(approxEqual(wrapcauchyPPF(0.2, 0.2, loc: 1.0, scale: 2.0),     2.804235486010719,  atol:1e-13))
    }
}