stack_exchange_randoms_2.jl

# This time with a circumference
using ProgressMeter
using DataFrames
using Gadfly

Gadfly.push_theme(:default)


############################
#    DATAPOINTS TO TAKE    #
############################
statistics = Int64(5e6)
Nbins = 200


############################################
# Seed random generators, wrap C generator #
############################################

srand(42) # Mersenne twister (Julia's default)

ccall(("srand","libc"),Void,(Int64,),42) # C randoms
crand() = ccall(("rand","libc"),Int64,())/2^31

# rand() as LCG
STATE=42
function lcgrand()
    next = ((STATE * 1103515245) + 12345) & 0x7fffffff
    global STATE = next
    return STATE/2^31
end


############################################
#        Unit Circle w sub               #
############################################
"Return a tuple (x,y) in the unit circumference"
function random_in_circle(gen)
    # Gen is supposed to return r∈[0,1)
    # Adapt it to r∈[-1,1)
    θ = 2π*gen()
    return cos(θ),sin(θ)
end

"""
Return the pairwise distances distribution of
`pointlist::Vector{Tuple}` in linspace(0,2,bins).
"""
function distance_distribution(pointlist,bins)
    N = length(pointlist)
    result = zeros(bins)
    # Dont pick the N² points. Too slow. Pick N distances aprox.
    # The results are the same picking all the distances, so I'm
    # choosing this path.
    rand_is = shuffle(1:N)[1:floor(Int64,sqrt(N))]
    rand_js = shuffle(1:N)[1:floor(Int64,sqrt(N))]
    P = Progress(N,desc="Binning distances...")
    for i ∈ rand_is, j ∈ rand_js
        if i != j
            x1,y1 = pointlist[i]
            x2,y2 = pointlist[j]
            d = sqrt((x1-x2)^2+(y1-y2)^2)
            for k in 0:bins-1
                minbound = 0 + k*2/bins
                maxbound = minbound + 1/bins
                if minbound ≤ d < maxbound
                    result[k+1] += 1
                    break
                end
            end
        end
        next!(P)
    end
    area = sum(result)*2/bins
    return result/area # normalize!
    println()
end

analytical_distance_distribution(d) = d/π*(4*atan(√(4-d^2)/d)-d*√(4-d^2))

"Substitute a point in the circle for another one"
function substitute_in_circle!(points, gen)
    N = length(points)
    idx = floor(Int64,gen()*N) + 1
    neu = random_in_circle(gen)
    points[idx] = neu
end

# Generate the lists of random points and substitute every one of
# them. Then, bin.
df_uniformly = DataFrame()
df_uniformly[:r] = linspace(0,2,Nbins)

info("\nComputing for glibc randoms.")
list_crand = [random_in_circle(crand) for i in 1:statistics]
for i in 1:statistics
    substitute_in_circle!(list_crand,crand)
end
df_uniformly[Symbol("glibc rand")]  = distance_distribution(list_crand,Nbins)

info("Computing for the LCG.")
list_lcgrand = [random_in_circle(lcgrand) for i in 1:statistics]
for i in 1:statistics
    substitute_in_circle!(list_lcgrand,lcgrand)
end
df_uniformly[Symbol("LCG")] = distance_distribution(list_lcgrand,Nbins)

info("Computing for the Mersenne Twister.")
list_mersenne = [random_in_circle(rand) for i in 1:statistics]
for i in 1:statistics
    substitute_in_circle!(list_mersenne,rand)
end
df_uniformly[Symbol("Mersenne Twister")]  = distance_distribution(list_mersenne,Nbins)

writetable("unit_circumference_substitution.csv",df_uniformly)

# Just plot.
info("Plotting!")
p1 = plot(melt(df_uniformly,:r)
          ,x=:r
          ,y=:value
          ,color=:variable
          ,Geom.step
          ,Guide.xlabel("Distance between points")
          ,Guide.ylabel("PDF")
          ,Guide.title("Empirical distance between points (circumference)")
          )


draw(PNG("unit_circumference_substitution.png",15cm,10cm),p1)
draw(PDF("unit_circumference_substitution.pdf",15cm,10cm),p1)