text:
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE DataKinds #-}

-- | Internal representation of Feldspar programs

module Feldspar.Representation where

import Control.Monad.Reader
import Data.Complex
import Data.Hash (hashInt)
import Data.Int
import Data.List (genericTake)
import Data.Proxy
import Data.Typeable (Typeable)
import Data.Word

import Data.Constraint (Dict (..))

import Language.Syntactic
import Language.Syntactic.Functional
import Language.Syntactic.Functional.Tuple
import Language.Syntactic.TH

import qualified Control.Monad.Operational.Higher as Operational

import qualified Language.Embedded.Expression as Imp
import qualified Language.Embedded.Imperative.CMD as Imp
import qualified Language.Embedded.Backend.C.Expression as Imp

import Data.Inhabited
import Data.Selection
import Data.TypedStruct
import Feldspar.Primitive.Representation
import Feldspar.Primitive.Backend.C ()

--------------------------------------------------------------------------------
-- * Object-language types
--------------------------------------------------------------------------------

-- | Representation of all supported types
type TypeRep = Struct PrimType' PrimTypeRep

-- | Supported types
class (Eq a, Show a, Typeable a, Inhabited a) => Type a
  where
    -- | Reify a type
    typeRep :: TypeRep a

instance Type Bool             where typeRep = Single BoolT
instance Type Int8             where typeRep = Single Int8T
instance Type Int16            where typeRep = Single Int16T
instance Type Int32            where typeRep = Single Int32T
instance Type Int64            where typeRep = Single Int64T
instance Type Word8            where typeRep = Single Word8T
instance Type Word16           where typeRep = Single Word16T
instance Type Word32           where typeRep = Single Word32T
instance Type Word64           where typeRep = Single Word64T
instance Type Float            where typeRep = Single FloatT
instance Type Double           where typeRep = Single DoubleT
instance Type (Complex Float)  where typeRep = Single ComplexFloatT
instance Type (Complex Double) where typeRep = Single ComplexDoubleT
instance (Type a, Type b) => Type (a,b) where typeRep = Two typeRep typeRep

