Haskell repa - how to reduce array and return index?

Question

In GNU Octave this code -

[e, ix] = min(X);

will return minimum element and it's location. How do you this in repa for arbitrary binary function?

This is what I came up with:

min x = z $ foldl' f (e,0,0) es
  where
    (e:es) = toList x
    f (a,ix,r) b = let ix' = ix+1 in if a < b then (a,ix',r) else (b,ix',ix')
    z (a,ix,r) = (a,r)

In above example we convert repa 1D matrix to list and use foldl' (from Data.List) with two accumulators - one for counting iterations (ix) and other to save position of min element (r). But the whole point of using repa is to use arrays, not lists!

In repa there are two folds for Array type (foldS and foldP) - but they can only take function of type (a -> a -> a) - meaning, I cannot pass tuple with accumulators to it. There is also traverse, which can, in principle, reduce 1D array to a scalar array:

min x = traverse x to0D min
  where
    to0D (Z:.i) = Z
    min f (Z) = ??? -- how to get elements for comparison?

The first thing that comes to mind is

[f (Z:.i) | i <- [1..n]], where n = (\(Z:.i) -> i) $ extent x

But this will also convert array to list, instead of doing computation on the array.

Not really familiar with repa, but can't you just map each item to a tuple before using `foldP?`. For example, you can describe a subarray using a tuple of the minimum element, the index of the minimum and the length of the subarray. So you map each element `x` to `(x, 0, 1)` and then parallel fold using `f (x, i, n) (y, j, m) = if x < y then (x, i, n+m) else (y, n+j, n+m)`. For an identity element, you can probably get away with `(minBound, 0, 0)`. — hammar, Nov 26 '12 at 22:10
In order to use foldP with a tuple you will have to have array with those tuples. It can be done, but this is kind of wasteful for image processing (millions of elements). — EvgenijM86, Nov 27 '12 at 15:05
@EvgenjiM86: The idea was that the intermediate array of tuples would be fused away. — hammar, Nov 27 '12 at 15:15

Itai Zukerman · Answer 1 · 2012-11-29T21:25:45.013

I'm no expert on Repa, but this seems to work for 1-D arrays. It can probably be adapted for other dimensions.

import Data.Array.Repa

indexed arr = traverse arr id (\src idx@(Z :. i) -> (src idx, i))

minimize arr = foldP f h t
  where
    (Z :. n) = extent arr
    arr' = indexed arr
    h = arr' ! (Z :. 0)
    t = extract (Z :. 1) (Z :. (n-1)) arr'
    f min@(valMin, iMin) x@(val, i) | val < valMin = x
                                    | otherwise = min

score 0 · Answer 2 · answered Jun 29 '16 at 21:55

At the risk of reviving a zombie post, here is a arbitrary dimension solution I worked out (which, reviewing the comments, is exactly what @hammar suggested). Because of the weird nature of foldAllP, which I use as an identity element when merging tuples, you also need to provide an upper bound for the minimum of the array you are seeking to find.

import Data.Array.Repa
import Data.Array.Repa.Eval
import Data.Array.Repa.Repr.Unboxed

minimize :: (Shape sh, Source r a, Elt a, Unbox a, Monad m, Ord a) => Array r sh a -> a -> m (a,sh)
minimize arr upperBound = do
      -- Zip the array with its (Int representation of) indices
      let zippedWithIndex = Data.Array.Repa.traverse arr id (\f idx -> (f idx, toIndex (extent arr) idx))

      -- Define a function that compares two tuple of (a,Int) on the value of the first entry
      let minimumIndex t@(v,_) t'@(v',_) = if v < v' then t else t'

      -- Do a parallel fold of all elements in the array 
      (min,idx) <- foldAllP minimumIndex (upperBound, error "Empty array") zippedWithIndex

      -- convert the indice back from its Int representation
      return (min, fromIndex (extent arr) idx)

Naturally, if your array contains a type that has an instance of Bounded, you can make the more convenient function

 minimize' arr = minimize arr maxBound

Some tests you can do at the prompt:

 λ> let x = fromListUnboxed (ix2 2 2) [20, 30, 10, 40] :: Array U DIM2 Int
 λ> minimize' x
(10,(Z :. 1) :. 0)

Haskell repa - how to reduce array and return index?

2 Answers2