R/kinship2_alignped4.R

In ggpedigree: Visualizing Pedigrees with 'ggplot2' and 'plotly'

# Automatically generated from all.nw using noweb

#' Compute optimal horizontal spacing for pedigree alignment
#'
#' This is an internal helper function for pedigree alignment. It uses quadratic
#' programming to find optimal horizontal positions for subjects that minimize
#' the distance between parents and children while keeping spouses together and
#' respecting spacing constraints. Requires the quadprog package.
#'
#' @param rval Aligned pedigree structure from previous alignment steps
#' @param spouse Logical matrix indicating spouse connections
#' @param level Integer vector of generation levels
#' @param width Numeric, maximum width of the pedigree plot
#' @param align Logical or numeric vector. If logical, uses default alignment parameters.
#'   If numeric, should be a vector c(a1, a2) where a1 controls parent-child penalties
#'   and a2 controls spouse penalties
#' @return Matrix of optimized horizontal positions for each subject
#' @keywords internal
kinship2_alignped4 <- function(rval, spouse, level, width, align) {
  ## Doc: alignped4 -part1, spacing across page
  if (is.logical(align)) align <- c(1.5, 2) # defaults
  maxlev <- nrow(rval$nid)
  width <- max(width, rval$n + .01) # width must be > the longest row

  n <- sum(rval$n) # total number of subjects
  myid <- matrix(0, maxlev, ncol(rval$nid)) # number the plotting points
  for (i in 1:maxlev) {
    myid[i, rval$nid[i, ] > 0] <- cumsum(c(0, rval$n))[i] + 1:rval$n[i]
  }

  # There will be one penalty for each spouse and one for each child
  npenal <- sum(spouse[rval$nid > 0]) + sum(rval$fam > 0)
  pmat <- matrix(0., nrow = npenal + 1, ncol = n)

  ## Doc: alignped4 -part2
  indx <- 0
  # Penalties to keep spouses close
  for (lev in 1:maxlev) {
    if (any(spouse[lev, ])) {
      who <- which(spouse[lev, ])
      indx <- max(indx) + 1:length(who)
      pmat[cbind(indx, myid[lev, who])] <- sqrt(align[2])
      pmat[cbind(indx, myid[lev, who + 1])] <- -sqrt(align[2])
    }
  }

  # Penalties to keep kids close to parents
  for (lev in (1:maxlev)[-1]) { # no parents at the top level
    families <- unique(rval$fam[lev, ])
    families <- families[families != 0] # 0 is the 'no parent' marker
    for (i in families) { # might be none
      who <- which(rval$fam[lev, ] == i)
      k <- length(who)
      indx <- max(indx) + 1:k # one penalty per child
      penalty <- sqrt(k^(-align[1]))
      pmat[cbind(indx, myid[lev, who])] <- -penalty
      pmat[cbind(indx, myid[lev - 1, rval$fam[lev, who]])] <- penalty / 2
      pmat[cbind(indx, myid[lev - 1, rval$fam[lev, who] + 1])] <- penalty / 2
    }
  }
  maxrow <- min(which(rval$n == max(rval$n)))
  pmat[nrow(pmat), myid[maxrow, 1]] <- 1e-5
  ncon <- n + maxlev # number of constraints
  cmat <- matrix(0., nrow = ncon, ncol = n)
  coff <- 0 # cumulative constraint lines so var
  dvec <- rep(1., ncon)
  for (lev in 1:maxlev) {
    nn <- rval$n[lev]
    if (nn > 1) {
      for (i in 1:(nn - 1)) {
        cmat[coff + i, myid[lev, i + 0:1]] <- c(-1, 1)
      }
    }

    cmat[coff + nn, myid[lev, 1]] <- 1 # first element >=0
    dvec[coff + nn] <- 0
    cmat[coff + nn + 1, myid[lev, nn]] <- -1 # last element <= width-1
    dvec[coff + nn + 1] <- 1 - width
    coff <- coff + nn + 1
  }

  if (requireNamespace("quadprog", quietly = TRUE)) {
    pp <- t(pmat) %*% pmat + 1e-8 * diag(ncol(pmat))
    fit <- quadprog::solve.QP(pp, rep(0., n), t(cmat), dvec)
  } else {
    stop("Need the quadprog package")
  }

  newpos <- rval$pos
  # fit <- lsei(pmat, rep(0, nrow(pmat)), G=cmat, H=dvec)
  # newpos[myid>0] <- fit$X[myid]
  newpos[myid > 0] <- fit$solution[myid]
  newpos
}