From 87b066f724b80f6be1ab342b8410d36e77dd43bc Mon Sep 17 00:00:00 2001 From: Adriano Meligrana <68152031+Tortar@users.noreply.github.com> Date: Sat, 16 Aug 2025 21:58:44 +0200 Subject: [PATCH 1/5] Improve performance of choose --- src/UnweightedSamplingMulti.jl | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/src/UnweightedSamplingMulti.jl b/src/UnweightedSamplingMulti.jl index 29979477..63f7843b 100644 --- a/src/UnweightedSamplingMulti.jl +++ b/src/UnweightedSamplingMulti.jl @@ -165,26 +165,27 @@ end @inline function choose(rng, n, p) z = exp(n*log1p(-p)) t = rand(rng, Uniform(z, 1.0)) + nt = t/z s = n*p q = 1-p - x = z + z*s/q - x > t && return 1 + x = 1 + s/q + x > nt && return 1 s *= (n-1)*p q *= 1-p - x += (s*z/q)/2 - x > t && return 2 + x += s/(q*2) + x > nt && return 2 s *= (n-2)*p q *= 1-p - x += (s*z/q)/6 - x > t && return 3 + x += s/(q*6) + x > nt && return 3 s *= (n-3)*p q *= 1-p - x += (s*z/q)/24 - x > t && return 4 + x += s/(q*24) + x > nt && return 4 s *= (n-4)*p q *= 1-p - x += (s*z/q)/120 - x > t && return 5 + x += s/(q*120) + x > nt && return 5 return quantile(Binomial(n, p), t) end @@ -274,3 +275,4 @@ function ordvalue(s::MultiOrdAlgRSWRSKIPSampler) return s.value[sortperm(s.ord)] end end + From 5c4245e12e090d5ebc3dae0969a5b25f3f0d148c Mon Sep 17 00:00:00 2001 From: Adriano Meligrana <68152031+Tortar@users.noreply.github.com> Date: Sun, 17 Aug 2025 00:51:24 +0200 Subject: [PATCH 2/5] Update UnweightedSamplingMulti.jl --- src/UnweightedSamplingMulti.jl | 43 ++++++++++++++++++---------------- 1 file changed, 23 insertions(+), 20 deletions(-) diff --git a/src/UnweightedSamplingMulti.jl b/src/UnweightedSamplingMulti.jl index 63f7843b..c7833fce 100644 --- a/src/UnweightedSamplingMulti.jl +++ b/src/UnweightedSamplingMulti.jl @@ -162,30 +162,32 @@ function recompute_skip!(s::MultiAlgRSWRSKIPSampler, n) return s end +macro quantile_fast(k) + block = Expr(:block) + initial_code = quote + $(esc(:s)) = $(esc(:n)) * $(esc(:p)) + $(esc(:q)) = 1 - $(esc(:p)) + $(esc(:x)) = 1 + $(esc(:s)) / $(esc(:q)) + $(esc(:x)) > $(esc(:nt)) && return 1 + end + append!(block.args, initial_code.args) + for i in 2:k + iteration_code = quote + $(esc(:s)) *= ($(esc(:n)) - $i) * $(esc(:p)) + $(esc(:q)) *= 1 - $(esc(:p)) + $(esc(:x)) += $(esc(:s)) / ($(esc(:q)) * $(factorial(i))) + $(esc(:x)) > $(esc(:nt)) && return $i + end + append!(block.args, iteration_code.args) + end + return block +end + @inline function choose(rng, n, p) z = exp(n*log1p(-p)) t = rand(rng, Uniform(z, 1.0)) nt = t/z - s = n*p - q = 1-p - x = 1 + s/q - x > nt && return 1 - s *= (n-1)*p - q *= 1-p - x += s/(q*2) - x > nt && return 2 - s *= (n-2)*p - q *= 1-p - x += s/(q*6) - x > nt && return 3 - s *= (n-3)*p - q *= 1-p - x += s/(q*24) - x > nt && return 4 - s *= (n-4)*p - q *= 1-p - x += s/(q*120) - x > nt && return 5 + @quantile_fast(8) return quantile(Binomial(n, p), t) end @@ -276,3 +278,4 @@ function ordvalue(s::MultiOrdAlgRSWRSKIPSampler) end end + From 77cd0af448d128a591c3cfb1a199672216cb9135 Mon Sep 17 00:00:00 2001 From: Adriano Meligrana <68152031+Tortar@users.noreply.github.com> Date: Sun, 17 Aug 2025 00:51:44 +0200 Subject: [PATCH 3/5] Update UnweightedSamplingMulti.jl --- src/UnweightedSamplingMulti.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/UnweightedSamplingMulti.jl b/src/UnweightedSamplingMulti.jl index c7833fce..f16bf960 100644 --- a/src/UnweightedSamplingMulti.jl +++ b/src/UnweightedSamplingMulti.jl @@ -277,5 +277,3 @@ function ordvalue(s::MultiOrdAlgRSWRSKIPSampler) return s.value[sortperm(s.ord)] end end - - From ae6610747d035184735ce7b8d59e11fa51bcbe68 Mon Sep 17 00:00:00 2001 From: Adriano Meligrana <68152031+Tortar@users.noreply.github.com> Date: Sun, 17 Aug 2025 01:01:21 +0200 Subject: [PATCH 4/5] Update UnweightedSamplingMulti.jl --- src/UnweightedSamplingMulti.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/UnweightedSamplingMulti.jl b/src/UnweightedSamplingMulti.jl index f16bf960..514ee015 100644 --- a/src/UnweightedSamplingMulti.jl +++ b/src/UnweightedSamplingMulti.jl @@ -166,15 +166,15 @@ macro quantile_fast(k) block = Expr(:block) initial_code = quote $(esc(:s)) = $(esc(:n)) * $(esc(:p)) - $(esc(:q)) = 1 - $(esc(:p)) - $(esc(:x)) = 1 + $(esc(:s)) / $(esc(:q)) + $(esc(:q)) = 1. - $(esc(:p)) + $(esc(:x)) = 1. + $(esc(:s)) / $(esc(:q)) $(esc(:x)) > $(esc(:nt)) && return 1 end append!(block.args, initial_code.args) for i in 2:k iteration_code = quote $(esc(:s)) *= ($(esc(:n)) - $i) * $(esc(:p)) - $(esc(:q)) *= 1 - $(esc(:p)) + $(esc(:q)) *= 1. - $(esc(:p)) $(esc(:x)) += $(esc(:s)) / ($(esc(:q)) * $(factorial(i))) $(esc(:x)) > $(esc(:nt)) && return $i end @@ -277,3 +277,4 @@ function ordvalue(s::MultiOrdAlgRSWRSKIPSampler) return s.value[sortperm(s.ord)] end end + From f9e797920b650dc587e0e66fa95b7d732dcc5f01 Mon Sep 17 00:00:00 2001 From: Adriano Meligrana <68152031+Tortar@users.noreply.github.com> Date: Sun, 17 Aug 2025 01:21:52 +0200 Subject: [PATCH 5/5] Update UnweightedSamplingMulti.jl --- src/UnweightedSamplingMulti.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/UnweightedSamplingMulti.jl b/src/UnweightedSamplingMulti.jl index 514ee015..691d52a3 100644 --- a/src/UnweightedSamplingMulti.jl +++ b/src/UnweightedSamplingMulti.jl @@ -164,21 +164,21 @@ end macro quantile_fast(k) block = Expr(:block) - initial_code = quote + firstv = quote $(esc(:s)) = $(esc(:n)) * $(esc(:p)) $(esc(:q)) = 1. - $(esc(:p)) $(esc(:x)) = 1. + $(esc(:s)) / $(esc(:q)) $(esc(:x)) > $(esc(:nt)) && return 1 end - append!(block.args, initial_code.args) + append!(block.args, firstv.args) for i in 2:k - iteration_code = quote + nextv = quote $(esc(:s)) *= ($(esc(:n)) - $i) * $(esc(:p)) $(esc(:q)) *= 1. - $(esc(:p)) $(esc(:x)) += $(esc(:s)) / ($(esc(:q)) * $(factorial(i))) $(esc(:x)) > $(esc(:nt)) && return $i end - append!(block.args, iteration_code.args) + append!(block.args, nextv.args) end return block end @@ -278,3 +278,4 @@ function ordvalue(s::MultiOrdAlgRSWRSKIPSampler) end end +