package scientific.stats.continuous

import std.math.*
import std.unittest.*
import std.unittest.testmacro.*

import scientific.numbers.*
import scientific.stats.random.*
import scientific.stats.normal.*

foreign func erfc(a: Float64): Float64

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

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

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

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

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

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


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

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

    let res = unsafe { erfc(sqrt(0.5 / y)) }
    return res
}


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

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

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

    return log(temp)
}


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

    let temp = normalPPF(0.5 * q)
    let res = 1.0 / (temp * temp)
    return res * scale + loc
}


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


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

/*
 * compute the std
 */
public func levyStd(loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    return Float64.Inf
}

@Test
public class TestLevy {
    @TestCase
    func testLevy(): Unit {
        @Assert(approxEqual(levyLogPDF(3.0, loc: 2.0, scale: 1.0), -1.4189385332046727, atol:1e-13))
        @Assert(approxEqual(levyLogCDF(3.0, loc: 2.0, scale: 1.0), -1.147874464449318,  atol:1e-13))
        @Assert(approxEqual(levyPPF(0.7, loc: 2.0, scale: 1.0),     8.735282952993835,  atol:1e-7))
    }
}