1## wrapper function for couplers
2## author : alexandros karatzoglou
3
4couple <- function(probin, coupler = "minpair")
5{
6  if(is.vector(probin))
7    probin <- matrix(probin,1)
8    m <- dim(probin)[1]
9
10  coupler <- match.arg(coupler, c("minpair", "pkpd", "vote", "ht"))
11
12#  if(coupler == "ht")
13#    multiprob <- sapply(1:m, function(x) do.call(coupler, list(probin[x ,], clscnt)))
14#  else
15    multiprob <- sapply(1:m, function(x) do.call(coupler, list(probin[x ,])))
16
17  return(t(multiprob))
18}
19
20
21ht <- function(probin, clscnt, iter=1000)
22  {
23    nclass <- length(clscnt)
24    probim <- matrix(0, nclass, nclass)
25    for(i in 1:nclass)
26      for(j in 1:nclass)
27        if(j>i)
28          {
29            probim[i,j] <- probin[i]
30            probim[j,i] <- 1 - probin[i]
31          }
32
33    p <- rep(1/nclass,nclass)
34    u <- matrix((1/nclass)/((1/nclass)+(1/nclass)) ,nclass,nclass)
35    iter <- 0
36
37    while(TRUE)
38      {
39        iter <- iter + 1
40        stoperror <- 0
41
42        for(i in 1:nclass){
43          num <- den <- 0
44          for(j in 1:nclass)
45            {
46              if (j!=i)
47                {
48                  num <- num + (clscnt[i] + clscnt[j]) * probim[i,j]
49                  den <- den + (clscnt[i] + clscnt[j]) * u[i,j]
50                }
51            }
52          alpha <- num/(den + 1e-308)
53          p[i] <- p[i]*alpha
54          stoperror <- stoperror + (alpha -1)^2
55          if(0)
56            {
57              sum <- 0
58              sum <- sum(p) + sum
59              p <- p/sum
60              for(ui in 1:nclass)
61                for(uj in 1:nclass)
62                  u[ui, uj] <- p[ui]/(p[ui] + p[uj])
63            }
64          else
65            {
66              for(j in 1:nclass)
67                if (i!=j)
68                  {
69                    u[i,j] <- p[i]/(p[i] + p[j])
70                    u[j,i] <- 1 - u[i,j]
71                  }
72            }
73        }
74        if(stoperror < 1e-3)
75          break
76        if(iter > 400)
77          {
78            cat("Too many iterations: aborting", probin, iter, stoperror, p)
79            break
80          }
81      }
82    ## normalize prob.
83    p <- p/sum(p)
84    return(p)
85  }
86
87
88minpair <- function(probin)
89  {  ## Count number of classes and construct prob. matrix
90    nclass <- (1+sqrt(1 + 8*length(probin)))/2
91    if(nclass%%1 != 0) stop("Vector has wrong length only one against one problems supported")
92    probim <- matrix(0, nclass, nclass)
93    probim[upper.tri(probim)] <- probin
94    probim[lower.tri(probim)] <- 1 - probin
95
96    sum <- colSums(probim^2)
97    Q <- diag(sum)
98    Q[upper.tri(Q)] <- - probin*(1 - probin)
99    Q[lower.tri(Q)] <- - probin*(1 - probin)
100    SQ <- matrix(0,nclass +1, nclass +1)
101    SQ[1:(nclass+1) <= nclass, 1:(nclass+1) <= nclass] <- Q
102    SQ[1:(nclass+1) > nclass, 1:(nclass+1) <= nclass] <- rep(1,nclass)
103    SQ[1:(nclass+1) <= nclass, 1:(nclass+1) > nclass] <- rep(1,nclass)
104
105    rhs <- rep(0,nclass+1)
106    rhs[nclass + 1] <- 1
107
108    p <- solve(SQ,rhs)
109
110    p <- p[-(nclass+1)]/sum(p[-(nclass+1)])
111    return(p)
112  }
113
114
115pkpd <- function(probin)
116  {  ## Count number of classes and constuct prob. matrix
117    nclass <- k <- (1+sqrt(1 + 8*length(probin)))/2
118    if(nclass%%1 != 0) stop("Vector has wrong length only one against one problems supported")
119    probim <- matrix(0, nclass, nclass)
120    probim[upper.tri(probim)] <- probin
121    probim[lower.tri(probim)] <- 1 - probin
122
123    probim[probim==0] <- 1e-300
124    R <- 1/probim
125    diag(R)  <-  0
126    p <- 1/(rowSums(R) - (k-2))
127
128    p <- p/sum(p)
129    return(p)
130  }
131
132
133vote<- function(probin)
134{
135  nclass <- (1+sqrt(1 + 8*length(probin)))/2
136  if(nclass%%1 != 0) stop("Vector has wrong length only one against one problems supported")
137
138  votev <- rep(0,nclass)
139  p <- 0
140  for(i in 1:(nclass-1))
141    {
142      jj <- i+1
143      for(j in jj:nclass)
144        {
145          p <- p+1
146          votev[i][probin[i] >= 0.5] <- votev[i][probin[i] >= 0.5] + 1
147          votev[j][probin[j] < 0.5] <- votev[j][probin[j] < 0.5] + 1
148        }
149    }
150
151  p <- votev/sum(votev)
152  return(p)
153}
154
155
156