### Vyhledávání

library(CRANsearcher)
#CRANsearcher()

library(pkgsearch)
#pkg_search_addin(query = "", viewer = c("dialog", "browser"))


### Základní matematické funkce

x <- -7
x = -7

# absolutni hodnota (absolute value)
abs(x)

# nasobeni a deleni (multiplication and division)
3*3
x*3
9/3
x/3
(x-3)/(x/3)

1/0
log(0)
Inf*Inf
Inf/Inf
0/0

# celociselne deleni (integer division)
13%/%2
library(numbers)
div(13,2)  
# zbytek po deleni (reminder after division)
13%%2
mod(13,2) 

# Prevod desetinných císel na zlomky (Convert decimal numbers to fractions) ###

library(MASS)
fractions(-0.1666667)
fractions(0.14)
fractions(0.4)

library(schoolmath)
decimal2fraction(0.1,6)

library(fractional)
fractional(-0.1666667, eps = 1e-06, maxConv = 20, sync = FALSE)
numerators(-0.1666667) # vypise citatele (numerator)
denominators(-0.1666667) # vypise jmenovatele (denominator)
numerical(fractional(-0.1666667, eps = 1e-06, maxConv = 20, sync = FALSE))
unfractional(fractional(-0.1666667, eps = 1e-06, maxConv = 20, sync = FALSE))

# zjednodusení zlomku (fractions simplification) ###
library(schoolmath)
cancel.fraction(42, 56)
cancel.fraction(-167, 100)


# mocniny
x = 9
y=4

sqrt(x)
x^(1/2)
x^(-1/2)

x**2
x^2
x^(-2)

library(numbers) # Determine whether p is the power of an integer.
isIntpower(144)
isIntpower(8)


# goniometricke funkce
cos(x)
sin(x)
tan(x)
acos(x)
asin(x)
atan(x)
atan2(y, x)
cospi(x) # cos(pi*x)
sinpi(x) # sin(pi*x)
tanpi(x) # tan(pi*x)

# logaritmy
x = 10
log(x)  # base e
log(x, base = exp(1))  # base e
logb(x, base = exp(1)) # base e
log10(x) # base 10
log(x,base = 10)
log2(x)  # base 2
log(x,base = 2)
log1p(x) # computes log(1+x) 


10^x # exponent
10**x # exponent
exp(x)   # e^x


library(schoolmath)
x <- c(1.002, 2, 3.14, 4, 5, 6.31, 7) 
schoolmath::is.decimal(x)
schoolmath::is.whole(x)

is.integer(x)
library(ttutils)
isInteger(x)
isInteger(c("test", 1, 2, 2.1)) # vektor (vector)

library(litteR)
is_natural_number(12)
is_natural_number(1.2)

library(numbers)
numbers::isNatural(x)


x <- c(1, 2, 3, 4, 5, 6, 7) 
schoolmath::is.even(x)
schoolmath::is.odd(x)
x %% 2 == 0 
x %% 2 != 0
is.even2 <- function(x) x %% 2 == 0; is.even2(x) 
is.odd2 <- function(x) x %% 2 != 0; is.odd2(x) 

x <- c(-1, -2, 3.02, 4, -5.2, 6, -7, 0) 
schoolmath::is.negative(x)
schoolmath::is.positive(x)
schoolmath::is.real.positive(x)
x>0
x<0
x==0


# Prvocisla
x <- c(1, 2, 3, 4, 5, 6, 7) 
schoolmath::is.prim(x)
library(DescTools)
DescTools::IsPrime(x)
library(numbers)
numbers::isPrime(x)

### Nejvetsí spolecný delitel (greatest common divisor) ###
x = 12
y = 8

library(schoolmath)
schoolmath::gcd(x, y)

library(FRACTION)
FRACTION::gcd(x, y)

library(DescTools)
DescTools::GCD(x, y, na.rm = FALSE)
DescTools::GCD(c(10,12,8), na.rm = FALSE)

library(numbers)
numbers::GCD(12, 10) 
numbers::GCD(46368, 75025)
numbers::mGCD(c(2, 3, 5, 7) * 11)
numbers::mGCD(c(2*3, 3*5, 5*7))

numbers::coprime(12, 10) # maji/nemaji prvociselneho delitele
numbers::coprime(46368, 75025) 


### Nejmensí spolecný násobek (smallest common multiple) ###
x = 12
y = 8

library(schoolmath)
schoolmath::scm(x, y)

library(DescTools)
DescTools::LCM(x, y)
DescTools::LCM(c(12,10,8), na.rm = FALSE)

library(numbers)
numbers::LCM(12, 10) 
numbers::LCM(46368, 75025) 
numbers::mLCM(c(2, 3, 5, 7) * 11) 
numbers::mLCM(c(2*3, 3*5, 5*7))


### Vektory


1:6
c(2:12)
rev(c(2:12))
c(12:2)
c(5, -5, -1)
seq(from=2, to=12, by=1)
seq(from=1, to=20, by=2)
seq(5, -5, -1)
seq(from=5, by=-1, along.with = 1:20)
rep(1, len = 9)
rep(1:4, len = 9)
rep(1:4, 2)
rep(1:4, each = 2)       # not the same.
rep(1:4, c(2,2,2,2))     # same as second.
rep(1:4, c(2,1,2,1))
rep(1:4, each = 2, len = 4)    # first 4 only.
rep(1:4, each = 2, len = 10)   # 8 integers plus two recycled 1's.
rep(1:4, each = 2, times = 3)  


sample(12:99, 9, replace = FALSE) # bez opakovani
sort(sample(1:20, 9, replace = FALSE))
sample(12:99, 9, replace = TRUE) # s opakovanim
sort(sample(1:20, 9, replace = TRUE))

set.seed(123)
x = sample(12:99,19, replace = FALSE); sort(x)
is.vector(x)
x
x[1]
x[1:5]
x[c(1:5)]
x[length(x)]
x[(length(x)-1)]
x[-1]
x[-(1:5)]
x[-c(1:5)]
x[-length(x)]
x[length(x)-3:length(x)]#poslednich 6 hodnot
x[c(1,4,7)]#vybrane hodnoty
which.max(x)
which.min(x)
library(Rmpfr)
Rmpfr::which.min(x)
which(x<=50)
x[which(x<=50)]
x[x<=50]
which(x>=50)
x[which(x>=50)]
which(x == max(x))#maximalni hodnota
which(x == 10)
which(x != 10)
x[x > 40]#hodnoty vyssi nez 40.
x[x < 40]#hodnoty mensi nez 40
which(x > 10)#hodnoty vyssi nez 10.
which(x < 10)#hodnoty mensi nez 10.
seq_along(x)
x[x > 80 | x<30] # nebo
x[x < 80 & x>30] # a soucasne

rank(x)  # vhodne pro vektory, nevhodne pro matice 
is.unsorted(x)
sort(x)
sort(x, decreasing = FALSE)
sort(x, decreasing = TRUE)
sort(x[which(x>=50)])
order(x)  # poradi (indexy) pro sort
x[order(x)]

a = c(1, 3, 5, 7) 
b = c(1, 2, 4, 8) 
5 * a
a + b 
a - b 
a * b 
a / b

u = c(10, 20, 30) 
v = c(1, 2, 3, 4, 5, 6, 7, 8, 9) 
u + v 

max(x)
min(x)
range(x)
library(DescTools)
Range(x)
sum(x)
cumsum(x)
prod(x)
cumprod(x)
diff(x)
sign(x)


### Matice

x1 = sample(12:99,19, replace = FALSE)
x1
x2 = sample(12:99,19, replace = TRUE)
x2
x12 = rbind(x1,x2)
x12 = cbind(x1,x2)
is.matrix(x12)
nrow(x12)
ncol(x12)

x12[order(x1),]
x12[order(x2),]
x12[,1]
x12[,"x2"]

colnames(x12)
colnames(x12)=c("X","Y")
rownames(x12)
rownames(x12) = LETTERS[1:length(x1)] 
t(x12)


### Dataframes


x1 = sample(12:99,19, replace = FALSE)
x2 = sample(12:99,19, replace = TRUE)
z = LETTERS[1:length(x1)] 

