#---------------------------------------------------------------#
## relenar una matriz, pir línea, 4 líneas                     ##

M <- matrix (byrow = TRUE,              # por línea
             nrow = 4,                  # num de línea
             c( 40, 1 , 75000 ,
                30, 0 , 69000 ,
                20, 0 , 69000 ,
                35, 1 , 75000 )         
             )
M                                       
# segunda linea
# y tercera columna
M[2.3]

M[2, ] # segunda línea
is.matrix(M[2, ]) # ¡no es una matriz!

m2_ <- matrix(byrow = TRUE, # una matriz de una línea
              fila = 1,
              M[2, ])
m2_
is.matrix(m2_) # esta es una matriz
dim(m2_) # dimensión una línea, 3 columnas

m_3 <- matrix(byrow = TRUE, # una matriz de una columna
                nrow = 4,
              M[ ,3])
m_3
is.matrix(m_3) # esto es una matriz
dim(m_3) # una columna de 3 lineas

m_3 <- matrix(M[ ,3]) # col por defecto
m_3

m2_ <- matrix(M[2, ]) #! columna por defecto
m2_#PB obtenemos una columna

m2_ <- t(matrix(M[2, ])) # transponer
m2_

#------------------------------------------------- --------------#
## manipulación de vectores versus vectores y matrices ##

v <- c(1,2,7,10) # creando un vector
v
dim(v) # sin dimensión
length(v) # su longitud
str(v) # en detalle

2*v # producto en cada componente
2*v + 100# aún mejor
v^2 # cuadrado en cada componente

w <- c(-10,-20,-30,-100) # otro vector
v+w # añadir miembro a miembro
v*w # producto miembro a miembro atención

sum(v) # suma de elementos
mean(v) # media de los elementos

sum(v*w) # producto escalar v, w
v%*%w # producto escalar v, w

is.vector(sum(v*w)) # esto es un vector
is.vector(1789) # un escalar simple es un vector
is.vector(v%*%w) # esto no es un vector pero
is.matrix(v%*%w) # ¡es una matriz!


f <- function(x){x^(1.001)-0.001} # cualquier función
f(10) # prueba de función

u <- apply(matrix(v),1,f) # aplicación en un vector
u
is.vector(u) # esto es un vector

u <- apply(matrix(v),2,f) # aplicación en un vector
u
is.vector(u) # esto no es un vector pero
is.matrix(u) # esto es una matriz

N <- apply(M,c(1,2),f) # aplicación en una matriz
N
is.matrix(N) # esta es una matriz

t(matrix(apply(M,2,mean))) # promedio para columna

T <- t(M)
T
matrix(apply(T,1,mean)) # media por línea

t(matrix(apply(M,2,sum))) # suma de cada columna
t(matrix(apply(M,2,min ))) # min de cada columna
t(matrix(apply(M,2,max ))) # max de cada columna
t(matrix(apply(M,2,stdev)))# desviación típica de cada columna
t(matrix(apply(M,2,prod ))) # produce términos por columna

matrix(apply(T,1,sum )) # suma de cada línea
matrix(apply(T,1,min )) # min de cada línea
matrix(apply(T,1,max )) # max de cada línea
matrix(apply(T, 1, stdev)) # desviación estándar de cada línea
matrix(apply(T,1,prod )) # produce términos por línea



#---------------------------------------------------------------#

M-N # resta miembro a miembro
2*M # multiplicación por 2 de miembros
M[2:3,2:3] # extraer submatriz

M[2:3,2:3] <- 2* M[2:3,2:3] # reemplazando valores en una
M       # parte de la matriz 

#------------------------------------------------- --------------#
x <- matrix( seq( 0, 4 ) ) # primeros 5 enteros
print(x)
is.matrix(x) # de hecho es una matriz
y <- matrix( seq( 5, 9 ) ) # los 5 siguientes
print(y)
print(x+y) # funciona porque dim(x) = dim(y)

#------------------------------------------------- --------------#
# MATRICES, CÁLCULOS UN POCO + RAROS #

library(matrixcalc) # para hacer cálculos matriciales

#------------------------------------------------- --------------#
# matrices conocidas #

H <- creation.matrix( 5 ) # subdiagonal 1 a n-1
print(H)
is.matrix(H) # de hecho es una matriz

V <- vandermonde.matrix(c(seq(1,5)),5) # matriz de vandermonde
printV)
is.matrix(V) # de hecho es una matriz


#------------------------------------------------- --------------#
# producto Hadamard #
# producto miembro * miembro: zij = aij * bij #
Z <- hadamard.prod (V,H)
print(Z)

# es lo mismo que
Z <- V * H # producto de Hadamard
print(Z)


# ¡no debe confundirse con el producto VERDADERO de matrices!
# producto zij = sum_over_k(aik * bkj) #
Z <- V %*% H # líneas * columnas
is.matrix(Z) # de hecho es una matriz
print(Z)

x%*%t(y) # obtenemos una matriz
t(x)%*%y # producto escalar
t(y)%*%x # ídem
sqrt(t(x)%*%x) # norma euclidiana

#------------------------------------------------- --------------#
# operación por bloque en matrices #

#------------------------------------------------- --------------#
# suma directa de matrices #
# es una operación por bloques, a partir de las matrices A y B #
# obtenemos una diagonal A B y el resto con 0s #
Z <- direct.sum(x, y) # suma "directa"
print(Z) # 2 columnas de dim (10)

