---
title: "R Notebook"
output: html_notebook
---


Nejistoty a chyby měření

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(errors)
xe <- set_errors(3.602, 0.008)
options(errors.notation = "plus-minus") 
print(xe, digits = 2)
options(errors.notation = "parenthesis")
print(xe, digits = 2)

# elementarni naboj
e <- set_errors(1.6021766208e-19, 0.0000000098e-19) 
options(errors.notation = "plus-minus") 
print(e, digits = 2)
options(errors.notation = "parenthesis")
print(e, digits = 3)

```


Ze zásilky kaprolaktamu bylo odebráno 10 vzorku a byl u nich stanoven bod tání. Vypočítejte průměrnou hodnotu bodu tání v zásilce a její směrodatnou odchylku.

```{r echo=FALSE, message=FALSE, warning=FALSE}

xc = c(68.5, 68.7, 68.3, 68.8, 68.5, 68.2, 68.6, 68.4, 68.2, 68.7)
mean(xc)
sd(xc) # sd
sd(xc)/sqrt(length(xc)) # str chyba prumeru

e <- set_errors(mean(xc), sd(xc)/sqrt(length(xc))) 
options(errors.notation = "plus-minus"); print(e, digits = 2)

```


Šíření chyb

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(propagate)

# z prumeru a nejistoty, bez uvedení poctu stupnu volnosti, resp. poctu opakovani

x <- c(5, 0.01) 
y <- c(1, 0.01) 

EXPR1 <- expression(x/y) 
DF1 <- cbind(x, y)
RES1 <- propagate(expr = EXPR1, data = DF1, type = "stat", do.sim = TRUE, verbose = TRUE, nsim = 100000) 
RES1
RES1$data
RES1$prop
RES1$sim
summary(RES1)
plot(RES1)

EXPR <- expression(a^b*x)
a = c(5, 0.1)
b = c(10, 0.1)
x = c(1, 0.1)
DAT <- cbind(a, b, x)
(res <- propagate(EXPR, DAT))
res$data
res$prop
res$sim
summary(res)
plot(res)


# z prumeru a nejistoty, s uvedením poctu stupnu volnosti
EXPR2 <- expression(x/y) 
x <- c(5, 0.01, 12) 
y <- c(1, 0.01, 5) 
DF2 <- cbind(x, y) 
RES2 <- propagate(expr = EXPR2, data = DF2, type = "stat", do.sim = TRUE, verbose = TRUE, nsim = 100000)
RES2
RES2$data
RES2$prop
RES2$sim
summary(RES2)
plot(RES2)


# Výpocet pomocí intervalu 

EXPR3 <- expression(C * sqrt((520 * H * P)/(M *(t + 460))))
H <- c(64, 65)
M <- c(16, 16.2)
P <- c(361, 365)
t <- c(165, 170)
C <- c(38.4, 38.5)
DAT3 <- makeDat(EXPR3)
interval(DAT3, EXPR3, seq = 2)

EXPR5 <- expression(x^2 - x + 1)
x <- c(-2, 1)
curve(x^2 - x + 1, -2, 1)
DAT5 <- makeDat(EXPR5)
interval(DAT5, EXPR5, seq = 2)


library(metRology)

data(GUM.H.1)
GUM.H.1

## a simple uncertainty analysis for the product of two quantities
GUM(c("x1","x2"),c(2.3,1.1),c(0.030,0.015),c(5,9999),"x1*x2")

## a simple uncertainty analysis for the product of two quantities
GUM.validate(c("x1","x2"), c(2.3,1.1), c(0.030,0.015), c(5,9999), c("A","B"),c("Normal","Rectangular"),"x1*x2")


expr <- expression(a+b*2+c*3+d/2)
x <- list(a=1, b=3, c=2, d=11)
u <- lapply(x, function(x) x/10) # nejistota 10 %
u.expr<-uncert(expr, x, u, method="NUM")
u.expr
# function method
f <- function(a,b,c,d) a+b*2+c*3+d/2
u.fun<-uncert(f, x, u, method="NUM")
u.fun
# formula method
u.form<-uncert(~a+b*2+c*3+d/2, x, u, method="NUM")
u.form


# s korelaci
u.cor<-diag(1,4)
u.cor[3,4]<-u.cor[4,3]<-0.5
u.cor
# num
u.formc<-uncert(~a+b*2+c*3+d/2, x, u, method="NUM", cor=u.cor)
u.formc
# Monte Carlo
u.formc.MC<-uncert(~a+b*2+c*3+d/2, x, u, method="MC", cor=u.cor, B=200)
u.formc.MC


expr <- expression(a/(b-c))
x <- list(a=1, b=3, c=2)
u <- lapply(x, function(x) x/20)
set.seed(403)
u.invexpr<-uncertMC(expr, x, u, distrib=rep("norm", 3), B=999, keep.x=TRUE )
u.invexpr

```


Vypočítejte hustotu a její chybu pro látku, u níž byla opakovaným měřením stanovena hmotnost 6.824 (0.008) g a objem 3.03 (0.01) ml.

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(propagate)

EXPR2 <- expression(m/V) 
m = c(6.824, 0.008) #[g] 
V = c(3.03, 0.01) #[ml]
DF2 <- cbind(m, V) 
RES2 <- propagate(expr = EXPR2, data = DF2, type = "stat", do.sim = TRUE, verbose = TRUE, nsim = 100000)
RES2
RES2$data
RES2$prop
RES2$sim
library(errors)
e <- set_errors(RES2$sim[1], RES2$sim[2]) 
options(errors.notation = "plus-minus"); print(e, digits = 1)
options(errors.notation = "parenthesis"); print(e, digits = 1)

```


Na píst o průměru 200 (0.05) mm působí pára tlakem 8.2 (0.1) atm. Jakou silou působí pára na píst?

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(propagate)
library(measurements)

d = c(200, 0.05) # [mm] na [m]
d = as.vector(conv_unit(c(200, 0.05), from="mm", to="m"))
p = c(8.2, 0.1)  # [atm] na [Pa]
p = as.vector(conv_unit(c(8.2, 0.1), from="atm", to="Pa"))
pi
EXPR7 <- expression(p*3.141593*(d/2)^2)
DF7 <- cbind(d, p) 
RES7 <- propagate(expr = EXPR7, data = DF7, type = "stat", do.sim = TRUE, verbose = TRUE, nsim = 100000)
RES7
RES7$data
RES7$prop
RES7$sim
library(errors)
e <- set_errors(RES7$sim[1], RES7$sim[2]) 
options(errors.notation = "plus-minus"); print(e, digits = 1)
options(errors.notation = "parenthesis"); print(e, digits = 1)

```


Vypočítejte koeficient viskozity roztoku glycerinu Stokesovou metodou

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(propagate)

r = c(0.0112, 0.0001) # polomer kulicky [cm]
l = c(31.23, 0.05) # draha kulicky za cas t [cm]
t = c(62.1, 0.2) # cas [s]
g = c(980.1,0) # tihove zrychleni [cm/s2] z tabulek
d0 = c(13.55,0) # hustota kulicky [g/cm3] z tabulek
d = c(1.28,0) # hustota roztoku [g/cm3] z tabulek
EXPR3 <- expression((2/9)*g*(((d0-d)*r^2)/l)*t)
DF3 <- cbind(r, l, t, g, d0, d) 
RES3 <- propagate(expr = EXPR3, data = DF3, type = "stat", do.sim = TRUE, verbose = TRUE, nsim = 100000)
RES3
RES3$data
RES3$prop
RES3$sim
EXPR4 <- expression((2/9)*980.1*(((13.55-1.28)*r^2)/l)*t)
DF4 <- cbind(r, l, t) 
RES4 <- propagate(expr = EXPR4, data = DF4, type = "stat", do.sim = TRUE, verbose = TRUE, nsim = 100000)
RES4$data
RES4$prop
RES4$sim
library(errors)
e <- set_errors(RES4$sim[1], RES4$sim[2]) 
options(errors.notation = "plus-minus"); print(e, digits = 1)
options(errors.notation = "parenthesis"); print(e, digits = 1)

```


Vypočítejte koeficient viskozity roztoku glycerinu pomocí kapilárního viskozimetru.

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(propagate)
library(measurements)

p = c(20.12, 0.01) # tlak na vytoku z kapilary [mm Hg] to [Pa]
p = as.vector(conv_unit(c(20.12, 0.01), from="mmHg", to="Pa"))
r = c(0.570, 0.003) # polomer kapilary [mm] to [m]
r = conv_unit(c(0.570, 0.003), from="mm", to="m")
l = c(10.526, 0.005) # delka kapilary [mm] to [m]
l = conv_unit(c(10.526, 0.005), from="mm", to="m")
V = c(5.025, 0.001)# objem kapaliny vytekle za cas t [cm3] to [m3]
V = conv_unit(c(5.025, 0.001), from="cm3", to="m3")
t = c(27.34, 0.02) # cas [s]
pi
EXPR5 <- expression((3.141593*p*(r^4)*t)/(8*V*l))
DF5 <- cbind(p, r, l, V, t) 
RES5 <- propagate(expr = EXPR5, data = DF5, type = "stat", do.sim = TRUE, verbose = TRUE, nsim = 100000)
RES5 # [kg / m.s]
RES5$data
RES5$prop
RES5$sim
# conv_multiunit(x = RES5$sim[1:2], from="kg / m", to="g / cm")
library(errors)
e <- set_errors(RES5$sim[1], RES5$sim[2]) 
options(errors.notation = "plus-minus"); print(e, digits = 1)
options(errors.notation = "parenthesis"); print(e, digits = 1)

```


Součin rozpustnosti stříbrné soli AgX má hodnotu Ks = 4.0 (0.4) x 10^-8. Jaká je chyba vypočtené rovnovážné koncentrace stříbrných iontů ve vode?

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(propagate)
EXPR <- expression(x^0.5)
x = c(4.0, 0.4)
DAT <- data.frame(x)
res <- propagate(EXPR,DAT)
# NEFUNGUJE pro jednu promennou.


library(metRology)

GUM("Ks",4.0e-8,0.4e-8,1,"sqrt(Ks)",sig.digits.U = 2) 

expr <- expression(a^0.5)
x <- list(a=4.0e-8)
u <- lapply(x, function(x) x/10) # nejistota 10 %
u.expr<-uncert(expr, x, u, method="NUM")
u.expr
# function method
f <- function(a) a^0.5
u.fun<-uncert(f, x, u, method="NUM")
u.fun
# formula method
u.form<-uncert(~a^0.5, x, u, method="NUM")
u.form

```


Nejistoty a korelace

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(errors)

x = c(0.9719378, 1.9840006, 2.9961830, 4.0123346, 5.0012799)
errors(x) = c(0.01, 0.01, 0.01, 0.01, 0.01)
y = c(0.9748992, 1.9627805, 2.9935831, 3.9921237, 4.9612555)
errors(y) = c(0.02, 0.02, 0.02, 0.02, 0.02)

# bez korelace
correl(x, y) =  c(0, 0, 0, 0, 0) 
z <- x / y; z

# s korelací
correl(x, y) =  c(0.8864282, 0.9761841, 0.9140209, 0.9496266, 0.9911837)
z_correlp <- x / y; z_correlp

correl(x, y) =  c(-0.8864282, -0.9761841, -0.9140209, -0.9496266, -0.9911837)
z_correln <- x / y; z_correln

```


Import dat

```{r echo=FALSE, message=FALSE, warning=FALSE}


# Import Dataset v RStudio


# clipboard

cyclamate<- read.delim("clipboard", header=T)


# txt, ascii

cyclamate1<- read.table("c:\\Users\\lubop\\Dropbox\\kursR\\cyclamate.txt", header=T)

cyclamate2<- read.delim("c:\\Users\\lubop\\Dropbox\\kursR\\cyclamate.txt", header=T)


## csv

data1 <- read.csv("c:\\Users\\lubop\\Dropbox\\kursR\\Chemical Composion of Ceramic.csv", header=TRUE, stringsAsFactors=FALSE)

library(data.table)
data2 <- fread("c:\\Users\\lubop\\Dropbox\\kursR\\Chemical Composion of Ceramic.csv")

data3 <- read.csv("https://archive.ics.uci.edu/ml/machine-learning-databases/00583/Chemical Composion of Ceramic.csv")


library(readr)


## Excel

library(openxlsx)

getSheetNames("c:\\Users\\lubop\\Dropbox\\kursR\\zarodky.xlsx")

zarodky1 <- read.xlsx("c:\\Users\\lubop\\Dropbox\\kursR\\zarodky.xlsx", sheet = "zarodky", startRow = 1, colNames = TRUE, rowNames = FALSE,detectDates = FALSE, skipEmptyRows = TRUE,rows = NULL, cols = NULL,check.names = FALSE)
head(zarodky1)

library(readxl)
zarodky2 <- read_excel("c:\\Users\\lubop\\Dropbox\\kursR\\zarodky.xlsx", sheet = "zarodky")
head(zarodky2)


path = "c:\\Users\\lubop\\Dropbox\\kursR\\"
file.names <- dir(path, pattern =".txt"); file.names

RD <- function(X){M <- read.delim(paste("c:\\Users\\lubop\\Dropbox\\kursR\\",X,sep=""))
return(M)}
listt = lapply(file.names,RD)
listt[[1]]
listt[[2]]
listt[[4]]
listt[[3]]

```


Export dat

```{r echo=FALSE, message=FALSE, warning=FALSE}

# txt, ascii

write.table(listt[[3]], file="c:\\Users\\lubop\\Dropbox\\kursR\\slinutypr2.txt", sep='\t')


# csv

write.csv(listt[[3]], "c:\\Users\\lubop\\Dropbox\\kursR\\slinutypr.csv", row.names=FALSE)

library(readr)
write_csv(listt[[4]], "c:\\Users\\lubop\\Dropbox\\kursR\\zarodky.csv")
library(data.table)
fwrite(listt[[2]], "c:\\Users\\lubop\\Dropbox\\kursR\\pneu.csv")


# Excel

library(openxlsx)
write.xlsx(listt[[1]], "c:\\Users\\lubop\\Dropbox\\kursR\\cyklamaty1.xlsx", asTable = FALSE, overwrite = TRUE)

library(writexl)
write_xlsx(listt[[1]], "c:\\Users\\lubop\\Dropbox\\kursR\\cyklamaty2.xlsx")

```


Čísla dle Chemical Abstracts Service (CAS)

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(erah)

findComp(name = "caffeine", id.database = mslib, CAS = NULL,chem.form = NULL)
findComp(name = "proline", id.database = mslib, CAS = NULL,chem.form = NULL)

findComp(name = NULL, id.database = mslib, CAS = "58-08-2",chem.form = NULL)
findComp(name = NULL, id.database = mslib, CAS = "50-78-2",chem.form = NULL)


library(httk)

chem.physical_and_invitro.data

CAS.checksum(CAS.string="50-78-2") # Aspirin
CAS.checksum(CAS.string="58-08-2") # Caffeine

```


BCPC Compendium of Pesticide Common Names https://pesticidecompendium.bcpc.org
     formerly known as 
Alan Woods Compendium of Pesticide Common Names
http://www.alanwood.net/pesticides

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(webchem)

# use names

bcpc_query("DEET", type = 'commonname', verbose = TRUE)

bcpc_query("Fluazinam", type = 'commonname', verbose = TRUE)$fluazinam$formula[1]

bcpc_query('Parathion', from ='name')

out = bcpc_query(c('Fluazinam','Diclofop'), from ='name')
out


# use CAS-numbers

bcpc_query("79622-59-6", from ='cas')

bcpc_query("51-03-6", from = 'cas', verbose = TRUE)
# from = c("name", "cas")


```


ChemicalIdentifierResolver(CIR)
http://cactus.nci.nih.gov/chemical/

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(webchem)

cir_query("piperonyl butoxide", representation = "smiles", resolver = NULL, first = FALSE, verbose = TRUE)
cir_query("51-03-6", representation = "smiles", resolver = NULL, first = FALSE, verbose = TRUE)

# SMILES
cir_query('Imidacloprid', representation = 'smiles')
cir_query('Aspirin', 'smiles')
cir_query('Triclosan', 'smiles')
cir_query("3380-34-5", 'smiles')

# CAS
cir_query('Imidacloprid', representation = 'cas')
cir_query('DEET', 'cas', first = TRUE)
cir_query('Triclosan', 'cas')
cir_query("3380-34-5", 'cas', first = TRUE)
cir_query("3380-34-5", 'cas', resolver = 'cas_number')

# InChIKey
cir_query('Imidacloprid', representation = 'stdinchikey')

# Molecular weight
cir_query('Imidacloprid', representation = 'mw')

# number of rings
cir_query('Imidacloprid', representation = 'ring_count')

# formula
cir_query("3380-34-5", 'formula')

# name
cir_query("3380-34-5", 'names')


cir_query("3380-34-5", 'pubchem_sid')
cir_query("3380-34-5", 'chemnavigator_sid')


name = 'Triclosan'
cir_query(name, 'mw')
cir_query(name, 'formula')
cir_query(name, 'monoisotopic_mass')
cir_query(name, 'heteroatom_count')
cir_query(name, 'hydrogen_atom_count')


comp <- c('Triclosan', 'Aspirin')
cir_query(comp, 'cas')
cir_query(comp, 'cas', first = TRUE)
cir_query(comp, 'smiles')
cir_query(comp, 'mw')


cir_img("CCO", dir = "c:\\Users\\lubop\\Dropbox\\kursR\\") # SMILES

query = c("Glyphosate", "Isoproturon", "BSYNRYMUTXBXSQ-UHFFFAOYSA-N")
cir_img(query, dir = "c:\\Users\\lubop\\Dropbox\\kursR\\", bgcolor = "transparent", antialising = 0)

query  = "Triclosan"
cir_img(query,dir = "c:\\Users\\lubop\\Dropbox\\kursR\\",format = "png",width = 600,height = 600,linewidth = 2,symbolfontsize = 30,bgcolor = "white",antialiasing = FALSE,atomcolor = "black",bondcolor = "black",csymbol = "all",hsymbol = "all",hcolor = "black",header = "",footer = "",frame = 1,verbose = getOption("verbose"))

```


ChemicalIdentifierResolver(CIR)
https://cactus.nci.nih.gov/chemical/structure

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(MSbox)

describe('Triclosan', "formula", info = TRUE) 

describe('malic acid', "formula", info = TRUE) 
describe('malic acid', "mw", info = FALSE) 
describe('malic acid', "hydrogen_atom_count", info = FALSE) 
describe('malic acid', "heteroatom_count", info = FALSE) 

describe(c('malic acid', 'citric acid', 'tartaric acid'), "smiles")

describe('glyphosate', "formula", info = FALSE) 
describe('glyphosate', "mw", info = FALSE)
describe('glyphosate', "smiles", info = FALSE) 
describe('glyphosate', "stdinchikey", info = FALSE) 

```


Chemical Translation Service (CTS)
http://cts.fiehnlab.ucdavis.edu/

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(webchem)

out <- cts_compinfo("XEFQLINVKFYRCS-UHFFFAOYSA-N", verbose = TRUE)
str(out) 
out[[1]][1:5]

inchikeys <- c("XEFQLINVKFYRCS-UHFFFAOYSA-N","BSYNRYMUTXBXSQ-UHFFFAOYSA-N" ) 
out2 <- cts_compinfo(inchikeys) 
str(out2)

sapply(out2, function(y) c(y$formula,y$molweight))


cts_convert('XEFQLINVKFYRCS-UHFFFAOYSA-N', 'inchikey', 'Chemical Name')
# 'Chemical Name','InChIKey','PubChem CID', 'ChemSpider', 'CAS'

comp <- c('XEFQLINVKFYRCS-UHFFFAOYSA-N', 'BSYNRYMUTXBXSQ-UHFFFAOYSA-N') 
cts_convert(comp, 'inchikey', 'CAS')

# cts_convert(query, from, to, first = FALSE, verbose = TRUE)

```


PubChem 
https://pubchem.ncbi.nlm.nih.gov/
CompoundID (CID) for a search query using PUG-REST https://pubchem.ncbi. nlm.nih.gov/

```{r echo=FALSE, message=FALSE, warning=FALSE}

# get_cid(query, from = "name", first = FALSE, verbose = TRUE, arg = NULL)

comp <- c('Triclosan', 'Aspirin') 
get_cid(comp)
# from = 'name'(default),'cid','sid','aid','smiles', 'inchi', 'inchikey'

pc_prop(5564, properties = NULL, verbose = TRUE)


pc_synonyms('Aspirin') 
pc_synonyms(c('Aspirin', 'Triclosan')) 
pc_synonyms(5564, from = 'cid') 
# from = 'name'(default),'cid','sid','aid','smiles', 'inchi', 'inchikey'



```


ETOX: Information System Ecotoxicology and Environmental Quality Targets 
https:// webetox.uba.de/webETOX/index.do 

```{r echo=FALSE, message=FALSE, warning=FALSE}

comps <- c('Triclosan','Glyphosate', 'DEET') 
get_etoxid(comps, match = 'all')

ids <- c("20179", "9051", "2001")
etox_basic(ids)

comp <- c('Triclosan', 'Aspirin') # nefunguje
etox_basic(comp)


id <- get_etoxid("fluazinam", match = 'best') 
etox_tests(id)
etox_basic("98976")


```


Wikidata Item ID

```{r echo=FALSE, message=FALSE, warning=FALSE}

get_wdid('Triclosan', language = 'de') 
get_wdid('DDT') 
get_wdid('DDT', match = 'all')
# match = c("best", "first", "all", "ask", "na")

comps <- c('Triclosan', 'Glyphosate') 
get_wdid(comps)

id <- c("Q408646")
wd_ident(id, verbose = TRUE)
# 'smiles', 'cas', 'cid', 'einecs', 'csid', 'inchi', 'inchikey', 'drugbank', 'zvg', 'chebi', 'chembl', 'unii' and source_url

```


Flavor percepts 
http://www.flavornet.org

```{r echo=FALSE, message=FALSE, warning=FALSE}

fn_percept("123-32-0", verbose = TRUE)

CASs <- c("75-07-0", "64-17-5", "109-66-0", "78-94-4", "78-93-3") 
fn_percept(CASs, verbose = TRUE)

```


ChemIDPlus
http://chem.sis.nlm.nih.gov/chemidplus
http://cts.fiehnlab.ucdavis.edu/

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(webchem)

ci_query('WSFSSNUMVMOOMR-UHFFFAOYSA-N', from ='inchikey')

y2 <- ci_query('WSFSSNUMVMOOMR-UHFFFAOYSA-N', from ='inchikey')
y2[[1]]$name

y2 <- ci_query('50-00-0', from ='rn')
y2[[1]]$name

comps <- c("50-00-0", "64-17-5")
ci_query(comps, from = "rn")

y2 <- ci_query('50-00-0', from = 'rn') 
y2[['50-00-0']]$inchikey
y2[['50-00-0']]$name
y2[['50-00-0']]$physprop
y2[['50-00-0']]$smiles
y2[['50-00-0']]$cas
y2[['50-00-0']]$synonyms
y2[['50-00-0']]$toxicity

y1 <- ci_query('50-00-0', from ='rn')
y1[['50-00-0']]$name
# extract log-P
sapply(y1, function(y){if (length(y) == 1 && is.na(y))
  return(NA)
  y$physprop$Value[y$physprop$`Physical Property`=='log P (octanol-water)']})

```


Query the OPSIN (Open Parser for Systematic IUPAC nomenclature) web service 
http://opsin.ch.cam.ac.uk/instructions.html

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(webchem)

opsin_query('Cyclopropane', verbose = TRUE) 

opsin_query(c('Cyclopropane', 'Octane'), verbose = TRUE) 

```


Acute toxicity data from U.S. EPA ECOTOX

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(webchem)

lc50
lc50[,1]

tnm = "67485-29-4"
nam = cir_query(tnm, 'names', first = FALSE)
lc50[which(lc50[,1]==tnm),2]

```



```{r echo=FALSE, message=FALSE, warning=FALSE}

library(PesticideLoadIndicator)

products.path()
products = products.load()
check_products_column_names(products)

substances = substances.load()

compute_pesticide_load_indicator(substances, products)


# Organic plant protection products in the river Jagst (Germany) in 2013
library(webchem)
jagst
unique(jagst[,"substance"])


library(ChemmineR)
pubchemSmilesSearch(smile)


library(CHNOSZ)
iCH <- info("SO2")
info(iCH)

```


```{r echo=FALSE, message=FALSE, warning=FALSE}

# library(rJava)
library(rcdk)

formula <- get.formula('NH4', charge = 1)
formula
formula@mass
formula@charge
formula@isotopes
formula@string


anle138b = parse.smiles("C1OC2=C(O1)C=C(C=C2)C3=CC(=NN3)C4=CC(=CC=C4)Br")[[1]]
anle138b = parse.smiles("c1ccccc1CC(=O)C(N)CC1CCCCOC1")[[1]]


get.depictor(width = 500, height = 500, zoom = 1.3, style = "cow", annotate = "off", abbr = "on", suppressh = TRUE, showTitle = FALSE, smaLimit = 100, sma = NULL) 

rcdkplot = function(molecule){
  par(mar=c(0,0,0,0)) # set margins to zero since this isn't a real plot
  temp1 = view.image.2d(molecule, depictor = NULL) # get Java representation into an image matrix. set number of pixels you want horiz and vertical
  plot(NA,NA,xlim=c(1,10),ylim=c(1,10),xaxt='n',yaxt='n',xlab='',ylab='') # create an empty plot
  rasterImage(temp1,1,1,10,10) # boundaries of raster: xmin, ymin, xmax, ymax. here i set them equal to plot boundaries
}
rcdkplot(anle138b)


todepict <- function(mol,pathsd){# raalizadeh, https://github.com/CDK-R/cdkr/issues/61
  result = tryCatch({
    factory <- .jnew("org.openscience.cdk.depict.DepictionGenerator")$withAtomColors()
    factory$withSize(1000,1000)#$getStyle("cow")
    temp1 <- paste0(pathsd)
    result<-factory$depict(mol)$writeTo(temp1)
  }, warning = function(w) {
    result=NULL
  }, error = function(e) {
    result=NULL
  })
  return(result)
}

library(MSbox)
smile <- as.character(describe('camptothecin', "smiles", info = FALSE) )
mol<-parse.smiles(smile)[[1]]

todepict(mol,'c:\\Users\\lubop\\Dropbox\\kursR\\camptothecin.png')

```





