Массивы и broadcast

1.8. Массивы и broadcast#

В Julia для массивов существует абстрактный тип AbstractArray{T,N}, где T тип элемента массива (Int64, String, Complex{Int64}, Number, Any…), а N размерность. При написании собственных типов с поведением массива в объявлении следует указать, что ваш тип является подтипом AbstractArray.

Встроенным типом для массивом является Array{T,N}. Это изменяемый массив динамического размера с column-major хранением элементов.

Примечание

Поскольку в Julia массивы хранятся постолбично, то для быстродействия вложенный цикл должен проходится по строке, а внешний по столбцу. Вообще, если есть массив A[i, j, k], то цикл должен выглядить так:

for k in eachindex(A, 3)
    for j in eachindex(A, 2)
        for i in eachindex(A, 1)
            ...
        end
    end
end

Существуют два часто используемых синонима

  • Vector{T} == Array{T,1} для одномерных вектор-столбцов;

  • Matrix{T} == Array{T,2} для матриц и вектор-строк.

Литералом Array{T,N} являются квадратные скобки [], при этом можно явно указать тип T[] (Int[1, 2], Float64[1, 2]). Если тип элементов массива не указывать, то Julia его вычислит.

Для создания массива через литерал используется вертикальная и горизонтальная конкатенация.

Для вертикальной конкатенации используется ; или новая строка.

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

julia> [[1, 2, 3]; [7, 8]]
5-element Vector{Int64}:
 1
 2
 3
 7
 8

Для горизонтальной конкатенации используется пробел или табуляция. Запятая , означает отсутствие конкатенации.

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

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

Если конкатенация не нужна, используйтся запятая ,. В примере ниже показано создание вектора векторов.

julia> [[1, 2, 3], [7, 8]]
2-element Vector{Vector{Int64}}:
 [1, 2, 3]
 [7, 8]

На первый взгляд синтаксис может показаться запутанным. Однако, с ним легко создавать блочные матрицы и тензоры.

Демонстрация 1.3 (Создание блочной матрицы)

Скажем, задана блочная матрица

\[\begin{split}\mathcal{M} = \left[\begin{array}{ccc|ccc} 1 & 2 & 3 & 1 & 0 & 0 \\ 4 & 5 & 6 & 0 & 1 & 0 \\ 7 & 8 & 9 & 0 & 0 & 1 \\ \hline 0 & 0 & 0 & 10 & 11 & 12 \\ \end{array}\right] = \left[\begin{array}{c|c} A & I \\ \hline 0^\top & b^\top \end{array}\right]\end{split}\]

Соответствующий код будет

julia> A = [1 2 3; 4 5 6; 7 8 9]   # или reshape(collect(1:9), 3, 3)'
3×3 Matrix{Int64}:
 1  2  3
 4  5  6
 7  8  9

julia> I = [1 0 0; 0 1 0; 0 0 1]
3×3 Matrix{Int64}:
 1  0  0
 0  1  0
 0  0  1

julia> zero = [0, 0, 0]  # или zeros(Int, 3)
3-element Vector{Int64}
 0
 0
 0

julia> b = [10, 11, 12]  # или collect(10:12)
3-element Vector{Int64}:
 10
 11
 12

julia> M = [A I; zero' b']
4×6 Matrix{Int64}:
 1  2  3   1   0   0
 4  5  6   0   1   0
 7  8  9   0   0   1
 0  0  0  10  11  12

Оператор ' (одинарная кавычка) выполняет транспонирование.

1.8.1. Индексирование#

В Julia встроено два вида индексов.

Синтаксис первого A[i], это линейный индекс LinearIndex. Он позволяет получить i-ый по счёту в памяти элемент массива. Например, для вектора A[4] это \(A_4\), а для матрицы размера \(2 \times 3\) это \(A_{2, 2}\).

Синтаксис второго A[i, j, k,..], это декартов индекс CartesianIndex. Данный индекс учитывает размерность массива. Например, для матрицы элемент A[i, j] это элемент \(A_{i,j}\) (\(i\)-ая строка, \(j\)-ый столбец).

Индексом может быть и коллекция элементов, что позволяет обращаться к целым строкам, столцам, их пересечению, или просто некоторому набору элементов. Часто для этого используется range operator start:[step:]stop.

Ниже показано несколько примеров для матрицы.

julia> A = [1 2 3; 4 5 6; 7 8 9]
3×3 Matrix{Int64}:
 1  2  3
 4  5  6
 7  8  9

julia> A[6]     # линейный индекс
8

julia> A[3, 2]  # декартов индекс
8

julia> A[end, end]
9

Строка.

julia> A[1, :]  # A[строка 1, все столбцы]
3-element Vector{Int64}:
 1
 2
 3

Столбец.

julia> A[:, 2]  # A[все строки, столбец 2]
3-element Vector{Int64}:
 2
 5
 8

«Пересечение» строк 1 и 2 со столбцами 2 и 3.

