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 powerlognormLogPDF(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("powerlognormLogCDF: shape parameter out of bound.")
    }
    if (y <= 0.0) {
        throw IllegalArgumentException("powerlognormLogCDF: input value out of bound.")
    }

    let t = log(y) / c

    let res = log(a) - log(y) - log(c) + log(normalPDF(t)) + (a - 1.0) * log(normalCDF(-t))
    return res - log(scale)
}

/*
 * Probability density function
 */
public func powerlognormPDF(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("powerlognormPDF: shape parameter out of bound.")
    }
    if (y <= 0.0) {
        throw IllegalArgumentException("powerlognormPDF: input value out of bound.")
    }

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

/*
 * Cumulative probability density function
 */
public func powerlognormCDF(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("powerlognormLogCDF: shape parameter out of bound.")
    }
    if (y <= 0.0) {
        throw IllegalArgumentException("powerlognormLogCDF: input value out of bound.")
    }

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

/*
 * Log of Cumulative probability density function
 */
public func powerlognormLogCDF(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("powerlognormCDF: shape parameter out of bound.")
    }
    if (y <= 0.0) {
        throw IllegalArgumentException("powerlognormCDF: input value out of bound.")
    }

    let temp = powerlognormCDF(x, a, c, loc:loc, scale: scale)
    if (temp < 0.000001) {
        throw IllegalArgumentException("powerlognormCDF: return-value too small.")
    }

    return log(temp)
}

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

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

@Test
public class TestPowerLognorm {
    @TestCase
    func testPowerlognorm(): Unit {
        @Assert(approxEqual(powerlognormLogPDF(2.0, 2.0, 1.0, loc: 1.0, scale: 2.0), -0.7458754179434196, atol:1e-13))
        @Assert(approxEqual(powerlognormLogCDF(2.0, 2.0, 1.0, loc: 1.0, scale: 2.0), -0.8471654374077482, atol:1e-13))
        @Assert(approxEqual(powerlognormPPF(0.2, 2.0, 1.0, loc: 1.0, scale: 2.0),     1.5727681262250766, atol:1e-6))
    }
}