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

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

    let res = log(2.0) - log(Float64.getPI()) + 0.5 * log(1.0 - y * y)
    return res - log(scale)
}

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

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

    let temp = semicircularLogPDF(x, loc: loc, scale: scale)
    return exp(temp)
}


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

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

    let res = 0.5 + 1.0 / Float64.getPI() * (y * sqrt(1.0 - y * y) + asin(y))
    return res
}


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

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

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

    return log(temp)
}


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

    let res = rdistPPF(q, 3.0)
    return res * scale + loc
}


/*
 * compute the mean
 */
public func semicircularMean(loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    return loc
}


/*
 * compute the var
 */
public func semicircularVar(loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    return 0.25 * scale * scale
}

/*
 * compute the std
 */
public func semicircularStd(loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let temp = halfnormVar(loc: loc, scale: scale)

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

    return sqrt(temp)
}


@Test
public class TestSemicircular {
    @TestCase
    func testSemicircular(): Unit {
        @Assert(approxEqual(semicircularLogPDF(2.0, loc: 2.0, scale: 1.0), -0.4515827052894548, atol:1e-13))
        @Assert(approxEqual(semicircularLogCDF(2.0, loc: 2.0, scale: 1.0), -0.6931471805599453, atol:1e-13))
        @Assert(approxEqual(semicircularPPF(0.7, loc: 2.0, scale: 1.0),     2.3196915097905038, atol:1e-6))
    }
}