julia> A[1:2, 2:3]
2×2 Matrix{Int64}:
 2  3
 5  6

Угольные элементы.

julia> A[[1, end], [1, end]]
2×2 Matrix{Int64}:
 1  3
 7  9

Главная диагональ.

julia> A[1:4:end]
3-element Vector{Int64}:
 1
 5
 9

Побочная диагональ.

julia> A[3:2:7]
3-element Vector{Int64}:
 7
 5
 3

1.8.2. Присваивание значений#

Массив Array{T,N} изменяемый тип, можно перезаписывать его элементы. Это делается с помощью оператора присваивания =, при этом по обе стороны должны находиться объекты одинаковых размерностей.

julia> A = [1 2 3; 4 5 6; 7 8 9];

julia> A[:, 1] = [7, 4, 1];

julia> A
3×3 Matrix{Int64}:
 7  2  3
 4  5  6
 1  8  9

julia> A[1, :] = [-1, 0, 1];

julia> A
3×3 Matrix{Int64}:
 -1  0  1
  4  5  6
  1  8  9

julia> A[1, 2] = 10;

julia> A
3×3 Matrix{Int64}:
 -1  10  1
  4   5  6
  1   8  9

При присваивании можно использовать векторизацию.

julia> A[1, :] .= 11;

julia> A
3×3 Matrix{Int64}:
 11  11  11
  4   5   6
  1   8   9

1.8.3. Вспомогательные функции#

Интроспекция массива

  • length(A): количество элементов;

  • ndims(A): размерность;

  • size(A): кортеж из размеров каждой размерности (число_строк, число_столбцов);

  • size(A, n): размер размерности n;

  • eachindex(A): итератор из индексов массива.

Создание массивов

  • collect(c): создание массива из коллекции;

  • zeros([T=Float64,] dims...): массив из нулей типа T размерности dims...;

  • ones([T=Float64,] dims...): аналогично zeros, но массив из единиц;

  • copy(A): нерекурсивная копия массива;

  • similar(A): неинициализированный массив того же типа и размера, что и A;

  • empty(A): пустой массив того же вида, что и A;

  • rand([T=Float64,] dims...): массив из равномерно-распределённых случайных чисел;

  • fill(x, dims...): массив размера dims... из x;

  • map(f, c): создание массива на основе применения функции f к каждому элементу коллекции c.

Другие операции

  • push!(A, x...): добавление в конец A элемента x;

  • pop!(A): извлечение элемента с конца A.

Для чтения и записи в файл массивов в форматах с сепаратором (.csv, .tsv,..) есть встроенная библиотека DelimitedFiles и сторонний пакет CSV.jl.

1.8.4. Векторизация функций#

Совет

Соответствующий раздел мануала [url].

В Julia любая функция может быть векторизована. Т.е. вы можете создать функцию для скаляров и автоматически получаете её векторизованный аналог. Векторизация осуществляется через механизм broadcasting. При этом достаточно в конце имени функции дописать точку . (dot-syntax).

julia> x = [0, π/6, π/3];

julia> sin(x)
ERROR: MethodError: no method matching sin(::Vector{Float64})
  ...

julia> sin.(x)
3-element Vector{Float64}:
 0.0
 0.49999999999999994
 0.8660254037844386

Механизм векторизации «насыщен» (fused). Для программиста это означает, что сложный вызов, например, sin.(cos.(x)) не медленнее цикла, в котором к каждому элементу вектора x применили sin(cos(x[i])).

Автоматически векторизуются и пользовательские функции над стандартными типами.

julia> f(x) = x^2 + 1
f (generic function with 1 method)

julia> f.([1, 2, 3])
3-element Vector{Int64}:
  2
  5
 10

julia> g(x, y) = f(x) + f(y)
g (generic function with 1 method)

julia> g.([1, 2], [3, 4])
2-element Vector{Int64}:
 12
 22

julia> g.([1, 2], 3)
2-element Vector{Int64}:
 12
 15

julia> parse.(Float64, ["1", "2.5", "-1.2e-9"])
3-element Vector{Float64}:
  1.0
  2.5
 -1.2e-9

Для собственных типов нужно объявлять правила векторизации. Однако, если ваш тип должен вести себя как скаляр, можно воспользоваться Ref(x).

Пример 1.4

Ниже определяется модель идеального газа и его давление. Затем, с помощью Ref и векторизации вычисляется давление кислорода при фиксированной плотности для ряда температур.

julia> struct IdealGasComponent
           name::String
           molarmass::Float64
       end

julia> pressure(c::IdealGasComponent, ρ, T) = ρ * 8.314 * T / c.molarmass
pressure (generic function with 1 method)

julia> oxygen = IdealGasComponent("oxygen", 0.016)
IdealGasComponent("oxygen", 0.016)

julia> pressure(oxygen, 500, 300)
7.794375e7

julia> pressure.(Ref(oxygen), 500, [300, 350, 400])
3-element Vector{Float64}:
 7.794375e7
 9.0934375e7
 1.03925e8

