### kvadraticke a kubicke rovnice ###

dc = c(1,14,-480)

# Descartova veta: Pocet kladnych korenu polynomu a0*x^n + a1*x^n-1 + ... + an-1*x + an = 0 je bud roven poctu znamenkovych zmen v posloupnosti a0, a1, ..., an-1, an jeho koeficientu, nebo je o sudy pocet mensi. Pocet zapornych korenu se urci obdobne po dosazeni -x za x.
sum(diff(sign(dc)) != 0)

library(sfsmisc)
nr.sign.chg(dc)

# Horner
library(sfsmisc)
polyn.eval(dc, x)
library(cmna)
horner(x, coefs)
rhorner(x, coefs)
naivepoly(x, coefs)
betterpoly(x, coefs)

### 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
F2 = F1 + 14
F2

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)$root 
F1
F2 = F1 + 14
F2

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

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

# pro kontrolu lze pouzit Vietovy vzorce
uniroot(fun, c(-1, 0),extendInt = "yes", tol = 1e-9)$root
(-dc[2]/dc[1]) - F1
(dc[3]/dc[1])/F1
# nebo
library(pracma)
Poly(c(16, -30))

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


library(ssanv)
uniroot.integer(f, interval, lower = min(interval), upper = max(interval),step.power = 6, step.up = TRUE, pos.side = FALSE, print.steps = FALSE,maxiter = 1000)


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

library(animation) # animation
bisection.method(FUN = function(x) x^2 - 4, rg = c(-1, 10), tol = 0.001, interact = FALSE)
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=.01, x=10, tol = 0.001, m = 100)



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

dc = c(1,14,-480)

# Descartova veta: Pocet kladnych korenu polynomu a0*x^n + a1*x^n-1 + ... + an-1*x + an = 0 je bud roven poctu znamenkovych zmen v posloupnosti a0, a1, ..., an-1, an jeho koeficientu, nebo je o sudy pocet mensi. Pocet zapornych korenu se urci obdobne po dosazeni -x za x.
sum(diff(sign(dc)) != 0)

library(sfsmisc)
nr.sign.chg(dc)

# Horner
library(sfsmisc)
polyn.eval(dc, x)
library(cmna)
horner(x, coefs)
rhorner(x, coefs)
naivepoly(x, coefs)
betterpoly(x, coefs)

### 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
F2 = F1 + 14
F2

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)$root 
F1
F2 = F1 + 14
F2

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

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

# pro kontrolu lze pouzit Vietovy vzorce
uniroot(fun, c(-1, 0),extendInt = "yes", tol = 1e-9)$root
(-dc[2]/dc[1]) - F1
(dc[3]/dc[1])/F1
# nebo
library(pracma)
Poly(c(16, -30))

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


library(ssanv)
uniroot.integer(f, interval, lower = min(interval), upper = max(interval),step.power = 6, step.up = TRUE, pos.side = FALSE, print.steps = FALSE,maxiter = 1000)



### 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]
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
# cgmmatrix(A, B, tol = 1e-06, maxiter = 100) # nefunguje
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)")


library(rSymPy)
### 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


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



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


# Taylor expansion:
yac_str("texp := Taylor(x,0,5) Exp(x)")
yac_str("texp := Taylor(x,0,3) 1/Cos(x)")
yac_str("texp := Taylor(x,0,6) Ln(1+x)")
yac_str("texp := Taylor(x,0,4) x/(4-x*x)")

eval(yac_expr("texp := Taylor(x,0,6) Ln(1+x)"), list(x = 2))


### Mame smes kysliku a oxidu dusnateho. Pri jake koncentraci kysliku (v %) ve smesi se oxid dusnaty oxiduje nejrychleji. 
# 2 NO + O2 = 2 NO2
# v = k*(x^2)*y = k*(x^2)*(100-x)
library(Ryacas)
L <- ysym("(x^2)*(100-x)"); L
der <- deriv(L, "x"); der
as_r(der)
solve(der, "x")
round(200/3,1) # [%]

### Elektricky obvod se sklada ze zdroje stejnosmerneho elektromotorickeho napeti Ue s vnitrnim odporem Rv a seriove pripojeneho spotrebice o odporu R.Pri jake hodnote odporu R bude prikon spotrebice maximalni?    
# P = R*I^2 = R*(Ue/(R+Rv))^2
library(Ryacas)
L <- ysym("R/(R+Rv)^2"); L
der <- deriv(L, "R"); der
solve(der, "R")
yac_str(paste0("Simplify(",sqrt(abs(-R^2)),")"))  # Simplify(x): Simplify an expression


### Urcete velikost hydrostaticke tlakove sily, ktera pusobi na plast valce vysky v a polomeru r, zcela naplneneho kapalinou o hustote rho. 
# dF = p*dS = p*2*pi*r*dh = rho*g*h*2*pi*r*dh
library(Ryacas)
L <- ysym("2*pi*r*rho*g*h*dh"); L
int = integrate(L, "h"); int
as_r(int)


############### vektory ###############################################################

x = sample(12:99, 9, replace = TRUE)


### norma vektoru

library(pracma)
Norm(x, p = 2) 

library(fBasics)
norm2(x, p = 2)

library(cmna)
vecnorm(x)

# Euklidovska norma (p = 2)
nv <- sum(abs(x)^2)^(1/2)
nv <-sqrt(t(x)%*%x)
norm_vec <- function(x) sqrt(sum(x^2))
norm_vec <- function(x){sqrt(crossprod(x))}

# Pythagorova veta - delka prepony pravouhleho trojuhelnika s odvesnami dilky 3 a 4, tez Euklidovska vzdalenost v CA.
norm(as.matrix(c(3, 4)), type = "F")
Norm(c(3, 4))          
norm_vec(c(3, 4))


### uhel 2 vektoru

library(matlib)
matA <- matrix(c(3, 1), nrow = 2)  ##column vectors
matB <- matrix(c(5, 5), nrow = 2)
angle(as.vector(matA), as.vector(matB),degree = FALSE) # radiany
angle(as.vector(matA), as.vector(matB),degree = TRUE)  # stupne
corner(p1=matrix(c(0,0), nrow = 2), p2=matA, p3=matB, d = 0.1, absolute = TRUE)

printMatrix(M, parent = TRUE, fractions = FALSE, latex = FALSE, tol = sqrt(.Machine$double.eps) )
printMatrix(M, parent = TRUE, fractions = TRUE, latex = TRUE, tol = sqrt(.Machine$double.eps) )



angle2v <- function(x,y){
  dot.prod <- x%*%y 
  norm.x <- norm(x,type="2")
  norm.y <- norm(y,type="2")
  theta <- acos(dot.prod / (norm.x * norm.y))
  as.numeric(theta)
}
ang = angle2v(t(matA),matB); ang  # radiany
library(pracma)
rad2deg(ang)                      # stupne


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

# delka prostorove diagonaly v jednotkovi krychli
sqrt(dot(c(1,1,1), c(1,1,1)))  #=> 1.732051
# v}kon: soucin smly a rychlosti
# mechanicka prace: soucin smly a drahy


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

A <- matrix(c(1,0,0, 0,1,0), nrow=2, ncol=3, byrow=TRUE)
crossn(A)  #=> 0 0 1

# https://stackoverflow.com/questions/15162741/what-is-rs-crossproduct-function
# https://stackoverflow.com/questions/36798301/r-compute-cross-product-of-vectors-physics


# Lorentzova sila
# Pohyb tuheho telesa
# uhlovy moment hybnosti



############### matice ##################################################################

#define a matrix for this example
M <- t(matrix(data = rnorm(12), ncol = 3))

matrix(1,3,3)
matrix(0,3,3)
matrix(1:12,ncol = 4)
matrix(1:12,ncol = 4, byrow=TRUE)
matrix(1:12,ncol = 4, byrow=FALSE)
matrix(1:12,nrow = 3, byrow=TRUE)
matrix(1:12,nrow = 3, byrow=FALSE)
mm = matrix(1:12,nrow = 3); mm
mm[,c(4:1)]
mm[c(3:1),]
mm[,1:3] = mm[,4]; mm
t(matrix(1:12,nrow = 3)) # transponovana matice

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)
as.vector(M)
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 konstitucnich koeficientu
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


### Urcete pocet linearne nezavislych reakci v soustave obsahujici dane slozky.
# Gibbsovo stechiometricke pravidlo: maximalni pocet linearne nezavislych reakci r je mensi nebo roven rozdilu poctu slozek v systemu (n) a hodnosti matice konstitucnich koeficientu (h), v n-slozkovem systemu existuje n - h stechiometrickych vztahu.

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

# spocitat pomoci vyse uvedene procedury
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)


### Vycislete rovnici: Fe2O3 + Al = Al2O3 + Fe
Fe = c(2, 0, 0, -1)
O = c(3, 0, -3, 0)
Al = c(0, 1, -2, 0)
# http://ndu2009algebra.blogspot.com/2011/04/as-we-struggle-to-increase-number-of.html
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)
#limSolve::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(Rlinsolve)
ls = lsolve.gs(A, B, xinit = NA, reltol = 1e-05, maxiter = 1000, adjsym = TRUE, verbose = TRUE)
ls$x # nelze pouzit (underdetermined system)

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


### TNT vznika reakci toluenu s kyselinou dusicnou: xC7H8 + yHNO3 = zC7H5O6N3 + wH2O. Vycislete rovnici.

### Nitroglycerin (trinitroglycerol) C3H5N3O9 se explozi rozklada na dusik, oxid uhlicity, vodu a kyslik, za vzniku velkeho mnozstvi energie. Vycislete rovnici:
# vC3H5N3O9 = wN2 + xCO2 + yH2O + zO2

### Vycislete rovnici fotosyntezy: xCO2 + yH2O = zC6H12O6 + wO2


### Vycislete rovnici: C3H8 + O2 = CO2 + H2O
C = c(3, 0, -1, 0)
H = c(8, 0, 0, -2)
O = c(0, 2, -2, -1)
A = as.matrix(rbind(C, H, O))
colnames(A) = c("C3H8","O2","CO2","H2O")
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)
#limSolve::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))]
mx = cmna::rrefmatrix(cbind(A,-B)) #  returns the reduced row echelon matrix
rownames(NN) = rownames(A)
colnames(NN) = colnames(A)[-length(colnames(A))]
apply(NN,2,sum)
library(Rlinsolve)
ls = lsolve.gs(A, B, xinit = NA, reltol = 1e-05, maxiter = 1000, adjsym = TRUE, verbose = TRUE)
ls$x # nelze pouzit (underdetermined system)


### Vycislete rovnici: KMnO4 + MnSO4 + H2O = MnO2 +  K2SO4 + H2SO4
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)
# https://www.physicsforums.com/threads/balancing-chemical-equations-with-linear-algebra.374530/
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)
#limSolve::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)
library(Rlinsolve)
ls = lsolve.gs(A, B, xinit = NA, reltol = 1e-05, maxiter = 1000, adjsym = TRUE, verbose = TRUE)
ls$x # nelze pouzit (underdetermined system)


### Kolik 96% kyseliny sirove a kolik vody je potreba na pripravu 1 l 78% kyseliny sirove?
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 kazdi slitiny A, B, C je treba na v}robu 100 kg slitiny obsahujmcm 2 % Si, 1 % Mn, 0.4 % P a 0.15 % S?
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


### Jedno kapalni hnojivo obsahuje 20 % dusmku, 10 % dusmku a 5 % draslmku, druhi kapalni hnojivo obsahuje 15 % dusmku a 20 % fosforu. Kolik prvnmho a kolik druhiho hnojiva je treba na prmpravu 10 l hnojiva, kteri by obsahovalo
# a) 19 % dusmku, 12 % fosforu a 4 % draslmku
# b) 28 % dusmku, 20 % fosforu a 4 % draslmku
hn1 = c(0.2, 0.1, 0.05)
hn2 = c(0.15, 0.2, 0)
A = cbind(hn1,hn2); rownames(A) = c("N","P","K")
B = c(0.19, 0.12, 0.04); names(B) = c("N","P","K")
B = c(0.28, 0.20, 0.04); names(B) = c("N","P","K")
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[,3]*100 # [kg]
library(Rlinsolve)
ls = lsolve.gs(A, B, xinit = NA, reltol = 1e-05, maxiter = 1000, adjsym = TRUE, verbose = TRUE)
ls$x


### Surovina 1 pro magnezitovi stavivo obsahuje 80 % MgO, surovina 2 pro chromitovi stavivo obsahuje 50 % Cr2O3. Kolik kg surovin 1 a 2 potrebujeme na v}robu 100 kg magnezitovo-chromitoviho staviva, kteri by obsahovalo 48 % MgO a 20 % Cr2O3?
sur1 = c(0.8, 0)
sur2 = c(0, 0.5)
A = cbind(sur1,sur2); rownames(A) = c("MgO","Cr2O3")
B = c(0.48, 0.2); names(B) = c("MgO","Cr2O3")
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)
GE1 = gaussianElimination(A, B, tol = sqrt(.Machine$double.eps), verbose = FALSE, latex = FALSE, fractions = FALSE)
GE1[,3]*100 # [kg]
library(Rlinsolve)
ls = lsolve.gs(A, B, xinit = NA, reltol = 1e-05, maxiter = 1000, adjsym = TRUE, verbose = TRUE)
ls$x


### Doporucena dennm davka vitammnu B2 je 1.4 mg, vitammnu C 80 mg a vitammnu K 0.075 mg. Kolik porcm mrkve, spenatu a rajcat je treba snmst, abyste dodrzeli doporucenou dennm davku vitammnu B2, C a K. V jedni porci mrkve je 0.059 mg B2. 1.5 mg C a 0 mg K, v jedni porci spenatu je 0.236 mg B2, 9.8 mg C a 0.145 mg K a v jedni porci rajcat je 0.148 mg B2, 19.1 mg C a 0 mg K. 
B = c(1.4,80,0.075) # [mg]
names(B) = c("B2","C","K")
mrkev = c(0.059,1.5,0) # [mg]
spenat = c(0.236,9.8,0.145) # [mg]
rajcata = c(0.148,19.1,0) # [mg]
A = cbind(mrkev,spenat,rajcata)
rownames(A) = c("B2","C","K")
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)
# jen nulove a/nebo kladne koeficienty
library(limSolve)
G <-matrix(ncol = 3, byrow = TRUE, data = c(1,0,0, 0,1,0, 0,0,1)) 
H <- c(0, 0, 0)
ldei(A, B, G = G, H = H)$X


### Kadernice potrebuje pripravit 12% roztok peroxidu vodiku. V jakem pomeru musime smichat 3% a 30% roztok peroxidu?  dvojnasobek 3%
B = c(0.12,1) # [mg]
names(B) = c("H2O2","H2O")
r1 = c(0.3,1) # [mg]
r2 = c(0.03,1) # [mg]
A = cbind(r1,r2)
colnames(A) = c("H2O2","H2O")
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


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

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

A <- matrix(c(1,-1,-5, 8,6,3, 3,-7,4), nrow=3, ncol=3, byrow=TRUE)
B <-  c(0, 2, 3)

### Gaussova eliminace
library(matlib)
matlib::gaussianElimination(A, B, tol = sqrt(.Machine$double.eps), verbose = FALSE, latex = FALSE, fractions = TRUE)
matlib::Solve(A, B)

library(cmna)
cmna::refmatrix(cbind(A,-B)) # reduces a matrix to row echelon form
cmna::rrefmatrix(cbind(A,-B)) #  returns the reduced row echelon matrix
cmna::solvematrix(A, B)

library(Rlinsolve)
out1 = lsolve.cg(A,B)
out2 = lsolve.bicg(A,B)
out3 = lsolve.bicgstab(A,B)
out4 = lsolve.sor(A,B,w=0.5)
out5 = lsolve.cgs(A,B)
out6 = lsolve.cheby(A,B)
matout = cbind(out1$x, out2$x, out3$x, out4$x, out5$x, out6$x)
colnames(matout) = c("CG result", "BiCG result", "BiCGSTAB result","SSOR result","CGS result", "Chebyshev result") 
print(matout)

library(limSolve)
limSolve::Solve(A, B)

library(RandomFieldsUtils)
RandomFieldsUtils::solvex(A, B, logdeterminant=FALSE)


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

# Linear algebra 
A <- outer(0:3, 1:4, "-") + diag(2:5)  # numericka matice
a <- 1:4
B <- ysym(A)                           # yacas matice
B
b <- ysym(a)
b

as_r(B) # objekt yacas na objekt R 
as_y(A)


dim(A)
dim(B)

cbind(b,b)
rbind(b,b)

A[, 2:3]
B[, 2:3]

diag(A)
diag(B)
lower.tri(A, diag = FALSE)
lower.tri(B, diag = FALSE)
as_y(lower.tri(B, diag = FALSE))
upper.tri(A, diag = FALSE)
upper.tri(B, diag = FALSE)
as_y(upper.tri(B, diag = FALSE))
A[upper.tri(A)] <- 1
B[upper.tri(B)] <- 1

C <- B
C[lower.tri(C)] <- "x"
C[upper.tri(C)] <- "1"
diag(C) <- 2
C[2, 3] <- "ou"
C

exp(A)
exp(B)

2*A - A
2*B - B

A %*% a
B %*% b
B * B 
B %*% B

# transpose of a matrix
t(A)
y_fn(B, "Transpose")
t(B)

# inverse of a matrix
library(cmna)
invmatrix(A)

y_fn(B, "Inverse")
solve(B)
A %*% solve(A)
round(A %*% solve(A),0)
B %*% solve(B)
solve(A %*% t(A))
solve(B %*% t(B))
solve(A, a)
solve(B, b)

# trace of a matrix, soucet hodnot na hlavni dagonale
y_fn(B, "Trace")


### Kadernice potrebuje pripravit 12% roztok peroxidu vodíku. V jakém pomeru musí smíchat 3% a 30% roztok peroxidu?  dvojnásobek 3%
r1 = c(0.03, 1)
r2 = c(0.3, 1)
A = ysym(cbind(r1, r2))
B = ysym(c(0.12, 1))
solve(A,B)
as_r(solve(A,B))


### Ze dvou kovu o hustotách 7.4 g/cm3 a 8.2 g/cm3 je treba pripravit 0.5 kg slitiny o hustote 7.6 g/cm3. Kolik g kazdého z kovu je k tomu potreba? 375 g lehcího a 125 g tezsího
r1 = c(7.4, 1)
r2 = c(8.2, 1)
A = ysym(rbind(r1, r2))
B = ysym(c(7.6, 1))
solve(A,B)


### Kolik g 60% a kolik g 30% roztoku NaCl je treba smíchat pri príprave 100 g 40% roztoku? 20 g 60% a a 80 g 35%
r1 = c(0.60, 1)
r2 = c(0.30, 1)
A = rbind(r1, r2); rownames(A) = c("NaCl","H2O")
A = ysym(A)
B = c(0.40,100); names(B) = c("NaCl","H2O")
B = ysym(B)
solve(A,B)

r1 = c(0.3, 1)
r2 = c(0.6, 1)
A = ysym(cbind(r1, r2))
B = ysym(c(0.4, 1))
solve(A,B)
as_r(solve(A,B))


############# Linearni programovani ##############################################

library(lpSolveAPI)

### Farmar ma 75 akru pudy na nichz pestuje psenici a jecmen. Naklady na pestovanm (semena, hnojiva, apod.) jsou 120 $ na akr pro psenici a 210 $ na akr pro jecmen. Farmar ma k dispozici 15 000 $ a zlozn} prostor na 4000 bushelu obilm. Na 1 akru vypestuje v prumeru 110 bushelu psenice resp. 30 bushelu jecmene. Pokud je v}nos z 1 bushelu u psenice (po odectenm vsech nakladu) 1.30 $ a u jecmene 2.00 $, jak by mel farmar osmt 75 akru pudy, aby dosahl maximalnmho zisku?
# maximize   P = (110)(1.30)x + (30)(2.00)y = 143x + 60y
# subject to: 120x + 210y <= 15000, 110x + 30y <= 4000, x + y <= 75, x >= 0, y >= 0
lprec <- make.lp(0, 2)
lp.control(lprec, sense="max")
set.objfn(lprec, c(143, 60))
add.constraint(lprec, c(120, 210), "<=", 15000)
add.constraint(lprec, c(110, 30), "<=", 4000)
add.constraint(lprec, c(1, 1), "<=", 75)
lprec
solve(lprec)
get.objective(lprec) # Get maximum profit
get.variables(lprec) # Get the solution

#### Kombinatorika ####################################################################################

choose(n, k)   # kombinace
choose(3,2)
lchoose(n, k)  # log kombinace
factorial(x)   # faktorial x!
factorial(5)
lfactorial(x)  # faktorial log x!

perm = function(n, x){ factorial(n) / factorial(n-x)}
comb = function(n, x){ factorial(n) / factorial(n-x) / factorial(x)}


## kartezsky soucin mnozin
expand.grid(c(1, 2, 3), c(4, 5, 6))
# https://stackoverflow.com/questions/4309217/cartesian-product-data-frame


my_list <- c('one','two','three','four','five')
combn(my_list, 4, simplify=F)
combn(my_list, 4, simplify=T)
combn(c(1, 2, 3, 4, 5, 6), m=2, FUN = NULL, simplify = TRUE)
combn(5,3)


library(combinat)
permn(3)
length(permn(3))
combn(3, 2)
dim(combn(3,2))[2]


library(gtools)
n=5
r=3
combinations(n, r, v=1:n, set=TRUE, repeats.allowed=FALSE)
permutations(n, r, v=1:n, set=TRUE, repeats.allowed=FALSE)


### Urna obsahuje 3 koule ruzných barev, kolika zpusoby (s vracením) lze vybrat 2 koule z urny ?
x <- c('red', 'blue', 'black')
permutations(n=3,r=2,v=x,repeats.allowed=T)
nrow(permutations(n=3,r=2,v=x,repeats.allowed=T))
library(prob)
nsamp(3, 2, replace = TRUE, ordered = TRUE)
urnsamples(1:3, size = 2, replace = TRUE, ordered = TRUE)
nrow(urnsamples(1:3, size = 2, replace = TRUE, ordered = TRUE))


### Je treba zvolit ctyrmístný PIN. Kolik je mozno vytvorit ruznych kombinací císlic za predpokladu, ze se císlice nebudou opakovat?
n=10
r=4
permutations(n=10,r=4,v=1:n,repeats.allowed=F)
nrow(permutations(n=10,r=4,v=1:n,repeats.allowed=F))


### Kolik je trojciferných císel v nichz se zádná císlice neopakuje ?
n=10
r=3
mpm = permutations(n=10,r=3,v=0:9,repeats.allowed=F)
nrow(mpm)
which(mpm[,1]==0)
mpm = mpm[-which(mpm[,1]==0),]
nrow(mpm)


### Urcete celkový pocet kladných trojciferných císel, která obsahují císlici 6.
n=10
r=3
mpm = permutations(n=10,r=3,v=0:9,repeats.allowed=T)
nrow(mpm)
which(mpm[,1]==0)
mpm = mpm[-which(mpm[,1]==0),]
nrow(mpm)
length(unique(unlist(lapply(c(1:3),function(xx) which(mpm[,xx]==6)))))


### Kolik ruzných tetiv muze být vedeno 20 body po obvodu kruhu?
n=20
r=2
combinations(20, r=2, v=1:n, set=TRUE, repeats.allowed=FALSE)
choose(20, 2) 

### V sérii 12 výrobku jsou práve 3 vadné. Kolika zpusoby z nich lze vybrat 
# 6 libovolných výrobku
n=12
r=6
c1 = combinations(12, r=6, v=1:n, set=TRUE, repeats.allowed=FALSE)
cc1 = length(c1[,1]); cc1
choose(12, 6) 
# 6 výrobku bezvadných
c2 = combinations(9, r=6, v=1:n, set=TRUE, repeats.allowed=FALSE)
cc2 = length(c2[,1])
choose(9, 6) 
# 6 výrobku z nichz práve 1 je vadný
c3a = combinations(3, r=1, v=1:n, set=TRUE, repeats.allowed=FALSE)
length(c3a[,1])
c3b = combinations(9, r=5, v=1:n, set=TRUE, repeats.allowed=FALSE)
length(c3b[,1])
cc3 = length(c3a[,1])*length(c3b[,1])
choose(3, 1) 
choose(9, 5) 
# 6 výrobku z nichz práve 2 jsou vadné
c4a = combinations(3, r=2, v=1:n, set=TRUE, repeats.allowed=FALSE)
length(c4a[,1])
c4b = combinations(9, r=4, v=1:n, set=TRUE, repeats.allowed=FALSE)
length(c4b[,1])
cc4 = length(c4a[,1])*length(c4b[,1])
choose(3, 2) 
choose(9, 4) 
# 6 výrobku z nichz práve 3 jsou vadné
c5a = combinations(3, r=3, v=1:n, set=TRUE, repeats.allowed=FALSE)
length(c5a[,1])
c5b = combinations(9, r=3, v=1:n, set=TRUE, repeats.allowed=FALSE)
length(c5b[,1])
cc5 = length(c5a[,1])*length(c5b[,1])
choose(3, 3) 
choose(9, 3) 
# check - dle kombinatorického pravidla souctu
(cc2 + cc3 + cc4 + cc5) == cc1



### Pravdepodobnost
library(prob)
S <- rolldie(times = 3, makespace = TRUE)
Prob(S, X1+X2 > 9 )
Prob(S, X1+X2 > 9, given = X1+X2+X3 > 7)


###### STATISTIKA ##########################################################

library(openxlsx)
pisky <- read.xlsx("c:\\Users\\lubop\\Documents\\kurs\\DVpisky1.xlsx", sheet = "Pisky", startRow = 1, colNames = TRUE, rowNames = TRUE,detectDates = FALSE, skipEmptyRows = TRUE,rows = NULL, cols = NULL,check.names = FALSE)
summary(pisky)
pisky = pisky[,-1]
colnames(pisky)
x =  pisky[,"Sr"]

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

### arithmetic mean

mean(x, na.rm=TRUE)
apply(pisky,2,mean)

# trimming
mean(x,trim=0.1)
apply(pisky,2,function(xx) mean(xx,trim=0.1))
library(asbio)
trim.me(x, trim = 0.1) # trimmovany soubor


### rozptyl (variance)
# populace n / vyber n-1 viz Excel
var(x, na.rm=TRUE)
library(asbio)
var(trim.me(x, trim = 0.1)) # var trim. souboru

#### SD
sd(x, na.rm=TRUE)
sqrt(var(x, na.rm=TRUE))
sd(trim.me(x, trim = 0.1)) # sd trim souboru


### geometric mean
prod(x)^(1/length(x)) #compute the geometric mean
library(FSA)
geomean(x, na.rm = FALSE, zneg.rm = FALSE)

# lognormal mean
library(fuel)
lognormalmean(x, estimator="ml", base = exp(1), n = length(x), m = n - 1, d = 1/n)


#### geometric SD
library(FSA)
sg = geosd(x, na.rm = FALSE, zneg.rm = FALSE)
10^(sg)

#### lognormal SD
library(fuel)
lognormalsd(x, estimator="ml", base = exp(1),n = length(x),m = n - 1, d = 1/n)

### variacni rozpeti
range(x)
max(x)
min(x)
max(x)-min(x)

max(x)/min(x) 
# ratio exceeds 2 orders of magnitude, logarithmic plots will be informative
# if the ratio is between 1.5 and 2 orders of magnitude, logarithmically scaled plots will likely be informative;


