text:
{-# LANGUAGE FlexibleContexts, FlexibleInstances, TypeFamilies #-}

module Tests.ApproxEq
    (
      ApproxEq(..)
    ) where

import Data.Complex (Complex(..), realPart)
import Data.List (intersperse)
import Data.Maybe (catMaybes)
import Numeric.MathFunctions.Constants (m_epsilon)
import Statistics.Matrix hiding (map, toList)
import Test.QuickCheck
import qualified Data.Vector as V
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U
import qualified Statistics.Matrix as M

class (Eq a, Show a) => ApproxEq a where
    type Bounds a

    eq   :: Bounds a -> a -> a -> Bool
    eql  :: Bounds a -> a -> a -> Property
    eql eps a b = counterexample (show a ++ " /=~ " ++ show b) (eq eps a b)

    (=~)  :: a -> a -> Bool

    (==~) :: a -> a -> Property
    a ==~ b = counterexample (show a ++ " /=~ " ++ show b) (a =~ b)

instance ApproxEq Double where
    type Bounds Double = Double

    eq eps a b
      | a == 0 && b == 0 = True
      | otherwise        = abs (a - b) <= eps * max (abs a) (abs b)
    (=~)  = eq m_epsilon

instance ApproxEq (Complex Double) where
    type Bounds (Complex Double) = Double

    eq eps a@(ar :+ ai) b@(br :+ bi)
      | a == 0 && b == 0 = True
      | otherwise        = abs (ar - br) <= eps * d
                        && abs (ai - bi) <= eps * d
      where
        d = max (realPart $ abs a) (realPart $ abs b)

    (=~)  = eq m_epsilon

instance ApproxEq [Double] where
    type Bounds [Double] = Double

    eq  eps (x:xs) (y:ys) = eq eps x y && eq eps xs ys
    eq  _   []     []     = True
    eq  _   _      _      = False

    eql   = eqll length id id
    (=~)  = eq m_epsilon
    (==~) = eql m_epsilon

instance ApproxEq (U.Vector Double) where
    type Bounds (U.Vector Double) = Double

    eq    = eqv
    (=~)  = eq m_epsilon
    eql   = eqlv
    (==~) = eqlv m_epsilon

instance ApproxEq (V.Vector Double) where
    type Bounds (V.Vector Double) = Double

    eq    = eqv
    (=~)  = eq m_epsilon
    eql   = eqlv
    (==~) = eqlv m_epsilon

instance ApproxEq Matrix where
    type Bounds Matrix = Double

    eq eps (Matrix r1 c1 e1 v1) (Matrix r2 c2 e2 v2) =
      (r1,c1,e1) == (r2,c2,e2) && eq eps v1 v2
    (=~)  = eq m_epsilon
    eql eps a b = eqll dimension M.toList (`quotRem` cols a) eps a b
    (==~) = eql m_epsilon

eqv :: (ApproxEq a, G.Vector v Bool, G.Vector v a) =>
       Bounds a -> v a -> v a -> Bool
eqv eps a b = G.length a == G.length b && G.and (G.zipWith (eq eps) a b)

eqlv :: (ApproxEq [a], G.Vector v a) => Bounds [a] -> v a -> v a -> Property
eqlv eps a b = eql eps (G.toList a) (G.toList b)

eqll :: (ApproxEq l, ApproxEq a, Show c, Show d, Eq d, Bounds l ~ Bounds a) =>
        (l -> d) -> (l -> [a]) -> (Int -> c) -> Bounds l -> l -> l -> Property
eqll dim toList coord eps a b = counterexample fancy $ eq eps a b
  where
    fancy
      | la /= lb  = "size mismatch: " ++ show la ++ " /= " ++ show lb
      | length summary < length full = summary
      | otherwise = full
    summary = concat . intersperse ", " . catMaybes $
              zipWith3 whee (map coord [(0::Int)..]) xs ys
    full | '\n' `elem` sa = sa ++ "  /=~\n" ++ sb
         | otherwise      = sa ++ " /=~" ++ sb
    (sa, sb) = (show a, show b)
    (xs, ys) = (toList a, toList b)
    (la, lb) = (dim a, dim b)
    whee i x y | eq eps x y = Nothing
               | otherwise  = Just $ show i ++ ": " ++ show x ++ " /=~ " ++ show y

max_stars_repo_name: harendra-kumar/statistics

max_issues_repo_name: harendra-kumar/statistics

max_forks_repo_name: harendra-kumar/statistics


text:
module HNNWrapper
  (

  ) where

import AI.HNN.FF.Network
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.HMatrix
import Foreign.Storable
import System.Random
import qualified Data.Vector as V

randomizeMatrixValue :: V.Vector (Matrix a) -> IO (V.Vector (Matrix a))
randomizeMatrixValue = undefined

setMatrixValue :: (Element a) => (Int, Int) -> a -> Matrix a -> Matrix a
setMatrixValue (i,j) a matrix = fromRows newRows
  where
    oldRows = toRows matrix
    oldRow  = oldRows !! i
    newRow  = fromList $ replaceElement j a (toList oldRow)
    newRows = replaceElement i newRow oldRows

(!!+) :: Maybe [a] -> Int -> Maybe a
(!!+) Nothing     _ = Nothing
(!!+) (Just list) i = if i >= length list || i < 0 then Nothing else Just (list !! i)

replaceElement :: Int -> a -> [a] -> [a]
replaceElement i x xs = let (ys,zs) = splitAt i xs
                        in if length zs > 0 then ys ++ [x] ++ (tail zs)
                           else ys ++ [x]

ext: hs

max_stars_repo_name: eklinkhammer/local-dpp

max_issues_repo_name: eklinkhammer/local-dpp

max_forks_repo_name: eklinkhammer/local-dpp


text:
module Grammar where

import Data.Complex
import Data.Matrix

data Command =
  --  let a =              [1, 0]
  InitQ { qName :: String, qVal :: QBit} |
  InitG { gName :: String, gVal :: Gate} |
  Measure QBit |
  Return QBit

data QBit =
  QRef String | QArr [Complex Float] | App Gate QBit

data Gate =
  GRef String | GMatrix (Matrix (Complex Float)) | Tensor Gate Gate | Product Gate Gate

max_stars_repo_name: lukasberglund/quami

max_issues_repo_name: lukasberglund/quami

max_forks_repo_name: lukasberglund/quami


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE TypeSynonymInstances #-}

-- |
-- Module      :  Spiral.Util.Pretty.Maple
-- Copyright   :  (c) 2017-2020 Drexel University
-- License     :  BSD-style
-- Maintainer  :  mainland@drexel.edu

module Spiral.Util.Pretty.Maple (
    Pretty(..)
  ) where

import Data.Ratio
import Data.Complex (Complex(..))
import Data.Monoid ((<>))
import Data.Set (Set)
import qualified Data.Set as Set
import qualified Data.Vector as V
import Text.PrettyPrint.Mainland
import qualified Text.PrettyPrint.Mainland.Class as Pretty

import Spiral.Array
import Spiral.Exp
import Spiral.SPL
import Spiral.Util.Pretty (Assoc(..),
                           Fixity(..),
                           HasFixity(..),
                           addPrec,
                           mulPrec,
                           mulPrec1,
                           infixl_,
                           precOf)

class Pretty a where
    {-# MINIMAL pprPrec | ppr #-}
    ppr     :: a -> Doc
    pprPrec :: Int -> a -> Doc

    pprList     :: [a] -> Doc
    pprPrecList :: Int -> [a] -> Doc

    ppr        = pprPrec 0
    pprPrec _  = ppr

    pprPrecList _  = pprList
    pprList xs     = list (map ppr xs)

instance Pretty a => Pretty [a] where
    ppr     = pprList
    pprPrec = pprPrecList

instance Pretty Bool where
    ppr True  = text "true"
    ppr False = text "false"

instance Pretty Char where
    ppr = char

instance Pretty Int where
    ppr = int

instance Pretty Integer where
    ppr = integer

pprRealFrac :: (RealFrac a, Pretty.Pretty a) => a -> Doc
pprRealFrac x
    | isIntegral = Pretty.ppr (ceiling x :: Integer)
    | otherwise  = Pretty.ppr x
  where
    x' :: Integer
    x' = ceiling x

    isIntegral :: Bool
    isIntegral = fromIntegral x' == x

instance Pretty Float where
    ppr = pprRealFrac

instance Pretty Double where
    ppr = pprRealFrac

instance Pretty a => Pretty (Ratio a)  where
    ppr x = text "Fraction" <> parens (commasep [ppr (numerator x), ppr (denominator x)])

instance Pretty a => Pretty (Set a) where
    ppr xs = enclosesep lbrace rbrace comma (map ppr (Set.toList xs))

instance (Pretty e, IArray r DIM2 e) => Pretty (Matrix r e) where
    ppr m = text "Matrix" <> (parens . ppr . toLists) m

instance Pretty (Const a) where
    pprPrec _ (BoolC x)     = ppr x
    pprPrec _ (IntC x)      = ppr x
    pprPrec _ (IntegerC x)  = ppr x
    pprPrec _ (FloatC x)    = ppr x
    pprPrec _ (DoubleC x)   = ppr x
    pprPrec _ (RationalC x) = ppr x

    pprPrec _ (ComplexC r i)
        | r == 0 && i == 0    = char '0'
        | r == 0 && i == 1    = char 'I'
        | r == 0 && i == (-1) = text "-I"
        | r == 0              = ppr i <> char '*' <> char 'I'
        | i == 0              = ppr r
        | otherwise           = text "Complex" <> pprArgs [ppr r, ppr i]

    pprPrec _ (W _ 0 _) = text "1"
    pprPrec _ (W n 1 _) = text "w_" <> Pretty.ppr n
    pprPrec _ (W n k _) = text "w_" <> Pretty.ppr n <> char '^' <> Pretty.ppr k

    pprPrec p (CycC x) = text (showsPrec p x "")

    pprPrec p (PiC x) = parensIf (p > mulPrec) $
                        Pretty.pprPrec mulPrec1 x <> char '*' <> text "Pi"

    pprPrec p (ModularC x) = Pretty.pprPrec p x

instance Pretty Var where
    ppr = Pretty.ppr

instance Pretty (Exp a) where
    pprPrec p (ConstE c) = pprPrec p c
    pprPrec p (VarE v)   = pprPrec p v

    pprPrec p (UnopE op@Neg e) =
        unop p op e

    pprPrec _ (UnopE op e) =
        Pretty.ppr op <> parens (Pretty.ppr e)

    pprPrec _ (BinopE Quot e1 e2) =
        text "iquo" <> pprArgs [ppr e1, ppr e2]

    pprPrec _ (BinopE Rem e1 e2) =
        text "irem" <> pprArgs [ppr e1, ppr e2]

    pprPrec p (BinopE op e1 e2) =
        infixop p op e1 e2

    pprPrec _ (IdxE ev eis) =
        ppr ev <> (brackets . commasep) [ppr (ix+1) | ix <- eis]

    pprPrec p (ComplexE er ei) =
        parensIf (p > addPrec) $
        pprComplex (er :+ ei)

    pprPrec _ (ReE e) =
        text "Re" <> parens (ppr e)

    pprPrec _ (ImE e) =
        text "Im" <> parens (ppr e)

    pprPrec p (BBinopE op e1 e2) =
        infixop p op e1 e2

    pprPrec _ (IfE e1 e2 e3) =
        text "if" <+> ppr e1 <+>
        text "then" <+> ppr e2 <+>
        text "else" <+> ppr e3 <+>
        text "end if"

instance Pretty Unop where
    ppr Neg    = char '-'
    ppr Abs    = text "abs"
    ppr Signum = text "signum"
    ppr Exp    = text "exp"
    ppr Log    = text "log"
    ppr Sqrt   = text "sqrt"
    ppr Sin    = text "sin"
    ppr Cos    = text "cos"
    ppr Asin   = text "asin"
    ppr Acos   = text "acos"
    ppr Atan   = text "atan"
    ppr Sinh   = text "sinh"
    ppr Cosh   = text "cosh"
    ppr Asinh  = text "asinh"
    ppr Acosh  = text "acosh"
    ppr Atanh  = text "atanh"

instance Pretty Binop where
    ppr Add  = char '+'
    ppr Sub  = char '-'
    ppr Mul  = char '*'
    ppr Quot = text "`quot`"
    ppr Rem  = text "`rem`"
    ppr Div  = text "`div`"
    ppr Mod  = text "`mod`"
    ppr FDiv = char '/'

instance Pretty BBinop where
    ppr Eq = text "="
    ppr Ne = text "<>"
    ppr Lt = text "<"
    ppr Le = text "<="
    ppr Ge = text ">="
    ppr Gt = text ">"

instance (Num e, Pretty e) => Pretty (SPL e) where
    pprPrec _ (I n)      = text "IdentityMatrix" <> parens (ppr n)
    pprPrec _ (Diag xs)  = text "DiagonalMatrix" <> parens (list (map ppr (V.toList xs)))
    pprPrec _ (Kron a b) = text "KroneckerProduct" <> pprArgs [ppr a, ppr b]
    pprPrec _ (DSum a b) = text "DiagonalMatrix" <> parens (list [ppr a, ppr b])
    pprPrec p (Prod a b) = infixop p POp a b
    --pprPrec _ (Prod a b) = text "MatrixMatrixMultiply" <> pprArgs [ppr a, ppr b]
    pprPrec p m          = pprPrec p (toMatrix m)

data MatrixBinop = POp
  deriving (Eq, Ord, Show)

instance HasFixity MatrixBinop where
    fixity POp = infixl_ 7

instance Pretty MatrixBinop where
    ppr POp = char '.'

pprArgs :: [Doc] -> Doc
pprArgs = parens . commasep

pprComplex :: (Eq a, Num a, Pretty a) => Complex a -> Doc
pprComplex (r :+ 0) = ppr r
pprComplex (0 :+ i) = ppr i <> char '*' <> char 'I'
pprComplex (r :+ i) = ppr r <+> char '+' <+> ppr i <> char '*' <> char 'I'

unop :: (Pretty a, Pretty op, HasFixity op)
     => Int -- ^ precedence of context
     -> op  -- ^ operator
     -> a
     -> Doc
unop prec op x =
    parensIf (prec > precOf op) $
    ppr op <> pprPrec (precOf op) x

infixop :: (Pretty a, Pretty b, Pretty op, HasFixity op)
        => Int -- ^ precedence of context
        -> op  -- ^ operator
        -> a   -- ^ left argument
        -> b   -- ^ right argument
        -> Doc
infixop prec op l r =
    parensIf (prec > opPrec) $
    pprPrec leftPrec l <+> ppr op <+/> pprPrec rightPrec r
  where
    leftPrec | opAssoc == RightAssoc = opPrec + 1
             | otherwise             = opPrec

    rightPrec | opAssoc == LeftAssoc = opPrec + 1
              | otherwise            = opPrec

    Fixity opAssoc opPrec = fixity op

max_stars_repo_name: mainland/hspiral

max_issues_repo_name: mainland/hspiral

max_forks_repo_name: mainland/hspiral


text:
{-# LANGUAGE BinaryLiterals #-}
{-# LANGUAGE OverloadedLists #-}

module Main (main) where

-- import Data.Yaml.Aeson
import qualified Data.Aeson as Aeson
import Data.Bits
import Data.Complex
import Data.Ratio ((%))
import Data.Type.Equality
import qualified Data.Vector as V
import qualified Data.Vector.Storable as S
import qualified Data.Vector.Unboxed as U
import Foreign.Storable
import GHC.Exts (IsList (..))
-- import LatticeSymmetries
import LatticeSymmetries.Algebra
import LatticeSymmetries.Basis
import LatticeSymmetries.Benes
-- import qualified LatticeSymmetries.CSR as CSR
import LatticeSymmetries.ComplexRational
import LatticeSymmetries.Dense
import LatticeSymmetries.Generator
import LatticeSymmetries.Group
-- import LatticeSymmetries.IO
import LatticeSymmetries.NonbranchingTerm
import LatticeSymmetries.Operator
import LatticeSymmetries.Parser
-- import LatticeSymmetries.Sparse
-- import LatticeSymmetries.Types
import Test.Hspec
import Text.Parsec (parse)
import Text.PrettyPrint.ANSI.Leijen (hPutDoc, hardline, pretty, putDoc)
import Prelude hiding (Product, Sum, toList)

-- anySpinEDException :: Selector SpinEDException
-- anySpinEDException = const True

-- anyLatticeSymmetriesException :: Selector LatticeSymmetriesException
-- anyLatticeSymmetriesException = const True

-- parseLines :: FromJSON a => [Text] -> Either ParseException a
-- parseLines xs = decodeEither' . encodeUtf8 . unlines $ xs

-- checkRight s y = case s of
--   Left _ -> s `shouldSatisfy` isRight
--   Right x -> x `shouldBe` y

-- checkJust s y = case s of
--   Just x -> x `shouldBe` y
--   Nothing -> s `shouldSatisfy` isJust

type SpinPolynomial =
  Sum (Scaled ComplexRational (Product (Generator Int SpinGeneratorType)))

type FermionicPolynomial =
  Sum (Scaled ComplexRational (Product (Generator Int FermionGeneratorType)))

main :: IO ()
main = hspec $ do
  describe "ToJSON/FromJSON" $ do
    it "parses Permutation" $ do
      Aeson.decode "[3, 0, 2, 1]" `shouldBe` Just (mkPermutation [3, 0, 2, 1])
      Aeson.decode "[3, 0, 2, 1, 4, 5, 6]" `shouldBe` Just (mkPermutation [3, 0, 2, 1, 4, 5, 6])
      Aeson.decode "[]" `shouldBe` Just (mkPermutation [])
      Aeson.decode "{ \"oops\": [2, 1, 0] }" `shouldBe` (Nothing :: Maybe Permutation)
    it "encodes Permutation" $ do
      Aeson.encode (mkPermutation [0, 1, 2, 3]) `shouldBe` "[0,1,2,3]"
      Aeson.encode (mkPermutation []) `shouldBe` "[]"
      Aeson.encode (mkPermutation [3, 5, 0, 1, 2, 4]) `shouldBe` "[3,5,0,1,2,4]"
    it "parses Symmetry" $ do
      Aeson.decode "{ \"permutation\": [1, 2, 3, 0], \"sector\": 2 }" `shouldBe` Just (mkSymmetry [1, 2, 3, 0] 2)
      Aeson.decode "{ \"permutation\": [0, 1], \"sector\": 0 }" `shouldBe` Just (mkSymmetry [0, 1] 0)
      Aeson.decode "{ \"permutation\": [], \"sector\": 0 }" `shouldBe` Just (mkSymmetry [] 0)
    -- it "encodes Symmetry" $ do
    --   Aeson.encode (mkSymmetry [0, 1, 2, 3, 4] 0) `shouldBe` "{\"sector\":0,\"permutation\":[0,1,2,3,4]}"
    --   Aeson.encode (mkSymmetry [] 0) `shouldBe` "{\"sector\":0,\"permutation\":[]}"
    --   Aeson.encode (mkSymmetry [1, 2, 0] 2) `shouldBe` "{\"sector\":2,\"permutation\":[1,2,0]}"
    it "parses Symmetries" $ do
      Aeson.decode "[{\"permutation\": [1, 2, 0], \"sector\": 1}]" `shouldBe` Just (mkSymmetries [mkSymmetry [1, 2, 0] 1])
      Aeson.decode
        "[{\"permutation\": [1, 2, 3, 0], \"sector\": 0}, \
        \ {\"permutation\": [3, 2, 1, 0], \"sector\": 0}]"
        `shouldBe` Just (mkSymmetries [mkSymmetry [1, 2, 3, 0] 0, mkSymmetry [3, 2, 1, 0] 0])
    -- it "encodes Symmetries" $ do
    --   Aeson.encode (mkSymmetries [mkSymmetry [1, 2, 0] 1])
    --     `shouldBe` "[{\"sector\":0,\"permutation\":[0,1,2]},{\"sector\":2,\"permutation\":[1,2,0]},{\"sector\":1,\"permutation\":[2,0,1]}]"
    it "parses/encodes ParticleTy" $ do
      forM_ ([SpinTy, SpinlessFermionTy, SpinfulFermionTy] :: [ParticleTy]) $ \tp ->
        Aeson.decode (Aeson.encode tp) `shouldBe` Just tp
      Aeson.decode "\"spin\"" `shouldBe` Just SpinTy
      Aeson.decode "\"spinless-fermion\"" `shouldBe` Just SpinlessFermionTy
      Aeson.decode "\"spinful-fermion\"" `shouldBe` Just SpinfulFermionTy
    it "parses/encodes Basis 'SpinTy" $ do
      let s₁ = mkSpinBasis 10 (Just 5) Nothing emptySymmetries
          s₂ = mkSpinBasis 4 (Just 2) (Just 1) (mkSymmetries [mkSymmetry [1, 2, 3, 0] 0])
      print $ Aeson.encode s₂
      forM_ ([s₁, s₂] :: [Basis 'SpinTy]) $ \s ->
        Aeson.decode (Aeson.encode s) `shouldBe` Just s
    it "parses/encodes Basis 'SpinfulFermionTy" $ do
      let s₁ = mkSpinfulFermionicBasis 5 (SpinfulTotalParticles 4)
          s₂ = mkSpinfulFermionicBasis 3 (SpinfulPerSector 2 1)
          s₃ = mkSpinfulFermionicBasis 30 SpinfulNoOccupation
      forM_ ([s₁, s₂, s₃] :: [Basis 'SpinfulFermionTy]) $ \s -> do
        print $ Aeson.encode s
        Aeson.decode (Aeson.encode s) `shouldBe` Just s
    it "parses/encodes Basis 'SpinlessFermionTy" $ do
      let s₁ = mkSpinlessFermionicBasis 5 Nothing
          s₂ = mkSpinlessFermionicBasis 10 (Just 5)
      forM_ ([s₁, s₂] :: [Basis 'SpinlessFermionTy]) $ \s -> do
        print $ Aeson.encode s
        Aeson.decode (Aeson.encode s) `shouldBe` Just s

  --   (mkdecode "[{\"permutation\": [1, 2, 0], \"sector\": 1}]" `shouldBe` Just
  --   `shouldBe` Just (mkSymmetries [mkSymmetry [1, 2, 3, 0] 0, mkSymmetry [3, 2, 1, 0] 0])
  -- describe "SymmetrySpec" $ do
  --   it "parses JSON specifications" $ do
  --     let s₁ =
  --           parseLines @SymmetrySpec
  --             [ "permutation: [3, 1, 2, 4]",
  --               "sector: 8"
  --             ]
  --         s₂ =
  --           parseLines @SymmetrySpec
  --             [ "permutation: [1, 1, 2, 4]",
  --               "sector: -5"
  --             ]
  --         s₃ =
  --           parseLines @SymmetrySpec
  --             [ "permutation: []",
  --               "sector: 5"
  --             ]
  --         s₄ =
  --           parseLines @SymmetrySpec
  --             [ "permutation: [1]",
  --               "sector: null"
  --             ]
  --     case s₁ of
  --       Left _ -> s₁ `shouldSatisfy` isRight
  --       Right x -> x `shouldBe` SymmetrySpec (fromList [3, 1, 2, 4]) 8
  --     case s₂ of
  --       Left _ -> s₂ `shouldSatisfy` isRight
  --       Right x -> x `shouldBe` SymmetrySpec (fromList [1, 1, 2, 4]) (-5)
  --     s₃ `shouldSatisfy` isLeft
  --     s₄ `shouldSatisfy` isLeft
  -- describe "denseMatrixFromList" $ do
  --   it "handles empty lists" $ do
  --     denseMatrixFromList []
  --       `shouldBe` Right (DenseMatrix S.empty :: DenseMatrix S.Vector Double)
  --     denseMatrixFromList [[], []]
  --       `shouldBe` Right (DenseMatrix S.empty :: DenseMatrix S.Vector Double)
  --   it "handles square matrices" $ do
  --     denseMatrixFromList [[1, 2], [3, 4]]
  --       `shouldBe` Right (DenseMatrix (fromList [1, 2, 3, 4]) :: DenseMatrix S.Vector Double)
  --   it "handles rectangular matrices" $ do
  --     denseMatrixFromList [[1, 2, 3], [4, 5, 6]]
  --       `shouldBe` Right (DenseMatrix (fromList [1, 2, 3, 4, 5, 6]) :: DenseMatrix S.Vector Double)
  --     denseMatrixFromList [[1, 2], [3, 4], [5, 6]]
  --       `shouldBe` Right (DenseMatrix (fromList [1, 2, 3, 4, 5, 6]) :: DenseMatrix S.Vector 3 2 Double)
  -- describe "DenseMatrixSpec" $ do
  --   it "parses JSON specifications" $ do
  --     parseLines @DenseMatrixSpec
  --       [ "[[1, 2,  1],",
  --         " [4, 2,  1],",
  --         " [1, -2, 3]]"
  --       ]
  --       `checkRight` DenseMatrixSpec [[1, 2, 1], [4, 2, 1], [1, -2, 3]]
  --     parseLines @DenseMatrixSpec
  --       [ "[[1, 2],",
  --         " [4, 2],",
  --         " [1, [-2, 3]]]"
  --       ]
  --       `checkRight` DenseMatrixSpec [[1, 2], [4, 2], [1, (-2) :+ 3]]
  --     parseLines @DenseMatrixSpec ["[]"] `checkRight` DenseMatrixSpec []
  --     parseLines @DenseMatrixSpec ["[[], []]"] `checkRight` DenseMatrixSpec [[], []]
  -- describe "denseToCSRSquare" $ do
  --   it "works" $ do
  --     denseToSparse (fromList [[1, 0, 2], [0, 0, -3]]) `shouldBe` Nothing
  --     print $ denseToSparse (fromList [[1, 0, 2], [0, 4, 0], [0, 0, -3]])
  --     True `shouldBe` True
  --   it ".." $ do
  --     let m =
  --           fromList
  --             [ [1, 0, 0.2, 0],
  --               [0, 1, (-2) :+ 8, 0],
  --               [0, -2.1, 1, 0],
  --               [0 :+ 0.1, 0, 0, 3]
  --             ]
  --         expr = denseToSparse m
  --      in case expr of
  --           Just m' -> sparseToDense m' `shouldBe` m
  --           Nothing -> expr `shouldSatisfy` isJust
  --     let m =
  --           fromList
  --             [[3]]
  --         expr = denseToSparse m
  --      in case expr of
  --           Just m' -> sparseToDense m' `shouldBe` m
  --           Nothing -> expr `shouldSatisfy` isJust
  --     let m =
  --           fromList
  --             []
  --         expr = denseToSparse m
  --      in case expr of
  --           Just m' -> sparseToDense m' `shouldBe` m
  --           Nothing -> expr `shouldSatisfy` isJust
  -- describe "Storable instances" $ do
  --   it "has correct sizeOf" $ do
  --     sizeOf (undefined :: Csparse_matrix) `shouldBe` trueCsparse_matrixSizeOf
  --     sizeOf (undefined :: Cbit_index) `shouldBe` trueCbit_indexSizeOf
  --     sizeOf (undefined :: Cterm) `shouldBe` trueCtermSizeOf
  --     sizeOf (undefined :: Coutput_buffer) `shouldBe` trueCoutput_bufferSizeOf
  --     sizeOf (undefined :: Csparse_operator) `shouldBe` trueCsparse_operatorSizeOf
  --   it "has correct alignment" $ do
  --     alignment (undefined :: Csparse_matrix) `shouldBe` trueCsparse_matrixAlignment
  --     alignment (undefined :: Cbit_index) `shouldBe` trueCbit_indexAlignment
  --     alignment (undefined :: Cterm) `shouldBe` trueCtermAlignment
  --     alignment (undefined :: Coutput_buffer) `shouldBe` trueCoutput_bufferAlignment
  --     alignment (undefined :: Csparse_operator) `shouldBe` trueCsparse_operatorAlignment
  -- describe "applyOperatorTerm'" $ do
  --   it ".." $ do
  --     let spec =
  --           InteractionSpec
  --             [ [1, 0, 0, 0],
  --               [0, -1, 2, 0],
  --               [0, 2, -1, 0],
  --               [0, 0, 0, 1]
  --             ]
  --             [[0, 1]]
  --         (Right term) = toOperatorTerm spec
  --     print =<< term `applyOperatorTerm'` [0x1]
  --     True `shouldBe` True
  -- describe "binarySearch" $ do
  --   it ".." $ do
  --     binarySearch (S.fromList [1, 3, 4, 5, 7, 9]) (4 :: Int) `shouldBe` (Just 2)
  --     binarySearch (S.fromList [1, 3, 4, 5, 7, 9]) (2 :: Int) `shouldBe` Nothing
  --     binarySearch (S.fromList []) (2 :: Int) `shouldBe` Nothing
  --     binarySearch (S.fromList [1]) (1 :: Int) `shouldBe` Just 0
  -- describe "cooToCSR" $ do
  --   it ".." $ do
  --     let (m :: CSR S.Vector Int 3 4 Int) = cooToCsr @U.Vector @S.Vector . fromList $ [(2, 2, -3), (1, 1, 1), (2, 3, 8), (1, 1 :: Int, 3 :: Int)]
  --     print m
  --     m `csrIndex` (1, 1) `shouldBe` 4
  -- describe "denseToCSR" $ do
  --   it ".." $ do
  --     let (m :: CSR S.Vector Int 3 4 Int) = denseToCsr (fromList [[1, 0, 2, 0], [0, 4, 0, 0], [0, 0, -3, 8]])
  --     print m
  --     m `csrIndex` (0, 0) `shouldBe` 1
  --     m `csrIndex` (0, 1) `shouldBe` 0
  --     m `csrIndex` (0, 2) `shouldBe` 2
  --     m `csrIndex` (2, 3) `shouldBe` 8
  -- describe "Num CSR" $ do
  --   it "(+)" $ do
  --     let (a :: CSR S.Vector Int 3 4 Int) =
  --           denseToCsr (fromList [[1, 0, 2, 0], [0, 4, 0, 0], [0, 0, -3, 8]])
  --         (b :: CSR S.Vector Int 3 4 Int) =
  --           denseToCsr (fromList [[3, -8, 0, 0], [1, 5, 0, 0], [0, 0, -3, 0]])
  --         (c :: CSR S.Vector Int 3 4 Int) =
  --           denseToCsr (fromList [[4, -8, 2, 0], [1, 9, 0, 0], [0, 0, -6, 8]])
  --     print a
  --     print b
  --     print c
  --     a + b `shouldBe` c
  --   it "(+)" $ do
  --     let a₁' =
  --           fromList @(SomeCSR S.Vector Int Double) $
  --             [ (0, 0, 0.05052698254161314),
  --               (0, 1, 0.34388186562432865),
  --               (0, 2, 0.3736125370797053),
  --               (0, 4, 0.37422323470438656),
  --               (1, 1, 0.08956717735368347),
  --               (1, 3, 0.20028803087704727),
  --               (2, 0, 0.010520730797154432),
  --               (2, 1, 0.7913535778080222),
  --               (2, 3, 0.6557777513743904),
  --               (2, 4, 0.9609587020099947)
  --             ]
  --         b₁' =
  --           fromList @(SomeCSR S.Vector Int Double) $
  --             [ (2, 3, 0.3942046813633675),
  --               (2, 4, 0.0008746925379253812)
  --             ]
  --         (c₁ :: [(Int, Int, Double)]) =
  --           [ (0, 0, 0.05052698254161314),
  --             (0, 1, 0.34388186562432865),
  --             (0, 2, 0.3736125370797053),
  --             (0, 4, 0.37422323470438656),
  --             (1, 1, 0.08956717735368347),
  --             (1, 3, 0.20028803087704727),
  --             (2, 0, 0.010520730797154432),
  --             (2, 1, 0.7913535778080222),
  --             (2, 3, 1.049982432737758),
  --             (2, 4, 0.96183339454792)
  --           ]
  --         a₂ =
  --           fromList @(CSR S.Vector Int 5 2 Double) $
  --             [(1, 1, 0.9913900269818051), (2, 0, 0.701927418412059)]
  --         b₂ =
  --           fromList @(CSR S.Vector Int 5 2 Double) $
  --             [(2, 0, 0.5563116167541139), (2, 1, 0.5507459794849707), (3, 1, 0.761606773358582)]
  --         (c₂ :: [(Int, Int, Double)]) =
  --           [ (1, 1, 0.9913900269818051),
  --             (2, 0, 1.2582390351661727),
  --             (2, 1, 0.5507459794849707),
  --             (3, 1, 0.761606773358582)
  --           ]
  --     withSomeCsr a₁' $ \a₁ -> withSomeCsr b₁' $ \b₁ ->
  --       case sameShape a₁ b₁ of
  --         Just Refl -> toList (a₁ + b₁) `shouldBe` c₁
  --         x -> x `shouldSatisfy` isJust
  --     toList (a₂ + b₂) `shouldBe` c₂
  -- describe "pPrimitiveOperator" $ do
  --   it ".." $ do
  --     parse pPrimitiveOperator "" ("c†↓₁₀" :: Text)
  --       `shouldBe` Right (SpinfulFermionicOperator FermionicCreationOperator 'c' SpinDown 10)
  --     parse pPrimitiveOperator "" ("f₀" :: Text)
  --       `shouldBe` Right (SpinlessFermionicOperator FermionicAnnihilationOperator 'f' 0)
  --     parse pPrimitiveOperator "" ("n↑₃₈" :: Text)
  --       `shouldBe` Right (SpinfulFermionicOperator FermionicNumberCountingOperator 'n' SpinUp 38)
  --     parse pPrimitiveOperator "" ("Sᶻ₃₈" :: Text)
  --       `shouldBe` Right (SpinOperator SpinZOperator 'S' 38)
  --     parse pOperatorString "" ("n₃₈ f↓₁₅" :: Text)
  --       `shouldBe` Right
  --         ( SpinlessFermionicOperator FermionicNumberCountingOperator 'n' 38
  --             :| [SpinfulFermionicOperator FermionicAnnihilationOperator 'f' SpinDown 15]
  --         )
  --     parse pOperatorString "" ("n↑₃₈f↓₁₅" :: Text)
  --       `shouldBe` Right
  --         ( SpinfulFermionicOperator FermionicNumberCountingOperator 'n' SpinUp 38
  --             :| [SpinfulFermionicOperator FermionicAnnihilationOperator 'f' SpinDown 15]
  --         )
  -- describe "sortByWithParity" $ do
  --   it ".." $ do
  --     sortByWithParity compare [3 :: Int, 1, 4, 8, 0] `shouldBe` (Odd, [0, 1, 3, 4, 8])
  --     sortByWithParity compare [3 :: Int, 1] `shouldBe` (Odd, [1, 3])
  --     sortByWithParity compare [1 :: Int, 3] `shouldBe` (Even, [1, 3])
  --     sortByWithParity compare ([] :: [Int]) `shouldBe` (Even, [])
  -- describe "csrMatMul" $ do
  --   it ".." $ do
  --     let (a :: StorableDenseMatrix 3 4 Float) =
  --           fromList $
  --             [ [0.0, 0.0, 0.0, 0.87383316],
  --               [0.0, 0.0, 0.0, 0.13542224],
  --               [0.24292645, 0.0, 0.0, 0.0]
  --             ]
  --         (b :: StorableDenseMatrix 4 5 Float) =
  --           fromList $
  --             [ [0.0, 0.0, 0.0, 0.67190817, 0.0],
  --               [0.80930562, 0.0, 0.0, 0.94956449, 0.51363394],
  --               [0.63021572, 0.94316889, 0.0, 0.89156391, 0.0],
  --               [0.0, 0.0, 0.64573977, 0.05604065, 0.28140917]
  --             ]
  --         (c :: StorableDenseMatrix 3 5 Float) =
  --           fromList $
  --             [ [0.0, 0.0, 0.56426882, 0.048970178, 0.24590467],
  --               [0.0, 0.0, 0.08744753, 0.00758915, 0.03810906],
  --               [0.0, 0.0, 0.0, 0.16322428, 0.0]
  --             ]
  --     csrMatMul (denseToCsr @S.Vector @Int a) (denseToCsr b) `shouldBe` (denseToCsr c)
  --   it "multiplies matrices" $ do
  --     let (a :: CSR S.Vector Int 5 2 Double) = fromList $ [(0, 0, 0.8213116953652885), (0, 1, 0.2721720739350296), (1, 1, 0.0004725096934642403), (2, 0, 0.03350833939255393), (4, 0, 0.2598645975559749)]
  --         (b :: CSR S.Vector Int 2 6 Double) = fromList $ [(0, 2, 0.7175422685574337), (0, 3, 0.2603509882468078), (0, 5, 0.20233694825026882), (1, 0, 0.9695204042199307), (1, 1, 0.25913885822625227), (1, 3, 0.3241512404858632)]
  --         (c :: CSR S.Vector Int 5 6 Double) = fromList $ [(0, 1, 0.07053036048059469), (0, 0, 0.2638763791388668), (0, 5, 0.16618170200246693), (0, 3, 0.3020542269386639), (0, 2, 0.589325857085161), (1, 3, 0.00015316460327802846), (1, 1, 0.00012244562246515968), (1, 0, 0.0004581077890052858), (2, 5, 0.006779975133623629), (2, 3, 0.008723929275360856), (2, 2, 0.02404364986332557), (4, 5, 0.05258020962776023), (4, 3, 0.06765600478405706), (4, 2, 0.1864638328480788)]
  --     csrMatMul a b `shouldBe` c

  -- describe "combineNeighbors" $ do
  --   it ".." $ do
  --     combineNeighbors (==) (+) (S.fromList [1 :: Int, 2, 3]) `shouldBe` (S.fromList [1, 2, 3])
  --     combineNeighbors (==) (+) (S.fromList [1 :: Int, 1, 3]) `shouldBe` (S.fromList [2, 3])
  --     combineNeighbors (==) (const) (S.fromList [1 :: Int, 1, 2, 1, 1, 1, 3, 3])
  --       `shouldBe` (S.fromList [1, 2, 1, 3])
  --     combineNeighbors (==) (+) (S.fromList ([] :: [Int])) `shouldBe` (S.fromList [])
  describe "pSpinOperator" $ do
    it "parses simple operators" $ do
      parse pSpinOperator "" ("σᶻ₁₀" :: Text)
        `shouldBe` Right (fromList [Scaled 1 (Generator 10 SpinZ)])
      parse pSpinOperator "" ("S⁺₀" :: Text)
        `shouldBe` Right (fromList [Scaled (fromRational (1 % 2)) (Generator 0 SpinPlus)])
    it "parses composite operators" $ do
      parse pSpinOperator "" ("σˣ₅" :: Text)
        `shouldBe` Right
          ( fromList
              [ Scaled 1 (Generator 5 SpinPlus),
                Scaled 1 (Generator 5 SpinMinus)
              ]
          )
      parse pSpinOperator "" ("σʸ₈₉₄₃" :: Text)
        `shouldBe` Right
          ( fromList
              [ Scaled (ComplexRational 0 (-1)) (Generator 8943 SpinPlus),
                Scaled (ComplexRational 0 1) (Generator 8943 SpinMinus)
              ]
          )
  describe "pOperatorString" $ do
    it ".." $ do
      parse (pOperatorString pSpinOperator) "" ("σᶻ₁₀σ⁺₁₁" :: Text)
        `shouldBe` Right [Scaled 1 [Generator 10 SpinZ, Generator 11 SpinPlus]]
      parse (pOperatorString pSpinOperator) "" ("σᶻ₁₀ σ⁺₁₁" :: Text)
        `shouldBe` Right [Scaled 1 [Generator 10 SpinZ, Generator 11 SpinPlus]]
      parse (pOperatorString pSpinOperator) "" ("σᶻ₁₀ Sˣ₅ σ⁺₁₁" :: Text)
        `shouldBe` Right
          ( [ Scaled
                (fromRational (1 % 2))
                [Generator 10 SpinZ, Generator 5 SpinPlus, Generator 11 SpinPlus],
              Scaled
                (fromRational (1 % 2))
                [Generator 10 SpinZ, Generator 5 SpinMinus, Generator 11 SpinPlus]
            ]
          )
  describe "sumToCanonical" $ do
    it ".." $ do
      -- let (Right [x₀]) =
      --       toList <$> parse (pOperatorString pSpinOperator) "" ("σᶻ₁₀ σ⁺₅ σ⁺₁₁" :: Text)
      -- productToCanonical x₀

      let (Right x₀) = parse (pOperatorString pSpinOperator) "" ("σᶻ₁₀ Sˣ₅ σ⁺₁₁" :: Text)
          y₀ :: SpinPolynomial
          y₀ =
            [ Scaled
                (fromRational (1 % 2))
                [Generator 5 SpinPlus, Generator 10 SpinZ, Generator 11 SpinPlus],
              Scaled
                (fromRational (1 % 2))
                [Generator 5 SpinMinus, Generator 10 SpinZ, Generator 11 SpinPlus]
            ]
          x₁ :: SpinPolynomial
          x₁ =
            [ Scaled
                5
                [Generator 5 SpinMinus, Generator 10 SpinZ, Generator 11 SpinPlus],
              Scaled
                3
                [Generator 11 SpinMinus, Generator 5 SpinPlus, Generator 10 SpinZ],
              Scaled
                (-1)
                [Generator 11 SpinMinus, Generator 5 SpinZ]
            ]
          y₁ :: SpinPolynomial
          y₁ =
            [ Scaled
                (-1)
                [Generator 5 SpinZ, Generator 11 SpinMinus],
              Scaled
                3
                [Generator 5 SpinPlus, Generator 10 SpinZ, Generator 11 SpinMinus],
              Scaled
                5
                [Generator 5 SpinMinus, Generator 10 SpinZ, Generator 11 SpinPlus]
            ]
      simplify x₀ `shouldBe` y₀
      simplify x₁ `shouldBe` y₁

  describe "simplify (spin)" $ do
    let parseString ::
          Text ->
          Polynomial ComplexRational (Generator Int SpinGeneratorType)
        -- Sum (Scaled ComplexRational (Product (Generator Int SpinGeneratorType)))
        parseString s = case parse (pOperatorString pSpinOperator) "" s of
          Left e -> error (show e)
          Right x -> x
    it "computes spin commutators" $ do
      simplify (parseString "σ⁺₁₀ σᶻ₁₀" - parseString "σᶻ₁₀ σ⁺₁₀")
        `shouldBe` [Scaled (-2) [Generator 10 SpinPlus]]
      simplify (parseString "σ⁺₁₀ σ⁻₁₀" - parseString "σ⁻₁₀ σ⁺₁₀")
        `shouldBe` [Scaled 1 [Generator 10 SpinZ]]
      simplify (parseString "σ⁻₁₀ σᶻ₁₀" - parseString "σᶻ₁₀ σ⁻₁₀")
        `shouldBe` [Scaled 2 [Generator 10 SpinMinus]]
    it "simplifies products of primitive operators" $ do
      simplify (parseString "σᶻ₁₀ σᶻ₁₀") `shouldBe` [Scaled 1 [Generator 10 SpinIdentity]]
      simplify (parseString "σ⁺₁₀ σ⁺₁₀") `shouldBe` []
      simplify (parseString "σ⁻₁₀ σ⁻₁₀") `shouldBe` []
      simplify (parseString "σ⁺₁₀ σᶻ₁₀") `shouldBe` [Scaled (-1) [Generator 10 SpinPlus]]
      simplify (parseString "σ⁻₁₀ σᶻ₁₀") `shouldBe` [Scaled 1 [Generator 10 SpinMinus]]
      simplify (parseString "σᶻ₁₀ σ⁺₁₀") `shouldBe` [Scaled 1 [Generator 10 SpinPlus]]
      simplify (parseString "σᶻ₁₀ σ⁻₁₀") `shouldBe` [Scaled (-1) [Generator 10 SpinMinus]]
      simplify (parseString "σ⁺₁₀ σ⁻₁₀")
        `shouldBe` [ Scaled (fromRational (1 % 2)) [Generator 10 SpinIdentity],
                     Scaled (fromRational (1 % 2)) [Generator 10 SpinZ]
                   ]
      simplify (parseString "σ⁻₁₀ σ⁺₁₀")
        `shouldBe` [ Scaled (fromRational (1 % 2)) [Generator 10 SpinIdentity],
                     Scaled (fromRational (-1 % 2)) [Generator 10 SpinZ]
                   ]
    it "simplifies products with identities" $ do
      forM_ ([SpinIdentity, SpinPlus, SpinMinus, SpinZ] :: [SpinGeneratorType]) $ \g -> do
        simplify [Scaled 1 [Generator 5 SpinIdentity, Generator 5 g]]
          `shouldBe` [Scaled (1 :: Rational) [Generator (5 :: Int) g]]
        simplify [Scaled 1 [Generator 5 g, Generator 5 SpinIdentity]]
          `shouldBe` [Scaled (1 :: Rational) [Generator (5 :: Int) g]]

  -- describe "CSR.csrMatrixFromDense" $ do
  --   it "converts dense matrices to sparse form" $ do
  --     CSR.csrMatrixFromDense ([[1, 2], [3, 4]] :: Dense.DenseMatrix S.Vector (Complex Double))
  --       `shouldBe` (CSR.CsrMatrix [0, 1, 2] [1, 0] [2, 3] [1, 4])
  -- describe "Num CSR" $ do
  --   it "(+)" $
  --     do
  --       let a = CSR.csrMatrixFromDense ([[1, 2], [0, 2]] :: Dense.DenseMatrix S.Vector (Complex Double))
  --           b = CSR.csrMatrixFromDense ([[0, 0], [1, 3]] :: Dense.DenseMatrix S.Vector (Complex Double))
  --           c = CSR.csrMatrixFromDense ([[1, 2], [1, 5]] :: Dense.DenseMatrix S.Vector (Complex Double))
  --       (a + b) `shouldBe` c

  -- describe "builds operators" $ do
  --   it ".." $ do
  --     print $ mkSpinOperator "σ⁺₁ σ⁻₂" [[0, 1], [1, 2], [2, 0]]
  --     True `shouldBe` True

  describe "NonbranchingTerm" $ do
    it ".." $ do
      let a = nonbranchingRepresentation (Generator (3 :: Int) FermionCreate)
          b = nonbranchingRepresentation (Generator (4 :: Int) FermionAnnihilate)
          c = nonbranchingRepresentation (Generator (3 :: Int) FermionAnnihilate)
          d = nonbranchingRepresentation (Generator (4 :: Int) FermionIdentity)
      print a
      print b
      case a <> a of
        (NonbranchingTerm v _ _ _ _ _) -> v `shouldBe` 0
      case b <> b of
        (NonbranchingTerm v _ _ _ _ _) -> v `shouldBe` 0
      (a <> c) `shouldBe` (nonbranchingRepresentation (Generator (3 :: Int) FermionCount))
      (d <> d) `shouldBe` d
      (a <> d) `shouldBe` a
      (d <> a) `shouldBe` a
      (b <> d) `shouldBe` b
      (d <> b) `shouldBe` b
      print (a <> b)
      print (b <> a)

  -- describe "binomialCoefficient" $ do
  --   it "computes binomial(n, k)" $ do
  --     binomialCoefficient 10 5 `shouldBe` 252
  --     binomialCoefficient 10 0 `shouldBe` 1
  --     binomialCoefficient 10 24 `shouldBe` 0
  --     binomialCoefficient 4 1 `shouldBe` 4
  describe "stateIndex" $ do
    it "computes indices of representatives" $ do
      -- isStateIndexIdentity (basisHeader (mkSpinBasis 10 Nothing)) `shouldBe` True
      -- isStateIndexIdentity (basisHeader (mkSpinBasis 10 (Just 8))) `shouldBe` False
      -- stateIndex (mkSpinBasis 10 Nothing) "|0000000101⟩" `shouldBe` Just 5
      -- stateIndex (mkSpinBasis 10 (Just 2)) "|0000000101⟩" `shouldBe` Just 1
      stateIndex (mkSpinfulFermionicBasis 3 (SpinfulPerSector 2 1)) "|001101⟩" `shouldBe` Just 1
      stateIndex (mkSpinfulFermionicBasis 3 (SpinfulPerSector 2 1)) "|100101⟩" `shouldBe` Just 7
  describe "BenesNetwork" $ do
    it "should permute bits" $ do
      let p₁ = toBenesNetwork (mkPermutation [0, 1, 2])
      permuteBits p₁ 0b100 `shouldBe` 0b100
      let p₂ = toBenesNetwork (mkPermutation [1, 2, 0])
      permuteBits p₂ 0b100 `shouldBe` 0b010
      permuteBits p₂ 0b101 `shouldBe` 0b110
      let v₃ = mkPermutation [1, 2, 3, 0, 5, 6, 7, 4, 9, 10, 11, 8]
          p₃ = toBenesNetwork v₃
      permuteBits p₃ 0b100111010010 `shouldBe` 0b110011100001
      permuteBits p₃ 0b100111010010 `shouldBe` permuteBits' v₃ 0b100111010010
  describe "mkSymmetries" $ do
    it "should build symmetry groups" $ do
      print $
        mkSymmetries
          [ mkSymmetry (mkPermutation [1, 2, 3, 0]) 0,
            mkSymmetry (mkPermutation [3, 2, 1, 0]) 1
          ]
      True `shouldBe` True

-- describe "applyOperator" $ do
--   it ".." $ do
--     let basis1 = LatticeSymmetries.Basis.SpinBasis 10 Nothing
--         operator1 = operatorFromString basis1 "σ⁺₀" [[5]]
--     -- hPutDoc stdout (pretty (opTerms operator1) <> hardline)
--     applyOperator operator1 "|0000000000⟩" `shouldBe` [(1, "|0000100000⟩")]
--     applyOperator operator1 "|0000100000⟩" `shouldBe` []
--     applyOperator operator1 "|0111011111⟩" `shouldBe` [(1, "|0111111111⟩")]
--     let basis2 = LatticeSymmetries.Basis.SpinBasis 2 Nothing
--         operator2 =
--           operatorFromString basis2 "σˣ₀ σˣ₁" [[0, 1]]
--             + operatorFromString basis2 "σʸ₀ σʸ₁" [[0, 1]]
--             + operatorFromString basis2 "σᶻ₀ σᶻ₁" [[0, 1]]
--     -- hPutDoc stdout (pretty (opTerms operator2) <> hardline)
--     -- print $ getNonbranchingTerms operator2
--     applyOperator operator2 "|00⟩" `shouldBe` [(1, "|00⟩")]
--     applyOperator operator2 "|01⟩" `shouldBe` [(2, "|10⟩"), (-1, "|01⟩")]
--     applyOperator operator2 "|10⟩" `shouldBe` [(2, "|01⟩"), (-1, "|10⟩")]
--     applyOperator operator2 "|11⟩" `shouldBe` [(1, "|11⟩")]

--     let basis3 = LatticeSymmetries.Basis.SpinfulFermionicBasis 4 SpinfulNoOccupation
--         operator3 =
--           operatorFromString basis3 "c†↑₀ c↑₁" [[0, 1], [3, 2]]
--     hPutDoc stdout (pretty (opTerms operator3) <> hardline)
--     applyOperator operator3 "|00000010⟩" `shouldBe` [(1, "|00000001⟩")]
--     applyOperator operator3 "|00000100⟩" `shouldBe` [(-1, "|00001000⟩")]

max_stars_repo_name: twesterhout/lattice-symmetries-haskell

max_issues_repo_name: twesterhout/lattice-symmetries-haskell

max_forks_repo_name: twesterhout/lattice-symmetries-haskell


text:
{-# LANGUAGE FlexibleContexts #-}
module Workflow.Toolset
    ( Array1D
    , mean, var, characteristicF, distributionF
    , fft, ifft, distToChar, charToDist, linSpace, asort, azipWith
    ) where

import Data.Array.Unboxed
import Data.Complex
import Data.Bits as Bit
import Data.Array.ST
import Control.Monad.ST
import Data.Array.Base

-- | Array indexed by integers
type Array1D = UArray Int

-- | mean estimation.
-- Complexity O(n).
mean :: (IArray UArray a, Fractional a) => Array1D a -> a
mean arr = s / size where
    (i1,i2) = bounds arr
    size = fromIntegral (i2 - i1 + 1)
    s = afoldl (+) 0 arr

-- | Unbiased variance estimation.
-- Complexity O(n).
var :: (IArray UArray a, Fractional a) => Array1D a -> a
var arr = s / size where
    (i1,i2) = bounds arr
    size = fromIntegral (i2 - i1)
    me = mean arr
    s = afoldl (\acc x -> acc + (x - me)^(2::Int)) 0 arr

-- |
-- Empirical characteristic function of the distribuition represented by array (sample).
-- Complexity O(n).
characteristicF :: (IArray UArray a, RealFloat a) => Array1D a -> [a] -> [Complex a]
characteristicF _ [] = []
characteristicF arr ts = size `seq` ts `seq` map ((/size) . f) ts
  where
    (i1,i2) = bounds arr
    size = fromIntegral (i2 - i1 + 1)
    f t = t `seq` afoldl g 0 arr where
        g acc x = x `seq` acc `seq` acc + (cis $ x*t)

-- |
-- Empirical cumulative distribuition function estimation
-- of the distribuition represented by !sorted! array (sample).
-- Uses linear approximation between points.
-- Complexity O(n).
distributionF :: (IArray UArray a, RealFloat a) => Array1D a -> [a] -> [a]
distributionF _ [] = []
distributionF arr tt@(t:ts) = if t <= (arr ! i1)
        then 0:(distributionF arr ts)
        else go (i1+1) tt
  where
    (i1,i2) = bounds arr
    n = fromIntegral (i2 - i1 + 1)
    go _ [] = []
    go i xt@(x:xs) | i == i2     = take (length xt) . repeat $ 1
                   | arr ! i > x = ( (fromIntegral (i - i1) +
                     (x - (arr ! (i-1))) / ((arr ! i) - (arr ! (i-1)))
                     )/ n):(go i xs)
                   | otherwise   = go (i+1) tt

-- | Converts CDF into Char Func., given grid (x,t)
-- Note, length of t is n/2, and length of x is n
-- resulting array also has length of n/2, but CDF array has length of n
distToChar :: ([Double],[Double]) -> [Double] -> [Complex Double]
distToChar (x,t) f = zipWith (\t' f' -> f' * (0 :+ t'*c) * (cis $ t'*a)) t phi'
  where
    (a,b) = (head x, last x)
    c = (a-b)/(n-1)
    n = fromIntegral . flip Bit.shiftL 1 . length $ t  :: Double
    phi' = reverse . fft True $ zipWith (\f' x' -> (f' - (x'-a)/(b-a)) :+ 0) f x

-- | Converts Char Func. into CDF, given grid (x,t)
-- Note, length of t is n/2, and length of x is n
-- resulting array also has length of n, but Char Func. array has length of n/2
charToDist :: ([Double],[Double]) -> [Complex Double] -> [Double]
charToDist (x,t) phi = postPass
    . shiftUp
    $ zipWith (\x' ps -> (realPart ps) + (x'-a)/(b-a)) x $ ifft False psi
  where
    (a,b) = (head x, last x)
    c = (n-1)/(a-b)
    n' = length $ t  :: Int
    n = fromIntegral . flip Bit.shiftL 1 $ n' :: Double
    phi' = zipWith (\t' ph -> conjugate ph * (0 :+ c/t') * (cis $ t'*a)) t phi
    psi = (:) 0 $ (++) phi' . tail . reverse . map conjugate $ phi'
    shiftUp xxs'@(x':_) = map (\z -> z-x') xxs'
    shiftUp [] = []
    postPass xs = map (min 1 . max 0) xs

-- | Fast Fourier Transformation.
-- Only accepts sizes n = 2^k.
-- Boolean parameter is whether to return full array, or only second (right) part of it.
fft :: RealFloat a => Bool -> [Complex a] -> [Complex a]
fft halve xs = fft' xs w halve
  where
    n = length $ xs
    w = cis $ -2 * pi / fromIntegral n

-- | Inverse Fast Fourier Transformation.
-- Only accepts sizes n = 2^k.
-- Boolean parameter is whether to return full array, or only second (right) part of it.
ifft :: RealFloat a => Bool -> [Complex a] -> [Complex a]
ifft halve xs = map (/ fromIntegral n) $ fft' xs w halve
  where
    n = length $ xs
    w = cis $ 2 * pi / fromIntegral n

-- Fourier transform (both, forward and inverse)
fft' :: Num t => [t] -> t -> Bool -> [t]
fft' []  _ _ = []
fft' [x] _ _ = [x]
fft' xs  w halve = xs `seq` ys `seq` zs `seq` w `seq` halve `seq`
    if halve then zipWith (-) ys ts
             else zipWith (+) ys ts ++ zipWith (-) ys ts
    where wx  = w*w
          ys = fft' evens wx False
          zs = fft' odds wx False
          (evens, odds)  = split xs
          split []       = ([], [])
          split [x]       = ([x], [])
          split (x:y:ss) = (x:xt, y:yt) where (xt, yt) = split ss
          ts = zipWith (*) zs . iterate (w*) $ 1

-- |
-- create x and t grids from bounds on x;
-- Arguments: k power of 2 - number of points: n = 2^k; bounds
linSpace :: Int -> (Double,Double) -> ([Double],[Double])
linSpace k (a,b) = (map (\i' -> a + i' * (b-a) / (n-1)) i, map (\j' -> 2*pi*j'*(n-1)/n/(b-a)) j)
  where
    n = fromIntegral $ (Bit.shiftL 1 k :: Int ):: Double
    i = [0..n-1] :: [Double]
    j = [1..n/2] :: [Double]

-- import Data.Foldable
-- | Left fold on an array
afoldl :: (IArray UArray a) => (b -> a -> b) -> b -> Array1D a -> b
afoldl f x0 arr = f `seq` x0 `seq` arr `seq` go x0 start where
    (start, end) = bounds arr
    g x y = x `seq` y `seq` f x y
    go x i | i > end   = x
           | otherwise = x `seq` go (g x (arr ! i)) (i+1)

---- | Right fold on an array
--afoldr :: (IArray UArray a) => (a -> b -> b) -> b -> Array1D a -> b
--afoldr f x0 arr = f `seq` x0 `seq` arr `seq` go x0 end where
--    (start, end) = bounds arr
--    go x i | i < start = x
--           | otherwise = x `seq` go (f (arr ! i) x) (i-1)

-- | zipWith for arrays
azipWith :: (IArray UArray a, IArray UArray b, IArray UArray c) => (a -> b -> c) -> Array1D a -> Array1D b -> Array1D c
azipWith f a b = a `seq` b `seq` listArray bnds [f (a ! i) (b ! i) | i <- range bnds] where
    bnds = bounds a

-- | working implementation of mergesort; uses the same array of additional memory, slow just as original Data.List.sort
asort :: Array1D Double -> Array1D Double
asort arr = if size < 2 then arr
                        else runSTUArray result
  where
  (start,end) = bounds arr :: (Int, Int)
  size = end-start+1 :: Int
  result = do
    xRight <- (thaw :: Array1D Double -> ST s (STUArray s Int Double)) arr
    xLeft <- (thaw :: Array1D Double -> ST s (STUArray s Int Double)) arr
    let f _ 1 _ = return ()
        f i 2 right = let x = if right then xLeft else xRight
                      in do
          x1 <- unsafeRead x i
          x2 <- unsafeRead x (i+1)
          if x1 > x2 then do
            unsafeWrite x i x2
            unsafeWrite x (i+1) x1
          else return ()
        f i 3 right = let x = if right then xLeft else xRight in do
          x1 <- unsafeRead x i
          x2 <- unsafeRead x (i+1)
          x3 <- unsafeRead x (i+2)
          case (x1 > x2, x1 > x3, x2 > x3) of
            (True,False,False) -> do
                unsafeWrite x i x2
                unsafeWrite x (i+1) x1
            (True,True,True) -> do
                unsafeWrite x i x3
                unsafeWrite x (i+2) x1
            (False,False,True) -> do
                unsafeWrite x (i+1) x3
                unsafeWrite x (i+2) x2
            (False,True,True) -> do
                unsafeWrite x i x3
                unsafeWrite x (i+1) x1
                unsafeWrite x (i+2) x2
            (True,True,False) -> do
                unsafeWrite x i x2
                unsafeWrite x (i+1) x3
                unsafeWrite x (i+2) x1
            _ -> return ()
        f i n right = do
            f i n1 (not right)
            f (i+n1) n2 (not right)
            merge 0 0
          where
            x = if right then xRight else xLeft
            y = if right then xLeft else xRight
            n1 = shiftR n 1
            n2 = n - n1
            i2 = i+n1
            merge k j | k == n1 && j < n2  = sequence_
                          [unsafeRead x (i2+m) >>= unsafeWrite y (i+k+m) | m <- range (j,n2-1)]
                      | k < n1 && j == n2  = sequence_
                          [unsafeRead x (i +m) >>= unsafeWrite y (i+m+j) | m <- range (k,n1-1)]
                      | k < n1 && j < n2 = do
               x1 <- unsafeRead x (i+k)
               x2 <- unsafeRead x (i2+j)
               if x1 > x2 then unsafeWrite y (i+k+j) x2 >> merge k (j+1)
                          else unsafeWrite y (i+k+j) x1 >> merge (k+1) j
                      | otherwise = return ()
    f 0 size True
    return xLeft

max_stars_repo_name: achirkin/haskell-estimatewf

max_issues_repo_name: achirkin/haskell-estimatewf

max_forks_repo_name: achirkin/haskell-estimatewf


text:
module KAB where

import Control.Applicative
import Control.Monad
import Statistics.Distribution (genContVar)
import Statistics.Distribution.Normal (normalDistr)
import Statistics.Distribution.Uniform (uniformDistr)
import System.Random.MWC (withSystemRandom, GenIO)

testComp = "."

normalDist = normalDistr 0 1
uniformDist = uniformDistr 0 1

numActions = 10

runLength = 10000

alpha = 0.1
epsilon = 0.1

valueBias = 0.0

type ActionMeans = [Double]

sampleNormal :: IO Double
sampleNormal = withSystemRandom (genContVar normalDist :: (GenIO -> IO Double))

actionMeans :: IO ActionMeans
actionMeans = replicateM numActions sampleNormal

sampleValue :: Double -> IO Double
sampleValue v = fmap (+ v) sampleNormal

sampleNonstationaryStep :: IO Double
sampleNonstationaryStep = fmap (/ 100) sampleNormal

stepValues :: [Double] -> IO [Double]
stepValues values = sequenceA $ fmap takeStep values
  where
    takeStep val = do
      step <- sampleNonstationaryStep
      return $ val + step

--progress :: [Double] -> Int -> IO (Double, [Double])
--progress values move = (,) <$> (sampleValue $ values !! move) <*> (stepActions values <$> (replicateM numActions sampleNonstationaryStep))

type SampleAverageStats = [(Double, Int)]

type AlphaValEstimates = [Double]

getReward :: ActionMeans -> Int -> IO Double
getReward means actionNum =
  if (length means <= actionNum)
    then error $ "getReward: index " ++ (show actionNum) ++" is too large for list of length " ++ (show $ length means)
    else ((+) (means !! actionNum)) <$> sampleNormal

observeSA :: SampleAverageStats -> Int -> Double -> SampleAverageStats
observeSA oldStats actionNum latestSample =
  if (length oldStats <= actionNum)
    then error $ "action number " ++ (show actionNum) ++ " does not exist in the list " ++ (show oldStats)
    else let (avg, count) = (oldStats !! actionNum) in take actionNum oldStats ++ ((avg + ((latestSample - avg) / (fromIntegral $ count + 1)), count + 1) : (drop (actionNum + 1) oldStats))

observeAlp :: AlphaValEstimates -> Int -> Double -> AlphaValEstimates
observeAlp oldStats actionNum latestSample =
  if (length oldStats <= actionNum)
    then error $ "observeAlp does not contain index " ++ (show actionNum)
    else let oldAvg = (oldStats !! actionNum) in take actionNum oldStats ++ (oldAvg + (alpha * (latestSample - oldAvg)) : (drop (actionNum + 1) oldStats))

argmax :: Ord b => (a -> b) -> [a] -> a
argmax f (a:as) = go a (f a) as
  where
    go amx _ [] = amx
    go amx maxSoFar (b:bs) = if (f b > maxSoFar)
      then go b (f b) bs
      else go amx maxSoFar bs

actGreedySA :: SampleAverageStats -> Int
actGreedySA stats = argmax (saValue stats) [0 .. (length stats - 1)]

actGreedyAlpha :: AlphaValEstimates -> Int
actGreedyAlpha stats = argmax ((!!) stats) [0 .. (length stats - 1)]

actEpsilonGreedy :: Int -> IO Int
actEpsilonGreedy greedyChoice = do
  randomSample <- withSystemRandom (genContVar uniformDist :: (GenIO -> IO Double))
  if (randomSample <= epsilon)
    then randomAction
    else return greedyChoice

randomAction :: IO Int
randomAction = do
  randomSample <- withSystemRandom (genContVar uniformDist :: (GenIO -> IO Double))
  return $ floor $ randomSample * (fromIntegral numActions)

saValue :: SampleAverageStats -> Int -> Double
saValue stats actionNum =
  if (length stats <= actionNum)
    then error "saValue index too large"
    else fst $ stats !! actionNum

initialValsAlp :: AlphaValEstimates
initialValsAlp = take numActions $ repeat valueBias

initialValsSA :: SampleAverageStats
initialValsSA = take numActions $ repeat (valueBias, 0)

initializeScenario :: IO (AlphaValEstimates, SampleAverageStats, ActionMeans)
initializeScenario = (,,) initialValsAlp initialValsSA <$> actionMeans

runScenario :: IO [[Double]]
runScenario = go runLength initialValsAlp initialValsSA =<< actionMeans
  where
    go :: Int -> AlphaValEstimates -> SampleAverageStats -> ActionMeans -> IO [[Double]]
    go 0 _ _ _ = return []
    go n alphaEsts averageEsts means = do
      --putStrLn $ "Time Remaining: " ++ (show n)

      alphaAction <- actEpsilonGreedy $ actGreedyAlpha alphaEsts
      alphaReward <- getReward means alphaAction

      saAction <- actEpsilonGreedy $ actGreedySA averageEsts
      saReward <- getReward means saAction

      newMeans <- stepValues means

      ((:) [alphaReward, saReward]) <$> (go (n - 1) (observeAlp alphaEsts alphaAction alphaReward) (observeSA averageEsts saAction saReward) newMeans)

average :: [[Double]] -> [Double]
average nums = fmap (flip (/) (fromIntegral $ length nums)) $ foldr (zipWith (+)) (head nums) (tail nums)

reportPerformance :: IO [Double]
reportPerformance = fmap average runScenario

f 0 = 0
f n = nextDown + alpha * (1 - nextDown)
  where
    nextDown = f (n - 1)


text:
-- |
-- Module    : Statistics.Test.KruskalWallis
-- Copyright : (c) 2014 Danny Navarro
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
module Statistics.Test.KruskalWallis
  ( -- * Kruskal-Wallis test
    kruskalWallisTest
    -- ** Building blocks
  , kruskalWallisRank
  , kruskalWallis
  , module Statistics.Test.Types
  ) where

import Data.Ord (comparing)
import qualified Data.Vector.Unboxed as U
import Statistics.Function (sort, sortBy, square)
import Statistics.Distribution (complCumulative)
import Statistics.Distribution.ChiSquared (chiSquared)
import Statistics.Types
import Statistics.Test.Types
import Statistics.Test.Internal (rank)
import Statistics.Sample
import qualified Statistics.Sample.Internal as Sample(sum)

-- | Kruskal-Wallis ranking.
--
-- All values are replaced by the absolute rank in the combined samples.
--
-- The samples and values need not to be ordered but the values in the result
-- are ordered. Assigned ranks (ties are given their average rank).
kruskalWallisRank :: (U.Unbox a, Ord a) => [U.Vector a] -> [U.Vector Double]
kruskalWallisRank samples = groupByTags
                          . sortBy (comparing fst)
                          . U.zip tags
                          $ rank (==) joinSample
  where
    (tags,joinSample) = U.unzip
                      . sortBy (comparing snd)
                      $ foldMap (uncurry tagSample) $ zip [(1::Int)..] samples
    tagSample t = U.map (\x -> (t,x))

    groupByTags xs
        | U.null xs = []
        | otherwise = sort (U.map snd ys) : groupByTags zs
      where
        (ys,zs) = U.span ((==) (fst $ U.head xs) . fst) xs

-- | The Kruskal-Wallis Test.
--
-- In textbooks the output value is usually represented by 'K' or 'H'. This
-- function already does the ranking.
kruskalWallis :: (U.Unbox a, Ord a) => [U.Vector a] -> Double
kruskalWallis samples = (nTot - 1) * numerator / denominator
  where
    -- Total number of elements in all samples
    nTot    = fromIntegral $ sumWith rsamples U.length
    -- Average rank of all samples
    avgRank = (nTot + 1) / 2
    --
    numerator = sumWith rsamples $ \sample ->
        let n = fromIntegral $ U.length sample
        in  n * square (mean sample - avgRank)
    denominator = sumWith rsamples $ \sample ->
        Sample.sum $ U.map (\r -> square (r - avgRank)) sample

    rsamples = kruskalWallisRank samples

-- | Perform Kruskal-Wallis Test for the given samples and required
-- significance. For additional information check 'kruskalWallis'. This is just
-- a helper function.
--
-- It uses /Chi-Squared/ distribution for approximation as long as the sizes are
-- larger than 5. Otherwise the test returns 'Nothing'.
kruskalWallisTest :: (Ord a, U.Unbox a) => [U.Vector a] -> Maybe (Test ())
kruskalWallisTest []      = Nothing
kruskalWallisTest samples
  -- We use chi-squared approximation here
  | all (>4) ns = Just Test { testSignificance = mkPValue $ complCumulative d k
                            , testStatistics   = k
                            , testDistribution = ()
                            }
  | otherwise   = Nothing
  where
    k  = kruskalWallis samples
    ns = map U.length samples
    d  = chiSquared (length ns - 1)

-- * Helper functions

sumWith :: Num a => [Sample] -> (Sample -> a) -> a
sumWith samples f = Prelude.sum $ fmap f samples
{-# INLINE sumWith #-}

max_stars_repo_name: vaerksted/statistics

max_issues_repo_name: vaerksted/statistics

max_forks_repo_name: vaerksted/statistics


text:
-----------------------------------------------------------------------------
--
-- Module      :  Drool.Utils.Conversions
-- Copyright   :
-- License     :  AllRightsReserved
--
-- Maintainer  :
-- Stability   :
-- Portability :
--
-- |
--
-----------------------------------------------------------------------------

module Drool.Utils.Conversions (
    gtkColorToGLColor3,
    gtkColorToGLColor4,
    glColor3ToGtkColor,
    glColor4ToGtkColor,
    freqToMs,
    msToFreq,
    floatToComplexDouble,
    floatsToComplexDoubles,
    complexDoubleToFloat,
    complexDoublesToFloats,
    listToCArray,
    listFromCArray,
    interleave,
    interleaveArrays,
    aMax,
    aLength,
    aZip,
    adjustBufferSize,
    adjustBufferSizeBack,
    blendModeSourceFromIndex,
    blendModeFrameBufferFromIndex,
    blendModeSourceIndex,
    blendModeFrameBufferIndex,
) where

import qualified Graphics.UI.Gtk as Gtk
import Graphics.Rendering.OpenGL
import Data.Complex
import Data.Array ( Array )
import Data.Array.IArray ( (!), bounds, listArray, ixmap, amap, rangeSize )
import Data.Array.CArray ( createCArray, elems )
import Foreign.Marshal.Array ( pokeArray, peekArray )
import GHC.Float
import Data.Word ( Word16 )

gtkColorToGLColor3 :: Gtk.Color -> Color3 GLfloat
gtkColorToGLColor3 (Gtk.Color r g b) = Color3 r' g' b'
  where r' = ((fromIntegral r) / 65535.0) :: GLfloat
        g' = ((fromIntegral g) / 65535.0) :: GLfloat
        b' = ((fromIntegral b) / 65535.0) :: GLfloat

gtkColorToGLColor4 :: Gtk.Color -> Word16 -> Color4 GLfloat
gtkColorToGLColor4 (Gtk.Color r g b) a = Color4 r' g' b' a'
  where r' = ((fromIntegral r) / 65535.0) :: GLfloat
        g' = ((fromIntegral g) / 65535.0) :: GLfloat
        b' = ((fromIntegral b) / 65535.0) :: GLfloat
        a' = ((fromIntegral a) / 65535.0) :: GLfloat

glColor3ToGtkColor :: Color3 GLfloat -> Gtk.Color
glColor3ToGtkColor (Color3 r g b) = Gtk.Color r' g' b'
  where r' = round(r * 65535.0)
        g' = round(g * 65535.0)
        b' = round(b * 65535.0)

glColor4ToGtkColor :: Color4 GLfloat -> ( Gtk.Color, Word16 )
glColor4ToGtkColor (Color4 r g b a) = (Gtk.Color r' g' b', a')
  where r' = round(r * 65535.0)
        g' = round(g * 65535.0)
        b' = round(b * 65535.0)
        a' = round(a * 65535.0)

freqToMs :: Int -> Int
freqToMs f = round(1000.0 / fromIntegral f)

msToFreq :: Int -> Int
msToFreq ms = round(1000.0 / fromIntegral ms)

floatToComplexDouble :: Float -> Complex Double
floatToComplexDouble f = (float2Double f :+ 0.0) :: Complex Double

floatsToComplexDoubles :: [Float] -> [Complex Double]
floatsToComplexDoubles fs = (map (\x -> floatToComplexDouble x) fs)

-- Returns magnitude of complex double as float, with magnitude = (re^2 + im^2)^0.5
complexDoubleToFloat :: Complex Double -> Float
complexDoubleToFloat cd = realToFrac(sqrt (realPart cd * realPart cd + imagPart cd * imagPart cd)) :: Float
-- complexDoubleToFloat cd = realToFrac(sqrt (realPart cd * realPart cd)) :: Float

complexDoublesToFloats :: [Complex Double] -> [Float]
complexDoublesToFloats cds = map (\x -> (complexDoubleToFloat x) / n) cds
  where n = fromIntegral $ length cds

listToCArray lst = createCArray (0,(length lst)-1) ( \ptr -> do { pokeArray ptr lst } )

listFromCArray carr = elems carr

interleave :: [a] -> [a] -> [a]
interleave xsA xsB = foldl (\acc (a,b) -> acc ++ [a,b] ) [] (zip xsA xsB)

aLength a = rangeSize $ bounds a

aMax a = snd $ bounds a

interleaveArrays :: Array Int e -> Array Int e -> Array Int e
interleaveArrays a b = listArray (0,(aLength a) + (aMax b)) (concat [ [(a ! i),(b ! i)] | i <- [0..(min (aMax a) (aMax b))] ])

aZip a b = [ ((a ! i),(b ! i)) | i <- [0..(min (aMax a) (aMax b))] ]

adjustBufferSize :: [a] -> Int -> [a]
adjustBufferSize buf maxLen = drop ((length buf)-maxLen) buf -- drop does not alter list for values <= 0

adjustBufferSizeBack :: [a] -> Int -> [a]
adjustBufferSizeBack buf maxLen = take maxLen buf -- drop does not alter list for values <= 0

blendSourceModes :: [ BlendingFactor ]
blendSourceModes = [ Zero, One, DstColor, OneMinusDstColor, SrcAlpha,
                     OneMinusSrcAlpha, DstAlpha, OneMinusDstAlpha,
                     SrcAlphaSaturate ]
blendFrameBufferModes :: [ BlendingFactor ]
blendFrameBufferModes = [ Zero, One, SrcColor, OneMinusSrcColor, SrcAlpha,
                          OneMinusSrcAlpha, DstAlpha, OneMinusDstAlpha ]

blendModeSourceFromIndex :: Int -> BlendingFactor
blendModeSourceFromIndex i = blendSourceModes !! i

blendModeFrameBufferFromIndex :: Int -> BlendingFactor
blendModeFrameBufferFromIndex i = blendFrameBufferModes !! i

blendModeSourceIndex :: BlendingFactor -> Int
blendModeSourceIndex bm = blendModeSourceIndex' bm 0 blendSourceModes
blendModeSourceIndex' :: BlendingFactor -> Int -> [ BlendingFactor ] -> Int
blendModeSourceIndex' bm i (mode:modes) = index
  where index = if bm == mode then i else blendModeSourceIndex' bm (i+1) modes
blendModeSourceIndex' _ _ [] = 0

blendModeFrameBufferIndex :: BlendingFactor -> Int
blendModeFrameBufferIndex bm = blendModeFrameBufferIndex' bm 0 blendFrameBufferModes
blendModeFrameBufferIndex' :: BlendingFactor -> Int -> [ BlendingFactor ] -> Int
blendModeFrameBufferIndex' bm i (mode:modes) = index
  where index = if bm == mode then i else blendModeFrameBufferIndex' bm (i+1) modes
blendModeFrameBufferIndex' _ _ [] = 0

max_stars_repo_name: fuchsto/drool

max_issues_repo_name: fuchsto/drool

max_forks_repo_name: fuchsto/drool


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE Strict           #-}
module FourierPinwheel.GaussianEnvelopePinwheel where

import           Data.Array.Repa             as R
import           Data.Complex
import           Data.Vector.Generic         as VG
import           Data.Vector.Unboxed         as VU
import           FourierPinwheel.Hypergeo1F1
import           Math.Gamma
import           Pinwheel.FourierSeries2D
import           Utils.Distribution
import           Utils.Parallel

-- The Fourier coefficients of e^{-a^2 r^2} r^\mu
{-# INLINE gaussianPinwheelFourierCoefficients #-}
gaussianPinwheelFourierCoefficients ::
     (RealFloat a, Gamma (Complex a), Enum a)
  => Int
  -> a
  -> a
  -> a
  -> Int
  -> Int
  -> a
  -> a
  -> a
  -> Complex a
gaussianPinwheelFourierCoefficients numR2Freqs periodR2 a sigma angularFreq radialFreq periodEnv phi rho =
  let periodConst = (-2) * pi / log periodEnv
      piRhoPConst = pi * rho / periodR2
      real = pi / periodR2 ^ 2 * piRhoPConst ^ abs angularFreq
      mu =
        (2 + sigma + fromIntegral (abs angularFreq)) :+
        (periodConst * fromIntegral radialFreq)
      alpha = mu / 2
      beta = fromIntegral (1 + abs angularFreq) :+ 0
      z = ((-1) * (piRhoPConst / a) ^ 2) :+ 0
      img =
        (0 :+ (-1)) ^ abs angularFreq * cis (fromIntegral (-angularFreq) * phi) *
        gamma alpha /
        gamma beta *
        (a :+ 0) ** (-mu) *
        hypergeom alpha beta z
   in (real :+ 0) * img

gaussianPinwheel ::
     ( RealFloat a
     , Gamma (Complex a)
     , Enum a
     , Unbox a
     , VG.Vector vector (Complex a)
     , NFData (vector (Complex a))
     )
  => Int
  -> a
  -> a
  -> a
  -> Int
  -> Int
  -> a
  -> a
  -> a
  -> vector (Complex a)
gaussianPinwheel numR2Freqs periodR2 stdR2 sigma thetaFreq rFreq periodEnv stdTheta stdR =
  let zeroVec = VG.replicate (numR2Freqs ^ 2) 0
      a = 1 / (stdR2 * sqrt 2)
   in VG.concat .
      parMap
        rdeepseq
        (\(radialFreq, angularFreq) ->
           if angularFreq == 0
             then let pinwheel =
                        centerHollowArray numR2Freqs $
                        createFrequencyArray
                          numR2Freqs
                          (gaussianPinwheelFourierCoefficients
                             numR2Freqs
                             periodR2
                             a
                             sigma
                             angularFreq
                             radialFreq
                             periodEnv)
                      arr =
                        R.map
                          (* ((gaussian1DFourierCoefficients
                                 (fromIntegral radialFreq)
                                 (log periodEnv)
                                 stdR) :+
                              0)) $
                        centerHollowArray numR2Freqs pinwheel
                   in VG.convert . toUnboxed . computeS $ arr
             else zeroVec) $
      [ (radialFreq, angularFreq)
      | radialFreq <- [-rFreq .. rFreq]
      , angularFreq <- [-thetaFreq .. thetaFreq]
      ]

-- This one has a orientation preference at 0 degree. It is used for Koffka cross problem.
gaussianPinwheel1 ::
     ( RealFloat a
     , Gamma (Complex a)
     , Enum a
     , Unbox a
     , VG.Vector vector (Complex a)
     , NFData (vector (Complex a))
     )
  => Int
  -> a
  -> a
  -> a
  -> Int
  -> Int
  -> a
  -> a
  -> a
  -> vector (Complex a)
gaussianPinwheel1 numR2Freqs periodR2 stdR2 sigma thetaFreq rFreq periodEnv stdTheta stdR =
  let a = 1 / (stdR2 * sqrt 2)
   in VG.concat .
      parMap
        rdeepseq
        (\(radialFreq, angularFreq) ->
           let pinwheel =
                 centerHollowArray numR2Freqs $
                 createFrequencyArray
                   numR2Freqs
                   (gaussianPinwheelFourierCoefficients
                      numR2Freqs
                      periodR2
                      a
                      sigma
                      0
                      radialFreq
                      periodEnv)
               arr =
                 R.map
                   (* ((gaussian1DFreq (fromIntegral angularFreq) stdTheta *
                        gaussian1DFourierCoefficients
                          (fromIntegral radialFreq)
                          (log periodEnv)
                          stdR) :+
                       0)) $
                 centerHollowArray numR2Freqs pinwheel
            in VG.convert . toUnboxed . computeS $ arr) $
      [ (radialFreq, angularFreq)
      | radialFreq <- [-rFreq .. rFreq]
      , angularFreq <- [-thetaFreq .. thetaFreq]
      ]

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 BangPatterns #-}
module Main where

import Data.Vector.Storable as V
import Data.Vector.Storable.Mutable as M
import Numeric.FFT.Vector.Unnormalized as U
import Numeric.FFT.Vector.Plan
import Data.Numskell.Vector as N
import Criterion.Main
import Data.Complex

import System.Environment

main = do
    [n] <- fmap (fmap Prelude.read) getArgs
    let numIters = 1000 * 10
    vIn <- M.unsafeNew n
    vOut <- M.unsafeNew n
    N.write vIn =: pure n 17
    let !p = plan dft n
    N.sequence_ $ pure numIters $ executeM p vIn vOut

max_stars_repo_name: TravisWhitaker/vector-fftw

max_issues_repo_name: TravisWhitaker/vector-fftw

max_forks_repo_name: TravisWhitaker/vector-fftw


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE ScopedTypeVariables #-}
#if __GLASGOW_HASKELL__ >= 702
{-# LANGUAGE Trustworthy #-}
#endif

-----------------------------------------------------------------------------
-- |
-- Module      :  Data.Distributive
-- Copyright   :  (C) 2011-2016 Edward Kmett
-- License     :  BSD-style (see the file LICENSE)
--
-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
-- Stability   :  provisional
-- Portability :  portable
--
----------------------------------------------------------------------------
module Data.Distributive
  ( Distributive(..)
  , cotraverse
  , comapM
  , fmapCollect
  ) where

import Control.Applicative
import Control.Applicative.Backwards
import Control.Monad (liftM)
#if __GLASGOW_HASKELL__ < 707
import Control.Monad.Instances ()
#endif
import Control.Monad.Trans.Identity
import Control.Monad.Trans.Reader
import Data.Coerce
import Data.Functor.Compose
import Data.Functor.Identity
import Data.Functor.Product
import Data.Functor.Reverse
import qualified Data.Monoid as Monoid
import Data.Orphans ()

#if MIN_VERSION_base(4,4,0)
import Data.Complex
#endif
#if __GLASGOW_HASKELL__ >= 707 || defined(MIN_VERSION_tagged)
import Data.Proxy
#endif
#if __GLASGOW_HASKELL__ >= 800 || defined(MIN_VERSION_semigroups)
import qualified Data.Semigroup as Semigroup
#endif
#ifdef MIN_VERSION_tagged
import Data.Tagged
#endif
#if __GLASGOW_HASKELL__ >= 702
import GHC.Generics (U1(..), (:*:)(..), (:.:)(..), Par1(..), Rec1(..), M1(..))
#endif

#ifdef HLINT
{-# ANN module "hlint: ignore Use section" #-}
#endif

-- | This is the categorical dual of 'Traversable'.
--
-- Due to the lack of non-trivial comonoids in Haskell, we can restrict
-- ourselves to requiring a 'Functor' rather than
-- some Coapplicative class. Categorically every 'Distributive'
-- functor is actually a right adjoint, and so it must be 'Representable'
-- endofunctor and preserve all limits. This is a fancy way of saying it
-- isomorphic to @(->) x@ for some x.
--
-- To be distributable a container will need to have a way to consistently
-- zip a potentially infinite number of copies of itself. This effectively
-- means that the holes in all values of that type, must have the same
-- cardinality, fixed sized vectors, infinite streams, functions, etc.
-- and no extra information to try to merge together.
--
class Functor g => Distributive g where
#if __GLASGOW_HASKELL__ >= 707
  {-# MINIMAL distribute | collect #-}
#endif
  -- | The dual of 'Data.Traversable.sequenceA'
  --
  -- >>> distribute [(+1),(+2)] 1
  -- [2,3]
  --
  -- @
  -- 'distribute' = 'collect' 'id'
  -- 'distribute' . 'distribute' = 'id'
  -- @
  distribute  :: Functor f => f (g a) -> g (f a)
  distribute  = collect id

  -- |
  -- @
  -- 'collect' f = 'distribute' . 'fmap' f
  -- 'fmap' f = 'runIdentity' . 'collect' ('Identity' . f)
  -- 'fmap' 'distribute' . 'collect' f = 'getCompose' . 'collect' ('Compose' . f)
  -- @

  collect     :: Functor f => (a -> g b) -> f a -> g (f b)
  collect f   = distribute . fmap f

  -- | The dual of 'Data.Traversable.sequence'
  --
  -- @
  -- 'distributeM' = 'fmap' 'unwrapMonad' . 'distribute' . 'WrapMonad'
  -- @
  distributeM :: Monad m => m (g a) -> g (m a)
  distributeM = fmap unwrapMonad . distribute . WrapMonad

  -- |
  -- @
  -- 'collectM' = 'distributeM' . 'liftM' f
  -- @
  collectM    :: Monad m => (a -> g b) -> m a -> g (m b)
  collectM f  = distributeM . liftM f

-- | The dual of 'Data.Traversable.traverse'
--
-- @
-- 'cotraverse' f = 'fmap' f . 'distribute'
-- @
cotraverse :: (Distributive g, Functor f) => (f a -> b) -> f (g a) -> g b
cotraverse f = fmap f . distribute

-- | The dual of 'Data.Traversable.mapM'
--
-- @
-- 'comapM' f = 'fmap' f . 'distributeM'
-- @
comapM :: (Distributive g, Monad m) => (m a -> b) -> m (g a) -> g b
comapM f = fmap f . distributeM

instance Distributive Identity where
  collect = coerce (fmap :: (a -> b) -> f a -> f b)
    :: forall a b f . Functor f => (a -> Identity b) -> f a -> Identity (f b)
  distribute = Identity . fmap runIdentity

#if __GLASGOW_HASKELL__ >= 707 || defined(MIN_VERSION_tagged)
instance Distributive Proxy where
  collect _ _ = Proxy
  distribute _ = Proxy
#endif

#if defined(MIN_VERSION_tagged)
instance Distributive (Tagged t) where
  collect = coerce (fmap :: (a -> b) -> f a -> f b)
    :: forall a b f . Functor f => (a -> Tagged t b) -> f a -> Tagged t (f b)
  distribute = Tagged . fmap unTagged
#endif

instance Distributive ((->)e) where
  distribute a e = fmap ($e) a
  collect f q e = fmap (flip f e) q

instance Distributive g => Distributive (ReaderT e g) where
  distribute a = ReaderT $ \e -> collect (flip runReaderT e) a
  collect f x = ReaderT $ \e -> collect (\a -> runReaderT (f a) e) x

instance Distributive g => Distributive (IdentityT g) where
  collect = coerce (collect :: (a -> g b) -> f a -> g (f b))
            :: forall a b f . Functor f => (a -> IdentityT g b) -> f a -> IdentityT g (f b)

instance (Distributive f, Distributive g) => Distributive (Compose f g) where
  distribute = Compose . fmap distribute . collect getCompose
  collect f = Compose . fmap distribute . collect (coerce f)

instance (Distributive f, Distributive g) => Distributive (Product f g) where
  -- It might be tempting to write a 'collect' implementation that
  -- composes the passed function with fstP and sndP. This could be bad,
  -- because it would lead to the passed function being evaluated twice
  -- for each element of the underlying functor.
  distribute wp = Pair (collect fstP wp) (collect sndP wp) where
    fstP (Pair a _) = a
    sndP (Pair _ b) = b

instance Distributive f => Distributive (Backwards f) where
  distribute = Backwards . collect forwards
  collect = coerce (collect :: (a -> f b) -> g a -> f (g b))
    :: forall g a b . Functor g
    => (a -> Backwards f b) -> g a -> Backwards f (g b)

instance Distributive f => Distributive (Reverse f) where
  distribute = Reverse . collect getReverse
  collect = coerce (collect :: (a -> f b) -> g a -> f (g b))
    :: forall g a b . Functor g
    => (a -> Reverse f b) -> g a -> Reverse f (g b)

instance Distributive Monoid.Dual where
  collect = coerce (fmap :: (a -> b) -> f a -> f b)
    :: forall f a b . Functor f
    => (a -> Monoid.Dual b) -> f a -> Monoid.Dual (f b)
  distribute = Monoid.Dual . fmap Monoid.getDual

instance Distributive Monoid.Product where
  collect = coerce (fmap :: (a -> b) -> f a -> f b)
    :: forall f a b . Functor f
    => (a -> Monoid.Product b) -> f a -> Monoid.Product (f b)
  distribute = Monoid.Product . fmap Monoid.getProduct

instance Distributive Monoid.Sum where
  collect = coerce (fmap :: (a -> b) -> f a -> f b)
    :: forall f a b . Functor f
    => (a -> Monoid.Sum b) -> f a -> Monoid.Sum (f b)
  distribute = Monoid.Sum . fmap Monoid.getSum

#if __GLASGOW_HASKELL__ >= 800 || defined(MIN_VERSION_semigroups)
instance Distributive Semigroup.Min where
  collect = coerce (fmap :: (a -> b) -> f a -> f b)
    :: forall f a b . Functor f
    => (a -> Semigroup.Min b) -> f a -> Semigroup.Min (f b)
  distribute = Semigroup.Min . fmap Semigroup.getMin

instance Distributive Semigroup.Max where
  collect = coerce (fmap :: (a -> b) -> f a -> f b)
    :: forall f a b . Functor f
    => (a -> Semigroup.Max b) -> f a -> Semigroup.Max (f b)
  distribute = Semigroup.Max . fmap Semigroup.getMax

instance Distributive Semigroup.First where
  collect = coerce (fmap :: (a -> b) -> f a -> f b)
    :: forall f a b . Functor f
    => (a -> Semigroup.First b) -> f a -> Semigroup.First (f b)
  distribute = Semigroup.First . fmap Semigroup.getFirst

instance Distributive Semigroup.Last where
  collect = coerce (fmap :: (a -> b) -> f a -> f b)
    :: forall f a b . Functor f
    => (a -> Semigroup.Last b) -> f a -> Semigroup.Last (f b)
  distribute = Semigroup.Last . fmap Semigroup.getLast
#endif

#if MIN_VERSION_base(4,4,0)
instance Distributive Complex where
  distribute wc = fmap realP wc :+ fmap imagP wc where
    -- Redefine realPart and imagPart to avoid incurring redundant RealFloat
    -- constraints on older versions of base
    realP (r :+ _) = r
    imagP (_ :+ i) = i
#endif

-- | 'fmapCollect' is a viable default definition for 'fmap' given
-- a 'Distributive' instance defined in terms of 'collect'.
fmapCollect :: forall f a b . Distributive f => (a -> b) -> f a -> f b
fmapCollect = coerce (collect :: (a -> Identity b) -> f a -> Identity (f b))

#if __GLASGOW_HASKELL__ >= 702
instance Distributive U1 where
  distribute _ = U1

instance (Distributive a, Distributive b) => Distributive (a :*: b) where
  -- It might be tempting to write a 'collect' implementation that
  -- composes the passed function with fstP and sndP. This could be bad,
  -- because it would lead to the passed function being evaluated twice
  -- for each element of the underlying functor.
  distribute f = collect fstP f :*: collect sndP f where
    fstP (l :*: _) = l
    sndP (_ :*: r) = r

instance (Distributive a, Distributive b) => Distributive (a :.: b) where
  distribute = Comp1 . fmap distribute . collect unComp1
  collect f = Comp1 . fmap distribute . collect (coerce f)

instance Distributive Par1 where
  distribute = Par1 . fmap unPar1
  collect = coerce (fmap :: (a -> b) -> f a -> f b)
    :: forall f a b . Functor f => (a -> Par1 b) -> f a -> Par1 (f b)

instance Distributive f => Distributive (Rec1 f) where
  distribute = Rec1 . collect unRec1
  collect = coerce (collect :: (a -> f b) -> g a -> f (g b))
    :: forall g a b . Functor g
    => (a -> Rec1 f b) -> g a -> Rec1 f (g b)

instance Distributive f => Distributive (M1 i c f) where
  distribute = M1 . collect unM1
  collect = coerce (collect :: (a -> f b) -> g a -> f (g b))
    :: forall g a b . Functor g
    => (a -> M1 i c f b) -> g a -> M1 i c f (g b)
#endif

max_stars_repo_name: aaronvargo/distributive

max_issues_repo_name: aaronvargo/distributive

max_forks_repo_name: aaronvargo/distributive


text:
{-# LANGUAGE DeriveFoldable    #-}
{-# LANGUAGE DeriveFunctor     #-}
{-# LANGUAGE DeriveGeneric     #-}
{-# LANGUAGE DeriveTraversable #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE TemplateHaskell   #-}

module Bench where

import           Control.Lens                    (set)
import           Control.Lens.Operators
import           Control.Lens.TH
import           Criterion
import           Criterion.Types                 (Regression (..), Report (..),
                                                  SampleAnalysis (..))
import qualified Data.Binary                     as Bin
import qualified Data.ByteString.Lazy            as B
import           Data.Csv
import qualified Data.Map.Lazy                   as Map
import           Data.Maybe                      (fromJust)
import           Data.Monoid                     ((<>))
import qualified Data.Text                       as T
import qualified Data.Text.Encoding              as T
import           GHC.Generics                    (Generic)
import           Statistics.Resampling.Bootstrap (Estimate (..))

data Bench a = Bench {
      _benchDescription :: T.Text
    , _benchData        :: [(Double, a)]
    }
    deriving (Functor, Generic, Show, Foldable, Traversable)

makeLenses ''Bench

data BenchGroup a = BenchGroup {
      _bgId          :: T.Text
    , _bgDescription :: T.Text
    , _bgBenches     :: [Bench a]
    , _bgXAxisName   :: T.Text
    }
    deriving (Functor, Generic, Show, Foldable, Traversable)

makeLenses ''BenchGroup

instance Bin.Binary a => Bin.Binary (Bench a)
instance Bin.Binary a => Bin.Binary (BenchGroup a)

runBenchGroup :: BenchGroup Benchmarkable -> IO (BenchGroup Report)
runBenchGroup group = do
    putStrLn $ "Benchmarking group: " <> T.unpack (group ^. bgDescription) <>
               " (" <> T.unpack (group ^. bgId) <> ")"
    results <- mapM describeAndRun (group ^. bgBenches)
    return $ set bgBenches results group

describeAndRun :: Bench Benchmarkable -> IO (Bench Report)
describeAndRun bench = do
    putStrLn $ "Running benchmark: " <> T.unpack (bench ^. benchDescription)
    mapM benchAndRetry bench
  where
    benchAndRetry singleBench = do
        result <- benchmark' singleBench
        if getRSquared result < 0.99
            then putStrLn "R^2 < 0.99 => RETRY!" >> benchAndRetry singleBench
            else return result

getOLS :: Report -> Double
getOLS = estPoint . fromJust . Map.lookup "iters" . regCoeffs . head . anRegress . reportAnalysis

getRSquared :: Report -> Double
getRSquared = estPoint . regRSquare . head . anRegress . reportAnalysis

getMean :: Report -> Double
getMean = estPoint . anMean . reportAnalysis

getHeader :: BenchGroup a -> [T.Text]
getHeader group = ["framework", group ^. bgXAxisName, "time"]

getMeasurements :: Bench a -> [Measurement a]
getMeasurements bench =
    map (\(x, y) -> Measurement (bench ^. benchDescription) x y)
        (bench ^. benchData)

data Measurement a = Measurement {
      _mLabel  :: T.Text
    , _mXValue :: Double
    , _mYValue :: a
    }
    deriving (Functor, Generic, Show)

makeLenses ''Measurement

instance ToField a => ToRecord (Measurement a)

toCsv :: ToField a => BenchGroup a -> B.ByteString
toCsv group =
       B.intercalate "," (map (B.fromStrict . T.encodeUtf8) (getHeader group))
    <> "\n"
    <> encode (concatMap getMeasurements (group ^. bgBenches))

max_stars_repo_name: jaaket/dippa-bench

max_issues_repo_name: jaaket/dippa-bench

max_forks_repo_name: jaaket/dippa-bench


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE TupleSections #-}

import Control.Monad
import Criterion
import Criterion.IO
import Criterion.Types
import Data.Aeson
import qualified Data.ByteString.Lazy as B
import qualified Data.HashMap as HM
import Data.List
import qualified Data.Vector as V
import Statistics.Types (confIntLDX, confIntUDX, estError, estPoint)
import System.Environment
import System.Exit
import System.IO
import System.Process as P

toScriptFormat :: [(Word, [Report])] -> Value
toScriptFormat = toJSON . concatMap mapRun
  where
    mapRun (cores, reps) =
        map
            (\rep ->
                 object
                     [ "config" .=
                       object ["experiment" .= reportName rep, "cores" .= cores]
                     , "mean" .=
                       let meanEstimate = anMean $ reportAnalysis rep
                        in object
                               [ "value" .= estPoint meanEstimate
                               , "upper" .= confIntLDX (estError meanEstimate)
                               , "lower" .= confIntUDX (estError meanEstimate)
                               ]
                     ])
            reps

reportFileName = "results/reports.json"

runExperiment lo hi benchmark variants = do
    let fname = benchmark ++ "-" ++ show lo ++ ":" ++ show hi ++ "-results"
        runFor c name = do
            callProcess
                "./bin/benchmarks"
                [ "--json"
                , reportFileName
                , "--match"
                , "prefix"
                , benchmark ++ "/" ++ name
                , "+RTS"
                , "-N" ++ show (c :: Word)
                ]
            (_, _, recs) <- either error pure =<< readJSONReports reportFileName
            pure recs
    let multiVars = filter (/= "sequential") variants
    sequential <- runFor 1 "sequential"
    multiRecs <-
        forM [lo .. hi] $ \c -> do
            recs <- concat <$> mapM (runFor c) multiVars
            pure (c, recs)
    let recs = (1, sequential) : multiRecs
    B.writeFile ("results/" ++ fname ++ ".json") $ encode $ toScriptFormat recs

setUpExperiment = do
    callProcess "stack" ["build"]
    callProcess
        "stack"
        [ "ghc"
        , "--"
        , "src/benchmarks.hs"
        , "-o"
        , "./bin/benchmarks"
        , "-outputdir"
        , "./bin/"
        , "-rtsopts"
        , "-threaded"
        ]

experiments =
    [ ("ohua-bench", ["sequential", "pure", "reading", "writing"])
    , ("comp-bench", ["sequential", "ohua", "par1", "par2"])
    , ("app-bench", ["sequential", "ohua", "par", "ohua-noinline", "par-noinline"])
    , ("cond-bench", ["sequential", "ohua", "par"])
    , ("matmult-bench", ["sequential", "ohua", "par"])
    , ("bs-bench", ["sequential", "ohua", "par"])
    , ("mandel-bench", ["sequential", "ohua", "par"])
    , ("sumeuler-bench", ["sequential", "ohua", "sum_ohua", "par"])
    ]

runAllExperiments lo hi = do
    forM experiments $ \(name, variants) -> runExperiment lo hi name variants

main = do
    [name, cs] <- getArgs
    let (lo, hi) = read cs
    setUpExperiment
    if name == "all"
        then runAllExperiments lo hi
        else (\x -> return [x]) =<<
             (runExperiment lo hi name $ getVariants name)
  where
    getVariants benchmark =
        let variants = HM.lookup benchmark (HM.fromList experiments)
         in case variants of
                (Just v) -> v
                (Nothing) -> error $ "No such benchmark: " ++ benchmark

max_stars_repo_name: ohua-dev/statefulness

max_issues_repo_name: ohua-dev/statefulness

max_forks_repo_name: ohua-dev/statefulness


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE OverloadedLists #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE Strict #-}
{-# LANGUAGE BangPatterns #-}
module SDLib where

import qualified SDL
import qualified Graphics.Rendering.OpenGL.GL as GL
import           Graphics.Rendering.OpenGL.GL (($=))
import           Control.Monad
import qualified Numeric.LinearAlgebra as L
import           Numeric.LinearAlgebra.Data (flatten)
import qualified Numeric.LinearAlgebra.Data as LD
import           Control.Lens      hiding (indices)
import           Debug.Trace
import           Control.Arrow
import qualified Data.Vector as V
import Data.Vector.Strategies
import Control.Concurrent (forkOS, killThread)
import Control.Concurrent.MVar

import           Utils (step, projectionMatrix,normalize,clampSpeed,norm)

import           OBJReader
import           Data.Shaders (loadShaders)

import           Data.Mesh    ( loadTexture
                              , createMesh
                              , flatQuad
                              , quaternion)
import           Data.GameObject
import           Physics.Collision
import           Data.Scene
import           Data.PhysicsData (acc,speed)
import           Data.SkyBox
import           Physics.Mechanics
import           Tasks.Swarm
import           EventHandler
import           Utils.Quaternions

frameTime = round $ 1000/60

-- drawThread ::

debugPrint (GL.DebugMessage _ GL.DebugTypeError _ _ err) = putStrLn err
debugPrint _ = return ()

setup = do
  SDL.swapInterval $= SDL.SynchronizedUpdates
  GL.depthFunc $= Just GL.Less
  -- GL.cullFace $= Just GL.Front
  -- GL.debugMessageCallback $= Just debugPrint
  -- GL.debugOutput $= GL.Enabled

drawThread mv window = forkOS $ do
  context <- SDL.glCreateContext window
  setup
  let uniforms = ["camera","proj","pos","rot","texSampler","normSampler"] :: [String]
  shad <- loadShaders "shaders/vert.glsl" "shaders/frag.glsl" uniforms
  shadNoTex <- loadShaders "shaders/vert.glsl" "shaders/frag_no_text.glsl" uniforms
  ultraMesh <- readOBJ "media/sphere.obj"
  ultraVeh  <- readOBJ "media/cube.obj"
  vehPrim   <- readBoxFromObj "media/cube.obj"
  let miniVeh = ultraMesh
  meh       <- createMesh ultraVeh shadNoTex Nothing Nothing
  veh       <- createMesh ultraMesh shadNoTex Nothing Nothing
  -- crvMesh   <- createMesh curveVertices shadNoTex Nothing
  tex       <- loadTexture "media/bricks_dif.png"
  normalTex <- loadTexture "media/bricks_norm.png"
  sbox      <- loadSkyBox
  qd        <- flatQuad shad tex normalTex
  let quads =
        V.fromList .
        map (\(x,z) -> GameObject qd [x, -2.5, z] Nothing [] (Plane [0,1,0] (-2.5))) $
        [(x,z) | x <- [-30,-26..30], z <- [-30,-26..30]]
      camPos = [ 0, 10,20]
      vehOb pos = initializePhysics (1/10) $
        GameObject meh pos Nothing [] vehPrim
      sphOb pos = initializePhysics (1/10) .
                  GameObject veh pos Nothing [] $
                  Sphere 1 (LD.ident 4)
      dim = 1
      sep = 3
      vehs = -- V.fromList [vehOb [0.5,1,-5]]
        -- V.fromList []
                    -- initializePhysics (1/10) .
                    -- GameObject veh [0.5,1,-5] Nothing [] $
                    -- Sphere 1 (LD.ident 4)]
        V.fromList [sphOb [0,2,sep*z] | z <- [0..dim]]
      mehOb =
        -- initializePhysics (1/100) .
        --             GameObject veh [0, 0, 10] Nothing [] $
        --             Sphere 1 (LD.ident 4)
        initializePhysics (1/10) $
        GameObject meh [0,0,10] Nothing [] vehPrim
      -- curve = GameObject crvMesh [0,0,0] Nothing []
      scene = Scene camPos (Objects mehOb vehs quads) (flatten . L.tr $ projectionMatrix) (Just sbox)
  putMVar mv scene
  drawLoop
  where drawLoop = takeMVar mv >>= draw window >> drawLoop

mainLoop :: MVar Scene -> SDL.Window -> Scene -> IO ()
mainLoop mv window scene = do
  begTime <- SDL.ticks
  events <- SDL.pollEvents
  kState <- SDL.getKeyboardState
  let qpressed   = any handleQuit events
      -- newTime    = time + (0.1*step)
      -- newTime2   = if newTime > 1 then 0 else newTime
      handleKeys = foldr ((.).keyHandler) id $ filter kState keys
      loc        = scene^.objects.player.location
      -- updateNPCS = applyLinear
      --            . over (physics._Just.speed) (clamp 20)
      --            . addGravityForce loc
      !newScene  = -- set objects.player.location) (countPosition newTime2)
                 adjustCamera .
                 over objects resolveObjectCollisions .
                 traverseObjects (floorReaction . applyForce . friction) .
                 traverseObjects gravityF $
                 -- . over (objects.npc) (computeForces updateNPCS)
                 handleKeys scene
  -- print $ newScene^.campos
  endTime <- SDL.ticks
  when ((endTime - begTime) < frameTime)
    $ SDL.delay (frameTime - (endTime - begTime))
  print $ endTime - begTime
  print $ scene^.objects.player
  let cp = prepareCPrim $ scene^.objects.player
      p  = V.toList $ _cpoints cp
      d  = [norm (a- b) | a <- p, b <- p]
  print $ p
  print d
  putMVar mv newScene
  unless qpressed $ mainLoop mv window newScene

sdlMain :: IO ()
sdlMain = do
  SDL.initializeAll
  window <- SDL.createWindow "yolo" SDL.defaultWindow {
      SDL.windowOpenGL = Just SDL.defaultOpenGL
    , SDL.windowInitialSize = SDL.V2 1366 768 }
  mv <- newEmptyMVar
  dthr <- drawThread mv window
  scene <- takeMVar mv
  mainLoop mv window scene
  killThread dthr
  GL.finish
  -- SDL.glDeleteContext context
  SDL.quit

max_stars_repo_name: Antystenes/CPG

max_issues_repo_name: Antystenes/CPG

max_forks_repo_name: Antystenes/CPG