myprojects-hg-reborn

Fork

rm(list=ls())

library(fitdistrplus)


library(tidyverse)
library(magrittr)
library(reshape2)

library(tikzDevice)
library(scales)

source("/home/lorenzo/myprojects-hg/R-codes/stat_lib.R")


dbeta2 <- function(x, shape, ...)
       dbeta(x, shape, shape, ...)
pbeta2 <- function(q, shape, ...)
       pbeta(q, shape, shape, ...)     
rbeta2 <- function(n, shape, ...)
       rbeta(n, shape, shape, ...)     
qbeta2 <- function(p, shape, ...)
       pbeta(p, shape, shape, ...)     
       


prepare_for_fitting <- function(df21, df_name, threshold){

data <- df21  %>% mutate(df=df_name) %>% select(c(X1,X2,X4,df)) %>% set_colnames(c("n", "i", "j", "df")) %>% mutate(i_norm=i/n, j_norm=j/n) %>% filter(i!=0, j!=0) %>% arrange(n) %>% group_by(n) %>% filter(n()>=threshold) %>% ungroup

return(data)
    
}




prepare_for_fitting_large <- function(df21, df_name, threshold){

data <- df21  %>% mutate(j=n-i) %>%  mutate(df=df_name) %>% select(c(n,i,j,df)) %>% set_colnames(c("n", "i", "j", "df")) %>% mutate(i_norm=i/n, j_norm=j/n) %>% filter(i!=0, j!=0) %>% arrange(n) %>% group_by(n) %>% filter(n()>=threshold) %>% ungroup

return(data)
    
}



##########################################################################
##########################################################################
##########################################################################
##########################################################################
##########################################################################
##########################################################################


threshold=500


dname <-  "beta2"  ## "beta"



shapevec <- c()
shapevec2 <- c()
shapevec3 <- c()


d1<-read_csv("df21kf13Lib2Gen4lib1.csv", col_names=F)
d2<-read_csv("df21kf13Lib2Gen5lib1.csv", col_names=F)
d3<-read_csv("df21kf13Lib2Gen6lib1.csv", col_names=F)
d4<-read_csv("df21kf13Lib2Gen7lib1.csv", col_names=F)

## df21 <- rbind(d1,d2,d3,d4) %>% mutate(df="df21") %>% select(c(X1,X2,X4,df)) %>% set_colnames(c("n", "i", "j", "df")) %>% mutate(i_norm=i/n, j_norm=j/n) %>% filter(i!=0, j!=0) %>% arrange(n) %>% group_by(n) %>% filter(n()>=threshold) %>% ungroup

data <- rbind(d1,d2,d3,d4)


df21 <- prepare_for_fitting(data, "df_21", threshold)

data_large <- readRDS("large_fragments_df21.RDS")  %>%  as_tibble

df21_large <- prepare_for_fitting_large(data_large, "df_21", threshold)



df21 <- rbind(df21, df21_large)


###############################################################
###############################################################
###############################################################
###############################################################
###############################################################


## Now I work on another fractal dimension



d1 <- read_csv("df18kf13Lib2Gen4lib2.csv", col_names=F)
d2 <- read_csv("df18kf13Lib2Gen4lib3.csv", col_names=F)
d3 <- read_csv("df18kf13Lib2Gen5lib2.csv", col_names=F)
d4 <- read_csv("df18kf13Lib2Gen5lib3.csv", col_names=F)
d5 <- read_csv("df18kf13Lib2Gen6lib2.csv", col_names=F)
d6 <- read_csv("df18kf13Lib2Gen6lib3.csv", col_names=F)
d7 <- read_csv("df18kf13Lib2Gen7lib2.csv", col_names=F)
d8 <- read_csv("df18kf13Lib2Gen7lib3.csv", col_names=F)



data <- rbind(d1,d2,d3,d4,d5,d6,d7,d8)




df18 <- prepare_for_fitting(data, "df_18", threshold)