#### Variacni koeficient
# CV > 100 % = plots on a logarithmic scale
# CV = 70% - 100%, the inspection of logarithmically scaled plots will likely be informative. 

sd(x)/mean(x)
library(DescTools)
CoefVar(x, weights = NULL, unbiased = FALSE, conf.level = 0,95, na.rm = FALSE)


### sikmost
library(moments)
skewness(x)# sikmost
skewness(x)/(sqrt(6/length(x))) # sikmost delena odhadem str. chyby
#hodnoty pod -2 a nad 2 indikují nenormalitu (Havranek)


### spicatost
library(moments)
kurtosis(x)# spicatost
kurtosis(x)-3
(kurtosis(x)-3)/(sqrt(24/length(x))) # sikmost delena odhadem str. chyby
#hodnoty pod -2 a nad 2 indikují nenormalitu (Havranek)


#### kvantily ####################################################################

### kvartily
Q1 = quantile(x,0.25)
Q2 = quantile(x,0.5)
Q3 = quantile(x,0.75)

quantile(x,c(0.25,0.5,0.75))


### median
median(x)
quantile(x, 0.5)


### interkvartilove rozpeti
IQR(x)
Q3-Q1


### absolutni medianova odchylka (MAD)
mad(x)



### SUMMARY ###########

# pro dany sloupec
summary(x)

library(pastecs)
stat.desc(x)

library(fields)
stats(x)

library(s20x)
summaryStats(x)

library(epiR)
epi.descriptives(x, conf.level = 0.95)


# pro celou matici
summary(pisky)

library(psych)
describe(pisky)

library(Hmisc)
describe(pisky)

library(pastecs)
stat.desc(pisky)

library(fields)
stats(pisky)


### konfidencni intervaly
t.test(x,conf.level=0.95)$conf.int

library(DescTools)
MeanCI(x,conf.level=0.95)

library(Rmisc)
CI(x,ci=0.95)


### bootstrap CI
library(DescTools)
BootCI(x, y = NULL, mean, bci.method = "norm", conf.level = 0.95, sides = "two.sided", R = 999)
# bci.method = c("norm", "basic", "stud", "perc", "bca")
# sides = c("two.sided", "left", "right")

library(boot)
Mboot = boot(x,function(x,i) mean(x[i]), R=10000)
mean(Mboot$t[,1])
boot.ci(Mboot,conf = 0.95,type = "norm")


#### Predikcni intervaly
library(predictionInterval)
pi.m(M=mean(x),SD=sd(x),n=length(x),rep.n=length(x))


#### Tolerancni intervaly
library(tolerance)
nptol.int(y,alpha=0.05,P=0.95,side=2,method="WALD",upper=NULL,lower=NULL)
# method = c("WILKS", "WALD", "HM","YM")


####################################################################################
######  GRAFY EDA ###################################################################
####################################################################################

library(MASS)
data(chem)
x = chem
# A numeric vector of 24 determinations of copper in wholemeal flour, in parts per million.

pisky <- read.xlsx("c:\\Users\\lubop\\Documents\\kurs\\DVpisky1.xlsx", sheet = "Pisky", startRow = 1, colNames = TRUE, rowNames = TRUE,detectDates = FALSE, skipEmptyRows = TRUE,rows = NULL, cols = NULL,check.names = FALSE)
summary(pisky)
pisky = pisky[,-1]
colnames(pisky)
x =  pisky[,"Ba"]

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

# Cullen Frey plot
library(fitdistrplus)
descdist(x, discrete=FALSE)
descdist(x, discrete=FALSE, boot=900)
descdist(trim.me(x, trim = 0.1), discrete=FALSE, boot=900) # asbio


### stripchart ############################################################################################################

stripchart(x,method="overplot",pch=1,col=1,cex=1.5)         # vertical=FALSE
stripchart(x,method="overplot",pch=1,col=1,cex=1.5,log="x")

stripchart(x,method="jitter",jitter=0.1, pch=1,col=1,cex=1.5)
stripchart(x,method="jitter",jitter=0.1, pch=1,col=1,cex=1.5,log="x")

stripchart(x,method="stack",pch=1,col=1,cex=1.5)
stripchart(x,method="stack",pch=1,col=1,cex=1.5,log="x")


library(beeswarm)
beeswarm(x, log = FALSE,pch = 1, col = 1, cex = 1.2, vertical=FALSE)
beeswarm(x, log = TRUE,pch = 1, col = 1, cex = 1.2, vertical=FALSE)



### boxplot #####################################################################

boxplot(x, horizontal=TRUE, notch = FALSE, xlab = "x")

boxplot(x, horizontal=TRUE, notch = TRUE, xlab = "x")

library(beeswarm)
boxplot(x, horizontal=TRUE, notch = FALSE, xlab = "x", log = "x", cex = 1.2)
beeswarm(x, log = TRUE,pch = 1, col = 1, cex = 1.2, vertical=FALSE, add=TRUE)

boxplot(x,xlab="x",las=2, horizontal=TRUE, cex=1.2,log="x") 
stripchart(x,pch=1, col=1, vert=FALSE, add=TRUE, cex=1.2,log="x")

boxplot.stats(x)
bs=boxplot.stats(x)
bs$out



### histogram #####################################################################

hist(x, freq=TRUE, breaks = 20, plot = TRUE,main = NULL) # pro counts

hist(x, freq=FALSE, breaks = 10, plot = TRUE,main = NULL)
rug(x)


#### Jadrovy odhad hustoty pravdepodobnosti (Kernel density estimator, KDE) ####################################################

den1 = density(x)
den1$bw
den1$x
den1$y
plot(den1$x, den1$y, type="l", ylab = "density", xlab = "x")
rug(x, lwd=2)

plot(density(x), main="")  #A bit bumpy
rug(x, lwd=2)

density(x,adjust = 10)
plot(density(x,adjust = 10), main="") #Very smooth
rug(x, lwd=2)

density(x,adjust = 0.1)
plot(density(x,adjust = 0.1), main="") #crazy bumpy
rug(x, lwd=2)


#### houslovy graf ##########################################################################################

library(vioplot)

vioplot(x, range=1.5, horizontal=TRUE, col=0)
stripchart(x, pch=1, col=2, vert=FALSE, add=TRUE, cex=1.2)

vioplot(x, range=1.5, horizontal=TRUE, col=0)
stripchart(x, method="jitter", pch=1, col=2, vert=FALSE, add=TRUE, cex=1.2)

vioplot(x, range=1.5, horizontal=TRUE, col=0)
beeswarm(x, log = FALSE,pch = 16, col = 2, cex = 1.2, vertical=FALSE, add=TRUE)


#### QQ graf ############################################################################################################################

my.qq <- qqnorm(sort(x), type = "p", xlab = "x", ylab = "Quantiles of the Normal Distribution",main=NULL)
qqline(x, col = "red")

library(car)
qqPlot(x, dist="norm",line="robust",envelope=.95)
# line=c("quartiles", "robust", "none")


#### EDA plots ################################################################################################

library(UsingR)
simple.eda(x)

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


### Testy na odlehlé body ###########################################################################

library(outliers)
dixon.test(x, type=22, opposite=FALSE, two.sided=TRUE)
# 10 pro 3-7, 11 pro 8-10, 21 pro 11-13, 22 pro 14 a vice hodnot
grubbs.test(x, type=10, opposite=FALSE, two.sided=TRUE)
grubbs.test(x, type=11, opposite=FALSE, two.sided=TRUE)
grubbs.test(x, type=20, opposite=FALSE, two.sided=TRUE)
# 10 pro 1 odl. bod, 11 pro 1 odl. bod na kazdem konci, 20 pro 2 odl. body na 1 konci

library(climtrends)
FindOutliersTietjenMooreTest(x,k=1,alpha=0.05)

library(EnvStats)
rosnerTest(x, k = 2, alpha = 0.05, warn = TRUE)


#### TESTY NORMALITY #####################################################################################

library(s20x)
normcheck(x, xlab = c("Theoretical Quantiles", ""),
          ylab = c("Sample Quantiles", ""), main = c("", ""), col = "light blue",
          bootstrap = TRUE, B = 5, bpch = 3, bcol = "lightgrey",
          shapiro.wilk = FALSE, whichPlot = 1:2, usePar = TRUE)


# Shapiro Francia test
library(nortest)
ShapiroFranciaTest(x)

# Shapiro Wilk test
shapiro.test(x)

# Anderson Darling test
library(nortest)
ad.test(x)

# DAgostinuv test sikmosti
library(moments)
agostino.test(x, alternative = "two.sided")
# alternative = c("two.sided", "less", "greater")

# Anscombe-Glynnuv test spicatosti
library(moments)
anscombe.test(x, alternative = "two.sided")
# alternative = c("two.sided", "less", "greater")


#### TESTY GOODNESS OF FIT #####################################################################################

# Anderson Darlinguv test
library(DescTools)
AndersonDarlingTest(x, "pnorm", mean=mean(x), sd=sd(x))

# Cramer von Mises
library(goftest)
cvm.test(x, "pnorm", mean=mean(x), sd=sd(x), estimated=TRUE)

# Kolmogorov Smirnov test
library(KScorrect)
LcKS(x, cdf="pnorm", nreps = 4999, G = 1:9, varModel = "E", parallel = FALSE, cores = NULL)
# varModel = c("E", "V")
# "pnorm" for normal,"pmixnorm" for (univariate) normal mixture,"plnorm" for lognormal (log-normal, log normal), "punif" for uniform, "plunif" for loguniform (log-uniform, log uniform),
library(dgof)
ks.test(x, y="pnorm",alternative = "two.sided", exact = NULL, tol=1e-8, simulate.p.value=FALSE, B=2000)
# alternative = c("two.sided", "less", "greater")


### Fitting distribution / Goodness of fit #################################################################

### normal
library(MASS)
dfp <- fitdistr(x, densfun="normal")
(dffl <- dfp$loglik)
dff <- dfp$estimate
mean <- as.numeric(dff[1]); mean
sd <- as.numeric(dff[2]); sd
mean(x)
sd(x)

ncd = seq(min(x), max(x), 0.05)
npd = dnorm(ncd, mean = mean, sd = sd)
hist(x,freq = FALSE)
lines(ncd,npd, type="l",col=2)
rug(x)

library(EnvStats)
norm.est <- enorm(x, method = "mvue", ci = TRUE, ci.type = "two-sided",ci.method = "exact", conf.level = 0.95, ci.param = "mean")
norm.m <- as.vector(norm.est[[3]])[1]
norm.sd <- as.vector(norm.est[[3]])[2]
###
sw.norm <- gofTest(x, test="sw", dist = "norm"); sw.norm
plot(sw.norm)
norm.W <- sw.norm[[6]]
norm.p <- sw.norm[[10]]


### lognormal
library(MASS)
dfp <- fitdistr(x, densfun="lognormal")
(dffl <- dfp$loglik)
dff <- dfp$estimate
meanlog <- as.numeric(dff[1])
sdlog <- as.numeric(dff[2])

library(EnvStats)
lnorm.est <- elnorm(x, method = "mvue", ci = TRUE, ci.type = "two-sided",ci.method = "exact", conf.level = 0.95)
lnorm.mlog <- as.vector(lnorm.est[[3]])[1]
lnorm.sdlog <- as.vector(lnorm.est[[3]])[2]
#
sw.lnorm <- gofTest(x, test="sw", dist = "lnorm",param.list = list(meanlog=lnorm.mlog, sdlog = lnorm.sdlog)); sw.lnorm
plot(sw.lnorm)
lnorm.W <- sw.lnorm[[6]]
lnorm.p <- sw.lnorm[[10]]


#### MULTIMODALITA #####################################################################################

##### modelovani kontaminovanych distribuci

# multimodalita 1
xfit<-seq(-20, 20, length=10000)
yfit1<-dnorm(xfit, mean=-5, sd=3)
yfit2<-dnorm(xfit, mean=5, sd=3)
yfits<-yfit1+yfit2
plot(xfit, yfit1, col=2, type="l",lty=2, ylim=c(0,(1.1*max(yfits))),xlab="x",ylab="pravdepodobnost")
lines(xfit, yfit2, col=4, lty=2)
lines(xfit, yfits, col=1)

# multimodalita 2
xfit<-seq(-10, 20, length=10000)
yfit1<-dnorm(xfit, mean=-2, sd=1)
yfit2<-dnorm(xfit, mean=5, sd=4)
yfits<-yfit1+yfit2
plot(xfit, yfit1, col=2, type="l",lty=2, ylim=c(0,(1.1*max(yfits))),xlab="x",ylab="pravdepodobnost")
lines(xfit, yfit2, col=4, lty=2)
lines(xfit, yfits, col=1)

# zmena sikmosti  1
xfit<-seq(-10, 20, length=10000)
yfit1<-dnorm(xfit, mean=0, sd=2)
yfit2<-dnorm(xfit, mean=4, sd=4)
yfits<-yfit1+yfit2
plot(xfit, yfit1, col=2, type="l",lty=2, ylim=c(0,(1.1*max(yfits))),xlab="x",ylab="pravdepodobnost")
lines(xfit, yfit2, col=4, lty=2)
lines(xfit, yfits, col=1)

# zmena sikmosti  2
xfit<-seq(-35,50, length=10000)
yfit1<-dnorm(xfit, mean=-6, sd=10)
yfit2<-dnorm(xfit, mean=12, sd=15)
yfits<-yfit1+yfit2
plot(xfit, yfit1, col=2, type="l",lty=2, ylim=c(0,(1.1*max(yfits))),xlab="x",ylab="pravdepodobnost")
lines(xfit, yfit2, col=4, lty=2)
lines(xfit, yfits, col=1)

