-
Notifications
You must be signed in to change notification settings - Fork 0
/
GC.IDX.R
111 lines (106 loc) · 3.78 KB
/
GC.IDX.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
GC.IDX <- function(x, cmax, cmin = 2, indexlist = "all", method = 'FCM', fzm = 2,
iter = 100, nstart = 20){
# Defined vector
IDX = c("gc1","gc2","gc3","gc4")
invisible(lapply(IDX, function(x) assign(x, vector(), envir = .GlobalEnv)))
distance =dist(x,diag = TRUE,upper= TRUE)
dd = as.matrix(distance)
dd[lower.tri(dd)] = 0
diag(dd) = 0
# For checking IDX
Check_GC1 = sum(indexlist %in% c("all","GC1"))>=1
Check_GC2 = sum(indexlist %in% c("all","GC2"))>=1
Check_GC3 = sum(indexlist %in% c("all","GC3"))>=1
Check_GC4 = sum(indexlist %in% c("all","GC4"))>=1
# Start k loop
for(k in cmin:cmax){
if(method == "EM"){ # EM Algorithm
EM.model <- Mclust(x,G=k,verbose=FALSE)
assign("m",EM.model$z)
assign("c",t(EM.model$parameters$mean))
}else if(method == "FCM"){ # FCM Algorithm
wd = Inf
# cm.out = list()
for (nr in 1:nstart){
FCM.model = cmeans(x,k,iter,verbose=FALSE,method="cmeans",m=fzm)
if (FCM.model$withinerror < wd){
wd = FCM.model$withinerror
FCM.model2 =FCM.model
}
}
assign("m",FCM.model2$membership)
assign("c",FCM.model2$centers)
}
# Defined variable
n = nrow(x)
mt = t(m)
NI = colSums(m)
nws = floor(sum(NI*(NI-1)/2))
# GC1: sum product
if(Check_GC1){
PD1 = m%*%mt
G1 = sum(PD1*dd)
Gmax1 = sum(sort(dd,decreasing = T)[1:nws]*sort(PD1[lower.tri(PD1)],decreasing = T)[1:nws])
Gmin1 = sum(sort(dd,decreasing = T)[1:nws]*sort(PD1[lower.tri(PD1)],decreasing = F)[1:nws])
gc1[k-cmin+1] = (G1-Gmin1)/(Gmax1-Gmin1)
}
if(sum(indexlist %in% c("all","GC2","GC3","GC4"))>=1){
PD2 = matrix(rep(0,n^2),n,n)
PD3 = matrix(rep(0,n^2),n,n)
PD4 = matrix(rep(0,n^2),n,n)
for (s in 1:n){
#GC2: sum-min
if(Check_GC2){
PD2[s,] = colSums(pmin(mt,m[s,]))
}
# GC3: max-product
if(Check_GC3){
PD3[s,] = apply(m[s,]*mt,2,max)
}
# GC4: max-min
if(Check_GC4){
PD4[s,] = apply(pmin(mt,m[s,]),2,max)
}
} # end loop s
if(Check_GC2){
G2 = sum(PD2*dd)
Gmax2 = sum(sort(dd,decreasing = T)[1:nws]*sort(PD2[lower.tri(PD2)],decreasing = T)[1:nws])
Gmin2 = sum(sort(dd,decreasing = T)[1:nws]*sort(PD2[lower.tri(PD2)],decreasing = F)[1:nws])
gc2[k-cmin+1] = (G2-Gmin2)/(Gmax2-Gmin2)
}
if(Check_GC3){
G3 = sum(PD3*dd)
Gmax3 = sum(sort(dd,decreasing = T)[1:nws]*sort(PD3[lower.tri(PD3)],decreasing = T)[1:nws])
Gmin3 = sum(sort(dd,decreasing = T)[1:nws]*sort(PD3[lower.tri(PD3)],decreasing = F)[1:nws])
gc3[k-cmin+1] = (G3-Gmin3)/(Gmax3-Gmin3)
}
if(Check_GC4){
G4 = sum(PD4*dd)
Gmax4 = sum(sort(dd,decreasing = T)[1:nws]*sort(PD4[lower.tri(PD4)],decreasing = T)[1:nws])
Gmin4 = sum(sort(dd,decreasing = T)[1:nws]*sort(PD4[lower.tri(PD4)],decreasing = F)[1:nws])
gc4[k-cmin+1] = (G4-Gmin4)/(Gmax4-Gmin4)
}
}# END GC index
}#end k loop
# Create data frame for IDX result
IDX.value = c("gc1","gc2","gc3","gc4")
IDX.name = c("GC1","GC2","GC3","GC4")
IDX.val_data = lapply(seq_along(IDX.value), function(i){
if(is.logical(get(IDX.value[i]))){
value = 0
}else{
value = get(IDX.value[i])
}
IDX.data = data.frame(K = cmin:cmax, value = value)
colnames(IDX.data) = c("c",paste0(IDX.name[i]))
return(IDX.data)
})
names(IDX.val_data) = IDX.name
invisible(list2env(IDX.val_data, envir = .GlobalEnv))
GC.list = list("GC1"= GC1,"GC2"= GC2,"GC3"= GC3,"GC4"= GC4)
if (sum(indexlist == "all")==1){
return(GC.list)
} else {
return(GC.list[indexlist])
}
}