text:
module HashedExpression.Interp
  ( Evaluable (..),
    Approximable (..),
  )
where

import Data.Array
import Data.Complex
import qualified Data.IntMap.Strict as IM
import Data.List (intercalate)
import Data.Map (Map, fromList)
import qualified Data.Map as Map
import Debug.Trace (traceId, traceShowId)
import GHC.TypeLits (KnownNat)
import HashedExpression.Internal.Expression
  ( C,
    ET (..),
    Expression (..),
    ExpressionMap,
    Node (..),
    NodeID,
    NodeID,
    R,
    Scalar,
  )
import HashedExpression.Internal.Node
import HashedExpression.Internal.Utils
import HashedExpression.Prettify (prettify, showExp)
import HashedExpression.Value
import Text.Printf

-- | This operation emulates the mathematical operation
-- | Turn expression to the right type
expZeroR :: ExpressionMap -> NodeID -> Expression Scalar R
expZeroR = flip Expression

expZeroC :: ExpressionMap -> NodeID -> Expression Scalar C
expZeroC = flip Expression

-- | Choose branch base on condition value
chooseBranch :: [Double] -> Double -> [a] -> a
chooseBranch marks val branches
  | val < head marks = head branches
  | otherwise =
    snd . last . filter ((val >=) . fst) $ zip marks (tail branches)

-- | Approximable class
class
  Show a =>
  Approximable a where
  (~=) :: a -> a -> Bool
  prettifyShow :: a -> String

infix 4 ~=

-- |
relativeError :: Double -> Double -> Double
relativeError a b = abs (a - b) / max (abs a) (abs b)

instance Approximable Double where
  (~=) :: Double -> Double -> Bool
  a ~= b
    | abs (a - b) < 1.0e-5 = True
    | a == b = True
    | otherwise = relativeError a b < 0.01
  prettifyShow a
    | abs a < 1e-10 = "0"
    | otherwise = printf "%.2f" a

instance Approximable (Complex Double) where
  (~=) :: Complex Double -> Complex Double -> Bool
  a ~= b = (realPart a ~= realPart b) && (imagPart a ~= imagPart b)
  prettifyShow a =
    prettifyShow (realPart a) ++ " + " ++ prettifyShow (imagPart a) ++ "i"

instance Approximable (Array Int Double) where
  (~=) :: Array Int Double -> Array Int Double -> Bool
  a ~= b = (bounds a == bounds b) && and (zipWith (~=) (elems a) (elems b))
  prettifyShow a =
    "[" ++ (intercalate ", " . map prettifyShow . elems $ a) ++ "]"

instance Approximable (Array Int (Complex Double)) where
  (~=) :: Array Int (Complex Double) -> Array Int (Complex Double) -> Bool
  a ~= b = (bounds a == bounds b) && and (zipWith (~=) (elems a) (elems b))
  prettifyShow a =
    "[" ++ (intercalate ", " . map prettifyShow . elems $ a) ++ "]"

instance Approximable (Array (Int, Int) Double) where
  (~=) :: Array (Int, Int) Double -> Array (Int, Int) Double -> Bool
  a ~= b = (bounds a == bounds b) && and (zipWith (~=) (elems a) (elems b))
  prettifyShow a =
    "[" ++ (intercalate ", " . map prettifyShow . elems $ a) ++ "]"

instance Approximable (Array (Int, Int) (Complex Double)) where
  (~=) ::
    Array (Int, Int) (Complex Double) ->
    Array (Int, Int) (Complex Double) ->
    Bool
  a ~= b = (bounds a == bounds b) && and (zipWith (~=) (elems a) (elems b))
  prettifyShow a =
    "[" ++ (intercalate ", " . map prettifyShow . elems $ a) ++ "]"

instance Approximable (Array (Int, Int, Int) Double) where
  (~=) :: Array (Int, Int, Int) Double -> Array (Int, Int, Int) Double -> Bool
  a ~= b = (bounds a == bounds b) && and (zipWith (~=) (elems a) (elems b))
  prettifyShow a =
    "[" ++ (intercalate ", " . map prettifyShow . elems $ a) ++ "]"

instance Approximable (Array (Int, Int, Int) (Complex Double)) where
  (~=) ::
    Array (Int, Int, Int) (Complex Double) ->
    Array (Int, Int, Int) (Complex Double) ->
    Bool
  a ~= b = (bounds a == bounds b) && and (zipWith (~=) (elems a) (elems b))
  prettifyShow a =
    "[" ++ (intercalate ", " . map prettifyShow . elems $ a) ++ "]"

-- | These should be commented properly.
class Evaluable d rc output | d rc -> output where
  eval :: ValMaps -> Expression d rc -> output

-- |
instance Evaluable Scalar R Double where
  eval :: ValMaps -> Expression Scalar R -> Double
  eval valMap e@(Expression n mp)
    | [] <- retrieveShape n mp =
      case retrieveNode n mp of
        Var name ->
          case Map.lookup name valMap of
            Just (VScalar val) -> val
            _ -> error "no value associated with the variable"
        Const val -> val
        Sum R args -> sum . map (eval valMap . expZeroR mp) $ args
        Mul R args -> product . map (eval valMap . expZeroR mp) $ args
        Neg R arg -> - (eval valMap $ expZeroR mp arg)
        Scale R arg1 arg2 ->
          eval valMap (expZeroR mp arg1)
            * eval valMap (expZeroR mp arg2)
        Power x arg -> eval valMap (expZeroR mp arg) ^ x
        Div arg1 arg2 ->
          eval valMap (expZeroR mp arg1)
            / eval valMap (expZeroR mp arg2)
        Sqrt arg -> sqrt (eval valMap (expZeroR mp arg))
        Sin arg -> sin (eval valMap (expZeroR mp arg))
        Cos arg -> cos (eval valMap (expZeroR mp arg))
        Tan arg -> tan (eval valMap (expZeroR mp arg))
        Exp arg -> exp (eval valMap (expZeroR mp arg))
        Log arg -> log (eval valMap (expZeroR mp arg))
        Sinh arg -> sinh (eval valMap (expZeroR mp arg))
        Cosh arg -> cosh (eval valMap (expZeroR mp arg))
        Tanh arg -> tanh (eval valMap (expZeroR mp arg))
        Asin arg -> asin (eval valMap (expZeroR mp arg))
        Acos arg -> acos (eval valMap (expZeroR mp arg))
        Atan arg -> atan (eval valMap (expZeroR mp arg))
        Asinh arg -> asinh (eval valMap (expZeroR mp arg))
        Acosh arg -> acosh (eval valMap (expZeroR mp arg))
        Atanh arg -> atanh (eval valMap (expZeroR mp arg))
        RealPart arg -> realPart (eval valMap (expZeroC mp arg))
        ImagPart arg -> imagPart (eval valMap (expZeroC mp arg))
        InnerProd R arg1 arg2 ->
          case retrieveShape arg1 mp of
            [] ->
              eval valMap (expZeroR mp arg1)
                * eval valMap (expZeroR mp arg2)
            [size] ->
              let res1 = evaluate1DReal valMap $ (mp, arg1)
                  res2 = evaluate1DReal valMap $ (mp, arg2)
               in sum
                    [ x * y
                      | i <- [0 .. size - 1],
                        let x = res1 ! i,
                        let y = res2 ! i
                    ]
            [size1, size2] ->
              let res1 = evaluate2DReal valMap $ (mp, arg1)
                  res2 = evaluate2DReal valMap $ (mp, arg2)
               in sum
                    [ x * y
                      | i <- [0 .. size1 - 1],
                        j <- [0 .. size2 - 1],
                        let x = res1 ! (i, j),
                        let y = res2 ! (i, j)
                    ]
            [size1, size2, size3] ->
              let res1 = evaluate3DReal valMap $ (mp, arg1)
                  res2 = evaluate3DReal valMap $ (mp, arg2)
               in sum
                    [ x * y
                      | i <- [0 .. size1 - 1],
                        j <- [0 .. size2 - 1],
                        k <- [0 .. size3 - 1],
                        let x = res1 ! (i, j, k),
                        let y = res2 ! (i, j, k)
                    ]
            _ -> error "4D shape?"
        Piecewise marks conditionArg branchArgs ->
          let cdt = eval valMap $ expZeroR mp conditionArg
              branches = map (eval valMap . expZeroR mp) branchArgs
           in chooseBranch marks cdt branches
        _ ->
          error
            ("expression structure Scalar R is wrong " ++ prettify e)
    | otherwise = error "one r but shape is not [] ??"

instance Evaluable Scalar C (Complex Double) where
  eval :: ValMaps -> Expression Scalar C -> Complex Double
  eval valMap e@(Expression n mp)
    | [] <- retrieveShape n mp =
      case retrieveNode n mp of
        Sum C args -> sum . map (eval valMap . expZeroC mp) $ args
        Mul C args -> product . map (eval valMap . expZeroC mp) $ args
        Power x arg -> eval valMap (expZeroC mp arg) ^ x
        Neg C arg -> - (eval valMap $ expZeroC mp arg)
        Scale C arg1 arg2 ->
          case retrieveElementType arg1 mp of
            R ->
              fromR (eval valMap (expZeroR mp arg1))
                * eval valMap (expZeroC mp arg2)
            C ->
              eval valMap (expZeroC mp arg1)
                * eval valMap (expZeroC mp arg2)
        RealImag arg1 arg2 ->
          eval valMap (expZeroR mp arg1)
            :+ eval valMap (expZeroR mp arg2)
        InnerProd C arg1 arg2 ->
          case retrieveShape arg1 mp of
            [] ->
              eval valMap (expZeroC mp arg1)
                * conjugate (eval valMap (expZeroC mp arg2))
            [size] ->
              let res1 = evaluate1DComplex valMap $ (mp, arg1)
                  res2 = evaluate1DComplex valMap $ (mp, arg2)
               in sum
                    [ x * conjugate y
                      | i <- [0 .. size - 1],
                        let x = res1 ! i,
                        let y = res2 ! i
                    ]
            [size1, size2] ->
              let res1 = evaluate2DComplex valMap $ (mp, arg1)
                  res2 = evaluate2DComplex valMap $ (mp, arg2)
               in sum
                    [ x * conjugate y
                      | i <- [0 .. size1 - 1],
                        j <- [0 .. size2 - 1],
                        let x = res1 ! (i, j),
                        let y = res2 ! (i, j)
                    ]
            [size1, size2, size3] ->
              let res1 = evaluate3DComplex valMap $ (mp, arg1)
                  res2 = evaluate3DComplex valMap $ (mp, arg2)
               in sum
                    [ x * conjugate y
                      | i <- [0 .. size1 - 1],
                        j <- [0 .. size2 - 1],
                        k <- [0 .. size3 - 1],
                        let x = res1 ! (i, j, k),
                        let y = res2 ! (i, j, k)
                    ]
            _ -> error "4D shape?"
        Piecewise marks conditionArg branchArgs ->
          let cdt = eval valMap $ expZeroR mp conditionArg
              branches = map (eval valMap . expZeroC mp) branchArgs
           in chooseBranch marks cdt branches
        _ ->
          error
            ("expression structure Scalar C is wrong " ++ prettify e)
    | otherwise = error "One C but shape is not [] ??"

-- |
zipWithA :: Ix x => (a -> b -> c) -> Array x a -> Array x b -> Array x c
zipWithA f xs ys = listArray (bounds xs) $ zipWith f (elems xs) (elems ys)

foldrElementwise :: Ix ix => (a -> a -> a) -> [Array ix a] -> Array ix a
foldrElementwise f [x] = x
foldrElementwise f (x : xs) = zipWithA f x (foldrElementwise f xs)

-- |
evaluate1DReal :: ValMaps -> (ExpressionMap, NodeID) -> Array Int Double
evaluate1DReal valMap (mp, n)
  | [size] <- retrieveShape n mp =
    case retrieveNode n mp of
      Var name ->
        case Map.lookup name valMap of
          Just (V1D val) -> val
          _ -> error "no value associated with the variable"
      Const val -> listArray (0, size - 1) $ replicate size val
      Sum R args ->
        foldrElementwise (+) . map (evaluate1DReal valMap . (mp,)) $
          args
      Mul R args ->
        foldrElementwise (*) . map (evaluate1DReal valMap . (mp,)) $
          args
      Power x arg -> fmap (^ x) (evaluate1DReal valMap $ (mp, arg))
      Neg R arg -> fmap negate . evaluate1DReal valMap $ (mp, arg)
      Scale R arg1 arg2 ->
        let scalar = eval valMap $ expZeroR mp arg1
         in fmap (scalar *) . evaluate1DReal valMap $ (mp, arg2)
      Div arg1 arg2 ->
        zipWithA
          (/)
          (evaluate1DReal valMap $ (mp, arg2))
          (evaluate1DReal valMap $ (mp, arg2))
      Sqrt arg -> fmap sqrt . evaluate1DReal valMap $ (mp, arg)
      Sin arg -> fmap sin . evaluate1DReal valMap $ (mp, arg)
      Cos arg -> fmap cos . evaluate1DReal valMap $ (mp, arg)
      Tan arg -> fmap tan . evaluate1DReal valMap $ (mp, arg)
      Exp arg -> fmap exp . evaluate1DReal valMap $ (mp, arg)
      Log arg -> fmap log . evaluate1DReal valMap $ (mp, arg)
      Sinh arg -> fmap sinh . evaluate1DReal valMap $ (mp, arg)
      Cosh arg -> fmap cosh . evaluate1DReal valMap $ (mp, arg)
      Tanh arg -> fmap tanh . evaluate1DReal valMap $ (mp, arg)
      Asin arg -> fmap asin . evaluate1DReal valMap $ (mp, arg)
      Acos arg -> fmap acos . evaluate1DReal valMap $ (mp, arg)
      Atan arg -> fmap atan . evaluate1DReal valMap $ (mp, arg)
      Asinh arg -> fmap asinh . evaluate1DReal valMap $ (mp, arg)
      Acosh arg -> fmap acosh . evaluate1DReal valMap $ (mp, arg)
      Atanh arg -> fmap atanh . evaluate1DReal valMap $ (mp, arg)
      RealPart arg -> fmap realPart . evaluate1DComplex valMap $ (mp, arg)
      ImagPart arg -> fmap imagPart . evaluate1DComplex valMap $ (mp, arg)
      -- Rotate rA arg ->
      Piecewise marks conditionArg branchArgs ->
        let cdt = evaluate1DReal valMap $ (mp, conditionArg)
            branches = map (evaluate1DReal valMap . (mp,)) branchArgs
         in listArray
              (0, size - 1)
              [ chosen ! i
                | i <- [0 .. size - 1],
                  let chosen = chooseBranch marks (cdt ! i) branches
              ]
      Rotate [amount] arg ->
        rotate1D size amount (evaluate1DReal valMap $ (mp, arg))
      TwiceReFT arg ->
        let innerRes = evaluate1DReal valMap $ (mp, arg)
            scaleFactor = fromIntegral size / 2
         in listArray
              (0, size - 1)
              [ scaleFactor
                  * (innerRes ! i + innerRes ! ((size - i) `mod` size))
                | i <- [0 .. size - 1]
              ]
      TwiceImFT arg ->
        let innerRes = evaluate1DReal valMap $ (mp, arg)
            scaleFactor = fromIntegral size / 2
         in listArray
              (0, size - 1)
              [ scaleFactor
                  * (innerRes ! i - innerRes ! ((size - i) `mod` size))
                | i <- [0 .. size - 1]
              ]
      ReFT arg ->
        case retrieveElementType arg mp of
          R ->
            let inner =
                  fmap (:+ 0) . evaluate1DReal valMap $ (mp, arg)
                ftResult = fourierTransform1D size inner
             in fmap realPart ftResult
          C ->
            let inner = evaluate1DComplex valMap $ (mp, arg)
                ftResult = fourierTransform1D size inner
             in fmap realPart ftResult
      ImFT arg ->
        case retrieveElementType arg mp of
          R ->
            let inner =
                  fmap (:+ 0) . evaluate1DReal valMap $ (mp, arg)
                ftResult = fourierTransform1D size inner
             in fmap imagPart ftResult
          C ->
            let inner = evaluate1DComplex valMap $ (mp, arg)
                ftResult = fourierTransform1D size inner
             in fmap imagPart ftResult
      _ -> error "expression structure One R is wrong"
  | otherwise = error "one r but shape is not [size] ??"

instance (KnownNat n) => Evaluable n R (Array Int Double) where
  eval :: ValMaps -> Expression n R -> Array Int Double
  eval valMap (Expression n mp) = evaluate1DReal valMap (mp, n)

-- |
evaluate1DComplex ::
  ValMaps -> (ExpressionMap, NodeID) -> Array Int (Complex Double)
evaluate1DComplex valMap (mp, n)
  | [size] <- retrieveShape n mp =
    case retrieveNode n mp of
      Sum C args ->
        foldrElementwise (+) . map (evaluate1DComplex valMap . (mp,)) $
          args
      Mul C args ->
        foldrElementwise (*) . map (evaluate1DComplex valMap . (mp,)) $
          args
      Power x arg -> fmap (^ x) (evaluate1DComplex valMap $ (mp, arg))
      Neg C arg -> fmap negate . evaluate1DComplex valMap $ (mp, arg)
      Scale C arg1 arg2 ->
        case retrieveElementType arg1 mp of
          R ->
            let scalar = fromR . eval valMap $ expZeroR mp arg1
             in fmap (scalar *) . evaluate1DComplex valMap $
                  (mp, arg2)
          C ->
            let scalar = eval valMap $ expZeroC mp arg1
             in fmap (scalar *) . evaluate1DComplex valMap $
                  (mp, arg2)
      RealImag arg1 arg2 ->
        zipWithA
          (:+)
          (evaluate1DReal valMap $ (mp, arg1))
          (evaluate1DReal valMap $ (mp, arg2))
      Piecewise marks conditionArg branchArgs ->
        let cdt = evaluate1DReal valMap $ (mp, conditionArg)
            branches =
              map (evaluate1DComplex valMap . (mp,)) branchArgs
         in listArray
              (0, size - 1)
              [ chosen ! i
                | i <- [0 .. size - 1],
                  let chosen = chooseBranch marks (cdt ! i) branches
              ]
      Rotate [amount] arg ->
        rotate1D size amount (evaluate1DComplex valMap $ (mp, arg))
      _ -> error "expression structure One C is wrong"
  | otherwise = error "one C but shape is not [size] ??"

--                         in chooseBranch marks cdt branches
instance (KnownNat n) => Evaluable n C (Array Int (Complex Double)) where
  eval :: ValMaps -> Expression n C -> Array Int (Complex Double)
  eval valMap (Expression n mp) = evaluate1DComplex valMap (mp, n)

-- |
evaluate2DReal :: ValMaps -> (ExpressionMap, NodeID) -> Array (Int, Int) Double
evaluate2DReal valMap (mp, n)
  | [size1, size2] <- retrieveShape n mp =
    case retrieveNode n mp of
      Var name ->
        case Map.lookup name valMap of
          Just (V2D val) -> val
          _ -> error $ "no value associated with the variable" ++ name
      Const val ->
        listArray ((0, 0), (size1 - 1, size2 - 1)) $
          replicate (size1 * size2) val
      Sum R args ->
        foldrElementwise (+) . map (evaluate2DReal valMap . (mp,)) $
          args
      Mul R args ->
        foldrElementwise (*) . map (evaluate2DReal valMap . (mp,)) $
          args
      Power x arg -> fmap (^ x) (evaluate2DReal valMap $ (mp, arg))
      Neg R arg -> fmap negate . evaluate2DReal valMap $ (mp, arg)
      Scale R arg1 arg2 ->
        let scalar = eval valMap $ expZeroR mp arg1
         in fmap (scalar *) . evaluate2DReal valMap $ (mp, arg2)
      Div arg1 arg2 ->
        zipWithA
          (/)
          (evaluate2DReal valMap $ (mp, arg2))
          (evaluate2DReal valMap $ (mp, arg2))
      Sqrt arg -> fmap sqrt . evaluate2DReal valMap $ (mp, arg)
      Sin arg -> fmap sin . evaluate2DReal valMap $ (mp, arg)
      Cos arg -> fmap cos . evaluate2DReal valMap $ (mp, arg)
      Tan arg -> fmap tan . evaluate2DReal valMap $ (mp, arg)
      Exp arg -> fmap exp . evaluate2DReal valMap $ (mp, arg)
      Log arg -> fmap log . evaluate2DReal valMap $ (mp, arg)
      Sinh arg -> fmap sinh . evaluate2DReal valMap $ (mp, arg)
      Cosh arg -> fmap cosh . evaluate2DReal valMap $ (mp, arg)
      Tanh arg -> fmap tanh . evaluate2DReal valMap $ (mp, arg)
      Asin arg -> fmap asin . evaluate2DReal valMap $ (mp, arg)
      Acos arg -> fmap acos . evaluate2DReal valMap $ (mp, arg)
      Atan arg -> fmap atan . evaluate2DReal valMap $ (mp, arg)
      Asinh arg -> fmap asinh . evaluate2DReal valMap $ (mp, arg)
      Acosh arg -> fmap acosh . evaluate2DReal valMap $ (mp, arg)
      Atanh arg -> fmap atanh . evaluate2DReal valMap $ (mp, arg)
      RealPart arg -> fmap realPart . evaluate2DComplex valMap $ (mp, arg)
      ImagPart arg -> fmap imagPart . evaluate2DComplex valMap $ (mp, arg)
      Piecewise marks conditionArg branchArgs ->
        let cdt = evaluate2DReal valMap $ (mp, conditionArg)
            branches = map (evaluate2DReal valMap . (mp,)) branchArgs
         in listArray
              ((0, 0), (size1 - 1, size2 - 1))
              [ chosen ! (i, j)
                | i <- [0 .. size1 - 1],
                  j <- [0 .. size2 - 1],
                  let chosen =
                        chooseBranch marks (cdt ! (i, j)) branches
              ]
      Rotate [amount1, amount2] arg ->
        rotate2D
          (size1, size2)
          (amount1, amount2)
          (evaluate2DReal valMap $ (mp, arg))
      TwiceReFT arg ->
        let innerRes = evaluate2DReal valMap $ (mp, arg)
            scaleFactor = fromIntegral size1 * fromIntegral size2 / 2
         in listArray
              ((0, 0), (size1 - 1, size2 - 1))
              [ scaleFactor
                  * ( innerRes ! (i, j)
                        + innerRes
                        ! ((size1 - i) `mod` size1, (size2 - j) `mod` size2)
                    )
                | i <- [0 .. size1 - 1],
                  j <- [0 .. size2 - 1]
              ]
      TwiceImFT arg ->
        let innerRes = evaluate2DReal valMap $ (mp, arg)
            scaleFactor = fromIntegral size1 * fromIntegral size2 / 2
         in listArray
              ((0, 0), (size1 - 1, size2 - 1))
              [ scaleFactor
                  * ( innerRes ! (i, j)
                        - innerRes
                        ! ((size1 - i) `mod` size1, (size2 - j) `mod` size2)
                    )
                | i <- [0 .. size1 - 1],
                  j <- [0 .. size2 - 1]
              ]
      ReFT arg ->
        case retrieveElementType arg mp of
          R ->
            let inner =
                  fmap (:+ 0) . evaluate2DReal valMap $ (mp, arg)
                ftResult = fourierTransform2D (size1, size2) inner
             in fmap realPart ftResult
          C ->
            let inner = evaluate2DComplex valMap $ (mp, arg)
                ftResult = fourierTransform2D (size1, size2) inner
             in fmap realPart ftResult
      ImFT arg ->
        case retrieveElementType arg mp of
          R ->
            let inner =
                  fmap (:+ 0) . evaluate2DReal valMap $ (mp, arg)
                ftResult = fourierTransform2D (size1, size2) inner
             in fmap imagPart ftResult
          C ->
            let inner = evaluate2DComplex valMap $ (mp, arg)
                ftResult = fourierTransform2D (size1, size2) inner
             in fmap imagPart ftResult
      _ -> error "expression structure Two R is wrong"
  | otherwise = error "Two r but shape is not [size1, size2] ??"

instance
  (KnownNat m, KnownNat n) =>
  Evaluable '(m, n) R (Array (Int, Int) Double)
  where
  eval :: ValMaps -> Expression '(m, n) R -> Array (Int, Int) Double
  eval valMap (Expression n mp) = evaluate2DReal valMap (mp, n)

-- |
evaluate2DComplex ::
  ValMaps -> (ExpressionMap, NodeID) -> Array (Int, Int) (Complex Double)
evaluate2DComplex valMap (mp, n)
  | [size1, size2] <- retrieveShape n mp =
    case retrieveNode n mp of
      Sum C args ->
        foldrElementwise (+) . map (evaluate2DComplex valMap . (mp,)) $
          args
      Mul C args ->
        foldrElementwise (*) . map (evaluate2DComplex valMap . (mp,)) $
          args
      Power x arg -> fmap (^ x) (evaluate2DComplex valMap $ (mp, arg))
      Neg C arg -> fmap negate . evaluate2DComplex valMap $ (mp, arg)
      Scale C arg1 arg2 ->
        case retrieveElementType arg1 mp of
          R ->
            let scalar = fromR . eval valMap $ expZeroR mp arg1
             in fmap (scalar *) . evaluate2DComplex valMap $
                  (mp, arg2)
          C ->
            let scalar = eval valMap $ expZeroC mp arg1
             in fmap (scalar *) . evaluate2DComplex valMap $
                  (mp, arg2)
      RealImag arg1 arg2 ->
        zipWithA
          (:+)
          (evaluate2DReal valMap $ (mp, arg1))
          (evaluate2DReal valMap $ (mp, arg2))
      Piecewise marks conditionArg branchArgs ->
        let cdt = evaluate2DReal valMap $ (mp, conditionArg)
            branches =
              map (evaluate2DComplex valMap . (mp,)) branchArgs
         in listArray
              ((0, 0), (size1 - 1, size2 - 1))
              [ chosen ! (i, j)
                | i <- [0 .. size1 - 1],
                  j <- [0 .. size2 - 1],
                  let chosen =
                        chooseBranch marks (cdt ! (i, j)) branches
              ]
      Rotate [amount1, amount2] arg ->
        rotate2D
          (size1, size2)
          (amount1, amount2)
          (evaluate2DComplex valMap $ (mp, arg))
      _ -> error "expression structure Two C is wrong"
  | otherwise = error "Two C but shape is not [size1, size2] ??"

