How to translate this list-based code into using mutable arrays?

Question

EDIT3: I'm writing a code to process very long input list of Ints with only few hundred non-duplicates. I use two auxiliary lists to maintain cumulative partial sums to calculate some accumulator value, the how's and why's are non-important. I want to ditch all lists here and turn it into nice destructive loop, and I don't know how. I don't need the whole code, just a skeleton code would be great, were read/write is done to two auxiliary arrays and some end result is returned. What I have right now would run 0.5 hour for the input. I've coded this now in C++, and it runs in 90 seconds for the same input.

I can't understand how to do this, at all. This is the list-based code that I have right now:^{(but the Map-based code below is clearer)}

ins :: (Num b, Ord a) => a -> b -> [(a, b)] -> ([(a, b)], b)
ins n x [] = ( [(n,x)], 0) 
ins n x l@((v, s):t) = 
  case compare n v of
    LT -> ( (n,s+x) : l , s )
    EQ -> ( (n,s+x) : t , if null t then 0 else snd (head t))
    GT -> let (u,z) = ins n x t
          in  ((v,s+x):u,z)

This is used in a loop, to process a list of numbers of known length, ^{(changed it to foldl now)}

scanl g (0,([],[])) ns  -- ns :: [Int]
g ::
  (Num t, Ord t, Ord a) =>
  (t, ([(a, t)], [(a, t)])) -> a -> (t, ([(a, t)], [(a, t)])) 
g (c,( a, b)) n = 
    let
      (a2,x) = ins n 1 a
      (b2,y) = if x>0 then ins n x b else (b,0)
      c2     = c + y
    in
      (c2,( a2, b2))

This works, but I need to speed it up. In C, I would keep the lists (a,b) as arrays; use binary search to find the element with the key just above or equal to n (instead of the sequential search used here); and use in-place update to change all the preceding entries.

I'm only really interested in final value. How is this done in Haskell, with mutable arrays?

I tried something, but I really don't know what I'm doing here, and am getting strange and very long error messages (like "can not deduce ... from context ..."):

goarr top = runSTArray $ do
  let sz = 10000
  a <- newArray (1,sz) (0,0) :: ST s (STArray s Int (Integer,Integer))
  b <- newArray (1,sz) (0,0) :: ST s (STArray s Int (Integer,Integer))
  let p1 = somefunc 2 -- somefunc :: Integer -> [(Integer, Int)]
  go1 p1 2 0 top a b

go1 p1 i c top a b = 
    if i >= top
     then 
      do
       return c
     else
       go2 p1 i c top a b

go2 p1 i c top a b =
  do
   let p2 = somefunc (i+1)  -- p2 :: [(Integer, Int)]
   let n  = combine p1 p2   -- n :: Int
   -- update arrays and calc new c 
   -- like the "g" function is doing:
   --    (a2,x) = ins n 1 a
   --    (b2,y) = if x>0 then ins n x b else (b,0)
   --    c2     = c + y
   go1 p2 (i+1) c2 top a b  -- a2 b2??

This doesn't work at all. I don't even know how to encode loops in do notation. Please help.

UPD: the Map based code that runs 3 times slower:

ins3 :: (Ord k, Num a) => k -> a -> Map.Map k a -> (Map.Map k a, a)
ins3 n x a | Map.null a = (Map.insert n x a , 0)
ins3 n x a = let (p,q,r) = Map.splitLookup n a in
  case q of 
    Nothing -> (Map.union (Map.map (+x) p) 
                 (Map.insert n (x+leftmost r) r) , leftmost r)
    Just s -> (Map.union (Map.map (+x) p) 
                 (Map.insert n (x+s) r) , leftmost r)

leftmost r | Map.null r = 0
           | otherwise = snd . head $ Map.toList r

UPD2: The error message is " Could not deduce (Num (STArray s1 i e)) from the context () arising from the literal `0' at filename.hs:417:11"

that's where it says return c in go1 function. Perhaps c is expected to be an array, but I want to return the accumulator value that is built while using the two auxiliary arrays.

EDIT3: I've replaced scanl and (!!) with foldl and take as per Chris's advice, and now it runs in constant space with sane empirical complexity and is actually projected to finish in under 0.5 hour - a.o.t. ... 3 days ! I knew about it of course but was so sure GHC optimizes the stuff away for me, surely it wouldn't make that much of a difference, I thought! And so felt only mutable arrays could help... Bummer.

Still, C++ does same in 90 sec, and I would very much appreciate help in learning how to code this with mutable arrays, in Haskell.

