text:
{-# language GADTs                 #-}
{-# language StandaloneDeriving    #-}
{-# language TypeOperators         #-}
{-# language FlexibleInstances     #-}
{-# language FlexibleContexts      #-}
{-# language UndecidableInstances  #-}
{-# language MultiParamTypeClasses #-}
{-# language TypeFamilies          #-}

{-# options_ghc -fwarn-incomplete-patterns #-}

module Feldspar.Software.Primitive where

import Feldspar.Representation
import Data.Struct

import Data.Array ((!))
import Data.Bits (Bits)
import Data.Complex
import Data.Int
import Data.Word
import Data.List (genericTake)
import Data.Typeable hiding (TypeRep)
import Data.Constraint hiding (Sub)
import qualified Data.Bits as Bits

-- syntactic.
import Language.Syntactic
import Language.Syntactic.Functional
import Language.Syntactic.Functional.Tuple
import qualified Language.Syntactic as Syn

-- imperative-edsl.
import Language.Embedded.Expression
import qualified Language.Embedded.Imperative.CMD as Imp (IArr(..))

--------------------------------------------------------------------------------
-- * Software primitives.
--------------------------------------------------------------------------------

--------------------------------------------------------------------------------
-- ** Software primitive types.

-- | Representation of supported, primitive software types.
data SoftwarePrimTypeRep a
  where
    -- booleans
    BoolST   :: SoftwarePrimTypeRep Bool
    -- signed numbers.
    Int8ST   :: SoftwarePrimTypeRep Int8
    Int16ST  :: SoftwarePrimTypeRep Int16
    Int32ST  :: SoftwarePrimTypeRep Int32
    Int64ST  :: SoftwarePrimTypeRep Int64
    -- unsigned numbers.
    Word8ST  :: SoftwarePrimTypeRep Word8
    Word16ST :: SoftwarePrimTypeRep Word16
    Word32ST :: SoftwarePrimTypeRep Word32
    Word64ST :: SoftwarePrimTypeRep Word64
    -- floating point numbers.
    FloatST  :: SoftwarePrimTypeRep Float
    DoubleST :: SoftwarePrimTypeRep Double
    -- complex numbers.
    ComplexFloatST  :: SoftwarePrimTypeRep (Complex Float)
    ComplexDoubleST :: SoftwarePrimTypeRep (Complex Double)

deriving instance Eq       (SoftwarePrimTypeRep a)
deriving instance Show     (SoftwarePrimTypeRep a)
deriving instance Typeable (SoftwarePrimTypeRep a)

instance Inhabited (Complex Float)
  where
    reset = 0 :+ 0

instance Inhabited (Complex Double)
  where
    reset = 0 :+ 0

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

-- | Class of supported, primitive software types.
class (Eq a, Show a, Typeable a, Inhabited a) => SoftwarePrimType a
  where
    softwareRep :: SoftwarePrimTypeRep a

instance SoftwarePrimType Bool   where softwareRep = BoolST
instance SoftwarePrimType Int8   where softwareRep = Int8ST
instance SoftwarePrimType Int16  where softwareRep = Int16ST
instance SoftwarePrimType Int32  where softwareRep = Int32ST
instance SoftwarePrimType Int64  where softwareRep = Int64ST
instance SoftwarePrimType Word8  where softwareRep = Word8ST
instance SoftwarePrimType Word16 where softwareRep = Word16ST
instance SoftwarePrimType Word32 where softwareRep = Word32ST
instance SoftwarePrimType Word64 where softwareRep = Word64ST
instance SoftwarePrimType Float  where softwareRep = FloatST
instance SoftwarePrimType Double where softwareRep = DoubleST
instance SoftwarePrimType (Complex Float)  where softwareRep = ComplexFloatST
instance SoftwarePrimType (Complex Double) where softwareRep = ComplexDoubleST

-- | Compare two primitive software types for equality.
softwarePrimTypeEq :: SoftwarePrimTypeRep a -> SoftwarePrimTypeRep b -> Maybe (Dict (a ~ b))
softwarePrimTypeEq (BoolST)   (BoolST)   = Just Dict
softwarePrimTypeEq (Int8ST)   (Int8ST)   = Just Dict
softwarePrimTypeEq (Int16ST)  (Int16ST)  = Just Dict
softwarePrimTypeEq (Int32ST)  (Int32ST)  = Just Dict
softwarePrimTypeEq (Int64ST)  (Int64ST)  = Just Dict
softwarePrimTypeEq (Word8ST)  (Word8ST)  = Just Dict
softwarePrimTypeEq (Word16ST) (Word16ST) = Just Dict
softwarePrimTypeEq (Word32ST) (Word32ST) = Just Dict
softwarePrimTypeEq (Word64ST) (Word64ST) = Just Dict
softwarePrimTypeEq (FloatST)  (FloatST)  = Just Dict
softwarePrimTypeEq (DoubleST) (DoubleST) = Just Dict
softwarePrimTypeEq (ComplexFloatST)  (ComplexFloatST)  = Just Dict
softwarePrimTypeEq (ComplexDoubleST) (ComplexDoubleST) = Just Dict
softwarePrimTypeEq _          _          = Nothing

-- | Construct the primitive software type representation of 'a'.
softwarePrimTypeOf :: SoftwarePrimType a => a -> SoftwarePrimTypeRep a
softwarePrimTypeOf _ = softwareRep

-- | Construct a primitive software type witness from its representation.
softwarePrimWitType :: SoftwarePrimTypeRep a -> Dict (SoftwarePrimType a)
softwarePrimWitType BoolST   = Dict
softwarePrimWitType Int8ST   = Dict
softwarePrimWitType Int16ST  = Dict
softwarePrimWitType Int32ST  = Dict
softwarePrimWitType Int64ST  = Dict
softwarePrimWitType Word8ST  = Dict
softwarePrimWitType Word16ST = Dict
softwarePrimWitType Word32ST = Dict
softwarePrimWitType Word64ST = Dict
softwarePrimWitType FloatST  = Dict
softwarePrimWitType DoubleST = Dict
softwarePrimWitType ComplexFloatST  = Dict
softwarePrimWitType ComplexDoubleST = Dict

--------------------------------------------------------------------------------
-- ** Software primitive expressions.

-- | Software primitive symbols.
data SoftwarePrim sig
  where
    -- free variables and literals.
    FreeVar :: (SoftwarePrimType a) => String -> SoftwarePrim (Full a)
    Lit     :: (Show a, Eq a)       => a      -> SoftwarePrim (Full a)
    -- numerical operations.
    Neg  :: (SoftwarePrimType a, Num a) => SoftwarePrim (a :-> Full a)
    Abs  :: (SoftwarePrimType a, Num a) => SoftwarePrim (a :-> Full a)
    Sign :: (SoftwarePrimType a, Num a) => SoftwarePrim (a :-> Full a)
    Add  :: (SoftwarePrimType a, Num a) => SoftwarePrim (a :-> a :-> Full a)
    Sub  :: (SoftwarePrimType a, Num a) => SoftwarePrim (a :-> a :-> Full a)
    Mul  :: (SoftwarePrimType a, Num a) => SoftwarePrim (a :-> a :-> Full a)
    -- integral operations.
    Div  :: (SoftwarePrimType a, Integral a) => SoftwarePrim (a :-> a :-> Full a)
    Mod  :: (SoftwarePrimType a, Integral a) => SoftwarePrim (a :-> a :-> Full a)
    Quot :: (SoftwarePrimType a, Integral a) => SoftwarePrim (a :-> a :-> Full a)
    Rem  :: (SoftwarePrimType a, Integral a) => SoftwarePrim (a :-> a :-> Full a)
    --
    FDiv :: (SoftwarePrimType a, Fractional a) => SoftwarePrim (a :-> a :-> Full a)
    -- floating point operators.
    Pi    :: (SoftwarePrimType a, Floating a) => SoftwarePrim (Full a)
    Exp   :: (SoftwarePrimType a, Floating a) => SoftwarePrim (a :-> Full a)
    Log   :: (SoftwarePrimType a, Floating a) => SoftwarePrim (a :-> Full a)
    Sqrt  :: (SoftwarePrimType a, Floating a) => SoftwarePrim (a :-> Full a)
    Pow   :: (SoftwarePrimType a, Floating a) => SoftwarePrim (a :-> a :-> Full a)
    Sin   :: (SoftwarePrimType a, Floating a) => SoftwarePrim (a :-> Full a)
    Cos   :: (SoftwarePrimType a, Floating a) => SoftwarePrim (a :-> Full a)
    Tan   :: (SoftwarePrimType a, Floating a) => SoftwarePrim (a :-> Full a)
    Asin  :: (SoftwarePrimType a, Floating a) => SoftwarePrim (a :-> Full a)
    Acos  :: (SoftwarePrimType a, Floating a) => SoftwarePrim (a :-> Full a)
    Atan  :: (SoftwarePrimType a, Floating a) => SoftwarePrim (a :-> Full a)
    Sinh  :: (SoftwarePrimType a, Floating a) => SoftwarePrim (a :-> Full a)
    Cosh  :: (SoftwarePrimType a, Floating a) => SoftwarePrim (a :-> Full a)
    Tanh  :: (SoftwarePrimType a, Floating a) => SoftwarePrim (a :-> Full a)
    Asinh :: (SoftwarePrimType a, Floating a) => SoftwarePrim (a :-> Full a)
    Acosh :: (SoftwarePrimType a, Floating a) => SoftwarePrim (a :-> Full a)
    Atanh :: (SoftwarePrimType a, Floating a) => SoftwarePrim (a :-> Full a)
    -- complex operators.
    Complex   :: (SoftwarePrimType a, SoftwarePrimType (Complex a), Num a) =>
      SoftwarePrim (a :-> a :-> Full (Complex a))
    Real      :: (SoftwarePrimType a, SoftwarePrimType (Complex a)) =>
      SoftwarePrim (Complex a :-> Full a)
    Imag      :: (SoftwarePrimType a, SoftwarePrimType (Complex a)) =>
      SoftwarePrim (Complex a :-> Full a)
    Polar     :: (SoftwarePrimType a, SoftwarePrimType (Complex a), Floating a) =>
      SoftwarePrim (a :-> a :-> Full (Complex a))
    Magnitude :: (SoftwarePrimType a, SoftwarePrimType (Complex a), RealFloat a) =>
      SoftwarePrim (Complex a :-> Full a)
    Phase     :: (SoftwarePrimType a, SoftwarePrimType (Complex a), RealFloat a) =>
      SoftwarePrim (Complex a :-> Full a)
    Conjugate :: (SoftwarePrimType a, SoftwarePrimType (Complex a), Num a) =>
      SoftwarePrim (Complex a :-> Full (Complex a))
    -- type casting.
    I2N   :: (SoftwarePrimType a, Integral a, SoftwarePrimType b, Num b) =>
      SoftwarePrim (a :-> Full b)
    I2B   :: (SoftwarePrimType a, Integral a) =>
      SoftwarePrim (a :-> Full Bool)
    B2I   :: (SoftwarePrimType a, Integral a) =>
      SoftwarePrim (Bool :-> Full a)
    Round :: (SoftwarePrimType a, RealFrac a, SoftwarePrimType b, Num b) =>
      SoftwarePrim (a :-> Full b)
    -- logical operations.
    Not     :: SoftwarePrim (Bool :-> Full Bool)
    And     :: SoftwarePrim (Bool :-> Bool :-> Full Bool)
    Or      :: SoftwarePrim (Bool :-> Bool :-> Full Bool)
    -- bitwise logical operations.
    BitAnd   :: (SoftwarePrimType a, Bits a) => SoftwarePrim (a :-> a :-> Full a)
    BitOr    :: (SoftwarePrimType a, Bits a) => SoftwarePrim (a :-> a :-> Full a)
    BitXor   :: (SoftwarePrimType a, Bits a) => SoftwarePrim (a :-> a :-> Full a)
    BitCompl :: (SoftwarePrimType a, Bits a) => SoftwarePrim (a :-> Full a)
    ShiftL   :: (SoftwarePrimType a, Bits a, SoftwarePrimType b, Integral b) =>
      SoftwarePrim (a :-> b :-> Full a)
    ShiftR   :: (SoftwarePrimType a, Bits a, SoftwarePrimType b, Integral b) =>
      SoftwarePrim (a :-> b :-> Full a)
    RotateL  :: (SoftwarePrimType a, Bits a, SoftwarePrimType b, Integral b) =>
      SoftwarePrim (a :-> b :-> Full a)
    RotateR  :: (SoftwarePrimType a, Bits a, SoftwarePrimType b, Integral b) =>
      SoftwarePrim (a :-> b :-> Full a)
    -- relational operations.
    Eq  :: (SoftwarePrimType a, Eq a)  => SoftwarePrim (a :-> a :-> Full Bool)
    Neq :: (SoftwarePrimType a, Eq a)  => SoftwarePrim (a :-> a :-> Full Bool)
    Lt  :: (SoftwarePrimType a, Ord a) => SoftwarePrim (a :-> a :-> Full Bool)
    Lte :: (SoftwarePrimType a, Ord a) => SoftwarePrim (a :-> a :-> Full Bool)
    Gt  :: (SoftwarePrimType a, Ord a) => SoftwarePrim (a :-> a :-> Full Bool)
    Gte :: (SoftwarePrimType a, Ord a) => SoftwarePrim (a :-> a :-> Full Bool)
    -- conditional.
    Cond :: SoftwarePrim (Bool :-> a :-> a :-> Full a)
    -- array indexing.
    ArrIx :: (SoftwarePrimType a) => Imp.IArr Index a ->
      SoftwarePrim (Index :-> Full a)

deriving instance Show     (SoftwarePrim a)
deriving instance Typeable (SoftwarePrim a)

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

-- | Software primitive symbols.
type SoftwarePrimConstructs = SoftwarePrim

-- | Software primitive symbols tagged with their type representation.
type SoftwarePrimDomain = SoftwarePrimConstructs :&: SoftwarePrimTypeRep

-- | Software primitive expressions.
newtype Prim a = Prim { unPrim :: ASTF SoftwarePrimDomain a }

-- | Evaluate a closed, software primitive expression.
evalPrim :: Prim a -> a
evalPrim = go . unPrim
  where
    go :: AST SoftwarePrimDomain sig -> Denotation sig
    go (Sym (s :&: _)) = evalSym s
    go (f :$ a)        = go f $ go a

-- | Sugar a software primitive symbol as a smart constructor.
sugarSymPrim
  :: ( Signature sig
     , fi  ~ SmartFun dom sig
     , sig ~ SmartSig fi
     , dom ~ SmartSym fi
     , dom ~ SoftwarePrimDomain
     , SyntacticN f fi
     , sub :<: SoftwarePrimConstructs
     , SoftwarePrimType (DenResult sig)
     )
  => sub sig -> f
sugarSymPrim = sugarSymDecor softwareRep

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

instance Syntactic (Prim a)
  where
    type Domain   (Prim a) = SoftwarePrimDomain
    type Internal (Prim a) = a
    desugar = unPrim
    sugar   = Prim

instance FreeExp Prim
  where
    type FreePred Prim = SoftwarePrimType
    constExp = sugarSymPrim . Lit
    varExp   = sugarSymPrim . FreeVar

instance EvalExp Prim
  where
    evalExp = evalPrim

--------------------------------------------------------------------------------
-- front-end.

instance (SoftwarePrimType a, Num a) => Num (Prim a)
  where
    fromInteger = constExp . fromInteger
    (+)         = sugarSymPrim Add
    (-)         = sugarSymPrim Sub
    (*)         = sugarSymPrim Mul
    negate      = sugarSymPrim Neg
    abs         = error "Num (Prim a): abs."
    signum      = error "Num (Prim a): signum."

--------------------------------------------------------------------------------
-- syntactic instances.

instance Eval SoftwarePrim
  where
    evalSym (FreeVar v) = error $ "evaluating free variable " ++ show v
    evalSym (Lit a)     = a
    evalSym Cond        = \c t f -> if c then t else f
    evalSym Neg         = negate
    evalSym Abs         = abs
    evalSym Sign        = signum
    evalSym Add         = (+)
    evalSym Sub         = (-)
    evalSym Mul         = (*)
    evalSym Div         = div
    evalSym Mod         = mod
    evalSym Quot        = quot
    evalSym Rem         = rem
    evalSym FDiv        = (/)
    evalSym Pi          = pi
    evalSym Exp         = exp
    evalSym Log         = log
    evalSym Sqrt        = sqrt
    evalSym Pow         = (**)
    evalSym Sin         = sin
    evalSym Cos         = cos
    evalSym Tan         = tan
    evalSym Asin        = asin
    evalSym Acos        = acos
    evalSym Atan        = atan
    evalSym Sinh        = sinh
    evalSym Cosh        = cosh
    evalSym Tanh        = tanh
    evalSym Asinh       = asinh
    evalSym Acosh       = acosh
    evalSym Atanh       = atanh
    evalSym Complex     = (:+)
    evalSym Polar       = mkPolar
    evalSym Real        = realPart
    evalSym Imag        = imagPart
    evalSym Magnitude   = magnitude
    evalSym Phase       = phase
    evalSym Conjugate   = conjugate
    evalSym I2N         = fromIntegral
    evalSym I2B         = (/=0)
    evalSym B2I         = \a -> if a then 1 else 0
    evalSym Round       = fromInteger . round
    evalSym Not         = not
    evalSym And         = (&&)
    evalSym Or          = (||)
    evalSym BitAnd      = (Bits..&.)
    evalSym BitOr       = (Bits..|.)
    evalSym BitXor      = Bits.xor
    evalSym BitCompl    = Bits.complement
    evalSym ShiftL      = \b i -> Bits.shiftL  b (fromIntegral i)
    evalSym ShiftR      = \b i -> Bits.shiftR  b (fromIntegral i)
    evalSym RotateL     = \b i -> Bits.rotateL b (fromIntegral i)
    evalSym RotateR     = \b i -> Bits.rotateR b (fromIntegral i)
    evalSym Eq          = (==)
    evalSym Neq         = (/=)
    evalSym Lt          = (<)
    evalSym Lte         = (<=)
    evalSym Gt          = (>)
    evalSym Gte         = (>=)
    evalSym (ArrIx (Imp.IArrRun arr)) = \i -> arr ! i
    evalSym (ArrIx _)   = error "eval of variable."

instance Symbol SoftwarePrim
  where
    symSig (FreeVar v) = signature
    symSig (Lit a)     = signature
    symSig Cond        = signature
    symSig Neg         = signature
    symSig Abs         = signature
    symSig Sign        = signature
    symSig Add         = signature
    symSig Sub         = signature
    symSig Mul         = signature
    symSig Div         = signature
    symSig Mod         = signature
    symSig Quot        = signature
    symSig Rem         = signature
    symSig FDiv        = signature
    symSig Pi          = signature
    symSig Exp         = signature
    symSig Log         = signature
    symSig Sqrt        = signature
    symSig Pow         = signature
    symSig Sin         = signature
    symSig Cos         = signature
    symSig Tan         = signature
    symSig Asin        = signature
    symSig Acos        = signature
    symSig Atan        = signature
    symSig Sinh        = signature
    symSig Cosh        = signature
    symSig Tanh        = signature
    symSig Asinh       = signature
    symSig Acosh       = signature
    symSig Atanh       = signature
    symSig Complex     = signature
    symSig Real        = signature
    symSig Imag        = signature
    symSig Polar       = signature
    symSig Magnitude   = signature
    symSig Phase       = signature
    symSig Conjugate   = signature
    symSig I2N         = signature
    symSig I2B         = signature
    symSig B2I         = signature
    symSig Round       = signature
    symSig Not         = signature
    symSig And         = signature
    symSig Or          = signature
    symSig BitAnd      = signature
    symSig BitOr       = signature
    symSig BitXor      = signature
    symSig BitCompl    = signature
    symSig ShiftL      = signature
    symSig ShiftR      = signature
    symSig RotateL     = signature
    symSig RotateR     = signature
    symSig Eq          = signature
    symSig Neq         = signature
    symSig Lt          = signature
    symSig Lte         = signature
    symSig Gt          = signature
    symSig Gte         = signature
    symSig (ArrIx a)   = signature

instance Render SoftwarePrim
  where
    renderSym  = show
    renderArgs = renderArgsSmart

instance Equality SoftwarePrim
  where
    equal s1 s2 = show s1 == show s2

instance StringTree SoftwarePrim
instance EvalEnv SoftwarePrim env

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

max_stars_repo_name: markus-git/co-feldspar

max_issues_repo_name: markus-git/co-feldspar

max_forks_repo_name: markus-git/co-feldspar


text:
module GPLVM.PPCA
       ( PPCA (..)
       , makePPCA
       ) where

import Prelude (isNaN, log)

import Universum hiding (All, Any, Vector, map, toList, transpose, (++))

import qualified Universum as U

import GPLVM.Types
import GPLVM.Util

import Data.Array.Repa
import Data.Array.Repa.Repr.ForeignPtr
import Data.List ((\\))
import Numeric.LinearAlgebra.Repa hiding (Matrix, Vector, diag)
import System.Random

data PPCA = PPCA
  {  _noMissedData        :: Bool
   , _learningData        :: Matrix D Double
   , desiredDimension     :: Int
   , stopParameter        :: Either Int Double
   , _variance            :: Double
   , _W                   :: Matrix D Double
   , _finalExpLikelihood  :: Double
   , _restoredMatrix      :: Maybe (Matrix D Double)
   }

makePPCA
  :: (HasCallStack) => RandomGen gen
  => Matrix D Double   -- ^ learning vectors
  -> Int               -- ^ desired dimention
  -> Either Int Double -- ^ number of iterations or
                       --   the value of stop difference of the likelihood expectation between interations of EM.
  -> gen
  -> PPCA
makePPCA leaningVectors desiredDimension stopParameter generator =
  let _learningData@(ADelayed (Z :. n :. _) _) = leaningVectors
      initMatrix = randomMatrixD generator (n, desiredDimension)
      initDispersion = fromIntegral $ fst . next $ generator
      _noMissedData = not $ isNaN $ runIdentity $ sumAllP _learningData
      (_W, _variance, _finalExpLikelihood, _restoredMatrix) =
        if _noMissedData
        then emStepsFast _learningData initMatrix initDispersion stopParameter
        else emStepsMissed _learningData initMatrix initDispersion stopParameter
  in PPCA{..}

-- TODO: add comments to the most of following local binds
-- Most of them does not have informative names because they are just intermediate calculations.

emStepsFast
  :: (HasCallStack) => Matrix D Double
  -> Matrix D Double
  -> Double
  -> Either Int Double
  -> (Matrix D Double, Double, Double, Maybe (Matrix D Double))
emStepsFast learnMatrix@(ADelayed (Z :. _ :. n) _) initMatrix@(ADelayed (Z :. d :. _) _) initVariance stopParam =
  let learnMatrixCentered = transpose $ substractMean (transpose learnMatrix)

      stepOne :: (HasCallStack) => Array D DIM2 Double -> Double -> Int -> (Matrix D Double, Double, Double, Maybe (Matrix D Double))
      stepOne (!oldW) oldVariance iteration =
        let m = delay $ mapDiagonal (oldVariance +) $ (transpose oldW) `mulS` oldW
            invM = invS m
            expEznZ = delay $ invM `mul` (transpose oldW `mulS` learnMatrixCentered)
            expXnEZnZ = learnMatrixCentered `mulS` (transpose expEznZ)
            expZtrZ = computeS $ (map (*((fromIntegral n)*oldVariance)) $ delay invM) +^ (expEznZ `mulS` (transpose $ expEznZ)) :: Matrix F Double
        in stepTwo (expEznZ, expZtrZ, expXnEZnZ) oldW oldVariance iteration

      stepTwo :: (HasCallStack) => (Matrix D Double, Matrix F Double, Matrix F Double) -> Matrix D Double -> Double -> Int -> (Matrix D Double, Double, Double, Maybe (Matrix D Double))
      stepTwo (expEznZ, expZtrZ, expXnEZnZ) oldW oldVariance iteration =
        let newW = delay $ expXnEZnZ `mul` (inv expZtrZ)
            u = chol $ trustSym expZtrZ
            wr = newW `mulS` (transpose u)
            totalSum = sumAllS $ map (^ (2 :: Int)) learnMatrixCentered
            secondDenominatorInVariance = (*2) $ sumAllS $ expEznZ *^ (delay $ (transpose newW) `mulS` learnMatrixCentered)
            thirdDenominatorInVariance  = sumAllS $ map (^(2 :: Int)) wr
            newVariance = (totalSum - secondDenominatorInVariance + thirdDenominatorInVariance)/(fromIntegral $ n*d)
            maxDiffNewOldW = foldAllS max 0.0 $ map abs $ oldW -^ newW
            diffVariance = abs $ newVariance - oldVariance
            newC = mapDiagonal ((+) newVariance) $ newW `mulS` (transpose newW)
            newInvM = delay $ invS $ mapDiagonal (newVariance +) $ (transpose newW) `mulS` newW
            newInvC = (mapDiagonal (+(1/newVariance))
                        (map (\x -> -x/newVariance) $ (delay $ newW `mulS` newInvM) `mulS` (transpose newW)))
            expLogLikelihood = (-1.0)*((fromIntegral n)/2.0)*((fromIntegral d)*(log $ 2*pi)
              + (log $ detS newC)
              + (trace2S $ computeS $ map (/(fromIntegral $ n-1)) $ newInvC `mulS` (delay $ learnMatrixCentered `mulS` (transpose learnMatrixCentered))))
        in case stopParam of
              Left maxIteration -> if maxIteration > iteration then stepOne newW newVariance (iteration + 1) else (newW,newVariance, expLogLikelihood, Nothing)
              Right stopDiffBetweenIterations -> if (max diffVariance maxDiffNewOldW) > stopDiffBetweenIterations then stepOne newW newVariance (iteration + 1) else (newW,newVariance, expLogLikelihood, Nothing)

  in stepOne initMatrix initVariance 0

emStepsMissed
  :: (HasCallStack) => Matrix D Double
  -> Matrix D Double
  -> Double
  -> Either Int Double
  -> (Matrix D Double, Double, Double, Maybe (Matrix D Double))
emStepsMissed learnMatrix@(ADelayed (Z :. _ :. n) _) initMatrix@(ADelayed (Z :. d :. _) _) initVariance stopParam =
  let initMu = delay $ fromListUnboxed (Z:.d:.1) $ replicate d 0.0 -- extend (Any :. (1 :: Int)) $ meanColumnWithNan learnMatrix

      stepOne :: (HasCallStack) => Array D DIM2 Double -> Array D DIM2 Double -> Double -> Int -> (Matrix D Double, Double, Double, Maybe (Matrix D Double))
      stepOne oldMu (!oldW) oldVariance iteration =
        let yi i = extend (Any :. (1 :: Int)) $ slice learnMatrix (Any :. (i :: Int))
            yList i = toList (yi i)
            unknownIndices i =
              let zipped = zip [0..] (yList i)
              in U.map fst $ filter (isNaN . snd) zipped
            yPi i  = deleteRows (unknownIndices i) (yi i)
            muP i =  deleteRows (unknownIndices i) oldMu
            oldWP i = deleteRows (unknownIndices i) oldW
            mP i = mapDiagonal (+ oldVariance) $ (transpose (oldWP i)) `mulS` (oldWP i)
            invMP i = delay $ invS $ mP i
            expXi i = delay $ (delay $ (delay $ invMP i) `mulS` (transpose (oldWP i))) `mulS` (yPi i -^ muP i)
        in stepTwo ([yPi, expXi, invMP], unknownIndices, oldW) oldVariance iteration

      stepTwo :: (HasCallStack) => ([Int -> Matrix D Double], Int -> [Int], Matrix D Double) -> Double -> Int -> (Matrix D Double, Double, Double, Maybe (Matrix D Double))
      stepTwo ([yP, expXi, invMP], unknownIndices, oldW) oldVariance iteration =
        let expX = foldl1 (\acc i -> acc ++ i) $ U.map expXi [0..(n-1)]
            newMu = extend (Any :. (1 :: Int)) $ meanColumnWithNan $ (learnMatrix -^ (oldW `mulS` expX))

            newW =
              let yj j = extend (Any :. (1 :: Int) :. All) $ slice learnMatrix (Any :. (j :: Int) :. All)
                  yJList j = toList (yj j)
                  unknownColumns j =
                    let zipped = zip [0..] (yJList j)
                    in U.map fst $ filter (isNaN . snd) zipped
                  knownColumns j = [0..(n-1)] \\ (unknownColumns j)
                  expXPj j = deleteColumns (unknownColumns j) expX
                  sumInvMP j = sumListMatrices $ U.map invMP (knownColumns j)
                  gj j = ((expXPj j) `mulS` (transpose $ expXPj j)) +^ (map (*oldVariance) (sumInvMP j))
                  yPj j = deleteColumns (unknownColumns j) (yj j)
                  expXtrXj j = (expXPj j) `mulS` (transpose (map (\x -> x - (newMu ! (Z:.j:.0))) (yPj j)))
                  newWj j = (gj j) `solveS` (delay $ expXtrXj j)
              in transpose $ foldl1 (\acc j -> acc ++ j) $ U.map (\x -> delay $ newWj x) [0..(d-1)]

            newVariance =
              let newWPi i = deleteRows (unknownIndices i) newW
                  newMuP i = deleteRows (unknownIndices i) newMu
                  varianceStep i = sumAllS $
                    (map (^(2 :: Int)) (yP i -^ ((newWPi i) `mulS` (expXi i)) -^ (newMuP i)))
                    +^ (map (*oldVariance) $ diag $ delay $ (delay $ (delay $ newWPi i) `mulS` (invMP i)) `mulS` (transpose $ newWPi i))
                  knownVars = foldAllS (\acc x -> if isNaN x then acc else acc + 1.0) 0.0 learnMatrix
              in (sum $ U.map varianceStep [0..(n-1)])/knownVars

            expLogLikelihood =
              let newMuP i = deleteRows (unknownIndices i) newMu
                  newWPi i = deleteRows (unknownIndices i) newW
                  yCenteredP i = (yP i) -^ (newMuP i)
                  yCenteredProd i = delay $ (yCenteredP i) `mulS` (transpose $ yCenteredP i)
                  invMY i = mapDiagonal (+newVariance) $ (newWPi i) `mulS` (transpose $ newWPi i)
                  knownIndices i = [0..(d-1)] \\ (unknownIndices i)
                  expLogLikelihoodStep i = (fromIntegral $ length $ knownIndices i)*(log $ 2*pi) + (log $ detS $ invMY i) +
                    (trace2S $ computeS $ delay $ (invMY i) `solveS` (yCenteredProd i))

              in (-1.0)*(sum $ U.map expLogLikelihoodStep [0..(n-1)])/2.0

            restoredData =
              let muX = extend (Any :. (1 :: Int)) $ meanColumn (transpose expX)
                  finalMu = slice (newMu +^ (newW `mulS` muX)) (Any :. (0 :: Int))
                  wtrw = (transpose newW) `mulS` newW
                  factor1 = newW `mulS` (delay $ inv wtrw)
                  factor2 = computeS $ mapDiagonal (+newVariance) wtrw
                  restoredCentered =  factor1 `mul` factor2 `mul` (computeS $ expX)
              in (delay restoredCentered) +^ (extend (Any :. (n :: Int)) finalMu)

            diffVariance = abs $ newVariance - oldVariance
            maxDiffNewOldW = foldAllS max 0.0 $ map abs $ oldW -^ newW

        in case stopParam of
              Left maxIteration ->
                if maxIteration > iteration
                then stepOne newMu newW newVariance (iteration + 1)
                else (newW,newVariance, expLogLikelihood, Just restoredData)
              Right stopDiffBetweenIterations ->
                if (max diffVariance maxDiffNewOldW) > stopDiffBetweenIterations
                then stepOne newMu newW newVariance (iteration + 1)
                else (newW,newVariance, expLogLikelihood, Just restoredData)
  in stepOne initMu initMatrix initVariance 0

ext: hs

max_stars_repo_name: serokell/GPLVMHaskell

max_issues_repo_name: serokell/GPLVMHaskell

max_forks_repo_name: serokell/GPLVMHaskell


text:
module Mand where
import Data.List.Split
import Data.Scientific as Scientific
import Data.Complex

type CS = Complex Scientific
type S = Scientific

makeSci :: String -> String -> Scientific
makeSci a b = read (a ++ "e" ++ b) :: Scientific

csqrt :: (Floating a, RealFloat a) => Complex a -> Complex a
csqrt x = (m * cos (t) :+ m * sin(t))
  where
    p = polar x
    m = sqrt $ fst p
    t = (snd p)/2

--cubed_iteration c z  = c + z * cos (z ** z)
cubed_iteration :: RealFloat a => Complex a -> Complex a -> Complex a
cubed_iteration c z = c + z * z * z

new_iteration :: RealFloat a => Complex a -> Complex a -> Complex a -> Complex a
new_iteration c n_z o_z = c + n_z * o_z

zzcosMask :: RealFloat a => Complex a -> Bool
zzcosMask (x :+y) = ((x > -0.95 && x < 0.1) && (y < 0.47 && y > -0.47)) || (( x <= -0.95 && x > -1.25) && (abs y < 0.42)) || (abs y < 0.2 && x <= -1.25 && x > -1.45) || (abs y < 0.4 && x >= 0.1 && x < 0.2) || (abs y < 0.17 && (x > 0.6 && x < 0.78)) || (abs y < 0.11 && x > 0.4 && x <= 0.85)

zzcos_iteration :: RealFloat a => Complex a -> Complex a -> Complex a
zzcos_iteration c z = c + z * z * (cos z)

mand_iteration :: RealFloat a => Complex a -> Complex a -> Complex a
mand_iteration c z  = c + z*z

mand_pow_iteration :: RealFloat a => Int -> Complex a -> Complex a -> Complex a
mand_pow_iteration pow c z  = c + z^pow

mand_iterationA :: CS -> CS -> CS
mand_iterationA (c:+b) (z:+x)  = (c + z*z - x* x) :+ (b + 2*z*x)

cexp :: Floating a => Complex a -> Complex a
cexp (a :+ b) = (s * cos b :+ s * sin b)
  where
    s = exp a

csin :: RealFloat a => Complex a -> Complex a
csin x = ((realPart c1 + realPart c2)/2 :+ (imagPart c1 + imagPart c2)/2)
 where
  c1 = exp x
  c2 = exp $ negate x

cpow :: RealFloat a => Complex a -> Complex a -> Complex a
cpow x y = cexp $ y * log  x

mand :: (RealFloat a) => Either (Complex a) CS -> Int -> Int
mand (Left c) max = general mand_iteration (Just mandSive) c max
mand (Right c) _ = generalA mand_iterationA c

sciMagnitude :: Complex Scientific  -> Scientific
sciMagnitude (a :+b) = a * a + b * b

mandSive :: RealFloat a => Complex a -> Bool
mandSive c@(x:+y) = ( let b = c + (1:+0) in realPart(b*(conjugate b)) < 0.05) ||((x < 0) && ( realPart (c * (conjugate c))  < 0.4)) ||(x < 0.25 && x  >= 0 && abs y  < 0.5)  -- takes 16 seconds

--( let b = a + (1:+0) in realPart(b*(conjugate b)) < 0.05) ||((realPart a < 0) && ( realPart (a * (conjugate a))  < 0.4)) ||(realPart a < 0.25 && realPart a >= 0 && abs (imagPart a) < 0.5) =  max -1  -- takes 16 seconds

general :: RealFloat a => (Complex a -> Complex a -> Complex a) -> Maybe (Complex a -> Bool)  -> Complex a -> Int -> Int
general g Nothing a max = count_iterations 0 (0 :+ 0) a
  where
    count_iterations n e x
      | n == max = max -1
      | realPart (x * (conjugate x)) > 4 = n
      | otherwise = count_iterations (n +1) e $ g a x
general g (Just sive) a max
--  | ((realPart a) < 0.25) && ((realPart a) > (-0.5)) && (abs (imagPart a) < 0.5) = 255 takes 21 secods
-- removed real check: (imagPart a == 0 && realPart a < 0.25 && realPart a > -2) ||
  | sive a = max - 1
  | otherwise = count_iterations 0 (0 :+ 0) a
  where
    count_iterations n e x
      | n == max = max -1
      | realPart (x * (conjugate x)) > 4 = n
      | otherwise = count_iterations (n +1) e $ g a x

generalA :: (CS -> CS -> CS) -> CS -> Int
generalA g a = count_iterations 0 (0 :+ 0) a
  where
    count_iterations n e x
      -- | n >= 255  = 255
      | n >= 16581375 = 16581375
      -- if using sci
      -- | (sciMagnitude x ) > (makeSci "4" "0") = n
      -- | otherwise = count_iterations (n + 1) e  $ g a $ roundComplex x
      -- if using float
      | (realPart x)^2 + (imagPart x)^2 > 4 = n
      | otherwise = count_iterations (n +1) e $ g a $ roundComplex x

general_list_R ::  (Complex Scientific -> Complex Scientific -> Complex Scientific) -> Complex Scientific -> [Complex Scientific]
general_list_R g a = get_list 0 a (0 :+ 0)
  where
    get_list 1000 _ _  = []
    get_list n (c :+d) (e :+ f)
      | (sciMagnitude (e :+ f)) > (scientific 2 0) = []
      | otherwise = (e:+f):  (get_list (n +1) (c :+d) ( g (c :+ d) (roundComplex (e :+ f))))
-- takes 1:13 to do 1000 iterations using roundComplex

general_list ::  (Complex Scientific -> Complex Scientific -> Complex Scientific) -> Complex Scientific -> [Complex Scientific]
general_list g (a :+ b) = get_list 0  (a :+ b) (0 :+ 0)
  where
    get_list 40  _ _  = []
    get_list n (c :+d) (e :+ f)
      | (sciMagnitude (e :+ f)) > (makeSci "2" "2") = []
      | otherwise = (e:+f):  (get_list (n +1) (c :+d) ( g (c :+ d) (e :+ f)))

general_list_D :: RealFloat a => (Complex a -> Complex a -> Complex a) -> Complex a -> [Complex a]
general_list_D g x = getlist 0 x (0 :+ 0)
  where
    getlist 1000 _ _ = []
    getlist n c z
      | realPart (z * conjugate z)> 4 = []
      | otherwise = z: (getlist (n +1) c (g c z))

roundSci' :: Scientific -> Scientific
roundSci' s = makeSci (getSign:cut) getExp
  where
    num = splitOn "e" $ show s
    removeSign = if s < 0 then tail (num !! 0) else (num !! 0)
    getSign = if s < 0 then '-'  else ' '
    getExp = if (length num) == 1 then "0" else num !! 1
    cut = take 100  removeSign

--rounds t 20 dp
roundSci :: Scientific -> Scientific
roundSci s
  | d <  0  = s
  | otherwise = scientific (div (coefficient s) (10^d)) $ d + base10Exponent s
  where
    d = (digitCount $ coefficient s) - 20

digitCount :: Integer -> Int
digitCount = go 1 . abs
    where
        go ds n = if n >= 10 then go (ds + 1) (n `div` 10) else ds
--(fromInteger $ round $ f * (10^n)) / (10.0^^n)
--Integer
roundComplex :: Complex Scientific -> Complex Scientific
roundComplex (a :+ b) = (roundSci' a :+ roundSci' b)

max_stars_repo_name: AlwinHughes/fractal

max_issues_repo_name: AlwinHughes/fractal

max_forks_repo_name: AlwinHughes/fractal


text:
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE MagicHash #-}
{-# LANGUAGE ParallelListComp #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE ViewPatterns #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}

module Data.Store.TH.Internal
    (
    -- * TH functions for generating Store instances
      deriveManyStoreFromStorable
    , deriveTupleStoreInstance
    , deriveGenericInstance
    , deriveManyStorePrimVector
    , deriveManyStoreUnboxVector
    , deriveStore
    , makeStore
    -- * Misc utilties used in Store test
    , getAllInstanceTypes1
    , isMonoType
    ) where

import           Control.Applicative
import           Data.Complex ()
import           Data.Generics.Aliases (extT)
import           Data.Generics.Schemes (listify, everywhere)
import           Data.List (find, nub)
import qualified Data.Map as M
import           Data.Maybe (fromMaybe)
import           Data.Primitive.ByteArray
import           Data.Primitive.Types
import           Data.Store.Core
import           Data.Store.Impl
import qualified Data.Text as T
import           Data.Traversable (forM)
import qualified Data.Vector.Primitive as PV
import qualified Data.Vector.Unboxed as UV
import           Data.Word
import           Foreign.Storable (Storable)
import           GHC.Types (Int(..))
import           Language.Haskell.TH
import           Language.Haskell.TH.ReifyMany.Internal (TypeclassInstance(..), getInstances, unAppsT)
import           Language.Haskell.TH.Syntax (lift)
import           Prelude
import           Safe (headMay)
import           TH.Derive (Deriver(..))
import           TH.ReifySimple
import           TH.Utilities (expectTyCon1, dequalify, plainInstanceD, appsT)

instance Deriver (Store a) where
    runDeriver _ preds ty = do
        argTy <- expectTyCon1 ''Store ty
        dt <- reifyDataTypeSubstituted argTy
        (:[]) <$> deriveStore preds argTy (dtCons dt)

-- | Given the name of a type, generate a Store instance for it,
-- assuming that all type variables also need to be Store instances.
--
-- Note that when used with datatypes that require type variables, the
-- ScopedTypeVariables extension is required.
makeStore :: Name -> Q [Dec]
makeStore name = do
    dt <- reifyDataType name
    let preds = map (storePred . VarT) (dtTvs dt)
        argTy = appsT (ConT name) (map VarT (dtTvs dt))
    (:[]) <$> deriveStore preds argTy (dtCons dt)

deriveStore :: Cxt -> Type -> [DataCon] -> Q Dec
deriveStore preds headTy cons0 =
    makeStoreInstance preds headTy
        <$> sizeExpr
        <*> peekExpr
        <*> pokeExpr
  where
    cons :: [(Name, [(Name, Type)])]
    cons =
      [ ( dcName dc
        , [ (mkName ("c" ++ show ixc ++ "f" ++ show ixf), ty)
          | ixf <- ints
          | (_, ty) <- dcFields dc
          ]
        )
      | ixc <- ints
      | dc <- cons0
      ]
    -- NOTE: tag code duplicated in th-storable.
    (tagType, _, tagSize) =
        fromMaybe (error "Too many constructors") $
        find (\(_, maxN, _) -> maxN >= length cons) tagTypes
    tagTypes :: [(Name, Int, Int)]
    tagTypes =
        [ ('(), 1, 0)
        , (''Word8, fromIntegral (maxBound :: Word8), 1)
        , (''Word16, fromIntegral (maxBound :: Word16), 2)
        , (''Word32, fromIntegral (maxBound :: Word32), 4)
        , (''Word64, fromIntegral (maxBound :: Word64), 8)
        ]
    fName ix = mkName ("f" ++ show ix)
    ints = [0..] :: [Int]
    fNames = map fName ints
    sizeNames = zipWith (\_ -> mkName . ("sz" ++) . show) cons ints
    tagName = mkName "tag"
    valName = mkName "val"
    sizeExpr =
        caseE (tupE (concatMap (map sizeAtType . snd) cons))
              (if null sizeNames then [matchConstSize] else [matchConstSize, matchVarSize])
      where
        sizeAtType :: (Name, Type) -> ExpQ
        sizeAtType (_, ty) = [| size :: Size $(return ty) |]
        matchConstSize :: MatchQ
        matchConstSize = do
            let sz0 = VarE (mkName "sz0")
                sizeDecls =
                    if null sizeNames
                        then [valD (varP (mkName "sz0")) (normalB [| 0 |]) []]
                        else zipWith constSizeDec sizeNames cons
            sameSizeExpr <-
                case sizeNames of
                    (_ : tailSizeNames) ->
                        foldl (\l r -> [| $(l) && $(r) |]) [| True |] $
                        map (\szn -> [| $(return sz0) == $(varE szn) |]) tailSizeNames
                    [] -> [| True |]
            result <- [| ConstSize (tagSize + $(return sz0)) |]
            match (tupP (map (\(n, _) -> conP 'ConstSize [varP n])
                             (concatMap snd cons)))
                  (guardedB [return (NormalG sameSizeExpr, result)])
                  sizeDecls
        constSizeDec :: Name -> (Name, [(Name, Type)]) -> DecQ
        constSizeDec szn (_, []) =
            valD (varP szn) (normalB [| 0 |]) []
        constSizeDec szn (_, fields) =
            valD (varP szn) body []
          where
            body = normalB $
                foldl1 (\l r -> [| $(l) + $(r) |]) $
                map (\(sizeName, _) -> varE sizeName) fields
        matchVarSize :: MatchQ
        matchVarSize = do
            match (tupP (map (\(n, _) -> varP n) (concatMap snd cons)))
                  (normalB [| VarSize $ \x -> tagSize +
                                  $(caseE [| x |] (map matchVar cons)) |])
                  []
        matchVar :: (Name, [(Name, Type)]) -> MatchQ
        matchVar (cname, []) =
            match (conP cname []) (normalB [| 0 |]) []
        matchVar (cname, fields) =
            match (conP cname (zipWith (\_ fn -> varP fn) fields fNames))
                  body
                  []
          where
            body = normalB $
                foldl1 (\l r -> [| $(l) + $(r) |])
                (zipWith (\(sizeName, _) fn -> [| getSizeWith $(varE sizeName) $(varE fn) |])
                         fields
                         fNames)
    -- Choose a tag size large enough for this constructor count.
    -- Expression used for the definition of peek.
    peekExpr = case cons of
        [] -> [| error ("Attempting to peek type with no constructors (" ++ $(lift (show headTy)) ++ ")") |]
        [con] -> peekCon con
        _ -> doE
            [ bindS (varP tagName) [| peek |]
            , noBindS (caseE (sigE (varE tagName) (conT tagType))
                      (map peekMatch (zip [0..] cons) ++ [peekErr]))
            ]
    peekMatch (ix, con) = match (litP (IntegerL ix)) (normalB (peekCon con)) []
    peekErr = match wildP (normalB
        [| peekException $ T.pack $ "Found invalid tag while peeking (" ++ $(lift (show headTy)) ++ ")" |]) []
    peekCon (cname, fields) =
        case fields of
            [] -> [| pure $(conE cname) |]
            _ -> doE $
                map (\(fn, _) -> bindS (varP fn) [| peek |]) fields ++
               [noBindS $ appE (varE 'return) $ appsE $ conE cname : map (\(fn, _) -> varE fn) fields]
    pokeExpr = lamE [varP valName] $ caseE (varE valName) $ zipWith pokeCon [0..] cons
    pokeCon :: Int -> (Name, [(Name, Type)]) -> MatchQ
    pokeCon ix (cname, fields) =
        match (conP cname (map (\(fn, _) -> varP fn) fields)) body []
      where
        body = normalB $
            case cons of
                (_:_:_) -> doE (pokeTag ix : map pokeField fields)
                _ -> doE (map pokeField fields)
    pokeTag ix = noBindS [| poke (ix :: $(conT tagType)) |]
    pokeField (fn, _) = noBindS [| poke $(varE fn) |]

-- FIXME: make this work even when there are too many fields

-- FIXME: make the ConstSize stuff explicit in the API. Have an option
-- that always errors at runtime if it isn't ConstSize?

-- TODO: It would be really awesome, though a bit tricky, to know at
-- compile time if we have a static size.

-- TODO: make sure that this tends to optimize even with tons of fields.
-- It should also optimize when some fields are known to be , but others
-- are unknown (determined by polymorphic var)

{- What the generated code looks like

data Foo = Foo Int Double Float

instance Store Foo where
    size =
        case (size :: Size Int, size :: Size Double, size :: Size Float) of
            (ConstSize c0f0, ConstSize c0f1, ConstSize c0f2) -> ConstSize (0 + sz0)
              where
                sz0 = c0f0 + c0f1 + c0f2
            (c0f0, c0f1, c0f2)
                VarSize $ \(Foo f0 f1 f2) -> 0 +
                    getSizeWith c0f0 f0 + getSizeWith c0f1 f1 + getSizeWith c0f2 f2
    peek = do
        f0 <- peek
        f1 <- peek
        f2 <- peek
        return (Foo f0 f1 f2)
    poke (Foo f0 f1 f2) = do
        poke f0
        poke f1
        poke f2

data Bar = Bar Int | Baz Double

instance Store Bar where
    size =
        case (size :: Size Int, size :: Size Double) of
            (ConstSize c0f0, ConstSize c1f0) | sz0 == sz1 -> ConstSize (1 + sz0)
              where
                sz0 = c0f0
                sz1 = c1f0
            (c0f0, c1f0) -> VarSize $ \x -> 1 +
                case x of
                    Bar f0 -> getSizeWith c0f0 f0
                    Baz f0 -> getSizeWith c1f0 f0
    peek = do
        tag <- peek
        case (tag :: Word8) of
            0 -> do
                f0 <- peek
                return (Bar f0)
            1 -> do
                f0 <- peek
                return (Baz f0)
            _ -> peekException "Found invalid tag while peeking (Bar)"
    poke (Bar f0) = do
        poke 0
        poke f0
    poke (Bar f0) = do
        poke 1
        poke f0
-}

------------------------------------------------------------------------
-- Generic

deriveTupleStoreInstance :: Int -> Dec
deriveTupleStoreInstance n =
    deriveGenericInstance (map storePred tvs)
                          (foldl1 AppT (TupleT n : tvs))
  where
    tvs = take n (map (VarT . mkName . (:[])) ['a'..'z'])

deriveGenericInstance :: Cxt -> Type -> Dec
deriveGenericInstance cs ty = plainInstanceD cs (AppT (ConT ''Store) ty) []

------------------------------------------------------------------------
-- Storable

-- TODO: Generate inline pragmas? Probably not necessary

deriveManyStoreFromStorable :: (Type -> Bool) -> Q [Dec]
deriveManyStoreFromStorable p = do
    storables <- postprocess . instancesMap <$> getInstances ''Storable
    stores <- postprocess . instancesMap <$> getInstances ''Store
    return $ M.elems $ flip M.mapMaybe (storables `M.difference` stores) $
        \(TypeclassInstance cs ty _) ->
        let argTy = head (tail (unAppsT ty))
            tyNameLit = LitE (StringL (pprint ty)) in
        if p argTy
            then Just $ makeStoreInstance cs argTy
                (AppE (VarE 'sizeStorableTy) tyNameLit)
                (AppE (VarE 'peekStorableTy) tyNameLit)
                (VarE 'pokeStorable)
            else Nothing

------------------------------------------------------------------------
-- Vector

deriveManyStorePrimVector :: Q [Dec]
deriveManyStorePrimVector = do
    prims <- postprocess . instancesMap <$> getInstances ''PV.Prim
    stores <- postprocess . instancesMap <$> getInstances ''Store
    let primInsts =
            M.mapKeys (map (AppT (ConT ''PV.Vector))) prims
            `M.difference`
            stores
    forM (M.toList primInsts) $ \primInst -> case primInst of
        ([_], TypeclassInstance cs ty _) -> do
            let argTy = head (tail (unAppsT ty))
            sizeExpr <- [|
                VarSize $ \x ->
                    I# $(primSizeOfExpr (ConT ''Int)) +
                    I# $(primSizeOfExpr argTy) * PV.length x
                |]
            peekExpr <- [| do
                len <- peek
                let sz = I# $(primSizeOfExpr argTy)
                array <- peekToByteArray $(lift ("Primitive Vector (" ++ pprint argTy ++ ")"))
                                         (len * sz)
                return (PV.Vector 0 len array)
                |]
            pokeExpr <- [| \(PV.Vector offset len (ByteArray array)) -> do
                let sz = I# $(primSizeOfExpr argTy)
                poke len
                pokeFromByteArray array (offset * sz) (len * sz)
                |]
            return $ makeStoreInstance cs (AppT (ConT ''PV.Vector) argTy) sizeExpr peekExpr pokeExpr
        _ -> fail "Invariant violated in derivemanyStorePrimVector"

primSizeOfExpr :: Type -> ExpQ
primSizeOfExpr ty = [| $(varE 'sizeOf#) (error "sizeOf# evaluated its argument" :: $(return ty)) |]

deriveManyStoreUnboxVector :: Q [Dec]
deriveManyStoreUnboxVector = do
    unboxes <- getUnboxInfo
    stores <- postprocess . instancesMap <$> getInstances ''Store
    let unboxInsts =
            M.fromList (map (\(preds, ty, cons) -> ([AppT (ConT ''UV.Vector) ty], (preds, cons))) unboxes)
            `M.difference`
            stores
    -- TODO: ideally this would use a variant of 'deriveStore' which
    -- assumes VarSize.
    forM (M.toList unboxInsts) $ \case
        ([ty], (preds, cons)) -> do
            {-
            -- While this approach is reasonable-ish, it ends up
            -- requiring UndecidableInstances.
            let extraPreds =
                    map (AppT (ConT ''Store)) $
                    filter (not . isMonoType) $
                    concatMap (map snd . dcFields) cons
            -}
            let extraPreds = map (storePred . AppT (ConT ''UV.Vector)) $ listify isVarT ty
            deriveStore (nub (preds ++ extraPreds)) ty cons
        _ -> fail "impossible case in deriveManyStoreUnboxVector"

-- TODO: Add something for this purpose to TH.ReifyDataType

getUnboxInfo :: Q [(Cxt, Type, [DataCon])]
getUnboxInfo = do
    FamilyI _ insts <- reify ''UV.Vector
    return (map (everywhere (id `extT` dequalVarT) . go) insts)
  where
#if MIN_VERSION_template_haskell(2,11,0)
    go (NewtypeInstD preds _ [ty] _ con _) = (preds, ty, conToDataCons con)
    go (DataInstD preds _ [ty] _ cons _) = (preds, ty, concatMap conToDataCons cons)
#else
    go (NewtypeInstD preds _ [ty] con _) = (preds, ty, conToDataCons con)
    go (DataInstD preds _ [ty] cons _) = (preds, ty, concatMap conToDataCons cons)
#endif
    go x = error ("Unexpected result from reifying Unboxed Vector instances: " ++ pprint x)
    dequalVarT (VarT n) = VarT (dequalify n)
    dequalVarT ty = ty

------------------------------------------------------------------------
-- Utilities

-- Filters out overlapping instances and instances with more than one
-- type arg (should be impossible).
postprocess :: M.Map [Type] [a] -> M.Map [Type] a
postprocess =
    M.mapMaybeWithKey $ \tys xs ->
        case (tys, xs) of
            ([_ty], [x]) -> Just x
            _ -> Nothing

makeStoreInstance :: Cxt -> Type -> Exp -> Exp -> Exp -> Dec
makeStoreInstance cs ty sizeExpr peekExpr pokeExpr =
    plainInstanceD
        cs
        (AppT (ConT ''Store) ty)
        [ ValD (VarP 'size) (NormalB sizeExpr) []
        , ValD (VarP 'peek) (NormalB peekExpr) []
        , ValD (VarP 'poke) (NormalB pokeExpr) []
        ]

-- TODO: either generate random types that satisfy instances with
-- variables in them, or have a check that there's at least a manual
-- check for polymorphic instances.

getAllInstanceTypes :: Name -> Q [[Type]]
getAllInstanceTypes n =
    map (\(TypeclassInstance _ ty _) -> drop 1 (unAppsT ty)) <$>
    getInstances n

getAllInstanceTypes1 :: Name -> Q [Type]
getAllInstanceTypes1 n =
    fmap (fmap (fromMaybe (error "getAllMonoInstances1 expected only one type argument") . headMay))
         (getAllInstanceTypes n)

isMonoType :: Type -> Bool
isMonoType = null . listify isVarT

isVarT :: Type -> Bool
isVarT VarT{} = True
isVarT _ = False

-- TOOD: move these to th-reify-many

-- | Get a map from the 'getTyHead' type of instances to
-- 'TypeclassInstance'.
instancesMap :: [TypeclassInstance] -> M.Map [Type] [TypeclassInstance]
instancesMap =
    M.fromListWith (++) .
    map (\ti -> (map getTyHead (instanceArgTypes ti), [ti]))

instanceArgTypes :: TypeclassInstance -> [Type]
instanceArgTypes (TypeclassInstance _ ty _) = drop 1 (unAppsT ty)

getTyHead :: Type -> Type
getTyHead (SigT x _) = getTyHead x
getTyHead (ForallT _ _ x) = getTyHead x
getTyHead (AppT l _) = getTyHead l
getTyHead x = x

storePred :: Type -> Pred
storePred ty =
#if MIN_VERSION_template_haskell(2,10,0)
        AppT (ConT ''Store) ty
#else
        ClassP ''Store [ty]
#endif

max_stars_repo_name: TravisWhitaker/store

max_issues_repo_name: TravisWhitaker/store

max_forks_repo_name: TravisWhitaker/store


text:
module Main where

import Data.Ord (comparing)
import qualified Data.Vector.Algorithms.Intro as Intro
import qualified Data.Vector.Storable as VS
import Data.Vector.Storable.Mutable (IOVector)
import Numeric.LinearAlgebra hiding (eig, vector)
import Numeric.LinearAlgebra.Arnoldi
import Test.Hspec
import Test.QuickCheck

main :: IO ()
main = hspec $ do
  describe "Numeric.LinearAlgebra.Arnoldi.eig" $ do

    it "returns the diagonal of an arbitrary diagonal matrix" $ property $ do
      dim <- arbitrary `suchThat` (> 3)
      diagonal <- VS.fromList <$> vector dim
      let
        _matrix = diag diagonal :: Matrix Double
        nev = dim - 3
        options = Options { which = SR
                          , number = nev
                          , maxIterations = Nothing
                          }
        sort = VS.modify Intro.sort
        actual = (sort . fst) (eig options dim (multiply _matrix))
        expected = (VS.take nev . sort) diagonal
        relative = VS.maximum (VS.zipWith (%) expected actual)
      pure (counterexample (show (expected, actual)) (relative < 1E-4))

    it "returns the diagonal of an arbitrary diagonal matrix" $ property $ do
      dim <- arbitrary `suchThat` (> 3)
      diagonal <- VS.fromList <$> vector dim
      let
        _matrix = diag diagonal :: Matrix (Complex Double)
        nev = dim - 3
        options = Options { which = SR
                          , number = nev
                          , maxIterations = Nothing
                          }
        sort = VS.modify (Intro.sortBy (comparing realPart))
        actual = (sort . fst) (eig options dim (multiply _matrix))
        expected = (VS.take nev . sort) diagonal
        relative = VS.maximum (VS.zipWith (%%) expected actual)
      pure (counterexample (show (expected, actual)) (relative < 1E-4))

multiply :: Numeric a => Matrix a -> IOVector a -> IOVector a -> IO ()
multiply _matrix dst src = do
  x <- VS.freeze src
  let y = _matrix #> x
  VS.copy dst y

(%) :: Double -> Double -> Double
(%) a b =
  let a_plus_b = a + b
  in if a_plus_b /= 0
     then abs ((a - b) / a_plus_b)
     else a - b

(%%) :: Complex Double -> Complex Double -> Double
(%%) a b =
  let a_plus_b = a + b
  in if a_plus_b /= 0
     then magnitude ((a - b) / a_plus_b)
     else magnitude (a - b)

max_stars_repo_name: ttuegel/arpack

max_issues_repo_name: ttuegel/arpack

max_forks_repo_name: ttuegel/arpack


text:
module Main where

import Parser
import Analyzer
import qualified Data.MultiSet as M
import Control.Lens hiding ((??),(|>),(<.>))
import qualified Data.Text.Lazy.IO as I
import Data.List (sortBy)
import Data.Function (on)
import Text.Megaparsec (parseMaybe)
import Numeric.LinearAlgebra
import Data.Complex (realPart)

file :: String
file = "/home/bertram/Downloads/la famiglia.txt"

messages :: IO [Chunk]
messages = filter isMessage <$> ((parseMaybe chunksP <$> I.readFile file)
           >>= maybe (putStrLn "Parse failed" >> return []) return)

main' :: IO ()
main' = do
  msg <- messages
  let users = map (toText . (^.user)) msg
  let counts = M.fromList users
  mapM_ (\(a,b)->I.putStr a >> putStr " | " >> print b) . sortBy (compare `on` snd) . M.toOccurList $ counts

main :: IO ()
main = do
  input <- I.readFile file
  let (matrix,users) = maybe undefined id (parseMaybe matrixP input)
  let n = length users
  let probs =  cmap realPart . flatten . (?? (All, Take 1))
               . snd . eig . tr $ matrix
  let probs' = scale (1 / (probs <.> (n |> repeat 1))) probs
  mapM_ print . sortBy (compare `on` snd). zip users $ toList probs'
  print $ sum . toList $ probs'

max_stars_repo_name: bnarnold/WAAnalyzer

max_issues_repo_name: bnarnold/WAAnalyzer

max_forks_repo_name: bnarnold/WAAnalyzer


text:
{-# LANGUAGE Safe #-}
{-# OPTIONS_GHC -interactive-print=Text.Pretty.Simple.pPrint #-}
{-# LANGUAGE ConstrainedClassMethods #-}
{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE DeriveFoldable #-}
{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE DeriveTraversable #-}
{-# LANGUAGE EmptyDataDecls #-}
{-# LANGUAGE ExistentialQuantification #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE ImplicitParams #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE LiberalTypeSynonyms #-}
{-# LANGUAGE MagicHash #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE PackageImports #-}
{-# LANGUAGE ParallelListComp #-}
{-# LANGUAGE PatternGuards #-}
{-# LANGUAGE PolymorphicComponents #-}
{-# LANGUAGE PostfixOperators #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE UnicodeSyntax #-}
module L where
import Control.Applicative
import Control.Arrow
import Control.Monad
import Control.Monad.Cont
import Control.Monad.Identity
import Control.Monad.Reader
import Control.Monad.ST
import Control.Monad.State
import Control.Monad.Writer
import Data.Bits
import Data.Bool
import Data.Char
import Data.Complex
import Data.Dynamic
import Data.Either
import Data.Eq
import Data.Function
import Data.Functor
import Data.Int
import Data.List
import Data.Maybe
import Data.Monoid
import Data.Ord
import Data.Ratio
import Data.STRef
import Data.Tree
import Data.Tuple
import Data.Typeable
import Data.Word
import Numeric
import ShowFun
import System.Random
import Lambdabot.Plugin.Haskell.Eval.Trusted
import Text.PrettyPrint.HughesPJ hiding (empty)
import Text.Printf
import qualified Data.ByteString as BS
import qualified Data.ByteString.Char8 as BSC
import qualified Data.ByteString.Lazy as BSL
import qualified Data.ByteString.Lazy.Char8 as BSLC
import qualified Data.Foldable
import qualified Data.IntMap as IM
import qualified Data.IntSet as IS
import qualified Data.Map as M
import qualified Data.Sequence
import qualified Data.Set as S
import qualified Data.Traversable

{-# LINE 1 "<local>" #-}

max_stars_repo_name: swamp-agr/lambdabot-telegram-plugins

max_issues_repo_name: swamp-agr/lambdabot-telegram-plugins

max_forks_repo_name: swamp-agr/lambdabot-telegram-plugins


text:
{-# LANGUAGE FlexibleContexts #-}

module Spectral(layout,default_layout) where

import Numeric.LinearAlgebra.HMatrix
import Numeric.LinearAlgebra.Data

degree_matrix m =
  let a = foldl1 (+) (toRows m)
      b = takeDiag m
      in
    a - b

compute_laplacian adj_matrix node_weights =
  let (n,m) = size adj_matrix
      weighted_adjency = adj_matrix
      mat_D = diag $ degree_matrix weighted_adjency
      mat_L = mat_D - weighted_adjency
      in
      (mat_L,mat_D)

symmetrize m = m + tr m

layout adjacent_matrix node_weights =
  let sym_adjacent_matrix = symmetrize adjacent_matrix
      (mat_L,mat_D) = compute_laplacian sym_adjacent_matrix node_weights
      eigenVectors = toColumns $ snd $ geigSH' mat_L mat_D
      v2 = eigenVectors !! 1
      v3 = eigenVectors !! 2
  in
    [(v2  ! i, v3 ! i) | i <- [0..size v2-1] ]

default_layout adjacent_matrix =
  let (n,_) = size adjacent_matrix
      ones = matrix n [1 | i <- [1..n]]
  in
    layout adjacent_matrix ones

max_stars_repo_name: Davidbrcz/renderer

max_issues_repo_name: Davidbrcz/renderer

max_forks_repo_name: Davidbrcz/renderer


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE DeriveTraversable #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE RankNTypes #-}
#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 702
{-# LANGUAGE Trustworthy #-}
{-# LANGUAGE DeriveGeneric #-}
#endif
{-# LANGUAGE DeriveDataTypeable #-}
#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ < 708
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE ScopedTypeVariables #-}
#endif
-----------------------------------------------------------------------------
-- |
-- License     :  BSD-style (see the file LICENSE)
-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
-- Stability   :  provisional
-- Portability :  portable
--
-- Operations on affine spaces.
-----------------------------------------------------------------------------
module Linear.Affine where

import Control.Applicative
import Control.DeepSeq
import Control.Monad (liftM)
import Control.Lens
import Data.Binary as Binary
import Data.Bytes.Serial
import Data.Complex (Complex)
import Data.Data
import Data.Distributive
import Data.Foldable as Foldable
import Data.Functor.Bind
import Data.Functor.Classes
import Data.Functor.Rep as Rep
import Data.HashMap.Lazy (HashMap)
import Data.Hashable
import Data.IntMap (IntMap)
import Data.Ix
import Data.Map (Map)
import Data.Serialize as Cereal
import Data.Vector (Vector)
import Foreign.Storable
#if __GLASGOW_HASKELL__ >= 702
import GHC.Generics (Generic)
#endif
#if __GLASGOW_HASKELL__ >= 706
import GHC.Generics (Generic1)
#endif
import qualified Data.Vector.Generic.Mutable as M
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed.Base as U
import Linear.Epsilon
import Linear.Metric
import Linear.Plucker
import Linear.Quaternion
import Linear.V
import Linear.V0
import Linear.V1
import Linear.V2
import Linear.V3
import Linear.V4
import Linear.Vector

#ifdef HLINT
{-# ANN module "HLint: ignore Unused LANGUAGE pragma" #-}
#endif

-- | An affine space is roughly a vector space in which we have
-- forgotten or at least pretend to have forgotten the origin.
--
-- > a .+^ (b .-. a)  =  b@
-- > (a .+^ u) .+^ v  =  a .+^ (u ^+^ v)@
-- > (a .-. b) ^+^ v  =  (a .+^ v) .-. q@
class Additive (Diff p) => Affine p where
  type Diff p :: * -> *

  infixl 6 .-.
  -- | Get the difference between two points as a vector offset.
  (.-.) :: Num a => p a -> p a -> Diff p a

  infixl 6 .+^
  -- | Add a vector offset to a point.
  (.+^) :: Num a => p a -> Diff p a -> p a

  infixl 6 .-^
  -- | Subtract a vector offset from a point.
  (.-^) :: Num a => p a -> Diff p a -> p a
  p .-^ v = p .+^ negated v
  {-# INLINE (.-^) #-}

-- | Compute the quadrance of the difference (the square of the distance)
qdA :: (Affine p, Foldable (Diff p), Num a) => p a -> p a -> a
qdA a b = Foldable.sum (fmap (join (*)) (a .-. b))
{-# INLINE qdA #-}

-- | Distance between two points in an affine space
distanceA :: (Floating a, Foldable (Diff p), Affine p) => p a -> p a -> a
distanceA a b = sqrt (qdA a b)
{-# INLINE distanceA #-}

#define ADDITIVEC(CTX,T) instance CTX => Affine T where type Diff T = T ; \
  (.-.) = (^-^) ; {-# INLINE (.-.) #-} ; (.+^) = (^+^) ; {-# INLINE (.+^) #-} ; \
  (.-^) = (^-^) ; {-# INLINE (.-^) #-}
#define ADDITIVE(T) ADDITIVEC((), T)

ADDITIVE([])
ADDITIVE(Complex)
ADDITIVE(ZipList)
ADDITIVE(Maybe)
ADDITIVE(IntMap)
ADDITIVE(Identity)
ADDITIVE(Vector)
ADDITIVE(V0)
ADDITIVE(V1)
ADDITIVE(V2)
ADDITIVE(V3)
ADDITIVE(V4)
ADDITIVE(Plucker)
ADDITIVE(Quaternion)
ADDITIVE(((->) b))
ADDITIVEC(Ord k, (Map k))
ADDITIVEC((Eq k, Hashable k), (HashMap k))
ADDITIVEC(Dim n, (V n))

-- | A handy wrapper to help distinguish points from vectors at the
-- type level
newtype Point f a = P (f a)
  deriving ( Eq, Ord, Show, Read, Monad, Functor, Applicative, Foldable
           , Eq1, Ord1, Show1, Read1
           , Traversable, Apply, Additive, Metric
           , Fractional , Num, Ix, Storable, Epsilon
           , Hashable
#if __GLASGOW_HASKELL__ >= 702
           , Generic
#endif
#if __GLASGOW_HASKELL__ >= 706
           , Generic1
#endif
#if __GLASGOW_HASKELL__ >= 708
           , Typeable, Data
#endif
           )

instance NFData (f a) => NFData (Point f a) where
  rnf (P x) = rnf x

instance Serial1 f => Serial1 (Point f) where
  serializeWith f (P p) = serializeWith f p
  deserializeWith m = P `liftM` deserializeWith m

instance Serial (f a) => Serial (Point f a) where
  serialize (P p) = serialize p
  deserialize = P `liftM` deserialize

instance Binary (f a) => Binary (Point f a) where
  put (P p) = Binary.put p
  get = P `liftM` Binary.get

instance Serialize (f a) => Serialize (Point f a) where
  put (P p) = Cereal.put p
  get = P `liftM` Cereal.get

#if __GLASGOW_HASKELL__ < 708
instance forall f. Typeable1 f => Typeable1 (Point f) where
  typeOf1 _ = mkTyConApp (mkTyCon3 "linear" "Linear.Affine" "Point") [] `mkAppTy`
              typeOf1 (undefined :: f a)

deriving instance (Data (f a), Typeable1 f, Typeable a) => Data (Point f a)
#endif

lensP :: Lens' (Point g a) (g a)
lensP afb (P a) = P <$> afb a
{-# INLINE lensP #-}

_Point :: Iso' (Point f a) (f a)
_Point = iso (\(P a) -> a) P
{-# INLINE _Point #-}

instance (t ~ Point g b) => Rewrapped (Point f a) t
instance Wrapped (Point f a) where
  type Unwrapped (Point f a) = f a
  _Wrapped' = _Point
  {-# INLINE _Wrapped' #-}

instance Bind f => Bind (Point f) where
  join (P m) = P $ join $ fmap (\(P m')->m') m

instance Distributive f => Distributive (Point f) where
  distribute = P . collect (\(P p) -> p)

instance Representable f => Representable (Point f) where
  type Rep (Point f) = Rep f
  tabulate f = P (tabulate f)
  {-# INLINE tabulate #-}
  index (P xs) = Rep.index xs
  {-# INLINE index #-}

type instance Index (Point f a) = Index (f a)
type instance IxValue (Point f a) = IxValue (f a)

instance Ixed (f a) => Ixed (Point f a) where
  ix l = lensP . ix l
  {-# INLINE ix #-}

instance Traversable f => Each (Point f a) (Point f b) a b where
  each = traverse
  {-# INLINE each #-}

instance R1 f => R1 (Point f) where
  _x = lensP . _x
  {-# INLINE _x #-}

instance R2 f => R2 (Point f) where
  _y = lensP . _y
  {-# INLINE _y #-}
  _xy = lensP . _xy
  {-# INLINE _xy #-}

instance R3 f => R3 (Point f) where
  _z = lensP . _z
  {-# INLINE _z #-}
  _xyz = lensP . _xyz
  {-# INLINE _xyz #-}

instance R4 f => R4 (Point f) where
  _w = lensP . _w
  {-# INLINE _w #-}
  _xyzw = lensP . _xyzw
  {-# INLINE _xyzw #-}

instance Additive f => Affine (Point f) where
  type Diff (Point f) = f
  P x .-. P y = x ^-^ y
  {-# INLINE (.-.) #-}
  P x .+^ v = P (x ^+^ v)
  {-# INLINE (.+^) #-}
  P x .-^ v = P (x ^-^ v)
  {-# INLINE (.-^) #-}

-- | Vector spaces have origins.
origin :: (Additive f, Num a) => Point f a
origin = P zero

-- | An isomorphism between points and vectors, given a reference
--   point.
relative :: (Additive f, Num a) => Point f a -> Iso' (Point f a) (f a)
relative p0 = iso (.-. p0) (p0 .+^)
{-# INLINE relative #-}

data instance U.Vector    (Point f a) =  V_P !(U.Vector    (f a))
data instance U.MVector s (Point f a) = MV_P !(U.MVector s (f a))
instance U.Unbox (f a) => U.Unbox (Point f a)

instance U.Unbox (f a) => M.MVector U.MVector (Point f a) where
  {-# INLINE basicLength #-}
  {-# INLINE basicUnsafeSlice #-}
  {-# INLINE basicOverlaps #-}
  {-# INLINE basicUnsafeNew #-}
  {-# INLINE basicUnsafeRead #-}
  {-# INLINE basicUnsafeWrite #-}
  basicLength (MV_P v) = M.basicLength v
  basicUnsafeSlice m n (MV_P v) = MV_P (M.basicUnsafeSlice m n v)
  basicOverlaps (MV_P v) (MV_P u) = M.basicOverlaps v u
  basicUnsafeNew n = MV_P `liftM` M.basicUnsafeNew n
  basicUnsafeRead (MV_P v) i = P `liftM` M.basicUnsafeRead v i
  basicUnsafeWrite (MV_P v) i (P x) = M.basicUnsafeWrite v i x
#if MIN_VERSION_vector(0,11,0)
  basicInitialize (MV_P v) = M.basicInitialize v
  {-# INLINE basicInitialize #-}
#endif

instance U.Unbox (f a) => G.Vector U.Vector (Point f a) where
  {-# INLINE basicUnsafeFreeze #-}
  {-# INLINE basicUnsafeThaw   #-}
  {-# INLINE basicLength       #-}
  {-# INLINE basicUnsafeSlice  #-}
  {-# INLINE basicUnsafeIndexM #-}
  basicUnsafeFreeze (MV_P v) = V_P `liftM` G.basicUnsafeFreeze v
  basicUnsafeThaw   ( V_P v) = MV_P `liftM` G.basicUnsafeThaw   v
  basicLength       ( V_P v) = G.basicLength v
  basicUnsafeSlice m n (V_P v) = V_P (G.basicUnsafeSlice m n v)
  basicUnsafeIndexM (V_P v) i = P `liftM` G.basicUnsafeIndexM v i

max_stars_repo_name: bgamari/linear

max_issues_repo_name: bgamari/linear

max_forks_repo_name: bgamari/linear


text:
{-# LANGUAGE DataKinds
           , GADTs
           , KindSignatures
           , OverloadedLists
           , ScopedTypeVariables
           , TypeApplications #-}

module MorphGine.Data.Primitives where

import qualified Data.Vector.Storable  as VS
import qualified Numeric.LinearAlgebra as L
import           Numeric.LinearAlgebra ((><))
import qualified Numeric.LinearAlgebra.Data as DL --(flatten, reshape,cmap)
import           GHC.Ptr (Ptr(..))
import           GHC.TypeLits
import qualified Graphics.Rendering.OpenGL.GL as GL
import           Graphics.Rendering.OpenGL.GL (($=))
import qualified Graphics.GL as GLRaw
import           Graphics.GL (GLfloat)

data Matrix :: * -> Nat -> Nat -> * where
  M :: VS.Storable a => L.Matrix a -> Matrix a m n

withPtr :: (VS.Storable a, L.Element a) => Matrix a m n -> (Ptr a -> IO b) -> IO b
withPtr (M x) = VS.unsafeWith . L.flatten $ x

projection :: GLfloat -> GLfloat -> GLfloat -> GLfloat -> Matrix GLfloat 4 4
projection fov aspect near far = M . (4><4) $ [ t/aspect, 0, 0, 0
                                               , 0,        t, 0, 0
                                               , 0,        0, (far+near)/(near-far), (2*far*near)/(near-far)
                                               , 0,        0,-1, 0]
  where t = 1 / tan (fov/2)

degProjection :: GLfloat -> GLfloat -> GLfloat -> GLfloat -> Matrix  GLfloat 4 4
degProjection = projection . (/180) . (*pi)

max_stars_repo_name: Antystenes/CPG

max_issues_repo_name: Antystenes/CPG

max_forks_repo_name: Antystenes/CPG


text:
{-# LANGUAGE FlexibleContexts #-}
module R2Z2T0S0ToR2S1RPT0S0 where

import           Control.Monad.Parallel    as MP
import           Data.Array.Repa           as R
import           Data.Binary               (decodeFile)
import           Data.Complex
import           Data.List                 as L
import           FokkerPlanck.DomainChange
import           FokkerPlanck.MonteCarlo
import           FokkerPlanck.Pinwheel
import           Image.IO
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Utils.Array
import           Utils.Time

main = do
  args@(numPointStr:numOrientationStr:numScaleStr:thetaSigmaStr:scaleSigmaStr:maxScaleStr:taoStr:lenStr:initStr:numTrailStr:maxTrailStr:theta0FreqsStr:thetaFreqsStr:scale0FreqsStr:scaleFreqsStr:histPath:alphaStr:pinwheelFlagStr:numThreadStr:_) <-
    getArgs
  print args
  let numPoint = read numPointStr :: Int
      numOrientation = read numOrientationStr :: Int
      numScale = read numScaleStr :: Int
      thetaSigma = read thetaSigmaStr :: Double
      scaleSigma = read scaleSigmaStr :: Double
      maxScale = read maxScaleStr :: Double
      tao = read taoStr :: Double
      len = read lenStr :: Int
      init = read initStr :: (Double, Double, Double, Double, Double, Double)
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      theta0Freq = read theta0FreqsStr :: Double
      thetaFreq = read thetaFreqsStr :: Double
      scale0Freq = read scale0FreqsStr :: Double
      scaleFreq = read scaleFreqsStr :: Double
      theta0Freqs = [-theta0Freq .. theta0Freq]
      thetaFreqs = [-thetaFreq .. thetaFreq]
      scale0Freqs = [-scale0Freq .. scale0Freq]
      scaleFreqs = [-scaleFreq .. scaleFreq]
      alpha = read alphaStr :: Double
      pinwheelFlag = read pinwheelFlagStr :: Bool
      numThread = read numThreadStr :: Int
      folderPath = "output/test/R2Z2T0S0ToR2S1RPT0S0"
  flag <- doesFileExist histPath
  -- arrR2Z2T0S0 <-
  --   -- if flag
  --   --   then getNormalizedHistogramArr <$> decodeFile histPath
  --   --   else
  --     solveMonteCarloR2Z2T0S0
  --            numThread
  --            numTrail
  --            maxTrail
  --            numPoint
  --            numPoint
  --            thetaSigma
  --            scaleSigma
  --            maxScale
  --            tao
  --            len
  --            theta0Freqs
  --            thetaFreqs
  --            scale0Freqs
  --            scaleFreqs
  --            histPath
  --            (emptyHistogram
  --               [ numPoint
  --               , numPoint
  --               , L.length scale0Freqs
  --               , L.length theta0Freqs
  --               , L.length scaleFreqs
  --               , L.length thetaFreqs
  --               ]
  --               0)
  radialArr <-
    if flag
      then R.map magnitude . getNormalizedHistogramArr <$> decodeFile histPath
      else do
        putStrLn "Couldn't find a Green's function data. Start simulation..."
        printCurrentTime
        solveMonteCarloR2Z2T0S0Radial
          numThread
          numTrail
          maxTrail
          numPoint
          numPoint
          thetaSigma
          scaleSigma
          maxScale
          tao
          theta0Freqs
          thetaFreqs
          scale0Freqs
          scaleFreqs
          ""
          (emptyHistogram
             [ (round . sqrt . fromIntegral $ 2 * (div numPoint 2) ^ 2)
             , L.length scale0Freqs
             , L.length theta0Freqs
             , L.length scaleFreqs
             , L.length thetaFreqs
             ]
             0)
  printCurrentTime
  arrR2Z2T0S0 <-
    computeUnboxedP $
    computeR2Z2T0S0ArrayRadial
      radialArr
      numPoint
      numPoint
      1
      maxScale
      thetaFreqs
      scaleFreqs
      theta0Freqs
      scale0Freqs
  arrR2Z2 <- R.sumP . R.sumS . rotateR2Z2T0S0Array $ arrR2Z2T0S0
      -- arr =
      --   r2z2t0s0Tor2s1rpt0s0
      --     numOrientation
      --     thetaFreqs
      --     numScale
      --     scaleFreqs
      --     arrR2Z2T0S0
      -- arr4d =
      --   R.slice arr $
      --   (Z :. All :. All :. (L.length theta0Freqs - 1) :.
      --    (L.length scale0Freqs - 1) :.
      --    All :.
      --    All)
  let arr4d = r2z2Tor2s1rp numOrientation thetaFreqs numScale scaleFreqs arrR2Z2
  createDirectoryIfMissing True folderPath
  printCurrentTime
  MP.mapM_
    (\(i, j) ->
       plotImageRepaComplex
         (folderPath </> show (i + 1) L.++ "_" L.++ show (j + 1) L.++ ".png") .
       ImageRepa 8 .
       computeS . R.extend (Z :. (1 :: Int) :. All :. All) . R.slice arr4d $
       (Z :. i :. j :. All :. All))
    [(i, j) | i <- [0 .. numOrientation - 1], j <- [0 .. numScale - 1]]
  arr2D <- R.sumP . R.sumS . rotate4D . rotate4D $ arr4d
  plotImageRepaComplex (folderPath </> "sum.png") .
    ImageRepa 8 . computeS . R.extend (Z :. (1 :: Int) :. All :. All) $
    arr2D

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-|
Module: MachineLearning.Types
Description: Common Types
Copyright: (c) Alexander Ignatyev, 2016
License: BSD-3
Stability: experimental
Portability: POSIX

Common type definitions used in all modules.
-}

module MachineLearning.Types
(
    R
  , Vector
  , Matrix
)

where

import qualified Numeric.LinearAlgebra.Data as LAD

-- | Scalar Type (Double)
type R = LAD.R

-- | Vector Types (of Doubles)
type Vector = LAD.Vector R

-- | Matrix Types (of Doubles)
type Matrix = LAD.Matrix R

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
-----------------------------------------------------------------------------
-- |
-- Module     : Numeric.LinearAlgebra
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- Linear algebra types and operations
--

module Numeric.LinearAlgebra (
    -- * Vector types
    Vector,
    RVector,
    STVector,
    IOVector,

    -- * Matrix types
    Matrix,
    RMatrix,
    STMatrix,
    IOMatrix,

    -- * Packed matrix types
    Packed,
    RPacked,
    STPacked,
    IOPacked,

    module Numeric.LinearAlgebra.Types,
    ) where

import Numeric.LinearAlgebra.Types
import Numeric.LinearAlgebra.Vector( Vector, RVector, STVector, IOVector )
import Numeric.LinearAlgebra.Matrix( Matrix, RMatrix, STMatrix, IOMatrix )
import Numeric.LinearAlgebra.Packed( Packed, RPacked, STPacked, IOPacked )

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
module PCA
    ( fitPCA, transformPCA
    ) where

import PCA.QModels
import Numeric.LinearAlgebra.Data
       (Matrix, Vector, toRows, fromRows)
import qualified Numeric.LinearAlgebra.HMatrix as H

-- | Theta parameteres which we will you for predict in model
data PCAModel = PCAModel
    { scale :: Matrix Double  -- ^ W
    , shift :: Vector Double  -- ^ \mu
    }

fitPCA :: Matrix Double -> PCAModel
fitPCA t =
    let (w,mu) = calculateQ t
    in PCAModel
       { scale = w
       , shift = mu
       }

transformPCA :: PCAModel -> Matrix Double -> Matrix Double
transformPCA (PCAModel w mu) t =
    fromRows
        (map
             (\tn ->
                   tn + mu)
             (toRows (t H.<> w)))

max_stars_repo_name: DbIHbKA/vbpca

max_issues_repo_name: DbIHbKA/vbpca

max_forks_repo_name: DbIHbKA/vbpca


text:
module Raytracer.Geometry where

import Numeric.LinearAlgebra.Data (Vector, (|>), toList, (!), fromList, toList)
import Numeric.LinearAlgebra.Data (fromColumns, asColumn, toColumns)
import Data.Maybe.HT (toMaybe)
import Data.Monoid (First(First), Monoid, getFirst, mempty, mappend, mconcat)
import Numeric.LinearAlgebra (luSolve, luPacked, pinv, (<>))
import Control.Applicative ((<*>), pure)
import Control.Monad (join)
import Data.Foldable (find)
import Data.List (sort)
import Data.Function (on)
import Data.Maybe (isJust)
import Data.Bool.HT (if')

(|*) = (*).(3|>).repeat

all_of :: [(a -> Bool)] -> a -> Bool
all_of conditions value = and $ conditions <*> (pure value)

type Pos = Vector Double
type Dir = Vector Double

-- This represents a ray of the form r = ax + b
data Ray = Ray Dir Pos deriving (Show)
calc_ray (Ray a b) n1 = round_to $ (n1 |* a) + b
  where
  round_to = fromList . map (epsilon_round 6) . toList
  epsilon_round n f = (fromInteger $ round $ f * (10^n)) / (10.0^^n)

-- This computes the ray that goes from the first point to the second
rayTo :: Pos -> Pos -> Ray
p1 `rayTo` p2 = Ray (p2 - p1) p1

-- This wraps a collision with a ray and an object
-- The `a` is what the object was holding, and the double is how far down the
-- ray the collision occured (used to find the "first" collision)
data Collision a = Collision Double a deriving (Show)

dist (Collision d _) = d

instance Eq (Collision a) where
  (==) = (==) `on` dist

instance Ord (Collision a) where
  compare = compare `on` dist

class Intersectable a where
	intersection :: Ray -> a b -> Maybe (Collision b)

ray_plane_intersect :: Ray -> Plane a -> (Vector Double)
ray_plane_intersect (Ray r1 r2) (Plane a p1 p2 p3) = head.toColumns.(\m -> (luSolve.luPacked) m (asColumn (r2 - p3))) $ fromColumns [p1, p2, -1 |* r1]

-- This represents a plane of the form p = ax + by + c
data Plane a = Plane a Dir Dir Pos deriving (Show)
calc_plane (Plane _ a b c) n1 n2 = (n1 |* a) + (n2 |* b) + c

instance Intersectable Plane where
	intersection ray plane@(Plane a _ _ _) = find (all_of [not.isInfinite, (>0)] . dist) $ Just $ Collision n a
		where
		n = (ray_plane_intersect ray plane) ! 2

-- Each of these Vectors represent a point in 3-space
data Triangle a = Triangle a Pos Pos Pos deriving (Show)

instance Intersectable Triangle where
	intersection ray (Triangle a p1 p2 p3) = if' (not $ all_of [all (>0), (1>).sum] [n1,n2]) Nothing $ find (all_of [not.isInfinite, (>0)] . dist) $ Just $ Collision n a
		where
		plane = Plane a (p2 - p1) (p3 - p1) p1
		(n1:n2:n:_) = toList $ ray_plane_intersect ray plane

-- We'll call a collection of triangles a mesh...
data Mesh a = Mesh [Triangle a] deriving (Show)

-- Most of this is just proxying the list of triangles.
instance Monoid (Mesh a) where
	mempty = Mesh []
	mappend (Mesh l1) (Mesh l2) = Mesh $ mappend l1 l2

instance Intersectable Mesh where
	-- We'll return, for the mesh, the first one we hit
	intersection ray (Mesh triangles) = join $ find isJust $ sort $ fmap (intersection ray) triangles

-- This function takes a point and two vectors and generates a mesh representing that square
square a v1 v2 p1 = Mesh [Triangle a p1 (p1 + v1) (p1 + v1 + v2), Triangle a (p1 + v1 + v2) (p1 + v2) p1]

-- This function takes a point and three vectors and generates a mesh representing that cube
cube a v1 v2 v3 p1 = mconcat $ fmap (uncurry3 $ square a) [(v1, v2, p1), (v2, v3, p1), (v3, v1, p1), (neg v1, neg v2, p2), (neg v2, neg v3, p2), (neg v3, neg v1, p2)]
	where
	p2 = p1 + v1 + v2 + v3
	neg = ((-1) |*)
	uncurry3 func (a, b, c) = func a b c


text:
module Fourier.FFT where

import Control.Parallel.Strategies
import Data.Complex(Complex,Complex((:+)),realPart,imagPart)
import Data.Matrix(matrix,fromList)

-- Cooley & Tukey Fast Fourier transformation of a sequence
fft :: [Double] -> [Complex Double]
fft xs = dft xs

-- Typical Discrete Fourier transformation of a sequence
dft :: [Double] -> [Complex Double]
dft [] = []
dft xs = do { k <- [0..length xs - 1]
  ; let re = xk k xs cos
  ; let im = -1 * (xk k xs sin)
  ; return (re :+ im)}

-- Calculate the sum of coefficience at [k]
-- Formula: xk = sum x[n] f(2pi * k * n/N)
xk :: Int -> [Double] -> (Double -> Double) -> Double
xk k xs f = let ys = do { (i,x) <- zip [0..] xs
  ; let n = fromIntegral $ length xs
  ; return (x * f( 2.0 * pi * (fromIntegral k) * (fromIntegral i) / n)) }
    in sum ys


text:
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications    #-}
module ArrayFire.GraphicsSpec where

import Control.Exception
import Data.Complex
import Data.Word
import Foreign.C.Types
import GHC.Int
import Test.Hspec

import ArrayFire

spec :: Spec
spec =
  describe "Graphics tests" $ do
    it "Should create window" $ do
      (1 + 1) `shouldBe` 2

max_stars_repo_name: danielkroeni/arrayfire-haskell

max_issues_repo_name: danielkroeni/arrayfire-haskell

max_forks_repo_name: danielkroeni/arrayfire-haskell


text:
{-# LANGUAGE GADTs, RankNTypes, MultiParamTypeClasses, FlexibleInstances, UndecidableInstances, KindSignatures, ScopedTypeVariables, ConstraintKinds, TemplateHaskell, GeneralizedNewtypeDeriving, TypeFamilies #-}

module Numeric.Trainee.Types (
	Gradee(..),
	Ow(..),
	params, forwardPass, Learnee(..),
	Loss, Regularization, HasNorm(..),
	Sample(..), (⇢),
	Samples, samples,

	module Numeric.Trainee.Params
	) where

import Prelude hiding (id, (.))
import Prelude.Unicode

import Control.Category
import Control.DeepSeq
import Control.Lens
import qualified Data.Vector as V
import Numeric.LinearAlgebra (Normed(norm_1), R, Vector)

import Numeric.Trainee.Params

newtype Gradee a b = Gradee {
	runGradee ∷ Lens' a b }

instance Category Gradee where
	id = Gradee id
	Gradee f . Gradee g = Gradee (g . f)

class Category a ⇒ Ow a where
	first ∷ a b c → a (b, d) (c, d)
	first = flip stars id

	second ∷ a b c → a (d, b) (d, c)
	second = stars id

	stars ∷ a b c → a b' c' → a (b, b') (c, c')
	vstars ∷ a b c → a (V.Vector b) (V.Vector c)

instance Ow Gradee where
	stars (Gradee f) (Gradee g) = Gradee $ lens g' s' where
		g' (x, y) = (view f x, view g y)
		s' (x, y) (x', y') = (set f x' x, set g y' y)

	vstars (Gradee f) = Gradee $ lens g' s' where
		g' = V.map (view f)
		s' = V.zipWith (flip (set f))

data Learnee a b = Learnee {
	_params ∷ Params,
	_forwardPass ∷ Params → a → (b, b → (Params, a)) }

makeLenses ''Learnee

instance NFData (Learnee a b) where
	rnf (Learnee ws _) = rnf ws

instance Show (Learnee a b) where
	show (Learnee ws _) = show ws

instance Category Learnee where
	id = Learnee (Params NoParams) fn where
		fn _ x = (x, const (Params NoParams, x))
	Learnee rws g . Learnee lws f = lws `deepseq` rws `deepseq` Learnee (Params (lws, rws)) (h ∘ castParams) where
		h (lws', rws') x = x `seq` y `seq` lws' `deepseq` rws' `deepseq` (z, up) where
			(y, f') = f lws' x
			(z, g') = g rws' y
			up dz = dz `seq` dy `seq` lws'' `deepseq` rws'' `deepseq` (Params (lws'', rws''), dx) where
				(rws'', dy) = g' dz
				(lws'', dx) = f' dy

type Loss b = b → b → (b, b)

type Regularization b = b → (b, b)

class HasNorm a where
	type Norm a
	norm ∷ a → Norm a

instance HasNorm Float where
	type Norm Float = Float
	norm = id

instance HasNorm Double where
	type Norm Double = Double
	norm = id

instance Normed (Vector a) ⇒ HasNorm (Vector a) where
	type Norm (Vector a) = R
	norm = norm_1

data Sample a b = Sample {
	sampleInput ∷ a,
	sampleOutput ∷ b }
		deriving (Eq, Ord)

(⇢) ∷ a → b → Sample a b
(⇢) = Sample

instance Bifunctor Sample where
	bimap f g (Sample x y) = Sample (f x) (g y)

instance (Show a, Show b) ⇒ Show (Sample a b) where
	show (Sample xs ys) = show xs ++ " => " ++ show ys

instance (NFData a, NFData b) ⇒ NFData (Sample a b) where
	rnf (Sample i o) = rnf i `seq` rnf o

type Samples a b = V.Vector (Sample a b)

samples ∷ [Sample a b] → Samples a b
samples = V.fromList

max_stars_repo_name: mvoidex/trainee

max_issues_repo_name: mvoidex/trainee

max_forks_repo_name: mvoidex/trainee


text:
module Dungeon
  ( Level
  , emptyLvl
  , execLevelGen
  , rndLvl
  , getCell
  , findRandomCell
  , allCells
  , allCellsBetween
  , printRndLvl
  ) where

import           Control.Monad.Random
import           Control.Monad.Random.Class
import           Control.Monad.State.Strict
import           Data.List
import qualified Math.Geometry.Grid           as Grid
import           Math.Geometry.Grid.Octagonal (RectOctGrid, rectOctGrid)
import qualified Math.Geometry.GridMap        as GridMap
import           Math.Geometry.GridMap.Lazy   (LGridMap, lazyGridMap)
import           Safe                         (headMay)
import           Numeric.Vector

data Object = At | Enemy | Loot deriving (Show, Eq)
type Cell = Maybe [Object]
type Level = LGridMap RectOctGrid Cell
type LevelGen g a = RandT g (State Level) a

data Rect
  = Rect
  { left   :: Int
  , top    :: Int
  , width  :: Int
  , height :: Int
  }

getLeft :: Rect -> Int
getLeft = left

getTop :: Rect -> Int
getTop = top

getCell_ :: Level -> (Int, Int) -> Cell
getCell_ lvl idx =
  if lvl `Grid.contains` idx
  then lvl GridMap.! idx
  else Nothing

getCell :: Level -> Vec2i -> Cell
getCell l (Vec2 x y) = getCell_ l (x, y)

findRandomCell :: MonadRandom m => Level -> (Cell -> Bool) -> m (Vec2i, Cell)
findRandomCell lvl approve = do
  let (w, h) = Grid.size lvl
      go = do
        x <- getRandomR (0, w-1)
        y <- getRandomR (0, h-1)
        let cell = getCell_ lvl (x, y)
        if approve cell
          then return (Vec2 x y, cell)
          else go
  go

allCells :: Level -> [(Vec2i, Cell)]
allCells lvl = flip map (GridMap.toList lvl) $ \((x, y), c) ->
  (Vec2 x y, c)

allCellsBetween :: Level -> Vec2i -> Vec2i -> [(Vec2i, Cell)]
allCellsBetween lvl (Vec2 left top) (Vec2 right bottom) =
  let idxs = [Vec2 x y | y <- [top..bottom], x <- [left..right]]
  in flip map idxs $ \pos@(Vec2 x y) ->
    (pos, join $ GridMap.lookup (x, y) lvl)

emptyLvl :: Int -> Int -> Level
emptyLvl w h = lazyGridMap (rectOctGrid h w) (replicate (w*h) Nothing)

execLevelGen :: Int -> Int -> LevelGen StdGen a -> IO Level
execLevelGen w h gen = do
  stc <- evalRandT gen <$> newStdGen
  return $ execState stc (emptyLvl w h)

center :: Rect -> (Double, Double)
center Rect { .. } =
  (fromIntegral left + fromIntegral width / 2,
    fromIntegral top + fromIntegral height / 2)

anyCell :: MonadRandom m => Rect -> m (Int, Int)
anyCell Rect { .. } = do
  x <- getRandomR (left, left + width - 1)
  y <- getRandomR (top, top + height - 1)
  return (x, y)

dist :: Rect -> Rect -> Double
dist ra rb =
  let (xa, ya) = center ra
      (xb, yb) = center rb
  in sqrt ((xa - xb)**2 + (ya - yb)**2)

lvlTxt :: Level -> String
lvlTxt lvl = unlines rows where
  (w, h) = Grid.size lvl
  rows = map renderLine [0 .. h-1]
  renderLine y = intersperse ' ' $ map (renderChar y) [0 .. w-1]
  renderChar y x =
    case lvl GridMap.! (x, y) of
      Nothing -> '.'
      Just l ->
        if length l < 10
        then last $ show (length l)
        else '*'

eatCells :: [(Int, Int)] -> LevelGen g ()
eatCells cells = do
  lvl <- get
  put $ foldl' (\l cell -> GridMap.insert cell (Just []) l) lvl cells

rndLvl :: RandomGen g => Int -> Int -> Double -> LevelGen g [Rect]
rndLvl minRoomEdge maxRoomEdge_ touchThresh = do
  (w, h) <- Grid.size <$> get
  descend (Rect 1 1 (w-2) (h-2))
    where
      -- something bigger than maxRoomEdge needs to be divisible without violating
      -- minRoomEdge
      maxRoomEdge = max maxRoomEdge_ (2 * minRoomEdge - 1)
      -- implicit bsp tree by recursion
      descend rect@Rect{..} = do
        let needVert = width > maxRoomEdge
            needHoriz = height > maxRoomEdge
        if not needVert && not needHoriz
          then do
            allowTouchChance :: Double <- getRandomR (0, 1)
            let minRoomWidth = max (width `div` 2) minRoomEdge
                minRoomHeight = max (height `div` 2) minRoomEdge
                space = if allowTouchChance <= touchThresh then 0 else 1
            if width - space < minRoomWidth || height - space < minRoomHeight
              then return []
              else do
                roomWidth <- getRandomR (minRoomWidth, width - space)
                roomHeight <- getRandomR (minRoomHeight, height - space)
                widthOffset <- getRandomR (0, width - roomWidth - space)
                heightOffset <- getRandomR (0, height - roomHeight - space)
                let roomTop = top + heightOffset + space
                    roomLeft = left + widthOffset + space
                eatCells [(roomLeft + x, roomTop + y)
                         | x <- [0 .. roomWidth-1], y <- [0 .. roomHeight-1]]
                return [Rect roomLeft roomTop roomWidth roomHeight]
          else do
            doVertical <- getRandom
            if not needHoriz || needVert && doVertical
              then do
                leftWidth <- getRandomR (minRoomEdge, width - minRoomEdge)
                rectsLeft <- descend rect { width = leftWidth }
                rectsRight <- descend rect { left = left + leftWidth, width = width - leftWidth }
                sequence_ $ do
                  rectA <- headMay $ sortOn getLeft rectsRight
                  rectB <- headMay $ sortOn (dist rectA) rectsLeft
                  return $ makeCorridor rectA rectB
                return (rectsLeft ++ rectsRight)
              else do
                topHeight <- getRandomR (minRoomEdge, height - minRoomEdge)
                rectsTop <- descend rect { height = topHeight }
                rectsBottom <- descend rect { top = top + topHeight, height = height - topHeight }
                sequence_ $ do
                  rectA <- headMay $ sortOn getTop rectsBottom
                  rectB <- headMay $ sortOn (dist rectA) rectsTop
                  return $ makeCorridor rectA rectB
                return (rectsTop ++ rectsBottom)

range :: Int -> Int -> [Int]
range a b =
  if a <= b
  then [a .. b]
  else [b .. a]

makeCorridor :: RandomGen g => Rect -> Rect -> LevelGen g ()
makeCorridor a b = do
  (xa, ya) <- anyCell a
  (xb, yb) <- anyCell b
  let horiz y = eatCells [(x, y) | x <- range xa xb]
      vert x = eatCells [(x, y) | y <- range ya yb]
  horizFirst <- getRandom
  if horizFirst
    then do
      horiz ya
      vert xb
    else do
      vert xa
      horiz yb

-- printRndLvl 100 100 3 20 0.7
printRndLvl :: Int -> Int -> Int -> Int -> Double -> IO ()
printRndLvl w h minRoomEdge maxRoomEdge touchThresh = do
  lvl <- execLevelGen w h (rndLvl minRoomEdge maxRoomEdge touchThresh)
  putStrLn $ lvlTxt lvl

max_stars_repo_name: cjay/some-roguelike

max_issues_repo_name: cjay/some-roguelike

max_forks_repo_name: cjay/some-roguelike


text:
module EquivarianceDiagram where

import           Data.Array.Repa     as R
import           Data.Complex
import           Data.List           as L
import           Data.Vector.Unboxed as VU
import           DFT.Plan
import           Filter.Gaussian
import           Image.IO
import           System.Directory
import           System.Environment
import           System.FilePath

main = do
  (inputPath:sigmaStr:_) <- getArgs
  (ImageRepa _ img) <- readImageRepa inputPath True
  let (Z :. channels :. cols :. rows) = extent img
      sigma = read sigmaStr :: Double
      folderPath = "output/test/EquivarianceDiagram"
  createDirectoryIfMissing True folderPath
  lock <- getFFTWLock
  (plan1, imgF) <-
    dft1dGPlan lock emptyPlan [channels, cols, rows] [1, 2] .
    VU.convert . VU.map (:+ 0) . toUnboxed $
    img
  (plan2, _) <- idft1dGPlan lock plan1 [channels, cols, rows] [1, 2] imgF
  (plan, gaussianFilterF) <-
    gaussian2DFilter plan2 (Gaussian2DParams sigma rows cols)
  convolvedImg <- convolveGaussian2D plan gaussianFilterF . R.map (:+ 0) $ img
  plotImageRepa (folderPath </> (takeBaseName inputPath) L.++ "_out.png") .
    ImageRepa 8 . computeS . R.map magnitude $
    convolvedImg