text:
-- monadic computations
-- (contributed by Vivian McPhail)

import Numeric.LinearAlgebra
import Control.Monad.State.Strict
import Control.Monad.Maybe
import Foreign.Storable(Storable)
import System.Random(randomIO)

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

-- an instance of MonadIO, a monad transformer
type VectorMonadT = StateT Int IO

test1 :: Vector Int -> IO (Vector Int)
test1 = mapVectorM $ \x -> do
    putStr $ (show x) ++ " "
    return (x + 1)

-- we can have an arbitrary monad AND do IO
addInitialM :: Vector Int -> VectorMonadT ()
addInitialM = mapVectorM_ $ \x -> do
    i <- get
    liftIO $ putStr $ (show $ x + i) ++ " "
    put $ x + i

-- sum the values of the even indiced elements
sumEvens :: Vector Int -> Int
sumEvens = foldVectorWithIndex (\x a b -> if x `mod` 2 == 0 then a + b else b) 0

-- sum and print running total of evens
sumEvensAndPrint :: Vector Int -> VectorMonadT ()
sumEvensAndPrint = mapVectorWithIndexM_ $ \ i x -> do
    when (i `mod` 2 == 0) $ do
        v <- get
        put $ v + x
        v' <- get
        liftIO $ putStr $ (show v') ++ " "

indexPlusSum :: Vector Int -> VectorMonadT ()
indexPlusSum v' = do
    let f i x = do
            s <- get
            let inc = x+s
            liftIO $ putStr $ show (i,inc) ++ " "
            put inc
            return inc
    v <- mapVectorWithIndexM f v'
    liftIO $ do
        putStrLn ""
        putStrLn $ show v

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

-- short circuit
monoStep :: Double -> MaybeT (State Double) ()
monoStep d = do
    dp <- get
    when (d < dp) (fail "negative difference")
    put d
{-# INLINE monoStep #-}

isMonotoneIncreasing :: Vector Double -> Bool
isMonotoneIncreasing v =
    let res = evalState (runMaybeT $ (mapVectorM_ monoStep v)) (v @> 0)
     in case res of
        Nothing -> False
        Just _  -> True

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

-- | apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs
successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool
successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ step (subVector 1 (dim v - 1) v))) (v @> 0)
   where step e = do
                  ep <- lift $ get
                  if t e ep
                     then lift $ put e
                     else (fail "successive_ test failed")

-- | operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input
successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b
successive f v = evalState (mapVectorM step (subVector 1 (dim v - 1) v)) (v @> 0)
   where step e = do
                  ep <- get
                  put e
                  return $ f ep e

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

v :: Vector Int
v = 10 |> [0..]

w = fromList ([1..10]++[10,9..1]) :: Vector Double