the second half is mostly gibberish as I really dont know what I'm doing. The first half is a working code, it just calculates something in a loop, while maintaining two auxiliary lists - which I want to turn into arrays, for speed. — darveter
– darveter, Commented Apr 2, 2012 at 21:46
The first half is still too difficult follow. Could we get some type signatures, maybe? Or some comments? — Louis Wasserman
– Louis Wasserman, Commented Apr 2, 2012 at 21:48
I endorse Louis' request. On a different note, inserting into a sorted array may not be the most efficient way. Perhaps you can use Maps for it? — Daniel Fischer
– Daniel Fischer, Commented Apr 2, 2012 at 21:51
with maps it runs 3 times slower, because it not just inserts, but also must update all entries before the insertion point, and so must union the updated first half with the second half. — darveter
– darveter, Commented Apr 2, 2012 at 22:08

Chris Kuklewicz · Accepted Answer · 2012-04-03 20:00:29Z

3

Are the input values ever EQ? If they are not EQ then the way scanl g (0,([],[])) ns is used means that the first [(,)] array, call it a always has map snd a == reverse [1..length a] at each stage of g. For example, in a length 10 list the value of snd (a !! 4) is going to be 10-4. Keeping these reversed index values by mutating the second value of each preceding entry in a is quite wasteful. If you need speed then this is one place to make a better algorithm.

None of this applies to the second [(,)] whose purpose is still mysterious to me. It records all insertions that were not done at the end of a, so perhaps it allows one to reconstruct the initial sequence of values.

You said "I'm only really interested in final value." Do you mean you only care about the last value in list output by the scanl .. line? If so then you need a foldl instead of scanl.

Edit: I am adding a non-mutable solution using a custom Finger Tree. It passes my ad hoc testing (at bottom of code):