# zmena spicatosti  1
xfit<-seq(-20, 20, length=10000)
yfit1<-dnorm(xfit, mean=-2.5, sd=4)
yfit2<-dnorm(xfit, mean=5, sd=4)
yfits<-yfit1+yfit2
plot(xfit, yfit1, col=2, type="l",lty=2, ylim=c(0,(1.1*max(yfits))),xlab="x",ylab="pravdepodobnost")
lines(xfit, yfit2, col=4, lty=2)
lines(xfit, yfits, col=1)

# zmena spicatosti 2
xfit<-seq(-8, 8, length=10000)
yfit1<-dnorm(xfit, mean=0, sd=0.8)
yfit2<-dnorm(xfit, mean=0, sd=2.4)
yfits<-yfit1+yfit2
plot(xfit, yfit1, col=2, type="l",lty=2, ylim=c(0,(1.1*max(yfits))),xlab="x",ylab="pravdepodobnost")
lines(xfit, yfit2, col=4, lty=2)
lines(xfit, yfits, col=1)

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

kyjov<- read.delim("c:\\Users\\lubop\\Documents\\kurs\\kyjov.txt", header=T)
x = kyjov[which(kyjov[,3]=="1034"),]
x = x[,6]

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

hist(x, main="")
rug(x)

library(diptest)
dtb = dip(x, full.result = TRUE, min.is.0 = TRUE, debug = FALSE); dtb
dip.test(x, simulate.p.value = FALSE, B = 2000)
plot(dtb)

library(Rfolding)
folding.test(as.matrix(x))


#### Kernel density estimate ########################################################################################

library(multimode)
locmodes(x,mod0=1,display=TRUE,xlab="x")
locmodes(x,mod0=2,display=TRUE,xlab="x")
locmodes(x,mod0=3,display=TRUE,xlab="x")
locmodes(x,mod0=4,display=TRUE,xlab="x")

modetree(x,bws=NULL,gridsize=NULL,cbw1=NULL,cbw2=NULL,display=TRUE,logbw=FALSE,logbw.regulargrid=NULL,xlab="x")
modetree(x,bws=c(0.1,1),gridsize=NULL,cbw1=NULL,cbw2=NULL,display=TRUE,logbw=FALSE,logbw.regulargrid=NULL,xlab="x")

bw1 = bw.crit(x,mod0=1,lowsup=-Inf,uppsup=Inf,n=2^15,tol=10^(-5),full.result=F)
bw2 = bw.crit(x,mod0=2,lowsup=-Inf,uppsup=Inf,n=2^15,tol=10^(-5),full.result=F)
bw3 = bw.crit(x,mod0=3,lowsup=-Inf,uppsup=Inf,n=2^15,tol=10^(-5),full.result=F)
bw4 = bw.crit(x,mod0=4,lowsup=-Inf,uppsup=Inf,n=2^15,tol=10^(-5),full.result=F)


library(silvermantest) # https://www.mathematik.uni-marburg.de/~stochastik/R_packages/
modes=1:4
densities.plot(x,modes=modes,mark_modes=TRUE,in_one=TRUE*(length(modes)<4))
#  null hypothesis: underlying density has at most k modes
silverman.test(x, k=1, M = 999, adjust = FALSE, digits = 6)
silverman.test(x, k=2, M = 999, adjust = FALSE, digits = 6)
silverman.test(x, k=3, M = 999, adjust = FALSE, digits = 6)
silverman.plot(x, kmin=1, kmax=5, alpha=0.05, adjust=FALSE)


library(mixtools)
fit <- normalmixEM(x, k = 2) #try to fit k Gaussians
summary(fit)



###################################################################################################################
######### Testy shody ###############################################################################################
##########################################################################################################

b314<- read.table("c:\\Users\\lubop\\Documents\\kurs\\Data1\\ascii03\\B\\B314.txt", header=T)
# dekl obsah 2.4 mg/l
is.list(b314)
x = unlist(b314)

#### t-test ############################################################################################
# one sample t-test: t.test(x,mu=3, alternative = "two.sided") 

