Commit 837e2ee6 authored by Dahua Lin's avatar Dahua Lin

new benchmark script for samplers (based on BenchmarkLite)

parent c8f01195
# Benchmarking samplers
using BenchmarkLite
using Distributions
### define benchmark tasks
type UnivariateSamplerRun{Spl} <: Proc end
const batch_unit = 1000
Base.isvalid(::UnivariateSamplerRun, cfg) = true
Base.length(p::UnivariateSamplerRun, cfg) = batch_unit
Base.string{Spl}(p::UnivariateSamplerRun{Spl}) = getname(Spl)
Base.start{Spl}(p::UnivariateSamplerRun{Spl}, cfg) = getsampler(Spl, cfg)
Base.done(p::UnivariateSamplerRun, cfg, s) = nothing
getsampler{Spl<:Sampleable}(::Type{Spl}, cfg) = Spl(cfg...)
function Base.run(p::UnivariateSamplerRun, cfg, s)
for i = 1:batch_unit
rand(s)
end
end
make_procs(spltypes...) = Proc[UnivariateSamplerRun{T}() for T in spltypes]
### specific benchmarking program
## categorical
import Distributions: AliasTable, CategoricalDirectSampler
getsampler(::Type{CategoricalDirectSampler}, k::Int) = CategoricalDirectSampler(fill(1/k, k))
getsampler(::Type{AliasTable}, k::Int) = AliasTable(fill(1/k, k))
getname(::Type{CategoricalDirectSampler}) = "direct"
getname(::Type{AliasTable}) = "alias"
benchmark_categorical() = (
make_procs(CategoricalDirectSampler, AliasTable),
"K", 2 .^ (1:12))
## binomial
import Distributions: BinomialRmathSampler, BinomialAliasSampler
import Distributions: BinomialGeomSampler, BinomialTPESampler, BinomialPolySampler
getname(::Type{BinomialRmathSampler}) = "rmath"
getname(::Type{BinomialAliasSampler}) = "alias"
getname(::Type{BinomialGeomSampler}) = "Geom"
getname(::Type{BinomialTPESampler}) = "BTPE"
getname(::Type{BinomialPolySampler}) = "Poly"
benchmark_binomial() = (
make_procs(BinomialRmathSampler,
BinomialAliasSampler,
BinomialGeomSampler,
BinomialTPESampler,
BinomialPolySampler),
"(n, p)", vec([(n, p) for p in [0.3, 0.5, 0.9], n in 2.^(1:12)]))
## poisson
import Distributions: PoissonRmathSampler, PoissonCountSampler, PoissonADSampler
getname(::Type{PoissonRmathSampler}) = "rmath"
getname(::Type{PoissonCountSampler}) = "count"
getname(::Type{PoissonADSampler}) = "AD"
Base.isvalid(::UnivariateSamplerRun{PoissonADSampler}, mu) = (mu >= 5.0)
benchmark_poisson() = (
make_procs(PoissonRmathSampler, PoissonCountSampler, PoissonADSampler),
"μ", [0.2, 0.5, 1.0, 2.0, 5.0, 10.0, 15.0, 20.0, 30.0])
## exponential
import Distributions: ExponentialSampler, ExponentialLogUSampler
getname(::Type{ExponentialSampler}) = "base"
getname(::Type{ExponentialLogUSampler}) = "logu"
benchmark_exponential() = (
make_procs(ExponentialSampler, ExponentialLogUSampler),
"scale", [1.0])
## gamma
import Distributions: GammaRmathSampler, GammaMTSampler
getname(::Type{GammaRmathSampler}) = "rmath"
getname(::Type{GammaMTSampler}) = "MT"
benchmark_gamma() = (
make_procs(GammaRmathSampler, GammaMTSampler),
"(α, scale)", [(α, 1.0) for α in [0.5, 1.0, 2.0, 5.0, 20.0]])
### main
const dnames = ["categorical",
"binomial",
"poisson",
"exponential",
"gamma"]
function printhelp()
println("Require exactly one argument. Usage:")
println()
println(" julia <path>/sampler.jl <distrname>")
println()
println(" <distrname> can be:")
println()
for dname in dnames
println(" - $dname")
end
println()
end
function getarg(args)
if length(args) == 1
dname = args[1]
if dname in dnames
return dname
else
printhelp()
exit(1)
end
else
printhelp()
exit(1)
end
end
dname = getarg(ARGS)
function do_benchmark(dname; verbose::Int=2)
(procs, cfghead, cfgs) =
dname == "categorical" ? benchmark_categorical() :
dname == "binomial" ? benchmark_binomial() :
dname == "poisson" ? benchmark_poisson() :
dname == "exponential" ? benchmark_exponential() :
dname == "gamma" ? benchmark_gamma() :
error("benchmarking function for $dname has not been implemented.")
r = run(procs, cfgs; duration=0.2, verbose=verbose)
println()
show(r; unit=:mps, cfghead=cfghead)
end
do_benchmark(dname)
println()
# BinomialRmathSampler is non-exported, which is mainly used
# for benchmarking or testing purpose
#
immutable BinomialRmathSampler <: Sampleable{Univariate,Discrete}
n::Int
prob::Float64
......
......@@ -7,3 +7,9 @@ end
rand(s::ExponentialSampler) = s.scale * randexp()
immutable ExponentialLogUSampler <: Sampleable{Univariate,Continuous}
scale::Float64
end
rand(s::ExponentialLogUSampler) = -s.scale * log(rand())
immutable GammaRmathSampler <: Sampleable{Univariate,Continuous}
α::Float64
scale::Float64
end
rand(s::GammaRmathSampler) =
ccall((:rgamma, "libRmath-julia"), Float64, (Float64, Float64), s.α, s.scale)
# A simple method for generating gamma variables - Marsaglia and Tsang (2000)
# http://www.cparity.com/projects/AcmClassification/samples/358414.pdf
# Page 369
......
immutable PoissonRmathSampler <: Sampleable{Univariate,Discrete}
mu::Float64
end
rand(s::PoissonRmathSampler) =
int(ccall((:rpois, "libRmath-julia"), Float64, (Float64,), s.mu))
function poissonpvec(μ::Float64, n::Int)
# Poisson probabilities, from 0 to n
pv = Array(Float64, n+1)
......
using Distributions
# config
Ks = [3, 10, 30, 100, 1000, 10000]
n = 10^5
function run_sample(sampler, n::Int)
for i = 1:n
rand(sampler)
end
end
macro bench_sampler(construct_expr)
quote
s = $(esc(construct_expr)) # construct sampler
tname = typeof(s)
# warm up
run_sample(s, 100)
et = @elapsed run_sample(s, n)
mps = n * 1.0e-6 / et
@printf(" %-28s : elapsed = %12.6f sec | %8.4f MPS\n",
tname, et, mps)
end
end
for K in Ks
println("K = $K")
println("-----------------------")
p = fill(1/K, K)
pu = fill(uint64(1), K)
@bench_sampler DiscreteUniform(1, K) # samples generated by rand(1:K)
@bench_sampler Categorical(p)
# @bench_sampler Distributions.DiscreteDistributionTable(p)
# @bench_sampler Distributions.huffman(1:K, pu)
@bench_sampler Distributions.AliasTable(p)
println()
end
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment