Arrays multidimensionales en

Acerca de Arrays multidimensionales

Mucho antes, en el concepto Vectores, señalamos que «los arrays pueden ser de tamaño arbitrario (sujetos únicamente a las limitaciones de memoria de tu hardware) y pueden tener tantas dimensiones como se quiera».

Desde entonces, hemos ignorado en gran medida los arrays con más de una dimensión, solo para mantener las cosas sencillas. Esta decisión tendrá más sentido si intentas leer los documentos de referencia de Julia en toda su complejidad.

Sin embargo, los arrays de dimensiones superiores son muy, muy importantes en la computación científica, así que necesitamos entenderlos.

Nomenclatura: Siguiendo siglos de precedente matemático, nos referimos a los arrays unidimensionales como Vectors y a los bidimensionales como Matrices.

Los ejemplos de este documento serán en su mayoría matrices. Trabajar con 3 o más dimensiones es sintácticamente casi idéntico, pero la salida es difícil y confusa de leer (en una pantalla bidimensional).

El type de un array N-dimensional de eltype T es Array{T, N}.

Por comodidad, y para mantener la coherencia con la nomenclatura matemática, Julia define algunos alias de tipo: Vector{T} para Array{T, 1} y Matrix{T} para Array{T, 2}.

# 3-D array
julia> m3 = ones(2, 3, 4);

julia> typeof(m3)
Array{Float64, 3}

# Vector
julia> v = [1, 2]
2-element Vector{Int64}:
 1
 2

julia> typeof(v) == Array{Int64, 1}
true

# Matrix
julia> m
2×3 Matrix{Int64}:
 1  2  3
 4  5  6

julia> typeof(m)
Matrix{Int64} (alias for Array{Int64, 2})

julia> typeof(m) == Array{Int64, 2}
true

Construir arrays

Hemos creado muchos vectores colocando una lista separada por comas entre corchetes. También se pueden usar los puntos y coma como separador.

julia> v = [1, 2, 3]
3-element Vector{Int64}:
 1
 2
 3

julia> w = [1; 2; 3]
3-element Vector{Int64}:
 1
 2
 3

julia> v == w
true

Si usamos el espacio (o el tabulador) como separador, el resultado es diferente.

julia> u = [1 2 3]
1×3 Matrix{Int64}:
 1  2  3

Julia ahora define esto como una matriz 1×3 (en otros contextos, lo llamaríamos row vector).

En general, los espacios unen cosas horizontalmente y los puntos y coma (o los saltos de línea) las unen verticalmente.

La referencia a «cosas» es deliberadamente vaga, ya que Julia intentará trabajar con lo que le des.

julia> [v 2v]
3×2 Matrix{Int64}:
 1  2
 2  4
 3  6

julia> [v; 2v]
6-element Vector{Int64}:
 1
 2
 3
 2
 4
 6

Existen las funciones hcat() y vcat(), que hacen lo mismo, dejando más explícito que se trata de concatenaciones horizontales y verticales. La generalización a más dimensiones es la función cat().

jjulia> hcat(v, 2v)
3×2 Matrix{Int64}:
 1  2
 2  4
 3  6

julia> vcat(v, 2v)
6-element Vector{Int64}:
 1
 2
 3
 2
 4
 6

Escribir matrices explícitas se hace cómodamente en orden por filas, porque encaja con la intuición humana (es más fácil de mirar, para culturas con texto horizontal):

julia> m = [1 2 3; 4 5 6]
2×3 Matrix{Int64}:
 1  2  3
 4  5  6

Sin embargo, ten en cuenta que Julia (como Fortran, R y Matlab, pero a diferencia de C/C++ o NumPy) almacena los arrays de N dimensiones en orden por columnas, y esto puede suponer una enorme diferencia de rendimiento si recorres los elementos en un bucle. ¡Ayuda a la caché de tu CPU para que te ayude!

# put these integers in 2 rows and 3 columns

julia> reshape(collect(1:6), 2, 3)
2×3 Matrix{Int64}:
 1  3  5
 2  4  6

El ejemplo anterior toma los enteros del 1 al 6 y rellena con ellos una matriz 2×3 por columnas.

Hay varias funciones de utilidad para construir tipos comunes de array (uniformes o aleatorios).

julia> zeros(2, 3)  # see also ones()
2×3 Matrix{Float64}:
 0.0  0.0  0.0
 0.0  0.0  0.0

julia> falses(2, 2)  # booleans, see also trues()
2×2 BitMatrix:
 0  0
 0  0

julia> rand(Float32, 2, 3)  # random numbers in the interval [0, 1)
2×3 Matrix{Float32}:
 0.768823  0.169633  0.632565
 0.388451  0.109176  0.850381

Generar valores equiespaciados

