LTS-17.9. Add Cyclic (from astro), doctests

bjornbm · May 1, 2021 · 6dbba52 · 6dbba52
1 parent c366183
commit 6dbba52
Show file tree

Hide file tree

Showing 12 changed files with 234 additions and 14 deletions.
diff --git a/dimensional-experimental.cabal b/dimensional-experimental.cabal
@@ -19,25 +19,45 @@ Description:
 
 Category:            Math, Physics
 Build-Type:          Simple
-Build-Depends:
+Extra-source-files:  README, LICENSE
+cabal-version:       >=1.10
+
+library
+  default-language:    Haskell2010
+  hs-source-dirs:      src
+  Build-Depends:
 
     base,
     time,
+    ieee754,
     dimensional,
     dimensional-dk-experimental,
     dimensional-homo,
     numtype-dk,
     numeric-quest,
     ad
 
-Exposed-Modules:
+  Exposed-Modules:
 
     Numeric.Units.Dimensional.AD
     Numeric.Units.Dimensional.Constants
+    Numeric.Units.Dimensional.Cyclic
     Numeric.Units.Dimensional.Formulae
     Numeric.Units.Dimensional.LinearAlgebra.VectorAD
     Numeric.Units.Dimensional.LinearAlgebra.Rotation
     Numeric.Units.Dimensional.LinearAlgebra.PosVel
 
-Extra-source-files:  README, LICENSE
 
+test-suite doctests
+  type:           exitcode-stdio-1.0
+  ghc-options:    -threaded
+  default-language:    Haskell2010
+  hs-source-dirs: test
+  main-is:        Doctests.hs
+  build-depends:  base
+                , doctest >= 0.8
+                , Glob
+                , ieee754
+                , QuickCheck
+                , dimensional
+                , dimensional-homo
diff --git a/Numeric/Units/Dimensional/AD.hs → src/Numeric/Units/Dimensional/AD.hs b/Numeric/Units/Dimensional/AD.hs → src/Numeric/Units/Dimensional/AD.hs
diff --git a/Numeric/Units/Dimensional/Constants.hs → src/Numeric/Units/Dimensional/Constants.hs b/Numeric/Units/Dimensional/Constants.hs → src/Numeric/Units/Dimensional/Constants.hs
diff --git a/src/Numeric/Units/Dimensional/Cyclic.hs b/src/Numeric/Units/Dimensional/Cyclic.hs
@@ -0,0 +1,152 @@
+-- | This module contains various function for working with functions
+-- that are is cyclic in meaning (but not in value). An obvious example
+-- are angles, where in many cases there is no value in distinguishing
+-- between x and x + n2π.
+
+module Numeric.Units.Dimensional.Cyclic where
+
+import qualified Prelude
+import Numeric.Units.Dimensional.Prelude
+
+import Data.AEq
+
+
+-- Adjusting values with respect to a reference
+-- --------------------------------------------
+
+-- | Returns the integral part of a value, leaving only a fractional
+-- part in the interval (-1,1).
+-- (@/= snd . properFraction@).
+fractionalPart :: RealFrac a => Dimensionless a -> Dimensionless a
+fractionalPart = fmap (snd . properFraction)
+
+-- | Removes the integral part of a value so that a fractional part
+-- in the interval [0,1) is returned. This differs in behavior from
+-- @(snd . properFraction)@ for negative values:
+--   @(snd . properFraction)   1.2@  =  0.2
+--   @adjustZeroOne            1.2@  =  0.2
+--   @(snd . properFraction) (-1.2)@ = -0.2
+--   @adjustZeroOne          (-1.2)@ =  0.8
+adjustZeroOne :: RealFrac a => Dimensionless a -> Dimensionless a
+--adjustZeroOne = fmap (\x -> x Prelude.- fromIntegral (floor x))
+adjustZeroOne x = x - fmap (fromIntegral . floor) x
+
+-- | @adjustAbove1 min x@ adjusts x by an integer so that it lies in the
+-- interval [min,min+1).
+adjust1Above :: RealFrac a
+             => Dimensionless a -> Dimensionless a -> Dimensionless a
+adjust1Above min x = adjustZeroOne (x - min) + min
+
+-- | @adjustAbove cycle min x@ adjusts x by an integer number or cycles
+-- so that it lies in the interval [min,min+cycle).
+adjustAbove :: RealFrac a
+            => Quantity d a -> Quantity d a -> Quantity d a -> Quantity d a
+adjustAbove cycle min x = adjust1Above (min/cycle) (x/cycle) * cycle
+
+-- | @adjustAbout1 center x@ adjusts x by an integer so that it lies in the
+-- interval [center-1/2,center+1/2).
+adjust1About :: RealFrac a
+             => Dimensionless a -> Dimensionless a -> Dimensionless a
+adjust1About center = adjust1Above (center - _1/_2)
+
+-- | @adjustAbout center cycle x@ adjusts x by an integer number or cycles
+-- so that it lies in the interval [min-cycle/2,min+cycle/2).
+adjustAbout :: RealFrac a
+            => Quantity d a -> Quantity d a -> Quantity d a -> Quantity d a
+adjustAbout center cycle x = adjust1About (center/cycle) (x/cycle) * cycle
+
+
+-- Adjusting angles
+-- ----------------
+
+-- | @adjustAngle center a@ adjusts the angle @a@ to
+-- be within ±π of @center@.
+adjustAngle :: RealFloat a => Angle a -> Angle a -> Angle a
+adjustAngle center = adjustAbout center tau
+
+-- | Adjusts an angle to the range [-pi,pi).
+plusMinusPi :: RealFloat a => Angle a -> Angle a
+plusMinusPi = adjustAngle _0
+
+-- | Adjusts an angle to the range [0,2pi).
+zeroTwoPi :: RealFloat a => Angle a -> Angle a
+zeroTwoPi = adjustAngle pi
+
+-- | Adjusts an angle to the range [0,tau).
+zeroTau :: RealFloat a => Angle a -> Angle a
+zeroTau = zeroTwoPi
+
+
+-- Consistency metric for cyclic values
+-- ------------------------------------
+
+-- | Assume that y(x) is cyclic in meaning (but not in value, as
+-- for e.g. angles) with a periodicity @period@ and cycle length @cycle@.
+-- Then @cyclesOff (x,y) period cycle@ compute the approximate (closest)
+-- integral number of cycles that @y@ differs from what one would expect
+-- given @x@ and an assumption that y = cycle * x / period.
+cyclesOff :: RealFrac a
+          => (Quantity d a, Quantity dy a)
+          -> Quantity d a -> Quantity dy a -> Dimensionless a
+cyclesOff (x,y) period cycle = cyclesOff1 (x/period, y/cycle)
+
+-- | Same as 'cyclesOff' but assumes that both the period and the cycle
+-- of y(x) is one.
+cyclesOff1 :: RealFrac a => (Dimensionless a, Dimensionless a) -> Dimensionless a
+cyclesOff1 (x, y) = fmap (fromIntegral . round) (x - y)
+
+
+-- Adjusting for relative total rotation
+-- -------------------------------------
+
+-- | Assume that y(x) is cyclic in meaning (but not in value, as
+-- for e.g. angles). Then @adjustCyclic (x0,y0) (x1,y1) period cycle@
+-- returns a new @y1@ adjusted so that the difference @y1 - y0@
+-- corresponds roughly to the difference @x1 - x0@.
+-- (See also adjustCyclic1.)
+adjustCyclic :: RealFrac a
+             => (Quantity d a, Quantity dy a) -> (Quantity d a, Quantity dy a)
+             -> Quantity d a -> Quantity dy a -> Quantity dy a
+adjustCyclic (x0,y0) (x1,y1) period cycle =
+  y1 + cyclesOff (x1-x0, y1-y0) period cycle * cycle
+  -- Could be defined as:
+  --   adjustCyclic1 (x0/period,y0/cycle) (x1/period,y1/cycle) * cycle
+  -- but that has worse numerical properties!
+
+-- | Assume that y(x) is cyclic in meaning (but not in value, as
+-- for e.g. angles) where the meaning has a cycle (in y) of 1 with
+-- a period (in x) of 1. Then @adjustCyclic1 (x0,y0) (x1,y1)@
+-- returns a new @y1@ adjusted so that @y1 - y0@ is roughly the same
+-- as @x1 - x0@.
+--
+-- Property of returned y1:
+--   | (y1 - y0) - (x1 - x0) | < 0.5
+--
+-- (This is a "normalized" version of 'adjustCyclic' for cycle and
+-- period of 1.)
+adjustCyclic1 :: RealFrac a
+              => (Dimensionless a, Dimensionless a)
+              -> (Dimensionless a, Dimensionless a)
+              -> Dimensionless a
+adjustCyclic1 (x0,y0) (x1,y1) = y1 + cyclesOff1 (x1-x0, y1-y0)
+
+-- | Adjust assuming x0 = 0 and y0 = 0
+adjustCyclic0 :: RealFrac a
+              => (Quantity d a, Quantity dy a)
+              -> Quantity d a -> Quantity dy a -> Quantity dy a
+adjustCyclic0 (x1,y1) period cycle = y1 + cyclesOff (x1,y1) period cycle * cycle
+
+
+-- Angle comparisons
+-- -----------------
+
+-- | Compares two angles for cyclic equality.
+(==~) :: (RealFloat a, Eq a) => Angle a -> Angle a -> Bool
+x ==~ y = plusMinusPi x == plusMinusPi y
+
+-- | Compares two angles for approximate cyclic equality.
+(~==~) :: (RealFloat a, AEq a) => Angle a -> Angle a -> Bool
+x ~==~ y = plusMinusPi x ~== plusMinusPi y
+        ||   zeroTwoPi x ~==   zeroTwoPi y  -- move the boundaries.
+
+infixl 4 ==~, ~==~
diff --git a/Numeric/Units/Dimensional/Formulae.hs → src/Numeric/Units/Dimensional/Formulae.hs b/Numeric/Units/Dimensional/Formulae.hs → src/Numeric/Units/Dimensional/Formulae.hs
diff --git a/...Units/Dimensional/LinearAlgebra/PosVel.hs → ...Units/Dimensional/LinearAlgebra/PosVel.hs b/...Units/Dimensional/LinearAlgebra/PosVel.hs → ...Units/Dimensional/LinearAlgebra/PosVel.hs
@@ -8,15 +8,21 @@
 {-# LANGUAGE StandaloneDeriving #-}
 {-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE DeriveAnyClass #-}
 
 module Numeric.Units.Dimensional.LinearAlgebra.PosVel where
 
 import qualified Prelude
+import Data.AEq
 import Numeric.NumType.DK.Integers (TypeInt (Pos2))
+import Numeric.Units.Dimensional.Cyclic
 import Numeric.Units.Dimensional.Prelude
 import Numeric.Units.Dimensional.LinearAlgebra
 import Numeric.Units.Dimensional.LinearAlgebra.VectorAD (applyLinear)
 
+-- $setup
+-- >>> :set -XDataKinds
+-- >>> import Data.AEq
 
 -- Type synonyms for clearer documentation.
 
@@ -40,9 +46,17 @@ data Sph d a = Sph { magnitude      :: Quantity d a
                    , rightAscension :: PlaneAngle a
                    } deriving (Eq)
 deriving instance (KnownDimension d, Show a, Real a) => Show (Sph d a)  -- Needed since d unknown.
+
+instance (RealFloat a, AEq a) => AEq (Sph d a) where
+  s1 ~== s2 =
+    magnitude s1 ~== magnitude s2 &&
+    zenith s1 ~==~ zenith s2 &&
+    rightAscension s1 ~==~ rightAscension s2
+
 type SPos = Sph DRadius
+
 data SVel a = SVel (Velocity a) (AngularVelocity a) (AngularVelocity a)
-  deriving (Eq, Show)
+  deriving (Eq, Show, AEq)
 
 
 
@@ -75,8 +89,13 @@ hourAngle = rightAscension
 
 -- Converting
 -- ----------
--- Converts a cartesian position vector into a spherical position vector.
 
+-- | Converts a cartesian vector into a spherical vector.
+--
+-- >>> let v = _1 <: _2 <:. _3 in s2c (c2s v) ~== v
+-- True
+-- >>> let v = fromListErr ([1, 2, 3] *~~ meter) in s2c (c2s v) ~== v
+-- True
 c2s :: RealFloat a => Car d a -> Sph d a
 c2s c = Sph r zen az
   where
@@ -86,8 +105,14 @@ c2s c = Sph r zen az
     az  = atan2 y x
 
 
--- Converts a spherical position vector into a cartesian position vector.
+-- | Converts a spherical vector into a cartesian vector.
-
+--
+-- >>> let v = Sph _4 _2 _3 in c2s (s2c v) ~== v
+-- True
+-- >>> let v = Sph (4 *~ meter) _2 _3 in c2s (s2c v) ~== v
+-- True
+--
+-- prop> let v = Sph _0 _3 _2 in s2c v == (nullVector :: CVel Double)
 s2c :: Floating a => Sph d a -> Car d a
 s2c s = x <: y <:. z
   where
@@ -96,8 +121,7 @@ s2c s = x <: y <:. z
     y = r * sin zen * sin az
     z = r * cos zen
 
--- Convert declination/latitude to zenith.
+-- | Convert declination/latitude to zenith.
-
 toZenith :: Floating a => Angle a -> Zenith a
 toZenith decl = tau / _4 - decl
 
@@ -123,6 +147,9 @@ type SPosVel a = (SPos a, SVel a)
 -- Conversions
 -- -----------
 
+-- |
+-- >>> let pv = (fromListErr ([1, 2, 3] *~~ meter), fromListErr ([2, 2, 3] *~~ (meter / second))) in s2cEphem (c2sEphem pv) ~== pv
+-- True
 c2sEphem :: RealFloat a => CPosVel a -> SPosVel a
 c2sEphem = pvf . applyLinear (fp . c2s) -- unlinearize (c2s . linearize c :: RealFloat b => Time b -> SPos b)
   where

diff --git a/...its/Dimensional/LinearAlgebra/Rotation.hs → ...its/Dimensional/LinearAlgebra/Rotation.hs b/...its/Dimensional/LinearAlgebra/Rotation.hs → ...its/Dimensional/LinearAlgebra/Rotation.hs
@@ -7,6 +7,15 @@ import Numeric.Units.Dimensional.Prelude
 import Numeric.Units.Dimensional.LinearAlgebra
 import qualified Prelude
 
+-- $setup
+-- >>> :set -XTypeSynonymInstances
+-- >>> :set -XFlexibleInstances
+-- >>> :set -XFlexibleContexts
+-- >>> :set -XDataKinds
+-- >>> import Data.AEq
+-- >>> import Test.QuickCheck
+-- >>> instance (HasDimension (Data.Proxy.Proxy d), Arbitrary a, Num a) => Arbitrary (Quantity d a) where arbitrary = (*~ siUnit) <$> arbitrary
+
 -- | The type of homogenous vectors with three elements.
 type Homo3 d = Vec d 3
 
@@ -26,16 +35,22 @@ type Homo33 d = Mat d 3 3
 -- Rotation matrices. Rotates a vector by the given angle (analogous
 -- to rotating the coordinate system in opposite direction).
 
+-- |
+-- prop> rotX a ~== rotV unit_x (a :: Angle Double)
 rotX :: Floating a => PlaneAngle a -> Homo33 DOne a
 rotX a =   (_1 <:    _0 <:.            _0 )
        |:  (_0 <: cos a <:. negate (sin a))
        |:. (_0 <: sin a <:.         cos a )
 
+-- |
+-- prop> rotY a ~== rotV unit_y (a :: Angle Double)
 rotY :: Floating a => PlaneAngle a -> Homo33 DOne a
 rotY a =   (        cos a  <: _0 <:. sin a)
        |:  (           _0  <: _1 <:.    _0)
        |:. (negate (sin a) <: _0 <:. cos a)
 
+-- |
+-- prop> rotZ a ~== rotV unit_z (a :: Angle Double)
 rotZ :: Floating a => PlaneAngle a -> Homo33 DOne a
 rotZ a =   (cos a <: negate (sin a) <:. _0)
        |:  (sin a <:         cos a  <:. _0)

diff --git a/...its/Dimensional/LinearAlgebra/VectorAD.hs → ...its/Dimensional/LinearAlgebra/VectorAD.hs b/...its/Dimensional/LinearAlgebra/VectorAD.hs → ...its/Dimensional/LinearAlgebra/VectorAD.hs
diff --git a/stack.yaml b/stack.yaml
@@ -15,7 +15,7 @@
 # resolver:
 #  name: custom-snapshot
 #  location: "./custom-snapshot.yaml"
-resolver: lts-13.19
+resolver: lts-17.9
 allow-newer: true
 
 # User packages to be built.
@@ -38,15 +38,13 @@ allow-newer: true
 # will not be run. This is useful for tweaking upstream packages.
 packages:
 - '.'
-- '../dimensional'
-- '../dimensional-dk-experimental'
 - '../dimensional-homo'
 # Dependency packages to be pulled from upstream that are not in the resolver
 # (e.g., acme-missiles-0.3)
 extra-deps:
 - numeric-quest-0.2.0.2
-- HaTeX-3.20.0.0
+- git: https://github.com/bjornbm/dimensional-dk-experimental.git
-- wl-pprint-extras-3.5.0.5
+  commit: 0516f3f10f2589b414a8ede6b949e44dbc1ceebc  # origin/dimensional-head
 
 # Override default flag values for local packages and extra-deps
 flags: {}

diff --git a/test/Doctests.hs b/test/Doctests.hs
@@ -0,0 +1,8 @@
+import System.FilePath.Glob (glob)
+import Test.DocTest (doctest)
+
+main = glob "src/**/*.*hs" >>= doctest
+-- main = do
+--     src_hs <- glob "src/**/*.hs"
+--     src_lhs <- glob "src/**/*.lhs"
+--     doctest (src_hs ++ src_lhs ++ ["test/TestInstances.hs", "test/TestUtil.hs"])
diff --git a/tests/PosVelTest.hs → test/PosVelTest.hs b/tests/PosVelTest.hs → test/PosVelTest.hs
diff --git a/tests/VectorADTest.hs → test/VectorADTest.hs b/tests/VectorADTest.hs → test/VectorADTest.hs