Skip to content

Instantly share code, notes, and snippets.

@stla
Last active March 31, 2024 11:50
Show Gist options
  • Save stla/ed755691ffcb6dec4846a72f18b51cba to your computer and use it in GitHub Desktop.
Save stla/ed755691ffcb6dec4846a72f18b51cba to your computer and use it in GitHub Desktop.
Kostka numbers with R
rg <- function(a, b) {
if(a <= b) a:b else integer(0L)
}
KostkaNumber <- function(lambda, mu) {
lambda <- as.integer(lambda)
mu <- as.integer(mu)
wmu <- sum(mu)
wlam <- sum(lambda)
if(wlam == 0L) {
return(as.integer(wmu == 0L))
}
if(wmu != wlam || !jack:::isDominated(mu, lambda)) {
return(0L)
}
nlam <- length(lambda)
nmu <- length(mu)
n <- max(nlam, nmu)
lambda <- c(lambda, rep(0L, n - nlam))
mu <- c(mu, rep(0L, n - nmu))
revlam <- rev(lambda)
boundedNonIncrSeqs <- function(h0, aas, bbs) {
if(length(aas) == 0L || length(bbs) == 0L) {
list(integer(0L))
} else {
a <- aas[1L]
b <- bbs[1L]
as <- aas[-1L]
bs <- bbs[-1L]
h_ <- rg(max(0L, a), min(h0, b))
do.call(c, lapply(h_, function(h) {
lapply(boundedNonIncrSeqs(h, as, bs), function(hs) {
c(h, hs)
})
}))
}
}
worker <- function(rlrls, smusmus, aacc, lastx0lastrowt) {
if(length(rlrls) <= 1L) {
1L
} else {
x0 <- smusmus[1L] - aacc[1L]
rl <- rlrls[1L]
rls <- rlrls[-1L]
smus <- smusmus[-1L]
acc <- aacc[-1L]
nacc <- length(acc)
lastx0 <- lastx0lastrowt[1L]
lastrowt <- lastx0lastrowt[-1L]
aas <- vapply(c(max(lastx0, x0), lastrowt), function(i) {
max(rl, i)
}, integer(1L))
rows <- boundedNonIncrSeqs(x0, aas, lambda)
sum(vapply(rows, function(row) {
l <- length(row) - 1L
trow <- tail(row, l)
irow <- head(row, l)
m <- min(nacc, l)
worker(rls, smus, head(acc, m) + head(trow, m), irow)
}, integer(1L)))
}
}
worker(revlam, cumsum(mu), rep(0L, n-1L), rep(0L, n))
}
KostkaNumber(c(4,1,1), c(2,1,1,1,1))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment