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 bradfordLogPDF(x: Float64, c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {

    let y = (x - loc) / scale

    if(y < 0.0 || y > 1.0) {
        throw IllegalArgumentException("bradfordLogPDF: input value x out of bound.")
    }

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

    return log(c) - log(log(1.0 + c)) - log(1.0 + c * y) - log(scale)
}

/*
 * Probability density function
 */
public func bradfordPDF(x: Float64, c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {

    let y = (x - loc) / scale

    if(y < 0.0 || y > 1.0) {
        throw IllegalArgumentException("bradfordLogPDF: input value x out of bound.")
    }

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

    return exp(bradfordLogPDF(x, c, loc: loc, scale: scale))
}


/*
 * Cumulative probability density function
 */
public func bradfordCDF(x: Float64, c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {

    let y = (x - loc) / scale

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

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

    return  log(1.0 + c * y) / log(1.0 + c) 
}


/*
 * Log of cumulative probability density function
 */
public func bradfordLogCDF(x: Float64, c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {

    let y = (x - loc) / scale

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

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

    let temp = bradfordCDF(x, c, loc: loc, scale: scale)

    if(temp < 0.000001) {
        throw IllegalArgumentException("bradfordCDF: return-value too small.")
    }

    return  log(temp)
}


/*
 * ppf
 */
public func bradfordPPF(q: Float64, c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {

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

    let temp = (exp(q * log(1.0 + c)) - 1.0) / c

    return  temp * scale + loc
}


/*
 * mean
 */
public func bradfordMean(c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {

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

    let temp = log(1.0 + c)
    
    return  (c - temp) / (c * temp) * scale + loc
}


/*
 * var
 */
public func bradfordVar(c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {

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

    let temp = log(1.0 + c)
    
    return ((c + 2.0) * temp - 2.0 * c) / (2.0 * c * temp * temp) * scale * scale
}


@Test
public class TestBradford {
    @TestCase
    func testBradford(): Unit {
        @Assert(approxEqual(bradfordLogPDF(1.2, 2.0, loc: 1.0, scale: 2.0), -0.2763693844106537,  atol:1e-13))
        @Assert(approxEqual(bradfordLogCDF(1.2, 2.0, loc: 1.0, scale: 2.0), -1.7960311828981994,  atol:1e-13))
        @Assert(approxEqual(bradfordPPF(0.2, 2.0, loc: 1.0, scale: 2.0),     1.2457309396155174,  atol:1e-13)) 
        @Assert(approxEqual(bradfordMean(2.0, loc: 1.0, scale: 2.0),         1.8204784532536746,  atol:1e-13)) 
        @Assert(approxEqual(bradfordVar(2.0, loc: 1.0, scale: 2.0),          0.32681510774645767, atol:1e-13)) 
    }
}