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

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

    var res = log(0.5 * k)
    if (0.0 < y && y < 1.0) {
        res = res + (k - 1.0) * log(y)
    } else if (y >= 1.0) {
        res = res + (-k - 1.0) * log(y)
    } else {
        throw IllegalArgumentException("loglaplaceLogPDF: input value out of bound.")
    }

    return res - log(scale)
}

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

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

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

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

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

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

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

    var res = 0.0
    if (y < 1.0) {
        res = 0.5 * pow(y, k)
    } else {
        res = 1.0 - 0.5 * pow(y, -k)
    }

    return res
}

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

    if (k <= 0.0) {
        throw IllegalArgumentException("loglaplaceLogCDF: shape parameter out of bound.")
    }
    if (y <= 0.0) {
        throw IllegalArgumentException("loglaplaceLogCDF: input value out of bound.")
    }

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

    return log(temp)
}


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

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

    var res = 0.0
    if (q < 0.5) {
        res = pow(2.0 * q, 1.0 / k)
    } else {
        res = pow(2.0 - 2.0 * q, -1.0/ k)
    }
    return res * scale + loc
}

@Test
public class TestLogLaplace {
    @TestCase
    func testLoglaplace(): Unit {
        @Assert(approxEqual(loglaplaceLogPDF(2.0, 2.0, loc: 1.0, scale: 2.0), -1.3862943611198906, atol:1e-13))
        @Assert(approxEqual(loglaplaceLogCDF(2.0, 2.0, loc: 1.0, scale: 2.0), -2.0794415416798357, atol:1e-13))
        @Assert(approxEqual(loglaplacePPF(0.2, 2.0, loc: 1.0, scale: 2.0),     2.264911064067352,  atol:1e-13))
    }
}