text:
-- Copyright © 2016 Bart Massey
-- [This work is made available under the "MIT License".
-- Please see the file LICENSE in this distribution for
-- license details.]

-- Analysis tool for a Hardware Random Number Generator

import Control.Lens
import Data.Char (toLower)
import Data.Colour
import Data.Colour.Names
import Data.Colour.SRGB
import Data.Default.Class
import Data.List.Split (splitOn)
import Data.Maybe
import Graphics.Rendering.Chart
import Graphics.Rendering.Chart.Backend.Cairo

import Data.Complex
import qualified Data.Vector.Unboxed as V
import Numeric.FFT.Vector.Unnormalized

import Control.Monad
import Data.Bits
import qualified Data.ByteString as B
import Data.List
import System.Console.ParseArgs
import System.Directory
import System.Exit (exitSuccess)
import System.FilePath
import System.Random
import Text.Printf

defaultAnalysisDir :: FilePath
defaultAnalysisDir = "analysis"

wWindowSize :: Int
wWindowSize = 4096

analysisFile :: FilePath -> FilePath -> String -> String-> FilePath
analysisFile dir what suff ext =
    joinPath [dir, addExtension (what ++ "-" ++ suff) ext]

plotSuffix :: Maybe FileFormat -> String
plotSuffix (Just PDF) = "pdf"
plotSuffix (Just SVG) = "svg"
plotSuffix _ = error "no suffix for (non)plot format"

plotRender :: FileFormat -> String -> Renderable a -> IO ()
plotRender ff name r = do
    _ <- renderableToFile (fo_format .~ ff $ def) name r
    return ()

readSamples :: Int -> B.ByteString -> [Int]
readSamples sampleCount stuff =
    take sampleCount $ makeInts $ B.unpack stuff
    where
      makeInts [] = []
      makeInts [_] = error "odd byte array"
      makeInts (b2 : b1 : bs) =
          let i1 = fromIntegral b1
              i2 = fromIntegral b2
          in
          ((i1 `shiftL` 8) .|. i2) : makeInts bs

takeSpan :: Int -> Int -> [a] -> [a]
takeSpan start len xs =
    take len $ drop start xs

rawHist :: [Int] -> [(Int, Int)]
rawHist samples =
    map histify $ group $ sort $ samples ++ [smallest..largest]
    where
      smallest = minimum samples
      largest = maximum samples
      histify [] = error "bad hist group"
      histify (x : xs) = (x, length xs)

showHist :: [(Int, Int)] -> String
showHist bins =
    unlines $ map showBin bins
    where
      showBin (x, y) = unwords [show x, show y]

type ColourPair = (Colour Double, Colour Double)

plotTimeSeries :: ColourPair -> [Int]
               -> Renderable (LayoutPick Int Int Int)
plotTimeSeries (color, color2) samples = do
  let nSamples = length samples
  let samplePoints = zip [(1::Int)..] samples
  let zoomedPoints =
          zip [0, nSamples `div` 100 ..] $
              takeSpan (nSamples `div` 3) 100 samples
  let allPlot =
          toPlot $
          plot_points_style .~ filledCircles 0.5 (opaque color2) $
          plot_points_values .~ samplePoints $
          plot_points_title .~ "all" $
          def
  let zoomedPlot =
          toPlot $
          plot_lines_style . line_color .~ opaque color $
          plot_lines_values .~ [zoomedPoints] $
          plot_lines_title .~ "zoomed" $
          def
  let tsLayout =
          layout_x_axis .~ (laxis_title .~ "sample" $ def) $
          layout_y_axis .~ (laxis_title .~ "value" $ def) $
          layout_plots .~ [allPlot, zoomedPlot] $
          def
  layoutToRenderable tsLayout

plotSampleHist :: ColourPair -> Int -> [Int]
               -> Renderable (LayoutPick Int Int Int)
plotSampleHist (color, _) nBins samples = do
  let histPlot =
          plotBars $
          plot_bars_values .~  sampleBars ++ [(nBins, [0])] $
          plot_bars_item_styles .~ [barStyle] $
          plot_bars_spacing .~ BarsFixGap 0 0 $
          plot_bars_alignment .~ BarsLeft $
          def
  let histLayout =
          layout_x_axis .~ (laxis_title .~ "value" $ def) $
          layout_y_axis .~ (laxis_title .~ "frequency" $ def) $
          layout_plots .~ [histPlot] $
          def
  layoutToRenderable histLayout
  where
    sampleBars =
        map (\(x, y) -> (x, [y])) $ rawHist samples
    barStyle = (FillStyleSolid (opaque gray),
                Just (line_width .~ 0.05 $ line_color .~ opaque color $ def))

plotSampleDFT :: ColourPair -> [Double]
              -> Renderable (LayoutPick Double Double Double)
plotSampleDFT (color, _) dftBins = do
  let dftPlot =
          plotBars $
          plot_bars_spacing .~ BarsFixWidth 0.3 $
          plot_bars_values .~ dftBars $
          plot_bars_item_styles .~ [barStyle] $
          def
  let dftLayout =
          layout_x_axis .~ (laxis_title .~ "frequency" $ def) $
          layout_y_axis .~ (laxis_title .~ "amplitude" $ def) $
          layout_plots .~ [dftPlot] $
          def
  layoutToRenderable dftLayout
  where
    dftBars = map (\(x, y) -> (x, [y])) $ zip [0..] dftBins
    barStyle = (FillStyleSolid (opaque color), Nothing)

data EntropyMode = EntropyModeRaw | EntropyModeNormalized

entropy :: Real a => EntropyMode -> [a] -> Double
entropy mode samples =
    negate $ sum $ map binEntropy samples
    where
      weight =
          case mode of
            EntropyModeRaw -> 1.0
            EntropyModeNormalized -> realToFrac $ sum samples
      binEntropy 0 = 0
      binEntropy count =
          p * logBase 2 p
          where
            p = realToFrac count / weight

hannWindow :: Int -> [Double]
hannWindow nSamples =
    map hannFunction [0 .. nSamples - 1]
    where
      hannFunction n =
          0.5 * (1.0 - cos(2 * pi * fromIntegral n / nn))
          where
            nn = fromIntegral (nSamples - 1)

sampleDFT :: [Double] -> [Double]
sampleDFT samples =
    map magnitude $ V.toList $ run dftR2C $ V.fromList samples

data Bias = BiasDebiased | BiasNominal Int
data DFTMode = DFTModeRaw |
               DFTModeProper Int (Int -> [Double])

splitSamples :: Int -> [a] -> [[a]]
splitSamples n xs
    | length first < n = []
    | otherwise =
        first : splitSamples n rest
    where
      (first, rest) = splitAt n xs

processDFT :: Bias -> DFTMode -> [Int] -> [Double]
processDFT bias dftMode samples =
    case dftMode of
      DFTModeRaw ->
          sampleDFT dftSamples
          where
            dftStart = (nSamples - dftLength) `div` 3
            dftLength = min 10000 nSamples
            dftSamples = takeSpan dftStart dftLength normedSamples
      DFTModeProper windowSize window ->
        avgBins $ map (sampleDFT . applyWindow) $
          splitSamples windowSize normedSamples
        where
          applyWindow xs =
              zipWith (*) xs $ window windowSize
          avgBins bins =
              map average $ transpose bins
    where
      nSamples = length samples
      normedSamples =
          map norm samples
          where
            norm sample =
                (fromIntegral sample - dftDC) / fromIntegral nSamples
                where
                  dftDC =
                      case bias of
                        BiasNominal nBits ->
                            fromIntegral (2 ^ (nBits - 1) - 1 :: Integer)
                        BiasDebiased ->
                            average samples

average :: Real a => [a] -> Double
average samples =
    realToFrac (sum samples) / fromIntegral (length samples)

stdDeviation :: Real a => [a] -> Double
stdDeviation samples =
    sqrt (sum devs / fromIntegral (length devs))
    where
      devs = map (square . (`subtract` mean) . realToFrac) samples
      square x = x * x
      mean = average samples

spectralFlatness :: [Double] -> Double
spectralFlatness rDFT =
    10.0 * (gMeanDB - aMeanDB)
    where
      xDFT = tail rDFT
      nxDFT = fromIntegral (length xDFT)
      gMeanDB = sum (map (logBase 10) xDFT) / nxDFT
      aMeanDB = logBase 10 (sum xDFT) - logBase 10 nxDFT

secondMoment64 :: [Int] -> (Double, Double)
secondMoment64 samples =
    (minimum devs, maximum devs)
    where
      devs = map stdDeviation $ splitSamples 64 samples

showStats :: Int -> [Int] -> [Double] -> [Double] -> String
showStats nBits samples rDFT wDFT = unlines [
  printf "min: %d" (minimum samples),
  printf "max: %d" (maximum samples),
  printf "mean: %0.3g" (average samples),
  printf "byte-entropy: %0.3g"
      (entropyAdj * entropy EntropyModeNormalized hist),
  printf "second-moment-64s-min: %0.3g" smMin,
  printf "second-moment-64s-max: %0.3g" smMax,
  printf "spectral-entropy: %0.3g"
      (spectralEntropyAdj * entropy EntropyModeNormalized rDFT2),
  printf "spectral-flatness-db: %0.3g" (spectralFlatness rDFT),
  printf "avg-spectral-flatness-db: %0.3g" (spectralFlatness wDFT) ]
  where
    (smMin, smMax) = secondMoment64 samples
    hist = map snd $ rawHist samples
    rDFT2 = map (**2.0) $ tail rDFT
    entropyAdj =
        max 1.0 $ 8.0 / fromIntegral nBits
    spectralEntropyAdj =
        fromIntegral (max 8 nBits)  / logBase 2.0 (fromIntegral (length rDFT2))

analyze :: Maybe FileFormat -> ColourPair
        -> String -> String -> Int -> [Int]
        -> IO ()
