library(gplots)


# get averages and standard deviations for phage timecourse data

dataTablePhage <- read.table( paste('experimentalData','phageTimecourse20111110-11_M9succinate_M9glucose.txt', sep= .Platform$file.sep), sep='\t', header = TRUE)
dataTablePhageAcetate <- read.table( paste('experimentalData','phageTimecourse20120223_Acetate.txt', sep= .Platform$file.sep), sep='\t', header = TRUE)
dataTablePhageRich <- read.table( paste('experimentalData','phageTimecourse_20110826_TB.txt', sep= .Platform$file.sep), sep='\t', header = TRUE)
dataTablePhageRich <- subset( dataTablePhageRich, time < 0.28  )
dataTablePhage$G1[2]

sdPhageGlucose <- 0
avePhageGlucose <- 0
for (i in 1:length(dataTablePhage$G1)) {
	sdPhageGlucose[i] <- sd( c(dataTablePhage$G1[i], dataTablePhage$G2[i], dataTablePhage$G3[i]) )
	avePhageGlucose[i] <- mean( c(dataTablePhage$G1[i], dataTablePhage$G2[i], dataTablePhage$G3[i]) )
}

sdPhageSuccinate <- 0
avePhageSuccinate <- 0
for (i in 1:length(dataTablePhage$S1)) {
	sdPhageSuccinate[i] <- sd( c(dataTablePhage$S1[i], dataTablePhage$S2[i], dataTablePhage$S3[i]) )
	avePhageSuccinate[i] <- mean( c(dataTablePhage$S1[i], dataTablePhage$S2[i], dataTablePhage$S3[i]) )
}

sdPhageTryptone <- 0
avePhageTryptone <- 0
for (i in 1:length(dataTablePhageRich$phage1)) {
	sdPhageTryptone[i] <- sd( c(dataTablePhageRich$phage1[i], dataTablePhageRich$phage2[i], dataTablePhageRich$phage3[i]) )
	avePhageTryptone[i] <- mean( c(dataTablePhageRich$phage1[i], dataTablePhageRich$phage2[i], dataTablePhageRich$phage3[i]) )
}

sdPhageAcetate <- 0
avePhageAcetate <- 0
for (i in 1:length(dataTablePhageAcetate$t)) {
	sdPhageAcetate[i] <- sd( c(dataTablePhageAcetate$A1[i], dataTablePhageAcetate$A2[i], dataTablePhageAcetate$A3[i]) )
	avePhageAcetate[i] <- mean( c(dataTablePhageAcetate$A1[i], dataTablePhageAcetate$A2[i], dataTablePhageAcetate$A3[i]) )
}


tComp <- 17

muData <- c(0.27,0.45,0.66,1.5)
phageData <- c(0,0,0,0)
phageSD <- c(0,0,0,0)

# Acetate for bar plot
idxMedia <- 1
idxAcetateData <- which.min( abs(dataTablePhageAcetate$t-tComp) )
phageData[idxMedia] <- avePhageAcetate[idxAcetateData]
phageSD[idxMedia] <- sdPhageAcetate[idxAcetateData]

# Succinate for bar plot
idxMedia <- 2
idxSuccinateData <- which.min( abs(dataTablePhage$time-tComp) )
phageData[idxMedia] <- avePhageSuccinate[idxSuccinateData]
phageSD[idxMedia] <- sdPhageSuccinate[idxSuccinateData]


# GLucose for bar plot
idxMedia <- 3
idxGlucoseData <- which.min( abs(dataTablePhage$time-tComp) )
phageData[idxMedia] <- avePhageGlucose[idxGlucoseData]
phageSD[idxMedia] <- sdPhageGlucose[idxGlucoseData]

# Tryptone for bar plot
idxMedia <- 4
idxTryptoneData <- which.min( abs(dataTablePhageRich$time-tComp) )
phageData[idxMedia] <- avePhageTryptone[idxTryptoneData]
phageSD[idxMedia] <- sdPhageTryptone[idxTryptoneData]

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
folderString <- paste( 'simulationsText', 'fluxBoundsSet_metabolicLimit' , sep= .Platform$file.sep)
files <- list.files( path = folderString, pattern = ".txt" )
nFile <- length( files )

