text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE RankNTypes #-}
{- |
Module      : Language.Scheme.Parser
Copyright   : Justin Ethier
Licence     : MIT (see LICENSE in the distribution)

Maintainer  : github.com/justinethier
Stability   : experimental
Portability : portable

This module implements parsing of Scheme code.
-}

module Language.Scheme.Parser
    (
      lispDef
    -- *Higher level parsing
    , mainParser
    , readOrThrow
    , readExpr
    , readExprList
    -- *Low level parsing
    , parseExpr
    , parseAtom
    , parseBool
    , parseChar
    , parseOctalNumber
    , parseBinaryNumber
    , parseHexNumber
    , parseDecimalNumber
    , parseNumber
    , parseRealNumber
    , parseRationalNumber
    , parseComplexNumber
    , parseEscapedChar
    , parseString
    , parseVector
    , parseByteVector
    , parseHashTable
    , parseList
    , parseDottedList
    , parseQuoted
    , parseQuasiQuoted
    , parseUnquoted
    , parseUnquoteSpliced
    ) where
import Language.Scheme.Types
import Control.Monad.Except
import Control.Monad.Fail
import Data.Array
import qualified Data.ByteString as BS
import qualified Data.Char as DC
import Data.Complex
import qualified Data.Map
import Data.Ratio
import Data.Word
import Numeric
import Text.ParserCombinators.Parsec hiding (spaces)
import Text.Parsec.Language
import qualified Text.Parsec.Token as P
#if __GLASGOW_HASKELL__ >= 702
import Data.Functor.Identity (Identity)
import Text.Parsec.Prim (ParsecT)
#endif

-- This was added by pull request #63 as part of a series of fixes
-- to get husk to build on ghc 7.2.2
--
-- For now this has been removed to allow husk to support the older
-- GHC 6.x.x series.
--
--import Data.Functor.Identity (Identity)

-- |Language definition for Scheme
lispDef :: LanguageDef ()
lispDef
  = emptyDef
  { P.commentStart   = "#|"
  , P.commentEnd     = "|#"
  , P.commentLine    = ";"
  , P.nestedComments = True
  , P.identStart     = letter <|> symbol
  , P.identLetter    = letter <|> digit <|> symbol
  , P.reservedNames  = []
  , P.caseSensitive  = True
  }

#if __GLASGOW_HASKELL__ >= 702
lexer :: P.GenTokenParser String () Data.Functor.Identity.Identity
#endif
lexer = P.makeTokenParser lispDef

#if __GLASGOW_HASKELL__ >= 702
dot :: ParsecT String () Identity String
#endif
dot = P.dot lexer

#if __GLASGOW_HASKELL__ >= 702
parens :: ParsecT String () Identity a -> ParsecT String () Identity a
#endif
parens = P.parens lexer

#if __GLASGOW_HASKELL__ >= 702
brackets :: ParsecT String () Identity a -> ParsecT String () Identity a
#endif
brackets = P.brackets lexer

#if __GLASGOW_HASKELL__ >= 702
identifier :: ParsecT String () Identity String
#endif
identifier = P.identifier lexer

#if __GLASGOW_HASKELL__ >= 702
whiteSpace :: ParsecT String () Identity ()
#endif
whiteSpace = P.whiteSpace lexer

#if __GLASGOW_HASKELL__ >= 702
lexeme :: ParsecT String () Identity a -> ParsecT String () Identity a
#endif
lexeme = P.lexeme lexer

-- |Match a special character
symbol :: Parser Char
symbol = oneOf "!$%&|*+-/:<=>?@^_~."

-- |Parse an atom (scheme symbol)
parseAtom :: Parser LispVal
parseAtom = do
  atom <- identifier
  if atom == "."
     then pzero -- Do not match this form
     else return $ Atom atom

-- |Parse a boolean
parseBool :: Parser LispVal
parseBool = do _ <- string "#"
               x <- oneOf "tf"
               return $ case x of
                          't' -> Bool True
                          'f' -> Bool False
                          _ -> Bool False

-- |Parse a character
parseChar :: Parser LispVal
parseChar = do
  _ <- try (string "#\\")
  c <- anyChar
  r <- many (letter <|> digit)
  let pchr = c : r
  case pchr of
    "space"     -> return $ Char ' '
    "newline"   -> return $ Char '\n'
    "alarm"     -> return $ Char '\a'
    "backspace" -> return $ Char '\b'
    "delete"    -> return $ Char '\DEL'
    "escape"    -> return $ Char '\ESC'
    "null"      -> return $ Char '\0'
    "return"    -> return $ Char '\n'
    "tab"       -> return $ Char '\t'
    _ -> case (c : r) of
        [ch] -> return $ Char ch
        ('x' : hexs) -> do
            rv <- parseHexScalar hexs
            return $ Char rv
        _ -> pzero

-- |Parse an integer in octal notation, base 8
parseOctalNumber :: Parser LispVal
parseOctalNumber = do
  _ <- try (string "#o")
  sign <- many (oneOf "-")
  num <- many1 (oneOf "01234567")
  case (length sign) of
     0 -> return $ Number $ fst $ head (Numeric.readOct num)
     1 -> return $ Number $ fromInteger $ (*) (-1) $ fst $ head (Numeric.readOct num)
     _ -> pzero

-- |Parse an integer in binary notation, base 2
parseBinaryNumber :: Parser LispVal
parseBinaryNumber = do
  _ <- try (string "#b")
  sign <- many (oneOf "-")
  num <- many1 (oneOf "01")
  case (length sign) of
     0 -> return $ Number $ fst $ head (Numeric.readInt 2 (`elem` "01") DC.digitToInt num)
     1 -> return $ Number $ fromInteger $ (*) (-1) $ fst $ head (Numeric.readInt 2 (`elem` "01") DC.digitToInt num)
     _ -> pzero

-- |Parse an integer in hexadecimal notation, base 16
parseHexNumber :: Parser LispVal
parseHexNumber = do
  _ <- try (string "#x")
  sign <- many (oneOf "-")
  num <- many1 (digit <|> oneOf "abcdefABCDEF")
  case (length sign) of
     0 -> return $ Number $ fst $ head (Numeric.readHex num)
     1 -> return $ Number $ fromInteger $ (*) (-1) $ fst $ head (Numeric.readHex num)
     _ -> pzero

-- |Parser for Integer, base 10
parseDecimalNumber :: Parser LispVal
parseDecimalNumber = do
  _ <- try (many (string "#d"))
  sign <- many (oneOf "-")
  num <- many1 digit
  if (length sign) > 1
     then pzero
     else return $ (Number . read) $ sign ++ num

-- |Parser for a base 10 Integer that will also
--  check to see if the number is followed by
--  an exponent (scientific notation). If so,
--  the integer is converted to a float of the
--  given magnitude.
parseDecimalNumberMaybeExponent :: Parser LispVal
parseDecimalNumberMaybeExponent = do
  num <- parseDecimalNumber
  parseNumberExponent num

-- |Parse an integer in any base
parseNumber :: Parser LispVal
parseNumber = parseDecimalNumberMaybeExponent <|>
              parseHexNumber <|>
              parseBinaryNumber <|>
              parseOctalNumber <?>
              "Unable to parse number"

-- |Parse a floating point number
parseRealNumber :: Parser LispVal
parseRealNumber = do
  sign <- many (oneOf "-+")
  num <- many digit
  _ <- char '.'
  frac <- many1 digit
  let dec = if not (null num)
               then num ++ "." ++ frac
               else "0." ++ frac
  f <- case (length sign) of
     0 -> return $ Float $ fst $ head (Numeric.readFloat dec)
          -- Bit of a hack, but need to support the + sign as well as the minus.
     1 -> if sign == "-"
             then return $ Float $ (*) (-1.0) $ fst $ head (Numeric.readFloat dec)
             else return $ Float $ fst $ head (Numeric.readFloat dec)
     _ -> pzero
  parseNumberExponent f

-- | Parse the exponent section of a floating point number
--   in scientific notation. Eg "e10" from "1.0e10"
parseNumberExponent :: LispVal -> Parser LispVal
parseNumberExponent n = do
  expnt <- many $ oneOf "Ee"
  case (length expnt) of
    0 -> return n
    1 -> do
      num <- try parseDecimalNumber
      case num of
        Number nexp -> buildResult n nexp
        _ -> pzero
    _ -> pzero
 where
  buildResult (Number num) nexp = return $ Float $ (fromIntegral num) * (10 ** (fromIntegral nexp))
  buildResult (Float num) nexp = return $ Float $ num * (10 ** (fromIntegral nexp))
  buildResult _ _ = pzero

-- |Parse a rational number
parseRationalNumber :: Parser LispVal
parseRationalNumber = do
  pnumerator <- parseDecimalNumber
  case pnumerator of
    Number n -> do
      _ <- char '/'
      sign <- many (oneOf "-")
      num <- many1 digit
      if (length sign) > 1
         then pzero
         else do
             let pdenominator = read $ sign ++ num
             if pdenominator == 0
                then return $ Number 0 -- TODO: Prevents a div-by-zero error, but not really correct either
                else return $ Rational $ n % pdenominator
    _ -> pzero

-- |Parse a complex number
parseComplexNumber :: Parser LispVal
parseComplexNumber = do
  lispreal <- (try parseRealNumber <|> try parseRationalNumber <|> parseDecimalNumber)
  let real = case lispreal of
                  Number n -> fromInteger n
                  Rational r -> fromRational r
                  Float f -> f
                  _ -> 0
  _ <- char '+'
  lispimag <- (try parseRealNumber <|> try parseRationalNumber <|> parseDecimalNumber)
  let imag = case lispimag of
                  Number n -> fromInteger n
                  Rational r -> fromRational r
                  Float f -> f
                  _ -> 0 -- Case should never be reached
  _ <- char 'i'
  return $ Complex $ real :+ imag

-- |Parse an escaped character
parseEscapedChar :: forall st .
                    GenParser Char st Char
parseEscapedChar = do
  _ <- char '\\'
  c <- anyChar
  case c of
    'a' -> return '\a'
    'b' -> return '\b'
    'n' -> return '\n'
    't' -> return '\t'
    'r' -> return '\r'
    'x' -> do
        num <- many $ letter <|> digit
        _ <- char ';'
        parseHexScalar num
    _ -> return c

-- |Parse a hexidecimal scalar
parseHexScalar :: MonadFail m => String -> m Char
parseHexScalar num = do
    let ns = Numeric.readHex num
    case ns of
        [] -> fail $ "Unable to parse hex value " ++ show num
        _ -> return $ DC.chr $ fst $ head ns

-- |Parse a string
parseString :: Parser LispVal
parseString = do
  _ <- char '"'
  x <- many (parseEscapedChar <|> noneOf "\"")
  _ <- char '"'
  return $ String x

-- |Parse a vector
parseVector :: Parser LispVal
parseVector = do
  vals <- sepBy parseExpr whiteSpace
  return $ Vector (listArray (0, (length vals - 1)) vals)

-- |Parse a bytevector
parseByteVector :: Parser LispVal
parseByteVector = do
  ns <- sepBy parseNumber whiteSpace
  return $ ByteVector $ BS.pack $ map conv ns
 where
   conv (Number n) = fromInteger n :: Word8
   conv _ = 0 :: Word8

-- |Parse a hash table. The table is either empty or is made up of
--  an alist (associative list)
parseHashTable :: Parser LispVal
parseHashTable = do
  -- This function uses explicit recursion to loop over the parsed list:
  -- As long as it is an alist, the members are appended to an accumulator
  -- so they can be added to the hash table. However, if the input list is
  -- determined not to be an alist, Nothing is returned, letting the parser
  -- know that a valid hashtable was not read.
  let f :: [(LispVal, LispVal)] -> [LispVal] -> Maybe [(LispVal, LispVal)]
      f acc [] = Just acc
      f acc (List [a, b] :ls) = f (acc ++ [(a, b)]) ls
      f acc (DottedList [a] b :ls) = f (acc ++ [(a, b)]) ls
      f _ (_:_) = Nothing
  vals <- sepBy parseExpr whiteSpace
  let mvals = f [] vals
  case mvals of
    Just m -> return $ HashTable $ Data.Map.fromList m
    Nothing -> pzero

-- |Parse a list
parseList :: Parser LispVal
parseList = liftM List $ sepBy parseExpr whiteSpace
-- TODO: wanted to use endBy (or a variant) above, but it causes an error such that dotted lists are not parsed

