Draw curveball trades among non-empty columns only

Project.toml

@@ -2,7 +2,7 @@ name = "RandomBooleanMatrices"
 uuid = "9ae346a0-3d16-5633-ad70-ddb60ab77eac"
 license = "MIT"
 authors = ["mkborregaard <mkborregaard@snm.ku.dk"]
-version = "0.2.0"
+version = "0.2.1"
 
 [deps]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

README.md

@@ -51,6 +51,9 @@ m3 = rand!(rmg)
     `sis` draw, so even the first sample is decorrelated from the input matrix, and
     the `trades` keyword sets how many trades separate successive draws — raise it
     to thin out the correlation between samples.
+    Trades are drawn among the non-empty columns only (empty columns can never
+    trade), and `trades` defaults to five per non-empty column, so matrices with
+    many empty columns mix as fast as the same matrix without them.
   * `sis` — sequential importance sampling (Harrison & Miller 2013). Each call
     draws an **independent** matrix from a fast, near-uniform approximation of the
     fixed-margin distribution. Scales to large matrices and is the natural choice

src/RandomBooleanMatrices.jl

@@ -14,7 +14,7 @@ include("sis.jl")
 @enum matrixrandomizations curveball exact sis
 
 """
-    randomize_matrix!(m [,rng]; method = curveball, trades = 5 * size(m, 2))
+    randomize_matrix!(m [,rng]; method = curveball, trades)
 
 Randomize the sparse boolean Matrix `m` in place while maintaining its row and
 column sums. The algorithm is chosen with `method`:
@@ -30,9 +30,12 @@ column sums. The algorithm is chosen with `method`:
     step is repeated on every call, so prefer [`matrixrandomizer`](@ref), which
     counts once and reuses the result.
 
-`trades` applies only to `curveball` and sets how far the chain moves per call.
+`trades` applies only to `curveball` and sets how far the chain moves per call. It
+defaults to five per non-empty column: empty columns can never trade, so trades
+are drawn among the non-empty columns only, and a sparse matrix with many empty
+columns mixes as fast as the same matrix without them.
 """
-function randomize_matrix!(m::SparseMatrixCSC{Bool, Int}, rng = Random.GLOBAL_RNG; method::matrixrandomizations = curveball, trades::Int = 5 * size(m, 2))
+function randomize_matrix!(m::SparseMatrixCSC{Bool, Int}, rng = Random.GLOBAL_RNG; method::matrixrandomizations = curveball, trades::Int = _defaulttrades(m))
     method == curveball && return _curveball!(m, rng, trades)
     method == sis       && return _sis!(m, rng)
     method == exact     && return _exact!(m, rng)
@@ -61,13 +64,13 @@ function _sampler!(method::matrixrandomizations, sm::SparseMatrixCSC{Bool, Int},
     method == exact && return _exact_counts(_margins(sm)...)
     if method == curveball
         _sis!(sm, rng)
-        return CurveballSampler(something(trades, 5 * size(sm, 2)))
+        return CurveballSampler(something(trades, _defaulttrades(sm)))
     end
     error("undefined method")
 end
 
 """
-    matrixrandomizer(m [,rng]; method = curveball, trades = 5 * size(m, 2))
+    matrixrandomizer(m [,rng]; method = curveball, trades)
 
 Create a matrix generator that returns a random boolean matrix every time it is
 called, maintaining row and column sums. Non-boolean input matrices are
@@ -78,7 +81,8 @@ For `method = exact` the matrix count is computed once here and reused by every
 draw. For `method = curveball` the chain is warm-started from an independent `sis`
 draw, so even the first sample is decorrelated from `m`; `trades` then sets how
 many curveball trades separate successive draws — raise it to reduce correlation
-between samples.
+between samples. It defaults to five per non-empty column, as for
+[`randomize_matrix!`](@ref).
 
 # Examples
 ```

src/curveball.jl

@@ -53,14 +53,25 @@ function _interdif!(v1, v2, inter, dif)
    ndiff, nshared
 end
 
-function _curveball!(m::SparseMatrixCSC{Bool, Int}, rng = Random.GLOBAL_RNG, trades::Int = 5 * size(m, 2))
-   R, C = size(m)
+# The columns that can take part in a trade. A pair involving an empty column is a
+# no-op and column sums never change, so this set is fixed along the chain; drawing
+# pairs only from it leaves the target distribution unchanged, and stops a sparse
+# matrix with many empty columns from spending nearly all its trades on nothing.
+_tradingcolumns(m::SparseMatrixCSC) = findall(>(0), diff(SparseArrays.getcolptr(m)))
+
+# The default number of trades per call: five per column that can trade.
+_defaulttrades(m::SparseMatrixCSC) = 5 * length(_tradingcolumns(m))
+
+function _curveball!(m::SparseMatrixCSC{Bool, Int}, rng = Random.GLOBAL_RNG, trades::Int = _defaulttrades(m))
+   cols = _tradingcolumns(m)
+   C = length(cols)
+   C < 2 && return m
    mcs = min(2maximum(diff(m.colptr)), size(m, 1))
    not_shared, shared = Vector{Int}(undef, mcs), Vector{Int}(undef, mcs)
    newa, newb = Vector{Int}(undef, mcs), Vector{Int}(undef, mcs)
 
    for rep ∈ 1:trades
-	   A, B = rand(rng, 1:C,2)
+	   A, B = cols[rand(rng, 1:C)], cols[rand(rng, 1:C)]
 
       # use views directly into the sparse matrix to avoid copying
 	   a, b = view(m.rowval, m.colptr[A]:m.colptr[A+1]-1), view(m.rowval, m.colptr[B]:m.colptr[B+1]-1)

test/runtests.jl

@@ -66,6 +66,25 @@ end
     # the generator warm-starts from an independent draw and preserves margins
     csm5, rsm5 = margins(m5)
     @test margins(rand(matrixrandomizer(m5, method = curveball, trades = 10))) == (csm5, rsm5)
+
+    # a sparse matrix with many empty columns still mixes: trades are drawn among
+    # the non-empty columns, so the empty ones do not use them up
+    sp = spzeros(Bool, 6, 20_000)
+    for (k, j) in enumerate(randperm(MersenneTwister(5), 20_000)[1:60])
+        sp[mod1(k, 6), j] = true
+        sp[mod1(k + 2, 6), j] = true
+    end
+    @test matrixrandomizer(sp, method = curveball).state.trades == 5 * 60
+    spcsm, sprsm = margins(sp)
+    rmg = matrixrandomizer(sp, MersenneTwister(6), method = curveball)
+    draws = [rand(rmg) for _ in 1:20]
+    @test all(d -> margins(d) == (spcsm, sprsm), draws)
+    @test minimum(i -> count(draws[i] .& .!draws[i + 1]), 1:19) >= 10
+
+    # with fewer than two non-empty columns nothing can trade
+    single = spzeros(Bool, 4, 50)
+    single[2, 7] = true
+    @test randomize_matrix!(copy(single), method = curveball) == single
 end
 
 @testset "sis" begin