package scientific.stats.discrete

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

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

/*
 * Log of Probability density function
 */
public func boltzmannLogPDF(x: Int64, n: Int64, a: Float64): Float64 {
    if (n <= 0 || x < 0 || x >= n || a <= 0.0) {
        throw IllegalArgumentException("boltzmannLogPDF: parameters out of bound.")
    }

    let res = boltzmannPDF(x, n, a)
    if (res < 0.000001) {
        throw IllegalArgumentException("boltzmannLogPDF: return-value too small.")
    }

    return log(res)
}


/*
 * Probability density function
 */
public func boltzmannPDF(x: Int64, n: Int64, a: Float64): Float64 {
    if (n <= 0 || x < 0 || x >= n || a <= 0.0) {
        throw IllegalArgumentException("boltzmannPDF: parameters out of bound.")
    }

    let x_f = Float64(x)
    let n_f = Float64(n)
    let t1 = 1.0 - exp(-a)
    let t2 = exp(-a * x_f)
    let t3 = 1.0 - exp(-a * n_f)
    return t1 * t2 / t3
}


/*
 * Cumulative probability density function
 */
public func boltzmannCDF(x: Int64, n: Int64, a: Float64): Float64 {
    if (n <= 0 || x < 0 || x >= n || a <= 0.0) {
        throw IllegalArgumentException("boltzmannCDF: parameters out of bound.")
    }

    let x_f = Float64(x)
    let n_f = Float64(n)

    let t1 = 1.0 - exp(-a * (x_f + 1.0))
    let t2 = 1.0 - exp(-a * n_f)
    return t1 / t2
}


/*
 * Log of cumulative probability density function
 */
public func boltzmannLogCDF(x: Int64, n: Int64, a: Float64): Float64 {
    if (n <= 0 || x < 0 || x >= n || a <= 0.0) {
        throw IllegalArgumentException("boltzmannCDF: parameters out of bound.")
    }

    let res = boltzmannCDF(x, n, a)
    if (res < 0.000001) {
        throw IllegalArgumentException("boltzmannCDF: return-value too small.")
    }

    return log(res)
}


@Test
public class TestBoltzmann {
    @TestCase
    func testBoltzmann(): Unit {
        @Assert(approxEqual(boltzmannPDF(2, 5, 0.2), 0.1922234742163609, atol:1e-13))
        @Assert(approxEqual(boltzmannCDF(2, 5, 0.2), 0.7137694821097315, atol:1e-8))
    }
}