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 mielkeLogPDF(x: Float64, a: Float64, c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale
    if (a <= 0.0 || c <= 0.0) {
        throw IllegalArgumentException("mielkeLogCDF: shape parameter out of bound.")
    }
    if (y <= 0.0) {
        throw IllegalArgumentException("mielkeLogCDF: input value out of bound.")
    }

    let res = log(a) + (a - 1.0) * log(y) - (1.0 + a / c) * log(1.0 + pow(y, c))
    return res - log(scale)
}

/*
 * Probability density function
 */
public func mielkePDF(x: Float64, a: Float64, c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale
    if (a <= 0.0 || c <= 0.0) {
        throw IllegalArgumentException("mielkePDF: shape parameter out of bound.")
    }
    if (y <= 0.0) {
        throw IllegalArgumentException("mielkePDF: input value out of bound.")
    }

    let res = exp(mielkeLogPDF(x, a, c, loc: loc, scale: scale))
    return res
}

/*
 * Cumulative probability density function
 */
public func mielkeCDF(x: Float64, a: Float64, c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale
    if (a <= 0.0 || c <= 0.0) {
        throw IllegalArgumentException("mielkeLogCDF: shape parameter out of bound.")
    }
    if (y <= 0.0) {
        throw IllegalArgumentException("mielkeLogCDF: input value out of bound.")
    }

    let res = exp(mielkeLogCDF(x, a, c, loc: loc, scale: scale))
    return res
}

/*
 * Log of Cumulative probability density function
 */
public func mielkeLogCDF(x: Float64, a: Float64, c: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale
    if (a <= 0.0 || c <= 0.0) {
        throw IllegalArgumentException("mielkeCDF: shape parameter out of bound.")
    }
    if (y <= 0.0) {
        throw IllegalArgumentException("mielkeCDF: input value out of bound.")
    }

    let res = a * log(y) - (a * 1.0 / c) * log(1.0 + pow(y, c))
    return res 
}

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

    let temp = pow(q, c * 1.0 / a)
    let res = pow(temp / (1.0 - temp), 1.0 / c)
    return res * scale + loc
}

@Test
public class TestMielke {
    @TestCase
    func testMielke(): Unit {
        @Assert(approxEqual(mielkeLogPDF(2.0, 10.0, 4.0, loc: 1.0, scale: 2.0), -4.8410728889629295, atol:1e-13))
        @Assert(approxEqual(mielkeLogCDF(2.0, 10.0, 4.0, loc: 1.0, scale: 2.0), -7.08303336014054,   atol:1e-13))
        @Assert(approxEqual(mielkePPF(0.2, 10.0, 4.0, loc: 1.0, scale: 2.0),     3.0513013186026394, atol:1e-6))
    }
}