Created
October 27, 2016 06:32
-
-
Save RedPointyJackson/b11bc010e9d41ae1fb5ac9eec762fb37 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
using ProgressMeter | |
using DataFrames | |
using Gadfly | |
############################ | |
# DATAPOINTS TO TAKE # | |
############################ | |
statistics = Int64(5e6) | |
Nbins = 100 | |
############################################ | |
# 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 # | |
############################################ | |
"True if p=(x,y) is in the unit circle" | |
is_in_circle(p) = p[1]^2+p[2]^2 < 1.0 | |
"Return a tuple (x,y) in the unit circle" | |
function random_in_circle(gen) | |
# Gen is supposed to return r∈[0,1) | |
# Adapt it to r∈[-1,1) | |
while true | |
x,y = 2*(gen()-0.5),2*(gen()-0.5) | |
is_in_circle((x,y)) && return (x,y) | |
end | |
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! | |
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_circle_substitution.csv",df_uniformly) | |
# Just plot. | |
info("Plotting!") | |
p1 = plot( | |
layer( | |
melt(df_uniformly,:r) | |
,x=:r | |
,y=:value | |
,color=:variable | |
,Geom.step | |
) | |
,layer(x->analytical_distance_distribution(x) | |
,0,2 | |
,style(line_width=10pt | |
,default_color=colorant"#DDDDDD")) | |
,Guide.xlabel("Distance between points") | |
,Guide.ylabel("Probability") | |
,Guide.title("Expected distance between points vs analytic") | |
) | |
df_deviation = deepcopy(df_uniformly) | |
expected = analytical_distance_distribution.(df_deviation[:r]) | |
df_deviation[:,2] = (df_deviation[:,2] - expected)./expected |> abs | |
df_deviation[:,3] = (df_deviation[:,3] - expected)./expected |> abs | |
df_deviation[:,4] = (df_deviation[:,4] - expected)./expected |> abs | |
p2 = plot(melt(df_deviation,:r) | |
,x=:r | |
,y=:value | |
,xgroup=:variable | |
,Geom.subplot_grid(Geom.line) | |
,Scale.y_log10 | |
,Guide.xlabel("Distance between points") | |
,Guide.ylabel("Deviation") | |
,Guide.title("Deviation from the expected distance (relative)") | |
) | |
draw(PNG("unit_circle_substitution.png",15cm,20cm),vstack(p1,p2)) | |
draw(PDF("unit_circle_substitution.pdf",15cm,20cm),vstack(p1,p2)) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment