source: branches/2886_SymRegGrammarEnumeration/ExpressionClustering_R/FindClustersForGPLog.R @ 16174

# find nearest neighbors for GP indiviuals

library(dplyr)
library(hexbin)
library(ggplot2);
#library(doSNOW); #dopar
#library(foreach); # dopar

#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,20500,500),1]
popSize <- 500;

# check qualities
# ggplot(gp_log, aes(x=idx,y=gp_log$V1)) + geom_line();

# generations <- seq(1,34000/popSize,1);
generations <- seq(1,41,1);

numClusters <- max(m$c);

gp_evals <- gp_log[,seq(3,202,2)];
all_evals <- m[,5:104];
gen_i=1
stats <- data.frame(gen=generations,
                    numDiffClusters = -generations,
                    medianRank=-generations)
for(gen_i in generations) {
 #gen_i <- 1;
  selectedRows <- seq((gen_i - 1)*popSize + 1, gen_i * popSize,1); 
  #min(select2edRows)
  #max(selectedRows)
  xcorrel <- cor(t(all_evals[,]), t(gp_evals[selectedRows,]))^2
  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) + 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)

  #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)