using Random
using HaltonSequences, Primes

function rand_points(n::Int, d::Int=2)
    x = rand(n, d)
    x = 2 .* x .- 1
    return x
end

function halton_points(n::Int, d::Int=2, P = primes(d^2+1)[1:d])
    P = shuffle!(P)[1:d]
    x = hcat( (Halton(P[k])[1:n] for k=1:d)... )
    x = 2 .* x .- 1
    return x
end

function grid_points(n::Int, d::Int=2)
    m = ceil(Int, n^(1/d))
    x = zeros(m^d, d)
    for i = 0:m^d-1
        j = [ (i÷m^(k-1)) % m for k=1:d ] 
        x[i+1,:] = j ./ m 
    end
    x = 2 .* x .- 1
    return x
end

function as_tuples(z::Matrix)
    n = size(z, 1)
    return [ tuple(z[i,:] ... ) for i=1:n]
end

function gauss_proxy(z::Matrix; μ::Number=0.0, σ::Number=1.0)
    z2 = sum(z .^ 2, dims=2)[:, 1]
    f = x -> - ( (sign(x) * sqrt( 4/pi * mean( z2 .< x^2 ) ) - μ ) / σ )^2 / 2
    return LogDensity(f)
end

function gauss_proxy(z::Matrix, nsamples::Vector{Int}; μ::Number=0.0, σ::Number=1.0)
    fs = [ gauss_proxy(z[1:n, :]; μ, σ) for n=nsamples ]  
    return MultilevelLogDensity(fs)
end


function gaussmix_proxy(z::Matrix, μ::Vector=[-.5, .5], σ::Vector=[.1, .1])
    fs = [ sample_gauss_proxy(z, m, s) for (m,s) in zip(μ, σ) ]
    return LogDensity( x -> log( sum(exp.([logdensity(f, x) for f in fs]))) )
end