x12 = rbind(x1,x2,z)
x12 = cbind(x1,x2,z)
nrow(x12)
ncol(x12)
is.data.frame(x12)
is.matrix(x12)

x12 = as.data.frame(x12)
is.data.frame(x12)
colnames(x12)
rownames(x12)
as.matrix(x12)


mean(x1)
mean(as.numeric(x12[,1]))


as.character(x1)

m = cbind(x1,x2,x1,x2)
apply(m,1,mean)


### Seznamy (lists)

x1 = sample(12:99,19, replace = FALSE)
x2 = sample(12:99,19, replace = TRUE)
x3 = sample(12:99,19, replace = TRUE)
z = LETTERS[1:length(x1)] 

x123 = list(x1,x2,x3,z)
is.list(x123)
names(x123) = c("x1","x2","x3","z") 
x123$x1
x123[[1]]


x1 = sample(12:99,19, replace = FALSE)
x2 = sample(12:99,19, replace = TRUE)
x3 = sample(12:99,19, replace = TRUE)
LL = list(x1,x2,x3)
names(LL) = c("x1","x2","x3")

lapply(LL, mean)
sapply(LL, mean)


unlist(LL)
do.call(rbind,LL)
do.call(cbind,LL)

stack(LL)


### Apply funkce
# https://www.r-bloggers.com/r-tutorial-on-the-apply-family-of-functions/
# http://www.datasciencemadesimple.com/apply-function-r/

### lapply a sapply (vektory, seznamy) ###

### Vypocitejte molekulove hmotnosti alkanu C1 - C12.
# lapply
Mhc = lapply(c(1:12), function(n){n*12.011 + (2*n + 2)*1.008})
Mhc = unlist(Mhc)
names(Mhc) = c(1:12)
Mhc
# sapply
Mhc = sapply(c(1:12), function(n){n*12.011 + (2*n + 2)*1.008})
names(Mhc) = c(1:12)
Mhc

library(live)
data(wine)
wine
head(wine)

# mean pomoci sapply a lapply
wn = sapply(c(1:length(wine[1,])), function(x){mean(wine[,x])})
names(wn) = colnames(wine)
wn
#
wn = lapply(c(1:length(wine[1,])), function(x){mean(wine[,x])})
names(wn) = colnames(wine)
wn

### tapply (matice, dataframes) ###

tapply(wine$alcohol,wine$quality,mean)
tapply(wine$alcohol,wine$quality,median)
tapply(wine$alcohol,wine$quality,length)
tapply(wine$alcohol,wine$quality,boxplot.stats)

### apply (matice, dataframes) ###

apply(wine,1,length)
apply(wine,2,length)
apply(wine,2,mean)
apply(wine,2,median)


### Matematické a fyzikální konstanty


### Ludolfovo cislo ####
pi
print(pi, digits=20)

### Eulerovo cislo ####
exp(1)
print(exp(1), digits=20) 

library(marelac)
data(Constants)
Constants

Constants$gasCt1
Constants$gasCt2
Constants$gasCt3
R = as.numeric(Constants$gasCt2[1]); R # hodnota
Constants$gasCt2[2] # jednotka


library(constants)
codata
lookup("planck", ignore.case=TRUE)
lookup("avogadro", ignore.case=TRUE)
lookup("faraday", ignore.case=TRUE)
lookup("molar", ignore.case=TRUE)
lookup("molar mass", ignore.case=TRUE)
lookup("light", ignore.case=TRUE)
lookup("boltzmann", ignore.case=TRUE)
lookup("mass", ignore.case=TRUE)
lookup("mass constant", ignore.case=TRUE)
lookup("proton mass", ignore.case=TRUE)
lookup("electron mass", ignore.case=TRUE)

# Univerzalni plynova konstanta
lookup("gas", ignore.case=TRUE)
codata[134,]
R = as.numeric(codata[134,5]) # hodnota
Ru = as.numeric(codata[134,6]) # nejistota
codata[134,7] # jednotka



### Prevody jednotek

library(scales) # %
percent(0.05)

library(pracma) # uhly vs radiany
deg2rad(120) 
rad2deg(3.14)


library(measurements)
conv_unit_options # pouzivane jednotky

# length
conv_unit(1, from="mm", to="m")
conv_unit(1, from="nm", to="angstrom")

# pressure
conv_unit(101, from="Pa", to="mmHg")
conv_unit(1, from="atm", to="Pa")

# temperature
conv_unit(100, from="C", to="K")

# energy
conv_unit(500, from="cal", to="J")
conv_unit(1, from="erg", to="J")

#volume
conv_unit(100, from="l", to="m3")

# mass
conv_unit(100, from="g", to="mg")
conv_unit(100, from="ug", to="g")

# speed
conv_unit(1, from="km_per_hr", to="m_per_sec")
conv_unit(1, from="light", to="km_per_hr")
conv_unit(1, from="light", to="m_per_sec")


### multiunits
conv_multiunit(x = 1, from="ug / l", to="g / m3")
conv_multiunit(x = 1, from="mg / l", to="kg / m3")
conv_multiunit(x = 1, from="cal / kg", to="J / g")


#### nasobky jednotek SI
library(sitools)
f2si(10000)
f2si(0.023, unit="l")
f2si(3.5e-5, unit="mol")
numbers <- c(1e5, 3.5e-12, 0.004) 
f2si(numbers, unit="g")

sapply(numbers,function(x){f2si(x, unit="g")})


# Prevod casovych jednotek

library(datetime)
xx = as.minute(60) # 60 min
as.second(xx)
as.hour(xx)
as.day(xx)

dur = duration(week=0,day=0.04166667,hours=1,minutes=0,seconds=3600)
as.numeric(dur, "minutes")


### Plynný systém mení svuj objem o 1200 ml za konstantního vnejsího tlaku 30 atm.Jakou práci vykoná plyn pri expanzi? Vyjádrete v ruzných energetických jednotkách.
# W = p . dV
dV = 1200 # [ml]
dV = conv_unit(dV, from="ml", to="m3"); dV
p = 30 # [atm] to [Pa]
p = conv_unit(p, from="atm", to="Pa"); p
W = p*dV; W # [J]
W1 = conv_unit(W, from="J", to="cal"); W1 # [cal]
W2 = conv_unit(W, from="J", to="Wsec"); W2 # [W.s]
W3 = conv_unit(W, from="J", to="erg"); W3 # [erg]


### Kolik tepla se uvolní pri pruchodu 26.43 coulombu elektrického proudu vodicem pri potenciálovém spádu 2.432 V.
# 64.33 J, 64.33e7 erg, 15.37 cal
# Q = U . I . t = U . q 
q = 26.43 # [C]
U = 2.432 # [V]
Q = U*q; Q # [J]
Q1 = conv_unit(Q, from="J", to="cal"); Q1 # [cal]
Q3 = conv_unit(Q, from="J", to="erg"); Q3 # [erg]


### Pokus trval 7658 minut. Kolik je to dní, hodin a minut?
tt = 7658 # [min]

library(lubridate)
seconds_to_period(tt*60)
library(datetime)
seconds_to_period(as.second(as.minute(tt)))


##### Prvky

library(PeriodicTable)
data(periodicTable)
periodicTable

library(CHNOSZ)
# atomova hmotnost
AW = mass("Ca"); AW
# molarni hmotnost
MW = mass("CaCO3"); MW
# pocet atomu prvku ve vzorci
ce = count.elements("CaCO3"); ce
as.chemical.formula(ce, drop.zero = TRUE)

library(IsoSpecR)
data(isotopicData)
isotopicData

library(loon.data)
data(elements)
elements


#### Zaokrouhlování na daný pocet desetinných míst

round(1234.567,2)
round(123.456,digits=2)
ceiling(1234.567)
floor(1234.567)

library(guf)
round_something(123.456, decimals = 2)
round_something("123.456", decimals = 2)

library(PKNCA)
roundString(3141.59265, digits = 3, sci_range = 0, sci_sep = "e")
roundString(3141.59265, digits = 3, sci_range = 0, sci_sep = "x10^")
roundString(3141.59265, digits = 3, sci_range = Inf, sci_sep = "e")


### Zaokrouhlování na daný pocet platných císlic

signif(1234.567,4)
signif(123.456,digits=4)