t.test(x, y = NULL, alternative = "two.sided", mu = 2.4, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
# alternative = c("two.sided", "less", "greater")

library(pwr)
p.t.one <- pwr.t.test(d=(abs(mean(x)-2.4)), power=0.9, type="one.sample", alternative="two.sided")
plot(p.t.one)


#### medianovy test ############################################################################################

library(BSDA)
SIGN.test(x, y = NULL, md = 2.4, alternative = "two.sided", conf.level = 0.95)


#### Wilcoxonuv (Mann-Whitney) test ############################################################################################

wilcox.test(x, y = NULL, alternative = "two.sided", mu = 2.4, paired = FALSE, exact = NULL, correct = TRUE, conf.int = TRUE, conf.level = 0.95) 

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


##### Dvouvyberove testy  ###########################################################################################################################################


###### Scale #######################################################################################

### F test
var.test(x, y, ratio = 1, alternative = "two.sided", conf.level = 0.95)
# alternative = c("two.sided", "less", "greater")

### Ansari Bradley two sample test
ansari.test(x, y, alternative = "two.sided", exact = FALSE, conf.int = FALSE, conf.level = 0.95)


###### Location ###################################################################################

b325ab<- read.table("c:\\Users\\lubop\\Documents\\kurs\\Data1\\ascii03\\B\\B325ab.txt", header=T)
# dekl obsah 2.4 mg/l
is.list(b325ab)
X = as.data.frame(b325ab)
# x = unlist(b325ab)
x= unlist(X[,1])
y= unlist(X[,2])


#### t-test ############################################################################################

# t.test(x,mu=3, alternative = "two.sided") 

sd(x)
sd(y)
# srovnani pomoci F-testu nebo jineho testu shody rozptylu

boxplot(x,y)

t.test(x, y, alternative = "two.sided", mu = 0, paired = FALSE, var.equal = TRUE, conf.level = 0.95)
t.test(x, y, alternative = "two.sided", mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) # Welchuv test
# alternative = c("two.sided", "less", "greater")

library(pwr)
p.t.two <- pwr.t.test(d=(abs(mean(x)-mean(y))), power=0.9, type="two.sample", alternative="two.sided")
plot(p.t.two)
plot(p.t.two, xlab="sample size per group")

p.t.two2 <- pwr.t2n.test(n1 = length(x), n2= length(y), d = (abs(mean(x)-mean(y))), sig.level = 0.05, power = NULL,
                         alternative ="two.sided")
# alternative = c("two.sided","less","greater")
plot(p.t.two2)
plot(p.t.two2, xlab="sample size per group")


#### 2 sample medianovy test ############################################################################################

# one sample sign test 
library(BSDA)
SIGN.test(x, y, md = 0, alternative = "two.sided", conf.level = 0.95)


#### 2 sample Wilcoxonuv (Mann-Whitney) test ############################################################################################


wilcox.test(x, y,alternative = "two.sided", mu = 0, paired = FALSE, exact = FALSE, correct = TRUE,conf.int = FALSE, conf.level = 0.95)
# alternative = c("two.sided", "less", "greater")

library(exactRankTests)
wilcox.exact(x, y, mu = 0,alternative = "two.sided",paired = FALSE, exact = FALSE,conf.int = TRUE, conf.level = 0.95)

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

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

### Srovnani 2 distribuci ##########################################################################################


### Kolmogorov Smirnov test
ks.test(x,y, alternative = "two.sided", exact = FALSE)
# alternative = c("two.sided", "less", "greater")

library(KScorrect)
ks_test_stat(x, y)

library(dgof)
ks.test(x, y,alternative = "two.sided", exact = NULL, tol=1e-8, simulate.p.value=TRUE, B=2000)
# alternative = c("two.sided", "less", "greater")

library(twosamples)
ks_test(x, y, nboots = 2000, p = 0.025)


### Cramer von Mises
library(CDFt)
CramerVonMisesTwoSamples(x, y)

library(cramer)
cramer.test(x,y,conf.level=0.95,replicates=1000,sim="ordinary",just.statistic=FALSE,kernel="phiCramer", maxM=2^14, K=160)



###################################################################################################################
######### Testy shody ###############################################################################################
##########################################################################################################

b314<- read.table("c:\\Users\\lubop\\Documents\\kurs\\Data1\\ascii03\\B\\B314.txt", header=T)
# dekl obsah 2.4 mg/l
is.list(b314)
x = unlist(b314)


#### t-test ############################################################################################
# one sample t-test: t.test(x,mu=3, alternative = "two.sided") 

t.test(x, y = NULL, alternative = "two.sided", mu = 2.4, paired = FALSE, var.equal = FALSE, conf.level = 0.95)
# alternative = c("two.sided", "less", "greater")

library(pwr)
p.t.one <- pwr.t.test(d=(abs(mean(x)-2.4)), power=0.9, type="one.sample", alternative="two.sided")
plot(p.t.one)


#### medianovy test ############################################################################################

library(BSDA)
SIGN.test(x, y = NULL, md = 2.4, alternative = "two.sided", conf.level = 0.95)


#### Wilcoxonuv (Mann-Whitney) test ############################################################################################

wilcox.test(x, y = NULL, alternative = "two.sided", mu = 2.4, paired = FALSE, exact = NULL, correct = TRUE, conf.int = TRUE, conf.level = 0.95) 

##### Dvouvyberove testy  ###########################################################################################################################################


###### Scale #######################################################################################

### F test
var.test(x, y, ratio = 1, alternative = "two.sided", conf.level = 0.95)
# alternative = c("two.sided", "less", "greater")

### Ansari Bradley two sample test
ansari.test(x, y, alternative = "two.sided", exact = FALSE, conf.int = FALSE, conf.level = 0.95)


###### Location ###################################################################################

b325ab<- read.table("c:\\Users\\lubop\\Documents\\kurs\\Data1\\ascii03\\B\\B325ab.txt", header=T)
# dekl obsah 2.4 mg/l
is.list(b325ab)
X = as.data.frame(b325ab)
# x = unlist(b325ab)
x= unlist(X[,1])
y= unlist(X[,2])


#### t-test ############################################################################################

# t.test(x,mu=3, alternative = "two.sided") 

sd(x)
sd(y)
# srovnani pomoci F-testu nebo jineho testu shody rozptylu

boxplot(x,y)

t.test(x, y, alternative = "two.sided", mu = 0, paired = FALSE, var.equal = TRUE, conf.level = 0.95)
t.test(x, y, alternative = "two.sided", mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95) # Welchuv test
# alternative = c("two.sided", "less", "greater")

library(pwr)
p.t.two <- pwr.t.test(d=(abs(mean(x)-mean(y))), power=0.9, type="two.sample", alternative="two.sided")
plot(p.t.two)
plot(p.t.two, xlab="sample size per group")

p.t.two2 <- pwr.t2n.test(n1 = length(x), n2= length(y), d = (abs(mean(x)-mean(y))), sig.level = 0.05, power = NULL,
                         alternative ="two.sided")
# alternative = c("two.sided","less","greater")
plot(p.t.two2)
plot(p.t.two2, xlab="sample size per group")


#### 2 sample medianovy test ############################################################################################

# one sample sign test 
library(BSDA)
SIGN.test(x, y, md = 0, alternative = "two.sided", conf.level = 0.95)


#### 2 sample Wilcoxonuv (Mann-Whitney) test ############################################################################################


wilcox.test(x, y,alternative = "two.sided", mu = 0, paired = FALSE, exact = FALSE, correct = TRUE,conf.int = FALSE, conf.level = 0.95)
# alternative = c("two.sided", "less", "greater")

library(exactRankTests)
wilcox.exact(x, y, mu = 0,alternative = "two.sided",paired = FALSE, exact = FALSE,conf.int = TRUE, conf.level = 0.95)


### Srovnani 2 distribuci ##########################################################################################


### Kolmogorov Smirnov test
ks.test(x,y, alternative = "two.sided", exact = FALSE)
# alternative = c("two.sided", "less", "greater")

library(KScorrect)
ks_test_stat(x, y)

library(dgof)
ks.test(x, y,alternative = "two.sided", exact = NULL, tol=1e-8, simulate.p.value=TRUE, B=2000)
# alternative = c("two.sided", "less", "greater")

library(twosamples)
ks_test(x, y, nboots = 2000, p = 0.025)


### Cramer von Mises
library(CDFt)
CramerVonMisesTwoSamples(x, y)

library(cramer)
cramer.test(x,y,conf.level=0.95,replicates=1000,sim="ordinary",just.statistic=FALSE,kernel="phiCramer", maxM=2^14, K=160)



lxy = list(x,y); names(lxy)= c("x","y")
data = stack(lxy)
data = as.data.frame(data); colnames(data)=c("xy","z")
xy = data[,1]
z = data[,2]
stripchart(xy~z, method="jitter",pch=1,col=1,cex=1.5) 
boxplot(xy~z, pch=1,col=5,cex=1.5) 
stripchart(xy~z, method="jitter",pch=1,col=1,cex=1.5) 

### Anderson Darling
library(kSamples)
ad.test(xy~z, data = data, method = "simulated", dist = FALSE, Nsim = 10000)
# method = c("asymptotic", "simulated", "exact")


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


E503 <- read.table("C:/Users/lubop/Documents/kurs/Data1/ascii05/E/E503.txt")[-1,]
E503
data = stack(E503)
data = as.data.frame(data); colnames(data)=c("elem","lom")
elem = data[,1]
lom = data[,2]


### Distribucni grafy

boxplot(elem~lom,ylab="conc",xlab="lab",las=1,cex=1,outline=TRUE) 
stripchart(elem~as.factor(lom), pch=1, col=1,vert=TRUE,cex=1,add=TRUE)

library(beeswarm)
boxplot(elem~lom,ylab="",xlab="",las=1,cex=1,outline=TRUE,cex.lab=1) 
beeswarm(elem~as.factor(lom), col = 1, pch = 1, add = TRUE,horizontal=FALSE,pwcol = NULL,cex=1)

library(beanplot)
beanplot(elem ~ lom, data = data, col = "bisque", ylim = "")



#### Global test of model assumptions
library(gvlma)
gvmodel <- gvlma(lm(elem~lom)) 
summary(gvmodel)



# Bonferroni Outlier Test
library(car)
ou = outlierTest(lm(elem~lom), cutoff=0.05, n.max=15, order=FALSE); ou
as.numeric(names(ou$rstudent))

# Lund Outlier Test
# modified https://stackoverflow.com/questions/1444306/how-to-use-outlier-tests-in-r-code/1444548#1444548
lundout<-function(x1,x2,a) {
  testoutlm<-lm(x1~x2)
  standardresid<-rstandard(testoutlm)
  nl<-length(x1)
  q<-length(testoutlm$coefficients)
  F<-qf(c(1-(a/nl)),df1=1,df2=nl-q-1,lower.tail=TRUE)
  crit<-((nl-q)*F/(nl-q-1+F))^0.5
  oul = which(abs(standardresid)>crit)
  return(oul)
}
lundout(elem,lom,0.05)

### outliers
stresid<-rstandard(lm(elem~lom))
library(referenceIntervals)
vanderLoo.outliers(stresid)
cook.outliers(stresid)
dixon.outliers(stresid)
horn.outliers(stresid)


#### shoda rozptylu 

## Bartlett
bartlett.test(elem~as.factor(lom))
library(RVAideMemoire)
pairwise.var.test(elem,as.factor(lom))

library(RVAideMemoire)
perm.bartlett.test(elem~lom, nperm = 999, progress = TRUE)
pairwise.perm.var.test(elem,as.factor(lom),nperm=999)


## Fligner Kileen
fligner.test(elem~lom)

library(coin)
fligner_test(elem~lom, data=data, subset = NULL, weights = NULL,  conf.int = TRUE, conf.level = 0.95, ties.method = "average-scores")
fligner_test(elem~lom, data=data, subset = NULL, weights = NULL,  conf.int = TRUE, conf.level = 0.95, ties.method = "average-scores", distribution = approximate(nresample = 1000))
# ties.method = c("mid-ranks", "average-scores"),


## Conover, Johnson, Johnson
library(coin)
conover_test(elem~lom, data=data, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95,ties.method = "average-scores")
conover_test(elem~lom, data=data, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95,ties.method = "average-scores", distribution = approximate(nresample = 1000))


## Levene
library(lawstat)
lawstat::levene.test(elem, group = as.factor(lom), location="mean", trim.alpha=0.25,bootstrap = TRUE, num.bootstrap=999, kruskal.test=TRUE,correction.method="none")
# location=c("median", "mean", "trim.mean")
# correction.method=c("none","correction.factor","zero.removal","zero.correction")
library(car)
car::leveneTest(elem~as.factor(lom),center=mean)
car::leveneTest(elem~as.factor(lom),center=median)
# center = c("median", "mean")


#### ANOVA homoskedastic

fit <- aov(elem~lom)
summary(fit) # display Type I ANOVA table
summary.lm(fit)
anova(fit)
model.tables(fit,se=T)
model.tables(fit,"means",se=T)
library(car)
Anova(fit, type="II")
### If you use type="III", you need the following line before the analysis
### options(contrasts = c("contr.sum", "contr.poly"))



library(FSA)   
Summarize(elem ~ lom, data=data, digits=3)

library(multcompView)
multcompBoxplot(elem ~ lom, data = data, horizontal = FALSE, compFn = "TukeyHSD",sortFn = "mean", decreasing = TRUE, plotList = list(boxplot = list(fig = c(0, 0.75, 0, 1)),                                                                                                                                multcompTs = list(fig = c(0.7, 0.85, 0, 1)),multcompLetters = list(fig = c(0.87, 0.97, 0.03, 0.98), fontsize = 20,fontface = "bold")))


## Scheffe test
library(agricolae)
out <- scheffe.test(aov(elem~lom),"lom",alpha = 0.95, group=TRUE,console=TRUE,main="")
plot(out,las=2)

library(PMCMRplus)
out = scheffeTest(elem~lom, subset=NULL, na.action=NULL)
summary(out)
summaryGroup(out)

# least significant difference all-pairs comparisons test for normally distributed data with equal group variances.
library(PMCMRplus)
out = lsdTest(elem,lom, subset=NULL, na.action=NULL); out
summary(out)
summaryGroup(out)

library(agricolae)
out <- LSD.test(aov(elem~lom),"lom", alpha=0.95, group=TRUE,main=NULL,console=TRUE,p.adj="bonferroni")
# p.adj=c("none","holm","hommel","hochberg", "bonferroni", "BH", "BY", "fdr")
plot(out,las=2)


#### ANOVA heteroskedastic


## Alexander-Govern test.
library(onewaytests)
out = ag.test(elem ~ lom, data = data, alpha = 0.05, na.rm = TRUE, verbose = TRUE)
out
out =paircomp(out, adjust.method = "fdr")
# p.adjust.methods: "holm", "hochberg","hommel","bonferroni","BH","BY","fdr","none","tukey","games-howell"


## Welch test

library(car)
mod = aov(elem~lom, data=data, na.action=na.omit)
# mod = lm(elem~lom, data=as.data.frame(cbind(as.numeric(elem),lom)), na.action=na.omit)
Anova(mod, white.adjust=TRUE)
# white.adjust=c(FALSE, TRUE, "hc3", "hc0", "hc1", "hc2", "hc4")

oneway.test(elem~as.factor(lom), subset=NULL, na.action=NULL, var.equal = FALSE)

library(welchADF)
summary(welchADF.test(elem~lom))


library(asbio)
pairw.oneway(elem,as.factor(lom), conf = 0.95, digits = 5, method = "fdr")

library(userfriendlyscience)
oneway(elem,lom, fullDescribe = TRUE, posthoc = 'fdr',corrections=TRUE,plot=TRUE)
# p.adjust.methods: "holm", "hochberg","hommel","bonferroni","BH","BY","fdr","none","tukey","games-howell"

library(coin)
oneway_test(elem~lom, data=data, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95, ties.method = "average-scores")
oneway_test(elem~lom, data=data, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95, ties.method = "average-scores", distribution = approximate(nresample = 1000))


## Games-Howell
library(userfriendlyscience)
posthocTGH(elem,lom, method= "games-howell", conf.level = 0.95, digits=2, p.adjust="holm",formatPvalue = TRUE)
# method=c("games-howell", "tukey")
# p.adjust= c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none")

library(PMCMRplus)
out = gamesHowellTest(elem~lom,alternative = "two.sided",subset=NULL, p.adjust="BH")
summary(out)
summaryGroup(out)


## Tamhane T2
library(PMCMRplus)
out = tamhaneT2Test(elem,lom, subset=NULL, na.action=NULL,welch = TRUE,alternative = "two.sided")
summary(out)
summaryGroup(out)


#### shoda distribuce #########################################################################################

library(PMCMRplus)
adKSampleTest(elem~lom, subset=NULL, na.action=NULL)
out = adAllPairsTest(elem~lom, subset=NULL, na.action=NULL, p.adjust.method = "BH"); out
# p.adjust.methods = c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none")
summary(out)
summaryGroup(out)
barPlot(out, alpha = 0.05)
plot(out, alpha = 0.05)

library(PMCMRplus) # Murakami k-sample Baumgartner Weiss Schindler Test of Equal Distributions
bwsKSampleTest(elem,lom, subset=NULL, na.action,nperm = 1000)
out <- bwsAllPairsTest(elem, lom, p.adjust="holm"); out
summary(out)
summaryGroup(out)


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

#### Non-parametic

## Kruskal-Wallis test  # homoscedastic
kruskal.test(elem~lom)

library(agricolae)
kruskal(elem,lom, alpha = 0.05, p.adj="holm",group=TRUE, main = NULL,console=TRUE)
# p.adj=c("none","holm","hommel", "hochberg", "bonferroni", "BH", "BY", "fdr")

library(asbio)
pairw.kw(elem, lom, conf=0.05)
plot(pairw.kw(elem, lom, conf=0.05),las=2)


library(NSM3)
pKW(elem,lom,method="Monte Carlo",n.mc=500) # Kruskal-Wallis
# method=c("Exact", "Monte Carlo", or "Asymptotic")

library(coin)
kruskal_test(elem~lom, data=data, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95, ties.method = "average-scores") 
kruskal_test(elem~lom, data=data, subset = NULL, weights = NULL, conf.int = FALSE, conf.level = 0.95, ties.method = "average-scores", distribution = approximate(nresample = 1000)) 

library(onewaytests)
out = kw.test(elem ~ lom, data = data, alpha = 0.05, na.rm = TRUE, verbose = TRUE)
# paircomp(out, adjust.method = "bonferroni")
# adjust.method = c("bonferroni", "holm", "hochberg", "hommel", "BH","BY", "fdr", "none")

library(PMCMRplus)
kruskalTest(elem, lom, subset=NULL, na.action=NULL, dist = "KruskalWallis")
# dist = c("Chisquare", "KruskalWallis", "FDist")


## Nemenyi
library(PMCMR)
posthoc.kruskal.nemenyi.test(elem,lom, method="Tukey")
# method="Tukey", "Chisq"

library(PMCMRplus)
out = kwAllPairsNemenyiTest(elem, lom, subset=NULL, na.action=NULL,dist = "Tukey"); out
# dist = c("Tukey", "Chisquare")
summary(out)
summaryGroup(out)

library(PMCMRplus)
out=kwAllPairsDunnTest(elem, lom, subset=NULL, na.action=NULL, p.adjust.method ="fdr"); out
summary(out)
summaryGroup(out)

library(PMCMR)
out = posthoc.kruskal.dunn.test(elem,lom,p.adjust.method ="holm"); out
summary(out)


## Conover's non-parametric all-pairs comparison test for Kruskal-type ranked data.
library(conover.test)
conover.test(elem,lom, method="holm", kw=TRUE, label=TRUE,wrap=FALSE, table=TRUE, list=FALSE, rmc=FALSE, alpha=0.95)
# method=c("none", "bonferroni", "sidak", "holm", "hs", "hochberg", "bh", "by")

library(PMCMRplus)
out = kwAllPairsConoverTest(elem, lom, subset=NULL, na.action=NULL,p.adjust.method = "holm"); out
# p.adjust.method = c("single-step", p.adjust.methods); single-step = Tukey's p-adjustment
summary(out)
summaryGroup(out)


library(PMCMR)
posthoc.kruskal.conover.test(elem, lom,p.adjust.method ="holm")


## Dwass, Steel, Critchlow and Fligner
library(PMCMRplus)
out = dscfAllPairsTest(elem, lom, subset=NULL, na.action=NULL, p.adjust.method = "holm", alternative = "two.sided")
summary(out)
summaryGroup(out)


E518 <- read.table("C:/Users/lubop/Documents/kurs/Data1/ascii05/E/E518.txt")
data = stack(E518[-1,])
data = as.data.frame(data); colnames(data)=c("elem","lom")
elem = data[,1]
lom = data[,2]


# Jonckheere-Terpstra (testuje trend v poradi souboru)
library(DescTools)
lomc =  ordered(sapply(lom,function(x) sub("V", "", x)))   
JonckheereTerpstraTest(elem, lomc, alternative = "increasing",nperm = 500)
# alternative = c("two.sided", "increasing", "decreasing")


#### normal scores

## van-der-Waerden normal scores test
library(agricolae)
out<-waerden.test(elem,lom,alpha=0.95,group=TRUE,console=TRUE,main="")
plot(out,las=2)

library(PMCMRplus)
vanWaerdenTest(elem,as.factor(lom), subset=NULL, na.action=NULL,alternative = "two.sided")
out = vanWaerdenAllPairsTest(elem,as.factor(lom), subset=NULL, na.action=NULL,alternative = "two.sided",p.adjust.method = "fdr")
#  p.adjust.method = c("single-step", p.adjust.methods)
summary(out)
summaryGroup(out)

library(PMCMR)
vanWaerden.test(elem,as.factor(lom), subset=NULL, na.action=NULL,alternative = "two.sided")
out = posthoc.vanWaerden.test(elem,as.factor(lom),p.adjust.method ="holm")
summary(out)

library(coin)
normal_test(as.numeric(elem)~as.factor(lom), ties.method = "mid-ranks", conf.int = FALSE, conf.level = 0.95)
normal_test(as.numeric(elem)~as.factor(lom), ties.method = "mid-ranks", conf.int = FALSE, conf.level = 0.95, distribution = approximate(B = 1000))
# method = c("mid-ranks", "average-scores")


# Lu-Smith normal scores test, transformace dat na normalitu
library(PMCMRplus)
normalScoresTest(elem,lom, subset=NULL, na.action=NULL,alternative = "two.sided")

out = normalScoresAllPairsTest(elem,lom, subset=NULL, na.action=NULL,alternative = "two.sided", p.adjust.method = "holm"); out
# p.adjust.method = c("single-step", p.adjust.methods)
summary(out)
summaryGroup(out)



## Median test
library(agricolae)
out<-Median.test(elem,lom,correct=TRUE,simulate.p.value = TRUE,alpha=0.95,group=TRUE,console=TRUE,main="")
plot(out,las=2)

library(coin)
median_test(elem~lom, mid.score = c("0", "0.5", "1"), conf.int = FALSE, conf.level = 0.95)
median_test(elem~lom, mid.score = c("0", "0.5", "1"), conf.int = FALSE, conf.level = 0.95, distribution = approximate(nresample = 1000))

library(RVAideMemoire)
mood.medtest(elem, lom, exact = TRUE)
pairwise.mood.medtest(elem, lom, exact = TRUE, p.method = "fdr")



##### 2-way ANOVA ####################################################################################################################################

library(openxlsx)
shn = getSheetNames("c:\\Users\\lubop\\Documents\\kurs\\DVpisky1.xlsx")
shn
sal <- read.xlsx("c:\\Users\\lubop\\Documents\\kurs\\DVpisky1.xlsx", sheet = "2wANOVA", startRow = 1, colNames = TRUE, rowNames = FALSE,detectDates = FALSE, skipEmptyRows = TRUE,rows = NULL, cols = NULL,check.names = FALSE)
data = as.data.frame(sal)
nitr = data[,1]
subst = data[,2]
hnoj = data[,3]

### Global test of model assumptions
library(gvlma)
gvmodel <- with(data,gvlma(lm(nitr~subst*hnoj)))
summary(gvmodel)

fit = lm(nitr~subst*hnoj,data=data)
fit = lm(nitr~subst*hnoj-1,data=data)
summary(fit)
anova(fit)

fit2 = aov(nitr~subst*hnoj,data=data)
anova(fit2)
summary(fit2)



library(Hmisc)
with(data,summary(nitr~subst*hnoj))
with(data,summary(nitr~subst*hnoj-1))

library(car)
(my_anova <- aov(nitr~subst*hnoj,data=data))
anova(my_anova)
Anova(my_anova, type = "II") # preferable
# Anova(my_anova, type = "III")

with(data,interaction.plot(subst,hnoj,nitr,ylab="",type ="l"))

summary(aov(nitr~subst*hnoj,data=data))
with(data,interaction.plot(subst,hnoj,nitr,ylab=""))
anova(aov(nitr~subst*hnoj,data=data))

library(ExpDes)
with(data, fat2.crd(subst, hnoj, nitr, quali = c(TRUE, TRUE), mcomp = "tukey",
                    fac.names = c("substrat", "hnojeni"), sigT = 0.05, sigF = 0.05))

library(welchADF)
summary(welchADF.test(lm(nitr~subst*hnoj,data=data)))


library(lmPerm)
# perm "Exact", "Prob", "SPR" will produce permutation probabilities. Anything else, such as "", will produce F-test probabilities.
out = aovp(nitr~subst*hnoj,data=data,perm="Prob",nCycle = 5000); out
summary(out)
anova(lmp(nitr~subst*hnoj,data=data))

library(permuco)
out = aovperm(nitr~subst*hnoj,data=data, np = 5000, method = "freedman_lane"); out
summary(out)

library(robustbase)
model = lmrob(nitr~subst*hnoj-1,data=data) 
summary(model)
plot(fitted(model), residuals(model))
abline(h=0, col=2, lty=2)


library(robustbase)
model = lmrob(nitr~subst*hnoj,data=data) 
summary(model)
plot(fitted(model), residuals(model))
abline(h=0, col=2, lty=2)


###### Regrese #####################################################################

library(investr)
data(arsenic)

library(ACSWR)
data(viscos)
data=viscos
x = data[,1]
y = data[,2]

library(openxlsx)
sal <- read.xlsx("c:\\Users\\lubop\\Documents\\kurs\\DVpisky1.xlsx", sheet = "regrese2", startRow = 1, colNames = TRUE, rowNames = FALSE,detectDates = FALSE, skipEmptyRows = TRUE,rows = NULL, cols = NULL,check.names = FALSE)
data = as.data.frame(sal)
x = data[,1]
y = data[,2]
z = data[,3]

plot(x,y)


lmMod <- lm(y~x, data=data, weights=NULL, na.action = na.pass)
summary(lmMod)

plot(x,y)
abline(lm(y~x, data=data), col=2, lty=2)
abline(lm(y~x-1, data=data), col=4, lty=2)

cp = predict(lmMod, newdata = 0, interval = 'confidence')
lines(x,cp[,2],col=6,lty=3)
lines(x,cp[,3],col=6,lty=3)

pp = predict(lmMod, newdata = 0, interval = 'prediction')
lines(x,pp[,2],col=3,lty=3)
lines(x,pp[,3],col=3,lty=3)

confint(lmMod,level = 0.95)
confint(lmMod,level = 0.99)

lm.out <- lm(y ~ x)
newx = seq(min(x),max(x),by = 0.05)

conf_interval <- predict(lm.out, newdata=data.frame(x=newx), interval="confidence", level = 0.95)
plot(x, y, xlab="x", ylab="y", main=" ")
abline(lm.out, col="lightblue")
lines(newx, conf_interval[,2], col="blue", lty=2)
lines(newx, conf_interval[,3], col="blue", lty=2)

pred_interval <- predict(lm.out, newdata=data.frame(x=newx), interval="prediction", level = 0.95)
lines(newx, pred_interval[,2], col="green", lty=2)
lines(newx, pred_interval[,3], col="green", lty=2)


# test koeficientu
library(lmtest)
coeftest(lmMod, vcov. = NULL, df = NULL)
coefci(lmMod, parm = NULL, level = 0.95, vcov. = NULL, df = NULL)

library(car)
confidenceEllipse(lmMod)


# fitted
fitted.values(lmMod)

library(visreg)
visreg(lmMod, "x")


# residuals
resid(lmMod) # klasicka
plot(fitted.values(lmMod), resid(lmMod)); abline(h=0, col=2, lty=2)
lines(lowess(fitted.values(lmMod), resid(lmMod)), col = 4)

rstandard(lmMod) # standardizovana
plot(fitted.values(lmMod), rstandard(lmMod)); abline(h=0, col=2, lty=2)
lines(lowess(fitted.values(lmMod), rstandard(lmMod)), col = 4)

rstudent(lmMod) # studentizovana
plot(fitted.values(lmMod), rstudent(lmMod)); abline(h=0, col=2, lty=2)
lines(lowess(fitted.values(lmMod), rstudent(lmMod)), col = 4)

library(s20x)
ciReg(lmMod, conf.level = 0.95, print.out = TRUE)

cooks20x(lmMod, main = "Cook's Distance plot", xlab = "observation number",ylab = "Cook's distance", line = c(0.5, 1.2, 2), cex.labels = 1,axisOpts = list(xAxis = TRUE, yAxisTight = FALSE))

eovcheck(lmMod, xlab = "Fitted values",ylab = "Residuals", col = NULL, smoother = FALSE, twosd = FALSE,levene = FALSE)

modcheck(lmMod, plotOrder = 1:4, args = list(eovcheck = list(smoother = FALSE, twosd = FALSE, levene = FALSE), normcheck = list(xlab = c("Theoretical Quantiles", ""), ylab = c("Sample Quantiles",""), main = c("", ""), col = "light blue", bootstrap = FALSE, B = 5, bpch = 3, bcol = "lightgrey", shapiro.wilk = FALSE, whichPlot = 1:2, usePar =   TRUE), cooks20x = list(main = "Cook's Distance plot", xlab = "observation number", ylab = "Cook's distance", line = c(0.5, 0.1, 2), cex.labels = 1, axisOpts = list(xAxis = TRUE))), parVals = list(mfrow = c(2, 2), xaxs = "r", yaxs = "r", pty = "s", mai= c(0.2, 0.2, 0.05, 0.05)))

normcheck(lmMod, xlab = c("Theoretical Quantiles", ""),ylab = c("Sample Quantiles", ""), main = c("", ""), col = "light blue", bootstrap = FALSE, B = 5, bpch = 3, bcol = "lightgrey", shapiro.wilk = FALSE, whichPlot = 1:2, usePar = TRUE)

residPlot(lmMod, f = 1)


library(car)
residualPlots(lmMod, terms = ~., layout = NULL, main = "", fitted = TRUE, AsIs=TRUE, plot = TRUE,tests = TRUE)
influencePlot(lmMod)
infIndexPlot(lmMod)


library(UsingR)
simple.lm(x, y, show.residuals=TRUE, show.ci=TRUE, conf.level=0.95,pred=NULL) 

par(mfrow=c(2,2))
plot(lmMod)
par(mfrow=c(1,1))

library(lindia)
gg_boxcox(lmMod, showlambda = TRUE, lambdaSF = 3, scale.factor = 1)
gg_cooksd(lmMod, label = TRUE, show.threshold = TRUE,threshold = "convention", scale.factor = 1)
gg_diagnose(lmMod, theme = NULL, ncol = NA, plot.all = TRUE,scale.factor = 1, boxcox = FALSE, max.per.page = NA)
gg_qqplot(lmMod, scale.factor = 1)
gg_resfitted(lmMod, scale.factor = 1)
gg_reshist(lmMod, bins = 20)
gg_resleverage(lmMod, method = "loess", se = FALSE, scale.factor = 1)
gg_resX(lmMod, plot.all = TRUE, scale.factor = 1, max.per.page = NA, ncol = NA)
gg_scalelocation(lmMod, method = "loess", scale.factor = 1, se = FALSE)


# Bonferroni Outlier Test
library(car)
ou = outlierTest(lm(y~x), cutoff=0.05, n.max=15, order=TRUE); ou
as.numeric(names(ou$rstudent))

# Lund Outlier Test
lundout<-function(x1,x2,a) {
  testoutlm<-lm(x1~x2)
  standardresid<-rstandard(testoutlm)
  nl<-length(x1)
  q<-length(testoutlm$coefficients)
  F<-qf(c(1-(a/nl)),df1=1,df2=nl-q-1,lower.tail=TRUE)
  crit<-((nl-q)*F/(nl-q-1+F))^0.5
  oul = which(abs(standardresid)>crit)
  return(oul)
}
lundout(y,x,0.05)

# Shapiro-Wilk test rezidui
shapiro.test(resid(lmMod))
shapiro.test(rstandard(lmMod))
shapiro.test(rstudent(lmMod))


# Breusch-Pagan test na heteroskedasticitu
library(lmtest)
lmtest::bptest(y~x, varformula = NULL, studentize = TRUE, data = data) 

# Harrison-McCabe test for heteroskedasticity
library(lmtest)
lmtest::hmctest(y~x, point = 0.5, order.by = NULL, simulate.p = TRUE, nsim = 1000, plot = TRUE, data = data)

# Score Test for Non-Constant Error Variance
library(car)
car::ncvTest(lmMod)

# Goldfeld-Quandt test against heteroskedasticity
library(lmtest)
lmtest::gqtest(y~x, point = 0.5, fraction = 0, alternative = "two.sided",order.by = NULL, data = data)
# alternative = c("greater", "two.sided", "less")


# Durbin-Watson test for autocorrelation (zavislost rezidui na x)
library(lmtest)
lmtest::dwtest(y~x, order.by = NULL, alternative = "two.sided", iterations = 15, exact = NULL, tol = 1e-10, data = data)
# alternative = c("greater", "two.sided", "less")
library(car)
durbinWatsonTest(lmMod)

### Testy linearity
# Rainbow test linearity
library(lmtest)
lmtest::raintest(y~x, fraction = 0.5, order.by = NULL, center = NULL, data=data)

# Harvey-Collier test for linearity
library(lmtest)
lmtest::harvtest(y~x, order.by = NULL, data = data)



###### regrese pocatkem ############################

s.lm <- lm(y ~ x, data = data)
summary(s.lm)

s.lm2 <- lm(y ~ 0 + x, data = data) # Adding the 0 term tells the lm() to fit the line through the origin
summary(s.lm2)

s.lm3 <- lm(y ~ x-1, data = data) # Adding the 0 term tells the lm() to fit the line through the origin
summary(s.lm3)

plot(x,y) # regrese bez useku
abline(lm(y ~ x-1, data = data), col=2, lty=2)


### Box-Cox #############################################################
library(AID)
lmbc = boxcoxlm(as.matrix(x), y, method = "lse", lambda = seq(-3,3,0.01), lambda2 = NULL, plot = TRUE,alpha = 0.05, verbose = TRUE)
lmbc$lambda.hat
caret::BoxCoxTrans(y)
library(forecast)
ynew = BoxCox(y, lambda=2)
lmMod_bc <- lm(ynew ~ x, data=data)
bptest(lmMod_bc)
plot(lmMod_bc)


### regrese polynomem
s.lmq <- lm(y ~ x + poly(x,2), data = data)
summary(s.lmq)

plot(x,y) # regrese polynomem
lines(lm(y ~ poly(x^2), data = data), col=2, lty=2)
abline(lm(y ~ x + poly(x,2), data = data), col=4, lty=2)

lmMod2 = lm(z ~ x, data = data)

### Srovnani 2 modelu
anova(lmMod,lmMod2)

library(performance)
compare_performance(lmMod, lmMod2, rank = TRUE)

# Srovnani 2 primek
library(gap)
chow.test(y1=y,x1=x,y2=z,x2=x,x=NULL)


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

### Resistant Regression
library(MASS)
lmM1 = lqs(y ~ x, data,method = "lts")
lmM1
# method = c("lts", "lqs", "lms", "S", "model.frame")

lmM2 = rlm(y ~ x, data, weights=NULL, method = "M", wt.method = "inv.var")
lmM2 
# method = c("M", "MM", "model.frame"),
# wt.method = c("inv.var", "case")

library(LearnEDA) # Fits Tukey resistant line of form a + b (x - xC).
lmM3 = rline(x,y,iter=1)


### regrese s LOD
library(lodr)
x2 = 
  fitl <- lod_lm(data=lod_data, frmla=y~x, lod=0.1,var_LOD="x2", nSamples=100, boots=500)
summary(fitl)
coef(fit)
residuals(fitl)


library(truncreg)
truncreg(y~x, data=data, point = 2, direction = "left",model = TRUE, y = FALSE, x = FALSE, scaled = FALSE)


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

library(chemCal)
data(din32645)
din32645
data(massart97ex1)
massart97ex1
data(massart97ex3)
massart97ex3
data(utstats14)


m <- lm(y ~ x, data = din32645)
calplot(m,xlim = c("auto", "auto"), ylim = c("auto", "auto"),xlab = "Concentration", ylab = "Response", alpha=0.05, varfunc = NULL)
plot(m, which=1)
plot(m, which=2)
plot(m, which=3)


## Prediction of x with confidence interval
prediction <- inverse.predict(m, 3500, alpha = 0.05)

## Critical value (decision limit)
crit <- lod(m, alpha = 0.05, beta = 0.05)

m <- lm(y ~ x, data = din32645)
lod(m)
lod(m, alpha = 0.05, beta = 0.05)

m <- lm(y ~ x, data = din32645)
loq(m)
loq(m, n = 3) # better by using replicate measurements



library(EnvStats)
calibrate(y ~ x, data=din32645, test.higher.orders = TRUE, max.order = 4, p.crit = 0.05, F.test = "partial", method = "qr", model = TRUE, x = FALSE, y = FALSE, contrasts = NULL, warn = TRUE)


### nonlinear regression ##################################################################################################

library(investr)
data(Puromycin, package = "datasets")
Puromycin2 <- Puromycin[Puromycin$state == "treated", ][, 1:2]


detach(openxlsx)
write.xlsx(Puromycin, "c:\\Users\\lubop\\Documents\\kurs\\DVpiskyXX.xlsx", col.names=TRUE, row.names=TRUE)

write.table(Puromycin, file = "c:\\Users\\lubop\\Documents\\kurs\\DVpiskyXX.txt", append = FALSE, quote = TRUE, sep = " ",
            eol = "\n", na = "NA", dec = ".", row.names = TRUE,col.names = TRUE)



Puro.nls <- nls(rate ~ Vm * conc/(K + conc), data = Puromycin2,start = c(Vm = 200, K = 0.05))
plotFit(Puro.nls, interval = "both", pch = 19, shade = TRUE, col.conf = "skyblue4", col.pred = "lightskyblue2")


################################################################################################################################3


library(deming)
data(arsenate)
data(ferritin)

# Deming
fit <- deming(aes ~ aas, data=arsenate, xstd=se.aas, ystd=se.aes, conf=.95, jackknife=TRUE, dfbeta=FALSE,x=FALSE, y=FALSE, model=TRUE)
print(fit)
fit$coefficients
fit$coefficients[1]
fit$coefficients[2]
fit$residuals
mean(fit$residuals)


# Passing Bablok
afit1 <- pbreg(aes ~ aas, data= arsenate,conf=.95,nboot = 0, method=1, eps=sqrt(.Machine$double.eps), x = FALSE, y = FALSE, model = TRUE)
print(afit1)
# Theil Sen
afit2 <- theilsen(aes ~ aas, data= arsenate,conf=.95, nboot = 0, symmetric=FALSE, eps=sqrt(.Machine$double.eps), x = FALSE, y = FALSE, model = TRUE)
print(afit2)

data(ferritin)
temp <- ferritin[ferritin$period <4,]
plot(temp$old.lot, temp$new.lot, type='n', log='xy', xlab="Old lot", ylab="New Lot")
text(temp$old.lot, temp$new.lot, temp$period,col=temp$period)


library(mcr)
data(creatinine)
fit.lr <- mcreg(as.matrix(creatinine), method.reg="LinReg", na.rm=TRUE)
fit.wlr <- mcreg(as.matrix(creatinine), method.reg="WLinReg", na.rm=TRUE)
compareFit( fit.lr, fit.wlr )



#### Bland-Altman plot


library(BlandAltmanLeh)
bland.altman.PEFR  # data
bland.altman.plot(group1=arsenate$aes, group2=arsenate$aas, two = 1.96, mode = 1, graph.sys = "base", conf.int = 0, silent = TRUE, sunflower = FALSE, geom_count = FALSE)
bland.altman.stats(group1=arsenate$aes, group2=arsenate$aas, two = 1.96, mode = 1, conf.int = 0.95)

library(MethComp)
library(blandr)


# Linuv koeficient konkordance
library(DescTools)
CCC(arsenate$aes, arsenate$aas, ci = "z-transform", conf.level = 0.95, na.rm = FALSE)
library(epiR)
epi.ccc(arsenate$aes, arsenate$aas, ci = "z-transform", conf.level = 0.95, rep.measure = FALSE, subjectid=NULL)
library(agRee)
ratings = cbind(arsenate$aes, arsenate$aas)
agree.ccc(ratings, conf.level=0.95,method="bootstrap",nboot=999, nmcmc=10000,mvt.para=list(prior=list(lower.v=4, upper.v=25,Mu0=rep(0, ncol(ratings)),Sigma0=diag(10000, ncol(ratings)),p=ncol(ratings),V=diag(1, ncol(ratings))),initial=list(v=NULL, Sigma=NULL)),NAaction="omit")
# method=c("jackknifeZ", "jackknife","bootstrap","bootstrapBC","mvn.jeffreys", "mvn.conjugate","mvt", "lognormalNormal", "mvsn", "mvst")
# NAaction=c("fail", "omit")


############################################################################################################
###### Regrese #####################################################################
############################################################################################################

library(investr)
data(arsenic)

library(ACSWR)
data(viscos)
data=viscos
x = data[,1]
y = data[,2]

library(openxlsx)
sal <- read.xlsx("c:\\Users\\lubop\\Documents\\kurs\\DVpisky1.xlsx", sheet = "regrese2", startRow = 1, colNames = TRUE, rowNames = FALSE,detectDates = FALSE, skipEmptyRows = TRUE,rows = NULL, cols = NULL,check.names = FALSE)
data = as.data.frame(sal)
x = data[,1]
y = data[,2]
z = data[,3]

plot(x,y)


lmMod <- lm(y~x, data=data, weights=NULL, na.action = na.pass)
summary(lmMod)

plot(x,y)
abline(lm(y~x, data=data), col=2, lty=2)
abline(lm(y~x-1, data=data), col=4, lty=2)

cp = predict(lmMod, newdata = 0, interval = 'confidence')
lines(x,cp[,2],col=6,lty=3)
lines(x,cp[,3],col=6,lty=3)

pp = predict(lmMod, newdata = 0, interval = 'prediction')
lines(x,pp[,2],col=3,lty=3)
lines(x,pp[,3],col=3,lty=3)

confint(lmMod,level = 0.95)
confint(lmMod,level = 0.99)

lm.out <- lm(y ~ x)
newx = seq(min(x),max(x),by = 0.05)

conf_interval <- predict(lm.out, newdata=data.frame(x=newx), interval="confidence", level = 0.95)
plot(x, y, xlab="x", ylab="y", main=" ")
abline(lm.out, col="lightblue")
lines(newx, conf_interval[,2], col="blue", lty=2)
lines(newx, conf_interval[,3], col="blue", lty=2)

pred_interval <- predict(lm.out, newdata=data.frame(x=newx), interval="prediction", level = 0.95)
lines(newx, pred_interval[,2], col="green", lty=2)
lines(newx, pred_interval[,3], col="green", lty=2)


# test koeficientu
library(lmtest)
coeftest(lmMod, vcov. = NULL, df = NULL)
coefci(lmMod, parm = NULL, level = 0.95, vcov. = NULL, df = NULL)

library(car)
confidenceEllipse(lmMod)


# fitted
fitted.values(lmMod)

library(visreg)
visreg(lmMod, "x")


# residuals
resid(lmMod) # klasicka
plot(fitted.values(lmMod), resid(lmMod)); abline(h=0, col=2, lty=2)
lines(lowess(fitted.values(lmMod), resid(lmMod)), col = 4)

rstandard(lmMod) # standardizovana
plot(fitted.values(lmMod), rstandard(lmMod)); abline(h=0, col=2, lty=2)
lines(lowess(fitted.values(lmMod), rstandard(lmMod)), col = 4)

rstudent(lmMod) # studentizovana
plot(fitted.values(lmMod), rstudent(lmMod)); abline(h=0, col=2, lty=2)
lines(lowess(fitted.values(lmMod), rstudent(lmMod)), col = 4)

library(s20x)
ciReg(lmMod, conf.level = 0.95, print.out = TRUE)

cooks20x(lmMod, main = "Cook's Distance plot", xlab = "observation number",ylab = "Cook's distance", line = c(0.5, 1.2, 2), cex.labels = 1,axisOpts = list(xAxis = TRUE, yAxisTight = FALSE))

eovcheck(lmMod, xlab = "Fitted values",ylab = "Residuals", col = NULL, smoother = FALSE, twosd = FALSE,levene = FALSE)

modcheck(lmMod, plotOrder = 1:4, args = list(eovcheck = list(smoother = FALSE, twosd = FALSE, levene = FALSE), normcheck = list(xlab = c("Theoretical Quantiles", ""), ylab = c("Sample Quantiles",""), main = c("", ""), col = "light blue", bootstrap = FALSE, B = 5, bpch = 3, bcol = "lightgrey", shapiro.wilk = FALSE, whichPlot = 1:2, usePar =   TRUE), cooks20x = list(main = "Cook's Distance plot", xlab = "observation number", ylab = "Cook's distance", line = c(0.5, 0.1, 2), cex.labels = 1, axisOpts = list(xAxis = TRUE))), parVals = list(mfrow = c(2, 2), xaxs = "r", yaxs = "r", pty = "s", mai= c(0.2, 0.2, 0.05, 0.05)))

normcheck(lmMod, xlab = c("Theoretical Quantiles", ""),ylab = c("Sample Quantiles", ""), main = c("", ""), col = "light blue", bootstrap = FALSE, B = 5, bpch = 3, bcol = "lightgrey", shapiro.wilk = FALSE, whichPlot = 1:2, usePar = TRUE)

residPlot(lmMod, f = 1)


library(car)
residualPlots(lmMod, terms = ~., layout = NULL, main = "", fitted = TRUE, AsIs=TRUE, plot = TRUE,tests = TRUE)
influencePlot(lmMod)
infIndexPlot(lmMod)


library(UsingR)
simple.lm(x, y, show.residuals=TRUE, show.ci=TRUE, conf.level=0.95,pred=NULL) 

par(mfrow=c(2,2))
plot(lmMod)
par(mfrow=c(1,1))

library(lindia)
gg_boxcox(lmMod, showlambda = TRUE, lambdaSF = 3, scale.factor = 1)
gg_cooksd(lmMod, label = TRUE, show.threshold = TRUE,threshold = "convention", scale.factor = 1)
gg_diagnose(lmMod, theme = NULL, ncol = NA, plot.all = TRUE,scale.factor = 1, boxcox = FALSE, max.per.page = NA)
gg_qqplot(lmMod, scale.factor = 1)
gg_resfitted(lmMod, scale.factor = 1)
gg_reshist(lmMod, bins = 20)
gg_resleverage(lmMod, method = "loess", se = FALSE, scale.factor = 1)
gg_resX(lmMod, plot.all = TRUE, scale.factor = 1, max.per.page = NA, ncol = NA)
gg_scalelocation(lmMod, method = "loess", scale.factor = 1, se = FALSE)


# Bonferroni Outlier Test
library(car)
ou = outlierTest(lm(y~x), cutoff=0.05, n.max=15, order=TRUE); ou
as.numeric(names(ou$rstudent))

# Lund Outlier Test
lundout<-function(x1,x2,a) {
  testoutlm<-lm(x1~x2)
  standardresid<-rstandard(testoutlm)
  nl<-length(x1)
  q<-length(testoutlm$coefficients)
  F<-qf(c(1-(a/nl)),df1=1,df2=nl-q-1,lower.tail=TRUE)
  crit<-((nl-q)*F/(nl-q-1+F))^0.5
  oul = which(abs(standardresid)>crit)
  return(oul)
}
lundout(y,x,0.05)


# Shapiro-Wilk test rezidui
shapiro.test(resid(lmMod))
shapiro.test(rstandard(lmMod))
shapiro.test(rstudent(lmMod))


# Breusch-Pagan test na heteroskedasticitu
library(lmtest)
lmtest::bptest(y~x, varformula = NULL, studentize = TRUE, data = data) 

# Harrison-McCabe test for heteroskedasticity
library(lmtest)
lmtest::hmctest(y~x, point = 0.5, order.by = NULL, simulate.p = TRUE, nsim = 1000, plot = TRUE, data = data)

# Score Test for Non-Constant Error Variance
library(car)
car::ncvTest(lmMod)

# Goldfeld-Quandt test against heteroskedasticity
library(lmtest)
lmtest::gqtest(y~x, point = 0.5, fraction = 0, alternative = "two.sided",order.by = NULL, data = data)
# alternative = c("greater", "two.sided", "less")


# Durbin-Watson test for autocorrelation (zavislost rezidui na x)
library(lmtest)
lmtest::dwtest(y~x, order.by = NULL, alternative = "two.sided", iterations = 15, exact = NULL, tol = 1e-10, data = data)
# alternative = c("greater", "two.sided", "less")
library(car)
durbinWatsonTest(lmMod)

### Testy linearity
# Rainbow test linearity
library(lmtest)
lmtest::raintest(y~x, fraction = 0.5, order.by = NULL, center = NULL, data=data)

# Harvey-Collier test for linearity
library(lmtest)
lmtest::harvtest(y~x, order.by = NULL, data = data)



###### regrese pocatkem ############################

s.lm <- lm(y ~ x, data = data)
summary(s.lm)

s.lm2 <- lm(y ~ 0 + x, data = data) # Adding the 0 term tells the lm() to fit the line through the origin
summary(s.lm2)

s.lm3 <- lm(y ~ x-1, data = data) # Adding the 0 term tells the lm() to fit the line through the origin
summary(s.lm3)

plot(x,y) # regrese bez useku
abline(lm(y ~ x-1, data = data), col=2, lty=2)


### Box-Cox #############################################################
library(AID)
lmbc = boxcoxlm(as.matrix(x), y, method = "lse", lambda = seq(-3,3,0.01), lambda2 = NULL, plot = TRUE,alpha = 0.05, verbose = TRUE)
lmbc$lambda.hat
caret::BoxCoxTrans(y)
library(forecast)
ynew = BoxCox(y, lambda=2)
lmMod_bc <- lm(ynew ~ x, data=data)
bptest(lmMod_bc)
plot(lmMod_bc)


### regrese polynomem
s.lmq <- lm(data$y ~ data$x + poly(data$x,2), data = data)
summary(s.lmq)

plot(x,y) # regrese polynomem
lines(lm(y ~ poly(x^2), data = data), col=2, lty=2)
abline(lm(y ~ x + poly(x,2), data = data), col=4, lty=2)

lmMod2 = lm(z ~ x, data = data)

### Srovnani 2 modelu
anova(lmMod,lmMod2)

library(performance)
compare_performance(lmMod, lmMod2, rank = TRUE)

# Srovnani 2 primek
library(gap)
chow.test(y1=y,x1=x,y2=z,x2=x,x=NULL)


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

### Resistant Regression
library(MASS)
lmM1 = lqs(y ~ x, data,method = "lts")
lmM1
# method = c("lts", "lqs", "lms", "S", "model.frame")

lmM2 = rlm(y ~ x, data, weights=NULL, method = "M", wt.method = "inv.var")
lmM2 
# method = c("M", "MM", "model.frame"),
# wt.method = c("inv.var", "case")

library(LearnEDA) # Fits Tukey resistant line of form a + b (x - xC).
lmM3 = rline(x,y,iter=1)


### regrese s LOD
library(lodr)
x2 = 
  fitl <- lod_lm(data=lod_data, frmla=y~x, lod=0.1,var_LOD="x2", nSamples=100, boots=500)
summary(fitl)
coef(fit)
residuals(fitl)


library(truncreg)
truncreg(y~x, data=data, point = 2, direction = "left",model = TRUE, y = FALSE, x = FALSE, scaled = FALSE)


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

library(chemCal)
data(din32645)
din32645
data(massart97ex1)
massart97ex1
data(massart97ex3)
massart97ex3
data(utstats14)


m <- lm(y ~ x, data = din32645)
calplot(m,xlim = c("auto", "auto"), ylim = c("auto", "auto"),xlab = "Concentration", ylab = "Response", alpha=0.05, varfunc = NULL)
plot(m, which=1)
plot(m, which=2)
plot(m, which=3)


## Prediction of x with confidence interval
prediction <- inverse.predict(m, 3500, alpha = 0.05)

## Critical value (decision limit)
crit <- lod(m, alpha = 0.05, beta = 0.05)

m <- lm(y ~ x, data = din32645)
lod(m)
lod(m, alpha = 0.05, beta = 0.05)

m <- lm(y ~ x, data = din32645)
loq(m)
loq(m, n = 3) # better by using replicate measurements



library(EnvStats)
calibrate(y ~ x, data=din32645, test.higher.orders = TRUE, max.order = 4, p.crit = 0.05, F.test = "partial", method = "qr", model = TRUE, x = FALSE, y = FALSE, contrasts = NULL, warn = TRUE)


### nonlinear regression ##################################################################################################

library(investr)
data(Puromycin, package = "datasets")
Puromycin2 <- Puromycin[Puromycin$state == "treated", ][, 1:2]


detach(openxlsx)
write.xlsx(Puromycin, "c:\\Users\\lubop\\Documents\\kurs\\DVpiskyXX.xlsx", col.names=TRUE, row.names=TRUE)

write.table(Puromycin, file = "c:\\Users\\lubop\\Documents\\kurs\\DVpiskyXX.txt", append = FALSE, quote = TRUE, sep = " ",
            eol = "\n", na = "NA", dec = ".", row.names = TRUE,col.names = TRUE)



Puro.nls <- nls(rate ~ Vm * conc/(K + conc), data = Puromycin2,start = c(Vm = 200, K = 0.05))
plotFit(Puro.nls, interval = "both", pch = 19, shade = TRUE, col.conf = "skyblue4", col.pred = "lightskyblue2")


################################################################################################################################3


library(deming)
data(arsenate)
data(ferritin)

# Deming
fit <- deming(aes ~ aas, data=arsenate, xstd=se.aas, ystd=se.aes, conf=.95, jackknife=TRUE, dfbeta=FALSE,x=FALSE, y=FALSE, model=TRUE)
print(fit)
fit$coefficients
fit$coefficients[1]
fit$coefficients[2]
fit$residuals
mean(fit$residuals)


# Passing Bablok
afit1 <- pbreg(aes ~ aas, data= arsenate,conf=.95,nboot = 0, method=1, eps=sqrt(.Machine$double.eps), x = FALSE, y = FALSE, model = TRUE)
print(afit1)
# Theil Sen
afit2 <- theilsen(aes ~ aas, data= arsenate,conf=.95, nboot = 0, symmetric=FALSE, eps=sqrt(.Machine$double.eps), x = FALSE, y = FALSE, model = TRUE)
print(afit2)

data(ferritin)
temp <- ferritin[ferritin$period <4,]
plot(temp$old.lot, temp$new.lot, type='n', log='xy', xlab="Old lot", ylab="New Lot")
text(temp$old.lot, temp$new.lot, temp$period,col=temp$period)


library(mcr)
data(creatinine)
fit.lr <- mcreg(as.matrix(creatinine), method.reg="LinReg", na.rm=TRUE)
fit.wlr <- mcreg(as.matrix(creatinine), method.reg="WLinReg", na.rm=TRUE)
compareFit( fit.lr, fit.wlr )



#### Bland-Altman plot


library(BlandAltmanLeh)
bland.altman.PEFR  # data
bland.altman.plot(group1=arsenate$aes, group2=arsenate$aas, two = 1.96, mode = 1, graph.sys = "base", conf.int = 0, silent = TRUE, sunflower = FALSE, geom_count = FALSE)
bland.altman.stats(group1=arsenate$aes, group2=arsenate$aas, two = 1.96, mode = 1, conf.int = 0.95)

library(MethComp)
library(blandr)


# Linuv koeficient konkordance
library(DescTools)
CCC(arsenate$aes, arsenate$aas, ci = "z-transform", conf.level = 0.95, na.rm = FALSE)
library(epiR)
epi.ccc(arsenate$aes, arsenate$aas, ci = "z-transform", conf.level = 0.95, rep.measure = FALSE, subjectid=NULL)
library(agRee)
ratings = cbind(arsenate$aes, arsenate$aas)
agree.ccc(ratings, conf.level=0.95,method="bootstrap",nboot=999, nmcmc=10000,mvt.para=list(prior=list(lower.v=4, upper.v=25,Mu0=rep(0, ncol(ratings)),Sigma0=diag(10000, ncol(ratings)),p=ncol(ratings),V=diag(1, ncol(ratings))),initial=list(v=NULL, Sigma=NULL)),NAaction="omit")
# method=c("jackknifeZ", "jackknife","bootstrap","bootstrapBC","mvn.jeffreys", "mvn.conjugate","mvt", "lognormalNormal", "mvsn", "mvst")
# NAaction=c("fail", "omit")


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


library(EBImage) # Bioc
Eg1 <- readImage("c:\\Users\\lubop\\Documents\\image\\Egypt1.jpg")
display(Eg1)
print(Eg1)
img <- Eg1

# Brightness
img1 = img+0.4
display(img1)
img1 = img-0.4
display(img1)

# Contrast
img2 = 2*img
display(img2)
img2 = .5*img
display(img2)

# Gamma correction
img3 = img^3
img3 = img^.7
img3 = (0.2+img)^3
display(img3)

# Treshold
img4 = img>0.5
display(img4)

imgc=img
colorMode(imgc) = Grayscale
display(imgc)

colorMode(imgc) = Color
display(imgc)

imgb = rgbImage(red=img, green=img*0, blue=img*0)
display(imgb)
imgb = rgbImage(red=img*0, green=img, blue=img*0)
display(imgb)
imgb = rgbImage(red=img*0, green=img*0, blue=img)
display(imgb)
imgb = rgbImage(red=img>0.6, green=img*0.5, blue=img*1.5)
display(imgb)

hist(img)


## Low-pass disc-shaped filter
flo = makeBrush(11,shape='disc', step=FALSE)^2
flo = flo/sum(flo)
imgflo = filter2(imgc, flo)
display(imgflo)

## High-pass Laplacian filter
fhi = matrix(1, nc=3, nr=3)
fhi[2,2] = -8
imgfhi = filter2(imgc, fhi)
display(imgfhi)
display(1-imgfhi)
display(.8-imgfhi)

## Low-pass Gaussian filter
sigma=6
imgflg = gblur(img, sigma, radius = 2 * ceiling(3 * sigma) + 1)
display(imgflg)


## Median filter
imgm = medianFilter(img, 3)
display(imgm)


red=as.vector(Eg1[,,1])
green=as.vector(Eg1[,,2])
blue=as.vector(Eg1[,,3])
plot(density(red),col=2,ylim=c(0,6.5),lwd=2,main="",sub="",cex.axis=1.5,cex.lab=1.5)
points(density(green),col=3,lty=1,type="l",lwd=2)
points(density(blue),col=4,lty=1,type="l",lwd=2)


library(squash)
Mtr1 <- xyzmat2xyz(img[,,1])
Mtr2 <- xyzmat2xyz(img[,,2])
Mtr3 <- xyzmat2xyz(img[,,3])
Mtr <- cbind(Mtr1$z,Mtr2$z,Mtr3$z)


library(e1071)
N <- 4
fc<-cmeans(Mtr,centers=N,iter.max=100,verbose=FALSE,dist="euclidean",method="cmeans",m=1.1,rate.par=NULL,weights=1)
fcl<-fc$cluster
xyz <- cbind(Mtr1$x,Mtr1$y[length(Mtr1$y):1],fcl)

library(raster)
Eg1tr <- rasterFromXYZ(xyz)

image(Eg1tr,col=c(1:N),add=FALSE,xaxs="i",yaxs="i",axes=FALSE,xlab="",ylab="")

image(Eg1tr,col=c(1,0,0,0,0,0,0),add=FALSE,xaxs="i",yaxs="i",axes=FALSE,xlab="",ylab="")
image(Eg1tr,col=c(0,2,0,0,0,0,0),add=FALSE,xaxs="i",yaxs="i",axes=FALSE,xlab="",ylab="")
image(Eg1tr,col=c(0,0,0,4,0,0,0),add=FALSE,xaxs="i",yaxs="i",axes=FALSE,xlab="",ylab="")
image(Eg1tr,col=c(0,0,3,0,0,0,0),add=FALSE,xaxs="i",yaxs="i",axes=FALSE,xlab="",ylab="")
image(Eg1tr,col=c(0,0,0,0,5,0,0),add=FALSE,xaxs="i",yaxs="i",axes=FALSE,xlab="",ylab="")
image(Eg1tr,col=c(0,0,0,0,0,6,0),add=FALSE,xaxs="i",yaxs="i",axes=FALSE,xlab="",ylab="")
image(Eg1tr,col=c(0,0,0,0,0,0,7),add=FALSE,xaxs="i",yaxs="i",axes=FALSE,xlab="",ylab="")

fclm = fcl*0
fclm[which(fcl==1)]= "black"
fclm[which(fcl==2)]= "red"
fclm[which(fcl==3)]= "green3"
fclm[which(fcl==4)]= "blue"
fclm[which(fcl==5)]= "cyan"
fclm[which(fcl==6)]= "magenta"
fclm[which(fcl==7)]= "yellow"
qwfc = Image(fclm,dim=c(dim(img)[1],dim(img)[2]))
display(qwfc)           

fclm = rep("white",length(fcl))
fclm[which(fcl==1)]= "black"
qwfc1 = Image(fclm,dim=c(dim(img)[1],dim(img)[2]))
display(qwfc1)

CL = 3
fcmem<-fc$membership
fcm<-fcmem[,CL]
xyz <- cbind(Mtr1$x,Mtr1$y[length(Mtr1$y):1],fcm)
library(raster)
Eg1tM <- rasterFromXYZ(xyz)
image(Eg1tM,col=gray.colors(12, start = 0.1, end = 0.9, gamma = 2.2, alpha = NULL),add=FALSE,xaxs="i",yaxs="i",axes=FALSE,xlab="",ylab="")
qwfcm = Image(fcm,dim=c(dim(img)[1],dim(img)[2]))
display(qwfcm)



#####################################################################################################################################]

rgb(0, 1, 0)
rgb(0.3, 0.7, 0.9)
col2rgb(paste0("gold", 1:4))
col2rgb("peachpuff")
col2rgb(c(blu = "royalblue", reddish = "tomato")) 


library(plotrix)
color.id("#cc00cc")
sapply(rainbow(4), color.id)


library(ggtern)
rgb2hex(crgb[,1])
rgb2hex(255,0,0)

library(gplots)
col2hex(c("red","yellow","lightgrey"))


library(jjb)
rgb_to_hex(255,255,255) # Hexadecimal with pound sign
rgb_to_hex(255,255,255,FALSE) # Heaxadecimal without pound sign
shade(c(22, 150, 230), shade_factor = 0.5)
tint(rgb_value, tint_factor = 0.2)


library(aqp)
# Napthol Red PR112 (8R-5-14)
munsell2rgb(the_hue="8R", the_value=5, the_chroma=14, return_triplets=TRUE)
munsell2rgb(the_hue=c('10YR', '2.5YR'), the_value=c(3, 5), the_chroma=c(5, 6), return_triplets=TRUE)


library(munsellinterpol)
# Napthol Red PR112 (8R-5-14)
# Ultramarine Blue Reddish PB29 (7PB-2-12)
hns = HueNumberFromString( '8R' )
(mrgb = MunsellToRGB( c(hns,2,12), space='sRGB', maxSignal=255, adapt='Bradford'))
rgbc = as.numeric(unlist(mrgb)[5:7])
library(jjb)
col= rgb_to_hex(rgbc[1],rgbc[2],rgbc[3])
barplot(5,col=col, axes = FALSE)


# Transmitance to absorbance
library(photobiology)
T2A(0.99, action = NULL, byref = FALSE, clean = TRUE) # transmittance into absorbance

(crgb = col2rgb(w_length2rgb(c(400, 500, 600, 700), color.name=c("a","b","c","d"))))

col2rgb(w_length_range2rgb(500:600))
col2rgb(w_length_range2rgb(c(500,600)))