lkb-pipeline/core/lkb.R

source("boot.R")
box_use(./laurent2)

pair_idx <- function(a, b, pairs) {
  if (a > b) { tmp <- a; a <- b; b <- tmp }
  for (i in seq_along(pairs)) {
    if (pairs[[i]][1] == a && pairs[[i]][2] == b) return(i)
  }
  stop(sprintf("Pair (%d,%d) not found", a, b))
}

pairs_list <- function(n) {
  result <- list()
  for (j in 0:(n-2)) for (k in (j+1):(n-1)) result[[length(result) + 1]] <- c(j, k)
  result
}

lkb_generator_matrix <- function(k, n = 8, positive = TRUE) {
  pairs <- pairs_list(n)
  N <- length(pairs)
  m <- laurent2$l2_mat_zero(N)
  k0 <- k - 1L
  q <- laurent2$l2_q
  q2 <- laurent2$l2_mul(q, q)
  qm1 <- laurent2$l2_qinv
  qm2 <- laurent2$l2_mul(qm1, qm1)
  t <- laurent2$l2_t
  one <- laurent2$l2_one
  zero <- laurent2$l2_zero
  qsq_minus_q <- laurent2$laurent2(c("2,0" = 1, "1,0" = -1))
  one_minus_q <- laurent2$laurent2(c("0,0" = 1, "1,0" = -1))
  one_minus_qinv <- laurent2$laurent2(c("0,0" = 1, "-1,0" = -1))
  qm1_minus_qm2 <- laurent2$laurent2(c("-2,0" = 1, "-1,0" = -1))
  tqm1_minus_tqm2 <- laurent2$laurent2(c("-1,-1" = 1, "-2,-1" = -1))
  neg_q2_t <- laurent2$laurent2(c("2,1" = -1))
  neg_tinv_qinv2 <- laurent2$laurent2(c("-2,-1" = -1))
  neg_t_q2 <- laurent2$laurent2(c("2,1" = -1))

  for (idx in seq_len(N)) {
    j <- pairs[[idx]][1]
    kk <- pairs[[idx]][2]

    if (positive) {
      if (k0 == j - 1) {
        m[[pair_idx(k0, kk, pairs), idx]] <- q
        m[[pair_idx(k0, j, pairs), idx]] <- qsq_minus_q
        m[[pair_idx(j, kk, pairs), idx]] <- one_minus_q
      } else if (k0 == j && j != kk - 1) {
        m[[pair_idx(j + 1L, kk, pairs), idx]] <- one
        m[[idx, idx]] <- zero
      } else if (kk - 1 == k0 && kk - 1 != j) {
        m[[pair_idx(j, k0, pairs), idx]] <- q
        m[[pair_idx(j, kk, pairs), idx]] <- one_minus_q
        m[[pair_idx(k0, kk, pairs), idx]] <- laurent2$l2_mul(one_minus_q, laurent2$l2_mul(q, t))
      } else if (k0 == kk) {
        m[[pair_idx(j, kk + 1L, pairs), idx]] <- one
        m[[idx, idx]] <- zero
      } else if (k0 == j && j == kk - 1) {
        # Sage: A[idx, m] = -t*q*q = -q^2*t
        m[[idx, idx]] <- neg_q2_t
      } else {
        m[[idx, idx]] <- one
      }
    } else {
      if (k0 == j - 1) {
        m[[pair_idx(j - 1L, kk, pairs), idx]] <- one
      } else if (k0 == j && j != kk - 1) {
        m[[pair_idx(j + 1L, kk, pairs), idx]] <- qm1
        m[[pair_idx(j, kk, pairs), idx]] <- one_minus_qinv
        m[[pair_idx(j, j + 1L, pairs), idx]] <- tqm1_minus_tqm2
      } else if (kk - 1 == k0 && kk - 1 != j) {
        m[[pair_idx(j, kk - 1L, pairs), idx]] <- one
      } else if (k0 == kk) {
        m[[pair_idx(j, kk + 1L, pairs), idx]] <- qm1
        m[[pair_idx(j, kk, pairs), idx]] <- one_minus_qinv
        m[[pair_idx(kk, kk + 1L, pairs), idx]] <- qm1_minus_qm2
      } else if (k0 == j && j == kk - 1) {
        # Sage: A[idx, m] = -t^{-1}*q^{-2}
        m[[idx, idx]] <- neg_tinv_qinv2
      } else {
        m[[idx, idx]] <- one
      }
    }
  }
  m
}

#' @export
lkb_pair_to_idx <- function(i, j, n) {
  if (i > j) { tmp <- i; i <- j; j <- tmp }
  as.integer((i - 1) * (2 * n - i) / 2 + (j - i))
}

#' @export
lkb_idx_to_pair <- function(idx, n) {
  i <- 1
  while (idx > (n - i)) { idx <- idx - (n - i); i <- i + 1 }
  c(i, i + idx)
}

#' @export
lkb_generator <- function(k, n = 8) {
  lkb_generator_matrix(k, n, TRUE)
}

#' @export
lkb_generator_inverse <- function(k, n = 8) {
  lkb_generator_matrix(k, n, FALSE)
}

#' @export
lkb_generators <- function(n = 8) {
  lapply(1:(n - 1), function(k) lkb_generator_matrix(k, n, TRUE))
}

#' @export
lkb_generators_inverse <- function(n = 8) {
  lapply(1:(n - 1), function(k) lkb_generator_matrix(k, n, FALSE))
}

#' @export
lkb_matrix <- function(word, n = 8, gens = NULL, gens_inv = NULL) {
  if (is.null(gens)) gens <- lkb_generators(n)
  if (is.null(gens_inv)) gens_inv <- lkb_generators_inverse(n)
  N <- n * (n - 1) / 2
  result <- laurent2$l2_mat_identity(N)
  for (s in unclass(word)) {
    g <- abs(s)
    if (g < 1 || g > n - 1) stop(sprintf("Invalid generator %d for B%d", g, n))
    if (s > 0) {
      result <- laurent2$l2_mat_mul(result, gens[[g]])
    } else {
      result <- laurent2$l2_mat_mul(result, gens_inv[[g]])
    }
  }
  result
}

#' @export
lkb_basis_label <- function(idx, n) {
  p <- lkb_idx_to_pair(idx, n)
  paste0("e_{", p[1], ",", p[2], "}")
}

#' @export
lkb_matrix_str <- function(m, n = 8, threshold = 8) {
  N <- nrow(m)
  lines <- character(0)
  for (i in 1:min(N, threshold)) {
    row_parts <- character(min(N, threshold))
    for (j in 1:min(N, threshold)) {
      s <- laurent2$l2_str(m[[i, j]])
      row_parts[j] <- if (nchar(s) > 12) paste0(substr(s, 1, 10), "..") else s
    }
    lines <- c(lines, paste0(lkb_basis_label(i, n), ": ", paste(row_parts, collapse = " | ")))
  }
  if (N > threshold) lines <- c(lines, paste0("... (", N, "x", N, " matrix)"))
  paste(lines, collapse = "\n")
}