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