text:
{-# language GADTs #-}
{-# language QuasiQuotes #-}
{-# language ScopedTypeVariables #-}
{-# language FlexibleContexts #-}

module Feldspar.Software.Primitive.Backend where

import Feldspar.Software.Primitive

import Data.Complex (Complex (..))
import Data.Constraint (Dict (..))
import Data.Proxy

-- syntactic.
import Language.Syntactic

-- language-c-quote.
import Language.C.Quote.C
import qualified Language.C.Syntax as C

--imperative-edsl.
import Language.C.Monad
import Language.Embedded.Backend.C

--------------------------------------------------------------------------------
-- * Compilation of software primitives.
--------------------------------------------------------------------------------

viewLitPrim :: ASTF SoftwarePrimDomain a -> Maybe a
viewLitPrim (Sym (Lit a :&: _)) = Just a
viewLitPrim _                   = Nothing

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

instance CompTypeClass SoftwarePrimType
  where
    compType _ (_ :: proxy a) = case softwareRep :: SoftwarePrimTypeRep a of
      BoolST   -> addInclude "<stdbool.h>" >> return [cty| typename bool |]
      Int8ST   -> addInclude "<stdint.h>"  >> return [cty| typename int8_t |]
      Int16ST  -> addInclude "<stdint.h>"  >> return [cty| typename int16_t |]
      Int32ST  -> addInclude "<stdint.h>"  >> return [cty| typename int32_t |]
      Int64ST  -> addInclude "<stdint.h>"  >> return [cty| typename int64_t |]
      Word8ST  -> addInclude "<stdint.h>"  >> return [cty| typename uint8_t |]
      Word16ST -> addInclude "<stdint.h>"  >> return [cty| typename uint16_t |]
      Word32ST -> addInclude "<stdint.h>"  >> return [cty| typename uint32_t |]
      Word64ST -> addInclude "<stdint.h>"  >> return [cty| typename uint64_t |]
      FloatST  -> return [cty| float |]
      DoubleST -> return [cty| double |]
      ComplexFloatST  ->
        addInclude "<tgmath.h>" >> return [cty| float  _Complex |]
      ComplexDoubleST ->
        addInclude "<tgmath.h>" >> return [cty| double _Complex |]

    compLit _ a = case softwarePrimTypeOf a of
      BoolST  ->
        do addInclude "<stdbool.h>"
           return $ if a then [cexp| true |] else [cexp| false |]
      Int8ST   -> return [cexp| $a |]
      Int16ST  -> return [cexp| $a |]
      Int32ST  -> return [cexp| $a |]
      Int64ST  -> return [cexp| $a |]
      Word8ST  -> return [cexp| $a |]
      Word16ST -> return [cexp| $a |]
      Word32ST -> return [cexp| $a |]
      Word64ST -> return [cexp| $a |]
      FloatST  -> return [cexp| $a |]
      DoubleST -> return [cexp| $a |]
      ComplexFloatST  -> return $ compComplexLit a
      ComplexDoubleST -> return $ compComplexLit a

instance CompExp Prim
  where
    compExp = compPrim

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

compUnOp
  :: MonadC m
  => C.UnOp
  -> ASTF SoftwarePrimDomain a
  -> m C.Exp
compUnOp op a = do
  a' <- compPrim $ Prim a
  return $ C.UnOp op a' mempty

compBinOp
  :: MonadC m
  => C.BinOp
  -> ASTF SoftwarePrimDomain a
  -> ASTF SoftwarePrimDomain b
  -> m C.Exp
compBinOp op a b = do
  a' <- compPrim $ Prim a
  b' <- compPrim $ Prim b
  return $ C.BinOp op a' b' mempty

compCast
  :: MonadC m
  => SoftwarePrimTypeRep a
  -> ASTF SoftwarePrimDomain b
  -> m C.Exp
compCast t a = do
  p <- compPrim $ Prim a
  compCastExp t p

compCastExp
  :: MonadC m
  => SoftwarePrimTypeRep a
  -> C.Exp
  -> m C.Exp
compCastExp t p = case softwarePrimWitType t of
  Dict -> do
    typ <- compType (Proxy :: Proxy SoftwarePrimType) t
    return [cexp|($ty:typ) $p|]

compFun
  :: MonadC m
  => String
  -> Args (AST SoftwarePrimDomain) sig
  -> m C.Exp
compFun fun args = do
  as <- sequence $ listArgs (compPrim . Prim) args
  return [cexp| $id:fun($args:as) |]

compRotateL_def = [cedecl|
unsigned int feld_rotl(const unsigned int value, int shift) {
    if ((shift &= sizeof(value)*8 - 1) == 0)
      return value;
    return (value << shift) | (value >> (sizeof(value)*8 - shift));
} |]

compRotateR_def = [cedecl|
unsigned int feld_rotr(const unsigned int value, int shift) {
    if ((shift &= sizeof(value)*8 - 1) == 0)
      return value;
    return (value >> shift) | (value << (sizeof(value)*8 - shift));
} |]

compComplexLit :: (Eq a, Num a, ToExp a) => Complex a -> C.Exp
compComplexLit (r :+ 0) = [cexp| $r |]
compComplexLit (0 :+ i) = [cexp| $i * I |]
compComplexLit (r :+ i) = [cexp| $r + $i * I |]

compAbs
  :: MonadC m
  => SoftwarePrimTypeRep a
  -> ASTF SoftwarePrimDomain a
  -> m C.Exp
compAbs t a
 | boolType t    = error "compAbs: type BoolT not supported"
 | integerType t = addInclude "<stdlib.h>" >> compFun "abs" (a :* Nil)
 | wordType t    = compPrim $ Prim a
 | otherwise     = addInclude "<tgmath.h>" >> compFun "fabs" (a :* Nil)

compSign
  :: MonadC m
  => SoftwarePrimTypeRep a
  -> ASTF SoftwarePrimDomain a
  -> m C.Exp
compSign t a | boolType t = do
  error "compSign: type BoolST not supported"
compSign t a | integerType t = do
  addTagMacro
  a' <- compPrim $ Prim a
  return [cexp| TAG("signum", ($a' > 0) - ($a' < 0)) |]
compSign t a | wordType t = do
  addTagMacro
  a' <- compPrim $ Prim a
  return [cexp| TAG("signum", $a' > 0) |]
compSign FloatST a = do
  addTagMacro
  a' <- compPrim $ Prim a
  return [cexp| TAG("signum", (float) (($a' > 0) - ($a' < 0))) |]
compSign DoubleST a = do
  addTagMacro
  a' <- compPrim $ Prim a
  return [cexp| TAG("signum", (double) (($a' > 0) - ($a' < 0))) |]
compSign ComplexFloatST a = do
  addInclude "<complex.h>"
  addGlobal complexSignf_def
  a' <- compPrim $ Prim a
  return [cexp| feld_complexSignf($a') |]
compSign ComplexDoubleST a = do
  addInclude "<tgmath.h>"
  addGlobal complexSign_def
  a' <- compPrim $ Prim a
  return [cexp| feld_complexSign($a') |]
  -- todo: The floating point cases give `sign (-0.0) = 0.0`, which is (slightly)
  -- wrong. They should return -0.0. I don't know whether it's correct for other
  -- strange values.

complexSignf_def = [cedecl|
float _Complex feld_complexSignf(float _Complex c) {
    float z = cabsf(c);
    if (z == 0)
      return 0;
    else
      return (crealf(c)/z + I*(cimagf(c)/z));
} |]

complexSign_def = [cedecl|
double _Complex feld_complexSign(double _Complex c) {
    double z = cabs(c);
    if (z == 0)
      return 0;
    else
      return (creal(c)/z + I*(cimag(c)/z));
} |]

compDiv
  :: MonadC m
  => SoftwarePrimTypeRep a
  -> ASTF SoftwarePrimDomain a
  -> ASTF SoftwarePrimDomain b
  -> m C.Exp
compDiv Int64ST a b = do
  addGlobal ldiv_def
  compFun "feld_ldiv" (a :* b :* Nil)
compDiv t a b | integerType t = do
  addGlobal div_def
  compFun "feld_div" (a :* b :* Nil)
compDiv t a b | wordType t = compBinOp C.Div a b
compDiv t a b = error $ "compDiv: type " ++ show t ++ " not supported"

ldiv_def = [cedecl|
long int feld_ldiv(long int x, long int y) {
    int q = x/y;
    int r = x%y;
    if ((r!=0) && ((r<0) != (y<0))) --q;
    return q;
} |]

div_def = [cedecl|
int feld_div(int x, int y) {
    int q = x/y;
    int r = x%y;
    if ((r!=0) && ((r<0) != (y<0))) --q;
    return q;
} |]

compMod
  :: MonadC m
  => SoftwarePrimTypeRep a
  -> ASTF SoftwarePrimDomain a
  -> ASTF SoftwarePrimDomain b
  -> m C.Exp
compMod Int64ST a b = do
  addGlobal lmod_def
  compFun "feld_lmod" (a :* b :* Nil)
compMod t a b | integerType t = do
  addGlobal div_def
  compFun "feld_mod" (a :* b :* Nil)
compMod t a b | wordType t = compBinOp C.Mod a b
compMod t a b = error $ "compMod: type " ++ show t ++ " not supported"

lmod_def = [cedecl|
long int feld_lmod(long int x, long int y) {
    int r = x%y;
    if ((r!=0) && ((r<0) != (y<0))) { r += y; }
    return r;
} |]

mod_def = [cedecl|
int feld_mod(int x, int y) {
    int r = x%y;
    if ((r!=0) && ((r<0) != (y<0))) { r += y; }
    return r;
} |]

compRound :: (SoftwarePrimType a, Num a, RealFrac b, MonadC m)
  => SoftwarePrimTypeRep a
  -> ASTF SoftwarePrimDomain b
  -> m C.Exp
compRound t a | integerType t || wordType t = do
  addInclude "<tgmath.h>"
  p <- compFun "lround" (a :* Nil)
  compCastExp t p
compRound t a | floatingType t || complexType t = do
  addInclude "<tgmath.h>"
  p <- compFun "round"  (a :* Nil)
  compCastExp t p
compRound t a = do
  error $ "compSign: type " ++ show t ++ " not supported"

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

compPrim :: MonadC m => Prim a -> m C.Exp
compPrim = simpleMatch (\(s :&: t) -> go t s) . unPrim
  where
    go :: forall m sig
       .  MonadC m
       => SoftwarePrimTypeRep (DenResult sig)
       -> SoftwarePrimConstructs sig
       -> Args (AST SoftwarePrimDomain) sig
       -> m C.Exp
    go _ (FreeVar v) Nil = touchVar v >> return [cexp| $id:v |]
    go t (Lit a)     Nil | Dict <- softwarePrimWitType t
                         = compLit (Proxy :: Proxy SoftwarePrimType) a

    go _ Neg  (a :* Nil)      = compUnOp  C.Negate a
    go _ Add  (a :* b :* Nil) = compBinOp C.Add a b
    go _ Sub  (a :* b :* Nil) = compBinOp C.Sub a b
    go _ Mul  (a :* b :* Nil) = compBinOp C.Mul a b
    go t Div  (a :* b :* Nil) = compDiv t a b
    go _ Quot (a :* b :* Nil) = compBinOp C.Div a b
    go _ Rem  (a :* b :* Nil) = compBinOp C.Mod a b
    go t Mod  (a :* b :* Nil) = compMod t a b
    go t Abs  (a :* Nil)      = compAbs t a
    go t Sign (a :* Nil)      = compSign t a
    go _ FDiv (a :* b :* Nil) = compBinOp C.Div a b

    go _ Exp   args = addInclude "<tgmath.h>" >> compFun "exp" args
    go _ Log   args = addInclude "<tgmath.h>" >> compFun "log" args
    go _ Sqrt  args = addInclude "<tgmath.h>" >> compFun "sqrt" args
    go _ Pow   args = addInclude "<tgmath.h>" >> compFun "pow" args

    go _ Sin   args = addInclude "<tgmath.h>" >> compFun "sin" args
    go _ Cos   args = addInclude "<tgmath.h>" >> compFun "cos" args
    go _ Tan   args = addInclude "<tgmath.h>" >> compFun "tan" args
    go _ Asin  args = addInclude "<tgmath.h>" >> compFun "asin" args
    go _ Acos  args = addInclude "<tgmath.h>" >> compFun "acos" args
    go _ Atan  args = addInclude "<tgmath.h>" >> compFun "atan" args
    go _ Sinh  args = addInclude "<tgmath.h>" >> compFun "sinh" args
    go _ Cosh  args = addInclude "<tgmath.h>" >> compFun "cosh" args
    go _ Tanh  args = addInclude "<tgmath.h>" >> compFun "tanh" args
    go _ Asinh args = addInclude "<tgmath.h>" >> compFun "asinh" args
    go _ Acosh args = addInclude "<tgmath.h>" >> compFun "acosh" args
    go _ Atanh args = addInclude "<tgmath.h>" >> compFun "atanh" args

    go t I2N   (a :* Nil) = compCast t a
    go t I2B   (a :* Nil) = compCast t a
    go t B2I   (a :* Nil) = compCast t a
    go t Round (a :* Nil) = compRound t a

    go _ Complex (a :* b :* Nil) = do
        addInclude "<tgmath.h>"
        a' <- compPrim $ Prim a
        b' <- compPrim $ Prim b
        return $ case (viewLitPrim a, viewLitPrim b) of
          (Just 0, _) -> [cexp| I*$b'       |]
          (_, Just 0) -> [cexp| $a'         |]
          _           -> [cexp| $a' + I*$b' |]

    go _ Polar (m :* p :* Nil)
        | Just 0 <- viewLitPrim m = do
            return [cexp| 0 |]
        | Just 0 <- viewLitPrim p = do
            m' <- compPrim $ Prim m
            return [cexp| $m' |]
        | Just 1 <- viewLitPrim m = do
            p' <- compPrim $ Prim p
            return [cexp| exp(I*$p') |]
        | otherwise = do
            m' <- compPrim $ Prim m
            p' <- compPrim $ Prim p
            return [cexp| $m' * exp(I*$p') |]

    go _ Real      args = addInclude "<tgmath.h>" >> compFun "creal" args
    go _ Imag      args = addInclude "<tgmath.h>" >> compFun "cimag" args
    go _ Magnitude args = addInclude "<tgmath.h>" >> compFun "cabs"  args
    go _ Phase     args = addInclude "<tgmath.h>" >> compFun "carg"  args
    go _ Conjugate args = addInclude "<tgmath.h>" >> compFun "conj"  args

    go _ Not (a :* Nil)      = compUnOp  C.Lnot a
    go _ And (a :* b :* Nil) = compBinOp C.Land a b
    go _ Or  (a :* b :* Nil) = compBinOp C.Lor  a b
    go _ Eq  (a :* b :* Nil) = compBinOp C.Eq   a b
    go _ Neq (a :* b :* Nil) = compBinOp C.Ne   a b
    go _ Lt  (a :* b :* Nil) = compBinOp C.Lt   a b
    go _ Lte (a :* b :* Nil) = compBinOp C.Le   a b
    go _ Gt  (a :* b :* Nil) = compBinOp C.Gt   a b
    go _ Gte (a :* b :* Nil) = compBinOp C.Ge   a b

    go _ Pi Nil = addGlobal pi_def >> return [cexp| FELD_PI |]
      where pi_def = [cedecl|$esc:("#define FELD_PI 3.141592653589793")|]

    go _ BitAnd   (a :* b :* Nil) = compBinOp C.And a b
    go _ BitOr    (a :* b :* Nil) = compBinOp C.Or  a b
    go _ BitXor   (a :* b :* Nil) = compBinOp C.Xor a b
    go _ BitCompl (a :* Nil)      = compUnOp  C.Not a
    go _ ShiftL   (a :* b :* Nil) = compBinOp C.Lsh a b
    go _ ShiftR   (a :* b :* Nil) = compBinOp C.Rsh a b

    go _ RotateL  (a :* b :* Nil) = do
      addGlobal compRotateL_def
      a' <- compPrim $ Prim a
      b' <- compPrim $ Prim b
      return [cexp| feld_rotl($a', $b') |]
    go _ RotateR  (a :* b :* Nil) = do
      addGlobal compRotateR_def
      a' <- compPrim $ Prim a
      b' <- compPrim $ Prim b
      return [cexp| feld_rotr($a', $b') |]

    go _ (ArrIx arr) (i :* Nil) = do
      i' <- compPrim $ Prim i
      touchVar arr
      return [cexp| $id:arr[$i'] |]

    go _ Cond (c :* t :* f :* Nil) = do
      c' <- compPrim $ Prim c
      t' <- compPrim $ Prim t
      f' <- compPrim $ Prim f
      return $ C.Cond c' t' f' mempty

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

addTagMacro :: MonadC m => m ()
addTagMacro = addGlobal [cedecl| $esc:("#define TAG(tag,exp) (exp)") |]

boolType :: SoftwarePrimTypeRep a -> Bool
boolType BoolST     = True
boolType _          = False

integerType :: SoftwarePrimTypeRep a -> Bool
integerType Int8ST  = True
integerType Int16ST = True
integerType Int32ST = True
integerType Int64ST = True
integerType _       = False

wordType :: SoftwarePrimTypeRep a -> Bool
wordType Word8ST    = True
wordType Word16ST   = True
wordType Word32ST   = True
wordType Word64ST   = True
wordType _          = False

floatingType :: SoftwarePrimTypeRep a -> Bool
floatingType FloatST  = True
floatingType DoubleST = True
floatingType _        = False

complexType :: SoftwarePrimTypeRep a -> Bool
complexType ComplexFloatST  = True
complexType ComplexDoubleST = True
complexType _               = False

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


text:
module Main where
import System.Environment (getArgs)
import Control.Concurrent.Async
import Control.Concurrent hiding (yield)
import Control.Monad
import System.Posix.Signals
import Pipes
import Pipes.Concurrent as PC
import qualified Pipes.Group as PG
import qualified Pipes.Prelude as P
import Sound.Pulse.Simple
import Data.WAVE
import qualified Data.List.Extra as E
import qualified Numeric.FFT as F
import Data.Complex
import Data.Monoid
import Data.Int
import Data.Maybe
import qualified UI.HSCurses.CursesHelper as CH
import qualified UI.HSCurses.Curses as C

takeStep :: Int -> [a] -> [a]
takeStep _ [] = []
takeStep n (x:xs)
  | n >= 0 = x : takeStep n (drop (n-1) xs)
  | otherwise = takeStep (-n) (reverse xs)

rjust :: Int -> a -> [a] -> [a]
rjust size filler inp
  | length inp >= size = inp
  | otherwise = inp ++ replicate (size - length inp) filler

fft :: Int -> Pipe [WAVESample] [Double] IO()
fft chunkSize = forever $ do
    let step = 18
    let fftNum = 128
    a <- await
    let vec = rjust fftNum 0 $takeStep step $ takeStep 2 a
    let spec = map magnitude $ F.dft(map fromIntegral $ vec)
    let ratio = 3e-9
    yield $ map (\x -> x * ratio) spec

convSpectrumMat :: [Double] -> Int -> Int -> [(Int, String)]
convSpectrumMat input offsetY height =
    serializer 0
    where
        arr = rjust 1 0 $ map (\x -> floor x) input
        serializer y
            | height - offsetY <= y= []
            | otherwise  = [(abs(height - y), map (\i -> if i>y then '#' else ' ') arr)]++serializer (y+1)

spectrogram :: Int -> Int -> String-> Consumer [Double] IO ()
spectrogram chunkNo chunkNum filename = forever $ do
    (height, width) <-lift $ C.scrSize
    a <- await
    let arr = rjust 1 0 a
    let offsetX = 0
    let offsetY = 4
    let mat = convSpectrumMat arr offsetY (height-2)
    let cprintLn y str = C.mvWAddStr C.stdScr y offsetX str
    lift $ mapM (\(y, s) -> cprintLn y $take width s) mat
    lift $ cprintLn 0 ("     chunkNo:"++(show chunkNo) ++ "/" ++ show chunkNum)
    lift $ cprintLn 1 ("max spectrum:"++(show $ maximum arr ))
    lift $ cprintLn 2 ("num spectrum:"++(show $ length arr ))
    lift $ cprintLn 3 ("    filename:"++filename)
    lift $ C.refresh
    spectrogram (chunkNo+1) chunkNum filename

concat' :: [[a]] -> [a]
concat' xss = foldr (++) [] xss

main = do
    tid <- myThreadId
    let handler = do
            putStrLn "keybord interupted"
            killThread tid
    installHandler keyboardSignal (Catch handler) Nothing
    (waveFileName:_) <- getArgs
    input <- getWAVEFile waveFileName
    let header = waveHeader input
    let resolution = 2^64 :: Double
    let channels = waveNumChannels header
    let samples = waveSamples input
    let sampleLength = fromMaybe 0 $waveFrames header
    let chunkSize = 2^11
    let chunkNum = (sampleLength `div` chunkSize) + 1

    s<-simpleNew Nothing "ninja" Play Nothing "youjo"
        (SampleSpec (S32 LittleEndian) (waveFrameRate header) 2) Nothing Nothing
        -- WAVESample is 32 bit, so you must set S32 in pulseaudio or casting.

    --runEffect $ each (E.chunksOf 2048 samples) >-> P.map concat' >-> P.mapM_ (simpleWrite s)

    (plotterOutput, plotterInput) <- spawn unbounded
    player <- async $ do
        runEffect $ each (((E.chunksOf (chunkSize)).(map$ map (`div` 3))) samples) >-> P.map concat' >-> P.chain(simpleWrite s) >-> toOutput plotterOutput
        performGC
    plotter <- async $ do
        CH.start
        C.refresh
        runEffect $ fromInput plotterInput >-> fft chunkSize>-> spectrogram 0 chunkNum waveFileName
        C.endWin
        performGC
    mapM_ wait [player, plotter]

    simpleDrain s
    simpleFree s

ext: hs

max_stars_repo_name: cocuh/hspectrogram

max_issues_repo_name: cocuh/hspectrogram

max_forks_repo_name: cocuh/hspectrogram


text:
-- |
-- Module      : Occlusion.Window
-- Description :
-- Copyright   : (c) Jonatan H Sundqvist, 2015
-- License     : MIT
-- Maintainer  : Jonatan H Sundqvist
-- Stability   : experimental|stable
-- Portability : POSIX (not sure)
--

-- Created September 20 2015

-- TODO | -
--        -

-- SPEC | -
--        -

--------------------------------------------------------------------------------------------------------------------------------------------
-- GHC Pragmas
--------------------------------------------------------------------------------------------------------------------------------------------

--------------------------------------------------------------------------------------------------------------------------------------------
-- API
--------------------------------------------------------------------------------------------------------------------------------------------
module Occlusion.Window where

--------------------------------------------------------------------------------------------------------------------------------------------
-- We'll need these
--------------------------------------------------------------------------------------------------------------------------------------------
import Data.Complex

import Graphics.UI.Gtk
import qualified Graphics.Rendering.Cairo as Cairo

import Occlusion.Types
import qualified Occlusion.Core   as Core
import qualified Occlusion.Events as Events

--------------------------------------------------------------------------------------------------------------------------------------------
-- Functions
--------------------------------------------------------------------------------------------------------------------------------------------
-- |
create :: Complex Double -> IO (Window, DrawingArea)
create (winx:+winy) = do
  initGUI
  window <- windowNew
  frame  <- frameNew

  set window [windowTitle := "Occlusion"]

  canvas <- drawingAreaNew
  containerAdd frame canvas
  set window [ containerChild := frame ]
  windowSetDefaultSize window (round winx) (round winy)
  windowSetIconFromFile window "C:/Users/Jonatan/Desktop/Haskell/projects/Occlusion/assets/images/TreeAtDusk.png"

  widgetAddEvents canvas [PointerMotionMask, ButtonMotionMask, ButtonPressMask] -- MouseButton1Mask
  widgetShowAll window

  return (window, canvas)

max_stars_repo_name: SwiftsNamesake/Occlusion

max_issues_repo_name: SwiftsNamesake/Occlusion

max_forks_repo_name: SwiftsNamesake/Occlusion


text:
{-# LANGUAGE CPP #-}
--------------------------------------------------------------------------------
{- |
Module      :  Numeric.LinearAlgebra.HMatrix
Copyright   :  (c) Alberto Ruiz 2006-14
License     :  BSD3
Maintainer  :  Alberto Ruiz
Stability   :  provisional

compatibility with previous version, to be removed

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

module Numeric.LinearAlgebra.HMatrix (
    module Numeric.LinearAlgebra,
    (¦),(——),ℝ,ℂ,(<·>),app,mul, cholSH, mbCholSH, eigSH', eigenvaluesSH', geigSH'
) where

import Numeric.LinearAlgebra
import Internal.Util
import Internal.Algorithms(cholSH, mbCholSH, eigSH', eigenvaluesSH', geigSH')
#if MIN_VERSION_base(4,11,0)
import Prelude hiding ((<>))
#endif

infixr 8 <·>
(<·>) :: Numeric t => Vector t -> Vector t -> t
(<·>) = dot

app :: Numeric t => Matrix t -> Vector t -> Vector t
app m v = m #> v

mul :: Numeric t => Matrix t -> Matrix t -> Matrix t
mul a b = a <> b

max_stars_repo_name: schnecki/hmatrix-float

max_issues_repo_name: schnecki/hmatrix-float

max_forks_repo_name: schnecki/hmatrix-float


text:
{-# LANGUAGE NoImplicitPrelude #-}
import Protolude
import Data.Vector as V
import Numeric.LinearAlgebra as H

-- |
-- >>> let x = V.fromList [0..5]
-- >>> V.length x
-- 6
--
-- >>> V.null x
-- False

-- | hmatrix examples
-- >>> let x = (2><2) [0..3] :: Matrix Double
-- >>> x
-- (2><2)
--  [ 0.0, 1.0
--  , 2.0, 3.0 ]
--
-- Matrix multiplication
-- >>> let y = fromLists [[0, 1], [2, 3]] :: Matrix Double
-- >>> x H.<> y
-- (2><2)
--  [ 2.0,  3.0
--  , 6.0, 11.0 ]

main :: IO ()
main = do
    -- vector examples
    let x = V.fromList [0..5]

    print $ V.length x -- 6
    print $ V.null x -- False

    -- indexing
    print $ x V.! 1 -- 1
    print $ V.head x -- 0
    print $ last x -- 5

    -- Slicing
    print $ slice 2 3 x -- [2, 3, 4]
    print $ V.splitAt 2 x -- ([0, 1], [2, 3, 4, 5])

    -- Prepending and Appending
    print $ cons (-1) x -- [-1, 0, 1, 2, 3, 4, 5]
    print $ snoc x 6 -- [0, 1, 2, 3, 4, 5, 6]

    -- Concatenation
    print $ x V.++ x -- [0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5]
    print $ V.concat [x, x] -- [0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5]

    -- Update
    print $ x // [(0, 1), (2, 6)] -- [1, 1, 6, 3, 4, 5]
    print $ V.map (+2) x -- [2, 3, 4, 5, 6, 7]

    -- hmatrix examples
    -- Creating matrices
    let x = (2><2) [0..3] :: Matrix Double

    -- [ 0.0, 1.0
    -- , 2.0, 3.0 ]
    let y = fromLists [[0, 1], [2, 3]] :: Matrix Double

    -- Random matrices
    r <- randn 2 3
    -- [ 0.7764496757867578,    1.246311658930589,    -0.684233085372455
    -- , -2.540045307941425, -0.20975584071908912, -9.039537343065803e-3 ]

    -- Matrix multiplication
    print $ x H.<> y -- [ 2.0,  3.0
                   -- , 6.0, 11.0 ]

    -- Transpose
    print $ tr x -- [ 0.0, 2.0
                 -- , 1.0, 3.0 ]

    -- Matrix slicing
    print $ r ?? (H.All, Take 2) -- [ 0.7764496757867578,    1.246311658930589
                               -- , -2.540045307941425, -0.20975584071908912 ]

    -- Mapping over matrices
    print $ cmap ((+ 2) . (*2)) x -- [ 2.0, 4.0
                                  -- , 6.0, 8.0 ]

    -- Flatten
    print $ flatten x -- [0.0, 1.0, 2.0, 3.0]

max_stars_repo_name: DataHaskell/data-haskell-examples

max_issues_repo_name: DataHaskell/data-haskell-examples

max_forks_repo_name: DataHaskell/data-haskell-examples


text:
module Mbrot where
import           Codec.Picture.Types
import           Coloring
import           Control.Arrow       ((***))
import           Data.Complex
import           ImageOutput

-- |Size parameters
size :: (Int,Int)
size = (1400,800)

-- |Floating point size, calculated
sizef :: (Float,Float)
sizef = (fromIntegral *** fromIntegral) size

-- |Maximum iterations
maxIter :: Int
maxIter = 1000

-- |Color step for integer esacpe coloring
mult :: Float
mult = 255 / fromIntegral maxIter

-- |Color step for continuously colored
multg :: Float
multg = 400 / fromIntegral maxIter

-- |Function for rendering pixel in greyscale
drawGreyscale :: Bool    -- ^Continuously colored
              -> Int
              -> Int
              -> Pixel8
drawGreyscale c x y
    |c =  fromIntegral (round(float * mult))
    |otherwise = fromIntegral(int * round mult)
    where int = mbrot maxIter (((x1/x2)-0.5) :+ (y1/y2))
          float = mbrotContinuous maxIter (((x1/x2)-0.5) :+ (y1/y2))
          x1 = 2 * (fromIntegral x - x2)
          y1 = 2 * (fromIntegral y - y2)
          x2 = fst sizef /2
          y2 = snd sizef /2

-- |Deprecated. Use 'drawColorZoom' instead
drawColor :: Int -> Int -> PixelRGB8
drawColor x y =  colorToPixel $ color (fromIntegral (int*400)/1000)
    where int = mbrot maxIter ((x1/x2) :+ (y1/y2))
          x1 = 1.75 * (fromIntegral x - x2)
          y1 = fromIntegral y - y2
          x2 = fst sizef /2
          y2 = snd sizef /2

-- |Function for rendering pixel in color, with zoom level and focus
drawColorZoom :: Float         -- ^ Zoom level
              -> (Float,Float) -- ^ Focus point (x,y)
              -> Int
              -> Int
              -> PixelRGB8
drawColorZoom zoom (fx, fy) x y =  colorToPixel $ color (float*multg)
    where float = mbrotContinuous maxIter (((x1/x2)+fx) :+ ((y1/y2)+fy))
          x1 = 1.75 * (fromIntegral x - x2) / zoom
          y1 = (fromIntegral y - y2) / zoom
          x2 = fst sizef /2
          y2 = snd sizef /2

-- |Function for drawing color gradient. Tests for correct linear interpolation of color
drawTestGradient :: Int -> Int -> PixelRGB8
drawTestGradient x _ = colorToPixel $ color (fromIntegral x *400/1000)

-- |Function for integer escape time of mbrot set of Re(x) + Im(y)
mbrot :: RealFloat x
     => Int           -- ^ Maximum iterations
     -> Complex x     -- ^ Re(x) + Im(y)
     -> Int
mbrot n = mbrot' maxIter (0 :+ 0)
    where
        mbrot' :: RealFloat x => Int -> Complex x -> Complex x -> Int
        mbrot' n a c
            |magnitude newIter < 65535 && n > 0 = mbrot' (n-1) newIter c
            |otherwise = maxIter - n
                where
                    newIter = a*a + c

-- |Function for real floating escape time of mbrot set of Re(x) + Im(y)
mbrotContinuous :: RealFloat x
                => Int           -- ^ Maximum iterations
                -> Complex x     -- ^ Re(x) + Im(y)
                -> x
mbrotContinuous n = mbrot' maxIter (0 :+ 0)
    where
        mbrot' :: RealFloat x => Int -> Complex x -> Complex x -> x
        mbrot' n a c
            |magnitude newIter < 65535 && n > 0 = mbrot' (n-1) newIter c
            |otherwise = normalized
                where
                    newIter = a*a + c
                    normalized =  if n > 0
                                  then
                                    fromIntegral(maxIter - n +1) - logBase 2 ( logBase 2 (magnitude newIter ^ 2) / 2 )
                                    --fromIntegral (maxIter - n +1) - (log (log (magnitude newIter)))/ log (2.0)
                                  else fromIntegral maxIter

max_stars_repo_name: veniversum/fractal-haskell

max_issues_repo_name: veniversum/fractal-haskell

max_forks_repo_name: veniversum/fractal-haskell


text:
{-# LANGUAGE NegativeLiterals #-}

module Main (main) where

import Criterion.Types
import qualified Data.Aeson as Aeson
import qualified Data.Vector as V
import Properties
import Statistics.Resampling.Bootstrap (Estimate(..))
import Test.Framework (defaultMain)
import Test.Framework.Providers.HUnit
import Test.HUnit

r1 :: Report
r1 = Report 0 "" [] v1 s1 (Outliers 0 0 0 0 0) []
 where
  m1 = Measured 4.613000783137977e-05 3.500000000000378e-05 31432 1 0 0 0 0.0 0.0 0.0 0.0
  v1 = V.fromList [m1]
  est1 = Estimate 0.0 0.0 0.0 0.0
  s1 = SampleAnalysis [] 0.0 est1 est1 (OutlierVariance Unaffected "" 0.0)

m2 :: Measured
m2 = Measured {measTime = 1.1438998626545072e-5
              , measCpuTime = 1.2000000001677336e-5
              , measCycles = 6208
              , measIters = 1

              , measAllocated = -9223372036854775808
              , measNumGcs = -9223372036854775808
              , measBytesCopied = -9223372036854775808

              , measMutatorWallSeconds = -1/0
              , measMutatorCpuSeconds = -1/0
              , measGcWallSeconds = -1/0
              , measGcCpuSeconds = -1/0}

main :: IO ()
main = defaultMain
       [ Properties.tests
       , testCase "json-roundtrip1"
           (assertEqual "round trip simple Measured"
              (Right m2) (Aeson.eitherDecode (Aeson.encode m2)))
       , testCase "json-roundtrip2"
           (assertEqual "round trip simple Report"
              (Right r1) (Aeson.eitherDecode (Aeson.encode r1)))
       ]

max_stars_repo_name: Shimuuar/criterion

max_issues_repo_name: Shimuuar/criterion

max_forks_repo_name: Shimuuar/criterion


text:
module STCR2Z1T0Image where

import           Control.Monad                as M
import           Control.Monad.Parallel       as MP
import           Data.Array.Repa              as R
import           Data.Binary                  (decodeFile)
import           Data.Complex
import           Data.List                    as L
import           Data.Vector.Generic          as VG
import           Data.Vector.Storable         as VS
import           DFT.Plan
import           Filter.Gaussian
import           Filter.Pinwheel
import           FokkerPlanck.DomainChange
import           FokkerPlanck.MonteCarlo
import           FokkerPlanck.Pinwheel
import           Image.IO
import           Image.Transform
import           STC
import           STC.OrientationScaleAnalysis
import           System.Directory
import           System.Environment
import           System.FilePath
import           System.Random
import           Types
import           Utils.Array

main = do
  args@(numPointStr:numOrientationStr:sigmaStr:taoStr:lenStr:initialScaleStr:numTrailStr:maxTrailStr:theta0FreqsStr:thetaFreqsStr:histFilePath:alphaStr:pinwheelFlagStr:imagePath:numIterationStr:writeSourceFlagStr:numThreadStr:_) <-
    getArgs
  print args
  let numPoint = read numPointStr :: Int
      numOrientation = read numOrientationStr :: Int
      sigma = read sigmaStr :: Double
      tao = read taoStr :: Double
      len = read lenStr :: Int
      initialScale = read initialScaleStr :: Double
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      theta0Freq = read theta0FreqsStr :: Double
      theta0Freqs = [-theta0Freq .. theta0Freq]
      thetaFreq = read thetaFreqsStr :: Double
      thetaFreqs = [-thetaFreq .. thetaFreq]
      alpha = read alphaStr :: Double
      pinwheelFlag = read pinwheelFlagStr :: Bool
      numIteration = read numIterationStr :: Int
      writeSourceFlag = read writeSourceFlagStr :: Bool
      numThread = read numThreadStr :: Int
      folderPath = "output/test/STCR2Z1T0Image"
  createDirectoryIfMissing True folderPath
  flag <- doesFileExist histFilePath
  (ImageRepa _ img'') <- readImageRepa imagePath False
  let img'
        -- resize25D (numPoint, numPoint) (0, 255) .
        -- R.backpermute (Z :. nf :. numPoint :. numPoint) id .
        -- upsample [4, 4, 1] . downsample [4, 4, 1] $
       = img''
      (Z :. nf :. _ :. _) = extent img''
      (Z :. _ :. cols :. rows) = extent img'
      pinwheelParams = PinwheelParams rows cols alpha (exp 1) theta0Freqs [0]
  radialArr <-
    if flag
      then R.map magnitude . getNormalizedHistogramArr <$>
           decodeFile histFilePath
      else do
        putStrLn "Couldn't find a Green's function data. Start simulation..."
        solveMonteCarloR2Z1T0Radial
          numThread
          numTrail
          maxTrail
          numPoint
          numPoint
          sigma
          tao
          1
          theta0Freqs
          thetaFreqs
          histFilePath
          (emptyHistogram
             [ (round . sqrt . fromIntegral $ 2 * (div numPoint 2) ^ 2)
             , L.length theta0Freqs
             , L.length thetaFreqs
             ]
             0)
  arrR2Z1T0 <-
    computeUnboxedP $
    computeR2Z1T0ArrayRadial
      radialArr
      numPoint
      numPoint
      1.5
      thetaFreqs
      theta0Freqs
  plan <- makeR2Z1T0Plan emptyPlan arrR2Z1T0
  -- (plan1, imgF) <- makeImagePlan plan0 . computeS $ img'
  -- (plan2, filterF, filterPIF) <- pinwheelFilter plan1 pinwheelParams
  -- (plan, gaussianFilterF) <-
  --   gaussian2DFilter plan2 (Gaussian2DParams 1 rows cols)
  -- -- img <- convolveGaussian2D plan gaussianFilterF . R.map (:+ 0) $ img'
  -- -- imgVec <- dropPixel (2 / 3) . toUnboxed . computeS . filterImage $ img'
  let img = R.map (:+ 0) img' -- = . fromUnboxed (extent img') $ imgVec
  plotImageRepa
    (folderPath </> "input.png") -- (ImageRepa 8 img)
    (ImageRepa 8 . computeS . R.map magnitude $ img)
  -- convolvedImg' <-
  --   convolvePinwheel
  --     plan
  --     filterPIF
  --     (R.slice imgF (Z :. (0 :: Int) :. All :. All))
  -- let convolvedImg =
  --       computeS . R.slice convolvedImg' $
  --       (Z :. All :. (0 :: Int) :. All :. All)
  -- let (initialDistSource, initialDistSink) =
  --       analyzeOrientation numOrientation theta0Freqs convolvedImg
  powerMethod1
    plan
    folderPath
    cols
    rows
    numOrientation
    thetaFreqs
    theta0Freqs
    arrR2Z1T0
    numIteration
    writeSourceFlag
    ""
    0.1
    (R.traverse img (const (Z :. (L.length theta0Freqs) :. cols :. rows)) $ \f (Z :. _ :. i :. j) ->
       if magnitude (f (Z :. (0 :: Int) :. i :. j)) > 0
         then 1
         else 0)
    (R.traverse
       img
       (const
          (Z :. (L.length thetaFreqs) :. (L.length theta0Freqs) :. cols :. rows)) $ \f (Z :. k :. _ :. i :. j) ->
       if k == div (L.length thetaFreqs) 2
         then f (Z :. (0 :: Int) :. i :. j) /
              (fromIntegral $ L.length theta0Freqs)
         else 0)

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-|
Module      : Grenade.Layers.Crop
Description : Cropping layer
Copyright   : (c) Huw Campbell, 2016-2017
License     : BSD2
Stability   : experimental
-}
module Grenade.Layers.Crop (
    Crop (..)
  ) where

import           Data.Maybe
import           Data.Proxy
import           Data.Singletons.TypeLits
import           GHC.TypeLits

import           Grenade.Core
import           Grenade.Layers.Internal.Pad

import           Numeric.LinearAlgebra (konst, subMatrix, diagBlock)
import           Numeric.LinearAlgebra.Static (extract, create)

-- | A cropping layer for a neural network.
data Crop :: Nat
          -> Nat
          -> Nat
          -> Nat -> * where
  Crop :: Crop cropLeft cropTop cropRight cropBottom

instance Show (Crop cropLeft cropTop cropRight cropBottom) where
  show Crop = "Crop"

instance UpdateLayer (Crop l t r b) where
  type Gradient (Crop l t r b) = ()
  runUpdate _ x _ = x
  createRandom = return Crop

-- | A two dimentional image can be cropped.
instance ( KnownNat cropLeft
         , KnownNat cropTop
         , KnownNat cropRight
         , KnownNat cropBottom
         , KnownNat inputRows
         , KnownNat inputColumns
         , KnownNat outputRows
         , KnownNat outputColumns
         , (outputRows + cropTop + cropBottom) ~ inputRows
         , (outputColumns + cropLeft + cropRight) ~ inputColumns
         ) => Layer (Crop cropLeft cropTop cropRight cropBottom) ('D2 inputRows inputColumns) ('D2 outputRows outputColumns) where
  type Tape (Crop cropLeft cropTop cropRight cropBottom) ('D2 inputRows inputColumns) ('D2 outputRows outputColumns) = ()
  runForwards Crop (S2D input) =
    let cropl = fromIntegral $ natVal (Proxy :: Proxy cropLeft)
        cropt = fromIntegral $ natVal (Proxy :: Proxy cropTop)
        nrows = fromIntegral $ natVal (Proxy :: Proxy outputRows)
        ncols = fromIntegral $ natVal (Proxy :: Proxy outputColumns)
        m  = extract input
        r  = subMatrix (cropt, cropl) (nrows, ncols) m
    in  ((), S2D . fromJust . create $ r)
  runBackwards _ _ (S2D dEdy) =
    let cropl = fromIntegral $ natVal (Proxy :: Proxy cropLeft)
        cropt = fromIntegral $ natVal (Proxy :: Proxy cropTop)
        cropr = fromIntegral $ natVal (Proxy :: Proxy cropRight)
        cropb = fromIntegral $ natVal (Proxy :: Proxy cropBottom)
        eo    = extract dEdy
        vs    = diagBlock [konst 0 (cropt,cropl), eo, konst 0 (cropb,cropr)]
    in  ((), S2D . fromJust . create $ vs)

-- | A two dimentional image can be cropped.
instance ( KnownNat cropLeft
         , KnownNat cropTop
         , KnownNat cropRight
         , KnownNat cropBottom
         , KnownNat inputRows
         , KnownNat inputColumns
         , KnownNat outputRows
         , KnownNat outputColumns
         , KnownNat channels
         , KnownNat (inputRows * channels)
         , KnownNat (outputRows * channels)
         , (outputRows + cropTop + cropBottom) ~ inputRows
         , (outputColumns + cropLeft + cropRight) ~ inputColumns
         ) => Layer (Crop cropLeft cropTop cropRight cropBottom) ('D3 inputRows inputColumns channels) ('D3 outputRows outputColumns channels) where
  type Tape (Crop cropLeft cropTop cropRight cropBottom) ('D3 inputRows inputColumns channels) ('D3 outputRows outputColumns channels)  = ()
  runForwards Crop (S3D input) =
    let padl  = fromIntegral $ natVal (Proxy :: Proxy cropLeft)
        padt  = fromIntegral $ natVal (Proxy :: Proxy cropTop)
        padr  = fromIntegral $ natVal (Proxy :: Proxy cropRight)
        padb  = fromIntegral $ natVal (Proxy :: Proxy cropBottom)
        inr   = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        inc   = fromIntegral $ natVal (Proxy :: Proxy inputColumns)
        outr  = fromIntegral $ natVal (Proxy :: Proxy outputRows)
        outc  = fromIntegral $ natVal (Proxy :: Proxy outputColumns)
        ch    = fromIntegral $ natVal (Proxy :: Proxy channels)
        m     = extract input
        cropped = crop ch padl padt padr padb outr outc inr inc m
    in  ((), S3D . fromJust . create $ cropped)

  runBackwards Crop () (S3D gradient) =
    let padl  = fromIntegral $ natVal (Proxy :: Proxy cropLeft)
        padt  = fromIntegral $ natVal (Proxy :: Proxy cropTop)
        padr  = fromIntegral $ natVal (Proxy :: Proxy cropRight)
        padb  = fromIntegral $ natVal (Proxy :: Proxy cropBottom)
        inr   = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        inc   = fromIntegral $ natVal (Proxy :: Proxy inputColumns)
        outr  = fromIntegral $ natVal (Proxy :: Proxy outputRows)
        outc  = fromIntegral $ natVal (Proxy :: Proxy outputColumns)
        ch    = fromIntegral $ natVal (Proxy :: Proxy channels)
        m     = extract gradient
        padded = pad ch padl padt padr padb outr outc inr inc m
    in  ((), S3D . fromJust . create $ padded)


text:
{-# LANGUAGE NamedFieldPuns, TemplateHaskell #-}

module School.Unit.Test.UnitBackward
( unitBackwardTest) where

import Conduit ((.|), sinkList, yield)
import Control.Monad.IO.Class (liftIO)
import Data.Either (isLeft)
import Numeric.LinearAlgebra (ident)
import School.TestUtils (assertRight, fromRight, randomAffineParams, randomMatrix, testState)
import School.Train.TrainState (TrainState(..), def)
import School.Types.FloatEq ((~=))
import School.Types.PingPong (pingPongSingleton, reversePingPong, toPingPong)
import School.Unit.Affine (affine)
import School.Unit.RecLin (recLin)
import School.Unit.Unit (Unit(..))
import School.Unit.UnitActivation (UnitActivation(..))
import School.Unit.UnitBackward
import School.Unit.UnitGradient (UnitGradient(..), isGradientFail)
import School.Unit.UnitParams (UnitParams(..))
import School.Utils.Tuple (snd3)
import Test.Tasty (TestTree)
import Test.Tasty.QuickCheck
import Test.Tasty.TH
import Test.QuickCheck.Monadic (assert, monadicIO)

prop_affine_input_fail :: Property
prop_affine_input_fail = monadicIO $ do
  let backward = unitBackward affine
  let acts = [ApplyFail "init"]
  let inGrad = BatchGradient $ ident 1
  let network =  yield (acts, inGrad, 0)
              .| backward
              .| sinkList
  result <- testState network def
  assert $ isLeft result

prop_affine_gradient_fail :: Property
prop_affine_gradient_fail = monadicIO $ do
  let backward = unitBackward affine
  let acts = [BatchActivation $ ident 1]
  let inGrad = GradientFail "init"
  let network =  yield (acts, inGrad, 0)
              .| backward
              .| sinkList
  result <- testState network def
  assert $ isLeft result

prop_affine_param_fail :: Positive Int -> Positive Int -> Property
prop_affine_param_fail (Positive fSize) (Positive oSize) = monadicIO $ do
  let backward = unitBackward affine
  actMat <- liftIO $ randomMatrix oSize fSize
  gradMat <- liftIO $ randomMatrix oSize fSize
  let acts = [BatchActivation actMat]
  let fStack = (acts, BatchGradient gradMat, 0)
  let network =  yield fStack
              .| backward
              .| sinkList
  result <- testState network def
  assertRight (isGradientFail . snd3 . head . fst)
              result

prop_reclin_gradient :: Positive Int -> Positive Int -> Property
prop_reclin_gradient (Positive bSize) (Positive fSize) = monadicIO $ do
  let backward = unitBackward recLin
  actMat <- liftIO $ randomMatrix bSize fSize
  let input = BatchActivation actMat
  gradMat <- liftIO $ randomMatrix bSize fSize
  let inGrad = BatchGradient gradMat
  let fStack = ([input], inGrad, 0)
  let network =  yield fStack
              .| backward
              .| sinkList
  result <- testState network def
  let (check, _) = deriv recLin EmptyParams inGrad input
  assertRight (((~=) check) . snd3 . head . fst)
              result

prop_affine_derivs :: Positive Int -> Positive Int -> Positive Int -> Property
prop_affine_derivs (Positive bSize) (Positive fSize) (Positive oSize) = monadicIO $ do
  let backward = unitBackward affine
  actMat <- liftIO $ randomMatrix bSize fSize
  let input = BatchActivation actMat
  gradMat <- liftIO $ randomMatrix bSize oSize
  let inGrad = BatchGradient gradMat
  let fStack = ([input], inGrad, 0)
  let network =  yield fStack
              .| backward
              .| sinkList
  params <- liftIO $ randomAffineParams fSize oSize
  let paramList = pingPongSingleton params
  let initState = def { paramList }
  result <- testState network initState
  let (_, check) = deriv affine params inGrad input
  assertRight (((~=) check) . head . paramDerivs . snd)
              result

prop_deriv_aff_rl_aff_rl :: Positive Int
                         -> Positive Int
                         -> Positive Int
                         -> Positive Int
                         -> Property
prop_deriv_aff_rl_aff_rl (Positive b)
                         (Positive f)
                         (Positive h)
                         (Positive o) = monadicIO $ do
  inMat <- liftIO $ randomMatrix b f
  act1Mat <- liftIO $ randomMatrix b h
  act2Mat <- liftIO $ randomMatrix b h
  act3Mat <- liftIO $ randomMatrix b o
  let acts = BatchActivation <$> [ act3Mat
                                 , act2Mat
                                 , act1Mat
                                 , inMat
                                 ]
  gradMat <- liftIO $ randomMatrix b o
  let inGrad = BatchGradient gradMat
  let network =  yield (acts, inGrad, 0)
              .| unitBackward recLin
              .| unitBackward affine
              .| unitBackward recLin
              .| unitBackward affine
              .| sinkList
  params1 <- liftIO $ randomAffineParams f h
  params2 <- liftIO $ randomAffineParams h o
  let allParams = toPingPong [ params1
                             , EmptyParams
                             , params2
                             , EmptyParams
                             ]
  let paramList = reversePingPong $ fromRight (pingPongSingleton EmptyParams) allParams
  let initState = def { paramList }
  result <- testState network initState
  let (grad1, _) = deriv recLin EmptyParams inGrad (head acts)
  let (grad2, dParams2) = deriv affine params2 grad1 (acts!!1)
  let (grad3, _) = deriv recLin EmptyParams grad2 (acts!!2)
  let (_, dParams1) = deriv affine params1 grad3 (last acts)
  let check = [ dParams1
              , EmptyParams
              , dParams2
              , EmptyParams
              ]
  assertRight (((~=) check) . paramDerivs . snd)
              result

unitBackwardTest :: TestTree
unitBackwardTest = $(testGroupGenerator)

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


text:
{-# LANGUAGE FlexibleInstances #-}
-----------------------------------------------------------------------------
-- |
-- Module     : Data.Elem.LAPACK.C
-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@stanford.edu>
-- Stability  : experimental
--
-- Low-level interface to LAPACK.
--

module Data.Elem.LAPACK.C
    where

import Control.Exception( assert )
import Control.Monad
import Data.Complex( Complex )
import Data.Elem.BLAS.Level3
import Data.Matrix.Class( TransEnum(..), SideEnum(..) )
import Foreign
import LAPACK.CTypes

import Data.Elem.LAPACK.Double
import Data.Elem.LAPACK.Zomplex

-- | The LAPACK typeclass.
class (BLAS3 e) => LAPACK e where
    geqrf :: Int -> Int -> Ptr e -> Int -> Ptr e -> IO ()
    gelqf :: Int -> Int -> Ptr e -> Int -> Ptr e -> IO ()
    unmqr :: SideEnum -> TransEnum -> Int -> Int -> Int -> Ptr e -> Int -> Ptr e -> Ptr e -> Int -> IO ()
    unmlq :: SideEnum -> TransEnum -> Int -> Int -> Int -> Ptr e -> Int -> Ptr e -> Ptr e -> Int -> IO ()
    larfg :: Int -> Ptr e -> Ptr e -> Int -> IO e

callWithWork :: (Storable e) => (Ptr e -> Int -> IO a) -> IO a
callWithWork call =
    alloca $ \pQuery -> do
        call pQuery (-1)
        ldWork <- peek (castPtr pQuery) :: IO Double
        let lWork = max 1 $ ceiling ldWork
        allocaArray lWork $ \pWork -> do
            call pWork lWork

checkInfo :: Int -> IO ()
checkInfo info = assert (info == 0) $ return ()

instance LAPACK Double where
    geqrf m n pA ldA pTau =
        checkInfo =<< callWithWork (dgeqrf m n pA ldA pTau)
    gelqf m n pA ldA pTau =
        checkInfo =<< callWithWork (dgelqf m n pA ldA pTau)
    unmqr s t m n k pA ldA pTau pC ldC =
        checkInfo =<< callWithWork (dormqr (cblasSide s) (cblasTrans t) m n k pA ldA pTau pC ldC)
    unmlq s t m n k pA ldA pTau pC ldC =
        checkInfo =<< callWithWork (dormlq (cblasSide s) (cblasTrans t) m n k pA ldA pTau pC ldC)
    larfg n alpha x incx = with 0 $ \pTau ->
        dlarfg n alpha x incx pTau >> peek pTau

instance LAPACK (Complex Double) where
    geqrf m n pA ldA pTau =
        checkInfo =<< callWithWork (zgeqrf m n pA ldA pTau)
    gelqf m n pA ldA pTau =
        checkInfo =<< callWithWork (zgelqf m n pA ldA pTau)
    unmqr s t m n k pA ldA pTau pC ldC =
        checkInfo =<< callWithWork (zunmqr (cblasSide s) (cblasTrans t) m n k pA ldA pTau pC ldC)
    unmlq s t m n k pA ldA pTau pC ldC =
        checkInfo =<< callWithWork (zunmlq (cblasSide s) (cblasTrans t) m n k pA ldA pTau pC ldC)
    larfg n alpha x incx = with 0 $ \pTau ->
        zlarfg n alpha x incx pTau >> peek pTau

max_stars_repo_name: patperry/lapack

max_issues_repo_name: patperry/lapack

max_forks_repo_name: patperry/lapack


text:
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE PatternSynonyms #-}
{-# LANGUAGE ViewPatterns #-}
{-# LANGUAGE OverloadedStrings #-}

{-|
Module      : Balanced
Description : Amplitude-balanced path sums
Copyright   : (c) Matthew Amy, 2020
Maintainer  : matt.e.amy@gmail.com
Stability   : experimental
Portability : portable
-}

module Feynman.Algebra.Pathsum.Balanced where

import Data.List
import qualified Data.Set as Set
import Data.Ratio
import Data.Semigroup
import Control.Monad (mzero, msum)
import Data.Maybe (maybeToList)
import Data.Complex (Complex, mkPolar)
import Data.Bits (shiftL)
import Data.Map (Map, (!))
import qualified Data.Map as Map
import Data.String (IsString(..))
import Data.Tuple (swap)

import qualified Feynman.Util.Unicode as U
import Feynman.Algebra.Base
import Feynman.Algebra.Polynomial (degree)
import Feynman.Algebra.Polynomial.Multilinear

{-----------------------------------
 Variables
 -----------------------------------}

-- | Variables are either input variables or path variables. The distinction
--   is due to the binding structure of our pathsum representation, and moreover
--   improves readability
data Var = IVar !Int | PVar !Int | FVar !String deriving (Eq, Ord)

instance Show Var where
  show (IVar i) = U.sub "x" $ fromIntegral i
  show (PVar i) = U.sub "y" $ fromIntegral i
  show (FVar x) = x

instance IsString Var where
  fromString = FVar

-- | Convenience function for the string representation of the 'i'th input variable
ivar :: Int -> String
ivar = show . IVar

-- | Convenience function for the string representation of the 'i'th path variable
pvar :: Int -> String
pvar = show . PVar

-- | Construct an integer shift for input variables
shiftI :: Int -> (Var -> Var)
shiftI i = shiftAll i 0

-- | Construct an integer shift for path variables
shiftP :: Int -> (Var -> Var)
shiftP j = shiftAll 0 j

-- | General shift. Constructs a substitution from shift values for I and P
shiftAll :: Int -> Int -> (Var -> Var)
shiftAll i j = go
  where go (IVar i') = IVar (i + i')
        go (PVar j') = PVar (j + j')

-- | Check if a variable is an input
isI :: Var -> Bool
isI (IVar _) = True
isI _        = False

-- | Check if a variable is a path variable
isP :: Var -> Bool
isP (PVar _) = True
isP _        = False

-- | Check if a variable is a free variable
isF :: Var -> Bool
isF (FVar _) = True
isF _        = False

-- | Get the string of a free variable
unF :: Var -> String
unF (FVar s) = s
unF _        = error "Not a free variable"

{-----------------------------------
 Path sums
 -----------------------------------}

-- | Path sums of the form
--   \(\frac{1}{\sqrt{2}^k}\sum_{y\in\mathbb{Z}_2^m}e^{i\pi P(x, y)}|f(x, y)\rangle\)
data Pathsum g = Pathsum {
  sde       :: !Int,
  inDeg     :: !Int,
  outDeg    :: !Int,
  pathVars  :: !Int,
  phasePoly :: !(PseudoBoolean Var g),
  outVals   :: ![SBool Var]
  } deriving (Eq)

instance (Show g, Eq g, Periodic g, Real g) => Show (Pathsum g) where
  show sop = inputstr ++ scalarstr ++ sumstr ++ amplitudestr ++ statestr
    where inputstr = case inDeg sop of
            0 -> ""
            1 -> U.ket (ivar 0) ++ " " ++ U.mapsto ++ " "
            2 -> U.ket (ivar 0) ++ U.ket (ivar 1) ++ " " ++ U.mapsto ++ " "
            j -> U.ket (ivar 0) ++ U.dots ++ U.ket (ivar (j-1)) ++ " " ++ U.mapsto ++ " "
          scalarstr = case compare (sde sop) 0 of
            LT -> U.sup ("(" ++ U.rt2 ++ ")") (fromIntegral . abs $ sde sop)
            EQ -> if (inDeg sop == 0 && outDeg sop == 0 && phasePoly sop == 0) then "1" else ""
            GT -> U.sup ("1/(" ++ U.rt2 ++ ")") (fromIntegral $ sde sop)
          sumstr = case pathVars sop of
            0 -> ""
            1 -> U.sum ++ "[" ++ pvar 0 ++ "]"
            2 -> U.sum ++ "[" ++ pvar 0 ++ pvar 1 ++ "]"
            j -> U.sum ++ "[" ++ pvar 0 ++ U.dots ++ pvar (j-1) ++ "]"
          amplitudestr = case order (phasePoly sop) of
            0 -> U.e ++ "^" ++ U.i ++ U.pi ++ "{" ++ show (phasePoly sop) ++ "}"
            1 -> ""
            2 -> "(-1)^{" ++ show (makeIntegral 1 $ phasePoly sop) ++ "}"
            4 -> U.i ++ "^{" ++ show (makeIntegral 2 $ phasePoly sop) ++ "}"
            8 -> U.omega ++ "^{" ++ show (makeIntegral 4 $ phasePoly sop) ++ "}"
            j -> U.sub U.zeta j ++ "^{" ++ show (makeIntegral j $ phasePoly sop) ++ "}"
          statestr = concatMap (U.ket . show) $ outVals sop

-- | Convenience function for pretty printing
makeIntegral :: Real g => Integer -> PseudoBoolean v g -> PseudoBoolean v Integer
makeIntegral i = cast (\a -> numerator $ toRational a * toRational i)

-- | Retrieve the internal path variables
internalPaths :: Pathsum g -> [Var]
internalPaths sop = [PVar i | i <- [0..pathVars sop - 1]] \\ outVars
  where outVars = Set.toList . Set.unions . map vars $ outVals sop

-- | Retrieve the free variables
freeVars :: Pathsum g -> [String]
freeVars sop = map unF . Set.toList . Set.filter isF . foldr (Set.union) Set.empty $ xs
  where xs = (vars $ phasePoly sop):(map vars $ outVals sop)

-- | Checks if the path sum is (trivially) the identity
isTrivial :: (Eq g, Num g) => Pathsum g -> Bool
isTrivial sop = sop == identity (inDeg sop)

-- | (To be deprecated) Drops constant term from the phase polynomial
dropGlobalPhase :: (Eq g, Num g) => Pathsum g -> Pathsum g
dropGlobalPhase sop = sop { phasePoly = dropConstant $ phasePoly sop }

-- | (To be deprecated) Drops the phase polynomial
dropPhase :: (Eq g, Num g) => Pathsum g -> Pathsum g
dropPhase sop = sop { phasePoly = 0 }

-- | (To be deprecated) Drops normalization
dropAmplitude :: Pathsum g -> Pathsum g
dropAmplitude sop = sop { sde = 0 }

{----------------------------
 Constructors
 ----------------------------}

-- | Construct an 'n'-qubit identity operator
identity :: (Eq g, Num g) => Int -> Pathsum g
identity n = Pathsum 0 n n 0 0 [ofVar (IVar i) | i <- [0..n-1]]

-- | Construct a (symbolic) state
ket :: (Eq g, Num g) => [SBool String] -> Pathsum g
ket xs = Pathsum 0 0 (fromIntegral $ length xs) 0 0 $ map (rename FVar) xs

-- | Initialize a fresh ancilla
initialize :: (Eq g, Num g) => FF2 -> Pathsum g
initialize b = ket [constant b]

{-# INLINE initialize #-}

-- | Construct a uniform superposition of classical states of a given form
superposition :: (Ord v, Eq g, Num g) => [SBool v] -> Pathsum g
superposition xs = Pathsum k 0 n k 0 $ map (rename sub) xs
  where n   = length xs
        k   = Set.size fv
        fv  = Set.unions $ map vars xs
        sub = ((Map.fromList [(v, PVar i) | (v, i) <- zip (Set.toList fv) [0..]])!)

-- | Construct a classical transformation from the free variables of 'xs' to 'xs'
--   Effectively binds the free variables in 'state xs'
compute :: (Ord v, Eq g, Num g) => [SBool v] -> Pathsum g
compute xs = Pathsum 0 (Set.size fv) (length xs) 0 0 $ map (rename sub) xs
  where fv  = Set.unions $ map vars xs
        sub = ((Map.fromList [(v, IVar i) | (v, i) <- zip (Set.toList fv) [0..]])!)

-- | Breaks the connection between inputs and outputs by mapping
--   any input basis state to the maximally mixed state. Non-unitary
disconnect :: (Eq g, Num g) => Int -> Pathsum g
disconnect n = Pathsum n n n n 0 [ofVar $ PVar i | i <- [0..n-1]]

-- | Construct a permutation
permutation :: (Eq g, Num g) => [Int] -> Pathsum g
permutation xs
  | all (\i -> i >= 0 && i <= n) xs = Pathsum 0 n n 0 0 (map (ofVar . IVar) xs)
  | otherwise = error "permutation: Input not a permutation"
  where n = length xs

-- | Construct an n-ary unit
etaN :: (Eq g, Num g) => Int -> Pathsum g
etaN n = Pathsum 0 0 (2*n) n 0 $ xs ++ xs
  where xs = map (ofVar . PVar) [0..n-1]

{----------------------------
 Dual constructors
 ----------------------------}

-- | Construct a (symbolic) state destructor
bra :: (Eq g, Abelian g) => [SBool String] -> Pathsum g
bra xs = Pathsum (2*m) m 0 m (lift $ p) []
  where m         = fromIntegral $ length xs
        p         = foldr (+) 0 . map go $ zip [0..] xs
        go (i, v) = ofVar (PVar i) * (ofVar (IVar i) + (rename fromString v))

-- | Alternate state destructor with fewer paths but a more complicated polynomial
unstateAlt :: (Eq g, Abelian g) => [SBool String] -> Pathsum g
unstateAlt xs = Pathsum 2 (fromIntegral $ length xs) 0 1 (lift $ y*(1 + p)) []
  where y = ofVar (PVar 0)
        p = foldr (*) 1 . map valF $ zip xs [0..]
        valF (val, i) = 1 + (rename fromString val) + ofVar (IVar i)

-- | Dagger of initialize -- i.e. unnormalized post-selection
postselect :: (Eq g, Abelian g) => FF2 -> Pathsum g
postselect b = bra [constant b]

{-# INLINE postselect #-}

-- | Select on a classical state of a given form
unsuper :: (Ord v, Eq g, Abelian g) => [SBool v] -> Pathsum g
unsuper xs = Pathsum (2*m + n) m 0 (m+n) poly []
  where m    = length xs
        n    = Set.size fv
        fv   = Set.unions $ map vars xs
        sub  = ((Map.fromList [(v, PVar (m + i)) | (v, i) <- zip (Set.toList fv) [0..]])!)
        poly = foldr (+) zero $ map constructTerm [0..m-1]
        constructTerm i = lift $ ofVar (PVar i) * (ofVar (IVar i) + rename sub (xs!!i))

-- | Invert a classical transformation
uncompute :: (Ord v, Eq g, Abelian g) => [SBool v] -> Pathsum g
uncompute xs = Pathsum (2*m) m n (m+n) poly [ofVar (PVar $ m + i) | i <- [0..n-1]]
  where m    = length xs
        n    = Set.size fv
        fv   = Set.unions $ map vars xs
        sub  = ((Map.fromList [(v, PVar (m + i)) | (v, i) <- zip (Set.toList fv) [0..]])!)
        poly = foldr (+) zero $ map constructTerm [0..m-1]
        constructTerm i = lift $ ofVar (PVar i) * (ofVar (IVar i) + rename sub (xs!!i))

-- | Construct an inverse permutation
unpermutation :: (Eq g, Num g) => [Int] -> Pathsum g
unpermutation = permutation . snd . unzip . sort . map swap . zip [0..]

-- | Construct an n-ary co-unit
epsilonN :: (Eq g, Abelian g) => Int -> Pathsum g
epsilonN n = Pathsum (2*n) (2*n) 0 n (lift poly) []
  where poly = sum $ map f [0..n-1]
        f i  = ofVar (PVar i) * (ofVar (IVar i) + ofVar (IVar $ n+i))

{----------------------------
 Constants & gates
 ----------------------------}

-- | \(\sqrt{2}\)
root2 :: (Eq g, Abelian g, Dyadic g) => Pathsum g
root2 = Pathsum 0 0 0 1 ((-constant (half * half)) + scale half (lift $ ofVar (PVar 0))) []

-- | \(1/\sqrt{2}\)
roothalf :: (Eq g, Abelian g, Dyadic g) => Pathsum g
roothalf = Pathsum 1 0 0 0 0 []

-- | \(i\)
iunit :: (Eq g, Abelian g, Dyadic g) => Pathsum g
iunit = Pathsum 0 0 0 0 (constant half) []

-- | \(e^{i\pi/4}\)
omega :: (Eq g, Abelian g, Dyadic g) => Pathsum g
omega = Pathsum 0 0 0 0 (constant (half * half)) []

-- | A fresh, 0-valued ancilla
fresh :: (Eq g, Num g) => Pathsum g
fresh = Pathsum 0 0 1 0 0 [0]

-- | The dagger of fresh
unfresh :: (Eq g, Abelian g) => Pathsum g
unfresh = dagger fresh

-- | The unit, \(\eta\)
eta :: (Eq g, Num g) => Pathsum g
eta = Pathsum 0 0 2 1 0 [ofVar (PVar 0), ofVar (PVar 0)]

-- | The co-unit, \(\epsilon\)
epsilon :: (Eq g, Abelian g) => Pathsum g
epsilon = Pathsum 2 2 0 1 p []
  where p = lift $ ofVar (PVar 0) * (ofVar (IVar 0) + ofVar (IVar 1))

-- | X gate
xgate :: (Eq g, Num g) => Pathsum g
xgate = Pathsum 0 1 1 0 0 [1 + ofVar (IVar 0)]

-- | Z gate
zgate :: (Eq g, Abelian g) => Pathsum g
zgate = Pathsum 0 1 1 0 p [ofVar (IVar 0)]
  where p = lift $ ofVar (IVar 0)

-- | Y gate
ygate :: (Eq g, Abelian g, Dyadic g) => Pathsum g
ygate = Pathsum 0 1 1 0 p [1 + ofVar (IVar 0)]
  where p = constant half + (lift $ ofVar (IVar 0))

-- | S gate
sgate :: (Eq g, Abelian g, Dyadic g) => Pathsum g
sgate = Pathsum 0 1 1 0 p [ofVar (IVar 0)]
  where p = scale half (lift $ ofVar (IVar 0))

-- | S* gate
sdggate :: (Eq g, Abelian g, Dyadic g) => Pathsum g
sdggate = Pathsum 0 1 1 0 p [ofVar (IVar 0)]
  where p = scale (-half) (lift $ ofVar (IVar 0))

-- | T gate
tgate :: (Eq g, Abelian g, Dyadic g) => Pathsum g
tgate = Pathsum 0 1 1 0 p [ofVar (IVar 0)]
  where p = scale (half*half) (lift $ ofVar (IVar 0))

-- | T* gate
tdggate :: (Eq g, Abelian g, Dyadic g) => Pathsum g
tdggate = Pathsum 0 1 1 0 p [ofVar (IVar 0)]
  where p = scale (-half*half) (lift $ ofVar (IVar 0))

-- | R_k gate
rkgate :: (Eq g, Abelian g, Dyadic g) => Int -> Pathsum g
rkgate k = Pathsum 0 1 1 0 p [ofVar (IVar 0)]
  where p = scale (fromDyadic $ dyadic 1 k) (lift $ ofVar (IVar 0))

-- | R_z gate
rzgate :: (Eq g, Abelian g, Dyadic g) => DyadicRational -> Pathsum g
rzgate theta = Pathsum 0 1 1 0 p [ofVar (IVar 0)]
  where p = scale (fromDyadic theta) (lift $ ofVar (IVar 0))

-- | H gate
hgate :: (Eq g, Abelian g, Dyadic g) => Pathsum g
hgate = Pathsum 1 1 1 1 p [ofVar (PVar 0)]
  where p = lift $ (ofVar $ IVar 0) * (ofVar $ PVar 0)

-- | CH gate
chgate :: (Eq g, Abelian g, Dyadic g) => Pathsum g
chgate = Pathsum 1 2 2 1 p [x1, x2 + x1*x2 + x1*y]
  where p = lift $ x1 * x2 * y
        x1 = ofVar $ IVar 0
        x2 = ofVar $ IVar 1
        y = ofVar $ PVar 0

-- | CNOT gate
cxgate :: (Eq g, Num g) => Pathsum g
cxgate = Pathsum 0 2 2 0 0 [x0, x0+x1]
  where x0 = ofVar $ IVar 0
        x1 = ofVar $ IVar 1

-- | Toffoli gate
ccxgate :: (Eq g, Num g) => Pathsum g
ccxgate = Pathsum 0 3 3 0 0 [x0, x1, x2 + x0*x1]
  where x0 = ofVar $ IVar 0
        x1 = ofVar $ IVar 1
        x2 = ofVar $ IVar 2

-- | k-control Toffoli gate
mctgate :: (Eq g, Num g) => Int -> Pathsum g
mctgate k = Pathsum 0 (k+1) (k+1) 0 0 (controls ++ [t + foldr (*) 1 controls])
  where controls = [ofVar (IVar i) | i <- [0..k-1]]
        t        = ofVar $ IVar k

-- | SWAP gate
swapgate :: (Eq g, Num g) => Pathsum g
swapgate = Pathsum 0 2 2 0 0 [x1, x0]
  where x0 = ofVar $ IVar 0
        x1 = ofVar $ IVar 1

{----------------------------
 Channels
 ----------------------------}

-- | Choi matrix of computational basis measurement
measureChoi :: (Eq g, Abelian g) => Pathsum g
measureChoi = Pathsum 2 2 2 1 (lift $ y * (x0 + x1)) [x0, x1]
  where x0 = ofVar $ IVar 0
        x1 = ofVar $ IVar 1
        y  = ofVar $ PVar 0

-- | CPM operator of computational basis measurement
measure :: (Eq g, Abelian g) => Pathsum g
measure = unChoi measureChoi

{----------------------------
 Bind, unbind, and subst
 ----------------------------}

-- | Bind some collection of free variables in a path sum
bind :: (Foldable f, Eq g, Abelian g) => f String -> Pathsum g -> Pathsum g
bind = flip (foldr go)
  where go x sop =
          let v = IVar $ inDeg sop in
            sop { inDeg = (inDeg sop) + 1,
                  phasePoly = subst (FVar x) (ofVar v) (phasePoly sop),
                  outVals = map (subst (FVar x) (ofVar v)) (outVals sop) }

-- | Close a path sum by binding all free variables
close :: (Eq g, Abelian g) => Pathsum g -> Pathsum g
close sop = bind (freeVars sop) sop

-- | Unbind (instantiate) some collection of inputs
unbind :: (Foldable f, Eq g, Abelian g) => f Int -> Pathsum g -> Pathsum g
unbind xs (Pathsum a b c d e f) = Pathsum a (b - length xs) c d e' f' where
  e'  = substMany sub e
  f'  = map (substMany sub) f
  sub = \v -> Map.findWithDefault (ofVar v) v tmp
  tmp = snd $ foldr buildMap (0, Map.empty) [0..b-1]
  buildMap i (j, acc) = case i `elem` xs of
    True  -> (j, Map.insert (IVar i) (ofVar . FVar $ "#" ++ show (IVar i)) acc)
    False -> (j+1, Map.insert (IVar i) (ofVar . IVar $ j) acc)

-- | Open a path sum by instantiating all inputs
open :: (Eq g, Abelian g) => Pathsum g -> Pathsum g
open sop = unbind [0..(inDeg sop) - 1] sop

-- | Substitute a monomial with a symbolic Boolean expression throughout
--
--   This is generally not a very safe thing to do. Convenience for certain
--   local transformations
substitute :: (Eq g, Abelian g) => [Var] -> SBool Var -> Pathsum g -> Pathsum g
substitute xs p (Pathsum a b c d e f) = Pathsum a b c d e' f' where
  e' = substMonomial xs p e
  f' = map (substMonomial xs p) f

{----------------------------
 Operators
 ----------------------------}

-- | Return the dual of a path sum
dualize :: (Eq g, Abelian g) => Pathsum g -> Pathsum g
dualize sop@(Pathsum a b c d e f) = inSOP .> midSOP .> outSOP
  where inSOP  = tensor (identity c) (etaN b)
        midSOP = tensor (tensor (identity c) sop) (identity b)
        outSOP = tensor (epsilonN c) (identity b)

-- | Return the (column) vectorized path sum. By convention we place the inputs
--   first (i.e. f : A -> B becomes vectorize f : A* \otimes B)
vectorize :: (Eq g, Abelian g) => Pathsum g -> Pathsum g
vectorize sop@(Pathsum a b c d e f) = etaN b .> tensor (identity b) sop

-- | Return the (row) vectorized path sum. By convention we place the outputs
--   first (i.e. f : A -> B becomes vectorize f : A* \otimes B)
covectorize :: (Eq g, Abelian g) => Pathsum g -> Pathsum g
covectorize sop@(Pathsum a b c d e f) = tensor (identity c) sop .> epsilonN c

-- | Take the dagger of a path sum
dagger :: (Eq g, Abelian g) => Pathsum g -> Pathsum g
dagger (Pathsum a b c d e f) = dualize $ Pathsum a b c d (-e) f

-- | Take the conjugate (c.f., lower star) of a path sum
conjugate :: (Eq g, Abelian g) => Pathsum g -> Pathsum g
conjugate (Pathsum a b c d e f) = Pathsum a b c d (-e) f

-- | Trace a square morphism. Throws an error if the input and outputs are not
--   the same size
trace :: (Eq g, Abelian g) => Pathsum g -> Pathsum g
trace sop@(Pathsum a b c d e f)
  | b /= c = error "Can't trace a non-square operator"
  | otherwise = etaN b .> tensor (identity b) sop .> epsilonN b

-- | Trace out the first qubit. Throws an error if the input is a vector
ptrace :: (Eq g, Abelian g) => Pathsum g -> Pathsum g
ptrace sop@(Pathsum a b c d e f)
  | b < 1 || c < 1 = error "Can't partial trace a vector"
  | otherwise = tensor eta (identity $ b-1) .>
                tensor (identity 1) sop .>
                tensor epsilon (identity $ c-1)

-- | Turn a pure state into a density matrix. Throws an error if the input is not
--   a column vector
densify :: (Eq g, Abelian g) => Pathsum g -> Pathsum g
densify sop@(Pathsum a b c d e f)
  | b /= 0 = error "Can't densify an operator "
  | otherwise = dagger sop .> sop

-- | Turn a (single-qubit) Choi matrix into a CPM-style linear operator
unChoi :: (Eq g, Abelian g) => Pathsum g -> Pathsum g
unChoi sop@(Pathsum a b c d e f)
  | b /= 2 || c /= 2 = error "Only single-qubit channels currently supported"
  | otherwise        = tensor eta swapgate .>
                       tensor (identity 1) (tensor sop $ identity 1) .>
                       tensor swapgate epsilon

-- | Turn a unitary operator into a channel. That is, f becomes f_* \otimes f
channelize :: (Eq g, Abelian g) => Pathsum g -> Pathsum g
channelize sop = tensor (conjugate sop) sop

-- | Construct a controlled path sum
controlled :: (Eq g, Abelian g) => Pathsum g -> Pathsum g
controlled sop@(Pathsum a b c d e f) = Pathsum a (b+1) (c+1) d e' f' where
  shift   = shiftI 1
  x       = ofVar $ IVar 0
  e'      = (lift x)*(renameMonotonic shift e)
  f'      = [lift x] ++ (map g . zip [1..] . map (renameMonotonic shift) $ f)
  g (i,y) = (ofVar $ IVar i) + x*((ofVar $ IVar i) + y)

-- | Attempt to add two path sums. Only succeeds if the resulting sum is balanced
--   and the dimensions match.
plusMaybe :: (Eq g, Abelian g) => Pathsum g -> Pathsum g -> Maybe (Pathsum g)
plusMaybe sop sop'
  | inDeg sop  /= inDeg sop'                                       = Nothing
  | outDeg sop /= outDeg sop'                                      = Nothing
  | (sde sop) + 2*(pathVars sop') /= (sde sop') + 2*(pathVars sop) = Nothing
  | otherwise = Just $ Pathsum sde' inDeg' outDeg' pathVars' phasePoly' outVals'
  where sde'       = (sde sop) + 2*(pathVars sop')
        inDeg'     = inDeg sop
        outDeg'    = outDeg sop
        pathVars'  = (pathVars sop) + (pathVars sop') + 1
        y          = ofVar $ PVar (pathVars' - 1)
        phasePoly' = (lift y)*(phasePoly sop) +
                     (lift (1+y))*(renameMonotonic shift $ phasePoly sop')
        outVals'   = map (\(a,b) -> b + y*(a + b)) $
                       zip (outVals sop) (map (renameMonotonic shift) $ outVals sop')
        shift x    = case x of
          PVar i -> PVar $ i + (pathVars sop)
          _      -> x

-- | Construct the sum of two path sums. Raises an error if the sums are incompatible
plus :: (Eq g, Abelian g) => Pathsum g -> Pathsum g -> Pathsum g
plus sop sop' = case plusMaybe sop sop' of
  Nothing    -> error "Incompatible path sums"
  Just sop'' -> sop''

-- | Compose two path sums in parallel
tensor :: (Eq g, Num g) => Pathsum g -> Pathsum g -> Pathsum g
tensor sop sop' = Pathsum sde' inDeg' outDeg' pathVars' phasePoly' outVals'
  where sde'       = (sde sop) + (sde sop')
        inDeg'     = (inDeg sop) + (inDeg sop')
        outDeg'    = (outDeg sop) + (outDeg sop')
        pathVars'  = (pathVars sop) + (pathVars sop')
        phasePoly' = (phasePoly sop) + (renameMonotonic shift $ phasePoly sop')
        outVals'   = (outVals sop) ++ (map (renameMonotonic shift) $ outVals sop')
        shift x    = case x of
          IVar i -> IVar $ i + (inDeg sop)
          PVar i -> PVar $ i + (pathVars sop)
          _      -> x

-- | Attempt to compose two path sums in sequence. Only succeeds if the dimensions
--   are compatible (i.e. if the out degree of the former is the in degree of the
--   latter)
timesMaybe :: (Eq g, Abelian g) => Pathsum g -> Pathsum g -> Maybe (Pathsum g)
timesMaybe sop sop'
  | outDeg sop /= inDeg sop' = Nothing
  | otherwise = Just $ Pathsum sde' inDeg' outDeg' pathVars' phasePoly' outVals'
  where sde'       = (sde sop) + (sde sop')
        inDeg'     = inDeg sop
        outDeg'    = outDeg sop'
        pathVars'  = (pathVars sop) + (pathVars sop')
        phasePoly' = (phasePoly sop) +
                     (substMany sub . renameMonotonic shift $ phasePoly sop')
        outVals'   = (map (substMany sub . renameMonotonic shift) $ outVals sop')
        shift x    = case x of
          PVar i -> PVar $ i + (pathVars sop)
          _      -> x
        sub x      = case x of
          IVar i -> (outVals sop)!!i
          _      -> ofVar x

-- | Compose two path sums in sequence. Throws an error if the dimensions are
--   not compatible
times :: (Eq g, Abelian g) => Pathsum g -> Pathsum g -> Pathsum g
times sop sop' = case timesMaybe sop sop' of
  Nothing    -> error "Incompatible path sum dimensions"
  Just sop'' -> sop''

-- | Left-to-right composition
(.>) :: (Eq g, Abelian g) => Pathsum g -> Pathsum g -> Pathsum g
(.>) = times

infixr 5 .>

-- | Scale the normalization factor
renormalize :: Int -> Pathsum g -> Pathsum g
renormalize k (Pathsum a b c d e f) = Pathsum (a + k) b c d e f

-- | Embed a path sum into a larger space with a specified input and
--   output embedding.
embed :: (Eq g, Abelian g) => Pathsum g -> Int -> (Int -> Int) -> (Int -> Int) -> Pathsum g
embed sop n embedIn embedOut
  | n < 0     = error "Can't embed in smaller space"
  | otherwise = inPerm .> tensor (identity n) sop .> outPerm where
      mIn = inDeg sop
      ins = map embedIn [0..mIn-1]
      inPerm  = permutation $ ([0..mIn+n-1] \\ ins) ++ ins
      mOut = outDeg sop
      outs = map embedOut [0..mOut-1]
      outPerm = unpermutation $ ([0..mOut+n-1] \\ outs) ++ outs

-- | Drop a qubit
discard :: Eq g => Int -> Pathsum g -> Pathsum g
discard i sop@(Pathsum a b c d e f) = Pathsum a b' c' d e f' where
  b' = if i < b then b-1 else b
  c' = if i < c then c-1 else c
  f' = snd . unzip . filter (\(j,_) -> i /= j) $ zip [0..] f

{--------------------------
 Type class instances
 --------------------------}

instance (Eq g, Num g) => Semigroup (Pathsum g) where
  (<>) = tensor

instance (Eq g, Num g) => Monoid (Pathsum g) where
  mempty  = Pathsum 0 0 0 0 0 []
  mappend = tensor

instance (Eq g, Abelian g) => Num (Pathsum g) where
  (+)                          = plus
  (*)                          = (flip times)
  negate (Pathsum a b c d e f) = Pathsum a b c d (lift 1 + e) f
  abs (Pathsum a b c d e f)    = Pathsum a b c d (dropConstant e) f
  signum sop                   = sop
  fromInteger                  = identity . fromInteger

instance Functor Pathsum where
  fmap g (Pathsum a b c d e f) = Pathsum a b c d (cast g e) f

{--------------------------
 Reduction rules
 --------------------------}

{-
class RewriteRule rule g where
  matchAll :: Pathsum g -> [rule]
  matchOne :: Pathsum g -> rule
  apply :: rule -> Pathsum g -> Pathsum g

data E = E !Var {-# UNBOX #-}

instance (Eq g, Periodic g) => RewriteRule E g where
  match = matchElim
  apply = applyElim
-}

-- | Maps the order 1 and order 2 elements of a group to FF2
injectFF2 :: Periodic g => g -> Maybe FF2
injectFF2 a = case order a of
  1 -> Just 0
  2 -> Just 1
  _ -> Nothing

-- | Gives a Boolean polynomial equivalent to the current polynomial, if possible
toBooleanPoly :: (Eq g, Periodic g) => PseudoBoolean v g -> Maybe (SBool v)
toBooleanPoly = castMaybe injectFF2

-- | Elim rule. \(\dots(\sum_y)\dots = \dots 2 \dots\)
matchElim :: (Eq g, Periodic g) => Pathsum g -> [Var]
matchElim sop = msum . (map go) $ internalPaths sop
  where go v = if Set.member v (vars $ phasePoly sop) then [] else [v]

-- | Generic HH rule. \(\dots(\sum_y (-1)^{y\cdot f})\dots = \dots|_{f = 0}\)
matchHH :: (Eq g, Periodic g) => Pathsum g -> [(Var, SBool Var)]
matchHH sop = msum . (map (maybeToList . go)) $ internalPaths sop
  where go v = toBooleanPoly (quotVar v $ phasePoly sop) >>= \p -> return (v, p)

-- | Solvable instances of the HH rule.
--   \(\dots(\sum_y (-1)^{y(z \oplus f)})\dots = \dots[z \gets f]\)
matchHHSolve :: (Eq g, Periodic g) => Pathsum g -> [(Var, Var, SBool Var)]
matchHHSolve sop = do
  (v, p)   <- matchHH sop
  (v', p') <- solveForX p
  case v' of
    PVar _ -> return (v, v', p')
    _      -> mzero

-- | Instances of the HH rule with a linear substitution
matchHHLinear :: (Eq g, Periodic g) => Pathsum g -> [(Var, Var, SBool Var)]
matchHHLinear sop = do
  (v, p)   <- filter (\(_, p) -> degree p <= 1) $ matchHH sop
  (v', p') <- solveForX p
  return (v, v', p')

-- | Instances of the (\omega\) rule
matchOmega :: (Eq g, Periodic g, Dyadic g) => Pathsum g -> [(Var, SBool Var)]
matchOmega sop = do
  v <- internalPaths sop
  p <- maybeToList . toBooleanPoly . addFactor v $ phasePoly sop
  return (v, p)
  where addFactor v p = constant (fromDyadic $ dyadic 3 1) + quotVar v p

{--------------------------
 Pattern synonyms
 --------------------------}

-- | Pattern synonym for Elim
pattern Triv :: (Eq g, Num g) => Pathsum g
pattern Triv <- (isTrivial -> True)

-- | Pattern synonym for Elim
pattern Elim :: (Eq g, Periodic g) => Var -> Pathsum g
pattern Elim v <- (matchElim -> (v:_))

-- | Pattern synonym for HH
pattern HH :: (Eq g, Periodic g) => Var -> SBool Var -> Pathsum g
pattern HH v p <- (matchHH -> (v, p):_)

-- | Pattern synonym for solvable HH instances
pattern HHSolved :: (Eq g, Periodic g) => Var -> Var -> SBool Var -> Pathsum g
pattern HHSolved v v' p <- (matchHHSolve -> (v, v', p):_)

-- | Pattern synonym for linear HH instances
pattern HHLinear :: (Eq g, Periodic g) => Var -> Var -> SBool Var -> Pathsum g
pattern HHLinear v v' p <- (matchHHLinear -> (v, v', p):_)

-- | Pattern synonym for HH instances where the polynomial is strictly a
--   function of input variables
pattern HHKill :: (Eq g, Periodic g) => Var -> SBool Var -> Pathsum g
pattern HHKill v p <- (filter (all (not . isP) . vars . snd) . matchHH -> (v, p):_)

-- | Pattern synonym for Omega instances
pattern Omega :: (Eq g, Periodic g, Dyadic g) => Var -> SBool Var -> Pathsum g
pattern Omega v p <- (matchOmega -> (v, p):_)

{--------------------------
 Applying reductions
 --------------------------}

-- | Apply an elim rule. Does not check if the instance is valid
applyElim :: Var -> Pathsum g -> Pathsum g
applyElim (PVar i) (Pathsum a b c d e f) = Pathsum (a-2) b c (d-1) e' f'
  where e' = renameMonotonic varShift e
        f' = map (renameMonotonic varShift) f
        varShift (PVar j)
          | j > i     = PVar $ j - 1
          | otherwise = PVar $ j
        varShift v = v

-- | Apply a (solvable) HH rule. Does not check if the instance is valid
applyHHSolved :: (Eq g, Abelian g) => Var -> Var -> SBool Var -> Pathsum g -> Pathsum g
applyHHSolved (PVar i) v p (Pathsum a b c d e f) = Pathsum a b c (d-1) e' f'
  where e' = renameMonotonic varShift . subst v p . remVar (PVar i) $ e
        f' = map (renameMonotonic varShift . subst v p) f
        varShift (PVar j)
          | j > i     = PVar $ j - 1
          | otherwise = PVar $ j
        varShift v = v

-- | Apply an (\omega\) rule. Does not check if the instance is valid
applyOmega :: (Eq g, Abelian g, Dyadic g) => Var -> SBool Var -> Pathsum g -> Pathsum g
applyOmega (PVar i) p (Pathsum a b c d e f) = Pathsum (a-1) b c (d-1) e' f'
  where e' = renameMonotonic varShift $ p' + remVar (PVar i) e
        f' = map (renameMonotonic varShift) f
        p' = constant (fromDyadic $ dyadic 1 2) + distribute (fromDyadic $ dyadic 3 1) (lift p)
        varShift (PVar j)
          | j > i     = PVar $ j - 1
          | otherwise = PVar $ j
        varShift v = v

-- | Finds and applies the first elimination instance
rewriteElim :: (Eq g, Periodic g) => Pathsum g -> Pathsum g
rewriteElim sop = case sop of
  Elim v -> applyElim v sop
  _      -> sop

-- | Finds and applies the first hh instance
rewriteHH :: (Eq g, Periodic g) => Pathsum g -> Pathsum g
rewriteHH sop = case sop of
  HHSolved v v' p -> applyHHSolved v v' p sop
  _               -> sop

-- | Finds and applies the first omega instance
rewriteOmega :: (Eq g, Periodic g, Dyadic g) => Pathsum g -> Pathsum g
rewriteOmega sop = case sop of
  Omega v p -> applyOmega v p sop
  _         -> sop

{--------------------------
 Reduction procedures
 --------------------------}

-- | Performs basic simplifications
simplify :: (Eq g, Periodic g, Dyadic g) => Pathsum g -> Pathsum g
simplify sop = case sop of
  Elim y         -> grind $ applyElim y sop
  HHLinear y z p -> grind $ applyHHSolved y z p sop
  Omega y p      -> grind $ applyOmega y p sop
  _              -> sop

-- | A complete normalization procedure for Clifford circuits. Originally described in
--   the paper M. Amy,
--   / Towards Large-Scaled Functional Verification of Universal Quantum Circuits /, QPL 2018.
grind :: (Eq g, Periodic g, Dyadic g) => Pathsum g -> Pathsum g
grind sop = case sop of
  Elim y         -> grind $ applyElim y sop
  HHSolved y z p -> grind $ applyHHSolved y z p sop
  Omega y p      -> grind $ applyOmega y p sop
  _              -> sop

-- | A single step of 'grind'
grindStep :: (Eq g, Periodic g, Dyadic g) => Pathsum g -> Pathsum g
grindStep sop = case sop of
  Elim y         -> applyElim y sop
  HHSolved y z p -> applyHHSolved y z p sop
  Omega y p      -> applyOmega y p sop
  _              -> sop

{--------------------------
 Simulation
 --------------------------}

-- | Gets the cofactors of some path variable
expand :: (Eq g, Abelian g) => Pathsum g -> Var -> (Pathsum g, Pathsum g)
expand (Pathsum a b c d e f) v = (p0, p1) where
  p0  = Pathsum a b c (d-1) (subst v 0 e) (map (subst v 0) f)
  p1  = Pathsum a b c (d-1) (subst v 1 e) (map (subst v 1) f)

-- | Simulates a pathsum on a given input
simulate :: (Eq g, Periodic g, Dyadic g, Real g, RealFloat f) => Pathsum g -> [FF2] -> Map [FF2] (Complex f)
simulate sop xs = go $ sop * ket (map constant xs)
  where go      = go' . grind
        go' ps  = case ps of
          (Pathsum k 0 _ 0 p xs) ->
            let phase     = fromRational . toRational $ getConstant p
                base      = case k `mod` 2 of
                  0 -> fromInteger $ 1 `shiftL` (abs k)
                  1 -> sqrt(2.0) * (fromInteger $ 1 `shiftL` (abs (k-1)))
                magnitude = base**(fromIntegral $ signum k)
            in
              Map.singleton (map getConstant xs) (mkPolar magnitude (pi * phase))
          (Pathsum k 0 n i p xs) ->
            let v     = PVar $ i-1
                left  = go (Pathsum k 0 n (i-1) (subst v zero p) (map (subst v zero) xs))
                right = go (Pathsum k 0 n (i-1) (subst v one p) (map (subst v one) xs))
            in
              Map.unionWith (+) left right
          _                      -> error "Incompatible dimensions"

-- | Evaluates a pathsum on a given input and output
amplitude :: (Eq g, Periodic g, Dyadic g, Real g, RealFloat f) => [FF2] -> Pathsum g -> [FF2] -> Complex f
amplitude o sop i = (simulate (bra (map constant o) * sop) i)![]

-- | Checks identity by checking inputs iteratively
isIdentity :: (Eq g, Periodic g, Dyadic g) => Pathsum g -> Bool
isIdentity sop
  | isTrivial sop = True
  | otherwise     = case inDeg sop of
      0 -> False
      i -> let p0 = (grind $ identity (i-1) <> ket [0] .> sop .> identity (i-1) <> bra [0])
               p1 = (grind $ identity (i-1) <> ket [1] .> sop .> identity (i-1) <> bra [1])
           in
             isIdentity p0 && isIdentity p1

{--------------------------
 Examples
 --------------------------}

-- | A symbolic state |x>
sstate :: Pathsum DMod2
sstate = open $ identity 1

-- | A bell state
bellstate :: Pathsum DMod2
bellstate = fresh <> fresh .> hgate <> (identity 1) .> cxgate

-- | Teleportation circuit
teleport :: Pathsum DMod2
teleport = (identity 1) <> bellstate .>
           cxgate <> (identity 1) .>
           hgate <> cxgate .>
           swapgate <> hgate .>
           (identity 1) <> cxgate .>
           swapgate <> hgate

-- | Teleportation channel
teleportChannel :: Pathsum DMod2
teleportChannel = channelize ((identity 1) <> bellstate) .>
                  channelize (cxgate <> (identity 1)) .>
                  channelize (hgate <> cxgate) .>
                  embed measure 4 (* 3) (* 3) .>
                  embed measure 4 (\i -> i*3 + 1) (\j -> j*3 + 1) .>
                  channelize (swapgate <> hgate) .>
                  channelize ((identity 1) <> cxgate) .>
                  channelize (swapgate <> hgate) .>
                  embed epsilon 4 (* 3) (* 3) .> -- trace out first qubit
                  embed epsilon 2 (* 2) (* 2)    -- trace out second

-- | Verify teleportation
verifyTele :: () -> IO ()
verifyTele _ = case (densify sstate == grind (ptrace . ptrace . densify $ sstate .> teleport)) of
  True -> putStrLn "Identity"
  False -> putStrLn "Not identity"

-- | Verify teleportation channel
verifyTeleC :: () -> IO ()
verifyTeleC _ = case (rho == grind (rho .> teleportChannel)) of
  True -> putStrLn "Identity"
  False -> putStrLn "Not identity"
  where rho = grind $ vectorize $ densify sstate

-- | The |A> = T|+> state
aState :: Pathsum DMod2
aState = fresh .> hgate .> tgate

-- | T gate teleportation channel
teleportTChannel :: Pathsum DMod2
teleportTChannel = channelize ((identity 1) <> aState) .>
                   channelize (cxgate) .>
                   embed measure 2 (\i -> i*2 + 1) (\j -> j*2 + 1) .>
                   channelize (controlled sgate) .>
                   embed epsilon 2 (\i -> i*2 + 1) (\j -> j*2 + 1) -- trace out the resource state

-- | Verify teleportation channel
verifyTeleT :: () -> IO ()
verifyTeleT _ = case (channelize tgate == grind (teleportTChannel)) of
  True -> putStrLn "Identity"
  False -> putStrLn "No identity"

max_stars_repo_name: PariaNaghavi/feynman

max_issues_repo_name: PariaNaghavi/feynman

max_forks_repo_name: PariaNaghavi/feynman


text:
{-# LANGUAGE FlexibleInstances, GeneralizedNewtypeDeriving #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}
-----------------------------------------------------------------------------
-- |
-- Module     : Test.QuickCheck.LinearAlgebra
-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@stanford.edu>
-- Stability  : experimental
--
--
-- Test generators for linear algebra types.
--

module Test.QuickCheck.LinearAlgebra (
    -- Testable element types
    TestElem(..),

    -- * Generating random objects
    -- ** Dimensions
    dim,
    dim2,
    Dim(..),
    Dim2(..),

    -- ** Indices
    index,
    index2,
    Index(..),
    Index2(..),

    -- ** Elements
    elem,
    elems,
    -- realElem,
    -- realElems,

    -- ** Association lists
    assocs,
    assocs2,
    Assocs(..),
    Assocs2(..),
    -- bandedAssocs,

    -- ** Vectors
    vector,
    VectorPair(..),
    VectorTriple(..),
    VectorList(..),
    WeightedVectorList(..),
    NonEmptyVectorList(..),
    NonEmptyWeightedVectorList(..),

    --  ** Matrices
    matrix,
    MatrixPair(..),
    MatrixTriple(..),
    -- hermMatrix,
    -- triMatrix,

    --  Banded matrices
    -- bandwidths,
    -- banded,
    -- bandedWith,
    -- hermBanded,
    -- triBanded,

    -- BandedAssocs(..),

    ) where

import Prelude hiding ( elem )

-- import BLAS.Types( UpLoEnum(..), DiagEnum(..) )
import Control.Monad
import Data.Complex( magnitude )
import Data.Maybe( fromJust )

import Test.QuickCheck hiding ( vector )
import qualified Test.QuickCheck as QC

import Numeric.LinearAlgebra.Types
import Numeric.LinearAlgebra.Vector( Vector )
import qualified Numeric.LinearAlgebra.Vector as V
import Numeric.LinearAlgebra.Matrix( Matrix )
import qualified Numeric.LinearAlgebra.Matrix as M
-- import Data.Matrix.Banded( Banded, maybeBandedFromMatrixStorage )
-- import Data.Matrix.Banded.ST( runSTBanded, unsafeThawBanded,
--     diagViewBanded )

class TestElem e where
    maybeToReal :: e -> Maybe Double
    toReal :: e -> Double
    norm1 :: e -> Double
    norm :: e -> Double
    conjugate :: e -> e

    toReal = fromJust . maybeToReal

instance (TestElem Double) where
    maybeToReal = Just
    toReal = id
    norm1 = abs
    norm = abs
    conjugate = id

instance (TestElem (Complex Double)) where
    maybeToReal (x :+ y) | y == 0    = Just x
                         | otherwise = Nothing
    norm1 (x :+ y) = abs x + abs y
    norm z = magnitude z
    conjugate (x :+ y) = x :+ (-y)

{-
-- | Element types that can be tested with QuickCheck properties.
class (Storable e, Arbitrary e, CoArbitrary e) => TestElem e where
    -- | Inicates whether or not the value should be used in tests.  For
    -- 'Double's, @isTestElem e@ is defined as
    -- @not (isNaN e || isInfinite e || isDenormalized e)@.
    isTestElem :: e -> Bool

instance TestElem Double where
    isTestElem e = not (isNaN e || isInfinite e || isDenormalized e)
    {-# INLINE isTestElem #-}

instance TestElem (Complex Double) where
    isTestElem (x :+ y) = isTestElem x && isTestElem y
    {-# INLINE isTestElem #-}
-}

minit:: Int ->Int
minit  n=  n `mod` 20

-- | Generate a random element.
elem :: (Arbitrary e) => Gen e
elem = arbitrary

-- | Generate a list of elements suitable for testing with.
elems :: (Arbitrary e) => Int -> Gen [e]
elems n = replicateM n elem

-- Generate a random element that has no imaginary part.
-- realElem :: (TestElem e) => Gen e
-- realElem = liftM fromReal elem

-- Generate a list of elements for testing that have no imaginary part.
-- realElems :: (TestElem e) => Int -> Gen [e]
-- realElems n = replicateM n realElem

-- | Get an appropriate dimension for a random vector
dim :: Gen Int
dim = sized $ \s -> do
    (NonNegative n) <- resize (s `div` 4) $ arbitrary
    return (minit n)

-- | Get an appropriate dimension for a random matrix
dim2 :: Gen (Int,Int)
dim2 = sized $ \s -> do
    m <- resize (s `div` 2) $ dim
    n <- resize (s `div` 2) $ dim
    return (m,n)

-- | A vector dimension.
newtype Dim = Dim Int
  deriving (Eq, Ord, Num, Integral, Real, Enum, Show, Read)

instance Arbitrary Dim where
    arbitrary = sized $ \s -> do
        (NonNegative n) <- resize (s `div` 4) $ arbitrary
        return $ Dim n

    shrink (Dim a) =
        [ Dim a' | (NonNegative a') <- shrink (NonNegative a) ]

-- | A matrix dimension.
newtype Dim2 = Dim2 (Int,Int)
  deriving (Eq, Show, Read)

instance Arbitrary Dim2 where
    arbitrary = do
        mn <- dim2
        return $ Dim2 mn

-- | Given a dimension generate a valid index.  The dimension must be positive.
index :: Int -> Gen Int
index n | n <= 0 =
            error $ "index " ++ (show n) ++ ":"
                  ++ " dimension must be positive (QuickCheck error)"
        | otherwise =
            choose (0,n-1)

-- | Given a matrix dimension generate a valid index.
index2 :: (Int,Int) -> Gen (Int,Int)
index2 (m,n) = do
    i <- index m
    j <- index n
    return (i,j)

-- | A dimension and a valid index for it.
data Index = Index Int Int deriving (Eq,Show)
instance Arbitrary Index where
    arbitrary = do
        n <- (1+) `fmap` dim
        i <- index n
        return $ Index n i

-- | A matrix dimension and a valid index for it.
data Index2 = Index2 (Int,Int) (Int,Int) deriving (Eq,Show)
instance Arbitrary Index2 where
    arbitrary = do
        (m',n') <- dim2
        let mn = (m' + 1, n' + 1)
        ij <- index2 mn
        return $ Index2 mn ij

-- | Generate an associations list for a vector of the given dimension.
assocs :: (Arbitrary e) => Int -> Gen [(Int,e)]
assocs n | n == 0    = return []
         | otherwise = do
    l <- choose(0, 2*n)
    is <- replicateM l $ index n
    es <- elems l
    return $ zip is es

-- | Generate an associations list for a matrix of the given shape.
assocs2 :: (Arbitrary e) => (Int,Int) -> Gen [((Int,Int),e)]
assocs2 (m,n) | m*n == 0  = return []
              | otherwise = do
    l <- choose(0, 2*m*n)
    is <- replicateM l $ index2 (m,n)
    es <- elems l
    return $ zip is es

-- | A dimension and an associations list.
data Assocs e = Assocs Int [(Int,e)] deriving (Eq,Show)

instance (Arbitrary e) => Arbitrary (Assocs e) where
    arbitrary = do
        n   <- dim
        ies <- assocs n
        return $ Assocs n ies

    shrink (Assocs n ies) =
        [ Assocs n' $ filter ((< n') . fst) ies
        | n' <- shrink n
        ] ++
        [ Assocs n ies'
        | ies' <- shrink ies
        ]

-- | A shape and an associations list.
data Assocs2 e = Assocs2 (Int,Int) [((Int,Int),e)] deriving (Eq,Show)
instance (Arbitrary e) => Arbitrary (Assocs2 e) where
    arbitrary = do
        mn   <- dim2
        ies  <- assocs2 mn
        return $ Assocs2 mn ies

-- | Generate a random vector of the given size.
vector :: (Arbitrary e, Storable e) => Int -> Gen (Vector e)
vector n = do
    es <- elems n
    return $ V.fromList n es

instance (Arbitrary e, Storable e) => Arbitrary (Vector e) where
    arbitrary = dim >>= vector

    shrink x =
        [ V.slice 0 n x
        | (NonNegative n) <- shrink (NonNegative $ V.dim x)
        ]

-- | Two vectors with the same dimension.
data VectorPair e f =
    VectorPair (Vector e) (Vector f) deriving (Eq, Show)
instance (Arbitrary e, Storable e, Arbitrary f, Storable f) =>
    Arbitrary (VectorPair e f) where
        arbitrary = do
            x <- arbitrary
            y <- vector (V.dim x)
            return $ VectorPair x y

        shrink (VectorPair x y) =
            [ VectorPair (V.slice 0 n' x) (V.slice 0 n' y)
            | n' <- shrink (V.dim x)
            ]

-- | Three vectors with the same dimension.
data VectorTriple e f g =
    VectorTriple (Vector e) (Vector f) (Vector g) deriving (Eq, Show)
instance (Arbitrary e, Storable e, Arbitrary f, Storable f,
          Arbitrary g, Storable g) =>
    Arbitrary (VectorTriple e f g) where
        arbitrary = do
            x <- arbitrary
            y <- vector (V.dim x)
            z <- vector (V.dim x)
            return $ VectorTriple x y z

-- | A nonempty list of vectors with the same dimension.
data NonEmptyVectorList e = NonEmptyVectorList Int [Vector e] deriving (Eq, Show)
instance (Arbitrary e, Storable e) => Arbitrary (NonEmptyVectorList e) where
    arbitrary = do
        x <- arbitrary
        n <- choose (0,20)
        let p = V.dim x
        xs <- replicateM n $ vector p
        return $ NonEmptyVectorList p $ x:xs

-- | A nonempty list of (weight, vector) pairs, with the weights all non-negative
-- and the vectors all having the same dimension.
data NonEmptyWeightedVectorList e = NonEmptyWeightedVectorList Int [(e, Vector e)]
    deriving (Eq, Show)
instance (Arbitrary e, Storable e, Num e) => Arbitrary (NonEmptyWeightedVectorList e) where
    arbitrary = do
        (NonEmptyVectorList p xs) <- arbitrary
        ws <- replicateM (length xs) $ fmap abs arbitrary
        return $ NonEmptyWeightedVectorList p $ zip ws xs

-- | A list of vectors with the same dimension.
data VectorList e = VectorList Int [Vector e] deriving (Eq, Show)
instance (Arbitrary e, Storable e) => Arbitrary (VectorList e) where
    arbitrary = do
        (NonEmptyVectorList p (_:xs)) <- arbitrary
        return $ VectorList p xs

-- | A list of (weight, vector) pairs, with the weights all non-negative
-- and the vectors all having the same dimension.
data WeightedVectorList e = WeightedVectorList Int [(e, Vector e)]
    deriving (Eq, Show)
instance (Arbitrary e, Storable e, Num e) => Arbitrary (WeightedVectorList e) where
    arbitrary = do
        (NonEmptyWeightedVectorList p (_:wxs)) <- arbitrary
        return $ WeightedVectorList p wxs

-- | Generate a random matrix of the given size.
matrix :: (Arbitrary e, Storable e) => (Int,Int) -> Gen (Matrix e)
matrix (m,n) =
    oneof [ raw, sub ]
  where
    raw = do
        es <- elems $ m * n
        return $ M.fromList (m,n) es
    sub = do
        m' <- choose (m, 2*m)
        es <- elems $ m' * n
        return $ M.slice (0,0) (m,n) (M.fromList (m',n) es)

instance (Arbitrary e, Storable e) => Arbitrary (Matrix e) where
    arbitrary = dim2 >>= matrix

-- | Two matrices with the same dimension.
data MatrixPair e f =
    MatrixPair (Matrix e) (Matrix f) deriving (Eq, Show)
instance (Arbitrary e, Storable e, Arbitrary f, Storable f) =>
    Arbitrary (MatrixPair e f) where
        arbitrary = do
            x <- arbitrary
            y <- matrix (M.dim x)
            return $ MatrixPair x y

-- | Three matrices with the same dimension.
data MatrixTriple e f g =
    MatrixTriple (Matrix e) (Matrix f) (Matrix g) deriving (Eq, Show)
instance (Arbitrary e, Storable e, Arbitrary f, Storable f,
          Arbitrary g, Storable g) =>
    Arbitrary (MatrixTriple e f g) where
        arbitrary = do
            x <- arbitrary
            y <- matrix (M.dim x)
            z <- matrix (M.dim x)
            return $ MatrixTriple x y z

instance Arbitrary Trans where
    arbitrary = elements [ NoTrans, Trans, ConjTrans ]
{-

-- | Generate a triangular dense matrix.
triMatrix :: (TestElem e) => (Int,Int) -> Gen (Tri Matrix e)
triMatrix (m,n) = do
    a <- matrix (m,n)
    u <- QC.elements [ Lower, Upper  ]
    d <- QC.elements [ Unit, NonUnit ]
    return $ Tri u d a

-- | Generate a Hermitian dense matrix.
hermMatrix :: (TestElem e) => Int -> Gen (Herm Matrix e)
hermMatrix n = do
    a <- matrix (n,n)
    d <- realElems n
    let a' = runSTMatrix $ do
                 ma <- unsafeThawMatrix a
                 setElems (diagView ma 0) d
                 return ma
    u <- QC.elements [ Lower, Upper ]
    return $ Herm u a'

rawMatrix :: (TestElem e) => (Int,Int) -> Gen (Matrix e)
rawMatrix (m,n) = do
    es <- elems (m*n)
    return $ M.fromList (m,n) es

subMatrix :: (TestElem e) => (Int,Int) -> Gen (SubMatrix e)
subMatrix (m,n) =
    oneof [ rawSubMatrix (m,n)
          , rawSubMatrix (n,m) >>= \(SubMatrix a (i,j) (m',n')) ->
                return $ SubMatrix (herm a) (j,i) (n',m')
          ]

rawSubMatrix :: (TestElem e) => (Int,Int) -> Gen (SubMatrix e)
rawSubMatrix (m,n) = do
    i <- choose (0,5)
    j <- choose (0,5)
    e <- choose (0,5)
    f <- choose (0,5)
    x <- rawMatrix (i+m+e, j+n+f)
    return $ SubMatrix x (i,j) (m,n)

instance (TestElem e) => Arbitrary (SubMatrix e) where
    arbitrary = do
        (m,n) <- shape
        (SubMatrix a ij mn) <- subMatrix (m,n)
        return $ SubMatrix a ij mn

-- | Generate valid bandwidth for a given matrix dimension size
bandwidth :: Int -> Gen Int
bandwidth n = if n == 0 then return 0 else choose (0,n-1)

-- | Generate valid bandwidths for the given matrix shape.
bandwidths :: (Int,Int) -> Gen (Int,Int)
bandwidths (m,n) = liftM2 (,) (bandwidth m) (bandwidth n)

-- | Generate a random banded matrix of the given shape.
banded :: (TestElem e) => (Int,Int) -> Gen (Banded e)
banded mn = do
    lu <- bandwidths mn
    bandedWith lu mn

-- | Generate a random banded matrix with the given bandwidths.
bandedWith :: (TestElem e)
           => (Int,Int) -> (Int,Int) -> Gen (Banded e)
bandedWith lu mn = frequency [ (3, rawBanded mn lu)
                             , (2, hermedBanded mn lu)
                             ]

-- | Generate a triangular banded matrix.
triBanded :: (TestElem e) => Int -> Gen (Tri Banded e)
triBanded n = do
    a <- banded (n,n)
    u <- QC.elements [ Lower, Upper  ]
    d <- QC.elements [ Unit, NonUnit ]
    return $ Tri u d a

-- | Generate a Hermitian banded matrix.
hermBanded :: (TestElem e) => Int -> Gen (Herm Banded e)
hermBanded n = do
    a <- banded (n,n)
    d <- realElems n
    let a' = runSTBanded $ do
                 ma <- unsafeThawBanded a
                 setElems (diagViewBanded ma 0) d
                 return ma
    u <- QC.elements [ Lower, Upper ]
    return $ Herm u a'

rawBanded :: (TestElem e) =>
    (Int,Int) -> (Int,Int) -> Gen (Banded e)
rawBanded (m,n) (kl,ku) =
    let bw = kl+ku+1
    in do
        a <- frequency [ (2, rawMatrix (bw,n))
                       , (1, rawSubMatrix (bw,n) >>= \(SubMatrix b ij _) ->
                                 return $ submatrix b ij (bw,n))
                       ]
        return $ fromJust (maybeBandedFromMatrixStorage (m,n) (kl,ku) a)

hermedBanded :: (TestElem e) =>
    (Int,Int) -> (Int,Int) -> Gen (Banded e)
hermedBanded (m,n) (kl,ku) = do
    x <- rawBanded (n,m) (ku,kl)
    return $ herm x

-- | Generate an associations list for a banded matrix of the given shape
-- and bandwidths.
bandedAssocs :: (TestElem e) => (Int,Int) -> (Int,Int) -> Gen [((Int,Int),e)]
bandedAssocs (m,n) (kl,ku) | m*n == 0  = return []
                           | otherwise = do
    (Nat l) <- arbitrary
    ijs     <- replicateM l $ index2 (kl+1+ku,n)
    let ijs' = mapMaybe (\(i,j) -> let i' = i - j - ku in
                                   if 0 <= i' && i' < m then Just (i',j)
                                                        else Nothing    ) ijs
    es      <- replicateM l elem
    return $ zip ijs' es

-- | A shape, bandwidths, and an associations list.
data BandedAssocs e = BandedAssocs (Int,Int) (Int,Int) [((Int,Int),e)] deriving (Eq,Show)
instance (TestElem e) => Arbitrary (BandedAssocs e) where
    arbitrary = do
        mn  <- shape
        bw  <- bandwidths mn
        ies <- bandedAssocs mn bw
        return $ BandedAssocs mn bw ies
-}

max_stars_repo_name: cartazio/hs-cblas

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas