text:
import Control.Applicative ((<$>), (<*>))
import Control.Monad (forM, forM_, replicateM)
import Control.Monad.Random (Rand, RandomGen, getRandom, runRand)
import Data.List (genericLength, transpose)
import Numeric.LinearAlgebra
import System.Random (mkStdGen)
import Text.Printf (printf)

-- Exercises 1 - 5
getX :: Matrix Double -> Matrix Double
getX = takeColumns 2

getY :: Matrix Double -> Vector Double
getY = head . toColumns . dropColumns 2

nonLinear :: Matrix Double -> Matrix Double
nonLinear = fromRows . map (fromList . nonLinearRow . toList) . toRows
  where
    nonLinearRow [x1, x2] = [1, x1, x2, x1^2, x2^2, x1 * x2, abs (x1 - x2), abs (x1 + x2)]

classificationError :: Vector Double -> Vector Double -> Double
classificationError v y = let sv = mapVector signum v
                              diffs = zipVectorWith (((fromIntegral . fromEnum).) . (/=)) sv y
                          in  foldVector (+) 0 diffs / fromIntegral (dim y)

regression :: Matrix Double -> Vector Double -> Vector Double
regression x y = let xt = trans x
                     xtx = multiply xt x
                     z = multiply (pinv xtx) xt
                 in  z `mXv` y

experiment :: Matrix Double -> Vector Double -> Matrix Double -> Vector Double -> Double
experiment xtrain ytrain xval yval = let z = regression xtrain ytrain
                                     in  classificationError (xval `mXv` z) yval
-- Exercise 6
minimumEstimator :: RandomGen g => Rand g Double
minimumEstimator = min <$> getRandom <*> getRandom

minimumExperiments :: RandomGen g => Rand g Double
minimumExperiments = replicateM iterations minimumEstimator >>= average
  where
    iterations = 10000
    average xs = return $ sum xs / genericLength xs

-- Exercise 7
type Point = (Double, Double)
type RegressionModel = Point -> Point -> (Double -> Double)

loo :: [a] -> [(a, [a])]
loo [] = []
loo (x:xs) = loo' [] x xs
  where
    loo' xs y [] = [(y, xs)]
    loo' xs y (z:zs) = (y, xs ++ z:zs) : loo' (xs ++ [y]) z zs

constantModel :: RegressionModel
constantModel (_, py) (_, qy) = const a
  where
    a = (qy + py) / 2

linearModel :: RegressionModel
linearModel (px, py) (qx, qy) x = m * x + b
  where
    m = dy / dx
    dy = py - qy
    dx = px - qx
    b = py - m * px

squaredError :: (Double -> Double) -> Point -> Double
squaredError f (x, y) = (f x - y)^2

trainAndValidate :: [Point] -> RegressionModel -> Double
trainAndValidate ps m = average
  where
    trainLOO (a, [p, q]) = squaredError (m p q) a
    errors = map trainLOO (loo ps)
    average = sum errors / genericLength errors

