### 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(-Rv^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)
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
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)


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

library(rSymPy)

# 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

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