# ________________create the final matrix with a LOOP ____________
iter <- nrow(t)
Mfinal <- matrix(NA, nrow=iter, ncol=m+4)
#use for 1D and 2D
for (i in 1:iter){
if (m==1){ Mfinal[i,] <- c(t[i,m], mu[i], sigma[i], lower[i], upper[i])}
if (m==2){ Mfinal[i,] <- c(t[i,m-1],t[i,m], mu[i], sigma[i])} #, lower[i], upper[i])}
}
#-------------PLOT IF NEEDED-----------------------------
a=0
if (a == 1){
xx <- seq(1, length(t), 1)
plot(XDoE, YDoE, type = "p", ylim=c(min(lower),max(upper)))
lines (Mfinal[ , 1], Mfinal[ , 2], type = "l")
lines(Mfinal[ , 1], Mfinal[ , 4], type = "l", col='red')
lines(Mfinal[ , 1], Mfinal[ , 5], type = "l", col='red')
}
#____________________ SAVE in EXCEL to use in the MATLAB _________________
#Mfinal
#Mfinal <- m[ ]
#write.xlsx(Mfinal, "14KrigingR_dataresults.xls", sheetName="Sheet1")
#print to CSV
write.table(Mfinal, file = "all_predictor_data", col.names = FALSE, row.names = FALSE, sep=",")
#close sink
## reset message sink and close the file connection
#sink(type="message")
#close(zz)
#clean space of variables to avoid problems
rm(list = ls())
rm(.Random.seed, envir=globalenv())
closeAllConnections()
.libPaths("C:/Users/ruidr/Documents/R/win-library/3.4")
.libPaths()
#install.packages("XLConnect", "C:/Users/Rui Teixeira/Documents/R/win-library/3.4")
rm(list=ls())
#____load all packages_____
#
library(DiceKriging)
library(XLConnect)
library(xlsx)
#
#capture all error messages to a file in order to evaluate in MATLAB
#zz <- file("all.Rout", open="wt")
#sink(zz, type="message")
# Half the Inputs from Matlab
# These are all included in an Excel file
#________________________________________________________________
#read current directory from txt file
# do not forget to change the pat here -------------------------  INPUT !!!!!
setwd("C:\\Users\\ruidr\\Google Drive\\Trinity Sync\\PhD\\021_DICE")
getwd()
#________________________________________________________________
#load data from Worksheet of Matlab
#
workbookfile1 <- XLConnect::loadWorkbook("01data_input_R.xls")
alldata <- readWorksheet(workbookfile1, sheet = "Sheet1");
alldata
# assign the variables
n = ncol(alldata)
len = nrow(alldata)
#
XDoE <- alldata[ , c(1:(n-2))]
YDoE <- alldata[ , c(n-1)]
sigmaDoE<- alldata[ , c(n)]
#alocate variables
x <- (XDoE)
x <- mapply(x, FUN=as.numeric) #transform in numeric
y <- c(YDoE)
y <- mapply(y, FUN=as.numeric) #transform in numeric
#------------------------------------------------
#___________________ BUILD THE KRIGING _______________________________________________
# Heteregeneous noise variances
sigmaDoE <- mapply(sigmaDoE, FUN=as.numeric)
noise.var <- (sigmaDoE^2)
#------------------------------------------------
workbookfile2 <- XLConnect::loadWorkbook("02control_data_R.xls")
controldata <- readWorksheet(workbookfile2, sheet = "Sheet1");
workbookfile3 <- XLConnect::loadWorkbook("03control_convergence.xls")
bounds <- readWorksheet(workbookfile3, sheet = "Sheet1");
#DEGREE
degreeindex <- controldata[ (1), ]
degreeindex <- mapply(degreeindex, FUN=as.numeric)
#
#degreeindex <- 5
############### select degree
# trend (beware to the usual I operator) ------ First for the case of 1D
if ((n-2)==1){if (degreeindex == 0){formula <- ~1}
if (degreeindex == 1){formula <- ~x}
if (degreeindex == 2){formula <-~polym(x, degree=2, raw=TRUE)}
if (degreeindex == 3){formula <- ~polym(x, degree=3, raw=TRUE)}
if (degreeindex == 4){formula <- ~polym(x, degree=4, raw=TRUE) }
if (degreeindex == 5){formula <- ~polym(x, degree=5, raw=TRUE)}
if (degreeindex == 6){formula <- ~polym(x, degree=6, raw=TRUE)}
}
#~x+I(x^2)+I(x^3)+I(x^4)+I(x^5)+I(x^6)}
############### select degree
# First for the case of 2D
if ((n-2)==2){if (degreeindex == 0){formula <- ~1}
if (degreeindex == 1){formula <- ~polym(x1,x2, degree=1, raw=TRUE)}
if (degreeindex == 2){formula <- ~polym(x1,x2, degree=2, raw=TRUE)}
if (degreeindex == 3){formula <- ~polym(x1,x2, degree=3, raw=TRUE)}
if (degreeindex == 4){formula <- ~polym(x1,x2, degree=4, raw=TRUE) }
if (degreeindex == 5){formula <- ~polym(x1,x2, degree=5, raw=TRUE)}
if (degreeindex == 6){formula <- ~polym(x1,x2, degree=6, raw=TRUE)}
}
#x1+x2+I(x1*x2)+I(x1^2)+I(x2^2)}
#formula1 <- formula1~polym(x1, x2, degree=2, raw=TRUE) # is equivalent to
#formula <- ~x1 + x2 + I(x1^2) + I(x2^2) + x1:x2
#CORRELATION
elmcovtype <- controldata[ (2), ]
#BOUNDS
upperlim <- (bounds[ (2), 1:(n-2)])
upperlim <- mapply(upperlim, FUN=as.numeric)
lowerlim <- (bounds[ (1), 1:(n-2)])
lowerlim <- mapply(lowerlim, FUN=as.numeric)
# Create the Kriging model
# check manual to see the different options to fill in !!!
#is.numeric(data.frame(y=y))
#is.numeric(data.frame(x=x))
#rename the cpolumns in data frame --- very important
if ((n-2) == 1){x <- data.frame (x=x)}
if ((n-2) == 2){
colnames(x) <- c("x1","x2")}
#VV = matrix(NA, nrow=len, ncol=1)
#for (i in 1:len) {VV[i,1] = var[i]}
#noise.var = data.frame(s2=VV)
#}
ini_estimator <- sample(50:300, 1) #esseential for convergence for now!!!
#myKriging <- km(formula=formula, design = x, response = data.frame(y=y),
#               covtype=elmcovtype, noise.var=noise.var, estim.method="MLE",
#              optim.method = "BFGS", control=list(maxit=100, popsize=ini_estimator),
#             lower = lowerlim, upper = upperlim, multistart = 1)
myKriging <- km(formula=formula, design = x, response = data.frame(y=y),
covtype=elmcovtype, estim.method="MLE",
optim.method = "BFGS", control=list(maxit=100, popsize=ini_estimator),
lower = lowerlim, upper = upperlim, multistart = 1)
# Be carefull with the optimization boundaries
#"BFGS", "gen" , "L-BFGS-B"
#control=list(maxit=500,trace=T,popsize=50)
#get the theta, sigma and beta parameters from the Kriging.
theta.est <- covparam2vect(myKriging@covariance)
sigma2.est <- myKriging@covariance@sd2
beta.est <- myKriging@trend.coef
# save these in a file - a file for each
write.xlsx(theta.est, "11theta_res.xls", sheetName="Sheet1")
write.xlsx(sigma2.est, "12sigma2_res.xls", sheetName="Sheet1")
write.xlsx(beta.est, "13beta_res.xls", sheetName="Sheet1")
#------------------------------------------------
#covtype="matern5_2" or "gauss"
#_____________________________________________________________
# Other half of the inputs from the matlab
# prediction
library(XLConnect)
workbookfile4 <- XLConnect::loadWorkbook("04Points_DoE_predictions.xls")
predictionpoints <- readWorksheet(workbookfile4, sheet = "Sheet1");
m = ncol(predictionpoints)
t <- predictionpoints[ , c(1:m)]
if (m == 1){}
if (m == 2){ t = expand.grid( t[ , c(1)], t[ , c(2)] ) }
#t <- seq(min(x), max(x), by=0.1)
t <- mapply(t, FUN=as.numeric)
#rename the cpolumns in data frame --- very important for 2D
if (m == 1){}
if (m == 2){colnames(t) <- c("x1","x2")}
#--------------------------------------------
#---- PREDICCTIONS --------------------------
#
#introduce an excel file that will control the prediction points
# the prediction matrix needs to be a data frame!!!! Very important
if (m == 1){t <- data.frame (x=t)}
if (m == 2){t <- data.frame(t)}
p <- predict(myKriging, newdata=t, type="UK", se.compute=TRUE)
mu <- p$mean
sigma <- p$sd
# Intervals
#lower <- p$lower95; upper <- p$upper95
# ________________create the final matrix with a LOOP ____________
iter <- nrow(t)
Mfinal <- matrix(NA, nrow=iter, ncol=m+4)
#use for 1D and 2D
for (i in 1:iter){
if (m==1){ Mfinal[i,] <- c(t[i,m], mu[i], sigma[i], lower[i], upper[i])}
if (m==2){ Mfinal[i,] <- c(t[i,m-1],t[i,m], mu[i], sigma[i])} #, lower[i], upper[i])}
}
#-------------PLOT IF NEEDED-----------------------------
a=0
if (a == 1){
xx <- seq(1, length(t), 1)
plot(XDoE, YDoE, type = "p", ylim=c(min(lower),max(upper)))
lines (Mfinal[ , 1], Mfinal[ , 2], type = "l")
lines(Mfinal[ , 1], Mfinal[ , 4], type = "l", col='red')
lines(Mfinal[ , 1], Mfinal[ , 5], type = "l", col='red')
}
#____________________ SAVE in EXCEL to use in the MATLAB _________________
#Mfinal
#Mfinal <- m[ ]
#write.xlsx(Mfinal, "14KrigingR_dataresults.xls", sheetName="Sheet1")
#print to CSV
write.table(Mfinal, file = "all_predictor_data", col.names = FALSE, row.names = FALSE, sep=",")
#close sink
## reset message sink and close the file connection
#sink(type="message")
#close(zz)
y
#clean space of variables to avoid problems
rm(list = ls())
rm(.Random.seed, envir=globalenv())
closeAllConnections()
.libPaths("C:/Users/ruidr/Documents/R/win-library/3.4")
.libPaths()
#install.packages("XLConnect", "C:/Users/Rui Teixeira/Documents/R/win-library/3.4")
rm(list=ls())
#____load all packages_____
#
library(DiceKriging)
library(XLConnect)
library(xlsx)
#
#capture all error messages to a file in order to evaluate in MATLAB
#zz <- file("all.Rout", open="wt")
#sink(zz, type="message")
# Half the Inputs from Matlab
# These are all included in an Excel file
#________________________________________________________________
#read current directory from txt file
# do not forget to change the pat here -------------------------  INPUT !!!!!
setwd("C:\\Users\\ruidr\\Google Drive\\Trinity Sync\\PhD\\021_DICE")
getwd()
#________________________________________________________________
#load data from Worksheet of Matlab
#
workbookfile1 <- XLConnect::loadWorkbook("01data_input_R.xls")
alldata <- readWorksheet(workbookfile1, sheet = "Sheet1");
alldata
# assign the variables
n = ncol(alldata)
len = nrow(alldata)
#
XDoE <- alldata[ , c(1:(n-2))]
YDoE <- alldata[ , c(n-1)]
sigmaDoE<- alldata[ , c(n)]
#alocate variables
x <- (XDoE)
x <- mapply(x, FUN=as.numeric) #transform in numeric
y <- c(YDoE)
y <- mapply(y, FUN=as.numeric) #transform in numeric
#------------------------------------------------
#___________________ BUILD THE KRIGING _______________________________________________
# Heteregeneous noise variances
sigmaDoE <- mapply(sigmaDoE, FUN=as.numeric)
noise.var <- (sigmaDoE^2)
#------------------------------------------------
workbookfile2 <- XLConnect::loadWorkbook("02control_data_R.xls")
controldata <- readWorksheet(workbookfile2, sheet = "Sheet1");
workbookfile3 <- XLConnect::loadWorkbook("03control_convergence.xls")
bounds <- readWorksheet(workbookfile3, sheet = "Sheet1");
#DEGREE
degreeindex <- controldata[ (1), ]
degreeindex <- mapply(degreeindex, FUN=as.numeric)
#
#degreeindex <- 5
############### select degree
# trend (beware to the usual I operator) ------ First for the case of 1D
if ((n-2)==1){if (degreeindex == 0){formula <- ~1}
if (degreeindex == 1){formula <- ~x}
if (degreeindex == 2){formula <-~polym(x, degree=2, raw=TRUE)}
if (degreeindex == 3){formula <- ~polym(x, degree=3, raw=TRUE)}
if (degreeindex == 4){formula <- ~polym(x, degree=4, raw=TRUE) }
if (degreeindex == 5){formula <- ~polym(x, degree=5, raw=TRUE)}
if (degreeindex == 6){formula <- ~polym(x, degree=6, raw=TRUE)}
}
#~x+I(x^2)+I(x^3)+I(x^4)+I(x^5)+I(x^6)}
############### select degree
# First for the case of 2D
if ((n-2)==2){if (degreeindex == 0){formula <- ~1}
if (degreeindex == 1){formula <- ~polym(x1,x2, degree=1, raw=TRUE)}
if (degreeindex == 2){formula <- ~polym(x1,x2, degree=2, raw=TRUE)}
if (degreeindex == 3){formula <- ~polym(x1,x2, degree=3, raw=TRUE)}
if (degreeindex == 4){formula <- ~polym(x1,x2, degree=4, raw=TRUE) }
if (degreeindex == 5){formula <- ~polym(x1,x2, degree=5, raw=TRUE)}
if (degreeindex == 6){formula <- ~polym(x1,x2, degree=6, raw=TRUE)}
}
#x1+x2+I(x1*x2)+I(x1^2)+I(x2^2)}
#formula1 <- formula1~polym(x1, x2, degree=2, raw=TRUE) # is equivalent to
#formula <- ~x1 + x2 + I(x1^2) + I(x2^2) + x1:x2
#CORRELATION
elmcovtype <- controldata[ (2), ]
#BOUNDS
upperlim <- (bounds[ (2), 1:(n-2)])
upperlim <- mapply(upperlim, FUN=as.numeric)
lowerlim <- (bounds[ (1), 1:(n-2)])
lowerlim <- mapply(lowerlim, FUN=as.numeric)
# Create the Kriging model
# check manual to see the different options to fill in !!!
#is.numeric(data.frame(y=y))
#is.numeric(data.frame(x=x))
#rename the cpolumns in data frame --- very important
if ((n-2) == 1){x <- data.frame (x=x)}
if ((n-2) == 2){
colnames(x) <- c("x1","x2")}
#VV = matrix(NA, nrow=len, ncol=1)
#for (i in 1:len) {VV[i,1] = var[i]}
#noise.var = data.frame(s2=VV)
#}
ini_estimator <- sample(50:300, 1) #esseential for convergence for now!!!
#myKriging <- km(formula=formula, design = x, response = data.frame(y=y),
#               covtype=elmcovtype, noise.var=noise.var, estim.method="MLE",
#              optim.method = "BFGS", control=list(maxit=100, popsize=ini_estimator),
#             lower = lowerlim, upper = upperlim, multistart = 1)
myKriging <- km(formula=formula, design = x, response = data.frame(y=y),
covtype=elmcovtype, estim.method="MLE",
optim.method = "BFGS", control=list(maxit=100, popsize=ini_estimator),
lower = lowerlim, upper = upperlim, multistart = 1)
# Be carefull with the optimization boundaries
#"BFGS", "gen" , "L-BFGS-B"
#control=list(maxit=500,trace=T,popsize=50)
#get the theta, sigma and beta parameters from the Kriging.
theta.est <- covparam2vect(myKriging@covariance)
sigma2.est <- myKriging@covariance@sd2
beta.est <- myKriging@trend.coef
# save these in a file - a file for each
write.xlsx(theta.est, "11theta_res.xls", sheetName="Sheet1")
write.xlsx(sigma2.est, "12sigma2_res.xls", sheetName="Sheet1")
write.xlsx(beta.est, "13beta_res.xls", sheetName="Sheet1")
#------------------------------------------------
#covtype="matern5_2" or "gauss"
#_____________________________________________________________
# Other half of the inputs from the matlab
# prediction
library(XLConnect)
workbookfile4 <- XLConnect::loadWorkbook("04Points_DoE_predictions.xls")
predictionpoints <- readWorksheet(workbookfile4, sheet = "Sheet1");
m = ncol(predictionpoints)
t <- predictionpoints[ , c(1:m)]
if (m == 1){}
if (m == 2){ t = expand.grid( t[ , c(1)], t[ , c(2)] ) }
#t <- seq(min(x), max(x), by=0.1)
t <- mapply(t, FUN=as.numeric)
#rename the cpolumns in data frame --- very important for 2D
if (m == 1){}
if (m == 2){colnames(t) <- c("x1","x2")}
#--------------------------------------------
#---- PREDICCTIONS --------------------------
#
#introduce an excel file that will control the prediction points
# the prediction matrix needs to be a data frame!!!! Very important
if (m == 1){t <- data.frame (x=t)}
if (m == 2){t <- data.frame(t)}
p <- predict(myKriging, newdata=t, type="UK", se.compute=TRUE)
mu <- p$mean
sigma <- p$sd
# Intervals
#lower <- p$lower95; upper <- p$upper95
# ________________create the final matrix with a LOOP ____________
iter <- nrow(t)
Mfinal <- matrix(NA, nrow=iter, ncol=m+4)
#use for 1D and 2D
for (i in 1:iter){
if (m==1){ Mfinal[i,] <- c(t[i,m], mu[i], sigma[i], lower[i], upper[i])}
if (m==2){ Mfinal[i,] <- c(t[i,m-1],t[i,m], mu[i], sigma[i])} #, lower[i], upper[i])}
}
#-------------PLOT IF NEEDED-----------------------------
a=0
if (a == 1){
xx <- seq(1, length(t), 1)
plot(XDoE, YDoE, type = "p", ylim=c(min(lower),max(upper)))
lines (Mfinal[ , 1], Mfinal[ , 2], type = "l")
lines(Mfinal[ , 1], Mfinal[ , 4], type = "l", col='red')
lines(Mfinal[ , 1], Mfinal[ , 5], type = "l", col='red')
}
#____________________ SAVE in EXCEL to use in the MATLAB _________________
#Mfinal
#Mfinal <- m[ ]
#write.xlsx(Mfinal, "14KrigingR_dataresults.xls", sheetName="Sheet1")
#print to CSV
write.table(Mfinal, file = "all_predictor_data", col.names = FALSE, row.names = FALSE, sep=",")
#close sink
## reset message sink and close the file connection
#sink(type="message")
#close(zz)
#clean space of variables to avoid problems
rm(list = ls())
closeAllConnections()
.libPaths("C:/Users/ruidr/Documents/R/win-library/3.4")
.libPaths()
#install.packages("XLConnect", "C:/Users/Rui Teixeira/Documents/R/win-library/3.4")
rm(list=ls())
#____load all packages_____
#
library(DiceKriging)
library(XLConnect)
library(xlsx)
library(WriteXLS)
#
# do not forget to change the pat here -------------------------  INPUT !!!!!
setwd("C:\\Users\\ruidr\\Google Drive\\Trinity Sync\\PhD\\021_DICE\\Custom_fitR")
getwd()
#load data from Worksheet of Matlab
#
workbookfile1 <- XLConnect::loadWorkbook("01data_input_R.xls")
alldata <- readWorksheet(workbookfile1, sheet = "Sheet1");
alldata
# assign the variables
n = ncol(alldata)
len = nrow(alldata)
#
XDoE <- alldata[ , c(1:(n-2))]
YDoE <- alldata[ , c(n-1)]
#alocate variables
x <- (XDoE)
x <- mapply(x, FUN=as.numeric) #transform in numeric
y <- c(YDoE)
y <- mapply(y, FUN=as.numeric) #transform in numeric
#----------------------------------------------------
#------ Kriging Parameters ----------------
#----------------------------------------------------
workbookfile2 <- XLConnect::loadWorkbook("11theta_res.xls")
thetadata <- readWorksheet(workbookfile2, sheet = "Sheet1")
aux <- thetadata
thetadata <- mapply(aux, FUN=as.numeric)
workbookfile2 <- XLConnect::loadWorkbook("12sigma2_res.xls")
sigma2data <- readWorksheet(workbookfile2, sheet = "Sheet1")
aux <- sigma2data
sigma2data <- mapply(aux, FUN=as.numeric)
workbookfile3 <- XLConnect::loadWorkbook("13beta_res.xls")
betadata <- readWorksheet(workbookfile3, sheet = "Sheet1")
aux <- betadata
betadata <- mapply(aux, FUN=as.numeric)
# -----------------------------------
#Create the custom KRIGING model trend
#------------------------------------
workbookfile2 <- XLConnect::loadWorkbook("02control_data_R.xls")
controldata <- readWorksheet(workbookfile2, sheet = "Sheet1");
#DEGREE
degreeindex <- controldata[ (1), ]
degreeindex <- mapply(degreeindex, FUN=as.numeric)
############### select degree
# First for the case of 2D
if ((n-2)==2){if (degreeindex == 0){formula <- ~1}
if (degreeindex == 1){formula <- ~polym(x1,x2, degree=1, raw=TRUE)}
if (degreeindex == 2){formula <- ~polym(x1,x2, degree=2, raw=TRUE)}
if (degreeindex == 3){formula <- ~polym(x1,x2, degree=3, raw=TRUE)}
if (degreeindex == 4){formula <- ~polym(x1,x2, degree=4, raw=TRUE) }
if (degreeindex == 5){formula <- ~polym(x1,x2, degree=5, raw=TRUE)}
if (degreeindex == 6){formula <- ~polym(x1,x2, degree=6, raw=TRUE)}
}
#CORRELATION
elmcovtype <- controldata[ (2), ]
# -----------------------------------
#Create the custom KRIGING model trend
#------------------------------------
#rename the cpolumns in data frame --- very important
if ((n-2) == 1){x <- data.frame (x=x)}
if ((n-2) == 2){
colnames(x) <- c("x1","x2")}
#treat data
theta = thetadata [ , 2]
#theta = data.frame(theta = theta)
sigma = sigma2data [ 2]
#sigma = data.frame(sigma2 = sigma)
beta = betadata [ ,2]
#beta = data.frame(beta = beta)
#Kriging
myKriging_reinter <- km(formula=formula, design = x, response = data.frame(y=y),
covtype=elmcovtype, coef.trend = beta, coef.cov = theta , coef.var = sigma)
#coef.trend = beta,
# ------------------------------------
#USE the results to do the predictions
#-------------------------------------
# Other half of the inputs from the matlab
# prediction
library(XLConnect)
workbookfile4 <- XLConnect::loadWorkbook("04Points_search_prediction.xls")
predictionpoints <- readWorksheet(workbookfile4, sheet = "Sheet1");
m = ncol(predictionpoints)
t <- predictionpoints[ , c(1:m)]
if (m == 1){}
if (m == 2){ t = expand.grid( t[ , c(1)], t[ , c(2)] ) }
#t <- seq(min(x), max(x), by=0.1)
t <- mapply(t, FUN=as.numeric)
#rename the cpolumns in data frame --- very important for 2D
if (m == 1){}
if (m == 2){colnames(t) <- c("x1","x2")}
#--------------------------------------------
#---- PREDICCTIONS --------------------------
#
#introduce an excel file that will control the prediction points
# the prediction matrix needs to be a data frame!!!! Very important
if (m == 1){t <- data.frame (x=t)}
if (m == 2){t <- data.frame(t)}
p <- predict(myKriging_reinter, newdata=t, type="UK", se.compute=TRUE)
mu <- p$mean
sigma <- p$sd
# Intervals
lower <- p$lower95; upper <- p$upper95
# ________________create the final matrix with a LOOP ____________
iter <- nrow(t)
Mfinal <- matrix(NA, nrow=iter, ncol=m+4)
#use for 1D and 2D
for (i in 1:iter){
if (m==1){ Mfinal[i,] <- c(t[i,m], mu[i], sigma[i], lower[i], upper[i])}
if (m==2){ Mfinal[i,] <- c(t[i,m-1],t[i,m], mu[i], sigma[i], lower[i], upper[i])}
}
#write.xlsx(Mfinal, "99KrigingR_searchfield.xlsx", sheetName="Sheet1")
#Mfinal = data.frame(x = Mfinal)
#colnames(Mfinal) <- c("x1","x2", 'y', 'sigma', 'up', 'low')
#WriteXLS("Mfinal",ExcelFileName="99KrigingR_searchfield.xlsx")
#print to CSV
write.table(Mfinal, file = "allsearch_data", col.names = FALSE, row.names = FALSE, sep=",")
#wirte table does not write tables and columns !!
theta
theta