-- | Alias for the conjunction of 'PrimType'' and 'Type'
class    (PrimType' a, Type a) => PrimType a
instance (PrimType' a, Type a) => PrimType a

instance Imp.CompTypeClass PrimType
  where
    compType _ = Imp.compType (Proxy :: Proxy PrimType')
    compLit _  = Imp.compLit (Proxy :: Proxy PrimType')
  -- This instance is used by <https://github.com/kmate/raw-feldspar-mcs>

-- | Convert any 'Struct' with a 'PrimType' constraint to a 'TypeRep'
toTypeRep :: Struct PrimType' c a -> TypeRep a
toTypeRep = mapStruct (const primTypeRep)

-- | Check whether two type representations are equal
typeEq :: TypeRep a -> TypeRep b -> Maybe (Dict (a ~ b))
typeEq (Single t) (Single u) = primTypeEq t u
typeEq (Two t1 t2) (Two u1 u2) = do
    Dict <- typeEq t1 u1
    Dict <- typeEq t2 u2
    return Dict
typeEq _ _ = Nothing

-- | Reflect a 'TypeRep' to a 'Typeable' constraint
witTypeable :: TypeRep a -> Dict (Typeable a)
witTypeable (Single t)
    | Dict <- witPrimType t
    = Dict
witTypeable (Two ta tb)
    | Dict <- witTypeable ta
    , Dict <- witTypeable tb
    = Dict

-- | Representation of supported value types + N-ary functions over such types
data TypeRepFun a
  where
    ValT :: TypeRep a -> TypeRepFun a
    FunT :: TypeRep a -> TypeRepFun b -> TypeRepFun (a -> b)
  -- Another option would have been to make `FunT` a constructor in `TypeRep`.
  -- That would have got rid of the extra layer at the expense of less accurate
  -- types (functions would be allowed in pairs, etc.). The current design was
  -- chosen in order to be able to reuse `Struct` instead of making `TypeRep` a
  -- new data type.

-- | Check whether two type representations are equal
typeEqFun :: TypeRepFun a -> TypeRepFun b -> Maybe (Dict (a ~ b))
typeEqFun (ValT t)     (ValT u)     = typeEq t u
typeEqFun (FunT ta tb) (FunT ua ub) = do
    Dict <- typeEq ta ua
    Dict <- typeEqFun tb ub
    return Dict
typeEqFun _ _ = Nothing

-- | Mutable variable
newtype Ref a = Ref { unRef :: Struct PrimType' Imp.Ref (Internal a) }
  -- A reference to a tuple is a struct of smaller references. This means that
  -- creating a reference to a tuple will generate several calls to generate new
  -- references. This must be done already in the front end, which means that
  -- the work in the back end becomes simpler.
  --
  -- Another option would be to allow a single reference to refer to a tuple,
  -- and then turn that into smaller references in the back end. However, this
  -- would complicate the back end for no obvious reason. (One way to do it
  -- would be to run the back end in a store that maps each front end reference
  -- to a struct of small references. Among other things, that would require
  -- dynamic typing.)

-- | Reference specialized to 'Data' elements
type DRef a = Ref (Data a)

instance Eq (Ref a)
  where
    Ref a == Ref b = and $ zipListStruct (==) a b

-- | Mutable array
data Arr a = Arr
    { arrOffset :: Data Index
    , arrLength :: Data Length
    , unArr     :: Struct PrimType' (Imp.Arr Index) (Internal a)
    }
  -- `arrOffset` gives the offset into the internal arrays. The user should not
  -- be able to access the offset or the internal arrays.
  --
  -- `arrLength` gives the maximum index+1 according to the view provided by
  -- `arrIx`. Thus, the maximum index in the internal arrays is
  -- `arrLength + arrOffset - 1`.
  -- `arrLength` should not be exported directly to prevent the user from using
  -- record syntax to change the length of an array. But an equivalent function
  -- can be exported.

  -- An array of tuples is represented as a struct of smaller arrays. See
  -- comment to `Ref`.

-- | Array specialized to 'Data' elements
type DArr a = Arr (Data a)

-- | '==' checks if two 'Arr' use the same physical array. The length and offset
-- are ignored.
instance Eq (Arr a)
  where
    Arr _ _ arr1 == Arr _ _ arr2 = and (zipListStruct (==) arr1 arr2)

-- | Immutable array
data IArr a = IArr
    { iarrOffset :: Data Index
    , iarrLength :: Data Length
    , unIArr     :: Struct PrimType' (Imp.IArr Index) (Internal a)
    }

-- | Immutable array specialized to 'Data' elements
type DIArr a = IArr (Data a)

-- | Check if an 'Arr' and and 'IArr' use the same physical array. The length
-- and offset are ignored. This operation may give false negatives, but never
-- false positives. Whether or not false negatives occur may also depend on the
-- interpretation of the program.
--
-- Due to this unreliability, the function should only be used to affect the
-- non-functional properties of a program (e.g. to avoid unnecessary array
-- copying).
unsafeEqArrIArr :: Arr a -> IArr a -> Bool
unsafeEqArrIArr (Arr _ _ arr1) (IArr _ _ arr2) =
    and (zipListStruct sameId arr1 arr2)
  where
    sameId :: Imp.Arr Index a -> Imp.IArr Index a -> Bool
    sameId (Imp.ArrComp a) (Imp.IArrComp i) = a==i
    sameId _ _ = False

--------------------------------------------------------------------------------
-- * Pure expressions
--------------------------------------------------------------------------------

-- | Assertion labels
data AssertionLabel
    = InternalAssertion
        -- ^ Internal assertion to guarantee meaningful results
    | LibraryAssertion String
        -- ^ Assertion related to a specific library
    | UserAssertion String
        -- ^ Assertion in user code. The default label for user assertions is
        --   @`UserAssertion` ""@
  deriving (Eq, Show)

-- | A selection that includes all labels defined as 'UserAssertion'
onlyUserAssertions :: Selection AssertionLabel
onlyUserAssertions = selectBy $ \l -> case l of
    UserAssertion _ -> True
    _ -> False

data ExtraPrimitive sig
  where
    -- Integer division assuming `divBalanced x y * y == x` (i.e. no remainder).
    -- The purpose of the assumption is to use it in simplifications.
    DivBalanced :: (Integral a, PrimType' a) =>
        ExtraPrimitive (a :-> a :-> Full a)

    -- Guard a value by an assertion
    GuardVal ::
        AssertionLabel -> String -> ExtraPrimitive (Bool :-> a :-> Full a)

instance Eval ExtraPrimitive
  where
    evalSym DivBalanced = \a b -> if a `mod` b /= 0
      then error $ unwords ["divBalanced", show a, show b, "is not balanced"]
      else div a b
    evalSym (GuardVal _ msg) = \cond a ->
        if cond then a else error $ "Feldspar assertion failure: " ++ msg

instance Equality ExtraPrimitive
  where
    equal DivBalanced    DivBalanced    = True
    equal (GuardVal _ _) (GuardVal _ _) = True
    equal _ _ = False

    hash DivBalanced    = hashInt 1
    hash (GuardVal _ _) = hashInt 2

-- | For loop
data ForLoop sig
  where
    ForLoop :: Type st =>
        ForLoop (Length :-> st :-> (Index -> st -> st) :-> Full st)

instance Eval ForLoop
  where
    evalSym ForLoop = \len init body ->
        foldl (flip body) init $ genericTake len [0..]

-- | Interaction with the IO layer
data Unsafe sig
  where
    -- Turn a program into a pure value
    UnsafePerform :: Comp (Data a) -> Unsafe (Full a)

instance Render Unsafe
  where
    renderSym (UnsafePerform _) = "UnsafePerform ..."

instance Eval Unsafe
  where
    evalSym s = error $ "eval: cannot evaluate unsafe operation " ++ renderSym s

-- | 'equal' always returns 'False'
instance Equality Unsafe
  where
    equal _ _ = False

type FeldConstructs
    =   BindingT
    :+: Let
    :+: Tuple
    :+: Primitive
    :+: ExtraPrimitive
    :+: ForLoop
    :+: Unsafe

type FeldDomain = FeldConstructs :&: TypeRepFun

newtype Data a = Data { unData :: ASTF FeldDomain a }

-- | Declaring 'Data' as syntactic sugar
instance Syntactic (Data a)
  where
    type Domain (Data a)   = FeldDomain
    type Internal (Data a) = a
    desugar = unData
    sugar   = Data

instance Syntactic (Struct PrimType' Data a)
    -- Note that this instance places no constraints on `a`. This is crucial in
    -- the way it is used in the rest of the code. It would be possible to
    -- define `desugar` and `sugar` in terms of the instance for pairs; however,
    -- that would require constraining `a`.
  where
    type Domain   (Struct PrimType' Data a) = FeldDomain
    type Internal (Struct PrimType' Data a) = a

    desugar (Single a) = unData a
    desugar (Two a b)  = sugarSymDecor (ValT $ Two ta tb) Pair a' b'
      where
        a' = desugar a
        b' = desugar b
        ValT ta = getDecor a'
        ValT tb = getDecor b'

    sugar a = case getDecor a of
        ValT (Single _)  -> Single $ Data a
        ValT (Two ta tb) ->
            Two (sugarSymDecor (ValT ta) Fst a) (sugarSymDecor (ValT tb) Snd a)

-- | Specialization of the 'Syntactic' class for the Feldspar domain
class    (Syntactic a, Domain a ~ FeldDomain, Type (Internal a)) => Syntax a
instance (Syntactic a, Domain a ~ FeldDomain, Type (Internal a)) => Syntax a

-- | Make a smart constructor for a symbol
sugarSymFeld
    :: ( Signature sig
       , fi         ~ SmartFun FeldDomain sig
       , sig        ~ SmartSig fi
       , FeldDomain ~ SmartSym fi
       , SyntacticN f fi
       , sub :<: FeldConstructs
       , Type (DenResult sig)
       )
    => sub sig -> f
sugarSymFeld = sugarSymDecor $ ValT typeRep

-- | Make a smart constructor for a symbol
sugarSymFeldPrim
    :: ( Signature sig
       , fi         ~ SmartFun FeldDomain sig
       , sig        ~ SmartSig fi
       , FeldDomain ~ SmartSym fi
       , SyntacticN f fi
       , sub :<: FeldConstructs
       , PrimType' (DenResult sig)
       )
    => sub sig -> f
sugarSymFeldPrim = sugarSymDecor $ ValT $ Single primTypeRep

instance Imp.FreeExp Data
  where
    type FreePred Data = PrimType'
    constExp = sugarSymFeldPrim . Lit
    varExp   = sugarSymFeldPrim . FreeVar

--------------------------------------------------------------------------------
-- * Monadic computations
--------------------------------------------------------------------------------

data AssertCMD fs a
  where
    Assert
        :: AssertionLabel
        -> exp Bool
        -> String
        -> AssertCMD (Operational.Param3 prog exp pred) ()

instance Operational.HFunctor AssertCMD
  where
    hfmap _ (Assert c cond msg) = Assert c cond msg

instance Operational.HBifunctor AssertCMD
  where
    hbimap _ g (Assert c cond msg) = Assert c (g cond) msg

instance Operational.InterpBi AssertCMD IO (Operational.Param1 PrimType')
  where
    interpBi (Assert _ cond msg) = do
        cond' <- cond
        unless cond' $ error $ "Assertion failed: " ++ msg

type CompCMD
  =               Imp.RefCMD
  Operational.:+: Imp.ArrCMD
  Operational.:+: Imp.ControlCMD
  Operational.:+: AssertCMD

-- | Monad for computational effects: mutable data structures and control flow
newtype Comp a = Comp
    { unComp ::
        Operational.Program CompCMD (Operational.Param2 Data PrimType') a
    }
  deriving (Functor, Applicative, Monad)

--------------------------------------------------------------------------------
-- Template Haskell instances
--------------------------------------------------------------------------------

deriveSymbol    ''ExtraPrimitive
deriveRender id ''ExtraPrimitive

deriveSymbol    ''ForLoop
deriveRender id ''ForLoop
deriveEquality  ''ForLoop

deriveSymbol ''Unsafe

--------------------------------------------------------------------------------
-- * Interpretation
--------------------------------------------------------------------------------

-- The stuff below depends on the TH-generated instances. From GHC 9, it seems
-- that generated instances are only in scope after the splice. The alternative
-- to move the splices up into the code is not possible (at least not easily)
-- because of the mutual dependencies between the different types. (Definitions
-- that occur after a splice are not in scope before the splice.)

instance EvalEnv ExtraPrimitive env
instance StringTree ExtraPrimitive

instance EvalEnv ForLoop env
instance StringTree ForLoop

instance EvalEnv Unsafe env
instance StringTree Unsafe

-- | Evaluate a closed expression
eval :: (Syntactic a, Domain a ~ FeldDomain) => a -> Internal a
eval = evalClosed . desugar
  -- Note that a `Syntax` constraint would rule out evaluating functions

instance Imp.EvalExp Data
  where
    evalExp = eval

max_stars_repo_name: Feldspar/raw-feldspar

max_issues_repo_name: Feldspar/raw-feldspar

max_forks_repo_name: Feldspar/raw-feldspar


text:
{-# OPTIONS_GHC -Wall #-}
{-# LANGUAGE Trustworthy #-}
{-# LANGUAGE CPP #-}

{- |
Module      :  Physics.Learn.BeamStack
Copyright   :  (c) Scott N. Walck 2016-2018
License     :  BSD3 (see LICENSE)
Maintainer  :  Scott N. Walck <walck@lvc.edu>
Stability   :  experimental

Splitters, recombiners, and detectors for Stern-Gerlach
experiments.
-}

-- Spin-1/2 mixed states.

module Physics.Learn.BeamStack
    (
    -- * Core laboratory components
      BeamStack()
    , randomBeam
    , split
    , recombine
    , applyBField
    , dropBeam
    , flipBeams
    , numBeams
    , detect
    -- * Standard splitters
    , splitX
    , splitY
    , splitZ
    -- * Standard magnetic fields
    , applyBFieldX
    , applyBFieldY
    , applyBFieldZ
    -- * Standard combiners
    , recombineX
    , recombineY
    , recombineZ
    -- * Filters
    , xpFilter
    , xmFilter
    , zpFilter
    , zmFilter
    )
    where

import Physics.Learn.QuantumMat
    ( zp
    , zm
    , nm
    , np
    , couter
    , oneQubitMixed
    )
import Numeric.LinearAlgebra
    ( C
    , Vector
    , Matrix
    , iC
    , (<>)
    , kronecker
    , fromLists
    , toList
    , toLists
    , scale
    , size
    , takeDiag
    , ident
    , tr
    )
import Data.Complex
    ( Complex(..)
    , realPart
    , imagPart
    )
import Data.List
    ( intercalate
    )
#if MIN_VERSION_base(4,11,0)
import Prelude hiding ((<>))
#endif

data BeamStack = BeamStack (Matrix C)

showOneBeam :: Double -> String
showOneBeam r = "Beam of intensity " ++ show r

instance Show BeamStack where
    show b = intercalate "\n" $ map showOneBeam (detect b)

{-
unBeamStack :: BeamStack -> Matrix C
unBeamStack (BeamStack m) = m
-}

--------------------
-- Core functions --
--------------------

-- | A beam of randomly oriented spin-1/2 particles.
randomBeam :: BeamStack
randomBeam = BeamStack oneQubitMixed

extendWithZeros :: Matrix C -> Matrix C
extendWithZeros m
    = let (_,q) = size m
          ml = toLists m
      in fromLists $ map (++ [0,0]) ml
             ++ [replicate (q+2) 0, replicate (q+2) 0]

-- reduce row and column size by 2
reduceMat :: Matrix C -> Matrix C
reduceMat m
    = let (p,q) = size m
          ml = toLists m
      in fromLists $ take (p-2) $ map (take (q-2)) ml

checkedRealPart :: C -> Double
checkedRealPart c
    = let eps = 1e-14
      in if imagPart c < eps
         then realPart c
         else error $ "checkRealPart: imagPart = " ++ show (imagPart c)

-- | Return the intensities of a stack of beams.
detect :: BeamStack -> [Double]
detect (BeamStack m)
    = addAlternate $ toList $ takeDiag m

addAlternate :: [C] -> [Double]
addAlternate [] = []
addAlternate [_] = error "addAlternate needs even number of elements"
addAlternate (x:y:xs) = checkedRealPart (x+y) : addAlternate xs

-- | Remove the most recent beam from the stack.
dropBeam :: BeamStack -> BeamStack
dropBeam (BeamStack m) = BeamStack (reduceMat m)

-- | Return the number of beams in a 'BeamStack'.
numBeams :: BeamStack -> Int
numBeams (BeamStack m)
    = let (p,_) = size m
      in p `div` 2

-- | Interchange the two most recent beams on the stack.
flipBeams :: BeamStack -> BeamStack
flipBeams (BeamStack m)
    = let (d,_) = size m
          fl = flipMat d
      in BeamStack $ fl <> m <> tr fl

flipMat :: Int -> Matrix C
flipMat d = bigM d (fromLists [[0,0,1,0]
                              ,[0,0,0,1]
                              ,[1,0,0,0]
                              ,[0,1,0,0]])

-- Turn a 2x2 into a dxd.
bigM2 :: Int -> Matrix C -> Matrix C
bigM2 d m
    | d < 2      = error "bigM2 requires d >= 2"
    | odd d      = error "bigM2 requires even d"
    | otherwise  = fromLists $ map (++ [0,0]) (toLists (ident (d-2)))
                   ++ map (replicate (d-2) 0 ++) (toLists m)

-- Turn a 4x4 into a dxd.
bigM :: Int -> Matrix C -> Matrix C
bigM d m
    | d < 4      = error "bigM requires d >= 4"
    | odd d      = error "bigM requires even d"
    | otherwise  = fromLists $ map (++ [0,0,0,0]) (toLists (ident (d-4)))
                   ++ map (replicate (d-4) 0 ++) (toLists m)

s :: Double -> Double -> Matrix C
s theta phi = kronecker (u `couter` u) (np theta phi `couter` np theta phi)
            + kronecker (l `couter` u) (nm theta phi `couter` nm theta phi)
            + kronecker (u `couter` l) (nm theta phi `couter` nm theta phi)
            + kronecker (l `couter` l) (np theta phi `couter` np theta phi)

u :: Vector C
u = zp

l :: Vector C
l = zm

-- | Given angles describing the orientation of the splitter,
--   removes an incoming beam from the stack and replaces
--   it with two beams, a spin-up and a spin-down beam.
--   The spin-down beam is the most recent beam on the stack.
split :: Double -> Double -> BeamStack -> BeamStack
split theta phi (BeamStack m)
    = let m' = extendWithZeros m
          (p,_) = size m'
          ss = bigM p (s theta phi)
      in BeamStack $ ss <> m' <> tr ss

-- | Given angles describing the orientation of the recombiner,
--   returns a single beam from an incoming pair of beams.
recombine :: Double -> Double -> BeamStack -> BeamStack
recombine theta phi (BeamStack m)
    = let (d,_) = size m
          ss = bigM d (s theta phi)
      in dropBeam $ BeamStack $ ss <> m <> tr ss

mag2x2 :: Double -> Double -> Double -> Matrix C
mag2x2 theta phi omegaT
    = let z = iC * (omegaT :+ 0) / 2
          np' = np theta phi
          nm' = nm theta phi
      in scale (exp   z ) (np' `couter` np')
       + scale (exp (-z)) (nm' `couter` nm')

-- | Given angles describing the direction of a
--   uniform magnetic field, and given an angle
--   describing the product of the Larmor frequency
--   and the time, return an output beam from an
--   input beam.
applyBField :: Double -> Double -> Double -> BeamStack -> BeamStack
applyBField theta phi omegaT (BeamStack m)
    = let (d,_) = size m
          uu = bigM2 d (mag2x2 theta phi omegaT)
      in BeamStack $ uu <> m <> tr uu

-----------------------
-- Derived functions --
-----------------------

-- | A Stern-Gerlach splitter in the x direction.
splitX :: BeamStack -> BeamStack
splitX = split (pi/2) 0

-- | A Stern-Gerlach splitter in the y direction.
splitY :: BeamStack -> BeamStack
splitY = split (pi/2) (pi/2)

-- | A Stern-Gerlach splitter in the z direction.
splitZ :: BeamStack -> BeamStack
splitZ = split 0 0

-- | Given an angle in radians
--   describing the product of the Larmor frequency
--   and the time, apply a magnetic in the x direction
--   to the most recent beam on the stack.
applyBFieldX :: Double -> BeamStack -> BeamStack
applyBFieldX = applyBField (pi/2) 0

-- | Given an angle in radians
--   describing the product of the Larmor frequency
--   and the time, apply a magnetic in the y direction
--   to the most recent beam on the stack.
applyBFieldY :: Double -> BeamStack -> BeamStack
applyBFieldY = applyBField (pi/2) (pi/2)

-- | Given an angle in radians
--   describing the product of the Larmor frequency
--   and the time, apply a magnetic in the z direction
--   to the most recent beam on the stack.
applyBFieldZ :: Double -> BeamStack -> BeamStack
applyBFieldZ = applyBField 0 0

-- | A Stern-Gerlach recombiner in the x direction.
recombineX :: BeamStack -> BeamStack
recombineX = recombine (pi/2) 0

-- | A Stern-Gerlach recombiner in the y direction.
recombineY :: BeamStack -> BeamStack
recombineY = recombine (pi/2) (pi/2)

-- | A Stern-Gerlach recombiner in the z direction.
recombineZ :: BeamStack -> BeamStack
recombineZ = recombine 0 0

-- | Filter for spin-up particles in the x direction.
xpFilter :: BeamStack -> BeamStack
xpFilter = dropBeam . splitX

-- | Filter for spin-down particles in the x direction.
xmFilter :: BeamStack -> BeamStack
xmFilter = dropBeam . flipBeams . splitX

-- | Filter for spin-up particles in the z direction.
zpFilter :: BeamStack -> BeamStack
zpFilter = dropBeam . splitZ

-- | Filter for spin-down particles in the z direction.
zmFilter :: BeamStack -> BeamStack
zmFilter = dropBeam . flipBeams . splitZ

ext: hs

max_stars_repo_name: walck/learn-physics

max_issues_repo_name: walck/learn-physics

max_forks_repo_name: walck/learn-physics


text:
module Isumi.Math.FFT.Internal.Bruun2D
  ( fft2dBruun
  ) where

import           Data.Complex
import qualified Data.Vector.Generic           as GV
import           Isumi.Math.FFT.Internal.Bruun (fftBruun, fftBruunC)
import qualified Numeric.LinearAlgebra         as LA
import           Numeric.LinearAlgebra.Data    (Matrix)

fft2dBruun :: Matrix Double -> Maybe (Matrix (Complex Double))
fft2dBruun m = do
  rows' <- traverse (fftBruun . GV.convert) (LA.toRows m)
  cols' <- traverse (fftBruunC . GV.convert) . LA.toColumns . LA.fromRows
           $ fmap GV.convert rows'
  pure $ LA.fromColumns (fmap GV.convert cols')

max_stars_repo_name: IsumiF/fft-bruun

max_issues_repo_name: IsumiF/fft-bruun

max_forks_repo_name: IsumiF/fft-bruun


text:
import Data.Complex

main :: IO ()
main = do let res = dft $ enumerate $ map (\x -> x :+ 0) [1..100]
       	  putStrLn $ show res

enumerate :: Integral b => [a] -> [(a, b)]
enumerate l = zip l [0..]

dft :: (RealFloat a, Integral b) => [(Complex a, b)] -> [Complex a]
dft x = map (\(_, k) -> singleDFT k x) x

singleDFT :: (RealFloat a, Integral b) => b -> [(Complex a, b)] -> Complex a
singleDFT k x = sum $ map (\(e, n) -> (e *) $ fCoeff k n $ length x) $ x

fCoeff :: (RealFloat a, Integral b) => b -> b -> b -> Complex a
fCoeff k n n_tot = exp (0.0 :+ ((-2.0) * pi * (fromIntegral $ k * n) /
       	   	       	       	       	      (fromIntegral n_tot)))

max_stars_repo_name: fredmorcos/attic

max_issues_repo_name: fredmorcos/attic

max_forks_repo_name: fredmorcos/attic


text:
{-# LANGUAGE OverloadedStrings #-}

module Simulation where

import DigitalAGC
import LutLogTests (sfToDouble, ufToDouble)

import Clash.Prelude (
  UFixed, Unsigned, Signed, SFixed(..), Vec(..), Nat, KnownDomain(..), HiddenClockResetEnable, SNat(..),
  createDomain, vInitBehavior, vName, vSystem, InitBehavior(..), Signal,
  simulate, System,
  shiftR,
  register,
  sf, bundle, df,
  riseEvery,
  fLitR, KnownNat
  )

import Prelude
import Data.Maybe (fromMaybe)
import Data.List (group)

import Lucid
import Lucid.Html5
import Graphics.Plotly hiding (sort)
import qualified Graphics.Plotly.Base as GPB
import Graphics.Plotly.Simple
import Graphics.Plotly.Lucid
import Lens.Micro
import System.Process
import Data.Array
import Data.Foldable

import qualified Data.Complex as DC
import Numeric.Transform.Fourier.FFT
import DSP.Basic (linspace)

import DataFlow.Extra
import Test.QuickCheck

import Statistics.Sample
import qualified Data.Vector.Unboxed as V
--import Data.Fixed

--dec :: Int -> [a] -> [a]
--dec n [] = []
--dec n (x : xs) = x : dec n (drop (n-1) xs)

interp :: Int -> [a] -> [a]
interp n [] = []
interp n (x : xs) = replicate n x ++ interp n xs

chunks :: Int -> [a] -> [[a]]
chunks n [] = []
chunks n xs = take n xs : chunks n (drop n xs)

quantiseS :: forall i f . (KnownNat i, KnownNat f) => SNat i -> SNat f -> Double -> Double
quantiseS _ _ x = sfToDouble (fLitR x :: SFixed i f)

quantiseU :: forall i f . (KnownNat i, KnownNat f) => SNat i -> SNat f -> Double -> Double
quantiseU _ _ x = ufToDouble (fLitR x :: UFixed i f)

d0 = SNat :: SNat 0

iSig   = SNat :: SNat 26
iGain  = SNat :: SNat 10
fGain  = SNat :: SNat 10
iAlpha = SNat :: SNat 1
fAlpha = SNat :: SNat 6
iRef   = SNat :: SNat 3
fRef   = SNat :: SNat 12
iLog   = SNat :: SNat 4
fLog   = SNat :: SNat 22

simAgc :: Double -> Double -> Double -> [(Double, Double)] -> [(Double, Double,Double)]
simAgc window ref alpha iqs =
  let window' = 2**window
      powers blk = map (\(i,q)-> quantiseU iSig d0 . sqrt $ i^2 + q^2) blk
      smooth blk = quantiseU iSig d0 $ (/window') $ sum  blk
      alg state x = quantiseS iLog fLog (
                      (quantiseS iLog fLog state) -
                      ((quantiseS iLog fLog $ logBase 10 (smooth $ powers x)) - (quantiseU iRef fRef ref)) *
                      (quantiseU iAlpha fAlpha alpha))
      gains = scanl (\state x ->
                quantiseS iLog fLog $
                alg state (map (\(i,q)->(
                  quantiseS iSig d0 $ (i*) $ quantiseU iGain fGain $ 10**state,
                  quantiseS iSig d0 $ (q*) $ quantiseU iGain fGain $ 10**state)
                  ) x)
              ) 0
              (init $ chunks (round window') iqs)
  in zipWith (\(i,q) s -> let g = quantiseU iGain fGain $ 10**s
                          in (g, quantiseS iSig d0 $ i*g, quantiseS iSig d0 $ q*g))
     iqs (interp (round window') gains)

-- My floating point simulation shows us that we _should_ be ok with our
-- approach... make some tests to verify the log and exp DF behaviour.

createDomain vSystem{vName="SyncDefined", vInitBehavior=Defined}

simInput = let iqs = map (\(i,q)->(shiftR i 1, shiftR q 1)) $ sinInputComplex 1 0.01
           in  take 3000 iqs ++ map (\(i,q)->(shiftR i 4, shiftR q 4)) iqs

-- from 1.25 to 0.75... our ref is basically divided by 4, so this makes sense for a 100x range
-- we have 0.5 range with 6 frac bits => resolution of 1/64, so 32 sensible settings for our signal height.
-- This is all with alpha = 0.9375 and a window of 7
-- Window doesn't really affect this too much
-- Alpha seems to?

-- we have a 2.6 log word ( really 4.4)

-- a = 0.9375 and window = 128 ~ 1.5k         1.5k
-- a = 0.5                       3k           3.5k
-- a = 0.25                      4.5k         5k
-- a = 0.125                     7k           8k
-- a = 0.0625                                 10k

--   alpha             recovery samples   diff      log10  log diff
--        1; 5         2                  ----      0.3    --------
--      0.5; 7.3       5.3                3.3       0.72   0.42
--     0.25; 12        9                  3.7       0.95   0.23
--    0.125; 20        17                 8         1.23   0.27
--   0.0625; 33        31                 14        1.49   0.26

-- ref = 1.13 for 32k
-- for 1000, log 10
--sim = let ref = 1.05 :: UFixed 2 10
--          window = 7 :: Unsigned 5
--          alpha = 2.0 :: UFixed 1 6
--          fLog = Clash.d6
--          fGain = Clash.d6
--          out_gain = Clash.simulate @System (uncurry (digiAgcMult (pure window) (pure ref) (pure alpha)). unbundle) simInput
--      in map (zip [1..]) [
--                           map (fromIntegral . (\(_,x,_)->x)) out_gain
--                         , map (fromIntegral . (\(_,_,x)->x)) out_gain
--                         , map (ufToDouble   . (\(x,_,_)->x)) out_gain
--                         ]

{-
Just thinking about sensible wordlengths for RFSoC...

Our input is always 16 bits.

==========

This will likely be after decimation, used inside the user's RX digital logic so we can expect deal with sampling rates well below the 4 GSamples/s... say probably a max of 512 MHz.

How many window bits are needed for a period of 1 ms?

1e-3 * 512e6 = 18.9 ... let's just say 19

==========

Fullscale signal after log will be...

after ID; we need (nWindow + nSig) bits (35!!!)
after log10; we need (1+log2 35) bits (9) and fLog fractional bits, let's say 9.fLog
after sub with ref; we need 10.fLog
after mul with alpha we need (10+ia).(fLog+fa)

{ What should alpha's range be to allow a good selection of response times? }

After the antilog, our wordlength increases exponentially! We should be super
careful about growing our wordlengths in the log domain.

==========

nSig = 16
nWindow = 19
-}

-- Experimentally found recovery cycles for every tenth decimal step in alpha
-- alphas = [0.1,0.2,...1.9]
recoveryCycles :: [(Int, Int)]
recoveryCycles = [(1,56),(2,40),(3,31),(4,26),(5,22),(6,19),(7,17),(8,15),(9,14),(10,13),(11,12),(12,11),(13,10),(14,10),(15,9),(16,8),(17,8),(18,7),(19,7)]

getRecoveryCycles :: Double -> Int
getRecoveryCycles = fromMaybe (error "No entry found for that alpha") . flip lookup recoveryCycles. round . (*10)

newtype InGain = InGain Double deriving Show
newtype InStepTime = InStepTime Int deriving Show

instance Arbitrary InStepTime where
  arbitrary = fmap InStepTime $ choose (1000,10000)

instance Arbitrary InGain where
  arbitrary = fmap (InGain . recip . fromIntegral) $ choose (1::Int,63)

steppedInput :: InGain -> InGain -> InStepTime -> [(Signed 16,Signed 16)]
steppedInput (InGain g1) (InGain g2) (InStepTime n) =
  let a = map (\(i,q)->(multD i g1, multD q g1)) . take n $ sinInputComplex 1 0.01
      b = map (\(i,q)->(multD i g2, multD q g2))           $ sinInputComplex 1 0.01
  in a ++ b
  where
  multD s d = fromIntegral . round $ fromIntegral s * d

splitInto n [] = []
splitInto n xs = let (a,b) = splitAt n xs
                 in a : splitInto n b

isSteady :: Int -> Double -> Double -> [Double] -> Bool
isSteady n ref percent = (<=n) . maximum . map length . filter (\a->False == a!!0) . group . map (\x->abs (x-ref)/ref < (percent/100))

simDfLogErr ref alpha xs =
  simulate @System (\x ->
    bundle $ df (dfLogErr (pure ref) (pure alpha)) x (riseEvery (SNat :: SNat 3)) (pure True) )
  (cycle xs) :: [(SFixed 4 22, Bool, Bool)]

{--

1) Can I reproduce the variation in simulation? Maybe generating a repeating pattern with a non-integer multiple of the window size

2) Does this also appear in the floating point simulation, because that doesn't have the delays from cordic stuff.

If no to 1) maybe it's to do with my DataFlow implementation?

I'm seeing spikes for 32 cycles in the output. Why is this?! We should only be seeing 512 sample changes because of the window of 2^9

I've fixed it... Two things now:

  1) Are the log10 and pow10 df units *really* returning consistient values?

  2) We can add a feature to set the error state! We can use that in the python instead of feeding it extra samples. Is there anything else we should flush? Maybe just pad with 32ish samples. Wait. isn't that just the reset functionality that should already exist? I've just mapped the areset pin to axi reg 0, bit 1.

--}

simOutPower g1 g2 n =
  let ref = 4.0 :: UFixed 3 12
      alpha = 1.0 :: UFixed 1 6
      window = 9 :: Unsigned 5
      --inputSig = take (10000 + rec_time*(2^window)) $ steppedInput g1 g2 n
      inputSig = take (6000*5) $ steppedInput g1 g2 n
      rec_time = (2+) . getRecoveryCycles $ ufToDouble alpha
      --out_gain = take 15000 $ simAgc (fromIntegral window) (ufToDouble ref) (ufToDouble alpha) (map (\(i,q)->(fromIntegral i, fromIntegral q)) inputSig)
      --out_pow = map (\(_,i,q)-> sqrt $ (i)**2 + (q)**2) out_gain
      ip x =
               let ip = df (testBufferDF (SNat :: SNat 6000) `seqDF` throttleDF (SNat :: SNat 4) `seqDF` dfAgc (pure window) (pure ref) (pure alpha) (pure True))
                   oR = (riseEvery (SNat :: SNat 3))
                   (y, oV, iR) = ip x (pure True) oR
               in bundle (y, oV, oR) :: Signal System ((UFixed 10 15, (Signed 16, Signed 16)), Bool, Bool)
      --outs = drop 1 . take (10000 + rec_time*(2^window))
      outs = drop 1 . take 6000 . filter (\(_,v,r)->v&&r)
             $ simulate @System ip inputSig
      out_gain = map (\((g,(i,q)), v,r)->(g,i,q)) outs
      out_pow = map (\(_,i,q)-> sqrt $ (fromIntegral i)**2 + (fromIntegral q)**2) out_gain
      out_pow_block = map ((/(2^window)) . sum) $ splitInto (2^window) out_pow
      expected_pow = 10**(ufToDouble ref)
  in (inputSig, out_gain, out_pow, out_pow_block, rec_time, expected_pow)

prop_OutPower :: InGain -> InGain -> InStepTime -> Property
prop_OutPower g1 g2 (InStepTime n) =
  let (_, _, _, out_pow_block, rec_time, expected_pow) = simOutPower g1 g2 (InStepTime n)
  in property . isSteady rec_time expected_pow 1 $ drop (ceiling $ (fromIntegral n) / 2**7) out_pow_block

showTest g1 g2 n =
    renderToFile "/tmp/clash/sim.html" $ doctypehtml_ $ do
    head_ $ do meta_ [charset_ "utf-8"]
               plotlyCDN
               reloadCDN
               styleSheet
    body_ $ do
               toHtml $ plotly "time_iq_in" (traceTime simInputTrace)
                          & layout . title ?~ "Time domain I/Q Input"
               toHtml $ plotly "time_iq_out" (traceTime tData)
                          & layout . title ?~ "Time domain I/Q Output"
               toHtml $ plotly "constl_iq_out" (traceConstl tData)
                          & layout . title ?~ "Constellation I/Q"
                          & layout . width ?~ 600
                          & layout . height ?~ 600
               toHtml $ plotly "time_iq_ctrl" (traceTime ctrlData)
                          & layout . title ?~ "Time domain Control"
  where
  (simInput, sim, out_pow, _, _, _) = simOutPower g1 g2 n
  tData = map (zip [1..]) $ [map (\(_,x,_)->fromIntegral x) sim
                            ,map (\(_,_,x)->fromIntegral x) sim
                            ,map (\(x,_,_)->ufToDouble x)   sim
                            --[map (\(_,x,_)-> x) sim
                            --,map (\(_,_,x)-> x) sim
                            --,map (\(x,_,_)->x)   sim
                            ] :: [[(Double, Double)]]
  ctrlData = [tData !! 2, tData !! 2]
  simInputTrace = map (zip [(1::Double)..] . map (sfToDouble . sf (SNat :: SNat 0)))  [map fst simInput, map snd simInput]

corr :: [(Double, Double)] -> Double
corr = correlation . V.fromList

--prop_corrTest g1 g2 n =
--  let (_, outsig, _, _, recTime, _) = simOutPower g1 g2 n
--      insig = steppedInput (InGain 0.3125) (InGain 0.3125) (InStepTime 0)
--      inI  = map (fromIntegral . fst) insig
--      inQ  = map (fromIntegral . snd) insig
--      outI = map (fromIntegral . (\(_,i,_)->i)) outsig
--      outQ = map (fromIntegral . (\(_,_,q)->q)) outsig
--      correlation = corr $ zip (drop recTime inI) (drop recTime outI)
--  in property $ correlation > 0.9

bl10 n = fromIntegral . ceiling $ logBase 2 (logBase 10 (2**n))
be10 n = fromIntegral . ceiling $ logBase 2 (10 ** (2**n))

bitsCalc :: Double -> Double -> Double -> (Double, Double) -> Double
         -> (Double, (Double, Double), (Double, Double))
bitsCalc window sig fLog (iAlpha, fAlpha) fGain =
  let intgInternal = 2**window + sig
      iLog = bl10 sig
      iGain = be10 iLog - fGain
  in (intgInternal, (iLog,fLog), (iGain, fGain))

sinInput :: Double -> Double -> [Signed 16]
sinInput fs dt = map (fromInteger . round . ((2**15-1) * ) .  sin . (2*pi*fs * )) [0.0, dt..]

sinInputComplex :: Double -> Double -> [(Signed 16, Signed 16)]
sinInputComplex fs dt = zip sins coss
  where sins = map (fromInteger . round . ((2**15-1) * ) .  sin . (2*pi*fs * )) [0.0, dt..]
        coss = map (fromInteger . round . ((2**15-1) * ) .  cos . (2*pi*fs * )) [0.0, dt..]

reloadCDN :: Monad m => HtmlT m ()
reloadCDN = script_ [src_ "http://livejs.com/live.js"] $ toHtml (""::String)

styleSheet :: Monad m => HtmlT m ()
styleSheet = style_ [type_ "text/css" ] $
             toHtml (".svg-container {margin: 0 auto !important;}"::String)

traceTime a =  [trace "I" $ a !! 0
               ,trace "Q" $ a !! 1]
  where trace label points = linePlot points
                               & name ?~ label
                               & (GPB.line ?~ (defLine & lineshape ?~ Hv))

traceConstl points = [scatterPlot $ zip (map snd (points!!0)) (map snd (points!!1))]

startServer = callCommand "mkdir -p /tmp/clash; cd /tmp/clash/; python -m http.server > /dev/null 2>&1 &"

return []
runTests = $quickCheckAll

{- TODO

Include bin files in IP make script

Try find what resolution we need to increase in order to better recover very very low signals.

-}


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

-- |
-- Module    : Statistics.Sample.Histogram
-- Copyright : (c) 2011 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Functions for computing histograms of sample data.

module Statistics.Sample.Histogram
    (
      histogram
    -- * Building blocks
    , histogram_
    , range
    ) where

import Numeric.MathFunctions.Constants (m_epsilon,m_tiny)
import Statistics.Function (minMax)
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Generic.Mutable as GM

-- | /O(n)/ Compute a histogram over a data set.
--
-- The result consists of a pair of vectors:
--
-- * The lower bound of each interval.
--
-- * The number of samples within the interval.
--
-- Interval (bin) sizes are uniform, and the upper and lower bounds
-- are chosen automatically using the 'range' function.  To specify
-- these parameters directly, use the 'histogram_' function.
histogram :: (G.Vector v0 Double, G.Vector v1 Double, Num b, G.Vector v1 b) =>
             Int                -- ^ Number of bins (must be positive).
          -> v0 Double          -- ^ Sample data (cannot be empty).
          -> (v1 Double, v1 b)
histogram numBins xs = (G.generate numBins step, histogram_ numBins lo hi xs)
    where (lo,hi)    = range numBins xs
          step i     = lo + d * fromIntegral i
          d          = (hi - lo) / fromIntegral numBins
{-# INLINE histogram #-}

-- | /O(n)/ Compute a histogram over a data set.
--
-- Interval (bin) sizes are uniform, based on the supplied upper
-- and lower bounds.
histogram_ :: (Num b, RealFrac a, G.Vector v0 a, G.Vector v1 b) =>
              Int
           -- ^ Number of bins.  This value must be positive.  A zero
           -- or negative value will cause an error.
           -> a
           -- ^ Lower bound on interval range.  Sample data less than
           -- this will cause an error.
           -> a
           -- ^ Upper bound on interval range.  This value must not be
           -- less than the lower bound.  Sample data that falls above
           -- the upper bound will cause an error.
           -> v0 a
           -- ^ Sample data.
           -> v1 b
histogram_ numBins lo hi xs0 = G.create (GM.replicate numBins 0 >>= bin xs0)
  where
    bin xs bins = go 0
     where
       go i | i >= len = return bins
            | otherwise = do
         let x = xs `G.unsafeIndex` i
             b = truncate $ (x - lo) / d
         write' bins b . (+1) =<< GM.read bins b
         go (i+1)
       write' bins' b !e = GM.write bins' b e
       len = G.length xs
       d = ((hi - lo) / fromIntegral numBins) * (1 + realToFrac m_epsilon)
{-# INLINE histogram_ #-}

-- | /O(n)/ Compute decent defaults for the lower and upper bounds of
-- a histogram, based on the desired number of bins and the range of
-- the sample data.
--
-- The upper and lower bounds used are @(lo-d, hi+d)@, where
--
-- @d = (maximum sample - minimum sample) / ((bins - 1) * 2)@
--
-- If all elements in the sample are the same and equal to @x@ range
-- is set to @(x - |x|/10, x + |x|/10)@. And if @x@ is equal to 0 range
-- is set to @(-1,1)@. This is needed to avoid creating histogram with
-- zero bin size.
range :: (G.Vector v Double) =>
         Int                    -- ^ Number of bins (must be positive).
      -> v Double               -- ^ Sample data (cannot be empty).
      -> (Double, Double)
range numBins xs
    | numBins < 1 = error "Statistics.Histogram.range: invalid bin count"
    | G.null xs   = error "Statistics.Histogram.range: empty sample"
    | lo == hi    = case abs lo / 10 of
                      a | a < m_tiny -> (-1,1)
                        | otherwise  -> (lo - a, lo + a)
    | otherwise   = (lo-d, hi+d)
  where
    d | numBins == 1 = 0
      | otherwise    = (hi - lo) / ((fromIntegral numBins - 1) * 2)
    (lo,hi)          = minMax xs
{-# INLINE range #-}

max_stars_repo_name: vaerksted/statistics

max_issues_repo_name: vaerksted/statistics

max_forks_repo_name: vaerksted/statistics


text:
{-# LANGUAGE StaticPointers #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE UndecidableSuperClasses #-}

module Control.Distributed.Closure.Static where

import Control.Distributed.Closure hiding (Static(..))
import Data.Binary
import Type.Reflection

import qualified Control.Concurrent
import qualified Data.ByteString
import qualified Data.ByteString.Lazy
import qualified Data.ByteString.Short
import qualified Data.Complex
import qualified Data.Int
import qualified Data.IntSet
import qualified Data.Maybe
import qualified Data.Monoid
import qualified Data.Semigroup
import qualified Data.Version
import qualified Data.Void
import qualified Data.Word
import qualified GHC.Exts
import qualified GHC.Fingerprint.Type
import qualified Numeric.Natural
import qualified Type.Reflection.Unsafe

class (c, Typeable c) => Static c where
  -- dict :: Dict c
  -- dict = Dict
  closureDict :: Closure (Dict c)
  closureDictStatic :: Closure (Dict (Static c))

instance (Static c1, Static c2, Typeable c1, Typeable c2, (c1, c2)) =>
         Static (c1, c2) where
  closureDict = static dict `cap` closureDict `cap` closureDict
    where dict :: Dict c3 -> Dict c4 -> Dict (c3, c4)
          dict Dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic `cap` closureDictStatic
    where dict :: Dict (Static c3) -> Dict (Static c4) -> Dict (Static (c3, c4))
          dict Dict Dict = Dict

-- Static Binary and Static Typeable instances for basic types

instance Static (Binary ()) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Data.ByteString.ByteString) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Data.ByteString.Lazy.ByteString) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Data.ByteString.Short.ShortByteString) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Data.Int.Int16) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Data.Int.Int32) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Data.Int.Int64) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Data.Int.Int8) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Data.IntSet.IntSet) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Data.Monoid.All) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Data.Monoid.Any) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Data.Word.Word16) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Data.Word.Word32) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Data.Word.Word64) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Data.Word.Word8) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Data.Version.Version) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Data.Void.Void) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary GHC.Exts.RuntimeRep) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary GHC.Exts.VecCount) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary GHC.Exts.VecElem) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary GHC.Fingerprint.Type.Fingerprint) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Numeric.Natural.Natural) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Prelude.Bool) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Prelude.Char) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Prelude.Double) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Prelude.Float) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Prelude.Int) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Prelude.Integer) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Prelude.Ordering) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Prelude.Word) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Type.Reflection.SomeTypeRep) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Type.Reflection.TyCon) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Type.Reflection.Unsafe.KindRep) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Binary Type.Reflection.Unsafe.TypeLitSort) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable ()) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Data.ByteString.ByteString) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Data.ByteString.Lazy.ByteString) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Data.ByteString.Short.ShortByteString) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Data.Int.Int16) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Data.Int.Int32) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Data.Int.Int64) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Data.Int.Int8) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Data.IntSet.IntSet) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Data.Monoid.All) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Data.Monoid.Any) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Data.Word.Word16) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Data.Word.Word32) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Data.Word.Word64) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Data.Word.Word8) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Data.Version.Version) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Data.Void.Void) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable GHC.Exts.RuntimeRep) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable GHC.Exts.VecCount) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable GHC.Exts.VecElem) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable GHC.Fingerprint.Type.Fingerprint) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Numeric.Natural.Natural) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Prelude.Bool) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Prelude.Char) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Prelude.Double) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Prelude.Float) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Prelude.Int) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Prelude.Integer) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Prelude.Ordering) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Prelude.Word) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Type.Reflection.SomeTypeRep) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Type.Reflection.TyCon) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Type.Reflection.Unsafe.KindRep) where
  closureDict = static Dict
  closureDictStatic = static Dict

