myprojects-hg-reborn

Fork

rm(list=ls())



#library(Cairo)
#library(igraph)
#library(lattice)
#library(MASS)
library(ggplot2)
library(tikzDevice)


log_binning_3 <- function(x_min,x_max,epsi,n_bin){
#epsi <- 0.001
x_max <- x_max+epsi
m <- n_bin-1
r <- (x_max/x_min)^(1/m)
my_seq <- seq(0,m,by=1)
grid <- x_min*r^my_seq

}


slice <- 20

cut <- 3600*12 #a 12 hour cutoff


#Impose a cutoff and get a discrete distr varying in units of 1


lyon <- read.table("delta_tab_tac-lyon.dat")

sel <- which(lyon<=cut)

lyon <- lyon[sel, ]/slice

obg <- read.table("delta_tab_tac-obg.dat")

sel <- which(obg<=cut)

obg <- obg[sel,]/slice

ht <- read.table("delta_tab_tac-ht09.dat")

sel <- which(ht<=cut)

ht <- ht[sel,]/slice


sg <- read.table("delta_tab_tac-sg.dat")


sel <- which(sg<=cut)

sg <- sg[sel,]/slice


eswc <- read.table("delta_tab_tac-eswc.dat")

sel <- which(eswc<=cut)

eswc <- eswc[sel,]/slice


#generate grid for rescaled distributions

my_grid <- log_binning_3(1, cut/slice, 0,8)


#take into account the rescaled nature of the distr

my_seq <- seq(1,cut/slice)

h2<-hist(my_seq, plot=F, breaks=my_grid)



#Measure the discrete width of each bin, i.e. how many
#points can fall in each bin.

norm_count <- h2$counts





#take a few stats

h_lyon<-hist(lyon, plot=F, breaks=my_grid)

sel <- which(h_lyon$counts>0)


y_max_lyon <- h_lyon$counts[sel]/sum(h_lyon$count)/norm_count[sel]


mid_vec_lyon <-  h_lyon$mid[sel]




h_obg<-hist(obg, plot=F, breaks=my_grid)

sel <- which(h_obg$counts>0)


y_max_obg <- h_obg$counts[sel]/sum(h_obg$count)/norm_count[sel]


mid_vec_obg <-  h_obg$mid[sel]


h_ht<-hist(ht, plot=F, breaks=my_grid)

sel <- which(h_ht$counts>0)

y_max_ht <- h_ht$counts[sel]/sum(h_ht$count)/norm_count[sel]


mid_vec_ht <-  h_ht$mid[sel]



h_sg<-hist(sg, plot=F, breaks=my_grid)

sel <- which(h_sg$counts>0)

y_max_sg <- h_sg$counts[sel]/sum(h_sg$count)/norm_count[sel]



mid_vec_sg <-  h_sg$mid[sel]


h_eswc<-hist(eswc, plot=F, breaks=my_grid)

sel <- which(h_eswc$counts>0)


y_max_eswc <- h_eswc$counts[sel]/sum(h_eswc$count)/norm_count[sel]


mid_vec_eswc <-  h_eswc$mid[sel]


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


data_lyon <- as.data.frame(cbind(y_max_lyon,mid_vec_lyon, rep("Lyon", length(y_max_lyon))))

names(data_lyon) <- c("y","x","set")
data_lyon$set <- factor(data_lyon$set)




data_obg <- as.data.frame(cbind(y_max_obg,mid_vec_obg, rep("OBG", length(y_max_obg))))

names(data_obg) <- c("y","x","set")
data_obg$set <- factor(data_obg$set)


#write.table(rbind(data_obg$x*slice,data_obg$y), "data_obg_t_ab_ac.dat", row.names=FALSE, col.names=FALSE, quote=FALSE)


data_ht <- as.data.frame(cbind(y_max_ht,mid_vec_ht, rep("HT09", length(y_max_ht))))

names(data_ht) <- c("y","x","set")

data_ht$set <- factor(data_ht$set)

#write.table(rbind(data_ht$x*slice,data_ht$y), "data_ht_t_ab_ac.dat", row.names=FALSE, col.names=FALSE, quote=FALSE)


data_sg <- as.data.frame(cbind(y_max_sg,mid_vec_sg, rep("SG", length(y_max_sg))))

names(data_sg) <- c("y","x","set")

data_sg$set <- factor(data_sg$set)


#write.table(rbind(data_sg$x*slice,data_sg$y), "data_sg_t_ab_ac.dat", row.names=FALSE, col.names=FALSE, quote=FALSE)



data_eswc <- as.data.frame(cbind(y_max_eswc,mid_vec_eswc, rep("ESWC10", length(y_max_eswc))))

names(data_eswc) <- c("y","x","set")

data_eswc$set <- factor(data_eswc$set)

#write.table(rbind(data_eswc$x*slice,data_eswc*y), "data_eswc_t_ab_ac.dat", row.names=FALSE, col.names=FALSE, quote=FALSE)