data_large <- readRDS("large_fragments_df18.RDS")  %>%  as_tibble

df18_large <- prepare_for_fitting_large(data_large, "df_18", threshold)



df18 <- rbind(df18, df18_large)





### and now with another one


d1 <- read_csv("df16kf13Lib2Gen4lib1.csv", col_names=F)
d2 <- read_csv("df16kf13Lib2Gen5lib1.csv", col_names=F)
d3 <- read_csv("df16kf13Lib2Gen6lib1.csv", col_names=F)
d4 <- read_csv("df16kf13Lib2Gen7lib1.csv", col_names=F)



data3 <- rbind(d1,d2,d3,d4)




df16 <- prepare_for_fitting(data3, "df_16", threshold)






########################################








sizes <- df21$n %>% unique

for (mysize in sizes){


temp <- df21 %>% filter(n==mysize) %>% select(c(n,i_norm, j_norm)) %>% melt(id.var=c("n"))

tt <- temp$value
    

## fitd <- fitdist( tt , dname, start=list(shape=0.5) )


## see http://bit.ly/2iGIrEB
    
fitd <- fitdistr( tt , dbeta2, start=list(shape=0.5) , method="Brent", lower=c(0), upper=c(1))

    

    
res <- coef(fitd) %>% as.numeric

shapevec <- c(shapevec, res[1])



    
}



output_df21 <- cbind(sizes, shapevec)  %>% as_tibble %>% set_colnames(c("N", "alpha")) %>% mutate(df="df_21")



#########################################################################



sizes2 <- df18$n %>% unique

for (mysize in sizes2){


temp <- df18 %>% filter(n==mysize) %>% select(c(n,i_norm, j_norm)) %>% melt(id.var=c("n"))

tt <- temp$value
    

## fitd <- fitdist( tt , dname, start=list(shape=0.5) )

fitd <- fitdistr( tt , dbeta2, start=list(shape=0.5) , method="Brent", lower=c(0), upper=c(1))

    
    
res <- coef(fitd) %>% as.numeric

shapevec2 <- c(shapevec2, res[1])



    
}



output_df18 <- cbind(sizes2, shapevec2)  %>% as_tibble %>% set_colnames(c("N", "alpha")) %>% mutate(df="df_18")




#########################################################################



sizes3 <- df16$n %>% unique

for (mysize in sizes3){


temp <- df16 %>% filter(n==mysize) %>% select(c(n,i_norm, j_norm)) %>% melt(id.var=c("n"))

tt <- temp$value
    

## fitd <- fitdist( tt , dname, start=list(shape=0.5) )

fitd <- fitdistr( tt , dbeta2, start=list(shape=0.5), method="Brent", lower=c(0), upper=c(1)) 

    
    
    
res <- coef(fitd) %>% as.numeric

shapevec3 <- c(shapevec3, res[1])



    
}



output_df16 <- cbind(sizes3, shapevec3)  %>% as_tibble %>% set_colnames(c("N", "alpha")) %>% mutate(df="df_16")











output_tot <- rbind(output_df18, output_df21,output_df16)



print("OK output_tot")

output_tot$df <- as.factor(output_tot$df)

lbls <- levels(output_tot$df)   ## output_tot$df %>% unique

col_seq <- seq(length(lbls))

gpl <- ggplot(output_tot, aes(x=N, y=alpha
                       , color=df, shape=df)) +

    
my_ggplot_theme(c(0.53, 0.8))+
## theme(legend.position = 'right')+
geom_point(size=3) + 

                        
#scale_y_continuous(limits=c(0.3,0.9),breaks=seq(0.3, 0.9, by=0.3))+
scale_y_continuous(breaks=pretty_breaks(n=5))+
scale_x_continuous(breaks=pretty_breaks(n=5))+

    scale_color_manual(NULL, breaks=lbls,   labels=c("$d_f=1.6$","$d_f=1.8$", "$d_f=2.1$"      ), values=col_seq) +

    scale_shape_manual(NULL, breaks=lbls,  labels= c("$d_f=1.6$","$d_f=1.8$", "$d_f=2.1$"  )   
                 ,values=col_seq) +


    