library(PKNCA)
signifString(3141.59265, digits = 1, sci_range = 0, sci_sep = "e")
signifString(3141.59265, digits = 1, sci_range = 0,sci_sep = "x10^")
signifString(3141.59265, digits = 1, sci_range = Inf, sci_sep = "e")


### NA and data imputation 1

### vektory

library(datasets)
head(airquality)
summary(airquality)
airquality[,"Ozone"]
airquality[,"Solar.R"]

x = airquality[,"Ozone"]
length(x)
any(is.na(x))
is.na(x[5])
is.na(x[6])

xx[6]<- NA
xx[6]

library(jwutil) # rel. podil hodnot ve vektoru
propIsNa(x)

xx = na.omit(x)
length(xx)

xx[6]<- NA
xx[6]


### Matice a dataframes

x = airquality
any(is.na(x))
is.na(x)

library(splus2R)
splus2R::which.na(x)

library(GLDEX)
GLDEX::which.na(x)

apply(x,2,splus2R::which.na)


library(ryouready)
count_na(x, along = 2)
# along: 1 = rows, 2=columns

xl = apply(x,2,splus2R::which.na)
sapply(xl,length)

library(jwutil) # rel. podil hodnot ve sloupcich
propNaPerField(x)

library(na.tools)
all_na(x)
any_na(x)
is_na(x)
which_na(x[,1])
which_na(x[1,])


library(agricolae)
delete.na(x, alternative="less")
delete.na(x, alternative="greater")
# "less"    = vynechava sloupce
# "greater" = vynechava radky

xx = delete.na(x, alternative="greater")
xx
xx[5,5]<- NA


library(jwutil) # NA na 0
zero_na(x)
zero_na(x, na_ish = TRUE)

library(DescTools)
ZeroIfNA(x)
NAIfZero(x)
#
BlankIfNA(x, blank="")
BlankIfNA(x, blank="UDL")
NAIfBlank(x)

# nahrazeni medianem sloupce
library(DescTools) 
x1 = Impute(x[,1], FUN = function(y) mean(y, na.rm = TRUE))
x2 = Impute(x[,1], FUN = function(y) median(y, na.rm = TRUE))
cbind(x[,1],x1,x2)

x[,1] = x1 # nahrazeni

library(na.tools)
na.constant(x[,1], 5)
na.max(x[,1])
na.min(x[,1])
na.mean(x[,1])
na.median(x[,1])
na.quantile(x[,1], prob=0.4)
na.mode(x[,1])
na.most_freq(x[,1])
na.inf(x[,1])
na.neginf(x[,1])
na.true(x[,1])
na.false(x[,1])
na.zero(x[,1])

set.seed(123)
xb = na.bootstrap(x[,5])
xr = na.resample(x[,5])
rbind(x[,5],xb,xr)


### opakujici se hodnoty (ties)

duplicated(c("a","b","c", "a"), incomparables = FALSE)

library(guf)
count(c("a","b","c", "a"), group.by = NULL, split.by = NULL, row.total = FALSE)


### Nejistoty a chyby mereni

### Zapis chyby (nejistoty)
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)

# elementární náboj
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 stanoven jejich bod tání. Vypocítejte hodnotu bodu tání v zásilce a její smerodatnou odchylku.
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)
options(errors.notation = "parenthesis"); print(e, digits = 2)

### Teplota vody byla merena dvema teplomery s ruznou presností
t1 = c(20.8, 21.1, 20.9)
t2 = c(20.95, 20.92, 20.94)
e1 <- set_errors(mean(t1), sd(t1)) 
options(errors.notation = "plus-minus"); print(e1, digits = 1)
options(errors.notation = "parenthesis"); print(e1, digits = 1)
round(100*sd(t1)/mean(t1),3) # relativni chyba [%]
e2 <- set_errors(mean(t2), sd(t2)) 
options(errors.notation = "plus-minus"); print(e2, digits = 1)
options(errors.notation = "parenthesis"); print(e2, digits = 1)
round(100*sd(t2)/mean(t2),3) # relativni chyba [%]


#### Sireni chyb ###################################

##### mean + S.E. 
## Error propagation using Taylor expansion of first- (Mean.1, sd.1) and second-order (Mean.2, sd.2)
## Monte Carlo (simulace)

## Priklad (bez uvedení poctu stupnu volnosti, resp. poctu opakovani)
library(propagate)
EXPR1 <- expression(x/y) 
x <- c(5, 0.01) 
y <- c(1, 0.01) 
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)


### Vypocítejte hustotu a její chybu pro látku, u níz byla opakovaným merením stanovena hmotnost 6.824 (0.008) g a objem 3.03 (0.01) ml.
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)


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


### Na píst o prumeru (200 +- 0.05) mm pusobí pára tlakem (8.2 +- 0.1) atm. Jakou silou pusobí pára na píst? Urcete také relativní chybu výsledku.
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)


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


### Vypocítejte koeficient viskozity roztoku glycerinu Stokesovou metodou
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)


library(measurements) # prevod jednotek
### Vypocítejte koeficient viskozity roztoku glycerinu pomoci kapilarniho viskozimetru
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)


## NEFUNGUJE pro jednu promennou, viz "metRology"
### Soucin rozpustnosti stríbrné soli AgX má hodnotu Ks = (4.0 +- 0.4) x 10-8. Jaká je chyba vypoctené rovnovázné koncentrace stríbrných iontu ve vode?
EXPR <- expression(sqrt(Ks))
Ks = c(4.0, 0.4)
DAT <- data.frame(Ks)
res <- propagate(EXPR,DAT)


##### Metoda 2: Výpocet pomocí intervalu 

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


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


# Soucin rozpustnosti stríbrné soli AgX má hodnotu Ks = (4.0 +- 0.4) x 10-8. Jaká je chyba vypoctené rovnovázné koncentrace stríbrných iontu ve vode?
GUM("Ks",4.0e-8,0.4e-8,1,"sqrt(Ks)",sig.digits.U = 2) # tohle funguje
GUM.validate("Ks",4.0e-8,0.4e-8,1,"A","Normal","sqrt(Ks)")
GUM.validate("Ks",4.0e-8,0.4e-8,1,"B","Normal","sqrt(Ks)")

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

# Priklad (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


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


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
z <- x / y; z
#(z <- x / y)

# s korelací
correl(x, y) =  c(0.8864282, 0.4761841, -0.4140209, -0.4496266, 0.2911837) # zavedení korekace chyb
covar(x, y)
(z_correl <- x / y)


################################################################################################################


### vkladani hodnot z konzole

widmark <- function(){ 
  print("Výpocet hladiny alkoholu v krvi podle Widmarka")
  print("Data zadavejte pres konzoli !!!")
  sex<-readline(prompt="Jste muz (m) nebo zena(z)? " )
  if (sex == "m") {BW = 0.58; MR = 0.015} 
  if (sex == "z") {BW = 0.49; MR = 0.017}
  Wt <-readline(prompt="Jaka je Vase hmotnost (v kg)? " )
  Wt = as.numeric(Wt)
  # print("Jeden standardni alkoholicky napoj odpovida 10 g resp. 12.7 ml cisteho ethanolu")
  # SD <-readline(prompt="Kolik jste vypil standardnich alkoholickych napoju? " )
  # prepocitat na piva, vina a panaky
  p0 <-readline(prompt="Kolik sklenic nealkoholickeho piva (0.5 l) jste vypil ? " )
  p10 <-readline(prompt="Kolik sklenic 10 stupnoveho piva (0.5 l) jste vypil ? " )
  p12 <-readline(prompt="Kolik sklenic 12 stupnoveho piva (0.5 l) jste vypil ? " )
  vb <-readline(prompt="Kolik sklenicek bileho vina (0.2 l) jste vypil ? " )
  vc <-readline(prompt="Kolik sklenicek cerveneho vina (0.2 l) jste vypil ? " )
  p30 <-readline(prompt="Kolik panaku 30% alkoholu (0.05 l) jste vypil ? " )
  p45 <-readline(prompt="Kolik panaku 45% alkoholu (0.05 l) jste vypil ? " )
  # SD = as.numeric(SD)
  p0 = as.numeric(p0) * 2 # [g 100% ethanolu]
  p10 = as.numeric(p10) * 16 # [g 100% ethanolu]
  p12 = as.numeric(p12) * 20 # [g 100% ethanolu]
  vb = as.numeric(vb) * 16 # [g 100% ethanolu]
  vc = as.numeric(vc) * 19 # [g 100% ethanolu]
  p30 = as.numeric(p30) * 12 # [g 100% ethanolu]
  p45 = as.numeric(p45) * 18 # [g 100% ethanolu]
  SD = sum(c(p0,p10,p12,vb,vc,p30,p45))/10
  DP <-readline(prompt="Pred kolika hodinami? " )
  DP = as.numeric(DP)
  EBAC = ((0.806*SD*1.2/(BW*Wt)) - MR*DP)*10   # [promile, resp. g/l]
  # if (EBAC <= 0) {EBAC = 0}
  cat("Mate v krvi", round(EBAC,2) , "promile (g/l) alkoholu.")
  # print(paste("Mate v krvi", round(EBAC,2) , "promile (g/l) alkoholu."))
}
widmark()

# https://en.wikipedia.org/wiki/Blood_alcohol_content
# https://en.wikipedia.org/wiki/Standard_drink
# https://www.bezpecnecesty.cz/cz/alkohol-kalkulacka/algoritmus-vypoctu
# https://autobild.pluska.sk/poradca/chcete-vypocitat-hladinu-alkoholu-poradime-vam
#

################################################################################################################

# f <- function(<arguments>) {## Do something interesting}

my_function <- function(country = "Czech Republic") {
  paste("I am from", country)
}

Nested_function <- function(x, y) {
  a <- x + y
  return(a)
}
Nested_function(2,2)
Nested_function(3,3))

Outer_func <- function(x) {
  Inner_func <- function(y) {
    a <- x + y
    return(a)
  }
  return (Inner_func)
}
output <- Outer_func(3) # To call the Outer_func
output(5)

############################################################################

# pipe operator in R example (alternative - intermediate variables, r pipe)
a = 3.14159
b = seq(from = a,10,3)
round(b,3)

# pipe operator in R example (alternative - nested syntax) 
a = 3.14159
seq(round(a,3),10,3)

# pipe operator in R - example (assignment pipe r)
a = 3.14159
a %>% seq(10,3) %>% round(3)



###############################################################################################################################################

### if - else ###

### Vypocítejte vazebnou energii atomového jádra pomoci Bethe - Weiszackerovy rovnice
#     B = 14.0 - 13.1*A^(2/3) + 0.585*Z*(Z-1)/A^(1/3) - (18.1*(A-2*Z)^2)/A + C/A
# A - nukleonove cislo, Z - atomove cislo, C - konstanta, pro jadra se sudym poctem protonu a neutronu je rovna 132, 
# pro jadra s lichym poctem protonu a neutronu je rovna -132, pro lichy pocet protonu a sudy pocet neutronu a naopak je rovna 0.

# alpha particle
A = 4 
Z = 2

# U-235
A = 235 
Z = 92

# funkce
N = A-Z
if (N%%2==0 & Z%%2==0){C = 132} else if (N%%2!=0 & Z%%2!=0){C = -132} else {C = 0}
B = 14.0 - 13.1*A^(2/3) + 0.585*Z*(Z-1)/A^(1/3) - (18.1*(A-2*Z)^2)/A + C/A
B # [MeV]
B/A # energie na 1 nukleon [MeV]

# polomer jadra: R = R0 * A^(1/3)
R0 = 1.4e-13 # [cm]
R = R0 * A^(1/3)
R = R*10e-15 
R   # [m]


### Anorganická kvalitativní analýza - 1. skupina

KAA1 <- function(){ 
  print("Kvalitativni analyza - 1. skupina iontu (Ag, Pb, Hg)")
  print("Odpovidejte a (ano) nebo n (ne)")
  print("Pozor. Nez zacnete, prepnete se do konzole.")
  an <- readline(prompt="Chcete pokracovat? " )
  if (an == "n") {print("Tak snad nekdy jindy ...")
  } else {print("Takze zaciname ...")  
    an <- readline(prompt="Vznikla po pridani HCl bila srazenina? " )
    if (an == "n") {print("Neni pritomen zadny prvek 1. skupiny.")
    } else {print("Je pritomno olovo, stribro, rtut, nebo jejich kombinace.")
      an <- readline(prompt="Rozpousti se srazenina zcela v horke vode? " )
      if (an == "a") {print("Je pritomno pouze olovo.")
      } else {
        an <- readline(prompt="Rozpousti se srazenina alespon castecne? " )  
        if (an == "a") {print("Je pritomno olovo, ale take stribro a/nebo rtut.")
        } else {print("Olovo neni pritomno, jsou pritomny stribro a/nebo rtut.")}
        an <- readline(prompt="Rozpousti se srazenina zcela v amoniaku? " )
        if (an == "a") {print("Je pritomno stribro, neni pritomna rtut.")
        } else {
          an <- readline(prompt="Rozpousti se srazenina alespon castecne? " )  
          if (an == "a") {print("Je pritomno stribro, a take rtut.")
          } else {print("Neni pritomno stribro, je pritomna rtut.")}
          an <- readline(prompt="Doslo pridavkem amoniaku ke zcernani srazeniny? " )
          if (an == "a") {print("Rtut je potvrzena.")
          } else {print("To je divne. Nekde se stala chyba ...")}
        } 
      }   
    }
    print("Hotovo. Nashledanou priste.")
  }
}
KAA1()



### while cyklus (while loop) ###

### Vypocitejte molekulove hmotnosti alkanu C1 - C12
n = 1
Mhc = NULL
while(n <= 12){
  hc = n*12.011 + (2*n + 2)*1.008
  Mhc = c(Mhc,hc)
  n = n+1
}
names(Mhc) = c(1:12)
Mhc


### switch
switch(2,"red","green","blue")
switch(1,"red","green","blue")

1

### for cyklus (for loop) ###


### Vypocitejte molekulove hmotnosti alkanu C1 - C12

Mhc = NULL
for(n in c(1:12)){
  hc = n*12.011 + (2*n + 2)*1.008
  Mhc = c(Mhc,hc) 
}
names(Mhc) = c(1:12)
Mhc


### Vypocitejte stredni kvadratickou rychlost molekul danych plynu v rozmezi teplot od 150 K do 500 K s krokem 50 K.

el = c("H2","He", "O2", "Kr")
M = c(2.0159, 4.0026, 31.9988, 83.8) # [g / mol]
Ts = seq(from=150, to=500, by=50) # [K]
R = 8314.34 # [J / kmol K]
MKV = NULL
for(tt in Ts){
  mkv= rep(0, length(M))
  for(ii in c(1:length(M))){
    mkv[ii] = (3*R*tt/M[ii])^(1/2)
  }
  MKV = rbind(MKV, mkv)   
}
colnames(MKV) = el
rownames(MKV) = Ts
MKV
plot(0,0,xlim=c(min(Ts),max(Ts)),ylim = c(min(MKV),max(MKV)), type ="n",xlab="T [K]",ylab="str. kv. rychlost [m/s]")
for(ii in c(1:length(M))){points(Ts,MKV[,ii], type="p",col=ii,pch=16)}
legend("topleft",legend=el,col=c(1:length(M)),pch=16,cex=.8)


#### mnoziny ####################################################################################

(x <- c(sort(sample(1:20, 3))))
(y <- c(sort(sample(3:23, 7))))

# x je prvkem mnoziny y
is.element(x, y) 
x %in% y 
library(prob) # pro cely vektor
isin(x,y, ordered = FALSE)
all(x %in% y)

library(prob)
x <- c(3,3,2,2,3,3,4,4)
isrep(x)  # one pair each of 2s and 4s
countrep(x)
isrep(x, nrep = 4)
countrep(x, nrep = 4)
isrep(x, vals = 4) # one pair of 4s
countrep(x, vals = 4) 

union(x, y) # sjednoceni
intersect(x, y) # prunik
setdiff(x, y) # rozdil mnozin
setdiff(y, x) # rozdil mnozin
setequal(x, y) # shoda mnozin
# https://stat.ethz.ch/R-manual/R-devel/library/base/html/sets.html
subset(x, .)


