thomson <- function (b, uold = rep (0, nrow (b)), p = 2, itmax = 1000, eps = 1e-10, verbose = TRUE) { bold <- b - diag (uold) cold <- lrank (bold, p) fold <- sum ((bold - cold) ^ 2) eold <- Inf itel <- 1 repeat { unew <- diag (b - cold) enew <- mnorm (uold - unew) rnew <- enew / eold bnew <- b - diag (unew) cnew <- lrank (bnew, p) fnew <- sum ((bnew - cnew) ^ 2) if (verbose) { cat ( formatC (itel, width = 4, format = "d"), formatC ( fold, digits = 10, width = 15, format = "f" ), formatC ( fnew, digits = 10, width = 15, format = "f" ), formatC ( enew, digits = 10, width = 15, format = "f" ), formatC ( rnew, digits = 10, width = 15, format = "f" ), "\n" ) } if ((itel == itmax) || (fold - fnew) < eps) break itel <- itel + 1 uold <- unew fold <- fnew cold <- cnew bold <- bnew eold <- enew } return (list ( u = unew, f = fnew, e = enew, r = rnew, itel = itel )) } tmap <- function (b, u, p = 2) { c <- lrank (b - diag (u)) return (diag (b - c)) } tderivative <- function (b, u, p = 2) { n <- length (u) e <- eigen (b - diag (u)) k <- e\$vectors l <- e\$values m <- matrix (0, n, n) for (i in 1:n) for (j in 1:n) for (s in 1:p) { m[i,j] <- m[i,j] + (k[i,s] ^ 2) * (k[j, s] ^ 2) for (t in 1:n) { if (s == t) next xi <- l[s] / (l[t] - l[s]) m[i, j] <- m[i, j] - 2 * xi * k[i, s] * k[i, t] * k[j, s] * k[j, t] } } return (m) }