text:
module Parser (readExpr, readExprList) where

import Control.Monad.Except
import Data.Char (digitToInt)
import Data.Maybe (maybeToList)
import Data.Ratio
import Data.Complex
import Numeric (readHex, readOct, readInt, readFloat, readDec)
import Text.ParserCombinators.Parsec
import Text.ParserCombinators.Parsec.Expr
import Text.ParserCombinators.Parsec.Language
import qualified Text.ParserCombinators.Parsec.Token as T

import Definition
import LispVector (vector)

------- The Exposed Parsing Function -------

readOrThrow :: Parser a -> String -> Evaled a
readOrThrow parser input = case parse parser "lisp" input of
  Left err  -> throwError $ ParserErr err
  Right val -> return val

readExpr :: String -> Evaled LispVal
readExpr = readOrThrow $ parseExpr <* eof

readExprList :: String -> Evaled [LispVal]
readExprList = readOrThrow $ endBy parseExpr spaces

------- Parsers -------

parseExpr :: Parser LispVal
parseExpr = parseNumber
        <|> parseAtom
        <|> parseChar
        <|> parseBool
        <|> parseString
        <|> parseAnyQuoted
        <|> parseAnyList
        <|> parseVector

------- Any Atom Parser -------

parseAtom :: Parser LispVal
parseAtom = try $ do
  atom <- identifier
  if atom == "."
    then pzero
    else return $ LAtom atom

------- Bool Parser -------

parseBool :: Parser LispVal
parseBool = try $ do
  char '#'
  boolChar <- oneOf "tf"
  return . LBool $ case boolChar of
    't' -> True
    _   -> False

------- Character Literal Parser -------

parseChar :: Parser LispVal
parseChar = try (string "#\\") >> fmap LChar character
  where
    character :: Parser Char
    character = do
      ch   <- anyChar
      rest <- many letter
      case ch:rest of
        "altmode"   -> return '\ESC'
        "backnext"  -> return '\US'
        "backspace" -> return '\BS'
        "call"      -> return '\SUB'
        "delete"    -> return '\DEL'
        "escape"    -> return '\ESC'
        "linefeed"  -> return '\LF'
        "newline"   -> return '\n'
        "null"      -> return '\NUL'
        "page"      -> return '\FF'
        "return"    -> return '\CR'
        "rubout"    -> return '\DEL'
        "space"     -> return ' '
        "tab"       -> return '\t'
        [c]         -> return c
        _           -> pzero

------- String Parser -------

parseString :: Parser LispVal
parseString = do
  char '"'
  chars <- many $ validStringChar <|> noneOf "\""
  char '"'
  return $ LString chars

validStringChar :: Parser Char
validStringChar = do
  char '\\'
  c <- oneOf "\\\"abnrt"
  return $ case c of
    'a' -> '\a'
    'b' -> '\b'
    'n' -> '\n'
    'r' -> '\r'
    't' -> '\t'
    _   -> c

------- Number Parsers -------

parseNumber :: Parser LispVal
parseNumber = fmap LNumber . try $ parseSNumber <* notFollowedBy (letter <|> digit <|> symbol)

parseSNumber :: Parser SchemeNumber
parseSNumber = tryRational <|> tryComplex <|> noBase <|> try withBase
  where
    maybeNeg :: Num a => Parser (a -> a)
    maybeNeg = maybe id (const negate) <$> optionMaybe (char '-')

    toInt :: String -> Integer
    toInt = fst . head . readDec

    toDouble :: String -> Double
    toDouble = fst . head . readFloat

    getNumInt :: Parser Integer
    getNumInt = do
      mbNeg  <- maybeNeg
      digits <- many1 digit
      return . mbNeg . toInt $ digits

    getNumDouble :: Parser Double
    getNumDouble = do
      mbNeg  <- maybeNeg
      (a, b) <- many1 digit `sepByOnce` char '.'
      return . mbNeg . toDouble $ a ++ '.':b

    noBase :: Parser SchemeNumber
    noBase = fmap SDouble (try getNumDouble) <|> fmap SInt (try getNumInt)

    withBase :: Parser SchemeNumber
    withBase = char '#' >> oneOf baseChars >>= getNumFromBaseChar

    baseChars :: String
    baseChars = "bodx"

    getNumFromBaseChar :: Char -> Parser SchemeNumber
    getNumFromBaseChar 'b' = SInt . fst . head . readBin <$> many1 (oneOf binChars)
    getNumFromBaseChar 'o' = SInt . fst . head . readOct <$> many1 octDigit
    getNumFromBaseChar 'd' = noBase
    getNumFromBaseChar 'x' = SInt . fst . head . readHex <$> many1 hexDigit

    tryRational :: Parser SchemeNumber
    tryRational = try $ do
      (a, b) <- getNumInt `sepByOnce` char '/'
      return . SRational $ a % b

    tryComplex :: Parser SchemeNumber
    tryComplex = fmap SComplex $ rectComplex <|> polarComplex

    rectComplex :: Parser (Complex Double)
    rectComplex = try $ do
      rPart    <- getNum
      iPartNeg <- oneOf "+-" >>= \op -> return $ if op == '-'
                                          then negate
                                          else id

      iPart    <- iPartNeg <$> getNum
      char 'i'
      return $ rPart :+ iPart

    polarComplex :: Parser (Complex Double)
    polarComplex = try $ do
      mag   <- maybeNeg <*> getNumInt
      char '@'
      angle <- getNum
      return $ mkPolar (fromIntegral mag) angle

    getNum :: Parser Double
    getNum = try getNumDouble <|> fmap fromIntegral getNumInt

------- Quoted Parsers -------

parseAnyQuoted :: Parser LispVal
parseAnyQuoted =   parseQuasiQuoted
               <|> parseUnquoteSpliced
               <|> parseUnquoted
               <|> parseQuoted

parseQuoted :: Parser LispVal
parseQuoted = do
  lexeme $ char '\''
  x <- parseExpr
  return $ LList [LAtom "quote", x]

parseQuasiQuoted :: Parser LispVal
parseQuasiQuoted = do
  lexeme $ char '`'
  expr <- parseExpr
  return $ LList [LAtom "quasiquote", expr]

parseUnquoted :: Parser LispVal
parseUnquoted = do
  lexeme $ char ','
  expr <- parseExpr
  return $ LList [LAtom "unquote", expr]

parseUnquoteSpliced :: Parser LispVal
parseUnquoteSpliced = do
  try . lexeme $ string ",@"
  expr <- parseExpr
  return $ LList [LAtom "unquote-splicing", expr]

------- List Parsers -------

parseAnyList :: Parser LispVal
parseAnyList = parens anyList
  where
    anyList :: Parser LispVal
    anyList = try parseDottedList <|> parseList

parseList :: Parser LispVal
parseList = LList <$> parseExpr `sepEndBy` whiteSpace

parseDottedList :: Parser LispVal
parseDottedList = do
  vals <- parseExpr `endBy` whiteSpace
  val  <- dot >> parseExpr
  return $ case val of
    LList vs -> LList $ vals ++ vs
    _        -> LDottedList vals val

------- Vector Parser -------

parseVector :: Parser LispVal
parseVector = do
  lexeme $ string "#("
  vals <- parseExpr `sepEndBy` whiteSpace
  lexeme $ char ')'
  return . LVector $ vector vals

------- Helper Parsers -------

binChars :: String
binChars = "01"

-- http://stackoverflow.com/questions/5921573/convert-a-string-representing-a-binary-number-to-a-base-10-string-haskell
readBin :: String -> [(Integer, String)]
readBin = readInt 2 (`elem` binChars) digitToInt

symbolChars :: String
symbolChars = "!$%&|*+/:<=>?@^_~-."

symbol :: Parser Char
symbol = oneOf symbolChars

spaces1 :: Parser ()
spaces1 =  skipMany1 space

notSpace :: Parser Char
notSpace = satisfy (/=' ')

sepByOnce :: Parser a -> Parser sep -> Parser (a, a)
sepByOnce p sep = (,) <$> p <* sep <*> p

------- Language Definition -------

lispDef :: LanguageDef ()
lispDef = emptyDef {
  T.commentStart    = "#|"
, T.commentEnd      = "|#"
, T.commentLine     = ";"
, T.nestedComments  = True
, T.identStart      = letter <|> symbol
, T.identLetter     = letter <|> symbol <|> digit
, T.caseSensitive   = True
}

lexer :: T.TokenParser ()
lexer = T.makeTokenParser lispDef

identifier :: Parser String
identifier = T.identifier lexer

parens :: Parser a -> Parser a
parens = T.parens lexer

whiteSpace :: Parser ()
whiteSpace = T.whiteSpace lexer

dot :: Parser String
dot = T.dot lexer

lexeme :: Parser a -> Parser a
lexeme = T.lexeme lexer

max_stars_repo_name: comraq/scheme-interpreter

max_issues_repo_name: comraq/scheme-interpreter

max_forks_repo_name: comraq/scheme-interpreter


text:
-- file: Astro/PostNewtonian.hs
-- Post-Newtonian approximations for two-body problem

module Astro.PostNewtonian
	where

import Astro.Misc
import Numeric.LinearAlgebra