library(prob)
union(x, y)
setdiff(x, y)
intersect(x, y)
subset(x, ...)

table()


### Vennovy diagramy
library(gplots)
venn( list(A=1:5,B=4:6,C=c(4,8:10)) )
v.table<-venn( list(A=1:5,B=4:6,C=c(4,8:10),D=c(4:12)) ) 
print(v.table)




### numericke metody ###########################################################################################################################

library(animation) # animation
bisection.method(FUN = function(x) x^2 - 4, rg = c(-1, 10), tol = 0.001, interact = FALSE)

library(animation) # animation
newton.method(FUN = function(x) x^2 - 4, init = 10, rg = c(-1, 10), tol = 0.001, interact = FALSE, col.lp = c("blue", "red", "red"))


library(cmna) # numerical
bisection(f = function(x) x^2 - 4, a=-1, b=10, tol = 0.001, m = 100)
newton(f = function(x)x^2 - 4, fp=fp, x=10, tol = 0.001, m = 100) # nefunguje



### polynomy ####

library(polynom)
(p <- polynomial(c(1,0,0,3,4)))
str(p)
p + polynomial(1:2)
p*p
f1 <- as.function(p)
f1(pi)
f1(matrix(1:6,2,3))


library(multipol)
(a <- as.multipol(matrix(1:10,nrow=2)))
b <-  as.multipol(matrix(1:10,ncol=2))
a+b
a*b


### kvadraticke a kubicke rovnice ###

# Velikost vyslednice dvou navzajem kolmych sil je 34 N. Jaka je velikost skladanych sil, je-li jedna z nich o 14 N vetsi nez druha?  
# 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) {print("Kvadraticka rovnice ma dva sobe rovne realne koreny (dvojnasobny koren).")}
if (D > 0) {print("Kvadraticka rovnice ma dva ruzne realne koreny.")}
if (D < 0) {print("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


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)
F1 <- uniroot(fun, c(0, 100),extendInt = "yes", tol = 1e-9)$root 
F1


fun <- function (x) x^2 + 14*x - 480
# uniroot(fun, c(-1, 1),extendInt = "yes", tol = 1e-9)
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)$root 
F1
F2 = F1 + 14
F2

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

library(pracma)
rr = roots(c(1, 14, -480))
rr



## graficke reseni
xx = seq(-50,50,1)
yy = xx^2
zz = -14*xx + 480
plot(xx,yy,type="l")
abline(480,-14,type="l",col=2)
plot(xx,yy,type="l",xlim = c(10,20),ylim=c(0,500))
abline(480,-14,type="l",col=2)

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


### Jake pH ma roztok kyseliny mravenci o koncentraci 8.5 x 10-4 mol/l?
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 



### Tlakova lahev s oxidem uhlicitym obsahuje 10.0 kg plynu. Jaký objem zaujíma stlaceny plyn, kdyz pri teplote 30 °C je tlak v lahvi 13.17e6 Pa? Vypocet provedte pomoci van der Waalsovy rovnice
# (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 neidealni 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)


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



### soustavy linearnich rovnic ###

### Ze dvou slitin s 60% a 80% obsahem medi se ma ziskat 40 kg slitiny se 75% obsahem medi. Kolik kg kayde slitiny je treba pouzit?  [10 a 30 kg]

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

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)


### Ze dvou druhu caje v cene 160 Kc a 220 Kc za 1 kg je treba pripravit 20 kg smesi v cene 205 Kc za 1 kg. Kolik kg kazdeho caje je treba smichat?   5 kg levnejsiho a 15 kg drazsiho 

library(rootSolve)   ### opravit zadani
model <- function(x){
  F1 <- 160*x[1] + 220*x[2] - 4100
  F2 <- x[1] + x[2] - 20
  c(F1 = F1, F2 = F2)}
ss <- multiroot(f = model, start = c(1, 1))
ss$root

### graficke reseni
# 160*m1 + 220*m2 = 4100 => m1 = 4100/160 - 220/160 * m2  => m1 = 25.625 - 1.375 * m2
# m1 + m2 = 20 => m1 = 20 - m2
xx = seq(0,20,0.1)
yy = 25.625 - 1.375*xx
zz = 20 - xx
plot(xx, yy, type="l",col=2)
points(xx, zz, type="l",col=4)
plot(xx, yy, type="l",col=2,xlim=c(10,20),ylim=c(0,10))
points(xx, zz, type="l",col=4)


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

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

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]



### Vodní nádrz by se naplnila prvním prívodem za 36 minut, druhým za 45 minut. Za jak dlouho se nádrz naplní, pritéká-li voda nejprve 9 minut prvním prívodem a pak obema soucasne?
# x/36 + (x - 9)/45 = 1
vr <- ysym("x/36 + (x-9)/45 - 1") 
solve(vr, "x") 

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


### 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?
# x/24 + x/30 + x/20 = 1
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


### Míc byl vyhozen svisle vzhuru rychlostí 25 m/s. Za jak dlouho bude míc 20 m nad zemí? Odpor vzduchu zanedbejte.
# h = v t - 1/2 g t^2
g = 10 # [m/s^2]
h = 20 # [m]
v = 25 # [m/s]
t <- ysym("t")
g <- ysym("g")
h <- ysym("h")
v <- ysym("v")
poly <- ysym("g*t^2 - 2*v*t + 2*h") 
ko = solve(poly, "t")   
ko[1]
ko[2]
eval(as_r((2*v+sqrt(-((-4)*v^2+8*g*h)))/(2*g)), list(g = 10, v = 25, h = 20))
eval(as_r((2*v-sqrt(-((-4)*v^2+8*g*h)))/(2*g)), list(g = 10, v = 25, h = 20))

poly <- ysym("10*t^2 - 50*t + 40") 
ko = solve(poly, "t")   
ko[1]
ko[2]

library(pracma)
rr = roots(c(10, -50, 40))
rr
F1 = rr[rr>=0]
F1

rr = polyroot(c(40, -50, 10)) # base
Re(rr)

g = 10 # [m/s^2]
h = 20 # [m]
v = 25 # [m/s]
fun <- function (t) g*t^2 - 2*v*t + 2*h
#uniroot(fun, c(-1, 0),extendInt = "yes", tol = 1e-9)  # base
uniroot(fun, c(0.1, 5),extendInt = "yes", tol = 1e-9)
F1 <- uniroot(fun, c(0, 1),extendInt = "yes", tol = 1e-9)$root 
F1 # udava jen jeden koren

library(rootSolve)
F1 <- rootSolve::uniroot.all(fun, c(0, 10), tol = 1e-9)
F1 

library(Rmpfr)
F1 <- Rmpfr::unirootR(fun, lower=0.1, upper=10, tol = 1e-9)
F1 # udava jen jeden koren

# Descartes rule
cat("Pocet kladnych korenu:", sum(diff(sign(c(10, -50, 40))) != 0)) # pocet kladnych korenu


###### derivace ####################################################################################################

D() 
deriv() 

derivative = deriv(~ x^3, "x") 
x <- 2
eval(derivative )

f = expression(x^3)
dx2x <- D(f,"x")
dx2x 


library(pracma)
f <- sin
xs <- seq(-pi, pi, length.out = 100)
ys <- f(xs)
y1 <- fderiv(f, xs, n = 1, method = "backward")
y2 <- fderiv(f, xs, n = 2, method = "backward")
y3 <- fderiv(f, xs, n = 3, method = "backward")
y4 <- fderiv(f, xs, n = 4, method = "backward")
plot(xs, ys, type = "l", col = "gray", lwd = 2,xlab = "", ylab = "", main = "Sinus and its Derivatives")
lines(xs, y1, col=1, lty=2)
lines(xs, y2, col=2, lty=3)
lines(xs, y3, col=3, lty=4)
lines(xs, y4, col=4, lty=5)
grid()

f <- sin
xs <- seq(-pi, pi, length.out = 100)
numdiff(f, xs, maxiter = 16, h = 1/2)


# derivace v bode
library(pracma)
f <- function(x) sin(x)*sqrt(1+sin(x))
F <- function(x)integrate(f, 0, x, rel.tol = 1e-12)$value
x0 <- 1
(dF0 <- numderiv(F, x0, tol = 6.5e-15))      
f(x0)
x=x0
eval(sin(x)*sqrt(1+sin(x)))