Z <- direct.sum(V, H) # suma "directa"
print(Z) # 10 x 10

#------------------------------------------------- --------------#
# producto de Kronecker = producto directo (tensorial) #
# se reemplaza cada valor aij de la primera matriz A #
# por el producto de este valor y toda la matriz B #
# dejar [aij * B] #

Z <- direct.prod( x, y) # producto "directo"
print(Z) # 1 columnas de dim (25)
#x11*y
#x21*y...


Z <- direct.prod( V, H ) # producto "directo"
print(Z) # 25 x 25
#v11*H, v12*H...
#v21*H, v22*H...
#v31*H...
#------------------------------------------------- --------------#
#estándares vectoriales#

z <- matrix(c(1,1,1)) # vector diagonal 3D
z
n1 <- entrywise.norm( z, 1 ) # p= 1
n2 <- entrywise.norm( z, 2 ) # norma euclidiana p= 2
n10 <- entrywise.norm(z, 10) # p= 10
n_i <- maximum.norm(z) # max norma p= infinito
print(c(n1,n2,n10,n_i))

z <- matrix(c(0,1,0)) # vector base
z
n1 <- entrywise.norm( z, 1 ) # p= 1
n2 <- entrywise.norm( z, 2 ) # norma euclidiana p= 2
n10 <- entrywise.norm(z, 10) # p= 10
n_i <- maximum.norm(z) # max norma p= infinito
print(c(n1,n2,n10,n_i))

#---------------------------------------------------------------#
# normas de matriz #

# Frobenius (generalización: Hilbert-Schmidt)
print(frobenius.norm(V)) # produce: 696.084

# que es un caso especial de
print(entrywise.norm(V, 2)) # que también da 696.084

# norma máxima
print(maximum.norm(V)) # da 625

# que es un caso especial de
print(entrywise.norm(V, Inf)) # que también da 625

spectral.norm(V) # da 695
#------------------------------------------------- --------------#
# operaciones matriciales #

# transformar una matriz en un vector en el sentido de R
c <- as.vector(H)
print(c) # perdemos las dimensiones y ya no es una matriz
c <- matrix(H)
head(c) # los primeros 5 elementos del vector columna

# transposición de matriz, matrices, dataframes
Tv <- aperm(V, c(2, 1)) # usamos la transposición de dimensiones
Tv <- t(V) # ídem específico de matrices o dataframes
print(Tv)

# diagonal (¡matrices cuadradas!)

d <- diag(V) # diagonal de la matriz
print(d) # línea vectorial en el sentido de R
is.matrix(d) # ¡¡¡esto no es una matriz!!!

d <- matrix(d) #
print(d) # vector de columna como en matemáticas
is.matrix(d) # esta es una matriz y un vector-columna

Dv <- diag(diag(V)) # matriz con diagonal de V
print(Dv) # ¡¡Por fin!!

V0 <- V-Dv # reemplazando el diag de V con 0s
print(V0)

# obtener la matriz identidad
I5 <- diag(1,5) # matriz de 5 x 5 con 1 en la diagonal
print(I5) # en la dimensión 5

# triángulos superior e inferior de la matriz

Vl <- lower.triangle(V) # triángulo inferior incluida la diagonal
print(VI)
Vu <- upper.triangle(V) # triángulo superior incluida la diagonal
print(Vu)
Vbis <- Vl + Vu - Dv # reconstrucción de matriz
print(Vbis)


# calculo de la traza
matrix.trace(V) # suma de términos diagonales

# determinante de una matriz
print(det(V)) # no nulo

#valores propios y vectores propios de una matriz
V_eigen<- eigen(V) # la función de vector propio
print(V_eigen) # resultado a descomponer
V_valp <- diag(V_eigen$values) # matriz diagonal con
print(V_valp) # valores propios clasificados
V_vectp <- V_eigen$vectors # matriz de vectores propios
print(V_vectp)

# inversa de una matriz no singular
Vi <- matrix.inversa(V)
print(Vi)
solve(Vi) # si no usa el paquete matrixcalc
# todavía tenemos una función equivalente
# considerada como menos fiable: solve
# aquí el inverso del inverso es
# la matriz original
print(V %*% Vi) # cerca de la identidad
# vamos a comprobar la desviación t'ipica con Ident nulo en 11 dígitos
print(round(frobenius.norm(V %*% Vi - I5),digits = 11)) # 0

#------------------------------------------------- --------------#
# aplicar una función en una matriz

m <- matrix( 1:12, nrow = 3) # una matriz de prueba
m

# aplicación de una función en cada término
h <- apply(m,c(1,2),exp) # aplicación exponencial en todas partes
print(h) # c(1,2) significa en línea y columna

# suma de cada línea
h <- matrix(apply(m,1,sum)) # vector-columna, suma por fila
print(h) # "1" significa "por línea"

# suma de cada columna
h <- aperm(matrix(apply(m,2,sum)), # línea-vector
           c(2, 1)) # transponemos dim 2 y 1
# con suma por columna
print(h) # "2" significa "por línea"

# aplicación de una función lambda en cada término
# ejemplo, transformar en 1 donde sea >= 3, 0 de lo contrario
h <- apply(m, c(1,2), function(x) {ifelse(x >= 3, 1, 0)})
print(h)

#------------------------------------------------- --------------#
# ¿Esta matriz es singular?
is.singular.matrix(H) # sí
is.singular.matrix( V ) # no










