package scientific.stats.continuous

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

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

/*
 * Log of Probability density function
 */
public func johnsonsuLogPDF(x: Float64, a: Float64, c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale

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

    let temp = a + c * log(y + sqrt(1.0 + y * y))
    return log(c) - 0.5 * log(1.0 + y * y) + log(normalPDF(temp)) - log(scale)
}

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

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

    return exp(johnsonsuLogPDF(x, a, c, loc: loc, scale: scale))
}

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

    if (c <= 0.0) {
        throw IllegalArgumentException("johnsonsuCDF: shape parameter out of bound.")
    }
    
    let temp = a + c * log(y + sqrt(1.0 + y * y))
    return normalCDF(temp)
}

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

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

    let temp = johnsonsuCDF(x, a, c, loc: loc, scale: scale)

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

    return log(temp)
}

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

    let r1 = (normalPPF(q) - a) / c
    let res = 0.5 * (exp(r1) - exp(-r1))
    return res * scale + loc
}

@Test
public class TestJohnsonsu {
    @TestCase
    func testJohnsonsu(): Unit {
        @Assert(approxEqual(johnsonsuLogPDF(2.0, 2.0, 3.0, loc: 1.0, scale: 2.0), -6.554357843708609,     atol:1e-13))
        @Assert(approxEqual(johnsonsuLogCDF(2.0, 2.0, 3.0, loc: 1.0, scale: 2.0), -0.0002870156673103043, atol:1e-13)) 
        @Assert(approxEqual(johnsonsuPPF(0.2, 2.0, 3.0, loc: 1.0, scale: 2.0),    -1.1906753781779535,    atol:1e-6))
    }
}