### 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)