# turn on  X11 output
X11()
# this should make pause (wait for user input) but doesn't work probably because we redirect the script to the stdin
# par(ask = TRUE)

# read both data files
P1_data=read.table("./noncompatible_2_P1/summary",head=T)
P0_data=read.table("./noncompatible_2_P0/summary",head=T)

# add factors for methods and connect both data frames
method_factor=factor(c("method P0"," method P1"))
P1_data$method=rep(levels(method_factor)[2],84)
P0_data$method=rep(levels(method_factor)[1],84)
dat=rbind(P0_data,P1_data)

# add factors for sigma and ratio (h1/h2)
dat$sigma_fac=as.factor(dat$sigma)
dat$ratio_fac=as.factor(dat$ratio)

# compute log of condition number
dat$log_cond=log10(dat$eig_max/dat$eig_min)

# exclude smallest ratio, P0 behaves wrong for two small 1d elements
dat=subset(dat, ratio!=0.2)

require(ggplot2)


# for sigma==100 plot all variables  in all ratios for both methods
dat_100=subset(dat, sigma == 100)


p<-ggplot(dat_100, aes(group=ratio, color=ratio_fac)) 
p<-p+scale_y_continuous(expression(paste(log[2],"(error)"))) + scale_x_continuous(expression(paste(-log[2],h[paste(2,d)])))
p<-p+scale_linetype_manual(name="domain", values=c("dashed","solid"),labels=c("1d","2d") )
p<-p+scale_size_manual(name="variable",values=c(1,0.3),labels=c("press.","vel."))
p<-p+scale_color_manual(name=expression( bold( rho == delta[2] / delta[1] )) ,values=c("grey80","grey50","grey0"))
p<-p+geom_line( aes(exp, log2(p1d), linetype="1" , size="1" )) 
p<-p+geom_line( aes(exp, log2(p2d), linetype="2", size="1" ))
p<-p+geom_line( aes(exp, log2(v1d), linetype="1", size="2" )) 
p<-p+geom_line( aes(exp, log2(v2d), linetype="2", size="2" ))
p<-p+facet_grid(.~method, )+opts(panel.background = theme_rect(), panel.grid.minor = theme_blank(), panel.grid.major = theme_blank()  )
p
ggsave(filename="ratio_plot.pdf", width=12)


library(tcltk)
tk_messageBox(message="Press a key")

# plot condition numbers
qplot(exp,log_cond,data=dat, shape=dat$method)

q()

order_fit <- function(data, quantity) {
  q=data[[quantity]]
  l<-lm( log2(q) ~ exp, data=data, subset=( q>0 & sigma==1 & ratio==0.5) )
  print( c( coef(l)[2], vcov(l)[2,2] ) )
}
                          #  est.        std. err.
order_fit(P0_data, "p1d") # -0.32324    0.03016
order_fit(P0_data, "p2d") # -0.29207    0.03394
order_fit(P0_data, "v1d") #
order_fit(P0_data, "v2d") #

order_fit(P1_data, "p1d") #
order_fit(P1_data, "p2d") #
order_fit(P1_data, "v1d") #
order_fit(P1_data, "v2d") #

res=lm( log2(v1d) ~ exp, data=dat, subset=( sigma==100 & ratio==1 & method=="P1") )


ggplot(df, aes(x, y)) +
   geom_smooth(method = 'lm') +
   geom_text(aes(x = 2, y = 10, label = rout[[1]]), hjust = 0) +
   geom_text(aes(x = 2, y = 9.5, label = "sum(x)"), hjust = 0, parse = TRUE )
   
library(ggplot2); pdf('test.pdf'); qplot(mpg, wt, data=mtcars);
dev.off()    

ggplot(data=dat, aes(group=ratio)) +geom_line(aes(exp, log2(p1d), color="light_red") + scale_hue(h="blueratio, c=ratio , group=ratio, color=ratio))+stat_smooth(method=lm, se=F, color="red")+geom_line()+facet_grid(method~sigma)
ggplot(data=dat, aes(exp, log2(p1d), group=ratio, color=ratio)) +geom_line() + scale_color_hue(c=c(0,100))+facet_grid(method~sigma)