text:
{-# LANGUAGE CPP                       #-}
{-# LANGUAGE DataKinds                 #-}
{-# LANGUAGE DeriveAnyClass            #-}
{-# LANGUAGE DeriveGeneric             #-}
{-# LANGUAGE ExistentialQuantification #-}
{-# LANGUAGE FlexibleContexts          #-}
{-# LANGUAGE FlexibleInstances         #-}
{-# LANGUAGE GADTs                     #-}
{-# LANGUAGE KindSignatures            #-}
{-# LANGUAGE RankNTypes                #-}
{-# LANGUAGE ScopedTypeVariables       #-}
{-# LANGUAGE TypeOperators             #-}
{-# LANGUAGE UndecidableInstances      #-}

{-|
Module      : Grenade.Dynamic.Specification
Description : Specifying network architectures for dynamically creating nets.
Copyright   : (c) Manuel Schneckenreither, 2020
License     : BSD2
Stability   : experimental
-}
module Grenade.Dynamic.Specification
  ( FromDynamicLayer (..)
  , ToDynamicLayer (..)
  , SpecNetwork (..)
  , SpecNet (..)
  , GNum (..)
  -- Layer specific data types
  , Dimensions
  , SpecFullyConnected(..)
  , SpecConvolution(..)
  , SpecDeconvolution(..)
  , SpecDropout(..)
  , SpecElu(..)
  , SpecGelu(..)
  , SpecLeakyRelu(..)
  , SpecLogit(..)
  , SpecRelu(..)
  , SpecReshape(..)
  , SpecSinusoid(..)
  , SpecSoftmax(..)
  , SpecTanh(..)
  , SpecTrivial(..)
  ) where

import           Control.DeepSeq
import           Data.Constraint                       (Dict (..))
import           Data.Proxy
import           Data.Reflection                       (reifyNat)
import           Data.Serialize
import           Data.Singletons
import           Data.Singletons.Prelude.List
import           Data.Typeable                         as T
import           GHC.Generics                          hiding (R)
import           GHC.TypeLits
import           Unsafe.Coerce                         (unsafeCoerce)
#if MIN_VERSION_base(4,9,0)
import           Data.Kind                             (Type)
#endif
#ifndef FLYCHECK
import {-# SOURCE #-} Grenade.Dynamic.Layers.Convolution    ()
import {-# SOURCE #-} Grenade.Dynamic.Layers.Deconvolution  ()
import {-# SOURCE #-} Grenade.Dynamic.Layers.Dropout        ()
import {-# SOURCE #-} Grenade.Dynamic.Layers.Elu            ()
import {-# SOURCE #-} Grenade.Dynamic.Layers.FullyConnected ()
import {-# SOURCE #-} Grenade.Dynamic.Layers.Gelu           ()
import {-# SOURCE #-} Grenade.Dynamic.Layers.LeakyRelu      ()
import {-# SOURCE #-} Grenade.Dynamic.Layers.Logit          ()
import {-# SOURCE #-} Grenade.Dynamic.Layers.Relu           ()
import {-# SOURCE #-} Grenade.Dynamic.Layers.Reshape        ()
import {-# SOURCE #-} Grenade.Dynamic.Layers.Sinusoid       ()
import {-# SOURCE #-} Grenade.Dynamic.Layers.Softmax        ()
import {-# SOURCE #-} Grenade.Dynamic.Layers.Tanh           ()
import {-# SOURCE #-} Grenade.Dynamic.Layers.Trivial        ()
#endif

import           Control.Monad.Primitive               (PrimBase, PrimState)
import           Control.Parallel.Strategies
import           Numeric.LinearAlgebra.Static          (L, R)
import qualified Numeric.LinearAlgebra.Static          as H
import           System.Random.MWC                     (Gen)

import           Grenade.Core
import           Grenade.Types
import           Grenade.Utils.ListStore

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

-- | Grenade Num class.
--
-- This allows for instance scalar multiplication of the weights, which is useful for slowly adapting networks, e.g. NN'
-- <- \tau * NN' + (1-\tau) * NN. Or one could sum up some gradients in parallel and apply them at once: @applyUpdate lp
-- net $ foldl1 (|+) ...@aq.
class GNum a where
  (|*) :: Rational -> a -> a
  (|+) :: a -> a -> a
  gFromRational :: Rational -> a

infixl 7 |*
infixr 5 |+

instance (SingI i) => GNum (Network '[] '[ i]) where
  _ |* NNil = NNil
  _ |+ NNil = NNil
  gFromRational _ = NNil

instance (SingI i, SingI o, Layer x i o, NFData (Network xs (o ': rs)), GNum x, GNum (Network xs (o ': rs))) => GNum (Network (x ': xs) (i ': o ': rs)) where
  s |* (x :~> xs) =
    let x' = (s |* x)
        xs' = (s |* xs) `using` rdeepseq
     in x' :~> xs'
  (x :~> xs) |+ (y :~> ys) =
    let x' = (x |+ y)
        xs' = (xs |+ ys) `using` rdeepseq
     in x' :~> xs'
  gFromRational r = gFromRational r :~> gFromRational r

instance GNum (Gradients '[]) where
  _ |* GNil = GNil
  _ |+ GNil = GNil
  gFromRational _ = GNil

instance (GNum a) => GNum [a] where
  r |* xs = fmap (r |*) xs
  xs |+ ys = zipWith (|+) xs ys
  gFromRational _ = error "Cannot create a list of elements using gFromRational"

instance (UpdateLayer x, GNum (Gradient x), GNum (Gradients xs), NFData (Gradients xs)) => GNum (Gradients (x ': xs)) where
  s |* (x :/> xs) =
    let x' = (s |* x)
        xs' = (s |* xs) `using` rdeepseq
     in x' :/> xs'
  (x :/> xs) |+ (y :/> ys) =
    let x' = (x |+ y)
        xs' = (xs |+ ys) `using` rdeepseq
     in x' :/> xs'
  gFromRational r = gFromRational r :/> gFromRational r

instance GNum () where
  _ |* () = ()
  _ |+ () = ()
  gFromRational _ = ()

instance (GNum a, GNum b) => GNum (a, b) where
  s |* (a, b) = (s |* a, s |* b)
  (a1, b1) |+ (a2, b2) = (a1 |+ a2, b1 |+ b2)
  gFromRational v = (gFromRational v, gFromRational v)

instance (KnownNat m) => GNum (R m) where
  s |* vec = H.dvmap (fromRational s *) vec
  (|+) = (+)
  gFromRational = fromRational

instance (KnownNat m, KnownNat n) => GNum (L m n) where
  s |* mat = H.dmmap (fromRational s *) mat
  (|+) = (+)
  gFromRational = fromRational

instance (GNum a) => GNum (ListStore a) where
  s |* (ListStore n xs) = ListStore n $ map (fmap (s |*)) xs
  (ListStore n1 xs) |+ (ListStore n2 ys) = ListStore (max n1 n2) $ zipWith (\x y -> (|+) <$> x <*> y) xs ys
  gFromRational _ = mkListStore

-- | Create a runtime dynamic specification of a network. Dynamic layers (and networks), for storing and restoring specific network structures (e.g. in saving the network structures to a DB and
-- restoring it from there) or simply generating them at runtime. This does not store the weights and biases! They have to be handled separately (see Serialize)!
class FromDynamicLayer x where
  fromDynamicLayer :: SomeSing Shape -> SomeSing Shape -> x -> SpecNet

-- | Class for generating layers from a specification. See the function @networkFromSpecification@.
class (Show spec) => ToDynamicLayer spec where
  toDynamicLayer :: (PrimBase m) => WeightInitMethod -> Gen (PrimState m) -> spec -> m SpecNetwork

----------------------------------------
-- Return value of toDynamicLayer

-- | Specification of a network or layer.
data SpecNetwork :: Type where
  SpecNetwork
    :: ( FromDynamicLayer (Network layers shapes)
       , FoldableGradient (Gradients layers)
       , GNum (Gradients layers)
       , GNum (Network layers shapes)
       , Layer (Network layers shapes) (Head shapes) (Last shapes)
       , NFData (Gradients layers)
       , NFData (Network layers shapes)
       , NFData (Tapes layers shapes)
       , RandomLayer (Network layers shapes)
       , Serialize (Network layers shapes)
       , Serialize (Gradients layers)
       , Show (Network layers shapes)
       , SingI (Head shapes)
       , SingI (Last shapes)
       , SingI shapes
       , Typeable (Head shapes)
       , Typeable layers
       , Typeable shapes
       , CreatableNetwork layers shapes
       , RunnableNetwork layers shapes
       )
    => !(Network layers shapes)
    -> SpecNetwork
  SpecLayer
    :: ( FromDynamicLayer x
       , FoldableGradient (Gradient x)
       , GNum (Gradient x)
       , GNum x
       , Layer x i o
       , NFData (Gradient x)
       , NFData (Tape x i o)
       , NFData x
       , RandomLayer x
       , Serialize x
       , Serialize (Gradient x)
       , Show x
       , Typeable i
       , Typeable x
       , CreatableNetwork '[x] '[i,o]
       )
    => !x
    -> !(Sing i)
    -> !(Sing o)
    -> SpecNetwork

instance Show SpecNetwork where
  show (SpecNetwork net) = show net
  show (SpecLayer x _ _) = show x

instance FromDynamicLayer SpecNetwork where
  fromDynamicLayer inp out (SpecNetwork net) = fromDynamicLayer inp out net
  fromDynamicLayer inp out (SpecLayer x _ _) = fromDynamicLayer inp out x

-- | Data structure for holding specifications for networks. Networks can be built dynamically with @toDynamicLayer@. Further see the functions @|=>@, @specNil1D@, @specNil2D@, @specNil3D@, and
-- possibly any other layer implementation of @ToDynamicLayer@ for building specifications. See functions @networkFromSpecification@ to create a network from a specification and
-- @networkToSpecification@ to create a specification from a network.
data SpecNet
  = SpecNNil1D !Integer                   -- ^ 1D network output
  | SpecNNil2D !Integer !Integer          -- ^ 2D network output
  | SpecNNil3D !Integer !Integer !Integer -- ^ 3D network output
  | SpecNCons !SpecNet !SpecNet           -- ^ x :~> xs, where x can also be a network
  | forall spec . (ToDynamicLayer spec, Typeable spec, Ord spec, Eq spec, Show spec, Serialize spec, NFData spec) => SpecNetLayer !spec -- ^ Specification of a layer

instance Eq SpecNet where
  SpecNNil1D r1 == SpecNNil1D r2 = r1 == r2
  SpecNNil2D r1 c1 == SpecNNil2D r2 c2 = r1 == r2 && c1 == c2
  SpecNNil3D r1 c1 d1 == SpecNNil3D r2 c2 d2 = r1 == r2 && c1 == c2 && d1 == d2
  SpecNCons l1 r1 == SpecNCons l2 r2 = l1 == l2 && r1 == r2
  SpecNetLayer spec1 == SpecNetLayer spec2 = T.typeOf spec1 == T.typeOf spec2 && Just spec1 == cast spec2
  _ == _ = False

instance Ord SpecNet where
  SpecNNil1D r1 `compare` SpecNNil1D r2 = r1 `compare` r2
  SpecNNil2D r1 c1 `compare` SpecNNil2D r2 c2 = compare (r1, c1) (r2, c2)
  SpecNNil3D r1 c1 d1 `compare` SpecNNil3D r2 c2 d2 = compare (r1, c1, d1) (r2, c2, d2)
  SpecNCons l1 r1 `compare` SpecNCons l2 r2 = compare (l1, r1) (l2, r2)
  SpecNetLayer spec1 `compare` SpecNetLayer spec2
    | T.typeOf spec1 == T.typeOf spec2 = compare (Just spec1) (cast spec2)
    | otherwise = LT
  x `compare` y = compare x y

instance NFData SpecNet where
  rnf (SpecNNil1D r1)       = rnf r1
  rnf (SpecNNil2D r1 c1)    = rnf r1 `seq` rnf c1
  rnf (SpecNNil3D r1 c1 d1) = rnf r1 `seq` rnf c1 `seq` rnf d1
  rnf (SpecNCons l1 r1)     = rnf l1 `seq` rnf r1
  rnf (SpecNetLayer spec1)  = rnf spec1

instance Serialize SpecNet where
  put (SpecNNil1D r1)       = put (0 :: Int) >> put r1
  put (SpecNNil2D r1 c1)    = put (1 :: Int) >> put r1 >> put c1
  put (SpecNNil3D r1 c1 d1) = put (2 :: Int) >> put r1 >> put c1 >> put d1
  put (SpecNCons l1 r1)     = put (3 :: Int) >> put l1 >> put r1
  put (SpecNetLayer spec1)  = put (4 :: Int) >> put (show $ typeOf spec1) >> put spec1
  get = do
    (nr :: Int) <- get
    case nr of
      0 -> SpecNNil1D <$> get
      1 -> SpecNNil2D <$> get <*> get
      2 -> SpecNNil3D <$> get <*> get <*> get
      3 -> SpecNCons <$> get <*> get
      4 -> get >>= parseSpecDataConstructor
      _ -> error "unexpected input in Serialize instance of DynamicNetwork"

instance Show SpecNet where
  show (SpecNNil1D o) = "SpecNNil1D " ++ show o
  show (SpecNNil2D r c) = "SpecNNil2D " ++ show r ++ "x" ++ show c
  show (SpecNNil3D x y z) = "SpecNNil3D " ++ show x ++ "x" ++ show y ++ "x" ++ show z
  show (SpecNCons layer rest) = case layer of
    SpecNCons{} -> "[ " ++ show layer ++ " ] :=> " ++ show rest
    _           -> show layer ++ " :=> " ++ show rest
  show (SpecNetLayer layer) = show layer

-- | Function to parse the known implemented types. Can only deserialise these!
parseSpecDataConstructor :: String -> Get SpecNet
#ifdef FLYCHECK
parseSpecDataConstructor = error "disable -DFLYCHECK while compiling!"
#else
parseSpecDataConstructor repStr
  | repStr == show (typeRep (Proxy :: Proxy SpecFullyConnected)) = SpecNetLayer <$> (get :: Get SpecFullyConnected)
  | repStr == show (typeRep (Proxy :: Proxy SpecConvolution))    = SpecNetLayer <$> (get :: Get SpecConvolution)
  | repStr == show (typeRep (Proxy :: Proxy SpecDeconvolution))  = SpecNetLayer <$> (get :: Get SpecDeconvolution)
  | repStr == show (typeRep (Proxy :: Proxy SpecDropout))        = SpecNetLayer <$> (get :: Get SpecDropout)
  | repStr == show (typeRep (Proxy :: Proxy SpecElu))            = SpecNetLayer <$> (get :: Get SpecElu)
  | repStr == show (typeRep (Proxy :: Proxy SpecGelu))           = SpecNetLayer <$> (get :: Get SpecGelu)
  | repStr == show (typeRep (Proxy :: Proxy SpecLeakyRelu))      = SpecNetLayer <$> (get :: Get SpecLeakyRelu)
  | repStr == show (typeRep (Proxy :: Proxy SpecLogit))          = SpecNetLayer <$> (get :: Get SpecLogit)
  | repStr == show (typeRep (Proxy :: Proxy SpecRelu))           = SpecNetLayer <$> (get :: Get SpecRelu)
  | repStr == show (typeRep (Proxy :: Proxy SpecReshape))        = SpecNetLayer <$> (get :: Get SpecReshape)
  | repStr == show (typeRep (Proxy :: Proxy SpecSinusoid))       = SpecNetLayer <$> (get :: Get SpecSinusoid)
  | repStr == show (typeRep (Proxy :: Proxy SpecSoftmax))        = SpecNetLayer <$> (get :: Get SpecSoftmax)
  | repStr == show (typeRep (Proxy :: Proxy SpecTanh))           = SpecNetLayer <$> (get :: Get SpecTanh)
  | repStr == show (typeRep (Proxy :: Proxy SpecTrivial))        = SpecNetLayer <$> (get :: Get SpecTrivial)
  | repStr == show (typeRep (Proxy :: Proxy SpecNet))            = SpecNetLayer <$> (get :: Get SpecNet)
  | otherwise = error $ "unexpected input parseSpecDataConstructor: " ++ repStr
#endif

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

-- Network instances

instance (KnownNat rows) => FromDynamicLayer (Network '[] '[ 'D1 rows ]) where
  fromDynamicLayer _ _ _ = SpecNNil1D nr
    where nr = natVal (Proxy :: Proxy rows)
instance (KnownNat rows, KnownNat cols) => FromDynamicLayer (Network '[] '[ 'D2 rows cols ]) where
  fromDynamicLayer _ _ _ = SpecNNil2D (natVal (Proxy :: Proxy rows)) (natVal (Proxy :: Proxy cols))

instance (KnownNat rows, KnownNat cols, KnownNat depth) => FromDynamicLayer (Network '[] '[ 'D3 rows cols depth ]) where
  fromDynamicLayer _ _ _ = SpecNNil3D (natVal (Proxy :: Proxy rows)) (natVal (Proxy :: Proxy cols)) (natVal (Proxy :: Proxy depth))

instance (FromDynamicLayer x, FromDynamicLayer (Network xs (h : rs)), SingI h) => FromDynamicLayer (Network (x ': xs) (i ': h ': rs)) where
  fromDynamicLayer inp out ((x :: x) :~> (xs :: Network xs (h ': rs))) = SpecNCons (fromDynamicLayer inp hShape x) (fromDynamicLayer hShape out xs)
    where
      hShape = SomeSing (sing :: Sing h)

instance ToDynamicLayer SpecNet where
  toDynamicLayer _ _ (SpecNNil1D nrOut) =
    reifyNat nrOut $ \(_ :: (KnownNat o') => Proxy o') ->
      case (unsafeCoerce (Dict :: Dict ()) :: Dict (layers ~ '[]), unsafeCoerce (Dict :: Dict ()) :: Dict (shapes ~ '[ 'D1 o'])) of
        (Dict, Dict) -> return $ SpecNetwork (NNil :: Network '[] '[ 'D1 o'])

  toDynamicLayer _ _ (SpecNNil2D rows cols) =
    reifyNat rows $ \(_ :: (KnownNat r) => Proxy r) -> reifyNat cols $ \(_ :: (KnownNat c) => Proxy c) ->
      case (unsafeCoerce (Dict :: Dict ()) :: Dict (layers ~ '[]), unsafeCoerce (Dict :: Dict ()) :: Dict (shapes ~ '[ 'D2 r c ])) of
        (Dict, Dict) -> return $ SpecNetwork (NNil :: Network '[] '[ 'D2 r c])

  toDynamicLayer _ _ (SpecNNil3D rows cols depth) =
    reifyNat rows $ \(_ :: (KnownNat r) => Proxy r) -> reifyNat cols $ \(_ :: (KnownNat c) => Proxy c) -> reifyNat depth $ \(_ :: (KnownNat d) => Proxy d) ->
      case ( unsafeCoerce (Dict :: Dict ()) :: Dict (layers ~ '[])
           , unsafeCoerce (Dict :: Dict ()) :: Dict (shapes ~ '[ 'D3 r c d ])
           , unsafeCoerce (Dict :: Dict ()) :: Dict (KnownNat (r GHC.TypeLits.* d))) of
        (Dict, Dict, Dict) -> return $ SpecNetwork (NNil :: Network '[] '[ 'D3 r c d])

  toDynamicLayer wInit gen (SpecNetLayer layer) = toDynamicLayer wInit gen layer

  toDynamicLayer wInit gen (SpecNCons left right) = do
    restNet <- toDynamicLayer wInit gen right
    case restNet of
      SpecLayer{} -> error "Expected a network, but received a layer in Specification"
      SpecNetwork (xs :: Network restLayers restShapes) ->
        case (sing :: Sing restShapes) of
          SNil -> error "unexpected empty network (shapes where SNil) in Specification"
          (SCons (h :: Sing h) (hs :: Sing hs)) -> withSingI h $ withSingI hs $ do
            layer <- toDynamicLayer wInit gen left
            case layer of
              SpecLayer (x :: xType) singIn singOut ->
                case (singIn, singOut) of
                  (pxIn :: Sing inShape, pxOut :: Sing outShape) ->
                    withSingI pxIn $ withSingI h $ withSingI pxOut $
                     case (unsafeCoerce (Dict :: Dict ()) :: Dict (outShape ~ h)) of
                       Dict -> return $ SpecNetwork (x :~> xs :: Network (xType ': restLayers) (inShape ': restShapes))
              SpecNetwork (x :: Network xLayers xShapes) ->
                case (sing :: Sing xShapes) of
                  SNil -> error "unexpected empty network (SpecNNil) as layer in specification. Cannot proceed."
                  SCons (i :: Sing i) (xHs :: Sing xHs) -> withSingI i $ withSingI xHs $
                    case ( unsafeCoerce (Dict :: Dict ()) :: Dict (Head restShapes ~ Last xShapes)) of
                      Dict -> return $ SpecNetwork (x :~> xs :: Network (Network xLayers xShapes ': restLayers) (i ': restShapes))

-- Data structures stances for Layers (needs to be defined here)

type Dimensions = (Integer, Integer, Integer)

-- | Data Structure to save a fully connected layer. Saves number of input and output nodes.
data SpecFullyConnected = SpecFullyConnected !Integer !Integer
  deriving (Show, Eq, Ord, Serialize, Generic, NFData)

-- data SpecConcat = SpecConcat SpecNet SpecNet
--   deriving (Show, Eq, Ord, Serialize, Generic, NFData)

-- | Specifiation of a convolutional layer. Saves input and channels, filters, kernelRows, kernelColumns, strideRows, strideColumns, and kernelFlattened.
data SpecConvolution =
  SpecConvolution !Dimensions !Integer !Integer !Integer !Integer !Integer !Integer
  deriving (Show, Read, Eq, Ord, Serialize, Generic, NFData)

-- | Specifiation of a deconvolutional layer. Saves input and channels, filters, kernelRows, kernelColumns, strideRows, strideColumns, and kernelFlattened.
data SpecDeconvolution =
  SpecDeconvolution !Dimensions !Integer !Integer !Integer !Integer !Integer !Integer
  deriving (Show, Read, Eq, Ord, Serialize, Generic, NFData)

-- | Specification of a dropout filter. Saves number of input nodes, ratio and maybe a seed.
data SpecDropout = SpecDropout !Integer !RealNum !(Maybe Int)
  deriving (Show, Read, Eq, Ord, Serialize, Generic, NFData)

-- | Specification of Elu, saves input dimensions as triple, where every element >= 1.
newtype SpecElu = SpecElu Dimensions
  deriving (Show, Read, Eq, Ord, Serialize, Generic, NFData)

-- | Specification of Leaky Relu (with 0.01), saves input dimensions as triple, where every element >= 1.
newtype SpecGelu = SpecGelu Dimensions
  deriving (Show, Read, Eq, Ord, Serialize, Generic, NFData)

-- | Specification of Leaky Relu, saves alpha input dimensions as triple, where every element >= 1.
data SpecLeakyRelu = SpecLeakyRelu Dimensions RealNum
  deriving (Show, Read, Eq, Ord, Serialize, Generic, NFData)

-- | Specification of Logit, saves input dimensions as triple, where every element >= 1.
newtype SpecLogit = SpecLogit Dimensions
  deriving (Show, Read, Eq, Ord, Serialize, Generic, NFData)

-- | Specification of Relu, saves input dimensions as triple, where every element >= 1.
newtype SpecRelu = SpecRelu Dimensions
  deriving (Show, Read, Eq, Ord, Serialize, Generic, NFData)

-- | Specification of a reshape layer, saves input and output dimensions, where every element >= 1.
data SpecReshape = SpecReshape !Dimensions !Dimensions
  deriving (Show, Read, Eq, Ord, Serialize, Generic, NFData)

-- | Specification of Sinusoid, saves input dimensions as triple, where every element >= 1.
newtype SpecSinusoid = SpecSinusoid Dimensions
  deriving (Show, Read, Eq, Ord, Serialize, Generic, NFData)

-- | Specification of Softmax, saves input dimensions as triple, where every element >= 1.
newtype SpecSoftmax = SpecSoftmax Integer
  deriving (Show, Read, Eq, Ord, Serialize, Generic, NFData)

-- | Specification of Tanh, saves input dimensions as triple, where every element >= 1.
newtype SpecTanh = SpecTanh Dimensions
  deriving (Show, Read, Eq, Ord, Serialize, Generic, NFData)

-- | Specification of Trivial, saves input dimensions as triple, where every element >= 1.
newtype SpecTrivial = SpecTrivial Dimensions
  deriving (Show, Read, Eq, Ord, Serialize, Generic, NFData)

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-# LANGUAGE UnicodeSyntax, DataKinds, TypeOperators, KindSignatures,
             TypeInType, GADTs, MultiParamTypeClasses, FunctionalDependencies,
             TypeFamilies, AllowAmbiguousTypes, FlexibleInstances,
             UndecidableInstances, InstanceSigs, TypeApplications,
             ScopedTypeVariables, EmptyCase, FlexibleContexts,
             RankNTypes, LambdaCase
#-}

-- Implementation of density matrix interpretation of quantum computation
module LinTrans2 where

import Prelim

import qualified Data.Complex as C
import Data.Complex (Complex(..))
import Control.Monad.State.Lazy
import Data.Singletons
import Data.Maybe (fromJust)
--import Data.Tuple hiding (swap)
import Unsafe.Coerce
import Data.Constraint (Dict(..))
import Data.List ( (\\) )
import qualified Debug.Trace as Tr

-- I don't like the show instance for ℂ
newtype ℂ = ℂ (Complex Double)

liftC :: (Complex Double -> Complex Double) -> ℂ -> ℂ
liftC f (ℂ c) = ℂ $ f c

i :: ℂ
i = ℂ $ 0 :+ 1
conjugate = liftC C.conjugate
instance Num ℂ where
  ℂ m + ℂ n = ℂ $ m + n
  ℂ m - ℂ n = ℂ $ m - n
  ℂ m * ℂ n = ℂ $ m * n
  abs = liftC abs
  signum = liftC signum
  fromInteger = ℂ . fromInteger
instance Fractional ℂ where
  fromRational = ℂ . fromRational
  ℂ m / ℂ n    = ℂ $ m / n
instance Floating ℂ where
  pi    = ℂ $ pi
  exp   = liftC exp
  log   = liftC log
  sin   = liftC sin
  cos   = liftC cos
  asin  = liftC asin
  acos  = liftC acos
  atan  = liftC atan
  sinh  = liftC sinh
  cosh  = liftC cosh
  asinh = liftC asinh
  acosh = liftC acosh
  atanh = liftC atanh
instance Show ℂ where
  show (ℂ (α :+ β)) =   if β == 0 then show α
                        else if α == 0 then show β ++ "i"
                        else show α ++ " + " ++ show β ++ "i"

--------------
-- Matrices --
--------------

-- A matrix Matrix (m,n) f of dimension 2^m × 2^n is a function f from pairs of
-- nats to ℂ
data Matrix = Matrix { dim :: (Int,Int)
                     , idx :: (Int, Int) -> ℂ
                     }
dom :: Matrix -> [Int]
dom (Matrix (m,_) _) = [0..(2^m)-1]

cod :: Matrix -> [Int]
cod (Matrix (_,n) _) = [0..(2^n)-1]

mkMatrix :: (Int,Int) -> [[ℂ]] -> Matrix
mkMatrix (m,n) ls = Matrix (m,n) $ \(i,j) -> (ls !! i) !! j

--------------
-- Printing --
--------------

rows :: Matrix -> [[ℂ]]
rows mat = f <$> dom mat
  where
    f i = (\j -> idx mat (i,j)) <$> cod mat
showRows :: Matrix -> [String]
showRows mat = show <$> rows mat
instance Show Matrix where
  show mat = unlines $ showRows mat

---------------------------------------
-- Primitive matrices and operations --
---------------------------------------

ident :: Int -> Matrix
ident n = Matrix (n,n) $ \(i,j) -> if i == j then 1 else 0

transpose :: Matrix -> Matrix
transpose (Matrix (m,n) f) = Matrix (n,m) $ \(i,j) -> f (j,i)

dot :: Matrix -> Matrix -> ℂ
dot v1 v2 | snd (dim v1) == fst (dim v2) =
            foldl (+) 0 $ (\x -> idx v1 (0,x) * idx v2 (x,0)) <$> dom v2
          | otherwise = error "Mismatched dimensions when taking dot product"

row :: Int -> Matrix -> Matrix
row i (Matrix (_,n) f) = Matrix (0,n) $ \(_,j) -> f(i,j)

col :: Int -> Matrix -> Matrix
col j (Matrix (m,_) f) = Matrix (m,0) $ \(i,_) -> f(i,j)

plusM :: Matrix -> Matrix -> Matrix
plusM (Matrix dim1 f1) (Matrix dim2 f2) | dim1 == dim2 =
       Matrix dim1 $ \(i,j) -> f1(i,j) + f2(i,j)
plusM _ _ | otherwise    = error "Mismatched dimensions when adding matrices"

multM :: Matrix -> Matrix -> Matrix
multM mat1 mat2 | n1 == m2 =
    Matrix (m1,n2) $ \(i,j) -> dot (row i mat1) (col j mat2)
  where
    (m1,n1) = dim mat1
    (m2,n2) = dim mat2
multM _ _ | otherwise = error "Mismatched dimensions when multiplying matrices"

kron :: Matrix -> Matrix -> Matrix
kron (Matrix (m1,n1) f1) (Matrix (m2,n2) f2) =
    Matrix (m1*n1, m2*n2) $ \(i,j) ->
           f1 (i `div` 2^m2, j `div` 2^n2) * f2 (i `mod` 2^m2, j `mod` 2^n2)

zero n = Matrix (n,n) $ \_ -> 0

scale :: ℂ -> Matrix -> Matrix
scale c (Matrix (m,n) f) = Matrix (m,n) $ \(i,j) -> c * f(i,j)

----------------------
-- Quantum Matrices --
----------------------

ket0 = mkMatrix (1,0) [[1],[0]]
ket1 = mkMatrix (1,0) [[0],[1]]
density0 = ket0 * transpose ket0
density1 = ket1 * transpose ket1
newD True  = density1
newD False = density0

hadamard = mkMatrix (1,1) [[1/sqrt 2, 1/sqrt 2]
                          ,[1/sqrt 2, -1/sqrt 2]]

pauliX = mkMatrix (1,1) [[0,1]
                       ,[1,0]]
pauliY = mkMatrix (1,1) [[0,-1]
                       ,[i,0]]
pauliZ = mkMatrix (1,1) [[1,0]
                       ,[0,-1]]

cnot = mkMatrix (2,2) [[1,0,0,0]
                      ,[0,1,0,0]
                      ,[0,0,0,1]
                      ,[0,0,1,0]]

------------------
-- Permutations --
------------------

-- encode and decode modulo a key k
encode :: Int -> [Int] -> Int
encode k [] = 0
encode k (i : ls) = i + k * encode k ls

decode :: Int -> Int -> [Int]
decode k i | i == 0 = []
decode k i | otherwise = i `mod` k : decode k (i `div` k)

fromAssocList :: [(Int,Int)] -> Int -> Int
fromAssocList []           = id
fromAssocList ((a,b) : ls) = assocFun a b . fromAssocList ls
  where
    assocFun a b i = if i == a then b else i

-- Check if two lists are equal up to the permutation specified, i.e.
-- if permuting ls1 by f is equal to ls2
isPermutation :: Eq a => (Int -> Int) -> [a] -> [a] -> Bool
isPermutation f ls1 ls2 = all (\i -> ls1 !! i == ls2 !! f i) [0..len]
  where
    len = length ls1

swapFun :: Int -> (Int -> Int) -> Matrix
swapFun k f = Matrix (k,k) $ \(i,j) ->
    if isPermutation f (decode k i) (decode k j) then 1 else 0

-- swap takes a list of qubit/bit variables and produces a matrix that permutes
-- those qubits.
-- i.e. swap ls |φ_0 ⋯ φ_n⟩ = |φ_f(0) ⋯ ρ_f(n)⟩
-- where f = fromAssocList ls
swap :: Int -> [Int] -> Matrix
swap k ls = swapFun k $ fromAssocList $ zip [0..] ls

-------------------
-- Density Monad --
-------------------

-- A density monad is a nondeterminism state monad:
-- Density -> [Density]
-- An element op of type DensityMonad corresponds to the superoperator
-- \ρ -> ∑ (op ρ)
type DensityMonad = StateT Matrix []

newM :: Bool -> DensityMonad Int
newM b = do
    ρ ← get
    put $ ρ `kron` newD b
    return $ fst (dim ρ)

-- runQ applies the superoperator to the identity density matrix of size 1.
runQ :: DensityMonad a -> [(a,Matrix)]
runQ m = runStateT m (ident 0)

-- getDensity combines the result of runQ into a single density matrix
getDensity :: DensityMonad a -> Matrix
getDensity m = foldr f (zero n) ls
  where
    ls = runQ m
    n = fst . dim . snd $ head ls
    f (_,ρ) ρ0 = ρ + ρ0

-----------------
-- Application --
-----------------

super :: Matrix -> DensityMonad ()
super mat = do ρ ← get
               put $ mat * ρ * transpose mat

applyMatrix :: Matrix -> [Int] -> DensityMonad ()
applyMatrix mat ls = do super $ swap n ls  -- moves the relevant qubits to the front
                        super mat -- applies operation
                        super . transpose $ swap n ls -- moves the qubits back to where they were
  where
    (n,_) = dim mat

branch :: DensityMonad a -> DensityMonad a -> DensityMonad a
branch m0 m1 = StateT $ \ρ -> [runStateT m0 ρ, runStateT m1 ρ] >>= id

measM :: Int -> DensityMonad Bool
measM i = branch (applyMatrix density0 [i] >> return False) (applyMatrix density1 [i] >> return True)

------------------
-- Num instance --
------------------

instance Num Matrix where
  (+) = plusM
  mat1 - mat2 = mat1 `plusM` scale (-1) mat2
  (*) = multM
  abs = undefined
  signum mat = undefined -- trace?
  fromInteger n = ident $ fromInteger n

-------------------
-- Show Instance --
-------------------

instance Show (DensityMonad a) where
  show = show . getDensity

{-

------------------
-- Num instance --
------------------

instance (SingI m, SingI n) => Num (Matrix m n) where
  (+) mat1 mat2 (i,j) = mat1(i,j) + mat2(i,j)
  (-) mat1 mat2 (i,j) = mat1(i,j) - mat2(i,j)
  (*) mat1 mat2 (i,j) = mat1(i,j) * mat2(i,j)
  abs mat (i,j) = abs $ mat(i,j)
  signum mat = undefined -- trace mat
  fromInteger n = matrix . repeat $ fromIntegral n
instance Num Density where
  Density n1 mat1 + Density n2 mat2 =
    withSingI (two `raiseToSNat` n1) $
    withSingI (two `raiseToSNat` n2) $
    case eqSNat n1 n2 of
      Left Dict -> Density n1 $ mat1 + mat2
      Right _   -> error $ "Cannot add mismatched matrices "
                            ++ show mat1 ++ " and " ++ show mat2
  _ - _ = undefined
  _ * _ = undefined
  abs _ = undefined
  signum = undefined
  fromInteger n = undefined

-- Debugging help
m0 = square @Two [1,0,0,0
                 ,0,1,0,0
                 ,0,0,0,1
                 ,0,0,1,0]
rho = Density three $ \case
  (0,0) -> -0.5
  (0,2) ->  0.5
  (3,0) -> -0.5
  (3,2) ->  0.5
  (_,_) -> 0

-}

ext: hs

max_stars_repo_name: jpaykin/LNLHaskell

max_issues_repo_name: jpaykin/LNLHaskell

max_forks_repo_name: jpaykin/LNLHaskell


text:
module Number
    (
    parseNumber
    ) where

import Text.ParserCombinators.Parsec hiding (spaces)
import Control.Monad
import Numeric
import Data.Ratio
import Data.Complex

import LispType

-- Parse Numbers
-- |
-- >>> readExpr "123"
-- "Found value"
parseNumber :: Parser LispVal
parseNumber = parseDecimal
          <|> parseHex
          <|> parseOct
          <|> parseBin
          <|> parseRatio
          <|> parseComplex

parseDecimal = parseDecimal1 <|> parseDecimal2

parseDecimal1 :: Parser LispVal
parseDecimal1 = liftM (Number . read) $ many1 digit

parseDecimal2 :: Parser LispVal
parseDecimal2 = do
  _ <- try $ string "#d"
  x <- many1 digit
  return $ Number $ read x

parseHex :: Parser LispVal
parseHex = do
  _ <- try $ string "#x"
  x <- many1 hexDigit
  return $ Number $ hex2dec x

-- |
-- >>> readExpr "#o12345"
-- "Found value"
parseOct :: Parser LispVal
parseOct = do
  _ <- try $ string "#o"
  x <- many1 octDigit
  return $ Number $ oct2dec x

-- |
-- >>> readExpr "#b1010"
-- "Found value"
parseBin :: Parser LispVal
parseBin = do
  _ <- try $ string "#b"
  x <- many1 $ oneOf "01"
  return $ Number $ bin2dec x

hex2dec :: String -> Integer
hex2dec x = fst $ head $ readHex x

oct2dec :: String -> Integer
oct2dec x = fst $ head $ readOct x

bin2dec :: String -> Integer
bin2dec = bin2dec' 0
  where
    bin2dec' digint "" = digint
    bin2dec' digint (x:xs) =
      let old = 2 * digint + (if x == '0' then 0 else 1) in
        bin2dec' old xs

-- |
-- >>> readExpr "10.123"
-- "Found value"
parseFloat :: Parser LispVal
parseFloat = do
  x <- many1 digit
  _ <- char '.'
  y <- many1 digit
  return $ Float $ fst . head $ readFloat $ x ++ "." ++ y

-- |
-- >> readExpr "/ 3 4"
-- "Found value"
parseRatio :: Parser LispVal
parseRatio = do
  x <- many1 digit
  _ <- char '/'
  y <- many1 digit
  return $ Ratio $ read x % read y

{- https://hackage.haskell.org/package/base-4.8.1.0/docs/Data-Complex.html -}

-- |
-- >> readExpr "44.2+33i"
-- "Found value"

{- on lisp
> (+ 12+1i 12+1i)
24+2i
>
-}

toDouble :: LispVal -> Double
toDouble (Float f)  = realToFrac f
toDouble (Number n) = realToFrac n

parseComplex :: Parser LispVal
parseComplex = do
  x <- try parseFloat <|> parseDecimal
  _ <- char '+'
  y <- try parseFloat <|> parseDecimal
  _ <- char 'i'
  return $ Complex (toDouble x :+ toDouble y)

max_stars_repo_name: bokuo-okubo/scheme-of-bko

max_issues_repo_name: bokuo-okubo/scheme-of-bko

max_forks_repo_name: bokuo-okubo/scheme-of-bko


text:
{-# LANGUAGE BangPatterns          #-}
{-# LANGUAGE DeriveFoldable        #-}
{-# LANGUAGE DeriveTraversable     #-}
{-# LANGUAGE EmptyCase             #-}
{-# LANGUAGE EmptyDataDeriving     #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE LambdaCase            #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE StandaloneDeriving    #-}
{-# LANGUAGE TemplateHaskell       #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeInType            #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE UndecidableInstances  #-}
{-# OPTIONS_GHC -fno-warn-orphans  #-}

-- |
-- Module      : Data.Mutable.Instances
-- Copyright   : (c) Justin Le 2020
-- License     : BSD3
--
-- Maintainer  : justin@jle.im
-- Stability   : experimental
-- Portability : non-portable
--
-- Exports 'Ref' data types for various common data types, and also the
-- tools for automatic derivation of instances.  See "Data.Mutable" for
-- more information.
module Data.Mutable.Instances (
    RecRef(..)
  , HListRef(..)
  , UnitRef(..)
  , VoidRef
  -- * Generic
  , GRef(..)
  , gThawRef, gFreezeRef
  , gCopyRef, gMoveRef, gCloneRef
  , gUnsafeThawRef, gUnsafeFreezeRef
  , GMutable (GRef_)
  -- * Higher-Kinded Data Pattern
  , thawHKD, freezeHKD
  , copyHKD, moveHKD, cloneHKD
  , unsafeThawHKD, unsafeFreezeHKD
  -- * Coercible
  , CoerceRef(..)
  , thawCoerce, freezeCoerce
  , copyCoerce, moveCoerce, cloneCoerce
  , unsafeThawCoerce, unsafeFreezeCoerce
  -- * Traversable
  , TraverseRef(..)
  , thawTraverse, freezeTraverse
  , copyTraverse, moveTraverse, cloneTraverse
  , unsafeThawTraverse, unsafeFreezeTraverse
  -- * Immutable
  , ImmutableRef(..), thawImmutable, freezeImmutable, copyImmutable
  -- * Instances for Generics combinators themselves
  , GMutableRef(..)
  , MutSumF(..)
  -- * Utility
  , MapRef
  ) where

import           Control.Applicative
import           Data.Complex
import           Data.Functor.Compose
import           Data.Functor.Identity
import           Data.Functor.Product
import           Data.Functor.Sum
import           Data.Generics.Product.Internal.HList      (HList(..))
import           Data.Kind
import           Data.Mutable.Internal
import           Data.Mutable.Internal.TH
import           Data.Ord
import           Data.Primitive.Array
import           Data.Primitive.ByteArray
import           Data.Primitive.PrimArray
import           Data.Primitive.SmallArray
import           Data.Primitive.Types
import           Data.Ratio
import           Data.Vinyl                    as V hiding (HList)
import           Data.Void
import           Data.Word
import           Foreign.C.Types
import           Foreign.Storable
import           Numeric.Natural
import qualified Data.Monoid                               as M
import qualified Data.Vector                               as V
import qualified Data.Vector.Generic                       as VG
import qualified Data.Vector.Generic.Mutable               as MVG
import qualified Data.Vector.Mutable                       as MV
import qualified Data.Vector.Primitive                     as VP
import qualified Data.Vector.Primitive.Mutable             as MVP
import qualified Data.Vector.Storable                      as VS
import qualified Data.Vector.Storable.Mutable              as MVS
import qualified Data.Vector.Unboxed                       as VU
import qualified Data.Vector.Unboxed.Mutable               as MVU
import qualified Data.Vinyl.ARec                           as V
import qualified Data.Vinyl.Functor                        as V
import qualified Data.Vinyl.TypeLevel                      as V

instance Mutable s Int
instance Mutable s Integer
instance Mutable s Natural
instance Mutable s (Ratio a)
instance Mutable s Float
instance Mutable s Double
instance Mutable s (Complex a)
instance Mutable s Bool
instance Mutable s Char

instance Mutable s Word
instance Mutable s Word8
instance Mutable s Word16
instance Mutable s Word64

instance Mutable s CChar
instance Mutable s CSChar
instance Mutable s CUChar
instance Mutable s CShort
instance Mutable s CUShort
instance Mutable s CInt
instance Mutable s CUInt
instance Mutable s CLong
instance Mutable s CULong
instance Mutable s CPtrdiff
instance Mutable s CSize
instance Mutable s CWchar
instance Mutable s CSigAtomic
instance Mutable s CLLong
instance Mutable s CULLong
instance Mutable s CBool
instance Mutable s CIntPtr
instance Mutable s CUIntPtr
instance Mutable s CIntMax
instance Mutable s CUIntMax
instance Mutable s CClock
instance Mutable s CTime
instance Mutable s CUSeconds
instance Mutable s CSUSeconds
instance Mutable s CFloat
instance Mutable s CDouble

instance Mutable s a => Mutable s (Identity a) where
    type Ref s (Identity a) = CoerceRef s (Identity a) a

instance Mutable s a => Mutable s (Const a b) where
    type Ref s (Const a b) = CoerceRef s (Const a b) a

instance Mutable s a => Mutable s (V.Const a b) where
    type Ref s (V.Const a b) = CoerceRef s (V.Const a b) a

instance Mutable s a => Mutable s (M.Product a) where
    type Ref s (M.Product a) = CoerceRef s (M.Product a) a

instance Mutable s a => Mutable s (M.Sum a) where
    type Ref s (M.Sum a) = CoerceRef s (M.Sum a) a

instance Mutable s a => Mutable s (Down a) where
    type Ref s (Down a) = CoerceRef s (Down a) a

instance Mutable s a => Mutable s (M.Dual a) where
    type Ref s (M.Dual a) = CoerceRef s (M.Dual a) a

instance Mutable s a => Mutable s (Maybe a) where
    type Ref s (Maybe a) = GRef s (Maybe a)

instance (Mutable s a, Mutable s b) => Mutable s (Either a b) where
    type Ref s (Either a b) = GRef s (Either a b)

instance (Mutable s (f a), Mutable s (g a)) => Mutable s (Product f g a) where
    type Ref s (Product f g a) = GRef s (Product f g a)

instance (Mutable s (f a), Mutable s (g a)) => Mutable s (Sum f g a) where
    type Ref s (Sum f g a) = GRef s (Sum f g a)

instance (Mutable s (f (g a))) => Mutable s (Compose f g a) where
    type Ref s (Compose f g a) = CoerceRef s (Compose f g a) (f (g a))

-- | Mutable linked list with mutable references in each cell.  See
-- 'Data.Mutable.MutBranch' documentation for an example of using this as
-- a mutable linked list.l
instance Mutable s a => Mutable s [a] where
    type Ref s [a] = GRef s [a]

-- | Meant for usage with higher-kinded data pattern (See 'X.HKD')
instance Mutable s a => Mutable s (V.Identity a) where
    type Ref s (V.Identity a) = RefFor s a
    thawRef (V.Identity x) = RefFor <$> thawRef x
    freezeRef (RefFor r) = V.Identity <$> freezeRef r
    copyRef (RefFor r) (V.Identity x) = copyRef r x
    moveRef (RefFor r) (RefFor v) = moveRef r v
    cloneRef = fmap RefFor . cloneRef . getRefFor
    unsafeThawRef (V.Identity x) = RefFor <$> unsafeThawRef x
    unsafeFreezeRef (RefFor r) = V.Identity <$> unsafeFreezeRef r

-- | Mutable reference is 'MV.MVector'.
instance Mutable s (V.Vector a) where
    type Ref s (V.Vector a) = MV.MVector s a
    thawRef         = VG.thaw
    freezeRef       = VG.freeze
    copyRef         = VG.copy
    moveRef         = MVG.move
    cloneRef        = MVG.clone
    unsafeThawRef   = VG.unsafeThaw
    unsafeFreezeRef = VG.unsafeFreeze

-- | Mutable reference is 'MVS.MVector'.
instance Storable a => Mutable s (VS.Vector a) where
    type Ref s (VS.Vector a) = MVS.MVector s a
    thawRef         = VG.thaw
    freezeRef       = VG.freeze
    copyRef         = VG.copy
    moveRef         = MVG.move
    cloneRef        = MVG.clone
    unsafeThawRef   = VG.unsafeThaw
    unsafeFreezeRef = VG.unsafeFreeze

-- | Mutable reference is 'MVU.MVector'.
instance VU.Unbox a => Mutable s (VU.Vector a) where
    type Ref s (VU.Vector a) = MVU.MVector s a
    thawRef         = VG.thaw
    freezeRef       = VG.freeze
    copyRef         = VG.copy
    moveRef         = MVG.move
    cloneRef        = MVG.clone
    unsafeThawRef   = VG.unsafeThaw
    unsafeFreezeRef = VG.unsafeFreeze

-- | Mutable reference is 'MVP.MVector'.
instance Prim a => Mutable s (VP.Vector a) where
    type Ref s (VP.Vector a) = MVP.MVector s a
    thawRef         = VG.thaw
    freezeRef       = VG.freeze
    copyRef         = VG.copy
    moveRef         = MVG.move
    cloneRef        = MVG.clone
    unsafeThawRef   = VG.unsafeThaw
    unsafeFreezeRef = VG.unsafeFreeze

instance Mutable s (Array a) where
    type Ref s (Array a) = MutableArray s a

    thawRef xs = thawArray xs 0 (sizeofArray xs)
    freezeRef rs = freezeArray rs 0 (sizeofMutableArray rs)
    copyRef rs xs = copyArray rs 0 xs 0 l
      where
        l = sizeofArray xs `min` sizeofMutableArray rs
    moveRef rs vs = copyMutableArray rs 0 vs 0 l
      where
        l = sizeofMutableArray vs `min` sizeofMutableArray rs
    cloneRef rs = cloneMutableArray rs 0 (sizeofMutableArray rs)
    unsafeThawRef   = unsafeThawArray
    unsafeFreezeRef = unsafeFreezeArray

instance Mutable s (SmallArray a) where
    type Ref s (SmallArray a) = SmallMutableArray s a

    thawRef xs = thawSmallArray xs 0 (sizeofSmallArray xs)
    freezeRef rs = freezeSmallArray rs 0 (sizeofSmallMutableArray rs)
    copyRef rs xs = copySmallArray rs 0 xs 0 l
      where
        l = sizeofSmallArray xs `min` sizeofSmallMutableArray rs
    moveRef rs vs = copySmallMutableArray rs 0 vs 0 l
      where
        l = sizeofSmallMutableArray vs `min` sizeofSmallMutableArray rs
    cloneRef rs = cloneSmallMutableArray rs 0 (sizeofSmallMutableArray rs)
    unsafeThawRef   = unsafeThawSmallArray
    unsafeFreezeRef = unsafeFreezeSmallArray

instance Mutable s ByteArray where
    type Ref s ByteArray = MutableByteArray s

    thawRef xs = do
        rs <- newByteArray (sizeofByteArray xs)
        copyByteArray rs 0 xs 0 (sizeofByteArray xs)
        pure rs
    freezeRef rs = do
        xs <- newByteArray (sizeofMutableByteArray rs)
        copyMutableByteArray xs 0 rs 0 (sizeofMutableByteArray rs)
        unsafeFreezeByteArray xs
    copyRef rs xs = copyByteArray rs 0 xs 0 l
      where
        l = sizeofByteArray xs `min` sizeofMutableByteArray rs
    moveRef rs vs = copyMutableByteArray rs 0 vs 0 l
      where
        l = sizeofMutableByteArray vs `min` sizeofMutableByteArray rs
    cloneRef rs = do
        vs <- newByteArray (sizeofMutableByteArray rs)
        copyMutableByteArray vs 0 rs 0 (sizeofMutableByteArray rs)
        pure vs
    unsafeThawRef   = unsafeThawByteArray
    unsafeFreezeRef = unsafeFreezeByteArray

instance Prim a => Mutable s (PrimArray a) where
    type Ref s (PrimArray a) = MutablePrimArray s a

    thawRef xs = do
        rs <- newPrimArray (sizeofPrimArray xs)
        copyPrimArray rs 0 xs 0 (sizeofPrimArray xs)
        pure rs
    freezeRef rs = do
        xs <- newPrimArray (sizeofMutablePrimArray rs)
        copyMutablePrimArray xs 0 rs 0 (sizeofMutablePrimArray rs)
        unsafeFreezePrimArray xs
    copyRef rs xs = copyPrimArray rs 0 xs 0 l
      where
        l = sizeofPrimArray xs `min` sizeofMutablePrimArray rs
    moveRef rs vs = copyMutablePrimArray rs 0 vs 0 l
      where
        l = sizeofMutablePrimArray vs `min` sizeofMutablePrimArray rs
    cloneRef rs = do
        vs <- newPrimArray (sizeofMutablePrimArray rs)
        copyMutablePrimArray vs 0 rs 0 (sizeofMutablePrimArray rs)
        pure vs
    unsafeThawRef   = unsafeThawPrimArray
    unsafeFreezeRef = unsafeFreezePrimArray

data VoidRef s
  deriving (Show, Read, Eq, Ord, Functor, Traversable, Foldable)

instance Mutable s Void where
    type Ref s Void = VoidRef s
    thawRef         = \case {}
    freezeRef       = \case {}
    copyRef         = \case {}
    moveRef         = \case {}
    cloneRef        = \case {}
    unsafeThawRef   = \case {}
    unsafeFreezeRef = \case {}

data UnitRef s = UnitRef
  deriving (Show, Read, Eq, Ord, Functor, Traversable, Foldable)

instance Applicative UnitRef where
    pure _  = UnitRef
    _ <*> _ = UnitRef

instance Monad UnitRef where
    return   = pure
    _ >>= _ = UnitRef

instance Mutable s () where
    type Ref s () = UnitRef s
    thawRef   _       = pure UnitRef
    freezeRef _       = pure ()
    copyRef _ _       = pure ()
    moveRef _ _       = pure ()
    cloneRef _        = pure UnitRef
    unsafeThawRef _   = pure UnitRef
    unsafeFreezeRef _ = pure ()

-- | A 'Ref' of a tuple is a tuple of 'Ref's, for easy accessing.
--
-- @
-- Ref s (Int, 'V.Vector' Double) = ('Data.Primitive.MutVar.MutVar' s Int, 'MV.MVector' s Double)
-- @
instance (Mutable s a, Mutable s b) => Mutable s (a, b) where
    type Ref s (a, b) = (Ref s a, Ref s b)
    thawRef   (!x, !y) = (,) <$> thawRef x   <*> thawRef y
    freezeRef (u , v ) = (,) <$> freezeRef u <*> freezeRef v
    copyRef   (u , v ) (!x, !y) = copyRef u x *> copyRef v y
    moveRef   (u , v ) ( x,  y) = moveRef u x *> moveRef v y
    cloneRef  (u , v ) = (,) <$> cloneRef u   <*> cloneRef v
    unsafeThawRef   (!x, !y) = (,) <$> unsafeThawRef x   <*> unsafeThawRef y
    unsafeFreezeRef (u , v ) = (,) <$> unsafeFreezeRef u <*> unsafeFreezeRef v

mutableTuples [3..12]

-- | 'Ref' for components in a vinyl 'Rec'.
newtype RecRef s f a = RecRef { getRecRef :: Ref s (f a) }

deriving instance Eq (Ref s (f a)) => Eq (RecRef s f a)
deriving instance Ord (Ref s (f a)) => Ord (RecRef s f a)

instance Mutable s (Rec f '[]) where
    type Ref s (Rec f '[]) = Rec (RecRef s f) '[]
    thawRef   _       = pure RNil
    freezeRef _       = pure RNil
    copyRef _ _       = pure ()
    moveRef _ _       = pure ()
    cloneRef _        = pure RNil
    unsafeThawRef _   = pure RNil
    unsafeFreezeRef _ = pure RNil

instance ( Mutable s (f a)
         , Mutable s (Rec f as)
         , Ref s (Rec f as) ~ Rec (RecRef s f) as
         ) => Mutable s (Rec f (a ': as)) where
    type Ref s (Rec f (a ': as)) = Rec (RecRef s f) (a ': as)
    thawRef   = \case
      x :& xs -> (:&) <$> (RecRef <$> thawRef x) <*> thawRef xs
    freezeRef = \case
      RecRef v :& vs -> (:&) <$> freezeRef v <*> freezeRef vs
    copyRef = \case
      RecRef v :& vs -> \case
        x :& xs -> copyRef v x >> copyRef vs xs
    moveRef = \case
      RecRef v :& vs -> \case
        RecRef r :& rs ->
          moveRef v r >> moveRef vs rs
    cloneRef = \case
      RecRef v :& rs -> (:&) <$> (RecRef <$> cloneRef v) <*> cloneRef rs
    unsafeThawRef   = \case
      x :& xs -> (:&) <$> (RecRef <$> unsafeThawRef x) <*> unsafeThawRef xs
    unsafeFreezeRef = \case
      RecRef v :& vs -> (:&) <$> unsafeFreezeRef v <*> unsafeFreezeRef vs

instance ( RecApplicative as
         , V.NatToInt (V.RLength as)
         , RPureConstrained (V.IndexableField as) as
         , Mutable s (Rec f as)
         , Ref s (Rec f as) ~ Rec (RecRef s f) as
         ) => Mutable s (ARec f as) where
    type Ref s (ARec f as) = ARec (RecRef s f) as

    thawRef         = fmap toARec . thawRef   . fromARec
    freezeRef       = fmap toARec . freezeRef . fromARec
    copyRef r x     = copyRef (fromARec r) (fromARec x)
    moveRef r v     = moveRef (fromARec r) (fromARec v)
    cloneRef        = fmap toARec . cloneRef . fromARec
    unsafeThawRef   = fmap toARec . unsafeThawRef   . fromARec
    unsafeFreezeRef = fmap toARec . unsafeFreezeRef . fromARec

-- | The mutable reference of the 'HList' type from generic-lens.
data HListRef :: Type -> [Type] -> Type where
    NilRef :: HListRef s '[]
    (:!>)  :: Ref s a -> HListRef s as -> HListRef s (a ': as)
infixr 5 :!>

instance Mutable s (HList '[]) where
    type Ref s (HList '[]) = HListRef s '[]
    thawRef   _       = pure NilRef
    freezeRef _       = pure Nil
    copyRef _ _       = pure ()
    moveRef _ _       = pure ()
    cloneRef _        = pure NilRef
    unsafeThawRef _   = pure NilRef
    unsafeFreezeRef _ = pure Nil

instance (Mutable s a, Mutable s (HList as), Ref s (HList as) ~ HListRef s as) => Mutable s (HList (a ': as)) where
    type Ref s (HList (a ': as)) = HListRef s (a ': as)
    thawRef   = \case
      x :> xs -> (:!>) <$> thawRef x <*> thawRef xs
    freezeRef = \case
      v :!> vs -> (:>) <$> freezeRef v <*> freezeRef vs
    copyRef = \case
      v :!> vs -> \case
        x :> xs -> copyRef v x >> copyRef vs xs
    moveRef = \case
      v :!> vs -> \case
        r :!> rs ->
          moveRef v r >> moveRef vs rs
    cloneRef = \case
      v :!> rs -> (:!>) <$> cloneRef v <*> cloneRef rs
    unsafeThawRef   = \case
      x :> xs -> (:!>) <$> unsafeThawRef x <*> unsafeThawRef xs
    unsafeFreezeRef = \case
      v :!> vs -> (:>) <$> unsafeFreezeRef v <*> unsafeFreezeRef vs

-- ListRefTuple instances

-- this one instance is the reason why we have to define a new typeclass
-- instead of using 'ListRef' -- because of @'ListRef' s () '[]@.
instance ListRefTuple s (UnitRef s) '[] where
    tupleToListRef _ = Nil
    listRefToTuple _ = UnitRef

instance (Ref s a ~ ra) => ListRefTuple s ra '[a] where
    tupleToListRef x = x :> Nil
    listRefToTuple (x :> _) = x

instance (Ref s a ~ ra, Ref s b ~ rb) => ListRefTuple s (ra, rb) '[a, b] where
    tupleToListRef (x, y) = x :> y :> Nil
    listRefToTuple (x :> y :> _) = (x, y)

listRefTuples [3..12]

max_stars_repo_name: alex-rozenshteyn-leapyear/mutable

max_issues_repo_name: alex-rozenshteyn-leapyear/mutable

max_forks_repo_name: alex-rozenshteyn-leapyear/mutable


text:
{-# OPTIONS_GHC -fno-warn-unused-do-bind #-}

module Wyas.Parser.Primitives where

import Wyas.Types

import Text.ParserCombinators.Parsec hiding (spaces)
import Numeric (readOct, readHex)
import Data.Ratio ((%))
import Data.Complex (Complex (..))
import Control.Monad (liftM)

parseBool :: Parser LispVal
parseBool = do char '#'
               b <- parseTrue <|> parseFalse
               return $ Bool b
               where parseTrue  = char 't' >> return True
                     parseFalse = char 'f' >> return False

parseAtom :: Parser LispVal
parseAtom = do first <- letter <|> symbol
               rest <- many (letter <|> digit <|> symbol)
               return $ Atom (first:rest)

--
-- Char & String
--

parseChar :: Parser LispVal
parseChar = liftM Character (string "#\\" >> parseChar')
          where parseWS = (string " " <|> string "space") >> return ' '
                parseNewline = string "newline" >> return '\n'
                parseChar' = parseWS <|> parseNewline <|> anyChar

parseString :: Parser LispVal
parseString = do char '"'
                 x <- many (escapedChars <|> noneOf "\\\"")
                 char '"'
                 (return . String) x

--
-- Parsing Integers
--

parseNumber :: Parser LispVal
parseNumber = parsePlainNumber <|> parseRadixNumber

parsePlainNumber :: Parser LispVal
parsePlainNumber = liftM (Number . read) (many1 digit)

--
-- Radix Numbers
--

parseRadixNumber :: Parser LispVal
parseRadixNumber = char '#' >>
                   (
                        parseDecimal
                        <|> parseBinary
                        <|> parseOct
                        <|> parseHex
                   )

parseDecimal :: Parser LispVal
parseDecimal = do char 'd'
                  n <- many1 digit
                  (return . Number . read) n

parseBinary :: Parser LispVal
parseBinary = do char 'b'
                 n <- many $ oneOf "01"
                 (return . Number . bin2int) n

parseHex :: Parser LispVal
parseHex = do char 'x'
              n <- many $ oneOf "0123456789abcdefABCDEF"
              return . Number . readWith readHex $ n
              where readWith f s = fst $ head (f s)

parseOct :: Parser LispVal
parseOct = do char 'o'
              n <- many $ oneOf "01234567"
              return . Number . readWith readOct $ n
              where readWith f s = fst $ head (f s)

--
-- Numeric Tower
--

parseFloat :: Parser LispVal
parseFloat = do x <- many1 digit
                char '.'
                y <- many1 digit
                return $ Float (read (x ++ "." ++y))

parseRational :: Parser LispVal
parseRational = do n <- many1 digit
                   char '/'
                   d <- many1 digit
                   return $ Ratio $ read n % read d

parseComplex :: Parser LispVal
parseComplex = do x <- try parseFloat <|> parsePlainNumber
                  char '+'
                  y <- try parseFloat <|> parsePlainNumber
                  char 'i'
                  return $ Complex (toDouble x :+ toDouble y)

--
-- Helpers
--

escapedChars :: Parser Char
escapedChars = char '\\' >>
               (
                    char '\\' <|> char '"' <|>
                    parseNewline <|>
                    parseReturn <|>
                    parseTab
               )
               where parseNewline = char 'n' >> return '\n'
                     parseReturn  = char 'r' >> return '\r'
                     parseTab     = char 't' >> return '\t'

symbol :: Parser Char
symbol = oneOf "!$%|*+-/:<=>?@^_~"

spaces :: Parser ()
spaces = skipMany1 space

bin2int :: String -> Integer
bin2int = bin2int' 0
        where bin2int' digint "" = digint
              bin2int' digint (x:xs) = let old = 2 * digint + (if x == '0' then 0 else 1)
                                       in bin2int' old xs

-- TODO clean up
toDouble :: LispVal -> Double
toDouble (Float f) = f
toDouble (Number n) = fromIntegral n
toDouble _ = error "ERROR: not a double"

max_stars_repo_name: grtlr/wyas

max_issues_repo_name: grtlr/wyas

max_forks_repo_name: grtlr/wyas


text:
{-# LANGUAGE DeriveGeneric     #-}
{-# LANGUAGE FlexibleContexts  #-}
{-# LANGUAGE ViewPatterns      #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE InstanceSigs      #-}

module Main where

import Numeric.LinearAlgebra.Data
import Numeric.LinearAlgebra.Devel          ( MatrixOrder(..)
                                            , matrixFromVector)
import Numeric.LinearAlgebra.HMatrix
-- import Control.Lens
import Control.Applicative
import qualified Data.ByteString.Lazy as BL
import Data.Csv
import Data.List
import qualified Data.Vector          as V
import qualified Data.Vector.Storable as VS
import GHC.Generics (Generic)
-- import Unsafe.Coerce

import TinyMath             as Tiny

matrixFromVector' :: VS.Storable a => MatrixOrder -> Int -> Int -> V.Vector [a] -> Matrix a
matrixFromVector' o r c = matrixFromVector o r c . VS.fromList . concat . V.toList

identity :: Int -> Matrix Double
identity = ident

pause :: IO ()
pause = do
  putStrLn "\n Program paused. Press enter to continue."
  _ <- getChar
  return ()

class Decodable a where
  decode' :: HasHeader -> BL.ByteString -> V.Vector a

instance Decodable [Double] where
  decode' :: HasHeader -> BL.ByteString -> V.Vector [Double]
  decode' h f = d
    where
      d = case decode NoHeader f :: Either String (V.Vector [Double]) of
        Left err -> V.empty
        Right xs -> xs

instance Decodable (Double, Double) where
  decode' :: HasHeader -> BL.ByteString -> V.Vector (Double, Double)
  decode' h f = d
    where
      d = case decode NoHeader f :: Either String (V.Vector (Double, Double)) of
        Left err -> V.empty
        Right xs -> xs

-- instance VS.Vector Integer where
--   konst :: Integer -> (Int, Int) -> Matrix Integer
--   konst = undefined

--computeCost :: Matrix a -> Matrix a -> a
computeCost :: Matrix Double -> Matrix Double -> Double
computeCost xs theta = cost
  where
    m    = fromIntegral . fst . size $ xs   :: Double
    x    = xs ¿ [0,1]
    y    = xs ¿ [2]
    hyp  = mul x theta                      :: Matrix Double
    temp = (hyp - y)^2                      :: Matrix Double
    cost = ((2*m)**(-1) * sumElements temp ) :: Double

-- app xs (VS.fromList [1.0,1.0,1.0])
-- ((2*m)**(-1)) * (sumElements $ (xs :: Matrix Double) ¿ [1])

gradientDescent :: Matrix a -> Matrix a -> Double -> Integer -> Double
gradientDescent = undefined

-- ===================== Main                   ====================

main :: IO ()
main = do
--  ==================== Part 1: Basic Function ====================
  putStrLn "Running warmUpExercise ..."
  putStrLn "5x5 Identity Matrix:"
  print (identity 5)
  pause

-- ======================= Part 2: Plotting =======================
  putStrLn "Plotting Data..."
  f <- BL.readFile "./data/ex1data1.txt"
  let
    d = decode' NoHeader f
    in Tiny.plot $ V.toList d

  pause

-- =================== Part 3: Cost and Gradient descent ===================
  let
    -- | Some gradient descent settings
    iterations = 1500 :: Integer
    alpha      = 0.1  :: Double

    d     = decode' NoHeader f :: V.Vector [Double]
    mtx   = matrixFromVector' ColumnMajor 2 (V.length d) d
    m     = V.length d
    -- | a Matrix with a column of 1's preadded.
    -- | initialize fittig parameters
    xs    = konst 1 (m, 1) ||| tr mtx      -- | Add a column of ones to x
    theta = konst 0 (2,1)      :: Matrix Double -- | initialize fitting parameters

    -- | compute and display initial cost
    cost  = computeCost xs theta
  print cost

  let
    theta'= gradientDescent xs theta alpha iterations

    -- m =
  return ()


text:
{-# LANGUAGE ForeignFunctionInterface #-}

module Numeric.LinearAlgebra.HArpack
       ( eigs
       , eigsN
       , eigsS
       , linearOperator
       , SparseMatrix (..)
       , Problem (..)
       , SWhich (..)
       , NWhich (..)
       , Result (..)
       , Dimension
       , HowMany
       , Tolerance
       , Iterations ) where

import           Control.Monad
import           Control.Monad.Loops
import           Data.Complex
import           Data.Vector.Storable         (Vector)
import qualified Data.Vector.Storable         as V
import           Data.Vector.Storable.Mutable (IOVector)
import qualified Data.Vector.Storable.Mutable as MV
import           Foreign
import           Foreign.C.String
import           Foreign.C.Types
import           Foreign.ForeignPtr.Unsafe    as FPUnsafe
import           System.IO.Unsafe

-- | Functions of type 'LinearOperator' represent a matrix vector product.
-- The values in the second vector should be overwritten by operator applied to
-- the first.

data SparseMatrix =
  SparseMatrix { dim ::Int
               , indexes :: [((Int,Int),Double)] }
  deriving Show

data LinearOperator =
  LinearOperator
    { applyOperator :: IOVector CDouble -> IOVector CDouble -> IO () }

linearOperator :: [((Int,Int),Double)] -> LinearOperator
linearOperator l = LinearOperator $ \x y ->
  do MV.set y 0
     forM_ l $ \((i,j),v) -> do
       v' <- (* realToFrac v) <$> MV.read x j
       MV.modify y (+ v') i

--------------------------------------------------------------------------------
-- | Which eigenvalues to calculate:
-- - The first character refers to 'S'ymmetric or 'N'on-Symmetric;
-- - The second refers to 'L'argest or 'S'mallest;
-- - The third refers to 'M'agnitude, 'A'lgebraic, 'R'eal, 'I'mmaginary;
-- - The option 'BE' means 'BE'tween, for central eigenvalues.

data SWhich = SLM | SSM | SLA | SSA | SBE

instance Show SWhich where
  show SLM = "LM"
  show SSM = "SM"
  show SLA = "LA"
  show SSA = "SA"
  show SBE = "BE"

----------

data NWhich = NLM | NSM | NLR | NSR | NLI | NSI

instance Show NWhich where
  show NLM = "LM"
  show NSM = "SM"
  show NLR = "LR"
  show NSR = "SR"
  show NLI = "LI"
  show NSI = "SI"

----------

data Problem = Symmetric SWhich | NonSymmetric NWhich
     deriving (Show)

data Eigenpairs = Error String
                | EReal [(Double,Vector Double)]
                | EComplex [(Complex Double, Vector (Complex Double))]

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

-- | Encapsulates the arrays and data needed to invoke _naupd routines in ARPACK

data ArpackSetup =
  ArpackSetup { aIDO    :: Ptr CInt
              , aBMAT   :: Ptr CChar
              , aN      :: Ptr CInt
              , aWHICH  :: Ptr CChar
              , aNEV    :: Ptr CInt
              , aTOL    :: Ptr CDouble
              , aRESID  :: Ptr CDouble
              , aNCV    :: Ptr CInt
              , aV      :: Ptr CDouble
              , aLDV    :: Ptr CInt
              , aIPARAM :: Ptr CInt
              , aIPNTR  :: Ptr CInt
              , aWORKD  :: Ptr CDouble
              , aWORKL  :: Ptr CDouble
              , aLWORKL :: Ptr CInt
              , aINFO   :: Ptr CInt }
  deriving (Show)

type Dimension = Int
type HowMany = Int
type Tolerance = Double
type Iterations = Int

-- | The 'setupArpack' function allocates the arrays and data
-- needed by ARPACK to solve an eigenvalue problem

setupArpack :: Problem
            -> Dimension
            -> HowMany
            -> Tolerance
            -> Iterations
            -> IO ArpackSetup

setupArpack (NonSymmetric which') n nev tol maxIter =
  do let ncv = min n (2 * nev + 1)
         lWorkl = 3 * ncv * ncv + 6 * ncv
         ldv = n

     idoPtr <- ptrInt 0
     nPtr   <- ptrInt n
     nevPtr <- ptrInt nev
     ncvPtr <- ptrInt ncv
     ldvPtr <- ptrInt ldv

     bmatPtr  <- ptrCharArray "I"
     whichPtr <- ptrCharArray (show which')

     tolPtr   <- ptrDouble tol

     residPtr <- mallocDoubleArray n
     vPtr     <- mallocDoubleArray (n * ncv)
     workdPtr <- mallocDoubleArray (3 * n)
     worklPtr <- mallocDoubleArray lWorkl

     iParamPtr <- mallocIntArray 11

     pokeElemOff iParamPtr 0 1
     pokeElemOff iParamPtr 2 (fromIntegral maxIter)
     pokeElemOff iParamPtr 3 1
     pokeElemOff iParamPtr 6 1

     iPntrPtr  <- mallocIntArray 14

     lworklPtr <- ptrInt lWorkl
     infoPtr   <- ptrInt 0

     setup <- return $
              ArpackSetup idoPtr bmatPtr nPtr whichPtr nevPtr tolPtr
                          residPtr ncvPtr vPtr ldvPtr iParamPtr iPntrPtr
                          workdPtr worklPtr lworklPtr infoPtr

     dnaupd setup
     return setup

setupArpack (Symmetric which) n nev tol maxIter =
  do let ncv = min n (2 * nev + 1)
         lWorkl = ncv * ncv + 8 * ncv
         ldv = n

     idoPtr <- ptrInt 0
     nPtr   <- ptrInt n
     nevPtr <- ptrInt nev
     ncvPtr <- ptrInt ncv
     ldvPtr <- ptrInt ldv

     bmatPtr  <- ptrCharArray "I"
     whichPtr <- ptrCharArray (show which)

     tolPtr   <- ptrDouble tol

     residPtr <- mallocDoubleArray n
     vPtr     <- mallocDoubleArray (n * ncv)
     workdPtr <- mallocDoubleArray (3 * n)
     worklPtr <- mallocDoubleArray lWorkl

     iParamPtr <- mallocIntArray 11

     pokeElemOff iParamPtr 0 1
     pokeElemOff iParamPtr 2 (fromIntegral maxIter)
     pokeElemOff iParamPtr 3 1
     pokeElemOff iParamPtr 6 1

     iPntrPtr  <- mallocIntArray 11

     lworklPtr <- ptrInt lWorkl
     infoPtr   <- ptrInt 0

     setup <- return $
              ArpackSetup idoPtr bmatPtr nPtr whichPtr nevPtr tolPtr
                          residPtr ncvPtr vPtr ldvPtr iParamPtr iPntrPtr
                          workdPtr worklPtr lworklPtr infoPtr

     dsaupd setup

     return setup

-- |Invoke ARPACK's iteration through the reverse communication interface.
iterateArpack :: Problem -> LinearOperator -> ArpackSetup -> IO ()
iterateArpack (NonSymmetric _) f setup = do
  workdXElem <- peekElemOff (aIPNTR setup) 0
  workdYElem <- peekElemOff (aIPNTR setup) 1

  n <- peek (aN setup)

  xPtr <- newForeignPtr_ $ advancePtr (aWORKD setup) (fromIntegral workdXElem - 1)
  yPtr <- newForeignPtr_ $ advancePtr (aWORKD setup) (fromIntegral workdYElem - 1)

  let x = MV.unsafeFromForeignPtr xPtr 0 (fromIntegral n)
      y = MV.unsafeFromForeignPtr yPtr 0 (fromIntegral n)

  applyOperator f x y

  dnaupd setup

iterateArpack (Symmetric _) f setup = do
  workdXElem <- peekElemOff (aIPNTR setup) 0
  workdYElem <- peekElemOff (aIPNTR setup) 1

  n <- peek (aN setup)

  xPtr <- newForeignPtr_ $ advancePtr (aWORKD setup) (fromIntegral workdXElem - 1)
  yPtr <- newForeignPtr_ $ advancePtr (aWORKD setup) (fromIntegral workdYElem - 1)

  let x = MV.unsafeFromForeignPtr xPtr 0 (fromIntegral n)
      y = MV.unsafeFromForeignPtr yPtr 0 (fromIntegral n)

  applyOperator f x y

  dsaupd setup

arpack :: LinearOperator -> Problem -> Dimension -> HowMany -> Tolerance ->
          Iterations -> IO ArpackSetup
arpack op problem n nev tol mxItr = do
  ar <- setupArpack problem n nev tol mxItr
  iterateArpack problem op ar
  whileM_ ((== 1) <$> peek (aIDO ar)) (iterateArpack problem op ar)
  return ar

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

data D = DReal (IOVector CDouble)
       | DComplex (IOVector CDouble) (IOVector CDouble)

data Sigma = SReal (Ptr CDouble)
           | SComplex (Ptr CDouble) (Ptr CDouble)

data ArpackResults =
  ArpackResults { eRVEC    :: Ptr CInt
                , eHOWMNY  :: Ptr CChar
                , eSELECT  :: Ptr CInt
                , eD       :: D
                , eZ       :: IOVector CDouble
                , eLDZ     :: Ptr CInt
                , eSIGMA   :: Sigma
                , eWORKEV  :: Maybe (Ptr CDouble)
                , eSetup   :: ArpackSetup }

errCheck :: ArpackResults -> IO Bool
errCheck ar = do
  info <- (peek . aINFO . eSetup) ar
  print info
  return (info < 0)

parseArpackOutput :: Problem -> ArpackSetup -> IO ArpackResults
parseArpackOutput (NonSymmetric _) setup = do
  ncv <- peek (aNCV setup)
  nev <- peek (aNEV setup)
  n <- peek (aN setup)

  rvecPtr <- ptrInt 1
  howmnyPtr <- ptrCharArray "A"
  drVec <- newVectorDouble (nev + 1)
  diVec <- newVectorDouble (nev + 1)
  selectPtr <- mallocIntArray (fromIntegral ncv)
  zVec <- newVectorDouble (n*(nev+1))
  ldzPtr <- ptrInt (fromIntegral n)
  sigmarPtr <- mallocDouble
  sigmaiPtr <- mallocDouble
  workevPtr <- mallocDoubleArray (fromIntegral (3*ncv))

  let (drfPtr,_,_) = MV.unsafeToForeignPtr drVec
      drPtr = FPUnsafe.unsafeForeignPtrToPtr drfPtr
      (difPtr,_,_) = MV.unsafeToForeignPtr diVec
      diPtr = FPUnsafe.unsafeForeignPtrToPtr difPtr
      (zfPtr,_,_) = MV.unsafeToForeignPtr zVec
      zPtr = FPUnsafe.unsafeForeignPtrToPtr zfPtr

  c_dneupd rvecPtr howmnyPtr selectPtr drPtr diPtr zPtr ldzPtr sigmarPtr
    sigmaiPtr workevPtr (aBMAT setup) (aN setup) (aWHICH setup) (aNEV setup)
    (aTOL setup) (aRESID setup) (aNCV setup) (aV setup) (aLDV setup)
    (aIPARAM setup) (aIPNTR setup) (aWORKD setup) (aWORKL setup) (aLWORKL setup)
    (aINFO setup)

  return $ ArpackResults rvecPtr howmnyPtr selectPtr (DComplex drVec diVec) zVec
               ldzPtr (SComplex sigmarPtr sigmaiPtr) (Just workevPtr) setup

parseArpackOutput (Symmetric _) setup = do
  ncv <- peek (aNCV setup)
  nev <- peek (aNEV setup)
  n <- peek (aN setup)

  rvecPtr <- ptrInt 1
  howmnyPtr <- ptrCharArray "A"
  dVec <- newVectorDouble nev
  selectPtr <- mallocIntArray (fromIntegral ncv)
  zVec <- newVectorDouble (n*(nev+1))
  ldzPtr <- ptrInt (fromIntegral n)
  sigmaPtr <- mallocDouble

  let (drfPtr,_,_) = MV.unsafeToForeignPtr dVec
      dPtr = FPUnsafe.unsafeForeignPtrToPtr drfPtr
      (zfPtr,_,_) = MV.unsafeToForeignPtr zVec
      zPtr = FPUnsafe.unsafeForeignPtrToPtr zfPtr

  c_dseupd rvecPtr howmnyPtr selectPtr dPtr zPtr ldzPtr
    sigmaPtr (aBMAT setup) (aN setup) (aWHICH setup) (aNEV setup)
    (aTOL setup) (aRESID setup) (aNCV setup) (aV setup) (aLDV setup)
    (aIPARAM setup) (aIPNTR setup) (aWORKD setup) (aWORKL setup) (aLWORKL setup)
    (aINFO setup)

  return $ ArpackResults rvecPtr howmnyPtr selectPtr (DReal dVec) zVec
               ldzPtr (SReal sigmaPtr) Nothing setup

getComplexEigenvalue :: ArpackResults -> Int -> IO (Complex Double)
getComplexEigenvalue results index = do
  let DComplex reEigs imEigs = eD results
  re <- MV.read reEigs index
  im <- MV.read imEigs index
  return (realToFrac re :+ realToFrac im)

getRealEigenvector :: ArpackResults -> Int -> IO (Vector (Complex Double))
getRealEigenvector results index = do
  n <- (peek . aN . eSetup) results
  ldz <- peek $ eLDZ results
  let z = eZ results
  evec <- V.freeze $ MV.slice (index * fromIntegral ldz) (fromIntegral n) z
  return $ V.map ((:+0) . realToFrac) evec

getComplexEigenvector
  :: ArpackResults -> Int -> Bool -> IO (Vector (Complex Double))
getComplexEigenvector results index conjugate = do
  n <- (peek . aN . eSetup) results
  ldz <- peek $ eLDZ results
  let z = eZ results
  reEvec <- V.freeze $ MV.slice (index * fromIntegral ldz) (fromIntegral n) z
  imEvec <- V.freeze $ MV.slice ((index+1) * fromIntegral ldz) (fromIntegral n) z
  if conjugate
     then return $
          V.zipWith (:+) (V.map realToFrac reEvec)
                         (V.map (negate . realToFrac) imEvec)
     else return $
          V.zipWith (:+) (V.map realToFrac reEvec)
                         (V.map realToFrac imEvec)

getEigenpair
  :: ArpackResults -> Int -> IO (Complex Double, V.Vector (Complex Double))
getEigenpair results index = do
  eig <- getComplexEigenvalue results index
  let evec = case compare (imagPart eig) 0 of
        EQ -> getRealEigenvector results index
        LT -> getComplexEigenvector results (index - 1) True
        GT -> getComplexEigenvector results index False
  evec' <- evec
  return (eig,evec')

-- | The function 'eigs' computes a few eigenvalues of a real valued
-- unsymmetric matrix by invoking the ARPACK library. First element is true if
-- ARPACK converged. If converged the second element is a list of eigenpairs.

eigsN :: LinearOperator
      -> NWhich
      -> Dimension
      -> HowMany
      -> Tolerance
      -> Iterations
      -> Maybe [(Complex Double, Vector (Complex Double))]
eigsN f which' n' nev' tol' iters = unsafePerformIO $ do
  setup <- arpack f (NonSymmetric which') n' nev' tol' iters
  results <- parseArpackOutput (NonSymmetric which') setup
  err <- errCheck results
  if err
    then return Nothing
    else do
      pairs <- mapM (getEigenpair results) [0..(nev'-1)]
      return (Just pairs)

data Result = RError
            | RReal [(Double,Vector Double)]
            | RComplex [(Complex Double, Vector (Complex Double))]
              deriving Show

eigs :: SparseMatrix
     -> Problem
     -> HowMany
     -> Result
eigs m (NonSymmetric w) h =
  case eigsN (linearOperator (indexes m)) w (dim m) h (1e-10) 6000 of
  Nothing -> RError
  Just x -> RComplex x
eigs m (Symmetric w) h =
  case eigsS (linearOperator (indexes m)) w (dim m) h (1e-10) 6000 of
  Nothing -> RError
  Just x -> RReal x

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

getSymmetricEigenvalue :: ArpackResults -> Int -> IO Double
getSymmetricEigenvalue results index = do
  let DReal eigV = eD results
  eig <- MV.read eigV index
  return (realToFrac eig)

getSymmetricEigenvector :: ArpackResults -> Int -> IO (Vector Double)
getSymmetricEigenvector results index = do
  n <- fromIntegral <$> (peek . aN . eSetup) results
  ldz <- peek $ eLDZ results
  let z = eZ results
  evec <- V.freeze $ MV.slice (index * fromIntegral ldz) n z
  return $ V.map realToFrac evec

getSymmetricEigenpair :: ArpackResults -> Int -> IO (Double, Vector Double)
getSymmetricEigenpair results index = do
  eig <- getSymmetricEigenvalue results index
  evec <- getSymmetricEigenvector results index
  return (eig,evec)

eigsS :: LinearOperator
      -> SWhich
      -> Dimension
      -> HowMany
      -> Tolerance
      -> Iterations
      -> Maybe [(Double,Vector Double)]
eigsS f which n nev tol' iters = unsafePerformIO $ do
  setup <- arpack f (Symmetric which) n nev tol' iters
  results <- parseArpackOutput (Symmetric which) setup
  err <- errCheck results
  if err
    then return Nothing
    else do
      pairs <- mapM (getSymmetricEigenpair results) [0..(nev-1)]
      return (Just pairs)

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

mallocInt :: IO (Ptr CInt)
mallocInt = malloc

ptrInt :: Int -> IO (Ptr CInt)
ptrInt k =
  do ptr <- malloc
     poke ptr (fromIntegral k)
     return ptr

mallocIntArray :: Int -> IO (Ptr CInt)
mallocIntArray = mallocArray

----------

mallocDouble :: IO (Ptr CDouble)
mallocDouble = malloc

ptrDouble :: Double -> IO (Ptr CDouble)
ptrDouble x =
  do ptr <- malloc
     poke ptr (realToFrac x)
     return ptr

mallocDoubleArray :: Int -> IO (Ptr CDouble)
mallocDoubleArray = mallocArray

newVectorDouble :: CInt -> IO (MV.IOVector CDouble)
newVectorDouble = MV.new . fromIntegral

----------

mallocChar :: IO (Ptr CChar)
mallocChar = malloc

mallocCharArray :: Int -> IO (Ptr CChar)
mallocCharArray = mallocArray

ptrCharArray :: String -> IO (Ptr CChar)
ptrCharArray s =
  do ptr <- mallocCharArray (length s + 1)
     pokeArray ptr (map castCharToCChar s)
     return ptr

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

foreign import ccall unsafe "arpack.h dnaupd_"
  c_dnaupd :: Ptr CInt    -- ^ reverse communication flag, to be started with 0
           -> Ptr CChar   -- ^ 'I' for standard problem, 'G' for generalized
           -> Ptr CInt    -- ^ dimension of the eigenproblem
           -> Ptr CChar   -- ^ which eigenpairs to compute
           -> Ptr CInt    -- ^ N : number of eigenpairs to compute
           -> Ptr CDouble -- ^ tolerance
           -> Ptr CDouble -- ^ residual vector
           -> Ptr CInt    -- ^ NCV : number of columns of V matrix
           -> Ptr CDouble -- ^ N x NCV output array
           -> Ptr CInt    -- ^ leading dimension of V
           -> Ptr CInt    -- ^ array for internal parameters, length 11
           -> Ptr CInt    -- ^ pointer of starting locations, length 14
           -> Ptr CDouble -- ^ workd : work array
           -> Ptr CDouble -- ^ workl : work array
           -> Ptr CInt    -- ^ length of workl
           -> Ptr CInt    -- ^ info
           -> IO ()

dnaupd :: ArpackSetup -> IO ()
dnaupd setup =
  c_dnaupd (aIDO setup) (aBMAT setup) (aN setup) (aWHICH setup) (aNEV setup)
           (aTOL setup) (aRESID setup) (aNCV setup) (aV setup) (aLDV setup)
           (aIPARAM setup) (aIPNTR setup) (aWORKD setup) (aWORKL setup)
           (aLWORKL setup) (aINFO setup)

foreign import ccall unsafe "arpack.h dneupd_"
  c_dneupd :: Ptr CInt    -- ^ whether to calculate eigenvectors
           -> Ptr CChar   -- ^ whether to calculate Ritz or Schur vectors
           -> Ptr CInt    -- ^ selects the vectors to be computed
           -> Ptr CDouble -- ^ real part of eigenvalues
           -> Ptr CDouble -- ^ imaginary part of eigenvalues
           -> Ptr CDouble -- ^ approximate eigenvectors
           -> Ptr CInt    -- ^ leading dimension of former array
           -> Ptr CDouble -- ^ real part of the shift
           -> Ptr CDouble -- ^ imaginary part of the shift
           -> Ptr CDouble -- ^ work array
           -> Ptr CChar -> Ptr CInt -> Ptr CChar -> Ptr CInt -> Ptr CDouble
           -> Ptr CDouble -> Ptr CInt -> Ptr CDouble -> Ptr CInt -> Ptr CInt
           -> Ptr CInt -> Ptr CDouble -> Ptr CDouble -> Ptr CInt -> Ptr CInt
           -> IO ()

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

foreign import ccall unsafe "arpack.h dsaupd_"
  c_dsaupd :: Ptr CInt    -- ^ reverse communication flag, to be started with 0
           -> Ptr CChar   -- ^ 'I' for standard problem, 'G' for generalized
           -> Ptr CInt    -- ^ dimension of the eigenproblem
           -> Ptr CChar   -- ^ which eigenpairs to compute
           -> Ptr CInt    -- ^ N : number of eigenpairs to compute
           -> Ptr CDouble -- ^ tolerance
           -> Ptr CDouble -- ^ residual vector
           -> Ptr CInt    -- ^ NCV : number of columns of V matrix
           -> Ptr CDouble -- ^ N x NCV output array
           -> Ptr CInt    -- ^ leading dimension of V
           -> Ptr CInt    -- ^ array for internal parameters, length 11
           -> Ptr CInt    -- ^ pointer of starting locations, length 14
           -> Ptr CDouble -- ^ workd : work array
           -> Ptr CDouble -- ^ workl : work array
           -> Ptr CInt    -- ^ length of workl
           -> Ptr CInt    -- ^ info
           -> IO ()

dsaupd :: ArpackSetup -> IO ()
dsaupd setup =
  c_dsaupd (aIDO setup) (aBMAT setup) (aN setup) (aWHICH setup) (aNEV setup)
           (aTOL setup) (aRESID setup) (aNCV setup) (aV setup) (aLDV setup)
           (aIPARAM setup) (aIPNTR setup) (aWORKD setup) (aWORKL setup)
           (aLWORKL setup) (aINFO setup)

foreign import ccall unsafe "arpack.h dseupd_"
  c_dseupd :: Ptr CInt    -- ^ whether to calculate eigenvectors
           -> Ptr CChar   -- ^ whether to calculate Ritz or Schur vectors
           -> Ptr CInt    -- ^ selects the vectors to be computed
           -> Ptr CDouble -- ^ eigenvalues
           -> Ptr CDouble -- ^ approximate eigenvectors
           -> Ptr CInt    -- ^ leading dimension of former array
           -> Ptr CDouble -- ^ shift
           -> Ptr CChar -> Ptr CInt -> Ptr CChar -> Ptr CInt -> Ptr CDouble
           -> Ptr CDouble -> Ptr CInt -> Ptr CDouble -> Ptr CInt -> Ptr CInt
           -> Ptr CInt -> Ptr CDouble -> Ptr CDouble -> Ptr CInt -> Ptr CInt
           -> IO ()

max_stars_repo_name: guaraqe/harpack

max_issues_repo_name: guaraqe/harpack

max_forks_repo_name: guaraqe/harpack


text:
{-# LANGUAGE Rank2Types #-}
{-# LANGUAGE TypeFamilies #-}

module Plotting where

import Data.List (sort)
import qualified Data.Vector as Vector
import Graphics.Rendering.Chart.Easy
import Statistics.Sample (meanVariance)

{-# ANN module "HLint: ignore Reduce duplication" #-}

errBars :: Double -> Vector.Vector Double -> (Double, Double)
errBars k xs = (mean, k * stdDev)
  where
    (mean, var) = meanVariance xs
    stdDev = sqrt var

anytimeEval :: ([(a, Double)] -> Double) -> [(a, Double)] -> [Int] -> [Double]
anytimeEval f xs ns = map eval (filter (<= length xs) ns)
  where
    eval n = f (take n xs)

anytimePlot :: String -> String -> [Int] -> [(String, [Double])] -> EC (Layout LogValue LogValue) ()
anytimePlot x_title y_title ns inputs = do
  let xLabelShow = map (show . ceiling)
  let yLabelShow = map show
  let generatePlot (algName, ys) = plot $ line algName [zip (map fromIntegral ns) (map LogValue ys)]
  mapM_ generatePlot inputs
  layout_x_axis . laxis_generate .= autoScaledLogAxis (loga_labelf .~ xLabelShow $ def)
  layout_y_axis . laxis_generate .= autoScaledLogAxis (loga_labelf .~ yLabelShow $ def)
  layout_x_axis . laxis_title .= x_title
  layout_y_axis . laxis_title .= y_title

oneShotPlot :: String -> String -> [(String, [(Double, (Double, Double))])] -> EC (Layout LogValue LogValue) ()
oneShotPlot x_title y_title inputs = do
  let xLabelShow = map (show . ceiling)
  let yLabelShow = map show
  let processPoint (x, (y, dy)) = symErrPoint (LogValue x) (LogValue y) 0 (LogValue dy)
  let generatePlot (algName, ps) = do
        plot $ points algName $ map (\(x, (y, _)) -> (LogValue x, LogValue y)) ps
        plot $ return (def & plot_errbars_values .~ map processPoint ps)
  mapM_ generatePlot inputs
  layout_x_axis . laxis_generate .= autoScaledLogAxis (loga_labelf .~ xLabelShow $ def)
  layout_y_axis . laxis_generate .= autoScaledLogAxis (loga_labelf .~ yLabelShow $ def)
  layout_x_axis . laxis_title .= x_title
  layout_y_axis . laxis_title .= y_title

errorbarPlot ::
  String ->
  String ->
  [(String, [(Double, [Double])])] ->
  EC (Layout LogValue LogValue) ()
errorbarPlot x_title y_title inputs = do
  let xLabelShow = map (show . ceiling)
  let yLabelShow = map show
  let makeBars (x, ys) = ErrPoint (ErrValue x' x' x') (ErrValue ymin y ymax)
        where
          x' = LogValue x
          y = LogValue $ sum ys / fromIntegral (length ys)
          sorted = sort ys
          n = length ys
          ymin = LogValue $ sorted !! (n `div` 4)
          ymax = LogValue $ sorted !! (n * 3 `div` 4)
  let generatePlot (algName, ps) = do
        plot $ liftEC $ do
          -- ensure error bars are the same color as the points
          (color : _) <- liftCState $ use colors
          plot_errbars_values .= map makeBars ps
          plot_errbars_line_style . line_color .= color
        plot $ points algName $
          map (\(x, ys) -> (LogValue x, LogValue (sum ys / fromIntegral (length ys)))) ps
  mapM_ generatePlot inputs
  layout_x_axis . laxis_generate .= autoScaledLogAxis (loga_labelf .~ xLabelShow $ def)
  layout_y_axis . laxis_generate .= autoScaledLogAxis (loga_labelf .~ yLabelShow $ def)
  layout_x_axis . laxis_title .= x_title
  layout_y_axis . laxis_title .= y_title

max_stars_repo_name: saeedhadikhanloo/monad-bayes

max_issues_repo_name: saeedhadikhanloo/monad-bayes

max_forks_repo_name: saeedhadikhanloo/monad-bayes


text:
module Numeric.LinearAlgebra.Matrix.Mat44 where

import Numeric.LinearAlgebra.Matrix.Class
import Numeric.LinearAlgebra.Matrix.Mat33 ( det33 )

import Data.List ( foldl1' )
import Data.NumInstances() -- Add function instances to Num

data Mat44 a = Mat44
  { m00 :: !a, m01 :: !a, m02 :: !a, m03 :: !a
  , m10 :: !a, m11 :: !a, m12 :: !a, m13 :: !a
  , m20 :: !a, m21 :: !a, m22 :: !a, m23 :: !a
  , m30 :: !a, m31 :: !a, m32 :: !a, m33 :: !a
  }
  deriving (Read, Show, Eq, Ord)

instance Matrix Mat44 where
  {-# INLINE mDim #-}
  mDim _ = 4
  {-# INLINE mElement #-}
  mElement m 0 0 = m00 m
  mElement m 0 1 = m01 m
  mElement m 0 2 = m02 m
  mElement m 0 3 = m03 m
  mElement m 1 0 = m10 m
  mElement m 1 1 = m11 m
  mElement m 1 2 = m12 m
  mElement m 1 3 = m13 m
  mElement m 2 0 = m20 m
  mElement m 2 1 = m21 m
  mElement m 2 2 = m22 m
  mElement m 2 3 = m23 m
  mElement m 3 0 = m30 m
  mElement m 3 1 = m31 m
  mElement m 3 2 = m32 m
  mElement m 3 3 = m33 m
  mElement _ i j = error ("Index " ++ show i ++ ", " ++ show j ++ ": out of range, must be 0,0 to 3,3")
  {-# INLINE mZip #-}
  mZip f m n = Mat44 (f (m00 m) (m00 n)) (f (m01 m) (m01 n)) (f (m02 m) (m02 n)) (f (m03 m) (m03 n))
                     (f (m10 m) (m10 n)) (f (m11 m) (m11 n)) (f (m12 m) (m12 n)) (f (m13 m) (m13 n))
                     (f (m20 m) (m20 n)) (f (m21 m) (m21 n)) (f (m22 m) (m22 n)) (f (m23 m) (m23 n))
                     (f (m30 m) (m30 n)) (f (m31 m) (m31 n)) (f (m32 m) (m32 n)) (f (m33 m) (m33 n))
  {-# INLINE mFold #-}
  mFold f m = foldl1' f [ mElement m i j | i <- [ 0 .. 3 ], j <- [ 0 .. 3 ] ]
  {-# INLINE mIndexOf #-}
  mIndexOf p m = fst (foldl1' p' [ ((i,j), mElement m i j) | j <- [ 3, 2 .. 0 ], i <- [ 3, 2 .. 0 ] ])
    where
    p' x@(_, a) y@(_, a') | a `p` a'  = x
                          | otherwise = y
  {-# INLINE det #-}
  det = det44

instance Functor Mat44 where
  {-# INLINE fmap #-}
  fmap f m = Mat44 (f (m00 m)) (f (m01 m)) (f (m02 m)) (f (m03 m))
                   (f (m10 m)) (f (m11 m)) (f (m12 m)) (f (m13 m))
                   (f (m20 m)) (f (m21 m)) (f (m22 m)) (f (m23 m))
                   (f (m30 m)) (f (m31 m)) (f (m32 m)) (f (m33 m))

{-# SPECIALIZE det44 :: Mat44 Double -> Double #-}
{-# SPECIALIZE det44 :: Mat44 Float  -> Float #-}
det44 :: Num a => Mat44 a -> a
det44 = m00 * det33 m11 m12 m13
                    m21 m22 m23
                    m31 m32 m33
      - m01 * det33 m10 m12 m13
                    m20 m22 m23
                    m30 m32 m33
      + m02 * det33 m10 m11 m13
                    m20 m21 m23
                    m30 m31 m33
      - m03 * det33 m10 m11 m12
                    m20 m21 m22
                    m30 m31 m32

{-# SPECIALIZE inv44 :: Mat44 Double -> Mat44 Double #-}
{-# SPECIALIZE inv44 :: Mat44 Float  -> Mat44 Float  #-}
inv44 :: Fractional a => Mat44 a -> Mat44 a
inv44 m = mApply m' m
  where
  d  = const $ det44 m
  m' = Mat44 (det33 m11 m12 m13
                    m21 m22 m23
                    m31 m32 m33 / d)
           (-(det33 m01 m02 m03
                    m21 m22 m23
                    m31 m32 m33 / d))
             (det33 m01 m02 m03
                    m11 m12 m13
                    m31 m32 m33 / d)
           (-(det33 m01 m02 m03
                    m11 m12 m13
                    m21 m22 m23 / d))

           (-(det33 m10 m12 m13
                    m20 m22 m23
                    m30 m32 m33 / d))
             (det33 m00 m02 m03
                    m20 m22 m23
                    m30 m32 m33 / d)
           (-(det33 m00 m02 m03
                    m10 m12 m13
                    m30 m32 m33 / d))
             (det33 m00 m02 m03
                    m10 m12 m13
                    m20 m22 m23 / d)

             (det33 m10 m11 m13
                    m20 m21 m23
                    m30 m31 m33 / d)
           (-(det33 m00 m01 m03
                    m20 m21 m23
                    m30 m31 m33 / d))
             (det33 m00 m01 m03
                    m10 m11 m13
                    m30 m31 m33 / d)
           (-(det33 m00 m01 m03
                    m10 m11 m13
                    m20 m21 m23 / d))

           (-(det33 m10 m11 m12
                    m20 m21 m22
                    m30 m31 m32 / d))
             (det33 m00 m01 m02
                    m20 m21 m22
                    m30 m31 m32 / d)
           (-(det33 m00 m01 m02
                    m10 m11 m12
                    m30 m31 m32 / d))
             (det33 m00 m01 m02
                    m10 m11 m12
                    m20 m21 m22 / d)

max_stars_repo_name: dagit/lin-alg

max_issues_repo_name: dagit/lin-alg

max_forks_repo_name: dagit/lin-alg


text:
{-# LANGUAGE DataKinds                                #-}
{-# LANGUAGE DeriveAnyClass                           #-}
{-# LANGUAGE DeriveDataTypeable                       #-}
{-# LANGUAGE DeriveGeneric                            #-}
{-# LANGUAGE DerivingVia                              #-}
{-# LANGUAGE FlexibleInstances                        #-}
{-# LANGUAGE GADTs                                    #-}
{-# LANGUAGE KindSignatures                           #-}
{-# LANGUAGE MultiParamTypeClasses                    #-}
{-# LANGUAGE PatternSynonyms                          #-}
{-# LANGUAGE RankNTypes                               #-}
{-# LANGUAGE RecordWildCards                          #-}
{-# LANGUAGE ScopedTypeVariables                      #-}
{-# LANGUAGE StandaloneDeriving                       #-}
{-# LANGUAGE TemplateHaskell                          #-}
{-# LANGUAGE TypeApplications                         #-}
{-# LANGUAGE TypeFamilies                             #-}
{-# LANGUAGE TypeInType                               #-}
{-# LANGUAGE TypeOperators                            #-}
{-# LANGUAGE UndecidableInstances                     #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}

module Backprop.Learn.Model.Neural.LSTM (
    -- * LSTM
    lstm
  , LSTMp(..), lstmForget, lstmInput, lstmUpdate, lstmOutput
  , reshapeLSTMpInput
  , reshapeLSTMpOutput
  , lstm'
    -- * GRU
  , gru
  , GRUp(..), gruMemory, gruUpdate, gruOutput
  , gru'
  ) where

import           Backprop.Learn.Initialize
import           Backprop.Learn.Model.Function
import           Backprop.Learn.Model.Neural
import           Backprop.Learn.Model.Regression
import           Backprop.Learn.Model.State
import           Backprop.Learn.Model.Types
import           Control.DeepSeq
import           Control.Monad
import           Control.Monad.Primitive
import           Data.Type.Tuple
import           Data.Typeable
import           GHC.Generics                          (Generic)
import           GHC.TypeNats
import           Lens.Micro
import           Lens.Micro.TH
import           Numeric.Backprop
import           Numeric.LinearAlgebra.Static.Backprop
import           Numeric.OneLiner
import           Numeric.Opto.Ref
import           Numeric.Opto.Update
import           Statistics.Distribution
import qualified Data.Binary                           as Bi
import qualified Numeric.LinearAlgebra.Static          as H
import qualified System.Random.MWC                     as MWC

-- TODO: allow parameterize internal activation function?
-- TODO: Peepholes

-- | 'LSTM' layer parmateters
data LSTMp (i :: Nat) (o :: Nat) =
    LSTMp { _lstmForget :: !(FCp (i + o) o)
          , _lstmInput  :: !(FCp (i + o) o)
          , _lstmUpdate :: !(FCp (i + o) o)
          , _lstmOutput :: !(FCp (i + o) o)
          }
  deriving stock     (Generic, Typeable, Show)
  deriving anyclass  (NFData, Linear Double, Metric Double, Bi.Binary, Regularize, Backprop)

deriving via (GNum (LSTMp i o)) instance (KnownNat i, KnownNat o) => Num (LSTMp i o)
deriving via (GNum (LSTMp i o)) instance (KnownNat i, KnownNat o) => Fractional (LSTMp i o)
deriving via (GNum (LSTMp i o)) instance (KnownNat i, KnownNat o) => Floating (LSTMp i o)

makeLenses ''LSTMp

instance (PrimMonad m, KnownNat i, KnownNat o) => Mutable m (LSTMp i o) where
    type Ref m (LSTMp i o) = GRef m (LSTMp i o)
    thawRef   = gThawRef
    freezeRef = gFreezeRef
    copyRef   = gCopyRef
instance (PrimMonad m, KnownNat i, KnownNat o) => LinearInPlace m Double (LSTMp i o)

instance (PrimMonad m, KnownNat i, KnownNat o) => Learnable m (LSTMp i o)

-- | Stateless version of 'lstm' that takes the "previous input" as a part
-- of the input vector.
lstm'
    :: (KnownNat i, KnownNat o)
    => Model ('Just (LSTMp i o)) ('Just (R o)) (R (i + o)) (R o)
lstm' = modelD $ \(PJust p) x (PJust s) ->
    let forget = logistic $ runLRp (p ^^. lstmForget) x
        input  = logistic $ runLRp (p ^^. lstmInput ) x
        update = tanh     $ runLRp (p ^^. lstmUpdate) x
        s'     = forget * s + input * update
        o      = logistic $ runLRp (p ^^. lstmOutput) x
        h      = o * tanh s'
    in  (h, PJust s')

-- | Long-term short-term memory layer
--
-- <http://colah.github.io/posts/2015-08-Understanding-LSTMs/>
--
lstm
    :: (KnownNat i, KnownNat o)
    => Model ('Just (LSTMp i o)) ('Just (R o :# R o)) (R i) (R o)
lstm = recurrent H.split (H.#) id lstm'

reshapeLSTMpInput
    :: (ContGen d, PrimMonad m, KnownNat i, KnownNat i', KnownNat o)
    => d
    -> MWC.Gen (PrimState m)
    -> LSTMp i o
    -> m (LSTMp i' o)
reshapeLSTMpInput d g (LSTMp forget input update output) =
    LSTMp <$> reshaper forget
          <*> reshaper input
          <*> reshaper update
          <*> reshaper output
  where
    reshaper = reshapeLRpInput d g

reshapeLSTMpOutput
    :: (ContGen d, PrimMonad m, KnownNat i, KnownNat o, KnownNat o')
    => d
    -> MWC.Gen (PrimState m)
    -> LSTMp i o
    -> m (LSTMp i o')
reshapeLSTMpOutput d g (LSTMp forget input update output) =
    LSTMp <$> reshaper forget
          <*> reshaper input
          <*> reshaper update
          <*> reshaper output
  where
    reshaper = reshapeLRpInput  d g
           <=< reshapeLRpOutput d g

-- | Forget biases initialized to 1
instance (KnownNat i, KnownNat o) => Initialize (LSTMp i o) where
    initialize d g = LSTMp <$> set (mapped . fcBias) 1 (initialize d g)
                           <*> initialize d g
                           <*> initialize d g
                           <*> initialize d g

-- | 'GRU' layer parmateters
data GRUp (i :: Nat) (o :: Nat) =
    GRUp { _gruMemory :: !(FCp (i + o) o)
         , _gruUpdate :: !(FCp (i + o) o)
         , _gruOutput :: !(FCp (i + o) o)
         }
  deriving stock     (Generic, Typeable, Show)
  deriving anyclass  (NFData, Linear Double, Metric Double, Bi.Binary, Initialize, Regularize, Backprop)

deriving via (GNum (GRUp i o)) instance (KnownNat i, KnownNat o) => Num (GRUp i o)
deriving via (GNum (GRUp i o)) instance (KnownNat i, KnownNat o) => Fractional (GRUp i o)
deriving via (GNum (GRUp i o)) instance (KnownNat i, KnownNat o) => Floating (GRUp i o)

makeLenses ''GRUp

instance (PrimMonad m, KnownNat i, KnownNat o) => Mutable m (GRUp i o) where
    type Ref m (GRUp i o) = GRef m (GRUp i o)
    thawRef   = gThawRef
    freezeRef = gFreezeRef
    copyRef   = gCopyRef
instance (KnownNat i, KnownNat o, Mutable m (GRUp i o)) => LinearInPlace m Double (GRUp i o)

instance (KnownNat i, KnownNat o, PrimMonad m) => Learnable m (GRUp i o)

-- | Stateless version of 'gru' that takes the "previous input" as a part
-- of the input vector.
gru'
    :: forall i o. (KnownNat i, KnownNat o)
    => Model ('Just (GRUp i o)) 'Nothing (R (i + o)) (R o)
gru' = modelStatelessD $ \(PJust p) x ->
    let z      = logistic $ runLRp (p ^^. gruMemory) x
        r      = logistic $ runLRp (p ^^. gruUpdate) x
        r'     = 1 # r
        h'     = tanh     $ runLRp (p ^^. gruOutput) (r' * x)
    in  (1 - z) * snd (split @i x) + z * h'

-- | Gated Recurrent Unit
--
-- <http://colah.github.io/posts/2015-08-Understanding-LSTMs/>
--
gru :: (KnownNat i, KnownNat o)
    => Model ('Just (GRUp i o)) ('Just (R o)) (R i) (R o)
gru = recurrent H.split (H.#) id gru'

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
module Gates where

import           Protolude

import           Control.Exception (assert)
import           Data.Complex      (Complex (..), magnitude)
import           Data.List         ((!!))
import qualified Data.Matrix       as M
import qualified Data.Vector       as V
import           System.Random     (randomIO)

import           Instruction       (Operation, Qubit, Qubits)
import           State             (QuantumState)

-- Dimension in qubits by taking the log base 2
dimensionQubits :: Integral a => a -> Int
dimensionQubits = ceiling . logBase 2 . fromIntegral

i :: Operation
i = M.identity 2

x :: Operation
x = M.fromLists [[0, 1], [1, 0]]

h :: Operation
h = M.scaleMatrix (1 / sqrt 2) $ M.fromLists [[1, 1], [1, -1]]

cnot :: Operation
cnot = M.fromLists [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]]

swap :: Operation
swap = M.fromLists [[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]]

-- Apply an operation to a quantum state by matrix multiplication
applyOperation :: Operation -> QuantumState -> QuantumState
applyOperation g s = M.getCol 1 $ g * M.colVector s

-- Compose operators by doing matrix multiplication
-- NOTE: (*) operator chosen over multStd for possible speed improvements for larger matrices
composeOperations :: Operation -> Operation -> Operation
composeOperations = (*)

-- Sample quantum state, randomly selecting a state, weighted by |phi_i|^2
sample :: MonadIO m => QuantumState -> m Int
sample s = do
  r <- liftIO $ randomIO
  return $ sampleIter r 0
  where
    sampleIter :: Double -> Int -> Int
    sampleIter r t =
      let newR = r - (magnitude (s V.! t)) ** 2
       in if newR < 0
            then t
            else sampleIter newR (t + 1)

-- Collapse the quantum state into the state at the given index
collapse :: QuantumState -> Int -> QuantumState
collapse s i =
  V.generate
    (V.length s)
    (\x ->
       if x == i
         then (1 :+ 0)
         else (0 :+ 0))

-- Observe the quantum state and collapse into that state
observe :: MonadIO m => QuantumState -> m QuantumState
observe s = do
  i <- sample s
  return $ collapse s i

-- Oh, to have dependent types
kroneckerMultiply :: (Num a) => M.Matrix a -> M.Matrix a -> M.Matrix a
kroneckerMultiply a b = M.flatten $ M.mapPos (\_ -> flip M.scaleMatrix b) a

-- TODO: re-eval constraints
kroneckerExpt :: (Num a) => M.Matrix a -> Int -> M.Matrix a
kroneckerExpt a n
  | n < 1 = M.fromLists [[1]]
  | n == 1 = a
  | otherwise = kroneckerMultiply (kroneckerExpt a (n - 1)) a

-- Lift an operation by padding with identity
liftOperation :: Operation -> Qubit -> Int -> Operation
liftOperation a index n =
  assert (index < n) kroneckerMultiply left (kroneckerMultiply a right)
  where
    left = kroneckerExpt i (n - index - (dimensionQubits $ M.nrows a))
    right = kroneckerExpt i index

-- Calculate needed transpositions from given permutation
permutationToTranspositions :: [Int] -> [(Int, Int)]
permutationToTranspositions p = permutationToTranspositionsIter 0
  where
    permutationToTranspositionsIter dest
      | dest == (length p) = []
      | otherwise =
        (case compare src dest of
           LT -> [(src, dest)]
           GT -> [(dest, src)]
           _  -> []) ++
        permutationToTranspositionsIter (dest + 1)
      where
        src = nextSrc (p !! dest) dest
    nextSrc s d
      | s < d = nextSrc (p !! s) d
      | otherwise = s

-- Convert transpositions to list of adjacent transpositions
transpositionToAdjacentTranspositions :: [(Int, Int)] -> [Int]
transpositionToAdjacentTranspositions ts = foldMap expandCons ts
  where
    expandCons (x, y)
      | y - x == 1 = [x]
      | otherwise =
        let trans = [x .. y - 1]
         in trans ++ (reverse (initSafe trans))

-- Apply 1 qubit gate
apply1QGate :: QuantumState -> Operation -> Qubit -> QuantumState
apply1QGate s a q =
  applyOperation (liftOperation a q (dimensionQubits (length s))) s

-- Apply multiple qubit gate
applynQGate :: QuantumState -> Operation -> Qubits -> QuantumState
applynQGate s a qubits = do
  applyOperation upq s
  where
    n :: Int
    n = dimensionQubits (length s)
    transpositionsToOperation :: [Int] -> Operation
    transpositionsToOperation ts =
      foldr composeOperations (M.identity (2 ^ n)) $
      map (\x -> liftOperation Gates.swap x n) ts
    u01 :: Operation
    u01 = liftOperation a 0 n
    fromSpace :: [Int]
    fromSpace =
      (reverse qubits) ++ (filter (\x -> not (elem x qubits)) [0 .. n - 1])
    trans :: [Int]
    trans =
      transpositionToAdjacentTranspositions . permutationToTranspositions $
      fromSpace
    toFrom :: Operation
    toFrom = transpositionsToOperation trans
    fromTo :: Operation
    fromTo = transpositionsToOperation (reverse trans)
    upq :: Operation
    upq = composeOperations toFrom (composeOperations u01 fromTo)

-- Apply gate, picking apply1QGate or applynQGate based on size
applyGate :: QuantumState -> Operation -> [Int] -> QuantumState
applyGate s a qubits =
  assert
    (length qubits == dimensionQubits (M.nrows a))
    (if length qubits == 1
       then apply1QGate s a (qubits !! 0)
       else applynQGate s a qubits)

max_stars_repo_name: colescott/quantum-interpreter

max_issues_repo_name: colescott/quantum-interpreter

max_forks_repo_name: colescott/quantum-interpreter


text:
{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE CPP                 #-}
{-# LANGUAGE TemplateHaskell     #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE KindSignatures      #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# OPTIONS_GHC -fno-warn-missing-signatures #-}

module Test.Grenade.Layers.PadCrop where

import           Grenade

import           Hedgehog

import           Numeric.LinearAlgebra.Static ( norm_Inf )

import           Test.Hedgehog.Hmatrix

prop_pad_crop :: Property
prop_pad_crop =
  let net :: Network '[Pad 2 3 4 6, Crop 2 3 4 6] '[ 'D3 7 9 5, 'D3 16 15 5, 'D3 7 9 5 ]
      net = Pad :~> Crop :~> NNil
  in  property $
    forAll genOfShape >>= \(d :: S ('D3 7 9 5)) ->
      let (tapes, res)  = runForwards  net d
          (_    , grad) = runBackwards net tapes d
      in  do assert $ d ~~~ res
             assert $ grad ~~~ d

prop_pad_crop_2d :: Property
prop_pad_crop_2d =
  let net :: Network '[Pad 2 3 4 6, Crop 2 3 4 6] '[ 'D2 7 9, 'D2 16 15, 'D2 7 9 ]
      net = Pad :~> Crop :~> NNil
  in  property $
    forAll genOfShape >>= \(d :: S ('D2 7 9)) ->
      let (tapes, res)  = runForwards  net d
          (_    , grad) = runBackwards net tapes d
      in  do assert $ d ~~~ res
             assert $ grad ~~~ d

(~~~) :: S x -> S x -> Bool
(S1D x) ~~~ (S1D y) = norm_Inf (x - y) < 0.00001
(S2D x) ~~~ (S2D y) = norm_Inf (x - y) < 0.00001
(S3D x) ~~~ (S3D y) = norm_Inf (x - y) < 0.00001

tests :: IO Bool
tests = checkParallel $$(discover)

max_issues_repo_name: Alien-Inc/grenade

max_forks_repo_name: Alien-Inc/grenade