```{r}


as.character(3.14)
as.character(pi)

data <- c(-11, 21, 1.5, -31)
as.character(data)


## number of strings
length("BaCrO4")
length("How many characters?")
length(c("How", "many", "characters?"))

## number of charaters
nchar("BaCrO4")
nchar("How many characters?")
nchar(c("How", "many", "characters?"))

library(stringr)
str_length("BaCrO4")
str_length("How many characters?")
str_length(c("How", "many", "characters?"))

some_text = c("one", "two", "three", NA, "five")
str_length(some_text)
nchar(some_text)

some_factor = factor(c(1, 1, 1, 2, 2, 2), labels = c("good", "bad"))
some_factor
str_length(some_factor)
nchar(some_factor)


# Usporadani

set11 = c("today", "produced", "example", "beautiful", "a", "nicely")
# sort (decreasing order)
sort(set11)
# sort (increasing order)
sort(set11, decreasing = TRUE)


library(stringi)
stri_reverse("dna")
stri_reverse(set11)


### Spojovani retezcu

toString(17.04)
toString(c(17.04, 1978))
toString(c("Bonjour", 123, TRUE, NA, log(exp(1))))


## paste

paste("This is",  "out of",   "examples.")
paste0("This is",  "out of",  "examples.") 
paste0("This is",  " out of",  " examples.")

a <- "acetic"
b <- "acid"
c <- "anhydride"

d1 <- paste(a, b, c); d1
d2 <- paste(a, b, c, sep = "-"); d2

paste("I", "love", "R", sep = "-")


str <- paste(c(1:3), "4", sep = ":")
print (str)

str <- paste(c(1:4), c(5:8), sep = "--")
print (str)

paste("X", 1:5, sep = ".")

paste(1:3, c("!", "?", "+"))
paste(1:3, c("!", "?", "+"), sep = "")
paste0(1:3, c("!", "?", "+"))
paste(1:3, c("!", "?", "+"), sep = "", collapse = "")
paste0(1:3, c("!", "?", "+"), collapse = "")

paste("The value of pi is", round(pi,3), "and value of e is", round(exp(1),3),".")


df <- data.frame(cation=c('sodium','mercury','iron','calcium'),
                 anion=c('acetate','sulphate','dioxide','oxalate'),
                 purity=c(0.99, 0.9, 0.999, 0.985))
df
df$name = paste(df$cation, df$anion)
df

library(stringr)
str_c("May", "The", "Force", "Be", "With", "You")


library(stringr)
str_c("May", "The", "Force", NULL, "Be", "With", "You", character(0))
paste("May", "The", "Force", NULL, "Be", "With", "You", character(0))

str_c("May", "The", "Force", "Be", "With", "You", sep = "_")

# str_join("May", "The", "Force", "Be", "With", "You", sep = "-") 
# analogicka fce k predchozi


## cat - pouze vypisuje retezec

cat("learn", "code", "tech", sep = ":")
str <- cat("learn", "code", "tech", sep = ":") # nelze ulozit do promenne
print (str) 

my_string = c("learn", "code", "tech")
cat(my_string, "with R")
cat(my_string, "with R", sep = " =) ")

cat(1:10, sep = "-")
cat(month.name[1:4], sep = " ")

cat(c(1:5), file ='sample.txt')
cat(my_string, "with R", file = "output.txt")
cat(11:20, sep = '\n', file = 'temp.csv')
readLines('temp.csv') # read the file temp.csv

cat("The value of pi is", round(pi,3), "and value of e is", round(exp(1),3),".")


## print - vypisuje vsechny formaty

print(df)
df

my_value <- 8235.675324 
my_value
print(my_value) 
print(my_value, digits = 5) 

print(c("learn", "code", "tech"))
paste(c("learn", "code", "tech"))
paste(c("learn", "code", "tech"), collapse=" ")
cat(c("learn", "code", "tech"))

my_string <- "This is \nthe example string" 
print(my_string)  
cat(my_string) 


### 

library(stringr)

str_dup("hola", 3)
str_dup("adios", 1:3)

words = c("lorem", "ipsum", "dolor", "sit", "amet")
str_dup(words, 2)
str_dup(words, 1:5)


####

format(c("A", "BB", "CCC"), width = 5, justify = "centre")
format(c("A", "BB", "CCC"), width = 5, justify = "left")
format(c("A", "BB", "CCC"), width = 5, justify = "right")
format(c("A", "BB", "CCC"), width = 5, justify = "none")

library(stringr)
str_pad("hola", width = 7)
str_pad("adios", width = 7, side = "both")
str_pad("hashtag", width = 8, pad = "#")
str_pad("hashtag", width = 9, side = "both", pad = "-")


##### editace ciselnych retezcu

## sprintf

# '%f' indicates 'fixed point' decimal notation
sprintf("%f", pi)
sprintf("%f", 123.45)
sprintf("%f", 123.456789)
# decimal notation with 3 decimal digits
sprintf("%.3f", pi)
sprintf("%.3f", 123.456789)
# 1 integer and 0 decimal digits
sprintf("%1.0f", pi)
sprintf("%1.0f", 123.456789)
# decimal notation with 3 decimal digits
sprintf("%5.1f", pi)
sprintf("%05.1f", pi)
sprintf("%6.1f", 123.456789)
sprintf("%06.1f", 123.456789)
# print with sign (positive)
sprintf("%+f", pi)
sprintf("%+f", 123.456789)
# prefix a space
sprintf("% f", pi)
sprintf("% f", 123.456789)
# left adjustment
sprintf("%-10f", pi)
sprintf("%-10f", 123.456)
# exponential decimal notation 'e'
sprintf("%e", pi)
sprintf("%e", 123.456789)
# exponential decimal notation 'E'
sprintf("%E", pi)
sprintf("%E", 123.456789)
# number of significant digits (6 by default)
sprintf("%g", pi)
sprintf("%g", 123.456789)


## format

format(pi)
format(123.45)
format(123.456789)

format(123.45, nsmall = 5)
format(12.3456789, nsmall=2)
format(12.3456789, nsmall=7)
format(12.3, nsmall=3)

format(12.3456789, digits=2)
format(12.3456789, digits=5)
format(c(6, 13.1), digits = 2)

format(12.3456789, digits=5, nsmall = 6)
format(c(6, 13.1), digits = 2, nsmall = 2)

format(1/1:5, digits = 2)
format(format(1/1:5, digits = 2), width = 6, justify = "centre")


format(12345678, big.mark = ",") # oddeleni tisicu
format(1230, big.mark = ",")


## uvozovky ve vypisu retezce

my_string = "programming with data is fun"
print(my_string)
print(my_string, quote = FALSE)

noquote(my_string)

no_quotes = noquote(c("some", "quoted", "text", "!%^(&="))
no_quotes
no_quotes[2:3]

noquote(paste("I", "love", "R", sep = "-"))
noquote(cat("The value of pi is", round(pi,3), "and value of e is", round(exp(1),3),".")
)


dQuote(my_string)
sQuote(my_string)

x <- "2020-05-29 19:18:05"
dQuote(x)
sQuote(x)


########

substr(month.name, 1, 3)

substring(month.name, 1, 3)

strtrim(month.name, 3)


rs <- ("This is First R String Example")
strsplit(rs, split = " ")

rs <- ("This&is&First&R&String&Example")
strsplit(rs, split = "&")

a <- "Alabama-Alaska-Arizona-Arkansas-California"
strsplit(a, split = "-")

str = "Splitting sentence into words"
strsplit(str, " ")
unlist(strsplit(str, " "))


rs <- ("C21H22N2O2") # strychnin

strsplit(rs, split = "[0-9]+")[[1]]

strsplit(rs, split = "[A-Z]+")[[1]][-1]
strsplit(rs, "\\D+")[[1]][-1]
regmatches(rs, gregexpr("[[:digit:]]+", rs))
library(stringr)
str_extract_all(rs,"[0-9]+")[[1]]
library(strex)
str_extract_numbers(rs,decimals = FALSE)

rs <- ("C21H22N2O2")
strsplit(rs, split = "")

string_date<-c("2-07-2020","5-07-2020","6-07-2020",
               "7-07-2020","8-07-2020")
ssp = strsplit(string_date,split = "-"); ssp


# extract 'bcd'
substr("abcdef", start=2, stop=4)

substring("ABCDEF", 2, 4)

# extract each letter
substring("ABCDEF", 1:6, 1:6)


library(stringr)
lorem = "Lorem Ipsum"
substring(lorem, first = 1, last = 5)
str_sub(lorem, start = 1, end = 5)

str_sub("adios", 1:3)

resto = c("brasserie", "bistrot", "creperie", "bouchon")
substring(resto, first = -4, last = -1)
str_sub(resto, start = -4, end = -1)

str_sub(lorem, seq_len(nchar(lorem)))
substring(lorem, seq_len(nchar(lorem)))

str_sub(lorem,-2) 


### orezani mezer na koncich retezce

trimws(" This has trailing spaces.  ")

bad_text = c("This", " example ", "has several   ", "   whitespaces ")
trimws(bad_text)
library(stringr)
str_trim(bad_text, side = "left")
str_trim(bad_text, side = "right")
str_trim(bad_text, side = "both")


## replace
chartr(old = "All", new = "aLL", "All ChaRacterS in Upper Case") 
chartr("a", "A", "This is a boring string")
chartr("a", "0", "This is a bad example")

crazy = c("Here's to the crazy ones", "The misfits", "The rebels")
chartr("aei", "#!?", crazy)


x = c("may", "the", "force", "be", "with", "you")
substr(x, 2, 2) <- "#"
x
y = c("may", "the", "force", "be", "with", "you")
substr(y, 2, 3) <- ":)"
y

z = c("may", "the", "force", "be", "with", "you")
substr(z, 2, 3) <- c("#", "@")
z

text = c("more", "emotions", "are", "better", "than", "less")
substring(text, 1:3) <- c(" ", "zzz")
text


library(stringr)
lorem = "Lorem Ipsum"
str_sub(lorem, -1) <- "Nullam"
lorem

lorem = "Lorem Ipsum"
str_sub(lorem, 1, 5) <- "Nullam"
lorem

lorem = "Lorem Ipsum"
str_sub(lorem, c(1, 7), c(5, 8)) <- c("Nullam", "Enim")
lorem


## to lower case
tolower("BaCrO4")
tolower(c("aLL ChaRacterS in LoweR caSe", "ABCDE"))
casefold("aLL ChaRacterS in LoweR caSe")

## to upper case
toupper("BaCrO4")
toupper(c("All ChaRacterS in Upper Case", "abcde"))
casefold("All ChaRacterS in Upper Case", upper = TRUE)


### filtrovani

startsWith(month.name, "J")
endsWith(month.name, "ember")

myStrings <- paste(1:3, month.name, sep = ". ")
myStrings
# Is a pattern present (returns a logical vector)?
grepl("ember", myStrings)
# In which elements is a pattern present (returns indices)?
grep("ember", myStrings)
# In which elements is a pattern present (returns the values)?
grep("ember", myStrings, value = TRUE)

x1<-c("R is a programming language and programming software environment",
      "R is freely available under the GNU General Public License",
      "This programming language was named R")
grep("programming",x1,fixed=TRUE)


# Srovnani 2 retezcu
message1 <- "Pro"
message2 <- "Pro"
message3 <- "pRO"
message1 == message2
message1 == message3

# Hledani retezce v jinem retezci
str <- "Hello World!"
grepl("H", str)
grepl("Hello", str)
grepl("X", str)
grepl(message1, message2)
grepl(message1, message3)

identical(message1, message2)
identical(message1, message3)
identical(tolower(message1), tolower(message3))

message1[message1 %in% message2]
message1[message1 %in% message3]  
message1[tolower(message1) %in% tolower(message3)] 

vector1 <- c("hey", "hello", "greetings")
vector2 <- c("hey", "hello", "hi")
vector1[vector1 %in% vector2]


###

library(stringr)

change = c("Be the change", "you want to be")
# extract first word
word(change, 1)
# extract second word
word(change, 2)
# extract last word
word(change, -1)
# extract all but the first words
word(change, 2, -1)

word(change,start=1,end=2,sep=fixed(" "))


data <-  c('Ab_Cd-001234.txt','Ab_Cd-001234.txt')

gsub('.*-([0-9]+).*','\\1','Ab_Cd-001234.txt')
x <- c('Ab_Cd-001234.txt','Ab_Cd-001234.txt')
sub('.*-([0-9]+).*','\\1',x)

library(stringr)
regexp <- "[[:digit:]]+"
str_extract(data, regexp)

library(qdap)
genXtract("Ab_Cd-001234.txt", "-", ".txt")
x <- c('Ab_Cd-001234.txt','Ab_Cd-001234.txt')
genXtract(x, "-", ".txt")

library(tools)
sub(".*-", "", file_path_sans_ext(x))

library(gsubfn)
strapplyc(x, "-(\\d+)\\.", simplify = TRUE)
strapply(x, "-(\\d+)\\.", as.numeric, simplify = TRUE)


x <- c("a very nice character string")  
library(stringr)
str_replace(x, "c", "xxx")
str_replace_all(x, "c", "xxx")

x <- "aaabbb" 
sub("a", "c", x)  
gsub("a", "c", x)
sub("a|b", "c", x)
gsub("a|b", "c", x)  

x <- c("d", "a", "c", "abba") 
grep("a", x)
grepl("a", x)
grep("a|c", x)
grepl("a|c", x)

regexpr("a", x)
gregexpr("a", x)
regexec("a", x) 

x <- "example_xxx_string"
library(stringr)
str_sub(string = x, start = 8, end = 12)
str_sub(string = x, start = 8, end = 12) <- " character "  # Replace substring
x  

x <- "xxxxyxxyxaaaaaay"
x 
sub("y", "NEW", x) 
gsub("y", "NEW", x)

library(stringr)
str_replace(x, "y", "NEW")
str_replace_all(x, "y", "NEW")

```


