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

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

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

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

    return log(temp)
}

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

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

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

    let t1 = 1.0 / (Float64.getPI() * (1.0 + (y - k) * (y - k)))
    let t2 = 1.0 / (Float64.getPI() * (1.0 + (y + k) * (y + k)))
    return (t1 + t2) / scale
}


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

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

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

    return log(temp)
}

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

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

    return 1.0 / Float64.getPI() * (atan(y - k) + atan(y + k))
}

@Test
public class TestFoldCauchy {
    @TestCase
    func testFoldcauchy(): Unit {
        @Assert(approxEqual(foldcauchyLogPDF(2.0, 2.0, loc: 1.0, scale: 2.0), -2.6461582744540975, atol:1e-13))
        @Assert(approxEqual(foldcauchyLogCDF(2.0, 2.0, loc: 1.0, scale: 2.0), -2.717372011180594,  atol:1e-13))
    }
}