metLimitMus <- rep(0,nFile)
metLimitPhage <- rep(0,nFile)

for( iFile in 1:nFile ){
	
	dataTableTemp <- read.table(paste(folderString,files[iFile],sep= .Platform$file.sep), sep='\t', header = TRUE)
	dataTableTemp <- subset( dataTableTemp, hostFlag == 1 )
	idxDataTemp <- which.min( abs(dataTableTemp$t*60-tComp) )
	
	metLimitMus[iFile] <- dataTableTemp$bBIOMASS[2]
	metLimitPhage[iFile] <- dataTableTemp$CONCvT7[idxDataTemp]
	
}

# corners of region

# # if convex
# metLimitMus <- c( metLimitMus, min(metLimitMus), max(metLimitMus) )
# metLimitPhage <- c( metLimitPhage, 0, 0 )
# metLimitRegion <- chull( x = metLimitMus, y = metLimitPhage )
# xMetLimit <- metLimitMus[metLimitRegion]
# yMetLimit <- metLimitPhage[metLimitRegion]

# smoothing
fSmooth <- 0.1
smoothed <- lowess(x=metLimitMus,y=metLimitPhage,f=fSmooth)
xMetLimit <- c(min(metLimitMus), smoothed$x ,max(metLimitMus))
yMetLimit <- c(0, smoothed$y ,0)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
folderString <- paste( 'simulationsText', 'fluxBoundsSet_mechineryLimit' , sep= .Platform$file.sep)
files <- list.files( path = folderString, pattern = ".txt" )
nFile <- length( files )

intLimitMus <- rep(0,nFile)
intLimitPhage <- rep(0,nFile)
mechLimitPhage <- rep(0,nFile)

for( iFile in 1:nFile ){
	
	dataTableTemp <- read.table(paste(folderString,files[iFile],sep= .Platform$file.sep), sep='\t', header = TRUE)
	dataTableTemp <- subset( dataTableTemp, hostFlag == 1 )
	idxDataTemp <- which.min( abs(dataTableTemp$t*60-tComp) )
	
	intLimitMus[iFile] <- dataTableTemp$bBIOMASS[2]
	intLimitPhage[iFile] <- dataTableTemp$CONCvT7[idxDataTemp]
	mechLimitPhage[iFile] <- dataTableTemp$CONCvT7Ref[idxDataTemp]
	
}

# corners of region
mechLimitSortIdx <- order( intLimitMus )
xMechLimitRegion <- c( min(intLimitMus), intLimitMus[mechLimitSortIdx] ,max(intLimitMus) )
yMechLimitRegion <- c( 0, mechLimitPhage[mechLimitSortIdx] ,0)
xInfeasRegion <- c( min(intLimitMus), intLimitMus[mechLimitSortIdx] )
yInfeasRegion <- c( max(mechLimitPhage), mechLimitPhage[mechLimitSortIdx])

# sort integrated data
intSortIdx <- order( intLimitMus )
intLimitMus <- intLimitMus[intSortIdx]
intLimitPhage <- intLimitPhage[intSortIdx]

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -	

		
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

figureName <- paste('figure7_floatLegend_8o7cm_width_',tComp,'min.eps', sep = '' )

# postscript(figureName, width = 7.0, height = 6.0, horizontal = FALSE, onefile = FALSE, paper = "special")
postscript(paste('figureOutput',figureName, sep= .Platform$file.sep), width = 3.27, height = 3.27, horizontal = FALSE, onefile = FALSE, paper = "special", pointsize = 8)

par( mar = (c(5, 5, 0, 0) + 0.1))

xAxis <- c(min(intLimitMus), 1.65)
yAxis <- c(0, 240)

