CLUSTER LAYOUTS

*1.3 Producing the space filling layouts *

# Cluster space filling layouts_________________________________

source("spaceFillingLayout.r")
library(MASS)
n = 50
m0 = rep(0, 8)
samp = mvrnorm(n, mu = m0, Sigma = diag(8))

for (i in 1:8) {
    m = m0  # All zeros
    m[i] = 8  # ith component 8
    samp = rbind(samp, mvrnorm(n, mu = m, diag(8)))
}

nineClust = hclust(dist(samp), method = "ward")
library(RColorBrewer)
colorCodes = brewer.pal(9, name = "Set1")
colorCodes[6] = "#000000"
nineClustLayout = spaceFillLayout(nineClust)
plot(nineClustLayout$points, type = "n", main = "Space Filling Layout")
colorIndex = (1:450 - 1)%/%50 + 1  # Many way to compute this
points(nineClustLayout$points, col = colorCodes[colorIndex], pch = 16)
spaceFillDrawBox(nineClustLayout, nclust = 9)

plot(nineClustLayout$points, type = "n", main = "Space Filling Layout")
spaceFillDrawBox(nineClustLayout, nclust = 80, pointFill = "#A0E0FF", pointBorder = "black")
points(nineClustLayout$points, col = colorCodes[colorIndex], pch = 16)

LJ_layout = function(position, radii, centricFactor = 1, separate = 4, slow = 10) {
    n = nrow(position)
    containRad = max(sqrt((position * position) %*% c(1, 1)) + radii)

    # central force on each child
    dxy = -position * cbind(radii, radii) * centricFactor/containRad

    # LJ force due to all pairs
    for (i in 1:(n - 1)) {
        for (j in (i + 1):n) {
            vec = position[i, ] - position[j, ]
            curDist = sqrt(sum(vec * vec))
            unitVec = vec/curDist
            desiredDist = radii[i] + radii[j] + separate  # variation 
            ratio = desiredDist/curDist
            expForce = (36 * ratio^12 - 6 * ratio^6)/curDist
            newDxy = expForce * unitVec
            dxy[i, ] = dxy[i, ] + newDxy
            dxy[j, ] = dxy[j, ] - newDxy
        }
    }
    # restrict maximum movement to 1/slow of radius
    dxyLength = sqrt((dxy * dxy) %*% c(1, 1))
    denom = ifelse(slow * dxyLength > radii, slow * dxyLength * radii, radii)
    dxy = dxy/cbind(denom, denom)
    position = position + dxy
    return(position)
}

# LJ_draw: Draws the clusters

LJ_draw = function(position, radii, mycolors = "#0080FF", dtime = 0.05) {
    z1 = Sys.time()
    z2 = Sys.time()
    while (difftime(z2, z1, units = "secs") < dtime) z2 = Sys.time()
    containRad = max(sqrt((position * position) %*% c(1, 1)) + radii)
    rb = c(-containRad, containRad)
    dev.hold(1)
    plot(rb, rb, type = "n")
    x = position[, 1]
    y = position[, 2]
    dev.hold(1)
    symbols(x, y, circles = radii, inches = FALSE, bg = mycolors, add = TRUE)
    dev.hold(1)
    symbols(0, 0, circles = containRad, inches = FALSE, add = TRUE)
    dev.hold(1)
    text(x, y, 1:length(x))
    dev.flush(4)
}

LJ_start = function(clusterSize) {
    radii = sqrt(clusterSize)
    nClust = length(radii)
    positionLength = cumsum(2 * radii) - radii - radii[1]
    positionLength[1] = 0
    angles = 2 * pi * runif(nClust)
    x = positionLength * cos(angles)
    y = positionLength * sin(angles)
    position = cbind(x, y)
    return(list(position = position, radii = radii))
}

# LJ_shiftCenters

LJ_shiftCenters = function(position, radii) {
    len = sqrt((position * position) %*% c(1, 1)) + radii
    containRad = max(len)
    ind = which(len == containRad)
    if (length(ind) > 1) 
        return(position)  # two or more touching circles
    touchPT = position[ind, ]  # circle center of touching circle 
    moveVec = -touchPT/sqrt(sum(touchPT^2))  # unit vector

    # Let the target position length be curRad-radii
    target = containRad - radii
    # Conceptually length(position + k*moveVec) = length(target)

    # Use quadratic formula
    quadC = position[, 1]^2 + position[, 2]^2 - target^2
    quadB = 2 * (position[, 1] * moveVec[1] + position[, 2] * moveVec[2])
    # quadA = sum(moveVec^2) = 1 since a is a unit vector

    tmp = sqrt(quadB^2 - 4 * quadC)
    kLarge = (-quadB + tmp)/2  # negative values in wrong direction
    # kSmall = (-quadB - tmp)/2 # for the ind circle this is 0
    k = min(kLarge)/2
    return(position + matrix(rep(k * moveVec, nrow(position)), ncol = ncol(position), 
        byrow = T))
}

# Example

clusterSize = c(60, 50, 40, 30, 20, 10)

tmp = LJ_start(clusterSize)
position = tmp$position
radii = tmp$radii
mycolors = rainbow(length(clusterSize))

par(pty = "s")
LJ_draw(position, radii, mycolors)

## [1] 0

centricFactor = 1
slow = 2  # Don't reduce the step size too much 
skip = 1
for (i in 1:400) {
    position = LJ_layout(position, radii, centricFactor, separate = 2, slow = slow)
    if (i%%skip == 0) 
        LJ_draw(position, radii, mycolors)
}

tmp = LJ_shiftCenters(position, radii)
mycolors = rainbow(length(clusterSize))
par(pty = "s")
LJ_draw(tmp, radii, mycolors)

## [1] 0