diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2015, Stephen Dekker + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Stephen Dekker nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/fractals.cabal b/fractals.cabal new file mode 100644 --- /dev/null +++ b/fractals.cabal @@ -0,0 +1,38 @@ +name: fractals +version: 0.1.0.0 +synopsis: A collection of useful fractal curve encoders +description: + A collection of efficient fractal curve encoders that are of general use + for creating spatial data structures. + . + Currently, the only encoder included is an implementation Butz's algorithm + for generating N-dimensional space-filling Hilbert curves. + . + An encoder for Morton (Z-order) curves is planned for a future release. + +license: BSD3 +license-file: LICENSE +copyright: Copyright (c) 2015, Stephen Dekker +author: Stephen Dekker +maintainer: Stephen Dekker +tested-with: GHC==7.10.2 +category: Math +build-type: Simple +cabal-version: >=1.10 + +library + exposed-modules: Data.SpaceFillingCurve.Hilbert.Integer + Data.SpaceFillingCurve.Hilbert.Integer.Internal + build-depends: base >= 4.8 && < 4.9 + hs-source-dirs: src + default-language: Haskell2010 + +test-suite property-tests + type: exitcode-stdio-1.0 + hs-source-dirs: test, src + main-is: TestSuite.hs + other-modules: PropertyTests + build-depends: base >= 4.8 && < 4.9, + QuickCheck >= 2.8.0 && < 2.9, + integer-gmp >= 1.0.0.0 && < 1.1 + default-language: Haskell2010 diff --git a/src/Data/SpaceFillingCurve/Hilbert/Integer.hs b/src/Data/SpaceFillingCurve/Hilbert/Integer.hs new file mode 100644 --- /dev/null +++ b/src/Data/SpaceFillingCurve/Hilbert/Integer.hs @@ -0,0 +1,51 @@ +-- | +-- Module : Data.SpaceFillingCurve.Hilbert.Integer +-- Copyright : (c) 2015 Stephen Dekker +-- License : BSD3 +-- +-- Maintainer : steve.dekk@gmail.com +-- Stability : experimental +-- Portability : portable +-- +-- An implementation of Butz's classic (and rather beautiful) algorithm for +-- computing the discrete Hilbert index of an N-dimensional point in +-- Cartesian space (A. R. Butz. "Alternative algorithm for Hilbert’s +-- space-filling curve. IEEE Transactions on Computers", pages 424–426, +-- April 1971). +-- +-- This particular implementation relies upon the 'Integer' numeric type in +-- order to handle unbounded input point coordinates. A version not built +-- around the 'Integer' type could offer improved performance, as the +-- algorithm essentially boils down to the repeated application of bitwise +-- operations. +-- +-- The specific algorithm used is the uncompact Hilbert indexing algorithm +-- described by Chris Hamilton (Hamilton, C. "Compact Hilbert Indices", +-- Dalhousie University, Faculty of Computer Science, Technical Report +-- CS-2006-07, July 2006). +-- +-- Hamilton's paper provides a thorough overview of the mathematics behind +-- the algorithm and also extends it to handle variable encoding widths for +-- the different Cartesian axes. The compact Hilbert indexing scheme +-- described in the technical report is not implemented in this module. +-- +-- The encoding function is written to accept a list of 'Bits' instances +-- for the input point and to produce a 'Num' instance for the output index +-- and is capable of handling unbounded 'Bits' instances such as 'Integer'. +-- +-- Similarly, the decoding function will take an unbounded 'Num' instance +-- and produce a point consisting of unbounded 'Bits' components with the +-- desired dimensionality. +-- +-- Lastly, the functions exported by this module will accept negative +-- inputs, but the behaviour of the functions for negative Hilbert indices +-- or point coordinates is undefined. These assumptions are made explicit +-- in the included QuickCheck property tests. + +module Data.SpaceFillingCurve.Hilbert.Integer ( + -- * Encoding and decoding the Hilbert curve + hilbert, -- :: (Bits a, Bits b, Num b) => Int -> [a] -> b + hilbertInverse, -- :: (Bits a, Bits b) => Int -> Int -> a -> [b] + ) where + +import Data.SpaceFillingCurve.Hilbert.Integer.Internal (hilbert, hilbertInverse) diff --git a/src/Data/SpaceFillingCurve/Hilbert/Integer/Internal.hs b/src/Data/SpaceFillingCurve/Hilbert/Integer/Internal.hs new file mode 100644 --- /dev/null +++ b/src/Data/SpaceFillingCurve/Hilbert/Integer/Internal.hs @@ -0,0 +1,177 @@ +-- | +-- Module : Data.SpaceFillingCurve.Hilbert.Integer.Internal +-- Copyright : (c) 2015 Stephen Dekker +-- License : BSD3 +-- +-- Maintainer : steve.dekk@gmail.com +-- Stability : experimental +-- Portability : portable +-- +-- This modules contains the implementation of Butz's Hilbert curve +-- encoding algorithm. Of these, the hilbert and hilbertInverse function +-- are exposed through the "Data.SpaceFillingCurve.Hilbert.Integer" module. + +module Data.SpaceFillingCurve.Hilbert.Integer.Internal ( + -- * Encoding and decoding the Hilbert curve + hilbert, -- :: (Bits a, Bits b, Num b) => Int -> [a] -> b + hilbertInverse, -- :: (Bits a, Bits b) => Int -> Int -> a -> [b] + + -- * Internal helper functions for the Hilbert transformations + bitAt, -- :: (Bits a, Bits b) => a -> Int -> b + trailingSetBits, -- :: (Bits a, Num b) => a -> b + mask, -- :: Num a => Int -> a + rotR, -- :: (Num a, Bits a) => Int -> a -> Int -> a + rotL, -- :: (Num a, Bits a) => Int -> a -> Int -> a + grayCode, -- :: Bits a => a -> a + grayCodeInverse, -- :: Bits a => a -> a + entryPoint, -- :: (Num a, Bits a) => a -> a + direction, -- :: (Num a, Bits a) => Int -> a -> Int + transform, -- :: (Num a, Bits a) => Int -> a -> Int -> a -> a + transformInverse, -- :: (Num a, Bits a) => Int -> a -> Int -> a -> a + pivot -- :: (Bits a, Bits b) => Int -> a -> [Int] -> [b] + ) where + +import Data.Bits (Bits, bit, clearBit, setBit, shiftL, shiftR, + testBit, xor, zeroBits, (.&.), (.|.)) + +------------------------------------------ +-- Encoding and decoding the Hilbert curve + +-- The variable names (symbols) used in Hamilton's paper are reproduced in +-- the bodies of these two functions. These names are overly terse, but are +-- useful when comparing the implementation side-by-side with Hamilton's +-- report. More descriptive names were chosen for the helper functions and +-- the sequences described in the original paper. + +-- | Given the number of bits required to represent the largest value in +-- the given input list (which represents a point in an N-dimensional +-- Cartesian space), returns the Hilbert index of the point. + +hilbert :: (Bits a, Bits b, Num b) => Int -> [a] -> b +hilbert precision ps = hilbertIndex + where (_, _, hilbertIndex) = foldr step start [0..precision-1] + n = length ps + start = (zeroBits, zeroBits, zeroBits) + step i (e, d, h) = seq e' $ seq d' $ seq h' (e', d', h') + where l = foldr (\x acc -> (acc `shiftL` 1) .|. (x `bitAt` i)) + zeroBits (reverse ps) + t = transform n e d l + w = grayCodeInverse t + h' = (h `shiftL` n) .|. w + e' = e `xor` rotL n (entryPoint w) (d+1) + d' = (d + direction n w + 1) `mod` n + +-- | Given the number of bits required to represent the largest value in +-- the output vector, the number of dimensions in the output space and the +-- Hilbert index of the output point, returns a list of values representing +-- the point in Cartesian space. + +hilbertInverse :: (Num a, Bits a, Bits b) => Int -> Int -> a -> [b] +hilbertInverse precision n hilbertIndex = points + where (_, _, points) = foldr step start [0..precision-1] + start = (0::Integer, zeroBits, replicate n zeroBits) + step i (e, d, ps) = seq e' $ seq d' $ seq ps' (e', d', ps') + where w = foldr (\x acc -> (acc `shiftL` 1) .|. + (hilbertIndex `bitAt` (i*n + x))) + zeroBits [0..n-1] + t = grayCode w + l = transformInverse n e d t + ps' = zipWith (.|.) ps (pivot i l (reverse [0..n-1])) + e' = e `xor` rotL n (entryPoint w) (d+1) + d' = (d + direction n w + 1) `mod` n + +------------------------------------------------------------ +-- Internal helper functions for the Hilbert transformations + +-- | Returns the value of the given bit in the source bit string. Note that +-- if the bit was set, the returned value will be of the output type with +-- only the first bit set. + +bitAt :: (Bits a, Bits b) => a -> Int -> b +bitAt x i = if x `testBit` i then bit 0 else zeroBits + +-- | Counts the number of trailing set bits in the given bit string. + +trailingSetBits :: (Bits a, Num b) => a -> b +trailingSetBits i = go i 0 + where go j acc = if not (testBit j 0) + then acc + else go (j `shiftR` 1) (acc+1) + +-- | Creates a bit mask extending the range of bits from [0, 'width' - 1]. + +mask :: Num a => Int -> a +mask width = 2^width - 2 + fromIntegral (signum width) + +-- | Performs a windowed right rotate by 'i' within a window from bit 0 to +-- bit 'width' on a number 'x'. + +rotR :: (Num a, Bits a) => Int -> a -> Int -> a +rotR width x i = trunc ((trunc x `shiftR` s) .|. (x `shiftL` (width - s))) + where s = i `mod` width + trunc = (.&.) (mask width) + +-- | Performs a windowed left rotate by 'i' within a window from bit 0 to +-- bit 'width' on a number 'x'. + +rotL :: (Num a, Bits a) => Int -> a -> Int -> a +rotL width x i = trunc ((x `shiftL` s) .|. (trunc x `shiftR` (width - s))) + where s = i `mod` width + trunc = (.&.) (mask width) + +-- | Returns the 'i'-th binary-reflected Gray code. + +grayCode :: Bits a => a -> a +grayCode i = i `xor` (i `shiftR` 1) + +-- | Returns the enumeration index of a given binary-reflected Gray code, +-- inverting the Gray code transform. + +grayCodeInverse :: Bits a => a -> a +grayCodeInverse g = go g g 1 + where go i acc j = if acc == zeroBits + then i + else go (i `xor` (g `shiftR` j)) + (acc `shiftR` 1) (j+1) + +-- | Returns the 'i'-th element in the sequence of entry points. + +entryPoint :: (Num a, Bits a) => a -> a +entryPoint i | signum i == -1 = error "Input must be positive" + | i == zeroBits = zeroBits + | otherwise = grayCode ((i-1) `clearBit` 0) + +-- | Given the dimensionality of the Hilbert curve and an index 'i', +-- returns the 'i'-th element in the sequence of directions. + +direction :: (Num a, Bits a) => Int -> a -> Int +direction n i | signum i == -1 = error "Input must be positive" + | i == zeroBits = zeroBits + | testBit i 0 = trailingSetBits i `mod` n + | otherwise = trailingSetBits(i-1) `mod` n + +-- | Given a dimensionality, an entry point, a direction and a Gray code +-- representing a canonical, unrotated sub-hypercube path we wish to +-- transform, returns the path rotated so that it is correctly oriented +-- within its quadrant. + +transform :: (Num a, Bits a) => Int -> a -> Int -> a -> a +transform n e d l = rotR n (l `xor` e) (d+1) + +-- | Given a dimensionality, an entry point, a direction and the Gray code +-- representing the rotated sub-hypercube path in a particular quadrant, +-- returns the path rotated and transformed back into its canonical form. + +transformInverse :: (Num a, Bits a) => Int -> a -> Int -> a -> a +transformInverse n e d l = e `xor` rotL n l (d+1) + +-- | Given a position, 'i', a bit-array 'l' and a list of positions to +-- test, returns a list containing either the value (2^'i') or 0 depending +-- on whether the bits in 'l' at the positions in the input list are set or +-- not. + +pivot :: (Bits a, Bits b) => Int -> a -> [Int] -> [b] +pivot i l = map (\j -> setBitIf (testBit l j) zeroBits i) + where setBitIf True = setBit + setBitIf False = const + diff --git a/test/PropertyTests.hs b/test/PropertyTests.hs new file mode 100644 --- /dev/null +++ b/test/PropertyTests.hs @@ -0,0 +1,178 @@ +{-# LANGUAGE MagicHash #-} + +-- | +-- Module : PropertyTests +-- Copyright : (c) 2015 Stephen Dekker +-- License : BSD3 +-- +-- Maintainer : steve.dekk@gmail.com +-- Stability : experimental +-- Portability : non-portable (MagicHash) + +module PropertyTests ( + runTests, -- :: IO () + prop_grayCodeInversion, -- :: Property + prop_hilbertInversion, -- :: Property + prop_hilbertLocality, -- :: Property + prop_rotRBitCountInvariance, -- :: Property + prop_rotLBitCountInvariance, -- :: Property + prop_rotRInteger, -- :: Property + prop_rotLInteger, -- :: Property + prop_rotLRInversion -- :: Property + ) where + +import Data.SpaceFillingCurve.Hilbert.Integer.Internal + +import GHC.Exts (Int(I#)) +import GHC.Integer.Logarithms (integerLog2#) + +import Data.Bits (clearBit, popCount, setBit, shiftL, shiftR, + testBit, (.&.), (.|.)) +import Test.QuickCheck (Gen, Property, arbitrary, choose, forAll, + listOf1, quickCheck, suchThat) + +------------------------ +-- Test value generators + +positiveInts :: Gen Int +positiveInts = (arbitrary :: Gen Int) `suchThat` (> 0) + +nonNegativeInts :: Gen Int +nonNegativeInts = (arbitrary :: Gen Int) `suchThat` (>= 0) + +nonNegativeIntegers :: Gen Integer +nonNegativeIntegers = (arbitrary :: Gen Integer) `suchThat` (>= 0) + +nonNegativeVectors :: Gen [Integer] +nonNegativeVectors = listOf1 ((arbitrary :: Gen Integer) `suchThat` (>= 0)) + + +-------------------------------------- +-- Some property tests for this module + +-- | Runs the QuickCheck property tests for the Hilbert encoder/decoder as +-- well as the internal helper functions. + +runTests :: IO () +runTests = do + quickCheck prop_rotLBitCountInvariance + quickCheck prop_rotRBitCountInvariance + quickCheck prop_grayCodeInversion + quickCheck prop_hilbertInversion + quickCheck prop_hilbertLocality + quickCheck prop_rotLRInversion + quickCheck prop_rotLInteger + quickCheck prop_rotRInteger + +-- | The Gray encoding is invertible. + +prop_grayCodeInversion :: Property +prop_grayCodeInversion = forAll nonNegativeInts $ \x -> + grayCodeInverse (grayCode x) == x + +-- | The Hilbert curve encoding is invertible. + +prop_hilbertInversion :: Property +prop_hilbertInversion = forAll nonNegativeVectors check + where check ps = hilbertInverse m' n (hilbert m' ps :: Integer) == ps + where m = unboundedBitSize (maximum ps) + n = length ps + m' = m + n - (m `mod` n) + +-- | Two points co-located on the Hilbert curve should be within a unit +-- step in Cartesian space. + +prop_hilbertLocality :: Property +prop_hilbertLocality = forAll nonNegativeVectors check + where check ps = norm (diff ps ps') `near` 1 + where h = hilbert m' ps :: Integer + ps' = hilbertInverse m' n (h+1) + n = length ps + m = unboundedBitSize (maximum ps) + m' = m + n - (m `mod` n) + +-- | The total population count of set bits must not change after a right +-- windowed rotate. + +prop_rotRBitCountInvariance :: Property +prop_rotRBitCountInvariance = forAll positiveInts $ \width -> + forAll nonNegativeIntegers $ \x -> + forAll nonNegativeInts $ \i -> check width x i + where check width x i = popCount (mask width .&. x) == popCount (rotR width x i) + +-- | The total population count of set bits must not change after a left +-- windowed rotate. + +prop_rotLBitCountInvariance :: Property +prop_rotLBitCountInvariance = forAll positiveInts $ \width -> + forAll nonNegativeIntegers $ \x -> + forAll nonNegativeInts $ \i -> check width x i + where check width x i = popCount (mask width .&. x) == popCount (rotL width x i) + +-- | Rotating a value right by one is equivalent to halving the value and +-- setting the last bit to 1 if the LSB in the initial value was set. + +prop_rotRInteger :: Property +prop_rotRInteger = forAll positiveInts $ \width -> + forAll nonNegativeIntegers $ \x -> + forAll positiveInts $ \i -> check width x i + where check width x i = rotR width x i == value + where y = rotR width x (i-1) + lsBit = testBit y 0 + value | lsBit = (y `shiftR` 1) `setBit` (width - 1) + | otherwise = y `shiftR` 1 + +-- | Rotating a value left by one is equivalent to doubling the value after +-- clearing the MSB and setting the first bit to 1 if the MSB in the +-- initial value was set. + +prop_rotLInteger :: Property +prop_rotLInteger = forAll nonNegativeIntegers $ \x -> + forAll (choose (1, unboundedBitSize x)) $ \width -> + forAll nonNegativeInts $ \i -> check width x i + where check width x i = rotL width x i == value + where y = rotL width x (i-1) + msBit = y `testBit` (width - 1) + value | msBit = ((y `clearBit` (width-1)) `shiftL` 1) .|. 1 + | otherwise = y `shiftL` 1 + +-- | Check that the windowed rotate functions are the inverse of each +-- other, taking into account the fact that both functions truncate the +-- input value to the window width. + +prop_rotLRInversion :: Property +prop_rotLRInversion = forAll positiveInts $ \width -> + forAll nonNegativeIntegers $ \x -> + forAll positiveInts $ \i -> check width x i + where check width x i = rotR width (rotL width x i) i == mask width .&. x + + +------------------------------------------- +-- Utility functions for the property tests + +-- | Calculates the number of bits required to represent an unbounded, +-- positive Integer. + +unboundedBitSize :: Integer -> Int +unboundedBitSize i | i == 0 = 1 + | otherwise = I# (integerLog2# i) + 1 + +-- | Calculates the norm of a given vector, returning the result as +-- a floating point number. + +norm :: (Integral a, Floating b) => [a] -> b +norm a = sqrt (sum (map (fromIntegral . (^(2 :: Int))) a)) + +-- | Calculates the difference between two vectors. + +diff :: Num a => [a] -> [a] -> [a] +diff = zipWith (-) + +-- | Determines whether or not two 'Double' precision numbers are near +-- enough to be considered equal (that is, the displacement between them is +-- less than the machine epsilon). + +near :: Double -> Double -> Bool +near a b = abs (a - b) < epsilon + where epsilon = 2**(-53) :: Double + diff --git a/test/TestSuite.hs b/test/TestSuite.hs new file mode 100644 --- /dev/null +++ b/test/TestSuite.hs @@ -0,0 +1,6 @@ +module Main where + +import PropertyTests + +main :: IO () +main = runTests