text:
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Data.Network.Utils
 (repeatM, scanl, randomInput, rowNorm, rowNorm', (.*), (*.),
  Differentiable, derivate, squeeze, squish, toGaussian,
  makeBatchP, unzipP, integersP, averageP,
  accuracy, accuracyBatch)
where

import Prelude hiding (scanl, zip)
import Data.Bifunctor as Bf
import Data.Foldable as F
import Numeric.LinearAlgebra as LA
import Control.Monad.Random as R
import qualified Data.Vector.Storable as V
import Control.Applicative as A
import Pipes as P
import Pipes.Core as PC
import qualified Pipes.Prelude as PP

scanl :: Foldable f => (b -> a -> b) -> b -> f a -> [b]
scanl f b a = reverse $ foldl (\(i:is) x -> (f i x):i:is) [b] a

scanr :: Foldable f => (a -> b -> b) -> b -> f a -> [b]
scanr f b a = foldr (\x (i:is) -> (f x i):i:is) [b] a

randomInput :: MonadRandom m => RandDist -> Int -> m (LA.Vector Double)
randomInput dist i = do
  seed :: Int <- getRandom
  return $ LA.randomVector seed dist i

a .* b = LA.scale a b
a *. b = LA.scale b a

rowNorm :: Matrix Double -> Vector R
rowNorm m = squeeze $ sqrt (m**2)

rowNorm' :: (Element a, Normed (Vector a)) => Matrix a -> Vector R
rowNorm' = LA.fromList . (fmap norm_2) . toRows

class Differentiable a b c | c -> a b where
  derivate :: c -> (a -> b)

repeatM :: Monad m => Int -> (a -> m a) -> a -> m a
repeatM n f i
  | n <= 0 = return i
  | otherwise = do
      o <- f i
      repeatM (n-1) f o

-- flatten laterally(summing over rows)
squeeze :: (V.Storable a, Num a, Numeric a) => Matrix a -> Vector a
squeeze m = m #> (V.generate (cols m) (const 1))

-- flatten vertically(summing over columns)
squish :: (V.Storable a, Num a, Numeric a) => Matrix a -> Vector a
squish m = (V.generate (rows m) (const 1)) <# m

--- Pipe utils
-----------------------------------------------------------------------------------------

-- Collect chunks of values of a particular size and yields those chunks as lists
-- Note: the elements of each list are in reverse order of arrival.
-- Compose with (PP.map reverse) to change that.
chunkP :: (Monad m) => Int -> Pipe a [a] m r
chunkP n = collect n [] where
  collect 0 xs = (yield xs) >> (collect n [])
  collect n xs = do
    x <- await
    collect (n-1) (x:xs)

-- transforms a stream of vectors into a stream of matrix(i.e. batch of vectors)
makeBatchP :: (Element a, Monad m) => Int -> Pipe (Vector a) (Matrix a) m r
makeBatchP n = chunkP n >-> PP.map (LA.fromRows)

unzipP :: (Monad m) => Producer (a,b) m r -> (Producer a m r, Producer b m r)
unzipP p = (p >-> (PP.map fst), p >-> (PP.map snd))

integersP :: (Monad m, Enum a) => a -> Producer a m r
integersP n = yield n >> integersP (succ n)

rangeP :: (Monad m, Enum a, Ord a) => a -> a -> Producer a m ()
rangeP !a !b
  | a < b = yield a >> rangeP (succ a) b
  | otherwise = return ()

whileMaybe :: (Monad m) => Proxy a' a b' b m r -> Proxy a' a b' (Maybe b) m r
whileMaybe p = for p (respond . Just) <* respond Nothing

averageP :: (Monad m, Num a, Floating a) => Producer a m r -> m (a, r)
averageP = PP.fold' (\x r -> bimap (+r) (+1) x) (0,0) (uncurry (/))

-----------------------------------------------------------------------------------
-- Random stream transformation
-----------------------------------------------------------------------------------

toGaussian :: [Double] -> [Double]
toGaussian (a:b:cs) = (r * cos θ):(r * sin θ):toGaussian cs
  where r = sqrt (-2 * log a)
        θ = 2*pi*b

----------------------------------------
accuracy :: [(Vector R, Vector R)] -> Double
accuracy outs = (/ (fromIntegral $ length outs)) $ fromIntegral $ length $ filter (\(y,y') -> (maxIndex y == maxIndex y')) outs

