diff --git a/lib/ChebFun.hs b/lib/ChebFun.hs index 49b799d..287ab89 100644 --- a/lib/ChebFun.hs +++ b/lib/ChebFun.hs @@ -1,11 +1,13 @@ {-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE UndecidableSuperClasses #-} module ChebFun ( ChebVal() , ChebFun(..) - , zeroCF, constCF, idCF - , ChebProd(..), ChebUnit(..) + , compactify, uncompactify, clamp + , law_ChebFun_compact_inv + , law_ChebFun_compact_inv' ) where import Prelude hiding ( id, (.), curry, uncurry @@ -13,211 +15,78 @@ import Prelude hiding ( id, (.), curry, uncurry , Foldable(..) , Applicative(..), (<$>) ) --- import qualified Prelude as P +import Control.Exception (assert) import Data.Constraint -import qualified Data.Vector.Unboxed as U -import Data.Vector.Unboxed.Deriving import GHC.Generics -import GHC.TypeLits import qualified Test.QuickCheck as QC import Applicative import Category import Chebyshev --- import Foldable +import Foldable +import qualified Function as F import Functor --- import Unfoldable import Unbox +import Unfoldable -class (RealFloat a, Unbox a) => ChebVal a -instance (RealFloat a, Unbox a) => ChebVal a +-- | ChebVal +class (RealFloat a, Obj (Dom v) a) => ChebVal v a +instance (RealFloat a, Obj (Dom v) a) => ChebVal v a -newtype ChebFun (n :: Nat) a b = ChebFun (NUVector n b) +instance F.MCategory (ChebVal v) + + + +-- | ChebFun +newtype ChebFun v a b = ChebFun { getChebFun :: v b } deriving ( Eq, Ord, Show, Generic - , QC.Arbitrary, QC.CoArbitrary -- , QC.Function - -- , P.Functor, P.Foldable - -- , Functor - , Applicative + , QC.Arbitrary, QC.CoArbitrary ) -instance Functor (ChebFun n a) where - type Dom (ChebFun n a) = Dom (NUVector n) - type Cod (ChebFun n a) = Cod (NUVector n) +instance (Functor v, Cod v ~ (->)) => Functor (ChebFun v a) where + type Dom (ChebFun v a) = Dom v + type Cod (ChebFun v a) = Cod v proveCod = Sub Dict fmap f (ChebFun xs) = ChebFun (fmap f xs) --- instance Foldable (ChebFun n a) where --- foldMap f (ChebFun xs) = foldMap f xs --- foldr f z (ChebFun xs) = foldr f z xs --- foldl f z (ChebFun xs) = foldl f z xs --- toList (ChebFun xs) = toList xs - - +instance (Applicative v, Cod v ~ (->)) => Applicative (ChebFun v a) where + pure = ChebFun . pure + ChebFun fs <*> ChebFun xs = ChebFun (fs <*> xs) + liftA2 f (ChebFun xs) (ChebFun ys) = ChebFun (liftA2 f xs ys) + liftA2' f (ChebFun xs, ChebFun ys) = ChebFun (liftA2' f (xs, ys)) -zeroCF :: forall n a b. (KnownNat n, Unbox b, Num b) => ChebFun n a b -zeroCF = ChebFun (WithSize (U.replicate (intVal @n) 0)) +instance Foldable v => F.Morphism (ChebFun v) where + type Cat (ChebFun v) = ChebVal v + chase (ChebFun cs) = chebyshev cs . realToFrac -constCF :: forall n a b. (KnownNat n, Unbox b, Num b, 1 <= n) - => b -> ChebFun n a b -constCF x = ChebFun (WithSize (U.fromListN (intVal @n) - ([x] ++ replicate (intVal @n - 1) 0))) +instance (Foldable v, Sized v, Unfoldable v) + => F.Discretization (ChebFun v) where + discretize = ChebFun . chebyshevApprox (size @v) -idCF :: forall n a. (KnownNat n, Unbox a, Num a, 2 <= n) => ChebFun n a a -idCF = ChebFun (WithSize (U.fromListN (intVal @n) - ([0, 1] ++ replicate (intVal @n - 2) 0))) +clamp :: Ord a => a -> a -> a -> a +clamp lo hi = max lo . min hi -instance (KnownNat n, 2 <= n) => Category (ChebFun n) where - type Obj (ChebFun n) = ChebVal - id = idCF - g . f = unapply (apply g . apply f) -- TODO: do this better! - -- g . f = ChebFun (chebyshevApprox' n2 n (apply g . apply f)) - -- where n = intVal @n - -- n2 = n ^ (2 :: Int) - apply (ChebFun cs) = chebyshev cs . realToFrac - unapply = ChebFun . chebyshevApprox (intVal @n) +-- | input: \( [-\infty, \infty] \), +-- output: \( [-1; 1] \) +compactify :: RealFloat a => a -> a +compactify x = let y = 2 / pi * atan x + in assert (-1 <= y && y <= 1) y +uncompactify :: RealFloat a => a -> a +uncompactify y = + assert (-1 <= y && y <= 1) $ + tan (pi / 2 * clamp (-1) 1 y) +-- | prop> uncompactify . compactify = id +law_ChebFun_compact_inv :: RealFloat a => FnEqual a a +law_ChebFun_compact_inv = (uncompactify . compactify) `fnEqual` id @(->) -newtype ChebProd (n :: Nat) a b = ChebProd { getChebProd :: (a, b) } - deriving ( Eq, Ord, Read, Show, Generic - , QC.Arbitrary, QC.CoArbitrary - -- , Functor, Foldable, Unfoldable, Applicative - ) - -derivingUnbox "ChebProd" - [t| forall n a b. (Unbox a, Unbox b) => ChebProd n a b -> (a, b) |] - [| getChebProd |] - [| ChebProd |] - -bimap :: (a -> c) -> (b -> d) -> ChebProd n a b -> ChebProd n c d -bimap f g (ChebProd (x, y)) = ChebProd (f x, g y) - -bipure :: a -> b -> ChebProd n a b -bipure x y = ChebProd (x, y) - -biliftA2 :: (a -> c -> e) -> (b -> d -> f) -> - ChebProd n a b -> ChebProd n c d -> ChebProd n e f -biliftA2 f g (ChebProd (x1, y1)) (ChebProd (x2, y2)) = - ChebProd (f x1 x2, g y1 y2) - -bifoldMap :: (a -> c) -> (b -> d) -> (c -> d -> e) -> ChebProd n a b -> e -bifoldMap f g op (ChebProd (x, y)) = f x `op` g y - -instance (Num a, Num b) => Num (ChebProd n a b) where - (+) = biliftA2 (+) (+) - (*) = biliftA2 (*) (*) - negate = bimap negate negate - abs = bimap abs abs - signum = bimap signum signum - fromInteger i = bipure (fromInteger i) (fromInteger i) - -instance (Fractional a, Fractional b) => Fractional (ChebProd n a b) where - recip = bimap recip recip - fromRational r = bipure (fromRational r) (fromRational r) - -instance (Real a, Real b) => Real (ChebProd n a b) where - toRational (ChebProd (x, _)) = toRational x - -instance (RealFrac a, RealFrac b) => RealFrac (ChebProd n a b) where - properFraction (ChebProd (x, _)) = undefined -- properFraction x - -instance (Floating a, Floating b) => Floating (ChebProd n a b) where - pi = bipure pi pi - exp = bimap exp exp - log = bimap log log - sin = bimap sin sin - cos = bimap cos cos - asin = bimap asin asin - acos = bimap acos acos - atan = bimap atan atan - sinh = bimap sinh sinh - cosh = bimap cosh cosh - asinh = bimap asinh asinh - acosh = bimap acosh acosh - atanh = bimap atanh atanh - -instance (RealFloat a, RealFloat b) => RealFloat (ChebProd n a b) where - floatRadix _ = 0 - floatDigits _ = 0 - floatRange _ = (0, 0) - decodeFloat _ = (0, 0) - encodeFloat i j = bipure (encodeFloat i j) (encodeFloat i j) - isNaN = bifoldMap isNaN isNaN (||) - isInfinite = bifoldMap isInfinite isInfinite (||) - isDenormalized = bifoldMap isDenormalized isDenormalized (||) - isNegativeZero = bifoldMap isNegativeZero isNegativeZero (||) - isIEEE _ = False - - - -newtype ChebUnit n = ChebUnit () - deriving ( Eq, Ord, Read, Show, Generic - , QC.Arbitrary, QC.CoArbitrary - ) - -derivingUnbox "ChebUnit" - [t| forall n. ChebUnit n -> () |] - [| const () |] - [| ChebUnit |] - -instance Num (ChebUnit n) where - _ + _ = ChebUnit () - _ * _ = ChebUnit () - negate _ = ChebUnit () - abs _ = ChebUnit () - signum _ = ChebUnit () - fromInteger _ = ChebUnit () - -instance Fractional (ChebUnit n) where - recip _ = ChebUnit () - fromRational _ = ChebUnit () - -instance Real (ChebUnit n) where - toRational _ = 0 - -instance RealFrac (ChebUnit n) where - properFraction _ = (0, ChebUnit ()) - -instance Floating (ChebUnit n) where - pi = ChebUnit () - exp _ = ChebUnit () - log _ = ChebUnit () - sin _ = ChebUnit () - cos _ = ChebUnit () - asin _ = ChebUnit () - acos _ = ChebUnit () - atan _ = ChebUnit () - sinh _ = ChebUnit () - cosh _ = ChebUnit () - asinh _ = ChebUnit () - acosh _ = ChebUnit () - atanh _ = ChebUnit () - -instance RealFloat (ChebUnit n) where - floatRadix _ = 0 - floatDigits _ = 0 - floatRange _ = (0, 0) - decodeFloat _ = (0, 0) - encodeFloat _ _ = ChebUnit () - isNaN _ = False - isInfinite _ = False - isDenormalized _ = False - isNegativeZero _ = False - isIEEE _ = False - - - -instance CatProd (ChebProd n) where - type Unit (ChebProd n) = ChebUnit n - punit = ChebUnit () - prod = ChebProd - unprod = getChebProd - -instance (KnownNat n, 2 <= n) => Cartesian (ChebFun n) where - type Prod (ChebFun n) = ChebProd n - proveCartesian = Sub Dict +-- | prop> compactify . uncompactify = id +law_ChebFun_compact_inv' :: RealFloat a => FnEqual a a +law_ChebFun_compact_inv' = + (compactify . uncompactify) `fnEqual` id @(->) diff --git a/lib/ChebFun1.hs b/lib/ChebFun1.hs new file mode 100644 index 0000000..7540133 --- /dev/null +++ b/lib/ChebFun1.hs @@ -0,0 +1,264 @@ +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE UndecidableInstances #-} + +module ChebFun1 + ( ChebVal() + , ChebFun(..) + -- , zeroCF, constCF, idCF + , compactify, uncompactify, clamp + , law_ChebFun_compact_inv + , law_ChebFun_compact_inv' + , ChebProd(..), ChebUnit(..) + ) where + +import Prelude hiding ( id, (.), curry, uncurry + , Functor(..) + , Foldable(..) + , Applicative(..), (<$>) + ) +-- import qualified Prelude as P + +import Control.Exception (assert) +import Data.Constraint +-- import qualified Data.Vector.Unboxed as U +import Data.Vector.Unboxed.Deriving +import GHC.Generics +import GHC.TypeLits +import qualified Test.QuickCheck as QC + +import Applicative +import Category +import Chebyshev +-- import Foldable +import Functor +-- import Unfoldable +import Unbox + + + +class (RealFloat a, Unbox a) => ChebVal a +instance (RealFloat a, Unbox a) => ChebVal a + +-- TODO: Make NUVector a parameter; use "Dom v in ChebVal definition" +-- TODO: Support rational numbers somehow +newtype ChebFun (n :: Nat) a b = ChebFun (NUVector n b) + deriving ( Eq, Ord, Show, Generic + , QC.Arbitrary, QC.CoArbitrary -- , QC.Function + -- , P.Functor, P.Foldable + -- , Functor + , Applicative + ) + +instance Functor (ChebFun n a) where + type Dom (ChebFun n a) = Dom (NUVector n) + type Cod (ChebFun n a) = Cod (NUVector n) + proveCod = Sub Dict + fmap f (ChebFun xs) = ChebFun (fmap f xs) + +-- instance Foldable (ChebFun n a) where +-- foldMap f (ChebFun xs) = foldMap f xs +-- foldr f z (ChebFun xs) = foldr f z xs +-- foldl f z (ChebFun xs) = foldl f z xs +-- toList (ChebFun xs) = toList xs + + + +-- zeroCF :: forall n a b. (KnownNat n, Unbox b, Num b) => ChebFun n a b +-- zeroCF = ChebFun (WithSize (U.replicate (intVal @n) 0)) +-- +-- constCF :: forall n a b. (KnownNat n, Unbox b, Num b, 1 <= n) +-- => b -> ChebFun n a b +-- constCF x = ChebFun (WithSize (U.fromListN (intVal @n) +-- ([x] ++ replicate (intVal @n - 1) 0))) +-- +-- idCF :: forall n a. (KnownNat n, Unbox a, Num a, 2 <= n) => ChebFun n a a +-- idCF = ChebFun (WithSize (U.fromListN (intVal @n) +-- ([0, 1] ++ replicate (intVal @n - 2) 0))) + + + +clamp :: Ord a => a -> a -> a -> a +clamp lo hi = max lo . min hi + +-- | input: \( [-\infty, \infty] \), +-- output: \( [-1; 1] \) +-- +-- > compactify x +-- +-- >>> uncompactivy . compactify +-- +-- prop> uncompactivy . compactify = id +compactify :: RealFloat a => a -> a +compactify x = let y = 2 / pi * atan x + in assert (-1 <= y && y <= 1) y + +uncompactify :: RealFloat a => a -> a +uncompactify y = + assert (-1 <= y && y <= 1) $ + -- assert (-1.1 <= y && y <= 1.1) $ + tan (pi / 2 * clamp (-1) 1 y) + +law_ChebFun_compact_inv :: RealFloat a => FnEqual a a +law_ChebFun_compact_inv = (uncompactify . compactify) `fnEqual` id @(->) + +law_ChebFun_compact_inv' :: RealFloat a => FnEqual a a +law_ChebFun_compact_inv' = + (compactify . uncompactify) `fnEqual` id @(->) + + + +instance KnownNat n => Category (ChebFun n) where + type Obj (ChebFun n) = ChebVal + -- TODO: do these better! + id = unapply id + g . f = unapply (apply g . apply f) + -- | [-1; 1] -> [-\infty; \infty] + apply (ChebFun cs) = chebyshev cs . realToFrac + unapply = ChebFun . chebyshevApprox (intVal @n) + -- | [-\infty; \infty] -> [-\infty; \infty] + -- apply (ChebFun cs) = chebyshev cs . realToFrac . compactify + -- unapply f = ChebFun . chebyshevApprox (intVal @n) $ f . uncompactify + -- | [-1; 1] -> [-1; 1] + -- apply (ChebFun cs) = compactify . chebyshev cs . realToFrac + -- unapply f = ChebFun . chebyshevApprox (intVal @n) $ uncompactify . f + + + +newtype ChebProd (n :: Nat) a b = ChebProd { getChebProd :: (a, b) } + deriving ( Eq, Ord, Read, Show, Generic + , QC.Arbitrary, QC.CoArbitrary + -- , Functor, Foldable, Unfoldable, Applicative + ) + +derivingUnbox "ChebProd" + [t| forall n a b. (Unbox a, Unbox b) => ChebProd n a b -> (a, b) |] + [| getChebProd |] + [| ChebProd |] + +bimap :: (a -> c) -> (b -> d) -> ChebProd n a b -> ChebProd n c d +bimap f g (ChebProd (x, y)) = ChebProd (f x, g y) + +bipure :: a -> b -> ChebProd n a b +bipure x y = ChebProd (x, y) + +biliftA2 :: (a -> c -> e) -> (b -> d -> f) -> + ChebProd n a b -> ChebProd n c d -> ChebProd n e f +biliftA2 f g (ChebProd (x1, y1)) (ChebProd (x2, y2)) = + ChebProd (f x1 x2, g y1 y2) + +bifoldMap :: (a -> c) -> (b -> d) -> (c -> d -> e) -> ChebProd n a b -> e +bifoldMap f g op (ChebProd (x, y)) = f x `op` g y + +instance (Num a, Num b) => Num (ChebProd n a b) where + (+) = biliftA2 (+) (+) + (*) = biliftA2 (*) (*) + negate = bimap negate negate + abs = bimap abs abs + signum = bimap signum signum + fromInteger i = bipure (fromInteger i) (fromInteger i) + +instance (Fractional a, Fractional b) => Fractional (ChebProd n a b) where + recip = bimap recip recip + fromRational r = bipure (fromRational r) (fromRational r) + +instance (Real a, Real b) => Real (ChebProd n a b) where + toRational (ChebProd (x, _)) = toRational x + +instance (RealFrac a, RealFrac b) => RealFrac (ChebProd n a b) where + properFraction (ChebProd (x, _)) = undefined -- properFraction x + +instance (Floating a, Floating b) => Floating (ChebProd n a b) where + pi = bipure pi pi + exp = bimap exp exp + log = bimap log log + sin = bimap sin sin + cos = bimap cos cos + asin = bimap asin asin + acos = bimap acos acos + atan = bimap atan atan + sinh = bimap sinh sinh + cosh = bimap cosh cosh + asinh = bimap asinh asinh + acosh = bimap acosh acosh + atanh = bimap atanh atanh + +instance (RealFloat a, RealFloat b) => RealFloat (ChebProd n a b) where + floatRadix _ = 0 + floatDigits _ = 0 + floatRange _ = (0, 0) + decodeFloat _ = (0, 0) + encodeFloat i j = bipure (encodeFloat i j) (encodeFloat i j) + isNaN = bifoldMap isNaN isNaN (||) + isInfinite = bifoldMap isInfinite isInfinite (||) + isDenormalized = bifoldMap isDenormalized isDenormalized (||) + isNegativeZero = bifoldMap isNegativeZero isNegativeZero (||) + isIEEE _ = False + + + +newtype ChebUnit n = ChebUnit () + deriving ( Eq, Ord, Read, Show, Generic + , QC.Arbitrary, QC.CoArbitrary + ) + +derivingUnbox "ChebUnit" + [t| forall n. ChebUnit n -> () |] + [| const () |] + [| ChebUnit |] + +instance Num (ChebUnit n) where + _ + _ = ChebUnit () + _ * _ = ChebUnit () + negate _ = ChebUnit () + abs _ = ChebUnit () + signum _ = ChebUnit () + fromInteger _ = ChebUnit () + +instance Fractional (ChebUnit n) where + recip _ = ChebUnit () + fromRational _ = ChebUnit () + +instance Real (ChebUnit n) where + toRational _ = 0 + +instance RealFrac (ChebUnit n) where + properFraction _ = (0, ChebUnit ()) + +instance Floating (ChebUnit n) where + pi = ChebUnit () + exp _ = ChebUnit () + log _ = ChebUnit () + sin _ = ChebUnit () + cos _ = ChebUnit () + asin _ = ChebUnit () + acos _ = ChebUnit () + atan _ = ChebUnit () + sinh _ = ChebUnit () + cosh _ = ChebUnit () + asinh _ = ChebUnit () + acosh _ = ChebUnit () + atanh _ = ChebUnit () + +instance RealFloat (ChebUnit n) where + floatRadix _ = 0 + floatDigits _ = 0 + floatRange _ = (0, 0) + decodeFloat _ = (0, 0) + encodeFloat _ _ = ChebUnit () + isNaN _ = False + isInfinite _ = False + isDenormalized _ = False + isNegativeZero _ = False + isIEEE _ = False + + + +instance CatProd (ChebProd n) where + type Unit (ChebProd n) = ChebUnit n + punit = ChebUnit () + prod = ChebProd + unprod = getChebProd + +instance (KnownNat n, 2 <= n) => Cartesian (ChebFun n) where + type Prod (ChebFun n) = ChebProd n + proveCartesian = Sub Dict diff --git a/lib/Chebyshev.hs b/lib/Chebyshev.hs index c2e97fa..8e521fc 100644 --- a/lib/Chebyshev.hs +++ b/lib/Chebyshev.hs @@ -1,5 +1,8 @@ +{-# LANGUAGE UndecidableInstances #-} + -- | Adapted from [math-functions-0.2.1.0] Numeric.Polynomial.Chebyshev module Chebyshev ( chebyshev + -- , Floating'(..) , chebyshevApprox , chebyshevApprox' ) where @@ -9,7 +12,9 @@ import Prelude hiding ( id, (.), curry, uncurry , Foldable(..) ) +-- import Data.Fixed import Data.Maybe +-- import Data.Ratio import Category import Foldable @@ -55,23 +60,43 @@ chebyshev cs' x = -- | Approximate a given function via Chebyshev polynomials. -- See . chebyshevApprox :: ( Unfoldable v, Obj (Dom v) b - , RealFloat a, Fractional b + , RealFrac a, Floating' a, Fractional b ) => Int -> (a -> b) -> v b chebyshevApprox n = chebyshevApprox' (2 * n) n chebyshevApprox' :: forall f a b. ( Unfoldable f, Obj (Dom f) b - , RealFloat a, Fractional b + , RealFrac a, Floating' a, Fractional b ) => Int -> Int -> (a -> b) -> f b chebyshevApprox' np nc f = (fromJust . fromList) [coeff i | i <- [0..nc-1]] where coeff j = rf ((if j == 0 then 1 else 2) / fi np) * sum [f (x k) * rf (c j k) | k <- [0..np-1]] t :: Int -> a - t k = pi * (fi k + 0.5) / fi np + t k = fi (2 * k + 1) / fi (2 * np) x :: Int -> a - x k = cos (t k) + x k = cospi' (t k) c :: Int -> Int -> a - c j k = cos (fi j * t k) + -- c j k = cos (pi * fi j * t k) + c j k = cospi' (u j k) + u :: Int -> Int -> a + u j k = fi jk / fi (2 * np) + where jk = (j * (2 * k + 1)) `mod` (4 * np) fi = fromIntegral rf = realToFrac + +class Floating' a where cospi' :: a -> a +instance Floating a => Floating' a where cospi' x = cos (pi * x) + +-- class Floating' a where +-- cospi' :: a -> a +-- +-- instance Floating a => Floating' a where +-- cospi' x = cos (pi * x) +-- +-- instance Integral a => Floating' (Ratio a) where +-- cospi' x = +-- let x1 = x `mod'` 2 +-- x2 = if x1 <= 1 then x1 else 2 - x1 +-- in if x2 <= 1/2 then c x2 else - c (1 - x2) +-- where c y = (1 - 4 * y^2) / (1 + y^2) diff --git a/lib/Domains.hs b/lib/Domains.hs new file mode 100644 index 0000000..66132f6 --- /dev/null +++ b/lib/Domains.hs @@ -0,0 +1,56 @@ +{-# LANGUAGE UndecidableInstances #-} + +module Domains where + +import Data.Validity +import Numeric.IEEE +import GHC.Generics + + + +-- x: (-\infty; +\infty) +-- sqrt: [0; +\infty) +-- sin: [-1; +1] +-- log: (0; +\infty) +-- atan: (-\pi/2; +\pi/2) +-- acos: [0; \pi) + +data Constant + = NegInf + | Neg2Pi + | NegPi + | NegOne + | NegPi2 + | Zero + | PosInf + | Pos2Pi + | PosPi + | PosOne + | PosPi2 + +class Value (v :: Constant) where value :: IEEE a => a +instance Value 'NegInf where value = -infinity +instance Value 'Neg2Pi where value = -2*pi +instance Value 'NegPi where value = -pi +instance Value 'NegOne where value = -1 +instance Value 'NegPi2 where value = -pi/2 +instance Value 'Zero where value = 0 +instance Value 'PosPi2 where value = pi/2 +instance Value 'PosOne where value = 1 +instance Value 'PosPi where value = pi +instance Value 'Pos2Pi where value = 2*pi +instance Value 'PosInf where value = infinity + + + +newtype BoundedReal (min :: Constant) (max :: Constant) a = B { getB :: a } + deriving ( Eq, Ord, Read, Show, Generic + , Floating, Fractional, IEEE, Num, Real, RealFloat, RealFrac) + +instance (Value min, Value max, IEEE a) => Bounded (BoundedReal min max a) where + minBound = B (value @min) + maxBound = B (value @max) + +instance (Bounded (BoundedReal min max a), IEEE a, Validity a) + => Validity (BoundedReal min max a) where + isValid x = (isNaN x || x >= minBound && x <= maxBound) && isValid (getB x) diff --git a/lib/Foldable.hs b/lib/Foldable.hs index d863151..d1a7c1b 100644 --- a/lib/Foldable.hs +++ b/lib/Foldable.hs @@ -24,7 +24,6 @@ import Data.Functor.Sum as F import Data.List.NonEmpty hiding (length, toList) import Data.Monoid hiding (First(..)) import Data.Proxy -import Data.Semigroup import GHC.Base (build) import Category diff --git a/lib/Function.hs b/lib/Function.hs index 56c4b2f..74f8584 100644 --- a/lib/Function.hs +++ b/lib/Function.hs @@ -16,11 +16,6 @@ import Data.Constraint import Data.Kind import qualified Data.Vector.Unboxed as U () import GHC.Generics -import GHC.TypeLits - -import ChebFun -import Chebyshev -import Unbox @@ -124,15 +119,3 @@ instance Morphism (-.>) where instance Discretization (-.>) where discretize = NFun - - - --- | ChebVal -instance MCategory ChebVal - -instance Morphism (ChebFun n) where - type Cat (ChebFun n) = ChebVal - chase (ChebFun cs) = chebyshev cs . realToFrac - -instance KnownNat n => Discretization (ChebFun n) where - discretize = ChebFun . chebyshevApprox (intVal @n) diff --git a/lib/Sized.hs b/lib/Sized.hs new file mode 100644 index 0000000..ce51ce4 --- /dev/null +++ b/lib/Sized.hs @@ -0,0 +1,19 @@ +module Sized + ( intVal + , Sized(..) + ) where + +import Data.Kind +import Data.Proxy +import GHC.TypeLits + + + +-- | Use as 'intVal @n' +intVal :: forall n. KnownNat n => Int +intVal = fromInteger (natVal (Proxy @n)) + +class KnownNat (Size v) => Sized (v :: Type -> Type) where + type Size v :: Nat + size :: Int + size = intVal @(Size v) diff --git a/lib/Unbox.hs b/lib/Unbox.hs index c887ca6..221e108 100644 --- a/lib/Unbox.hs +++ b/lib/Unbox.hs @@ -27,10 +27,8 @@ import Prelude hiding ( id, (.), curry, uncurry import qualified Prelude import Data.Constraint -import Data.Kind import Data.Maybe import Data.Monoid -import Data.Proxy import qualified Data.Vector as V import Data.Vector.Unboxed (Unbox) import qualified Data.Vector.Unboxed as U @@ -45,6 +43,7 @@ import Category import Comonad import Foldable import Functor +import Sized import Unfoldable @@ -179,15 +178,6 @@ instance Indexed U.Vector where --- | Use as 'intVal @n' -intVal :: forall n. KnownNat n => Int -intVal = fromInteger (natVal (Proxy @n)) - -class KnownNat (Size v) => Sized (v :: Type -> Type) where - type Size v :: Nat - size :: Int - size = intVal @(Size v) - newtype WithSize (n :: Nat) v a = WithSize (v a) deriving Generic diff --git a/package.yaml b/package.yaml index 27ddc9e..f04fa01 100644 --- a/package.yaml +++ b/package.yaml @@ -77,8 +77,9 @@ library: # - bifunctors - constraints # - data-default + - ieee754 - quickcheck-instances - # - validity + - validity # - validity-vector - vector # # - vector-space diff --git a/test/ChebFun1Spec.hs b/test/ChebFun1Spec.hs new file mode 100644 index 0000000..7f2a595 --- /dev/null +++ b/test/ChebFun1Spec.hs @@ -0,0 +1,131 @@ +{-# OPTIONS_GHC -Wno-incomplete-patterns #-} + +module ChebFun1Spec where + +import Prelude hiding ( id, (.), curry, uncurry + , Functor(..) + , Foldable(..) + , concat, concatMap, sum, product, maximum, minimum + , and, or, all, any + , Applicative(..), (<$>) + ) + +import Data.Fixed +import Test.QuickCheck hiding (scale) + +import Applicative +import Category +import ChebFun1 +import Foldable +import Functor +import Unbox + + + +maxabs :: (Foldable f, k ~ Dom f, Obj k a, Num a, Ord a) => f a -> a +maxabs = foldl (\r -> max r . abs) 0 + +sumabs :: (Foldable f, k ~ Dom f, Obj k a, Num a) => f a -> a +sumabs = foldl (\r -> (r +) . abs) 0 + +-- | Approximate floating-point equality +approx :: (RealFrac a, Show a) => a -> a -> a -> Property +approx delta x y = counterexample + (show x ++ " /= " ++ show y ++ "; " ++ + "absdiff = " ++ show (y - x) ++ ", " ++ + "reldiff = " ++ show ((y - x) / maxabs [x, y]) ++ "; " ++ + "delta = " ++ show delta) + (abs (x - y) <= delta) + +limitArg :: Real a => a -> a +limitArg x = mod' (x + 1) 2 - 1 +limitFun :: (Unbox b, Fractional b, Ord b) => ChebFun n a b -> ChebFun n a b +limitFun f = let ChebFun cs = f + scale = max 1 (sumabs cs) + in ChebFun (fmap (unapply (/ scale)) cs) + + + +type N = 10 +type N' = 2 + +type CA = Double +type CB = Double +type CC = Double + +eps :: Double +eps = 1.0e-12 + + + +prop_ChebFun1_compact_inv :: CA -> Property +prop_ChebFun1_compact_inv x' = + uncurry (approx eps) (getFnEqual law_ChebFun_compact_inv x) + where x = compactify x' + +prop_ChebFun1_compact_inv' :: CA -> Property +prop_ChebFun1_compact_inv' x' = + uncurry (approx eps) (getFnEqual law_ChebFun_compact_inv x) + where x = mod' x' 2 - 1 + + + +prop_ChebFun1_Category_comp_id_left :: ChebFun N CA CB -> CA -> Property +prop_ChebFun1_Category_comp_id_left f' x' = + uncurry (approx eps) (getFnEqual (law_Category_comp_id_left f) x) + where x = limitArg x' + f = limitFun f' + +prop_ChebFun1_Category_comp_id_right :: ChebFun N CA CB -> CA -> Property +prop_ChebFun1_Category_comp_id_right f' x' = + uncurry (approx eps) (getFnEqual (law_Category_comp_id_right f) x) + where x = limitArg x' + f = limitFun f' + -- where x = mod' x' 2 - 1 + -- eps = let (ChebFun cs) = f in eps1 * max 1 (maxabs cs) + +-- TODO: Make this work for larger N' +prop_ChebFun1_Category_comp_assoc :: + ChebFun N' CA CC -> ChebFun N' CB CA -> ChebFun N' CA CB -> CA -> Property +prop_ChebFun1_Category_comp_assoc h' g' f' x' = + uncurry (approx eps) (getFnEqual (law_Category_comp_assoc h g f) x) + where x = limitArg x' + f = limitFun f' + g = limitFun g' + h = limitFun h' + + + +prop_ChebFun1_Functor_id :: ChebFun N CB CA -> Property +prop_ChebFun1_Functor_id xs = + uncurry (===) (getFnEqual law_Functor_id xs) + +prop_ChebFun1_Functor_assoc :: + Fun CB CC -> Fun CA CB -> ChebFun N CB CA -> Property +prop_ChebFun1_Functor_assoc (Fn g) (Fn f) xs = + uncurry (===) (getFnEqual (law_Functor_assoc (UFun g) (UFun f)) xs) + + + +prop_ChebFun1_Applicative_id' :: Fun CA CB -> ChebFun N CB CA -> Property +prop_ChebFun1_Applicative_id' (Fn f) xs = + uncurry (===) (getEqual (law_Applicative_id' (UFun f) xs)) + +prop_ChebFun1_Applicative_id_left' :: + Fun (CA *#* CB) CC -> CA -> ChebFun N CB CB -> Property +prop_ChebFun1_Applicative_id_left' (Fn f) x ys = + uncurry (===) (getEqual (law_Applicative_id_left' (UFun f) x ys)) + +prop_ChebFun1_Applicative_id_right' :: + Fun (CA *#* CB) CC -> ChebFun N CB CA -> CB -> Property +prop_ChebFun1_Applicative_id_right' (Fn f) xs y = + uncurry (===) (getEqual (law_Applicative_id_right' (UFun f) xs y)) + +prop_ChebFun1_Applicative_assoc' :: + Fun CA CA -> Fun CB CB -> Fun CC CC -> Fun (CA *#* (CB *#* CC)) CA -> + ChebFun N CB CA -> ChebFun N CB CB -> ChebFun N CB CC -> Property +prop_ChebFun1_Applicative_assoc' + (Fn f) (Fn g) (Fn h) (Fn i) xs ys zs = + uncurry (===) (getEqual (law_Applicative_assoc' + (UFun f) (UFun g) (UFun h) (UFun i) + xs ys zs)) diff --git a/test/ChebFunSpec.hs b/test/ChebFunSpec.hs index 252871a..c79d879 100644 --- a/test/ChebFunSpec.hs +++ b/test/ChebFunSpec.hs @@ -11,12 +11,13 @@ import Prelude hiding ( id, (.), curry, uncurry ) import Data.Fixed -import Test.QuickCheck +import Test.QuickCheck hiding (scale) import Applicative import Category import ChebFun import Foldable +import qualified Function as F import Functor import Unbox @@ -25,84 +26,105 @@ import Unbox maxabs :: (Foldable f, k ~ Dom f, Obj k a, Num a, Ord a) => f a -> a maxabs = foldl (\r -> max r . abs) 0 +sumabs :: (Foldable f, k ~ Dom f, Obj k a, Num a) => f a -> a +sumabs = foldl (\r -> (r +) . abs) 0 + -- | Approximate floating-point equality approx :: (RealFrac a, Show a) => a -> a -> a -> Property -approx eps x y = counterexample +approx delta x y = counterexample (show x ++ " /= " ++ show y ++ "; " ++ "absdiff = " ++ show (y - x) ++ ", " ++ "reldiff = " ++ show ((y - x) / maxabs [x, y]) ++ "; " ++ - "eps = " ++ show eps) - (abs (x - y) <= eps) + "delta = " ++ show delta) + (abs (x - y) <= delta) + +limitArg :: Real a => a -> a +limitArg x = mod' (x + 1) 2 - 1 +limitFun :: ( Functor v, Cod v ~ (->), Foldable v + , Obj (Dom v) b, Fractional b, Ord b) + => ChebFun v a b -> ChebFun v a b +limitFun f = let scale = max 1 (sumabs (getChebFun f)) + in ChebFun (fmap (unapply (/ scale)) (getChebFun f)) type N = 10 -type N' = 2 +type V = NUVector N type CA = Double type CB = Double type CC = Double -eps1 :: Double -eps1 = 1.0e-13 +eps :: Double +eps = 1.0e-12 -prop_ChebFun_Category_comp_id_left :: ChebFun N CA CB -> CA -> Property -prop_ChebFun_Category_comp_id_left f x' = - let x = mod' (x' + 1) 2 - 1 - in uncurry (approx eps) (getFnEqual (law_Category_comp_id_left f) x) - where eps = let (ChebFun cs) = f in eps1 * max 1 (maxabs cs) +prop_ChebFun_compact_inv :: CA -> Property +prop_ChebFun_compact_inv x' = + uncurry (approx eps) (getFnEqual law_ChebFun_compact_inv x) + where x = compactify x' -prop_ChebFun_Category_comp_id_right :: ChebFun N CA CB -> CA -> Property -prop_ChebFun_Category_comp_id_right f x' = - let x = mod' x' 2 - 1 - in uncurry (approx eps) (getFnEqual (law_Category_comp_id_right f) x) - where eps = let (ChebFun cs) = f in eps1 * max 1 (maxabs cs) +prop_ChebFun_compact_inv' :: CA -> Property +prop_ChebFun_compact_inv' x' = + uncurry (approx eps) (getFnEqual law_ChebFun_compact_inv x) + where x = mod' x' 2 - 1 --- TODO: Make this work for larger N' -prop_ChebFun_Category_comp_assoc :: - ChebFun N' CA CC -> ChebFun N' CB CA -> ChebFun N' CA CB -> CA -> Property -prop_ChebFun_Category_comp_assoc h g f x' = - let x = mod' x' 2 - 1 - in uncurry (approx eps) (getFnEqual (law_Category_comp_assoc h g f) x) - where eps = let (ChebFun csf) = f - (ChebFun csg) = g - (ChebFun csh) = h - in eps1 * maxabs csf * maxabs csg * maxabs csh - -prop_ChebFun_Functor_id :: ChebFun N CB CA -> Property +prop_ChebFun_Functor_id :: ChebFun V CB CA -> Property prop_ChebFun_Functor_id xs = uncurry (===) (getFnEqual law_Functor_id xs) prop_ChebFun_Functor_assoc :: - Fun CB CC -> Fun CA CB -> ChebFun N CB CA -> Property + Fun CB CC -> Fun CA CB -> ChebFun V CB CA -> Property prop_ChebFun_Functor_assoc (Fn g) (Fn f) xs = uncurry (===) (getFnEqual (law_Functor_assoc (UFun g) (UFun f)) xs) -prop_ChebFun_Applicative_id' :: Fun CA CB -> ChebFun N CB CA -> Property +prop_ChebFun_Applicative_id' :: Fun CA CB -> ChebFun V CB CA -> Property prop_ChebFun_Applicative_id' (Fn f) xs = uncurry (===) (getEqual (law_Applicative_id' (UFun f) xs)) prop_ChebFun_Applicative_id_left' :: - Fun (CA *#* CB) CC -> CA -> ChebFun N CB CB -> Property + Fun (CA *#* CB) CC -> CA -> ChebFun V CB CB -> Property prop_ChebFun_Applicative_id_left' (Fn f) x ys = uncurry (===) (getEqual (law_Applicative_id_left' (UFun f) x ys)) prop_ChebFun_Applicative_id_right' :: - Fun (CA *#* CB) CC -> ChebFun N CB CA -> CB -> Property + Fun (CA *#* CB) CC -> ChebFun V CB CA -> CB -> Property prop_ChebFun_Applicative_id_right' (Fn f) xs y = uncurry (===) (getEqual (law_Applicative_id_right' (UFun f) xs y)) prop_ChebFun_Applicative_assoc' :: Fun CA CA -> Fun CB CB -> Fun CC CC -> Fun (CA *#* (CB *#* CC)) CA -> - ChebFun N CB CA -> ChebFun N CB CB -> ChebFun N CB CC -> Property + ChebFun V CB CA -> ChebFun V CB CB -> ChebFun V CB CC -> Property prop_ChebFun_Applicative_assoc' (Fn f) (Fn g) (Fn h) (Fn i) xs ys zs = uncurry (===) (getEqual (law_Applicative_assoc' (UFun f) (UFun g) (UFun h) (UFun i) xs ys zs)) + + + +prop_ChebVal_MCategory_comp_id_left :: ChebFun V CA CB -> CA -> Property +prop_ChebVal_MCategory_comp_id_left f' x' = + uncurry (===) (F.getFnEqual (F.law_MCategory_comp_id_left f) x) + where x = limitArg x' + f = limitFun f' + +prop_ChebVal_MCategory_comp_id_right :: ChebFun V CA CB -> CA -> Property +prop_ChebVal_MCategory_comp_id_right f' x' = + uncurry (===) (F.getFnEqual (F.law_MCategory_comp_id_right f) x) + where x = limitArg x' + f = limitFun f' + +prop_ChebVal_MCategory_comp_assoc :: + ChebFun V CA CC -> ChebFun V CB CA -> ChebFun V CA CB -> CA -> Property +prop_ChebVal_MCategory_comp_assoc h' g' f' x' = + uncurry (===) (F.getFnEqual (F.law_MCategory_comp_assoc h g f) x) + where x = limitArg x' + f = limitFun f' + g = limitFun g' + h = limitFun h' diff --git a/test/FunctionSpec.hs b/test/FunctionSpec.hs index 796d6c2..b7aba5d 100644 --- a/test/FunctionSpec.hs +++ b/test/FunctionSpec.hs @@ -53,32 +53,3 @@ prop_Num_MCategory_comp_assoc :: prop_Num_MCategory_comp_assoc (Fn h) (Fn g) (Fn f) x = uncurry (===) (getFnEqual (law_MCategory_comp_assoc (NFun h) (NFun g) (NFun f)) x) - - - -type N = 5 -- 10 - -type CA = Double -type CB = Double -type CC = Double - -eps1 :: Double -eps1 = 1.0e-13 - - - -prop_ChebVal_MCategory_comp_id_left :: ChebFun N CA CB -> CA -> Property -prop_ChebVal_MCategory_comp_id_left f x' = - let x = mod' (x' + 1) 2 - 1 - in uncurry (===) (getFnEqual (law_MCategory_comp_id_left f) x) - -prop_ChebVal_MCategory_comp_id_right :: ChebFun N CA CB -> CA -> Property -prop_ChebVal_MCategory_comp_id_right f x' = - let x = mod' (x' + 1) 2 - 1 - in uncurry (===) (getFnEqual (law_MCategory_comp_id_right f) x) - -prop_ChebVal_MCategory_comp_assoc :: - ChebFun N CA CC -> ChebFun N CB CA -> ChebFun N CA CB -> CA -> Property -prop_ChebVal_MCategory_comp_assoc h g f x' = - let x = mod' (x' + 1) 2 - 1 - in uncurry (===) (getFnEqual (law_MCategory_comp_assoc h g f) x)