text:
-- This file is part of Quipper. Copyright (C) 2011-2014. Please see the
-- file COPYRIGHT for a list of authors, copyright holders, licensing,
-- and other details. All rights reserved.
--
-- ======================================================================

{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE UndecidableInstances #-}

-- ----------------------------------------------------------------------
-- | This program reads an execution trace produced by QCL, and turns
-- it into a Quipper circuit.

module Programs.QCLParser.Main where

import Quipper hiding (cnot, Format)

import Data.IntMap (IntMap)
import qualified Data.IntMap as IntMap
import qualified Data.Traversable as Trav
import Control.Monad.State
import Prelude hiding (not)
import qualified Prelude
import System.IO
import Data.Complex
import Data.Maybe

import Text.ParserCombinators.ReadP hiding (get)
import Data.Char (isAlpha, isAlphaNum, isDigit)

-- ----------------------------------------------------------------------
-- * A monad for a QCL state

-- | In QCL, qubits are identified by integers. We have to map these
-- to Quipper's native qubits. A 'QCLState' holds such a map.
-- Implicitly, it also holds the set of qubits currently defined.
type QCLState = IntMap Qubit

-- | The 'QCLCirc' monad is like the 'Circ' monad, except that it also
-- keeps track of an additional 'QCLState'. The 'lift' function must
-- be used to lift any command for the 'Circ' monad to the 'QCLCirc'
-- monad.
type QCLCirc a = StateT QCLState Circ a

-- ----------------------------------------------------------------------
-- * Auxiliary state functions

-- | Look up the qubit corresponding to a QCL register, or allocate a
-- new qubit if it doesn't already exist.
provide :: Int -> QCLCirc Qubit
provide r = do
  s <- get
  case IntMap.lookup r s of
    Just q -> return q
    Nothing -> do
      q <- lift $ qinit False
      let s' = IntMap.insert r q s
      put s'
      return q

-- ----------------------------------------------------------------------
-- * Implementation of the QCL primitives

-- | Reset all qubits to state 0.
qcl_reset :: QCLCirc ()
qcl_reset = do
  s <- get
  lift $ Trav.mapM qdiscard s
  let s' = IntMap.empty
  put s'
  return ()

-- | Apply a controlled-not operation to the first argument.
qcl_cnot :: Int -> [Int] -> QCLCirc ()
qcl_cnot r ctrls = do
  q <- provide r
  cs <- Trav.mapM provide ctrls
  lift $ qnot_at q `controlled` cs
  return ()

-- | Apply an uncontrolled not operation.
qcl_not :: Int -> QCLCirc ()
qcl_not r = qcl_cnot r []

-- | @'qcl_fanout' ins outs ctrls@: Copy the qubit register /ins/ to
-- the qubit register /outs/ by means of c-not operations, provided
-- that /outs/ was previously 0. The whole operation is controlled by
-- /ctrls/.
qcl_fanout :: [Int] -> [Int] -> [Int] -> QCLCirc ()
qcl_fanout ins outs ctrls = do
  qins <- Trav.mapM provide ins
  qouts <- Trav.mapM provide outs
  qctrls <- Trav.mapM provide ctrls
  lift $ with_controls qctrls $ do
    let zips = zip qins qouts
    Trav.mapM (\(x,y) -> qnot_at y `controlled` x) zips
  return ()

-- | Calculate the square distance between two vectors, which must be
-- of the same length.
sqdist :: [Complex Double] -> [Complex Double] -> Double
sqdist v1 v2 = d2 where
  w = zipWith (-) v1 v2
  w2 = map (\x -> x * x) $ map magnitude w
  d2 = sum w2

-- | If the matrix looks like a /W/-gate, return 'True'.
matrix_w :: [Complex Double] -> Bool
matrix_w amps | length amps == 16 =
  let v = sqrt 0.5
      m = [ 1, 0, 0, 0,
            0, v, v, 0,
            0, v, -v, 0,
            0, 0, 0, 1]
  in
   if sqdist amps m <= 0.001 then
     True
   else
     False
matrix_w _ = False

-- | If the matrix looks like a controlled /e/^/tiZ/-gate, return the
-- angle /t/.
matrix_exp :: [Complex Double] -> Maybe Double
matrix_exp amps | length amps == 16 =
  let t = phase ((amps !! 10) + conjugate (amps !! 15))
      u = cis t
      m = [ 1, 0, 0, 0,
            0, 1, 0, 0,
            0, 0, u, 0,
            0, 0, 0, conjugate u]
  in
   if sqdist amps m <= 0.001 then
     Just t
   else
     Nothing
matrix_exp _ = Nothing

-- | @'qcl_matrix' n amps regs@: Apply the /n/-by-/n/ unitary gate
-- whose matrix is given in /amps/, to the qubit list /regs/.  This
-- function must guess, based on the complex entries of the matrix,
-- what the name of the gate should be. This guessing is crude at the
-- moment, and must be extended to include additional gates as
-- required by each algorithm. If the first argument is 'True', invert
-- the matrix.
qcl_matrix :: Bool -> Int -> [Complex Double] -> [Int] -> QCLCirc ()
qcl_matrix inv n amps [x,y] | matrix_w amps = do
  [q,r] <- Trav.mapM provide [x,y]
  lift $ gate_W_at q r  -- this gate is self-inverse
qcl_matrix inv n amps [x,y] | isJust t = do
  [q,r] <- Trav.mapM provide [x,y]
  lift $ expZt_at (sign * fromJust t) r `controlled` q
    where
      t = matrix_exp amps
      sign = if inv then -1 else 1
qcl_matrix inv n amps regs = do
  qregs <- Trav.mapM provide regs
  lift $ named_gate_at (if inv then "Gate*" else "Gate") qregs

-- | The inverse of 'qcl_cnot'.
qcl_cnot_inv = qcl_cnot

-- | The inverse of 'qcl_not'.
qcl_not_inv = qcl_not

-- | The inverse of 'qcl_fanout'.
qcl_fanout_inv = qcl_fanout

-- | A sample circuit to illustrate how to use the primitives.
testcircuit1 :: QCLCirc ()
testcircuit1 = do
  qcl_reset
  qcl_not 5
  qcl_not 0
  qcl_cnot 31 [4]
  qcl_cnot 31 [4,30]
  qcl_cnot 29 [31]
  qcl_cnot_inv 31 [4,30]
  qcl_cnot_inv 31 [4]

-- ----------------------------------------------------------------------
-- * Unpacking QCLCirc

-- | Run function for the 'QCLCirc' monad: execute the actions and
-- produce a circuit.
run :: QCLCirc a -> Circ a
run f = do
  (x,_) <- runStateT f IntMap.empty
  return x

-- ----------------------------------------------------------------------
-- * An abstract syntax for QCL output

-- | A data type to hold a QCL gate.
data QCLGate =
  Comment
  | Reset
  | Not Int
  | XNot Int
  | CNot Int [Int]
  | XCNot Int [Int]
  | Fanout [Int] [Int] [Int]
  | XFanout [Int] [Int] [Int]
  | Matrix Int [Complex Double] [Int]
  | XMatrix Int [Complex Double] [Int]
  deriving (Show)

-- | Take a gate from the abstract syntax and execute it in the
-- 'QCLCirc' monad.
do_qcl_gate :: QCLGate -> QCLCirc ()
do_qcl_gate Comment = return ()
do_qcl_gate Reset = qcl_reset
do_qcl_gate (Not x) = qcl_not x
do_qcl_gate (XNot x ) = qcl_not_inv x
do_qcl_gate (CNot x ctrl) = qcl_cnot x ctrl
do_qcl_gate (XCNot x ctrl) = qcl_cnot_inv x ctrl
do_qcl_gate (Fanout a b c) = qcl_fanout a b c
do_qcl_gate (XFanout a b c) = qcl_fanout_inv a b c
do_qcl_gate (Matrix n amps regs) = qcl_matrix False n amps regs
do_qcl_gate (XMatrix n amps regs) = qcl_matrix True n amps regs

-- ----------------------------------------------------------------------
-- * Parsing

-- $ The output of QCL consists of lines of the following forms. Lines
-- not starting with \"\@\" are comments or other non-circuit output.
--
-- > @ RESET
-- > @ NOT(qureg q=<5>)
-- > @ CNOT(qureg q=<31>,quconst c=<4,30>)
-- > @ !CNOT(qureg q=<31>,quconst c=<4,30>)
-- > @ Fanout(quconst a=<47,48,49,50,51,52>,quvoid b=<40,41,42,43,44,45>;cond=<>)
-- > @ !Fanout(quconst a=<47,48,49,50,51,52>,quvoid b=<40,41,42,43,44,45>;cond=<>)
-- > @ Matrix4x4(complex u00=1,complex u01=0,complex u02=0,complex u03=0,complex u10=0,complex u11=1,complex u12=0,complex u13=0,complex u20=0,complex u21=0,complex u22=(0.995004,-0.0998334),complex u23=0,complex u30=0,complex u31=0,complex u32=0,complex u33=(0.995004,0.0998334),qureg q=<12,13>)
-- > @ !Matrix4x4(complex u00=1,complex u01=0,complex u02=0,complex u03=0,complex u10=0,complex u11=0.707107,complex u12=0.707107,complex u13=0,complex u20=0,complex u21=0.707107,complex u22=-0.707107,complex u23=0,complex u30=0,complex u31=0,complex u32=0,complex u33=1,qureg q=<0,6>)
--
-- We use Koen Claessen's parser combinators (see
-- "Text.ParserCombinators.ReadP") to implement the parser.

-- | Parse a QCL identifier, which we take to be a non-empty string of
-- alphanumeric characters, starting with a letter
identifier :: ReadP String
identifier = do
  satisfy isAlpha
  munch isAlphaNum

-- | Parse a signless integer. We avoid the usual trick ('readS_to_P'
-- 'reads'), because this introduces backtracking errors.
int :: ReadP Int
int = do
  s <- munch1 isDigit
  return $ (read s :: Int)

-- | Parse a floating point number. We avoid the usual trick
-- ('readS_to_P' 'reads'), because this introduces backtracking
-- errors.
double :: ReadP Double
double = do
  (s, _) <- gather parse_double
  return $ (read s :: Double)
  where
    parse_double = do
      option '+' (char '+' +++ char '-')
      munch isDigit
      optional (char '.' >> munch1 isDigit)

-- | Parse a comma separated list.
commalist :: ReadP a -> ReadP [a]
commalist elt = sepBy elt (skipSpaces >> char ',' >> skipSpaces)

-- | Parse a QCL quantum register of the form
--
-- > q=<31,32>
-- > c=<4,31>
-- > b=<40,41,42,43,44,45>.
qureg :: ReadP (String, [Int])
qureg = do
  id <- identifier
  skipSpaces
  char '='
  skipSpaces
  rs <- between (char '<' >> skipSpaces) (skipSpaces >> char '>') $ do
    commalist int
  return (id, rs)

-- | Consume an optional \"!\". Return 'True' if consumed, and 'False'
-- otherwise.
inversechar :: ReadP Bool
inversechar = do
  c <- option '+' (char '!')
  return (c == '!')

-- | Parse a complex number declaration of the format
--
-- > complex u00=1
--
-- or
--
-- > complex u22=(0.995004,-0.0998334).
complex :: ReadP (Complex Double)
complex = do
  string "complex"
  skipSpaces
  identifier
  skipSpaces
  char '='
  skipSpaces
  choice [
    do
      x <- double
      return (x :+ 0)
    ,
    do
      char '('
      skipSpaces
      x <- double
      skipSpaces
      char ','
      skipSpaces
      y <- double
      skipSpaces
      char ')'
      return (x :+ y)
    ]

-- | Parse a single line of QCL output into a 'QCLGate'.
qcl_line :: ReadP QCLGate
qcl_line = choice [
  do -- @ RESET
    char '@'
    skipSpaces
    string "RESET"
    return Reset
  ,
  do -- @ NOT(qureg q=<5>)
    char '@'
    skipSpaces
    inv <- inversechar
    skipSpaces
    string "NOT"
    skipSpaces
    char '('
    skipSpaces
    string "qureg"
    skipSpaces
    (_,[r]) <- qureg
    skipSpaces
    char ')'
    return (if inv then XNot r else Not r)
  ,
  do -- @ CNOT(qureg q=<31>,quconst c=<4,30>)
    char '@'
    skipSpaces
    inv <- inversechar
    skipSpaces
    string "CNOT"
    skipSpaces
    char '('
    skipSpaces
    string "qureg"
    skipSpaces
    (_,[r]) <- qureg
    skipSpaces
    char ','
    skipSpaces
    string "quconst"
    skipSpaces
    (_,ctrls) <- qureg
    skipSpaces
    char ')'
    return (if inv then XCNot r ctrls else CNot r ctrls)
  ,
  do -- @ Fanout(quconst a=<47,48>,quvoid b=<40,41>;cond=<>)
    char '@'
    skipSpaces
    inv <- inversechar
    skipSpaces
    string "Fanout"
    skipSpaces
    char '('
    skipSpaces
    string "quconst"
    skipSpaces
    (_,ins) <- qureg
    skipSpaces
    char ','
    skipSpaces
    string "quvoid"
    skipSpaces
    (_,outs) <- qureg
    skipSpaces
    char ';'
    skipSpaces
    (_,ctrls) <- qureg
    skipSpaces
    char ')'
    return (if inv then XFanout ins outs ctrls else Fanout ins outs ctrls)
  ,
  do -- @ Matrix4x4(complex u00=1,complex u01=0,complex u02=0,complex u03=0,complex u10=0,complex u11=1,complex u12=0,complex u13=0,complex u20=0,complex u21=0,complex u22=(0.995004,-0.0998334),complex u23=0,complex u30=0,complex u31=0,complex u32=0,complex u33=(0.995004,0.0998334),qureg q=<12,13>)
    char '@'
    skipSpaces
    inv <- inversechar
    skipSpaces
    string "Matrix"
    skipSpaces
    dim1 <- int
    skipSpaces
    char 'x'
    skipSpaces
    dim2 <- int
    skipSpaces
    char '('
    skipSpaces
    amps <- commalist complex
    skipSpaces
    char ','
    skipSpaces
    string "qureg"
    skipSpaces
    (_,r) <- qureg
    skipSpaces
    char ')'
    when (dim1 /= dim2) $ do
      error "Non-square matrix"
    return (if inv then XMatrix dim1 amps r else Matrix dim1 amps r)
  ,
  do -- any line not starting with '@' is a comment
    satisfy ((/=) '@')
    munch (\x -> True)
    return Comment
  ,
  do -- empty lines are comments
    eof
    return Comment
  ]