import Data.List (foldl')

-- Calculate Post-Newtonian acceleration functions for non-spin
-- dependent terms
mkNonSpinPNs :: Double->Double->Double->Double->
	[(Vector Double->Vector Double->Vector Double)]
mkNonSpinPNs gc c m1 m2 = 	[ a1
				, a2
				, a2_5
				]
			where
	-- Keplerian contribution
	--a0 :: Vector Double -> Vector Double -> Vector Double
	--a0 rv vv = scale (-gm / r^3) rv where r = norm rv

	a1 :: Vector Double -> Vector Double -> Vector Double
	a1 rv vv = {-# SCC "a1" #-}
		scale (-gmperc2 / r^2) $
			scale nbrac nv + scale vbrac vv
		where
			nbrac = -2*(2+eta)*gm/r + (1+3*eta)*v^2 - 3/2*eta*rd^2
			vbrac = -2*(2-eta)*rd
			r	= norm rv
			v	= norm vv
			nv	= unitV rv
			rd	= (rv <.> vv)/r

	a2 :: Vector Double -> Vector Double -> Vector Double
	a2 rv vv = {-# SCC "a2" #-}
		scale (- gmperc2 / (c^2 * r^2)) $
			scale (nbrac/r) rv + scale ((-rd/2) * vbrac) vv
		where
			nbrac = 3/4*(12+29*eta)*(gm/r)^2 + eta*(3-4*eta)*v^4
				+ 15/8*eta*(1-3*eta)*rd^4
				- 3/2*eta*(3-4*eta)*v^2*rd^2
				- 1/2*eta*(13-4*eta)*(gm/r)*v^2
				- (2+25*eta+2*eta^2)*(gm/r)*rd^2
			vbrac = eta*(15+4*eta)*v^2 - (eta+41*eta+8*eta^2)*(gm/r)
				- 3*eta*(3+2*eta)*rd^2
			r	= norm rv
			v	= norm vv
			rd	= (rv <.> vv)/r

	a2_5 :: Vector Double -> Vector Double -> Vector Double
	a2_5 rv vv = {-# SCC "a2_5" #-}
		scale (8/15*(gmperc2)^2*eta/(c*r^2)) $
			scale (nbrac*rd/r) rv + scale (-vbrac) vv
		where
			nbrac 	= 9*v^2 + 17*gm/r
			vbrac 	= 3*v^2 + 9*gm/r
			r	= norm rv
			v	= norm vv
			rd	= (rv <.> vv)/r
	gm 	= gc * m
	gmperc2 = gm / c ^2
	m 	= m1+m2
	eta 	= m1*m2/m^2

-- Calculate Post-Newtonian acceleration functions for spin
-- dependent terms
--mkSpinPNs :: Double->Double->Double->Double->Double->Double->
--	[(Vector Double->Vector Double->Vector Double->Vector Double)]
mkSpinPNs gc c m1 m2 chi q = 	[ spinOrbit
				, quadMonoInteraction
				] where
	--spinOrbit :: Vector Double -> Vector Double -> Vector Double
	--		-> Vector Double
	spinOrbit rv vv s = {-# SCC "SO" #-}
		scale (gm^2/(r^3*c^3)*(1+sqrt (1-4*eta))/4*chi) $
				scale nbrac nv + scale (nsbrac) (cross nv s')
					+ scale (-vsbrac) (cross vv s')
		where
			nbrac = 12 * s' <.> (cross nv vv)
			nsbrac = (9 + 3 * sqrt (1-4*eta))*rd
			vsbrac = 7 + sqrt(1-4*eta)
			r	= norm rv
			v	= norm vv
			nv	= unitV rv
			rd	= (rv <.> vv)/r
			s'	= spinCylToSpher s

	--quadMonoInteraction :: Vector Double -> Vector Double -> Vector Double
	--			-> Vector Double
	quadMonoInteraction rv vv s = {-# SCC "Q" #-}
		-- FIXME: Minus or no minus? Depends on whether
		-- coordinates fixed on m1 or m2, see OJ287 paper, minus
		-- in paper.
		scale (-q*chi^2*3*gc^3*m1^2*m/(2*c^4*r^4)) $
			scale nbrac nv + scale (-2 * nv <.> s') s'
		where
			nbrac = 5 * (nv <.> s')^2 - 1
			r	= norm rv
			v	= norm vv
			nv	= unitV rv
			s'	= spinCylToSpher s
	gm 	= gc * m
	m 	= m1+m2
	eta 	= m1*m2/m^2

-- Time derivative of spin unit vector, in spherical coordinates
spinDt :: Double->Double->Double->Double->Double->Double->
	Vector Double->Vector Double->Vector Double->Vector Double
spinDt gc c m1 m2 chi q rv vv s = cross omega s
	where
	omega = scale (gm*eta/(2*c^2*r^2)*(7+sqrt (1-4*eta))
				/(1+sqrt (1-4*eta))) (cross nv vv)
	r = norm rv
	nv = scale (1/r) rv
	gm 	= gc * m
	m 	= m1+m2
	eta 	= m1*m2/m^2

-- Time derivatives of spin unit vector components in cylindrical
-- coordinates, with x as the axisymmery axis
spinDtCyl :: Double->Double->Double->Double->Double->Double->
	Vector Double->Vector Double->[(Double,Double)->Double]
spinDtCyl gc c m1 m2 chi q rv vv = [dtheta, dxi]
	where
	dtheta (th,xi) =
		-xiPerRho*(dhs@>1*cos th + dhs@>2*sin th) + dhs@>0
		where
		xiPerRho = if abs xi >= 1 then 0 else xi/sqrt (1-xi^2)
	dxi (th,xi) =
		rho * (dhs@>1*sin th - dhs@>2 * cos th)
		where
		rho = if abs xi >= 1 then 0 else sqrt $ 1 - xi^2
	r = norm rv
	nv = scale (1/r) rv
	dhs = omega
	omega = scale (gm*eta/(2*c^2*r^2)*(7+sqrt (1-4*eta))
				/(1+sqrt (1-4*eta))) (cross nv vv)
	gm 	= gc * m
	m 	= m1+m2
	eta 	= m1*m2/m^2

-- Convert spin from cylindrical to cartesian representation
-- Axisymmetry axis is x-axis
spinCylToSpher :: (Double,Double)->Vector Double
spinCylToSpher (theta,xi) = --3|>[rho * cos theta, rho*sin theta, xi]
	3|>[xi, rho * cos theta, rho*sin theta]
	where rho = if abs xi >= 1 then 0 else sqrt $ 1 - xi^2

-- Post-Newtonian acceleration functions for highest supported degree
nonSpinPNAcc gc c m1 m2 rv vv =
	sum [f rv vv | f <- mkNonSpinPNs gc c m1 m2]
spinPNAcc gc c m1 m2 chi q rv vv s =
	sum [f rv vv s | f <- mkSpinPNs gc c m1 m2 chi q]

totalPNAcc gc c m1 m2 chi q rv vv s = ((f1 rv vv) + (f2 rv vv s))
	where (f1, f2) = (nonSpinPNAcc gc c m1 m2,
				spinPNAcc gc c m1 m2 chi q)

-- Explicitly calculated total acceleration
totalPNAccExp gc c m1 m2 chi q rv vv s =
	a1 + spinOrbit + quadMono + a2 + a2_5
	where
	a1 = {-# SCC "a1" #-}
		scale (-gmperc2 / r^2) $
			scale nbrac nv + scale vbrac vv
		where
			nbrac = -2*(2+eta)*gm/r + (1+3*eta)*v^2 - 3/2*eta*rd^2
			vbrac = -2*(2-eta)*rd

	a2 = {-# SCC "a2" #-}
		scale (- gmperc2 / (c^2 * r^2)) $
			scale nbrac nv + scale ((-rd/2) * vbrac) vv
		where
			nbrac = 3/4*(12+29*eta)*(gm/r)^2 + eta*(3-4*eta)*v^4
				+ 15/8*eta*(1-3*eta)*rd^4
				- 3/2*eta*(3-4*eta)*v^2*rd^2
				- 1/2*eta*(13-4*eta)*(gm/r)*v^2
				- (2+25*eta+2*eta^2)*(gm/r)*rd^2
			vbrac = eta*(15+4*eta)*v^2 - (eta+41*eta+8*eta^2)*(gm/r)
				- 3*eta*(3+2*eta)*rd^2

	a2_5 = {-# SCC "a2_5" #-}
		scale (8/15*(gmperc2)^2*eta/(c*r^2)) $
			scale (nbrac*rd) nv + scale (-vbrac) vv
		where
			nbrac 	= 9*v^2 + 17*gm/r
			vbrac 	= 3*v^2 + 9*gm/r

	spinOrbit = {-# SCC "SO" #-}
		scale (gm^2/(r^3*c^3)*(1+sqrt (1-4*eta))/4*chi) $
				scale nbrac nv + scale (nsbrac) (cross nv s')
					+ scale (-vsbrac) (cross vv s')
		where
			nbrac = 12 * s' <.> (cross nv vv)
			nsbrac = (9 + 3 * sqrt (1-4*eta))*rd
			vsbrac = 7 + sqrt(1-4*eta)

	quadMono = {-# SCC "Q" #-}
		scale (-q*chi^2*3*gc^3*m1^2*m/(2*c^4*r^4)) $
			scale nbrac nv + scale (-2 * nv <.> s') s'
		where
			nbrac = 5 * (nv <.> s')^2 - 1
	r	= norm rv
	v	= norm vv
	s'	= spinCylToSpher s
	nv	= scale (1/r) rv
	rd	= (rv <.> vv)/r
	gm 	= gc * m
	gmperc2 = gm / c ^2
	m 	= m1+m2
	eta 	= m1*m2/m^2

ext: hs

max_stars_repo_name: sageh/AstroHaskell

max_issues_repo_name: sageh/AstroHaskell

max_forks_repo_name: sageh/AstroHaskell


text:
{-# LANGUAGE FlexibleContexts #-}
-- {-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE LambdaCase #-}
module Main where

import           Common                         ( Analysis(..)
                                                , Path(..)
                                                )
import           Control.Monad                  ( foldM
                                                , forM
                                                , forM_
                                                , replicateM
                                                , zipWithM_
                                                )
import           Control.Monad.Except           ( runExceptT )
import           Data.Either                    ( rights )
import           Data.List                      ( unzip4 )
import qualified Data.List                     as L
import           Data.Maybe                     ( catMaybes
                                                , listToMaybe
                                                , mapMaybe
                                                )
import qualified Data.Vector                   as V
import           GHC.Float                      ( int2Double )
import qualified Graphics.Matplotlib           as Plt
import qualified GreedyParser                  as Greedy
import           Inference.Conjugate            ( Hyper
                                                , HyperRep
                                                , Prior(sampleProbs)
                                                , Probs
                                                , ProbsRep
                                                , Trace
                                                , evalTraceLogP
                                                , evalTracePredLogP
                                                , getPosterior
                                                , jeffreysPrior
                                                , traceTrace
                                                , uniformPrior
                                                )
import           Musicology.Pitch               ( SPitch )
import qualified Numeric.Log                   as Log
import           PVGrammar                      ( Edge
                                                , PVAnalysis
                                                , loadAnalysis
                                                , loadInput
                                                , topEdges
                                                )
import           PVGrammar.Parse                ( protoVoiceEvaluator
                                                , pvCountUnrestricted
                                                )
import           PVGrammar.Prob.Simple          ( PVParams(PVParams)
                                                , observeDerivation'
                                                , sampleDerivation'
                                                )
import qualified Parser
import qualified Statistics.Sample             as Stats
import           System.FilePath                ( (<.>)
                                                , (</>)
                                                )
import qualified System.FilePattern            as FP
import qualified System.FilePattern.Directory  as FP
import qualified System.Random.MWC.Probability as MWC
import           System.Random.Stateful         ( StatefulGen
                                                , initStdGen
                                                , newIOGenM
                                                )
import           Text.Pretty.Simple             ( CheckColorTty(..)
                                                , OutputOptions(..)
                                                , defaultOutputOptionsNoColor
                                                , pPrintOpt
                                                )
import           Text.Printf                    ( printf )

prettyPrint :: (Show a) => a -> IO ()
prettyPrint = pPrintOpt
  NoCheckColorTty
  (defaultOutputOptionsNoColor { outputOptionsCompact = True })

dataDir :: FilePath
dataDir = "data/"

main = mainLearn

nBaselines :: Int
nBaselines = 100 -- default: 100

nSamples :: Int
nSamples = 250 -- default: 250

mainLearn :: IO ()
mainLearn = do
  -- initialize
  genPure <- initStdGen
  gen     <- newIOGenM genPure
  genMWC  <- MWC.create -- uses a fixed seed
  let prior = uniformPrior @PVParams
  -- load data
  articleExamples <- loadDir
    (dataDir </> "theory-article")
    ["05b_cello_prelude_1-4", "09a_hinunter", "03_bwv784_pattern"]
  Just bwv939 <- loadItem (dataDir </> "bach" </> "fünf-kleine-präludien")
                          "BWV_0939"
  Just bwv940 <- loadItem (dataDir </> "bach" </> "fünf-kleine-präludien")
                          "BWV_0940"
  let dataset = bwv939 : bwv940 : articleExamples
  -- let dataset = take 3 articleExamples
  putStrLn "list of pieces:"
  forM_ dataset $ \(name, _ana, _trace, _surface) -> do
    putStrLn $ "  " <> name
  -- compute overall posterior
  posteriorTotal <- learn prior dataset
  prettyPrint posteriorTotal
  -- cross validation
  let splits = leaveOneOut dataset
  crossPerps <- forM splits (comparePerNote gen genMWC prior)
  let (logpPriors, logpPosts, counts, basemeans, baselines) =
        L.unzip5 crossPerps
      count           = sum counts
      logppnPrior     = sum logpPriors / count
      logppnTrained   = sum logpPosts / count
      logppnBaseline  = sum basemeans / count
      -- look at each split
      logppns         = zipWith (/) logpPosts counts -- fmap (\(_, logp, count, _, _) -> logp / count) crossPerps
      baselinelogppns = zipWith (\bs n -> (/ n) <$> bs) baselines counts
      testpieces      = fst <$> splits
      testnames       = (\(name, _, _, _) -> name) <$> testpieces
  showSplits logppns baselinelogppns counts testnames
      -- nsplits         = int2Double $ length splits
      -- meanlogppn      = sum logppns / nsplits
  putStrLn $ "prior logppn: " <> show logppnPrior
  putStrLn $ "prior perppn: " <> show (exp $ negate logppnPrior)
  putStrLn $ "overall trained logppn: " <> show logppnTrained
  putStrLn $ "overall trained perppn: " <> show (exp $ negate logppnTrained)
  putStrLn $ "baseline logppn: " <> show logppnBaseline
  putStrLn $ "baseline perppn: " <> show (exp $ negate logppnBaseline)
  putStrLn $ "prior-posterior Δlogppn:" <> show (logppnPrior - logppnTrained)
  putStrLn $ printf "logppn (nats) & %.3f & %.3f & %.3f\\\\"
                    logppnPrior
                    logppnTrained
                    logppnBaseline
  let natsbits = log 2
  putStrLn $ printf "logppn (bits) & %.3f & %.3f & %.3f\\\\"
                    (logppnPrior / natsbits)
                    (logppnTrained / natsbits)
                    (logppnBaseline / natsbits)
  putStrLn $ printf "perppn & %.2f & %.2f & %.2f\\\\"
                    (exp $ negate logppnPrior)
                    (exp $ negate logppnTrained)
                    (exp $ negate logppnBaseline)

testEstimator = do
  (long, short)   <- loadExamples
  articleExamples <- loadDir
    (dataDir </> "theory-article")
    ["05b_cello_prelude_1-4", "09a_hinunter", "03_bwv784_pattern"]
  Just bwv939 <- loadItem (dataDir </> "bach" </> "fünf-kleine-präludien")
                          "BWV_0939"
  Just bwv940 <- loadItem (dataDir </> "bach" </> "fünf-kleine-präludien")
                          "BWV_0940"
  let dataset = bwv939 : bwv940 : articleExamples
  genMWC <- MWC.createSystemRandom
  let prior = uniformPrior @PVParams
  posteriorTotal <- learn prior dataset
  estimates      <- derivationLogProb genMWC posteriorTotal long
  -- let means = meanOfLogs <$> tail (L.inits estimates)
  -- -- print means
  -- Plt.onscreen $ Plt.line [1 .. length means] means
  print $ meanOfLogs estimates

countRareIntDerivs = do
  Just (_, _, _, surface) <- loadItem (dataDir </> "theory-article")
                                      "10c_rare_int"
  count <- Parser.parseSize pvCountUnrestricted surface
  print count

-- loading data
-- ------------

type Piece
  = (String, PVAnalysis SPitch, Trace PVParams, Path [SPitch] [Edge SPitch])

loadItem :: FilePath -> FilePath -> IO (Maybe Piece)
loadItem dir name = do
  ana <- loadAnalysis (dir </> name <.> "analysis.json")
  case ana of
    Left  _err -> pure Nothing
    Right a    -> if anaTop a == PathEnd topEdges
      then do
        surface <- loadInput (dir </> name <.> "musicxml")
        case observeDerivation' (anaDerivation a) of
          Left _err -> do
            putStrLn $ "could not observe trace for " <> name <> ", skipping."
            pure Nothing
          Right trace -> pure $ Just (name, a, trace, surface)
      else do
        putStrLn $ "derivation for " <> name <> " is incomplete, skipping."
        pure Nothing

loadDir :: FilePath -> [String] -> IO [Piece]
loadDir dir exclude = do
  files <- FP.getDirectoryFiles dir ["*.analysis.json"]
  let getName file = FP.match "*.analysis.json" file >>= listToMaybe
      names =
        -- exclude duplicats
        filter (`L.notElem` exclude) $ catMaybes $ getName <$> files
  -- print names
  items <- mapM (loadItem dir) names
  pure $ catMaybes items

loadExamples :: IO (Trace PVParams, Trace PVParams)
loadExamples = do
  Just (_, _, bwv940Trace, _) <- loadItem
    (dataDir </> "bach" </> "fünf-kleine-präludien")
    "BWV_0940"
  Just (_, _, rareTrace, _) <- loadItem (dataDir </> "theory-article")
                                        "10c_rare_int"
  pure (bwv940Trace, rareTrace)

-- learning
-- --------

learn :: Hyper PVParams -> [Piece] -> IO (Hyper PVParams)
learn = foldM train
 where
  train prior (name, _, trace, _) =
    case getPosterior prior trace sampleDerivation' of
      Nothing -> do
        putStrLn $ "couldn't compute posterior for " <> name <> ", skipping."
        pure prior
      Just post -> do
        -- putStrLn $ "learned from " <> name <> "."
        pure post

-- evaluating
-- ----------

leaveOneOut :: [a] -> [(a, [a])]
leaveOneOut dataset = go dataset [] []
 where
  go []       _    splits = splits
  go (x : xs) done splits = go xs (x : done) ((x, xs <> done) : splits)

derivationLogProb
  :: MWC.GenIO -> Hyper PVParams -> Trace PVParams -> IO [Double]
derivationLogProb gen hyper trace = do
  probs <- replicateM nSamples $ MWC.sample (sampleProbs @PVParams hyper) gen
  let estimates = mapMaybe
        (\params -> snd <$> evalTraceLogP params trace sampleDerivation')
        probs
  pure $! estimates

derivationLogProb'
  :: MWC.GenIO -> Hyper PVParams -> Trace PVParams -> IO Double
derivationLogProb' gen hyper trace = do
  estimates <- derivationLogProb gen hyper trace
  pure $! meanOfLogs estimates

meanOfLogs :: (RealFloat b, Foldable t, Functor t) => t b -> b
meanOfLogs logs =
  Log.ln $! Log.sum (Log.Exp <$> logs) / fromIntegral (length logs)

countNotes :: Path [a] b -> Int
countNotes (PathEnd notes       ) = length notes
countNotes (Path notes edges rst) = length notes + countNotes rst

sampleBaselines
  :: (StatefulGen g IO)
  => g
  -> MWC.GenIO
  -> Hyper PVParams
  -> Piece
  -> IO [Double]
sampleBaselines gen genMWC posterior (name, _, _, surface) = do
  derivsTry <- replicateM nBaselines $ runExceptT $ Greedy.parseRandom'
    gen
    protoVoiceEvaluator
    surface
  let derivs = rights derivsTry
  putStrLn $ "  collected " <> show (length derivs) <> " samples for " <> name
  let baselines = flip fmap derivs $ \ana -> do
        trace <- observeDerivation' (anaDerivation ana)
        pure $ derivationLogProb' genMWC posterior trace
  baselines <- traverse sequence baselines -- maaaaagic
  fmap catMaybes $ forM baselines $ \case
    Left  err -> putStrLn err >> pure Nothing
    Right val -> pure $ Just val

summarize :: [Double] -> (Double, Double)
summarize xs = (Stats.mean sample, Stats.stdDev sample)
  where sample = V.fromList xs

comparePerNote
  :: (StatefulGen g IO)
  => g
  -> MWC.GenIO
  -> Hyper PVParams
  -> (Piece, [Piece])
  -> IO (Double, Double, Double, Double, [Double])
comparePerNote gen genMWC prior (test@(tstName, _, tstTrace, tstSurface), train)
  = do
    putStrLn $ "testing on " <> tstName
    posterior <- learn prior train
    -- compute normalized logprobs
    let nnotes = int2Double $ countNotes tstSurface
    logpPrior <- derivationLogProb' genMWC prior tstTrace
    logpPost  <- derivationLogProb' genMWC posterior tstTrace
    let logppnPrior = logpPrior / nnotes
        logppnPost  = logpPost / nnotes
    baselines <- sampleBaselines gen genMWC posterior test
    let (basemean, basestd) = summarize baselines
        basemeanpn          = basemean / nnotes
        basestdpn           = basestd / nnotes
    putStrLn $ "  nnotes: " <> show nnotes
    putStrLn $ "  logppn prior: " <> show logppnPrior <> " nats"
    putStrLn $ "  logppn posterior: " <> show logppnPost <> " nats"
    putStrLn
      $  "  logppn baseline: mean="
      <> show basemeanpn
      <> ", std="
      <> show basestdpn
    putStrLn $ "  baselines top 5: " <> show
      (drop 95 $ (/ nnotes) <$> L.sort baselines)
    putStrLn $ "  perplexity prior: " <> show (exp $ negate logppnPrior)
    putStrLn $ "  perplexity posterior: " <> show (exp $ negate logppnPost)
    putStrLn $ "  Δlogppn: " <> show (logppnPost - logppnPrior) <> " nats"
    pure (logpPrior, logpPost, nnotes, basemean, baselines)

-- plotting
-- --------

a % b = a Plt.% Plt.mp Plt.# b
infixl 5 %

showSplits
  :: [Double] -> [[Double]] -> [Double] -> [String] -> IO (Either String String)
showSplits logppns blogppns counts testpieces =
  Plt.file "splits.svg"
    $ Plt.readData (logppns, blogppns, counts, testpieces)
    % "import numpy as np"
    % "import pandas as pd"
    % "import seaborn as sns"
    % "(logppns, blogppns, counts, pieces) = tuple(map(np.array, data))"
    % "baselines = pd.concat([pd.DataFrame({'logppn': bppns, 'piece': piece}) for bppns, piece in zip(blogppns, pieces)])"
    % "testscores = pd.DataFrame({'logppn': logppns, 'piece': pieces})"
    % "colors = sns.color_palette()"
    % "fig, ax = plot.subplots(figsize=(10,5))"
    % "sns.boxplot(ax=ax, y='piece', x='logppn', data=baselines, whis=(2,98), color=colors[0], fliersize=4)"
    % "sns.scatterplot(ax=ax, y='piece', x='logppn', data=testscores, color=colors[1])"
    % "ax.invert_yaxis()"
    % "fig.tight_layout()"
    % "fig.savefig('splits.pdf')"
    % "fig.savefig('splits.png')"

max_stars_repo_name: DCMLab/protovoices-haskell

max_issues_repo_name: DCMLab/protovoices-haskell

max_forks_repo_name: DCMLab/protovoices-haskell


text:
module Statistics.Quantile.Bench.Gen(
    gauss
  , gamma
  , sampleGaussian
  , sampleGamma
) where

import Control.Applicative
import Control.Monad

import Data.Random (runRVar)
import Data.Random.Distribution
import Data.Random.Distribution.Normal
import qualified Data.Random.Distribution.Gamma as G

import Pipes

import Statistics.Quantile.Types

import System.Random.MWC

gauss :: Double -> Double -> Producer Double IO ()
gauss = runSampler sampleGaussian

gamma :: Double -> Double -> Producer Double IO ()
gamma = runSampler sampleGamma

runSampler :: (GenIO -> Double -> Double -> IO Double)
           -> Double
           -> Double
           -> Producer Double IO ()
runSampler sample a b = do
  gen <- liftIO createSystemRandom
  forever $ do
    x <- liftIO $ sample gen a b
    yield x

sampleGaussian :: GenIO -> Double -> Double -> IO Double
sampleGaussian gen mu sigma = runRVar (normal mu sigma) gen

sampleGamma :: GenIO -> Double -> Double -> IO Double
sampleGamma gen alpha beta = runRVar (G.gamma alpha beta) gen

max_stars_repo_name: fractalcat/slides

max_issues_repo_name: fractalcat/slides

max_forks_repo_name: fractalcat/slides


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

module School.FileIO.Source
( Source(..)
, pool
, source
) where

import Conduit ((.|), ConduitM, MonadResource, mapMC, nullC, sourceFileBS, takeCE)
import Control.Monad ((>=>), unless)
import Control.Monad.Except (MonadError)
import Data.ByteString (ByteString)
import Numeric.LinearAlgebra (Matrix, R)
import School.FileIO.BinConversion (binConversion)
import School.FileIO.ConduitHeader (conduitHeader)
import School.FileIO.FileHeader (FileHeader(..))
import School.FileIO.FilePath (FilePath)
import School.Types.DataType (DataType(..), getSize)
import School.Types.Decoding (binToDouble, binToListInt, binToMatrixDouble)
import School.Types.Error (Error)
import School.Types.LiftResult (LiftResult(..))

class Source a where
  decodeByteString :: FileHeader -> ByteString -> Either Error a

pool :: (Monad m)
     => Int
     -> (ByteString -> m a)
     -> ConduitM ByteString a m ()
pool chunkSize transformer = loop where
  loop = do
    takeCE chunkSize .| mapMC transformer
    isEmpty <- nullC
    unless isEmpty loop

source :: (LiftResult m, MonadError Error m, MonadResource m, Source a)
       => FileHeader
       -> FilePath
       -> ConduitM () a m ()
source header@FileHeader { dataType, cols, rows } path =
  let size = getSize dataType * cols * rows
      byteSource = sourceFileBS path
      cHeader = conduitHeader header
      trans = liftResult . decodeByteString header
  in byteSource .| cHeader .| pool size trans

instance Source (Matrix R) where
  decodeByteString FileHeader{ dataType, cols, rows } =
    binToMatrixDouble dataType rows cols

instance Source Double where
  decodeByteString FileHeader{ dataType } =
    binConversion dataType DBL64B >=> binToDouble

instance Source [Int] where
  decodeByteString FileHeader{ dataType } = binToListInt dataType

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


text:
module Prelude.Containers.Algo (module X) where

import Data.Vector.Lens as X
import Data.IntSet.Lens as X

import Data.Array.Lens as X
import Data.Bits.Lens as X
import Data.ByteString.Lens as X
-- import Data.ByteString.Lazy.Lens as X
-- import Data.ByteString.Strict.Lens as X
import Data.Complex.Lens as X
import Data.Data.Lens as X
import Data.Dynamic.Lens as X
-- import Data.HashSet.Lens as X hiding (setOf, setmapped)  -- conlicts with Data.IntSet.Lens
import Data.List.Lens as X
import Data.Map.Lens as X
import Data.Sequence.Lens as X hiding (sliced)  -- conlicts with Data.Vector.Lens
-- import Data.Set.Lens as X hiding (setOf, setmapped)    -- conlicts with Data.IntSet.Lens
import Data.Text.Lens as X
-- import Data.Text.Lazy.Lens as X
-- import Data.Text.Strict.Lens as X
import Data.Tree.Lens as X
import Data.Typeable.Lens as X
-- import Data.Vector.Generic.Lens as X

max_stars_repo_name: jackblades/lens-prelude

max_issues_repo_name: jackblades/lens-prelude

max_forks_repo_name: jackblades/lens-prelude


text:
{-# LANGUAGE CPP #-}
module Main (main) where

import Prelude hiding (sum)

#if !MIN_VERSION_base(4, 8, 0)
import Control.Applicative ((<$>))
#endif
import Control.Monad (when)
import Data.Foldable (sum)
import Data.Maybe (fromJust)

import System.Environment (getArgs, getProgName)
import System.Random (getStdGen, randoms)

import Text.Printf (printf)

import qualified Data.Vector.Generic as V

import qualified Criterion as C
import qualified Criterion.Types as C
import qualified Statistics.Resampling.Bootstrap as S

import qualified Data.ReedSolomon as RS

runBench :: RS.Encoder -> RS.Matrix -> IO Double
runBench encoder shards = do
    let bench = C.nf (fromJust . RS.encode RS.defaultBackend encoder) shards

    r <- C.benchmark' bench

    let mean = S.estPoint $ C.anMean $ C.reportAnalysis r
        iters = 1.0 / mean
        dataSize = sum $ V.map V.length shards
        throughput = iters * fromIntegral dataSize

    return throughput

main :: IO ()
main = do
    args <- getArgs
    when (length args /= 3) $ do
        progName <- getProgName
        error $ "Usage: " ++ progName ++ " dataShards parityShards dataSize"

    let [dataShards, parityShards, dataSize] = fmap read args

    instructionSet <- RS.simdInstructions
    printf "Using CPU instructions: %s\n" (maybe "None" show instructionSet)

    putStrLn "generating encoder..."
    encoder <- RS.new dataShards parityShards

    putStrLn "generating data..."
    vs <- RS.split encoder =<< V.fromListN dataSize . randoms <$> getStdGen

    throughput <- runBench encoder vs

    printf "Configuration: dataShards %d, parityShards %d, dataSize %d, shardSize %d\n"
        dataShards parityShards dataSize (V.length $ V.head vs)
    printf "Mean throughput: %.2f MB/s\n" (throughput / (1024 * 1024))

max_stars_repo_name: tetcoin-staging/reed-solomon-erasure

max_issues_repo_name: tetcoin-staging/reed-solomon-erasure

max_forks_repo_name: tetcoin-staging/reed-solomon-erasure


text:
{-# LANGUAGE BinaryLiterals            #-}
{-# LANGUAGE CPP                       #-}
{-# LANGUAGE FlexibleContexts          #-}
{-# LANGUAGE FlexibleInstances         #-}
{-# LANGUAGE MultiParamTypeClasses     #-}
{-# LANGUAGE NegativeLiterals          #-}
{-# LANGUAGE NoMonomorphismRestriction #-}
{-# LANGUAGE ScopedTypeVariables       #-}

-- | Tests for the flat module
module Main where

import           Control.Monad
import           Data.Bits
import qualified Data.ByteString          as B
import qualified Data.ByteString.Lazy     as L
import qualified Data.ByteString.Short    as SBS
import           Data.Char
import           Data.Either
import           Data.Flat
import           Data.Flat.Bits
import           Data.Flat.Decoder
import qualified Data.Flat.Encoder        as E
import qualified Data.Flat.Encoder.Prim   as E
import qualified Data.Flat.Encoder.Strict as E
import           Data.Int
import           Data.Proxy
import qualified Data.Sequence            as Seq
import qualified Data.Text                as T
import           Data.Word
import           Numeric.Natural
import           System.Exit
import           Test.Data
import           Test.Data.Arbitrary()
import           Test.Data.Flat
import           Test.Data.Values         hiding (lbs, ns)
import           Test.E
import           Test.E.Arbitrary()
import           Test.E.Flat
import           Test.Tasty
import           Test.Tasty.HUnit
import           Test.Tasty.QuickCheck    as QC hiding (getSize)
import           Data.Flat.Endian
import           Data.FloatCast
import Data.Text.Arbitrary
-- import Test.QuickCheck.Arbitrary
import qualified Data.Complex as B
import qualified Data.Ratio as B

import qualified Data.Map                 as C
import qualified Data.Map.Strict                 as CS
import qualified Data.Map.Lazy                 as CL
import qualified Data.IntMap.Strict as CS
import qualified Data.IntMap.Lazy  as CL

-- import           Data.List
-- import           Data.Ord

#if MIN_VERSION_base(4,9,0)
import qualified Data.List.NonEmpty as BI
#endif

instance Arbitrary UTF8Text where
  arbitrary = UTF8Text <$> arbitrary
  shrink t = UTF8Text <$> shrink (unUTF8 t)

instance Arbitrary UTF16Text where
    arbitrary = UTF16Text <$> arbitrary
    shrink t = UTF16Text <$> shrink (unUTF16 t)

-- instance Flat [Int16]
-- instance Flat [Word8]
-- instance Flat [Bool]

main = do
-- #ifdef ghcjs_HOST_OS
--   print "GHCJS"
-- #endif

  -- printInfo

  -- print $ flat asciiStrT

  mainTest
  -- print $ flatRaw 18446744073709551615::Word64
  -- print $ B.unpack . flat $ (True,0::Word64,18446744073709551615::Word64)
  -- print (2^56::Word64,fromIntegral (1::Word8) `shiftL` 56 :: Word64,(18446744073709551615::Word64) `shiftR` 1)
  -- mainShow
  -- eWord64E id 0b

mainShow = do
  mapM_ (\_ -> generate (arbitrary :: Gen Int) >>= print) [1..10]
  exitFailure

mainTest = defaultMain tests

tests :: TestTree
tests = testGroup "Tests" [
  testPrimitives

  ,testEncDec

  ,testFlat
 ]

testPrimitives = testGroup "conversion/memory primitives" [
   testEndian
  ,testFloatingConvert
  --,testShifts -- ghcjs fails this
  ]

testEncDec = testGroup "encode/decode primitives" [
   testEncodingPrim
  ,testDecodingPrim
#ifdef TEST_DECBITS
  ,testDecBits
#endif
  ]

testFlat = testGroup "flat/unflat" [
  testSize
  ,testLargeEnum
  ,testContainers
  ,flatUnflatRT
  ,flatTests
  ]

-- Data.Flat.Endian tests (to run, need to modify imports and cabal file)
testEndian = testGroup "Endian" [
   conv toBE16 (2^10 + 3)  (2^9+2^8+4)
  ,conv toBE32 (2^18 + 3)  50332672
  ,conv toBE64 (2^34 + 3)  216172782180892672
  ,conv toBE16 0x1234 0x3412
  ,conv toBE32 0x11223344 0x44332211
  ,conv toBE64 0x0123456789ABCDEF 0xEFCDAB8967452301
  ]

testFloatingConvert = testGroup "Floating conversions" [
   conv floatToWord -0.15625 3189768192
  ,conv wordToFloat 3189768192 -0.15625
  ,conv doubleToWord -0.15625 13818169556679524352
  ,conv wordToDouble 13818169556679524352 -0.15625
  ,rt "floatToWord" (prop_float_conv :: RT Float)
  ,rt "doubleToWord" (prop_double_conv :: RT Double)
 ]

-- ghcjs bug on shiftR 0, see: https://github.com/ghcjs/ghcjs/issues/706
testShifts = testGroup "Shifts" $ map tst [0..33]
   where
 tst n = testCase ("shiftR " ++ show n) $
  let val = 4294967295::Word32
      s = val `shift` (-n)
      r = val `shiftR` n
  in r @?= s

-- shR = shiftR
-- shR = unsafeShiftR
shR val 0 = val
shR val n = shift val (-n)

testEncodingPrim = testGroup "Encoding Primitives" [
   encRawWith 1 E.eTrueF [0b10000001]
  ,encRawWith 3 (E.eTrueF >=> E.eFalseF >=> E.eTrueF) [0b10100001]

  -- Depends on endianess
  --,encRawWith 32 (E.eWord32E id $ 2^18 + 3) [3,0,4,0,1]
  -- ,encRawWith 64 (E.eWord64E id $ 0x1122334455667788) [0x88,0x77,0x66,0x55,0x44,0x33,0x22,0x11,1]
  --,encRawWith 65 (E.eTrueF >=> E.eWord64E id (2^34 + 3)) [1,0,0,0,2,0,0,128,129]
  --,encRawWith 65 (E.eFalseF >=> E.eWord64E id (2^34 + 3)) [1,0,0,0,2,0,0,0,129]

  -- Big Endian
  ,encRawWith 32 (E.eWord32BEF  $ 2^18 + 3) [0,4,0,3,1]
  ,encRawWith 64 (E.eWord64BEF  $ 2^34 + 3) [0,0,0,4,0,0,0,3,1]
  ,encRawWith 65 (E.eTrueF >=> E.eWord64BEF (2^34 + 3)) [128,0,0,2,0,0,0,1,129]
  ,encRawWith 65 (E.eFalseF >=> E.eWord64BEF (2^34 + 3))  [0,0,0,2,0,0,0,1,129]
  ]
  where
    encRawWith sz enc exp = testCase (unwords ["encode raw with size",show sz]) $ flatRawWith sz enc @?= exp

conv f v e = testCase (unwords ["conv",sshow v,showB . flat $ v,"to",sshow e]) $ f v @?= e

testDecodingPrim = testGroup "Decoding Primitives" [
   dec ((,,,) <$> dropBits 13 <*> dBool <*> dBool <*> dBool) [0b10111110,0b10011010] ((),False,True,False)
  ,dec ((,,,) <$> dropBits 1 <*> dBE16 <*> dBool <*> dropBits 6) [0b11000000
                                                                 ,0b00000001
                                                                 ,0b01000000] ((),2^15+2,True,())
  ,dec ((,,,) <$> dropBits 1 <*> dBE32 <*> dBool <*> dropBits 6) [0b11000000
                                                                 ,0b00000000
                                                                 ,0b00000000
                                                                 ,0b00000001
                                                                 ,0b01000000] ((),2^31+2,True,())
  ,dec dBE64 [0b10000000
                                                                 ,0b00000000
                                                                 ,0b00000000
                                                                 ,0b00000000
                                                                 ,0b00000000
                                                                 ,0b00000000
                                                                 ,0b00000000
                                                                 ,0b00000010
                                                                 ] (2^63+2)

  ,dec ((,,,) <$> dropBits 1 <*> dBE64 <*> dBool <*> dropBits 6) [0b11000000
                                                                 ,0b00000000
                                                                 ,0b00000000
                                                                 ,0b00000000
                                                                 ,0b00000000
                                                                 ,0b00000000
                                                                 ,0b00000000
                                                                 ,0b00000001
                                                                 ,0b01000000] ((),2^63+2,True,())

  ]
    where
      dec decOp v e = testCase (unwords ["decode",sshow v]) $ unflatRawWith decOp (B.pack v) @?= Right e

testDecBits = testGroup "Decode Bits" $ concat [
   decBitsN dBEBits8
  ,decBitsN dBEBits16
  ,decBitsN dBEBits32
  ,decBitsN dBEBits64
 ]
  where
          -- Test dBEBits8/16/32/64, extraction of up to 8/16/32/bits from various positions
          decBitsN :: forall a. (Num a,FiniteBits a,Show a,Flat a) => (Int -> Get a) -> [TestTree]
          decBitsN dec = let s = finiteBitSize (undefined::a)
                         in [decBits_ dec val numBitsToTake pre | numBitsToTake <- [0 .. s], val <- [0::a ,1+2^(s - 2)+2^(s - 5) ,fromIntegral $ (2^s::Integer) - 1],pre <- [0,1,7]]

          decBits_ :: forall a. (FiniteBits a,Show a,Flat a) => (Int -> Get a) -> a -> Int -> Int -> TestTree
          decBits_ deco val numBitsToTake pre =
            -- a sequence composed by pre zero bits followed by the val and zero bits till the next byte boundary
            let vs = B.pack . asBytes . fromBools $ replicate pre False ++ toBools (asBits val)
                len = B.length vs
                sz = finiteBitSize (undefined::a)
                dec :: Get a
                dec = do
                  dropBits pre
                  r <- deco numBitsToTake
                  dropBits (len*8-numBitsToTake-pre)
                  return r
                -- we expect the first numBitsToTake bits of the value
                expectedD@(Right expected) :: Decoded a = Right $ val `shR` (sz - numBitsToTake) -- ghcjs: shiftR fails, see: https://github.com/ghcjs/ghcjs/issues/706
                actualD@(Right actual) :: Decoded a = unflatRawWith dec vs
            in testCase (unwords ["take",show numBitsToTake,"bits from",show val,"of size",show sz,"with prefix",show pre,"sequence",showB vs,show expected,show actual,show $ val == actual,show $ expected == actual,show $ expected /= actual,show $ show expected == show actual,show $ flat expected == flat actual])
                $ actualD @?= expectedD

testSize = testGroup "Size" $ concat [
  sz () 0
  ,sz True 1
  ,sz One 2
  ,sz Two 2
  ,sz Three 2
  ,sz Four 3
  ,sz Five 3
  ,sz 'a' 8
  ,sz 'à' 16
  ,sz '经' 24
  ,sz (0::Word8) 8
  ,sz (1::Word8) 8
  ,concatMap (uncurry sz) ns
  ,concatMap (uncurry sz) nsI
  ,concatMap (uncurry sz) nsII
  ,sz (1.1::Float) 32
  ,sz (1.1::Double) 64
  ,sz "" 1
  ,sz "abc" (4+3*8)
  ,sz ((),(),Unit) 0
  ,sz (True,False,One,Five) 7
  ,sz map1 7
  ,sz bs (4+3*8)
  ,sz stBS bsSize
  ,sz lzBS bsSize
#ifndef ghcjs_HOST_OS
  ,sz shBS bsSize
#endif
  ,sz tx utf8Size
  ,sz (UTF8Text tx) utf8Size
  ,sz (UTF16Text tx) utf16Size
 ]
   where
    tx = T.pack "txt"
    utf8Size = 8+8+3*32+8
    utf16Size = 8+8+3*16+8
    bsSize = 8+8+3*8+8

sz v e = [testCase (unwords ["size of",sshow v]) $ getSize v @?= e]

-- E258_256 = 11111110 _257 = 111111110 _258 = 111111111
testLargeEnum = testGroup "test enum with more than 256 constructors" $ concat
  [
#ifdef ENUM_LARGE
      sz E258_256 8
    , sz E258_257 9
    , sz E258_258 9

    -- As encodes are inlined, this is going to take for ever if this is compiled with -O1 or -O2
    -- , encRaw (E258_256) [0b11111110]
    -- , encRaw (E258_257) [0b11111111,0b00000000]
    -- , encRaw (E258_258) [0b11111111,0b10000000]
    -- , encRaw (E258_256,E258_257,E258_258) [0b11111110,0b11111111,0b01111111,0b11000000]

    , map trip [E258_1,E258_256,E258_257,E258_258]
    , map trip [E256_1,E256_134,E256_256]
#endif
  ]

testContainers = testGroup "containers" [
    trip longSeq
    , trip dataMap
    , trip listMap
    -- , trip intMap
    ]

flatUnflatRT = testGroup "unflat (flat v) == v"
  [  rt "()" (prop_Flat_roundtrip:: RT ())
    ,rt "Bool" (prop_Flat_roundtrip::RT Bool)
    ,rt "Char" (prop_Flat_roundtrip:: RT Char)
    ,rt "Complex" (prop_Flat_roundtrip:: RT (B.Complex Float))
    ,rt "Either N Bool" (prop_Flat_roundtrip:: RT (Either N Bool))
    ,rt "Either Int Char" (prop_Flat_roundtrip:: RT (Either Int Char))
    ,rt "Int8" (prop_Flat_Large_roundtrip:: RTL Int8)
    ,rt "Int16" (prop_Flat_Large_roundtrip:: RTL Int16)
    ,rt "Int32" (prop_Flat_Large_roundtrip:: RTL Int32)
    ,rt "Int64" (prop_Flat_Large_roundtrip:: RTL Int64)
    ,rt "Int" (prop_Flat_Large_roundtrip:: RTL Int)
    ,rt "[Int16]" (prop_Flat_roundtrip:: RT [Int16])
    ,rt "String" (prop_Flat_roundtrip:: RT String)
#if MIN_VERSION_base(4,9,0)
    ,rt "NonEmpty" (prop_Flat_roundtrip:: RT (BI.NonEmpty Bool))
#endif
    ,rt "Maybe N" (prop_Flat_roundtrip:: RT (Maybe N))
    ,rt "Ratio" (prop_Flat_roundtrip:: RT (B.Ratio Int32))
    ,rt "Word8" (prop_Flat_Large_roundtrip:: RTL Word8)
    ,rt "Word16" (prop_Flat_Large_roundtrip:: RTL Word16)
    ,rt "Word32" (prop_Flat_Large_roundtrip:: RTL Word32)
    ,rt "Word64" (prop_Flat_Large_roundtrip:: RTL Word64)
    ,rt "Word" (prop_Flat_Large_roundtrip:: RTL Word)
    ,rt "Natural" (prop_Flat_roundtrip:: RT Natural)
    ,rt "Integer" (prop_Flat_roundtrip:: RT Integer)
    ,rt "Float" (prop_Flat_roundtrip:: RT Float)
    ,rt "Double" (prop_Flat_roundtrip:: RT Double)

    ,rt "Text" (prop_Flat_roundtrip:: RT T.Text)
    ,rt "UTF8 Text" (prop_Flat_roundtrip:: RT UTF8Text)
    ,rt "UTF16 Text" (prop_Flat_roundtrip:: RT UTF16Text)

    ,rt "ByteString" (prop_Flat_roundtrip:: RT B.ByteString)
    ,rt "Lazy ByteString" (prop_Flat_roundtrip:: RT L.ByteString)
#ifndef ghcjs_HOST_OS
    ,rt "Short ByteString" (prop_Flat_roundtrip:: RT SBS.ShortByteString)
#endif

    ,rt "Map.Strict" (prop_Flat_roundtrip:: RT (CS.Map Int Bool))
    ,rt "Map.Lazy" (prop_Flat_roundtrip:: RT (CL.Map Int Bool))
    ,rt "IntMap.Strict" (prop_Flat_roundtrip:: RT (CS.IntMap Bool))
    ,rt "IntMap.Lazy" (prop_Flat_roundtrip:: RT (CL.IntMap Bool))

    ,rt "Unit" (prop_Flat_roundtrip:: RT Unit)
    ,rt "Un" (prop_Flat_roundtrip:: RT Un )
    ,rt "N" (prop_Flat_roundtrip:: RT N )
    ,rt "E2" (prop_Flat_roundtrip:: RT E2 )
    ,rt "E3" (prop_Flat_roundtrip:: RT E3 )
    ,rt "E4" (prop_Flat_roundtrip:: RT E4 )
    ,rt "E8" (prop_Flat_roundtrip:: RT E8 )
    ,rt "E16" (prop_Flat_roundtrip:: RT E16 )
    ,rt "E17" (prop_Flat_roundtrip:: RT E17 )
    ,rt "E32" (prop_Flat_roundtrip:: RT E32 )
    ,rt "A" (prop_Flat_roundtrip:: RT A )
    ,rt "B" (prop_Flat_roundtrip:: RT B )
    -- ,rt "Tree Bool" (prop_Flat_roundtrip:: RT (Tree Bool))
    -- ,rt "Tree N" (prop_Flat_roundtrip:: RT (Tree N))
    ,rt "List N" (prop_Flat_roundtrip:: RT (List N))
    ]

rt n = QC.testProperty (unwords ["round trip",n])

flatTests = testGroup "flat/unflat Unit tests" $ concat [

  -- Expected errors
   errDec (Proxy::Proxy Bool) [] -- no data
  ,errDec (Proxy::Proxy Bool) [128] -- no filler
  ,errDec (Proxy::Proxy Bool) [128+1,1,2,4,8] -- additional bytes

  ,encRaw () []
  ,encRaw ((),(),Unit) []
  ,encRaw (Unit,'a',Unit,'a',Unit,'a',Unit) [97,97,97]
  ,a () [1]
  ,a True [128+1]
  ,a (True,True) [128+64+1]
  ,a (True,False,True) [128+32+1]
  ,a (True,False,True,True) [128+32+16+1]
  ,a (True,False,True,True,True) [128+32+16+8+1]
  ,a (True,False,True,True,True,True) [128+32+16+8+4+1]
  ,a (True,False,True,True,True,True,True) [128+32+16+8+4+2+1]
  ,a (True,False,True,True,(True,True,True,True)) [128+32+16+8+4+2+1,1]
  ,encRaw (True,False,True,True) [128+32+16]
  ,encRaw ((True,True,False,True,False),(False,False,True,False,True,True)) [128+64+16+1,64+32]
  ,encRaw ('\0','\1','\127') [0,1,127]
  ,encRaw (33::Word32,44::Word32) [33,44]
    --,s (Elem True) [64]
    --,s (NECons True (NECons False (Elem True))) [128+64+32+4]
  ,encRaw (0::Word8) [0]
  ,encRaw (1::Word8) [1]
  ,encRaw (255::Word8) [255]
  ,encRaw (0::Word16) [0]
  ,encRaw (1::Word16) [1]
  ,encRaw (255::Word16) [255,1]
  ,encRaw (256::Word16) [128,2]
  ,encRaw (65535::Word16) [255,255,3]
  ,encRaw (127::Word32) [127]
  ,encRaw (128::Word32) [128,1]
  ,encRaw (129::Word32) [129,1]
  ,encRaw (255::Word32) [255,1]
  ,encRaw (16383::Word32) [255,127]
  ,encRaw (16384::Word32) [128,128,1]
  ,encRaw (16385::Word32) [129,128,1]
  ,encRaw (32767::Word32) [255,255,1]
  ,encRaw (32768::Word32) [128,128,2]
  ,encRaw (32769::Word32) [129,128,2]
  ,encRaw (65535::Word32) [255,255,3]
  ,encRaw (2097151::Word32) [255,255,127]
  ,encRaw (2097152::Word32) [128,128,128,1]
  ,encRaw (2097153::Word32) [129,128,128,1]
  ,encRaw (4294967295::Word32) [255,255,255,255,15]
  ,encRaw (255::Word64) [255,1]
  ,encRaw (65535::Word64) [255,255,3]
  ,encRaw (4294967295::Word64) [255,255,255,255,15]
  ,encRaw (18446744073709551615::Word64)       [255,255,255,255,255,255,255,255,255,1]
  ,encRaw (False,18446744073709551615::Word64) [127,255,255,255,255,255,255,255,255,128,128]
  ,encRaw (255::Word) [255,1]
  ,encRaw (65535::Word) [255,255,3]
  ,encRaw (4294967295::Word) [255,255,255,255,15]
  ,tstI [0::Int8,2,-2]
  ,encRaw (127::Int8) [254]
  ,encRaw (-128::Int8) [255]
  ,tstI [0::Int16,2,-2,127,-128]
  ,tstI [0::Int32,2,-2,127,-128]
  ,tstI [0::Int64,2,-2,127,-128]
  ,encRaw (-1024::Int64) [255,15]
  ,encRaw (maxBound::Word8)       [255]
  ,encRaw (True,maxBound::Word8)  [255,128]
  ,encRaw (maxBound::Word16)      [255,255,3]
  ,encRaw (True,maxBound::Word16) [255,255,129,128]
  ,encRaw (maxBound::Word32)      [255,255,255,255,15]
  ,encRaw (True,maxBound::Word32) [255,255,255,255,135,128]
  ,encRaw (maxBound::Word64)      [255,255,255,255,255,255,255,255,255,1]
  ,encRaw (True,maxBound::Word64) [255,255,255,255,255,255,255,255,255,128,128]
  ,encRaw (minBound::Int64) [255,255,255,255,255,255,255,255,255,1]
  ,encRaw (maxBound::Int64) [254,255,255,255,255,255,255,255,255,1]
  ,tstI [0::Int,2,-2,127,-128]
  ,tstI [0::Integer,2,-2,127,-128,-256,-512]
  ,encRaw (-1024::Integer) [255,15]
  ,encRaw (0::Float)  [0,0,0,0]
  ,encRaw (-2::Float) [0b11000000,0,0,0]
  ,encRaw (0.085::Float) [0b00111101,0b10101110,0b00010100,0b01111011]
  ,encRaw (0::Double)  [0,0,0,0,0,0,0,0]
  ,encRaw (-2::Double) [0b11000000,0,0,0,0,0,0,0]
  ,encRaw (23::Double) [0b01000000,0b00110111,0,0,0,0,0,0]
  ,encRaw (-0.15625::Float)  [0b10111110,0b00100000,0,0]
  ,encRaw (-0.15625::Double) [0b10111111,0b11000100,0,0,0,0,0,0]
  ,encRaw (-123.2325E-23::Double) [0b10111011,0b10010111,0b01000111,0b00101000,0b01110101,0b01111011,0b01000111,0b10111010]
  ,encRaw (Left True :: Either Bool (Double, Double)) [0b01000000]
  ,encRaw (-2.1234E15 :: Double) [195,30,44,226,90,221,64,0]
  ,encRaw (1.1234E-22 :: Double) [59,96,249,241,120,219,249,174]
  ,encRaw ((False,-2.1234E15) :: (Bool,Double)) [97,143,22,113,45,110,160,0,0]
  ,encRaw ((True,-2.1234E15) :: (Bool,Double)) [225,143,22,113,45,110,160,0,0]
  ,encRaw ((-2.1234E15 , 1.1234E-22) :: (Double, Double)) $ [0b11000011,30,44,226,90,221,64,0] ++ [59,96,249,241,120,219,249,174]
  ,encRaw ((True,-2.1234E15 , 1.1234E-22) :: (Bool,Double, Double)) [0b11100001,143,22,113,45,110,160,0,29,176,124,248,188,109,252,215,0]
  ,encRaw (Right (-2.1234E15 , 1.1234E-22) :: Either Bool (Double, Double)) [0b11100001,143,22,113,45,110,160,0,29,176,124,248,188,109,252,215,0]
  ,encRaw (Left True:: Either Bool Direction) [0b01000000]
  ,encRaw (Right West :: Either Bool Direction) [0b11110000]

  ,map trip [minBound,maxBound::Word8]
  ,map trip [minBound,maxBound::Word16]
  ,map trip [minBound,maxBound::Word32]
  ,map trip [minBound,maxBound::Word64]
  ,map trip [minBound::Int8,maxBound::Int8]
  ,map trip [minBound::Int16,maxBound::Int16]
  ,map trip [minBound::Int32,maxBound::Int32]
  ,map trip [minBound::Int64,maxBound::Int64]
  ,map tripShow [0::Float,-0::Float,0/0::Float,1/0::Float]
  ,map tripShow [0::Double,-0::Double,0/0::Double,1/0::Double]
  ,encRaw '\0' [0]
  ,encRaw '\1' [1]
  ,encRaw '\127' [127]
  ,encRaw 'a' [97]
  ,encRaw 'à' [224,1]
  ,encRaw '经' [207,253,1]
  ,[trip [chr 0x10FFFF]]
  ,encRaw Unit []
  ,encRaw (Un False) [0]
  ,encRaw (One,Two,Three) [16+8]
  ,encRaw (Five,Five,Five) [255,128]
    --,s (NECons True (Elem True)) [128+64+16]
  ,encRaw "" [0]
#ifdef LIST_BIT
  ,encRaw "abc" [176,216,172,96]
  ,encRaw [False,True,False,True] [128
                               +32+16
                               +8
                               +2+1,0]
#elif defined(LIST_BYTE)
  ,s "abc" s3
  ,s (cs 600) s600
#endif
    -- Aligned structures
    --,s (T.pack "") [1,0]
    --,s (Just $ T.pack "abc") [128+1,3,97,98,99,0]
    --,s (T.pack "abc") (al s3)
    --,s (T.pack $ cs 600) (al s600)
  ,encRaw map1 [0b10111000]
  ,encRaw (B.pack $ csb 3) (bsl c3)
  ,encRaw (B.pack $ csb 600) (bsl s600)
  ,encRaw (L.pack $ csb 3) (bsl c3)
   -- Long LazyStrings can have internal sections shorter than 255
   --,s (L.pack $ csb 600) (bsl s600)
  ,[trip [1..100::Int16]]

-- See https://github.com/typelead/eta/issues/901
#ifndef ETA_VERSION
  ,[trip longAsciiStrT]
  ,[trip longBoolListT]
#endif

  ,[trip asciiTextT]
  ,[trip english]
  ,[trip "维护和平正"]
  ,[trip (T.pack "abc")]
  ,[trip unicodeText]
  ,[trip unicodeTextUTF8T]

  ,[trip longBS,trip longLBS]
#ifndef ghcjs_HOST_OS
  ,[trip longSBS]
  ,[trip unicodeTextUTF16T]
#endif
  ]
    where

      --al = (1:) -- prealign
      bsl = id -- noalign

      tstI = map ti

      ti v | v >= 0    = testCase (unwords ["Int",show v]) $ teq v (2 * fromIntegral v ::Word64)
           | otherwise = testCase (unwords ["Int",show v]) $ teq v (2 * fromIntegral (-v) - 1 ::Word64)

      teq a b = ser a @?= ser b

              --,testCase (unwords ["unflat raw",sshow v]) $ desRaw e @?= Right v]

      -- Aligned values unflat to the original value, modulo the added filler.
      a v e = [testCase (unwords ["flat",sshow v]) $ ser v @?= e
              ,testCase (unwords ["unflat",sshow v]) $ let Right v' = des e in v @?= v']
      -- a v e = [testCase (unwords ["flat postAligned",show v]) $ ser (postAligned v) @?= e
      --         ,testCase (unwords ["unflat postAligned",show v]) $ let Right (PostAligned v' _) = des e in v @?= v']

encRaw :: forall a. (Show a, Flat a) => a -> [Word8] -> [TestTree]
encRaw v e = [testCase (unwords ["flat raw",sshow v,show . B.unpack . flat $ v]) $ serRaw v @?= e]

trip :: forall a .(Show a,Flat a,Eq a) => a -> TestTree
trip v = testCase (unwords ["roundtrip",sshow v]) $
  -- direct comparison
  (unflat (flat v::B.ByteString)::Decoded a) @?= (Right v::Decoded a)

tripShow :: forall a .(Show a,Flat a,Eq a) => a -> TestTree
tripShow v = testCase (unwords ["roundtrip",sshow v]) $
  -- we use show to get Right NaN == Right NaN
  show (unflat (flat v::B.ByteString)::Decoded a) @?= show (Right v::Decoded a)

-- Test Data
lzBS = L.pack bs
stBS = B.pack bs
bs = [32,32,32::Word8]
s3 = [3,97,98,99,0]
c3a = [3,99,99,99,0] -- Array Word8
c3 = pre c3a
s600 = pre s600a
pre = (1:)
s600a = concat [[255],csb 255,[255],csb 255,[90],csb 90,[0]]
s600B = concat [[55],csb 55,[255],csb 255,[90],csb 90,[200],csb 200,[0]]
longSeq :: Seq.Seq Word8
longSeq = Seq.fromList lbs
longBS = B.pack lbs
longLBS = L.concat $ concat $ replicate 10 [L.pack lbs]
lbs = concat $ replicate 100 [234,123,255,0]
cs n = replicate n 'c' -- take n $ cycle ['a'..'z']
csb = map (fromIntegral . ord) . cs
map1 = C.fromList [(False,True),(True,False)]

ns :: [(Word64, Int)]
ns =  [( (-) (2 ^(i*7)) 1,fromIntegral (8*i)) | i <- [1 .. 10]]

nsI :: [(Int64, Int)]
nsI = nsI_
nsII :: [(Integer, Int)]
nsII = nsI_
nsI_ =  [( (-) (2 ^(((-) i 1)*7)) 1,fromIntegral (8*i)) | i <- [1 .. 10]]

#ifndef ghcjs_HOST_OS
shBS = SBS.toShort stBS
longSBS = SBS.toShort longBS
#endif

sshow = take 80 . show

showB = show . B.unpack

errDec :: forall a . (Flat a, Eq a, Show a) => Proxy a -> [Word8] -> [TestTree]
--errDec _ bs = [testCase "bad decode" $ let ev = (des bs::Decoded a) in ev @?= Left ""]
errDec _ bs = [testCase "bad decode" $ let ev = (des bs::Decoded a) in isRight ev @?= False]

ser :: Flat a => a -> [Word8]
ser = B.unpack . flat

des :: Flat a => [Word8] -> Decoded a
des = unflat

flatRawWith sz enc = B.unpack $ E.strictEncoder (sz+8) (E.Encoding $ enc >=> E.eFillerF)

serRaw :: Flat a => a -> [Word8]
-- serRaw = B.unpack . flatRaw
-- serRaw = L.unpack . flatRaw
serRaw = asBytes . bits

--desRaw :: Flat a => [Word8] -> Decoded a
--desRaw = unflatRaw . L.pack

type RT a = a -> Bool
type RTL a = Large a -> Bool

prop_Flat_roundtrip :: (Flat a, Eq a) => a -> Bool
prop_Flat_roundtrip = roundTripExt

prop_Flat_Large_roundtrip :: (Eq b, Flat b) => Large b -> Bool
prop_Flat_Large_roundtrip (Large x) = roundTripExt x

roundTrip x = unflat (flat x::B.ByteString) == Right x

-- Test roundtrip for both the value and the value embedded between bools
roundTripExt x = roundTrip x && roundTrip (True,x,False)

prop_double_conv d = wordToDouble (doubleToWord d) == d

prop_float_conv d = wordToFloat (floatToWord d) == d

{-
prop_common_unsigned :: (Num l,Num h,Flat l,Flat h) => l -> h -> Bool
prop_common_unsigned n _ = let n2 :: h = fromIntegral n
                           in flat n == flat n2
-}

-- e :: Stream Bool
-- e = unflatIncremental . flat $ stream1

-- el :: List Bool
-- el = unflatIncremental . flat $ infList

-- deflat = unflat

-- b1 :: BLOB UTF8
-- b1 = BLOB UTF8 (preAligned (List255 [97,98,99]))
-- -- b1 = BLOB (preAligned (UTF8 (List255 [97,98,99])))

max_stars_repo_name: AnisimoffNikita/flat

max_issues_repo_name: AnisimoffNikita/flat

max_forks_repo_name: AnisimoffNikita/flat


text:
-- comment
module Main where

import qualified Graphics.Rendering.OpenGL as GL
import Numeric.LinearAlgebra.HMatrix as H hiding (reshape)
import Prelude hiding (init)
import GLCode
import Points

import Sketch

rrand :: Int -> IO (Matrix Double)
rrand n = randn n n

crand :: Int -> IO (Matrix (Complex Double))
crand n = curry toComplex <$> randn n n <*> randn n n

unit :: Floating a => Int -> [a]
unit n = [fromIntegral j/fromIntegral n | j <- [0..n-1]]

initialPoints :: Int -> Vector (Complex Double)
initialPoints num = fromList (n ++ s) where
    n = [(-12.5+25*j) :+ (-20) | j <- unit num]
    s = [(-12.5+25*j) :+ 20 | j <- unit num]

finalPoints :: Int -> Vector (Complex Double)
finalPoints num = fromList (e ++ w) where
    e = [(-20) :+ (-12.5+25*j) | j <- unit num]
    w = [20 :+ (-12.5+25*j) | j <- unit num]

main :: IO ()
main = do
    let num = 20

    u' <- crand (num * 2)
    let u = scale 1.5 u'

    let d = initialPoints num
    let d' = finalPoints num

    let p = diag d
    let p' = u `mul` diag d' `mul` inv u

    mainGifLoop "xxx" 12.0 512 512 1 0 100 $ \time -> do
        let t = time/200

        let eigs = H.toList $ eigenvalues $ scale (1-realToFrac t) p + scale (realToFrac t) p'

        GL.clearColor GL.$= GL.Color4 0.0 0.0 0.0 1
        io $ GL.clear [GL.ColorBuffer]

        let points = [GL.Vertex2 (0.04*realToFrac x) (0.04*realToFrac y) |
                        x :+ y <- eigs] :: [GL.Vertex2 Float]
        setUniform "shader" "pointSize" (8.0 :: Float)
        drawPoint "shader" "vPosition" points

        io GL.flush

max_stars_repo_name: dpiponi/SketchBox

max_issues_repo_name: dpiponi/SketchBox

max_forks_repo_name: dpiponi/SketchBox