Features Download
From: <oleg <at> pobox.com>
Subject: Eliminating Multiple-Array Bound Checking through Non-dependent types
Newsgroups: gmane.comp.lang.haskell.general
Date: Saturday 11th February 2006 06:05:04 UTC (over 12 years ago)
An earlier message showed an example of writing code with non-trivial
static guarantees in the present-day Haskell (i.e., Haskell98 + rank-2


Although this approach obviously falls short of the genuine dependent
typing, we can use it in the existing, well-maintained functional
language systems. Thus we may begin to acquire taste for dependent
typing -- which, McBride at al. point out, matter.

The present message shows a more advanced example -- eliminating array
bound checking when processing several arrays at a time. The number of
arrays to process is not statically known. Furthermore, the arrays may
have different sizes and bounds -- potentially, empty and
non-overlapping bounds. The example involves native Haskell arrays,
index computations, and general recursion. All array indexing
operations are statically guaranteed to be safe -- and so we can
safely use an efficient unsafeAt provided by GHC seemingly for that
purpose. The code is efficient; the static assurances in the main loop
over arrays cost us no run-time overhead. The example uses only
Haskell98 + higher-ranked types. No new type classes are
introduced. The safety is based on: Haskell type system, quantified
type variables, and a compact general-purpose trusted kernel.
I thank Daniel Yokomizo for the challenge.

Our example is folding over multiple, variously-sized arrays.  This
is like a fold over an array -- generalized to an arbitrary number of
arrays, whose index ranges do not have to be the same (and do not have
to overlap). Typing this example in a genuinely dependent type system
is probably going to be quite challenging.

Our goal is to implement a function

 marray_fold :: (Ix i, Integral i) =>
	         (seed -> [e] -> seed) -> seed -> [Array i e] -> seed

Its third argument is a list of arrays; the arrays have all the same
element and index types; the actual sizes (that is, the lower and
upper bounds) may differ. Some arrays in the list may even be empty
(with the lower bound higher than the upper one). The function
marray_fold, like left fold, applies its left argument to the
values extracted from the corresponding arrays. Because arrays may
have different sizes and bounds, marray_fold operates over the range
of indices that is the intersection of the ranges of the argument

For example:

  dot a1 a2 = marray_fold (\seed l -> product l + seed) 0 [a1,a2]

computes the dot products of two arrays. 

Granted our example is a bit contrived: in practice one would probably
use ``fixed arity'' multiple array fold (similar to zipWith, zipWith3,
etc -- which are more efficient and simpler to implement).  Also
simple is the processing of arrays of the same bounds (and invoking an
error continuations if the bounds turn out different).  OTH, our
present example is more challenging and thus more interesting.

Here's the complete implementation (assuming the short trusted
branding library explained below).

