persistent-spatial (empty) → 0.1.0.0

raw patch · 8 files changed

+647/−0 lines, 8 filesdep +QuickCheckdep +aesondep +basesetup-changed

Dependencies added: QuickCheck, aeson, base, hspec, http-api-data, integer-logarithms, lens, persistent, persistent-spatial, text

Files

• CHANGELOG.md +5/−0
• LICENSE +21/−0
• README.md +4/−0
• Setup.hs +2/−0
• persistent-spatial.cabal +53/−0
• src/Data/LatLong.hs +179/−0
• src/Data/Morton.hs +224/−0
• test/Spec.hs +159/−0

@@ -0,0 +1,5 @@+# Changlog for persistent-spatial++## 0.1.0.0++Initial release, start of versioning.

@@ -0,0 +1,21 @@+MIT License++Copyright © 2018-2019 Satsuma Labs++Permission is hereby granted, free of charge, to any person obtaining a copy+of this software and associated documentation files (the "Software"), to deal+in the Software without restriction, including without limitation the rights+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell+copies of the Software, and to permit persons to whom the Software is+furnished to do so, subject to the following conditions:++The above copyright notice and this permission notice shall be included in all+copies or substantial portions of the Software.++THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE+SOFTWARE.

@@ -0,0 +1,4 @@+# persistent-spatial++This package implements a type for storing and indexing geographic coordinates which can be used with any database which supports Word64. +See the haddocks for details.

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

@@ -0,0 +1,53 @@+name:                persistent-spatial+version:             0.1.0.0+synopsis:            Database agnostic, spatially indexed type for geographic points.+description:         Defines type for storing geographic coordinates that can be spatially indexed by any database which supports Word64.+                     This inxeding is implemented using a normal integer index on points repersented using a Morton Z-Order curve.+                     Geographic regions are transformed into a covering set of tiles (contigious ranges) which can be used in a single query.+homepage:            https://github.com/george-steel/persistent-spatial#readme+license:             MIT+license-file:        LICENSE+author:              George Steel+maintainer:          george.steel@gmail.com+copyright:           2019 Satsuma Labs+category:            Database, Geography+build-type:          Simple+cabal-version:       >=1.10+extra-source-files:  README.md+                   , CHANGELOG.md++library+  hs-source-dirs:      src+  ghc-options:         -Wall -Wcompat -Wincomplete-record-updates -Wincomplete-uni-patterns -Wredundant-constraints+  exposed-modules:     Data.Morton+                       Data.LatLong+  build-depends:       base >= 4.7 && < 5+                     , integer-logarithms+                     , aeson+                     , lens+                     , http-api-data+                     , text+                     , persistent+  default-language:    Haskell2010+++test-suite tests+  type:                exitcode-stdio-1.0+  hs-source-dirs:      test+  ghc-options:         -Wall -Wcompat -Wincomplete-record-updates -Wincomplete-uni-patterns -Wredundant-constraints -Wno-unused-imports+  main-is:             Spec.hs+  build-depends:       persistent-spatial+                     , base >= 4.7 && < 5+                     , aeson+                     , http-api-data+                     , text+                     , persistent+                     , hspec+                     , QuickCheck+  default-language:    Haskell2010++++source-repository head+  type:     git+  location: https://github.com/george-steel/persistent-spatial

@@ -0,0 +1,179 @@+{-# LANGUAGE OverloadedStrings, ScopedTypeVariables, ViewPatterns, PatternSynonyms #-}++{-|+Module: Data.LatLong+Description: Spatially indexed type for geographic coordinates.+Copyright: © 2018-2019 Satsuma labs++Defines a type for georgraphic coordinates that can be spatially indexed by any database supporting 64 bit integer values.+This indexing works by reperesenting points using a 'Morton' Z-Order curce, with each coordinate reperesented as a 32-bit fixed-point value+which then have their bits interleaved into a 64-bit integer to the internal reperesentation.++Taking binary prefixes of these values divides the globe into a hierarchy of rectangular tiles (repereseteh here as 'LatLongTile' objects),+each of which is a contiguous interval when points are ordered according to their integer reperesentations.+As any geographic region can be covered by a small number of tiles of simillar size, this provides an easy to loop up data for specific reguions.+Instances and a filter for persistent are provided for this purpose.++-}++module Data.LatLong (+    LatLong(LatLongZ, LatLong), lat, long,+    earthRadius, geoDistance, geoSquare,+    -- * Tiles+    LatLongTile, latLongTileInterval,+    latLongTileCover, latLongTileCoverSquare, tileSetElem, withinTileSet+) where++import Data.Morton+import Data.Aeson+import Data.Proxy+import qualified Data.Text as T+import Control.Monad+import Control.Lens (Lens')+import Data.Word+import Numeric+import Web.HttpApiData+import Database.Persist.Sql++-- | Type for storing geographic coordinates that can be spatially indexed ('Morton' ordering).+-- Each coordinate is reperesented as as 32-bit fixed point value and is also accessible as a Double through a pattern synonym.+-- Order follows a Morton Z-order curve which can be used to search a database by tiles.+-- This works with any database capable of storing and indexing 'Word64' (although this type only uses those values fitting in a 64 bit signed integer)+newtype LatLong =+    -- | Underlying reperesentation and source of ordering for indexing+    LatLongZ Morton+    deriving (Eq, Ord)++two32f :: Double+two32f = 2 ^ (32 :: Int)+++makeLatLong :: Double -> Double -> LatLong+makeLatLong theta phi = LatLongZ (MortonPair theta' phi') where+    theta' = clampLat . floor $ (theta + 90) / 360 * two32f+    phi' = wrapLong . floor $ (phi + 180) / 360 * two32f+    clampLat (x::Int) = fromIntegral (max 0 (min 0x7fffffff x))+    wrapLong (x::Int) = fromIntegral (x `mod` 0x100000000)++latLongCoords :: LatLong -> (Double, Double)+latLongCoords (LatLongZ (MortonPair theta' phi')) = (theta,phi) where+    theta = (fromIntegral theta' * 360 / two32f) - 90+    phi = (fromIntegral phi' * 360 / two32f) - 180++-- | Pattern for accessing latitide and longitude coordinates as 'Double' values.+-- This is not fully isomoprphic as latitude is clipped to ±90, longitude is wrapped mod 360 ±180,+-- and rounding error exists due to the internal fixed-point reperesentation.+pattern LatLong :: Double -> Double -> LatLong+pattern LatLong theta phi <- (latLongCoords -> (theta,phi)) where+    LatLong theta phi = makeLatLong theta phi+{-# COMPLETE LatLong #-}+++instance ToJSON LatLong where+    toJSON (LatLong theta phi) = object ["lat" .= theta, "long" .= phi]+    toEncoding (LatLong theta phi) = pairs $ "lat" .= theta <> "long" .= phi+instance FromJSON LatLong where+    parseJSON = withObject "LatLong" $ \o -> do+        theta <- o .: "lat"+        phi <- o .: "long"+        guard $ theta > -90 && theta < 90+        guard $ phi >= -180 && phi < 180+        return $ LatLong theta phi++instance FromHttpApiData LatLong where+    parseUrlPiece s = maybe (Left "malformed coordinate pair") Right $ do+        [theta',phi'] <- return $ T.splitOn "," s+        Right theta <- return $ parseUrlPiece theta'+        Right phi <- return $ parseUrlPiece phi'+        guard $ theta >= -90 && theta <= 90+        guard $ phi >= -180 && phi <= 180+        return $ LatLong theta phi+instance ToHttpApiData LatLong where+    toUrlPiece (LatLong theta phi) = toUrlPiece theta <> "," <> toUrlPiece phi++instance Show LatLong where+    show (LatLong theta phi) = join+        [ showFFloat (Just 5) (abs theta) []+        , if theta >= 0 then " N " else " S "+        , showFFloat (Just 5) (abs phi) []+        , if phi >= 0 then " E" else " W" ]++instance PersistField LatLong where+    toPersistValue (LatLongZ (Morton x)) = toPersistValue x+    fromPersistValue = fmap (LatLongZ . Morton) . fromPersistValue+instance PersistFieldSql LatLong where+    sqlType _ = sqlType (Proxy :: Proxy Word64)+++-- | Lens for latitude.+lat :: Lens' LatLong Double+lat f (LatLong theta phi) = fmap (\theta' -> LatLong theta' phi) (f theta)+-- | Lens for longitude.+long :: Lens' LatLong Double+long f (LatLong theta phi) = fmap (\phi' -> LatLong theta phi') (f phi)++-- | Earth's average radius in meters+earthRadius :: Double+earthRadius = 6371.2e3++rads :: Double->Double+rads x = x / 180 * pi+degs :: Double->Double+degs x = x * 180 / pi+sindeg :: Double->Double+sindeg = sin . rads+cosdeg :: Double->Double+cosdeg = cos . rads++-- | Calculate distance between two points using the Haversine formula (up to 0.5% due to the assumption of a spherical Earth).+-- Distance is returned in meters.+geoDistance :: LatLong -> LatLong -> Double+geoDistance (LatLong theta1 phi1) (LatLong theta2 phi2) = earthRadius * sigma where+    havdelta x y = sindeg ((x-y)/2) ^ (2::Int)+    hav = havdelta theta1 theta2 + (cosdeg theta1 * cosdeg theta2 * havdelta phi1 phi2)+    sigma = 2 * asin (sqrt hav)++-- | Calculates the corner coordinates of a square with a given center and radius (in meters).+-- Based on the Mercator projection thus has distortion near the poles+-- (within 5% for a radius at most 200km and latitude within ±70).+geoSquare :: LatLong -> Double -> (LatLong, LatLong)+geoSquare (LatLong theta phi) r = let+    dtheta = degs (r / earthRadius)+    dphi = degs (r / earthRadius / cosdeg theta)+    se = LatLong (theta - dtheta) (phi - dphi)+    nw = LatLong (theta + dtheta) (phi + dphi)+    in (se,nw)++++++-- | Represents a LatLong tile, which is both a rectangle and a contoguous interval in the ordering.+newtype LatLongTile = LatLongTile MortonTile deriving (Eq, Read, Show)++-- | Gets the corners of a tile, which are also the bounds of its interval in sort order.+latLongTileInterval :: LatLongTile -> Interval LatLong+latLongTileInterval (LatLongTile t) = fmap LatLongZ (mortonTileBounds t)++-- | Covers a rectangle (defined by its corners) tiles of at most its size.+latLongTileCover :: LatLong -> LatLong -> [LatLongTile]+latLongTileCover se nw = let+    LatLong s _ = se+    LatLong n _ = nw+    LatLongZ y = se+    LatLongZ x = nw+    in if s > n then [] else fmap LatLongTile (mortonTileCoverTorus x y)++-- | Covers a square (defined by its center and radius) by tiles.+latLongTileCoverSquare :: LatLong -> Double -> [LatLongTile]+latLongTileCoverSquare c r = uncurry latLongTileCover $ geoSquare c r++-- | Tests whether a point is contasined in a tile set.+tileSetElem :: LatLong -> [LatLongTile] -> Bool+tileSetElem p ts = or [intervalElem p (latLongTileInterval t) | t <- ts]++-- | Persistent filter producing the SQL equiveland ot 'tileSetElem'.+withinTileSet :: (EntityField row LatLong) -> [LatLongTile] -> Filter row+withinTileSet field tiles = let+    tfilter tile = let Interval a b = latLongTileInterval tile in FilterAnd [field >=. a, field <=. b]+    in FilterOr $ fmap tfilter tiles

@@ -0,0 +1,224 @@+{-# LANGUAGE BangPatterns, ScopedTypeVariables, GeneralizedNewtypeDeriving, DeriveFunctor, ViewPatterns, PatternSynonyms #-}++{- |+Module: Data.Morton+Description: Morton reperesention of integer pairs+Copyright: © 2018-2019 Satsuma labs++Morton reperesentation of integer pairs (interleaved bits) used for creating spatial indexes.++Bit interleaving code is originally from the documentation of Data.Sparse by Edward Kmett at+<https://www.schoolofhaskell.com/user/edwardk/revisiting-matrix-multiplication/part-1>+-}++module Data.Morton (+    Morton(Morton,MortonPair),+    -- * Intervals+    Interval(..), intersectInterval, intervalElem, intervalSize, intervalSizeMorton,+    -- * Rectangles+    MortonRect(MortonRect,MortonRectSides), mortonRectBounds, intersectMortonRect, mortonRectSize,+    -- * Tiles+    MortonTile(..), mortonTileBounds, mortonTileRect, enclosingMortonTile, splitMortonTile, trimMortonTile,+    mortonTileCoverSized, mortonTileCover, mortonTileCoverTorus,+) where++import Data.Bits+import Data.Word+import Data.Monoid ((<>))+import Data.Maybe+import Data.Ord+import Control.Monad+import Numeric+import Text.Read+import Text.Read.Lex+import Text.ParserCombinators.ReadPrec (readP_to_Prec)+import Text.ParserCombinators.ReadP+import Math.NumberTheory.Logarithms++zeropad :: Int -> String -> String+zeropad n s = replicate (n - length s) '0' ++ s++-- | Type implementing a Morton Z-Order Curve.+-- Stores two 'Word32' values with bits interleaved.+-- This allows for spatial indexing by rectangular tiles which form contiguous intervals.+newtype Morton = Morton Word64 deriving (Eq, Ord, Enum)+-- Shows value in hex.+instance (Show Morton) where+    show (Morton m) = "Z" <> zeropad 16 (showHex m [])+instance (Read Morton) where+    readPrec = readP_to_Prec (const readMorton)++readMorton :: ReadP Morton+readMorton = char 'Z' >> fmap Morton readHexP+++-- interleaves the bits of two integers by performing AH,AL,BH,BL -> AH,BH,AL,BL then recursing on H and L portions in SIMD fashion+interleaveM :: Word32 -> Word32 -> Morton+interleaveM !x !y = Morton k5 where+    k0 = unsafeShiftL (fromIntegral x) 32 .|. fromIntegral y+    k1 = unsafeShiftL (k0 .&. 0x00000000FFFF0000) 16  .|. unsafeShiftR k0 16 .&. 0x00000000FFFF0000  .|. k0 .&. 0xFFFF00000000FFFF+    k2 = unsafeShiftL (k1 .&. 0x0000FF000000FF00) 8   .|. unsafeShiftR k1 8  .&. 0x0000FF000000FF00  .|. k1 .&. 0xFF0000FFFF0000FF+    k3 = unsafeShiftL (k2 .&. 0x00F000F000F000F0) 4   .|. unsafeShiftR k2 4  .&. 0x00F000F000F000F0  .|. k2 .&. 0xF00FF00FF00FF00F+    k4 = unsafeShiftL (k3 .&. 0x0C0C0C0C0C0C0C0C) 2   .|. unsafeShiftR k3 2  .&. 0x0C0C0C0C0C0C0C0C  .|. k3 .&. 0xC3C3C3C3C3C3C3C3+    k5 = unsafeShiftL (k4 .&. 0x2222222222222222) 1   .|. unsafeShiftR k4 1  .&. 0x2222222222222222  .|. k4 .&. 0x9999999999999999+++uninterleaveM :: Morton -> (Word32,Word32)+uninterleaveM (Morton k0) = (fromIntegral (unsafeShiftR k5 32), fromIntegral (k5 .&. 0x00000000FFFFFFFF)) where+    k5 = unsafeShiftL (k4 .&. 0x00000000FFFF0000) 16  .|. unsafeShiftR k4 16 .&. 0x00000000FFFF0000  .|. k4 .&. 0xFFFF00000000FFFF+    k4 = unsafeShiftL (k3 .&. 0x0000FF000000FF00) 8   .|. unsafeShiftR k3 8  .&. 0x0000FF000000FF00  .|. k3 .&. 0xFF0000FFFF0000FF+    k3 = unsafeShiftL (k2 .&. 0x00F000F000F000F0) 4   .|. unsafeShiftR k2 4  .&. 0x00F000F000F000F0  .|. k2 .&. 0xF00FF00FF00FF00F+    k2 = unsafeShiftL (k1 .&. 0x0C0C0C0C0C0C0C0C) 2   .|. unsafeShiftR k1 2  .&. 0x0C0C0C0C0C0C0C0C  .|. k1 .&. 0xC3C3C3C3C3C3C3C3+    k1 = unsafeShiftL (k0 .&. 0x2222222222222222) 1   .|. unsafeShiftR k0 1  .&. 0x2222222222222222  .|. k0 .&. 0x9999999999999999++-- | Construct a Morton value from its two coordinates.+pattern MortonPair :: Word32 -> Word32 -> Morton+pattern MortonPair x y <- (uninterleaveM -> (x,y)) where+    MortonPair x y = interleaveM x y+{-# COMPLETE MortonPair #-}+++-- | Type for closed intervals. The second field should be greater than the first.+data Interval a = Interval a a deriving (Eq, Show, Read, Functor)++-- | Returns intersection of two intervals, or Nothing if they do not overlap.+intersectInterval :: Ord a => Interval a -> Interval a -> Maybe (Interval a)+intersectInterval (Interval a b) (Interval a' b') = do+    guard (a <= b)+    guard (a' <= b')+    let x = max a a'+        y = min b b'+    guard (x <= y)+    return $ Interval x y++-- | Tests whether an element is contained within a given Interval.+intervalElem :: (Ord a) => a -> Interval a -> Bool+intervalElem x (Interval a b) = a <= x && x <= b++-- | Returns the size of an integer interval.+intervalSize :: (Integral a) => Interval a -> Integer+intervalSize (Interval a b) = fromIntegral b - fromIntegral a + 1++-- | Returns the size of a Morton interval.+intervalSizeMorton :: Interval Morton -> Integer+intervalSizeMorton (Interval (Morton a) (Morton b)) = fromIntegral b - fromIntegral a + 1+++-- | Type for retangles in Morton space reperesented by upper-left and lower-right corners+data MortonRect = MortonRect {-# UNPACK #-} !Morton {-# UNPACK #-} !Morton deriving (Eq, Show, Read)++-- | Returns x,y bounds of a rectangle+mortonRectBounds :: MortonRect -> (Interval Word32, Interval Word32)+mortonRectBounds (MortonRect (MortonPair x y) (MortonPair x' y')) = (Interval x x', Interval y y')++-- | Construct/match rectangles by their sides.+pattern MortonRectSides :: Interval Word32 -> Interval Word32 -> MortonRect+pattern MortonRectSides xs ys <- (mortonRectBounds -> (xs,ys)) where+    MortonRectSides (Interval x x') (Interval y y') = MortonRect (MortonPair x y) (MortonPair x' y')+{-# COMPLETE MortonRectSides #-}++-- | Rerurns intersection of two rectangles.+intersectMortonRect :: MortonRect -> MortonRect -> Maybe MortonRect+intersectMortonRect (MortonRect a1 b1) (MortonRect a2 b2) = mint where+    selx (Morton m) = m .&. 0xAAAAAAAAAAAAAAAA+    sely (Morton m) = m .&. 0x5555555555555555+    -- interleaving bits with 0 preserves ordering so bit shifts are not necesary for comparisons+    ax = max (selx a1) (selx a2)+    ay = max (sely a1) (sely a2)+    bx = min (selx b1) (selx b2)+    by = min (sely b1) (sely b2)+    mint = if (ax <= bx) && (ay <= by) then Just (MortonRect (Morton $ ax .|. ay) (Morton $ bx .|. by)) else Nothing++-- | Returns the area of a rectangle+mortonRectSize :: MortonRect -> Integer+mortonRectSize (MortonRect (MortonPair ax ay) (MortonPair bx by)) =+    (fromIntegral bx - fromIntegral ax + 1) * (fromIntegral by - fromIntegral ay + 1)+++-- | Type for a tile in Morton space, a special type of rectangle which is the set of all points sharing a common binary prefex.+-- Reperesented as a point and mask length simillarly to a CIDR subnet.+data MortonTile = MortonTile {-# UNPACK #-} !Morton {-# UNPACK #-} !Int+instance (Show MortonTile) where+    show t@(MortonTile _ n) = let Interval m _ = mortonTileBounds t in show m <> "/" <> show n+instance (Read MortonTile) where+    readPrec = readP_to_Prec . const $ do+        m <- readMorton+        char '/'+        n <- readDecP+        guard (n >= 0 && n <= 64)+        return $ MortonTile m n++-- | Values which reperesent the same tile compare equal even if the reperesentative points differ.+instance (Eq MortonTile) where+    a == b = mortonTileBounds a == mortonTileBounds b++-- | A tile sorts immediately before its subtiles, i.e. x sorts before 0 and 1.+instance (Ord MortonTile) where+    compare = comparing (\x -> let Interval a b = mortonTileBounds x in (a, Down b))++-- | Returns a tile as an 'Interval'.+mortonTileBounds :: MortonTile -> Interval Morton+mortonTileBounds (MortonTile (Morton x) n) = let+    mask = shiftR 0xFFFFFFFFFFFFFFFF n+    a = x .&. complement mask+    b = x .|. mask+    in Interval (Morton a) (Morton b)++-- | Returns a tile as a 'MortonRect'.+mortonTileRect :: MortonTile -> MortonRect+mortonTileRect t = let+    Interval a b = mortonTileBounds t+    in MortonRect a b++-- | Finds the smallest tile completely enclosing a rectangle.+--  This can be arbitrarily large if the rectangle crosses a seam.+enclosingMortonTile :: MortonRect -> MortonTile+enclosingMortonTile (MortonRect (Morton a) (Morton b)) =+    let n = countLeadingZeros $ a `xor` b in MortonTile (Morton a) n++-- | Splits a 'MortonTile' in half. Does not split tiles containing a single value.+splitMortonTile :: MortonTile -> [MortonTile]+splitMortonTile t@(MortonTile _ 64) = [t]+splitMortonTile (MortonTile (Morton m) n) = let+    mask = shiftR 0xFFFFFFFFFFFFFFFF n+    low = m .&. complement mask+    high = m .|. mask+    in [MortonTile (Morton low) (n+1), MortonTile (Morton high) (n+1)]++-- | Trims a @MortonTile@ to its subtile overlapping a given rectangle.+-- Returns 'Nothing' if the rectabgle and tile do not intersect.+trimMortonTile :: MortonRect -> MortonTile -> Maybe MortonTile+trimMortonTile rect t = let+    Interval a b = mortonTileBounds t+    trect = MortonRect a b+    irect = intersectMortonRect rect trect+    in fmap enclosingMortonTile irect++-- | Covers a rectangle using tiles within a range of sizes (specified by their mask values).+mortonTileCoverSized :: Int -> Maybe Int -> MortonRect -> [MortonTile]+mortonTileCoverSized big small rect = let+    expandTile tile@(MortonTile p m) = case small of+        Just s | m > s -> MortonTile p s+        _              -> tile+    subdiv tile@(MortonTile _ m)+        | m >= big = [expandTile tile]+        | otherwise = (>>= subdiv) . mapMaybe (trimMortonTile rect) . splitMortonTile $ tile+    in subdiv (enclosingMortonTile rect)++-- | Covers a rectangle with tiles no larger then the area to be covered (no lower size limit).+-- The total area coverd by these tiles bas a trivial upper bound of 8 tiles the rectangle's area plus the area of its enclosing square+-- and the actual performance is usually significantly better (possibly always, although I have not proven so).+mortonTileCover :: MortonRect -> [MortonTile]+mortonTileCover rect = let+    big = max 0 $ 64 - integerLog2 (mortonRectSize rect)+    in mortonTileCoverSized big Nothing rect++-- | Version of 'mortonTileCover' which allows the rectangle to wrap around the maximum x/y coordinates (as if the space were a torus).+mortonTileCoverTorus :: Morton -> Morton -> [MortonTile]+mortonTileCoverTorus (MortonPair ax ay) (MortonPair bx by) = let+    rsize = (fromIntegral (bx - ax) + 1) * (fromIntegral (by - ay) + 1) -- integer overflows handle wraparound case+    big = max 0 $ 64 - integerLog2 rsize+    initrects = MortonRectSides+        <$> (if bx >= ax then [Interval ax bx] else [Interval ax maxBound, Interval minBound bx])+        <*> (if by >= ay then [Interval ay by] else [Interval ay maxBound, Interval minBound by])+    in initrects >>= mortonTileCoverSized big Nothing

@@ -0,0 +1,159 @@+{-# OPTIONS_GHC -Wno-orphans #-}+{-# LANGUAGE ScopedTypeVariables #-}++module Main where++import Test.Hspec+import Test.QuickCheck+import Test.Hspec.QuickCheck+import Data.Word+import Data.List+import Data.Maybe+import Data.Aeson+import Text.Read+import Numeric+import Web.HttpApiData+import Database.Persist++import Data.Morton+import Data.LatLong++instance (Arbitrary Morton) where+    arbitrary = fmap Morton arbitrary+instance (Arbitrary a, Ord a) => (Arbitrary (Interval a)) where+    arbitrary = do+        x <- arbitrary+        y <- arbitrary+        return $ Interval (min x y) (max x y)+instance (Arbitrary MortonRect) where+    arbitrary = MortonRectSides <$> arbitrary <*> arbitrary++prop_morton_isom :: Word32 -> Word32 -> Bool+prop_morton_isom x y = let+    z = MortonPair x y+    MortonPair a b = z+    in (x == a) && (y == b)++prop_morton_isom_rev :: Morton -> Bool+prop_morton_isom_rev z = let+    MortonPair x y = z+    z' = MortonPair x y+    in z == z'++prop_morton_parse :: Morton -> Bool+prop_morton_parse x = readMaybe (show x) == Just x++prop_morton_intersect :: Interval Word32 -> Interval Word32 -> Interval Word32 -> Interval Word32 -> Property+prop_morton_intersect a b c d = let+    rx = MortonRectSides a b+    ry = MortonRectSides c d+    rz = intersectMortonRect rx ry+    rz' = MortonRectSides <$> (intersectInterval a c) <*> (intersectInterval b d)+    in classify (isJust rz) "overlapping" $ rz == rz'++prop_morton_tile_parse :: Morton -> Property+prop_morton_tile_parse m = forAll (choose (0,64)) $ \n -> let+    t = MortonTile m n+    mt = readMaybe (show t)+    in mt == Just t++prop_morton_tile_bounds :: MortonRect -> Bool+prop_morton_tile_bounds rect = let+    t@(MortonTile _ n) = enclosingMortonTile rect+    s = intervalSizeMorton (mortonTileBounds t)+    sr = mortonRectSize (mortonTileRect t)+    s' = 2 ^ (64 - n)+    in s == s' && sr == s'+++mortonRectAspect :: MortonRect -> Double+mortonRectAspect (MortonRectSides ix iy) = let+    dx = intervalSize ix+    dy = intervalSize iy+    in fromIntegral (max dx dy) / fromIntegral (min dx dy)++mortonRectExpansion :: MortonRect -> Integer+mortonRectExpansion rect = let+    tiles = mortonTileCover rect+    rsize = mortonRectSize rect+    tsize = sum . fmap (intervalSizeMorton . mortonTileBounds) $ tiles+    in (tsize + 1) `div` rsize++prop_morton_split_8 :: MortonRect -> Property+prop_morton_split_8 rect = let+    ratio = mortonRectExpansion rect+    in collect ratio $ mortonRectAspect rect < 8 ==> (ratio >= 1 && ratio <= 16)++badSquare = MortonRect (Morton 0x0fffFfffFfffFfff) (Morton 0xc000000000000000)++++instance (Arbitrary LatLong) where+    arbitrary = LatLong <$> choose (-90,90) <*> choose (-180,180)++prop_latlong_json :: LatLong -> Bool+prop_latlong_json p = let+    mp = decode (encode p)+    in mp == Just p++prop_latlong_http :: LatLong -> Bool+prop_latlong_http p = let+    mp = parseUrlPiece (toUrlPiece p)+    in mp == Right p++prop_latlong_persist :: LatLong -> Bool+prop_latlong_persist p = let+    mp = fromPersistValue (toPersistValue p)+    in mp == Right p++prop_geo_triangle :: LatLong -> LatLong -> LatLong -> Bool+prop_geo_triangle a b c = let+    ttest x y z = geoDistance x y + geoDistance y z >= geoDistance x z+    in ttest a b c && ttest b c a && ttest c a b++pctError :: Double -> Double -> Double+pctError ref samp = abs (samp - ref) / ref++prop_square_corner_dist :: Property+prop_square_corner_dist = let+    gen = (,) <$> choose (10,200000) <*> (LatLong <$> choose (-70,70) <*> choose (-180,180))+    test (r,p) = let+        d = 2*r+        (LatLong s w, LatLong n e) = geoSquare p r+        laterr =  max (pctError d (geoDistance (LatLong n w) (LatLong s w)))+                      (pctError d (geoDistance (LatLong n e) (LatLong s e)))+        --laterrbin :: Int = max (-9) (ceiling (logBase 10 laterr))+        longerr = max (pctError d (geoDistance (LatLong n w) (LatLong n e)))+                      (pctError d (geoDistance (LatLong s w) (LatLong s e)))+        longerrbin = minimum [x | x <- [100,10,5,2,1,0.5,0.1,0.01], x > longerr * 100]+        in collect longerrbin . counterexample (show (laterr,longerr)) $ laterr < 1e-3 && longerr < 0.1+    in forAll gen test++prop_tiles_cover_points :: Property+prop_tiles_cover_points = let+    gen = (,) <$> choose (10,1000000) <*> (LatLong <$> choose (-70,70) <*> choose (-180,180))+    test (r,p) = let+        (sw@(LatLong s w), ne@(LatLong n e)) = geoSquare p r+        tiles = latLongTileCover sw ne+        pgen = LatLong <$> choose (s,n) <*> choose (w,e)+        in forAll pgen (`tileSetElem` tiles)+    in forAll gen test++main :: IO ()+main = hspec . modifyMaxSuccess (const 10000) $ do+    describe "Data.Morton" $ do+        prop "interleaving is reversible" prop_morton_isom+        prop "deinterleaving is reversible" prop_morton_isom_rev+        prop "point Read instance works" prop_morton_parse+        prop "rectangle intersection mathes interval intersection" prop_morton_intersect+        prop "tile Read instance works" prop_morton_tile_parse+        prop "tile size (both definitions) matches mask value" prop_morton_tile_bounds+        it   "pathological square expands at most 6-fold" $ mortonRectExpansion badSquare `shouldSatisfy` (< 6)+        modifyMaxSuccess (const 100000) $ prop "rectangle expansion tile cover (eccentricity limit 8, floored expansion collected)" prop_morton_split_8+    describe "Data.LatLong" $ do+        prop "Aeson instances work" prop_latlong_json+        prop "HttpApiData instances work" prop_latlong_http+        prop "PersistField instance works" prop_latlong_persist+        prop "geoDistance obeys triangle inequality" prop_geo_triangle+        prop "geoSquare corners within tolerance (% error collected), 10m < r < 200km, 70S < lat < 70N" prop_square_corner_dist+        prop "tile covers contain all points in square" prop_tiles_cover_points