package scientific.stats.continuous

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

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

foreign func fdtr(a: Int32, b: Int32, t: Float64): Float64

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

    if (dfn <= 0 || dfd <= 0) {
        throw IllegalArgumentException("fLogPDF: freedom parameter should be positive.")
    }

    if (y <= 0.0) {
        throw IllegalArgumentException("fLogPDF: input value x out of bound.")
    }
    
    let d1 = Float64(dfn)
    let d2 = Float64(dfd)

    return  0.5 * d2 * log(d2) + 0.5 * d1 * log(d1) + (0.5 * d1 - 1.0) * log(y) -
            0.5 * (d1 + d2) * log(d2 + d1 * y) - betaLog(0.5 * d1, 0.5 * d2) - log(scale)
} 


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

    if (dfn <= 0 || dfd <= 0) {
        throw IllegalArgumentException("fPDF: freedom parameter should be positive.")
    }

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

    return exp(fLogPDF(x, dfn, dfd, loc: loc, scale: scale))
}  


/*
/*
 * Cumulative probability density function
 */
public func fCDF(x: Float64, dfn: Int64, dfd: Int64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {

    let y = (x - loc) / scale

    if(dfn <= 0 || dfd <= 0) {
        throw IllegalArgumentException("fCDF: freedom parameter should be positive.")
    }

    if(y <= 0.0) {
        throw IllegalArgumentException("fCDF: input value x out of bound.")
    }
    
    return unsafe { fdtr(Int32(dfn), Int32(dfd), y) }
}

/*
 * Log of Cumulative probability density function
 */
public func fLogCDF(x: Float64, dfn: Int64, dfd: Int64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {

    let y = (x - loc) / scale

    if(dfn <= 0 || dfd <= 0) {
        throw IllegalArgumentException("fLogCDF: freedom parameter should be positive.")
    }

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

    let temp = fCDF(x, dfn, dfd, loc: loc, scale: scale)
    print("${temp}\n")
    if(temp < 0.000001) {
        throw IllegalArgumentException("fLogCDF: return-value too small.")
    }

    return log(temp)
}
*/

/*
 * Mean of the distribution.
 */
public func fMean(dfn: Int64, dfd: Int64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    if (dfd <= 2 || scale == 0.0) {
        return Float64.Inf    //return infinite
    }

    let df2 = Float64(dfd)

    return df2 / (df2 - 2.0) * scale + loc
}  


/*
 * Var of the distribution.
 */
public func fVar(dfn: Int64, dfd: Int64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    if (dfd <= 4) {
        return Float64.Inf    //return infinite
    }

    let df1 = Float64(dfn)
    let df2 = Float64(dfd)

    return 2.0 * df2 * df2 * (df1 + df2 - 2.0) / (df1 * (df2 - 2.0) * (df2 - 2.0) * (df2 - 4.0)) * scale * scale
}  

@Test
public class TestF {
    @TestCase
    func testF(): Unit {
        @Assert(approxEqual(fLogPDF(2.0, 2, 3, loc: 1.0, scale: 2.0), -1.4123523616893974,   atol:1e-13))
        //@Assert(approxEqual(fLogCDF(2.0, 2, 3, loc: 1.0, scale: 2.0), -1.0484489330101874,   atol:1e-13))
        @Assert(approxEqual(fMean(2, 3, loc: 1.0, scale: 2.0),        7.0,                   atol:1e-13))
        @Assert(approxEqual(fVar(2, 6, loc: 1.0, scale: 2.0),         27.0,                  atol:1e-13))
    }
}