accuracyBatch :: (Matrix R, Matrix R) -> Double
accuracyBatch out = let n = fromIntegral $ rows $ fst out
                        hits = let (idx, idx') = bimap (fmap maxIndex . toRows) (fmap maxIndex . toRows) out
                               in fromIntegral $ length $ filter id (zipWith (==) idx idx')
                     in  hits/n

hitBatchP :: (Floating a, Monad m) => Pipe (Matrix R, Matrix R) a m ()
hitBatchP = do
  (o,t) <- await
  let (idx, idx') = bimap (fmap maxIndex . toRows) (fmap maxIndex . toRows) (o,t)
      hits = filter id (zipWith (==) idx idx')
  yield $ (fromIntegral $ length hits) / (fromIntegral $ rows o)

max_stars_repo_name: DrPyser/NeuralNetwork

max_issues_repo_name: DrPyser/NeuralNetwork

max_forks_repo_name: DrPyser/NeuralNetwork


text:
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ExplicitNamespaces #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE InstanceSigs #-}
{-# LANGUAGE RecordWildCards #-}
module Data.Matrix.Static.LinearAlgebra
    ( module Data.Matrix.Static.LinearAlgebra.Types
    , Arithmetic(..)
    , Factorization(..)
    , LinearAlgebra(..)
    , EigenOptions(..)
    , defaultEigenOptions
    , SortRule(..)

      -- * Dense matrix operation
    , zeros
    , ones
    , inverse
    , eig
    , svd
    , cond
    ) where

import qualified Data.Vector.Storable as VS
import Data.Vector.Storable (Vector)
import System.IO.Unsafe (unsafePerformIO)
import Data.Complex (Complex)
import Data.Singletons.Prelude hiding ((@@), type (==), type (-), type (<=))
import Data.Type.Equality
import GHC.TypeLits
import Data.Maybe (fromMaybe)

import qualified Data.Matrix.Static.Dense as D
import qualified Data.Matrix.Static.Sparse as S
import qualified Data.Matrix.Static.Generic.Mutable as CM
import qualified Data.Matrix.Static.Generic as C
import qualified Data.Matrix.Static.Internal as Internal
import Data.Matrix.Static.LinearAlgebra.Types
import Data.Matrix.Static.LinearAlgebra.Internal

class Arithmetic (mat1 :: C.MatrixKind) (mat2 :: C.MatrixKind) where
    -- | Matrix multiplication between different types of matrices.
    (@@) :: ( Numeric a, SingI n, SingI m
            , If (mat1 == mat2) mat1 D.Matrix ~ mat3 )
         => mat1 n p Vector a
         -> mat2 p m Vector a
         -> mat3 n m Vector a
    infixr 8 @@

    -- | Element-wise addition between different types of matrices.
    (%+%) :: ( Numeric a, SingI n, SingI m
             , If (mat1 == mat2) mat1 D.Matrix ~ mat3 )
          => mat1 n m Vector a
          -> mat2 n m Vector a
          -> mat3 n m Vector a
    infixr 8 %+%

    -- | Element-wise substraction between different types of matrices.
    (%-%) :: ( Numeric a, SingI n, SingI m
             , If (mat1 == mat2) mat1 D.Matrix ~ mat3 )
          => mat1 n m Vector a
          -> mat2 n m Vector a
          -> mat3 n m Vector a
    infixr 8 %-%

    -- | Element-wise multiplication between different types of matrices.
    (%*%) :: ( Numeric a, SingI n, SingI m
             , If (mat1 == mat2) mat1 S.SparseMatrix ~ mat3 )
          => mat1 n m Vector a
          -> mat2 n m Vector a
          -> mat3 n m Vector a
    infixr 8 %*%

instance Arithmetic D.Matrix D.Matrix where
    (@@) = withFun2 Internal.c_dd_mul
    (%+%) = (+)
    (%-%) = (-)
    (%*%) = (*)

instance Arithmetic D.Matrix S.SparseMatrix where
    (@@) = withDS Internal.c_ds_mul
    (%+%) = flip (%+%)
    (%-%) a b = a %+% C.map negate b
    (%*%) = undefined

instance Arithmetic S.SparseMatrix D.Matrix where
    (@@) = withSD Internal.c_sd_mul
    (%+%) = withSD Internal.c_sd_plus
    (%-%) a b = a %+% C.map negate b
    (%*%) = undefined

instance Arithmetic S.SparseMatrix S.SparseMatrix where
    (@@) = withSS Internal.c_ss_mul
    (%+%) = withSS Internal.c_ss_plus
    (%-%) a b = a %+% C.map negate b
    (%*%) = withSS Internal.c_ss_cmul

class LinearAlgebra (mat :: C.MatrixKind) where
    ident :: (Numeric a, SingI n) => mat n n Vector a

    colSum :: (Numeric a, SingI n, C.Matrix mat Vector a)
           => mat m n Vector a
           -> Matrix 1 n a
    colSum mat = D.create $ do
        m <- CM.replicate 0
        flip C.imapM_ mat $ \(_,j) v -> CM.unsafeModify m (+v) (0, j)
        return m
    {-# INLINE colSum #-}

    rowSum :: (Numeric a, SingI m, C.Matrix mat Vector a)
           => mat m n Vector a
           -> Matrix m 1 a
    rowSum mat = D.create $ do
        m <- CM.replicate 0
        flip C.imapM_ mat $ \(i,_) x -> CM.unsafeModify m (+x) (i, 0)
        return m
    {-# INLINE rowSum #-}

instance LinearAlgebra D.Matrix where
    ident = D.diag 0 $ D.replicate 1

instance LinearAlgebra S.SparseMatrix where
    ident = S.diag $ D.replicate 1

data EigenOptions = EigenOptions
    { _sigma :: Maybe Double -- ^ Find eigenvalues near sigma using shift-invert mode.
    , _ncv :: Maybe Int      -- ^ Parameter that controls the convergence
                             -- speed of the algorithm. Typically a larger ncv
                             -- means faster convergence, but it may also result
                             -- in greater memory use and more matrix operations
                             -- in each iteration. This parameter must satisfy
                             -- k+2 <= ncv <= n, and is advised to take ncv > 2k+1.
                             -- Default: min(n, max(2*k + 1, 20))
    , _maxit :: Maybe Int          -- ^ Maximum number of iterations allowed in the algorithm.
    , _tol   :: Double       -- Precision parameter for the calculated eigenvalues.
    , _sort_rule :: SortRule
    }

defaultEigenOptions :: EigenOptions
defaultEigenOptions = EigenOptions
    { _sigma = Nothing
    , _ncv = Nothing
    , _maxit = Nothing
    , _tol = 2.2204460492503131e-16
    , _sort_rule = LM }

data SortRule = LM
              | LR
              | LI
              | LA
              | SM
              | SR
              | SI
              | SA
              | BE
              deriving (Enum)

class Factorization mat where
    -- | Eigenvalues (from largest to smallest) and
    -- eigenvectors (as columns) of a general square matrix.
    eigS :: (SingI k, SingI n, k <= n - 2)
         => Sing k
         -> mat n n Vector Double
         -> (Matrix k 1 (Complex Double), Matrix n k (Complex Double))
    eigS a b = eigS' a b defaultEigenOptions

    eigS' :: (SingI k, SingI n, k <= n - 2)
         => Sing k
         -> mat n n Vector Double
         -> EigenOptions
         -> (Matrix k 1 (Complex Double), Matrix n k (Complex Double))

    -- | Eigenvalues (from largest to smallest) and
    -- eigenvectors (as columns) of a symmetric square matrix.
    eigSH :: (SingI k, SingI n, k <= n - 1)
          => Sing k
          -> mat n n Vector Double
          -> (Matrix k 1 Double, Matrix n k Double)
    eigSH a b = eigSH' a b defaultEigenOptions

    eigSH' :: (SingI k, SingI n, k <= n - 1)
           => Sing k
           -> mat n n Vector Double
           -> EigenOptions
           -> (Matrix k 1 Double, Matrix n k Double)

    -- | Generalized eigen solver for real symmetric matrices
    geigSH :: (SingI k, SingI n, k <= n - 1)
           => Sing k
           -> mat n n Vector Double
           -> SparseMatrix n n Double
           -> (Matrix k 1 Double, Matrix n k Double)
    geigSH a b c = geigSH' a b c defaultEigenOptions

    geigSH' :: (SingI k, SingI n, k <= n - 1)
            => Sing k
            -> mat n n Vector Double
            -> SparseMatrix n n Double
            -> EigenOptions
            -> (Matrix k 1 Double, Matrix n k Double)

    -- | Cholesky decomposition
    cholesky :: (Numeric a, SingI n) => mat n n Vector a -> mat n n Vector a

instance Factorization D.Matrix where
    eigS' s mat EigenOptions{..}
        | D.all (==0) mat = ( D.replicate 0, D.replicate 1)
        | otherwise = unsafePerformIO $ do
            m1 <- CM.new
            m2 <- CM.new
            info <- fmap Internal.computationInfo $ unsafeWith' m1 $ \v1 _ _ ->
                unsafeWith' m2 $ \v2 _ _ -> unsafeWith mat $ \v n _ ->
                    Internal.c_eigs k v1 v2 v n ncv' maxit _tol mode sigma $ fromIntegral $ fromEnum _sort_rule
            case info of
                Nothing -> do
                    m1' <- C.unsafeFreeze m1
                    m2' <- C.unsafeFreeze m2
                    return (m1', m2')
                Just err -> error err
      where
        ncv' = case _ncv of
            Just x -> fromIntegral x
            Nothing -> min d $ max 20 $ 2 * k + 1
        maxit = case _maxit of
            Just x -> fromIntegral x
            Nothing -> d * 10
        mode = case _sigma of
            Nothing -> 0
            _ -> 1
        sigma = fromMaybe 0 _sigma
        k = fromIntegral $ fromSing s
        d = fromIntegral $ D.rows mat
    {-# INLINE eigS' #-}

    eigSH' s mat EigenOptions{..}
        | D.all (==0) mat = (D.replicate 0, D.replicate 1)
        | otherwise = unsafePerformIO $ do
            m1 <- CM.new
            m2 <- CM.new
            ( fmap Internal.computationInfo $ unsafeWith' m1 $ \v1 _ _ ->
                unsafeWith' m2 $ \v2 _ _ -> unsafeWith mat $ \v n _ ->
                    Internal.c_eigsh k v1 v2 v n ncv' maxit _tol mode sigma $
                    fromIntegral $ fromEnum _sort_rule ) >>= maybe
                    ( do
                        m1' <- C.unsafeFreeze m1
                        m2' <- C.unsafeFreeze m2
                        return (m1', m2') ) error
      where
        ncv' = case _ncv of
            Just x -> fromIntegral x
            Nothing -> min d $ max 20 $ 2 * k + 1
        maxit = case _maxit of
            Just x -> fromIntegral x
            Nothing -> d * 10
        mode = case _sigma of
            Nothing -> 0
            _ -> 1
        sigma = fromMaybe 0 _sigma
        k = fromIntegral $ fromSing s
        d = fromIntegral $ D.rows mat
    {-# INLINE eigSH' #-}

    geigSH' s matA matB EigenOptions{..} = unsafePerformIO $ do
        m1 <- CM.new
        m2 <- CM.new
        ( fmap Internal.computationInfo $ unsafeWith' m1 $ \v1 _ _ ->
            unsafeWith' m2 $ \v2 _ _ -> unsafeWith matA $ \v n _ ->
            unsafeWithS matB $ \pv pin po _ _ size ->
            Internal.c_geigsh k v1 v2 v n pv po pin size
            ncv' maxit _tol $ fromIntegral $ fromEnum _sort_rule ) >>= maybe
            ( do
                m1' <- C.unsafeFreeze m1
                m2' <- C.unsafeFreeze m2
                return (m1', m2') ) error
      where
        ncv' = case _ncv of
            Just x -> fromIntegral x
            Nothing -> min d $ max 20 $ 2 * k + 1
        maxit = case _maxit of
            Just x -> fromIntegral x
            Nothing -> d * 10
        k = fromIntegral $ fromSing s
        d = fromIntegral $ D.rows matA
    {-# INLINE geigSH' #-}

    cholesky mat = flip withFun1 mat $
        \code p1 c1 _ p2 _ _ -> Internal.c_cholesky code p1 p2 c1
    {-# INLINE cholesky #-}

instance Factorization S.SparseMatrix where
    eigS' s mat EigenOptions{..} = unsafePerformIO $ do
        m1 <- CM.new
        m2 <- CM.new
        ( fmap Internal.computationInfo $ unsafeWith' m1 $ \v1 _ _ ->
            unsafeWith' m2 $ \v2 _ _ -> unsafeWithS mat $ \pv pin po n _ size ->
                Internal.c_seigs k v1 v2 pv po pin n size
                ncv' maxit _tol mode sigma $ fromIntegral $ fromEnum _sort_rule
                ) >>= maybe
            ( do
                m1' <- C.unsafeFreeze m1
                m2' <- C.unsafeFreeze m2
                return (m1', m2') ) error
      where
        ncv' = case _ncv of
            Just x -> fromIntegral x
            Nothing -> min d $ max 20 $ 2 * k + 1
        maxit = case _maxit of
            Just x -> fromIntegral x
            Nothing -> d * 10
        mode = case _sigma of
            Nothing -> 0
            _ -> 1
        sigma = fromMaybe 0 _sigma
        k = fromIntegral $ fromSing s
        d = fromIntegral $ D.rows mat
    {-# INLINE eigS' #-}

    eigSH' s mat EigenOptions{..} = unsafePerformIO $ do
        m1 <- CM.new
        m2 <- CM.new
        ( fmap Internal.computationInfo $ unsafeWith' m1 $ \v1 _ _ ->
            unsafeWith' m2 $ \v2 _ _ -> unsafeWithS mat $ \pv pin po n _ size ->
                Internal.c_seigsh k v1 v2 pv po pin n size
                ncv' maxit _tol mode sigma $ fromIntegral $ fromEnum _sort_rule
                ) >>= maybe
                ( do
                    m1' <- C.unsafeFreeze m1
                    m2' <- C.unsafeFreeze m2
                    return (m1', m2') ) error
      where
        ncv' = case _ncv of
            Just x -> fromIntegral x
            Nothing -> min d $ max 20 $ 2 * k + 1
        maxit = case _maxit of
            Just x -> fromIntegral x
            Nothing -> d * 10
        mode = case _sigma of
            Nothing -> 0
            _ -> 1
        sigma = fromMaybe 0 _sigma
        k = fromIntegral $ fromSing s
        d = fromIntegral $ D.rows mat
    {-# INLINE eigSH' #-}

    cholesky = undefined

type family R a where
    R Float = Float
    R Double = Double
    R (Complex Double) = Double
    R (Complex Float) = Float

zeros :: (SingI m, SingI n) => Matrix m n Double
zeros = D.replicate 0
{-# INLINE zeros #-}

ones :: (SingI m, SingI n) => Matrix m n Double
ones = D.replicate 1
{-# INLINE ones #-}

-- | The inverse of a dense matrix.
inverse :: (SingI n, Numeric a) => Matrix n n a -> Matrix n n a
inverse = withFun1 Internal.c_inverse
{-# INLINE inverse #-}

-- | Compute the full eigendecomposition for dense matrix.
eig :: forall n . SingI n
    => Matrix n n Double
    -> (Matrix n 1 (Complex Double), Matrix n n (Complex Double))
eig mat = unsafePerformIO $ do
    m1 <- CM.new
    m2 <- CM.new
    _ <- unsafeWith' m1 $ \v1 _ _ -> unsafeWith' m2 $ \v2 _ _ ->
        unsafeWith mat $ \v n _ -> Internal.c_eig v1 v2 v n
    m1' <- C.unsafeFreeze m1
    m2' <- C.unsafeFreeze m2
    return (m1', m2')
{-# INLINE eig #-}

-- | Compute the full singular value decomposition for dense matrix.
svd :: forall n p a m. (Numeric (R a), Numeric a, SingI n, SingI p, SingI m, m ~ Min n p)
    => Matrix n p a
    -> (Matrix n m a, Matrix m 1 (R a), Matrix p m a)
svd mat = unsafePerformIO $ do
    mu <- CM.new
    ms <- CM.new
    mv <- CM.new
    checkResult $ unsafeWith' mu $ \pu _ _ -> unsafeWith' ms $ \ps _ _ ->
        unsafeWith' mv $ \pv _ _ -> unsafeWith mat $ \px r c ->
            Internal.c_bdcsvd (foreignType (undefined :: a))
                pu ps pv px r c
    u <- C.unsafeFreeze mu
    s <- C.unsafeFreeze ms
    v <- C.unsafeFreeze mv
    return (u, s, v)
{-# INLINE svd #-}

-- | Condition number.
cond :: ( Numeric a, Numeric (R a), Ord (R a), Fractional (R a)
        , SingI n, SingI m, SingI (Min n m))
     => Matrix n m a -> R a
cond mat = VS.maximum val / VS.minimum val
  where
    val = VS.filter (/=0) $ D.flatten s
    (_,s,_) = svd mat
{-# INLINE cond #-}

ext: hs

max_stars_repo_name: kaizhang/matrix-sized

max_issues_repo_name: kaizhang/matrix-sized

max_forks_repo_name: kaizhang/matrix-sized


text:
module PCA.QModels where

import PCA.Distribution
import Numeric.LinearAlgebra.Data
       (Matrix, Vector, ident, konst, size, toRows, fromRows, diag, diagl,
        takeDiag)
import qualified Numeric.LinearAlgebra.HMatrix as H

initQtau :: Distr Double Double
initQtau = Gamma 1.0e-3 1.0e-3

initQmu :: Int -> Distr (Vector Double) (Matrix Double)
initQmu d = MNormal (konst 0 d) (ident d)

initQalpha :: Int -> Distr (Vector Double) (Vector Double)
initQalpha d = MGamma (konst 1.0e-3 d) (konst 1.0e-3 d)

initQW :: Int -> Distr (Matrix Double) (Matrix Double)
initQW d = MatrixMNormal (ident d) (ident d)

initQX :: Matrix Double -> Distr (Matrix Double) (Matrix Double)
initQX trainT = let (_, d) = size trainT
                    sigmaX = calcSigmaX d initQtau (initQW d)
                    mX = calcMX trainT initQtau sigmaX (initQW d) (initQmu d)
                in MatrixMNormal mX sigmaX

type QModel = (Matrix Double, Vector Double)  -- ^ Mode of W and \mu

calculateQ :: Matrix Double -> QModel
calculateQ trainT =
    go 5 initQtau (initQmu d) (initQalpha d) (initQW d) (initQX trainT) trainT
  where
    (n,d) = size trainT
    go 0 _ mu _ w _ _ = (mean w, mean mu)
    go k (Gamma aTau bTau) mu (MGamma aAlpha bAlpha) w x t =
        go (k - 1) ntau nmu nalpha nw nx t
      where
        ntau = Gamma (calcAtau n d aTau) (calcBtau bTau t mu w x)
        nalpha = MGamma (calcAalpha d aAlpha) (calcBalpha bAlpha w)
        nSigmaMu = calcSigmaMu n ntau mu
        nmu = MNormal (calcMMu t ntau nSigmaMu w x) nSigmaMu
        nSigmaW = calcSigmaW ntau nalpha x
        nw = MatrixMNormal (calcMW ntau nSigmaW x t nmu) nSigmaW
        nSigmaX = calcSigmaX d ntau nw
        nx = MatrixMNormal (calcMX t ntau nSigmaX nw nmu) nSigmaX

calcSigmaX
    :: Int
    -> Distr Double Double
    -> Distr (Matrix Double) (Matrix Double)
    -> Matrix Double
calcSigmaX d tau w = H.inv (ident d + mean tau `H.scale` wSquareMean)
  where
    mW = mean w
    wSquareMean = variance w + (H.tr mW H.<> mW)

calcMX
    :: Matrix Double
    -> Distr Double Double
    -> Matrix Double
    -> Distr (Matrix Double) (Matrix Double)
    -> Distr (Vector Double) (Matrix Double)
    -> Matrix Double
calcMX trainT tau sigmaX w mu =
    H.tr
        (mean tau `H.scale` sigmaX H.<> H.tr (mean w) H.<>
         H.tr
             (fromRows
                  (map
                       (\r ->
                             r - mean mu)
                       (toRows trainT))))

calcSigmaMu
    :: Int
    -> Distr Double Double
    -> Distr (Vector Double) (Matrix Double)
    -> Matrix Double
calcSigmaMu n tau mu =
    vMu + 1 / (fromIntegral n * mean tau) `H.scale` ident d
  where
    vMu = variance mu
    (_,d) = size vMu

calcMMu
    :: Matrix Double
    -> Distr Double Double
    -> Matrix Double
    -> Distr (Matrix Double) (Matrix Double)
    -> Distr (Matrix Double) (Matrix Double)
    -> Vector Double
calcMMu trainT tau sigmaMu w x =
    (mean tau `H.scale` sigmaMu) H.#>
    foldr
         (flip (+))
         (konst 0 d)
         (toRows (trainT - H.tr (mW H.<> H.tr mX)))
  where
    mW = mean w
    mX = mean x
    (_, d) = size trainT

calcSigmaW
    :: Distr Double Double
    -> Distr (Vector Double) (Vector Double)
    -> Distr (Matrix Double) (Matrix Double)
    -> Matrix Double
calcSigmaW tau alpha x =
    H.inv
        (diag (mean alpha) +
         mean tau `H.scale`
         (fromIntegral n `H.scale` variance x + H.tr mX H.<> mX))
  where
    mX = mean x
    (n,_) = size mX

calcMW
    :: Distr Double Double
    -> Matrix Double
    -> Distr (Matrix Double) (Matrix Double)
    -> Matrix Double
    -> Distr (Vector Double) (Matrix Double)
    -> Matrix Double
calcMW tau sigmaW x trainT mu =
    (mean tau `H.scale` sigmaW) H.<>
    (H.tr mX H.<>
     (fromRows
          (map
               (\tn ->
                     tn - mMu)
               (toRows trainT))))
  where
    mX = mean x
    mMu = mean mu

calcAalpha :: Int -> Vector Double -> Vector Double
calcAalpha d = H.cmap (+ fromIntegral d / 2)

calcBalpha :: Vector Double
           -> Distr (Matrix Double) (Matrix Double)
           -> Vector Double
calcBalpha b w = b + 0.5 `H.scale` (takeDiag (H.tr mW H.<> mW))
  where
    mW = mean w

calcAtau :: Int -> Int -> Double -> Double
calcAtau n d a = a + fromIntegral (n * d) / 2

-- | TODO: Formula is not correct.
-- Skip part Tr (\langle W^TW \rangle \langle x_n x^T_n 'rangle')
-- because i think it is not correct in that result must be number
calcBtau
    :: Double
    -> Matrix Double
    -> Distr (Vector Double) (Matrix Double)
    -> Distr (Matrix Double) (Matrix Double)
    -> Distr (Matrix Double) (Matrix Double)
    -> Double
calcBtau b trainT mu w x =
    b +
    0.5 *
    (fromIntegral n * norm mMu +
     H.sumElements (takeDiag (trainT H.<> H.tr trainT)) +
     2 * (H.sumElements (H.tr (mW H.<> H.tr mX) H.#> mMu)) -
     2 * (H.sumElements (takeDiag (trainT H.<> mW H.<> H.tr mX))) -
     2 * (H.sumElements (trainT H.#> mMu)))
  where
    (n,_) = size trainT
    norm v = v `H.dot` v
    mMu = mean mu
    mW = mean w
    mX = mean x

max_stars_repo_name: DbIHbKA/vbpca

max_issues_repo_name: DbIHbKA/vbpca

max_forks_repo_name: DbIHbKA/vbpca


text:
{-# LANGUAGE CPP #-}
module Driver (
    E,
    V,
    M,

    AEq,
    (===),
    (~==),

    field,

    mytest,
    mycheck,
    mytests,
    done,

    ) where

import Data.AEq( AEq )
import qualified Data.AEq as AEq
import Data.Complex
import Data.List
import Debug.Trace
import System.IO
import System.Random
import Test.QuickCheck
import Text.Printf
import Text.Show.Functions

import Data.Vector.Dense( Vector )
import Data.Matrix.Dense( Matrix )

#ifdef COMPLEX
field = "Complex Double"
type E = Complex Double
#else
field = "Double"
type E = Double
#endif

type V = Vector Int E
type M = Matrix (Int,Int) E

infix 4 ===, ~==

x === y | (AEq.===) x y = True
        | otherwise = trace (printf "expected `%s', but got `%s'" (show y) (show x)) False

x ~== y | (AEq.~==) x y = True
        | otherwise = trace (printf "expected `%s', but got `%s'" (show y) (show x)) False

------------------------------------------------------------------------
--
-- QC driver ( taken from xmonad-0.6 )
--

debug = False

mytest :: Testable a => a -> Int -> IO (Bool, Int)
mytest a n = mycheck defaultConfig
    { configMaxTest=n
    , configEvery   = \n args -> let s = show n in s ++ [ '\b' | _ <- s ] } a
 -- , configEvery= \n args -> if debug then show n ++ ":\n" ++ unlines args else [] } a

mycheck :: Testable a => Config -> a -> IO (Bool, Int)
mycheck config a = do
    rnd <- newStdGen
    mytests config (evaluate a) rnd 0 0 []

mytests :: Config -> Gen Result -> StdGen -> Int -> Int -> [[String]] -> IO (Bool, Int)
mytests config gen rnd0 ntest nfail stamps
    | ntest == configMaxTest config = done "OK," ntest stamps >> return (True, ntest)
    | nfail == configMaxFail config = done "Arguments exhausted after" ntest stamps >> return (True, ntest)
    | otherwise               =
      do putStr (configEvery config ntest (arguments result)) >> hFlush stdout
         case ok result of
           Nothing    ->
             mytests config gen rnd1 ntest (nfail+1) stamps
           Just True  ->
             mytests config gen rnd1 (ntest+1) nfail (stamp result:stamps)
           Just False ->
             putStr ( "Falsifiable after "
                   ++ show ntest
                   ++ " tests:\n"
                   ++ unlines (arguments result)
                    ) >> hFlush stdout >> return (False, ntest)
     where
      result      = generate (configSize config ntest) rnd2 gen
      (rnd1,rnd2) = split rnd0

done :: String -> Int -> [[String]] -> IO ()
done mesg ntest stamps = putStr ( mesg ++ " " ++ show ntest ++ " tests" ++ table )
  where
    table = display
            . map entry
            . reverse
            . sort
            . map pairLength
            . group
            . sort
            . filter (not . null)
            $ stamps

    display []  = ".\n"
    display [x] = " (" ++ x ++ ").\n"
    display xs  = ".\n" ++ unlines (map (++ ".") xs)

    pairLength xss@(xs:_) = (length xss, xs)
    entry (n, xs)         = percentage n ntest
                       ++ " "
                       ++ concat (intersperse ", " xs)

    percentage n m        = show ((100 * n) `div` m) ++ "%"

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

max_stars_repo_name: patperry/lapack

max_issues_repo_name: patperry/lapack

max_forks_repo_name: patperry/lapack


text:
{-# LANGUAGE CPP, TypeOperators #-}

#if __GLASGOW_HASKELL__ >= 702
{-# LANGUAGE Trustworthy #-}
#endif

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

-----------------------------------------------------------------------------
-- |
-- Copyright   :  (C) 2011-2015 Edward Kmett
-- License     :  BSD-style (see the file LICENSE)
--
-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
-- Stability   :  provisional
-- Portability :  portable
--
----------------------------------------------------------------------------
module Data.Semigroup.Foldable.Class
  ( Foldable1(..)
  , Bifoldable1(..)
  ) where

import Control.Applicative
import Control.Applicative.Backwards
import Control.Applicative.Lift
import Control.Monad.Trans.Identity
import Data.Bifoldable
import Data.Bifunctor.Biff
import Data.Bifunctor.Clown
import Data.Bifunctor.Flip
import Data.Bifunctor.Join
import Data.Bifunctor.Product as Bifunctor
import Data.Bifunctor.Joker
import Data.Bifunctor.Tannen
import Data.Bifunctor.Wrapped
import Data.Foldable

import Data.Functor.Identity
import Data.Functor.Product as Functor
import Data.Functor.Reverse
import Data.Functor.Sum as Functor
import Data.Functor.Compose
import Data.List.NonEmpty (NonEmpty(..))

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

#ifdef MIN_VERSION_tagged
import Data.Tagged
#endif

import Data.Traversable.Instances ()

#ifdef MIN_VERSION_containers
import Data.Tree
#endif

import Data.Monoid as Monoid hiding ((<>))
import Data.Semigroup as Semigroup hiding (Product, Sum)
import Data.Orphans ()
-- import Data.Ord -- missing Foldable, https://ghc.haskell.org/trac/ghc/ticket/15098#ticket

#ifdef MIN_VERSION_generic_deriving
import Generics.Deriving.Base
#else
import GHC.Generics
#endif

import Prelude hiding (foldr)

class Foldable t => Foldable1 t where
  fold1 :: Semigroup m => t m -> m
  foldMap1 :: Semigroup m => (a -> m) -> t a -> m
  toNonEmpty :: t a -> NonEmpty a

  foldMap1 f = maybe (error "foldMap1") id . getOption . foldMap (Option . Just . f)
  fold1 = foldMap1 id
  toNonEmpty = foldMap1 (:|[])

instance Foldable1 Monoid.Sum where
  foldMap1 f (Monoid.Sum a) = f a

instance Foldable1 Monoid.Product where
  foldMap1 f (Monoid.Product a) = f a

instance Foldable1 Monoid.Dual where
  foldMap1 f (Monoid.Dual a) = f a

#if MIN_VERSION_base(4,8,0)
instance Foldable1 f => Foldable1 (Monoid.Alt f) where
  foldMap1 g (Alt m) = foldMap1 g m
#endif

instance Foldable1 Semigroup.First where
  foldMap1 f (Semigroup.First a) = f a

instance Foldable1 Semigroup.Last where
  foldMap1 f (Semigroup.Last a) = f a

instance Foldable1 Semigroup.Min where
  foldMap1 f (Semigroup.Min a) = f a

instance Foldable1 Semigroup.Max where
  foldMap1 f (Semigroup.Max a) = f a

instance Foldable1 f => Foldable1 (Rec1 f) where
  foldMap1 f (Rec1 as) = foldMap1 f as

instance Foldable1 f => Foldable1 (M1 i c f) where
  foldMap1 f (M1 as) = foldMap1 f as

instance Foldable1 Par1 where
  foldMap1 f (Par1 a) = f a

instance (Foldable1 f, Foldable1 g) => Foldable1 (f :*: g) where
  foldMap1 f (as :*: bs) = foldMap1 f as <> foldMap1 f bs

instance (Foldable1 f, Foldable1 g) => Foldable1 (f :+: g) where
  foldMap1 f (L1 as) = foldMap1 f as
  foldMap1 f (R1 bs) = foldMap1 f bs

instance Foldable1 V1 where
  foldMap1 _ v = v `seq` undefined

instance (Foldable1 f, Foldable1 g) => Foldable1 (f :.: g) where
  foldMap1 f (Comp1 m) = foldMap1 (foldMap1 f) m

class Bifoldable t => Bifoldable1 t where
  bifold1 :: Semigroup m => t m m -> m
  bifold1 = bifoldMap1 id id
  {-# INLINE bifold1 #-}

  bifoldMap1 :: Semigroup m => (a -> m) -> (b -> m) -> t a b -> m
  bifoldMap1 f g = maybe (error "bifoldMap1") id . getOption . bifoldMap (Option . Just . f) (Option . Just . g)
  {-# INLINE bifoldMap1 #-}

#if MIN_VERSION_semigroups(0,16,2)
instance Bifoldable1 Arg where
  bifoldMap1 f g (Arg a b) = f a <> g b
#endif

instance Bifoldable1 Either where
  bifoldMap1 f _ (Left a) = f a
  bifoldMap1 _ g (Right b) = g b
  {-# INLINE bifoldMap1 #-}

instance Bifoldable1 (,) where
  bifoldMap1 f g (a, b) = f a <> g b
  {-# INLINE bifoldMap1 #-}

instance Bifoldable1 ((,,) x) where
  bifoldMap1 f g (_,a,b) = f a <> g b
  {-# INLINE bifoldMap1 #-}

instance Bifoldable1 ((,,,) x y) where
  bifoldMap1 f g (_,_,a,b) = f a <> g b
  {-# INLINE bifoldMap1 #-}

instance Bifoldable1 ((,,,,) x y z) where
  bifoldMap1 f g (_,_,_,a,b) = f a <> g b
  {-# INLINE bifoldMap1 #-}

instance Bifoldable1 Const where
  bifoldMap1 f _ (Const a) = f a
  {-# INLINE bifoldMap1 #-}

#ifdef MIN_VERSION_tagged
instance Bifoldable1 Tagged where
  bifoldMap1 _ g (Tagged b) = g b
  {-# INLINE bifoldMap1 #-}
#endif

instance (Bifoldable1 p, Foldable1 f, Foldable1 g) => Bifoldable1 (Biff p f g) where
  bifoldMap1 f g = bifoldMap1 (foldMap1 f) (foldMap1 g) . runBiff
  {-# INLINE bifoldMap1 #-}

instance Foldable1 f => Bifoldable1 (Clown f) where
  bifoldMap1 f _ = foldMap1 f . runClown
  {-# INLINE bifoldMap1 #-}

instance Bifoldable1 p => Bifoldable1 (Flip p) where
  bifoldMap1 f g = bifoldMap1 g f . runFlip
  {-# INLINE bifoldMap1 #-}

instance Bifoldable1 p => Foldable1 (Join p) where
  foldMap1 f (Join a) = bifoldMap1 f f a
  {-# INLINE foldMap1 #-}

instance Foldable1 g => Bifoldable1 (Joker g) where
  bifoldMap1 _ g = foldMap1 g . runJoker
  {-# INLINE bifoldMap1 #-}

instance (Bifoldable1 f, Bifoldable1 g) => Bifoldable1 (Bifunctor.Product f g) where
  bifoldMap1 f g (Bifunctor.Pair x y) = bifoldMap1 f g x <> bifoldMap1 f g y
  {-# INLINE bifoldMap1 #-}

instance (Foldable1 f, Bifoldable1 p) => Bifoldable1 (Tannen f p) where
  bifoldMap1 f g = foldMap1 (bifoldMap1 f g) . runTannen
  {-# INLINE bifoldMap1 #-}

instance Bifoldable1 p => Bifoldable1 (WrappedBifunctor p) where
  bifoldMap1 f g = bifoldMap1 f g . unwrapBifunctor
  {-# INLINE bifoldMap1 #-}

#if MIN_VERSION_base(4,4,0)
instance Foldable1 Complex where
  foldMap1 f (a :+ b) = f a <> f b
  {-# INLINE foldMap1 #-}
#endif

#ifdef MIN_VERSION_containers
instance Foldable1 Tree where
  foldMap1 f (Node a []) = f a
  foldMap1 f (Node a (x:xs)) = f a <> foldMap1 (foldMap1 f) (x :| xs)
#endif

instance Foldable1 Identity where
  foldMap1 f = f . runIdentity

#ifdef MIN_VERSION_tagged
instance Foldable1 (Tagged a) where
  foldMap1 f (Tagged a) = f a
#endif

instance Foldable1 m => Foldable1 (IdentityT m) where
  foldMap1 f = foldMap1 f . runIdentityT

instance Foldable1 f => Foldable1 (Backwards f) where
  foldMap1 f = foldMap1 f . forwards

instance (Foldable1 f, Foldable1 g) => Foldable1 (Compose f g) where
  foldMap1 f = foldMap1 (foldMap1 f) . getCompose

instance Foldable1 f => Foldable1 (Lift f) where
  foldMap1 f (Pure x)  = f x
  foldMap1 f (Other y) = foldMap1 f y

instance (Foldable1 f, Foldable1 g) => Foldable1 (Functor.Product f g) where
  foldMap1 f (Functor.Pair a b) = foldMap1 f a <> foldMap1 f b

instance Foldable1 f => Foldable1 (Reverse f) where
  foldMap1 f = getDual . foldMap1 (Dual . f) . getReverse

instance (Foldable1 f, Foldable1 g) => Foldable1 (Functor.Sum f g) where
  foldMap1 f (Functor.InL x) = foldMap1 f x
  foldMap1 f (Functor.InR y) = foldMap1 f y

instance Foldable1 NonEmpty where
  foldMap1 f (a :| []) = f a
  foldMap1 f (a :| b : bs) = f a <> foldMap1 f (b :| bs)
  toNonEmpty = id

instance Foldable1 ((,) a) where
  foldMap1 f (_, x) = f x

instance Foldable1 g => Foldable1 (Joker g a) where
  foldMap1 g = foldMap1 g . runJoker
  {-# INLINE foldMap1 #-}

max_stars_repo_name: quasicomputational/semigroupoids

max_issues_repo_name: quasicomputational/semigroupoids

max_forks_repo_name: quasicomputational/semigroupoids


text:
import Paraiso
import Data.Complex

main = do
  args <- getArgs
  let arch = if "--cuda" `elem` args then
               CUDA 128 128
             else
               X86
  putStrLn $ compile arch code
    where
      code = do
         parallel 16384 $ do
           r <- allocate
           x <- allocate
           r =$ Rand (0.0::Double) 4.0
           x =$ Rand 0 1
           cuda $ do
             sequential 65536 $ do
               r =$ r * x * (1-x)
           output [r,x]


text:
{-# LANGUAGE OverloadedLists#-}
{-# LANGUAGE RankNTypes#-}
{-# LANGUAGE ScopedTypeVariables #-}
--{-# LANGUAGE Strict#-}
module Lib where

import Control.Monad.ST
import Debug.Trace
import qualified Numeric.AD as AD
import Numeric.AD.Internal.Reverse
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Data
import qualified Data.Vector.Storable as VS
import System.Random
import Debug.Trace
import Text.Printf
import Control.Monad
import Control.Concurrent
import Data.List (permutations)

someFunc :: IO ()
someFunc = putStrLn "someFunc"

data Neuron = Neuron {
                     w  :: Vector Double,
                     f  :: Double -> Double }

data MNeuron = MNeuron {
                       mw :: Matrix Double,
                       mf :: Double -> Double }

instance Show Neuron where
  show = show . w

instance Show MNeuron where
  show = show . mw

applyMNeuron :: MNeuron -> Vector Double -> Vector Double
applyMNeuron (MNeuron wagi funkcja) = VS.map funkcja . (wagi #>) . VS.cons 1

applyNeuron :: Neuron -> Vector Double -> Double
applyNeuron (Neuron w f) = f . (w <.>) . VS.cons 1

acFunc :: Double -> Double
acFunc x | x >= 0 = 1
          | otherwise = 0

notN = Neuron [1,-1.5] acFunc

andN = Neuron [-2, 1.5, 1.5] acFunc

nandN = Neuron [2, -1.5, -1.5] acFunc

orN = Neuron [-1, 1.5, 1.5] acFunc

-- 2. Picture classification

trainSet :: [(Vector Double, Double)]
trainSet = [([ 0, 0, 0, 0, 0
              , 0, 1, 1, 0, 0
              , 0, 0, 1, 0, 0
              , 0, 0, 1, 0, 0
              , 0, 0, 1, 0, 0], 1)

           , ([ 0, 0, 1, 1, 0
              , 0, 0, 0, 1, 0
              , 0, 0, 0, 1, 0
              , 0, 0, 0, 0, 0
              , 0, 0, 0, 0, 0], 1)

           , ([ 0, 0, 0, 0, 0
              , 1, 1, 0, 0, 0
              , 0, 1, 0, 0, 0
              , 0, 1, 0, 0, 0
              , 0, 1, 0, 0, 0], 1)

           , ([ 0, 0, 0, 0, 0
              , 0, 1, 1, 1, 0
              , 0, 1, 0, 1, 0
              , 0, 1, 1, 1, 0
              , 0, 0, 0, 0, 0], 0)

           , ([ 0, 0, 0, 0, 0
              , 0, 0, 0, 0, 0
              , 1, 1, 1, 0, 0
              , 1, 0, 1, 0, 0
              , 1, 1, 1, 0, 0], 0)]

test :: [Vector Double]
test = [[ 0, 1, 1, 0, 0
        , 0, 0, 1, 0, 0
        , 0, 0, 1, 0, 0
        , 0, 0, 0, 0, 0
        , 0, 0, 0, 0, 0
        ]
       ,[ 1, 1, 1, 0, 0
        , 1, 0, 1, 0, 0
        , 1, 0, 1, 0, 0
        , 1, 1, 1, 0, 0
        , 0, 0, 0, 0, 0
        ]
       ,[ 0, 0, 1, 0, 0
        , 0, 0, 1, 0, 0
        , 0, 0, 1, 0, 0
        , 0, 0, 1, 0, 0
        , 0, 0, 0, 0, 0
        ]
       ,[ 0, 0, 0, 0, 0
        , 0, 0, 0, 0, 0
        , 1, 1, 1, 0, 0
        , 1, 0, 1, 0, 0
        , 1, 1, 1, 0, 0
        ]]

train st perc trSet = go st 0 perc (cycle trSet) (length trSet)
  where
    go st iters n@(Neuron w f) ((inp,z):xs) c =
      if c <=0 then
        (iters, n)
      else
        go st (iters+1) (Neuron newW f) xs newC
      where
        y = applyNeuron n inp
        newC = if z == y then c - 1 else 5
        newW = w + VS.map (\x -> x * st * (z - y)) (VS.cons 1 inp)

perceptron :: Neuron
perceptron = Neuron (VS.fromList . replicate 26 $ 1) acFunc

-- 3. Picture recognition

prPrint :: Vector Double -> String
prPrint = morph (\x (acc, count) -> ((if x <= 0 then " " else "■") ++ (if count `mod` 5 == 0 then "\n" else []) ++ acc, count+1)) []
  where
    morph f a = fst . VS.foldr f (a,0)

readPic :: String -> Vector Double
readPic = VS.fromList . foldr (\x acc -> case x of
                                  ' ' -> (-1):acc
                                  '*' -> 1:acc
                                  '\n'-> acc) []
pics :: [String]
pics = [ "     \n *** \n * * \n *** \n     "
       , "     \n **  \n  *  \n  *  \n     "]

sgn :: Double -> Double
sgn x
     | x >= 0 = 1
     | x < 0 = -1

picTaker :: MNeuron
picTaker = MNeuron weights sgn
  where
    [v1, v2] = map readPic pics
    w1       = cmap (/25) (v1 `outer` v1)
    w2       = cmap (/25) (v2 `outer` v2)
    weights = w1 + w2

apMNeu2 (MNeuron weights f) = VS.map f . (weights #>)

pics2 :: [String]
pics2 = [ " *** \n * * \n * * \n *** \n     "
        , "  *  \n  *  \n  *  \n  *  \n  *  "]

-- 4. Gradient Descent
gradStep :: Floating a => a
gradStep = 0.1

eps :: Floating a => a
eps = 0.000001

gradDesc :: (Floating a, Ord a, Show a) => (forall b. (Floating b, Ord b, Show b) => [b] -> b) -> [a] -> ([a], a)
gradDesc f x = if diff < eps then (newX, f newX) else gradDesc f newX
  where gr   = AD.grad f x
        newX = zipWith (-) x . map (gradStep *) $ gr
        diff = abs . foldr (+) 0 . zipWith (-) x $ newX

funkcja1 :: Num a => [a] -> a
funkcja1 [x,y,z] = (2 * x * x) + (2 * y * y) + (z * z) - (2 * x * y) - (2 * y * z) - (2 * x) + 3

dfunkcja1 :: Num a => [a] -> [a]
dfunkcja1 [x,y,z] = [(4*x) - (2*y) - 2, 4*y - 2*x - 2*z, 2*z - 2*y]

funkcja2 :: Num a => [a] -> a
funkcja2 [x,y] = (3 * x * x * x * x) + (4 * x * x* x) - (12 * x * x) + (12 * y * y) - (27 * y)

dfunkcja2 :: Num a => [a] -> [a]
dfunkcja2 [x,y] = [12* x* x* x + 12* x* x - 24* x, 24 * y - 24]

-- 5. Backpropagation

beta :: Double
beta = 1.25

krok :: Double
krok = 1

sig :: Double -> Double
sig = (1/).(1+).exp.negate

dsig :: Double -> Double
dsig x = sig x * (1 - sig x)

trainSetXOR :: [(Vector Double, Double)]
trainSetXOR =
  [([0,0],0)
  ,([1,0],1)
  ,([0,1],1)
  ,([1,1],0)]

fstLayer :: IO (Neuron, Neuron)
fstLayer = do
  randVec1 <- VS.replicateM 3 randomIO
  randVec2 <- VS.replicateM 3 randomIO
  return (Neuron randVec1 sig, Neuron randVec2 sig)

outLayer :: IO Neuron
outLayer = do
  randVec <- VS.replicateM 3 randomIO
  return (Neuron randVec sig)

beginning :: ((Neuron, Neuron), Neuron)
beginning = ((Neuron [2,0,1] sig, Neuron [2,0,1] sig), Neuron [2,0,1] sig)

applyXor :: ((Neuron,Neuron), Neuron) -> Vector Double -> Double
applyXor ((h1, h2), o) u = applyNeuron o out
  where
    x1 = applyNeuron h1 u
    x2 = applyNeuron h2 u
    out= [x1, x2]

infixr 3 #*
(#*) :: Double -> Vector Double -> Vector Double
x #* v = VS.map (*x) v

zipWith4 :: (a -> b -> c -> d -> e) -> [a] -> [b] -> [c] -> [d] -> [e]
zipWith4 f x y z = zipWith ($) (zipWith3 f x y z)

trainXOR perc trSet = go 0 perc
  where
    go iters n@((n1@(Neuron w1 f1), n2@(Neuron w2 f2)), n3@(Neuron s f3)) =
      if eps > (maximum . VS.toList .  VS.map abs $ difference) then
        (iters, newN)
      else
--        (if iters `mod` 10000 == 0 then trace (show iters ++ " iterations. Gradients:\n" ++ unlines [show gradS,show gradW1,show gradW2] ++ "\nDifference:\n" ++ show difference ++ "\n") else id ) $
        go (iters+1) newN
      where
        difference = (w1 VS.++ w2 VS.++s) - (newW1 VS.++ newW2 VS.++ newS)
        newN =  ((Neuron newW1 f1, Neuron newW2 f2), Neuron newS f3)
        u   = map (VS.cons 1 . fst) trSet
        x1s = map (applyNeuron n1 . fst) trSet
        x2s = map (applyNeuron n2 . fst) trSet
        x :: [Vector Double]
        x   = zipWith (\x1 x2 -> [1,x1,x2]) x1s x2s
        y   = map (applyXor n . fst) trSet
        z   = map snd trSet
        grs x y z = dsig (x<.>s) * (y - z) #* x
        gradsS= zipWith3 grs x y z
        gradS = foldr (+) [0,0,0] gradsS
        newS = s - (VS.map (*gradStep) gradS)
        [_,s1,s2] = s
        grw w sx x y z u= (y-z) * sx * dsig (x<.>s) * dsig (w<.>u) #* u
        gradsW w s = zipWith4 (grw w s) x y z u
        gradW w s = foldr (+) [0,0,0] $ gradsW w s
        gradW1 = gradW w1 s1
        gradW2 = gradW w2 s2
        newW1 = w1 - (VS.map (*gradStep) $ gradW w1 s1)
        newW2 = w2 - (VS.map (*gradStep) $ gradW w2 s2)

prettyPrintVals :: [Double] -> String
prettyPrintVals = unwords . map (printf "%.4f")

-- Zad. 6 Encoder - Decoder

autoEncTrain = map readPic [ "     \n **  \n  *  \n  *  \n  *  "
                            , "  ** \n   * \n   * \n     \n     "
                            , "     \n**   \n *   \n *   \n *   "]

type EncDec = (MNeuron, MNeuron)

makeDec (_, MNeuron w _) = MNeuron w acFunc

zeroM inp out = (out><inp) $ replicate (out*inp) 0

--startEncDec :: IO EncDec
startEncDec = do --(MNeuron (zeroM 26 16) sig, MNeuron (zeroM 17 25) sig) --do
  m1 <- map (\x -> x - 0.5) <$> replicateM (25*26) (randomIO :: IO Double)
  m2 <- map (\x -> x - 0.5) <$> replicateM (25*17) (randomIO :: IO Double)
  return (MNeuron ((16><26) m1) sig, MNeuron ((25><17) m2) sig)

decode (_,a) = applyMNeuron a

encode (a,_) = applyMNeuron a

ende a = applyMNeuron (makeDec a) . encode a

alfa = 0.8

trainED trSet = trainHelp 0
  where trainHelp old ed@(MNeuron w f1, MNeuron w' f2) =
          if abs (old - rmse) < 0.000001 then ed
          else trace ("RMSE = "++ show rmse) $ trainHelp rmse (MNeuron nW f1, MNeuron nW' f2)
          where nW            = w - cmap (*alfa) gradW
                nW'           = w'- cmap (*alfa) gradW'
                diffDecode :: Vector Double -> Vector Double
                diffDecode    = cmap dsig . (#>) w' . VS.cons 1
                diffEncode    = cmap dsig . (#>) w  . VS.cons 1
                gradDec x y x'= let dx = (x' - x) * diffDecode y
                                in dx `outer` VS.cons 1 y
                gradEnc x y x'= let dx = (x' - x) * diffDecode y
                                    dy = VS.tail (tr w' #> dx) * diffEncode x
                                in dy `outer` VS.cons 1 x
                triples       = map (\x -> (x, encode ed x,decode ed . encode ed $ x )) trSet
                gradsW'       = map (\(x,y,x') -> gradDec x y x') triples
                gradsW        = map (\(x,y,x') -> gradEnc x y x') triples
                rmse          = sqrt . VS.foldr (+) 0 . foldr (\(x,_,x') acc -> (x-x')^2 +acc ) (fromList . replicate 25 $ 0) $ triples
                gradW'        = foldr (+) (zeroM 17 25) gradsW'
                gradW         = foldr (+) (zeroM 26 16) gradsW

-- Zad 7. Sieci Hopfielda

readPic2 :: String -> Vector Double
readPic2 = VS.fromList . foldr (\x acc -> case x of
                                  ' ' -> 0:acc
                                  '*' -> 1:acc
                                  '\n'-> acc) []

jedynka :: Vector Double
jedynka = readPic2 "     \n **  \n  *  \n  *  \n  *  "

wagi :: Matrix Double
wagi = cmap (*2) $ (xs `outer` xs) * (cmap ((+1).negate) $ ident 25)
  where xs = cmap (flip (-) 0.5) jedynka

bias :: Vector Double
bias = cmap (pred.(/2)) (wagi #> fromList (replicate 25 1))

iteration :: Vector Double -> Vector Double
iteration = (flip (-) bias).(wagi #>) >>= VS.zipWith (\u x -> if u > 0 then 1 else if u < 0 then 0 else x)

inp :: Vector Double
inp = readPic2 " *   \n   * \n* *  \n    *\n *   "

process = process . iteration . (prPrint >>= trace)

randInp :: IO (Vector Double)
randInp = VS.fromList . map (\x -> if x > 0.5 then 1 else 0) <$> replicateM 25 (randomIO :: IO Double)

zeroV :: Vector Double
zeroV = readPic2 " *** \n * * \n * * \n * * \n *** "

wagi2 :: Matrix Double
wagi2 = cmap (*2) (((cx `outer` cx) + (dx `outer` dx)) * (cmap ((+1).negate) $ ident 25))
  where dx = cmap (flip (-) 0.5) zeroV
        cx = cmap (flip (-) 0.5) jedynka

bias2 :: Vector Double
bias2 = cmap (pred.pred.(/2)) (wagi2 #> fromList (replicate 25 1))

iterate2 = (flip (-) bias2).(wagi2 #>) >>= VS.zipWith (\u x -> if u > 0 then 1 else if u < 0 then 0 else x)

inp2 :: Vector Double
inp2 = readPic2 "  *  \n * * \n     \n* *  \n    *"

-- Boltzmann

prob :: Double -> Double -> Double
prob t x = (1/).(1+).exp.negate $ (x/t)

boltzStep :: Double -> Vector Double -> IO (Vector Double)
boltzStep temp = VS.mapM (\x -> randomIO >>= \y -> if y <= prob temp x then return 1 else return 0).(flip (-) bias).(wagi#>)

--replicateIO :: (a -> IO a) -> a -> IO [a]
--replicateIO f a = [ a : rest | rest <- replicateIO f nA , nA <- f a ]

annealing :: Double -> Vector Double -> IO ()
annealing = aSt 0

aSt :: Double -> Double -> Vector Double -> IO ()
aSt x inT input = do
  nv <- boltzStep (tann inT x) input
  threadDelay 300000
  putStrLn $ "Step number: " ++ show x
  putStrLn $ "Temperature: " ++ show (tann inT x)
  putStrLn . prPrint $ input
  aSt (x+1) inT nv

bSt :: Double -> Vector Double -> IO ()
bSt temp input = do
  nv <- boltzStep temp input
  threadDelay 300000
  putStrLn . prPrint $ input
  bSt temp nv

tann :: Double -> Double -> Double
tann t = (t/).(1+).log.(1+)

-- TSP

randomPermutation :: Int -> IO [Int]
randomPermutation n = do
  let ps= permutations [0..n-1]
      x = length ps
  r <- randomRIO (0,x-1) :: IO Int
  return $ ps !! r

odleglosci1 = cmap abs (row [1..10] - col [1..10]) - (col (8:replicate 9 0) * row (replicate 9 0 ++ [1])) - (col (replicate 9 0++[8]) * row (1 : replicate 9 0))

max_stars_repo_name: Antystenes/Sieci

max_issues_repo_name: Antystenes/Sieci

max_forks_repo_name: Antystenes/Sieci


text:
module EvalSpec (spec) where

import           Control.Monad.Error
import           Data.Complex
import           Data.Either
import           Data.IORef
import           Data.Ratio
import           Env
import           Eval
import           Read
import           System.IO.Unsafe
import           Test.Hspec
import           Types

-- same as Main.hs
primitiveBindings :: IO Env
primitiveBindings = nullEnv >>= flip bindVars (map makePrimitiveFunc primitives)
  where makePrimitiveFunc (var, func) = (var, PrimitiveFunc func)

runIOThrowsToLispVal :: IOThrowsError LispVal -> IO LispVal
runIOThrowsToLispVal action = fmap extractValue (runErrorT action)

exprToVal :: String -> IOThrowsError LispVal
exprToVal expr = liftThrows $ readExpr expr

evalExpr :: Env -> String -> IO LispVal
evalExpr env expr = runIOThrowsToLispVal $ exprToVal expr >>= eval env

re expr = unsafePerformIO $ do
  env <- primitiveBindings
  evalExpr env expr

le expr = unsafePerformIO $ do
  env <- primitiveBindings
  evalString env expr

spec :: Spec
spec = do
  describe "eval" $ do
    it "evaluates functions" $ do
      re "(+ 1 2 3)" `shouldBe` Number 6
      re "(- 3 1)" `shouldBe` Number 2
      re "(- 10 1 2)" `shouldBe` Number 7
      re "(* 1 2 3 4 5)" `shouldBe` Number 120
      re "(/ 23 7)" `shouldBe` Number 3
      re "(/ 51 7 3)" `shouldBe` Number 2
      re "(mod 7 3)" `shouldBe` Number 1
      re "(quotient 7 3)" `shouldBe` Number 2
      re "(remainder 7 3)" `shouldBe` Number 1
      re "(symbol? 'a)" `shouldBe` Bool True
      re "(symbol? 1)" `shouldBe` Bool False
      re "(string? \"foo\")" `shouldBe` Bool True
      re "(string? 1)" `shouldBe` Bool False
      re "(number? 1)" `shouldBe` Bool True
      re "(number? 'a)" `shouldBe` Bool False
      re "(boolean? #t)" `shouldBe` Bool True
      re "(boolean? 'f)" `shouldBe` Bool False

      re "(= 1 1)" `shouldBe` Bool True
      re "(= 1 2)" `shouldBe` Bool False
      re "(> 1 0)" `shouldBe` Bool True
      re "(> 1 2)" `shouldBe` Bool False
      re "(< 1 0)" `shouldBe` Bool False
      re "(< 1 2)" `shouldBe` Bool True
      re "(/= 1 1)" `shouldBe` Bool False
      re "(/= 1 2)" `shouldBe` Bool True
      re "(>= 1 2)" `shouldBe` Bool False
      re "(>= 1 1)" `shouldBe` Bool True
      re "(<= 1 0)" `shouldBe` Bool False
      re "(<= 1 1)" `shouldBe` Bool True
      re "(&& #t #t)" `shouldBe` Bool True
      re "(&& #t #f)" `shouldBe` Bool False
      re "(&& #f #f)" `shouldBe` Bool False
      re "(|| #t #t)" `shouldBe` Bool True
      re "(|| #t #f)" `shouldBe` Bool True
      re "(|| #f #f)" `shouldBe` Bool False

    it "string?" $ do
      re "(string=? \"a\" \"a\")" `shouldBe` Bool True
      re "(string=? \"a\" \"b\")" `shouldBe` Bool False
      re "(string<? \"a\" \"a\")" `shouldBe` Bool False
      re "(string<? \"a\" \"b\")" `shouldBe` Bool True
      re "(string>? \"a\" \"a\")" `shouldBe` Bool False
      re "(string>? \"c\" \"b\")" `shouldBe` Bool True
      re "(string<=? \"b\" \"a\")" `shouldBe` Bool False
      re "(string<=? \"b\" \"b\")" `shouldBe` Bool True
      re "(string>=? \"a\" \"b\")" `shouldBe` Bool False
      re "(string>=? \"b\" \"b\")" `shouldBe` Bool True

    it "if" $ do
      re "(if (> 2 3) \"foo\" \"bar\")" `shouldBe` String "bar"
      re "(if (= 3 3) (+ 2 3 (- 5 1)) \"unequal\")" `shouldBe` Number 9

    it "car, cdr" $ do
      re "(cdr '(a simple test))" `shouldBe` List [Atom "simple", Atom "test"]
      re "(car (cdr '(a simple test)))" `shouldBe` Atom "simple"
      re "(car '((this is) a test))" `shouldBe` List [Atom "this", Atom "is"]
      re "(cons '(this is) 'test)" `shouldBe` DottedList [List [Atom "this", Atom "is"]] (Atom "test")
      re "(cons '(this is) '())" `shouldBe` List [List [Atom "this", Atom "is"]]

    it "eqv?" $ do
      re "(eqv? 1 3)" `shouldBe` Bool False
      re "(eqv? 3 3)" `shouldBe` Bool True
      re "(eqv? 'atom 'atom)" `shouldBe` Bool True

    it "equal?" $ do
      re "(equal? 123 123)" `shouldBe` Bool True
      re "(equal? 123 \"123\")" `shouldBe` Bool True
      re "(equal? 123 \"foo\")" `shouldBe` Bool False

    it "throws" $ do
      le "(symbol? 1 2)" `shouldStartWith` "Expected 1 args;"
      le "(string? 1 2)" `shouldStartWith` "Expected 1 args;"
      le "(number? 1 2)" `shouldStartWith` "Expected 1 args;"
      le "(boolean? 1 2)" `shouldStartWith` "Expected 1 args;"
      le "(+ 2 \"two\")" `shouldStartWith` "Invalid type:"
      le "(+ 2)" `shouldStartWith` "Expected 2 args;"
      le "(what? 2)" `shouldBe` "Getting an unbound variable: : what?"
      le "(if 1 2 3)" `shouldStartWith` "Invalid type:"

  describe "env" $ do
    it "defines functions" $ do
      env <- primitiveBindings

      evalExpr env "(define (f x y) (+ x y))"
      unsafePerformIO (evalExpr env "(f 1 2)") `shouldBe` Number 3

      evalExpr env "(define (factorial x) (if (= x 1) 1 (* x (factorial (- x 1)))))"
      unsafePerformIO (evalExpr env "(factorial 10)") `shouldBe` Number 3628800

      evalExpr env "(define (counter inc) (lambda (x) (set! inc (+ x inc)) inc))"
      evalExpr env "(define my-count (counter 5))"
      unsafePerformIO (evalExpr env "(my-count 3)") `shouldBe` Number 8
      unsafePerformIO (evalExpr env "(my-count 6)") `shouldBe` Number 14
      unsafePerformIO (evalExpr env "(my-count 5)") `shouldBe` Number 19

    it "fails if function is not defined" $
      le "(f 1 2)" `shouldBe` "Getting an unbound variable: : f"

    it "defines functions" $ do
      env <- primitiveBindings
      evalExpr env "(define (f x y) (+ x y))"
      unsafePerformIO (evalExpr env "(f 1 2)") `shouldBe` Number 3

    it "fails if function is not defined" $
      le "(f 1 2)" `shouldBe` "Getting an unbound variable: : f"


text:
{-# LANGUAGE PatternSynonyms, ViewPatterns #-}
module ParsersTest (parsersTests) where

import Prelude hiding (getChar)

import Test.Tasty
import Test.Tasty.HUnit
import Test.Tasty.QuickCheck

import Data.Complex (Complex(..))
import Data.Maybe
import Data.Either
import qualified Data.Vector as V

import Control.Arrow (first)

import Parsers
import Val
import LispValTest ()
import EvaluationTest ((?=))

pattern IList xs <- (unrollList -> Just xs)
  where IList xs = makeImmutableList xs

pattern IDottedList xs x <- (unrollDList -> Just (xs, x))
  where IDottedList xs x = makeImproperImmutableList xs x

parsersTests :: TestTree
parsersTests = testGroup "Parser tests" [unitTests, propTests]

unitTests :: TestTree
unitTests = testGroup "Unit tests"
    [ stringParserTests
    , atomParserTests
    , numberParserTests
    , charParserTests
    , listParserTests
    , dotListParserTests
    , vectorParserTests
    , endToEndTests
    ]

propTests :: TestTree
propTests = testGroup "Property tests" [qcTests]

qcTests :: TestTree
qcTests = testGroup "(QuickCheck)"
    [ prop_QuoteParser
    , prop_QuasiquoteParser
    , prop_UnquoteParser
    , prop_UnquoteSplicingParser
    ]

-- END TO END TESTS
endToEndTests = testGroup "End to End" $ map testParseExpr
    [ mkE2Etest
        { testName = "#d Number"
        , input = "12"
        , expectedContents = Just $ Number 12
        }
    , mkE2Etest
        { testName = "#x, negative number"
        , input = "#x-ef"
        , expectedContents = Just $ Number (-239)
        }
    , mkE2Etest
        { testName = "#o, number with +"
        , input = "#o+55"
        , expectedContents = Just $ Number 45
        }
    , mkE2Etest
        { testName = "Char"
        , input = "#\\f"
        , expectedContents = Just $ Char 'f'
        }
    , mkE2Etest
        { testName = "Regular Symbol"
        , input = "symbol"
        , expectedContents = Just $ Symbol "symbol"
        }
    , mkE2Etest
        { testName = "#t"
        , input = "#t"
        , expectedContents = Just $ Bool True
        }
    , mkE2Etest
        { testName = "#f"
        , input = "#f"
        , expectedContents = Just $ Bool False
        }
    , mkE2Etest
        { testName = "Simple string"
        , input = "\"string\""
        , expectedContents = Just $ IString "string"
        }
    , mkE2Etest
        { testName = "Quoted atom"
        , input = "'$name"
        , expectedContents = Just $ IList [Symbol "quote", Symbol "$name"]
        }
    , mkE2Etest
        { testName = "Quoted char"
        , input = "'#\\*"
        , expectedContents = Just $ IList [Symbol "quote", Char '*']
        }
    , mkE2Etest
        { testName = "Immediately nested quote forms"
        , input = "`,4"
        , expectedContents = Just $ IList [Symbol "quasiquote", IList
                                          [Symbol "unquote", Number 4]]
        }
    , mkE2Etest
        { testName = "Quasiquoted mess"
        , input = "`(,(foo) ,@(bar) x 1)"
        , expectedContents = Just $
            IList [Symbol "quasiquote", IList
                  [ IList [Symbol "unquote", IList [Symbol "foo"]]
                  , IList [Symbol "unquote-splicing", IList [Symbol "bar"]]
                  , Symbol "x"
                  , Number 1
                  ]]
        }
    , mkE2Etest
        { testName = "Nested List"
        , input = "(f \"x\" #\\5 (+ 3))"
        , expectedContents = Just $
            IList [Symbol "f", IString "x", Char '5', IList
                  [Symbol "+", Number 3]
            ]
        }
    , mkE2Etest
        { testName = "List with []"
        , input = "[+ 1 2]"
        , expectedContents = Just $
            IList [Symbol "+", Number 1, Number 2]
        }
    , mkE2Etest
        { testName = "List with {}"
        , input = "{test-func #\\a 0}"
        , expectedContents = Just $
            IList [Symbol "test-func", Char 'a', Number 0]
        }
    , mkE2Etest
        { testName = "Regular dotted List"
        , input = "(f 5 #\\$ . (1 2 \"test\"))"
        , expectedContents = Just $
            IList [Symbol "f", Number 5, Char '$'
                  , Number 1, Number 2, IString "test"]
        }
    , mkE2Etest
        { testName = "Dotted with nil"
        , input = "(+ 0 . ())"
        , expectedContents = Just $ IList [Symbol "+", Number 0]
        }
    , mkE2Etest
        { testName = "Irregular dotted List"
        , input = "(f 0 . m)"
        , expectedContents = Just $ IDottedList [Symbol "f", Number 0] (Symbol "m")
        }
    , mkE2Etest
        { testName = "Degenerate dotted List"
        , input = "(. xs)"
        , expectedContents = Just $ Symbol "xs"
        }
    , mkE2Etest
        { testName = "Empty vector"
        , input = "#()"
        , expectedContents = Just $ IVector V.empty
        }
    , mkE2Etest
        { testName = "regular vector"
        , input = "#(\"test\" 0 atom)"
        , expectedContents = Just $ IVector $ V.fromList
            [IString "test", Number 0, Symbol "atom"]
        }
    , mkE2Etest
        { testName = "Irregular spacing"
        , input = "(f 5 #\\$(1 2 #x45) #b10(test))"
        , expectedContents = Just $
            IList
                [Symbol "f", Number 5, Char '$', IList
                    [Number 1, Number 2, Number 69]
                , Number 2, IList
                    [Symbol "test"]
                ]
        }
    , mkE2Etest
        { testName = "Invalid #d number and invalid atom"
        , input = "(+ 1 5b)"
        }
    , mkE2Etest
        { testName = "Invalid symbol"
        , input = "symb`ol"
        }
    , mkE2Etest
        { testName = "Too few parens"
        , input = "(test-func 1 (cons (car xs) (cdr xs))"
        }
    , mkE2Etest
        { testName = "Too many parens"
        , input = "(test-func \"t\" car xs))"
        }
    , mkE2Etest
        { testName = "Extra input after list expr"
        , input = "(test-func + 0) randomwords"
        }
    , mkE2Etest
        { testName = "Extra input after normal expr"
        , input = "#\\c (test-func 0 1)"
        }
    ]

mkE2Etest :: TestBuilder Val
mkE2Etest = TB End2End "" "" Nothing

testParseExpr :: TestBuilder Val -> TestTree
testParseExpr testBuilder =
    testCaseSteps (testName testBuilder) $ \step -> do
        step "Parsing"
        let pInput = input testBuilder
        parse <- run pInput

        if isJust $ expectedContents testBuilder
        then do
            step "Verifying success"
            isRight parse @? "Parse failed on: " ++ pInput

            let val = fromRight undefined parse

            step "Verifying parse"
            -- TODO: we should also verify that the result is
            -- appropriately mutable/immutable.
            Right val ?= Right (fromJust $ expectedContents testBuilder)
        else do
            step "Verifying failure"
            isLeft parse @? "Parse succeeded on: " ++ pInput

-- UNIT TESTS
stringParserTests = testGroup "Parsing Strings" $ map testStringParser
    [ mkStringTest
        { testName = "Simple string"
        , input = "\"simple test string\""
        , expectedContents = Just "simple test string"
        }
    , mkStringTest
        { testName = "Missing start quote"
        , input = "simple test string\""
        }
    , mkStringTest
        { testName = "Missing end quote"
        , input = "\"simple test string"
        }
    , mkStringTest
        { testName = "Unquoted"
        , input = "simple test string"
        }
    , mkStringTest
        { testName = "String with \\"
        , input = "\"test\\\\string\""
        , expectedContents = Just "test\\string"
        }
    , mkStringTest
        { testName = "String with \""
        , input =  "\"test \\\"string in a string\\\" string\""
        , expectedContents = Just "test \"string in a string\" string"
        }
    , mkStringTest
        { testName = "String with newline"
        , input = "\"test string\\nwith two lines\""
        , expectedContents = Just "test string\nwith two lines"
        }
    , mkStringTest
        { testName = "String with carriage return"
        , input = "\"test string \\r with carriage ret\""
        , expectedContents = Just "test string \r with carriage ret"
        }
    , mkStringTest
        { testName = "String with tab"
        , input = "\"test string\\t with tab\""
        , expectedContents = Just "test string\t with tab"
        }
    ]

atomParserTests :: TestTree
atomParserTests = testGroup "Parsing Symbols" $ map testSymbolParser
    [ let symbol = "atom" in mkSymbolTest
        { testName = "Simple atom"
        , input = symbol
        , expectedContents = Just $ pack symbol
        }
    , let symbol = "#name" in mkSymbolTest
        { testName = "Symbol with #"
        , input = symbol
        }
    , let symbol = "!@#$%&*-_=+|:\\/?<>~" in mkSymbolTest
        { testName = "Symbol with all symbols"
        , input = symbol
        , expectedContents = Just $ pack symbol
        }
    , let symbol = "t1" in mkSymbolTest
        { testName = "Symbol with digit"
        , input = symbol
        , expectedContents = Just $ pack symbol
        }
    , mkSymbolTest
        { testName = "Symbol starts with digit"
        , input = "5a"
        }
    , mkSymbolTest
        { testName = "Symbol with brackets"
        , input = "na[]me"
        }
    , mkSymbolTest
        { testName = "Symbol with braces"
        , input = "na{}me"
        }
    , mkSymbolTest
        { testName = "Symbol with invalid symbols"
        , input = "AlmostValidSymbol`name"
        }
    , let symbol = "Complex#Valid15Symbol\\Name8?" in mkSymbolTest
        { testName = "Complex valid atom"
        , input = symbol
        , expectedContents = Just $ pack symbol
        }
    , mkSymbolTest
        { testName = "vertical bar symbol"
        , input = "|H\\x65;llo|"
        , expectedContents = Just "Hello"
        }
    ]

numberParserTests :: TestTree
numberParserTests = testGroup "Parsing Numbers" $ map testNumberParser
    [ mkNumTest
        { testName = "25"
        , input = "25"
        , expectedContents = Just 25
        }
   , mkNumTest
        { testName = "-2"
        , input = "-2"
        , expectedContents = Just (-2)
        }
    , mkNumTest
        { testName = "+15"
        , input = "+15"
        , expectedContents = Just 15
        }
    , mkNumTest
        { testName = "#d80915"
        , input = "#d80915"
        , expectedContents = Just 80915
        }
    , mkNumTest
        { testName = "#d-100"
        , input = "#d-100"
        , expectedContents = Just (-100)
        }
    , mkNumTest
        { testName = "#d+06"
        , input = "#d+06"
        , expectedContents = Just 6
        }
    , mkNumTest
        { testName = "#b1010"
        , input = "#b1010"
        , expectedContents = Just 10
        }
    , mkNumTest
        { testName = "#b-11101"
        , input = "#b-11101"
        , expectedContents = Just (-29)
        }
    , mkNumTest
        { testName = "#b+1011"
        , input = "#b+1011"
        , expectedContents = Just 11
        }
    , mkNumTest
        { testName = "#o517"
        , input = "#o517"
        , expectedContents = Just 335
        }
    , mkNumTest
        { testName = "#o-24"
        , input = "#o-24"
        , expectedContents = Just (-20)
        }
    , mkNumTest
        { testName = "#o+36"
        , input = "#o+36"
        , expectedContents = Just 30
        }
    , mkNumTest
        { testName = "#xf4"
        , input = "#xf4"
        , expectedContents = Just 244
        }
    , mkNumTest
        { testName = "#x-cb8"
        , input = "#x-cb8"
        , expectedContents = Just (-3256)
        }
    , mkNumTest
        { testName = "#x+aa1"
        , input = "#x+aa1"
        , expectedContents = Just 2721
        }
    , mkNumTest
        { testName = "0.5"
        , input = "0.5"
        , expectedContents = Just (Real (Flonum 0.5))
        }
    , mkNumTest
        { testName = "#e0.5"
        , input = "#e0.5"
        , expectedContents = Just (Real (Ratnum 0.5))
        }
    , mkNumTest
        { testName = "infinity"
        , input = "+inf.0"
        , expectedContents = Just (Real (Flonum (1/0)))
        }
    , mkNumTest
        { testName = "rational"
        , input = "3/5"
        , expectedContents = Just (Real (Ratnum (3/5)))
        }
    , mkNumTest
        { testName = "rectangular complex"
        , input = "1+1i"
        , expectedContents = Just (Complex (1 :+ 1))
        }
    , mkNumTest
        { testName = "rectangular with scientific notation"
        , input = "2e+5+3i"
        , expectedContents = Just (Complex (2e+5 :+ 3))
        }
    , mkNumTest
        -- this is _almost_ a peculiar identifier, but the dot subsequent
        -- cannot be a number.
        { testName = "werid rectangular 0"
        , input = "-.0i"
        , expectedContents = Just (Complex (0 :+ (-0.0)))
        }
    , mkNumTest
        { testName = "15a3"
        , input = "15a3"
        }
    , mkNumTest
        { testName = "#b1201"
        , input = "#b1201"
        }
    , mkNumTest
        { testName = "#bnx"
        , input = "#bnx"
        }
    , mkNumTest
        { testName = "#o358"
        , input = "#o358"
        }
    , mkNumTest
        { testName = "#oa"
        , input = "#oa"
        }
    , mkNumTest
        { testName = "#x54mk"
        , input = "#x54mk"
        }
    , mkNumTest
        { testName = "Empty string"
        , input = ""
        }
    ]

charParserTests :: TestTree
charParserTests = testGroup "Parsing Chars" $ map testCharParser
    [ mkCharTest
        { testName = "Basic single character"
        , input = "#\\k"
        , expectedContents = Just 'k'
        }
    , mkCharTest
        { testName = "Basic single character"
        , input = "#\\$"
        , expectedContents = Just '$'
        }
    , mkCharTest
        { testName = "Basic single character"
        , input = "#\\`"
        , expectedContents = Just '`'
        }
    , mkCharTest
        { testName = "Too many characters"
        , input = "#\\ab"
        }
    , mkCharTest
        { testName = "No \\"
        , input = "#8"
        }
    , mkCharTest
        { testName = "Space"
        , input = "#\\space"
        , expectedContents = Just ' '
        }
    , mkCharTest
        { testName = "Newline"
        , input = "#\\newline"
        , expectedContents = Just '\n'
        }
    , mkCharTest
        { testName = "Carriage return"
        , input = "#\\return"
        , expectedContents = Just '\r'
        }
    , mkCharTest
        { testName = "Tab"
        , input = "#\\tab"
        , expectedContents = Just '\t'
        }
    , mkCharTest
        { testName = "escaped unicode"
        , input = "#\\x3bb"
        , expectedContents = Just 'λ'
        }
    ]

listParserTests :: TestTree
listParserTests = testGroup "Parsing Lists"
    [ testCase "\"(5 a #\\a)\"" $ do
        p <- run "(5 a #\\a)"
        isRight p @? "Parse failed."
        let val = unrollList $ fromRight undefined p
        isJust val @? "Parse did not return a list."
        let correct = case fromJust val of
                [Number 5, Symbol "a", Char 'a'] -> True
                _ -> False
        correct @? "Contents of the list are wrong: " ++ show (fromJust val)
    ]

dotListParserTests :: TestTree
dotListParserTests = testGroup "Parsing Dotted Lists"
    [ testCase "\"(5 a . #\\a)\"" $ do
        p <- run "(5 a . #\\a)"
        isRight p @? "Parse failed."
        let Right val = p
            mpair = unrollDList val
        isJust mpair @? "Parse did not return a dotted list."
        let correct = case fromJust mpair of
              ([Number 5, Symbol "a"], Char 'a') -> True
              _ -> False
        correct @? "Contents of the dotted list are wrong: " ++ show val
    ]

vectorParserTests :: TestTree
vectorParserTests = testGroup "Parsing vectors"
    [ testCase "\"#(1 #\\x ())\"" $ do
        p <- run "#(1 #\\x ())"
        val <- verify p
        fromJust val @?= V.fromList [Number 1, Char 'x', Nil]
    , testCase "\"#()\"" $ do
        p <- run "#()"
        val <- verify p
        fromJust val @?= V.empty
    ]
  where verify p = do
            isRight p @? "Parse failed."
            let val = case fromRight undefined p of
                    IVector v -> Just v
                    _ -> Nothing
            isJust val @? "Parse did not return an IVector."
            return val

-- PROPERTY TESTS
prop_QuoteParser :: TestTree
prop_QuoteParser = testProperty "tick always results in quote" $
    withMaxSuccess 50 $ \input -> ioProperty $ do
        result <- run ('\'' : show (input :: Val))
        return $ isRight result ==> case result of
            Right (IPair (Symbol "quote") _) -> True
            _ -> False

prop_QuasiquoteParser :: TestTree
prop_QuasiquoteParser = testProperty "backtick always results in quasiquote" $
    withMaxSuccess 50 $ \input -> ioProperty $ do
        result <- run ('`' : show (input :: Val))
        return $ isRight result ==> case result of
            Right (IPair (Symbol "quasiquote") _)  -> True
            _ -> False

prop_UnquoteParser :: TestTree
prop_UnquoteParser = testProperty "comma always results in unquote" $
    withMaxSuccess 50 $ \input -> ioProperty $ do
        result <- run (',' : clean (show (input :: Val)))
        return $ isRight result ==> case result of
            Right (IPair (Symbol "unquote") _)  -> True
            _ -> False
  -- if we don't do this, a LispString/Symbol that starts with '@' will turn it into
  -- an unquote splicing, which happened in a test at least once.
  where clean = dropWhile (== '@')

prop_UnquoteSplicingParser :: TestTree
prop_UnquoteSplicingParser = testProperty "comma@ always results in unquote-splicing" $
    withMaxSuccess 50 $ \input -> ioProperty $ do
        result <- run (",@" ++ show (input :: Val))
        return $ isRight result ==> case result of
            Right (IPair (Symbol "unquote-splicing") _)  -> True
            _ -> False

run :: String -> IO (Either LispErr Val)
run = labeledReadExpr ""

parseSucceeds :: String -> IO Bool
parseSucceeds s = isRight <$> run s

testStringParser :: TestBuilder Text -> TestTree
testStringParser tb = testParser tb getString

testSymbolParser :: TestBuilder Symbol -> TestTree
testSymbolParser tb = testParser tb getSymbol

testNumberParser :: TestBuilder Number -> TestTree
testNumberParser tb = testParser tb getNumber

testCharParser :: TestBuilder Char -> TestTree
testCharParser tb = testParser tb getChar

data TestType = TString | TSymbol | TNumber | TChar | End2End
instance Show TestType where
    show TString = "String"
    show TSymbol = "Symbol"
    show TNumber = "Number"
    show TChar   = "Char"
    show End2End = "End to End"

data TestBuilder a = TB { testType :: TestType
                        , testName :: String
                        , input :: String
                        , expectedContents :: Maybe a
                        }

mkStringTest, mkSymbolTest, mkNumTest, mkCharTest :: TestBuilder a
mkStringTest = TB TString "String Parser Test" "" Nothing
mkSymbolTest = TB TSymbol "Symbol Parser Test" "" Nothing
mkNumTest = TB TNumber "Number Parser Test" "" Nothing
mkCharTest = TB TChar "Char Parser Test" "" Nothing

testParser :: (Eq a, Show a)
           => TestBuilder a
           -> (Val -> Maybe a)
           -> TestTree
testParser testBuilder decons =
    testCaseSteps (testName testBuilder) $ \step -> do
        step "Parsing"
        parse <- run $ input testBuilder

        if isJust $ expectedContents testBuilder
        then do
            step "Verifying success"
            isRight parse @? "Parse failed on: " ++ input testBuilder

            step "Verifying type"
            let val = (decons . fromRight undefined) parse

            isJust val @? "Parse did not produce a "
                          ++ show (testType testBuilder) ++ "."

            step "Verifying contents"
            fromJust val @?= fromJust (expectedContents testBuilder)
        else do
            step "Verifying failure"
            isLeft parse @? "Parse succeeded on: " ++ input testBuilder

getString :: Val -> Maybe Text
getString (IString s) = Just s
getString _          = Nothing

getSymbol :: Val -> Maybe Symbol
getSymbol (Symbol s) = Just s
getSymbol _        = Nothing

getNumber :: Val -> Maybe Number
getNumber (Number n) = Just n
getNumber _ = Nothing

getChar :: Val -> Maybe Char
getChar (Char c) = Just c
getChar _        = Nothing

unrollList :: Val -> Maybe [Val]
unrollList (IPair c d) =
    (c:) <$> unrollList d
unrollList Nil = Just []
unrollList _ = Nothing

unrollDList :: Val -> Maybe ([Val], Val)
unrollDList (IPair c d) =
    first (c:) <$> unrollDList d
unrollDList Nil = Nothing -- not dotted
unrollDList obj = Just ([], obj)

max_stars_repo_name: JKTKops/Thcheme

max_issues_repo_name: JKTKops/Thcheme

max_forks_repo_name: JKTKops/Thcheme


text:
{-# LANGUAGE Strict #-}

module HEP.Util.Polynomial (quarticEqSol, quadEqSol) where

import Numeric.GSL.Polynomials (polySolve)

import Data.Complex            (imagPart, realPart)
import Data.List               (nub, sort)

-- | returns the real roots of quartic polynomial equation in ascending order.
--
-- It does not return zero roots. Thus, if the roots are empty, it implies
-- that the polynomial has only zero roots.
quarticEqSol :: (Double -> Double)
             -> [Double]               -- ^ four distinct nonzero inputs
             -> Double                 -- ^ the cut value for complex roots
             -> Maybe (Int, [Double])  -- ^ (the number of real roots, roots)
quarticEqSol f xs eps
    | length (filter (/= 0) xs) /= 4 || xs /= nub xs = Nothing
    | otherwise = do
          let roots = filter (\r -> abs (imagPart r) < eps) $
                      polySolve (coeffQuartic f xs)
              realroots = (sort . filter ((> eps) . abs)) $ realPart <$> roots
          return (length realroots, realroots)

-- | returns the coefficients of quartic polynomial
--   using four distinct nonzero inputs.
--
--   a * x^4 + b * x^3 + c * x^2 + d * x + e
--
--   the output is the list of the coefficients [e, d, c, b, a].
coeffQuartic :: (Double -> Double) -> [Double] -> [Double]
coeffQuartic f xs =
    let [x1, x2, x3, x4] = xs
        [f1, f2, f3, f4] = f <$> xs

        d1 = x1 * (x2 - x1) * (x3 - x1) * (x4 - x1)
        d2 = x2 * (x2 - x1) * (x3 - x2) * (x4 - x2)
        d3 = x3 * (x3 - x1) * (x3 - x2) * (x4 - x3)
        d4 = x4 * (x4 - x1) * (x4 - x2) * (x4 - x3)

        p12 = x1 * x2
        p13 = x1 * x3
        p14 = x1 * x4
        p23 = x2 * x3
        p24 = x2 * x4
        p34 = x3 * x4
        p1234 = p12 * p34

        e = f 0

        a = - f1 / d1 + f2 / d2 - f3 / d3 + f4 / d4 + e / p1234

        b =   f1 * (x2 + x3 + x4) / d1
            - f2 * (x4 + x3 + x1) / d2
            + f3 * (x4 + x2 + x1) / d3
            - f4 * (x3 + x2 + x1) / d4
            - e * (x1 + x2 + x3 + x4) / p1234

        c = - f1 * (p34 + p24 + p23) / d1
            + f2 * (p34 + p14 + p13) / d2
            - f3 * (p24 + p14 + p12) / d3
            + f4 * (p23 + p13 + p12) / d4
            + e * (p34 + p24 + p14 + p23 + p13 + p12) / p1234

        d =   f1 * x2 * p34 / d1
            - f2 * x3 * p14 / d2
            + f3 * x4 * p12 / d3
            - f4 * x1 * p23 / d4
            - e * (x2 * p34 + x3 * p14 + x4 * p12 + x1 * p23) / p1234
    in
        [e, d, c, b, a]

-- | the real roots of quadratic equation: a x^2 + b x + c = 0.
--
-- If the root x0 is complex and |Im(x0)| < cut, it returns Re(x0).
quadEqSol :: Double  -- ^ the coefficient a
          -> Double  -- ^ the coefficient b
          -> Double  -- ^ the coefficient c
          -> Double  -- ^ the cut value for complex roots
          -> Maybe (Double, Double)
quadEqSol a b c epsScale = do
    let d = b * b - 4 * a * c
        eps = 1.0e-12
    if d >= 0
        -- then return ((-b + sqrt d) / (2*a), (-b - sqrt d) / (2*a))
        then do let q = -0.5 * (b + signum b * sqrt d)
                return (q / (a + eps), c / (q  + eps))
        else do let r = sqrt $ c / (a + eps)
                    th = acos $ - b / (2 * sqrt (a * c) + eps)
                    x = r * cos th
                if abs (r * sin th) < epsScale
                    then return (x, x)
                    else Nothing

max_stars_repo_name: cbpark/hep-vector

max_issues_repo_name: cbpark/hep-vector

max_forks_repo_name: cbpark/hep-vector


text:
{-# LANGUAGE NamedFieldPuns #-}
{-# LANGUAGE BlockArguments #-}
module AmpOp where

import Control.Monad.Fix

import Signal
import Circuit
import Data.Complex

type Input = Complex Double
type Output = Complex Double
type Value = Double
type Impedance = Signal (Complex Double)

dc12 :: Signal Input
dc12 = Signal $ const 12

dc14 :: Signal Input
dc14 = Signal $ const 14

dc100 :: Signal Input
dc100 = Signal $ const 100

sen :: Signal Input
sen = Signal $ \time -> 12 * sin ((pi / 2) * time)

outlet :: Signal Input
outlet = Signal $ \time -> 220 * sqrt 2 * sin (2 * pi * 60 * time)

ground :: Signal Input
ground = Signal $ const 0

data AmpOp = AmpOp { name :: String,
                     openLoopGain :: Signal (Complex Double),
                     manufacturer :: String
                }

data State =
  OneNode (Complex Double) |
  TwoNodes (Complex Double, Complex Double) |
  ThreeNodes (Complex Double, Complex Double, Complex Double) |
  FourNodes (Complex Double, Complex Double, Complex Double, Complex Double)

getOutput :: State -> Output
getOutput (OneNode x) = x
getOutput (TwoNodes (_,x)) = x
getOutput (ThreeNodes (_,_,x)) = x
getOutput (FourNodes (_,_,_,x)) = x

data SignalState =
  SignalOneNode (Signal (Complex Double)) |
  SignalTwoNodes (Signal (Complex Double), Signal (Complex Double)) |
  SignalThreeNodes (Signal (Complex Double), Signal (Complex Double), Signal (Complex Double)) |
  SignalFourNodes (Signal (Complex Double), Signal (Complex Double), Signal (Complex Double), Signal (Complex Double))
  deriving (Show)

getSignalOutput :: SignalState -> Signal Output
getSignalOutput (SignalOneNode x) = x
getSignalOutput (SignalTwoNodes (_,x)) = x
getSignalOutput (SignalThreeNodes (_,_,x)) = x
getSignalOutput (SignalFourNodes (_,_,_,x)) = x

instance Show AmpOp where
    show ap = Prelude.show "AmpOp => Name: " ++ Prelude.show (name ap)

lm741 :: AmpOp
lm741 = AmpOp "LM741" 200000 "Texas. Instruments"

r1 :: Impedance
r1 = 10000

r2 :: Impedance
r2 = 50000

r3 :: Impedance
r3 = 20000

r4 :: Impedance
r4 = 30000

eps = 0.000001

ampOpBuffer :: AmpOp -> SignalState -> Either String (Circuit Input SignalState)
ampOpBuffer model (SignalOneNode initial) =
    Right $ ($ initial) <$> mfix
    (\f -> Circuit $ \vIn -> Signal $ \time vOutOld ->
        let vOut = (vIn * openLoopGain model) / (1 + openLoopGain model)
        in if abs((realPart <$> vOut) - (realPart <$> vOutOld)) <= eps
            then SignalOneNode vOut
            else f vOut)
ampOpBuffer _ _ = Left "Wrong type for initial state!"

ampOpNonInverting :: AmpOp -> Impedance -> Impedance -> SignalState -> Either String (Circuit Input SignalState)
ampOpNonInverting model r1 r2 (SignalTwoNodes initial) =
    Right $ ($ initial) <$> mfix
    (\f -> Circuit $ \vIn -> Signal $ \time (vXOld, vOutOld) ->
        let vX = ((vIn * openLoopGain model) - vOutOld) / openLoopGain model
            vOut = (vXOld * (r1 + r2)) / r1
        in if abs((realPart <$> vOut) - (realPart <$> vOutOld)) <= eps && abs((realPart <$> vX) - (realPart <$> vXOld)) <= eps
            then SignalTwoNodes (vX, vOut)
            else f (vX, vOut))
ampOpNonInverting _ _ _ _ = Left "Wrong type for initial state!"

ampOpInverting :: AmpOp -> Impedance -> Impedance -> SignalState -> Either String (Circuit Input SignalState)
ampOpInverting model r1 r2 (SignalTwoNodes initial) =
    Right $ ($ initial) <$> mfix
    (\f -> Circuit $ \vIn -> Signal $ \time (vXOld, vOutOld) ->
        let vX = (- vOutOld) / openLoopGain model
            vOut = ((vXOld * (r1 + r2)) - (vIn * r2)) / r1
        in if abs((realPart <$> vOutOld) - (realPart <$> vOut)) <= eps && abs((realPart <$> vXOld) - (realPart <$> vX)) <= eps
            then SignalTwoNodes (vX, vOut)
            else f (vX, vOut))
ampOpInverting _ _ _ _ = Left "Wrong type for initial state!"

-- ampOpSpecial :: AmpOp -> Impedance -> Impedance -> Impedance -> Impedance -> SignalState -> Either String (Circuit Input SignalState)
-- ampOpSpecial model z1 z2 z3 z4 (SignalThreeNodes initial) =
--   Right $ ($ initial) <$> mfix
--   (\f -> Circuit $ \vIn -> Signal $ \time (vXOld, vKOld, vOutOld) ->
--       let vK = vOutOld + (vOutOld / openLoopGain model)
--           vX = ((z2 + z3) * vKOld) / z3
--           vOut = (((z2 * z4 + z1 * z4 + z1 * z2) * vXOld) - z1 * z4 * vKOld - z2 * z4 * vIn) / (z1 * z2)
--       in if abs((realPart <$> vOutOld) - (realPart <$> vOut)) <= eps
--           && abs((realPart <$> vKOld) - (realPart <$> vK)) <= eps
--           && abs((realPart <$> vXOld) - (realPart <$> vX)) <= eps
--           then SignalThreeNodes (vX, vK, vOut)
--           else f (vX, vK, vOut))

ampOpInvertingTest :: AmpOp -> Impedance -> Impedance -> SignalState -> Circuit Input SignalState
ampOpInvertingTest model r1 r2 (SignalTwoNodes initial) =
    ($ initial) <$> mfix
    (\f -> Circuit $ \vIn -> Signal $ \time (vXOld, vOutOld) ->
        let vX = (- vOutOld) / openLoopGain model
            vOut = ((vXOld * (r1 + r2)) - (vIn * r2)) / r1
        in if abs((realPart <$> vOutOld) - (realPart <$> vOut)) <= eps && abs((realPart <$> vXOld) - (realPart <$> vX)) <= eps
            then SignalTwoNodes (vX, vOut)
            else f (vX, vOut))

-- Calculating the output given a fixed input, to demonstrate the fix point only across the iteration axis
ampOpInverting' :: AmpOp -> Impedance -> Impedance -> Input -> State -> Either String Output
ampOpInverting' _ _ _ _ (OneNode _) = Left "Wrong type for initial state!"
ampOpInverting' model r1 r2 input (TwoNodes initial) = Right $ fix calculate initial
    where calculate f (vXOld, vOutOld) = if abs(realPart vOutOld - realPart vOut) <= 0.000001 && abs(realPart vXOld - realPart vX) <= 0.000001 then vOut else f (vX, vOut)
            where vX = (- vOutOld) / (openLoopGain model `at` 0)
                  vOut = (vXOld * (r1 + r2) `at` 0 - (input * r2 `at` 0)) / r1 `at` 0

max_stars_repo_name: FP-Modeling/fixingAnalog

max_issues_repo_name: FP-Modeling/fixingAnalog

max_forks_repo_name: FP-Modeling/fixingAnalog