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 rayleighLogPDF(x: Float64, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let y = (x - loc) / scale

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

    let res = log(y) - 0.5 * y * y
    return res - log(scale)
}

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

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

    let temp = rayleighLogPDF(x, loc: loc, scale: scale)
    return exp(temp)
}


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

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

    let res = 1.0 - exp(-0.5 * y * y)
    return res
}


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

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

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

    return log(temp)
}


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

    let res = sqrt(-2.0 * log(1.0 - q))
    return res * scale + loc
}


/*
 * compute the mean
 */
public func rayleighMean(loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let res = sqrt(0.5 * Float64.getPI())
    return res * scale + loc
}


/*
 * compute the var
 */
public func rayleighVar(loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let res = 0.5 * (4.0 - Float64.getPI())
    return res * scale * scale
}

/*
 * compute the std
 */
public func rayleighStd(loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let temp = halfnormVar(loc: loc, scale: scale)

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

    return sqrt(temp)
}


/*
 * sample
 */
public func rayleighSample(r: Random, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let res = chiSample(r, 2.0)
    return res * scale + loc
}

@Test
public class TestRayleigh {
    @TestCase
    func testRayleigh(): Unit {
        @Assert(approxEqual(rayleighLogPDF(3.0, loc: 2.0, scale: 1.0), -0.5,                atol:1e-13))
        @Assert(approxEqual(rayleighLogCDF(3.0, loc: 2.0, scale: 1.0), -0.9327521295671886, atol:1e-13))
        @Assert(approxEqual(rayleighPPF(0.7, loc: 2.0, scale: 1.0),     3.551755653655521,  atol:1e-13))
    }
}