-- |Parse a dotted list (scheme pair)
parseDottedList :: Parser LispVal
parseDottedList = do
  phead <- endBy parseExpr whiteSpace
  case phead of
    [] -> pzero -- car is required; no match
    _ -> do
      ptail <- dot >> parseExpr
      case ptail of
        DottedList ls l -> return $ DottedList (phead ++ ls) l
        -- Issue #41
        -- Improper lists are tricky because if an improper list ends in a
        -- proper list, then it becomes proper as well. The following cases
        -- handle that, as well as preserving necessary functionality when
        -- appropriate, such as for unquoting.
        --
        -- FUTURE: I am not sure if this is complete, in fact the "unquote"
        -- seems like it could either be incorrect or one special case among
        -- others. Anyway, for the 3.3 release this is good enough to pass all
        -- test cases. It will be revisited later if necessary.
        --
        List (Atom "unquote" : _) -> return $ DottedList phead ptail
        List ls -> return $ List $ phead ++ ls
        {- Regarding above, see
           http://community.schemewiki.org/?scheme-faq-language#dottedapp

           Note, however, that most Schemes expand literal lists occurring in
           function applications, e.g. (foo bar . (1 2 3)) is expanded into
           (foo bar 1 2 3) by the reader. It is not entirely clear whether this
           is a consequence of the standard - the notation is not part of the
           R5RS grammar but there is strong evidence to suggest a Scheme
           implementation cannot comply with all of R5RS without performing this
           transformation. -}
        _ -> return $ DottedList phead ptail

-- |Parse a quoted expression
parseQuoted :: Parser LispVal
parseQuoted = do
  _ <- lexeme $ char '\''
  x <- parseExpr
  return $ List [Atom "quote", x]

-- |Parse a quasi-quoted expression
parseQuasiQuoted :: Parser LispVal
parseQuasiQuoted = do
  _ <- lexeme $ char '`'
  x <- parseExpr
  return $ List [Atom "quasiquote", x]

-- |Parse an unquoted expression (a quasiquotated expression preceded
--  by a comma)
parseUnquoted :: Parser LispVal
parseUnquoted = do
  _ <- try (lexeme $ char ',')
  x <- parseExpr
  return $ List [Atom "unquote", x]

-- |Parse an unquote-spliced expression
parseUnquoteSpliced :: Parser LispVal
parseUnquoteSpliced = do
  _ <- try (lexeme $ string ",@")
  x <- parseExpr
  return $ List [Atom "unquote-splicing", x]

