packages feed

numhask-space 0.4.0 → 0.13.3.0

raw patch · 12 files changed

Files

+ ChangeLog.md view
@@ -0,0 +1,9 @@+# Changelog++## 0.13.3.0++- Fixed `Show`/`Read` mismatch for `Point`, `Rect` and `Range` by adding+  `showsPrec` implementations that parenthesise values when precedence is+  greater than function-application precedence.+- Added `Read` instances for `Point` and `Rect`.+- Relaxed the upper bound on `time` to `<1.17`.
LICENSE view
@@ -1,4 +1,4 @@-Copyright Tony Day (c) 2016+Copyright (c) 2016, Tony Day  All rights reserved.
− Setup.hs 
@@ -1,2 +0,0 @@-import Distribution.Simple-main = defaultMain
numhask-space.cabal view
@@ -1,79 +1,89 @@ cabal-version: 3.0 name: numhask-space-version: 0.4.0-synopsis:-  numerical spaces+version: 0.13.3.0+license: BSD-3-Clause+license-file: LICENSE+copyright: Tony Day (c) 2016+category: math+author: Tony Day+maintainer: tonyday567@gmail.com+homepage: https://github.com/tonyday567/numhask-space#readme+bug-reports: https://github.com/tonyday567/numhask-space/issues+synopsis: Numerical spaces. description:-  Spaces and the numerical elements that inhabit them.-category:-  mathematics-homepage:-  https://github.com/tonyday567/numhask-space#readme-bug-reports:-  https://github.com/tonyday567/numhask-space/issues-author:-  Tony Day-maintainer:-  tonyday567@gmail.com-copyright:-  Tony Day-license:-  BSD-3-Clause-license-file:-  LICENSE-build-type:-  Simple+  @numhask-space@ provides support for spaces where [space](https://en.wikipedia.org/wiki/Space_(mathematics\)) is defined as a set of numbers with a lower and upper bound.++  == Usage++  >>> {-# LANGUAGE RebindableSyntax #-}+  >>> import NumHask.Prelude+  >>> import NumHask.Space++build-type: Simple+tested-with:+  ghc ==9.10.2+  ghc ==9.12.2+  ghc ==9.14.1+  ghc ==9.6.7+  ghc ==9.8.4++extra-doc-files: ChangeLog.md+ source-repository head-  type:-    git-  location:-    https://github.com/tonyday567/numhask-space+  type: git+  location: https://github.com/tonyday567/numhask-space -library-  hs-source-dirs:-    src+common ghc-options-stanza   ghc-options:     -Wall     -Wcompat+    -Widentities     -Wincomplete-record-updates     -Wincomplete-uni-patterns+    -Wpartial-fields     -Wredundant-constraints++common ghc2024-additions   default-extensions:-    NegativeLiterals-    NoImplicitPrelude-    OverloadedStrings-    UnicodeSyntax+    DataKinds+    DerivingStrategies+    DisambiguateRecordFields+    ExplicitNamespaces+    GADTs+    LambdaCase+    MonoLocalBinds+    RoleAnnotations++common ghc2024-stanza+  if impl(ghc >=9.10)+    default-language:+      GHC2024+  else+    import: ghc2024-additions+    default-language:+      GHC2021++library+  import: ghc-options-stanza+  import: ghc2024-stanza+  hs-source-dirs: src   build-depends:-      base >=4.7 && <5-    , adjunctions >=4.0 && <5-    , semigroupoids >=5 && <6-    , distributive >=0.2.2 && <1-    , lattices >= 2.0.1 && <2.1-    , time >= 1.8.0.2 && <2-    , text >= 1.2.3.1 && <2-    , foldl >= 1.4.5 && <2-    , containers >= 0.6 && < 0.7-    , tdigest >= 0.2.1 && < 0.3-    , protolude >= 0.3.0 && < 0.4.0+    adjunctions >=4.0 && <5,+    base >=4.14 && <5,+    containers >=0.6 && <0.9,+    distributive >=0.4 && <0.7,+    numhask >=0.10 && <0.14,+    semigroupoids >=5.3 && <6.1,+    tdigest >=0.2.1 && <0.4,+    text >=1.2 && <2.2,+    time >=1.9.1 && <1.17,+    vector >=0.12.3 && <0.14,+   exposed-modules:     NumHask.Space-    NumHask.Space.Types+    NumHask.Space.Histogram+    NumHask.Space.Point     NumHask.Space.Range     NumHask.Space.Rect-    NumHask.Space.Point     NumHask.Space.Time-    NumHask.Space.Histogram-  other-modules:-  default-language: Haskell2010--test-suite test-  type: exitcode-stdio-1.0-  main-is: test.hs-  hs-source-dirs:-      test-  ghc-options: -Wall -Wcompat -Wincomplete-record-updates -Wincomplete-uni-patterns -Wredundant-constraints-  build-depends:-      base >=4.7 && <5-    , doctest-  default-language: Haskell2010-+    NumHask.Space.Types
src/NumHask/Space.hs view
@@ -1,18 +1,27 @@-{-# OPTIONS_GHC -Wall #-}+{-# LANGUAGE RebindableSyntax #-} --- | A continuous set of numbers.------ Mathematics does not define a space, leaving library devs to experiment.+-- | Mathematics does not rigorously define a [space](https://en.wikipedia.org/wiki/Space_(mathematics\)), leaving library devs free to explore. ----- https://en.wikipedia.org/wiki/Space_(mathematics)+-- “But who can quantify+--  the algebra of space,+--  or weigh those worlds that swim+--  each in its place?+--  Who can outdo the dark?+--  And what computer knows+--  how beauty comes to birth -+--  shell star and rose? --+-- ~ Technicians by Jean Kenward”+-- ~ John Foster module NumHask.Space-  ( -- * Space+  ( -- * Usage+    -- $setup++    -- * Space     -- $space     module NumHask.Space.Types,      -- * Instances-    -- $instances     module NumHask.Space.Point,     module NumHask.Space.Range,     module NumHask.Space.Rect,@@ -21,23 +30,36 @@   ) where +import NumHask.Space.Histogram hiding () import NumHask.Space.Point hiding () import NumHask.Space.Range hiding () import NumHask.Space.Rect hiding () import NumHask.Space.Time hiding ()-import NumHask.Space.Histogram hiding () import NumHask.Space.Types hiding () --- $space--- The final frontier.+-- $setup+--+-- >>> :m -Prelude+-- >>> :set -XRebindableSyntax+-- >>> import NumHask.Prelude+-- >>> import NumHask.Space+-- >>> Point 1 1+-- Point 1 1+--+-- >>> one :: Range Double+-- Range -0.5 0.5+--+-- >>> grid OuterPos (Range 0 50 :: Range Double) 5+-- [0.0,10.0,20.0,30.0,40.0,50.0] --- $instances+-- $space+-- -- Space is an interesting cross-section of many programming domains. ----- - A Range is a Space of numbers.+-- - A 'Range' is a 'Space' of numbers. ----- - A Rect is a Space of Points.+-- - A 'Rect' is a 'Space' of 'Point's. It can also be a 'Space' of 'Rect's (but this is not yet coded up here). ----- - A time span is a space containing moments of time.+-- - A time span is a 'Space' containing moments of time. ----- - A histogram is a divided Range with a count of elements within each division.+-- - A 'Histogram' is a divided 'Range' with a count of elements within each division.
src/NumHask/Space/Histogram.hs view
@@ -1,9 +1,10 @@-{-# LANGUAGE DataKinds #-}-{-# OPTIONS_GHC -Wall #-}+{-# LANGUAGE RebindableSyntax #-}+{-# OPTIONS_GHC -fno-warn-name-shadowing #-} --- | A histogram, if you squint, is a series of contiguous ranges, annotated with values.+-- | A histogram, if you squint, is a series of contiguous 'Range's, annotated with values. module NumHask.Space.Histogram   ( Histogram (..),+    emptyHistogram,     DealOvers (..),     fill,     cutI,@@ -11,119 +12,160 @@     makeRects,     regularQuantiles,     quantileFold,-    fromQuantiles,     freq,+    average,+    quantiles,+    quantile,   ) where -import qualified Control.Foldl as L-import Data.Bool-import Data.Foldable-import qualified Data.List-import qualified Data.Map as Map-import Data.Maybe-import Data.TDigest+import Data.Map qualified as Map+import Data.TDigest qualified as TD+import Data.Vector qualified as V+import NumHask.Prelude import NumHask.Space.Range import NumHask.Space.Rect import NumHask.Space.Types-import Prelude +-- $setup+-- >>> :m -Prelude+-- >>> :set -XRebindableSyntax+-- >>> import NumHask.Prelude+-- >>> import NumHask.Space+ -- | This Histogram is a list of contiguous boundaries (a boundary being the lower edge of one bucket and the upper edge of another), and a value (usually a count) for each bucket, represented here as a map -- -- Overs and Unders are contained in key = 0 and key = length cuts-data Histogram-  = Histogram-      { cuts :: [Double], -- bucket boundaries-        values :: Map.Map Int Double -- bucket counts-      }+-- Intervals are defined as (l,u]+data Histogram = Histogram+  { cuts :: V.Vector Double, -- bucket boundaries+    values :: Map.Map Int Double -- bucket counts+  }   deriving (Show, Eq) +-- | A histogram with no cuts nor data.+emptyHistogram :: Histogram+emptyHistogram = Histogram V.empty Map.empty+ -- | Whether or not to ignore unders and overs.  If overs and unders are dealt with, IncludeOvers supplies an assumed width for the outer buckets. data DealOvers = IgnoreOvers | IncludeOvers Double  -- | Fill a Histogram using pre-specified cuts ----- >>> fill [0,50,100] [1..100]+-- >>> fill [0,50,100] [0..99] -- Histogram {cuts = [0.0,50.0,100.0], values = fromList [(1,50.0),(2,50.0)]} fill :: (Foldable f) => [Double] -> f Double -> Histogram-fill cs xs = Histogram cs (foldl' (\x a -> Map.insertWith (+) (cutI cs a) 1 x) Map.empty xs)+fill cs xs =+  Histogram+    (V.fromList cs)+    (foldl' (\x a -> Map.insertWith (+) (cutI (V.fromList cs) a) 1 x) Map.empty xs)  -- | find the index of the bucket the value is contained in.-cutI :: (Ord a) => [a] -> a -> Int-cutI bs n = go bs 0+cutI :: (Ord a) => V.Vector a -> a -> Int+cutI cs a = go (Range zero (V.length cs))   where-    go [] i = i-    go (x:xs) i = bool i (go xs (i+1)) (n>x)+    go (Range l u) =+      let k = (u + l) `div` 2+       in case compare a (cs V.! k) of+            EQ -> k + 1+            LT -> bool (go (Range l k)) k (l == k)+            GT ->+              bool+                ( case compare a (cs V.! (k + one)) of+                    EQ -> k + 2+                    LT -> k + 1+                    GT -> go (Range k u)+                )+                (k + 1)+                (k >= u - one)  -- | Make a histogram using n equally spaced cuts over the entire range of the data -- -- >>> regular 4 [0..100]--- Histogram {cuts = [0.0,25.0,50.0,75.0,100.0], values = fromList [(0,1.0),(1,25.0),(2,25.0),(3,25.0),(4,25.0)]}+-- Histogram {cuts = [0.0,25.0,50.0,75.0,100.0], values = fromList [(1,25.0),(2,25.0),(3,25.0),(4,25.0),(5,1.0)]} regular :: Int -> [Double] -> Histogram+regular _ [] = emptyHistogram regular n xs = fill cs xs   where-    cs = grid OuterPos (space1 xs :: Range Double) n+    cs = grid OuterPos (unsafeSpace1 xs :: Range Double) n  -- | Transform a Histogram to Rects -- -- >>> makeRects IgnoreOvers (regular 4 [0..100]) -- [Rect 0.0 25.0 0.0 0.25,Rect 25.0 50.0 0.0 0.25,Rect 50.0 75.0 0.0 0.25,Rect 75.0 100.0 0.0 0.25] makeRects :: DealOvers -> Histogram -> [Rect Double]-makeRects o (Histogram cs counts) = Data.List.zipWith4 Rect x z y w'+makeRects o (Histogram cs counts) = V.toList $ V.zipWith3 (\x z w' -> Rect x z zero w') x z w'   where-    y = repeat 0     w =-      zipWith+      V.zipWith         (/)-        ((\x' -> Map.findWithDefault 0 x' counts) <$> [f .. l])-        (zipWith (-) z x)+        ((\x' -> Map.findWithDefault 0 x' counts) <$> V.enumFromN f (l - f + one))+        (V.zipWith (-) z x)     f = case o of-      IgnoreOvers -> 1-      IncludeOvers _ -> 0+      IgnoreOvers -> one+      IncludeOvers _ -> zero     l = case o of-      IgnoreOvers -> length cs - 1+      IgnoreOvers -> length cs - one       IncludeOvers _ -> length cs     w' = (/ sum w) <$> w     x = case o of       IgnoreOvers -> cs       IncludeOvers outw ->-        [Data.List.head cs - outw]+        V.singleton (V.head cs - outw)           <> cs-          <> [Data.List.last cs + outw]-    z = drop 1 x+          <> V.singleton (V.last cs + outw)+    z = V.drop one x  -- | approx regular n-quantiles -- -- >>> regularQuantiles 4 [0..100] -- [0.0,24.75,50.0,75.25,100.0] regularQuantiles :: Double -> [Double] -> [Double]-regularQuantiles n = L.fold (quantileFold qs)+regularQuantiles n xs = quantileFold qs xs   where     qs = ((1 / n) *) <$> [0 .. n]  -- | one-pass approximate quantiles fold-quantileFold :: [Double] -> L.Fold Double [Double]-quantileFold qs = L.Fold step begin done+quantileFold :: [Double] -> [Double] -> [Double]+quantileFold qs xs = done $ foldl' step begin xs   where-    step x a = Data.TDigest.insert a x-    begin = tdigest ([] :: [Double]) :: TDigest 25-    done x = fromMaybe (0 / 0) . (`quantile` compress x) <$> qs+    step x a = TD.insert a x+    begin = TD.tdigest ([] :: [Double]) :: TD.TDigest 25+    done x = fromMaybe (0 / 0) . (`TD.quantile` TD.compress x) <$> qs --- | take a specification of quantiles and make a Histogram+-- | normalize a histogram ----- >>> fromQuantiles [0,0.25,0.5,0.75,1] (regularQuantiles 4 [0..100])--- Histogram {cuts = [0.0,24.75,50.0,75.25,100.0], values = fromList [(1,0.25),(2,0.25),(3,0.25),(4,0.25)]}-fromQuantiles :: [Double] -> [Double] -> Histogram-fromQuantiles qs xs = Histogram xs (Map.fromList $ zip [1 ..] (diffq qs))-  where-    diffq [] = []-    diffq [_] = []-    diffq (x : xs') = L.fold (L.Fold step (x, []) (reverse . snd)) xs'-    step (a0, xs') a = (a, (a - a0) : xs')---- | normalize a histogram so that sum values = one+-- > \h -> sum (values $ freq h) == one ----- >>> freq $ fill [0,50,100] [1..100]+-- >>> freq $ fill [0,50,100] [0..99] -- Histogram {cuts = [0.0,50.0,100.0], values = fromList [(1,0.5),(2,0.5)]} freq :: Histogram -> Histogram freq (Histogram cs vs) = Histogram cs $ Map.map (* recip (sum vs)) vs++-- | average+--+-- >>> average [0..1000]+-- 500.0+average :: (Foldable f) => f Double -> Double+average xs = sum xs / fromIntegral (length xs)++-- | Regularly spaced (approx) quantiles+--+-- >>> quantiles 5 [1..1000]+-- [1.0,200.5,400.5,600.5000000000001,800.5,1000.0]+quantiles :: (Foldable f) => Int -> f Double -> [Double]+quantiles n xs =+  ( \x ->+      fromMaybe 0 $+        TD.quantile x (TD.tdigest xs :: TD.TDigest 25)+  )+    . (/ fromIntegral n)+    . fromIntegral+    <$> [0 .. n]++-- | single (approx) quantile+--+-- >>> quantile 0.1 [1..1000]+-- 100.5+quantile :: (Foldable f) => Double -> f Double -> Double+quantile p xs = fromMaybe 0 $ TD.quantile p (TD.tdigest xs :: TD.TDigest 25)
src/NumHask/Space/Point.hs view
@@ -1,31 +1,43 @@-{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE RebindableSyntax #-} {-# LANGUAGE TypeFamilies #-}-{-# OPTIONS_GHC -Wall #-}  -- | A 2-dimensional point. module NumHask.Space.Point   ( Point (..),-    rotate,+    rotateP,     gridP,+    dotP,+    (<.>),+    crossP,+    flipY,+    Line (..),+    lineSolve,+    lineDistance,+    closestPoint,+    lineIntersect,+    translate,+    scaleT,+    skew,   ) where -import Algebra.Lattice-import Data.Distributive as D+import Data.Data+import Data.Distributive import Data.Functor.Classes import Data.Functor.Rep-import GHC.Generics (Generic)+import NumHask.Prelude hiding (Distributive) import NumHask.Space.Range import NumHask.Space.Types-import Text.Show-import Prelude  -- $setup--- +-- >>> :m -Prelude+-- >>> :set -XRebindableSyntax+-- >>> import NumHask.Prelude+-- >>> import NumHask.Space --- | A 2-dim point of a's+-- | A 2-dimensional Point of a's ----- A Point is functorial over both arguments, and is a Num instance.+-- In contrast with a tuple, a Point is functorial over both arguments. -- -- >>> let p = Point 1 1 -- >>> p + p@@ -33,30 +45,57 @@ -- >>> (2*) <$> p -- Point 2 2 ----- A major reason for this bespoke treatment of a point is that Points do not have maximums and minimums but they form a lattice, and this is useful for folding points to find out the (rectangular) Space they occupy.+-- A major reason for this bespoke treatment (compared to just using linear, say) is that Points do not have maximums and minimums but they do form a lattice, and this is useful for folding sets of points to find out the (rectangular) Space they occupy. -- -- >>> Point 0 1 /\ Point 1 0 -- Point 0 0 -- >>> Point 0 1 \/ Point 1 0 -- Point 1 1-data Point a-  = Point a a-  deriving (Eq, Generic)--instance (Show a) => Show (Point a) where-  show (Point a b) = "Point " <> Text.Show.show a <> " " <> Text.Show.show b--instance Functor Point where-  fmap f (Point a b) = Point (f a) (f b)+--+-- This is used extensively in [chart-svg](https://hackage.haskell.org/package/chart-svg) to ergonomically obtain chart areas.+--+-- > unsafeSpace1 [Point 1 0, Point 0 1] :: Rect Double+-- Rect 0.0 1.0 0.0 1.0+data Point a = Point+  { _x :: a,+    _y :: a+  }+  deriving (Eq, Generic, Data)  instance Eq1 Point where   liftEq f (Point a b) (Point c d) = f a c && f b d -instance Show1 Point where-  liftShowsPrec sp _ d (Point a b) = showsBinaryWith sp sp "Point" d a b+instance (Ord a, Additive a, Show a) => Show (Point a) where+  show (Point a b) = "Point " <> wrap a <> " " <> wrap b+    where+      wrap x = bool (show x) ("(" <> show x <> ")") (x < zero)+  showsPrec d p = showParen (d > app_prec) (showString (show p))+    where+      app_prec = 10 -instance Applicative Point where+instance (Read a) => Read (Point a) where+  readsPrec d s = readParen (d > app_prec) parsePoint s+    where+      app_prec = 10+      parsePoint r = do+        ("Point", r1) <- lex r+        (x, r2) <- readsMaybeNeg r1+        (y, r3) <- readsMaybeNeg r2+        pure (Point x y, r3)+      readsMaybeNeg t = case reads t of+        [(v, rest)] -> [(v, rest)]+        [] -> case t of+          '(' : r0 -> case reads r0 of+            [(v, rest)] -> case rest of+              ')' : r -> [(v, r)]+              _ -> []+            _ -> []+          _ -> [] +instance Functor Point where+  fmap f (Point a b) = Point (f a) (f b)++instance Applicative Point where   pure a = Point a a    (Point fa fb) <*> Point a b = Point (fa a) (fb b)@@ -76,36 +115,28 @@ instance (Semigroup a) => Semigroup (Point a) where   (Point a0 b0) <> (Point a1 b1) = Point (a0 <> a1) (b0 <> b1) -instance (Semigroup a, Monoid a) => Monoid (Point a) where-+instance (Monoid a) => Monoid (Point a) where   mempty = Point mempty mempty    mappend = (<>)  instance (Bounded a) => Bounded (Point a) where-   minBound = Point minBound minBound    maxBound = Point maxBound maxBound -instance (Num a) => Num (Point a) where-+instance (Additive a) => Additive (Point a) where   (Point a0 b0) + (Point a1 b1) = Point (a0 + a1) (b0 + b1)+  zero = Point zero zero +instance (Subtractive a) => Subtractive (Point a) where   negate = fmap negate +instance (Multiplicative a) => Multiplicative (Point a) where   (Point a0 b0) * (Point a1 b1) = Point (a0 * a1) (b0 * b1)--  signum = fmap signum--  abs = fmap abs--  fromInteger x = Point (fromInteger x) (fromInteger x)--instance (Fractional a) => Fractional (Point a) where--  fromRational x = Point (fromRational x) (fromRational x) +  one = Point one one +instance (Divisive a) => Divisive (Point a) where   recip = fmap recip  instance Distributive Point where@@ -114,8 +145,21 @@       getL (Point l _) = l       getR (Point _ r) = r -instance Representable Point where+instance (Additive a) => AdditiveAction (Point a) where+  type AdditiveScalar (Point a) = a+  (|+) (Point x y) a = Point (a + x) (a + y) +instance (Subtractive a) => SubtractiveAction (Point a) where+  (|-) (Point x y) a = Point (x - a) (y - a)++instance (Multiplicative a) => MultiplicativeAction (Point a) where+  type Scalar (Point a) = a+  (|*) (Point x y) a = Point (a * x) (a * y)++instance (Divisive a) => DivisiveAction (Point a) where+  (|/) (Point x y) a = Point (x / a) (y / a)++instance Representable Point where   type Rep Point = Bool    tabulate f = Point (f False) (f True)@@ -123,24 +167,129 @@   index (Point l _) False = l   index (Point _ r) True = r -instance (Ord a) => Lattice (Point a) where-+instance (Ord a) => JoinSemiLattice (Point a) where   (\/) (Point x y) (Point x' y') = Point (max x x') (max y y') +instance (Ord a) => MeetSemiLattice (Point a) where   (/\) (Point x y) (Point x' y') = Point (min x x') (min y y') --- | rotate a point by x degrees relative to the origin------ >>> rotate 90 (Point 0 1)--- Point 1.0 6.123233995736766e-17-rotate :: (Floating a) => a -> Point a -> Point a-rotate d (Point x y) = Point (x * cos d' + y * sin d') (y * cos d' - x * sin d')+instance+  (ExpField a, Eq a) =>+  Basis (Point a)   where-    d' = d * pi / 180+  type Mag (Point a) = a+  type Base (Point a) = Point a +  magnitude (Point x y) = sqrt (x * x + y * y)+  basis p = let m = magnitude p in bool (p |/ m) zero (m == zero)++-- | angle formed by a vector from the origin to a Point and the x-axis (Point 1 0). Note that an angle between two points p1 & p2 is thus angle p2 - angle p1+instance (TrigField a) => Direction (Point a) where+  type Dir (Point a) = a+  angle (Point x y) = atan2 y x+  ray x = Point (cos x) (sin x)++instance (Multiplicative a, Additive a) => Affinity (Point a) a where+  transform (Transform a b c d e f) (Point x y) =+    Point (a * x + b * y + c) (d * x + e * y + f)++-- | move an 'Affinity' by a 'Point'+translate :: (TrigField a) => Point a -> Transform a+translate (Point x y) = Transform one zero x zero one y++-- | scale an 'Affinity' by a 'Point'+scaleT :: (TrigField a) => Point a -> Transform a+scaleT (Point x y) = Transform x zero zero y zero zero++-- | Skew transform+--+-- x-axis skew+--+-- > skew (Point x 0)+skew :: (TrigField a) => Point a -> Transform a+skew (Point x y) = Transform one (tan x) zero (tan y) one zero++-- | rotate a point by x relative to the origin+--+-- >>> rotateP (pi/2) (Point 1 0)+-- Point 6.123233995736766e-17 1.0+rotateP :: (TrigField a) => a -> Point a -> Point a+rotateP d p = rotate d |. p+ -- | Create Points for a formulae y = f(x) across an x range ----- >>> gridP (**2) (Range 0 4) 4+-- >>> gridP (^^2) (Range 0 4) 4 -- [Point 0.0 0.0,Point 1.0 1.0,Point 2.0 4.0,Point 3.0 9.0,Point 4.0 16.0]-gridP :: (Ord a, Fractional a) => (a -> a) -> Range a -> Int -> [Point a]+gridP :: (FieldSpace (Range a)) => (a -> a) -> Range a -> Grid (Range a) -> [Point a] gridP f r g = (\x -> Point x (f x)) <$> grid OuterPos r g++-- | dot product+dotP :: (Multiplicative a, Additive a) => Point a -> Point a -> a+dotP (Point x y) (Point x' y') = x * x' + y * y'++infix 4 <.>++-- | dot product operator+(<.>) :: (Multiplicative a, Additive a) => Point a -> Point a -> a+(<.>) = dotP++-- | cross product+crossP :: (Multiplicative a, Subtractive a) => Point a -> Point a -> a+crossP (Point x y) (Point x' y') = x * y' - y * x'++-- | reflect on x-axis+flipY :: (Subtractive a) => Point a -> Point a+flipY (Point x y) = Point x (-y)++-- | A line is a composed of 2 'Point's+data Line a = Line+  { lineStart :: Point a,+    lineEnd :: Point a+  }+  deriving (Show, Eq, Functor, Foldable, Traversable)++instance (Multiplicative a, Additive a) => Affinity (Line a) a where+  transform t (Line s e) = Line (transform t s) (transform t e)++-- | Return the parameters (a, b, c) for the line equation @a*x + b*y + c = 0@.+lineSolve :: (ExpField a, Eq a) => Line a -> (a, a, a)+lineSolve (Line p1 p2) = (-my, mx, c)+  where+    m@(Point mx my) = basis (p2 - p1)+    c = crossP p1 m++-- | Return the signed distance from a point to the line.  If the+-- distance is negative, the point lies to the right of the line+lineDistance :: (ExpField a) => Line a -> Point a -> a+lineDistance (Line (Point x1 y1) (Point x2 y2)) =+  let dy = y1 - y2+      dx = x2 - x1+      d = sqrt (dx * dx + dy * dy)+   in dy `seq`+        dx `seq`+          d `seq`+            \(Point x y) -> (x - x1) * dy / d + (y - y1) * dx / d++-- | Return the point on the line closest to the given point.+closestPoint :: (Field a) => Line a -> Point a -> Point a+closestPoint (Line p1 p2) p3 = Point px py+  where+    d@(Point dx dy) = p2 - p1+    u = dy * _y p3 + dx * _x p3+    v = _x p1 * _y p2 - _x p2 * _y p1+    m = d <.> d+    px = (dx * u + dy * v) / m+    py = (dy * u - dx * v) / m++-- | Calculate the intersection of two lines.  If the determinant is+-- less than tolerance (parallel or coincident lines), return Nothing.+lineIntersect :: (Ord a, Epsilon a, Absolute a, Field a) => Line a -> Line a -> Maybe (Point a)+lineIntersect (Line p1 p2) (Line p3 p4)+  | abs det <= epsilon = Nothing+  | otherwise = Just $ (a *| d2 - b *| d1) |/ det+  where+    d1 = p1 - p2+    d2 = p3 - p4+    det = crossP d1 d2+    a = crossP p1 p2+    b = crossP p3 p4
src/NumHask/Space/Range.hs view
@@ -1,27 +1,30 @@-{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE RebindableSyntax #-} {-# LANGUAGE TypeFamilies #-}-{-# OPTIONS_GHC -Wall #-}+{-# LANGUAGE UndecidableInstances #-}  -- | A Space containing numerical elements module NumHask.Space.Range   ( Range (..),     gridSensible,+    stepSensible,   ) where -import Algebra.Lattice-import Data.Bool (bool)+import Data.Data import Data.Distributive as D import Data.Functor.Apply (Apply (..)) import Data.Functor.Classes import Data.Functor.Rep-import Data.Semigroup.Foldable (Foldable1 (..))-import Data.Semigroup.Traversable (Traversable1 (..))-import GHC.Generics (Generic)+import GHC.Show (show)+import NumHask.Prelude hiding (show) import NumHask.Space.Types as S-import Prelude  -- $setup+--+-- >>> :m -Prelude+-- >>> :set -XFlexibleContexts+-- >>> import NumHask.Prelude+-- >>> import NumHask.Space  -- | A continuous range over type a --@@ -29,12 +32,12 @@ -- >>> a -- Range -1 1 ----- Num instance based on interval arithmetic (with Ranges normalising to lower ... upper)--- -- >>> a + a -- Range -2 2+-- -- >>> a * a--- Range -1 1+-- Range -2.0 2.0+-- -- >>> (+1) <$> (Range 1 2) -- Range 2 3 --@@ -46,24 +49,24 @@ -- -- Create an equally spaced grid including outer bounds over a Range ----- >>> grid OuterPos (Range 0 10) 5+-- >>> grid OuterPos (Range 0.0 10.0) 5 -- [0.0,2.0,4.0,6.0,8.0,10.0] -- -- divide up a Range into equal-sized sections ----- >>> gridSpace (Range 0 1) 4+-- >>> gridSpace (Range 0.0 1.0) 4 -- [Range 0.0 0.25,Range 0.25 0.5,Range 0.5 0.75,Range 0.75 1.0] data Range a = Range a a-  deriving (Eq, Generic)--instance (Show a) => Show (Range a) where-  show (Range a b) = "Range " <> show a <> " " <> show b+  deriving (Eq, Read, Generic, Data)  instance Eq1 Range where   liftEq f (Range a b) (Range c d) = f a c && f b d -instance Show1 Range where-  liftShowsPrec sp _ d (Range a b) = showsBinaryWith sp sp "Range" d a b+instance (Show a) => Show (Range a) where+  show (Range a b) = "Range " <> show a <> " " <> show b+  showsPrec d p = showParen (d > app_prec) (showString (show p))+    where+      app_prec = 10  instance Functor Range where   fmap f (Range a b) = Range (f a) (f b)@@ -72,7 +75,6 @@   Range fa fb <.> Range a b = Range (fa a) (fb b)  instance Applicative Range where-   pure a = Range a a    (Range fa fb) <*> Range a b = Range (fa a) (fb b)@@ -80,14 +82,9 @@ instance Foldable Range where   foldMap f (Range a b) = f a `mappend` f b -instance Foldable1 Range- instance Traversable Range where   traverse f (Range a b) = Range <$> f a <*> f b -instance Traversable1 Range where-  traverse1 f (Range a b) = Range <$> f a Data.Functor.Apply.<.> f b- instance D.Distributive Range where   collect f x = Range (getL . f <$> x) (getR . f <$> x)     where@@ -95,7 +92,6 @@       getR (Range _ r) = r  instance Representable Range where-   type Rep Range = Bool    tabulate f = Range (f False) (f True)@@ -103,14 +99,13 @@   index (Range l _) False = l   index (Range _ r) True = r -instance (Ord a) => Lattice (Range a) where-+instance (Ord a) => JoinSemiLattice (Range a) where   (\/) = liftR2 min +instance (Ord a) => MeetSemiLattice (Range a) where   (/\) = liftR2 max -instance (Eq a, Ord a) => Space (Range a) where-+instance (Ord a) => Space (Range a) where   type Element (Range a) = a    lower (Range l _) = l@@ -119,11 +114,10 @@    (>.<) = Range -instance (Ord a, Fractional a) => FieldSpace (Range a) where-+instance (Field a, Ord a, FromIntegral a Int) => FieldSpace (Range a) where   type Grid (Range a) = Int -  grid o s n = (+ bool 0 (step / 2) (o == MidPos)) <$> posns+  grid o s n = (+ bool zero (step / two) (o == MidPos)) <$> posns     where       posns = (lower s +) . (step *) . fromIntegral <$> [i0 .. i1]       step = (/) (width s) (fromIntegral n)@@ -139,30 +133,45 @@       ps = grid OuterPos r n  -- | Monoid based on convex hull union-instance (Eq a, Ord a) => Semigroup (Range a) where+instance (Ord a) => Semigroup (Range a) where   (<>) a b = getUnion (Union a <> Union b) --- | Numeric algebra based on Interval arithmetic-instance (Num a, Eq a, Ord a) => Num (Range a) where--  (Range l u) + (Range l' u') = space1 [l + l', u + u']+instance (Additive a, Ord a) => Additive (Range a) where+  (Range l u) + (Range l' u') = unsafeSpace1 [l + l', u + u']+  zero = Range zero zero +instance (Subtractive a, Ord a) => Subtractive (Range a) where   negate (Range l u) = negate u ... negate l -  (Range l u) * (Range l' u') =-    space1 [l * l', l * u', u * l', u * u']+instance (Field a, Ord a) => Multiplicative (Range a) where+  a * b = bool (Range (m - r / (one + one)) (m + r / (one + one))) zero (a == zero || b == zero)+    where+      m = mid a + mid b+      r = width a * width b -  signum (Range l u) = bool (negate 1) 1 (u >= l)+  one = Range (negate one / (one + one)) (one / (one + one)) -  abs (Range l u) = bool (u ... l) (l ... u) (u >= l)+instance (Ord a, Field a) => Divisive (Range a) where+  recip a = bool (Range (-m - one / (two * r)) (-m + one / (two * r))) zero (r == zero)+    where+      m = mid a+      r = width a -  fromInteger x = fromInteger x ... fromInteger x+instance (Field a, Ord a) => Basis (Range a) where+  type Mag (Range a) = Range a+  type Base (Range a) = a+  basis (Range l u) = bool (negate one) one (u >= l)+  magnitude (Range l u) = bool (u ... l) (l ... u) (u >= l) -stepSensible :: (Fractional a, RealFrac a, Floating a, Integral b) => Pos -> a -> b -> a+-- | Find a step that feels pleasent for a 10 digit species.+--+-- >>> stepSensible OuterPos 35 6+-- 5.0+stepSensible :: Pos -> Double -> Int -> Double stepSensible tp span' n =   step + bool 0 (step / 2) (tp == MidPos)   where-    step' = 10.0 ^^ (floor (logBase 10 (span' / fromIntegral n)) :: Integer)+    step' = 10.0 ^^ floor (logBase 10 (span' / fromIntegral n))     err = fromIntegral n / span' * step'     step       | err <= 0.15 = 10.0 * step'@@ -170,30 +179,33 @@       | err <= 0.75 = 2.0 * step'       | otherwise = step' --- | a grid with human sensible (rounded) values+-- | a grid for five-digits per limb species -- -- >>> gridSensible OuterPos False (Range (-12.0) 23.0) 6 -- [-15.0,-10.0,-5.0,0.0,5.0,10.0,15.0,20.0,25.0] gridSensible ::-  (Ord a, RealFrac a, Floating a, Integral b) =>   Pos ->   Bool ->-  Range a ->-  b ->-  [a]+  Range Double ->+  Int ->+  [Double] gridSensible tp inside r@(Range l u) n =-  bool id (filter (`memberOf` r)) inside $-    (+ bool 0 (step / 2) (tp == MidPos)) <$> posns+  bool+    ( bool id (filter (`memberOf` r)) inside $+        (+ bool 0 (step / 2) (tp == MidPos)) <$> posns+    )+    [l - 0.5, l + 0.5]+    (span' == zero)   where     posns = (first' +) . (step *) . fromIntegral <$> [i0 .. i1]     span' = u - l     step = stepSensible tp span' n-    first' = step * fromIntegral (floor (l / step + 1e-6) :: Integer)-    last' = step * fromIntegral (ceiling (u / step - 1e-6) :: Integer)+    first' = step * fromIntegral (floor (l / step + 1e-6))+    last' = step * fromIntegral (ceiling (u / step - 1e-6))     n' = round ((last' - first') / step)     (i0, i1) =       case tp of-        OuterPos -> (0 :: Integer, n')+        OuterPos -> (0, n')         InnerPos -> (1, n' - 1)         LowerPos -> (0, n' - 1)         UpperPos -> (1, n')
src/NumHask/Space/Rect.hs view
@@ -1,69 +1,84 @@-{-# LANGUAGE DeriveGeneric #-}-{-# LANGUAGE DeriveTraversable #-}-{-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE RebindableSyntax #-} {-# LANGUAGE TypeFamilies #-}-{-# OPTIONS_GHC -Wall #-}-{-# OPTIONS_GHC -Wincomplete-patterns #-}+{-# LANGUAGE UndecidableInstances #-} --- | a two-dimensional plane, implemented as a composite of a 'Point' of 'Range's.+-- | A (finite) two-dimensional plane, implemented as a composite of a 'Point' of 'Range's. module NumHask.Space.Rect   ( Rect (..),     pattern Rect,     pattern Ranges,+    rx,+    rz,+    ry,+    rw,+    flipAxes,     corners,     corners4,     projectRect,-    addRect,-    multRect,-    unitRect,     foldRect,+    foldRectUnsafe,     addPoint,-    rotateRect,+    rotationBound,     gridR,     gridF,     aspect,     ratio,+    projectOnR,+    projectOnP,   ) where -import Algebra.Lattice-import Data.Bool (bool)+import Data.Data import Data.Distributive as D import Data.Functor.Compose import Data.Functor.Rep import Data.List.NonEmpty-import Data.Semigroup-import GHC.Exts-import GHC.Generics (Generic)+import NumHask.Prelude hiding (Distributive) import NumHask.Space.Point import NumHask.Space.Range import NumHask.Space.Types-import Prelude  -- $setup+--+-- >>> :m -Prelude+-- >>> :set -XRebindableSyntax+-- >>> import NumHask.Prelude+-- >>> import NumHask.Space --- | a rectangular space often representing a 2-dimensional or XY plane.+-- | a rectangular space often representing a finite 2-dimensional or XY plane. ----- >>> let a = Rect (-1) 1 (-2) 4+-- >>> one :: Rect Double+-- Rect (-0.5) 0.5 (-0.5) 0.5+--+-- >>> zero :: Rect Double+-- Rect 0.0 0.0 0.0 0.0+--+-- >>> one + one :: Rect Double+-- Rect (-1.0) 1.0 (-1.0) 1.0+--+-- >>> let a = Rect (-1.0) 1.0 (-2.0) 4.0 -- >>> a--- Rect -1 1 -2 4+-- Rect (-1.0) 1.0 (-2.0) 4.0+--+-- >>> a * one+-- Rect (-1.0) 1.0 (-2.0) 4.0+-- -- >>> let (Ranges x y) = a -- >>> x--- Range -1 1+-- Range -1.0 1.0 -- >>> y--- Range -2 4+-- Range -2.0 4.0 -- >>> fmap (+1) (Rect 1 2 3 4) -- Rect 2 3 4 5 -- -- as a Space instance with Points as Elements ----- >>> project (Rect 0 1 (-1) 0) (Rect 1 4 10 0) (Point 0.5 1)--- Point 2.5 -10.0--- >>> gridSpace (Rect 0 10 0 1) (Point 2 2)+-- >>> project (Rect 0.0 1.0 (-1.0) 0.0) (Rect 1.0 4.0 10.0 0.0) (Point 0.5 1.0)+-- Point 2.5 (-10.0)+-- >>> gridSpace (Rect 0.0 10.0 0.0 1.0) (Point (2::Int) (2::Int)) -- [Rect 0.0 5.0 0.0 0.5,Rect 0.0 5.0 0.5 1.0,Rect 5.0 10.0 0.0 0.5,Rect 5.0 10.0 0.5 1.0]--- >>> grid MidPos (Rect 0 10 0 1) (Point 2 2)+-- >>> grid MidPos (Rect 0.0 10.0 0.0 1.0) (Point (2::Int) (2::Int)) -- [Point 2.5 0.25,Point 2.5 0.75,Point 7.5 0.25,Point 7.5 0.75] newtype Rect a   = Rect' (Compose Point Range a)@@ -73,23 +88,52 @@       Applicative,       Foldable,       Traversable,-      Generic+      Generic,+      Data     )  -- | pattern of Rect lowerx upperx lowery uppery pattern Rect :: a -> a -> a -> a -> Rect a pattern Rect a b c d = Rect' (Compose (Point (Range a b) (Range c d)))+ {-# COMPLETE Rect #-}  -- | pattern of Ranges xrange yrange pattern Ranges :: Range a -> Range a -> Rect a pattern Ranges a b = Rect' (Compose (Point a b))+ {-# COMPLETE Ranges #-} -instance (Show a) => Show (Rect a) where+instance (Ord a, Additive a, Show a) => Show (Rect a) where   show (Rect a b c d) =-    "Rect " <> show a <> " " <> show b <> " " <> show c <> " " <> show d+    "Rect " <> wrap a <> " " <> wrap b <> " " <> wrap c <> " " <> wrap d+    where+      wrap x = bool (show x) ("(" <> show x <> ")") (x < zero)+  showsPrec d p = showParen (d > app_prec) (showString (show p))+    where+      app_prec = 10 +instance (Read a) => Read (Rect a) where+  readsPrec d s = readParen (d > app_prec) parseRect s+    where+      app_prec = 10+      parseRect r = do+        ("Rect", r1) <- lex r+        (a, r2) <- readsMaybeNeg r1+        (b, r3) <- readsMaybeNeg r2+        (c, r4) <- readsMaybeNeg r3+        (d, r5) <- readsMaybeNeg r4+        pure (Rect a b c d, r5)+      readsMaybeNeg t = case reads t of+        [(v, rest)] -> [(v, rest)]+        [] -> case t of+          '(' : r0 -> case reads r0 of+            [(v, rest)] -> case rest of+              ')' : r -> [(v, r)]+              _ -> []+            _ -> []+          _ -> []+ instance Distributive Rect where   collect f x =     Rect (getA . f <$> x) (getB . f <$> x) (getC . f <$> x) (getD . f <$> x)@@ -100,7 +144,6 @@       getD (Rect _ _ _ d) = d  instance Representable Rect where-   type Rep Rect = (Bool, Bool)    tabulate f =@@ -115,7 +158,6 @@   (<>) = union  instance (Ord a) => Space (Rect a) where-   type Element (Rect a) = Point a    union (Ranges a b) (Ranges c d) = Ranges (a `union` c) (b `union` d)@@ -141,11 +183,10 @@    (|<|) s0 s1 = lower s1 `joinLeq` upper s0 -instance (Ord a, Fractional a, Num a) => FieldSpace (Rect a) where-+instance (FromIntegral a Int, Field a, Ord a) => FieldSpace (Rect a) where   type Grid (Rect a) = Point Int -  grid o s n = (+ bool 0 (step / 2) (o == MidPos)) <$> posns+  grid o s n = (+ bool zero (step / (one + one)) (o == MidPos)) <$> posns     where       posns =         (lower s +) . (step *) . fmap fromIntegral@@ -153,32 +194,67 @@       step = (/) (width s) (fromIntegral <$> n)       (Point x0 y0, Point x1 y1) =         case o of-          OuterPos -> (0, n)-          InnerPos -> (1, n - 1)-          LowerPos -> (0, n - 1)-          UpperPos -> (1, n)-          MidPos -> (0, n - 1)+          OuterPos -> (zero, n)+          InnerPos -> (one, n - one)+          LowerPos -> (zero, n - one)+          UpperPos -> (one, n)+          MidPos -> (zero, n - one)    gridSpace (Ranges rX rY) (Point stepX stepY) =     [ Rect x (x + sx) y (y + sy)-      | x <- grid LowerPos rX stepX,-        y <- grid LowerPos rY stepY+    | x <- grid LowerPos rX stepX,+      y <- grid LowerPos rY stepY     ]     where       sx = width rX / fromIntegral stepX       sy = width rY / fromIntegral stepY +-- | The first X coordinate of a Rect+--+-- >>> rx one+-- -0.5+rx :: Rect a -> a+rx (Rect x _ _ _) = x++-- | The second X coordinate of a Rect+--+-- >>> rz one+-- 0.5+rz :: Rect a -> a+rz (Rect _ z _ _) = z++-- | The first Y coordinate of a Rect+--+-- >>> ry one+-- -0.5+ry :: Rect a -> a+ry (Rect _ _ y _) = y++-- | The second Y coordinate of a Rect+--+-- >>> rw one+-- 0.5+rw :: Rect a -> a+rw (Rect _ _ _ w) = w++-- | flip axes+--+-- >>> flipAxes (Rect 1 2 3 4)+-- Rect 3 4 1 2+flipAxes :: Rect a -> Rect a+flipAxes (Rect x z y w) = Rect y w x z+ -- | create a list of points representing the lower left and upper right corners of a rectangle. ----- >>> corners unitRect--- [Point -0.5 -0.5,Point 0.5 0.5]+-- >>> corners one+-- [Point (-0.5) (-0.5),Point 0.5 0.5] corners :: (Ord a) => Rect a -> [Point a] corners r = [lower r, upper r]  -- | the 4 corners ----- >>> corners4 unitRect--- [Point -0.5 -0.5,Point -0.5 0.5,Point 0.5 -0.5,Point 0.5 0.5]+-- >>> corners4 one+-- [Point (-0.5) (-0.5),Point (-0.5) 0.5,Point 0.5 (-0.5),Point 0.5 0.5] corners4 :: Rect a -> [Point a] corners4 (Rect x z y w) =   [ Point x y,@@ -194,7 +270,7 @@ -- >>> projectRect (Rect 0 1 (-1) 0) (Rect 0 4 0 8) (Rect 0.25 0.75 (-0.75) (-0.25)) -- Rect 1.0 3.0 2.0 6.0 projectRect ::-  (Ord a, Fractional a) =>+  (Field a, Ord a) =>   Rect a ->   Rect a ->   Rect a ->@@ -204,80 +280,67 @@     (Point a' c') = project r0 r1 (Point a c)     (Point b' d') = project r0 r1 (Point b d) --- | Numeric algebra based on interval arithmetioc for addition and unitRect and projection for multiplication-instance (Fractional a, Num a, Eq a, Ord a) => Num (Rect a) where--  (+) = addRect+-- | Numeric algebra based on interval arithmetic for addition and unitRect and projection for multiplication+-- >>> one + one :: Rect Double+-- Rect (-1.0) 1.0 (-1.0) 1.0+instance (Additive a) => Additive (Rect a) where+  (+) (Rect a b c d) (Rect a' b' c' d') =+    Rect (a + a') (b + b') (c + c') (d + d')+  zero = Rect zero zero zero zero +instance (Subtractive a) => Subtractive (Rect a) where   negate = fmap negate -  (*) = multRect--  signum (Rect x z y w) = bool (negate 1) 1 (z >= x && (w >= y))--  abs (Ranges x y) = Ranges (norm x) (norm y)--  fromInteger x = fromInteger x ... fromInteger x+instance (Ord a, Field a) => Multiplicative (Rect a) where+  (*) (Ranges x0 y0) (Ranges x1 y1) =+    Ranges (x0 * x1) (y0 * y1) --- | Rect addition------ >>> unitRect `addRect` unitRect--- Rect -1.0 1.0 -1.0 1.0-addRect :: (Num a) => Rect a -> Rect a -> Rect a-addRect (Rect a b c d) (Rect a' b' c' d') =-  Rect (a + a') (b + b') (c + c') (d + d')+  one = Ranges one one --- | Rect multiplication------ >>> unitRect `multRect` Rect 0 2 0 4--- Rect 0.0 2.0 0.0 4.0-multRect :: (Ord a, Fractional a) => Rect a -> Rect a -> Rect a-multRect (Ranges x0 y0) (Ranges x1 y1) =-  Ranges (x0 `rtimes` x1) (y0 `rtimes` y1)-  where-    rtimes a b = bool (Range (m - r / 2) (m + r / 2)) 0 (a == 0 || b == 0)-      where-        m = mid a + mid b-        r = width a * width b+instance (Ord a, Field a) => Divisive (Rect a) where+  recip (Ranges x y) = Ranges (recip x) (recip y) --- | a unit Rectangle, with values chosen so that width and height are one and mid is zero------ >>> unitRect :: Rect Double--- Rect -0.5 0.5 -0.5 0.5-unitRect :: (Fractional a) => Rect a-unitRect = Ranges rone rone-  where-    rone = Range (-0.5) 0.5+instance (Ord a, Field a) => Basis (Rect a) where+  type Mag (Rect a) = Rect a+  type Base (Rect a) = a+  basis (Rect x z y w) = bool (negate one) one (z >= x && (w >= y))+  magnitude (Ranges x y) = Ranges (magnitude x) (magnitude y)  -- | convex hull union of Rect's ----- >>> foldRect [Rect 0 1 0 1, unitRect]--- Just Rect -0.5 1.0 -0.5 1.0+-- >>> foldRect [Rect 0 1 0 1, one]+-- Just (Rect (-0.5) 1.0 (-0.5) 1.0) foldRect :: (Ord a) => [Rect a] -> Maybe (Rect a) foldRect [] = Nothing foldRect (x : xs) = Just $ sconcat (x :| xs) +-- | convex hull union of Rect's applied to a non-empty structure+--+-- >>> foldRectUnsafe [Rect 0 1 0 1, one]+-- Rect (-0.5) 1.0 (-0.5) 1.0+foldRectUnsafe :: (Foldable f, Ord a) => f (Rect a) -> Rect a+foldRectUnsafe = foldr1 (<>)+ -- | add a Point to a Rect ----- >>> addPoint (Point 0 1) unitRect--- Rect -0.5 0.5 0.5 1.5-addPoint :: (Num a) => Point a -> Rect a -> Rect a+-- >>> addPoint (Point 0 1) one+-- Rect (-0.5) 0.5 0.5 1.5+addPoint :: (Additive a) => Point a -> Rect a -> Rect a addPoint (Point x' y') (Rect x z y w) = Rect (x + x') (z + x') (y + y') (w + y')  -- | rotate the corners of a Rect by x degrees relative to the origin, and fold to a new Rect ----- >>> rotateRect 45 unitRect--- Rect -0.7071067811865475 0.7071067811865475 -5.551115123125783e-17 5.551115123125783e-17-rotateRect :: (Floating a, Ord a) => a -> Rect a -> Rect a-rotateRect d r =-  space1 $ rotate d <$> corners r+-- >>> rotationBound (pi/4) one+-- Rect (-0.7071067811865475) 0.7071067811865475 (-0.7071067811865475) 0.7071067811865475+rotationBound :: (TrigField a, Ord a) => a -> Rect a -> Rect a+rotationBound d = unsafeSpace1 . fmap (rotate d |.) . corners4  -- | Create Rects for a formulae y = f(x) across an x range where the y range is Range 0 y ----- >>> gridR (**2) (Range 0 4) 4+-- >>> gridR (^2) (Range 0 4) 4 -- [Rect 0.0 1.0 0.0 0.25,Rect 1.0 2.0 0.0 2.25,Rect 2.0 3.0 0.0 6.25,Rect 3.0 4.0 0.0 12.25]-gridR :: (Ord a, Fractional a) => (a -> a) -> Range a -> Int -> [Rect a]-gridR f r g = (\x -> Rect (x - tick / 2) (x + tick / 2) 0 (f x)) <$> grid MidPos r g+gridR :: (Field a, FromIntegral a Int, Ord a) => (a -> a) -> Range a -> Int -> [Rect a]+gridR f r g = (\x -> Rect (x - tick / two) (x + tick / two) zero (f x)) <$> grid MidPos r g   where     tick = width r / fromIntegral g @@ -285,19 +348,46 @@ -- -- >>> gridF (\(Point x y) -> x * y) (Rect 0 4 0 4) (Point 2 2) -- [(Rect 0.0 2.0 0.0 2.0,1.0),(Rect 0.0 2.0 2.0 4.0,3.0),(Rect 2.0 4.0 0.0 2.0,3.0),(Rect 2.0 4.0 2.0 4.0,9.0)]-gridF :: (Ord a, Fractional a) => (Point a -> b) -> Rect a -> Grid (Rect a) -> [(Rect a, b)]+gridF :: (Point Double -> b) -> Rect Double -> Grid (Rect Double) -> [(Rect Double, b)] gridF f r g = (\x -> (x, f (mid x))) <$> gridSpace r g  -- | convert a ratio (eg x:1) to a Rect with a height of one. -- -- >>> aspect 2--- Rect -1.0 1.0 -0.5 0.5-aspect :: (Fractional a) => a -> Rect a+-- Rect (-1.0) 1.0 (-0.5) 0.5+aspect :: Double -> Rect Double aspect a = Rect (a * (-0.5)) (a * 0.5) (-0.5) 0.5  -- | convert a Rect to a ratio --+-- >>> :set -XNegativeLiterals -- >>> ratio (Rect (-1) 1 (-0.5) 0.5) -- 2.0-ratio :: (Fractional a) => Rect a -> a+ratio :: (Field a) => Rect a -> a ratio (Rect x z y w) = (z - x) / (w - y)++-- | project a Rect from one Rect to another, preserving relative position, with guards for singleton Rects.+--+-- >>> projectOnR one (Rect 0 1 0 1) (Rect 0 0.5 0 0.5)+-- Rect (-0.5) 0.0 (-0.5) 0.0+projectOnR :: Rect Double -> Rect Double -> Rect Double -> Rect Double+projectOnR new old@(Rect x z y w) ao@(Rect ox oz oy ow)+  | x == z && y == w = ao+  | x == z = Rect ox oz ny nw+  | y == w = Rect nx nz oy ow+  | otherwise = a+  where+    a@(Rect nx nz ny nw) = projectRect old new ao++-- | project a Point from one Rect to another, preserving relative position, with guards for singleton Rects.+--+-- >>> projectOnP one (Rect 0 1 0 1) zero+-- Point (-0.5) (-0.5)+projectOnP :: Rect Double -> Rect Double -> Point Double -> Point Double+projectOnP new old@(Rect x z y w) po@(Point px py)+  | x == z && y == w = po+  | x == z = Point px py'+  | y == w = Point px' py+  | otherwise = Point px' py'+  where+    (Point px' py') = project old new po
src/NumHask/Space/Time.hs view
@@ -1,104 +1,118 @@-{-# LANGUAGE DeriveGeneric #-}-{-# OPTIONS_GHC -Wall #-}-{-# OPTIONS_GHC -Wno-unused-top-binds #-}+{-# LANGUAGE RebindableSyntax #-}+{-# OPTIONS_GHC -fno-warn-name-shadowing #-}  -- | data algorithms related to time (as a Space) module NumHask.Space.Time-  ( parseUTCTime,-    TimeGrain (..),+  ( TimeGrain (..),     floorGrain,     ceilingGrain,+    addGrain,     sensibleTimeGrid,     PosDiscontinuous (..),     placedTimeLabelDiscontinuous,+    placedTimeLabelContinuous,+    fromNominalDiffTime,+    toNominalDiffTime,+    fromDiffTime,+    toDiffTime,   ) where -import qualified Control.Foldl as L-import qualified Data.Text as Text-import Data.Text (Text)-import Data.Time-import GHC.Generics+import Data.Containers.ListUtils (nubOrd)+import Data.Data+import Data.Fixed (Fixed (MkFixed))+import Data.Sequence qualified as Seq+import Data.Text (Text, pack, unpack)+import Data.Time hiding (Hours, Minutes, Seconds)+import NumHask.Prelude+import NumHask.Space.Range import NumHask.Space.Types-import Prelude --- | parse text as per iso8601+-- $setup ----- >>> :set -XOverloadedStrings--- >>> let t0 = parseUTCTime ("2017-12-05" :: Text)--- >>> t0--- Just 2017-12-05 00:00:00 UTC-parseUTCTime :: Text -> Maybe UTCTime-parseUTCTime =-  parseTimeM False defaultTimeLocale (iso8601DateFormat Nothing) . Text.unpack+-- >>> :m -Prelude+-- >>> :set -XRebindableSyntax+-- >>> import NumHask.Prelude+-- >>> import NumHask.Space+-- >>> import NumHask.Space.Time (TimeGrain(..))+-- >>> import Data.Text (Text, pack)+-- >>> import Data.Time hiding (Hours, Minutes, Seconds)  -- | a step in time data TimeGrain-  = Years Integer+  = Years Int   | Months Int   | Days Int   | Hours Int   | Minutes Int   | Seconds Double-  deriving (Show, Eq, Generic)+  deriving (Eq, Show, Generic, Data)  grainSecs :: TimeGrain -> Double-grainSecs (Years n) = fromIntegral n * 365.0 * toDouble nominalDay-grainSecs (Months n) = fromIntegral n * 365.0 / 12 * toDouble nominalDay-grainSecs (Days n) = fromIntegral n * toDouble nominalDay+grainSecs (Years n) = fromIntegral n * 365.0 * fromNominalDiffTime nominalDay+grainSecs (Months n) = fromIntegral n * 365.0 / 12 * fromNominalDiffTime nominalDay+grainSecs (Days n) = fromIntegral n * fromNominalDiffTime nominalDay grainSecs (Hours n) = fromIntegral n * 60 * 60 grainSecs (Minutes n) = fromIntegral n * 60 grainSecs (Seconds n) = n -toDouble :: NominalDiffTime -> Double-toDouble t =-  (/ 1000000000000.0) $-    fromIntegral (floor $ t * 1000000000000 :: Integer)--toDouble' :: DiffTime -> Double-toDouble' =-  (\x -> x / ((10 :: Double) ^ (12 :: Integer))) . fromIntegral . fromEnum+-- | convenience conversion to Double+fromNominalDiffTime :: NominalDiffTime -> Double+fromNominalDiffTime t = fromInteger i * 1e-12+  where+    (MkFixed i) = nominalDiffTimeToSeconds t -fromDouble :: Double -> NominalDiffTime-fromDouble x =+-- | convenience conversion from Double+toNominalDiffTime :: Double -> NominalDiffTime+toNominalDiffTime x =   let d0 = ModifiedJulianDay 0-      days = floor (x / toDouble nominalDay)-      secs = x - fromIntegral days * toDouble nominalDay+      days = floor (x / fromNominalDiffTime nominalDay)+      secs = x - fromIntegral days * fromNominalDiffTime nominalDay       t0 = UTCTime d0 (picosecondsToDiffTime 0)-      t1 = UTCTime (addDays days d0) (picosecondsToDiffTime $ floor (secs / 1.0e-12))+      t1 = UTCTime (addDays (fromIntegral days) d0) (picosecondsToDiffTime . fromIntegral $ floor (secs / 1.0e-12))    in diffUTCTime t1 t0 -fromDouble' :: Double -> DiffTime-fromDouble' d = toEnum . fromEnum $ d * ((10 :: Double) ^ (12 :: Integer))+-- | Convert from 'DiffTime' to seconds (as a Double)+--+-- >>> fromDiffTime $ toDiffTime 1+-- 1.0+fromDiffTime :: DiffTime -> Double+fromDiffTime = (* 1e-12) . fromInteger . diffTimeToPicoseconds +-- | Convert from seconds (as a Double) to 'DiffTime'+-- >>> toDiffTime 1+-- 1s+toDiffTime :: Double -> DiffTime+toDiffTime = picosecondsToDiffTime . fromIntegral . floor . (* 1e12)+ -- | add a TimeGrain to a UTCTime ----- >>> addGrain (Years 1) 5 (UTCTime (fromGregorian 2015 2 28) 0)+-- >>> addGrain (Years 1) 5 (UTCTime (fromGregorian 2015 2 28) (toDiffTime 0)) -- 2020-02-29 00:00:00 UTC ----- >>> addGrain (Months 1) 1 (UTCTime (fromGregorian 2015 2 28) 0)+-- >>> addGrain (Months 1) 1 (UTCTime (fromGregorian 2015 2 28) (toDiffTime 0)) -- 2015-03-31 00:00:00 UTC ----- >>> addGrain (Hours 6) 5 (UTCTime (fromGregorian 2015 2 28) 0)+-- >>> addGrain (Hours 6) 5 (UTCTime (fromGregorian 2015 2 28) (toDiffTime 0)) -- 2015-03-01 06:00:00 UTC ----- >>> addGrain (Seconds 0.001) (60*1000+1) (UTCTime (fromGregorian 2015 2 28) 0)+-- >>> addGrain (Seconds 0.001) (60*1000+1) (UTCTime (fromGregorian 2015 2 28) (toDiffTime 0)) -- 2015-02-28 00:01:00.001 UTC addGrain :: TimeGrain -> Int -> UTCTime -> UTCTime addGrain (Years n) x (UTCTime d t) =-  UTCTime (addDays (-1) $ addGregorianYearsClip (n * fromIntegral x) (addDays 1 d)) t+  UTCTime (addDays (-1) $ addGregorianYearsClip (fromIntegral n * fromIntegral x) (addDays 1 d)) t addGrain (Months n) x (UTCTime d t) =   UTCTime (addDays (-1) $ addGregorianMonthsClip (fromIntegral (n * x)) (addDays 1 d)) t addGrain (Days n) x (UTCTime d t) = UTCTime (addDays (fromIntegral x * fromIntegral n) d) t-addGrain g@(Hours _) x d = addUTCTime (fromDouble (fromIntegral x * grainSecs g)) d-addGrain g@(Minutes _) x d = addUTCTime (fromDouble (fromIntegral x * grainSecs g)) d-addGrain g@(Seconds _) x d = addUTCTime (fromDouble (fromIntegral x * grainSecs g)) d+addGrain g@(Hours _) x d = addUTCTime (toNominalDiffTime (fromIntegral x * grainSecs g)) d+addGrain g@(Minutes _) x d = addUTCTime (toNominalDiffTime (fromIntegral x * grainSecs g)) d+addGrain g@(Seconds _) x d = addUTCTime (toNominalDiffTime (fromIntegral x * grainSecs g)) d  addHalfGrain :: TimeGrain -> UTCTime -> UTCTime addHalfGrain (Years n) (UTCTime d t) =   UTCTime     ( addDays (-1) . (if m' == 1 then addGregorianMonthsClip 6 else id) $-        addGregorianYearsClip d' (addDays 1 d)+        addGregorianYearsClip (fromIntegral d') (addDays 1 d)     )     t   where@@ -113,116 +127,124 @@   where     (d', m') = divMod 2 n addHalfGrain (Days n) (UTCTime d t) =-  (if m' == 1 then addUTCTime (fromDouble (0.5 * grainSecs (Days 1))) else id) $+  (if m' == 1 then addUTCTime (toNominalDiffTime (0.5 * grainSecs (Days 1))) else id) $     UTCTime (addDays (fromIntegral d') d) t   where     (d', m') = divMod 2 n-addHalfGrain g@(Hours _) d = addUTCTime (fromDouble (0.5 * grainSecs g)) d-addHalfGrain g@(Minutes _) d = addUTCTime (fromDouble (0.5 * grainSecs g)) d-addHalfGrain g@(Seconds _) d = addUTCTime (fromDouble (0.5 * grainSecs g)) d+addHalfGrain g@(Hours _) d = addUTCTime (toNominalDiffTime (0.5 * grainSecs g)) d+addHalfGrain g@(Minutes _) d = addUTCTime (toNominalDiffTime (0.5 * grainSecs g)) d+addHalfGrain g@(Seconds _) d = addUTCTime (toNominalDiffTime (0.5 * grainSecs g)) d  -- | compute the floor UTCTime based on the timegrain ----- >>> floorGrain (Years 5) (UTCTime (fromGregorian 1999 1 1) 0)+-- >>> floorGrain (Years 5) (UTCTime (fromGregorian 1999 1 1) (toDiffTime 0)) -- 1995-12-31 00:00:00 UTC ----- >>> floorGrain (Months 3) (UTCTime (fromGregorian 2016 12 30) 0)+-- >>> floorGrain (Months 3) (UTCTime (fromGregorian 2016 12 30) (toDiffTime 0)) -- 2016-09-30 00:00:00 UTC ----- >>> floorGrain (Days 5) (UTCTime (fromGregorian 2016 12 30) 1)+-- >>> floorGrain (Days 5) (UTCTime (fromGregorian 2016 12 30) (toDiffTime 1)) -- 2016-12-30 00:00:00 UTC ----- >>> floorGrain (Minutes 15) (UTCTime (fromGregorian 2016 12 30) (fromDouble' $ 15*60+1))+-- >>> floorGrain (Minutes 15) (UTCTime (fromGregorian 2016 12 30) (toDiffTime $ 15*60+1)) -- 2016-12-30 00:15:00 UTC ----- >>> floorGrain (Seconds 0.1) (UTCTime (fromGregorian 2016 12 30) 0.12)+-- >>> floorGrain (Seconds 0.1) (UTCTime (fromGregorian 2016 12 30) ((toDiffTime 0.12))) -- 2016-12-30 00:00:00.1 UTC floorGrain :: TimeGrain -> UTCTime -> UTCTime-floorGrain (Years n) (UTCTime d _) = UTCTime (addDays (-1) $ fromGregorian y' 1 1) 0+floorGrain (Years n) (UTCTime d _) = UTCTime (addDays (-1) $ fromGregorian y' 1 1) (secondsToDiffTime 0)   where     (y, _, _) = toGregorian (addDays 1 d)-    y' = fromIntegral $ 1 + n * floor (fromIntegral (y - 1) / fromIntegral n :: Double)-floorGrain (Months n) (UTCTime d _) = UTCTime (addDays (-1) $ fromGregorian y m' 1) 0+    y' = fromIntegral $ 1 + n * fromIntegral (floor (fromIntegral (y - 1) / fromIntegral n :: Double))+floorGrain (Months n) (UTCTime d _) = UTCTime (addDays (-1) $ fromGregorian y m' 1) (secondsToDiffTime 0)   where     (y, m, _) = toGregorian (addDays 1 d)-    m' = fromIntegral (1 + fromIntegral n * floor (fromIntegral (m - 1) / fromIntegral n :: Double) :: Integer)-floorGrain (Days _) (UTCTime d _) = UTCTime d 0+    m' = fromIntegral (1 + fromIntegral n * floor (fromIntegral (m - 1) / fromIntegral n :: Double))+floorGrain (Days _) (UTCTime d _) = UTCTime d (secondsToDiffTime 0) floorGrain (Hours h) u@(UTCTime _ t) = addUTCTime x u   where-    s = toDouble' t-    x = fromDouble $ fromIntegral (h * 3600 * fromIntegral (floor (s / (fromIntegral h * 3600)) :: Integer)) - s+    s = fromDiffTime t+    x = toNominalDiffTime $ fromIntegral (h * 3600 * fromIntegral (floor (s / (fromIntegral h * 3600)))) - s floorGrain (Minutes m) u@(UTCTime _ t) = addUTCTime x u   where-    s = toDouble' t-    x = fromDouble $ fromIntegral (m * 60 * fromIntegral (floor (s / (fromIntegral m * 60)) :: Integer)) - s+    s = fromDiffTime t+    x = toNominalDiffTime $ fromIntegral (m * 60 * fromIntegral (floor (s / (fromIntegral m * 60)))) - s floorGrain (Seconds secs) u@(UTCTime _ t) = addUTCTime x u   where-    s = toDouble' t-    x = fromDouble $ (secs * fromIntegral (floor (s / secs) :: Integer)) - s+    s = fromDiffTime t+    x = toNominalDiffTime $ (secs * fromIntegral (floor (s / secs))) - s  -- | compute the ceiling UTCTime based on the timegrain ----- >>> ceilingGrain (Years 5) (UTCTime (fromGregorian 1999 1 1) 0)+-- >>> ceilingGrain (Years 5) (UTCTime (fromGregorian 1999 1 1) (toDiffTime 0)) -- 2000-12-31 00:00:00 UTC ----- >>> ceilingGrain (Months 3) (UTCTime (fromGregorian 2016 12 30) 0)+-- >>> ceilingGrain (Months 3) (UTCTime (fromGregorian 2016 12 30) (toDiffTime 0)) -- 2016-12-31 00:00:00 UTC ----- >>> ceilingGrain (Days 5) (UTCTime (fromGregorian 2016 12 30) 1)+-- >>> ceilingGrain (Days 5) (UTCTime (fromGregorian 2016 12 30) (toDiffTime 1)) -- 2016-12-31 00:00:00 UTC ----- >>> ceilingGrain (Minutes 15) (UTCTime (fromGregorian 2016 12 30) (fromDouble' $ 15*60+1))+-- >>> ceilingGrain (Minutes 15) (UTCTime (fromGregorian 2016 12 30) (toDiffTime $ 15*60+1)) -- 2016-12-30 00:30:00 UTC ----- >>> ceilingGrain (Seconds 0.1) (UTCTime (fromGregorian 2016 12 30) 0.12)+-- >>> ceilingGrain (Seconds 0.1) (UTCTime (fromGregorian 2016 12 30) (toDiffTime 0.12)) -- 2016-12-30 00:00:00.2 UTC ceilingGrain :: TimeGrain -> UTCTime -> UTCTime-ceilingGrain (Years n) (UTCTime d _) = UTCTime (addDays (-1) $ fromGregorian y' 1 1) 0+ceilingGrain (Years n) (UTCTime d _) = UTCTime (addDays (-1) $ fromGregorian y' 1 1) (secondsToDiffTime 0)   where     (y, _, _) = toGregorian (addDays 1 d)-    y' = fromIntegral $ 1 + n * ceiling (fromIntegral (y - 1) / fromIntegral n :: Double)-ceilingGrain (Months n) (UTCTime d _) = UTCTime (addDays (-1) $ fromGregorian y' m'' 1) 0+    y' = fromIntegral $ 1 + n * fromIntegral (ceiling (fromIntegral (y - 1) / fromIntegral n :: Double))+ceilingGrain (Months n) (UTCTime d _) = UTCTime (addDays (-1) $ fromGregorian y' m'' 1) (secondsToDiffTime 0)   where     (y, m, _) = toGregorian (addDays 1 d)     m' = (m + n - 1) `div` n * n     (y', m'') = fromIntegral <$> if m' == 12 then (y + 1, 1) else (y, m' + 1)-ceilingGrain (Days _) (UTCTime d t) = if t == 0 then UTCTime d 0 else UTCTime (addDays 1 d) 0+ceilingGrain (Days _) (UTCTime d t) = if t == secondsToDiffTime 0 then UTCTime d (secondsToDiffTime 0) else UTCTime (addDays 1 d) (secondsToDiffTime 0) ceilingGrain (Hours h) u@(UTCTime _ t) = addUTCTime x u   where-    s = toDouble' t-    x = fromDouble $ fromIntegral (h * 3600 * fromIntegral (ceiling (s / (fromIntegral h * 3600)) :: Integer)) - s+    s = fromDiffTime t+    x = toNominalDiffTime $ fromIntegral (h * 3600 * fromIntegral (ceiling (s / (fromIntegral h * 3600)))) - s ceilingGrain (Minutes m) u@(UTCTime _ t) = addUTCTime x u   where-    s = toDouble' t-    x = fromDouble $ fromIntegral (m * 60 * fromIntegral (ceiling (s / (fromIntegral m * 60)) :: Integer)) - s+    s = fromDiffTime t+    x = toNominalDiffTime $ fromIntegral (m * 60 * fromIntegral (ceiling (s / (fromIntegral m * 60)))) - s ceilingGrain (Seconds secs) u@(UTCTime _ t) = addUTCTime x u   where-    s = toDouble' t-    x = fromDouble $ (secs * fromIntegral (ceiling (s / secs) :: Integer)) - s+    s = fromDiffTime t+    x = toNominalDiffTime $ (secs * fromIntegral (ceiling (s / secs))) - s  -- | whether to include lower and upper times data PosDiscontinuous = PosInnerOnly | PosIncludeBoundaries --- | dates used for time series analysis or attached to charts are often discontinuous, but we want to smooth that reality over and show a continuous range on the axis+-- | Dates used for time series analysis or attached to charts are often discontinuous, but we want to smooth that reality over and show a continuous range on the axis.+-- -- The assumption with getSensibleTimeGrid is that there is a list of discountinuous UTCTimes rather than a continuous range.  Output is a list of index points for the original [UTCTime] and label tuples, and a list of unused list elements. ----- >>> placedTimeLabelDiscontinuous PosIncludeBoundaries (Just "%d %b") 2 [UTCTime (fromGregorian 2017 12 6) 0, UTCTime (fromGregorian 2017 12 29) 0, UTCTime (fromGregorian 2018 1 31) 0, UTCTime (fromGregorian 2018 3 3) 0]+-- >>> placedTimeLabelDiscontinuous PosIncludeBoundaries (Just (pack "%d %b")) 2 [UTCTime (fromGregorian 2017 12 6) (toDiffTime 0), UTCTime (fromGregorian 2017 12 29) (toDiffTime 0), UTCTime (fromGregorian 2018 1 31) (toDiffTime 0), UTCTime (fromGregorian 2018 3 3) (toDiffTime 0)] -- ([(0,"06 Dec"),(1,"31 Dec"),(2,"28 Feb"),(3,"03 Mar")],[]) placedTimeLabelDiscontinuous :: PosDiscontinuous -> Maybe Text -> Int -> [UTCTime] -> ([(Int, Text)], [UTCTime])+placedTimeLabelDiscontinuous _ _ _ [] = ([], []) placedTimeLabelDiscontinuous posd format n ts = (zip (fst <$> inds') labels, rem')   where-    l = minimum ts-    u = maximum ts-    (grain, tps) = sensibleTimeGrid InnerPos n (l, u)+    r@(Range l u) = unsafeSpace1 ts+    (grain, tps) = sensibleTimeGrid InnerPos n r     tps' = case posd of       PosInnerOnly -> tps       PosIncludeBoundaries -> [l] <> tps <> [u]-    (rem', inds) = L.fold (matchTimes tps') ts+    begin = (tps', Seq.empty, zero :: Int)+    done (p, x, _) = (p, toList x)+    step ([], xs, n) _ = ([], xs, n)+    step (p : ps, xs, n) a+      | p == a = step (ps, xs Seq.:|> (n, p), n) a+      | p > a = (p : ps, xs, n + 1)+      | otherwise = step (ps, xs Seq.:|> (n - 1, p), n) a+    (rem', inds) = done $ foldl' step begin ts     inds' = laterTimes inds     fmt = case format of-      Just f -> Text.unpack f+      Just f -> unpack f       Nothing -> autoFormat grain-    labels = Text.pack . formatTime defaultTimeLocale fmt . snd <$> inds'+    labels = pack . formatTime defaultTimeLocale fmt . snd <$> inds'  autoFormat :: TimeGrain -> String autoFormat (Years x)@@ -234,92 +256,90 @@   | x > 3 = "%d/%m/%y %R"   | otherwise = "%R" autoFormat (Minutes _) = "%R"-autoFormat (Seconds _) = "%R%Q"--matchTimes :: [UTCTime] -> L.Fold UTCTime ([UTCTime], [(Int, UTCTime)])-matchTimes ticks = L.Fold step begin (\(p, x, _) -> (p, reverse x))-  where-    begin = (ticks, [], 0)-    step ([], xs, n) _ = ([], xs, n)-    step (p : ps, xs, n) a-      | p == a = step (ps, (n, p) : xs, n) a-      | p > a = (p : ps, xs, n + 1)-      | otherwise = step (ps, (n - 1, p) : xs, n) a+autoFormat (Seconds _) = "%T%Q"  laterTimes :: [(Int, a)] -> [(Int, a)] laterTimes [] = [] laterTimes [x] = [x]-laterTimes (x : xs) = L.fold (L.Fold step (x, []) (\(x0, x1) -> reverse $ x0 : x1)) xs+laterTimes (x : xs) =+  (\(x, xs) -> toList $ xs Seq.:|> x) $+    foldl' step (x, Seq.empty) xs   where-    step ((n, a), rs) (na, aa) = if na == n then ((na, aa), rs) else ((na, aa), (n, a) : rs)+    step ((n, a), rs) (na, aa) =+      bool ((na, aa), rs Seq.:|> (n, a)) ((na, aa), rs) (na == n) +-- | A sensible time grid between two dates, projected onto (0,1) with no attempt to get finnicky.+--+-- >>> placedTimeLabelContinuous PosIncludeBoundaries (Just (pack "%d %b")) 2 (Range (UTCTime (fromGregorian 2017 12 6) (toDiffTime 0)) (UTCTime (fromGregorian 2017 12 29) (toDiffTime 0)))+-- [(0.0,"06 Dec"),(0.4347826086956521,"16 Dec"),(0.8695652173913042,"26 Dec"),(0.9999999999999999,"29 Dec")]+placedTimeLabelContinuous :: PosDiscontinuous -> Maybe Text -> Int -> Range UTCTime -> [(Double, Text)]+placedTimeLabelContinuous posd format n r@(Range l u) = zip tpsd labels+  where+    (grain, tps) = sensibleTimeGrid InnerPos n r+    tps' = case posd of+      PosInnerOnly -> tps+      PosIncludeBoundaries -> nubOrd $ [l] <> tps <> [u]+    fmt = case format of+      Just f -> unpack f+      Nothing -> autoFormat grain+    labels = pack . formatTime defaultTimeLocale fmt <$> tps'+    r' = fromNominalDiffTime $ diffUTCTime u l+    tpsd = (/ r') . fromNominalDiffTime . flip diffUTCTime l <$> tps'+ -- | compute a sensible TimeGrain and list of UTCTimes ----- >>> sensibleTimeGrid InnerPos 2 (UTCTime (fromGregorian 2016 12 31) 0, UTCTime (fromGregorian 2017 12 31) 0)+-- >>> sensibleTimeGrid InnerPos 2 (Range (UTCTime (fromGregorian 2016 12 31) (toDiffTime 0)) (UTCTime (fromGregorian 2017 12 31) (toDiffTime 0))) -- (Months 6,[2016-12-31 00:00:00 UTC,2017-06-30 00:00:00 UTC,2017-12-31 00:00:00 UTC]) ----- >>> sensibleTimeGrid InnerPos 2 (UTCTime (fromGregorian 2017 1 1) 0, UTCTime (fromGregorian 2017 12 30) 0)+-- >>> sensibleTimeGrid InnerPos 2 (Range (UTCTime (fromGregorian 2017 1 1) (toDiffTime 0)) (UTCTime (fromGregorian 2017 12 30) (toDiffTime 0))) -- (Months 6,[2017-06-30 00:00:00 UTC]) ----- >>>  sensibleTimeGrid UpperPos 2 (UTCTime (fromGregorian 2017 1 1) 0, UTCTime (fromGregorian 2017 12 30) 0)+-- >>> sensibleTimeGrid UpperPos 2 (Range (UTCTime (fromGregorian 2017 1 1) (toDiffTime 0)) (UTCTime (fromGregorian 2017 12 30) (toDiffTime 0))) -- (Months 6,[2017-06-30 00:00:00 UTC,2017-12-31 00:00:00 UTC]) ----- >>>sensibleTimeGrid LowerPos 2 (UTCTime (fromGregorian 2017 1 1) 0, UTCTime (fromGregorian 2017 12 30) 0)+-- >>> sensibleTimeGrid LowerPos 2 (Range (UTCTime (fromGregorian 2017 1 1) (toDiffTime 0)) (UTCTime (fromGregorian 2017 12 30) (toDiffTime 0))) -- (Months 6,[2016-12-31 00:00:00 UTC,2017-06-30 00:00:00 UTC])-sensibleTimeGrid :: Pos -> Int -> (UTCTime, UTCTime) -> (TimeGrain, [UTCTime])-sensibleTimeGrid p n (l, u) = (grain, ts)+sensibleTimeGrid :: Pos -> Int -> Range UTCTime -> (TimeGrain, [UTCTime])+sensibleTimeGrid p n (Range l u) = (grain, ts)   where     span' = u `diffUTCTime` l     grain = stepSensibleTime p span' n     first' = floorGrain grain l     last' = ceilingGrain grain u-    n' = round $ toDouble (diffUTCTime last' first') / grainSecs grain :: Integer+    n' =+      round $+        fromNominalDiffTime (diffUTCTime last' first')+          / grainSecs grain     posns = case p of       OuterPos -> take (fromIntegral $ n' + 1)-      InnerPos -> drop (if first' == l then 0 else 1) . take (fromIntegral $ n' + if last' == u then 1 else 0)+      InnerPos ->+        drop (bool one zero (first' == l))+          . take (fromIntegral $ n' + bool zero one (last' == u))       UpperPos -> drop 1 . take (fromIntegral $ n' + 1)       LowerPos -> take (fromIntegral n')       MidPos -> take (fromIntegral n')     ts = case p of-      MidPos -> take (fromIntegral n') $ addHalfGrain grain . (\x -> addGrain grain x first') <$> [0 ..]-      _ -> posns $ (\x -> addGrain grain x first') <$> [0 ..]---- come up with a sensible step for a grid over a Field-stepSensible ::-  (Fractional a, RealFrac a, Floating a) =>-  Pos ->-  a ->-  Int ->-  a-stepSensible tp span' n =-  step-    + if tp == MidPos-      then step / 2-      else 0-  where-    step' = 10 ^^ (floor (logBase 10 (span' / fromIntegral n)) :: Integer)-    err = fromIntegral n / span' * step'-    step-      | err <= 0.15 = 10 * step'-      | err <= 0.35 = 5 * step'-      | err <= 0.75 = 2 * step'-      | otherwise = step'+      MidPos ->+        take (fromIntegral n') $+          addHalfGrain grain+            . (\x -> addGrain grain x first')+            <$> [0 ..]+      _notMid -> posns $ (\x -> addGrain grain x first') <$> [0 ..]  -- come up with a sensible step for a grid over a Field, where sensible means the 18th century -- practice of using multiples of 3 to round stepSensible3 ::-  (Fractional a, Floating a, RealFrac a) =>   Pos ->-  a ->+  Double ->   Int ->-  a+  Double stepSensible3 tp span' n =   step     + if tp == MidPos       then step / 2       else 0   where-    step' = 10 ^^ (floor (logBase 10 (span' / fromIntegral n)) :: Integer)+    step' = 10 ^^ floor (logBase 10 (span' / fromIntegral n))     err = fromIntegral n / span' * step'     step       | err <= 0.05 = 12 * step'@@ -331,18 +351,18 @@ stepSensibleTime :: Pos -> NominalDiffTime -> Int -> TimeGrain stepSensibleTime tp span' n   | yearsstep >= 1 = Years (floor yearsstep)-  | monthsstep >= 1 = Months (fromIntegral (floor monthsstep :: Integer))-  | daysstep >= 1 = Days (fromIntegral (floor daysstep :: Integer))-  | hoursstep >= 1 = Hours (fromIntegral (floor hoursstep :: Integer))-  | minutesstep >= 1 = Minutes (fromIntegral (floor minutesstep :: Integer))+  | monthsstep >= 1 = Months (fromIntegral (floor monthsstep))+  | daysstep >= 1 = Days (fromIntegral (floor daysstep))+  | hoursstep >= 1 = Hours (fromIntegral (floor hoursstep))+  | minutesstep >= 1 = Minutes (fromIntegral (floor minutesstep))   | secondsstep >= 1 = Seconds secondsstep3   | otherwise = Seconds secondsstep   where-    sp = toDouble span'+    sp = fromNominalDiffTime span'     minutes = sp / 60     hours = sp / (60 * 60)-    days = sp / toDouble nominalDay-    years = sp / 365 / toDouble nominalDay+    days = sp / fromNominalDiffTime nominalDay+    years = sp / 365 / fromNominalDiffTime nominalDay     months' = years * 12     yearsstep = stepSensible tp years n     monthsstep = stepSensible3 tp months' n
src/NumHask/Space/Types.hs view
@@ -1,5 +1,6 @@ {-# LANGUAGE DefaultSignatures #-}-{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE RebindableSyntax #-} {-# LANGUAGE TypeFamilies #-} {-# OPTIONS_HADDOCK hide #-} @@ -10,9 +11,11 @@     Intersection (..),     FieldSpace (..),     mid,+    interpolate,     project,     Pos (..),     space1,+    unsafeSpace1,     memberOf,     contains,     disjoint,@@ -24,13 +27,25 @@     widenEps,     scale,     move,+    Transform (..),+    inverseTransform,+    Affinity (..),+    (|.),+    rotate,   ) where -import Protolude-import Data.Foldable+import Control.Monad+import NumHask.Prelude+import Prelude qualified as P --- | Space is a continuous range of numbers that contains elements and has an upper and lower value.+-- $setup+-- >>> :m -Prelude+-- >>> :set -XRebindableSyntax+-- >>> import NumHask.Prelude+-- >>> import NumHask.Space++-- | A 'Space' is a continuous set of numbers. Continuous here means that the set has an upper and lower bound, and an element that is between these two bounds is a member of the 'Space'. -- -- > a `union` b == b `union` a -- > a `intersection` b == b `intersection` a@@ -39,9 +54,7 @@ -- > norm (norm a) = norm a -- > a |>| b == b |<| a -- > a |.| singleton a- class Space s where-   -- | the underlying element in the space   type Element s :: Type @@ -57,7 +70,6 @@    -- | the intersection of two spaces   intersection :: s -> s -> s-   default intersection :: (Ord (Element s)) => s -> s -> s   intersection a b = l >.< u     where@@ -66,7 +78,6 @@    -- | the union of two spaces   union :: s -> s -> s-   default union :: (Ord (Element s)) => s -> s -> s   union a b = l >.< u     where@@ -74,16 +85,16 @@       u = upper a `max` upper b    -- | Normalise a space so that+  --   -- > lower a \/ upper a == lower a   -- > lower a /\ upper a == upper a-  norm :: s -> s-  norm s = lower s ... upper s+  normalise :: s -> s+  normalise s = lower s ... upper s    -- | create a normalised space from two elements   infix 3 ...    (...) :: Element s -> Element s -> s-   default (...) :: (Ord (Element s)) => Element s -> Element s -> s   (...) a b = (a `min` b) >.< (a `max` b) @@ -96,7 +107,6 @@   infixl 7 |.|    (|.|) :: Element s -> s -> Bool-   default (|.|) :: (Ord (Element s)) => Element s -> s -> Bool   (|.|) a s = (a >= lower s) && (upper s >= a) @@ -104,7 +114,6 @@   infixl 7 |>|    (|>|) :: s -> s -> Bool-   default (|>|) :: (Ord (Element s)) => s -> s -> Bool   (|>|) s0 s1 =     lower s0 >= upper s1@@ -113,7 +122,6 @@   infixl 7 |<|    (|<|) :: s -> s -> Bool-   default (|<|) :: (Ord (Element s)) => s -> s -> Bool   (|<|) s0 s1 =     lower s1 <= upper s0@@ -136,13 +144,13 @@ memberOf = (|.|)  -- | distance between boundaries-width :: (Space s, Num (Element s)) => s -> Element s+width :: (Space s, Subtractive (Element s)) => s -> Element s width s = upper s - lower s  -- | create a space centered on a plus or minus b infixl 6 +/- -(+/-) :: (Space s, Num (Element s)) => Element s -> Element s -> s+(+/-) :: (Space s, Subtractive (Element s)) => Element s -> Element s -> s a +/- b = a - b ... a + b  -- | a convex hull@@ -159,10 +167,10 @@  -- | a space that can be divided neatly ----- > space1 (grid OuterPos s g) == s+-- > unsafeSpace1 (grid OuterPos s g) == s -- > getUnion (sconcat (Union <$> (gridSpace s g))) == s-class (Space s, Num (Element s)) => FieldSpace s where-+class (Space s, Field (Element s)) => FieldSpace s where+  -- | the type that divides or quotients the space   type Grid s :: Type    -- | create equally-spaced elements across a space@@ -172,44 +180,60 @@   gridSpace :: s -> Grid s -> [s]  -- | Pos suggests where points should be placed in forming a grid across a field space.-data Pos =-  -- | include boundaries-  OuterPos |-  -- | don't include boundaries-  InnerPos |-  -- | include the lower boundary-  LowerPos |-  -- | include the upper boundary-  UpperPos |-  -- | use the mid-point of the space-  MidPos deriving (Show, Eq)+data Pos+  = -- | include boundaries+    OuterPos+  | -- | don't include boundaries+    InnerPos+  | -- | include the lower boundary+    LowerPos+  | -- | include the upper boundary+    UpperPos+  | -- | use the mid-point of the space+    MidPos+  deriving (Show, Eq)  -- | middle element of the space-mid :: (Space s, Fractional (Element s)) => s -> Element s-mid s = (lower s + upper s) / 2.0+mid :: (Space s, Field (Element s)) => s -> Element s+mid s = (lower s + upper s) / (one + one) --- | project a data point from one space to another, preserving relative position+-- | interpolate a space --+-- > interpolate s x == project s (zero ... one) x+interpolate :: (Space s, Ring (Element s)) => s -> Element s -> Element s+interpolate s x = lower s + x * width s++-- | project an element from one space to another, preserving relative position.+-- -- > project o n (lower o) = lower n -- > project o n (upper o) = upper n+-- > project o n (mid o) = mid n -- > project a a x = x-project :: (Space s, Fractional (Element s)) => s -> s -> Element s -> Element s+project :: (Space s, Field (Element s)) => s -> s -> Element s -> Element s project s0 s1 p =   ((p - lower s0) / (upper s0 - lower s0)) * (upper s1 - lower s1) + lower s1 --- | the containing space of a non-empty Traversable--- > all $ space1 a `contains` <$> a-space1 :: (Space s, Traversable f) => f (Element s) -> s-space1 = foldr1 union . fmap singleton+-- | the containing space of a non-empty Traversable.+--+-- partial function.+--+-- > all $ unsafeSpace1 a `contains` <$> a+unsafeSpace1 :: (Space s, Traversable f) => f (Element s) -> s+unsafeSpace1 = P.foldr1 union . fmap singleton +-- | Maybe containing space of a traversable.+space1 :: (Space s, Traversable f) => f (Element s) -> Maybe s+space1 s = bool (Just $ unsafeSpace1 s) Nothing (null s)+ -- | lift a monotone function (increasing or decreasing) over a given space monotone :: (Space a, Space b) => (Element a -> Element b) -> a -> b-monotone f s = space1 [f (lower s), f (upper s)]+monotone f s = unsafeSpace1 [f (lower s), f (upper s)]  -- | a small space eps ::   ( Space s,-    Fractional (Element s)+    FromRational (Element s),+    Field (Element s)   ) =>   Element s ->   Element s ->@@ -219,7 +243,7 @@ -- | widen a space widen ::   ( Space s,-    Num (Element s)+    Ring (Element s)   ) =>   Element s ->   s ->@@ -229,7 +253,8 @@ -- | widen by a small amount widenEps ::   ( Space s,-    Fractional (Element s)+    FromRational (Element s),+    Ring (Element s)   ) =>   Element s ->   s ->@@ -237,9 +262,81 @@ widenEps accuracy = widen (accuracy * 1e-6)  -- | Scale a Space. (scalar multiplication)-scale :: (Num (Element s), Space s) => Element s -> s -> s+scale :: (Multiplicative (Element s), Space s) => Element s -> s -> s scale e s = (e * lower s) ... (e * upper s)  -- | Move a Space. (scalar addition)-move :: (Num (Element s), Space s) => Element s -> s -> s+move :: (Additive (Element s), Space s) => Element s -> s -> s move e s = (e + lower s) ... (e + upper s)++-- | linear transform + translate of a point-like number+--+-- > (x, y) -> (ax + by + c, dx + ey + d)+--+-- or+--+-- \[+-- \begin{pmatrix}+-- a & b & c\\+-- d & e & f\\+-- 0 & 0 & 1+-- \end{pmatrix}+-- \begin{pmatrix}+-- x\\+-- y\\+-- 1+-- \end{pmatrix}+-- \]+data Transform a = Transform+  { ta :: !a,+    tb :: !a,+    tc :: !a,+    td :: !a,+    te :: !a,+    tf :: !a+  }+  deriving (Eq, Show, Functor, Foldable, Traversable)++-- | Calculate the inverse of a transformation.+inverseTransform :: (Eq a, Field a) => Transform a -> Maybe (Transform a)+inverseTransform (Transform a b c d e f) =+  let det = a * e - b * d+   in bool+        ( Just+            ( Transform+                (a / det)+                (d / det)+                (-(a * c + d * f) / det)+                (b / det)+                (e / det)+                (-(b * c + e * f) / det)+            )+        )+        Nothing+        (det == zero)++-- | An 'Affinity' is something that can be subjected to an affine transformation in 2-dimensional space, where affine means a linear matrix operation or a translation (+).+--+-- https://en.wikipedia.org/wiki/Affine_transformation+class Affinity a b | a -> b where+  transform :: Transform b -> a -> a++infix 3 |.++-- | Apply a 'Transform' to an 'Affinity'+(|.) :: (Affinity a b) => Transform b -> a -> a+(|.) = transform++instance (Multiplicative a, Additive a) => Affinity (Transform a) a where+  transform (Transform a' b' c' d' e' f') (Transform a b c d e f) =+    Transform+      (a * a' + b' * d)+      (a' * b + b' * e)+      (a' * c + b' * f + c')+      (d' * a + e' * d)+      (d' * b + e' * e)+      (d' * c + e' * f + f')++-- | Rotate an 'Affinity' (counter-clockwise)+rotate :: (TrigField a) => a -> Transform a+rotate a = Transform (cos a) (-sin a) zero (sin a) (cos a) zero
− test/test.hs 
@@ -1,16 +0,0 @@-{-# OPTIONS_GHC -Wall #-}--module Main where--import Test.DocTest-import Prelude--main :: IO ()-main =-  doctest-    [ "src/NumHask/Space/Histogram.hs",-      "src/NumHask/Space/Point.hs",-      "src/NumHask/Space/Range.hs",-      "src/NumHask/Space/Rect.hs",-      "src/NumHask/Space/Time.hs"-    ]