fderiv(F, x0)

library(numDeriv)
numDeriv::grad(F, x0)



### Hmotny bod kona harmonicky pohyb, zavislost vychylky y na case t je dana vztahem y = A * sin(om*t + phi), kde A, om a phi jsou konstanty. Urcete rychlost hmotneho bodu v case t = 0.
library(Deriv)
der = Deriv("A*sin(om*t+phi)", "t")
der
der0 = gsub("t", "0", der, fixed = TRUE)
der0
Simplify(der0)




library(mosaic)

f <- makeFun(m * x + b ~ x, m = 3.5, b = 10)
f(x = 2)

g <- makeFun(A * x * cos(pi * x * y) ~ x + y, A = 3)
# function (x, y, A = 3) A * x * cos(pi * x * y)
g
g(x = 1, y = 2)

plotFun(A * exp(k * t) * sin(2 * pi * t/P) ~ t + k, t.lim = range(0, 10), k.lim = range(-0.3,0), A = 10, P = 4)

library(manipulate)
plotFun(A * exp(k * t) * sin(2 * pi * t/P) ~ t + k, t.lim = range(0, 10),k.lim = range(-0.3,0), A = 10, P = 4, surface = TRUE)



library(mosaic)
library(mosaicCalc)
D(sin(x) ~ x)   # derivace
D(A * x^2 * sin(y) ~ x)  # 2. derivace
D(A * x^2 * sin(y) ~ x + y)  # 3. derivace


findZeros(sin(t) ~ t, t.lim = range(-5, 1))
findZeros(sin(t) ~ t, nearest = 5, near = 10)
findZeros(sin(t) ~ t, near = 0, within = 8)


D = function(f,delta=.000001){function(x){ (f(x+delta) - f(x-delta))/(2*delta)} }
f = function(x){ x^2 + 2*x }
plot(f, 0, 10)



###### integrace ####################################################################################################

# integrate() 

# 1/((x+1)*sqrt(x))
integrand <- function(x) {1/((x+1)*sqrt(x))}
integrate(integrand, lower = 0, upper = Inf)

# 1/sqrt(2*pi)*exp(-x^2/2)
f <- function(x) {1/sqrt(2*pi)*exp(-x^2/2)}
integrate(f, lower = -1.96, upper = 1.96)


### mnohorozmerny integral

library(cubature)
f <- function(x) { 2/3 * (x[1] + x[2] + x[3]) }  
adaptIntegrate(f, lowerLimit = c(0, 0, 0), upperLimit = c(0.5, 0.5, 0.5))



library(sfsmisc)
integrate.xy(x, fx=y, a, b, use.spline=TRUE, xtol=2e-08)


### geometrie ##########################################################################################################################################

library(retistruct)
line.line.intersection(c(0, 0), c(1, 1), c(0, 1), c(1, 0))
line.line.intersection(c(0, 0), c(0, 1), c(1, 0), c(1, 1))


library(PBSmapping)
p1 <- data.frame(PID=rep(1, 4), POS=1:4, X=c(1,1,6,6), Y=c(1,3,3,1))
p2 <- data.frame(PID=rep(2, 5), POS=1:5, X=c(4,4,8,8,6), Y=c(2,4,4,2,1))
p3 <- joinPolys(p1,p2)
x11()
par(mar=c(3,3,1,1))
plot(1,1,ylim=c(0,5),xlim=c(0,9), t="n", xlab="", ylab="")
polygon(p1$X, p1$Y, border=2)
polygon(p2$X, p2$Y)
polygon(p3$X, p3$Y, col=rgb(0,0,1,0.2))


### kuzelosecky 

library(sfsmisc)
ellipsePoints(a, b, alpha = 0, loc = c(0, 0), n = 201)

library(conics)
# 2x2^2 + 2x1x2 + 2x2^2 - 20x1 - 28x2 + 10 = 0 
v <- c(2,2,2,-20,-28,10)
A <- conicMatrix(v)
conicCenter(v)
conicAxes(v)
conicCenter(A)
conicAxes(A)
conicPlot(v, main="conicPlot(v)", xlab="", ylab="")
conicPlot(v, col="red", lty=3)
conicPlot(v, asp=1)
conicPlot(v, center=T, sym.axes=T)

v <- c(2,2,2,-20,-28,10)
conicPlot(v, center=T, lty=1)
v[6] <- 30 
conicPlot(v, add=T, lty=2)
v[6] <- 50 
conicPlot(v, add=T, lty=4)

v <- c(2,2,-2,-20,20,10) 
conicPlot(v, xlim=c(-10,10), ylim=c(-5,18))
conicPlot(v, asymptotes=T, sym.axes=T,as.col="red", as.lty=2, ax.col="blue", ax.lty=4, xlim=c(-10,10), ylim=c(-5,18))
conicPlot(v, asymptotes=T, sym.axes=T, xlim=c(-10,10), ylim=c(-5,18), asp=1, col="blue", main="Hyperbola", bty="n")

v <- c(-4,0,1,0,0,1) 
cp <- conicPlot(v, sym.axes=TRUE, asymptote=TRUE, asp=1, ax.lty=2, as.col="gray") 
points(cp$foci$x,cp$foci$y,col="red",pch=19) 
text(cp$foci$x,cp$foci$y+0.1,paste("F",2:1,sep="")) 
points(cp$vertices$x,cp$vertices$y,col="blue",pch=19)



#### SYMBOLICKA MATEMATIKA ##################################################################


library(rSymPy)

# výpocet na daný pocet platných císlic
sympy("sqrt(8).evalf()")  # default
sympy("sqrt(8).evalf(50)")
sympy("pi.evalf()")
sympy("pi.evalf(3)")
sympy("pi.evalf(120)")


## Recognizing numbers
sympy("nsimplify(4.242640687119286)")
sympy("nsimplify(cos(pi/6))")
sympy("nsimplify(6.28, [pi], tolerance=0.01)")
sympy("nsimplify(pi, tolerance=1e-5)")
sympy("nsimplify(pi, tolerance=1e-6)")
sympy("nsimplify(29.60881, constants=[pi,E], tolerance=1e-5)")


# 
x <- Var("x")
x+x
x*x/x
y <- Var("x**3")
x/y

z <- sympy("2.5*x**2")
z + y

sympy("one = cos(1)**2 + sin(1)**2")
sympy("(one - 1).evalf()")  # chyba zaokrouhlení
sympy("(one - 1).evalf(chop=True)")  # rouding this type of roundoff errors


# FUNKCE
sympy("sqrt(8)")
sympy("sqrt(9)")  
sympy("sqrt(3.14)")  
sympy("sin(1)")  
# create function 
Exp <- function(x) Sym("exp(", x, ")") 
Exp(-x) * Exp(x) 


# rovnice
sympy("Eq(x**2+2*x+1,(x+1)**2)") # vytvoreni rovnice
# jsou obe strany rovnice ekvivalentni?
sympy("a = x**2+2*x+1")
sympy("b = (x+1)**2")
"0" == sympy("simplify(a-b)")  # if they are equal, the result is zero

# urceni korenu rovnice
sympy("solve(x**2 - 2, x)")  # solve x^2-2=0


# zjednoduseni vyrazu
sympy("simplify((x**3 + x**2 - x - 1)/(x**2 + 2*x + 1))")
sympy("factor(x**3 - x**2 + x - 1)")
sympy("(x*y**2 - 2*x*y*z + x*z**2 + y**2 - 2*y*z + z**2)/(x**2 - 1)")

x <- Var("x")
y <- Var("y")
z <- Var("z")
sympy("collect(x*y + x - 3 + 2*x**2 - z*x**2 + x**3, x)")  # organize equation around var 'x'
sympy("cancel((x*y**2 - 2*x*y*z + x*z**2 + y**2 - 2*y*z + z**2)/(x**2 - 1))")
sympy("simplify((x*y**2 - 2*x*y*z + x*z**2 + y**2 - 2*y*z + z**2)/(x**2 - 1))")