labs(title=NULL)+
theme(plot.title = element_text(lineheight=.8, size=24, face="bold", vjust=1))+
theme(legend.text = element_text(vjust=1,lineheight=1 ))+
theme(legend.title = element_text(vjust=1,lineheight=1, size=12 ))+

    xlab("$N$")+
 ylab("$\\alpha$")


fn <- paste("alpha_N_df_all.tex")

#use tikzDevice to generate and the compile the latex file ultimately
#leading to the figures.

tikz(fn, standAlone = TRUE, width=10,height=5)
print(gpl)
dev.off() 

tools::texi2dvi(fn,pdf=T)





df_plot2 <- output_tot %>% mutate(generation = ifelse(N<150, 4, ifelse( N>=150 & N<300,5, ifelse(N>=300 & N<600, 6, ifelse(N>=600 & N<1400,7, ifelse(N>=1400 & N<2800, 8, ifelse( N>=2800 & N<5500, 9 , ifelse(N>=5500 & N<8000, 10, ifelse(N>=8000, 11, NA))  )))))))


barycentre <- df_plot2 %>% group_by(generation, df) %>% summarise(mean_N=mean(N), mean_alpha=mean(alpha))



write_csv(barycentre, "barycentre.csv")

barycentre$df <- as.factor(barycentre$df)


lbls <- levels(barycentre$df)   ## unique(barycentre$df)

col_seq <- seq(length(lbls))

gpl <- ggplot(barycentre, aes(x=mean_N, y=mean_alpha
                       , color=df, shape=df)) +

    
my_ggplot_theme(c(0.53, 0.8))+
## theme(legend.position = 'right')+
geom_point(size=3) + 

                        
#scale_y_continuous(limits=c(0.3,0.9),breaks=seq(0.3, 0.9, by=0.3))+
scale_y_continuous(breaks=pretty_breaks(n=5))+
scale_x_continuous(breaks=pretty_breaks(n=5))+

    scale_color_manual(NULL, breaks=lbls,   labels=c("$d_f=1.6$","$d_f=1.8$", "$d_f=2.1$"
                                                     ), values=col_seq) +

    scale_shape_manual(NULL, breaks=lbls,  labels= c("$d_f=1.6$", "$d_f=1.8$", "$d_f=2.1$"
                                                     )   
                 ,values=col_seq) +

## coord_cartesian(xlim = c(0,800)) +
    
labs(title=NULL)+
theme(plot.title = element_text(lineheight=.8, size=24, face="bold", vjust=1))+
theme(legend.text = element_text(vjust=1,lineheight=1 ))+
theme(legend.title = element_text(vjust=1,lineheight=1, size=12 ))+

    xlab("$\\langle N \\rangle$")+
 ylab("$\\langle \\alpha  \\rangle$")


fn <- paste("alpha_N_df_barycentre.tex")

#use tikzDevice to generate and the compile the latex file ultimately
#leading to the figures.

tikz(fn, standAlone = TRUE, width=10,height=5)
print(gpl)
dev.off() 

tools::texi2dvi(fn,pdf=T)


###########################################################



df_tot <- rbind(df16,df18,df21) %>% mutate(generation = ifelse(n<150, 4, ifelse( n>=150 & n<300,5, ifelse(n>=300 & n<600, 6, ifelse(n>=600 & n<1400,7, ifelse(n>=1400 & n<2800, 8, ifelse( n>=2800 & n<5500, 9 , ifelse(n>=5500 & n<8000, 10, ifelse(n>=8000, 11, NA))  ))))))) %>% filter(i_norm>=0.5)


saveRDS(df_tot,"filtered_data.RDS")









print("So far so good")