# pobieranie R
# R download

# https://cran.r-project.org/


# RStudio IDE

# https://rstudio.com/products/rstudio/


# utworz macierz 2x3 gdzie elementy sa wypisane wierszami 
# create 3x3 matrix where elements are listed by rows

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

# utworz macierz 2x2 gdzie elementy sa wypisane wierszami 
# create 2x2 matrix where elements are listed by rows

B <- matrix(c(1,2,-1,4),nrow=2,ncol=2,byrow=TRUE)
B

# iloczyn macierzy B i A
# the product of B by A

B %*% A


# blad - liczba kolumn macierzy A nie jest rowna liczbie wierszy macierzy B
# error - the number of rows of matrix A is different from 
# the number of rows of matrix B

A %*% B


# zainstaluj pakiet "expm"
# install library "expm"

install.packages("expm")


# zaladuj pakiet "expm"
# load library "expm"

library(expm)


# macierz B do potegi 2 
# (potrzebuje pakietu "expm")
# matrix B raised to the power of 2
# (requires library "expm")

B %^% 2


# macierz A transponowana
# matrix A transposed 

t(A)


# pierwszy wiersz macierzy A
# the first row of matrix A

A[1,]


# trzecia kolumna macierzy A
# the third column of matrix A

A[,3]


# dwie ostatnie kolumny macierzy A
# the two last columns of matrix A

A[,2:3]


# macierze A i B polaczone kolumnami
# matrix A joined with B by their columns

cbind(A,B)


# macierz A polaczona ze swoja kopia wierszami
# matrix A joined with itself by rows

rbind(A,A)


# odwrotnosc macierzy B (e-16 znaczy 10^(-16))
# the inverse of matrix B (e-16 stands for 10^(-16))

solve(B)
solve(B) %*% B
B %*% solve(B) 


# zaladuj pakiet "pracma"
# load library "pracma"

library(pracma)


# liczba 1 powtorzona 2 razy
# the number 1 repeated twice

rep(1,2)


# macierz diagonalna 3x3 z liczbami 1,2,3 na przekatnej
# the diagonal matrix with number 1,2,3 on the diagonal

diag(c(1,2,3))


# macierz jednostkowa 3x3
# the 3x3 unit matrix 

diag(3)
diag(rep(1,3))


# postac schodkowa zredukowana macierzy A
# (potrzebuje pakietu "pracma")
# the (row) reduced echelon form of matrix A
# (requires library "pracma")

rref(A)


# znajdz wspolrzedne wektora v=(1, 8, 10, 10) 
# w bazie  (1, 2, 3, 1), (2, 1, 3, 3),(-1, 1, 0, -1), (0, 0, 1, 2)
# find coordinates of the vector v=(1, 8, 10, 10) 
# relative to the basis (1, 2, 3, 1), (2, 1, 3, 3),(-1, 1, 0, -1), (0, 0, 1, 2)

M <- matrix(c(1,2,3,1,2,1,3,3,-1,1,0,-1,0,0,1,2),nrow=4,ncol=4,byrow=FALSE)
v <- c(1,8,10,10)
coords <- rref(cbind(M,v))[,5]
coords


# znajdz baze podprzestrzeni zadanej jednorodnym ukladem rownan Ax=0
# find basis of a subspace given by a homogeneous system of linear equations Ax=0

remFirstCol <- function(M) {
	if (ncol(M)==1) { return(matrix(,nrow=nrow(M),ncol=1)) } else if (ncol(M)==2) {
	return(as.matrix(M[,-1]))} else { return(M[,-1]) }
}

getFirstCol <- function(M) {
	return(as.matrix(M[,1]))
}

doShift <- function(M) {
	if (nrow(M)==1) return(M) else {
		m <- nrow(M)
		return(as.matrix(M[c(m,1:(m-1)),]))}
}

if (rankMatrix(A)[1]==nrow(A)) {print("zero subspace has empty basis")
	} else { 	rrefA <- rref(A)[1:rankMatrix(A)[1],]
				pivots <- apply(rrefA,1,function(row) which.max(row==1))
				params <- -t(rrefA[,-pivots]) # remaining columns
				d <- ncol(A)-nrow(rrefA)
				colIns <- zeros(d,1)
				colIns[1,1] <- 1
				ans <- colIns
				for (j in 1:ncol(A)) {
						if (j %in% pivots){	ans <- cbind(ans,getFirstCol(params))
											params <- remFirstCol(params)
							} else { 	ans <- cbind(ans,colIns)
										colIns <- doShift(colIns)
									}	
				}
				ans[,-1]
				# rows give a basis, number of rows is equal to the dimension
	}
# try on A <- matrix(c(1,-2,0,4,0,0,1,7,1,-2,1,11),nrow=3,ncol=4,byrow=TRUE)
# try on A <- matrix(c(1,-2,0,4,11,0,0,1,7,-12,1,-2,1,11,-1),nrow=3,ncol=5,byrow=TRUE)


# odwrotnosc macierzy B obliczona metoda I
# (potrzebuje pakietu "pracma")
# the inverse of matrix B computed with method I
# (requires library "pracma")

rref(cbind(B,diag(2)))[,3:4]


# wyznacznik macierzy B
# the determinant of matrix B

det(B)


# macierz exp(B)=I+A+A^2/2!+A^3/3!+..
# (potrzebuje pakietu "expm")
# the matrix exp(B)=I+A+A^2/2!+A^3/3!+..
# (requires library "expm")

expm(B)


# rzad macierzy A
# (potrzebuje pakietu "Matrix")
# the rank of matrix A
# (requires library "Matrix")

r <- rankMatrix(A)
r[1]


# rozklad QR macierzy A
# the QR decomposition of matrix A

q <- qr(A)
Q <- qr.Q(q)
R <- qr.R(q)
Q %*% R 	# gives back A
Q %*% t(Q)	# matrix Q is orthogonal 


# rozklad macierzy A wedlug wartosci osobliwych (rozklad SVD)
# the SVD decomposition of matrix A

s <- svd(A)
names(s)
u <- s$u
v <- s$v
sig <- cbind(diag(s$d),rep(0,2))
v <- cbind(v,rep(0,3))
u %*% sig %*% t(v) # gives back A


# najlepsze przyblizenie macierzy A macierzami rangi 1 
# w normie Frobeniusa tzn. ||A-A1|| -> min, 
# gdzie ||A||^2=tr(A^T A)=\sum_{i,j} a_{ij}^2
# the best rank 1 approximation of matrix A 
# with respect to the Frobenius norm

sig1 <- cbind(diag(c(s$d[-2],0)),rep(0,2))
u %*% sig1 %*% t(v)
rankMatrix(u %*% sig1 %*% t(v))


# macierz dolaczona macierzy B
# the adjoint matrix of matrix B

adj <- function(B)
{
if (ncol(B)!=nrow(B)) 
	{	
		print("nie jest macierza kwadratowa/not a square matrix")
	} else
	{
		temp <- array(0,dim=c(ncol(B),ncol(B)))
		for(i in 1:ncol(B))
		{
			for(j in 1:ncol(B))
			{
				temp[i,j] <- (-1)^(i+j)*det(as.matrix(B[-i,-j]))
			}
		}
		t(temp)
	}
}

adj(B)


# macierz odwrotna macierzy B z macierzy dolaczonej, metoda II
# the inverse matrix of matrix B by using the formula with adjoint matrix, method II

(1/det(B))*adj(B)
solve(B)
(1/det(B))*adj(B) %*% B


# wektory i wartosci wlasne macierzy B
# eigenvectors and eigenvalues of matrix B

e <- eigen(B)
names(e)
e$values
e$vectors

# wektory wlasne w kolumnach C=M(id)_A^st
# eigenvectors in columns C=M(id)_A^st

C <- e$vectors
D <- diag(e$values)
C %*% D %*% solve(C)	# daje macierz B/recovers matrix B
B


# macierz B^10
# matrix B^10

C <- e$vectors
D <- diag(e$values)
C %*% diag(c(3^10,2^10)) %*% solve(C)	
B %^% 10