main = do
    v' <- test1 v
    putStrLn ""
    putStrLn $ show v'
    evalStateT (addInitialM v) 0
    putStrLn ""
    putStrLn $ show (sumEvens v)
    evalStateT (sumEvensAndPrint v) 0
    putStrLn ""
    evalStateT (indexPlusSum v) 0
    putStrLn "-----------------------"
    mapVectorM_ print v
    print =<< (mapVectorM (const randomIO) v :: IO (Vector Double))
    print =<< (mapVectorM (\a -> fmap (+a) randomIO) (5|>[0,100..1000]) :: IO (Vector Double))
    putStrLn "-----------------------"
    print $ isMonotoneIncreasing w
    print $ isMonotoneIncreasing (subVector 0 7 w)
    print $ successive_ (>) v
    print $ successive_ (>) w
    print $ successive (+) v

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
{-# LANGUAGE CPP                   #-}
{-# LANGUAGE ConstraintKinds       #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RebindableSyntax      #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeApplications      #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeSynonymInstances  #-}
{-# LANGUAGE UndecidableInstances  #-}
{-# LANGUAGE ViewPatterns          #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}
-- |
-- Module      : Data.Array.Accelerate.Data.Complex
-- Copyright   : [2015..2017] Trevor L. McDonell
-- License     : BSD3
--
-- Maintainer  : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>
-- Stability   : experimental
-- Portability : non-portable (GHC extensions)
--
-- Complex numbers, stored in the usual C-style array-of-struct representation,
-- for easy interoperability.
--
module Data.Array.Accelerate.Data.Complex (

  -- * Rectangular from
  Complex(..),
  real,
  imag,

  -- * Polar form
  mkPolar,
  cis,
  polar,
  magnitude,
  phase,

  -- * Conjugate
  conjugate,

) where

import Data.Array.Accelerate.Array.Sugar
import Data.Array.Accelerate.Classes
import Data.Array.Accelerate.Data.Functor
import Data.Array.Accelerate.Prelude
import Data.Array.Accelerate.Product
import Data.Array.Accelerate.Smart
import Data.Array.Accelerate.Type

import Prelude                                                      ( ($) )
import Data.Complex                                                 ( Complex(..) )
import qualified Data.Complex                                       as C
import qualified Prelude                                            as P

-- Use an array-of-structs representation for complex numbers. This matches the
-- standard C-style layout, but means that we can define instances only at
-- specific types (not for any type 'a') as we can only have vectors of
-- primitive type.
--

instance Elt (Complex Half) where
  type EltRepr (Complex Half) = V2 Half
  {-# INLINE eltType     #-}
  {-# INLINE [1] toElt   #-}
  {-# INLINE [1] fromElt #-}
  eltType          = TypeRscalar scalarType
  toElt (V2 r i)   = r :+ i
  fromElt (r :+ i) = V2 r i

instance Elt (Complex Float) where
  type EltRepr (Complex Float) = V2 Float
  {-# INLINE eltType     #-}
  {-# INLINE [1] toElt   #-}
  {-# INLINE [1] fromElt #-}
  eltType          = TypeRscalar scalarType
  toElt (V2 r i)   = r :+ i
  fromElt (r :+ i) = V2 r i

instance Elt (Complex Double) where
  type EltRepr (Complex Double) = V2 Double
  {-# INLINE eltType     #-}
  {-# INLINE [1] toElt   #-}
  {-# INLINE [1] fromElt #-}
  eltType          = TypeRscalar scalarType
  toElt (V2 r i)   = r :+ i
  fromElt (r :+ i) = V2 r i

instance Elt (Complex CFloat) where
  type EltRepr (Complex CFloat) = V2 Float
  {-# INLINE eltType     #-}
  {-# INLINE [1] toElt   #-}
  {-# INLINE [1] fromElt #-}
  eltType                        = TypeRscalar scalarType
  toElt (V2 r i)                 = CFloat r :+ CFloat i
  fromElt (CFloat r :+ CFloat i) = V2 r i

instance Elt (Complex CDouble) where
  type EltRepr (Complex CDouble) = V2 Double
  {-# INLINE eltType     #-}
  {-# INLINE [1] toElt   #-}
  {-# INLINE [1] fromElt #-}
  eltType                          = TypeRscalar scalarType
  toElt (V2 r i)                   = CDouble r :+ CDouble i
  fromElt (CDouble r :+ CDouble i) = V2 r i

instance cst a => IsProduct cst (Complex a) where
  type ProdRepr (Complex a) = ProdRepr (a,a)
  fromProd (r :+ i) = fromProd @cst (r,i)
  toProd p          = let (r,i) = toProd @cst p in (r :+ i)
  prod              = prod @cst @(a,a)

instance (Lift Exp a, Elt (Plain a), Elt (Complex (Plain a))) => Lift Exp (Complex a) where
  type Plain (Complex a) = Complex (Plain a)
  lift (r :+ i)          = Exp $ Tuple (NilTup `SnocTup` lift r `SnocTup` lift i)

instance (Elt a, Elt (Complex a)) => Unlift Exp (Complex (Exp a)) where
  unlift e
    = let r     = Exp $ SuccTupIdx ZeroTupIdx `Prj` e
          i     = Exp $ ZeroTupIdx `Prj` e
      in
      r :+ i

instance (Eq a, Elt (Complex a)) => Eq (Complex a) where
  x == y = let r1 :+ c1 = unlift x
               r2 :+ c2 = unlift y
           in  r1 == r2 && c1 == c2
  x /= y = let r1 :+ c1 = unlift x
               r2 :+ c2 = unlift y
           in  r1 /= r2 || c1 /= c2

instance (RealFloat a, Elt (Complex a)) => P.Num (Exp (Complex a)) where
  (+)           = lift2 ((+) :: Complex (Exp a) -> Complex (Exp a) -> Complex (Exp a))
  (-)           = lift2 ((-) :: Complex (Exp a) -> Complex (Exp a) -> Complex (Exp a))
  (*)           = lift2 ((*) :: Complex (Exp a) -> Complex (Exp a) -> Complex (Exp a))
  negate        = lift1 (negate :: Complex (Exp a) -> Complex (Exp a))
  signum z      = if z == 0
                    then z
                    else let x :+ y = unlift z
                             r      = magnitude z
                         in
                         lift (x/r :+ y/r)
  abs z         = lift (magnitude z :+ 0)
  fromInteger n = lift (fromInteger n :+ 0)

instance (RealFloat a, Elt (Complex a)) => P.Fractional (Exp (Complex a)) where
  fromRational x  = lift (fromRational x :+ 0)
  z / z'          = lift ((x*x''+y*y'') / d :+ (y*x''-x*y'') / d)
    where
      x  :+ y   = unlift z
      x' :+ y'  = unlift z'
      --
      x'' = scaleFloat k x'
      y'' = scaleFloat k y'
      k   = - max (exponent x') (exponent y')
      d   = x'*x'' + y'*y''

instance (RealFloat a, Elt (Complex a)) => P.Floating (Exp (Complex a)) where
  pi                      = lift $ pi :+ 0

  exp (unlift -> x :+ y)  = let expx = exp x
                            in  complex $ expx * cos y :+ expx * sin y

  log z                   = lift $ log (magnitude z) :+ phase z

  sqrt z@(unlift -> x :+ y) =
    if z == 0
      then 0
      else lift $ u :+ (y < 0 ? (-v, v))
    where
      (u,v) = unlift (x < 0 ? (lift (v',u'), lift (u',v')))
      v'    = abs y / (u'*2)
      u'    = sqrt ((magnitude z + abs x) / 2)

  x ** y =
    if y == 0 then 1 else
    if x == 0 then if exp_r > 0 then 0 else
                   if exp_r < 0 then lift (inf :+ 0)
                                else lift (nan :+ nan)
              else if isInfinite r || isInfinite i
                     then if exp_r > 0 then lift (inf :+ 0) else
                          if exp_r < 0 then 0
                                       else lift (nan :+ nan)
                     else exp (log x * y)
    where
      r     :+ i  = unlift x
      exp_r :+ _  = unlift y
      --
      inf = 1 / 0
      nan = 0 / 0

  sin (unlift -> x :+ y)  = complex $ sin x * cosh y :+ cos x * sinh y
  cos (unlift -> x :+ y)  = complex $ cos x * cosh y :+ (- sin x * sinh y)
  tan (unlift -> x :+ y)  = (complex $ sinx*coshy :+ cosx*sinhy) / (complex $ cosx*coshy :+ (-sinx*sinhy))
    where
      sinx  = sin x
      cosx  = cos x
      sinhy = sinh y
      coshy = cosh y

  sinh (unlift -> x :+ y) = complex $ cos y * sinh x :+ sin  y * cosh x
  cosh (unlift -> x :+ y) = complex $ cos y * cosh x :+ sin y * sinh x
  tanh (unlift -> x :+ y) = (complex $ cosy*sinhx :+ siny*coshx) / (complex $ cosy*coshx :+ siny*sinhx)
    where
      siny  = sin y
      cosy  = cos y
      sinhx = sinh x
      coshx = cosh x

  asin z@(unlift -> x :+ y) = complex $ y' :+ (-x')
    where
      x' :+ y' = unlift $ log ((complex ((-y):+x)) + sqrt (1 - z*z))

  acos z                    = complex $ y'' :+ (-x'')
    where
      x'' :+ y''  = unlift $ log (z + (complex ((-y') :+ x')))
      x'  :+ y'   = unlift $ sqrt (1 - z*z)

  atan z@(unlift -> x :+ y) = complex $ y' :+ (-x')
    where
      x' :+ y' = unlift $ log ((complex ((1-y):+x)) / sqrt (1+z*z))

  asinh z =  log (z + sqrt (1+z*z))
  acosh z =  log (z + (z+1) * sqrt ((z-1)/(z+1)))
  atanh z =  0.5 * log ((1.0+z) / (1.0-z))

instance (FromIntegral a b, Num b, Elt (Complex b)) => FromIntegral a (Complex b) where
  fromIntegral x = lift (fromIntegral x :+ 0)

-- | @since 1.2.0.0
instance Functor Complex where
  fmap f (unlift -> r :+ i) = lift (f r :+ f i)

-- Helper function to fix the types for lift (ugh)
--
complex :: (Elt a, Elt (Complex a)) => Complex (Exp a) -> Exp (Complex a)
complex = lift

-- | The non-negative magnitude of a complex number
--
magnitude :: (RealFloat a, Elt (Complex a)) => Exp (Complex a) -> Exp a
-- magnitude (unlift -> r :+ i) = sqrt (r*r + i*i)
magnitude (unlift -> r :+ i) = scaleFloat k (sqrt (sqr (scaleFloat mk r) + sqr (scaleFloat mk i)))
  where
    k     = max (exponent r) (exponent i)
    mk    = -k
    sqr z = z * z

-- | The phase of a complex number, in the range @(-'pi', 'pi']@. If the
-- magnitude is zero, then so is the phase.
--
phase :: (RealFloat a, Elt (Complex a)) => Exp (Complex a) -> Exp a
phase z@(unlift -> r :+ i) =
  if z == 0
    then 0
    else atan2 i r

-- | The function 'polar' takes a complex number and returns a (magnitude,
-- phase) pair in canonical form: the magnitude is non-negative, and the phase
-- in the range @(-'pi', 'pi']@; if the magnitude is zero, then so is the phase.
--
polar :: (RealFloat a, Elt (Complex a)) => Exp (Complex a) -> Exp (a,a)
polar z =  lift (magnitude z, phase z)

-- | Form a complex number from polar components of magnitude and phase.
--
#if __GLASGOW_HASKELL__ <= 708
mkPolar :: forall a. (RealFloat a, Elt (Complex a)) => Exp a -> Exp a -> Exp (Complex a)
#else
mkPolar :: forall a. (Floating a,  Elt (Complex a)) => Exp a -> Exp a -> Exp (Complex a)
#endif
mkPolar = lift2 (C.mkPolar :: Exp a -> Exp a -> Complex (Exp a))

-- | @'cis' t@ is a complex value with magnitude @1@ and phase @t@ (modulo
-- @2*'pi'@).
--
#if __GLASGOW_HASKELL__ <= 708
cis :: forall a. (RealFloat a, Elt (Complex a)) => Exp a -> Exp (Complex a)
#else
cis :: forall a. (Floating a,  Elt (Complex a)) => Exp a -> Exp (Complex a)
#endif
cis = lift1 (C.cis :: Exp a -> Complex (Exp a))

-- | Return the real part of a complex number
--
real :: (Elt a, Elt (Complex a)) => Exp (Complex a) -> Exp a
real (unlift -> r :+ _) = r

-- | Return the imaginary part of a complex number
--
imag :: (Elt a, Elt (Complex a)) => Exp (Complex a) -> Exp a
imag (unlift -> _ :+ i) = i

-- | Return the complex conjugate of a complex number, defined as
--
-- > conjugate(Z) = X - iY
--
conjugate :: (Num a, Elt (Complex a)) => Exp (Complex a) -> Exp (Complex a)
conjugate z = lift $ real z :+ (- imag z)

ext: hs

max_stars_repo_name: ocramz/accelerate

max_issues_repo_name: ocramz/accelerate

max_forks_repo_name: ocramz/accelerate


text:
{-# LANGUAGE RecordWildCards #-}

module Data.Matrix.CSC where

import Control.Arrow (second)
import Control.Exception
import Data.Maybe
import qualified Data.Vector.Storable as VS
import qualified Data.Vector as VB
import Foreign.C.Types
import Numeric.LinearAlgebra

data CSC a = CSC
    { cscVals :: VS.Vector a
    , cscCols :: VS.Vector CInt
    , cscRows :: VS.Vector CInt
    , cscNRows :: CInt
    , cscNCols :: CInt
    } deriving (Eq, Show)

fromMatrix :: (Element a, VS.Storable a, Eq a, Num a) => Matrix a -> CSC a
fromMatrix mat = CSC{..}
  where
    cscNRows = fromIntegral $ rows mat
    cscNCols = fromIntegral $ cols mat
    cscRows = VS.fromList $ concat irs
    cscVals = VS.fromList $ concat vals
    cscCols = VS.fromList $ scanl (\a b -> a + fromIntegral (length b)) 0 vals
    (irs, vals) = unzip $ map (go . VS.toList) (toColumns mat)
    go col = unzip $ catMaybes $ zipWith (\idx val -> if val == 0 then Nothing else Just (fromIntegral idx, val)) [0..] col

toMatrix :: (Element a, Num a, VS.Storable a) => CSC a -> Matrix a
toMatrix CSC {..} = fromColumns cols
  where
    cols = zipWith go rows vals
    go rs vs = VS.replicate (fromIntegral cscNRows) 0 VS.//
        (zipWith (\r v -> (fromIntegral r, v)) (VS.toList rs) (VS.toList vs))
    rows = split cscRows cscCols
    vals = split cscVals cscCols
    split v = snd . VS.foldr splitter (v, []) . VS.init
      where
        splitter idx (v, xs) = second (:xs) $ VS.splitAt (fromIntegral idx) v

(===) :: (VS.Storable a) => CSC a -> CSC a -> CSC a
a === b = assert (cscNCols a == cscNCols b) $ CSC
    { cscNCols = cscNCols a
    , cscNRows = cscNRows a + cscNRows b
    , cscRows = rows
    , cscVals = vals
    , cscCols = cols
    }
  where
    cols = VS.zipWith (+) (cscCols a) (cscCols b)
    rows = VS.concat $ concat $
        zipWith (\a b -> [a, b]) (split (cscRows a) (cscCols a))
        (split (VS.map (+ cscNRows a) $ cscRows b) (cscCols b))
    vals = VS.concat $ concat $
        zipWith (\a b -> [a, b]) (split (cscVals a) (cscCols a))
        (split (cscVals b) (cscCols b))
    split v = snd . VS.foldr splitter (v, [])
      where
        splitter idx (v, xs) = second (:xs) $ VS.splitAt (fromIntegral idx) v

(|||) :: (VS.Storable a) => CSC a -> CSC a -> CSC a
a ||| b = assert (cscNRows a == cscNRows b) $ CSC
    { cscNCols = cscNCols a + cscNCols b
    , cscNRows = cscNRows a
    , cscRows = cscRows a VS.++ cscRows b
    , cscVals = cscVals a VS.++ cscVals b
    , cscCols = cscCols a VS.++ (VS.map (+ VS.last (cscCols a)) $ VS.tail (cscCols b))
    }

diagBlock :: (VS.Storable a) => CSC a -> CSC a -> CSC a
diagBlock a b = CSC
    { cscNCols = cscNCols a + cscNCols b
    , cscNRows = cscNRows a + cscNRows b
    , cscRows = cscRows a VS.++ VS.map (+ cscNRows a) (cscRows b)
    , cscVals = cscVals a VS.++ cscVals b
    , cscCols = cscCols a VS.++ (VS.map (+ VS.last (cscCols a)) $ VS.tail (cscCols b))
    }

zero :: (Integral a, VS.Storable b) => a -> a -> CSC b
zero rows cols = CSC
    { cscNCols = fromIntegral cols
    , cscNRows = fromIntegral rows
    , cscRows = VS.empty
    , cscVals = VS.empty
    , cscCols = VS.replicate (fromIntegral cols + 1) 0
    }

fromBlocks :: (VS.Storable a) => [[CSC a]] -> CSC a
fromBlocks = foldr1 (===) . map (foldr1 (|||))

scale :: (VS.Storable a, Num a) => a -> CSC a -> CSC a
scale mult csc = csc { cscVals = VS.map (* mult) (cscVals csc) }


text:
{-|
Module: Numeric.Morpheus
Description: Low-level machine learning auxiliary functions
Copyright: (c) Alexander Ignatyev, 2017
License: BSD-3
Stability: experimental
Portability: POSIX

Purely functional interface to morpheus based on hmatrix. Morpheus library contains a bunch of cache line aware numerical algorithms suitable for using as a low-level primitives to build machine learning applications.

Naive implementations, which do not take into account cache lines, would suffer order of magnitude performance degradation if tried to traverse a row major matrix in column order.

The library functions take into account CPU cache and work well for both row major and column major matrices, as you can see from the benchmark below which measure time for finding indices of max elements for every column and every row in comparison with naive approach:

@
| Benchmark                           | Morheus   | Naive    |
+-------------------------------------+-----------+----------+
| max index in columns - row major    | 0.9077 ms | 12.93 ms |
+-------------------------------------+-----------+----------+
| max index in columns - column major | 0.8778 ms | 1.490 ms |
+-------------------------------------+-----------+----------+
| max index in rows - row major       | 0.8622 ms | 1.504 ms |
+-------------------------------------+-----------+----------+
| max index in rows - column major    | 0.8913 ms | 12.80 ms |
+-------------------------------------+-----------+----------+
@

== Examples

@
import Numeric.LinearAlgebra
import Numeric.Morpheus

a = matrix 5 [
  71, 11, 3,  -9, -7,
  21, -7,  -2,  23, 11,
  -11, 32, 53, -49, 37,
  1,  -24, 78, 90, 17
  ]

main :: IO ()
main = do
  putStrLn "\nNumeric.Morpheus.MatrixReduce functions: "

  putStr "sum of elements of rows: "
  print $ rowSum a

  putStr "product of elements of rows: "
  print $ rowPredicate (*) a

  putStr "max values of columns: "
  print $ fst $ columnMaxIndex a

  putStr "indices of max values of columns: "
  print $ snd $ columnMaxIndex a

  putStrLn "\n\nNumeric.Morpheus.Activation functions: "

  putStrLn "\nSigmoid:"
  disp 3 (sigmoid a)

  putStrLn "ReLu:"
  disp 3 (relu a)

  putStrLn "ReLu gradient:"
  disp 3 (reluGradient a)

  putStrLn "\n\nNumeric.Morpheus.Statistics functions: "

  putStrLn "\nColumn Means:"
  let colmeans = columnMean a
  print colmeans

  putStrLn "Column Standard Deviations: "
  print (columnStddev_m colmeans a)
@
-}

module Numeric.Morpheus
(
  module Numeric.Morpheus.Activation
  , module Numeric.Morpheus.MatrixReduce
  , module Numeric.Morpheus.Statistics
)

where

import Numeric.Morpheus.Activation
import Numeric.Morpheus.MatrixReduce
import Numeric.Morpheus.Statistics

max_stars_repo_name: Alexander-Ignatyev/morpheus

max_issues_repo_name: aligusnet/morpheus

max_forks_repo_name: aligusnet/morpheus


text:
{-# LANGUAGE CPP #-}

{-# LANGUAGE ConstraintKinds #-} -- 'Functor' class.
{-# LANGUAGE TypeFamilies    #-} -- 'Functor' class.

{-# LANGUAGE UndecidableInstances #-} -- Required for some of the 'transformer' instances.

{-# LANGUAGE TypeOperators #-} -- For ':*:' instance and others.

-- | Definition of constrained functors as they are required to work with
--   constrained monads and constrained supermonads.
module Control.Super.Monad.Constrained.Functor
  ( Functor(..)
  ) where

import Prelude
  ( Ord
  , (.), ($), const
  )

import GHC.Exts ( Constraint )

import qualified Prelude as P

-- To define instances:
import Data.Functor.Identity ( Identity )

import qualified Data.Monoid as Mon
import qualified Data.Proxy as Proxy
import qualified Data.Functor.Product as Product
import qualified Data.Functor.Compose as Compose
#if MIN_VERSION_GLASGOW_HASKELL(8,0,0,0)
import qualified Data.Semigroup as Semigroup
import qualified Data.Complex as Complex
import qualified Data.List.NonEmpty as NonEmpty
#endif

import qualified Control.Arrow as Arrow
import qualified Control.Applicative as App
import qualified Control.Monad.ST as ST
import qualified Control.Monad.ST.Lazy as STL

import qualified Text.ParserCombinators.ReadP as Read
import qualified Text.ParserCombinators.ReadPrec as Read

import qualified GHC.Conc as STM
#if MIN_VERSION_GLASGOW_HASKELL(8,0,0,0)
import qualified GHC.Generics as Generics
#endif

-- To defined constrained instances:
import qualified Data.Set as S

-- To define "transformers" instances:
import qualified Control.Monad.Trans.Cont     as Cont
import qualified Control.Monad.Trans.Except   as Except
import qualified Control.Monad.Trans.Identity as Identity
import qualified Control.Monad.Trans.Maybe    as Maybe
import qualified Control.Monad.Trans.RWS.Lazy      as RWSL
import qualified Control.Monad.Trans.RWS.Strict    as RWSS
import qualified Control.Monad.Trans.Reader        as Reader
import qualified Control.Monad.Trans.State.Lazy    as StateL
import qualified Control.Monad.Trans.State.Strict  as StateS
import qualified Control.Monad.Trans.Writer.Lazy   as WriterL
import qualified Control.Monad.Trans.Writer.Strict as WriterS

infixl 4  <$

-- -----------------------------------------------------------------------------
-- Constrained functor class
-- -----------------------------------------------------------------------------

-- | Class for constrained functors. Obeys all of the same laws as the standard
--   'Prelude.Functor' class, but allows to constrain the functors result type.
class Functor f where
  type FunctorCts f (a :: *) (b :: *) :: Constraint
  type FunctorCts f a b = ()

  fmap :: (FunctorCts f a b) => (a -> b) -> f a -> f b
  (<$) :: (FunctorCts f b a) => a -> f b -> f a
  (<$) = fmap . const

-- Unconstrained instances -----------------------------------------------------

instance Functor ((->) r) where
  fmap = P.fmap
  (<$) = (P.<$)
instance Functor Identity where
  fmap = P.fmap
  (<$) = (P.<$)
instance Functor [] where
  fmap = P.fmap
  (<$) = (P.<$)
instance Functor P.Maybe where
  fmap = P.fmap
  (<$) = (P.<$)
instance Functor P.IO where
  fmap = P.fmap
  (<$) = (P.<$)
instance Functor (P.Either e) where
  fmap = P.fmap
  (<$) = (P.<$)

instance Functor Mon.First where
  fmap = P.fmap
  (<$) = (P.<$)
instance Functor Mon.Last where
  fmap = P.fmap
  (<$) = (P.<$)
#if MIN_VERSION_GLASGOW_HASKELL(8,0,0,0)
instance Functor Mon.Sum where
  fmap = P.fmap
  (<$) = (P.<$)
instance Functor Mon.Product where
  fmap = P.fmap
  (<$) = (P.<$)
instance Functor Mon.Dual where
  fmap = P.fmap
  (<$) = (P.<$)
#endif
instance (Functor f) => Functor (Mon.Alt f) where
  type FunctorCts (Mon.Alt f) a b = FunctorCts f a b
  fmap f (Mon.Alt ma) = Mon.Alt $ fmap f ma
  a <$ (Mon.Alt mb) = Mon.Alt $ a <$ mb

#if MIN_VERSION_GLASGOW_HASKELL(8,0,0,0)
instance Functor Semigroup.Min where
  fmap = P.fmap
  (<$) = (P.<$)
instance Functor Semigroup.Max where
  fmap = P.fmap
  (<$) = (P.<$)
instance Functor Semigroup.Option where
  fmap = P.fmap
  (<$) = (P.<$)
instance Functor Semigroup.First where
  fmap = P.fmap
  (<$) = (P.<$)
instance Functor Semigroup.Last where
  fmap = P.fmap
  (<$) = (P.<$)
#endif

instance Functor Proxy.Proxy where
  fmap = P.fmap
  (<$) = (P.<$)
#if MIN_VERSION_GLASGOW_HASKELL(8,0,0,0)
instance Functor Complex.Complex where
  fmap = P.fmap
  (<$) = (P.<$)
instance Functor NonEmpty.NonEmpty where
  fmap = P.fmap
  (<$) = (P.<$)
#endif

instance (Functor f, Functor g) => Functor (Product.Product f g) where
  type FunctorCts (Product.Product f g) a b = (FunctorCts f a b, FunctorCts g a b)
  fmap f (Product.Pair fa fb) = Product.Pair (fmap f fa) (fmap f fb)

instance (Functor f, Functor g) => Functor (Compose.Compose f g) where
  type FunctorCts (Compose.Compose f g) a b = (FunctorCts f (g a) (g b), FunctorCts g a b)
  fmap f (Compose.Compose x) = Compose.Compose (fmap (fmap f) x)

instance Functor Read.ReadP where
  fmap = P.fmap
  (<$) = (P.<$)
instance Functor Read.ReadPrec where
  fmap = P.fmap
  (<$) = (P.<$)

instance Functor (ST.ST s) where
  fmap = P.fmap
  (<$) = (P.<$)
instance Functor (STL.ST s) where
  fmap = P.fmap
  (<$) = (P.<$)
instance (Arrow.ArrowApply a) => Functor (Arrow.ArrowMonad a) where
  fmap = P.fmap
  (<$) = (P.<$)
instance (Functor m) => Functor (App.WrappedMonad m) where
  type FunctorCts (App.WrappedMonad m) a b = FunctorCts m a b
  fmap f (App.WrapMonad ma) = App.WrapMonad $ fmap f ma
  a <$ (App.WrapMonad mb) = App.WrapMonad $ a <$ mb

instance Functor STM.STM where
  fmap = P.fmap
  (<$) = (P.<$)

#if MIN_VERSION_GLASGOW_HASKELL(8,0,0,0)
instance Functor Generics.U1 where
  fmap = P.fmap
  (<$) = (P.<$)
instance (Functor f) => Functor (Generics.Rec1 f) where
  type FunctorCts (Generics.Rec1 f) a b = FunctorCts f a b
  fmap f (Generics.Rec1 ma) = Generics.Rec1 $ fmap f ma
  a <$ (Generics.Rec1 mb) = Generics.Rec1 $ a <$ mb

instance (Functor f, Functor g) => Functor (f Generics.:*: g) where
  type FunctorCts (f Generics.:*: g) a b = (FunctorCts f a b, FunctorCts g a b)
  fmap f (a Generics.:*: b) = (fmap f a) Generics.:*: (fmap f b)
instance (Functor f, Functor g) => Functor (f Generics.:.: g) where
  type FunctorCts (f Generics.:.: g) a b = (FunctorCts f (g a) (g b), FunctorCts g a b)
  fmap f (Generics.Comp1 ma) = Generics.Comp1 $ fmap (fmap f) ma
instance Functor f => Functor (Generics.M1 i c f) where
  type FunctorCts (Generics.M1 i c f) a b = FunctorCts f a b
  fmap f (Generics.M1 ma) = Generics.M1 $ fmap f ma
#endif

-- Constrained instances -------------------------------------------------------

instance Functor S.Set where
  type FunctorCts S.Set a b = Ord b
  fmap = S.map

-- "transformers" package instances: -------------------------------------------

-- Continuations are so wierd...
-- | TODO / FIXME: Still need to figure out how and if we can generalize the continuation implementation.
instance Functor (Cont.ContT r m) where
  fmap = P.fmap
  (<$) = (P.<$)

instance Functor m => Functor (Except.ExceptT e m) where
  type FunctorCts (Except.ExceptT e m) a b = FunctorCts m (P.Either e a) (P.Either e b)
  fmap f = Except.ExceptT . fmap (fmap f) . Except.runExceptT
  {-# INLINE fmap #-}

instance (Functor m) => Functor (Identity.IdentityT m) where
  type FunctorCts (Identity.IdentityT m) a b = FunctorCts m a b
  fmap f = Identity.mapIdentityT (fmap f)
  {-# INLINE fmap #-}

instance (Functor m) => Functor (Maybe.MaybeT m) where
  type FunctorCts (Maybe.MaybeT m) a b = FunctorCts m (P.Maybe a) (P.Maybe b)
  fmap f = Maybe.mapMaybeT (fmap (fmap f))
  {-# INLINE fmap #-}

instance (Functor m) => Functor (RWSL.RWST r w s m) where
  type FunctorCts (RWSL.RWST r w s m) a b = FunctorCts m (a, s, w) (b, s, w)
  fmap f m = RWSL.RWST $ \ r s ->
      fmap (\ ~(a, s', w) -> (f a, s', w)) $ RWSL.runRWST m r s
  {-# INLINE fmap #-}

instance (Functor m) => Functor (RWSS.RWST r w s m) where
  type FunctorCts (RWSS.RWST r w s m) a b = FunctorCts m (a, s, w) (b, s, w)
  fmap f m = RWSS.RWST $ \ r s ->
      fmap (\ (a, s', w) -> (f a, s', w)) $ RWSS.runRWST m r s
  {-# INLINE fmap #-}

instance (Functor m) => Functor (Reader.ReaderT r m) where
  type FunctorCts (Reader.ReaderT r m) a b = FunctorCts m a b
  fmap f  = Reader.mapReaderT (fmap f)
  {-# INLINE fmap #-}

instance (Functor m) => Functor (StateL.StateT s m) where
  type FunctorCts (StateL.StateT s m) a b = FunctorCts m (a, s) (b, s)
  fmap f m = StateL.StateT $ \ s ->
      fmap (\ ~(a, s') -> (f a, s')) $ StateL.runStateT m s
  {-# INLINE fmap #-}

instance (Functor m) => Functor (StateS.StateT s m) where
  type FunctorCts (StateS.StateT s m) a b = FunctorCts m (a, s) (b, s)
  fmap f m = StateS.StateT $ \ s ->
      fmap (\ (a, s') -> (f a, s')) $ StateS.runStateT m s
  {-# INLINE fmap #-}

instance (Functor m) => Functor (WriterL.WriterT w m) where
  type FunctorCts (WriterL.WriterT w m) a b = FunctorCts m (a, w) (b, w)
  fmap f = WriterL.mapWriterT $ fmap $ \ ~(a, w) -> (f a, w)
  {-# INLINE fmap #-}

instance (Functor m) => Functor (WriterS.WriterT w m) where
  type FunctorCts (WriterS.WriterT w m) a b = FunctorCts m (a, w) (b, w)
  fmap f = WriterS.mapWriterT $ fmap $ \ (a, w) -> (f a, w)
  {-# INLINE fmap #-}

max_stars_repo_name: jbracker/supermonad-plugin

max_issues_repo_name: jbracker/supermonad

max_forks_repo_name: jbracker/supermonad-plugin


text:
{-# OPTIONS_GHC -fno-warn-type-defaults #-}
{-# LANGUAGE FlexibleInstances    #-}
{-# LANGUAGE TypeSynonymInstances #-}
-- ProcessingTime.hs ---
--
-- Filename: ProcessingTime.hs
-- Description:
-- Author: Manuel Schneckenreither
-- Maintainer:
-- Created: Tue Nov 21 10:16:08 2017 (+0100)
-- Version:
-- Package-Requires: ()
-- Last-Updated:
--           By:
--     Update #: 69
-- URL:
-- Doc URL:
-- Keywords:
-- Compatibility:
--
--

-- Commentary:
--
--
--
--

-- Change Log:
--
--
--
--
--
--
--

-- Code:

module SimSim.ProcessingTime.Type where

import           ClassyPrelude
import qualified Data.Map.Strict                     as M
import           Statistics.Distribution
import           Statistics.Distribution.Exponential
import           System.Random.MWC

import           SimSim.Block
import           SimSim.ProductType
import           SimSim.Time

-- ^ Use for instance Statistics.Distribution.Exponential and Statistics.Distribution to generate random numbers.
type ProcessingTime = GenIO -> IO Time

type ProcessingTimes = M.Map Block (M.Map ProductType ProcessingTime)

type ProcTimes = [(Block, [(ProductType, ProcessingTime)])]

fromProcTimes :: ProcTimes -> ProcessingTimes
fromProcTimes xs = M.fromList $ fmap (second M.fromList) xs

-- d :: ExponentialDistribution
-- d = exponential (1/0.8)

-- test :: IO [Double]
-- test = do
--   g <- createSystemRandom
--   mapM (\_ -> genContVar d g) [1..7000]

--
-- ProcessingTime.hs ends here

max_stars_repo_name: schnecki/simsim

max_issues_repo_name: schnecki/simsim

max_forks_repo_name: schnecki/simsim


text:
{-# LANGUAGE Safe #-}
{-# OPTIONS_GHC -fno-warn-deprecations #-}
-- | Import all modules from prelude that should be safe
module Main where

import Numeric
import Prelude

import Control.Applicative
import Control.Arrow
import Control.Category

import Control.Concurrent
import Control.Concurrent.Chan
import Control.Concurrent.MVar
import Control.Concurrent.QSem
import Control.Concurrent.QSemN

import Control.Exception
import Control.Exception.Base

import Control.Monad
import Control.Monad.Fix
import Control.Monad.Instances
import Control.Monad.Zip

import Control.Monad.ST
import Control.Monad.ST.Lazy
import Control.Monad.ST.Strict

import Data.Bifunctor
import Data.Bits
import Data.Bool
import Data.Char
-- import Data.Coerce
import Data.Complex
import Data.Data
import Data.Dynamic
import Data.Either
import Data.Eq
import Data.Fixed
import Data.Foldable
import Data.Function
import Data.Functor.Identity
import Data.Functor
import Data.IORef
import Data.Int
import Data.Ix
import Data.List
import Data.Maybe
import Data.Monoid
import Data.Ord
import Data.Proxy
import Data.Ratio
import Data.String
import Data.Traversable
import Data.Tuple
import Data.Typeable
import Data.Unique
import Data.Version
import Data.Word

import Data.STRef
import Data.STRef.Lazy
import Data.STRef.Strict

import Data.Type.Bool
-- import Data.Type.Coercion
-- import Data.Type.Equality

-- import Debug.Trace

import Foreign

import Foreign.Concurrent
import Foreign.ForeignPtr
import Foreign.Ptr
import Foreign.StablePtr
import Foreign.Storable

import Foreign.C
import Foreign.C.Error
import Foreign.C.String
import Foreign.C.Types

import Foreign.Marshal
import Foreign.Marshal.Alloc
import Foreign.Marshal.Array
import Foreign.Marshal.Error
import Foreign.Marshal.Pool
import Foreign.Marshal.Utils

import System.CPUTime
import System.Console.GetOpt
import System.Environment
import System.Exit
import System.IO
import System.IO.Error
import System.Info
import System.Mem
import System.Mem.StableName
import System.Mem.Weak
import System.Posix.Internals
import System.Posix.Types
import System.Timeout

import Text.ParserCombinators.ReadP
import Text.ParserCombinators.ReadPrec
import Text.Printf
import Text.Read
import Text.Read.Lex
import Text.Show
import Text.Show.Functions

import Type.Reflection

-- import Unsafe.Coerce

f :: Int
f = 2

main :: IO ()
main = putStrLn $ "X is: " ++ show f

max_stars_repo_name: janpath/ghc

max_issues_repo_name: janpath/ghc

max_forks_repo_name: janpath/ghc


text:
{-# LANGUAGE OverloadedStrings, DuplicateRecordFields, RecordWildCards, NamedFieldPuns #-}
{-|
 - Module      : Numeric.GLP
 - Description : Specifiy a General Linear Problem and transform to a Standardized LP
 - Copyright   : (c) Torsten Kemps-Benedix, 2019
 - This module 'Numeric.GLP' specifiies a [Generalized Linear Problem](https://en.wikipedia.org/wiki/Linear_programming) (or LP in short) and provides means to tranform this into a Standardized LP (see module 'Numeric.LP') that may be solved with e.g. the simplex solver in module 'Numeric.LP.Simplex'.
-}
module Numeric.GLP
-- (GLP, Numeric.GLP.a, c, m1, m2, m3, dir, mkGLP, standardizeGLP)
where

import qualified Data.List                   as L
import           Control.Monad
import qualified Data.Set                    as S
import           Data.Set                    (member, isSubsetOf)
import           Numeric.LinearAlgebra       (I, Element, Matrix, Vector, ident, diagl,
                                              cols, rows, find, toColumns,
                                              fromColumns, cmap, fromList, toList,
                                              (!), (><),
                                              )
import qualified Numeric.LinearAlgebra       as LA
import           Numeric.LinearAlgebra.Devel (modifyMatrix, runSTMatrix,
                                              thawMatrix, mapVectorWithIndex, fi)
import           Numeric.LinearAlgebra       (vjoin, (|||), (===))

-- | This record defines a General Linear Problem containing ≤
-- and ≥ inequalities and == equalities. The sum of the numbers of
-- inequalities ('m1' times ≤ and 'm2' times ≥) and equalities ('m3'
-- times ==) must match one less than the number of
-- rows of the matrix 'a'. The first (i.e. the zero-th) row of the matrix
-- contains the coefficients of the objective function.
--
-- We assume that \(a_{00}\) is zero. The matrix \(a\) shall have \(m+1\) rows
-- and \(n+1\) columns.
--
-- The formal specification of the problem we intend to solve is as follows:
--
-- Maximize resp. minimize (according to 'dir') the following linear objective
-- function:
--
-- \[ z_0 = \sum_{i=1}^{n} a_{0,i} x_i \]
--
-- Subject to the following linear constraints:
--
-- \[ z_j = a_{j,0} + \sum_{i=1}^{n} a_{j,i} x_i,\, j=1,\ldots,m\]
-- \[ z_j \leq 0,\, j=1,\ldots,m_1\]
-- \[ z_j \geq 0,\, j=m_1+1,\ldots,m_1+m_2\]
-- \[ z_j = 0,\, j=m_1+m_2+1,\ldots,m=m1+m_2+m_3\]
--
-- The original variables \(x_1, \ldots, x_n\) may be constraint to be non-negative,
-- non-positive or they may be unconstrained according to the n-vector
-- \(c\in\{-1,0,+1\}^n\).
--
-- \[ c_j = -1 \Rightarrow x_j \leq 0\]
-- \[ c_j = 1 \Rightarrow x_j \geq 0\]
-- \[ c_j = 0 \Rightarrow x_j \text{ unrestricted} \]
data GLP = GLP {
    a   :: Matrix Double, -- ^ Tableau of coefficients, see above for a formal
                          -- description of their meaning
    c   :: Vector I,      -- ^ Type of constraint for \(x_j\), \(j=1,\ldots,n\),
                          -- see 'GLP' for further details. The element \(c_i\)
                          -- corresponds to column \(i+1\) of the matrix \(a\).
    n :: Int,             -- ^ Number of original variables \(x_1, \ldots, x_n\)
    m1  :: Int,           -- ^ Number of ≤ inequalities
    m2  :: Int,           -- ^ Number of ≥ inequalities
    m3  :: Int,           -- ^ Number of = equalities
    dir :: Dir,           -- ^ Whether the objective function is to be maximized
                          -- or minimzed
    lhs, rhs :: Vector I  -- ^ Indices of left and right hand variables respectively.
                          -- The original variables \(x_1,\ldots,x_n\) have indices
                          -- \(1,\ldots,n\), the slack variables get indices
                          -- \(n+1, \ldots, n+m_1+m_2+1\), and the auxiliary variables
                          -- finally get indices \(n+m_1+m_2+2, \ldots, n+m_1+m_2+m+2\).
    } deriving (Show)

-- | This smart constructor checks all the pre-conditions for 'GLP' stated in the
-- description of 'GLP' and initializes the variables 'lhs' and 'rhs'. Initially only
-- the \(m\) auxiliary variables \(z_1,\ldots,z_m\) are on the left hand side and all the
-- original variables \(x_1,\ldots,x_n\) are on the right hand side. Initially there are
-- no slack variables \(y_1,\ldots,y_{m_1+m_2+1}\). These will be introdced by the function
-- 'introduceSlackVariables' below.
mkGLP :: Matrix Double -> Vector I -> (Int, Int, Int) -> Dir -> Either String GLP
mkGLP a c (m1,m2,m3) dir =
    let m = rows a - 1
        n = cols a - 1
    in if (m /= m1+m2+m3)
       then Left $ "rows a - 1 == "++show m++" /= m1+m2+m3 == "++show (m1+m2+m3)
       else if (m <= 0)
            then Left "The number of constraints is m <= 0"
            else if (n <= 0)
                 then Left "The number of variables n is <= 0"
                 else if not (S.fromList (LA.toList c) `isSubsetOf` (S.fromList [-1,0,1]))
                      then Left "Vector c has elements other than -1, 0, +1"
                      else if LA.size c /= n
                           then Left $ "Vector c has not the proper length "++show n
                           else Right (GLP a c n m1 m2 m3 dir
                                           (LA.fromList [fi (n+m1+m2+2)..fi (n+m1+m2+m+2)])
                                           (LA.fromList [1..fi n]))

data Dir = Maximize | Minimize deriving Show

-- | Transform the Generalized Linear Problem to a Standardized Linear Problem
-- by means of the following steps:
--
-- 1. If the objective function is to be minimized multiply it by -1.
-- 2. Turn around the ≥-inequalities by negation of the corresponding rows. After this
-- step there are 'm1'+'m2' inequalities of the <= type.
-- 3. For those variables \(x_i\) that have a constraint \(x_i \leq 0\) substitute
-- \(x_i \rightarrow -x_i\). This means multiplying the column \(i\) of the matrix
-- \(a\) by \(-1\).
-- 4. For those variables \(x_i\) that are unconstrained introduce new non-negative
-- variables \(x_i = x_i⁺ - x_i⁻\). We give these the indices \(x_i⁺ = x_i\) and
-- \(x_⁻ = x_{i+1}\). The indices of all other variables \(x_k\), \(k>i\) have to be
-- increased by \(1\). The matrix \(a\) therefore grows by one column.
-- 5. For the 'm1'+'m2' ≤-inequalities introduce as many slack variables on the right
-- hand side and make these equalities instead. After this step there are only
-- equalities left as constraints. We do not represent these slack variables in the
-- matrix 'a' and therefore nothing really changes in this step.
-- 6. Introduce artificial variables \(z_j\) ( \(j=1,\ldots,m\) ) for each equation.
-- This changes the optimization problem and the solver has to take this into account.
--
-- The function returns a modified and standardized GLP and a matrix \(b\) that maps
-- solution vectors (including the objective functions value) of the modified GLP back
-- to solution vectors of the original GLP.
standardizeGLP :: GLP -> (GLP, Matrix Double)
standardizeGLP glp =
    let (glp1, b1) = resolveMin2Max glp
        (glp2, b2) = resolveGeqConstraints glp1
        (glp3, b3) = resolveNonPositiveVars glp2
        (glp4, b4) = resolveUnconstrainedVars glp3
        glp5 = introduceSlackVariables glp4
    in (glp5, b1 <> b2 <> b3 <> b4)

resolveMin2Max :: GLP -> (GLP, Matrix Double)
resolveMin2Max glp@GLP{..} = case dir of
        Maximize -> (glp, ident (1+m1+m2+m3))
        Minimize -> (glp{a = mapRow a 0 negate, dir = Maximize},
                     diagl (-1:(replicate (m1+m2+m3) 1)))

resolveGeqConstraints :: GLP -> (GLP, Matrix Double)
resolveGeqConstraints glp@GLP{..} =
    let a' = foldl (\a j -> mapRow a j negate) a [m1+1..m1+m2] -- negate rows from m1+1 to m1+m2
    in (glp{a = a', m1 = m1+m2, m2 = 0},
        ident (1+m1+m2+m3))

resolveNonPositiveVars :: GLP -> (GLP, Matrix Double)
resolveNonPositiveVars glp@GLP{..} =
    let varsToBeNegated = find (== -1) c
        varsToBeNegatedS = S.fromList varsToBeNegated
        (a', c') = foldl (\(a,c) i ->
                            (mapCol a (i+1) negate,
                             mapVectorWithIndex (\k x -> if k==i
                                                         then 1
                                                         else c!k)
                                                c))
                         (a, c)
                         (find (== -1) c)
    in (glp{a = a', c = c'},
        diagl $ 1:[if i `member` varsToBeNegatedS
                   then -1
                   else 1 | i <- [0..m1+m2+m3-1]])

resolveUnconstrainedVars :: GLP -> (GLP, Matrix Double)
resolveUnconstrainedVars glp@GLP{..} =
    let m = 1+m1+m2+m3
        zCol = vec1 m 0 1
        (a', c', b) =
            let listOfColsAndCs =
                   fmap (\(c1,ci,i1) ->
                            if ci==0
                            then ([c1, cmap negate c1],
                                  [1, 1],
                                  [vec1 m (i1+1) 1,
                                   vec1 m (i1+1) (-1)]) -- duplicate and adjust
                            else ([c1], [ci],
                                  [vec1 m (i1+1) 1])) -- return unchanged
                        (L.zip3 (toColumns a) (toList c) [0..(cols a)-1])
                columns = L.concat $ fmap (\(x,_,_) -> x) listOfColsAndCs
                cs = L.concat $ fmap (\(_,x,_) -> x) listOfColsAndCs
                bs = L.concat $ fmap (\(_,_,x) -> x) listOfColsAndCs
            in (fromColumns columns, fromList cs, fromColumns $ zCol:bs)
    in (glp{a=a', c=c'}, b)

introduceSlackVariables :: GLP -> GLP
introduceSlackVariables glp@GLP{..} =
    let m3' = m1+m2+m3
        n = cols a - 1
        a' = a ||| ((1><(m1+m2)) (replicate (m1+m2) 0) ===
                    diagl ((replicate m1 1)++(replicate m2 (-1))) ===
                    (m3><(m1+m2)) (replicate (m3*(m1+m2)) 0))
        rhs' = vjoin [rhs, LA.fromList [fi (n+1)..fi (n+1+m1+m2)]]
    in glp{a=a', m1=0, m2=0, m3=m3', rhs=rhs'}

vec1 :: Int -> Int -> Double -> Vector Double
vec1 n i x = fromList $ L.replicate i 0++x:L.replicate (n-i-1) 0

mapRow :: (Element t) => Matrix t -> Int -> (t -> t) -> Matrix t
mapRow a k f = runSTMatrix $ do
                    a' <- thawMatrix a
                    forM_ [0..cols a-1] (\i -> modifyMatrix a' k i f)
                    return a'

mapCol :: (Element t) => Matrix t -> Int -> (t -> t) -> Matrix t
mapCol a k f = runSTMatrix $ do
                    a' <- thawMatrix a
                    forM_ [0..rows a-1] (\i -> modifyMatrix a' i k f)
                    return a'

prettyPrintAsFormulas :: GLP -> String
prettyPrintAsFormulas glp@GLP{..} =
    let rs = LA.toRows a
        m = cols a - 1
        n = rows a - 1
        objFun = (show dir)++
                    " z_0 = "++linExp 0
        zRows = fmap (\k -> "z_"++show k++" = "++linExp k) [1..n]
        linExp k = L.intercalate " + " (
                    L.filter (/= "") $
                    fmap (\(a,i) ->
                            if a /= 0
                            then show a++(if i>0 then " "++varName i else "")
                            else "") (L.zip (toList $ rs!!k) [0..m]))
        varName i = if i <= n then "x_"++show i  else "y_"++show (i-n)
        zConstr  = L.intercalate ", " $ filter (/= "") $
                   [if m1 == 1
                    then "z_1 ≤ 0"
                    else if m1 == 2
                         then "z_1, z_2 ≤ 0"
                         else if m1 > 2
                              then "z_1, ..., z_"++show m1++" ≤ 0"
                              else "",
                    if m2 == 1
                    then "z_"++show (m1+1)++" ≥ 0"
                    else if m2 > 1
                         then "z_"++show (m1+1)++", ..., z_"++show (m1+m2)++" ≥ 0"
                         else "",
                    if m3 == 1
                    then "z_"++show (m1+m2+1)++" = 0"
                    else if m3 > 1
                         then "z_"++show (m1+m2+1)++", ..., z_"++show (m1+m2+m3)++" = 0"
                         else ""]
        xConstr ck k = "x_"++show k++if ck == -1
                                     then " ≤ 0"
                                     else if ck == 1
                                          then " ≥ 0"
                                          else " unconstrained"
        xConstrs = fmap (uncurry xConstr) (L.zip (toList c) [1..m])
    in L.intercalate "\n" (
            objFun:
            "Subject to":
            (fmap (\line -> "\t"++line) (zRows++[""]++
                                         [zConstr, ""]++
                                         xConstrs)))
            ++"\n"

(Right ex1') = ex1
ex1 = let a = (5><5)
                 [   0,    1,  1,  3, -0.5
                 , 740,   -1,  0, -2,  0
                 ,   0,    0, -2,  0,  7
                 ,   0.5,  0, -1,  1, -2
                 ,   9,   -1, -1, -1, -1]
          c = LA.fromList [-1, 1, 0, -1]
      in mkGLP a c (2, 1, 1) Minimize

allEx1 = do
    let g1 = ex1'
    putStr $ prettyPrintAsFormulas g1
    putStrLn $ "lhs = "++show (lhs g1)
    putStrLn $ "rhs = "++show (rhs g1)
    putStrLn "\nResolve min/max"
    let (g2, b1) = resolveMin2Max g1
    putStr $ prettyPrintAsFormulas g2
    putStrLn $ "lhs = "++show (lhs g2)
    putStrLn $ "rhs = "++show (rhs g2)
    print b1
    putStrLn "\nResove geq constraints"
    let (g3, b2) = resolveGeqConstraints g2
    putStr $ prettyPrintAsFormulas g3
    putStrLn $ "lhs = "++show (lhs g3)
    putStrLn $ "rhs = "++show (rhs g3)
    print b2
    putStrLn "\nResove non-positive variables"
    let (g4, b3) = resolveNonPositiveVars g3
    putStr $ prettyPrintAsFormulas g4
    putStrLn $ "lhs = "++show (lhs g4)
    putStrLn $ "rhs = "++show (rhs g4)
    print b3
    putStrLn "\nResove unconstrained variables"
    let (g5, b4) = resolveUnconstrainedVars g4
    putStr $ prettyPrintAsFormulas g5
    putStrLn $ "lhs = "++show (lhs g5)
    putStrLn $ "rhs = "++show (rhs g5)
    print b4
    putStrLn "\nIntroduce slack variables"
    let g6 = introduceSlackVariables g5
    putStr $ prettyPrintAsFormulas g6
    putStrLn $ "lhs = "++show (lhs g6)
    putStrLn $ "rhs = "++show (rhs g6)
    putStrLn "\n"
    return (g6, b1 <> b2 <> b3 <> b4)

max_stars_repo_name: tkemps/int-prog

max_issues_repo_name: tkemps/int-prog

max_forks_repo_name: tkemps/int-prog


text:
import Data.Word
import Data.List(elemIndex)
import Data.Maybe
import Data.Complex
import qualified Data.Array.Repa as R
import Data.Array.Repa.IO.BMP
import System.Environment

type Number = Double
type ComplexNumber = Complex Number
type IterationCount = Integer
type PointData = (ComplexNumber, IterationCount)
type RootData = (Int, IterationCount)
type Bounds = (Number, Number, Number, Number)
type Step = Number
type Color = (Word8, Word8, Word8)

e :: ComplexNumber
e = (exp 1) :+ 0

data Function = Constant ComplexNumber
              | Z
              | Sum Function Function
              | Product Function Function
              | Quotient Function Function
              | Power ComplexNumber Function
              | Sin Function
              | Cos Function
              | Exponential ComplexNumber Function
              | Log ComplexNumber Function
              | Undefined
              | Indeterminate
              deriving (Show)

realConstant :: Number -> Function
realConstant n = Constant (n :+ 0)

derivative :: Function -> Function
derivative (Constant _) = realConstant 0
derivative Z = realConstant 1
derivative (Sum u v) = Sum (derivative u) (derivative v)
derivative (Product u v) = Sum (Product (derivative u) v) (Product u (derivative v))
derivative (Quotient u v) = (Quotient (Sum (Product v (derivative u)) (Product (Product (realConstant (-1)) u) (derivative v))) (Power 2 v))
derivative (Power p f) = Product (Product (Constant p) (Power (p - 1) f)) (derivative f)
derivative (Sin f) = Product (Cos f) (derivative f)
derivative (Cos f) = Product (Product (realConstant (-1)) (Sin f)) (derivative f)
derivative (Exponential b f) = Product (Product (Exponential b f) (Log e (Constant b))) (derivative f)
derivative (Log b f) = Quotient (derivative f) (Product f (Log e (Constant b)))
derivative Undefined = Undefined
derivative Indeterminate = Indeterminate

evaluate :: Function -> ComplexNumber -> ComplexNumber
evaluate (Constant a) z = a
evaluate Z z = z
evaluate (Sum u v) z = (evaluate u z) + (evaluate v z)
evaluate (Product u v) z = (evaluate u z) * (evaluate v z)
evaluate (Quotient u v) z = (evaluate u z) / (evaluate v z)
evaluate (Power p f) z = (evaluate f z) ** p
evaluate (Sin f) z = sin (evaluate f z)
evaluate (Cos f) z = cos (evaluate f z)
evaluate (Exponential b f) z = b ** (evaluate f z)
evaluate (Log b f) z = logBase b (evaluate f z)
evaluate Undefined z = undefined
evaluate Indeterminate z = undefined

maxIterations :: IterationCount
maxIterations = 100

epsilon :: Number
epsilon = 1e-12

newton :: Function -> Function -> ComplexNumber -> PointData
newton f f' z0 = newton' z1 z0 0
    where approximate z = z - (evaluate f z) / (evaluate f' z)
          z1 = approximate z0
          newton' current previous n
              | n > maxIterations = (current, n)
              | magnitude (current - previous) < epsilon = (current, n)
              | otherwise = newton' (approximate current) current (n + 1)

windowApply :: (ComplexNumber -> PointData) -> Bounds -> Step -> [PointData]
windowApply f (top, left, bottom, right) step = [f (x :+ y ) | y <- ys, x <- xs]
    where xs = [left, left + step .. right]
          ys = [top, top + step .. bottom]

approximateIndex :: ComplexNumber -> [ComplexNumber] -> Maybe Int
approximateIndex z zs = elemIndex True [magnitude (z - zi) < epsilon * 2 | zi <- zs]

closestIndex :: ComplexNumber -> [ComplexNumber] -> Int
closestIndex z zs = snd . minimum $ zip [magnitude (z - zi) | zi <- zs] [0..]

closeEnough :: ComplexNumber -> [ComplexNumber] -> Bool
closeEnough z zs = case approximateIndex z zs of
    Just _ -> True
    Nothing ->  False

roots :: [PointData] -> [ComplexNumber]
roots ps = foldr checkRoot [] zs
    where checkRoot z roots
              | z `closeEnough` roots = roots
              | otherwise = z:roots
          zs = fst . unzip $ ps

assignRoot :: [ComplexNumber] -> PointData -> RootData
assignRoot rs p@(z, i) = (closestIndex z rs, i)

colorRoot :: RootData -> Color
colorRoot (n, i) = (ra, ga, ba)
    where colors = [ (0, 0, 0)
                   , (255, 0, 0)
                   , (0, 255, 0)
                   , (0, 0, 255)
                   , (255, 255, 0)
                   , (255, 0, 255)
                   , (0, 255, 255)
                   , (255, 255, 255)
                   ]
          (r, g, b) = colors !! n
          ii = max i 1
          ra = fromIntegral . min 255  $ 4 * r `div` ii
          ga = fromIntegral . min 255  $ 4 * g `div` ii
          ba = fromIntegral . min 255  $ 4 * b `div` ii

colorRoots :: Function -> Bounds -> Step -> [Color]
colorRoots f bounds step = map colorRoot rds
    where ps = windowApply (newton f (derivative f)) bounds step
          rs = roots ps
          rds = map (assignRoot rs) ps

arrayify :: Int -> Int -> [Color] -> R.Array R.U R.DIM2 Color
arrayify width height =
    R.fromListUnboxed (R.Z R.:. height R.:. width :: R.DIM2)

f = Sum (Power 3 Z) (realConstant (-1))

main = do
    args <- getArgs
    let [fileName, stringTop, stringLeft, stringBottom, stringRight, stringStep] = args
        top = read stringTop
        left = read stringLeft
        bottom = read stringBottom
        right = read stringRight
        step = read stringStep
        width = floor ((right - left) / step) + 1
        height = floor ((bottom - top) / step) + 1
        colorList = colorRoots f (top, left, bottom, right) step
        colorArray = arrayify width height colorList
    writeImageToBMP fileName colorArray


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

module Main where
import qualified Data.Foldable as Foldable

-- vector
import Data.Vector (Vector)
import qualified Data.Vector as Vector

-- Matrix
import Data.Matrix (Matrix)
import qualified Data.Matrix as Matrix

-- HMatrix
import qualified Numeric.LinearAlgebra.Data as HMatrix
import qualified Numeric.LinearAlgebra.HMatrix as HMatrix

-- bytestring
import Data.ByteString.Lazy (ByteString)
import qualified Data.ByteString.Lazy as ByteString

-- cassava
import Data.Csv
  ( DefaultOrdered(headerOrder)
  , FromField(parseField)
  , FromNamedRecord(parseNamedRecord)
  , Header
  , ToField(toField)
  , ToNamedRecord(toNamedRecord)
  , (.:)
  , (.=)
  )
import qualified Data.Csv as Cassava

-- text
import Data.Text (Text)
import qualified Data.Text.Encoding as Text

-- Exception
import Control.Exception

-- base
import qualified Control.Monad as Monad
import qualified System.Exit as Exit
import qualified Debug.Trace as Trace

data ComponentData =
  ComponentData {
    componentType :: ComponentType,
    nodeK :: Int,
    nodeM :: Int,
    magnitude :: Double,
    param1 :: Double,
    param2 :: Double,
    plot :: Int
     }
  deriving (Eq, Show)

data SimulationData =
  SimulationData {
    nodes :: Int,
    voltageSources :: Int,
    stepSize :: Double,
    tmax :: Double
     }
  deriving (Eq, Show)

data ComponentType = Resistor | Capacitor | Inductor | EAC | EDC | Other Text deriving (Eq, Show)
type SimulationResults = (Vector Double, Matrix Double)

instance FromNamedRecord SimulationData where
  parseNamedRecord m =
    SimulationData
      <$> m .: "Number of Nodes"
      <*> m .: "Number of Voltages Sources"
      <*> m .: "Step Size"
      <*> m .: "Maximum time for simulation"

instance FromNamedRecord ComponentData where
  parseNamedRecord m =
    ComponentData
      <$> m .: "Element Type"
      <*> m .: "Node K"
      <*> m .: "Node M"
      <*> m .: "Value"
      <*> m .: "Source param 1"
      <*> m .: "Source param 2"
      <*> m .: "Plot"

instance FromField ComponentType where
  parseField "R" =
    pure Resistor

  parseField "L" =
    pure Inductor

  parseField "C" =
    pure Capacitor

  parseField "EDC" =
    pure EDC

  parseField "EAC" =
    pure EAC

  parseField otherType =
    Other <$> parseField otherType

decodeItems :: ByteString -> Either String (Vector ComponentData)
decodeItems =
  fmap snd . Cassava.decodeByName

decodeItemsFromFile :: FilePath -> IO (Either String (Vector ComponentData))
decodeItemsFromFile filePath =
  catchShowIO (ByteString.readFile filePath)
    >>= return . either Left decodeItems

decodeSimulation :: ByteString -> Either String (Vector SimulationData)
decodeSimulation =
  fmap snd . Cassava.decodeByName

decodeSimulationFromFile :: FilePath -> IO (Either String (Vector SimulationData))
decodeSimulationFromFile filePath =
  catchShowIO (ByteString.readFile filePath)
    >>= return . either Left decodeSimulation

getSingleSimulationLine :: Vector SimulationData -> SimulationData
getSingleSimulationLine =
    Vector.head

nhComponents :: [ComponentData] -> [ComponentData]
nhComponents =
  filter (\r -> (componentType r == Capacitor) || (componentType r == Inductor))

filterEnergyStorageComponent :: Vector ComponentData -> Vector ComponentData
filterEnergyStorageComponent =
  Vector.filter (\r -> (componentType r == Capacitor) || (componentType r == Inductor))

nh :: Vector ComponentData -> Int
nh components =
  length $ filterEnergyStorageComponent components

-- filter ((== Zaal) . reviewLocation) reviews

filterSources :: Vector ComponentData -> Vector ComponentData
filterSources =
  Vector.filter (\r -> (componentType r == EDC) || (componentType r == EAC))

condutance :: ComponentData -> Double -> Double
condutance component dt =
  case componentType component of
    Resistor -> 1.0 / (magnitude component)
    Capacitor -> (magnitude component) * 0.000001 * 2 / dt
    Inductor -> dt / (2 * 0.001 * (magnitude component))
    _ -> 0.0

gkm :: Vector ComponentData -> Double -> Vector Double
gkm components dt =
  Vector.map (\c -> condutance c dt) components
  -- Matrix.colVector (Vector.map (\c -> condutance c dt) components)
  -- Matrix.fromList (length components) 1 (Vector.toList (Vector.map (\c -> condutance c dt) components))

buildCompactGMatrix :: Double -> [ComponentData] -> Matrix Double -> Matrix Double
buildCompactGMatrix dt [] buffer = buffer
buildCompactGMatrix dt (component:cs) buffer =
  case (nodeK component, nodeM component) of (0, m) -> buildCompactGMatrix dt cs (Matrix.setElem (Matrix.getElem m m buffer + condutance component dt) (m, m) buffer)
                                             (k, 0) -> buildCompactGMatrix dt cs (Matrix.setElem (Matrix.getElem k k buffer + condutance component dt) (k, k) buffer)
                                             (k, m) -> buildCompactGMatrix dt cs (Matrix.setElem (Matrix.getElem k k buffer + condutance component dt) (k, k) (Matrix.setElem (Matrix.getElem m m buffer + condutance component dt) (m, m) (Matrix.setElem (Matrix.getElem k m buffer - condutance component dt) (k, m) (Matrix.setElem (Matrix.getElem m k buffer - condutance component dt) (m, k) buffer))))
                                             (_, _) -> buildCompactGMatrix dt cs buffer

buildGMatrixFromVector :: SimulationData -> Vector ComponentData -> Matrix Double
buildGMatrixFromVector simulation components =
  buildCompactGMatrix (stepSize simulation) (Vector.toList components) (Matrix.zero (nodes simulation) (nodes simulation))

buildIhVector :: [ComponentData] -> Double -> Int -> [Double] -> [Double] -> Matrix Double -> Vector Double
buildIhVector [] _ _ _ ihnew _ = Vector.fromList ihnew
buildIhVector (component:cs) dt n (hold:ihold) ihnew vMatrix =
  case (componentType component, nodeK component, nodeM component) of (Inductor, 0, m) -> buildIhVector cs dt n ihold (ihnew ++ [(2*(condutance component dt)*(Matrix.getElem m n vMatrix) + hold)]) vMatrix
                                                                      (Inductor, k, 0) -> buildIhVector cs dt n ihold (ihnew ++ [(-2*(condutance component dt)*(Matrix.getElem k n vMatrix) + hold)]) vMatrix
                                                                      (Inductor, k, m) -> buildIhVector cs dt n ihold (ihnew ++ [(-2*(condutance component dt)*((Matrix.getElem k n vMatrix) - (Matrix.getElem m n vMatrix)) + hold)]) vMatrix
                                                                      (Capacitor, 0, m) -> buildIhVector cs dt n ihold (ihnew ++ [(-2*(condutance component dt)*(Matrix.getElem m n vMatrix) - hold)]) vMatrix
                                                                      (Capacitor, k, 0) -> buildIhVector cs dt n ihold (ihnew ++ [(2*(condutance component dt)*(Matrix.getElem k n vMatrix) - hold)]) vMatrix
                                                                      (Capacitor, k, m) -> buildIhVector cs dt n ihold (ihnew ++ [(2*(condutance component dt)*((Matrix.getElem k n vMatrix) - (Matrix.getElem m n vMatrix)) - hold)]) vMatrix
                                                                      (_, _, _) -> buildIhVector cs dt n ihold ihnew vMatrix

buildVBVector :: [ComponentData] -> Double -> [Double] -> Vector Double
buildVBVector [] _ buffer = Vector.fromList buffer
buildVBVector (c:components) time buffer =
  case (componentType c) of EDC -> buildVBVector components time ((magnitude c) : buffer)
                            EAC -> buildVBVector components time (((magnitude c * cos (2 * pi * param2 c * time + (param1 c * (pi/180))))) : buffer)
                            _   -> buildVBVector components time buffer

buildIVector :: [ComponentData] -> [Double] -> Vector Double -> Vector Double
buildIVector [] _ iVector = iVector
buildIVector (component:cs) (ihEl:ih) iVector =
  case (componentType component, nodeK component, nodeM component) of (Inductor, k, 0) -> buildIVector cs ih (iVector Vector.// [((k - 1), ((iVector Vector.! (k-1)) + ihEl))])
                                                                      (Inductor, 0, m) -> buildIVector cs ih (iVector Vector.// [((m - 1), ((iVector Vector.! (m-1)) - ihEl))])
                                                                      (Inductor, k, m) -> buildIVector cs ih (iVector Vector.// [((m - 1), ((iVector Vector.! (m-1)) - ihEl)), ((k-1), ((iVector Vector.! (k-1)) + ihEl))])
                                                                      (Capacitor, k, 0) -> buildIVector cs ih (iVector Vector.// [((k - 1), ((iVector Vector.! (k-1)) + ihEl))])
                                                                      (Capacitor, 0, m) -> buildIVector cs ih (iVector Vector.// [((m - 1), ((iVector Vector.! (m-1)) - ihEl))])
                                                                      (Capacitor, k, m) -> buildIVector cs ih (iVector Vector.// [((m - 1), ((iVector Vector.! (m-1)) - ihEl)), ((k-1), ((iVector Vector.! (k-1)) + ihEl))])
                                                                      (_, _, _) -> buildIVector cs ih iVector

thtaControl :: Int -> Double -> Vector Double -> Vector Double -> SimulationData -> (Int, Vector Double, Double)
thtaControl thtactl time ihnew ih simulation
  | thtactl <= 0 = (thtactl, ihnew, (stepSize simulation + time))
  | thtactl < 3 = (thtactl + 1, (Vector.map (\i -> i/2) $ Vector.zipWith (+) ih ihnew), (time + (stepSize simulation/2)))
  | otherwise = (0, ihnew, (stepSize simulation + time))

fromHMatrixTransformer :: HMatrix.Matrix Double -> Matrix Double
fromHMatrixTransformer matrix =
  Matrix.fromLists $ HMatrix.toLists matrix

toHMatrixTransformer :: Matrix Double -> HMatrix.Matrix Double
toHMatrixTransformer matrix =
  HMatrix.fromLists $ Matrix.toLists matrix

fromHMatrixVectorTransformer :: HMatrix.Vector Double -> Vector Double
fromHMatrixVectorTransformer vec =
  Vector.fromList $ HMatrix.toList vec

toHMatrixVectorTransformer :: Vector Double -> HMatrix.Vector Double
toHMatrixVectorTransformer vec =
  HMatrix.fromList $ Vector.toList vec

solver :: HMatrix.Vector Double -> HMatrix.Matrix Double -> HMatrix.Matrix Double -> HMatrix.Matrix Double -> HMatrix.Matrix Double -> HMatrix.Vector Double -> SimulationData -> (Vector Double, Vector Double)
solver iVector gaa gab gba gbb vb simulation =
  let ia = HMatrix.subVector 0 ((nodes simulation) - (voltageSources simulation)) iVector
      rhsa = ia - (gab HMatrix.#> vb)
      va = gaa HMatrix.<\> rhsa
      ib = (gba HMatrix.#> va) + (gbb HMatrix.#> vb)
      iVec = HMatrix.vjoin [ia, ib]
      vVec = HMatrix.vjoin [va, vb]
  in
    ((fromHMatrixVectorTransformer iVec), (fromHMatrixVectorTransformer vVec))

thtaSimulationStep :: [ComponentData] -> Matrix Double -> SimulationData -> Int -> Int -> Double  -> Vector Double -> Matrix Double -> Vector Double -> Vector Double -> SimulationResults
thtaSimulationStep _ _ _ _ 1 _ _ vMatrix _ iVector = (iVector, vMatrix)
thtaSimulationStep components condutances simulation thtactl n time ih vMatrix vbVector iVector =
  let (gaa, gab, gba, gbb) = Matrix.splitBlocks (nodes simulation - voltageSources simulation) (nodes simulation - voltageSources simulation) condutances
      ihBuffer = buildIhVector (nhComponents components) (stepSize simulation) n (Vector.toList ih) [] vMatrix
      (thta, ihThta, timeThta) = thtaControl thtactl time ihBuffer ih simulation
      vbVec = buildVBVector components timeThta []
      iVec = buildIVector (nhComponents components) (Vector.toList ihThta) (Vector.replicate (nodes simulation) 0)
      -- (iVecCalc, vVec) = Trace.trace ("Solver = \n" ++ show (solver (toHMatrixVectorTransformer iVec) (toHMatrixTransformer gaa) (toHMatrixTransformer gab) (toHMatrixTransformer gba) (toHMatrixTransformer gbb) (toHMatrixVectorTransformer vbVec) simulation)) solver (toHMatrixVectorTransformer iVec) (toHMatrixTransformer gaa) (toHMatrixTransformer gab) (toHMatrixTransformer gba) (toHMatrixTransformer gbb) (toHMatrixVectorTransformer vbVec) simulation
      (iVecCalc, vVec) = solver (toHMatrixVectorTransformer iVec) (toHMatrixTransformer gaa) (toHMatrixTransformer gab) (toHMatrixTransformer gba) (toHMatrixTransformer gbb) (toHMatrixVectorTransformer vbVec) simulation
      vMatr = Matrix.mapCol (\r _ -> vVec Vector.! (r - 1)) (n-1) vMatrix
  in
      thtaSimulationStep components condutances simulation thta (n-1) timeThta ihThta vMatr vbVec iVecCalc

thtaSimulation :: Vector ComponentData -> SimulationData -> SimulationResults
thtaSimulation components simulation =
  thtaSimulationStep (Vector.toList components) (buildGMatrixFromVector simulation components) simulation 1 (npoints simulation) 0.0 (Vector.replicate (nh components) 0) (Matrix.zero (nodes simulation) (npoints simulation)) (Vector.replicate (voltageSources simulation) 0) (Vector.replicate (nodes simulation) 0)

npoints :: SimulationData -> Int
npoints sim =
  round ((tmax sim)/(stepSize sim)) + 1

catchShowIO :: IO a -> IO (Either String a)
catchShowIO action =
  fmap Right action
    `catch` handleIOException
  where
    handleIOException :: IOException -> IO (Either String a)
    handleIOException =
      return . Left . show

main :: IO ()
main = do

  eitherSimulation <-
      fmap getSingleSimulationLine
        <$> decodeSimulationFromFile "data/simulation.csv"

  case eitherSimulation of
    Left reason ->
      Exit.die reason

    Right simulation -> do
      components_list <- decodeItemsFromFile "data/components.csv"
      case components_list of
        Left reason -> Exit.die reason
        Right components -> do
          let gmt = buildGMatrixFromVector simulation components
          putStr "GMatrix: \n"
          print (gmt)
          let results = thtaSimulation components simulation
          putStr "Simulation: \n"
          print (results)

max_stars_repo_name: hannelita/thtahs

max_issues_repo_name: hannelita/thtahs

max_forks_repo_name: hannelita/thtahs


text:
module NNSpec where

import Numeric.LinearAlgebra as LA
import Test.Hspec
import NN

spec :: Spec
spec =
  describe "Toy Neural Network" $

--    it "Train Point" $ do
--        let wT1 = (3 >< 2) $ repeat 1.0
--            b1 = (3 >< 1) $ repeat 0.0
--            vT1 = (1 >< 3) $ repeat 1.0
--            learningRate = 0.001
--            res = trainPoint wT1 b1 vT1 learningRate 2 1 10
--        res `shouldBe` (wT1, b1, vT1)

    it "Test" $ do
        let wT = (3 >< 2) $ repeat 1.0
            b = (3 >< 1) $ repeat 0.0
            vT = (1 >< 3) $ repeat 1.0
            learningRate = 0.001
            epochs = 10
            xyPairs = [((2.0, 1.0), 10.0), ((3.0, 3.0), 21.0), ((4.0, 5.0), 32.0), ((6.0, 6.0), 42.0)]
            (wT', b', vT') = trainNeuralNetwork xyPairs wT b vT learningRate epochs
            x = [(1.0,1.0), (2.0,2.0), (3.0,3.0), (5.0,5.0), (10.0,10.0)]
            expected = fmap (snd . fmap (* 7)) x
            res = testNeuralNetwork x wT' b' vT' learningRate
            res' = map ((< 0.1) . abs) $ zipWith (-) res expected
        res' `shouldBe` [True, True, True, True, True]


text:
module LogHarmonicsNormalization where

import Data.Complex
import Data.List as L

main = do
  let delta = 0.00001 :: Double
      maxScale = 16 :: Double
      freq = 5
      xs =
        L.map
          (\x -> cis $ (-pi) * freq * (log x) /  (log maxScale))
          [1 / maxScale,1 / maxScale + delta .. maxScale]
      y =
        (exp (log maxScale :+ (-pi) * freq) -
         exp (((-1) * log maxScale) :+ pi * freq)) / (1 :+ (-pi) * freq / (log maxScale))
  print . polar $ ((delta :+ 0) * (L.sum xs))
  print . polar $ y

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE DeriveFunctor #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}

module Tests.Matrix.Types
    (
      Mat(..)
    , fromMat
    , toMat
    ) where

import Control.Monad (join)
import Control.Applicative ((<$>), (<*>))
import Statistics.Matrix (Matrix(..), fromList)
import Test.QuickCheck
import Tests.Helpers (shrinkFixedList, small)
import qualified Data.Vector.Unboxed as U

data Mat a = Mat { mrows :: Int , mcols :: Int
                 , asList :: [[a]] }
              deriving (Eq, Ord, Show, Functor)

fromMat :: Mat Double -> Matrix
fromMat (Mat r c xs) = fromList r c (concat xs)

toMat :: Matrix -> Mat Double
toMat (Matrix r c _ v) = Mat r c . split . U.toList $ v
  where split xs@(_:_) = let (h,t) = splitAt c xs
                         in h : split t
        split []       = []

instance (Arbitrary a) => Arbitrary (Mat a) where
    arbitrary = small $ join (arbMat <$> arbitrary <*> arbitrary)
    shrink (Mat r c xs) = Mat r c <$> shrinkFixedList (shrinkFixedList shrink) xs

arbMat :: (Arbitrary a) => Positive (Small Int) -> Positive (Small Int)
       -> Gen (Mat a)
arbMat (Positive (Small r)) (Positive (Small c)) =
    Mat r c <$> vectorOf r (vector c)

instance Arbitrary Matrix where
    arbitrary = fromMat <$> arbitrary
    -- shrink    = map fromMat . shrink . toMat

max_stars_repo_name: StefanHubner/statistics

max_issues_repo_name: StefanHubner/statistics

max_forks_repo_name: StefanHubner/statistics


text:
-- |
-- = Nonlinear transient computing
--
-- This module was developed as a part of author's PhD project:
-- https://www.researchgate.net/project/Theory-and-Modeling-of-Complex-Nonlinear-Delay-Dynamics-Applied-to-Neuromorphic-Computing
--

{-# LANGUAGE BangPatterns #-}

module RC.NTC
  ( new
  , learn
  , learnClassifier
  , predict
  , par0
  , NTCParameters (..)
  ) where

import           Numeric.LinearAlgebra
import           System.Random ( StdGen
                               , split
                               )
import qualified Data.List as List
import qualified Data.Vector.Storable as V
import qualified Data.Vector.Storable.Mutable as VM
import           System.IO.Unsafe ( unsafePerformIO )
import qualified Learning
import qualified Numeric.DDE.Model as DDEModel

import           RC.NTC.Types
import qualified RC.NTC.Reservoir as Reservoir
import qualified RC.Helpers as H

-- | An untrained NTC network
new
  :: StdGen
  -> NTCParameters
  -> (Int, Int, Int)
  -- ^ Input dimension, network nodes, and output dimension
  -> NTC
new g par (ind, nodes, out) =
  let iwgen = _inputWeightsGenerator par
      iw = iwgen g (nodes, ind) (_inputWeightsRange par)
      ntc = NTC { _inputWeights = iw
                , _reservoir = _reservoirModel par
                , _outputWeights = Nothing
                , _par = par
                }
  in ntc

-- | Default NTC parameters
par0 :: NTCParameters
par0 = Par
  { _preprocess = id
  , _inputWeightsGenerator = H.randMatrix
  , _postprocess = H.addBiases  -- Usually `id` will work
  , _inputWeightsRange = undefined  -- To be manually set, e.g. (-1, 1)
  , _reservoirModel = Reservoir.genReservoir p
  }
  where
    filt' = DDEModel.BandpassFiltering {
              DDEModel._tau = 7.8125e-3
            , DDEModel._theta = recip 0.34375
            }

    p = DDEModel.RC { DDEModel._filt = filt'
                    , DDEModel._rho = 3.25
                    , DDEModel._fnl = H.hsigmoid (1.09375, 1.5, 0.0)
                    }

-- | Nonlinear transformation performed by an NTC network
forwardPass :: NTC  -- ^ NTC network
            -> Matrix Double  -- ^ Input information
            -> Matrix Double
forwardPass NTC { _par = Par { _preprocess = prep, _postprocess = post }
                , _inputWeights = iw
                , _reservoir = Reservoir res
                } !sample =
  let pipeline = post. res. (iw <>). prep
  in pipeline sample

-- | NTC training: learn the readout weights offline
learn
  :: NTC
  -> Int
  -- ^ Discard the first N points
  -> Matrix Double
  -- ^ Input matrix of features rows and observations columns
  -> Matrix Double
  -- ^ Desired output matrix of observations columns
  -> Either String NTC
learn ntc forgetPts inp out = ntc'
  where
    state' = forwardPass ntc inp ?? (All, Drop forgetPts)
    teacher' = out ?? (All, Drop forgetPts)
    ntc' = case Learning.learn' state' teacher' of
      Nothing -> Left "Cannot create a readout matrix"
      w -> Right $ ntc { _outputWeights = w }

_concatForwardPass :: NTC
                   -> Int
                   -- ^ Number of samples in the list (for
                   -- memory allocation)
                   -> [Matrix Double]
                   -- ^ List of samples
                   -> Matrix Double
_concatForwardPass ntc m (state0_:samples') = unsafePerformIO $ do
      -- Detect the postprocessing output dimension `nodes`.
      -- Alternatively, use Static shapes from hmatrix.
      let state0 = forwardPass ntc state0_
          nodes = rows state0
          state0' = H.flatten' state0

      -- Allocate memory for a mutable vector
      v <- VM.new (nodes * m)

      -- Copy the first computed state0'
      copy state0' v 0 0 nodes

      let foldA _ [] = return ()
          foldA f ((y, i):ys) = do
            let v0 = f y
            copy v0 v 0 (i * nodes) ((i + 1) * nodes)
            foldA f ys

      -- NB: Does compiler know (H.flatten'. H.unflatten') == id?
      -- If not, consider refactoring genReservoir and forwardPass
      foldA (H.flatten'. forwardPass ntc) $ zip samples' [1..m - 1]

      processed' <- V.unsafeFreeze v
      let processed = H.unflatten' nodes processed'
      return processed
  where
    copy !v0 !v !i0 !i !k
      | i < k = do
        VM.unsafeWrite v i (v0 V.! i0)
        copy v0 v (i0 + 1) (i + 1) k
      | otherwise = return ()

-- | NTC training specific to classification task.
-- The readout weights are learned offline.
--
-- Alternatively, one could use `learn` function. However,
-- to make sure the training samples do not mix in the reservoir RNN,
-- a significant padding of zeros would be needed. Although that
-- way directly corresponds to a physical experiment, that would be not
-- memory-efficient on a computer.
learnClassifier
  :: NTC
  -- ^ NTC network
  -> [Int]
  -- ^ Target labels
  -> Int
  -- ^ Number of inputs (for memory allocation)
  -> [Matrix Double]
  -- ^ Training inputs
  -> Either String (Learning.Classifier Int)
learnClassifier ntc labels m samples =
  case Learning.learn' state teacher' of
    Nothing -> Left "Cannot create a readout matrix"
    Just w -> let clf = Learning.winnerTakesAll w klasses. forwardPass ntc
           in Right (Learning.Classifier clf)
  where
    klasses = fromList. List.sort. List.nub $ labels
    klassesNo = V.length klasses
    teacher' = concatHor. map (flip (Learning.teacher klassesNo) 1) $ labels

    -- Alternatively, a streaming interface might be a solution
    -- https://hackage.haskell.org/package/pipes-4.3.7/docs/Pipes-Tutorial.html
    state = _concatForwardPass ntc m samples

-- | Horizontal matrix concatenation, alternative to fromBlocks
concatHor :: Element a => [Matrix a] -> Matrix a
concatHor ms@(m:_) = foldr (|||) (zeroCols (rows m)) ms
-- concatHor = concatMapHor id
{-# SPECIALIZE concatHor :: [Matrix Double] -> Matrix Double #-}

concatMapHor
  :: Element b =>
     (Matrix a -> Matrix b) -> [Matrix a] -> Matrix b
concatMapHor f ms@(m:_) = foldr (\a b -> f a ||| b) (zeroCols (rows m)) ms
{-# SPECIALIZE concatMapHor :: (Matrix Double -> Matrix Double) -> [Matrix Double] -> Matrix Double #-}

-- | Matrix with zero elements
zeroCols :: V.Storable a => Int -> Matrix a
zeroCols rows' = (rows'><0) []
{-# SPECIALIZE zeroCols :: Int -> Matrix Double #-}

-- | Run prediction using a "clean" (uninitialized) reservoir and then
-- forget the reservoir's state.
-- This can be used for both forecasting and classification tasks.
predict :: NTC  -- ^ Trained network
        -> Int  -- ^ Washout (forget) points
        -> Matrix Double
        -- ^ Input matrix where measurements are columns and features are rows
        -> Either String (Matrix Double)
        -- ^ Either error string or predicted output
predict ntc@NTC { _outputWeights = ow
                } forgetPts inp =
  case ow of
    Nothing -> Left "Please train the NTC first"
    Just w -> let y = forwardPass ntc inp
                  y2 = y ?? (All, Drop forgetPts)
                  prediction = w <> y2
              in Right prediction

max_stars_repo_name: masterdezign/rc

max_issues_repo_name: masterdezign/rc

max_forks_repo_name: masterdezign/rc


text:
-- | Tests for Kernel density estimates.
module Tests.KDE (
  tests
  )where

import Data.Vector.Unboxed             ((!))
import Numeric.Sum                     (kbn, sumVector)
import Statistics.Sample.KernelDensity
import Test.Tasty                      (TestTree, testGroup)
import Test.Tasty.QuickCheck           (testProperty)
import Test.QuickCheck                 (Property, (==>), counterexample)
import Text.Printf                     (printf)
import qualified Data.Vector.Unboxed as U

tests :: TestTree
tests = testGroup "KDE"
  [ testProperty "integral(PDF) == 1" t_densityIsPDF
  ]

t_densityIsPDF :: [Double] -> Property
t_densityIsPDF vec
  = not (null vec) ==> test
  where
    (xs,ys)  = kde 4096 (U.fromList vec)
    step     = (xs ! 1) - (xs ! 0)
    integral = integratePDF step ys
    --
    test = counterexample (printf "Integral %f" integral)
         $ abs (1 - integral) <= 1e-3

integratePDF :: Double -> U.Vector Double -> Double
integratePDF step vec
  = step * sumVector kbn (U.zipWith (*) vec weights)
  where
    n       = U.length vec
    weights = U.generate n go
      where
        go i | i == 0    = 0.5
             | i == n-1  = 0.5
             | otherwise = 1


text:
{-# LANGUAGE ScopedTypeVariables #-}

module Isumi.Math.FFT.SpecUtil where

import Test.QuickCheck
import Data.Complex
import qualified Data.Vector.Unboxed as UV

newtype VectorOf2ExpSize = VectorOf2ExpSize (UV.Vector Double)

instance Show VectorOf2ExpSize where
  show (VectorOf2ExpSize v) = show v

instance Arbitrary VectorOf2ExpSize where
  arbitrary = do
    e :: Int <- choose (0, 10)
    let n = 2^e
    xs <- sequenceA (replicate n arbitrary)
    pure . VectorOf2ExpSize . UV.fromList $ xs

realVToComplexV :: UV.Vector Double -> UV.Vector (Complex Double)
realVToComplexV = UV.map (\x -> x :+ 0)

max_stars_repo_name: IsumiF/fft-bruun

max_issues_repo_name: IsumiF/fft-bruun

max_forks_repo_name: IsumiF/fft-bruun


text:
{-# LANGUAGE OverloadedStrings #-}

module Data.Wavelets.Reconstruction where
import Prelude hiding (map,maximum,minimum,init)
import Data.Wavelets
import Numeric.Statistics
import Statistics.Sample
import Data.Vector.Storable -- hiding (map)
-- import Linear

-- |explicitly name a Fractional and Ord instance on the idwt function
fidwt ::  Int -> WaveletFilter Double -> WaveletPacker Double c -> c -> [Double]
fidwt = idwt

-- | reconstruct the time series raw, without worrying about scaling
reconstructTimeSeries :: Int -> WaveletFilter Double -> WaveletPacker Double c -> c -> Vector Double
reconstructTimeSeries i wft wptc c = vRslt
    where rslt = fidwt i wft wptc c
          vRslt = init.fromList $ rslt


text:
module NeuralNet.Matrix (
  sumRows,
  dims,
  dimsEq,
  zero
) where

import Numeric.LinearAlgebra

zero :: Int -> Int -> Matrix Double
zero r c = (r><c) (repeat 0)

-- TODO: Check if a function available for this
sumRows :: Matrix Double -> Matrix Double
sumRows m = col (map sum (toLists m))

-- TODO: Have a look at size
dims :: Matrix a -> (Int, Int)
dims m = (rows m, cols m)

dimsEq :: Matrix a -> Matrix b -> Bool
dimsEq x y = dims x == dims y

max_stars_repo_name: danielholmes/neural-net

max_issues_repo_name: danielholmes/neural-net

max_forks_repo_name: danielholmes/neural-net


text:
module Neural
  ( Layer(Layer)
  , Network
  , evalLayer
  , evalLayerM
  , evalNetwork
  , evalNetworkM
  , evalOutputError
  , evalOutputErrorM
  , evalLayerError
  , evalLayerErrorM
  , evalNetworkError
  , evalNetworkErrorM
  , gradientDescent
  , gradientDescentM
  , sumColumns
  , vecToMat
  )
where

import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Data

-- Dependent types would be great to represent some of these constraints.

-- A layer in the neural network. The weights are represented by a matrix with
-- 1 row for each neuron, and 1 column for each weight. The biases are represented
-- by a vector. The number of biases must match the number of weights.
data Layer = Layer (Matrix R) (Vector R) deriving Show

-- A network is a list of layers, where each layer has the same number of weights as
-- the previous one has neurons. Otherwise, evalNetwork will fail.
type Network = [Layer]

-- A utility function that converts a vector to a row matrix
vecToMat :: Vector R -> Matrix R
vecToMat = fromRows . (:[])

-- A utility function that sums the colums in a matrix to produce a row vector
sumColumns :: Matrix R -> Vector R
sumColumns m = (vector . replicate (rows m) $ 1) <# m

-- A sigmoid function which is at 0.5 for the input 0, and tends to 1 as the input
-- increases, and to 0 as the input decreases.
sigmoid :: Floating f => f -> f
sigmoid z = 1 / (1 + exp (-z))

-- Differential of the sigmoid fuunction
sigmoid' :: Floating f => f -> f
sigmoid' z = sigmoid z * (1 - sigmoid z)

-- Evaluates a layer of the neural network given a set of inputs. The number of inputs
-- must match the number of weights in the input network or this function will fail
-- with an exception. Returns the outputs as well as the weighted inputs.
evalLayer :: Layer -> Vector R -> (Vector R, Vector R)
evalLayer layer inputs = (flatten output, flatten weightedInput)
  where
    (output, weightedInput) = evalLayerM layer (fromRows [inputs])

-- Evaluates a layer of the neural network given a set of inputs. Contrary to evalLayer,
-- the set of inputs is provided as a matrix where each row is a vector of inputs as
-- given in evalLayer.
evalLayerM :: Layer -> Matrix R -> (Matrix R, Matrix R)
evalLayerM (Layer weights biases) inputs = (sigmoid weightedInputs, weightedInputs)
  where
    inputRows = rows inputs
    biasMatrix = fromColumns . replicate inputRows $ biases
    weightedInputs = tr' $ (weights <> tr' inputs) + biasMatrix

-- Feed forward the given function (e.g. evalLayer) over the network
feedForward :: (Layer -> a -> (a, a)) -> a -> Network -> [(a, a)]
feedForward f inputs network = tail . reverse $ feedForward' f inputs network
  where
    feedForward' f inputs network = foldl (\acc@((a,_):_) l -> (f l a) : acc) [(inputs, undefined)] network

-- Evaluates a neural network by folding a set of inputs through each layer.
-- The number of inputs must match the number of weights in the first layer,
-- and the number of elements in the output vector depends on the number of neurons
-- in the final layer of the network. Returns the output and weighted input of each
-- layer in the network.
evalNetwork :: Network -> Vector R -> [(Vector R, Vector R)]
evalNetwork network inputs = feedForward evalLayer inputs network

-- Evaluates a neural network. Contrary to evalNetwork, the input is taken in the same form
-- as evalLayerM, i.e. a matrix where the rows are a single input vector.
evalNetworkM :: Network -> Matrix R -> [(Matrix R, Matrix R)]
evalNetworkM network inputs = feedForward evalLayerM inputs network

-- Evaluates the error in the output layer of a network
evalOutputError :: Vector R -> Vector R -> Vector R -> Vector R -> Vector R
evalOutputError output weightedInput inputs training = (output - training) * sigmoid' weightedInput

-- Evaluates the error in the output layer of a network. Input is a matrix where rows are inputs.
evalOutputErrorM :: Matrix R -> Matrix R -> Matrix R -> Matrix R -> Matrix R
evalOutputErrorM output weightedInput inputs training = (output - training) * sigmoid' weightedInput

-- Evaluates the error in a layer of a network given the weighted cost, and the next layer's weight
-- and error.
evalLayerError :: Vector R -> Matrix R -> Vector R -> Vector R
evalLayerError weightedCost nextLayerWeights nextLayerError =
  ((tr' nextLayerWeights) #> nextLayerError) * (sigmoid' weightedCost)

-- Evaluates the error in a layer like evalLayerError, but takes a matrix of inputs as rows.
evalLayerErrorM :: Matrix R -> Matrix R -> Matrix R -> Matrix R
evalLayerErrorM weightedCost nextLayerWeights nextLayerError =
  (nextLayerError <> nextLayerWeights) * (sigmoid' weightedCost)

-- Evaluates the error for each layer in the network by means of backpropogation.
evalNetworkError :: Network -> [(Vector R, Vector R)] -> Vector R -> Vector R -> [Vector R]
evalNetworkError network networkValues inputs training =
  map flatten $ evalNetworkErrorM network networkValuesM (vecToMat inputs) (vecToMat training)
    where
      networkValuesM = map (\(a,b) -> (vecToMat a, vecToMat b)) networkValues

-- Evaluates the error for each layer in the network by means of backpropogation.
evalNetworkErrorM :: Network -> [(Matrix R, Matrix R)] -> Matrix R -> Matrix R -> [Matrix R]
evalNetworkErrorM network networkValues inputs training = map fst backPropogation
  where
    -- The layers of the network in reverse order
    reverseNetwork = reverse network
    -- The activations of the layers of the network in reverse order
    reverseLayerVals = reverse networkValues
    -- (reverse layers, reverse activations)
    reverseLayerData = zip reverseNetwork reverseLayerVals
    -- The activations of the output layer of the network
    outputLayer = tail reverseLayerVals
    -- The weights of the output layer
    Layer outputWeights _ = head reverseNetwork
    -- The activations and weighted inputs of the output layer
    (outputValues, outputLayerWeightedInput) = head reverseLayerVals
    -- The error of the output layer
    outputError = evalOutputErrorM outputValues outputLayerWeightedInput inputs training
    -- An accumulator for backpropogation, started with the output layer's error
    acc = [(outputError, outputWeights)]
    -- The backpropogation algorithhm (apply evalOutputError to each layer in reverse order)
    backPropogate = \acc@((loe,lw):_) (Layer w _, (_,z)) -> (evalLayerErrorM z lw loe, w):acc
    backPropogation = foldl backPropogate acc (tail reverseLayerData)

-- Gradient descent function. Given a network, its inputs, its activations, and its error,
-- produces a new network corrected through gradient descent.
gradientDescent :: Network -> Vector R -> [Vector R] -> [Vector R] -> R -> Network
gradientDescent network input activations error learnRate =
  gradientDescentM network (vecToMat input) (map vecToMat activations) (map vecToMat error) learnRate

-- Gradient descent function. Given a network, its inputs, its activations, and its error,
-- produces a new network corrected through gradient descent.
gradientDescentM :: Network -> Matrix R -> [Matrix R] -> [Matrix R] -> R -> Network
gradientDescentM network input activations error learnRate =
  map gradientDescent' layerData
    where
      -- Get the activations from the previous layer of neurons
      -- by offsetting the activations with the input at the start
      -- and discarding the final element
      previousActivations = input : (init activations)
      -- Layers zipped with their errors and the previous layer's activation
      layerData = zip3 network error previousActivations
      -- Gradient descent for one layer
      gradientDescent' (Layer weight bias, error, previousActivation) =
        -- New layer is evaluated by calculating the gradient for the weight and bias and
        -- subtracting gradient * learnRate from the current value
        Layer
          (weight - (weightGradients * realToFrac learnRate))
          (bias - (biasGradients * realToFrac learnRate))
            where
              -- Calculate gradient for each rate
              -- Calculated using the kronecker product, a generalisation of tensor product to matrices
              --weightGradients = kronecker previousActivation (tr' error)
              --weightGradients = (previousActivation <> (tr' error)) -- / (fromIntegral $ cols previousActivation)
              weightGradients = (tr' error <> previousActivation) / (fromIntegral . cols $ previousActivation)
              -- Sum the bias gradient components for each neuron (the columns)
              biasGradients = sumColumns error / (fromIntegral . size $ bias)

max_stars_repo_name: Catchouli/neuron

max_issues_repo_name: Catchouli/neuron

max_forks_repo_name: Catchouli/neuron