text:
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Lib.Utility.NeuralNet where

import qualified Numeric.LinearAlgebra as LA
import Numeric.LinearAlgebra.Static
import Numeric.LinearAlgebra.Devel
import Control.Monad.State
import Control.Arrow
import GHC.TypeLits
import Data.Proxy

data Layer (il :: Nat) (is :: Nat) (wc :: Nat) (state :: *) where
    Layer :: state -> (X il is -> (Y il wc, state)) -> (state -> Delta il wc -> (DX il is, state)) -> Layer il is wc state

type ReluState il wc = Maybe (L il wc)
reluLayer :: forall il wc. (KnownNat il, KnownNat wc) =>
   Layer il wc wc (ReluState il wc)
reluLayer = Layer Nothing forwardF backwardF
    where forwardF x = (lt0, Just lt0)
              where (Just lt0) = create $ LA.cmap (\ e -> if e <= 0 then 0 else e) $ extract x
          backwardF Nothing _ = error "State is Nothing"
          backwardF (Just lt0) d = (zipWithM' (\ e1 e2 -> if e1 == 0 then 0 else e2) lt0 d, Just lt0)

type SigmoidState il wc = Maybe (L il wc)
sigmoidLayer :: forall il wc. (KnownNat il, KnownNat wc) =>
   Layer il wc wc (ReluState il wc)
sigmoidLayer = Layer Nothing forwardF backwardF
    where forwardF x = (y, Just y)
              where y = 1 / (1 + exp (-x))
          backwardF Nothing _ = error "State Y is Nothing"
          backwardF (Just s) d = (d * (1 - s) * s, Just s)

type AffineState il is wc = (NeuralParam is wc, Maybe (X il is), Maybe (NeuralParam is wc))
affineLayer :: forall il is wc. (KnownNat il, KnownNat is, KnownNat wc) =>
    NeuralParam is wc -> Layer il is wc (AffineState il is wc)
affineLayer param@(w, b) = Layer (param, Nothing, Nothing) forwardF backwardF
    where forwardF x = (x <> w + v2m b, (param, Just x, Nothing))
          backwardF (_, Nothing, _) _ = error "State X is Nothing"
          backwardF (param@(w, b), Just x, _) d = (dx, (param, Just x, Just (dw, db)))
              where dx = d <> tr w
                    dw = tr x <> d
                    db = unrow $ (1 :: L 1 il) <> d

type SoftmaxWithLossState il os = (T il os, Maybe (Y il os), Maybe (Loss il 1))
softmaxWithLossLayer :: forall il os. (KnownNat il, KnownNat os) =>
    T il os -> Layer il os 1 (SoftmaxWithLossState il os)
softmaxWithLossLayer t = Layer (t, Nothing, Nothing) forwardF backwardF
    where forwardF x = (loss, (t, Just y, Just loss))
              where y = softmax x
                    loss = crossEntropyError t y
          backwardF (_, Nothing, _) d = error "State Y is Nothing"
          backwardF state@(t, Just y, _) d = ((y - t) / fromIntegral il, state)
              where (il, _) = size d

class NeuralNet (a :: Nat -> Nat -> Nat -> *) where
    predict :: forall il is os. (KnownNat il, KnownNat is, KnownNat os)
         => a il is os -> X il is -> Prediction il os
    loss :: forall il is os. (KnownNat il, KnownNat is, KnownNat os)
         => a il is os -> X il is -> Loss il os
    accuracy :: forall il is os. (KnownNat il, KnownNat is, KnownNat os)
         => a il is os -> X il is -> T il os -> Double
    gradients :: forall il is os. (KnownNat il, KnownNat is, KnownNat os)
         => a il is os -> X il is -> T il os -> ()
    training :: forall il is os. (KnownNat il, KnownNat is, KnownNat os)
         => a il is os -> X il is -> T il os -> a il is os

    accuracy a x t = LA.sumElements (extract t1f0) / fromIntegral (size t1f0)
        where y = predict a x
              (y', t') = (maxIndex' y, maxIndex' t)
              t1f0 = uncol $ zipWithM' (\ e1 e2 -> if e1 == e2 then 1 else 0) y' t'

data MnistTwoLayerNet (wc :: Nat) (il :: Nat)  (is :: Nat) (os :: Nat) where
    MnistTwoLayerNet ::
        Layer il is wc (AffineState il is wc)
     -> Layer il wc wc (ReluState il wc)
     -> Layer il wc os (AffineState il wc os)
     -> AdaGrad wc il is os
     -> MnistTwoLayerNet wc il is os
instance NeuralNet (MnistTwoLayerNet 50) where
    predict (MnistTwoLayerNet af1 relu af2 _) x = y3
        where (y1, af1s) = case af1 of (Layer _ ff _) -> ff x
              (y2, relus) = case relu of (Layer _ ff _) -> ff y1
              (y3, af2s) = case af2 of (Layer _ ff _) -> ff y2
    loss = undefined
    gradients = undefined
    training = undefined

class Updater (a :: Nat -> Nat -> Nat -> Nat -> *) where
    update :: forall wc il is os. (KnownNat wc, KnownNat il, KnownNat is, KnownNat os)
         => a wc il is os
         -> (NeuralParam is wc, NeuralParam wc os)
         -> (NeuralParam is wc, NeuralParam wc os)
         -> ((NeuralParam is wc, NeuralParam wc os), a wc il is os)

data SGD (wc :: Nat) (il :: Nat) (is :: Nat) (os :: Nat) where
    SGD :: Double -> SGD wc il is os
instance Updater SGD where
    update (SGD lr) (p1, p2) (dp1, dp2) = ((newParam lr p1 dp1, newParam lr p2 dp2), SGD lr)

newParam :: forall is wc. (KnownNat is, KnownNat wc) => Double
    -> NeuralParam is wc -> NeuralParam is wc -> NeuralParam is wc
newParam lr (w, b) (dw, db) = (w - lr_w * dw, b - lr_b * db)
    where (lr_w, lr_b) = (realToFrac lr, realToFrac lr)

data AdaGrad (wc :: Nat) (il :: Nat) (is :: Nat) (os :: Nat) where
    AdaGrad :: Double -> (NeuralParam is wc, NeuralParam wc os) -> AdaGrad wc il is os
instance Updater AdaGrad where
    update (AdaGrad lr (h1, h2)) (p1, p2) (dp1, dp2) = (result, AdaGrad lr newH)
        where result = (newParam lr p1 dp1', newParam lr p2 dp2')
              (dp1', dp2') = (f' newH1 dp1, f' newH2 dp2)
                  where f h dx = dx / sqrt (h + 1e-7)
                        f' (hw, hb) (dw, db) = (f hw dw, f hb db)
              newH@(newH1, newH2) = (f h1 dp1, f h2 dp2)
                  where f (hw, hb) (dw, db) = (dw ** 2 + hw, db ** 2 + hb)
newAdaGrad lr = AdaGrad lr ((0, 0), (0, 0))

initTestNet :: forall wc il is os. (KnownNat wc, KnownNat il, KnownNat is, KnownNat os) => Double -> IO (MnistTwoLayerNet wc il is os)
initTestNet wis = result
    where initParams = do
              w1 <- randn'
              let b1 = konst 0
              w2 <- randn'
              let b2 = konst 0
              return (w1 * realToFrac wis, b1, w2 * realToFrac wis, b2)
          result = do
              (w1, b1, w2, b2) <- initParams
              let affine1 = affineLayer (w1, b1)
              let relu = reluLayer
              let affine2 = affineLayer (w2, b2)
              return $ MnistTwoLayerNet affine1 relu affine2 $ newAdaGrad 0.1

type LastLayer il hs os = X il hs -> T il os

type X il is = L il is
type Y il wc = L il wc
type T il os = Y il os
type W is wc = L is wc
type B wc = R wc
type NeuralParam is wc = (W is wc, B wc)
type Prediction il wc = Y il wc
type Loss il wc = Y il wc
type Delta il is = X il is
type DX il is = X il is
type ForwardF il is wc = X il is -> Y il wc
type BackwardF il is wc = Y il wc -> X il is

v2m :: forall m n. (KnownNat m, KnownNat n) => R n -> L m n
v2m v = fromList $ (>>= id) $ replicate mVal $ LA.toList $ extract v
    where mVal = fromIntegral $ natVal (Proxy :: Proxy m)

softmax :: forall m n. (KnownNat m, KnownNat n) => L m n -> L m n
softmax a = expA / sumExpA
    where maxs = foldAndExtractRow' LA.maxElement a
          a' = a - maxs
          expA = exp a'
          sumExpA = foldAndExtractRow' LA.sumElements expA

crossEntropyError :: forall m n. (KnownNat m, KnownNat n) => L m n -> L m n -> L m 1
crossEntropyError t y = negate $ foldRow' LA.sumElements $ t * log (y + delta)
    where delta = 1e-7

foldRow' :: forall m n. (KnownNat m, KnownNat n) => (LA.Vector Double -> Double) -> L m n -> L m 1
foldRow' f m = fromList $ map (f . extract) $ toRows m

foldAndExtractRow' :: forall m n. (KnownNat m, KnownNat n) => (LA.Vector Double -> Double) -> L m n -> L m n
foldAndExtractRow' f m = fromList $ toRows m >>= (replicate nVal . f . extract)
    where nVal = fromIntegral $ natVal (Proxy :: Proxy n)

randn' :: forall m n. (KnownNat m, KnownNat n) => IO (L m n)
randn' = result
    where mVal = fromIntegral $ natVal (Proxy :: Proxy m)
          nVal = fromIntegral $ natVal (Proxy :: Proxy n)
          result = do
              m <- LA.randn mVal nVal
              return $ case create m of
                  Nothing -> error "Something wrong"
                  Just x -> x

dmmap' :: forall m n. (KnownNat m, KnownNat n) => (Double -> Double) -> L m n -> L m n
dmmap' f m = mapped
    where m' = extract m
          mapped = case create $ LA.cmap f m' of
              Nothing -> error "Something wrong"
              Just x -> x

zipWithM' :: forall m n. (KnownNat m, KnownNat n) => (Double -> Double -> Double) -> L m n -> L m n -> L m n
zipWithM' f m1 m2 = result
    where m1' = extract m1
          m2' = extract m2
          result = build (\ r c -> f (atIndex m1' r c) (atIndex m2' r c))
              where atIndex m r c = LA.atIndex m (truncate r, truncate c)

maxIndex' :: forall m n. (KnownNat m, KnownNat n) => L m n -> L m 1
maxIndex' = foldRow' (fromIntegral . LA.maxIndex)


text:
{-# LANGUAGE Rank2Types #-}
-----------------------------------------------------------------------------
-- |
-- Module     : Numeric.LinearAlgebra.Matrix.Herm
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- Hermitian views of matrices.
--

module Numeric.LinearAlgebra.Matrix.Herm (
    -- * Immutable interface

    -- ** Vector multiplication
    hermMulVector,
    hermMulVectorWithScale,
    addHermMulVectorWithScales,

    -- ** Matrix  multiplication
    hermMulMatrix,
    hermMulMatrixWithScale,
    addHermMulMatrixWithScales,

    -- ** Updates
    hermRank1Update,
    hermRank2Update,
    hermRankKUpdate,
    hermRank2KUpdate,

    -- * Mutable interface
    hermCreate,

    -- ** Vector multiplication
    hermMulVectorTo,
    hermMulVectorWithScaleTo,
    addHermMulVectorWithScalesM_,

    -- ** Matrix multiplication
    hermMulMatrixTo,
    hermMulMatrixWithScaleTo,
    addHermMulMatrixWithScalesM_,

    -- ** Updates
    hermRank1UpdateM_,
    hermRank2UpdateM_,
    hermRankKUpdateM_,
    hermRank2KUpdateM_,

    ) where

import Control.Monad( when )
import Control.Monad.ST( ST, runST, unsafeIOToST )
import Text.Printf( printf )

import Numeric.LinearAlgebra.Vector( Vector, RVector, STVector )
import qualified Numeric.LinearAlgebra.Vector as V
import Numeric.LinearAlgebra.Matrix.Base( Matrix )
import Numeric.LinearAlgebra.Matrix.STBase( STMatrix, RMatrix )
import qualified Numeric.LinearAlgebra.Matrix.STBase as M
import Numeric.LinearAlgebra.Types
import qualified Foreign.BLAS as BLAS

-- | A safe way to create and work with a mutable Herm Matrix before returning
-- an immutable one for later perusal.
hermCreate :: (Storable e)
           => (forall s. ST s (Herm (STMatrix s) e))
           -> Herm Matrix e
hermCreate mh = runST $ do
    (Herm u ma) <- mh
    a <- M.unsafeFreeze ma
    return $ Herm u a

-- | @hermRank1Update alpha x a@ returns
-- @alpha * x * x^H + a@.
hermRank1Update :: (BLAS2 e)
                => Double -> Vector e -> Herm Matrix e -> Herm Matrix e
hermRank1Update alpha x (Herm uplo a) = runST $ do
    ma' <- M.newCopy a
    hermRank1UpdateM_ alpha x (Herm uplo ma')
    a' <- M.unsafeFreeze ma'
    return $ Herm uplo a'

-- | @hermRank2Update alpha x y a@ returns
-- @alpha * x * y^H + conj(alpha) * y * x^H + a@.
hermRank2Update :: (BLAS2 e)
                => e -> Vector e -> Vector e -> Herm Matrix e
                -> Herm Matrix e
hermRank2Update alpha x y (Herm uplo a) = runST $ do
    ma' <- M.newCopy a
    hermRank2UpdateM_ alpha x y (Herm uplo ma')
    a' <- M.unsafeFreeze ma'
    return $ Herm uplo a'

-- | @hermRankKUpdate alpha trans a beta c@ returns
-- @c := alpha * a * a^H + beta * c@ when @trans@ is @NoTrans@ and
-- @c := alpha * a^H * a + beta * c@ when @trans@ is @ConjTrans@.  The
-- function signals an error when @trans@ is @Trans@.
hermRankKUpdate :: (BLAS3 e)
                => e -> Trans -> Matrix e -> e -> Herm Matrix e
                -> Herm Matrix e
hermRankKUpdate alpha trans a beta (Herm uplo c) = runST $ do
    mc' <- M.newCopy c
    hermRankKUpdateM_ alpha trans a beta (Herm uplo mc')
    c' <- M.unsafeFreeze mc'
    return $ Herm uplo c'

-- | @hermRank2KUpdate alpha trans a b beta c@ returns
-- @c := alpha * a * b^H + conj(alpha) * b * a^H + beta * c@ when @trans@ is
-- @NoTrans@ and @c := alpha * b^H * a + conj(alpha) * a^H * b + beta * c@
-- when @trans@ is @ConjTrans@.  The function signals an error when @trans@
-- is @Trans@.
hermRank2KUpdate :: (BLAS3 e)
                 => e -> Trans -> Matrix e -> Matrix e -> e -> Herm Matrix e
                 -> Herm Matrix e
hermRank2KUpdate alpha trans a b beta (Herm uplo c) = runST $ do
    mc' <- M.newCopy c
    hermRank2KUpdateM_ alpha trans a b beta (Herm uplo mc')
    c' <- M.unsafeFreeze mc'
    return $ Herm uplo c'

-- | @hermRank1UpdateM_ alpha x a@ sets
-- @a := alpha * x * x^H + a@.
hermRank1UpdateM_ :: (RVector v, BLAS2 e)
                  => Double -> v e -> Herm (STMatrix s) e -> ST s ()
hermRank1UpdateM_ alpha x (Herm uplo a) = do
    nx <- V.getDim x
    (ma,na) <- M.getDim a
    let n = nx

    when ((not . and) [ nx == n, (ma,na) == (n,n) ]) $ error $
        printf ("hermRank1UpdateM_ _ <vector with dim %d>"
                 ++ " (Herm _ <matrix with dim (%d,%d)>):"
                 ++ " invalid dimensions") nx ma na

    unsafeIOToST $
        V.unsafeWith x $ \px ->
        M.unsafeWith a $ \pa lda ->
            BLAS.her uplo n alpha px 1 pa lda

-- | @hermRank2UpdateM_ alpha x y a@ sets
-- @a := alpha * x * y^H + conj(alpha) * y * x^H + a@.
hermRank2UpdateM_ :: (RVector v1, RVector v2, BLAS2 e)
                  => e -> v1 e -> v2 e -> Herm (STMatrix s) e -> ST s ()
hermRank2UpdateM_ alpha x y (Herm uplo a) = do
    nx <- V.getDim x
    ny <- V.getDim y
    (ma,na) <- M.getDim a
    let n = nx

    when ((not . and) [ nx == n, ny == n, (ma,na) == (n,n) ]) $ error $
        printf ("hermRank2UpdateM_ _ <vector with dim %d>"
                 ++ " <vector with dim %d>"
                 ++ " (Herm _ <matrix with dim (%d,%d)>):"
                 ++ " invalid dimensions") nx ny ma na

    unsafeIOToST $
        V.unsafeWith x $ \px ->
        V.unsafeWith y $ \py ->
        M.unsafeWith a $ \pa lda ->
            BLAS.her2 uplo n alpha px 1 py 1 pa lda

-- | @hermRankKUpdateM_ alpha trans a beta c@ sets
-- @c := alpha * a * a^H + beta * c@ when @trans@ is @NoTrans@ and
-- @c := alpha * a^H * a + beta * c@ when @trans@ is @ConjTrans@.  The
-- function signals an error when @trans@ is @Trans@.
hermRankKUpdateM_ :: (RMatrix m, BLAS3 e)
                  => e -> Trans -> m e -> e -> Herm (STMatrix s) e
                  -> ST s ()
hermRankKUpdateM_ alpha trans a beta (Herm uplo c) = do
    (ma,na) <- M.getDim a
    (mc,nc) <- M.getDim c
    let (n,k) = if trans == NoTrans then (ma,na) else (na,ma)

    when (trans == Trans) $ error $
        printf ("hermRankKUpdateM_ _ %s:"
                 ++ " trans argument must be NoTrans or ConjTrans")
               (show trans)

    when ((not . and) [ (mc,nc) == (n,n)
                      , case trans of NoTrans -> (ma,na) == (n,k)
                                      _       -> (ma,na) == (k,n)
                      ]) $ error $
            printf ("hermRankKUpdateM_ _ %s <matrix with dim (%d,%d)> _"
                    ++ " (Herm _ <matrix with dim (%d,%d)>):"
                    ++ " invalid dimensions") (show trans) ma na mc nc

    unsafeIOToST $
        M.unsafeWith a $ \pa lda ->
        M.unsafeWith c $ \pc ldc ->
            BLAS.herk uplo trans n k alpha pa lda beta pc ldc

-- | @hermRank2KUpdateM_ alpha trans a b beta c@ sets
-- @c := alpha * a * b^H + conj(alpha) * b * a^H + beta * c@ when @trans@ is
-- @NoTrans@ and @c := alpha * b^H * a + conj(alpha) * a^H * b + beta * c@
-- when @trans@ is @ConjTrans@.  The function signals an error when @trans@
-- is @Trans@.
hermRank2KUpdateM_ :: (RMatrix m1, RMatrix m2, BLAS3 e)
                   => e -> Trans -> m1 e -> m2 e -> e -> Herm (STMatrix s) e
                   -> ST s ()
hermRank2KUpdateM_ alpha trans a b beta (Herm uplo c) = do
    (ma,na) <- M.getDim a
    (mb,nb) <- M.getDim b
    (mc,nc) <- M.getDim c
    let (n,k) = if trans == NoTrans then (ma,na) else (na,ma)

    when (trans == Trans) $ error $
        printf ("hermRank2KUpdateM_ _ %s:"
                 ++ " trans argument must be NoTrans or ConjTrans")
               (show trans)

    when ((not . and) [ (mc,nc) == (n,n)
                      , (mb,nb) == (ma,na)
                      , case trans of NoTrans -> (ma,na) == (n,k)
                                      _       -> (ma,na) == (k,n)
                      ]) $ error $
            printf ("hermRank2KUpdateM_ _ %s <matrix with dim (%d,%d)>"
                    ++ " <matrix with dim (%d,%d)> _"
                    ++ " (Herm _ <matrix with dim (%d,%d)>):"
                    ++ " invalid dimensions") (show trans) ma na mb nb mc nc

    unsafeIOToST $
        M.unsafeWith a $ \pa lda ->
        M.unsafeWith b $ \pb ldb ->
        M.unsafeWith c $ \pc ldc ->
            BLAS.her2k uplo trans n k alpha pa lda pb ldb beta pc ldc

-- | @hermMulVector a x@ returns @a * x@.
hermMulVector :: (BLAS2 e)
              => Herm Matrix e
              -> Vector e
              -> Vector e
hermMulVector a x =
    V.create $ do
        n <- V.getDim x
        y <- V.new_ n
        hermMulVectorTo y a x
        return y

-- | @hermMulVectorWithScale alpha a x@ retunrs @alpha * a * x@.
hermMulVectorWithScale :: (BLAS2 e)
                       => e
                       -> Herm Matrix e
                       -> Vector e
                       -> Vector e
hermMulVectorWithScale alpha a x =
    V.create $ do
        n <- V.getDim x
        y <- V.new_ n
        hermMulVectorWithScaleTo y alpha a x
        return y

-- | @addHermMulVectorWithScales alpha a x y@
-- returns @alpha * a * x + beta * y@.
addHermMulVectorWithScales :: (BLAS2 e)
                           => e
                           -> Herm Matrix e
                           -> Vector e
                           -> e
                           -> Vector e
                           -> Vector e
addHermMulVectorWithScales alpha a x beta y =
    V.create $ do
        y' <- V.newCopy y
        addHermMulVectorWithScalesM_ alpha a x beta y'
        return y'

-- | @hermMulMatrix side a b@
-- returns @alpha * a * b@ when @side@ is @LeftSide@ and
-- @alpha * b * a@ when @side@ is @RightSide@.
hermMulMatrix :: (BLAS3 e)
              => Side -> Herm Matrix e
              -> Matrix e
              -> Matrix e
hermMulMatrix side a b =
    M.create $ do
        mn <- M.getDim b
        c <- M.new_ mn
        hermMulMatrixTo c side a b
        return c

-- | @hermMulMatrixWithScale alpha side a b@
-- returns @alpha * a * b@ when @side@ is @LeftSide@ and
-- @alpha * b * a@ when @side@ is @RightSide@.
hermMulMatrixWithScale :: (BLAS3 e)
                       => e
                       -> Side -> Herm Matrix e
                       -> Matrix e
                       -> Matrix e
hermMulMatrixWithScale alpha side a b =
    M.create $ do
        mn <- M.getDim b
        c <- M.new_ mn
        hermMulMatrixWithScaleTo c alpha side a b
        return c

-- | @addHermMulMatrixWithScales alpha side a b beta c@
-- returns @alpha * a * b + beta * c@ when @side@ is @LeftSide@ and
-- @alpha * b * a + beta * c@ when @side@ is @RightSide@.
addHermMulMatrixWithScales :: (BLAS3 e)
                           => e
                           -> Side -> Herm Matrix e
                           -> Matrix e
                           -> e
                           -> Matrix e
                           -> Matrix e
addHermMulMatrixWithScales alpha side a b beta c =
    M.create $ do
        c' <- M.newCopy c
        addHermMulMatrixWithScalesM_ alpha side a b beta c'
        return c'

-- | @hermMulVectorTo dst a x@ sets @dst := a * x@.
hermMulVectorTo :: (RMatrix m, RVector v, BLAS2 e)
                => STVector s e
                -> Herm m e
                -> v e
                -> ST s ()
hermMulVectorTo dst = hermMulVectorWithScaleTo dst 1

-- | @hermMulVectorWithScaleTo dst alpha a x@
-- sets @dst := alpha * a * x@.
hermMulVectorWithScaleTo :: (RMatrix m, RVector v, BLAS2 e)
                         => STVector s e
                         -> e
                         -> Herm m e
                         -> v e
                         -> ST s ()
hermMulVectorWithScaleTo dst alpha a x =
    addHermMulVectorWithScalesM_ alpha a x 0 dst

-- | @addHermMulVectorWithScalesM_ alpha a x beta y@
-- sets @y := alpha * a * x + beta * y@.
addHermMulVectorWithScalesM_ :: (RMatrix m, RVector v, BLAS2 e)
                             => e
                             -> Herm m e
                             -> v e
                             -> e
                             -> STVector s e
                             -> ST s ()
addHermMulVectorWithScalesM_ alpha (Herm uplo a) x beta y = do
    (ma,na) <- M.getDim a
    nx <- V.getDim x
    ny <- V.getDim y
    let n = ny

    when (ma /= na) $ error $
        printf ("addHermMulVectorWithScalesM_ _"
                ++ " (Herm %s <matrix with dim (%d,%d)>)"
                ++ " %s <vector with dim %d>"
                ++ " _"
                ++ " <vector with dim %d>: Herm matrix is not square")
               (show uplo) ma na
               nx ny

    when ((not . and) [ (ma,na) == (n,n)
                      , nx == n
                      , ny == n
                      ]) $ error $
        printf ("addHermMulVectorWithScalesM_ _"
                ++ " (Herm %s <matrix with dim (%d,%d)>)"
                ++ " %s <vector with dim %d>"
                ++ " _"
                ++ " <vector with dim %d>: dimension mismatch")
               (show uplo) ma na
               nx ny

    unsafeIOToST $
        M.unsafeWith a $ \pa lda ->
        V.unsafeWith x $ \px ->
        V.unsafeWith y $ \py ->
            BLAS.hemv uplo n alpha pa lda px 1 beta py 1

-- | @hermMulMatrixTo dst side a b@
-- sets @dst := a * b@ when @side@ is @LeftSide@ and
-- @dst := b * a@ when @side@ is @RightSide@.
hermMulMatrixTo :: (RMatrix m1, RMatrix m2, BLAS3 e)
                => STMatrix s e
                -> Side -> Herm m1 e
                -> m2 e
                -> ST s ()
hermMulMatrixTo dst = hermMulMatrixWithScaleTo dst 1

-- | @hermMulMatrixWithScaleTo dst alpha side a b@
-- sets @dst := alpha * a * b@ when @side@ is @LeftSide@ and
-- @dst := alpha * b * a@ when @side@ is @RightSide@.
hermMulMatrixWithScaleTo :: (RMatrix m1, RMatrix m2, BLAS3 e)
                         => STMatrix s e
                         -> e
                         -> Side -> Herm m1 e
                         -> m2 e
                         -> ST s ()
hermMulMatrixWithScaleTo dst alpha side a b =
    addHermMulMatrixWithScalesM_ alpha side a b 0 dst

-- | @addHermMulMatrixWithScalesM_ alpha side a b beta c@
-- sets @c := alpha * a * b + beta * c@ when @side@ is @LeftSide@ and
-- @c := alpha * b * a + beta * c@ when @side@ is @RightSide@.
addHermMulMatrixWithScalesM_ :: (RMatrix m1, RMatrix m2, BLAS3 e)
                             => e
                             -> Side -> Herm m1 e
                             -> m2 e
                             -> e
                             -> STMatrix s e
                             -> ST s ()
addHermMulMatrixWithScalesM_ alpha side (Herm uplo a) b beta c = do
    (ma,na) <- M.getDim a
    (mb,nb) <- M.getDim b
    (mc,nc) <- M.getDim c
    let (m,n) = (mc,nc)

    when (ma /= na) $ error $
        printf ("addHermMulMatrixWithScalesM_ _"
                ++ " %s (Herm %s <matrix with dim (%d,%d)>)"
                ++ " <matrix with dim (%d,%d)>"
                ++ " _"
                ++ " <matrix with dim (%d,%d)>: Herm matrix is not square")
               (show side) (show uplo) ma na
               mb nb
               mc nc
    when ((not . and) [ case side of LeftSide  -> (ma,na) == (m,m)
                                     RightSide -> (ma,na) == (n,n)
                      , (mb, nb ) == (m,n)
                      , (mc, nc ) == (m,n)
                      ]) $ error $
        printf ("addHermMulMatrixWithScalesM_ _"
                ++ " %s (Herm %s <matrix with dim (%d,%d)>)"
                ++ " <matrix with dim (%d,%d)>"
                ++ " _"
                ++ " <matrix with dim (%d,%d)>: dimension mismatch")
               (show side) (show uplo) ma na
               mb nb
               mc nc

    unsafeIOToST $
        M.unsafeWith a $ \pa lda ->
        M.unsafeWith b $ \pb ldb ->
        M.unsafeWith c $ \pc ldc ->
            BLAS.hemm side uplo m n alpha pa lda pb ldb beta pc ldc

ext: hs

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
-- ch1 figure 1.1
-- https://stackoverflow.com/a/48212667

import Statistics.Sample

import qualified Statistics.Sample.Histogram as S
import qualified Data.Vector as V

vs :: V.Vector Double
vs = V.fromList xs

xs :: [Double]
xs = [2.1, 2.4, 2.2, 2.3, 2.7, 2.5, 2.4, 2.6, 2.6, 2.9]

hs :: V.Vector Int
hs = S.histogram_ bins lowerBound upperBound vs
  where
    bins = 5
    lowerBound = 2.05
    upperBound = 3.05

main :: IO ()
main = do
   putStrLn "Let us show the histgram, which corresponds to Figure1.1:"
   putStrLn $ "The date is:           " ++ show vs
   putStrLn $ "The number of smple is " ++ show (V.length vs)
   putStrLn $ "The histogram is:      " ++ show hs

   let m = mean vs
       s = stdDev vs
   putStrLn $ "The mean is: " ++ show m ++ " and the standard deviation is: " ++ show s

max_stars_repo_name: raydsameshima/Stat

max_issues_repo_name: raydsameshima/Stat

max_forks_repo_name: raydsameshima/Stat


text:
{-# LANGUAGE OverloadedStrings   #-}
{-# LANGUAGE ScopedTypeVariables #-}
-----------------------------------------------------------------------------
-- |
-- Module :
-- Copyright : (c) 2013 Boyun Tang
-- License : BSD-style
-- Maintainer : tangboyun@hotmail.com
-- Stability : experimental
-- Portability : ghc
--
--
--
-----------------------------------------------------------------------------
module Main where

import           Control.Arrow
import           Data.Colour.Names
import           Data.Colour.Palette.BrewerSet
import qualified Data.HashMap.Strict           as H
import qualified Data.List                     as L
import qualified Data.Text                     as T
import qualified Data.Text.IO                  as T
import qualified Data.Text.Read                as T
import qualified Data.Vector                   as V
import qualified Data.Vector.Unboxed           as UV
import           Diagrams.Backend.Cairo
import           Diagrams.HeatMap
import           Diagrams.HeatMap.Type
import           Diagrams.Prelude
import           Statistics.Quantile
import           Statistics.Sample
import           System.Environment
import           System.FilePath

cluDiffOpt :: ClustOpt
cluDiffOpt = ClustOpt
  { colorOpt = color2
  , rowCluster = Nothing -- Just (eucDis,UPGMA,LeftTree)
  , colCluster = Just (eucDis,UPGMA,BottomTree)}

color1 = let cSet = brewerSet RdBu 11
         in Three (cSet !! 9) white (cSet !! 1)

color2 = Three green black red

mkPara w h i j (a,b,c) = Para
  { clustOpt = cluDiffOpt
  , colorVal = ColorVal a b c
  , matrixHeight = h
  , matrixWidth = w
  , rowTreeHeight = 0.2 * min w h
  , colTreeHeight = 0.2 * min w h
  , rowFontSize = 0.8 * h / i
  , colFontSize = 0.6 * w / j
  , legendFontSize = 0.5 * w * 0.1
  , fontName = "Arial"
  , colorBarPos = Horizontal
  , tradeOff = Quality
  , colTreeLineWidth = 0.5
  , rowTreeLineWidth = 0.5
  }

eucDis :: (UV.Unbox a, Num a) => UV.Vector a -> UV.Vector a -> a
eucDis vec1 vec2 = UV.sum $ UV.zipWith (\v1 v2 -> (v1-v2)*(v1-v2)) vec1 vec2

toZscore :: UV.Vector Double -> UV.Vector Double
toZscore vec =
    let (m,v) = meanVarianceUnb vec
    in UV.map (\e -> (e - m) / sqrt v ) vec

parseTSV :: Bool -> FilePath -> FilePath -> IO Dataset
parseTSV doNormalization datFile labelFile = do
    sampleHash <- T.readFile labelFile >>=
                  return . H.fromList .
                  map (((!! 0) &&& (!! 1)) . T.split (== '\t')) .
                  T.lines .
                  T.filter (/= '\r')
    T.readFile datFile >>=
        return .
        (\(h:ts) ->
          let sIDs = V.fromList $ tail h
              grIDs = V.map (sampleHash `myIdx`) sIDs
              gIDs = V.fromList $ map head ts
              j = V.length sIDs
              i = V.length gIDs
              func = if doNormalization
                     then toZscore
                     else id
              datum = UV.concat $ map (func . UV.fromList . map (fst . fromRight . T.double) . tail) ts
              m = Matrix i j RowMajor datum
          in Dataset (Just gIDs) (Just sIDs) (Just grIDs) m
        ) . map (T.split (== '\t')) . T.lines . T.filter (/= '\r')
  where
    myIdx h k = if k `H.member` h
                then h H.! k
                else error $ show k
    fromRight :: Show a => Either a b -> b
    fromRight (Right b) = b
    fromRight (Left a) = error $ show a

main :: IO ()
main = do
   doNormal:w:h:datFile:labelFile:_ <- getArgs
   dataset <- parseTSV (read doNormal) datFile labelFile
   let vMin = continuousBy medianUnbiased 1 100 $ dat $ datM $ dataset
       vMean = if (read doNormal)
               then 0
               else continuousBy medianUnbiased 50 100 $ dat $ datM $ dataset
       vMax = continuousBy medianUnbiased 99 100 $ dat $ datM $ dataset
       j = fromIntegral $ nCol $ datM dataset
       i = fromIntegral $ nRow $ datM dataset
       para = mkPara (read w) (read h) i j (vMin,vMean,vMax)
   renderCairo (replaceExtension datFile "pdf") (mkWidth 600) $ pad 1.03 $ fst $ plotHeatMap para dataset

max_stars_repo_name: tangboyun/diagrams-heatmap

max_issues_repo_name: tangboyun/diagrams-heatmap

max_forks_repo_name: tangboyun/diagrams-heatmap


text:
{-# LANGUAGE OverloadedStrings #-}

import Criterion
import Criterion.Main (defaultConfig)
import Criterion.Types (anMean, reportAnalysis, timeLimit, anOutlierVar, ovEffect, OutlierEffect(Moderate))
import Statistics.Types (Estimate(..), ConfInt(..))

import Control.Monad (filterM)
import Control.Exception (catch, throwIO)
import Data.Foldable (for_)
import Data.Traversable (for)
import GHC.Exts (fromString)

import Language.Rust.Data.InputStream (InputStream)
import Language.Rust.Syntax (SourceFile)
import Language.Rust.Parser (readInputStream, Span, parse')

import System.Directory (getCurrentDirectory, listDirectory, createDirectoryIfMissing, doesFileExist, removeFile)
import System.FilePath ((</>), (<.>), takeFileName)
import System.Process (proc, readCreateProcess)
import System.IO.Error (isDoesNotExistError)

import Data.Aeson
import qualified Data.ByteString.Lazy as BL

main :: IO ()
main = do
  -- Open the output log file
  status <- readCreateProcess (proc "git" ["status", "--porcelain"]) ""
  logFileName <- case status of
                   "" -> init <$> readCreateProcess (proc "git" ["rev-parse", "HEAD"]) ""
                   _ -> pure "WIP"

  -- Get the test cases
  workingDirectory <- getCurrentDirectory
  let sampleSources = workingDirectory </> "sample-sources"
      benchIgnore = sampleSources </> ".benchignore"
  benchIgnoreExists <- doesFileExist benchIgnore
  ignore <- if benchIgnoreExists
              then (\f -> map (sampleSources </>) (lines f)) <$> readFile benchIgnore
              else pure []
  entries <- map (sampleSources </>) <$> listDirectory sampleSources
  files <- filterM doesFileExist (filter (`notElem` ignore) entries)

  -- Clear out previous WIP (if there is one)
  catch (removeFile (workingDirectory </> "bench" </> "timings" </> "WIP" <.> "json"))
        (\e -> if isDoesNotExistError e then pure () else throwIO e)

  -- Run 'criterion' tests
  reports <- for files $ \f -> do
    let name = takeFileName f
    putStrLn name
    is <- readInputStream f
    bnch <- benchmarkWith' defaultConfig{ timeLimit = 20 } (nf (parse' :: InputStream -> SourceFile Span) is)
    pure (name, bnch)
  let results = object [ fromString name .= object [ "mean" .= m
                                                   , "lower bound" .= l
                                                   , "upper bound" .= u
                                                   ]
                       | (name,report) <- reports
                       , let Estimate m (ConfInt l u _) = anMean (reportAnalysis report)
                       , ovEffect (anOutlierVar (reportAnalysis report)) < Moderate
                       ]
  for_ [ name | (name,report) <- reports, ovEffect (anOutlierVar (reportAnalysis report)) >= Moderate ] $ \n ->
    putStrLn $ "Benchmark for `" ++ n ++ "' will not be considered since it was inflated"

  -- Save the output to JSON
  createDirectoryIfMissing False (workingDirectory </> "bench" </> "timings")
  let logFile = workingDirectory </> "bench" </> "timings" </> logFileName <.> "json"
  putStrLn $ "writing results to: " ++ logFile
  logFile `BL.writeFile` encode results

max_stars_repo_name: cronburg/language-rust

max_issues_repo_name: cronburg/language-rust

max_forks_repo_name: cronburg/language-rust


text:
module FFT.DFTData(
    dftInnerMap,
    dftFoldMap,
    dftOuterMap,
    dftInnerMapReduce,
    dftInnerParMapReduce
) where
import Data.Complex
import FFT.Samples
import Control.Parallel
import Strategies

-- twiddle factors
tw :: Int -> Int -> Complex Float
tw n k = cis (-2 * pi * fromIntegral k / fromIntegral n)

dftInnerMap xs = [ sum (innerMap k) | k <- [0..n']]
  where
    innerMap k = parmap (\j -> xs!!j * tw n (j*k)) [0..n']
    n = length xs
    n' = n-1

dftFoldMap xs = [ (parfold (+) (0 :+ 0) (innerMap k)) | k <- [0..n']]
  where
    innerMap k = parmap (\j -> xs!!j * tw n (j*k)) [0..n']
    n = length xs
    n' = n-1

dftOuterMap xs = parmap innerFunc [0..n']
    where
        innerFunc k = sum [ xs!!j * tw n (j*k) | j <- [0..n']]
        n = length xs
        n' = n-1

dftInnerMapReduce xs = map (\k -> parMapReduceSimple rdeepseq (mapper k) rdeepseq reducer [0..n']) [0..n']
    where
        mapper k j = xs!!j * tw n (j*k)
        reducer = sum
        n = length xs
        n' = n-1

dftInnerParMapReduce xs = parmap (\k -> parMapReduceSimple rdeepseq (mapper k) rdeepseq reducer [0..n']) [0..n']
    where
        mapper k j = xs!!j * tw n (j*k)
        reducer = sum
        n = length xs
        n' = n-1

max_stars_repo_name: BNJHope/parallel-fourier-transforms

max_issues_repo_name: BNJHope/parallel-fourier-transforms

max_forks_repo_name: BNJHope/parallel-fourier-transforms


text:
import Data.Complex

main = do
  print $ 0 ^ 0
  print $ 0.0 ^ 0
  print $ 0 ^^ 0
  print $ 0 ** 0
  print $ (0 :+ 0) ^ 0
  print $ (0 :+ 0) ** (0 :+ 0)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE Arrows, FlexibleContexts #-}

module Examples.SingleTargetTracking where

import Control.Arrow (returnA)
import Control.Monad.Trans (liftIO, lift)
import Control.Monad.Bayes.Class (MonadSample)
import Control.Monad.Bayes.Sampler (SamplerIO, sampleIO)

import Data.Aeson (encode)
import Data.Maybe (fromJust)

import qualified Data.ByteString.Lazy.Char8 as BS (putStrLn)

import Numeric.LinearAlgebra.Static

import Inference (zdsparticles, zunheap, zparticles)
import DelayedSampling (DelayedInfer, Result (..))
import qualified SymbolicDistr as DS (sample, mvNormal)
import DSProg (DeepForce (..), Expr' (..), Expr, marginal, zdeepForce, deepForce', zdeepForce')
import Util.ZStream (ZStream)
import qualified Util.ZStream as ZS
import Util.Ref (MonadState, Heap)

import qualified Metaprob as MP

model :: MonadState Heap m => MonadSample m => Bool -> ZStream (MP.Gen m) () (Expr (R 6), Expr (R 3))
model delay = ZS.fromStep step initPosVel
  where
  initPosVel :: Expr (R 6)
  initPosVel = Const (konst 0 :: R 6)
  step posvel () = do
    posvel' <- lift (force (DS.sample (DS.mvNormal (MVMul (Const motionMatrix) posvel) motionCov)))
    observation <- "obs" MP.~~ (MP.dsPrim (DS.mvNormal (MVMul (Const posFromPosVel) posvel') (10 * sym eye)))
    return (posvel', (posvel', observation))

  force = if delay then id else (>>= deepForce')
  -- Constants
  posCov = 0.01 * sym eye
  velCov = 0.1 * sym eye
  tdiff = 1
  motionCov :: Sym 6
  motionCov =
    sym $ konst tdiff * ((unSym posCov ||| (konst 0 :: Sq 3))
          ===
          ((konst 0 :: Sq 3) ||| unSym velCov))
  motionMatrix :: Sq 6
  motionMatrix =
      ((eye :: Sq 3) ||| (konst tdiff * eye :: Sq 3))
                    ===
      ((konst 0 :: Sq 3) ||| (eye :: Sq 3))
  posFromPosVel :: L 3 6
  posFromPosVel = (eye :: Sq 3) ||| (konst 0 :: Sq 3)

processObservationStream :: DelayedInfer m => Bool -> ZStream m (R 3) (Result (R 6))
processObservationStream delay = proc observations -> do
  (posvel, _) <- MP.zobserving (model delay) -< ((), "obs" MP.|-> MP.obs observations)
  ZS.run -< fromJust <$> marginal posvel

runInference :: Bool -> Int -> ZStream SamplerIO () (R 6, [Result (R 6)], R 3)
runInference delay numParticles = proc () -> do
  (groundTruth, obs) <- simulate (model True) -< ()
  particles <- zdsparticles numParticles (processObservationStream delay) -< obs
  returnA -< (groundTruth, particles, obs)

runExample :: Bool -> Int -> IO ()
runExample delay n = sampleIO $ ZS.runStream (liftIO . BS.putStrLn . encode . f) (runInference delay n)
  where f (gt, particles, obs) = ([(0 :: Int, gt)], map (\x -> [(0 :: Int, x)]) particles, [obs])

simulate = zdeepForce . ZS.liftM MP.sim

max_stars_repo_name: psg-mit/probzelus-haskell

max_issues_repo_name: psg-mit/probzelus-haskell

max_forks_repo_name: psg-mit/probzelus-haskell


text:
{-# LANGUAGE Arrows #-}

import Control.Arrow
import Data.Complex

-- The DFT forms the specification for the FFT.
dft :: RealFloat a => [Complex a] -> [Complex a]
dft xs = [ sum [ x * exp((0.0 :+ (-1.0)) * 2 * pi / lenXs * fromIntegral(k * n))
                 | (x, n) <- zip xs [0..]
               ]
           | k <- [0..(length xs - 1)]
         ]
    where lenXs = fromIntegral $ length xs

-- FFT implementation
-- Note my reversal of Ross' even/odd convention, below.
-- I did this, to remain consistent with the explanation, above, which came from a pre-existing work.
fft :: RealFloat b => Hom (Complex b) (Complex b)
fft = id :&: proc (e, o) -> do
                e' <- fft             -< e
                o' <- fft >>> twiddle -< o
                unriffle -< f (e', o')
    where f (x, y) = (x + y, x - y)

twiddle :: RealFloat b => Hom (Complex b) (Complex b)
twiddle = id :&: twiddle' 1

twiddle' :: RealFloat b => Int -> Hom (Pair (Complex b)) (Pair (Complex b))
twiddle' n = (id *** (* phi)) :&: ((twiddle' n') *** ((twiddle' n') >>> (arr ((* phi) `prod` (* phi)))))
    where phi = cis (-pi / (fromIntegral n))
          n'  = n + 1

-- Testing machinery
-- Conversion to list, for comparison to specification (i.e. - DFT)
btToList :: BalTree b -> [b]
btToList (Zero x)      = [x]
btToList (Succ (x, y)) = btToList x ++ (btToList y)

-- Build a tree of depth 'n', taking values from supplied list, sequentially.
buildTree :: Int -> [b] -> BalTree b
buildTree = invertTree . buildTree'

data BalTree' b = Zero' b
                | Succ' (Pair (BalTree' b))

buildTree' :: Int -> [b] -> BalTree' b
buildTree' _       []     = error "Exhausted input!"
buildTree' 0       (x:xs) = Zero' x
buildTree' n       xs     = Succ' ((buildTree' (n-1) (evens xs)), (buildTree' (n-1) (odds xs)))

invertTree :: BalTree' b -> BalTree b
invertTree Zero' x            = Zero x
invertTree Succ' (bt'1, bt'2) = Succ (joinTrees (invertTree bt'1, invertTree bt'2))

-- joinTrees :: Pair (BalTree b) -> BalTree (Pair b)
joinTrees :: Pair (BalTree b) -> BalTree b
joinTrees (Zero x,   Zero y)   = Succ (Zero (x, y))
joinTrees (Succ bt1, Succ bt2) = Succ (joinTrees (bt1, bt2))


text:
{-# LANGUAGE PackageImports,BangPatterns,StandaloneDeriving,TypeFamilies,FlexibleContexts,FlexibleInstances,TypeOperators,ExistentialQuantification,NoMonomorphismRestriction,DeriveTraversable,DeriveFunctor,DeriveFoldable #-}
module Kalman where
import Space.Class hiding(dim)
import Solver.RungeKutta
import qualified Space.Class as M
import Data.Typeable
import Data.Distributive
import Data.Functor.Compose
import Linear.Matrix hiding(trace)
import qualified Linear.Metric as Metric
import Vectorization
import Control.Applicative
import Data.Maybe
import Data.Foldable (Foldable,foldl')
import qualified  Data.Foldable as F
import Data.Traversable (Traversable,fmapDefault,mapAccumL)
import qualified Data.Traversable as T
import Linear.V0
import Linear.V1
import Foreign.Storable
import Foreign.FStorable
import Foreign.Cholesky
import qualified Numeric.LinearAlgebra as M
import qualified Data.Packed.Matrix as M
import Local
import Linear.Vector
import CartesianProduct
import Control.Monad.Trans.State

-- Covariance Version

fullM
  :: (RealFloat a, FStorable (Local b),
      Vectorize (Local b), Vectorize (Local b1), Space b,
      Space b1, Cholesky a, M.Field a) =>
     (b a -> b a)
     -> Local b (Local b a)
     -> (b a -> b1 a)
     -> Local b1 (Local b1 a)
     -> b1 a
     -> (b a, Local b (Local b a))
     -> (b a, Local b (Local b a))
fullM f r h q  m s= execState  ( do
      modify (prediction f r)
      modify (measure h q m)
      --modify (smoothing f r s)
      --modify (prediction f r)
      )  s

full
  :: (RealFloat a, FStorable (Local b2), FStorable (Local b),
      Vectorize (Local b), Vectorize (Local b1), Space b2, Space b,
      Space b1, Cholesky a, M.Field a) =>
     (b2 a -> b a)
     -> Local b (Local b a)
     -> (b a -> b1 a)
     -> Local b1 (Local b1 a)
     -> b1 a
     -> (b2 a, Local b2 (Local b2 a))
     -> (b a, Local b (Local b a))
full f r h q  m s= sm
  where sp = prediction f r s
        sm = measure h q m sp
        ss = smoothing f r s sm
        ssp = prediction f r ss
        ssm = measure h q m ssp

predictionM f = modify . prediction f
measureM h q  =  modify . measure  h q
smoothingM f r  = modify . smoothing f r

prediction
  :: (M.Field a , Vectorize (Local b) ,RealFloat a,
      FStorable (Local b1), Space b1, Space b, Cholesky a) =>
     (b1 a -> b a)
     -> Local b (Local b a)
     -> (b1 a, Local b1 (Local b1 a))
     -> (b a, Local b (Local b a))
prediction f r (x ,p)= (mean ,pxx )
    where
        (mean,pxx,_,_) = sigmaSamplingCovGain f r (x,p)

measure
  :: (RealFloat a, FStorable (Local b),
      Vectorize (Local b1), Space b, Space b1, Cholesky a, M.Field a) =>
     (b a -> b1 a)
     -> Local b1 (Local b1 a)
     -> b1 a
     -> (b a, Local b (Local b a))
     -> (b a, Local b (Local b a))
measure h q m (x ,p) = ( xb  , pb )
    where
        (mmean,pxx,pxy,k) =sigmaSamplingCovGain h q (x,p)
        xb = x |+| ( k !*  ( m |-| mmean ) )
        pb = p !-! ( k !*! pxx !*! distribute k )

smoothing
  :: (RealFloat a, FStorable (Local b), Vectorize (Local b1),
      Space b, Space b1, Cholesky a,
      M.Field a) =>
     (b a -> b1 a)
     -> Local b1 (Local b1 a)
     -> (b a, Local b (Local b a))
     -> (b1 a, Local b1 (Local b1 a))
     -> (b a, Local b (Local b a))
smoothing f r (x,p) (x1,p1) = (xb ,pb)
    where
        (xmean ,pxx,pxy,k) = sigmaSamplingCovGain f r (x,p)
        xb = x |+| (k !* ( x1 |-| xmean ) )
        pb = p !+! (k !*! ( p1 !-! pxx ) !*! distribute k )

sigmaSamplingCovGain
  :: (RealFloat a,M.Field a ,Vectorize (Local b),
      FStorable (Local b1), Space b1, Space b, Cholesky a) =>
     (b1 a -> b a)
     ->  (Local b (Local b a))
     -> (b1 a, Local b1 (Local b1 a))
     -> (b a, Local b (Local b a), Local b (Local b1 a),Local b1 (Local b a ))
sigmaSamplingCovGain f r (x,p) = (mean , pxx,pxy,k)
    where
        spoints = (sigmaPoints1 x p)
        points = fmap f spoints
        mean = sigmaMean1 points
        pxx = sigmaCov mean points !+! r
        pxy = sigmaCovXY x mean spoints  points
        k =  distribute pxy  !*! inverseM pxx

sigmaSamplingCov
  :: (RealFloat a, Show (Local b1 (Local b1 a)),
      FStorable (Local b1), Space b1, Space b, Cholesky a) =>
     (b1 a -> b a)
     ->  (Local b (Local b a))
     -> (b1 a, Local b1 (Local b1 a))
     -> (b a, Local b (Local b a), Local b (Local b1 a))

sigmaSamplingCov f r (x,p) = (mean , pxx,pxy)
    where
        spoints = (sigmaPoints1 x p)
        points = fmap f spoints
        mean = sigmaMean1 points
        pxx = sigmaCov mean points !+! r
        pxy = sigmaCovXY x mean spoints  points

inverseM = fromLists .  M.toLists . M.inv  . M.fromLists . toLists

{-
measureDelayed h m (x,p) (xlag,plag) =  (xnew,pxk)
    where
        -- Xlag
        sxlag = sigmaPoints1 xlag plag
        -- X
        sx = sigmaPoints1 x p
        -- XXlag
        pxxlag = sigmaCovXY x xlag sx sxlag
        -- Xk
        xk = x :|: xlag
        pxk = prodSym p plag pxxlag
        sxk = sigmaPoints1 xk pxk
        -- Y
        (sy,ym,py) = sigmaTransform h xlag plag
        pxlagy = sigmaCovXY xk ym sxk sy
        k = (distribute pxlagy) !*! (inverseM py)
        xnew = xk |+| (k !* (m |-| ym) )
-}

prodSym x y xy = prod x y xy (distribute xy)

prod x y xy yx = (liftA2 (:|:) x xy ) :|: (liftA2 (:|:) yx  y)

-- Helpers

square :: (Num a, R f) => f a -> f (f a)
square v = mult v v

squareT v = (distribute v) !*! v

mult x m = fmap (x ^*) m

sigmaMean
  :: (RealFloat a, Traversable t, FStorable t, Space b) =>
       t (b a) -> b a -> b a
sigmaMean p p0= case sigmaMean' p p0 of
                  Just x -> x
                  Nothing -> error $ "sigmaMean nonConvergence"

sigmaMean1 (V1 p0 :|: p )= sigmaMean p p0

sigmaMean'
  :: (RealFloat a, Traversable t, FStorable t, Space b) =>
       t (b a) -> b a -> Maybe (b a)
sigmaMean' p p0 = go 1000 p0
    where
        go ! i ! mik
            | i == 0 = Nothing
            | Metric.quadrance m > 1.0e-14 = go (i-1) mik1
            | otherwise = Just mik1
            where
                mik1 = mik |+| m
                s = F.foldr1 (^+^) $  fmap (|-|mik) p
                m =  s ^/ n
                n = 2 * (fromIntegral $ fsize p)

sigmaCovXY xmean ymean samplex sampley  =  covxy
    where
        xstd = fmap (|-| xmean) samplex
        ystd = fmap (|-| ymean) sampley
        covxy =  fmap (^/ 2) $ F.foldr1 (!+!) (liftA2 mult xstd ystd)

sigmaCov mean sample = cov
    where
        std  = fmap (|-| mean) sample
        cov  =  fmap (^/2) $ F.foldr1 (!+!) (fmap square std)

sigmaPoints1 x p = V1 x :|: points
  where
    m = case potrf p of
            Right x -> x
            Left i -> error $ "cholesky factorization error: sigmapoints1 : "
    points = fmap (x|+|) $  m  :|: fmap ((-1) *^) m

sqrtPrediction
  :: (Applicative g1, Applicative g, Traversable g, Traversable g1,
      FStorable g, FStorable g1, FStorable (Local b1), Space b,
      Space b1) =>
     (b Double -> b1 Double)
     -> g1 (Local b1 Double)
     -> (b Double, g (Local b Double))
     -> (b1 Double, Local b1 (Local b1 Double))
sqrtPrediction f r (x,s)= (xmean,sk'')
    where
        (y,w0m,w0c,w1) = weigthsGamma (fromIntegral $ fsize s)
        sp = fmap (x|+|) $ fmap (fmap (*y)) s :|: fmap (fmap (*(-y))) s
        sp' = fmap f sp
        x'  = f x
        xmean = sigmaMean1 $ V1 x' :|: sp'
        std  =fmap (|-| xmean) sp'
        sk' = dgeqrf $ distribute $ (fmap (fmap (*sqrt w1)) std):|: r
        sk'' = distribute $  ch1upJ sk' $ fmap (*sqrt w0c) $ x' |-| xmean

weigthsGamma l =  weigths 1 2 0 l

weigths a b k l= (y,w0m,w0c,w1) where
    lam =  a^2 * (l + k ) - l
    y =  sqrt ( l + lam)
    w0m = lam /(l + lam)
    w0c = w0m  +  (1 - a^2 + b)
    w1 = 1 / (2*(l+ lam))

ch1dnR x = unRight . ch1dn  x
sqrtMeasure
  :: (Show (Local b Double),Show (b Double),Applicative g, Traversable g, FStorable g, FStorable (Local b),
      FStorable (Local b1), Space b1, Space b) =>
     (b1 Double -> b Double)
     -> g (Local b Double)
     -> b Double
     -> (b1 Double, Local b1 (Local b1 Double))
     -> (b1 Double, Local b1 (Local b1 Double))
sqrtMeasure h r m (x,s) =  (x'',s''')
    where
{-
        spx = V1 x :|: (fmap (x|+|) $ fmap (fmap (*y)) s :|: fmap (fmap (*(-y))) s )
        wc = V1 w0c :|: pure w1
        --wm = V1 w0m :|: pure w1
        V1 y' :|: _ =  spy
        xmean = sigmaMean1  spx
        stdx  =  liftA2 (*^) wc $ fmap (|-| xmean ) spx
        spy = fmap h spx
        ymean = sigmaMean1  spy
        stdy  =  fmap (|-| ymean) spy
        pxy =   distribute stdx !*!  stdy
        sk' =  dgeqrf $ distribute $ fmap (fmap (*sqrt w1))  stdy :|: r
        sk''= distribute .fromJust  $   ch1up sk' (fmap (*w0c) $ ymean |-| y')
-}
        (ymean,pxx,pxy) = sigmaSampling h r (x,s)
        k  =  pseudoInverse pxy pxx
        -- x' = xmean |+| (k !* (m |-| y')  )
        x' = x |+| (k !* (m |-| ymean)  )
        u  = pxx !*! distribute k
        s' = distribute $ F.foldl ch1dnR  (distribute s)  u
        (x'',s''',_) = sigmaSampling id V0 (x',s')
{-      spx' = fmap (x' |+| ) $ fmap (fmap (*y))  s' :|: fmap (fmap (*(-y))) s'
        x'' =  sigmaMean1 $ V1 x' :|: spx'
        stdx' =  fmap (|-| x'') spx'
        s'' = dgeqrf $ distribute $ fmap (fmap (*sqrt w1))  stdx'
        s''' = distribute $ fromJust $ ch1up s'' $ fmap (*sqrt w0c) $ x'' |-| x'
-}

pseudoInverse
  :: (Applicative nhrs, Applicative m, Applicative m1,
      Traversable nhrs, Traversable m, Traversable m1, Distributive m1,
      FStorable nhrs, FStorable m, FStorable m1) =>
     nhrs (m1 Double) -> m (m1 Double) -> nhrs (m1 Double)
pseudoInverse y x  = dgelssR (distribute x )  (dgelssR x  y )

sqrtBackward f p (x,s) (x1,s1) = (xb ,sb)
    where
        (xmean ,pxx,pxy) = sigmaSampling f p (x,s)
        k =  pseudoInverse pxy pxx
        xb = x |+| (k !* (x1 |-| xmean ))
        s1pxx = distribute $ F.foldl (\y-> ch1dnR y)  (distribute s1 )  (distribute pxx )
        u  = s1pxx !*! distribute k
        sb = distribute $ F.foldl  ch1upJ  (distribute s)  u

sigmaSampling h r (x,s) = (ymean,sxx,sxy)
    where
        (y,w0m,w0c,w1) = weigthsGamma (fromIntegral $ fsize s)
        spx = V1 x :|: (fmap (x|+|) $ fmap (fmap (*y)) s :|: fmap (fmap (*(-y))) s )
        wc = V1 w0c :|: pure w1
        --wm = V1 w0m :|: pure w1
        V1 y' :|: _ =  spy
        xmean = sigmaMean1  spx
        stdx  =  liftA2 (*^) wc $ fmap (|-| xmean ) spx
        spy = fmap h spx
        ymean = sigmaMean1  spy
        stdy  =  fmap (|-| ymean) spy
        sxy =   distribute stdx !*!  stdy
        sk' =  dgeqrf $ distribute $ fmap (fmap (*sqrt w1))  stdy :|: r
        sxx = distribute $ ch1upJ sk' (fmap (*w0c) $ ymean |-| y')

ch1upJ x = fromJust . ch1up x
dgelssR x = unRight . dgelss x

unRight (Right x) = x
unRight (Left x) = error $ "no right: "

max_stars_repo_name: massudaw/mtk

max_issues_repo_name: massudaw/mtk

max_forks_repo_name: massudaw/mtk


text:
{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE PatternGuards #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE TypeFamilies #-}
#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 702
{-# LANGUAGE Trustworthy #-}
{-# LANGUAGE DeriveGeneric #-}
#endif

#ifndef MIN_VERSION_vector
#define MIN_VERSION_vector(x,y,z) 1
#endif
-----------------------------------------------------------------------------
-- |
-- Copyright   :  (C) 2012-2015 Edward Kmett
-- License     :  BSD-style (see the file LICENSE)
--
-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
-- Stability   :  experimental
-- Portability :  non-portable
--
-- Quaternions
----------------------------------------------------------------------------
module Linear.Quaternion
  ( Quaternion(..)
  , Complicated(..)
  , Hamiltonian(..)
  , ee, ei, ej, ek
  , slerp
  , asinq
  , acosq
  , atanq
  , asinhq
  , acoshq
  , atanhq
  , absi
  , pow
  , rotate
  , axisAngle
  ) where

import Control.Applicative
import Control.DeepSeq (NFData(rnf))
import Control.Monad (liftM)
import Control.Monad.Fix
import Control.Monad.Zip
import Control.Lens hiding ((<.>))
import Data.Binary as Binary
import Data.Bytes.Serial
import Data.Complex (Complex((:+)))
import Data.Data
import Data.Distributive
import Data.Foldable
import Data.Functor.Bind
import Data.Functor.Classes
import Data.Functor.Rep
import Data.Hashable
import Data.Serialize as Cereal
import GHC.Arr (Ix(..))
import qualified Data.Foldable as F
import Data.Monoid
import Foreign.Ptr (castPtr, plusPtr)
import Foreign.Storable (Storable(..))
#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 702
import GHC.Generics (Generic)
#endif
#if defined(__GLASGOW_HASKELL__) && __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.Conjugate
import Linear.Metric
import Linear.V3
import Linear.Vector
import Prelude hiding (any)

{-# ANN module "HLint: ignore Reduce duplication" #-}

-- | Quaternions
data Quaternion a = Quaternion !a {-# UNPACK #-}!(V3 a)
                    deriving (Eq,Ord,Read,Show,Data,Typeable
#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 702
                             ,Generic
#endif
#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 706
                             ,Generic1
#endif
                             )

instance Functor Quaternion where
  fmap f (Quaternion e v) = Quaternion (f e) (fmap f v)
  {-# INLINE fmap #-}
  a <$ _ = Quaternion a (V3 a a a)
  {-# INLINE (<$) #-}

instance Apply Quaternion where
  Quaternion f fv <.> Quaternion a v = Quaternion (f a) (fv <.> v)
  {-# INLINE (<.>) #-}

instance Applicative Quaternion where
  pure a = Quaternion a (pure a)
  {-# INLINE pure #-}
  Quaternion f fv <*> Quaternion a v = Quaternion (f a) (fv <*> v)
  {-# INLINE (<*>) #-}

instance Additive Quaternion where
  zero = pure 0
  {-# INLINE zero #-}
  liftU2 = liftA2
  {-# INLINE liftU2 #-}
  liftI2 = liftA2
  {-# INLINE liftI2 #-}

instance Bind Quaternion where
  Quaternion a (V3 b c d) >>- f = Quaternion a' (V3 b' c' d') where
    Quaternion a' _          = f a
    Quaternion _ (V3 b' _ _) = f b
    Quaternion _ (V3 _ c' _) = f c
    Quaternion _ (V3 _ _ d') = f d
  {-# INLINE (>>-) #-}

instance Monad Quaternion where
  return = pure
  {-# INLINE return #-}
  -- the diagonal of a sedenion is super useful!
  Quaternion a (V3 b c d) >>= f = Quaternion a' (V3 b' c' d') where
    Quaternion a' _          = f a
    Quaternion _ (V3 b' _ _) = f b
    Quaternion _ (V3 _ c' _) = f c
    Quaternion _ (V3 _ _ d') = f d
  {-# INLINE (>>=) #-}

instance Ix a => Ix (Quaternion a) where
    {-# SPECIALISE instance Ix (Quaternion Int) #-}

    range (Quaternion l1 l2, Quaternion u1 u2) =
      [ Quaternion i1 i2 | i1 <- range (l1,u1), i2 <- range (l2,u2) ]
    {-# INLINE range #-}

    unsafeIndex (Quaternion l1 l2, Quaternion u1 u2) (Quaternion i1 i2) =
      unsafeIndex (l1,u1) i1 * unsafeRangeSize (l2,u2) + unsafeIndex (l2,u2) i2
    {-# INLINE unsafeIndex #-}

    inRange (Quaternion l1 l2, Quaternion u1 u2) (Quaternion i1 i2) =
      inRange (l1,u1) i1 && inRange (l2,u2) i2
    {-# INLINE inRange #-}

instance Representable Quaternion where
  type Rep Quaternion = E Quaternion
  tabulate f = Quaternion (f ee) (V3 (f ei) (f ej) (f ek))
  {-# INLINE tabulate #-}
  index xs (E l) = view l xs
  {-# INLINE index #-}

instance FunctorWithIndex (E Quaternion) Quaternion where
  imap f (Quaternion a (V3 b c d)) = Quaternion (f ee a) $ V3 (f ei b) (f ej c) (f ek d)
  {-# INLINE imap #-}

instance FoldableWithIndex (E Quaternion) Quaternion where
  ifoldMap f (Quaternion a (V3 b c d)) = f ee a `mappend` f ei b `mappend` f ej c `mappend` f ek d
  {-# INLINE ifoldMap #-}

instance TraversableWithIndex (E Quaternion) Quaternion where
  itraverse f (Quaternion a (V3 b c d)) = Quaternion <$> f ee a <*> (V3 <$> f ei b <*> f ej c <*> f ek d)
  {-# INLINE itraverse #-}

type instance Index (Quaternion a) = E Quaternion
type instance IxValue (Quaternion a) = a

instance Ixed (Quaternion a) where
  ix = el
  {-# INLINE ix #-}

instance Each (Quaternion a) (Quaternion b) a b where
  each = traverse
  {-# INLINE each #-}

instance Foldable Quaternion where
  foldMap f (Quaternion e v) = f e `mappend` foldMap f v
  {-# INLINE foldMap #-}
  foldr f z (Quaternion e v) = f e (F.foldr f z v)
  {-# INLINE foldr #-}

instance Traversable Quaternion where
  traverse f (Quaternion e v) = Quaternion <$> f e <*> traverse f v
  {-# INLINE traverse #-}

instance Storable a => Storable (Quaternion a) where
  sizeOf _ = 4 * sizeOf (undefined::a)
  {-# INLINE sizeOf #-}
  alignment _ = alignment (undefined::a)
  {-# INLINE alignment #-}
  poke ptr (Quaternion e v) = poke (castPtr ptr) e >>
                              poke (castPtr (ptr `plusPtr` sz)) v
    where sz = sizeOf (undefined::a)
  {-# INLINE poke #-}
  peek ptr = Quaternion <$> peek (castPtr ptr)
                        <*> peek (castPtr (ptr `plusPtr` sz))
    where sz = sizeOf (undefined::a)
  {-# INLINE peek #-}

instance RealFloat a => Num (Quaternion a) where
  {-# SPECIALIZE instance Num (Quaternion Float) #-}
  {-# SPECIALIZE instance Num (Quaternion Double) #-}
  (+) = liftA2 (+)
  {-# INLINE (+) #-}
  (-) = liftA2 (-)
  {-# INLINE (-) #-}
  negate = fmap negate
  {-# INLINE negate #-}
  Quaternion s1 v1 * Quaternion s2 v2 = Quaternion (s1*s2 - (v1 `dot` v2)) $
                                        (v1 `cross` v2) + s1*^v2 + s2*^v1
  {-# INLINE (*) #-}
  fromInteger x = Quaternion (fromInteger x) 0
  {-# INLINE fromInteger #-}
  abs z = Quaternion (norm z) 0
  {-# INLINE abs #-}
  signum q@(Quaternion e (V3 i j k))
    | m == 0.0 = q
    | not (isInfinite m || isNaN m) = q ^/ sqrt m
    | any isNaN q = qNaN
    | not (ii || ij || ik) = Quaternion 1 (V3 0 0 0)
    | not (ie || ij || ik) = Quaternion 0 (V3 1 0 0)
    | not (ie || ii || ik) = Quaternion 0 (V3 0 1 0)
    | not (ie || ii || ij) = Quaternion 0 (V3 0 0 1)
    | otherwise = qNaN
    where
      m = quadrance q
      ie = isInfinite e
      ii = isInfinite i
      ij = isInfinite j
      ik = isInfinite k
  {-# INLINE signum #-}

instance Hashable a => Hashable (Quaternion a) where
  hashWithSalt s (Quaternion a b) = s `hashWithSalt` a `hashWithSalt` b
  {-# INLINE hashWithSalt #-}

qNaN :: RealFloat a => Quaternion a
qNaN = Quaternion fNaN (V3 fNaN fNaN fNaN) where fNaN = 0/0
{-# INLINE qNaN #-}

-- {-# RULES "abs/norm" abs x = Quaternion (norm x) 0 #-}
-- {-# RULES "signum/signorm" signum = signorm #-}

-- this will attempt to rewrite calls to abs to use norm intead when it is available.

instance RealFloat a => Fractional (Quaternion a) where
  {-# SPECIALIZE instance Fractional (Quaternion Float) #-}
  {-# SPECIALIZE instance Fractional (Quaternion Double) #-}
  Quaternion q0 (V3 q1 q2 q3) / Quaternion r0 (V3 r1 r2 r3) =
    Quaternion (r0*q0+r1*q1+r2*q2+r3*q3)
               (V3 (r0*q1-r1*q0-r2*q3+r3*q2)
                   (r0*q2+r1*q3-r2*q0-r3*q1)
                   (r0*q3-r1*q2+r2*q1-r3*q0))
               ^/ (r0*r0 + r1*r1 + r2*r2 + r3*r3)
  {-# INLINE (/) #-}
  recip q = q ^/ quadrance q
  {-# INLINE recip #-}
  fromRational x = Quaternion (fromRational x) 0
  {-# INLINE fromRational #-}

instance Metric Quaternion where
  Quaternion e v `dot` Quaternion e' v' = e*e' + (v `dot` v')
  {-# INLINE dot #-}

-- | A vector space that includes the basis elements '_e' and '_i'
class Complicated t where
  _e, _i :: Lens' (t a) a

ee, ei :: Complicated t => E t
ee = E _e
ei = E _i

instance Complicated Complex where
  _e f (a :+ b) = (:+ b) <$> f a
  {-# INLINE _e #-}
  _i f (a :+ b) = (a :+) <$> f b
  {-# INLINE _i #-}

instance Complicated Quaternion where
  _e f (Quaternion a v) = (`Quaternion` v) <$> f a
  {-# INLINE _e #-}
  _i f (Quaternion a v) = Quaternion a <$> _x f v
  {-# INLINE _i #-}

-- | A vector space that includes the basis elements '_e', '_i', '_j' and '_k'
class Complicated t => Hamiltonian t where
  _j, _k :: Lens' (t a) a
  _ijk :: Lens' (t a) (V3 a)

ej, ek :: Hamiltonian t => E t
ej = E _j
ek = E _k

instance Hamiltonian Quaternion where
  _j f (Quaternion a v) = Quaternion a <$> _y f v
  {-# INLINE _j #-}
  _k f (Quaternion a v) = Quaternion a <$> _z f v
  {-# INLINE _k #-}
  _ijk f (Quaternion a v) = Quaternion a <$> f v
  {-# INLINE _ijk #-}

instance Distributive Quaternion where
  distribute f = Quaternion (fmap (\(Quaternion x _) -> x) f) $ V3
    (fmap (\(Quaternion _ (V3 y _ _)) -> y) f)
    (fmap (\(Quaternion _ (V3 _ z _)) -> z) f)
    (fmap (\(Quaternion _ (V3 _ _ w)) -> w) f)
  {-# INLINE distribute #-}

instance (Conjugate a, RealFloat a) => Conjugate (Quaternion a) where
  conjugate (Quaternion e v) = Quaternion (conjugate e) (negate v)
  {-# INLINE conjugate #-}

reimagine :: RealFloat a => a -> a -> Quaternion a -> Quaternion a
reimagine r s (Quaternion _ v)
  | isNaN s || isInfinite s = let aux 0 = 0
                                  aux x = s * x
                              in Quaternion r (aux <$> v)
  | otherwise = Quaternion r (v^*s)
{-# INLINE reimagine #-}

-- | quadrance of the imaginary component
qi :: Num a => Quaternion a -> a
qi (Quaternion _ v) = quadrance v
{-# INLINE qi #-}

-- | norm of the imaginary component
absi :: Floating a => Quaternion a -> a
absi = sqrt . qi
{-# INLINE absi #-}

-- | raise a 'Quaternion' to a scalar power
pow :: RealFloat a => Quaternion a -> a -> Quaternion a
pow q t = exp (t *^ log q)
{-# INLINE pow #-}

sqrte2pqiq :: (Floating a, Ord a) => a -> a -> a
sqrte2pqiq e qiq -- = sqrt (e*e) + qiq
  | e < - 1.5097698010472593e153 = -(qiq/e) - e
  | e < 5.582399551122541e57      = sqrt ((e*e) + qiq) -- direct definition
  | otherwise                     = (qiq/e) + e
-- {-# SPECIALIZE sqrte2pqiq :: Double -> Double -> Double #-}
-- {-# SPECIALIZE sqrte2pqiq :: Float -> Float -> Float #-}
#ifdef HERBIE
{-# ANN sqrte2pqiq "NoHerbie" #-}
#endif

tanrhs :: (Floating a, Ord a) => a -> a -> a -> a
tanrhs sai ai d -- = cosh ai * (sai / ai) / d -- improved from 6.04 bits of error to 0.19 bits
  | sai < -4.618902267687042e-52 = (sai / d / ai) * cosh ai
  | sai < 1.038530535935153e-39 = (cosh ai * sai) / ai / d
  | otherwise = (sai / d / ai) * cosh ai
-- {-# SPECIALIZE tanrhs :: Double -> Double -> Double -> Double #-}
-- {-# SPECIALIZE tanrhs :: Float -> Float -> Float -> Float #-}
#ifdef HERBIE
{-# ANN tanrhs "NoHerbie" #-}
#endif

-- ehh..
instance RealFloat a => Floating (Quaternion a) where
  {-# SPECIALIZE instance Floating (Quaternion Float) #-}
  {-# SPECIALIZE instance Floating (Quaternion Double) #-}
  pi = Quaternion pi 0
  {-# INLINE pi #-}
  exp q@(Quaternion e v)
    | qiq == 0 = Quaternion (exp e) v
    | ai <- sqrt qiq, exe <- exp e = reimagine (exe * cos ai) (exe * (sin ai / ai)) q
    where qiq = qi q
  {-# INLINE exp #-}
  log q@(Quaternion e v@(V3 _i j k))
    | qiq == 0 = if e >= 0
                 then Quaternion (log e) v
                 else Quaternion (log (negate e)) (V3 pi j k) -- mmm, pi
    | ai <- sqrt qiq = reimagine (log m) (atan2 m e / ai) q
    where qiq = qi q
          m = sqrte2pqiq e qiq
  {-# INLINE log #-}

  x ** y = exp (y * log x)
  {-# INLINE (**) #-}

  sqrt q@(Quaternion e v)
    | m   == 0 = q
    | qiq == 0 = if e > 0
                 then Quaternion (sqrt e) 0
                 else Quaternion 0 (V3 (sqrt (negate e)) 0 0)
    | im <- sqrt (0.5*(m-e)) / sqrt qiq = Quaternion (0.5*(m+e)) (v^*im)
    where qiq = qi q
          m = sqrte2pqiq e qiq
  {-# INLINE sqrt #-}

  cos q@(Quaternion e v)
    | qiq == 0 = Quaternion (cos e) v
    | ai <- sqrt qiq = reimagine (cos e * cosh ai) (- sin e / ai / sinh ai) q -- 0.15 bits error
                    -- reimagine (cos e * cosh ai) (- sin e * sinh ai / ai) q -- 13.5 bits worse
    where qiq = qi q
  {-# INLINE cos #-}

  sin q@(Quaternion e v)
    | qiq == 0 = Quaternion (sin e) v
    | ai <- sqrt qiq = reimagine (sin e * cosh ai) (cos e * sinh ai / ai) q
    where qiq = qi q
  {-# INLINE sin #-}

  tan q@(Quaternion e v)
    | qiq == 0 = Quaternion (tan e) v
    | ai <- sqrt qiq, ce <- cos e, sai <- sinh ai, d <- ce*ce + sai*sai =
      reimagine (ce * sin e / d) (tanrhs sai ai d) q
    where qiq = qi q
  {-# INLINE tan #-}

  sinh q@(Quaternion e v)
    | qiq == 0 = Quaternion (sinh e) v
    | ai <- sqrt qiq = reimagine (sinh e * cos ai) (cosh e * sin ai / ai) q
    where qiq = qi q
  {-# INLINE sinh #-}

  cosh q@(Quaternion e v)
    | qiq == 0 = Quaternion (cosh e) v
    | ai <- sqrt qiq = reimagine (cosh e * cos ai) (sin ai * (sinh e / ai)) q
    where qiq = qi q
  {-# INLINE cosh #-}

  tanh q@(Quaternion e v)
    | qiq == 0 = Quaternion (tanh e) v
    | ai <- sqrt qiq, se <- sinh e, cai <- cos ai, d <- se*se + cai*cai =
      reimagine (cosh e * se / d) (tanhrhs cai ai d) q
    where qiq = qi q
  {-# INLINE tanh #-}

  asin = cut asin
  {-# INLINE asin #-}
  acos = cut acos
  {-# INLINE acos #-}
  atan = cut atan
  {-# INLINE atan #-}

  asinh = cut asinh
  {-# INLINE asinh #-}
  acosh = cut acosh
  {-# INLINE acosh #-}
  atanh = cut atanh
  {-# INLINE atanh #-}

tanhrhs :: (Floating a, Ord a) => a -> a -> a -> a
tanhrhs cai ai d -- = cai * (sin ai / ai) / d
  | d >= -4.2173720203427147e-29 && d < 4.446702369113811e64 = cai / (d * (ai / sin ai))
  | otherwise                                                = cai * (1 / ai / sin ai) / d
-- {-# SPECIALIZE tanhrhs :: Double -> Double -> Double -> Double #-}
-- {-# SPECIALIZE tanhrhs :: Float -> Float -> Float -> Float #-}
#ifdef HERBIE
{-# ANN tanhrhs "NoHerbie" #-}
#endif

-- | Helper for calculating with specific branch cuts
cut :: RealFloat a => (Complex a -> Complex a) -> Quaternion a -> Quaternion a
cut f q@(Quaternion e (V3 _ y z))
  | qiq == 0 = Quaternion a (V3 b y z)
  | otherwise = reimagine a (b / ai) q
  where qiq = qi q
        ai = sqrt qiq
        a :+ b = f (e :+ ai)
{-# INLINE cut #-}

-- | Helper for calculating with specific branch cuts
cutWith :: RealFloat a => Complex a -> Quaternion a -> Quaternion a
cutWith (r :+ im) q@(Quaternion e v)
  | e /= 0 || qiq == 0 || isNaN qiq || isInfinite qiq = error "bad cut"
  | s <- im / sqrt qiq = Quaternion r (v^*s)
  where qiq = qi q
{-# INLINE cutWith #-}

-- | 'asin' with a specified branch cut.
asinq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a
asinq q@(Quaternion e _) u
  | qiq /= 0.0 || e >= -1 && e <= 1 = asin q
  | otherwise = cutWith (asin (e :+ sqrt qiq)) u
  where qiq = qi q
{-# INLINE asinq #-}

-- | 'acos' with a specified branch cut.
acosq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a
acosq q@(Quaternion e _) u
  | qiq /= 0.0 || e >= -1 && e <= 1 = acos q
  | otherwise = cutWith (acos (e :+ sqrt qiq)) u
  where qiq = qi q
{-# INLINE acosq #-}

-- | 'atan' with a specified branch cut.
atanq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a
atanq q@(Quaternion e _) u
  | e /= 0.0 || qiq >= -1 && qiq <= 1 = atan q
  | otherwise = cutWith (atan (e :+ sqrt qiq)) u
  where qiq = qi q
{-# INLINE atanq #-}

-- | 'asinh' with a specified branch cut.
asinhq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a
asinhq q@(Quaternion e _) u
  | e /= 0.0 || qiq >= -1 && qiq <= 1 = asinh q
  | otherwise = cutWith (asinh (e :+ sqrt qiq)) u
  where qiq = qi q
{-# INLINE asinhq #-}

-- | 'acosh' with a specified branch cut.
acoshq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a
acoshq q@(Quaternion e _) u
  | qiq /= 0.0 || e >= 1 = asinh q
  | otherwise = cutWith (acosh (e :+ sqrt qiq)) u
  where qiq = qi q
{-# INLINE acoshq #-}

-- | 'atanh' with a specified branch cut.
atanhq :: RealFloat a => Quaternion a -> Quaternion a -> Quaternion a
atanhq q@(Quaternion e _) u
  | qiq /= 0.0 || e > -1 && e < 1 = atanh q
  | otherwise = cutWith (atanh (e :+ sqrt qiq)) u
  where qiq = qi q
{-# INLINE atanhq #-}

-- | Spherical linear interpolation between two quaternions.

slerp :: RealFloat a => Quaternion a -> Quaternion a -> a -> Quaternion a
slerp q p t
  | 1.0 - cosphi < 1e-8 = q
  | otherwise           = ((sin ((1-t)*phi) *^ q) + sin (t*phi) *^ f p) ^/ sin phi
  where
    dqp = dot q p
    (cosphi, f) = if dqp < 0 then (-dqp, negate) else (dqp, id)
    phi = acos cosphi
{-# SPECIALIZE slerp :: Quaternion Float -> Quaternion Float -> Float -> Quaternion Float #-}
{-# SPECIALIZE slerp :: Quaternion Double -> Quaternion Double -> Double -> Quaternion Double #-}

-- | Apply a rotation to a vector.
rotate :: (Conjugate a, RealFloat a) => Quaternion a -> V3 a -> V3 a
rotate q v = ijk where
  Quaternion _ ijk = q * Quaternion 0 v * conjugate q
{-# SPECIALIZE rotate :: Quaternion Float -> V3 Float -> V3 Float #-}
{-# SPECIALIZE rotate :: Quaternion Double -> V3 Double -> V3 Double #-}

instance (RealFloat a, Epsilon a) => Epsilon (Quaternion a) where
  nearZero = nearZero . quadrance
  {-# INLINE nearZero #-}

-- | @'axisAngle' axis theta@ builds a 'Quaternion' representing a
-- rotation of @theta@ radians about @axis@.
axisAngle :: (Epsilon a, Floating a) => V3 a -> a -> Quaternion a
axisAngle axis theta = Quaternion (cos half) (sin half *^ normalize axis)
  where half = theta / 2
{-# INLINE axisAngle #-}

data instance U.Vector    (Quaternion a) =  V_Quaternion !Int (U.Vector    a)
data instance U.MVector s (Quaternion a) = MV_Quaternion !Int (U.MVector s a)
instance U.Unbox a => U.Unbox (Quaternion a)

instance U.Unbox a => M.MVector U.MVector (Quaternion a) where
  basicLength (MV_Quaternion n _) = n
  basicUnsafeSlice m n (MV_Quaternion _ v) = MV_Quaternion n (M.basicUnsafeSlice (4*m) (4*n) v)
  basicOverlaps (MV_Quaternion _ v) (MV_Quaternion _ u) = M.basicOverlaps v u
  basicUnsafeNew n = liftM (MV_Quaternion n) (M.basicUnsafeNew (4*n))
  basicUnsafeRead (MV_Quaternion _ v) i =
    do let o = 4*i
       x <- M.basicUnsafeRead v o
       y <- M.basicUnsafeRead v (o+1)
       z <- M.basicUnsafeRead v (o+2)
       w <- M.basicUnsafeRead v (o+3)
       return (Quaternion x (V3 y z w))
  basicUnsafeWrite (MV_Quaternion _ v) i (Quaternion x (V3 y z w)) =
    do let o = 4*i
       M.basicUnsafeWrite v o     x
       M.basicUnsafeWrite v (o+1) y
       M.basicUnsafeWrite v (o+2) z
       M.basicUnsafeWrite v (o+3) w
#if MIN_VERSION_vector(0,11,0)
  basicInitialize (MV_Quaternion _ v) = M.basicInitialize v
#endif

instance U.Unbox a => G.Vector U.Vector (Quaternion a) where
  basicUnsafeFreeze (MV_Quaternion n v) = liftM ( V_Quaternion n) (G.basicUnsafeFreeze v)
  basicUnsafeThaw   ( V_Quaternion n v) = liftM (MV_Quaternion n) (G.basicUnsafeThaw   v)
  basicLength       ( V_Quaternion n _) = n
  basicUnsafeSlice m n (V_Quaternion _ v) = V_Quaternion n (G.basicUnsafeSlice (4*m) (4*n) v)
  basicUnsafeIndexM (V_Quaternion _ v) i =
    do let o = 4*i
       x <- G.basicUnsafeIndexM v o
       y <- G.basicUnsafeIndexM v (o+1)
       z <- G.basicUnsafeIndexM v (o+2)
       w <- G.basicUnsafeIndexM v (o+3)
       return (Quaternion x (V3 y z w))

instance MonadZip Quaternion where
  mzipWith = liftA2

instance MonadFix Quaternion where
  mfix f = Quaternion (let Quaternion a _ = f a in a)
                      (V3 (let Quaternion _ (V3 a _ _) = f a in a)
                          (let Quaternion _ (V3 _ a _) = f a in a)
                          (let Quaternion _ (V3 _ _ a) = f a in a))

instance NFData a => NFData (Quaternion a) where
  rnf (Quaternion a b) = rnf a `seq` rnf b

instance Serial1 Quaternion where
  serializeWith f (Quaternion a b) = f a >> serializeWith f b
  deserializeWith f = Quaternion <$> f <*> deserializeWith f

instance Serial a => Serial (Quaternion a) where
  serialize = serializeWith serialize
  deserialize = deserializeWith deserialize

instance Binary a => Binary (Quaternion a) where
  put = serializeWith Binary.put
  get = deserializeWith Binary.get

instance Serialize a => Serialize (Quaternion a) where
  put = serializeWith Cereal.put
  get = deserializeWith Cereal.get

instance Eq1 Quaternion where eq1 = (==)
instance Ord1 Quaternion where compare1 = compare
instance Show1 Quaternion where showsPrec1 = showsPrec
instance Read1 Quaternion where readsPrec1 = readsPrec

max_stars_repo_name: bgamari/linear

max_issues_repo_name: bgamari/linear

max_forks_repo_name: bgamari/linear


text:
{-# LANGUAGE RankNTypes #-}

-- |
--  Author: Benedikt Spies
--
--  Utility functions to generate simple Images
--  Image rows are generated in parallel
module Utils.ImageGenerator (generateRainbowImage, generateMandelbrotImage) where

import Codec.Picture
  ( Image,
    Pixel8,
    PixelRGB8 (PixelRGB8),
  )
import Data.Complex (Complex ((:+)), magnitude)
import Utils.ParallelImageUtils (parGenerateImage)

generateRainbowImage :: Image PixelRGB8
generateRainbowImage = parGenerateImage f 256 256
  where
    f :: Int -> Int -> PixelRGB8
    f x y = PixelRGB8 r g b
      where
        r = fromIntegral $ 255 - x
        g = fromIntegral $ 255 - y
        b = fromIntegral x

generateMandelbrotImage :: Int -> Int -> Image PixelRGB8
generateMandelbrotImage width height = parGenerateImage (generateMandelbrotPixel width height) width height

generateMandelbrotPixel :: Int -> Int -> Int -> Int -> PixelRGB8
generateMandelbrotPixel imageWidth imageHeight x y = color $ mandelbrot (x' / w' * 3 - 2) (y' / h' * 3 - 1.5)
  where
    x' = fromIntegral x
    y' = fromIntegral y
    w' = fromIntegral imageWidth
    h' = fromIntegral imageHeight
    maxIter = 300
    color :: Int -> PixelRGB8
    color iter = PixelRGB8 colorVal colorVal colorVal
      where
        colorVal = floor $ 255 - amt * 255 :: Pixel8
        amt = fromIntegral iter / fromIntegral maxIter :: Float
    mandelbrot :: Float -> Float -> Int
    mandelbrot r i =
      length . takeWhile (\z -> magnitude z <= 2)
        . take maxIter
        $ iterate (\z -> z ^ (2 :: Int) + (r :+ i)) 0

max_stars_repo_name: birneee/haskell-spotify-tui

max_issues_repo_name: birneee/haskell-spotify-tui

max_forks_repo_name: birneee/haskell-spotify-tui


text:
{-# LANGUAGE RecordWildCards, BangPatterns #-}

module NeuralNetwork where

import qualified  Statistics.Matrix as M
import            Control.Monad.Random

-- hyperparameters

data NeuralNetwork =
     NeuralNetwork { inputLayerSize  :: Int
                   , outputLayerSize :: Int
                   , hiddenLayerSize :: Int
                   , weight1Matrix   :: M.Matrix
                   , weight2Matrix   :: M.Matrix
                   }

initNetwork :: MonadRandom r => Int -> Int -> Int -> r NeuralNetwork
initNetwork inputSize outputSize hiddenSize = do
    rs  <- getRandomRs (0.0, 1.0)
    rs' <- getRandomRs (0.0, 1.0)
    return $ NeuralNetwork { inputLayerSize  = inputSize
                           , outputLayerSize = outputSize
                           , hiddenLayerSize = hiddenSize
                           , weight1Matrix   = M.fromList inputSize  hiddenSize $ take (inputSize  * hiddenSize) rs
                           , weight2Matrix   = M.fromList hiddenSize outputSize $ take (outputSize * hiddenSize) rs'
                           }

forwardRs :: NeuralNetwork -> M.Matrix -> (M.Matrix, M.Matrix, M.Matrix, M.Matrix)
forwardRs NeuralNetwork{..} matrix =
    let !z2 = matrix * weight1Matrix
        !a2 = M.map sigmoid z2
        !z3 = a2 * weight2Matrix
    in  (M.map sigmoid z3, z2, a2, z3)

forward :: NeuralNetwork -> M.Matrix -> M.Matrix
forward nn matrix =
    let (r, _, _, _) = forwardRs nn matrix
    in r

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

-- | Derivative of sigmoid function
sigmoidPrime :: Double -> Double
sigmoidPrime z = ez / ((1 + ez) ** 2)
    where ez = exp (-z)

costFunction :: NeuralNetwork -> M.Matrix -> M.Matrix -> Double
costFunction nn@NeuralNetwork{..} x y =
    let yHat = nn `forward` x
        err  = sum . map (**2) . M.toList $ y - yHat
    in  0.5 * err

costFunctionPrime :: NeuralNetwork -> M.Matrix -> M.Matrix -> (M.Matrix, M.Matrix)
costFunctionPrime nn@NeuralNetwork{..} x y =
    let (yHat, z2, a2, z3) = nn `forwardRs` x
        !delta3            =  (- (y - yHat)) `pointMultiply` (M.map sigmoidPrime z3)
        !dJdW2             = (M.transpose a2) * delta3

        !delta2            = delta3 * (M.transpose weight2Matrix) * (M.map sigmoidPrime z2)
        !dJdW1             = M.transpose x * delta2

    in (dJdW1, dJdW2)

pointMultiply :: M.Matrix -> M.Matrix -> M.Matrix
pointMultiply m1 m2 =
    let (r, c) = (M.rows m1, M.cols m1)
        (l1, l2) = (M.toList m1, M.toList m2)
    in M.fromList r c (foldr (\(a,b) xs -> (a * b) : xs) [] (zip l1 l2))

instance Num M.Matrix where
    m1 + m2  = let (r, c) = (M.rows m1, M.cols m1)
                   (l1, l2) = (M.toList m1, M.toList m2)
               in M.fromList r c (foldr (\(a,b) xs -> (a + b) : xs) [] (zip l1 l2))

    negate m      = let (r, c) = (M.rows m, M.cols m)
                        l      = M.toList m
                    in M.fromList r c (map negate l)

    m1 * m2  = m1 `M.multiply` m2

    abs m         = undefined
    signum m      = undefined
    fromInteger m = undefined

max_stars_repo_name: cirquit/nn-playground

max_issues_repo_name: cirquit/nn-playground

max_forks_repo_name: cirquit/nn-playground


text:
import Graphics.Gloss
import Graphics.Gloss.Juicy
import Codec.Picture
import Codec.Picture.RGBA8
import Data.Complex

width = 512
height = width `div` 2

mandelbrotness :: Complex Float -> Complex Float -> Float
mandelbrotness c s = (fromIntegral $ length $ takeWhile (\c -> magnitude c <= 2) $ take maxIter $ z c) / (fromIntegral maxIter) where
  z c = iterate (\x -> x^2 + c) s
  maxIter = 32

image :: Float -> Image PixelRGB8
image time = generateImage (\x y -> let [x',y',w,h] = map fromIntegral [x,y,width,height] in mandelbrotColor ((4*x'/w - 2.5) :+ (2*y'/h - 1)) (0.1*time :+ 0) ) width height where
  mandelbrotColor c s = PixelRGB8 (floor $ 255 * mandelbrotness c s) 120 120

main :: IO ()
main = animate (InWindow "Mandelbrot" (width,height) (0,0)) black (fromImageRGB8 . image)

max_stars_repo_name: cirquit/random-code

max_issues_repo_name: cirquit/Personal-Repository

max_forks_repo_name: cirquit/Personal-Repository


text:
{-# LANGUAGE BangPatterns, GeneralizedNewtypeDeriving, StandaloneDeriving #-}

import           Prelude hiding (mapM)

import           Options.Applicative
import           Data.Monoid ((<>))
import           Control.Monad.Trans.Class

import           Data.Vector (Vector)
import qualified Data.Vector.Generic as V
import           Statistics.Sample (mean)

import           Data.Traversable (mapM)
import qualified Data.Set as S
import           Data.Set (Set)
import qualified Data.Map.Strict as M

import           ReadData
import           SerializeText
import qualified RunSampler as Sampler
import           BayesStack.DirMulti
import           BayesStack.Models.Topic.LDA
import           BayesStack.UniqueKey

import qualified Data.Text as T
import qualified Data.Text.IO as TIO

import           System.FilePath.Posix ((</>))
import           Data.Binary
import qualified Data.ByteString as BS
import           Text.Printf

import           Data.Random
import           System.Random.MWC

data RunOpts = RunOpts { nodesFile       :: FilePath
                       , stopwords       :: Maybe FilePath
                       , nTopics         :: Int
                       , samplerOpts     :: Sampler.SamplerOpts
                       , hyperParams     :: HyperParams
                       }

runOpts :: Parser RunOpts
runOpts = RunOpts
    <$> strOption  ( long "nodes"
                  <> short 'n'
                  <> metavar "FILE"
                  <> help "File containing nodes and their associated items"
                   )
    <*> nullOption ( long "stopwords"
                  <> short 's'
                  <> metavar "FILE"
                  <> reader (pure . Just)
                  <> value Nothing
                  <> help "Stop word list"
                   )
    <*> option     ( long "topics"
                  <> short 't'
                  <> metavar "N"
                  <> value 20
                  <> help "Number of topics"
                   )
    <*> Sampler.samplerOpts
    <*> hyperOpts

hyperOpts = HyperParams
    <$> option     ( long "prior-theta"
                  <> value 1
                  <> help "Dirichlet parameter for prior on theta"
                   )
    <*> option     ( long "prior-phi"
                  <> value 0.1
                  <> help "Dirichlet parameter for prior on phi"
                   )

mapMKeys :: (Ord k, Ord k', Monad m, Applicative m)
         => (a -> m a') -> (k -> m k') -> M.Map k a -> m (M.Map k' a')
mapMKeys f g x = M.fromList <$> (mapM (\(k,v)->(,) <$> g k <*> f v) $ M.assocs x)

termsToItems :: M.Map NodeName [Term]
             -> (M.Map Node [Item], (M.Map Item Term, M.Map Node NodeName))
termsToItems nodes =
    let ((d', nodeMap), itemMap) =
            runUniqueKey' [Item i | i <- [0..]] $
            runUniqueKeyT' [Node i | i <- [0..]] $ do
                mapMKeys (mapM (lift . getUniqueKey)) getUniqueKey nodes
    in (d', (itemMap, nodeMap))

netData :: HyperParams -> M.Map Node [Item] -> Int -> NetData
netData hp nodeItems nTopics =
    NetData { dHypers           = hp
            , dItems            = S.unions $ map S.fromList $ M.elems nodeItems
            , dTopics           = S.fromList [Topic i | i <- [1..nTopics]]
            , dNodeItems        = M.fromList
                                  $ zip [NodeItem i | i <- [0..]]
                                  $ do (n,items) <- M.assocs nodeItems
                                       item <- items
                                       return (n, item)
            , dNodes            = M.keysSet nodeItems
            }

opts :: ParserInfo RunOpts
opts = info runOpts (  fullDesc
                    <> progDesc "Learn LDA model"
                    <> header "run-lda - learn LDA model"
                    )

instance Sampler.SamplerModel MState where
    estimateHypers = reestimate
    modelLikelihood = modelLikelihood
    summarizeHypers ms =
        "  phi  : "++show (dmAlpha $ snd $ M.findMin $ stPhis ms)++"\n"++
        "  theta: "++show (dmAlpha $ snd $ M.findMin $ stThetas ms)++"\n"

main :: IO ()
main = do
    args <- execParser opts
    stopWords <- case stopwords args of
                     Just f  -> S.fromList . T.words <$> TIO.readFile f
                     Nothing -> return S.empty
    printf "Read %d stopwords\n" (S.size stopWords)

    (nodeItems, (itemMap, nodeMap)) <- termsToItems
                            <$> readNodeItems stopWords (nodesFile args)

    Sampler.createSweeps $ samplerOpts args
    let sweepsDir = Sampler.sweepsDir $ samplerOpts args
    encodeFile (sweepsDir </> "item-map") itemMap
    encodeFile (sweepsDir </> "node-map") nodeMap

    let termCounts = V.fromListN (M.size nodeItems)
                     $ map length $ M.elems nodeItems :: Vector Int
    printf "Read %d nodes\n" (M.size nodeItems)
    printf "Mean items per node:  %1.2f\n" (mean $ V.map realToFrac termCounts)

    withSystemRandom $ \mwc->do
    let nd = netData (hyperParams args) nodeItems (nTopics args)
    encodeFile (sweepsDir </> "data") nd
    mInit <- runRVar (randomInitialize nd) mwc
    let m = model nd mInit
    Sampler.runSampler (samplerOpts args) m (updateUnits nd)
    return ()

max_stars_repo_name: bgamari/bayes-stack

max_issues_repo_name: bgamari/bayes-stack

max_forks_repo_name: bgamari/bayes-stack


text:
module Main where

import Quadratic
import Data.Complex

main :: IO ()
main = do
  putStrLn $ show $ roots (Quadratic 0 1 2)
  putStrLn $ show $ roots (Quadratic 1 3 4)
  putStrLn $ show $ roots (Quadratic 1 3 4)
  putStrLn $ show $ roots (Quadratic 1 4 4)
  putStrLn $ show $ roots (Quadratic 1 0 4)

max_stars_repo_name: subhashischatterjee/HuskellCookbook

max_issues_repo_name: subhashischatterjee/HuskellCookbook

max_forks_repo_name: subhashischatterjee/HuskellCookbook


text:
{-# LANGUAGE FlexibleContexts #-}
module Sound.Hommage.Sample
 (
 -- * Sample Datastructure
   Sample
 , arraySample
 , sizeSample
 , newSample
 , newSampleFromList

 -- * File Operations
 , openWavSample
 , saveWavSampleMono
 , saveWavSampleStereo

 -- * Pitching, Playing and Scratching
 , pitchWavSignal
 , scratchWavSignal
 , playSample
 , pitchSample
 , scratchSample
 , playCoeffs

 -- * Transforming Samples
 , mapSample
 , foldlSample
 , updateSample
 , normaliseSample
 , fadeinoutSample
 , stretchSample

 -- * Filter and Fourier Coeffs
 , Coeffs
 , arrayCoeffs
 , sizeCoeffs
 , wavesizeCoeffs
 , syntheseSample
 , analyseSample
 , filterCoeffs
 , mkPlayCoeffs
 , readCoeffs
 , writeCoeffs

 -- * Playing a Sample
 , mkTriggerSample
 , mkPlaySample
 , mkScratchSample
 , mkPitchSample

 -- * Array functions
 , interpolArray
 , mkLoopArray
 , mkPlayArray

 )
 where

import Foreign.Storable
import Data.Array.Storable
import Data.Array.IO
import Data.Int
import Data.IORef
import Data.Complex

import Sound.Hommage.Misc
import Sound.Hommage.WavFile
import Sound.Hommage.DFTFilter
import Sound.Hommage.Signal

---------------------------------------------------------------------------------------------------
-- | Open a Wav-File, but play it with variable (positive) speed\/ frequency.
pitchWavSignal :: FilePath -> [Double] -> Signal
pitchWavSignal fp = either Mono Stereo . inListE''
 (openWavSample fp >>= either (fmap Left . mkPitchSample) (fmap Right . mkPitchSample) )

-- | Open a Wav-File, but play it with variable (positive and\/or negative) speed\/ frequency.
--   Backward playing is possible.
scratchWavSignal :: FilePath -> [Double] -> Signal
scratchWavSignal fp = either Mono Stereo . inListE'
 (openWavSample fp >>= either (fmap Left . mkScratchSample) (fmap Right . mkScratchSample) )
---------------------------------------------------------------------------------------------------
-- | Creates a finite list with the content of the Sample.
playSample :: Sample a -> [a]
playSample = toList' . mkPlaySample

-- | Creates a list with the content of the Sample, but played with variable speed\/frequency.
pitchSample :: Fractional a => Sample a -> [Double] -> [a]
pitchSample s = inList'' $ mkPitchSample s

-- | Creates a list with the content of the Sample, but played with variable speed\/frequency,
--   backward playing possible.
scratchSample :: Fractional a => Sample a -> [Double] -> [a]
scratchSample s = inList' $ mkScratchSample s

-- | Creates a list with the content of 'Coeffs' (a representation of the fourier coefficients of
--   a 'Sample').
playCoeffs :: Coeffs -> [Complex Double]
playCoeffs = toList' . mkPlayCoeffs
---------------------------------------------------------------------------------------------------

---------------------------------------------------------------------------------------------------
--interpolArray :: MArray array Double IO => array Int Double -> Int -> Double -> IO Double
interpolArray :: (Fractional a, MArray array a IO) => array Int a -> Int -> Double -> IO a
interpolArray arr size =
 let f  = fromRational . toRational
 in \pos -> do
  let p0 = floor pos
      p1 = p0 + 1
      d :: Double
      d  = pos - fromIntegral p0
      d' = 1.0 - d
  v1 <- readArray arr $ mod p0 size
  v2 <- readArray arr $ mod p1 size
  let r = v1 * f d' + v2 * f d
  return r

mkLoopArray :: (Fractional a, MArray arr a IO) => arr Int a -> Int -> IO (Double -> IO a)
mkLoopArray arr size = do
 r <- newIORef 0.0
 let f  = fromRational . toRational
 return $ \dp ->
  do pos <- readIORef r
     let p  = floor pos
         p0 = mod p size
         p1 = mod (p + 1) size
         d :: Double
         d  = pos - fromIntegral p
         pos' = d + fromIntegral p0
         d' = 1.0 - d
     v1 <- readArray arr p0
     v2 <- readArray arr p1
     writeIORef r (pos' + dp)
     return (v1 * f d' + v2 * f d)

mkPlayArray :: (Fractional a, MArray arr a IO) => arr Int a -> Int -> IO (Double -> IO (Maybe a))
mkPlayArray arr size = do
 r <- newIORef 0.0
 let f  = fromRational . toRational
 return $ \dp ->
  do pos <- readIORef r
     let p0 = floor pos
         p1 = p0 + 1
         d :: Double
         d  = pos - fromIntegral p0
         d' = 1.0 - d
     if p1 >= size then return Nothing else do
      v1 <- readArray arr p0
      v2 <- readArray arr p1
      writeIORef r (pos + abs dp)
      return $ Just (v1 * f d' + v2 * f d)

---------------------------------------------------------------------------------------------------
data Sample a = SAMPLE
 { arraySample :: !(IOArray Int a)
 , sizeSample  :: !Int
 }
---------------------------------------------------------------------------------------------------
newSample :: Int -> a -> IO (Sample a)
newSample n a = newArray (0, n-1) a >>= \arr -> return $ SAMPLE arr n

newSampleFromList :: [a] -> IO (Sample a)
newSampleFromList xs = do
 let n = length xs
 arr <- newListArray (0, n - 1) xs
 return $ SAMPLE arr n

mapSample :: (a -> b) -> Sample a -> IO (Sample b)
mapSample f (SAMPLE arr n) = do
 brr <- newArray_ (0, n - 1)
 for 0 (<n) (+1) (\i -> readArray arr i >>= writeArray brr i . f)
 return $ SAMPLE brr n
---------------------------------------------------------------------------------------------------
openWavSample :: FilePath -> IO (Either (Sample Mono) (Sample Stereo))
openWavSample filepath = do
 (close, read, len, stereo) <- openInputWavFile filepath
 let bsize = 4000
 buf <- newArray_ (0, bsize-1)
 if stereo
  then do let lenh = div len 2
              hsize = div bsize 2
          arr <- newArray_ (0, lenh - 1)
          let step o i = do let i1=i*2                  -- i: maximal 1999
                            x1 <- readArray buf i1      --
                            x2 <- readArray buf (i1+1)
                            writeArray arr (o + i) (wavInt16ToDouble x1 :><: wavInt16ToDouble x2)
              loop o k | k < bsize = do read buf k
                                        for 0 (<(div k 2)) (+1) (step o)
                                        close
                                        return $ Right $ SAMPLE arr lenh
                       | otherwise = do read buf bsize
                                        for 0 (<hsize) (+1) (step o)
                                        loop (o+hsize) (k-bsize)
          loop 0 len
  else do arr <- newArray_ (0, len - 1)
          let loop o k | k < bsize = do read buf k
                                        for 0 (<k) (+1) (\i -> readArray buf i >>= writeArray arr (o + i) . wavInt16ToDouble )
                                        close
                                        return $ Left $ SAMPLE arr len
                       | otherwise = do read buf bsize
                                        for 0 (<bsize) (+1) (\i -> readArray buf i >>= writeArray arr (o + i) . wavInt16ToDouble )
                                        loop (o+bsize) (k-bsize)
          loop 0 len
---------------------------------------------------------------------------------------------------
saveWavSampleMono :: FilePath -> Sample Double -> IO ()
saveWavSampleMono filepath (SAMPLE arr len) = do
 (close, write) <- openOutputWavFileMono filepath
 let bsize = 4000
 buf <- newArray_ (0, bsize-1)
 let step o i = do x <- readArray arr (o+i)
                   writeArray buf i (wavDoubleToInt16 x)
     loop o k | k < bsize = do for 0 (<k) (+1) (step o)
                               write buf k
                               close
              | otherwise = do for 0 (<bsize) (+1) (step o)
                               write buf bsize
                               loop (o + bsize) (k - bsize)
 loop 0 len
---------------------------------------------------------------------------------------------------
saveWavSampleStereo :: FilePath -> Sample Stereo -> IO ()
saveWavSampleStereo filepath (SAMPLE arr len) = do
 (close, write) <- openOutputWavFileStereo filepath
 let bsize = 4000
     hsize = div bsize 2
 buf <- newArray_ (0, bsize-1)
 let step o i = do let i1 = i * 2
                       i2 = i1 + 1
                   x <- readArray arr (o+i)
                   writeArray buf i1 (wavDoubleToInt16 $ leftStereo x)
                   writeArray buf i2 (wavDoubleToInt16 $ rightStereo x)
     loop o k | k < hsize = do for 0 (<k) (+1) (step o)
                               write buf (2*k)
                               close
              | otherwise = do for 0 (<hsize) (+1) (step o)
                               write buf bsize
                               loop (o + hsize) (k - hsize)
 loop 0 len

--12000    348

---------------------------------------------------------------------------------------------------
mkTriggerSample :: Sample a -> a -> IO (Bool -> IO a)
mkTriggerSample s a = do
 rp <- newIORef 0
 let n = sizeSample s
     arr = arraySample s
     fun b = if b
              then readIORef rp >>= \p ->
                   if p < n
                    then writeIORef rp (p+1) >>
                         readArray arr p
                    else return a
              else writeIORef rp 0 >>
                   return a
 return fun

mkPlaySample :: Sample a -> IO (IO (Maybe a))
mkPlaySample (SAMPLE arr n) = do
 rpos <- newIORef 0
 let read = readIORef rpos >>= \pos ->
            if pos >= n then return Nothing
                        else writeIORef rpos (pos+1) >> readArray arr pos >>= return . Just
 return read
---------------------------------------------------------------------------------------------------
foldlSample :: Sample a -> (b -> a -> b) -> b -> IO b
foldlSample (SAMPLE arr n) f x =
 for' 0 (<n) (+1) x (\i x -> readArray arr i >>= \y -> return (f x y))

updateSample :: Sample a -> (a -> a) -> IO ()
updateSample (SAMPLE arr n) f =
 for 0 (<n) (+1) (\i -> readArray arr i >>= writeArray arr i . f)
---------------------------------------------------------------------------------------------------
normaliseSample :: (Ord a, Fractional a) => Sample a -> IO ()
normaliseSample s = do
 m <- foldlSample s (\b a -> max b (abs a)) 0.001
 updateSample s ((*)(1.0 / m))
---------------------------------------------------------------------------------------------------
fadeinoutSample :: Fractional a => (Int, Int) -> Sample a -> IO ()
fadeinoutSample (li,lo) sa = do
 let n = sizeSample sa
     arr = arraySample sa
     li' = min li n
     lo' = min lo n
     ol  = n - lo'
     si  = 1.0 / fromIntegral li'
     so  = 1.0 / fromIntegral lo'
 for' 0 (<li') (+1) 0.0 (\i s -> readArray arr i >>= \v -> writeArray arr i (v * s) >> return (s+si))
 for' ol (<n) (+1) 1.0 (\i s -> readArray arr i >>= \v -> writeArray arr i (v * s) >> return (s-so))
 return ()

---------------------------------------------------------------------------------------------------
stretchSample :: Fractional a => Int -> Sample a -> IO (Sample a)
stretchSample n sam = do
 let s   = sizeSample sam
     arr = arraySample sam
     d   = fromIntegral s / fromIntegral n
 ram <- newSample n 0.0
 for 0 (<n) (+1) (\i -> interpolArray arr s (d * fromIntegral i) >>= writeArray (arraySample ram) i)
 return ram

mkScratchSample :: Fractional a => Sample a -> IO (Double -> IO a)
mkScratchSample (SAMPLE arr n) = mkLoopArray arr n

mkPitchSample :: Fractional a => Sample a -> IO (Double -> IO (Maybe a))
mkPitchSample (SAMPLE arr n) = mkPlayArray arr n
---------------------------------------------------------------------------------------------------
data Coeffs = COEFFS
 { arrayCoeffs    :: IOArray Int CoeffArr
 , sizeCoeffs     :: Int
 , wavesizeCoeffs :: Int
 }

mkPlayCoeffs :: Coeffs -> IO (IO (Maybe (Complex Double)))
mkPlayCoeffs (COEFFS arr siz wsiz) = do
 r <- newIORef (0,0)
 let read = readIORef r >>= \(n,k) ->
            if n >= siz then return Nothing else
            if k >= 512 then writeIORef r (n+1,0) >> read else
            readArray arr n >>= \cs ->
            readCoeffArr cs k >>= \c ->
            writeIORef r (n,k+1) >>
            return (Just c)
 return read

analyseSample :: Sample Mono -> IO Coeffs
analyseSample s = do
 let arr = arraySample s
     siz   = sizeSample s
     n     = div siz 512
     r     = mod siz 512
     n'    = if r == 0 then n-1 else n
     n''   = n' - 1
     r'    = if r == 0 then 1024 else r+512
     p'    = n'' * 512
 dummy <- newArray_ (0,1023)
 brr <- newArray (0, n'') undefined
 for 0 (<n'') (+1) (\i -> do let p = i * 512
                             for 0 (<1024) (+1) (\i -> readArray arr (i+p) >>= writeArray dummy i)
                             crr <- newArray_ (0, 1023)
                             writeArray brr i crr
                             analyseDFT dummy crr)
 crr <- newArray_ (0, 1023)
 writeArray brr n'' crr
 for 0 (<r') (+1) (\i -> readArray arr (i+p') >>= writeArray dummy i)
 for r' (<1024) (+1) (\i -> writeArray dummy i 0.0)
 analyseDFT dummy crr
 return $ COEFFS brr n' siz

syntheseSample :: Coeffs -> IO (Sample Mono)
syntheseSample (COEFFS arr n siz) = do
 so <- newSample siz 0.0
 let n' = n-1
 dummy <- newArray (0, 1023) 0.0
-- sa :: StoraArray Int Double) <- newArray (0, 1023) 0.0
 for 0 (<n') (+1) (\i -> do crr <- readArray arr i
                            syntheseDFT crr dummy
                            kurveDFT dummy
                            let pl = i * 512
                            for 0 (<512) (+1) (\j -> let jpl = j + pl
                                                         jpr = jpl + 512
                                                     in
                                                     readArray (arraySample so) jpl >>= \x ->
                                                     readArray dummy j >>= \y ->
                                                     writeArray (arraySample so) jpl (x+y) >>
                                                     readArray dummy (j+512) >>= \z ->
                                                     writeArray (arraySample so) jpr z
                                              )
                  )
 crr <- readArray arr n'
 syntheseDFT crr dummy
 kurveDFT dummy
 let pl = n' * 512
     pm = siz - (pl + 512)
 for 0 (<512) (+1) (\j -> let jpl = j + pl
                          in
                          readArray (arraySample so) jpl >>= \x ->
                          readArray dummy j >>= \y ->
                          writeArray (arraySample so) jpl (x+y) )
 for 0 (<pm) (+1) (\j -> readArray dummy (j+512) >>= \z ->
                          writeArray (arraySample so) (j+pl+512) z )
 return so

filterCoeffs :: IO (CoeffMap ()) -> Coeffs -> IO Coeffs
filterCoeffs mkf (COEFFS arrci n siz) = do
 arrco <- newArray (0, n-1) undefined
 f   <- mkf
 for 0 (<n) (+1) (\i -> do ci <- readArray arrci i
                           co <- newArray_ (0, 1023)
                           f ci co
                           writeArray arrco i co)
 return $ COEFFS arrco n siz

-- | usage: @ readCoeffs obj nrOfCoeffarr nrOfCoefficient @
readCoeffs :: Coeffs -> Int -> Int -> IO (Complex Double)
readCoeffs cf a c = readArray (arrayCoeffs cf) a >>= \ca ->
                    readCoeffArr ca c

writeCoeffs :: Coeffs -> Int -> Int -> (Complex Double) -> IO ()
writeCoeffs cf a c v = readArray (arrayCoeffs cf) a >>= \ca ->
                       writeCoeffArr ca c v

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

max_stars_repo_name: aische/hommage

max_issues_repo_name: aische/hommage

max_forks_repo_name: aische/hommage