text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE ScopedTypeVariables #-}
module Data.Lisp (Number(..), SourceRange(..), Lisp(..), parseLisp,
                  parseLispFile, parseLispExpr, showLispPos, CharParser,
                  lispParser) where
import qualified Data.Text as Text
import qualified Data.Text.IO as Text
import Data.Text(Text)
import Control.Applicative hiding (some, many)
import Data.Ratio
import Data.Char
import Control.Monad
import Data.Complex
import Data.Void
import Text.Megaparsec
import Text.Megaparsec.Char

data Number =
  Integer Integer |
  SingleFloat Float |
  DoubleFloat Double |
  NumRatio (Ratio Integer) |
  ComplexDouble (Complex Double)

replaceChar :: Char -> Char -> String -> String
replaceChar _ _ [] = []
replaceChar from to (c:cs)
  | c == from = to:replaceChar from to cs
  | otherwise = c:replaceChar from to cs

specialChars :: String
specialChars = "()#\"\\,'`| ;"

instance Show Number where
  show (Integer i) = show i
  show (SingleFloat f) = replaceChar 'e' 's' $ show f
  show (DoubleFloat f)
    | 'd' `elem` str = str
    | otherwise = str ++ "d0"
    where str = replaceChar 'e' 'd' $ show f
  show (NumRatio r) = show (numerator r) ++ "/" ++ show (denominator r)
  show (ComplexDouble (a :+ b)) = "#(" ++ show a ++ " " ++ show b ++ ")"

-- | A position range in the Lisp source file or string.
data SourceRange = SourceRange
  { sourceFrom :: SourcePos
  , sourceTo :: SourcePos
  }

instance Show SourceRange where
  show (SourceRange from to) =
    "<" ++ showPos from ++ "," ++ showPos to ++ ">"
    where showPos (SourcePos _ l c) =
            show (unPos l) ++ ":" ++ show (unPos c)

data Lisp =
  LispString Text SourceRange |
  LispNumber Number SourceRange |
  LispSymbol Text SourceRange |
  LispVector [Lisp] SourceRange |
  LispList [Lisp] SourceRange |
  LispDotList [Lisp] Lisp SourceRange

instance Show Lisp where
  show (LispString t _) = show t
  show (LispNumber n _) = show n
  show (LispSymbol s _)
    | Text.null s = "||"
    | Text.any (`elem` specialChars) s = '|': Text.unpack s ++ "|"
    | otherwise = Text.unpack s

  show (LispVector l _) =
    "#(" ++ unwords (map show l) ++ ")"
  show (LispList l _) = "(" ++ unwords (map show l) ++ ")"
  show (LispDotList l e _) =
    "(" ++ unwords (map show l) ++ " . " ++ show e ++ ")"

-- | show the lisp with position info
showLispPos :: Lisp -> String
showLispPos (LispString t p) = show t ++ show p
showLispPos (LispNumber n p) = show n ++ show p
showLispPos (LispSymbol s p)
  | Text.null s = "||" ++ show p
  | Text.any (`elem` specialChars) s = '|': Text.unpack s ++ "|" ++ show p
  | otherwise = Text.unpack s ++ show p
showLispPos (LispVector l p) =
  "#(" ++ unwords (map showLispPos l) ++ ")" ++ show p
showLispPos (LispList l p) = "(" ++ unwords (map showLispPos l) ++ ")" ++ show p
showLispPos (LispDotList l e p) =
  "(" ++ unwords (map showLispPos l) ++ " . " ++ showLispPos e ++ ")" ++ show p

instance Read Lisp where
  readsPrec _ input =
    case runParser' (whiteSpace >> withSourceRange lispExprP) $
         State input 0 (PosState input 0 (initialPos "read") (mkPos 0) []) [] of
      (_, Left _) -> []
      (rest, Right r) -> [(r, stateInput rest)]

dummyRange :: SourceRange
dummyRange = SourceRange (initialPos "dummy") (initialPos "dummy")

-- | A megaparsec parser that has characters as tokens.
type CharParser t a = (Stream t, Token t ~ Char)
                    => Parsec Void t a

-- | A megaparsec parser for lisp expressions
lispParser :: CharParser t Lisp
lispParser = withSourceRange lispExprP

-- | Parse a lisp file
parseLispFile :: String -> IO (Either (ParseErrorBundle Text Void) [Lisp])
parseLispFile file =
  runParser (many lispParser <* whiteSpace <* eof) file
  <$> Text.readFile file

-- | @parse source text@: parse the text into a list of lisp
-- expressions.  Source is used for the error messages, and in the
-- `SourceRanges`.
parseLisp :: String -> Text -> Either (ParseErrorBundle Text Void) [Lisp]
parseLisp = runParser (sepEndBy lispParser whiteSpace <* eof)

-- | parse a single expression
parseLispExpr :: String -> Text -> Either (ParseErrorBundle Text Void) Lisp
parseLispExpr = runParser lispParser

signP :: CharParser t String
signP = option "" $ ("" <$ char '+') <|> ("-" <$ char '-')

withSourceRange :: CharParser t (SourceRange -> a) -> CharParser t a
withSourceRange p = do
  startRange <- getSourcePos
  mkParser <- p
  endRange <- getSourcePos
  pure $ mkParser $ SourceRange startRange endRange

-- numbers not starting with #
numP :: CharParser t (SourceRange -> Lisp)
numP = label "number" $ do
  sign <- signP
  -- 'try' the dot, because it could be a single dot, and then we need
  -- to backtrack
  let decimalP :: CharParser t String
      decimalP = some digitChar

      -- number starting with number
      numNumP :: CharParser t (SourceRange -> Lisp)
      numNumP = do
        decimal <- decimalP
        choice [ ratioP decimal
               , try (floatP decimal)
               , do _ <- optional (char '.')
                    pure $ LispNumber $ Integer $ read (sign++decimal)
               ]

      ratioP :: String -> CharParser t (SourceRange -> Lisp)
      ratioP d = do
        _ <- char '/'
        denom <- decimalP
        pure $ LispNumber $ NumRatio $ read (sign++d) % read denom

      floatP :: String -> CharParser t (SourceRange -> Lisp)
      floatP d =
        exptP sign d "0"
        <|> do _ <- char '.'
               exptP sign d "0" <|> do
                 fract <- decimalP
                 exptP sign d fract <|>
                   do pure $ LispNumber $ DoubleFloat $
                        read (sign++d ++ "." ++ fract)
      dotNumP = do
        _ <- char '.'
        fract <- decimalP
        exptP sign "0" fract <|>
          pure (LispNumber $ DoubleFloat $ read $ sign++"0." ++ fract)

  (numNumP <|> try dotNumP) <* notFollowedBy identifierBlocksP

exptP :: String -> String -> String -> CharParser t (SourceRange -> Lisp)
exptP sign num fract = do
  -- e would be context dependend, but I am defaulting it to Double here
  e <- oneOf ("esd" :: String)
  eSign <- option '+' $ char '+' <|> char '-'
  expt <- some digitChar
  let toFloat :: (Read a, Num a) => a
      toFloat = read $ sign ++ num ++ "." ++ fract ++ "e" ++ eSign:expt
  pure $ LispNumber $ case e of
    's' -> SingleFloat toFloat
    _   -> DoubleFloat toFloat

quoteAnyChar :: CharParser t Char
quoteAnyChar = char '\\' >> anySingle

stringP :: CharParser t (SourceRange -> Lisp)
stringP =
  label "string" $ do
  str <- between (char '"') (char '"') $
         Text.pack <$> many (quoteAnyChar <|> noneOf ("\\\"" :: String))
  pure $ LispString str

identifierP :: CharParser t (SourceRange -> Lisp)
identifierP =
  label "identifier" $ do
  str <- fmap Text.pack $ (++) <$> (firstBlock <|> quotedBlockP) <*> moreBlocksP
  if str == "."
    then fail ("dot" :: String)
    else pure $ LispSymbol str

  where firstBlock :: CharParser t String
        firstBlock = (:) <$> (notSpecial <|> quoteAnyChar) <*> many blockCharP

        moreBlocksP :: CharParser t String
        moreBlocksP = concat <$> many (some blockCharP <|> quotedBlockP)

quotedBlockP :: CharParser t String
quotedBlockP = between (char '|') (char '|') $
               many (noneOf ("|\\" :: String) <|> quoteAnyChar)

notSpecial :: CharParser t Char
notSpecial = toUpper <$> noneOf specialChars

blockCharP :: CharParser t Char
blockCharP = notSpecial <|> char '#' <|> char '`' <|> quoteAnyChar

identifierBlocksP :: CharParser t String
identifierBlocksP = concat <$> some (some blockCharP <|> quotedBlockP)

lispExprP :: CharParser t (SourceRange -> Lisp)
lispExprP = choice [ stringP
                   , listP
                   , try numP
                   , try identifierP
                   , quoteP
                   , readersP
                   ]

listP :: CharParser t (SourceRange -> Lisp)
listP =
  label "list" $
  between (char '(') (char ')') $ do
  elems <- lispParser `sepEndBy` whiteSpace
  dotElem <- optional $
    char '.' *> whiteSpace *>
    lispParser <* whiteSpace
  pure $ case dotElem of
    Nothing -> LispList elems
    Just (LispList l _) ->  LispList $ elems ++ l
    Just (LispDotList l el _) -> LispDotList (elems ++ l) el
    Just el -> LispDotList elems el

commentP :: CharParser t ()
commentP =
  label "comment" $
  char ';' >>
  many (noneOf ("\r\n" :: String)) >>
  (void eol <|> eof)

whiteSpace :: CharParser t ()
whiteSpace = () <$ many (space1 <|> commentP)

quoteSymbol :: String -> Text -> CharParser t (SourceRange -> Lisp)
quoteSymbol symbol quoteExpansion  = do
  foldMap (void <$> char) symbol
  pure $ \(SourceRange from _) ->
    LispSymbol quoteExpansion $ SourceRange from $
    from {sourceColumn = mkPos $ length symbol + unPos (sourceColumn from)}

quoteP :: CharParser t (SourceRange -> Lisp)
quoteP = do
  quote <- (quoteSymbol "'" "quote" <|>
            try (quoteSymbol ",@" "unquote-splicing") <|>
            quoteSymbol "," "unquote" <|>
            quoteSymbol "`" "quasiquote")
           <* whiteSpace
  expr <- lispParser
  pure $ \range -> LispList [quote range, expr] range

readersP :: CharParser t (SourceRange -> Lisp)
readersP = do
  _ <- char '#'
  (vectorReaderP <|> complexReaderP) <|>
    (octalReaderP <|> hexReaderP <|> binaryReaderP)
    <* notFollowedBy identifierBlocksP

vectorReaderP :: CharParser t (SourceRange -> Lisp)
vectorReaderP =
  between (char '(') (char ')') $
  LispVector <$> (lispParser `sepEndBy` whiteSpace)

octalReaderP :: CharParser t (SourceRange -> Lisp)
octalReaderP = do
  _ <- char 'o' <|> char 'O'
  sign <- signP
  digits <- some octDigitChar
  pure $ LispNumber $ Integer $ read $ sign ++ "0o" ++ digits

binaryReaderP :: CharParser t (SourceRange -> Lisp)
binaryReaderP = do
  _ <- char 'b' <|> char 'B'
  sign <- signP
  digits <- some binDigitChar
  let digitSum = foldl (\tot dig -> tot*2 + if dig == '1' then 1 else 0)
                 0 digits
      signedSum | sign == "-" = negate digitSum
                | otherwise = digitSum
  pure $ LispNumber $ Integer signedSum

hexReaderP :: CharParser t (SourceRange -> Lisp)
hexReaderP = do
  _ <- char 'x' <|> char 'X'
  sign <- signP
  digits <- some hexDigitChar
  pure $ LispNumber $ Integer $ read $ sign ++ "0x" ++ digits

convertToDouble :: Number -> Double
convertToDouble l = case l of
  Integer i -> realToFrac i
  SingleFloat f -> realToFrac f
  DoubleFloat f -> realToFrac f
  NumRatio r -> realToFrac r
  ComplexDouble _ -> error "convertToDouble"

complexReaderP :: CharParser t (SourceRange -> Lisp)
complexReaderP = do
  _ <- char 'c' <|> char 'C'
  between (char '(') (char ')') $ do
    _ <- many whiteSpace
    LispNumber rl _ <- ($ dummyRange) <$> numP
    _ <- some whiteSpace
    LispNumber imag _ <- ($ dummyRange) <$> numP
    _ <- many whiteSpace
    pure $ LispNumber $ ComplexDouble $
      convertToDouble rl :+ convertToDouble imag

max_stars_repo_name: kuribas/megalisp