-- | String version of 'qcl_line': parse a string and turn it into a
-- 'QCLGate'.
parse_qcl_line :: String -> QCLGate
parse_qcl_line s =
  case readP_to_S readline s of
    (h, _):_ -> h
    _ -> error ("Unrecognized line: " ++ s)
  where
    readline = do
      skipSpaces
      l <- qcl_line
      skipSpaces
      eof
      return l

-- | Monad version of 'parse_qcl_line': parse a string and execute the
-- resulting gate directly in the 'QCLCirc' monad.
run_qcl_line :: String -> QCLCirc ()
run_qcl_line = do_qcl_gate . parse_qcl_line

-- | Parse a stream consisting of many lines of QCL output and execute
-- the parsed gates in the 'QCLCirc' monad.
run_qcl_lines :: String -> QCLCirc ()
run_qcl_lines s =
  mapM_ run_qcl_line (lines s)

-- | A sample circuit to illustrate the parser.
testcircuit2 :: QCLCirc ()
testcircuit2 = do
  run_qcl_line "@ RESET"
  run_qcl_line "@ NOT(qureg q=<5>)"
  run_qcl_line "@ NOT(qureg q=<0>)"
  run_qcl_line "@ CNOT(qureg q=<31>,quconst c=<4>)"

-- ----------------------------------------------------------------------
-- * Main function

-- | Main function: read a circuit in QCL format from 'stdin', and
-- preview the translated Quipper circuit.
main :: IO ()
main = do
  lines <- hGetContents stdin
  print_simple Preview (run (run_qcl_lines lines))

max_stars_repo_name: fritzo/quipper

max_issues_repo_name: fritzo/quipper

max_forks_repo_name: fritzo/quipper


text:
module Common where

import Numeric
import Data.Ratio
import Data.Complex
import Data.Array
import Text.ParserCombinators.Parsec hiding (spaces)
import Data.IORef
import Control.Monad.Except
import System.IO

data LispVal
    = Atom String
    | List [LispVal]
    | DottedList [LispVal] LispVal
    | Number Integer
    | String String
    | Bool Bool
    | Character Char
    | Float Double
    | Ratio Rational
    | Complex (Complex Double)
    | Vector (Array Int LispVal)
    | PrimitiveFunc ([LispVal] -> ThrowsError LispVal)
    | Func {params :: [String], vararg :: (Maybe String),
            body :: [LispVal], closure :: Env}
    | IOFunc ([LispVal] -> IOThrowsError LispVal)
    | Port Handle
-- Constructors and types have different namespaces, so you can have both a constructor named String and a type named String. Both types and constructor tags always begin with capital letters.
instance Eq LispVal where
    (Atom x) == (Atom y) = x == y
    (List xs) == (List ys) = xs == ys
    (DottedList xs x) == (DottedList ys y) = xs == ys && x == y
    (Number x) == (Number y) = x == y
    (String x) == (String y) = x == y
    (Bool x) == (Bool y) = x == y
    (Character x) == (Character y) = x == y
    (Float x) == (Float y) = x == y
    (Ratio x) == (Ratio y) = x == y
    (Complex x) == (Complex y) = x == y
    (Vector xa) == (Vector ya) =  xa == ya
    (PrimitiveFunc x) == (PrimitiveFunc y) = False
    (Func {params=paramsx, vararg=varargx, body=bodyx, closure=closurex}) == (Func {params=paramsy, vararg=varargy, body=bodyy, closure=closurey}) = paramsx == paramsy && varargx == varargy && bodyx == bodyy && closurex == closurey

unwordsList :: [LispVal] -> String
unwordsList = unwords . map showVal

showVal :: LispVal -> String
showVal (String contents) = "\"" ++ contents ++ "\""
showVal (Atom name) = name
showVal (Number contents) = show contents
showVal (Bool True) = "#t"
showVal (Bool False) = "#f"
showVal (List contents) = "(" ++ unwordsList contents ++ ")"
showVal (DottedList head tail) = "(" ++ unwordsList head ++ " . " ++ showVal tail ++ ")"
showVal (PrimitiveFunc _) = "<primitive>"
showVal (Func {params = args, vararg = varargs, body = body, closure = env}) =
    "(lambda (" ++ unwords (map show args) ++
        (case varargs of
            Nothing -> ""
            Just arg -> " . " ++ arg) ++ ") ...)"
showVal (Port _) = "<IO port>"
showVal (IOFunc _) = "<IO primitive>"

instance Show LispVal where show = showVal

data LispError  = NumArgs Integer [LispVal]
                | TypeMismatch String LispVal
                | Parser ParseError
                | BadSpecialForm String LispVal
                | NotFunction String String
                | UnboundVar String String
                | Default String

type ThrowsError = Either LispError
type Env = IORef [(String, IORef LispVal)]
type IOThrowsError = ExceptT LispError IO

showError :: LispError -> String
showError (UnboundVar message varname)  = message ++ ": " ++ varname
showError (BadSpecialForm message form) = message ++ ": " ++ show form
showError (NotFunction message func)    = message ++ ": " ++ show func
showError (NumArgs expected found)      = "Expected" ++ show expected ++ "args: found values " ++ unwordsList found
showError (TypeMismatch expected found) = "Invalid type: expected " ++ expected ++ ", found " ++ show found
showError (Parser parseErr)             = "Parse error at " ++ show parseErr

instance Show LispError where show = showError

ext: hs


text:
module Ch20 where

import Data.Complex

dft :: RealFloat a => [Complex a] -> [Complex a]
dft xs =
    let lenI = length xs
        lenR = fromIntegral lenI
        lenC = lenR :+ 0
    in [let i = negate 2 * pi * fromIntegral k/lenR
        in (1/lenC) * sum [(xs !! n) * exp (0 :+ i * fromIntegral n)
            | n <- [0, 1 .. lenI - 1]]
        | k <- [0, 1 .. lenI - 1]]

-- This seems to work, but I don't truly understand it
idft :: RealFloat a => [Complex a] -> [Complex a]
idft xs =
    let lenI = length xs
        lenR = fromIntegral lenI
        lenC = lenR :+ 0
    in [let i = 2 * pi * fromIntegral k/lenR
        in sum [(xs !! n) * exp (0 :+ i * fromIntegral n)
            | n <- [0, 1 .. lenI - 1]]
        | k <- [0, 1 .. lenI - 1]]

printComplexL :: [Complex Double] -> IO ()
printComplexL xs =
    let f (i, rl :+ im) =
            do putStr (spaces (3 - length (show i)))
               putStr (show i ++ ": (" )
               putStr (niceNum rl ++ ", " )
               putStr (niceNum im ++ ")\n" )
    in mapM_ f (zip [0 .. length xs - 1] xs)