# macierz B^100 
# matrix B^100

library(gmp) 	# library for big rational numbers
power <- 100
gD <- matrix(c(as.bigq(3)^power,0,0,as.bigq(2)^power),nrow=2,ncol=2)
gC <- as.bigq(C)
gB <- as.bigq(B)
gC %*% gD %*% solve(gC)
B%^%power

gA <- as.bigq(diag(2))
for (i in 1:power) { gA <- gA %*% gB }
gA


# macierz, ktora nie jest diagonalizowalna nad liczbami rzeczywistym
# ale jest diagonalizowalna nad liczbami zespolonymi
# matrix which is not diagonalizable over the real numbers
# but it is diagonalizable over the complex numbers

B1 <- matrix(c(0,-1,1,0),nrow=2,ncol=2,byrow=TRUE)
e <- eigen(B1)
e

# macierz, ktora nie jest diagonalizowalna nad liczbami rzeczywistym
# ani nad liczbami zespolonymi
# matrix which is not diagonalizable over the real numbers
# and over the complex numbers

B2 <- matrix(c(1,-1,0,1),nrow=2,ncol=2,byrow=TRUE)
e <- eigen(B2)
e
rankMatrix(e$vectors)


# losowa macierz jest diagonalizowalna nad liczbami zespolonymi
# random matrix is diagonalizable over the complex numbers

n <- 3
M <- 100
range <- 10 
# M losowych macierzy nxn o elementach od -range do range
# M random nxn matrices with entries from -range to range
eigenvNumber <- replicate(M,rankMatrix(eigen(matrix(runif(n*n,-range,range),ncol=n,nrow=n))$vectors)) 
# zawiera liczby wektorow wlasnych M losowych macierzy
# contains the number of eigenvectors of M random matrices
isDiag <- eigenvNumber==n
isDiag


# losowa macierz ma maksymalny rzad
# random matrix has maximal rank

m <- 3
n <- 5
M <- 100
range <- 10
# M losowych macierzy mxn o elementach od -range do range
# M random mxn matrices with entries from -range to range
randMatrices <- replicate(M,matrix(runif(m*n,-range,range),ncol=n,nrow=m))
# randMatrices[,,1]
# randMatrices[,,1:3]
# zastosowanie funkcji rankMatrix do 3-ciej wspolrzednej randMatrices
# apply function rankMatrix to the third coordinate of randMatrices
ranks <- apply(randMatrices,3,rankMatrix)
isMaximal <- ranks == min(n,m)
isMaximal


# iloczyn skalarny
# scalar/dot product
# (needs library "pracma")

dot(c(1,1,1),c(1,2,3))

# wektor unormowany
# normed vector
# (needs library "pracma")

x <- c(1,2,2)
x/sqrt(dot(x,x))


# ortogonalizacja Grama-Schmidta
# Gram-Schmidt process
# V=lin((1, 0, 1, 0),(0, 1, 0, 2),(2, -2, 2, -4))
# (needs library "pracma")

A <- matrix(c(1,0,1,0,0,1,0,2,2,-2,2,-4),nrow=3,ncol=4,byrow=TRUE)
A <- t(rref(A)[1:2,])
A
# kolumny A tworza baze podprzestrzeni V
# columns of A gives a basis of subspace V
g <- gramSchmidt(A)
g$Q 
# kolumny Q tworza baze ortonormalna podprzestrzeni V
# columns of Q give an orthonormal basis of subspace V


# rzut prostopadly wektora w=(1,1,1,1) na podprzestrzen V
# the (linear) orthogonal projection of vector w=(1,1,1,1) onto the subspace V
#V=lin((1, 0, 1, 0),(0, 1, 0, 2),(2, -2, 2, -4))

A <- matrix(c(1,0,1,0,0,1,0,2,2,-2,2,-4),nrow=3,ncol=4,byrow=TRUE)
A <- t(rref(A)[1:2,])
A
# kolumny A tworza baze podprzestrzeni V
# columns of A give a basis of subspace V
w <- rep(1,4)
A%*%solve(t(A)%*%A)%*%t(A)%*%w


