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 triangLogPDF(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("triangLogPDF: shape parameter out of bound.")
    }

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

    return log(temp)
}

/*
 * Probability density function
 */
public func triangPDF(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("triangPDF: shape parameter out of bound.")
    }

    var res = 0.0
    if (k == 0.0) {
        res = 2.0 - 2.0 * y
    } else if (y < k) {
        res = 2.0 * y / k
    } else if (y >= k && k != 1.0) {
        res = 2.0 * (1.0 - y) / (1.0 - k)
    } else {
        res = 2.0 * y
    }
    return res / scale
}

/*
 * Cumulative probability density function
 */
public func triangCDF(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("triangCDF: shape parameter out of bound.")
    }

    var res = 0.0
    if (k == 0.0) {
        res = 2.0 * y - y * y
    } else if (y < k) {
        res = y * y / k
    } else if (y >= k && k != 1.0) {
        res = (y * y - 2.0 * y + k) / (k - 1.0)
    } else {
        res = y * y
    }
    return res
}

/*
 * Log of Cumulative probability density function
 */
public func triangLogCDF(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("triangLogCDF: shape parameter out of bound.")
    }

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

    return log(temp)
}

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

    var res = 0.0
    if (q < k) {
        res = sqrt(k * q)
    } else {
        res = 1.0 - sqrt((1.0 - k) * (1.0 - q))
    }

    return res * scale + loc
}

@Test
public class TestTriang {
    @TestCase
    func testTriang(): Unit {
        @Assert(approxEqual(triangLogPDF(2.0, 0.2, loc: 1.0, scale: 2.0), -0.4700036292457355, atol:1e-13))
        @Assert(approxEqual(triangLogCDF(2.0, 0.2, loc: 1.0, scale: 2.0), -0.3746934494414107, atol:1e-13))
        @Assert(approxEqual(triangPPF(0.2, 0.2, loc: 1.0, scale: 2.0),     1.4,                atol:1e-13))  
    }
}