1## manifold ranking 2## author: alexandros 3 4setGeneric("ranking",function(x, ...) standardGeneric("ranking")) 5setMethod("ranking",signature(x="matrix"), 6 function (x, 7 y, 8 kernel = "rbfdot", 9 kpar = list(sigma = 1), 10 scale = FALSE, 11 alpha = 0.99, 12 iterations = 600, 13 edgegraph = FALSE, 14 convergence = FALSE, 15 ...) 16 { 17 m <- dim(x)[1] 18 d <- dim(x)[2] 19 if(length(y) != m) 20 { 21 ym <- matrix(0,m,1) 22 ym[y] <- 1 23 y <- ym 24 } 25 if (is.null(y)) 26 y <- matrix(1, m, 1) 27 labelled <- y != 0 28 if (!any(labelled)) stop("no labels sublied") 29 30 if(is.character(kernel)) 31 kernel <- match.arg(kernel,c("rbfdot","polydot","tanhdot","vanilladot","besseldot","laplacedot")) 32 33 34 if(!is(kernel,"kernel")) 35 { 36 if(is(kernel,"function")) kernel <- deparse(substitute(kernel)) 37 kernel <- do.call(kernel, kpar) 38 } 39 40 if(!is(kernel,"kernel")) stop("kernel must inherit from class `kernel'") 41 42 if(scale) 43 x <- scale(x) 44 ## scaling from ksvm 45 ## normalize ? 46 47 if (is(kernel)[1]=='rbfkernel' && edgegraph){ 48 sigma = kpar(kernel)$sigma 49 n <- dim(x)[1] 50 dota <- rowSums(x*x)/2 51 sed <- crossprod(t(x)) 52 for (i in 1:n) 53 sed[i,] <- - 2*(sed[i,] - dota - rep(dota[i],n)) 54 diag(sed) <- 0 55 K <- exp(- sigma * sed) 56 57 mst <- minimum.spanning.tree(sed) 58 algo.mst <- mst$E 59 max.squared.edge.length <- mst$max.sed.in.tree 60 edgegraph <- (sed <= max.squared.edge.length) 61 K[!edgegraph] <- 0 62 ##algo.edge.graph <- sparse(algo.edge.graph) 63 rm(sed) 64 gc() 65 } 66 else { 67 if(edgegraph && is(kernel)[1]!="rbfkernel") 68 warning('edge graph is only implemented for use with the RBF kernel') 69 edgegraph <- matrix() 70 K <- kernelMatrix(kernel,x) 71 } 72 73 diag(K) <- 0 74 ##K <- sparse(K) 75 cs <- colSums(K) 76 ##cs[cs <= 10e-6] <- 1 77 78 D <- 1/sqrt(cs) 79 K <- D * K %*% diag(D) 80 81 if(sum(labelled)==1) 82 y <- K[, labelled,drop = FALSE] 83 else 84 y <- as.matrix(colSums(K[, labelled])) 85 K <- alpha * K[, !labelled] 86 ym <- matrix(0,m,iterations) 87 ym[,1] <- y 88 for (iteration in 2:iterations) 89 ym[, iteration] <- ym[, 1] + K %*% ym[!labelled, iteration-1] 90 91 ym[labelled,] <- NA 92 r <- ym 93 r[!labelled,] <- compute.ranks(-r[!labelled, ]) 94 if(convergence) 95 convergence <- (r - rep(r[,dim(r)[2]],iterations))/(m-sum(labelled)) 96 else 97 convergence <- matrix() 98 res <- cbind(t(t(1:m)), ym[,iterations], r[,iterations]) 99 return(new("ranking", .Data=res, convergence = convergence, edgegraph = edgegraph)) 100 }) 101 102 103## kernelMatrix interface 104setMethod("ranking",signature(x="kernelMatrix"), 105 function (x, 106 y, 107 alpha = 0.99, 108 iterations = 600, 109 convergence = FALSE, 110 ...) 111 { 112 m <- dim(x)[1] 113 114 if(length(y) != m) 115 { 116 ym <- matrix(0,m,1) 117 ym[y] <- 1 118 y <- ym 119 } 120 if (is.null(y)) 121 y <- matrix(1, m, 1) 122 labelled <- y != 0 123 if (!any(labelled)) stop("no labels sublied") 124 125 diag(x) <- 0 126 ##K <- sparse(K) 127 cs <- colSums(x) 128 ##cs[cs <= 10e-6] <- 1 129 130 D <- 1/sqrt(cs) 131 x <- D * x %*% diag(D) 132 133 if(sum(labelled)==1) 134 y <- x[, labelled,drop = FALSE] 135 else 136 y <- as.matrix(colSums(x[, labelled])) 137 x <- alpha * x[, !labelled] 138 ym <- matrix(0,m,iterations) 139 ym[,1] <- y 140 for (iteration in 2:iterations) 141 ym[, iteration] <- ym[, 1] + x %*% ym[!labelled, iteration-1] 142 143 ym[labelled,] <- NA 144 r <- ym 145 r[!labelled,] <- compute.ranks(-r[!labelled, ]) 146 if(convergence) 147 convergence <- (r - rep(r[,dim(r)[2]],iterations))/(m-sum(labelled)) 148 else 149 convergence <- matrix() 150 res <- cbind(t(t(1:m)), ym[,iterations], r[,iterations]) 151 return(new("ranking", .Data=res, convergence = convergence)) 152 }) 153 154 155## list interface 156setMethod("ranking",signature(x="list"), 157 function (x, 158 y, 159 kernel = "stringdot", 160 kpar = list(length = 4, lambda = 0.5), 161 alpha = 0.99, 162 iterations = 600, convergence = FALSE, ...) 163 { 164 m <- length(x) 165 166 if(length(y) != m) 167 { 168 ym <- matrix(0,m,1) 169 ym[y] <- 1 170 y <- ym 171 } 172 173 if (is.null(y)) 174 y <- matrix(1, m, 1) 175 labelled <- y != 0 176 if (!any(labelled)) stop("no labels sublied") 177 178 if(is.character(kernel)) 179 kernel <- match.arg(kernel,c("rbfdot","polydot","tanhdot","vanilladot","besseldot","laplacedot")) 180 181 if(!is(kernel,"kernel")) 182 { 183 if(is(kernel,"function")) kernel <- deparse(substitute(kernel)) 184 kernel <- do.call(kernel, kpar) 185 } 186 187 if(!is(kernel,"kernel")) stop("kernel must inherit from class `kernel'") 188 189 edgegraph <- matrix() 190 K <- kernelMatrix(kernel,x) 191 192 diag(K) <- 0 193 ##K <- sparse(K) 194 cs <- colSums(K) 195 ##cs[cs <= 10e-6] <- 1 196 197 D <- 1/sqrt(cs) 198 K <- D * K %*% diag(D) 199 200 if(sum(labelled)==1) 201 y <- K[, labelled,drop = FALSE] 202 else 203 y <- as.matrix(colSums(K[, labelled])) 204 K <- alpha * K[, !labelled] 205 ym <- matrix(0,m,iterations) 206 ym[,1] <- y 207 for (iteration in 2:iterations) 208 ym[, iteration] <- ym[, 1] + K %*% ym[!labelled, iteration-1] 209 210 ym[labelled,] <- NA 211 r <- ym 212 r[!labelled,] <- compute.ranks(-r[!labelled, ]) 213 if(convergence) 214 convergence <- (r - rep(r[,dim(r)[2]],iterations))/(m-sum(labelled)) 215 else 216 convergence <- matrix() 217 res <- cbind(t(t(1:m)), ym[,iterations], r[,iterations]) 218 return(new("ranking", .Data=res, convergence = convergence, edgegraph = NULL)) 219 }) 220 221minimum.spanning.tree <- function(sed) 222 { 223 max.sed.in.tree <- 0 224 E <- matrix(0,dim(sed)[1],dim(sed)[2]) 225 n <- dim(E)[1] 226 C <- logical(n) 227 cmp <- sed 228 diag(cmp) <- NA 229 ans <- min(cmp, na.rm = TRUE) 230 i <- which.min(cmp) 231 j <- i%/%n + 1 232 i <- i%%n +1 233 234 for (nC in 1:n) { 235 cmp <- sed 236 cmp[C,] <- NA 237 cmp[,!C] <- NA 238 if(nC == 1) 239 { 240 ans <- 1 241 i <- 1 242 } 243 else{ 244 ans <- min(cmp, na.rm=TRUE) 245 i <- which.min(cmp)} 246 j <- i%/%n + 1 247 i <- i%%n + 1 248 E[i, j] <- nC 249 E[j, i] <- nC 250 C[i] <- TRUE 251 max.sed.in.tree <- max(max.sed.in.tree, sed[i, j]) 252 } 253 ## E <- sparse(E) 254 res <- list(E=E, max.sed.in.tree=max.sed.in.tree) 255 } 256 257compute.ranks <- function(am) { 258 259 rm <- matrix(0,dim(am)[1],dim(am)[2]) 260 for (j in 1:dim(am)[2]) 261 { 262 a <- am[, j] 263 sort <- sort(a, index.return = TRUE) 264 sorted <- sort$x 265 r <- sort$ix 266 r[r] <- 1:length(r) 267 268 while(1) 269 { 270 if(sum(na.omit(diff(sorted) == 0)) == 0) 271 break 272 tied <- sorted[min(which(diff(sorted) == 0))] 273 sorted[sorted==tied] <- NA 274 r[a==tied] <- mean(r[a==tied]) 275 } 276 rm[, j] <- r 277 } 278 return(rm) 279} 280 281setMethod("show","ranking", 282 function(object) 283 { cat("Ranking object of class \"ranking\"","\n") 284 cat("\n") 285 show(object@.Data) 286 cat("\n") 287 if(!any(is.na(convergence(object)))) 288 cat("convergence matrix included.","\n") 289 if(!any(is.na(edgegraph(object)))) 290 cat("edgegraph matrix included.","\n") 291 }) 292