5

I built a function for finding the determinant of a matrix using the ST-Monad and unboxed STArrays (STUArray). The type for a matrix is the following:

newtype Matrix e = Matrix (Array Int (UArray Int e))

that is, an immutable array containing immutable unboxed arrays containing the elements. This will require me to add the Predicate IArray UArray e to functions dealing with Matrix, which in turn requires FlexibleContexts. Okay, done.

The function used to calculate the determinant has the following signature:

detST :: (IArray UArray e, MArray (STUArray s) e (ST s),
          Num e, Eq e, Division e)
   => Array Int (UArray Int e) -> ST s e

I am required to also add the Predicate MArray (STUArray s) e (ST s) since internally the arrays are converted into mutable arrays (the outer boxed, the inner unboxed).

This function can be used like so:

main = do
    let m@(Matrix x) = matrix [ [1,-2,3,234]
                              , [5,2,3,-3]
                              , [7,18,3,40]
                              , [2,9,71,0] ]
        d = runST (detST x) :: Int -- needed for type check, ambiguous otherwise

print d

Works, fine. But look how ugly it is! Of course I do not want to give away the internals of of Matrix (at least not further than the predicates attached to my functions already make me to). I would like to define the following function:

det :: Matrix e -> e

And I can't.

I tried without a proper signature:

det (Matrix arr) = runST (detST arr)

fails. Fair enough, I'll put my brain to work: detST requires IArray UArray e, MArray (STUArray s) e (ST s), Num e, Eq e, Division e, so does det, doesn't it?

det :: (IArray UArray e, MArray (STUArray s) e (ST s),
          Num e, Eq e, Division e) => Matrix e -> e

fails. But I don't see how. The message that GHC (7.4.2) gives me is:

Could not deduce (MArray (STUArray s) t (ST s))
  arising from a use of `detST'

but that exact term is in the predicates!

GHC suggests:

add (MArray (STUArray s) t (ST s)) to the context of
  a type expected by the context: ST s t
  or the inferred type of
     det :: (Eq t, Num t, IArray UArray t, Division t) => Matrix t -> t
or add an instance declaration for (MArray (STUArray s) t (ST s))

Okay. As to my understanding I have done that first thing. Also there does exist an instance for that (MArray ...) (otherwise, how could I have used it successfully in main?!).

I am not sure about what is wrong here. I believe it has something to do with the "hidden" ST state in s, and that the s of detST is some other s than the s in det would be, but I don't know how to write a type for this.

What is the proper definition of det - and why?!

The program without det compiles fine by using only FlexibleContexts, no warnings with -Wall. The complete source code can be found as a gist here.

scravy
  • 11,904
  • 14
  • 72
  • 127
  • 1
    To begin with - why not just use an UArray indexed with (Int, Int) instead of an array of arrays? – Mikhail Glushenkov Apr 05 '13 at 07:54
  • So that I can easily switch the rows. Anyways, there are of course other possible versions of this program. But I want to know why I can't define that particular function and what is wrong with it. – scravy Apr 05 '13 at 07:56
  • Can you show us the definition of `detST` so that we can test your code? I don't quite understand why it needs to be in `ST`. – Mikhail Glushenkov Apr 05 '13 at 08:01
  • The complete code is linked as a gist, see the last four words of the question. It needs to be in ST since it works with an STUArray internally. i.e. detST is the function which should be the internal function, and det the one I would publicly export. – scravy Apr 05 '13 at 08:08

1 Answers1

4

I managed to get this to work using the trick described by Keegan McAllister in this article:

{-# LANGUAGE FlexibleContexts, ScopedTypeVariables, RankNTypes, GADTs #-}

data Evidence s e where
  Evidence :: (MArray (STUArray s) e (ST s)) => Evidence s e

data ElemType e = ElemType (forall s. Evidence s e)

det :: forall e . (IArray UArray e, Num e, Eq e, Division e)
       => ElemType e -> Matrix e -> e
det (ElemType e) mat = runST (f e mat)
  where
    f :: Evidence s e -> Matrix e -> ST s e
    f Evidence (Matrix arr) = detST arr

Usage:

main :: IO ()
main = do
    let m = matrix [ [1,-2,3,234]
                   , [5,2,3,-3]
                   , [7,18,3,40]
                   , [2,9,71,0] ]
    print $ det (ElemType Evidence) (m :: Matrix Int)

The problem stems from the lack of higher-rank constraints - runST has type (forall s. ST s a) -> a, so you need a constraint like forall s . MArray (STUArray s) e (ST s), which are not supported by GHC. This trick allows you to convince the type checker that the constraint actually holds. A more in-depth discussion of this problem is available in the article I linked above.

Mikhail Glushenkov
  • 14,928
  • 3
  • 52
  • 65