instance
  (KnownNat m, KnownNat n) =>
  Evaluable '(m, n) C (Array (Int, Int) (Complex Double))
  where
  eval ::
    ValMaps -> Expression '(m, n) C -> Array (Int, Int) (Complex Double)
  eval valMap (Expression n mp) = evaluate2DComplex valMap (mp, n)

evaluate3DReal ::
  ValMaps -> (ExpressionMap, NodeID) -> Array (Int, Int, Int) Double
evaluate3DReal valMap (mp, n)
  | [size1, size2, size3] <- retrieveShape n mp =
    case retrieveNode n mp of
      Var name ->
        case Map.lookup name valMap of
          Just (V3D val) -> val
          _ -> error "no value associated with the variable"
      Const val ->
        listArray ((0, 0, 0), (size1 - 1, size2 - 1, size3 - 1)) $
          replicate (size1 * size2 * size3) val
      Sum R args ->
        foldrElementwise (+) . map (evaluate3DReal valMap . (mp,)) $
          args
      Mul R args ->
        foldrElementwise (*) . map (evaluate3DReal valMap . (mp,)) $
          args
      Power x arg -> fmap (^ x) (evaluate3DReal valMap $ (mp, arg))
      Neg R arg -> fmap negate . evaluate3DReal valMap $ (mp, arg)
      Scale R arg1 arg2 ->
        let scalar = eval valMap $ expZeroR mp arg1
         in fmap (scalar *) . evaluate3DReal valMap $ (mp, arg2)
      Div arg1 arg2 ->
        zipWithA
          (/)
          (evaluate3DReal valMap $ (mp, arg2))
          (evaluate3DReal valMap $ (mp, arg2))
      Sqrt arg -> fmap sqrt . evaluate3DReal valMap $ (mp, arg)
      Sin arg -> fmap sin . evaluate3DReal valMap $ (mp, arg)
      Cos arg -> fmap cos . evaluate3DReal valMap $ (mp, arg)
      Tan arg -> fmap tan . evaluate3DReal valMap $ (mp, arg)
      Exp arg -> fmap exp . evaluate3DReal valMap $ (mp, arg)
      Log arg -> fmap log . evaluate3DReal valMap $ (mp, arg)
      Sinh arg -> fmap sinh . evaluate3DReal valMap $ (mp, arg)
      Cosh arg -> fmap cosh . evaluate3DReal valMap $ (mp, arg)
      Tanh arg -> fmap tanh . evaluate3DReal valMap $ (mp, arg)
      Asin arg -> fmap asin . evaluate3DReal valMap $ (mp, arg)
      Acos arg -> fmap acos . evaluate3DReal valMap $ (mp, arg)
      Atan arg -> fmap atan . evaluate3DReal valMap $ (mp, arg)
      Asinh arg -> fmap asinh . evaluate3DReal valMap $ (mp, arg)
      Acosh arg -> fmap acosh . evaluate3DReal valMap $ (mp, arg)
      Atanh arg -> fmap atanh . evaluate3DReal valMap $ (mp, arg)
      RealPart arg -> fmap realPart . evaluate3DComplex valMap $ (mp, arg)
      ImagPart arg -> fmap imagPart . evaluate3DComplex valMap $ (mp, arg)
      Piecewise marks conditionArg branchArgs ->
        let cdt = evaluate3DReal valMap $ (mp, conditionArg)
            branches = map (evaluate3DReal valMap . (mp,)) branchArgs
         in listArray
              ((0, 0, 0), (size1 - 1, size2 - 1, size3 - 1))
              [ chosen ! (i, j, k)
                | i <- [0 .. size1 - 1],
                  j <- [0 .. size2 - 1],
                  k <- [0 .. size3 - 1],
                  let chosen =
                        chooseBranch marks (cdt ! (i, j, k)) branches
              ]
      Rotate [amount1, amount2, amount3] arg ->
        rotate3D
          (size1, size2, size3)
          (amount1, amount2, amount3)
          (evaluate3DReal valMap $ (mp, arg))
      TwiceReFT arg ->
        let innerRes = evaluate3DReal valMap $ (mp, arg)
            scaleFactor =
              fromIntegral size1 * fromIntegral size2
                * fromIntegral size3
                / 2
         in listArray
              ((0, 0, 0), (size1 - 1, size2 - 1, size3 - 1))
              [ scaleFactor
                  * ( innerRes ! (i, j, k)
                        + innerRes
                        ! ( (size1 - i) `mod` size1,
                            (size2 - j) `mod` size2,
                            (size3 - k) `mod` size3
                          )
                    )
                | i <- [0 .. size1 - 1],
                  j <- [0 .. size2 - 1],
                  k <- [0 .. size3 - 1]
              ]
      TwiceImFT arg ->
        let innerRes = evaluate3DReal valMap $ (mp, arg)
            scaleFactor =
              fromIntegral size1 * fromIntegral size2
                * fromIntegral size3
                / 2
         in listArray
              ((0, 0, 0), (size1 - 1, size2 - 1, size3 - 1))
              [ scaleFactor
                  * ( innerRes ! (i, j, k)
                        - innerRes
                        ! ( (size1 - i) `mod` size1,
                            (size2 - j) `mod` size2,
                            (size3 - k) `mod` size3
                          )
                    )
                | i <- [0 .. size1 - 1],
                  j <- [0 .. size2 - 1],
                  k <- [0 .. size3 - 1]
              ]
      ReFT arg ->
        case retrieveElementType arg mp of
          R ->
            let inner =
                  fmap (:+ 0) . evaluate3DReal valMap $ (mp, arg)
                ftResult =
                  fourierTransform3D (size1, size2, size3) inner
             in fmap realPart ftResult
          C ->
            let inner = evaluate3DComplex valMap $ (mp, arg)
                ftResult =
                  fourierTransform3D (size1, size2, size3) inner
             in fmap realPart ftResult
      ImFT arg ->
        case retrieveElementType arg mp of
          R ->
            let inner =
                  fmap (:+ 0) . evaluate3DReal valMap $ (mp, arg)
                ftResult =
                  fourierTransform3D (size1, size2, size3) inner
             in fmap imagPart ftResult
          C ->
            let inner = evaluate3DComplex valMap $ (mp, arg)
                ftResult =
                  fourierTransform3D (size1, size2, size3) inner
             in fmap imagPart ftResult
      _ -> error "expression structure Three R is wrong"
  | otherwise = error "Three r but shape is not [size1, size2, size3] ??"

instance
  (KnownNat m, KnownNat n, KnownNat p) =>
  Evaluable '(m, n, p) R (Array (Int, Int, Int) Double)
  where
  eval :: ValMaps -> Expression '(m, n, p) R -> Array (Int, Int, Int) Double
  eval valMap (Expression n mp) = evaluate3DReal valMap (mp, n)

evaluate3DComplex ::
  ValMaps -> (ExpressionMap, NodeID) -> Array (Int, Int, Int) (Complex Double)
evaluate3DComplex valMap (mp, n)
  | [size1, size2, size3] <- retrieveShape n mp =
    case retrieveNode n mp of
      Sum C args ->
        foldrElementwise (+) . map (evaluate3DComplex valMap . (mp,)) $
          args
      Mul C args ->
        foldrElementwise (*) . map (evaluate3DComplex valMap . (mp,)) $
          args
      Power x arg -> fmap (^ x) (evaluate3DComplex valMap $ (mp, arg))
      Neg C arg -> fmap negate . evaluate3DComplex valMap $ (mp, arg)
      Scale C arg1 arg2 ->
        case retrieveElementType arg1 mp of
          R ->
            let scalar = fromR . eval valMap $ expZeroR mp arg1
             in fmap (scalar *) . evaluate3DComplex valMap $
                  (mp, arg2)
          C ->
            let scalar = eval valMap $ expZeroC mp arg1
             in fmap (scalar *) . evaluate3DComplex valMap $
                  (mp, arg2)
      RealImag arg1 arg2 ->
        zipWithA
          (:+)
          (evaluate3DReal valMap $ (mp, arg1))
          (evaluate3DReal valMap $ (mp, arg2))
      Piecewise marks conditionArg branchArgs ->
        let cdt = evaluate3DReal valMap $ (mp, conditionArg)
            branches =
              map (evaluate3DComplex valMap . (mp,)) branchArgs
         in listArray
              ((0, 0, 0), (size1 - 1, size2 - 1, size3 - 1))
              [ chosen ! (i, j, k)
                | i <- [0 .. size1 - 1],
                  j <- [0 .. size2 - 1],
                  k <- [0 .. size3 - 1],
                  let chosen =
                        chooseBranch marks (cdt ! (i, j, k)) branches
              ]
      Rotate [amount1, amount2, amount3] arg ->
        rotate3D
          (size1, size2, size3)
          (amount1, amount2, amount3)
          (evaluate3DComplex valMap $ (mp, arg))
      _ -> error "expression structure Three C is wrong"
  | otherwise = error "Three C but shape is not [size1, size2, size3] ??"

instance
  (KnownNat m, KnownNat n, KnownNat p) =>
  Evaluable '(m, n, p) C (Array (Int, Int, Int) (Complex Double))
  where
  eval ::
    ValMaps ->
    Expression '(m, n, p) C ->
    Array (Int, Int, Int) (Complex Double)
  eval valMap (Expression n mp) = evaluate3DComplex valMap (mp, n)

-- NOTE: `mod` in Haskell works with negative number, e.g, (-5) `mod` 3 = 1

-- | One dimension rotation
rotate1D ::
  -- | Size of the input array
  Int ->
  -- | amount of rotation
  Int ->
  -- | Input array
  Array Int a ->
  -- | Rotated array
  Array Int a
rotate1D size amount arr =
  listArray
    (0, size - 1)
    [arr ! ((i - amount) `mod` size) | i <- [0 .. size - 1]]

-- | Two dimension rotation
rotate2D ::
  -- | Size of the 2d input array
  (Int, Int) ->
  -- | amount of rotation for 2d array
  (Int, Int) ->
  -- | Input 2d array
  Array (Int, Int) a ->
  -- | Rotated 2d array
  Array (Int, Int) a
rotate2D (size1, size2) (amount1, amount2) arr =
  listArray
    ((0, 0), (size1 - 1, size2 - 1))
    [ arr ! ((i - amount1) `mod` size1, (j - amount2) `mod` size2)
      | i <- [0 .. size1 - 1],
        j <- [0 .. size2 - 1]
    ]

-- | Three dimension rotation
rotate3D ::
  -- | Size of 3d input array
  (Int, Int, Int) ->
  -- | Amount of Rotation for 3d array
  (Int, Int, Int) ->
  -- | Input 3d array
  Array (Int, Int, Int) a ->
  -- | Rotated 3d array
  Array (Int, Int, Int) a
rotate3D (size1, size2, size3) (amount1, amount2, amount3) arr =
  listArray
    ((0, 0, 0), (size1 - 1, size2 - 1, size3 - 1))
    [ arr
        ! ( (i - amount1) `mod` size1,
            (j - amount2) `mod` size2,
            (k - amount3) `mod` size3
          )
      | i <- [0 .. size1 - 1],
        j <- [0 .. size2 - 1],
        k <- [0 .. size3 - 1]
    ]

-- |
fourierTransform1D ::
  Int -> Array Int (Complex Double) -> Array Int (Complex Double)
fourierTransform1D size arr =
  listArray (0, size - 1) [computeX i | i <- [0 .. size - 1]]
  where
    computeX i = sum $ zipWithA (*) arr (fourierBasis i)
    fourierBasis i =
      let frequency n = 2 * pi * fromIntegral (i * n) / fromIntegral size
       in listArray
            (0, size - 1)
            [ cos (frequency n) :+ (- sin (frequency n))
              | n <- [0 .. size - 1]
            ]

-- |
fourierTransform2D ::
  (Int, Int) ->
  Array (Int, Int) (Complex Double) ->
  Array (Int, Int) (Complex Double)
fourierTransform2D (size1, size2) arr =
  listArray
    ((0, 0), (size1 - 1, size2 - 1))
    [computeX i j | i <- [0 .. size1 - 1], j <- [0 .. size2 - 1]]
  where
    computeX i j = sum $ zipWithA (*) arr (fourierBasis i j)
    fourierBasis i j =
      let frequency m n =
            2 * pi * fromIntegral (i * m) / fromIntegral size1
              + 2 * pi * fromIntegral (j * n) / fromIntegral size2
       in listArray
            ((0, 0), (size1 - 1, size2 - 1))
            [ cos (frequency m n) :+ (- sin (frequency m n))
              | m <- [0 .. size1 - 1],
                n <- [0 .. size2 - 1]
            ]

-- |
fourierTransform3D ::
  (Int, Int, Int) ->
  Array (Int, Int, Int) (Complex Double) ->
  Array (Int, Int, Int) (Complex Double)
fourierTransform3D (size1, size2, size3) arr =
  listArray
    ((0, 0, 0), (size1 - 1, size2 - 1, size3 - 1))
    [ computeX i j k
      | i <- [0 .. size1 - 1],
        j <- [0 .. size2 - 1],
        k <- [0 .. size3 - 1]
    ]
  where
    computeX i j k = sum $ zipWithA (*) arr (fourierBasis i j k)
    fourierBasis i j k =
      let frequency m n p =
            2 * pi * fromIntegral (i * m) / fromIntegral size1
              + 2 * pi * fromIntegral (j * n) / fromIntegral size2
              + 2 * pi * fromIntegral (k * p) / fromIntegral size3
       in listArray
            ((0, 0, 0), (size1 - 1, size2 - 1, size3 - 1))
            [ cos (frequency m n p) :+ (- sin (frequency m n p))
              | m <- [0 .. size1 - 1],
                n <- [0 .. size2 - 1],
                p <- [0 .. size3 - 1]
            ]

max_stars_repo_name: Turboscient/HashedExpression

max_issues_repo_name: Turboscient/HashedExpression

max_forks_repo_name: Turboscient/HashedExpression