# roznasobeni zavorek
sympy("(x + 2)*(x - 3)")
sympy("expand((x + 2)*(x - 3))")
sympy("expand_func(gamma(x + 3))") # pro funkci


# limita funkce
sympy("limit(1/x, x, oo)")   # limit eg, x -> Inf
sympy("limit(1/x, x, 0)") 
sympy("limit(sin(2*x)/(2*x), x, 0)") 
sympy("limit((1+(1/(3*x)))**(3*x), x, oo)") 
sympy("limit((exp(x)-E)/(x-1), x, 1)") 
sympy("limit((ln(x)-1)/(x-E), x, E)") 
sympy("limit((exp(x)-exp(-x))/sin(x), x, 0)") 


### Pri redení daného mnozství kyseliny sírové pridáním vody o hmotnosti x vznikne mnozství tepla y = Ax / x+B , kde A, B jsou kladné konstanty. 
# Jaké mnozství tepla se uvolní pro x -> oo ? Kolik vody se musí pridat aby vzniklo 80 % tohoto mnozství tepla ?
x = Var("x")
A = Var("A")
B = Var("B")
sympy("limit(A*x/(x+B), x, oo)")


### Okamzitá koncentrace látky C, vzniklé pri nevratné chemické reakci A + B = C, je c(t) = ab(1-e^(b-a)kt) / a-bee^(b-a)kt , kde a, b, k jsou kladné konstanty, t > 0.
# Jaké hodnoty nabývá koncentrace c pro t -> oo v prípadech a > b, b > a. 
a = Var("a")
b = Var("b")
k = Var("k")
t = Var("t")
sympy("limit(a*b*(1-exp((b-a)*k*t))/(a-b*exp((b-a)*k*t)), t, oo)") 


### Molární tepelná kapacita cV dvouatomového ideálního plynu pri konstantním objemu V je dána vztahem cV(T) = 5R/2 + (hv/kT)^2 (R exp(hv/kT))/(exp(hv/kT)-1)^2
R = Var("R")
h = Var("h")
v = Var("v")
k = Var("k")
Tk = Var("Tk")
sympy("limit(R*5/2 + ((h*v/(k*Tk))**2)*(R*exp(h*v/(k*Tk))/(exp(h*v/(k*Tk))-1)**2), Tk, oo)") 


### Molární tepelná kapacita cV závisí na teplote T podle vztahu cV(T) = 3R(hv/kT)^2 (exp(hv/kT))/(exp(hv/kT)-1)^2
R = Var("R")
h = Var("h")
v = Var("v")
k = Var("k")
Tk = Var("Tk")
sympy("limit(3*R*((h*v/(k*Tk))**2)*(exp(h*v/(k*Tk))/(exp(h*v/(k*Tk))-1)**2), Tk, oo)") 




sympy("y = x*x") # create a variable 'y' in the SymPy persistant state

sympy("A = Matrix([[1,x], [y,1]])")
sympy("A**2")

sympy("B = A.subs(x,1.1)")  # replace x by 1.1 (btw, SymPy objects are immutable)
sympy("B**2")

# more replacement, a subexpression by another:
sympy("expr = sin(2*x) + cos(2*x)")
sympy("expr.subs(sin(2*x), 2*sin(x)*cos(x))")
sympy("expr.subs(x,pi/2)")

x <- Var("x")
y <- Var("y")
(x+y)*(x+y)



### Derivace ###

D(expression(sin(2*x)), "x") # also possible in base R

sympy("diff(sin(2*x), x, 1)")  # first derivative
sympy("diff(sin(2*x), x, 2)")  # second derivative
sympy("diff(sin(2*x), x, 3)")  # third  derivative
sympy("diff(exp(x*y*z), x, y, y)") # d^3/dxdy^2
sympy("diff(exp(x*y*z), x, z, 3)") # d^4/dxdz^3

library(Deriv)
Deriv("sin(2*x)", "x", cache.exp=FALSE) # differentiate only by x 


### Závislost mnozství m látky získané chemickou reakcí v case t je dána vztahem m = A(1-exp(-kt)), kde A a k jsou konstanty. Jak závisí rychlost reakce na case t ?
D(expression(A*(1-exp(-k*t))), "t") # base R
library(rSymPy)
sympy("diff(A*(1-exp(-k*t)), t, 1)")
library(Deriv)
Deriv("A*(1-exp(-k*t))", "t", cache.exp=FALSE)  


### Pri chemické reakci A + B = C v níz obe pocátecní koncentrace výchozích látek nabývají stejné hodnoty a, je závislost koncentrace x vznikající látky na case dána vztahem x(t) = (a^2)kt / 1+akt. Jak se mení koncentrace s casem, je-li k > 1 ? kt . 
D(expression((a^2)*k*t/(1+a*k*t)), "t") # base R
library(rSymPy)
sympy("diff((a**2)*k*t/(1+a*k*t), t, 1)")
library(Deriv)
Deriv("(a^2)*k*t/(1+a*k*t)", "t", cache.exp=FALSE)  



library(Deriv)
Deriv("sin(x^2) * y", "x") # differentiate only by x 
Deriv("sin(x^2) * y", "y") # differentiate only by x 
Deriv("sin(x^2) * y", cache.exp=FALSE) # differentiate by all variables (here by x and y) 
Deriv("sin(x^2) * y", cache.exp=TRUE) # differentiate by all variables (here by x and y) 
Deriv("sin(cos(x + y^2))", cache.exp=FALSE) # differentiate by all variables (here by x and y) 
Deriv("sin(cos(x + y^2))", "x")
Deriv("sin(cos(x + y^2))", "y")
# Compound function example (here abs(x) smoothed near 0) 
fc <- function(x, h=0.1) if (abs(x) < h) 0.5*h*(x/h)**2 else abs(x)-0.5*h 
Deriv("fc(x)", "x", cache.exp=FALSE) 

Simplify("x^3 - x^2 + x - 1")
Simplify("x*y + x - 3 + 2*x^2 - z*x^2 + x^3")
Simplify("(x*y**2 - 2*x*y*z + x*z**2 + y**2 - 2*y*z + z**2)/(x**2 - 1)")



### Integral ###

integrate(function(x) exp(-x), 0, Inf)  # integration is possible in R

sympy("integrate(exp(-x))")    # indefinite integral
sympy("integrate(exp(-x*y),x)")  # indefinite integral  
sympy("integrate(exp(-x), (x, 0, oo))") # definite integral
sympy("integrate(x**2 - y, (x, -5, 5), (y, -pi, pi))")  # definite integral


### Mnozství tepla Q (v J) potrebné na zahrátí 10 kg zeleza o 1 K je pro T = <273.15, 473.15> vyjádrena funkcí Q = 4.1686 (105.3T + 0.42T^2). Vypocítejte strední mnozství tepla (v J / kg K) potrebného na zahrátí zeleza z teploty 293.15 K na 373.15 K.  
qm = integrate(function(Tk) 4.1868*(105.3*Tk + 0.142*Tk^2), 293.15, 373.15)  # integration is possible in R
unlist(qm)
0.1*unlist(qm)$value/(373.15-293.15) # [J / kg K]

Tk = Var("Tk")
smp = sympy("integrate(4.1868*(105.3*Tk + 0.142*Tk**2), (Tk, 293.15, 373.15))")  # definite integral
0.1*as.numeric(smp)/(373.15-293.15) # [J / kg K]


### Merné teplo c benzínu (pri konstantním tlaku) v závislosti na teplote T je vyjádrené funkcí c = -0.06 + 0.0010228T. Vypocítejte strední merné teplo benzínu pro T = <389.15, 491.15>.
qm = integrate(function(Tk) -0.06 + 0.0010228*Tk, 389.15, 491.15)  # integration is possible in R
unlist(qm)
0.1*unlist(qm)$value/(491.15-389.15) # [J / kg K]

Tk = Var("Tk")
smp = sympy("integrate(-0.06 + 0.0010228*Tk, (Tk, 389.15, 491.15))") 
as.numeric(smp)/(491.15-389.15) # [J / K]


##################################################################################################################

library(Ryacas)
# https://www.andrewheiss.com/blog/2019/02/16/algebra-calculus-r-yacas/

yac_str("N(Pi, 50)")