Kvadraticke a kubicke rovnice


Velikost výslednice dvou navzajem kolmých sil je 34 N. Jaká je velikost skládaných sil, je-li jedna z nich o 14 N větší než druha?

```{r echo=FALSE, message=FALSE, warning=FALSE}

# Pythagorova veta: x^2 + (x + 14)^2 = 34^2      
library(Ryacas)
eq <- "x^2 + (x + 14)^2 - 34^2" 
yac_str(paste0("Simplify(", eq, ")"))  # zjednodusení výrazu

# diskriminant
dc = c(1,14,-480)
D = dc[2]^2 - 4*dc[1]*dc[3]
if (D == 0) {cat("Kvadraticka rovnice ma dva sobe rovne realne koreny (dvojnasobny koren).")}
if (D > 0) {cat("Kvadraticka rovnice ma dva ruzne realne koreny.")}
if (D < 0) {cat("Kvadraticka rovnice nema zadny koren v oboru realnych cisel. V oboru komplexnich cisel ma dva imaginarni komplexne sdruzene koreny.")}

# Descartes rule
sum(diff(sign(dc)) != 0)
cat("Pocet kladnych korenu:", sum(diff(sign(dc)) != 0)) # pocet kladnych korenu

library(sfsmisc)
nr.sign.chg(dc)
cat("Pocet kladnych korenu:", nr.sign.chg(dc)) # pocet kladnych korenu

# Hornerovo schema
library(sfsmisc)
polyn.eval(dc, x)


fun <- function (x) {x^2 + (x + 14)^2 - 34^2}
uniroot(fun, c(-1, 0),extendInt = "yes", tol = 1e-9)  # base
uniroot(fun, c(0, 1),extendInt = "yes", tol = 1e-9)
# uniroot(fun, c(-1, 1),extendInt = "yes", tol = 1e-9)
F1 <- uniroot(fun, c(0, 100),extendInt = "yes", tol = 1e-9) 
F1$root

fun <- function(x){x^2 + 14*x - 480}
uniroot(fun, c(-1, 0),extendInt = "yes", tol = 1e-9)  # base
uniroot(fun, c(0, 1),extendInt = "yes", tol = 1e-9)
F1 <- uniroot(fun, c(0, 100),extendInt = "yes", tol = 1e-9)
F1$root 


rr = polyroot(c(-480, 14, 1)) # base
Re(rr)


library(pracma)
rr = roots(c(1, 14, -480))
rr
F1 = rr[rr>=0]
F1
F2 = rr[rr<=0]
F2

# Vietovy vzorce
uniroot(fun, c(-1, 0),extendInt = "yes", tol = 1e-9)$root
(-dc[2]/dc[1]) - F1
(dc[3]/dc[1])/F1

# koeficienty kvadraticke rovnice z korenu
library(pracma)
Poly(c(16, -30))


## graficke reseni
xx = seq(-50,50,1)
yy = abs(xx^2 + (xx + 14)^2 - 34^2)
plot(xx,yy,type="l")
abline(h=0,col=2)
xx[which(yy==0)]
xx[yy==0]

## graficke reseni
xx = seq(-50,50,1)
yy = xx^2
zz = -14*xx + 480
plot(xx,yy,type="l")
abline(480,-14,col=2)

plot(xx,yy,type="l",xlim = c(10,20),ylim=c(0,500))
abline(480,-14,col=2)


```


Tlaková láhev s oxidem uhličitým obsahuje 10.0 kg plynu. Jaký objem zaujímá stlačený plyn, když při teplotě 30 °C je tlak v láhvi 13.17e-6 Pa? Výpočet proved'te pomocí van der Waalsovy rovnice.

```{r echo=FALSE, message=FALSE, warning=FALSE}

# (p + a/Vm^2)(Vm - b) = R*T     
# [Vm = 0.075 m3/kmol, n = 0.2273 kmol, V = 0.0171 m3]

library(measurements)
p = 13.17e6 # [Pa]
m = conv_unit(10.0, from="kg", to="g")
t = conv_unit(30, from="C", to="K") 
R = 8.3141  # [J / mol K]
R = R*1000  # [J / kmol K]
Mr = 44.01  # [g/mol]
n = m/Mr    # [mol]
n = n/100   # [kmol]

# vypocet pro idealni plyn
# p*Vm - R*T = 0  
Vm = R*t/p
V = Vm*n; V  # [m3]

# vypocet pro realny plyn
a = 0.365e6 # [m6 Pa / kmol2]
b = 0.0428  # [m3 / kmol]
# Vm**3 - (b + R*t/p)*Vm**2 + (a/p)*Vm - a*b/p 

# realne i komplexni koreny
library(RConics) 
CE = c(1,-(b + R*t/p),a/p,-(a*b/p))
x0 = cubic(CE); x0
Vm = Re(x0[which(Im(x0)==0)])  # [m3 / kmol]
Vm
V = Vm*n; V  # [m3]

# jen realne koreny
fun <- function (x) CE[1]*x^3 + CE[2]*x^2 + CE[3]*x + CE[4]
#curve(fun(x),-4,4)
#abline(h = 0, lty = 3)
Vm <- uniroot(fun, c(-4, 4))$root  # base
V = Vm*n; V  # [m3]

library(rootSolve)
Vm <- rootSolve::uniroot.all(fun, c(-4, 4))
V = Vm*n; V  # [m3]

library(Rmpfr)
Vm <- Rmpfr::unirootR(fun, lower=-4, upper=4, tol = 1e-9)$root
V = Vm*n; V  # [m3]

## graficke reseni
# CE = c(1,-(b + R*t/p),a/p,-(a*b/p))
xx <- seq(-4,4, by=0.01)
yy <- CE[1]*xx^3 
zz <- -CE[2]*xx^2 - CE[3]*xx - CE[4]
plot(xx,yy,type="l", col=2)
points(xx,zz,type="l", col=4)
plot(xx, yy, type="l",col=2,xlim=c(0.05,0.1),ylim=c(0,0.001))
points(xx, zz, type="l",col=4)

```


Jake pH má roztok kyseliny mravenčí o koncentraci 8.5 x 10-4 mol/l?

```{r echo=FALSE, message=FALSE, warning=FALSE}

pKA = 3.752
KA = 10^-pKA; KA
Kw = 10^-14
cA = 8.5e-4 # [mol/l ] 


# H = sqrt(KA*cA)
H <- sqrt(KA*cA)
pH = -log10(H); pH 


# [H+]^2 + KA*[H+] + KA*cA = 0
CE = c(1,KA,-KA*cA)
fun <- function (x) CE[1]*x^2 + CE[2]*x + CE[3]
uniroot(fun, c(-1e-1, 0),tol = 1e-9)
uniroot(fun, c(0, 1e-1),tol = 1e-9)
H <- uniroot(fun, c(0, 1e-1),tol = 1e-9)$root  # base
pH = -log10(H); pH

library(rootSolve)
rootSolve::uniroot.all(fun, c(-1e-1, 0), tol = 1e-9)
rootSolve::uniroot.all(fun, c(0, 1e-1), tol = 1e-9)
H <- rootSolve::uniroot.all(fun, c(-1e-1, 1e-1), tol = 1e-9)
pH = -log10(H); pH 

library(Rmpfr)
Rmpfr::unirootR(fun, lower=-1e-1, upper=0, tol = 1e-9)
Rmpfr::unirootR(fun, lower=0, upper=1e-1, tol = 1e-9)
H <- Rmpfr::unirootR(fun, lower=0, upper=1e-1, tol = 1e-9)$root
pH = -log10(H); pH 


# [H+]^3 + KA*[H+]^2 - [H+]*(KA*cA + Kw) - KA*Kw = 0
CE = c(1,KA,-(KA*cA + Kw),-KA*Kw)

library(RConics) 
x0 = cubic(CE); x0
H = Re(x0[which(Im(x0)==0)])
H = Re(x0[which(Re(x0)>=0)])
pH = -log10(H); pH 


fun <- function (x) CE[1]*x^3 + CE[2]*x^2 + CE[3]*x + CE[4]
uniroot(fun, c(-1e-1,0),tol = 1e-9)
uniroot(fun, c(0, 1e-1),tol = 1e-9)
H <- uniroot(fun, c(0, 1e-1),tol = 1e-9)$root  # base
pH = -log10(H); pH 


library(rootSolve)
rootSolve::uniroot.all(fun, c(-1e-1, 0),tol = 1e-9)
rootSolve::uniroot.all(fun, c(0,1e-1),tol = 1e-9)
H <- rootSolve::uniroot.all(fun, c(-1e-1, 1e-1),tol = 1e-9)
pH = -log10(H); pH 


library(Rmpfr)
Rmpfr::unirootR(fun, lower=-1e-1, upper=0, tol = 1e-9)
Rmpfr::unirootR(fun, lower=0, upper=1e-1, tol = 1e-9)
H <- Rmpfr::unirootR(fun, lower=0, upper=1e-1, tol = 1e-9)$root
pH = -log10(H); pH 

```



Soustavy lineárních rovnic


Ze dvou slitin s 60% a 80% obsahem mědi se ma získat 40 kg slitiny se 75% obsahem mědi. Kolik kg každé slitiny je třeba použít?  [10 a 30 kg]

```{r}

library(rootSolve)
model <- function(x){
  F1 <- 0.6*x[1] + 0.8*x[2] - 30
  F2 <- x[1] + x[2] - 40
  c(F1 = F1, F2 = F2)}
ss <- multiroot(f = model, start = c(1, 1))
ss$root

### graficke reseni
# 0.6*m1 + 0.8*m2 = 40*0.75 => m1 = 40*0.75/0.6 - 0.8/0.6 * m2  => m1 = 50 - 1.33 * m2 
# m1 + m2 = 40  => m1 = 40 - m2
xx = seq(0,50,0.1)
yy = 50 - 1.33*xx
zz = 40 - xx
plot(xx, yy, type="l",col=2)
points(xx, zz, type="l",col=4)
plot(xx, yy, type="l",col=2,xlim=c(25,35),ylim=c(0,20))
points(xx, zz, type="l",col=4)

```


Do bazenu natece prutokem A za 3 h a prutokem B za 4 h celkem 2150 hl vody. Prutokem A za 4 h a prutokem B za 2 h by nateklo 1700 hl vody. Kolik hl vody natece prutokem A a kolik prutokem B za 1 hodinu?      A 250 hl, B 350 hl

```{r}

library(rootSolve)
model <- function(x){
  F1 <- 3*x[1] + 4*x[2] - 2150
  F2 <- 4*x[1] + 2*x[2] - 1700
  c(F1 = F1, F2 = F2)}
ss <- multiroot(f = model, start = c(1, 1))
ss$root

### graficke reseni
# 3*m1 + 4*m2 = 2150 => m2 = 2150/4 - 3/4 * m1  => m2 = 537.5 - 0.75 * m1
# 4*m1 + 2*m2 = 1700 => m2 = 1700/2 - 4/2 * m1 => m2 = 850 - 2 * m1
xx = seq(0,500,1)
yy = 537.5 - 0.75*xx
zz = 850 - 2*xx
plot(xx, yy, type="l",col=2)
points(xx, zz, type="l",col=4)
plot(xx, yy, type="l",col=2,xlim=c(200,300),ylim=c(300,400))
points(xx, zz, type="l",col=4)

```





```{r echo=FALSE, message=FALSE, warning=FALSE}

### skalarni soucin (inner product, dot product, scalar product)

u <- rep(3,3)
v <- 1:3
u%*%v # the inner product
library(pracma)
dot(u, v)
dot(c(1, 2, 3), c(4, 5, 6))


### vektorovy soucin (outer product, cross product, vector product)

u <- rep(3,3)
v <- 1:3
u%o%v # The outer product
library(pracma)
cross(u, v)
cross(c(1, 2, 3), c(4, 5, 6))  
A <- matrix(c(1,0,0, 0,1,0), nrow=2, ncol=3, byrow=TRUE)
crossn(A)  


### Kartezsky soucin

expand.grid(u, v)
expand.grid(c(1, 2, 3), c(4, 5, 6))
expand.grid(x=c("a","b","c"),y=c(1,2,3))
library(data.table)
CJ(u, v)
CJ(c(1, 2, 3), c(4, 5, 6))
CJ(x=c("a","b","c"),y=c(1,2,3))
library(tidyr)
x <- data.frame(x=c("a","b","c"))
y <- data.frame(y=c(1,2,3))
crossing(x, y)

```



