text:
{-# LANGUAGE GADTs             #-}
{-# LANGUAGE OverloadedStrings #-}

module MatrixMarket (Matrix(..), readMatrix) where

import Control.Applicative                      hiding ( many )

import Data.Int
import Data.Complex
import Data.Attoparsec.ByteString.Char8
import Data.ByteString.Lex.Fractional
import qualified Data.Attoparsec.Lazy           as L
import qualified Data.ByteString.Lazy           as L

-- | Specifies the element type.  Pattern matrices do not have any elements,
-- only indices, and only make sense for coordinate matrices and vectors.
--
data Field  = Real | Complex | Integer | Pattern
    deriving (Eq, Show)

-- | Specifies either sparse or dense storage.  In sparse (\"coordinate\")
-- storage, elements are given in (i,j,x) triplets for matrices (or (i,x) for
-- vectors).  Indices are 1-based, so that A(1,1) is the first element of a
-- matrix, and x(1) is the first element of a vector.
--
-- In dense (\"array\") storage, elements are given in column-major order.
--
-- In both cases, each element is given on a separate line.
--
data Format = Coordinate | Array
    deriving (Eq, Show)

-- | Specifies any special structure in the matrix.  For symmetric and hermition
-- matrices, only the lower-triangular part of the matrix is given. For skew
-- matrices, only the entries below the diagonal are stored.
--
data Structure = General | Symmetric | Hermitian | Skew
    deriving (Eq, Show)

-- We really want a type parameter to Matrix, but I think that requires some
-- kind of dynamic typing so that we can determine (a ~ Integral) or (a ~
-- RealFloat), and so forth, depending on the file being read. This will do for
-- our purposes...
--
-- Format is: (rows,columns) nnz [(row,column,value)]
--
data Matrix
  = PatternMatrix (Int,Int) Int [(Int32,Int32)]
  | IntMatrix     (Int,Int) Int [(Int32,Int32,Int)]
  | RealMatrix    (Int,Int) Int [(Int32,Int32,Float)]
  | ComplexMatrix (Int,Int) Int [(Int32,Int32,Complex Float)]
  deriving Show

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

comment :: Parser ()
comment = char '%' *> skipWhile (not . eol) *> endOfLine
  where
    eol w = w `elem` ("\n\r" :: String)

floating :: Fractional a => Parser a
floating = do
  mv <- readDecimal <$> (skipSpace *> takeTill isSpace)  -- readDecimal does the fancy stuff
  case mv of
       Just (v,_) -> return v
       Nothing    -> fail "floating-point number"

integral :: Integral a => Parser a
integral = skipSpace *> decimal

format :: Parser Format
format =  string "coordinate" *> pure Coordinate
      <|> string "array"      *> pure Array
      <?> "matrix format"

field :: Parser Field
field =  string "real"    *> pure Real
     <|> string "complex" *> pure Complex
     <|> string "integer" *> pure Integer
     <|> string "pattern" *> pure Pattern
     <?> "matrix field"

structure :: Parser Structure
structure =  string "general"        *> pure General
         <|> string "symmetric"      *> pure Symmetric
         <|> string "hermitian"      *> pure Hermitian
         <|> string "skew-symmetric" *> pure Skew
         <?> "matrix structure"

header :: Parser (Format,Field,Structure)
header =  string "%%MatrixMarket matrix"
       >> (,,) <$> (skipSpace *> format)
               <*> (skipSpace *> field)
               <*> (skipSpace *> structure)
               <*  endOfLine
               <?> "MatrixMarket header"

extent :: Parser (Int,Int,Int)
extent = do
  [m,n,l] <- skipWhile isSpace *> count 3 integral <* endOfLine
  return (m,n,l)

line :: Integral i => Parser a -> Parser (i,i,a)
line f = (,,) <$> integral
              <*> integral
              <*> f
              <*  endOfLine

--------------------------------------------------------------------------------
-- Matrix Market
--------------------------------------------------------------------------------

matrix :: Parser Matrix
matrix = do
  (_,t,_) <- header
  (m,n,l) <- skipMany comment *> extent
  case t of
    Real    -> RealMatrix    (m,n) l `fmap` many1 (line floating)
    Complex -> ComplexMatrix (m,n) l `fmap` many1 (line ((:+) <$> floating <*> floating))
    Integer -> IntMatrix     (m,n) l `fmap` many1 (line integral)
    Pattern -> PatternMatrix (m,n) l `fmap` many1 ((,) <$> integral <*> integral)

readMatrix :: FilePath -> IO Matrix
readMatrix file = do
  chunks <- L.readFile file
  case L.parse matrix chunks of
    L.Fail _ _ msg      -> error $ file ++ ": " ++ msg
    L.Done _ mtx        -> return mtx

max_stars_repo_name: jmatsushita/accelerate-examples

max_issues_repo_name: gckeller/accelerate-examples

max_forks_repo_name: gckeller/accelerate-examples


text:
{-# OPTIONS_GHC -fno-warn-type-defaults #-}
{-# LANGUAGE DataKinds, FlexibleContexts, FlexibleInstances, GADTs #-}
{-# LANGUAGE MultiParamTypeClasses, NoMonomorphismRestriction      #-}
{-# LANGUAGE OverloadedStrings, PolyKinds, TypeFamilies            #-}
{-# LANGUAGE UndecidableInstances                                  #-}
{-# OPTIONS_GHC -fno-warn-type-defaults -fno-warn-orphans -Wall #-}
module Main where
import           Algebra.Bridge.Singular
import           Algebra.Prelude.Core
import           Control.Monad                   (void)
import           Control.Parallel.Strategies
import qualified Data.Text                       as T
import qualified Data.Vector.Unboxed             as V
import           Numeric.Field.Fraction          (Fraction)
import           Prelude                         (read)
import           Statistics.Resampling
import           Statistics.Resampling.Bootstrap
import           Statistics.Types
import qualified System.Random.MWC               as Rand

i2 :: [OrderedPolynomial (Fraction Integer) Grevlex 5]
i2 =  [35 * y^4 - 30*x*y^2 - 210*y^2*z + 3*x^2 + 30*x*z - 105*z^2 +140*y*t - 21*u
      ,5*x*y^3 - 140*y^3*z - 3*x^2*y + 45*x*y*z - 420*y*z^2 + 210*y^2*t -25*x*t + 70*z*t + 126*y*u
      ]
     where [t,u,x,y,z] = vars

i3 :: [OrderedPolynomial (Fraction Integer) Grevlex 4]
i3 = [ x^31 - x^6 - x- y, x^8 - z, x^10 -t]
  where
    [t,x,y,z] = vars

i4 :: [OrderedPolynomial (Fraction Integer) Grevlex 4]
i4 = [ w+x+y+z, w*x+x*y+y*z+z*w, w*x*y + x*y*z + y*z*w + z*w*x, w*x*y*z - 1]
  where
    [x,y,z,w] = vars

benchIdeal :: IsSingularPolynomial poly => Text -> Ideal poly -> IO Double
benchIdeal fun i = fmap ((/1000) . read . T.unpack) $ singular $
  toProg fun i

toProg :: (SingularOrder (Arity poly) (MOrder poly), SingularCoeff (Coefficient poly), IsOrderedPolynomial poly) => Text -> Ideal poly -> Text
toProg fun i =
  let expr = funE fun [idealE' i]
  in prettySingular $ do
    directC "system(\"--ticks-per-sec\", 1000000)"
    void $ ringC "R" expr
    declOnlyC IdealT "G"
    directC "timer=1"
    directC "int t=rtimer"
    letC "G" expr
    directC "print(rtimer-t)"
    directC "exit"

analyse :: (IsSingularPolynomial poly) => String -> Text -> Ideal poly -> IO ()
analyse lab fun ideal = do
  gen <- Rand.create
  i2Gr <- V.replicateM 200 $ benchIdeal fun ideal
  res  <- resample gen [Mean, StdDev] 1000 i2Gr
  let [Estimate mn (ConfInt lbmn ubmn _) ,Estimate dv _]
        = bootstrapBCA cl95 i2Gr res
  putStrLn lab
  putStrLn $ unlines $ map ('\t':)
    ["Mean:\t" ++ show mn ++ "(ms)"
    ,"MeanLB:\t" ++ show lbmn
    ,"MeanUB:\t" ++ show ubmn
    ,"StdDev:\t" ++ show dv
    ]

runTestCases :: (IsMonomialOrder n o, SingularCoeff k, KnownNat n)
             => String -> Ideal (OrderedPolynomial k o n) -> IO ()
runTestCases lab i = do
  analyse (lab ++ " (Lex, Sing(gr))") "groebner" $ fmap (changeOrder Lex) i
  analyse (lab ++ " (Lex, Sing(sba))") "sba" $ fmap (changeOrder Lex) i
  analyse (lab ++ " (Grevlex, Sing(gr))") "groebner" $ fmap (changeOrder Grevlex) i
  analyse (lab ++ " (Grevlex, Sing(sba))") "sba" $ fmap (changeOrder Grevlex) i

main :: IO ()
main = do
  ideal2 <- return $! (toIdeal i2 `using` rdeepseq)
  ideal3 <- return $! (toIdeal i3 `using` rdeepseq)
  ideal4 <- return $! (toIdeal i4 `using` rdeepseq)

  runTestCases "I1" ideal2
  runTestCases "Cyclic4" ideal4
  analyse "I3 (Grevlex, Sing(gr))" "groebner" ideal3
  analyse "I3 (Grevlex, Sing(sba))" "sba" ideal3
  -- analyse "groebner" ideal4
  -- analyse "sba" ideal4
  -- analyse "groebner" ideal3
  -- analyse "sba" ideal3
  -- defaultMainWith defaultConfig $
  --      mkTestCases (Left 1) ideal2
  --   ++ mkTestCases (Right "Cyclic-4") ideal4
  --   ++ [mkTC "grevlex03" ideal3]

ext: hs

max_stars_repo_name: hangingman/computational-algebra

max_issues_repo_name: hangingman/computational-algebra

max_forks_repo_name: hangingman/computational-algebra


text:
-- code from Edward Kmett's blog post
-- Mirrored Lenses
-- http://comonad.com/reader/2012/mirrored-lenses/

-- Slightly adapted since Data.PolyLens doesn't use Rank 2 types.
-- Mostly this means adding "Functor f =>" everywhere

import Data.Complex
import Data.PolyLens

realLens :: (RealFloat a, Functor f) => Lens (Complex a) a f
realLens f (r :+ i) = fmap (:+ i) (f r)

imagLens :: (RealFloat a, Functor f) => Lens (Complex a) a f
imagLens f (r :+ i) = fmap (r :+) (f i)

fstLens :: Functor f => PolyLens (a,c) (b,c) a b f
fstLens f (a, b) = fmap (\x -> (x, b)) (f a)

sndLens :: Functor f => PolyLens (a,b) (a,c) b c f
sndLens f (a, b) = fmap ((,) a) (f b)

swap :: (a, b) -> (b, a)
swap (a, b) = (b, a)

swapped :: Functor f => PolyLens (a, b) (c, d) (b, a) (d, c) f
swapped = iso swap swap

negated :: (Num a, Functor f) => Lens a a f
negated = iso negate negate

traverseLens :: ((c -> c) -> a -> b) -> a -> b
traverseLens f = f id

getFst :: Getter (a, b) a r
getFst = getting fst

getSnd :: Getter (a, b) b r
getSnd = getting snd

getPhase :: RealFloat a => Getter (Complex a) a r
getPhase = getting phase

getAbs, getSignum :: Num a => Getter a a r
getAbs = getting abs
getSignum = getting signum

plus, times :: Num a => Setter a a a
plus = setting (+)
times = setting (*)

reals :: (RealFloat a, RealFloat b) => Modifier (Complex a) (Complex b) a b
reals = modifying (\f (r :+ i) -> f r :+ f i)

max_stars_repo_name: DanBurton/polylens

max_issues_repo_name: DanBurton/polylens

max_forks_repo_name: DanBurton/polylens


text:
module Deanie.Measure (
    measure

  -- * queries

  , integrate
  , expectation
  , variance
  , mgf
  , cgf
  , cdf
  ) where

import Control.Monad
import Data.List (foldl')
import Deanie.Language
import Control.Foldl (Fold)
import Numeric.Integration.TanhSinh
import Numeric.SpecFunctions

newtype Measure a = Measure ((a -> Double) -> Double)

integrate :: (a -> Double) -> Measure a -> Double
integrate f (Measure nu) = nu f

expectation :: Measure Double -> Double
expectation = integrate id

variance :: Measure Double -> Double
variance nu = integrate (^ 2) nu - expectation nu ^ 2

mgf :: Measure Double -> Double -> Double
mgf mu t = integrate (\x -> exp (t * x)) mu

cgf :: Measure Double -> Double -> Double
cgf mu = log . mgf mu

cdf :: Measure Double -> Double -> Double
cdf nu x = integrate (negativeInfinity `to` x) nu where
  negativeInfinity :: Double
  negativeInfinity = negate (1 / 0)

  to :: (Num a, Ord a) => a -> a -> a -> a
  to a b x
    | x >= a && x <= b = 1
    | otherwise        = 0

instance Functor Measure where
  fmap f nu = Measure $ \g ->
    integrate (g . f) nu

instance Applicative Measure where
  pure x = Measure (\f -> f x)
  Measure h <*> Measure g = Measure $ \f ->
    h (\k -> g (f . k))

instance Monad Measure where
  return x  = Measure (\f -> f x)
  rho >>= g = Measure $ \f ->
    integrate (\nu -> integrate f (g nu)) rho

fromMassFunction :: Foldable f => (a -> Double) -> f a -> Measure a
fromMassFunction f support = Measure $ \g ->
  foldl' (\acc x -> acc + f x * g x) 0 support

fromDensityFunction :: (Double -> Double) -> Measure Double
fromDensityFunction d = Measure $ \f ->
    quadratureTanhSinh (\x -> f x * d x)
  where
    quadratureTanhSinh = result . last . everywhere trap

mbernoulli :: Double -> Measure Bool
mbernoulli p = fromMassFunction (pmf p) [False, True] where
  pmf p x
    | p < 0 || p > 1 = 0
    | otherwise      = if x then p else 1 - p

mbeta :: Double -> Double -> Measure Double
mbeta a b = fromDensityFunction (density a b) where
  density a b p
    | p < 0 || p > 1 = 0
    | otherwise      = 1 / exp (logBeta a b) * p ** (a - 1) * (1 - p) ** (b - 1)

mgamma :: Double -> Double -> Measure Double
mgamma a b = fromDensityFunction (density a b) where
  density a b x
    | a < 0 || b < 0 = 0
    | otherwise  =
       b ** a / exp (logGamma a) * x ** (a - 1) * exp (negate (b * x))

mgaussian :: Double -> Double -> Measure Double
mgaussian m s = fromDensityFunction (density m s) where
  density m s x
    | s <= 0    = 0
    | otherwise =
        1 / (s * sqrt (2 * pi)) *
          exp (negate ((x - m) ^^ 2) / (2 * (s ^^ 2)))

measure :: Program a -> Measure a
measure = iterM $ \case
    ProgramF (InL term) -> evalAlg term
    ProgramF (InR term) -> join (runAp measure term)
  where
    evalAlg = \case
      BernoulliF p k  -> mbernoulli p >>= k
      BetaF a b k     -> mbeta a b >>= k
      GammaF a b k    -> mgamma a b >>= k
      GaussianF m s k -> mgaussian m s >>= k

max_stars_repo_name: jtobin/deanie

max_issues_repo_name: silky/deanie

max_forks_repo_name: silky/deanie


text:
{-# LANGUAGE LambdaCase, DoAndIfThenElse, GeneralizedNewtypeDeriving #-}
module Truthiness where

import Data.Text (Text)
import qualified Data.Text as Text
import Data.ByteString (ByteString)
import qualified Data.ByteString as ByteString

import Data.Map (Map)
import qualified Data.Map as Map
import Data.Set (Set)
import qualified Data.Set as Set
import Data.Monoid (Any(..), All(..))
import System.Exit (ExitCode(..))
import Numeric.Natural (Natural)
import Data.Complex (Complex(..))
import Data.Ratio (Ratio)
import Control.Applicative (Alternative(..), Alternative(..))
import Control.Monad (MonadPlus(..))
import Data.Foldable (Foldable(..))

{-| Pythonic truthiness/falsines. see <https://docs.python.org/2.4/lib/truth.html>.

keep around information (e.g. use @Maybe@, a custom flag type, etc.)
to avoid <https://existentialtype.wordpress.com/2011/03/15/boolean-blindness/ Boolean Blindness>,
then drop that information with 'Boolean', when necessary for convenience.

e.g. 'when'

when 'Foldable', instances should satisfy the following law:

@
'falsy' = 'null' . 'toList'
@

when a 'Monoid', instances should satisfy the following law:

@
'falsy' = ('mempty' ==)
@

when a 'Num', instances should satisfy the following law:

@
'falsy' = (0 ==)
@

non-instances:

* @()@
* @Dual@

-}
class Boolean a where

 falsy :: a -> Bool

 truthy :: a -> Bool
 truthy = not . falsy

-- booleans
instance Boolean Bool         where falsy = (== False)
instance Boolean Any          where falsy = falsy . getAny
instance Boolean All          where falsy = falsy . getAll

-- errors
instance Boolean (Maybe a)    where falsy = (\case Nothing -> True; _ -> False) -- not (== Nothing)
instance Boolean (Either e a) where falsy = (\case Left{} -> True; _ -> False)
-- | (includes 'String's)
instance Boolean [a]          where falsy = (\case [] -> True; _ -> False) -- not (== [])

-- numbers
instance Boolean Int          where falsy = (== 0)
instance Boolean Word         where falsy = (== 0)
instance Boolean Integer      where falsy = (== 0)
instance Boolean Natural      where falsy = (== 0)
instance (Integral a) => Boolean (Ratio a)    where falsy = (== 0)
instance (RealFloat a) => Boolean (Complex a)  where falsy = (== 0)

-- containers
instance Boolean (Set a)      where falsy = Set.null
instance Boolean (Map k v)    where falsy = Map.null

-- strings
instance Boolean Char         where falsy = (== '\NUL')
instance Boolean Text         where falsy = Text.null
instance Boolean ByteString   where falsy = ByteString.null

-- et cetera
instance Boolean ExitCode      where falsy = (\case ExitFailure{} -> True; _ -> False)

instance (Foldable t)       => Boolean (WrappedFoldable t a) where falsy = null . toList
instance (Monoid a, Eq a)   => Boolean (WrappedMonoid a)     where falsy = (mempty ==)
instance (Num a,    Eq a)   => Boolean (WrappedNum a)        where falsy = (0 ==)

newtype WrappedFoldable t a = WrappedFoldable { getWrappedFoldable :: t a } deriving (Foldable)
newtype WrappedMonoid   a   = WrappedMonoid   { getWrappedMonoid   :: a } deriving (Monoid,Eq)
newtype WrappedNum      a   = WrappedNum      { getWrappedNum      :: a } deriving (Num,Eq)

-- control

whenB :: (Applicative m, Boolean b) => b -> m () -> m ()
whenB condition action = ifB condition action nothing
{-# INLINEABLE whenB #-}

whenM :: (Monad m, Boolean b) => m b -> m () -> m ()
whenM condition action = ifM condition action nothing
{-# INLINEABLE whenM #-}

guardB :: (Alternative m, Boolean b) => b -> m ()
guardB condition = unlessB condition empty
{-# INLINEABLE guardB #-}

guardM :: (MonadPlus m, Boolean b) => m b -> m ()
guardM condition = unlessM condition mzero
{-# INLINEABLE guardM #-}

unlessB :: (Applicative m, Boolean b) => b -> m () -> m ()
unlessB condition action = ifB condition nothing action
{-# INLINEABLE unlessB #-}

unlessM :: (Monad m, Boolean b) => m b -> m () -> m ()
unlessM condition action = ifM condition nothing action
{-# INLINEABLE unlessM #-}

boolB :: (Boolean b) => a -> a -> b -> a
boolB x y c = ifB c x y
{-# INLINEABLE boolB #-}

ifB :: (Boolean b) => b -> a -> a -> a
ifB c x y =
 if   truthy c
 then x
 else y
{-# INLINEABLE ifB #-}

ifM :: (Monad m, Boolean b) => m b -> m a -> m a -> m a
ifM condition a b = (\c -> ifB c a b) =<< condition
{-# INLINEABLE ifM #-}

nothing :: (Applicative m) => m ()
nothing = pure()

{-$ alternatives

too explicit: <https://hackage.haskell.org/package/base-4.8.2.0/docs/Control-Monad.html>

fewer instances: <https://hackage.haskell.org/package/cond>

-}

max_stars_repo_name: sboosali/truthiness

max_issues_repo_name: sboosali/truthiness

max_forks_repo_name: sboosali/truthiness


text:
{-# LANGUAGE GADTs #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE DataKinds, PolyKinds #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE FlexibleContexts, FlexibleInstances #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE MultiParamTypeClasses #-}

module TensorHMatrix where

import Nats
import Numeric.LinearAlgebra hiding ((!))
import Data.Vector.Storable
import Data.Default
import Data.Singletons.Prelude

import Data.Vinyl

import Prelude hiding (zipWith)

type Usable a = (Element a, Num a, Numeric a, Num (Vector a), Container Vector a)

type IntegralK (p :: KProxy k) = (SingKind p, Integral (DemoteRep p))
type IntegralN (n :: k) = IntegralK ('KProxy :: KProxy k)

natVal' :: (Num a, IntegralN n) => Sing n -> a
natVal' = fromIntegral . fromSing

data Tensor (dims :: [k]) a where
  Scalar :: a -> Tensor '[] a
  Vector :: Vector a -> Tensor '[n] a
  Matrix :: Matrix a -> Tensor '[n, m] a

deriving instance (Show a, Element a) => Show (Tensor dims a)

instance (Default a) => Default (Tensor '[] a) where
  def = Scalar def

fill1 :: forall a n. (SingI n, IntegralN n, Usable a) => a -> Tensor '[n] a
fill1 a = Vector $ konst a (natVal' (sing::Sing n))

instance (SingI n, IntegralN n, Default a, Usable a) => Default (Tensor '[n] a) where
  def = fill1 def

instance (SingI n, IntegralN n, SingI m, Default a, Usable a) => Default (Tensor '[n, m] a) where
  def = Matrix $ konst def (natVal' (sing::Sing n), natVal' (sing::Sing m))

--type instance IndexOf (Tensor l) = Rec LT l
--type instance IndexOf (Tensor l) = Rec (Const Int) l
--type instance IndexOf (Tensor l) = Vec (Length l) Int

{- Fails backwards fundep :(
instance Indexable (Tensor (n ': l) a) (Tensor l a) where
  Matrix m ! i = Vector (m ! i)
  Vector v ! i = Scalar (v ! i)
-}

select :: Storable a => Int -> Tensor '[n] a -> a
select i (Vector v) = v ! i

-- pretty much all of these instances should be automatically derivable

instance Num a => Num (Tensor '[] a) where
  Scalar a + Scalar b = Scalar (a + b)
  Scalar a - Scalar b = Scalar (a - b)
  Scalar a * Scalar b = Scalar (a * b)
  negate (Scalar a) = Scalar (negate a)
  abs (Scalar a) = Scalar (abs a)
  signum (Scalar a) = Scalar (signum a)
  fromInteger n = Scalar (fromInteger n)

instance (SingI n, IntegralN n, Usable a) => Num (Tensor '[n] a) where
  Vector a + Vector b = Vector (a + b)
  Vector a - Vector b = Vector (a - b)
  Vector a * Vector b = Vector (a * b)
  negate (Vector a) = Vector (negate a)
  abs (Vector a) = Vector (abs a)
  signum (Vector a) = Vector (signum a)
  fromInteger n = Vector $ konst (fromInteger n) (natVal' (sing::Sing n))

instance (SingI n, IntegralN n, SingI m, Usable a) => Num (Tensor '[n, m] a) where
  Matrix a + Matrix b = Matrix (a + b)
  Matrix a - Matrix b = Matrix (a - b)
  Matrix a * Matrix b = Matrix (a * b)
  negate (Matrix a) = Matrix (negate a)
  abs (Matrix a) = Matrix (abs a)
  signum (Matrix a) = Matrix (signum a)
  fromInteger n = Matrix $ konst (fromInteger n) (natVal' (sing::Sing n), natVal' (sing::Sing m))

instance Fractional a => Fractional (Tensor '[] a) where
  Scalar a / Scalar b = Scalar (a / b)
  recip (Scalar a) = Scalar (recip a)
  fromRational r = Scalar (fromRational r)

instance (SingI n, IntegralN n, Numeric a, Fractional a, Fractional (Vector a)) => Fractional (Tensor '[n] a) where
  Vector a / Vector b = Vector (a / b)
  recip (Vector a) = Vector (recip a)
  fromRational r = Vector $ konst (fromRational r) (natVal' (sing::Sing n))

instance Floating a => Floating (Tensor '[] a) where
  pi = Scalar pi
  exp (Scalar a) = Scalar (exp a)
  log (Scalar a) = Scalar (log a)
  sqrt (Scalar a) = Scalar (sqrt a)
  Scalar a ** Scalar b = Scalar (a ** b)
  logBase (Scalar a) (Scalar b) = Scalar (logBase a b)
  sin (Scalar a) = Scalar (sin a)
  cos (Scalar a) = Scalar (cos a)
  tan (Scalar a) = Scalar (tan a)
  asin (Scalar a) = Scalar (asin a)
  acos (Scalar a) = Scalar (acos a)
  atan (Scalar a) = Scalar (atan a)
  sinh (Scalar a) = Scalar (sinh a)
  cosh (Scalar a) = Scalar (cosh a)
  tanh (Scalar a) = Scalar (tanh a)
  asinh (Scalar a) = Scalar (asinh a)
  acosh (Scalar a) = Scalar (acosh a)
  atanh (Scalar a) = Scalar (atanh a)

instance (SingI n, IntegralN n, Usable a, Floating a, Floating (Vector a)) => Floating (Tensor '[n] a) where
  pi = fill1 pi
  exp (Vector a) = Vector (exp a)
  log (Vector a) = Vector (log a)
  sqrt (Vector a) = Vector (sqrt a)
  Vector a ** Vector b = Vector (a ** b)
  logBase (Vector a) (Vector b) = Vector (logBase a b)
  sin (Vector a) = Vector (sin a)
  cos (Vector a) = Vector (cos a)
  tan (Vector a) = Vector (tan a)
  asin (Vector a) = Vector (asin a)
  acos (Vector a) = Vector (acos a)
  atan (Vector a) = Vector (atan a)
  sinh (Vector a) = Vector (sinh a)
  cosh (Vector a) = Vector (cosh a)
  tanh (Vector a) = Vector (tanh a)
  asinh (Vector a) = Vector (asinh a)
  acosh (Vector a) = Vector (acosh a)
  atanh (Vector a) = Vector (atanh a)

tmap :: (Container Vector a, Num a, Element b) => (a -> b) -> Tensor dims a -> Tensor dims b
tmap f (Scalar a) = Scalar $ f a
tmap f (Vector v) = Vector $ cmap f v
tmap f (Matrix m) = Matrix $ cmap f m

tZipWith :: (Storable a, Storable b, Storable c) => (a -> b -> c) -> Tensor '[n] a -> Tensor '[n] b -> Tensor '[n] c
tZipWith f (Vector a) (Vector b) = Vector (zipWith f a b)

transpose :: (Numeric a) => Tensor '[n, m] a -> Tensor '[m, n] a
transpose (Matrix m) = Matrix (tr' m)

asCol :: Storable a => Tensor '[n] a -> Tensor '[n, FromInteger 1] a
asCol (Vector v) = Matrix $ asColumn v

asRow' :: Storable a => Tensor '[n] a -> Tensor '[FromInteger 1, n] a
asRow' (Vector v) = Matrix $ asRow v

dot :: (Storable a, Numeric a) => Tensor '[n] a -> Tensor '[n] a -> a
dot (Vector a) (Vector b) = a <.> b

mv :: Numeric a => Tensor '[n, m] a -> Tensor '[m] a -> Tensor '[n] a
mv (Matrix m) (Vector v) = Vector $ m #> v

mm :: Numeric a => Tensor '[n, m] a -> Tensor '[m, k] a -> Tensor '[n, k] a
mm (Matrix m1) (Matrix m2) = Matrix $ m1 <> m2

gradDot :: Numeric a => Tensor '[n] a -> Tensor '[n] a -> a -> (Tensor '[n] a, Tensor '[n] a)
gradDot (Vector a) (Vector b) g = (Vector $ scale g b, Vector $ scale g a)

gradMV :: Numeric a => Tensor '[n, m] a -> Tensor '[m] a -> Tensor '[n] a -> (Tensor '[n, m] a, Tensor '[m] a)
gradMV m v g = (mm (asCol g) (asRow' v), mv (transpose m) g)

gradMM :: Numeric a => Tensor '[n, m] a -> Tensor '[m, k] a -> Tensor '[n, k] a -> (Tensor '[n, m] a, Tensor '[m, k] a)
gradMM a b g = (mm g (transpose b), mm (transpose a) g)

oneHot :: forall a n. (SingI n, IntegralN n, Usable a) => Int -> Tensor '[n] a
oneHot m = Vector $ (konst 0 (natVal' (sing::Sing n))) // [(m, 1)]

gradSelect :: forall a n. (SingI n, IntegralN n, Usable a) => Int -> Tensor '[n] a -> a -> Tensor '[n] a
gradSelect i _ a = Vector $ (konst 0 (natVal' (sing::Sing n))) // [(i, a)]


text:
{-# LANGUAGE OverloadedStrings, NoImplicitPrelude, BangPatterns, DeriveGeneric #-}
module NLP.Albemarle.Sparse
    ( dense
    , fromDocuments
    , transpose
    , shift
    , size
    , sparse
    , mult
    , multCol
    , SparseMatrix(..)
    , Coord
    , Coords
    ) where
import ClassyPrelude
import NLP.Albemarle hiding (SparseMatrix)
import GHC.Generics (Generic)
import Data.Vector.Binary
import qualified Data.Vector.Unboxed as Vec
import qualified Data.Vector.Storable as SVec
import qualified Data.Vector.Algorithms.Intro as Introsort
import qualified Numeric.LinearAlgebra as HMatrix
import qualified Data.Binary as Bin
import Control.Parallel.Strategies (rpar, rseq, parList, withStrategy, dot, rdeepseq)
type Coord = (Int, Int, Double)
type Coords = Vec.Vector Coord
data SparseMatrix = RowMatrix !Coords | ColMatrix !Coords
  deriving (Show, Eq, Generic)
instance Bin.Binary SparseMatrix

-- | Create a sparse matrix from rows
fromDocuments :: [SparseVector] -> SparseMatrix
fromDocuments rows = RowMatrix $ Vec.concat [ docToSparseRow docid doc | (docid, doc) <- zip [0..] rows ]
  where docToSparseRow docid = Vec.map (\(colid, val) -> (docid, colid, fromIntegral val)) . Vec.fromList

-- | Find the size of a sparse matrix
size :: SparseMatrix -> (Int, Int)
size (RowMatrix mat) = (maxx+1, maxy+1)
  where (maxx, maxy) = Vec.foldr' (\(r, c, _) (mr, mc) -> (max mr r, max mc c)) (-1, -1) mat
size (ColMatrix mat) = (maxx+1, maxy+1)
  where (maxx, maxy) = Vec.foldr' (\(r, c, _) (mr, mc) -> (max mr r, max mc c)) (-1, -1) mat

-- | Convert a sparse matrix into a dense one
dense :: SparseMatrix -> HMatrix.Matrix Double
-- dense (RowMatrix mat) = HMatrix.toDense $ (\(row, col, val) -> ((row, col), val)) <$> Vec.toList mat
-- dense (ColMatrix mat) = HMatrix.toDense $ (\(row, col, val) -> ((row, col), val)) <$> Vec.toList mat
dense (RowMatrix mat) =
  HMatrix.reshape width $ Vec.convert flat_mat
  where
    (width, height) = size $ RowMatrix mat
    flat_mat = Vec.update
      (Vec.replicate (width*height) 0.0) $
      Vec.map (\(r,c,v) -> (r*width+c, v)) mat
dense (ColMatrix mat) = dense $ RowMatrix mat -- a little fib but it works

-- | Convert a dense matrix into a sparse one
sparse :: HMatrix.Matrix Double -> SparseMatrix
sparse mat =
  RowMatrix $ Vec.imap (\ix v -> (
    ix `div` width,
    ix `mod` width,
    v)) $ Vec.convert $ HMatrix.flatten mat
  where (width, height) = HMatrix.size mat
  --RowMatrix $ Vec.fromList [ (ri, ci, val)
  --  | (ri, row) <- zip [0..] lmat,
  --    (ci, val) <- zip [0..] row, abs val > 1e-6 ]
  --where lmat = HMatrix.toLists mat

transpose :: SparseMatrix -> SparseMatrix
transpose (RowMatrix mat) = RowMatrix $ Vec.map (\(row, col, val) -> (col, row, val)) mat
transpose (ColMatrix mat) = ColMatrix $ Vec.map (\(col, row, val) -> (row, col, val)) mat

shift :: SparseMatrix -> SparseMatrix
shift (RowMatrix mat) = ColMatrix $ Vec.modify (Introsort.sortBy $ comparing extractCol) mat
shift (ColMatrix mat) = RowMatrix $ Vec.modify (Introsort.sortBy $ comparing extractRow) mat

extractRow (row, col, val) = row
extractCol (row, col, val) = col

-- | Multiplies a RowMatrix by a ColMatrix (converting as necessary beforehand)
-- and generates a RowMatrix as a result
mult :: SparseMatrix -> SparseMatrix -> SparseMatrix
mult left@(ColMatrix _) right@(ColMatrix _) = mult (shift left) right
mult left@(ColMatrix _) right@(RowMatrix _) = mult (shift left) (shift right)
mult left@(RowMatrix _) right@(RowMatrix _) = mult left (shift right)
mult (RowMatrix !left) (ColMatrix !right) =
  RowMatrix $ Vec.fromList [ (rowid, colid, multLines 0.0 row col)
    | (rowid, row) <- rows,
      (colid, col) <- cols]
  where
    rows = zip [0..] (getLines extractRow left)
    cols = zip [0..] (getLines extractCol right)

-- | Multiplies a RowMatrix by a ColMatrix (converting as necessary beforehand)
-- and generates a ColMatrix as a result
multCol :: SparseMatrix -> SparseMatrix -> SparseMatrix
multCol left@(ColMatrix _) right@(ColMatrix _) = multCol (shift left) right
multCol left@(ColMatrix _) right@(RowMatrix _) = multCol (shift left) (shift right)
multCol left@(RowMatrix _) right@(RowMatrix _) = multCol left (shift right)
multCol (RowMatrix left) (ColMatrix right) =
  ColMatrix $ Vec.fromList [ (rowid, colid, multLines 0.0 row col)
    | (colid, col) <- zip [0..] (getLines extractCol right),
      (rowid, row) <- zip [0..] (getLines extractRow left)]

--
--
--   INTERNAL METHODS
--
--

multLines :: Double -> [Coord] -> [Coord] -> Double
multLines x [] _ = x
multLines x _ [] = x
multLines !total xs@((_, xcol, xval):xr) ys@((yrow, _, yval):yr) =
  case compare xcol yrow of
    LT -> multLines total xr ys
    GT -> multLines total xs yr
    EQ -> multLines (total + xval * yval) xr yr

getLines :: (Coord -> Int) -> Coords -> [[Coord]]
getLines extract mat
  | Vec.null mat = []
  | otherwise = Vec.toList thisrow : getLines extract rest
  where
    (thisrow, rest) = Vec.span (\c -> extract c == target) mat
    target = extract $ Vec.head mat

-- Gets the index that a row starts at (but it works with cols too)
csrRowBeginsAt :: (Coord -> Int) -> Coords -> Vec.Vector Int
csrRowBeginsAt extract mat =
  Vec.cons 0 $ Vec.accumulate max zeros $ Vec.imap (\idx coord -> (extract coord, idx+1)) mat
  where
    zeros = Vec.replicate (end+1) 0
    end = extract $ Vec.last mat

max_stars_repo_name: SeanTater/albemarle

max_issues_repo_name: SeanTater/albemarle

max_forks_repo_name: SeanTater/albemarle


text:
--
-- Copyright (c) 2009-2011, 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.
--

module Feldspar.Core.Frontend.Fractional
where

import Data.Complex

import Feldspar.Core.Constructs
import Feldspar.Core.Constructs.Fractional

import Feldspar.Core.Frontend.Literal
import Feldspar.Core.Frontend.Num

-- | Fractional types. The relation to the standard 'Fractional' class is
--
-- @instance `Fraction` a => `Fractional` (`Data` a)@
class (Fractional a, Numeric a) => Fraction a
  where
    fromRationalFrac :: Rational -> Data a
    fromRationalFrac = value . fromRational

    divFrac :: Data a -> Data a -> Data a
    divFrac = sugarSymF DivFrac

instance Fraction Float
instance Fraction Double

instance (Fraction a, RealFloat a) => Fraction (Complex a)

instance (Fraction a) => Fractional (Data a)
  where
    fromRational = fromRationalFrac
    (/)          = divFrac

max_stars_repo_name: kffaxen/feldspar-language

max_issues_repo_name: kffaxen/feldspar-language

max_forks_repo_name: kffaxen/feldspar-language


text:
{-# LANGUAGE FlexibleContexts       #-}
{-# LANGUAGE FlexibleInstances      #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE MultiParamTypeClasses  #-}
{-# LANGUAGE OverloadedStrings      #-}
{-# LANGUAGE TemplateHaskell        #-}

module DBPnet.ReadCount
    ( readCount
    , ReadCountOpt
    , binSize
    , pValue
    , chromSize
    , varFilter
    ) where

import           Bio.Data.Bed
import Bio.Data.Bed.Utils
import           Bio.Utils.Misc                    (readDouble, readInt)
import           Conduit
import Lens.Micro ((.~), (^.))
import Lens.Micro.TH (makeFields)
import           Control.Monad                     (forM, forM_)
import           Control.Monad.Base                (liftBase)
import           Control.Monad.Morph               (hoist)
import qualified Data.ByteString.Char8             as B
import qualified Data.Conduit.Zlib                 as Zlib
import           Data.Default
import           Data.Double.Conversion.ByteString (toFixed)
import           Data.Int                          (Int32)
import           Data.List
import qualified Data.Matrix.Unboxed               as MU
import qualified Data.Text                         as T
import qualified Data.Vector.Unboxed               as U
import           Shelly                            hiding (FilePath, withTmpDir)
import           Statistics.Distribution           (complCumulative)
import           Statistics.Distribution.Poisson   (poisson)
import           Statistics.Sample
import           System.IO

import           DBPnet.Type
import           DBPnet.Utils

data ReadCountOpt = ReadCountOpt
    { readCountOptBinSize   :: !Int
    , readCountOptPValue    :: !Double
    , readCountOptChromSize :: ![(B.ByteString, Int)]
    , readCountOptVarFilter :: !Double   -- ^ remove constant signal
    } deriving (Show, Read)

makeFields ''ReadCountOpt

instance Default ReadCountOpt where
    def = ReadCountOpt
        { readCountOptBinSize = 1000
        , readCountOptPValue = 1e-2
        , readCountOptChromSize = []
        , readCountOptVarFilter  = 0.1
        }

readCount :: [Experiment] -> FilePath -> ReadCountOpt -> IO [(String, FilePath)]
readCount es outDir opt = withTmpDir outDir $ \tmp -> do
    -- count reads for each experiment
    rs <- forM es $ \e -> do
        let fls = e^.files
            targetName = e^.eid
            output = tmp ++ "/" ++ targetName
        shelly $ mkdir_p $ fromText $ T.pack output
        rc fls output
        return (targetName, output)

    combineAndFilter rs outDir
  where
    rc inputs output = do
        readcounts <- forM inputs $ \input -> do
            let fl = input^.location
                fileFormat = input^.format
            case fileFormat of
                Bed -> runResourceT $ runConduit $ streamBed fl .| countTagsBinBed (opt^.binSize) regions
                BedGZip -> runResourceT $ runConduit $ streamBedGzip fl .| countTagsBinBed (opt^.binSize) regions
                _ -> undefined

        forM_ (zip (opt^.chromSize) $ merge readcounts) $ \((chr,_), v) -> do
            let outFile = output ++ "/" ++ B.unpack chr
                l = U.length v
                bg_glob = estimateBG v
            withFile outFile WriteMode $ \handle ->
                forM_ [0..l-1] $ \i -> do
                    let getBG x = let len = x `div` 2
                                      vec = U.map (v U.!) $ U.fromList $
                                            filter (\x -> x >= 0 && x < l) $
                                            [i-len..i-1] ++ [i+1,i+len]
                                  in estimateBG vec
                        bg_local = maximum $ bg_glob : map getBG [14, 24]
                        c = v U.! i
                    case () of
                        _ | c == 0 -> B.hPutStrLn handle "0"
                          | complCumulative (poisson bg_local) c <= (opt^.pValue) ->
                                B.hPutStrLn handle $ toFixed 4 $ c / bg_local
                          | otherwise -> B.hPutStrLn handle "0"
      where
        regions = map (\(chr,s) -> BED3 chr 0 s) $ opt^.chromSize
        estimateBG xs = let (m, var) = meanVarianceUnb xs
                        in mean $ U.filter (< (m + 2 * sqrt var)) xs
        merge rs = map f $ transpose counts
          where
            f = foldl1' (\acc x -> U.zipWith (+) acc x)
            counts = flip map rs $ \(values, _) ->
                     map (U.map (fromIntegral :: Int32 -> Double)) values
-- we do not need to normalize or average because we don't use input/control
--            n = fromIntegral $ length rs

    combineAndFilter xs output = do
        shelly $ mkdir_p $ fromText $ T.pack output
        let (ids, dirs)  = unzip xs

        dat <- forM (opt^.chromSize) $ \(chr', _) -> do
            let chr = B.unpack chr'
            rs <- forM dirs $ \dir -> do
                c <- B.readFile $ dir ++ "/" ++ chr
                return $ U.fromList $ map readDouble $ B.lines c
            let (idx, rc) = unzip $ filter (highCV . snd) $
                            zip [0, opt^.binSize ..] $
                            MU.toRows $ MU.fromColumns rs
            return (zip (repeat chr') idx, rc)

        let idx = concat $ fst $ unzip dat
            dat' = MU.fromColumns $ concat $ snd $ unzip dat
        forM (zip ids $ MU.toRows dat') $ \(i, v) -> do
            let outFile = output ++ "/" ++ i ++ ".bed"
            writeBed outFile $ zipWith toBed idx $ U.toList v
            return (i, outFile)
      where
        highCV v = let (_,var) = meanVarianceUnb v in sqrt var >= opt^.varFilter
        toBed (chr,i) x = BEDGraph chr i (i+(opt^.binSize)) x

max_stars_repo_name: kaizhang/DBPnet

max_issues_repo_name: kaizhang/DBPnet

max_forks_repo_name: kaizhang/DBPnet


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE TypeSynonymInstances #-}

#if __GLASGOW_HASKELL__ >= 701
{-# LANGUAGE DefaultSignatures #-}
#endif

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

#if __GLASGOW_HASKELL__ >= 710
{-# LANGUAGE Safe #-}
#elif __GLASGOW_HASKELL__ >= 701
{-# LANGUAGE Trustworthy #-}
#endif

module Generics.Deriving.Foldable (
  -- * Generic Foldable class
    GFoldable(..)

  -- * Default method
  , gfoldMapdefault

  -- * Derived functions
  , gtoList
  , gconcat
  , gconcatMap
  , gand
  , gor
  , gany
  , gall
  , gsum
  , gproduct
  , gmaximum
  , gmaximumBy
  , gminimum
  , gminimumBy
  , gelem
  , gnotElem
  , gfind

  -- * Internal Foldable class
  , GFoldable'(..)
  ) where

import           Control.Applicative (Const, ZipList)

import           Data.Maybe
import qualified Data.Monoid as Monoid (First, Last, Product(..), Sum(..))
import           Data.Monoid (All(..), Any(..), Dual(..), Endo(..))
#if !(MIN_VERSION_base(4,8,0))
import           Data.Monoid (Monoid(..))
#endif

import           Generics.Deriving.Base

#if MIN_VERSION_base(4,4,0)
import           Data.Complex (Complex)
#endif

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

#if MIN_VERSION_base(4,7,0)
import           Data.Proxy (Proxy)
#endif

#if MIN_VERSION_base(4,8,0)
import           Data.Functor.Identity (Identity)
#endif

#if MIN_VERSION_base(4,9,0)
import qualified Data.Functor.Product as Functor (Product)
import qualified Data.Functor.Sum as Functor (Sum)
import           Data.List.NonEmpty (NonEmpty)
import qualified Data.Semigroup as Semigroup (First, Last)
import           Data.Semigroup (Arg, Max, Min, WrappedMonoid)
#endif

--------------------------------------------------------------------------------
-- Generic fold
--------------------------------------------------------------------------------

class GFoldable' t where
  gfoldMap' :: Monoid m => (a -> m) -> t a -> m

instance GFoldable' V1 where
  gfoldMap' _ _ = mempty

instance GFoldable' U1 where
  gfoldMap' _ U1 = mempty

instance GFoldable' Par1 where
  gfoldMap' f (Par1 a) = f a

instance GFoldable' (K1 i c) where
  gfoldMap' _ (K1 _) = mempty

instance (GFoldable f) => GFoldable' (Rec1 f) where
  gfoldMap' f (Rec1 a) = gfoldMap f a

instance (GFoldable' f) => GFoldable' (M1 i c f) where
  gfoldMap' f (M1 a) = gfoldMap' f a

instance (GFoldable' f, GFoldable' g) => GFoldable' (f :+: g) where
  gfoldMap' f (L1 a) = gfoldMap' f a
  gfoldMap' f (R1 a) = gfoldMap' f a

instance (GFoldable' f, GFoldable' g) => GFoldable' (f :*: g) where
  gfoldMap' f (a :*: b) = mappend (gfoldMap' f a) (gfoldMap' f b)

instance (GFoldable f, GFoldable' g) => GFoldable' (f :.: g) where
  gfoldMap' f (Comp1 x) = gfoldMap (gfoldMap' f) x

instance GFoldable' UAddr where
  gfoldMap' _ (UAddr _) = mempty

instance GFoldable' UChar where
  gfoldMap' _ (UChar _) = mempty

instance GFoldable' UDouble where
  gfoldMap' _ (UDouble _) = mempty

instance GFoldable' UFloat where
  gfoldMap' _ (UFloat _) = mempty

instance GFoldable' UInt where
  gfoldMap' _ (UInt _) = mempty

instance GFoldable' UWord where
  gfoldMap' _ (UWord _) = mempty

class GFoldable t where
  gfoldMap :: Monoid m => (a -> m) -> t a -> m
#if __GLASGOW_HASKELL__ >= 701
  default gfoldMap :: (Generic1 t, GFoldable' (Rep1 t), Monoid m)
                   => (a -> m) -> t a -> m
  gfoldMap = gfoldMapdefault
#endif

  gfold :: Monoid m => t m -> m
  gfold = gfoldMap id

  gfoldr :: (a -> b -> b) -> b -> t a -> b
  gfoldr f z t = appEndo (gfoldMap (Endo . f) t) z

  gfoldr' :: (a -> b -> b) -> b -> t a -> b
  gfoldr' f z0 xs = gfoldl f' id xs z0
    where f' k x z = k $! f x z

  gfoldl :: (a -> b -> a) -> a -> t b -> a
  gfoldl f z t = appEndo (getDual (gfoldMap (Dual . Endo . flip f) t)) z

  gfoldl' :: (a -> b -> a) -> a -> t b -> a
  gfoldl' f z0 xs = gfoldr f' id xs z0
    where f' x k z = k $! f z x

  gfoldr1 :: (a -> a -> a) -> t a -> a
  gfoldr1 f xs = fromMaybe (error "gfoldr1: empty structure")
                   (gfoldr mf Nothing xs)
    where
      mf x Nothing = Just x
      mf x (Just y) = Just (f x y)

  gfoldl1 :: (a -> a -> a) -> t a -> a
  gfoldl1 f xs = fromMaybe (error "foldl1: empty structure")
                   (gfoldl mf Nothing xs)
    where
      mf Nothing y = Just y
      mf (Just x) y = Just (f x y)

gfoldMapdefault :: (Generic1 t, GFoldable' (Rep1 t), Monoid m)
                => (a -> m) -> t a -> m
gfoldMapdefault f x = gfoldMap' f (from1 x)

-- Base types instances
instance GFoldable ((,) a) where
  gfoldMap = gfoldMapdefault

instance GFoldable [] where
  gfoldMap = gfoldMapdefault

#if MIN_VERSION_base(4,9,0)
instance GFoldable (Arg a) where
  gfoldMap = gfoldMapdefault
#endif

#if MIN_VERSION_base(4,4,0)
instance GFoldable Complex where
  gfoldMap = gfoldMapdefault
#endif

instance GFoldable (Const m) where
  gfoldMap = gfoldMapdefault

instance GFoldable Down where
  gfoldMap = gfoldMapdefault

instance GFoldable Dual where
  gfoldMap = gfoldMapdefault

instance GFoldable (Either a) where
  gfoldMap = gfoldMapdefault

instance GFoldable Monoid.First where
  gfoldMap = gfoldMapdefault

#if MIN_VERSION_base(4,9,0)
instance GFoldable (Semigroup.First) where
  gfoldMap = gfoldMapdefault
#endif

#if MIN_VERSION_base(4,8,0)
instance GFoldable Identity where
  gfoldMap = gfoldMapdefault
#endif

instance GFoldable Monoid.Last where
  gfoldMap = gfoldMapdefault

#if MIN_VERSION_base(4,9,0)
instance GFoldable Semigroup.Last where
  gfoldMap = gfoldMapdefault

instance GFoldable Max where
  gfoldMap = gfoldMapdefault
#endif

instance GFoldable Maybe where
  gfoldMap = gfoldMapdefault

#if MIN_VERSION_base(4,9,0)
instance GFoldable Min where
  gfoldMap = gfoldMapdefault

instance GFoldable NonEmpty where
  gfoldMap = gfoldMapdefault
#endif

instance GFoldable Monoid.Product where
  gfoldMap = gfoldMapdefault

#if MIN_VERSION_base(4,9,0)
instance (GFoldable f, GFoldable g) => GFoldable (Functor.Product f g) where
  gfoldMap = gfoldMapdefault
#endif

#if MIN_VERSION_base(4,7,0)
instance GFoldable Proxy where
  gfoldMap = gfoldMapdefault
#endif

instance GFoldable Monoid.Sum where
  gfoldMap = gfoldMapdefault

#if MIN_VERSION_base(4,9,0)
instance (GFoldable f, GFoldable g) => GFoldable (Functor.Sum f g) where
  gfoldMap = gfoldMapdefault

instance GFoldable WrappedMonoid where
  gfoldMap = gfoldMapdefault
#endif

instance GFoldable ZipList where
  gfoldMap = gfoldMapdefault

gtoList :: GFoldable t => t a -> [a]
gtoList = gfoldr (:) []

gconcat :: GFoldable t => t [a] -> [a]
gconcat = gfold

gconcatMap :: GFoldable t => (a -> [b]) -> t a -> [b]
gconcatMap = gfoldMap

gand :: GFoldable t => t Bool -> Bool
gand = getAll . gfoldMap All

gor :: GFoldable t => t Bool -> Bool
gor = getAny . gfoldMap Any

gany :: GFoldable t => (a -> Bool) -> t a -> Bool
gany p = getAny . gfoldMap (Any . p)

gall :: GFoldable t => (a -> Bool) -> t a -> Bool
gall p = getAll . gfoldMap (All . p)

gsum :: (GFoldable t, Num a) => t a -> a
gsum = Monoid.getSum . gfoldMap Monoid.Sum

gproduct :: (GFoldable t, Num a) => t a -> a
gproduct = Monoid.getProduct . gfoldMap Monoid.Product

gmaximum :: (GFoldable t, Ord a) => t a -> a
gmaximum = gfoldr1 max

gmaximumBy :: GFoldable t => (a -> a -> Ordering) -> t a -> a
gmaximumBy cmp = gfoldr1 max'
  where max' x y = case cmp x y of
                        GT -> x
                        _  -> y

gminimum :: (GFoldable t, Ord a) => t a -> a
gminimum = gfoldr1 min

gminimumBy :: GFoldable t => (a -> a -> Ordering) -> t a -> a
gminimumBy cmp = gfoldr1 min'
  where min' x y = case cmp x y of
                        GT -> y
                        _  -> x

gelem :: (GFoldable t, Eq a) => a -> t a -> Bool
gelem = gany . (==)

gnotElem :: (GFoldable t, Eq a) => a -> t a -> Bool
gnotElem x = not . gelem x

gfind :: GFoldable t => (a -> Bool) -> t a -> Maybe a
gfind p = listToMaybe . gconcatMap (\ x -> if p x then [x] else [])

max_stars_repo_name: treeowl/generic-deriving

max_issues_repo_name: treeowl/generic-deriving

max_forks_repo_name: treeowl/generic-deriving


text:
{-# LANGUAGE ViewPatterns #-}
module EigenDFT where

import           Control.Monad                      as M
import           Control.Monad.Parallel             as MP
import qualified Data.Array.Accelerate              as A
import qualified Data.Array.Accelerate.Data.Complex as A
import           Data.Array.Accelerate.LLVM.PTX
import           Data.Array.Repa                    as R
import           Data.List                          as L
import           Foreign.CUDA.Driver                as CUDA
import           Image.IO
import           Numeric.LinearAlgebra              as NL
import           Pinwheel.Base
import           Pinwheel.FourierSeries2D
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf

main = do
  args@(nStr:numFreqStr:deltaStr:_) <- getArgs
  let n = read nStr :: Int
      numFreq = read numFreqStr :: Int
      delta = read deltaStr :: Double
      len = n ^ 2
      center = div n 2
      centerFreq = div numFreq 2
      constant = (-2) * pi * delta / (fromIntegral n) :: Double
      matArr' =
        fromFunction (Z :. n :. n :. numFreq :. numFreq) $ \(Z :. x' :. y' :. wx' :. wy') ->
          let wx = wx' - centerFreq
              wy = wy' - centerFreq
              x = x' - center
              y = y' - center
          in cis $ constant * fromIntegral (wx * x + wy * y)
      pinwheelVec =
        NL.fromList . R.toList . fromFunction (Z :. n :. n) $ \(Z :. x' :. y') ->
          let x = x' - center
              y = y' - center
          in fourierMellin (0.5 :: Double) 5 5 (fromIntegral x, fromIntegral y)
      folderPath = "output/test/EigenDFT"
  when (numFreq > n) (error "error: numFreq > n")
  removePathForcibly folderPath
  createDirectoryIfMissing True folderPath
  matArr <- computeUnboxedP matArr'
  let mat = (len >< (numFreq ^ 2)) . R.toList $ matArr
      coefficients = pinwheelVec <# mat
  --     xs =
  --       if numFreq == n
  --         then let (eigVal, eigVec) = eig mat
  --              in snd . L.unzip . L.reverse . L.sortOn fst $
  --                 L.zip
  --                   (L.map magnitude . NL.toList $ eigVal)
  --                   (L.map NL.toList . NL.toColumns $ eigVec)
  --         else let (u, s, v) = thinSVD mat
  --              in snd . L.unzip . L.reverse . L.sortOn fst $
  --                 L.zip (NL.toList $ s) (L.map NL.toList . NL.toColumns $ u)
  -- M.mapM_
  --   (\(i, ys) ->
  --      plotImageRepaComplex (folderPath </> printf "%d.png" i) .
  --      ImageRepa 8 . fromListUnboxed (Z :. (1 :: Int) :. n :. n) $
  --      ys) $
  --   L.zip [0 :: Int .. (numFreq ^ 2 - 1)] xs
  plotImageRepaComplex (folderPath </> "Pinwheel.png" ) . ImageRepa 8 . fromListUnboxed (Z :. (1 :: Int) :. n :. n) . NL.toList $ pinwheelVec
  plotImageRepaComplex (folderPath </> "coefficients.png" ) . ImageRepa 8 . fromListUnboxed (Z :. (1 :: Int) :. numFreq :. numFreq) . NL.toList $ coefficients
  -- let centerAcc = A.constant center
  --     centerFreqAcc = A.constant centerFreq
  --     constantAcc = A.constant constant
  --     matArrAcc' =
  --       A.generate (A.constant (A.Z A.:. numFreq A.:. numFreq A.:. n A.:. n)) $ \(A.unlift -> (A.Z A.:. wx' A.:. wy' A.:. x' A.:. y')) ->
  --         let wx = wx' - centerFreqAcc
  --             wy = wy' - centerFreqAcc
  --             x = x' - centerAcc
  --             y = y' - centerAcc
  --         in A.cis $ constantAcc * A.fromIntegral (wx * x + wy * y)
  -- initialise []
  -- dev <- device 7
  -- ctx <- CUDA.create dev []
  -- ptx <- createTargetFromContext ctx
  -- let matArrAcc = runWith ptx matArrAcc'
  --     mat = ((numFreq^2) >< len) . A.toList $ matArrAcc

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:
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE BangPatterns #-}

module Statistics.Quantile.Exact (
    external
  , external'
  , countOut
) where

import Control.Applicative
import Control.Monad.IO.Class

import Data.Monoid

import Debug.Trace

import Numeric

import Pipes
import qualified Pipes.Prelude as P

import Statistics.Quantile.Types
import Statistics.Quantile.Util

import System.Directory
import System.Exit
import System.IO
import System.IO.Temp
import System.FilePath
import System.Process

external :: Selector IO
external =  Selector external'

external' :: Quantile
          -> Stream IO
          -> IO Double
external' q src = do
  (fp, fh) <- openTempFile tmpDir "quantile.txt"
  n <- countOut fh src
  hClose fh
  fp' <- externalMergeSort fp
  fh' <- openFile fp' ReadMode
  quant <- select (qidx q n) $
                P.fromHandle fh'
            >-> P.map realOrDie
  removeFile fp >> removeFile fp'
  pure quant
  where realOrDie :: String -> Double
        realOrDie = read

countOut :: Handle
         -> Stream IO
         -> IO Int
countOut fh (Stream src) =
  P.foldM (incrementWrite fh) (pure 0) pure src
  where incrementWrite h n v = do
          -- Can't use the Show instance because scientific notation
          -- would break sorting.
          hPutStrLn h (showFFloat Nothing v "")
          pure $! (n + 1)

select :: Monad m
       => Int
       -> Producer Double m ()
       -> m Double
select n p0 = next p0 >>= \case
  (Left _)       -> error $ show EmptySample
  (Right (v, p)) -> go 1 v p
   where go i v p
           | i == n    = pure v
           | otherwise =
               next p >>= \case
                 (Left _) -> error "universe broke"
                 (Right (v', p')) -> go (i + 1) v' p'

externalMergeSort :: FilePath -> IO FilePath
externalMergeSort fp =
  let fp' = fp <> ".sorted" in do
  fh <- openFile fp' WriteMode
  ph <- runProcess "/usr/bin/sort" sortOpts Nothing Nothing Nothing (Just fh) Nothing
  st <- waitForProcess ph
  hClose fh
  case st of
    ExitSuccess -> pure fp'
    _ -> error $ "sort -n " <> fp <> " failed"
  where
    -- Numeric, TMPDIR=/var/tmp so we don't fill up a tmpfs.
    sortOpts = [ "-n", "-T", "/var/tmp", fp]

max_stars_repo_name: fractalcat/slides

max_issues_repo_name: fractalcat/slides

max_forks_repo_name: fractalcat/slides


text:
module Fungine.Graphics.Texture where
import Fungine.Prelude
import Fungine.Graphics.RGBA
import Numeric.LinearAlgebra
data Texture = StaticImage (Matrix RGBA)

max_stars_repo_name: jmdurr/fungine

max_issues_repo_name: jmdurr/fungine

max_forks_repo_name: jmdurr/fungine


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

-- Created September 20 2015

-- TODO | -
--        -

-- SPEC | -
--        -

--------------------------------------------------------------------------------------------------------------------------------------------
-- GHC Pragmas
--------------------------------------------------------------------------------------------------------------------------------------------
{-# LANGUAGE Rank2Types #-}

--------------------------------------------------------------------------------------------------------------------------------------------
-- API
--------------------------------------------------------------------------------------------------------------------------------------------
module Occlusion.Events where

--------------------------------------------------------------------------------------------------------------------------------------------
-- We'll need these
--------------------------------------------------------------------------------------------------------------------------------------------
import Data.IORef
import Data.Complex
import Data.Functor
import Data.Maybe
import qualified Data.Text as T
import qualified Data.Set  as S
import Control.Lens
import Control.Monad (when, unless, forM, mapM)

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

import Occlusion.Types
import Occlusion.Lenses
import Occlusion.Vector
import qualified Occlusion.Behaviours as Behaviours
import qualified Occlusion.Render     as Render

--------------------------------------------------------------------------------------------------------------------------------------------
-- Event listeners
--------------------------------------------------------------------------------------------------------------------------------------------
-- |
onrender :: AppState -> Cairo.Render ()
onrender appstate = do
  Render.background appstate
  Render.scene (appstate^.scene)

-- |
-- TODO: Simulation speed
-- TODO: Factour out 'tick'
-- TODO: Use state monad
onanimate :: IORef AppState -> IO Bool
onanimate stateref = do
  --
  appstate <- readIORef stateref

  --
  widgetQueueDraw (appstate^.gui.canvas)
  modifyIORef stateref (animation.frame %~ (+1)) -- Increment frame count

  -- TODO: Refactor (cf. State, traverse, each)
  -- modifyIORef stateref $ \appstate -> forM [scene.player, scene.npcs.ix 0] $ \char -> Behaviours.run char
  modifyIORef stateref $ Behaviours.run (scene.player)
  modifyIORef stateref $ Behaviours.run (scene.npcs.nth 0)
  return True

-- |
onmousemoves :: IORef AppState -> EventM EMotion Bool
onmousemoves stateref = do
  p <- uncurry (:+) <$> eventCoordinates
  Cairo.liftIO $ do
    modifyIORef stateref (input.mouse .~ p)
  return False

-- |
onmousepressed :: IORef AppState -> EventM EButton Bool
onmousepressed stateref = do
  p <- uncurry (:+) <$> eventCoordinates
  Cairo.liftIO $ do
    modifyIORef stateref (input.click .~ Just p)
  return False

-- |
onmousereleased :: IORef AppState -> EventM EButton Bool
onmousereleased stateref = do
  p <- uncurry (:+) <$> eventCoordinates
  Cairo.liftIO $ do
    modifyIORef stateref (input.click .~ Nothing)
  return False

-- |
onkeypressed :: IORef AppState -> EventM EKey Bool
onkeypressed stateref = do
  key <- T.unpack <$> eventKeyName
  Cairo.liftIO $ do
    appstate <- readIORef stateref
    unless (S.member key $ appstate^.input.keyboard) $ do
      modifyIORef stateref (input.keyboard %~ S.insert key)
  return False

-- |
onkeyreleased :: IORef AppState -> EventM EKey Bool
onkeyreleased stateref = do
  key <- T.unpack <$> eventKeyName
  Cairo.liftIO $ do
    appstate <- readIORef stateref
    modifyIORef stateref (input.keyboard %~ S.delete key)
  return False

-- |
-- onwindowresized
-- onwindowminimised

--------------------------------------------------------------------------------------------------------------------------------------------
-- Utilities
--------------------------------------------------------------------------------------------------------------------------------------------
-- |
attach :: Window -> DrawingArea -> IORef AppState -> IO ()
attach window canvas stateref = do
  appstate <- readIORef stateref

  canvas `on` draw $ Cairo.liftIO (readIORef stateref) >>= onrender

  canvas `on` motionNotifyEvent  $ onmousemoves    stateref
  canvas `on` buttonPressEvent   $ onmousepressed  stateref
  canvas `on` buttonReleaseEvent $ onmousereleased stateref
  window `on` keyPressEvent      $ onkeypressed    stateref
  window `on` keyReleaseEvent    $ onkeyreleased   stateref

  window `on` deleteEvent $ Cairo.liftIO mainQuit >> return False

  timeoutAdd (onanimate stateref) (round $ 1000.0 / appstate^.animation.fps)

  pass

--------------------------------------------------------------------------------------------------------------------------------------------
-- |
velocityFromKey :: Double -> String -> Complex Double
velocityFromKey speed key = case key of
      "Left"  -> (-speed):+0
      "Right" -> ( speed):+0
      "Up"    ->       0 :+(-speed)
      "Down"  ->       0 :+( speed)
      _       ->       0:+0

max_stars_repo_name: SwiftsNamesake/Occlusion

max_issues_repo_name: SwiftsNamesake/Occlusion

max_forks_repo_name: SwiftsNamesake/Occlusion


text:
module Net
  (
    repTrain
  , buildNetwork
  , forward
  , mse
  , sigmoid
  , sigmoid'
  , tanh'
  , printNet
  , NN
  ) where

import Numeric.LinearAlgebra
import Data.Traversable

type ActivationF = Double -> Double
type ActivationD = Double -> Double
type Activation = (ActivationF, ActivationD)
type Bias = Matrix R
type Weights = Matrix R
type Layer = (Weights, Bias, Activation)
type NN = [Layer]

type Input = Matrix R
type Output = Matrix R
type Deriv = Matrix R
type Error = Matrix R
type LearningRate = Double

repTrain :: Int -> LearningRate -> Input -> Output -> NN -> NN
repTrain 0 _ _ _ nn = nn
repTrain n lr x y nn = repTrain (n-1) lr x y (backprop lr x y nn)

buildNetwork :: Activation -> [Int] -> IO NN
buildNetwork act (a:b:xs) = do
  weights <- randn a b
  bias    <- randn 1 b
  rest    <- buildNetwork act (b:xs)
  return $ (weights, bias, act):rest
buildNetwork _ _ = return []

forward :: NN -> Input -> Output
forward = (fst .) . flip feedForward

feedForward :: Input -> NN -> (Output, [(Input, Deriv)])
feedForward = mapAccumL activateLayer

activateLayer :: Input -> Layer -> (Output, (Input, Deriv))
activateLayer input (wh, bh, act) = (output, (input, deriv))
  where
    output = cmap (fst act) ((input <> wh + bh))
    deriv  = cmap (snd act) output

backprop :: LearningRate -> Input -> Output -> NN -> NN
backprop lr input y net = backpropLayers lr error inpDs net
   where
     (output, inpDs) = feedForward input net
     error = y - output

mse :: NN -> Input -> Output -> Error
mse net x y = let (output, _) = feedForward x net in y - output

backpropLayers :: LearningRate -> Error -> [(Input, Deriv)] -> NN -> NN
backpropLayers lr error inpDs = snd . mapAccumR (backpropMatrix lr) error . zip inpDs
--  = ((snd . mapAccumR (backpropMatrix lr) error .) . zip) inpDs

backpropMatrix :: LearningRate -> Matrix R -> ((Matrix R, Matrix R), Layer) -> (Error, Layer)
backpropMatrix lr error (inputInfo, layer) = (prev_error, new_layer)
  where
    (layer_weights, layer_bias, act) = layer
    (layer_inp, layer_deriv) = inputInfo
    new_layer = (new_weights, new_bias, act)

    d_layer = layer_deriv * error
    new_weights = (+) layer_weights $ cmap (* lr) $ (tr layer_inp) <> d_layer
    new_bias = (+) layer_bias $ cmap (* lr) $ asRow $ vector $ map sumElements $ toColumns d_layer

    prev_error = d_layer <> (tr layer_weights)

sigmoid :: Double -> Double
sigmoid x = 1.0 / (1 + exp (-x))

sigmoid' :: Double -> Double
sigmoid' x = x * (1 - x)

tanh' :: Double -> Double
tanh' x = let fx = tanh x in 1 - fx * fx

printNet :: NN -> IO ()
printNet layers = do
  mapM_ (\(w,b,_) -> do
            disp 3 w
            disp 3 b) layers


text:
module InterpSpec where

import Commons
import Data.Complex (Complex (..))
import Data.Map.Strict (union)
import Data.Maybe (fromJust)
import Debug.Trace (traceShowId)
import HashedExpression.Internal.Expression
import HashedExpression.Internal.Normalize
import HashedExpression.Internal.Utils
import HashedExpression.Interp
import HashedExpression.Operation hiding (product, sum)
import qualified HashedExpression.Operation
import HashedExpression.Prettify
import Test.Hspec
import Test.QuickCheck (property)
import Var
import Prelude hiding ((^))

-- |
prop_AddScalarR :: SuiteScalarR -> SuiteScalarR -> Bool
prop_AddScalarR (Suite exp1 valMaps1) (Suite exp2 valMaps2) =
  eval valMaps (exp1 + exp2) == eval valMaps exp1 + eval valMaps exp2
  where
    valMaps = valMaps1 `union` valMaps2

-- |
prop_MultiplyScalarR :: SuiteScalarR -> SuiteScalarR -> Bool
prop_MultiplyScalarR (Suite exp1 valMaps1) (Suite exp2 valMaps2) =
  eval valMaps (exp1 * exp2) == eval valMaps exp1 * eval valMaps exp2
  where
    valMaps = valMaps1 `union` valMaps2

-- |
prop_AddScalarC :: SuiteScalarC -> SuiteScalarC -> Bool
prop_AddScalarC (Suite exp1 valMaps1) (Suite exp2 valMaps2) =
  eval valMaps (exp1 + exp2) == eval valMaps exp1 + eval valMaps exp2
  where
    valMaps = valMaps1 `union` valMaps2

-- |
prop_MultiplyScalarC :: SuiteScalarC -> SuiteScalarC -> Bool
prop_MultiplyScalarC (Suite exp1 valMaps1) (Suite exp2 valMaps2) =
  eval valMaps (exp1 * exp2) == eval valMaps exp1 * eval valMaps exp2
  where
    valMaps = valMaps1 `union` valMaps2

-- |
prop_RotateOneR1 :: SuiteOneR -> Bool
prop_RotateOneR1 (Suite exp valMaps) =
  eval valMaps (rotate 0 exp) == eval valMaps exp

-- |
prop_RotateOneR2 :: SuiteOneR -> Int -> Bool
prop_RotateOneR2 (Suite exp valMaps) amount =
  eval valMaps (f exp) == eval valMaps exp
  where
    f = rotate amount . rotate (- amount)

-- |
prop_RotateOneR3 :: SuiteOneR -> Int -> Int -> Bool
prop_RotateOneR3 (Suite exp valMaps) amount1 amount2 =
  eval valMaps (f1 exp) == eval valMaps (f2 exp)
  where
    f1 = rotate amount1 . rotate amount2
    f2 = rotate (amount1 + amount2)

-- |
prop_RotateTwoR1 :: SuiteTwoR -> Bool
prop_RotateTwoR1 (Suite exp valMaps) =
  eval valMaps (rotate (0, 0) exp) == eval valMaps exp

-- |
prop_RotateTwoR2 :: SuiteTwoR -> (Int, Int) -> Bool
prop_RotateTwoR2 (Suite exp valMaps) (offset1, offset2) =
  eval valMaps (f exp) == eval valMaps exp
  where
    f = rotate (offset1, offset2) . rotate (- offset1, - offset2)

-- |
prop_RotateTwoR3 :: SuiteTwoR -> (Int, Int) -> (Int, Int) -> Bool
prop_RotateTwoR3 (Suite exp valMaps) amount1 amount2 =
  eval valMaps (f1 exp) == eval valMaps (f2 exp)
  where
    f1 = rotate amount1 . rotate amount2
    f2 = rotate (fst amount1 + fst amount2, snd amount1 + snd amount2)

spec :: Spec
spec =
  describe "Interp spec" $ do
    specify "prop_Add Scalar R" $ property prop_AddScalarR
    specify "prop_Multiply Scalar R" $ property prop_MultiplyScalarR
    specify "prop_Add Scalar C" $ property prop_AddScalarC
    specify "prop_Multiply Scalar C" $ property prop_MultiplyScalarC
    specify "prop_Rotate One R rotate 0 should stay the same" $
      property prop_RotateOneR1
    specify "prop_Rotate One R rotate a and -a should stay the same" $
      property prop_RotateOneR2
    specify
      "prop_Rotate One R rotate a then rotate b should equal rotate (a + b)"
      $ property prop_RotateOneR3
    specify "prop_Rotate Two R rotate (0, 0) should stay the same" $
      property prop_RotateTwoR1
    specify "prop_Rotate Two R rotate a and -a should stay the same" $
      property prop_RotateTwoR2
    specify
      "prop_Rotate Two R rotate a then rotate b should equal rotate (a + b)"
      $ property prop_RotateTwoR3

max_stars_repo_name: Turboscient/HashedExpression

max_issues_repo_name: Turboscient/HashedExpression

max_forks_repo_name: Turboscient/HashedExpression


text:
{-# LANGUAGE CPP               #-}
{-# LANGUAGE DeriveAnyClass    #-}
{-# LANGUAGE FlexibleContexts  #-}
{-# LANGUAGE OverloadedStrings #-}

-- | Marshal a limited subset of J arrays into Repa arrays.
--
-- = Tutorial
--
-- Suppose we wish to perform linear regression. In J we could do:
--
-- @
-- xs := 1 2 3
-- ys := 2 4 6
--
-- reg_result =: ys %. xs ^/ i.2
-- @
--
-- To do this with Haskell data:
--
-- @
-- do
--    jenv <- 'jinit' 'libLinux'
--
--    let hsArr0 = R.fromListUnboxed (R.ix1 3) [1.0,2.0,3.0]
--        hsArr1 = R.fromListUnboxed (R.ix1 3) [2.0,4.0,6.0]
--        jArr0 = 'JDoubleArr' $ R.copyS $ R.map (realToFrac :: Double -> CDouble) hsArr0
--        jArr1 = 'JDoubleArr' $ R.copyS $ R.map (realToFrac :: Double -> CDouble) hsArr1
--
--    'setJData' jenv "xs" jArr0
--    'setJData' jenv "ys" jArr1
--
--    'bsDispatch' jenv "reg_result =: ys %. xs ^/ i.2"
--
--    'JDoubleArr' res <- 'getJData' jenv "reg_result"
--    R.toList res
-- @
--
-- There are three steps to do the calculation, plus one to get a J environment.
--
--     (1) Use 'jinit' with the appropriate file path for your platform
--
--     2. Marshal Haskell values and send them to the J environment. To do so, we
--     use 'setJData', which takes a 'JData' containing a repa array or
--     a string.
--
--     3. Perform calculations within the J environment. Here, we use
--     'bsDispatch' to compute some results and assign them within J
--
--     4. Marshal J values back to Haskell. We use 'getJData'.
--
--
--  Since marshaling data between J and Haskell is expensive, it's best to do as
--  much computation as possible in J.
--
--  == Loading Profile
--
--  If you would like to use user libraries, you need to use 'jLoad' on the
--  'JEnv'. As an example:
--
--  @
--  do
--      jenv <- 'jinit' 'libLinux'
--      'jLoad' jenv ('linuxProfile' "9.02")
--      'bsDispatch' 'jenv' "load'tables/csv'"
--  @
--
--  This will load the CSV addon, assuming it is installed.
--
--  = FFI
--
--  If you want to marshal data yourself, say to use a @Vector@, look at 'JEnv'.
module Language.J ( -- * Environment
                    JEnv (..)
                  , jinit
                  , jLoad
                  , Profile (..)
                  , linuxProfile
                  , macProfile
                  , windowsProfile
#ifndef mingw32_HOST_OS
                  , libLinux
                  , libMac
                  , profLinux
#else
                  , libWindows
#endif
                  , bsDispatch
                  , bsOut
                  , JVersion
                  -- * Repa
                  , JData (..)
                  , getJData
                  , setJData
                  -- * Exceptions
                  , JErr (..)
                  -- * Vector helpers
                  , tryIntVect
                  , copyIntVect
                  -- * FFI
                  , J
                  , JDoType
                  , JGetMType
                  , JGetRType
                  , JSetAType
                  ) where

import           Control.Applicative             (pure, (<$>), (<*>))
import           Control.Exception               (Exception, throw)
import qualified Data.Array.Repa                 as R
import qualified Data.Array.Repa.Repr.ForeignPtr as RF
import qualified Data.ByteString                 as BS
import qualified Data.ByteString.Char8           as ASCII
import qualified Data.ByteString.Internal        as BS
import           Data.Complex                    (Complex (..))
import           Data.Functor                    (void)
import           Data.Semigroup                  ((<>))
import qualified Data.Vector.Unboxed             as V
import           Foreign.C.String                (CString)
import           Foreign.C.Types                 (CChar, CDouble, CInt (..), CLLong (..))
import           Foreign.ForeignPtr              (ForeignPtr, castForeignPtr, mallocForeignPtrBytes, withForeignPtr)
import           Foreign.Marshal                 (alloca, copyArray, mallocBytes, peekArray, pokeArray)
import           Foreign.Ptr                     (FunPtr, Ptr, castPtrToFunPtr, plusPtr)
import           Foreign.Storable                (Storable, peek, pokeByteOff, sizeOf)
import           System.Info                     (arch)
#ifndef mingw32_HOST_OS
import           System.Posix.ByteString         (RTLDFlags (RTLD_LAZY), RawFilePath, dlopen, dlsym)
#else
import           System.Win32.DLL                (getProcAddress, loadLibrary)
#endif

-- Upstream reference
-- https://github.com/jsoftware/stats_jserver4r/blob/4c94fc6df351fab34791aa9d78d158eaefd33b17/source/lib/j2r.c
-- https://github.com/jsoftware/stats_jserver4r/blob/4c94fc6df351fab34791aa9d78d158eaefd33b17/source/lib/r2j.c

-- | Abstract context
data J

data JEnv = JEnv { context   :: Ptr J
                 , evaluator :: JDoType
                 , reader    :: JGetMType
                 , out       :: JGetRType
                 , setter    :: JSetAType
                 }

type JDoType = Ptr J -> CString -> IO CInt
type JGetMType = Ptr J -> CString -> Ptr CLLong -> Ptr CLLong -> Ptr (Ptr CLLong) -> Ptr (Ptr CChar) -> IO CInt
type JGetRType = Ptr J -> IO CString
type JSetAType = Ptr J -> CLLong -> CString -> CLLong -> Ptr () -> IO CInt

foreign import ccall "dynamic" mkJDo :: FunPtr JDoType -> JDoType
foreign import ccall "dynamic" mkJInit :: FunPtr (IO (Ptr J)) -> IO (Ptr J)
foreign import ccall "dynamic" mkJGetM :: FunPtr JGetMType -> JGetMType
foreign import ccall "dynamic" mkJGetR :: FunPtr JGetRType -> JGetRType
foreign import ccall "dynamic" mkJSetA :: FunPtr JSetAType -> JSetAType

type JVersion = [Int]

squashVersion :: JVersion -> String
squashVersion = concatMap show

squashVersionBS :: JVersion -> BS.ByteString
squashVersionBS = ASCII.pack . squashVersion

#ifndef mingw32_HOST_OS
-- | Expected 'RawFilePath' to the library on a Linux machine.
libLinux :: RawFilePath
libLinux = "/usr/lib/" <> ASCII.pack arch <> "-linux-gnu/libj.so"

-- | Expected 'RawFilePath' to the library on Mac.
libMac :: JVersion -> RawFilePath
libMac v = "/Applications/j64-" <> squashVersionBS v <> "/bin/libj.dylib"
#else
-- | @since 0.1.1.0
libWindows :: JVersion -> FilePath
libWindows v = "C:\\Program Files\\J" <> squashVersion v <> "\\bin\\j.dll"
#endif

profLinux :: BS.ByteString -> BS.ByteString
profLinux v = "/etc/j/" <> v <> "/profile.ijs"

binpathLinux :: BS.ByteString
binpathLinux = "/usr/bin"

dllLinux :: BS.ByteString -> BS.ByteString
dllLinux v = "libj.so." <> v

-- | @since 0.1.2.0
linuxProfile :: BS.ByteString -- ^ J version, e.g. @"9.02"@
             -> Profile
linuxProfile ver = Profile (profLinux ver) binpathLinux (dllLinux ver)

-- | @since 0.1.2.0
macProfile :: JVersion
           -> Profile
macProfile v =
    let binPathMac = "/Applications/j64-" <> squashVersionBS v <> "/bin"
        in Profile (binPathMac <> "/profile.ijs") binPathMac (binPathMac <> "/libj.dylib")

-- | @since 0.1.2.0
windowsProfile :: JVersion
               -> Profile
windowsProfile v =
    let binPathWindows = "C:\\Program Files\\J" <> squashVersionBS v <> "\\bin"
        in Profile (binPathWindows <> "\\profile.ijs") binPathWindows (binPathWindows <> "j.dll")

-- | @since 0.1.2.0
data Profile = Profile { profPath :: BS.ByteString -- ^ @profile.ijs@
                       , binPath  :: BS.ByteString
                       , dllName  :: BS.ByteString
                       }

-- | Load user profile.
--
-- @since 0.1.2.0@
jLoad :: JEnv
      -> Profile
      -> IO ()
jLoad jenv (Profile fp bin dll) = bsDispatch jenv ("(3 : '0!:0 y')<'"<> fp <> "'[BINPATH_z_=:'" <> bin <> "'[LIBFILE_z_=:'" <> dll <> "'[ARGV_z_=:''")

-- | Get a J environment
--
-- Passing the resultant 'JEnv' between threads can cause unexpected bugs.
#ifndef mingw32_HOST_OS
jinit :: RawFilePath -- ^ Path to J library
      -> IO JEnv
jinit libFp = do
    libj <- dlopen libFp [RTLD_LAZY]
    jt <- mkJInit =<< dlsym libj "JInit"
    let jeval = mkJDo <$> dlsym libj "JDo"
    let jread = mkJGetM <$> dlsym libj "JGetM"
    let jOut = mkJGetR <$> dlsym libj "JGetR"
    let jSet = mkJSetA <$> dlsym libj "JSetA"
    JEnv jt <$> jeval <*> jread <*> jOut <*> jSet
#else
jinit :: FilePath
      -> IO JEnv
jinit libFp = do
    libj <- loadLibrary libFp
    jt <- mkJInit . castPtrToFunPtr =<< getProcAddress libj "JInit"
    let jeval = mkJDo . castPtrToFunPtr <$> getProcAddress libj "JDo"
    let jread = mkJGetM . castPtrToFunPtr <$> getProcAddress libj "JGetM"
    let jOut = mkJGetR . castPtrToFunPtr <$> getProcAddress libj "JGetR"
    let jSet = mkJSetA . castPtrToFunPtr <$> getProcAddress libj "JSetA"
    JEnv jt <$> jeval <*> jread <*> jOut <*> jSet
#endif

-- | Send some J code to the environment.
bsDispatch :: JEnv -> BS.ByteString -> IO ()
bsDispatch (JEnv ctx jdo _ _ _) bs =
    void $ BS.useAsCString bs $ jdo ctx

-- | Read last output
--
-- For debugging
bsOut :: JEnv -> IO BS.ByteString
bsOut (JEnv ctx _ _ jout _) = BS.packCString =<< jout ctx

-- | \( O(n) \) in the array size
--
-- Throws 'JErr' for unspported conversion types.
getJData :: R.Shape sh
         => JEnv -> BS.ByteString -- ^ Name of the value in question
         -> IO (JData sh)
getJData jenv bs = jData <$> getAtomInternal jenv bs

getAtomInternal :: JEnv -> BS.ByteString -- ^ Name of the value in question
                -> IO JAtom
getAtomInternal (JEnv ctx _ jget _ _) bs = do
    BS.useAsCString bs $ \name ->
        alloca $ \t ->
        alloca $ \s ->
        alloca $ \r ->
        alloca $ \d -> do
            jget ctx name t r s d
            ty' <- intToJType <$> peek t
            rank' <- peek r
            let intRank = fromIntegral rank'
            shape' <- peekArray intRank =<< peek s
            let mult = jTypeWidth ty'
            let resBytes = mult * intRank
            res <- mallocForeignPtrBytes resBytes
            let arrSz = mult * fromIntegral (product shape')
            withForeignPtr res $ \r' -> do
                d' <- peek d
                copyArray r' d' arrSz
            pure $ JAtom ty' shape' res

data JAtom = JAtom !JType ![CLLong] !(ForeignPtr CChar)

-- | J data backed by repa array
data JData sh = JIntArr !(R.Array RF.F sh CLLong)
              | JDoubleArr !(R.Array RF.F sh CDouble)
              | JComplexArr !(R.Array RF.F sh (Complex CDouble))
              | JBoolArr !(R.Array RF.F sh CChar)
              | JString !BS.ByteString

data JErr = TypeError
          | UnsupportedType
          deriving (Show, Exception)

-- | Copy into a 'V.Vector' 'Int', if possible.
--
-- Throws a 'JErr' on type error.
--
-- \( O(n) \)
--
-- @since 0.2.2.0
tryIntVect :: JData R.DIM1 -> V.Vector Int
tryIntVect (JIntArr arr) = R.toUnboxed (R.copyS $ R.map fromIntegral arr)
tryIntVect _             = throw TypeError

-- | \( O(n) \)
--
-- @since 0.2.2.0
copyIntVect :: V.Vector Int -> JData R.DIM1
copyIntVect = JIntArr . R.copyS . R.map fromIntegral . (\v -> R.fromUnboxed (R.ix1 $ V.length v) v)

-- | \( O(n) \) in the array size
setJData :: (R.Shape sh) => JEnv -> BS.ByteString -- ^ Name
                         -> JData sh -> IO CInt
setJData (JEnv ctx _ _ _ jset) name (JIntArr iarr) = BS.useAsCStringLen name $ \(n, sz) -> do
    (ds, d) <- repaArr JInteger iarr
    jset ctx (fromIntegral sz) n ds d
setJData (JEnv ctx _ _ _ jset) name (JDoubleArr iarr) = BS.useAsCStringLen name $ \(n, sz) -> do
    (ds, d) <- repaArr JDouble iarr
    jset ctx (fromIntegral sz) n ds d
setJData (JEnv ctx _ _ _ jset) name (JComplexArr iarr) = BS.useAsCStringLen name $ \(n, sz) -> do
    (ds, d) <- repaArr JComplex iarr
    jset ctx (fromIntegral sz) n ds d
setJData (JEnv ctx _ _ _ jset) name (JBoolArr iarr) = BS.useAsCStringLen name $ \(n, sz) -> do
    (ds, d) <- repaArr JBool iarr
    jset ctx (fromIntegral sz) n ds d
setJData (JEnv ctx _ _ _ jset) name (JString bs) = BS.useAsCStringLen name $ \(n, sz) -> do
    (ds, d) <- strArr bs
    jset ctx (fromIntegral sz) n ds d

-- | Return a @'Ptr' ()@ suitable to be passed to @JSetA@
--
-- To be used on integer, double, and complex arrays
repaArr :: (R.Shape sh, Storable e) => JType -> R.Array RF.F sh e -> IO (CLLong, Ptr ())
repaArr jty arr = do
    let (rank', sh) = repaSize arr
        sz = product sh
    let wid = 32 + fromIntegral (jTypeWidth jty) * (rank' + sz)
    ptr <- mallocBytes (fromIntegral wid)
    pokeByteOff ptr 0 (227 :: CLLong) -- I think this is because it's non-boxed
    pokeByteOff ptr (sizeOf (undefined :: CLLong)) (jTypeToInt jty)
    pokeByteOff ptr (2 * sizeOf (undefined :: CLLong)) sz
    pokeByteOff ptr (3 * sizeOf (undefined :: CLLong)) rank'
    let dimOff = 4 * sizeOf (undefined :: CLLong)
    pokeArray (ptr `plusPtr` dimOff) sh
    let dataOff = dimOff + fromIntegral rank' * sizeOf (undefined :: CLLong)
    withForeignPtr (RF.toForeignPtr arr) $ \src ->
        copyArray (ptr `plusPtr` dataOff) src (fromIntegral sz)
    pure (wid, ptr)

strArr :: BS.ByteString -> IO (CLLong, Ptr ())
strArr bs = do
    let len = BS.length bs
        wid = 40 + 8 * (1 + len `div` 8)
        len' = fromIntegral len :: CLLong
    ptr <- mallocBytes wid
    pokeByteOff ptr 0 (227 :: CLLong)
    pokeByteOff ptr (sizeOf (undefined :: CLLong)) (jTypeToInt JChar)
    pokeByteOff ptr (2 * sizeOf (undefined :: CLLong)) len'
    pokeByteOff ptr (3 * sizeOf (undefined :: CLLong)) (1 :: CLLong)
    pokeByteOff ptr (4 * sizeOf (undefined :: CLLong)) len'
    let dataOff = 5 * sizeOf (undefined :: CLLong)
    BS.useAsCString bs $ \pSrc ->
        copyArray (ptr `plusPtr` dataOff) pSrc len
    pure (fromIntegral wid, ptr)

repaSize :: (R.Source r e, R.Shape sh) => R.Array r sh e -> (CLLong, [CLLong])
repaSize arr = let sh = R.extent arr in (fromIntegral $ R.rank sh, fromIntegral <$> R.listOfShape sh)

-- | J types
data JType = JBool
           | JChar
           | JInteger
           | JDouble
           | JComplex

intToJType :: CLLong -> JType
intToJType 1  = JBool
intToJType 2  = JChar
intToJType 4  = JInteger
intToJType 8  = JDouble
intToJType 16 = JComplex
intToJType _  = throw UnsupportedType

jTypeToInt :: JType -> CLLong
jTypeToInt JBool    = 1
jTypeToInt JChar    = 2
jTypeToInt JInteger = 4
jTypeToInt JDouble  = 8
jTypeToInt JComplex = 16

jTypeWidth :: JType -> Int
jTypeWidth JBool    = sizeOf (undefined :: CChar)
jTypeWidth JChar    = sizeOf (undefined :: CChar)
jTypeWidth JInteger = sizeOf (undefined :: CLLong)
jTypeWidth JDouble  = sizeOf (undefined :: CDouble)
jTypeWidth JComplex = sizeOf (undefined :: Complex CDouble)

jData :: R.Shape sh => JAtom -> JData sh
jData (JAtom JInteger sh fp) = JIntArr $ RF.fromForeignPtr (R.shapeOfList $ fmap fromIntegral sh) (castForeignPtr fp)
jData (JAtom JDouble sh fp)  = JDoubleArr $ RF.fromForeignPtr (R.shapeOfList $ fmap fromIntegral sh) (castForeignPtr fp)
jData (JAtom JComplex sh fp) = JComplexArr $ RF.fromForeignPtr (R.shapeOfList $ fmap fromIntegral sh) (castForeignPtr fp)
jData (JAtom JBool sh fp)    = JBoolArr $ RF.fromForeignPtr (R.shapeOfList $ fmap fromIntegral sh) (castForeignPtr fp)
jData (JAtom JChar [l] fp)   = JString $ BS.fromForeignPtr (castForeignPtr fp) 0 (fromIntegral l)
jData (JAtom JChar _ _)      = error "Not supported."

max_stars_repo_name: vmchale/j-hs

max_issues_repo_name: vmchale/j-hs

max_forks_repo_name: vmchale/j-hs


text:
{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE CPP                 #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE Strict              #-}
{-# LANGUAGE TemplateHaskell     #-}
{-# LANGUAGE TypeOperators       #-}

module Grenade.Layers.Internal.CUDA
  ( setCudaTriggerSize
  , useCuda
  , cudaDescendUnsafeAdamGPU
  , cudaSumVectorsGPU
  ) where

import           Control.Concurrent
import           Control.Concurrent.Lock      (Lock)
import qualified Control.Concurrent.Lock      as Lock
import           Control.Concurrent.MVar
import           Control.Exception            (IOException, handle)
import           Control.Monad
import qualified Data.Map.Strict              as M
import           Data.Maybe                   (fromMaybe)
import           Data.Proxy
import           Data.String
import qualified Data.Vector.Storable         as V
import qualified Data.Vector.Storable         as U (unsafeFromForeignPtr0,
                                                    unsafeToForeignPtr0)
import           Foreign                      (withForeignPtr)
import           Foreign.C.Types
-- import qualified Foreign.CUDA.BLAS            as BLAS
import           Control.Monad.Primitive      (PrimBase)
import qualified Data.ByteString
import           Data.FileEmbed               (embedFile)
import qualified Foreign.CUDA.Driver          as CUDA
import           Foreign.Ptr
import           Foreign.Storable             (sizeOf)
import           GHC.IO.Handle.Text           (memcpy)
import           GHC.TypeLits
-- import qualified Numeric.LinearAlgebra        as LA
import qualified Numeric.LinearAlgebra.Static as LAS
import           System.IO.Unsafe             (unsafePerformIO)

import           Grenade.Types
import           Grenade.Utils.Vector

import           Debug.Trace

-- | Cuda context
cudaTriggerSize :: MVar (Maybe Int)
cudaTriggerSize = unsafePerformIO $ newMVar Nothing

-- | Set the minimum vector/matrix size (nr of elements) of when CUDA is used.
setCudaTriggerSize :: (PrimBase m) => Maybe Int -> m ()
setCudaTriggerSize sz = unsafePerformIO (tryPutMVar cudaTriggerSize sz) `seq` return ()

useCuda :: Int -> Bool
useCuda sz = unsafePerformIO $ do
  minSz <- tryReadMVar cudaTriggerSize
  return $ maybe False (sz >=) (join minSz)

-- Load the ptx file into the binary at compile-time!
bsGPUGradientDescent :: Data.ByteString.ByteString
bsGPUGradientDescent = $(embedFile "cbits/gpu_gradient_descent.ptx")

-- Load the ptx file into the binary at compile-time!
bsGPUSumVectors :: Data.ByteString.ByteString
bsGPUSumVectors = $(embedFile "cbits/gpu_sum_vectors.ptx")

-- | Cuda function for Adam
cudaDescendAdamGPUFun :: MVar (Maybe CUDA.Fun)
cudaDescendAdamGPUFun = unsafePerformIO $ newMVar Nothing
{-# NOINLINE cudaDescendAdamGPUFun #-}

-- | Cuda function for sumVectors
cudaSumVectorsGPUFun :: MVar (Maybe CUDA.Fun)
cudaSumVectorsGPUFun = unsafePerformIO $ newMVar Nothing
{-# NOINLINE cudaSumVectorsGPUFun #-}

-- | Cuda function for sumVectors
lock :: Lock
lock = unsafePerformIO Lock.new
{-# NOINLINE lock #-}

writeMVar :: MVar a -> a -> IO ()
writeMVar var x = modifyMVar var (\_ -> return (x, ()))

-- cudaInitialise :: IO ()
-- cudaInitialise = do
--   CUDA.initialise []
--   count <- CUDA.count
--   devs <- mapM CUDA.device [0 .. count - 1]
--   putStrLn $ "Initialising CUDA. Found " ++ show count ++ " devices: "
--   mapM_ (CUDA.props >=> print) devs
--   putStrLn ""
--   unless (null devs) $ do
--     props <- CUDA.props (head devs)
--     let threadsPerBlock = CUDA.maxThreadsPerBlock props
--     putStrLn $ "Using CUDA with device: " ++ show (CUDA.deviceName props) ++ "; threads per Block: " ++ show threadsPerBlock ++ "\n"
--     writeIORef cudaThreadsPerBlock threadsPerBlock
--     device <- CUDA.device 0
--     writeIORef cudaDevice (Just device)
--     ctx <- CUDA.create device []
--     writeIORef cudaCtx (Just ctx)
--     -- Compile to ptx?
--     mdl <- CUDA.loadData bsGPUGradientDescent
--     descendAdamGpuFun <- CUDA.getFun mdl (fromString "descend_adam_gpu")
--     writeIORef cudaDescendAdamGPUFun (Just descendAdamGpuFun)

cudaInitialise :: IO ()
cudaInitialise = do
  threadId <- myThreadId
  Lock.acquire lock
  -- m <- readIORef cudaCtx
  -- case M.lookup threadId m of
  --   Just{} -> return ()
  --   Nothing -> do
  do
      CUDA.initialise []
      count <- CUDA.count
      devs <- mapM CUDA.device [0 .. count - 1]
      if null devs
        then return ()
        else do
          putStrLn $ "Initialising CUDA. Found " ++ show count ++ " devices: "
          mapM_ (CUDA.props >=> print) devs
          putStrLn ""
          props <- CUDA.props (head devs)
          let threadsPerBlock = CUDA.maxThreadsPerBlock props
          putStrLn $ "Using CUDA with device: " ++ show (CUDA.deviceName props) ++ "; threads per Block: " ++ show threadsPerBlock ++ "\n"
          writeMVar cudaThreadsPerBlock threadsPerBlock
          device <- CUDA.device 0
          writeMVar cudaDevice (Just device)
          -- ctx <- CUDA.create device []
          -- -- atomicModifyIORef cudaCtx (\m' -> (M.insert threadId ctx m', ()))
          -- -- writeIORef cudaCtx (M.Just ctx)
          --      -- Inline compile to ptx?
          -- flip mapM_ functions $ \(file, functionName, ioRef) -> do
          --   mdl <- CUDA.loadData file
          --   fun <- CUDA.getFun mdl (fromString functionName)
          --   modifyMVar ioRef (\_ -> return (Just fun, ()))
  Lock.release lock
  where
    functions =
      [ (bsGPUGradientDescent, "descend_adam_gpu", cudaDescendAdamGPUFun) -- ^ Adam
      , (bsGPUSumVectors, "sum_vectors_gpu", cudaSumVectorsGPUFun) -- ^ Sum a list of vectors
      ]

-- | Cuda context
cudaThreadsPerBlock :: MVar Int
cudaThreadsPerBlock = unsafePerformIO $ newMVar 128

-- | Cuda device
cudaDevice :: MVar (Maybe CUDA.Device)
cudaDevice = unsafePerformIO $ newMVar Nothing

withCudaCtx :: CUDA.Device -> (CUDA.Fun -> IO b) -> IO b
withCudaCtx device f = do
  -- device <- fromMaybe (error "empty device") . join <$> tryReadMVar cudaDevice -- device 0
  -- myThreadId >>= print
  !ctx <- CUDA.create device []
  mdl <- CUDA.loadData bsGPUSumVectors
  fun <- CUDA.getFun mdl (fromString "sum_vectors_gpu")
  !b <- f fun
  CUDA.destroy ctx
  return b

-- | To inform the user when it works again
cudaErrored :: MVar Bool
cudaErrored = unsafePerformIO $ newMVar False
{-# NOINLINE cudaErrored #-}

cudaDescendUnsafeAdamGPU ::
     Int -- Len
  -> Int -- Step
  -> RealNum -- Alpha
  -> RealNum -- Beta1
  -> RealNum -- Beta2
  -> RealNum -- Epsilon
  -> RealNum -- Lambda
  -> V.Vector RealNum -- Weights
  -> V.Vector RealNum -- Gradient
  -> V.Vector RealNum -- M
  -> V.Vector RealNum -- V
  -> Maybe (V.Vector RealNum, V.Vector RealNum, V.Vector RealNum)
cudaDescendUnsafeAdamGPU len step alpha beta1 beta2 epsilon lambda weights gradient m v =
  unsafePerformIO $ do
    mFun <- join <$> tryReadMVar cudaDescendAdamGPUFun
    case mFun of
      Nothing -> cudaInitialise >> tryReadMVar cudaDescendAdamGPUFun >>= maybe (return Nothing) cudaDescendAdamGPU' -- Ensure it was initialised
      x -> handle handler (cudaDescendAdamGPU' x)
  where
    handler (e :: CUDA.CUDAException) = writeMVar cudaErrored True >> print e >> return Nothing
    cudaDescendAdamGPU' mFun =
      case mFun of
        Nothing -> return Nothing
        Just fun -> do
          threadsPerBlock <- fromMaybe 128 <$> tryReadMVar cudaThreadsPerBlock
          V.unsafeWith weights $ \wPtr' ->
            V.unsafeWith gradient $ \gPtr' ->
              V.unsafeWith m $ \mPtr' ->
                V.unsafeWith v $ \vPtr' -> do
                  wPtrDev <- CUDA.mallocArray len :: IO (CUDA.DevicePtr RealNum)
                  gPtrDev <- CUDA.mallocArray len :: IO (CUDA.DevicePtr RealNum)
                  mPtrDev <- CUDA.mallocArray len :: IO (CUDA.DevicePtr RealNum)
                  vPtrDev <- CUDA.mallocArray len :: IO (CUDA.DevicePtr RealNum)
                  CUDA.pokeArray len wPtr' wPtrDev
                  CUDA.pokeArray len gPtr' gPtrDev
                  CUDA.pokeArray len mPtr' mPtrDev
                  CUDA.pokeArray len vPtr' vPtrDev
                  let tx = ceiling (fromIntegral len / (fromIntegral threadsPerBlock :: Float)) :: Int
                  -- print (len, step+1)
                  -- print (alpha, beta1, beta2, epsilon, lambda)
                  CUDA.launchKernel
                    fun
                    (tx, 1, 1 :: Int)
                    (threadsPerBlock, 1, 1)
                    0
                    Nothing
                    [ CUDA.IArg (fromIntegral len)
                    , CUDA.IArg (fromIntegral $ step + 1)
                    , CUDA.VArg (alpha :: RealNum)
                    , CUDA.VArg (beta1 :: RealNum)
                    , CUDA.VArg (beta2 :: RealNum)
                    , CUDA.VArg (epsilon :: RealNum)
                    , CUDA.VArg (lambda :: RealNum)
                    , CUDA.VArg wPtrDev
                    , CUDA.VArg gPtrDev
                    , CUDA.VArg mPtrDev
                    , CUDA.VArg vPtrDev
                    ]
                  CUDA.sync
                  CUDA.peekArray len wPtrDev wPtr'
                  CUDA.peekArray len mPtrDev mPtr'
                  CUDA.peekArray len vPtrDev vPtr'
                  mapM_ CUDA.free [wPtrDev,gPtrDev,mPtrDev,vPtrDev] -- put pointers back for reuse
                  err <- fromMaybe False <$> tryReadMVar cudaErrored
                  when err $ writeMVar cudaErrored False >> putStrLn "Cuda works again"
                  return $ Just (weights, m, v)
{-# NOINLINE cudaDescendUnsafeAdamGPU #-}

cudaSumVectorsGPU :: [V.Vector RealNum] -> Maybe (V.Vector RealNum) -- vectors to sum, all of same length
cudaSumVectorsGPU [] = Nothing
cudaSumVectorsGPU vs
  | useCuda (length vs * V.length (head vs)) =
    unsafePerformIO $! do
      mDev <- join <$> tryReadMVar cudaDevice
      case mDev of
        Nothing -> cudaInitialise >> tryReadMVar cudaDevice >>= maybe (return Nothing) cudaSumVectorsGPU' -- Ensure it was initialised
        dev -> handle handler (cudaSumVectorsGPU' dev)
  | otherwise = Nothing
  where
    handler (e :: CUDA.CUDAException) = writeMVar cudaErrored True >> print e >> return Nothing
    cudaSumVectorsGPU' (mDev :: Maybe CUDA.Device) =
      case mDev of
        Nothing -> return Nothing
        Just dev ->
          withCudaCtx dev $ \fun -> do
            let outVec = createVectorUnsafe (V.length $ head vs)
            let outLen = V.length (head vs)
                inLen = outLen * length vs
            let inVec = V.generate inLen genFun
                genFun idx =
                  let vIdx = idx `div` outLen
                      vec = vs !! vIdx
                   in vec V.! (idx - vIdx * outLen)
            threadsPerBlock <- fromMaybe 128 <$> tryReadMVar cudaThreadsPerBlock
            V.unsafeWith inVec $ \inVecPtr' -> do
              V.unsafeWith outVec $ \outVecPtr' -> do
                inVecPtrDev <- CUDA.mallocArray inLen :: IO (CUDA.DevicePtr RealNum)
                outVecPtrDev <- CUDA.mallocArray outLen :: IO (CUDA.DevicePtr RealNum)
                CUDA.pokeArray inLen inVecPtr' inVecPtrDev
                let tx = ceiling (fromIntegral outLen / (fromIntegral threadsPerBlock :: Float)) :: Int
                CUDA.launchKernel
                  fun
                  (tx, 1, 1 :: Int)
                  (threadsPerBlock, 1, 1)
                  0
                  Nothing
                  [ CUDA.IArg (fromIntegral inLen)
                  , CUDA.IArg (fromIntegral outLen)
                  , CUDA.VArg (inVecPtrDev :: CUDA.DevicePtr RealNum)
                  , CUDA.VArg (outVecPtrDev :: CUDA.DevicePtr RealNum)
                  ]
                CUDA.sync
                CUDA.peekArray outLen outVecPtrDev outVecPtr'
                mapM_ CUDA.free [inVecPtrDev, outVecPtrDev] -- put pointers back for reuse
                err <- fromMaybe False <$> tryReadMVar cudaErrored
                when err $ writeMVar cudaErrored False >> putStrLn "Cuda works again"
                return $ Just outVec
{-# NOINLINE cudaSumVectorsGPU #-}

-- test =
--   cudaSumVectorsGPU (map V.fromList (replicate 32 [1..10]))

max_stars_repo_name: schnecki/grenade

max_issues_repo_name: schnecki/grenade

max_forks_repo_name: schnecki/grenade