instance Static (Typeable Type.Reflection.Unsafe.TypeLitSort) where
  closureDict = static Dict
  closureDictStatic = static Dict

-- Static Binary and Static Typeable instances for basic type constructors

instance (Static (Typeable a), Static (Binary a)) => Static (Binary [a]) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Binary b) -> Dict (Binary [b])
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic `cap` closureDictStatic
    where dict :: Dict (Static (Typeable b))
               -> Dict (Static (Binary b))
               -> Dict (Static (Binary [b]))
          dict Dict Dict = Dict

instance (Static (Typeable a), Static (Binary a)) =>
         Static (Binary (Data.Complex.Complex a)) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Binary b) -> Dict (Binary (Data.Complex.Complex b))
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic `cap` closureDictStatic
    where dict :: Dict (Static (Typeable b))
               -> Dict (Static (Binary b))
               -> Dict (Static (Binary (Data.Complex.Complex b)))
          dict Dict Dict = Dict

instance (Static (Typeable a), Static (Binary a)) =>
         Static (Binary (Data.Maybe.Maybe a)) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Binary b) -> Dict (Binary (Data.Maybe.Maybe b))
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic `cap` closureDictStatic
    where dict :: Dict (Static (Typeable b))
               -> Dict (Static (Binary b))
               -> Dict (Static (Binary (Data.Maybe.Maybe b)))
          dict Dict Dict = Dict

instance (Static (Typeable a), Static (Binary a)) =>
         Static (Binary (Data.Monoid.Dual a)) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Binary b) -> Dict (Binary (Data.Monoid.Dual b))
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic `cap` closureDictStatic
    where dict :: Dict (Static (Binary b))
               -> Dict (Static (Typeable b))
               -> Dict (Static (Binary (Data.Monoid.Dual b)))
          dict Dict Dict = Dict

instance (Static (Typeable a), Static (Binary a)) =>
         Static (Binary (Data.Monoid.First a)) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Binary b) -> Dict (Binary (Data.Monoid.First b))
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic `cap` closureDictStatic
    where dict :: Dict (Static (Binary b))
               -> Dict (Static (Typeable b))
               -> Dict (Static (Binary (Data.Monoid.First b)))
          dict Dict Dict = Dict

instance (Static (Typeable a), Static (Binary a)) =>
         Static (Binary (Data.Monoid.Last a)) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Binary b) -> Dict (Binary (Data.Monoid.Last b))
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic `cap` closureDictStatic
    where dict :: Dict (Static (Binary b))
               -> Dict (Static (Typeable b))
               -> Dict (Static (Binary (Data.Monoid.Last b)))
          dict Dict Dict = Dict

instance (Static (Typeable a), Static (Binary a)) =>
         Static (Binary (Data.Monoid.Product a)) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Binary b) -> Dict (Binary (Data.Monoid.Product b))
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic `cap` closureDictStatic
    where dict :: Dict (Static (Binary b))
               -> Dict (Static (Typeable b))
               -> Dict (Static (Binary (Data.Monoid.Product b)))
          dict Dict Dict = Dict

instance (Static (Typeable a), Static (Binary a)) =>
         Static (Binary (Data.Monoid.Sum a)) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Binary b) -> Dict (Binary (Data.Monoid.Sum b))
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic `cap` closureDictStatic
    where dict :: Dict (Static (Binary b))
               -> Dict (Static (Typeable b))
               -> Dict (Static (Binary (Data.Monoid.Sum b)))
          dict Dict Dict = Dict

instance (Static (Typeable a), Static (Binary a)) =>
         Static (Binary (Data.Semigroup.Max a)) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Binary b) -> Dict (Binary (Data.Semigroup.Max b))
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic `cap` closureDictStatic
    where dict :: Dict (Static (Typeable b))
               -> Dict (Static (Binary b))
               -> Dict (Static (Binary (Data.Semigroup.Max b)))
          dict Dict Dict = Dict

instance (Static (Typeable a), Static (Binary a)) =>
         Static (Binary (Data.Semigroup.Min a)) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Binary b) -> Dict (Binary (Data.Semigroup.Min b))
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic `cap` closureDictStatic
    where dict :: Dict (Static (Typeable b))
               -> Dict (Static (Binary b))
               -> Dict (Static (Binary (Data.Semigroup.Min b)))
          dict Dict Dict = Dict

instance Static (Typeable a) => Static (Typeable [a]) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Typeable b) -> Dict (Typeable [b])
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic
    where dict :: Dict (Static (Typeable b))
               -> Dict (Static (Typeable [b]))
          dict Dict = Dict

instance Static (Typeable a) =>
         Static (Typeable (Data.Complex.Complex a)) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Typeable b) -> Dict (Typeable (Data.Complex.Complex b))
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic
    where dict :: Dict (Static (Typeable b))
               -> Dict (Static (Typeable (Data.Complex.Complex b)))
          dict Dict = Dict

instance Static (Typeable a) =>
         Static (Typeable (Data.Maybe.Maybe a)) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Typeable b) -> Dict (Typeable (Data.Maybe.Maybe b))
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic
    where dict :: Dict (Static (Typeable b))
               -> Dict (Static (Typeable (Data.Maybe.Maybe b)))
          dict Dict = Dict

