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

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

7 голосов
/ 29 мая 2011

Я пытаюсь реализовать функцию кумулятивной суммы, используя 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 ]

5 голосов
/ 31 мая 2011

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

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

3 голосов
/ 30 мая 2011

На самом деле, я думаю, что главная проблема в том, что у Repa нет примитива сканирования.(Однако очень похожая библиотека Accelerate делает.) Существует два варианта сканирования: сканирование префикса и сканирование суффикса.При заданном одномерном массиве

[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.

...