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

    let y = (x - loc) / scale   

    if(y < (-Float64.getPI() * 0.25) || y > Float64.getPI() * 0.25) {
        throw IllegalArgumentException("anglitLogPDF: input value out of bound.")
    }

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

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

    return log(anglitPDF(x, loc: loc, scale: scale))
}

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

    let y = (x - loc) / scale

    if(y < (-Float64.getPI() * 0.25) || y > Float64.getPI() * 0.25) {
        throw IllegalArgumentException("anglitPDF: input value out of bound.")
    }

    return cos(2.0 * y) / scale
}


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

    let y = (x - loc) / scale

    if(y < (-Float64.getPI() * 0.25) || y > Float64.getPI() * 0.25) {
        throw IllegalArgumentException("anglitCDF: input value out of bound.")
    }

    return 0.5 + 0.5 * sin(2.0 * y)
}


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

    let y = (x - loc) / scale

    if(y < (-Float64.getPI() * 0.25) || y > Float64.getPI() * 0.25) {
        throw IllegalArgumentException("anglitLogCDF: input value out of bound.")
    }

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

    return log(temp)
}


/*
 * compute the mean
 */
public func anglitMean(loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    return loc
}


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

/*
 * compute the std
 */
public func anglitStd(loc!: Float64 = 0.0, scale!: Float64 = 1.0): Float64 {
    return sqrt(anglitVar(loc: loc, scale: scale))
}

@Test
public class TestAnglit {
    @TestCase
    func testAnglit(): Unit {
        @Assert(approxEqual(anglitLogPDF(1.0, loc: 2.0, scale: 3.0), -1.3395542196347185, atol:1e-13))
        @Assert(approxEqual(anglitLogCDF(1.0, loc: 2.0, scale: 3.0), -1.656450390681924,  atol:1e-13))
        @Assert(approxEqual(anglitMean(loc: 2.0, scale: 3.0),         2.0,                atol:1e-13))
        @Assert(approxEqual(anglitVar(loc: 2.0, scale: 3.0),          1.0516524756127639, atol:1e-13))
    }
}