instance Static (Typeable a) =>
         Static (Typeable (Data.Monoid.Dual a)) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Typeable b) -> Dict (Typeable (Data.Monoid.Dual b))
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic
    where dict :: Dict (Static (Typeable b))
               -> Dict (Static (Typeable (Data.Monoid.Dual b)))
          dict Dict = Dict

instance Static (Typeable a) =>
         Static (Typeable (Data.Monoid.First a)) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Typeable b) -> Dict (Typeable (Data.Monoid.First b))
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic
    where dict :: Dict (Static (Typeable b))
               -> Dict (Static (Typeable (Data.Monoid.First b)))
          dict Dict = Dict

instance Static (Typeable a) =>
         Static (Typeable (Data.Monoid.Last a)) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Typeable b) -> Dict (Typeable (Data.Monoid.Last b))
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic
    where dict :: Dict (Static (Typeable b))
               -> Dict (Static (Typeable (Data.Monoid.Last b)))
          dict Dict = Dict

instance Static (Typeable a) =>
         Static (Typeable (Data.Monoid.Product a)) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Typeable b) -> Dict (Typeable (Data.Monoid.Product b))
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic
    where dict :: Dict (Static (Typeable b))
               -> Dict (Static (Typeable (Data.Monoid.Product b)))
          dict Dict = Dict

instance Static (Typeable a) =>
         Static (Typeable (Data.Monoid.Sum a)) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Typeable b) -> Dict (Typeable (Data.Monoid.Sum b))
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic
    where dict :: Dict (Static (Typeable b))
               -> Dict (Static (Typeable (Data.Monoid.Sum b)))
          dict Dict = Dict

instance Static (Typeable a) =>
         Static (Typeable (Data.Semigroup.Max a)) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Typeable b) -> Dict (Typeable (Data.Semigroup.Max b))
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic
    where dict :: Dict (Static (Typeable b))
               -> Dict (Static (Typeable (Data.Semigroup.Max b)))
          dict Dict = Dict

instance Static (Typeable a) =>
         Static (Typeable (Data.Semigroup.Min a)) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Typeable b) -> Dict (Typeable (Data.Semigroup.Min b))
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic
    where dict :: Dict (Static (Typeable b))
               -> Dict (Static (Typeable (Data.Semigroup.Min b)))
          dict Dict = Dict

-- Static Typeable without Static Binary instances for basic type constructors

instance Static (Typeable a) =>
         Static (Typeable (Control.Concurrent.MVar a)) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Typeable b)
               -> Dict (Typeable (Control.Concurrent.MVar b))
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic
    where dict :: Dict (Static (Typeable b))
               -> Dict (Static (Typeable (Control.Concurrent.MVar b)))
          dict Dict = Dict

instance Static (Typeable a) => Static (Typeable (Data.Monoid.Endo a)) where
  closureDict = static dict `cap` closureDict
    where dict :: Dict (Typeable b) -> Dict (Typeable (Data.Monoid.Endo b))
          dict Dict = Dict
  closureDictStatic =
    static dict `cap` closureDictStatic
    where dict :: Dict (Static (Typeable b))
               -> Dict (Static (Typeable (Data.Monoid.Endo b)))
          dict Dict = Dict