Quienes sean aficionados a NumPy sabrán que np.linspace es una función muy utilizada para generar un array de longitud especificada, con extremos especificados y valores equiespaciados. Normalmente se usa como el eje x de una gráfica, o como la variable independiente en modelos lineales.

Julia no tiene un equivalente exacto, pero la muy flexible función range() puede imitarla especificando el límite inferior y el superior, además del argumento de palabra clave length.

Cualquiera que trabaje en local, con acceso al paquete Plots, puede ejecutar este código:

julia> using Plots

julia> x = range(0, 2π; length=100)
0.0:0.06346651825433926:6.283185307179586

julia> typeof(x)
StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}

julia> vals = [x sin.(x) cos.(x)]
100×3 Matrix{Float64}:
 0.0         0.0           1.0
 0.0634665   0.0634239     0.997987
 0.126933    0.126592      0.991955
 0.1904      0.189251      0.981929
(...truncated)

julia> plot(vals[:, 1], vals[:, 2:3])

El StepRangeLen se puede usar como un vector, incluida su conversión a una columna de una matriz.

Consulta también la función equivalente logrange(), para generar valores equiespaciados en un eje logarítmico.

Indexación

Para un array bidimensional, generalmente usamos dos índices en el orden [row, col].

julia> m
2×3 Matrix{Int64}:
 1  2  3
 4  5  6

julia> m[1, 2] # row 1, col 2
2

# Stay within bounds! There is no row 3.
julia> m[3, 1]
ERROR: BoundsError: attempt to access 2×3 Matrix{Int64} at index [3, 1]

julia> m[3]
2

El último ejemplo es quizá sorprendente: un índice único no es un error, y devuelve un solo elemento.

La explicación se remonta al comentario sobre el orden por columnas: Julia recorre la columna 1 hacia abajo, luego la columna 2, hasta encontrar el tercer elemento en memoria.

Ten cuidado: esto tiene algunos usos al escribir bibliotecas de propósito general, pero lo más probable es que resulte confuso.

Surge un problema similar al consultar el tamaño de un array. length() da el número total de elementos y size() da una tupla de ndims() elementos con la longitud de cada dimensión.

julia> m
2×3 Matrix{Int64}:
 1  2  3
 4  5  6

julia> length(m)
6

julia> size(m)  # 2 rows, 3 cols
(2, 3)

julia> ndims(m)  # how many dimensions? Like `m |> size |> length`
2

Indexar con rangos y arrays

Podemos copiar fácilmente una submatriz.

En el ejemplo siguiente, la función reshape() fuerza un array a las dimensiones especificadas rellenándolo por columnas, y luego lo cortamos.

julia> m12 = reshape(collect(1:12), 4, 3)
4×3 Matrix{Int64}:
 1  5   9
 2  6  10
 3  7  11
 4  8  12

julia> m12[2:4, 1:2]
3×2 Matrix{Int64}:
 2  6
 3  7
 4  8

# some rows, all columns
julia> m12[2:4, :]
3×3 Matrix{Int64}:
 2  6  10
 3  7  11
 4  8  12

Un : por sí solo significa «copiar todo en esta dimensión».

Usa vectores para filas o columnas no consecutivas:

julia> m12[[1, 3], :]  # rows 1 and 3
2×3 Matrix{Int64}:
 1  5   9
 3  7  11

Un : por sí solo aplanará cualquier array a un vector, por columnas.

julia> m
2×3 Matrix{Int64}:
 1  2  3
 4  5  6

julia> m[:]
6-element Vector{Int64}:
 1
 4
 2
 5
 3
 6

Aplicar funciones a un array

Ya hemos visto funciones de agregación como sum() y maximum() aplicadas a colecciones unidimensionales, donde operan sobre todos los elementos y devuelven un resultado escalar.

Esto también funciona en dimensiones superiores. Sin embargo, puede que queramos aplicar la función a una sola dimensión; por ejemplo, sumando hacia abajo o hacia los lados para devolver un array con una singleton dimension de tamaño 1.

Para ello, existe un argumento de palabra clave opcional dims.

julia> m
2×3 Matrix{Int64}:
 1  2  3
 4  5  6

julia> sum(m)  # sum everything
21

julia> sum(m; dims=1)  # sum down
1×3 Matrix{Int64}:
 5  7  9

julia> sum(m; dims=2)  # sum across
2×1 Matrix{Int64}:
  6
 15

# other aggregation functions are similar
julia> maximum(m; dims=2)
2×1 Matrix{Int64}:
 3
 6

El valor de dims es la dimensión que pasa a ser singleton, así que dims=1 => 1×3 y dims=2 => 2×1 en los ejemplos anteriores.

Por extensión, los arrays de dimensiones superiores pueden reducir varias dimensiones si dims se establece en un array o un rango (de nuevo, entender la salida puede requerir algo de reflexión).

Escribir funciones que tienen en cuenta las dimensiones

