### R code from vignette source 'xadditivity.Rnw'
### Encoding: ISO8859-1

###################################################
### code chunk number 1: setup
###################################################

require(CHNOSZ)

options(width=65,continue=" ")



###################################################
### code chunk number 2: SolveFour
###################################################
A <- matrix(c(2,2,1,1, 0,1,0,0, 0,0,1,0, 0,0,0,1), 4)
b <- c(88.3,110.6,62.2,40.56)
x <- solve(A,b)
names(x) <- c("CH3","CH2","CH2OH","COOH")
x


###################################################
### code chunk number 3: CpAmendHelgeson
###################################################
Cp.AH97 <- c(CH3=47.00,CH2=20.70,CH2OH=17.28,COOH=-5.94,CHCH3=39.92)
Cp.AH97


###################################################
### code chunk number 4: BigMatrixAdd
###################################################
# xtable for making nice LaTeX tables
require(xtable)
# read the "big group" specification
f <- system.file("extdata/thermo/groups_big.csv",package="CHNOSZ")
A.big <- read.csv(f,check.names=FALSE,row.names=1)
# only take selected rows
imod <- c("ethane","propane","butane","ethanol","1-propanol","1-butanol","1-pentanol",
  "acetone","butanone","3-pentanone","2-heptanone",
  "acetic acid","propanoic acid","butanoic acid")
A.big <- A.big[imod,]
# get the heat capacities
is1 <- info(rownames(A.big),quiet=TRUE)
mydf <- data.frame(Cp=info(is1,quiet=TRUE)$Cp,A.big,check.names=FALSE)
print(xtable(mydf,align=c("l","r","r","r","r","r","r")))


###################################################
### code chunk number 5: QRexample
###################################################
set.seed(24)
print(A <- matrix(runif(12), 4))
b <- 1:4
x <- qr.solve(A,b)
as.numeric(A %*% x)


###################################################
### code chunk number 6: ReadBigAdd
###################################################
file.big <- system.file("extdata/thermo/groups_big.csv",package="CHNOSZ")
A.big <- read.csv(file.big,check.names=FALSE,row.names=1)
iuse <- c("ethane","propane","butane","ethanol","1-propanol","1-butanol","1-pentanol",
  "acetone","butanone","3-pentanone","2-heptanone",
  "acetic acid","propanoic acid","butanoic acid")
A.big <- A.big[iuse,]
A.big[is.na(A.big)] <- 0


###################################################
### code chunk number 7: CpSpeciesAdd
###################################################
ispecies <- info(rownames(A.big),quiet=TRUE)
cp.species <- info(ispecies)$Cp


###################################################
### code chunk number 8: SolveBigAdd
###################################################
cp.big <- qr.solve(A.big,cp.species)
cp.big


###################################################
### code chunk number 9: RMSDBigAdd
###################################################
pred.big <- as.matrix(A.big) %*% cp.big
rmsd <- rmsd(pred.big,cp.species)
rmsd


###################################################
### code chunk number 10: ResidualsBigAdd (eval = FALSE)
###################################################
## residuals <- pred.big - cp.species
## names(residuals) <- rownames(A.big)
## residuals <- residuals
## residualsplot(residuals,"Cp","big groups (added)")


###################################################
### code chunk number 11: ResidualsPlotBigAdd
###################################################
residuals <- pred.big - cp.species
names(residuals) <- rownames(A.big)
residuals <- residuals
residualsplot(residuals,"Cp","big groups")


###################################################
### code chunk number 12: BigMatrixSub
###################################################
# xtable for making nice LaTeX tables
require(xtable)
# read the "big group" specification
f <- system.file("extdata/thermo/groups_big.csv",package="CHNOSZ")
A.big <- read.csv(f,check.names=FALSE,row.names=1)
# get the heat capacities
is1 <- info(rownames(A.big),quiet=TRUE)
mydf <- data.frame(Cp=info(is1,quiet=TRUE)$Cp,A.big,check.names=FALSE)
rownames(mydf) <- paste(paste("(",1:nrow(mydf),")",sep=""),rownames(mydf))
print(xtable(mydf,align=c("l","r","r","r","r","r","r")))


###################################################
### code chunk number 13: ReadBigSub
###################################################
file.big <- system.file("extdata/thermo/groups_big.csv",package="CHNOSZ")
A.big <- read.csv(file.big,check.names=FALSE,row.names=1)
A.big[is.na(A.big)] <- 0