niceNum :: Double -> String
niceNum d =
    let d' = fromIntegral (round (1e10 * d))/1e10
        (dec, fra) = break (== '.') (show d')
        (fra', exp) = break (== 'e') fra
    in spaces (3 - length dec) ++ dec ++ take 11 fra'
        ++ exp ++ spaces (12 - length fra' - length exp)

spaces :: Int -> String
spaces n = take n (repeat ' ')

mkTerm :: Int -> Double -> [Complex Double]
mkTerm num n = let f = 2 * pi/fromIntegral num
    in [sin (n * f * fromIntegral i)/n :+ 0
        | i <- [0, 1 .. num - 1]]

mkxa, mkxb, mkxc :: Int -> [Complex Double]
mkxa num = mkTerm num 1
mkxb num = zipWith (+) (mkxa num) (mkTerm num 3)
mkxc num = zipWith (+) (mkxb num) (mkTerm num 5)

testXc n = printComplexL (dft (mkxc n))

mkPulse :: Int -> [Complex Double]
mkPulse n = 100 : take (n - 1) (repeat 0)

testPulse n = printComplexL (dft (mkPulse n))

max_stars_repo_name: sackyhack/HSoM-Exercises

max_issues_repo_name: sackyhack/HSoM-Exercises

max_forks_repo_name: sackyhack/HSoM-Exercises


text:
module CommonTypes
  (LispVal(..),
   LispError(..),
   ThrowsError,
   trapError,
   extractValue,
   Env,
   IOThrowsError,
   nullEnv,
   liftThrows,
   runIOThrows)
  where

import qualified Data.Vector as V
import Data.Ratio
import Data.Complex
import Control.Monad.Except
import Text.ParserCombinators.Parsec
import Data.IORef
import System.IO

data LispVal  = Atom String
              | List [LispVal]
              | DottedList [LispVal] LispVal
              | Vector (V.Vector LispVal)
              | Number Integer
              | String String
              | Bool Bool
              | Character Char
              | Float Double
              | Ratio Rational
              | Complex (Complex Double)
              | PrimitiveFunc ([LispVal] -> ThrowsError LispVal)
              | Func { params :: [String], vararg :: (Maybe String),
                     body :: [LispVal], closure :: Env }
              | IOFunc ([LispVal] -> IOThrowsError LispVal)
              | Port Handle

instance Show LispVal where show = showVal

showVal :: LispVal -> String
showVal (String string) = "\"" ++ string ++ "\""
showVal (Atom name) = name
showVal (Number num) = show num
showVal (Bool value) = if value then "#t" else "#f"
showVal (List contents) = "(" ++ unwordsList contents ++ ")"
showVal (DottedList head tail) = "(" ++ unwordsList head ++ " . " ++ showVal tail ++ ")"
showVal (Vector vector) = "#(" ++ (unwordsList $ V.toList vector) ++ ")"
showVal (Float value) = show value
showVal (Ratio value) = show (numerator value) ++ "/" ++ show (denominator value)
showVal (Complex value) = show (realPart value) ++ "+" ++ show (imagPart value) ++ "i"
showVal (Character char) | char == '\n' = "#\\newline"
                         | char == ' '  = "#\\space"
                         | otherwise    = "#\\" ++ [char]
showVal (PrimitiveFunc _) = "<primitive>"
showVal (Func {params = args, vararg = varargs, body = body, closure = env}) =
     "(lambda (" ++ unwords (map show args) ++
             (case varargs of
                          Nothing -> ""
                          Just arg -> " . " ++ arg) ++ ") ...)"
showVal (Port _) = "<IO port>"
showVal (IOFunc _) = "<IO primitve>"

unwordsList = unwords . map showVal

data LispError = NumArgs Integer [LispVal]
               | TypeMismatch String LispVal
               | Parser ParseError
               | BadSpecialForm String LispVal
               | NotFunction String String
               | UnboundVar String String
               | Default String

instance Show LispError where show = showError

type ThrowsError = Either LispError

showError :: LispError -> String
showError (UnboundVar message varname)  = message ++ ": " ++ varname
showError (BadSpecialForm message form) = message ++ ": " ++ show form
showError (NotFunction message func)    = message ++ ": " ++ show func
showError (NumArgs expected found)      = "Expected " ++ show expected
                                       ++ " args; found values " ++ unwordsList found
showError (TypeMismatch expected found) = "Invalid type: expected " ++ expected
                                       ++ ", found " ++ show found
showError (Parser parseErr)             = "Parse error at " ++ show parseErr

trapError action = catchError action (return . show)

extractValue :: ThrowsError a -> a
extractValue (Right val) = val

type Env = IORef [(String, IORef LispVal)]
type IOThrowsError = ExceptT LispError IO

nullEnv :: IO Env
nullEnv = newIORef []

liftThrows :: ThrowsError a -> IOThrowsError a
liftThrows (Left err) = throwError err
liftThrows (Right val) = return val

runIOThrows :: IOThrowsError String -> IO String
runIOThrows action = runExceptT (trapError action) >>= return . extractValue


text:
{-# LANGUAGE DefaultSignatures #-}
{-# LANGUAGE ExistentialQuantification #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-|
Module      : Text.Loquacious
Description : Overridable Pretty Printed Values with Generic Annotations
Stability   : experimental

Loquacious enables pretty printing of values with the ability to override
individual typeclass instances in particular context, for example to provide
locali[sz]ation.
-}

module Text.Loquacious
  ( -- * Loquacious class and existential type.
    Loquacious(..)
  , Loq(..)
    -- * Handle Loqs, in particular to override them.
  , defLoq
  , LoqHandler(..)
  , loqHandle
  , loqHandles
  , overLoq
    -- * Helpers for polymorphic Loquacious instances
  , toLoqFunctor
  , fromLoqTraversable
  , toLoqBifunctor
  , fromLoqBitraversable
  , toLoqFirst
  , fromLoqFirst
  ) where

-- For core functionality.
import           Data.Maybe (fromMaybe)
import           Data.Text.Prettyprint.Doc (Doc, Pretty)
import           Data.Typeable (Typeable, cast)
import           Text.Loquacious.Annotations (Ann)

-- For helpers.
import           Data.Bifunctor (Bifunctor, bimap, first)
import           Data.Bitraversable (Bitraversable, bitraverse)
import           Data.Foldable (toList)
import           Data.Text.Prettyprint.Doc (align, colon, list, parens, pretty, tupled, viaShow,
                     vsep, (<+>))

-- For instances.
import qualified Control.Exception as E
import qualified Control.Exception.Base as EB
import           Data.Array.Unboxed (Array, IArray, Ix, UArray, elems)
import           Data.Complex (Complex((:+)))
import           Data.Functor.Const (Const, getConst)
import           Data.Functor.Identity (Identity, runIdentity)
import           Data.Int (Int16, Int32, Int64, Int8)
import           Data.List.NonEmpty (NonEmpty)
import           Data.Ratio (Ratio, denominator, numerator)
import           Data.String (fromString)
import qualified Data.Text as Text
import qualified Data.Text.Lazy as LazyText
import           Data.Typeable (Proxy(Proxy), TypeRep, typeOf, typeRep, (:~:), (:~~:))
import           Data.Version (Version, showVersion)
import           Data.Void (Void)
import           Data.Word (Word16, Word32, Word64, Word8)
import           GHC.Stack
import           Numeric.Natural (Natural)

-- | The root of the Loquacious Type Heirarchy, compare 'E.SomeException'.
data Loq = forall t. Loquacious t => Loq t

-- | Construct a 'Doc' with the default loq handlers (those specified in
-- the typeclass bodies).
defLoq :: Loquacious s => s -> Doc Ann
defLoq = loq defLoq

instance Show Loq where
  show = show . defLoq

instance Loquacious Loq where
  loq m (Loq x) = m x
  toLoq = id
  fromLoq = pure

-- | Types than can be turned into Loquaciouss.
class Typeable t => Loquacious t where
  -- | Produce a 'Doc'ument ready for pretty-printing that may or may not
  -- include annotations.
  --
  -- For base types this can simply be 'pretty', but higher kinded types should
  -- implement this manually.
  loq :: (forall s. Loquacious s => s -> Doc Ann) -> t -> Doc Ann
  default loq :: Pretty t => (forall s. Loquacious s => s -> Doc Ann) -> t -> Doc Ann
  loq _ = pretty

  -- | Embed a 'Loquacious' as a 'Loq'.
  --
  -- Higher Kinded Types should recursively encode their arguments! This allows
  -- polymorphic deconstruction of 'Loq'.
  --
  -- For 'Functor's and 'Bifunctor's use 'toLoqFunctor' and 'toLoqBifunctor',
  -- respectively.
  toLoq :: t -> Loq
  toLoq = Loq

  -- | Attempt to extract some 'Loquacious' from a 'Loq'.
  --
  -- For 'Data.Traversable.Traverable's and 'Data.Bitraversable.Bitraverable's
  -- use 'fromLoqTraversable' and 'fromLoqBitraversable', respectively.
  fromLoq :: Loq -> Maybe t
  fromLoq (Loq sm) = cast sm

-- | Wrap a function from some 'Loquacious' type.
data LoqHandler t = forall m. Loquacious m => LoqHandler (m -> t)

-- | Given some loq and a handler for some loq, apply the handler if
-- possible, otherwise return a default value.
loqHandle :: Loq -> LoqHandler t -> t -> t
loqHandle m (LoqHandler handle) def = fromMaybe def (handle <$> fromLoq m)

-- | Given some loq and a bunch of handlers for some loq, apply each
-- handler in turn if possible, otherwise return a default value.
loqHandles :: Foldable f => Loq -> t -> f (LoqHandler t) -> t
loqHandles m def hs = foldr (loqHandle m) def hs

-- | Override a Loquacious Handler!
overLoq :: (Foldable f, Loquacious s) => f (LoqHandler (Doc Ann)) -> s -> Doc Ann
overLoq hs x = loqHandles (Loq x) (loq (overLoq hs) x) hs

-- | Recursively encode a 'Functor' instance to 'Loq'.
--
-- Useful for types of kind '* -> *'.
toLoqFunctor :: (Loquacious (f Loq), Functor f, Loquacious a) => f a -> Loq
toLoqFunctor = Loq . fmap toLoq

-- | Decode a recursively encoded 'Loq' via a 'Traversable' instance.
--
-- Useful for types of kind '* -> *'.
fromLoqTraversable :: (Traversable t, Loquacious a, Typeable t) => Loq -> Maybe (t a)
fromLoqTraversable (Loq sm) = cast sm >>= traverse fromLoq

-- | Recursively encode a 'Bifunctor' instance to 'Loq'.
--
-- Useful for types of kind '* -> * -> *'.
toLoqBifunctor :: (Loquacious (f Loq Loq), Bifunctor f, Loquacious a, Loquacious b) => f a b -> Loq
toLoqBifunctor = Loq . bimap toLoq toLoq

-- | Decode a recursively encoded 'Loq' via a 'Traversable' instance.
--
-- Useful for types of kind '* -> *'.
fromLoqBitraversable :: (Bitraversable t, Loquacious a, Loquacious b, Typeable t) => Loq -> Maybe (t a b)
fromLoqBitraversable (Loq sm) = cast sm >>= bitraverse fromLoq fromLoq

-- | Recursively encode a 'Bifunctor' instance to 'Loq', but only on the first type argument.
--
-- Useful for types of kind '* -> k -> *', for example 'Const'.
toLoqFirst :: (Loquacious (f Loq b), Bifunctor f, Loquacious a) => f a b -> Loq
toLoqFirst = Loq . first toLoq

-- | Decode a recursively encoded 'Loq' via a 'Traversable' instance.
--
-- Useful for types of kind '* -> *'.
fromLoqFirst :: (Bitraversable t, Loquacious a, Typeable t, Typeable b) => Loq -> Maybe (t a b)
fromLoqFirst (Loq sm) = cast sm >>= bitraverse fromLoq pure

-- basic base types
instance Loquacious ()
instance Loquacious Bool
instance Loquacious Version where loq _ = pretty . showVersion
instance Loquacious Void

-- basic number types
instance Loquacious Integer
instance Loquacious Natural

instance Loquacious Int
instance Loquacious Int8
instance Loquacious Int16
instance Loquacious Int32
instance Loquacious Int64

instance Loquacious Word
instance Loquacious Word8
instance Loquacious Word16
instance Loquacious Word32
instance Loquacious Word64

instance Loquacious Double
instance Loquacious Float

-- polymorphic types
instance Loquacious a => Loquacious (Identity a) where
  loq m = m . runIdentity
  toLoq = toLoqFunctor
  fromLoq = fromLoqTraversable

instance (Loquacious a, Typeable b) => Loquacious (Const a b) where
  loq m = m . getConst
  toLoq = toLoqFirst
  fromLoq = fromLoqFirst

instance Loquacious a => Loquacious (Maybe a) where
  loq _ Nothing  = mempty
  loq m (Just x) = m x
  toLoq = toLoqFunctor
  fromLoq = fromLoqTraversable

instance Loquacious a => Loquacious (NonEmpty a) where
  loq m = align . list . fmap m . toList
  toLoq = toLoqFunctor
  fromLoq = fromLoqTraversable

instance Loquacious a => Loquacious [a] where
  loq m = align . list . fmap m
  toLoq = toLoqFunctor
  fromLoq = fromLoqTraversable

-- | The first paramater of 'Array' cannot be matched generically.
instance (Ix i, Typeable i, Loquacious a) => Loquacious (Array i a) where
  loq m = align . list . fmap m . toList
  toLoq = toLoqFunctor
  fromLoq = fromLoqTraversable

-- | 'UArray' cannot be matched generically due to boxing requirements.
instance (Ix i, Typeable i, IArray UArray a, Loquacious a) => Loquacious (UArray i a) where
  loq m = align . list . fmap m . elems

instance (Loquacious a, Loquacious b) => Loquacious (Either a b) where
  loq m (Left x)  = m x
  loq m (Right x) = m x
  toLoq = toLoqBifunctor
  fromLoq = fromLoqBitraversable

instance (Loquacious a, Loquacious b) => Loquacious (a, b) where
  loq m (a,b) = align $ tupled [m a, m b]
  toLoq = toLoqBifunctor
  fromLoq = fromLoqBitraversable

-- | FIXME The first paramater of cannot be matched generically, no @Tritraversable@.
instance (Loquacious a, Loquacious b, Loquacious c) => Loquacious (a, b, c) where
  loq m (a,b,c) = align $ tupled [m a, m b, m c]
  toLoq = toLoqBifunctor
  fromLoq = fromLoqBitraversable

instance Loquacious t => Loquacious (Complex t) where
  loq m (a :+ b) = m a <+> "+" <+> m b <> "i"
  toLoq = toLoqFunctor
  fromLoq = fromLoqTraversable

-- | FIXME Cannot match generic 'Ratio's due to lack of 'Functor' instance.
instance Loquacious t => Loquacious (Ratio t) where
  loq m x = m (numerator x) <> "/" <> m (denominator x)
  -- toLoq = toLoqFunctor
  -- fromLoq = fromLoqTraversable

-- basic string types
instance Loquacious Char
instance {-# OVERLAPS #-} Loquacious String
instance Loquacious Text.Text
instance Loquacious LazyText.Text
-- TODO instance Loquacious ByteString.ByteString
-- TODO instance Loquacious LazyByteString.ByteString

-- Call Stack Support
instance Loquacious CallStack where
  loq m = withFrozenCallStack ( align . vsep . fmap messageStack . getCallStack )
    where messageStack (func, srcLoc) = fromString func <+> "←" <+> m srcLoc

instance Loquacious SrcLoc where
  loq m SrcLoc{..} = msgFile <+> parens msgPackage
    where
      msgFile = fromString srcLocFile <> colon <> m srcLocStartLine <> colon <> m srcLocStartCol
      msgPackage = fromString srcLocPackage <> colon <> fromString srcLocModule

-- Proxy and other Type Representations
instance Loquacious TypeRep where loq _ = viaShow

instance (Typeable a) => Loquacious (Proxy a) where
  loq m _ = m $ typeRep (Proxy @a)

instance (Typeable a, Typeable b) => Loquacious (a -> b) where
  loq m = m . typeOf

instance (Typeable a, Typeable b) => Loquacious (a :~: b) where
  loq m _ = m (typeRep (Proxy @a)) <+> ":~:" <+> m (typeRep (Proxy @b))

instance (Typeable a, Typeable b) => Loquacious (a :~~: b) where
  loq m _ = m (typeRep (Proxy @a)) <+> ":~~:" <+> m (typeRep (Proxy @b))

-- Exception Types
instance Loquacious E.AllocationLimitExceeded where loq _ = viaShow
instance Loquacious E.ArithException where loq _ = viaShow
instance Loquacious E.ArrayException where loq _ = viaShow
instance Loquacious E.AssertionFailed where loq _ = viaShow
instance Loquacious E.AsyncException where loq _ = viaShow
instance Loquacious E.BlockedIndefinitelyOnMVar where loq _ = viaShow
instance Loquacious E.BlockedIndefinitelyOnSTM where loq _ = viaShow
instance Loquacious E.CompactionFailed where loq _ = viaShow
instance Loquacious E.Deadlock where loq _ = viaShow
instance Loquacious E.IOException where loq _ = viaShow
instance Loquacious E.NestedAtomically where loq _ = viaShow
instance Loquacious E.NoMethodError where loq _ = viaShow
instance Loquacious E.NonTermination where loq _ = viaShow
instance Loquacious E.PatternMatchFail where loq _ = viaShow
instance Loquacious E.RecConError where loq _ = viaShow
instance Loquacious E.RecSelError where loq _ = viaShow
instance Loquacious E.RecUpdError where loq _ = viaShow
instance Loquacious E.SomeAsyncException where loq _ = viaShow
instance Loquacious E.SomeException where loq _ = viaShow
instance Loquacious E.TypeError where loq _ = viaShow
instance Loquacious EB.FixIOException where loq _ = viaShow

max_stars_repo_name: imuli/loquacious

max_issues_repo_name: imuli/loquacious

max_forks_repo_name: imuli/loquacious


text:
{-# LANGUAGE FlexibleContexts #-}

module Stastics where

import           Data.Heap
import           Data.Maybe                  (fromJust)

import qualified Data.Vector.Unboxed         as U
import           Graphics.EasyPlot
import           Statistics.LinearRegression

-- Moving/Running Mean

-- convert raw data to Double
clean raw = map (\s -> read s :: Double) (lines raw)

-- Calculate moving average -- hardcoded smoothening factor = 0.95
avg :: [Double] -> Double
avg xs = case xs of
    (x:xs) -> a * x + (1 - a) * (avg xs)
    []     -> 0
    where a = 0.95

movingMean = do
  rawInput <- readFile "numbers.txt"
  let input = avg . clean $ rawInput
  print input

-- Calculating a Moving Median
strToInt raw = map (\s -> read s :: Int) (lines raw)
-- put smaller nums in maxheap and big nums in minheap
median (x:xs) maxheap minheap = case viewHead maxheap of
    Just theMax -> if x < theMax
            then median xs (insert x maxheap) minheap
            else median xs maxheap (insert x minheap)
    Nothing -> median xs (insert x maxheap) minheap

median [] maxheap minheap
    | size maxheap + 1 < size minheap =
        median [] (insert minelem maxheap) $ (snd.fromJust.view) minheap
    | size minheap + 1 < size maxheap =
        median [] ((snd.fromJust.view) maxheap) $ insert maxelem minheap
    |  size maxheap == size minheap =
        (fromIntegral maxelem + fromIntegral minelem) / 2.0
    | size maxheap > size minheap = fromIntegral maxelem
    | otherwise = fromIntegral minelem
    where
        maxelem = fromJust (viewHead maxheap)
        minelem = fromJust (viewHead minheap)

movingMedian = do
    raw <- readFile "numbers.txt"
    let input = strToInt raw
    print $ median input (empty :: MaxHeap Int) (empty :: MinHeap Int)

-- APPROXIMATING A LINEAR REGRESSION
linearR =do
    let xs = U.fromList [1.0, 2.0, 3.0, 4.0, 5.0] :: U.Vector Double
    let ys = U.fromList [1.0, 2.0, 1.3, 3.75, 2.25] :: U.Vector Double
    let (b,m) = linearRegression xs ys
    print $ concat ["y = ", show m, " x + ", show b]

--
main = do
        let values = [4,5,16,15,14,13,13,17]
        let values2 = [2,5,7,10,21, 18, 15]
        -- plot X11 $
        --    Data2D [ Title "Line Graph"
        --           , Style Linespoints
        --           , Color Blue]
        --    [] (zip [1..] values)

        -- plot X11 $ [ Data2D [Color Red] [] (zip [1..] values)
        --             , Data2D [Color Blue] [] (zip [1..] values2) ]

        -- 3D
        -- plot' [Interactive] X11
        --                [ Data3D [Color Red] [] [(1,1,1), (1,2,1), (0,1,1), (1,1,0)]
        --                , Data3D [Color Blue] [] [(4,3,2), (3,3,2), (3,2,3), (4,4,3)] ]

        -- regression

        let winEstimate = map (\x -> x * 0.425 + 0.7850000000000001) [1.3, 2.4 .. 4.7]
        let regressionLine = zip [1.3, 2.4 .. 4.7] winEstimate
        plot X11 [Data2D [Title "Runs Per Game VS Win % in 2014" ] [] (zip [1.0, 2.0, 3.0, 4.0, 5.0] [1.0, 2.0, 1.3, 3.75, 2.25]), Data2D [Title "Regression Line", Style Lines, Color Blue] [] regressionLine]

max_stars_repo_name: alokpndy/haskell-learn

max_issues_repo_name: alokpndy/haskell-learn

max_forks_repo_name: alokpndy/haskell-learn


text:
module EndPointMarkovChain where

import           Control.Monad           as M
import           Data.Array.Repa         as R
import           Data.List               as L
import           Graphics.Gnuplot.Simple
import           Numeric.LinearAlgebra   as NL
import           System.Directory
import           System.Environment
import           System.FilePath         ((</>))
import           Text.Printf

{-# INLINE splitList #-}
splitList :: Int ->  [a] -> [[a]]
splitList _ [] = []
splitList n xs = (L.take n xs) : splitList n (L.drop n xs)

{-# INLINE makeMatrix #-}
makeMatrix :: Int -> Double -> R.Array D DIM2 Double
makeMatrix n p =
  fromFunction (Z :. (2 * n) :. (2 * n)) $ \(Z :. i :. j) ->
    let a = div i 2
        b = div j 2
        c = mod i 2
        d = mod j 2
     in if abs (a - b) == 1
          then if a == 0
                 then if d == 1
                        then if c == 0
                               then p
                               else 1
                        else 0
                 else if a == n - 1
                        then if d == 0
                               then if c == 0
                                      then 1
                                      else p
                               else 0
                        else if b < a
                               then if d == 0
                                      then if c == 0
                                             then 1
                                             else p
                                      else if c == 0 -- 0
                                              then p
                                              else 0
                               else if d == 0
                                      then  if c == 0 -- 0
                                               then 0
                                               else p
                                      else if c == 0
                                             then p
                                             else 1
          else 0

{-# INLINE powerMethodIO #-}
powerMethodIO :: Matrix Double -> Vector Double -> Int -> IO (Vector Double)
powerMethodIO mat vec n = do
  let (m,_) = NL.size mat
      x =  ((inv (mat - (0.1) * ident m)) #> vec)
      maxV = sumElements x -- maxElement x
      eigenVec = x / scalar maxV
  printf "%d: %f\n" n (rayleighQuotient mat eigenVec)
  printf "%03d: %s" n (printList . sumVec  . NL.toList $ eigenVec)
  return eigenVec

{-# INLINE powerMethod #-}
powerMethod :: Matrix Double -> Vector Double -> (Double, Vector Double)
powerMethod mat vec =
  let x = mat #> vec
      maxV = maxElement x
      eigenVec = x / scalar maxV
      rq = rayleighQuotient mat eigenVec
   in (rq, eigenVec)

{-# INLINE sumVec #-}
sumVec :: [Double] -> [Double]
sumVec []       = []
sumVec (x:y:xs) = (x + y) : sumVec xs

{-# INLINE printList #-}
printList :: [Double] -> String
printList []     = "\n"
printList (x:xs) = printf "%.3f " x L.++ printList xs

{-# INLINE rayleighQuotient #-}
rayleighQuotient :: Matrix Double -> Vector Double -> Double
rayleighQuotient mat vec = vec <.> (mat #> vec) / vec <.> vec

{-# INLINE inverse #-}
inverse :: Vector Double -> Vector Double
inverse vec = cmap (\x -> sumElements vec - x ) vec

{-# INLINE timeReverse #-}
timeReverse :: [a] -> [a]
timeReverse [] = []
timeReverse (x:y:xs) = y : (x : timeReverse xs)

main = do
  args <- getArgs
  let (nStr:aStr:mStr:_) = args
      n = read nStr :: Int
      a = read aStr :: Double
      m = read mStr :: Int
      mat = (tr . ((2 * n) >< (2 * n)) . R.toList . makeMatrix n $ a) + (0.0) * ident (2 * n)
      -- mat = NL.fromLists . timeReverse . L.map timeReverse . NL.toLists $ mat'
      vec = vector . L.take (2 * n) $ [1,1 ..]
      (eigVal', eigVec') = eig mat
      (maxEigVal, maxEigVec) =
        (\(x, y) -> (realPart x, sumVec $ L.map realPart y)) .
        L.maximumBy (\x y -> compare (realPart $ fst x) (realPart $ fst y)) .
        L.zip (NL.toList eigVal') . L.map NL.toList . NL.toColumns $
        eigVec'
      (maxEigVal', maxEigVec') =
        (\(x, y) -> (realPart x, L.map realPart y)) .
        L.maximumBy (\x y -> compare (realPart $ fst x) (realPart $ fst y)) .
        L.zip (NL.toList eigVal') . L.map NL.toList . NL.toColumns $
        eigVec'
      maxV = L.maximumBy (\x y -> compare (abs x) (abs y)) maxEigVec
      (rqs, eigVecs) =
        L.unzip $ L.scanl' (\(_, v) _ -> powerMethod mat v) (0, vec) [1 .. m]
      -- eigVec = L.last eigVecs
      folderPath = "output/test/EndPointMarkovChain"
  createDirectoryIfMissing True folderPath
  -- eigVec <- M.foldM (powerMethodIO mat) vec [1 .. m]
  -- -- print . rayleighQuotient mat $ eigVec
  -- -- putStr . printList . sumVec . NL.toList $ eigVec
  -- -- print maxEigVal
  -- -- putStr . printList . L.map (/ maxV) $ maxEigVec
  -- printf
  --   "%.3f: %s"
  --   (rayleighQuotient mat eigVec)
  --   (printList . sumVec . NL.toList $ eigVec)
  -- printf
  --   "%.3f: %s"
  --   maxEigVal
  --   -- (L.sum . L.map (/ maxV) $ maxEigVec)
  --   (printList . L.map (/ maxV) $ maxEigVec)
  -- plotPathStyle
  --   [ PNG (folderPath </> "RayleighQuotient.png")
  --   , Title ("Rayleigh Quotient")
  --   , XLabel "Iteration"
  --   ]
  --   (defaultStyle {plotType = Lines, lineSpec = CustomStyle [LineTitle ""]}) .
  --   L.zip [1 ..] . L.tail $
  --   rqs
  -- let m = 10
  --     croppedVec =
  --       L.take (2 * m) . L.drop (div (2 * n - 2 * m) 2) $
  --       L.map (/ maxV) $ maxEigVec'
  -- print . L.sum $ croppedVec
  M.mapM_ print .
    L.map
      (\(x, ys) ->
         let v = L.maximumBy (\a b -> compare (abs a) (abs b)) ys
          in (x, L.map (/ v) ys)) .
    L.map (\(x, y) -> (realPart x, sumVec $ L.map realPart y)) .
    L.sortBy (\x y -> compare (realPart $ fst x) (realPart $ fst y)) .
    L.zip (NL.toList eigVal') . L.map NL.toList . NL.toColumns $
    eigVec'
  -- M.mapM_  print . NL.toLists $ (mat + 1 * ident 20)
  -- M.mapM_  print . NL.toLists . inv $ (mat + 1 * ident 20)
  -- M.mapM_ print . NL.toLists $ mat - (0.017) * ident (2 * n)
  -- putStrLn ""
  M.mapM_ (putStr . printList) . NL.toLists $ (inv $ (mat - (0.6) * ident (2*n))) / 59.592--2261.598 --494.987 -- 4962.166
  --   (mat + 0.001 * ident 20) NL.<> (tr $ (mat + 0.001 * ident 20))

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
module Ch02.Perceptron where

import Prelude hiding (and, or)
import Numeric.LinearAlgebra

and_ :: (Ord a, Fractional a, Num b) => a -> a -> b
and_ x1 x2 = if (x1*w1 + x2*w2) <= theta then 0 else 1
  where
    w1    = 0.5
    w2    = 0.5
    theta = 0.7
-- λ> and_ 0 0
-- 0
-- λ> and_ 1 0
-- 0
-- λ> and_ 0 1
-- 0
-- λ> and_ 1 1
-- 1

perceptron :: R -> R -> R -> R -> R -> R
perceptron w1 w2 b = \x1 x2 ->
  if (sum $ toList $ (vector [x1,x2]) * vector [w1,w2]) + b <= 0 then 0 else 1

and       = perceptron   0.5    0.5  (-0.7)
nand      = perceptron (-0.5) (-0.5)   0.7
or        = perceptron   0.5    0.5  (-0.2)
xor x1 x2 = and (nand x1 x2) (or x1 x2)


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

import Data.Complex

import Debug.Trace

data Nav = N Float
         | S Float
         | E Float
         | W Float
         | L Int
         | R Int
         | F Float
 deriving (Show)

data Change = Dir (Complex Float)
            | Rot (Complex Float)
            | Move (Complex Float)

 deriving (Show)

degToRot :: Int -> Float
degToRot d = fromIntegral (d `div` 90)*pi/2

navToChange :: Nav -> Change
navToChange (N m) = Dir  $ 0 :+ m
navToChange (S m) = Dir  $ 0 :+ (-m)
navToChange (E m) = Dir  $ m :+ 0
navToChange (W m) = Dir  $ (-m) :+ 0
navToChange (L d) = Rot  $ mkPolar 1.0 (degToRot d)
navToChange (R d) = Rot  $ mkPolar 1.0 (-(degToRot d))
navToChange (F m) = Move $ mkPolar m 0

instance Read Nav where
  readsPrec _ ('N':r) = [( N $ read r, "")]
  readsPrec _ ('S':r) = [( S $ read r, "")]
  readsPrec _ ('E':r) = [( E $ read r, "")]
  readsPrec _ ('W':r) = [( W $ read r, "")]
  readsPrec _ ('L':r) = [( L $ read r, "")]
  readsPrec _ ('R':r) = [( R $ read r, "")]
  readsPrec _ ('F':r) = [( F $ read r, "")]

getInput :: FilePath -> IO [Nav]
getInput = fmap (map read . lines) . readFile

dirToComplex :: Nav -> Complex Float
dirToComplex (N n) = 0 :+ n
dirToComplex (S s) = 0 :+ (-s)
dirToComplex (E e) = e :+ 0
dirToComplex (W w) = (-w) :+ 0
dirToComplex _ = undefined

solution :: [Nav] -> Int
solution navs = round $ abs (dx + ix) + abs (dy + iy)
  where (dirs, instr) = (filter isDir navs, filter (not . isDir) navs)
        (dx :+ dy) = sum $ map dirToComplex dirs
        (ix :+ iy) = instrsToComplex 0 (0 :+ 0) instr

        isDir :: Nav -> Bool
        isDir (N _) = True
        isDir (S _) = True
        isDir (E _) = True
        isDir (W _) = True
        isDir _ = False

        instrsToComplex :: Float -> Complex Float -> [Nav] -> Complex Float
        instrsToComplex _ curLoc [] = curLoc
        instrsToComplex curPhase curLoc ((L d):rest)
            = instrsToComplex (curPhase + degToRot d) curLoc rest
        instrsToComplex curPhase curLoc ((R d):rest)
            = instrsToComplex (curPhase - degToRot d) curLoc rest
        instrsToComplex curPhase ship ((F m):rest)
            = instrsToComplex curPhase (ship+rot) rest
            where rot = mkPolar m curPhase

solution2 :: [Nav] -> Int
solution2 navs = round $ abs ix + abs iy
  where (ix :+ iy) = run 0 (10 :+ 1) $ map navToChange navs
        run :: Complex Float -> Complex Float -> [Change] -> Complex Float
        run ship _ [] = ship
        run ship wp (Dir d:rest) = run ship (wp + d) rest
        run ship wp (Rot r:rest) = run ship (wp*r) rest
        run ship wp (Move m:rest)= run (ship + m*wp) wp rest

main :: IO ()
main = do getInput "test-input" >>= print . solution
          getInput "test-input" >>= print . solution2
          getInput "input" >>= print . solution
          getInput "input" >>= print . solution2


text:
{-# LANGUAGE RankNTypes #-}

module Sim
  ( State (..),
    Plant (..),
    Measurement (..),
    simulate,
  )
where

import Numeric.LinearAlgebra

type Measurement = Vector Double

type State = Vector Double

data Plant = Plant {fA, fC, fB, fD :: forall a. Floating a => ([a] -> [a])}

simulate :: Plant -> State -> Int -> [(State, Measurement)]
simulate (Plant fA fB fC fD) x n = (x', y') : simulate (Plant fA fB fC fD) xh (n + 1)
  where
    x' = vector $ fA (toList x)
    y = vector $ fC (toList x')
    xh = x' + scale 0.0001 (randomVector (n + 1000) Gaussian (size x))
    y' = y + scale 0.5 (randomVector n Gaussian (size y))


text:
module Signal where

import Control.Monad.Fix
import Data.Complex

type Time = Double
type Metric = Complex Double
newtype Signal a = Signal {at :: Metric -> a}

instance Functor Signal where
    fmap f (Signal a) = Signal $ \time -> f (a time)

instance Applicative Signal where
    pure x = Signal $ const x
    sf <*> (Signal a) = Signal $ \time -> (sf `at` time) (a time)

instance Monad Signal where
    return = pure
    sa >>= f = Signal $ \time -> f (sa `at` time) `at` time

instance MonadFix Signal where
    mfix f = Signal $ \time -> fix ((`at` time) . f)

instance Eq a => Eq (Signal a) where
    (Signal x) == (Signal y) = x 0 == y 0

instance Ord a => Ord (Signal a) where
    (Signal x) <= (Signal y) = x 0 <= y 0

instance Num a => Num (Signal a) where
    sx + sy = Signal $ \time -> (sx `at` time) + (sy `at` time)
    sx * sy = Signal $ \time -> (sx `at` time) * (sy `at` time)
    abs sx = Signal $ \time -> abs $ sx `at` time
    fromInteger i = Signal $ const $ fromInteger i
    signum sx = Signal $ \time -> signum $ sx `at` time
    negate sx = Signal $ \time -> negate $ sx `at` time

instance Fractional a => Fractional (Signal a) where
    sx / sy = Signal $ \time -> (sx `at` time) / (sy `at` time)
    fromRational i = Signal $ const $ fromRational i

instance Show a => Show (Signal a) where
    show (Signal a) = show $ a 0

max_stars_repo_name: FP-Modeling/fixingAnalog

max_issues_repo_name: FP-Modeling/fixingAnalog

max_forks_repo_name: FP-Modeling/fixingAnalog


text:
-- example 3.11

import Statistics.Distribution
import Statistics.Distribution.Geometric

g :: GeometricDistribution
g = geometric 0.02

main = do
  putStrLn $ "the rate of engines which survive two hours: "
          ++ show (1 - cumulative g 2.0)
  putStrLn $ "the rate of engines which survive two hours: "
          ++ show (complCumulative g 2.0)
  putStrLn $ "the rate of engines which survive two hours: "
          ++ show (1 - sum (map (probability g) [0 .. 2]))

max_stars_repo_name: raydsameshima/Stat

max_issues_repo_name: raydsameshima/Stat

max_forks_repo_name: raydsameshima/Stat


text:
{-# LANGUAGE TemplateHaskell, RecordWildCards, ScopedTypeVariables #-}
module TestEvolution where

import qualified Data.ByteString.Lazy          as B
import           Data.MonoTraversable           ( MonoFoldable(olength) )
import           Data.Ord                       ( Down(Down) )
import qualified Data.Set                      as S
import           Data.Typeable                  ( Typeable )
import           Data.Word                      ( Word8 )
import           Instances                      ( )
import           Internal.Evolution
import           Statistics.ConfidenceInt       ( binomialCI )
import           Statistics.Types               ( ConfInt(ConfInt)
                                                , Estimate(Estimate)
                                                , cl95
                                                )
import           Test.QuickCheck                ( Arbitrary(arbitrary)
                                                , CoArbitrary
                                                , Gen
                                                , NonNegative(getNonNegative)
                                                , NonPositive(getNonPositive)
                                                , Positive(getPositive)
                                                , Property
                                                , choose
                                                , cover
                                                , elements
                                                , listOf1
                                                , quickCheckAll
                                                , suchThat
                                                , vectorOf
                                                )
import           Util

newtype ReproduceArgs a = ReproduceArgs (EvolutionParams, [a] -> [Float], [Agent a]) deriving (Show)

instance (Show a, Typeable a, Eq a, Arbitrary a, CoArbitrary a) => Arbitrary (ReproduceArgs a) where
    arbitrary = do
        params@EvolutionParams {..} <- arbitrary
        agents                      <- genPop evolvePopSize =<< arbitrary
        -- generate a fitness fn that strictly increases fitness
        -- fitness starts at 0 so this also guarantees fitness will be > 0
        -- necessary because fitnesses are passed as weights to `frequency`
        f                           <- (arbitrary :: Gen (a -> Positive Float))
        return $ ReproduceArgs (params, map (getPositive . f), agents)

newtype MutateArgs = MutateArgs (Double, [Word8], B.ByteString) deriving (Show, Eq)

instance Arbitrary MutateArgs where
    arbitrary = do
        p     <- choose (0, 1)
        genes <- listOf1 arbitrary
            `suchThat` combineWith (&&) [unique, (> 1) . length]
        genome <- B.pack <$> vectorOf 1000 (elements genes)
        return $ MutateArgs (p, genes, genome)

matchup :: [Agent a] -> [[Agent a]]
--matchup = matchWithBest 2 10 (Down . agentFitness) . S.fromDistinctAscList
matchup = matchups 2 . S.fromDistinctAscList

genPop :: (Arbitrary a) => Int -> Agent a -> Gen [Agent a]
genPop n Agent {..} =
    genPopulation n agentGenes agentEncoder agentDecoder arbitrary

prop_mergeAgentsUnique :: [Agent a] -> Bool
prop_mergeAgentsUnique = unique . mergeAgents

prop_mergeAgentsKeepAll :: [Agent a] -> Bool
prop_mergeAgentsKeepAll xs = all (`elem` mergeAgents xs) xs

prop_mergeAgentsNoNew :: [Agent a] -> Bool
prop_mergeAgentsNoNew xs = all (`elem` xs) (mergeAgents xs)

prop_populationLength :: Arbitrary a => NonNegative Int -> Agent a -> Gen Bool
prop_populationLength n agent = do
    pop <- genPop (getNonNegative n) agent
    return $ length pop == getNonNegative n

prop_populationIds :: Arbitrary a => NonNegative Int -> Agent a -> Gen Bool
prop_populationIds n agent = do
    pop <- genPop (getNonNegative n) agent
    return . and $ zipWith (==) (map agentId pop) [0 ..]

prop_populationUniform
    :: (Eq a, Arbitrary a) => NonNegative Int -> Agent a -> a -> Gen Bool
prop_populationUniform n agent c = do
    pop <- genPop (getNonNegative n) agent
    return $ uniform pop True
  where
    uniform []       acc = acc
    uniform (x : xs) acc = uniform xs (acc && all (agentCommon x) xs)
    agentCommon x y = fEq enc && fEq dec
      where
        fEq f = f x == f y
        enc z = agentEncoder z c
        dec z = agentDecoder z $ enc z

prop_mutate :: MutateArgs -> Gen Property
prop_mutate (MutateArgs (p, genes, genome)) = do
    mutated <- mutate p genes genome
    let Estimate p' (ConfInt lower upper _) = binomialCI
            cl95
            (olength genome)
            (olength $ filter id (B.zipWith (/=) genome mutated))
    -- proportion of mutations == p with 95% confidence
    return $ cover 95 ((p > p' - lower) && (p < p' + upper)) "near p" True

prop_reproduceLength :: ReproduceArgs a -> Gen Bool
prop_reproduceLength (ReproduceArgs (params, f, pop)) = do
    pop' <- reproduce params f matchup pop
    return $ length pop == length pop'

prop_reproduceIds :: ReproduceArgs a -> Gen Bool
prop_reproduceIds (ReproduceArgs (params, f, pop)) =
    fIncreasing True agentId <$> reproduce params f matchup pop
  where
    fIncreasing acc _ []           = acc
    fIncreasing acc _ [_         ] = acc
    fIncreasing acc g (x : y : xs) = fIncreasing (acc && g x < g y) g (y : xs)

prop_evolveId :: ReproduceArgs a -> NonPositive Int -> Gen Bool
prop_evolveId (ReproduceArgs (params, f, pop)) n =
    (pop ==)
        <$> evolve params { evolveGenerations = getNonPositive n } f matchup pop

prop_evolvePreserve :: ReproduceArgs a -> Gen Bool
prop_evolvePreserve (ReproduceArgs (params, f, pop)) = do
    pop' <- evolve params f matchup pop
    let args' = ReproduceArgs (params, f, pop')
    len <- prop_reproduceLength args'
    ids <- prop_reproduceIds args'
    return $ len && ids

return []
runTests :: IO Bool
runTests = $quickCheckAll

max_stars_repo_name: desophos/gamer-evolution

max_issues_repo_name: desophos/gamer-evolution

max_forks_repo_name: desophos/gamer-evolution


text:
module Main where

import Test.Hspec
-- import Test.Hspec.QuickCheck

-- import qualified Statistics.Matrix as M
-- import qualified Statistics.Matrix.Types as T
-- import qualified Statistics.Matrix.Algorithms as Alg

import qualified Fixtures as F
import AlgorithmsSpec
-- import Utils

spec :: Spec
spec = describe "Q-R Decomposition" $ do
        it "Q x R returns the original matrix" $ do
          qrDecompositionInvariant F.matB `shouldBe` True
        it "Matrix Q is orthogonal for the QR factorization of an invertible matrix" $ do
          qrFirstOrthoInvariant F.matD`shouldBe` True
        it "Matrix R is triangular" $ do
          qrSecondTriInvariant F.matB `shouldBe` True

main :: IO ()
main = do
  hspec spec

max_stars_repo_name: Magalame/dh-core

max_issues_repo_name: Magalame/dh-core

max_forks_repo_name: Magalame/dh-core


text:
module Sondages where

import qualified Data.Text as T
import Data.Time (UTCTime)
import qualified Data.Map as M
import Data.Time.Segment
import Statistics.LinearRegression (linearRegression)
import qualified Data.Vector.Unboxed as U

data Candidat = Macron | LePen | Clinton | Trump
    deriving (Show, Read, Eq, Ord, Enum, Bounded)

data VoteIntention = VoteIntention { unCandidat :: Candidat
                                   , unIntention :: Double
             } deriving (Show, Read, Eq)

data Sondage = Sondage { unSondeur   :: T.Text
                       , unDateBegin :: UTCTime
                       , unDateEnd   :: Maybe UTCTime
                       , echantillon :: Maybe Int
                       , margeErreur :: Maybe Double
                       , abstention  :: Maybe Double
                       , indecis     :: Maybe Double
                       , intentions  :: [(Candidat, Double)]
                       } deriving (Show, Read, Eq)

sondages :: [Sondage]
sondages = [ Sondage "Harris"     (jour 2017 04 23) Nothing                  (Just 2684)   Nothing   Nothing     Nothing            [(Macron, 0.64)]
           , Sondage "Ipsos"      (jour 2017 04 23) Nothing                  (Just 1379)   Nothing   Nothing     Nothing            [(Macron, 0.62)]
           , Sondage "OpinionWay" (jour 2017 04 23) (Just $ jour 2017 04 24) (Just 1461)   Nothing   Nothing     Nothing     [(Macron, 0.61)]
           , Sondage "Ifop"       (jour 2017 04 23) (Just $ jour 2017 04 24) (Just 846)    Nothing  (Just 0.26) (Just 0.12)  [(Macron, 0.60)]
           , Sondage "Harris"     (jour 2017 04 25) (Just $ jour 2017 04 27) (Just 1016)   Nothing   Nothing     Nothing     [(Macron, 0.61)]
           , Sondage "OpinionWay" (jour 2017 04 25) (Just $ jour 2017 04 27) (Just 1500)   Nothing  (Just 0.25)  Nothing     [(Macron, 0.60)]
           , Sondage "Odoxa"      (jour 2017 04 26) (Just $ jour 2017 04 27) (Just 1003)   Nothing  (Just 0.20)   Nothing     [(Macron, 0.59)]
           , Sondage "BVA"        (jour 2017 04 26) (Just $ jour 2017 04 28) (Just 1506)   Nothing  (Just 0.245) (Just 0.15) [(Macron, 0.59)]
           , Sondage "Ipsos"      (jour 2017 04 28) (Just $ jour 2017 04 29) (Just 1504)   Nothing  (Just 0.25)   Nothing     [(Macron, 0.60)]
           , Sondage "Kantar"     (jour 2017 04 28) (Just $ jour 2017 04 30) (Just 1539)   Nothing   Nothing      Nothing     [(Macron, 0.59)]
           , Sondage "Elabe"      (jour 2017 04 28) (Just $ jour 2017 05 02) (Just 3956)   Nothing   Nothing     (Just 0.11)  [(Macron, 0.59)]
           , Sondage "Ifop"       (jour 2017 04 28) (Just $ jour 2017 05 02) (Just 1388)   Nothing  (Just 0.27)   Nothing     [(Macron, 0.595)]
           , Sondage "Cevifop"    (jour 2017 04 30) (Just $ jour 2017 05 01) (Just 13742)  Nothing  (Just 0.24)  (Just 0.15)  [(Macron, 0.59)]
           , Sondage "BVA"        (jour 2017 05 01) (Just $ jour 2017 05 02) (Just 1435)   Nothing  (Just 0.22)   Nothing     [(Macron, 0.60)]
           , Sondage "OpinionWay" (jour 2017 05 02) (Just $ jour 2017 05 04) (Just 1500)   Nothing  (Just 0.25)  (Just 0.15)  [(Macron, 0.62)]
           , Sondage "Elabe"      (jour 2017 05 04)  Nothing                 (Just 1009)   Nothing   Nothing      Nothing     [(Macron, 0.62)]
           , Sondage "Ipsos"      (jour 2017 05 04)  Nothing                 (Just 2632)   Nothing  (Just 0.24)   Nothing     [(Macron, 0.615)]
           , Sondage "Odoxa"      (jour 2017 05 04)  Nothing                 (Just 998)    Nothing  (Just 0.25)  (Just 0.21)  [(Macron, 0.62)]
           , Sondage "Harris"     (jour 2017 05 02) (Just $ jour 2017 05 05) (Just 4991)   Nothing  (Just 0.25)   Nothing     [(Macron, 0.62)]
           , Sondage "Ifop"       (jour 2017 05 02) (Just $ jour 2017 05 05) (Just 1861)   Nothing  (Just 0.245)  Nothing     [(Macron, 0.63)]
           , Sondage "Ipsos"      (jour 2017 05 05)  Nothing                 (Just 5331)   Nothing  (Just 0.24)  Nothing      [(Macron, 0.63)]
           ]

polls :: [Sondage]
polls = [ Sondage "YouGov" (jour 2016 10 22) (Just $ jour 2016 10 26) (Just 1376) (Just 0.031) Nothing Nothing [(Clinton, 0.49), (Trump, 0.46)]
        , Sondage "ABC"    (jour 2016 10 23) (Just $ jour 2016 10 26) (Just 1150) (Just 0.03)  Nothing Nothing [(Clinton, 0.50), (Trump, 0.45)]
        , Sondage "USC"    (jour 2016 10 21) (Just $ jour 2016 10 27) (Just 3248) (Just 0.045) Nothing Nothing [(Clinton, 0.44), (Trump, 0.46)]
        , Sondage "Ipsos"  (jour 2016 10 21) (Just $ jour 2016 10 27) (Just 1627) (Just 0.03)  Nothing Nothing [(Clinton, 0.42), (Trump, 0.36)]
        , Sondage "IBD"    (jour 2016 10 22) (Just $ jour 2016 10 27) (Just 973)  (Just 0.033) Nothing Nothing [(Clinton, 0.45), (Trump, 0.42)]
        , Sondage "ABC"    (jour 2016 10 24) (Just $ jour 2016 10 27) (Just 1148) (Just 0.03)  Nothing Nothing [(Clinton, 0.49), (Trump, 0.46)]
        , Sondage "IBD"    (jour 2016 10 23) (Just $ jour 2016 10 28) (Just 1013) (Just 0.033) Nothing Nothing [(Clinton, 0.46), (Trump, 0.41)]
        , Sondage "ABC"    (jour 2016 10 25) (Just $ jour 2016 10 28) (Just 1160) (Just 0.03)  Nothing Nothing [(Clinton, 0.46), (Trump, 0.45)]
        ]

pollsResult :: [Sondage]
pollsResult = [Sondage "REAL" (jour 2016 11 08) Nothing Nothing Nothing Nothing Nothing [(Clinton, 0.481), (Trump, 0.46)]]

normalizePolls :: [Sondage] -> [Sondage]
normalizePolls xs = Prelude.map normalizePolls' xs
    where
        normalizePolls' (Sondage a b c d e f g is) = Sondage a b c d e f g is'
            where
                is' = Prelude.map (\(m, n) -> (m, n/total)) is
                total = sum $ Prelude.map snd is

sondage2date :: [Sondage] -> [(UTCTime, [(Candidat, Double)])]
sondage2date xs = concat ( map (dureeSondage) xs )
    where
        dureeSondage (Sondage _ begin end _ _ _ _ is) = case end of
                    Nothing   -> zip (timesBetween begin begin D) (repeat is)
                    Just end' -> zip (timesBetween begin end'  D) (repeat is)

toMap :: forall k a. Ord k => [(k, [a])] -> M.Map k [a]
toMap = Prelude.foldr (\s -> M.insertWith (++) (fst s) (snd s)) M.empty

toMap' :: forall k a. Ord k => [(k, a)] -> M.Map k [a]
toMap' = Prelude.foldr (\s -> M.insertWith (++) (fst s) ([snd s])) M.empty

mapSondages :: [Sondage] -> M.Map UTCTime (M.Map Candidat Double)
mapSondages sondage = M.map (M.map mean) $ M.map toMap' $ toMap (sondage2date sondage)
    where
        mean :: [Double] -> Double
        mean xs = (sum xs) / (fromIntegral $ length xs)

predict :: (Double -> Double) -> Int -> [(UTCTime, Double)]
predict f days = predict' sondages f days

predict' :: [Sondage] -> (Double -> Double) -> Int -> [(UTCTime, Double)]
predict' sondes f days = zip dates $ Prelude.map (\x -> f $ a' + b' * x) [1.. (fromIntegral n)]
    where
        sondes' = M.toAscList (mapSondages sondes)
        firstDate = fst $ head sondes'
        n = length sondes' + days
        dates = timesAfter n D firstDate
        (a', b') = linearRegression (U.fromList a) (U.fromList b)
        (a, b) = unzip $ zip ([1..] :: [Double]) $  Prelude.map snd $ concat $ Prelude.map M.toList $ Prelude.map snd sondes'

max_stars_repo_name: delanoe/vote-abstention

max_issues_repo_name: delanoe/vote-abstention

max_forks_repo_name: delanoe/vote-abstention


text:
{-# LANGUAGE FlexibleContexts #-}

module WeightDecay
  ( input
  , testdata
  , traindata
  , testY
  , testX
  , trainLinRegUnkownTarget
  , trainAndTestWithData
  , trainAndTestWithRegularization
  , makeTransformedMatrix
  , listToTuple'
  ) where

import           Data.Maybe
import           LinearRegression
import           Network.HTTP
import           NonlinearTransform
import           Numeric.LinearAlgebra.Data
import           Numeric.LinearAlgebra.HMatrix as HMatrix

--------------------------------------------------------------
-- Solutions to homework 6 of "Learning from data" question 5-6
--------------------------------------------------------------
input http = do
  response <- simpleHTTP $ getRequest http
  let body = fmap rspBody response
  let rows = fmap (map (map readDouble . words) . lines) body
  return rows

readDouble :: String -> Double
readDouble = read

listToTuple' :: [a] -> Maybe (a, a)
listToTuple' (x:y:xs) = Just (x, y)
listToTuple' _        = Nothing

-- | downloads the test data
testdata = input "http://work.caltech.edu/data/out.dta"

-- | downloads the training data
traindata = input "http://work.caltech.edu/data/in.dta"

testY = fmap (fmap makeVectorY) testdata
  where
    makeVectorY lists = vector $ map last lists

testX = makeDataMatrix testdata

trainY = fmap (fmap makeVectorY) traindata
  where
    makeVectorY lists = vector $ map last lists

trainX = makeDataMatrix traindata

makeLabelVector inputdata = fmap (fmap makeVectorY) inputdata
  where
    makeVectorY lists = vector $ map last lists

makeDataMatrix inputdata = fmap (fmap makeMatrixX) inputdata
  where
    makeMatrixX lists = matrix 3 $ concatMap xvectors lists
    xvectors list = 1 : init list

makeTransformedMatrix2 inputdata = fmap (fmap makeMatrixX) inputdata
  where
    makeMatrixX lists =
      createTransformedX2 $ map (fromMaybe (0, 0) . listToTuple') lists

makeTransformedMatrix k inputdata = fmap (fmap makeMatrixX) inputdata
  where
    makeMatrixX lists =
      createTransformedXk k $ map (fromMaybe (0, 0) . listToTuple') lists

trainXTransformed2 = makeTransformedMatrix2 traindata

testXTransformed2 = makeTransformedMatrix2 testdata

-- | Trains linear Regression from the data on the course website
trainLinRegUnkownTarget = do
  matrixXTrain <- trainXTransformed2
  vectorYTrain <- trainY
  let weights = linearRegressionWeight <$> matrixXTrain <*> vectorYTrain
  let inSampleError =
        linRegClassificationError <$> matrixXTrain <*> vectorYTrain <*> weights
  return (weights, inSampleError)

-- | Tests linear Regression from the data on the course website
testLinRegUnkownTarget weights = do
  matrixXTest <- testXTransformed2
  vectorYTest <- testY
  let outOfSampleError =
        linRegClassificationError <$> matrixXTest <*> vectorYTest <*> weights
  return outOfSampleError

-- | solution for question 2: train and test with given data, with nonlinear transformation but without regularization and print the result
trainAndTestWithData :: IO ()
trainAndTestWithData = do
  (weights, inSampleError) <- trainLinRegUnkownTarget
  outOfSampleError <- testLinRegUnkownTarget weights
  putStr "In Sample Error: "
  print inSampleError
  putStr "Out of Sample Error: "
  print outOfSampleError

-- | returns the weights determined by linear regression on the training data with regularization
linearRegressionWeightReg :: R -> Matrix R -> Vector R -> Vector R
linearRegressionWeightReg regParam matrixOfInputData vectorOfLabels =
  inv
    (tr matrixOfInputData HMatrix.<> matrixOfInputData +
     scalar regParam * (ident $fst $size $tr matrixOfInputData)) HMatrix.#>
  (tr matrixOfInputData #> vectorOfLabels)

-- | Trains linear Regression from the data on the course website with regularization
trainLinRegWithRegularization regParam = do
  matrixXTrain <- trainXTransformed2
  vectorYTrain <- trainY
  let weights =
        linearRegressionWeightReg regParam <$> matrixXTrain <*> vectorYTrain
  let inSampleError =
        linRegClassificationError <$> matrixXTrain <*> vectorYTrain <*> weights
  return (weights, inSampleError)

-- | performs training and testing with the regularization parameter provided as an input and prints the in-sample and out-of-sample error
trainAndTestWithRegularization :: R -> IO ()
trainAndTestWithRegularization regParam = do
  (weights, inSampleError) <- trainLinRegWithRegularization regParam
  outOfSampleError <- testLinRegUnkownTarget weights
  putStr "In Sample Error: "
  print inSampleError
  putStr "Out of Sample Error: "
  print outOfSampleError


text:
{-# LANGUAGE DataKinds         #-}
{-# LANGUAGE BangPatterns      #-}
{-# LANGUAGE RankNTypes        #-}
{-# LANGUAGE FlexibleContexts  #-}
{-# LANGUAGE TypeFamilies      #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE TypeOperators     #-}
-- {-# LANGUAGE UndecidableInstances #-}

module GrenadeExtras (
  epochTraining,
  trainOnBatchesEpochs,
  binaryNetError,
  normalize,
  hotvector,
  applyUpdates
) where

import           Shuffle (shuffle)
import           Control.Monad.Random (MonadRandom)

import           Data.List (foldl')

import           Lens.Micro (over, both)--Lens', set)

import           Data.Singletons (SingI)
import           Data.Singletons.Prelude (Head, Last)

import           Numeric.LinearAlgebra.Static (extract, R, ℝ, toRows, eye)
import           Numeric.LinearAlgebra.Data (Vector, (!))

import           GHC.TypeLits (KnownNat)

import           Grenade (Network((:~>), NNil), Shape(D1), LearningParameters, runNet, S(S1D),
                          Gradients((:/>)), backPropagate, train, runUpdates, Gradient, applyUpdate)

import           GrenadeExtras.OrphanNum()
import           GrenadeExtras.GradNorm (GradNorm(normSquared))

trainOnBatchesEpochs :: (SingI (Last shapes), MonadRandom m, Num (Gradients layers), GradNorm (Gradients layers))
                 => Network layers shapes
                 -> LearningParameters
                 -> [(S (Head shapes), S (Last shapes))]
                 -> Int
                 -> m [(Double, Network layers shapes)]
trainOnBatchesEpochs net0 rate input_data batchSize =

    foldMeOutList (0, net0) [(1::Int)..] $ \(_,net) _-> do
      shuffledInput <- shuffle input_data
      -- traning net (an epoch) with the input shuffled
      let batches = splitInBatches shuffledInput
          (gradientNorm, newNet) = foldl' trainBatch (0, net) batches
      return (gradientNorm, newNet)

  where
    trainBatch :: (SingI (Last shapes), Num (Gradients layers), GradNorm (Gradients layers))
               => (Double, Network layers shapes)
               -> [(S (Head shapes), S (Last shapes))]
               -> (Double, Network layers shapes)
    trainBatch (accNorm, !network) ios =
      let grads = fmap (uncurry $ backPropagate network) ios
          grad = sum grads
          norm = sqrt $ normSquared grad
       in (accNorm+norm, applyUpdate rate network grad)
       --in applyUpdates rate network grads

    --len = length input_data

    splitInBatches :: [a] -> [[a]]
    splitInBatches [] = []
    splitInBatches xs =
      let (start, finish) = splitAt batchSize xs
       in start : splitInBatches finish

applyUpdates :: LearningParameters
             -> Network layers shapes
             -> [Gradients layers]
             -> Network layers shapes
applyUpdates rate (layer :~> rest) gradients
  = runUpdates rate layer layerGradients :~> applyUpdates rate rest restLayersGradients
    where headTailGrad :: Gradients (layer ': layers) -> (Gradient layer, Gradients layers)
          headTailGrad (gradient :/> grest) = (gradient, grest)
          (layerGradients, restLayersGradients) = unzip $ fmap headTailGrad gradients

applyUpdates _ NNil _
  = NNil

epochTraining :: (SingI (Last shapes), MonadRandom m) =>
                 Network layers shapes
                 -> LearningParameters
                 -> [(S (Head shapes), S (Last shapes))]
                 -> m [Network layers shapes]
epochTraining net0 rate input_data =

    foldMeOutList net0 [(1::Int)..] $ \net _-> do
      shuffledInput <- shuffle input_data
      -- traning net (an epoch) with the input shuffled
      let newNet = foldl' trainEach net shuffledInput
      return newNet

  where
    trainEach !network (i,o) = train rate network i o

foldMeOutList :: Monad m => a -> [b] -> (a -> b -> m a) -> m [a]
foldMeOutList z xs_ op = f' z xs_
  where f' _    []     = return []
        f' zero (x:xs) = do
          zero' <- op zero x
          rec   <- f' zero' xs
          return $ zero':rec

binaryNetError :: (Last shapes ~ 'D1 1, Foldable t) =>
  Network layers shapes -> t (S (Head shapes), S ('D1 1)) -> (Double, Double)

binaryNetError net test = (fromIntegral errors / fromIntegral total, distance)
  where
    total, errors :: Integer
    {-(total, errors) = foldl' step (0,0) test-}
    distance :: Double
    (total, errors, distance) = foldl' step (0,0,0) test

    {-step (t, e) song = (t+1, e')-}
    step (t, e, d) song = (t+1, e', d')
      where
        (label, netOut) = valueFromDataAndNet song
        -- the predictions from the network come with numbers between 0 and 1,
        -- everything above .5 is considered 1 and below 0
        e' = if (label > 0.5) == (netOut > 0.5)
                then e
                else e+1
        d' = d + (label - netOut)^(2::Int)

    valueFromDataAndNet (input, label) = over both sd1toDouble (label, runNet net input)
      where
        sd1toDouble :: S ('D1 1) -> Double
        sd1toDouble (S1D r) = (extract :: R 1 -> Vector ℝ) r ! 0

-- taken from https://en.wikipedia.org/wiki/Normalization_(statistics)
-- normalization method: Student's t-statistic
normalize :: KnownNat n => [R n] -> [R n]
normalize features = fmap (\x -> (x-mean)/stdDeviation ) features
  where
    len = length features

    --mean :: R n
    mean = sum features / fromIntegral len

    --stdDeviation :: R n
    stdDeviation = sqrt $ (sum . fmap (\x-> (x-mean)^(2::Int)) $ features)
                          / fromIntegral (len-1)

hotvector :: KnownNat n => Int -> Maybe (R n)
hotvector m = toRows eye ~!! m
  where
    (~!!) :: [a] -> Int -> Maybe a
    []     ~!! _ = Nothing
    (x:xs) ~!! n
      | n == 0 = Just x
      | n < 0  = Nothing
      | otherwise = xs ~!! (n-1)

max_stars_repo_name: helq/haskell-binary-classification

max_issues_repo_name: helq/haskell-binary-classification

max_forks_repo_name: helq/haskell-binary-classification


text:
{-# LANGUAGE GADTs #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# OPTIONS_GHC -Wall #-}
-- This file takes 30+ seconds to compile with pattern match warnings
-- in GHC 8.4. https://gitlab.haskell.org/ghc/ghc/-/issues/14987
{-# OPTIONS_GHC -Wno-overlapping-patterns -Wno-incomplete-record-updates #-}
{-# OPTIONS_GHC -Wno-incomplete-patterns -Wno-incomplete-uni-patterns #-}
{-# OPTIONS_GHC -Wno-inaccessible-code #-}

--
-- Copyright (c) 2019, ERICSSON AB
-- All rights reserved.
--
-- Redistribution and use in source and binary forms, with or without
-- modification, are permitted provided that the following conditions are met:
--
--     * Redistributions of source code must retain the above copyright notice,
--       this list of conditions and the following disclaimer.
--     * Redistributions in binary form must reproduce the above copyright
--       notice, this list of conditions and the following disclaimer in the
--       documentation and/or other materials provided with the distribution.
--     * Neither the name of the ERICSSON AB nor the names of its contributors
--       may be used to endorse or promote products derived from this software
--       without specific prior written permission.
--
-- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-- AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
-- DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
-- FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
-- DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
-- SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
-- CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
-- OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
-- OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
--

-- | Translate from the typed to an untyped functional representation
module Feldspar.Core.Middleend.FromTyped
  ( toU
  ) where

import qualified Feldspar.Core.UntypedRepresentation as U
import qualified Feldspar.Core.NestedTuples as T
import Feldspar.Core.Types (TypeRep(..), TypeF(..), (:>)(..))
import qualified Feldspar.Core.Types as T
import Feldspar.Core.UntypedRepresentation hiding (Type(..), ScalarType(..))
import Feldspar.Core.ValueInfo (ValueInfo(..), ValueKind(..))
import Feldspar.Range (Range(..))
import qualified Feldspar.Core.Representation as R
import Feldspar.Core.Representation (AExpr(..), Expr(..))
import Data.Complex (Complex(..))
import Feldspar.Lattice (Lattice, universal)

toU :: AExpr a -> UntypedFeld ValueInfo
toU (((R.Info i) :: R.Info a) :& e)
  | Sym (R.Variable (R.Var n s)) <- e
  = i2 $ Variable $ Var n (untypeType tr i) s
  | Sym (R.Literal v) <- e
  = i2 $ Literal $ literal tr i v
  | Sym op <- e
  = i2 $ App (trOp op) (untypeType tr i) []
  | Sym (R.Lambda (R.Var n s)) :@ e' <- e
  , FunType b _ <- tr
  = i2 $ Lambda (Var n (untypeType b (fst i)) s) $ toU e'
  | Sym R.Cons :@ a1 :@ a2 <- e
  , e' <- toU a1
  , In _ (App Tup (U.TupType ts) es) <- toU a2
  = i2 $ App Tup (U.TupType $ typeof e' : ts) $ e' : es
  | Sym R.Car :@ a <- e
  , In _ (App (Drop n) (U.TupType (t:_)) es) <- addDrop $ toU a
  = i2 $ App (Sel n) t es
  | Sym R.Cdr :@ a <- e
  , In _ (App (Drop n) (U.TupType (_:ts)) es) <- addDrop $ toU a
  = i2 $ App (Drop $ n + 1) (U.TupType ts) es
  | Sym R.Tup :@ a <- e
  , In _ e' <- toU a
  = i2 e'
  | (op, es) <- go e []
  = i2 $ App op (untypeType tr i) es
  where tr = typeRepF :: TypeRep a
        i2 = In $ toValueInfo tr i
        go :: forall a' . Expr a' -> [UntypedFeld ValueInfo] -> (Op, [UntypedFeld ValueInfo])
        go (Sym op) es = (trOp op, es)
        go (f :@ e') es = go f $ toU e' : es
        addDrop e'@(In _ (App (Drop _) _ _)) = e'
        addDrop e' = In (error "FromTyped: temporary drop")
                         (App (Drop 0) (typeof e') [e'])

-- | Translate a Typed operator to the corresponding untyped one
trOp :: R.Op a -> Op
trOp R.GetLength       = GetLength
trOp R.Parallel        = Parallel
trOp R.Append          = Append
trOp R.GetIx           = GetIx
trOp R.SetLength       = SetLength
trOp R.Sequential      = Sequential
trOp R.SetIx           = SetIx
trOp R.Let             = Let
trOp R.Bit             = Bit
trOp R.Complement      = Complement
trOp R.ReverseBits     = ReverseBits
trOp R.BitScan         = BitScan
trOp R.BitCount        = BitCount
trOp R.BAnd            = BAnd
trOp R.BOr             = BOr
trOp R.BXor            = BXor
trOp R.SetBit          = SetBit
trOp R.ClearBit        = ClearBit
trOp R.ComplementBit   = ComplementBit
trOp R.TestBit         = TestBit
trOp R.ShiftLU         = ShiftLU
trOp R.ShiftRU         = ShiftRU
trOp R.ShiftL          = ShiftL
trOp R.ShiftR          = ShiftR
trOp R.RotateLU        = RotateLU
trOp R.RotateRU        = RotateRU
trOp R.RotateL         = RotateL
trOp R.RotateR         = RotateR
trOp R.RealPart        = RealPart
trOp R.ImagPart        = ImagPart
trOp R.Conjugate       = Conjugate
trOp R.Magnitude       = Magnitude
trOp R.Phase           = Phase
trOp R.Cis             = Cis
trOp R.MkComplex       = MkComplex
trOp R.MkPolar         = MkPolar
trOp R.Condition       = Condition
trOp R.ConditionM      = ConditionM
trOp R.F2I             = F2I
trOp R.I2N             = I2N
trOp R.B2I             = B2I
trOp R.Round           = Round
trOp R.Ceiling         = Ceiling
trOp R.Floor           = Floor
trOp R.ESkip           = ESkip
trOp R.EMaterialize    = EMaterialize
trOp R.EWrite          = EWrite
trOp R.EPar            = EPar
trOp R.EparFor         = EparFor
trOp R.Equal           = Equal
trOp R.NotEqual        = NotEqual
trOp R.Undefined       = Undefined
trOp (R.Assert s)      = Assert s
trOp (R.ForeignImport s _) = ForeignImport s
trOp R.Exp             = Exp
trOp R.Sqrt            = Sqrt
trOp R.Log             = Log
trOp R.Sin             = Sin
trOp R.Tan             = Tan
trOp R.Cos             = Cos
trOp R.Asin            = Asin
trOp R.Atan            = Atan
trOp R.Acos            = Acos
trOp R.Sinh            = Sinh
trOp R.Tanh            = Tanh
trOp R.Cosh            = Cosh
trOp R.Asinh           = Asinh
trOp R.Atanh           = Atanh
trOp R.Acosh           = Acosh
trOp R.Pow             = Pow
trOp R.LogBase         = LogBase
trOp R.Pi              = Pi
trOp R.DivFrac         = DivFrac
trOp R.MkFuture        = MkFuture
trOp R.Await           = Await
trOp R.Quot            = Quot
trOp R.Rem             = Rem
trOp R.Div             = Div
trOp R.Mod             = Mod
trOp R.IExp            = IExp
trOp R.Not             = Not
trOp R.And             = And
trOp R.Or              = Or
trOp R.ForLoop         = ForLoop
trOp R.WhileLoop       = WhileLoop
trOp R.While           = While
trOp R.For             = For
trOp R.Run             = Run
trOp R.Return          = Return
trOp R.Bind            = Bind
trOp R.Then            = Then
trOp R.When            = When
trOp R.NewArr_         = NewArr_
trOp R.ArrLength       = ArrLength
trOp R.NewArr          = NewArr
trOp R.GetArr          = GetArr
trOp R.SetArr          = SetArr
trOp R.RunMutableArray = RunMutableArray
trOp R.WithArray       = WithArray
trOp R.NewRef          = NewRef
trOp R.GetRef          = GetRef
trOp R.SetRef          = SetRef
trOp R.ModRef          = ModRef
trOp R.Nil             = Tup
trOp R.NoInline        = NoInline
trOp R.Abs             = Abs
trOp R.Sign            = Sign
trOp R.Add             = Add
trOp R.Sub             = Sub
trOp R.Mul             = Mul
trOp R.ParRun          = ParRun
trOp R.ParGet          = ParGet
trOp R.ParFork         = ParFork
trOp R.ParNew          = ParNew
trOp R.ParYield        = ParYield
trOp R.ParPut          = ParPut
trOp R.LTH             = LTH
trOp R.GTH             = GTH
trOp R.LTE             = LTE
trOp R.GTE             = GTE
trOp R.Min             = Min
trOp R.Max             = Max
trOp R.Atan2           = Atan2
trOp R.Save            = Save
trOp (R.PropSize _)    = PropSize
-- trOp R.SourceInfo      = SourceInfo
trOp R.Switch          = Switch
-- trOp R.Call            = Call
trOp (R.Variable _)    = error "trOp: Variable"
trOp (R.Lambda _)      = error "trOp: Lambda"
trOp (R.Literal _)     = error "trOp: Literal"
trOp R.Cons            = error "trOp: Cons"
trOp R.Car             = error "trOp: Car"
trOp R.Cdr             = error "trOp: Cdr"
trOp R.Tup             = error "trOp: Tup"

-- | Convert a BitWidth to an untyped Size
convSize :: T.BitWidth a -> U.Size
convSize T.N8      = U.S8
convSize T.N16     = U.S16
convSize T.N32     = U.S32
convSize T.N64     = U.S64

-- | Default size for numeric types
defaultNumSize :: Lattice (Range a) => TypeRep a -> T.Size a
defaultNumSize IntType{} = universal
defaultNumSize FloatType = universal
defaultNumSize DoubleType = universal
defaultNumSize t = error $ "defaultNumSize: missing numeric type: " ++ show t

untypeType :: TypeRep a -> T.Size a -> U.Type
untypeType UnitType _               = U.TupType []
untypeType BoolType _               = 1 U.:# U.BoolType
untypeType (IntType s n) _          = 1 U.:# U.IntType (convSign s) (convSize n)
untypeType FloatType _              = 1 U.:# U.FloatType
untypeType DoubleType _             = 1 U.:# U.DoubleType
untypeType (ComplexType t) _        = 1 U.:# U.ComplexType
  (case t of
    FloatType -> untypeType t (defaultNumSize t)
    DoubleType -> untypeType t (defaultNumSize t)
    _ -> error $ "untypeType: Missing complex type:" ++ show t)
untypeType (Tup2Type t) sz          = U.TupType $ untypeTup t sz
untypeType (Tup3Type t) sz          = U.TupType $ untypeTup t sz
untypeType (Tup4Type t) sz          = U.TupType $ untypeTup t sz
untypeType (Tup5Type t) sz          = U.TupType $ untypeTup t sz
untypeType (Tup6Type t) sz          = U.TupType $ untypeTup t sz
untypeType (Tup7Type t) sz          = U.TupType $ untypeTup t sz
untypeType (Tup8Type t) sz          = U.TupType $ untypeTup t sz
untypeType (Tup9Type t) sz          = U.TupType $ untypeTup t sz
untypeType (Tup10Type t) sz         = U.TupType $ untypeTup t sz
untypeType (Tup11Type t) sz         = U.TupType $ untypeTup t sz
untypeType (Tup12Type t) sz         = U.TupType $ untypeTup t sz
untypeType (Tup13Type t) sz         = U.TupType $ untypeTup t sz
untypeType (Tup14Type t) sz         = U.TupType $ untypeTup t sz
untypeType (Tup15Type t) sz         = U.TupType $ untypeTup t sz
untypeType (MutType a) sz              = U.MutType (untypeType a sz)
untypeType (RefType a) sz              = U.RefType (untypeType a sz)
untypeType (ArrayType a) (rs :> es)    = U.ArrayType rs (untypeType a es)
untypeType (MArrType a) (rs :> es)     = U.MArrType rs (untypeType a es)
untypeType (ParType a) sz              = U.ParType (untypeType a sz)
untypeType (ElementsType a) (_ :> es)  = U.ElementsType (untypeType a es)
untypeType c@ConsType{} sz             = U.TupType $ untypeTup c sz
untypeType NilType        _            = U.TupType []
untypeType (IVarType a) sz             = U.IVarType $ untypeType a sz
untypeType (FunType a b) (sa, sz)      = U.FunType (untypeType a sa) (untypeType b sz)
untypeType (FValType a) sz             = U.FValType (untypeType a sz)

untypeTup :: TypeRep (T.Tuple a) -> T.Size (T.Tuple a) -> [U.Type]
untypeTup (ConsType a b) (sa,sb) = untypeType a sa : untypeTup b sb
untypeTup NilType        _       = []

convSign :: T.Signedness a -> U.Signedness
convSign T.U       = Unsigned
convSign T.S       = Signed

literal :: TypeRep a -> T.Size a -> a -> Lit
literal UnitType        _ () = LTup []
literal BoolType        _  a = LBool a
literal (IntType s n)   _  a = LInt (convSign s) (convSize n) (toInteger a)
literal FloatType       _  a = LFloat a
literal DoubleType      _  a = LDouble a
literal (ArrayType t) (_ :> sz) a = LArray t' $ map (literal t sz) a
  where t' = untypeType t sz
literal (ComplexType t) _ (r :+ i)
  = case t of
     FloatType -> LComplex (literal t (defaultNumSize t) r)
                           (literal t (defaultNumSize t) i)
     DoubleType -> LComplex (literal t (defaultNumSize t) r)
                            (literal t (defaultNumSize t) i)
     _ -> error $ "literal: Missing complex type: " ++ show t
literal _ _ _ = error "Missing pattern: FromTyped.hs: literal"

-- | Construct a ValueInfo from a TypeRep and a Size
toValueInfo :: TypeRep a -> T.Size a -> ValueInfo
toValueInfo UnitType          _             = VIProd ProdKind []
-- FIXME: No range for boolean types yet.
toValueInfo BoolType                r       = VIBool    r
toValueInfo (IntType T.U T.N8)      r       = VIWord8   r
toValueInfo (IntType T.S T.N8)      r       = VIInt8    r
toValueInfo (IntType T.U T.N16)     r       = VIWord16  r
toValueInfo (IntType T.S T.N16)     r       = VIInt16   r
toValueInfo (IntType T.U T.N32)     r       = VIWord32  r
toValueInfo (IntType T.S T.N32)     r       = VIInt32   r
toValueInfo (IntType T.U T.N64)     r       = VIWord64  r
toValueInfo (IntType T.S T.N64)     r       = VIInt64   r
-- FIXME: No range for FP types and ComplexType yet.
toValueInfo FloatType         _             = VIFloat
toValueInfo DoubleType        _             = VIDouble
toValueInfo (ComplexType _)   _             = VIProd ProdKind []
toValueInfo (Tup2Type t) sz                 = toValueInfo t sz
toValueInfo (Tup3Type t) sz                 = toValueInfo t sz
toValueInfo (Tup4Type t) sz                 = toValueInfo t sz
toValueInfo (Tup5Type t) sz                 = toValueInfo t sz
toValueInfo (Tup6Type t) sz                 = toValueInfo t sz
toValueInfo (Tup7Type t) sz                 = toValueInfo t sz
toValueInfo (Tup8Type t) sz                 = toValueInfo t sz
toValueInfo (Tup9Type t) sz                 = toValueInfo t sz
toValueInfo (Tup10Type t) sz                = toValueInfo t sz
toValueInfo (Tup11Type t) sz                = toValueInfo t sz
toValueInfo (Tup12Type t) sz                = toValueInfo t sz
toValueInfo (Tup13Type t) sz                = toValueInfo t sz
toValueInfo (Tup14Type t) sz                = toValueInfo t sz
toValueInfo (Tup15Type t) sz                = toValueInfo t sz
toValueInfo (MutType a) sz                  = toValueInfo a sz
toValueInfo (RefType a) sz                  = toValueInfo a sz
toValueInfo (ArrayType a) (Range l r :> es)
  = VIProd RangeKind [VIWord32 (Range l r), toValueInfo a es]
toValueInfo (MArrType a) (Range l r :> es)
  = VIProd RangeKind [VIWord32 (Range l r), toValueInfo a es]
toValueInfo (ParType a) sz                  = toValueInfo a sz
toValueInfo (ElementsType a) (Range l r :> es)
  = VIProd RangeKind [VIWord32 (Range l r), toValueInfo a es]
toValueInfo (ConsType a b) (sa,sb) = VIProd TupleKind $ toValueInfo a sa : ss
  where VIProd _ ss = toValueInfo b sb
toValueInfo NilType _ = VIProd TupleKind []
toValueInfo (IVarType a) sz                 = toValueInfo a sz
-- TODO: Maybe keep argument information for FunType.
toValueInfo (FunType a _) (sa, _)           = toValueInfo a sa
toValueInfo (FValType a) sz                 = toValueInfo a sz

max_stars_repo_name: Feldspar/feldspar-language

max_issues_repo_name: Feldspar/feldspar-language

max_forks_repo_name: Feldspar/feldspar-language


text:
{- |
Module:      : ML.CostFunctions
Description  : Cost function derivatives for use in machine learning algorithms.
-}
module ML.CostFunctions
       (
           -- * Data types
           CostDerivative

           -- * Cost function derivatives
       ,   mse') where

import Numeric.LinearAlgebra (R, Vector)

-- | A vectorized function which returns ∂Cₓ/∂a.
--
-- The first parameter is the output activation, the second parameter
-- is the expected output.
type CostDerivative = Vector R -> Vector R -> Vector R

-- | The derivative of the mean squared error cost function.
mse' :: CostDerivative
mse' outputActivations y = outputActivations - y

max_stars_repo_name: m-renaud/ML

max_issues_repo_name: m-renaud/ML

max_forks_repo_name: m-renaud/ML


text:
import Control.Parallel

import Numeric.Container hiding (linspace)
import Numeric.LinearAlgebra.Util
import Numeric.LinearAlgebra.Data hiding (linspace)
import Numeric.LinearAlgebra.Algorithms
import Numeric.LinearAlgebra.Devel
import Data.List hiding (find)
import DSP.Basic

--gridsize
n = 75

--index function
index :: (Int, Int) -> Int -> Double
index (i,j) n
	| i == j       =  6
	| i == j + n*n = -1
	| i == j - n*n = -1
	| i == j + n   = -1
	| i == j - n   = -1
	| i == j + 1   = -1
	| i == j - 1   = -1
	| otherwise    =  0

-- Laplace Matrix in 3 dimensions
t = buildMatrix (n*n*n) (n*n*n) (\(i,j) -> index (i,j) n)
eigen = eigenvalues t

--Calculate corresponding Eigenfunction to n-th energy-eigenvalue
psi :: Int -> Vector Double
psi x = nullVector (a x)
	where a x = t - (mapMatrix (* (realPart (eigen @> x))) (ident (n*n*n)) )

--For further image processing
v = asColumn (psi 56)

main :: IO ()
main = do
-- write out matrix to visualize with python
saveMatrix "vmat.txt" "%f" v

max_stars_repo_name: davidschlegel/hartree-fock

max_issues_repo_name: davidschlegel/hartree-fock

max_forks_repo_name: davidschlegel/hartree-fock


text:
{-# LANGUAGE FlexibleInstances #-}
{-|
Module      : Math
Description : Auxiliary math functions
Copyright   : (c) Julian Kopka Larsen, 2015
Stability   : experimental

Definition of math functions like 'scatterMatrix' and the multivariate log Gamma function 'mlgamma'
-}
{-# LANGUAGE TypeSynonymInstances #-}
module Math (
    X,
    dimension,
    Sstat,
    add_means,
    sub_means,
    add_scatter,
    sub_scatter,
    fromX,
    AbelianGroup((<+>), (<->)),
    mlgamma,
    scatterMatrix,
    scatterMatrixAlt,
    mean,
    meanv,
    lnDeterminant
) where

import Numeric.LinearAlgebra (Matrix, Vector, trans, (<>), fromRows, dim, meanCov, invlndet, outer, konst, scale)
import Numeric.LinearAlgebra.Data (size)
import Debug.Trace
import Data.List ((\\))
import Partition (group)

qtr x = trace ("value: " ++ show x) x

-- | Class for the
class AbelianGroup a where
    (<+>) :: a -> a -> a
    (<->) :: a -> a -> a

-- | Type synonym for data objects that can be anything as long as it is compatible with the likelihood
--   could be extended to a type class with sufficient statistics member functions.
type X = [Vector Double]

-- | Computes the scatter matrix defined as:
scatterMatrixAlt :: X -> Matrix Double
scatterMatrixAlt x = sum $ map (\v -> outer (v-xm) (v-xm)) x
                where xm = meanv x

scatterMatrix :: X -> Matrix Double
scatterMatrix x = sum (map (\ v -> outer v v) x) - scale n (outer xm xm)
                    where xm = meanv x
                          n = fromIntegral $ length x

-- | Mean of a dataset
mean :: [Vector Double] -> Matrix Double
mean x = fromRows [meanv x]

meanv :: [Vector Double] -> Vector Double
meanv x = scale (1/(fromIntegral $ length x))(sum x)

dimension :: [Vector Double] -> Int
dimension = dim . head

type Sstat = (Int, Matrix Double, Vector Double)

fromX :: X -> Sstat
fromX x = (length x, scatterMatrix x, meanv x)

instance AbelianGroup Sstat where
    s <+> (0,_,_) = s
    (n1, s1, m1) <+> (n2, s2, m2) = (n1+n2, add_scatter n1 s1 m1 n2 s2 m2 , add_means n1 m1 n2 m2)
    s <-> (0,_,_) = s
    (n1, s1, m1) <-> (n2, s2, m2) = (n1-n2, sub_scatter n1 s1 m1 n2 s2 m2, sub_means n1 m1 n2 m2)

instance (AbelianGroup a) => AbelianGroup [a] where
    x <+> y = map (\(a,b) -> a <+> b) $ zip x y
    x <-> y = map (\(a,b) -> a <-> b) $ zip x y

add_means :: Int -> Vector Double -> Int -> Vector Double -> Vector Double
add_means _ m1 0 _ = m1
add_means 0 _ _ m2 = m2
add_means i1 m1 i2 m2 = scale w1 m1 + scale w2 m2
        where n1 = fromIntegral i1
              n2 = fromIntegral i2
              w1 = (n1/(n1+n2))
              w2 = (n2/(n1+n2))

sub_means :: Int -> Vector Double -> Int -> Vector Double -> Vector Double
sub_means _ m1 0 _ = m1
sub_means i1 _ i2 _ | i1==i2 = meanv []
sub_means i1 m1 i2 m2 = scale w1 m1 - scale w2 m2
        where n1 = fromIntegral i1
              n2 = fromIntegral i2
              w1 = (n1/(n1-n2))
              w2 = (n2/(n1-n2))

add_scatter :: Int -> Matrix Double -> Vector Double -> Int -> Matrix Double -> Vector Double -> Matrix Double
add_scatter _ sa _ 0 _ _ = sa
add_scatter 0 _ _ _ sb _ = sb
add_scatter ia sa ma ib sb mb = sa + sb + (scale na xaa) + (scale nb xbb) - (scale n xx)
        where na = fromIntegral ia
              nb = fromIntegral ib
              xaa = outer ma ma
              xbb = outer mb mb
              mx = add_means ia ma ib mb
              xx = outer mx mx
              n = na + nb

sub_scatter :: Int -> Matrix Double -> Vector Double -> Int -> Matrix Double -> Vector Double -> Matrix Double
sub_scatter 0 _ _ n _ _ = error "Nothing to subtract from."
sub_scatter _ sa _ 0 _ _ = sa
sub_scatter ia sa ma ib sb mb | ia==ib = scatterMatrix [ma]
sub_scatter ia sa ma ib sb mb | ia-ib==1 = konst 0 (size sa)
sub_scatter ia sa ma ib sb mb = sa - sb + (scale na xaa) - (scale nb xbb) - (scale n xx)
        where na = fromIntegral ia
              nb = fromIntegral ib
              xaa = outer ma ma
              xbb = outer mb mb
              mx = sub_means ia ma ib mb
              xx = outer mx mx
              n = na - nb

lnDeterminant :: Matrix Double -> Double
lnDeterminant x = absLnDet * signDet
        where (_,(absLnDet,signDet)) = invlndet x

cofG :: [Double]
cofG = [76.18009172947146,-86.50532032941677,24.01409824083091,-1.231739572450155,0.001208650973866179,-0.000005395239384953]

serG :: Double
serG = 1.000000000190015

lgamma :: Double -> Double
lgamma xx = let tmp' = (xx+5.5) - (xx+0.5)*log(xx+5.5)
                ser' = serG + (sum $ zipWith (/) cofG [xx+1..])
            in -tmp' + log(2.5066282746310005 * ser' / xx)

mlgamma :: Double -> Double -> Double
mlgamma p a = (log pi)*(p*(p-1)/4) + sum (map (\j -> lgamma (a + ((1-j)/2))) [1..p])

max_stars_repo_name: juliankopkalarsen/FpStats

max_issues_repo_name: juliankopkalarsen/FpStats

max_forks_repo_name: juliankopkalarsen/FpStats


text:
{-# LANGUAGE FlexibleContexts #-}
-- | Mixed-radix FFT calculation.
--
-- Arbitrary input vector lengths are handled using a mixed-radix
-- Cooley-Tukey decimation in time algorithm with residual prime
-- length vectors being treated using Rader's algorithm or hand-coded
-- codelets for small primes.
module Numeric.FFT
       ( fft, ifft, fftWith, ifftWith
       , plan, planFromFactors, execute
       , Plan (..), Direction (..), BaseTransform (..)
       ) where

import Prelude hiding (length, map, sum, zipWith)
import Data.Vector.Generic
import Data.Complex

import Numeric.FFT.Types
import Numeric.FFT.Plan
import Numeric.FFT.Execute

-- | Forward FFT with embedded plan calculation.  For an input vector
-- /h/ of length /N/, with entries numbered from 0 to /N - 1/,
-- calculates the entries in /H/, the discrete Fourier transform of
-- /h/, as:
--
-- <<doc-formulae/fft-formula.svg>>
fft :: Vector v (Complex Double) =>
       v (Complex Double) -> IO (v (Complex Double))
fft xs = do
  p <- plan $ length xs
  return $ fftWith p xs

-- | Inverse FFT with embedded plan calculation.  For an input vector
-- /H/ of length /N/, with entries numbered from 0 to /N - 1/,
-- representing Fourier amplitudes of a signal, calculates the entries
-- in /h/, the inverse discrete Fourier transform of /H/, as:
--
-- <<doc-formulae/ifft-formula.svg>>
ifft :: Vector v (Complex Double) =>
        v (Complex Double) -> IO (v (Complex Double))
ifft xs = do
  p <- plan $ length xs
  return $ ifftWith p xs

-- | Forward FFT with pre-computed plan.
fftWith :: Vector v (Complex Double) =>
           Plan -> v (Complex Double) -> v (Complex Double)
fftWith p = convert . execute p Forward . convert

-- | Inverse FFT with pre-computed plan.
ifftWith :: Vector v (Complex Double) =>
            Plan -> v (Complex Double) -> v (Complex Double)
ifftWith p = convert . execute p Inverse . convert

max_stars_repo_name: ian-ross/arb-fft

max_issues_repo_name: ian-ross/arb-fft

max_forks_repo_name: ian-ross/arb-fft