Table of Contents
library(ggplot2)
args <- commandArgs(trailingOnly = TRUE)
rep = paste("./",args[1],sep="")
V_REP = c(paste(rep,"/ITE/A0/1E2/WgivenTrue/",sep=""), paste(rep,"/ITE/A1/1E2/WgivenTrue/",sep=""),
paste(rep,"/ITE/A2/1E2/WgivenTrue/",sep=""), paste(rep,"/ITE/A3/1E2/WgivenTrue/",sep=""), paste(rep,"/ITE/A4/1E2/WgivenTrue/",sep=""),
paste(rep,"/ITE/A2/1E4/WgivenTrue/",sep=""), paste(rep,"/ITE/A3/1E4/WgivenTrue/",sep=""), paste(rep,"/ITE/A4/1E4/WgivenTrue/",sep=""),
paste(rep,"/ITE/A2/1E6/WgivenTrue/",sep=""), paste(rep,"/ITE/A3/1E6/WgivenTrue/",sep=""), paste(rep,"/ITE/A4/1E6/WgivenTrue/",sep=""))
ARG1 = c("S0", "S1", "S2", "S3")
ARG1 = c("S1", "S2", "S3")
ARG2 = c("LD", "MD", "LMD")
ARG3 = c(args[2])
V_FILE <- ""
V_OUT <- ""
V_PNG <- ""
list <- ""
t = 1; z = 1
for(z in 1:length(V_REP)){
t = 1
for(i in ARG1){
a = paste("*", i, sep="_")
a2 = i
if(i == "S3"){ARG22 = ARG2}
if(i != "S3"){ARG22 = c("LD")}
for(j in ARG22){
b = paste(a, j, sep="-")
b2 = paste(a2, j, sep="-")
for(k in ARG3){
c = paste(b, k, sep="-")
c2 = paste(b2, k, sep="-")
V_FILE[t] = paste(c, ".csv", sep="")
V_OUT[t] = paste(c2, ".csv", sep="")
V_PNG[t] = paste(c2, ".pdf", sep="")
t = t+1
system(command = paste(paste(paste("cat", V_REP[z],sep=" "), V_FILE[t-1], sep=""), paste(">",paste(V_REP[z],V_OUT[t-1],sep=""),sep=" "), sep=" "), intern=TRUE)
system(command = paste(paste("rm", V_REP[z],sep=" "), V_FILE[t-1], sep=""), intern = TRUE)
}
}
}
}
iliste = 0
for(IV_REP in V_REP){
ip = 0; iliste = iliste + 1; dall <- 0
for(IV_OUT in V_OUT){
ip = ip + 1
# Read Data:
filen = paste(IV_REP, IV_OUT, sep="")
df3 = read.csv(file = filen)
ite <- 0; ba <- 0.0; r <- 0.0; de <- 0.0; kw <- 0; s <- 0; enum <- 0; Si <-0; Swhich <- 0; Stheta <- 0; cpt <- 0
dimi1 = 1; dimi2 = length(df3$iterations); dimj = dimi2; dim = dimi2
ns = max(df3$enum)
z = 0; t = 1; i = 1; maxrun = 3; val = 0
for(k in seq(0, ns, by=maxrun)){
z = 0
ite <- 0; ba <- 0.0; r <- 0.0; de <- 0.0; kw <- 0; s <- 0; enum <- 0; Si <-0; Swhich <- 0; Stheta <- 0
for(i in dimi1:dimi2){
if( (df3$s[i] == k & k==0) | (df3$enum[i] %in% seq(k-maxrun+1,k,by=1)) ){
ite[z+1] = df3$iterations[i]; ba[z+1] = df3$backward_error[i]
de[z+1] = df3$direct_error[i]; r[z+1] = df3$approx_backward[i]
kw[z+1] = df3$k[i]; s[z+1] = df3$s[i]; enum[z+1] = val #df3$enum[i]
Si[z+1] = df3$Si[i]; Swhich[z+1] = df3$which[i]; Stheta[z+1] = df3$theta[i]
z = z + 1
}
}# end i
val = val + 1
z = length(s)
itec = 0 * seq(1, z, 1) - 1
bac = 0.0*itec; rc = bac; dec = bac; kwc = itec; sc = itec; enumc = itec
Sic = itec; Swhichc = itec; Sthetac = bac
cpt = 0*seq(1, z, 1)
# *Mean in Group(enum)
t = 1; i = 1
while(i <= z){
if((ite[i] %in% itec) == FALSE){
for(j in i:z){
if(ite[j] == ite[i]){
itec[t] = ite[i]
bac[t] = bac[t] + ba[j]; rc[t] = rc[t] + r[j]; dec[t] = dec[t] + de[j]
kwc[t] = kw[i]; sc[t] = s[i]; enumc[t] = enum[i]
Sic[t] = Si[i]; Swhichc[t] = Swhich[i]; Sthetac[t] = Stheta[i]
cpt[t] = cpt[t] + 1
}#end if
}#end j
t = t + 1
}#end if
i = i + 1
}#end while
ierr = 1
while(itec[ierr] > -1){
ierr = ierr + 1
}#end while
ierr = ierr - 1
# **Delate if -1
itec = itec[1:ierr]
bac = bac[1:ierr]; rc = rc[1:ierr]; dec = dec[1:ierr]
kwc = kwc[1:ierr]; sc = sc[1:ierr]; enumc = enumc[1:ierr]
Sic = Sic[1:ierr]; Swhichc = Swhichc[1:ierr]; Sthetac = Sthetac[1:ierr]
cpt = cpt[1:ierr]
# **Compute mean
bac = bac/cpt; rc = rc/cpt; dec = dec/cpt
# *Create data frame with group(s)
#dff <- data.frame(itec, bac, rc, dec, kwc, sc)
if(k == 0){
df = data.frame(itec, bac, rc, dec, kwc, sc, enumc, Sic, Swhichc, Sthetac)
}else{df = rbind(df, data.frame(itec, bac, rc, dec, kwc, sc, enumc, Sic, Swhichc, Sthetac))}
}#end k
#pname = paste(IV_REP, V_PNG[ip], sep="")
if(ip==1){
dall = df
}else{dall = rbind(dall, df)}
#print(dall)
#write.csv(df, filen)
list[iliste] = filen
#p <- ggplot(data = df, aes(x=itec, y=bac,color=kwc, group=enumc))+geom_line()+scale_color_gradient(low = "blue", high = "red")+scale_y_continuous(trans='log2')
#ggsave(pname, plot = p, width = 11, height = 8)
}
write.csv(dall, paste(IV_REP,"R.csv",sep=""))
}
library(ggplot2)
args <- commandArgs(trailingOnly = TRUE)
rep = paste("./",args[1],sep="")
rep2 = paste(paste("./",args[1],sep=""),"/ITE/",sep="")
V_REP1 = c(paste(rep,"/ITE/A0/1E2/WgivenTrue/",sep=""), paste(rep,"/ITE/A1/1E2/WgivenTrue/",sep=""))
V_REP = matrix(c(V_REP1), nrow=1, ncol=2, byrow=TRUE)
A_FILE = c("A0", "A1")
dall <- 0
i = 1
for(j in 1:2){
r = V_REP[i,j]; fn = paste(r, "R.csv", sep="")
dall = read.csv(fn)
ik=1;for(ij in dall$Sic){if(ij==1){dall$Sic[ik]="S1"};if(ij==2){dall$Sic[ik]="S2"};if(ij==3){dall$Sic[ik]="S3"};ik = ik + 1}
ik=1;for(ij in dall$Swhichc){if(ij==1){dall$Swhichc[ik]="LD"};if(ij==2){dall$Swhichc[ik]="MD"};if(ij==3){dall$Swhichc[ik]="LMD"};ik = ik + 1}
dall <- subset(dall, Sic!="S1")
#BACKWARD ERROR:
p <- ggplot(data = dall, aes(x=itec, y=bac, color=kwc, group=enumc))+geom_line()+scale_y_continuous(trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15,1E-16))
p <- p+facet_grid(. ~ Sic + Swhichc)+labs(y="Backward error")+labs(x="Solver iter.")+guides(color=guide_legend(title="k"))+scale_color_gradient(low = "blue", high = "red")
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
r = paste(paste(paste(rep2,A_FILE[j],sep=""),"/BA",sep=""),paste(A_FILE[j],".pdf",sep=""),sep="_")
ggsave(r, plot = p, width = 8, height = 6)
r = paste(paste(paste(rep2,A_FILE[j],sep=""),"/BA",sep=""),paste(A_FILE[j],".png",sep=""),sep="_")
ggsave(r, plot = p, width = 8, height = 6)
#DIRECT ERROR:
p <- ggplot(data = dall, aes(x=itec, y=dec, color=kwc, group=enumc))+geom_line()+scale_y_continuous(trans='log2',breaks=c(1,1E-4,1E-8,1E-13))
p <- p+facet_grid(. ~ Sic + Swhichc)+labs(y="Direct error (A-norm)")+labs(x="Solver iter.")+guides(color=guide_legend(title="k"))+scale_color_gradient(low = "blue", high = "red")
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
r = paste(paste(paste(rep2,A_FILE[j],sep=""),"/DE",sep=""),paste(A_FILE[j],".pdf",sep=""),sep="_")
ggsave(r, plot = p, width = 8, height = 6)
r = paste(paste(paste(rep2,A_FILE[j],sep=""),"/DE",sep=""),paste(A_FILE[j],".png",sep=""),sep="_")
ggsave(r, plot = p, width = 8, height = 6)
#ITERATED BACKWARD ERROR:
p <- ggplot(data = dall, aes(x=itec, y=rc, color=kwc, group=enumc))+geom_line()+scale_y_continuous(trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15,1E-16))
p <- p+facet_grid(. ~ Sic + Swhichc)+labs(y="Iterated backward error")+labs(x="Solver iter.")+guides(color=guide_legend(title="k"))+scale_color_gradient(low = "blue", high = "red")
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
r = paste(paste(paste(rep2,A_FILE[j],sep=""),"/IBA",sep=""),paste(A_FILE[j],".pdf",sep=""),sep="_")
ggsave(r, plot = p, width = 8, height = 6)
r = paste(paste(paste(rep2,A_FILE[j],sep=""),"/IBA",sep=""),paste(A_FILE[j],".png",sep=""),sep="_")
ggsave(r, plot = p, width = 8, height = 6)
}
V_REP1 = c(paste(rep,"/ITE/A2/1E2/WgivenTrue/",sep=""),paste(rep,"/ITE/A2/1E4/WgivenTrue/",sep=""),paste(rep,"/ITE/A2/1E6/WgivenTrue/",sep=""))
V_REP2 = c(paste(rep,"/ITE/A3/1E2/WgivenTrue/",sep=""), paste(rep,"/ITE/A3/1E4/WgivenTrue/",sep=""),paste(rep,"/ITE/A3/1E6/WgivenTrue/",sep=""))
V_REP3 = c(paste(rep,"/ITE/A4/1E2/WgivenTrue/",sep=""),paste(rep,"/ITE/A4/1E4/WgivenTrue/",sep=""),paste(rep,"/ITE/A4/1E6/WgivenTrue/",sep=""))
V_REP = matrix(c(V_REP1, V_REP2, V_REP3), nrow=3, ncol=3, byrow=TRUE)
A_FILE = c("A2", "A3", "A4")
dall <- 0
for(i in 1:3){
for(j in 1:3 ){
r = V_REP[i, j]; fn = paste(r, "R.csv",sep=""); dfff = read.csv(fn)
if(j==1){dall = read.csv(fn)}else{dall = rbind(dall, dfff)}
}
ik=1;for(ij in dall$Sic){if(ij==1){dall$Sic[ik]="S1"};if(ij==2){dall$Sic[ik]="S2"};if(ij==3){dall$Sic[ik]="S3"};ik = ik + 1}
ik=1;for(ij in dall$Swhichc){if(ij==1){dall$Swhichc[ik]="LD"};if(ij==2){dall$Swhichc[ik]="MD"};if(ij==3){dall$Swhichc[ik]="LMD"};ik = ik + 1}
ik=1;for(ij in dall$Sthetac){if(ij==1){dall$Sthetac[ik]="1E2"};if(ij==2){dall$Sthetac[ik]="1E4"};if(ij==3){dall$Sthetac[ik]="1E6"};ik = ik + 1}
xlim = 0
for(ik in 1:length(dall$kwc)){if(dall$kwc[ik]==0){xlim =max(xlim,dall$itec[ik])}}
xlim = xlim + 2
dall$lt <- 0
dall <- subset(dall, Sic!="S1")
for(ik in 1:length(dall$kwc)){if(dall$kwc[ik]==1){dall$lt[ik]=1};if(dall$kwc[ik]==2){dall$lt[ik]=2};if(dall$kwc[ik]==3){dall$lt[ik]=3};if(dall$kwc[ik]==4){dall$lt[ik]=4}}
for(iz in 1:2){
#p <- ggplot(data = dall, aes(x=itec, y=bac, color=factor(kwc), group=enumc))+geom_line()+scale_y_continuous(limits=c(1E-16, 1),trans='log2')
#p <- p+facet_grid(Sthetac ~ Sic + Swhichc)+labs(y="backward error")+labs(x="solver iter.")+guides(color=guide_legend(title="k"))
#BACKWARD ERROR:
p <- ggplot(data = dall, aes(x=itec, y=bac, color=kwc, group=enumc))+geom_line(aes(linetype=factor(lt)))+scale_y_continuous(trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15))+ scale_x_continuous(limits=c(0, xlim))
p <- p+facet_grid(Sthetac ~ Sic + Swhichc)+labs(y="Backward error")+labs(x="Solver iter.")+guides(color=guide_legend(title="k"), linetype=guide_legend(title="k"))+scale_color_gradient(low = "blue", high = "red")
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
if(iz==2){
r = paste(paste(paste(rep2,A_FILE[i],sep=""),paste("/ZOOM","BA",sep="_"),sep=""),paste(A_FILE[i],".pdf",sep=""),sep="_")
p <- p + scale_x_continuous(limits=c(0, 5))+ scale_y_continuous(trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15))
}else{r = paste(paste(paste(rep2,A_FILE[i],sep=""),"/BA",sep=""),paste(A_FILE[i],".pdf",sep=""),sep="_")}
ggsave(r, plot = p, width = 8, height = 6)
if(iz==2){
r = paste(paste(paste(rep2,A_FILE[i],sep=""),paste("/ZOOM","BA",sep="_"),sep=""),paste(A_FILE[i],".png",sep=""),sep="_")
p <- p + scale_x_continuous(limits=c(0, 5))+ scale_y_continuous(trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15))
}else{r = paste(paste(paste(rep2,A_FILE[i],sep=""),"/BA",sep=""),paste(A_FILE[i],".png",sep=""),sep="_")}
ggsave(r, plot = p, width = 8, height = 6)
#DIRECT ERROR (A-NORM):
p <- ggplot(data = dall, aes(x=itec, y=dec, color=kwc, group=enumc))+geom_line(aes(linetype=factor(lt)))+scale_y_continuous(trans='log2',breaks=c(1,1E-4,1E-8,1E-13))+ scale_x_continuous(limits=c(0, xlim))
p <- p+facet_grid(Sthetac ~ Sic + Swhichc)+labs(y="Direct error (A-norm)")+labs(x="Solver iter.")+guides(color=guide_legend(title="k"), linetype=guide_legend(title="k"))+scale_color_gradient(low = "blue", high = "red")
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
if(iz==2){
r = paste(paste(paste(rep2,A_FILE[i],sep=""),paste("/ZOOM","DE",sep="_"),sep=""),paste(A_FILE[i],".pdf",sep=""),sep="_")
p <- p + scale_x_continuous(limits=c(0, 5))+ scale_y_continuous(trans='log2',breaks=c(1,1E-4,1E-8,1E-13))
}else{r = paste(paste(paste(rep2,A_FILE[i],sep=""),"/DE",sep=""),paste(A_FILE[i],".pdf",sep=""),sep="_")}
ggsave(r, plot = p, width = 8, height = 6)
if(iz==2){
r = paste(paste(paste(rep2,A_FILE[i],sep=""),paste("/ZOOM","DE",sep="_"),sep=""),paste(A_FILE[i],".png",sep=""),sep="_")
p <- p + scale_x_continuous(limits=c(0, 5))+ scale_y_continuous(trans='log2',breaks=c(1,1E-4,1E-8,1E-13))
}else{r = paste(paste(paste(rep2,A_FILE[i],sep=""),"/DE",sep=""),paste(A_FILE[i],".png",sep=""),sep="_")}
ggsave(r, plot = p, width = 8, height = 6)
#ITERATED BACKWARD ERROR:
p <- ggplot(data = dall, aes(x=itec, y=rc, color=kwc, group=enumc))+geom_line(aes(linetype=factor(lt)))+scale_y_continuous(trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15))+ scale_x_continuous(limits=c(0, xlim))
p <- p+facet_grid(Sthetac ~ Sic + Swhichc)+labs(y="Iterated bacward error")+labs(x="Solver iter.")+guides(color=guide_legend(title="k"), linetype=guide_legend(title="k"))+scale_color_gradient(low = "blue", high = "red")
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
if(iz==2){
r = paste(paste(paste(rep2,A_FILE[i],sep=""),paste("/ZOOM","IBA",sep="_"),sep=""),paste(A_FILE[i],".pdf",sep=""),sep="_")
p <- p + scale_x_continuous(limits=c(0, 5))+ scale_y_continuous(trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15))
}else{r = paste(paste(paste(rep2,A_FILE[i],sep=""),"/IBA",sep=""),paste(A_FILE[i],".pdf",sep=""),sep="_")}
ggsave(r, plot = p, width = 8, height = 6)
if(iz==2){
r = paste(paste(paste(rep2,A_FILE[i],sep=""),paste("/ZOOM","IBA",sep="_"),sep=""),paste(A_FILE[i],".png",sep=""),sep="_")
p <- p + scale_x_continuous(limits=c(0, 5))+ scale_y_continuous(trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15))
}else{r = paste(paste(paste(rep2,A_FILE[i],sep=""),"/IBA",sep=""),paste(A_FILE[i],".png",sep=""),sep="_")}
ggsave(r, plot = p, width = 8, height = 6)
}
}
library(ggplot2)
args <- commandArgs(trailingOnly = TRUE)
rep = paste("./",args[1],sep="")
V_REP = c(paste(rep,"/ITE/A2/1E2/WgivenTrue/",sep=""), paste(rep,"/ITE/A2/1E4/WgivenTrue/",sep=""), paste(rep,"/ITE/A2/1E6/WgivenTrue/",sep=""),
paste(rep,"/ITE/A3/1E2/WgivenTrue/",sep=""), paste(rep,"/ITE/A3/1E4/WgivenTrue/",sep=""), paste(rep,"/ITE/A3/1E6/WgivenTrue/",sep=""))
ARG1 = c("S3")
ARG2 = c("LD", "MD")
ARG3 = c(args[2])
V_FILE <- ""
V_OUT <- ""
V_PNG <- ""
list <- ""
t = 1; z = 1
for(z in 1:length(V_REP)){
t = 1
for(i in ARG1){
a = paste("*", i, sep="_")
a2 = i
if(i == "S3"){ARG22 = ARG2}
if(i != "S3"){ARG22 = c("LD")}
for(j in ARG22){
b = paste(a, j, sep="-")
b2 = paste(a2, j, sep="-")
for(k in ARG3){
c = paste(b, k, sep="-")
c2 = paste(b2, k, sep="-")
V_FILE[t] = paste(c, ".csv", sep="")
V_OUT[t] = paste(c2, ".csv", sep="")
V_PNG[t] = paste(c2, ".pdf", sep="")
t = t+1
system(command = paste(paste(paste("cat", V_REP[z],sep=" "), V_FILE[t-1], sep=""), paste(">",paste(V_REP[z],V_OUT[t-1],sep=""),sep=" "), sep=" "), intern=TRUE)
system(command = paste(paste("rm", V_REP[z],sep=" "), V_FILE[t-1], sep=""), intern = TRUE)
}
}
}
}
ext = paste(args[2],".csv",sep="")
F_REP = c(paste(rep,paste("/ITE/A2/1E2/WgivenTrue/S3-LD-",ext,sep=""),sep=""), paste(rep,paste("/ITE/A2/1E4/WgivenTrue/S3-LD-"),ext,sep=""), paste(rep,paste("/ITE/A2/1E6/WgivenTrue/S3-LD-",ext,sep=""),sep=""),
paste(rep,paste("/ITE/A3/1E2/WgivenTrue/S3-MD-",ext,sep=""),sep=""), paste(rep,paste("/ITE/A3/1E4/WgivenTrue/S3-MD-",ext,sep=""),sep=""), paste(rep,paste("/ITE/A3/1E6/WgivenTrue/S3-MD-",ext,sep=""),sep=""))
ip = 0
for(fn in F_REP){
ip = ip + 1
df = read.csv(fn)
dred <- subset(df, k == 0)
df <- subset(df, k != 0)
dba <- aggregate(x=dred$backward_error, by=list(dred$iterations), FUN=mean)
dr <- aggregate(x=dred$approx_backward, by=list(dred$iterations), FUN=mean)
dde <- aggregate(x=dred$direct_error, by=list(dred$iterations), FUN=mean)
dred = dred[,-2]
t = max(dred$iterations)
t = t + 1
dred = dred[1:t,]
dred$backward_error = dba$x
dred$approx_backward = dr$x
dred$direct_error = dde$x
for(iz in dred$iterations){
if(iz==max(dred$iterations)){dred = dred[1:iz,]}
}
df = rbind(df,dred)
df = df[1:length(df$iterations)-1,]
df$Sin = df$enum
df = df[,-7]
df$lt <- 0
for(ik in 1:length(df$k)){if(df$k[ik]==0){df$lt[ik]=4};if(df$k[ik]==1){df$lt[ik]=1}}
write.csv(df, fn)
if(ip<=3){val=c(1e-6,1e-5,1e-4,1e-3,1e-2,1e-1,0.25,0.5,0.75,0.90,1.0)}
else{val=c(1e-6,1e-5,1e-4,1e-3,1e-2,1e-1,0.25,0.5,0.75,0.90,0.95,1.0)}
p <- ggplot(data = df, aes(x=iterations, y=backward_error, color=Sin, group=Sin))+geom_line(linetype=df$lt)+scale_y_continuous(trans='log', breaks=c(1,1E-4,1E-8,1E-12,1E-15))
p <- p+labs(y="Backward error")+labs(x="Solver iter.")+scale_color_gradient(low = "red", high = "blue", trans="log", breaks = val)
p <- p + guides(color=guide_legend(title="Sin<(w,u)"))
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
r = paste(fn,"-BA.pdf",sep="")
ggsave(r, plot = p, width = 8, height = 6)
r = paste(fn,"-BA.png",sep="")
ggsave(r, plot = p, width = 8, height = 6)
p <- ggplot(data = df, aes(x=iterations, y=direct_error, color=Sin, group=Sin))+geom_line(linetype=df$lt)+scale_y_continuous(trans='log', breaks=c(1,1E-4,1E-8,1E-13))
p <- p+labs(y="Direct error (A-norm)")+labs(x="Solver iter.")+scale_color_gradient(low = "red", high = "blue", trans="log", breaks = val)
p <- p + guides(color=guide_legend(title="Sin<(w,u)"))
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
r = paste(fn,"-DE.pdf",sep="")
ggsave(r, plot = p, width = 8, height = 6)
r = paste(fn,"-DE.png",sep="")
ggsave(r, plot = p, width = 8, height = 6)
p <- ggplot(data = df, aes(x=iterations, y=approx_backward, color=Sin, group=Sin))+geom_line(linetype=df$lt)+scale_y_continuous(trans='log', breaks=c(1,1E-4,1E-8,1E-12,1E-15))
p <- p+labs(y="Iterated backward error")+labs(x="Solver iter.")+scale_color_gradient(low = "red", high = "blue", trans="log", breaks = val)
p <- p + guides(color=guide_legend(title="Sin<(w,u)"))
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
r = paste(fn,"-IBA.pdf",sep="")
ggsave(r, plot = p, width = 8, height = 6)
r = paste(fn,"-IBA.png",sep="")
ggsave(r, plot = p, width = 8, height = 6)
if(ip == 1 | ip == 4){dall = df}else{dall = rbind(dall,df)}
if(ip == 3 | ip == 6){
p <- ggplot(data = dall, aes(x=iterations, y=backward_error, color=Sin, group=Sin))+geom_line(linetype=dall$lt)+scale_y_continuous(trans='log', breaks=c(1,1E-4,1E-8,1E-12,1E-15))
p <- p + facet_grid(~ theta)
p <- p+labs(y="Backward error")+labs(x="Solver iter.")+scale_color_gradient(low = "red", high = "blue", trans="log", breaks = val)
p <- p + guides(color=guide_legend(title="Sin<(w,u)"))
r = paste(fn,"-all-BA.pdf",sep="")
ggsave(r, plot = p, width = 8, height = 6)
r = paste(fn,"-all-BA.png",sep="")
ggsave(r, plot = p, width = 8, height = 6)
p <- ggplot(data = dall, aes(x=iterations, y=direct_error, color=Sin, group=Sin))+geom_line(linetype=dall$lt)+scale_y_continuous(trans='log', breaks=c(1,1E-4,1E-8,1E-13))
p <- p + facet_grid(~ theta)
p <- p+labs(y="Direct error (A-norm)")+labs(x="Solver iter.")+scale_color_gradient(low = "red", high = "blue", trans="log", breaks = val)
p <- p + guides(color=guide_legend(title="Sin<(w,u)"))
r = paste(fn,"-all-DE.pdf",sep="")
ggsave(r, plot = p, width = 8, height = 6)
r = paste(fn,"-all-DE.png",sep="")
ggsave(r, plot = p, width = 8, height = 6)
p <- ggplot(data = dall, aes(x=iterations, y=approx_backward, color=Sin, group=Sin))+geom_line(linetype=dall$lt)+scale_y_continuous(trans='log', breaks=c(1,1E-4,1E-8,1E-12,1E-15))
p <- p + facet_grid(~ theta)
p <- p+labs(y="Iterated backward error")+labs(x="Solver iter.")+scale_color_gradient(low = "red", high = "blue", trans="log", breaks = val)
p <- p + guides(color=guide_legend(title="Sin<(w,u)"))
r = paste(fn,"-all-IBA.pdf",sep="")
ggsave(r, plot = p, width = 8, height = 6)
r = paste(fn,"-all-IBA.png",sep="")
ggsave(r, plot = p, width = 8, height = 6)
}
}
library(ggplot2)
args <- commandArgs(trailingOnly = TRUE)
rep = paste("./",args[3],sep="")
rep2 = paste(paste("./",args[3],sep=""),"/ITE/",sep="")
V_REP1 = c(paste(rep,"/ITE/A0/1E2/WgivenFalse/",sep=""), paste(rep,"/ITE/A1/1E2/WgivenFalse/",sep=""))
V_REP = matrix(c(V_REP1), nrow=1, ncol=2, byrow=TRUE)
A_FILE = c("A0", "A1")
dall <- 0
outfn1 = paste(paste("L_",args[1],sep=""),args[2],sep="_")
i = 1
for(j in 1:2){
r = V_REP[i,j]; fn = paste(r, "R.csv", sep="")
dall = read.csv(fn)
ik=1;for(ij in dall$Sic){if(ij==1){dall$Sic[ik]="S1"};if(ij==2){dall$Sic[ik]="S2"};if(ij==3){dall$Sic[ik]="S3"};ik = ik + 1}
ik=1;for(ij in dall$Swhichc){if(ij==1){dall$Swhichc[ik]="LD"};if(ij==2){dall$Swhichc[ik]="MD"};if(ij==3){dall$Swhichc[ik]="LMD"};ik = ik + 1}
#BACKWARD ERROR:
p <- ggplot(data = dall, aes(x=itec, y=bac, color=kwc, group=enumc))+geom_line()+scale_y_continuous(limits=c(1E-16, 1),trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15))
p <- p+facet_grid(. ~ Sic + Swhichc)+labs(y="Backward error")+labs(x="Solver iter.")+guides(color=guide_legend(title="k"))+scale_color_gradient(low = "blue", high = "red")
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
r = paste(paste(paste(rep2,A_FILE[j],sep=""),"/BA",sep=""),paste(paste(A_FILE[j],outfn1,sep="_"),".pdf",sep=""),sep="_")
ggsave(r, plot = p, width = 8, height = 6)
r = paste(paste(paste(rep2,A_FILE[j],sep=""),"/BA",sep=""),paste(paste(A_FILE[j],outfn1,sep="_"),".png",sep=""),sep="_")
ggsave(r, plot = p, width = 8, height = 6)
#DIRECT ERROR:
p <- ggplot(data = dall, aes(x=itec, y=dec, color=kwc, group=enumc))+geom_line()+scale_y_continuous(trans='log2', breaks=c(1,1E-4,1E-8,1E-13))
p <- p+facet_grid(. ~ Sic + Swhichc)+labs(y="Direct error (A-norm)")+labs(x="Solver iter.")+guides(color=guide_legend(title="k"))+scale_color_gradient(low = "blue", high = "red")
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
r = paste(paste(paste(rep2,A_FILE[j],sep=""),"/DE",sep=""),paste(paste(A_FILE[j],outfn1,sep="_"),".pdf",sep=""),sep="_")
ggsave(r, plot = p, width = 8, height = 6)
r = paste(paste(paste(rep2,A_FILE[j],sep=""),"/DE",sep=""),paste(paste(A_FILE[j],outfn1,sep="_"),".png",sep=""),sep="_")
ggsave(r, plot = p, width = 8, height = 6)
#ITERATED BACKWARD ERROR:
p <- ggplot(data = dall, aes(x=itec, y=rc, color=kwc, group=enumc))+geom_line()+scale_y_continuous(trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15))
p <- p+facet_grid(. ~ Sic + Swhichc)+labs(y="Iterated backward error")+labs(x="Solver iter.")+guides(color=guide_legend(title="k"))+scale_color_gradient(low = "blue", high = "red")
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
r = paste(paste(paste(rep2,A_FILE[j],sep=""),"/IBA",sep=""),paste(paste(A_FILE[j],outfn1,sep="_"),".pdf",sep=""),sep="_")
ggsave(r, plot = p, width = 8, height = 6)
r = paste(paste(paste(rep2,A_FILE[j],sep=""),"/IBA",sep=""),paste(paste(A_FILE[j],outfn1,sep="_"),".png",sep=""),sep="_")
ggsave(r, plot = p, width = 8, height = 6)
}
V_REP1 = c(paste(rep,"/ITE/A2/1E2/WgivenFalse/",sep=""), paste(rep,"/ITE/A2/1E4/WgivenFalse/",sep=""), paste(rep,"/ITE/A2/1E6/WgivenFalse/",sep=""))
V_REP2 = c(paste(rep,"/ITE/A3/1E2/WgivenFalse/",sep=""), paste(rep,"/ITE/A3/1E4/WgivenFalse/",sep=""), paste(rep,"/ITE/A3/1E6/WgivenFalse/",sep=""))
V_REP3 = c(paste(rep,"/ITE/A4/1E2/WgivenFalse/",sep=""), paste(rep,"/ITE/A4/1E4/WgivenFalse/",sep=""), paste(rep,"/ITE/A4/1E6/WgivenFalse/",sep=""))
V_REP = matrix(c(V_REP1, V_REP2, V_REP3), nrow=3, ncol=3, byrow=TRUE)
A_FILE = c("A2", "A3", "A4")
dall <- 0
for(i in 1:3){
for(j in 1:3 ){
r = V_REP[i, j]; fn = paste(r, "R.csv",sep=""); dfff = read.csv(fn)
if(j==1){dall = read.csv(fn)}else{dall = rbind(dall, dfff)}
}
ik=1;for(ij in dall$Sic){if(ij==1){dall$Sic[ik]="S1"};if(ij==2){dall$Sic[ik]="S2"};if(ij==3){dall$Sic[ik]="S3"};ik = ik + 1}
ik=1;for(ij in dall$Swhichc){if(ij==1){dall$Swhichc[ik]="LD"};if(ij==2){dall$Swhichc[ik]="MD"};if(ij==3){dall$Swhichc[ik]="LMD"};ik = ik + 1}
ik=1;for(ij in dall$Sthetac){if(ij==1){dall$Sthetac[ik]="1E2"};if(ij==2){dall$Sthetac[ik]="1E4"};if(ij==3){dall$Sthetac[ik]="1E6"};ik = ik + 1}
outfn = paste(paste(paste("L_",args[1],sep=""),args[2],sep="_"),".csv",sep="")
print(outfn)
dall$l = paste(args[1],args[2],sep=" - ")
for(ik in 1:length(dall$l)){if(dall$Sic[ik]=="S1" | dall$Sic[ik]=="S2"){dall$l[ik]="-"}}
write.csv(dall, paste(paste(rep2,A_FILE[i],sep=""),outfn,sep="/L/"))
for(ik in 1:length(dall$kwc)){if(dall$kwc[ik]==0 & dall$sc[ik]>1){dall$kwc[ik]=dall$kwc[ik-1]}}
xlim = 0
for(ik in 1:length(dall$kwc)){if(dall$kwc[ik]==0){xlim =max(xlim,dall$itec[ik])}}
xlim = xlim + 2
dall$lt <- 0
for(ik in 1:length(dall$kwc)){if(dall$kwc[ik]==1){dall$lt[ik]=1};if(dall$kwc[ik]==2){dall$lt[ik]=2};if(dall$kwc[ik]==3){dall$lt[ik]=3};if(dall$kwc[ik]==4){dall$lt[ik]=4}}
dcopy <- dall
for(isi in 1:2){
dall <- dcopy
if(isi==1){dall <- subset(dcopy, Sic=="S1");dall <- rbind(dall,subset(dcopy, Sic=="S2")); isep = paste("_","1",sep="")}
if(isi==2){dall <- subset(dcopy, Sic!="S1");dall <- subset(dall, Sic!="S2"); isep = paste("_","1",sep=""); isep = paste("_","2",sep="")}
for(iz in 1:2){
#p <- ggplot(data = dall, aes(x=itec, y=bac, color=factor(kwc), group=enumc))+geom_line()+scale_y_continuous(trans='log2')
#p <- p+facet_grid(Sthetac ~ Sic + Swhichc)+labs(y="backward error")+labs(x="solver iter.")+guides(color=guide_legend(title="k"))
#BACKWARD ERROR:
p <- ggplot(data = dall, aes(x=itec, y=bac, color=kwc, group=enumc))+geom_line(aes(linetype=factor(lt)))+scale_y_continuous(limits=c(1E-16, 1),trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15))+ scale_x_continuous(limits=c(0, xlim))
p <- p+facet_grid(Sthetac ~ Sic + Swhichc+ l)+labs(y="Backward error")+labs(x="Solver iter.")+guides(color=guide_legend(title="k"), linetype=guide_legend(title="k"))+scale_color_gradient(low = "blue", high = "red")
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
if(iz==2){
r = paste(paste(paste(rep2,A_FILE[i],sep=""),paste("/ZOOM","BA",sep="_"),sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".pdf",sep=isep),sep="_")
p <- p + scale_x_continuous(limits=c(0, 7))+ scale_y_continuous(trans='log2')
}else{r = paste(paste(paste(rep2,A_FILE[i],sep=""),"/BA",sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".pdf",sep=isep),sep="_")}
ggsave(r, plot = p, width = 8, height = 6)
if(iz==2){
r = paste(paste(paste(rep2,A_FILE[i],sep=""),paste("/ZOOM","BA",sep="_"),sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".pdf",sep=isep),sep="_")
p <- p + scale_x_continuous(limits=c(0, 7))+ scale_y_continuous(trans='log2')
}else{r = paste(paste(paste(rep2,A_FILE[i],sep=""),"/BA",sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".png",sep=isep),sep="_")}
ggsave(r, plot = p, width = 8, height = 6)
#DIRECT ERROR (A-NORM):
p <- ggplot(data = dall, aes(x=itec, y=dec, color=kwc, group=enumc))+geom_line(aes(linetype=factor(lt)))+scale_y_continuous(trans='log2', breaks=c(1,1E-4,1E-8,1E-13))+ scale_x_continuous(limits=c(0, xlim))
p <- p+facet_grid(Sthetac ~ Sic + Swhichc + l)+labs(y="Direct error (A-norm)")+labs(x="Solver iter.")+guides(color=guide_legend(title="k"), linetype=guide_legend(title="k"))+scale_color_gradient(low = "blue", high = "red")
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
if(iz==2){
r = paste(paste(paste(rep2,A_FILE[i],sep=""),paste("/ZOOM","DE",sep="_"),sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".pdf",sep=isep),sep="_")
p <- p + scale_x_continuous(limits=c(0, 5))+ scale_y_continuous(trans='log2')
}else{r = paste(paste(paste(rep2,A_FILE[i],sep=""),"/DE",sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".pdf",sep=isep),sep="_")}
ggsave(r, plot = p, width = 8, height = 6)
if(iz==2){
r = paste(paste(paste(rep2,A_FILE[i],sep=""),paste("/ZOOM","DE",sep="_"),sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".pdf",sep=isep),sep="_")
p <- p + scale_x_continuous(limits=c(0, 5))+ scale_y_continuous(trans='log2')
}else{r = paste(paste(paste(rep2,A_FILE[i],sep=""),"/DE",sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".png",sep=isep),sep="_")}
ggsave(r, plot = p, width = 8, height = 6)
#ITERATED BACKWARD ERROR:
p <- ggplot(data = dall, aes(x=itec, y=rc, color=kwc, group=enumc))+geom_line(aes(linetype=factor(lt)))+scale_y_continuous(trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15))+ scale_x_continuous(limits=c(0, xlim))
p <- p+facet_grid(Sthetac ~ Sic + Swhichc + l)+labs(y="Iterated bacward error")+labs(x="Solver iter.")+guides(color=guide_legend(title="k"), linetype=guide_legend(title="k"))+scale_color_gradient(low = "blue", high = "red")
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
if(iz==2){
r = paste(paste(paste(rep2,A_FILE[i],sep=""),paste("/ZOOM","IBA",sep="_"),sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".pdf",sep=isep),sep="_")
p <- p + scale_x_continuous(limits=c(0, 5))+ scale_y_continuous(trans='log2')
}else{r = paste(paste(paste(rep2,A_FILE[i],sep=""),"/IBA",sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".pdf",sep=""),sep=isep)}
ggsave(r, plot = p, width = 8, height = 6)
if(iz==2){
r = paste(paste(paste(rep2,A_FILE[i],sep=""),paste("/ZOOM","IBA",sep="_"),sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".pdf",sep=isep),sep="_")
p <- p + scale_x_continuous(limits=c(0, 5))+ scale_y_continuous(trans='log2')
}else{r = paste(paste(paste(rep2,A_FILE[i],sep=""),"/IBA",sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".png",sep=isep),sep="_")}
ggsave(r, plot = p, width = 8, height = 6)
}
}
}
library(ggplot2)
args <- commandArgs(trailingOnly = TRUE)
rep = paste("./",args[1],sep="")
V_REP1 = c(paste(rep,"/ITE/A2/L/",sep=""), paste(rep,"/ITE/A3/L/",sep=""), paste(rep,"/ITE/A4/L/",sep=""))
V_REP = matrix(c(V_REP1), nrow=1, ncol=3, byrow=TRUE)
A_FILE = c("A2", "A3", "A4")
pp <- 0
print(getwd())
for(j in 1:3){
setwd(V_REP[1,j])
print(getwd())
liste = dir()
print(liste)
ns = length(liste)
for(ifn in 1:ns){
print(liste[ifn])
dff=read.csv(liste[ifn])
if(ifn == 1){df = read.csv(liste[ifn])}
else{df = rbind(df, dff)}
}
df <- subset(df, Sic != "S1")
df <- subset(df, Sic != "S2")
xlim = 0
for(ik in 1:length(df$kwc)){if(df$kwc[ik]==0){xlim =max(xlim,df$itec[ik])}}
xlim = xlim + 2
#df$lt <- 0
#for(ik in 1:length(df$kwc)){if(df$kwc[ik]==1){df$lt[ik]=1};if(df$kwc[ik]==1){df$lt[ik]=2};if(df$kwc[ik]==3){df$lt[ik]=3};if(df$kwc[ik]==4){df$lt[ik]=4}}
for(z in 1:2){
if(z==1){df1 <- subset(df, Swhichc != "LMD");df1 <- subset(df1, Swhichc != "MD");dff <- df1}
if(z==2){df2 <- subset(df, Swhichc != "LMD");df2 <- subset(df2, Swhichc != "LD");dff <- df2}
#BACKWARD ERROR
p <- ggplot(data = dff, aes(x=itec, y=bac, color=kwc, group=enumc))+geom_line()+scale_y_continuous(limits=c(1E-16, 1),trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15))+ scale_x_continuous(limits=c(0, xlim))
p <- p+facet_grid(Sthetac ~ Sic+Swhichc+l)+labs(y="Backward error")+labs(x="Solver iter.")+guides(color=guide_legend(title="k"))+scale_color_gradient(low = "blue", high = "red")
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
if(z==1){ggsave(paste("../LD",paste(j,"BA.pdf",sep=""),sep=""), plot = p, width = 8, height = 6)}
if(z==2){ggsave(paste("../MD",paste(j,"BA.pdf",sep=""),sep=""), plot = p, width = 8, height = 6)}
if(z==1){ggsave(paste("../LD",paste(j,"BA.png",sep=""),sep=""), plot = p, width = 8, height = 6)}
if(z==2){ggsave(paste("../MD",paste(j,"BA.png",sep=""),sep=""), plot = p, width = 8, height = 6)}
#DIRECT ERROR
p <- ggplot(data = dff, aes(x=itec, y=dec, color=kwc, group=enumc))+geom_line()+scale_y_continuous(limits=c(1E-16, 1),trans='log2', breaks=c(1,1E-4,1E-8,1E-13))+ scale_x_continuous(limits=c(0, xlim))
p <- p+facet_grid(Sthetac ~ Sic+Swhichc+l)+labs(y="Direct error (A-norm)")+labs(x="Solver iter.")+guides(color=guide_legend(title="k"))+scale_color_gradient(low = "blue", high = "red")
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
if(z==1){ggsave(paste("../LD",paste(j,"DE.pdf",sep=""),sep=""), plot = p, width = 8, height = 6)}
if(z==2){ggsave(paste("../MD",paste(j,"DE.pdf",sep=""),sep=""), plot = p, width = 8, height = 6)}
if(z==1){ggsave(paste("../LD",paste(j,"DE.png",sep=""),sep=""), plot = p, width = 8, height = 6)}
if(z==2){ggsave(paste("../MD",paste(j,"DE.png",sep=""),sep=""), plot = p, width = 8, height = 6)}
#ITERATED BACKWARD ERROR
p <- ggplot(data = dff, aes(x=itec, y=rc, color=kwc, group=enumc))+geom_line()+scale_y_continuous(limits=c(1E-16, 1),trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15))+ scale_x_continuous(limits=c(0, xlim))
p <- p+facet_grid(Sthetac ~ Sic+Swhichc+l)+labs(y="Iterated backward error")+labs(x="Solver iter.")+guides(color=guide_legend(title="k"))+scale_color_gradient(low = "blue", high = "red")
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
if(z==1){ggsave(paste("../LD",paste(j,"IBA.pdf",sep=""),sep=""), plot = p, width = 8, height = 6)}
if(z==2){ggsave(paste("../MD",paste(j,"IBA.pdf",sep=""),sep=""), plot = p, width = 8, height = 6)}
if(z==1){ggsave(paste("../LD",paste(j,"IBA.png",sep=""),sep=""), plot = p, width = 8, height = 6)}
if(z==2){ggsave(paste("../MD",paste(j,"IBA.png",sep=""),sep=""), plot = p, width = 8, height = 6)}
}
setwd("../../../../")
print(getwd())
}
library(ggplot2)
args <- commandArgs(trailingOnly = TRUE)
rep = paste("./",args[1],sep="")
V_REP = c(paste(rep,"/ITE/A0/1E2/WgivenFalse/",sep=""), paste(rep,"/ITE/A1/1E2/WgivenFalse/",sep=""),
paste(rep,"/ITE/A2/1E2/WgivenFalse/",sep=""), paste(rep,"/ITE/A3/1E2/WgivenFalse/",sep=""), paste(rep,"/ITE/A4/1E2/WgivenFalse/",sep=""),
paste(rep,"/ITE/A2/1E4/WgivenFalse/",sep=""), paste(rep,"/ITE/A3/1E4/WgivenFalse/",sep=""), paste(rep,"/ITE/A4/1E4/WgivenFalse/",sep=""),
paste(rep,"/ITE/A2/1E6/WgivenFalse/",sep=""), paste(rep,"/ITE/A3/1E6/WgivenFalse/",sep=""), paste(rep,"/ITE/A4/1E6/WgivenFalse/",sep=""))
ARG1 = c("S0", "S1", "S2", "S3")
ARG1 = c("S1", "S2", "S3")
ARG2 = c("LD", "MD", "LMD")
ARG3 = c(args[2])
V_FILE <- ""
V_OUT <- ""
V_PNG <- ""
list <- ""
t = 1; z = 1
for(z in 1:length(V_REP)){
t = 1
for(i in ARG1){
a = paste("*", i, sep="_")
a2 = i
if(i == "S3"){ARG22 = ARG2}
if(i != "S3"){ARG22 = c("LD")}
for(j in ARG22){
b = paste(a, j, sep="-")
b2 = paste(a2, j, sep="-")
for(k in ARG3){
c = paste(b, k, sep="-")
c2 = paste(b2, k, sep="-")
V_FILE[t] = paste(c, ".csv", sep="")
V_OUT[t] = paste(c2, ".csv", sep="")
V_PNG[t] = paste(c2, ".pdf", sep="")
t = t+1
system(command = paste(paste(paste("cat", V_REP[z],sep=" "), V_FILE[t-1], sep=""), paste(">",paste(V_REP[z],V_OUT[t-1],sep=""),sep=" "), sep=" "), intern=TRUE)
system(command = paste(paste("rm", V_REP[z],sep=" "), V_FILE[t-1], sep=""), intern = TRUE)
}
}
}
}
iliste = 0
for(IV_REP in V_REP){
ip = 0; iliste = iliste + 1; dall <- 0
for(IV_OUT in V_OUT){
ip = ip + 1
# Read Data:
filen = paste(IV_REP, IV_OUT, sep="")
df3 = read.csv(file = filen)
ite <- 0; ba <- 0.0; r <- 0.0; de <- 0.0; kw <- 0; s <- 0; enum <- 0; Si <-0; Swhich <- 0; Stheta <- 0; cpt <- 0
dimi1 = 1; dimi2 = length(df3$iterations); dimj = dimi2; dim = dimi2
ns = max(df3$enum)
z = 0; t = 1; i = 1; maxrun = 3; val = 0
for(k in seq(0, ns, by=maxrun)){
z = 0
ite <- 0; ba <- 0.0; r <- 0.0; de <- 0.0; kw <- 0; s <- 0; enum <- 0; Si <-0; Swhich <- 0; Stheta <- 0
for(i in dimi1:dimi2){
if( (df3$s[i] == k & k==0) | (df3$enum[i] %in% seq(k-maxrun+1,k,by=1)) ){
ite[z+1] = df3$iterations[i]; ba[z+1] = df3$backward_error[i]
de[z+1] = df3$direct_error[i]; r[z+1] = df3$approx_backward[i]
kw[z+1] = df3$k[i]; s[z+1] = df3$s[i]; enum[z+1] = val #df3$enum[i]
Si[z+1] = df3$Si[i]; Swhich[z+1] = df3$which[i]; Stheta[z+1] = df3$theta[i]
z = z + 1
}
}# end i
val = val + 1
z = length(s)
itec = 0 * seq(1, z, 1) - 1
bac = 0.0*itec; rc = bac; dec = bac; kwc = itec; sc = itec; enumc = itec
Sic = itec; Swhichc = itec; Sthetac = bac
cpt = 0*seq(1, z, 1)
# *Mean in Group(enum)
t = 1; i = 1
while(i <= z){
if((ite[i] %in% itec) == FALSE){
for(j in i:z){
if(ite[j] == ite[i]){
itec[t] = ite[i]
bac[t] = bac[t] + ba[j]; rc[t] = rc[t] + r[j]; dec[t] = dec[t] + de[j]
kwc[t] = kw[i]; sc[t] = s[i]; enumc[t] = enum[i]
Sic[t] = Si[i]; Swhichc[t] = Swhich[i]; Sthetac[t] = Stheta[i]
cpt[t] = cpt[t] + 1
}#end if
}#end j
t = t + 1
}#end if
i = i + 1
}#end while
ierr = 1
while(itec[ierr] > -1){
ierr = ierr + 1
}#end while
ierr = ierr - 1
# **Delate if -1
itec = itec[1:ierr]
bac = bac[1:ierr]; rc = rc[1:ierr]; dec = dec[1:ierr]
kwc = kwc[1:ierr]; sc = sc[1:ierr]; enumc = enumc[1:ierr]
Sic = Sic[1:ierr]; Swhichc = Swhichc[1:ierr]; Sthetac = Sthetac[1:ierr]
cpt = cpt[1:ierr]
# **Compute mean
bac = bac/cpt; rc = rc/cpt; dec = dec/cpt
# *Create data frame with group(s)
#dff <- data.frame(itec, bac, rc, dec, kwc, sc)
if(k == 0){
df = data.frame(itec, bac, rc, dec, kwc, sc, enumc, Sic, Swhichc, Sthetac)
}else{df = rbind(df, data.frame(itec, bac, rc, dec, kwc, sc, enumc, Sic, Swhichc, Sthetac))}
}#end k
#pname = paste(IV_REP, V_PNG[ip], sep="")
if(ip==1){
dall = df
}else{dall = rbind(dall, df)}
#print(dall)
#write.csv(df, filen)
list[iliste] = filen
#p <- ggplot(data = df, aes(x=itec, y=bac,color=kwc, group=enumc))+geom_line()+scale_color_gradient(low = "blue", high = "red")+scale_y_continuous(trans='log2')
#ggsave(pname, plot = p, width = 11, height = 8)
}
write.csv(dall, paste(IV_REP,"R.csv",sep=""))
}
library(ggplot2)
args <- commandArgs(trailingOnly = TRUE)
rep = paste("./",args[1],sep="")
V_REP = c(paste(rep,"/ITE/A2/1E2/WgivenFalse/",sep=""), paste(rep,"/ITE/A3/1E2/WgivenFalse/",sep=""), paste(rep,"/ITE/A4/1E2/WgivenFalse/",sep=""),
paste(rep,"/ITE/A2/1E4/WgivenFalse/",sep=""), paste(rep,"/ITE/A3/1E4/WgivenFalse/",sep=""), paste(rep,"/ITE/A4/1E4/WgivenFalse/",sep=""),
paste(rep,"/ITE/A2/1E6/WgivenFalse/",sep=""), paste(rep,"/ITE/A3/1E6/WgivenFalse/",sep=""), paste(rep,"/ITE/A4/1E6/WgivenFalse/",sep=""))
ARG1 = c("S3")
ARG2 = c("LD", "MD", "LMD")
ARG3 = c(args[2])
V_FILE <- ""
V_OUT <- ""
V_PNG <- ""
list <- ""
t = 1; z = 1
for(z in 1:length(V_REP)){
t = 1
for(i in ARG1){
a = paste("*", i, sep="_")
a2 = i
if(i == "S3"){ARG22 = ARG2}
if(i != "S3"){ARG22 = c("LD")}
for(j in ARG22){
b = paste(a, j, sep="-")
b2 = paste(a2, j, sep="-")
for(k in ARG3){
c = paste(b, k, sep="-")
c2 = paste(b2, k, sep="-")
V_FILE[t] = paste(c, ".csv", sep="")
V_OUT[t] = paste(c2, ".csv", sep="")
V_PNG[t] = paste(c2, ".pdf", sep="")
t = t+1
system(command = paste(paste(paste("cat", V_REP[z],sep=" "), V_FILE[t-1], sep=""), paste(">",paste(V_REP[z],V_OUT[t-1],sep=""),sep=" "), sep=" "), intern=TRUE)
system(command = paste(paste("rm", V_REP[z],sep=" "), V_FILE[t-1], sep=""), intern = TRUE)
}
}
}
}
iliste = 0
for(IV_REP in V_REP){
ip = 0; iliste = iliste + 1; dall <- 0
for(IV_OUT in V_OUT){
ip = ip + 1
# Read Data:
filen = paste(IV_REP, IV_OUT, sep="")
df3 = read.csv(file = filen)
ite <- 0; ba <- 0.0; r <- 0.0; de <- 0.0; kw <- 0; s <- 0; enum <- 0; Si <-0; Swhich <- 0; Stheta <- 0; cpt <- 0
dimi1 = 1; dimi2 = length(df3$iterations); dimj = dimi2; dim = dimi2
ns = max(df3$enum)
z = 0; t = 1; i = 1; maxrun = 3; val = 0
for(k in seq(0, ns, by=maxrun)){
z = 0
ite <- 0; ba <- 0.0; r <- 0.0; de <- 0.0; kw <- 0; s <- 0; enum <- 0; Si <-0; Swhich <- 0; Stheta <- 0
for(i in dimi1:dimi2){
if( (df3$s[i] == k & k==0) | (df3$enum[i] %in% seq(k-maxrun+1,k,by=1)) ){
ite[z+1] = df3$iterations[i]; ba[z+1] = df3$backward_error[i]
de[z+1] = df3$direct_error[i]; r[z+1] = df3$approx_backward[i]
kw[z+1] = df3$k[i]; s[z+1] = df3$s[i]; enum[z+1] = val #df3$enum[i]
Si[z+1] = df3$Si[i]; Swhich[z+1] = df3$which[i]; Stheta[z+1] = df3$theta[i]
z = z + 1
}
}# end i
val = val + 1
z = length(s)
itec = 0 * seq(1, z, 1) - 1
bac = 0.0*itec; rc = bac; dec = bac; kwc = itec; sc = itec; enumc = itec
Sic = itec; Swhichc = itec; Sthetac = bac
cpt = 0*seq(1, z, 1)
# *Mean in Group(enum)
t = 1; i = 1
while(i <= z){
if((ite[i] %in% itec) == FALSE){
for(j in i:z){
if(ite[j] == ite[i]){
itec[t] = ite[i]
bac[t] = bac[t] + ba[j]; rc[t] = rc[t] + r[j]; dec[t] = dec[t] + de[j]
kwc[t] = kw[i]; sc[t] = s[i]; enumc[t] = enum[i]
Sic[t] = Si[i]; Swhichc[t] = Swhich[i]; Sthetac[t] = Stheta[i]
cpt[t] = cpt[t] + 1
}#end if
}#end j
t = t + 1
}#end if
i = i + 1
}#end while
ierr = 1
while(itec[ierr] > -1){
ierr = ierr + 1
}#end while
ierr = ierr - 1
# **Delate if -1
itec = itec[1:ierr]
bac = bac[1:ierr]; rc = rc[1:ierr]; dec = dec[1:ierr]
kwc = kwc[1:ierr]; sc = sc[1:ierr]; enumc = enumc[1:ierr]
Sic = Sic[1:ierr]; Swhichc = Swhichc[1:ierr]; Sthetac = Sthetac[1:ierr]
cpt = cpt[1:ierr]
# **Compute mean
bac = bac/cpt; rc = rc/cpt; dec = dec/cpt
# *Create data frame with group(s)
#dff <- data.frame(itec, bac, rc, dec, kwc, sc)
if(k == 0){
df = data.frame(itec, bac, rc, dec, kwc, sc, enumc, Sic, Swhichc, Sthetac)
}else{df = rbind(df, data.frame(itec, bac, rc, dec, kwc, sc, enumc, Sic, Swhichc, Sthetac))}
}#end k
#pname = paste(IV_REP, V_PNG[ip], sep="")
if(ip==1){
dall = df
}else{dall = rbind(dall, df)}
#print(dall)
#write.csv(df, filen)
list[iliste] = filen
#p <- ggplot(data = df, aes(x=itec, y=bac,color=kwc, group=enumc))+geom_line()+scale_color_gradient(low = "blue", high = "red")+scale_y_continuous(trans='log2')
#ggsave(pname, plot = p, width = 11, height = 8)
}
write.csv(dall, paste(IV_REP,"R.csv",sep=""))
}
library(ggplot2)
args <- commandArgs(trailingOnly = TRUE)
rep = paste("./",args[1],sep="")
V_REP1 = c(paste(rep,"/ITE/A2/L/",sep=""), paste(rep,"/ITE/A3/L/",sep=""), paste(rep,"/ITE/A4/L/",sep=""))
V_REP = matrix(c(V_REP1), nrow=1, ncol=3, byrow=TRUE)
A_FILE = c("A2", "A3", "A4")
pp <- 0
print(getwd())
for(j in 1:3){
setwd(V_REP[1,j])
print(getwd())
liste = dir()
print(liste)
ns = length(liste)
for(ifn in 1:ns){
print(liste[ifn])
dff=read.csv(liste[ifn])
if(ifn == 1){df = read.csv(liste[ifn])}
else{df = rbind(df, dff)}
}
df <- subset(df, Sic != "S1")
df <- subset(df, Sic != "S2")
xlim = 0
for(ik in 1:length(df$kwc)){if(df$kwc[ik]==0){xlim =max(xlim,df$itec[ik])}}
xlim = xlim + 2
#df$lt <- 0
#for(ik in 1:length(df$kwc)){if(df$kwc[ik]==1){df$lt[ik]=1};if(df$kwc[ik]==1){df$lt[ik]=2};if(df$kwc[ik]==3){df$lt[ik]=3};if(df$kwc[ik]==4){df$lt[ik]=4}}
for(z in 1:2){
if(z==1){df1 <- subset(df, Swhichc != "LMD");df1 <- subset(df1, Swhichc != "MD");dff <- df1}
if(z==2){df2 <- subset(df, Swhichc != "LMD");df2 <- subset(df2, Swhichc != "LD");dff <- df2}
dff$s =dff$sc
#BACKWARD ERROR
p <- ggplot(data = dff, aes(x=itec, y=bac, color=s, group=enumc))+geom_line()+scale_y_continuous(limits=c(1E-16, 1),trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15))+ scale_x_continuous(limits=c(0, xlim))
p <- p+facet_grid(Sthetac ~ Sic+Swhichc+l)+labs(y="Backward error")+labs(x="Solver iter.")+scale_color_gradient(low = "blue", high = "red")
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
if(z==1){ggsave(paste("../LD",paste(j,"BA.pdf",sep=""),sep=""), plot = p, width = 8, height = 6)}
if(z==2){ggsave(paste("../MD",paste(j,"BA.pdf",sep=""),sep=""), plot = p, width = 8, height = 6)}
if(z==1){ggsave(paste("../LD",paste(j,"BA.png",sep=""),sep=""), plot = p, width = 8, height = 6)}
if(z==2){ggsave(paste("../MD",paste(j,"BA.png",sep=""),sep=""), plot = p, width = 8, height = 6)}
#DIRECT ERROR
p <- ggplot(data = dff, aes(x=itec, y=dec, color=s, group=enumc))+geom_line()+scale_y_continuous(limits=c(1E-16, 1),trans='log2', breaks=c(1,1E-4,1E-8,1E-13))+ scale_x_continuous(limits=c(0, xlim))
p <- p+facet_grid(Sthetac ~ Sic+Swhichc+l)+labs(y="Direct error (A-norm)")+labs(x="Solver iter.")+scale_color_gradient(low = "blue", high = "red")
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
if(z==1){ggsave(paste("../LD",paste(j,"DE.pdf",sep=""),sep=""), plot = p, width = 8, height = 6)}
if(z==2){ggsave(paste("../MD",paste(j,"DE.pdf",sep=""),sep=""), plot = p, width = 8, height = 6)}
if(z==1){ggsave(paste("../LD",paste(j,"DE.png",sep=""),sep=""), plot = p, width = 8, height = 6)}
if(z==2){ggsave(paste("../MD",paste(j,"DE.png",sep=""),sep=""), plot = p, width = 8, height = 6)}
#ITERATED BACKWARD ERROR
p <- ggplot(data = dff, aes(x=itec, y=rc, color=s, group=enumc))+geom_line()+scale_y_continuous(limits=c(1E-16, 1),trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15))+ scale_x_continuous(limits=c(0, xlim))
p <- p+facet_grid(Sthetac ~ Sic+Swhichc+l)+labs(y="Iterated backward error")+labs(x="Solver iter.")+scale_color_gradient(low = "blue", high = "red")
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
if(z==1){ggsave(paste("../LD",paste(j,"IBA.pdf",sep=""),sep=""), plot = p, width = 8, height = 6)}
if(z==2){ggsave(paste("../MD",paste(j,"IBA.pdf",sep=""),sep=""), plot = p, width = 8, height = 6)}
if(z==1){ggsave(paste("../LD",paste(j,"IBA.png",sep=""),sep=""), plot = p, width = 8, height = 6)}
if(z==2){ggsave(paste("../MD",paste(j,"IBA.png",sep=""),sep=""), plot = p, width = 8, height = 6)}
}
setwd("../../../../")
print(getwd())
}
library(ggplot2)
args <- commandArgs(trailingOnly = TRUE)
rep = paste("./",args[3],sep="")
rep2 = paste(paste("./",args[3],sep=""),"/ITE/",sep="")
outfn1 = paste(paste("L_",args[1],sep=""),args[2],sep="_")
V_REP1 = c(paste(rep,"/ITE/A2/1E2/WgivenFalse/",sep=""), paste(rep,"/ITE/A2/1E4/WgivenFalse/",sep=""), paste(rep,"/ITE/A2/1E6/WgivenFalse/",sep=""))
V_REP2 = c(paste(rep,"/ITE/A3/1E2/WgivenFalse/",sep=""), paste(rep,"/ITE/A3/1E4/WgivenFalse/",sep=""), paste(rep,"/ITE/A3/1E6/WgivenFalse/",sep=""))
V_REP3 = c(paste(rep,"/ITE/A4/1E2/WgivenFalse/",sep=""), paste(rep,"/ITE/A4/1E4/WgivenFalse/",sep=""), paste(rep,"/ITE/A4/1E6/WgivenFalse/",sep=""))
V_REP = matrix(c(V_REP1, V_REP2, V_REP3), nrow=3, ncol=3, byrow=TRUE)
A_FILE = c("A2", "A3", "A4")
dall <- 0
for(i in 1:3){
for(j in 1:3 ){
r = V_REP[i, j]; fn = paste(r, "R.csv",sep=""); dfff = read.csv(fn)
if(j==1){dall = read.csv(fn)}else{dall = rbind(dall, dfff)}
}
ik=1;for(ij in dall$Sic){if(ij==1){dall$Sic[ik]="S1"};if(ij==2){dall$Sic[ik]="S2"};if(ij==3){dall$Sic[ik]="S3"};ik = ik + 1}
ik=1;for(ij in dall$Swhichc){if(ij==1){dall$Swhichc[ik]="LD"};if(ij==2){dall$Swhichc[ik]="MD"};if(ij==3){dall$Swhichc[ik]="LMD"};ik = ik + 1}
ik=1;for(ij in dall$Sthetac){if(ij==1){dall$Sthetac[ik]="1E2"};if(ij==2){dall$Sthetac[ik]="1E4"};if(ij==3){dall$Sthetac[ik]="1E6"};ik = ik + 1}
outfn = paste(paste(paste("L_",args[1],sep=""),args[2],sep="_"),".csv",sep="")
print(outfn)
dall$l = paste(args[1],args[2],sep=" - ")
for(ik in 1:length(dall$l)){if(dall$Sic[ik]=="S1" | dall$Sic[ik]=="S2"){dall$l[ik]="-"}}
write.csv(dall, paste(paste(rep2,A_FILE[i],sep=""),outfn,sep="/L/"))
xlim = 0
for(ik in 1:length(dall$kwc)){if(dall$kwc[ik]==0){xlim =max(xlim,dall$itec[ik])}}
xlim = xlim + 2
dall$lt <- 0
for(ik in 1:length(dall$kwc)){if(dall$kwc[ik]==1){dall$lt[ik]=1};if(dall$kwc[ik]==2){dall$lt[ik]=2};if(dall$kwc[ik]==3){dall$lt[ik]=3};if(dall$kwc[ik]==4){dall$lt[ik]=4}}
dall$s = dall$sc
for(iz in 1:2){
#p <- ggplot(data = dall, aes(x=itec, y=bac, color=factor(kwc), group=enumc))+geom_line()+scale_y_continuous(trans='log2')
#p <- p+facet_grid(Sthetac ~ Sic + Swhichc)+labs(y="backward error")+labs(x="solver iter.")+guides(color=guide_legend(title="k"))
#BACKWARD ERROR:
p <- ggplot(data = dall, aes(x=itec, y=bac, color=s, group=enumc))+geom_line()+scale_y_continuous(limits=c(1E-16, 1),trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15))+ scale_x_continuous(limits=c(0, xlim))
p <- p+facet_grid(Sthetac ~ Sic + Swhichc+ l)+labs(y="Backward error")+labs(x="Solver iter.")+scale_color_gradient(low = "blue", high = "red")
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
if(iz==2){
r = paste(paste(paste(rep2,A_FILE[i],sep=""),paste("/ZOOM","BA",sep="_"),sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".pdf",sep=""),sep="_")
p <- p + scale_x_continuous(limits=c(0, 7))+ scale_y_continuous(trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15))
}else{r = paste(paste(paste(rep2,A_FILE[i],sep=""),"/BA",sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".pdf",sep=""),sep="_")}
ggsave(r, plot = p, width = 8, height = 6)
if(iz==2){
r = paste(paste(paste(rep2,A_FILE[i],sep=""),paste("/ZOOM","BA",sep="_"),sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".png",sep=""),sep="_")
p <- p + scale_x_continuous(limits=c(0, 7))+ scale_y_continuous(trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15))
}else{r = paste(paste(paste(rep2,A_FILE[i],sep=""),"/BA",sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".png",sep=""),sep="_")}
ggsave(r, plot = p, width = 8, height = 6)
#DIRECT ERROR (A-NORM):
p <- ggplot(data = dall, aes(x=itec, y=dec, color=s, group=enumc))+geom_line()+scale_y_continuous(trans='log2', breaks=c(1,1E-4,1E-8,1E-13))+ scale_x_continuous(limits=c(0, xlim))
p <- p+facet_grid(Sthetac ~ Sic + Swhichc + l)+labs(y="Direct error (A-norm)")+labs(x="Solver iter.")+scale_color_gradient(low = "blue", high = "red")
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
if(iz==2){
r = paste(paste(paste(rep2,A_FILE[i],sep=""),paste("/ZOOM","DE",sep="_"),sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".pdf",sep=""),sep="_")
p <- p + scale_x_continuous(limits=c(0, 5))+ scale_y_continuous(trans='log2', breaks=c(1,1E-4,1E-8,1E-13))
}else{r = paste(paste(paste(rep2,A_FILE[i],sep=""),"/DE",sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".pdf",sep=""),sep="_")}
ggsave(r, plot = p, width = 8, height = 6)
if(iz==2){
r = paste(paste(paste(rep2,A_FILE[i],sep=""),paste("/ZOOM","DE",sep="_"),sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".png",sep=""),sep="_")
p <- p + scale_x_continuous(limits=c(0, 5))+ scale_y_continuous(trans='log2', breaks=c(1,1E-4,1E-8,1E-13))
}else{r = paste(paste(paste(rep2,A_FILE[i],sep=""),"/DE",sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".png",sep=""),sep="_")}
ggsave(r, plot = p, width = 8, height = 6)
#ITERATED BACKWARD ERROR:
p <- ggplot(data = dall, aes(x=itec, y=rc, color=s, group=enumc))+geom_line()+scale_y_continuous(trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15))+ scale_x_continuous(limits=c(0, xlim))
p <- p+facet_grid(Sthetac ~ Sic + Swhichc + l)+labs(y="Iterated bacward error")+labs(x="Solver iter.")+scale_color_gradient(low = "blue", high = "red")
p <- p + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
if(iz==2){
r = paste(paste(paste(rep2,A_FILE[i],sep=""),paste("/ZOOM","IBA",sep="_"),sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".pdf",sep=""),sep="_")
p <- p + scale_x_continuous(limits=c(0, 5))+ scale_y_continuous(trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15))
}else{r = paste(paste(paste(rep2,A_FILE[i],sep=""),"/IBA",sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".pdf",sep=""),sep="_")}
ggsave(r, plot = p, width = 8, height = 6)
if(iz==2){
r = paste(paste(paste(rep2,A_FILE[i],sep=""),paste("/ZOOM","IBA",sep="_"),sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".png",sep=""),sep="_")
p <- p + scale_x_continuous(limits=c(0, 5))+ scale_y_continuous(trans='log2', breaks=c(1,1E-4,1E-8,1E-12,1E-15))
}else{r = paste(paste(paste(rep2,A_FILE[i],sep=""),"/IBA",sep=""),paste(paste(A_FILE[i],outfn1,sep="_"),".png",sep=""),sep="_")}
ggsave(r, plot = p, width = 8, height = 6)
}
}
import numpy as np import scipy.linalg as la import os
def csv_cg(ite, data, filen, tags=None): """Used to store in a csv format data during CG. data : iteration backward_error approx_backward direct_error W s Wgiven Si """ tags = ["A", "0", "N", "EXP", "1", "1", "0", "0"] if tags is None else tags header = 'iterations, backward_error, approx_backward, direct_error, k, s, enum, Si, which, theta' REP = "./"+tags[3]+"/ITE/"+tags[0]+"/"+tags[2]+"/Wgiven"+str(data[6])+"/" os.makedirs(REP, exist_ok=True) os.makedirs("./"+tags[3]+"/ITE/"+tags[0]+"/L/", exist_ok=True) nk = 0 if data[4] is None else data[4].shape[1] T1=0;T2=0 if(tags[6]=='S0'):T1=0 if(tags[6]=='S1'):T1=1 if(tags[6]=='S2'):T1=2 if(tags[6]=='S3'):T1=3 if(tags[7]=='LD'):T2=1 if(tags[7]=='MD'):T2=2 if(tags[7]=='LMD'):T2=3 if ite == 0 and data[5] == 1 and tags[1] == "1" and tags[4] == "1": F = open(REP+tags[1]+"_"+filen,"w",newline="\n"); strF = str(data[0])+", "+str(data[1])+", "+str(data[2])+", "+str(data[3])+", "+str(nk)+", "+str(data[5])+", "+tags[5]+", "+str(T1)+", "+str(T2)+", "+tags[2]; F.write(header+"\n"); F.write(strF+"\n"); F.close(); else: F = open(REP+tags[1]+"_"+filen,"a",newline="\n"); strF = str(data[0])+", "+str(data[1])+", "+str(data[2])+", "+str(data[3])+", "+str(nk)+", "+str(data[5])+", "+tags[5]+", "+str(T1)+", "+str(T2)+", "+tags[2]; F.write(strF+"\n"); F.close(); return None
def csv_eigen(A, W, TEval, TEvec, strat, tags=None): """Used to strore in csv format eigen residual Parameters: ----------- A: Matrix-like object, the linear operator W: Matrix-like object, Approximate eigenspace TEval: List, true eigen values of A TEvec: Matrix-like obkect, True eigen vectors of A strat: List, parameters of cg-def Tags: List, [str(Matrix-name),str(int)] Header: '# cols number of W corresponding one in TEvec s Rayleigh Quotient Eigen Residual Sin<(wi,TEveci) ' """ if W is None: return None tags = ["A", "0", "N", "EXP", "1", "1", "0", "0"] if tags is None else tags header = 'i, a, s, WRQ, WEQ, SinW, A, SI, part, L, k, theta' REP = "./"+tags[3]+"/EIGEN/"+tags[0]+"/"+tags[2]+"/Wgiven"+str(strat[0])+"/" os.makedirs(REP, exist_ok=True) filen = REP+strat[5] + "-" + strat[7]+"-"+"dense"+"-"+str(strat[9])+".csv" if strat[6] == 1: if tags[1] == "1": os.makedirs("./"+tags[3]+"/EIGEN/"+tags[0]+"/"+tags[2]+"/", exist_ok=True) F = open("./"+tags[3]+"/EIGEN/"+tags[0]+"/"+tags[2]+"/"+"TrueEigenVals_n_"+str(len(TEval))+".csv", "w", newline="\n") F.write("i, len, eval"+"\n") for i in range(len(TEval)): F.write(str(i)+", "+str(len(TEval))+", "+str(TEval[i])+"\n") F.close() if tags[1] == "1" and tags[6] == "S1" and tags[7] == "LD": F = open(filen, "w", newline="\n") F.write(header+"\n") F.close() F = open(filen, "a", newline="\n") for i in range(W.shape[1]): WRQ = np.asscalar( (W[...,i:i+1].T@A@W[...,i:i+1])/(W[...,i:i+1].T@W[...,i:i+1]) ) WEQ = la.norm(A@W[...,i:i+1]-WRQ*W[...,i:i+1]) S1 = 10.0; S2 = 0.0 ; a = 0 for k in range(TEvec.shape[1]): C1 = np.asscalar( ( TEvec[...,k:k+1].T@W[...,i:i+1] )/( la.norm(W[...,i:i+1])*la.norm(TEvec[...,k:k+1]) ) ) S2 = np.sqrt(np.abs(1-C1*C1)) if S2 <= S1: a = k; S1 = S2 F.write( str(i)+", "+str(a)+", "+str(strat[6])+", "+str(WRQ)+", "+str(WEQ)+", "+str(S1)+", "+tags[0]+", "+tags[6]+", "+tags[7]+", "+str(strat[2])+"-"+str(strat[3])+", "+str(W.shape[1])+", "+tags[2]+"\n" ) F.close() return None
def lanczos_RR(A, teve, tever, nk, Aname, Stheta, L=None, b=None, maxiter=None, reortho=None): """Compute and return Vm and Tm such that Tm = Vm.T@A@Vm via Lanczos iteration. A: a matrix object-like b: a vector object-like randomly set if is None maxiter: integer, the number of iterations to perform, set to A.shape[0] if is None reortho: default false, if True compute a reorthogonalization of Vm This function return Vm,Tm """ header = 'iterations, WQRR, WERR, WSRR, WQHR, WEHR, WSHR, A, theta' header = 'iterations, VAL, WVAL, RHR, A, theta' REP = "./LANCZOS/" if L is None else "./LANCZOS/L/" os.makedirs(REP, exist_ok=True) #WVAL = ["WQRR", "WERR", "WSRR", "WQHR", "WEHR", "WSHR"] WVAL = ["Q", "E", "S", "Q", "E", "S"] WVAL2 = [3, 3, 3, 4, 4, 4] maxiter = A.shape[1] if maxiter is None else maxiter b = A@np.random.rand(1,maxiter).T if b is None else b reortho = False if reortho is None else reortho filen = REP+"n_"+str(A.shape[1])+"_k_"+str(nk)+"_" if reortho == False: filen = filen+"NoReortho"+".csv" else: filen = filen+"Reortho"+".csv" if Aname == "A0" and Stheta == "1E2": F = open(filen, "w", newline="\n") F.write(header+"\n") else: F = open(filen, "a", newline="\n") m = maxiter Ta = np.zeros(m) Tb = np.zeros(m) v = b / la.norm(b) Vm = np.vstack((v)) r = A @ v alpha = v.T@r Ta[0] = alpha r = r - alpha * v beta = la.norm(r) Tm = ( np.diag(Ta[0:1], 0) ) #Compute Eigenvalues lT and Eigenvectors of Tm: lT, gT = la.eigh(Tm) #Compute Eigenpair (RR): RR = Vm @ gT WQRR = np.asscalar( (RR.T@A@RR)/(RR.T@RR) ) WERR = la.norm(A @ RR - WQRR * RR) WSRR = 11.25 for ik in range(tever.shape[1]): CC = np.asscalar((tever[...,ik:ik+1].T@RR)/(la.norm(RR)*la.norm(tever[...,ik:ik+1]))) WSRR = min(np.sqrt(np.abs(1.0-CC*CC)), WSRR) #Solve Harmonic-Rayleigh Ritz Eigenvalue Pb: G = ((A@Vm).T)@(A@Vm) Lambda, tildW = la.eigh(a=G, b=Tm, lower=True, eigvals_only=False, overwrite_a=False, overwrite_b=False, type=1, check_finite=True) #Compute Eigenpair (HR): HR = Vm @ tildW WQHR = np.asscalar( (HR.T@A@HR)/(HR.T@HR) ) WEHR = la.norm(A @ HR - WQHR * HR) WSHR = 11.25 for ik in range(tever.shape[1]): CC = np.asscalar((tever[...,ik:ik+1].T@HR)/(la.norm(HR)*la.norm(tever[...,ik:ik+1]))) WSHR = min(np.sqrt(np.abs(1.0-CC*CC)), WSHR) #F.write(str(1)+", "+str(WQRR)+", "+str(WERR)+", "+str(WSRR)+", "+str(WQHR)+", "+str(WEHR)+", "+str(WSHR)+", "+Aname+", "+Stheta+"\n") VAL = [WQRR,WERR,WSRR,WQHR,WEHR, WSHR] for I_WVAL in range(len(WVAL)): F.write(str(1)+", "+str(VAL[I_WVAL])+", "+WVAL[I_WVAL]+", "+str(WVAL2[I_WVAL])+", "+Aname+", "+Stheta+"\n") for j in range(2, m+1, 1): Tb[j-2] = beta w = v v = r / beta Vm = np.hstack((Vm,v)) r = A @ v - beta * w alpha = v.T @ r Ta[j-1] = alpha Tm = (np.diag(Ta[0:j],0) + np.diag(Tb[0:j-1],-1) + np.diag(Tb[0:j-1],+1)) r = r - alpha * v if reortho == True: r = r - Vm @ (Vm.T @ r) beta = la.norm(r) if L is None: #Compute Eigenvalues lT and Eigenvectors of Tm: lT, gT = la.eigh(Tm) #Solve Harmonic-Rayleigh Ritz Eigenvalue Pb: G = ((A@Vm).T)@(A@Vm) Lambda, tildW = la.eigh(a=G, b=Tm, lower=True, eigvals_only=False, overwrite_a=False, overwrite_b=False, type=1, check_finite=True) for im in range(Tm.shape[1]): #Compute Eigenpair (RR): RR = Vm @ gT[...,im:im+1] WQRR = np.asscalar( (RR.T@A@RR)/(RR.T@RR) ) WERR = la.norm(A @ RR - WQRR * RR) WSRR = 11.25 for ik in range(tever.shape[1]): CC = np.asscalar((tever[...,ik:ik+1].T@RR)/(la.norm(RR)*la.norm(tever[...,ik:ik+1]))) WSRR = min(np.sqrt(np.abs(1.0-CC*CC)), WSRR) if WSRR == 0: WSRR = 1e-8 #Compute Eigenpair (HR): HR = Vm @ tildW[...,im:im+1] WQHR = np.asscalar( (HR.T@A@HR)/(HR.T@HR) ) WEHR = la.norm(A @ HR - WQHR * HR) WSHR = 11.25 for ik in range(tever.shape[1]): CC = np.asscalar((tever[...,ik:ik+1].T@HR)/(la.norm(HR)*la.norm(tever[...,ik:ik+1]))) WSHR = min(np.sqrt(np.abs(1.0-CC*CC)), WSHR) if WSHR == 0: WSHR = 1e-8 #F.write(str(j)+", "+str(WQRR)+", "+str(WERR)+", "+str(WSRR)+", "+str(WQHR)+", "+str(WEHR)+", "+str(WSHR)+", "+Aname+", "+Stheta+"\n") VAL = [WQRR,WERR,WSRR,WQHR,WEHR, WSHR] for I_WVAL in range(len(WVAL)): F.write(str(j)+", "+str(VAL[I_WVAL])+", "+WVAL[I_WVAL]+", "+str(WVAL2[I_WVAL])+", "+Aname+", "+Stheta+"\n") elif j%L==0: T = Tm.copy(); V = Vm.copy() Tm = Tm[j-L:j,j-L:j]; Vm = Vm[...,j-L:j] #Compute Eigenvalues lT and Eigenvectors of Tm: lT, gT = la.eigh(Tm) #Solve Harmonic-Rayleigh Ritz Eigenvalue Pb: G = ((A@Vm).T)@(A@Vm) Lambda, tildW = la.eigh(a=G, b=Tm, lower=True, eigvals_only=False, overwrite_a=False, overwrite_b=False, type=1, check_finite=True) for im in range(Tm.shape[1]): #Compute Eigenpair (RR): RR = Vm @ gT[...,im:im+1] WQRR = np.asscalar( (RR.T@A@RR)/(RR.T@RR) ) WERR = la.norm(A @ RR - WQRR * RR) WSRR = 11.25 for ik in range(tever.shape[1]): CC = np.asscalar((tever[...,ik:ik+1].T@RR)/(la.norm(RR)*la.norm(tever[...,ik:ik+1]))) WSRR = min(np.sqrt(np.abs(1.0-CC*CC)), WSRR) if WSRR == 0: WSRR = 1e-8 #Compute Eigenpair (HR): HR = Vm @ tildW[...,im:im+1] WQHR = np.asscalar( (HR.T@A@HR)/(HR.T@HR) ) WEHR = la.norm(A @ HR - WQHR * HR) WSHR = 11.25 for ik in range(tever.shape[1]): CC = np.asscalar((tever[...,ik:ik+1].T@HR)/(la.norm(HR)*la.norm(tever[...,ik:ik+1]))) WSHR = min(np.sqrt(np.abs(1.0-CC*CC)), WSHR) if WSHR == 0: WSHR = 1e-8 #F.write(str(j)+", "+str(WQRR)+", "+str(WERR)+", "+str(WSRR)+", "+str(WQHR)+", "+str(WEHR)+", "+str(WSHR)+", "+Aname+", "+Stheta+"\n") VAL = [WQRR,WERR,WSRR,WQHR,WEHR, WSHR] for I_WVAL in range(len(WVAL)): F.write(str(j)+", "+str(VAL[I_WVAL])+", "+WVAL[I_WVAL]+", "+str(WVAL2[I_WVAL])+", "+Aname+", "+Stheta+"\n") Tm = T.copy(); Vm = V.copy() if beta == 0.0: for ida in range(len(teve)): F.write(str(j+2)+", "+str(teve[ida])+", "+WVAL[0]+", "+str(20)+", "+Aname+", "+Stheta+"\n") F.close() return Tm, Vm, filen for ida in range(len(teve)): F.write(str(j+2)+", "+str(teve[ida])+", "+WVAL[0]+", "+str(20)+", "+Aname+", "+Stheta+"\n") F.close() return Tm, Vm, filen
library(ggplot2)
args <- commandArgs(trailingOnly = TRUE)
fn = args[1]
#df = read.csv(fn)
dff = read.csv(fn)
#A_FILE = c(" A0", " A1", " A2", " A3", " A4")
#for(AF in A_FILE){
# file = paste(fn, paste("_facet_",AF, sep=""), sep="_")
# dA <- subset(df, A==AF)
# #dA <- subset(dA, theta==1E6)
# p1 <- ggplot(data = dA, aes(x=iterations, y=VAL, color=VAL, group=RHR))+geom_point(shape=dA$RHR)#+scale_y_continuous(trans='log')
# p1 <- p1 + facet_grid(theta ~ WVAL)+ labs(y=" ")+labs(x="Lanczos iter.")+scale_color_gradient(low = "blue", high = "red")
# p1 <- p1 + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
# r = paste(file,".pdf",sep="")
# ggsave(r, plot = p1, width = 11, height = 8)
#}
#A_FILE = c(" A0", " A1", " A2", " A3", " A4")
#ARG = c(" Q", " E", " S")
#for(AF in A_FILE){
# file = paste(fn, AF, sep="_")
# dA <- subset(dff, A==AF)
# print(head(dA))
# for(i in ARG){
# file2 = paste(file, i, sep="_")
# dAA <- subset(dA, WVAL==i)
# p1 <- ggplot(data = dAA, aes(x=iterations, y=VAL, color=VAL, group=RHR))+geom_point(shape=dAA$RHR)
# if((AF != " A0" & AF != " A1") | ((AF == " A0" | AF == " A1") & (i != " Q")) ){p1 <- p1 + scale_y_continuous(trans='log')}
# p1 <- p1 + facet_grid(~ theta)+ labs(y=" ")+labs(x="Lanczos iter.")+scale_color_gradient(low = "blue", high = "red")#+ scale_y_continuous(limits=c(ymin, ymax))
# p1 <- p1 + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
# r = paste(file2,".pdf",sep="")
# ggsave(r, plot = p1, width = 11, height = 8)
# }
#}
A_FILE = c(" A0", " A1", " A2", " A3", " A4")
ARG = c(" Q", " E", " S")
RHR = c("RR", "HR")
for(AF in A_FILE){
file = paste(fn, AF, sep="_")
dA <- subset(dff, A==AF)
if(AF==" A0" | AF==" A1"){dA <- subset(dA, theta==1E2)}
dA$Eigenvalue_Cluster <- 10
n = length(dA$VAL)
if(AF==" A0"){c1 = 5; c2 = max(dA$iterations)-5}
else if(AF==" A1"){c1 = 0.25; c2 = 1.25}
else{c1 = 0.5; c2 = 1.50}
for(ic in seq(1, n, by=1)){
if(dA$WVAL[ic]==" Q"){
if(dA$VAL[ic]>=c1 & dA$VAL[ic]<=c2){dA$Eigenvalue_Cluster[ic]="Centered"}
else if(dA$VAL[ic]<c1){dA$Eigenvalue_Cluster[ic]="Least dominant"}
else{dA$Eigenvalue_Cluster[ic]="Most dominant"}
#dA$color[ic]=dA$VAL[ic]
#;dA$color[ic+1]=dA$VAL[ic];dA$color[ic+2]=dA$VAL[ic]
#dA$color[ic]=1;dA$color[ic+1]=1;dA$color[ic+2]=1
}
else if(dA$WVAL[ic]==" E"){
dA$Eigenvalue_Cluster[ic]=dA$Eigenvalue_Cluster[ic-1]
#dA$color[ic]=dA$VAL[ic-1]
}
else if(dA$WVAL[ic]==" S"){
dA$Eigenvalue_Cluster[ic]=dA$Eigenvalue_Cluster[ic-2]
#dA$color[ic]=dA$VAL[ic-2]
}
}
YTAGS=c("Rayleigh quotient", "Eigen residual", "Sin<(wi, ui)")
it = 0
for(i in ARG){
it = it + 1
file2 = paste(file, i, sep="_")
print(file2)
dAA <- subset(dA, WVAL==i)
for(k in 3:4){
dAAA <- subset(dAA, RHR == k)
if(i==" Q"){dAAA = rbind(dAAA, subset(dAA, RHR == 20))}
file3 = paste(file2, RHR[k-2], sep="_")
write.csv(dAAA, paste(file3,".csv",sep=""))
p1 <- ggplot(data = dAAA, aes(x=iterations, y=VAL, color=Eigenvalue_Cluster, group=RHR))+geom_point(shape=dAAA$RHR)
if((AF != " A0" & AF != " A1") | ((AF == " A0" | AF == " A1") & (i != " Q")) ){p1 <- p1 + scale_y_continuous(trans='log')}
p1 <- p1 + labs(y=YTAGS[it])+labs(x="Lanczos iter.")
if(AF!=" A0" & AF!=" A1"){p1 <- p1 + facet_grid(~ theta)}
#if(i == " S"){p1 <- p1 + ylim(max(c(min(dAAA$VAL),1e-8)), min(c(max(dAAA$VAL),1e8)))}
#p1 <- p1 + scale_color_gradient(low = "blue", high = "red")#+ scale_y_continuous(limits=c(ymin, ymax))
#p1 <- p1 + theme(legend.position="bottom",legend.box="horizontal",legend.direction="horizontal")+theme(text = element_text(size = 10))+theme(axis.text.x = element_text(angle = 45, hjust = 1))
r = paste(file3,".pdf",sep="")
ggsave(r, plot = p1, width = 11, height = 8)
}
}
}
1 Test Cases
Will performed test cases in order to determine which part of the spectrum need to be deflated to accelerate convergence. Multiple strategies to extract and compress eigenvectors approximations will be set up and convergence of the method will be study under the quality of those approximations. We first show how to build matrices with a desired spectrum and them we explain how to build the deflation subspace.
1.1 Matrices
The \textbf{Least Dominant} (LD) refers to the eigenpair having the smallest eigenvalue and \textbf{Most Dominant} refers to the eigenpair having the largest eigenvalue from the matrix \(A\). If the spectrum of \(A\) has some clusters, then \(\theta\) will refer to the distance between them and \textbf(Least Dominant Cluster} (LDC) is the cluster contening the smallest Eigenvalues close to \(\frac{1}{\theta}\) \textbf{Most Dominant Cluster} (MDC) is the one containing the Largest Eigenvlaues close to \(\theta\). A \textbf{Center Cluster} (CC) is used to define Eigenvalues between \(a\) and \(b\). To Construct matrices having a specific spectrum will use the formation based on EigenDecompositon: \(A=Q^{T}.\Sigma.Q\) with \(\Sigma\) a diagonal matrix such that \(A.Q^{T}_{j}=Q^{T}_{j}.\Sigma_{jj}\). Will construct five matrices \(A_{i=0,4}\) such that \(A_{0}\) has \(n\) eigen values from \(1\) to \(n\) linearly spaced by \(1\), \(A_{1}\) has \(n\) eigen values randomly pick into \([0.5,1.5]\). \(A_{2}\) has \(k\) eigen values randomly pick into \([0.5,1.5]*\theta^{-1}\) and \(n-k\) eigen values into \([0.5,1.5]\) randomly pick. \(A_{3}\) has \(k\) eigen values randomly pick into \([0.5,1.5]*\theta\) and \(n-k\) eigen values into \([0.5,1.5]\) randomly pick. $A4$= has \(k\) eigen values randomly pick into \([0.5,1.5]*\theta^{-1}\) and \(k\) eigen values randomly pick into \([0.5,1.5]*\theta\) and \(n-2*k\) eigen values into \([0.5,1.5]\) randomly pick. To construct the Eigen space \(Q^{T}\) will perform a \(=QR=\) decomposition onto a random matrice \(B\).
def setmatrix4exp(which="A0", n=40, k=1, theta=1e4): """ Set some matrices A, DA and Q such that A= Q.T@DA@Q to be used during deflation experimentations. INPUT: ------ which: str('[A0|A1|A2|A3|A4]') indicates which matrix to construct over: A0 --> diag(1,2...,N) a diagonal matrix with integer entries A1 --> a diagonal matrix with entries randomly (normal law) pick up in [0.5,1.5] A2 --> AL+A1 where AL+ means that k of A1 entries are randomly pick up in (1.0/theta)*[1/2,3/2] A3 --> A1+AR where AR+ means that k percent of A1 entries are randomly pick up in theta*[1/2,3/2] A4 --> A2 and A3 effects over A1 default is A0 n: the row number of the A matrix, default is 40 theta: parameter to set AL (1/theta) and AR (theta), default is 1e4 k: the number of elements of A to set in AL and AR, default is 1. OUTPUT: ------- A matrix-like object, DA vector-like object, Q matrix-like object """ #Compute DA the diagonal matrix involved in A = Q.T@DA@Q: if which== "A0" : R = np.arange(1,n+1,1); DA = 0*np.arange(1,n+1,1) for i in range(n): alpha = np.random.randint(low=0, high=len(R), size=None) DA[i] = R[alpha] RR = 0*np.arange(1,n-i,1) jj = 0 for j in range(len(R)): if j != alpha : RR[jj] = R[j]; jj+=1 R = RR DA = np.diag(DA,0) elif which == "A1" or "A2" or "A3" or "A4": DA = 1.0*np.random.sample((1,n))+0.5 DA = DA[0,0:n] if which == "A2": AL = (1.0/theta)*np.random.sample((1,k))+(1.0/(2.0*theta)) AL = AL[0,0:k] for i in range(k): DA[i] = AL[i] elif which == "A3": AR = theta*np.random.sample((1,k))+(0.5*theta) AR = AR[0,0:k] for j in range(k): DA[(n-1)-j] = AR[j] elif which == "A4": if(k%2==0): hk = int(k/2) else: hk = int((k+1)/2) AL = (1.0/theta)*np.random.sample((1,hk))+(1.0/(2.0*theta)) AL = AL[0,0:hk] AR = theta*np.random.sample((1,hk))+(0.5*theta) AR = AR[0,0:hk] for i in range(hk): DA[i] = AL[i];DA[(n-1)-i] = AR[i] DA = np.diag(DA,0) #Compute Q involved in A = Q.T@DA@Q like a Q-matrix from the QR decomposition # of a random matrix B: B = 1.0*np.random.rand(n,n).T Q,R = la.qr(B,mode='full') del R; del B #Finalized by setting A=Q.T@DA@Q #Return A Matrix like object, DA vector like object, Q.T matrix like object DA = np.hstack((DA.diagonal())) #Sort DA, Q.T in ascending order on DA: i=0 QT = Q.T.copy() while i<len(DA)-1: for j in range(i+1, len(DA)): if DA[j]<DA[i]: sw=DA[j]; DA[j]=DA[i]; DA[i]=sw vw=QT[...,j:j+1].copy(); QT[...,j:j+1]=QT[...,i:i+1].copy() QT[...,i:i+1]=vw.copy() i+=1 del Q return QT@np.diag(DA,0)@QT.T,DA,QT
1.2 Deflation Subspace
In practise the Deflated subspace W is build using approximate eigenvectors and
search vectors from the previous solve. But for this study will used also an existing
W. Here we solve \(A.x^{s}=b^{s},\, \for s=1,2\) and the first system (s=1) is solved
by CG and \(l\) search vectors are stored in order to compute \(k\) approximate eigen vectors \(w_{i=1,..,k}\)
such that during the second solve (s=2) the deflated subspace is spanned by \(W=span(w_{1},...,w_{k})\).
The part of the spectrum to be deflated is specified by arguments:
\begin{itemize} \item $LD-->MD$: $W_{i=1,k}$ are the $k$ eigen vectors associated to the $k$ smallest eigen values \item $MD-->LD$: $W_{i=1,k}$ are the $k$ eigen vectors associated to the $k$ Largest eigen values \item $LMD$: Combine the two arguments \end{itemize}
In this case will study the convergence of DEF-CG when \(l\), \(k\) vary and define the \(S_{i=0,1,2,3}\) test cases:
The stopping criterion used for those tests is based on the iterated residual such that DEF-CG stop when \(\frac{||r_{i}||_{2}}{||b||_{2}}<1E-16\).