From 194e88fa1b0f9f06659245be5569939b1738206d Mon Sep 17 00:00:00 2001
From: Guenther Raidl <raidl@ac.tuwien.ac.at>
Date: Mon, 15 Dec 2025 19:22:41 +0100
Subject: [PATCH 01/10] Fix: don't sample elements with weight 0

---
 src/sampling.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/sampling.jl b/src/sampling.jl
index 358092e9a..c07bbd470 100644
--- a/src/sampling.jl
+++ b/src/sampling.jl
@@ -596,6 +596,10 @@ function sample(rng::AbstractRNG, wv::AbstractWeights)
         i += 1
         cw += wv[i]
     end
+    if cw < t
+        # may happen with floating point weights due to numerical inaccuracies
+        return findlast(!iszero, wv) 
+    end
     return i
 end
 sample(wv::AbstractWeights) = sample(default_rng(), wv)

From 069b79f4fe11b03cb792ef366e4d3368e8c2de42 Mon Sep 17 00:00:00 2001
From: Guenther Raidl <raidl@ac.tuwien.ac.at>
Date: Tue, 16 Dec 2025 09:53:30 +0100
Subject: [PATCH 02/10] Alternative: replace sum(wv)

---
 src/sampling.jl | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/src/sampling.jl b/src/sampling.jl
index c07bbd470..1a7d8f4ea 100644
--- a/src/sampling.jl
+++ b/src/sampling.jl
@@ -586,7 +586,7 @@ Optionally specify a random number generator `rng` as the first argument
 function sample(rng::AbstractRNG, wv::AbstractWeights)
     1 == firstindex(wv) ||
         throw(ArgumentError("non 1-based arrays are not supported"))
-    wsum = sum(wv)
+    wsum = foldl(+, wv)  # instead of sum(wv) for avoiding numerical discrepancies with cw
     isfinite(wsum) || throw(ArgumentError("only finite weights are supported"))
     t = rand(rng) * wsum
     n = length(wv)
@@ -596,10 +596,6 @@ function sample(rng::AbstractRNG, wv::AbstractWeights)
         i += 1
         cw += wv[i]
     end
-    if cw < t
-        # may happen with floating point weights due to numerical inaccuracies
-        return findlast(!iszero, wv) 
-    end
     return i
 end
 sample(wv::AbstractWeights) = sample(default_rng(), wv)

From 6582fa0d501fa882f814c8f73cee67e33caecd14 Mon Sep 17 00:00:00 2001
From: Guenther Raidl <raidl@ac.tuwien.ac.at>
Date: Tue, 16 Dec 2025 11:49:38 +0100
Subject: [PATCH 03/10] Revert "Alternative: replace sum(wv)"

This reverts commit 069b79f4fe11b03cb792ef366e4d3368e8c2de42.
---
 src/sampling.jl | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/sampling.jl b/src/sampling.jl
index 1a7d8f4ea..c07bbd470 100644
--- a/src/sampling.jl
+++ b/src/sampling.jl
@@ -586,7 +586,7 @@ Optionally specify a random number generator `rng` as the first argument
 function sample(rng::AbstractRNG, wv::AbstractWeights)
     1 == firstindex(wv) ||
         throw(ArgumentError("non 1-based arrays are not supported"))
-    wsum = foldl(+, wv)  # instead of sum(wv) for avoiding numerical discrepancies with cw
+    wsum = sum(wv)
     isfinite(wsum) || throw(ArgumentError("only finite weights are supported"))
     t = rand(rng) * wsum
     n = length(wv)
@@ -596,6 +596,10 @@ function sample(rng::AbstractRNG, wv::AbstractWeights)
         i += 1
         cw += wv[i]
     end
+    if cw < t
+        # may happen with floating point weights due to numerical inaccuracies
+        return findlast(!iszero, wv) 
+    end
     return i
 end
 sample(wv::AbstractWeights) = sample(default_rng(), wv)

From c78063e84bff797dce9b273c19b4bb0c7322624c Mon Sep 17 00:00:00 2001
From: Guenther Raidl <raidl@ac.tuwien.ac.at>
Date: Tue, 16 Dec 2025 12:42:09 +0100
Subject: [PATCH 04/10] Test included

---
 sample.jl        | 52 ++++++++++++++++++++++++++++++++++++++++++++++++
 test/runtests.jl | 30 ++++++++++++++--------------
 test/sampling.jl |  9 +++++++++
 3 files changed, 76 insertions(+), 15 deletions(-)
 create mode 100644 sample.jl

