Массивы

Фикус разработан специально для числовой обработки, поэтому в нём обеспечена хорошая поддержка многомерных массивов.

Многомерный массив – встроенный тип в Фикус, обозначаемый как 't [, ...], где количество запятых соответствует размерности массива минус единица. Многомерные массивы хранятся в памяти подряд в порядке строк (аналогично C/C++ и в отличие от Fortran и MATLAB, где используется порядок столбцов), и они не являются массивами массивов.

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

// Создадим массив 480x640 восьмибитных чисел,
// пригодный для представления серого изображения.
// Инициализируем элементы нулями.
val image = array((480, 640), 0u8)

// Построим другое изображение, выполнив арифметическую операцию
val lighter_image = image .+ 100u8

val phi = M_PI/3
val (a, b) = (cos(phi), sin(phi))
// Создадим матрицу поворота 2x2, явно указав все её элементы
val R = [a, -b; b, a]
val RR = R*R // произведём умножение матриц
val shift = (100., 10.)
// Составим аффинную матрицу 2x3:
// левый блок 2x2 будет инициализирован R,
// правый блок 2x1 — смещающим вектором
val affine = [\R, $$shift.0; shift.1]]

// Используя comprehension
val gradient = [for y <- 0:256 for x <- 0:256
                  { sat_uint8((x, y, (x+y)/2)) }]

Создав массив, мы можем легко читать и изменять его:

fun sumpixels(img: uint8 [,]): int =
    // Итерация по массиву с помощью цикла for
    fold sum=0 for x <- img {sum + x}

// Выделим интересующую область (ROI);
// данные не копируются, а совместно используются
val roi = image[50:100, 300:400]
// Подсчитаем сумму пикселов в ROI
val sum_roi = sumpixels(roi)

// Изменим ROI; родительское изображение изменится тоже
val (m, n) = size(roi)
for y <- 0:m
  for x <- 0:n {
     roi[y, x] ^= 255u8
  }

Отметим, что элементы roi могут быть изменены, даже если сам roi является значением. Это связано с тем, что массивы — изменяемые структуры.

Массивы индексируются с помощью нотации arr[idx_or_range1, idx_or_range2, ...], где каждый idx_or_rangej — это либо целое число (индекс), либо диапазон [start]:[end]:[step], где любой из компонентов start, end или step может отсутствовать, либо двоеточие :, означающее “взять все элементы вдоль этой оси (всё поле, весь столбец и т.д.)”:

// создадим ненормализованную матрицу Адамара
fun Hadamard(n: int, one: 't): 't [,] {
    assert(n & (n-1) == 0)
    if n == 1 { [for i <- 0:1 for j <- 0:1 {one}] }
    else { val h = Hadamard(n/2, one); [\h, \h; \h, $-h)] }
}

val mtx = Hadamard(8, 1.)

// возьмём пятую строку матрицы, умножим все элементы на 2
// и прибавим полученный продукт к первой строке той же матрицы
mtx[1, :] += mtx[5, :]*2

Такие матричные преобразования возможны, поскольку операция выделения подмассива не копирует никаких данных. Вместо этого формируется новая заголовочная структура, инициализируется её размер, тип и размерность, а затем она указывает на оригинальную память массива с некоторым смещением. Счётчик ссылок оригинальной памяти увеличивается. Есть два исключения из этого правила:

  • если в диапазоне подмассива использован отрицательный шаг или ненулевой шаг в последней размерности, то подмассив нужно скопировать в новое место.
  • выполнена пперация flatten (arr[:]), а массив не непрерывен. В этом случае данные также подлежат копированию.

Когда подмассив больше не нужен, он разрушается, счётчик ссылок уменьшается.

Обращение к элементам массива

Когда в программе встречается обращение к элементу массива A[i, j, k, ...], компилятор Фикус генерирует код, аналогичный следующему псевдокоду:

if (i < 0 || i >= A.size[0] ||
   j < 0 || j >= A.size[1] ||
   k < 0 || k >= A.size[2] ||
   ...) {
    __throw__(OutOfRangeError);
}
elemtype* ptr = (elemtype*)(A.data + A.step[0]*i +
                   A.step[1]*j + A.step[2]*k + ... );

То есть каждая операция обращения к элементу проверяет, находится ли индекс элемента в границах массива. Если индекс выходит за пределы, возникает исключение OutOfRangeError. Хотя это повышает безопасность и помогает быстро выявить ошибки, это может замедлить доступ к массивам. В случае выбора диапазонов вместо конкретных индексов границы также проверяются, но накладные расходы заметно ниже.

Для повышения производительности и сохранения безопасности компилятор Фикус применяет особые оптимизации, уменьшая количество проверок границ:

  • если внутри цикла for обращение к массиву выполняется безусловным образом (то есть не внутри вложенных операторов if, match, try и т.д.)
  • и если индекс ij является линейной комбинацией индекса цикла и некоторых инвариантов цикла: ij == loop_idx*loop_inv + another_loop_inv

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

