text:
-- | Type classes for random generation of values.
--
-- __Note__: the contents of this module are re-exported by
-- "Test.QuickCheck". You do not need to import it directly.
{-# LANGUAGE CPP #-}
{-# LANGUAGE FlexibleContexts #-}
#ifndef NO_GENERICS
{-# LANGUAGE DefaultSignatures, FlexibleContexts, TypeOperators #-}
{-# LANGUAGE FlexibleInstances, KindSignatures, ScopedTypeVariables #-}
{-# LANGUAGE MultiParamTypeClasses #-}
#if __GLASGOW_HASKELL__ >= 710
#define OVERLAPPING_ {-# OVERLAPPING #-}
#else
{-# LANGUAGE OverlappingInstances  #-}
#define OVERLAPPING_
#endif
#endif
#ifndef NO_POLYKINDS
{-# LANGUAGE PolyKinds #-}
#endif
#ifndef NO_SAFE_HASKELL
{-# LANGUAGE Trustworthy #-}
#endif
#ifndef NO_NEWTYPE_DERIVING
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
#endif
module Test.QuickCheck.Arbitrary
  (
  -- * Arbitrary and CoArbitrary classes
    Arbitrary(..)
  , CoArbitrary(..)

  -- ** Unary and Binary classes
  , Arbitrary1(..)
  , arbitrary1
  , shrink1
  , Arbitrary2(..)
  , arbitrary2
  , shrink2

  -- ** Helper functions for implementing arbitrary
  , applyArbitrary2
  , applyArbitrary3
  , applyArbitrary4
  , arbitrarySizedIntegral        -- :: Integral a => Gen a
  , arbitrarySizedNatural         -- :: Integral a => Gen a
  , arbitraryBoundedIntegral      -- :: (Bounded a, Integral a) => Gen a
  , arbitrarySizedBoundedIntegral -- :: (Bounded a, Integral a) => Gen a
  , arbitrarySizedFractional      -- :: Fractional a => Gen a
  , arbitraryBoundedRandom        -- :: (Bounded a, Random a) => Gen a
  , arbitraryBoundedEnum          -- :: (Bounded a, Enum a) => Gen a
  -- ** Generators for various kinds of character
  , arbitraryUnicodeChar   -- :: Gen Char
  , arbitraryASCIIChar     -- :: Gen Char
  , arbitraryPrintableChar -- :: Gen Char
  -- ** Helper functions for implementing shrink
#ifndef NO_GENERICS
  , RecursivelyShrink
  , GSubterms
  , genericShrink      -- :: (Generic a, Arbitrary a, RecursivelyShrink (Rep a), GSubterms (Rep a) a) => a -> [a]
  , subterms           -- :: (Generic a, Arbitrary a, GSubterms (Rep a) a) => a -> [a]
  , recursivelyShrink  -- :: (Generic a, RecursivelyShrink (Rep a)) => a -> [a]
  , genericCoarbitrary -- :: (Generic a, GCoArbitrary (Rep a)) => a -> Gen b -> Gen b
#endif
  , shrinkNothing            -- :: a -> [a]
  , shrinkList               -- :: (a -> [a]) -> [a] -> [[a]]
  , shrinkMap                -- :: Arbitrary a -> (a -> b) -> (b -> a) -> b -> [b]
  , shrinkMapBy              -- :: (a -> b) -> (b -> a) -> (a -> [a]) -> b -> [b]
  , shrinkIntegral           -- :: Integral a => a -> [a]
  , shrinkRealFrac           -- :: RealFrac a => a -> [a]
  , shrinkDecimal            -- :: RealFrac a => a -> [a]
  -- ** Helper functions for implementing coarbitrary
  , coarbitraryIntegral      -- :: Integral a => a -> Gen b -> Gen b
  , coarbitraryReal          -- :: Real a => a -> Gen b -> Gen b
  , coarbitraryShow          -- :: Show a => a -> Gen b -> Gen b
  , coarbitraryEnum          -- :: Enum a => a -> Gen b -> Gen b
  , (><)

  -- ** Generators which use arbitrary
  , vector       -- :: Arbitrary a => Int -> Gen [a]
  , orderedList  -- :: (Ord a, Arbitrary a) => Gen [a]
  , infiniteList -- :: Arbitrary a => Gen [a]
  )
 where

--------------------------------------------------------------------------
-- imports

import Control.Applicative
import Data.Foldable(toList)
import System.Random(Random)
import Test.QuickCheck.Gen
import Test.QuickCheck.Random
import Test.QuickCheck.Gen.Unsafe

{-
import Data.Generics
  ( (:*:)(..)
  , (:+:)(..)
  , Unit(..)
  )
-}

import Data.Char
  ( ord
  , isLower
  , isUpper
  , toLower
  , isDigit
  , isSpace
  , isPrint
  , generalCategory
  , GeneralCategory(..)
  )

#ifndef NO_FIXED
import Data.Fixed
  ( Fixed
  , HasResolution
  )
#endif

import Data.Ratio
  ( Ratio
  , (%)
  , numerator
  , denominator
  )

import Data.Complex
  ( Complex((:+)) )

import Data.List
  ( sort
  , nub
  )

import Data.Version (Version (..))

#if defined(MIN_VERSION_base)
#if MIN_VERSION_base(4,2,0)
import System.IO
  ( Newline(..)
  , NewlineMode(..)
  )
#endif
#endif

import Control.Monad
  ( liftM
  , liftM2
  , liftM3
  , liftM4
  , liftM5
  )

import Data.Int(Int8, Int16, Int32, Int64)
import Data.Word(Word, Word8, Word16, Word32, Word64)
import System.Exit (ExitCode(..))
import Foreign.C.Types

#ifndef NO_GENERICS
import GHC.Generics
#endif

import qualified Data.Set as Set
import qualified Data.Map as Map
import qualified Data.IntSet as IntSet
import qualified Data.IntMap as IntMap
import qualified Data.Sequence as Sequence
import qualified Data.Tree as Tree
import Data.Bits

import qualified Data.Monoid as Monoid

#ifndef NO_TRANSFORMERS
import Data.Functor.Identity
import Data.Functor.Constant
import Data.Functor.Compose
import Data.Functor.Product
#endif

--------------------------------------------------------------------------
-- ** class Arbitrary

-- | Random generation and shrinking of values.
--
-- QuickCheck provides @Arbitrary@ instances for most types in @base@,
-- except those which incur extra dependencies.
-- For a wider range of @Arbitrary@ instances see the
-- <http://hackage.haskell.org/package/quickcheck-instances quickcheck-instances>
-- package.
class Arbitrary a where
  -- | A generator for values of the given type.
  --
  -- It is worth spending time thinking about what sort of test data
  -- you want - good generators are often the difference between
  -- finding bugs and not finding them. You can use 'sample',
  -- 'label' and 'classify' to check the quality of your test data.
  --
  -- There is no generic @arbitrary@ implementation included because we don't
  -- know how to make a high-quality one. If you want one, consider using the
  -- <http://hackage.haskell.org/package/testing-feat testing-feat> or
  -- <http://hackage.haskell.org/package/generic-random generic-random> packages.
  --
  -- The <http://www.cse.chalmers.se/~rjmh/QuickCheck/manual.html QuickCheck manual>
  -- goes into detail on how to write good generators. Make sure to look at it,
  -- especially if your type is recursive!
  arbitrary :: Gen a

  -- | Produces a (possibly) empty list of all the possible
  -- immediate shrinks of the given value.
  --
  -- The default implementation returns the empty list, so will not try to
  -- shrink the value. If your data type has no special invariants, you can
  -- enable shrinking by defining @shrink = 'genericShrink'@, but by customising
  -- the behaviour of @shrink@ you can often get simpler counterexamples.
  --
  -- Most implementations of 'shrink' should try at least three things:
  --
  -- 1. Shrink a term to any of its immediate subterms.
  --    You can use 'subterms' to do this.
  --
  -- 2. Recursively apply 'shrink' to all immediate subterms.
  --    You can use 'recursivelyShrink' to do this.
  --
  -- 3. Type-specific shrinkings such as replacing a constructor by a
  --    simpler constructor.
  --
  -- For example, suppose we have the following implementation of binary trees:
  --
  -- > data Tree a = Nil | Branch a (Tree a) (Tree a)
  --
  -- We can then define 'shrink' as follows:
  --
  -- > shrink Nil = []
  -- > shrink (Branch x l r) =
  -- >   -- shrink Branch to Nil
  -- >   [Nil] ++
  -- >   -- shrink to subterms
  -- >   [l, r] ++
  -- >   -- recursively shrink subterms
  -- >   [Branch x' l' r' | (x', l', r') <- shrink (x, l, r)]
  --
  -- There are a couple of subtleties here:
  --
  -- * QuickCheck tries the shrinking candidates in the order they
  --   appear in the list, so we put more aggressive shrinking steps
  --   (such as replacing the whole tree by @Nil@) before smaller
  --   ones (such as recursively shrinking the subtrees).
  --
  -- * It is tempting to write the last line as
  --   @[Branch x' l' r' | x' <- shrink x, l' <- shrink l, r' <- shrink r]@
  --   but this is the /wrong thing/! It will force QuickCheck to shrink
  --   @x@, @l@ and @r@ in tandem, and shrinking will stop once /one/ of
  --   the three is fully shrunk.
  --
  -- There is a fair bit of boilerplate in the code above.
  -- We can avoid it with the help of some generic functions.
  -- The function 'genericShrink' tries shrinking a term to all of its
  -- subterms and, failing that, recursively shrinks the subterms.
  -- Using it, we can define 'shrink' as:
  --
  -- > shrink x = shrinkToNil x ++ genericShrink x
  -- >   where
  -- >     shrinkToNil Nil = []
  -- >     shrinkToNil (Branch _ l r) = [Nil]
  --
  -- 'genericShrink' is a combination of 'subterms', which shrinks
  -- a term to any of its subterms, and 'recursivelyShrink', which shrinks
  -- all subterms of a term. These may be useful if you need a bit more
  -- control over shrinking than 'genericShrink' gives you.
  --
  -- A final gotcha: we cannot define 'shrink' as simply @'shrink' x = Nil:'genericShrink' x@
  -- as this shrinks @Nil@ to @Nil@, and shrinking will go into an
  -- infinite loop.
  --
  -- If all this leaves you bewildered, you might try @'shrink' = 'genericShrink'@ to begin with,
  -- after deriving @Generic@ for your type. However, if your data type has any
  -- special invariants, you will need to check that 'genericShrink' can't break those invariants.
  shrink :: a -> [a]
  shrink _ = []

-- | Lifting of the 'Arbitrary' class to unary type constructors.
class Arbitrary1 f where
  liftArbitrary :: Gen a -> Gen (f a)
  liftShrink    :: (a -> [a]) -> f a -> [f a]
  liftShrink _ _ = []

arbitrary1 :: (Arbitrary1 f, Arbitrary a) => Gen (f a)
arbitrary1 = liftArbitrary arbitrary

shrink1 :: (Arbitrary1 f, Arbitrary a) => f a -> [f a]
shrink1 = liftShrink shrink

-- | Lifting of the 'Arbitrary' class to binary type constructors.
class Arbitrary2 f where
  liftArbitrary2 :: Gen a -> Gen b -> Gen (f a b)
  liftShrink2    :: (a -> [a]) -> (b -> [b]) -> f a b -> [f a b]
  liftShrink2 _ _ _ = []

arbitrary2 :: (Arbitrary2 f, Arbitrary a, Arbitrary b) => Gen (f a b)
arbitrary2 = liftArbitrary2 arbitrary arbitrary

shrink2 :: (Arbitrary2 f, Arbitrary a, Arbitrary b) => f a b -> [f a b]
shrink2 = liftShrink2 shrink shrink

#ifndef NO_GENERICS
-- | Shrink a term to any of its immediate subterms,
-- and also recursively shrink all subterms.
genericShrink :: (Generic a, RecursivelyShrink (Rep a), GSubterms (Rep a) a) => a -> [a]
genericShrink x = subterms x ++ recursivelyShrink x

-- | Recursively shrink all immediate subterms.
recursivelyShrink :: (Generic a, RecursivelyShrink (Rep a)) => a -> [a]
recursivelyShrink = map to . grecursivelyShrink . from

class RecursivelyShrink f where
  grecursivelyShrink :: f a -> [f a]

instance (RecursivelyShrink f, RecursivelyShrink g) => RecursivelyShrink (f :*: g) where
  grecursivelyShrink (x :*: y) =
    [x' :*: y | x' <- grecursivelyShrink x] ++
    [x :*: y' | y' <- grecursivelyShrink y]

instance (RecursivelyShrink f, RecursivelyShrink g) => RecursivelyShrink (f :+: g) where
  grecursivelyShrink (L1 x) = map L1 (grecursivelyShrink x)
  grecursivelyShrink (R1 x) = map R1 (grecursivelyShrink x)

instance RecursivelyShrink f => RecursivelyShrink (M1 i c f) where
  grecursivelyShrink (M1 x) = map M1 (grecursivelyShrink x)

instance Arbitrary a => RecursivelyShrink (K1 i a) where
  grecursivelyShrink (K1 x) = map K1 (shrink x)

instance RecursivelyShrink U1 where
  grecursivelyShrink U1 = []

instance RecursivelyShrink V1 where
  -- The empty type can't be shrunk to anything.
  grecursivelyShrink _ = []

-- | All immediate subterms of a term.
subterms :: (Generic a, GSubterms (Rep a) a) => a -> [a]
subterms = gSubterms . from

class GSubterms f a where
  -- | Provides the immediate subterms of a term that are of the same type
  -- as the term itself.
  --
  -- Requires a constructor to be stripped off; this means it skips through
  -- @M1@ wrappers and returns @[]@ on everything that's not `(:*:)` or `(:+:)`.
  --
  -- Once a `(:*:)` or `(:+:)` constructor has been reached, this function
  -- delegates to `gSubtermsIncl` to return the immediately next constructor
  -- available.
  gSubterms :: f a -> [a]

instance GSubterms V1 a where
  -- The empty type can't be shrunk to anything.
  gSubterms _ = []

instance GSubterms U1 a where
  gSubterms U1 = []

instance (GSubtermsIncl f a, GSubtermsIncl g a) => GSubterms (f :*: g) a where
  gSubterms (l :*: r) = gSubtermsIncl l ++ gSubtermsIncl r

instance (GSubtermsIncl f a, GSubtermsIncl g a) => GSubterms (f :+: g) a where
  gSubterms (L1 x) = gSubtermsIncl x
  gSubterms (R1 x) = gSubtermsIncl x

instance GSubterms f a => GSubterms (M1 i c f) a where
  gSubterms (M1 x) = gSubterms x

instance GSubterms (K1 i a) b where
  gSubterms (K1 _) = []

class GSubtermsIncl f a where
  -- | Provides the immediate subterms of a term that are of the same type
  -- as the term itself.
  --
  -- In contrast to `gSubterms`, this returns the immediate next constructor
  -- available.
  gSubtermsIncl :: f a -> [a]

instance GSubtermsIncl V1 a where
  -- The empty type can't be shrunk to anything.
  gSubtermsIncl _ = []

instance GSubtermsIncl U1 a where
  gSubtermsIncl U1 = []

instance (GSubtermsIncl f a, GSubtermsIncl g a) => GSubtermsIncl (f :*: g) a where
  gSubtermsIncl (l :*: r) = gSubtermsIncl l ++ gSubtermsIncl r

instance (GSubtermsIncl f a, GSubtermsIncl g a) => GSubtermsIncl (f :+: g) a where
  gSubtermsIncl (L1 x) = gSubtermsIncl x
  gSubtermsIncl (R1 x) = gSubtermsIncl x

instance GSubtermsIncl f a => GSubtermsIncl (M1 i c f) a where
  gSubtermsIncl (M1 x) = gSubtermsIncl x

-- This is the important case: We've found a term of the same type.
instance OVERLAPPING_ GSubtermsIncl (K1 i a) a where
  gSubtermsIncl (K1 x) = [x]

instance OVERLAPPING_ GSubtermsIncl (K1 i a) b where
  gSubtermsIncl (K1 _) = []

#endif

-- instances

instance (CoArbitrary a) => Arbitrary1 ((->) a) where
  liftArbitrary arbB = promote (`coarbitrary` arbB)

instance (CoArbitrary a, Arbitrary b) => Arbitrary (a -> b) where
  arbitrary = arbitrary1

instance Arbitrary () where
  arbitrary = return ()

instance Arbitrary Bool where
  arbitrary = chooseEnum (False,True)
  shrink True = [False]
  shrink False = []

instance Arbitrary Ordering where
  arbitrary = elements [LT, EQ, GT]
  shrink GT = [EQ, LT]
  shrink LT = [EQ]
  shrink EQ = []

instance Arbitrary1 Maybe where
  liftArbitrary arb = frequency [(1, return Nothing), (3, liftM Just arb)]

  liftShrink shr (Just x) = Nothing : [ Just x' | x' <- shr x ]
  liftShrink _   Nothing  = []

instance Arbitrary a => Arbitrary (Maybe a) where
  arbitrary = arbitrary1
  shrink = shrink1

instance Arbitrary2 Either where
  liftArbitrary2 arbA arbB = oneof [liftM Left arbA, liftM Right arbB]

  liftShrink2 shrA _ (Left x)  = [ Left  x' | x' <- shrA x ]
  liftShrink2 _ shrB (Right y) = [ Right y' | y' <- shrB y ]

instance Arbitrary a => Arbitrary1 (Either a) where
  liftArbitrary = liftArbitrary2 arbitrary
  liftShrink = liftShrink2 shrink

instance (Arbitrary a, Arbitrary b) => Arbitrary (Either a b) where
  arbitrary = arbitrary2
  shrink = shrink2

instance Arbitrary1 [] where
  liftArbitrary = listOf
  liftShrink = shrinkList

instance Arbitrary a => Arbitrary [a] where
  arbitrary = arbitrary1
  shrink = shrink1

-- | Shrink a list of values given a shrinking function for individual values.
shrinkList :: (a -> [a]) -> [a] -> [[a]]
shrinkList shr xs = concat [ removes k n xs | k <- takeWhile (>0) (iterate (`div`2) n) ]
                 ++ shrinkOne xs
 where
  n = length xs

  shrinkOne []     = []
  shrinkOne (x:xs) = [ x':xs | x'  <- shr x ]
                  ++ [ x:xs' | xs' <- shrinkOne xs ]

  removes k n xs
    | k > n     = []
    | null xs2  = [[]]
    | otherwise = xs2 : map (xs1 ++) (removes k (n-k) xs2)
   where
    xs1 = take k xs
    xs2 = drop k xs

{-
  -- "standard" definition for lists:
  shrink []     = []
  shrink (x:xs) = [ xs ]
               ++ [ x:xs' | xs' <- shrink xs ]
               ++ [ x':xs | x'  <- shrink x ]
-}

instance Integral a => Arbitrary (Ratio a) where
  arbitrary = arbitrarySizedFractional
  shrink    = shrinkRealFrac

#if defined(MIN_VERSION_base) && MIN_VERSION_base(4,4,0)
instance Arbitrary a => Arbitrary (Complex a) where
#else
instance (RealFloat a, Arbitrary a) => Arbitrary (Complex a) where
#endif
  arbitrary = liftM2 (:+) arbitrary arbitrary
  shrink (x :+ y) = [ x' :+ y | x' <- shrink x ] ++
                    [ x :+ y' | y' <- shrink y ]

#ifndef NO_FIXED
instance HasResolution a => Arbitrary (Fixed a) where
  arbitrary = arbitrarySizedFractional
  shrink    = shrinkDecimal
#endif

instance Arbitrary2 (,) where
  liftArbitrary2 = liftM2 (,)
  liftShrink2 shrA shrB (x, y) =
       [ (x', y) | x' <- shrA x ]
    ++ [ (x, y') | y' <- shrB y ]

instance (Arbitrary a) => Arbitrary1 ((,) a) where
  liftArbitrary = liftArbitrary2 arbitrary
  liftShrink = liftShrink2 shrink

instance (Arbitrary a, Arbitrary b) => Arbitrary (a,b) where
  arbitrary = arbitrary2
  shrink = shrink2

instance (Arbitrary a, Arbitrary b, Arbitrary c)
      => Arbitrary (a,b,c)
 where
  arbitrary = liftM3 (,,) arbitrary arbitrary arbitrary

  shrink (x, y, z) =
    [ (x', y', z')
    | (x', (y', z')) <- shrink (x, (y, z)) ]

instance (Arbitrary a, Arbitrary b, Arbitrary c, Arbitrary d)
      => Arbitrary (a,b,c,d)
 where
  arbitrary = liftM4 (,,,) arbitrary arbitrary arbitrary arbitrary

  shrink (w, x, y, z) =
    [ (w', x', y', z')
    | (w', (x', (y', z'))) <- shrink (w, (x, (y, z))) ]

instance (Arbitrary a, Arbitrary b, Arbitrary c, Arbitrary d, Arbitrary e)
      => Arbitrary (a,b,c,d,e)
 where
  arbitrary = liftM5 (,,,,) arbitrary arbitrary arbitrary arbitrary arbitrary

  shrink (v, w, x, y, z) =
    [ (v', w', x', y', z')
    | (v', (w', (x', (y', z')))) <- shrink (v, (w, (x, (y, z)))) ]

instance ( Arbitrary a, Arbitrary b, Arbitrary c, Arbitrary d, Arbitrary e
         , Arbitrary f
         )
      => Arbitrary (a,b,c,d,e,f)
 where
  arbitrary = return (,,,,,)
          <*> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary
          <*> arbitrary <*> arbitrary

  shrink (u, v, w, x, y, z) =
    [ (u', v', w', x', y', z')
    | (u', (v', (w', (x', (y', z'))))) <- shrink (u, (v, (w, (x, (y, z))))) ]

instance ( Arbitrary a, Arbitrary b, Arbitrary c, Arbitrary d, Arbitrary e
         , Arbitrary f, Arbitrary g
         )
      => Arbitrary (a,b,c,d,e,f,g)
 where
  arbitrary = return (,,,,,,)
          <*> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary
          <*> arbitrary <*> arbitrary <*> arbitrary

  shrink (t, u, v, w, x, y, z) =
    [ (t', u', v', w', x', y', z')
    | (t', (u', (v', (w', (x', (y', z')))))) <- shrink (t, (u, (v, (w, (x, (y, z)))))) ]

instance ( Arbitrary a, Arbitrary b, Arbitrary c, Arbitrary d, Arbitrary e
         , Arbitrary f, Arbitrary g, Arbitrary h
         )
      => Arbitrary (a,b,c,d,e,f,g,h)
 where
  arbitrary = return (,,,,,,,)
          <*> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary
          <*> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary

  shrink (s, t, u, v, w, x, y, z) =
    [ (s', t', u', v', w', x', y', z')
    | (s', (t', (u', (v', (w', (x', (y', z')))))))
      <- shrink (s, (t, (u, (v, (w, (x, (y, z))))))) ]

instance ( Arbitrary a, Arbitrary b, Arbitrary c, Arbitrary d, Arbitrary e
         , Arbitrary f, Arbitrary g, Arbitrary h, Arbitrary i
         )
      => Arbitrary (a,b,c,d,e,f,g,h,i)
 where
  arbitrary = return (,,,,,,,,)
          <*> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary
          <*> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary
          <*> arbitrary

  shrink (r, s, t, u, v, w, x, y, z) =
    [ (r', s', t', u', v', w', x', y', z')
    | (r', (s', (t', (u', (v', (w', (x', (y', z'))))))))
      <- shrink (r, (s, (t, (u, (v, (w, (x, (y, z)))))))) ]

instance ( Arbitrary a, Arbitrary b, Arbitrary c, Arbitrary d, Arbitrary e
         , Arbitrary f, Arbitrary g, Arbitrary h, Arbitrary i, Arbitrary j
         )
      => Arbitrary (a,b,c,d,e,f,g,h,i,j)
 where
  arbitrary = return (,,,,,,,,,)
          <*> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary
          <*> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary
          <*> arbitrary <*> arbitrary

  shrink (q, r, s, t, u, v, w, x, y, z) =
    [ (q', r', s', t', u', v', w', x', y', z')
    | (q', (r', (s', (t', (u', (v', (w', (x', (y', z')))))))))
      <- shrink (q, (r, (s, (t, (u, (v, (w, (x, (y, z))))))))) ]

-- typical instance for primitive (numerical) types

instance Arbitrary Integer where
  arbitrary = arbitrarySizedIntegral
  shrink    = shrinkIntegral

instance Arbitrary Int where
  arbitrary = arbitrarySizedIntegral
  shrink    = shrinkIntegral

instance Arbitrary Int8 where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink    = shrinkIntegral

instance Arbitrary Int16 where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink    = shrinkIntegral

instance Arbitrary Int32 where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink    = shrinkIntegral

instance Arbitrary Int64 where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink    = shrinkIntegral

instance Arbitrary Word where
  arbitrary = arbitrarySizedNatural
  shrink    = shrinkIntegral

instance Arbitrary Word8 where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink    = shrinkIntegral

instance Arbitrary Word16 where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink    = shrinkIntegral

instance Arbitrary Word32 where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink    = shrinkIntegral

instance Arbitrary Word64 where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink    = shrinkIntegral

instance Arbitrary Char where
  arbitrary =
    frequency
      [(3, arbitraryASCIIChar),
       (1, arbitraryUnicodeChar)]

  shrink c = filter (<. c) $ nub
            $ ['a','b','c']
            ++ [ toLower c | isUpper c ]
            ++ ['A','B','C']
            ++ ['1','2','3']
            ++ [' ','\n']
     where
      a <. b  = stamp a < stamp b
      stamp a = ( (not (isLower a)
                , not (isUpper a)
                , not (isDigit a))
                , (not (a==' ')
                , not (isSpace a)
                , a)
                )

instance Arbitrary Float where
  arbitrary = oneof
    -- generate 0..1 numbers with full precision
    [ genFloat
    -- generate integral numbers
    , fromIntegral <$> (arbitrary :: Gen Int)
    -- generate fractions with small denominators
    , smallDenominators
    -- uniform -size..size with with denominators ~ size
    , uniform
    -- and uniform -size..size with higher precision
    , arbitrarySizedFractional
    ]
    where
      smallDenominators = sized $ \n -> do
        i <- chooseInt (0, n)
        pure (fromRational (streamNth i rationalUniverse))

      uniform = sized $ \n -> do
        let n' = toInteger n
        b <- chooseInteger (1, max 1 n')
        a <- chooseInteger ((-n') * b, n' * b)
        return (fromRational (a % b))

  shrink    = shrinkDecimal

instance Arbitrary Double where
  arbitrary = oneof
    -- generate 0..1 numbers with full precision
    [ genDouble
    -- generate integral numbers
    , fromIntegral <$> (arbitrary :: Gen Int)
    -- generate fractions with small denominators
    , smallDenominators
    -- uniform -size..size with with denominators ~ size
    , uniform
    -- and uniform -size..size with higher precision
    , arbitrarySizedFractional
    ]
    where
      smallDenominators = sized $ \n -> do
        i <- chooseInt (0, n)
        pure (fromRational (streamNth i rationalUniverse))

      uniform = sized $ \n -> do
        let n' = toInteger n
        b <- chooseInteger (1, max 1 n')
        a <- chooseInteger ((-n') * b, n' * b)
        return (fromRational (a % b))

  shrink    = shrinkDecimal

instance Arbitrary CChar where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink = shrinkIntegral

instance Arbitrary CSChar where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink = shrinkIntegral

instance Arbitrary CUChar where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink = shrinkIntegral

instance Arbitrary CShort where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink = shrinkIntegral

instance Arbitrary CUShort where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink = shrinkIntegral

instance Arbitrary CInt where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink = shrinkIntegral

instance Arbitrary CUInt where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink = shrinkIntegral

instance Arbitrary CLong where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink = shrinkIntegral

instance Arbitrary CULong where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink = shrinkIntegral

instance Arbitrary CPtrdiff where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink = shrinkIntegral

instance Arbitrary CSize where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink = shrinkIntegral

instance Arbitrary CWchar where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink = shrinkIntegral

instance Arbitrary CSigAtomic where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink = shrinkIntegral

instance Arbitrary CLLong where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink = shrinkIntegral

instance Arbitrary CULLong where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink = shrinkIntegral

instance Arbitrary CIntPtr where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink = shrinkIntegral

instance Arbitrary CUIntPtr where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink = shrinkIntegral

instance Arbitrary CIntMax where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink = shrinkIntegral

instance Arbitrary CUIntMax where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink = shrinkIntegral

#ifndef NO_CTYPES_CONSTRUCTORS
-- The following four types have no Bounded instance,
-- so we fake it by discovering the bounds at runtime.
instance Arbitrary CClock where
  arbitrary = fmap CClock arbitrary
  shrink (CClock x) = map CClock (shrink x)

instance Arbitrary CTime where
  arbitrary = fmap CTime arbitrary
  shrink (CTime x) = map CTime (shrink x)

#ifndef NO_FOREIGN_C_USECONDS
instance Arbitrary CUSeconds where
  arbitrary = fmap CUSeconds arbitrary
  shrink (CUSeconds x) = map CUSeconds (shrink x)

instance Arbitrary CSUSeconds where
  arbitrary = fmap CSUSeconds arbitrary
  shrink (CSUSeconds x) = map CSUSeconds (shrink x)
#endif
#endif

instance Arbitrary CFloat where
  arbitrary = arbitrarySizedFractional
  shrink = shrinkDecimal

instance Arbitrary CDouble where
  arbitrary = arbitrarySizedFractional
  shrink = shrinkDecimal

-- Arbitrary instances for container types
instance (Ord a, Arbitrary a) => Arbitrary (Set.Set a) where
  arbitrary = fmap Set.fromList arbitrary
  shrink = map Set.fromList . shrink . Set.toList
instance (Ord k, Arbitrary k) => Arbitrary1 (Map.Map k) where
  liftArbitrary = fmap Map.fromList . liftArbitrary . liftArbitrary
  liftShrink shr = map Map.fromList . liftShrink (liftShrink shr) . Map.toList
instance (Ord k, Arbitrary k, Arbitrary v) => Arbitrary (Map.Map k v) where
  arbitrary = arbitrary1
  shrink = shrink1
instance Arbitrary IntSet.IntSet where
  arbitrary = fmap IntSet.fromList arbitrary
  shrink = map IntSet.fromList . shrink . IntSet.toList
instance Arbitrary1 IntMap.IntMap where
  liftArbitrary = fmap IntMap.fromList . liftArbitrary . liftArbitrary
  liftShrink shr = map IntMap.fromList . liftShrink (liftShrink shr) . IntMap.toList
instance Arbitrary a => Arbitrary (IntMap.IntMap a) where
  arbitrary = arbitrary1
  shrink = shrink1
instance Arbitrary1 Sequence.Seq where
  liftArbitrary = fmap Sequence.fromList . liftArbitrary
  liftShrink shr = map Sequence.fromList . liftShrink shr . toList
instance Arbitrary a => Arbitrary (Sequence.Seq a) where
  arbitrary = arbitrary1
  shrink = shrink1
instance Arbitrary1 Tree.Tree where
    liftArbitrary arb = sized $ \n -> do
        k <- chooseInt (0, n)
        go k
      where
        go n = do -- n is the size of the trees.
            value <- arb
            pars <- arbPartition (n - 1) -- can go negative!
            forest <- mapM go pars
            return $ Tree.Node value forest

        arbPartition :: Int -> Gen [Int]
        arbPartition k = case compare k 1 of
            LT -> pure []
            EQ -> pure [1]
            GT -> do
                first <- chooseInt (1, k)
                rest <- arbPartition $ k - first
                shuffle (first : rest)

    liftShrink shr = go
      where
        go (Tree.Node val forest) = forest ++
            [ Tree.Node e fs
            | (e, fs) <- liftShrink2 shr (liftShrink go) (val, forest)
            ]
instance Arbitrary a => Arbitrary (Tree.Tree a) where
  arbitrary = arbitrary1
  shrink = shrink1

-- Arbitrary instance for Ziplist
instance Arbitrary1 ZipList where
  liftArbitrary = fmap ZipList . liftArbitrary
  liftShrink shr = map ZipList . liftShrink shr . getZipList
instance Arbitrary a => Arbitrary (ZipList a) where
  arbitrary = arbitrary1
  shrink = shrink1

#ifndef NO_TRANSFORMERS
-- Arbitrary instance for transformers' Functors
instance Arbitrary1 Identity where
  liftArbitrary = fmap Identity
  liftShrink shr = map Identity . shr . runIdentity
instance Arbitrary a => Arbitrary (Identity a) where
  arbitrary = arbitrary1
  shrink = shrink1

instance Arbitrary2 Constant where
  liftArbitrary2 arbA _ = fmap Constant arbA
  liftShrink2 shrA _ = fmap Constant . shrA . getConstant
instance Arbitrary a => Arbitrary1 (Constant a) where
  liftArbitrary = liftArbitrary2 arbitrary
  liftShrink = liftShrink2 shrink
-- Have to be defined explicitly, as Constant is kind polymorphic
instance Arbitrary a => Arbitrary (Constant a b) where
  arbitrary = fmap Constant arbitrary
  shrink = map Constant . shrink . getConstant

instance (Arbitrary1 f, Arbitrary1 g) => Arbitrary1 (Product f g) where
  liftArbitrary arb = liftM2 Pair (liftArbitrary arb) (liftArbitrary arb)
  liftShrink shr (Pair f g) =
    [ Pair f' g | f' <- liftShrink shr f ] ++
    [ Pair f g' | g' <- liftShrink shr g ]
instance (Arbitrary1 f, Arbitrary1 g, Arbitrary a) => Arbitrary (Product f g a) where
  arbitrary = arbitrary1
  shrink = shrink1

instance (Arbitrary1 f, Arbitrary1 g) => Arbitrary1 (Compose f g) where
  liftArbitrary = fmap Compose . liftArbitrary . liftArbitrary
  liftShrink shr = map Compose . liftShrink (liftShrink shr) . getCompose
instance (Arbitrary1 f, Arbitrary1 g, Arbitrary a) => Arbitrary (Compose f g a) where
  arbitrary = arbitrary1
  shrink = shrink1
#endif

-- Arbitrary instance for Const
instance Arbitrary2 Const where
  liftArbitrary2 arbA _ = fmap Const arbA
  liftShrink2 shrA _ = fmap Const . shrA . getConst
instance Arbitrary a => Arbitrary1 (Const a) where
  liftArbitrary = liftArbitrary2 arbitrary
  liftShrink = liftShrink2 shrink
-- Have to be defined explicitly, as Const is kind polymorphic
instance Arbitrary a => Arbitrary (Const a b) where
  arbitrary = fmap Const arbitrary
  shrink = map Const . shrink . getConst

instance Arbitrary (m a) => Arbitrary (WrappedMonad m a) where
  arbitrary = WrapMonad <$> arbitrary
  shrink (WrapMonad a) = map WrapMonad (shrink a)

instance Arbitrary (a b c) => Arbitrary (WrappedArrow a b c) where
  arbitrary = WrapArrow <$> arbitrary
  shrink (WrapArrow a) = map WrapArrow (shrink a)

-- Arbitrary instances for Monoid
instance Arbitrary a => Arbitrary (Monoid.Dual a) where
  arbitrary = fmap Monoid.Dual arbitrary
  shrink = map Monoid.Dual . shrink . Monoid.getDual

instance (Arbitrary a, CoArbitrary a) => Arbitrary (Monoid.Endo a) where
  arbitrary = fmap Monoid.Endo arbitrary
  shrink = map Monoid.Endo . shrink . Monoid.appEndo

instance Arbitrary Monoid.All where
  arbitrary = fmap Monoid.All arbitrary
  shrink = map Monoid.All . shrink . Monoid.getAll

instance Arbitrary Monoid.Any where
  arbitrary = fmap Monoid.Any arbitrary
  shrink = map Monoid.Any . shrink . Monoid.getAny

instance Arbitrary a => Arbitrary (Monoid.Sum a) where
  arbitrary = fmap Monoid.Sum arbitrary
  shrink = map Monoid.Sum . shrink . Monoid.getSum

instance Arbitrary a => Arbitrary (Monoid.Product a) where
  arbitrary = fmap Monoid.Product  arbitrary
  shrink = map Monoid.Product  . shrink . Monoid.getProduct

#if defined(MIN_VERSION_base)
#if MIN_VERSION_base(3,0,0)
instance Arbitrary a => Arbitrary (Monoid.First a) where
  arbitrary = fmap Monoid.First arbitrary
  shrink = map Monoid.First . shrink . Monoid.getFirst

instance Arbitrary a => Arbitrary (Monoid.Last a) where
  arbitrary = fmap Monoid.Last arbitrary
  shrink = map Monoid.Last . shrink . Monoid.getLast
#endif

#if MIN_VERSION_base(4,8,0)
instance Arbitrary (f a) => Arbitrary (Monoid.Alt f a) where
  arbitrary = fmap Monoid.Alt arbitrary
  shrink = map Monoid.Alt . shrink . Monoid.getAlt
#endif
#endif

-- | Generates 'Version' with non-empty non-negative @versionBranch@, and empty @versionTags@
instance Arbitrary Version where
  arbitrary = sized $ \n ->
    do k <- chooseInt (0, log2 n)
       xs <- vectorOf (k+1) arbitrarySizedNatural
       return (Version xs [])
    where
      log2 :: Int -> Int
      log2 n | n <= 1 = 0
             | otherwise = 1 + log2 (n `div` 2)

  shrink (Version xs _) =
    [ Version xs' []
    | xs' <- shrink xs
    , length xs' > 0
    , all (>=0) xs'
    ]

instance Arbitrary QCGen where
  arbitrary = MkGen (\g _ -> g)

instance Arbitrary ExitCode where
  arbitrary = frequency [(1, return ExitSuccess), (3, liftM ExitFailure arbitrary)]

  shrink (ExitFailure x) = ExitSuccess : [ ExitFailure x' | x' <- shrink x ]
  shrink _        = []

#if defined(MIN_VERSION_base)
#if MIN_VERSION_base(4,2,0)
instance Arbitrary Newline where
  arbitrary = elements [LF, CRLF]

  -- The behavior of code for LF is generally simpler than for CRLF
  -- See the documentation for this type, which states that Haskell
  -- Internally always assumes newlines are \n and this type represents
  -- how to translate that to and from the outside world, where LF means
  -- no translation.
  shrink LF = []
  shrink CRLF = [LF]

instance Arbitrary NewlineMode where
  arbitrary = NewlineMode <$> arbitrary <*> arbitrary

  shrink (NewlineMode inNL outNL) = [NewlineMode inNL' outNL' | (inNL', outNL') <- shrink (inNL, outNL)]
#endif
#endif

-- ** Helper functions for implementing arbitrary

-- | Apply a binary function to random arguments.
applyArbitrary2 :: (Arbitrary a, Arbitrary b) => (a -> b -> r) -> Gen r
applyArbitrary2 f = liftA2 f arbitrary arbitrary

-- | Apply a ternary function to random arguments.
applyArbitrary3
  :: (Arbitrary a, Arbitrary b, Arbitrary c)
  => (a -> b -> c -> r) -> Gen r
applyArbitrary3 f = liftA3 f arbitrary arbitrary arbitrary

-- | Apply a function of arity 4 to random arguments.
applyArbitrary4
  :: (Arbitrary a, Arbitrary b, Arbitrary c, Arbitrary d)
  => (a -> b -> c -> d -> r) -> Gen r
applyArbitrary4 f = applyArbitrary3 (uncurry f)

-- | Generates an integral number. The number can be positive or negative
-- and its maximum absolute value depends on the size parameter.
arbitrarySizedIntegral :: Integral a => Gen a
arbitrarySizedIntegral =
  sized $ \n ->
  inBounds fromIntegral (chooseInt (-n, n))

-- | Generates a natural number. The number's maximum value depends on
-- the size parameter.
arbitrarySizedNatural :: Integral a => Gen a
arbitrarySizedNatural =
  sized $ \n ->
  inBounds fromIntegral (chooseInt (0, n))

inBounds :: Integral a => (Int -> a) -> Gen Int -> Gen a
inBounds fi g = fmap fi (g `suchThat` (\x -> toInteger x == toInteger (fi x)))

-- | Uniformly generates a fractional number. The number can be positive or negative
-- and its maximum absolute value depends on the size parameter.
arbitrarySizedFractional :: Fractional a => Gen a
arbitrarySizedFractional =
  sized $ \n -> do
    denom <- chooseInt (1, max 1 n)
    numer <- chooseInt (-n*denom, n*denom)
    pure $ fromIntegral numer / fromIntegral denom

-- Useful for getting at minBound and maxBound without having to
-- fiddle around with asTypeOf.
{-# INLINE withBounds #-}
withBounds :: Bounded a => (a -> a -> Gen a) -> Gen a
withBounds k = k minBound maxBound

-- | Generates an integral number. The number is chosen uniformly from
-- the entire range of the type. You may want to use
-- 'arbitrarySizedBoundedIntegral' instead.
arbitraryBoundedIntegral :: (Bounded a, Integral a) => Gen a
arbitraryBoundedIntegral = chooseBoundedIntegral (minBound, maxBound)

-- | Generates an element of a bounded type. The element is
-- chosen from the entire range of the type.
arbitraryBoundedRandom :: (Bounded a, Random a) => Gen a
arbitraryBoundedRandom = choose (minBound,maxBound)

-- | Generates an element of a bounded enumeration.
arbitraryBoundedEnum :: (Bounded a, Enum a) => Gen a
arbitraryBoundedEnum = chooseEnum (minBound, maxBound)

-- | Generates an integral number from a bounded domain. The number is
-- chosen from the entire range of the type, but small numbers are
-- generated more often than big numbers. Inspired by demands from
-- Phil Wadler.
arbitrarySizedBoundedIntegral :: (Bounded a, Integral a) => Gen a
-- INLINEABLE so that this combinator gets specialised at each type,
-- which means that the constant 'bits' in the let-block below will
-- only be computed once.
{-# INLINEABLE arbitrarySizedBoundedIntegral #-}
arbitrarySizedBoundedIntegral =
  withBounds $ \mn mx ->
  let ilog2 1 = 0
      ilog2 n | n > 0 = 1 + ilog2 (n `div` 2)

      -- How many bits are needed to represent this type?
      -- (This number is an upper bound, not exact.)
      bits = ilog2 (toInteger mx - toInteger mn + 1) in
  sized $ \k ->
    let
      -- Reach maximum size by k=80, or quicker for small integer types
      power = ((bits `max` 40) * k) `div` 80

      -- Bounds should be 2^power, but:
      --   * clamp the result to minBound/maxBound
      --   * clamp power to 'bits', in case k is a huge number
      lo = toInteger mn `max` (-1 `shiftL` (power `min` bits))
      hi = toInteger mx `min` (1 `shiftL` (power `min` bits)) in
    fmap fromInteger (chooseInteger (lo, hi))

-- ** Generators for various kinds of character

-- | Generates any Unicode character (but not a surrogate)
arbitraryUnicodeChar :: Gen Char
arbitraryUnicodeChar =
  arbitraryBoundedEnum `suchThat` isValidUnicode
  where
    isValidUnicode c = case generalCategory c of
      Surrogate -> False
      NotAssigned -> False
      _ -> True

-- | Generates a random ASCII character (0-127).
arbitraryASCIIChar :: Gen Char
arbitraryASCIIChar = chooseEnum ('\0', '\127')

-- | Generates a printable Unicode character.
arbitraryPrintableChar :: Gen Char
arbitraryPrintableChar = arbitrary `suchThat` isPrint

-- ** Helper functions for implementing shrink

-- | Returns no shrinking alternatives.
shrinkNothing :: a -> [a]
shrinkNothing _ = []

-- | Map a shrink function to another domain. This is handy if your data type
-- has special invariants, but is /almost/ isomorphic to some other type.
--
-- @
-- shrinkOrderedList :: (Ord a, Arbitrary a) => [a] -> [[a]]
-- shrinkOrderedList = shrinkMap sort id
--
-- shrinkSet :: (Ord a, Arbitrary a) => Set a -> [Set a]
-- shrinkSet = shrinkMap fromList toList
-- @
shrinkMap :: Arbitrary a => (a -> b) -> (b -> a) -> b -> [b]
shrinkMap f g = shrinkMapBy f g shrink

-- | Non-overloaded version of `shrinkMap`.
shrinkMapBy :: (a -> b) -> (b -> a) -> (a -> [a]) -> b -> [b]
shrinkMapBy f g shr = map f . shr . g

-- | Shrink an integral number.
shrinkIntegral :: Integral a => a -> [a]
shrinkIntegral x =
  nub $
  [ -x
  | x < 0, -x > x
  ] ++
  [ x'
  | x' <- takeWhile (<< x) (0:[ x - i | i <- tail (iterate (`quot` 2) x) ])
  ]
 where
   -- a << b is "morally" abs a < abs b, but taking care of overflow.
   a << b = case (a >= 0, b >= 0) of
            (True,  True)  -> a < b
            (False, False) -> a > b
            (True,  False) -> a + b < 0
            (False, True)  -> a + b > 0

-- | Shrink a fraction, preferring numbers with smaller
-- numerators or denominators. See also 'shrinkDecimal'.
shrinkRealFrac :: RealFrac a => a -> [a]
shrinkRealFrac x
  | not (x == x)  = 0 : take 10 (iterate (*2) 0) -- NaN
  | not (2*x+1>x) = 0 : takeWhile (<x) (iterate (*2) 0) -- infinity
  | x < 0 = negate x:map negate (shrinkRealFrac (negate x))
  | otherwise =
    -- To ensure termination
    filter (\y -> abs y < abs x) $
      -- Try shrinking to an integer first
      map fromInteger (shrink (truncate x) ++ [truncate x]) ++
      -- Shrink the numerator
      [fromRational (num' % denom) | num' <- shrink num] ++
      -- Shrink the denominator, and keep the fraction as close
      -- to the original as possible, rounding towards zero
      [fromRational (truncate (num * denom' % denom) % denom')
      | denom' <- shrink denom, denom' /= 0 ]
  where
    num = numerator (toRational x)
    denom = denominator (toRational x)

-- | Shrink a real number, preferring numbers with shorter
-- decimal representations. See also 'shrinkRealFrac'.
shrinkDecimal :: RealFrac a => a -> [a]
shrinkDecimal x
  | not (x == x)  = 0 : take 10 (iterate (*2) 0)        -- NaN
  | not (2*abs x+1>abs x) = 0 : takeWhile (<x) (iterate (*2) 0) -- infinity
  | otherwise =
    -- e.g. shrink pi =
    --   shrink 3 ++ map (/ 10) (shrink 31) ++
    --   map (/ 100) (shrink 314) + ...,
    -- where the inner calls to shrink use integer shrinking.
    [ y
    | precision <- take 6 (iterate (*10) 1),
      let m = round (toRational x * precision),
      precision == 1 || m `mod` 10 /= 0, -- don't allow shrinking to increase digits
      n <- m:shrink m,
      let y = fromRational (fromInteger n / precision),
      abs y < abs x ]

--------------------------------------------------------------------------
-- ** CoArbitrary

#ifndef NO_GENERICS
-- | Used for random generation of functions.
-- You should consider using 'Test.QuickCheck.Fun' instead, which
-- can show the generated functions as strings.
--
-- If you are using a recent GHC, there is a default definition of
-- 'coarbitrary' using 'genericCoarbitrary', so if your type has a
-- 'Generic' instance it's enough to say
--
-- > instance CoArbitrary MyType
--
-- You should only use 'genericCoarbitrary' for data types where
-- equality is structural, i.e. if you can't have two different
-- representations of the same value. An example where it's not
-- safe is sets implemented using binary search trees: the same
-- set can be represented as several different trees.
-- Here you would have to explicitly define
-- @coarbitrary s = coarbitrary (toList s)@.
#else
-- | Used for random generation of functions.
#endif
class CoArbitrary a where
  -- | Used to generate a function of type @a -> b@.
  -- The first argument is a value, the second a generator.
  -- You should use 'variant' to perturb the random generator;
  -- the goal is that different values for the first argument will
  -- lead to different calls to 'variant'. An example will help:
  --
  -- @
  -- instance CoArbitrary a => CoArbitrary [a] where
  --   coarbitrary []     = 'variant' 0
  --   coarbitrary (x:xs) = 'variant' 1 . coarbitrary (x,xs)
  -- @
  coarbitrary :: a -> Gen b -> Gen b
#ifndef NO_GENERICS
  default coarbitrary :: (Generic a, GCoArbitrary (Rep a)) => a -> Gen b -> Gen b
  coarbitrary = genericCoarbitrary

-- | Generic CoArbitrary implementation.
genericCoarbitrary :: (Generic a, GCoArbitrary (Rep a)) => a -> Gen b -> Gen b
genericCoarbitrary = gCoarbitrary . from

class GCoArbitrary f where
  gCoarbitrary :: f a -> Gen b -> Gen b

instance GCoArbitrary U1 where
  gCoarbitrary U1 = id

instance (GCoArbitrary f, GCoArbitrary g) => GCoArbitrary (f :*: g) where
  -- Like the instance for tuples.
  gCoarbitrary (l :*: r) = gCoarbitrary l . gCoarbitrary r

instance (GCoArbitrary f, GCoArbitrary g) => GCoArbitrary (f :+: g) where
  -- Like the instance for Either.
  gCoarbitrary (L1 x) = variant 0 . gCoarbitrary x
  gCoarbitrary (R1 x) = variant 1 . gCoarbitrary x

instance GCoArbitrary f => GCoArbitrary (M1 i c f) where
  gCoarbitrary (M1 x) = gCoarbitrary x

instance CoArbitrary a => GCoArbitrary (K1 i a) where
  gCoarbitrary (K1 x) = coarbitrary x
#endif

{-# DEPRECATED (><) "Use ordinary function composition instead" #-}
-- | Combine two generator perturbing functions, for example the
-- results of calls to 'variant' or 'coarbitrary'.
(><) :: (Gen a -> Gen a) -> (Gen a -> Gen a) -> (Gen a -> Gen a)
(><) = (.)

instance (Arbitrary a, CoArbitrary b) => CoArbitrary (a -> b) where
  coarbitrary f gen =
    do xs <- arbitrary
       coarbitrary (map f xs) gen

instance CoArbitrary () where
  coarbitrary _ = id

instance CoArbitrary Bool where
  coarbitrary False = variant 0
  coarbitrary True  = variant 1

instance CoArbitrary Ordering where
  coarbitrary GT = variant 0
  coarbitrary EQ = variant 1
  coarbitrary LT = variant 2

instance CoArbitrary a => CoArbitrary (Maybe a) where
  coarbitrary Nothing  = variant 0
  coarbitrary (Just x) = variant 1 . coarbitrary x

instance (CoArbitrary a, CoArbitrary b) => CoArbitrary (Either a b) where
  coarbitrary (Left x)  = variant 0 . coarbitrary x
  coarbitrary (Right y) = variant 1 . coarbitrary y

instance CoArbitrary a => CoArbitrary [a] where
  coarbitrary []     = variant 0
  coarbitrary (x:xs) = variant 1 . coarbitrary (x,xs)

instance (Integral a, CoArbitrary a) => CoArbitrary (Ratio a) where
  coarbitrary r = coarbitrary (numerator r,denominator r)

#ifndef NO_FIXED
instance HasResolution a => CoArbitrary (Fixed a) where
  coarbitrary = coarbitraryReal
#endif

#if defined(MIN_VERSION_base) && MIN_VERSION_base(4,4,0)
instance CoArbitrary a => CoArbitrary (Complex a) where
#else
instance (RealFloat a, CoArbitrary a) => CoArbitrary (Complex a) where
#endif
  coarbitrary (x :+ y) = coarbitrary x . coarbitrary y

instance (CoArbitrary a, CoArbitrary b)
      => CoArbitrary (a,b)
 where
  coarbitrary (x,y) = coarbitrary x
                    . coarbitrary y

instance (CoArbitrary a, CoArbitrary b, CoArbitrary c)
      => CoArbitrary (a,b,c)
 where
  coarbitrary (x,y,z) = coarbitrary x
                      . coarbitrary y
                      . coarbitrary z

instance (CoArbitrary a, CoArbitrary b, CoArbitrary c, CoArbitrary d)
      => CoArbitrary (a,b,c,d)
 where
  coarbitrary (x,y,z,v) = coarbitrary x
                        . coarbitrary y
                        . coarbitrary z
                        . coarbitrary v

instance (CoArbitrary a, CoArbitrary b, CoArbitrary c, CoArbitrary d, CoArbitrary e)
      => CoArbitrary (a,b,c,d,e)
 where
  coarbitrary (x,y,z,v,w) = coarbitrary x
                          . coarbitrary y
                          . coarbitrary z
                          . coarbitrary v
                          . coarbitrary w

-- typical instance for primitive (numerical) types

instance CoArbitrary Integer where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Int where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Int8 where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Int16 where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Int32 where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Int64 where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Word where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Word8 where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Word16 where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Word32 where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Word64 where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Char where
  coarbitrary = coarbitrary . ord

instance CoArbitrary Float where
  coarbitrary = coarbitraryReal

instance CoArbitrary Double where
  coarbitrary = coarbitraryReal

-- Coarbitrary instances for container types
instance CoArbitrary a => CoArbitrary (Set.Set a) where
  coarbitrary = coarbitrary. Set.toList
instance (CoArbitrary k, CoArbitrary v) => CoArbitrary (Map.Map k v) where
  coarbitrary = coarbitrary . Map.toList
instance CoArbitrary IntSet.IntSet where
  coarbitrary = coarbitrary . IntSet.toList
instance CoArbitrary a => CoArbitrary (IntMap.IntMap a) where
  coarbitrary = coarbitrary . IntMap.toList
instance CoArbitrary a => CoArbitrary (Sequence.Seq a) where
  coarbitrary = coarbitrary . toList
instance CoArbitrary a => CoArbitrary (Tree.Tree a) where
  coarbitrary (Tree.Node val forest) = coarbitrary val . coarbitrary forest

-- CoArbitrary instance for Ziplist
instance CoArbitrary a => CoArbitrary (ZipList a) where
  coarbitrary = coarbitrary . getZipList

#ifndef NO_TRANSFORMERS
-- CoArbitrary instance for transformers' Functors
instance CoArbitrary a => CoArbitrary (Identity a) where
  coarbitrary = coarbitrary . runIdentity

instance CoArbitrary a => CoArbitrary (Constant a b) where
  coarbitrary = coarbitrary . getConstant
#endif

-- CoArbitrary instance for Const
instance CoArbitrary a => CoArbitrary (Const a b) where
  coarbitrary = coarbitrary . getConst

-- CoArbitrary instances for Monoid
instance CoArbitrary a => CoArbitrary (Monoid.Dual a) where
  coarbitrary = coarbitrary . Monoid.getDual

instance (Arbitrary a, CoArbitrary a) => CoArbitrary (Monoid.Endo a) where
  coarbitrary = coarbitrary . Monoid.appEndo

instance CoArbitrary Monoid.All where
  coarbitrary = coarbitrary . Monoid.getAll

instance CoArbitrary Monoid.Any where
  coarbitrary = coarbitrary . Monoid.getAny

instance CoArbitrary a => CoArbitrary (Monoid.Sum a) where
  coarbitrary = coarbitrary . Monoid.getSum

instance CoArbitrary a => CoArbitrary (Monoid.Product a) where
  coarbitrary = coarbitrary . Monoid.getProduct

#if defined(MIN_VERSION_base)
#if MIN_VERSION_base(3,0,0)
instance CoArbitrary a => CoArbitrary (Monoid.First a) where
  coarbitrary = coarbitrary . Monoid.getFirst

instance CoArbitrary a => CoArbitrary (Monoid.Last a) where
  coarbitrary = coarbitrary . Monoid.getLast
#endif

#if MIN_VERSION_base(4,8,0)
instance CoArbitrary (f a) => CoArbitrary (Monoid.Alt f a) where
  coarbitrary = coarbitrary . Monoid.getAlt
#endif
#endif

instance CoArbitrary Version where
  coarbitrary (Version a b) = coarbitrary (a, b)

#if defined(MIN_VERSION_base)
#if MIN_VERSION_base(4,2,0)
instance CoArbitrary Newline where
  coarbitrary LF = variant 0
  coarbitrary CRLF = variant 1

instance CoArbitrary NewlineMode where
  coarbitrary (NewlineMode inNL outNL) = coarbitrary inNL . coarbitrary outNL
#endif
#endif

-- ** Helpers for implementing coarbitrary

-- | A 'coarbitrary' implementation for integral numbers.
coarbitraryIntegral :: Integral a => a -> Gen b -> Gen b
coarbitraryIntegral = variant

-- | A 'coarbitrary' implementation for real numbers.
coarbitraryReal :: Real a => a -> Gen b -> Gen b
coarbitraryReal x = coarbitrary (toRational x)

-- | 'coarbitrary' helper for lazy people :-).
coarbitraryShow :: Show a => a -> Gen b -> Gen b
coarbitraryShow x = coarbitrary (show x)

-- | A 'coarbitrary' implementation for enums.
coarbitraryEnum :: Enum a => a -> Gen b -> Gen b
coarbitraryEnum = variant . fromEnum

--------------------------------------------------------------------------
-- ** arbitrary generators

-- these are here and not in Gen because of the Arbitrary class constraint

-- | Generates a list of a given length.
vector :: Arbitrary a => Int -> Gen [a]
vector k = vectorOf k arbitrary

-- | Generates an ordered list.
orderedList :: (Ord a, Arbitrary a) => Gen [a]
orderedList = sort `fmap` arbitrary

-- | Generates an infinite list.
infiniteList :: Arbitrary a => Gen [a]
infiniteList = infiniteListOf arbitrary

--------------------------------------------------------------------------
-- ** Rational helper

infixr 5 :<
data Stream a = !a :< Stream a

streamNth :: Int -> Stream a -> a
streamNth n (x :< xs) | n <= 0    = x
                      | otherwise = streamNth (n - 1) xs

-- We read into this stream only with ~size argument,
-- so it's ok to have it as CAF.
--
rationalUniverse :: Stream Rational
rationalUniverse = 0 :< 1 :< (-1) :< go leftSideStream
  where
    go (x :< xs) =
      let nx = -x
          rx = recip x
          nrx = -rx
      in nx `seq` rx `seq` nrx `seq` (x :< rx :< nx :< nrx :< go xs)

-- All the rational numbers on the left side of the Calkin-Wilf tree,
-- in breadth-first order.
leftSideStream :: Stream Rational
leftSideStream = (1 % 2) :< go leftSideStream
  where
    go (x :< xs) =
        lChild `seq` rChild `seq`
        (lChild :< rChild :< go xs)
      where
        nd = numerator x + denominator x
        lChild = numerator x % nd
        rChild = nd % denominator x

--------------------------------------------------------------------------
-- the end.


text:
{-# LANGUAGE GADTs,RebindableSyntax,CPP,FlexibleContexts,FlexibleInstances,ConstraintKinds #-}
{-# LANGUAGE StandaloneDeriving,DeriveDataTypeable #-}
{-# OPTIONS_GHC -dcore-lint #-}
{-
 - This test suite ensures that the rewrites that HerbiePlugin performs
 - give the correct results.
 -}

module Main
    where

import SubHask

import System.IO
-- import Data.Complex
-- import Linear.Quaternion
-- import Linear.V3
-- import Linear.Vector

--------------------------------------------------------------------------------

test1a :: Double -> Double -> Double
test1a far near = -(2 * far * near) / (far - near)

{-# ANN test1b "NoHerbie" #-}
test1b :: Double -> Double -> Double
test1b far near = -(2 * far * near) / (far - near)

{-# ANN test1c "NoHerbie" #-}
test1c :: Double -> Double -> Double
test1c far near = if far < -1.7210442634149447e81
    then ((-2 * far) / (far - near)) * near
    else if far < 8.364504563556443e16
        then -2 * far * (near / (far - near))
        else ((-2 * far) / (far - near)) * near

{-
--------------------

test2a :: Double -> Double -> Double
test2a a b = a + ((b - a) / 2)

{-# ANN test2b "NoHerbie" #-}
test2b :: Double -> Double -> Double
test2b a b = a + ((b - a) / 2)

--------------------

-- test3a :: Quaternion Double -> Quaternion Double -> Quaternion Double
-- test3a (Quaternion q0 (V3 q1 q2 q3)) (Quaternion r0 (V3 r1 r2 r3)) =
--     Quaternion (r0*q0+r1*q1+r2*q2+r3*q3)
--                (V3 (r0*q1-r1*q0-r2*q3+r3*q2)
--                    (r0*q2+r1*q3-r2*q0-r3*q1)
--                    (r0*q3-r1*q2+r2*q1-r3*q0))
--                ^/ (r0*r0 + r1*r1 + r2*r2 + r3*r3)
--
-- {-# ANN test3b "NoHerbie" #-}
-- test3b :: Quaternion Double -> Quaternion Double -> Quaternion Double
-- test3b (Quaternion q0 (V3 q1 q2 q3)) (Quaternion r0 (V3 r1 r2 r3)) =
--     Quaternion (r0*q0+r1*q1+r2*q2+r3*q3)
--                (V3 (r0*q1-r1*q0-r2*q3+r3*q2)
--                    (r0*q2+r1*q3-r2*q0-r3*q1)
--                    (r0*q3-r1*q2+r2*q1-r3*q0))
--                ^/ (r0*r0 + r1*r1 + r2*r2 + r3*r3)

--------------------

data Yo a = Yo
    { yo_x2y :: a
    , yo_y2x :: a
    }

test4 :: Real a => a -> a -> Yo a
test4 x y = Yo
    { yo_x2y = x * x * y
    , yo_y2x = y * y * x
    }

test5 :: Float -> Float -> Float
test5 x y = (x * x) + (2 * x * y) + (y * y)

--------------------------------------------------------------------------------

-- asinh_ :: Complex Double -> Complex Double
-- asinh_ x = log (x + sqrt (1.0+x*x))
--
-- acosh_ :: Complex Double -> Complex Double
-- acosh_ x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))

atanh_ :: Double -> Double
atanh_ x = 0.5 * log ((1.0+x) / (1.0-x))

--------------------------------------------------------------------------------
-}

#define mkTest(f1,f2,a,b) \
    putStrLn $ "mkTest: " ++ show (f1 (a) (b)); \
    putStrLn $ "mkTest: " ++ show (f2 (a) (b)); \
    putStrLn ""

#define mkTestB(f1,f2,a) \
    putStrLn $ "mkTest: " ++ show (f1 (a)); \
    putStrLn $ "mkTest: " ++ show (f2 (a)); \
    putStrLn ""

main = do
    mkTest(test1a,test1b,-2e90,6)
    mkTest(test1a,test1b,3,4)
    mkTest(test1a,test1b,2e90,6)

    mkTest(test1a,test1c,-2e90,6)
    mkTest(test1a,test1c,3,4)
    mkTest(test1a,test1c,2e90,6)

{-
    mkTest(test2a,test2b,1,2)

--     mkTest(test3a,test3b,(Quaternion 1 (V3 1 2 3)),(Quaternion 2 (V3 2 3 4)))

--     mkTestB(asinh,asinh_,5e-17::Complex Double)
--     mkTestB(acosh,acosh_,5e-17::Complex Double)
    mkTestB(atanh,atanh_,5e-17::Double)
-}

    putStrLn "done"

ext: hs

max_stars_repo_name: mikeizbicki/herbie-haskell

max_issues_repo_name: mikeizbicki/herbie-haskell

max_forks_repo_name: mikeizbicki/herbie-haskell


text:
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE OverlappingInstances #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}

module Numeric.LinearAlgebra.Dimensional.DK.Internal
{-
(
   -- * Data.Packed.Vector
   (@>),
   -- * Data.Packed.Matrix
   -- ** dimension
   cols, rows,
   colsNT, rowsNT,
   hasRows, hasCols,
   -- (><),
   trans,
   -- reshape, flatten, fromLists, toLists, buildMatrix,
   (@@>),

   -- asRow, asColumn, fromRows, toRows, fromColumns, toColumns
   -- fromBlocks
   -- diagBlock,
   -- toBlocks, toBlocksEvery, repmat, flipud, fliprl
   -- subMatrix, takeRows, dropRows, takeColumns, dropColumns,
   -- extractRows, diagRect, takeDiag, mapMatrix,
   -- mapMatrixWithIndexM, mapMatrixWithIndexM_, liftMatrix,
   -- liftMatrix2, liftMatrix2Auto, fromArray2D,

   ident, -- where to put this?
   -- * Numeric.Container
   -- constant, linspace,
   diag,
   ctrans,
   -- ** Container class
   scalar,
   conj,
   scale, scaleRecip,
   recipMat,
   add,
   sub,
   mul,
   divide,
   equal,
   arctan2,
   hconcat,
   vconcat,
   concat,
   konst,
   zeroes,
   -- build, atIndex, minIndex, maxIndex, minElement, maxElement,
   -- sumElements, prodElements, step, cond, find, assoc, accum,
   -- Convert
   -- ** Product class
   multiply,
   -- dot, absSum, norm1, norm2, normInf,
   -- norm1, normInf,
   -- optimiseMult, mXm, mXv, vXm, (<.>),
   -- (<>), (<\>), outer, kronecker,
   (<>),
   -- ** Random numbers
   -- ** Element conversion
   -- ** Input/Output
   -- ** Experimental

   -- * Numeric.LinearAlgebra.Algorithms
   -- | incomplete wrapper for "Numeric.LinearAlgebra.Algorithms"

   -- ** Linear Systems
   -- linearSolve, luSolve, cholSolve, linearSolveLS, linearSolveSVD,
   inv,
   PInv(pinv),
   pinvTol,
   det,
   -- invlndet,
   rank,
   -- rcond,
   -- ** Matrix factorizations

   -- *** Singular value decomposition
   -- *** Eigensystems
   -- $eigs
   --wrapEig, wrapEigOnly,
   --EigCxt,
   -- **** eigenvalues and eigenvectors
   --eig,
   --eigC,
   --eigH,
   --eigH',
   --eigR,
   --eigS,
   --eigS',
   --eigSH,
   --eigSH',

   -- **** eigenvalues
   --eigOnlyC,
   --eigOnlyH,
   --eigOnlyR,
   --eigOnlyS,
   --eigenvalues,
   --eigenvaluesSH,
   --eigenvaluesSH',

   -- *** QR
   -- *** Cholesky
   -- *** Hessenberg
   -- *** Schur
   -- *** LU

   -- ** Matrix functions
   -- sqrtm, matFunc
   expm,

   -- ** Nullspace
   -- ** Norms
   -- ** Misc
   -- ** Util

   -- * actually internal
   toDM,
   DimMat(..),
  )
  -}
   where
import Foreign.Storable (Storable)
import GHC.Exts (Constraint)
import Numeric.Units.Dimensional.DK.Prelude hiding (concat)
import qualified Prelude as P
import qualified Numeric.NumType.DK.Integers as N

import Numeric.LinearAlgebra.Dimensional.DK.Shapes as S
import Data.Proxy

import qualified Numeric.Matrix as M
import Numeric.Matrix ((<|>), (<->))

data DimMat (shape :: Shape) a where
  DimMat :: M.Matrix a -> DimMat ('MatrixShape g rs cs) a
  DimVec :: M.Matrix a -> DimMat ('VectorShape d ds) a -- as a column vector

type ValidElement = M.MatrixElement

deriving instance (Show a, ValidElement a) => Show (DimMat s a)

(@>) :: (NN.KnownNat n, KnownDimension (VectorElement ('VectorShape d ds) n), Fractional a, ValidElement a)
     => DimMat ('VectorShape d ds) a
     -> Proxy n
     -> Quantity (VectorElement ('VectorShape d ds) n) a
(DimVec v) @> n = (v `M.at` (asIntFrom1 n,1)) *~ siUnit

(@@>) :: (NN.KnownNat nr, NN.KnownNat nc, KnownDimension (MatrixElement ('MatrixShape g rs cs) nr nc), Fractional a, ValidElement a)
    => DimMat ('MatrixShape g rs cs) a
    -> (Proxy nr, Proxy nc)
    -> Quantity (MatrixElement ('MatrixShape g rs cs) nr nc) a
DimMat m @@> (nr,nc) = (m `M.at` (asIntFrom1 nr, asIntFrom1 nc)) *~ siUnit

asInt, asIntFrom1 :: (NN.KnownNat n) => Proxy n -> Int
asInt = fromInteger . NN.natVal
asIntFrom1 = (P.+ 1) . asInt

{-
norm1 :: (sh ~ [r11 ': rs,ci], rs ~ MapConst r11 rs, ci ~ MapConst DOne ci, a ~ H.RealOf a)
         => DimMat sh a
         -> Quantity r11 a
norm1 (DimMat a) = Dimensional (H.pnorm H.PNorm1 a)

normInf :: (sh ~ [r11 ': rs,ci], rs ~ MapConst r11 rs, ci ~ MapConst DOne ci, a ~ H.RealOf a)
           => DimMat sh a
           -> Quantity r11 a
normInf (DimMat a) = Dimensional (H.pnorm H.Infinity a)
-}

{- | does H.'H.mXm' and H.'H.mXv'.

vXm and vXv (called dot) might be supported in the future too
-}
multiply :: (HasProduct s1 s2, ValidElement a)
    => DimMat s1 a -> DimMat s2 a
    -> DimMat (ShapeProduct s1 s2) a
multiply (DimMat a) (DimMat b) = DimMat (M.times a b)
multiply (DimMat a) (DimVec b) = DimVec (M.times a b)

infixl 7 <>
(<>) :: (HasProduct s1 s2, ValidElement a)
    => DimMat s1 a -> DimMat s2 a
    -> DimMat (ShapeProduct s1 s2) a
(<>) = multiply

trans :: (ValidElement a) => DimMat s a -> DimMat (ShapeTranspose s) a
trans (DimMat m) = DimMat $ M.transpose m

{-
pinvTol :: (PInv sh sh',
            a ~ Double,
           sh' ~ [ri2 ': _1 , DOne ': ci2]) => Double -> DimMat sh a -> DimMat sh' a
pinvTol tol (DimMat a) = DimMat (H.pinvTol tol a)

-}

det :: (Square s, KnownDimension (ShapeDeterminant s), Fractional a, ValidElement a) => DimMat s a -> Quantity (ShapeDeterminant s) a
det (DimMat m) = M.det m *~ siUnit

expm :: (s ~ ShapeProduct s s, HasProduct s s) => DimMat s a -> DimMat s a
expm = undefined

scale :: (Fractional a, ValidElement a, KnownDimension d)
         => Quantity d a -> DimMat s a -> DimMat (ShapeScale d s) a
scale x (DimMat m) = DimMat $ M.scale m (x /~ siUnit)
scale x (DimVec v) = DimVec $ M.scale v (x /~ siUnit)

add :: (ValidElement a) => DimMat s a -> DimMat s a -> DimMat s a
add (DimMat x) (DimMat y) = DimMat (M.plus x y)

sub :: (ValidElement a) => DimMat s a -> DimMat s a -> DimMat s a
sub (DimMat x) (DimMat y) = DimMat (M.minus x y)

equal :: (Eq a, ValidElement a) => DimMat s a -> DimMat s a -> Bool
equal (DimMat m1) (DimMat m2) = m1 == m2
equal (DimVec v1) (DimVec v2) = v1 == v2

hconcat :: (HorizontallyConcatenable s1 s2, ValidElement a) => DimMat s1 a -> DimMat s2 a -> DimMat (HorizontalConcatenation s1 s2) a
hconcat (DimMat m1) (DimMat m2) = DimMat (m1 <|> m2)
hconcat (DimMat m1) (DimVec v2) = DimMat (m1 <|> v2)
hconcat (DimVec v1) (DimMat m2) = DimMat (v1 <|> m2)
hconcat (DimVec v1) (DimVec v2) = DimMat (v1 <|> v2)

vconcat :: (VerticallyConcatenable s1 s2, ValidElement a) => DimMat s1 a -> DimMat s2 a -> DimMat (VerticalConcatenation s1 s2) a
vconcat (DimMat m1) (DimMat m2) = DimMat (m1 <-> m2)
vconcat (DimMat m1) (DimVec v2) = undefined --DimMat (m1 <-> M.transpose v2)
vconcat (DimVec v1) (DimMat m2) = undefined --DimMat (M.transpose v1 <-> m2)
vconcat (DimVec v1) (DimVec v2) = DimMat (M.transpose v1 <-> M.transpose v2)

fromRowVector :: (ValidElement a) => DimMat (VectorShape d ds) a -> DimMat (MatrixShape d '[] (MapDiv d ds)) a
fromRowVector (DimVec v) = DimMat (M.transpose v)

vconcat'  :: (ValidElement a,
              MapMulEq ds d cs
              )
             => DimMat (VectorShape d ds) a
             -> DimMat (MatrixShape g rs cs) a
             -> DimMat (MatrixShape d ((g/d) ': (MapMul (g/d) rs)) cs) a
vconcat' (DimVec v1) (DimMat m2) = DimMat (M.transpose v1 <-> m2)

-- I can't figure out why this is needed. But I also can't make a version of vconcat' that seems to work when used with fromRowVector. I don't get it.
vconcat'' :: (ValidElement a,
              MapMulEq ds2 (d2/d1) ds1
             )
             => DimMat (VectorShape d1 ds1) a
             -> DimMat (VectorShape d2 ds2) a
             -> DimMat (MatrixShape d1 '[d2/d1] (MapDiv d1 ds1)) a
vconcat'' (DimVec v1) (DimVec v2) = DimMat (M.transpose v1 <-> M.transpose v2)

concat :: (ValidElement a) => DimMat s1 a -> DimMat s2 a -> DimMat (VectorConcatenation s1 s2) a
concat (DimVec v1) (DimVec v2) = DimVec (v1 <-> v2)

--vecSingleton :: (KnownDimension d, Fractional a, ValidElement a) => Quantity d a -> DimMat (VectorShape d '[]) a
--vecSingleton x = DimVec . M.fromList $ [[x /~ siUnit]]

--vecCons :: (KnownDimension d', Fractional a, ValidElement a) => Quantity d' a -> DimMat (VectorShape d ds) a -> DimMat (VectorShape d' (d ': ds)) a
--vecCons x = concat (vecSingleton x)

vecSingleton :: (Fractional a, ValidElement a) => Quantity d a -> DimMat (VectorShape d '[]) a
vecSingleton x = DimVec . M.fromList $ [[12345]]

vecCons :: (Fractional a, ValidElement a) => Quantity d' a -> DimMat (VectorShape d ds) a -> DimMat (VectorShape d' (d ': ds)) a
vecCons x = concat (vecSingleton x)

rank :: DimMat s a -> Integer
rank = undefined

rows :: forall s a.(NN.KnownNat (ShapeRows s)) => DimMat s a -> Integer
rows _ = NN.natVal (Proxy :: Proxy (ShapeRows s))

cols :: forall s a.(NN.KnownNat (ShapeCols s)) => DimMat s a -> Integer
cols _ = NN.natVal (Proxy :: Proxy (ShapeCols s))

-- TODO: add a constraint that the row exists for better error message?
row :: forall n s a d ds.(MatrixRow s n ~ 'VectorShape d ds) => Proxy n -> DimMat s a -> DimMat ('VectorShape d ds) a
row _ _ = DimVec undefined

-- TODO: add a constraint that the column exists for better error message?
col :: forall n s a d ds.(MatrixColumn s n ~ 'VectorShape d ds) => Proxy n -> DimMat s a -> DimMat ('VectorShape d ds) a
col _ _ = DimVec undefined

{-
scalar :: (H.Field a,
          sh ~ ['[u], '[DOne]]) => Quantity u a -> DimMat sh a
scalar (Dimensional a) = DimMat (H.scalar a)
-}

{- | Numeric.Container.'H.konst', but the size is determined by the type.

>>> let n = hSucc (hSucc hZero) -- 2
>>> konst ((1::Double) *~ second) `hasRows` n `hasCols` n
2><2 1   1
s    1.0 1.0
s    1.0 1.0

-}
{-
konst :: forall u us ones a _1.
    (H.Field a,
     HNat2Integral (HLength ones),
     HNat2Integral (HLength us),
     ones ~ (DOne ': _1),
     AllEq DOne _1,
     AllEq u us)
    => Quantity u a -> DimMat [us, ones] a
konst (Dimensional a) = DimMat (H.konst a
    (hNat2Integral (proxy :: Proxy (HLength us)),
     hNat2Integral (proxy :: Proxy (HLength ones))))

-}

-- | identity matrix. The size is determined by the type.
ident :: forall g rs cs s a.(s ~ 'MatrixShape g rs cs, HasIdentity s, NN.KnownNat (ShapeRows s), ValidElement a) => DimMat ('MatrixShape g rs cs) a
ident = DimMat $ M.unit $ asInt (Proxy :: Proxy (ShapeRows s))

-- | zero matrix. The size and dimension is determined by the type.
zeroes :: forall g rs cs s a.(s ~ 'MatrixShape g rs cs, NN.KnownNat (ShapeRows s), NN.KnownNat (ShapeCols s), ValidElement a) => DimMat ('MatrixShape g rs cs) a
zeroes = DimMat $ M.matrix (r,c) (\(i,j) -> 0)
           where
             r = asInt (Proxy :: Proxy (ShapeRows s))
             c = asInt (Proxy :: Proxy (ShapeCols s))

trace :: (HasTrace s, Fractional a, ValidElement a, KnownDimension (ShapeTrace s)) => DimMat s a -> Quantity (ShapeTrace s) a
trace (DimMat m) = (P.sum $ M.trace m) *~ siUnit

conj :: DimMat s a -> DimMat s a
conj = undefined

-- | conjugate transpose
ctrans :: DimMat s a -> DimMat (ShapeTranspose s) a
ctrans = undefined

{-

diag :: (MapConst DOne v ~ c,
        c ~ (DOne ': _1)
        ) => DimMat '[v] t -> DimMat '[v,c] t
diag (DimVec a) = DimMat (H.diag a)
-}

{- $eigs

The Hmatrix eig factors A into P and D where A = P D inv(P) and D is diagonal.

The units for eigenvalues can be figured out:

>               _____
>      -1       |  c
> P D P  = A =  |r
>               |

>       _______
>       |   d
> P   = |c
>       |

>       _______
>       |   -1
>       |  c
>  -1   |
> P   = | -1
>       |d

So we can see that the dimension labeled `d-1` in P inverse is actually the
same `c` in `A`. The actual units of `d` don't seem to matter because the
`inv(d)` un-does any units that the `d` adds. So `d` can be all DOne. But
another choice, such as 1/c would be more appropriate, since then you can
expm your eigenvectors (not that that seems to be something people do)?

To get the row-units of A to match up, sometimes `D` will have units.
The equation ends up as D/c = r

Please ignore the type signatures on 'eig' 'eigC' etc. instead look at the type of
'wrapEig' 'wrapEigOnly' together with the hmatrix documentation (linked).

Perhaps the convenience definitions `eig m = wrapEig H.eig m` should be in
another module.
-}

{-

-- | 'wrapEig' H.'H.eig'
eig m = wrapEig H.eig m
-- | 'wrapEig' H.'H.eigC'
eigC m = wrapEig H.eigC m
-- | 'wrapEig' H.'H.eigH'
eigH m = wrapEig H.eigH m
-- | 'wrapEig' H.'H.eigH''
eigH' m = wrapEig H.eigH' m
-- | 'wrapEig' H.'H.eigR'
eigR m = wrapEig H.eigR m
-- | 'wrapEig' H.'H.eigS'
eigS m = wrapEig H.eigS m
-- | 'wrapEig' H.'H.eigS''
eigS' m = wrapEig H.eigS' m
-- | 'wrapEig' H.'H.eigSH'
eigSH m = wrapEig H.eigSH m
-- | 'wrapEig' H.'H.eigSH''
eigSH' m = wrapEig H.eigSH' m

-- | 'wrapEigOnly' H.'H.eigOnlyC'
eigOnlyC m = wrapEigOnly H.eigOnlyC m
-- | 'wrapEigOnly' H.'H.eigOnlyH'
eigOnlyH m = wrapEigOnly H.eigOnlyH m
-- | 'wrapEigOnly' H.'H.eigOnlyR'
eigOnlyR m = wrapEigOnly H.eigOnlyR m
-- | 'wrapEigOnly' H.'H.eigOnlyS'
eigOnlyS m = wrapEigOnly H.eigOnlyS m
-- | 'wrapEigOnly' H.'H.eigenvalues'
eigenvalues m = wrapEigOnly H.eigenvalues m
-- | 'wrapEigOnly' H.'H.eigenvaluesSH'
eigenvaluesSH m = wrapEigOnly H.eigenvaluesSH m
-- | 'wrapEigOnly' H.'H.eigenvaluesSH''
eigenvaluesSH' m = wrapEigOnly H.eigenvaluesSH' m

wrapEig :: (ones ~ (DOne ': _1), EigCxt [r,c] eigVal [cinv,ones],
    H.Field y, H.Field z)
    => (H.Matrix x -> (H.Vector y, H.Matrix z)) ->
    DimMat [r,c] x ->
    (DimMat '[eigVal] y, DimMat [cinv,ones] z)
wrapEig hmatrixFun (DimMat a) = case hmatrixFun a of
    (e,v) -> (DimVec e, DimMat v)

wrapEigOnly :: (EigCxt [r,c] eigVal '(), H.Field y)
    => (H.Matrix x -> H.Vector y) ->
    DimMat [r,c] x -> DimMat '[eigVal] y
wrapEigOnly hmatrixFun (DimMat a) = case hmatrixFun a of
    (e) -> DimVec e

-}

max_stars_repo_name: dmcclean/dimensional-dk-linalg

max_issues_repo_name: dmcclean/dimensional-dk-linalg

max_forks_repo_name: dmcclean/dimensional-dk-linalg


text:
{-# LANGUAGE Rank2Types, FlexibleContexts #-}

module Numeric.AD.Lagrangian.Internal where

import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Storable as S
import qualified Data.Packed.Vector as V
import qualified Data.Packed.Matrix as M

import GHC.IO (unsafePerformIO)

import Numeric.AD
import Numeric.Optimization.Algorithms.HagerZhang05
import Numeric.LinearAlgebra.Algorithms

-- | An equality constraint of the form @g(x, y, ...) = c@. Use '<=>' to
-- construct a 'Constraint'.
newtype Constraint = Constraint
  {unConstraint :: forall a. (Floating a) => ([a] -> a, a)}

infixr 1 <=>
-- | Build a 'Constraint' from a function and a constant
(<=>) :: (forall a. (Floating a) => [a] -> a)
      -> (forall b. (Floating b) => b)
      -> Constraint
f <=> c = Constraint (f, c)

-- | Numerically minimize the Langrangian. The objective function and each of
-- the constraints must take the same number of arguments.
minimize :: (forall a. (Floating a) => [a] -> a)
         -- ^ The objective function to minimize
         -> [Constraint]
         -- ^ A list of constraints @g \<=\> c@ corresponding to equations of
         -- the form @g(x, y, ...) = c@
         -> Double
         -- ^ Stop iterating when the largest component of the gradient is
         -- smaller than this value
         -> Int
         -- ^ The arity of the objective function, which must equal the arity of
         -- the constraints
         -> Either (Result, Statistics) (V.Vector Double, V.Vector Double)
         -- ^ Either a 'Right' containing the argmin and the Lagrange
         -- multipliers, or a 'Left' containing an explanation of why the
         -- gradient descent failed
minimize f constraints tolerance argCount = result where
    -- At a constrained minimum of `f`, the gradient of the Lagrangian must be
    -- zero. So we square the Lagrangian's gradient (making it non-negative) and
    -- minimize it.
    (sqGradLgn, gradSqGradLgn) = (fst . g, snd . g) where
        g = grad' $ squaredGrad $ lagrangian f constraints argCount

    -- Perhaps this should be exposed. ...
    guess = U.replicate (argCount + length constraints) 1

    result = case unsafePerformIO $
                    optimize
                        (defaultParameters {printFinal = False})
                        tolerance
                        guess
                        (VFunction (sqGradLgn . U.toList))
                        (VGradient (U.fromList . gradSqGradLgn . U.toList))
                        Nothing of
       (vs, ToleranceStatisfied, _) -> Right (S.take argCount vs,
                                              S.drop argCount vs)
       (_, x, y) -> Left (x, y)

-- | Numerically maximize the Langrangian. The objective function and each of
-- the constraints must take the same number of arguments.
maximize :: (forall a. (Floating a) => [a] -> a)
         -- ^ The objective function to minimize
         -> [Constraint]
         -- ^ A list of constraints @g \<=\> c@ corresponding to equations of
         -- the form @g(x, y, ...) = c@
         -> Double
         -- ^ Stop iterating when the largest component of the gradient is
         -- smaller than this value
         -> Int
         -- ^ The arity of the objective function, which must equal the arity of
         -- the constraints
         -> Either (Result, Statistics) (V.Vector Double, V.Vector Double)
         -- ^ Either a 'Right' containing the argmax and the Lagrange
         -- multipliers, or a 'Left' containing an explanation of why the
         -- gradient ascent failed
maximize f = minimize $ negate . f

lagrangian :: (Floating a)
           => (forall b. (Floating b) => [b] -> b)
           -> [Constraint]
           -> Int
           -> [a]
           -> a
lagrangian f constraints argCount argsAndLams = result where
    args = take argCount argsAndLams
    lams = drop argCount argsAndLams

    -- g(x, y, ...) = c <=> g(x, y, ...) - c = 0
    appliedConstraints = fmap (\(Constraint (g, c)) -> g args - c) constraints

    -- L(x, y, ..., lam0, ...) = f(x, y, ...) + lam0 * (g0 - c0) ...
    result = (f args) + (sum . zipWith (*) lams $ appliedConstraints)

squaredGrad :: (Floating a)
            => (forall b. (Floating b) => [b] -> b)
            -> [a]
            -> a
squaredGrad f = sum . fmap square . grad f where
    square x = x * x

-- | WARNING: Experimental.
--   This is not a true feasibility test for the function. I am not sure
--   exactly how to implement that. This just checks the feasiblility at a
--   point. If this ever returns false, 'solve' can fail.
feasible :: (Floating a, Field a, M.Element a)
         => (forall b. (Floating b) => [b] -> b)
         ->[Constraint]
         -> [a]
         -> Bool
feasible f constraints points = result where
    sqGradLgn :: (Floating a) => [a] -> a
    sqGradLgn = squaredGrad $ lagrangian f constraints $ length points

    hessianMatrix = M.fromLists . hessian sqGradLgn $ points

    -- make sure all of the eigenvalues are positive
    result = all (>0) . V.toList . eigenvaluesSH $ hessianMatrix

max_stars_repo_name: alexkalderimis/lagrangian

max_issues_repo_name: alexkalderimis/lagrangian

max_forks_repo_name: alexkalderimis/lagrangian


text:
{-# LANGUAGE BangPatterns #-}
module Main (main) where

import Data.Vector.Storable as V
import Numeric.FFT.Vector.Unnormalized as U
import Numeric.FFT.Vector.Invertible as I
import Numeric.FFT.Vector.Unitary as O
import Criterion.Main
import Data.Complex

main = do
    let mkSize planner k = do
            let !v = V.replicate k 17 :: Vector (Complex Double)
            let !p = plan planner k
            return $ bench (show k) $ whnf (\w -> execute p w) v
    let mkSize1 k = do
            let !v = V.replicate k 17 :: Vector (Complex Double)
            let !p = plan U.dft k
            return $ bench (show k) $ whnf (\w -> execute p w) v
    let mkSize2 k = do
            let !v = V.replicate k 17 :: Vector (Complex Double)
            let !p = plan O.dft k
            return $ bench (show k) $ whnf (\w -> execute p w) v

    let mkSize3 k = do
            let !v = V.replicate k 17 :: Vector (Complex Double)
            return $ bench (show k) $ whnf (\w -> I.run O.dft w) v

    let sizes = [64,256,2048,8192]
    let sizedGroup name f = fmap (bgroup name) $ Prelude.mapM f sizes
    benches <- sizedGroup "U.idft_1" (mkSize U.idft)
    benches1 <- sizedGroup "U.dft" mkSize1
    benches2 <- sizedGroup "O.dft" mkSize2
    benches3 <- sizedGroup "run_O.dft" mkSize3
    defaultMain $ [benches, benches1, benches2, benches3]

max_stars_repo_name: TravisWhitaker/vector-fftw

max_issues_repo_name: TravisWhitaker/vector-fftw

max_forks_repo_name: TravisWhitaker/vector-fftw


text:
{-# LANGUAGE RecordWildCards #-}

module NeuralNetwork.BackForwardPropagation
  ( genRandomThetas
  , computeCost
  , minimizeBFGS2
  ) where

import Control.Monad (mapM)
import Numeric.GSL.Minimization (MinimizeMethodD(VectorBFGS2), minimizeVD)
import Numeric.LinearAlgebra (( #> ), (<.>), (<>), sumElements)
import Numeric.LinearAlgebra.Data
  ( Matrix
  , R
  , Vector
  , (===)
  , (?)
  , asColumn
  , atIndex
  , cmap
  , flatten
  , fromList
  , fromLists
  , matrix
  , reshape
  , size
  , takesV
  , tr'
  , vjoin
  )
import Numeric.LinearAlgebra.Devel (mapMatrixWithIndex)
import Numeric.Natural (Natural)
import Prelude hiding ((<>))
import System.Random (randomRIO)

import Common
  ( MinimizationOpts(..)
  , addOnesColumn
  , getDimensions
  , sigmoid
  , sigmoidMatrix
  )
import NeuralNetwork.Types (Layer(numUnits), Network(..))

genRandomThetas :: R -> Network -> IO [Matrix R]
genRandomThetas epsilon Network {..} =
  let allLayers = inputLayer : hiddenLayer ++ [outputLayer]
   in mapM mkTheta (zip allLayers (tail allLayers))
  where
    mkTheta :: (Layer, Layer) -> IO (Matrix R)
    mkTheta (l1, l2) = do
      let numRows = numUnits l2
          numCols = numUnits l1 + 1
          numElts = fromIntegral $ numRows * numCols
      elts <- sequence $ replicate numElts (randomRIO (0.0, 1.0))
      return $
        matrix
          (fromIntegral numCols)
          ((\r -> r * (2 * epsilon) - epsilon) <$> elts)

computeCost :: Int -> Matrix R -> Vector R -> [Matrix R] -> R
computeCost numLabels x y thetas = cost
  where
    m = size y
    --
    hypotheses =
      foldl (\an th -> sigmoidMatrix ((addOnesColumn an) <> (tr' th))) x thetas
    --
    cost =
      let rowElemsCost i k r =
            let yik =
                  if (round $ y `atIndex` i) == k
                    then 1.0
                    else 0.0
             in (-yik * (log (r))) - ((1 - yik) * log (1 - r))
          summation =
            sumElements $
            mapMatrixWithIndex (\(i, k) r -> rowElemsCost i k r) hypotheses
       in summation / fromIntegral m

minimizeBFGS2 ::
     Matrix R
  -> Vector R
  -> Int
  -> [Matrix R]
  -> Int
  -> MinimizationOpts
  -> Int
  -> [Matrix R]
minimizeBFGS2 x y numLabels thetas numIters MinimizationOpts {..} regFactor =
  let (flattenedParams, dimensions) = flattenParameters thetas
      (minParams, _) =
        minimizeVD
          VectorBFGS2
          precision
          numIters
          sizeOfFirstTrialStep
          tolerance
          (computeCost' dimensions)
          (computeGradients' dimensions)
          flattenedParams
   in unflattenParameters dimensions minParams
  where
    m = size y
    --
    flattenParameters :: [Matrix R] -> (Vector R, [(Int, Int)])
    flattenParameters mtxs =
      let (vecs, dims) =
            unzip $ (\mx -> (flatten mx, getDimensions mx)) <$> mtxs
       in (vjoin vecs, dims)
    --
    unflattenParameters :: [(Int, Int)] -> Vector R -> [Matrix R]
    unflattenParameters dims params =
      let sizes = uncurry (*) <$> dims
       in zipWith
            (\vec (ncols, _) -> reshape ncols vec)
            (takesV sizes params)
            dims
    --
    computeCost' :: [(Int, Int)] -> Vector R -> R
    computeCost' dims params =
      computeCost numLabels x y (unflattenParameters dims params)
    --
    computeGradients' :: [(Int, Int)] -> Vector R -> Vector R
    computeGradients' dims gradients =
      fst $
      flattenParameters $
      computeGradients 0 (unflattenParameters dims gradients)
    --
    computeActivationVals :: Matrix R -> [Matrix R]
    computeActivationVals a0 =
      let f vals@[] th = sigmoidMatrix (th * (addBiasUnit a0)) : vals
          f vals@(an:_) th = sigmoidMatrix (th * (addBiasUnit an)) : vals
       in foldl f [] thetas
    --
    computeErrors :: Matrix R -> [Matrix R] -> [Matrix R]
    computeErrors dL activationVals =
      let f vals@(dl:_) (an, th) = ((tr' th) * dl) .* an .* (1 - an) : vals
       in foldl f [dL] (zip activationVals (reverse thetas))
    --
    computeGradients :: Int -> [Matrix R] -> [Matrix R]
    computeGradients i gradients
      | i >= m = gradients
      | otherwise =
        let xi = (asColumn . flatten) $ x ? [i]
            a0 = xi
            activationVals@(aL:_) = computeActivationVals a0
            yi = round $ y `atIndex` i
            binEq :: Int -> R
            binEq n =
              if yi == n
                then 1.0
                else 0.0
            yik = asColumn $ fromList (binEq <$> [0 .. (numLabels - 1)])
            errorVals = computeErrors (aL - yik) (tail activationVals)
            accumGradient grad errVal actVal = grad + (errVal * (tr' actVal))
         in computeGradients (i + 1) $
            zipWith3
              accumGradient
              gradients
              errorVals
              (a0 : reverse activationVals)

elemWiseMult :: Matrix R -> Matrix R -> Matrix R
elemWiseMult mx1 mx2 =
  mapMatrixWithIndex (\idx r -> (mx1 `atIndex` idx) * r) mx2

(.*) :: Matrix R -> Matrix R -> Matrix R
(.*) = elemWiseMult

infixr 8 .*

-- Helpers
addBiasUnit :: Matrix R -> Matrix R
addBiasUnit mx = fromLists [[1]] === mx

addBiasUnit' :: Vector R -> Vector R
addBiasUnit' vec = vjoin [fromList [1], vec]

max_stars_repo_name: sebashack/mlToolBox

max_issues_repo_name: sebashack/mlToolBox

max_forks_repo_name: sebashack/mlToolBox


text:
module Main where
-- package の import
-- CSVの読み込み
import qualified    Text.CSV                                as CSV
import qualified    Data.Text                               as T
import              Data.Attoparsec.Text                    as DAT
-- 統計のためのパッケージ
import qualified    Statistics.Correlation                  as SC
-- その他
import qualified    Control.Monad                           as CM
import qualified    Data.List                               as L
import qualified    Data.Vector.Unboxed                     as VU
import              Data.IORef

main = do
    -- csvファイルの読み込み
    xs <- read_csv "Data/multiTimeline.csv"
    -- データの表示
    print_data (L.head xs) (L.tail xs)
    -- Headerの抽出
    let headers = L.map T.unpack $ L.head xs
    -- 値のDoubleへの変換
    let body    = L.transpose
                $ L.map  (L.map parseDouble)
                $ L.tail xs
    -- 分析するデータの選択
    -- xとyのHeader名を書き換えよう
    let x = column body "AI" headers
    let y = column body "Python" headers
    -- 相関係数を求める
    let res = SC.pearson
            $ VU.zip (VU.fromList x) (VU.fromList y)
    -- 相関係数の表示
    putStrLn $ "相関係数:" ++ show res

-- ** 作業用関数
-- 別のモジュールにまとめておくのが良いが複雑になるのでここにまとめておく
-- | Doubleのパーサ
{-# INLINE parseDouble #-}
parseDouble :: T.Text -> Double
parseDouble tx = case DAT.parseOnly DAT.double tx of
        Right r   -> r
        Left  l   -> error $ "Error on parseDouble : " ++ show l

-- | データの読み込み
read_csv :: String -> IO [[T.Text]]
read_csv file
    = CSV.parseCSVFromFile file >>= \res
    -> case res of
        Right xs -> return $ L.map (L.map T.pack) xs
        Left  xs -> error $ "error on parseCSVFromFile: " ++ show xs

-- | data set っぽく表示する
print_data :: (Show a)  => [T.Text] {- header -}
                        -> [[a]]    {- body -}
                        -> IO ()
print_data hs bs = do
    let max_length  =  let hs_max = L.maximum (L.map T.length hs)
                    in let bs_max = L.maximum
                                  $ (flip L.map) bs
                                  $ \x  -> L.maximum
                                  $ L.map (L.length . show) x
                    in case compare hs_max bs_max  of
                        LT -> bs_max
                        _  -> hs_max

    let bs_digit_num =  L.length (show (L.length bs))
    putStr $ " " ++ L.replicate bs_digit_num ' '
    CM.forM_ hs $ \x -> let tx = T.unpack x
                     in case compare (L.length tx) max_length of
                        LT -> putStr $ tx ++ (L.replicate  (max_length -  (L.length tx)) ' ')  ++ " "
                        _  -> putStr $ tx ++ " "
    putStrLn ""

    currentNum <- newIORef 0 :: IO (IORef Int)

    let  total_col = L.length bs
    case total_col  > 10 of
        False -> CM.forM_ bs $  \xs -> readIORef currentNum >>= \num
                                    -> putStr ((show num)   ++ " "
                                                            ++ L.replicate
                                                                (bs_digit_num - (L.length (show num))) ' ')
                                    >> modifyIORef currentNum (\x -> x + 1)
                                    >> CM.forM_ xs  ( \x -> let tx = show x
                                                         in putStr  $ tx
                                                                    ++ (L.replicate  (max_length -  (L.length tx)) ' ')
                                                                    ++ " ")
                                    >> putStrLn ""
        True ->  CM.forM_ bs $  \xs -> readIORef currentNum >>= \num
                                    -> case (num <= 4 || (total_col - num) <= 4) of
                                        True  -> putStr ((show num)   ++ " " ++ L.replicate (bs_digit_num - (L.length (show num))) ' ')
                                              >> modifyIORef currentNum (\x -> x + 1)
                                              >> CM.forM_ xs  ( \x -> let tx = show x
                                                         in putStr  $ tx
                                                                    ++ (L.replicate  (max_length -  (L.length tx)) ' ')
                                                                    ++ " ")
                                              >> putStrLn ""
                                        False -> modifyIORef currentNum (\x -> x + 1)

-- | headerによる情報の抜き出し
column :: [[Double]] -> String -> [String] -> [Double]
column xs name hs = case (L.elemIndex name hs) of
                    Just i  -> xs !! i
                    Nothing -> error $ "error at column: " ++ name ++ "doesn't exest."


text:
{-# OPTIONS_GHC -fno-warn-orphans #-}
module Main where

import MNIST.Prelude
import MNIST.DataSet
import Numeric.LinearAlgebra.Static hiding (dot)
import Numeric.Backprop

import qualified Data.Vector.Generic             as VG
import qualified Data.Vector                     as V
import qualified Generics.SOP                    as SOP
import qualified Numeric.LinearAlgebra           as HM
import qualified System.Random.MWC               as MWC
import qualified System.Random.MWC.Distributions as MWC

-- | type-combinators alias for the terminal constructor
nil :: Prod f '[]
nil = Ø

-- ========================================================================= --
-- Internal structure of our FF neural net

data Layer i o = Layer
  { _lWeights :: !(L o i)
  , _lBiases  :: !(R o)
  } deriving (Show, Generic)

instance SOP.Generic (Layer i o)
instance NFData (Layer i o)

data Network i h1 h2 o = Net
  { _nLayer1 :: !(Layer i  h1)
  , _nLayer2 :: !(Layer h1 h2)
  , _nLayer3 :: !(Layer h2 o)
  } deriving (Show, Generic)

instance SOP.Generic (Network i h1 h2 o)
instance NFData (Network i h1 h2 o)

-- ========================================================================= --
-- type classes

type KnownNat4 i h1 h2 o = (KnownNat i, KnownNat h1, KnownNat h2, KnownNat o)
type KnownNat2 i o = (KnownNat i, KnownNat o)

instance KnownNat2 i o => Num (Layer i o) where
  Layer w1 b1 + Layer w2 b2 = Layer (w1 + w2) (b1 + b2)
  Layer w1 b1 - Layer w2 b2 = Layer (w1 - w2) (b1 - b2)
  Layer w1 b1 * Layer w2 b2 = Layer (w1 * w2) (b1 * b2)
  abs    (Layer w b)        = Layer (abs    w) (abs    b)
  signum (Layer w b)        = Layer (signum w) (signum b)
  negate (Layer w b)        = Layer (negate w) (negate b)
  fromInteger x = Layer (fromInteger x) (fromInteger x)

instance KnownNat4 i h1 h2 o => Num (Network i h1 h2 o) where
  Net a b c + Net d e f = Net (a + d) (b + e) (c + f)
  Net a b c - Net d e f = Net (a - d) (b - e) (c - f)
  Net a b c * Net d e f = Net (a * d) (b * e) (c * f)
  abs    (Net a b c)    = Net (abs    a) (abs    b) (abs    c)
  signum (Net a b c)    = Net (signum a) (signum b) (signum c)
  negate (Net a b c)    = Net (negate a) (negate b) (negate c)
  fromInteger x         = Net (fromInteger x) (fromInteger x) (fromInteger x)

instance KnownNat2 i o => Fractional (Layer i o) where
  Layer w1 b1 / Layer w2 b2 = Layer (w1 / w2) (b1 / b2)
  recip (Layer w b)         = Layer (recip w) (recip b)
  fromRational x            = Layer (fromRational x) (fromRational x)

instance KnownNat4 i h1 h2 o => Fractional (Network i h1 h2 o) where
  Net a b c / Net d e f = Net (a / d) (b / e) (c / f)
  recip (Net a b c)     = Net (recip a) (recip b) (recip c)
  fromRational x        = Net (fromRational x) (fromRational x) (fromRational x)

instance KnownNat n => MWC.Variate (R n) where
  uniform g = randomVector <$> MWC.uniform g <*> pure Uniform
  uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g

instance KnownNat2 m n => MWC.Variate (L m n) where
  uniform g = uniformSample <$> MWC.uniform g <*> pure 0 <*> pure 1
  uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g

instance KnownNat2 i o => MWC.Variate (Layer i o) where
  uniform g = Layer <$> MWC.uniform g <*> MWC.uniform g
  uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g

instance KnownNat4 i h1 h2 o => MWC.Variate (Network i h1 h2 o) where
  uniform g = Net <$> MWC.uniform g <*> MWC.uniform g <*> MWC.uniform g
  uniformR (l, h) g = (\x -> x * (h - l) + l) <$> MWC.uniform g

-- ========================================================================= --
-- Basic math functions with back propagation

matVec :: forall m n . (KnownNat m, KnownNat n) => Op '[ L m n, R n ] (R m)
matVec = op2' $ \m v -> (forward m v, backward m v)
  where
    forward :: L m n -> R n -> R m
    forward m v = m #> v

    backward :: L m n -> R n -> Maybe (R m) -> (L m n, R n)
    backward m v (fromMaybe 1 -> g) = (g `outer` v, tr m #> g)

dot :: forall n . KnownNat n => Op '[ R n, R n ] Double
dot = op2' $ \x y -> (forward x y, backward x y)
  where
    forward :: R n -> R n -> Double
    forward x y = x <.> y

    backward :: R n -> R n -> Maybe Double -> (R n, R n)
    backward x y = \case
      Nothing -> (y, x)
      Just g  -> (konst g * y, x * konst g)

scale :: forall n . KnownNat n => Op '[ Double, R n ] (R n)
scale = op2' $ \a x -> (forward a x, backward a x)
  where
    forward :: Double -> R n -> R n
    forward a x = konst a * x

    backward :: Double -> R n -> Maybe (R n) -> (Double, R n)
    backward a x = \case
      Nothing -> (HM.sumElements (extract x      ), konst a    )
      Just g  -> (HM.sumElements (extract (x * g)), konst a * g)

vsum :: forall n . KnownNat n => Op '[ R n ] Double
vsum = op1' $ \x -> (forward x, backward)
  where
    forward :: R n -> Double
    forward = HM.sumElements . extract

    backward :: Maybe Double -> R n
    backward = maybe 1 konst

logistic :: Floating a => a -> a
logistic x = 1 / (1 + exp (-x))

-- ========================================================================= --
-- run backpropagation
type LayerCtx i o = '[R i, Layer i o]
type NetCtx i h1 h2 o = '[ R i, Network i h1 h2 o ]

runLayer :: forall i o s . (KnownNat i, KnownNat o) => BPOp s '[ R i, Layer i o ] (R o)
runLayer = withInps $ \(decombinate -> (x, l)) -> do
  w :< b :< _ <- partsVar gTuple l
  y <- opVar matVec (w :< x :< nil)
  return $ y + b
  where
    decombinate
      :: Prod (BVar s (LayerCtx i o)) '[R i, Layer i o]
      -> (BVar s (LayerCtx i o) (R i), BVar s (LayerCtx i o) (Layer i o))
    decombinate (x :< l :< _) = (x, l)

runNetwork :: KnownNat4 i h1 h2 o => BPOp s (NetCtx i h1 h2 o) (R o)
runNetwork = withInps $ \(x :< n :< _) -> do
  l1 :< l2 :< l3 :< _ <- partsVar gTuple n

  y <- bindVar $ liftB (bpOp runLayer) (x :< l1 :< nil)
  -- or
  z <- bindVar $ (bpOp runLayer) .$ (logistic y :< l2 :< nil)
  -- or
  r <- (bpOp runLayer) ~$ (logistic z :< l3 :< nil)

  bpOp softmax ~$ only r
  where
    softmax :: KnownNat n => BPOp s '[ R n ] (R n)
    softmax = withInps $ \(x :< _) -> do
      expX <- bindVar (exp x)
      totX <- vsum ~$ (expX :< nil)
      scale        ~$ (1 / totX :< expX :< nil)

runNetOnInp :: KnownNat4 i h1 h2 o => Network i h1 h2 o -> R i -> R o
runNetOnInp n x = evalBPOp runNetwork (x ::< n ::< nil)

gradNet :: KnownNat4 i h1 h2 o => Network i h1 h2 o -> R i -> Network i h1 h2 o
gradNet n x = case gradBPOp runNetwork (x ::< n ::< nil) of
    _gradX ::< gradN ::< nil -> gradN

-- ========================================================================= --

-- crossEntropy :: forall s n . KnownNat n => R n -> BPOp s '[ R n ] Double
-- crossEntropy targ = withInps $ \(r :< _) ->
--   negate (dot ~$ (log r :< only t))
--   where
--     t :: BVar s '[R n] (R n)
--     t = constVar targ

crossEntropyI :: forall s n . KnownNat n => R n -> BPOpI s '[ R n ] Double
crossEntropyI targ (r :< _) = negate (dot .$ (log r :< only t))
  where
    t :: BVar s '[R n] (R n)
    t = constVar targ

softMaxCrossEntropy :: forall s n . KnownNat n => R n -> BPOp s '[ R n ] Double
softMaxCrossEntropy targ = withInps $ \(r :< Ø) -> do
  bindVar $ realToFrac tsum * log (vsum .$ (only r)) - (dot .$ (r :< t :< nil))
  where
    tsum :: Double
    tsum = HM.sumElements . extract $ targ

    t :: BVar s '[R n] (R n)
    t = constVar targ

softMaxCrossEntropyI :: forall s n . KnownNat n => R n -> BPOpI s '[ R n ] Double
softMaxCrossEntropyI targ (r :< Ø) =
  realToFrac tsum * log (vsum .$ (only r)) - (dot .$ (r :< t :< nil))
  where
    tsum :: Double
    tsum = HM.sumElements . extract $ targ

    t :: BVar s '[R n] (R n)
    t = constVar targ

trainStep
  :: forall i h1 h2 o. KnownNat4 i h1 h2 o
  => Double
  -> R i
  -> R o
  -> Network i h1 h2 o
  -> Network i h1 h2 o
trainStep r !x !t !n =
  case gradBPOp o (x ::< n ::< nil) of
    (_ :< I gN :< _) -> n - (realToFrac r * gN)
  where
    o :: BPOp s '[ R i, Network i h1 h2 o ] Double
    o = do
      y <- runNetwork
      implicitly (crossEntropyI t) -$ (y :< nil)

trainList
  :: KnownNat4 i h1 h2 o
  => Double
  -> [(R i, R o)]
  -> Network i h1 h2 o
  -> Network i h1 h2 o
trainList r = flip $ foldl' (\n (x,y) -> trainStep r x y n)

testNet
  :: forall i h1 h2 o. KnownNat4 i h1 h2 o
  => [(R i, R o)]
  -> Network i h1 h2 o
  -> Double
testNet xs n = sum (map (\(i,o) -> test i o) xs) / fromIntegral (length xs)
  where
    test :: R i -> R o -> Double
    test x (extract->t) = fromIntegral . fromEnum $
      HM.maxIndex t == HM.maxIndex (extract r)
      where
        r :: R o
        r = evalBPOp runNetwork (x ::< n ::< nil)

main :: IO ()
main = MWC.withSystemRandom $ \g -> do
  -- initialize data and network
  !trainingSet <- trainingDataBp
  !testSet     <- testDataBp
  !net0        <- MWC.uniformR @(Network 784 300 100 9) (-0.5, 0.5) g

  flip evalStateT net0 . forM_ [1..100] $ \e -> do
    trainingSet' <- liftIO . fmap V.toList $ MWC.uniformShuffle (V.fromList trainingSet) g
    liftIO $ printf "[Epoch %d]\n" (e :: Int)

    forM_ ([1..] `zip` chunksOf batch trainingSet') $ \(b, chnk) -> StateT $ \n0 -> do
      printf "(Batch %d)\n" (b :: Int)

      -- t0 <- getCurrentTime
      n' <- evaluate . force $ trainList rate chnk n0
      -- t1 <- getCurrentTime
      -- printf "Trained on %d points in %s.\n" batch (show (t1 `diffUTCTime` t0))

      let trainScore = testNet chnk    n'
          testScore  = testNet testSet n'
      printf "Training error:   %.2f%%\n" ((1 - trainScore) * 100)
      -- printf "Validation error: %.2f%%\n" ((1 - testScore ) * 100)

      return ((), n')
  where
    rate = 0.0001
    batch = 100

  --  go :: StateT (Network 784 300 100 9) IO ()
  --  go = do
  --    e <- get

--
--   -- Test
--   let testPreds = map (take 3 testSet) $ \(x, _) -> evalBPOp runNetwork (x ::< net0 ::< nil)
--
--   liftIO $ forM_ ([0..3] :: [Int]) $ \i -> do
--       -- T.putStrLn $ drawMNIST $ testImages !! i
--       putStrLn $ "\n" ++ "expected " ++ show (testLabels !! i)
--       putStrLn $         "     got " ++ show (testPreds !! i)
--
--

max_stars_repo_name: stites/haskell-mnist-benchmarks

max_issues_repo_name: stites/haskell-mnist-benchmarks

max_forks_repo_name: stites/haskell-mnist-benchmarks