text:
{-# LANGUAGE TypeApplications #-}
module Main where

import Config
import File

import Statistics.Distribution
import Statistics.Distribution.Empirical
import Statistics.Distribution.Normal

import Data.List
import Data.Yaml

import Text.Printf

import Debug.Trace

main :: IO ()
main = do
  -- YAML setup
  yml <- decodeFileEither "splits.yaml"

  f <- case yml of
         Right (Config _ levels) -> load levels
         Left  _ -> load'

  let l = 1
      (Just distrs) = lookup l $ levelEmpiricalDistribution f
      (Just level) = fmap (sort . map realToFrac) . lookup l $ onlyValidSplits $ levelData f
      elems = [0, 0.01..0.99]
      dens = 1 / (fromIntegral $ length level)
      domain = reverse [1.0, 1.0 - dens..dens]

      vec = zipWith Vector2 level domain

      vecs = zip4 vec (drop 1 vec) (drop 2 vec) (drop 3 vec)

      vals = concat $ [ [ let (Vector2 x y) = centripetalInterpolate v m in (x, y) | m <- elems ] | v <- vecs ]

      (xs, ys) = unzip vals

      test = zipWith (\(x1, y1) (x2, y2) -> (abs $ 1 / (y2 - y1)) * (x2 - x1)) vals $ tail vals
      (dx, dy) = unzip $ zipWith (\(x1, y1) (x2, y2) -> (,) (x2 - x1) (y2 - y1)) vals $ tail vals

  mapM_ putStrLn $ zipWith5 (printf "%f %f %f %f %f") xs ys dx dy test

linearInterpolate :: (Double, Double) -> Double -> Double
linearInterpolate (lRange, uRange) m =
  let notM = 1 - m
  in lRange * notM + uRange * m

hermiteInterpolate :: (Double, Double, Double, Double) -> Double -> Double
hermiteInterpolate (y0, y1, y2, y3) mu =
  let mu2 = mu * mu
      a0 = (-0.5) * y0 + 1.5 * y1 - 1.5 * y2 + 0.5 * y3
      a1 = y0 - 2.5 * y1 + 2 * y2 - 0.5 * y3
      a2 = (-0.5) * y0 + 0.5 * y2
      a3 = y1
  in a0 * mu * mu2 + a1 * mu2 + a2 * mu + a3

centripetalInterpolate (p0, p1, p2, p3) mu =
  let t0 = 0
      t1 = getT t0 p0 p1
      t2 = getT t1 p1 p2
      t3 = getT t2 p2 p3

      t = linearInterpolate (t1, t2) mu

      a1 = vscale ((t1 - t) / (t1 - t0)) p0 + vscale ((t - t0) / (t1 - t0)) p1
      a2 = vscale ((t2 - t) / (t2 - t1)) p1 + vscale ((t - t1) / (t2 - t1)) p2
      a3 = vscale ((t3 - t) / (t3 - t2)) p2 + vscale ((t - t2) / (t3 - t2)) p3

      b1 = vscale ((t2 - t) / (t2 - t0)) a1 + vscale ((t - t0) / (t2 - t0)) a2
      b2 = vscale ((t3 - t) / (t3 - t1)) a2 + vscale ((t - t1) / (t3 - t1)) a3
  in vscale ((t2 - t) / (t2 - t1)) b1 + vscale ((t - t1) / (t2 - t1)) b2
  where getT tn (Vector2 x1 y1) (Vector2 x2 y2) = tn + ((sqrt $ (x2 - x1) ^ 2 + (y2 - y1) ^ 2) ** alpha)
        alpha = 0.5

type Scalar = Double

class Vector v where
  vmap  :: (Scalar -> Scalar) -> v -> v
  vzip  :: (Scalar -> Scalar -> Scalar) -> v -> v -> v
  vfold :: (x -> Scalar -> x) -> x -> v -> x

vdot :: Vector v => v -> v -> Scalar
vdot v0 v1 = vfold (+) 0 $ vzip (*) v0 v1

vmag_sqr :: Vector v => v -> Scalar
vmag_sqr v = v `vdot` v

vmag :: Vector v => v -> Scalar
vmag = sqrt . vmag_sqr

vscale :: Vector v => Scalar -> v -> v
vscale s = vmap (s*)

vunit :: Vector v => v -> v
vunit v =
  if vmag v == 0
    then v
    else vscale (1 / vmag v) v

data Vector2 = Vector2 {v2x, v2y :: Scalar} deriving (Eq)

instance Show Vector2 where
  show (Vector2 x y) = "<" ++ (show x) ++ ", " ++ (show y) ++ ">"

instance Vector Vector2 where
  vmap  f   (Vector2 x y) = Vector2 (f x) (f y)
  vfold f i (Vector2 x y) = (i `f` x) `f` y
  vzip  f   (Vector2 x0 y0) (Vector2 x1 y1) = Vector2 (f x0 x1) (f y0 y1)

instance Num Vector2 where
  (+) = vzip (+)
  (-) = vzip (-)
  (*) = vzip (*)
  negate = vmap negate
  fromInteger s = Vector2 (fromInteger s) (fromInteger s)

instance Fractional Vector2 where
  (/) = vzip (/)
  fromRational s = let r = realToFrac s in Vector2 r r

ext: hs

max_stars_repo_name: argiopetech/timer

max_issues_repo_name: argiopetech/timer

max_forks_repo_name: argiopetech/timer


text:
{-# LANGUAGE FlexibleContexts, Rank2Types, ScopedTypeVariables #-}
-- | Student's T-test is for assessing whether two samples have
--   different mean. This module contain several variations of
--   T-test. It's a parametric tests and assumes that samples are
--   normally distributed.
module Statistics.Test.StudentT
    (
      studentTTest
    , welchTTest
    , pairedTTest
    , module Statistics.Test.Types
    ) where

import Statistics.Distribution hiding (mean)
import Statistics.Distribution.StudentT
import Statistics.Sample (mean, varianceUnbiased)
import Statistics.Test.Types
import Statistics.Types    (mkPValue,PValue)
import Statistics.Function (square)
import qualified Data.Vector.Generic  as G
import qualified Data.Vector.Unboxed  as U
import qualified Data.Vector.Storable as S
import qualified Data.Vector          as V

-- | Two-sample Student's t-test. It assumes that both samples are
--   normally distributed and have same variance. Returns @Nothing@ if
--   sample sizes are not sufficient.
studentTTest :: (G.Vector v Double)
             => PositionTest  -- ^ one- or two-tailed test
             -> v Double      -- ^ Sample A
             -> v Double      -- ^ Sample B
             -> Maybe (Test StudentT)
studentTTest test sample1 sample2
  | G.length sample1 < 2 || G.length sample2 < 2 = Nothing
  | otherwise                                    = Just Test
      { testSignificance = significance test t ndf
      , testStatistics   = t
      , testDistribution = studentT ndf
      }
  where
    (t, ndf) = tStatistics True sample1 sample2
{-# INLINABLE  studentTTest #-}
{-# SPECIALIZE studentTTest :: PositionTest -> U.Vector Double -> U.Vector Double -> Maybe (Test StudentT) #-}
{-# SPECIALIZE studentTTest :: PositionTest -> S.Vector Double -> S.Vector Double -> Maybe (Test StudentT) #-}
{-# SPECIALIZE studentTTest :: PositionTest -> V.Vector Double -> V.Vector Double -> Maybe (Test StudentT) #-}

-- | Two-sample Welch's t-test. It assumes that both samples are
--   normally distributed but doesn't assume that they have same
--   variance. Returns @Nothing@ if sample sizes are not sufficient.
welchTTest :: (G.Vector v Double)
           => PositionTest  -- ^ one- or two-tailed test
           -> v Double      -- ^ Sample A
           -> v Double      -- ^ Sample B
           -> Maybe (Test StudentT)
welchTTest test sample1 sample2
  | G.length sample1 < 2 || G.length sample2 < 2 = Nothing
  | otherwise                                    = Just Test
      { testSignificance = significance test t ndf
      , testStatistics   = t
      , testDistribution = studentT ndf
      }
  where
    (t, ndf) = tStatistics False sample1 sample2
{-# INLINABLE  welchTTest #-}
{-# SPECIALIZE welchTTest :: PositionTest -> U.Vector Double -> U.Vector Double -> Maybe (Test StudentT) #-}
{-# SPECIALIZE welchTTest :: PositionTest -> S.Vector Double -> S.Vector Double -> Maybe (Test StudentT) #-}
{-# SPECIALIZE welchTTest :: PositionTest -> V.Vector Double -> V.Vector Double -> Maybe (Test StudentT) #-}

-- | Paired two-sample t-test. Two samples are paired in a
-- within-subject design. Returns @Nothing@ if sample size is not
-- sufficient.
pairedTTest :: forall v. (G.Vector v (Double, Double), G.Vector v Double)
            => PositionTest          -- ^ one- or two-tailed test
            -> v (Double, Double)    -- ^ paired samples
            -> Maybe (Test StudentT)
pairedTTest test sample
  | G.length sample < 2 = Nothing
  | otherwise           = Just Test
      { testSignificance = significance test t ndf
      , testStatistics   = t
      , testDistribution = studentT ndf
      }
  where
    (t, ndf) = tStatisticsPaired sample
{-# INLINABLE  pairedTTest #-}
{-# SPECIALIZE pairedTTest :: PositionTest -> U.Vector (Double,Double) -> Maybe (Test StudentT) #-}
{-# SPECIALIZE pairedTTest :: PositionTest -> V.Vector (Double,Double) -> Maybe (Test StudentT) #-}

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

significance :: PositionTest    -- ^ one- or two-tailed
             -> Double          -- ^ t statistics
             -> Double          -- ^ degree of freedom
             -> PValue Double   -- ^ p-value
significance test t df =
  case test of
    -- Here we exploit symmetry of T-distribution and calculate small tail
    SamplesDiffer -> mkPValue $ 2 * tailArea (negate (abs t))
    AGreater      -> mkPValue $ tailArea (negate t)
    BGreater      -> mkPValue $ tailArea  t
  where
    tailArea = cumulative (studentT df)

-- Calculate T statistics for two samples
tStatistics :: (G.Vector v Double)
            => Bool               -- variance equality
            -> v Double
            -> v Double
            -> (Double, Double)
{-# INLINE tStatistics #-}
tStatistics varequal sample1 sample2 = (t, ndf)
  where
    -- t-statistics
    t = (m1 - m2) / sqrt (
      if varequal
        then ((n1 - 1) * s1 + (n2 - 1) * s2) / (n1 + n2 - 2) * (1 / n1 + 1 / n2)
        else s1 / n1 + s2 / n2)

    -- degree of freedom
    ndf | varequal  = n1 + n2 - 2
        | otherwise = square (s1 / n1 + s2 / n2)
                    / (square s1 / (square n1 * (n1 - 1)) + square s2 / (square n2 * (n2 - 1)))
    -- statistics of two samples
    n1 = fromIntegral $ G.length sample1
    n2 = fromIntegral $ G.length sample2
    m1 = mean sample1
    m2 = mean sample2
    s1 = varianceUnbiased sample1
    s2 = varianceUnbiased sample2

-- Calculate T-statistics for paired sample
tStatisticsPaired :: (G.Vector v (Double, Double), G.Vector v Double)
                  => v (Double, Double)
                  -> (Double, Double)
{-# INLINE tStatisticsPaired #-}
tStatisticsPaired sample = (t, ndf)
  where
    -- t-statistics
    t = let d    = G.map (uncurry (-)) sample
            sumd = G.sum d
        in sumd / sqrt ((n * G.sum (G.map square d) - square sumd) / ndf)
    -- degree of freedom
    ndf = n - 1
    n   = fromIntegral $ G.length sample


text:
{-# LANGUAGE ForeignFunctionInterface, GeneralizedNewtypeDeriving #-}

module Numerical.HBLAS.BLAS.FFI.Level2  where

import Foreign.Ptr
import Foreign()
import Foreign.C.Types
import Data.Complex
import Numerical.HBLAS.BLAS.FFI

type GbmvFunFFI sc el =
       CBLAS_ORDERT -> CBLAS_TRANSPOSET -> CInt -> CInt
    -> CInt -> CInt -> sc -> Ptr el -> CInt -> Ptr el -> CInt -> sc -> Ptr el -> CInt -> IO ()

foreign import ccall unsafe "cblas_sgbmv"
    cblas_sgbmv_unsafe :: GbmvFunFFI Float Float
foreign import ccall unsafe "cblas_dgbmv"
    cblas_dgbmv_unsafe :: GbmvFunFFI Double Double
foreign import ccall unsafe "cblas_cgbmv"
    cblas_cgbmv_unsafe :: GbmvFunFFI (Ptr (Complex Float)) (Complex Float)
foreign import ccall unsafe "cblas_zgbmv"
    cblas_zgbmv_unsafe :: GbmvFunFFI (Ptr (Complex Double)) (Complex Double)

foreign import ccall "cblas_sgbmv"
    cblas_sgbmv_safe :: GbmvFunFFI Float Float
foreign import ccall "cblas_dgbmv"
    cblas_dgbmv_safe :: GbmvFunFFI Double Double
foreign import ccall "cblas_cgbmv"
    cblas_cgbmv_safe :: GbmvFunFFI (Ptr (Complex Float)) (Complex Float)
foreign import ccall "cblas_zgbmv"
    cblas_zgbmv_safe :: GbmvFunFFI (Ptr (Complex Double)) (Complex Double)
--void cblas_sgbmv(  enum CBLAS_ORDER order,   enum CBLAS_TRANSPOSE TransA,   CInt M,   CInt N,
--                   CInt KL,   CInt KU,   Float alpha,   Float *A,   CInt lda,   Float *X,   CInt incX,   Float beta, Float *Y,   CInt incY);
--void cblas_dgbmv(  enum CBLAS_ORDER order,   enum CBLAS_TRANSPOSE TransA,   CInt M,   CInt N,
--                   CInt KL,   CInt KU,   Double alpha,   Double *A,   CInt lda,   Double *X,   CInt incX,   Double beta, Double *Y,   CInt incY);
--void cblas_cgbmv(  enum CBLAS_ORDER order,   enum CBLAS_TRANSPOSE TransA,   CInt M,   CInt N,
--                   CInt KL,   CInt KU,   Float *alpha,   Float *A,   CInt lda,   Float *X,   CInt incX,   Float *beta, Float *Y,   CInt incY);
--void cblas_zgbmv(  enum CBLAS_ORDER order,   enum CBLAS_TRANSPOSE TransA,   CInt M,   CInt N,
--                   CInt KL,   CInt KU,   Double *alpha,   Double *A,   CInt lda,   Double *X,   CInt incX,   Double *beta, Double *Y,   CInt incY);

{-
matrix vector product for general matrices
 perform one of the matrix-vector operations   y :=
      alpha*A*x + beta*y, or y := alpha*A'*x + beta*y,
-}

type GemvFunFFI sc el =
       CBLAS_ORDERT -> CBLAS_TRANSPOSET -> CInt -> CInt
    -> sc -> Ptr el -> CInt -> Ptr el -> CInt -> sc -> Ptr el -> CInt -> IO ()

foreign import ccall unsafe "cblas_sgemv"
    cblas_sgemv_unsafe :: GemvFunFFI Float Float
foreign import ccall safe   "cblas_sgemv"
    cblas_sgemv_safe   :: GemvFunFFI Float Float

foreign import ccall unsafe "cblas_dgemv"
    cblas_dgemv_unsafe :: GemvFunFFI Double Double
foreign import ccall safe   "cblas_dgemv"
    cblas_dgemv_safe   :: GemvFunFFI Double Double

foreign import ccall unsafe "cblas_cgemv"
    cblas_cgemv_unsafe :: GemvFunFFI (Ptr (Complex Float)) (Complex Float)
foreign import ccall safe   "cblas_cgemv"
    cblas_cgemv_safe   :: GemvFunFFI (Ptr (Complex Float)) (Complex Float)

foreign import ccall unsafe "cblas_zgemv"
    cblas_zgemv_unsafe :: GemvFunFFI (Ptr (Complex Double)) (Complex Double)
foreign import ccall safe   "cblas_zgemv"
    cblas_zgemv_safe   :: GemvFunFFI (Ptr (Complex Double)) (Complex Double)

--void cblas_sgemv(  enum CBLAS_ORDER order,    enum CBLAS_TRANSPOSE trans,    CInt m,   CInt n,
--           Float alpha,   Float  *a,   CInt lda,    Float  *x,   CInt incx,    Float beta,  Float  *y,   CInt incy);
--void cblas_dgemv(  enum CBLAS_ORDER order,    enum CBLAS_TRANSPOSE trans,    CInt m,   CInt n,
--           Double alpha,   Double  *a,   CInt lda,    Double  *x,   CInt incx,    Double beta,  Double  *y,   CInt incy);
--void cblas_cgemv(  enum CBLAS_ORDER order,    enum CBLAS_TRANSPOSE trans,    CInt m,   CInt n,
--           Float *alpha,   Float  *a,   CInt lda,    Float  *x,   CInt incx,    Float *beta,  Float  *y,   CInt incy);
--void cblas_zgemv(  enum CBLAS_ORDER order,    enum CBLAS_TRANSPOSE trans,    CInt m,   CInt n,
--           Double *alpha,   Double  *a,   CInt lda,    Double  *x,   CInt incx,    Double *beta,  Double  *y,   CInt incy);

-- perform the rank 1 operation   A := alpha*x*y' + A,

type GerxFunFFI scale el = CBLAS_ORDERT -> CInt -> CInt -> scale -> Ptr el -> CInt -> Ptr el -> CInt -> Ptr el -> CInt -> IO ()

foreign import ccall unsafe "cblas_sger" cblas_sger_unsafe ::
        GerxFunFFI Float Float
foreign import ccall safe   "cblas_sger" cblas_sger_safe   ::
        GerxFunFFI Float Float

foreign import ccall unsafe "cblas_dger" cblas_dger_unsafe ::
        GerxFunFFI Double Double
foreign import ccall safe   "cblas_dger" cblas_dger_safe   ::
        GerxFunFFI Double Double

foreign import ccall unsafe "cblas_cgerc" cblas_cgerc_unsafe ::
        GerxFunFFI (Ptr (Complex Float)) (Complex Float)
foreign import ccall unsafe "cblas_zgerc" cblas_zgerc_unsafe ::
        GerxFunFFI (Ptr (Complex Double)) (Complex Double)

foreign import ccall safe   "cblas_cgerc" cblas_cgerc_safe   ::
        GerxFunFFI (Ptr (Complex Float)) (Complex Float)
foreign import ccall safe   "cblas_zgerc" cblas_zgerc_safe   ::
        GerxFunFFI (Ptr (Complex Double)) (Complex Double)

foreign import ccall unsafe "cblas_cgeru" cblas_cgeru_unsafe ::
        GerxFunFFI (Ptr (Complex Float)) (Complex Float)
foreign import ccall unsafe "cblas_zgeru" cblas_zgeru_unsafe ::
        GerxFunFFI (Ptr (Complex Double)) (Complex Double)

foreign import ccall safe   "cblas_cgeru" cblas_cgeru_safe   ::
        GerxFunFFI (Ptr (Complex Float)) (Complex Float)
foreign import ccall safe   "cblas_zgeru" cblas_zgeru_safe   ::
        GerxFunFFI (Ptr (Complex Double)) (Complex Double)

--void cblas_sger (  enum CBLAS_ORDER order,   CInt M,   CInt N,   Float   alpha,   Float  *X,   CInt incX,   Float  *Y,   CInt incY, Float  *A,   CInt lda);
--void cblas_dger (  enum CBLAS_ORDER order,   CInt M,   CInt N,   Double  alpha,   Double *X,   CInt incX,   Double *Y,   CInt incY, Double *A,   CInt lda);
--void cblas_cgeru(  enum CBLAS_ORDER order,   CInt M,   CInt N,   Float  *alpha,   Float  *X,   CInt incX,   Float  *Y,   CInt incY, Float  *A,   CInt lda);
--void cblas_cgerc(  enum CBLAS_ORDER order,   CInt M,   CInt N,   Float  *alpha,   Float  *X,   CInt incX,   Float  *Y,   CInt incY, Float  *A,   CInt lda);
--void cblas_zgeru(  enum CBLAS_ORDER order,   CInt M,   CInt N,   Double *alpha,   Double *X,   CInt incX,   Double *Y,   CInt incY, Double *A,   CInt lda);
--void cblas_zgerc(  enum CBLAS_ORDER order,   CInt M,   CInt N,   Double *alpha,   Double *X,   CInt incX,   Double *Y,   CInt incY, Double *A,   CInt lda);

type HbmvFunFFI sc el =
       CBLAS_ORDERT -> CBLAS_UPLOT -> CInt -> CInt
    -> sc -> Ptr el -> CInt -> Ptr el -> CInt -> sc -> Ptr el -> CInt -> IO ()
foreign import ccall unsafe "cblas_chbmv"
    cblas_chbmv_unsafe :: HbmvFunFFI (Ptr (Complex Float)) (Complex Float)
foreign import ccall safe   "cblas_chbmv"
    cblas_chbmv_safe   :: HbmvFunFFI (Ptr (Complex Float)) (Complex Float)

foreign import ccall unsafe "cblas_zhbmv"
    cblas_zhbmv_unsafe :: HbmvFunFFI (Ptr (Complex Double)) (Complex Double)
foreign import ccall safe   "cblas_zhbmv"
    cblas_zhbmv_safe   :: HbmvFunFFI (Ptr (Complex Double)) (Complex Double)
--void cblas_chbmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   CInt K,
--           Float *alpha,   Float *A,   CInt lda,   Float *X,   CInt incX,   Float *beta, Float *Y,   CInt incY);
--void cblas_zhbmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   CInt K,
--           Double *alpha,   Double *A,   CInt lda,   Double *X,   CInt incX,   Double *beta, Double *Y,   CInt incY);

type HemvFunFFI sc el =
       CBLAS_ORDERT -> CBLAS_UPLOT -> CInt
    -> sc -> Ptr el -> CInt -> Ptr el -> CInt -> sc -> Ptr el -> CInt -> IO ()
foreign import ccall unsafe "cblas_chemv"
    cblas_chemv_unsafe :: HemvFunFFI (Ptr (Complex Float)) (Complex Float)
foreign import ccall safe   "cblas_chemv"
    cblas_chemv_safe   :: HemvFunFFI (Ptr (Complex Float)) (Complex Float)

foreign import ccall unsafe "cblas_zhemv"
    cblas_zhemv_unsafe :: HemvFunFFI (Ptr (Complex Double)) (Complex Double)
foreign import ccall safe   "cblas_zhemv"
    cblas_zhemv_safe   :: HemvFunFFI (Ptr (Complex Double)) (Complex Double)
--void cblas_chemv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Float *alpha,   Float *A,
--                   CInt lda,   Float *X,   CInt incX,   Float *beta, Float *Y,   CInt incY);
--void cblas_zhemv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double *alpha,   Double *A,
--                   CInt lda,   Double *X,   CInt incX,   Double *beta, Double *Y,   CInt incY);

type HerFunFFI sc el =
      CBLAS_ORDERT -> CBLAS_UPLOT -> CInt
      -> sc -> Ptr el -> CInt -> Ptr el -> CInt -> IO ()
foreign import ccall unsafe "cblas_cher"
    cblas_cher_unsafe :: HerFunFFI Float (Complex Float)
foreign import ccall safe "cblas_cher"
    cblas_cher_safe :: HerFunFFI Float (Complex Float)

foreign import ccall unsafe "cblas_zher"
    cblas_zher_unsafe :: HerFunFFI Double (Complex Double)
foreign import ccall safe "cblas_zher"
    cblas_zher_safe :: HerFunFFI Double (Complex Double)
--void cblas_cher(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Float alpha,   Float *X,   CInt incX, Float *A,   CInt lda);
--void cblas_zher(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double alpha,   Double *X,   CInt incX, Double *A,   CInt lda);

type Her2FunFFI sc el =
      CBLAS_ORDERT -> CBLAS_UPLOT -> CInt
      -> sc -> Ptr el -> CInt -> Ptr el -> CInt -> Ptr el -> CInt -> IO ()
foreign import ccall unsafe "cblas_cher2"
    cblas_cher2_unsafe :: Her2FunFFI (Ptr (Complex Float)) (Complex Float)
foreign import ccall safe "cblas_cher2"
    cblas_cher2_safe :: Her2FunFFI (Ptr (Complex Float)) (Complex Float)

foreign import ccall unsafe "cblas_zher2"
    cblas_zher2_unsafe :: Her2FunFFI (Ptr (Complex Double)) (Complex Double)
foreign import ccall safe "cblas_zher2"
    cblas_zher2_safe :: Her2FunFFI (Ptr (Complex Double)) (Complex Double)
--void cblas_cher2(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Float *alpha,   Float *X,   CInt incX,
--                  Float *Y,   CInt incY, Float *A,   CInt lda);
--void cblas_zher2(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double *alpha,   Double *X,   CInt incX,
--                  Double *Y,   CInt incY, Double *A,   CInt lda);

type HpmvFunFFI sc el =
      CBLAS_ORDERT -> CBLAS_UPLOT -> CInt
      -> sc -> Ptr el -> Ptr el -> CInt -> sc -> Ptr el -> CInt -> IO ()
foreign import ccall unsafe "cblas_chpmv"
    cblas_chpmv_unsafe :: HpmvFunFFI (Ptr (Complex Float)) (Complex Float)
foreign import ccall safe "cblas_chpmv"
    cblas_chpmv_safe :: HpmvFunFFI (Ptr (Complex Float)) (Complex Float)

foreign import ccall unsafe "cblas_zhpmv"
    cblas_zhpmv_unsafe :: HpmvFunFFI (Ptr (Complex Double)) (Complex Double)
foreign import ccall safe "cblas_zhpmv"
    cblas_zhpmv_safe :: HpmvFunFFI (Ptr (Complex Double)) (Complex Double)
--void cblas_chpmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,
--           Float *alpha,   Float *Ap,   Float *X,   CInt incX,   Float *beta, Float *Y,   CInt incY);
--void cblas_zhpmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,
--           Double *alpha,   Double *Ap,   Double *X,   CInt incX,   Double *beta, Double *Y,   CInt incY);

type HprFunFFI sc el =
      CBLAS_ORDERT -> CBLAS_UPLOT -> CInt
      -> sc -> Ptr el -> CInt -> Ptr el -> IO ()
foreign import ccall unsafe "cblas_chpr"
    cblas_chpr_unsafe :: HprFunFFI Float (Complex Float)
foreign import ccall safe "cblas_chpr"
    cblas_chpr_safe :: HprFunFFI Float (Complex Float)

foreign import ccall unsafe "cblas_zhpr"
    cblas_zhpr_unsafe :: HprFunFFI Double (Complex Double)
foreign import ccall safe "cblas_zhpr"
    cblas_zhpr_safe :: HprFunFFI Double (Complex Double)
--void cblas_chpr(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Float alpha,   Float *X,   CInt incX, Float *A);
--void cblas_zhpr(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double alpha,   Double *X,  CInt incX, Double *A);

type Hpr2FunFFI sc el =
      CBLAS_ORDERT -> CBLAS_UPLOT -> CInt
      -> sc -> Ptr el -> CInt -> Ptr el -> CInt -> Ptr el -> IO ()
foreign import ccall unsafe "cblas_chpr2"
    cblas_chpr2_unsafe :: Hpr2FunFFI (Ptr (Complex Float)) (Complex Float)
foreign import ccall safe "cblas_chpr2"
    cblas_chpr2_safe :: Hpr2FunFFI (Ptr (Complex Float)) (Complex Float)

foreign import ccall unsafe "cblas_zhpr2"
    cblas_zhpr2_unsafe :: Hpr2FunFFI (Ptr (Complex Double)) (Complex Double)
foreign import ccall safe "cblas_zhpr2"
    cblas_zhpr2_safe :: Hpr2FunFFI (Ptr (Complex Double)) (Complex Double)
--void cblas_chpr2(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Float *alpha,   Float *X,   CInt incX,   Float *Y,   CInt incY, Float *Ap);
--void cblas_zhpr2(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double *alpha,   Double *X,   CInt incX,   Double *Y,   CInt incY, Double *Ap);

type SbmvFunFFI sc el =
      CBLAS_ORDERT -> CBLAS_UPLOT -> CInt -> CInt
      -> sc -> Ptr el -> CInt -> Ptr el -> CInt -> sc -> Ptr el -> CInt -> IO ()
foreign import ccall unsafe "cblas_ssbmv"
    cblas_ssbmv_unsafe :: SbmvFunFFI Float Float
foreign import ccall safe "cblas_ssbmv"
    cblas_ssbmv_safe :: SbmvFunFFI Float Float

foreign import ccall unsafe "cblas_dsbmv"
    cblas_dsbmv_unsafe :: SbmvFunFFI Double Double
foreign import ccall safe "cblas_dsbmv"
    cblas_dsbmv_safe :: SbmvFunFFI Double Double
--void cblas_ssbmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   CInt K,   Float alpha,   Float *A,
--                   CInt lda,   Float *X,   CInt incX,   Float beta, Float *Y,   CInt incY);
--void cblas_dsbmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   CInt K,   Double alpha,   Double *A,
--                   CInt lda,   Double *X,   CInt incX,   Double beta, Double *Y,   CInt incY);

---------------
--- | packed symmetric matrix * vector product  y:= alpha * Av  + beta * y
---------------
type SpmvFunFFI sc el =
      CBLAS_ORDERT -> CBLAS_UPLOT -> CInt
      -> sc -> Ptr el -> Ptr el -> CInt -> sc -> Ptr el -> CInt -> IO ()
foreign import ccall unsafe "cblas_sspmv"
    cblas_sspmv_unsafe :: SpmvFunFFI Float Float
foreign import ccall safe "cblas_sspmv"
    cblas_sspmv_safe :: SpmvFunFFI Float Float

foreign import ccall unsafe "cblas_dspmv"
    cblas_dspmv_unsafe :: SpmvFunFFI Double Double
foreign import ccall safe "cblas_dspmv"
    cblas_dspmv_safe :: SpmvFunFFI Double Double
--void cblas_sspmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Float alpha,   Float *Ap,
--                   Float *X,   CInt incX,   Float beta, Float *Y,   CInt incY);
--void cblas_dspmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double alpha,   Double *Ap,
--                   Double *X,   CInt incX,   Double beta, Double *Y,   CInt incY);

type SprFunFFI sc el =
      CBLAS_ORDERT -> CBLAS_UPLOT -> CInt
      -> sc -> Ptr el -> CInt -> Ptr el -> IO ()
foreign import ccall unsafe "cblas_sspr"
    cblas_sspr_unsafe :: SprFunFFI Float Float
foreign import ccall safe "cblas_sspr"
    cblas_sspr_safe :: SprFunFFI Float Float

foreign import ccall unsafe "cblas_dspr"
    cblas_dspr_unsafe :: SprFunFFI Double Double
foreign import ccall safe "cblas_dspr"
    cblas_dspr_safe :: SprFunFFI Double Double
--void cblas_sspr(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Float alpha,   Float *X,   CInt incX, Float *Ap);
--void cblas_dspr(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double alpha,   Double *X,   CInt incX, Double *Ap);

type Spr2FunFFI sc el =
      CBLAS_ORDERT -> CBLAS_UPLOT -> CInt
      -> sc -> Ptr el -> CInt -> Ptr el -> CInt -> Ptr el -> IO ()
foreign import ccall unsafe "cblas_sspr2"
    cblas_sspr2_unsafe :: Spr2FunFFI Float Float
foreign import ccall safe "cblas_sspr2"
    cblas_sspr2_safe :: Spr2FunFFI Float Float

foreign import ccall unsafe "cblas_dspr2"
    cblas_dspr2_unsafe :: Spr2FunFFI Double Double
foreign import ccall safe "cblas_dspr2"
    cblas_dspr2_safe :: Spr2FunFFI Double Double
--void cblas_sspr2(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Float alpha,   Float *X,   CInt incX,   Float *Y,   CInt incY, Float *A);
--void cblas_dspr2(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double alpha,   Double *X,   CInt incX,   Double *Y,   CInt incY, Double *A);

----------------------------------
---- |  (unpacked) symmetric matrix vector product    x:=Av, writes result x into v
---------------------------------

type SymvFunFFI el = CBLAS_ORDERT -> CBLAS_UPLOT -> CInt -> el -> Ptr el ->  CInt ->
                        Ptr el -> CInt -> el -> Ptr el -> CInt -> IO ()
foreign import ccall unsafe "cblas_ssymv"
    cblas_ssymv_unsafe :: SymvFunFFI Float
foreign import ccall safe "cblas_ssymv"
    cblas_ssymv_safe :: SymvFunFFI Float

foreign import ccall unsafe "cblas_dsymv"
    cblas_dsymv_unsafe :: SymvFunFFI Double
foreign import ccall safe "cblas_dsymv"
    cblas_dsymv_safe :: SymvFunFFI Double
--void cblas_ssymv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Float alpha,   Float *A,
--                   CInt lda,   Float *X,   CInt incX,   Float beta, Float *Y,   CInt incY);
--void cblas_dsymv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double alpha,   Double *A,
--                   CInt lda,   Double *X,   CInt incX,   Double beta, Double *Y,   CInt incY);

type SyrFunFFI el = CBLAS_ORDERT -> CBLAS_UPLOT -> CInt -> el -> Ptr el -> CInt -> Ptr el -> CInt -> IO ()
foreign import ccall unsafe "cblas_ssyr"
      cblas_ssyr_unsafe :: SyrFunFFI Float
foreign import ccall safe "cblas_ssyr"
      cblas_ssyr_safe :: SyrFunFFI Float

foreign import ccall unsafe "cblas_dsyr"
      cblas_dsyr_unsafe :: SyrFunFFI Double
foreign import ccall safe "cblas_dsyr"
      cblas_dsyr_safe :: SyrFunFFI Double
--void cblas_ssyr(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Float alpha,   Float *X,   CInt incX, Float *A,   CInt lda);
--void cblas_dsyr(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double alpha,   Double *X,   CInt incX, Double *A,   CInt lda);

type Syr2FunFFI el = CBLAS_ORDERT -> CBLAS_UPLOT -> CInt -> el -> Ptr el -> CInt -> Ptr el -> CInt -> Ptr el -> CInt -> IO ()
foreign import ccall unsafe "cblas_ssyr2"
      cblas_ssyr2_unsafe :: Syr2FunFFI Float
foreign import ccall safe "cblas_ssyr2"
      cblas_ssyr2_safe :: Syr2FunFFI Float

foreign import ccall unsafe "cblas_dsyr2"
      cblas_dsyr2_unsafe :: Syr2FunFFI Double
foreign import ccall safe "cblas_dsyr2"
      cblas_dsyr2_safe :: Syr2FunFFI Double
--void cblas_ssyr2(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,  CInt N,   Float alpha,   Float *X,
--                  CInt incX,   Float *Y,   CInt incY, Float *A,   CInt lda);
--void cblas_dsyr2(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double alpha,   Double *X,
--                  CInt incX,   Double *Y,   CInt incY, Double *A,   CInt lda);

type TbmvFunFFI el = CBLAS_ORDERT -> CBLAS_UPLOT -> CBLAS_TRANSPOSET -> CBLAS_DIAGT ->
                      CInt -> CInt -> Ptr el -> CInt -> Ptr el -> CInt -> IO ()
foreign import ccall unsafe "cblas_stbmv"
      cblas_stbmv_unsafe :: TbmvFunFFI Float
foreign import ccall safe "cblas_stbmv"
      cblas_stbmv_safe :: TbmvFunFFI Float

foreign import ccall unsafe "cblas_dtbmv"
      cblas_dtbmv_unsafe :: TbmvFunFFI Double
foreign import ccall safe "cblas_dtbmv"
      cblas_dtbmv_safe :: TbmvFunFFI Double

foreign import ccall unsafe "cblas_ctbmv"
      cblas_ctbmv_unsafe :: TbmvFunFFI (Complex Float)
foreign import ccall safe "cblas_ctbmv"
      cblas_ctbmv_safe :: TbmvFunFFI (Complex Float)

foreign import ccall unsafe "cblas_ztbmv"
      cblas_ztbmv_unsafe :: TbmvFunFFI (Complex Double)
foreign import ccall safe "cblas_ztbmv"
      cblas_ztbmv_safe :: TbmvFunFFI (Complex Double)
--void cblas_stbmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   CInt K,   Float *A,   CInt lda, Float *X,   CInt incX);
--void cblas_dtbmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   CInt K,   Double *A,   CInt lda, Double *X,   CInt incX);
--void cblas_ctbmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   CInt K,   Float *A,   CInt lda, Float *X,   CInt incX);
--void cblas_ztbmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   CInt K,   Double *A,   CInt lda, Double *X,   CInt incX);

----------------
--- | solves  Ax=v where A is k+1 banded triangular matrix, and x and
----------------
type TbsvFunFFI el = CBLAS_ORDERT -> CBLAS_UPLOT -> CBLAS_TRANSPOSET -> CBLAS_DIAGT ->
                      CInt -> CInt -> Ptr el -> CInt -> Ptr el -> CInt -> IO ()
foreign import ccall unsafe "cblas_stbsv"
      cblas_stbsv_unsafe :: TbsvFunFFI Float
foreign import ccall safe "cblas_stbsv"
      cblas_stbsv_safe :: TbsvFunFFI Float

foreign import ccall unsafe "cblas_dtbsv"
      cblas_dtbsv_unsafe :: TbsvFunFFI Double
foreign import ccall safe "cblas_dtbsv"
      cblas_dtbsv_safe :: TbsvFunFFI Double

foreign import ccall unsafe "cblas_ctbsv"
      cblas_ctbsv_unsafe :: TbsvFunFFI (Complex Float)
foreign import ccall safe "cblas_ctbsv"
      cblas_ctbsv_safe :: TbsvFunFFI (Complex Float)

foreign import ccall unsafe "cblas_ztbsv"
      cblas_ztbsv_unsafe :: TbsvFunFFI (Complex Double)
foreign import ccall safe "cblas_ztbsv"
      cblas_ztbsv_safe :: TbsvFunFFI (Complex Double)