{-# LANGUAGE MultiParamTypeClasses #-}
import Data.Monoid
import Data.FingerTree

data Entry a v = E !a !v deriving Show

data ME a v = NoF | F !(Entry a v) deriving Show

instance Num v => Monoid (ME a v) where
  mempty = NoF
  NoF `mappend` k = k
  k `mappend` NoF = k
  (F (E _a1 v1)) `mappend` (F (E a2 v2)) = F (E a2 (v1 + v2))

instance Num v => Measured (ME a v) (Entry a v) where
  measure = F

type M a v = FingerTree (ME a v) (Entry a v)

getV NoF = 0
getV (F (E _a v)) = v

expand :: Num v => M a v -> [(a, v)]
expand m = case viewl m of
             EmptyL -> []
             (E a _v) :< m' -> (a, getV (measure m)) : expand m'

ins :: (Ord a, Num v) => a -> v -> M a v -> (M a v, v)
ins n x m =
  let comp (F (E a _)) = n <= a
      comp NoF = False
      (lo, hi) = split comp m
  in case viewl hi of
       EmptyL -> (lo |> E n x, 0)
       (E v s) :< higher | n < v ->
         (lo >< (E n x <| hi), getV (measure hi))
                         | otherwise ->
         (lo >< (E n (s+x) <| higher), getV (measure higher))

g :: (Num t, Ord t, Ord a) =>
     (t, (M a t, M a t)) -> a -> (t, (M a t, M a t))
g (c, (a, b)) n =
  let (a2, x) = ins n 1 a
      (b2, y) = if x>0 then ins n x b else (b, 0)
  in (c+y, (a2, b2))

go :: (Ord a, Num v, Ord v) => [a] -> (v, ([(a, v)], [(a, v)]))
go ns = let (t, (a, b)) = foldl g (0, (mempty, mempty)) ns
        in (t, (expand a, expand b))

up = [1..6]
down = [5,4..1]
see'tests = map go [ up, down, up ++ down, down ++ up ]

main = putStrLn . unlines . map show $ see'test

edited Apr 3, 2012 at 20:00

answered Apr 3, 2012 at 13:16

Chris Kuklewicz

8,16325 silver badges33 bronze badges

Sign up to request clarification or add additional context in comments.

8 Comments

darveter Over a year ago

yes, of very long input list there only few hundred non-duplicates. The two lists just maintain cumulative partial sums to calculate some accumulator value, the hows and whys are non-important. In C each ins operation would be essentially O(1). Yes on a foldl, I used scanl for debugging purposes. Anyway I want to ditch all lists here and turn it into nice destructive loop, and I don't know how. I don't need the whole code, just a skeleton code would be great, were read/write is done to two auxiliary arrays and some end result is returned.

darveter Over a year ago

Thanks for your suggestion of foldl, it turned projected run time of 3 days to 0.5 hour!! (It would've been enough, yesterday, to run it for 0.5 hour. The irony!). In C++ it takes 90 sec though, and I'd still like to learn how to code this with mutable arrays in Haskell. It not nice to feel ignorant, and I can't make sense of anything on haskellwiki or learnyouahaskell etc. about this. I've updated the Q.

Chris Kuklewicz Over a year ago

@user1308992 : The fingertree code I just posted should be O(log n) for ins instead of updating O(n) entries each time.

darveter Over a year ago

A. It seems what I have is almost Fenwick tree; I have O(m) update, WP says it should be O(log m), m=number of keys; B. I don't know FingerTree module, it'll take me forever to get all the |> and ><s; I understand the concept, have read apfelmus' blog entry; C. I could easily have O(log n) upd at cost of storing individual frequencies instead of cumulative, but I need the sum to calc the increment on each step. Is there any summing going on deep in the M's guts on each increment in your code? And if not and it is all delayed and stored somewhere, will it create a massive space leak...?

darveter Over a year ago

btw I only need fst $ go ns; the two dictionaries just build/pass the info forward as foldl g rolls over, so there is no need for expand (except in data-debugging).

|

Chris Kuklewicz · Accepted Answer · 2012-04-04 09:39:28Z

2

Slightly unorthodox, I am adding a second answer using a mutable technique. Since user1308992 mentioned Fenwick trees, I have used them to implement the algorithm. Two STUArray are allocated and mutated during the run. The basic Fenwick tree keeps totals for all smaller indices and the algorithm here needs totals for all larger indices. This change is handled by the (sz-x) subtraction.

import Control.Monad.ST(runST,ST)
import Data.Array.ST(STUArray,newArray)
import Data.Array.Base(unsafeRead, unsafeWrite)
import Data.Bits((.&.))
import Debug.Trace(trace)
import Data.List(group,sort)

{-# INLINE lsb #-}
lsb :: Int -> Int
lsb i = (negate i) .&. i

go :: [Int] -> Int
go xs = compute (maximum xs) xs

-- Require "top == maximum xs" and "all (>=0) xs"
compute :: Int -> [Int] -> Int
compute top xs = runST mutating where
  -- Have (sz - (top+1)) > 0 to keep algorithm simple
  sz = top + 2

  -- Reversed Fenwick tree (no bounds checking)
  insert :: STUArray s Int Int -> Int -> Int -> ST s ()
  insert arr x v = loop (sz-x) where
    loop i | i > sz = return ()
           | i <= 0 = error "wtf"
           | otherwise = do
      oldVal <- unsafeRead arr i
      unsafeWrite arr i (oldVal + v)
      loop (i + lsb i)

  getSum :: STUArray s Int Int -> Int -> ST s Int
  getSum arr x = loop (sz - x) 0 where
     loop i acc | i <= 0 = return acc
                | otherwise = do
       val <- unsafeRead arr i
       loop (i - lsb i) $! acc + val

  ins n x arr = do
    insert arr n x
    getSum arr (succ n)

  mutating :: ST s Int
  mutating = do
    -- Start index from 0 to make unsafeRead, unsafeWrite easy
    a <- newArray (0,sz) 0 :: ST s (STUArray s Int Int)
    b <- newArray (0,sz) 0 :: ST s (STUArray s Int Int)
    let loop [] c = return c
        loop (n:ns) c = do
          x <- ins n 1 a
          y <- if x > 0
               then 
                 ins n x b
               else
                 return 0
          loop ns $! c + y
    -- Without debugging use the next line
    -- loop xs 0
    -- With debugging use the next five lines
    c <- loop xs 0
    a' <- see a
    b' <- see b
    trace (show (c,(a',b'))) $ do 
    return c

  -- see is only used in debugging
  see arr = do
    let zs = map head . group . sort $ xs
    vs <- sequence [ getSum arr z | z <- zs ]
    let ans = filter (\(a,v) -> v>0) (zip zs vs)
    return ans

up = [1..6]
down = [5,4..1]
see'tests = map go [ up, down, up ++ down, down ++ up ]

main = putStrLn . unlines . map show $ see'tests

answered Apr 4, 2012 at 9:39

Chris Kuklewicz

8,16325 silver badges33 bronze badges

3 Comments

darveter Over a year ago

Thank you so very much for your incredibly generous help! Not only I have a code to study mutable arrays now, but also a clear code for a Fenwick tree! One thing if I may: the input list is very long; your code is not "online" in two places. I can guesstimate the top value instead of calling maximum; but in see you use xs to find out all unique keys in the input. This info is available in the first tree, as it counts each incoming key. So all elts of the first tree with non-zero frequencies are exactly the keys that we need to see in the second tree. Thanks again!

Chris Kuklewicz Over a year ago

Well. Entries with non-zero keys in a Fenwick tree may not have been inserted directly, the insert may add to several entries. By computing all the running totals you can detect which ones are bigger than the previous one and this does indicate an inserted key.

darveter Over a year ago

That is what I meant. Individual frequency, not cumulative frequency. Fenwick tree ought to provide for querying of both, both in O(log n) time. I guess it's trivial with getFrq arr k = do {a<-getSum arr k; b<-getSum arr(k+1); return (b-a)}, right? Or here in see with using let ans = filter (\(a,v) -> v>0) (zipWith(\(a,v)(b,u)->(a,u-v)) vs (tail vs)) (or is it ->(b,u-v)?). Thanks again!

Collectives™ on Stack Overflow

How to translate this list-based code into using mutable arrays?

2 Answers 2

8 Comments

3 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

8 Comments

3 Comments

Your Answer

Sign up or log in

Post as a guest

Related