Операторы тоже функции. Так, broadcast порождает, например, поэлементные операции над векторами и матрицами (.+, .*, …).

julia> 2 * sin.(x) .* cos.(x) == sin.(2 * x)
true
Наблюдение (Полиморфизм)

В этом примере несколько поэлементых видов умножений.

  • 2 * x: умножение скаляра на вектор, оно и без броадкаста действует поэлементно;

  • x .* y: поэлементное умножение векторов.

Если вам надоедает ставить точки везде, воспользуйтесь макросом @.

julia> x = [0, π/6, π/3];

julia> @. 2 * sin(x) * cos(x)
3-element Vector{Float64}:
 0.0
 0.8660254037844386
 0.8660254037844388

Ещё одно важное применение броадкаста: in-place присваивание. Оно делает вашу программу менее прожорливой по памяти, а сборщик мусора не будет добавлять времени работы.

julia> x = rand(100); y = rand(100);

julia> @allocated z = x + y;    # первый вызов включает компиляцию

julia> @allocated z = x + y
896

julia> z = similar(x);          # выделили память под `z` один раз

julia> @allocated z .= x .+ y;  # первый вызов включает компиляцию

julia> @allocated z .= x .+ y
64
Демонстрация 1.5 ((Ещё один) бенчмарк)

На аллокации требуется время, а ещё время требуется для работы сборщика мусора. На данных небольшого размера время сборщика мусора существенно по сравнению с временем исполнения.

Так, сложение двух векторов размера 100 с аллокацией требует в тесте в среднем 85 нс, а с использованием in-place операций 22 нс. Причём разницу, судя по бенчмарку, занимает сборщик мусора (garbage collector, GC).

julia> using BenchmarkTools

julia> x = rand(100); y = rand(100); z = similar(x);

julia> @benchmark $z = $x + $y
BenchmarkTools.Trial: 10000 samples with 973 evaluations.
 Range (min … max):  75.194 ns … 619.990 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     78.805 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   85.932 ns ±  41.564 ns  ┊ GC (mean ± σ):  5.88% ± 9.87%

  █▄▁                                                          ▁
  █████▇█▇▅▅▅▅▄█▇▄▃▃▅▄▁▄▃▃▁▁▄▄▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▄▁▁▃▁▃▃▁▁▄▅▆▅▇▆▇▆ █
  75.2 ns       Histogram: log(frequency) by time       365 ns <

 Memory estimate: 896 bytes, allocs estimate: 1.

julia> @benchmark $z .= $x .+ $y
BenchmarkTools.Trial: 10000 samples with 997 evaluations.
 Range (min … max):  20.640 ns … 398.402 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     20.930 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   22.264 ns ±   8.109 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▃▁▁▁                                                        ▁
  █████▇▆▆▇▇▆▆▆▆▆▆▆▆▅▅▅▆▆▆▆▆▇▇▇▇▇▇▆▇▇▆▇▇▇▆▅▆▆▆▆▄▄▅▅▄▄▄▄▄▅▃▁▄▁▄ █
  20.6 ns       Histogram: log(frequency) by time        45 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

Повторим тест для больших векторов, возьмём их размером 1000000.

julia> x = rand(1_000_000); y = rand(1_000_000); z = similar(x);

julia> @benchmark $z = $x + $y
BenchmarkTools.Trial: 2902 samples with 1 evaluation.
 Range (min … max):  901.201 μs …   6.733 ms  ┊ GC (min … max):  0.00% … 46.36%
 Time  (median):       1.275 ms               ┊ GC (median):     0.00%
 Time  (mean ± σ):     1.713 ms ± 976.100 μs  ┊ GC (mean ± σ):  26.05% ± 25.78%

        ██▄▁                                            ▁▄▄▂▁   ▁
  ▃▁▁▁▁▅████▇▄▇▅▃▃▁▃▁▁▄▃▁▁▁▃▁▃▃▁▃▁▁▁▃▁▁▁▁▁▁▁▁▁▃▃▁▁▁▁▁▁▁▅██████▇ █
  901 μs        Histogram: log(frequency) by time        470 ms <

 Memory estimate: 7.63 MiB, allocs estimate: 2.

julia> @benchmark $z .= $x .+ $y
BenchmarkTools.Trial: 5360 samples with 1 evaluation.
 Range (min … max):  789.767 μs …   5.811 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     832.923 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   921.733 μs ± 295.907 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▅█▆▃▃▂▁▁▁▁▁▁▁▁                                                ▁
  ██████████████████▇█▇▇▆▆▇▆▅▆▅▆▆▅▆▅▅▅▄▄▄▄▄▄▄▃▄▅▅▄▅▅▁▄▃▄▅▄▄▄▄▄▅ █
  790 μs        Histogram: log(frequency) by time       2.42 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

На векторах размера 100 разница среднего времени исполнения составляла \(\approx 4\) раза, а на больших векторах уменьшилась до \(\approx 2\) раз.

Содержание