data_overall <- as.data.frame(rbind(data_lyon,data_obg,data_ht,data_sg,data_eswc))

data_overall$y <- as.numeric(levels(data_overall$y))[data_overall$y]  #.     as.matrix(data_overall$y)

data_overall$x  <- as.numeric(levels(data_overall$x))[data_overall$x] # as.matrix(data_overall$x)

data_overall$x  <- data_overall$x*slice #bring everything back to the right units! 

cols <- c("OBG" = "black","ESWC10" = "red","HT09" = "green","SG" = "blue", "Lyon"="orange")

cat <- c("OBG","ESWC10","HT09","SG", "Lyon")

mb <- c(seq(10,100, by=10), seq(200,1000, by=100), seq(2000, 1e4, by=1e3),
        seq(2e4, 1e5, by=1e4)    )

my_sel <- seq(1,length(mb), by=9)

my_label2 <- rep("",length(mb))
my_label2[my_sel] <- expression(10^1 ,10^2, 10^3, 10^4, 10^5)




mby <- c(seq(1e-8,1e-7,by=1e-8),seq(2e-7,1e-6,by=1e-7),seq(2e-6, 1e-5, by=1e-6),seq(2e-5,1e-4, by=1e-5),
         seq(2e-4,1e-3, by=1e-4), seq(2e-3, 1e-2, by=1e-3),
        seq(2e-2, 1e-1, by=1e-2), seq(2e-1,1,by=1e-1)    )


my_sel <- seq(1,length(mby), by=9)

my_label2y <- rep("",length(mby))
my_label2y[my_sel] <- expression(10^-8,10^-7,10^-6 ,10^-5, 10^-4, 10^-3, 10^-2, 10^-1,10^0)




gpl <- #ggplot(data_overall, aes(group = ds, fill = set)) +
#  ggplot(data_overall, aes(group = ds, color = set, shape=set)) +
   ggplot(data_overall, aes(group=set, color = set, shape=set)) +
#  geom_rect(aes(ymin = y_min, ymax = y_max,xmin=x_min,xmax=x_max), colour="black", alpha=0.3)+
  geom_line(aes(y = y,x=x), alpha=1, size=1.5)+
  geom_point(aes(y = y,x=x), alpha=1, size=4)+
  geom_point(aes(y = y,x=x), alpha=1, size=3.8)+
  geom_point(aes(y = y,x=x), alpha=1, size=3.6)+

    ## geom_vline(xintercept = seq(0,4,by=1), col="white", size=1)+
    ## geom_hline(yintercept = seq(-6,0,by=1), col="white", size=1)+
  ## geom_point(size=3.6) +
  ## geom_point(size=3.8) +
  ##   geom_point(size=4) +
    opts( panel.background=theme_rect(fill='white'))+
    opts(legend.text = theme_text(size = 17, vjust=0.4))+
  opts(legend.title = theme_text(size = 24,hjust=0))+
   opts(legend.position = c(0.73, 0.82),legend.title = theme_blank())+
  opts( legend.background=theme_rect(col=0),
       legend.key = theme_rect(colour = NA))+
  #extra stuff
  xlab(expression(paste("$\\Delta t_{AB-AC}$ (sec)")))+
  ylab(expression(paste("$P(\\Delta t_{AB-AC})$")))+
#    opts(title = nm)+
   opts( panel.background=theme_rect(fill='white', size=1.5))+
  opts(panel.grid.minor = theme_blank())+
opts(panel.grid.major = theme_blank())+
scale_x_continuous(trans = "log10",limits=c(10,1e5),breaks=mb, label=my_label2)+
scale_y_continuous(trans="log10", limits=c(1e-8, 1e0), breaks=mby,label=my_label2y )+
    opts(axis.ticks = theme_segment(colour = "black", size=1.5),axis.ticks.length = unit(0.15, "cm"))+
#    scale_fill_manual("",labels=cat_aggregate, breaks=cat_aggregate,  values=cols) +
    scale_colour_manual("",labels=cat, breaks=cat,  values=cols) +

  scale_shape_manual("", breaks=cat, labels=cat, values=c(21,22,23,24,25))+
opts(axis.title.x = theme_text(size = 24))+
opts(axis.title.y = theme_text(size = 24, angle=90,hjust=0.3))+
#  opts(title = expression(paste("OBG, all data")))+
  opts(plot.title = theme_text(size = 25))+
opts(axis.text.x = theme_text(size=22, colour="black",vjust=1))+
opts(axis.text.y = theme_text(size=22, colour="black", hjust=1))


tikz('comparison_t_ab_ac.tex', standAlone = TRUE, width=5,height=5)
print(gpl)
dev.off() 

tools::texi2dvi('comparison_t_ab_ac.tex',pdf=T)


write.table(data_overall, "data_comparison_t_ab_ac.dat", row.names=FALSE, col.names=FALSE, quote=FALSE)



print("So far so good")