> {-# OPTIONS -fglasgow-exts #-}
> module MDep where
> import Data.Array
> marray_fold :: (Ix i, Integral i) =>
>	         (seed -> [e] -> seed) -> seed -> [Array i e] -> seed
> marray_fold _ seed [] = seed		-- No arrays, no fold
> marray_fold f seed arrs = brands arrs (marray_fold' f seed) seed
> marray_fold' f seed (barrs, bl, bh) = loop (bl2i bl) seed
>   where loop bi seed = bsucc bh bi (\bi' -> loop bi' seed') seed'
>	   where seed' = f seed (map (!. bi) barrs)
> brands :: (Ix i, Integral i) => [Array i e]
>                   -> (forall s. ([BArray s i e], BLow s i, BHi s i) -> w)
>	            -> w -> w
> brands [arr] consumer onempty =
>     brand arr (\ (barr,bl,bh) -> consumer ([barr],bl,bh)) onempty
> brands (a:arrs) consumer onempty =
>     brands arrs (\bbars -> brand_merge bbars a consumer onempty) onempty

The function |brands| computes the intersection index range for the
given arrays, and brands each array and the intersection range. The
range is given as a pair (BLow s i, BHi s i). The quantified (and thus
unforgeable) type variable |s| represents the brand. The function
|brands| has a higher-ranked type and must have a
signature. Incidentally, functions like that are relatively rare;
processors of the list of arrays of the same brand do not require any
signature (cf. marray_fold').

The folding loop has as few index comparison operations as
algorithmically needed, i.e., only one per iteration. Operator (!.) is
a statically safe array indexing operator. The type system and the
trust properties of the kernel below guarantee that in the expression
"arr !. m" the index `m' is positively in range of the array `arr'

Here are a few test cases, of the dot product over several arrays

> dots ars = marray_fold (\seed l -> product l + seed) 0 ars
> test1 = dots [listArray (1,5) [1..]]
> test2 = dots [listArray (1,5) [1..],listArray (1,5) [1,1..]]
> test3 = dots [listArray (1,5) [1..],listArray (5,1) [1,1..]]
> test4 = dots [listArray (1,5) [1..],listArray (2,7) [1,1..]]
> test5 = dots [listArray (1,5) [1..],listArray (7,7) [1,1..]]
> test6 = dots [listArray (1,5) [1..],listArray (5,5) [1,1..],
>               listArray (3,7) [3..]]

The case test3 demonstrates the inclusion of an empty array; test5
shows non-overlapping arrays.

The code relies on the compact general-purpose trusted kernel
explained below. That code should be put into a separate module.

First we introduce tags for Branded arrays, Branded indices,
and Branded low and high bounds.

>     -- the data constructors of those four must *not* be exported!
> newtype BArray s i a = BArray (Array i a)
> newtype BIndex s i   = BIndex i
> newtype BLow s i     = BLow i
> newtype BHi  s i     = BHi i
> unbi (BIndex i) = i

These are `newtype's and so impose no run-time overhead. Of interest
is a phantom type variable 's', which marks a _brand_ of an array and
of an array index. An index or a bound are branded if they are
positively within the range of the array of their brand.  The type
variable 's' is similar to that in the ST monad. The latter relies on
's' to enforce serialization. We, OTH, do not impose any linearity
constraints on 's' -- it may be freely duplicated and discarded (see
`unbi'). It is created however under controlled conditions.  The
safety depends on the trusted way of creating branded types: the data
constructors BIndex, BLow, BHi, and BArray must be used in the trusted
kernel only, and should not be available anywhere else.  The
uniqueness of 's' afforded by the explicit universal quantification
prevents mixing up of different brands.

We must re-iterate that safety depends on assurances of the code that
constructs branded values. Because of the high assurance, we must
formulate the safety properties as propositions, and prove
them. Fortunately, the code below is compact and straightforward, as
well as general purpose.

We start with an introduction rule for BArray:

> brand:: (Ix i, Integral i) 
>  => Array i e 
>      -> (forall s. (BArray s i e, BLow s i, BHi s i) -> w)
>      -> w -> w
> brand (a::Array i e) k kempty = 
>     let (l,h) = bounds a
>     in if l <= h then k ((BArray a)::BArray () i e, BLow l, BHi h)
>     else kempty

The function takes an array and two continuations. If the array is
empty, kempty is returned. Otherwise, the continuation k is invoked,
which receives the branded array along with two branded indices, for
the lower and the upper bounds of the array.

Proposition: For the argument (a,l,h) of the continuation k, 
the branded integral indices l and h are within the bounds of the non-empty

branded array a.
Proof: from the semantics of the functions `bounds' and `inRange',
properties of integral numbers and ranges over integral numbers
(specifically, contiguity), and the code itself.

The function brand has a higher-rank type. It is the existential
quantification of 's' as well as the absence of BArray constructor
elsewhere guarantee that the same brand entails the same bounds.

For this example, we need another introduction rule for BArray:

> brand_merge :: (Ix i, Integral i) => 
>	 ([BArray s i e], BLow s i, BHi s i)
>	 -> Array i e
>        -> (forall s'. ([BArray s' i e], BLow s' i, BHi s' i) -> w)
>        -> w -> w
> brand_merge (bas,(BLow bl),(BHi bh)) (a :: Array i e) k kempty =
>     let (l,h) = bounds a
>	  l' = max l bl
>         h' = min h bh
>     in if l' <= h' then 
>	       k (((BArray a)::BArray () i e) :
>	           (map (\ (BArray a) -> BArray a) bas),
>	          BLow l', BHi h')
>     else kempty

The function |brand_merge| takes a list of branded arrays and their
greatest common bound (specified as a pair of branded indices). The
function also takes an unbranded array. We prepend the new array to
the existing ones, and adjust the greatest common bound. The result
will have a new brand, which signifies that the new |bl| and |bh|
indices are within the range of all of the arrays.

The mapped function  (\ (BArray a) -> BArray a) is operationally identity.
Hopefully the compiler recognizes this fact.

Proposition: Given (bas, l, h) where 
      [l,h] is an integral range,
      forall arr in bas. [l,h] in bounds(arr)
   if the continuation k is ever invoked, its argument (bas',l',h')
   has the property: 
      bas' is a non-empty list,
      [l',h'] is a non-empty linear (integral) range,
      forall arr in bas'. [l',h'] in bounds(arr)

Proof: from the semantics of the functions `bounds' and `inRange',
properties of integral numbers and ranges over integral numbers
(specifically, contiguity), and the code itself.

> bl2i:: BLow s i -> BIndex s i
> bl2i (BLow i) = BIndex i

Proposition: if bl is within the bounds of an array branded with s,
then the result is BIndex within the bounds of that array.
Proof: trivial.

> bsucc:: (Integral i, Ix i) 
> 	=> BHi s i -> BIndex s i -> (BIndex s i -> r) -> r -> r
> bsucc (BHi upb) (BIndex i) on_within on_out
> 	= let i'    = succ i
> 	  in if i' <= upb then (on_within (BIndex i')) else on_out

The function `bsucc' takes two branded indices that correspond to the
same bounds (see the variable 's'). The function also takes two
continuations, on_within and on_out. The first index is considered to
be an upper limit. The function increments the second index. If the
result does not exceed the upper limit, we invoke on_within and pass
it the result. Otherwise, we invoke `on_out'.

Safety Proposition: l <= upb <= h, l <= i <= h, (i+1) <= upb |-
		    l <= (i+1) <= h
Proof: from i < (i+1) and properties of inequalities.
The safety proposition justifies our use of the data constructor

> bpred:: (Ix i, Integral i)
> 	=> BLow s i -> BIndex s i -> (BIndex s i -> r) -> r -> r
> bpred (BLow lwb) (BIndex i) on_within on_out
> 	= let i'    = pred i
> 	  in if i' >= lwb then (on_within (BIndex i')) else on_out