max_issues_repo_name: kuribas/megalisp

max_forks_repo_name: kuribas/megalisp


text:
{-# LANGUAGE Rank2Types #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}
-----------------------------------------------------------------------------
-- |
-- Copyright   :  (c) Edward Kmett 2010-2014
-- License     :  BSD3
-- Maintainer  :  ekmett@gmail.com
-- Stability   :  experimental
-- Portability :  GHC only
--
-----------------------------------------------------------------------------

module Numeric.AD.Newton
  (
  -- * Newton's Method (Forward AD)
    findZero
  , inverse
  , fixedPoint
  , extremum
  -- * Gradient Ascent/Descent (Reverse AD)
  , gradientDescent
  , gradientAscent
  , conjugateGradientDescent
  , conjugateGradientAscent
  ) where

import Data.Foldable (all, sum)
import Data.Reflection (Reifies)
import Data.Traversable
import Numeric.AD.Internal.Combinators
import Numeric.AD.Internal.Forward (Forward)
import Numeric.AD.Internal.On
import Numeric.AD.Internal.Or
import Numeric.AD.Internal.Reverse (Reverse, Tape)
import Numeric.AD.Internal.Type (AD(..))
import Numeric.AD.Mode
import Numeric.AD.Mode.Reverse as Reverse (gradWith')
import Numeric.AD.Rank1.Kahn as Kahn (Kahn, grad)
import qualified Numeric.AD.Rank1.Newton as Rank1
import Prelude hiding (all, mapM, sum)

-- $setup
-- >>> import Data.Complex

-- | The 'findZero' function finds a zero of a scalar function using
-- Newton's method; its output is a stream of increasingly accurate
-- results.  (Modulo the usual caveats.) If the stream becomes constant
-- ("it converges"), no further elements are returned.
--
-- Examples:
--
-- >>> take 10 $ findZero (\x->x^2-4) 1
-- [1.0,2.5,2.05,2.000609756097561,2.0000000929222947,2.000000000000002,2.0]
--
-- >>> last $ take 10 $ findZero ((+1).(^2)) (1 :+ 1)
-- 0.0 :+ 1.0
findZero :: (Fractional a, Eq a) => (forall s. AD s (Forward a) -> AD s (Forward a)) -> a -> [a]
findZero f = Rank1.findZero (runAD.f.AD)
{-# INLINE findZero #-}

-- | The 'inverse' function inverts a scalar function using
-- Newton's method; its output is a stream of increasingly accurate
-- results.  (Modulo the usual caveats.) If the stream becomes
-- constant ("it converges"), no further elements are returned.
--
-- Example:
--
-- >>> last $ take 10 $ inverse sqrt 1 (sqrt 10)
-- 10.0
inverse :: (Fractional a, Eq a) => (forall s. AD s (Forward a) -> AD s (Forward a)) -> a -> a -> [a]
inverse f = Rank1.inverse (runAD.f.AD)
{-# INLINE inverse  #-}

-- | The 'fixedPoint' function find a fixedpoint of a scalar
-- function using Newton's method; its output is a stream of
-- increasingly accurate results.  (Modulo the usual caveats.)
--
-- If the stream becomes constant ("it converges"), no further
-- elements are returned.
--
-- >>> last $ take 10 $ fixedPoint cos 1
-- 0.7390851332151607
fixedPoint :: (Fractional a, Eq a) => (forall s. AD s (Forward a) -> AD s (Forward a)) -> a -> [a]
fixedPoint f = Rank1.fixedPoint (runAD.f.AD)
{-# INLINE fixedPoint #-}

-- | The 'extremum' function finds an extremum of a scalar
-- function using Newton's method; produces a stream of increasingly
-- accurate results.  (Modulo the usual caveats.) If the stream
-- becomes constant ("it converges"), no further elements are returned.
--
-- >>> last $ take 10 $ extremum cos 1
-- 0.0
extremum :: (Fractional a, Eq a) => (forall s. AD s (On (Forward (Forward a))) -> AD s (On (Forward (Forward a)))) -> a -> [a]
extremum f = Rank1.extremum (runAD.f.AD)
{-# INLINE extremum #-}

-- | The 'gradientDescent' function performs a multivariate
-- optimization, based on the naive-gradient-descent in the file
-- @stalingrad\/examples\/flow-tests\/pre-saddle-1a.vlad@ from the
-- VLAD compiler Stalingrad sources.  Its output is a stream of
-- increasingly accurate results.  (Modulo the usual caveats.)
--
-- It uses reverse mode automatic differentiation to compute the gradient.
gradientDescent :: (Traversable f, Fractional a, Ord a) => (forall s. Reifies s Tape => f (Reverse s a) -> Reverse s a) -> f a -> [f a]
gradientDescent f x0 = go x0 fx0 xgx0 0.1 (0 :: Int)
  where
    (fx0, xgx0) = Reverse.gradWith' (,) f x0
    go x fx xgx !eta !i
      | eta == 0     = [] -- step size is 0
      | fx1 > fx     = go x fx xgx (eta/2) 0 -- we stepped too far
      | zeroGrad xgx = [] -- gradient is 0
      | otherwise    = x1 : if i == 10
                            then go x1 fx1 xgx1 (eta*2) 0
                            else go x1 fx1 xgx1 eta (i+1)
      where
        zeroGrad = all (\(_,g) -> g == 0)
        x1 = fmap (\(xi,gxi) -> xi - eta * gxi) xgx
        (fx1, xgx1) = Reverse.gradWith' (,) f x1
{-# INLINE gradientDescent #-}

-- | Perform a gradient descent using reverse mode automatic differentiation to compute the gradient.
gradientAscent :: (Traversable f, Fractional a, Ord a) => (forall s. Reifies s Tape => f (Reverse s a) -> Reverse s a) -> f a -> [f a]
gradientAscent f = gradientDescent (negate . f)
{-# INLINE gradientAscent #-}

-- | Perform a conjugate gradient descent using reverse mode automatic differentiation to compute the gradient, and using forward-on-forward mode for computing extrema.
--
-- >>> let sq x = x * x
-- >>> let rosenbrock [x,y] = sq (1 - x) + 100 * sq (y - sq x)
-- >>> rosenbrock [0,0]
-- 1
-- >>> rosenbrock (conjugateGradientDescent rosenbrock [0, 0] !! 5) < 0.1
-- True
conjugateGradientDescent
  :: (Traversable f, Ord a, Fractional a)
  => (forall s. Chosen s => f (Or s (On (Forward (Forward a))) (Kahn a)) -> Or s (On (Forward (Forward a))) (Kahn a))
  -> f a -> [f a]
conjugateGradientDescent f = conjugateGradientAscent (negate . f)
{-# INLINE conjugateGradientDescent #-}

lfu :: Functor f => (f (Or F a b) -> Or F a b) -> f a -> a
lfu f = runL . f . fmap L

rfu :: Functor f => (f (Or T a b) -> Or T a b) -> f b -> b
rfu f = runR . f . fmap R

-- | Perform a conjugate gradient ascent using reverse mode automatic differentiation to compute the gradient.
conjugateGradientAscent
  :: (Traversable f, Ord a, Fractional a)
  => (forall s. Chosen s => f (Or s (On (Forward (Forward a))) (Kahn a)) -> Or s (On (Forward (Forward a))) (Kahn a))
  -> f a -> [f a]
conjugateGradientAscent f x0 = takeWhile (all (\a -> a == a)) (go x0 d0 d0 delta0)
  where
    dot x y = sum $ zipWithT (*) x y
    d0 = Kahn.grad (rfu f) x0
    delta0 = dot d0 d0
    go xi _ri di deltai = xi : go xi1 ri1 di1 deltai1
      where
        ai = last $ take 20 $ Rank1.extremum (\a -> lfu f $ zipWithT (\x d -> auto x + a * auto d) xi di) 0
        xi1 = zipWithT (\x d -> x + ai*d) xi di
        ri1 = Kahn.grad (rfu f) xi1
        deltai1 = dot ri1 ri1
        bi1 = deltai1 / deltai
        di1 = zipWithT (\r d -> r + bi1 * d) ri1 di
{-# INLINE conjugateGradientAscent #-}

ext: hs

max_stars_repo_name: silky/ad

max_issues_repo_name: silky/ad

max_forks_repo_name: silky/ad


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}

-- |
-- Module      :  Test.Factorization
-- Copyright   :  (c) 2016-2020 Drexel University
-- License     :  BSD-style
-- Maintainer  :  mainland@drexel.edu

module Test.Factorization (
    factorizationTests,
    f2Test,
    f4Test,
    f8Test,
    cooleyTukeyDITTest,
    cooleyTukeyDIFTest,
    goodThomasTest,
    raderTest,
    bluesteinTest,
    winogradSmallTest,
    ditTest,
    difTest,
    splitRadixTest,
    splitRadix8Test,
    conjPairSplitRadixTest,
    improvedSplitRadixTest,
    opcountSearchTest
  ) where

import Data.Complex
import Data.Modular
import Data.Proxy (Proxy(..))
import Test.HUnit ((@?=))
import Test.Hspec

import Spiral.Array (M,
                     Matrix)
import qualified Spiral.Array as A
import Spiral.Convolution
import Spiral.Driver
import Spiral.Exp
import Spiral.FFT.Bluestein
import Spiral.FFT.CooleyTukey
import Spiral.FFT.GoodThomas
import Spiral.FFT.Rader
import Spiral.FFT.Winograd
import Spiral.NumberTheory (factors)
import Spiral.RootOfUnity
import Spiral.SPL
import Spiral.Search.FFTBreakdowns (getCycs)
import Spiral.Search.OpCount

factorizationTests :: Spec
factorizationTests = do
    f2Test
    f4Test
    f8Test
    complexFactorizationTests
    describe "Modular DFT" $ do
        describe "ℤ/17" $
            sequence_ [ditTest (Proxy :: Proxy (ℤ/17)) i | i <- [1..4::Int]]
        describe "ℤ/2013265921" $
            sequence_ [ditTest (Proxy :: Proxy (ℤ/2013265921)) i | i <- [1..7::Int]]

complexFactorizationTests :: Spec
complexFactorizationTests = do
    describe "Cooley-Tukey" $ do
        cooleyTukeyDITTest p 5 7
        cooleyTukeyDIFTest p 5 7
    describe "Good-Thomas" $
        goodThomasTest p 5 7
    describe "Rader" $ do
        mapM_ (raderTest p "Rader" rader) [7, 23]
        mapM_ (raderTest p "RaderI" (\n w -> raderI n w (ConvolutionTheorem (n-1)))) [7, 17, 23]
        mapM_ (raderTest p "RaderII" raderII) [7, 17, 23]
        mapM_ (raderTest p "RaderIII" raderIII) [7, 17, 23]
        mapM_ (raderTest p "RaderIV" raderIV) [7, 17, 23]
        mapM_ (raderTest p "raderLuIII" raderLuIII) [7, 17, 23]
    describe "Bluestein" $ do
        mapM_ (bluesteinTest p "Bluestein" bluestein) [(3, 6), (4, 8), (4, 9), (9, 18)]
        mapM_ (bluesteinTest p "Bluestein'" bluestein') [(3, 6), (4, 8), (4, 9), (9, 18)]
    describe "Winograd Small" $ do
        winogradSmallTest 3 (ConvolutionTheorem 2)
        winogradSmallTest 5 (ConvolutionTheorem 4)
        winogradSmallTest 7 (SplitNesting [3, 2] [(Winograd 3 [Standard 1, Standard 2]), (Winograd 2 [Standard 1, Standard 1])])
        winogradSmallTest 7 (AgarwalCooley 3 2 (Winograd 3 [Standard 1, Lift 2 6 (Tensor [3,2] [Standard 3,Standard 2])]) (Winograd 2 [Standard 1,Standard 1]))
    describe "Winograd Large" $
        mapM_ winogradLargeTest [(2,3), (5, 2), (2,7), (3,4), (4,3), (5,4), (3,8), (5,8), (8,3), (3,5), (5,3)]
    describe "Winograd 2^n" $
        mapM_ winogradTwoTest [2^i | i <- [1..5::Int]]
    describe "Winograd Power" $ do
        winogradPowerTest 3 2
        winogradSquareTest 3
        winogradSquareTest 5
    describe "DIT" $
        sequence_ [ditTest p n | n <- [1..7]]
    describe "DIF" $
        sequence_ [difTest p n | n <- [1..7]]
    describe "Split Radix" $
        sequence_ [splitRadixTest p n | n <- [1..3]]
    describe "Split Radix 8" $
        sequence_ [splitRadix8Test p n | n <- [1..2]]
    describe "Conjugate Pair Split Radix" $
        sequence_ [conjPairSplitRadixTest p n | n <- [1..3]]
    describe "Improved Split Radix" $
        sequence_ [improvedSplitRadixTest p n | n <- [1..3]]
  where
    p :: Proxy (Complex Double)
    p = Proxy

-- $F_2$
f2Test :: Spec
f2Test = it "F_2" $ toMatrix (dit 2) @?= f2
  where
    f2 :: Matrix M (Exp (Complex Double))
    f2 = A.matrix [[1,  1],
                   [1, -1]]

-- $F_4$ calculated per "SPL: A Language and Compiler for DSP Algorithms"
f4Test :: Spec
f4Test = it "F_4" $ toMatrix (dit 4) @?= f4
  where
    f4 :: Matrix M (Exp (Complex Double))
    f4 = A.matrix [[1,  1,  1,  1],
                   [1, -i, -1,  i],
                   [1, -1,  1, -1],
                   [1,  i, -1, -i]]
      where
        i :: Exp (Complex Double)
        i = complexE (0 :+ 1)

-- $F_8$ calculated per "SPL: A Language and Compiler for DSP Algorithms"
-- See also:
--   https://en.wikipedia.org/wiki/DFT_matrix
f8Test :: Spec
f8Test = it "F_8" $ toMatrix (dit 8) @?= f8
  where
    f8 :: Matrix M (Exp (Complex Double))
    f8 = A.matrix [[1,  1,  1,  1, 1, 1, 1, 1],
                   [1,  w, -i,  -i*w, -1, -w, i, i*w],
                   [1, -i, -1, i, 1, -i, -1, i],
                   [1, -i*w, i, w, -1, i*w, -i, -w],
                   [1, -1, 1, -1, 1, -1, 1, -1],
                   [1, -w, -i, i*w, -1, w, i, -i*w],
                   [1, i, -1, -i, 1, i, -1, -i],
                   [1, i*w, i, -w, -1, -i*w, -i, w]]
      where
        i = complexE (0 :+ 1)

        w = omega 8

-- Test Cooley-Tukey DIT
cooleyTukeyDITTest :: forall p . RootOfUnity (Exp p)
                   => Proxy p
                   -> Int
                   -> Int
                   -> Spec
cooleyTukeyDITTest _ r s =
    it ("CooleyTukeyDIT(" ++ show r ++ "," ++ show s ++ ")") $
    toMatrix (cooleyTukeyDIT r s w) @?= toMatrix (F rs w)
  where
    w :: Exp p
    w = omega rs

    rs :: Int
    rs = r*s

-- Test Cooley-Tukey DIF
cooleyTukeyDIFTest :: forall p . RootOfUnity (Exp p)
                   => Proxy p
                   -> Int
                   -> Int
                   -> Spec
cooleyTukeyDIFTest _ r s =
    it ("CooleyTukeyDIF(" ++ show r ++ "," ++ show s ++ ")") $
    toMatrix (cooleyTukeyDIF r s w) @?= toMatrix (F rs w)
  where
    w :: Exp p
    w = omega rs

    rs :: Int
    rs = r*s

-- Test Good-Thomas with factors 5 and 7
goodThomasTest :: forall p . RootOfUnity (Exp p)
               => Proxy p
               -> Int
               -> Int
               -> Spec
goodThomasTest _ r s =
    it ("GoodThomas(" ++ show r ++ "," ++ show s ++ ")") $
    toMatrix (goodThomas r s w) @?= toMatrix (F rs w)
  where
    w :: Exp p
    w = omega rs

    rs :: Int
    rs = r*s

-- Test Rader for given prime
raderTest :: forall p . RootOfUnity (Exp p)
          => Proxy p
          -> String
          -> (Int -> Exp p -> SPL (Exp p))
          -> Int
          -> Spec
raderTest _ desc rader n = it (desc ++ "(" ++ show n ++ ")") $
    toMatrix (rader (fromIntegral n) w) @?= toMatrix (F n w)
  where
    w :: Exp p
    w = omega n

-- Test Bluestein for given n and m
bluesteinTest :: forall p . RootOfUnity (Exp p)
              => Proxy p
              -> String
              -> (Int -> Int -> Exp p -> SPL (Exp p))
              -> (Int, Int)
              -> Spec
bluesteinTest _ desc bluestein (n, m) = it (desc ++ "(" ++ show n ++ "," ++ show m ++ ")") $
    toMatrix (bluestein n m w) @?= toMatrix (DFT n)
  where
    w :: Exp p
    w = omega (2*n)

winogradSmallTest :: Int
                  -> CyclicConvolution (Exp (Complex Double))
                  -> Spec
winogradSmallTest n cyc =
    it ("WinogradSmall(" ++ show n ++ "--" ++ show cyc ++ ")") $
    toMatrix (winogradSmall n w cyc) @?= toMatrix (F n w)
  where
    w :: Exp (Complex Double)
    w = omega n

winogradLargeTest :: (Int, Int) -> Spec
winogradLargeTest (r, s) = it ("WinogradLarge(" ++ show r ++ ", " ++ show s ++ ")") $
    sequence_ [toMatrix (winogradLarge r s w cr cc) @?= toMatrix (F n w)
                | cr <- take 2 $ getWinogradTriple' r s w getCycs
                , cc <- take 2 $ getWinogradTriple' s r w getCycs]
  where
    n :: Int
    n = r*s

    w :: Exp (Complex Double)
    w = omega n

winogradTwoTest :: Int -> Spec
winogradTwoTest n@2 = describe ("Winograd N=" ++ show n) $ do
    sequence_ [it ("DIF(2,1): w^" ++ show i) $ toMatrix (winogradDIF 2 1 (w^i)) @?= toMatrix (F n (w^i)) | i <- [1..n-1], rem i 2 /= 0]
    sequence_ [it ("DIT(2,1): w^" ++ show i) $ toMatrix (winogradDIT 2 1 (w^i)) @?= toMatrix (F n (w^i)) | i <- [1..n-1], rem i 2 /= 0]
  where
    w :: Exp (Complex Double)
    w = omega n

winogradTwoTest n = describe ("Winograd N=" ++ show n) $ do
    sequence_ [it ("DIF(" ++ show r ++ "," ++ show s ++ "): w^" ++ show i) $ toMatrix (winogradDIF r s (w^i)) @?= toMatrix (F n (w^i)) | i <- [1,n-1], rem i 2 /= 0, (r,s) <- factors n]
    sequence_ [it ("DIT(" ++ show r ++ "," ++ show s ++ "): w^" ++ show i) $ toMatrix (winogradDIT r s (w^i)) @?= toMatrix (F n (w^i)) | i <- [1,n-1], rem i 2 /= 0, (r,s) <- factors n]
    sequence_ [it ("SplitRadix: w^" ++ show i) $ toMatrix (winogradSplitRadix n (w^i)) @?= toMatrix (F n (w^i)) | i <- [1,n-1], rem i 2 /= 0]
    sequence_ [it ("ConjPairSplitRadix: w^" ++ show i) $ toMatrix (winogradConjPairSplitRadix n (w^i)) @?= toMatrix (F n (w^i)) | i <- [1,n-1], rem i 2 /= 0]
    sequence_ [it ("SplitRadix8: w^" ++ show i) $ toMatrix (winogradSplitRadix8 n (w^i)) @?= toMatrix (F n (w^i)) | i <- [1,n-1], n >= 8, rem i 2 /= 0]
    sequence_ [it ("ImprovedSplitRadix: w^" ++ show i) $ toMatrix (winogradImprovedSplitRadix n (w^i)) @?= toMatrix (F n (w^i)) | i <- [1,n-1], rem i 2 /= 0]
  where
    w :: Exp (Complex Double)
    w = omega n

winogradPowerTest :: Int -> Int -> Spec
winogradPowerTest p k = describe ("WinogradPower(" ++ show p ++ "^" ++ show k ++ ")") $
    sequence_
    [it ("w^" ++ show i ++ " -- p': " ++ show cP ++ " -- s: " ++ show cS) $
        toMatrix (winogradPower p k (w^^i) cP cS) @?= toMatrix (F n (w^^i))
        | i  <- [1,8],
          cP <- take 2 $ getCycs p',
          cS <- take 2 $ getCycs s]
  where
    n, p', s :: Int
    n  = p^k
    p' = (p-1)
    s  = p * p'

    w :: Exp (Complex Double)
    w = omega n

winogradSquareTest :: Int -> Spec
winogradSquareTest p = describe ("WinogradSquare(" ++ show p ++ "^2)") $
    sequence_
    [it ("w^" ++ show i) $ toMatrix f @?= toMatrix (F n (w^^i))
        | i  <- [1,p],
          f  <- take 1 $ winogradSquare p k (w^^i) getCycs]
  where
    n, p', s, k :: Int
    k  = 2
    n  = p^k
    p' = (p-1)
    s  = p * p'

    w :: Exp (Complex Double)
    w = omega n

-- DIT test
ditTest :: forall p . RootOfUnity (Exp p)
        => Proxy p
        -> Int
        -> Spec
ditTest _ n = it ("DIT(2^" ++ show n ++ ")") $
    toMatrix fft_spl @?= toMatrix (DFT (2^n))
  where
    fft_spl :: SPL (Exp p)
    fft_spl = dit (2^n)

-- DIF test
difTest :: forall p . RootOfUnity (Exp p)
        => Proxy p
        -> Int
        -> Spec
difTest _ n = it ("DIF(2^" ++ show n ++ ")") $
    toMatrix fft_spl @?= toMatrix (DFT (2^n))
  where
    fft_spl :: SPL (Exp p)
    fft_spl = dif (2^n)

-- Split radix test
splitRadixTest :: forall p . RootOfUnity (Exp p)
               => Proxy p
               -> Int
               -> Spec
splitRadixTest _ n = it ("SplitRadix(4^" ++ show n ++ ")") $
    toMatrix (splitRadix (4^n) w) @?= toMatrix (DFT (4^n))
  where
    w :: Exp p
    w = omega (4^n)

-- Split radix 8 test
splitRadix8Test :: forall p . RootOfUnity (Exp p)
                => Proxy p
                -> Int
                -> Spec
splitRadix8Test _ n = it ("SplitRadix(8^" ++ show n ++ ")") $
    toMatrix (splitRadix8 (8^n) w) @?= toMatrix (DFT (8^n))
  where
    w :: Exp p
    w = omega (8^n)

-- Split radix test
conjPairSplitRadixTest :: forall p . RootOfUnity (Exp p)
                       => Proxy p
                       -> Int
                       -> Spec
conjPairSplitRadixTest _ n = it ("ConjPairSplitRadix(4^" ++ show n ++ ")") $
    toMatrix (conjPairSplitRadix (4^n) w) @?= toMatrix (DFT (4^n))
  where
    w :: Exp p
    w = omega (4^n)

-- Split radix test
improvedSplitRadixTest :: forall p . (Floating (Exp p), RootOfUnity (Exp p))
                       => Proxy p
                       -> Int
                       -> Spec
improvedSplitRadixTest _ n = it ("ImprovedSplitRadix(4^" ++ show n ++ ")") $
    toMatrix (improvedSplitRadix (4^n) w) @?= toMatrix (DFT (4^n))
  where
    w :: Exp p
    w = omega (4^n)

-- Lowest-opcount test
opcountSearchTest :: Int -> Spec
opcountSearchTest n = it ("OpcountSearch(" ++ show n ++ ")") $ do
    Re e <- runSpiralWith mempty $
            searchOpCount (Re (DFT n) :: SPL (Exp Double))
    toMatrix e @?= toMatrix (DFT n)

max_stars_repo_name: mainland/hspiral

max_issues_repo_name: mainland/hspiral

max_forks_repo_name: mainland/hspiral


text:
-- A simple example testing arithmetic functions.
{-# LANGUAGE FlexibleInstances #-}
import           Control.Monad
import           Data.Complex
import           Data.Either
import           Data.Matrix
import           Data.Ord
import           Data.Vector                (Vector)
import           Numeric.Natural
import           QuickSpec                  hiding (A)
import           QuickSpec.Internal         hiding (A)
import qualified QuickSpec.Internal.Haskell as Haskell
import           Test.QuickCheck
import           Test.QuickCheck.Instances
import           Twee.Pretty

dim = 5

type A = Rational
type T = Matrix A

instance Ord a => Ord (Matrix a) where
  compare = comparing toLists

instance (RealFrac a, Ord a) => Ord (Complex a) where
  compare = comparing (\(x :+ y) -> (truncate (x*1000), truncate (y*1000)))

--instance Ord Cyclotomic where
--  compare = comparing show
--
--instance Arbitrary Cyclotomic where
--  arbitrary = sized arb
--    where
--      arb n =
--        frequency $
--          [(5, fromInteger  <$> arbitrary),
--           (5, fromRational <$> arbRat),
--           (5, gaussianRat <$> arbRat <*> arbRat),
--           (5, polarRat <$> arbRat <*> arbRat),
--           (5, e <$> fmap (succ . abs) arbitrary),
--           (3, sqrtInteger <$> arbitrary),
--           (3, sqrtRat <$> arbRat),
--           (2, sinRev <$> arbRat),
--           (2, cosRev <$> arbRat),
--           (2, return i)] ++
--          concat
--            [ [(5, conj <$> arb (n-1)),
--               (3, real <$> arb (n-1)),
--               (3, imag <$> arb (n-1))]
--            | n > 0 ] ++
--          concat
--            [ let arb2 = arb (n `div` 2) in
--              [(n*5, (+) <$> arb2 <*> arb2),
--               (n*5, (*) <$> arb2 <*> arb2)]
--            | n > 0 ]
--      arbRat = sized $ \n -> do
--        x <- choose (-(n `min` 3), n `min` 3)
--        d <- fmap (succ . abs) (choose (0, n `min` 5))
--        n <- choose (0, d-1)
--        return (fromIntegral x + fromIntegral n / fromIntegral d)

instance Arbitrary T where
  arbitrary =
    makeHermitian <$> frequency [
--(1, return $ zero dim dim),
--      (1, return $ identity dim),
--      (1, return $ negate $ identity dim),
--      (1, genDiagonal),
      (5, fromList dim dim <$> infiniteListOf genA)]

genA :: Gen A
genA = fromInteger <$> arbitrary

genDiagonal :: Gen T
genDiagonal = diagonalList dim 0 <$> infiniteListOf genA

isDiagonal :: T -> Bool
isDiagonal m = m == diagonal 0 (getDiag m)

makeUpperTriangular :: T -> T
makeUpperTriangular = mapPos (\(i, j) x -> if i > j then 0 else x)

isUpperTriangular :: T -> Bool
isUpperTriangular m = makeUpperTriangular m == m

genUpperTriangular :: Gen T
genUpperTriangular = makeUpperTriangular <$> arbitrary

makeLowerTriangular :: T -> T
makeLowerTriangular = mapPos (\(i, j) x -> if i > j then 0 else x)

isLowerTriangular :: T -> Bool
isLowerTriangular m = makeLowerTriangular m == m

genLowerTriangular :: Gen T
genLowerTriangular = makeLowerTriangular <$> arbitrary

makeHermitian :: T -> T
makeHermitian m = m + {-fmap conjugate-} transpose m

isHermitian :: T -> Bool
isHermitian m = m == {-fmap conjugate-} transpose m

genHermitian :: Gen T
genHermitian = makeHermitian <$> arbitrary

instance Floating Rational where
  sqrt = error "sqrt"
instance RealFloat Rational

newtype RowV = RowV {unRowV :: T} deriving (Eq, Ord)
newtype ColV = ColV {unColV :: T} deriving (Eq, Ord)

instance Arbitrary RowV where
  arbitrary = RowV <$> oneof [return $ zero 1 dim, fromList 1 dim <$> infiniteList]
instance Arbitrary ColV where
  arbitrary = (ColV <$> transpose) . unRowV <$> arbitrary

times :: T -> T -> T
times x y = scaleMatrix (1/2) (x*y+y*x)

pow :: Natural -> T -> T
pow 1 t = t
pow n t | n > 0 = times t (pow (n-1) t)

main = quickSpec [
  series [
    background
    [withMaxTermSize 7,
     -- Good pruning settings for Hard Maths (TM)
     --withPruningTermSize 9,
     arith (Proxy :: Proxy A),
     arith (Proxy :: Proxy Natural) `without` ["0"],
     con "*" ((*) :: Natural -> Natural -> Natural),
     monoType (Proxy :: Proxy T),
     con "0" (zero dim dim :: T),
     con "I" (identity dim :: T),
     con "+" ((+) :: T -> T -> T),
     vars ["A", "B", "C"] (Proxy :: Proxy T),
     vars ["n", "m", "l"] (Proxy :: Proxy Natural),
     instFun (arbitrary `suchThat` (\n -> n > 0 && n < 4) :: Gen Natural)],
     --con "*" ((*) :: T -> T -> T),
     --con "-" (negate :: T -> T)],
    --[(con "&&" (&&))],
    --[predicateGen "hermitian" isHermitian (\() -> liftM2 (,) genHermitian (return ()))],
    con "∘" times,
    con "^" (flip pow)]]
{-    [--con "perm" (permMatrix dim :: Int -> Int -> T),
     con "transpose" (transpose :: T -> T),
     --con "inverse" (either (const Nothing) Just . (inverse :: T -> Either String T)),
     con "trace" (trace :: T -> A),
     con "diagProd" (diagProd :: T -> A),
     con "det" (detLU :: T -> A)],
     {-predicateGen "isDiagonal" isDiagonal (\() -> liftM2 (,) genDiagonal (return ())),
     predicateGen "isUpperTriangular" isUpperTriangular (\() -> liftM2 (,) genUpperTriangular (return ())),
     predicateGen "isLowerTriangular" isLowerTriangular (\() -> liftM2 (,) genLowerTriangular (return ()))],-}
    [monoType (Proxy :: Proxy RowV),
     con "0R" (RowV (zero 1 dim)),
     con "+R" (\(RowV x) (RowV y) -> RowV (x+y)),
     con "*R" (\(RowV x) y -> RowV (x*y)),
     con "-R" (\(RowV x) -> RowV (negate x))],
    [con "transposeR" (\(RowV x) -> ColV (transpose x)),
     con "detR" (detLU . unRowV)],
    [monoType (Proxy :: Proxy ColV),
     con "0C" (ColV (zero dim 1)),
     con "+C" (\(ColV x) (ColV y) -> ColV (x+y)),
     con "*C" (\(ColV x) y -> ColV (x*y)),
     con "-C" (\(ColV x) -> ColV (negate x))],
    [con "transposeC" (\(ColV x) -> RowV (transpose x)),
     con "detC" (detLU . unColV)]]]-}

max_stars_repo_name: danwdart/quickspec

max_issues_repo_name: danwdart/quickspec

max_forks_repo_name: danwdart/quickspec


text:
{-|
Module      : Data.Meta.Multiarm
Description : Matrix function for graphs
Copyright   : (c) Thodoris Papakonstantinou, 2019
License     : GPL-3
Maintainer  : mail@tpapak.com
Stability   : experimental
Portability : POSIX

-}

module Data.Meta.Matrices
  ( createBMatrix
  , createAMatrix
  , createCMatrix
  )
where

import           Control.Applicative
import           Data.List
import qualified Data.Map.Strict               as Map
import qualified Data.Set                      as Set
import           Data.Maybe
import           Data.Either
import           Data.Meta.Effects
import qualified Numeric.LinearAlgebra         as LA
import qualified Numeric.LinearAlgebra.Devel   as LAD
import qualified Data.Graph.AdjacencyList      as G
import qualified Data.Graph.AdjacencyList.BFS      as BFS

-- | create the edge vertex adjacency matrix given a graph
createBMatrix :: G.Graph -> LA.Matrix Double
createBMatrix g =
  let vs  = G.vertices g
      nvs = length vs
      es  = G.edges g
      nes = length es
      rowfromEdge e =
          let
            (u, v) = G.toTuple e
            r      = LA.fromList $ replicate nvs 0.0 :: LA.Vector Double
            r' =
              LAD.mapVectorWithIndex
                (\i _ ->
                  if i == u - 1 then 1.0 else if i == v - 1 then (-1.0) else 0.0
                )
                r :: LA.Vector Double
          in
            r'
      rows = map rowfromEdge es
      b    = LA.fromRows rows
  in  b

createAMatrix :: G.Graph -> LA.Matrix Double
createAMatrix g =
  let vs  = G.vertices g
      nvs = length vs
      es  = G.edges g
      nes = length es
      rowfromEdge e =
          let
            (u, v) = G.toTuple e
            r      = LA.fromList $ replicate (nvs - 1) 0.0 :: LA.Vector Double
            r' =
              LAD.mapVectorWithIndex
                (\i _ ->
                  if i == u - 2 then 1.0 else if i == v - 2 then (-1.0) else 0.0
                )
                r :: LA.Vector Double
           in r'
      rows = map rowfromEdge es
      a    = LA.fromRows rows
  in  a

-- | create the Vertex Edge adjacency matrix given a graph

-- | create the Vertex Edge adjacency matrix given a graph
createCMatrix :: G.Graph -> LA.Matrix Double
createCMatrix g =
  let vs  = G.vertices g
      nvs = length vs
      es  = G.edges g
      nes = length es
      rowfromVertex v =
        let r = map (\e ->
              let (u, v') = G.toTuple e
                  fillcell x | x==u = (-1.0)
                  fillcell x | x==v' = (1.0)
                  fillcell _ = (0.0)
               in  fillcell v
                  ) es
         in LA.fromList $ r :: LA.Vector Double
      rows = map rowfromVertex vs
      c    = LA.fromRows rows
  in  c

-- | create the Vertex Edge adjacency matrix given a graph
createXMatrix :: G.Graph -> LA.Matrix Double
createXMatrix g =
  let vs  = G.vertices g
      nvs = length vs
      es  = G.edges g
      nes = length es
      firstVertex = head vs
      bfsNet = BFS.bfs g firstVertex
      spanningtree = BFS.spanningTree bfsNet
      rowfromVertex v =
        let r = map (\e ->
              let (u, v') = G.toTuple e
                  fillcell x | x==u = (-1.0)
                  fillcell x | x==v' = (1.0)
                  fillcell _ = (0.0)
               in  fillcell v
                  ) es
         in LA.fromList $ r :: LA.Vector Double
      rows = map rowfromVertex vs
      c    = LA.fromRows rows
  in  c

max_stars_repo_name: tpapak/meta-analysis

max_issues_repo_name: tpapak/meta-analysis

max_forks_repo_name: tpapak/meta-analysis


text:
module Resultant where

import Numeric.LinearAlgebra

type Poly = [R]
  -- 多項式はこれで表す
  -- 主係数から定数までを，次数の大きい項から順に全て列挙
  -- 例 : x^2 - 2 -> [1, 0, (-2)]

deg :: Poly -> Int -- 次数計算
deg f = length f - 1

mj :: Poly -> Poly -> Int -> Matrix R -- M_j (f, g) の計算
mj f g j
  | j < 0 || k <= 0 = undefined
  | otherwise       = (k >< k) $ concat -- 行ごとの配列になっているものをくっつける
                      $  ((makeMsRow f k) <$> [0..(n - j - 1)]) -- 上側
                      ++ ((makeMsRow g k) <$> [0..(m - j - 1)]) -- 下側
  where
    (m, n) = (deg f, deg g)
    k = m + n - 2 * j

syl :: Poly -> Poly -> Matrix R
syl f g = mj f g 0 -- シルベスター行列

psc :: Poly -> Poly -> Int -> R -- PSC_j (f, g) の計算
psc f g = det . mj f g

res :: Poly -> Poly -> R
res f g = det $ syl f g -- シルベスター行列の PSC

degGcd :: Poly -> Poly -> Int
degGcd f g = head $ filter (\j -> psc f g j /= 0) [0..] -- PSC が最初に 0 じゃなくなる j

makeMsRow :: Poly -> Int -> Int -> Poly          -- M_j の各行を生成
makeMsRow f k = flip adjustLength k . add0Head f -- 右に n ずらしてから長さ k の配列を作る

add0Head :: Poly -> Int -> Poly
add0Head a 0 = a
add0Head a n    -- n 個の 0 を先頭に加える
  | n < 0     = undefined
  | otherwise = 0 : add0Head a (n - 1)

adjustLength :: Poly -> Int -> Poly
adjustLength a n -- 配列の長さを n にする
  | dif == 0  = a
  | dif <  0  = take n a             -- 長い分を削る
  | otherwise = a ++ replicate dif 0 -- 短い分 0 を後ろに加える
  where
    dif = n - length a

mat :: Matrix R
mat = (6 >< 6) [ 3, 1, 2, 5, 2, 2
               , 5, 2, 3, 3, 1, 2
               , 2, 0, 1, 5, 3, 4
               , 4, 3, 5, 2, 1, 2
               , 1, 2, 3, 5, 6, 4
               , 6, 4, 1, 4, 3, 3
               ]


text:
-- |
-- Copyright: (c) 2020 Tom Westerhout
-- SPDX-License-Identifier: BSD-3-Clause
-- Maintainer: Tom Westerhout <14264576+twesterhout@users.noreply.github.com>
--
-- High-level wrapper around [PRIMME C
-- library](https://github.com/primme/primme). /Quote from README/:
--
-- @
--   PRIMME, pronounced as prime, is a high-performance library for computing a few
--   eigenvalues/eigenvectors, and singular values/vectors. PRIMME is especially
--   optimized for large, difficult problems. Real symmetric and complex Hermitian
--   problems, standard @A x = λ x@ and generalized @A x = λ B x@, are supported.
--   Besides, standard eigenvalue problems with a normal matrix are supported. It
--   can find largest, smallest, or interior singular/eigenvalues, and can use
--   preconditioning to accelerate convergence.
-- @
module Numeric.PRIMME
  ( -- * Defining the matrix

    -- | One of the great things about PRIMME library is that it works with
    -- block matrix-vector products (i.e. matrix-matrix products). Following the
    -- example of [vector](https://hackage.haskell.org/package/vector) library,
    -- we define two types for mutable and immutable dense blocks.
    MBlock (..),
    Block (..),
    -- | Now we are can define the "operator" we want to diagonalize. Since this
    -- library is meant to be used with rather large matrices, it might very
    -- well be the case that the matrix does not fit into the memory of your
    -- computer (even in sparse format such as CSR). Sometimes, however, we can
    -- define the "operator" implicitly, namely, by defining its action on a
    -- vector (or in case of PRIMME on a block of vectors). This is done with
    -- 'PrimmeOperator' type.
    PrimmeOperator,

    -- * Choosing what to compute

    -- | Having defined a 'PrimmeOperator' we now need to tell PRIMME what to
    -- compute. This is done by constructing a 'PrimmeOptions' object.
    PrimmeOptions (..),
    primmeDefaults,
    finalizeOptions,
    PrimmeTarget (..),

    -- * Diagonalizing
    eigh,
    eigh',

    -- * Logging
    PrimmeMonitor (..),
    PrimmeInfo (..),
    PrimmeEventInfo (..),
    PrimmeStats (..),
    primmePrettyInfo,

    -- * Dense Matrices

    -- | The most trivial example of an "operator" is of course a square dense
    -- matrix. We thus provide a function which constructs an operator from a
    -- matrix.
    primmeFromDense,

    -- * Misc
    getPrimmeVersion,
    PrimmeException (..),
    BlasDatatype (blasTag),
    BlasRealPart,
    BlasDatatypeTag (..),
  )
where

import Control.Exception.Safe (throw)
import Control.Monad (forM_, unless)
import Control.Monad.ST (RealWorld)
import Data.Complex
import Data.Proxy
import Data.Vector.Storable (Vector)
import qualified Data.Vector.Storable as V
import Data.Vector.Storable.Mutable (MVector)
import qualified Data.Vector.Storable.Mutable as MV
import Foreign.C.Types (CInt (..))
import Foreign.Ptr (Ptr)
import Foreign.Storable
import qualified Language.C.Inline as C
import qualified Language.C.Inline.Unsafe as CU
import Numeric.PRIMME.Context
import Numeric.PRIMME.Dense
import Numeric.PRIMME.Monitor
import Numeric.PRIMME.Options
import Numeric.PRIMME.Types
import Prelude hiding (init)

C.context (C.baseCtx <> primmeCtx)
C.include "<primme.h>"

-- | Diagonalize the operator to find the first few eigenpairs.
eigh ::
  forall a.
  BlasDatatype a =>
  PrimmeOptions ->
  PrimmeOperator a ->
  IO (Vector (BlasRealPart a), Block a, Vector (BlasRealPart a))
eigh options matrix = do
  let dim = pDim options
      numEvals = pNumEvals options
  (evals :: MVector RealWorld (BlasRealPart a)) <- MV.new numEvals
  (evecs :: MVector RealWorld a) <- MV.new (dim * numEvals)
  (rnorms :: MVector RealWorld (BlasRealPart a)) <- MV.new numEvals
  status <-
    MV.unsafeWith evals $ \evalsPtr ->
      MV.unsafeWith evecs $ \evecsPtr ->
        MV.unsafeWith rnorms $ \rnormsPtr ->
          withOptions options matrix $ \optionsPtr ->
            cPrimme evalsPtr evecsPtr rnormsPtr optionsPtr
  unless (status == 0) $ throwPrimmeError status
  evals' <- V.unsafeFreeze evals
  evecs' <- Block (dim, numEvals) dim <$> V.unsafeFreeze evecs
  rnorms' <- V.unsafeFreeze rnorms
  return (evals', evecs', rnorms')

mkEmptyBlock :: Storable a => Int -> Int -> IO (MBlock RealWorld a)
mkEmptyBlock rows columns
  | rows <= 0 || columns <= 0 = error $ "invalid shape specified: " <> show (rows, columns)
  | otherwise = MBlock (rows, columns) rows <$> MV.new (rows * columns)

copyBlock :: Storable a => MBlock RealWorld a -> Block a -> IO ()
copyBlock dest@(MBlock (r, c) _ _) src@(Block (r', c') _ _)
  | r /= r' || c /= c' = error $ "block shape mismatch: " <> show (r, c) <> " != " <> show (r', c')
  | otherwise = do
    src' <- unsafeThaw src
    forM_ [0 .. (c - 1)] $ \i ->
      MV.copy (getColumn i dest) (getColumn i src')
  where
    getColumn :: Storable a => Int -> MBlock RealWorld a -> MVector RealWorld a
    getColumn i (MBlock (rows, _) stride v) = MV.slice (i * stride) rows v

unsafeThaw :: Storable a => Block a -> IO (MBlock RealWorld a)
unsafeThaw (Block shape stride v) = MBlock shape stride <$> V.unsafeThaw v

unsafeFreeze :: Storable a => MBlock RealWorld a -> IO (Block a)
unsafeFreeze (MBlock shape stride v) = Block shape stride <$> V.unsafeFreeze v

sliceBlock1 :: Storable a => Int -> Int -> MBlock s a -> MBlock s a
sliceBlock1 i n (MBlock (r, c) stride v)
  | i < 0 || n < 0 || i + n > c = error $ "invalid slice: [" <> show i <> ", " <> show (i + n) <> ")"
  | otherwise = MBlock (r, n) stride $ MV.slice (i * stride) (n * stride) v

eigh' ::
  forall a.
  BlasDatatype a =>
  PrimmeOptions ->
  MBlock RealWorld a ->
  PrimmeOperator a ->
  IO (Vector (BlasRealPart a), Block a, Vector (BlasRealPart a))
eigh' options init@(MBlock (dim', initSize) _ _) matrix
  | pDim options /= dim' =
    error $
      "'init' has wrong shape: "
        <> show (dim', initSize)
        <> "; expected a block with "
        <> show (pDim options)
        <> " rows"
  | otherwise = do
    let dim = pDim options
        numEvals = pNumEvals options
    evecs@(MBlock (_, _) evecsStride evecsData) <- case initSize < numEvals of
      True -> do
        temp <- mkEmptyBlock dim numEvals
        copyBlock (sliceBlock1 0 initSize temp) =<< unsafeFreeze init
        return temp
      False -> return init
    (evals :: MVector RealWorld (BlasRealPart a)) <- MV.new numEvals
    (rnorms :: MVector RealWorld (BlasRealPart a)) <- MV.new numEvals
    status <-
      MV.unsafeWith evals $ \evalsPtr ->
        MV.unsafeWith evecsData $ \evecsPtr ->
          MV.unsafeWith rnorms $ \rnormsPtr ->
            withOptions options matrix $ \optionsPtr -> do
              let c_initSize = fromIntegral initSize
                  c_evecsStride = fromIntegral evecsStride
              [C.block| void {
                $(primme_params* optionsPtr)->initSize = $(int c_initSize);
                $(primme_params* optionsPtr)->ldevecs = $(PRIMME_INT c_evecsStride);
              } |]
              cPrimme evalsPtr evecsPtr rnormsPtr optionsPtr
    unless (status == 0) $ throwPrimmeError status
    evals' <- V.unsafeFreeze evals
    evecs' <- unsafeFreeze $ sliceBlock1 0 numEvals evecs
    rnorms' <- V.unsafeFreeze rnorms
    return (evals', evecs', rnorms')

throwPrimmeError :: CInt -> IO a
throwPrimmeError c = case (- c) of
  0 -> error "no error"
  3 -> throw PrimmeMaximumIterationsReached
  40 -> throw PrimmeLapackFailure
  41 -> throw PrimmeUserFailure
  _ -> throw $ PrimmeOtherFailure c

foreign import ccall "sprimme"
  sprimme :: Ptr Float -> Ptr Float -> Ptr Float -> Ptr Cprimme_params -> IO CInt

foreign import ccall "dprimme"
  dprimme :: Ptr Double -> Ptr Double -> Ptr Double -> Ptr Cprimme_params -> IO CInt

foreign import ccall "cprimme"
  cprimme :: Ptr Float -> Ptr (Complex Float) -> Ptr Float -> Ptr Cprimme_params -> IO CInt

foreign import ccall "zprimme"
  zprimme :: Ptr Double -> Ptr (Complex Double) -> Ptr Double -> Ptr Cprimme_params -> IO CInt

cPrimme ::
  forall a.
  (BlasDatatype a) =>
  Ptr (BlasRealPart a) ->
  Ptr a ->
  Ptr (BlasRealPart a) ->
  Ptr Cprimme_params ->
  IO CInt
cPrimme = case blasTag (Proxy :: Proxy a) of
  FloatTag -> sprimme
  DoubleTag -> dprimme
  ComplexFloatTag -> cprimme
  ComplexDoubleTag -> zprimme

getPrimmeVersion :: (Int, Int)
getPrimmeVersion =
  ( fromIntegral [CU.pure| int { PRIMME_VERSION_MAJOR } |],
    fromIntegral [CU.pure| int { PRIMME_VERSION_MINOR } |]
  )

max_stars_repo_name: twesterhout/primme-hs

max_issues_repo_name: twesterhout/primme-hs

max_forks_repo_name: twesterhout/primme-hs


text:
-- {-# LANGUAGE DeriveDataTypeable #-}

module Main where

-- import Lib

import Control.Monad.State
import Data.List (transpose)
import Data.Vector (snoc)
import Numeric.LinearAlgebra
  ( Linear (scale),
    Matrix,
    Vector,
    cmap,
    matrix,
    maxIndex,
    outer,
    vector,
    (#>),
  )
import Text.CSV (parseCSVFromFile)

type WeightValue = Matrix Double

type WeightState = Matrix Double

predict ::
  Matrix Double -> -- curried - Weights
  Vector Double -> -- curried - Input vector
  Vector Double -- Guess
predict w x = cmap sigmoid w #> x

predictEpoch ::
  Matrix Double -> -- w - Weights
  [Vector Double] -> -- x - An epoch of input vectorors
  [Vector Double] -- A guess for each input in the epoch
predictEpoch w x = map (cmap sigmoid) (fmap (w #>) x)

sigmoid :: Double -> Double
sigmoid x = 1.0 / (1.0 + exp (negate x))

stocTrain ::
  Vector Double -> -- x - Single input vectoror
  Vector Double -> -- l - label
  State WeightState () -- Weights
stocTrain x l = do
  w <- get
  let error = predict w x - l
      sigDeriv = error * cmap (1 -) error
      delta = scale 0.01 $ sigDeriv `outer` x
      newWeights = w - delta
  put newWeights

train ::
  Int -> -- t - number of epochs to train on
  [Vector Double] -> -- x - An epoch of inputs
  [Vector Double] -> -- l - An epoch of labels
  State WeightState WeightValue -- Monad stuff
train t xs ls
  | t == 0 = get
  | otherwise = do
    trainEpoch xs ls
    train (t - 1) xs ls

trainEpoch ::
  [Vector Double] -> -- x - An epoch of inputs
  [Vector Double] -> -- l - An epoch of labels
  State WeightState WeightValue -- Weights
trainEpoch [] [] = get
trainEpoch (x : xs) (l : ls) = do
  stocTrain x l
  trainEpoch xs ls

test ::
  Matrix Double -> -- w - Weights
  [Vector Double] -> -- x - An epoch of input vectors
  [Vector Double] -> -- l - An epoch of labels
  Double -- A guess for each input in the epoch
test w x l = intDiv correct total
  where
    correct = sum (zipWith checkSame l (predictEpoch w x))
    total = length l

intDiv :: Int -> Int -> Double
intDiv n l = fromIntegral n / fromIntegral l

checkSame :: Vector Double -> Vector Double -> Int
checkSame x y
  | maxIndex x == maxIndex y = 1
  | otherwise = 0

run ::
  [[String]] -> -- tesCSV - read csv of test data
  [[String]] -> --  trainCSV - read csv of train data
  Double -- correct / total
run testCSV trainCSV = test w xTest lTest
  where
    (xTest, lTest) = csvToInputsLabels testCSV
    w = evalState (train 1 xTrain lTrain) initW
    (xTrain, lTrain) = csvToInputsLabels trainCSV
    initW = matrix 785 (replicate (785 * 10) 0)

csvToInputsLabels ::
  [[String]] -> -- read csv
  ([Vector Double], [Vector Double]) -- (x,l)
csvToInputsLabels csv = (finInputs, oneHotLabels)
  where
    oneHotLabels = fmap oneHotVector rawLabels
    finInputs = fmap vector scaledBiasedInputs
    scaledBiasedInputs = fmap (++ [1.0]) (scaleInput rawInputs)
    (rawLabels, rawInputs) = labelInputSplit (map (map read) (init csv))

scaleInput :: [[Double]] -> [[Double]]
scaleInput = map (map (/ 255))

labelInputSplit :: [[Double]] -> ([Double], [[Double]])
labelInputSplit xs = (head t, transpose (tail t))
  where
    t = transpose xs

oneHotVector :: Double -> Vector Double
oneHotVector x = vector $ replicate int 0 ++ [1] ++ replicate (9 - int) 0
  where
    int = round x

main :: IO ()
main = do
  -- testCSV <- parseCSVFromFile "./data/mnist_test.csv"
  -- trainCSV <- parseCSVFromFile "./data/mnist_train.csv"
  -- case testCSV of
  --   Left err -> print err
  --   Right testCSV ->
  --     case trainCSV of
  --       Left err -> print err
  --       Right trainCSV -> print (run testCSV testCSV)

  let inputs = map vector [[0, 0, 1], [0, 1, 1], [1, 0, 1], [1, 1, 1]]
  let labels = map vector [[0, 1], [0, 1], [0, 1], [1, 0]]
  let weights = matrix 3 (replicate 6 0)
  print "Initialized Weights"
  print weights
  print "target"
  print labels
  print "untrained prediction"
  print $ predictEpoch weights inputs
  print $ test weights inputs labels
  let trainedWeights = evalState (train 10 inputs labels) weights
  print "trained weights"
  print trainedWeights
  print "target"
  print labels
  print "trained prediciton"
  print $ predictEpoch trainedWeights inputs
  print "Test"
  print $ test trainedWeights inputs labels


text:
module HaskellWorks.Playground.Linear where

import Numeric.LinearAlgebra
import Prelude               hiding ((<>))

variance :: Matrix Double -> Matrix Double
variance m = x <> ((m - asRow rdMean) ** 2)
  where (rdMean, _) = meanCov m
        r = rows m
        x = konst (1 / fromIntegral r) (1, r)

max_stars_repo_name: haskell-works/hw-playground-linear

max_issues_repo_name: haskell-works/hw-playground-linear

max_forks_repo_name: haskell-works/hw-playground-linear


text:
{-# LANGUAGE BangPatterns           #-}
{-# LANGUAGE DataKinds              #-}
{-# LANGUAGE FlexibleContexts       #-}
{-# LANGUAGE GADTs                  #-}
{-# LANGUAGE LambdaCase             #-}
{-# LANGUAGE ScopedTypeVariables    #-}
{-# LANGUAGE TupleSections          #-}
{-# LANGUAGE TypeApplications       #-}
{-# LANGUAGE TypeFamilyDependencies #-}

module Backprop.Learn.Model (
    module M, Backprop(..)
  , runLearn_, runLearnStoch_, runLearnStateless_, runLearnStochStateless_
  , gradLearn, gradLearnStoch
  -- * Work with parameters
  , initParam, initParamMaybe, initParamNormal, initParamNormalMaybe
  , encodeParam, decodeParam, decodeParamOrFail, saveParam, loadParam, loadParamOrFail
  -- * Iterated runners
  , iterateLearn, iterateLearnM, iterateLearnStoch
  , scanLearn, scanLearnStoch
  -- * No final state
  , iterateLearn_, iterateLearnM_, iterateLearnStoch_
  , scanLearn_, scanLearnStoch_
  -- * "Prime" runners
  , primeLearn, primeLearnStoch, selfPrime, selfPrimeM
  ) where

import           Backprop.Learn.Initialize
import           Backprop.Learn.Model.Class       as M
import           Backprop.Learn.Model.Combinator  as M
import           Backprop.Learn.Model.Function    as M
import           Backprop.Learn.Model.Neural      as M
import           Backprop.Learn.Model.Neural.LSTM as M
import           Backprop.Learn.Model.Parameter   as M
import           Backprop.Learn.Model.Regression  as M
import           Backprop.Learn.Model.State       as M
import           Backprop.Learn.Model.Stochastic  as M
import           Control.Monad.Primitive
import           Control.Monad.ST
import           Control.Monad.Trans.State
import           Data.Bifunctor
import           Data.Foldable
import           Data.Functor.Identity
import           Data.Type.Mayb
import           Data.Word
import           Numeric.Backprop
import           Statistics.Distribution
import qualified Data.Binary                      as Bi
import qualified Data.ByteString.Lazy             as BSL
import qualified Data.Vector.Unboxed              as VU
import qualified System.Random.MWC                as MWC

-- TODO: this can be more efficient by breaking out into separate functions
runLearn_
    :: (Learn a b l, MaybeC Backprop (LStateMaybe l), Backprop b)
    => l
    -> LParam_ I l
    -> a
    -> LState_ I l
    -> (b, LState_ I l)
runLearn_ l mp x ms = case mp of
    N_ -> case ms of
      N_       -> (evalBP (runLearnStateless l N_) x, N_)
      J_ (I s) -> second (J_ . I)
                . evalBP2 (\x' s' -> uncurry T2
                                       . second fromJ_
                                       $ runLearn l N_ x' (J_ s')
                          ) x
                $ s
    J_ (I p) -> case ms of
      N_       -> (evalBP2 (runLearnStateless l . J_) p x, N_)
      J_ (I s) -> second (J_ . I)
                . evalBPN (\(p' :< x' :< s' :< Ø) ->
                                uncurry T2
                              . second fromJ_
                              $ runLearn l (J_ p') x' (J_ s')
                          )
                $ (p ::< x ::< s ::< Ø)

runLearnStoch_
    :: (Learn a b l, MaybeC Backprop (LStateMaybe l), Backprop b, PrimMonad m)
    => l
    -> MWC.Gen (PrimState m)
    -> LParam_ I l
    -> a
    -> LState_ I l
    -> m (b, LState_ I l)
runLearnStoch_ l g mp x ms = do
    -- TODO: is this the best way to handle this?
    seed <- MWC.uniformVector @_ @Word32 @VU.Vector g 2
    pure $ case mp of
      N_ -> case ms of
        N_       -> (, N_) $ evalBP (\x' -> runST $ do
            g' <- MWC.initialize seed
            runLearnStochStateless l g' N_ x'
          ) x
        J_ (I s) -> second (J_ . I) $ evalBP2 (\x' s' -> runST $ do
            g' <- MWC.initialize seed
            uncurry T2 . second fromJ_
              <$> runLearnStoch l g' N_ x' (J_ s')
          ) x s
      J_ (I p) -> case ms of
        N_       -> (, N_) $ evalBP2 (\p' x' -> runST $ do
            g' <- MWC.initialize seed
            runLearnStochStateless l g' (J_ p') x'
          ) p x
        J_ (I s) -> second (J_ . I) $ evalBPN (\(p' :< x' :< s' :< Ø) -> runST $ do
            g' <- MWC.initialize seed
            uncurry T2 . second fromJ_
              <$> runLearnStoch l g' (J_ p') x' (J_ s')
          ) (p ::< x ::< s ::< Ø)

runLearnStateless_
    :: (Learn a b l, NoState l)
    => l
    -> LParam_ I l
    -> a
    -> b
runLearnStateless_ l = \case
    N_       -> evalBP  (runLearnStateless l N_  )
    J_ (I p) -> evalBP2 (runLearnStateless l . J_) p

runLearnStochStateless_
    :: (Learn a b l, NoState l, PrimMonad m)
    => l
    -> MWC.Gen (PrimState m)
    -> LParam_ I l
    -> a
    -> m b
runLearnStochStateless_ l g mp x = do
    seed <- MWC.uniformVector @_ @Word32 @VU.Vector g 2
    pure $ case mp of
      N_       -> evalBP  (\x' -> runST $ do
          g' <- MWC.initialize seed
          runLearnStochStateless l g' N_ x'
        ) x
      J_ (I p) -> evalBP2 (\p' x' -> runST $ do
          g' <- MWC.initialize seed
          runLearnStochStateless l g' (J_ p') x'
        ) p x

gradLearn
    :: (Learn a b l, NoState l, Backprop a, Backprop b, MaybeC Backprop (LParamMaybe l))
    => l
    -> LParam_ I l
    -> a
    -> (LParam_ I l, a)
gradLearn l = \case
    N_       ->       (N_,)    . gradBP  (runLearnStateless l   N_)
    J_ (I p) -> first (J_ . I) . gradBP2 (runLearnStateless l . J_) p

gradLearnStoch
    :: (Learn a b l, NoState l, Backprop a, Backprop b, MaybeC Backprop (LParamMaybe l), PrimMonad m)
    => l
    -> MWC.Gen (PrimState m)
    -> LParam_ I l
    -> a
    -> m (LParam_ I l, a)
gradLearnStoch l g mp x = do
    seed <- MWC.uniformVector @_ @Word32 @VU.Vector g 2
    pure $ case mp of
      N_       ->       (N_,)    $ gradBP  (\x' -> runST $ do
          g' <- MWC.initialize seed
          runLearnStochStateless l g' N_ x'
        ) x
      J_ (I p) -> first (J_ . I) $ gradBP2 (\p' x' -> runST $ do
          g' <- MWC.initialize seed
          runLearnStochStateless l g' (J_ p') x'
        ) p x

iterateLearn
    :: (Learn a b l, Backprop b, MaybeC Backprop (LStateMaybe l))
    => (b -> a)         -- ^ loop
    -> Int              -- ^ num times
    -> l
    -> LParam_ I l
    -> a
    -> LState_ I l
    -> ([b], LState_ I l)
iterateLearn f n l p x = runIdentity . iterateLearnM (Identity . f) n l p x

iterateLearn_
    :: (Learn a b l, Backprop b, MaybeC Backprop (LStateMaybe l))
    => (b -> a)         -- ^ loop
    -> l
    -> LParam_ I l
    -> a
    -> LState_ I l
    -> [b]
iterateLearn_ f l p = go
  where
    go !x !s = y : go (f y) s'
      where
        (y, s') = runLearn_ l p x s

selfPrime
    :: (Learn a b l, Backprop b, MaybeC Backprop (LStateMaybe l))
    => (b -> a)         -- ^ loop
    -> l
    -> LParam_ I l
    -> a
    -> LState_ I l
    -> [LState_ I l]
selfPrime f l p = go
  where
    go !x !s = s' : go (f y) s'
      where
        (y, s') = runLearn_ l p x s

iterateLearnM
    :: (Learn a b l, Backprop b, MaybeC Backprop (LStateMaybe l), Monad m)
    => (b -> m a)           -- ^ loop
    -> Int                  -- ^ num times
    -> l
    -> LParam_ I l
    -> a
    -> LState_ I l
    -> m ([b], LState_ I l)
iterateLearnM f n l p = go 0
  where
    go !i !x !s
      | i <= n    = do
          let (y, s') = runLearn_ l p x s
          (ys, s'') <- flip (go (i + 1)) s' =<< f y
          pure (y : ys, s'')
      | otherwise = pure ([], s)

iterateLearnM_
    :: (Learn a b l, Backprop b, MaybeC Backprop (LStateMaybe l), Monad m)
    => (b -> m a)           -- ^ loop
    -> Int                  -- ^ num times
    -> l
    -> LParam_ I l
    -> a
    -> LState_ I l
    -> m [b]
iterateLearnM_ f n l p x = fmap fst . iterateLearnM f n l p x

selfPrimeM
    :: (Learn a b l, Backprop b, MaybeC Backprop (LStateMaybe l), Monad m)
    => (b -> m a)           -- ^ loop
    -> Int                  -- ^ num times
    -> l
    -> LParam_ I l
    -> a
    -> LState_ I l
    -> m (LState_ I l)
selfPrimeM f n l p x = fmap snd . iterateLearnM f n l p x

iterateLearnStoch
    :: (Learn a b l, Backprop b, MaybeC Backprop (LStateMaybe l), PrimMonad m)
    => (b -> m a)           -- ^ loop
    -> Int                  -- ^ num times
    -> l
    -> MWC.Gen (PrimState m)
    -> LParam_ I l
    -> a
    -> LState_ I l
    -> m ([b], LState_ I l)
iterateLearnStoch f n l g p = go 0
  where
    go !i !x !s
      | i <= n    = do
          (y , s' ) <- runLearnStoch_ l g p x s
          (ys, s'') <- flip (go (i + 1)) s' =<< f y
          pure (y : ys, s'')
      | otherwise = pure ([], s)

iterateLearnStoch_
    :: (Learn a b l, Backprop b, MaybeC Backprop (LStateMaybe l), PrimMonad m)
    => (b -> m a)           -- ^ loop
    -> Int                  -- ^ num times
    -> l
    -> MWC.Gen (PrimState m)
    -> LParam_ I l
    -> a
    -> LState_ I l
    -> m [b]
iterateLearnStoch_ f n l g p x = fmap fst . iterateLearnStoch f n l g p x

scanLearn
    :: (Learn a b l, Backprop b, MaybeC Backprop (LStateMaybe l), Traversable t)
    => l
    -> LParam_ I l
    -> t a
    -> LState_ I l
    -> (t b, LState_ I l)
scanLearn l p = runState . traverse (state . runLearn_ l p)

scanLearn_
    :: (Learn a b l, Backprop b, MaybeC Backprop (LStateMaybe l), Traversable t)
    => l
    -> LParam_ I l
    -> t a
    -> LState_ I l
    -> t b
scanLearn_ l p xs = fst . scanLearn l p xs

primeLearn
    :: (Learn a b l, Backprop b, MaybeC Backprop (LStateMaybe l), Foldable t)
    => l
    -> LParam_ I l
    -> t a
    -> LState_ I l
    -> LState_ I l
primeLearn l p = execState . traverse_ (state . runLearn_ l p)

scanLearnStoch
    :: (Learn a b l, Backprop b, MaybeC Backprop (LStateMaybe l), PrimMonad m, Traversable t)
    => l
    -> MWC.Gen (PrimState m)
    -> LParam_ I l
    -> t a
    -> LState_ I l
    -> m (t b, LState_ I l)
scanLearnStoch l g p = runStateT . traverse (StateT . runLearnStoch_ l g p)

scanLearnStoch_
    :: (Learn a b l, Backprop b, MaybeC Backprop (LStateMaybe l), PrimMonad m, Traversable t)
    => l
    -> MWC.Gen (PrimState m)
    -> LParam_ I l
    -> t a
    -> LState_ I l
    -> m (t b)
scanLearnStoch_ l g p xs = fmap fst . scanLearnStoch l g p xs

primeLearnStoch
    :: (Learn a b l, Backprop b, MaybeC Backprop (LStateMaybe l), PrimMonad m, Foldable t)
    => l
    -> MWC.Gen (PrimState m)
    -> LParam_ I l
    -> t a
    -> LState_ I l
    -> m (LState_ I l)
primeLearnStoch l g p = execStateT . traverse_ (StateT . runLearnStoch_ l g p)

initParamMaybe
    :: forall l m d proxy.
     ( MaybeC Initialize (LParamMaybe l)
     , ContGen d
     , PrimMonad m
     , KnownMayb (LParamMaybe l)
     )
    => proxy l                            -- ^ ignored
    -> d
    -> MWC.Gen (PrimState m)
    -> LParam_ m l
initParamMaybe _ d g = case knownMayb @(LParamMaybe l) of
                    N_            -> N_
                    J_ (_ :: P p) -> J_ $ initialize @p d g

initParam
    :: forall l m d proxy.
     ( Initialize (LParam l)
     , ContGen d
     , PrimMonad m
     )
    => proxy l                            -- ^ ignored
    -> d
    -> MWC.Gen (PrimState m)
    -> m (LParam l)
initParam _ = initialize

initParamNormalMaybe
    :: forall l m proxy.
     ( MaybeC Initialize (LParamMaybe l)
     , PrimMonad m
     , KnownMayb (LParamMaybe l)
     )
    => proxy l                            -- ^ ignored
    -> Double
    -> MWC.Gen (PrimState m)
    -> LParam_ m l
initParamNormalMaybe _ d g = case knownMayb @(LParamMaybe l) of
    N_            -> N_
    J_ (_ :: P p) -> J_ $ initializeNormal @p d g

initParamNormal
    :: forall l m proxy.
     ( Initialize (LParam l)
     , PrimMonad m
     )
    => proxy l                            -- ^ ignored
    -> Double
    -> MWC.Gen (PrimState m)
    -> m (LParam l)
initParamNormal _ = initializeNormal

encodeParam
    :: Bi.Binary (LParam l)
    => proxy l                              -- ^ ignored
    -> LParam l
    -> BSL.ByteString
encodeParam _ = Bi.encode

decodeParam
    :: Bi.Binary (LParam l)
    => proxy l                              -- ^ ignored
    -> BSL.ByteString
    -> LParam l
decodeParam _ = Bi.decode

decodeParamOrFail
    :: Bi.Binary (LParam l)
    => proxy l                              -- ^ ignored
    -> BSL.ByteString
    -> Either String (LParam l)
decodeParamOrFail _ = bimap thrd thrd . Bi.decodeOrFail

saveParam
    :: Bi.Binary (LParam l)
    => proxy l
    -> FilePath
    -> LParam l
    -> IO ()
saveParam p fp = BSL.writeFile fp . encodeParam p

loadParam
    :: Bi.Binary (LParam l)
    => proxy l
    -> FilePath
    -> IO (LParam l)
loadParam p fp = decodeParam p <$> BSL.readFile fp

loadParamOrFail
    :: Bi.Binary (LParam l)
    => proxy l
    -> FilePath
    -> IO (Either String (LParam l))
loadParamOrFail p fp = decodeParamOrFail p <$> BSL.readFile fp

thrd :: (a,b,c) -> c
thrd (_,_,z) = z

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE UndecidableSuperClasses #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE LambdaCase    #-}
{-# LANGUAGE EmptyCase #-}
{-# LANGUAGE TypeApplications #-}

{-# OPTIONS_GHC -Wall #-}

{-# OPTIONS_GHC -Wno-deprecations #-} -- for errorWithStackTrace

-- | Miscellany

module ConCat.Misc where

import GHC.Types (Constraint)

-- import Control.Arrow ((&&&))
-- import Data.Type.Equality

import Data.Typeable (Typeable,TypeRep,typeRep,Proxy(..))
import Data.Data (Data)
import Data.Monoid (Endo(..))
import Data.Semigroup (Semigroup(..))
import Data.Complex (Complex)
import GHC.Generics hiding (R)
-- import Unsafe.Coerce (unsafeCoerce)  -- for oops
import GHC.Stack (errorWithStackTrace)  -- for oops
import GHC.TypeLits

import Control.Newtype.Generics

{--------------------------------------------------------------------
    Type abbreviations
-------------------------------------------------------------------
-}

infixr 8 :^
infixl 7 :*
infixl 6 :+
infixr 1 :=>

type s :^ n = n -> s
type (:*)  = (,)
type (:+)  = Either
type (:=>) = (->)

{--------------------------------------------------------------------
    Helpers for GHC.Generics
--------------------------------------------------------------------}

-- | Operate inside a Generic1
inGeneric1 :: (Generic1 f, Generic1 g) => (Rep1 f a -> Rep1 g b) -> (f a -> g b)
inGeneric1 = to1 <~ from1
{-# INLINE inGeneric1 #-}

-- | Apply a unary function within the 'Comp1' constructor.
inComp :: (g (f a) -> g' (f' a')) -> ((g :.: f) a -> (g' :.: f') a')
inComp = Comp1 <~ unComp1
{-# INLINE inComp #-}

-- | Apply a binary function within the 'Comp1' constructor.
inComp2 :: (  g (f a)   -> g' (f' a')     -> g'' (f'' a''))
        -> ((g :.: f) a -> (g' :.: f') a' -> (g'' :.: f'') a'')
inComp2 = inComp <~ unComp1
{-# INLINE inComp2 #-}

-- TODO: phase out inComp and inComp2 in favor of inNew and inNew2.

absurdF :: V1 a -> b
absurdF = \ case
{-# INLINE absurdF #-}

-- infixr 1 +->
-- data (a +-> b) p = Fun1 { unFun1 :: a p -> b p }

-- -- TODO: resolve name conflict with tries. Using ":->:" for functors fits with
-- -- other type constructors in GHC.Generics.

-- instance Newtype ((a +-> b) t) where
--   type O ((a +-> b) t) = a t -> b t
--   pack = Fun1
--   unpack = unFun1

#if 0

{--------------------------------------------------------------------
    Evaluation
--------------------------------------------------------------------}

-- class Evalable e where
--   type ValT e
--   eval :: e -> ValT e

class PrimBasics p where
  unitP :: p ()
  pairP :: p (a :=> b :=> a :* b)

class Evalable p where eval :: p a -> a

-- TODO: Are we still using PrimBasics or Evalable?

#endif

{--------------------------------------------------------------------
    Other
--------------------------------------------------------------------}

type Unop   a = a -> a
type Binop  a = a -> Unop a
type Ternop a = a -> Binop a

infixl 1 <~
infixr 1 ~>

-- | Add pre- and post-processing
(~>) :: forall a b a' b'. (a' -> a) -> (b -> b') -> ((a -> b) -> (a' -> b'))
(f ~> h) g = h . g . f
-- (~>) = flip (<~)
{-# INLINE (~>) #-}

-- | Add post- and pre-processing
(<~) :: forall a b a' b'. (b -> b') -> (a' -> a) -> ((a -> b) -> (a' -> b'))
(h <~ f) g = h . g . f
{-# INLINE (<~) #-}

-- For SEC-style programming. I was using fmap instead, but my rules interfered.
result :: (b -> c) -> ((a -> b) -> (a -> c))
result = (.)
{-# INLINE result #-}

class    Yes0
instance Yes0

class    Yes1 a
instance Yes1 a

class    Yes2 a b
instance Yes2 a b

-- | Compose list of unary transformations
compose :: Foldable f => f (Unop a) -> Unop a
compose = appEndo . foldMap Endo
-- compose = foldr (.) id

infixr 3 `xor`

xor :: Binop Bool
xor = (/=)
{-# NOINLINE xor #-}

newtype Parity = Parity { getParity :: Bool }

instance Newtype Parity where
  type O Parity = Bool
  pack = Parity
  unpack (Parity x) = x

instance Semigroup Parity where
  Parity a <> Parity b = Parity (a `xor` b)

instance Monoid Parity where
  mempty = Parity False
  mappend = (<>)
  -- Parity a `mappend` Parity b = Parity (a `xor` b)

boolToInt :: Bool -> Int
boolToInt c = if c then 1 else 0
{-# INLINE boolToInt #-}

cond :: a -> a -> Bool -> a
cond t e i = if i then t else e
{-# INLINE cond #-}  -- later INLINE?

{--------------------------------------------------------------------
    Type level computations
--------------------------------------------------------------------}

infixr 3 &&

class    (a,b) => a && b
instance (a,b) => a && b

-- Saying (b,a) instead of (a,b) causes Oks k [a,b,c] to expand in order, oddly.
-- TODO: investigate.

infixr 3 &+&
class    (a t, b t) => (a &+& b) t
instance (a t, b t) => (a &+& b) t

class    f b a => Flip f a b
instance f b a => Flip f a b

-- • Potential superclass cycle for ‘&&’
--     one of whose superclass constraints is headed by a type variable: ‘a’
--   Use UndecidableSuperClasses to accept this

-- Same for Flip

type family FoldrC op b0 as where
  FoldrC op z '[]      = z
  FoldrC op z (a : as) = a `op` FoldrC op z as

type family MapC f us where
  MapC f '[]      = '[]
  MapC f (u : us) = f u : MapC f us

-- type Comp g f u = g (f u)
-- -- Operator applied to too few arguments: :
-- type MapC' f us = FoldrC (Comp (':) f) '[] us

type AndC   cs = FoldrC (&&) Yes0 cs
type AllC f us = AndC (MapC f us)

-- type family AndC' cs where
--   AndC' '[]      = Yes0
--   AndC' (c : cs) = c && AndC' cs

-- type family AllC f as where
--   AllC f '[]      = Yes0
--   AllC f (a : as) = f a && AllC f as

-- -- Operator applied to too few arguments: :
-- type as ++ bs = FoldrC (':) bs as

infixr 5 ++
type family as ++ bs where
  '[]      ++ bs = bs
  (a : as) ++ bs = a : as ++ bs

type family CrossWith f as bs where
  CrossWith f '[]      bs = '[]
  CrossWith f (a : as) bs = MapC (f a) bs ++ CrossWith f as bs

-- Illegal nested type family application ‘MapC (f a1) bs
--                                               ++ CrossWith f as bs’
--       (Use UndecidableInstances to permit this)

type AllC2 f as bs = AndC (CrossWith f as bs)

-- | Annotation for pseudo-function, i.e., defined by rules. During ccc
-- generation, don't split applications. TODO: maybe add an arity.
data PseudoFun = PseudoFun { pseudoArgs :: Int } deriving (Typeable,Data)

-- Alternatively, we could keep PseudoFun abstract:

-- pseudoFun :: Int -> PseudoFun
-- pseudoFun = PseudoFun

-- | Pseudo function to fool GHC's divergence checker.
oops :: String -> b
oops str = errorWithStackTrace ("Oops: "++str)
{-# NOINLINE oops #-}

--     In the use of ‘errorWithStackTrace’ (imported from GHC.Stack):
--     Deprecated: "'error' appends the call stack now"

-- When we use error, the divergence checker eliminates a lot of code early. An
-- alternative is unsafeCoerce, but it leads to terrible run-time errors. A safe
-- alternative seems to be errorWithStackTrace. Oddly, the doc for
-- errorWithStackTrace says "Deprecated: error appends the call stack now."

-- | Hack: delay inlining to thwart some of GHC's rewrites
delay :: a -> a
delay a = a
{-# INLINE [0] delay #-}

bottom :: a
-- bottom = error "bottom evaluated"
bottom = oops "bottom evaluated"
{-# NOINLINE bottom #-}

-- Convenient alternative to typeRep
typeR :: forall a. Typeable a => TypeRep
typeR = typeRep (Proxy :: Proxy a)

type R = Double -- Float

type C = Complex R

sqr :: Num a => a -> a
sqr a = a * a
{-# INLINE sqr #-}

magSqr :: Num a => a :* a -> a
magSqr (a,b) = sqr a + sqr b
{-# INLINE magSqr #-}

transpose :: (Traversable t, Applicative f) => t (f a) -> f (t a)
transpose = sequenceA

inTranspose :: (Applicative f, Traversable t, Applicative f', Traversable t')
            => (f (t a) -> t' (f' a)) -> (t (f a) -> f' (t' a))
inTranspose = transpose <~ transpose
{-# INLINE inTranspose #-}
-- inTranspose h = transpose . h . transpose

unzip :: Functor f => f (a :* b) -> f a :* f b
unzip ps = (fst <$> ps, snd <$> ps)
{-# INLINE unzip #-}

natValAt :: forall n. KnownNat n => Integer
natValAt = nat @n

-- Shorter name
nat :: forall n. KnownNat n => Integer
nat = natVal (Proxy @n)
{-# INLINE nat #-}

int :: forall n. KnownNat n => Int
int = fromIntegral (nat @n)
{-# INLINE int #-}

{--------------------------------------------------------------------
    Newtype
--------------------------------------------------------------------}

-- See <https://github.com/jcristovao/newtype-generics/pull/5>

-- Type generalization of underF from newtype-generics.
underF :: (Newtype n, Newtype n', o' ~ O n', o ~ O n, Functor f, Functor g)
       => (o -> n) -> (f n -> g n') -> (f o -> g o')
underF _ f = fmap unpack . f . fmap pack
{-# INLINE underF #-}

-- Type generalization of overF from newtype-generics.
overF :: (Newtype n, Newtype n', o' ~ O n', o ~ O n, Functor f, Functor g)
      => (o -> n) -> (f o -> g o') -> (f n -> g n')
overF _ f = fmap pack . f . fmap unpack
{-# INLINE overF #-}

inNew :: (Newtype p, Newtype q) =>
         (O p -> O q) -> (p -> q)
inNew = pack <~ unpack
{-# INLINE inNew #-}

inNew2 :: (Newtype p, Newtype q, Newtype r) =>
          (O p -> O q -> O r) -> (p -> q -> r)
inNew2 = inNew <~ unpack
{-# INLINE inNew2 #-}

-- TODO: use inNew and inNew2 in place of ad hoc versions throughout.

exNew :: (Newtype p, Newtype q) =>
         (p -> q) -> (O p -> O q)
exNew = unpack <~ pack
{-# INLINE exNew #-}

exNew2 :: (Newtype p, Newtype q, Newtype r) =>
          (p -> q -> r) -> (O p -> O q -> O r)
exNew2 = exNew <~ pack
{-# INLINE exNew2 #-}

{--------------------------------------------------------------------
    Constraint shorthands
--------------------------------------------------------------------}

#if 1
-- Experiment. Smaller Core?
type C1 (con :: u -> Constraint) a = con a
type C2 (con :: u -> Constraint) a b = (con a, con b)
type C3 (con :: u -> Constraint) a b c = (con a, con b, con c)
type C4 (con :: u -> Constraint) a b c d = (con a, con b, con c, con d)
type C5 (con :: u -> Constraint) a b c d e = (con a, con b, con c, con d, con e)
type C6 (con :: u -> Constraint) a b c d e f = (con a, con b, con c, con d, con e, con f)
#else
type C1 (con :: u -> Constraint) a = con a
type C2 con a b         = (C1 con a, con b)
type C3 con a b c       = (C2 con a b, con c)
type C4 con a b c d     = (C2 con a b, C2 con c d)
type C5 con a b c d e   = (C3 con a b c, C2 con d e)
type C6 con a b c d e f = (C3 con a b c, C3 con d e f)
#endif

max_stars_repo_name: kenranunderscore/concat

max_issues_repo_name: con-kitty/concat

max_forks_repo_name: con-kitty/concat


text:
module Main where

import Math.Probably.Sampler
import Math.Probably.RandIO
import Math.Probably.FoldingStats

import Numeric.LinearAlgebra

import Math.Probably.Unscented

import Graphics.Gnewplot.Exec
import Graphics.Gnewplot.Types
import Graphics.Gnewplot.Style
import Graphics.Gnewplot.Panels
import Graphics.Gnewplot.Instances
import Graphics.Gnewplot.Histogram

dt = 0.1

f v =  fromList [t1, x0 + dt * sin (t1) + noisex]
  where t0 = v @> 0
        x0 = v @> 1
        noiset = v @> 2
        noisex = v @> 3
        t1 = t0+dt + noiset

f' v =  fromList [t1, x0 + dt * sin (t1) ]
  where t0 = v @> 0
        x0 = v @> 1
        t1 = t0+dt

g v = fromList [t1+noiset, x1+noisex]
  where t1 = v @> 0
        x1 = v @> 1
        noiset = v @> 2
        noisex = v @> 3

g' v = fromList [t1, x1]
  where t1 = v @> 0
        x1 = v @> 1

procCov = (2><2) [0.001, 0,
                  0,     0.01]

obsCov = (2><2) [0.01, 0,
                  0,     0.01]

--obsCov = (1><1)  [0.01]

testUT = do
  let cov = (2><2) [0.1,0.08,0.08,0.1]
      mn = (2 |> [1,1])
      xform v = let [x,y] = toList v in fromList [x*x, y*y]
      mkpt = ((\[x,y] -> (x,y)) .toList)
  pts <- sample $ sequence $ replicate 1000 $ multiNormal mn cov
  let (mnU, covU, sigmapts) = unscentedTransform (mn,cov) xform
  let xformPts = map xform pts
  io $ gnuplotOnScreen $ map mkpt pts :+: map mkpt sigmapts --  :+: map mkpt xformPts
  io $ print mnU
  io $ print covU
  (mnE, covE) <- sample $ smellyTransform (multiNormal mn cov) 1000 xform
  io $ print mnE
  io $ print covE

main = runRIO $ do
     xyvs <- sample $ simDynamical 2000 f g procCov obsCov (fromList [0, 0.1])
     --let xyvs = flip map xyvs' $ \(xv,yv) -> (xv,negate yv)
--     pts <- sample $ sequence $ replicate 1000 $ lastSDESampler dt
     --w <- sample $ fmap (eulerMaruyama dt 0 a b) $ weiner dt 1
     let xhats = unscentedKalmanFilterAdditive f' g' procCov obsCov (fromList [0,0], ((2><2) [1,0,0,1])) $ map snd xyvs
     let get s f =  (s, Lines [] $ zip [(0::Double)..] $ map f xyvs)
     io $ gnuplotOnScreen $ get "y" ((@>1) . snd) :+: get "x" ((@>1) . fst)
                            :+: ("xhat", Lines [] $ zip [(0::Double)..] $ map ((@>1) . fst)  xhats)
     io $ print $ take 10 $ map snd xyvs

     testUT
--     io $ gnuplotOnScreen $ get "x" ((@>1) . fst)
--     io $ print $ runStat meanSDF $ pts

max_stars_repo_name: glutamate/probably-base

max_issues_repo_name: glutamate/probably-base

max_forks_repo_name: glutamate/probably-base


text:
module Main where

import Data.List (sort)
import Data.Maybe (isJust)
import Data.Primitive.Contiguous (PrimArray)
import Statistics.Array.Types (AscList,Quartiles(..))
import Test.Tasty (defaultMain,testGroup)
import Test.Tasty.HUnit (testCase,(@=?))
import Test.Tasty.QuickCheck ((===),(==>))
import Test.Tasty.QuickCheck (testProperty)

import qualified Data.Primitive.Contiguous as Arr
import qualified Statistics.Array as Stats
import qualified Statistics.Array.Types as Asc

main :: IO ()
main = defaultMain $ testGroup "stats"
  [ testGroup "quartiles"
    [ testCase "empty" $
      let arr = Arr.empty :: PrimArray Int
       in (Stats.quartiles <$> Asc.fromArray arr) @=? Nothing
    , testProperty "singleton" $ \x ->
      let arr = Arr.singleton x :: PrimArray Int
       in (Stats.quartiles <$> Asc.fromArray arr) === Just (Quartiles x x x)
    , testProperty "doubleton" $ \x0 y0 ->
      let (x, y) = if x0 <= y0 then (x0,y0) else (y0,x0)
          arr = Arr.doubleton x y :: PrimArray Int
       in (Stats.quartiles <$> Asc.fromArray arr) === Just (Quartiles x y y)
    , testProperty "tripleton" $ \x0 y0 z0 ->
      let (x1, y1) = if x0 <= y0 then (x0,y0) else (y0,x0)
          (x,y,z) = if z0 < x1 then (z0,x1,y1) else
                    if z0 < y1 then (x1,z0,y1) else
                      (x1,y1,z0)
          arr = Arr.tripleton x y z :: PrimArray Int
       in (Stats.quartiles <$> Asc.fromArray arr) === Just (Quartiles x y z)
    , testCase "bigun" $
      let arr = Arr.fromList nineDatums :: PrimArray Int
       in (Stats.quartiles <$> Asc.fromArray arr) @=? Just (Quartiles 3 4 8)
    ]
  , testProperty "median is second quartile" $ \xs ->
      let asc = Asc.fromList xs :: Maybe (AscList Int)
       in isJust asc ==>
          (Stats.median <$> asc) === (q2 . Stats.quartiles <$> asc)
  , testGroup "bowley skew" $
    [ testCase "smoke" $
      let arr = Arr.fromList nineDatums :: PrimArray Int
          qs = Stats.quartiles <$> Asc.fromArray arr
       in (Stats.bowleySkew <$> qs) @=? Just 0.6
    , testProperty "outlier" $ \x ->
      let arr = Arr.fromList (x:eightDatums) :: PrimArray Int
          qs = Stats.quartiles <$> Asc.fromArray arr
       in x > 9 ==>
          (Stats.bowleySkew <$> qs) === Just 0.6
    ]
  , testGroup "median of absolute deviations" $
    [ testCase "smoke" $
      let arr = Arr.fromList nineDatums :: PrimArray Int
       in (Stats.mad <$> Asc.fromArray arr) @=? Just 2
    ]
  , testGroup "list derivative" $
    [ testProperty "like naive" $ \unsortedXs ->
      let xs = sort unsortedXs
          lst' = zipWith (-) (tail xs) xs
          asc = Asc.unsafeFromAscendingList xs :: AscList Int
       in (not . null) xs ==>
          (Stats.listDerivative asc) === (Arr.fromList lst' :: PrimArray Int)
    ]
  ]

eightDatums :: [Int]
-- 0,2,3,3,4,6,8,9
eightDatums = [3,6,8,2,4,0,3,9]

nineDatums :: [Int]
-- 0,2,3,3,4,6,8,9,54
nineDatums = [3,6,8,2,4,0,3,54,9]

max_stars_repo_name: Zankoku-Okuno/array-statistics

max_issues_repo_name: Zankoku-Okuno/array-statistics

max_forks_repo_name: Zankoku-Okuno/array-statistics