main = do
  inputMatrix <- fromFile "in.dta" (35, 3)
  testMatrix <- fromFile "out.dta" (250, 3)

  let dx = nonLinear (getX inputMatrix)
      dy = getY inputMatrix
      xtest = nonLinear (getX testMatrix)
      ytest = getY testMatrix
      xtrain = takeRows 25 dx
      ytrain = subVector 0 25 dy
      xval = dropRows 25 dx
      yval = subVector 25 10 dy

  putStr "1 & 2)"
  firstTest <- forM [3..7] $ \k ->
    let trim = takeColumns (k + 1)
        run = experiment (trim xtrain) ytrain
        e_val = run (trim xval) yval
        e_test = run (trim xtest) ytest
        str = printf "\tk = %d: E_val = %.5f, E_out = %.5f" k e_val e_test
    in putStrLn str >> return (e_val, e_test)

  putStr "\n3 & 4)"
  secondTest <- forM [1..7] $ \k ->
    let trim = takeColumns (k + 1)
        run = experiment (trim xval) yval
        e_val = run (trim xtrain) ytrain
        e_test = run (trim xtest) ytest
        str = printf "\tk = %d: E_val = %.5f, E_out = %.5f" k e_val e_test
    in putStrLn str >> return (e_val, e_test)

  putStrLn $ printf "\n5)\t(%.5f, %.5f)" (snd $ minimum firstTest) (snd $ minimum secondTest)

  let g = mkStdGen 12345
  let (expectedMinimum, g') = runRand minimumExperiments g
  putStrLn $ printf "\n6)\t%.8f" expectedMinimum

  let ps = map sqrt [sqrt 3 + 4,
                     sqrt 3 - 1,
                     9 + 4 * sqrt 6,
                     9 - sqrt 6]
      fixed = [(-1, 0), (1, 0)]
  let triples = map (\p -> (p, 1):fixed) ps
  putStr "\n7)"
  forM_ (zip ['a'..'d'] triples) $ \(label, triple) ->
    let f = trainAndValidate triple
        ec = f constantModel
        el = f linearModel
    in putStrLn $ printf "\t%c) Constant: %.5f, Linear: %.5f" label ec el

ext: hs


text:
module TransferFunction where

import Numeric.GSL.ODE
import Numeric.LinearAlgebra
-- import Numeric.LinearAlgebra (linspace, toLists)

type Time = Double
type Freq = Double

solve :: (Time -> [Double] -> [Double]) -> Time -> [Double] -> Time -> [Double]
solve xdot t0 x0 t1 = last . toLists $ odeSolve xdot x0 ts
  where
    ts = linspace 2 (t0, t1)

diffStepFunc t = if 0.0 <= t && t <= trf then 1 / trf else 0
  where
    trf = 0.01
diffSinFunc fsig t = 2*pi*cos(2*pi*fsig*t)

lpf1 :: Freq -> (Time -> Double) -> Time -> [Double] -> [Double]
lpf1 fp xdot t [x, y] = [sx, sy]
  where
    sx = xdot t
    sy = (x - y) / t1 -- Y/X = 1 / (1 + s.t1) <=> Y + t1.sY = X
    t1 = 1 / 2 / pi / fp

lpf2 :: (Freq, Freq) -> (Time -> Double) -> Time -> [Double] -> [Double]
lpf2 (fp1, fp2) xdot t [x, y, sy] = [sx, sy, ssy]
  where
    sx  = xdot t
    ssy = x - y - (t1 + t2) * sy
    t1  = 1 / 2 / pi / fp1
    t2  = 1 / 2 / pi / fp2
    --     H(s) = 1 / (1 + s/wp1) * (1 + s/wp2)
    -- <=> Y + (t1+t2).sY + t1.t2.ssY = X
    -- <=> t1.t2.ssY = X - Y - (t1+t2).sY

hpf1 :: Freq -> (Time -> Double) -> Time -> [Double] -> [Double]
-- true 1-zero tranfer function can't be solved
hpf1 fz xdot t [x, y] = [sx, sy]
  where
    sx  = xdot t
    sy  = (tz*sx - y) / tz
    tz  = 1 / 2 / pi / fz

laglead :: (Freq, Freq) -> (Time -> Double) -> Time -> [Double] -> [Double]
laglead (fp, fz) xdot t [x, y] = [sx, sy]
  where
    sx  = xdot t
    sy  = (x + tz*sx - y) / tp
    tp  = 1 / 2 / pi / fp
    tz  = 1 / 2 / pi / fz
    -- H(x) = (1 + tz.s) / (1 + tp.s)
    -- <=> Y + tp.sY = X + tz.sX
    -- <=> sY = (X + tz.sX - Y) / tp

p2p :: [Double] -> Double
p2p xs = maximum xs - minimum xs


text:
{-# LANGUAGE NoMonomorphismRestriction #-}
-- Addapted from: "The MathematicalOrchid, 24 Feb 2007"
--  http://warp.povusers.org/MandScripts/haskell.html

import Data.Complex
import Diagrams.Prelude hiding (magnitude,image)
import Diagrams.Backend.Postscript.CmdLine

quadratic c z = z*z + c

critical_orbit z = iterate (quadratic z) 0

pixel = length . takeWhile (\z -> magnitude z <= 2) . take max_iter

side n v0 v1 =
   let sv = (v1 - v0) / fromIntegral n
   in  [v0, (v0 + sv) .. v1]

grid = map (\y -> map (:+ y) side_x) side_y

image = map (map (to_circle . pixel . critical_orbit)) grid

d :: Diagram Postscript
d = vcat . map hcat $ image

main = defaultMain d

side_x = side 64 (-2) 2
side_y = side 64 (-2) 2

max_iter = 32

minBox = centerXY (strutX w <> strutY w)
  where w = 2 * sqrt max_iter

to_circle :: Integral a => a -> Diagram Postscript
to_circle = (<> minBox) . centerXY . lwG 0 . fc blue . circle . sqrt . fromIntegral

max_stars_repo_name: diagrams/diagrams-postscript

max_issues_repo_name: diagrams/diagrams-postscript

max_forks_repo_name: diagrams/diagrams-postscript


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

module Backprop.Learn.Test (
  -- * Tests
    Test
  , maxIxTest, rmseTest
  , squaredErrorTest, absErrorTest, totalSquaredErrorTest, squaredErrorTestV
  , crossEntropyTest, crossEntropyTest1
  , boolTest
  -- ** Manipulate tests
  , lossTest, lmapTest
  -- * Run tests
  , testModel, testModelStoch, testModelAll, testModelStochAll
  -- ** Correlation tests
  , testModelCov, testModelCorr
  , testModelStochCov, testModelStochCorr
  ) where

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

-- TODO: support non-double results?

type Test o = o -> o -> Double

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

boolTest :: forall a. RealFrac a => Test a
boolTest x y
    | ri x == ri y = 1
    | otherwise    = 0
  where
    ri :: a -> Int
    ri = round

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

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

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

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

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

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

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

crossEntropyTest1 :: Test Double
crossEntropyTest1 targ res = -(log res * targ + log (1 - res) * (1 - targ))

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

testModel
    :: Test b
    -> Model p 'Nothing a b
    -> TMaybe p
    -> a
    -> b
    -> Double
testModel t f mp x y = t y $ runModelStateless f mp x

testModelStoch
    :: PrimMonad m
    => Test b
    -> Model p 'Nothing a b
    -> MWC.Gen (PrimState m)
    -> TMaybe p
    -> a
    -> b
    -> m Double
testModelStoch t f g mp x y = t y <$> runModelStochStateless f g mp x

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

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

testModelCov
    :: (Foldable t, Fractional b)
    => Model p 'Nothing a b
    -> TMaybe p
    -> t (a, b)
    -> b
testModelCov f p = L.fold $ (lmap . first) (runModelStateless f p) cov

testModelCorr
    :: (Foldable t, Floating b)
    => Model p 'Nothing a b
    -> TMaybe p
    -> t (a, b)
    -> b
testModelCorr f p = L.fold $ (lmap . first) (runModelStateless f p) corr

testModelAll
    :: Foldable t
    => Test b
    -> Model p 'Nothing a b
    -> TMaybe p
    -> t (a, b)
    -> Double
testModelAll t f p = L.fold $ lmap (uncurry (testModel t f p)) L.mean

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

testModelStochAll
    :: (Foldable t, PrimMonad m)
    => Test b
    -> Model p 'Nothing a b
    -> MWC.Gen (PrimState m)
    -> TMaybe p
    -> t (a, b)
    -> m Double
testModelStochAll t f g p = L.foldM $ L.premapM (uncurry (testModelStoch t f g p))
                                      (L.generalize L.mean)

testModelStochCov
    :: (Foldable t, PrimMonad m, Fractional b)
    => Model p 'Nothing a b
    -> MWC.Gen (PrimState m)
    -> TMaybe p
    -> t (a, b)
    -> m b
testModelStochCov f g p = L.foldM $ (L.premapM . flip bitraverse pure)
                                      (runModelStochStateless f g p)
                                      (L.generalize cov)

testModelStochCorr
    :: (Foldable t, PrimMonad m, Floating b)
    => Model p 'Nothing a b
    -> MWC.Gen (PrimState m)
    -> TMaybe p
    -> t (a, b)
    -> m b
testModelStochCorr f g p = L.foldM $ (L.premapM . flip bitraverse pure)
                                       (runModelStochStateless f g p)
                                       (L.generalize corr)

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
{-# LANGUAGE UnicodeSyntax
           , MultiParamTypeClasses
           , FunctionalDependencies
           , FlexibleInstances
           , FlexibleContexts
--           , OverlappingInstances
           , UndecidableInstances #-}

import Numeric.LinearAlgebra

class Scaling a b c | a b -> c where
 -- ^ 0x22C5	8901	DOT OPERATOR, scaling
 infixl 7 ⋅
 (⋅) :: a -> b -> c

class Contraction a b c | a b -> c where
 -- ^ 0x00D7	215	MULTIPLICATION SIGN	×, contraction
 infixl 7 ×
 (×) :: a -> b -> c

class Outer a b c | a b -> c where
 -- ^ 0x2297	8855	CIRCLED TIMES	⊗, outer product (not associative)
 infixl 7 ⊗
 (⊗) :: a -> b -> c

-------

instance (Num t) => Scaling t t t where
    (⋅) = (*)

instance Container Vector t => Scaling t (Vector t) (Vector t) where
    (⋅) = scale

instance Container Vector t => Scaling (Vector t) t (Vector t) where
    (⋅) = flip scale

instance Container Vector t => Scaling t (Matrix t) (Matrix t) where
    (⋅) = scale

instance Container Vector t => Scaling (Matrix t) t (Matrix t) where
    (⋅) = flip scale

instance Product t => Contraction (Vector t) (Vector t) t where
    (×) = dot

instance Product t => Contraction (Matrix t) (Vector t) (Vector t) where
    (×) = mXv

instance Product t => Contraction (Vector t) (Matrix t) (Vector t) where
    (×) = vXm

instance Product t => Contraction (Matrix t) (Matrix t) (Matrix t) where
    (×) = mXm

--instance Scaling a b c => Contraction a b c where
--    (×) = (⋅)

-----

instance Product t => Outer (Vector t) (Vector t) (Matrix t) where
    (⊗) = outer

instance Product t => Outer (Vector t) (Matrix t) (Matrix t) where
    v ⊗ m = kronecker (asColumn v) m

instance Product t => Outer (Matrix t) (Vector t) (Matrix t) where
    m ⊗ v = kronecker m (asRow v)

instance Product t => Outer (Matrix t) (Matrix t) (Matrix t) where
    (⊗) = kronecker

-----

v = 3 |> [1..] :: Vector Double

m = (3 >< 3) [1..] :: Matrix Double

s = 3 :: Double

a = s ⋅ v × m × m × v ⋅ s

b = (v ⊗ m) ⊗ (v ⊗ m)

c = v ⊗ m ⊗ v ⊗ m

d = s ⋅ (3 |> [10,20..] :: Vector Double)

main = do
    print $ scale s v <> m <.> v
    print $ scale s v <.> (m <> v)
    print $ s * (v <> m <.> v)
    print $ s ⋅ v × m × v
    print a
    print (b == c)
    print d

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
module Parser
    ( parseExpr
    , spaces
    ) where

import Text.ParserCombinators.Parsec hiding (spaces)
import Control.Monad
import Numeric
import Data.Ratio
import Data.Complex
import Data.Array
import Common (LispVal(..))

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

spaces :: Parser ()
spaces = skipMany1 space

parseString :: Parser LispVal
parseString = do
    char '"'
    x <- many $ escapeChars <|> noneOf "\"\\"
    char '"'
    return $ String x

escapeChars :: Parser Char
escapeChars = do
    char '\\' -- backslash
    x <- oneOf "\\\"nrt" -- either backslash or quotation mark
    return $ case x of -- return escaped character
        '\\' -> x
        '\"' -> x
        'n'  -> '\n'
        'r'  -> '\r'
        't'  -> '\t'

parseCharacter :: Parser LispVal
parseCharacter = do
    try $ string "#\\"
    value <- try (string "newline" <|> string "space")
             <|> do {x <- anyChar; notFollowedBy alphaNum; return [x] }
    return $ Character $ case value of
        "space" -> ' '
        "newline" -> '\n'
        otherwise -> (value !! 0)

parseAtom :: Parser LispVal
parseAtom = do
    first <- letter <|> symbol
    rest <- many (letter <|> digit <|> symbol)
    let atom = first:rest
    return $ Atom atom

parseBool :: Parser LispVal
parseBool = do
    char '#'
    (char 't' >> return (Bool True)) <|> (char 'f' >> return (Bool False))

parseNumber :: Parser LispVal
parseNumber = parseDecimal1 <|> parseDecimal2 <|> parseHex <|> parseOct <|> parseBin

parseDecimal1 :: Parser LispVal
parseDecimal1 = many1 digit >>= (return . Number . read)

parseDecimal2 :: Parser LispVal
parseDecimal2 = do  try $ string "#d"
                    x <- many1 digit
                    (return . Number . read) x

parseHex :: Parser LispVal
parseHex = do   try $ string "#x"
                x <- many1 hexDigit
                return $ Number (hex2dig x)

parseOct :: Parser LispVal
parseOct = do   try $ string "#o"
                x <- many1 octDigit
                return $ Number (oct2dig x)

parseBin :: Parser LispVal
parseBin = do   try $ string "#b"
                x <- many1 (oneOf "10")
                return $ Number (bin2dig x)

oct2dig x = fst $ readOct x !! 0
hex2dig x = fst $ readHex x !! 0
bin2dig = bin2dig' 0
bin2dig' digint "" = digint
bin2dig' digint (x:xs) = let old = 2 * digint + (if x == '0' then 0 else 1) in bin2dig' old xs

parseFloat :: Parser LispVal
parseFloat = do x <- many1 digit
                char '.'
                y <- many1 digit
                return $ Float (fst.head $ readFloat (x++"."++y))

parseRatio :: Parser LispVal
parseRatio = do x <- many1 digit
                char '/'
                y <- many1 digit
                return $ Ratio ((read x) % (read y))

toDouble :: LispVal -> Double
toDouble (Float f) = realToFrac f
toDouble (Number n) = fromIntegral n

parseComplex :: Parser LispVal
parseComplex = do   x <- (try parseFloat <|> parseDecimal1)
                    char '+'
                    y <- (try parseFloat <|> parseDecimal1)
                    char 'i'
                    return $ Complex (toDouble x :+ toDouble y)

parseExpr :: Parser LispVal
parseExpr = parseAtom
         <|> parseString
         <|> try parseComplex
         <|> try parseFloat
         <|> try parseRatio
         <|> try parseNumber -- try is there because of the hash char
         <|> try parseBool
         <|> try parseCharacter
         <|> parseQuoted
         <|> try (do    string "#("
                        x <- parseVector
                        char ')'
                        return x)
         <|> parseList
         <|> parseQuasiQuoted
         <|> parseUnQuote
         <|> parseUnQuoteSplicing

parseList :: Parser LispVal
parseList = between beg end parseList1
    where   beg = (char '(' >> skipMany space)
            end = (skipMany space >> char ')')

parseList1 :: Parser LispVal
parseList1 = do list <- sepEndBy parseExpr spaces
                maybeDatum <- optionMaybe (char '.' >> spaces >> parseExpr)
                return $ case maybeDatum of
                    Nothing -> List list
                    Just datum -> DottedList list datum

parseDottedList :: Parser LispVal
parseDottedList = do
    head <- endBy parseExpr spaces
    tail <- char '.' >> spaces >> parseExpr
    return $ DottedList head tail

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

parseQuasiQuoted :: Parser LispVal
parseQuasiQuoted = do
    char '`'
    x <- parseExpr
    return $ List [Atom "quasiquote", x]

parseUnQuote :: Parser LispVal
parseUnQuote = do
    char ','
    x <- parseExpr
    return $ List [Atom "unquote", x]

parseUnQuoteSplicing :: Parser LispVal
parseUnQuoteSplicing = do
    char ','
    char '@'
    x <- parseExpr
    return $ List [Atom "unquote-splicing", x]

parseVector :: Parser LispVal
parseVector = do    arrayValues <- sepBy parseExpr spaces
                    return $ Vector (listArray (0, (length arrayValues - 1)) arrayValues)


text:
{-|
Module: Numeric.Morpheus.Activation
Description: Activation Functions
Copyright: (c) Alexander Ignatyev, 2017
License: BSD-3
Stability: experimental
Portability: POSIX
-}

module Numeric.Morpheus.Activation
(
  sigmoid
  , sigmoidGradient
  , relu
  , reluGradient
  , tanh_
  , tanhGradient
)

where

import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Devel
import System.IO.Unsafe(unsafePerformIO)
import Foreign.C.Types
import Foreign.Ptr(Ptr)
import Foreign.Storable

createMatrixOfShape :: Storable a => Matrix a -> IO (Matrix a)
createMatrixOfShape m = createMatrix (orderOf m) (rows m) (cols m)

{- morpheus_sigmoid -}
foreign import ccall unsafe "morpheus_sigmoid"
  c_morpheus_sigmoid :: CInt -> Ptr Double -> Ptr Double -> IO ()

call_morpheus_sigmoid :: CInt -> CInt -> CInt -> CInt -> Ptr Double
                      -> CInt -> CInt -> CInt -> CInt -> Ptr Double
                      -> IO ()
call_morpheus_sigmoid rows cols _ _ xPtr _ _ _ _ yPtr = do
  c_morpheus_sigmoid (rows*cols) xPtr yPtr

-- | Calculates sigmoid.
sigmoid :: Matrix R -> Matrix R
sigmoid m = unsafePerformIO $ do
    y <- createMatrixOfShape m
    apply m (apply y id) call_morpheus_sigmoid
    return y

{- morpheus_sigmoid_gradient -}
foreign import ccall unsafe "morpheus_sigmoid_gradient"
  c_morpheus_sigmoid_gradient :: CInt -> Ptr Double -> Ptr Double -> IO ()

call_morpheus_sigmoid_gradient :: CInt -> CInt -> CInt -> CInt -> Ptr Double
                               -> CInt -> CInt -> CInt -> CInt -> Ptr Double
                               -> IO ()
call_morpheus_sigmoid_gradient rows cols _ _ xPtr _ _ _ _ yPtr = do
  c_morpheus_sigmoid_gradient (rows*cols) xPtr yPtr

-- | Calculates derivates of sigmoid.
sigmoidGradient :: Matrix R -> Matrix R
sigmoidGradient m = unsafePerformIO $ do
    y <- createMatrixOfShape m
    apply m (apply y id) call_morpheus_sigmoid_gradient
    return y

{- morpheus_relu -}
foreign import ccall unsafe "morpheus_relu"
  c_morpheus_relu :: CInt -> Ptr Double -> Ptr Double -> IO ()

call_morpheus_relu :: CInt -> CInt -> CInt -> CInt -> Ptr Double
                      -> CInt -> CInt -> CInt -> CInt -> Ptr Double
                      -> IO ()
call_morpheus_relu rows cols _ _ xPtr _ _ _ _ yPtr = do
  c_morpheus_relu (rows*cols) xPtr yPtr

-- | Calculates ReLu.
relu :: Matrix R -> Matrix R
relu m = unsafePerformIO $ do
    y <- createMatrixOfShape m
    apply m (apply y id) call_morpheus_relu
    return y

{- morpheus_relu_gradient -}
foreign import ccall unsafe "morpheus_relu_gradient"
  c_morpheus_relu_gradient :: CInt -> Ptr Double -> Ptr Double -> IO ()

call_morpheus_relu_gradient :: CInt -> CInt -> CInt -> CInt -> Ptr Double
                               -> CInt -> CInt -> CInt -> CInt -> Ptr Double
                               -> IO ()
call_morpheus_relu_gradient rows cols _ _ xPtr _ _ _ _ yPtr = do
  c_morpheus_relu_gradient (rows*cols) xPtr yPtr

-- | Calculates derivates of ReLu.
reluGradient :: Matrix R -> Matrix R
reluGradient m = unsafePerformIO $ do
    y <- createMatrixOfShape m
    apply m (apply y id) call_morpheus_relu_gradient
    return y

{- morpheus_tanh -}
foreign import ccall unsafe "morpheus_tanh"
  c_morpheus_tanh :: CInt -> Ptr Double -> Ptr Double -> IO ()

call_morpheus_tanh :: CInt -> CInt -> CInt -> CInt -> Ptr Double
                      -> CInt -> CInt -> CInt -> CInt -> Ptr Double
                      -> IO ()
call_morpheus_tanh rows cols _ _ xPtr _ _ _ _ yPtr = do
  c_morpheus_tanh (rows*cols) xPtr yPtr

-- | Calculates tanh.
tanh_ :: Matrix R -> Matrix R
tanh_ m = unsafePerformIO $ do
    y <- createMatrixOfShape m
    apply m (apply y id) call_morpheus_tanh
    return y

{- morpheus_tanh_gradient -}
foreign import ccall unsafe "morpheus_tanh_gradient"
  c_morpheus_tanh_gradient :: CInt -> Ptr Double -> Ptr Double -> IO ()

call_morpheus_tanh_gradient :: CInt -> CInt -> CInt -> CInt -> Ptr Double
                               -> CInt -> CInt -> CInt -> CInt -> Ptr Double
                               -> IO ()
call_morpheus_tanh_gradient rows cols _ _ xPtr _ _ _ _ yPtr = do
  c_morpheus_tanh_gradient (rows*cols) xPtr yPtr

-- | Calculates derivates of tanh.
tanhGradient :: Matrix R -> Matrix R
tanhGradient m = unsafePerformIO $ do
    y <- createMatrixOfShape m
    apply m (apply y id) call_morpheus_tanh_gradient
    return y

max_stars_repo_name: Alexander-Ignatyev/morpheus

max_issues_repo_name: aligusnet/morpheus

max_forks_repo_name: aligusnet/morpheus


text:
{-# LANGUAGE CPP                 #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE FlexibleInstances   #-}
{-# LANGUAGE ScopedTypeVariables #-}
#if __GLASGOW_HASKELL__ >= 706
{-# LANGUAGE PolyKinds           #-}
#endif
{-# OPTIONS_GHC -fno-warn-orphans -fno-warn-unused-imports #-}
module Data.Binary.Orphans () where

import Data.Binary
import Data.Binary.Put

import           Control.Applicative   (Alternative (..))
import           Control.Monad         (MonadPlus (..), liftM, liftM2)
import qualified Control.Monad.Fail    as Fail
import           Data.Bits             (Bits, shiftL, shiftR, (.|.))
import           Data.Complex          (Complex (..))
import qualified Data.Fixed            as Fixed
import           Data.Functor.Identity (Identity (..))
import           Data.List             (foldl', unfoldr)
import qualified Data.List.NonEmpty    as NE
import qualified Data.Monoid           as Monoid
import           Data.Semigroup        ((<>))
import qualified Data.Semigroup        as Semigroup
import           Data.Version          (Version (..))
import           Data.Void             (Void, absurd)
import           GHC.Fingerprint       (Fingerprint (..))
import           Numeric.Natural       (Natural)

-------------------------------------------------------------------------------
-- binary-0.7.1.0
-------------------------------------------------------------------------------

#if !(MIN_VERSION_binary(0,7,1))

-------------------------------------------------------------------------------
-- Add MonadPlus instance for Get.

-- Not implementable, as far as I can see

#endif

-------------------------------------------------------------------------------
-- binary-0.7.3.0
-------------------------------------------------------------------------------

#if !(MIN_VERSION_binary(0,7,3) )

-------------------------------------------------------------------------------
-- Add Binary instance for Natural (only with base > 4.8).

#ifndef MIN_VERSION_nats
#define MIN_VERSION_nats(x,y,z) 0
#endif

#if !(MIN_VERSION_nats(1,1,0))
-- Fixed-size type for a subset of Natural
type NaturalWord = Word64

-- | /Since: 0.7.3.0/
instance Binary Natural where
    {-# INLINE put #-}
    put n | n <= hi =
        putWord8 0
        >> put (fromIntegral n :: NaturalWord)  -- fast path
     where
        hi = fromIntegral (maxBound :: NaturalWord) :: Natural

    put n =
        putWord8 1
        >> put (unroll (abs n))         -- unroll the bytes

    {-# INLINE get #-}
    get = do
        tag <- get :: Get Word8
        case tag of
            0 -> liftM fromIntegral (get :: Get NaturalWord)
            _ -> do bytes <- get
                    return $! roll bytes

--
-- Fold and unfold an Integer to and from a list of its bytes
--
unroll :: (Integral a, Bits a) => a -> [Word8]
unroll = unfoldr step
  where
    step 0 = Nothing
    step i = Just (fromIntegral i, i `shiftR` 8)

roll :: (Integral a, Bits a) => [Word8] -> a
roll   = foldl' unstep 0 . reverse
  where
    unstep a b = a `shiftL` 8 .|. fromIntegral b
#endif

#endif

-------------------------------------------------------------------------------
-- binary-0.7.6.0
-------------------------------------------------------------------------------

#if !MIN_VERSION_binary(0,7,6)

-------------------------------------------------------------------------------
-- Added binary instance for GHC.Fingerprint (from GHC >= 7.4).

instance Binary Fingerprint where
    put (Fingerprint x1 x2) = put x1 <> put x2
    get = do
        x1 <- get
        x2 <- get
        return $! Fingerprint x1 x2

#endif

-------------------------------------------------------------------------------
-- binary-0.8.0.0
-------------------------------------------------------------------------------

#if !(MIN_VERSION_binary(0,8,0))

-------------------------------------------------------------------------------
-- Added binary instance for Version from Data.Version.

instance Binary Version where
    put (Version br tags) = put br <> put tags
    get = liftM2 Version get get

-------------------------------------------------------------------------------
-- Added binary instance for Void from GHC 7.10.1.

instance Binary Void where
    put     = absurd
    get     = fail "get @Void"

-------------------------------------------------------------------------------
-- Added binary instance for (Data.Fixed a) from GHC 7.8.1.

#if MIN_VERSION_base(4,7,0)
instance Binary (Fixed.Fixed a) where
  put (Fixed.MkFixed a) = put a
  get = Fixed.MkFixed `fmap` get
#else
instance Fixed.HasResolution a => Binary (Fixed.Fixed a) where
  -- Using undefined :: Maybe a as a proxy, as Data.Proxy is introduced only in base-4.7
  put x = put (truncate (x * fromInteger (Fixed.resolution (undefined :: Maybe a))) :: Integer)
  get = (\x -> fromInteger x / fromInteger (Fixed.resolution (undefined :: Maybe a))) `fmap` get
#endif

#endif

-------------------------------------------------------------------------------
-- Added semigroup instance for Data.Binary.Builder from GHC 8.0.

-- in semigroups

-------------------------------------------------------------------------------
-- binary-0.8.2.0
-------------------------------------------------------------------------------

-------------------------------------------------------------------------------
-- When using GHC >= 8, Data.Binary.Get.Get implements MonadFail and delegates its fail to MonadFail.fail.

#if !(MIN_VERSION_binary(0,8,2)) || !(MIN_VERSION_base(4,9,0))
instance Fail.MonadFail Get where
    -- this is ok, as if old base: Prelude.fail is Monad's fail
    fail = Prelude.fail
#endif

-------------------------------------------------------------------------------
-- binary-0.8.3
-------------------------------------------------------------------------------

-------------------------------------------------------------------------------
-- Add Monoid and Semigroup instance for Put.

#if !(MIN_VERSION_binary(0,8,3)) || !(MIN_VERSION_base(4,9,0))
instance Semigroup.Semigroup Put where
    (<>) = (>>)
#endif

#if !(MIN_VERSION_binary(0,8,3))
instance Monoid.Monoid Put where
    mempty = return ()
    mappend = (<>)

#endif

-------------------------------------------------------------------------------
-- Add Binary instance for Complex a.

#if !(MIN_VERSION_binary(0,8,3))

instance Binary a => Binary (Complex a) where
    {-# INLINE put #-}
    put (r :+ i) = put (r, i)
    {-# INLINE get #-}
    get = fmap (\(r,i) -> r :+ i) get

#endif

-------------------------------------------------------------------------------
-- binary-0.8.4
-------------------------------------------------------------------------------

-------------------------------------------------------------------------------
-- Add Binary instances for datatypes in Data.Monoid/Data.Semigroup
-- https://github.com/kolmodin/binary/pull/114

#if !MIN_VERSION_binary(0,8,4)

instance Binary a => Binary (Monoid.Dual a) where
  get = fmap Monoid.Dual get
  put = put . Monoid.getDual

instance Binary Monoid.All where
  get = fmap Monoid.All get
  put = put . Monoid.getAll

instance Binary Monoid.Any where
  get = fmap Monoid.Any get
  put = put . Monoid.getAny

instance Binary a => Binary (Monoid.Sum a) where
  get = fmap Monoid.Sum get
  put = put . Monoid.getSum

instance Binary a => Binary (Monoid.Product a) where
  get = fmap Monoid.Product get
  put = put . Monoid.getProduct

instance Binary a => Binary (Monoid.First a) where
  get = fmap Monoid.First get
  put = put . Monoid.getFirst

instance Binary a => Binary (Monoid.Last a) where
  get = fmap Monoid.Last get
  put = put . Monoid.getLast

#if MIN_VERSION_base(4,8,0)
instance Binary (f a) => Binary (Monoid.Alt f a) where
  get = fmap Monoid.Alt get
  put = put . Monoid.getAlt
#endif

#endif

------------------------------------------------------------------------
-- Data.Semigroup datatypes

#if !MIN_VERSION_binary(0,8,4) || !MIN_VERSION_base(4,9,0)

instance Binary a => Binary (Semigroup.Min a) where
  get = fmap Semigroup.Min get
  put = put . Semigroup.getMin

instance Binary a => Binary (Semigroup.Max a) where
  get = fmap Semigroup.Max get
  put = put . Semigroup.getMax

instance Binary a => Binary (Semigroup.First a) where
  get = fmap Semigroup.First get
  put = put . Semigroup.getFirst

instance Binary a => Binary (Semigroup.Last a) where
  get = fmap Semigroup.Last get
  put = put . Semigroup.getLast

instance Binary a => Binary (Semigroup.Option a) where
  get = fmap Semigroup.Option get
  put = put . Semigroup.getOption

instance Binary m => Binary (Semigroup.WrappedMonoid m) where
  get = fmap Semigroup.WrapMonoid get
  put = put . Semigroup.unwrapMonoid

instance (Binary a, Binary b) => Binary (Semigroup.Arg a b) where
  get                     = liftM2 Semigroup.Arg get get
  put (Semigroup.Arg a b) = put a <> put b

instance Binary a => Binary (NE.NonEmpty a) where
  get = do
      x <- get
      case x of
          []     -> fail "empty NonEmpty"
          (x:xs) -> return (x NE.:| xs)
  put = put . NE.toList

#endif

-------------------------------------------------------------------------------
-- binary-0.8.5.0
-------------------------------------------------------------------------------

#if !(MIN_VERSION_binary(8,6,0))

-------------------------------------------------------------------------------
-- Typeable TypeReps
-- Add Binary instances for Typeable TypeReps
-- https://github.com/kolmodin/binary/pull/131

-- TODO

#endif

-------------------------------------------------------------------------------
-- binary-0.8.6.0
-------------------------------------------------------------------------------

#if !(MIN_VERSION_binary(8,6,0))

-------------------------------------------------------------------------------
-- Add binary instance for Data.Functor.Identity from base,
-- https://github.com/kolmodin/binary/pull/146

#if !MIN_VERSION_base(4,8,0)
instance Binary a => Binary (Identity a) where
  put (Identity x) = put x
  get = fmap Identity get
#endif

#endif

max_stars_repo_name: haskellari/binary-orphans

max_issues_repo_name: haskellari/binary-orphans

max_forks_repo_name: haskellari/binary-orphans


text:
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE TypeApplications    #-}
-- |
-- Module      : Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Level2
-- Copyright   : [2017..2020] Trevor L. McDonell
-- License     : BSD3
--
-- Maintainer  : Trevor L. McDonell <trevor.mcdonell@gmail.com>
-- Stability   : experimental
-- Portability : non-portable (GHC extensions)
--

module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Level2
  where

import Data.Complex
import Data.Array.Accelerate.Representation.Array
import Data.Array.Accelerate.Representation.Shape
import Data.Array.Accelerate.Sugar.Elt

import Data.Array.Accelerate.LLVM.Native.Foreign
import Data.Array.Accelerate.Numeric.LinearAlgebra.Type
import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Base

import Foreign.Ptr
import qualified Blas.Primitive.Types                               as C
import qualified Blas.Primitive.Unsafe                              as C

gemv :: NumericR s e
     -> Transpose
     -> ForeignAcc ((((((), Scalar e), Matrix e), Vector e)) -> Vector e)
gemv nR opA = ForeignAcc "native.gemv" (gemv' nR opA)

gemv' :: NumericR s e
      -> Transpose
      -> ((((), Scalar e), Matrix e), Vector e)
      -> Par Native (Future (Vector e))
gemv' nR opA ((((), alpha), matA), vecx) = do
  let
      (((), rowsA), colsA) = shape matA

      sizeY   = case opA of
                  N -> rowsA
                  _ -> colsA

      opA'    = encodeTranspose opA
      alpha'  = indexArray (ArrayR dim0 eR) alpha ()

      aR      = ArrayR dim1 eR
      eR      = case nR of
                  NumericRfloat32   -> eltR @Float
                  NumericRfloat64   -> eltR @Double
                  NumericRcomplex32 -> eltR @(Complex Float)
                  NumericRcomplex64 -> eltR @(Complex Double)
  --
  future <- new
  vecy   <- allocateRemote aR ((), sizeY)
  ()     <- liftIO $ do
    withArray nR matA   $ \ptr_A -> do
     withArray nR vecx  $ \ptr_x -> do
      withArray nR vecy $ \ptr_y -> do
        case nR of
          NumericRfloat32   -> C.sgemv C.RowMajor opA' rowsA colsA alpha' ptr_A colsA ptr_x 1 0 ptr_y 1
          NumericRfloat64   -> C.dgemv C.RowMajor opA' rowsA colsA alpha' ptr_A colsA ptr_x 1 0 ptr_y 1
          NumericRcomplex32 -> C.cgemv C.RowMajor opA' rowsA colsA (toElt alpha') (castPtr ptr_A) colsA (castPtr ptr_x) 1 0 (castPtr ptr_y) 1
          NumericRcomplex64 -> C.zgemv C.RowMajor opA' rowsA colsA (toElt alpha') (castPtr ptr_A) colsA (castPtr ptr_x) 1 0 (castPtr ptr_y) 1
  --
  put future vecy
  return future

max_stars_repo_name: statusfailed/accelerate-blas

max_issues_repo_name: statusfailed/accelerate-blas

max_forks_repo_name: statusfailed/accelerate-blas


text:
{-# LANGUAGE FlexibleInstances, BangPatterns, Rank2Types, StandaloneDeriving #-}

module LinearTesting
    ( properties
    ) where

import Control.Applicative

import Data.List (transpose, intersperse)
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Storable as S
import qualified Data.Packed.Matrix as M
import Data.Traversable

import Debug.Trace

import Foreign.Storable

import Numeric (showFFloat)
import Numeric.AD
import Numeric.MaxEnt.Linear
import Numeric.LinearAlgebra

import Test.Tasty
import Test.Tasty.QuickCheck

--------------------------------------------------------------------------------
-- Test groups
--------------------------------------------------------------------------------

properties :: TestTree
properties = testGroup "Properties"
    [ testProperty "solvableSystemsAreSolvable" solvableSystemsAreSolvable
    , testProperty "probsSumToOne" probsSumToOne
    , testProperty "solutionFitsConstraints" solutionFitsConstraints
    ]

--------------------------------------------------------------------------------
-- QuickCheck and helper datatypes, instances, etc.
--------------------------------------------------------------------------------

class Approximate a where
    (=~=) :: a -> a -> Bool

data InvalidLinearConstraints a = InvalidLinearConstraints
  { imatrix :: [[a]]
  , ioutput :: [a]
  } deriving (Eq, Show)

instance Arbitrary (InvalidLinearConstraints Double) where
    arbitrary = sized $ \size -> do

        matrix <- vectorOf size (vector size)
        unnormalizedProbs <- vector size

        badSum <- suchThat arbitrary (/= 1.0)

        let probs    = let total = sum unnormalizedProbs in
                                map (\x -> badSum * (x/total)) unnormalizedProbs
            hmatrix     = M.fromLists matrix
            hprobs      = M.fromLists $ transpose [probs]
            inputVector = hmatrix `multiply` hprobs

        return . InvalidLinearConstraints matrix . (map head) . M.toLists $ inputVector

data ValidLinearConstraints = VLC
  { unVLC :: forall a. Floating a => ([[a]], [a]) }

toLC :: ValidLinearConstraints -> LinearConstraints
toLC vlc = LC $ unVLC vlc

deriving instance Eq ValidLinearConstraints

instance Approximate Double where
    x =~= y = abs (x - y) < 0.1

instance Approximate a => Approximate [a] where
    xs =~= ys = all id . zipWith (=~=) xs $ ys

instance (Approximate a, Storable a) => Approximate (S.Vector a) where
    xs =~= ys = S.all id . S.zipWith (=~=) xs $ ys

instance Approximate (ValidLinearConstraints) where
    (VLC (mx, ox)) =~= (VLC (my, oy)) =
        (mx :: [[Double]]) =~= my && (ox :: [Double]) =~= oy

traceIt x = trace (show x) x

printRow :: [Double] -> String
printRow xs = "[" ++
   concat (intersperse "," (map (\x -> showFFloat (Just 6) x "") xs)) ++ "]"

instance Show (ValidLinearConstraints) where
   show (VLC (xss, xs)) = "matrix = [" ++
       concat (intersperse "," $ map printRow xss) ++ "]\n output = " ++ printRow xs

normalize xs = let total = sum xs in map (/total) xs

instance Arbitrary (ValidLinearConstraints) where
    arbitrary = sized $ \size' -> do

        let size = size' + 1

        matrix <- vectorOf size (vectorOf size (suchThat arbitrary (>0.0))) :: Gen [[Double]]
        unnormalizedProbs <- vectorOf size (suchThat arbitrary (>0.0)) :: Gen [Double]

        let probs       = normalize unnormalizedProbs
            --matrix' :: (Floating a) => [a]
            matrix'     = map normalize matrix
            hmatrix     = M.fromLists matrix'
            hprobs      = M.fromLists $ transpose [probs]
            inputVector = hmatrix `multiply` hprobs
            output = map head $ M.toLists inputVector

        return $ VLC (map (map (fromRational . toRational)) matrix', map (fromRational . toRational) output)

--
--toPoly :: RealFloat a => LinearConstraints Double -> LinearConstraints a
--toPoly (LC x y) =
--     LC (map (map (fromRational . toRational)) x)
--                        (map (fromRational . toRational) y)

--------------------------------------------------------------------------------
-- QuickCheck properties
--------------------------------------------------------------------------------

solvableSystemsAreSolvable :: ValidLinearConstraints -> Bool
solvableSystemsAreSolvable vlc =
    case linear 0.0000005 (toLC vlc) of
        Right _ -> True
        Left  _ -> False

traceItNote msg x = trace (msg ++ " " ++ show x) x

probsSumToOne :: ValidLinearConstraints -> Bool
probsSumToOne vlc =
    case linear 0.000005 (toLC vlc) of
        Right ps -> case 1.0 =~= S.sum ps of
            True -> True
            False -> trace ("new probs" ++ show ps) False
        Left _   -> False

solutionFitsConstraints :: ValidLinearConstraints -> Bool
solutionFitsConstraints vlc = let lc = toLC vlc in
    case linear 0.000005 lc of
        Right ps -> result where
            (x, y) = unLC lc
            result = ((map head) . M.toLists $ inputVector) =~= y

            hmatrix     = M.fromLists x
            hprobs      = M.fromLists $ transpose [S.toList ps]
            inputVector = hmatrix <> hprobs

        Left _   -> False

--------------------------------------------------------------------------------
-- Unfinished stuff. After finishing, remember to move to the appropriate place
-- and add to an exported test group.
--------------------------------------------------------------------------------

--This is not the test I want ..  but it does seem to work
-- TODO you want to test against the original probs which should
-- have less then or equal to the estimated
entropyIsGreaterOrEqual :: ValidLinearConstraints -> Bool
entropyIsGreaterOrEqual vlc = let lc = toLC vlc in
    case linear 0.000005 lc of
        Right ps -> result where
            entropy xs = negate . sum . map (\x -> x * log x) $ xs

            yEntropy         = entropy . snd $ unLC lc
            estimatedEntropy = entropy $ S.toList ps

            result = yEntropy >= estimatedEntropy
        Left  _  -> error "failed!"

--probabilityNeighborhood ps =

-- TODO make this
entropyIsMaximum :: ValidLinearConstraints -> Bool
entropyIsMaximum vlc = let lc = toLC vlc in
    case linear 0.000005 lc of
        Right ps -> result where
            entropy xs = negate . sum . map (\x -> x * log x) $ xs

            yEntropy         = entropy . snd $ unLC lc
            estimatedEntropy = entropy $ S.toList ps

            result = yEntropy >= estimatedEntropy
        Left  _  -> error "failed!"

-- also if it is a maximum a small change in either direction that still fits the constraints
-- should not lower the entropy
--

max_stars_repo_name: jfischoff/maxent

max_issues_repo_name: jfischoff/maxent

max_forks_repo_name: jfischoff/maxent


text:
{-# LANGUAGE AllowAmbiguousTypes   #-}
{-# LANGUAGE ConstraintKinds       #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE InstanceSigs          #-}
{-# LANGUAGE LambdaCase            #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE StandaloneDeriving    #-}
{-# LANGUAGE TypeApplications      #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeInType            #-}
{-# LANGUAGE TypeOperators         #-}

module Learn.Neural.Layer.Mapping (
    Mapping
  , CommonMap(..)
  , IdentMap, ScaleMap, LogitMap, ReLUMap, ReLUpMap, ELUMap, ELUpMap
  , MapFunc(..)
  , PMapping
  , CommonPMap(..)
  , PReLUMap, PELUMap
  , PMapFunc(..)
  ) where

import           Control.Applicative
import           Data.Kind
import           Data.Proxy
import           Data.Reflection
import           Data.Singletons
import           Data.Traversable
import           Data.Type.Combinator
import           Data.Type.Conjunction
import           Data.Type.Util
import           Data.Type.Vector hiding         (head')
import           GHC.Generics                    (Generic)
import           Learn.Neural.Layer
import           Numeric.BLAS
import           Numeric.Backprop
import           Numeric.Backprop.Iso
import           Numeric.Backprop.Op
import           Statistics.Distribution
import           Statistics.Distribution.Uniform
import           Type.Class.Higher
import           Type.Class.Known
import           Type.Class.Witness
import qualified Data.Type.Nat                   as TCN
import qualified Generics.SOP                    as SOP
import qualified Type.Family.Nat                 as TCN

data Mapping :: k -> Type

newtype MapFunc :: Type where
    MF :: { runMapFunc :: (forall a. RealFloat a => a -> a)
          }
       -> MapFunc

instance Num (CParam (Mapping s) b i i) where
    _ + _         = MapP
    _ * _         = MapP
    _ - _         = MapP
    negate _      = MapP
    abs    _      = MapP
    signum _      = MapP
    fromInteger _ = MapP

instance Fractional (CParam (Mapping s) b i i) where
    _ / _          = MapP
    recip _        = MapP
    fromRational _ = MapP

instance Floating (CParam (Mapping s) b i i) where
    sqrt _ = MapP

instance Num (CState (Mapping s) b i i) where
    _ + _         = MapS
    _ * _         = MapS
    _ - _         = MapS
    negate _      = MapS
    abs    _      = MapS
    signum _      = MapS
    fromInteger _ = MapS

instance Fractional (CState (Mapping s) b i i) where
    _ / _          = MapS
    recip _        = MapS
    fromRational _ = MapS

instance Floating (CState (Mapping s) b i i) where
    sqrt _ = MapS

instance (BLAS b, Reifies s MapFunc, SingI i) => Component (Mapping s) b i i where
    data CParam (Mapping s) b i i = MapP
    data CState (Mapping s) b i i = MapS
    data CConf  (Mapping s) b i i = MapC

    componentOp = componentOpDefault

    initParam _ _ _ _ = return MapP
    initState _ _ _ _ = return MapS
    defConf           = MapC

instance (BLAS b, Reifies s MapFunc, SingI i) => ComponentFF (Mapping s) b i i where
    componentOpFF = bpOp . withInps $ \(x :< _ :< Ø) -> do
        y <- tmapOp (runMapFunc mf) ~$ (x :< Ø)
        return . only $ y
      where
        mf :: MapFunc
        mf = reflect (Proxy @s)

instance (BLAS b, Reifies s MapFunc, SingI i) => ComponentLayer r (Mapping s) b i i where
    componentRunMode = RMIsFF

data CommonMap :: Type where
    MF_Ident :: CommonMap
    MF_Scale :: a -> CommonMap
    MF_Logit :: CommonMap
    MF_Tanh  :: CommonMap
    MF_ReLU  :: CommonMap
    MF_ReLUp :: a -> CommonMap
    MF_ELU   :: CommonMap
    MF_ELUp  :: a -> CommonMap

instance Reifies 'MF_Ident MapFunc where
    reflect _ = MF id
instance Reifies s Double => Reifies ('MF_Scale s) MapFunc where
    reflect _ = MF $ (realToFrac α *)
      where
        α :: Double
        α = reflect (Proxy @s)
instance Reifies 'MF_Logit MapFunc where
    reflect _ = MF $ \x -> 1 / (1 + exp (-x))
instance Reifies 'MF_Tanh MapFunc where
    reflect _ = MF tanh
instance Reifies 'MF_ReLU MapFunc where
    reflect _ = MF $ \x -> if x < 0 then 0 else x
instance Reifies s Double => Reifies ('MF_ReLUp s) MapFunc where
    reflect _ = MF $ \x -> if x < 0 then realToFrac α * x else x
      where
        α :: Double
        α = reflect (Proxy @s)
instance Reifies 'MF_ELU MapFunc where
    reflect _ = MF $ \x -> if x < 0 then exp x - 1 else x
instance Reifies s Double => Reifies ('MF_ELUp s) MapFunc where
    reflect _ = MF $ \x -> if x < 0 then realToFrac α * (exp x - 1) else x
      where
        α :: Double
        α = reflect (Proxy @s)

type IdentMap   = Mapping 'MF_Ident
type ScaleMap s = Mapping ('MF_Scale s)
type LogitMap   = Mapping 'MF_Logit
type TanhMap    = Mapping 'MF_Tanh
type ReLUMap    = Mapping 'MF_ReLU
type ReLUpMap s = Mapping ('MF_ReLUp s)
type ELUMap     = Mapping 'MF_ELU
type ELUpMap s  = Mapping ('MF_ELUp s)

data PMapping :: k -> TCN.N -> Type

data PMapFunc :: TCN.N -> Type where
    PMF :: { runPMapFunc :: (forall a. RealFloat a => (I :&: Vec n) a -> a)
           , getPMapDef  :: Vec n (SomeC ContGen I)
           }
        -> PMapFunc n

instance (Tensor b, Known TCN.Nat n) => Num (CParam (PMapping s n) b i i) where
    PMapP x + PMapP y = PMapP $ x + y
    PMapP x - PMapP y = PMapP $ x - y
    PMapP x * PMapP y = PMapP $ x * y
    negate (PMapP x) = PMapP (negate x)
    abs    (PMapP x) = PMapP (abs    x)
    signum (PMapP x) = PMapP (signum x)
    fromInteger x         = PMapP (fromInteger x)

-- until Vec has a Fractional instance
instance (Tensor b, Known TCN.Nat n) => Fractional (CParam (PMapping s n) b i i) where
    PMapP x / PMapP y = PMapP $ liftA2 (/) x y
    recip (PMapP x)   = PMapP $ fmap recip x
    fromRational x    = PMapP $ pure (fromRational x)

instance (Tensor b, Known TCN.Nat n) => Floating (CParam (PMapping s n) b i i) where
    sqrt (PMapP x)   = PMapP $ fmap sqrt x

instance Num (CState (PMapping s n) b i i) where
    _ + _         = PMapS
    _ * _         = PMapS
    _ - _         = PMapS
    negate _      = PMapS
    abs    _      = PMapS
    signum _      = PMapS
    fromInteger _ = PMapS

instance Fractional (CState (PMapping s n) b i i) where
    _ / _          = PMapS
    recip _        = PMapS
    fromRational _ = PMapS

instance Floating (CState (PMapping s n) b i i) where
    sqrt _        = PMapS

pMapP :: Known TCN.Nat n => Iso' (CParam (PMapping s n) b i i) (Tuple (Replicate n (Scalar b)))
pMapP = gTuple . iso (vecToProd . getI . head') (only_ . prodToVec' known)

deriving instance Generic (CParam (PMapping s n) b i i)
instance SOP.Generic (CParam (PMapping s n) b i i)

instance (BLAS b, Reifies s (PMapFunc n), SingI i, Known TCN.Nat n) => Component (PMapping s n) b i i where
    data CParam (PMapping s n) b i i = PMapP { _getPMapP :: !(Vec n (Scalar b)) }
    data CState (PMapping s n) b i i = PMapS
    data CConf  (PMapping s n) b i i = PMapC { _getPMapC :: !(Vec n (SomeC ContGen I)) }

    componentOp = componentOpDefault

    initParam _ _ c g = do
        ps <- forM (_getPMapC c) $ \(SomeC (I d)) ->
          realToFrac <$> genContVar d g
        return $ PMapP ps
    initState _ _ _ _ = return PMapS

    defConf = PMapC (getPMapDef pmf)
      where
        pmf :: PMapFunc n
        pmf = reflect (Proxy @s)

instance (BLAS b, Reifies s (PMapFunc n), SingI i, Known TCN.Nat n) => ComponentFF (PMapping s n) b i i where
    componentOpFF
        :: forall q. Num (b i)
        => OpB q '[b i, CParam (PMapping s n) b i i] '[b i]
    componentOpFF = bpOp . withInps $ \(x :< mp :< Ø) -> replWit @n @Num @(Scalar b) n Wit //
                                                         replLen @n @(Scalar b) n          // do
        ps :: Prod (BVar q '[b i, CParam (PMapping s n) b i i]) (Replicate n (Scalar b))
          <- replWit @n @Num @(Scalar b) n Wit //
             replLen @n @(Scalar b) n //
               pMapP #<~ mp
        let psV :: VecT n (BVar q '[b i, CParam (PMapping s n) b i i]) (Scalar b)
            psV = prodToVec' n ps
        psTV :: VecT n (BVar q '[b i, CParam (PMapping s n) b i i]) (b i)
          <- vtraverse (\p -> tkonstOp sing ~$ (p :< Ø)) psV
        let psT :: Prod (BVar q '[b i, CParam (PMapping s n) b i i]) (Replicate n (b i))
            psT = vecToProd psTV
        y <- tzipNOp @_ @_ @('TCN.S n) (\case x' :* psT' -> runPMapFunc pmf (x' :&: psT'))
               ~$ (x :< psT)
        return $ only y
      where
        n :: TCN.Nat n
        n = known
        pmf :: PMapFunc n
        pmf = reflect (Proxy @s)

instance (BLAS b, Reifies s (PMapFunc n), SingI i, Known TCN.Nat n) => ComponentLayer r (PMapping s n) b i i where
    componentRunMode = RMIsFF

data CommonPMap :: Type where
    PMF_PReLU :: CommonPMap
    PMF_PELU  :: CommonPMap

instance Reifies 'PMF_PReLU (PMapFunc TCN.N1) where
    reflect _ = PMF (\case I x :&: (I α :* ØV) ->
                            if x < 0 then α * x else x)
                    (SomeC (I (uniformDistr 0 0.01)) :+ ØV)

instance Reifies 'PMF_PELU (PMapFunc TCN.N1) where
    reflect _ = PMF (\case I x :&: (I α :* ØV) ->
                            if x < 0 then α * (exp x - 1) else x)
                    (SomeC (I (uniformDistr 0 0.01)) :+ ØV)

type PReLUMap = PMapping 'PMF_PReLU
type PELUMap  = PMapping 'PMF_PELU


text:
module StatisticsTests
  ( statsSpecs
  )
where

import           Control.Exception   (evaluate)
import           Data.Complex        (Complex)
import qualified Data.HashMap.Strict as HM
import           Data.List           (nub, sort)
import           Data.Statistics
import           Test.Invariant
import           Test.QuickCheck
import           Test.Tasty
import           Test.Tasty.Hspec

statsSpecs :: IO TestTree
statsSpecs = testSpec "Statistics module specs" $ do
  describe "mean spec"        meanSpec
  describe "median spec"      medianSpec
  describe "counter spec"     counterSpec
  describe "quantile spec"    quantileSpec
  describe "mode spec"        modeSpec
  describe "deMean spec"      deMeanSpec
  describe "variance spec"    varianceSpec
  describe "variance spec"    stdDevSpec
  describe "correlation spec" correlationSpec

meanSpec :: Spec
meanSpec = do
  it "returns 0 for empty list" $ do
    mean ([] :: [Double]) `shouldBe` 0
  it "computes mean correctly" $ do
    mean [1 :: Double, 2, 3, 4] `shouldBe` 2.5
  it "mean of single element list is element" $ property $ \x ->
    mean [x] == (x :: Double)
  it "mean is always between min and max" $ property meanInValueRange

counterSpec :: Spec
counterSpec = do
  it "number of entries in counter equals unique elements" $ do
    property $ \xs ->
      let distinct    = sort $ nub (xs :: [Int])
          counterKeys = sort $ HM.keys $ counter id xs
      in  distinct == counterKeys
  it "sum of counted values equals input size" $ do
    property countedAllElements

medianSpec :: Spec
medianSpec = do
  it "throws exception on empty list" $ do
    evaluate (median ([] :: [Double])) `shouldThrow` anyException
  it "uses average of two inner elements on list with even length"
    $          median [3 :: Rational, 8, 17, 33]
    `shouldBe` 12.5
  it "is indifferent to sorting" $ property $ \xs ->
    (not $ null xs) ==> median xs `shouldBe` median (sort (xs :: [Double]))
  it "single element list should have that median" $ property $ \x ->
    median [x] == (x :: Double)
  it "median is always within the input range" $ property medianInValueRange
  it "median and mean are the same for two-element lists" $ property $ \x y ->
    mean [x, y] == median [x, (y :: Double)]

quantileSpec :: Spec
quantileSpec = do
  it "throws exception on empty list"
    $             evaluate (quantile ([] :: [Double]) 0)
    `shouldThrow` anyException
  it "produces the only value in a singleton list"
    $ forAll (choose (0, 1))
    $ \p x -> quantile [x] (p :: Double) == (x :: Int)
  it "produces largest element for p=1" $ property $ \xs ->
    (not $ null xs) ==> quantile xs 1 == maximum (xs :: [Double])
  it "produces smarllest element for p=0" $ property $ \xs ->
    (not $ null xs) ==> quantile xs 0 == minimum (xs :: [Double])
  it "is monotonically increasing in p" $ forAll (choose (0, 1)) $ \p1 ->
    forAll (choose (0, 1)) $ \p2 xs ->
      (not $ null xs) ==> monotonicIncreasing (quantile (xs :: [Double])) p1 p2

modeSpec :: Spec
modeSpec = do
  it "finds multiple modes"
    $          mode [1, 2, 2, 4, 4, 1, 3]
    `shouldBe` [1, 2, (4 :: Int)]

deMeanSpec :: Spec
deMeanSpec = do
  it "does not change the number of elements" $ property $ \xs ->
    (not $ null xs) ==> length (deMean xs) == length (xs :: [Double])
  it "doesn't alter the data range" $ property $ \xs ->
    (not $ null xs) ==> dataRange (deMean xs) `doubleEq` dataRange
      (xs :: [Double])
  it "resulting data has mean of zero" $ property $ \xs ->
    (not $ null xs) ==> mean (deMean xs) `doubleEq` (0 :: Double)

varianceSpec :: Spec
varianceSpec = do
  it "is insensitive to shifting" $ property $ \xs y ->
    (length xs >= 2) ==> variance xs `doubleEq` variance
      ((+ (y :: Double)) <$> xs)

-- | compare double values, accepting a small absolute and relative deviation
-- Equality is assumed if the relative *or* the absolute deviation is small
doubleEq' :: Double -- ^ epsilon: maximum relative and absolute deviation
          -> Double
          -> Double
          -> Bool
doubleEq' e x y =
  let mx = max x y
      mn = min x y
  in  mx - mn < e || mx / mn < (1 + e)

-- | the double comparison used in the local tests
doubleEq :: Double -> Double -> Bool
doubleEq = doubleEq' 0.0000000001

stdDevSpec :: Spec
stdDevSpec = do
  it "is insensitive to shifting" $
    property $ \xs y -> (length xs >= 2) ==>
      stdDev xs `doubleEq` stdDev ((+ (y :: Double)) <$> xs)
  it "scales with a scaling factor" $
    property $ \xs y -> (length xs >= 2) ==>
      (y * stdDev xs) `doubleEq` (stdDev ((* (y :: Double)) <$> xs))

correlationSpec :: Spec
correlationSpec = do
  it "is 1 for self-correlation" $
    property $ \xs -> (length xs >= 2) ==>
      correlation (xs :: [Double]) xs `doubleEq` 1

countedAllElements :: [Complex Double] -> Bool
countedAllElements xs = let c = counter id xs in sum (HM.elems c) == length xs

meanInValueRange :: [Double] -> Property
meanInValueRange xs =
  length xs >=  1
    ==> let minV = minimum xs
            maxV = maximum xs
            res  = mean xs
        in  res >= minV && res <= maxV

medianInValueRange :: [Double] -> Property
medianInValueRange xs =
  length xs >=  1
    ==> let minV = minimum xs
            maxV = maximum xs
        in  median xs >= minV && median xs <= maxV


text:
{-# LANGUAGE CPP                   #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DeriveAnyClass        #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE UndecidableInstances  #-}
{-|
Module      : Grenade.Layers.Reshape
Description : Multipurpose reshaping layer
Copyright   : (c) Huw Campbell, 2016-2017
License     : BSD2
Stability   : experimental
-}
module Grenade.Layers.Reshape
  ( Reshape(..)
  , SpecReshape (..)
  , specReshape
  , specReshape3D1D
  , specReshape3D2D
  , specReshape2D1D
  , specReshape2D3D
  , specReshape1D2D
  , specReshape1D3D
  , reshape
  ) where

#if MIN_VERSION_singletons(2,6,0)
import           Data.Singletons.TypeLits       (SNat (..))
#endif

import           Control.DeepSeq                (NFData (..))
import           Control.Monad                  (when)
import           Data.Constraint                (Dict (..))
import           Data.Reflection                (reifyNat)
import           Data.Serialize
import           Data.Singletons
import           Data.Singletons.Prelude.Num    ((%*))
import qualified Data.Vector.Storable           as V
import           GHC.Generics                   (Generic)
import           GHC.TypeLits
import           Numeric.LinearAlgebra.Data     as LA (flatten)
import           Numeric.LinearAlgebra.Static   hiding (toRows)
import           Unsafe.Coerce                  (unsafeCoerce)

import           Grenade.Core
import           Grenade.Dynamic
import           Grenade.Dynamic.Internal.Build
import           Grenade.Layers.Trivial
import           Grenade.Utils.Conversion

import           Debug.Trace

-- | Reshape Layer
--
-- The Reshape layer can flatten any 2D or 3D image to 1D vector with the
-- same number of activations, as well as cast up from 1D to a 2D or 3D
-- shape.
--
-- Can also be used to turn a 3D image with only one channel into a 2D image
-- or vice versa.
data Reshape = Reshape
  deriving (Show,Generic,NFData)

instance UpdateLayer Reshape where
  type Gradient Reshape = ()
  runUpdate _ _ _ = Reshape

instance RandomLayer Reshape where
  createRandomWith _ _ = return Reshape

instance (KnownNat a, KnownNat x, KnownNat y, a ~ (x * y)) => Layer Reshape ('D2 x y) ('D1 a) where
  type Tape Reshape ('D2 x y) ('D1 a) = ()
  runForwards _ (S2D y)  =  ((), fromJust' . fromStorable . flatten . extract $ y)
  runForwards _ (S2DV y) = ((), fromStorableV y)
  runBackwards _ _ (S1D y)  = ((), fromJust' . fromStorable . extract $ y)
  runBackwards _ _ (S1DV y) = ((), fromStorableV y)

instance (KnownNat a, KnownNat x, KnownNat y, KnownNat (x * z), KnownNat z, a ~ (x * y * z)) => Layer Reshape ('D3 x y z) ('D1 a) where
  type Tape Reshape ('D3 x y z) ('D1 a) = ()
  runForwards _ (S3D y)     = ((), fromJust' . fromStorable . flatten . extract $ y)
  runBackwards _ _ (S1D y)  = ((), fromJust' . fromStorable . extract $ y)
  runBackwards _ _ (S1DV _) = error "not yet implemented"

instance (KnownNat y, KnownNat x, KnownNat z, z ~ 1) => Layer Reshape ('D3 x y z) ('D2 x y) where
  type Tape Reshape ('D3 x y z) ('D2 x y) = ()
  runForwards _ (S3D y)    = ((), S2D y)
  runBackwards _ _ (S2D y) = ((), S3D y)
  runBackwards _ _ _       = error "S3DV not yet implemented"

instance (KnownNat y, KnownNat x, KnownNat z, z ~ 1) => Layer Reshape ('D2 x y) ('D3 x y z) where
  type Tape Reshape ('D2 x y) ('D3 x y z) = ()
  runForwards _ (S2D y)  = ((), S3D y)
  runForwards _ (S2DV _) = error "S3DV not yet implemented"
  runBackwards _ _ (S3D y) = ((), S2D y)

instance (KnownNat a, KnownNat x, KnownNat y, a ~ (x * y)) => Layer Reshape ('D1 a) ('D2 x y) where
  type Tape Reshape ('D1 a) ('D2 x y) = ()
  runForwards _ (S1D y)  = ((), fromJust' . fromStorable . extract $ y)
  runForwards _ (S1DV y) = ((), S2DV y) -- We use no conversion at all. This should not be a problem as long as the runBackwards does it in the same way.
  runBackwards _ _ (S2D y)    = ((), fromJust' . fromStorable . flatten . extract $ y)
  runBackwards _ _ x@(S2DV y) = ((), S1DV y)

instance (KnownNat a, KnownNat x, KnownNat y, KnownNat (x * z), KnownNat z, a ~ (x * y * z)) => Layer Reshape ('D1 a) ('D3 x y z) where
  type Tape Reshape ('D1 a) ('D3 x y z) = ()
  runForwards _ (S1D y)  = ((), fromJust' . fromStorable . extract $ y)
  runForwards _ (S1DV _) = error "not yet implemented"
  runBackwards _ _ (S3D y)  = ((), fromJust' . fromStorable . flatten . extract $ y)

instance Serialize Reshape where
  put _ = return ()
  get = return Reshape

fromJust' :: Maybe x -> x
fromJust' (Just x) = x
fromJust' Nothing  = error "Reshape error: data shape couldn't be converted."

-------------------- DynamicNetwork instance --------------------

instance FromDynamicLayer Reshape where
  fromDynamicLayer inShape outShape _ =
     SpecNetLayer $ SpecReshape (tripleFromSomeShape inShape) (tripleFromSomeShape outShape)

instance ToDynamicLayer SpecReshape where
  toDynamicLayer _ _ (SpecReshape inp@(rowsI, colsI, depthI) out@(rowsO, colsO, depthO)) =
    reifyNat rowsI $ \(pxRI :: (KnownNat rI) => Proxy rI) ->
    reifyNat colsI $ \(_ :: (KnownNat cI) => Proxy cI) ->
    reifyNat depthI $ \(pxDI :: (KnownNat dI) => Proxy dI) ->
    reifyNat rowsO $ \(pxRO:: (KnownNat rO) => Proxy rO) ->
    reifyNat colsO $ \(_ :: (KnownNat cO) => Proxy cO) ->
    reifyNat depthO $ \(pxDO :: (KnownNat dO) => Proxy dO) ->
    case (singByProxy pxRI %* singByProxy pxDI, singByProxy pxRO %* singByProxy pxDO) of
          (SNat, SNat) ->
            if (rowsI, colsI, depthI) == (rowsO, colsO, depthO)
            then case (rowsI, colsI, depthI) of
              (_,1,1) -> return $ SpecLayer Trivial (sing :: Sing ('D1 rI)) (sing :: Sing ('D1 rI))
              (_,_,1) -> return $ SpecLayer Trivial (sing :: Sing ('D2 rI cI)) (sing :: Sing ('D2 rI cI))
              (_,_,_) -> return $ SpecLayer Trivial (sing :: Sing ('D3 rI cI dI)) (sing :: Sing ('D3 rI cI dI))
            else
            case ((rowsI, colsI, depthI), (rowsO, colsO, depthO)) of
              ((_,1,1), (_, _, 1)) | rowsI /= rowsO * colsO -> err
              ((_,1,1), (_, _, 1)) ->
                case (unsafeCoerce (Dict::Dict()) :: Dict (rI ~ (rO * cO))) of
                  Dict -> return $ SpecLayer Reshape (sing :: Sing ('D1 rI)) (sing :: Sing ('D2 rO cO))
              ((_,1,1), (_, _, _)) | rowsI /= rowsO * colsO * depthO -> err
              ((_,1,1), (_, _, _)) ->
                case (unsafeCoerce (Dict::Dict()) :: Dict (rI ~ ((rO * cO) * dO))) of
                  Dict -> return $ SpecLayer Reshape (sing :: Sing ('D1 rI)) (sing :: Sing ('D3 rO cO dO))
              ((_,_,1), (_, 1, 1)) | rowsI * colsI /= rowsO -> err
              ((_,_,1), (_, 1, 1)) ->
                case (unsafeCoerce (Dict::Dict()) :: Dict ((rI * cI) ~ rO)) of
                  Dict -> return $ SpecLayer Reshape (sing :: Sing ('D2 rI cI)) (sing :: Sing ('D1 rO))
              ((_,_,1), (_, _, _)) | 1 /= depthO -> err
              ((_,_,1), (_, _, _)) ->
                case (unsafeCoerce (Dict::Dict()) :: Dict (dO ~ 1)) of
                  Dict -> return $ SpecLayer Reshape (sing :: Sing ('D2 rI cI)) (sing :: Sing ('D3 rI cI dO))
              ((_,_,_), (_, 1, 1)) | rowsI * colsI * depthO /= rowsO -> err
              ((_,_,_), (_, 1, 1)) ->
                case (unsafeCoerce (Dict::Dict()) :: Dict (((rI * cI) * dI) ~ rO)) of
                  Dict -> return $ SpecLayer Reshape (sing :: Sing ('D3 rI cI dI)) (sing :: Sing ('D1 rO))
              ((_,_,_), (_, _, 1)) | depthI == 1 -> err
              ((_,_,_), (_, _, 1)) ->
                case (unsafeCoerce (Dict::Dict()) :: Dict (dI ~ 1)) of
                  Dict -> return $ SpecLayer Reshape (sing :: Sing ('D3 rI cI dI)) (sing :: Sing ('D2 rI cI))
              (_, _) -> error $ "Reshaping using a specificaiton from " ++ show inp ++ " to " ++ show out ++ " is not possible!"
    where err = error $ "cannot reshape from " ++ show inp ++ " to " ++ show out ++ ". Sizes (number of elements) do not match or it is trivial."

specReshape :: (Integer, Integer, Integer) -> (Integer, Integer, Integer) -> SpecNet
specReshape inp out = SpecNetLayer $ SpecReshape inp out

specReshape3D2D :: (Integer, Integer, Integer) -> (Integer, Integer) -> SpecNet
specReshape3D2D inp (rows, cols) = specReshape inp (rows, cols, 1)

specReshape3D1D :: (Integer, Integer, Integer) -> Integer -> SpecNet
specReshape3D1D inp rows = specReshape inp (rows, 1, 1)

specReshape2D3D :: (Integer, Integer) -> (Integer, Integer, Integer) -> SpecNet
specReshape2D3D (rows, cols) = specReshape (rows, cols, 1)

specReshape2D1D :: (Integer, Integer) -> Integer -> SpecNet
specReshape2D1D (rowsI, colsI) rowsO = specReshape (rowsI, colsI, 1) (rowsO, 1, 1)

specReshape1D3D :: Integer -> (Integer, Integer, Integer) -> SpecNet
specReshape1D3D rows = specReshape (rows, 1, 1)

specReshape1D2D :: Integer -> (Integer, Integer) -> SpecNet
specReshape1D2D rowsI (rowsO, colsO) = specReshape (rowsI, 1, 1) (rowsO, colsO, 1)

-- | Reshape to the given dimensions. Input and output nodes must match. If the input and output dimensions are the same, then this layer is omitted.
reshape :: Dimensions -> BuildM ()
reshape out@(rIn, cIn, dIn) = do
  inp@(rOut, cOut, dOut) <- buildGetLastLayerOut
  when (rIn * cIn * dIn /= rOut * cOut * dOut) $ error $ "Number of input and output nodes do not match in reshape! From " ++ show inp ++ " to " ++ show out
  if rIn == rOut && cIn == cOut && dIn == dOut
    then return () -- ignore as same size
    else buildAddSpec (specReshape inp out) >> buildSetLastLayer out

-------------------- GNum instances --------------------

instance GNum Reshape where
  _ |* Reshape = Reshape
  _ |+ Reshape = Reshape

max_stars_repo_name: schnecki/grenade

max_issues_repo_name: schnecki/grenade

max_forks_repo_name: schnecki/grenade


text:
module Main where

import CV.ColourUtils
import CV.DFT
import CV.Drawing
import CV.Edges
import CV.Filters
import CV.Image
import CV.ImageOp
import CV.Pixelwise
import Control.Applicative hiding ((<**>))
import qualified CV.ImageMath as IM
import qualified CV.Transforms as T
import Data.Complex

main = do
    Just x <- loadImage "smallLena.jpg"
    let  y, u, z :: Image GrayScale D32
         y = T.flip T.Horizontal x
         u = T.flip T.Vertical   y
         z = T.flip T.Vertical   x
         d :: Pixelwise (Complex D32)
         d = fromImage $ dft x
    saveImage "PixelWise.png" $ montage (2,5) 5
        [x,y
        ,stretchHistogram . toImage $ (fromImage z + fromImage y)
        ,stretchHistogram . toImage $ fromImage y + fromImage u + fromImage z
        ,stretchHistogram . toImage $ fmap (sin . (*9)) $ fromImage y
        ,stretchHistogram . toImage $ fmap log $ fromImage x

        ,stretchHistogram . gaussian (3,3)
                          . toImage $ atan2 <$$> (sobel (1,0) s5 x)
                                            <+>  (sobel (0,1) s5 x)
        ,toImage $ fmap (\x -> if x > 0.5 then 0 else 1) $ fromImage x
        ,toImage $ fmap (\x -> if x > 0.5 && x < 0.6 then 0 else 1)  $ fromImage x
        ,logarithmicCompression $ fst $ dftSplit $ (toImage d) <# circleOp (0:+0) (60,60) 30 Filled -- * (fromFunction (getSize d) (\(x,y) -> (exp(-(fromIntegral (91-x)^2+fromIntegral (91-y)^2)/(1600)):+0)))
        ]

max_stars_repo_name: maaleske/CV

max_issues_repo_name: maaleske/CV

max_forks_repo_name: maaleske/CV


text:
{-# LANGUAGE ForeignFunctionInterface #-}
{-|
Module      : Grenade.Layers.Internal.Add
Description : Fast addition functions that call efficient C function
Maintainer  : Theo Charalambous
License     : BSD2
Stability   : experimental
-}
module Grenade.Layers.Internal.BatchNorm (
    batchnorm
  ) where

import qualified Data.Vector.Storable        as U (unsafeFromForeignPtr0,
                                                   unsafeToForeignPtr0)

import           Foreign                     (mallocForeignPtrArray, withForeignPtr)
import           Foreign.Ptr                 (Ptr)
import           Numeric.LinearAlgebra       (Vector, Matrix, flatten)
import qualified Numeric.LinearAlgebra.Devel as U
import           System.IO.Unsafe            (unsafePerformIO)

import           Grenade.Types

-- | Efficient implementation of the forward propogation of batch normalization in testing mode.
batchnorm :: Int -> Int -> Int -> RealNum -> Matrix RealNum
             -> Vector RealNum -> Vector RealNum -> Vector RealNum -> Vector RealNum
             -> Matrix RealNum
batchnorm channels rows cols epsilon m gamma beta running_mean running_var
  = let outMatSize      = rows * cols * channels
        vec             = flatten m
    in unsafePerformIO $ do
      outPtr <- mallocForeignPtrArray outMatSize
      let (inPtr, _)    = U.unsafeToForeignPtr0 vec
          (gammaPtr, _) = U.unsafeToForeignPtr0 gamma
          (betaPtr, _)  = U.unsafeToForeignPtr0 beta
          (meanPtr, _)  = U.unsafeToForeignPtr0 running_mean
          (varPtr, _)   = U.unsafeToForeignPtr0 running_var

      withForeignPtr inPtr $ \inPtr' ->
        withForeignPtr gammaPtr $ \gammaPtr' ->
          withForeignPtr betaPtr $ \betaPtr' ->
            withForeignPtr meanPtr $ \meanPtr' ->
              withForeignPtr varPtr $ \varPtr' ->
                withForeignPtr outPtr $ \outPtr' ->
                  batchnorm_forward_cpu inPtr' gammaPtr' betaPtr' meanPtr' varPtr' epsilon channels rows cols outPtr'

      let matVec = U.unsafeFromForeignPtr0 outPtr outMatSize
      return (U.matrixFromVector U.RowMajor (rows * channels) cols matVec)
{-# INLINE batchnorm #-}

foreign import ccall unsafe
    batchnorm_forward_cpu
      :: Ptr RealNum -> Ptr RealNum -> Ptr RealNum -> Ptr RealNum -> Ptr RealNum -> RealNum -> Int -> Int -> Int -> Ptr RealNum -> IO ()

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
-- |
-- Module    : Statistics.Regression
-- Copyright : 2014 Bryan O'Sullivan
-- License   : BSD3
--
-- Functions for regression analysis.

module Statistics.Regression
    (
      olsRegress
    , bootstrapRegress
    ) where

import Control.Applicative
import Control.Concurrent (forkIO)
import Control.Concurrent.Chan (newChan, readChan, writeChan)
import Control.DeepSeq (rnf)
import Control.Monad (forM_, replicateM)
import GHC.Conc (getNumCapabilities)
import Statistics.Function as F
import Statistics.Matrix
import Statistics.Matrix.Algorithms (qr)
import Statistics.Types      (Estimate(..),ConfInt,CL,estimateFromInterval,significanceLevel)
import Statistics.Sample (mean)
import Statistics.Sample.Internal (sum)
import System.Random.MWC (GenIO, uniformR, splitGen)
import qualified Data.Vector as V
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Unboxed.Mutable as M
import Prelude hiding (pred, sum)

-- | Perform an ordinary least-squares regression on a set of
-- predictors, and calculate the goodness-of-fit of the regression.
--
-- The returned pair consists of:
--
-- * A vector of regression coefficients.  This vector has /one more/
--   element than the list of predictors; the last element is the
--   /y/-intercept value.
--
-- * /R²/, the coefficient of determination (see 'rSquare' for
--   details).
olsRegress :: [Vector]
              -- ^ Non-empty list of predictor vectors.  Must all have
              -- the same length.  These will become the columns of
              -- the matrix /A/ solved by 'ols'.
           -> Vector
              -- ^ Responder vector.  Must have the same length as the
              -- predictor vectors.
           -> (Vector, Double)
olsRegress preds@(_:_) resps
  | any (/=n) ls        = error $ "predictor vector length mismatch " ++
                                  show lss
  | G.length resps /= n = error $ "responder/predictor length mismatch " ++
                                  show (G.length resps, n)
  | otherwise           = (coeffs, rSquare mxpreds resps coeffs)
  where
    coeffs    = ols mxpreds resps
    mxpreds   = transpose .
                fromVector (length lss + 1) n .
                G.concat $ preds ++ [G.replicate n 1]
    lss@(n:ls) = map G.length preds
olsRegress _ _ = error "no predictors given"

-- | Compute the ordinary least-squares solution to /A x = b/.
ols :: Matrix     -- ^ /A/ has at least as many rows as columns.
    -> Vector     -- ^ /b/ has the same length as columns in /A/.
    -> Vector
ols a b
  | rs < cs   = error $ "fewer rows than columns " ++ show d
  | otherwise = solve r (transpose q `multiplyV` b)
  where
    d@(rs,cs) = dimension a
    (q,r)     = qr a

-- | Solve the equation /R x = b/.
solve :: Matrix     -- ^ /R/ is an upper-triangular square matrix.
      -> Vector     -- ^ /b/ is of the same length as rows\/columns in /R/.
      -> Vector
solve r b
  | n /= l    = error $ "row/vector mismatch " ++ show (n,l)
  | otherwise = U.create $ do
  s <- U.thaw b
  rfor n 0 $ \i -> do
    si <- (/ unsafeIndex r i i) <$> M.unsafeRead s i
    M.unsafeWrite s i si
    for 0 i $ \j -> F.unsafeModify s j $ subtract (unsafeIndex r j i * si)
  return s
  where n = rows r
        l = U.length b

-- | Compute /R²/, the coefficient of determination that
-- indicates goodness-of-fit of a regression.
--
-- This value will be 1 if the predictors fit perfectly, dropping to 0
-- if they have no explanatory power.
rSquare :: Matrix               -- ^ Predictors (regressors).
        -> Vector               -- ^ Responders.
        -> Vector               -- ^ Regression coefficients.
        -> Double
rSquare pred resp coeff = 1 - r / t
  where
    r   = sum $ flip U.imap resp $ \i x -> square (x - p i)
    t   = sum $ flip U.map resp $ \x -> square (x - mean resp)
    p i = sum . flip U.imap coeff $ \j -> (* unsafeIndex pred i j)

-- | Bootstrap a regression function.  Returns both the results of the
-- regression and the requested confidence interval values.
bootstrapRegress
  :: GenIO
  -> Int         -- ^ Number of resamples to compute.
  -> CL Double   -- ^ Confidence level.
  -> ([Vector] -> Vector -> (Vector, Double))
     -- ^ Regression function.
  -> [Vector]    -- ^ Predictor vectors.
  -> Vector      -- ^ Responder vector.
  -> IO (V.Vector (Estimate ConfInt Double), Estimate ConfInt Double)
bootstrapRegress gen0 numResamples cl rgrss preds0 resp0
  | numResamples < 1   = error $ "bootstrapRegress: number of resamples " ++
                                 "must be positive"
  | otherwise = do
  caps <- getNumCapabilities
  gens <- splitGen caps gen0
  done <- newChan
  forM_ (zip gens (balance caps numResamples)) $ \(gen,count) -> forkIO $ do
      v <- V.replicateM count $ do
           let n = U.length resp0
           ixs <- U.replicateM n $ uniformR (0,n-1) gen
           let resp  = U.backpermute resp0 ixs
               preds = map (flip U.backpermute ixs) preds0
           return $ rgrss preds resp
      rnf v `seq` writeChan done v
  (coeffsv, r2v) <- (G.unzip . V.concat) <$> replicateM caps (readChan done)
  let coeffs  = flip G.imap (G.convert coeffss) $ \i x ->
                est x . U.generate numResamples $ \k -> (coeffsv G.! k) G.! i
      r2      = est r2s (G.convert r2v)
      (coeffss, r2s) = rgrss preds0 resp0
      est s v = estimateFromInterval s (w G.! lo, w G.! hi) cl
        where w  = F.sort v
              lo = round c
              hi = truncate (n - c)
              n  = fromIntegral numResamples
              c  = n * (significanceLevel cl / 2)
  return (coeffs, r2)

-- | Balance units of work across workers.
balance :: Int -> Int -> [Int]
balance numSlices numItems = zipWith (+) (replicate numSlices q)
                                         (replicate r 1 ++ repeat 0)
 where (q,r) = numItems `quotRem` numSlices

max_stars_repo_name: runeksvendsen/hs-gauge

max_issues_repo_name: runeksvendsen/hs-gauge

max_forks_repo_name: runeksvendsen/hs-gauge


text:
-- {-# LANGUAGE ExistentialQuantification, RankNTypes #-}

module Language.Ruby.AST (Term(..)) where

import Data.Ratio (Rational)
import Data.Complex

data Term
       = Begin [Term]
       -- | RComplex (forall a. Num a => Complex (a))
       | Alias Term Term
       | And Term Term
       | Anddot
       | Arg String
       | Args [Term]
       | BackRef String
       | BlockArg String
       | BlockPass Term
       | Break [Term]
       | Case [Term]
       | Casgn Term String (Maybe Term) -- constants
       | Cbase
       | Const Term String
       | Cvar String
       | Cvasgn String (Maybe Term) -- class variable
       | Def String Term Term
       | Defined [Term]
       | Defs Term String Term Term
       | Dot
       | Dstr [Term]
       | Dsym [Term]
       | Encoding
       | Ensure Term Term
       | ERange Term Term
       | File
       | For Term Term Term
       | Gvar String
       | Gvasgn String (Maybe Term) -- global variables
       | Hash [Term]
       | Index Term [Term]
       | IndexAsgn Term [Term]
       | If Term Term Term
       | IRange Term Term
       | Ivar String
       | Ivasgn String (Maybe Term) -- instance variables
       | KWArg String
       | KWBegin [Term]
       | KWOptArg String Term
       | KWRestArg (Maybe String)
       | KWSplat Term
       | Lambda
       | Line
       | Lvar String
       | Lvasgn String (Maybe Term) -- variables
       | Masgn Term Term
       | Mlhs [Term]
       | Next [Term]
       | Nil
       | NthRef Int
       | Pair Term Term
       | Postexe Term
       | Preexe Term
       | Or Term Term
       | OptArg String Term
       | RArray [Term]
       | RComplex (Complex Double)
       | Redo [Term]
       | Retry [Term]
       | Return [Term]
       | RFalse
       | RFloat Double
       | RInt Int
       | RRational Rational
       | RTrue
       | Self
       | Send Term String [Term]
       | Csend Term String [Term]
       | Resbody Term Term Term
       | Rescue Term [Term]
       | RestArg (Maybe String)
       | ShadowArg String
       | Splat (Maybe Term)
       | Str String
       | Super [Term]
       | Sym String
       | Undef [Term]
       | Until Term Term
       | UntilPost Term Term
       | When [Term]
       | While Term Term
       | WhilePost Term Term
       | Yield [Term]
       | Zsuper [Term]
       deriving (Eq, Show)

max_stars_repo_name: flbulgarelli/language-ruby

max_issues_repo_name: flbulgarelli/language-ruby

max_forks_repo_name: flbulgarelli/language-ruby


text:
module Statistics.QQplot where

import qualified Data.List as DL
import qualified Statistics.Distribution.Normal as SDnorm
import qualified Statistics.Distribution as SD
{-
Python scipy code:

import numpy as np
from scipy.stats import beta
from scipy import stats

a = range(1,5)
b = 4-a+1
>>> pmean = beta.mean(a,b)
array([0.2, 0.4, 0.6, 0.8])
>>> pmedia = beta.median(a,b)
array([0.15910358, 0.38572757, 0.61427243, 0.84089642])
>>> dist = stats.norm

>>> pmeanppf = dist.ppf(pmean)
>>> pmeanppf
array([-0.84162123, -0.2533471 ,  0.2533471 ,  0.84162123])

>>> pmedia = beta.median(a,b)
>>> pmediappf = dist.ppf(pmedia)
>>> pmediappf
array([-0.99814888, -0.29047204,  0.29047204,  0.99814888])
-}

{-
Haskell code:
import Playground.Stats

*Main Playground.Stats> let hmean = betaMean 4
*Main Playground.Stats> hmean
[0.2,0.4,0.6,0.8]

*Main Playground.Stats> let hmedian = betaMedian 4
*Main Playground.Stats> hmedian
[0.1591035847462855,0.3854524627720504,0.6145475372279496,0.8408964152537145]
Main Playground.Stats> zipWith (-) hmedian pmedian
[4.74628550350964e-9,-2.751072279496003e-4,2.751072279496558e-4,-4.746285475754064e-9]

*Main Playground.Stats> hmeanppf = normPPF hmean
*Main Playground.Stats> hmeanppf
[-0.8416212335729141,-0.2533471031357998,0.2533471031357998,0.8416212335729144]

*Main Playground.Stats> hmedianppf = normPPF hmedian
*Main Playground.Stats> hmedianppf
[-0.9981488825015569,-0.2911914182975527,0.29119141829755285,0.9981488825015569]

*Main Playground.Stats> zipWith (-) hmedianppf pmedianppf
[-2.5015568505537544e-9,-7.193782975526974e-4,7.19378297552864e-4,2.5015568505537544e-9]

*Main Playground.Stats> zipWith (-) hmeanppf pmeanppf
[-3.572914075711253e-9,-3.135799786502247e-9,3.135799786502247e-9,3.5729144087781606e-9]
-}
-- | About beta mean & beta median
-- | scipy.stats.morestats
-- | https://github.com/scipy/scipy/blob/abdab61d65dda1591f9d742230f0d1459fd7c0fa/scipy/stats/morestats.py#L523
-- | line: 393
-- | betaMean definition is the same as
-- | https://github.com/statsmodels/statsmodels/blob/66fc298c51dc323ce8ab8564b07b1b3797108dad/statsmodels/graphics/gofplots.py#L58
-- | line: 639, function: plotting_pos
-- | betaMedain approximation definition is
-- | https://github.com/scipy/scipy/blob/abdab61d65dda1591f9d742230f0d1459fd7c0fa/scipy/stats/morestats.py#L523
-- | line 361, function: _calc_uniform_order_statistic_medians
-- | This approximation also being introduced in doc scipy.stats.probplot:
-- | https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.stats.probplot.html

-- | statistic has function to compute mean for beta distribution but no median !!!
betaMean :: (Enum a, Fractional a) => a -> [a]
betaMean nobs = (\x-> x / (nobs+1)) <$> [1.0..nobs]

-- | Filliben’s estimate
betaMedian :: (Enum a, Floating a) => a -> [a]
betaMedian nobs = headElem : body ++ [lastElem]
                    where
                        headElem = 1 - 0.5**(1/nobs)
                        lastElem = 0.5**(1/nobs)
                        body = (\x-> (x-0.3175)/(nobs+0.365)) <$> [2.0..nobs-1]

-- | http://hackage.haskell.org/package/statistics-0.15.1.1/docs/src/Statistics.Distribution.Normal.html
-- | inverse CDF function of normal distribution
normPPF quantile =
    let norm = SDnorm.normalDistr 0 1
    in realToFrac . SD.quantile norm <$> quantile

sort :: Ord a => [a] -> [a]
sort = DL.sort

max_stars_repo_name: emmettng/hExample

max_issues_repo_name: emmettng/hExample

max_forks_repo_name: emmettng/hExample