```{r}

M <- matrix(data = rnorm(12), ncol = 3)

# transponovana matice
t(M) 

dim(M)
nrow(M) 
ncol(M)
length(M)
library(Matrix)
nnzero(M, na.counted = NA) # pocet nenulovych prvku matice

round(M,2)
library(Rfast)
Round(M,digit=2,na.rm = FALSE)


### prevod symbolicke matice na numerickou
mch = matrix("1",3,3)
mch
as.numeric(mch) # nefunguje
matrix(as.numeric(mch),ncol=ncol(mch))
library(gmp)
asNumeric(mch)


### Premena matice na vektor
library(gdata)
unmatrix(M, byrow=FALSE)

### Matice na x, y, z
library(squash)
xyzmat2xyz(M)
expand.grid(c(1, 2, 3), c(4, 5, 6))


library(Rfast2) # Split the matrix in lower,upper triangular and diagonal
lud(M)


# jednotkova matice, identita
diag(3)
MM = diag(x = 1, nrow = 3, ncol = 3, names = TRUE)
diag(MM) <- 2
library(Rfast)
Diag.matrix(3,v=1)
Diag.fill(MM,v=0) 
library(Matrix)
Diagonal(3, x = 1)

upper.tri(matrix(1, 3, 3))
round(upper.tri(matrix(1, 3, 3)))
upper.tri(MM)
MM[upper.tri(MM)] <- 0
lower.tri(matrix(1, 3, 3))
round(lower.tri(matrix(1, 3, 3)))
lower.tri(MM)
MM[lower.tri(MM)] <- 0
library(gdata)
upperTriangle(MM, diag=FALSE, byrow=FALSE)
upperTriangle(MM, diag=TRUE, byrow=FALSE)
upperTriangle(MM, diag=FALSE, byrow=FALSE) <- 1
MM
lowerTriangle(MM, diag=FALSE, byrow=FALSE) 
lowerTriangle(MM, diag=FALSE, byrow=FALSE) <- 3
MM
library(Rfast)
lower_tri(MM, suma = FALSE, diag = FALSE) 
upper_tri(MM, suma = FALSE, diag = FALSE) 
lower_tri.assign(MM, 1, diag = FALSE) 
upper_tri.assign(MM, 3, diag = FALSE)
library(Matrix) # v maticovem tvaru
tril(MM) # lower triangular
triu(MM) # upper triangular

# LU decomposition (decompose a matrix into lower- and upper-triangular matrices)
library(cmna)
lumatrix(MM)

### nasobeni matic
b <- matrix(nrow = 2, ncol = 2, c(1, 2, 3, 4)); b
a <- matrix(nrow = 2, ncol = 2, c(1, 0, 0, -1)); a
a%*%b
b%*%a
library(Rfast)
mat.mult(a, b)


### stopa (trace) - soucet prvku hlavni diagonaly
library(fBasics)
tr(MM)


## determinant
det(M[,c(1:3)])

library(matlib)
Det(M[,c(1:3)], method = "elimination", verbose = FALSE, fractions = FALSE)
# method = c("elimination", "eigenvalues", "cofactors")


### je matice pozitivne definitni
library(fBasics)
isPositiveDefinite(MM) 
makePositiveDefinite(MM)


### hodnost matice

# Frobeniova veta: soustava linearnich rovnic ma reseni tehdy a jen tehdy, maji-li matice a rozsirena matice soustavy stejnou hodnost.
qr(M)$rank

library(Matrix)
rankMatrix(M, tol = NULL, method = "qr", sval = svd(M, 0, 0)$d, warn.t = TRUE)[1]
# method = c("tolNorm2", "qr.R", "qrLINPACK", "qr", "useGrad", "maybeGrad")

library(fBasics)
rk(M, method = "qr")
# method = c("qr", "chol")

library(matrixcalc) # JEN PRO CTVERCOVOU MATICI
matrix.rank(M[,c(1:3)], method = "qr")
# method = c("qr", "chol")

library(limSolve)
resolution(M)  # $nsolvable = hodnost matice


### Inverzni matice

## ctvercova matice
M2 <- cbind(c(1,0,1),c(0,1,2),c(0,0,1))
solve(M2)
solve(M2)%*%M2  # check

library(fBasics)
inv(M2)
inv(M2)%*%M2  # check

library(cmna)
invmatrix(M2)
invmatrix(M2)%*%M2  # check

library(matlib)
Inverse(M2, verbose=FALSE, fractions=FALSE)
Inverse(M2)%*%M2 # check

## zobecnena inverze
library(MASS)  #  Moore-Penrose generalized inverse 
ginv(M2)
round(ginv(M2),0)
ginv(M2)%*%M2   # check
round(ginv(M2)%*%M2,0)

library(matlib) 
Ginv(M2, fractions=FALSE) 
Ginv(M2) %*% M2 # check

library(limSolve)
limSolve::Solve(M2)
Solve(M2) %*% M2 # check


### Eigenvalues, eigenvectors
eigen(M2)

library(cmna)
library(geigen)

library(matlib)
Eigen(M2,tol = sqrt(.Machine$double.eps), max.iter = 100, retain.zeroes = TRUE)
C <- matrix(c(1,2,3,2,5,6,3,6,10), 3, 3) # nonsingular, symmetric C 
EC <- Eigen(C) # eigenanalysis of C 
EC$vectors %*% diag(EC$values) %*% t(EC$vectors) # check




```


Sestavte matici konstitučních koeficientů

```{r}

form =c("CH4","H2O","CO2","CO")
form =c("CO","COS","CH4","S2","H2O")
form = c("KMnO4", "MnSO4", "H2O", "MnO2", "K2SO4", "H2SO4")
form = c("Fe2O3", "Al", "Al2O3", "Fe")

library(CHNOSZ)
chars = lapply(form,function(x){names(count.elements(x))})
nums = lapply(form,function(x){as.vector(count.elements(x))})
ulch = sort(unique(unlist(chars)))
Mrl = as.data.frame(matrix(0,length(chars),length(ulch)))
colnames(Mrl)=ulch
for(X in c(1:length(chars))){  
  for(i in c(1:length(chars[[X]]))){ 
    cn1 = which(ulch == chars[[X]][i])
    cn2 = which(chars[[X]] == chars[[X]][i])
    Mrl[X,cn1] = as.numeric(nums[[X]][cn2])
  }
}
rownames(Mrl)=form
Mrl

```


Určete počet lineárně nezávislých reakcí v soustavě obsahující dané složky.

Gibbsovo stechiometricke pravidlo: maximální počet lineárně nezávislých reakcí r je menší nebo roven rozdílu počtu složek v systému (n) a hodnosti matice konstitučních koeficientů (h), v n-složkovém systému existuje n - h stechiometrických vztahů.

```{r}

form = c("CH4","H2O","CO","CO2","H2")
form = c("NH3","N2","NO","NO2","H2O","O2")

library(CHNOSZ)
chars = lapply(form,function(x){names(count.elements(x))})
nums = lapply(form,function(x){as.vector(count.elements(x))})
ulch = sort(unique(unlist(chars)))
Mrl = as.data.frame(matrix(0,length(chars),length(ulch)))
colnames(Mrl)=ulch
for(X in c(1:length(chars))){  
  for(i in c(1:length(chars[[X]]))){ 
    cn1 = which(ulch == chars[[X]][i])
    cn2 = which(chars[[X]] == chars[[X]][i])
    Mrl[X,cn1] = as.numeric(nums[[X]][cn2])
  }
}
rownames(Mrl)=form
Mrl
M = Mrl

# hodnost matice
qr(M)$rank
library(Matrix)
rankMatrix(M)[1]

# pocet linearne nezavisl}ch reakcm
(r = length(form) - qr(M)$rank)

```


```{r}

### Vycislete rovnici: Fe2O3 + Al = Al2O3 + Fe
Fe = c(2, 0, 0, -1)
O = c(3, 0, -3, 0)
Al = c(0, 1, -2, 0)
A = as.matrix(rbind(Fe, O, Al))
colnames(A) = c("Fe2O3", "Al", "Al2O3", "Fe")
B = c(rep(0,length(A[,1])))
C = cbind(A,B)
library(matlib)
c(R(A), R(cbind(A,B)))          # show ranks
all.equal(R(A), R(cbind(A,B)))  # consistent?
showEqn(A, B)
matlib::Solve(A, B)
GE1 = gaussianElimination(A, B, tol = sqrt(.Machine$double.eps), verbose = FALSE, latex = FALSE, fractions = TRUE)
as.character(GE1)
# cim nasobit koeficienty
nnv1 = 1/min(abs(as.numeric(GE1[,length(A[1,])])))
nnv = as.numeric(GE1[,length(A[1,])])*nnv1
NN = -1*(as.matrix(GE1)*nnv)[,c(1:(length(A[1,])-1))]
rownames(NN) = rownames(A)
colnames(NN) = colnames(A)[-length(colnames(A))]
apply(NN,2,sum)

library(MASS)
nn = length(A[1,])
a = A[,-nn]
b = -A[,nn]
ge1 <- MASS::ginv(a) %*% b
ge1*2

library(limSolve)
Ls = limSolve::Solve(A, B = diag(nrow = nrow(A)), tol = sqrt(.Machine$double.eps))
library(MASS)
fractions(Ls)
fractions(Ls)*30

```


Vyčíslete rovnici: KMnO4 + MnSO4 + H2O = MnO2 +  K2SO4 + H2SO4