yac_str("x+x+x+x")
yac_expr("x+x+x+x")
eval(yac_expr("x+x+x+x"), list(x = 4))

yac_str("Factor(x^2+x-6)")
yac_expr("Factor(x^2+x-6)")

eq <- "x^2 + 4 + 2*x + 2*x"
yac_str(eq) # Evaluate yacas command x (a string) and get result as string/character. No task was given to yacas, so we simply get the same returned

# zjednodusení výrazu
yac_str(paste0("Simplify(", eq, ")"))  # Simplify(x): Simplify an expression

eq <- "x^2 + (x + 14)^2 - 34^2" 
rr = yac_str(paste0("Simplify(", eq, ")"))  # zjednodusení výrazu
x = 2
eval(parse(text="x^2+14*x-480"))

yac_str(paste0("Factor(", eq, ")"))    # Factor(x): Factorise an expression
yac_str(y_fn(eq, "Factor"))

y_fn(eq, "Factor")
yac_expr(paste0("Factor(", eq, ")"))

zz = y_fn(eq, "Factor")
yac_expr(paste0("Expand(", zz, ")"))

#
expr <- yac_expr(paste0("Factor(", eq, ")"))
expr
eval(expr, list(x = 2))


x <- ysym("x")
2*x^2 - 5
c(-2, 5)*x
c(2*x, -x^3)
as_r(c(-2, 5)*x) # or yac_expr(c(-2, 5)*x)


xs <- ysym("x")
poly <- xs^2 - xs - 6
poly
zeroes <- solve(poly, "x") # Solve(x^2 - x - 6 == 0, x)
zeroes

solve(poly, 3, "x") # Solve(x^2 - x - 6 == 3, x)


x <- ysym("x")
y <- ysym("y")
lhs <- c(3*x*y - y, x)
rhs <- c(-5*x, y+4)
sol <- solve(lhs, rhs, c("x", "y"))
sol
sol_vals <- lapply(seq_len(nrow(sol)), function(sol_no){y_rmvars(sol[sol_no, ])})
sol_vals
sol_envir <- lapply(sol_vals, function(l){list(x = as_r(l[1]), y = as_r(l[2]))})
sol_envir
do.call(rbind, lapply(seq_along(sol_envir), function(sol_no) {sol_val <- sol_envir[[sol_no]]
data.frame(sol_no = sol_no,
           eq_no = seq_along(sol_val),
           lhs = eval(as_r(lhs), sol_val),
           rhs = eval(as_r(rhs), sol_val))}))


# Rosenbrock function:
x <- ysym("x")
y <- ysym("y")
f <- (1 - x)^2 + 100*(y - x^2)^2
f

N <- 30
x <- seq(-1, 2, length=N)
y <- seq(-1, 2, length=N)
f_r <- as_r(f)
f_r
z <- outer(x, y, function(x, y) eval(f_r, list(x = x, y = y)))
levels <- c(0.001, .1, .3, 1:5, 10, 20, 30, 40, 50, 60, 80, 100, 500, 1000)
cols <- rainbow(length(levels))
contour(x, y, z, levels = levels, col = cols)
g <- deriv(f, c("x", "y")) # minimum
g
crit_sol_all <- solve(g, c("x", "y"))
crit_sol_all


Solve(expr, var) solve an equation (refer to the Ryacas function y_rmvars())
Variables(): List yacas variables
# yacas( "OldSolve({a*x+y==0,x+z==0},{x,y})" )
#  equation $a x = b$ for $x$.
solve(a,b)   # ???nd roots of a for variable b, i.e. yacas Solve(a == 0,b) 
solve(a,b,v) #  ???nd solutions to a == b for variable v, i.e. yacas Solve(a == b,v) 
#  also works for a system of equations (when a is a vector)
A <- outer(0:3, 1:4, "-") + diag(2:5) 
a <- 1:4 
B <- ysym(A) 
b <- ysym(a) solve(A) 
solve(B) 
solve(A, a) 
solve(B, b)
poly <- ysym("x^2 - x - 6") 
solve(poly, "x")   # Solve(poly == 0, x) 
solve(poly, 3, "x") # Solve(poly == 3, x)


sum(expr, var, lower, upper, ..., na.rm = FALSE)

## D(x) expr: Take the derivative of expr with respect to x

yac_str("D(x) x^2+x-6")
yac_expr("D(x) x^2+x-6")

L <- ysym("x^2 * (y/4) - a*(3*x + 3*y/2 - 45)")
L
deriv(L, "x")

my_log <- function(x){return(sin(log(2 + x)))}
deriv(my_log(x), x)
integrate(my_log(x), x)
integrate(f, var, lower, upper )

# Hessian(function, var) # Create the Hessian matrix
Hessian(L, "x")
deriv(L, c("x", "y", "a"))
H <- Hessian(L, c("x", "y", "a")); H
as_r(H)
eval(as_r(H), list(x = 2, y = 2, a = 2))
# Jacobian(function, var) # Create the Jacobian matrix
L2 <- ysym(c("x^2 * (y/4) - a*(3*x + 3*y/2 - 45)", "x^3 + 4*a^2")) # just some function
L2
Jacobian(L2, "x")
Jacobian(L2, c("x", "y", "a"))


# lim(f, var, val, from_left, from_right) # Limit of f(n) for n going towards a (e.g. Infinity or 0)
cmd <- "Limit(n, Infinity) (1+(1/n))^n"
yac_str(cmd)
yac_expr(cmd)

yac_str("Limit(h, 0) (Sin(x+h)-Sin(x))/h") 


# sum(k, a, b, f(k)) # Sum of f(k) for k from a to b # urcity integral
yac_str("Sum(k, 0, n, a^k)")
yac_expr("Sum(k, 0, n, a^k)")


yac_str("poly := (x-3)*(x+2)")
yac_silent("poly := (x-3)*(x+2)")
yac_str("Expand(poly)")


################################################################################################################


### Spojovani retezcu

paste("The value of pi:", pi)

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

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

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


cat(my_string)
cat(my_string, "with R")
cat(my_string, "with R", sep = " =) ")

cat(1:10, sep = "-")

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

cat("Loooooooooong strings", "can be displayed", "in a nice format","by using the 'fill' argument")

cat(my_string, "with R", file = "output.txt")


#####

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]


##### Format retezce #######

format(13.7)
format(13.12345678)
format(13.7, nsmall = 3)


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

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

format(123456789, big.mark = ",")


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


#####################

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

toString(c("one", "two", "3333333333"), width = 8)
toString(c("one", "two", "3333333333"), width = 12)


#####################

length("How many characters?")
length(c("How", "many", "characters?"))

nchar("How many characters?")
nchar(c("How", "many", "characters?"))


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

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


# replace 'a' by 'A'
chartr("a", "A", "This is a boring string")
chartr("ai", "X", "This is a bad example")

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

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

# replace 2nd letter with hash symbol
x = c("may", "the", "force", "be", "with", "you")
substr(x, 2, 2) <- "#"
x

# replace 2nd and 3rd letters with happy face
y = c("may", "the", "force", "be", "with", "you")
substr(y, 2, 3) <- ":)"
y

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


substring("ABCDEF", 2, 4)

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

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


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


###################################################

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

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

#######################################################

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)


#######################################################


startsWith(month.name, "J")

endsWith(month.name, "ember")


#######################################################

substr(month.name, 1, 3)

substring(month.name, 1, 3)

strtrim(month.name, 3)



trimws(" This has trailing spaces.  ")

#######################################################


library(stringr)

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

# removing zero length objects
paste("May", "The", "Force", NULL, "Be", "With", "You", character(0))
str_c("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


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


# some factor
some_factor = factor(c(1, 1, 1, 2, 2, 2), labels = c("good", "bad"))
some_factor

nchar(some_factor)
str_length(some_factor)


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, -seq_len(nchar(lorem)))
substring(lorem, -seq_len(nchar(lorem)))

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

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

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


##################

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)


##################

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 = "-")


###################

bad_text = c("This", " example ", "has several   ", "   whitespaces ")

str_trim(bad_text, side = "left")

str_trim(bad_text, side = "right")

str_trim(bad_text, side = "both")


###################

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)


#####################


library(stringi)

stri_reverse("dna")


library(Biostrings)

str_rev("dna")