using Random, Distributions
################33
# normal
##############3
function pdf_normal(x::Real)
    return pdf(Normal(zero.(x), one(x)), x)
end

function lpdf_normal(x::T) where T<:Real
    return logpdf(Normal(), x)
end
function lpdf_normal(x::AbstractVector{T}) where T<:Real
    return logpdf(MvNormal(zero.(x), one.(x)), x)
end

function ∇lpdf_normal(x)
    return -x 
end

function cdf_normal(x::Real) 
    z = x / sqrt(2.0)
    if abs(x) > 0.0
        return 1.0- 0.5*erfc(z)
    else
        return 0.5 + 0.5*erf(z)
    end
    # return cdf(Normal(zero(x),one(x)), x)
end 

function invcdf_normal(x)
    # if x > 0.5
    #     z = 2.0 - 2.0*x 
    #     return sqrt(2.0)*erfcinv(z)
    # else 
    #     z = 2.0*x - 1.0
    #     return sqrt(2.0)*erfinv(z)
    # end
    return quantile(Normal(zero(x), one(x)), x)
end 

function invcdf_normal(x::BigFloat)
    if x > 0.5
        z = 2.0 - 2.0*x 
        return sqrt(2.0)*erfcinv_finite(z)
    else 
        z = 2.0*x - 1.0
        return sqrt(2.0)*erfinv_finite(z)
    end
end 


############# reasonable erfinv and erfcinv ###########
function erfcinv_finite(y::BigFloat)
    yfloat = Float64(y)
    xfloat = erfcinv(yfloat)
    if isfinite(xfloat)
        x = BigFloat(xfloat)
    else
        # Float64 overflowed, use asymptotic estimate instead
        # from erfc(x) ≈ exp(-x²)/x√π ≈ y  ⟹  -log(yπ) ≈ x² + log(x) ≈ x²
        if yfloat < 1
            x = sqrt(-log(y*sqrtπ))
        else # y must be close to 2
            x = -sqrt(-log((2-y)*sqrtπ))
        end
        # TODO: Newton convergence is slow near y=0 singularity; accelerate?
        isfinite(x) || return x
    end
    sqrtπhalf = sqrtπ * BigFloat(0.5)
    tol = 2eps(abs(x))
    n = 1
    while n < 100000 # Newton iteration for maximal 100000 steps
        Δx = sqrtπhalf * (erfc(x) - y) * exp(x^2)
        x += Δx
        abs(Δx) < tol && break
        n += 1
    end
    return x
end

function erfinv_finite(y::BigFloat)
    xfloat = erfinv(Float64(y))
    if isfinite(xfloat)
        x = BigFloat(xfloat)
    else
        # Float64 overflowed, use asymptotic estimate instead
        # from erfc(x) ≈ exp(-x²)/x√π ≈ y  ⟹  -log(yπ) ≈ x² + log(x) ≈ x²
        x = copysign(sqrt(-log((1-abs(y))*sqrtπ)), y)
        isfinite(x) || return x
    end
    sqrtπhalf = sqrtπ * big(0.5)
    tol = 2eps(abs(x))
    n = 1
    while n < 100000 # Newton iterations
        Δx = sqrtπhalf * (erf(x) - y) * exp(x^2)
        x -= Δx
        abs(Δx) < tol && break
        n +=1
    end
    return x
end

############################
# Laplace
############################

function cdf_laplace_std(x) 
    return 0.5 - 0.5 * sign(x) * expm1(-abs(x))
end 

function invcdf_laplace_std(x)
    return sign(0.5 - x) * log1p(-2.0*abs(x - 0.5))
end 

function pdf_laplace_std(x)
    return 0.5 * exp(-abs(x))
end

function lpdf_laplace_std(x::T) where T<:Real
    return log(0.5) - abs(x)
end

function lpdf_laplace_std(x::AbstractVector{T}) where T<:Real
    d = size(x, 1)
    return d*log(0.5) - sum(abs, x)
end

function ∇lpdf_laplace_std(x::T) where T<:Real
    return -sign(x)
end

function ∇lpdf_laplace_std(x::AbstractVector{T}) where T<:Real
    return -sign.(x)
end

function randl(size)
    return rand(Laplace(), size)
end

randl()=rand(Laplace())