###################################################
### code chunk number 14: CpBigSub
###################################################
ispecies <- info(rownames(A.big),quiet=TRUE)
cp.species <- info(ispecies,quiet=TRUE)$Cp


###################################################
### code chunk number 15: SolveBigSub
###################################################
i.big <- 1:24
cp.big <- qr.solve(A.big[i.big,],cp.species[i.big])
cp.big


###################################################
### code chunk number 16: RMSDBigSub
###################################################
pred.big <- as.matrix(A.big) %*% cp.big
rmsd <- rmsd(pred.big[i.big],cp.species[i.big])
rmsd


###################################################
### code chunk number 17: ResidualsBigSub (eval = FALSE)
###################################################
## residuals <- pred.big - cp.species
## names(residuals) <- rownames(A.big)
## residuals <- residuals[i.big]
## residualsplot(residuals,"Cp","big groups")


###################################################
### code chunk number 18: ResidualsPlotBigSub
###################################################
residuals <- pred.big - cp.species
names(residuals) <- rownames(A.big)
residuals <- residuals[i.big]
residualsplot(residuals,"Cp","big groups")


###################################################
### code chunk number 19: SmallMatrix
###################################################
# read the "small group" specification
file.small <- system.file("extdata/thermo/groups_small.csv",package="CHNOSZ")
A.small <- read.csv(file.small,check.names=FALSE,row.names=1)
# add stars to indicate species from Aug.11
ispecies <- info(head(rownames(A.small),-1),quiet=TRUE)
inew <- grep("Aug.11",info(ispecies,quiet=TRUE)$date)
row.names(A.small)[inew] <- paste(row.names(A.small)[inew],"*")
xt <- xtable(A.small,align=c("p{4.5cm}",rep("p{0.6cm}",12)))
# now perform some work on the LaTeX code to rotate the column names
xp <- capture.output(print(xt))
# column names are in line 7
stuff <- xp[[7]]
# remove the \\\\ at end
stuff <- strsplit(stuff,"\\\\\\\\ ")[[1]]
# split the text at the latex separator
stuff <- strsplit(stuff,"& ")[[1]]
# remove the separating spaces
stuff <- gsub("\ $","",stuff)
# remember which ones are blank
iblank <- stuff==""
# add the latex rotation instructions
# have to escape the backslash
stuff <- paste("\\rotatebox{60}{",stuff,"}",sep="")
# restore the blanks
stuff[iblank] <- ""
# restore the separating latex separator character
stuff <- list(paste(stuff,"& "))
# unsplit it
stuff <- sapply(stuff,paste,collapse="")
# restore the \\\\ at end
stuff <- gsub("& $","\\\\\\\\ ",stuff)
# put in our modified code for column names
xp[[7]] <- stuff
cat(xp,sep="\n")


###################################################
### code chunk number 20: ReadSmall
###################################################
file.small <- system.file("extdata/thermo/groups_small.csv",package="CHNOSZ")
A.small <- read.csv(file.small,check.names=FALSE,row.names=1)
A.small[is.na(A.small)] <- 0


###################################################
### code chunk number 21: CpSpeciesSmall
###################################################
ispecies <- info(head(rownames(A.small),-1),quiet=TRUE)
cp.species <- info(ispecies,quiet=TRUE)$Cp
cp.species <- c(cp.species,NA)


###################################################
### code chunk number 22: SolveSmall
###################################################
i.small <- 1:24
cp.small <- qr.solve(A.small[i.small,],cp.species[i.small])
cp.small


###################################################
### code chunk number 23: RMSDSmall
###################################################
pred.small <- as.matrix(A.small) %*% cp.small
rmsd <- rmsd(pred.small[i.small],cp.species[i.small])
rmsd


###################################################
### code chunk number 24: ResidualsSmall (eval = FALSE)
###################################################
## residuals <- pred.small - cp.species
## names(residuals) <- rownames(A.small)
## residuals <- residuals[i.small]
## residualsplot(residuals,"Cp","small groups")


###################################################
### code chunk number 25: ResidualsPlotSmall
###################################################
residuals <- pred.small - cp.species
names(residuals) <- rownames(A.small)
residuals <- residuals[i.small]
residualsplot(residuals,"Cp","small groups")


###################################################
### code chunk number 26: LacCit
###################################################
ilac <- match("lactic acid",rownames(A.small))
icit <- match("citric acid",rownames(A.small))
iiso <- match("isocitric acid",rownames(A.small))
pb <- round(pred.big,2)
ps <- round(pred.small,2)