The dual of `bsucc'.
Safety Proposition: l <= lwb <= h, l <= i <= h, (i-1) >= lwb |-
		    l <= (i-1) <= h

Because a branded index is assuredly within the bounds of the array of
the same brand, we can write

> infixl 5 !.
> (!.):: (Ix i) => BArray s i e -> (BIndex s i) -> e
> (BArray a) !. (BIndex i) = a ! i

actually, we may _safely_ replace a ! i with `unsafeAt a i'

A reader may have noticed a seemingly artificial |Integral i|
constraint on indices. Why can't the mere |Ix i| constraint suffice?
That would let us potentially deal with multi-dimensional
arrays. Alas, the sole |Ix i| constraint does not suffice because the
Ix laws guaranteed in Haskell are lame. For example:

  inRange (l,h) l ?
   Does not always hold: inRange (2,1) 2 ---> False

  if rangeSize (l,h) >0, inRange (l,h) l ?
   does not seem to follow from the provided laws

  l <= h, inRange(l,h) l ?
   False: let l = (1,2); h = (2,1) in (l<=h, inRange (l,h) l)
          gives (True,False)

let (bl,bh) = ((1,2),(1,2)) -- non-empty range
    (l,h)   = ((1,1),(2,1)) -- non-empty range
    l' = max l bl -- (1,2)
    h' = min h bh -- (1,2)
 rangeSize (l',h') is positive
 but inRange (l,h) l' is False

Although Ord is a super-class of Ix, the comparison operation of Ord
is not required to have any relationship with ranges. And in fact, it
doesn't: the comparison operation on tuples is lexicographic, which
isn't the best choice for tuples representing multi-dimensional

Perhaps this message might provide the incentive for more precise
specification of the laws for various basic classes in future Haskell
CD: 4ms