fractals (empty) → 0.1.0.0

raw patch · 7 files changed

+482/−0 lines, 7 filesdep +QuickCheckdep +basedep +integer-gmpsetup-changed

Dependencies added: QuickCheck, base, integer-gmp

Files

• LICENSE +30/−0
• Setup.hs +2/−0
• fractals.cabal +38/−0
• src/Data/SpaceFillingCurve/Hilbert/Integer.hs +51/−0
• src/Data/SpaceFillingCurve/Hilbert/Integer/Internal.hs +177/−0
• test/PropertyTests.hs +178/−0
• test/TestSuite.hs +6/−0

@@ -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.

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

@@ -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 <steve.dekk@gmail.com>+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

@@ -0,0 +1,51 @@+-- |+-- Module      : Data.SpaceFillingCurve.Hilbert.Integer+-- Copyright   : (c) 2015 Stephen Dekker <steve.dekk@gmail.com>+-- 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)

@@ -0,0 +1,177 @@+-- |+-- Module      : Data.SpaceFillingCurve.Hilbert.Integer.Internal+-- Copyright   : (c) 2015 Stephen Dekker <steve.dekk@gmail.com>+-- 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+

@@ -0,0 +1,178 @@+{-# LANGUAGE MagicHash #-}++-- |+-- Module      : PropertyTests+-- Copyright   : (c) 2015 Stephen Dekker <steve.dekk@gmail.com>+-- 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+

@@ -0,0 +1,6 @@+module Main where++import PropertyTests++main :: IO ()+main = runTests