Как взять срез массива с Repa в диапазоне

Я пытаюсь реализовать функцию кумулятивного суммирования с использованием Repa для вычисления интегральных изображений. Моя текущая реализация выглядит следующим образом:

cumsum :: (Elt a, Num a) => Array DIM2 a -> Array DIM2 a
cumsum array = traverse array id cumsum' 
    where
        elementSlice inner outer = slice array (inner :. (0 :: Int))
        cumsum' f (inner :. outer) = Repa.sumAll $ elementSlice inner outer

Проблема в функции elementSlice. В matlab или скажем numpy это можно указать как массив [inner,0: external]. Так что я ищу что-то вроде:

slice array (inner :. (Range 0 outer))

Однако, похоже, что нельзя разрешать указывать срезы в диапазонах, которые в настоящее время находятся в Repa. Я рассмотрел использование разделов, как обсуждалось в Efficient Parallel Stencil Convolution в Haskell, но это кажется довольно тяжелым подходом, если они будут меняться при каждой итерации. Я также рассмотрел маскирование среза (умножение на двоичный вектор) - но опять-таки казалось, что он будет работать очень плохо на больших матрицах, так как я буду выделять вектор маски для каждой точки в матрице...

Мой вопрос - кто-нибудь знает, есть ли планы добавить поддержку для нарезки в диапазоне в Repa? Или есть эффективный способ, которым я уже могу атаковать эту проблему, может быть, с другим подходом?

Источник

2 ответа

Решение

На самом деле, я думаю, что главная проблема в том, что у Repa нет примитива сканирования. (Однако, очень похожая библиотека Accelerate делает.) Есть два варианта сканирования: сканирование префикса и сканирование суффикса. Учитывая 1D массив

[a_1, ..., a_n]

сканирование префикса возвращает

[0, a_0, a_0 + a_1, ..., a_0 + ... + a_{n-1} ]

в то время как сканирование суффикса производит

[a_0, a_0 + a_1, ..., a_0 + a_1 + ... + a_n ]

Я предполагаю, что это то, что вы собираетесь с вашей накопленной суммой (cumsum) функция.

Сканирование префиксов и суффиксов вполне естественно обобщает многомерные массивы и имеет эффективную реализацию, основанную на сокращении дерева. Относительно старая статья на эту тему - "Сканирование примитивов для векторных компьютеров". Кроме того, Конал Эллиотт недавно написал несколько постов в блоге о создании эффективных параллельных сканирований в Haskell.

Интегральное изображение (на двумерном массиве) можно рассчитать, выполнив два сканирования, одно по горизонтали и одно по вертикали. В отсутствие примитива сканирования я реализовал один, очень неэффективно.

horizScan :: Array DIM2 Int -> Array DIM2 Int
horizScan arr = foldl addIt arr [0 .. n - 1]
  where 
    addIt :: Array DIM2 Int -> Int -> Array DIM2 Int
    addIt accum i = accum +^ vs
       where 
         vs = toAdd (i+1) n (slice arr (Z:.All:.i))
    (Z:.m:.n) = arrayExtent arr

--
-- Given an @i@ and a length @len@ and a 1D array @arr@ 
-- (of length n) produces a 2D array of dimensions n X len.
-- The columns are filled with zeroes up to row @i@.
-- Subsequently they are filled with columns equal to the 
-- values in @arr.
--
toAdd :: Int -> Int -> Array DIM1 Int -> Array DIM2 Int
toAdd i len arr = traverse arr (\sh -> sh:.(len::Int)) 
               (\_ (Z:.n:.m) -> if m >= i then arr ! (Z:.n) else 0)

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

vertScan :: Array DIM2 Int -> Array DIM2 Int
vertScan = transpose . horizScan . transpose

integralImage = horizScan . vertScan

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

Источник

Извлечение поддиапазона - это манипулирование индексным пространством, которое достаточно легко выразить с помощью fromFunction, хотя, вероятно, нам следует добавить для него более приятную оболочку в API.

let arr = fromList (Z :. (5 :: Int)) [1, 2, 3, 4, 5 :: Int] 
in  fromFunction (Z :. 3) (\(Z :. ix) -> arr ! (Z :. ix + 1))

> [2,3,4]

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

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

Источник
Другие вопросы по тегам