Рассмотрим пример реализации функции умножения матриц, приведённой ранее:

fun matmul(A: double [,], B: double [,])
{
    val (ma, na) = size(A), (mb, nb) = size(B)
    assert(na == mb)
    [for i <- 0:ma for j <- 0:nb {
        fold s = 0. for k <- 0:na {s + A[i,k]*B[k,j]}
     }]
}

Генерируемый код (псевдокод на языке Фикус) выглядел бы так:

fun matmul(A: double [,], B: double [,])
{
    val (ma, na) = size(A), (mb, nb) = size(B)
    assert(na == mb)
    // Проверяем границы массива для индексов i и j
    check_range(0, ma, 1, 0, A, 0)
    check_range(0, nb, 1, 0, B, 1)
    [for i <- 0:ma for j <- 0:nb {
        // Аналогично проверяем индекс k
        check_range(0, na, 1, 0, A, 1)
        check_range(0, na, 1, 0, B, 0)
        // <.> обозначает быстрый доступ без проверки границ,
        // это не реальный синтаксис Фикус
        fold s = 0. for k <- 0:na {s + A[<i>,<k>]*B[<k>,<j>]}
     }]
}

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

Более быстрым методом доступа к элементам массива является их запрос прямо в заголовке цикла for, например:

// здесь x извлекается с помощью *(указатель_на_субмассив + индекс_цикла) на уровне C-кода
fold sum = 0. for x <- arr {s + x}

Преимущества такого подхода:

  • проверка границ не требуется, так как известно, что мы проходим по всему массиву.
  • адрес элемента не нужно вычислять заново, достаточно увеличить счетчик цикла и простым действием обратиться к элементу. Особенно полезно это в случае прохода по двумерной матрице, когда перед внутренним циклом устанавливается указатель на текущую строку матрицы, а затем он используется как базовый указатель, добавляясь к счетчику цикла для получения текущего элемента матрицы. Такой подход необходим, поскольку матрица не обязана быть смежной (может быть подматрицей большей матрицы).

Следовательно, если нам нужны как элемент, так и его индекс, эффективнее запрашивать их оба в заголовке цикла for, чем обращаться к элементу через оператор [].

Унарный оператор .-

В Python удобно получать последний элемент строки или списка с помощью отрицательных индексов, например, python_str[-1] возвращает последний символ строки.

В Фикус, как мы недавно узнали, подобная операция приведет к исключению OutOfRangeError, что зачастую и требуется. Однако в Фикус предусмотрен оператор .-, предназначенный для доступа к последнему, предпоследнему и т.д. элементам. То есть вместо использования выражений вида N1-i1, N2-i2, где Nj — размер соответствующего измерения массива, а ij — индекс, можно просто поставить точку вместо N1, N2 и т.д.:

  • A[0, .-1] — доступ к верхнему правому углу матрицы
  • A[.-1, .-1] — доступ к нижнему правому углу матрицы
  • A[:.-1] — взятие всех элементов одномерного массива, кроме последнего

Базовые операции с массивами

Ниже приведены некоторые базовые операции обработки массивов, собранные для удобства:

  1. Получить размер или форму массива: size(A). Узнать размер i-го измерения: __intrin_size__(A, i)
  2. Извлечь строку матрицы: A[i, :]
  3. Извлечь столбец матрицы: A[:, i]
  4. Извлечь четные строки: A[::2, :]. Это можно сделать без копирования элементов матрицы.
  5. Отразить матрицу вокруг горизонтальной оси, то есть перевернуть порядок строк: A[::-1,:]
  6. Отразить матрицу вокруг вертикальной оси, то есть перевернуть каждую строку: A[:,::-1]
  7. Отразить матрицу симметрично относительно центра: A[::-1,::-1]
  8. Сплющить массив, превратив его в одномерный: A[:]
  9. Применить поэлементную бинарную операцию к каждому соответствующему элементу пары матриц: A op B, где op — один из операторов .+, .-, .*, ./, .%, .**, &, |, ^.
  10. Горизонтально соединить несколько массивов: [\A1, \A2, ...]
  11. Вертикально соединить несколько массивов: [\A1; \A2; ...]
  12. Преобразовать список или строку в массив: [\source]
  13. Преобразовать массив обратно в список или строку: list(arr), string(arr). Во втором случае массив должен иметь тип char [], иначе функция string(arr) создаст текстовое представление массива.
  14. Объединить несколько массивов в многоканальный массив: [for x1 <- A1, x2 <- A2, ... {(x1, x2, ...)}] TODO
  15. Разделить многоканальный массив на несколько одноканальных массивов: [@unzip for x <- A {x}]
  16. Транспонировать матрицу с помощью апострофа после оператора: A'
  17. Умножить две матрицы: A*B
  18. Инвертировать матрицу (текущие реализации используют QR-алгоритм): A\1
  19. Решить систему уравнений вида Ax=b: A\b
  20. Вычислить детерминант или след матрицы: det(A), trace(A)
  21. Отсортировать элементы одномерного массива согласно пользовательскому критерию упорядоченности: sort(A, less_than)
  22. Вычислить сумму, среднее значение, нормы L∞, L₁ или L₂ (также известную как норма Фробениуса) массива: sum(A), mean(A), normInf(A), normL1(A), normL2(A). Для вариантов функций norm* с двумя параметрами вычисляется норма разницы между двумя массивами.

Экспраполяция границ

В обработке сигналов, изображений и глубоком обучении нередко возникает потребность обращаться к массиву за пределами его границ. Например, при применении фильтра к изображению, скажем, Гауссова размытия, для каждого пикселя, включая краевые и угловые, нужно вычислить взвешенную сумму окрестности размеров 2R+1 × 2R+1, где R — радиус фильтра. Возможные подходы к решению этой задачи:

  1. создавать на выходе массив меньшего размера, обрабатывая только те элементы входного массива, которые достаточно удалены от границы.
  2. Копировать массив в центр нового увеличенного массива, заполняя оставшиеся элементы новым массивом с помощью определенной формулы (например, распространение крайних элементов исходного массива на границу нового массива):
// Исходный массив:
a b c d
e f g h

// Расширенный массив с радиусом 2:
a a a b c d d d
a a a b c d d d
a a a b c d d d
e e e f g h h h
e e e f g h h h
e e e f g h h h

  1. Изменить алгоритм так, чтобы каждый доступ к элементам входного массива заменялся более сложной формулой, например, вместо A[i, j] использовать нечто вроде A[min(max(i, 0), M), min(max(j, 0), N)].

  2. Разделить цикл обработки массива на несколько частей (минимум на две), где основная “быстрая” часть обрабатывает внутреннюю область массива без учета границ (метод 1 в этом списке), а затем обработать граничные элементы с помощью метода 3.

Выбор конкретной стратегии зависит от приложения и ряда обстоятельств, таких как размер массива, радиус экстраполяции, сложность фильтра и т.д. Любой из четырех подходов может оказаться полезным.

Фикус пытается облегчить реализацию методов №3 и №4, вводя специальные разновидности операций доступа к элементам массива:

  • A.clip[i1, i2, ...] — доступ к элементам массива с отсечением. Каждый индекс ij ограничивается с помощью формулы min(max(ij,0),Nj), где Nj == __intrin_size__(A,j).
  • A.zero[i1, i2, ...] — доступ к элементам массива с дополнением нулями. Если любой из индексов ij выходит за пределы диапазона [0,Nj), возвращается ноль (того же типа, что и элементы массива; для многоканального массива возвращается нулевой кортеж).
  • A.wrap[i1, i2, ...] — доступ к элементам массива с цикличным повторением. Индекс ij обрабатывается по формуле ij mod Nj, где mod — операция взятия остатка от деления, возвращающая неотрицательные значения. Например, для одномерного массива A.wrap[N1] ~ A[0], A.wrap[N1+1] ~ A[1], A.wrap[-1] ~ A[N1-1].

Операции экстраполяции границ имеют ограничения и особенности:

  • их нельзя использовать с диапазонами; все индексы должны быть скалярными значениями.
  • эти операции применимы только к массивам, чьи элементы являются скалярными значениями (числа, bool, char) или кортежами скалярных значений.
  • ни одна из этих операций не генерирует исключение. .clip и .wrap всегда возвращают один из элементов массива, если массив не пуст. Если массив пустой, возвращается нулевой элемент (аналогично .zero).
  • данные операции могут применяться не только к массивам, но и к объектам типа vector и string.

Вот наивная, но рабочая реализация фильтра размытия 3x3 с использованием операции .clip:

fun blur3x3(img: uint8 [,])
{
   val (h, w) = size(img)
   [for y <- 0:h for x <- 0:w {
	  sat_uint8((img.clip[y-1, x-1] + img.clip[y-1, x] +
				 img.clip[y-1, x+1]+ img.clip[y, x-1] +
				 img[y, x] + img.clip[y, x+1]+
				 img.clip[y+1, x-1] + img.clip[y+1, x] +
				 img.clip[y+1, x+1])/9)
	 }]
}

В качестве упражнения попробуйте реализовать вариант этой функции, который:

  • создает выходной массив размера h×w с помощью array((h, w), 0u8)
  • вычисляет центральную часть результата, используя цикл по диапазону 1≤y<h-1, 1≤x<w-1. В этом цикле компилятор Фикус должен исключить проверку границ, так как мы гарантированно остаемся в пределах исходного изображения.
  • рассчитывает границы результата с помощью операции .clip. Вероятно, эту часть можно оформить в отдельную функцию, так как она должна вызываться четыре раза для каждой границы.