library(ggplot2)
library(dplyr,warn.conflicts=F, quietly=T)
library(plotly,warn.conflicts=F, quietly=T)source("kruskal_montone.r")
source("nonmetric_plot.r")
source("nonmetric_mds.r")Playing around with the interactive plot below
delta = as.matrix(read.table("munsell.dat"))
#nonmetric_mds = function(delta, k=2, config_plot = TRUE, Shepard_plot = TRUE)
myresult = nonmetric_mds(delta, k=3,config_plot = TRUE, Shepard_plot = TRUE , html=TRUE)
myresult$myresult## $Xhat
## [,1] [,2] [,3]
## [1,] -3.131514e+00 -2.3851464 1.3250389
## [2,] -3.360069e+00 0.2650006 0.7949932
## [3,] -2.187123e+00 -1.3250242 -1.5900157
## [4,] -2.717217e+00 2.3851153 -1.0600289
## [5,] 2.120092e+00 -2.9152077 0.5299940
## [6,] 8.930722e-06 1.8550952 -0.7950183
## [7,] 3.445393e+00 -1.8550995 -0.5300012
## [8,] 3.180296e+00 1.0600075 -0.2650072
## [9,] 2.650132e+00 2.9152592 1.5900451
##
## $dhat
## 1 2 3 4 5 6 7 8
## 1 0.000000 2.712280 3.242418 5.349354 5.337830 5.681613 6.854037 7.364528
## 2 2.712280 0.000000 3.097136 2.889507 6.341618 4.043101 7.250155 6.673232
## 3 3.242418 3.097136 0.000000 3.785105 5.057196 3.940650 5.755853 6.021064
## 4 5.349354 2.889507 3.785105 0.000000 7.349909 2.781091 7.499207 6.096608
## 5 5.337830 6.341618 5.057196 7.349909 0.000000 5.385741 2.000960 4.190274
## 6 5.681613 4.043101 3.940650 2.781091 5.385741 0.000000 5.070153 3.320739
## 7 6.854037 7.250155 5.755853 7.499207 2.000960 5.070153 0.000000 2.939106
## 8 7.364528 6.673232 6.021064 6.096608 4.190274 3.320739 2.939106 0.000000
## 9 7.848054 6.616533 7.175781 6.009358 5.949705 3.719628 5.280469 2.676612
## 9
## 1 7.848054
## 2 6.616533
## 3 7.175781
## 4 6.009358
## 5 5.949705
## 6 3.719628
## 7 5.280469
## 8 2.676612
## 9 0.000000
##
## $stress
## [1] 0.03545689htmltools::tagList(list(myresult$pplist$p1, myresult$pplist$p2,myresult$pplist$p3,myresult$pp))Shepard_mine## function (myresult, delta, html = FALSE)
## {
## require(ggplot2, warn.conflicts = F, quietly = T)
## require(plotly, warn.conflicts = F, quietly = T)
## dhatmatrix = as.matrix(myresult$dhat)
## dhat_arr = dhatmatrix[lower.tri(dhatmatrix)]
## dobs_arr = delta[lower.tri(delta)]
## data = data.frame(dhat_arr, dobs_arr)
## p = ggplot(data = data, aes(dhat_arr, dobs_arr)) + geom_point() +
## labs(x = "dhat", y = "Observations", title = "Shepard Plot")
## pp = ggplotly(p)
## if (!html) {
## show(pp)
## return(NULL)
## }
## else return(pp)
## }
## <bytecode: 0x55f5e7a50da0>configplot## function (myresult, html = FALSE)
## {
## require(ggplot2, warn.conflicts = F, quietly = T)
## require(plotly, warn.conflicts = F, quietly = T)
## X = data.frame(myresult$Xhat)
## k = dim(X)[2]
## X$label = paste("V", 1:dim(X)[1], sep = "")
## plist = list()
## for (i in 1:3) {
## if (i == 3 & k >= 3) {
## p1 = plot_ly(X, x = ~X1, y = ~X2, z = ~X3, text = ~label) %>%
## add_markers(size = 2) %>% add_text(textposition = "Top right") %>%
## layout(title = "3D configuration")
## if (!html)
## show(p1)
## else plist$p3 = p1
## }
## if (i == 2 & k >= 2) {
## p = ggplot(X, aes(X1, X2, label = label)) + geom_point() +
## geom_text(nudge_x = 0.2) + labs(title = "2D Configuration Plot",
## x = "Dimension 1", y = "Dimension 2")
## p2 = ggplotly(p)
## if (!html)
## show(p2)
## else plist$p2 = p2
## }
## if (i == 1) {
## tempdf = data.frame(x = rep(0, dim(X)[1]), y = X$X1,
## label = X$label)
## p = ggplot(data = tempdf, aes(x = x, y = y, label = label)) +
## geom_point() + geom_text(nudge_x = 0.04) + scale_x_continuous(name = "",
## breaks = c(0), expand = c(-1, 1)) + labs(title = "1D Configuration Plot",
## y = "Dimension 1")
## p3 = ggplotly(p)
## if (!html)
## show(p3)
## else plist$p1 = p3
## }
## }
## if (html)
## return(plist)
## else return(NULL)
## }
## <bytecode: 0x55f5e4d68d80>#R CMD SHLIB jbkPava.c
dyn.load("jbkPava.so")
jbkPava <- function (x, w = rep(1, length(x))) {
h <-
.C("jbkPava",
x = as.double(x),
w = as.double (w),
n = as.integer(length(x)))
return (h$x)
}
objf = function(X_org){
delta = as.matrix(dist(X_org))
delta_arr = delta[lower.tri(delta)]
dhat = jbkPava(delta_arr)#jbkPava(sort(delta_arr))
return(sqrt(sum((dhat-delta_arr)^2)/sum(delta_arr)^2 ))
}
X_org = cmdscale(delta)
myresult = optim(X_org, objf, NULL, method = "L-BFGS-B")
myresult## $par
## [,1] [,2]
## [1,] -3.4671589 -2.1083321
## [2,] -3.8364419 0.2342597
## [3,] -1.7361587 -0.7027776
## [4,] -2.6731948 2.1083332
## [5,] 1.6398143 -3.0453694
## [6,] -0.2342604 1.1712944
## [7,] 3.9824021 -1.1901549
## [8,] 3.5138859 0.6839178
## [9,] 2.5768524 3.0830889
##
## $value
## [1] 0.05189119
##
## $counts
## function gradient
## 5 5
##
## $convergence
## [1] 0
##
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"#use mass to do isomass
MASS::isoMDS(as.dist(delta),k=3)## initial value 2.122891
## iter 5 value 0.668595
## iter 10 value 0.450975
## iter 15 value 0.300896
## iter 20 value 0.228028
## iter 25 value 0.148885
## iter 30 value 0.127102
## iter 35 value 0.125328
## iter 40 value 0.123509
## final value 0.121339
## converged## $points
## [,1] [,2] [,3]
## V1 -3.4433930 -1.7216619 1.15357470
## V2 -3.8219055 0.2129022 0.38293357
## V3 -1.5583455 -1.1002060 -1.13000036
## V4 -2.7685513 1.6440282 -0.77405673
## V5 1.5675403 -2.6395698 -0.20721640
## V6 0.2399494 1.5923721 -0.92661027
## V7 4.3951461 -1.5942707 0.16326221
## V8 3.4830533 0.3350341 0.04515819
## V9 1.9065063 3.2713718 1.29295509
##
## $stress
## [1] 0.1213394kruskal_montone## function (data, debugnow = FALSE)
## {
## delta_arr = data
## n = length(delta_arr)
## block = data.frame(value = delta_arr, size_perb = rep(1,
## n), previousi = (1:n) - 1, nexti = (1:n) + 1)
## active = 1
## upsatisfied = FALSE
## downsatisfied = FALSE
## while (TRUE) {
## if (active > n)
## break
## if (debugnow)
## print("*******************start big loop *******************")
## if (debugnow)
## print("*****************************************************")
## if (debugnow)
## print(paste("Active: ", active))
## next_curr = block$nexti[active]
## previous_curr = block$previousi[active]
## if (active < n & next_curr <= n)
## upactive = TRUE
## else (upactive = FALSE)
## acloop = 1
## while (TRUE) {
## next_curr = block$nexti[active]
## previous_curr = block$previousi[active]
## acloop = acloop + 1
## if (debugnow)
## print("**************true loop start*****************")
## if (debugnow)
## print(paste("active:", upactive))
## if (upactive) {
## if (next_curr == (n + 1))
## break
## if (debugnow)
## print("\n active: up")
## tempupsat = (block$value[next_curr]/block$size_perb[next_curr]) >=
## (block$value[active]/block$size_perb[active])
## if (tempupsat) {
## if (debugnow)
## print("\n upsatisfy: true")
## if (active != 1 & previous_curr >= 1)
## upactive = FALSE
## else upactive = TRUE
## }
## else {
## if (debugnow)
## print("\n upsatisfy: false")
## if (debugnow)
## print(paste(" downvalue", block$value[previous_curr],
## "curvalue", block$value[active], "upvalue",
## block$value[next_curr]))
## nextnext = block$nexti[next_curr]
## block$value[active] = (block$value[next_curr] +
## block$value[active])
## block$size_perb[active] = block$size_perb[active] +
## block$size_perb[next_curr]
## block$nexti[active] = nextnext
## if (debugnow)
## print(paste("next_curr", next_curr))
## if (nextnext < n + 1)
## block$previousi[nextnext] = active
## block$size_perb[next_curr] = 0
## if (debugnow)
## print(paste("* downvalue", block$value[previous_curr],
## "curvalue", block$value[active], "upvalue",
## block$value[next_curr]))
## if (active > 1 & previous_curr >= 1)
## upactive = FALSE
## else upactive = TRUE
## next_curr = block$nexti[active]
## previous_curr = block$previousi[active]
## }
## }
## else {
## if (previous_curr < 1)
## break
## if (debugnow)
## print("\n active: down")
## if (downsatisfied) {
## if (debugnow)
## print("\n downsatisfy: true")
## if (active < n & next_curr <= n)
## upactive = TRUE
## else upactive = FALSE
## }
## else {
## if (debugnow)
## print("\n downsatisfy: false")
## if (debugnow)
## print(paste(" downvalue", block$value[previous_curr],
## "curvalue", block$value[active], "upvalue",
## block$value[next_curr]))
## prevprev = block$previousi[previous_curr]
## if (debugnow)
## print(paste("**prev", previous_curr))
## block$value[active] = (block$value[previous_curr] +
## block$value[active])
## block$size_perb[active] = block$size_perb[active] +
## block$size_perb[previous_curr]
## block$previousi[active] = prevprev
## if (prevprev > 0)
## block$nexti[prevprev] = nextnext
## block$size_perb[previous_curr] = 0
## if (debugnow)
## print(paste("* downvalue", block$value[previous_curr],
## "curvalue", block$value[active], "upvalue",
## block$value[next_curr]))
## if (active < n & next_curr <= n) {
## upactive = TRUE
## }
## else {
## upactive = FALSE
## }
## next_curr = block$nexti[active]
## previous_curr = block$previousi[active]
## }
## }
## if (active == 1 | previous_curr < 1)
## downsatisfied = TRUE
## else {
## if ((block$value[previous_curr]/block$size_perb[previous_curr]) <=
## (block$value[active]/block$size_perb[active]))
## downsatisfied = TRUE
## else downsatisfied = FALSE
## }
## if (active == n | next_curr > n)
## upsatisfied = TRUE
## else {
## if ((block$value[next_curr]/block$size_perb[next_curr]) >=
## (block$value[active]/block$size_perb[active])) {
## upsatisfied = TRUE
## }
## else upsatisfied = FALSE
## }
## if (debugnow)
## print(paste("preavg", (block$value[previous_curr]/block$size_perb[previous_curr]),
## "curavg", (block$value[active]/block$size_perb[active]),
## "nextavg", (block$value[next_curr]/block$size_perb[next_curr])))
## if (upsatisfied & downsatisfied)
## break
## }
## if (active == n)
## break
## else (active = block$nexti[active])
## }
## block = cbind(block, index = 1:n)
## sizeperb_list = block$size_perb[block$size_perb != 0]
## stored_place = block$index[block$size_perb != 0]
## prev_acc = 0
## istoreindex = 1
## store_accum = array()
## for (istore in stored_place) {
## store_accum[istoreindex] = block$size_perb[istore] +
## prev_acc
## prev_acc = store_accum[istoreindex]
## istoreindex = istoreindex + 1
## }
## store_accum
## x = array()
## for (i in 1:length(store_accum)) {
## if (i == 1)
## x[1:store_accum[i]] = block$value[stored_place[i]]/block$size_perb[stored_place[i]]
## else x[(store_accum[i - 1] + 1):store_accum[i]] = block$value[stored_place[i]]/block$size_perb[stored_place[i]]
## }
## return(x)
## }
## <bytecode: 0x55f5e4c36060>