-- FUTURE: should be able to use the grammar from R5RS
-- to make parsing more efficient (mostly by minimizing
-- or eliminating the number of try's below)

-- |Parse an expression
parseExpr :: Parser LispVal
parseExpr =
      try (lexeme parseComplexNumber)
  <|> try (lexeme parseRationalNumber)
  <|> try (lexeme parseRealNumber)
  <|> try (lexeme parseNumber)
  <|> lexeme parseChar
  <|> parseUnquoteSpliced
  <|> do _ <- try (lexeme $ string "#(")
         x <- parseVector
         _ <- lexeme $ char ')'
         return x
  <|> do _ <- try (lexeme $ string "#u8(")
         x <- parseByteVector
         _ <- lexeme $ char ')'
         return x
--  <|> do _ <- try (lexeme $ string "#hash(")
--         x <- parseHashTable
--         _ <- lexeme $ char ')'
--         return x
  <|> try parseAtom
  <|> lexeme parseString
  <|> lexeme parseBool
  <|> parseQuoted
  <|> parseQuasiQuoted
  <|> parseUnquoted
  <|> try (parens parseList)
  <|> parens parseDottedList
  <|> try (brackets parseList)
  <|> brackets parseDottedList
  <?> "Expression"

-- |Initial parser used by the high-level parse functions
mainParser :: Parser LispVal
mainParser = do
    _ <- whiteSpace
    parseExpr

-- |Use a parser to parse the given text, throwing an error
--  if there is a problem parsing the text.
readOrThrow :: Parser a -> String -> ThrowsError a
readOrThrow parser input = case parse parser "lisp" input of
  Left err -> throwError $ Parser err
  Right val -> return val

-- |Parse an expression from a string of text
readExpr :: String -> ThrowsError LispVal
readExpr = readOrThrow mainParser

-- |Parse many expressions from a string of text
readExprList :: String -> ThrowsError [LispVal]
readExprList = readOrThrow (endBy mainParser whiteSpace)

max_stars_repo_name: svenkeidel/husk-scheme

max_issues_repo_name: svenkeidel/husk-scheme

max_forks_repo_name: svenkeidel/husk-scheme


text:
module Math.IRT.Model.Rasch
  ( RaschModel (..)
  ) where

import Statistics.Distribution

import Math.IRT.Internal.Distribution
import Math.IRT.Internal.LogLikelihood
import Math.IRT.Model.FourPLM ( FourPLM(..) )
import Math.IRT.Model.Generic

data RaschModel = RaschModel { difficulty :: !Double
                             } deriving (Show)

instance Distribution RaschModel where
    cumulative = cumulative . toFourPLM

instance ContDistr RaschModel where
    density    = density . toFourPLM
    quantile _ = error "This shouldn't be needed"

instance DensityDeriv RaschModel where
    densityDeriv = densityDeriv . toFourPLM

instance GenericModel RaschModel where
    fromRasch           = RaschModel
    fromOnePLM          = RaschModel
    fromTwoPLM      _ b = RaschModel b
    fromThreePLM  _ b _ = RaschModel b
    fromFourPLM _ b _ _ = RaschModel b

instance LogLikelihood RaschModel where
    logLikelihood b = logLikelihood b . toFourPLM

toFourPLM :: RaschModel -> FourPLM
toFourPLM (RaschModel sb) = FourPLM 1.0 sb 0.0 1.0

ext: hs

max_stars_repo_name: argiopetech/irt

max_issues_repo_name: argiopetech/irt

max_forks_repo_name: argiopetech/irt


text:
{-# LANGUAGE BangPatterns #-}

module FFT where

import Prelude hiding ((++))
import qualified Data.ByteString.Lazy.Char8      as BS
import qualified Data.ByteString.Lex.Lazy.Double as BS
import qualified Data.Vector as V
import Data.Vector ((!),(++),Vector)
import Data.List hiding ((++))
import Data.Complex
import System.Environment

headerSize bs = (+) 4 (BS.length . BS.concat . Prelude.take 4 . BS.lines $ bs)

dropHeader bs = BS.drop (headerSize bs) $ bs

main = do
   [f] <- getArgs
   s   <- BS.readFile f
   print . V.map (round. realPart) . ifft_CT. fft_CT . parse $ (dropHeader s)

dft :: Vector (Double) -> Vector (Complex Double)
dft xr = V.map (\k -> V.sum (V.imap (arg k) xr)) (numvec n)
 where
   n = V.length xr
   nf = fromIntegral n
   numvec n = V.enumFromStepN 0 1 n
   arg k i x = (x:+0) * exp (-2 * pi * kf * ifl * (0:+1)/nf)
             where
               kf = fromIntegral k
               ifl = fromIntegral i

idft :: Vector (Complex Double) -> Vector (Complex Double)
idft xs = V.map (\k ->(/) (V.sum (V.imap (arg k) xs)) nf) (numvec n)
 where
   n = V.length xs
   nf = fromIntegral n
   numvec n = V.enumFromStepN 0 1 n
   arg k i x = x * exp (2 * pi * kf * ifl * (0:+1)/nf)
             where
               kf = fromIntegral k
               ifl = fromIntegral i

fft_CT xs = if V.length xs == 1 then
             dft xs
           else
             let xtop = V.zipWith3 fft_top (fft_CT xse) (fft_CT xso) (V.enumFromStepN 0 1 m)

                 xbottom = V.zipWith3 fft_bottom (fft_CT xse) (fft_CT xso) (V.enumFromStepN 0 1 m)
                 xse = V.map (xs !) (V.enumFromStepN 0 2 m)
                 xso = V.map (xs !) (V.enumFromStepN 1 2 m)
                 n = V.length xs
                 nf = fromIntegral n
                 m = n `div` 2
                 fft_top xe xo k = xe + xo * exp (-2 * (0:+1) * pi * k / nf)
                 fft_bottom xe xo k = xe - xo * exp (-2 * (0:+1) * pi * k / nf)
             in
                 xtop ++ xbottom

ifft_CT xs = if V.length xs == 1 then
              idft xs
           else
              let xtop = V.zipWith3 fft_top (ifft_CT xse) (ifft_CT xso) (V.enumFromStepN 0 1 m)
                  xtop_scaled = V.map (/2) xtop
                  xbottom = V.zipWith3 fft_bottom (ifft_CT xse) (ifft_CT xso) (V.enumFromStepN 0 1 m)
                  xbottom_scaled = V.map (/2) xbottom
                  xse = V.map (xs !) (V.enumFromStepN 0 2 m)
                  xso = V.map (xs !) (V.enumFromStepN 1 2 m)
                  n = V.length xs
                  nf = fromIntegral n
                  m = n `div` 2
                  fft_top xe xo k = xe + xo * exp (2 * (0:+1) * pi * k / nf)
                  fft_bottom xe xo k = xe - xo * exp (2 * (0:+1) * pi * k / nf)
             in
                xtop_scaled ++ xbottom_scaled

-- Fill a new vector from a file containing a list of numbers.
parse = V.unfoldr step
 where
   step !s = case BS.readDouble s of
     Nothing       -> Nothing
     Just (!k, !t) -> Just (k, BS.tail t)

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
{-# LANGUAGE ConstraintKinds  #-}
{-# LANGUAGE FlexibleContexts #-}

module Advection1D
where

import qualified Basis
import           Data.List
import qualified Element
import           Flux
import           Mesh
import           Numeric.LinearAlgebra
--import           Numeric.LinearAlgebra.HMatrix
import           Quadrature            (integrate)
import qualified ShapeFcns

-- Synonym for numeric/fractional constraints
type FrElNuFi a = (Fractional a,Element a,Numeric a,Field a)

-- Element stiffness integrand
stiffnessIntegrand :: (Element.Element e,ShapeFcns.ShapeFcn s,Basis.Basis b,FrElNuFi a) => e a -> s b -> [a] -> Matrix a
--stiffnessIntegrand elem shpFcn xi = fromLists [[
--  ShapeFcns.n shpFcn xi i [0] * atIndex (Element.dndx elem shpFcn xi j) 0 * Element.computeJacobianDet elem shpFcn xi
--  | i <- idxRange] | j <- idxRange]
stiffnessIntegrand elem shpFcn xi = fromLists [[
  (atIndex (Element.dndx elem shpFcn xi i) 0 * atIndex (Element.dndx elem shpFcn xi j) 0 + ShapeFcns.n shpFcn xi i [0] * ShapeFcns.n shpFcn xi j [0]) * Element.computeJacobianDet elem shpFcn xi
  | i <- idxRange] | j <- idxRange]
  where
    idxRange = [0..(Element.getNumNodes elem - 1)]

-- Element mass integrand
massIntegrand :: (Element.Element e,ShapeFcns.ShapeFcn s,Basis.Basis b,FrElNuFi a) => e a -> s b -> [a] -> Matrix a
massIntegrand elem shpFcn xi = fromLists [[
  ShapeFcns.n shpFcn xi i [0] * ShapeFcns.n shpFcn xi j [0] * Element.computeJacobianDet elem shpFcn xi
  | i <- idxRange] | j <- idxRange]
  where
    idxRange = [0..(ShapeFcns.getShapeFcnOrder shpFcn)]

-- Computes the element flux matrix assuming an upwind flux.
-- Assumes that advSpeed is positive.
elemFluxMatrix :: FrElNuFi a => a -> Matrix a
elemFluxMatrix advSpeed = scale advSpeed $ fromLists [[-1.0, 0.0, 0.0],[0.0, 0.0, 1.0]]

-- Element mass and stiffness matrices
-- Must multiply by the determinant of the Jacobian here!!!
elemMatrices :: (Element.Element e,ShapeFcns.ShapeFcn s,Basis.Basis b,FrElNuFi a) => e a -> s b -> Int -> (Matrix a,Matrix a)
elemMatrices elem shpFcn ngpts = (stiffnessMat, massMat)
  where
    stiffnessMat = integrate 1 ngpts (stiffnessIntegrand elem shpFcn)
    massMat      = integrate 1 ngpts (massIntegrand elem shpFcn)

-- Function which computes the dimensions of the global matrices
globalMatrixDimension :: Int -> Int
globalMatrixDimension numNodes = 2 * numNodes - 2

-- Function which takes a local element matrix along with an element and transforms it to global
localMatToGlobal :: (Element.Element e,FrElNuFi a) => e a -> Matrix a -> [((Int, Int), a)]
localMatToGlobal elem elemMat = concat [[((Mesh.globalNodeNum en i, Mesh.globalNodeNum en j), atIndex elemMat (i,j) ) | i <- [0..dim-1]] | j <- [0..dim-1]]
  where
    en = Element.getElementNumber elem
    dim = (fst . size) elemMat

-- Function which does the same as the above, but with contributions due to fluxes.
localFluxMatToGlobal :: (Element.Element e,FrElNuFi a) => e a -> Matrix a -> [((Int, Int), a)]
localFluxMatToGlobal elem elemMat = concat [[((Mesh.globalNodeNum en i, Mesh.globalNodeNum en j - 1), atIndex elemMat (i,j) ) | i <- [0..dim-1]] | j <- [0..dim]]
  where
    en = Element.getElementNumber elem
    dim = (fst . size) elemMat

-- Function which assembles global stiffness and mass matrices
-- This function takes a mesh and returns a tuple which corresponds to the
-- global stiffness and global mass matrices, respectively.
-- THIS ONLY WORKS WITH THE DISCONTINUOUS MESH...  IT DOES NOT SUM DUPLICATE ENTRIES.
assembleGlobalMatrices :: (Element.Element e, ShapeFcns.ShapeFcn s,Basis.Basis b, FrElNuFi a) => Mesh e a -> s b -> Int -> a -> [[((Int,Int),a)]]
assembleGlobalMatrices grid shpFcn ngpts advSpd = [globalK, globalM, globalF]
  where
    elems = Mesh.getMeshElements grid
    nnodes = length (concatMap Element.getElementNodes elems)
    elemMats = map (\elem -> elemMatrices elem shpFcn ngpts) elems
    elemK = map fst elemMats
    elemM = map snd elemMats
    elemF = elemFluxMatrix advSpd
    globalK = concat (zipWith localMatToGlobal elems elemK)
    globalM = concat (zipWith localMatToGlobal elems elemM)
    globalF = filter (\entry -> (snd . fst) entry >= 0 && (snd . fst) entry < nnodes) (concatMap (`localFluxMatToGlobal` elemF) elems)

-- This function computes the right hand side of the semi-discrete statement.
-- udot_j = (M_ij)^(-1) ( a K_ij - F_ij ) u_j.
-- This is what needs to be integrated point-by-point.
uDot :: (Floating a,Field a,Num (Vector a)) => Matrix a -> Matrix a -> Matrix a -> a -> Vector a -> Vector a
uDot globalM globalK globalF advSpd u = inv globalM #> ((scale advSpd globalK - globalF) #> u)

-- This function advances the solution one step in time.
advance :: (Floating a,Field a,Num (Vector a)) => Matrix a -> Matrix a -> Matrix a -> a -> Vector a -> a -> Vector a
advance globalM globalK globalF advSpd u dt = scale dt (uDot globalM globalK globalF advSpd u - u)

-- This function applies boundary conditions to the semi-discrete system.
applyBCs _ = []

-- This function computes the solution as a function of time.
{-
solveAdvection numElements icFunc tf = do {
  linBasis <- Lagrange 1;
  shpFcn   <- TensorProduct linBasis;
  grid     <- genDMesh 0.0 1.0 numElements StructElem;
  }
-}


text:
{-# LANGUAGE CPP              #-}
{-# LANGUAGE DataKinds        #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeFamilies     #-}

{-# OPTIONS_GHC -fno-warn-missing-signatures #-}

-----------------------------------------------------------------------------
{- |
Module      :  Numeric.LinearAlgebra.Tests.Properties
Copyright   :  (c) Alberto Ruiz 2008
License     :  BSD3
Maintainer  :  Alberto Ruiz
Stability   :  provisional

Testing properties.

-}

module Test.Numeric.LinearAlgebra.Tests.Properties (
    dist, (|~|), (~~), (~:), Aprox((:~)), (~=),
    zeros, ones,
    square,
    unitary,
    hermitian,
    wellCond,
    positiveDefinite,
    upperTriang,
    upperHessenberg,
    luProp,
    invProp,
    pinvProp,
    detProp,
    nullspaceProp,
--    bugProp,
    svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4,
    svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7,
    eigProp, eigSHProp, eigProp2, eigSHProp2,
    geigProp,
    qrProp, rqProp, rqProp1, rqProp2, rqProp3,
    hessProp,
    schurProp1, schurProp2,
    cholProp, exactProp,
    expmDiagProp,
    multProp1, multProp2,
    subProp,
    linearSolveProp, linearSolvePropH, linearSolveProp2,

    -- Binary properties
    vectorBinaryRoundtripProp
  , staticVectorBinaryRoundtripProp
  , matrixBinaryRoundtripProp
  , staticMatrixBinaryRoundtripProp
  , staticVectorBinaryFailProp
) where

import           Numeric.LinearAlgebra.HMatrix hiding (Testable)
import qualified Numeric.LinearAlgebra.Static  as Static
import           Test.QuickCheck

import           Data.Binary
import           Data.Binary.Get               (runGet)
import           Data.Either                   (isLeft)
#if MIN_VERSION_base(4,11,0)
import           Prelude                       hiding ((<>))
#endif

import           Debug.Trace

(~=) :: Float -> Float -> Bool
a ~= b = abs (a - b) < 1e-10

trivial :: Testable a => Bool -> a -> Property
trivial = (`classify` "trivial")

-- relative error
dist :: (Num a, Normed a) => a -> a -> Float
dist = relativeError norm_Inf

infixl 4 |~|
a |~| b = a :~1~: b
--a |~| b = dist a b < 10^^(-10)

a ~~ b = fromList a |~| fromList b

data Aprox a = (:~) a Int
-- (~:) :: (Normed a, Num a) => Aprox a -> a -> Bool
a :~n~: b = dist a b < 10^^(-n)

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

square m = rows m == cols m

-- orthonormal columns
orthonormal m = tr m <> m |~| ident (cols m)

unitary m = square m && orthonormal m

hermitian m = square m && m |~| tr m

wellCond m = rcond m > 1/100

positiveDefinite m = minimum (toList e) > 0
    where (e,_v) = eigSH m

upperTriang m = rows m == 1 || down == z
    where down = fromList $ concat $ zipWith drop [1..] (toLists (tr m))
          z = konst 0 (size down)

upperHessenberg m = rows m < 3 || down == z
    where down = fromList $ concat $ zipWith drop [2..] (toLists (tr m))
          z = konst 0 (size down)

zeros (r,c) = reshape c (konst 0 (r*c))

ones (r,c) = zeros (r,c) + 1

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

luProp m =
  -- trace ("x: " ++ show m)
  -- trace ("m: " ++ show m ++ " =~ " ++ show (p <> l <> u) ++ " ... " ++ show (m |~| p <> l <> u))
  -- trace ("f (det p) |~| f s: " ++ show (f (det p)) ++ " = " ++ show  (f s) ++ " ... " ++ show (f (det p) |~| f s))
  m |~| p <> l <> u && f (det p) |~| f s
    where (l,u,p,s) = lu m
          f x = fromList [x]

invProp m = m <> inv m |~| ident (rows m)

pinvProp m =  m <> p <> m |~| m
           && p <> m <> p |~| p
           && hermitian (m<>p)
           && hermitian (p<>m)
    where p = pinv m

detProp m =
  trace ("d1: " ++ show (s d1, s d2))
  s d1 |~| s d2
    where d1 = det m
          d2 = det' * det q
          det' = product $ toList $ takeDiag r
          (q,r) = qr m
          s x = fromList [x]

nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c)
                                     && orthonormal n)
    where n = nullspaceSVD (Left (1*peps)) m (rightSV m)
          nl = toColumns n
          r = rows m
          c = cols m - rank m

------------------------------------------------------------------
{-
-- testcase for nonempty fpu stack
-- uncommenting unitary' signature eliminates the problem
bugProp m = m |~| u <> real d <> tr v && unitary' u && unitary' v
    where (u,d,v) = svd m
          -- unitary' :: (Num (Vector t), Field t) => Matrix t -> Bool
          unitary' a = unitary a
-}
------------------------------------------------------------------

-- fullSVD
svdProp1 m = m |~| u <> real d <> tr v && unitary u && unitary v
  where
    (u,s,v) = svd m
    d = diagRect 0 s (rows m) (cols m)

svdProp1a svdfun m = m |~| u <> real d <> tr v && unitary u && unitary v
  where
    (u,s,v) = svdfun m
    d = diagRect 0 s (rows m) (cols m)

svdProp1b svdfun m = unitary u && unitary v
  where
    (u,_,v) = svdfun m

-- thinSVD
svdProp2 thinSVDfun m
    =  m |~| u <> diag (real s) <> tr v
    && orthonormal u && orthonormal v
    && size s == min (rows m) (cols m)
  where
    (u,s,v) = thinSVDfun m

-- compactSVD
svdProp3 m = (m |~| u <> real (diag s) <> tr v
             && orthonormal u && orthonormal v)
  where
    (u,s,v) = compactSVD m

svdProp4 m' = m |~| u <> real (diag s) <> tr v
           && orthonormal u && orthonormal v
           && (size s == r || r == 0 && size s == 1)
  where
    (u,s,v) = compactSVD m
    m = fromBlocks [[m'],[m']]
    r = rank m'

svdProp5a m = all (s1|~|) [s3,s5] where
    s1       = singularValues (m :: Matrix Float)
--  s2       = svRd m
    (_,s3,_) = svd m
--  (_,s4,_) = svdRd m
    (_,s5,_) = thinSVD m
--  (_,s6,_) = thinSVDRd m

svdProp5b m = all (s1|~|) [s3,s5] where
    s1       = singularValues (m :: Matrix (Complex Float))
--  s2       = svCd m
    (_,s3,_) = svd m
--  (_,s4,_) = svdCd m
    (_,s5,_) = thinSVD m
--  (_,s6,_) = thinSVDCd m

svdProp6a m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'
  where
    (u,s,v) = svd (m :: Matrix Float)
    (s',v') = rightSV m
    (u',s'') = leftSV m

svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'
  where
    (u,s,v) = svd (m :: Matrix (Complex Float))
    (s',v') = rightSV m
    (u',s'') = leftSV m

svdProp7 m = s |~| s' && u |~| u' && v |~| v' && s |~| s'''
  where
    (u,s,v) = svd m
    (s',v') = rightSV m
    (u',_s'') = leftSV m
    s''' = singularValues m

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

eigProp m = complex m <> v |~| v <> diag s
    where (s, v) = eig m

eigSHProp m = m <> v |~| v <> real (diag s)
              && unitary v
              && m |~| v <> real (diag s) <> tr v
    where (s, v) = eigSH' m

eigProp2 m = fst (eig m) |~| eigenvalues m

eigSHProp2 m = fst (eigSH' m) |~| eigenvaluesSH' m

geigProp a b = a' <> v <> diag betas' |~| b' <> v <> diag alphas
    where (alphas, betas, v) = geig a b
          betas' = complex betas
          a' = complex a
          b' = complex b

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

qrProp m = q <> r |~| m && unitary q && upperTriang r
    where (q,r) = qr m

rqProp m = r <> q |~| m && unitary q && upperTriang' r
    where (r,q) = rq m

rqProp1 m = r <> q |~| m
    where (r,q) = rq m

rqProp2 m = unitary q
    where (_r,q) = rq m

rqProp3 m = upperTriang' r
    where (r,_q) = rq m

upperTriang' r = upptr (rows r) (cols r) * r |~| r
    where upptr f c = build (f,c) $ \r' c' -> if r'-t > c' then 0 else 1
              where t = fromIntegral (f-c)

hessProp m = m |~| p <> h <> tr p && unitary p && upperHessenberg h
    where (p,h) = hess m

schurProp1 m = m |~| u <> s <> tr u && unitary u && upperTriang s
    where (u,s) = schur m

schurProp2 m = m |~| u <> s <> tr u && unitary u && upperHessenberg s -- fixme
    where (u,s) = schur m

cholProp m = m |~| tr c <> c && upperTriang c
    where c = chol (trustSym m)

exactProp m = chol (trustSym m) == chol (trustSym (m+0))

expmDiagProp m = expm (logm m) :~ 7 ~: complex m
    where logm = matFunc log

-- reference multiply
mulH a b = fromLists [[ doth ai bj | bj <- toColumns b] | ai <- toRows a ]
    where doth u v = sum $ zipWith (*) (toList u) (toList v)

multProp1 p (a,b) = (a <> b) :~p~: (mulH a b)

multProp2 p (a,b) = (tr (a <> b)) :~p~: (tr b <> tr a)

linearSolveProp f m = f m m |~| ident (rows m)

linearSolvePropH f m = f m (unSym m) |~| ident (rows (unSym m))

linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b)
    where q = min (rows a) (cols a)
          b = a <> x
          wc = rank a == q

subProp m = m == (conj . tr . fromColumns . toRows) m

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

vectorBinaryRoundtripProp :: Vector Float -> Bool
vectorBinaryRoundtripProp vec = decode (encode vec) == vec

staticVectorBinaryRoundtripProp :: Static.R 5 -> Bool
staticVectorBinaryRoundtripProp vec =
  let
    decoded = decode (encode vec) :: Static.R 500
  in
    Static.extract decoded == Static.extract vec

matrixBinaryRoundtripProp :: Matrix Float -> Bool
matrixBinaryRoundtripProp mat = decode (encode mat) == mat

staticMatrixBinaryRoundtripProp :: Static.L 100 200 -> Bool
staticMatrixBinaryRoundtripProp mat =
  let
    decoded = decode (encode mat) :: Static.L 100 200
  in
    (Static.extract decoded) == (Static.extract mat)

staticVectorBinaryFailProp :: Static.R 20 -> Bool
staticVectorBinaryFailProp vec =
  let
    decoded = runGet get (encode vec) :: Either String (Static.R 50)
  in
    isLeft decoded

max_stars_repo_name: schnecki/hmatrix-float

max_issues_repo_name: schnecki/hmatrix-float

max_forks_repo_name: schnecki/hmatrix-float


text:
module Haskov where

import System.Random
import Data.Tuple (swap)
import Data.List (foldl')
import Data.Sequence (Seq, (|>))
import qualified Data.Sequence as Seqq
import Data.Map.Strict (Map, keys, foldlWithKey', foldrWithKey', mapWithKey)
import qualified Data.Map.Strict as Map
import Numeric.LinearAlgebra (scale, luSolve, luPacked)
--import qualified Numeric.LinearAlgebra as Lin
import Numeric.LinearAlgebra.Data (Matrix, matrix, col, cols, ident, tr)
import qualified Numeric.LinearAlgebra.Data as Dat
import Numeric.LinearAlgebra.HMatrix (sumElements, cmap)
import qualified Numeric.LinearAlgebra.HMatrix as Hma

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

-- Markov Type --
type IndexMap a = Map a Int
type ReverseIndexMap a = Map Int a
type HMatrixMap = Map (Int, Int) Double

data Markov a = Markov {  imap :: IndexMap a
                        , hmap :: HMatrixMap }

instance (Show a, Ord a) => Show (Markov a) where
    show haskov = show $ toList haskov

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

-- Query --
-- Safe query
lookup :: (Ord a) => a -> a -> Markov a -> Maybe Double
lookup i j (Markov imap hmap)
    | (Map.member i imap) && (Map.member j imap) = Just (hmap Map.! ((imap Map.!i), (imap Map.! j)))
    | otherwise                                  = Nothing

-- Returns number of states
size :: (Ord a) => Markov a -> Int
size (Markov imap hmap) = Map.size imap

-- Returns true if the chain is empty
null :: (Ord a) => Markov a -> Bool
null (Markov imap hmap) = Map.null imap

-- Returns true if input is a state of the chain
member :: (Ord a) => a -> Markov a -> Bool
member i (Markov imap hmap)
    | (Map.member i imap) = True
    | otherwise = False

-- Returns true if input is not a state of the chain
notMember :: (Ord a) => a -> Markov a -> Bool
notMember i haskov = not (member i haskov)

-- States of the Markov chain
states :: (Ord a) => Markov a -> [a]
states (Markov imap hmap) = keys imap

statesI :: (Ord a) => Markov a -> IndexMap a
statesI (Markov imap hmap) = imap

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

-- Construction --
empty :: Markov a
empty = Markov (Map.empty) (Map.empty)

insert :: (Ord a) => a -> a -> Double -> Markov a -> Markov a
insert i j d (Markov imap hmap)
    | hasI && (not hasJ) = Markov (Map.insert j s imap) (Map.insert (imap Map.! i, s) d hmap)
    | (not hasI) && hasJ = Markov (Map.insert i s imap) (Map.insert (s, imap Map.! j) d hmap)
    | (not hasI) && (not hasJ) && (i /= j) = Markov (Map.insert j (s+1) (Map.insert i s imap)) (Map.insert (s, (s+1)) d hmap)
    | (not hasI) && (not hasJ) && (i == j) = Markov (Map.insert i s imap) (Map.insert (s, s) d hmap)
    | otherwise          = Markov imap (Map.insert (imap Map.! i, imap Map.! j) d hmap)
    where
        hasI = Map.member i imap
        hasJ = Map.member j imap
        s = Map.size imap

insertWith :: (Ord a) => (Double -> Double -> Double) -> a -> a -> Double -> Markov a -> Markov a
insertWith f i j d (Markov imap hmap)
    | hasI && (not hasJ) = Markov (Map.insert j s imap) (Map.insert (imap Map.! i, s) d hmap)
    | (not hasI) && hasJ = Markov (Map.insert i s imap) (Map.insert (s, imap Map.! j) d hmap)
    | (not hasI) && (not hasJ) && (i /= j) = Markov (Map.insert j (s+1) (Map.insert i s imap)) (Map.insert (s, (s+1)) d hmap)
    | (not hasI) && (not hasJ) && (i == j) = Markov (Map.insert i s imap) (Map.insert (s, s) d hmap)
    | otherwise          = Markov imap (Map.insertWith f (imap Map.! i, imap Map.! j) d hmap)
    where
        hasI = Map.member i imap
        hasJ = Map.member j imap
        s = Map.size imap

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

-- Lists --
fromList :: (Ord a) => [((a, a), Double)] -> Markov a
fromList list =
    let imap = foldl' toimap Map.empty list
        hmap = foldl' (tohmap imap) Map.empty list
    in Markov imap hmap

toList :: (Ord a) => Markov a -> [((a, a), Double)]
toList (Markov imap hmap) =
    let rIMap = reverseIMap imap
    in foldrWithKey' (combineImapHmap rIMap) [] hmap

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

-- Maps --
fromMap :: (Ord a) => Map (a, a) Double -> Markov a
fromMap mapp =
    let imap = foldlWithKey' toimap' Map.empty mapp
        hmap = foldlWithKey' (tohmap' imap) Map.empty mapp
    in Markov imap hmap

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

-- Matrices --
hmatrix :: (Ord a) => Markov a -> Matrix Double
hmatrix (Markov imap hmap) =
    let m = Map.size imap
        zeros = matrix m (replicate (m*m) 0)
        hlist = Map.toList hmap
    in  Dat.accum zeros (+) hlist

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

-- Chains --

walkFrom :: (Ord a, Show a) => a -> Int -> Markov a -> StdGen -> IO [a]
walkFrom start n haskov gen = do
    let mat = hmatrix haskov
        imap = statesI haskov
        rimap = transMap imap
        nextRow = mat Dat.? [(imap Map.! start)]
    return (start : (steps rimap mat nextRow n gen))

walk :: (Ord a, Show a) => Int -> Markov a -> StdGen -> IO [a]
walk n haskov gen = do
    let mat = hmatrix haskov
        ss = tr . steady $ mat
        imap = statesI haskov
        rimap = transMap imap
    return (steps rimap mat ss n gen)

steps :: (Ord a, Show a) => ReverseIndexMap a -> Matrix Double -> Matrix Double -> Int -> StdGen -> [a]
steps _ _ _ 0 _ = []
steps rimap mat row n gen
    | (sumElements row) == 0 = []
    | otherwise = choice : steps rimap mat newRow (n-1) (snd rand)
    where
        rand = random gen
        index = randomStep row (fst rand) 0.0 0
        choice = rimap Map.! index
        newRow = mat Dat.? [index]

randomStep :: Matrix Double -> Double -> Double -> Int -> Int
randomStep row rand total j
    | newTotal < rand = randomStep row rand newTotal newJ
    | otherwise       = j
    where
        newTotal = total + (row Dat.! 0 Dat.! j)
        newJ = (j + 1) `mod` (cols row)

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

-- Steady State --
steadyState :: (Ord a) =>  Markov a -> [(a, Double)]
steadyState haskov = zip (states haskov) (concat . Dat.toLists $ steady (hmatrix haskov))

steady :: Matrix Double -> Matrix Double
steady mat =
    let m = Dat.rows mat -- Number of states
        q = tr $ ident m - mat -- Transpose (I-matrix size m - hmatrix)
        e = col ((replicate (m-1) 0) ++ [1.0]) -- col of zeros size m, last element is 1.0
        x = luSolve (luPacked q) e -- Solve linear system of Q and e (Qx = e)
    in  scale (1 / (sumElements x)) x  -- Divide x by sum of elements

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

-- Normalize --
normalize :: (Ord a) => Markov a -> Markov a
normalize (Markov imap hmap) =
    let ss = foldlWithKey' sumSeq Seqq.empty hmap
    in Markov imap (mapWithKey (calcNorm ss) hmap)

calcNorm :: Seq Double -> (Int, Int) -> Double -> Double
calcNorm ss (i, j) n = n / (Seqq.index ss i)

sumSeq :: Seq Double -> (Int, Int) -> Double -> Seq Double
sumSeq s (i, j) n =
    if Seqq.length s < (i+1)
    then s |> n
    else Seqq.adjust (+n) i s

norm :: Matrix Double -> Matrix Double
norm row = cmap (/ (sumElements row)) row

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

---- Mapping --
--map :: (Ord a) => (Double -> Double) -> Markov a -> Markov a
--map f (Markov imap hmap) = Markov imap (fmap f hmap)

--mapWithState :: (Ord a) => (a -> Double -> Double) -> Markov a -> Markov a
--mapWithState f (Markov imap hmap) =
    --let rIMap = reverseIMap imap

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

-- Delete/Update --
--delete :: (Ord a) => a -> a -> Markov a -> Markov a
--delete i j (Markov imap hmap) =

adjust :: (Ord a) => (Double -> Double) -> a -> a -> Markov a -> Markov a
adjust f i j (Markov imap hmap) =
    if Map.member i imap && Map.member j imap
    then Markov imap (Map.adjust f (imap Map.! i, imap Map.! j) hmap)
    else Markov imap hmap

-- Helper Functions --
-- for fromList
toimap :: (Ord a) => IndexMap a -> ((a, a), Double) -> IndexMap a
toimap acc ((i, j), d)
    | i /= j && Map.notMember i acc && Map.notMember j acc = Map.insert j (Map.size acc+1) (Map.insert i (Map.size acc) acc)
    | Map.notMember i acc = Map.insert i (Map.size acc) acc
    | Map.notMember j acc = Map.insert j (Map.size acc) acc
    | otherwise = acc

-- For fromMap
toimap' :: (Ord a) => IndexMap a -> (a, a) -> Double -> IndexMap a
toimap' acc (i, j) d
    | i /= j && Map.notMember i acc && Map.notMember j acc = Map.insert j (Map.size acc+1) (Map.insert i (Map.size acc) acc)
    | Map.notMember i acc = Map.insert i (Map.size acc) acc
    | Map.notMember j acc = Map.insert j (Map.size acc) acc
    | otherwise = acc

-- For fromList
tohmap :: (Ord a) => IndexMap a -> HMatrixMap -> ((a, a), Double) -> HMatrixMap
tohmap imap acc ((i, j), n) = Map.insert (imap Map.! i, imap Map.! j) n acc

-- For fromMap
tohmap' :: (Ord a) => IndexMap a -> HMatrixMap -> (a, a) -> Double -> HMatrixMap
tohmap' imap acc (i, j) n = Map.insert (imap Map.! i, imap Map.! j) n acc

combineImapHmap :: (Ord a) => Map Int a -> (Int, Int) -> Double -> [((a, a), Double)] -> [((a, a), Double)]
combineImapHmap rIMap (i, j) d acc = ((rIMap Map.! i, rIMap Map.! j), d) : acc

reverseIMap :: (Ord a) => Map a Int -> Map Int a
reverseIMap imap = foldrWithKey' reverseIMap' Map.empty imap
    where reverseIMap' state i m = Map.insert i state m

-- swap map key and index
transMap :: (Ord k, Ord a) => Map k a -> Map a k
transMap = Map.fromList . map swap . Map.toList

-- Machine Epsilon --

machineE :: Double
machineE = until ((== 1) . (+1)) (/2) 1

max_stars_repo_name: francisdb/haskov

max_issues_repo_name: francisdb/haskov

max_forks_repo_name: francisdb/haskov


text:
module ForeignHABQT where

import Data.Bifunctor (bimap)
import Data.Complex (Complex(..))
import Data.List (unzip)
import qualified Data.Vector as V
import Foreign
import Foreign.C
import HABQTlib.Data
import HABQTlib.MeasurementProcessing (PurePOVM)
import qualified Numeric.LinearAlgebra as LA

tdm :: DensityMatrix
tdm = DensityMatrix $ LA.real $ (2 LA.>< 2) [1 ..]

type CMatrix = [[CDouble]]

convertDM :: DensityMatrix -> (CMatrix, CMatrix)
convertDM (DensityMatrix dm) =
  let (rp, ip) = LA.fromComplex dm
      conv = (map . map) realToFrac . LA.toLists
   in (conv rp, conv ip)

unmarshallSV :: Dim -> Ptr CDouble -> Ptr CDouble -> IO PureStateVector
unmarshallSV dim rPtr iPtr = do
  rl <- peekArray dim rPtr
  il <- peekArray dim iPtr
  let toHM = (dim LA.>< 1) . map realToFrac
  return . PureStateVector $ LA.toComplex (toHM rl, toHM il)

convertWPSV :: WeighedPureStateVector -> ([CDouble], [CDouble])
convertWPSV (WeighedPureStateVector (w, PureStateVector sv)) =
  let sv' = LA.scale (sqrt w :+ 0) sv
      (rp, ip) = LA.fromComplex sv'
      conv = map realToFrac . LA.toList . head . LA.toColumns
   in (conv rp, conv ip)

convertPOVM :: PurePOVM -> ([CDouble], [CDouble])
convertPOVM povm =
  let svs = V.toList povm
      csvs = map convertWPSV svs
   in bimap concat concat . unzip $ csvs

marshallPOVM :: PurePOVM -> Ptr CDouble -> Ptr CDouble -> IO ()
marshallPOVM povm rPtr iPtr = do
  let (r, i) = convertPOVM povm
  pokeArray rPtr r
  pokeArray iPtr i

marshallDM :: DensityMatrix -> Ptr CDouble -> Ptr CDouble -> IO ()
marshallDM dm rPtr iPtr = do
  let (r, i) = convertDM dm
  pokeArray rPtr $ concat r
  pokeArray iPtr $ concat i

max_stars_repo_name: Belinsky-L-V/HABQT

max_issues_repo_name: Belinsky-L-V/HABQT

max_forks_repo_name: Belinsky-L-V/HABQT


text:
module Tests.ComplexFromIntegral where
import Data.Complex

runTest = do
  putStrLn $ show $ p' where
    p' = zipWith (*) (map (fromIntegral) [1..3]) p
    p = [1 :+ 0, 2 :+ 0, 3 :+ 0]

max_stars_repo_name: freizl/haste-compiler

max_issues_repo_name: freizl/haste-compiler

max_forks_repo_name: freizl/haste-compiler


text:
-- |
-- Module      : SpriteClip.Events
-- Description :
-- Copyright   : (c) Jonatan H Sundqvist, 2015
-- License     : MIT
-- Maintainer  : Jonatan H Sundqvist
-- Stability   : experimental|stable
-- Portability : POSIX (not sure)
--

-- Created September 7 2015

-- TODO | - Encode all logic in types (state machine?)
--        - Preferences and keymaps
--        - Sophisticated command dispatch (eg. composite events, overlapping events)

-- SPEC | -
--        -

--------------------------------------------------------------------------------------------------------------------------------------------
-- GHC Pragmas
--------------------------------------------------------------------------------------------------------------------------------------------
{-# LANGUAGE OverloadedStrings #-}

--------------------------------------------------------------------------------------------------------------------------------------------
-- API
--------------------------------------------------------------------------------------------------------------------------------------------
module SpriteClip.Events (module SpriteClip.Events) where

--------------------------------------------------------------------------------------------------------------------------------------------
-- We'll need these
--------------------------------------------------------------------------------------------------------------------------------------------
import Data.Complex
import Data.Functor
import Data.Char (toLower)
import Data.IORef
import Data.Maybe (fromMaybe)

import qualified Data.Set  as S
import qualified Data.Text as T
import qualified Data.Map  as M

import Text.Printf
import Control.Monad (liftM, when, forM_, forM, unless)

import System.FilePath (takeExtension)

import qualified Graphics.Rendering.Cairo as Cairo --
import           Graphics.UI.Gtk          as Gtk   --

import qualified Southpaw.Interactive.Application     as App
import qualified Southpaw.Cartesian.Plane.BoundingBox as BBox
import Southpaw.Cartesian.Plane.Utilities

import SpriteClip.Types
import SpriteClip.Lenses
import SpriteClip.Render             as Render
import qualified SpriteClip.Commands as Commands

--------------------------------------------------------------------------------------------------------------------------------------------
-- Functions
--------------------------------------------------------------------------------------------------------------------------------------------
-- Event configurations --------------------------------------------------------------------------------------------------------------------
-- |
-- widgetAcceptsDNDEvents ::
widgetAcceptsDNDEvents :: WidgetClass self => self -> IO ()
widgetAcceptsDNDEvents widget = do
    dragDestSet widget [DestDefaultMotion, DestDefaultDrop] [ActionCopy]
    -- dragDestAddTextTargets  widget
    -- dragDestAddImageTargets widget
    dragDestAddURITargets   widget

-- |
-- widgetAddDNDListeners ::
widgetAddDNDListeners :: WidgetClass widget => App.App AppState -> widget -> IO ()
widgetAddDNDListeners app widget = do
  widget `on` dragBegin        $ ondragbegin    app
  widget `on` dragMotion       $ ondragmotion   app -- TODO: Why won't this work?
  widget `on` dragDrop         $ ondragdrop     app
  widget `on` dragDataReceived $ ondragreceived app
  return ()

-- Rendering -------------------------------------------------------------------------------------------------------------------------------
-- |
-- TODO: Figure out canvasSize
ondraw :: AppState -> Cairo.Render ()
ondraw = Render.render

-- Animation -------------------------------------------------------------------------------------------------------------------------------
-- |
onanimate :: Gtk.DrawingArea -> IORef AppState -> IO Bool
onanimate cvs stateref = Gtk.widgetQueueDraw cvs >> return True

-- Keyboard --------------------------------------------------------------------------------------------------------------------------------
-- |
onkeydown :: IORef AppState -> EventM EKey Bool
onkeydown stateref = do
  key <- eventKeyName
  appstate <- Cairo.liftIO $ readIORef stateref

  Cairo.liftIO $ unless (S.member key $ appstate^.keys) $ do
    -- TODO: Needs heavy refactoring
    -- TODO: Use foldM (?)
    forM (Commands.dispatch (T.unpack key) (S.map T.unpack . S.insert key $ appstate^.keys) Commands.commands) $ \command -> do
      putStrLn "Running command"
      appstate <- readIORef stateref
      newstate <- command appstate
      writeIORef stateref newstate
      -- newstate <- maybe (return appstate) ($ appstate) $ M.lookup (S.fromList [key]) Commands.commands -- TODO: Factour out
    modifyIORef stateref (keys %~ S.insert (T.toLower key))

  --
  w <- eventWindow
  Cairo.liftIO $ do
    cursor <- cursorNew $ case key of
      "a" -> Pencil
      "b" -> Hand1
      _  -> BlankCursor

    -- children <- containerGetChildren w
    -- widgetGrabFocus
    -- forM_ [w] (liftM (maybe (return ()) (const $ putStrLn "Setting cursor")) . widgetGetWindow)
    forM_ [w] (flip drawWindowSetCursor (Just cursor))

  return True

-- |
onkeyup :: IORef AppState -> EventM EKey Bool
onkeyup stateref = do
  key <- eventKeyName
  Cairo.liftIO $ do
    modifyIORef stateref (keys %~ S.delete (T.toLower key))
    printf "Key %s was released.\n" (T.unpack key)
  return False

-- Mouse -----------------------------------------------------------------------------------------------------------------------------------
-- |
onmousemotion :: IORef AppState -> EventM EMotion Bool
onmousemotion stateref = do
  (mx, my) <- eventCoordinates
  Cairo.liftIO $ modifyIORef stateref (mouse .~ (mx:+my))
  return False

-- |
onmousedown :: IORef AppState -> EventM EButton Bool
onmousedown stateref = do
  -- TODO: Check which button it is
  -- mouse <- Cairo.liftIO $ readIORef cutouts
  (mx, my) <- eventCoordinates
  Cairo.liftIO $ modifyIORef stateref (pin .~ Just (mx:+my))

  -- Flesh out cursor behaviour (ie. grab hand when hover over marker)
  -- Cairo.liftIO $ putStrLn "Setting icon"
  -- w <- eventWindow
  -- cursor <- Cairo.liftIO $ Gtk.cursorNew Gtk.Hand1
  -- Cairo.liftIO $ Gtk.drawWindowSetCursor w Nothing -- (Just cursor)

  return True -- Keep sending motion-notify events

-- |
onmouseup :: IORef AppState -> EventM EButton Bool
onmouseup stateref = do
  (mx, my) <- eventCoordinates
  -- TODO: Flip corners if necessary (no negative sizes)
  -- ((makebox topleft (cutouts ^. _2)) :)
  Cairo.liftIO $ do
    modifyIORef stateref (\s -> maybe s (\click -> s & cutouts %~ (BBox.fromCorners click (mx:+my) :)) (s ^. pin))
    s <- readIORef stateref
    modifyIORef stateref (pin .~ Nothing)

  return True

-- Drag'n'Drop -----------------------------------------------------------------------------------------------------------------------------
-- |
acceptURIDrops :: IO ()
acceptURIDrops = do
  return ()

-- |
acceptImageDrops :: IO ()
acceptImageDrops = do
  return ()

-- |
-- ondragbegin ::
ondragbegin :: App.App state -> DragContext -> IO ()
ondragbegin app context = do
  -- dragSetIconName
  Cairo.liftIO $ putStrLn "Drag begin"
  return ()

-- |
ondragstatus :: App.App state -> DragContext -> Maybe DragAction -> TimeStamp -> IO ()
ondragstatus app context maction timestamp = do
  Cairo.liftIO $ putStrLn "Drag status"
  return ()

-- |
-- ondragmotion ::
ondragmotion :: App.App AppState -> DragContext -> Point -> TimeStamp -> IO Bool
ondragmotion app context pt timestamp = do
  Cairo.liftIO $ putStrLn "Dragging..."
  atom <- atomNew ("String" :: String)
  s <- dragGetData (App._window app) context atom timestamp
  return True

-- |
-- ondragdrop ::
ondragdrop :: App.App state -> DragContext -> Point -> TimeStamp -> IO Bool
ondragdrop app context pt timestamp = do
  putStrLn "Dropping data. Catch!"
  return True

-- |
-- ondatareceived ::
ondragreceived :: App.App AppState -> DragContext -> Point -> InfoId -> TimeStamp -> SelectionDataM ()
ondragreceived app context pos infoid timestamp = do
  Cairo.liftIO $ putStrLn "Receing data..."
  s <- selectionDataGetURIs
  Cairo.liftIO $ print s
  Cairo.liftIO $ maybe
      (putStrLn "Didn't understand the drop")
      (\uris -> Cairo.liftIO $ when (takeExtension (uris !! 0) == ".png") $ do
         newsheet <- Cairo.imageSurfaceCreateFromPNG . drop (length ("file:///" :: String)) $ uris !! 0
         modifyIORef (App._state app) (sheet .~ newsheet)

         (cx:+cy) <- imageSurfaceSize newsheet >>= return . dotmap floor
         windowResize (App._window app) cx cy)
      -- (mapM_ (printf "Understood, here it is: <%s>.\n"))
      (s :: Maybe [String])

-- Size & Position -------------------------------------------------------------------------------------------------------------------------

-- Admin -----------------------------------------------------------------------------------------------------------------------------------
-- |
onrealised :: App.App AppState -> IO ()
onrealised app = do
  putStrLn "On realised"

max_stars_repo_name: SwiftsNamesake/SpriteClip

max_issues_repo_name: SwiftsNamesake/SpriteClip

max_forks_repo_name: SwiftsNamesake/SpriteClip


text:
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleContexts #-}
module Main where

import qualified AnalyticsLaop as LM
import qualified AnalyticsHMatrix as HM
import qualified AnalyticsMatrix as DM
import qualified LAoP.Matrix.Type as LM
import qualified Numeric.LinearAlgebra as HM
import qualified Numeric.LinearAlgebra.Data as HM
import qualified Numeric.LinearAlgebra.HMatrix as HM
import qualified Data.Matrix as DM
import GHC.TypeLits
import Data.Proxy
import Criterion.Main
import Test.QuickCheck

randomBDLaop ::
             forall jn en .
             ( KnownNat jn,
               KnownNat en,
               LM.FromLists Int (LM.FromNat (jn + 1)) (LM.FromNat 3),
               LM.FromLists Int (LM.FromNat (jn + 1)) (),
               LM.FromLists Int (LM.FromNat (en + 1)) (LM.FromNat 5),
               LM.FromLists Int (LM.FromNat (en + 1)) (LM.FromNat 3),
               LM.FromLists Int (LM.FromNat (en + 1)) (LM.FromNat 2)
             ) => Gen (LM.BD jn en)
randomBDLaop = do
  let jn' = fromInteger (natVal (Proxy :: Proxy jn)) + 1
      en' = fromInteger (natVal (Proxy :: Proxy en)) + 1
  descs <- vectorOf (jn' * 3) (elements [0, 1])
  salary <- vectorOf jn' (choose (0, 10000))
  code <- vectorOf (jn' * 3) (elements [0, 1])
  let jobs = LM.JT {
      LM.jDesc   = LM.fromLists (buildList descs jn'),
      LM.jSalary = LM.fromLists (buildList salary jn'),
      LM.jCode   = LM.fromLists (buildList code jn')
          }
  name <- vectorOf (en' * 5) (elements [0, 1])
  job <- vectorOf (en' * 3) (elements [0, 1])
  country <- vectorOf (en' * 2) (elements [0, 1])
  branch <- vectorOf (en' * 2) (elements [0, 1])
  let employees = LM.ET {
      LM.eName    = LM.fromLists (buildList name en'),
      LM.eJob     = LM.fromLists (buildList job en'),
      LM.eCountry = LM.fromLists (buildList country en'),
      LM.eBranch  = LM.fromLists (buildList branch en')
          }
  return LM.BD { LM.jobsTable = jobs, LM.employeesTable = employees }

randomBDHMatrix ::
             forall jn en .
             ( KnownNat jn,
               KnownNat en
             ) => Gen (HM.BD jn en)
randomBDHMatrix = do
  let jn' = fromInteger (natVal (Proxy :: Proxy jn)) + 1
      en' = fromInteger (natVal (Proxy :: Proxy en)) + 1
  descs <- vectorOf (jn' * 3) (elements [0, 1])
  salary <- vectorOf jn' (choose (0, 10000))
  code <- vectorOf (jn' * 3) (elements [0, 1])
  let jobs = HM.JT {
      HM.jDesc   = HM.fromLists (buildList descs jn'),
      HM.jSalary = HM.fromLists (buildList salary jn'),
      HM.jCode   = HM.fromLists (buildList code jn')
          }
  name <- vectorOf (en' * 5) (elements [0, 1])
  job <- vectorOf (en' * 3) (elements [0, 1])
  country <- vectorOf (en' * 2) (elements [0, 1])
  branch <- vectorOf (en' * 2) (elements [0, 1])
  let employees = HM.ET {
      HM.eName    = HM.fromLists (buildList name en'),
      HM.eJob     = HM.fromLists (buildList job en'),
      HM.eCountry = HM.fromLists (buildList country en'),
      HM.eBranch  = HM.fromLists (buildList branch en')
          }
  return HM.BD { HM.jobsTable = jobs, HM.employeesTable = employees }

randomBDMatrix ::
             forall jn en .
             ( KnownNat jn,
               KnownNat en
             ) => Gen (DM.BD jn en)
randomBDMatrix = do
  let jn' = fromInteger (natVal (Proxy :: Proxy jn)) + 1
      en' = fromInteger (natVal (Proxy :: Proxy en)) + 1
  descs <- vectorOf (jn' * 3) (elements [0, 1])
  salary <- vectorOf jn' (choose (0, 10000))
  code <- vectorOf (jn' * 3) (elements [0, 1])
  let jobs = DM.JT {
      DM.jDesc   = DM.fromLists (buildList descs jn'),
      DM.jSalary = DM.fromLists (buildList salary jn'),
      DM.jCode   = DM.fromLists (buildList code jn')
          }
  name <- vectorOf (en' * 5) (elements [0, 1])
  job <- vectorOf (en' * 3) (elements [0, 1])
  country <- vectorOf (en' * 2) (elements [0, 1])
  branch <- vectorOf (en' * 2) (elements [0, 1])
  let employees = DM.ET {
      DM.eName    = DM.fromLists (buildList name en'),
      DM.eJob     = DM.fromLists (buildList job en'),
      DM.eCountry = DM.fromLists (buildList country en'),
      DM.eBranch  = DM.fromLists (buildList branch en')
          }
  return DM.BD { DM.jobsTable = jobs, DM.employeesTable = employees }

buildList [] _ = []
buildList l r  = take r l : buildList (drop r l) r

setupEnv1 = do
  -- HMatrix
  bd1 <- generate (resize 1 (randomBDHMatrix @999 @999))
  -- Matrix
  bd2 <- generate (resize 1 (randomBDMatrix @999 @999))
  -- Laop
  bd3 <- generate (resize 1 (randomBDLaop @999 @999))
  return (bd1, bd2, bd3)

benchmark :: IO ()
benchmark = do
  print "Starting benchmarks..."
  defaultMain [
   env setupEnv1 $ \ ~(bd1, bd2, bd3) -> bgroup "DB query" [
   bgroup "1000 entries" [
     bench "hmatrix" $ nf HM.query bd1
   , bench "matrix" $ nf DM.query bd2
   , bench "laop" $ nf LM.query bd3
   ] ] ]

main :: IO ()
main = benchmark


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE OverlappingInstances #-}
{-# LANGUAGE EmptyDataDecls #-}

module TestBase(
    GEMMTestData(..),
    StrassenTestData(..),
    TRSMTestData(..),
    PotrfTestData(..),
    Tag2(),
    eqMatrix,
    handleResult,
    Truncable(..)

) where

import           Data.Complex (Complex ((:+)))
import Test.QuickCheck (Result(..), Arbitrary(arbitrary), Testable(property), Gen, choose, vectorOf, elements)
import Data.Tagged (Tagged(Tagged, unTagged))
import Data.Default (Default(def))
import System.Random (Random)
import Debug.Trace (trace)
import Control.Applicative (liftA)

import FPNLA.Matrix
       (Matrix(fromList_m, dim_m, elem_m),
        Matrix(generate_m, transpose_m, zipWith_m))
import FPNLA.Matrix.Instances.HMatrix ()
import FPNLA.Operations.Parameters (Elt, TriangType(..), TransType(..), UnitType(..), StratCtx)
import FPNLA.Operations.LAPACK.Strategies (HMatrixBindSeq, CholLLVSeq, CholLLVPar_Repa, CholLLVBlocksSeq, newNullContext, newSqrBlockContext)
import FPNLA.Operations.BLAS.Strategies (DefPar_ST, CBindSeq, DefSeq, BlocksSeq, BlocksByRows, BlocksByCols, DefPar_MP, BlocksPar_MP, StrassenPar_MP, ColsPar_MP, DefPar_ST, DefPar_Acc, BlocksPar_ST, StrassenPar_ST, ColsPar_ST, newBlockContext, newStrassenContext)
import FPNLA.Utils (iif)

-- Links interesantes:
-- http://www.haskell.org/haskellwiki/Scoped_type_variables
-- http://www.haskell.org/ghc/docs/7.4.1/html/users_guide/other-type-extensions.html seccion 7.11.7
-- http://hackage.haskell.org/trac/ghc/ticket/3485

blockSize :: Int
blockSize = 2

-- Valor de epsilon para hacer las comparaciones:
epsilon :: Double
epsilon = 1e-5

class Truncable a where
    inEpsilonRange :: Real e => e -> a -> a -> Bool

instance (Floating n, RealFrac n) => Truncable n where
    inEpsilonRange eps n1 n2 = expN1 == expN2 && abs(mantN1 - mantN2) <= epsFractional
        where
            epsFractional = fromRational $ toRational eps
            shift n d1 = n * (10 ** fromIntegral d1)
            getExponent r
                    | r == 0 = 0
                    | r < 0 = getExponent (abs r)
                    | r > 10 = 1 + getExponent (r / 10)
                    | r < 1 = (-1) + getExponent (r * 10)
                    | otherwise = 0
            expN1 :: Integer = getExponent n1
            expN2 :: Integer = getExponent n2
            mantN1 = shift n1 (-expN1)
            mantN2 = shift n2 (-expN2)

--
instance (RealFrac n, Truncable n) => Truncable (Complex n) where
    inEpsilonRange eps (r1:+i1) (r2:+i2) = inEpsilonRange eps r1 r2 && inEpsilonRange eps i1 i2

-- Compara una matriz elemento a elemento truncando los d√≠gitos. Si falla imprime los √≠ndices del primer elemento que fall√≥:
eqMatrix :: (Show e, Eq e, Truncable e, Matrix m e) => m e -> m e -> Bool
eqMatrix m1 m2 = dimm1 == dimm2 && eqElems 0 0
    where
        dimm1 = dim_m m1
        dimm2 = dim_m m2
        (m,n) = dimm1
        getElem = {--truncDigits truncValue $--} elem_m
        eqElems i j
          | j >= n = eqElems (i+1) 0
          | i >= m = True
          -- | getElem i j m1 /= getElem i j m2 = trace ("(i,j) = (" ++ show i ++ "," ++ show j ++ ")") False
          | not $ inEpsilonRange epsilon (getElem i j m1) (getElem i j m2) = trace ("(i,j) = (" ++ show i ++ "," ++ show j ++ ")" ++ " / m1 = " ++ show (getElem i j m1) ++ " / m2 = " ++ show (getElem i j m2)) False
          | otherwise =  eqElems i (j+1)

arbSDPMatrix :: (Elt e, Arbitrary e, Random e, Matrix m e) => Int -> Gen (m e)
arbSDPMatrix n = do
    la <- vectorOf (n*n) $ choose (0, 1)
    ld <- vectorOf n $ choose (100, 200)
    let a = fromList_m n n la
        d = generate_m n n (\i j -> iif (i == j) (ld !! i) 0)
        s = zipWith_m (+) a (transpose_m a)
        r = zipWith_m (+) s d
    return r

handleResult :: Monad m => Result -> m ()
handleResult (Success{}) = return ()
handleResult (GaveUp{}) = fail "GaveUp"
handleResult (Failure{}) = fail "Failure"
handleResult (NoExpectedFailure{}) = fail "NoExpectedFailure"

--
instance Arbitrary (TransType a) where
    arbitrary = elements [NoTrans undefined, Trans undefined, ConjTrans undefined]
--
instance Arbitrary (a -> TransType a) where
    arbitrary = elements [NoTrans, Trans, ConjTrans]
instance Arbitrary (a -> TriangType a) where
    arbitrary = elements [Lower, Upper]
instance Arbitrary (a -> UnitType a) where
    arbitrary = elements [Unit, NoUnit]

--
arbitraryIntBounded :: Gen Int
arbitraryIntBounded = choose (2, 10)

arbitraryMatrix :: (Arbitrary e, Matrix m e) => Int -> Int -> Gen (m e)
arbitraryMatrix m n = do
    l <- vectorOf (m*n) arbitrary
    return $ fromList_m m n l

arbitraryTransMatrix :: (Arbitrary e, Matrix m e) => TransType a -> Int -> Int -> Gen (TransType (m e))
arbitraryTransMatrix (NoTrans _) m n = liftA NoTrans $ arbitraryMatrix m n
arbitraryTransMatrix (Trans _) m n = liftA Trans $ arbitraryMatrix n m
arbitraryTransMatrix (ConjTrans _) m n = liftA ConjTrans $ arbitraryMatrix n m

data GEMMTestData m e = GEMMTestData (TransType (m e)) (TransType (m e)) e e (m e) deriving (Show)
data StrassenTestData m e = StrassenTestData (TransType (m e)) (TransType (m e)) e e (m e) deriving (Show)
data TRSMTestData m e = TRSMTestData e (TransType (TriangType (UnitType (m e)))) (m e) deriving (Show)

instance (Matrix m e, Arbitrary e) => Arbitrary (GEMMTestData m e) where
    arbitrary = do
        m <- arbitraryIntBounded
        n <- arbitraryIntBounded
        k <- arbitraryIntBounded
        ttmA <- arbitrary
        ttmB <- arbitrary
        pmA <- arbitraryTransMatrix ttmA m k
        pmB <- arbitraryTransMatrix ttmB k n
        alpha <- arbitrary
        beta <- arbitrary
        mC <- arbitraryMatrix m n
        return $ GEMMTestData pmA pmB alpha beta mC

instance (Matrix m e, Arbitrary e) => Arbitrary (StrassenTestData m e) where
    arbitrary = do
        e <- choose (2, 4) :: Gen Int
        let n = truncate $ (2 :: Double) ** fromIntegral e
        ttmA <- arbitrary
        ttmB <- arbitrary
        alpha <- arbitrary
        mA <- arbitraryTransMatrix ttmA n n
        mB <- arbitraryTransMatrix ttmB n n
        beta <- arbitrary
        mC <- arbitraryMatrix n n
        return $ StrassenTestData mA mB alpha beta mC
--
instance (Matrix m e, Arbitrary e) => Arbitrary (TRSMTestData m e) where
    arbitrary = do
        m <- arbitraryIntBounded
        n <- arbitraryIntBounded
        trans <- arbitrary :: Gen (a -> TransType a)
        triang <- arbitrary :: Gen (a -> TriangType a)
        unit <- arbitrary :: Gen (a -> UnitType a)
        alpha <- arbitrary
        mA <- arbitraryMatrix m m
        mB <- arbitraryMatrix m n
        return $ TRSMTestData alpha (trans . triang . unit $ mA) mB

data PotrfTestData m e = PotrfTestData (TriangType (m e)) deriving (Show)

instance (Matrix m e, Elt e, Arbitrary e, Random e) => Arbitrary (PotrfTestData m e) where
    arbitrary = do
        _ <- arbitrary :: Gen (a -> TriangType a)
        n <- choose (3, 5)
        m <- arbSDPMatrix (blockSize ^ (n :: Int))
        return $ PotrfTestData $ Lower m

-- Cualquier propiedad tageada es una propiedad:
instance (Testable p) => Testable (Tagged t p) where
    property p = property $ unTagged p

data Tag2 (s :: *) (v :: * -> *)

instance (StratCtx HMatrixBindSeq ~ ctx_Seq) => Default (Tagged HMatrixBindSeq ctx_Seq) where
    def = Tagged newNullContext

instance (StratCtx dots ~ ctx_dots, Default (Tagged dots ctx_dots), StratCtx gs ~ ctx_gs, Default (Tagged gs ctx_gs), StratCtx (CholLLVSeq dots gs) ~ ctx_bp) => Default (Tagged (CholLLVSeq dots gs) ctx_bp) where
    def = Tagged (unTagged (def :: Tagged dots ctx_dots), unTagged (def :: Tagged gs ctx_gs))

instance (StratCtx dots ~ ctx_dots, Default (Tagged dots ctx_dots), StratCtx gs ~ ctx_gs, Default (Tagged gs ctx_gs), StratCtx (CholLLVPar_Repa dots gs) ~ ctx_bp) => Default (Tagged (CholLLVPar_Repa dots gs) ctx_bp) where
    def = Tagged (unTagged (def :: Tagged dots ctx_dots), unTagged (def :: Tagged gs ctx_gs))

instance (StratCtx ps ~ ctx_ps, Default (Tagged ps ctx_ps), StratCtx ss ~ ctx_ss, Default (Tagged ss ctx_ss), StratCtx gs ~ ctx_gs, Default (Tagged gs ctx_gs), StratCtx ts ~ ctx_ts, Default (Tagged ts ctx_ts), StratCtx (CholLLVBlocksSeq ps ss gs ts) ~ ctx_cs) => Default (Tagged (CholLLVBlocksSeq ps ss gs ts) ctx_cs) where
    def = Tagged (newSqrBlockContext blockSize, unTagged (def :: Tagged ps ctx_ps), unTagged (def :: Tagged ss ctx_ss), unTagged (def :: Tagged gs ctx_gs), unTagged (def :: Tagged ts ctx_ts))

--
instance (StratCtx DefSeq ~ ctx_DefSeq) => Default (Tagged DefSeq ctx_DefSeq) where
    def = Tagged newNullContext
instance (StratCtx CBindSeq ~ ctx_DefSeq) => Default (Tagged CBindSeq ctx_DefSeq) where
    def = Tagged newNullContext
instance (StratCtx gs ~ ctx_gs, Default (Tagged gs ctx_gs), StratCtx (BlocksSeq gs) ~ ctx_bp) => Default (Tagged (BlocksSeq gs) ctx_bp) where
    def = Tagged (newBlockContext 2 2, unTagged (def :: Tagged gs ctx_gs))
instance (StratCtx gs ~ ctx_gs, Default (Tagged gs ctx_gs), StratCtx ts ~ ctx_ts, Default (Tagged ts ctx_ts), StratCtx (BlocksByCols gs ts) ~ ctx_bbcp) => Default (Tagged (BlocksByCols gs ts) ctx_bbcp) where
    def = Tagged (newSqrBlockContext 2, unTagged (def :: Tagged gs ctx_gs), unTagged (def :: Tagged ts ctx_ts))
instance (StratCtx gs ~ ctx_gs, Default (Tagged gs ctx_gs), StratCtx ts ~ ctx_ts, Default (Tagged ts ctx_ts), StratCtx (BlocksByRows gs ts) ~ ctx_bbcp) => Default (Tagged (BlocksByRows gs ts) ctx_bbcp) where
    def = Tagged (newSqrBlockContext 2, unTagged (def :: Tagged gs ctx_gs), unTagged (def :: Tagged ts ctx_ts))

instance (StratCtx DefPar_MP ~ ctx_DefPar) => Default (Tagged DefPar_MP ctx_DefPar) where
    def = Tagged newNullContext
instance (StratCtx gs ~ ctx_gs, Default (Tagged gs ctx_gs), StratCtx (BlocksPar_MP gs) ~ ctx_bp) => Default (Tagged (BlocksPar_MP gs) ctx_bp) where
    def = Tagged (newBlockContext 2 2, unTagged (def :: Tagged gs ctx_gs))
instance (StratCtx ts ~ ctx_ts, Default (Tagged ts ctx_ts), StratCtx (ColsPar_MP ts) ~ ctx_cp) => Default (Tagged (ColsPar_MP ts) ctx_cp) where
    def = Tagged $ unTagged (def :: Tagged ts ctx_ts)
instance (StratCtx gs ~ ctx_gs, Default (Tagged gs ctx_gs), StratCtx (StrassenPar_MP gs) ~ ctx_bp) => Default (Tagged (StrassenPar_MP gs) ctx_bp) where
   def = Tagged (newStrassenContext 2, unTagged (def :: Tagged gs ctx_gs))

instance (StratCtx DefPar_ST ~ ctx_DefPar) => Default (Tagged DefPar_ST ctx_DefPar) where
    def = Tagged newNullContext
instance (StratCtx gs ~ ctx_gs, Default (Tagged gs ctx_gs), StratCtx (BlocksPar_ST gs) ~ ctx_bp) => Default (Tagged (BlocksPar_ST gs) ctx_bp) where
    def = Tagged (newBlockContext 2 2, unTagged (def :: Tagged gs ctx_gs))
instance (StratCtx ts ~ ctx_ts, Default (Tagged ts ctx_ts), StratCtx (ColsPar_ST ts) ~ ctx_cp) => Default (Tagged (ColsPar_ST ts) ctx_cp) where
    def = Tagged $ unTagged (def :: Tagged ts ctx_ts)
instance (StratCtx gs ~ ctx_gs, Default (Tagged gs ctx_gs), StratCtx (StrassenPar_ST gs) ~ ctx_bp) => Default (Tagged (StrassenPar_ST gs) ctx_bp) where
   def = Tagged (newStrassenContext 2, unTagged (def :: Tagged gs ctx_gs))

instance (StratCtx DefPar_Acc ~ ctx_DefPar) => Default (Tagged DefPar_Acc ctx_DefPar) where
    def = Tagged newNullContext

max_stars_repo_name: mauroblanco/fpnla-examples

max_issues_repo_name: mauroblanco/fpnla-examples

max_forks_repo_name: mauroblanco/fpnla-examples


text:
{-# options_ghc -Wwarn #-}
{-# language TypeOperators #-}
{-# language FlexibleInstances #-}
module Main where

import Solitaire.Yukon

import System.Directory
import System.Environment
import System.Microtimer

import qualified Data.Vector as V

import Statistics.Sample
import Data.List (sort)
import Options.Generic
import LinearRegression

data CliArgs w = CliArgs
  { numTrials  :: w ::: Int <!> "10"  <?> "Number of trials to run"
  , limitGames :: w ::: Int <!> "100" <?> "Limit number of games to test"
  , step       :: w ::: Float <!> "0.1" <?> "Step size"
  }
  deriving stock (Generic)

deriving instance Show (CliArgs Unwrapped)
instance ParseRecord (CliArgs Wrapped)

main :: IO ()
main = do
  args :: CliArgs Unwrapped <- unwrapRecord "Optimize heuristics"

  let
    limitGames = args ^. #limitGames
    numTrials = args ^. #numTrials

  fileNames <- take limitGames <$> listDirectory "games/Solvable"

  games <- for fileNames $ \name -> do
    fileContents <- readFile $ "games/Solvable/" <> name
    pure $ read fileContents

  let
    optimizeConfig = OptimizeConfig { numTrials, games }

    initialWeights = mempty
      & at "numFaceUp" ?~ 2
      & at "numFaceDown" ?~ 5
      & at "totalRunScore" ?~ (-0.5)

    test = testWeights optimizeConfig

  runTime <- test initialWeights

  prettyPrint runTime
  pure ()

data PartialDArgs = PartialDArgs
  { step :: Float
  , field :: Text
  , runTime :: RunTime
  , weights :: Weights
  , function :: Weights -> IO RunTime
  }
  deriving stock (Generic)

partialDerivative :: PartialDArgs -> IO RunTime
partialDerivative args = do
  let
    w_1 :: Weights
    w_1 = args ^. #weights & ix (args ^. #field) -~ (args ^. #step)

    w_2 :: Weights
    w_2 = args ^. #weights & ix (args ^. #field) +~ (args ^. #step)

    r0 = args ^. #runTime

  r1 <- args ^. #function $ w_1
  r2 <- args ^. #function $ w_2

  pure $ undefined

data Uncertain a = a :+- a

data Line = Line
  { m :: Uncertain Float
  , b :: Uncertain Float
  }
  deriving stock (Generic)

linearRegression :: [(Weight, RunTime)] -> Line
linearRegression = undefined

data OptimizeConfig = OptimizeConfig
  { numTrials :: !Int
  , games :: ![Game]
  }
  deriving stock (Generic)

data RunTime = RunTime
  { mean :: !Double
  , standardDeviation :: !Double
  }
  deriving stock (Generic)

instance Pretty RunTime where
  prettyExpr run =
    prettyExpr $ (Empty :: Map Text Text)
      & at "mean" ?~ (run ^. #mean . to formatSeconds . packed)
      & at "std dev" ?~ (run ^. #standardDeviation . to formatSeconds . packed)

testWeights :: OptimizeConfig -> Weights -> IO RunTime
testWeights config weights = do
  let
    games = config ^. #games
    numTrials = config ^. #numTrials

  times <- for games $ \game -> do
    let
      run = runGame @Yukon game
        & (runReaderT ?? weights)
        & time_

    sequenceA $ replicate numTrials run

  let
    sortedData = join times & sort & V.fromList
    (mean, variance) = meanVariance sortedData
    standardDeviation = sqrt variance

  pure $ RunTime { mean, standardDeviation }

updateWeightsAtRandom
  :: forall k m. (Ord k, MonadRandom m)
  => Weights
  -> m Weights
updateWeightsAtRandom map = do
  let
    n = length map
    keys = (map ^.. ifolded . asIndex)

  delta <- randomDelta n

  pure $ flip execState map $
    ifor_ keys $ \i key -> do
      ix key += (delta ^?! ix i)


text:
module Numeric.Sundials.ODEOpts where

import           Data.Word (Word32)
import qualified Data.Vector.Storable as VS

import           Numeric.LinearAlgebra.HMatrix (Vector, Matrix)

type Jacobian = Double -> Vector Double -> Matrix Double

data ODEOpts = ODEOpts {
    maxNumSteps :: Word32
  , minStep     :: Double
  , relTol      :: Double
  , absTols     :: VS.Vector Double
  , initStep    :: Maybe Double
  , maxFail     :: Word32
  } deriving (Read, Show, Eq, Ord)

data SundialsDiagnostics = SundialsDiagnostics {
    aRKodeGetNumSteps               :: Int
  , aRKodeGetNumStepAttempts        :: Int
  , aRKodeGetNumRhsEvals_fe         :: Int
  , aRKodeGetNumRhsEvals_fi         :: Int
  , aRKodeGetNumLinSolvSetups       :: Int
  , aRKodeGetNumErrTestFails        :: Int
  , aRKodeGetNumNonlinSolvIters     :: Int
  , aRKodeGetNumNonlinSolvConvFails :: Int
  , aRKDlsGetNumJacEvals            :: Int
  , aRKDlsGetNumRhsEvals            :: Int
  } deriving Show

max_stars_repo_name: steinitznavican/hmatrix-sundials

max_issues_repo_name: steinitznavican/hmatrix-sundials

max_forks_repo_name: steinitznavican/hmatrix-sundials


text:
-- Modulo destinado a todas as funcoes relativas ao treino da Rede.

module Training
(train) where

import Execution
import InputOutput
import Types
import Data.List.Split
import System.Random
import System.Random.Shuffle
import Numeric.LinearAlgebra.HMatrix
import Numeric.LinearAlgebra.Data

type Image = [Double]
type Sample = (Int, Image)

-- Funcao de treino que executa uma
-- dada quantidade de epocas de treino
train :: Int -> IO String
train epochAmount = do
                    trainingSet <- getTraining
                    testSet <- getTest
                    network <- readIn
                    manageTrainingEpoch epochAmount trainingSet testSet network

--
manageTrainingEpoch :: Int -> [Sample] -> [Sample] -> Data -> IO String
manageTrainingEpoch 0 _ _ _ = return ""
manageTrainingEpoch epochAmount trainingSet testSet network = do
                                                newNetwork <- trainingEpoch trainingSet network
                                                let
                                                    correctCnt = testEpoch testSet newNetwork
                                                    totalAmount = length testSet
                                                do  return $ printEpoch epochAmount correctCnt totalAmount
                                                    writeIn newNetwork
                                                    manageTrainingEpoch (epochAmount - 1) trainingSet testSet newNetwork

trainingEpoch :: [Sample] -> Data -> IO Data
trainingEpoch trainingSet network = do
                                standardGenerator <- getStdGen
                                let trainingSize = length trainingSet
                                    shuffledTrainingSet = shuffle' trainingSet trainingSize standardGenerator
                                    minibatchAmount = 20 -- adequar quantidade
                                    minibatchSize = trainingSize `div` 20 -- adequar quantidade
                                    minibatches = chunksOf minibatchSize shuffledTrainingSet
                                    newNetwork = manageMinibatch minibatchAmount 0 network minibatches
                                return newNetwork

manageMinibatch :: Int -> Int -> Data -> [[Sample]] -> Data
manageMinibatch amount counter network minibatches = if amount /= (counter + 1)
                                                        then
                                                          let changes = minibatchEvaluation (minibatches !! counter) amount network
                                                          in plus (plus network changes) (manageMinibatch amount (counter + 1) network minibatches)
                                                        else
                                                          generateBasedOf network

minibatchEvaluation :: [Sample] -> Int -> Data -> Data
minibatchEvaluation minibatch amount network = let sumedDesiredChanges = manageSample minibatch amount network
                                                   averageDesiredChanges = divide sumedDesiredChanges (fromIntegral amount)
                                               in averageDesiredChanges

manageSample :: [Sample] -> Int -> Data -> Data
manageSample minibatch counter networkModel = if counter > 0
                                    then
                                        let representedInt = fst $ minibatch !! counter
                                            image = snd $ minibatch !! counter
                                            network = feedforward image networkModel
                                            expectedOutput = buildExpectedOutput representedInt
                                            desiredChanges = backpropagation network image expectedOutput
                                            sumChanges = generateBasedOf networkModel
                                        in plus (plus sumChanges desiredChanges) (manageSample minibatch counter networkModel)
                                    else
                                        generateBasedOf networkModel

buildExpectedOutput :: Int -> [Double]
buildExpectedOutput representedValue = let indexes = [0.0 .. 9.0]
                                       in [if x == (fromIntegral representedValue) then x else 0.0 | x <- indexes]

-- Recebe as informacoes da rede neural, o resultado
-- esperado e retorna um Data com as modificacoes necessarias
-- na rede
-- N = network, E = expected list, I = image values
backpropagation :: Data -> Image -> [Double] -> Data
backpropagation n i e
    | isEmpty n = error "Backpropagation error: Data is empty"
    | length e /= 10 = error "Backpropagation error: expectedOutput list is invalid"
    | otherwise = let outputError = computeOutputError (aOutput n) (fromList e) (zetaOutput n)
                      hiddenError = computeHiddenError (wOutput n) outputError (zetaHidden n)
                      outputDesiredChanges = computeOutputDesiredChanges outputError (aHidden n)
                      hiddenDesiredChanges = computeHiddenDesiredChanges hiddenError i
                      wH = fst hiddenDesiredChanges
                      bH = snd hiddenDesiredChanges
                      wO = fst outputDesiredChanges
                      bO = snd outputDesiredChanges
                      newNetwork = (Data wH bH (aHidden n) (zetaHidden n) wO bO (aOutput n) (zetaOutput n))
                  in newNetwork

-- Gera o vetor de erro da camada output, recebe
-- as ativacoes do output atual, as ativacoes esperadas
-- e a lista zeta do output
computeOutputError :: Vector Double -> Vector Double -> Vector Double -> Vector Double
computeOutputError a e z = fromList $ hadamardV (toList (add a (scale (-1) e))) (sigV' $ toList z)

computeHiddenError :: Matrix R -> Vector Double -> Vector Double -> Vector Double
computeHiddenError ow oe hz = fromList $ hadamardV (toList $ (tr' ow) #> oe) (sigV' (toList hz))

computeOutputDesiredChanges :: Vector Double -> Vector Double -> (Matrix R, Vector Double)
computeOutputDesiredChanges oe ah = let owDesired = oe `outer` ah
                                        ohDesired = oe
                                    in (owDesired, ohDesired)

computeHiddenDesiredChanges :: Vector Double -> Image -> (Matrix R, Vector Double)
computeHiddenDesiredChanges he image = let hwDesired = he `outer` (fromList image)
                                           hbDesired = he
                                       in (hwDesired, hbDesired)

testEpoch :: [Sample] -> Data -> Int
testEpoch testSet network = manageEpoch testSet network (length testSet)

manageEpoch :: [Sample] -> Data -> Int -> Int
manageEpoch testSet network counter = if counter > 0
                                        then
                                            let newNetwork = feedforward (snd $ testSet !! counter) network
                                            in if (toList $ aOutput newNetwork) == buildExpectedOutput (fst $ testSet !! counter)
                                                then
                                                    1 + manageEpoch testSet network (counter - 1)
                                                else
                                                    0 + manageEpoch testSet network (counter - 1)
                                        else
                                            0


text:
{-# LANGUAGE GeneralisedNewtypeDeriving #-}
-- |
-- Module Orbital
--
-- A probability distribution for driver laptimes, initialised from
-- Elo-like ratings.
--
-- In the PDF, t = 0 is taken to be the ideal laptime in a track.
--
-- The name "orbital" alludes to this PDF originally having been, modulo
-- constant factors, the radial probability density of the 1s hydrogen atom
-- orbital. The corresponding distribution can be obtained by setting the
-- gamma shape parameter to 3.
module Orbital
    ( deltaWP
    , eloWP
    , initialRating
    , kFromRating
    , ratingFromK
    , OrbitalDistribution
    , orbitalDistr
    , orbitalK
    ) where

import Numeric.SpecFunctions (choose)
import Statistics.Distribution
import Statistics.Distribution.Gamma
import Data.MemoTrie
import Data.List (foldl')

-- | The base implementation of victory probabilities for performances
-- modelled by gamma distributions with the same (integer) shape.
ratioWP :: Int -> Double -> Double
ratioWP gsh = \w -> w^gsh * evalPoly (memoWpCoefficient gsh) (gsh-1) w

-- | A strict pair type.
data SP a b = SP !a !b

-- | A polynomial evaluator. For our immediate purposes, it performs slightly
-- better than the one available in 'Numeric.Polynomial'.
evalPoly :: (Int -> Double) -> Int -> Double -> Double
evalPoly coef n x = y
    where
    SP xn y = foldl' alg (SP 1 0) [0..n]
    alg (SP xx s) m = SP (x * xx) (coef m * xx + s)

-- | Coefficients of the polynomial used in the gamma winning probability
-- calculation.
wpCoefficient :: Int -> Int -> Double
wpCoefficient gsh m =
    (-1)^m * choose gsh m * choose (2*gsh-1) (gsh-1)
        * fromIntegral (gsh-m) / fromIntegral (gsh+m)

-- | Memoised version of 'wpCoefficient'.
memoWpCoefficient :: Int -> Int -> Double
memoWpCoefficient = memo2 wpCoefficient

-- | Initial rating for new racers. Defined as a constant here for
-- expediteness, taking into account that making it configurable isn't as
-- essential as it is making other engine parameters configurable.
initialRating :: Double
initialRating = 1500

-- | The base Elo exponent. In essence, a logistic growth factor for expected
-- scores/winning probabilities. The denominator used here, 400, is standard
-- for chess rating computations. Using it, a 0.6 expected score corresponds
-- to a ~70.4 rating gap, and a 0.9 expected score, to a ~382 gap.
eloAlpha :: Double
eloAlpha = log 10 / 400

-- | Adjusts the k values obtained from ratings depending on the gamma shape so
-- that the same rating gaps correspond to approximately the same winning
-- probabilities.
--
-- The formula used here is made of two factors:
--
-- - A base correction, which on its own would adjust the derivative of
--   @deltaWP gsh@ to match that of @deltaWP 1@ at zero. This factor, which is
--   the best adjustment for rating gaps close to zero, slightly increases
--   the absolute excess winning probabilities over the entire range of gaps.
--
-- - An additional correction, which reduces the aforementioned bias. Its value
--   is chosen to (approximately) minimise the sum of the squares of the
--   differences between @ratioWP gsh@ and @ratioWP 1@ over the entire range
--   of ratios.
eloGammaCorrection :: Int -> Double
eloGammaCorrection gsh = extraAdjustment * baseCorrection
    where
    gsh' = fromIntegral gsh
    baseCorrection = 4^(gsh-1) / (gsh' * choose (2*gsh-1) (gsh-1))
    fittedConstant = (2*sqrt 29 - 3)*pi
    extraAdjustment = 1 - (gsh'-1) / (gsh' * fittedConstant)

-- | Elo conversion factor. Amounts to eloAlpha for shape 1, in which case the
-- gamma model coincides with the conventional Elo system.
eloFactor :: Int -> Double
eloFactor gsh = eloGammaCorrection gsh * eloAlpha

-- | Memoised Elo conversion factor.
memoEloFactor :: Int -> Double
memoEloFactor = memo eloFactor

-- | Winning probabilities for (scaled) rating differences. Ratings are
-- essentially logarithms of the gamma distribution rate parameters.
deltaWP :: Int -> Double -> Double
deltaWP gsh = \d -> ratioWP gsh (1 / (1 + exp(- memoEloFactor gsh * d)))
-- If ru = log u, rv = log v, and d = ru - rv, then
-- u / (u + v) = 1 / (1 + exp (-d))
-- Another way of expressing the ratio is given by
-- 1 / (1 + exp(-d)) = (1 + tanh (d/2)) / 2

-- | Elo-based winning probability. This version takes the ratings of the two
-- racers. Only here for reference, at least for now.
eloWP :: Double -> Double -> Double
eloWP ru rv = deltaWP 1 (ru-rv)

-- | Rating that corresponds to a gamma rate parameter of 1.
--
-- As far as the results go, the effect of 'referenceRating' is superficial.
-- Since it adjusts all ratings in the same manner, changing it does not affect
-- the overall results.
--
-- There used to be a counterpart @referenceK@ definition, which was removed
-- for being redundant. It is worth noting that scaling the conversions to make
-- the gamma rate ("k") values larger slows down the calculations, probably by
-- making it harder for the gamma distribution algorithms.
referenceRating :: Double
referenceRating = initialRating

-- | k-to-rating conversion. Only for reference, at least for now.
ratingFromK :: Int -> Double -> Double
ratingFromK gsh u = referenceRating + log u / eloFactor gsh

-- | rating-to-k conversion.
kFromRating :: Int -> Double -> Double
kFromRating gsh r = exp ((r - referenceRating) * eloFactor gsh)
-- Using memoEloFactor instead of eloFactor here doesn't seem to improve
-- performance.

-- | The model distribution is a special case of the gamma distribution, with
-- theta = 1/k.
newtype OrbitalDistribution = OrbitalDistribution
    { orbitalGamma :: GammaDistribution }
    deriving (Eq, Show, Distribution, ContDistr, ContGen, MaybeEntropy
        , Variance, MaybeVariance, Mean, MaybeMean)

-- | Sets up the model distribution given a rating.
orbitalDistr
    :: Int    -- ^ Shape parameter to be used.
    -> Double -- ^ Elo-like rating.
    -> OrbitalDistribution
orbitalDistr gsh r = OrbitalDistribution
    (gammaDistr (fromIntegral gsh) (1 / kFromRating gsh r))

-- | Recovers the k factor from the distribution.
orbitalK :: OrbitalDistribution -> Double
orbitalK distr = 1 / gdScale (orbitalGamma distr)

-- | Recovers the shape parameter from the distribution.
orbitalShape :: OrbitalDistribution -> Double
orbitalShape distr = gdShape (orbitalGamma distr)

max_stars_repo_name: duplode/elo-zs

max_issues_repo_name: duplode/elo-zs

max_forks_repo_name: duplode/elo-zs