xScale <- seq(from = 0.25, to = 1.5, by = 0.25)
yScale <- seq(from = 0, to = 200, by = 50)

yString <- paste('Phage T7 Produced (phage/infected host) \n at ', tComp, ' Minutes Post Infection', sep='')

plot(x=NULL, y=NULL, xlab = 'Growth Rate (/hour)', xlim = xAxis, ylab = yString, ylim = yAxis, axes = FALSE)

# print(xMechLimitRegion)
# print(xMetLimit)
polygon( x=xInfeasRegion, y =yInfeasRegion, col = 'grey95' , border = NA)
polygon( x=xMechLimitRegion, y =yMechLimitRegion, col = 'grey70' , border = NA)
# polygon( x=xMetLimit, y = yMetLimit , col = 'olivedrab1', border = NA )
# polygon( x=xMechLimitRegion, y =yMechLimitRegion, col = 'black', density = 10, lwd = 1, angle = 45 )
polygon( x=xMetLimit, y = yMetLimit , col = 'black' , density = 10, lwd = 1, angle = 135 )
 # points(x=xMechLimitRegion, y =yMechLimitRegion,pch = 4, col = 'grey30')
 # points( x=xMetLimit, y = yMetLimit,pch = 4, col = 'olivedrab3')
lines(x=intLimitMus,y=intLimitPhage, lwd = 4 , col =  'darkolivegreen3' ) # 'orangered1'

plotCI( x = muData, y = phageData, uiw = phageSD, col = 'cornflowerblue' , barcol = 'cornflowerblue' , add = TRUE, cex = 1.5, lwd = 1.5)
axis(1, at = xScale )
axis(2, at = yScale)

# legend(x=1.75, y=250, legend=c('Mechanistically Infeasible','Mechanistically Feasible','Metabolicly Feasible ','Integrated Simulation','Expirimental Data'),  
# pch = c(NA, NA, NA, NA, 1), lty = c(NA,NA,NA,1, NA), lwd = c(1,1,1,4,2),col=c(NA, NA ,"black",'grey20','black'),
# fill=c('grey90','grey50',NA,NA,NA), bty="n", border = c(NA,NA, 'black', NA, NA),  pt.cex = c(4,4,4,NA,1), seg.len = 1, 
# density = c(NA, NA, 10, NA, NA), angle = c(NA, NA, 135, NA, NA) )

wLeg <- 0.12
hLeg <- 15
yLeg <- max(yInfeasRegion) - 10 - hLeg#235
xLeg <- min(xInfeasRegion)+0.1 #0.25
sLeg <- 1.75

rect(xleft = (xLeg-0.05), ybottom = (yLeg-4*sLeg*hLeg-0.35*hLeg ) , xright= (xLeg+(0.8-0.2)), ytop = (yLeg+1.35*hLeg), col = 'white', border = NA)

yLegBot <- yLeg-(0:2)*(sLeg*hLeg) -2
yLegTop <- yLegBot + hLeg + 2

rect(xleft = xLeg, ybottom = yLegBot , xright= (xLeg+wLeg), ytop = yLegTop, density = c(NA, NA, 10),
col=c('grey95','grey70','black'),border = c(NA,NA, 'black'), angle = c(NA, NA, 135))

yLegBot <- yLeg-(0:4)*(sLeg*hLeg)
yLegTop <- yLegBot + hLeg

text(c('Machinery \nInfeasible','Machinery \nFeasible','Metabolically \nFeasible','Integrated \nSimulation','Experimental \nData'),
x= (xLeg+wLeg*1.2), y = (yLegBot+ yLegTop)/2,adj = c(0,0.5))
lines(x=c(xLeg,(xLeg+wLeg)),y=rep(((yLegBot[4]+ yLegTop[4])/2),2), lwd = 4 , col = 'darkolivegreen3' )		
points(x=c((xLeg+wLeg/2)),y=(yLegBot[5]+ yLegTop[5])/2, col =  'cornflowerblue', cex = 1.5, lwd =1.5 )		
dev.off()