R version 3.5.0 alpha (2018-03-28 r74481) Copyright (C) 2018 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > #### Interface to intermediate level Fortran Routines > #### Purpose : explore different versions of dysta(), dysta3(), .. in > #### Fortran code. > > library(cluster) > > dysta <- function(x, kind = c("euclidean","manhattan", "SqEuclidean"), + dystaK = "dysta") + { + ## Purpose: + ## ------------------------------------------------------------------------- + ## Arguments: + ## ------------------------------------------------------------------------- + ## Author: Martin Maechler, Date: 3 Sep 2002, 08:21 + + kind <- match.arg(kind) + ndyst <- which(kind == eval(formals()$kind))# 1, 2, or 3 + n <- nrow(x <- as.matrix(x)) + p <- ncol(x) + storage.mode(x) <- "double" + hasNA <- apply(is.na(x), 2, any) # == apply(x, 2, anyNA) + if(any(hasNA)) { + ina <- is.na(x) + x[ina] <- valmd <- -1.1*max(abs(range(x, na.rm = TRUE))) + valmd <- rep(valmd, p) + } else valmd <- 0. + + dys <- double(1 + n*(n-1)/2) + jtmd <- as.integer(ifelse(hasNA, -1, 1)) + + r <- + if(dystaK == "dysta3") { + .C(cluster:::dysta3, + n, + jp = p, + x, + dys= dys, + ndyst= ndyst, + jtmd= jtmd, + valmd= valmd, + jhalt= integer(1))[c("dys", "jhalt")] + } else { + .Fortran(cluster:::dysta, + n, + jp = p, + x, + dys= dys, + ndyst= ndyst, + jtmd= jtmd, + valmd= valmd, + jhalt= integer(1))[c("dys", "jhalt")] + } + if(r$jhalt) { + cat("'jhalt' was ", r$jhalt, + " -- some dissimilarities will be missing.\n") + r$dys[r$dys == -1.] <- NA + } + r$dys + } > > (x <- cbind(c(0:6,NA), c(1,2,NA,7,NA,8:9,8))) [,1] [,2] [1,] 0 1 [2,] 1 2 [3,] 2 NA [4,] 3 7 [5,] 4 NA [6,] 5 8 [7,] 6 9 [8,] NA 8 > dysta(x) 'jhalt' was 1 -- some dissimilarities will be missing. [1] 0.000000 1.414214 2.828427 1.414214 6.708204 5.385165 1.414214 [8] 5.656854 4.242641 2.828427 1.414214 8.602325 7.211103 4.242641 [15] 2.236068 1.414214 10.000000 8.602325 5.656854 3.605551 2.828427 [22] 1.414214 9.899495 8.485281 NA 1.414214 NA 0.000000 [29] 1.414214 > (d1 <- dysta(x, kind = "m")) 'jhalt' was 1 -- some dissimilarities will be missing. [1] 0 2 4 2 9 7 2 8 6 4 2 12 10 6 3 2 14 12 8 5 4 2 14 12 NA [26] 2 NA 0 2 > (d3 <- dysta(x, kind = "m",dystaK = "dysta3")) 'jhalt' was 1 -- some dissimilarities will be missing. [1] 2 4 9 8 12 14 14 2 7 6 10 12 12 2 4 6 8 NA 2 3 5 2 2 4 NA [26] 2 0 2 0 > > identical(sort(d1), sort(d3)) # TRUE [1] TRUE > cbind(d1=d1[-1], d3=d3[-length(d3)], + dist=dist(x,"manhattan"), daisy= daisy(x,"manhattan")) d1 d3 dist daisy [1,] 2 2 2 2 [2,] 4 4 4 4 [3,] 2 9 9 9 [4,] 9 8 8 8 [5,] 7 12 12 12 [6,] 2 14 14 14 [7,] 8 14 14 14 [8,] 6 2 2 2 [9,] 4 7 7 7 [10,] 2 6 6 6 [11,] 12 10 10 10 [12,] 10 12 12 12 [13,] 6 12 12 12 [14,] 3 2 2 2 [15,] 2 4 4 4 [16,] 14 6 6 6 [17,] 12 8 8 8 [18,] 8 NA NA NA [19,] 5 2 2 2 [20,] 4 3 3 3 [21,] 2 5 5 5 [22,] 14 2 2 2 [23,] 12 2 2 2 [24,] NA 4 4 4 [25,] 2 NA NA NA [26,] NA 2 2 2 [27,] 0 0 0 0 [28,] 2 2 2 2 > > identical(d3[-length(d3)], + c(dist(x,"manhattan")))# ! [1] TRUE > identical(c(daisy(x,"manhattan")), c(dist(x,"manhattan"))) [1] TRUE > identical(c(daisy(x,"euclidean")), c(dist(x,"euclidean"))) [1] TRUE > identical(dysta(x, dystaK="dysta3")[-length(d3)], + c(dist(x,"euclidean")))# ! 'jhalt' was 1 -- some dissimilarities will be missing. [1] TRUE > > proc.time() user system elapsed 0.171 0.038 0.255