text:
{-# OPTIONS_GHC -Wall #-}
{-# LANGUAGE Trustworthy #-}
{-# LANGUAGE CPP #-}

{- |
Module      :  Physics.Learn.QuantumMat
Copyright   :  (c) Scott N. Walck 2016-2018
License     :  BSD3 (see LICENSE)
Maintainer  :  Scott N. Walck <walck@lvc.edu>
Stability   :  experimental

This module contains state vectors and matrices
for quantum mechanics.
-}

-- Using only Complex Double here, no cyclotomic

module Physics.Learn.QuantumMat
    (
    -- * Complex numbers
      C
    -- * State Vectors
    , xp
    , xm
    , yp
    , ym
    , zp
    , zm
    , np
    , nm
    , dim
    , scaleV
    , inner
    , norm
    , normalize
    , probVector
    , gramSchmidt
    , conjV
    , fromList
    , toList
    -- * Matrices (operators)
    , sx
    , sy
    , sz
    , scaleM
    , (<>)
    , (#>)
    , (<#)
    , conjugateTranspose
    , fromLists
    , toLists
    , size
    , matrixFunction
    -- * Density matrices
    , couter
    , dm
    , trace
    , normalizeDM
    , oneQubitMixed
    -- * Quantum Dynamics
    , timeEvMat
    , timeEv
    , timeEvMatSpec
    -- * Composition
    , Kronecker(..)
    -- * Measurement
    , possibleOutcomes
    , outcomesProjectors
    , outcomesProbabilities
    -- * Vector and Matrix
    , Vector
    , Matrix
    )
    where

import Numeric.LinearAlgebra
    ( C
    , Vector
    , Matrix
    , Herm
    , iC        -- square root of negative one
    , (><)      -- matrix definition
    , ident
    , scale
    , norm_2
    , inv
    , (<\>)
    , sym
    , eigenvaluesSH
    , eigSH
    , cmap
    , takeDiag
    , conj
    , dot
    , tr
    )
--    , (<>)      -- matrix product (not * !!!!)
--    , (#>)      -- matrix-vector product
--    , fromList  -- vector definition

import qualified Numeric.LinearAlgebra as H
-- because H.outer does not conjugate
import Data.Complex
    ( Complex(..)
    , magnitude
    , realPart
    )
#if MIN_VERSION_base(4,11,0)
import Prelude hiding ((<>))
#endif

-- | The state resulting from a measurement of
--   spin angular momentum in the x direction
--   on a spin-1/2 particle
--   when the result of the measurement is hbar/2.
xp :: Vector C
xp = normalize $ fromList [1, 1]

-- | The state resulting from a measurement of
--   spin angular momentum in the x direction
--   on a spin-1/2 particle
--   when the result of the measurement is -hbar/2.
xm :: Vector C
xm = normalize $ fromList [1, -1]

-- | The state resulting from a measurement of
--   spin angular momentum in the y direction
--   on a spin-1/2 particle
--   when the result of the measurement is hbar/2.
yp :: Vector C
yp = normalize $ fromList [1, iC]

-- | The state resulting from a measurement of
--   spin angular momentum in the y direction
--   on a spin-1/2 particle
--   when the result of the measurement is -hbar/2.
ym :: Vector C
ym = normalize $ fromList [1, -iC]

-- | The state resulting from a measurement of
--   spin angular momentum in the z direction
--   on a spin-1/2 particle
--   when the result of the measurement is hbar/2.
zp :: Vector C
zp = normalize $ fromList [1, 0]

-- | The state resulting from a measurement of
--   spin angular momentum in the z direction
--   on a spin-1/2 particle
--   when the result of the measurement is -hbar/2.
zm :: Vector C
zm = normalize $ fromList [0, 1]

-- | The state resulting from a measurement of
--   spin angular momentum in the direction
--   specified by spherical angles theta (polar angle)
--   and phi (azimuthal angle)
--   on a spin-1/2 particle
--   when the result of the measurement is hbar/2.
np :: Double -> Double -> Vector C
np theta phi = fromList [ cos (theta/2) :+ 0
                        , exp(0 :+ phi) * (sin (theta/2) :+ 0) ]

-- | The state resulting from a measurement of
--   spin angular momentum in the direction
--   specified by spherical angles theta (polar angle)
--   and phi (azimuthal angle)
--   on a spin-1/2 particle
--   when the result of the measurement is -hbar/2.
nm :: Double -> Double -> Vector C
nm theta phi = fromList [ sin (theta/2) :+ 0
                        , -exp(0 :+ phi) * (cos (theta/2) :+ 0) ]

-- | Dimension of a vector.
dim :: Vector C -> Int
dim = H.size

-- | Scale a complex vector by a complex number.
scaleV :: C -> Vector C -> Vector C
scaleV = scale

-- | Complex inner product.  First vector gets conjugated.
inner :: Vector C -> Vector C -> C
inner = dot

-- | Length of a complex vector.
norm :: Vector C -> Double
norm = norm_2

-- | Return a normalized version of a given state vector.
normalize :: Vector C -> Vector C
normalize v = scale (1 / norm_2 v :+ 0) v

-- | Return a vector of probabilities for a given state vector.
probVector :: Vector C       -- ^ state vector
           -> Vector Double  -- ^ vector of probabilities
probVector = cmap (\c -> magnitude c**2)

-- | Conjugate the entries of a vector.
conjV :: Vector C -> Vector C
conjV = conj

-- | Construct a vector from a list of complex numbers.
fromList :: [C] -> Vector C
fromList = H.fromList

-- | Produce a list of complex numbers from a vector.
toList :: Vector C -> [C]
toList = H.toList

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

-- | The Pauli X matrix.
sx :: Matrix C
sx = (2><2) [ 0, 1
            , 1, 0 ]

-- | The Pauli Y matrix.
sy :: Matrix C
sy = (2><2) [  0, -iC
            , iC,   0 ]

-- | The Pauli Z matrix.
sz :: Matrix C
sz = (2><2) [ 1,  0
            , 0, -1 ]

-- | Scale a complex matrix by a complex number.
scaleM :: C -> Matrix C -> Matrix C
scaleM = scale

-- | Matrix product.
(<>) :: Matrix C -> Matrix C -> Matrix C
(<>) = (H.<>)

-- | Matrix-vector product.
(#>) :: Matrix C -> Vector C -> Vector C
(#>) = (H.#>)

-- | Vector-matrix product
(<#) :: Vector C -> Matrix C -> Vector C
(<#) = (H.<#)

-- | Conjugate transpose of a matrix.
conjugateTranspose :: Matrix C -> Matrix C
conjugateTranspose = tr

-- | Construct a matrix from a list of lists of complex numbers.
fromLists :: [[C]] -> Matrix C
fromLists = H.fromLists

-- | Produce a list of lists of complex numbers from a matrix.
toLists :: Matrix C -> [[C]]
toLists = H.toLists

-- | Size of a matrix.
size :: Matrix C -> (Int,Int)
size = H.size

-- | Apply a function to a matrix.
--   Assumes the matrix is a normal matrix (a matrix
--   with an orthonormal basis of eigenvectors).
matrixFunction :: (C -> C) -> Matrix C -> Matrix C
matrixFunction f m
    = let (valv,vecm) = H.eig m
          fvalv = fromList [f val | val <- toList valv]
      in vecm <> H.diag fvalv <> tr vecm

----------------------
-- Density Matrices --
----------------------

-- | Complex outer product
couter :: Vector C -> Vector C -> Matrix C
couter v w = v `H.outer` conj w

-- | Build a pure-state density matrix from a state vector.
dm :: Vector C -> Matrix C
dm cvec = cvec `couter` cvec

-- | Trace of a matrix.
trace :: Matrix C -> C
trace = sum . toList . takeDiag

-- | Normalize a density matrix so that it has trace one.
normalizeDM :: Matrix C -> Matrix C
normalizeDM rho = scale (1 / trace rho) rho

-- | The one-qubit totally mixed state.
oneQubitMixed :: Matrix C
oneQubitMixed = normalizeDM $ ident 2

----------------------
-- Quantum Dynamics --
----------------------

-- | Given a time step and a Hamiltonian matrix,
--   produce a unitary time evolution matrix.
--   Unless you really need the time evolution matrix,
--   it is better to use 'timeEv', which gives the
--   same numerical results without doing an explicit
--   matrix inversion.  The function assumes hbar = 1.
timeEvMat :: Double -> Matrix C -> Matrix C
timeEvMat dt h
    = let ah = scale (0 :+ dt / 2) h
          (l,m) = size h
          n = if l == m then m else error "timeEv needs square Hamiltonian"
          identity = ident n
      in inv (identity + ah) <> (identity - ah)

-- | Given a time step and a Hamiltonian matrix,
--   advance the state vector using the Schrodinger equation.
--   This method should be faster than using 'timeEvMat'
--   since it solves a linear system rather than calculating
--   an inverse matrix.  The function assumes hbar = 1.
timeEv :: Double -> Matrix C -> Vector C -> Vector C
timeEv dt h v
    = let ah = scale (0 :+ dt / 2) h
          (l,m) = size h
          n = if l == m then m else error "timeEv needs square Hamiltonian"
          identity = ident n
      in (identity + ah) <\> ((identity - ah) #> v)

-- | Given a Hamiltonian matrix, return a function from time
--   to evolution matrix.  Uses spectral decomposition.
--   Assumes hbar = 1.
timeEvMatSpec :: Matrix C -> Double -> Matrix C
timeEvMatSpec m t = matrixFunction (\h -> exp(-iC * h * (t :+ 0))) m

-----------------
-- Composition --
-----------------

class Kronecker a where
    kron :: a -> a -> a

instance H.Product t => Kronecker (Vector t) where
    kron v1 v2 = H.fromList [c1 * c2 | c1 <- H.toList v1, c2 <- H.toList v2]

instance H.Product t => Kronecker (Matrix t) where
    kron = H.kronecker

-----------------
-- Measurement --
-----------------

-- | The possible outcomes of a measurement
--   of an observable.
--   These are the eigenvalues of the matrix
--   of the observable.
possibleOutcomes :: Matrix C -> [Double]
possibleOutcomes observable
    = H.toList $ eigenvaluesSH (sym observable)

-- From a Hermitian matrix, a list of pairs of eigenvalues and eigenvectors.
valsVecs :: Herm C -> [(Double,Vector C)]
valsVecs h = let (valv,m) = eigSH h
                 vals = H.toList valv
                 vecs = map (conjV . fromList) $ toLists (conjugateTranspose m)
             in zip vals vecs

-- From a Hermitian matrix, a list of pairs of eigenvalues and projectors.
valsPs :: Herm C -> [(Double,Matrix C)]
valsPs h = [(val,couter vec vec) | (val,vec) <- valsVecs h]

combineFst :: (Eq a, Num b) => [(a,b)] -> [(a,b)]
combineFst [] = []
combineFst [p] = [p]
combineFst ((x1,m1):(x2,m2):ps)
    = if x1 == x2
      then combineFst ((x1,m1+m2):ps)
      else (x1,m1):combineFst ((x2,m2):ps)

-- | Given an obervable, return a list of pairs
--   of possible outcomes and projectors
--   for each outcome.
outcomesProjectors :: Matrix C -> [(Double,Matrix C)]
outcomesProjectors m = combineFst (valsPs (sym m))

-- | Given an observable and a state vector, return a list of pairs
--   of possible outcomes and probabilites
--   for each outcome.
outcomesProbabilities :: Matrix C -> Vector C -> [(Double,Double)]
outcomesProbabilities m v
    = [(a,realPart (inner v (p #> v))) | (a,p) <- outcomesProjectors m]

------------------
-- Gram-Schmidt --
------------------

-- | Form an orthonormal list of complex vectors
--   from a linearly independent list of complex vectors.
gramSchmidt :: [Vector C] -> [Vector C]
gramSchmidt [] = []
gramSchmidt (v:vs) = let nvs = gramSchmidt vs
                         nv = normalize (v - sum [scale (inner w v) w | w <- nvs])
                     in nv:nvs

-- To Do
--   Generate higher spin operators and state vectors
--   eigenvectors
--   projection operators


text:
{-# LANGUAGE CPP              #-}
{-# LANGUAGE DataKinds        #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE GADTs            #-}
{-# LANGUAGE TemplateHaskell  #-}
{-# LANGUAGE TypeOperators    #-}

module Test.Grenade.Sys.Utils where

import           Numeric.LinearAlgebra        hiding (R, konst, randomVector,
                                               uniformSample, (===))
import qualified Numeric.LinearAlgebra        as LA
import qualified Numeric.LinearAlgebra.Data   as D
import           Numeric.LinearAlgebra.Static (L, R)
import qualified Numeric.LinearAlgebra.Static as H

import           Hedgehog
import qualified Hedgehog.Gen                 as Gen
import qualified Hedgehog.Range               as Range

import           Test.Hedgehog.Compat
import           Test.Hedgehog.Hmatrix

import           GHC.TypeLits
import           Data.List
import           Data.Maybe
import           Data.Function
import           Data.Either
import           Control.Monad                (guard)
import           System.FilePath
import qualified Data.ByteString.Lazy as BS
import           Data.Binary

import           Grenade
import           Grenade.Utils.PascalVoc
import           Grenade.Utils.ImageNet

imagesDir :: FilePath
imagesDir = (takeDirectory __FILE__) </> "Images"

loadSerializedImage :: (KnownNat d, KnownNat c, KnownNat (d * c), c ~ 3) => FilePath -> IO (Maybe (S ('D3 d d c)))
loadSerializedImage path = do
  bs <- BS.readFile path
  let mat = decode bs :: LA.Matrix Double
  return $ S3D <$> (H.create (D.cmap doubleToRealNum mat))

loadSerializedChannel :: (KnownNat d, KnownNat c, KnownNat (d * c), c ~ 1) => FilePath -> IO (Maybe (S ('D3 d d c)))
loadSerializedChannel path = do
  bs <- BS.readFile path
  let mat = decode bs :: LA.Matrix Double
  return $ S3D <$> (H.create (D.cmap doubleToRealNum mat))

-- Performs linear regression on data points, where
-- the x-values consist of the data set [1..n]
linearRegression :: Int -> [Double] -> Double
linearRegression n ys = m
  where
    n'    = fromIntegral n :: Double
    xs    = [1..n] :: [Int]
    x     = (n' + 1) / 2
    x2    = (fromIntegral $ sum (map (^ 2) xs)) / n'
    xyNum = (sum $ zipWith (\l i -> l * (fromIntegral i)) ys xs) :: Double
    xy    = xyNum / n' :: Double
    y     = (sum ys) / n'
    m     = (xy - (x * y)) / (x2 - (x * x)) :: Double

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
module Optimization where

import Numeric.LinearAlgebra

class Optimize a where
  --optimizeInit :: a -> Double -> Double -> (Int, Int)  -> a
  paramUpdate :: a -> a

data SGD = SGD {
  sgdLearningRate :: Double,
  sgdDecay :: Double,
  sgdMomentum :: Double,
  sgdParams :: [Matrix R],
  sgdGparams :: [Matrix R],
  sgdDeltaPre :: [Matrix R],
  sgdBatchSize :: Double
                 }

instance Optimize SGD where
  paramUpdate sgd =
    let learningRate = sgdLearningRate sgd in
    let decay = sgdDecay sgd in
    let momentum = sgdMomentum sgd in
    let params = sgdParams sgd in
    let gparams = sgdGparams sgd in
    let deltaPre = sgdDeltaPre sgd in
    let batchSize = sgdBatchSize sgd in
    let _sgd (px:pxs) (gx: gxs) (dx:dxs) paUpdate deUpdate gaZeros  =
          let r = rows gx in
          let c = cols gx in
          let lrMat = (r >< c) [learningRate..]:: (Matrix R) in
          let batchSizeMat = (r >< c) [batchSize..] :: (Matrix R) in
          let momMat = (r >< c) [momentum..] :: (Matrix R) in
          let gUp = gx * lrMat / batchSizeMat in
          let ugd = momMat * dx - gx in
          let pu = px + ugd in _sgd pxs gxs dxs (paUpdate++[pu]) (deUpdate++[ugd]) (gaZeros++[(r >< c) [0..]])
        _sgd [] [] [] paUpdate deUpdate gaZeros = SGD {
          sgdLearningRate=learningRate,
          sgdDecay=decay,
          sgdMomentum=momentum,
          sgdParams=paUpdate,
          sgdGparams=gaZeros,
          sgdDeltaPre=deUpdate,
          sgdBatchSize=0
          } in
    _sgd params gparams deltaPre [] [] []

data Adadelta = Adadelta {
  adaDecay :: Double,
  adaEpsilon :: Double,
  adaBatchSize :: Double,
  adaParams :: [Matrix R],
  adaGparams :: [Matrix R],
  adaAccGrad :: [Matrix R],
  adaAccDelta :: [Matrix R]
                         }

instance Optimize Adadelta where
  paramUpdate ada =
    let decay = adaDecay ada in
    let epsilon = adaEpsilon ada in
    let batchSize = adaBatchSize ada in
    let params = adaParams ada in
    let gparams = adaGparams ada in
    let accGrad = adaAccGrad ada in
    let accDelta = adaAccDelta ada in
    let _ada (px:pxs) (gx:gxs) (accgx:accgxs) (accdx:accdxs) paUpdate gradUpdate deltaUpdate gaZeros =
          let r = rows px in
          let c = cols px in
          let decayMat = (r >< c) [decay ..] :: (Matrix R) in
          let epsilonMat = (r >< c) [epsilon ..] :: (Matrix R) in
          let batchSizeMat = (r >< c) [batchSize ..] :: (Matrix R) in
          let oneMat = (r >< c) [1.0 ..]:: (Matrix R) in
          let gUp = accgx / batchSizeMat in
          let gradUp = decayMat * accgx + (oneMat - decayMat) * gUp * gUp in
          let ugd = (cmap (\x -> -(sqrt x)) (accdx + epsilonMat)) /  (cmap (\x -> sqrt x) (accgx + epsilonMat)) * gUp in
          let deltaUp = decayMat * accdx + (oneMat - decayMat) * ugd * ugd in
          let pu = px + ugd in
          _ada pxs gxs accgxs accdxs (pu:paUpdate) (gradUp:gradUpdate) (deltaUp:deltaUpdate) (((r >< c)[0.0 ..]):gaZeros)
        _ada [] [] [] [] paUpdate gradUpdate deltaUpdate gaZeros = Adadelta {
          adaDecay=decay,
          adaEpsilon=epsilon,
          adaBatchSize=batchSize,
          adaParams=(reverse paUpdate),
          adaGparams=(reverse gaZeros),
          adaAccGrad=(reverse gradUpdate),
          adaAccDelta=(reverse deltaUpdate)
                                                                        }
    in _ada params gparams accGrad accDelta [] [] [] []

max_stars_repo_name: neutronest/vortex

max_issues_repo_name: neutronest/vortex

max_forks_repo_name: neutronest/vortex


text:
{-# OPTIONS_GHC -Wall #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE InstanceSigs #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE TypeFamilies #-}
module DataRep where

import Prelude hiding ((.), id, curry, uncurry)
import Control.Exception (assert)
import Data.List
import Data.Proxy
import qualified Data.Vector.Generic as VG
import Foreign.Storable
import qualified Numeric.LinearAlgebra as HM

import qualified Base
import Base hiding (Scalable (..))

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

class (Fractional (Scalar a), Storable (Scalar a), HM.Numeric (Scalar a), Additive a) => VectorSpace a where
  type Scalar a
  dim :: Proxy a -> Int
  toVector :: a -> HM.Vector (Scalar a)
  fromVector :: HM.Vector (Scalar a) -> a

  scale :: Scalar a -> a -> a
  scale a = fromVector . HM.scale a . toVector

linComb :: VectorSpace a => [(a, Scalar a)] -> a
linComb xs = foldl' (.+.) zero [scale s a | (a,s) <- xs]

instance VectorSpace Double where
  type Scalar Double = Double
  scale = (*)
  dim _ = 1
  toVector = VG.singleton
  fromVector = (VG.! 0)

instance (VectorSpace a, VectorSpace b, Scalar a ~ Scalar b) => VectorSpace (a, b) where
  type Scalar (a, b) = Scalar a
  scale s (a,b) = (scale s a, scale s b)
  dim _ = dim (Proxy :: Proxy a) + dim (Proxy :: Proxy b)
  toVector (a,b) = toVector a <> toVector b
  fromVector v =
    case VG.splitAt (dim (Proxy :: Proxy a)) v of
      (v1, v2) -> (fromVector v1, fromVector v2)

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

newtype LinMap s a b = LinMap (HM.Matrix s)

{-
infixr 0 ⊸

data a ⊸ b where
  LinMap :: (Scalar a ~ Scalar b) => !(HM.Matrix (Scalar a)) -> a ⊸ b

と書けると格好良いが、スカラー型の異なるベクトル空間の直積とかで困ったことになる。
-}

instance Category (LinMap s) where
  type Obj (LinMap s) a = (VectorSpace a, Scalar a ~ s)

  id :: forall a. Obj (LinMap s) a => LinMap s a a
  id = LinMap $ HM.ident $ dim (Proxy :: Proxy a)

  LinMap m1 . LinMap m2 = LinMap (m1 HM.<> m2)

instance ToFun (LinMap s) where
  toFun (LinMap m) a = fromVector (m HM.#> toVector a)

instance (VectorSpace a, VectorSpace b, Scalar a ~ s, Scalar b ~ s) => Additive (LinMap s a b) where
  zero = LinMap $ HM.konst 0 (m, n)
    where
      m = dim (Proxy :: Proxy a)
      n = dim (Proxy :: Proxy b)
  LinMap m1 .+. LinMap m2 = LinMap $ HM.add m1 m2

instance (VectorSpace a, VectorSpace b, Scalar a ~ s, Scalar b ~ s) => VectorSpace (LinMap s a b) where
  type Scalar (LinMap s a b) = s
  scale s (LinMap m) = LinMap (HM.scale s m)
  dim _ = dim (Proxy :: Proxy a) * dim (Proxy :: Proxy b)
  toVector (LinMap m) = HM.flatten m
  fromVector v = LinMap (HM.reshape (dim (Proxy :: Proxy b)) v)

instance Monoidal (LinMap s) where
  LinMap m1 >< LinMap m2 = LinMap $ HM.diagBlock [m1, m2]
  monObj = ObjR

instance Cartesian (LinMap s) where
  exl :: forall a b. (Obj (LinMap s) a, Obj (LinMap s) b) => LinMap s (a,b) a
  exl = LinMap $ HM.ident m HM.||| HM.konst 0 (m, n)
    where
      m = dim (Proxy :: Proxy a)
      n = dim (Proxy :: Proxy b)

  exr :: forall a b. (Obj (LinMap s) a, Obj (LinMap s) b) => LinMap s (a,b) b
  exr = LinMap $ HM.konst 0 (n, m) HM.||| HM.ident n
    where
      m = dim (Proxy :: Proxy a)
      n = dim (Proxy :: Proxy b)

  dup :: forall a. Obj (LinMap s) a => LinMap s a (a,a)
  dup = LinMap $
      HM.ident n
      HM.===
      HM.ident n
    where
      n = dim (Proxy :: Proxy a)

instance Cocartesian (LinMap s) where
  inl :: forall a b. (Obj (LinMap s) a, Obj (LinMap s) b) => LinMap s a (a,b)
  inl = LinMap $
      HM.ident m
      HM.===
      HM.konst 0 (n, m)
    where
      m = dim (Proxy :: Proxy a)
      n = dim (Proxy :: Proxy b)

  inr :: forall a b. (Obj (LinMap s) a, Obj (LinMap s) b) => LinMap s b (a,b)
  inr = LinMap $
      HM.konst 0 (n, m)
      HM.===
      HM.ident m
    where
      m = dim (Proxy :: Proxy a)
      n = dim (Proxy :: Proxy b)

  jam :: forall a. (Obj (LinMap s) a) => LinMap s (a,a) a
  jam = LinMap $ HM.ident n HM.||| HM.ident n
    where
      n = dim (Proxy :: Proxy a)

instance (VectorSpace s, Scalar s ~ s) => Base.Scalable (LinMap s) s where
  scale s = assert (dim (Proxy :: Proxy s) == 1) $ LinMap ((1 HM.>< 1) [s])

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

-- | Dual vector space
newtype Dual a = Dual (LinMap (Scalar a) a (Scalar a))

instance VectorSpace a => Additive (Dual a) where
  zero = Dual $ LinMap $ HM.konst 0 (n, 1)
    where
      n = dim (Proxy :: Proxy a)
  Dual (LinMap m1) .+. Dual (LinMap m2) = Dual $ LinMap $ HM.add m1 m2

instance VectorSpace a => VectorSpace (Dual a) where
  type Scalar (Dual a) = Scalar a
  scale s (Dual (LinMap m)) = Dual (LinMap (HM.scale s m))
  dim _ = dim (Proxy :: Proxy a)
  toVector (Dual (LinMap m)) = HM.flatten m
  fromVector v = Dual (LinMap (HM.reshape 1 v))

toDual :: VectorSpace a => a -> Dual a
toDual a = Dual $ LinMap $ HM.asRow (toVector a)

toDualMap :: forall a b s. (VectorSpace a, VectorSpace b, Scalar a ~ s, Scalar b ~ s) => LinMap s a b -> LinMap s (Dual b) (Dual a)
toDualMap (LinMap m) = LinMap (HM.tr' m)

fromDual :: forall a. VectorSpace a => Dual a -> a
fromDual (Dual (LinMap m)) = fromVector $ HM.flatten m

fromDualMap :: forall a b s. (VectorSpace a, VectorSpace b, Scalar a ~ s, Scalar b ~ s) => LinMap s (Dual b) (Dual a) -> LinMap s a b
fromDualMap (LinMap m) = LinMap (HM.tr' m)

instance (VectorSpace u, VectorSpace s, Scalar u ~ s, Scalar s ~ s) => Base.HasDot (LinMap s) s u where
  dot x = case toDual x of Dual f -> f
  undot f = fromDual (Dual f)

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

testDual = (toFun f (2,1) == 7, toFun f2 (2,1) == 7)
  where
    f :: LinMap Double (Double, Double) Double
    f = LinMap $ HM.asRow (VG.fromList [2,3])
    f2 = fromDualMap (toDualMap f)

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

data a :⊗ b where
  TensorProd :: (Scalar a ~ s, Scalar b ~ s) => HM.Matrix s -> a :⊗ b

instance (VectorSpace a, VectorSpace b, Scalar a ~ Scalar b) => Additive (a :⊗ b) where
  zero = TensorProd (HM.konst 0 (m, n))
    where
      m = dim (Proxy :: Proxy a)
      n = dim (Proxy :: Proxy b)
  TensorProd m1 .+. TensorProd m2 = TensorProd $ m1 `HM.add` m2

instance (VectorSpace a, VectorSpace b, Scalar a ~ Scalar b) => VectorSpace (a :⊗ b) where
  type Scalar (a :⊗ b) = Scalar a
  scale s (TensorProd m) = TensorProd $ HM.scale s m
  dim _ = dim (Proxy :: Proxy a) * dim (Proxy :: Proxy b)
  toVector (TensorProd m) = HM.flatten m
  fromVector v = TensorProd (HM.reshape (dim (Proxy :: Proxy b)) v)

curry
  :: forall a b c s. (VectorSpace a, VectorSpace b, VectorSpace c, Scalar a ~ s, Scalar b ~ s, Scalar c ~ s)
  => LinMap s (a :⊗ b) c -> LinMap s a (LinMap s b c)
curry (LinMap m) = LinMap $ HM.reshape (nb*nc) $ HM.flatten m
  where
    nb = dim (Proxy :: Proxy b)
    nc = dim (Proxy :: Proxy c)

uncurry
  :: forall a b c s. (VectorSpace a, VectorSpace b, VectorSpace c, Scalar a ~ s, Scalar b ~ s, Scalar c ~ s)
  => LinMap s a (LinMap s b c) -> LinMap s (a :⊗ b) c
uncurry (LinMap m) = LinMap $ HM.reshape nc $ HM.flatten m
  where
    nc = dim (Proxy :: Proxy c)

mapTensor
  :: forall a b c d s.
     (VectorSpace a, VectorSpace b, VectorSpace c, VectorSpace d, Scalar a ~ s, Scalar b ~ s, Scalar c ~ s, Scalar d ~ s)
  => LinMap s a b -> LinMap s c d -> LinMap s (a :⊗ c) (b :⊗ d)
mapTensor (LinMap m1) (LinMap m2) = LinMap $ ((na*nc) HM.>< (nb*nd)) [(m1 `HM.atIndex` (i, k)) * (m2 `HM.atIndex` (j, l)) | i <- [0..na-1], j <- [0..nc-1], k <- [0..nb-1], l <- [0..nd-1]]
  where
    na = dim (Proxy :: Proxy a)
    nb = dim (Proxy :: Proxy b)
    nc = dim (Proxy :: Proxy c)
    nd = dim (Proxy :: Proxy d)

TensorProd test_Tensor = toFun h t
  where
    f, g :: LinMap Double Double Double
    f = Base.scale 2
    g = Base.scale 3
    h :: LinMap Double (Double :⊗ Double) (Double :⊗ Double)
    h = mapTensor f g

    t :: Double :⊗ Double
    t = TensorProd $ (1 HM.>< 1) [4]

max_stars_repo_name: msakai/essence-of-ad

max_issues_repo_name: msakai/essence-of-ad

max_forks_repo_name: msakai/essence-of-ad