El argumento de palabra clave dims es habitual en funciones integradas como sum(), pero ¿cómo escribimos algo equivalente en nuestro propio código?

Una buena respuesta es usar funciones de orden superior como reduce(), algo que se tratará con cierto detalle en un concepto posterior.

Como alternativa, Julia proporciona varias funciones que nos permiten tratar los arrays de N dimensiones como si fueran vectores de vectores anidados.

Para las matrices, eachrow() y eachcol() son cómodas, pero la función más general es eachslice(), que puede manejar dimensiones arbitrarias.

# m is as in the previous examples

julia> eachrow(m)
2-element RowSlices{Matrix{Int64}, Tuple{Base.OneTo{Int64}}, SubArray{Int64, 1, Matrix{Int64}, Tuple{Int64, Base.Slice{Base.OneTo{Int64}}}, true}}:
 [1, 2, 3]
 [4, 5, 6]

julia> eachcol(m)
3-element ColumnSlices{Matrix{Int64}, Tuple{Base.OneTo{Int64}}, SubArray{Int64, 1, Matrix{Int64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}}:
 [1, 4]
 [2, 5]
 [3, 6]

El tipo es un poco alarmante, pero solo porque esto es una vista del array original, que nos permite trabajar sobre él sin copiarlo. Los arrays de Julia pueden llegar a ocupar terabytes, así que copiarlos puede ser una pesadilla de rendimiento.

Una vista como esta se puede usar para hacer bucles, broadcasting o cualquiera de las demás operaciones que hemos visto en el temario hasta ahora.

Además, las comprensiones de arrays pueden ser potentes y versátiles con un array como entrada. Se mencionaron casos sencillos en el concepto Bucles, pero habrá un análisis más amplio en un concepto posterior.

Broadcasting en múltiples dimensiones

Tratamos el caso unidimensional en el concepto Operaciones con vectores.

julia> v = [1.2, 1.5, 1.7]
3-element Vector{Float64}:
 1.2
 1.5
 1.7

julia> v .- 0.5
3-element Vector{Float64}:
 0.7
 1.0
 1.2

El operador con punto .- trata el 0.5 singleton como equivalente a [0.5, 0.5, 0.5] al hacerle broadcasting para que coincida en dimensiones, y luego realiza la resta elemento a elemento.

Extender esto a dimensiones superiores no es más que lo mismo otra vez.

Por ejemplo, al hacer broadcasting de la multiplicación de un vector fila 2x1 [1.0 1.5 2.0] con la matriz 2x3 [1 2 3; 4 5 6], equivale a hacer broadcasting de [1.0 1.5 2.0] a [1.0 1.5 2.0; 1.0 1.5 2.0] y luego realizar la multiplicación elemento a elemento.

julia> m
2×3 Matrix{Int64}:
 1  2  3
 4  5  6

julia> m .* [1.0 1.5 2.0]  # broadcast a 1x3 row vector down
2×3 Matrix{Float64}:
 1.0  3.0   6.0
 4.0  7.5  12.0

julia> m .* [1.0 1.5 2.0; 1.0 1.5 2.0]  # 2x3 elementwise multiplication
2×3 Matrix{Float64}:
 1.0  3.0   6.0
 4.0  7.5  12.0

julia> m ./ [1, 2]  # broadcast a 2x1 column vector across
2×3 Matrix{Float64}:
 1.0  2.0  3.0
 2.0  2.5  3.0

julia> m ./ [1 1 1; 2 2 2]  # 2x3 elementwise division
2×3 Matrix{Float64}:
 1.0  2.0  3.0
 2.0  2.5  3.0

Nota: Son las dimensiones no singleton las que deben coincidir en tamaño (por ejemplo, 2x3 Matrix .* 2x1 Matrix) para poder hacer broadcasting.

Además, una función se puede aplicar con broadcasting a cada elemento de una Matrix exactamente igual que como se hace con un Vector.

julia> m
2×3 Matrix{Int64}:
 1  2  3
 4  5  6

julia> (x -> x^2).(m)  # broadcast a function to all elements
2×3 Matrix{Int64}:
  1   4   9
 16  25  36

Puede que al principio te parezca que el margen de confusión del programador crece exponencialmente con el número de dimensiones, pero la práctica ayuda mucho. Además, la familiaridad con los arrays de NumPy se traslada bastante bien a Julia, a pesar de la sintaxis diferente.

Lecturas adicionales

Muchas operaciones habituales con vectores y matrices entran dentro del ámbito general del álgebra lineal.

Actualmente hay dos conceptos adicionales en el sitio web de Exercism, que son solo documentos hasta que podamos diseñar ejercicios que los acompañen.

1. Fundamentos de álgebra lineal trata principalmente varios tipos de multiplicación.
2. Resolución de ecuaciones lineales trata los sistemas de ecuaciones lineales y cómo resolverlos de forma eficiente.