diff --git a/sample.jl b/sample.jl
new file mode 100644
index 000000000..a21650996
--- /dev/null
+++ b/sample.jl
@@ -0,0 +1,52 @@
+using StatsBase
+using Random
+
+mask = Bool[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
+
+probs = Float32[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0437019, 0.0, 0.0, 0.04302464, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.039748967, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.040406376, 0.042578973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.040906563, 0.039586294, 0.0, 0.0, 0.0, 0.0, 0.0, 0.04302464, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.042357873, 0.04302464, 0.0, 0.039262936, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.040406376, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.040406376, 0.041919112, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.041484896, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.04057242, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.048750814, 0.0, 0.0, 0.04057242, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.042137925, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.040906563, 0.0, 0.0, 0.0, 0.0, 0.0, 0.04280123, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.041919112, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.04107466, 0.03942428, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0]
+
+
+function sample(rng::AbstractRNG, wv::AbstractWeights)
+    1 == firstindex(wv) ||
+        throw(ArgumentError("non 1-based arrays are not supported"))
+    wsum = sum(wv)  #foldl(+, wv) #sum(wv)
+    #println(wsum, " ",sum(wv))
+    isfinite(wsum) || throw(ArgumentError("only finite weights are supported"))
+    t = rand(rng) * wsum
+    n = length(wv)
+    i = 1
+    cw = wv[1]
+    while cw < t && i < n
+        i += 1
+        cw += wv[i]
+        # println("i=$i cw=$cw t=$t")
+    end
+    if cw < t && i >= n
+        println("catched $i")
+        # may happen with floating point weights due to numerical inaccuracies
+        return findlast(!iszero, wv) +10000
+    end
+    # println("selected $i")
+    return i
+end
+
+# probs = Float16[sqrt(2), π, 0.0]
+# probs = Float32[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0437019, 0.0, 0.0, 0.04302464, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.039748967, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.040406376, 0.042578973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.040906563, 0.039586294, 0.0, 0.0, 0.0, 0.0, 0.0, 0.04302464, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.042357873, 0.04302464, 0.0, 0.039262936, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.040406376, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.040406376, 0.041919112, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.041484896, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.04057242, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.048750814, 0.0, 0.0, 0.04057242, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.042137925, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.040906563, 0.0, 0.0, 0.0, 0.0, 0.0, 0.04280123, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.041919112, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.04107466, 0.03942428, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0]
+
+probs = Float32[ 0.0437019, 0.04302464, 0.039748967,  0.040406376, 0.042578973,  0.040906563, 0.039586294, 0.04302464, 0.042357873, 0.04302464, 0.039262936, 0.040406376, 0.040406376, 0.041919112, 0.041484896, 0.04057242, 0.048750814, 0.04057242, 0.042137925,  0.040906563, 0.04280123, 0.041919112, 0.04107466, 0.03942428, 0.0]
+
+probs = Float32[ 0.0437019, 0.04302464, 0.039748967,  0.040406376, 0.042578973,  0.040906563, 0.039586294, 0.04302464, 0.042357873, 0.04302464, 0.039262936, 0.040406376, 0.040406376, 0.041919112, 0.041484896, 0.04057242, 0.0]
+
+for i = 1:1000
+    i += 1
+    Random.seed!(i)
+    a = sample(Random.default_rng(), Weights(probs))
+    #println(a)
+    if a>=1000  #probs[a] == 0.0
+        println("Sampled $i index: ", a)
+    end
+end
+
+Random.seed!(10513378)
+a = sample(Random.default_rng(), Weights(probs))
+println("Result: ", a, length(probs))
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index da8212d9e..89ee47444 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,21 +4,21 @@ using LinearAlgebra
 using Random
 
 tests = ["ambiguous",
-         "weights",
-         "moments",
-         "scalarstats",
-         "deviation",
-         "cov",
-         "counts",
-         "ranking",
-         "empirical",
-         "hist",
-         "rankcorr",
-         "signalcorr",
-         "misc",
-         "pairwise",
-         "reliability",
-         "robust",
+        #  "weights",
+        #  "moments",
+        #  "scalarstats",
+        #  "deviation",
+        #  "cov",
+        #  "counts",
+        #  "ranking",
+        #  "empirical",
+        #  "hist",
+        #  "rankcorr",
+        #  "signalcorr",
+        #  "misc",
+        #  "pairwise",
+        #  "reliability",
+        #  "robust",
          "sampling",
          "wsampling",
          "statmodels",
diff --git a/test/sampling.jl b/test/sampling.jl
index cd2ff0960..9d29fa3b0 100644
--- a/test/sampling.jl
+++ b/test/sampling.jl
@@ -297,3 +297,12 @@ end
         end
     end
 end
+
+@testset "issue #982" begin
+    w = Float32[0.0437019, 0.04302464, 0.039748967,  0.040406376, 0.042578973, 
+        0.040906563, 0.039586294, 0.04302464, 0.042357873, 0.04302464, 0.039262936, 
+        0.040406376, 0.040406376, 0.041919112, 0.041484896, 0.04057242, 0.0]
+    Random.seed!(10513378)
+    s = sample(Random.default_rng(), Weights(w))
+    @test s == length(w) - 1
+end
\ No newline at end of file

From 54388a96b2ccb840b89de81f346a9671ec5ec15d Mon Sep 17 00:00:00 2001
From: Guenther Raidl <raidl@ac.tuwien.ac.at>
Date: Tue, 16 Dec 2025 12:45:59 +0100
Subject: [PATCH 05/10] Removed accidentally committed file

---
 sample.jl | 52 ----------------------------------------------------
 1 file changed, 52 deletions(-)
 delete mode 100644 sample.jl

diff --git a/sample.jl b/sample.jl
deleted file mode 100644
index a21650996..000000000
--- a/sample.jl
+++ /dev/null
@@ -1,52 +0,0 @@
-using StatsBase
-using Random
-
-mask = Bool[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
-
-probs = Float32[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0437019, 0.0, 0.0, 0.04302464, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.039748967, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.040406376, 0.042578973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.040906563, 0.039586294, 0.0, 0.0, 0.0, 0.0, 0.0, 0.04302464, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.042357873, 0.04302464, 0.0, 0.039262936, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.040406376, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.040406376, 0.041919112, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.041484896, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.04057242, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.048750814, 0.0, 0.0, 0.04057242, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.042137925, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.040906563, 0.0, 0.0, 0.0, 0.0, 0.0, 0.04280123, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.041919112, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.04107466, 0.03942428, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0]
-
-
-function sample(rng::AbstractRNG, wv::AbstractWeights)
-    1 == firstindex(wv) ||
-        throw(ArgumentError("non 1-based arrays are not supported"))
-    wsum = sum(wv)  #foldl(+, wv) #sum(wv)
-    #println(wsum, " ",sum(wv))
-    isfinite(wsum) || throw(ArgumentError("only finite weights are supported"))
-    t = rand(rng) * wsum
-    n = length(wv)
-    i = 1
-    cw = wv[1]
-    while cw < t && i < n
-        i += 1
-        cw += wv[i]
-        # println("i=$i cw=$cw t=$t")
-    end
-    if cw < t && i >= n
-        println("catched $i")
-        # may happen with floating point weights due to numerical inaccuracies
-        return findlast(!iszero, wv) +10000
-    end
-    # println("selected $i")
-    return i
-end
-
-# probs = Float16[sqrt(2), π, 0.0]
-# probs = Float32[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0437019, 0.0, 0.0, 0.04302464, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.039748967, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.040406376, 0.042578973, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.040906563, 0.039586294, 0.0, 0.0, 0.0, 0.0, 0.0, 0.04302464, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.042357873, 0.04302464, 0.0, 0.039262936, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.040406376, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.040406376, 0.041919112, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.041484896, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.04057242, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.048750814, 0.0, 0.0, 0.04057242, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.042137925, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.040906563, 0.0, 0.0, 0.0, 0.0, 0.0, 0.04280123, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.041919112, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.04107466, 0.03942428, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0]
-
-probs = Float32[ 0.0437019, 0.04302464, 0.039748967,  0.040406376, 0.042578973,  0.040906563, 0.039586294, 0.04302464, 0.042357873, 0.04302464, 0.039262936, 0.040406376, 0.040406376, 0.041919112, 0.041484896, 0.04057242, 0.048750814, 0.04057242, 0.042137925,  0.040906563, 0.04280123, 0.041919112, 0.04107466, 0.03942428, 0.0]
-
-probs = Float32[ 0.0437019, 0.04302464, 0.039748967,  0.040406376, 0.042578973,  0.040906563, 0.039586294, 0.04302464, 0.042357873, 0.04302464, 0.039262936, 0.040406376, 0.040406376, 0.041919112, 0.041484896, 0.04057242, 0.0]
-
-for i = 1:1000
-    i += 1
-    Random.seed!(i)
-    a = sample(Random.default_rng(), Weights(probs))
-    #println(a)
-    if a>=1000  #probs[a] == 0.0
-        println("Sampled $i index: ", a)
-    end
-end
-
-Random.seed!(10513378)
-a = sample(Random.default_rng(), Weights(probs))
-println("Result: ", a, length(probs))
\ No newline at end of file

From c3f6c74f1fe336f31311d4c1c744de8d01bdcc9b Mon Sep 17 00:00:00 2001
From: Guenther Raidl <raidl@ac.tuwien.ac.at>
Date: Tue, 16 Dec 2025 16:09:37 +0100
Subject: [PATCH 06/10] findlast avoided, test uses StableRNG

---
 src/sampling.jl  | 4 +++-
 test/sampling.jl | 4 ++--
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/sampling.jl b/src/sampling.jl
index c07bbd470..6f570678d 100644
--- a/src/sampling.jl
+++ b/src/sampling.jl
@@ -598,7 +598,9 @@ function sample(rng::AbstractRNG, wv::AbstractWeights)
     end
     if cw < t
         # may happen with floating point weights due to numerical inaccuracies
-        return findlast(!iszero, wv) 
+        while iszero(wv[i])
+            i -= 1
+        end
     end
     return i
 end
diff --git a/test/sampling.jl b/test/sampling.jl
index 9d29fa3b0..3307585d9 100644
--- a/test/sampling.jl
+++ b/test/sampling.jl
@@ -302,7 +302,7 @@ end
     w = Float32[0.0437019, 0.04302464, 0.039748967,  0.040406376, 0.042578973, 
         0.040906563, 0.039586294, 0.04302464, 0.042357873, 0.04302464, 0.039262936, 
         0.040406376, 0.040406376, 0.041919112, 0.041484896, 0.04057242, 0.0]
-    Random.seed!(10513378)
-    s = sample(Random.default_rng(), Weights(w))
+    rng = StableRNG(889858990530)
+    s = sample(rng, Weights(w))
     @test s == length(w) - 1
 end
\ No newline at end of file

From 56688510c44294cb04d34505f57519104491c240 Mon Sep 17 00:00:00 2001
From: Guenther Raidl <raidl@ac.tuwien.ac.at>
Date: Thu, 18 Dec 2025 17:26:27 +0100
Subject: [PATCH 07/10] Added warning in test if issue not triggered

---
 test/sampling.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/test/sampling.jl b/test/sampling.jl
index 3307585d9..ce008cde3 100644
--- a/test/sampling.jl
+++ b/test/sampling.jl
@@ -302,6 +302,10 @@ end
     w = Float32[0.0437019, 0.04302464, 0.039748967,  0.040406376, 0.042578973, 
         0.040906563, 0.039586294, 0.04302464, 0.042357873, 0.04302464, 0.039262936, 
         0.040406376, 0.040406376, 0.041919112, 0.041484896, 0.04057242, 0.0]
+    if sum(w) != 0.662413f0
+        @warn "Issue #982 not triggered and therefore resolution not tested \
+            as no appropriate SIMD support available."
+    end
     rng = StableRNG(889858990530)
     s = sample(rng, Weights(w))
     @test s == length(w) - 1

From c0b32f33e855aa04497c62c57db151b6cafcff42 Mon Sep 17 00:00:00 2001
From: Guenther Raidl <raidl@ac.tuwien.ac.at>
Date: Thu, 18 Dec 2025 17:38:32 +0100
Subject: [PATCH 08/10] Situation commented and warning adapted

---
 test/sampling.jl | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/test/sampling.jl b/test/sampling.jl
index ce008cde3..eb2e76ce5 100644
--- a/test/sampling.jl
+++ b/test/sampling.jl
@@ -302,7 +302,11 @@ end
     w = Float32[0.0437019, 0.04302464, 0.039748967,  0.040406376, 0.042578973, 
         0.040906563, 0.039586294, 0.04302464, 0.042357873, 0.04302464, 0.039262936, 
         0.040406376, 0.040406376, 0.041919112, 0.041484896, 0.04057242, 0.0]
-    if sum(w) != 0.662413f0
+    # It appears issue #982 only is triggered in conjunction with SIMD support
+    # under certain circumstances, in particular when sum(w) != 0.662413f0
+    # Without SIMD support, sum(w) == 0.66241294f0 and this test cannot check the
+    # resolution of the issue.
+    if sum(w) ∉ (0.662413f0, 0.66241294f0)
         @warn "Issue #982 not triggered and therefore resolution not tested \
             as no appropriate SIMD support available."
     end

From 011c073aced4800076ff11de2095c53a91df321c Mon Sep 17 00:00:00 2001
From: Guenther Raidl <raidl@ac.tuwien.ac.at>
Date: Thu, 18 Dec 2025 17:41:08 +0100
Subject: [PATCH 09/10] Warning message changed

---
 test/sampling.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/test/sampling.jl b/test/sampling.jl
index eb2e76ce5..e51418aa9 100644
--- a/test/sampling.jl
+++ b/test/sampling.jl
@@ -307,8 +307,7 @@ end
     # Without SIMD support, sum(w) == 0.66241294f0 and this test cannot check the
     # resolution of the issue.
     if sum(w) ∉ (0.662413f0, 0.66241294f0)
-        @warn "Issue #982 not triggered and therefore resolution not tested \
-            as no appropriate SIMD support available."
+        @warn "So far unrecognized value for sum(w) encountered."
     end
     rng = StableRNG(889858990530)
     s = sample(rng, Weights(w))

From 2393aac91628d1d5e9111566abdf65ef11c08383 Mon Sep 17 00:00:00 2001
From: Guenther Raidl <raidl@ac.tuwien.ac.at>
Date: Fri, 19 Dec 2025 10:59:47 +0100
Subject: [PATCH 10/10] Test simplified by providing sum

---
 src/sampling.jl  |  1 +
 test/runtests.jl | 30 +++++++++++++++---------------
 test/sampling.jl | 17 +++++++----------
 3 files changed, 23 insertions(+), 25 deletions(-)

diff --git a/src/sampling.jl b/src/sampling.jl
index 6f570678d..42a9fef00 100644
--- a/src/sampling.jl
+++ b/src/sampling.jl
@@ -1036,3 +1036,4 @@ wsample(rng::AbstractRNG, a::AbstractArray{T}, w::AbstractVector{<:Real}, dims::
 wsample(a::AbstractArray, w::AbstractVector{<:Real}, dims::Dims;
         replace::Bool=true, ordered::Bool=false) =
     wsample(default_rng(), a, w, dims; replace=replace, ordered=ordered)
+
diff --git a/test/runtests.jl b/test/runtests.jl
index 89ee47444..da8212d9e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,21 +4,21 @@ using LinearAlgebra
 using Random
 
 tests = ["ambiguous",
-        #  "weights",
-        #  "moments",
-        #  "scalarstats",
-        #  "deviation",
-        #  "cov",
-        #  "counts",
-        #  "ranking",
-        #  "empirical",
-        #  "hist",
-        #  "rankcorr",
-        #  "signalcorr",
-        #  "misc",
-        #  "pairwise",
-        #  "reliability",
-        #  "robust",
+         "weights",
+         "moments",
+         "scalarstats",
+         "deviation",
+         "cov",
+         "counts",
+         "ranking",
+         "empirical",
+         "hist",
+         "rankcorr",
+         "signalcorr",
+         "misc",
+         "pairwise",
+         "reliability",
+         "robust",
          "sampling",
          "wsampling",
          "statmodels",
diff --git a/test/sampling.jl b/test/sampling.jl
index e51418aa9..cc3fd3411 100644
--- a/test/sampling.jl
+++ b/test/sampling.jl
@@ -299,17 +299,14 @@ end
 end
 
 @testset "issue #982" begin
+    # Issue #982 was triggered in conjunction with SIMD support under certain circumstances, 
+    # in particular with the following weights when sum(w) == 0.662413f0.
+    # The following mimics this case.
     w = Float32[0.0437019, 0.04302464, 0.039748967,  0.040406376, 0.042578973, 
-        0.040906563, 0.039586294, 0.04302464, 0.042357873, 0.04302464, 0.039262936, 
-        0.040406376, 0.040406376, 0.041919112, 0.041484896, 0.04057242, 0.0]
-    # It appears issue #982 only is triggered in conjunction with SIMD support
-    # under certain circumstances, in particular when sum(w) != 0.662413f0
-    # Without SIMD support, sum(w) == 0.66241294f0 and this test cannot check the
-    # resolution of the issue.
-    if sum(w) ∉ (0.662413f0, 0.66241294f0)
-        @warn "So far unrecognized value for sum(w) encountered."
-    end
+    0.040906563, 0.039586294, 0.04302464, 0.042357873, 0.04302464, 0.039262936, 
+    0.040406376, 0.040406376, 0.041919112, 0.041484896, 0.04057242, 0.0]
     rng = StableRNG(889858990530)
-    s = sample(rng, Weights(w))
+    s = sample(rng, Weights(w, 0.662413f0))
     @test s == length(w) - 1
+    @test sample(rng, Weights([1, 2, 0, 0], 10000)) == 2  # another more trivial test
 end
\ No newline at end of file