# jak wyzej, przy uzyciu pakietu "matlib"
# as above using library "matlib" 
# (needs library "matlib")

library(matlib)
A <- matrix(c(1,0,1,0,0,1,0,2,2,-2,2,-4),nrow=3,ncol=4,byrow=TRUE)
w <- rep(1,4)
Proj(w,t(A))


# jak wyzej, przy uzyciu bazy ortonormalnej V
# as above using orthonormal basis of V
# (needs library "pracma")

A <- matrix(c(1,0,1,0,0,1,0,2,2,-2,2,-4),nrow=3,ncol=4,byrow=TRUE)
A <- t(rref(A)[1:2,])
g <- gramSchmidt(A)
w <- rep(1,4)
g$Q %*% t(g$Q)%*% w 
# macierz t(q$Q)%*% w zawiera iloczyny skalarne w z kolejnymi wektorami bazy ortonormalnej
# matrix t(q$Q)%*% w contains scalar products of vector w with vectors of orthonormal basis


# symetria prostopadla wektora w=(1,1,1,1) wzgledem podprzestrzeni V
# the (linear) orthogonal reflection of vector w=(1,1,1,1) about the subspace V
#V=lin((1, 0, 1, 0),(0, 1, 0, 2),(2, -2, 2, -4))

A <- matrix(c(1,0,1,0,0,1,0,2,2,-2,2,-4),nrow=3,ncol=4,byrow=TRUE)
A <- t(rref(A)[1:2,])
A
# kolumny A tworza baze podprzestrzeni V
# columns of A gives a basis of subspace V
w <- rep(1,4)
2*A%*%solve(t(A)%*%A)%*%t(A)%*%w-w


# zaladuj pakiet "matrixcalc"
# load library "matrixcalc"
# install.packages("matrixcalc")

library(matrixcalc)


# macierz symetryczna M jest ujemnie określona
# symmetric matrix M is negative definite
# (needs library "matrixcalc")

M=matrix(c(-1,2,2,-5),nrow=2,ncol=2,byrow=TRUE)
is.negative.definite(M)
all(eigen(M)$values<0)


# macierz A^T A jest dodatnio polokreslona
# matrix A^T A is positive semidefinite
# (needs library "matrixcalc")

A=matrix(c(1,2,3,4,5,6),nrow=2,ncol=3,byrow=TRUE)
is.positive.semi.definite(t(A) %*% A)
all(eigen(t(A) %*% A)$values>=0)


# macierz A^T A jest nie dodatnio okreslona
# matrix A^T A is positive definite
# (needs library "matrixcalc")

A=matrix(c(1,2,3,4,5,6),nrow=2,ncol=3,byrow=TRUE)
is.positive.definite(t(A) %*% A)
all(eigen(t(A) %*% A)$values>0)


# znajdz wszystkie zbiory i rozwiazania dopuszczalne
# find all basic feasible sets and corresponding feasible solutions

# problem 55d
# A <- matrix(c(1,2,3,1,,2,5,1,1,1,1,1,-1),nrow=3,byrow=TRUE)
# b <- as.matrix(c(0,1,2))
# problem 55c - one equation is redundant
A <- matrix(c(1,1,-1,2,-1,2,1,-2,2,3,-2,4),nrow=3,byrow=TRUE)
b <- as.matrix(c(1,2,3))
rrefAb <- rref(cbind(A,b))
A <- rrefAb[1:2,1:4]
b <- as.matrix(rrefAb[1:2,5])
# problem 55b is inconsistent
# Ab <- matrix(c(1,1,2,1,2,2,1,4,1,1,3,-1,6,-1,0),nrow=3,byrow=TRUE)
# rref(Ab)
# problem 55a
# A <- matrix(c(1,2,3,-1,1,2,5,6,1,-1),nrow=2,byrow=TRUE)
# b <- as.matrix(c(1,2))

n <- ncol(A)
m <- nrow(A)
rankMatrix(A)==rankMatrix(cbind(A,b))
# install.packages("combinat")
require(combinat)
sets <- combn(n,m)
for (i in 1:(ncol(sets))) 
	{ 
		H <- A[,sets[,i]]
		# print(H)
		if (det(H)!=0) {	if (all((solve(H) %*% b)>=0)) sFeas <- "feasible" else sFeas <- "infeasible"
							sSet <- paste(as.character(sets[,i]),collapse=", ")
							sSolution <- paste(as.character(solve(H) %*% b),collapse=", ")
							cat(paste("B={",sSet,"} is basic ",sFeas," with x_B=(",sSolution,")\n",sep="")) }
	}


# zaladuj pakiet "Rglpk"
# load library(Rglpk)
# install.packages("Rglpk")

library(Rglpk)


# rozwiaz problem programowania liniowego, zadanie 56a
# solve linear programming problem 56a

# funkcja celu/the objective function
f.obj <- c(-1,-2)
# macierz wiezow A/the constraints matrix A
f.con <- matrix(c(4,4),nrow=1,byrow=TRUE)
# nierownosci/inequalities
f.dir <- c("<=")
# prawa strona nierownosci/the right hand side of inequalities		   
f.rhs <- c(12)
# dolne i gorne ograniczenia na wartosci zmiennych
# lower and upper variables bounds
f.bounds <- list(lower = list(ind = c(1,2), val = c(0,0)),
upper = list(ind = c(1,2), val = c(2, 2)))
# minimum funkcji celu na zbiorze dopuszczalnym
# the minimal value of the objective function on the feasible set
Rglpk_solve_LP(f.obj,f.con,f.dir,f.rhs,f.bounds,max=FALSE)
# rozwiazanie optymalne - punkt/wierzcholek w ktorym przyjmowane jest minium
# optimal solution - point/vertex in which the minimum is attained
Rglpk_solve_LP(f.obj,f.con,f.dir,f.rhs,f.bounds,max=FALSE)$solution


# zadanie 56a w postaci standardowej
# problem 56a in standard form

f.obj <- c(-1,-2,0,0,0)
f.con <- matrix(c(4,4,1,0,0,
                  1,0,0,1,0,
				  0,1,0,0,1),nrow=3,byrow=TRUE)
f.dir <- c("==",
           "==",
           "==")
f.rhs <- c(12,
           2,
           2)
f.bounds <- list(lower = list(ind = 1:5, val = rep(0,5)))
Rglpk_solve_LP(f.obj,f.con,f.dir,f.rhs,f.bounds,max=FALSE)
Rglpk_solve_LP(f.obj,f.con,f.dir,f.rhs,f.bounds,max=FALSE)$solution
# rozwiazanie optymalne pokrywa sie z poprzednim 
# po opuszczeniu wartosci pomocnicznych zmiennych x_3,x_4,x_5
# one can recover original optimal solution by forgetting 
# the values of slack variables x_3,x_4,x_5


# zadanie 56b
# problem 56b

f.obj <- c(2,1)
f.con <- matrix(c(4,4),nrow=1,byrow=TRUE)
f.dir <- c("<=")
f.rhs <- c(12)
f.bounds <- list(lower = list(ind = c(1,2), val = c(0,0)),
upper = list(ind = c(1,2), val = c(2, 2)))
Rglpk_solve_LP(f.obj,f.con,f.dir,f.rhs,f.bounds,max=TRUE)
Rglpk_solve_LP(f.obj,f.con,f.dir,f.rhs,f.bounds,max=TRUE)$solution


# zadanie 56c
# problem 56c

f.obj <- c(3,-2)
f.con <- matrix(c(1,1,
                 -3,2,
				  1,-1),nrow=3,byrow=TRUE)
f.dir <- c(">=",
           ">=",
           "<=")
f.rhs <- c(4,
           8,
           0)
f.bounds <- list(lower = list(ind = c(1,2), val = c(0,0)))
Rglpk_solve_LP(f.obj,f.con,f.dir,f.rhs,f.bounds,max=TRUE)
Rglpk_solve_LP(f.obj,f.con,f.dir,f.rhs,f.bounds,max=TRUE)$solution

# eksportuj do LaTeXa
# export to LaTeX

library(xtable)
print(xtable(A))