text:
{-# LANGUAGE ConstraintKinds         #-}
{-# LANGUAGE DataKinds               #-}
{-# LANGUAGE DefaultSignatures       #-}
{-# LANGUAGE FlexibleContexts        #-}
{-# LANGUAGE FlexibleInstances       #-}
{-# LANGUAGE KindSignatures          #-}
{-# LANGUAGE MagicHash               #-}
{-# LANGUAGE MultiParamTypeClasses   #-}
{-# LANGUAGE PolyKinds               #-}
{-# LANGUAGE RankNTypes              #-}
{-# LANGUAGE ScopedTypeVariables     #-}
{-# LANGUAGE TypeApplications        #-}
{-# LANGUAGE TypeFamilies            #-}
{-# LANGUAGE TypeOperators           #-}
{-# LANGUAGE UndecidableInstances    #-}
{-# LANGUAGE UndecidableSuperClasses #-}
module Data.Vector.HFixed.Class (
    -- * Types and type classes
    -- ** N-ary functions
    Fn
  , Fun
  , TFun(..)
    -- ** Type functions
  , Proxy(..)
  , type (++)
  , Len
  , HomList
    -- ** Type classes
  , Arity(..)
  , ArityC(..)
  , (:&&:)
  , HVector(..)
  , tupleSize
  , HVectorF(..)
  , tupleSizeF
    -- *** Lookup in vector
  , Index(..)
  , TyLookup(..)
    -- ** CPS-encoded vector
  , ContVec
  , ContVecF(..)
  , cons
  , consF
    -- ** Interop with homogeneous vectors
  , HomArity(..)
  , homInspect
  , homConstruct
    -- * Operations of Fun
    -- ** Primitives for Fun
  , curryFun
  , uncurryFun
  , uncurryMany
  , curryMany
  , constFun
    -- ** Primitives for TFun
  , constTFun
  , curryTFun
  , uncurryTFun
  , shuffleTF
  , stepTFun
    -- ** More complicated functions
  , concatF
  , lensWorkerF
  , lensWorkerTF
    -- * Lens
  , Lens
  , Lens'
  ) where

import Data.Coerce
import Data.Complex          (Complex(..))
import Data.Functor.Identity (Identity(..))
import Data.Type.Equality    (type (==))

import           Data.Vector.Fixed.Cont   (Peano,PeanoNum(..),ArityPeano)
import qualified Data.Vector.Fixed                as F
import qualified Data.Vector.Fixed.Cont           as F (curryFirst)
import qualified Data.Vector.Fixed.Unboxed        as U
import qualified Data.Vector.Fixed.Primitive      as P
import qualified Data.Vector.Fixed.Storable       as S
import qualified Data.Vector.Fixed.Boxed          as B

import Unsafe.Coerce (unsafeCoerce)
import GHC.Exts (Proxy#,proxy#)
import GHC.TypeLits
import GHC.Generics hiding (S)

import Data.Vector.HFixed.TypeFuns

----------------------------------------------------------------
-- Types
----------------------------------------------------------------

-- | Type family for N-ary function. Types of function parameters are
--   encoded as the list of types.
type family Fn (f :: α -> *) (as :: [α]) b where
  Fn f '[]      b = b
  Fn f (a : as) b = f a -> Fn f as b

-- | Newtype wrapper for function where all type parameters have same
--   type constructor. This type is required for writing function
--   which works with monads, appicatives etc.
newtype TFun f as b = TFun { unTFun :: Fn f as b }

-- | Newtype wrapper to work around of type families' lack of
--   injectivity.
type Fun = TFun Identity

----------------------------------------------------------------
-- Generic operations
----------------------------------------------------------------

-- | Type class for combining two constraint constructors. Those are
--   required for 'ArityC' type class.
class (c1 a, c2 a) => (:&&:) c1 c2 a

instance (c1 a, c2 a) => (:&&:) c1 c2 a

-- | Type class for dealing with N-ary function in generic way. Both
--   'accum' and 'apply' work with accumulator data types which are
--   polymorphic. So it's only possible to write functions which
--   rearrange elements in vector using plain ADT. It's possible to
--   get around it by using GADT as accumulator (See 'ArityC' and
--   function which use it)
--
--   This is also somewhat a kitchen sink module. It contains
--   witnesses which could be used to prove type equalities or to
--   bring instance in scope.
class Arity (xs :: [α]) where
  -- | Fold over /N/ elements exposed as N-ary function.
  accum :: (forall a as. t (a : as) -> f a -> t as)
        -- ^ Step function. Applies element to accumulator.
        -> (t '[] -> b)
        -- ^ Extract value from accumulator.
        -> t xs
        -- ^ Initial state.
        -> TFun f xs b

  -- | Apply values to N-ary function
  apply :: (forall a as. t (a : as) -> (f a, t as))
        -- ^ Extract value to be applied to function.
        -> t xs
        -- ^ Initial state.
        -> ContVecF xs f

  -- | Size of type list as integer.
  arity :: p xs -> Int

class (Arity xs) => ArityC c xs where
  accumC :: proxy c
         -- ^
         -> (forall a as. (c a) => t (a : as) -> f a -> t as)
         -- ^ Step function. Applies element to accumulator.
         -> (t '[] -> b)
         -- ^ Extract value from accumulator.
         -> t xs
         -- ^ Initial state.
         -> TFun f xs b

  -- | Apply values to N-ary function
  applyC :: proxy c
         --
         -> (forall a as. (c a) => t (a : as) -> (f a, t as))
         -- ^ Extract value to be applied to function.
         -> t xs
         -- ^ Initial state.
         -> ContVecF xs f

instance Arity '[] where
  accum _ f t = TFun (f t)
  apply _ _   = ContVecF unTFun
  {-# INLINE accum #-}
  {-# INLINE apply #-}
  arity _     = 0
  {-# INLINE arity #-}

instance Arity xs => Arity (x : xs) where
  accum f g t = uncurryTFun (\a -> accum f g (f t a))
  apply f t   = case f t of (a,u) -> consF a (apply f u)
  {-# INLINE accum #-}
  {-# INLINE apply #-}
  arity _     = 1 + arity (Proxy :: Proxy xs)
  {-# INLINE arity        #-}

instance ArityC c '[] where
  accumC _ _ f t = TFun (f t)
  applyC _ _ _   = ContVecF unTFun
  {-# INLINE accumC #-}
  {-# INLINE applyC #-}

instance (c x, ArityC c xs) => ArityC c (x : xs) where
  accumC w f g t = uncurryTFun (\a -> accumC w f g (f t a))
  applyC w f t   = case f t of (a,u) -> consF a (applyC w f u)
  {-# INLINE accumC #-}
  {-# INLINE applyC #-}

-- |
-- Type class for product type. Any product type could have instance
-- of this type.  Its methods describe how to construct and
-- deconstruct data type. For example instance for simple data type
-- with two fields could be written as:
--
-- > data A a = A Int a
-- >
-- > instance HVector (A a) where
-- >   type Elems (A a) = '[Int,a]
-- >   construct = TFun $ \i a -> A i a
-- >   inspect (A i a) (TFun f) = f i a
--
-- Another equivalent description of this type class is descibes
-- isomorphism between data type and
-- 'Data.Vector.HFixed.Cont.ContVec', where @constuct@ implements
-- @ContVec → a@ (see 'Data.Vector.HFixed.Cont.vector') and @inspect@
-- implements @a → ContVec@ (see 'Data.Vector.HFixed.Cont.cvec')
--
-- Istances should satisfy one law:
--
-- > inspect v construct = v
--
-- Default implementation which uses 'Generic' is provided.
class Arity (Elems v) => HVector v where
  type Elems v :: [*]
  type Elems v = GElems (Rep v)
  -- | Function for constructing vector
  construct :: Fun (Elems v) v
  default construct :: (Generic v, GHVector (Rep v), GElems (Rep v) ~ Elems v)
                    => Fun (Elems v) v
  construct = fmap to gconstruct
  -- | Function for deconstruction of vector. It applies vector's
  --   elements to N-ary function.
  inspect :: v -> Fun (Elems v) a -> a
  default inspect :: (Generic v, GHVector (Rep v), GElems (Rep v) ~ Elems v)
                  => v -> Fun (Elems v) a -> a
  inspect v = ginspect (from v)
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

-- | Number of elements in product type
tupleSize :: forall v proxy. HVector v => proxy v -> Int
tupleSize _ = arity (Proxy :: Proxy (Elems v))

-- | Type class for partially homogeneous vector where every element
--   in the vector have same type constructor. Vector itself is
--   parametrized by that constructor
class Arity (ElemsF v) => HVectorF (v :: (α -> *) -> *) where
  -- | Elements of the vector without type constructors
  type ElemsF v :: [α]
  inspectF   :: v f -> TFun f (ElemsF v) a -> a
  constructF :: TFun f (ElemsF v) (v f)

-- | Number of elements in parametrized product type
tupleSizeF :: forall v f proxy. HVectorF v => proxy (v f) -> Int
tupleSizeF _a = arity (Proxy :: Proxy (ElemsF v))

----------------------------------------------------------------
-- Interop with homogeneous vectors
----------------------------------------------------------------

-- | Conversion between homogeneous and heterogeneous N-ary functions.
class (ArityPeano n, Arity (HomList n a)) => HomArity n a where
  -- | Convert n-ary homogeneous function to heterogeneous.
  toHeterogeneous :: F.Fun n a r -> Fun (HomList n a) r
  -- | Convert heterogeneous n-ary function to homogeneous.
  toHomogeneous   :: Fun (HomList n a) r -> F.Fun n a r

instance HomArity 'Z a where
  toHeterogeneous = coerce
  toHomogeneous   = coerce
  {-# INLINE toHeterogeneous #-}
  {-# INLINE toHomogeneous   #-}

instance HomArity n a => HomArity ('S n) a where
  toHeterogeneous f
    = coerce $ \a -> unTFun $ toHeterogeneous (F.curryFirst f a)
  toHomogeneous (f :: Fun (a : HomList n a) r)
    = coerce $ \a -> (toHomogeneous $ curryFun f a :: F.Fun n a r)
  {-# INLINE toHeterogeneous #-}
  {-# INLINE toHomogeneous   #-}

-- | Default implementation of 'inspect' for homogeneous vector.
homInspect :: (F.Vector v a, HomArity (Peano (F.Dim v)) a)
           => v a -> Fun (HomList (Peano (F.Dim v)) a) r -> r
homInspect v f = F.inspect v (toHomogeneous f)
{-# INLINE homInspect #-}

-- | Default implementation of 'construct' for homogeneous vector.
homConstruct :: forall v a.
                (F.Vector v a, HomArity (Peano (F.Dim v)) a)
             => Fun (HomList (Peano (F.Dim v)) a) (v a)
homConstruct = toHeterogeneous (F.construct :: F.Fun (Peano (F.Dim v)) a (v a))
{-# INLINE homConstruct #-}

instance ( HomArity (Peano n) a
         , KnownNat n
         , Peano (n + 1) ~ 'S (Peano n)
         ) => HVector (B.Vec n a) where
  type Elems (B.Vec n a) = HomList (Peano n) a
  inspect   = homInspect
  construct = homConstruct
  {-# INLINE inspect   #-}
  {-# INLINE construct #-}

instance ( U.Unbox n a
         , HomArity (Peano n) a
         , KnownNat n
         , Peano (n + 1) ~ 'S (Peano n)
         ) => HVector (U.Vec n a) where
  type Elems (U.Vec n a) = HomList (Peano n) a
  inspect   = homInspect
  construct = homConstruct
  {-# INLINE inspect   #-}
  {-# INLINE construct #-}

instance ( S.Storable a
         , HomArity (Peano n) a
         , KnownNat n
         , Peano (n + 1) ~ 'S (Peano n)
         ) => HVector (S.Vec n a) where
  type Elems (S.Vec n a) = HomList (Peano n) a
  inspect   = homInspect
  construct = homConstruct
  {-# INLINE inspect   #-}
  {-# INLINE construct #-}

instance ( P.Prim a
         , HomArity (Peano n) a
         , KnownNat n
         , Peano (n + 1) ~ 'S (Peano n)
         ) => HVector (P.Vec n a) where
  type Elems (P.Vec n a) = HomList (Peano n) a
  inspect   = homInspect
  construct = homConstruct
  {-# INLINE inspect   #-}
  {-# INLINE construct #-}

----------------------------------------------------------------
-- CPS-encoded vectors
----------------------------------------------------------------

--
-- newtype ContVec xs = ContVec { runContVec :: forall r. Fun xs r -> r }

instance Arity xs => HVector (ContVecF xs Identity) where
  type Elems (ContVecF xs Identity) = xs
  construct = accum
    (\(T_mkN f) (Identity x) -> T_mkN (f . cons x))
    (\(T_mkN f)              -> f (ContVecF unTFun))
    (T_mkN id)
  inspect (ContVecF cont) f = cont f
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

newtype T_mkN all xs = T_mkN (ContVec xs -> ContVec all)

-- | CPS-encoded heterogeneous vector.
type ContVec xs = ContVecF xs Identity

-- | CPS-encoded partially heterogeneous vector.
newtype ContVecF (xs :: [α]) (f :: α -> *) =
  ContVecF { runContVecF :: forall r. TFun f xs r -> r }

instance Arity xs => HVectorF (ContVecF xs) where
  type ElemsF (ContVecF xs) = xs
  inspectF (ContVecF cont) = cont
  constructF = constructFF
  {-# INLINE constructF #-}
  {-# INLINE inspectF   #-}

constructFF :: forall f xs. (Arity xs) => TFun f xs (ContVecF xs f)
{-# INLINE constructFF #-}
constructFF = accum (\(TF_mkN f) x -> TF_mkN (f . consF x))
                    (\(TF_mkN f)   -> f $ ContVecF unTFun)
                    (TF_mkN id)

newtype TF_mkN f all xs = TF_mkN (ContVecF xs f -> ContVecF all f)

-- | Cons element to the vector
cons :: x -> ContVec xs -> ContVec (x : xs)
cons x (ContVecF cont) = ContVecF $ \f -> cont $ curryFun f x
{-# INLINE cons #-}

-- | Cons element to the vector
consF :: f x -> ContVecF xs f -> ContVecF (x : xs) f
consF x (ContVecF cont) = ContVecF $ \f -> cont $ curryTFun f x
{-# INLINE consF #-}

----------------------------------------------------------------
-- Instances of Fun
----------------------------------------------------------------

instance (Arity xs) => Functor (TFun f xs) where
  fmap f (TFun g0)
    = accum (\(TF_fmap g) a -> TF_fmap (g a))
            (\(TF_fmap r)   -> f r)
            (TF_fmap g0)
  {-# INLINE fmap #-}

instance (Arity xs) => Applicative (TFun f xs) where
  pure r = accum (\Proxy _ -> Proxy)
                 (\Proxy   -> r)
                 (Proxy)
  {-# INLINE pure  #-}
  (TFun f0 :: TFun f xs (a -> b)) <*> (TFun g0 :: TFun f xs a)
    = accum (\(TF_ap f g) a -> TF_ap (f a) (g a))
            (\(TF_ap f g)   -> f g)
            ( TF_ap f0 g0 :: TF_ap f (a -> b) a xs)
  {-# INLINE (<*>) #-}

instance Arity xs => Monad (TFun f xs) where
  return  = pure
  f >>= g = shuffleTF g <*> f
  {-# INLINE return #-}
  {-# INLINE (>>=)  #-}

newtype TF_fmap f a   xs = TF_fmap (Fn f xs a)
data    TF_ap   f a b xs = TF_ap   (Fn f xs a) (Fn f xs b)

----------------------------------------------------------------
-- Operations on Fun
----------------------------------------------------------------

-- | Apply single parameter to function
curryFun :: Fun (x : xs) r -> x -> Fun xs r
curryFun = coerce
{-# INLINE curryFun #-}

-- | Uncurry N-ary function.
uncurryFun :: (x -> Fun xs r) -> Fun (x : xs) r
uncurryFun = coerce
{-# INLINE uncurryFun #-}

-- | Conversion function
uncurryMany :: forall xs ys r. Arity xs => Fun xs (Fun ys r) -> Fun (xs ++ ys) r
-- NOTE: GHC is not smart enough to figure out that:
--
--       > Fn xs (Fn ys) r ~ Fn (xs ++ ys) r
--
--       It's possible to construct type safe definition but it's
--       quite complicated and increase compile time and may hurrt
--       performance
{-# INLINE uncurryMany #-}
uncurryMany = unsafeCoerce

-- | Curry first /n/ arguments of N-ary function.
curryMany :: forall xs ys r. Arity xs => Fun (xs ++ ys) r -> Fun xs (Fun ys r)
-- NOTE: See uncurryMany
{-# INLINE curryMany #-}
curryMany = unsafeCoerce

-- | Add one parameter to function which is ignored.
constFun :: Fun xs r -> Fun (x : xs) r
constFun = uncurryFun . const
{-# INLINE constFun #-}

-- | Add one parameter to function which is ignored.
constTFun :: TFun f xs r -> TFun f (x : xs) r
constTFun = uncurryTFun . const
{-# INLINE constTFun #-}

-- | Transform function but leave outermost parameter untouched.
stepTFun :: (TFun f xs a       -> TFun f ys b)
         -> (TFun f (x : xs) a -> TFun f (x : ys) b)
stepTFun g = uncurryTFun . fmap g . curryTFun
{-# INLINE stepTFun #-}

-- | Concatenate n-ary functions. This function combine results of
--   both N-ary functions and merge their parameters into single list.
concatF :: (Arity xs, Arity ys)
        => (a -> b -> c) -> Fun xs a -> Fun ys b -> Fun (xs ++ ys) c
{-# INLINE concatF #-}
concatF f funA funB = uncurryMany $ fmap go funA
  where
    go a = fmap (\b -> f a b) funB

-- | Helper for lens implementation.
lensWorkerF :: forall f r x y xs. (Functor f, Arity xs)
            => (x -> f y) -> Fun (y : xs) r -> Fun (x : xs) (f r)
{-# INLINE lensWorkerF #-}
lensWorkerF g f
  = uncurryFun
  $ \x -> (\r -> fmap (r $) (g x)) <$> shuffleTF (curryFun f)

-- | Helper for lens implementation.
lensWorkerTF :: forall f g r x y xs. (Functor f, Arity xs)
             => (g x -> f (g y))
             -> TFun g (y : xs) r
             -> TFun g (x : xs) (f r)
{-# INLINE lensWorkerTF #-}
lensWorkerTF g f
  = uncurryTFun
  $ \x -> (\r -> fmap (r $) (g x)) <$> shuffleTF (curryTFun f)

----------------------------------------------------------------
-- Operations on TFun
----------------------------------------------------------------

-- | Apply single parameter to function
curryTFun :: TFun f (x : xs) r -> f x -> TFun f xs r
curryTFun = coerce
{-# INLINE curryTFun #-}

-- | Uncurry single parameter
uncurryTFun :: (f x -> TFun f xs r) -> TFun f (x : xs) r
uncurryTFun = coerce
{-# INLINE uncurryTFun #-}

-- | Move first argument of function to its result. This function is
--   useful for implementation of lens.
shuffleTF :: forall f x xs r. Arity xs
          => (x -> TFun f xs r) -> TFun f xs (x -> r)
{-# INLINE shuffleTF #-}
shuffleTF fun0 = accum
  (\(TF_shuffle f) a -> TF_shuffle (\x -> f x a))
  (\(TF_shuffle f)   -> f)
  (TF_shuffle (fmap unTFun fun0))

data TF_shuffle f x r xs = TF_shuffle (x -> Fn f xs r)

----------------------------------------------------------------
-- Indexing
----------------------------------------------------------------

-- | Indexing of vectors
class ArityPeano n => Index (n :: PeanoNum) (xs :: [*]) where
  -- | Type at position n
  type ValueAt n xs :: *
  -- | List of types with n'th element replaced by /a/.
  type NewElems n xs a :: [*]
  -- | Getter function for vectors
  getF :: proxy n -> Fun xs (ValueAt n xs)
  -- | Putter function. It applies value @x@ to @n@th parameter of
  --   function.
  putF :: proxy n -> ValueAt n xs -> Fun xs r -> Fun xs r
  -- | Helper for implementation of lens
  lensF   :: (Functor f, v ~ ValueAt n xs)
          => proxy n -> (v -> f v) -> Fun xs r -> Fun xs (f r)
  -- | Helper for type-changing lens
  lensChF :: (Functor f)
          => proxy n -> (ValueAt n xs -> f a) -> Fun (NewElems n xs a) r -> Fun xs (f r)

instance Arity xs => Index 'Z (x : xs) where
  type ValueAt  'Z (x : xs)   = x
  type NewElems 'Z (x : xs) a = a : xs
  getF  _     = TFun $ \(Identity x) -> unTFun (pure x :: Fun xs x)
  putF  _ x f = constFun $ curryFun f x
  lensF   _     = lensWorkerF
  lensChF _     = lensWorkerF
  {-# INLINE getF    #-}
  {-# INLINE putF    #-}
  {-# INLINE lensF   #-}
  {-# INLINE lensChF #-}

instance Index n xs => Index ('S n) (x : xs) where
  type ValueAt  ('S n) (x : xs)   = ValueAt n xs
  type NewElems ('S n) (x : xs) a = x : NewElems n xs a
  getF    _   = constFun $ getF    (Proxy @n)
  putF    _ x = stepTFun $ putF    (Proxy @n) x
  lensF   _ f = stepTFun $ lensF   (Proxy @n) f
  lensChF _ f = stepTFun $ lensChF (Proxy @n) f
  {-# INLINE getF    #-}
  {-# INLINE putF    #-}
  {-# INLINE lensF   #-}
  {-# INLINE lensChF #-}

----------------------------------------------------------------
-- Type lookup
----------------------------------------------------------------

-- | Type class to supporty looking up value in product type by its
--   type. Latter must not contain two elements of type @x@.
class Arity xs => TyLookup x xs where
  lookupTFun :: TFun f xs (f x)

-- Case analysis for type equality
class Arity xs => TyLookupCase (eq :: Bool) x xs where
  lookupTFunCase :: Proxy# eq -> TFun f xs (f x)

-- List xs does not contain type x
class NoType                  x xs
class NoTypeCase (eq :: Bool) x xs
instance                             NoType a '[]
instance NoTypeCase (a == x) a xs => NoType a (x ': xs)
instance ( TypeError ('Text "Duplicate type found: " ':$$: 'ShowType a)
         )           => NoTypeCase 'True  a xs
instance NoType a xs => NoTypeCase 'False a xs

instance ( TypeError ('Text "Cannot find type: " ':$$: 'ShowType a)
         ) => TyLookup a '[] where
  lookupTFun = error "Unreachable"

-- Case analysis of type equality
instance ( Arity xs
         , TyLookupCase (a == x) a (x ': xs)
         ) => TyLookup a (x ': xs) where
  lookupTFun = lookupTFunCase (proxy# :: Proxy# (a == x))
  {-# INLINE lookupTFun #-}

-- Found x
instance ( Arity xs
         , NoType x xs
         ) => TyLookupCase 'True x (x ': xs) where
  lookupTFunCase _ = uncurryTFun pure
  {-# INLINE lookupTFunCase #-}

-- Go deeper
instance ( Arity xs
         , TyLookup a xs
         ) => TyLookupCase 'False a (x ': xs) where
  lookupTFunCase _ = uncurryTFun $ const lookupTFun
  {-# INLINE lookupTFunCase #-}

----------------------------------------------------------------
-- Instances
----------------------------------------------------------------

-- | Unit is empty heterogeneous vector
instance HVector () where
  type Elems () = '[]
  construct = TFun ()
  inspect () (TFun f) = f

instance HVector (Complex a) where
  type Elems (Complex a) = '[a,a]
  construct = TFun $ \(Identity r) (Identity i) -> (:+) r i
  inspect (r :+ i) f = coerce f r i
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b) where
  type Elems (a,b) = '[a,b]
  construct = coerce ((,) :: a->b -> (a,b))
  inspect (a,b) f = coerce f a b
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c) where
  type Elems (a,b,c) = '[a,b,c]
  construct = coerce ((,,) :: a->b->c -> (a,b,c))
  inspect (a,b,c) f = coerce f a b c
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d) where
  type Elems (a,b,c,d) = '[a,b,c,d]
  construct = coerce ((,,,) :: a->b->c->d -> (a,b,c,d))
  inspect (a,b,c,d) f = coerce f a b c d
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e) where
  type Elems (a,b,c,d,e) = '[a,b,c,d,e]
  construct = coerce ((,,,,) :: a->b->c->d->e -> (a,b,c,d,e))
  inspect (a,b,c,d,e) f = coerce f a b c d e
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f) where
  type Elems (a,b,c,d,e,f) = '[a,b,c,d,e,f]
  construct = coerce ((,,,,,) :: a->b->c->d->e->f
                              -> (a,b,c,d,e,f))
  inspect (a,b,c,d,e,f) fun = coerce fun a b c d e f
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g) where
  type Elems (a,b,c,d,e,f,g) = '[a,b,c,d,e,f,g]
  construct = coerce ((,,,,,,) :: a->b->c->d->e->f->g
                               -> (a,b,c,d,e,f,g))
  inspect (a,b,c,d,e,f,g) fun = coerce fun a b c d e f g
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h) where
  type Elems (a,b,c,d,e,f,g,h) = '[a,b,c,d,e,f,g,h]
  construct = coerce ((,,,,,,,) :: a->b->c->d->e->f->g->h
                                -> (a,b,c,d,e,f,g,h))
  inspect (a,b,c,d,e,f,g,h) fun = coerce fun a b c d e f g h
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i) where
  type Elems (a,b,c,d,e,f,g,h,i) = '[a,b,c,d,e,f,g,h,i]
  construct = coerce ((,,,,,,,,) :: a->b->c->d->e->f->g->h->i
                                 -> (a,b,c,d,e,f,g,h,i))
  inspect (a,b,c,d,e,f,g,h,i) fun = coerce fun a b c d e f g h i
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j) where
  type Elems (a,b,c,d,e,f,g,h,i,j) = '[a,b,c,d,e,f,g,h,i,j]
  construct = coerce ((,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j
                                  -> (a,b,c,d,e,f,g,h,i,j))
  inspect (a,b,c,d,e,f,g,h,i,j) fun = coerce fun a b c d e f g h i j
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k) where
  type Elems (a,b,c,d,e,f,g,h,i,j,k) = '[a,b,c,d,e,f,g,h,i,j,k]
  construct = coerce ((,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k
                                   -> (a,b,c,d,e,f,g,h,i,j,k))
  inspect (a,b,c,d,e,f,g,h,i,j,k) fun = coerce fun a b c d e f g h i j k
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k,l) where
  type Elems (a,b,c,d,e,f,g,h,i,j,k,l) = '[a,b,c,d,e,f,g,h,i,j,k,l]
  construct = coerce ((,,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k->l
                                    -> (a,b,c,d,e,f,g,h,i,j,k,l))
  inspect (a,b,c,d,e,f,g,h,i,j,k,l) fun = coerce fun a b c d e f g h i j k l
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k,l,m) where
  type Elems (a,b,c,d,e,f,g,h,i,j,k,l,m) = '[a,b,c,d,e,f,g,h,i,j,k,l,m]
  construct = coerce ((,,,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k->l->m
                                     -> (a,b,c,d,e,f,g,h,i,j,k,l,m))
  inspect (a,b,c,d,e,f,g,h,i,j,k,l,m) fun = coerce fun a b c d e f g h i j k l m
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k,l,m,n) where
  type Elems (a,b,c,d,e,f,g,h,i,j,k,l,m,n) = '[a,b,c,d,e,f,g,h,i,j,k,l,m,n]
  construct = coerce ((,,,,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k->l->m->n
                                      -> (a,b,c,d,e,f,g,h,i,j,k,l,m,n))
  inspect (a,b,c,d,e,f,g,h,i,j,k,l,m,n) fun
    = coerce fun a b c d e f g h i j k l m n
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o) where
  type Elems (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o) = '[a,b,c,d,e,f,g,h,i,j,k,l,m,n,o]
  construct = coerce ((,,,,,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k->l->m->n->o
                                       -> (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o))
  inspect (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o) fun
    = coerce fun a b c d e f g h i j k l m n o
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p) where
  type Elems (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p) =
    '[a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p]
  construct = coerce ((,,,,,,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k->l->m->n->o->p
                                        -> (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p))
  inspect (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p) fun
    = coerce fun a b c d e f g h i j k l m n o p
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q) where
  type Elems (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q) =
    '[a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q]
  construct = coerce ((,,,,,,,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k->l->m->n->o->p->q
                                         -> (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q))
  inspect (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q) fun
    = coerce fun a b c d e f g h i j k l m n o p q
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r) where
  type Elems (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r) =
    '[a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r]
  construct = coerce ((,,,,,,,,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k->l->m->n->o->p->q->r
                                          -> (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r))
  inspect (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r) fun
    = coerce fun a b c d e f g h i j k l m n o p q r
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s) where
  type Elems (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s) =
    '[a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s]
  construct = coerce ((,,,,,,,,,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k->l->m->n->o->p->q->r->s
                                           -> (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s))
  inspect (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s) fun
    = coerce fun a b c d e f g h i j k l m n o p q r s
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t) where
  type Elems (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t) =
    '[a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t]
  construct = coerce ((,,,,,,,,,,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k->l->m->n->o->p->q->r->s->t
                                            -> (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t))
  inspect (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t) fun
    = coerce fun a b c d e f g h i j k l m n o p q r s t
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u) where
  type Elems (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u) =
    '[a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u]
  construct = coerce ((,,,,,,,,,,,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k->l->m->n->o->p->q->r->s->t->u
                                             -> (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u))
  inspect (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u) fun
    = coerce fun a b c d e f g h i j k l m n o p q r s t u
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v) where
  type Elems (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v) =
    '[a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v]
  construct = coerce ((,,,,,,,,,,,,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k->l->m->n->o->p->q->r->s->t->u->v
                                              -> (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v))
  inspect (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v) fun
    = coerce fun a b c d e f g h i j k l m n o p q r s t u v
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w) where
  type Elems (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w) =
    '[a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w]
  construct = coerce ((,,,,,,,,,,,,,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k->l->m->n->o->p->q->r->s->t->u->v->w
                                               -> (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w))
  inspect (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w) fun
    = coerce fun a b c d e f g h i j k l m n o p q r s t u v w
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x) where
  type Elems (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x) =
    '[a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x]
  construct = coerce ((,,,,,,,,,,,,,,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k->l->m->n->o->p->q->r->s->t->u->v->w->x
                                                -> (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x))
  inspect (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x) fun
    = coerce fun a b c d e f g h i j k l m n o p q r s t u v w x
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y) where
  type Elems (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y) =
    '[a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y]
  construct = coerce ((,,,,,,,,,,,,,,,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k->l->m->n->o->p->q->r->s->t->u->v->w->x->y
                                                 -> (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y))
  inspect (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y) fun
    = coerce fun a b c d e f g h i j k l m n o p q r s t u v w x y
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z) where
  type Elems (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z) =
    '[a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z]
  construct = coerce ((,,,,,,,,,,,,,,,,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k->l->m->n->o->p->q->r->s->t->u->v->w->x->y->z
                                                  -> (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z))
  inspect (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z) fun
    = coerce fun a b c d e f g h i j k l m n o p q r s t u v w x y z
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a') where
  type Elems (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a') =
    '[a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a']
  construct = coerce ((,,,,,,,,,,,,,,,,,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k->l->m->n->o->p->q->r->s->t->u->v->w->x->y->z->a'
                                                   -> (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a'))
  inspect (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a') fun
    = coerce fun a b c d e f g h i j k l m n o p q r s t u v w x y z a'
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b') where
  type Elems (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b') =
    '[a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b']
  construct = coerce ((,,,,,,,,,,,,,,,,,,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k->l->m->n->o->p->q->r->s->t->u->v->w->x->y->z->a'->b'
                                                    -> (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b'))
  inspect (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b') fun
    = coerce fun a b c d e f g h i j k l m n o p q r s t u v w x y z a' b'
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b',c') where
  type Elems (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b',c') =
    '[a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b',c']
  construct = coerce ((,,,,,,,,,,,,,,,,,,,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k->l->m->n->o->p->q->r->s->t->u->v->w->x->y->z->a'->b'->c'
                                                     -> (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b',c'))
  inspect (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b',c') fun
    = coerce fun a b c d e f g h i j k l m n o p q r s t u v w x y z a' b' c'
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b',c',d') where
  type Elems (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b',c',d') =
    '[a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b',c',d']
  construct = coerce ((,,,,,,,,,,,,,,,,,,,,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k->l->m->n->o->p->q->r->s->t->u->v->w->x->y->z->a'->b'->c'->d'
                                                      -> (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b',c',d'))
  inspect (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b',c',d') fun
    = coerce fun a b c d e f g h i j k l m n o p q r s t u v w x y z a' b' c' d'
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b',c',d',e') where
  type Elems (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b',c',d',e')
    = '[a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b',c',d',e']
  construct = coerce ((,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k->l->m->n->o->p->q->r->s->t->u->v->w->x->y->z->a'->b'->c'->d'->e'
                                                       -> (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b',c',d',e'))
  inspect (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b',c',d',e') fun
    = coerce fun a b c d e f g h i j k l m n o p q r s t u v w x y z a' b' c' d' e'
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

instance HVector (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b',c',d',e',f') where
  type Elems (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b',c',d',e',f')
    = '[a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b',c',d',e',f']
  construct = coerce ((,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,) :: a->b->c->d->e->f->g->h->i->j->k->l->m->n->o->p->q->r->s->t->u->v->w->x->y->z->a'->b'->c'->d'->e'->f'
                                                        -> (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b',c',d',e',f'))
  inspect (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,a',b',c',d',e',f') fun
    = coerce fun a b c d e f g h i j k l m n o p q r s t u v w x y z a' b' c' d' e' f'
  {-# INLINE construct #-}
  {-# INLINE inspect   #-}

-- | Copy of lens type definition from lens package
type Lens s t a b = forall f. Functor f => (a -> f b) -> s -> f t

-- | Copy of type preserving lens definition from lens package
type Lens' s a = Lens s s a a

----------------------------------------------------------------
-- Generics
----------------------------------------------------------------

class GHVector (v :: * -> *) where
  type GElems v :: [*]
  gconstruct :: Fun (GElems v) (v p)
  ginspect   :: v p -> Fun (GElems v) r -> r

-- We simply skip metadata
instance (GHVector f, Arity (GElems f)) => GHVector (M1 i c f) where
  type GElems (M1 i c f) = GElems f
  gconstruct = fmap M1 gconstruct
  ginspect v = ginspect (unM1 v)
  {-# INLINE gconstruct #-}
  {-# INLINE ginspect   #-}

instance ( GHVector f, GHVector g, Arity (GElems f), Arity (GElems g)
         ) => GHVector (f :*: g) where
  type GElems (f :*: g) = GElems f ++ GElems g
  gconstruct = concatF (:*:) gconstruct gconstruct
  ginspect (f :*: g) fun
    = ginspect g $ ginspect f $ curryMany fun
  {-# INLINE gconstruct #-}
  {-# INLINE ginspect   #-}

instance ( TypeError ('Text "It's impossible to derive HVector for type without constructors")
         ) => GHVector V1 where
  type GElems V1 = TypeError ('Text "It's impossible to derive HVector for type without constructors")
  gconstruct = error "Unreachable"
  ginspect   = error "Unreachable"

instance ( TypeError ('Text "It's impossible to derive HVector for sum types")
         ) => GHVector (f :+: g) where
  type GElems (f :+: g) = TypeError ('Text "It's impossible to derive HVector for sum types")
  gconstruct = error "Unreachable"
  ginspect   = error "Unreachable"

-- Recursion is terminated by simple field
instance GHVector (K1 R x) where
  type GElems (K1 R x) = '[x]
  gconstruct               = TFun (K1 . runIdentity)
  ginspect (K1 x) (TFun f) = f (Identity x)
  {-# INLINE gconstruct #-}
  {-# INLINE ginspect   #-}

-- Unit types are empty vectors
instance GHVector U1 where
  type GElems U1      = '[]
  gconstruct          = coerce U1
  ginspect _ (TFun f) = f
  {-# INLINE gconstruct #-}
  {-# INLINE ginspect   #-}

max_stars_repo_name: Shimuuar/fixed-vector-hetero

max_issues_repo_name: Shimuuar/fixed-vector-hetero

max_forks_repo_name: Shimuuar/fixed-vector-hetero


text:
{-# LANGUAGE FlexibleContexts       #-}
{-# LANGUAGE FlexibleInstances      #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE MultiParamTypeClasses  #-}
{-# LANGUAGE TypeFamilies           #-}
{-# LANGUAGE UndecidableInstances   #-}

-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Conversion
-- Copyright   :  (c) Alberto Ruiz 2010
-- License     :  BSD3
-- Maintainer  :  Alberto Ruiz
-- Stability   :  provisional
--
-- Conversion routines
--
-----------------------------------------------------------------------------

module Internal.Conversion (
    Complexable(..), RealElement,
    module Data.Complex
) where

import           Control.Arrow       ((***))
import           Data.Complex
import           Internal.Matrix
import           Internal.Vector
import           Internal.Vectorized

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

-- | Supported single-double precision type pairs
class (Element s, Element d) => Precision s d | s -> d, d -> s where
    double2FloatG :: Vector d -> Vector s
    float2DoubleG :: Vector s -> Vector d

-- instance Precision Float Double where
--     double2FloatG = double2FloatV
--     float2DoubleG = float2DoubleV

-- instance Precision (Complex Float) (Complex Double) where
--     double2FloatG = asComplex . double2FloatV . asReal
--     float2DoubleG = asComplex . float2DoubleV . asReal

instance Precision I Z where
    double2FloatG = long2intV
    float2DoubleG = int2longV

-- | Supported real types
class (Element t, Element (Complex t), RealFloat t)
    => RealElement t

-- instance RealElement Double
instance RealElement Float

-- | Structures that may contain complex numbers
class Complexable c where
    toComplex'   :: (RealElement e) => (c e, c e) -> c (Complex e)
    fromComplex' :: (RealElement e) => c (Complex e) -> (c e, c e)
    comp'        :: (RealElement e) => c e -> c (Complex e)
    single'      :: Precision a b => c b -> c a
    double'      :: Precision a b => c a -> c b

instance Complexable Vector where
    toComplex' = toComplexV
    fromComplex' = fromComplexV
    comp' v = toComplex' (v,constantD 0 (dim v))
    single' = double2FloatG
    double' = float2DoubleG

-- | creates a complex vector from vectors with real and imaginary parts
toComplexV :: (RealElement a) => (Vector a, Vector a) ->  Vector (Complex a)
toComplexV (r,i) = asComplex $ flatten $ fromColumns [r,i]

-- | the inverse of 'toComplex'
fromComplexV :: (RealElement a) => Vector (Complex a) -> (Vector a, Vector a)
fromComplexV z = (r,i) where
    [r,i] = toColumns $ reshape 2 $ asReal z

instance Complexable Matrix where
    toComplex' = uncurry $ liftMatrix2 $ curry toComplex'
    fromComplex' z = (reshape c *** reshape c) . fromComplex' . flatten $ z
        where c = cols z
    comp' = liftMatrix comp'
    single' = liftMatrix single'
    double' = liftMatrix double'

ext: hs

max_stars_repo_name: schnecki/hmatrix-float

max_issues_repo_name: schnecki/hmatrix-float

max_forks_repo_name: schnecki/hmatrix-float


text:
module Imports where

import           Control.Applicative
import           Control.Arrow
import           Control.Monad
import           Control.Monad.Cont
import           Control.Monad.Error
import           Control.Monad.Fix
import           Control.Monad.Identity
import           Control.Monad.Loops
import           Control.Monad.RWS
import           Control.Monad.Reader
import           Control.Monad.State
import           Control.Monad.Writer
import qualified Data.Aeson as A
import           Data.Array
import           Data.Bits
import           Data.Bool
import qualified Data.ByteString as B
import qualified Data.ByteString.Lazy as BL
import           Data.Char
import           Data.Complex
import           Data.Dynamic
import           Data.Either
import           Data.Eq
import           Data.Fixed
import           Data.Function
import           Data.Graph
import           Data.Int
import           Data.Ix
import           Data.List
import           Data.List.Split
import qualified Data.Map as M
import           Data.Maybe
import           Data.Monoid
import           Data.Ord
import           Data.Ratio
import qualified Data.Set as S
import           Data.Text (Text)
import qualified Data.Text as T
import qualified Data.Text.Encoding as T
import qualified Data.Text.Lazy as TL
import qualified Data.Text.Lazy.Encoding as TL
import           Data.Tree
import           Data.Tuple
import           Data.Typeable
import           Data.Word
import           Debug.SimpleReflect
import           Prelude hiding (IO,putStr,putStrLn,getLine,readLn,print,readIO,readFile,writeFile,appendFile)
import           Safe hiding (abort)
import           ShowFun
import           System.Random
import           Text.Printf

-- Convenience for golfing.
asS :: (String -> String) -> Text -> Text
asS f = T.pack . f . T.unpack

asS' = flip asS

max_stars_repo_name: robinp/repli

max_issues_repo_name: robinp/repli

max_forks_repo_name: robinp/repli


text:
{-# OPTIONS_GHC -Wall -Wno-partial-type-signatures #-}

{-# LANGUAGE QuasiQuotes #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE MultiWayIf #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE PartialTypeSignatures #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}

-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Sundials.ARKode.ODE
-- Copyright   :  Dominic Steinitz 2018,
--                Novadiscovery 2018
-- License     :  BSD
-- Maintainer  :  Dominic Steinitz
-- Stability   :  provisional
--
-- Solution of ordinary differential equation (ODE) initial value problems.
-- See <https://computation.llnl.gov/projects/sundials/sundials-software> for more detail.
--
-- A simple example:
--
-- <<diagrams/brusselator.png#diagram=brusselator&height=400&width=500>>
--
-- @
-- import           Numeric.Sundials.ARKode.ODE
-- import           Numeric.LinearAlgebra
--
-- import           Plots as P
-- import qualified Diagrams.Prelude as D
-- import           Diagrams.Backend.Rasterific
--
-- brusselator :: Double -> [Double] -> [Double]
-- brusselator _t x = [ a - (w + 1) * u + v * u * u
--                    , w * u - v * u * u
--                    , (b - w) / eps - w * u
--                    ]
--   where
--     a = 1.0
--     b = 3.5
--     eps = 5.0e-6
--     u = x !! 0
--     v = x !! 1
--     w = x !! 2
--
-- lSaxis :: [[Double]] -> P.Axis B D.V2 Double
-- lSaxis xs = P.r2Axis &~ do
--   let ts = xs!!0
--       us = xs!!1
--       vs = xs!!2
--       ws = xs!!3
--   P.linePlot' $ zip ts us
--   P.linePlot' $ zip ts vs
--   P.linePlot' $ zip ts ws
--
-- main = do
--   let res1 = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
--   renderRasterific "diagrams/brusselator.png"
--                    (D.dims2D 500.0 500.0)
--                    (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1))
-- @
--
-- With Sundials ARKode, it is possible to retrieve the Butcher
-- tableau for the solver. FIXME: Not available just now and hopefully
-- normal service will be resumed soon.
--
-- @
-- import           Numeric.Sundials.ARKode.ODE
-- import           Numeric.LinearAlgebra
--
-- import           Data.List (intercalate)
--
-- import           Text.PrettyPrint.HughesPJClass
--
--
-- butcherTableauTex :: ButcherTable -> String
-- butcherTableauTex (ButcherTable m c b b2) =
--   render $
--   vcat [ text ("\n\\begin{array}{c|" ++ (concat $ replicate n "c") ++ "}")
--        , us
--        , text "\\hline"
--        , text bs <+> text "\\\\"
--        , text b2s <+> text "\\\\"
--        , text "\\end{array}"
--        ]
--   where
--     n = rows m
--     rs = toLists m
--     ss = map (\r -> intercalate " & " $ map show r) rs
--     ts = zipWith (\i r -> show i ++ " & " ++ r) (toList c) ss
--     us = vcat $ map (\r -> text r <+> text "\\\\") ts
--     bs  = " & " ++ (intercalate " & " $ map show $ toList b)
--     b2s = " & " ++ (intercalate " & " $ map show $ toList b2)
--
-- main :: IO ()
-- main = do
--
--   let res = butcherTable (SDIRK_2_1_2 undefined)
--   putStrLn $ show res
--   putStrLn $ butcherTableauTex res
--
--   let resA = butcherTable (KVAERNO_4_2_3 undefined)
--   putStrLn $ show resA
--   putStrLn $ butcherTableauTex resA
--
--   let resB = butcherTable (SDIRK_5_3_4 undefined)
--   putStrLn $ show resB
--   putStrLn $ butcherTableauTex resB
-- @
--
--  Using the code above from the examples gives
--
-- KVAERNO_4_2_3
--
-- \[
-- \begin{array}{c|cccc}
-- 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\
-- 0.871733043 & 0.4358665215 & 0.4358665215 & 0.0 & 0.0 \\
-- 1.0 & 0.490563388419108 & 7.3570090080892e-2 & 0.4358665215 & 0.0 \\
-- 1.0 & 0.308809969973036 & 1.490563388254106 & -1.235239879727145 & 0.4358665215 \\
-- \hline
--  & 0.308809969973036 & 1.490563388254106 & -1.235239879727145 & 0.4358665215 \\
--  & 0.490563388419108 & 7.3570090080892e-2 & 0.4358665215 & 0.0 \\
-- \end{array}
-- \]
--
-- SDIRK_2_1_2
--
-- \[
-- \begin{array}{c|cc}
-- 1.0 & 1.0 & 0.0 \\
-- 0.0 & -1.0 & 1.0 \\
-- \hline
--  & 0.5 & 0.5 \\
--  & 1.0 & 0.0 \\
-- \end{array}
-- \]
--
-- SDIRK_5_3_4
--
-- \[
-- \begin{array}{c|ccccc}
-- 0.25 & 0.25 & 0.0 & 0.0 & 0.0 & 0.0 \\
-- 0.75 & 0.5 & 0.25 & 0.0 & 0.0 & 0.0 \\
-- 0.55 & 0.34 & -4.0e-2 & 0.25 & 0.0 & 0.0 \\
-- 0.5 & 0.2727941176470588 & -5.036764705882353e-2 & 2.7573529411764705e-2 & 0.25 & 0.0 \\
-- 1.0 & 1.0416666666666667 & -1.0208333333333333 & 7.8125 & -7.083333333333333 & 0.25 \\
-- \hline
--  & 1.0416666666666667 & -1.0208333333333333 & 7.8125 & -7.083333333333333 & 0.25 \\
--  & 1.2291666666666667 & -0.17708333333333334 & 7.03125 & -7.083333333333333 & 0.0 \\
-- \end{array}
-- \]
-----------------------------------------------------------------------------
module Numeric.Sundials.ARKode.ODE ( odeSolve
                                   , odeSolveV
                                   , odeSolveVWith
                                   , odeSolveVWith'
                                   , odeSolveWithEvents
                                   , ODEMethod(..)
                                   , StepControl(..)
                                   ) where

import qualified Language.C.Inline as C
import qualified Language.C.Inline.Unsafe as CU

import           Data.Monoid ((<>))
import           Data.Maybe (isJust)

import           Foreign.C.Types (CDouble, CInt)
import           Foreign.Ptr (Ptr)
import           Foreign.Storable (poke, peek)

import qualified Data.Vector.Storable as V

import           Data.Coerce (coerce)
import           System.IO.Unsafe (unsafePerformIO)
import           GHC.Generics (C1, Constructor, (:+:)(..), D1, Rep, Generic, M1(..),
                               from, conName)

import           Numeric.LinearAlgebra.Devel (createVector)

import           Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, rows,
                                                cols, toLists, size, reshape,
                                                (><))

import           Numeric.Sundials.Types
import qualified Numeric.Sundials.Arkode as T
import           Numeric.Sundials.Arkode (SunIndexType)
import           Numeric.Sundials.Arkode (sDIRK_2_1_2,
                                          bILLINGTON_3_3_2,
                                          tRBDF2_3_3_2,
                                          kVAERNO_4_2_3,
                                          aRK324L2SA_DIRK_4_2_3,
                                          cASH_5_2_4,
                                          cASH_5_3_4,
                                          sDIRK_5_3_4,
                                          kVAERNO_5_3_4,
                                          aRK436L2SA_DIRK_6_3_4,
                                          kVAERNO_7_4_5,
                                          aRK548L2SA_DIRK_8_4_5,
                                          hEUN_EULER_2_1_2,
                                          bOGACKI_SHAMPINE_4_2_3,
                                          aRK324L2SA_ERK_4_2_3,
                                          zONNEVELD_5_3_4,
                                          aRK436L2SA_ERK_6_3_4,
                                          sAYFY_ABURUB_6_3_4,
                                          cASH_KARP_6_4_5,
                                          fEHLBERG_6_4_5,
                                          dORMAND_PRINCE_7_4_5,
                                          aRK548L2SA_ERK_8_4_5,
                                          vERNER_8_5_6,
                                          fEHLBERG_13_7_8)

C.context (C.baseCtx <> C.vecCtx <> C.funCtx <> sunCtx)

C.include "<stdlib.h>"
C.include "<stdio.h>"
C.include "<math.h>"
C.include "<arkode/arkode_arkstep.h>"                 -- prototypes for ARKODE fcts., consts.
C.include "<arkode/arkode_erkstep.h>"
C.include "<nvector/nvector_serial.h>"        -- serial N_Vector types, fcts., macros
C.include "<sunmatrix/sunmatrix_dense.h>"     -- access to dense SUNMatrix
C.include "<sunlinsol/sunlinsol_dense.h>"     -- access to dense SUNLinearSolver
C.include "<sundials/sundials_types.h>"       -- definition of type realtype
C.include "<sundials/sundials_math.h>"
C.include "../../../helpers.h"
C.include "Numeric/Sundials/Arkode_hsc.h"

-- | Stepping functions
data ODEMethod = SDIRK_2_1_2            Jacobian
               | SDIRK_2_1_2'
               | BILLINGTON_3_3_2       Jacobian
               | BILLINGTON_3_3_2'
               | TRBDF2_3_3_2           Jacobian
               | TRBDF2_3_3_2'
               | KVAERNO_4_2_3          Jacobian
               | KVAERNO_4_2_3'
               | ARK324L2SA_DIRK_4_2_3  Jacobian
               | ARK324L2SA_DIRK_4_2_3'
               | CASH_5_2_4             Jacobian
               | CASH_5_2_4'
               | CASH_5_3_4             Jacobian
               | CASH_5_3_4'
               | SDIRK_5_3_4            Jacobian
               | SDIRK_5_3_4'
               | KVAERNO_5_3_4          Jacobian
               | KVAERNO_5_3_4'
               | ARK436L2SA_DIRK_6_3_4  Jacobian
               | ARK436L2SA_DIRK_6_3_4'
               | KVAERNO_7_4_5          Jacobian
               | KVAERNO_7_4_5'
               | ARK548L2SA_DIRK_8_4_5  Jacobian
               | ARK548L2SA_DIRK_8_4_5'
               | HEUN_EULER_2_1_2         Jacobian
               | HEUN_EULER_2_1_2'
               | BOGACKI_SHAMPINE_4_2_3   Jacobian
               | BOGACKI_SHAMPINE_4_2_3'
               | ARK324L2SA_ERK_4_2_3     Jacobian
               | ARK324L2SA_ERK_4_2_3'
               | ZONNEVELD_5_3_4          Jacobian
               | ZONNEVELD_5_3_4'
               | ARK436L2SA_ERK_6_3_4     Jacobian
               | ARK436L2SA_ERK_6_3_4'
               | SAYFY_ABURUB_6_3_4       Jacobian
               | SAYFY_ABURUB_6_3_4'
               | CASH_KARP_6_4_5          Jacobian
               | CASH_KARP_6_4_5'
               | FEHLBERG_6_4_5         Jacobian
               | FEHLBERG_6_4_5'
               | DORMAND_PRINCE_7_4_5     Jacobian
               | DORMAND_PRINCE_7_4_5'
               | ARK548L2SA_ERK_8_4_5     Jacobian
               | ARK548L2SA_ERK_8_4_5'
               | VERNER_8_5_6            Jacobian
               | VERNER_8_5_6'
               | FEHLBERG_13_7_8         Jacobian
               | FEHLBERG_13_7_8'
  deriving Generic

constrName :: (HasConstructor (Rep a), Generic a)=> a -> String
constrName = genericConstrName . from

class HasConstructor (f :: * -> *) where
  genericConstrName :: f x -> String

instance HasConstructor f => HasConstructor (D1 c f) where
  genericConstrName (M1 x) = genericConstrName x

instance (HasConstructor x, HasConstructor y) => HasConstructor (x :+: y) where
  genericConstrName (L1 l) = genericConstrName l
  genericConstrName (R1 r) = genericConstrName r

instance Constructor c => HasConstructor (C1 c f) where
  genericConstrName x = conName x

instance Show ODEMethod where
  show x = constrName x

-- FIXME: We can probably do better here with generics
getMethod :: ODEMethod -> Int
getMethod (SDIRK_2_1_2 _)            = sDIRK_2_1_2
getMethod (SDIRK_2_1_2')             = sDIRK_2_1_2
getMethod (BILLINGTON_3_3_2 _)       = bILLINGTON_3_3_2
getMethod (BILLINGTON_3_3_2')        = bILLINGTON_3_3_2
getMethod (TRBDF2_3_3_2 _)           = tRBDF2_3_3_2
getMethod (TRBDF2_3_3_2')            = tRBDF2_3_3_2
getMethod (KVAERNO_4_2_3  _)         = kVAERNO_4_2_3
getMethod (KVAERNO_4_2_3')           = kVAERNO_4_2_3
getMethod (ARK324L2SA_DIRK_4_2_3 _)  = aRK324L2SA_DIRK_4_2_3
getMethod (ARK324L2SA_DIRK_4_2_3')   = aRK324L2SA_DIRK_4_2_3
getMethod (CASH_5_2_4 _)             = cASH_5_2_4
getMethod (CASH_5_2_4')              = cASH_5_2_4
getMethod (CASH_5_3_4 _)             = cASH_5_3_4
getMethod (CASH_5_3_4')              = cASH_5_3_4
getMethod (SDIRK_5_3_4 _)            = sDIRK_5_3_4
getMethod (SDIRK_5_3_4')             = sDIRK_5_3_4
getMethod (KVAERNO_5_3_4 _)          = kVAERNO_5_3_4
getMethod (KVAERNO_5_3_4')           = kVAERNO_5_3_4
getMethod (ARK436L2SA_DIRK_6_3_4 _)  = aRK436L2SA_DIRK_6_3_4
getMethod (ARK436L2SA_DIRK_6_3_4')   = aRK436L2SA_DIRK_6_3_4
getMethod (KVAERNO_7_4_5 _)          = kVAERNO_7_4_5
getMethod (KVAERNO_7_4_5')           = kVAERNO_7_4_5
getMethod (ARK548L2SA_DIRK_8_4_5 _)  = aRK548L2SA_DIRK_8_4_5
getMethod (ARK548L2SA_DIRK_8_4_5')   = aRK548L2SA_DIRK_8_4_5
getMethod (HEUN_EULER_2_1_2 _)       = hEUN_EULER_2_1_2
getMethod (HEUN_EULER_2_1_2')        = hEUN_EULER_2_1_2
getMethod (BOGACKI_SHAMPINE_4_2_3 _) = bOGACKI_SHAMPINE_4_2_3
getMethod (BOGACKI_SHAMPINE_4_2_3')  = bOGACKI_SHAMPINE_4_2_3
getMethod (ARK324L2SA_ERK_4_2_3 _)   = aRK324L2SA_ERK_4_2_3
getMethod (ARK324L2SA_ERK_4_2_3')    = aRK324L2SA_ERK_4_2_3
getMethod (ZONNEVELD_5_3_4 _)        = zONNEVELD_5_3_4
getMethod (ZONNEVELD_5_3_4')         = zONNEVELD_5_3_4
getMethod (ARK436L2SA_ERK_6_3_4 _)   = aRK436L2SA_ERK_6_3_4
getMethod (ARK436L2SA_ERK_6_3_4')    = aRK436L2SA_ERK_6_3_4
getMethod (SAYFY_ABURUB_6_3_4 _)     = sAYFY_ABURUB_6_3_4
getMethod (SAYFY_ABURUB_6_3_4')      = sAYFY_ABURUB_6_3_4
getMethod (CASH_KARP_6_4_5 _)        = cASH_KARP_6_4_5
getMethod (CASH_KARP_6_4_5')         = cASH_KARP_6_4_5
getMethod (FEHLBERG_6_4_5 _)         = fEHLBERG_6_4_5
getMethod (FEHLBERG_6_4_5' )         = fEHLBERG_6_4_5
getMethod (DORMAND_PRINCE_7_4_5 _)   = dORMAND_PRINCE_7_4_5
getMethod (DORMAND_PRINCE_7_4_5')    = dORMAND_PRINCE_7_4_5
getMethod (ARK548L2SA_ERK_8_4_5 _)   = aRK548L2SA_ERK_8_4_5
getMethod (ARK548L2SA_ERK_8_4_5')    = aRK548L2SA_ERK_8_4_5
getMethod (VERNER_8_5_6 _)           = vERNER_8_5_6
getMethod (VERNER_8_5_6')            = vERNER_8_5_6
getMethod (FEHLBERG_13_7_8 _)        = fEHLBERG_13_7_8
getMethod (FEHLBERG_13_7_8')         = fEHLBERG_13_7_8

getJacobian :: ODEMethod -> Maybe Jacobian
getJacobian (SDIRK_2_1_2 j)            = Just j
getJacobian (BILLINGTON_3_3_2 j)       = Just j
getJacobian (TRBDF2_3_3_2 j)           = Just j
getJacobian (KVAERNO_4_2_3  j)         = Just j
getJacobian (ARK324L2SA_DIRK_4_2_3 j)  = Just j
getJacobian (CASH_5_2_4 j)             = Just j
getJacobian (CASH_5_3_4 j)             = Just j
getJacobian (SDIRK_5_3_4 j)            = Just j
getJacobian (KVAERNO_5_3_4 j)          = Just j
getJacobian (ARK436L2SA_DIRK_6_3_4 j)  = Just j
getJacobian (KVAERNO_7_4_5 j)          = Just j
getJacobian (ARK548L2SA_DIRK_8_4_5 j)  = Just j
getJacobian (HEUN_EULER_2_1_2 j)       = Just j
getJacobian (BOGACKI_SHAMPINE_4_2_3 j) = Just j
getJacobian (ARK324L2SA_ERK_4_2_3 j)   = Just j
getJacobian (ZONNEVELD_5_3_4 j)        = Just j
getJacobian (ARK436L2SA_ERK_6_3_4 j)   = Just j
getJacobian (SAYFY_ABURUB_6_3_4 j)     = Just j
getJacobian (CASH_KARP_6_4_5 j)        = Just j
getJacobian (FEHLBERG_6_4_5 j)         = Just j
getJacobian (DORMAND_PRINCE_7_4_5 j)   = Just j
getJacobian (ARK548L2SA_ERK_8_4_5 j)   = Just j
getJacobian (VERNER_8_5_6 j)           = Just j
getJacobian (FEHLBERG_13_7_8 j)        = Just j
getJacobian _                          = Nothing

-- | A version of 'odeSolveVWith' with reasonable default step control.
odeSolveV
    :: ODEMethod
    -> Maybe Double      -- ^ initial step size - by default, ARKode
                         -- estimates the initial step size to be the
                         -- solution \(h\) of the equation
                         -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
                         -- \(\ddot{y}\) is an estimated value of the
                         -- second derivative of the solution at \(t_0\)
    -> Double            -- ^ absolute tolerance for the state vector
    -> Double            -- ^ relative tolerance for the state vector
    -> (Double -> Vector Double -> Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
    -> Vector Double     -- ^ initial conditions
    -> Vector Double     -- ^ desired solution times
    -> Matrix Double     -- ^ solution
odeSolveV meth hi epsAbs epsRel f y0 ts =
  odeSolveVWith meth (X epsAbs epsRel) hi g y0 ts
    where
      g t x0 = coerce $ f t x0

-- | A version of 'odeSolveV' with reasonable default parameters and
-- system of equations defined using lists. FIXME: we should say
-- something about the fact we could use the Jacobian but don't for
-- compatibility with hmatrix-gsl.
odeSolve :: (Double -> [Double] -> [Double]) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
         -> [Double]                         -- ^ initial conditions
         -> Vector Double                    -- ^ desired solution times
         -> Matrix Double                    -- ^ solution
odeSolve f y0 ts =
  -- FIXME: These tolerances are different from the ones in GSL
  odeSolveVWith SDIRK_5_3_4' (XX' 1.0e-6 1.0e-10 1 1)  Nothing g (V.fromList y0) (V.fromList $ toList ts)
  where
    g t x0 = V.fromList $ f t (V.toList x0)

odeSolveVWith ::
  ODEMethod
  -> StepControl
  -> Maybe Double -- ^ initial step size - by default, ARKode
                  -- estimates the initial step size to be the
                  -- solution \(h\) of the equation
                  -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
                  -- \(\ddot{y}\) is an estimated value of the second
                  -- derivative of the solution at \(t_0\)
  -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
  -> V.Vector Double                     -- ^ Initial conditions
  -> V.Vector Double                     -- ^ Desired solution times
  -> Matrix Double                       -- ^ Error code or solution
odeSolveVWith method control initStepSize f y0 tt =
  case odeSolveVWith' opts method control initStepSize f y0 tt of
    Left  (c, _v) -> error $ show c -- FIXME
    Right (v, _d) -> v
  where
    opts = ODEOpts { maxNumSteps = 10000
                   , minStep     = 1.0e-12
                   , maxFail     = 10
                   , odeMethod   = error "ARKode: unexpected use of ODEOpts.odeMethod"
                   , stepControl = error "ARKode: unexpected use of ODEOpts.stepControl"
                   , initStep    = error "ARKode: unexpected use of ODEOpts.initStep"
                   }

odeSolveVWith' ::
  ODEOpts ODEMethod
  -> ODEMethod
  -> StepControl
  -> Maybe Double -- ^ initial step size - by default, ARKode
                  -- estimates the initial step size to be the
                  -- solution \(h\) of the equation
                  -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
                  -- \(\ddot{y}\) is an estimated value of the second
                  -- derivative of the solution at \(t_0\)
  -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
  -> V.Vector Double                     -- ^ Initial conditions
  -> V.Vector Double                     -- ^ Desired solution times
  -> Either (Matrix Double, Int) (Matrix Double, SundialsDiagnostics) -- ^ Error code or solution
odeSolveVWith' opts method control initStepSize f y0 tt =
  case solveOdeC (fromIntegral $ maxFail opts)
                 (fromIntegral $ maxNumSteps opts) (coerce $ minStep opts)
                 (fromIntegral $ getMethod method) (coerce initStepSize) jacH (scise control)
                 (coerce f) (coerce y0) (coerce tt) of
    Left  (v, c) -> Left  (reshape l (coerce v), fromIntegral c)
    Right (v, d)
      | V.null y0 -> Right ((V.length tt >< 0) [], emptyDiagnostics)
      | otherwise -> Right (reshape l (coerce v), d)
  where
    l = size y0
    scise (X aTol rTol)                          = coerce (V.replicate l aTol, rTol)
    scise (X' aTol rTol)                         = coerce (V.replicate l aTol, rTol)
    scise (XX' aTol rTol yScale _yDotScale)      = coerce (V.replicate l aTol, yScale * rTol)
    -- FIXME; Should we check that the length of ss is correct?
    scise (ScXX' aTol rTol yScale _yDotScale ss) = coerce (V.map (* aTol) ss, yScale * rTol)
    jacH = fmap (\g t v -> matrixToSunMatrix $ g (coerce t) (coerce v)) $
           getJacobian method
    matrixToSunMatrix m = T.SunMatrix { T.rows = nr, T.cols = nc, T.vals = vs }
      where
        nr = fromIntegral $ rows m
        nc = fromIntegral $ cols m
        -- FIXME: efficiency
        vs = V.fromList $ map coerce $ concat $ toLists m

-- | This function implements the same interface as
-- 'Numeric.Sundials.CVode.ODE.odeSolveWithEvents', although it does not
-- currently support events.
odeSolveWithEvents
  :: ODEOpts ODEMethod
  -> [EventSpec]
    -- ^ Event specifications
  -> Int
    -- ^ Maximum number of events
  -> (Double -> V.Vector Double -> V.Vector Double)
    -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
  -> Maybe (Double -> Vector Double -> Matrix Double)
    -- ^ The Jacobian (optional)
  -> V.Vector Double
    -- ^ Initial conditions
  -> V.Vector Double
    -- ^ Desired solution times
  -> Either Int SundialsSolution
    -- ^ Either an error code or a solution
odeSolveWithEvents opts events _ rhs _mb_jac y0 times
  | (not . null) events =
      -- Call error rather than return a Left because this is a programming
      -- error, not just a runtime issue.
      error $ "ARKode called with a non-empty list of events (" ++ show (length events) ++
      " in total).\
      \ ARKode does not support events at this point and should not be passed any."
  | otherwise =
      let
        result :: Either (Matrix Double, Int)
                         (Matrix Double, SundialsDiagnostics)
        result =
          odeSolveVWith' opts
            (odeMethod opts)
            (stepControl opts)
            (initStep opts)
            rhs y0 times
      in
        case result of
          Left (_, code) -> Left code
          Right (mx, diagn) ->
            Right $ SundialsSolution
                { actualTimeGrid = times
                , solutionMatrix =
                    -- Note: at this time, ARKode's output matrix does not
                    -- include the time column, so we're not dropping it
                    -- here unlike in CVode. If/when we add event support
                    -- to ARKode, this is going to change.
                    mx
                , eventInfo = []
                , diagnostics = diagn
                }

solveOdeC ::
  CInt ->
  SunIndexType ->
  CDouble ->
  CInt ->
  Maybe CDouble ->
  (Maybe (CDouble -> V.Vector CDouble -> T.SunMatrix)) ->
  (V.Vector CDouble, CDouble) ->
  (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
  -> V.Vector CDouble -- ^ Initial conditions
  -> V.Vector CDouble -- ^ Desired solution times
  -> Either (V.Vector CDouble, CInt) (V.Vector CDouble, SundialsDiagnostics) -- ^ Partial solution and error code or
                                                                             -- solution and diagnostics
solveOdeC maxErrTestFails maxNumSteps_ minStep_ method initStepSize
          jacH (aTols, rTol) fun f0 ts
  | V.null f0 = -- 0-dimensional (empty) system
    Right (V.empty, emptyDiagnostics)
  | otherwise =
  unsafePerformIO $ do
  let isInitStepSize :: CInt
      isInitStepSize = fromIntegral $ fromEnum $ isJust initStepSize
      ss :: CDouble
      ss = case initStepSize of
             -- It would be better to put an error message here but
             -- inline-c seems to evaluate this even if it is never
             -- used :(
             Nothing -> 0.0
             Just x  -> x

  let dim = V.length f0
      nEq :: SunIndexType
      nEq = fromIntegral dim
      nTs :: CInt
      nTs = fromIntegral $ V.length ts
  quasiMatrixRes <- createVector ((fromIntegral dim) * (fromIntegral nTs))
  qMatMut <- V.thaw quasiMatrixRes
  diagMut :: V.MVector _ SunIndexType <- V.thaw =<< createVector 10 -- FIXME
  -- We need the types that sundials expects. These are tied together
  -- in 'CLangToHaskellTypes'. FIXME: The Haskell type is currently empty!
  let funIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt
      funIO t y f _ptr = do
        sv <- peek y
        poke f $ T.SunVector { T.sunVecN = T.sunVecN sv
                             , T.sunVecVals = fun t (T.sunVecVals sv)
                             }
        [CU.exp| int{ 0 } |]
  let isJac :: CInt
      isJac = fromIntegral $ fromEnum $ isJust jacH
      jacIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunMatrix ->
               Ptr () -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunVector ->
               IO CInt
      jacIO t y _fy jacS _ptr _tmp1 _tmp2 _tmp3 = do
        case jacH of
          Nothing   -> error "Numeric.Sundials.ARKode.ODE: Jacobian not defined"
          Just jacI -> do j <- jacI t <$> (T.sunVecVals <$> peek y)
                          poke jacS j
                          -- FIXME: I don't understand what this comment means
                          -- Unsafe since the function will be called many times.
                          [CU.exp| int{ 0 } |]

  res <- [C.block| int {
                         /* general problem variables */

                         int flag;                  /* reusable error-checking flag                 */
                         int i, j;                  /* reusable loop indices                        */
                         N_Vector y = NULL;         /* empty vector for storing solution            */
                         N_Vector tv = NULL;        /* empty vector for storing absolute tolerances */
                         SUNMatrix A = NULL;        /* empty matrix for linear solver               */
                         SUNLinearSolver LS = NULL; /* empty linear solver object                   */
                         void *arkode_mem = NULL;   /* empty ARKode memory structure                */
                         realtype t;
                         long nst, nst_a, nfe, nfi, nsetups, nje, nfeLS, nni, ncfn, netf;

                         /* general problem parameters */

                         realtype T0 = RCONST(($vec-ptr:(double *ts))[0]); /* initial time              */
                         sunindextype NEQ = $(sunindextype nEq);             /* number of dependent vars. */

                         /* Initialize data structures */

                         y = N_VNew_Serial(NEQ); /* Create serial vector for solution */
                         if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1;
                         /* Specify initial condition */
                         for (i = 0; i < NEQ; i++) {
                           NV_Ith_S(y,i) = ($vec-ptr:(double *f0))[i];
                         };

                         tv = N_VNew_Serial(NEQ); /* Create serial vector for absolute tolerances */
                         if (check_flag((void *)tv, "N_VNew_Serial", 0)) return 1;
                         /* Specify tolerances */
                         for (i = 0; i < NEQ; i++) {
                           NV_Ith_S(tv,i) = ($vec-ptr:(double *aTols))[i];
                         };

                         /* Call ARKStepCreate to initialize the ARK timestepper module and */
                         /* specify the right-hand side function in y'=f(t,y), the inital time */
                         /* T0, and the initial dependent variable vector y.  Note: since this */
                         /* problem is fully implicit, we set f_E to NULL and f_I to f. */

                         /* Here we use the C types defined in helpers.h which tie up with */
                         /* the Haskell types defined in CLangToHaskellTypes               */
                         if ($(int method) < MIN_DIRK_NUM) {
                           arkode_mem = ARKStepCreate($fun:(int (* funIO) (double t, SunVector y[], SunVector dydt[], void * params)), NULL, T0, y);
                           if (check_flag((void *)arkode_mem, "ARKStepCreate", 0)) return 1;
                             } else {
                           arkode_mem = ARKStepCreate(NULL, $fun:(int (* funIO) (double t, SunVector y[], SunVector dydt[], void * params)), T0, y);
                           if (check_flag(&flag, "ARKStepCreate", 0)) return 1;
                         }

                         flag = ARKStepSetMinStep(arkode_mem, $(double minStep_));
                         if (check_flag(&flag, "ARKStepSetMinStep", 1)) return 1;
                         flag = ARKStepSetMaxNumSteps(arkode_mem, $(sunindextype maxNumSteps_));
                         if (check_flag(&flag, "ARKStepSetMaxNumSteps", 1)) return 1;
                         flag = ARKStepSetMaxErrTestFails(arkode_mem, $(int maxErrTestFails));
                         if (check_flag(&flag, "ARKStepSetMaxErrTestFails", 1)) return 1;

                         /* Set routines */
                         flag = ARKStepSVtolerances(arkode_mem, $(double rTol), tv);
                         if (check_flag(&flag, "ARKStepSVtolerances", 1)) return 1;

                         /* Initialize dense matrix data structure and solver */
                         A = SUNDenseMatrix(NEQ, NEQ);
                         if (check_flag((void *)A, "SUNDenseMatrix", 0)) return 1;
                         LS = SUNDenseLinearSolver(y, A);
                         if (check_flag((void *)LS, "SUNDenseLinearSolver", 0)) return 1;

                         /* Attach matrix and linear solver */
                         flag = ARKStepSetLinearSolver(arkode_mem, LS, A);
                         if (check_flag(&flag, "ARKStepSetLinearSolver", 1)) return 1;

                         /* Set the initial step size if there is one */
                         if ($(int isInitStepSize)) {
                           /* FIXME: We could check if the initial step size is 0 */
                           /* or even NaN and then throw an error                 */
                           flag = ARKStepSetInitStep(arkode_mem, $(double ss));
                           if (check_flag(&flag, "ARKStepSetInitStep", 1)) return 1;
                         }

                         /* Set the Jacobian if there is one */
                         if ($(int isJac)) {
                           flag = ARKStepSetJacFn(arkode_mem, $fun:(int (* jacIO) (double t, SunVector y[], SunVector fy[], SunMatrix Jac[], void * params, SunVector tmp1[], SunVector tmp2[], SunVector tmp3[])));
                           if (check_flag(&flag, "ARKStepSetJacFn", 1)) return 1;
                         }

                         /* Store initial conditions */
                         for (j = 0; j < NEQ; j++) {
                           ($vec-ptr:(double *qMatMut))[0 * $(int nTs) + j] = NV_Ith_S(y,j);
                         }

                         /* Explicitly set the method */
                         if ($(int method) >= MIN_DIRK_NUM) {
                           flag = ARKStepSetTableNum(arkode_mem, $(int method), -1);
                           if (check_flag(&flag, "ARKStepSetTableNum", 1)) return 1;
                         } else {
                           flag = ARKStepSetTableNum(arkode_mem, -1, $(int method));
                           if (check_flag(&flag, "ERKStepSetTableNum", 1)) return 1;
                         }

                         /* Main time-stepping loop: calls ARKStep to perform the integration */
                         /* Stops when the final time has been reached                       */
                         for (i = 1; i < $(int nTs); i++) {

                           flag = ARKStepEvolve(arkode_mem, ($vec-ptr:(double *ts))[i], y, &t, ARK_NORMAL); /* call integrator */
                           if (check_flag(&flag, "ARKStep solver failure, stopping integration", 1)) return 1;

                           /* Store the results for Haskell */
                           for (j = 0; j < NEQ; j++) {
                             ($vec-ptr:(double *qMatMut))[i * NEQ + j] = NV_Ith_S(y,j);
                           }
                         }

                         /* Get some final statistics on how the solve progressed */

                         flag = ARKStepGetNumSteps(arkode_mem, &nst);
                         check_flag(&flag, "ARKStepGetNumSteps", 1);
                         ($vec-ptr:(sunindextype *diagMut))[0] = nst;

                         flag = ARKStepGetNumStepAttempts(arkode_mem, &nst_a);
                         check_flag(&flag, "ARKStepGetNumStepAttempts", 1);
                         ($vec-ptr:(sunindextype *diagMut))[1] = nst_a;

                         flag = ARKStepGetNumRhsEvals(arkode_mem, &nfe, &nfi);
                         check_flag(&flag, "ARKStepGetNumRhsEvals", 1);
                         ($vec-ptr:(sunindextype *diagMut))[2] = nfe;
                         ($vec-ptr:(sunindextype *diagMut))[3] = nfi;

                         flag = ARKStepGetNumLinSolvSetups(arkode_mem, &nsetups);
                         check_flag(&flag, "ARKStepGetNumLinSolvSetups", 1);
                         ($vec-ptr:(sunindextype *diagMut))[4] = nsetups;

                         flag = ARKStepGetNumErrTestFails(arkode_mem, &netf);
                         check_flag(&flag, "ARKStepGetNumErrTestFails", 1);
                         ($vec-ptr:(sunindextype *diagMut))[5] = netf;

                         flag = ARKStepGetNumNonlinSolvIters(arkode_mem, &nni);
                         check_flag(&flag, "ARKStepGetNumNonlinSolvIters", 1);
                         ($vec-ptr:(sunindextype *diagMut))[6] = nni;

                         flag = ARKStepGetNumNonlinSolvConvFails(arkode_mem, &ncfn);
                         check_flag(&flag, "ARKStepGetNumNonlinSolvConvFails", 1);
                         ($vec-ptr:(sunindextype *diagMut))[7] = ncfn;

                         flag = ARKStepGetNumJacEvals(arkode_mem, &nje);
                         check_flag(&flag, "ARKStepGetNumJacEvals", 1);
                         ($vec-ptr:(sunindextype *diagMut))[8] = nje;

                         flag = ARKStepGetNumLinRhsEvals(arkode_mem, &nfeLS);
                         check_flag(&flag, "ARKStepGetNumLinRhsEvals", 1);
                         ($vec-ptr:(sunindextype *diagMut))[9] = nfeLS;

                         /* Clean up and return */
                         N_VDestroy(y);            /* Free y vector          */
                         N_VDestroy(tv);           /* Free tv vector         */
                         ARKStepFree(&arkode_mem);  /* Free integrator memory */
                         SUNLinSolFree(LS);        /* Free linear solver     */
                         SUNMatDestroy(A);         /* Free A matrix          */

                         return flag;
                       } |]
  preD <- V.freeze diagMut
  let d = SundialsDiagnostics (fromIntegral $ preD V.!0)
                              (fromIntegral $ preD V.!1)
                              (fromIntegral $ preD V.!2)
                              (fromIntegral $ preD V.!3)
                              (fromIntegral $ preD V.!4)
                              (fromIntegral $ preD V.!5)
                              (fromIntegral $ preD V.!6)
                              (fromIntegral $ preD V.!7)
                              (fromIntegral $ preD V.!8)
                              (fromIntegral $ preD V.!9)
  m <- V.freeze qMatMut
  if res == 0
    then do
      return $ Right (m, d)
    else do
      return $ Left  (m, res)

max_stars_repo_name: idontgetoutmuch/hmatrix-sundials

max_issues_repo_name: idontgetoutmuch/hmatrix-sundials

max_forks_repo_name: idontgetoutmuch/hmatrix-sundials


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ViewPatterns     #-}
{-# LANGUAGE TypeFamilies     #-}
{-# LANGUAGE ScopedTypeVariables #-}

-- |
-- Module      : Main
-- Description : Test runner
-- Copyright   : (c) Tom Westerhout, 2017
-- License     : BSD3
-- Maintainer  : t.westerhout@student.ru.nl
-- Stability   : experimental

module Main where

import           Prelude hiding (map, zipWithM)

import           Debug.Trace
import qualified System.Random.MWC as MWC
-- import           Control.Lens hiding((<.>))
import           Control.Monad.Reader hiding (zipWithM)
import           Control.Monad.Primitive
import qualified Data.List            as L
import           Data.Complex
import           Data.Semigroup
import           Data.Vector.Storable (Vector, (!))
import qualified Data.Vector.Storable as V
import           System.Exit
import           System.IO hiding (hGetLine)
import           System.Environment(getArgs, getProgName)
import           Foreign.Storable
import           Data.Text (Text)
import           Data.Text.IO (hGetLine)
import qualified Data.Text.IO as T

import           Lens.Micro
import           Lens.Micro.Extras

import           PSO.Random
import           PSO.Swarm
-- import           PSO.Heisenberg
-- import           PSO.Energy
-- import           PSO.Neural
import           PSO.FromPython
import           NQS.Rbm

meanVariance :: (RealFloat a, Storable a) => V.Vector a -> (a, a)
meanVariance xs
  | n > 1     = (m, sumVar m xs / fromIntegral (n - 1))
  | otherwise = (m, 0)
    where n = V.length xs
          m = V.sum xs / fromIntegral (V.length xs)
          sumVar m' = V.sum . V.map ((^2) . subtract m')

type R = Float
type C = Complex Float
type X = Rbm C

fromPyFile :: FilePath -> IO X
fromPyFile name = withFile name ReadMode toRbm
  where toRight :: Either String (V.Vector (Complex Float)) -> V.Vector (Complex Float)
        toRight (Right x) = x
        toRight (Left x)  = error x
        toRbm h = do
          a <- trace ("a...") $ toRight <$> readVector <$> hGetLine h
          b <- trace ("b...") $ toRight <$> readVector <$> hGetLine h
          s <- hGetLine h
          T.putStrLn s
          let !w = trace ("w...") $ toRight $ readMatrix s
          return $ mkRbm (a, b, w)

randomRbm ::
  ( UniformDist m R
  )
  => Int
  -> Int
  -> (R, R)
  -> (R, R)
  -> (R, R)
  -> m X
randomRbm n m (lowV, highV) (lowH, highH) (lowW, highW) =
  do
    visible <- uniformVector n (lowV :+ lowV, highV :+ highV)
    hidden <- uniformVector m (lowH :+ lowH, highH :+ highH)
    weights <- uniformVector (n * m) (lowW :+ lowW, highW :+ highW)
    return $ trace ("mkRbm...") (mkRbm (visible, hidden, weights))

rbms :: (UniformDist m R) => m [X]
rbms = trace ("Creating rbms...") $ replicateM 20 (randomRbm 10 20 (-0.1, 0.1) (-0.1, 0.1) (-0.1, 0.1))

rbms' :: IO [X]
rbms' = Prelude.mapM fromPyFile $
  [ "test.txt"
  -- "/home/tom/src/tcm-swarm/10_spins/Weights_10_spins_2_density_-1.23981661958_energy"
  -- "/home/tom/src/tcm-swarm/10_spins/Weights_10_spins_2_density_-1.89541865304_energy"
  -- , "/home/tom/src/tcm-swarm/10_spins/Weights_10_spins_2_density_0.0887681979723_energy"
  -- , "/home/tom/src/tcm-swarm/10_spins/Weights_10_spins_2_density_1.40938434022_energy"
  -- , "/home/tom/src/tcm-swarm/10_spins/Weights_10_spins_2_density_-3.61205064542_energy"
  -- , "/home/tom/src/tcm-swarm/10_spins/Weights_10_spins_2_density_-4.1917579599_energy"
  -- , "/home/tom/src/tcm-swarm/10_spins/Weights_10_spins_2_density_-5.41817224673_energy"
  -- , "/home/tom/src/tcm-swarm/10_spins/Weights_10_spins_2_density_-5.74064288804_energy"
  -- , "/home/tom/src/tcm-swarm/10_spins/Weights_10_spins_2_density_-6.33997058402_energy"
  -- , "/home/tom/src/tcm-swarm/10_spins/Weights_10_spins_2_density_-6.48266135073_energy"
  -- , "/home/tom/src/tcm-swarm/10_spins/Weights_10_spins_2_density_-7.29704001394_energy"
  -- , "/home/tom/src/tcm-swarm/10_spins/Weights_10_spins_2_density_-7.4024108663_energy"
  -- , "/home/tom/src/tcm-swarm/10_spins/Weights_10_spins_2_density_-7.83407453451_energy"
  -- , "/home/tom/src/tcm-swarm/10_spins/Weights_10_spins_2_density_-8.17005224883_energy"
  -- , "/home/tom/src/tcm-swarm/10_spins/Weights_10_spins_2_density_-8.18318151693_energy"
  -- , "/home/tom/src/tcm-swarm/10_spins/Weights_10_spins_2_density_-8.75410430005_energy"
  -- , "/home/tom/src/tcm-swarm/10_spins/Weights_10_spins_2_density_-9.04188789775_energy"
  -- , "/home/tom/src/tcm-swarm/10_spins/Weights_10_spins_2_density_-9.34302369964_energy"
  -- , "/home/tom/src/tcm-swarm/10_spins/Weights_10_spins_2_density_-9.5917604763_energy"
  -- , "/home/tom/src/tcm-swarm/10_spins/Weights_10_spins_2_density_2.01304877555_energy"
  ]

{-
  [ "/home/tom/src/tcm-swarm/8_spins/input_singleflip_0.txt"
  , "/home/tom/src/tcm-swarm/8_spins/input_singleflip_1.txt"
  , "/home/tom/src/tcm-swarm/8_spins/input_singleflip_2.txt"
  , "/home/tom/src/tcm-swarm/8_spins/input_singleflip_3.txt"
  , "/home/tom/src/tcm-swarm/8_spins/input_singleflip_4.txt"
  , "/home/tom/src/tcm-swarm/8_spins/input_singleflip_5.txt"
  , "/home/tom/src/tcm-swarm/8_spins/input_singleflip_6.txt"
  , "/home/tom/src/tcm-swarm/8_spins/input_singleflip_7.txt"
  , "/home/tom/src/tcm-swarm/8_spins/input_singleflip_8.txt"
  , "/home/tom/src/tcm-swarm/8_spins/input_singleflip_9.txt"
  , "/home/tom/src/tcm-swarm/8_spins/input_singleflip_10.txt"
  , "/home/tom/src/tcm-swarm/8_spins/input_singleflip_11.txt"
  , "/home/tom/src/tcm-swarm/8_spins/input_singleflip_12.txt"
  , "/home/tom/src/tcm-swarm/8_spins/input_singleflip_13.txt"
  , "/home/tom/src/tcm-swarm/8_spins/input_singleflip_14.txt"
  , "/home/tom/src/tcm-swarm/8_spins/input_singleflip_15.txt"
  , "/home/tom/src/tcm-swarm/8_spins/input_singleflip_16.txt"
  , "/home/tom/src/tcm-swarm/8_spins/input_singleflip_17.txt"
  , "/home/tom/src/tcm-swarm/8_spins/input_singleflip_18.txt"
  , "/home/tom/src/tcm-swarm/8_spins/input_singleflip_19.txt"
  ]
-}

function :: QMState X -> ReaderT g IO (Vector C)
function x = lift $
  sampleMoments (x ^. pos) Heisenberg1DPeriodic 4 (1000, 7000 + 1000, 7) (Just 0) 4

instance {-# OVERLAPS #-} Ord (Vector C) where
  a <= b = let k x = realPart (x ! 3) / (realPart (x ! 1))^2
            in k a <= k b

process :: Swarm m (SwarmGuide X (Vector C)) (BeeGuide X (Vector C)) (QMState X) (Vector C)
        -> IO ()
process swarm = do
  putStrLn $ "At iteration #" <> show (swarm ^. guide . iteration)
  withFile "Energies.dat" AppendMode $ \h ->
    do
      let xs = view bees swarm
      hPutStrLn h (L.concat $ L.intersperse "\t" $ show . realPart . (! 0) . view val <$> xs)
  withFile "Variances.dat" AppendMode $ \h ->
    do
      let xs = view bees swarm
      hPutStrLn h (L.concat $ L.intersperse "\t" $ show . realPart . (! 1) . view val <$> xs)
  withFile "Kurtosis.dat" AppendMode $ \h ->
    do
      let xs = view bees swarm
          k x = realPart (x ! 3) / (realPart (x ! 1))^2
      hPutStrLn h (L.concat $ L.intersperse "\t" $ show . k . view val <$> xs)

instance Scalable Float (Rbm C) where
  scale λ x = map (* (λ :+ 0)) x

instance DeltaWell m R R => DeltaWell m R X where
  upDeltaWell κ p x = zipWithM (upDeltaWell κ) p x

update ::
  ( m ~ ReaderT g IO
  , Randomisable m Float
  )
  => PhaseUpdater m (SwarmGuide X r) (BeeGuide X r) (QMState X) r
update = PhaseUpdater $ deltaUpdater (1.9 * log 2 :: R)

runHeisenbergFromList ::
  ( m ~ ReaderT g IO
  -- , r ~ MeanVar Float
  -- , χ ~ Rbm (Complex Float)
  -- , DeltaWell m Float χ
  -- , Scalable Float χ
  , UniformDist m Int
  , Randomisable m Float
  , Randomisable m Bool
  )
  => g -> IO ()
runHeisenbergFromList gen = runReaderT go gen
  where go =
          do
            -- states <- fmap QMState <$> rbms
            states <- fmap QMState <$> lift rbms'
            optimiseNDFromList
                  states
                  update
                  function
                  (\s -> (s^.guide.iteration == 100))
                  (lift . process)
            return ()

{-
runHeisenberg ::
  ( m ~ ReaderT g IO
  -- , r ~ MeanVar Float
  -- , χ ~ Rbm (Complex Float)
  -- , DeltaWell m Float χ
  -- , Scalable Float χ
  , UniformDist m Int
  , Randomisable m Float
  , Randomisable m Bool
  )
  => Int -> Int -> g -> IO ()
runHeisenberg spins count gen = do
  let h  = heisenberg1DOpen spins
      initBounds = ((-2.0E-1) :+ (-2.0E-1), 2.0E-1 :+ 2.0E-1)
      n = spins
      newState = QMState <$> uniformRbm n (2 * spins) initBounds
      func x = uncurry MV <$> energyHH1DOpenMKLC (x ^. pos) 100 10000
      xs = optimiseND
        newState
        update
        func
        count
        (\s -> (s^.guide.iteration == 100))
  swarms <- runReaderT xs gen
  writeEnergies2TSV "Energies.dat"  (view (val . mean)) swarms
  writeEnergies2TSV "Variances.dat" (view (val . var)) swarms
  -- writeVariances2TSV "Variances.dat" swarms
  let swarm = last swarms
      (m, v) = meanVariance
                . V.fromList
                . map (view (val . mean))
                $ swarm ^. bees
      min = minimum $ view mean <$> view val <$> swarm ^. bees
      max = maximum $ view mean <$> view val <$> swarm ^. bees
  -- mapM_ (print . (!!0) . (view bees)) $ swarms
  putStrLn ""
  putStrLn $ "[+] After " ++ show (swarm ^. guide . iteration)
    ++ " iterations: "
  putStrLn $ "[+] E[<H>]         = " ++ show m
  putStrLn $ "[+] Min[<H>]       = " ++ show min
  putStrLn $ "[+] Max[<H>]       = " ++ show max
  putStrLn $ "[+] StdDev[<H>]    = " ++ show (sqrt v)
  putStrLn $ "[+] Best[<H>]      = " ++ show (swarm ^. guide . val . mean)
  putStrLn $ "[+] Best[Var[<H>]] = " ++ show (swarm ^. guide . val . var)
  putStrLn $ "[+] Actual states:   " ++ show (fst . LA.eigSH $ h)
-}

main :: IO ()
main = do
  args <- getArgs
  case args of
    [ reads -> [(spins :: Int, _)], reads -> [(count :: Int, _)], reads -> [(seed, _)]
      ] -> do g <- mkMWCGen (Just seed)
              runHeisenbergFromList g
    _ -> do name <- getProgName
            hPutStrLn stderr $ "usage: " ++ name ++ " <#spins> <#bees> <seed>"
            exitFailure

max_stars_repo_name: twesterhout/tcm-swarm

max_issues_repo_name: twesterhout/tcm-swarm

max_forks_repo_name: twesterhout/tcm-swarm


text:
{-# LANGUAGE BangPatterns #-}
{- | Parser for MNIST data.
Format: n lines of 785 double precision numbers, first 784 are the image, last is the expected number.
-}
module Main where

import           Data.Attoparsec.ByteString.Char8 (Parser, char, double, parseOnly, sepBy)
import qualified Data.ByteString.Char8 as BS
import           Numeric.LinearAlgebra (R, Vector, assoc, vector)
import qualified Data.Vector as V

data Sample = Sample
              { sampleInput          :: {-# UNPACK #-} !(Vector R)
              , sampleExpectedOutput :: {-# UNPACK #-} !(Vector R)
              } deriving (Read, Show)

makeSample :: ([R], [R]) -> Sample
makeSample (image, [expected]) = Sample (vector image) (assoc 10 0 [(floor expected, 1)])

sample :: Parser Sample
sample = do
  x <- double `sepBy` (char ' ')
  let !s = makeSample $ splitAt 784 x
  pure s

toSampleSet :: V.Vector Sample -> (V.Vector Sample, V.Vector Sample, V.Vector Sample)
toSampleSet vec = (train, validation, test)
    where (train, rest) = V.splitAt 50000 vec
          (validation, test) = V.splitAt 10000 rest

parseMnist :: BS.ByteString -> (V.Vector Sample, V.Vector Sample, V.Vector Sample)
parseMnist = toSampleSet . V.unfoldr parseLines . BS.lines
    where parseLines :: [BS.ByteString] -> Maybe (Sample, [BS.ByteString])
          parseLines []     = Nothing
          parseLines (l:ls) = Just (fromRight $ parseOnly sample l, ls)

          fromRight (Right b) = b

main :: IO ()
main = do
    !contents <- BS.readFile "/home/matt/Downloads/mnist.dat"
    let (!train, !validation, !test) = parseMnist contents
    print $
        show (length train)      ++ " " ++
        show (length validation) ++ " " ++
        show (length test)

max_stars_repo_name: m-renaud/ML

max_issues_repo_name: m-renaud/ML

max_forks_repo_name: m-renaud/ML


text:
{-# LANGUAGE NoImplicitPrelude #-}
module Test where
import Data.Complex.Compat

max_stars_repo_name: RyanGlScott/cabal-sdist-bug

max_issues_repo_name: RyanGlScott/cabal-sdist-bug

max_forks_repo_name: RyanGlScott/cabal-sdist-bug


text:
{-# LANGUAGE RecordWildCards #-}
module Main where

import Control.Monad (when)
import qualified Data.Text.IO as TIO
import qualified Data.ByteString.Lazy as BL (readFile)
import Data.Csv (decodeByName)

import QuoteData
import Statistics
import StatReport
import Charts
import Params

main :: IO ()
main = cmdLineParser >>= work

work :: Params -> IO ()
work params = do
  csvData <- BL.readFile (fname params)
  case decodeByName csvData of
    Left err -> putStrLn err
    Right (_, quotes) -> generateReports params quotes

generateReports :: (Functor t, Foldable t) => Params -> t QuoteData -> IO ()
generateReports Params {..} quotes = do
  TIO.putStr $ statReport statInfo'
  when prices $ plotChart title quotes [Open, Close, High, Low] fname_prices
  when volumes $ plotChart title quotes [Volume] fname_volumes
 where
   statInfo' = statInfo quotes
   withCompany pref = if company /= "" then pref ++ company else ""
   img_suffix = withCompany "_" ++ ".svg"
   fname_prices = "prices" ++ img_suffix
   fname_volumes = "volumes" ++ img_suffix
   title = "Historical Quotes" ++ withCompany " for "

max_stars_repo_name: Francososa/stockell

max_issues_repo_name: Francososa/stockell

max_forks_repo_name: Francososa/stockell


text:
{-# LANGUAGE InstanceSigs #-}
{-# LANGUAGE FlexibleContexts #-}

module Matrix.Traversable
  (
    matrixTraverse
  ) where

import Numeric.LinearAlgebra.HMatrix hiding (corr)

import Data.Traversable

import Util.Tuples

-- Matrix cannot be an instance of Functor, so it cannot be an instance of traversable
-- instead, convert to list and use the list traverse
matrixTraverse :: (Element b, Element c) => (a -> b -> (a,c)) -> a -> Matrix b -> (a, Matrix c)
matrixTraverse f a m = applySnd (fromLists) $ mapAccumL (mapAccumL f) a (toLists m)

max_stars_repo_name: eklinkhammer/neural-algorithms

max_issues_repo_name: eklinkhammer/neural-algorithms

max_forks_repo_name: eklinkhammer/neural-algorithms


text:
module Ch03.Activation where

import Prelude hiding (and, or)
import Numeric.LinearAlgebra hiding (step)
import Graphics.Gnuplot.Simple

-- Activation functions
step :: [R] -> [R]
step = map (\x -> if x > 0.0 then 1.0 else 0.0)

sigmoid :: [R] -> [R]
sigmoid = map (\x -> 1 / (1 + exp (-x)))

relu :: [R] -> [R]
relu = map $ max 0

-- Graphs
step_graph :: IO ()
step_graph = plotPath  [Title "step", YRange (-0.2, 1.2)] $ zip xs ys
  where
    xs = linearScale 1000 (-5.0, 5.0)
    ys = step xs

sigmoid_graph :: IO ()
sigmoid_graph = plotPath  [Title "sigmoid", YRange (-0.2, 1.2)] $ zip xs ys
  where
    xs = linearScale 1000 (-5.0, 5.0)
    ys = sigmoid xs

sigmoid_step_graph :: IO ()
sigmoid_step_graph = plotPathsStyle
                     [Title "sigmoid and step", YRange (-0.2, 1.2)]
                     [(style1, zip xs ys1), (style2, zip xs ys2)]
  where
    xs     = linearScale 1000 (-5.0, 5.0)
    ys1    = sigmoid xs
    ys2    = step xs
    style1 = PlotStyle Lines $ CustomStyle [(LineTitle "sigmoid")]
    style2 = PlotStyle Lines $ CustomStyle [(LineTitle "step")]

relu_graph :: IO ()
relu_graph = plotPath [Title "ReLU", YRange (-1.0, 5.5)] $ zip xs ys
  where
    xs = linearScale 1000 (-5.0, 5.0)
    ys = relu xs


text:
--stolen from kmeans!

module Math.Probably.KMedoids
    where
import Data.List (transpose, sort, groupBy, minimumBy, nub)
import Data.Ord (comparing)
--import Numeric.LinearAlgebra
import qualified Data.Map.Strict as Map

import Numeric.LinearAlgebra
import Control.Monad
import Debug.Trace
import Data.Random

type DistMap = ([Int] , (Int -> Int -> Double))

kmedoids :: Int -> Int -> DistMap -> RVar (Map.Map Int [Int])
kmedoids nclust iters dist = do
  inits <- kmedoidInit nclust dist
  return $ kmedoidIter dist inits iters

kmedoidsMany :: Int -> Int -> Int -> DistMap -> RVar [Int]
kmedoidsMany nclust iters tries dist = do
  ress <- forM [1..tries] $ \i -> do
     clusters <- kmedoids nclust iters dist
     let cost = totalCost dist clusters
     return $ trace ("KMed try = "++show i++ " nclus="++show nclust++" cost= "++show cost) (clusters, cost)
  return $ Map.keys $ fst $ minimumBy (comparing snd) ress

kmedoidInit :: Int -> DistMap -> RVar [Int]
kmedoidInit nclust dist = fmap (sort . take nclust) $ shuffle $  fst dist

kmedoidIter :: DistMap -> [Int] -> Int -> Map.Map Int [Int]
kmedoidIter dists meds 0 = assignToClosestMedoid dists meds
kmedoidIter dists meds0 iter = next where
  clusters = assignToClosestMedoid dists meds0
  meds1 = sort $ nub $ map (\(_, pts)-> reassignMedoid dists pts) $ Map.toList clusters
  next = {-trace (show (meds0, totalCost dists clusters))
         $ -} if meds1 == meds0
                 then assignToClosestMedoid dists meds1
                 else kmedoidIter dists meds1 (iter-1)

totalCost :: DistMap -> Map.Map Int [Int] -> Double
totalCost dists  = sum . map ccost . Map.toList where
  ccost (med, pts) = sum $ map (getDist dists med) pts

getDist :: DistMap -> Int -> Int -> Double
getDist (_,f) from to
  | from == to = 0
  | otherwise = f from to

assignToClosestMedoid :: DistMap -> [Int] -> Map.Map Int [Int]
assignToClosestMedoid dists medoids = clusmap where
  assignmap = Map.fromList $ map (\k -> (k, findNearestMedoid dists medoids k)) $ fst dists
  clusmap = Map.fromList $ map (\med -> (med, keysWhereValIs assignmap med)) medoids

keysWhereValIs :: Eq a => Map.Map k a -> a -> [k]
keysWhereValIs mp val = Map.keys $ Map.filter (==val) mp

findNearestMedoid :: DistMap -> [Int] -> Int -> Int
findNearestMedoid dists meds k
  = fst $ minimumBy (comparing snd) $ map (\med-> (med, getDist dists k med)) meds

calcCost :: DistMap -> Int -> [Int] -> Double
calcCost dists med  = sum . map (getDist dists med)

reassignMedoid :: DistMap -> [Int] -> Int
reassignMedoid dists points
  = fst $ minimumBy (comparing snd) $ map (\pt-> (pt, calcCost dists pt points) ) points

max_stars_repo_name: glutamate/probably-baysig

max_issues_repo_name: silky/probably-baysig

max_forks_repo_name: silky/probably-baysig


text:
{-# OPTIONS_GHC -fno-warn-missing-signatures #-}

module HVX.BlackBoxTests.ManyvarTestSpec where

import Test.Hspec
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Util (ones, zeros)

import HVX
import HVX.Internal.TestUtil

-- Problem definition.
nx = 2
ny = 3
nz = 4
a = EConst $ (3><nx)
    ([-0.1022 , 0.3129
    , -0.2414 , -0.8649
    , 0.3192 , -0.0301 ] :: [Double])
b = EConst $ (3><ny)
    ([-0.1649 , 1.1093 , -1.2141
    , 0.6277 , -0.8637 , -1.1135
    , 1.0933 , 0.0774 , -0.0068 ] :: [Double])
c = EConst $ (3><nz)
    ([1.5326 , -0.2256 , 0.0326 , 1.5442
    , -0.7697 , 1.1174 , 0.5525 , 0.0859
    , 0.3714 , -1.0891 , 1.1006 , -1.4916 ] :: [Double])
d = EConst $ (3><1)
    ([-0.6156 , 0.7481 , -0.1924] :: [Double])

x = EVar "x"
y = EVar "y"
z = EVar "z"
const10nx = EConst $ scale 10 $ ones nx 1
const10ny = EConst $ scale 10 $ ones ny 1
const01nz = EConst $ scale 0.1 $ ones nz 1

subgradAns = subgradMaximize
  ( neg (norm 2 (a *~ x +~ b *~ y +~ c *~ z +~ d))
    +~ neg (norm 1 x)
    +~ neg (norm 4.2 y) )
  [ const01nz >=~ hexp z
    , powBaseP01 0.25 x >=~ const10nx
    , powBaseP1InfNotInt 1.75 y <=~ const10ny ]
  (decNonSumStep 100.0) 10
  [("x", zeros nx 1), ("y", zeros ny 1), ("z", zeros nz 1)]

ellipsoidAns = ellipsoidMaximize
  ( neg (norm 2 (a *~ x +~ b *~ y +~ c *~ z +~ d))
    +~ neg (norm 1 x)
    +~ neg (norm 4.2 y) )
  [ const01nz >=~ hexp z
    , powBaseP01 0.25 x >=~ const10nx
    , powBaseP1InfNotInt 1.75 y <=~ const10ny ]
  [("x", nx), ("y", ny), ("z", nz)]
  1e-16 1e10

(subgradVars, subgradOptval) = subgradAns

(ellipsoidVars, ellipsoidOptval, ellipsoidUBound) = ellipsoidAns

-- CVX's results.
cvxOptval = -29051

-- Verify that HVX matches CVX.
-- TODO(mh): This test fails because primitives that generate implicit
-- constraints are currently unsupported. When support is added for them, the
-- test will be added back in. (2014-06-04)
spec :: Spec
spec =
  describe "Placeholder test" $ do
    it "placeholder description" $
      True `shouldBe` True
--  describe "Verify that HVX matches CVX for triple variable stress test" $ do
--    it "Verify that HVX subgrad matches CVX for huber/berhu" $
--      subgradOptval `shouldSatisfy` fpequalsApprox cvxOptval
--    it "Verify that HVX ellipsoid matches CVX for huber/berhu" $
--      ellipsoidOptval `shouldSatisfy` fpequalsApprox cvxOptval
--
--main :: IO ()
--main = hspec spec

max_stars_repo_name: wellposed/hvx

max_issues_repo_name: wellposed/hvx

max_forks_repo_name: wellposed/hvx