analyze plotMode colors dir what nBits samples = do
  let rDFT = processDFT BiasDebiased DFTModeRaw samples
  let wDFT = processDFT BiasDebiased (DFTModeProper wWindowSize hannWindow) samples
  writeFile (af "stats" "txt") $
    showStats nBits samples rDFT wDFT
  case plotMode of
    Just ff -> do
      writeFile (af "hist" "txt") $ showHist $ rawHist samples
      let pr = plotRender ff
      pr (af' "ts") $ plotTimeSeries colors samples
      pr (af' "hist") $ plotSampleHist colors (2 ^ nBits) samples
      pr (af' "dft") $ plotSampleDFT colors rDFT
      pr (af' "wdft") $ plotSampleDFT colors wDFT
    Nothing -> return ()
  where
    af = analysisFile dir what
    af' name = af name (plotSuffix plotMode)

data ArgIndex = ArgIndexBitsFile
              | ArgIndexPlotFormat
              | ArgIndexAnalysisDir
              | ArgIndexSampleCount
              | ArgIndexColor
              | ArgIndexColor2
              | ArgIndexTests
                deriving (Eq, Ord, Show)

maybeReadS :: ReadS a -> String -> Maybe a
maybeReadS f s =
    case f s of
      [(v, "")] -> Just v
      _ -> Nothing

readColor :: String -> Colour Double
readColor name =
    head $ mapMaybe id [
               maybeReadS sRGB24reads name,
               readColourName name,
               error "illegal color" ]

argd :: [Arg ArgIndex]
argd = [
  Arg {
    argIndex = ArgIndexPlotFormat,
    argName = Just "plot-format",
    argAbbr = Just 'p',
    argData = argDataDefaulted "type" ArgtypeString "pdf",
    argDesc = "Plot format (\"pdf\", \"svg\", or \"none\" for just stats)" },
  Arg {
    argIndex = ArgIndexAnalysisDir,
    argName = Just "analysis-dir",
    argAbbr = Just 'a',
    argData = argDataDefaulted "dir" ArgtypeString defaultAnalysisDir,
    argDesc = "Directory for analysis results" },
  Arg {
    argIndex = ArgIndexSampleCount,
    argName = Just "sample-count",
    argAbbr = Just 's',
    argData = argDataOptional "count" ArgtypeInt,
    argDesc = "Number of samples to analyze" },
  Arg {
    argIndex = ArgIndexColor,
    argName = Just "color",
    argAbbr = Nothing,
    argData = argDataDefaulted "name/rgb" ArgtypeString "black",
    argDesc = "Plot color" },
  Arg {
    argIndex = ArgIndexColor2,
    argName = Just "secondary-color",
    argAbbr = Nothing,
    argData = argDataDefaulted "name/rgb" ArgtypeString "gray",
    argDesc = "Plot secondary color" },
  Arg {
    argIndex = ArgIndexTests,
    argName = Just "tests",
    argAbbr = Just 't',
    argData = argDataOptional "test-list" ArgtypeString,
    argDesc = "Comma-separated list of tests (\"help\" for help)" },
  Arg {
    argIndex = ArgIndexBitsFile,
    argName = Nothing,
    argAbbr = Nothing,
    argData = argDataOptional "bits-file" ArgtypeString,
    argDesc = "Random bits." } ]

testSet :: [String]
testSet = ["raw","prng","twobit","mid","mid7","low"]

testHelp :: IO ()
testHelp = do
  printf "Available Tests:\n"
  mapM_ (\t -> printf "  %s\n" t) testSet
  exitSuccess

main :: IO ()
main = do
  args <- parseArgsIO ArgsComplete argd

  let tests =
          case getArg args ArgIndexTests of
            Nothing -> testSet
            Just(ts) -> splitOn "," ts
  when (tests == ["help"]) testHelp

  let plotFormat =
          case map toLower $ getRequiredArg args ArgIndexPlotFormat of
            "none" -> Nothing
            "pdf" -> Just PDF
            "svg" -> Just SVG
            _ -> error "unrecognized plot format"

  bitsFile <- getArgStdio args ArgIndexBitsFile ReadMode
  rawSamples <- B.hGetContents bitsFile
  let nSamples = B.length rawSamples `div` 2
  let sampleCount = case getArg args ArgIndexSampleCount of
                      Nothing -> nSamples
                      Just n -> nSamples `min` n
  let samples = readSamples sampleCount rawSamples

  let analysisDir = getRequiredArg args ArgIndexAnalysisDir
  createDirectoryIfMissing True analysisDir

  let plotColor = readColor $ getRequiredArg args ArgIndexColor
  let plotColor2 = readColor $ getRequiredArg args ArgIndexColor2

  let aa name bits source =
          when (name `elem` tests) $
               analyze plotFormat (plotColor, plotColor2) analysisDir
                       name bits source

  aa "raw" 12 samples

  prng12Samples <- replicateM (length samples) (randomRIO (0, 4095) :: IO Int)
  aa "prng12" 12 prng12Samples

  prngSamples <- replicateM (length samples) (randomRIO (0, 255) :: IO Int)
  aa "prng" 8 prngSamples

  let twoBitSamples = map (.&. 0x03) samples
  aa "twobit" 2 twoBitSamples

  let midSamples = map ((.&. 0xff) . (`shiftR` 1)) samples
  aa "mid" 8 midSamples

  let mid7Samples = map ((.&. 0x7f) . (`shiftR` 2)) samples
  aa "mid7" 7 mid7Samples

  let lowSamples = map (.&. 0xff) samples
  aa "low" 8 lowSamples

  let hmSamples = map hmBits samples
                  where
                    hmBits sample =
                        let mid = (sample `shiftR` 1) .&. 0xff in
                        mid `xor` (sample .&. 0x01) `xor`
                          ((sample `shiftR` 9) .&. 0x01) `xor`
                          ((sample `shiftR` 10) .&. 0x01) `xor`
                          ((sample `shiftR` 11) .&. 0x01)
  aa "hm" 8 hmSamples

max_stars_repo_name: BartMassey/anrand

max_issues_repo_name: BartMassey/anrand

max_forks_repo_name: BartMassey/anrand


text:
module Isumi.Math.MatrixUtil
  ( imageToMatrix
  ) where

import           Codec.Picture
import           Numeric.LinearAlgebra (Matrix)
import qualified Numeric.LinearAlgebra as LA

-- | Convert image of Pixel8 to matrix
imageToMatrix :: Image Pixel8 -> Matrix Double
imageToMatrix =
  LA.fromColumns . fmap (LA.fromList . fmap fromIntegral) . imageToPixels

-- | Convert image to list of list of pixels, in column major order
imageToPixels :: Pixel a => Image a -> [[a]]
imageToPixels image =
  fmap (\c -> fmap (\r -> pixelAt image c r) [0.. height - 1]) [0..width - 1]
  where
  width = imageWidth image
  height = imageHeight image

ext: hs

max_stars_repo_name: IsumiF/fft-bruun

max_issues_repo_name: IsumiF/fft-bruun

max_forks_repo_name: IsumiF/fft-bruun


text:
module Main where

import Data.Complex

main :: IO ()
main = do
  putStrLn "hello world"

maxIter :: Int
maxIter = 10000

-- data Scale = {
--                ULCorner :: Point
--                DRCorner :: Point
--                XScale :: Double
--                YScale :: Double
--              }

-- Simple boolean test to see if a point is in the set
isInSet :: RealFloat a => a -> a -> Bool
isInSet r i = case escapeTime r i of
              Nothing -> True
              Just _ -> False

-- Escape Time of Nothing means the point never escaped
escapeTime :: RealFloat a => a -> a -> Maybe Int
escapeTime r i
    | numIter == maxIter = Nothing
    | otherwise = Just numIter
    where numIter = mandel r i

-- Figures out how many iterations a point takes to escape
-- Or stops after maxIter iterations
mandel :: RealFloat a => a -> a -> Int
mandel r i = length . takeWhile ((<= 2) . magnitude) .
    take maxIter $ iterate (\z -> z^2 + (r :+ i)) 0

max_stars_repo_name: IQubic/mandelbrot

max_issues_repo_name: IQubic/mandelbrot

max_forks_repo_name: IQubic/mandelbrot


text:
{-# LANGUAGE ApplicativeDo                        #-}
{-# LANGUAGE FlexibleContexts                     #-}
{-# LANGUAGE PartialTypeSignatures                #-}
{-# LANGUAGE RankNTypes                           #-}
{-# LANGUAGE ScopedTypeVariables                  #-}
{-# LANGUAGE TypeApplications                     #-}
{-# LANGUAGE TypeFamilies                         #-}
{-# OPTIONS_GHC -fno-warn-partial-type-signatures #-}

module Backprop.Learn.Test (
  -- * Tests
    Test
  , maxIxTest, rmseTest
  , squaredErrorTest, absErrorTest, totalSquaredErrorTest, squaredErrorTestV
  , crossEntropyTest
  -- ** Manipulate tests
  , lossTest, lmapTest
  -- * Run tests
  , testLearn, testLearnStoch, testLearnAll, testLearnStochAll
  -- ** Correlation tests
  , testLearnCov, testLearnCorr
  , testLearnStochCov, testLearnStochCorr
  ) where

import           Backprop.Learn.Loss
import           Backprop.Learn.Model
import           Control.Monad.Primitive
import           Data.Bifunctor
import           Data.Bitraversable
import           Data.Function
import           Data.Profunctor
import           Data.Proxy
import           GHC.TypeNats
import           Numeric.Backprop
import qualified Control.Foldl                as L
import qualified Numeric.LinearAlgebra        as HU
import qualified Numeric.LinearAlgebra.Static as H
import qualified System.Random.MWC            as MWC

-- TODO: support non-double results?

type Test o = o -> o -> Double

-- | Create a 'Test' from a 'Loss'
lossTest :: Loss a -> Test a
lossTest l x = evalBP (l x)

maxIxTest :: KnownNat n => Test (H.R n)
maxIxTest x y
    | match x y = 1
    | otherwise = 0
  where
    match = (==) `on` (HU.maxIndex . H.extract)

rmseTest :: forall n. KnownNat n => Test (H.R n)
rmseTest x y = H.norm_2 (x - y) / sqrt (fromIntegral (natVal (Proxy @n)))

squaredErrorTest :: Real a => Test a
squaredErrorTest x y = e * e
  where
    e = realToFrac (x - y)

absErrorTest :: Real a => Test a
absErrorTest x y = realToFrac . abs $ x - y

totalSquaredErrorTest :: (Applicative t, Foldable t, Real a) => Test (t a)
totalSquaredErrorTest x y = realToFrac (sum e)
  where
    e = do
      x' <- x
      y' <- y
      pure ((x' - y') ^ (2 :: Int))

squaredErrorTestV :: KnownNat n => Test (H.R n)
squaredErrorTestV x y = e `H.dot` e
  where
    e = x - y

crossEntropyTest :: KnownNat n => Test (H.R n)
crossEntropyTest targ res = -(log res H.<.> targ)

lmapTest
    :: (a -> b)
    -> Test b
    -> Test a
lmapTest f t x y = t (f x) (f y)

testLearn
    :: (Learn a b l, NoState l)
    => Test b
    -> l
    -> LParam_ I l
    -> a
    -> b
    -> Double
testLearn t l mp x y = t y $ runLearnStateless_ l mp x

testLearnStoch
    :: (Learn a b l, NoState l, PrimMonad m)
    => Test b
    -> l
    -> MWC.Gen (PrimState m)
    -> LParam_ I l
    -> a
    -> b
    -> m Double
testLearnStoch t l g mp x y = t y <$> runLearnStochStateless_ l g mp x

cov :: Fractional a => L.Fold (a, a) a
cov = do
    x  <- lmap fst           L.sum
    y  <- lmap snd           L.sum
    xy <- lmap (uncurry (*)) L.sum
    n  <- fromIntegral <$> L.length
    pure (xy / n - (x * y) / n / n)

corr :: Floating a => L.Fold (a, a) a
corr = do
    x  <- lmap fst           L.sum
    x2 <- lmap ((**2) . fst) L.sum
    y  <- lmap snd           L.sum
    y2 <- lmap ((**2) . snd) L.sum
    xy <- lmap (uncurry (*)) L.sum
    n  <- fromIntegral <$> L.length
    pure $ (xy / n - (x * y) / n / n)
         / sqrt ( x2 / n - (x / n)**2 )
         / sqrt ( y2 / n - (y / n)**2 )

testLearnCov
    :: (Learn a b l, NoState l, Foldable t, Fractional b)
    => l
    -> LParam_ I l
    -> t (a, b)
    -> b
testLearnCov l p = L.fold $ (lmap . first) (runLearnStateless_ l p) cov

testLearnCorr
    :: (Learn a b l, NoState l, Foldable t, Floating b)
    => l
    -> LParam_ I l
    -> t (a, b)
    -> b
testLearnCorr l p = L.fold $ (lmap . first) (runLearnStateless_ l p) corr

testLearnAll
    :: (Learn a b l, NoState l, Foldable t)
    => Test b
    -> l
    -> LParam_ I l
    -> t (a, b)
    -> Double
testLearnAll t l p = L.fold $ lmap (uncurry (testLearn t l p)) L.mean

-- newtype M m a = M { getM :: m a }
-- instance (Semigroup a, Applicative m) => Semigroup (M m a) where
--     M x <> M y = M $ liftA2 (<>) x y
-- instance (Monoid a, Applicative m) => Monoid (M m a) where
--     mappend = (<>)
--     mempty = M (pure mempty)

testLearnStochAll
    :: (Learn a b l, NoState l, PrimMonad m, Foldable t)
    => Test b
    -> l
    -> MWC.Gen (PrimState m)
    -> LParam_ I l
    -> t (a, b)
    -> m Double
testLearnStochAll t l g p = L.foldM $ L.premapM (uncurry (testLearnStoch t l g p))
                                      (L.generalize L.mean)

testLearnStochCov
    :: (Learn a b l, NoState l, PrimMonad m, Foldable t, Fractional b)
    => l
    -> MWC.Gen (PrimState m)
    -> LParam_ I l
    -> t (a, b)
    -> m b
testLearnStochCov l g p = L.foldM $ (L.premapM . flip bitraverse pure)
                                      (runLearnStochStateless_ l g p)
                                      (L.generalize cov)

testLearnStochCorr
    :: (Learn a b l, NoState l, PrimMonad m, Foldable t, Floating b)
    => l
    -> MWC.Gen (PrimState m)
    -> LParam_ I l
    -> t (a, b)
    -> m b
testLearnStochCorr l g p = L.foldM $ (L.premapM . flip bitraverse pure)
                                       (runLearnStochStateless_ l g p)
                                       (L.generalize corr)

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
module SchemeParser.ParserSpec where

import Test.Hspec

import Data.Ratio ((%))
import Data.Complex (Complex(..))

import qualified Data.Map as M
import Text.Parsec (runParser)

import SchemeParser.Types
import SchemeParser.Parser
import SchemeParser.Printer

main :: IO ()
main = hspec spec

parse :: String -> Either () LispVal
parse input = case runParser parseExpr M.empty "scheme-test" input of
  Left err -> Left ()
  Right val -> Right val

spec :: Spec
spec =
  describe "parser" $ do
    it "parses atoms" $ do
      parse "atom" `shouldBe` Right (LAtom "atom")
      parse "atom123" `shouldBe` Right (LAtom "atom123")
      parse "atom!$%" `shouldBe` Right (LAtom "atom!$%")
      parse "!atom" `shouldBe` Right (LAtom "!atom")
      parse "!123" `shouldBe` Right (LAtom "!123")
      parse "!atom123" `shouldBe` Right (LAtom "!atom123")
      parse "!123atom" `shouldBe` Right (LAtom "!123atom")
      parse "!%123$atom%$" `shouldBe` Right (LAtom "!%123$atom%$")
    it "parses strings" $ do
      parse "\"string\"" `shouldBe` Right (LString "string")
      parse "\"string123\"" `shouldBe` Right (LString "string123")
      parse "\"string!$%\"" `shouldBe` Right (LString "string!$%")
      parse "\"!string\"" `shouldBe` Right (LString "!string")
      parse "\"!123\"" `shouldBe` Right (LString "!123")
      parse "\"!string123\"" `shouldBe` Right (LString "!string123")
      parse "\"!123string\"" `shouldBe` Right (LString "!123string")
      parse "\"!%123$string%$\"" `shouldBe` Right (LString "!%123$string%$")
      parse "\"string\\\"nested\\\"123\"" `shouldBe` Right (LString "string\"nested\"123")
      parse "\"string'nested singlequote'123\"" `shouldBe` Right (LString "string'nested singlequote'123")
    it "parses bools" $ do
      parse "#t" `shouldBe` Right (LBool True)
      parse "#f" `shouldBe` Right (LBool False)
    it "parses chars" $ do
      parse "#\\a" `shouldBe` Right (LChar 'a')
      parse "#\\ " `shouldBe` Right (LChar ' ')
      parse "#\\space" `shouldBe` Right (LChar ' ')
      parse "#\\newline" `shouldBe` Right (LChar '\n')
      parse "#\\n" `shouldBe` Right (LChar 'n')
    it "parses numbers" $ do
      parse "123" `shouldBe` Right (LNumber 123)
      parse "#d123" `shouldBe` Right (LNumber 123)
      parse "#b1010" `shouldBe` Right (LNumber 10)
      parse "#o123" `shouldBe` Right (LNumber 83)
      parse "#xff" `shouldBe` Right (LNumber 255)
    it "parses floats" $ do
      parse "123.123" `shouldBe` Right (LFloat 123.123)
      parse "0.123" `shouldBe` Right (LFloat 0.123)
      parse "123.0" `shouldBe` Right (LFloat 123.0)
    it "parses ratios" $ do
      parse "10/3" `shouldBe` Right (LRatio (10 % 3))
      parse "25/5" `shouldBe` Right (LRatio (25 % 5))
      parse "25/5" `shouldBe` Right (LRatio (5 % 1))
    it "parses complex numbers" $ do
      parse "1+1i" `shouldBe` Right (LComplex (1 :+ 1))
      parse "123+456i" `shouldBe` Right (LComplex (123 :+ 456))
      -- TODO: test floats
    it "parses lists" $ do
      parse "()" `shouldBe` Right (LList [])
      parse "(atom)" `shouldBe` Right (LList [LAtom "atom"])
      parse "(one two three)" `shouldBe` Right (LList [LAtom "one", LAtom "two", LAtom "three"])
      parse "(\"foo\")" `shouldBe` Right (LList [LString "foo"])
      parse "(\"one\" \"two\" \"three\")" `shouldBe` Right (LList [LString "one", LString "two", LString "three"])
      parse "(#\\a)" `shouldBe` Right (LList [LChar 'a'])
      parse "(#\\a #\\  #\\space #\\newline #\\n)" `shouldBe` Right (LList [LChar 'a', LChar ' ', LChar ' ', LChar '\n', LChar 'n'])
      parse "(#t)" `shouldBe` Right (LList [LBool True])
      parse "(1)" `shouldBe` Right (LList [LNumber 1])
      parse "(1 2 3)" `shouldBe` Right (LList [LNumber 1, LNumber 2, LNumber 3])
      parse "(#d1)" `shouldBe` Right (LList [LNumber 1])
      parse "(#d1 #d2 #d3)" `shouldBe` Right (LList [LNumber 1, LNumber 2, LNumber 3])
      parse "(#o1)" `shouldBe` Right (LList [LNumber 1])
      parse "(#o1 #o2 #o3)" `shouldBe` Right (LList [LNumber 1, LNumber 2, LNumber 3])
      parse "(#b1)" `shouldBe` Right (LList [LNumber 1])
      parse "(#b1 #b10 #b11)" `shouldBe` Right (LList [LNumber 1, LNumber 2, LNumber 3])
      parse "(#x1)" `shouldBe` Right (LList [LNumber 1])
      parse "(#x1 #x2 #x3)" `shouldBe` Right (LList [LNumber 1, LNumber 2, LNumber 3])
      parse "(12.34)" `shouldBe` Right (LList [LFloat 12.34])
      parse "(12.34 56.78 90.12)" `shouldBe` Right (LList [LFloat 12.34, LFloat 56.78, LFloat 90.12])
      parse "(12/34)" `shouldBe` Right (LList [LRatio (12 % 34)])
      parse "(12/34 56/78 90/12)" `shouldBe` Right (LList [LRatio (12 % 34), LRatio (56 % 78), LRatio (90 % 12)])
      parse "(12+34i)" `shouldBe` Right (LList [LComplex (12 :+ 34)])
      parse "(12+34i 56+78i 90+12i)" `shouldBe` Right (LList [LComplex (12 :+ 34), LComplex (56 :+ 78), LComplex (90 :+ 12)])
      -- TODO: test complex with floats
      -- TODO: test more stuff

max_stars_repo_name: alexpeits/haskell-lab-old

max_issues_repo_name: alexpeits/haskell-lab-old

max_forks_repo_name: alexpeits/haskell-lab-old


text:
{-# LANGUAGE BangPatterns               #-}
{-# LANGUAGE CPP                        #-}
{-# LANGUAGE FlexibleInstances          #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE NoImplicitPrelude          #-}
{-# LANGUAGE ScopedTypeVariables        #-}
{-# LANGUAGE StandaloneDeriving         #-}
{-# LANGUAGE TypeOperators              #-}

#if __GLASGOW_HASKELL__ >=706
{-# LANGUAGE PolyKinds #-}
#endif

#if __GLASGOW_HASKELL__ >=702
{-# LANGUAGE Trustworthy                #-}
#endif

-- | A class of non-empty data structures that can be folded to a summary value.
module Data.Foldable1 (
    Foldable1(..),
    foldr1, foldr1',
    foldl1, foldl1',
    intercalate1,
    foldrM1,
    foldlM1,
    foldrMapM1,
    foldlMapM1,
    maximumBy,
    minimumBy,
    ) where

import Data.Foldable      (Foldable, foldlM, foldr)
import Data.List          (foldl, foldl')
import Data.List.NonEmpty (NonEmpty (..))
import Data.Semigroup
       (Dual (..), First (..), Last (..), Max (..), Min (..), Product (..),
       Semigroup (..), Sum (..))
import Prelude
       (Maybe (..), Monad (..), Ord, Ordering (..), id, seq, ($!), ($), (.),
       (=<<), flip, const)

import qualified Data.List.NonEmpty as NE

#if MIN_VERSION_base(4,4,0)
import Data.Complex (Complex (..))
import GHC.Generics
       (M1 (..), Par1 (..), Rec1 (..), V1, (:*:) (..), (:+:) (..), (:.:) (..))
import Prelude      (error)
#endif

#if MIN_VERSION_base(4,6,0)
import Data.Ord (Down (..))
#endif

#if MIN_VERSION_base(4,8,0)
import qualified Data.Monoid as Mon
#endif

#if !MIN_VERSION_base(4,12,0)
import Data.Orphans ()
#endif

#ifdef MIN_VERSION_tagged
import Data.Tagged (Tagged (..))
#endif

-- Instances
import Control.Applicative.Backwards (Backwards (..))
import Control.Applicative.Lift      (Lift (..))
import Control.Monad.Trans.Identity  (IdentityT (..))
import Data.Functor.Compose          (Compose (..))
import Data.Functor.Identity         (Identity (..))
import Data.Functor.Reverse          (Reverse (..))
import Data.Tree                     (Tree (..))

import qualified Data.Functor.Product as Functor
import qualified Data.Functor.Sum     as Functor

-- coerce
#if __GLASGOW_HASKELL__ <708
import Unsafe.Coerce (unsafeCoerce)
#else
import Data.Coerce (Coercible, coerce)
#endif

-- $setup
-- >>> import Prelude hiding (foldr1, foldl1, head, last, minimum, maximum)

-------------------------------------------------------------------------------
-- Foldable1 type class
-------------------------------------------------------------------------------

-- | Non-empty data structures that can be folded.
class Foldable t => Foldable1 t where
#if __GLASGOW_HASKELL__ >= 708
    {-# MINIMAL foldMap1 | foldrMap1 #-}
#endif

    -- At some point during design it was possible to define this class using
    -- only 'toNonEmpty'. But it seems a bad idea in general.
    --
    -- So currently we require either foldMap1 or foldrMap1
    --
    -- * foldMap1 defined using foldrMap1
    -- * foldrMap1 defined using foldMap1
    --
    -- One can alsays define instance using following pattern:
    --
    --     toNonEmpty = ...
    --     foldMap f     = foldMap f     . toNonEmpty
    --     foldrMap1 f g = foldrMap1 f g . toNonEmpty

    -- | Combine the elements of a structure using a semigroup.
    fold1 :: Semigroup m => t m -> m
    fold1 = foldMap1 id

    -- | Map each element of the structure to a semigroup,
    -- and combine the results.
    --
    -- >>> foldMap1 Sum (1 :| [2, 3, 4])
    -- Sum {getSum = 10}
    --
    foldMap1 :: Semigroup m => (a -> m) -> t a -> m
    foldMap1 f = foldrMap1 f (\a m -> f a <> m)

    -- | A variant of 'foldMap1' that is strict in the accumulator.
    --
    -- >>> foldMap1' Sum (1 :| [2, 3, 4])
    -- Sum {getSum = 10}
    --
    foldMap1' :: Semigroup m => (a -> m) -> t a -> m
    foldMap1' f = foldlMap1' f (\m a -> m <> f a)

    -- | List of elements of a structure, from left to right.
    --
    -- >>> toNonEmpty (Identity 2)
    -- 2 :| []
    --
    toNonEmpty :: t a -> NonEmpty a
    toNonEmpty = runNonEmptyDList . foldMap1 singleton

    -- | The largest element of a non-empty structure.
    --
    -- >>> maximum (32 :| [64, 8, 128, 16])
    -- 128
    --
    maximum :: Ord a => t a -> a
    maximum = getMax #. foldMap1' Max

    -- | The least element of a non-empty structure.
    --
    -- >>> minimum (32 :| [64, 8, 128, 16])
    -- 8
    --
    minimum :: Ord a => t a -> a
    minimum = getMin #. foldMap1' Min

    -- | The first element of a non-empty structure.
    --
    -- >>> head (1 :| [2, 3, 4])
    -- 1
    --
    head :: t a -> a
    head = getFirst #. foldMap1 First

    -- | The last element of a non-empty structure.
    --
    -- >>> last (1 :| [2, 3, 4])
    -- 4
    --
    last :: t a -> a
    last = getLast #. foldMap1 Last

    -- | Generalized 'foldr1'.
    foldrMap1 :: (a -> b) -> (a -> b -> b) -> t a -> b
    foldrMap1 f g xs =
        appFromMaybe (foldMap1 (FromMaybe #. h) xs) Nothing
      where
        h a Nothing  = f a
        h a (Just b) = g a b

    -- | Generalized 'foldl1''.
    foldlMap1' :: (a -> b) -> (b -> a -> b) -> t a -> b
    foldlMap1' f g xs =
        foldrMap1 f' g' xs SNothing
      where
        -- f' :: a -> SMaybe b -> b
        f' a SNothing  = f a
        f' a (SJust b) = g b a

        -- g' :: a -> (SMaybe b -> b) -> SMaybe b -> b
        g' a x SNothing  = x $! SJust (f a)
        g' a x (SJust b) = x $! SJust (g b a)

    -- | Generalized 'foldl1'.
    foldlMap1 :: (a -> b) -> (b -> a -> b) -> t a -> b
    foldlMap1 f g xs =
        appFromMaybe (getDual (foldMap1 ((Dual . FromMaybe) #. h) xs)) Nothing
      where
        h a Nothing  = f a
        h a (Just b) = g b a

    -- | Generalized 'foldr1''.
    foldrMap1' :: (a -> b) -> (a -> b -> b) -> t a -> b
    foldrMap1' f g xs =
        foldlMap1 f' g' xs SNothing
      where
        f' a SNothing  = f a
        f' a (SJust b) = g a b

        g' bb a SNothing  = bb $! SJust (f a)
        g' bb a (SJust b) = bb $! SJust (g a b)

-------------------------------------------------------------------------------
-- Combinators
-------------------------------------------------------------------------------

-- | Right-associative fold of a structure.
--
-- In the case of lists, 'foldr1', when applied to a binary operator,
-- and a list, reduces the list using the binary operator,
-- from right to left:
--
-- > foldr1 f [x1, x2, ..., xn] == x1 `f` (x2 `f` ... (xn1 `f` xn )...)
--
-- Note that, since the head of the resulting expression is produced by
-- an application of the operator to the first element of the list,
-- 'foldr1' can produce a terminating expression from an infinite list.
--
-- For a general 'Foldable1' structure this should be semantically identical
-- to,
--
-- @foldr1 f = foldr1 f . 'toNonEmpty'@
--
foldr1 :: Foldable1 t => (a -> a -> a) -> t a -> a
foldr1 = foldrMap1 id
{-# INLINE foldr1 #-}

-- | Right-associative fold of a structure, but with strict application of
-- the operator.
--
foldr1' :: Foldable1 t => (a -> a -> a) -> t a -> a
foldr1' = foldrMap1' id
{-# INLINE foldr1' #-}

-- | Left-associative fold of a structure.
--
-- In the case of lists, 'foldl1', when applied to a binary
-- operator, and a list, reduces the list using the binary operator,
-- from left to right:
--
-- > foldl1 f [x1, x2, ..., xn] == (...((x1 `f` x2) `f`...) `f` xn
--
-- Note that to produce the outermost application of the operator the
-- entire input list must be traversed. This means that 'foldl1' will
-- diverge if given an infinite list.
--
-- Also note that if you want an efficient left-fold, you probably want to
-- use 'foldl1'' instead of 'foldl1'. The reason for this is that latter does
-- not force the "inner" results (e.g. @x1 \`f\` x2@ in the above example)
-- before applying them to the operator (e.g. to @(\`f\` x3)@). This results
-- in a thunk chain \(\mathcal{O}(n)\) elements long, which then must be
-- evaluated from the outside-in.
--
-- For a general 'Foldable1' structure this should be semantically identical
-- to,
--
-- @foldl1 f z = foldl1 f . 'toNonEmpty'@
--
foldl1 :: Foldable1 t => (a -> a -> a) -> t a -> a
foldl1 = foldlMap1 id
{-# INLINE foldl1 #-}

-- | Left-associative fold of a structure but with strict application of
-- the operator.
--
-- This ensures that each step of the fold is forced to weak head normal
-- form before being applied, avoiding the collection of thunks that would
-- otherwise occur. This is often what you want to strictly reduce a finite
-- list to a single, monolithic result (e.g. 'length').
--
-- For a general 'Foldable1' structure this should be semantically identical
-- to,
--
-- @foldl1' f z = foldl1 f . 'toNonEmpty'@
--
foldl1' :: Foldable1 t => (a -> a -> a) -> t a -> a
foldl1' = foldlMap1' id
{-# INLINE foldl1' #-}

-- | Insert an @m@ between each pair of @t m@.
--
-- >>> intercalate1 ", " $ "hello" :| ["how", "are", "you"]
-- "hello, how, are, you"
--
-- >>> intercalate1 ", " $ "hello" :| []
-- "hello"
--
-- >>> intercalate1 mempty $ "I" :| ["Am", "Fine", "You?"]
-- "IAmFineYou?"
--
intercalate1 :: (Foldable1 t, Semigroup m) => m -> t m -> m
intercalate1 = flip intercalateMap1 id

intercalateMap1 :: (Foldable1 t, Semigroup m) => m -> (a -> m) -> t a -> m
intercalateMap1 j f = flip joinee j . foldMap1 (JoinWith . const . f)

-- | Monadic fold over the elements of a non-empty structure,
-- associating to the right, i.e. from right to left.
foldrM1 :: (Foldable1 t, Monad m) => (a -> a -> m a) -> t a -> m a
foldrM1 = foldrMapM1 return

-- | Map variant of 'foldrM1'.
foldrMapM1 :: (Foldable1 t, Monad m) => (a -> m b) -> (a -> b -> m b) -> t a -> m b
foldrMapM1 g f = go . toNonEmpty
  where
    go (e:|es) =
      case es of
        []   -> g e
        x:xs -> f e =<< go (x:|xs)

-- | Monadic fold over the elements of a non-empty structure,
-- associating to the left, i.e. from left to right.
foldlM1 :: (Foldable1 t, Monad m) => (a -> a -> m a) -> t a -> m a
foldlM1 = foldlMapM1 return

-- | Map variant of 'foldlM1'.
foldlMapM1 :: (Foldable1 t, Monad m) => (a -> m b) -> (b -> a -> m b) -> t a -> m b
foldlMapM1 g f t = g x >>= \y -> foldlM f y xs
  where x:|xs = toNonEmpty t

-- | The largest element of a non-empty structure with respect to the
-- given comparison function.

-- See Note [maximumBy/minimumBy space usage]
maximumBy :: Foldable1 t => (a -> a -> Ordering) -> t a -> a
maximumBy cmp = foldl1' max'
  where max' x y = case cmp x y of
                        GT -> x
                        _  -> y

-- | The least element of a non-empty structure with respect to the
-- given comparison function.

-- See Note [maximumBy/minimumBy space usage]
minimumBy :: Foldable1 t => (a -> a -> Ordering) -> t a -> a
minimumBy cmp = foldl1' min'
  where min' x y = case cmp x y of
                        GT -> y
                        _  -> x

-------------------------------------------------------------------------------
-- Auxiliary types
-------------------------------------------------------------------------------

-- | Used for default toNonEmpty implementation.
newtype NonEmptyDList a = NEDL { unNEDL :: [a] -> NonEmpty a }

instance Semigroup (NonEmptyDList a) where
  xs <> ys = NEDL (unNEDL xs . NE.toList . unNEDL ys)
  {-# INLINE (<>) #-}

-- | Create dlist with a single element
singleton :: a -> NonEmptyDList a
singleton = NEDL #. (:|)

-- | Convert a dlist to a non-empty list
runNonEmptyDList :: NonEmptyDList a -> NonEmpty a
runNonEmptyDList = ($[]) . unNEDL
{-# INLINE runNonEmptyDList #-}

-- | Used for foldrMap1 and foldlMap1 definitions
newtype FromMaybe b = FromMaybe { appFromMaybe :: Maybe b -> b }

instance Semigroup (FromMaybe b) where
    FromMaybe f <> FromMaybe g = FromMaybe (f . Just . g)

-- | Strict maybe, used to implement default foldlMap1' etc.
data SMaybe a = SNothing | SJust !a

-- | Used to implement intercalate1/Map
newtype JoinWith a = JoinWith {joinee :: (a -> a)}

instance Semigroup a => Semigroup (JoinWith a) where
  JoinWith a <> JoinWith b = JoinWith $ \j -> a j <> j <> b j

-------------------------------------------------------------------------------
-- Instances for misc base types
-------------------------------------------------------------------------------

instance Foldable1 NonEmpty where
    foldMap1 f (x :| xs) = go (f x) xs where
        go y [] = y
        go y (z : zs) = y <> go (f z) zs

    foldMap1' f (x :| xs) = foldl' (\m y -> m <> f y) (f x) xs

    toNonEmpty = id

    foldrMap1 g f (x :| xs) = go x xs where
        go y [] = g y
        go y (z : zs) = f y (go z zs)

    foldlMap1  g f (x :| xs) = foldl f (g x) xs
    foldlMap1' g f (x :| xs) = let gx = g x in gx `seq` foldl' f gx xs

    head = NE.head
    last = NE.last

#if MIN_VERSION_base(4,6,0)
instance Foldable1 Down where
    foldMap1 = coerce
#endif

#if MIN_VERSION_base(4,4,0)
instance Foldable1 Complex where
    foldMap1 f (x :+ y) = f x <> f y

    toNonEmpty (x :+ y) = x :| y : []
#endif

-------------------------------------------------------------------------------
-- Instances for tuples
-------------------------------------------------------------------------------

-- 3+ tuples are not Foldable/Traversable

instance Foldable1 ((,) a) where
    foldMap1 f (_, y) = f y
    toNonEmpty (_, x) = x :| []
    minimum (_, x) = x
    maximum (_, x) = x
    head (_, x) = x
    last (_, x) = x

-------------------------------------------------------------------------------
-- Monoid / Semigroup instances
-------------------------------------------------------------------------------

instance Foldable1 Dual where
    foldMap1 = coerce

instance Foldable1 Sum where
    foldMap1 = coerce

instance Foldable1 Product where
    foldMap1 = coerce

instance Foldable1 Min where
    foldMap1 = coerce

instance Foldable1 Max where
    foldMap1 = coerce

instance Foldable1 First where
    foldMap1 = coerce

instance Foldable1 Last where
    foldMap1 = coerce

#if MIN_VERSION_base(4,8,0)
deriving instance (Foldable1 f) => Foldable1 (Mon.Alt f)
#endif

#if MIN_VERSION_base(4,12,0)
deriving instance (Foldable1 f) => Foldable1 (Mon.Ap f)
#endif

-------------------------------------------------------------------------------
-- GHC.Generics instances
-------------------------------------------------------------------------------

#if MIN_VERSION_base(4,4,0)
instance Foldable1 V1 where
    foldMap1 _ x = x `seq` error "foldMap1 @V1"

instance Foldable1 Par1 where
    foldMap1 = coerce

deriving instance Foldable1 f => Foldable1 (Rec1 f)

deriving instance Foldable1 f => Foldable1 (M1 i c f)

instance (Foldable1 f, Foldable1 g) => Foldable1 (f :+: g) where
    foldMap1 f (L1 x) = foldMap1 f x
    foldMap1 f (R1 y) = foldMap1 f y

instance (Foldable1 f, Foldable1 g) => Foldable1 (f :*: g) where
    foldMap1 f (x :*: y) = foldMap1 f x <> foldMap1 f y

instance (Foldable1 f, Foldable1 g) => Foldable1 (f :.: g) where
    foldMap1 f = foldMap1 (foldMap1 f) . unComp1
#endif

-------------------------------------------------------------------------------
-- Extra instances
-------------------------------------------------------------------------------

instance Foldable1 Identity where
    foldMap1      = coerce

    foldrMap1  g _ = coerce g
    foldrMap1' g _ = coerce g
    foldlMap1  g _ = coerce g
    foldlMap1' g _ = coerce g

    toNonEmpty (Identity x) = x :| []

    last    = coerce
    head    = coerce
    minimum = coerce
    maximum = coerce

-- | It would be enough for either half of a product to be 'Foldable1'.
-- Other could be 'Foldable'.
instance (Foldable1 f, Foldable1 g) => Foldable1 (Functor.Product f g) where
    foldMap1 f (Functor.Pair x y)    = foldMap1 f x <> foldMap1 f y
    foldrMap1 g f (Functor.Pair x y) = foldr f (foldrMap1 g f y) x

    head (Functor.Pair x _) = head x
    last (Functor.Pair _ y) = last y

instance (Foldable1 f, Foldable1 g) => Foldable1 (Functor.Sum f g) where
    foldMap1 f (Functor.InL x) = foldMap1 f x
    foldMap1 f (Functor.InR y) = foldMap1 f y

    foldrMap1 g f (Functor.InL x) = foldrMap1 g f x
    foldrMap1 g f (Functor.InR y) = foldrMap1 g f y

    toNonEmpty (Functor.InL x) = toNonEmpty x
    toNonEmpty (Functor.InR y) = toNonEmpty y

    head (Functor.InL x) = head x
    head (Functor.InR y) = head y
    last (Functor.InL x) = last x
    last (Functor.InR y) = last y

    minimum (Functor.InL x) = minimum x
    minimum (Functor.InR y) = minimum y
    maximum (Functor.InL x) = maximum x
    maximum (Functor.InR y) = maximum y

instance (Foldable1 f, Foldable1 g) => Foldable1 (Compose f g) where
    foldMap1 f = foldMap1 (foldMap1 f) . getCompose

    foldrMap1 f g = foldrMap1 (foldrMap1 f g) (\xs x -> foldr g x xs) . getCompose

    head = head . head . getCompose
    last = last . last . getCompose

-------------------------------------------------------------------------------
-- containers
-------------------------------------------------------------------------------

instance Foldable1 Tree where
    foldMap1 f (Node x [])       = f x
    foldMap1 f (Node x (y : ys)) = f x <> foldMap1 (foldMap1 f) (y :| ys)

    foldMap1' f = go where
        go (Node x ys) =
            foldl' (\m zs -> let gozs = go zs in gozs `seq` m <> gozs) (f x) ys

    foldlMap1 f g (Node x xs) = goForest (f x) xs where
        goForest = foldl' go
        go y (Node z zs) = goForest (g y z) zs

    foldlMap1' f g (Node x xs) = goForest (f x) xs where
        goForest !y = foldl' go y
        go !y (Node z zs) = goForest (g y z) zs

    head (Node x _) = x

-------------------------------------------------------------------------------
-- transformers
-------------------------------------------------------------------------------

instance Foldable1 f => Foldable1 (Reverse f) where
    foldMap1 f = getDual . foldMap1 (Dual . f) . getReverse

    foldrMap1  f g (Reverse xs) = foldlMap1  f (flip g) xs
    foldlMap1  f g (Reverse xs) = foldrMap1  f (flip g) xs
    foldrMap1' f g (Reverse xs) = foldlMap1' f (flip g) xs
    foldlMap1' f g (Reverse xs) = foldrMap1' f (flip g) xs

    head = last . getReverse
    last = head . getReverse

deriving instance Foldable1 f => Foldable1 (IdentityT f)

instance Foldable1 f => Foldable1 (Backwards f) where
    foldMap1 f = foldMap1 f . forwards

instance Foldable1 f => Foldable1 (Lift f) where
    foldMap1 f (Pure x)  = f x
    foldMap1 f (Other y) = foldMap1 f y

-------------------------------------------------------------------------------
-- tagged
-------------------------------------------------------------------------------

#ifdef MIN_VERSION_tagged
instance Foldable1 (Tagged b) where
    foldMap1      = coerce

    foldrMap1  g _ = coerce g
    foldrMap1' g _ = coerce g
    foldlMap1  g _ = coerce g
    foldlMap1' g _ = coerce g

    toNonEmpty x = coerce x :| []

    last    = coerce
    head    = coerce
    minimum = coerce
    maximum = coerce
#endif

-------------------------------------------------------------------------------
-- coerce shim
-------------------------------------------------------------------------------

#if __GLASGOW_HASKELL__ <708
coerce :: a -> b
coerce = unsafeCoerce

(#.) :: (b -> c) -> (a -> b) -> a -> c
(#.) _f = coerce
#else
(#.) :: Coercible b c => (b -> c) -> (a -> b) -> a -> c
(#.) _f = coerce
#endif


text:
import Paraiso
import Data.Complex
import System.Environment

main = do
  args <- getArgs
  let arch = if "--cuda" `elem` args then
               CUDA 128 128
             else
               X86
  putStrLn $ compile arch code
    where
      code = do
         parallel 16384 $ do
           z <- allocate
           z0<- allocate
           z =$ (Rand (-2.0) 2.0) :+ (Rand (-2.0) 2.0)
           z0=$ (z:: Complex (Expr Double))
           cuda $ do
             sequential 65536 $ do
               z =$ z * z - 1
           output [realPart z0,imagPart z0, realPart z, imagPart z]


text:
{-# LANGUAGE DataKinds            #-}
{-# LANGUAGE FlexibleContexts     #-}
{-# LANGUAGE GADTs                #-}
{-# LANGUAGE InstanceSigs         #-}
{-# LANGUAGE KindSignatures       #-}
{-# LANGUAGE LambdaCase           #-}
{-# LANGUAGE RankNTypes           #-}
{-# LANGUAGE ScopedTypeVariables  #-}
{-# LANGUAGE TypeApplications     #-}
{-# LANGUAGE TypeFamilies         #-}
{-# LANGUAGE UndecidableInstances #-}

module TensorOps.BLAS.HMat
  ( HMat
  , HMatD
  ) where

import           Control.DeepSeq
import           Data.Kind
import           Data.Singletons
import           Data.Singletons.TypeLits
import           Data.Type.Combinator
import           Data.Type.Vector            (Vec, VecT(..))
import           Data.Type.Vector.Util       (curryV2', curryV3')
import           Numeric.LinearAlgebra
import           Numeric.LinearAlgebra.Data  as LA
import           Numeric.LinearAlgebra.Devel
import           TensorOps.BLAS
import           Type.Class.Higher
import           Type.Class.Higher.Util
import qualified Data.Finite                 as DF
import qualified Data.Finite.Internal        as DF
import qualified Data.Vector.Storable        as VS

type HMatD = HMat Double

data HMat :: Type -> BShape Nat -> Type where
    HMV :: { unHMV :: !(Vector a) } -> HMat a ('BV n)
    HMM :: { unHMM :: !(Matrix a) } -> HMat a ('BM n m)

instance (VS.Storable a, Show a, Element a) => Show (HMat a s) where
    showsPrec p = \case
      HMV x -> showParen (p > 10) $ showString "HMV "
                                  . showsPrec 11 x
      HMM x -> showParen (p > 10) $ showString "HMM "
                                  . showsPrec 11 x

instance (VS.Storable a, Show a, Element a) => Show1 (HMat a)

instance (VS.Storable a, NFData a) => NFData (HMat a s) where
    rnf = \case
      HMV xs -> rnf xs
      HMM xs -> rnf xs
    {-# INLINE rnf #-}

instance (VS.Storable a, NFData a) => NFData1 (HMat a)

instance (SingI s, Container Vector a, Container Matrix a, Num a) => Num (HMat a s) where
    (+) = unsafeZipH add add
    (*) = unsafeZipH (VS.zipWith (*)) (liftMatrix2 (VS.zipWith (*)))
    (-) = unsafeZipH (VS.zipWith (-)) (liftMatrix2 (VS.zipWith (-)))
    negate = unsafeMapH (scale (-1)) (scale (-1))
    abs    = unsafeMapH (cmap abs) (cmap abs)
    signum = unsafeMapH (cmap signum) (cmap signum)
    fromInteger = case (sing :: Sing s) of
        SBV n   -> HMV . flip konst (fromIntegral (fromSing n)) . fromInteger
        SBM n m -> HMM . flip konst (fromIntegral (fromSing n)
                                    ,fromIntegral (fromSing m)
                                    ) . fromInteger

-- | WARNING!! Functions should assume equal sized inputs and return
-- outputs of the same size!  This is not checked!!!
unsafeZipH
    :: (Vector a -> Vector a -> Vector a)
    -> (Matrix a -> Matrix a -> Matrix a)
    -> HMat a s -> HMat a s -> HMat a s
unsafeZipH f g = \case
    HMV x -> \case
      HMV y -> HMV $ f x y
    HMM x -> \case
      HMM y -> HMM $ g x y

-- | WARNING!! Functions should return outputs of the same size!  This is
-- not checked!!!
unsafeMapH
    :: (Vector a -> Vector a)
    -> (Matrix a -> Matrix a)
    -> HMat a s -> HMat a s
unsafeMapH f g = \case
    HMV x -> HMV $ f x
    HMM x -> HMM $ g x

liftB'
    :: (Numeric a)
    => Sing s
    -> (Vec n a -> a)
    -> Vec n (HMat a s)
    -> HMat a s
liftB' s f xs = bgen s $ \i -> f (indexB i <$> xs)
{-# INLINE liftB' #-}

instance (Container Vector a, Numeric a) => BLAS (HMat a) where
    type ElemB (HMat a) = a

    -- TODO: rewrite rules
    -- write in parallel?
    liftB
        :: forall n s. ()
        => Sing s
        -> (Vec n a -> a)
        -> Vec n (HMat a s)
        -> HMat a s
    liftB s f = \case
        ØV -> case s of
          SBV sN    -> HMV $ konst (f ØV) ( fromIntegral (fromSing sN) )
          SBM sN sM -> HMM $ konst (f ØV) ( fromIntegral (fromSing sN)
                                          , fromIntegral (fromSing sM)
                                          )
        I x :* ØV -> case x of
          HMV x' -> HMV (cmap (f . (:* ØV) . I) x')
          HMM x' -> HMM (cmap (f . (:* ØV) . I) x')
        I x :* I y :* ØV -> case x of
          HMV x' -> case y of
            HMV y' -> HMV $ VS.zipWith (curryV2' f) x' y'
          HMM x' -> case y of
            HMM y' -> HMM $ liftMatrix2 (VS.zipWith (curryV2' f)) x' y'
        xs@(I x :* I y :* I z :* ØV) -> case x of
          HMV x' -> case y of
            HMV y' -> case z of
              HMV z' -> HMV $ VS.zipWith3 (curryV3' f) x' y' z'
          _ -> liftB' s f xs
        xs@(_ :* _ :* _ :* _ :* _) -> liftB' s f xs

    axpy α (HMV x) my
        = HMV
        . maybe id (add . unHMV) my
        . scale α
        $ x
    {-# INLINE axpy #-}
    dot (HMV x) (HMV y)
        = x <.> y
    {-# INLINE dot #-}
    ger (HMV x) (HMV y)
        = HMM $ x `outer` y
    {-# INLINE ger #-}
    gemv α (HMM a) (HMV x) mβy
        = HMV
        . maybe id (\(β, HMV y) -> add (scale β y)) mβy
        . (a #>)
        . scale α
        $ x
    {-# INLINE gemv #-}
    gemm α (HMM a) (HMM b) mβc
        = HMM
        . maybe id (\(β, HMM c) -> add (scale β c)) mβc
        . (a <>)
        . scale α
        $ b
    {-# INLINE gemm #-}
    scaleB α = unsafeMapH (scale α) (scale α)
    {-# INLINE scaleB #-}
    addB = unsafeZipH add add
    {-# INLINE addB #-}
    indexB = \case
        PBV i -> \case
          HMV x -> x `atIndex` fromInteger (DF.getFinite i)
        PBM i j -> \case
          HMM x -> x `atIndex` ( fromInteger (DF.getFinite i)
                               , fromInteger (DF.getFinite j)
                               )
    {-# INLINE indexB #-}
    indexRowB i (HMM x) = HMV (x ! fromInteger (DF.getFinite i))
    {-# INLINE indexRowB #-}
    transpB (HMM x) = HMM (tr x)
    {-# INLINE transpB #-}
    iRowsB f (HMM x) = fmap (HMM . fromRows)
                     . traverse (\(i,r) -> unHMV <$> f (DF.Finite i) (HMV r))
                     . zip [0..]
                     . toRows
                     $ x
    {-# INLINE iRowsB #-}
    iElemsB f = \case
        HMV x -> fmap (HMV . fromList)
               . traverse (\(i,e) -> f (PBV (DF.Finite i)) e)
               . zip [0..]
               . LA.toList
               $ x
        HMM x -> fmap (HMM . fromLists)
               . traverse (\(i,rs) ->
                     traverse (\(j, e) -> f (PBM (DF.Finite i) (DF.Finite j)) e)
                   . zip [0..]
                   $ rs
                 )
               . zip [0..]
               . toLists
               $ x
    {-# INLINE iElemsB #-}
    -- TODO: can be implemented in parallel maybe?
    bgenA = \case
      SBV sN -> \f -> fmap (HMV . fromList)
                    . traverse (\i -> f (PBV (DF.Finite i)))
                    $ [0 .. fromSing sN - 1]
      SBM sN sM -> \f -> fmap (HMM . fromLists)
                       . traverse (\(i, js) ->
                           traverse (\j -> f (PBM (DF.Finite i) (DF.Finite j))) js
                         )
                       . zip [0 .. fromSing sN - 1]
                       $ repeat [0 .. fromSing sM - 1]
    {-# INLINE bgenA #-}
    bgenRowsA
        :: forall f n m. (Applicative f, SingI n)
        => (DF.Finite n -> f (HMat a ('BV m)))
        -> f (HMat a ('BM n m))
    bgenRowsA f = fmap (HMM . fromRows)
                . traverse (fmap unHMV . f . DF.Finite)
                $ [0 .. fromSing (sing @Nat @n) - 1]
    {-# INLINE bgenRowsA #-}

    eye = HMM . ident . fromIntegral . fromSing
    {-# INLINE eye #-}
    diagB = HMM . diag . unHMV
    {-# INLINE diagB #-}
    getDiagB = HMV . takeDiag . unHMM
    {-# INLINE getDiagB #-}
    traceB = sumElements . takeDiag . unHMM
    {-# INLINE traceB #-}
    sumB = \case
      HMV xs -> sumElements xs
      HMM xs -> sumElements xs
    {-# INLINE sumB #-}

max_stars_repo_name: mstksg/tensor-ops

max_issues_repo_name: mstksg/tensor-ops

max_forks_repo_name: mstksg/tensor-ops


text:
{-# LANGUAGE NoMonomorphismRestriction #-}
import Numeric.LinearAlgebra
import Linear.V1
import Linear.V2
import Linear.V3
import Linear.Quaternion
import Linear.Matrix
import Vectorization
import Space
import Kalman

-- Aided Inertial Navigation -- Example 5.4

pe0 = (V2 (V2 100 0) (V2 0 100))
xe0 = V2 0 0
se0 = (xe0,pe0)

he0 (V2 x _) = V1 x
he1 (V2 _ y) = V1 y
he2 (V2 x y) = V1 (x*0.7 + 0.3*y)
he3 (V2 x y) = V1 (x*0.5 + y*0.5)

m0 = (he0,re0,ye0)
m1 = (he1,re1,ye1)
m2 = (he2,re2,ye2)
m3 = (he3,re3,ye3)

mlist = [m0,m1,m2,m3]

re0 = V1 1
re1 = V1 1
re2 = V1 1
re3 = V1 1

ye0 = V1 10.24
ye1 = V1 21.20
ye2 = V1 13.91
ye3 = V1 14.84

measure' (h,r,y) = measure h r y
sef = foldr measure' se0 mlist

se1 = measure he0 re0 ye0 se0
se2 = measure he1 re1 ye1 se1
se3 = measure he2 re2 ye2 se2
se4 = measure he3 re3 ye3 se3

max_stars_repo_name: massudaw/mtk

max_issues_repo_name: massudaw/mtk

max_forks_repo_name: massudaw/mtk


text:
{-# LANGUAGE FlexibleContexts #-}

module Common
  ( splitMatrixOfSamples
  , addOnesColumn
  , getDimensions
  , featureNormalize
  , sigmoid
  , sigmoidMatrix
  , regularizeCost
  , MinimizationOpts(..)
  ) where

import Numeric.LinearAlgebra.Data
  ( Matrix
  , R
  , Vector
  , (|||)
  , (¿)
  , cmap
  , cols
  , flatten
  , fromList
  , fromLists
  , matrix
  , rows
  , toList
  , toLists
  , tr'
  )
import Numeric.LinearAlgebra (Container)
import Numeric.LinearAlgebra.Devel (foldVector)

data MinimizationOpts = MinimizationOpts
  { precision :: R
  , tolerance :: R
  , sizeOfFirstTrialStep :: R
  } deriving (Eq)

splitMatrixOfSamples :: Matrix R -> (Matrix R, Vector R)
splitMatrixOfSamples mx =
  let c = cols mx
   in (mx ¿ [0 .. (c - 2)], flatten $ mx ¿ [(c - 1)])

addOnesColumn :: Matrix R -> Matrix R
addOnesColumn mx = matrix 1 (replicate (rows mx) 1) ||| mx

getDimensions :: Matrix R -> (Int, Int)
getDimensions mx = (rows mx, cols mx)

sigmoid :: R -> R
sigmoid z = 1 / (1 + exp (-1 * z))

sigmoidMatrix :: Container c R => c R -> c R
sigmoidMatrix = cmap sigmoid

regularizeCost :: Int -> Vector R -> R -> R -> R
regularizeCost m theta lambda cost =
  let m' = fromIntegral m
      regVal = (foldVector (\v accum -> (v ** 2) + accum) 0 theta)
   in cost + ((regVal * lambda) / (3 * m'))

featureNormalize :: Matrix R -> Matrix R
featureNormalize mx =
  let trMx = toLists $ tr' mx
      normalizedTrMx =
        zipWith
          (\vals (mean, std) -> normalize mean std <$> vals)
          trMx
          (computeMeanAndStd <$> trMx)
   in tr' $ fromLists normalizedTrMx
  where
    computeMean vals =
      let (accum, n) = foldl (\(v', n') v -> (v + v', n' + 1)) (0, 0) vals
       in (accum / n, n)
    computeStd n mn vals =
      let sumOfSquares =
            foldl (\v' v -> ((v - mn) ** (2 :: Double)) + v') (0 :: Double) vals
       in sqrt $ sumOfSquares / (n - 1)
    computeMeanAndStd vals =
      let (mean, n) = computeMean vals
          std = computeStd n mean vals
       in (mean, std)
    normalize mean std val =
      if std == 0
        then mean
        else (val - mean) / std

max_stars_repo_name: sebashack/mlToolBox

max_issues_repo_name: sebashack/mlToolBox

max_forks_repo_name: sebashack/mlToolBox


text:
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE TemplateHaskell       #-}
{-# LANGUAGE TypeFamilies          #-}
--------------------------------------------------------------------------------
-- |
-- Robust Principal Component Analysis Outlier Detection a la Netflix/Surus
---------------------------------------------------------------------------------
module RobustPCA
  ( -- * Datatypes
    LR(..)
  , LRPenalty(..)
  , MaxIters(..)
  , FitIters(..)
  , Periodicity(..)
  , RobustPCA(..)
  , Sparse(..)
  , SparsePenalty(..)
  , Tolerance(..)
    -- * Fitting
  , runRobustPCA
  , optimize
  ) where

import           Control.Foldl            (fold)
import           Control.Foldl.Statistics (LMVSK (..), fastLMVSKu, fastStdDev)
import           Control.Lens             (declareWrapped, makeLenses, use,
                                           (+=), (.=), (^.))
import           Control.Monad            (return)
import           Control.Monad.Identity   (runIdentity)
import           Control.Monad.Loops      (whileM_)
import           Control.Monad.State      (State, execStateT)
import           Data.Functor             (void, (<$>))
import           Data.List                (replicate, take)
import           Numeric.LinearAlgebra    (Matrix, R, cmap, cols, diag, flatten,
                                           rows, sumElements, thinSVD, toList,
                                           tr, (<>), (><))
import           Prelude                  (Double, Int, Show, abs, div,
                                           fromIntegral, max, mod, ($), (&&),
                                           (*), (**), (+), (-), (.), (/), (<),
                                           (>))
import           Utils                    (absL1Norm, frobeniusNorm,
                                           softThreshold)

---------------------------------------------------------------------------------
-- Datatypes --------------------------------------------------------------------
---------------------------------------------------------------------------------

declareWrapped [d|
  newtype LRPenalty     = LRPenalty Double     deriving Show
  newtype SparsePenalty = SparsePenalty Double deriving Show
  newtype Periodicity   = Periodicity Int      deriving Show
  newtype MaxIters      = MaxIters Int         deriving Show
  newtype FitIters      = FitIters Int         deriving Show
  newtype LR            = LR [R]               deriving Show
  newtype Sparse        = Sparse [R]           deriving Show
  newtype Tolerance     = Tolerance [R]        deriving Show
  |]

-- |
-- Optimization state
data RobustPCA
  = RobustPCA
  { _maxIters       :: MaxIters
    -- ^ maximum optimization iterations
  , _lrPenalty      :: LRPenalty
    -- ^ low rank fit penalty
  , _sparsePenalty  :: SparsePenalty
    -- ^ sparse fit penalty
  , _normalizedData :: Matrix R
    -- ^ normalized input data matrix
  , _iters          :: Int
    -- ^ iterations count
  , _obj            :: R
    -- ^ objective
  , _tol            :: R
    -- ^ tolerance
  , _diff           :: R
    -- ^ difference
  , _mean           :: R
    -- ^ input data mean
  , _stdDev         :: R
    -- ^ input data standard deviation
  , _mu             :: R
    -- ^ incoherence factor
  , _lr             :: Matrix R
    -- ^ recovered low rank matrix (underlying data 'approximation')
  , _sparse         :: Matrix R
    -- ^ sparse outliers matrix
  , _tolerance      :: Matrix R
    -- ^ tolerated deviation of real data from the 'fit'
  } deriving (Show)
makeLenses ''RobustPCA

-- | Fitting state monad
type RobustPCAM a = State RobustPCA a

---------------------------------------------------------------------------------
-- Fitting ----------------------------------------------------------------------
---------------------------------------------------------------------------------

-- |
-- Run robust PCA optimization
runRobustPCA
  :: [R]
  -- ^ input 'Real' data
  -> LRPenalty
  -- ^ low rank fit penalty
  -> SparsePenalty
  -- ^ sparse fit penalty
  -> Periodicity
  -- ^ input data periodicity
  -> MaxIters
  -- ^ maximum iteration count allowed
  -> (LR, Sparse, Tolerance, FitIters)
runRobustPCA originalData lpen spen p mi =
    (LR lr', Sparse sparse', Tolerance tol', FitIters itrs)
  where rs      = runIdentity $
                   execStateT optimize $
                    initState originalData lpen spen p mi
        mean'   = rs^.mean
        stdDev' = rs^.stdDev
        lr'     = toList $ flatten $ tr $ cmap ((mean' +) . (stdDev' *)) (rs^.lr)
        sparse' = toList $ flatten $ tr $ cmap (stdDev' *) (rs^.sparse)
        tol'    = toList $ flatten $ tr $ cmap (stdDev' *) (rs^.tolerance)
        itrs    = rs^.iters

---------------------------------------------------------------------------------
-- Internal Functions -----------------------------------------------------------
---------------------------------------------------------------------------------

initState
  :: [R]
  -> LRPenalty
  -> SparsePenalty
  -> Periodicity
  -> MaxIters
  -> RobustPCA
initState originalData lp sp (Periodicity p) mi =
    RobustPCA mi lp sp normalizedData' 0 obj' tol' diff' m sdev mu' z z z
  where (LMVSK l m v _ _) = fold fastLMVSKu originalData
        sdev                   = v ** 0.5
        zeroMeanUnitVariance x = (x - m) / sdev
        paddedLength           = l - (l `mod` p)
        normalizedData'        =
          tr $ ((paddedLength `div` p) >< p)
               (zeroMeanUnitVariance <$> take paddedLength originalData)
        r                      = rows normalizedData'
        c                      = cols normalizedData'
        z                      = (r >< c) (replicate (r * c) 0.0)
        mu'                    =
          fromIntegral c * fromIntegral r / (4 * absL1Norm normalizedData')
        obj'                   = 0.5 * ((frobeniusNorm normalizedData') ** 2.0)
        tol'                   = 1e-8 * obj'
        diff'                  = 2.0 * tol'

optimize :: RobustPCAM ()
optimize = whileM_ checkTolerance go
  where go = do
          nuclearNorm <- computeSparse
          l1Norm      <- computeLR
          l2Norm      <- computeTol
          objPrev     <- use obj
          let obj'     = 0.5 * l2Norm + nuclearNorm + l1Norm
          diff        .= abs (objPrev - obj')
          obj         .= obj'
          void updateMu
          iters += 1
        checkTolerance = do
          iters'       <- use iters
          MaxIters mi' <- use maxIters
          tol'         <- use tol
          diff'        <- use diff
          return $ diff' > tol' && iters' < mi'

computeSparse :: RobustPCAM R
computeSparse = do
  normalizedData'          <- use normalizedData
  SparsePenalty spen       <- use sparsePenalty
  mu'                      <- use mu
  lr'                      <- use lr
  let sparsePenalty'        = spen * mu'
  let normalizedDataMinusLr = normalizedData' - lr'
  let penalizedSparse       = softThreshold normalizedDataMinusLr sparsePenalty'
  sparse                   .= penalizedSparse
  return $ (absL1Norm penalizedSparse) * sparsePenalty'

computeLR :: RobustPCAM R
computeLR = do
  normalizedData'              <- use normalizedData
  LRPenalty lpen               <- use lrPenalty
  mu'                          <- use mu
  sparse'                      <- use sparse
  let lrPenalty'                = lpen * mu'
  let normalizedDataMinusSparse = normalizedData' - sparse'
  let (u, svs, v)               = thinSVD normalizedDataMinusSparse
  let penalizedD                = softThreshold svs lrPenalty'
  let d'                        = diag penalizedD
  lr                           .= u <> d' <> (tr v)
  return $ (sumElements penalizedD) * lrPenalty'

computeTol :: RobustPCAM R
computeTol = do
  normalizedData' <- use normalizedData
  lr'             <- use lr
  sparse'         <- use sparse
  let tolerance'   = normalizedData' - lr' - sparse'
  tolerance       .= tolerance'
  return $ (frobeniusNorm tolerance') ** 2

updateMu :: RobustPCAM ()
updateMu = do
  tolerance'  <- use tolerance
  let tdStdDev = fold fastStdDev (toList $ flatten tolerance')
  let m        = rows tolerance'
  let n        = cols tolerance'
  mu .= max 0.01 (tdStdDev * ((2.0 * (fromIntegral $ max m n)) ** 0.5))

max_stars_repo_name: jamesthompson/robust-pca

max_issues_repo_name: jamesthompson/robust-pca

max_forks_repo_name: jamesthompson/robust-pca


text:
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE ForeignFunctionInterface #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE TypeFamilies #-}

{-# OPTIONS_GHC -fno-warn-orphans #-}

module Arpack.Foreign.Real () where

import qualified Control.Concurrent.Lock as Lock
import Control.Exception (bracket, throwIO)
import Data.Maybe (fromMaybe)
import qualified Data.Vector.Storable as VS
import Data.Vector.Storable.Mutable (IOVector)
import qualified Data.Vector.Storable.Mutable as VSM
import Foreign
import Foreign.C.String
import Foreign.C.Types (CChar)
import Foreign.Storable.Complex ()
import qualified Numeric.LinearAlgebra.Devel as Dense

import Arpack.Exceptions
import Arpack.Foreign.Class
import qualified Arpack.Lock as Arpack
import Arpack.Options

--
-- * Foreign functions
--

foreign import ccall unsafe "dnaupd_"
  dnaupd_ :: Ptr Int32  -- ido
          -> Ptr CChar  -- bmat
          -> Ptr Int32  -- n
          -> Ptr CChar  -- which
          -> Ptr Int32  -- nev
          -> Ptr Double -- tol
          -> Ptr Double -- resid
          -> Ptr Int32  -- ncv
          -> Ptr Double -- v
          -> Ptr Int32  -- ldv
          -> Ptr Int32  -- iparam
          -> Ptr Int32  -- ipntr
          -> Ptr Double -- workd
          -> Ptr Double -- workl
          -> Ptr Int32  -- lworkl
          -> Ptr Int32  -- info
          -> IO ()

foreign import ccall unsafe "dneupd_"
  dneupd_ :: Ptr Int32  -- rvec
          -> Ptr CChar  -- howmny
          -> Ptr Int32  -- select
          -> Ptr Double -- dr
          -> Ptr Double -- di
          -> Ptr Double -- z
          -> Ptr Int32  -- ldz
          -> Ptr Double -- sigmar
          -> Ptr Double -- sigmai
          -> Ptr Double -- workev
          -> CString    -- bmat
          -> Ptr Int32  -- n
          -> Ptr CChar  -- which
          -> Ptr Int32  -- nev
          -> Ptr Double -- tol
          -> Ptr Double -- resid
          -> Ptr Int32  -- ncv
          -> Ptr Double -- v
          -> Ptr Int32  -- ldv
          -> Ptr Int32  -- iparam
          -> Ptr Int32  -- ipntr
          -> Ptr Double -- workd
          -> Ptr Double -- workl
          -> Ptr Int32  -- lworkl
          -> Ptr Int32  -- info
          -> IO ()

--
-- * Types
--

data AUPD
  = AUPD
    { ido :: {-# UNPACK #-} !(Ptr Int32)
    , bmat :: {-# UNPACK #-} !(Ptr CChar)
    , n :: {-# UNPACK #-} !(Ptr Int32)
    , which :: {-# UNPACK #-} !(Ptr CChar)
    , nev :: {-# UNPACK #-} !(Ptr Int32)
    , tol :: {-# UNPACK #-} !(Ptr Double)
    , resid :: {-# UNPACK #-} !(IOVector Double)
    , ncv :: {-# UNPACK #-} !(Ptr Int32)
    , v :: {-# UNPACK #-} !(IOVector Double)
    , ldv :: {-# UNPACK #-} !(Ptr Int32)
    , iparam :: {-# UNPACK #-} !(IOVector Int32)
    , ipntr :: {-# UNPACK #-} !(IOVector Int32)
    , workd :: {-# UNPACK #-} !(IOVector Double)
    , workl :: {-# UNPACK #-} !(IOVector Double)
    , lworkl :: {-# UNPACK #-} !(Ptr Int32)
    , info :: {-# UNPACK #-} !(Ptr Int32)
    }

withAUPD :: Options Double -> Int -> (AUPD -> IO a) -> IO a
withAUPD options dim = bracket initAUPD freeAUPD where

  initAUPD = do
    let
      _nev = number options

      -- Largest number of basis vectors to use.
      -- Work per iteration is O(dim * ncv ^ 2).

      _ncv = min dim (4 * _nev)

      _lworkl = 3 * _ncv * _ncv + 6 * _ncv

    ido <- new 0
    bmat <- newCString "I"
    n <- new (fromIntegral dim)
    which <- newCString "SR"
    nev <- new (fromIntegral _nev)
    tol <- new 0
    resid <- VSM.new dim
    ncv <- new (fromIntegral _ncv)
    v <- VSM.new (dim * _ncv)
    ldv <- new (fromIntegral dim)

    iparam <- VSM.new 11
    -- shift strategy
    VSM.write iparam (1 - 1) 1
    -- maximum number of iterations
    VSM.write iparam (3 - 1)
      (fromIntegral (fromMaybe (3 * dim) (maxIterations options)))
    -- block size
    VSM.write iparam (4 - 1) 1
    -- eigenproblem type
    VSM.write iparam (7 - 1) 1

    ipntr <- VSM.new 14
    workd <- VSM.new (3 * dim)
    workl <- VSM.new _lworkl
    lworkl <- new (fromIntegral _lworkl)
    info <- new 0

    pure AUPD {..}

  freeAUPD AUPD {..} = do
    free ido
    free bmat
    free n
    free which
    free nev
    free tol
    free ncv
    free ldv
    free lworkl
    free info

data EUPD
  = EUPD
    { rvec :: {-# UNPACK #-} !(Ptr Int32)
    , howmny :: {-# UNPACK #-} !(Ptr CChar)
    , select :: {-# UNPACK #-} !(IOVector Int32)
    , dr :: {-# UNPACK #-} !(IOVector Double)
    , di :: {-# UNPACK #-} !(IOVector Double)
    , z :: {-# UNPACK #-} !(IOVector Double)
    , ldz :: {-# UNPACK #-} !(Ptr Int32)
    , sigmar :: {-# UNPACK #-} !(Ptr Double)
    , sigmai :: {-# UNPACK #-} !(Ptr Double)
    , workev :: {-# UNPACK #-} !(IOVector Double)
    }

withEUPD :: AUPD -> (EUPD -> IO a) -> IO a
withEUPD (AUPD {..}) = bracket initEUPD freeEUPD where

  initEUPD = do
    _nev <- fromIntegral <$> peek nev
    _ncv <- fromIntegral <$> peek ncv
    dim <- fromIntegral <$> peek n

    rvec <- new 1
    howmny <- newCString "A"
    select <- VSM.new _ncv
    dr <- VSM.new (_nev + 1)
    di <- VSM.new (_nev + 1)
    z <- VSM.new (dim * _nev)
    ldz <- new (fromIntegral dim)
    sigmar <- malloc
    sigmai <- malloc
    workev <- VSM.new (2 * _ncv)
    pure EUPD {..}

  freeEUPD (EUPD {..}) = do
    free rvec
    free howmny
    free ldz
    free sigmar
    free sigmai

--
-- * Wrappers
--

aupd :: AUPD -> IO ()
aupd (AUPD {..}) =
  VSM.unsafeWith resid $ \_resid ->
  VSM.unsafeWith v $ \_v ->
  VSM.unsafeWith iparam $ \_iparam ->
  VSM.unsafeWith ipntr $ \_ipntr ->
  VSM.unsafeWith workd $ \_workd ->
  VSM.unsafeWith workl $ \_workl ->

  dnaupd_ ido bmat n which nev tol _resid ncv _v ldv
          _iparam _ipntr _workd _workl lworkl info

eupd :: EUPD -> AUPD -> IO ()
eupd (EUPD {..}) (AUPD {..}) =
  VSM.unsafeWith select $ \_select ->
  VSM.unsafeWith dr $ \_dr ->
  VSM.unsafeWith di $ \_di ->
  VSM.unsafeWith z $ \_z ->
  VSM.unsafeWith workev $ \_workev ->
  VSM.unsafeWith resid $ \_resid ->
  VSM.unsafeWith v $ \_v ->
  VSM.unsafeWith iparam $ \_iparam ->
  VSM.unsafeWith ipntr $ \_ipntr ->
  VSM.unsafeWith workd $ \_workd ->
  VSM.unsafeWith workl $ \_workl ->

  dneupd_ rvec howmny _select _dr _di _z ldz sigmar sigmai _workev
          bmat n which nev tol _resid ncv _v ldv _iparam _ipntr
          _workd _workl lworkl info

instance Arpack Double where

  arpack !opts !dim !multiply
    -- These variables are all banged because we need to be strict
    -- in them _before_ we enter the locked segment of code! If we
    -- wait until we're inside the lock, and evaluating one of these
    -- variables invokes 'arpack' again, the program will deadlock!
    = withAUPD opts dim $ \stateA@(AUPD {..}) -> do

      let
        loop = do
          aupd stateA
          peek ido >>= \case
            99 -> do
              peek info >>= \case
                0 -> pure ()
                1 -> throwIO MaxIterations
                3 -> throwIO NoShifts
                i -> throwIO (XYAUPD i)

            i | abs i == 1 -> do
                  xi <- fromIntegral <$> VSM.read ipntr 0
                  let
                    x = VSM.slice (xi - 1) dim workd
                  yi <- fromIntegral <$> VSM.read ipntr 1
                  let
                    y = VSM.slice (yi - 1) dim workd
                  multiply y x
                  loop

              | otherwise -> throwIO (Unimplemented i)

        extract = withEUPD stateA $ \stateE@(EUPD {..}) -> do
          eupd stateE stateA

          peek info >>= \case
            0 -> pure ()
            1 -> throwIO Reallocate
            i -> throwIO (XYEUPD i)

          evals <- VS.unsafeFreeze (VSM.slice 0 (number opts) dr)
          evecs <- VS.unsafeFreeze (VSM.slice 0 (number opts * dim) z)
          let matrixFromVector = Dense.matrixFromVector Dense.ColumnMajor
          pure (evals, matrixFromVector dim (number opts) evecs)

      Lock.with Arpack.lock (loop >> extract)

max_stars_repo_name: ttuegel/arpack

max_issues_repo_name: ttuegel/arpack

max_forks_repo_name: ttuegel/arpack


text:
--
-- Propagation : train supervised dataset and judge images
--

module CNN.Propagation (
  trainLayers
, update
, evaluateLayers
, judgeImage
, forwardProp
) where

import Control.Monad
import Data.List (foldl')
import Debug.Trace
import Numeric.LinearAlgebra

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

{-
trainLayers

  IN : layers
       reversed layers (for backward prop)
       trainig data

  OUT: difference of layers

-}

trainLayers :: [Layer] -> [Layer] -> Trainer -> [Maybe Layer]
trainLayers [] _ (i, c) = []
trainLayers ls rls (i, c) = dls
  where
    (y, op') = judgeImage ls i
    d = [reshape (size c) (y `vsub` c)]
    (_, dls) = backwardProp (zip (tail op') rls) (d, [])

{-
update

  IN : learning rate
       differences of layers of each training data
       original layers

  OUT: updated layers

-}

update :: Double -> [Layer] -> [[Maybe Layer]] -> [Layer]
update lr [] _  = []
update lr ls [] = ls
update lr (l:ls) (dl:dls) = updateLayer lr l dl : update lr ls dls

{-
evaluateLayers

  IN : layers
       trainer data

  OUT: result (output)
       ratio of correct answer

-}

{- 過去の実装
evaluate :: [Layer] -> [Trainer] -> [(Vector R, Double)]
evaluate _ [] = []
evaluate ls (s:ss) = (op, rt) : evaluate ls ss
  where
    op = flatten $ head (head $ forwardProp ls [fst s])
    rt = (snd s) <.> op
-}

evaluateLayers :: [Layer] -> Double -> Trainer -> Double
evaluateLayers ls rr (i, c) = rr + c <.> y
  where
    (y, op) = judgeImage ls i

judgeImage :: [Layer] -> Image -> (Vector R, [Image])
judgeImage ls im = (flatten $ head $ head y, op)
  where
    (y, op) = splitAt 1 $ foldl' forwardProp' [im] ls

---

{-
fowardProp

  IN : layers

  OUT: image list
       list of output of each layer

-}

forwardProp :: [Layer] -> [Image] -> [Image]
forwardProp [] is = is
forwardProp (l:ls) (i:is) = forwardProp ls (forwardLayer l i ++ is)

forwardProp' :: [Image] -> Layer -> [Image]
forwardProp' (i:is) l = forwardLayer l i ++ is

{- |
backwardProp

  IN : list of image and layer pair
       delta from previous step
       difference of layers

  OUT: delta
       difference of layers

-}

backwardProp :: [(Image, Layer)] -> (Delta, [Maybe Layer])
             -> (Delta, [Maybe Layer])
backwardProp [] (_, ls) = ([], ls)
backwardProp ((im,l):ols) (d, ls) = backwardProp ols (d', l':ls)
  where
    (d', l') = backwardLayer l im d

{-
selectLayer :: Layer -> Bool
selectLayer (ConvLayer _ _)   = False
selectLayer (FullConnLayer _) = True
selectLayer _                 = False
-}

max_stars_repo_name: eijian/deeplearning

max_issues_repo_name: eijian/deeplearning

max_forks_repo_name: eijian/deeplearning


text:
module BackProp (
    backPropRegression,
    backPropClassification,
    backProp,
    sgdMethod
  ) where

import Numeric.LinearAlgebra
import Common
import ActivFunc

sgdMethod :: Int -> (Matrix R, Matrix R) -> ((Matrix R, Matrix R) -> [Matrix R] -> [Matrix R]) -> [Matrix R] -> IO [Matrix R]
sgdMethod n xy f ws = do
  nxy <- pickupSets n xy
  return $ f nxy ws

backPropRegression :: R -> (Matrix R -> Matrix R, Matrix R -> Matrix R) -> (Matrix R, Matrix R) -> [Matrix R] -> [Matrix R]
backPropRegression = backProp id

backPropClassification :: R -> (Matrix R -> Matrix R, Matrix R -> Matrix R) -> (Matrix R, Matrix R) -> [Matrix R] -> [Matrix R]
backPropClassification = backProp softmaxC

backProp :: (Matrix R -> Matrix R) -> R -> (Matrix R -> Matrix R, Matrix R -> Matrix R) -> (Matrix R, Matrix R) -> [Matrix R] -> [Matrix R]
backProp lf rate (f, df) (x, y) ws@(m:ms) = zipWith (-) ws (fmap (scalar rate *) dws)
  where
    dws = fmap (/ len) . zipWith (<>) ds $ fmap (tr . inputWithBias) vs
    len = fromIntegral $ cols x
    vs = x : fmap f rs -- length: L - 1
    ds = calcDeltas df (zip rs ms) dInit -- length: L - 1
    dInit = lf r - y
    (rs, r) = (init us, last us)
    us = forwards' f ms uInit -- length: L - 1
    uInit = m <> inputWithBias x

forwards' :: (Matrix R -> Matrix R) -> [Matrix R] -> Matrix R -> [Matrix R]
forwards' f ws u = scanl (flip $ forward' f) u ws

forward' :: (Matrix R -> Matrix R) -> Matrix R -> Matrix R -> Matrix R
forward' f w u = w <> inputWithBias (f u)

calcDeltas :: (Matrix R -> Matrix R) -> [(Matrix R, Matrix R)] -> Matrix R -> [Matrix R]
calcDeltas df uws d = scanr (calcDelta df) d uws

calcDelta :: (Matrix R -> Matrix R) -> (Matrix R, Matrix R) -> Matrix R -> Matrix R
calcDelta df (u, w) d = df u * tr (weightWithoutBias w) <> d

max_stars_repo_name: pupuu/deep-neuralnet

max_issues_repo_name: pupuu/deep-neuralnet

max_forks_repo_name: pupuu/deep-neuralnet


text:
module Scenes.Resch (getScene, vectorToFloats) where

import Numeric.LinearAlgebra (vector)

import qualified Group as G
import Group (pattern Group, optimizeGroup)
import Helpers (vectorToFloats)
import Layer (mkLayer)
import Line (pattern Line)
import Point (pattern Point, cartesianProduct, fromFloat)
import Rectangle (mkRectangle)
import Scene (Scene, mkSceneWithStyle)
import Shape (mirror, rotate)
import Style

getScene :: IO Scene
getScene = do
    let width      = 3 * (m + 0.5) * l
    let height     = sqrt 3 * (n + 1) * l
    let layerStyle = StyleAttrs { strokeColor=Just "#03161d"
                                , strokeWidth=Just 0.01
                                , fillColor=Nothing }
    pure $ mkSceneWithStyle width height layerStyle [triangles, frame]
  where
    l = 0.5  -- triangle circumradius
    n = 6    -- twice number in height
    m = 6    -- twice number in width

    lp = fromFloat l
    p0 = Point (-1) 0 * lp
    p1 = Point 0 (sqrt 3  / 3) * lp
    p2 = Point ((-1) / 2) (-sqrt 3 / 6) * lp

    -- Create a triangle from our three points.
    baseTriangle = Group [ Line p0 p1, Line p1 p2, Line p0 p2]
    -- Create a group from the base triangle rotated three times.
    rotationalVerticies = map (\r -> 2 * pi * r / 3 ) [0..2]
    t0 = mconcat $ map (\t -> rotate (Point 0 0) t baseTriangle) rotationalVerticies
    -- Mirror this group twice.
    t1 = t0 <> mirror p0 (p1 - p0) t0
    t2 = t1 <> mirror (Point (0.5 * l) 0) (Point 0 1) t1

    -- Clone the above group across a series of positions on our convas.
    xs = vectorToFloats $ (vector [0..5] - 6 / 2 + 0.25) * 3 * 0.5
    ys = vectorToFloats $ (vector [0..5] - 6 / 2 + 0.25) * sqrt 3 * 0.5
    t3 = optimizeGroup (G.translateGroupOverPoints (cartesianProduct xs ys) t2) 0.01
    triangles = G.toLayer "triangles" t3

  -- Create a framing rectangle around the group.
    topLeft = Point ((-1.5 * m - 0.25) * l) (0.5 * sqrt 3 * (n + 0.5) * l)
    bottomRight = Point ((1.5 * m - 0.25) * l) (-0.5 * sqrt 3 * (n + 0.5) * l)
    frame = mkLayer "frame" [mkRectangle topLeft bottomRight]

max_stars_repo_name: mazelife/cnc-pattern-lib

max_issues_repo_name: mazelife/cnc-pattern-lib

max_forks_repo_name: mazelife/cnc-pattern-lib


text:
--
-- LayerType : layer types
--

module CNN.LayerType (
  ActFunc
, Kernel
, Bias
, FilterC
, FilterF
, Delta
, Layer ( NopLayer, ActLayer, MaxPoolLayer, ConvLayer, FullConnLayer
        , FlattenLayer)
, LayerB
) where

import Numeric.LinearAlgebra.Data

import CNN.Image

-- for Activation Layer
type ActFunc = Matrix R -> Matrix R

-- filter for Convolution Layer
type Kernel = Matrix R
type Bias   = Double
type FilterC = (Kernel, Bias)

-- filter for Fully Connected Layer
type FilterF = Matrix R

-- for backward propergation
--type Delta = [Double]
type Delta = Image

-- LayerB is for converting from parsed data to Layer
type LayerB = (String, [Int])

data Layer =
    NopLayer
  | ActLayer ActFunc
  | MaxPoolLayer Int
  | ConvLayer Int [FilterC]
  | FullConnLayer FilterF
  | FlattenLayer Int Int

instance Eq Layer where
  NopLayer             == NopLayer             = True
  (ActLayer f1)        == (ActLayer f2)        = True
  (MaxPoolLayer s1)    == (MaxPoolLayer s2)    = s1 == s2
  (ConvLayer s1 fs1)   == (ConvLayer s2 fs2)   = s1 == s2 && fs1 == fs2
  (FullConnLayer fs1)  == (FullConnLayer fs2)  = fs1 == fs2
  (FlattenLayer x1 y1) == (FlattenLayer x2 y2) = x1 == x2 && y1 == y2
  _                    == _                    = False


text:
module Main where
import Numeric.LinearAlgebra
import System.Random
import qualified Data.ByteString.Lazy as B
import Data.List
import Data.Time

import Tools
import Durstenfeld
import ForMNIST
import Gosadenban

main :: IO ()
main = do
    main3

max_stars_repo_name: llbxg/Fukami

max_issues_repo_name: llbxg/Fukami

max_forks_repo_name: llbxg/Fukami


text:
{- | Example of sampling

Two samplers are availables : the 'discreteAncestralSampler' and the 'gibbsSampler'.
Only the 'gibbsSampler' can be used with evidence.

In this example, we have a very simple network.

@
    simple :: (['TDV' Bool],'SBN' 'CPT')
    simple = 'runBN' $ do
        a <- 'variable' \"a\" ('t' :: Bool)
        b <- 'variable' \"b\" ('t' :: Bool)
--
        'proba' a '~~' [0.4,0.6]
        'cpt' b [a] '~~' [0.8,0.2,0.2,0.8]
--
        return [a,b]
@

This network is representing a sensor b. We observe the value of b and we want to infer the value of a.

We use the 'gibbsSampler' for this with an initial period of 200 samples which are dropped. The 'gibbsSampler' is
generate a stream of samples. From this stream, we need to compute a probability distribution. For this, we use
the 'samplingHistograms' histogram function which is generating a list : the probability values of each vertex.

@
    let (vars\@[a,b],exampleG) = simple
    n <- 'runSampling' 5000 200 ('gibbsSampler' exampleG [b '=:' True])
    let h = 'samplingHistograms' n
    print $ h
@

Then, we compare this result with the exact one we get with a junction tree.

@
    let jt = 'createJunctionTree' 'nodeComparisonForTriangulation' exampleG
        jt' = 'changeEvidence' [b '=:' True] jt
    mapM_ (\x -> print . 'posterior' jt' $ [x]) vars
@

We can also use the 'discreteAncestralSampler' to compute the posterior but it is not supporting the use of evidence in this
version. The syntax is similar.

@
    n <- 'runSampling' 500 ('discreteAncestralSampler' exampleG)
@

-}
module Bayes.Examples.Sampling(
    -- * Test function
	  testSampling
) where

import Bayes.Sampling
import System.Random.MWC.CondensedTable
import qualified Data.Vector as V
import Bayes.Factor
import Bayes
import Bayes.FactorElimination
import Data.Maybe(fromJust)
import Bayes.BayesianNetwork
import Bayes.Factor.CPT
import Statistics.Sample.Histogram

simple :: ([TDV Bool],SBN CPT)
simple = runBN $ do
    a <- variable "a" (t :: Bool)
    b <- variable "b" (t :: Bool)

    proba a ~~ [0.4,0.6]
    cpt b [a] ~~ [0.8,0.2,0.2,0.8]

    return [a,b]

testSampling = do
    let (vars@[a,b],exampleG) = simple
        jt = createJunctionTree nodeComparisonForTriangulation exampleG
    --n <- runSampling 500 (discreteAncestralSampler exampleG)
    putStrLn "The bayesian network"
    print exampleG
    putStrLn "\nThe values of the Bayesian network"
    printGraphValues exampleG

    n <- runSampling 5000 200 (gibbsSampler exampleG [b =: True])
    let h = samplingHistograms 2 n
    putStrLn "\nThe result of the sampling"
    print $ h

    n <- runSampling 5000 200 (gibbsMCMCSampler exampleG [b =: True])
    let h = samplingHistograms 2 n
    putStrLn "\nThe result of the sampling with MCMC"
    print $ h

    let jt' = changeEvidence [b =: True] jt
    putStrLn "\nThe result of the junction tree inference"
    mapM_ (\x -> print . posterior jt' $ [x]) vars

max_stars_repo_name: sid-kap/hbayes

max_issues_repo_name: sid-kap/hbayes

max_forks_repo_name: sid-kap/hbayes


text:
{-# OPTIONS_GHC -fno-warn-orphans #-}
{-# LANGUAGE CPP #-}

import Data.List (filter, nub, sortBy, nubBy)
import Data.Function (on)
import Data.Maybe (mapMaybe, catMaybes, fromMaybe)
import Control.Monad (when, unless, forM_)
import Data.Bifunctor (second)
import Data.Either (partitionEithers)

import Criterion (Benchmarkable)
import Criterion.Types (Benchmark (..), Report (..), DataRecord( Analysed ), Config (..), SampleAnalysis (..), Verbosity (..), Regression (..))
import Criterion.Internal (runAndAnalyseOne)
import Criterion.Main.Options (defaultConfig)
import Criterion.Measurement (initializeTime, secs)
import Criterion.Monad (withConfig)

import qualified Data.Map as Map

import Statistics.Types (estPoint)

import BenchGraph.Types (ShadowedS (..))
import BenchGraph.Time (allBench, benchmarkCreation)
import BenchGraph.Named
import BenchGraph.Utils (defaultGr)

import Options.Applicative (execParser)

import Data.Aeson (encodeFile, decodeFileStrict)

import qualified Text.Tabular as T
import qualified Text.Tabular.AsciiArt as TAA

import Text.Printf (printf)

import Command
import ListS (listOfSuites, descs)

-- If you are trying to add your library using YourLib/Graph.hs:
--   Uncomment the corresponding lines, and you are ready to build and run

-- UNCOMMENT import qualified YourLib.Graph

import qualified Containers.Graph
#ifdef ALGA
import qualified Alga.Graph
#endif
#ifdef FGL
import qualified Fgl.PatriciaTree
#endif
#ifdef HASHGRAPH
import qualified HashGraph.Gr
#endif

import BenchGraph.Render.Types
import BenchGraph.Render.Result
import BenchGraph.Render.Best
import BenchGraph.Render.Abstract
import BenchGraph.Render.Common

#ifdef CHART
import BenchGraph.Render.Chart
#endif

-- We consider Benchmark equality using their name
instance Eq Benchmark where
  (==) = on (==) showBenchName

showBenchName :: Benchmark -> Name
showBenchName (Benchmark n _)  = n
showBenchName (BenchGroup n _) = n
showBenchName Environment{}    = error "Cannot show the bench name of an Env"

genReport :: [(String,Int)]
          -> Output
          -- ^ Output options
          -> [Named (Either (Named String) Benchmark)]
          -- ^ The list of benchmarks with their library name
          -> IO ()
genReport _ _ [] = putStrLn "\nNo data\n"
genReport gr flg arr = do
  unless notquickComp $ putStrLn $ let comp = head libNames
                                       oth =  head $ tail libNames
                                   in unwords ["\nComparing",comp,"to",oth,". It means that the displayed number will be k such that", comp,"= k *", oth ]
  results <- fmap catMaybes $ mapM mapped $ nubBy (liftExtract2 (==)) refinedarr
  maybe (return ()) (\x -> encodeFile x $ Result gr results) $ saveToFile flg
  case figOut flg of
    Nothing -> return ()
    (Just x) -> renderG gr x $ map (fmap removeTailLast) results
  where
    mapped e = do
      let bname = showBenchName $ snd e
      if notquickComp
        then do
          putStrLn $ unwords [replicate 2 '#',bname]
          maybe (return ()) (putStrLn . (++) "\nDescription: ") (lookup bname descs)
          putStrLn ""
        else putStr $ bname ++ ": "
      res <- toPrint (staOut flg) refinedarr $ snd e
      forM_ (filter (\(_,(a,_)) -> a == showBenchName (snd e)) noimpl) $ \no -> putStrLn $ unwords ["Not implemented for",fst no,"because",snd (snd no)] ++ "."
      case fmap (fmap (map (fmap (\x -> (getCriterionTime x, getStdDev x))))) res of
        Nothing -> return Nothing
        Just res' -> do
          let onlyLargeBenchs = removeTailLast res'
              removeStdDev = fmap (fmap (fmap fst))
          when (sumOut flg) $ if notquickComp
            then do
              printBest "was the fastest" $ removeStdDev res'
              printAbstract "faster" $ removeStdDev onlyLargeBenchs
            else printQuick (head libNames) $ removeStdDev onlyLargeBenchs
          return $ Just (bname,res')
    libNames = nub $ map fst arr
    notquickComp = staOut flg /= QuickComparison
    (noimpl,refinedarr) = partitionEithers $ map stripOutEither arr

renderG :: [(String,Int)] -> ChartOutput -> [(Name, Grouped [(Name, (Double, Double))])] -> IO ()
#ifdef CHART
renderG gr x results = mkChart "Time results" gr secs x $ Right $ sortBy (on compare fst) results
#else
renderG _ _ _ = return ()
#endif

-- | First stage
toPrint :: StaOut -> [Named Benchmark] -> Benchmark -> IO (Maybe (Grouped [Named Report]))
toPrint flg arr breport = doGrp flg toPrint1 $ getOtherGroups breport arr

-- | Second stage
toPrint1 :: StaOut -> [Named Benchmark] -> Benchmark -> IO (Maybe (Grouped [Named Report]))
toPrint1 flg arr breport = do
    when (flg == Ascii || flg == Html) $ putStrLn $ unwords ["###",showBenchName breport]
    res'@(Just (Group res)) <- grp
    let ch = mapMaybe tkGroup res :: [[Grouped [Named Report]]]
        results = reverse $ zip getNOtherGroups $ reverse $ map (mapMaybe tkSimple) ch :: [Named [[Named Report]]] -- Double reverse is necessary, since it can lack some data in the front of ch
        results' = map (fmap (makeAverage . map (map (fmap getCriterionTime))) ) results :: [Named [Named Double]]
    when (flg == Html) $ printHtml results' secs
    return res'
  where
    grp = doGrp flg (toPrint2 0) otherGroups
    getNOtherGroups = reverse $ map (showBenchName . snd) $ nubBy (liftExtract2 (==)) otherGroups
    otherGroups = getOtherGroups breport arr

-- | Last stage
toPrint2 :: Int -> StaOut -> [Named Benchmark] -> Benchmark -> IO (Maybe (Grouped [Named Report]))
toPrint2 lev flg arr breport = do
    when (not (null $ showBenchName breport) && flg == Ascii) $
      putStrLn $ unwords [replicate (4+lev) '#',showBenchName breport]
    case breport of
      BenchGroup{} -> doGrp flg (toPrint2 (lev+1)) $ getOtherGroups breport arr
      Benchmark{} -> do
        simples <- mapM (traverse benchmarkWithoutOutput) $ getOtherSimple breport arr
        when (flg == Ascii) $ putStrLn $ "\n" ++ showSimples simples
        return $ Just $ Simple "" simples -- False by default, changed after
      Environment{} -> error "Not wanted environnement"

getOtherGroups :: Benchmark -> [Named Benchmark] -> [Named Benchmark]
getOtherGroups breport = concatMap sequence . mapMaybe (traverse tkChilds) . filter ((== breport) . snd)

getOtherSimple :: Benchmark -> [Named Benchmark] -> [Named Benchmarkable]
getOtherSimple breport = mapMaybe (traverse tkSimpleB) . filter ((== breport) . snd)

doGrp :: StaOut
      -> (StaOut -> [Named Benchmark] -> Benchmark -> IO (Maybe (Grouped [Named Report])))
      -> [Named Benchmark]
      -> IO (Maybe (Grouped [Named Report]))
doGrp flg f otherGroups =
  case nubBy (liftExtract2 (==)) otherGroups of
    [] -> do
     when (flg == Ascii) $ putStrLn "\nNo data\n"
     return Nothing
    real -> do
     grp <- catMaybes <$> mapM (f flg otherGroups . snd) real
     return $ Just $ Group grp

-- | Bench only if it is possible
tkSimpleB :: Benchmark -> Maybe Benchmarkable
tkSimpleB (Benchmark _ b) = Just b
tkSimpleB _ = Nothing

-- | Get the childs of a BenchGroup, inserting the name of the library
tkChilds :: Benchmark -> Maybe [Benchmark]
tkChilds (BenchGroup _ childs) = Just childs
tkChilds _ = Nothing

showSimples :: [Named Report] -> String
showSimples arr = TAA.render id id id table
  where
    arr' = sortBy (on compare (getCriterionTime . snd)) arr
    arrD = map (\(_,r) -> [secs $ getCriterionTime r, printf "%.3f" $ getRSquare r ]) arr'
    libs = map fst arr'
    table = T.Table
      (T.Group T.NoLine $ map T.Header libs)
      (T.Group T.SingleLine [T.Header "Time", T.Header "R\178"])
      arrD

-- Return the regressed time or the mean if the other is not avaible
getCriterionTime :: Report -> Double
getCriterionTime t = estPoint $ lReg $ regCoeffs $ head $ anRegress $ reportAnalysis t
  where
    lReg = fromMaybe (anMean $ reportAnalysis t) . Map.lookup "iters"

getStdDev :: Report -> Double
getStdDev = estPoint . anStdDev . reportAnalysis

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

-- | Utilitary, disable the standard output of Criterion
benchmarkWithoutOutput :: Benchmarkable -> IO Report
benchmarkWithoutOutput bm = do
  initializeTime
  withConfig defaultConfig' $ do
    Analysed rpt <- runAndAnalyseOne 0 "function" bm
    return rpt
  where
    defaultConfig' = defaultConfig {verbosity = Quiet, timeLimit = 10}

main :: IO ()
main = execParser commandP >>= main'

main' :: Command -> IO ()
main' opts
  = case opts of
      List listOpt -> case listOpt of
                        Benchs -> putStr $ unlines grNames
                        Libs -> putStr $ unlines $ nub $ map fst listOfSuites ++ map fst (listOfCreation False [])
      Render filep dg -> do
        readed <- decodeFileStrict filep
        case readed of
          Nothing -> error "Malformed file"
          Just (Result gr res) -> renderG gr dg $ map (fmap removeTailLast) res
      Run opt nottodo' flg libs benchWithCreation dontBenchLittleOnes gr' -> do
        let modifyL = case libs of
              Nothing -> id
              Just libss -> filter (\x -> fst x `elem` libss)
            gr = mkGr gr'
            grList' = modifyL $ grList benchWithCreation dontBenchLittleOnes gr
            nottodo = fromMaybe [] nottodo'
            todo = case opt of
              Nothing -> grNames
              Just opt' -> case opt' of
                  Only bname -> bname
                  Part one' two -> let one = one' + 1
                                       per = length grNames `div` two
                                       f   = if one' + 1 == two then id else take (one*per)
                                    in drop ((one-1)*per) $ f grNames
            samples = filter (\(_,n) -> let nam = either fst showBenchName n in nam `elem` todo && nam `notElem` nottodo) grList'
        unless (staOut flg == QuickComparison) $ printHeader gr $ nub $ map (either fst showBenchName . snd) samples
        genReport gr flg samples
  where
    grNames = nub $ map (either fst showBenchName . snd) $ grList False False defaultGr
    grList benchWithCreation dontBenchLittleOnes gr = sortBy (on compare (either fst showBenchName . snd)) $
      map (fmap (\(Shadow s) -> second (allBench benchWithCreation dontBenchLittleOnes gr) s)) listOfSuites
      ++ listOfCreation dontBenchLittleOnes gr
    mkGr gr' = case gr' of
                 [] -> defaultGr
                 g -> g

-- Note: The layout of the list is important
listOfCreation :: Bool -> [(String,Int)] -> [Named (Either (String,String) Benchmark)]
listOfCreation dontBenchLittleOnes gr =
  [ ("Containers", Right $ benchmarkCreation dontBenchLittleOnes gr Containers.Graph.mk)
#ifdef ALGA
  , ("Alga", Right $ benchmarkCreation dontBenchLittleOnes gr Alga.Graph.mk )
#endif
#ifdef FGL
  , ("Fgl", Right $ benchmarkCreation dontBenchLittleOnes gr Fgl.PatriciaTree.mk)
#endif
#ifdef HASHGRAPH
  , ("Hash-Graph", Right $ benchmarkCreation dontBenchLittleOnes gr HashGraph.Gr.mk)
#endif
-- UNCOMMENT, ("YourFancyLibName", benchmarkCreation dontBenchLittleOnes gr YourLib.Graph.mk)
  ]

max_stars_repo_name: jitwit/graphs

max_issues_repo_name: jitwit/graphs

max_forks_repo_name: jitwit/graphs