Changeset 15927 for branches/2886_SymRegGrammarEnumeration/ExpressionClustering_R/ClusteringScript.R

Timestamp:

05/03/18 19:39:55 (6 years ago)

Author:

gkronber

Message:

#2886 made some changes to the R script for clustering and visualzation for the book chapter

File:

• branches/2886_SymRegGrammarEnumeration/ExpressionClustering_R/ClusteringScript.R (modified) (2 diffs)

branches/2886_SymRegGrammarEnumeration/ExpressionClustering_R/ClusteringScript.R

@@ -3 +3 @@
 library(dplyr)
 
+setwd("D:/heal/documents/trunk/Publications/2018/GPTP/data");
+
+sentenceFileName <- "evaluations_allSentences_2018-04-13_16-40_TreeSize-7_1d.csv.gz";
+
+# read from CSV and store as R binary (must be done once to produce the .rds file)
+#evalData <- read.csv(sentenceFileName,header = TRUE, sep = ";", dec=",");
+#saveRDS(evalData, "evaluations_allSentences_2018-04-13_16-40_TreeSize-7_1d.csv.rds");
+
+# read from R binary (faster)
+evalData <- readRDS("evaluations_allSentences_2018-04-13_16-40_TreeSize-7_1d.csv.rds");
+
+max(evalData$R2.keijzer4);
+max(evalData$R2.keijzer9);
+max(evalData$R2.pagie);
+max(evalData$R2.nguyen5);
+max(evalData$R2.nguyen6);
+max(evalData$R2.nguyen7);
+
+outputs <- evalData[,10:109];
+
+#check zero mean, unit variance
+#mean(t(outputs[2,]))
+#sd(t(outputs[2,]))
+
+# check
+# plot(t(outputs[4,]))
+
+#apprNN <- randomProjectionTreeSearch(t(outputs), K=100, n_trees=50, distance_method="Euclidean", verbose=TRUE)
+# check ANN
+#cluster_1 <- tidyr::gather(dplyr::tbl_df(t(outputs)[,apprNN[,5]]), "rowNum", "value");
+#xs <- rep(seq(1:100),100)
+#ggplot(cluster_1, aes(x=xs, y=value, c=rowNum)) + geom_line();
+
+#edgeMatrix <- buildEdgeMatrix(t(outputs),apprNN, verbose=TRUE);
+#clusters <- hdbscan(edgeMatrix, apprNN, minPts = 10, K = 5, verbose=TRUE);
+
+# check cluster
+#cluster_1 <- tidyr::gather(dplyr::tbl_df(t(outputs)[,!is.na(clusters$clusters) & clusters$clusters==3]), "rowNum", "value");
+#xs <- rep(seq(1:100),nrow(cluster_1)/100) # reps must be the number of functions in the cluster
+#ggplot(cluster_1, aes(x=xs, y=value, c=rowNum)) + geom_line();
 
 
-sentenceFileName <- "D:/heal/documents/trunk/Publications/2018/GPTP/data/evaluations_allSentences_2018-04-13_16-40_TreeSize-7_1d.csv.gz";
-evalData <- read.csv(sentenceFileName,header = TRUE, sep = ";", dec=",");
-qualities <- evalData$R2.keijzer4;
-outputs <- evalData[,6:105];
-
-lv <- largeVis(outputs, dim=2, K = 50, n_trees = 50) # TODO scale?
-clusters <- hdbscan(lv, minPts = 3, K=50);
+lv <- largeVis(t(outputs), dim=2, distance_method="Cosine",
+               perplexity=100, K = 100, n_trees = 150, threads=4, 
+               save_neighbors = TRUE, save_edges = TRUE, verbose=TRUE) ;
+clusters <- hdbscan(lv, verbose=TRUE, threads=4, minPts = 10, K = 20);
 
 
+# calculate quality distribution for each cluster
+qualities <- evalData$R2.keijzer4;
+clusterQualities <- data.frame(Qualities = qualities, Clusters = clusters$clusters, x=t(lv$coords)[,1], y=t(lv$coords)[,2] );
+
+clusterQualityAvg <- clusterQualities %>% group_by(Clusters) %>% summarize(AvgQuality = mean(Qualities)) ;
+clusterQualityStdDev <- clusterQualities %>% group_by(Clusters) %>% summarize(StdDevQuality = sd(Qualities));
+clusterQualityCount <- clusterQualities %>% group_by(Clusters) %>% summarize(Count = n());
+clusterXCenter <- clusterQualities %>% group_by(Clusters) %>% summarize(meanX = mean(x));
+clusterYCenter <- clusterQualities %>% group_by(Clusters) %>% summarize(meanY = mean(y));
+clusterStats <- clusterQualityAvg %>% full_join(clusterQualityStdDev, by="Clusters") %>% full_join(clusterQualityCount, by="Clusters") %>% full_join(clusterXCenter, by ="Clusters") %>% full_join(clusterYCenter, by="Clusters");
+clusterStats <- dplyr::arrange(clusterStats, desc(AvgQuality));
+clusterStats$Rank <- seq(1:nrow(clusterStats));
+ggplot(clusterStats, aes(x = Rank, y=AvgQuality)) + geom_point();
+
+write.csv2(clusters$clusters, "cluster_assignment_new.csv", sep = " ", dec = ".");
+
+#check clusters
+for(i in seq(1:nrow(clusterStats))) {
+  clusterNumber <- clusterStats$Clusters[i] # number of cluster with smallest quality (error!)
+  cluster_i <- tidyr::gather(dplyr::tbl_df(t(outputs)[,!is.na(clusters$clusters) & clusters$clusters==clusterNumber]), "rowNum", "value");
+  xs <- rep(seq(1:100),nrow(cluster_i)/100) # reps must be the number of functions in the cluster
+  ggplot(cluster_i, aes(x=xs, y=value, c=rowNum)) +
+    theme_void() +
+    geom_line(alpha=0.1);
+  
+  ggsave(paste(as.character(i), as.character(round(clusterStats$meanX[i], 3)), as.character(round(clusterStats$meanY[i], 3)), ".png"));
+}
+
+funs_in_cluster <- t(outputs)[,!is.na(clusters$clusters) & clusters$clusters==1748]
+cor(method="pearson", target_keijzer4, t(outputs[20281, ]))
+plot(funs_in_cluster[,4], target_keijzer4)
+
+xi <- seq(0,9.99,0.1);
+#  x³  * exp(-x) * cos(x) * sin(x) * (sin(x)² * cos(x) - 1)
+target_keijzer4 <- xi^3 * exp(-xi) * cos(xi) * sin(xi) * (sin(xi)*sin(xi) * cos(xi) - 1);
+plot(xi, target_keijzer4);
 
 m <- data.frame(x=t(lv$coords)[,1], y=t(lv$coords)[,2], c=clusters$clusters, q=qualities, outputs)
-# plot mapped points
-ggplot(data=m, aes(x=x, y=y)) + geom_point(aes(color=c))  + theme(legend.position = "none") ;
+m_sub <- m[m$q<1.0,];
+
+# plot mapped points (clusters)
+ggplot(data=m, aes(x=x, y=y)) +
+  geom_point(aes(color=c))  +
+  theme(legend.position = "none")
+#  scale_color_gradient(low = "red",high = "black") 
+;
+ggsave("phenotypic_clusters.png")
+
+
+# plot mapped points (qualities)
+ggplot(data=m, aes(x=x, y=y)) +
+  geom_point(aes(color=q))  +
+  scale_color_gradientn(colors=heat.colors(30)) 
+;
+
+
+write.csv2(m, "mapping_evaluations_allSentences_2018-04-13_16-40_TreeSize-7_1d.csv.gz");
 
 cluster_n <- dplyr::filter(m, c==5);
@@ -25 +115 @@
 p <- ggplot(evals_cluster_n, aes(x=x, y=fx,color=f)) + geom_line();
 p
-