```{r}

K = c(1, 0, 0, 0, -2, 0)
Mn = c(0, 1, 0, -1, 0, 0)
O = c(4, 4, 1, -2, -4, -4)
S = c(0, 1, 0, 0, -1, -1)
H = c(0, 0, 2, 0, 0, -2)

A = as.matrix(rbind(K, Mn, O, S, H))
colnames(A) = c("KMnO4", "MnSO4", "H2O", "MnO2", "K2SO4", "H2SO4")
B = c(rep(0,length(A[,1])))
C = cbind(A,B)
library(matlib)
c(R(A), R(cbind(A,B)))          # show ranks
all.equal(R(A), R(cbind(A,B)))  # consistent?
showEqn(A, B)
matlib::Solve(A, B)
GE2 = gaussianElimination(A, B, tol = sqrt(.Machine$double.eps), verbose = FALSE, latex = FALSE, fractions = TRUE)
as.character(GE2)
# cim nasobit koeficienty
nnv1 = 1/min(abs(as.numeric(GE2[,length(A[1,])])))
nnv = as.numeric(GE2[,length(A[1,])])*nnv1
NN = -1*(as.matrix(GE2)*nnv)[,c(1:(length(A[1,])-1))]
rownames(NN) = rownames(A)
colnames(NN) = colnames(A)[-length(colnames(A))]
apply(NN,2,sum)

```


Kolik 96% kyseliny sírové a kolik vody je potřeba na pripravu 1 l 78% kyseliny sírové?

```{r}

r1 = c(0, 1)
r2 = c(0.96, 0.04)
A = cbind(r1, r2); rownames(A) = c("H2SO4","H2O")
B = c(0.78, 0.22); names(B) = c("H2SO4","H2O")
library(matlib)
c(R(A), R(cbind(A,B)))          # show ranks, viz Frobeniova veta
all.equal(R(A), R(cbind(A,B)))  # consistent?
showEqn(A, B)
matlib::Solve(A, B)
limSolve::Solve(A, B)
GE1 = gaussianElimination(A, B, tol = sqrt(.Machine$double.eps), verbose = FALSE, latex = FALSE, fractions = FALSE)
NN = GE1[,3] # [l]
names(NN) = rownames(A)
NN
library(Rlinsolve)
ls = lsolve.gs(A, B, xinit = NA, reltol = 1e-05, maxiter = 1000, adjsym = TRUE, verbose = TRUE)
ls$x
library(cmna)
cgmmatrix(A, B, tol = 1e-06, maxiter = 100) # iterativematrix
solvematrix(A, B) # refmatrix

```


Slitina A obsahuje 1.5 % Si, 1.4 % Mn, 0.4 % P a 0.3 % S. Slitina B obsahuje 0.5 % Si, 1.6 % Mn, 0.2 % P a 0.2 % S. Slitina C obsahuje 3 % Si, 0.5 % Mn, 0.5 % P a 0.05 % S. Kolik každé slitiny A, B, C je třeba na výrobu 100 kg slitiny obsahující 2 % Si, 1 % Mn, 0.4 % P a 0.15 % S?

```{r}

sl1 = c(0.015, 0.014, 0.004, 0.003)
sl2 = c(0.005,0.016, 0.002, 0.002)
sl3 = c(0.03, 0.005, 0.005, 0.0005)
A = cbind(sl1, sl2, sl3); rownames(A) = c("Si","Mn","P","S")
B = c(0.02, 0.01, 0.004, 0.0015); names(B) = c("Si","Mn","P","S")
# B = c(2, 1, 0.4, 0.15); names(B) = c("Si","Mn","P","S")
library(matlib)
c(R(A), R(cbind(A,B)))          # show ranks, viz Frobeniova veta
all.equal(R(A), R(cbind(A,B)))  # consistent?
showEqn(A, B)
matlib::Solve(A, B)
limSolve::Solve(A, B)
GE1 = gaussianElimination(A, B, tol = sqrt(.Machine$double.eps), verbose = FALSE, latex = FALSE, fractions = FALSE)
GE1[,4]*100 # [kg]
library(Rlinsolve)
ls = lsolve.gs(A, B, xinit = NA, reltol = 1e-05, maxiter = 1000, adjsym = TRUE, verbose = TRUE)
ls$x

```


Kolik g 60% a kolik g 30% roztoku NaCl je treba smichat pri priprave 100 g 40% roztoku? 20 g 60% a a 80 g 35%

```{r}

library(rootSolve)
model <- function(x){
  F1 <- 0.6*x[1] + 0.3*x[2] - 40
  F2 <- x[1] + x[2] - 100
  c(F1 = F1, F2 = F2)}
ss <- multiroot(f = model, start = c(1, 1))
ss$root

### graficke reseni
# 0.6*m1 + 0.3*m2 = 100*0.4 => m1 = 40/0.6 - 0.3/0.6 * m2  => m1 = 66.67 - 0.5 * m2 
# m1 + m2 = 100  => m1 = 100 - m2
xx = seq(0,100,0.1)
yy = 66.67 - 0.5*xx
zz = 100 - xx
plot(xx, yy, type="l",col=2)
points(xx, zz, type="l",col=4)
plot(xx, yy, type="l",col=2,xlim=c(65,70),ylim=c(30,40))
points(xx, zz, type="l",col=4)

### matice

B = c(0.40*100,100) # [mg]
names(B) = c("60%","30%")
r1 = c(0.60,1) # [mg]
r2 = c(0.30,1) # [mg]
A = cbind(r1,r2)
colnames(A) = c("60%","30%")
rownames(A) = c("r30","r3")
det(A) # matice je regularni n = h
library(matlib)
c(R(A), R(cbind(A,B)))          # show ranks
all.equal(R(A), R(cbind(A,B)))  # consistent?
showEqn(A, B)
matlib::Solve(A, B)
limSolve::Solve(A, B)
#
library(limSolve)
G <-matrix(ncol = 2, byrow = TRUE, data = c(1, 0, 0, 1)) 
H <- c(0, 0)
ldei(A, B, G = G, H = H)$X
#
library(cmna) 
gdls(A, B, alpha = 0.05, tol = 1e-06, m = 1e+05) #  least squares with graident descent
jacobi(A, B, tol = 1e-06, maxiter = 100)  # iterativematrix
gaussseidel(A, B, tol = 1e-06, maxiter = 100) # iterativematrix
solvematrix(A, B) # refmatrix

```


Ze dvou kovu o hustotach 7.4 g/cm3 a 8.2 g/cm3 je treba pripravit 0.5 kg slitiny o hustote 7.6 g/cm3. Kolik g kazdiho z kovu je k tomu potreba? 375 g lehciho a 125 g tezsiho

```{r}

library(rootSolve)
model <- function(x){
  F1 <- 7.4*x[1] + 8.2*x[2] - 3800
  F2 <- x[1] + x[2] - 500
  c(F1 = F1, F2 = F2)}
ss <- multiroot(f = model, start = c(1, 1))
ss$root  # [kg]

### graficke reseni
# 7.4*m1 + 8.2*m2 = 0.5*7.6 => m1 = 3800/7.4 - 8.2/7.4 * m2  => m1 = 513.5 - 1.108 * m2 
# m1 + m2 = 500  => m1 = 500 - m2
xx = seq(0,500,1)
yy = 513.5 - 1.108*xx
zz = 500 - xx
plot(xx, yy, type="l",col=2)
points(xx, zz, type="l",col=4)
plot(xx, yy, type="l",col=2,xlim=c(100,150),ylim=c(350,400))
points(xx, zz, type="l",col=4)

### matice
B = c(7.6,1) # [mg]
names(B) = c("kov1","kov2")
r1 = c(7.4,1) # [mg]
r2 = c(8.2,1) # [mg]
A = cbind(r1,r2)
colnames(A) = c("kov1","kov2")
rownames(A) = c("7.4","8.2")
det(A) # matice je regularni n = h
library(matlib)
c(R(A), R(cbind(A,B)))          # show ranks
all.equal(R(A), R(cbind(A,B)))  # consistent?
showEqn(A, B)
rr = matlib::Solve(A, B)
read.table(text = rr[1], fill = TRUE)[[3]]*500 # [g]
read.table(text = rr[2], fill = TRUE)[[3]]*500 # [g]
rr = limSolve::Solve(A, B)
rr*500 # [g]
# 
library(limSolve)
G <-matrix(ncol = 2, byrow = TRUE, data = c(1, 0, 0, 1)) 
H <- c(0, 0)
rr = ldei(A, B, G = G, H = H)$X
rr*500 # [g]

```


V tepelné elektrárne mají zásobu uhlí, která vystací na 24 dní, bude-li v cinnosti pouze první blok, na 30 dní, bude-li v provozu pouze 2. blok a na 20 dní, bude-li v provozu pouze 3. blok. Na jak dloho vystací zásoba, budou-li v provozu vsechny bloky najednou?

```{r}

# x/24 + x/30 + x/20 = 1

library(Ryacas)
vr <- ysym("x/24 + x/30 + x/20 - 1") 
solve(vr, "x") 

fun <- function (x) x/24 + x/30 + x/20 - 1
uniroot(fun, c(0, 1),extendInt = "yes", tol = 1e-9)$root

```