Changeset 16174 for branches/2886_SymRegGrammarEnumeration/ExpressionClustering_R/FindClustersForGPLog.R

Timestamp:

09/22/18 07:12:57 (6 years ago)

Author:

gkronber

Message:

#2886 made several changes to clustering scripts for GPTP paper

File:

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

branches/2886_SymRegGrammarEnumeration/ExpressionClustering_R/FindClustersForGPLog.R

@@ -1 +1 @@
 # find nearest neighbors for GP indiviuals
 
+library(dplyr)
+library(hexbin)
 library(ggplot2);
+#library(doSNOW); #dopar
+#library(foreach); # dopar
 
-gp_log <- read.csv2("C:/Users/P24581/filebox/GPTP 2018/symbreg-models-13.05.2018-1418.csv",header = FALSE,sep='\t',dec=',');
+#cl<-makeCluster(2) #change the 2 to your number of CPU cores
+#registerDoSNOW(cl)
+
+setwd("C:/Users/P24581/filebox/GPTP 2018/Keijzer4");
+
+gp_log <- read.csv2("symbreg-models-13.05.2018-1412.csv.gz",header = FALSE,sep='\t',dec=',');
 idx <- seq(1:nrow(gp_log));
 # check popSize in gp_log
-gp_log[seq(1,34000,500),1]
+gp_log[seq(1,20500,500),1]
 popSize <- 500;
 
@@ -13 +22 @@
 
 # generations <- seq(1,34000/popSize,1);
-generations <- seq(1,15,1);
+generations <- seq(1,41,1);
+
 numClusters <- max(m$c);
 
 gp_evals <- gp_log[,seq(3,202,2)];
-all_evals <- m[,6:105];
+all_evals <- m[,5:104];
+gen_i=1
+stats <- data.frame(gen=generations,
+                    numDiffClusters = -generations,
+                    medianRank=-generations)
 for(gen_i in generations) {
-#gen_i <- 15;
+ #gen_i <- 1;
   selectedRows <- seq((gen_i - 1)*popSize + 1, gen_i * popSize,1); 
-  min(selectedRows)
-  max(selectedRows)
+  #min(select2edRows)
+  #max(selectedRows)
   xcorrel <- cor(t(all_evals[,]), t(gp_evals[selectedRows,]))^2
-  mapped_gp_log <- m[max.col(t(xcorrel)), 1:5]
+  xcorrel[is.na(xcorrel)] <- 0; # for constant expressions
+  
+  idxOfMaxCorrel <- data.frame(idx=max.col(t(xcorrel)));
+  
+  # number of different nearest neighbours
+  numDistinctNeighbors <- count(distinct(.data=idxOfMaxCorrel))
+  stats$numDiffClusters[gen_i] <- numDistinctNeighbors
+  
+  mapped_gp_log <- m[idxOfMaxCorrel$idx, 1:5]
 #check
 #cor(t(all_evals[128082,]), t(gp_evals[2,]))^2
 #max(cor(t(all_evals[,]), t(gp_evals[2,]))^2)
-  ggplot(mapped_gp_log, aes(x=x, y=y)) + xlim(-75,75) + ylim(-75,75) + geom_point();
+  ggplot(mapped_gp_log, aes(x=x, y=y)) + xlim(-75,75) + ylim(-75,75) + labs(caption=paste("# distinct neighbors",numDistinctNeighbors))+ geom_point();
   ggsave(paste("scatter",gen_i,".png"))
+  
+#  bestFunc <- data.frame(x = seq(0,9.8,0.1),
+#                         t(m[idxOfMaxCorrel$idx[1], 6:104]), 
+#                         t(gp_evals[selectedRows[1],1:99]));
+#  # best neighbor
+#  ggplot(bestFunc, aes(x=x)) + 
+#    geom_line(aes(y=scale.default(bestFunc[,2]), color=rep("Nearest Neighbor",nrow(bestFunc)))) +
+#    geom_line(aes(y=scale.default(bestFunc[,3]), color=rep("Best GP Solution",nrow(bestFunc)))) + 
+#    guides(fill="legend") + 
+#    labs(y="f(x)", color="Function");
+#  ggsave(paste("best_function",gen_i,".png"))
+ 
+# count number of unclustered    
+#  count(m, m$c==-1)
 
-  ggplot(mapped_gp_log, aes(x=c)) +xlim(0,numClusters+1) + geom_histogram(binwidth = 1);
-  ggsave(paste("cluster_freq",gen_i,".png"))
+  #plot expr freq by ranks
+  mapped_gp_log_with_clusters <- dplyr::inner_join(mapped_gp_log, clusterStats, by=c("c" = "Clusters"))
+  stats$medianRank[gen_i] <- median(mapped_gp_log_with_clusters$Rank)
+  
+  mapped_gp_log_grouped <- mapped_gp_log_with_clusters %>% group_by(c) %>% summarise(count = n())
+  mapped_gp_log_grouped <- dplyr::inner_join(mapped_gp_log_grouped, clusterStats, by=c("c"="Clusters"))
+  ggplot(mapped_gp_log_grouped, aes(x=Rank,y=count)) + 
+    scale_y_log10() + 
+    xlim(0,17540) +
+    geom_point()
+  ggsave(paste("cluster-freqs",gen_i,".pdf"),
+         device=pdf,width=3)
+  
+  
+  #plot hex map
+  ggplot(mapped_gp_log, aes(x=x, y=y)) + 
+    xlim(-60,60) +
+    ylim(-60,60) + 
+    theme_void() +
+    geom_hex()
+  
+  ggsave(paste("hexbin",gen_i,".pdf"),
+         device=pdf,width=3,height=2)
+  
+#    ggplot(mapped_gp_log, aes(x=c)) +
+#  xlim(-2,numClusters+1) +
+#  geom_histogram(binwidth = 1) + 
+#  ylim(0,500);
+#  ggsave(paste("cluster_freq",gen_i,".png"))
 }
+
+# plot line chart of stats
+ggplot(stats, aes(x=gen)) +
+  xlim(1,40) + 
+  geom_line(aes(y=medianRank)) +  
+  theme_classic() +
+  labs(x="Generations", y="Median cluster rank") + 
+  scale_y_log10()
+ggsave("gp-medianClusterRank.pdf", device=pdf)
+
+print(stats$numDiffClusters)
+
+stats$numDiffClusters <- unlist(stats$numDiffClusters)
+
+ggplot(stats, aes(x=gen, y=numDiffClusters)) +
+  xlim(1,40) + 
+  geom_line() +  
+  theme_classic() +
+  labs(x="Generations", y="Number of explored clusters")
+ggsave("gp-numClusters.pdf", device=pdf)