Skip to content

Exact and sequential-importance-sampling fixed-margin samplers#6

Open
mkborregaard wants to merge 7 commits into
mainfrom
feature/exact-and-sis-sampling
Open

Exact and sequential-importance-sampling fixed-margin samplers#6
mkborregaard wants to merge 7 commits into
mainfrom
feature/exact-and-sis-sampling

Conversation

@mkborregaard

Copy link
Copy Markdown
Member

Adds two new fixed-margin matrix samplers alongside curveball, plus a weighted importance sampler, and updates dependencies. The public API (randomize_matrix!, matrixrandomizer) is extended with a method keyword; weighted sampling gets its own importance_sampler entry point.

New methods

  • exact — exact uniform sampling and counting (Miller & Harrison 2013). Ported from the authors' reference exact_pure.py; counts with BigInt, samples row-by-row from the cached count. Feasible for small matrices (rows + cols ≲ 100, or very sparse); the generator counts once and reuses it.
  • sis — the sequential importance sampling heuristic (Harrison & Miller 2013), uniform case. Builds a matrix column-by-column via a banded dynamic program over per-row Canfield odds and Gale–Ryser feasibility. Independent draws, scales to large matrices.

Both are exposed through the existing generator, which now dispatches a single draw on a per-method sampler type (CurveballSampler / SISSampler / ExactCounts) rather than branching on the enum.

Weighted importance sampling

importance_sampler(m, w) targets P*(z) ∝ ∏ w[i,j]^z[i,j] over fixed-margin binary matrices (uniform sis is the w ≡ 1 case). Each rand returns (matrix, logweight), so estimates reweight by exp(logweight)mean(exp(logweight)) estimates the weighted count κ, and a weight-weighted average estimates E[h(Z)].

  • Per-row odds are scaled by a factor vᵢ built from the elementary symmetric polynomials of the row's remaining weights (carried in logs).
  • Structural zeros (w[i,j]=0, e.g. a zero diagonal for directed graphs) forbid a one there; an uncompletable draw returns logweight = -Inf and is discarded.
  • Weights are Sinkhorn-canonicalised by default (canonicalize=false to opt out); P* is invariant to it.
  • Verified: the mean importance weight recovers brute-force κ within noise (including zero-diagonal targets), the uniform case reduces to the plain matrix count, and cv² stays low (≈0.002–0.11).

Curveball improvements

The curveball generator now warm-starts its chain from an independent SIS draw, so the first sample is decorrelated from the (often atypical) input matrix, and a trades keyword controls how many trades separate successive draws (thinning). Empirically the running chain samples uniformly (χ² ≈ df).

Other

  • Dependency bump (StatsBase 0.34); resolves to the latest releases.
  • Fixed a latent bug: MatrixGenerator's show defined a shadowing local instead of extending Base.show.
  • Tests: curveball (9), sis (6), exact (11), weighted sis (10) — all green, with brute-force κ/count cross-checks.

Not included (efficiency-only)

The paper's variance-based column ordering for the weighted sampler; the importance weights are already correct and low-variance without it.

🤖 Generated with Claude Code

mkborregaard and others added 7 commits June 20, 2026 21:54
Add two independent-draw alternatives to the curveball sampler, both from
Miller & Harrison, selected via the `method` keyword on `randomize_matrix!`
and `matrixrandomizer`:

- `exact`: exact uniform sampling via memoized counting of fixed-margin
  binary matrices (Miller & Harrison 2013). The generator counts once and
  reuses the result across draws.
- `sis`: sequential importance sampling, a fast near-uniform proposal
  (Harrison & Miller 2013), suited to large matrices. Implements the uniform
  case; weighted sampling is left as a future extension.

The dynamic programs are expressed as small dispatch-based functions over
bespoke state types (ExactProblem/ExactCounts/Cursor, SISResidual/
SISWorkspace/ColumnScorer). Verified against brute-force counts and a
uniformity check; the SIS proposal probability is validated by recovering
the matrix count via importance sampling.

Also bump StatsBase compat to 0.34 and the package version to 0.2.0.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Replace the per-column full re-sort with an O(ones-placed) repair of the
descending row order (each placed row drops by one, so it only slides past
the block holding its old value), tracked via an added inverse-permutation
field. Also clear only the touched band of the DP buffer each row instead
of the whole vector, and normalise with a reciprocal multiply. ~20% faster
per draw with identical output; importance-weight count estimates unchanged.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The curveball generator now warm-starts its chain from one independent SIS
draw instead of the observed matrix, so the first sample is decorrelated from
the input. A new `trades` keyword on randomize_matrix! and matrixrandomizer
sets how many trades each draw performs (default 5*ncols), controlling thinning
between successive samples.

Refactor the generator to dispatch a single draw on a per-method sampler type
(CurveballSampler/SISSampler/ExactCounts) rather than branching on the method
enum, keeping each method's draw logic in its own file.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The bare `show` defined a new RandomBooleanMatrices.show instead of extending
Base.show, so the custom display never fired. Qualify it and use print (no
trailing newline) so it shows as intended.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Generalise the SIS sampler to the Harrison & Miller (2013) weighted target
P*(z) ∝ ∏ w[i,j]^z[i,j], of which the uniform sampler is the w≡1 special case.
A WeightModel (UniformWeights / WeightMatrix) multiplies each row's odds of a
one by a factor v_i built from the elementary symmetric polynomials of the row's
remaining weights (carried in logs), leaving the dynamic program otherwise
unchanged. Each draw now carries its log importance weight log(∏w^z / Q*), and
the public importance_sampler(m, w) returns (matrix, logweight) pairs whose mean
weight estimates the weighted count κ.

Verified: the mean importance weight recovers brute-force κ within noise, the
uniform case reduces to the matrix count, and cv² stays low. Scope: strictly
positive weights; structural zeros, weight canonicalisation, and variance-based
column ordering are deferred (they affect efficiency, not correctness).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Allow nonnegative weights. A zero weight is a structural zero forbidding a one
at that position: the log-space symmetric polynomials and the v_i factor drop
those positions automatically, and a draw the zeros leave no way to complete is
returned with logweight = -Inf (importance weight zero) rather than erroring.
The v_i binomial term now uses the count of available (positive) positions in
the row's remaining columns.

Build the proposal from the Sinkhorn-canonical w̄ (eq. 16) by default
(canonicalize=false to opt out); P* is invariant to it, so the importance
weights are unchanged. Verified on zero-diagonal targets (κ recovered, 0%
rejections) and that canonicalize on/off give the same estimate.

Also make the exact m3!=m4 test deterministic and robust (it relied on an
unseeded generator and a large enough support).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
@codecov

codecov Bot commented Jun 21, 2026

Copy link
Copy Markdown

Codecov Report

❌ Patch coverage is 97.90576% with 8 lines in your changes missing coverage. Please review.
✅ Project coverage is 97.31%. Comparing base (798b019) to head (8a1349e).

Files with missing lines Patch % Lines
src/RandomBooleanMatrices.jl 86.66% 4 Missing ⚠️
src/weights.jl 97.29% 2 Missing ⚠️
src/exact.jl 98.97% 1 Missing ⚠️
src/sis.jl 99.31% 1 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##             main       #6      +/-   ##
==========================================
+ Coverage   90.78%   97.31%   +6.52%     
==========================================
  Files           2        6       +4     
  Lines          76      447     +371     
==========================================
+ Hits           69      435     +366     
- Misses          7       12       +5     
Flag Coverage Δ
unittests 97.31% <97.90%> (+6.52%) ⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant