text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE IncoherentInstances #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}
module Data.C
  ( CVal (..)
  , CTy(..)
  , IsFun
  , eraseTy
  , getCTyR
  , getTy
  , eitherTy
  , pairTy
  , cTySpec
  , (<::)
  , tagFld
  , valFld
  , inlFld
  , inrFld
  , fstFld
  , sndFld
  , cCase
  , cTagl
  , cTagr
  , declVar
  , Int
  , Bool
  , Float
  , Double
  , String
  , Either
  ) where

import Control.Monad.CGen
import Type.Reflection ( Typeable, TypeRep, typeRep )
import Data.Complex

data CTy a where
  CUnit   :: CTy ()
  CInt    :: CTy Int
  CBool   :: CTy Bool
  CFlt    :: CTy Float
  CDbl    :: CTy Double
  CCplx   :: CTy (Complex Double)
  CStr    :: CTy String
  CPair   :: (CVal a, CVal b) => CTy a -> CTy b -> CTy (a,b)
  CEither :: (CVal a, CVal b) => CTy a -> CTy b -> CTy (Either a b)
  CVec    :: CVal a => CTy a -> CTy [a]

eraseTy :: CTy a -> ECTy
eraseTy CUnit         = ECUnit
eraseTy CInt          = ECInt
eraseTy CBool         = ECBool
eraseTy CFlt          = ECFlt
eraseTy CDbl          = ECDbl
eraseTy CCplx         = ECCplx
eraseTy CStr          = ECStr
eraseTy (CPair   l r) = ECPair (eraseTy l) (eraseTy r)
eraseTy (CEither l r) = ECEither (eraseTy l) (eraseTy r)
eraseTy (CVec    l  ) = ECVec (eraseTy l)

getTy :: forall a t. CVal a => t a -> ECTy
getTy _ = eraseTy (getCTy :: CTy a)

getCTyR :: CTy a -> TypeRep a
getCTyR (CPair _ _) = typeRep
getCTyR (CEither _ _) = typeRep
getCTyR (CVec _) = typeRep
getCTyR CUnit = typeRep
getCTyR CInt = typeRep
getCTyR CBool = typeRep
getCTyR CFlt = typeRep
getCTyR CDbl = typeRep
getCTyR CCplx = typeRep
getCTyR CStr = typeRep

eitherTy :: (CVal a, CVal b) => t a -> t b -> TypeRep (Either a b)
eitherTy _ _ = typeRep

pairTy :: (CVal a, CVal b) => t a -> t b -> TypeRep (a, b)
pairTy _ _ = typeRep

cTySpec :: forall a t st. CVal a => t a -> CGen st ([CDeclSpec], [CDerivedDeclr])
cTySpec _ = typeSpec $ eraseTy $ (getCTy :: CTy a)

type family IsFun f where
  IsFun (a -> b) = 'True
  IsFun a = 'False

class (Eq a, Ord a, Show a, Typeable a, IsFun a ~ 'False) => CVal a where
  getCTy  :: CTy a
  cVal    :: a -> CGen st CExpr

instance CVal () where
  getCTy = CUnit
  cVal () = cTySpec CUnit *> pure (cVar cUnit)

instance CVal Int where
  getCTy = CInt
  cVal i = pure $ cInt i

instance CVal Bool where
  getCTy = CBool
  cVal True = pure $ cInt 1
  cVal False = pure $ cInt 0

instance CVal Float where
  getCTy = CFlt
  cVal i = pure $ cFlt i

instance CVal Double where
  getCTy = CDbl
  cVal i = pure $ cDbl i

instance Ord (Complex Double) where
  compare i j =
    case compare (realPart i) (realPart j) of
      EQ -> compare (imagPart i) (imagPart j)
      o -> o

instance CVal (Complex Double) where
  getCTy = CCplx
  cVal i = pure $
    CBinary CAddOp (cDbl (realPart i))
    (CBinary CMulOp (cDbl (imagPart i)) (cVar $ internalIdent "I") undefNode)
    undefNode

instance (CVal a, CVal b) => CVal (a, b) where
  getCTy = CPair getCTy getCTy
  cVal (x, y) = do
    vx <- cVal x
    vy <- cVal y
    (t, q) <- cTySpec $ getCTy @(a,b)
    cVar <$> newVar t q (Just $ ini [mkInit vx, mkInit vy])
      where
        mkInit e = ([], CInitExpr e undefNode)

instance (CVal a, CVal b) => CVal (Either a b) where
  getCTy = CEither getCTy getCTy
  cVal (Left x) = do
    vx <- cVal x
    (t, q) <- cTySpec $ getCTy @(Either a b)
    cVar <$> newVar t q (Just $ ini [tagL, ([], ini [([mkInL], iniE vx)])])
      where
        tagL = ([], iniE $ cVar cTagl)
        mkInL = CMemberDesig inlFld undefNode
  cVal (Right x) = do
    vx <- cVal x
    (t, q) <- cTySpec $ getCTy @(Either a b)
    cVar <$> newVar t q (Just $ ini [tagR, ([], ini [([mkInR], iniE vx)])])
      where
        tagR = ([], iniE $ cVar cTagr)
        mkInR = CMemberDesig inrFld undefNode

ini :: CInitializerList NodeInfo -> CInitializer NodeInfo
ini = (`CInitList` undefNode)
iniE :: CExpression NodeInfo -> CInitializer NodeInfo
iniE = (`CInitExpr` undefNode)

instance (CVal a) => CVal [a] where
  getCTy = CVec getCTy
  cVal _xs = error "FIXME: unimplemented"

instance {-# OVERLAPPING #-} CVal String where
  getCTy = CStr
  cVal s = pure $ cStr s

(<::) :: forall a t st. CVal a => Ident -> t a -> CGen st [CBlockItem]
v <:: _ = do
  (t, q) <- cTySpec (getCTy :: CTy a)
  pure $ [CBlockDecl $ varDecl v t q Nothing]

cCase :: CExpr -> [CBlockItem] -> [CBlockItem] -> CStat
cCase e l r = CSwitch e cases undefNode
  where
    cases = cComp $ zipWith3 mkCase [True, False] [l, r] [cTagl, cTagr]
    mkCase bb ss tag = CBlockStmt $ CCase (cVar tag) (caseCode bb ss) undefNode
    caseCode b ss = cComp $ ss ++ mbreak
      where
        mbreak = if b then [CBlockStmt $ CBreak undefNode] else []

declVar :: forall a t st. CVal a => t a -> CGen st (CExpr, [CBlockItem])
declVar e
  | getTy e == ECUnit = pure (cVar cUnit, [])
  | otherwise = do
      v <- freshVar
      dv <- v <:: e
      pure (cVar v, dv)

max_stars_repo_name: session-arr/session-arr

max_issues_repo_name: session-arr/session-arr

max_forks_repo_name: dcastrop/SAlg


text:
{-|
Module: MachineLearning.NeuralNetwork.Regularization
Description: Regularization
Copyright: (c) Alexander Ignatyev, 2017
License: BSD-3
Stability: experimental
Portability: POSIX

Regularization.
-}

module MachineLearning.NeuralNetwork.Regularization
(
  Regularization(..)
  , forwardReg
  , backwardReg
)

where

import MachineLearning.Types (R, Matrix)
import qualified Numeric.LinearAlgebra as LA
import MachineLearning.Regularization (Regularization(..))

-- | Calcaulates regularization for forward propagation.
-- It takes regularization parameter and theta list.
forwardReg :: Regularization -> [(Matrix, Matrix)] -> R
forwardReg RegNone _ = 0
forwardReg (L2 lambda) thetaList = 0.5 * lambda * (sum $ map LA.norm_2 $ snd $ unzip thetaList)

-- | Calculates regularization for step of backward propagation.
-- It takes regularization parameter and theta.
backwardReg :: Regularization -> Matrix -> Matrix
backwardReg RegNone _ = 0
backwardReg (L2 lambda) w = w * (LA.scalar lambda)

ext: hs

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
-- | Type classes for random generation of values.
{-# LANGUAGE CPP #-}
#ifndef NO_GENERICS
{-# LANGUAGE DefaultSignatures, FlexibleContexts, TypeOperators #-}
#endif
module Test.QuickCheck.Arbitrary
  (
  -- * Arbitrary and CoArbitrary classes
    Arbitrary(..)
  , CoArbitrary(..)

  -- ** Helper functions for implementing arbitrary
  , arbitrarySizedIntegral        -- :: Integral a => Gen a
  , arbitraryBoundedIntegral      -- :: (Bounded a, Integral a) => Gen a
  , arbitrarySizedBoundedIntegral -- :: (Bounded a, Integral a) => Gen a
  , arbitrarySizedFractional      -- :: Fractional a => Gen a
  , arbitraryBoundedRandom        -- :: (Bounded a, Random a) => Gen a
  , arbitraryBoundedEnum          -- :: (Bounded a, Enum a) => Gen a
  -- ** Helper functions for implementing shrink
#ifndef NO_GENERICS
  , genericShrink      -- :: (Generic a, Typeable a, RecursivelyShrink (Rep a), Subterms (Rep a)) => a -> [a]
  , subterms           -- :: (Generic a, Subterms (Rep a)) => a -> [a]
  , recursivelyShrink  -- :: (Generic a, RecursivelyShrink (Rep a)) => a -> [a]
#endif
  , shrinkNothing            -- :: a -> [a]
  , shrinkList               -- :: (a -> [a]) -> [a] -> [[a]]
  , shrinkIntegral           -- :: Integral a => a -> [a]
  , shrinkRealFrac           -- :: RealFrac a => a -> [a]
  , shrinkRealFracToInteger  -- :: RealFrac a => a -> [a]
  -- ** Helper functions for implementing coarbitrary
  , coarbitraryIntegral      -- :: Integral a => a -> Gen b -> Gen b
  , coarbitraryReal          -- :: Real a => a -> Gen b -> Gen b
  , coarbitraryShow          -- :: Show a => a -> Gen b -> Gen b
  , coarbitraryEnum          -- :: Enum a => a -> Gen b -> Gen b
  , (><)

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

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

import System.Random(Random)
import Test.QuickCheck.Gen
import Test.QuickCheck.Gen.Unsafe

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

import Data.Char
  ( chr
  , ord
  , isLower
  , isUpper
  , toLower
  , isDigit
  , isSpace
  )

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

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

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

import Data.List
  ( sort
  , nub
  )

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

import Data.Int(Int8, Int16, Int32, Int64)
import Data.Word(Word, Word8, Word16, Word32, Word64)

#ifndef NO_GENERICS
import GHC.Generics
import Data.Typeable
#endif

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

-- | Random generation and shrinking of values.
class Arbitrary a where
  -- | A generator for values of the given type.
  arbitrary :: Gen a
  arbitrary = error "no default generator"

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

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

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

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

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

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

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

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

instance RecursivelyShrink U1 where
  grecursivelyShrink U1 = []

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

class Subterms f where
  gsubterms :: Typeable b => f a -> [b]

instance (Subterms f, Subterms g) => Subterms (f :*: g) where
  gsubterms (x :*: y) =
    gsubterms x ++ gsubterms y

instance (Subterms f, Subterms g) => Subterms (f :+: g) where
  gsubterms (L1 x) = gsubterms x
  gsubterms (R1 x) = gsubterms x

instance Subterms f => Subterms (M1 i c f) where
  gsubterms (M1 x) = gsubterms x

instance Typeable a => Subterms (K1 i a) where
  gsubterms (K1 x) =
    case cast x of
      Nothing -> []
      Just y -> [y]

instance Subterms U1 where
  gsubterms U1 = []

#endif

-- instances

instance (CoArbitrary a, Arbitrary b) => Arbitrary (a -> b) where
  arbitrary = promote (`coarbitrary` arbitrary)

instance Arbitrary () where
  arbitrary = return ()

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

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

instance Arbitrary a => Arbitrary (Maybe a) where
  arbitrary = frequency [(1, return Nothing), (3, liftM Just arbitrary)]

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

instance (Arbitrary a, Arbitrary b) => Arbitrary (Either a b) where
  arbitrary = oneof [liftM Left arbitrary, liftM Right arbitrary]

  shrink (Left x)  = [ Left  x' | x' <- shrink x ]
  shrink (Right y) = [ Right y' | y' <- shrink y ]

instance Arbitrary a => Arbitrary [a] where
  arbitrary = sized $ \n ->
    do k <- choose (0,n)
       sequence [ arbitrary | _ <- [1..k] ]

  shrink xs = shrinkList shrink xs

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

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

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

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

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

instance (RealFloat a, Arbitrary a) => Arbitrary (Complex a) where
  arbitrary = liftM2 (:+) arbitrary arbitrary
  shrink (x :+ y) = [ x' :+ y | x' <- shrink x ] ++
                    [ x :+ y' | y' <- shrink y ]

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

instance (Arbitrary a, Arbitrary b)
      => Arbitrary (a,b)
 where
  arbitrary = liftM2 (,) arbitrary arbitrary

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

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

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

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

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

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

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

-- typical instance for primitive (numerical) types

instance Arbitrary Integer where
  arbitrary = arbitrarySizedIntegral
  shrink    = shrinkIntegral

instance Arbitrary Int where
  arbitrary = arbitrarySizedIntegral
  shrink    = shrinkIntegral

instance Arbitrary Int8 where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink    = shrinkIntegral

instance Arbitrary Int16 where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink    = shrinkIntegral

instance Arbitrary Int32 where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink    = shrinkIntegral

instance Arbitrary Int64 where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink    = shrinkIntegral

instance Arbitrary Word where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink    = shrinkIntegral

instance Arbitrary Word8 where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink    = shrinkIntegral

instance Arbitrary Word16 where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink    = shrinkIntegral

instance Arbitrary Word32 where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink    = shrinkIntegral

instance Arbitrary Word64 where
  arbitrary = arbitrarySizedBoundedIntegral
  shrink    = shrinkIntegral

instance Arbitrary Char where
  arbitrary = chr `fmap` oneof [choose (0,127), choose (0,255)]
  shrink c  = filter (<. c) $ nub
            $ ['a','b','c']
           ++ [ toLower c | isUpper c ]
           ++ ['A','B','C']
           ++ ['1','2','3']
           ++ [' ','\n']
   where
    a <. b  = stamp a < stamp b
    stamp a = ( (not (isLower a)
              , not (isUpper a)
              , not (isDigit a))
              , (not (a==' ')
              , not (isSpace a)
              , a)
              )

instance Arbitrary Float where
  arbitrary = arbitrarySizedFractional
  shrink    = shrinkRealFrac

instance Arbitrary Double where
  arbitrary = arbitrarySizedFractional
  shrink    = shrinkRealFrac

-- ** Helper functions for implementing arbitrary

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

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

-- | Generates a fractional number. The number can be positive or negative
-- and its maximum absolute value depends on the size parameter.
arbitrarySizedFractional :: Fractional a => Gen a
arbitrarySizedFractional =
  sized $ \n ->
    let n' = toInteger n in
      do a <- choose ((-n') * precision, n' * precision)
         b <- choose (1, precision)
         return (fromRational (a % b))
 where
  precision = 9999999999999 :: Integer

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

-- | Generates an integral number. The number is chosen uniformly from
-- the entire range of the type. You may want to use
-- 'arbitrarySizedBoundedIntegral' instead.
arbitraryBoundedIntegral :: (Bounded a, Integral a) => Gen a
arbitraryBoundedIntegral =
  withBounds $ \mn mx ->
  do n <- choose (toInteger mn, toInteger mx)
     return (fromInteger n)

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

-- | Generates an element of a bounded enumeration.
arbitraryBoundedEnum :: (Bounded a, Enum a) => Gen a
arbitraryBoundedEnum =
  withBounds $ \mn mx ->
  do n <- choose (fromEnum mn, fromEnum mx)
     return (toEnum n)

-- | Generates an integral number from a bounded domain. The number is
-- chosen from the entire range of the type, but small numbers are
-- generated more often than big numbers. Inspired by demands from
-- Phil Wadler.
arbitrarySizedBoundedIntegral :: (Bounded a, Integral a) => Gen a
arbitrarySizedBoundedIntegral =
  withBounds $ \mn mx ->
  sized $ \s ->
    do let bits n | n `quot` 2 == 0 = 0
                  | otherwise = 1 + bits (n `quot` 2)
           k  = 2^(s*(bits mn `max` bits mx `max` 40) `div` 100)
       n <- choose (toInteger mn `max` (-k), toInteger mx `min` k)
       return (fromInteger n)

-- ** Helper functions for implementing shrink

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

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

-- | Shrink a fraction, but only shrink to integral values.
shrinkRealFracToInteger :: RealFrac a => a -> [a]
shrinkRealFracToInteger x =
  nub $
  [ -x
  | x < 0
  ] ++
  map fromInteger (shrinkIntegral (truncate x))

-- | Shrink a fraction.
shrinkRealFrac :: RealFrac a => a -> [a]
shrinkRealFrac x =
  nub $
  shrinkRealFracToInteger x ++
  [ x - x'
  | x' <- take 20 (iterate (/ 2) x)
  , (x - x') << x ]
 where
  a << b = abs a < abs b

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

-- | Used for random generation of functions.
class CoArbitrary a where
  -- | Used to generate a function of type @a -> b@.
  -- The first argument is a value, the second a generator.
  -- You should use 'variant' to perturb the random generator;
  -- the goal is that different values for the first argument will
  -- lead to different calls to 'variant'. An example will help:
  --
  -- @
  -- instance CoArbitrary a => CoArbitrary [a] where
  --   coarbitrary []     = 'variant' 0
  --   coarbitrary (x:xs) = 'variant' 1 . coarbitrary (x,xs)
  -- @

  coarbitrary :: a -> Gen b -> Gen b

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

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

instance CoArbitrary () where
  coarbitrary _ = id

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

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

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

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

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

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

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

instance (RealFloat a, CoArbitrary a) => CoArbitrary (Complex a) where
  coarbitrary (x :+ y) = coarbitrary x . coarbitrary y

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

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

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

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

-- typical instance for primitive (numerical) types

instance CoArbitrary Integer where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Int where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Int8 where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Int16 where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Int32 where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Int64 where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Word where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Word8 where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Word16 where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Word32 where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Word64 where
  coarbitrary = coarbitraryIntegral

instance CoArbitrary Char where
  coarbitrary = coarbitrary . ord

instance CoArbitrary Float where
  coarbitrary = coarbitraryReal

instance CoArbitrary Double where
  coarbitrary = coarbitraryReal

-- ** Helpers for implementing coarbitrary

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

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

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

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

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

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

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

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

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

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

max_stars_repo_name: ppetr/quickcheck

max_issues_repo_name: ppetr/quickcheck

max_forks_repo_name: ppetr/quickcheck


text:
module CategoricDefinitions where

import Prelude hiding (id, (.))
import qualified Prelude as P
import GHC.Exts
import Data.Kind (Type)
import Data.NumInstances.Tuple

import Numeric.LinearAlgebra.Array
import Numeric.LinearAlgebra.Array.Util

{-
Unfortunately, there is a known bug which prevents quantifying constraints involving type families.
https://gitlab.haskell.org/ghc/ghc/issues/14860
This results in all the functions below having a list of constraints AllAllowed k [...]
-}

-- Standard haskell way to define a category
class Category (k :: Type -> Type -> Type) where
    type Allowed k a :: Constraint
    type Allowed k a = ()

    id  :: Allowed k a => a `k` a
    (.) :: AllAllowed k [a, b, c] => b `k` c -> a `k` b -> a `k` c

-- By monoidal here we mean symmetric monoidal category
class Category k => Monoidal (k :: Type -> Type -> Type) where
    -- unit object in haskell is ()
    x :: AllAllowed k [a, b, c, d, (a, b), (c, d)]
      => (a `k` c) -> (b `k` d) -> ((a, b) `k` (c, d))
    assocL :: AllAllowed k [a, b, c, (a, b), ((a, b), c), (b, c), (a, (b, c))]
           => ((a, b), c) `k` (a, (b, c))
    assocR :: AllAllowed k [a, b, c, (a, b), ((a, b), c), (b, c), (a, (b, c))]
           => (a, (b, c)) `k` ((a, b), c)
    unitorL :: Allowed k a
            => ((), a) `k` a
    unitorL' :: Allowed k a -- inverse of unitorL
            => a `k` ((), a)
    -- could also potentially add right unitor and their inverses
    swap :: AllAllowed k [a, b, (a, b), (b, a)]
         => (a, b) `k` (b, a)

class Monoidal k => Cartesian k where
    type AllowedCar k a :: Constraint
    type AllowedCar k a = ()

    exl :: AllowedCar k b => (a, b) `k` a
    exr :: AllowedCar k a => (a, b) `k` b
    dup :: AllowedCar k a => a `k` (a, a)
    counit :: AllowedCar k a => a `k` ()

class Category k => Cocartesian k where
    type AllowedCoCar k a :: Constraint
    type AllowedCoCar k a = Allowed k a

    inl :: AllowedCoCar k b => a `k` (a, b)
    inr :: AllowedCoCar k a => b `k` (a, b)
    jam :: AllowedCoCar k a => (a, a) `k` a
    unit :: AllowedCoCar k a => () `k` a

{-
This is a hacky way of modelling a weak 2-category which is needed for Para.
Notice the tick' after class name
(.*) corresponds to id
(.-) corresponds to . (sequential comp)
(.|) corresponds to `x` (parallel comp)

Ideally, we'd like to specify in code that given any symmetric monoidal category C we can construct _another_ symmetric monoidal category called Para(C), but sometimes dreams will have to be dreams.
-}

class Category' (k :: Type -> Type -> Type -> Type) where
    type Allowed' k a :: Constraint
    type Allowed' k a = ()

    (.*) :: (Allowed' k a) => k () a a
    (.-) :: AllAllowed' k [p, q, a, b, c, (q, p), (p, q), ((p, q), a), ((q, p), a), (q, (p, a)), (q, b), (p, a)]
          => k q b c -> k p a b -> k (p, q) a c

-- the constraints are unfortunately very ugly
class Category' k => Monoidal' (k :: Type -> Type -> Type -> Type) where
    (.|) :: AllAllowed' k [a, b, c, d, p, q, ((p, q), (a, b)), ((p, a), (q, b)), (c, d), (p, a), (q, b), (p, (q, (a, b))), (p, ((q, a), b)), (p, ((a, q), b)), (p, (a, (q, b))), (a, (q, b)), ((a, q), b), (a, q), ((q, a), b), (q, a), (q, (a, b)), (a, b), (p, q)]
      => k p a c -> k q b d -> k (p, q) (a, b) (c, d)

-- Sequential composition of parametrized functions
(.--) :: (Monoidal k, _)
    => (q, b) `k` c
    -> (p, a) `k` b
    -> ((p, q), a) `k` c
g .-- f = g . (id `x` f) . assocL . (swap `x` id)

-- Parallel composition of parametrized functions
(.||) :: (Monoidal k, _)
    => (p, a) `k` b
    -> (q, c) `k` d
    -> ((p, q), (a, c)) `k` (b, d)
f .|| g = f `x` g . swapParam

{-
Swap map for monoidal product of parametrized functions, basically bracket bookkeeping.
Read from top to bottom
(a b) (c d)
a (b, (c, d))
a ((b, c), d)
a ((c, b), d)
a (c, (b, d))
(a c) (b d)
-}
swapParam :: (Monoidal k, _) => ((a, b), (c, d)) `k` ((a, c), (b, d))
swapParam = assocR . (id `x` assocL) . (id `x` (swap `x` id)) . (id `x` assocR) . assocL

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

class Additive a where
    zero :: a
    (^+) :: a -> a -> a

class NumCat (k :: Type -> Type -> Type) a where
    negateC :: a `k` a
    addC :: (a, a) `k` a
    mulC :: (a, a) `k` a
    increaseC :: a -> a `k` a -- curried add, add a single number

class FloatCat (k :: Type -> Type -> Type) a where
    expC :: a `k` a

class FractCat (k :: Type -> Type -> Type) a where
    recipC :: a `k` a

class Scalable (k :: Type -> Type -> Type) a where
    scale :: a -> (a `k` a)

type Tensor = NArray None Double

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

instance Category (->) where
    id    = \a -> a
    g . f = \a -> g (f a)

instance Monoidal (->) where
    f `x` g = \(a, b) -> (f a, g b)
    assocL = \((a, b), c) -> (a, (b, c))
    assocR = \(a, (b, c)) -> ((a, b), c)
    unitorL = \((), a) -> a
    unitorL' = \a -> ((), a)
    swap = \(a, b) -> (b, a)

instance Cartesian (->) where
    exl = \(a, _) -> a
    exr = \(_, b) -> b
    dup = \a -> (a, a)
    counit = \_ -> ()

instance Num a => NumCat (->) a where
    negateC = negate
    addC = uncurry (+)
    mulC = uncurry (*)
    increaseC a = (+a)

instance Floating a => FloatCat (->) a where
    expC = exp

instance Fractional a => FractCat (->) a where
    recipC = recip

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

instance Additive () where
    zero = ()
    () ^+ () = ()

instance {-# OVERLAPPABLE #-} Num a => Additive a where
    zero = 0
    (^+) = (+)

instance (Additive a, Additive b) => Additive (a, b) where
    zero = (zero, zero)
    (a1, b1) ^+ (a2, b2) = (a1 ^+ a2, b1 ^+ b2)

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

(/\) :: (Cartesian k, _) => b `k` c -> b `k` d -> b `k` (c, d)
f /\ g = (f `x` g) . dup

(\/) :: (Monoidal k, Cocartesian k, _) => a `k` c -> b `k` c -> (a, b) `k` c
f \/ g = jam . (f `x` g)

fork :: (Cartesian k, _) => (b `k` c,  b `k` d) -> b `k` (c, d)
fork (f, g) = f /\ g

unfork :: (Cartesian k, _) => b `k` (c, d) -> (b `k` c, b `k` d)
unfork h = (exl . h, exr . h)

join :: (Monoidal k, Cocartesian k, _) => (a `k` c, b `k` c) -> (a, b) `k` c
join (f, g) = f \/ g

unjoin :: (Cocartesian k, _) => (a, b) `k` c -> (a `k` c, b `k` c)
unjoin h = (h . inl, h . inr)

divide :: (Monoidal k, FractCat k a, _) => k (a, a) a
divide = mulC . (id `x` recipC)
-------------------------------------

type family AllAllowed k xs :: Constraint where
    AllAllowed k '[] = ()
    AllAllowed k (x : xs) = (Allowed k x, AllAllowed k xs)

type family AllAllowed' k xs :: Constraint where
    AllAllowed' k '[] = ()
    AllAllowed' k (x : xs) = (Allowed' k x, AllAllowed' k xs)


text:
module Sound (soundFFT, wavesToTime) where

import qualified Config       as Conf
import           Data.Complex
import           Data.WAVE
import           Numeric.FFT

-- | Perform Fast Fourier Transform on the input data, called the
-- "time" domain, in order to get its representation in the "frequency"
-- domain.
soundFFT :: [Double] -> [Double]
soundFFT timeDomain = take Conf.samplesPerFreq (drop Conf.samplesShift $ mfft $ map (\t -> t :+ 0) timeDomain)

mfft :: [Complex Double] -> [Double]
mfft p = map (\l -> realPart $ abs l) (fft p)

-- | Converts WAVESamples to [Double] using sampleToDouble
-- function. For each channel it takes the average.
wavesToTime :: WAVESamples -> [Double]
wavesToTime []     = []
wavesToTime (h:tl) = (average $ map sampleToDouble h) : wavesToTime tl

-- | Calculates the average value of a given list.
average :: [Double] -> Double
average list = (sum list) / realToFrac (length list)

max_stars_repo_name: adamszczerba/haskell-fft

max_issues_repo_name: adamszczerba/haskell-fft

max_forks_repo_name: adamszczerba/haskell-fft


text:
module PCA.PCA
       ( PCA (..)
       , covariance
       , eigenValues
       , eigenVectors
       , inputData
       , finalData
       , makePCA
       , restoredData
       , meanMatrix
       ) where

import Universum hiding (All, Vector, transpose)

import PCA.Types

import Control.Lens (makeLenses)
import Data.Array.Repa
import Data.Array.Repa.Repr.ForeignPtr
import Numeric.LinearAlgebra.Repa hiding (Matrix, Vector)

data PCA = PCA
  {   _inputData    :: Matrix D Double
    , _covariance   :: Covariance Double
    , _eigenVectors :: EigenVectors Double
    , _eigenValues  :: Vector F Double
    , _finalData    :: Matrix D Double
    , _restoredData :: Matrix D Double
    , _meanMatrix   :: Matrix D Double
  }

makeLenses ''PCA

makePCA
    :: Int
    -> Matrix D Double
    -> PCA
makePCA desiredDimensions input =
  let _inputData = input
         -- get dimension
      (Z :. yInp :. _) = extent input
      (meanVec, _covariance) = meanCovS input
      _meanMatrix = extend (Z :. yInp :. All) meanVec
      adjustInput'@(ADelayed (Z :. _ :. _) _) = input -^ _meanMatrix
      (_eigenValues, eigenVec) = eigSH _covariance
      eigenVecD@(ADelayed (Z :. y :. x) f) = transpose eigenVec
        -- leave only n desired eigenvectors
      eigenVectors' =
        if desiredDimensions >= y
        then eigenVecD
        else ADelayed (Z :. desiredDimensions :. x) f
        -- colunmns are eigenvectors
      _eigenVectors = transpose eigenVectors'
        -- data in new eigenvectors space
      finalData' = runIdentity $ eigenVectors' `mulP` transpose adjustInput'
      _finalData = transpose finalData'
        -- restore to the original space
      restoredDataWOMean = transpose $ pinvS eigenVectors' `mul` finalData'
        -- add mean
      _restoredData = restoredDataWOMean +^ _meanMatrix
  in PCA{..}

max_stars_repo_name: serokell/PCA

max_issues_repo_name: serokell/PCA

max_forks_repo_name: serokell/PCA


text:
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE RankNTypes #-}
module AI.Singularity.Data.Network where
import GHC.TypeLits
import qualified Numeric.LinearAlgebra.Data as LD
import qualified Numeric.LinearAlgebra.Devel as LDev
import qualified Numeric.LinearAlgebra as LA
import qualified Control.Monad.State as S
import Control.Lens
import Data.Proxy
import Data.List (foldl')
import Control.Monad.ST
import System.Random
import Control.Monad
import Control.Arrow
import qualified Data.Vector.Storable as VS

import AI.Singularity.Utils
import AI.Singularity.Data.Matrix
import AI.Singularity.Data.Vector

data GeneralNetwork vec mat inp out where
  FFLayer :: forall vec mat m n. (KnownNat m, KnownNat n) =>
    mat (m+1) n ->
    (forall a. Floating a => a -> a) ->
    GeneralNetwork vec mat (vec m) (vec n)
  FFSeq :: forall vec mat a b c.
    GeneralNetwork  vec mat a b ->
    GeneralNetwork  vec mat b c ->
    GeneralNetwork  vec mat a c
  FromFunc :: (Num a) =>
    (forall b. Num b => b -> b) ->
    GeneralNetwork vec mat a a
  FromBinFunc :: forall a vec mat. Fractional a =>
    (forall b. Fractional b => b -> b -> b) ->
    GeneralNetwork vec mat (a,a) a
  FFDiv :: forall vec mat a' a b b'.
    GeneralNetwork vec mat a b ->
    GeneralNetwork vec mat a' b' ->
    GeneralNetwork vec mat (a, a') (b, b')
  Split :: forall vec mat a b c. Num a =>
    GeneralNetwork vec mat a b ->
    GeneralNetwork vec mat a c ->
    GeneralNetwork vec mat a (b,c)
  SplitNet :: forall vec mat n k l c. (KnownNat n, KnownNat k, KnownNat l, (k + l) ~ n) =>
    GeneralNetwork vec mat (vec n) c ->
    GeneralNetwork vec mat (vec k, vec l) c
  Recurse :: forall vec mat a b. (Num b) =>
    GeneralNetwork vec mat (a,b) b ->
    GeneralNetwork vec mat [a] b
  First :: forall vec mat a b c. (Num c) =>
    GeneralNetwork vec mat a b ->
    GeneralNetwork vec mat (a,c) b
  Second :: forall vec mat a b c. (Num c) =>
    GeneralNetwork vec mat a b ->
    GeneralNetwork vec mat (c,a) b

type Network inp out = GeneralNetwork Vector Matrix inp out

type TrainedNetwork inp out = GeneralNetwork Vector (GeneralInTrainingMat LD.Matrix) inp out

type STNetwork s = GeneralNetwork Vector (STInTrainMat s)

changeMatType :: Applicative f => (forall m n. (KnownNat m, KnownNat n) => mat m n -> f (mat' m n)) ->
                                  GeneralNetwork vec mat inp out -> f (GeneralNetwork vec mat' inp out)
{-# INLINE changeMatType #-}
changeMatType conv (FFLayer m f)    = (\x -> FFLayer x f) <$> conv m
changeMatType conv (FFSeq n1 n2)    = FFSeq <$> changeMatType conv n1 <*> changeMatType conv n2
changeMatType conv (FFDiv n1 n2)    = FFDiv <$> changeMatType conv n1 <*> changeMatType conv n2
changeMatType conv (FromFunc f)     = pure $ FromFunc f
changeMatType conv (FromBinFunc f)  = pure $ FromBinFunc f
changeMatType conv (Split n1 n2)    = Split <$> changeMatType conv n1 <*> changeMatType conv n2
changeMatType conv (SplitNet n1)    = SplitNet <$> changeMatType conv n1
changeMatType conv (Recurse n)      = Recurse <$> changeMatType conv n
changeMatType conv (First n)        = First <$> changeMatType conv n
changeMatType conv (Second n)       = Second <$> changeMatType conv n

matToTrain :: (KnownNat m, KnownNat n) => Matrix m n -> Identity (InTrainMat m n)
matToTrain = pure . (\x -> GeneralInTrainingMat x 0 0)

saveMat :: Matrix m n -> S.State [LD.Matrix Float] (Matrix m n)
saveMat (Mat m) = do
  l <- S.get
  S.put (m:l)
  return (Mat m)

saveWeights = map LD.toLists . reverse . flip S.execState [] . changeMatType saveMat

loadMat :: Matrix m n -> S.State [LD.Matrix Float] (Matrix m n)
loadMat _ = do
  (a:l) <- S.get
  S.put l
  return (Mat a)

loadWeights w = flip S.evalState (map LD.fromLists w) . changeMatType loadMat

toTrain :: Network inp out -> TrainedNetwork inp out
toTrain = runIdentity . changeMatType matToTrain

fromTrain :: TrainedNetwork inp out -> Network inp out
fromTrain = runIdentity . changeMatType (Identity . view trMatr)

toSTTrain :: Network inp out -> ST s (STNetwork s inp out)
toSTTrain = changeMatType toSTInTrainM

instance Show (Network inp out) where
  show (FFLayer m _)     = "weights:\n" ++ show m
  show (FFSeq net1 net2) = show net1 ++ "\n" ++ show net2
  show (FromFunc _ )     = "FromFunc"
  show (FFDiv a b)       = show a ++ "\n\tdiv \n" ++ show b
  show (SplitNet n)      = "\tSplitNet \n" ++ show n
  show (Recurse n)       = "\tRec\n" ++ show n

createLayer :: forall m n. (KnownNat m, KnownNat n) => Matrix (m+1) n -> (forall a. Floating a => a -> a)  -> Network (Vector m) (Vector n)
createLayer weights activator = FFLayer weights activator

randomLayer :: forall m n k. (KnownNat m, KnownNat n, KnownNat k, k ~ (m + 1)) => (forall a. Floating a => a -> a) -> IO (Network (Vector m) (Vector n))
randomLayer activator = (\w -> createLayer w activator) <$> randomMat

identLayer :: Num a => Network a a
identLayer = FromFunc id

generalConductSignal :: forall mat a b. (forall m n. mat m n -> Matrix m n) -> a -> GeneralNetwork Vector mat a b -> b
generalConductSignal _ v (FromFunc f)        = f v
generalConductSignal _ v (FromBinFunc f)     = uncurry f v
generalConductSignal c v (FFSeq net1 net2)   = generalConductSignal c (generalConductSignal c v net1) net2
generalConductSignal c v (FFLayer weights f) = cmap f (c weights #> consV 1 v)
generalConductSignal c v (FFDiv n1 n2)       = (generalConductSignal c (fst v) n1, generalConductSignal c (snd v) n2)
generalConductSignal c v (Split n1 n2)       = (generalConductSignal c v n1, generalConductSignal c v n2)
generalConductSignal c v (SplitNet n)        = generalConductSignal c (uncurry appendV v) n
generalConductSignal c v (First n)           = generalConductSignal c (fst v) n
generalConductSignal c v (Second n)          = generalConductSignal c (snd v) n
generalConductSignal c v (Recurse n)         = foldl' (\b a -> generalConductSignal c (a,b) n) 0 v

conductSignalInTrain = generalConductSignal (view trMatr)

conductSignal :: forall a b. a -> Network a b -> b
conductSignal = generalConductSignal id

runNetwork :: Network a b -> a -> b
runNetwork = flip conductSignal

compose :: Network m n -> Network k m -> Network k n
compose net2 net1 = FFSeq net1 net2

sigmoid :: Floating a => a -> a
sigmoid = (1/) . (1+). exp . negate

softmax :: forall n. KnownNat n => Vector n -> Vector n
softmax = exp >>> id &&& (foldV (+) 0 >>> fromValV) >>> uncurry (/)

atFirst :: Num a => Network b c -> Network (b,a) (c,a)
atFirst = flip FFDiv identLayer
atSecond :: Num a => Network b c -> Network (a,b) (a,c)
atSecond = FFDiv identLayer

takeFst :: (Num c, Num a) => Network (a,c) a
takeFst = First identLayer
takeSnd :: (Num c, Num a) => Network (c,a) a
takeSnd = Second identLayer

prelstm :: forall m n k. (KnownNat m, KnownNat n, KnownNat k) => IO (Network [Vector m] (Vector n, Vector k))
prelstm = do
  let takeInp  = First identLayer
      takeOut  = Second . Second $ identLayer
      -- takeMem  = Second . First $ identLayer
      splitter = Split takeInp takeOut
  x <- SplitNet <$> randomLayer sigmoid
  y <- SplitNet <$> randomLayer sigmoid
  z <- SplitNet <$> randomLayer tanh
  k <- SplitNet <$> randomLayer sigmoid
  l <- randomLayer tanh
  let
    gateCell   = Split (FFSeq splitter x) (Second . First $ identLayer) `FFSeq` FromBinFunc (*)
    secondCell = Split gateCell (splitter `FFSeq` Split y z `FFSeq` FromBinFunc (*)) `FFSeq` FromBinFunc (+)
    thirdCell  = Split secondCell (splitter `FFSeq` k) `FFSeq` Split takeFst (atFirst l `FFSeq` FromBinFunc (*))
      -- get2 pre = Split pre (Split)
  return . Recurse $ thirdCell
-- MnistReader

gru :: forall m n. (KnownNat m, KnownNat n) => IO (Network [Vector m] (Vector n))
gru = do
  z <- SplitNet <$> randomLayer sigmoid
  r <- SplitNet <$> randomLayer sigmoid
  h <- SplitNet <$> randomLayer tanh
  let
    rxt = Split r takeSnd `FFSeq` FromBinFunc (*)
    ht' = Split rxt takeFst `FFSeq` h
    ht  = Split (z `FFSeq` Split (FromFunc (1-)) identLayer) (Split takeSnd ht') `FFSeq` FromBinFunc (*) `FFSeq` FromBinFunc (+)
  return $ Recurse ht
    --beg = Split (Split x takeFst `FFSeq` FromBinFunc (*)) takeFst `FFSeq` z

iden :: forall a. (Num a) => a -> a
iden x = x

max_stars_repo_name: Antystenes/Memetic-Predictor

max_issues_repo_name: Antystenes/Memetic-Predictor

max_forks_repo_name: Antystenes/Memetic-Predictor


text:
module Deconvolution where

import           Control.Monad             as M
import           Control.Monad.Parallel    as MP
import           Data.Array.Repa           as R
import           Data.Binary
import           Data.Complex
import           Data.List                 as L
import           DFT.Plan
import           Filter.Utils
import           FokkerPlanck.DomainChange
import           FokkerPlanck.MonteCarlo
import           FokkerPlanck.Pinwheel
import           Graphics.Gnuplot.Simple
import           Image.IO
import           STC
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Types
import           Utils.Array

main = do
  args@(numPointStr:numOrientationStr:numScaleStr:thetaSigmaStr:scaleSigmaStr:maxScaleStr:taoStr:numTrailStr:maxTrailStr:theta0FreqsStr:thetaFreqsStr:scale0FreqsStr:scaleFreqsStr:initDistStr:histFilePath:alphaStr:pinwheelFlagStr:numThreadStr:_) <-
    getArgs
  print args
  let numPoint = read numPointStr :: Int
      numOrientation = read numOrientationStr :: Int
      numScale = read numScaleStr :: Int
      thetaSigma = read thetaSigmaStr :: Double
      scaleSigma = read scaleSigmaStr :: Double
      maxScale = read maxScaleStr :: Double
      tao = read taoStr :: Double
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      theta0Freq = read theta0FreqsStr :: Double
      theta0Freqs = [-theta0Freq .. theta0Freq]
      thetaFreq = read thetaFreqsStr :: Double
      thetaFreqs = [-thetaFreq .. thetaFreq]
      initDist = read initDistStr :: [R2S1RPPoint]
      scale0Freq = read scale0FreqsStr :: Double
      scaleFreq = read scaleFreqsStr :: Double
      scale0Freqs = [-scale0Freq .. scale0Freq]
      scaleFreqs = [-scaleFreq .. scaleFreq]
      alpha = read alphaStr :: Double
      pinwheelFlag = read pinwheelFlagStr :: Bool
      numThread = read numThreadStr :: Int
      folderPath = "output/test/Deconvolution"
      (R2S1RPPoint (_, _, _, s)) = L.head initDist
  createDirectoryIfMissing True folderPath
  flag <- doesFileExist histFilePath
  radialArr <-
    if flag
      then R.map magnitude . getNormalizedHistogramArr <$>
           decodeFile histFilePath
      else do
        putStrLn "Couldn't find a Green's function data. Start simulation..."
        solveMonteCarloR2Z2T0S0Radial
          numThread
          numTrail
          maxTrail
          numPoint
          numPoint
          thetaSigma
          scaleSigma
          maxScale
          tao
          theta0Freqs
          thetaFreqs
          scale0Freqs
          scaleFreqs
          histFilePath
          (emptyHistogram
             [ (round . sqrt . fromIntegral $ 2 * (div numPoint 2) ^ 2)
             , L.length scale0Freqs
             , L.length theta0Freqs
             , L.length scaleFreqs
             , L.length thetaFreqs
             ]
             0)
  radialArrSink <-
    R.map magnitude . getNormalizedHistogramArr <$> decodeFile histFilePath
  arrR2Z2T0S0 <-
    computeUnboxedP $
    computeR2Z2T0S0ArrayRadial
      radialArr
      numPoint
      numPoint
      1
      maxScale
      thetaFreqs
      scaleFreqs
      theta0Freqs
      scale0Freqs
  arrR2Z2T0S0Sink <-
    computeUnboxedP $
    computeR2Z2T0S0ArrayRadial
      (cutoff 24 radialArrSink)
      numPoint
      numPoint
      1
      maxScale
      thetaFreqs
      scaleFreqs
      theta0Freqs
      scale0Freqs
  plan <- makeR2Z2T0S0Plan emptyPlan arrR2Z2T0S0
  sourceDistArr <-
    computeInitialDistributionR2T0S0
      plan
      numPoint
      numPoint
      theta0Freqs
      scale0Freqs
      maxScale
      initDist
  sourceDistArrR2S1RP <-
    r2z2Tor2s1rpP numOrientation thetaFreqs numScale scaleFreqs $ sourceDistArr
  let xIndex =
        [ (fromIntegral i) * 360 / (fromIntegral numOrientation)
        | i <- [0 .. numOrientation - 1]
        ] :: [Double]
      inputStyle =
        L.zipWith
          (\(R2S1RPPoint (x', y', _, _)) i ->
             let x = x' + center numPoint
                 y = y' + center numPoint
              in ( defaultStyle
                     { plotType = LinesPoints
                     , lineSpec =
                         CustomStyle
                           [LineTitle (printf "(%d,%d)" x y), PointType i]
                     }
                 , L.zip xIndex .
                   R.toList . R.slice (R.map magnitude sourceDistArrR2S1RP) $
                   (Z :. All :. (0 :: Int) :. x :. y)))
          initDist
          [1 ..]
  plotPathsStyle [PNG (folderPath </> "Input.png"), Title "Input"] inputStyle
  arrR2Z2T0S0F <-
    dftR2Z2T0S0 plan .
    computeS .
    makeFilter2D .
    R.traverse arrR2Z2T0S0
      -- (R.traverse (arrR2Z2T0S0) id $ \f idx@(Z :. _ :. _ :. _ :. _ :. a :. b) ->
      --    if a == center numPoint && b == center numPoint
      --      then 0
      --      else f idx)
      id $ \f (Z :. tf' :. sf' :. t0f :. s0f :. i :. j) ->
      let idx = (Z :. tf' :. sf' :. t0f :. s0f :. i :. j)
       in if i == center numPoint && j == center numPoint &&
             tf' == t0f
                    -- tf' == div (L.length thetaFreqs) 2 &&
                    -- sf' == div (L.length scaleFreqs) 2 &&
                    -- t0f == div (L.length theta0Freqs) 2 &&
                    -- s0f == div (L.length scale0Freqs) 2
            then f idx + (0.0 :+ 0)
            else f idx
  -- let arr' =
  --       r2z2t0s0Tor2s1rps1rp
  --         numOrientation
  --         thetaFreqs
  --         theta0Freqs
  --         numScale
  --         scaleFreqs
  --         scale0Freqs
  --         maxScale $
  --       arrR2Z2T0S0Sink
  --     arr'' =
  --       r2s1rps1rpTor2z2t0s0
  --         numOrientation
  --         thetaFreqs
  --         theta0Freqs
  --         numScale
  --         scaleFreqs
  --         scale0Freqs
  --         maxScale .
  --       computeUnboxedS .
  --       R.map
  --         (\x ->
  --            let (m, p) = polar x
  --             in if m == 0
  --                  then 0
  --                  else mkPolar (1 / m) p) $
  --       arr'
  arrR2Z2T0S0SinkF <-
    dftR2Z2T0S0 plan .
    computeS .
    makeFilter2D . computeSinkFromSourceR2Z2T0S0 thetaFreqs theta0Freqs
    $
    arrR2Z2T0S0Sink
  --   R.traverse
  --     (-- computeSinkFromSourceR2Z2T0S0 thetaFreqs theta0Freqs $
  --      arrR2Z2T0S0)
  --     id $ \f idx@(Z :. a :. b :. c :. d :. i :. j) ->
  --     let (m, p) =
  --           polar $
  --           conjugate
  --             (f (Z :. (L.length thetaFreqs - 1 - a) :.
  --                 (L.length scaleFreqs - 1 - b) :.
  --                 (L.length theta0Freqs - 1 - c) :.
  --                 (L.length scale0Freqs - 1 - d) :.
  --                 i :.
  --                 j))
  --      in if a == div (L.length thetaFreqs) 2 &&
  --            b == div (L.length scaleFreqs) 2 &&
  --            c == div (L.length theta0Freqs) 2 &&
  --            d == div (L.length scale0Freqs) 2
  --           then 1 + conjugate
  --                       (f (Z :. (L.length thetaFreqs - 1 - a) :.
  --                           (L.length scaleFreqs - 1 - b) :.
  --                           (L.length theta0Freqs - 1 - c) :.
  --                           (L.length scale0Freqs - 1 - d) :.
  --                           i :.
  --                           j))
  --           else conjugate
  --                  (f (Z :. (L.length thetaFreqs - 1 - a) :.
  --                      (L.length scaleFreqs - 1 - b) :.
  --                      (L.length theta0Freqs - 1 - c) :.
  --                      (L.length scale0Freqs - 1 - d) :.
  --                      i :.
  --                      j))
  -- Source field
  sourceArr <- convolveR2T0S0 plan arrR2Z2T0S0F sourceDistArr
  sourceR2Z2' <- R.sumP . R.sumS . rotateR2Z2T0S0Array $ sourceArr
  sourceR2S1RP <-
    r2z2Tor2s1rpP numOrientation thetaFreqs numScale scaleFreqs $ sourceR2Z2'
  let sourceR2Z2 -- = sourceR2Z2'
       =
        R.traverse2 sourceR2Z2' sourceDistArr const $ \f1 f2 idx@(Z :. a :. b :. c :. d) ->
          let x = div (L.length theta0Freqs) 2
              y = div (L.length scale0Freqs) 2
           in if f2 (Z :. x :. y :. c :. d) == 0
                then 0
                else f1 idx
      sourceStyle =
        L.zipWith
          (\(R2S1RPPoint (x', y', _, _)) i ->
             let x = x' + center numPoint
                 y = y' + center numPoint
              in ( defaultStyle
                     { plotType = LinesPoints
                     , lineSpec =
                         CustomStyle
                           [LineTitle (printf "(%d,%d)" x y), PointType i]
                     }
                 , L.zip xIndex .
                   R.toList . R.slice (R.map magnitude sourceR2S1RP) $
                   (Z :. All :. (0 :: Int) :. x :. y)))
          initDist
          [1 ..]
  plotPathsStyle [PNG (folderPath </> "Conv.png"), Title "Conv"] sourceStyle
  sourceField <-
    fmap (computeS . R.extend (Z :. (1 :: Int) :. All :. All)) .
    R.sumP .
    R.sumS .
    rotate4D .
    rotate4D . r2z2Tor2s1rp numOrientation thetaFreqs numScale scaleFreqs $
    sourceR2Z2
  plotImageRepaComplex (folderPath </> "Source.png") . ImageRepa 8 $ sourceField
  -- Deconvolution
  deconvSourceArr <- convolveR2T0S0 plan arrR2Z2T0S0SinkF . computeUnboxedS $ sourceR2Z2
  deconvSourceR2Z2 <-
    R.sumP . R.sumS . rotateR2Z2T0S0Array $ deconvSourceArr
  deconvSourceField <-
    fmap (computeS . R.extend (Z :. (1 :: Int) :. All :. All)) .
    R.sumP .
    R.sumS .
    rotate4D .
    rotate4D . r2z2Tor2s1rp numOrientation thetaFreqs numScale scaleFreqs $
    deconvSourceR2Z2
  plotImageRepaComplex (folderPath </> "DeconvSource.png") . ImageRepa 8 $
    deconvSourceField
  deconvSourceR2S1RP <-
    r2z2Tor2s1rpP numOrientation thetaFreqs numScale scaleFreqs $
    deconvSourceR2Z2
  let outputStyle =
        L.zipWith
          (\(R2S1RPPoint (x', y', _, _)) i ->
             let x = x' + center numPoint
                 y = y' + center numPoint
              in ( defaultStyle
                     { plotType = LinesPoints
                     , lineSpec =
                         CustomStyle
                           [LineTitle (printf "(%d,%d)" x y), PointType i]
                     }
                 , L.zip xIndex .
                   R.toList . R.slice (R.map magnitude deconvSourceR2S1RP) $
                   (Z :. All :. (0 :: Int) :. x :. y)
                     -- (R.sumS .
                   --    R.backpermute
                   --      (Z :. (L.length theta0Freqs) :. numPoint :. numPoint :.
                   --       (L.length scale0Freqs))
                   --      (\(Z :. a :. b :. c :. d) -> (Z :. a :. d :. b :. c)) .
                   --    R.map magnitude $
                   --    deconvSourceR2S1RP) $
                   -- (Z :. All :. x :. y)
                  ))
          initDist
          [1 ..]
  plotPathsStyle [PNG (folderPath </> "Output.png"), Title "Output"] outputStyle

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
module Neural.Util where

import Control.Monad
import Control.Monad.ST
import Control.Parallel.Strategies
import Data.Array.ST
import Data.STRef
import Numeric.LinearAlgebra
import System.Random

 -- shamelessly nicked from https://wiki.haskell.org/Random_shuffle

shuffle' :: [a] -> StdGen -> ([a],StdGen)
shuffle' xs gen = runST (do
        g <- newSTRef gen
        let randomRST lohi = do
              (a,s') <- liftM (randomR lohi) (readSTRef g)
              writeSTRef g s'
              return a
        ar <- newArray n xs
        xs' <- forM [1..n] $ \i -> do
                j <- randomRST (i,n)
                vi <- readArray ar i
                vj <- readArray ar j
                writeArray ar j vi
                return vj
        gen' <- readSTRef g
        return (xs',gen'))
  where
    n = length xs
    newArray :: Int -> [a] -> ST s (STArray s Int a)
    newArray n xs =  newListArray (1,n) xs

shuffleIO :: [a] -> IO [a]
shuffleIO xs = getStdRandom (shuffle' xs)

parZipWith :: Strategy c -> (a -> b -> c) -> [a] -> [b] -> [c]
parZipWith strat f x y = zipWith f x y `using` parList strat

parZipWith3 :: Strategy d -> (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
parZipWith3 strat f x y z = zipWith3 f x y z `using` parList strat

fullyFlatten :: [Matrix Double] -> Vector Double
fullyFlatten = vjoin . (parMap rdeepseq flatten)


text:
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE OverloadedStrings #-}

module MLP.DataSet (
   toPat,
   normPat
) where

import Numeric.LinearAlgebra.Data (cols, Matrix, R, takeColumns, dropColumns, toRows)
import Numeric.LinearAlgebra.Static (create, zipWithVector)
import MLP.Types (Pattern(..), Topology(..), zipWithPat)
import MLP.Network (scalarR)
import qualified Data.Vector.Storable as V
import qualified Data.Text as T
import Control.Monad (when)
import GHC.TypeLits (KnownNat)
import Data.Singletons ( SingI(sing), SingKind(fromSing) )
import Data.Maybe (mapMaybe)

toPat :: forall i o. (KnownNat i, KnownNat o) => Matrix R -> Either String [Pattern i o]
toPat m
   | cols m == (dimI + dimO) = Right $ zipWith Pattern is os
   | otherwise = Left "Mismatch on number of columns"

   where
      [dimI, dimO] = map fromIntegral . fromSing $ sing @'[i,o]
      is = mapMaybe create . toRows . takeColumns dimI $ m
      os = mapMaybe create . toRows . dropColumns dimI $ m

normPat :: forall i o. (KnownNat i, KnownNat o) => [Pattern i o] -> [Pattern i o]
normPat [] = []
normPat xxs@(x:xs) = map (\p -> (p - minp) / scalep) xxs
   where
      zeroPat = Pattern (scalarR 0) (scalarR 0)
      minp = foldl (zipWithPat min) x xs
      maxp = foldl (zipWithPat max) x xs
      scalep = zipWithPat max zeroPat (maxp - minp)

max_stars_repo_name: ovanr/Typed-MLP

max_issues_repo_name: ovanr/Typed-MLP

max_forks_repo_name: ovanr/Typed-MLP


text:
module UKF
  ( KF
  , kf_x
  , kf_p
  , kf_t

  , makeFilter
  , predict
  , update
  ) where

import Prelude hiding ((<>))

import App (Sensor (..))

import Numeric.LinearAlgebra

data KF = KF
  { kf_x :: Vector Double
  , kf_p :: Matrix Double
  , kf_t :: Word
  } deriving Show

makeFilter (Laser px py t) =
  KF { kf_x = 5 |> [px,
                    py,
                    0,
                    0,
                    0]
     , kf_p = (5><5) [0.15, 0.00, 0.0, 0.0, 0.0,
                      0.00, 0.15, 0.0, 0.0, 0.0,
                      0.00, 0.00, 1.0, 0.0, 0.0,
                      0.00, 0.00, 0.0, 1.0, 0.0,
                      0.00, 0.00, 0.0, 0.0, 1.0]
     , kf_t = t }

makeFilter (Radar rho phi rho' t) =
  KF { kf_x = 5 |> [rho * cos phi,
                    rho * sin phi,
                    rho',
                    0,
                    0]
     , kf_p = (5><5) [0.15, 0.00, 0.0, 0.0, 0.0,
                      0.00, 0.15, 0.0, 0.0, 0.0,
                      0.00, 0.00, 1.0, 0.0, 0.0,
                      0.00, 0.00, 0.0, 1.0, 0.0,
                      0.00, 0.00, 0.0, 0.0, 1.0]
     , kf_t = t }

-- Standard Deviation Tuning Parameters
std_a_     = 2.0
std_yawdd_ = 0.3

std_laspx_  = 0.15
std_laspy_  = 0.15

std_radr_   = 0.3
std_radphi_ = 0.03
std_radrd_  = 0.3

-- Unwind angle
normAngle θ
  | θ > pi  = normAngle (θ - 2*pi)
  | θ < -pi = normAngle (θ + 2*pi)
  | otherwise = θ

-- Calculate A st. m = AA^T
sqrtMat :: Matrix Double -> Matrix Double
sqrtMat = tr . chol . sym

{- constant turn rate and velocity magnitude (CTRV) motion model.

   x is the state augmented with the noise vector
   dt is a reasonable time step
 -}
process :: Vector Double -> Double -> Vector Double
process x dt =
  (5 |> [px, py, vel, yaw, yaw']) + f + v
  where
    -- State vector components
    px   = x ! 0
    py   = x ! 1
    vel  = x ! 2
    yaw  = x ! 3
    yaw' = x ! 4

    -- Noise vector components
    noiseA  = x ! 5
    noiseY' = x ! 6

    -- State model vector
    f = if abs yaw' > 0.001
        then
          5 |> [vel/yaw' * ((sin (yaw + yaw' * dt)) - (sin yaw)),
                vel/yaw' * ((cos yaw) - (cos (yaw + yaw' * dt))),
                0,
                yaw' * dt,
                0]
        else
          -- The special case where yaw rate is 0.
          -- Meaning that the motion is traveling
          -- along a straight line.
          5 |> [vel * dt * (cos yaw),
                vel * dt * (sin yaw),
                0,
                0,
                0]

    -- Noise model vector
    v = 5 |> [0.5 * noiseA * (dt ** 2) * (cos yaw),
              0.5 * noiseA * (dt ** 2) * (sin yaw),
              noiseA * dt,
              0.5 * noiseY' * (dt ** 2),
              noiseY' * dt]

{- This gives us KF(k+1 | k) for the filter.
 -}
predict :: Word -> KF -> KF
predict t kf =
  KF { kf_x = x'
     , kf_p = p'
     , kf_t = t }
  where
    -- Time inside the filter.
    filterTime =
      kf_t kf

    -- Delta-time. Essentially a timestep.
    dt =
      (fromIntegral (t - filterTime)) / 1000000.0

    -- n is the size of the state space with noise components
    n = 7

    -- l is a tuning parameter that controls the distance of the points.
    l = 3 - n

    -- Generate a list of augmented sigma points --
    sigmaPoints :: [Vector Double]
    sigmaPoints =
      let
        x = kf_x kf
        p = kf_p kf

        xAug = 7 |> [x ! 0,
                     x ! 1,
                     x ! 2,
                     x ! 3,
                     x ! 4,
                     0,
                     0]
        pAug = diagBlock [p, std_a_ ** 2, std_yawdd_ ** 2]

        sqrtP = toColumns $ sqrtMat $ (l+n) `scale` pAug
      in
        xAug : map (xAug +) sqrtP ++ map (xAug -) sqrtP

    -- Run the sigma points thru the motion model --
    predictedPoints =
      x0 : xs
      where
        x0 = process (head sigmaPoints) dt
        xs = map (\pnt -> process pnt dt)
                 (tail sigmaPoints)

    -- x(k+1 | k)
    x' :: Vector Double
    x' = foldr (\xn x -> x + (wn `scale` xn))
               (w0 `scale` x0)
               xs
      where
        x0 = head predictedPoints
        w0 = l / (l+n)
        wn = 1 / (2 * (l+n))
        xs = tail predictedPoints

    -- p(k+1 | k)
    p' :: Matrix Double
    p' = foldr (\xn p ->
                  let
                    xn_d  = xn - x'
                    xn_d' = 5 |> [xn_d ! 0,
                                  xn_d ! 1,
                                  xn_d ! 2,
                                  normAngle $ xn_d ! 3,
                                  xn_d ! 4]
                  in
                    p + (wn `scale` asColumn xn_d' <> asRow xn_d'))
               (w0 `scale` asColumn x0_d' <> asRow x0_d')
               xs
      where
        x0    = head predictedPoints
        w0    = l / (l+n)
        x0_d  = x0 - x'
        x0_d' = 5 |> [x0_d ! 0,
                      x0_d ! 1,
                      x0_d ! 2,
                      normAngle $ x0_d ! 3,
                      x0_d ! 4]

        wn = 1 / (2 * (l+n))
        xs = tail predictedPoints

{- For the case of laser measurments we don't need to
   use an unscented transformation for the measurment
   update. Standard kalman filter equations apply just
   fine here. However, we do need an unscented
   transformation for the model prediction. Because that
   is a non-linear model space. This gives us
   KF(k+1 | k+1) for the filter with a laser measurment.
 -}
update (Laser px py t) kf =
  KF { kf_x = x'
     , kf_p = p'
     , kf_t = t }
  where
    kf' = predict t kf

    -- Measurment transition matricies
    h = (2><5) [1, 0, 0, 0, 0,
                0, 1, 0, 0, 0]
    r = (2><2) [std_laspx_ ** 2, 0.00000000,
                0.00000000,      std_laspy_ ** 2]
    y = 2 |> [px, py] - (h #> kf_x kf')

    -- Kalman Filter Equations
    i  = ident 5
    ht = tr h
    s  = h <> kf_p kf' <> ht + r
    si = inv s
    k  = kf_p kf' <> ht <> si

    -- Estimation based on kalman filter gain kf(k+1 | k+1)
    x' = kf_x kf' + (k #> y)
    p' = (i - k <> h) <> kf_p kf'

{- For the case of radar measurments we do
   need to use an unscented transformation as
   the both the mapping to measurment space is
   non-linear. This gives us  KF(k+1 | k+1) for
   the filter
 -}
update (Radar rho phi rho' t) kf =
  kf -- skip
  where
    kf' = predict t kf

    n = 5
    l = 3 - n

    --Calculates the state space point as a radar space point.
    radarModel x =
      let
        _x  = x ! 0
        _y  = x ! 1
        _v  = x ! 2
        _θ  = x ! 3
        _vx = _v * (cos _θ)
        _vy = _v * (sin _θ)

        ρ  = sqrt ((_x ** 2) + (_y ** 2))
        φ  = atan2 _y _x
        ρ' = (_x*_vx + _y*_vy) / ρ
      in
        3 |> [ρ, φ, ρ']

    -- Generate a list of sigma points --
    sigmaPoints :: [Vector Double]
    sigmaPoints =
      let
        x = kf_x kf
        p = kf_p kf

        sqrtP = toColumns $ sqrtMat $ (l+n) `scale` p
      in
        x : map (x +) sqrtP ++ map (x -) sqrtP

    -- Run the sigma points thru the radar model --
    zs = map radarModel sigmaPoints

    -- z(k+1 | k)
    z' :: Vector Double
    z' = foldr (\zn z -> z + (wn `scale` zn))
               (w0 `scale` z0)
               (tail zs)
      where
        z0 = head zs
        w0 = l / (l+n)
        wn = 1 / (2 * (l+n))

    -- s(k+1 | k)
    s' :: Matrix Double
    s' = s + r
      where
        z0 = head zs
        w0 = l / (l+n)
        wn = 1 / (2 * (l+n))
        z0_d  = z0 - z'

        r = (3><3) [std_radr_ ** 2, 0, 0,
                    0, std_radphi_ ** 2, 0,
                    0, 0, std_radrd_ ** 2]

        s = foldr (\zn s ->
                     let
                       zn_d  = zn - z'
                     in
                       s + (wn `scale` asColumn zn_d <> asRow zn_d))
                  (w0 `scale` asColumn z0_d <> asRow z0_d)
                  (tail zs)


text:
{-# OPTIONS_HADDOCK hide #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}

{-# LANGUAGE CPP #-}
{-# LANGUAGE EmptyDataDecls  #-}
{-# LANGUAGE FlexibleInstances  #-}
{-# LANGUAGE ForeignFunctionInterface  #-}
{-# LANGUAGE FunctionalDependencies  #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ScopedTypeVariables  #-}

module Data.Eigen.Internal where

import Foreign.Ptr
import Foreign.ForeignPtr
import Foreign.Storable
import Foreign.C.Types
import Foreign.C.String
import Control.Monad
#if __GLASGOW_HASKELL__ >= 710
#else
import Control.Applicative
#endif
import System.IO.Unsafe
import Data.Complex
import Data.Bits
import Data.Binary
import Data.Binary.Get
import Data.Binary.Put
import qualified Data.Vector.Storable as VS
import qualified Data.ByteString as BS
import qualified Data.ByteString.Internal as BSI

class (Num a, Cast a b, Cast b a, Storable b, Code b) => Elem a b | a -> b where

instance Elem Float CFloat where
instance Elem Double CDouble where
instance Elem (Complex Float) (CComplex CFloat) where
instance Elem (Complex Double) (CComplex CDouble) where

class Cast a b where
    cast :: a -> b

instance Storable a => Binary (VS.Vector a) where
    put vs = put (BS.length bs) >> putByteString bs where
        (fp,fs) = VS.unsafeToForeignPtr0 vs
        es = sizeOf (VS.head vs)
        bs = BSI.fromForeignPtr (castForeignPtr fp) 0 (fs * es)

    get = get >>= getByteString >>= \bs -> let
        (fp,fo,fs) = BSI.toForeignPtr bs
        es = sizeOf (VS.head vs)
        vs = VS.unsafeFromForeignPtr0 (Data.Eigen.Internal.plusForeignPtr fp fo) (fs `div` es)
        in return vs

-- | Complex number for FFI with the same memory layout as std::complex\<T\>
data CComplex a = CComplex !a !a deriving Show

instance Storable a => Storable (CComplex a) where
    sizeOf _ = sizeOf (undefined :: a) * 2
    alignment _ = alignment (undefined :: a)
    poke p (CComplex x y) = do
        pokeElemOff (castPtr p) 0 x
        pokeElemOff (castPtr p) 1 y
    peek p = CComplex
        <$> peekElemOff (castPtr p) 0
        <*> peekElemOff (castPtr p) 1

data CTriplet a = CTriplet !CInt !CInt !a deriving Show

instance Storable a => Storable (CTriplet a) where
    sizeOf _ = sizeOf (undefined :: a) + sizeOf (undefined :: CInt) * 2
    alignment _ = alignment (undefined :: CInt)
    poke p (CTriplet row col val) = do
        pokeElemOff (castPtr p) 0 row
        pokeElemOff (castPtr p) 1 col
        pokeByteOff p (sizeOf (undefined :: CInt) * 2) val
    peek p = CTriplet
        <$> peekElemOff (castPtr p) 0
        <*> peekElemOff (castPtr p) 1
        <*> peekByteOff p (sizeOf (undefined :: CInt) * 2)

instance Cast CInt Int where; cast = fromIntegral
instance Cast Int CInt where; cast = fromIntegral
instance Cast CFloat Float where; cast (CFloat x) = x
instance Cast Float CFloat where; cast = CFloat
instance Cast CDouble Double where; cast (CDouble x) = x
instance Cast Double CDouble where; cast = CDouble
instance Cast (CComplex CFloat) (Complex Float) where; cast (CComplex x y) = cast x :+ cast y
instance Cast (Complex Float) (CComplex CFloat) where; cast (x :+ y) = CComplex (cast x) (cast y)
instance Cast (CComplex CDouble) (Complex Double) where; cast (CComplex x y) = cast x :+ cast y
instance Cast (Complex Double) (CComplex CDouble) where; cast (x :+ y) = CComplex (cast x) (cast y)

instance Cast a b => Cast (CTriplet a) (Int, Int, b) where; cast (CTriplet x y z) = (cast x, cast y, cast z)
instance Cast a b => Cast (Int, Int, a) (CTriplet b) where; cast (x,y,z) = CTriplet (cast x) (cast y) (cast z)

intSize :: Int
intSize = sizeOf (undefined :: CInt)

encodeInt :: CInt -> BS.ByteString
encodeInt x = BSI.unsafeCreate (sizeOf x) $ (`poke` x) . castPtr

decodeInt :: BS.ByteString -> CInt
decodeInt (BSI.PS fp fo fs)
    | fs == sizeOf x = x
    | otherwise = error "decodeInt: wrong buffer size"
    where x = performIO $ withForeignPtr fp $ peek . (`plusPtr` fo)

data CSparseMatrix a b
type CSparseMatrixPtr a b = Ptr (CSparseMatrix a b)

data CSolver a b
type CSolverPtr a b = Ptr (CSolver a b)

performIO :: IO a -> a
performIO = unsafeDupablePerformIO

plusForeignPtr :: ForeignPtr a -> Int -> ForeignPtr b
plusForeignPtr fp fo = castForeignPtr fp' where
    vs :: VS.Vector CChar
    vs = VS.unsafeFromForeignPtr (castForeignPtr fp) fo 0
    (fp', _) = VS.unsafeToForeignPtr0 vs

foreign import ccall "eigen-proxy.h free" c_freeString :: CString -> IO ()

call :: IO CString -> IO ()
call func = func >>= \c_str -> when (c_str /= nullPtr) $
    peekCString c_str >>= \str -> c_freeString c_str >> fail str

foreign import ccall "eigen-proxy.h free" free :: Ptr a -> IO ()

foreign import ccall "eigen-proxy.h eigen_setNbThreads" c_setNbThreads :: CInt -> IO ()
foreign import ccall "eigen-proxy.h eigen_getNbThreads" c_getNbThreads :: IO CInt

class Code a where; code :: a -> CInt
instance Code CFloat where; code _ = 0
instance Code CDouble where; code _ = 1
instance Code (CComplex CFloat) where; code _ = 2
instance Code (CComplex CDouble) where; code _ = 3

newtype MagicCode = MagicCode CInt deriving Eq

instance Binary MagicCode where
    put (MagicCode code) = putWord32be $ fromIntegral code
    get = MagicCode . fromIntegral <$> getWord32be

magicCode :: Code a => a -> MagicCode
magicCode x = MagicCode (code x `xor` 0x45696730)

#let api1 name, args = "foreign import ccall \"eigen_%s\" c_%s :: CInt -> %s\n%s :: forall b . Code b => %s\n%s = c_%s (code (undefined :: b))", #name, #name, args, #name, args, #name, #name

#api1 random,        "Ptr b -> CInt -> CInt -> IO CString"
#api1 identity,      "Ptr b -> CInt -> CInt -> IO CString"
#api1 add,           "Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"
#api1 sub,           "Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"
#api1 mul,           "Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"
#api1 diagonal,      "Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"
#api1 transpose,     "Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"
#api1 inverse,       "Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"
#api1 adjoint,       "Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"
#api1 conjugate,     "Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"
#api1 normalize,     "Ptr b -> CInt -> CInt -> IO CString"
#api1 sum,           "Ptr b -> Ptr b -> CInt -> CInt -> IO CString"
#api1 prod,          "Ptr b -> Ptr b -> CInt -> CInt -> IO CString"
#api1 mean,          "Ptr b -> Ptr b -> CInt -> CInt -> IO CString"
#api1 norm,          "Ptr b -> Ptr b -> CInt -> CInt -> IO CString"
#api1 trace,         "Ptr b -> Ptr b -> CInt -> CInt -> IO CString"
#api1 squaredNorm,   "Ptr b -> Ptr b -> CInt -> CInt -> IO CString"
#api1 blueNorm,      "Ptr b -> Ptr b -> CInt -> CInt -> IO CString"
#api1 hypotNorm,     "Ptr b -> Ptr b -> CInt -> CInt -> IO CString"
#api1 determinant,   "Ptr b -> Ptr b -> CInt -> CInt -> IO CString"
#api1 rank,          "CInt -> Ptr CInt -> Ptr b -> CInt -> CInt -> IO CString"
#api1 image,         "CInt -> Ptr (Ptr b) -> Ptr CInt -> Ptr CInt -> Ptr b -> CInt -> CInt -> IO CString"
#api1 kernel,        "CInt -> Ptr (Ptr b) -> Ptr CInt -> Ptr CInt -> Ptr b -> CInt -> CInt -> IO CString"
#api1 solve,         "CInt -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"
#api1 relativeError, "Ptr b -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"

#let api2 name, args = "foreign import ccall \"eigen_%s\" c_%s :: CInt -> %s\n%s :: forall a b . Code b => %s\n%s = c_%s (code (undefined :: b))", #name, #name, args, #name, args, #name, #name

#api2 sparse_new,           "CInt -> CInt -> Ptr (CSparseMatrixPtr a b) -> IO CString"
#api2 sparse_clone,         "CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"
#api2 sparse_fromList,      "CInt -> CInt -> Ptr (CTriplet b) -> CInt -> Ptr (CSparseMatrixPtr a b) -> IO CString"
#api2 sparse_toList,        "CSparseMatrixPtr a b -> Ptr (CTriplet b) -> CInt -> IO CString"
#api2 sparse_free,          "CSparseMatrixPtr a b -> IO CString"
#api2 sparse_makeCompressed,"CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"
#api2 sparse_uncompress,    "CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"
#api2 sparse_isCompressed,  "CSparseMatrixPtr a b -> Ptr CInt -> IO CString"
#api2 sparse_transpose,     "CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"
#api2 sparse_adjoint,       "CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"
#api2 sparse_pruned,        "CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"
#api2 sparse_prunedRef,     "CSparseMatrixPtr a b -> Ptr b -> Ptr (CSparseMatrixPtr a b) -> IO CString"
#api2 sparse_scale,         "CSparseMatrixPtr a b -> Ptr b -> Ptr (CSparseMatrixPtr a b) -> IO CString"
#api2 sparse_nonZeros,      "CSparseMatrixPtr a b -> Ptr CInt -> IO CString"
#api2 sparse_innerSize,     "CSparseMatrixPtr a b -> Ptr CInt -> IO CString"
#api2 sparse_outerSize,     "CSparseMatrixPtr a b -> Ptr CInt -> IO CString"
#api2 sparse_coeff,         "CSparseMatrixPtr a b -> CInt -> CInt -> Ptr b -> IO CString"
#api2 sparse_coeffRef,      "CSparseMatrixPtr a b -> CInt -> CInt -> Ptr (Ptr b) -> IO CString"
#api2 sparse_cols,          "CSparseMatrixPtr a b -> Ptr CInt -> IO CString"
#api2 sparse_rows,          "CSparseMatrixPtr a b -> Ptr CInt -> IO CString"
#api2 sparse_norm,          "CSparseMatrixPtr a b -> Ptr b -> IO CString"
#api2 sparse_squaredNorm,   "CSparseMatrixPtr a b -> Ptr b -> IO CString"
#api2 sparse_blueNorm,      "CSparseMatrixPtr a b -> Ptr b -> IO CString"
#api2 sparse_add,           "CSparseMatrixPtr a b -> CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"
#api2 sparse_sub,           "CSparseMatrixPtr a b -> CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"
#api2 sparse_mul,           "CSparseMatrixPtr a b -> CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"
#api2 sparse_block,         "CSparseMatrixPtr a b -> CInt -> CInt -> CInt -> CInt -> Ptr (CSparseMatrixPtr a b) -> IO CString"
#api2 sparse_fromMatrix,    "Ptr b -> CInt -> CInt -> Ptr (CSparseMatrixPtr a b) -> IO CString"
#api2 sparse_toMatrix,      "CSparseMatrixPtr a b -> Ptr b -> CInt -> CInt -> IO CString"
#api2 sparse_values,        "CSparseMatrixPtr a b -> Ptr CInt -> Ptr (Ptr b) -> IO CString"
#api2 sparse_outerStarts,   "CSparseMatrixPtr a b -> Ptr CInt -> Ptr (Ptr CInt) -> IO CString"
#api2 sparse_innerIndices,  "CSparseMatrixPtr a b -> Ptr CInt -> Ptr (Ptr CInt) -> IO CString"
#api2 sparse_innerNNZs,     "CSparseMatrixPtr a b -> Ptr CInt -> Ptr (Ptr CInt) -> IO CString"
#api2 sparse_setZero,       "CSparseMatrixPtr a b -> IO CString"
#api2 sparse_setIdentity,   "CSparseMatrixPtr a b -> IO CString"
#api2 sparse_reserve,       "CSparseMatrixPtr a b -> CInt -> IO CString"
#api2 sparse_resize,        "CSparseMatrixPtr a b -> CInt -> CInt -> IO CString"

#api2 sparse_conservativeResize,    "CSparseMatrixPtr a b -> CInt -> CInt -> IO CString"
#api2 sparse_compressInplace,       "CSparseMatrixPtr a b -> IO CString"
#api2 sparse_uncompressInplace,     "CSparseMatrixPtr a b -> IO CString"

#let api3 name, args = "foreign import ccall \"eigen_%s\" c_%s :: CInt -> CInt -> %s\n%s :: forall s a b . (Code s, Code b) => s -> %s\n%s s = c_%s (code (undefined :: b)) (code s)", #name, #name, args, #name, args, #name, #name

#api3 sparse_la_newSolver,          "Ptr (CSolverPtr a b) -> IO CString"
#api3 sparse_la_freeSolver,         "CSolverPtr a b -> IO CString"
#api3 sparse_la_factorize,          "CSolverPtr a b -> CSparseMatrixPtr a b -> IO CString"
#api3 sparse_la_analyzePattern,     "CSolverPtr a b -> CSparseMatrixPtr a b -> IO CString"
#api3 sparse_la_compute,            "CSolverPtr a b -> CSparseMatrixPtr a b -> IO CString"
#api3 sparse_la_tolerance,          "CSolverPtr a b -> Ptr CDouble -> IO CString"
#api3 sparse_la_setTolerance,       "CSolverPtr a b -> CDouble -> IO CString"
#api3 sparse_la_maxIterations,      "CSolverPtr a b -> Ptr CInt -> IO CString"
#api3 sparse_la_setMaxIterations,   "CSolverPtr a b -> CInt -> IO CString"
#api3 sparse_la_info,               "CSolverPtr a b -> Ptr CInt -> IO CString"
#api3 sparse_la_error,              "CSolverPtr a b -> Ptr CDouble -> IO CString"
#api3 sparse_la_iterations,         "CSolverPtr a b -> Ptr CInt -> IO CString"
#api3 sparse_la_solve,              "CSolverPtr a b -> CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"
-- #api3 sparse_la_solveWithGuess,     "CSolverPtr a b -> CSparseMatrixPtr a b -> CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"
#api3 sparse_la_matrixQ,            "CSolverPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"
#api3 sparse_la_matrixR,            "CSolverPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"
#api3 sparse_la_setPivotThreshold,  "CSolverPtr a b -> CDouble -> IO CString"
#api3 sparse_la_rank,               "CSolverPtr a b -> Ptr CInt -> IO CString"
#api3 sparse_la_matrixL,            "CSolverPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"
#api3 sparse_la_matrixU,            "CSolverPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"
#api3 sparse_la_setSymmetric,       "CSolverPtr a b -> CInt -> IO CString"
#api3 sparse_la_determinant,        "CSolverPtr a b -> Ptr b -> IO CString"
#api3 sparse_la_logAbsDeterminant,  "CSolverPtr a b -> Ptr b -> IO CString"
#api3 sparse_la_absDeterminant,     "CSolverPtr a b -> Ptr b -> IO CString"
#api3 sparse_la_signDeterminant,    "CSolverPtr a b -> Ptr b -> IO CString"

ext: hsc

max_stars_repo_name: osidorkin/haskell-eigen

max_issues_repo_name: osidorkin/haskell-eigen

max_forks_repo_name: osidorkin/haskell-eigen


text:
{-# language DataKinds, ScopedTypeVariables #-}
module Main where

import System.Environment
import GHC.TypeLits
import Data.Proxy
import Data.List
import Data.List.Split (chunksOf)
import Data.Maybe
import Numeric.LinearAlgebra (toList)

import GAPoly

type Args = (String, Int, Int, Int, Integer, Double)

validateArgs :: [String] -> Args
validateArgs args =
  case args of
       (fname:nGens:nPop:nTerms:maxK:pm:_) -> (fname, read nGens, read nPop, read nTerms, read maxK, read pm)
       _ -> error "missing arguments"

term2str :: [(Bool,Int)] -> String
term2str ts = intercalate "*" $ filter (/="") $ zipWith toStr ts [0..]
  where
    toStr (False,_) _  = ""
    toStr (True, k) ix = "x[:," ++ show ix ++ "]**(" ++ show k ++ ")"

poly2str :: Int -> Solution (Poly n) -> String
poly2str nVars (Sol cs betas _) = let cs'    = map _getPoly cs
                                      terms  = map term2str $ chunksOf nVars cs'
                                      (b:bs) = reverse $ toList $ fromJust betas
                                      bts    = intercalate " + " $ zipWith toStr bs terms
                                   in show b ++ " + " ++ bts
  where
    toStr beta term = show beta ++ "*" ++ term

main :: IO ()
main = do args <- getArgs
          let (dataname, nGens, nPop, nTerms, maxK, pm) = validateArgs args
          dat <- readFile dataname
          let (xss, ys, nVars) = parseFile dat
          case someNatVal maxK of
               Just (SomeNat (_ :: Proxy n)) ->
                 do (avgs, best) <- runGA nTerms nVars nGens nPop pm xss ys
                    print $ poly2str nVars (best :: Solution (Poly n))
               Nothing -> error "Negative maxK"

max_stars_repo_name: folivetti/gapoly

max_issues_repo_name: folivetti/gapoly

max_forks_repo_name: folivetti/gapoly