--void cblas_stbsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   CInt K,   Float *A,   CInt lda, Float *X,   CInt incX);
--void cblas_dtbsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   CInt K,   Double *A,   CInt lda, Double *X,   CInt incX);
--void cblas_ctbsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   CInt K,   Float *A,   CInt lda, Float *X,   CInt incX);
--void cblas_ztbsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   CInt K,   Double *A,   CInt lda, Double *X,   CInt incX);

-------------------------------------------------------------------------
-- | matrix vector product Av, writes result into v, where A is a packed triangular nxn matrix
-------------------------------------------------------------------------
type TpmvFunFFI el = CBLAS_ORDERT -> CBLAS_UPLOT -> CBLAS_TRANSPOSET -> CBLAS_DIAGT ->
                      CInt -> Ptr el -> Ptr el -> CInt -> IO ()
foreign import ccall unsafe "cblas_stpmv"
      cblas_stpmv_unsafe :: TpmvFunFFI Float
foreign import ccall safe "cblas_stpmv"
      cblas_stpmv_safe :: TpmvFunFFI Float

foreign import ccall unsafe "cblas_dtpmv"
      cblas_dtpmv_unsafe :: TpmvFunFFI Double
foreign import ccall safe "cblas_dtpmv"
      cblas_dtpmv_safe :: TpmvFunFFI Double

foreign import ccall unsafe "cblas_ctpmv"
      cblas_ctpmv_unsafe :: TpmvFunFFI (Complex Float)
foreign import ccall safe "cblas_ctpmv"
      cblas_ctpmv_safe :: TpmvFunFFI (Complex Float)

foreign import ccall unsafe "cblas_ztpmv"
      cblas_ztpmv_unsafe :: TpmvFunFFI (Complex Double)
foreign import ccall safe "cblas_ztpmv"
      cblas_ztpmv_safe :: TpmvFunFFI (Complex Double)
--void cblas_stpmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   Float *Ap, Float *X,   CInt incX);
--void cblas_dtpmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   Double *Ap, Double *X,   CInt incX);
--void cblas_ctpmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   Float *Ap, Float *X,   CInt incX);
--void cblas_ztpmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   Double *Ap, Double *X,   CInt incX);

--------------------------------------------------
---  | solve  Ax=v where A is a nxn packed triangular matrix, v vector input, writes the solution into x.
--------------------------------------------------
type TpsvFunFFI el = CBLAS_ORDERT -> CBLAS_UPLOT -> CBLAS_TRANSPOSET -> CBLAS_DIAGT ->
                      CInt -> Ptr el -> Ptr el -> CInt -> IO ()
foreign import ccall unsafe "cblas_stpsv"
      cblas_stpsv_unsafe :: TpsvFunFFI Float
foreign import ccall safe "cblas_stpsv"
      cblas_stpsv_safe :: TpsvFunFFI Float

foreign import ccall unsafe "cblas_dtpsv"
      cblas_dtpsv_unsafe :: TpsvFunFFI Double
foreign import ccall safe "cblas_dtpsv"
      cblas_dtpsv_safe :: TpsvFunFFI Double

foreign import ccall unsafe "cblas_ctpsv"
      cblas_ctpsv_unsafe :: TpsvFunFFI (Complex Float)
foreign import ccall safe "cblas_ctpsv"
      cblas_ctpsv_safe :: TpsvFunFFI (Complex Float)

foreign import ccall unsafe "cblas_ztpsv"
      cblas_ztpsv_unsafe :: TpsvFunFFI (Complex Double)
foreign import ccall safe "cblas_ztpsv"
      cblas_ztpsv_safe :: TpsvFunFFI (Complex Double)
--void cblas_stpsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   Float *Ap, Float *X,   CInt incX);
--void cblas_dtpsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   Double *Ap, Double *X,   CInt incX);
--void cblas_ctpsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   Float *Ap, Float *X,   CInt incX);
--void cblas_ztpsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   Double *Ap, Double *X,   CInt incX);

type TrmvFunFFI el = CBLAS_ORDERT -> CBLAS_UPLOT -> CBLAS_TRANSPOSET -> CBLAS_DIAGT ->
                      CInt -> Ptr el -> CInt -> Ptr el -> CInt -> IO ()
foreign import ccall unsafe "cblas_strmv"
      cblas_strmv_unsafe :: TrmvFunFFI Float
foreign import ccall safe "cblas_strmv"
      cblas_strmv_safe :: TrmvFunFFI Float

foreign import ccall unsafe "cblas_dtrmv"
      cblas_dtrmv_unsafe :: TrmvFunFFI Double
foreign import ccall safe "cblas_dtrmv"
      cblas_dtrmv_safe :: TrmvFunFFI Double

foreign import ccall unsafe "cblas_ctrmv"
      cblas_ctrmv_unsafe :: TrmvFunFFI (Complex Float)
foreign import ccall safe "cblas_ctrmv"
      cblas_ctrmv_safe :: TrmvFunFFI (Complex Float)

foreign import ccall unsafe "cblas_ztrmv"
      cblas_ztrmv_unsafe :: TrmvFunFFI (Complex Double)
foreign import ccall safe "cblas_ztrmv"
      cblas_ztrmv_safe :: TrmvFunFFI (Complex Double)
--void cblas_strmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,   CInt N,   Float *A,   CInt lda, Float *X,   CInt incX);
--void cblas_dtrmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,   CInt N,   Double *A,   CInt lda, Double *X,   CInt incX);
--void cblas_ctrmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,   CInt N,   Float *A,   CInt lda, Float *X,   CInt incX);
--void cblas_ztrmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,   CInt N,   Double *A,   CInt lda, Double *X,   CInt incX);

--STRSV - solve one of the systems of equations   A*x = b, or A'*x = b, where A is a (non)unit upper(/lower) triangular matrix
type TrsvFunFFI el = CBLAS_ORDERT -> CBLAS_UPLOT -> CBLAS_TRANSPOSET -> CBLAS_DIAGT ->
                      CInt -> Ptr el -> CInt -> Ptr el -> CInt -> IO ()

foreign import ccall unsafe "cblas_strsv"
  cblas_strsv_unsafe :: TrsvFunFFI Float
foreign import ccall safe   "cblas_strsv"
  cblas_strsv_safe   :: TrsvFunFFI Float

foreign import ccall unsafe "cblas_dtrsv"
  cblas_dtrsv_unsafe :: TrsvFunFFI Double
foreign import ccall safe   "cblas_dtrsv"
  cblas_dtrsv_safe   :: TrsvFunFFI Double

foreign import ccall unsafe "cblas_ctrsv"
  cblas_ctrsv_unsafe :: TrsvFunFFI (Complex Float)
foreign import ccall safe   "cblas_ctrsv"
  cblas_ctrsv_safe   :: TrsvFunFFI (Complex Float)

foreign import ccall unsafe "cblas_ztrsv"
  cblas_ztrsv_unsafe :: TrsvFunFFI (Complex Double)
foreign import ccall safe   "cblas_ztrsv"
  cblas_ztrsv_safe   :: TrsvFunFFI (Complex Double)

--void cblas_strsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,   CInt N,   Float *A,   CInt lda, Float *X,   CInt incX);
--void cblas_dtrsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,   CInt N,   Double *A,   CInt lda, Double *X,   CInt incX);
--void cblas_ctrsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,   CInt N,   Float *A,   CInt lda, Float *X,   CInt incX);
--void cblas_ztrsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,   CInt N,   Double *A,   CInt lda, Double *X,   CInt incX);

max_stars_repo_name: schnecki/hblas

max_issues_repo_name: schnecki/hblas

max_forks_repo_name: schnecki/hblas


text:
--
-- FullConnLayer : fully connected layer
--

module CNN.FullConnLayer (
  initFilterF
, zeroFilterF
, connect
, deconnect
, reverseFullConnFilter
, updateFullConnFilter
) where

import Control.Monad hiding (msum)
--import Data.List
import Data.Maybe
import Debug.Trace
import Numeric.LinearAlgebra
import System.Random.Mersenne as MT

import CNN.Algebra
import CNN.Image
import CNN.LayerType

{- |
initFilterF

  IN : #kernel
       #channel

  OUT: filter of fully connected layer

  OUT: filter of fully connected layer

>>> f <- initFilterF 3 20
>>> rows f
3
>>> cols f
21

-}

initFilterF :: Int -> Int -> IO FilterF
initFilterF k c = do
  rs <- forM [1..k] $ \i -> do
    w <- forM [1..c] $ \j -> do
      r <- MT.randomIO :: IO Double
      return ((r * 2.0 - 1.0) * a)
    return (0.0:w)
  return $ fromLists rs
  where
    --a = 1.0 / fromIntegral c
    a = 4.0 * sqrt (6.0 / fromIntegral (c + k))
    --a = sqrt (6.0 / fromIntegral (c + k))

zeroFilterF :: Int -> Int -> IO FilterF
zeroFilterF k c = return $ (k><(c+1)) $ repeat 0.0

--

{- |
connect

  IN : filter of fully connected layer
       image

  OUT: updated image

>>> let fs = fromLists [[0.5,1.0,2.0,3.0],[0.1,4.0,5.0,6.0]]
>>> let im = [fromLists [[9.0,8.0,7.0]]]
>>> connect fs im
[(1><2)
 [ 46.5, 118.1 ]]
>>> connect fs []
*** Exception: invalid Image
CallStack (from HasCallStack):
  error, called at CNN/FullConnLayer.hs:82:16 in main:CNN.FullConnLayer

-}

