Press n or j to go to the next uncovered block, b, p or k for the previous block.
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 | 2x 313710866x 313710866x 313710866x 4450x 4450x 4450x 4450x 4450x 2x 584x 584x 2225x 2225x 2225x 2225x 2225x 2225x 584x 584x 2x 2x 4659x 4659x 4659x 2x 4659x 4659x 18636x 18636x 18636x 18636x 18636x 18636x 18636x 4659x 4659x 4659x 4659x 4659x 4659x 2x 156855433x 156855433x 156855433x 156855433x 156855433x 156855433x 156855433x 156855433x 156855433x 156855433x 156855433x 156855433x 156855433x 156855433x 156855433x 156855433x 156855433x 156855433x 2x 77896432x 77896432x 77896432x 77896432x 2x 850126x 850126x 850124x 850124x 850124x 850124x 850124x 850124x 850124x 1062569x 1062569x 1062569x 850124x 1x 850126x 2x 2x 756634x 756634x 756634x 716593x 716593x 40038x 756634x 756634x 625718x 625718x 625718x 625718x 625718x 625718x 625718x 625718x 625718x 1594736x 1594736x 1594736x 1594736x 1594736x 1594736x 1594736x 625718x 625718x 7180x 7180x 7180x 7180x 130913x 130913x 130913x 130913x 130913x 130913x 130913x 130913x 130913x 130913x 130913x 130913x 130913x 130913x 130913x 130913x 130913x 130913x 130913x 130913x 160402x 160402x 160402x 160402x 160402x 42913x 42913x 42913x 37508x 160402x 5690x 5690x 5690x 11689x 5999x 5999x 5999x 5999x 49197x 160402x 26539x 26539x 26539x 26539x 13407x 26539x 13074x 13074x 26539x 13774x 13774x 22658x 22658x 22658x 160402x 12037x 1795x 1795x 1795x 1795x 1795x 1795x 1795x 1795x 1795x 1795x 1795x 1795x 1795x 718x 718x 130913x 2x 60x 60x 60x 2x 2x 2x 12562273x 6981806x 6981806x 2x 1368434x 1368434x 1364262x 1368434x 1368434x 120022x 120022x 120022x 120022x 120022x 120022x 120022x 120022x 120022x 780124x 780124x 780124x 780124x 120022x 120022x 120022x 2x 2x 1244239x 1244239x 1244239x 1244239x 1244239x 1244239x 1244239x 1244239x 1244239x 1244239x 1244239x 1244239x 1244239x 1244239x 1244239x 1244239x 1244239x 1244239x 1244239x 1244239x 1244239x 2162986x 2162986x 2162986x 2162986x 2162986x 1895754x 1895754x 1895754x 1895754x 1895754x 1895754x 1895754x 1895754x 724840x 161443x 2162986x 1244239x 1244239x 1244239x 1244239x | /**
* The draws a resampled map is made of (§6.5): a seeded generator, and the two
* discrete distributions that drawing contacts reduces to.
*
* **Subsampling** is choosing `N` of a file's `T` contacts, every subset equally
* likely — cooltools' `random-sample --count N --exact`, which draws `N` distinct
* indices out of `range(T)` and counts how many land in each pixel. The count a
* group of contacts receives is then **hypergeometric**, and a sample drawn group
* by group — sequence pairs, then pixels — with each group's share hypergeometric
* in what is left is the same distribution, drawn without ever holding `T`
* indices. **Oversampling** adds contacts drawn with replacement from a
* probability map, whose shares are **binomial** in the same way.
*
* Both samplers are numpy's, ported line for line (`random_binomial` with its
* BTPE and inversion arms, `random_hypergeometric` with HRUA and the direct
* draw), because they are the published, tested answers to how these two are
* drawn exactly at any size — a count of 190,000 on a diagonal pixel and a
* genome of 10⁹ contacts included — and a sampler written from a textbook is a
* sampler whose tails nobody has checked.
*
* The generator is **xoshiro128\*\***, not the analysis tasks' mulberry32
* (`analysis/random.ts`): that one's period is 2³², which is fine for resampling
* a few thousand rows and not for a stream that walks a hundred million contacts.
* Each group gets a generator of its own, seeded from the map's seed and the
* group's coordinates ({@link seedOf}), so what a pixel receives does not depend
* on which other pixels happened to be read beside it — the property that lets
* two workers, and the command line, arrive at the same map.
*/
/** 2⁵³, the denominator of a double drawn from 53 random bits. */
const TWO_53 = 9007199254740992
function rotl(x: number, k: number): number {
return (x << k) | (x >>> (32 - k))
}
/** A murmur3 finaliser: every input bit moves every output bit. */
function avalanche(h: number): number {
h = Math.imul(h ^ (h >>> 16), 0x85ebca6b)
h = Math.imul(h ^ (h >>> 13), 0xc2b2ae35)
return (h ^ (h >>> 16)) >>> 0
}
/**
* One 32-bit seed from several integers — a map's seed and a group's coordinates.
*
* Each part is folded in as two 32-bit halves, so a coordinate past 2³² (a pixel
* key on a large genome) is not truncated into a collision with a small one.
*/
export function seedOf(...parts: number[]): number {
let h = 0x811c9dc5
for (const part of parts) {
const value = Number.isFinite(part) ? Math.trunc(part) : 0
const lo = value >>> 0
const hi = Math.floor(value / 4294967296) >>> 0
h = avalanche((h ^ lo) + 0x9e3779b9)
h = avalanche((h ^ hi) + 0x7f4a7c15)
}
return h
}
/**
* xoshiro128\*\* (Blackman & Vigna), seeded through splitmix32.
*
* Four 32-bit words of state and a period of 2¹²⁸ − 1; every output is one
* `Math.imul` chain, which is what keeps a draw per contact affordable.
*/
export class Rng {
private s0: number
private s1: number
private s2: number
private s3: number
constructor(seed: number) {
let a = seed | 0
const splitmix = (): number => {
a = (a + 0x9e3779b9) | 0
let t = a ^ (a >>> 16)
t = Math.imul(t, 0x21f0aaad)
t ^= t >>> 15
t = Math.imul(t, 0x735a2d97)
return (t ^ (t >>> 15)) >>> 0
}
this.s0 = splitmix()
this.s1 = splitmix()
this.s2 = splitmix()
this.s3 = splitmix()
// The one state the generator cannot leave. Four zero words out of
// splitmix is not a thing that happens, and this makes it impossible.
if ((this.s0 | this.s1 | this.s2 | this.s3) === 0) this.s0 = 1
}
/** A uniform integer in `[0, 2³²)`. */
nextUint32(): number {
let s0 = this.s0
let s1 = this.s1
let s2 = this.s2
let s3 = this.s3
const result = Math.imul(rotl(Math.imul(s1, 5), 7), 9) >>> 0
const t = s1 << 9
s2 ^= s0
s3 ^= s1
s1 ^= s2
s0 ^= s3
s2 ^= t
s3 = rotl(s3, 11)
this.s0 = s0
this.s1 = s1
this.s2 = s2
this.s3 = s3
return result
}
/** A uniform double in `[0, 1)`, from 53 random bits. */
nextDouble(): number {
const hi = this.nextUint32() >>> 5
const lo = this.nextUint32() >>> 6
return (hi * 67108864 + lo) / TWO_53
}
/**
* A uniform integer in `[0, max]`, inclusive — numpy's `random_interval`.
*
* Masked rejection below 2³², which is exact; above it a 53-bit double is
* scaled, whose bias is one part in 2⁵³⁻ᵏ for a `max` of `k` bits — below
* anything a count could show.
*/
interval(max: number): number {
if (max <= 0) return 0
if (max <= 0xffffffff) {
let mask = max
mask |= mask >>> 1
mask |= mask >>> 2
mask |= mask >>> 4
mask |= mask >>> 8
mask |= mask >>> 16
for (;;) {
const value = (this.nextUint32() & mask) >>> 0
if (value <= max) return value
}
}
return Math.floor(this.nextDouble() * (max + 1))
}
}
// ---------------------------------------------------------------------------
// Binomial — numpy's `random_binomial`
// ---------------------------------------------------------------------------
/**
* The number of successes in `n` trials of probability `p`.
*
* numpy's dispatch: inversion where the mean is at most 30, BTPE above it, each
* on the smaller of `p` and `1 − p`.
*/
export function binomial(rng: Rng, n: number, p: number): number {
if (!(n > 0) || !(p > 0)) return 0
if (p >= 1) return n
if (p <= 0.5) {
return p * n <= 30 ? binomialInversion(rng, n, p) : binomialBtpe(rng, n, p)
}
const q = 1 - p
return n - (q * n <= 30 ? binomialInversion(rng, n, q) : binomialBtpe(rng, n, q))
}
/** `random_binomial_inversion`: walk the CDF from zero, restarting past the bound. */
function binomialInversion(rng: Rng, n: number, p: number): number {
const q = 1 - p
// `log1p` where numpy writes `log(q)`: the same number, without the digits
// `1 − p` loses when `p` is a few parts in a billion.
const qn = Math.exp(n * Math.log1p(-p))
const np = n * p
const bound = Math.min(n, Math.floor(np + 10 * Math.sqrt(np * q + 1)))
let x = 0
let px = qn
let u = rng.nextDouble()
while (u > px) {
x++
if (x > bound) {
x = 0
px = qn
u = rng.nextDouble()
} else {
u -= px
px = ((n - x + 1) * p * px) / (x * q)
}
}
return x
}
/** The Stirling correction BTPE's final acceptance test adds for one argument. */
function stirlingTail(v: number): number {
const v2 = v * v
return (13680 - (462 - (132 - (99 - 140 / v2) / v2) / v2) / v2) / v / 166320
}
/**
* `random_binomial_btpe` (Kachitvichyanukul & Schmeiser, 1988), for `p ≤ ½`.
*
* The four regions of the dominating function — triangle, two parallelograms,
* two exponential tails — and numpy's step labels in the comments, so the port
* can be read against its source.
*/
function binomialBtpe(rng: Rng, n: number, p: number): number {
const r = Math.min(p, 1 - p)
const q = 1 - r
const fm = n * r + r
const m = Math.floor(fm)
const p1 = Math.floor(2.195 * Math.sqrt(n * r * q) - 4.6 * q) + 0.5
const xm = m + 0.5
const xl = xm - p1
const xr = xm + p1
const c = 0.134 + 20.5 / (15.3 + m)
let a = (fm - xl) / (fm - xl * r)
const laml = a * (1 + a / 2)
a = (xr - fm) / (xr * q)
const lamr = a * (1 + a / 2)
const p2 = p1 * (1 + 2 * c)
const p3 = p2 + c / laml
const p4 = p3 + c / lamr
const nrq = n * r * q
const done = (y: number): number => (p > 0.5 ? n - y : y)
for (;;) {
// Step 10: the triangle.
const u = rng.nextDouble() * p4
let v = rng.nextDouble()
let y: number
if (u <= p1) return done(Math.floor(xm - p1 * v + u))
if (u <= p2) {
// Step 20: the parallelograms.
const x = xl + (u - p1) / c
v = v * c + 1 - Math.abs(m - x + 0.5) / p1
if (v > 1) continue
y = Math.floor(x)
} else if (u <= p3) {
// Step 30: the left tail.
y = Math.floor(xl + Math.log(v) / laml)
if (y < 0 || v === 0) continue
v = v * (u - p2) * laml
} else {
// Step 40: the right tail.
y = Math.floor(xr - Math.log(v) / lamr)
if (y > n || v === 0) continue
v = v * (u - p3) * lamr
}
// Step 50: an explicit ratio of probabilities, near the mode.
const k = Math.abs(y - m)
if (!(k > 20 && k < nrq / 2 - 1)) {
const s = r / q
const aa = s * (n + 1)
let f = 1
if (m < y) {
for (let i = m + 1; i <= y; i++) f *= aa / i - s
} else if (m > y) {
for (let i = y + 1; i <= m; i++) f /= aa / i - s
}
if (v > f) continue
return done(y)
}
// Step 52: squeeze on the log, then the Stirling bound.
const rho = (k / nrq) * ((k * (k / 3 + 0.625) + 0.16666666666666666) / nrq + 0.5)
const t = (-k * k) / (2 * nrq)
const A = Math.log(v)
if (A < t - rho) return done(y)
if (A > t + rho) continue
const x1 = y + 1
const f1 = m + 1
const z = n + 1 - m
const w = n - y + 1
const bound =
xm * Math.log(f1 / x1) +
(n - m + 0.5) * Math.log(z / w) +
(y - m) * Math.log((w * r) / (x1 * q)) +
stirlingTail(f1) +
stirlingTail(z) +
stirlingTail(x1) +
stirlingTail(w)
if (A > bound) continue
return done(y)
}
}
// ---------------------------------------------------------------------------
// Hypergeometric — numpy's `random_hypergeometric`
// ---------------------------------------------------------------------------
/** `log k!` below 126, where numpy tabulates it. */
const LOG_FACTORIAL = (() => {
const table = new Float64Array(126)
for (let k = 2; k < 126; k++) table[k] = table[k - 1] + Math.log(k)
return table
})()
const HALF_LN_2PI = 0.9189385332046728
/** `log k!`: the table, then Stirling's series to the `1/k³` term, as numpy has it. */
export function logFactorial(k: number): number {
if (k < 126) return LOG_FACTORIAL[k]
return (k + 0.5) * Math.log(k) - k + (HALF_LN_2PI + (1 / k) * (1 / 12 - 1 / (360 * k * k)))
}
/**
* How many of `sample` draws without replacement, from `good` marked and `bad`
* unmarked items, are marked.
*
* numpy's dispatch: the ratio-of-uniforms method once ten or more draws are
* taken and ten or more are left, the direct draw otherwise.
*/
export function hypergeometric(rng: Rng, good: number, bad: number, sample: number): number {
if (!(sample > 0) || !(good > 0)) return 0
if (!(bad > 0)) return Math.min(sample, good)
const total = good + bad
if (sample >= total) return good
if (sample >= 10 && sample <= total - 10) return hypergeometricHrua(rng, good, bad, sample)
return hypergeometricDirect(rng, good, bad, sample)
}
/** `hypergeometric_sample`: take the draws one at a time, from the smaller side. */
function hypergeometricDirect(rng: Rng, good: number, bad: number, sample: number): number {
const total = good + bad
const half = Math.floor(total / 2)
let left = sample > half ? total - sample : sample
let remainingTotal = total
let remainingGood = good
while (left > 0 && remainingGood > 0 && remainingTotal > remainingGood) {
remainingTotal--
if (rng.interval(remainingTotal) < remainingGood) remainingGood--
left--
}
// Only marked items remain: every draw left takes one.
if (remainingTotal === remainingGood) remainingGood -= left
return sample > half ? remainingGood : good - remainingGood
}
const HRUA_D1 = 1.7155277699214135
const HRUA_D2 = 0.8989161620588988
/**
* `hypergeometric_hrua` (Stadlober, 1989): ratio of uniforms under a table
* mountain, on the smaller of the two sides and the smaller of sample and
* complement, mapped back at the end.
*/
function hypergeometricHrua(rng: Rng, good: number, bad: number, sample: number): number {
const popsize = good + bad
const computedSample = Math.min(sample, popsize - sample)
const minGoodBad = Math.min(good, bad)
const maxGoodBad = Math.max(good, bad)
const p = minGoodBad / popsize
const q = maxGoodBad / popsize
const mu = computedSample * p
const a = mu + 0.5
const variance = ((popsize - computedSample) * computedSample * p * q) / (popsize - 1)
const c = Math.sqrt(variance + 0.5)
const h = HRUA_D1 * c + HRUA_D2
const m = Math.floor(((computedSample + 1) * (minGoodBad + 1)) / (popsize + 2))
const g =
logFactorial(m) +
logFactorial(minGoodBad - m) +
logFactorial(computedSample - m) +
logFactorial(maxGoodBad - computedSample + m)
// Sixteen standard deviations past the mean, or the end of the support.
const b = Math.min(Math.min(computedSample, minGoodBad) + 1, Math.floor(a + 16 * c))
let k: number
for (;;) {
const u = rng.nextDouble()
const v = rng.nextDouble()
// numpy divides by `U` unguarded; a zero here is one draw in 2⁵³, and it
// would turn `X` into a NaN that every comparison below lets through.
if (u === 0) continue
const x = a + (h * (v - 0.5)) / u
if (x < 0 || x >= b) continue
k = Math.floor(x)
const gp =
logFactorial(k) +
logFactorial(minGoodBad - k) +
logFactorial(computedSample - k) +
logFactorial(maxGoodBad - computedSample + k)
const t = g - gp
if (u * (4 - u) - 3 <= t) break
if (u * (u - t) >= 1) continue
if (2 * Math.log(u) <= t) break
}
if (good > bad) k = computedSample - k
if (computedSample < sample) k = good - k
return k
}
|