text:
--------------------------------------------------------------------------------

{-# OPTIONS_GHC -fno-warn-orphans #-}

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

{-# LANGUAGE AllowAmbiguousTypes       #-}
{-# LANGUAGE BangPatterns              #-}
{-# LANGUAGE CPP                       #-}
{-# LANGUAGE EmptyDataDecls            #-}
{-# LANGUAGE FlexibleContexts          #-}
{-# LANGUAGE FlexibleInstances         #-}
{-# LANGUAGE ForeignFunctionInterface  #-}
{-# LANGUAGE FunctionalDependencies    #-}
{-# LANGUAGE GADTs                     #-}
{-# LANGUAGE KindSignatures            #-}
{-# LANGUAGE LambdaCase                #-}
{-# LANGUAGE MagicHash                 #-}
{-# LANGUAGE MultiParamTypeClasses     #-}
{-# LANGUAGE Rank2Types                #-}
{-# LANGUAGE ScopedTypeVariables       #-}
{-# LANGUAGE StandaloneDeriving        #-}
{-# LANGUAGE TypeApplications          #-}
{-# LANGUAGE TypeFamilyDependencies    #-}
{-# LANGUAGE TypeInType                #-}
{-# LANGUAGE UnboxedTuples             #-}
{-# LANGUAGE UndecidableInstances      #-}

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

-- | Internal module to Eigen.
--   Here we define all foreign function calls,
--   and some typeclasses integral to the public and private interfaces
--   of the library.
module Eigen.Internal where --   FIXME: Explicit export list

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

--import           Control.Monad            (when)
import           Data.Binary              (Binary(put,get))
import           Data.Binary.Get          (getByteString, getWord32be)
import           Data.Binary.Put          (putByteString, putWord32be)
import           Data.Bits                (xor)
import           Data.Complex             (Complex((:+)))
import           Data.Kind                (Type)
import           Data.Proxy               (Proxy(Proxy))
--import           Foreign.C.String         (CString, peekCString)
import           Foreign.C.Types          (CInt(CInt), CFloat(CFloat), CDouble(CDouble), CChar)
import           Foreign.ForeignPtr       (ForeignPtr, castForeignPtr, withForeignPtr)
import           Foreign.Ptr              (Ptr, castPtr, plusPtr)
import           Foreign.Storable         (Storable(sizeOf, alignment, poke, peek, peekByteOff, peekElemOff, pokeByteOff, pokeElemOff))
import           GHC.TypeLits             (natVal, KnownNat, Nat)
import           System.IO.Unsafe         (unsafeDupablePerformIO)
import qualified Data.Vector.Storable     as VS
import qualified Data.ByteString          as BS
import qualified Data.ByteString.Internal as BSI

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

-- | Like 'Proxy', but specialised to 'Nat'.
data Row (r :: Nat) = Row
-- | Like 'Proxy', but specialised to 'Nat'.
data Col (c :: Nat) = Col

-- | Used internally. Given a 'KnownNat' constraint, turn the type-level 'Nat' into an 'Int'.
natToInt :: forall n. KnownNat n => Int
{-# INLINE natToInt #-}
natToInt = fromIntegral (natVal @n Proxy)

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

-- | Cast to and from a C-FFI type
--   'Cast' is a closed typeclass with an associated injective type family.
--   It is closed in the sense that we provide only four types
--   with instances for it; and intend for eigen to only be used
--   with those four types. The injectivity of the type family is
--   then useful for avoiding MPTCs. 'Cast' has two functions; 'toC'
--   and 'fromC', where 'toC' goes from a Haskell type to its associated
--   C type for internal use, with the C FFI, and 'fromC' goes from the
--   associated C type to the Haskell type.
class Cast (a :: Type) where
  type family C a = (result :: Type) | result -> a
  toC   :: a -> C a
  fromC :: C a -> a

instance Cast Int where
  type C Int = CInt
  toC = CInt . fromIntegral
  {-# INLINE toC #-}
  fromC (CInt x) = fromIntegral x
  {-# INLINE fromC #-}

instance Cast Float where
  type C Float = CFloat
  toC = CFloat
  {-# INLINE toC #-}
  fromC (CFloat x) = x
  {-# INLINE fromC #-}

instance Cast Double where
  type C Double = CDouble
  toC = CDouble
  {-# INLINE toC #-}
  fromC (CDouble x) = x
  {-# INLINE fromC #-}

instance Cast a => Cast (Complex a) where
  type C (Complex a) = CComplex (C a)
  toC (a :+ b) = CComplex (toC a) (toC b)
  {-# INLINE toC #-}
  fromC (CComplex a b) = (fromC a) :+ (fromC b)
  {-# INLINE fromC #-}

-- | WARNING! 'toC' is lossy for any Int greater than (maxBound :: Int32)!
instance Cast a => Cast (Int, Int, a) where
  type C (Int, Int, a) = CTriplet a
  {-# INLINE toC #-}
  toC (x, y, z) = CTriplet (toC x) (toC y) (toC z)
  {-# INLINE fromC #-}
  fromC (CTriplet x y z) = (fromC x, fromC y, fromC z)

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

-- | Complex number for FFI with the same memory layout as std::complex\<T\>
data CComplex a = CComplex !a !a deriving (Show)

instance Storable a => Storable (CComplex a) where
    sizeOf _ = sizeOf (undefined :: a) * 2
    alignment _ = alignment (undefined :: a)
    poke p (CComplex x y) = do
        pokeElemOff (castPtr p) 0 x
        pokeElemOff (castPtr p) 1 y
    peek p = CComplex
        <$> peekElemOff (castPtr p) 0
        <*> peekElemOff (castPtr p) 1

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

-- | FIXME: Doc
data CTriplet a where
  CTriplet :: Cast a => !CInt -> !CInt -> !(C a) -> CTriplet a

deriving instance (Show a, Show (C a)) => Show (CTriplet a)

instance (Storable a, Elem a) => Storable (CTriplet a) where
    sizeOf _ = sizeOf (undefined :: a) + sizeOf (undefined :: CInt) * 2
    alignment _ = alignment (undefined :: CInt)
    poke p (CTriplet row col val) = do
        pokeElemOff (castPtr p) 0 row
        pokeElemOff (castPtr p) 1 col
        pokeByteOff p (sizeOf (undefined :: CInt) * 2) val
    peek p = CTriplet
        <$> peekElemOff (castPtr p) 0
        <*> peekElemOff (castPtr p) 1
        <*> peekByteOff p (sizeOf (undefined :: CInt) * 2)

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

-- | `Elem` is a closed typeclass that encompasses the properties
--   eigen expects its values to possess, and simplifies the external
--   API quite a bit.
class (Num a, Cast a, Storable a, Storable (C a), Code (C a)) => Elem a

instance Elem Float
instance Elem Double
instance Elem (Complex Float)
instance Elem (Complex Double)

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

-- | Encode a C Type as a CInt
--
--   Hack used in FFI wrapper functions when constructing FFI calls
class Code a where; code :: a -> CInt
instance Code CFloat             where; code _ = 0
instance Code CDouble            where; code _ = 1
instance Code (CComplex CFloat)  where; code _ = 2
instance Code (CComplex CDouble) where; code _ = 3

-- | Hack used in constructing FFI calls.
newtype MagicCode = MagicCode CInt deriving Eq

instance Binary MagicCode where
    put (MagicCode _code) = putWord32be $ fromIntegral _code
    get = MagicCode . fromIntegral <$> getWord32be

-- | Hack used in constructing FFI calls.
magicCode :: Code a => a -> MagicCode
magicCode x = MagicCode (code x `xor` 0x45696730)

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

-- | Machine size of a 'CInt'.
intSize :: Int
intSize = sizeOf (undefined :: CInt)

-- | FIXME: Doc
encodeInt :: CInt -> BS.ByteString
encodeInt x = BSI.unsafeCreate (sizeOf x) $ (`poke` x) . castPtr

-- | FIXME: Doc
decodeInt :: BS.ByteString -> CInt
decodeInt (BSI.PS fp fo fs)
    | fs == sizeOf x = x
    | otherwise = error "decodeInt: wrong buffer size"
    where x = performIO $ withForeignPtr fp $ peek . (`plusPtr` fo)

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

-- | 'Binary' instance for 'Data.Vector.Storable.Mutable.Vector'
instance Storable a => Binary (VS.Vector a) where
    put vs = put (BS.length bs) >> putByteString bs where
        (fp,fs) = VS.unsafeToForeignPtr0 vs
        es = sizeOf (VS.head vs)
        bs = BSI.fromForeignPtr (castForeignPtr fp) 0 (fs * es)

    get = get >>= getByteString >>= \bs -> let
        (fp,fo,fs) = BSI.toForeignPtr bs
        es = sizeOf (VS.head vs)
        -- `plusForeignPtr` is used qualified here to just remind a reader
        -- that it is defined internally within eigen
        vs = VS.unsafeFromForeignPtr0 (Eigen.Internal.plusForeignPtr fp fo) (fs `div` es)
        in return vs

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

-- | FIXME: Doc
data CSparseMatrix a
-- | FIXME: Doc
type CSparseMatrixPtr a = Ptr (CSparseMatrix a)

-- | FIXME: Doc
data CSolver a
-- | FIXME: Doc
type CSolverPtr a = Ptr (CSolver a)

-- {-# INLINE unholyPerformIO #-}
-- unholyPerformIO :: IO a -> a
-- unholyPerformIO (IO m) = case m realWorld# of (# _, r #) -> r

-- | FIXME: replace with unholyPerformIO (?)
performIO :: IO a -> a
performIO = unsafeDupablePerformIO

-- | FIXME: Doc
plusForeignPtr :: ForeignPtr a -> Int -> ForeignPtr b
plusForeignPtr fp fo = castForeignPtr fp1 where
    vs :: VS.Vector CChar
    vs = VS.unsafeFromForeignPtr (castForeignPtr fp) fo 0
    (fp1, _) = VS.unsafeToForeignPtr0 vs

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

ext: hsc

max_stars_repo_name: hasanalrasyid/inline-fortran

max_issues_repo_name: hasanalrasyid/inline-fortran

max_forks_repo_name: hasanalrasyid/inline-fortran


text:
{-
module Solve
  ( Ord
  , DoubleMap
  , innerKeys
  , toMatrix
  , toMatrix'
  , trans
  , solve
  , smartSolve ) where
-}
module Solve where

import           Data.List             (elemIndex, nub, sort, union, (\\))
import           Data.List.Extra       (groupSortOn)
import           Data.Map.Strict       (Map)
import qualified Data.Map.Strict       as M
import           Data.Maybe            (catMaybes, mapMaybe)
import           Data.Tuple            (swap)
import           Numeric.LinearAlgebra

type DoubleMap k l a = Map k (Map l a)

pullMaybe :: Monad m => (m a, m b) -> m (a, b)
pullMaybe = (>>= (fmap swap . sequence . swap)) . sequence

innerKeys :: Eq l => DoubleMap k l a -> [l]
innerKeys = nub . concatMap M.keys . M.elems

keyToInt :: Ord k => k -> DoubleMap k l a -> Maybe Int
keyToInt n = elemIndex n . sort . M.keys

innerKeyToInt :: Ord l => l -> DoubleMap k l a -> Maybe Int
innerKeyToInt n = elemIndex n . sort . innerKeys

toAssocList :: DoubleMap k l a -> [((k, l), a)]
toAssocList = fmap f . concatMap sequence . M.toList . fmap M.toList
  where f (n, (i, a)) = ((n, i), a)

toAssocInt :: (Ord k, Ord l) =>  DoubleMap k l a -> [((Int, Int), a)]
toAssocInt b = mapMaybe f . toAssocList $ b
  where
    f ((n, i), a) = pullMaybe (pullMaybe (keyToInt n b, innerKeyToInt i b), Just a)

fromAssocList :: (Ord k, Ord l) => [((k, l), a)] -> DoubleMap k l a
fromAssocList = fmap M.fromList . M.fromList . fmap shift . groupSortOn fst . fmap f
  where
    shift xs = (fst . head $ xs, fmap snd xs)
    f ((k, l), a) = (k, (l, a))

trans :: (Ord k, Ord l) => DoubleMap k l a -> DoubleMap l k a
trans = fromAssocList . fmap f . toAssocList
  where f ((l, k), a) = ((k, l), a)

-- Outer keys are rows
toMatrix :: (Container Vector a, Num a, Ord k, Ord l)
  => DoubleMap k l a -> Matrix a
toMatrix b = assoc (M.size b, length $ innerKeys b) 0 (toAssocInt b)

-- Column matrix from a Map
toMatrix' :: (Container Vector a, Num a, Ord k) => Map k a -> Matrix a
toMatrix' = toMatrix . fmap (M.fromList . pure . ((),))

-- Are the rows independent?
isInd :: (Field a, Ord k, Ord l) => DoubleMap k l a -> Bool
isInd m = rank (toMatrix m) == M.size m

-- Choices of n elements from a list, fast algorithm from stack overflow.
-- https://stackoverflow.com/questions/21265454/subsequences-of-length-n-from-list-performance/21288092#21288092
choices :: Int -> [a] -> [[a]]
choices n xs =
  let l = length xs
  in if n>l then [] else subsequencesBySize xs !! (l-n)
  where
    subsequencesBySize [] = [[[]]]
    subsequencesBySize (y:ys) =
      let next = subsequencesBySize ys
      in zipWith (++) ([]:next) (map (map (y:)) next ++ [[]])

-- SubMaps with n keys, excluding some.
allSizeNSub :: Ord k => Int -> Map k v -> [Map k v]
allSizeNSub n m =
  [ M.filterWithKey (\k _ -> k `elem` keyList) m
  | keyList <- choices n . M.keys $ m ]

-- The max rank submatricies with independent rows, excluding some.
allMaxInd :: (Field a, Ord k, Ord l)
  => DoubleMap k l a -> [DoubleMap k l a]
allMaxInd m = filter isInd . allSizeNSub (rank $ toMatrix m) $ m

-- The max rank submatricies with independent rows,
-- including or excluding certain rows.
chooseInd :: (Field a, Ord k, Ord l)
  => [k] -> [k] -> DoubleMap k l a -> [DoubleMap k l a]
chooseInd fixedKeys free m = filter hasKeys . allMaxInd $ m'
  where
    hasKeys b = null $ fixedKeys \\ M.keys b
    m' = M.filterWithKey (\k _ -> k `notElem` free) m

-- Maximal sets of independent columns.
chooseCols :: (Field a, Ord k, Ord l) => DoubleMap k l a -> [DoubleMap k l a]
chooseCols = fmap trans . allMaxInd . trans

-- All submatricies with independent rows and columns, respecting row preferences.
chooseBoth :: (Field a, Ord k, Ord l)
  => [k] -> [k] -> DoubleMap k l a -> [DoubleMap k l a]
chooseBoth fixedKeys free m = concatMap chooseCols (chooseInd fixedKeys free m)

-- All legal constraints where extra terms to zero.
allExtraZero :: (Field a, Ord k, Num v, Ord l)
  => Map k v -> [k] -> DoubleMap k l a -> [Map k v]
allExtraZero fixed free = fmap (M.union fixed . fmap (const 0))
  . chooseInd (M.keys fixed) free

-- Solve Ax = B, where B is a column matrix.
innerSolver :: (Field a, Ord l, Ord k)
  => DoubleMap k l a -> Map k a -> Maybe (Map l a)
innerSolver b c
  | M.keys b == M.keys c = fromCol <$> linearSolve (toMatrix b) (toMatrix' c)
  | otherwise = Nothing
  where
    fromCol = M.fromList . zip (sort . innerKeys $ b) . concat . toLists

-- Takes a list of fixed values, a list of free values,
-- and a DoubleMap, and attempts to solve.
solve :: (Eq a, Field a, Ord k, Ord l)
  => Map k a -> [k] -> DoubleMap k l a -> [Map l a]
solve fixed free m = nub . catMaybes $
  [ innerSolver x y
  | x <- chooseBoth (M.keys fixed) free m
  , y <- allExtraZero fixed free m ]

-- Rows with a unique nonzero entry.
looseEnds :: (Eq a, Ord k, Num a) => DoubleMap k l a -> [k]
looseEnds m =
  [ k
  | k <- M.keys m
  , let rowK = M.findWithDefault M.empty k m
  , (== 1) . M.size . M.filter (/= 0) $ rowK ]

-- Tries setting loose ends to be free,
-- excluding those already declared to be fixed.
smartSolve :: (Field a, Ord k, Ord l, Ord a)
  => Map k a -> [k] -> DoubleMap k l a -> [Map l a]
smartSolve fixed free m = case filter (all (>= 0)) . solve fixed free' $ m of
  [] -> solve fixed free m
  xs -> xs
  where
    free' = union free $ looseEnds m \\ M.keys fixed

ext: hs

max_stars_repo_name: alexloomis/factory-planner

max_issues_repo_name: alexloomis/factory-planner

max_forks_repo_name: alexloomis/factory-planner


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE OverloadedLists #-}
{-# LANGUAGE CPP #-}

module Main where

import qualified Data.Vector.Unboxed         as U
import           Data.Vector.Unboxed         (Vector)
import qualified Data.Vector as V

-- DLA
import qualified Statistics.Matrix as M
import qualified Statistics.Matrix.Fast as MF
import qualified Statistics.Matrix.Fast.Algorithms as A

-- hmatrix
import qualified Numeric.LinearAlgebra as H

-- numhask
import qualified NumHask.Array as NH
import qualified NumHask.Prelude as NP

-- data.matrix from matrix
import qualified Data.Matrix as DMX

-- massiv
import qualified Data.Massiv.Array as MA

import qualified System.Random.MWC as Mwc

import qualified Criterion.Main as C

#define N 10
#define N2 100

n :: Int
n = N

vectorGen :: IO (Vector Double)
vectorGen =  do
    gen <- Mwc.create
    Mwc.uniformVector gen (n*n)

matrixDLA :: IO M.Matrix
matrixDLA = do
    vec <- vectorGen
    return $ M.Matrix n n vec

matrixH :: IO (H.Matrix H.R)
matrixH = do
    vec <- vectorGen
    return $ (n H.>< n) $ U.toList $ vec

identH :: Int -> H.Matrix Double
identH = H.ident

identDMX :: Int -> DMX.Matrix Double
identDMX = DMX.identity

elemZero :: Double -> Double
elemZero = const 0

elemSqr :: Double -> Double
elemSqr x = x*x

mapH :: (Double -> Double) -> H.Matrix Double -> H.Matrix Double
mapH = H.cmap

main :: IO ()
main = do

    vDLA' <- vectorGen
    uDLA' <- vectorGen

    let

    --
      subDLA' = U.take n vDLA'
      aDLA' = M.Matrix n n vDLA'
      bDLA' = M.Matrix n n uDLA'

    --
      vList = U.toList vDLA'
      uList = U.toList uDLA'

    --
      aH' = (n H.>< n) vList
      bH' = (n H.>< n) uList

      subH' = H.fromList . take n $ vList
      vH' = H.fromList vList

    --
      aNH' = NP.fromList vList :: NH.Array V.Vector '[N, N] Double
      bNH' = NP.fromList uList :: NH.Array V.Vector '[N, N] Double

      vNH' = NP.fromList vList :: NH.Array V.Vector '[N2] Double

    --
      aDMX' = DMX.fromList n n vList
      bDMX' = DMX.fromList n n uList

    --
      vMA' = MA.fromList MA.Seq vList :: MA.Array MA.P MA.Ix1 Double
      aMA' = MA.resize' (MA.Sz (n MA.:. n)) vMA' :: MA.Array MA.P MA.Ix2 Double
      bMA' = MA.resize' (MA.Sz (n MA.:. n)) $ MA.fromList MA.Seq uList :: MA.Array MA.P MA.Ix2 Double

    C.defaultMain [
        C.env (pure (aDLA', bDLA', subDLA', vDLA')) $ \ ~(aDLA, bDLA, subDLA, vDLA) ->
            C.bgroup "DLA" [
                         C.bench "multiplication" $ C.nf (MF.multiply aDLA) bDLA,
                         C.bench "repeated multiplication" $ C.nf (U.sum . (flip M.row) 1 . MF.multiply bDLA . MF.multiply aDLA . MF.multiply aDLA ) bDLA,
                         C.bench "multiplicationV" $ C.nf (MF.multiplyV aDLA) subDLA,
                         C.bench "qr factorization" $ C.nf A.qr aDLA,
                         C.bench "transpose" $ C.nf MF.transpose aDLA,
                         C.bench "norm" $ C.nf MF.norm vDLA,
                         C.bench "row" $ C.nf (M.row  aDLA) 0,
                         C.bench "column" $ C.nf (M.column  aDLA) 0,
                         C.bench "identity" $ C.nf M.ident n,
                         C.bench "diag" $ C.nf M.diag subDLA,
                         C.bench "map const 0" $ C.nf (M.map elemZero) aDLA,
                         C.bench "map sqr" $ C.nf (M.map elemSqr) aDLA
                       ],

        C.env (pure (aH', bH', subH', vH')) $ \ ~(aH, bH, subH, vH) ->
            C.bgroup "Hmatrix" [
                             C.bench "multiplication" $ C.nf ((<>) aH) bH,
                             C.bench "repeated multiplication" $ C.nf ( H.sumElements . flip (H.?) [1] . (<>) bH . (<>) aH . (<>) aH) bH,
                             C.bench "multiplicationV" $ C.nf ((H.#>) aH) subH,
                             C.bench "qr factorization" $ C.nf H.qr aH,
                             C.bench "transpose" $ C.nf H.tr aH,
                             C.bench "norm" $ C.nf H.norm_2 vH,
                             C.bench "row" $ C.nf ((H.?) aH) [0],
                             C.bench "column" $ C.nf ((H.¿) aH) [0],
                             C.bench "identity" $ C.nf identH n,
                             C.bench "diag" $ C.nf H.diag subH,
                             C.bench "map const 0" $ C.nf (mapH elemZero) aH,
                             C.bench "map sqr" $ C.nf (mapH elemSqr) aH
                           ],

        C.env (pure (aNH', bNH', vNH')) $ \ ~(aNH, bNH, vNH) ->
            C.bgroup "NumHask" [
                             C.bench "multiplication" $ C.nf (NH.mmult aNH) bNH,
                             C.bench "repeated multiplication" $ C.nf ( (\(NH.Array a) -> V.sum a) . NH.row (NP.Proxy :: NP.Proxy 1) . NH.mmult bNH . NH.mmult aNH . NH.mmult aNH ) bNH,
                             C.bench "transpose" $ C.nf NH.transpose aNH,
                             C.bench "norm" $ C.nf (sqrt . (NP.<.> vNH)) vNH,
                             C.bench "row" $ C.nf (NH.row (NP.Proxy :: NP.Proxy 0)) aNH,
                             C.bench "column" $ C.nf (NH.col (NP.Proxy :: NP.Proxy 0)) aNH
                           ],

        C.env (pure (aMA', bMA', vMA')) $ \ ~(aMA, bMA, vMA) ->
            C.bgroup "Massiv" [
                                 C.bench "multiplication" $ C.nf ((MA.|*|) aMA) bMA,
                                 C.bench "multiplication (Par)" $ C.nf ((MA.|*|) (MA.setComp MA.Par aMA)) bMA,
                                 C.bench "repeated multiplication" $ C.nf ( MA.foldlS (+) 0 .  flip (MA.!>) 1 . (MA.|*|) bMA . (MA.|*|) aMA . (MA.|*|) aMA) bMA,
                                 C.bench "repeated multiplication (Par)" $ C.nf ( MA.foldlS (+) 0 .  flip (MA.!>) 1 . (MA.|*|) bMA . (MA.|*|) aMA . (MA.|*|) (MA.setComp MA.Par aMA)) bMA,
                                 C.bench "norm" $ C.nf (sqrt . MA.foldlS (+) 0 . (MA.zipWith (*) vMA)) vMA,
                                 C.bench "transpose" $ C.nf (MA.computeAs MA.P . MA.transpose) aMA,
                                 C.bench "row" $ C.nf (MA.computeAs MA.P . (MA.!>) aMA) 0,
                                 C.bench "column" $ C.nf (MA.computeAs MA.P . (MA.<!) aMA) 0
                                 ],

        C.env (pure (aDMX', bDMX')) $ \ ~(aDMX, bDMX) ->
                C.bgroup "Matrix" [
                                 C.bench "multiplication" $ C.nf (DMX.multStrassenMixed aDMX) bDMX,
                                 C.bench "transpose" $ C.nf DMX.transpose aDMX,
                                 C.bench "row" $ C.nf (DMX.getRow 1) aDMX,
                                 C.bench "column" $ C.nf (DMX.getCol 1) aDMX,
                                 C.bench "identity" $ C.nf identDMX n
                                 ]
                  ]

max_stars_repo_name: ocramz/fastest-matrices

max_issues_repo_name: ocramz/fastest-matrices

max_forks_repo_name: ocramz/fastest-matrices


text:
module GPLVM.Types
       ( Covariance
       , Distribution
       , EigenVectors
       , InputObservations (..)
       , LatentSpacePoints (..)
       , Matrix (..)
       , Mean (..)
       , ObservedData (..)
       , Vector (..)
       , unLatentSpacePoints
       , unObservedData
       , unInputObs
       ) where

import Universum hiding (Vector, transpose)

import Control.Lens (makeLenses)

import Data.Array.Repa
import Data.Array.Repa.Repr.ForeignPtr
-- import Foreign.Storable
import Numeric.LinearAlgebra.Repa hiding (Matrix, Vector)

type Distribution r a = Matrix r a -> a

type Mean r a = Matrix r a -> a

type Covariance a = Herm a

type Matrix r a = Array r DIM2 a

type Vector r a = Array r DIM1 a

type EigenVectors a = Matrix D a

-- Real matrix, represented as unboxed vector.
-------

newtype LatentSpacePoints a = LatentSpacePoints { _unLatentSpacePoints :: Matrix D a }
    deriving Eq

makeLenses ''LatentSpacePoints

newtype InputObservations a = InputObservations { _unInputObs :: Vector D a }
    deriving Eq

newtype ObservedData a = ObservedData { _unObservedData :: Matrix D a }
    deriving Eq

makeLenses ''ObservedData
makeLenses ''InputObservations

max_stars_repo_name: serokell/GPLVMHaskell

max_issues_repo_name: serokell/GPLVMHaskell

max_forks_repo_name: serokell/GPLVMHaskell


text:
#!/usr/bin/env stack
-- stack --resolver lts-15.04 runghc --package reanimate
module Main where

import           Codec.Picture.Types
import           Control.Lens
import           Control.Monad
import           Data.List
import           Data.Maybe
import qualified Data.Text                     as T
import           Data.Tuple
import qualified Data.Vector                   as V
import           Linear.V2
import           Linear.Vector
import qualified Numeric.LinearAlgebra         as Matrix
import           Numeric.LinearAlgebra.HMatrix (Matrix, linearSolve, toLists,
                                                (><))
import           Reanimate
import           Reanimate.Math.Common
import           Reanimate.Math.Polygon

bgColor :: PixelRGBA8
bgColor = PixelRGBA8 252 252 252 0xFF

type Points = V.Vector (V2 Double)
type Edges = [(Int, Int, Int)]
data Mesh = Mesh { meshPoints :: Points, meshEdges :: Edges }
data MeshPair = MeshPair Points Points Edges
-- The points in a RelMesh are:
--   relMeshStatic ++ x where Ax = B
data RelMesh = RelMesh
  { relMeshStatic :: Points
  , relMeshEdges  :: Edges
  , relMeshA      :: Matrix Double
  , relMeshB      :: Matrix Double
  }
data RelMeshPair = RelMeshPair Points Edges (Matrix Double) (Matrix Double) (Matrix Double) (Matrix Double)
-- Linear interpolation on RelMesh gives smooth morph.
-- solveMesh :: RelMesh -> Mesh
-- mkRelative :: Mesh -> RelMesh
-- mkRelativePair :: MeshPair -> RelMeshPair
-- triangulate :: Polygon -> Polygon -> MeshPair ?
-- embed :: MeshPair -> MeshPair
-- compatible :: Mesh -> Mesh -> Maybe MeshPair
-- linearInterpolate :: MeshPair -> Double -> Mesh
-- convexInterpolate :: RelMeshPair -> Double -> RelMesh

meshToPolygon :: Mesh -> Polygon
meshToPolygon mesh = pScale 2 $ mkPolygon $ V.fromList
    [ realToFrac <$> (meshPoints mesh V.! (n-1))
    | n <- polygonNodes ]
  where
    polygonNodes :: [Int]
    polygonNodes = [5,8,9,7,4,6]

main :: IO ()
main = reanimate morphAnimation

genTrails :: (Double -> Mesh) -> SVG
genTrails mkMesh =
    withFillOpacity 0 $
    withStrokeWidth (defaultStrokeWidth*0.5) $
    withStrokeColor "black" $
    withStrokeDashArray [0.1,0.1] $
    mkGroup $ map mkTrail $ transpose
      [ V.toList $ V.map (fmap realToFrac) $ polygonPoints $ meshToPolygon mesh
      | n <- [0..steps]
      , let mesh = mkMesh (fromIntegral n / fromIntegral steps)
      ]
  where
    steps = 100 :: Int
    mkTrail lst =mkLinePath [ (x,y) | V2 x y <- lst ]

morphAnimation :: Animation
morphAnimation = playThenReverseA $ pauseAround 1 1 $ addStatic (mkBackgroundPixel bgColor) $
  signalA (curveS 2) $ mkAnimation 5 $ \t -> lowerTransformations $ scale 3 $ pathify $
  mkGroup
  [ translate (-1.5) 0 $ withFillColor "lightgreen" $ mkGroup
    [ withStrokeWidth defaultStrokeWidth $
      withStrokeColor "black" $
      polygonShape $ meshToPolygon $ linearInterpolate meshPair t
    , linearTrails
    , polygonNumDots (meshToPolygon $ linearInterpolate meshPair t) 0
    ]
  , translate 1 0 $ withFillColor "cyan" $ mkGroup
    [ withStrokeWidth defaultStrokeWidth $
      withStrokeColor "black" $
      polygonShape $ meshToPolygon $ solveMesh $ convexInterpolate relPair t
    , curveTrails
    , polygonNumDots (meshToPolygon $ solveMesh $ convexInterpolate relPair t) 0
    ]
  ]
  where
    linearTrails = genTrails (linearInterpolate meshPair)
    curveTrails = genTrails (solveMesh . convexInterpolate relPair)
    relPair = mkRelativePair meshPair
    meshPair = fromJust $ compatible example1 example2

mkLineP :: P -> P -> SVG
mkLineP (V2 x1 y1) (V2 x2 y2) = mkLine (x1,y1) (x2,y2)

-- FIXME: Check that the triangles are all anticlockwise.
-- FIXME: Check that the edges connect all the points.
-- FIXME: Check that the edgse leave no gaps.
compatible :: Mesh -> Mesh -> Maybe MeshPair
compatible a b =
  if meshEdges a == meshEdges b && V.length (meshPoints a) == V.length (meshPoints b)
    then Just $ MeshPair (meshPoints a) (meshPoints b) (meshEdges a)
    else Nothing

linearInterpolate :: MeshPair -> Double -> Mesh
linearInterpolate (MeshPair aP bP edges) t = Mesh
    { meshPoints = V.zipWith (lerp (1-t)) aP bP
    , meshEdges = edges }

example1 :: Mesh
example1 = Mesh points edges
  where
    points = V.fromList
      [ V2 1 0
      , V2 (-1/2) (sqrt 3 / 2)
      , V2 (-1/2) (-sqrt 3 / 2)
      , 3 * points V.! 0 ^/ 4
      , 3 * points V.! 1 ^/ 4
      , 3 * points V.! 2 ^/ 4
      , points V.! 0 ^/ 2
      , points V.! 1 ^/ 2
      , points V.! 2 ^/ 2
      ]
    edges =
      [ (1,5,4), (1,2,5), (2,6,5), (2,3,6), (3,4,6), (3,1,4)
      , (4,8,7), (4,5,8), (5,9,8), (5,6,9), (6,7,9), (6,4,7), (7,8,9)]

example2 :: Mesh
example2 = Mesh points edges
  where
    points = V.fromList
      [ V2 1 0
      , V2 (-1/2) (sqrt 3 / 2)
      , V2 (-1/2) (-sqrt 3 / 2)
      , 3 * points V.! 2 ^/ 4
      , 3 * points V.! 0 ^/ 4
      , 3 * points V.! 1 ^/ 4
      , points V.! 1  ^/ 2
      , points V.! 2 ^/ 2
      , points V.! 0 ^/ 2
      ]
    edges = meshEdges example1

-- T = (U, G)
-- G = [Polygon]
-- U = nub $ concat G

findStarNeighbours :: Eq a => [(a,a,a)] -> a -> [(a, a)]
findStarNeighbours allTrig self =
  [ (b,c)
  | (a,b,c) <- allTrig
  , self == a
  ] ++
  [ (c,a)
  | (a,b,c) <- allTrig
  , self == b
  ] ++
  [ (a,b)
  | (a,b,c) <- allTrig
  , self == c
  ]

isInterior :: Eq a => [(a, a)] -> Bool
isInterior = isJust . getExteriorPoly

getExteriorPoly :: Eq a => [(a, a)] -> Maybe [a]
getExteriorPoly [] = Nothing
getExteriorPoly ((a,b):rest) = worker [a] a b rest
  where
    worker acc start this [] = do
      guard (start == this)
      return (reverse acc)
    worker acc start this xs =
      case lookup this xs of
        Just next -> worker (this:acc) start next (delete (this,next) xs)
        Nothing   ->
          case lookup this (map swap xs) of
            Just next -> worker (this:acc) start next (delete (next, this) xs)
            Nothing   -> Nothing

convexInterpolate :: RelMeshPair -> Double -> RelMesh
convexInterpolate (RelMeshPair static edges leftM leftB rightM rightB) t =
  RelMesh
  { relMeshStatic = static
  , relMeshEdges  = edges
  , relMeshA      = Matrix.scale (1-t) leftM +
                    Matrix.scale t rightM
  , relMeshB      = Matrix.scale (1-t) leftB +
                    Matrix.scale t rightB
  }

solveMesh :: RelMesh -> Mesh
solveMesh (RelMesh static edges m b) =
  case linearSolve m b of
    Nothing -> error "Failed to solve mesh"
    Just ret ->
      Mesh (static <> V.fromList (worker (toLists ret))) edges
  where
    worker []             = []
    worker ([x]:[y]:rest) = V2 x y : worker rest
    worker _              = error "invalid result"

mkRelative :: Mesh -> RelMesh
mkRelative (Mesh points edges) = RelMesh (V.fromList exteriorPoints) edges mM bM
  where
    mM = (s><s) (concat m)
    bM = (s><1) b
    (s,exterior, (m, b)) = toParameters points edges
    exteriorPoints =
      [ points V.! (i-1)
      | i <- exterior
      ]

mkRelativePair :: MeshPair -> RelMeshPair
mkRelativePair (MeshPair p1 p2 edges) =
  let RelMesh static _ leftM leftB = mkRelative (Mesh p1 edges)
      RelMesh _ _ rightM rightB = mkRelative (Mesh p2 edges)
  in RelMeshPair static edges leftM leftB rightM rightB

toParameters :: (Ord b, Fractional b) => V.Vector (V2 b) -> [(Int, Int, Int)] -> (Int, [Int], ([[b]], [b]))
toParameters points groups = (length interior*2,exterior,unzip $ concat
  [ let lst = [(if i == j then -1 else t)
              | j <- interior
              , let t = fromMaybe 0 $ lookup (i,j) lam_ij_cache
              ]
        pos = negate $ sum
           [ pj ^* t
           | j <- exterior
           , let t = fromMaybe 0 $ lookup (i,j) lam_ij_cache
                 pj = points V.! (j-1)
           ]
    in [ (dupX lst, pos ^. _x)
       , (dupY lst, pos ^. _y)]
  | i <- interior ])
  where
    lam_ij_cache = lam_ij points groups
    dupX []     = []
    dupX (x:xs) = x:0:dupX xs
    dupY []     = []
    dupY (x:xs) = 0:x:dupY xs
    (interior, exterior) =
      partition (isInterior . findStarNeighbours groups) [1 .. length points]

lam_ij :: (Ord b, Fractional b) => V.Vector (V2 b) -> [(Int, Int, Int)] -> [((Int, Int), b)]
lam_ij points groups =
  [ ((i, j), t)
  | i <- [1..V.length points]
  , (j, t) <- lam_j points groups i
  ]

lam_j :: (Ord b, Fractional b) => V.Vector (V2 b) -> [(Int, Int, Int)] -> Int -> [(Int, b)]
lam_j points groups p =
  [ (nP, sum [ t | (j,_k,t) <- mu, j == nP ] / fromIntegral (length nPoints))
  | let n = findStarNeighbours groups p
        nPoints = fromMaybe [] $ getExteriorPoly n
        mu = calcMu points groups p
  , nP <- nPoints ]

calcMu :: (Ord c, Fractional c) => V.Vector (V2 c) -> [(Int, Int, Int)] -> Int -> [(Int, Int, c)]
calcMu points groups p = concat
    [ [ (nP, nP, t1)
      , (a, nP, t2)
      , (b, nP, t3) ]
      -- (nP, a, b)
    | {-p <- [1..length points]-}
      let selfVert = points V.! (p-1)
          n = findStarNeighbours groups p
          nPoints :: [Int]
          nPoints = fromMaybe [] $ getExteriorPoly n
    , nP <- nPoints
    , let vert = points V.! (nP-1)
    , let line = (vert, selfVert)
    , let (a,b,aP,bP) = head $
            [ (_a,_b,_aP,_bP)
            | (_a,_b) <- n
            , let _aP = points V.! (_a-1)
                  _bP = points V.! (_b-1)
                  segment = (points V.! (_a-1), points V.! (_b-1))
            , case rayIntersect line segment of
                Nothing -> False
                Just u  -> isBetween u segment
            , _a /= nP
            , _b /= nP ]
    -- , b == (nPoints ++ nPoints) !! i
    , let (t1,t2,t3) = barycentricCoords vert aP bP selfVert
    ]

genColor :: Int -> Int -> PixelRGBA8
genColor n m =
    promotePixel $ parula (fromIntegral (n+offset) / fromIntegral (m+offset))
  where
    offset = 5

polygonNumDots :: Polygon -> Double -> SVG
polygonNumDots p t = mkGroup $ reverse
    [ mkGroup
      [ colored n $ withStrokeWidth (defaultStrokeWidth*0.5) $ withStrokeColor "black" $
        translate x y $ mkCircle circR
      , withFillColor "black" $
        translate x y $ ppNum n ]
    | n <- [0..pSize p-1]
    , let a = realToFrac <$> pAccess p n
          b = realToFrac <$> pAccess p (pNext p n)
          V2 x y = lerp t b a ]
  where
    circR = 0.1
    colored n =
      let c = genColor n (pSize p-1)
      in withFillColorPixel c
    ppNum n = scaleToHeight (circR*1.5) $ center $ latex $ T.pack $ "\\texttt{" ++ show n ++ "}"

polygonShape :: Polygon -> SVG
polygonShape p = mkLinePathClosed
  [ (x,y) | V2 x y <- map (fmap realToFrac) $ V.toList (polygonPoints p) ]

max_stars_repo_name: TristanCacqueray/reanimate

max_issues_repo_name: TristanCacqueray/reanimate

max_forks_repo_name: TristanCacqueray/reanimate


text:
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE ScopedTypeVariables #-}

module LinearRegression
  ( trainLinReg
  , trainAndTestLinReg
  , trainLinRegWNoise
  , avError
  , linRegAsInitialForPLA
  , createVectorY
  , createRandomPoints
  , flipLabels
  , linearRegressionWeight
  , linRegClassificationError
  , testLinRegWNoise
  , createMatrixX
  , makeTargetFunction
  , y
  ) where

import           Control.Monad                 (replicateM)
import           Data.Maybe
import           Data.Traversable              (for)
import           Numeric.LinearAlgebra.Data
import           Numeric.LinearAlgebra.HMatrix
import           PLA
import           System.Random

--------------------------------------------------------------
-- Solutions to homework 2 of "Learning from data" question 5-6
--------------------------------------------------------------
---------------------------------------------
-- Creating the line for the target function
--------------------------------------------
-- | takes 2 points and returns the slope of the line connecting them
m :: (R, R) -> (R, R) -> R
m (x1, x2) (z1, z2) = (x2 - z2) / (x1 - z1)

-- | takes 2 points and returns the y-intercept of the line connecting them
b :: (R, R) -> (R, R) -> R
b (x1, x2) (z1, z2) = x2 - x1 * (x2 - z2) / (x1 - z1)

-- | Takes 2 Points and returns the function of a line connecting them
line point1 point2 = \x -> m point1 point2 * x + b point1 point2

-----------------------------------
-- Creating random data points
------------------------------------
-- |Creates a random Point between (-1,-1) and (1,1)
createRandomPoint :: IO (R, R)
createRandomPoint = do
  rand1 <- randomRIO (-1, 1)
  rand2 <- randomRIO (-1, 1)
  return (rand1, rand2)

-- |Creates a list of n random Points between (-1,-1) and (1,1)
createRandomPoints :: Int -> IO [(R, R)]
createRandomPoints n = do
  randList1 <- for [1 .. n] $ \_ -> randomRIO (-1, 1)
  randList2 <- for [1 .. n] $ \_ -> randomRIO (-1, 1)
  return $ zip randList1 randList2

-----------------------------------------------------
-- Preparing data for processing by Linear Regression
------------------------------------------------------
-- | looks up the label of the point by applying the target function
y :: (R -> R) -> (R, R) -> R
y yLine (x1, x2) =
  if yLine x1 < x2
    then (-1)
    else 1

-- | Creates the vector of labels using the target function 'yLine' and a list of data points
createVectorY :: ((R, R) -> R) -> [(R, R)] -> Vector R
createVectorY target listOfPoints = vector $ map target listOfPoints

-- | Creates the matrix X from a list of data points
createMatrixX :: [(R, R)] -> Matrix R
createMatrixX listOfPoints = matrix 3 listOfNumbers
  where
    listOfNumbers = concat [[1, a, b] | (a, b) <- listOfPoints]

flipLabels :: Float -> Vector R -> IO (Vector R)
flipLabels noiselevel labelvector = do
  gen <- getStdGen
  let numLabels = size labelvector
  let randnums :: [Int] = randomRs (1, numLabels) gen
  let toFlip = take (floor (noiselevel * fromIntegral numLabels)) randnums
  let flipVector =
        vector
          [ if x `elem` toFlip
            then (-1)
            else 1
          | x <- [1 .. numLabels]
          ]
  let flippedLabels = labelvector * flipVector
  return flippedLabels

--------------------------------------------------------
-- Linear Regression algorithm
--------------------------------------------------------
-- | returns the weights determined by linear regression on the training data
linearRegressionWeight :: Matrix R -> Vector R -> Vector R
linearRegressionWeight matrixOfInputData vectorOfLabels =
  pinv matrixOfInputData #> vectorOfLabels

-- | Squared error
squaredError :: Matrix R -> Vector R -> Vector R -> R
squaredError matrixOfInputData vectorOfLabels weight =
  1 / fromIntegral n * dot v v
  where
    v = signum (matrixOfInputData #> weight) - vectorOfLabels
    n = size vectorOfLabels

-- | Classification error for linear Regression (i.e. fraction of misclassified points)
linRegClassificationError :: Matrix R -> Vector R -> Vector R -> R
linRegClassificationError matrixOfInputData vectorOfLabels weight =
  numberMisclassified / totalDatapoints
  where
    numberMisclassified =
      norm_0 $ signum (matrixOfInputData #> weight) - vectorOfLabels
    totalDatapoints = fromIntegral $ size vectorOfLabels

--------------------------------------------------
-- training and testing the Linear Regression model
--------------------------------------------------
-- | create the target function
makeTargetFunction :: IO (R -> R)
makeTargetFunction = do
  point1 <- createRandomPoint
  point2 <- createRandomPoint
  let target = line point1 point2
  return target

-- | training linear Regression on random dataset of n noisy points and returning (weights, in-Sample error, trainingpoints)
trainLinRegWNoise ::
     ((R, R) -> R)
  -> Int
  -> Float
  -> ([(R, R)] -> Matrix R)
  -> IO (Vector R, R, [(R, R)])
trainLinRegWNoise target n noiselevel transformFunction = do
  trainpoints <- createRandomPoints n
  let trainX = transformFunction trainpoints
  let trainY = createVectorY target trainpoints
  trainYNoisy <- flipLabels noiselevel trainY
  let weights = linearRegressionWeight trainX trainYNoisy
  let inSampleError = linRegClassificationError trainX trainYNoisy weights
  return (weights, inSampleError, trainpoints)

-- | training linear Regression on random dataset of n noiseless points and returning (weights, in-Sample error, trainingpoints)
trainLinReg ::
     (R -> R) -> Int -> ([(R, R)] -> Matrix R) -> IO (Vector R, R, [(R, R)])
trainLinReg target n transformFunction =
  trainLinRegWNoise (y target) n 0 transformFunction

-- | evaluates linear Regression result on test data
testLinReg :: (R -> R) -> Int -> ([(R, R)] -> Matrix R) -> Vector R -> IO R
testLinReg target n transformFunction weights =
  testLinRegWNoise (y target) n 0 transformFunction weights

-- | training linear Regression on random dataset of n noisy points and returning (weights, in-Sample error, trainingpoints)
testLinRegWNoise ::
     ((R, R) -> R) -> Int -> Float -> ([(R, R)] -> Matrix R) -> Vector R -> IO R
testLinRegWNoise target n noiselevel transformFunction weights = do
  testpoints <- createRandomPoints n
  let testX = transformFunction testpoints
  let testY = createVectorY target testpoints
  testYNoisy <- flipLabels noiselevel testY
  let outOfSampleError = linRegClassificationError testX testYNoisy weights
  return outOfSampleError

-- | performs training and testing
trainAndTestLinReg :: IO (R, R)
trainAndTestLinReg = do
  target <- makeTargetFunction
  (weights, inSampleError, _) <- trainLinReg target 100 createMatrixX
  outOfSampleError <- testLinReg target 1000 createMatrixX weights
  return (inSampleError, outOfSampleError)

---------------------------------------------------------
-- averaging several runs of the algorithm
---------------------------------------------------------
-- | average in a list of numbers
average :: [R] -> R
average xs = sum xs / fromIntegral (length xs)

-- | average of in-sample and out-of-sample error, calculated from a list of pairs and returned as a pair.
averageErrors :: [(R, R)] -> (R, R)
averageErrors pairs = (avFst, avSnd)
  where
    avFst = average $ map fst pairs
    avSnd = average $ map snd pairs

-- |pair of average in-sample and out-of-sample error in n runs of linear Regression on n different sets of points
avError :: Int -> IO ()
avError n =
  putStrLn . ("In-sample error, out-of-sample error: " ++) . show =<<
  fmap averageErrors (replicateM n trainAndTestLinReg)

----------------------------------------------
-- Using linear regression to obtain initial weights for perceptron
--------------------------------------------------
listToTuple :: [a] -> Maybe (a, a, a)
listToTuple (x:y:z:xs) = Just (x, y, z)
listToTuple _          = Nothing

-- | trains a perceptron wiht initial weights obtained by linear Regression
linRegAsInitialForPLA :: IO ()
linRegAsInitialForPLA = do
  target <- makeTargetFunction
  (initialweightsVec, _, trainPoints) <- trainLinReg target 10 createMatrixX
  let maybeInitialWeights = listToTuple $ toList initialweightsVec
  let initialWeights = fromMaybe (0, 0, 0) maybeInitialWeights
  (finalWeights, epochs) <-
    trainPLAWithInitial initialWeights trainPoints target
  putStr "Weights: "
  print finalWeights
  putStr "Epochs: "
  print epochs


text:
{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE ConstraintKinds #-}
module DSLsofMath.Algebra where
import qualified Data.Ratio
import qualified Prelude
import Prelude (Double, Rational, Int, Integer, Bool(..), otherwise,
                Foldable(foldr), (.), const, (==), (<), error)
import Data.Complex

-------------------------------
-- Classes
infixl 6 -
infixl 6 +

infixl 7 *
infixl 7 /

class Additive a where
  zero :: a
  (+) :: a -> a -> a

sum :: (Foldable t, Additive a) => t a -> a
sum = foldr (+) zero

times :: Additive a => Integer -> a -> a
times n0 = if n0 < 0  then  error "Algebra.Classes.times: negative number of times"
                      else  go n0
  where go 0 _ = zero
        go 1 x = x
        go n x = if r == 0 then twoy else x + twoy
          where (m,r) = n `Prelude.divMod` 2
                y = go m x
                twoy = y+y

(-) :: AddGroup a => a -> a -> a
x - y           =  x + negate y

class Additive a => AddGroup a where
  negate :: a -> a

mult :: AddGroup a => Integer -> a -> a
mult n x = if n < 0 then negate (times (negate n) x) else times n x

class Multiplicative a where
  one :: a
  (*) :: a -> a -> a

two :: (Additive a, Multiplicative a) => a
two = one+one

(^+) :: Multiplicative a => a -> Int -> a

x0 ^+ n0 = if n0 < 0  then  error "Algebra.Classes.^: negative exponent"
                      else  go n0 x0
  where go 0 _ = one
        go 1 x = x
        go n x = if r == 0 then y2 else x * y2
          where  (m,r) = n `Prelude.divMod` 2
                 y = go m x
                 y2 = y * y

type Ring a = (AddGroup a, Multiplicative a)

fromInteger :: Ring a => Integer -> a
fromInteger n = mult n one

class Multiplicative a => MulGroup a where
  {-# MINIMAL (recip | (/)) #-}
  recip :: a -> a
  recip x         =  one / x

  (/) :: a -> a -> a
  x / y           =  x * recip y

(^) :: MulGroup a => a -> Int -> a
a ^ b | b < 0      = recip (a ^+ (negate b))
      | otherwise  = (a ^+ b)

type Field a = (Ring a, MulGroup a)

fromRational :: Field a => Data.Ratio.Ratio Integer -> a
fromRational x  =  fromInteger (Data.Ratio.numerator x) / fromInteger (Data.Ratio.denominator x)

class Field a => Algebraic a where
  sqrt :: a -> a

-- normally it should be "Algebraic" instead of "Field" but we're lazy like that.
--   (Also Transcendental is a terrible name; taken from the "numeric prelude".)
class Field a => Transcendental a where
  pi  :: a
  exp :: a -> a
  sin :: a -> a
  cos :: a -> a

cosh, sinh :: Transcendental a => a -> a
cosh x = (exp x + exp (negate x))/two
sinh x = (exp x - exp (negate x))/two

---------------------------------
-- Instances

instance Additive        Int       where  (+)  = (Prelude.+);  zero  = 0
instance Additive        Integer   where  (+)  = (Prelude.+);  zero  = 0
instance Additive        Rational  where  (+)  = (Prelude.+);  zero  = 0
instance Additive        Double    where  (+)  = (Prelude.+);  zero  = 0

instance AddGroup        Int       where  negate = Prelude.negate
instance AddGroup        Integer   where  negate = Prelude.negate
instance AddGroup        Rational  where  negate = Prelude.negate
instance AddGroup        Double    where  negate = Prelude.negate

instance Multiplicative  Int       where  (*)  = (Prelude.*);  one   = 1
instance Multiplicative  Integer   where  (*)  = (Prelude.*);  one   = 1
instance Multiplicative  Rational  where  (*)  = (Prelude.*);  one   = 1
instance Multiplicative  Double    where  (*)  = (Prelude.*);  one   = 1

instance MulGroup        Rational  where  (/)  = (Prelude./);  recip = Prelude.recip
instance MulGroup        Double    where  (/)  = (Prelude./);  recip = Prelude.recip

lift0 ::  a         ->  (x->a)
lift1 :: (a->b)     ->  (x->a) -> (x->b)
lift2 :: (a->b->c)  ->  (x->a) -> (x->b) -> (x->c)
lift0 = const
lift1 = (.)
lift2 op2 f g = \x -> op2 (f x) (g x)

instance Additive a        => Additive        (x -> a) where (+)  = lift2 (+);  zero  = lift0 zero
instance Multiplicative a  => Multiplicative  (x -> a) where (*)  = lift2 (*);  one   = lift0 one

instance AddGroup a   => AddGroup   (x -> a)  where   negate  =  lift1 negate
instance MulGroup a   => MulGroup   (x -> a)  where   recip   =  lift1 recip
instance Algebraic a  => Algebraic  (x -> a)  where   sqrt    =  lift1 sqrt

instance Transcendental a => Transcendental (x -> a) where
   pi = lift0 pi;  sin = lift1 sin;  cos = lift1 cos;  exp = lift1 exp

instance Algebraic       Double where   sqrt = Prelude.sqrt
instance Transcendental  Double where
   pi = Prelude.pi;  sin = Prelude.sin;  cos = Prelude.cos;  exp = Prelude.exp

instance Additive a  => Additive        (Complex a) where  (+) = addC;  zero  = zeroC
instance Ring a      => Multiplicative  (Complex a) where  (*) = mulC;  one   = oneC
instance AddGroup a  => AddGroup        (Complex a) where  negate  = negateC
instance Field a     => MulGroup        (Complex a) where  recip   = recipC

addC :: Additive a => Complex a -> Complex a -> Complex a
addC (x :+ y) (x' :+ y') = (x + x') :+ (y+y')

negateC :: AddGroup a => Complex a -> Complex a
negateC (a :+ b) = negate a :+ negate b

mulC :: Ring a => Complex a -> Complex a -> Complex a
mulC (a :+ b) (a' :+ b') = (a * a' - b * b') :+ (a * b' + b * a')

toC :: Additive a => a -> Complex a
toC x = x :+ zero

zeroC :: Additive a => Complex a
zeroC = toC zero

oneC :: (Additive a, Multiplicative a) => Complex a
oneC = toC one

recipC (a :+ b) = (a / m)  :+ (negate b / m)
  where m = a*a + b*b

instance (Algebraic a, Prelude.RealFloat a) =>  Algebraic (Complex a) where
  sqrt = Prelude.sqrt

instance (Transcendental a) => Transcendental (Complex a) where
  pi = piC;  exp = expC;  sin = sinC;  cos = cosC

piC :: Transcendental a => Complex a
piC = toC pi

expC, sinC, cosC, sinhC, coshC :: Transcendental a => Complex a -> Complex a
expC   (x:+y)  =  expx  * cos y  :+ expx * sin y  where expx = exp x
sinC   (x:+y)  =  sin x * cosh y :+ cos x * sinh y
cosC   (x:+y)  =  cos x * cosh y :+ negate (sin x * sinh y)

sinhC  (x:+y)  =  cos y * sinh x :+ sin y * cosh x
coshC  (x:+y)  =  cos y * cosh x :+ sin y * sinh x

---------------------------------
-- For RebindableSyntax

ifThenElse :: Bool -> p -> p -> p
ifThenElse c a b = if c then a else b

-------------------------------
-- Typesetting aliases.

neg :: AddGroup a => a -> a
neg = negate   -- |neg| is used to typeset unary minus as a shorter dash, closer to its argument

frac :: MulGroup a => a -> a -> a
frac = (/)     -- |frac| is used to typeset a fraction (more compactly than |x / y|)

max_stars_repo_name: nicolabotta/DSLsofMath

max_issues_repo_name: nicolabotta/DSLsofMath

max_forks_repo_name: nicolabotta/DSLsofMath


text:
-- ch1 ex 1.25

import qualified Data.Vector as V

import Statistics.Sample.Histogram (histogram)
import Statistics.Sample (mean, varianceUnbiased, stdDev, range)

samples :: V.Vector Double
samples = V.fromList
  [  16, 224,  16,  80,  96, 536, 400,  80
  , 392, 576, 128,  56, 656, 224,  40,  32
  , 358, 384, 256, 246, 328, 464, 448, 716
  , 304,  16,  72,   8,  80,  72,  56, 608
  , 108, 194, 136, 224,  80,  16, 424, 264
  , 156, 216, 168, 184, 552,  72, 184, 240
  , 438, 120, 308,  32, 272, 152, 328, 480
  ,  60, 208, 340, 104,  72, 168,  40, 152
  , 360, 232,  40, 112, 112, 288, 168, 352
  ,  56,  72,  64,  40, 184, 264,  96, 224
  , 168, 168, 114, 280, 152, 208, 160, 176
  ]

hist :: (V.Vector Double, V.Vector Int)
hist = histogram 15 samples

main = do
  putStrLn $ "the range/4: " ++ show (range samples / 4)

  let av = mean samples
      sd = stdDev samples

  putStrLn $ "the standard deviation: " ++ show sd
  putStrLn $ "the mean value: " ++ show av

  let (bins, hs) = hist
  print hs

  let f k = V.length . V.filter (\x -> (av-k*sd)<x && x< (av+k*sd)) $ samples
  putStrLn $ "the sample size:    " ++ show (V.length samples)
  putStrLn $ "within one sigma:   " ++ show (f 1)
  putStrLn $ "within two sigma:   " ++ show (f 2)
  putStrLn $ "within three sigma: " ++ show (f 3)

max_stars_repo_name: raydsameshima/Stat

max_issues_repo_name: raydsameshima/Stat

max_forks_repo_name: raydsameshima/Stat


text:
module Main where

import Criterion.Main (defaultMain, bgroup, bench, nf, nfIO, whnf)
import qualified Data.Vector.Storable as V
import qualified Numeric.LinearAlgebra as LA
import qualified Numeric.Morpheus as M

main :: IO ()
main = do
  !a <- LA.rand 1000 1000
  !b <- LA.rand 1000 1000
  !x <- LA.flatten <$> LA.rand 1000 1000
  !y <- LA.flatten <$> LA.rand 1000 1000
  let !ata = LA.mTm a

  defaultMain [
    bgroup "Linear Algebra" [
        bench "random matrix" $ nfIO (LA.rand 1000 1000)
        , bench "inv" $ nf LA.inv a
        , bench "Cholesky decomposition" $ nf LA.chol ata
        , bench "vector dot product" $ whnf (LA.dot x) y
        , bench "matrix transpose" $ nf LA.tr a
        , bench "matrix multiplication" $ nf (a LA.<>) b
        , bench "sigmoid" $ nf sigmoid a
        , bench "ReLU" $ nf relu a
        , bench "sum by rows" $ nf (reduceByRows V.sum) a
        , bench "sum by columns" $ nf (reduceByColumns V.sum) a
        , bench "max index in rows" $ nf (reduceByRowsV (fromIntegral . LA.maxIndex) ) a
        , bench "max index in columns" $ nf (reduceByColumnsV (fromIntegral . LA.maxIndex) ) a
        , bench "sigmoid (morpheus)" $ nf M.sigmoid a
        , bench "sigmoid gradient (morpheus)" $ nf M.sigmoidGradient a
        , bench "ReLU (morpheus)" $ nf M.relu a
        , bench "sum by rows (morpheus)" $ nf M.rowSum a
        , bench "sum by columns (morpheus)" $ nf M.columnSum a
        , bench "max index in rows (morpheus)" $ nf M.rowMaxIndex a
        , bench "max index in columns (morpheus)" $ nf M.columnMaxIndex a

        ]
    ]

type Matrix = LA.Matrix LA.R
type Vector = LA.Vector LA.R
type R = LA.R

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

relu :: Matrix -> Matrix
relu x = x * (LA.step x)

reduceByRowsV :: (Vector -> R) -> Matrix -> Vector
reduceByRowsV f = LA.vector . map f . LA.toRows

reduceByColumnsV :: (Vector -> R) -> Matrix -> Vector
reduceByColumnsV f = LA.vector . map f . LA.toColumns

reduceByRows :: (Vector -> R) -> Matrix -> Matrix
reduceByRows f = LA.asColumn . reduceByRowsV f

reduceByColumns :: (Vector -> R) -> Matrix -> Matrix
reduceByColumns f = LA.asRow . reduceByColumnsV f

max_stars_repo_name: Alexander-Ignatyev/linarg-benchmarks

max_issues_repo_name: Alexander-Ignatyev/linarg-benchmarks

max_forks_repo_name: Alexander-Ignatyev/linarg-benchmarks


text:
import           Control.DeepSeq
import           Data.Graph.Inductive
import           Data.Graph.Inductive.PatriciaTree as GP
import           GraphBuilder
import           InterpolationAlgorithms
import           Numeric.LinearAlgebra             (size)
import           Plots
import           System.Environment
import           System.Random                     (randomRIO)
import           Utils

main :: IO ()
main = do
  args <- getArgs
  file <- checkFile args
  let nbds = checkArgs args

  matrix' <- getReqMatrix file

  let (rsize,dims) = size matrix'

{-
The purpose of this algorithm is to find the Riemann normal
coordinates (a coordinate system permissible on Riemann manifolds) of
the tangent space at this point. This is built up by first finding the
K-nearest neighbors to this point, and then applying some linear
dimensionality reduction step to this collection of points (also
called a neighborhood). We perform principal component analysis,
however, other linear dimensionality reduction alogrithms are
permissible. The variable *BASEPOINT* is a randomly chosen point from
the dataset.
-}
  basePoint <- randomRIO (0,rsize-1)

  print $ "The base point is point " ++ show basePoint

{-
Find the M-nearest neighbors about each point (hereby referred to as
the "representative point" with respect to its neighbors) in the
dataset and determine the edges (their necessary existence and weight)
between each member of a neighorhood set and their shared
"representative" point as the Euclidean distance between the
representative point and every point in the neighborhood set.
-}

{-The variable *GCTX*, is a list of contexts for the graph. Essentially,
each context contains the nodes leaving and entering a given node. You
may wish to consult the FGL (Functional Graph Library) documentation
for more information. The variable *POINTS* is a list of the K nearest
neighbors from *BASEPOINT*, a randomnly chosen point in the
dataset. The variable *INDS* is a list of the indices of the points
contained in *POINTS*.
-}
  let (gctx,points,inds) = theContextsAndDists nbds basePoint matrix'

      -- Build the graph from the contexts *GCTX*!
      graph = force . foldr (joinContext . force) (empty :: GP.Gr Double Double) $ gctx

      -- Deterimne the sets of shortest distances from points *POINTS*
      graphShortDists = shortDistMat rsize $ shortestDists inds graph

      {- This matrix is used in the interpolation step for finding the Riemann
          normal coordinates. Please refer to the paper in the repository for
          more information.
          -}
      res = estimateGrad (snd nbds) (dims-1) points
      -- The Riemann normal coordinates.
      res' = riemannNormCoord res graphShortDists

  -- If the graph is connected, then perform and save the results of the analysis.
  if isConnected graph
     then plotAndSave2D' file rsize basePoint nbds res'
     else error "The graph is not connected."

max_stars_repo_name: emmanueldenloye/manifoldRNC

max_issues_repo_name: emmanueldenloye/manifoldRNC

max_forks_repo_name: emmanueldenloye/manifoldRNC


text:
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE OverloadedLists #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE DeriveAnyClass #-}
{-# LANGUAGE DataKinds #-}

module Myelin.SNN where

import Control.Monad.Trans.State
import Control.Monad.Trans.Class
import Control.Lens hiding ((.=), (*~))
import Control.Monad
import Control.Monad.Trans.Identity

import Data.Aeson
import Data.Aeson.Encode.Pretty
import Data.Aeson.TH

import Data.ByteString.Lazy.Char8 as B
import Data.Monoid
import Data.Text
import Data.Traversable

import Numeric.LinearAlgebra

import GHC.Generics

import Myelin.Model
import Myelin.Neuron

-- Builder

-- | The backend targets supported by Myelin
data ExecutionTarget =
    Nest {
        _minTimestep :: Float,
        _maxTimestep :: Float
    }
    | BrainScaleS {
        _wafer :: Int, -- ^ wafer to run the experiment on
        _hicann :: Int -- ^ hicann chip to use
    } -- ^ first generation (wafer) brainscales system
    -- TODO: Add support for other platforms
    -- | Spikey {
    --     _mappingOffset :: Int -- 0..192 (really only 0 and 192 are sensible)
    -- }
    | SpiNNaker -- ^ SpiNNaker neuromorphic platform
    -- | BrainScaleS2 -- ^ second generation brainscales system
    deriving (Eq, Show)

{--
An execution task specifies all information needed to execute a SNN
on a specific target.
--}
data Task = Task {
    _executionTarget :: ExecutionTarget, -- ^ which neuromorphic hardware or simulator to run on
    _network :: Network, -- ^ network that should be implemented
    _simulationTime :: Double -- ^ simulation time in milliseconds
} deriving (Eq, Show)

-- | A spiking neural network described as a list of inputs, nodes, edges and
-- outputs. Also keeps a counter to uniquely label items in the network.
data Network = Network {
    _nextId :: Int,
    _inputs :: [Node],
    _nodes :: [Node],
    _edges :: [Edge],
    _outputs :: [Node]
} deriving (Eq, Show)

-- | Information about the population in the context of data injection/retrieval
data PopulationVisibility = Input | Output | Hidden

makeLenses ''Network

type SNN a m = StateT Network m a

-- | The initial empty network
initialNetwork = Network 0 [] [] [] []

newId :: Monad m => SNN Int m
newId = do
    l <- use nextId
    nextId += 1
    return l

-- | Creates a spike source array from a list of spike times in milliseconds
spikeSourceArray :: Monad m => [Float] -> SNN Node m
spikeSourceArray spikeTimes = do
    id <- newId
    let spikeSource = SpikeSourceArray spikeTimes id
    inputs <>= [spikeSource]
    return spikeSource

-- | Creates a spike source from a poisson distribution
spikeSourcePoisson :: Monad m =>
  Float -> -- ^ The poisson rate
  Integer -> -- ^ The start time in milliseconds
  SNN Node m
spikeSourcePoisson rate start = do
    id <- newId
    let spikeSource = SpikeSourcePoisson rate start id
    inputs <>= [spikeSource]
    return spikeSource

-- | Creates a population of neurons, defined by 'Myelin.Neuron.NeuronType'
population :: Monad m =>
    String -- ^ label of the population (used for printing)
    -> Integer -- ^ size of the population
    -> NeuronType -- ^ type of neuron
    -> PopulationVisibility -- ^ Input, output or hidden population
    -> SNN Node m
population label i typ visibility = do
    l <- newId
    let pop = Population i typ label l
    case visibility of
      Input -> inputs <>= [pop]
      Output -> outputs <>= [pop]
      Hidden -> nodes <>= [pop]
    return pop

-- | Projects two 'Node's together by projecting the first node to the second
-- node, as prescribed by the 'ProjectionEffect'
projection :: Monad m => ProjectionEffect -> [Node] -> [Node] -> SNN () m
projection proj [p0] [p1] = edges <>= [DenseProjection proj p0 p1]
projection proj [p0, p1] [p2] = edges <>= [MergeProjection proj (p0, p1) p2]
projection proj [p0] [p1, p2] = edges <>= [ReplicateProjection proj p0 (p1, p2)]

max_stars_repo_name: volr/myelin

max_issues_repo_name: volr/myelin

max_forks_repo_name: volr/myelin


text:
{-# LANGUAGE FlexibleContexts #-}
-- |
-- Module      : PSO.NeuralSpec
-- Description : Tests for PSO.Neural
-- Copyright   : (c) Tom Westerhout, 2017
-- License     : BSD3
-- Maintainer  : t.westerhout@student.ru.nl
-- Stability   : experimental

module PSO.NeuralSpec where

import           Control.Monad
import           Control.Monad.ST
import           Control.Monad.Reader
import           Data.Complex
import qualified Data.Vector.Storable as V
import           Foreign.Storable
import qualified Numeric.LinearAlgebra as LA
import qualified Numeric.LinearAlgebra.Devel as LA.Devel

import           Test.Hspec

import           PSO.Random
import           PSO.Neural

main :: IO ()
main = hspec spec

roundTo :: (RealFrac a) => Int -> a -> a
roundTo n x = (fromInteger . round $ x * (10^n)) / (10.0^^n)

roundToC :: (RealFrac a) => Int -> Complex a -> Complex a
roundToC n (x :+ y) = (roundTo n x) :+ (roundTo n y)

class AlmostEq a where
  eq' :: a -> a -> Bool

instance AlmostEq Float where
  eq' x y = let epsilon = 1.0E-5
             in abs (x - y) <= epsilon * max (abs x) (abs y)

instance AlmostEq a => AlmostEq (Complex a) where
  eq' (xr :+ xi) (yr :+ yi) = xr `eq'` yr && xi `eq'` yi

instance (Storable a, AlmostEq a) => AlmostEq (V.Vector a) where
  eq' x y = V.and $ V.zipWith eq' x y

toW n m w = LA.Devel.matrixFromVector LA.Devel.ColumnMajor m n w

mkTheta a b w' σ = b + (toW (V.length a) (V.length b) w') LA.#> σ

logWF' :: (Storable a, RealFloat a, Num (V.Vector (Complex a)), LA.Numeric (Complex a))
       => V.Vector (Complex a)
       -> V.Vector (Complex a)
       -> V.Vector (Complex a)
       -> V.Vector (Complex a)
       -> Complex a
logWF' a b w' σ =
  let n = V.length a
      m = V.length b
      θ = mkTheta a b w' σ
   in V.sum (V.zipWith (*) a σ) + V.sum (V.map (log . cosh) θ)

logQuotient1' ::
     (Storable a, RealFloat a, Num (V.Vector (Complex a)), LA.Numeric (Complex a))
  => V.Vector (Complex a)
  -> V.Vector (Complex a)
  -> V.Vector (Complex a)
  -> V.Vector (Complex a)
  -> Int -> Complex a
logQuotient1' a b w' σ flip =
  let σ' = σ V.// [(flip, (-1) * (σ V.! flip))]
   in logWF' a b w' σ' - logWF' a b w' σ

logQuotient2' ::
     (Storable a, RealFloat a, Num (V.Vector (Complex a)), LA.Numeric (Complex a))
  => V.Vector (Complex a)
  -> V.Vector (Complex a)
  -> V.Vector (Complex a)
  -> V.Vector (Complex a)
  -> Int -> Int -> Complex a
logQuotient2' a b w' σ flip1 flip2 =
  let σ' = σ V.// [ (flip1, (-1) * (σ V.! flip1))
                  , (flip2, (-1) * (σ V.! flip2)) ]
   in logWF' a b w' σ' - logWF' a b w' σ

locEnergyHH1DOpen' ::
     (Storable α, RealFloat α, Eq α, Num (V.Vector (Complex α)), LA.Numeric (Complex α))
  => V.Vector (Complex α)
  -> V.Vector (Complex α)
  -> V.Vector (Complex α)
  -> V.Vector (Complex α)
  -> Complex α
locEnergyHH1DOpen' a b w' σ =
  let zipper i x y
        | x == y    = 1
        | otherwise = -1 + 2 * exp (logQuotient2' a b w' σ i (i + 1))
   in V.sum $ V.izipWith zipper σ (V.tail σ)

a1 = V.fromList [   3.7e-01  :+   1.5e-01
                ,   6.5e-01  :+   3.0e-01
                , (-2.8e-01) :+ (-3.0e-01)
                ] :: V.Vector (Complex Float)

ra1 :: UniformDist m (Complex Float) => m (V.Vector (Complex Float))
ra1 = uniformVector 10 ((-1.0) :+ (-1.0), 1.0 :+ 1.0)

b1 = V.fromList [   6.0e-01  :+   8.4e-02
                , (-2.6e-01) :+   9.2e-01
                ,   5.5e-01  :+ (-8.2e-02)
                ,   4.3e-01  :+   3.7e-01
                ,   9.3e-01  :+   5.1e-01
                ,   6.0e-02  :+   8.9e-01
                ] :: V.Vector (Complex Float)

rb1 :: UniformDist m (Complex Float) => m (V.Vector (Complex Float))
rb1 = uniformVector 20 ((-1.0) :+ (-1.0), 1.0 :+ 1.0)

w1 = V.fromList $
      [   0.22 :+ (-0.93),   (-0.16) :+ (-0.47), (-0.57) :+ (-0.83),   (-0.16) :+ 0.86,    (-0.44) :+ 0.16,       0.43 :+ 0.96
      , (-3.6e-2) :+ 0.36, (-1.5e-2) :+ (-0.14),    (-0.15) :+ 0.57,   0.79 :+ (-0.26),     0.95 :+ (-0.6),       0.56 :+ 0.98
      ,      0.82 :+ 0.34,      0.85 :+ (-0.35),    0.64 :+ (-0.24), (-0.73) :+ (-0.8), (-0.34) :+ (-0.24), (-0.94) :+ (-0.46)
      ] :: V.Vector (Complex Float)

rw1 :: UniformDist m (Complex Float) => m (V.Vector (Complex Float))
rw1 = uniformVector (10 * 20) ((-1.0) :+ (-1.0), 1.0 :+ 1.0)

σ1  = V.fromList [ -1, -1, -1 ] :: V.Vector (Complex Float)
σ2  = V.fromList [ 1,   1,  1 ] :: V.Vector (Complex Float)
σ3  = V.fromList [ 1,  -1,  1 ] :: V.Vector (Complex Float)

rσ1 :: Randomisable m Bool => m (V.Vector (Complex Float))
rσ1 = randomSpin 10

rbm1  = runST $ mkRbm a1 b1 w1

mcmc1 = runST $ newMcmc rbm1 σ1 >>= unsafeFreezeMcmc
mcmc2 = runST $ newMcmc rbm1 σ2 >>= unsafeFreezeMcmc
mcmc3 = runST $ newMcmc rbm1 σ3 >>= unsafeFreezeMcmc

spec :: Spec
spec = do
  describe "mkRbm" $ do
    it "Constructs a new RBM given a, b, and w" $
      debugPrintRbm rbm1
  describe "newMcmc" $ do
    it "Constructs a new MCMC given RBM and σ" $
      do debugPrintMcmc mcmc1
         print $ (b1 + (toW 3 6 w1) LA.#> σ1)
  describe "logWF" $ do
    it "1) Calculates ln(ψ(S)) " $
      do putStrLn $ "logWF:  " ++ show (logWF mcmc1)
         putStrLn $ "logWF': " ++ show (logWF' a1 b1 w1 σ1)
    it "2) Calculates ln(ψ(S)) " $
      do putStrLn $ "logWF:  " ++ show (logWF mcmc2)
         putStrLn $ "logWF': " ++ show (logWF' a1 b1 w1 σ2)
    it "3) Calculates ln(ψ(S)) " $
      do putStrLn $ "logWF:  " ++ show (logWF mcmc3)
         putStrLn $ "logWF': " ++ show (logWF' a1 b1 w1 σ3)
  describe "logQuotient1" $ do
    it "1) Calculates ln(ψ(S')/ψ(S))" $
      do putStrLn $ "logQuotient1:  " ++ show (logQuotient1 mcmc1 0)
         putStrLn $ "logQuotient1': " ++ show (logQuotient1' a1 b1 w1 σ1 0)
    it "2) Calculates ln(ψ(S')/ψ(S))" $
      do putStrLn $ "logQuotient1:  " ++ show (logQuotient1 mcmc1 1)
         putStrLn $ "logQuotient1': " ++ show (logQuotient1' a1 b1 w1 σ1 1)
    it "3) Calculates ln(ψ(S')/ψ(S))" $
      do putStrLn $ "logQuotient1:  " ++ show (logQuotient1 mcmc1 2)
         putStrLn $ "logQuotient1': " ++ show (logQuotient1' a1 b1 w1 σ1 2)
    it "4) Calculates ln(ψ(S')/ψ(S))" $
      do putStrLn $ "logQuotient1:  " ++ show (logQuotient1 mcmc2 0)
         putStrLn $ "logQuotient1': " ++ show (logQuotient1' a1 b1 w1 σ2 0)
    it "5) Calculates ln(ψ(S')/ψ(S))" $
      do putStrLn $ "logQuotient1:  " ++ show (logQuotient1 mcmc2 2)
         putStrLn $ "logQuotient1': " ++ show (logQuotient1' a1 b1 w1 σ2 2)
  describe "logQuotient2" $ do
    it "1) Calculates ln(ψ(S')/ψ(S))" $
      do putStrLn $ "logQuotient2:  " ++ show (logQuotient2 mcmc1 0 1)
         putStrLn $ "logQuotient2': " ++ show (logQuotient2' a1 b1 w1 σ1 0 1)
    it "2) Calculates ln(ψ(S')/ψ(S))" $
      do putStrLn $ "logQuotient2:  " ++ show (logQuotient2 mcmc1 0 2)
         putStrLn $ "logQuotient2': " ++ show (logQuotient2' a1 b1 w1 σ1 0 2)
  describe "locEnergyHH1DOpen'" $ do
    it "1) Calculates Eloc(ψ)" $
      do putStrLn $ "locEnergyHH1DOpen:  " ++ show (locEnergyHH1DOpen mcmc3)
         putStrLn $ "locEnergyHH1DOpen': " ++ show (locEnergyHH1DOpen' a1 b1 w1 σ3)
  describe "Random!" $ do
    it "1) Calculates Eloc(ψ)" $
      do let job = do
               a <- ra1
               b <- rb1
               w <- rw1
               σ <- rσ1
               rbm <- mkRbm a b w
               mcmc <- newMcmc rbm σ >>= unsafeFreezeMcmc
               lift $ putStrLn
                    $ "locEnergyHH1DOpen:  " ++ show (locEnergyHH1DOpen mcmc)
               lift $ putStrLn
                    $ "locEnergyHH1DOpen': " ++ show (locEnergyHH1DOpen' a b w σ)
         runReaderT job =<< mkMWCGen (Just 135)
         runReaderT job =<< mkMWCGen (Just 136)
         runReaderT job =<< mkMWCGen (Just 137)

max_stars_repo_name: twesterhout/tcm-swarm

max_issues_repo_name: twesterhout/tcm-swarm

max_forks_repo_name: twesterhout/tcm-swarm


text:
{-# LANGUAGE DataKinds, KindSignatures, TypeOperators #-}
{-# OPTIONS_GHC -Wno-missing-export-lists #-}
-- | Parsing and pre-processing MNIST data.
module HordeAd.Tool.MnistData where

import Prelude

import           Codec.Compression.GZip (decompress)
import           Control.Arrow (first)
import           Data.Array.Internal (valueOf)
import qualified Data.Array.Shaped as OSB
import qualified Data.Array.ShapedS as OS
import qualified Data.ByteString.Lazy as LBS
import           Data.IDX
import           Data.List (sortOn)
import           Data.Maybe (fromMaybe)
import qualified Data.Vector.Generic as V
import qualified Data.Vector.Unboxed
import           GHC.TypeLits (KnownNat, Nat)
import           Numeric.LinearAlgebra (Matrix, Numeric, Vector)
import qualified Numeric.LinearAlgebra as HM
import           System.IO (IOMode (ReadMode), withBinaryFile)
import           System.Random

-- until stylish-haskell accepts NoStarIsType
import qualified GHC.TypeLits

type SizeMnistWidth = 28 :: Nat
type SizeMnistHeight = SizeMnistWidth

sizeMnistWidth :: Int
sizeMnistWidth = valueOf @SizeMnistWidth

type SizeMnistGlyph = SizeMnistWidth GHC.TypeLits.* SizeMnistHeight

sizeMnistGlyph :: Int
sizeMnistGlyph = valueOf @SizeMnistGlyph

type SizeMnistLabel = 10 :: Nat

sizeMnistLabel :: Int
sizeMnistLabel = valueOf @SizeMnistLabel

type LengthTestData = 10000 :: Nat

-- Actually, a better representation, supported by @Data.IDX@,
-- is an integer label and a picture (the same vector as below).
-- Then we'd use @lossCrossEntropy@ that picks a component according
-- to the label instead of performing a dot product with scaling.
-- This results in much smaller Delta expressions.
-- Our library makes this easy to express and gradients compute fine.
-- OTOH, methods with only matrix operations and graphs can't handle that.
-- However, the goal of the exercise it to implement the same
-- neural net that backprop uses for comparative benchmarks.
-- Also, loss computation is not the bottleneck and the more general
-- mechanism that admits non-discrete target labels fuses nicely
-- with softMax. This also seems to be the standard or at least
-- a simple default in tutorial.
type MnistData r = (Vector r, Vector r)

type MnistData2 r = (Matrix r, Vector r)

type MnistDataS r =
  ( OS.Array '[SizeMnistHeight, SizeMnistWidth] r
  , OS.Array '[SizeMnistLabel] r )

type MnistDataBatchS batch_size r =
  ( OS.Array '[batch_size, SizeMnistHeight, SizeMnistWidth] r
  , OS.Array '[batch_size, SizeMnistLabel] r )

shapeBatch :: Numeric r => MnistData r -> MnistDataS r
shapeBatch (input, target) = (OS.fromVector input, OS.fromVector target)

packBatch :: forall batch_size r. (Numeric r, KnownNat batch_size)
          => [MnistDataS r] -> MnistDataBatchS batch_size r
packBatch l =
  let (inputs, targets) = unzip l
  in (OS.ravel $ OSB.fromList inputs, OS.ravel $ OSB.fromList targets)

readMnistData :: LBS.ByteString -> LBS.ByteString -> [MnistData Double]
readMnistData glyphsBS labelsBS =
  let glyphs = fromMaybe (error "wrong MNIST glyphs file")
               $ decodeIDX glyphsBS
      labels = fromMaybe (error "wrong MNIST labels file")
               $ decodeIDXLabels labelsBS
      intData = fromMaybe (error "can't decode MNIST file into integers")
                $ labeledIntData labels glyphs
      f :: (Int, Data.Vector.Unboxed.Vector Int) -> MnistData Double
      -- Copied from library backprop to enable comparison of results.
      -- I have no idea how this is different from @labeledDoubleData@, etc.
      f (labN, v) =
        ( V.map (\r -> fromIntegral r / 255) $ V.convert v
        , V.generate sizeMnistLabel (\i -> if i == labN then 1 else 0) )
  in map f intData

trainGlyphsPath, trainLabelsPath, testGlyphsPath, testLabelsPath :: FilePath
trainGlyphsPath = "samplesData/train-images-idx3-ubyte.gz"
trainLabelsPath = "samplesData/train-labels-idx1-ubyte.gz"
testGlyphsPath  = "samplesData/t10k-images-idx3-ubyte.gz"
testLabelsPath  = "samplesData/t10k-labels-idx1-ubyte.gz"

loadMnistData :: FilePath -> FilePath -> IO [MnistData Double]
loadMnistData glyphsPath labelsPath =
  withBinaryFile glyphsPath ReadMode $ \glyphsHandle ->
    withBinaryFile labelsPath ReadMode $ \labelsHandle -> do
      glyphsContents <- LBS.hGetContents glyphsHandle
      labelsContents <- LBS.hGetContents labelsHandle
      return $! readMnistData (decompress glyphsContents)
                              (decompress labelsContents)

loadMnistData2 :: FilePath -> FilePath -> IO [MnistData2 Double]
loadMnistData2 glyphsPath labelsPath = do
  ds <- loadMnistData glyphsPath labelsPath
  return $! map (first $ HM.reshape sizeMnistWidth) ds

-- Good enough for QuickCheck, so good enough for me.
shuffle :: RandomGen g => g -> [a] -> [a]
shuffle g l =
  let rnds = randoms g :: [Int]
  in map fst $ sortOn snd $ zip l rnds

chunksOf :: Int -> [e] -> [[e]]
chunksOf n = go where
  go [] = []
  go l = let (chunk, rest) = splitAt n l
         in chunk : go rest

max_stars_repo_name: Mikolaj/horde-ad

max_issues_repo_name: Mikolaj/horde-ad

max_forks_repo_name: Mikolaj/horde-ad


text:
module Distributions where

import Control.Monad (replicateM)
import qualified Control.Monad.Bayes.Class as B
import Control.Monad.Bayes.Class (MonadSample, MonadCond)
import Numeric.SpecFunctions (logGamma)

import qualified Data.Vector as V
import qualified Data.Map as M
import qualified Data.MultiSet as S

import Numeric.Log (Log (Exp))

import SymbolicArithmetic
import Util.Numeric (logFact)
import Util.Bijection

type Distr = Distr' Double

data Distr' d a = Distr
  { sample :: forall m. MonadSample m => m a
  , score :: a -> d
  }

mapOutputDiscrete :: Bijection a b -> Distr' d a -> Distr' d b
mapOutputDiscrete (Bijection to from) (Distr sample score) =
  Distr (to <$> sample) (score . from)

factor :: MonadCond m => Double -> m ()
factor = B.score . Numeric.Log.Exp

negInf :: Double
negInf = - 1 / 0

observe :: MonadCond m => Distr a -> a -> m ()
observe d x = factor (score d x)

data ExpFam a = ExpFam
  { dim :: Int
  , suffStats :: [Exp () Double]
  , logPartition :: [Double] -> Double
  , logBaseMeas :: Exp () Double
  , efSample :: forall m. MonadSample m => [Double] -> m a
  , efRepr :: a -> Double
  , mle :: [Double] -> [Double]
  , efDescribe :: [Double] -> String }

bernoulliEF :: ExpFam Bool
bernoulliEF = ExpFam
  { dim = 1
  , suffStats = [Var ()]
  , logPartition = undefined
  , logBaseMeas = undefined
  , efSample = \[p] -> B.bernoulli p
  , efRepr = \x -> if x then 1 else 0
  , mle = id
  , efDescribe = \[eta] -> let z = exp eta in "Bernoulli(" ++ show (z / (1 + z)) ++ ")"
  }

betaEF :: ExpFam Double
betaEF = ExpFam
  { dim = 2
  , suffStats = let x = Var () in  [log x, log (1 - x)]
  , logPartition = \[a, b] -> logBeta a b
  , logBaseMeas = let x = Var () in - log x - log (1 - x)
  , efSample = \[a, b] -> B.beta a b
  , efRepr = id
  , mle = \[lx, lox] -> undefined
  , efDescribe = \[a, b] -> "Beta(" ++ show a ++ ", " ++ show b ++ ")"
  }

logBeta :: Double -> Double -> Double
logBeta a b = logGamma a + logGamma b - logGamma (a + b)

logChoose :: Int -> Int -> Double
logChoose n k = logFact n - logFact k - logFact (n - k)

normalEF :: ExpFam Double
normalEF = ExpFam
  { dim = 2
  , suffStats = let x = Var () in [x, x^2]
  , logPartition = \[eta1, eta2] -> - eta1^2 / (4 * eta2) - 1 / 2 * log(-2 * eta2)
  , logBaseMeas = - 1 / 2 * log (2 * pi)
  , efSample = \eta -> let (mu, sigma2) = conv eta in
      B.normal mu (sqrt sigma2)
  , efRepr = id
  , mle = \[mu, ex2] -> let sigma2 = ex2 - mu^2 in normalNatParams mu sigma2
  , efDescribe = \eta -> let (mu, sigma2) = conv eta in
    "Normal(" ++ show mu ++ ", " ++ show sigma2 ++ ")"
  }
  where
  conv [eta1, eta2] = let sigma2 = - 1 / (2 * eta2) in (sigma2 * eta1, sigma2)
  conv _ = error "conv"

-- indep :: ExpFam a -> ExpFam b -> ExpFam (a, b)
-- indep efa efb = ExpFam
--   { dim = dim efa + dim efb
--   , suffStats = \(a, b) -> suffStats efa a ++ suffStats efb b
--   , logPartition = \eta -> let (eta1, eta2) = split eta in
--       logPartition efa eta1 + logPartition efb eta2
--   , logBaseMeas = \(a, b) -> logBaseMeas efa a + logBaseMeas efb b
--   , efSample = \eta -> let (eta1, eta2) = split eta in
--       (,) <$> efSample efa eta1 <*> efSample efb eta2
--   , mle = \ss -> let (ss1, ss2) = split ss in
--       mle efa ss1 ++ mle efb ss2
--   , efDescribe = \eta -> let (eta1, eta2) = split eta in
--     "Indep(" ++ efDescribe efa eta1 ++ ", " ++ efDescribe efb eta2 ++ ")"
--   }
--   where
--   split = splitAt (dim efa)

normalNatParams :: Fractional d => d -> d -> [d]
normalNatParams mu sigma2 = [mu / sigma2, -1 / (2 * sigma2)]

dirac :: Eq a => a -> Distr a
dirac x = Distr (return x) (\y -> if x == y then 0 else negInf)

dot :: Num d => [d] -> [d] -> d
dot xs ys = sum (zipWith (*) xs ys)

efToDistr :: ExpFam a -> [Double] -> Distr a
efToDistr ef naturalParams =
  Distr (efSample ef naturalParams)
  (\x -> let ev = eval (\_ -> efRepr ef x) in
    naturalParams `dot` map ev (suffStats ef) - logPartition ef naturalParams + ev (logBaseMeas ef))

bernoulli :: Double -> Distr Bool
bernoulli p = Distr (B.bernoulli p) (\b -> bernoulli_ll p (if b then 1 else 0))

bernoulli01 :: Double -> Distr Int
bernoulli01 p = Distr sam sco where
  sam :: MonadSample m => m Int
  sam = do
    b <- B.bernoulli p
    return (if b then 1 else 0)
  sco 1 = log p
  sco 0 = log (1 - p)
  sco _ = negInf

categorical :: [Double] -> Distr Int
categorical ps = Distr (B.categorical (V.fromList ps)) (\i -> log (ps !! i))

categoricalM :: Ord a => M.Map a Double -> Distr a
categoricalM ps = Distr (fmap (\i -> fst (M.elemAt i ps)) (B.categorical (V.fromList (M.elems ps))))
  (\k -> log (ps M.! k))

bernoulli_ll :: Floating a => a -> a -> a
bernoulli_ll p b = b * log p + (1 - b) * log (1 - p)

-- WARNING: NOT NORMALIZED!
beta_ll :: Floating a => a -> a -> a -> a
beta_ll a b p = a * log p + b * log (1 - p)

beta :: Double -> Double -> Distr Double
beta a b = Distr (B.beta a b) (\x -> beta_ll a b x + logPartition betaEF [a, b])

gaussian_ll :: Floating a => a -> a -> a -> a
gaussian_ll mu sigma2 x = - x^2 / sigma2 + 2 * x * mu / sigma2 - mu^2 / sigma2 - log (2 * pi) - log sigma2

gaussian_ll' :: Floating a => a -> a -> a -> a
gaussian_ll' mu sigma2 x = - (x - mu)^2 / sigma2 - log (2 * pi) - log sigma2

gaussian :: Double -> Double -> Distr Double
gaussian mu sigma2 = Distr (B.normal mu (sqrt sigma2)) (gaussian_ll' mu sigma2)

normal = gaussian

poisson_ll :: Double -> Int -> Double
poisson_ll lambda k = fromIntegral k * log lambda - lambda - logGamma (fromIntegral k + 1)

poisson :: Double -> Distr Int
poisson lambda = Distr (B.poisson lambda) (poisson_ll lambda)

uniform :: Double -> Double -> Distr Double
uniform a b = Distr (B.uniform a b) $ \x ->
  if a <= x && x <= b
    then - log range
    else negInf
  where
  range = b - a

replicateNIID :: Int -> Distr a -> Distr [a]
replicateNIID n d = Distr (replicateM n (sample d)) (\obs ->
  if length obs == n
    then sum (map (score d) obs)
    else (-1 / 0))

uniformIntRange :: Int -> Distr Int
uniformIntRange max = Distr (B.uniformD [0 .. max - 1])
  (\obs -> if 0 <= obs && obs < max then - log (fromIntegral max) else -1 / 0)

bind :: Distr a -> (a -> Distr b) -> Distr (a, b)
bind d f = Distr (do {x <- sample d; y <- sample (f x); pure (x, y)})
                 (\(x, y) -> score d x + score (f x) y)

replicateIID :: forall a. Distr Int -> Distr a -> Distr [a]
replicateIID howMany d = Distr sam obs where
  sam :: MonadSample m => m [a]
  sam = do
    n <- sample howMany
    replicateM n (sample d)
  obs xs = score howMany (length xs) + sum (map (score d) xs)

geometric' :: MonadSample m => m Int
geometric' = do
  x <- B.bernoulli (0.5 :: Double)
  if x then return 0 else fmap (+1) geometric'

exponential :: MonadSample m => Double -> m Double
exponential lambda = do
  u <- B.random
  pure (- log u / lambda)

passert :: MonadCond m => Bool -> m ()
passert True = pure ()
passert False = B.score 0

data Some f where
  Some :: f a -> Some f

exponentialFamilies :: [Some ExpFam]
exponentialFamilies = [Some betaEF, Some bernoulliEF, Some normalEF]

shuffleWithRepeats' :: forall a. Ord a => M.Map a (Distr Int) -> Distr [a]
shuffleWithRepeats' countDistrs = Distr sam sco where
  sam :: MonadSample m => m [a]
  sam = do
    counts <- mapM sample countDistrs
    randomlyInterleave [ replicate n k | (k, n) <- M.toList counts ]
  sco :: [a] -> Double
  sco idxes = let ns = mconcat [ S.singleton i | i <- idxes ] in
    sum (M.mapWithKey (\k d -> score d (length (filter (== k) idxes))) countDistrs)
    + randomlyInterleaveLogPDF (map snd (S.toOccurList ns))

shuffleWithRepeats :: [Distr Int] -> Distr [Int]
shuffleWithRepeats countDistrs = Distr sam sco where
  sam :: MonadSample m => m [Int]
  sam = do
    counts <- mapM sample countDistrs
    randomlyInterleave [ replicate n i | (i, n) <- zip [0..] counts ]
  sco :: [Int] -> Double
  sco idxes = let ns = [ length (filter (== i) idxes) | i <- take (length countDistrs) [0..] ] in
    sum (zipWith score countDistrs ns) + randomlyInterleaveLogPDF ns

randomlyInterleave :: MonadSample m => [[a]] -> m [a]
randomlyInterleave xs = if ntot == 0
  then pure []
  else do
    i <- sample (categorical [ fromIntegral n / fromIntegral ntot | n <- ns ])
    let (before, (v : vs) : after) = splitAt i xs
    (v :) <$> randomlyInterleave (before ++ vs : after)
  where
    ns = map length xs
    ntot = sum ns

logNumPermutationsWithRepeats :: [Int] -> Double
logNumPermutationsWithRepeats ns = logFact (sum ns) - sum [ logFact n | n <- ns ]

randomlyInterleaveLogPDF :: [Int] -> Double
randomlyInterleaveLogPDF ns = - logNumPermutationsWithRepeats ns

shuffleListLogPDF :: Int -> Double
shuffleListLogPDF n = - logFact n

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 PartialTypeSignatures #-}
module PMMHTest where

import           Control.Monad.Bayes.Class
import           Control.Monad.Bayes.Sampler
import           Control.Monad
import           Control.Monad.Bayes.Inference.PMMH
                                               as PMMH
import           Control.Monad.Bayes.Weighted
import           Control.Monad.State
import           Numeric.Log
import           DataParser
import           Statistics
import           Utils

generateData
  :: MonadSample m
  =>
  -- | T
     Int
  ->
  -- | list of latent and observable states from t=1
     m [(Double, Double)]
generateData t = do
  (sigmaX, sigmaY) <- param
  let sq x = x * x
      simulate 0 _ acc = return acc
      simulate k x acc = do
        let n = length acc
        x' <- normal (mean x n) sigmaX
        y' <- normal (sq x' / 20) sigmaY
        simulate (k - 1) x' ((x', y') : acc)
  x0  <- normal 0 (sqrt 5)
  xys <- simulate t x0 []
  return $ reverse xys

param :: MonadSample m => m (Double, Double)
param = do
  let a = 0.01
  let b = 0.01
  precX <- gamma a b
  let sigmaX = 1 / sqrt precX
  precY <- gamma a b
  let sigmaY = 1 / sqrt precY
  return (sigmaX, sigmaY)

mean :: Double -> Int -> Double
mean x n =
  let sq x = x * x
  in  0.5 * x + 25 * x / (1 + sq x) + 8 * cos (1.2 * fromIntegral n)

model
  :: (MonadInfer m)
  =>
  -- | observed data
     [Double]
  ->
  -- | prior on the parameters
     (Double, Double)
  ->
  -- | list of latent states from t=1
     m [Double]
model obs (sigmaX, sigmaY) = do
  let sq x = x * x
      simulate []       _ acc = return acc
      simulate (y : ys) x acc = do
        let n = length acc
        x' <- normal (mean x n) sigmaX
        factor $ normalPdf (sq x' / 20) sigmaY y
        simulate ys x' (x' : acc)
  x0 <- normal 0 (sqrt 5)
  xs <- simulate obs x0 []
  return $ reverse xs

test :: IO [[([Double], Numeric.Log.Log Double)]]
test = sampleIO $ do
  let t = 5
  dat <- generateData t
  let ys = map snd dat
  pmmhRes <- prior $ pmmh 2 t 3 param (model ys)
  liftIO $ print pmmhRes
  return pmmhRes

max_stars_repo_name: rossng/sir-monad

max_issues_repo_name: rossng/sir-monad

max_forks_repo_name: rossng/sir-monad


text:
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE FlexibleContexts #-}
#if __GLASGOW_HASKELL__ >= 800
  {-# OPTIONS_GHC -Wno-redundant-constraints #-}
#endif
-- |
-- Module      : Graphics.Image.ColorSpace.Complex
-- Copyright   : (c) Alexey Kuleshevich 2016
-- License     : BSD3
-- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
-- Stability   : experimental
-- Portability : non-portable
--
module Graphics.Image.ColorSpace.Complex (
  -- *** Rectangular form
  Complex(..), (+:), realPart, imagPart,
  -- *** Polar form
  mkPolar, cis, polar, magnitude, phase,
  -- *** Conjugate
  conjugate
  ) where

import Graphics.Image.Interface (Pixel)
import Control.Applicative
import Data.Complex (Complex(..))
import qualified Data.Complex as C hiding (Complex(..))

infix 6 +:

-- | Constrcut a complex pixel from two pixels representing real and imaginary parts.
--
-- @ PixelRGB 4 8 6 '+:' PixelRGB 7 1 1 __==__ PixelRGB (4 ':+' 7) (8 ':+' 1) (6 ':+' 1) @
--
(+:) :: Applicative (Pixel cs) => Pixel cs e -> Pixel cs e -> Pixel cs (Complex e)
(+:) = liftA2 (:+)
{-# INLINE (+:) #-}

-- | Extracts the real part of a complex pixel.
realPart :: (Applicative (Pixel cs), RealFloat e) => Pixel cs (Complex e) -> Pixel cs e
realPart = liftA C.realPart
{-# INLINE realPart #-}

-- | Extracts the imaginary part of a complex pixel.
imagPart :: (Applicative (Pixel cs), RealFloat e) => Pixel cs (Complex e) -> Pixel cs e
imagPart = liftA C.imagPart
{-# INLINE imagPart #-}

-- | Form a complex pixel from polar components of magnitude and phase.
mkPolar :: (Applicative (Pixel cs), RealFloat e) =>
           Pixel cs e -> Pixel cs e -> Pixel cs (Complex e)
mkPolar = liftA2 C.mkPolar
{-# INLINE mkPolar #-}

-- | @'cis' t@ is a complex pixel with magnitude 1 and phase t (modulo @2*'pi'@).
cis :: (Applicative (Pixel cs), RealFloat e) => Pixel cs e -> Pixel cs (Complex e)
cis = liftA C.cis
{-# INLINE cis #-}

-- | The function @'polar'@ takes a complex pixel and returns a (magnitude, phase)
-- pair of pixels in canonical form: the magnitude is nonnegative, and the phase
-- in the range @(-'pi', 'pi']@; if the magnitude is zero, then so is the phase.
polar :: (Applicative (Pixel cs), RealFloat e) => Pixel cs (Complex e) -> (Pixel cs e, Pixel cs e)
polar !zPx = (magnitude zPx, phase zPx)
{-# INLINE polar #-}

-- | The nonnegative magnitude of a complex pixel.
magnitude :: (Applicative (Pixel cs), RealFloat e) => Pixel cs (Complex e) -> Pixel cs e
magnitude = liftA C.magnitude
{-# INLINE magnitude #-}

-- | The phase of a complex pixel, in the range @(-'pi', 'pi']@. If the
-- magnitude is zero, then so is the phase.
phase :: (Applicative (Pixel cs), RealFloat e) => Pixel cs (Complex e) -> Pixel cs e
phase = liftA C.phase
{-# INLINE phase #-}

-- | The conjugate of a complex pixel.
conjugate :: (Applicative (Pixel cs), RealFloat e) => Pixel cs (Complex e) -> Pixel cs (Complex e)
conjugate = liftA C.conjugate
{-# INLINE conjugate #-}

max_stars_repo_name: khilanravani/hip

max_issues_repo_name: khilanravani/hip

max_forks_repo_name: khilanravani/hip


text:
module MachineLearning.Classification.BinaryTest
(
  tests
  , testOptPath
)

where

import Test.Framework (testGroup)
import Test.Framework.Providers.HUnit
import Test.HUnit
import Test.HUnit.Approx
import Test.HUnit.Plus

import MachineLearning.Types
import qualified Data.Vector as V
import qualified Numeric.LinearAlgebra as LA

import MachineLearning.DataSets (dataset2)

import qualified MachineLearning as ML
import MachineLearning.Classification.Binary

(x, y) = ML.splitToXY dataset2

processX muSigma x = ML.addBiasDimension $ ML.featureNormalization muSigma $ ML.mapFeatures 6 x

muSigma = ML.meanStddev (ML.mapFeatures 6 x)
x1 = processX muSigma x
zeroTheta = LA.konst 0 (LA.cols x1)

xPredict = LA.matrix 2 [ -0.5, 0.5
                       , 0.2, -0.2
                       , 1, 1
                       , 1, 0
                       , 0, 0
                       , 0, 1]
xPredict1 = processX muSigma xPredict
yExpected = LA.vector [1, 1, 0, 0, 1, 0]

eps = 0.0001

-- Binary

(thetaCGFR, optPathCGFR) = learn (ConjugateGradientFR 0.1 0.1) eps 50 (L2 0.5) x1 y zeroTheta
(thetaCGPR, optPathCGPR) = learn (ConjugateGradientPR 0.1 0.1) eps 50 (L2 0.5) x1 y zeroTheta
(thetaBFGS, optPathBFGS) = learn (BFGS2 0.1 0.1) eps 50 (L2 0.5) x1 y zeroTheta

isInDescendingOrder :: V.Vector Double -> Bool
isInDescendingOrder lst = V.and . snd . V.unzip $ V.scanl (\(prev, _) current -> (current, prev-current > (-0.001))) (1/0, True) lst

testOptPath optPath = do
  let js = V.convert $ (LA.toColumns optPath) !! 1
  assertBool ("non-increasing errors: " ++ show js) $ isInDescendingOrder js

testAccuracyBinary theta eps = do
  let yPredicted = predict x1 theta
      accuracy = calcAccuracy y yPredicted
  assertApproxEqual "" eps 1 accuracy

tests = [
  testGroup "learn" [
      testCase "Conjugate Gradient FR" $ assertVector "" 0.01 yExpected (predict xPredict1 thetaCGFR)
      , testCase "Conjugate Gradient PR" $ assertVector "" 0.01 yExpected (predict xPredict1 thetaCGPR)
      , testCase "BFGS" $ assertVector "" 0.01 yExpected (predict xPredict1 thetaBFGS)
      ]
  , testGroup "optPath" [
      testCase "Conjugate Gradient FR" $ testOptPath optPathCGFR
      , testCase "Conjugate Gradient PR" $ testOptPath optPathCGPR
      , testCase "BFGS" $ testOptPath optPathBFGS
      ]
    , testGroup "accuracy" [
        testCase "Conjugate Gradient FR" $ testAccuracyBinary thetaCGFR 0.2
        , testCase "Conjugate Gradient PR" $ testAccuracyBinary thetaCGPR 0.2
        , testCase "BFGS" $ testAccuracyBinary thetaBFGS 0.2
        ]
  ]

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool