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

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

    let res = 0.5 * log(2.0) - 0.5 * log(Float64.getPI()) - 0.5 * y * y 
    return res - log(scale)
}

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

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

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


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

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

    return normalCDF(y) * 2.0 - 1.0
}


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

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

    let temp = halfnormCDF(x, loc: loc, scale: scale)

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

    return log(temp)
}


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

    let res = normalPPF(0.5 * (1.0 + q))
    return res * scale + loc
}


/*
 * Sample
 */
public func halfnormSample(r: Random, loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    let res = abs(normalSampleFloat64(r))
    return res * scale + loc
}


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


/*
 * compute the var
 */
public func halfnormVar(loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    return  (1.0 - 2.0 / Float64.getPI()) * scale * scale
}

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

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

    return sqrt(temp)
}

@Test
public class TestHalfnorm {
    @TestCase
    func testHalfnorm(): Unit {
        @Assert(approxEqual(halfnormLogPDF(3.0, loc: 2.0, scale: 1.0), -0.7257913526447274,  atol:1e-13))
        @Assert(approxEqual(halfnormLogCDF(3.0, loc: 2.0, scale: 1.0), -0.38171514630212616, atol:1e-13))
        @Assert(approxEqual(halfnormPPF(0.7, loc: 2.0, scale: 1.0),     3.03643338949379,    atol:1e-8))
        @Assert(approxEqual(halfnormMean(loc: 2.0, scale: 1.0),         2.7978845608028653,  atol:1e-13))
        @Assert(approxEqual(halfnormVar(loc: 2.0, scale: 1.0),          0.3633802276324186,  atol:1e-13))
    }
}