connect :: FilterF -> Image -> Image
connect _ [] = error "invalid Image"
connect fs [im] = [reshape (rows fs) $ fs #> v]
  where
    v = vjoin [konst 1.0 1, flatten im]

-- back prop

{- |
deconnect

  IN : filter of fully connected layer
       image
       difference from previous layer

  OUT: difference and updated layer

>>> let im = [fromLists [[1.0, 2.0, 3.0]]]
>>> let delta = [fromLists [[1.0, 2.0]]]
>>> let fs = fromLists [[1.0, 2.0],[3.0,4.0],[5.0,6.0]]
>>> let (d,l) = deconnect fs im delta
>>> d
[(1><2)
 [ 11.0, 17.0 ]]
>>> l
Just FullConnLayer:(2><4)
 [ 1.0, 1.0, 2.0, 3.0
 , 2.0, 2.0, 4.0, 6.0 ]

-}

deconnect :: FilterF -> Image -> Delta -> (Delta, Maybe Layer)
deconnect fs im delta = ([reshape (rows fs') $ fs' #> dl]
                        , Just (FullConnLayer $ outer dl im'))
  where
    dl  = flatten $ head delta    -- to Vector
    fs' = dropRows 1 fs
    im' = vjoin [konst 1.0 1, flatten $ head im]

-- reverse

reverseFullConnFilter :: FilterF -> Layer
reverseFullConnFilter fs = FullConnLayer $ tr' fs

-- update filter

updateFullConnFilter :: FilterF -> Double -> [Maybe Layer] -> Layer
updateFullConnFilter fs lr dl
  | dl' == [] = FullConnLayer fs
  | otherwise = FullConnLayer fs'
  where
    dl' = catMaybes dl
    delta = mscale (lr / fromIntegral (length dl')) (msum $ strip dl')
    fs' = msub fs delta

strip :: [Layer] -> [FilterF]
strip [] = []
strip (FullConnLayer fs:ds) = fs:strip ds
strip (_:ds) = strip ds

max_stars_repo_name: eijian/deeplearning

max_issues_repo_name: eijian/deeplearning

max_forks_repo_name: eijian/deeplearning


text:
module Template where

import Data.Complex

-------------------------------------------------------------------------------
-- Utility functions

-- | Convert an 'Int' to a 'Double'
i2d :: Int -> Double
i2d = fromIntegral

-- | Convert a 'Double' to a 'Complex Double' with no imaginary part.
d2rc :: Double -> Complex Double
d2rc a = (a :+ 0.0)

-- | Convert a 'Double' to a 'Complex Double' with no real part.
d2ic :: Double -> Complex Double
d2ic a = (0.0 :+ a)

-- | A signal is a list of double precision samples
type Signal = [Double]

-- | Complete the definition of Fourier
data Fourier =  Fourier [Complex Double]

mkFourier :: [Complex Double] -> Fourier
mkFourier cl = (Fourier cl)

unFourier :: Fourier -> [Complex Double]
unFourier (Fourier cl) = cl

-------------------------------------------------------------------------------
-- Almost Equal

{-class AlmostEq a where
  (~=) :: a -> a -> Bool

infix 4 ~=

instance AlmostEq Double where
  a ~= b = abs(a-b) < 1e-14

instance AlmostEq (Complex Double) where
  a ~= b = magnitude (a-b) < 1e-14

instance AlmostEq a => AlmostEq [a] where
  xs ~= ys = (length xs == length ys) && equ xs ys where
    equ [] [] = True
    equ (x:xs) (y:ys) = x ~= y && equ xs ys

instance AlmostEq Fourier where
  xs ~= ys = equf xs ys where
    equf (Fourier a) (Fourier b) = a ~= b-}

-------------------------------------------------------------------------------
-- Naive DFT

-- | Naive discrete Fourier transform
dft :: Signal -> Fourier
dft [] = (Fourier [])
dft l = (Fourier (makehead : (maketail len))) where
  makehead = d2rc (sum l) --n=0
  len = (length l-1) --N-1
  maketail 0 = [] --n!=0
  maketail n = eoc : maketail (n-1) where
     eoc = d2rc (l !! 1) +  sum ([(d2rc (l !! k)) *  (exp (d2ic ((-2) * pi * (((i2d k) * (i2d n))/(i2d len))))) | k<-[1..len]])

makeTerm::Int->[Complex Double]->Int-> Double
makeTerm bn xs k= (1/(i2d bn)) * sum [realPart $ (xs !! n) * (exp $ d2ic $ pi * (i2d (2*n*k)) /(i2d bn)) | n <- [0..(bn-1)]]

-- | Naive inverse discrete Fourier transform
idft :: Fourier -> Signal
idft (Fourier list) = map (makeTerm (length list) list) [0..(length list -1)]

-------------------------------------------------------------------------------
-- Spectrum Analyser

rescale::[Double]->[Double]
rescale list = map (/maxi) list
  where maxi = maximum list

bar::Double->String
bar x = replicate (ceiling (x*40)) '#'

spectrumAnalyser :: Signal -> IO ()
spectrumAnalyser =   putStr . unlines . map bar . rescale . map magnitude . unFourier . dft

-------------------------------------------------------------------------------
-- Fast Fourier Transform

-- | Split a list in its even and oddly indexed elements.
scatter :: [a] -> ([a],[a])
scatter [] = ([],[])
scatter (x:y:xs) = ((x:l),(y:r)) where (l,r) = scatter xs

-- | Gather the Fourier coordinates into the correct order.
-- > gather [X0, X4, X1, X5, X2, X6, X3, X7] = [X0,X1,X2,X3,X4,X5,X6,X7]
gather :: [a] -> [a]
gather list = l ++ r where (l,r) = scatter list

-- | Combine even and odd Fourier coefficients with the twiddle factor.
twiddle :: Int -> Int -> Complex Double -> Complex Double -> [Complex Double]
twiddle m n c1 c2 = [c1 + number * c2 ,c1 - number * c2 ] where
  number = exp (d2ic ((-2)*pi*((i2d n)/(i2d m))))

-- | The butterfly step: construct Fourier coefficients from the
-- Fourier coefficients of the even an odd elements.
butterfly :: Int -> Fourier -> Fourier -> [Complex Double]
butterfly bn (Fourier even) (Fourier odd)= concat $ zipWith3 (twiddle bn)  [0..] even odd

-- | Cooley-Tukey Fast Fourier Transform algorithm.
fft :: Signal -> Fourier
fft xs
  | odd (length xs) = dft xs
  | otherwise = let (evens,odds) = scatter xs
                in Fourier $ gather $ butterfly (length xs) (fft evens) (fft odds)

max_stars_repo_name: SeraphinaLiang/LectureHaskell

max_issues_repo_name: SeraphinaLiang/LectureHaskell

max_forks_repo_name: SeraphinaLiang/LectureHaskell


text:
{-# LANGUAGE GADTs
           , DataKinds
           , TypeApplications
           , ScopedTypeVariables
           , ExplicitForAll
           , AllowAmbiguousTypes
           , KindSignatures
           , TypeOperators
           , RankNTypes
           , GeneralizedNewtypeDeriving
           , StrictData #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}

module Lib where

import Numeric.LinearAlgebra hiding (Vector, Matrix)
import qualified Numeric.LinearAlgebra as LA
import qualified Data.Vector.Storable as VS
import GHC.TypeLits
import Data.Proxy
import System.Random
import Control.Monad
import Numeric.AD

newtype Vector (n :: Nat) = Vector { _vec :: VS.Vector Float }
  deriving (Show, Floating)

instance KnownNat n => Num (Vector n) where
  (Vector v1) + (Vector v2) = Vector $ v1 + v2
  (Vector v1) - (Vector v2) = Vector $ v1 - v2
  (Vector v1) * (Vector v2) = Vector $ v1 * v2
  abs                       = vectorMap abs
  fromInteger               = vecFromList . (:[]) . fromInteger
  negate                    = vectorMap negate
  signum                    = vectorMap signum

instance KnownNat n => Fractional (Vector n) where
  (Vector v1) / (Vector v2) = Vector $ v1 / v2
  fromRational              = vecFromList . (:[]) . fromRational

natToInt :: forall n. KnownNat n => Int
natToInt = fromIntegral . natVal $ Proxy @n

vecFromList :: forall n. KnownNat n => [Float] -> Vector n
vecFromList = Vector . VS.fromList . take (natToInt @n) . cycle

vecLength :: forall n. KnownNat n => Vector n -> Int
vecLength _ = natToInt @n

vecTail :: forall n. KnownNat n => Vector (n+1) -> Vector n
vecTail (Vector v) = Vector $ VS.tail v

newtype Matrix (n :: Nat) (m :: Nat) = Matrix { _mat :: LA.Matrix Float }
  deriving Show

matrixWidth :: forall n m. (KnownNat n, KnownNat m) => Matrix n m -> Int
matrixWidth _ = natToInt @n

matrixHeight :: forall n m. (KnownNat n, KnownNat m) => Matrix n m -> Int
matrixHeight _ = natToInt @m

matFromList :: forall n m. (KnownNat n, KnownNat m) => [Float] -> Matrix n m
matFromList = Matrix . (h><w) . cycle
  where w = natToInt @n
        h = natToInt @m

randomMatrix :: forall n m. (KnownNat n, KnownNat m) => IO (Matrix n m)
randomMatrix = matFromList <$> replicateM (w*h) (randomRIO (0.001, 0.1))
    where w = natToInt @n
          h = natToInt @m

data Network :: Nat -> Nat -> * where
  Layer :: forall m n. (KnownNat n, KnownNat m) =>
                       Matrix (n+1) m ->
                       (forall a. Floating a => a -> a) ->
                       Network n m
  Composition :: (KnownNat n, KnownNat m, KnownNat k) =>
                 Network m k ->
                 Network n m ->
                 Network n k

conductSignal :: forall n m. (KnownNat n, KnownNat m) =>
                               Network n m ->
                               Vector n ->
                               Vector m
conductSignal (Layer mat f) = vectorMap f . mvMul mat . vectorCons 1
conductSignal (Composition network2 network1) = conductSignal network2 . conductSignal network1

vectorMap :: (Float -> Float) -> Vector n -> Vector n
vectorMap f (Vector v) = Vector $ VS.map f v

mvMul :: Matrix n m -> Vector n -> Vector m
mvMul (Matrix m) (Vector v) = Vector $ m #> v

vectorCons :: KnownNat n => Float -> Vector n -> Vector (n+1)
vectorCons n (Vector v) = Vector $ VS.cons n v

randomLayer :: forall n m. (KnownNat n, KnownNat m) => (forall a. Floating a => a -> a) -> IO (Network n m)
randomLayer f = (`Layer` f) <$> randomMatrix

scaleMatrix :: (KnownNat n, KnownNat m) => Float -> Matrix m n -> Matrix m n
scaleMatrix s (Matrix m) = Matrix $ LA.scale s m

transposeMat :: (KnownNat n, KnownNat m) => Matrix n m -> Matrix m n
transposeMat (Matrix mat) = Matrix $ LA.tr mat

outerProduct :: (KnownNat n, KnownNat m) => Vector n -> Vector m -> Matrix n m
outerProduct (Vector v1) (Vector v2) = Matrix $ v2 `outer` v1

vvMul :: (KnownNat n) => Vector n -> Vector n -> Vector n
vvMul (Vector v1) (Vector v2) = Vector $ v1 * v2

matrixSub :: (KnownNat m, KnownNat n) => Matrix n m -> Matrix n m -> Matrix n m
matrixSub (Matrix m1) (Matrix m2) = Matrix $ m1 - m2

useTrainingExample :: forall m n. (KnownNat n, KnownNat m) =>
                        Vector m ->
                        Network n m ->
                        Vector n ->
                        (Network n m, Vector n)
useTrainingExample diffLoss (Layer weights activation) input = (trainedLayer, inputGradient)
  where
    trainedLayer = Layer newWeights activation
    newWeights :: Matrix (n+1) m
    newWeights = matrixSub weights $ scaleMatrix 0.01 weightGradient
    consedInput = vectorCons 1 input
    preActivationOutput = mvMul weights consedInput
    output = vectorMap activation preActivationOutput
    diffActivation = vectorMap (diff activation) preActivationOutput
    weightGradient = outerProduct consedInput (diffLoss * diffActivation)
    inputGradient  = vecTail $ mvMul (transposeMat weights) (diffLoss * diffActivation)
useTrainingExample diffLoss (Composition network2 network1) input =
    (Composition trainedNet2 trainedNet1, inputGradient)
  where (trainedNet2,outputGrad1)   = useTrainingExample diffLoss network2 input2
        input2                      = conductSignal network1 input
        (trainedNet1,inputGradient) = useTrainingExample outputGrad1 network1 input

showExample :: (KnownNat n, KnownNat m) =>
                (Vector m -> Vector m -> Vector m) ->
                Network n m ->
                (Vector n, Vector m) ->
                Network n m
showExample lossGrad network (input,expected) = fst $ useTrainingExample diffLoss network input
  where diffLoss = lossGrad expected output
        output   = conductSignal network input

sigmoid input = 1/(1+exp(-input))
                -- (1/) . (1+) . exp . negate

epsilon :: Floating a => a
epsilon = 0.0000001

logLoss :: Floating a => [a] -> a
logLoss [expected, output] = negate $ expected * log (output + epsilon) + (1 - expected) * log (1 - output + epsilon)

gradLogLoss expected output = grad logLoss [expected, output] !! 1

-- "mnist/train-set"

getTrainSet :: IO [(Vector 784, Vector 10)]
getTrainSet =
  map (\(a,b) -> (vecFromList a, toOneHot b)) .
  map (read @([Float],Int)) .
  lines <$> readFile  "mnist/train-set"
  where
    toOneHot n = vecFromList $ map (\x -> if x == n then 1 else 0) [0..]

eval network (input,expected) = do
  putStrLn "EXPECTED:"
  print expected
  putStrLn "PREDICTED:"
  print $ conductSignal network input

trainTest k n = do
  trainSet <- take n <$> getTrainSet
  untrainedLayer1 <- randomLayer sigmoid
  untrainedLayer2 <- randomLayer @50 @10 sigmoid
  let untrainedNetwork = Composition untrainedLayer2 untrainedLayer1
      trainedNetwork = foldl (showExample gradLogLoss) untrainedNetwork . take (k*n) $ cycle trainSet
      evalSet = take 10 trainSet
  mapM_ (eval trainedNetwork) evalSet

{-

  diff (loss . activation . mvMul(weights, input))

   => diff loss (activation . mvMul (weights,input)) * diff (activation . mvMul (weights,input))

   => diff loss (output) * diff activation (mvMul (weights,input)) * diff mvMul(weights,input)
-}

max_stars_repo_name: Antystenes/MPNN

max_issues_repo_name: Antystenes/MPNN

max_forks_repo_name: Antystenes/MPNN