text:
{-# LANGUAGE CPP                   #-}
{-# LANGUAGE UndecidableInstances  #-}
{-# LANGUAGE PolyKinds             #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE ViewPatterns          #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# OPTIONS_GHC -Wno-incomplete-uni-patterns #-}

module Grenade.Recurrent.Layers.LSTM (
    LSTM (..)
  , LSTMWeights (..)
  , randomLSTM
  ) where

import           Control.Monad.Primitive           (PrimBase, PrimState)
import           System.Random.MWC hiding (create)
import GHC.TypeLits

-- import           Data.List ( foldl1' )
import           Data.Proxy
import           Data.Serialize

#if MIN_VERSION_base(4,9,0)
import           Data.Kind (Type)
#endif

import qualified Numeric.LinearAlgebra as LA
import           Numeric.LinearAlgebra.Static

import           Grenade.Core
import           Grenade.Utils.LinearAlgebra
import           Grenade.Recurrent.Core
import           Grenade.Layers.Internal.Update

-- | Long Short Term Memory Recurrent unit
--
--   This is a Peephole formulation, so the recurrent shape is
--   just the cell state, the previous output is not held or used
--   at all.
data LSTM :: Nat -> Nat -> Type where
  LSTM :: ( KnownNat input
          , KnownNat output
          ) => !(LSTMWeights input output) -- Weights
            -> !(LSTMWeights input output) -- Momentums
            -> LSTM input output

data LSTMWeights :: Nat -> Nat -> Type where
  LSTMWeights :: ( KnownNat input
                 , KnownNat output
                 ) => {
                   lstmWf :: !(L output input)  -- Weight Forget     (W_f)
                 , lstmUf :: !(L output output) -- Cell State Forget (U_f)
                 , lstmBf :: !(R output)        -- Bias Forget       (b_f)
                 , lstmWi :: !(L output input)  -- Weight Input      (W_i)
                 , lstmUi :: !(L output output) -- Cell State Input  (U_i)
                 , lstmBi :: !(R output)        -- Bias Input        (b_i)
                 , lstmWo :: !(L output input)  -- Weight Output     (W_o)
                 , lstmUo :: !(L output output) -- Cell State Output (U_o)
                 , lstmBo :: !(R output)        -- Bias Output       (b_o)
                 , lstmWc :: !(L output input)  -- Weight Cell       (W_c)
                 , lstmBc :: !(R output)        -- Bias Cell         (b_c)
                 } -> LSTMWeights input output

instance Show (LSTM i o) where
  show LSTM {} = "LSTM"

instance FoldableGradient (LSTMWeights input output) where
  mapGradient f (LSTMWeights wf uf bf wi ui bi wo uo bo wc bc) =
    LSTMWeights (dmmap f wf) (dmmap f uf) (dvmap f bf) (dmmap f wi) (dmmap f ui) (dvmap f bi) (dmmap f wo) (dmmap f uo) (dvmap f bo) (dmmap f wc) (dvmap f bc)
  squaredSums (LSTMWeights wf uf bf wi ui bi wo uo bo wc bc) =
    [ sumM . squareM $ wf
    , sumM . squareM $ uf
    , sumV . squareV $ bf
    , sumM . squareM $ wi
    , sumM . squareM $ ui
    , sumV . squareV $ bi
    , sumM . squareM $ wo
    , sumM . squareM $ uo
    , sumV . squareV $ bo
    , sumM . squareM $ wc
    , sumV . squareV $ bc
    ]

instance (KnownNat i, KnownNat o) => UpdateLayer (LSTM i o) where
  -- The gradients are the same shape as the weights and momentum
  -- This seems to be a general pattern, maybe it should be enforced.
  type Gradient (LSTM i o) = (LSTMWeights i o)

  -- Run the update function for each group matrix/vector of weights, momentums and gradients.
  -- Hmm, maybe the function should be used instead of passing in the learning parameters.
  runUpdate opt@OptSGD{} (LSTM w m) g =
    let MatrixResultSGD wf wf' = u lstmWf w m g
        MatrixResultSGD uf uf' = u lstmUf w m g
        VectorResultSGD bf bf' = v lstmBf w m g
        MatrixResultSGD wi wi' = u lstmWi w m g
        MatrixResultSGD ui ui' = u lstmUi w m g
        VectorResultSGD bi bi' = v lstmBi w m g
        MatrixResultSGD wo wo' = u lstmWo w m g
        MatrixResultSGD uo uo' = u lstmUo w m g
        VectorResultSGD bo bo' = v lstmBo w m g
        MatrixResultSGD wc wc' = u lstmWc w m g
        VectorResultSGD bc bc' = v lstmBc w m g
    in LSTM (LSTMWeights wf uf bf wi ui bi wo uo bo wc bc) (LSTMWeights wf' uf' bf' wi' ui' bi' wo' uo' bo' wc' bc')
      where
    -- Utility function for updating with the momentum, gradients, and weights.
    u :: forall x ix out. (KnownNat ix, KnownNat out) => (x -> L out ix) -> x -> x -> x -> MatrixResult out ix
    u e (e -> weights) (e -> momentum) (e -> gradient) =
      descendMatrix opt (MatrixValuesSGD weights gradient momentum)

    v :: forall x ix. (KnownNat ix) => (x -> R ix) -> x -> x -> x -> VectorResult ix
    v e (e -> weights) (e -> momentum) (e -> gradient) =
      descendVector opt (VectorValuesSGD weights gradient momentum)
  runUpdate _ d g = runUpdate defOptimizer d g

  -- There's a lot of updates here, so to try and minimise the number of data copies
  -- we'll create a mutable bucket for each.
  -- runUpdates rate lstm gs =
  --   let combinedGradient = foldl1' uu gs
  --   in  runUpdate rate lstm combinedGradient
  --     where
  --   uu :: (KnownNat i, KnownNat o) => LSTMWeights i o -> LSTMWeights i o -> LSTMWeights i o
  --   uu a b =
  --     let wf = u lstmWf a b
  --         uf = u lstmUf a b
  --         bf = v lstmBf a b
  --         wi = u lstmWi a b
  --         ui = u lstmUi a b
  --         bi = v lstmBi a b
  --         wo = u lstmWo a b
  --         uo = u lstmUo a b
  --         bo = v lstmBo a b
  --         wc = u lstmWc a b
  --         bc = v lstmBc a b
  --     in LSTMWeights wf uf bf wi ui bi wo uo bo wc bc
  --   u :: forall x ix out. (KnownNat ix, KnownNat out) => (x -> (L out ix)) -> x -> x -> L out ix
  --   u e (e -> a) (e -> b) = tr $ tr a + tr b

  --   v :: forall x ix. (x -> (R ix)) -> x -> x -> R ix
  --   v e (e -> a) (e -> b) = a + b

instance (KnownNat i, KnownNat o, KnownNat (i*o), KnownNat (o*o)) => RandomLayer (LSTM i o) where
  createRandomWith = randomLSTM

instance (KnownNat i, KnownNat o) => RecurrentUpdateLayer (LSTM i o) where
  -- The recurrent shape is the same size as the output.
  -- It's actually the cell state however, as this is a peephole variety LSTM.
  type RecurrentShape (LSTM i o) = S ('D1 o)

instance (KnownNat i, KnownNat o) => RecurrentLayer (LSTM i o) ('D1 i) ('D1 o) where

  -- The tape stores essentially every variable we calculate,
  -- so we don't have to run any forwards component again.
  type RecTape (LSTM i o) ('D1 i) ('D1 o) = (R o, R i, R o, R o, R o, R o, R o, R o, R o, R o, R o)
  -- Forward propagation for the LSTM layer.
  -- The size of the cell state is also the size of the output.
  runRecurrentForwards (LSTM lw _) (S1D cell) (S1D input) =
    let -- Forget state vector
        f_s = lstmBf lw + lstmWf lw #> input + lstmUf lw #> cell
        f_t = sigmoid f_s
        -- Input state vector
        i_s = lstmBi lw + lstmWi lw #> input + lstmUi lw #> cell
        i_t = sigmoid i_s
        -- Output state vector
        o_s = lstmBo lw + lstmWo lw #> input + lstmUo lw #> cell
        o_t = sigmoid o_s
        -- Cell input state vector
        c_s = lstmBc lw + lstmWc lw #> input
        c_x = tanh c_s
        -- Cell state
        c_t = f_t * cell + i_t * c_x
        -- Output (it's sometimes recommended to use tanh c_t)
        h_t = o_t * c_t
    in ((cell, input, f_s, f_t, i_s, i_t, o_s, o_t, c_s, c_x, c_t), S1D c_t, S1D h_t)

  -- Run a backpropogation step for an LSTM layer.
  -- We're doing all the derivatives by hand here, so one should
  -- be extra careful when changing this.
  --
  -- There's a test version using the AD library without hmatrix in the test
  -- suite. These should match always.
  runRecurrentBackwards (LSTM lw _) (cell, input, f_s, f_t, i_s, i_t, o_s, o_t, c_s, c_x, c_t) (S1D cellGrad) (S1D h_t') =
    let -- Reverse Mode AD Derivitives
        c_t' = h_t' * o_t + cellGrad

        f_t' = c_t' * cell
        f_s' = sigmoid' f_s * f_t'

        o_t' = h_t' * c_t
        o_s' = sigmoid' o_s * o_t'

        i_t' = c_t' * c_x
        i_s' = sigmoid' i_s * i_t'

        c_x' = c_t' * i_t
        c_s' = tanh' c_s * c_x'

        -- The derivatives to pass sideways (recurrent) and downwards
        cell'  = tr (lstmUf lw) #> f_s' + tr (lstmUo lw) #> o_s' + tr (lstmUi lw) #> i_s' + c_t' * f_t
        input' = tr (lstmWf lw) #> f_s' + tr (lstmWo lw) #> o_s' + tr (lstmWi lw) #> i_s' + tr (lstmWc lw) #> c_s'

        -- Calculate the gradient Matricies for the input
        lstmWf' = f_s' `outer` input
        lstmWi' = i_s' `outer` input
        lstmWo' = o_s' `outer` input
        lstmWc' = c_s' `outer` input

        -- Calculate the gradient Matricies for the cell
        lstmUf' = f_s' `outer` cell
        lstmUi' = i_s' `outer` cell
        lstmUo' = o_s' `outer` cell

        -- The biases just get the values, but we'll write it so it's obvious
        lstmBf' = f_s'
        lstmBi' = i_s'
        lstmBo' = o_s'
        lstmBc' = c_s'

        gradients = LSTMWeights lstmWf' lstmUf' lstmBf' lstmWi' lstmUi' lstmBi' lstmWo' lstmUo' lstmBo' lstmWc' lstmBc'
    in  (gradients, S1D cell', S1D input')

-- | Generate an LSTM layer with random Weights
--   one can also just call createRandom from UpdateLayer
--
--   Has forget gate biases set to 1 to encourage early learning.
--
--   https://github.com/karpathy/char-rnn/commit/0dfeaa454e687dd0278f036552ea1e48a0a408c9
--
randomLSTM :: forall m i o. (PrimBase m, KnownNat i, KnownNat o, KnownNat (i*o),KnownNat (o*o))
           => NetworkInitSettings -> Gen (PrimState m) -> m (LSTM i o)
randomLSTM (NetworkInitSettings m HMatrix _) gen = do
  let w = getRandomMatrix i o m gen
  let u = getRandomMatrix i o m gen
  let v = getRandomVector i o m gen

  let w0 = konst 0
      u0 = konst 0
      v0 = konst 0

  LSTM <$> (LSTMWeights <$> w <*> u <*> pure (konst 1) <*> w <*> u <*> v <*> w <*> u <*> v <*> w <*> v)
         <*> pure (LSTMWeights w0 u0 v0 w0 u0 v0 w0 u0 v0 w0 v0)

  where i = natVal (Proxy :: Proxy i)
        o = natVal (Proxy :: Proxy o)
randomLSTM (NetworkInitSettings _ cpu _) _ = error $ "CPU backend " ++ show cpu ++ " not supported by LSTM layer"

-- | Maths
--
-- TODO: Move to not here
--       Optimise backwards derivative
sigmoid :: Floating a => a -> a
sigmoid x = 1 / (1 + exp (-x))

sigmoid' :: Floating a => a -> a
sigmoid' x = logix * (1 - logix)
  where
    logix = sigmoid x

tanh' :: (Floating a) => a -> a
tanh' t = 1 - s ^ (2 :: Int)  where s = tanh t

instance (KnownNat i, KnownNat o) => Serialize (LSTM i o) where
  put (LSTM lw _) = do
      u (lstmWf lw)
      u (lstmUf lw)
      v (lstmBf lw)
      u (lstmWi lw)
      u (lstmUi lw)
      v (lstmBi lw)
      u (lstmWo lw)
      u (lstmUo lw)
      v (lstmBo lw)
      u (lstmWc lw)
      v (lstmBc lw)
    where
      u :: forall a b. (KnownNat a, KnownNat b) => Putter  (L b a)
      u = putListOf put . LA.toList . LA.flatten . extract
      v :: forall a. (KnownNat a) => Putter (R a)
      v = putListOf put . LA.toList . extract

  get = do
      w <- LSTMWeights <$> u <*> u <*> v <*> u <*> u <*> v <*> u <*> u <*> v <*> u <*> v
      return $ LSTM w (LSTMWeights w0 u0 v0 w0 u0 v0 w0 u0 v0 w0 v0)
    where
      u :: forall a b. (KnownNat a, KnownNat b) => Get  (L b a)
      u = let f = fromIntegral $ natVal (Proxy :: Proxy a)
          in  maybe (fail "Vector of incorrect size") return . create . LA.reshape f . LA.fromList =<< getListOf get
      v :: forall a. (KnownNat a) => Get (R a)
      v = maybe (fail "Vector of incorrect size") return . create . LA.fromList =<< getListOf get

      w0 = konst 0
      u0 = konst 0
      v0 = konst 0

-------------------- GNum instances --------------------

instance (KnownNat i, KnownNat o) => GNum (LSTM i o) where
  n |* (LSTM w m) = LSTM (n |* w) (n |* m)
  (LSTM w1 m1) |+ (LSTM w2 m2) = LSTM (w1 |+ w2) (m1 |+ m2)

instance (KnownNat i, KnownNat o) => GNum (LSTMWeights i o) where
  n |* (LSTMWeights wf uf bf wi ui bi wo uo bo wc bc) =
    LSTMWeights
      (fromRational n * wf)
      (fromRational n * uf)
      (fromRational n * bf)
      (fromRational n * wi)
      (fromRational n * ui)
      (fromRational n * bi)
      (fromRational n * wo)
      (fromRational n * uo)
      (fromRational n * bo)
      (fromRational n * wc)
      (fromRational n * bc)
  (LSTMWeights wf1 uf1 bf1 wi1 ui1 bi1 wo1 uo1 bo1 wc1 bc1) |+ (LSTMWeights wf2 uf2 bf2 wi2 ui2 bi2 wo2 uo2 bo2 wc2 bc2) =
    LSTMWeights (wf1 + wf2) (uf1 + uf2) (bf1 + bf2) (wi1 + wi2) (ui1 + ui2) (bi1 + bi2) (wo1 + wo2) (uo1 + uo2) (bo1 + bo2) (wc1 + wc2) (bc1 + bc2)

max_stars_repo_name: schnecki/grenade

max_issues_repo_name: schnecki/grenade

max_forks_repo_name: schnecki/grenade


text:
import Criterion.Main
import Criterion.Types
import Dimacs (parseDIMACSFromFile)
import Internal.Sat (solveDIMACS)
import Sat (solve)
import Statistics.Types
import Test.HUnit

solveFile :: String -> Bool -> IO ()
solveFile fileName isSat = do
  parsed <- parseDIMACSFromFile ("resources/problems/" ++ fileName)
  case parsed of
    Right dimacs -> do
      case solveDIMACS dimacs of
        Just _ -> assertBool "" isSat
        Nothing -> assertBool "" (not isSat)
      return ()
    left -> assertFailure $ "" ++ fileName ++ "cannot be parsed"

main :: IO ()
main =
  Criterion.Main.defaultMainWith
    (defaultConfig {timeLimit = 30})
    [ bgroup
        "sat tests"
        [ benchit "add4.cnf" False,
          benchit "diamond1.cnf" False,
          benchit "diamond2.cnf" False,
          benchit "diamond3.cnf" False,
          benchit "ph2.cnf" False,
          benchit "ph3.cnf" False,
          benchit "prime4.cnf" True,
          benchit "prime9.cnf" True,
          benchit "prime25.cnf" True,
          benchit "prime49.cnf" True,
          benchit "prime121.cnf" True
        ]
    ]
  where
    benchit fileName isSat = bench fileName $ whnfIO $ solveFile fileName isSat

ext: hs

max_stars_repo_name: ysnrkdm/frost

max_issues_repo_name: ysnrkdm/frost

max_forks_repo_name: ysnrkdm/frost


text:
{-# LANGUAGE FlexibleContexts #-}

module Data.AdaptiveCoordinateDescent.Internal
  ( normalize
  , perturb
  , makePerturbed
  , fromVector
  , generate
  , plus
  , minus
  , sphere
  , rosen )
  where

import Control.Monad
import Control.Monad.Primitive
import Control.Monad.State.Strict
import Control.Monad.ST
import Data.Foldable
import Data.Traversable
import qualified Data.Vector.Unboxed as V
import qualified Data.Vector.Unboxed.Mutable as VM
import Numeric.LinearAlgebra hiding ( toList )
import Pipes hiding ( for )
import System.Random.MWC

normalize :: Matrix Double -> Matrix Double
normalize mat = runST $ do
  stds <- VM.replicate columns (0, 0)
  for_ [0..columns-1] $ \column -> do
    for_ [0..rows-1] $ \row -> do
      (avg, std) <- VM.read stds column
      VM.write stds column (avg + mat `atIndex` (row, column),
                            std + (mat `atIndex` (row, column))**2)
    (avg, std) <- VM.read stds column
    VM.write stds column (avg / fromIntegral rows
                         ,sqrt $ (std / fromIntegral rows) - (avg / fromIntegral rows)**2)

  fstds <- V.unsafeFreeze stds

  return $ (rows><columns) $ generate (rows*columns) $ \offset ->
             let row = offset `div` columns
                 column = offset `mod` columns

                 (avg, _var) = fstds V.! column -- variance ignored for now

              in ((mat `atIndex` (row, column)) - avg)
 where
  (rows, columns) = size mat

perturb :: (PrimMonad m, Traversable f)
        => f Double
        -> Gen (PrimState m)
        -> m (f Double)
perturb thing rng = for thing $ \val -> do
  perturbance <- uniformR (-0.1, 0.1) rng
  return $ val + perturbance

makePerturbed :: (PrimMonad m, Traversable f)
              => f Double
              -> Gen (PrimState m)
              -> Int
              -> Producer (f Double) m ()
makePerturbed params rng num_items =
  replicateM_ num_items $ do
    new_values <- lift $ perturb params rng
    yield new_values

fromVector :: (Container Vector a, Traversable f) => Vector a -> f void -> f a
fromVector vec structure =
  flip evalState 0 $ for structure $ \_ -> do
    idx <- get
    put (idx+1)
    return $ vec `atIndex` idx
{-# INLINE fromVector #-}

generate' :: Int -> (Int -> a) -> [a]
generate' 0 _ = []
generate' n fun = fun (n-1):generate' (n-1) fun

generate :: Int -> (Int -> a) -> [a]
generate x fun = reverse $ generate' x fun

plus :: (Num a, Traversable f) => f a -> f a -> f a
plus tr1 tr2 = flip evalState (toList tr2) $ for tr1 $ \item -> do
  (x:rest) <- get
  put rest
  return $ item + x

minus :: (Num a, Traversable f) => f a -> f a -> f a
minus tr1 tr2 = flip evalState (toList tr2) $ for tr1 $ \item -> do
  (x:rest) <- get
  put rest
  return $ item - x

sphere :: [Double] -> Double
sphere [x, y] = sqrt $ x*x + y*y
sphere _ = error "sphere takes two arguments."

rosen :: [Double] -> Double
rosen [x, y] = (1 - x)**2 + 100*(y - x**2)**2
rosen _ = error "rosen takes two arguments."

max_stars_repo_name: Noeda/adaptivecoordinatedescent

max_issues_repo_name: Noeda/adaptivecoordinatedescent

max_forks_repo_name: Noeda/adaptivecoordinatedescent


text:
module Camera where

import Types
import Numeric.Vector
import Numeric.Scalar
import VecUtil (unpack3)

rayFrom :: Camera -> (Double, Double) -> Ray
rayFrom (OrthoCamera p forward up right size) (x, y) = let dir = forward
                                                           x' = right * fromScalar (scalar (x * size))
                                                           y' = up * fromScalar (scalar (y * size))
                                                           origin = p + x' + y'
                                                       in  Ray origin dir
rayFrom (Camera p forward up right fov) (x, y) = let mul = tan $ fov * 0.00872664625 -- Magic number is to convert degrees to radians and divide by two
                                                     (r, u, f) = unpack3 $ normalized $ vec3 (x * mul) (y * mul) 1
                                                     dir = forward * (fromScalar $ scalar f) + right * (fromScalar $ scalar r) + up * (fromScalar $ scalar u)
                                                 in  Ray p dir


text:
module AoC2021.Day07 where

import           Data.Bifunctor               (bimap)
import           Data.List                    (sort)
import qualified Numeric.LinearAlgebra.Static as N
import qualified Parsers                      as P
import           Solver                       (Solver)

mean :: (Fractional n, Foldable f) => f n -> n
mean ns
  | null ns = 0
  | otherwise = uncurry (/) . foldl f (0,0) $ ns
  where f (t,c) n = (t+n,c+1)

median :: Ord a => [a] -> a
median as = let m = length as `div` 2
            in  sort as!!m

minFuelCost :: [Int] -> Int
minFuelCost ns = sum . map (abs . (median ns-)) $ ns

day07 :: ([Int] -> Int) -> Solver
day07 algo input = case P.runParser P.parseIntegerList () "" input of
  Left e          -> show e
  Right positions -> show . algo $ positions

solve01 :: Solver
solve01 = day07 minFuelCost


text:
{-# LANGUAGE BangPatterns #-}
module ML.NN.Training
       (
           -- * Training networks
           TrainingConfig(..)
       ,   Sample(..)
       ,   Gradient(..)
       ,   sgd
       ,   gradientDescentCore

           -- * Internal (exposed for testing)
       ,   computeZsAndAs
       ) where

import           ML.CostFunctions (CostDerivative, mse')
import           ML.NN
import           ML.NN.ActivationFunction (ActivationFunctionDerivative)
import           ML.Sample (Sample(..))

import           Control.Monad.Random (evalRandIO)
import           Control.Monad.Random.Extended (vshuffle)
import           Control.Monad.State (State, evalState, get, put)
import           Data.Foldable (foldl')
import           Data.List (unfoldr)
import qualified Data.Vector as V
import           Numeric.LinearAlgebra (Matrix, R, Vector, (#>), (<>), asColumn, asRow, cmap, konst, outer, size, tr)

-- | Configuration for training a network.
data TrainingConfig = TrainingConfig
                      { trainingEta :: R  -- ^ η - Learning rate.
                      , trainingActivation :: ActivationFunction
                      , trainingActivationDerivative :: ActivationFunctionDerivative
                      , trainingCostDerivative :: CostDerivative
                      }

-- | The gradient of the network.
data Gradient = Gradient
                { gradientNablaB :: {-# UNPACK #-} !(V.Vector (Vector R))
                , gradientNablaW :: {-# UNPACK #-} !(V.Vector (Matrix R))
                } deriving (Read, Show)

sumGradients :: V.Vector Gradient -> Gradient
sumGradients = V.foldl1' addGradient
    where (Gradient !lb !lw) `addGradient` (Gradient !rb !rw) =
              Gradient (V.zipWith (+) lb rb) (V.zipWith (+) lw rw)

-- | Return z and activation values for each layer in the network.
--
-- The values are returned in reverse order for use by the
-- backpropogation algorithm.  We the State monad so it's more
-- explicit that the output activation of one layer is the input to
-- the next.
computeZsAndAs :: ActivationFunction -> Network -> Vector R -> ([Vector R], [Vector R])
computeZsAndAs act (Network layers) input = (reverse zs, reverse (input:as))
    where (zs, as) = unzip $ evalState (mapM computeZA layers) input

          computeZA :: Layer -> State (Vector R) (Vector R, Vector R)
          computeZA (Layer b w) = do
              x <- get
              let z = w #> x + b
                  a = cmap act z
              put a
              pure $ (z,a)

-- | Compute the gradient for the training example.
backpropogation :: ActivationFunction
                -> ActivationFunctionDerivative
                -> CostDerivative
                -> Network
                -> Sample
                -> Gradient
backpropogation act act' cost' !net@(Network !layers) (Sample !x !y) =
    Gradient (V.reverse $ V.cons nablaB_L nablaBs) (V.reverse $ V.cons nablaW_L nablaWs)
    where (z:zs, a:a':as) = computeZsAndAs act net x
          delta_L  = cost' a y * cmap act' z
          nablaB_L = delta_L
          nablaW_L = asColumn delta_L <> asRow a'
          (nablaBs, nablaWs) =
              backpropogationBackwardsPass act' delta_L zs as (reverse $ tail layers)

-- | Perform the backwards pass of the backpropogation algorithm.
--
-- Traverse the previously reverse z values, activations, and layers
-- from layer L to 2.  We thread the delta value along as state
backpropogationBackwardsPass :: ActivationFunctionDerivative
                             -> Vector R    -- ^ delta for the last layer.
                             -> [Vector R]  -- ^ zs in reverse order.
                             -> [Vector R]  -- ^ activations in reverse order.
                             -> [Layer]     -- ^ Layers L to 2.
                             -> (V.Vector (Vector R), V.Vector (Matrix R))  -- ^
backpropogationBackwardsPass act' delta_L zs as layers = V.unzip output
    where !output = V.fromList $ evalState (mapM processLayer (zip3 layers zs as)) delta_L

          processLayer :: (Layer, Vector R, Vector R) -> State (Vector R) (Vector R, Matrix R)
          processLayer ((Layer _bias weights), z, a) = do
              delta <- get
              let actPrime = cmap act' z
                  delta'   = (tr weights #> delta) * actPrime
              put delta'
              pure (delta', delta' `outer` a)

-- | Update the network's weights and biases by applying gradient
--   descent for the given sample input.
gradientDescentCore :: TrainingConfig -> V.Vector Sample -> Network -> Network
gradientDescentCore (TrainingConfig eta act act' cost') !trainingData !net@(Network !layers) =
    Network layers'
    where sampleGradients :: V.Vector Gradient
          !sampleGradients = V.map (backpropogation act act' cost' net) trainingData

          Gradient !nablaB !nablaW = sumGradients sampleGradients

          !layers' = V.toList $ V.zipWith3 updateWeightsAndBiases (V.fromList layers) nablaB nablaW

          updateWeightsAndBiases :: Layer -> Vector R -> Matrix R -> Layer
          updateWeightsAndBiases (Layer b w) nb nw =
              Layer (b-(konst (eta/numSamples) (size b)) * nb)
                    (w-(konst (eta/numSamples) (size w)) * nw)
              where numSamples = fromIntegral $ length trainingData

-- | Train the neural network using stochastic gradient descent.
sgd :: TrainingConfig
    -> Int  -- ^ epochs
    -> Int  -- ^ mini-batch size
    -> V.Vector Sample
    -> Network
    -> IO Network
sgd trainingConfig !epochs !miniBatchSize !trainingData !network = trainEpoch 0 network
    where trainEpoch epoch !net
              | epoch == epochs = putStrLn "done!" >> pure net
              | otherwise = do
                    !shuffledTrainingData <- evalRandIO $ vshuffle trainingData
                    let miniBatches = miniBatchSize `vChunksOf` shuffledTrainingData
                        !net' = foldl' (flip (gradientDescentCore trainingConfig)) net miniBatches
                    putStrLn $ "Completed epoch: " ++ show epoch
                    trainEpoch (epoch+1) net'

vChunksOf :: Int -> V.Vector a -> [V.Vector a]
vChunksOf n vec = Data.List.unfoldr makeChunk vec
    where makeChunk :: V.Vector a -> Maybe (V.Vector a, V.Vector a)
          makeChunk v | V.null v  = Nothing
                      | otherwise = Just $ V.splitAt n v

max_stars_repo_name: m-renaud/ML

max_issues_repo_name: m-renaud/ML

max_forks_repo_name: m-renaud/ML


text:
module Statistics.Classification (
    module Classification
  ) where

import Statistics.Classification.ConfusionMatrix as Classification
import Statistics.Classification.ROC as Classification
import Statistics.Classification.Types as Classification

max_stars_repo_name: tsbattman/rochs

max_issues_repo_name: tsbattman/rochs

max_forks_repo_name: tsbattman/rochs


text:
import Test.HUnit
import FFT.FFTData
import FFT.Orig
import FFT.Samples
import System.Environment
import Data.Complex

testfftbFlySParZip = TestCase (assertEqual "FFT-blfys-parzip" expected actual)
    where
        expected = (sum $ fftbFlySParZip samples512)
        actual = (sum $ fft samples512)
        samples512 = samples 1 512

testfftbFlySParZipMap = TestCase (assertEqual "FFT-blfys-parzipmap" expected actual)
    where
        expected = (sum $ fftbFlySParZipMap samples512)
        actual = (sum $ fft samples512)
        samples512 = samples 1 512

testfftMapReduce = TestCase (assertEqual "FFT-map-reduce" expected actual)
    where
        expected = (sum $ fftMapReduce samples512)
        actual = (sum $ fft samples512)
        samples512 = samples 1 512

tests = TestList [TestLabel "testfftFlySParZip" testfftbFlySParZip,
                TestLabel "testfftbFlySParZipMap" testfftbFlySParZipMap,
		TestLabel "testfftMapReduce" testfftMapReduce]

main = runTestTT tests

max_stars_repo_name: BNJHope/parallel-fourier-transforms

max_issues_repo_name: BNJHope/parallel-fourier-transforms

max_forks_repo_name: BNJHope/parallel-fourier-transforms


text:
{-# LANGUAGE FlexibleContexts, UndecidableInstances #-}

module School.Unit.UnitActivation
( UnitActivation(..)
, isApplyFail
) where

import Numeric.LinearAlgebra (Container, Matrix, Vector)

data UnitActivation a =
    BatchActivation (Matrix a)
 |  ApplyFail String deriving (Show)

instance (Container Vector a, Eq a, Num a) => Eq (UnitActivation a) where
  (BatchActivation a1) == (BatchActivation a2) =
     a1 == a2
  (ApplyFail _) == (ApplyFail _) = True
  _ == _ = False

isApplyFail :: UnitActivation a -> Bool
isApplyFail (ApplyFail _) = True
isApplyFail _ = False

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


text:
module EnvelopeTest where

import           Control.Monad                  as M
import           Control.Monad.Parallel         as MP
import qualified Data.Array.Accelerate          as A
import           Data.Array.Accelerate.LLVM.PTX
import           Data.Array.Repa                as R
import           Data.Binary
import           Data.Complex
import           Data.List                      as L
import           Data.Vector.Storable           as VS
import           Data.Vector.Unboxed            as VU
import           Filter.Utils
import           FokkerPlanck.FourierSeries
import           Foreign.CUDA.Driver            as CUDA
import           FourierMethod.BlockMatrixAcc
import           FourierMethod.FourierSeries2D
import           Image.IO
import           Pinwheel.FourierSeries2D
import           STC
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Utils.Array
import           Utils.List
import           Utils.Time

main = do
  args@(deviceIDsStr:numPointsStr:deltaStr:numPointsReconStr:deltaReconStr:numR2FreqStr:periodR2Str:phiFreqsStr:rhoFreqsStr:stdR2Str:stdStr:numBatchR2Str:numBatchR2FreqsStr:sStr:numThreadStr:_) <-
    getArgs
  let deviceIDs = read deviceIDsStr :: [Int]
      numPoints = read numPointsStr :: Int
      delta = read deltaStr :: Double
      numPointsRecon = read numPointsReconStr :: Int
      deltaRecon = read deltaReconStr :: Double
      numR2Freq = read numR2FreqStr :: Int
      periodR2 = read periodR2Str :: Double
      phiFreq = read phiFreqsStr :: Int
      rhoFreq = read rhoFreqsStr :: Int
      numThread = read numThreadStr :: Int
      stdR2 = read stdR2Str :: Double
      std = read stdStr :: Double
      numBatchR2 = read numBatchR2Str :: Int
      numBatchR2Freqs = read numBatchR2FreqsStr :: Int
      s = read sStr :: Double
      folderPath = "output/test/EnvelopeTest"
      func a =
        R.fromFunction (Z :. (1 :: Int) :. numPointsRecon :. numPointsRecon) $ \(Z :. _ :. i :. j) ->
          let r =
                sqrt
                  (fromIntegral $
                   (i - div numPointsRecon 2) ^ 2 +
                   (j - div numPointsRecon 2) ^ 2) :: Double
          in if r == 0
               then 0
               else r ** a :+ 0
  createDirectoryIfMissing True folderPath
  plotImageRepaComplex (folderPath </> "Reference(-0.5).png") .
    ImageRepa 8 . computeS . func $
    (-0.5)
  plotImageRepaComplex (folderPath </> "Reference(-1).png") .
    ImageRepa 8 . computeS . func $
    (-1)
  plotImageRepaComplex (folderPath </> "Reference(-1.5).png") .
    ImageRepa 8 . computeS . func $
    (-1.5)
  initialise []
  devs <- M.mapM device deviceIDs
  ctxs <- M.mapM (\dev -> CUDA.create dev []) devs
  ptxs <- M.mapM createTargetFromContext ctxs
  let envelope' =
        VU.convert . toUnboxed . computeS $
        analyticalFourierCoefficients2 numR2Freq 1 phiFreq rhoFreq (-s) periodR2
      envelope =
        VS.map
          (\x -> x - VS.sum envelope' / (fromIntegral $ numR2Freq ^ 2))
          envelope'
      envelopeMat =
        A.use .
        A.fromList (A.Z A.:. (numR2Freq ^ 2) A.:. (1 :: Int)) . VS.toList $
        envelope
  plotImageRepaComplex (folderPath </> "EnvelopeFreq.png") .
    ImageRepa 8 .
    fromUnboxed (Z :. (1 :: Int) :. numR2Freq :. numR2Freq) . VS.convert $
    envelope
  envelopeR2 <-
    computeFourierSeriesR2StreamAcc
      ptxs
      numR2Freq
      numPointsRecon
      1
      periodR2
      deltaRecon
      numBatchR2
      envelopeMat
  plotImageRepaComplex (folderPath </> "Envelope.png") .
    ImageRepa 8 .
    computeS . extend (Z :. (1 :: Int) :. All :. All) . sumS . rotate3D $
    envelopeR2
      -- envelope2 = VS.map (^ 2) envelope
  let envelope2' =
        VS.map (^ 2) . VU.convert . toUnboxed . computeS $
        analyticalFourierCoefficients2 numR2Freq 1 phiFreq rhoFreq (-s) periodR2
      envelope2 =
        VS.map
          (\x -> x - VS.sum envelope2' / (fromIntegral $ numR2Freq ^ 2))
          envelope2'
      envelope2Mat =
        A.use .
        A.fromList (A.Z A.:. (numR2Freq ^ 2) A.:. (1 :: Int)) . VS.toList $
        envelope2
  plotImageRepaComplex (folderPath </> "Envelope2Freq.png") .
    ImageRepa 8 .
    fromUnboxed (Z :. (1 :: Int) :. numR2Freq :. numR2Freq) . VS.convert $
    envelope2
  envelope2R2 <-
    computeFourierSeriesR2StreamAcc
      ptxs
      numR2Freq
      numPointsRecon
      1
      periodR2
      deltaRecon
      numBatchR2
      envelope2Mat
  plotImageRepaComplex (folderPath </> "Envelope2.png") .
    ImageRepa 8 .
    computeS . extend (Z :. (1 :: Int) :. All :. All) . sumS . rotate3D $
    envelope2R2
  plan <-
    makePlan
      folderPath
      emptyPlan
      numPointsRecon
      numPointsRecon
      numR2Freq
      (2 * phiFreq + 1)
      (2 * rhoFreq + 1)
  let r2Freqs = L.map fromIntegral . getListFromNumber $ numR2Freq
      -- bias =
      --   VU.convert . toUnboxed . computeUnboxedS $
      --   analyticalFourierCoefficients2 numR2Freq 1 0 0 stdR2 periodR2
      bias' =
        computeUnboxedS $
        analyticalFourierCoefficients2 numR2Freq 1 0 0 stdR2 periodR2
      bias =
        VU.convert .
        toUnboxed .
        computeUnboxedS .
        R.map (\x -> x - sumAllS bias' / (fromIntegral $ numR2Freq ^ 2)) $
        bias' :: VS.Vector (Complex Double)
      biasMat =
        A.use .
        A.fromList (A.Z A.:. (numR2Freq ^ 2) A.:. (1 :: Int)) . VS.toList $
        bias
  plotImageRepaComplex (folderPath </> "BiasFreq.png") .
    ImageRepa 8 .
    fromUnboxed (Z :. (1 :: Int) :. numR2Freq :. numR2Freq) . VS.convert $
    bias
  biasR2 <-
    computeFourierSeriesR2StreamAcc
      ptxs
      numR2Freq
      numPointsRecon
      1
      periodR2
      deltaRecon
      numBatchR2
      biasMat
  plotImageRepaComplex (folderPath </> "Bias.png") .
    ImageRepa 8 .
    computeS . extend (Z :. (1 :: Int) :. All :. All) . sumS . rotate3D $
    biasR2
  let forwardPlanID = DFTPlanID DFT1DG [numR2Freq, numR2Freq] [0, 1]
      backwardPlanID = DFTPlanID IDFT1DG [numR2Freq, numR2Freq] [0, 1]
      forwardPlan = getDFTPlan plan forwardPlanID
      backwardPlan = getDFTPlan plan backwardPlanID
      biasVec' =
        VU.convert .
        toUnboxed .
        computeUnboxedS .
        makeFilter2D . fromUnboxed (Z :. numR2Freq :. numR2Freq) . VS.convert $
        bias
      biasVec = VS.map (\x -> x - VS.sum biasVec' / (fromIntegral $ numR2Freq^2)) biasVec'
  biasF <- dftExecuteWithPlan forwardPlanID forwardPlan biasVec
  envelope2F <- dftExecuteWithPlan forwardPlanID forwardPlan envelope2
  biased' <-
    dftExecuteWithPlan backwardPlanID backwardPlan $
    VS.zipWith (*) biasF envelope2F
  let biased =
        VS.map
          (\x -> x - VS.sum biased' / (fromIntegral $ numR2Freq ^ 2))
          biased'
      biasedMat =
        A.use .
        A.fromList (A.Z A.:. (numR2Freq ^ 2) A.:. (1 :: Int)) . VS.toList $
        biased
  plotImageRepaComplex (folderPath </> "BiasedFreq.png") .
    ImageRepa 8 .
    fromUnboxed (Z :. (1 :: Int) :. numR2Freq :. numR2Freq) . VS.convert $
    biased
  biasedR2 <-
    computeFourierSeriesR2StreamAcc
      ptxs
      numR2Freq
      numPointsRecon
      1
      periodR2
      deltaRecon
      numBatchR2
      biasedMat
  plotImageRepaComplex (folderPath </> "Biased.png") .
    ImageRepa 8 .
    computeS . extend (Z :. (1 :: Int) :. All :. All) . sumS . rotate3D $
    biasedR2

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:
import Data.Complex (Complex, realPart)

type CD = Complex Double

quadraticRoots :: (CD, CD, CD) -> (CD, CD)
quadraticRoots (a, b, c) =
  if realPart b > 0
    then ((2 * c) / (-b - d), (-b - d) / (2 * a))
    else ((-b + d) / (2 * a), (2 * c) / (-b + d))
  where
    d = sqrt $ b ^ 2 - 4 * a * c

main :: IO ()
main =
  mapM_
    (print . quadraticRoots)
    [(3, 4, 4 / 3), (3, 2, -1), (3, 2, 1), (1, -10e5, 1), (1, -10e9, 1)]

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
{-# LANGUAGE BangPatterns #-}
module Main where

import qualified Data.Vector                   as V
import qualified Data.Vector.Mutable           as VM
import qualified Data.Vector.Unboxed           as VU
import qualified Data.Vector.Unboxed.Mutable   as VUM
import qualified Data.Vector.Generic           as VG
import qualified Data.Vector.Generic.Mutable   as VGM
import qualified Data.ByteString.Char8         as BS
import qualified Data.ByteString.Lazy.Char8    as BSL
import qualified Data.ByteString.Builder       as BSB
import           Data.Char                      ( isSpace )
import           Data.List                      ( unfoldr )
import           Data.Function
import           Data.Complex
import           Control.Monad
import           Control.Monad.Primitive
import           Control.Monad.ST

main :: IO ()
main = do
  n   <- readLn :: IO Int
  abs <-
    VU.replicateM n
    $   (\v -> (v VU.! 0, v VU.! 1))
    .   VU.unfoldr (BS.readInt . BS.dropWhile isSpace)
    <$> BS.getLine
  let v = uncurry multiplyCT $ VU.unzip abs
  print 0
  VU.mapM_ print v

-- O(n^2) multiply polynomials
naiveMultiply :: VU.Vector Int -> VU.Vector Int -> VU.Vector Int
naiveMultiply f g = VU.create $ do
  let len = VU.length f + VU.length g - 1
  h <- VUM.new len
  flip VU.imapM_ f $ \i x ->
    flip VU.imapM_ g $ \j y -> VUM.modify h (\z -> z + x * y) (i + j)
  return h

-- O(n log n) multiply polynomials by FFT （Immutable Vector)
-- 859 ms in atc001-c
-- c.f. https://qiita.com/ageprocpp/items/0d63d4ed80de4a35fe79
multiply :: VU.Vector Int -> VU.Vector Int -> VU.Vector Int
{-# INLINE multiply #-}
multiply !f !g =
  let !len = VU.length f + VU.length g - 1
      !sz  = fix (\rec !k -> if k >= len then k else rec (k * 2)) 1
      !f'  = VU.map fromIntegral f
        VU.++ VU.replicate (sz - VU.length f) (0 :: Complex Double)
      !g' = VU.map fromIntegral g
        VU.++ VU.replicate (sz - VU.length g) (0 :: Complex Double)
  in  VU.map (round . magnitude . (/ fromIntegral sz))
        . VU.take len
        . dft (-1)
        $ VU.zipWith (*) (dft 1 f') (dft 1 g')

-- DFT algorithm
-- if sign == 1 then DFT
-- if sign == -1 then inverse DFT
dft :: Double -> VU.Vector (Complex Double) -> VU.Vector (Complex Double)
{-# INLINE dft #-}
dft !sign !f
  | VU.length f == 1
  = f
  | otherwise
  = let !n    = VU.length f
        !nd2  = n `div` 2
        !f0   = dft sign $ VU.map (f VU.!) $ VU.enumFromStepN 0 2 nd2
        !f1   = dft sign $ VU.map (f VU.!) $ VU.enumFromStepN 1 2 nd2
        -- f0 = dft sign $ VU.ifilter (\i _ -> i `mod` 2 == 0) f
        -- f1 = dft sign $ VU.ifilter (\i _ -> i `mod` 2 == 1) f
        !zeta = cis (sign * 2 * acos (-1) / fromIntegral n) :: Complex Double
    in  VU.imap (\i z -> (f0 VU.! (i `mod` nd2)) + z * (f1 VU.! (i `mod` nd2)))
          $ VU.iterateN n (* zeta) 1

-- O(n log n) multiply polynomials by FFT
-- c.f. https://qiita.com/ageprocpp/items/0d63d4ed80de4a35fe79
-- multiply :: VU.Vector Int -> VU.Vector Int -> VU.Vector Int
-- multiply f g = runST $ do
--   let len = VU.length f + VU.length g - 1
--   let sz = fix (\rec k -> if k >= len then k else rec (k*2)) 1
--   tf <- VU.thaw (VU.map fromIntegral f :: VU.Vector (Complex Double))
--   nf <- VUM.grow tf (sz - VU.length f)
--   tg <- VU.thaw (VU.map fromIntegral g :: VU.Vector (Complex Double))
--   ng <- VUM.grow tg (sz - VU.length g)
--   dft nf 1 sz
--   dft ng 1 sz
--   forM_ [0..sz-1] $ \i -> do
--     y <- VUM.read ng i
--     VUM.modify nf (* y) i
--   dft nf (-1) sz
--   ret <- VU.freeze nf
--   return $ VU.map ((`div` sz) . round . realPart) $ VU.take len ret

-- DFT algorithm
-- if sign == 1 then DFT
-- if sign == -1 then inverse DFT
-- dft :: PrimMonad m => VUM.MVector (PrimState m) (Complex Double) -> Double -> Int -> m ()
-- {-# INLINE dft #-}
-- dft f sign 1 = return ()
-- dft f sign n = do
--   let nd2 = n `div` 2
--   f0 <- VUM.new nd2
--   f1 <- VUM.new nd2
--   forM_ [0..nd2-1] $ \i -> do
--     x0 <- VUM.read f (2*i)
--     VUM.write f0 i x0
--     x1 <- VUM.read f (2*i+1)
--     VUM.write f1 i x1
--   dft f0 sign nd2
--   dft f1 sign nd2
--   let zeta = cis (sign * 2 * acos (-1) / fromIntegral n)
--   let zetas = VU.iterateN n (*zeta) 1
--   flip VU.imapM_ zetas $ \i z -> do
--     x0 <- VUM.read f0 (i `mod` nd2)
--     x1 <- VUM.read f1 (i `mod` nd2)
--     VUM.write f i (x0 + z * x1)

-- O(n log n) multiply polynomials by Cooley-Tukey FFT （Mutable Vector)
-- 448 ms in atc001-c
-- c.f. http://wwwa.pikara.ne.jp/okojisan/stockham/cooley-tukey.html

-- This function computes polynomial coefficient of (f * g)(x)
-- f & g are polynomial coefficient array of f(x) & g(x)
multiplyCT :: VU.Vector Int -> VU.Vector Int -> VU.Vector Int
multiplyCT f g = runST $ do
  let len = VU.length f + VU.length g - 1
  let sz  = fix (\rec k -> if k >= len then k else rec (k * 2)) 1
  tf <- VU.thaw (VU.map fromIntegral f :: VU.Vector (Complex Double))
  nf <- VUM.grow tf (sz - VU.length f)
  tg <- VU.thaw (VU.map fromIntegral g :: VU.Vector (Complex Double))
  ng <- VUM.grow tg (sz - VU.length g)
  fft sz nf
  fft sz ng
  forM_ [0 .. sz - 1] $ \i -> do
    y <- VUM.read ng i
    VUM.modify nf (* y) i
  ifft sz nf
  ret <- VU.freeze nf
  return $ VU.map (round . realPart) $ VU.take len ret
{-# INLINE multiplyCT #-}

fft
  :: PrimMonad m
  => Int -- length
  -> VUM.MVector (PrimState m) (Complex Double) -- target data of FFT
  -> m ()
fft n = cooleyTukey n 1 0 0
{-# INLINE fft #-}

ifft
  :: PrimMonad m
  => Int -- length
  -> VUM.MVector (PrimState m) (Complex Double) -- target data of FFT
  -> m ()
ifft n x = do
  forM_ [0 .. n - 1] $ VUM.modify x conjugate
  cooleyTukey n 1 0 0 x
  forM_ [0 .. n - 1] $ VUM.modify x ((/ fromIntegral n) . conjugate)
{-# INLINE ifft #-}

cooleyTukey
  :: PrimMonad m
  => Int -- length
  -> Int -- stride
  -> Int -- start position of block
  -> Int -- temp var for bit reverse position in start position of block
  -> VUM.MVector (PrimState m) (Complex Double) -- target data of FFT
  -> m ()
cooleyTukey n s q d x
  | n > 1 = do
    let m  = n `div` 2
    let θ0 = 2 * pi / fromIntegral n
    forM_ [0 .. m - 1] $ \p -> do
      let wp = cos (fromIntegral p * θ0) :+ (-sin (fromIntegral p * θ0))
      a <- VUM.read x (q + p + 0)
      b <- VUM.read x (q + p + m)
      VUM.write x (q + p + 0) (a + b)
      VUM.write x (q + p + m) ((a - b) * wp)
    cooleyTukey m (2 * s) (q + 0) (d + 0) x
    cooleyTukey m (2 * s) (q + m) (d + s) x
  | q > d = VUM.swap x q d
  | otherwise = return ()
{-# INLINE cooleyTukey #-}

max_stars_repo_name: matonix/atcoder

max_issues_repo_name: matonix/atcoder

max_forks_repo_name: matonix/atcoder


text:
{-# LANGUAGE PatternGuards #-}

module Math
  (getInfo, predict, subTrend, randoms)
  where

import Control.Monad
import System.Random hiding (randoms)
import Data.List (transpose,minimumBy)
import Data.Function (on)
import Data.Packed.Matrix
import Numeric.LinearAlgebra.Algorithms

-- import Data.Random.RVar
-- import Data.Random.Distribution.Normal
-- import Data.Random.Source.DevRandom

import Types

solve :: [[Double]] -> [Double] -> [Double]
solve a b = head $ transpose $ toLists $ linearSolve (fromLists a) (fromLists $ map (:[]) b)

type SystemBuilder = Int -> [AnyNumber] -> [AnyNumber] -> ([[AnyNumber]], [AnyNumber])

systemLinear :: SystemBuilder
systemLinear n xs ys = ([[sx2,sx],[sx,fromIntegral n]], [sxy, sy])
  where
    sx = sum xs
    sx2 = sum $ map (^2) xs
    sy = sum ys
    sxy = sum $ zipWith (*) xs ys

systemSquare :: SystemBuilder
systemSquare n xs ys = ([[sx4, sx3, sx2], [sx3, sx2, sx], [sx2, sx, fromIntegral n]],
                        [sxy2, sxy, sy])
  where
    sx = sum xs
    sy = sum ys
    sx2 = sum $ map (^2) xs
    sx3 = sum $ map (^3) xs
    sx4 = sum $ map (^4) xs
    sxy = sum $ zipWith (*) xs ys
    sxy2 = sum $ zipWith (*) (map (^2) xs) ys

systemExponent :: SystemBuilder
systemExponent n xs ys = systemLinear n xs (map log ys)

mls :: SystemBuilder -> Int -> [AnyNumber] -> [AnyNumber] -> (AnyNumber, AnyNumber,AnyNumber)
mls system n xs ys =
  let (ma,mb) = system n xs ys
      ans = solve (map (map toDouble) ma) (map toDouble mb)
      (a,b,c) = case ans of
                  (u:v:[]) -> (u,v,0)
                  (u:v:w:[]) -> (u,v,w)
                  _ -> error "Unexpected number of coefficients!"
  in  (Number a,Number b, Number c)

getInfo :: Int -> Formula -> [AnyNumber] -> [AnyNumber] -> Info
getInfo n f xs ys | Auto <- f =
  let infoL = getInfo' n Linear xs ys
      infoS = getInfo' n Square xs ys
      infoE = getInfo' n Exponent xs ys
  in  minimumBy (compare `on` stdDev) [infoL, infoS, infoE]
                  | otherwise = getInfo' n f xs ys

getInfo' :: Int -> Formula -> [AnyNumber] -> [AnyNumber] -> Info
getInfo' n f xs ys = Info xs ys ts a b c formula
  where
    (a,b,c) = mls system n xs ys
    ts      = map formula xs
    (system, formula) =
      case f of
        Linear -> (systemLinear, linear a b)
        Square -> (systemSquare, square a b c)
        Exponent -> (systemExponent, exponential a b)

avgStep :: Info -> AnyNumber
avgStep info =
  let s = zipWith (-) (tail $ xvals info) (xvals info)
  in  ((sum s)/(fromIntegral $ length s))

stdDev :: Info -> AnyNumber
stdDev info =
  let ts = trend info
      ys = yvals info
      n = length ts
  in  sqrt $ (sum $ map (^2) $ zipWith (-) ys ts)/(fromIntegral n-1)

-- randomsD :: Double -> Int -> IO [Double]
-- randomsD s n = do
--   replicateM n (runRVar (normal 0 s) DevRandom)
--
randoms :: AnyNumber -> Int -> IO [AnyNumber]
randoms s n = do
  lst <- randomsD (toDouble s) n
  return $ map Number lst

randomsD :: Double -> Int -> IO [Double]
randomsD s n = replicateM n $ do
  r1 <- randomIO
  r2 <- randomIO
  return $ s * cos(2*pi*r1)*sqrt(-2*log(r2))

predict :: Int -> Bool -> Info -> IO Info
predict n r info@(Info xs ys ts a b c func) = do
  let st = avgStep info
      xs' = take n $ tail $ iterate (+st) (last xs)
      ys' = map func xs'
  rs <- if r
         then randoms (stdDev info) n
         else return $ replicate n 0
  let ys'' = zipWith (+) ys' rs
  return $ Info (xs++xs') (ys++ys'') (ts++ys') a b c func

subTrend :: Info -> Info
subTrend (Info xs ys ts a b c f) = Info xs ys' ts a b c f
  where
    ys' = zipWith (-) ys ts

linear :: AnyNumber -> AnyNumber -> AnyNumber -> AnyNumber
linear a b x = a*x + b

square :: AnyNumber -> AnyNumber -> AnyNumber -> AnyNumber -> AnyNumber
square a b c x = a*x^2 + b*x + c

exponential :: AnyNumber -> AnyNumber -> AnyNumber -> AnyNumber
exponential a b x = exp (a*x + b)

max_stars_repo_name: portnov/trend

max_issues_repo_name: portnov/trend

max_forks_repo_name: portnov/trend


text:
{-| Filter design and plotting of frequency responses. -}

module SDR.FilterDesign (
    -- * Sinc Function
    sinc,

    -- * Root raised cosine
    srrc,

    -- * Windows
    hanning,
    hamming,
    blackman,

    -- * Convenience Functions
    windowedSinc,

    -- * Frequency Response Plot
    plotFrequency
    ) where

import Graphics.Rendering.Chart.Easy
import Graphics.Rendering.Chart.Backend.Cairo
import Data.Complex

import qualified Data.Vector.Generic as VG

-- | Compute a sinc function
sinc :: (Floating n, VG.Vector v n)
     => Int -- ^ The length. Must be odd.
     -> n   -- ^ The cutoff frequency (from 0 to 1)
     -> v n
sinc size cutoff  = VG.generate size (func . (-) ((size - 1) `quot` 2))
    where
    func 0   = cutoff
    func idx = sin (pi * cutoff * fromIntegral idx) / (fromIntegral idx * pi)

-- | Compute a Hanning window.
hanning :: (Floating n, VG.Vector v n)
        => Int -- ^ The length of the window
        -> v n
hanning size = VG.generate size func
    where
    func idx = 0.5 * (1 - cos((2 * pi * fromIntegral idx) / (fromIntegral size - 1)))

-- | Compute a Hamming window.
hamming :: (Floating n, VG.Vector v n)
        => Int -- ^ The length of the window
        -> v n
hamming size = VG.generate size func
    where
    func idx = 0.54 - 0.46 * cos((2 * pi * fromIntegral idx) / (fromIntegral size - 1))

-- | Compute a Blackman window.
blackman :: (Floating n, VG.Vector v n)
        => Int -- ^ The length of the window
        -> v n
blackman size = VG.generate size func
    where
    func idx = 0.42 - 0.5 * cos((2 * pi * fromIntegral idx) / (fromIntegral size - 1)) + 0.08 * cos((4 * pi * fromIntegral idx) / (fromIntegral size - 1))

-- | Compute a windowed sinc function
windowedSinc :: (Floating n, VG.Vector v n)
             => Int          -- ^ The length
             -> n            -- ^ The cutoff frequency (from 0 to 1)
             -> (Int -> v n) -- ^ The window function
             -> v n
windowedSinc size cutoff window = VG.zipWith (*) (sinc size cutoff) (window size)

signal :: [Double] -> [Double] -> [(Double, Double)]
signal coeffs xs = [ (x / pi, func x) | x <- xs ]
    where
    func phase = magnitude $ sum $ zipWith (\index mag -> mkPolar mag (phase * (- index))) (iterate (+ 1) (- ((fromIntegral (length coeffs) - 1) / 2))) coeffs

-- | Given filter coefficients, plot their frequency response and save the graph as a png file
plotFrequency :: [Double] -- ^ The filter coefficients
              -> FilePath -- ^ The filename
              -> IO ()
plotFrequency coeffs fName = toFile def fName $ do
    layout_title .= "Frequency Response"
    plot (line "Frequency Response" [signal coeffs $ takeWhile (< pi) $ iterate (+ 0.01) 0])

--ts is really the ratio ts / sampling period
-- | Square root raised cosine
srrc :: (Ord a, Floating a)
     => Int -- ^ size: from [-n .. n]
     -> Int -- ^ sampling period
     -> a   -- ^ beta
     -> [a]
srrc n ts beta = map func [(-n) .. n]
    where
    func x
        | x == 0                                                 = 1 - beta + 4 * beta / pi
        | abs (fromIntegral x) ~= (fromIntegral ts / (4 * beta)) = (beta / sqrt 2) * ((1 + 2/pi) * sin (pi / (4 * beta)) + (1 - 2/pi) * cos (pi / (4 * beta)))
        | otherwise                                              = (sin (pi * xdivts * (1 - beta)) + 4 * beta * xdivts * cos (pi * xdivts * (1 + beta))) / (pi * xdivts * (1 - (4 * beta * xdivts) ** 2))
            where
            xdivts = fromIntegral x / fromIntegral ts
    x ~= y = abs (x - y) < 0.001

max_stars_repo_name: adamwalker/sdr

max_issues_repo_name: peixian/sdr

max_forks_repo_name: adamwalker/sdr


text:
-- ch1 ex 1.27

import Statistics.Distribution.Normal

import Statistics.Distribution ( cumulative )

normalD :: NormalDistribution
normalD = normalDistr mean stdD
  where
    mean = 72
    stdD = 8

sampleSize = 340

main = do
  let f low high = (cumulative normalD high) - (cumulative normalD low)
      f1 = f 64 80

  putStrLn $ "The percentage: " ++ show f1
  putStrLn $ "The number of students in the range [64, 80] is approximately: " ++ show (sampleSize * f1)

  let f2 = f 56 88

  putStrLn $ "The percentage: " ++ show f2
  putStrLn $ "The number of students in the range [56, 88] is approximately: " ++ show (sampleSize * f2)

max_stars_repo_name: raydsameshima/Stat

max_issues_repo_name: raydsameshima/Stat

max_forks_repo_name: raydsameshima/Stat


text:
import Control.Applicative
import Data.Generics.Random.Boltzmann.GeneratingFunctions
import Numeric.LinearAlgebra
import Test.HUnit

main = runTestTT $ TestList
  [ test_solve
  ]

test_solve = "F(x) = 1 + x F(x)^2" ~: do
    print xs
    print (evalDeltas <$> xs <*> pure es)
  where
    x : fx : f'x : _ = fmap X [0 ..]
    expectedSize x n = (n * fx, x * f'x)
    es =
      [ expectedSize x 10000
      , (fx, 1 + x * fx * fx)
      , (f'x, fx * fx + 2 * x * f'x * fx)
      ]
    xs = solveEquations defSolveArgs es (vector [0, 1, 1])

max_stars_repo_name: blackheaven/generic-random

max_issues_repo_name: blackheaven/generic-random

max_forks_repo_name: blackheaven/generic-random


text:
module Main where

import           Numeric.LinearAlgebra
import           System.Environment (getArgs)
import           System.Random
import           System.ProgressBar
import           Control.Monad (replicateM)
import           Control.Monad.IO.Class
import           Data.Maybe
import           Data.List
import qualified Data.Vector as V
import           Polynomials
import           Genetic
import           Formulas

main :: IO ()
main = do
    (nGenerations:nPop:goal:nElit:_) <- getArgs

    ------------ Information ------------
    putStrLn ""
    putStrLn $ "No. of generatoins: "  ++ nGenerations
    putStrLn $ "Population Size: "     ++ nPop
    putStrLn $ "Score of goal: "       ++ goal
    putStrLn $ "Elite size (N-elit): " ++ nElit
    putStrLn ""

    ------------ Progress Bar -----------
    putStrLn "Progress:"
    (pb, asy) <- startProgress (msg "") (msg "") 10 (Progress 0 (read nGenerations))

    ------------ Algorithm --------------
    powers <- randomPowerGen 6 2
    let
        ------ Inps and outps of function ------
        inps  = inpsRosen
        outps = fromList $ V.toList yRosen
        ----------------------------------------
        individual = Article powers inps outps
        init = replicate (read nPop) individual

    algo <- runGeneticAlgorithm (read goal) (read nGenerations) (read nElit) init pb
    let
        finalPowers = _articlePowers algo
        x = calcX inps finalPowers
        y = outps
        coeffs = normEq x y
        ----------------------------------------
    grade <- show . fromJust <$> fitness algo

    putStrLn ""
    putStrLn $ showPolynomial coeffs finalPowers
    putStrLn $ "Fitness of final Polynomial is: " ++ grade
    putStrLn ""
    putStrLn "Attention: component with coefficient smaller tan 0.001 are rounded to nothing"

openVLists = V.map V.toList

-- Check x + 1:
simple x = x + 1.0

inpsSimple :: V.Vector Double
inpsSimple = V.fromList
       [ 1
       , 2
       , 3
       , 4]

powsSimple :: V.Vector Double
powsSimple = V.fromList
            [1.0,
             0]

ySimple = V.map simple inpsSimple

-- Check x**2 + y**2 -1:
circle [x, y] = x**2 + y**2 --  -1

inpsCircle :: Inplist
inpsCircle = V.fromList
     . map V.fromList
     $ [ [0.0001, 1.1]
       , [2, 0.001]
       , [2, 1.05]
       , [1, 2.11]
       , [0.1, 2]]

powsCircle :: Powers
powsCircle = V.fromList
          . map V.fromList
          $ [[2, 0],
             [0, 2],
             [0, 0]]

yCircle = V.map circle $ openVLists inpsCircle

-- Check Rosenbrock:
rosen [x, y] =
    1 - 2*x + x**2 + 100*y**2 + 200*x**2*y + 100*x**4

inpsRosen :: Inplist
inpsRosen = V.fromList
     . map V.fromList
     $ [[i, j] | i <- [-0.9, -0.8 .. -0.2], j <- [0.2, 0.3 .. 0.9]]

powsRosen :: Powers
powsRosen = V.fromList
        . map V.fromList
        $ [[0, 0],
           [1, 0],
           [2, 0],
           [0, 2],
           [2, 1],
           [4, 0]]

yRosen = V.map rosen $ openVLists inpsRosen


text:
module Trade.Algorithm.Bollinger where

import qualified Data.Vector as Vec
import Data.Vector (Vector)

import Statistics.Sample (mean, stdDev)

-- import Trade.Render.Svg.Plot

import Debug.Trace

-- number stddev
newtype BolK = BolK Int deriving (Show)

-- rolling window
newtype BolWin = BolWin Int deriving (Show)

meanStDev :: BolWin -> Vector Double -> (Vector (Double, Double), Vector (Double, Double))
meanStDev (BolWin l) v
  | l > Vec.length v = error $ "meanStDev: input vector to short: " ++ show l ++ "/" ++ show (Vec.length v)
meanStDev (BolWin l) v =
  let len = Vec.length v
      f a =
        let u = Vec.slice a l v
        in (mean u, stdDev u)
      idx = Vec.generate (len-l) (fromIntegral . (l+))
      (ms, ss) = Vec.unzip $ Vec.generate (len-l) f
  in (Vec.zip idx ms, Vec.zip idx ss)

data Bollinger = Bollinger {
  n :: BolWin
  , k :: BolK
  , bMean :: Vector (Double, Double)
  , bStDev :: Vector (Double, Double)
  } deriving (Show)

bollinger :: BolK -> BolWin -> Vector Double -> Bollinger
bollinger k n = uncurry (Bollinger n k) . meanStDev n

{-
bollinger2Lines :: String -> BolK -> BolWin -> Vector Double -> [PlotItem Vector Double]
bollinger2Lines str k@(BolK j) n v =
  let Bollinger _ _ ms ss = bollinger k n v
      k' = fromIntegral j
      f j (i, m) (_, s) = (i, m + j*s)
  in [ Line (str ++ " μ") ms
     , Line (str ++ " " ++ show j ++ "σ") (Vec.zipWith (f k') ms ss)
     , Line (str ++ " " ++ show j ++ "σ") (Vec.zipWith (f (-k')) ms ss) ]
-}

max_stars_repo_name: fphh/trade

max_issues_repo_name: fphh/trade

max_forks_repo_name: fphh/trade


text:
{-# LANGUAGE RankNTypes, BangPatterns, GADTs #-}
{-# OPTIONS -Wall #-}

module Language.Hakaru.Distribution where

import System.Random
import Language.Hakaru.Mixture
import Language.Hakaru.Types
import Data.Ix
import Data.Maybe (fromMaybe)
import Data.List (findIndex, foldl')
import Numeric.SpecFunctions
import qualified Data.Map.Strict as M
import qualified Data.Number.LogFloat as LF

mapFst :: (t -> s) -> (t, u) -> (s, u)
mapFst f (a,b) = (f a, b)

dirac :: (Eq a) => a -> Dist a
dirac theta = Dist {logDensity = (\ (Discrete x) -> if x == theta then 0 else log 0),
                    distSample = (\ g -> (Discrete theta,g))}

bern :: Double -> Dist Bool
bern p = Dist {logDensity = (\ (Discrete x) -> log (if x then p else 1 - p)),
               distSample = (\ g -> case randomR (0, 1) g of
                                  (t, g') -> (Discrete $ t <= p, g'))}

uniform :: Double -> Double -> Dist Double
uniform lo hi =
    let uniformLogDensity lo' hi' x | lo' <= x && x <= hi' = log (recip (hi' - lo'))
        uniformLogDensity _ _ _ = log 0
    in Dist {logDensity = (\ (Lebesgue x) -> uniformLogDensity lo hi x),
             distSample = (\ g -> mapFst Lebesgue $ randomR (lo, hi) g)}

uniformD :: (Ix a, Random a) => a -> a -> Dist a
uniformD lo hi =
    let uniformLogDensity lo' hi' x | lo' <= x && x <= hi' = log density
        uniformLogDensity _ _ _ = log 0
        density = recip (fromInteger (toInteger (rangeSize (lo,hi))))
    in Dist {logDensity = (\ (Discrete x) -> uniformLogDensity lo hi x),
             distSample = (\ g -> mapFst Discrete $ randomR (lo, hi) g)}

marsaglia :: (RandomGen g, Random a, Ord a, Floating a) => g -> ((a, a), g)
marsaglia g0 = -- "Marsaglia polar method"
  let (x, g1) = randomR (-1,1) g0
      (y, g ) = randomR (-1,1) g1
      s       = x * x + y * y
      q       = sqrt ((-2) * log s / s)
  in if 1 >= s && s > 0 then ((x * q, y * q), g) else marsaglia g

choose :: (RandomGen g) => Mixture k -> g -> (k, Prob, g)
choose (Mixture m) g0 =
  let peak = maximum (M.elems m)
      unMix = M.map (LF.fromLogFloat . (/peak)) m
      total = M.foldl' (+) (0::Double) unMix
      (p, g) = randomR (0, total) g0
      f !k !v b !p0 = let p1 = p0 + v in if p <= p1 then k else b p1
      err p0 = error ("choose: failure p0=" ++ show p0 ++
                      " total=" ++ show total ++
                      " size=" ++ show (M.size m))
  in (M.foldrWithKey f err unMix 0, LF.logFloat total * peak, g)

chooseIndex :: (RandomGen g) => [Double] -> g -> (Int, g)
chooseIndex probs g0 =
  let (p, g) = random g0
      k = fromMaybe (error ("chooseIndex: failure p=" ++ show p))
                    (findIndex (p <=) (scanl1 (+) probs))
  in (k, g)

normal_rng :: (Real a, Floating a, Random a, RandomGen g) =>
              a -> a -> g -> (a, g)
normal_rng mu sd g | sd > 0 = case marsaglia g of
                                ((x, _), g1) -> (mu + sd * x, g1)
normal_rng _ _ _ = error "normal: invalid parameters"

normalLogDensity :: Floating a => a -> a -> a -> a
normalLogDensity mu sd x = (-tau * square (x - mu)
                            + log (tau / pi / 2)) / 2
  where square y = y * y
        tau = 1 / square sd

normal :: Double -> Double -> Dist Double
normal mu sd = Dist {logDensity = normalLogDensity mu sd . fromLebesgue,
                     distSample = mapFst Lebesgue . normal_rng mu sd}

categoricalLogDensity :: (Eq b, Floating a) => [(b, a)] -> b -> a
categoricalLogDensity list x = log $ fromMaybe 0 (lookup x list)
categoricalSample :: (Num b, Ord b, RandomGen g, Random b) =>
    [(t,b)] -> g -> (t, g)
categoricalSample list g = (elem', g1)
    where
      (p, g1) = randomR (0, total) g
      elem' = fst $ head $ filter (\(_,p0) -> p <= p0) sumList
      sumList = scanl1 (\acc (a, b) -> (a, b + snd(acc))) list
      total = sum $ map snd list

categorical :: Eq a => [(a,Double)] -> Dist a
categorical list = Dist {logDensity = categoricalLogDensity list . fromDiscrete,
                         distSample = mapFst Discrete . categoricalSample list}

lnFact :: Integer -> Double
lnFact = logFactorial

-- Makes use of Atkinson's algorithm as described in:
-- Monte Carlo Statistical Methods pg. 55
--
-- Further discussion at:
-- http://www.johndcook.com/blog/2010/06/14/generating-poisson-random-values/
poisson_rng :: (RandomGen g) => Double -> g -> (Integer, g)
poisson_rng lambda g0 = make_poisson g0
   where smu = sqrt lambda
         b  = 0.931 + 2.53*smu
         a  = -0.059 + 0.02483*b
         vr = 0.9277 - 3.6224/(b - 2)
         arep  = 1.1239 + 1.1368/(b-3.4)
         lnlam = log lambda

         make_poisson :: (RandomGen g) => g -> (Integer,g)
         make_poisson g = let (u, g1) = randomR (-0.5,0.5) g
                              (v, g2) = randomR (0,1) g1
                              us = 0.5 - abs u
                              k = floor $ (2*a / us + b)*u + lambda + 0.43 in
                          case () of
                            () | us >= 0.07 && v <= vr -> (k, g2)
                            () | k < 0 -> make_poisson g2
                            () | us <= 0.013 && v > us -> make_poisson g2
                            () | accept_region us v k -> (k, g2)
                            _  -> make_poisson g2

         accept_region :: Double -> Double -> Integer -> Bool
         accept_region us v k = log (v * arep / (a/(us*us)+b)) <=
                                -lambda + (fromIntegral k)*lnlam - lnFact k

poisson :: Double -> Dist Integer
poisson l =
    let poissonLogDensity l' x | l' > 0 && x> 0 = (fromIntegral x)*(log l') - lnFact x - l'
        poissonLogDensity l' x | x==0 = -l'
        poissonLogDensity _ _ = log 0
    in Dist {logDensity = poissonLogDensity l . fromDiscrete,
             distSample = mapFst Discrete . poisson_rng l}

-- Direct implementation of  "A Simple Method for Generating Gamma Variables"
-- by George Marsaglia and Wai Wan Tsang.
gamma_rng :: (RandomGen g) => Double -> Double -> g -> (Double, g)
gamma_rng shape _   _ | shape <= 0.0  = error "gamma: got a negative shape paramater"
gamma_rng _     scl _ | scl <= 0.0  = error "gamma: got a negative scale paramater"
gamma_rng shape scl g | shape <  1.0  = (gvar2, g2)
                      where (gvar1, g1) = gamma_rng (shape + 1) scl g
                            (w,  g2) = randomR (0,1) g1
                            gvar2 = scl * gvar1 * (w ** recip shape)
gamma_rng shape scl g =
    let d = shape - 1/3
        c = recip $ sqrt $ 9*d
        -- Algorithm recommends inlining normal generator
        n = normal_rng 1 c
        (v, g2) = until (\y -> fst y > 0.0) (\ (_, g') -> normal_rng 1 c g') (n g)
        x = (v - 1) / c
        sqr = x * x
        v3 = v * v * v
        (u, g3) = randomR (0.0, 1.0) g2
        accept  = u < 1.0 - 0.0331*(sqr*sqr) || log u < 0.5*sqr + d*(1.0 - v3 + log v3)
    in case accept of
         True -> (scl*d*v3, g3)
         False -> gamma_rng shape scl g3

gammaLogDensity :: Double -> Double -> Double -> Double
gammaLogDensity shape scl x | x>= 0 && shape > 0 && scl > 0 =
     scl * log shape - scl * x + (shape - 1) * log x - logGamma shape
gammaLogDensity _ _ _ = log 0

gamma :: Double -> Double -> Dist Double
gamma shape scl = Dist {logDensity = gammaLogDensity shape scl . fromLebesgue,
                          distSample = mapFst Lebesgue . gamma_rng shape scl}

beta_rng :: (RandomGen g) => Double -> Double -> g -> (Double, g)
beta_rng a b g | a <= 1.0 && b <= 1.0 =
                 let (u, g1) = randomR (0.0, 1.0) g
                     (v, g2) = randomR (0.0, 1.0) g1
                     x = u ** (recip a)
                     y = v ** (recip b)
                 in  case (x+y) <= 1.0 of
                       True -> (x / (x + y), g2)
                       False -> beta_rng a b g2
beta_rng a b g = let (ga, g1) = gamma_rng a 1 g
                     (gb, g2) = gamma_rng b 1 g1
                 in (ga / (ga + gb), g2)

betaLogDensity :: Double -> Double -> Double -> Double
betaLogDensity _ _ x | x < 0 || x > 1 = error "beta: value must be between 0 and 1"
betaLogDensity a b _ | a <= 0 || b <= 0 = error "beta: parameters must be positve"
betaLogDensity a b x = (logGamma (a + b)
                        - logGamma a
                        - logGamma b
                        + (a - 1) * log x
                        + (b - 1) * log (1 - x))

beta :: Double -> Double -> Dist Double
beta a b = Dist {logDensity = betaLogDensity a b . fromLebesgue,
                 distSample = mapFst Lebesgue . beta_rng a b}

laplace_rng :: (RandomGen g) => Double -> Double -> g -> (Double, g)
laplace_rng mu sd g = sample (randomR (0.0, 1.0) g)
   where sample (u, g1) = case u < 0.5 of
                            True  -> (mu + sd * log (u + u), g1)
                            False -> (mu - sd * log (2.0 - u - u), g1)

laplaceLogDensity :: Floating a => a -> a -> a -> a
laplaceLogDensity mu sd x = - log (2 * sd) - abs (x - mu) / sd

laplace :: Double -> Double -> Dist Double
laplace mu sd = Dist {logDensity = laplaceLogDensity mu sd . fromLebesgue,
                      distSample = mapFst Lebesgue . laplace_rng mu sd}

-- Consider having dirichlet return Vector
-- Note: This is acutally symmetric dirichlet
dirichlet_rng :: (RandomGen g) => Int ->  Double -> g -> ([Double], g)
dirichlet_rng n' a g' = normalize (gammas g' n')
  where gammas g 0 = ([], 0, g)
        gammas g n = let (xs, total, g1) = gammas g (n-1)
                         ( x, g2) = gamma_rng a 1 g1
                     in ((x : xs), x+total, g2)
        normalize (b, total, h) = (map (/ total) b, h)

dirichletLogDensity :: [Double] -> [Double] -> Double
dirichletLogDensity a x | all (> 0) x = sum' (zipWith logTerm a x) + logGamma (sum a)
  where sum' = foldl' (+) 0
        logTerm b y = (b-1) * log y - logGamma b
dirichletLogDensity _ _ = error "dirichlet: all values must be between 0 and 1"

max_stars_repo_name: zaxtax/hakaru-old

max_issues_repo_name: zaxtax/hakaru-old

max_forks_repo_name: zaxtax/hakaru-old


text:
module Tests.Matrix (tests) where

import Statistics.Matrix hiding (map)
import Statistics.Matrix.Algorithms
import Test.Tasty (TestTree, testGroup)
import Test.Tasty.QuickCheck (testProperty)
import Test.QuickCheck
import Tests.ApproxEq (ApproxEq(..))
import Tests.Matrix.Types
import qualified Data.Vector.Unboxed as U

t_row :: Mat Double -> Gen Property
t_row ms@(Mat r _ xs) = do
  i <- choose (0,r-1)
  return $ row (fromMat ms) i === U.fromList (xs !! i)

t_column :: Mat Double -> Gen Property
t_column ms@(Mat _ c xs) = do
  i <- choose (0,c-1)
  return $ column (fromMat ms) i === U.fromList (map (!! i) xs)

t_center :: Mat Double -> Property
t_center ms@(Mat r c xs) =
  (xs !! (r `quot` 2)) !! (c `quot` 2) === center (fromMat ms)

t_transpose :: Matrix -> Property
t_transpose m = U.concat (map (column n) [0..rows m-1]) === toVector m
  where n = transpose m

t_qr :: Matrix -> Property
t_qr a = hasNaN p .||. eql 1e-10 a p
  where p = uncurry multiply (qr a)

tests :: TestTree
tests = testGroup "Matrix"
  [ testProperty "t_row" t_row
  , testProperty "t_column" t_column
  , testProperty "t_center" t_center
  , testProperty "t_transpose" t_transpose
  , testProperty "t_qr" t_qr
  ]


text:
import Control.Monad
import qualified Data.Vector.Storable as VS
import Numeric.LinearAlgebra
    ( (><), ident, Additive(add), Linear(scale) )
import Numeric.Minimization.QuadProgPP ( solveQuadProg )
import Test.Tasty ( defaultMain, testGroup, TestTree )
import Test.Tasty.HUnit ( testCase, assertEqual, assertFailure )
import Test.HUnit.Approx ( assertApproxEqual )
import Debug.Trace
import Numeric.LinearAlgebra.HMatrix (tr)

main :: IO ()
main = defaultMain tests

tests :: TestTree
tests = testGroup "small"
    [ problem0
    , problem1
    , problemQuadprogppDemo
    ]

-- Problem 0
-- minimize: 2(x-1)^2 + (y-4)^2
-- s.t.: x >= 0, y >= 0, x + 2y <= 2
problem0 :: TestTree
problem0 = testCase "Problem 0" $ do
    let answer = solveQuadProg
            ( (2><2)
                [ 4, 0
                , 0, 2
                ]
            , VS.fromList [-4, -8]
            )
            Nothing
            ( Just
                ( (3><2)
                    [ 1, 0
                    , 0, 1
                    ,-1,-2
                    ]
                , VS.fromList [0, 0, 2]
                )
            )
    (actual, best) <- case answer of
        Left err -> assertFailure $ show err
        Right x -> return x
    let expected = VS.fromList [2/9, 8/9] :: VS.Vector Double
    assertEqual "The answer should have the expected length"
        (VS.length expected)
        (VS.length actual)
    VS.zipWithM_ (assertApproxEqual "" 1e-12) (traceShowId expected) (traceShowId actual)

-- Problem 1
-- Same as Problem 0, but substitute:
--  s = x
--  t = x + y
--  u = x - y
-- and add a small positive constant to the diagonal elements to keep
-- the problem positive definite.
problem1 :: TestTree
problem1 = testCase "Problem 1" $ do
    let answer = solveQuadProg
            ((3><3)
                [   1, 2/3, 1/3
                , 2/3, 2/3,   0
                , 1/3,   0, 1/3
                ] `add` (1e-12 `scale` ident 3)
            , VS.fromList [-2, -4, 2]
            )
            (Just
                ( (1><3) [-3, 2, 1]
                , VS.fromList [0]
                )
            )
            (Just
                ((3><3)
                    [ 1,    0,    0
                    , 0,  1/3, -1/3
                    , 0, -4/3,  1/3
                    ]
                , VS.fromList [0, 0, 2]
                )
            )
    (actual, best) <- case answer of
        Left err -> assertFailure $ show err
        Right x -> return x
    let expected = VS.fromList [2/9, 10/9, -14/9] :: VS.Vector Double
    assertEqual "The answer should have the expected length"
        (traceShowId $ VS.length expected)
        (traceShowId $ VS.length actual)
    VS.zipWithM_ (assertApproxEqual "" 1e-5)
        (traceShowId expected)
        (traceShowId actual)

problemQuadprogppDemo :: TestTree
problemQuadprogppDemo = testCase "Problem quadpropp demo" $ do
    let answer = solveQuadProg
            ((2><2)
                [  4, -2
                , -2,  4 ]
            , VS.fromList [6, 0]
            )
            (Just
                ( (1><2) [1, 1]
                , VS.fromList [-3]
                )
            )
            (Just
                ((3><2)
                    [ 1, 0
                    , 0, 1
                    , 1, 1 ]
                , VS.fromList [0, -2, 0]
                )
            )
    (actual, best) <- case answer of
        Left err -> assertFailure $ show err
        Right x -> return x
    let epsilon = 1e-5
    assertApproxEqual "best" epsilon 12 best
    let expected = VS.fromList [1, 2] :: VS.Vector Double
    assertEqual "The answer should have the expected length"
        (traceShowId $ VS.length expected)
        (traceShowId $ VS.length actual)
    VS.zipWithM_ (assertApproxEqual "" epsilon)
        (traceShowId expected)
        (traceShowId actual)

max_forks_repo_name: tsurucapital/hmatrix-quadprogpp


text:
{-# LANGUAGE TypeOperators, CPP #-}
{-# LANGUAGE FlexibleInstances  #-}
{-# LANGUAGE FlexibleContexts   #-}
{-# LANGUAGE DefaultSignatures   #-}
{-# LANGUAGE ScopedTypeVariables #-}
----------------------------------------------------------------------
-- |
-- Module      :   Data.AdditiveGroup
-- Copyright   :  (c) Conal Elliott and Andy J Gill 2008
-- License     :  BSD3
--
-- Maintainer  :  conal@conal.net, andygill@ku.edu
-- Stability   :  experimental
--
-- Groups: zero, addition, and negation (additive inverse)
----------------------------------------------------------------------

module Data.AdditiveGroup
  (
    AdditiveGroup(..), sumV
  , Sum(..), inSum, inSum2
  ) where

import Prelude hiding (foldr)

import Control.Applicative
#if !(MIN_VERSION_base(4,8,0))
import Data.Monoid (Monoid(..))
import Data.Foldable (Foldable)
#endif
import Data.Foldable (foldr)
import Data.Complex hiding (magnitude)
import Data.Ratio
#if !(MIN_VERSION_base(4,11,0))
import Data.Semigroup (Semigroup(..))
#endif
import Foreign.C.Types (CSChar, CInt, CShort, CLong, CLLong, CIntMax, CFloat, CDouble)

import Data.MemoTrie

import Data.VectorSpace.Generic
import qualified GHC.Generics as Gnrx
import GHC.Generics (Generic, (:*:)(..))

infixl 6 ^+^, ^-^

-- | Additive group @v@.
class AdditiveGroup v where
  -- | The zero element: identity for '(^+^)'
  zeroV :: v
  default zeroV :: (Generic v, AdditiveGroup (VRep v)) => v
  zeroV = Gnrx.to (zeroV :: VRep v)
  -- | Add vectors
  (^+^) :: v -> v -> v
  default (^+^) :: (Generic v, AdditiveGroup (VRep v)) => v -> v -> v
  v ^+^ v' = Gnrx.to (Gnrx.from v ^+^ Gnrx.from v' :: VRep v)
  -- | Additive inverse
  negateV :: v -> v
  default negateV :: (Generic v, AdditiveGroup (VRep v)) => v -> v
  negateV v = Gnrx.to (negateV $ Gnrx.from v :: VRep v)
  -- | Group subtraction
  (^-^) :: v -> v -> v
  v ^-^ v' = v ^+^ negateV v'

-- | Sum over several vectors
sumV :: (Foldable f, AdditiveGroup v) => f v -> v
sumV = foldr (^+^) zeroV

instance AdditiveGroup () where
  zeroV     = ()
  () ^+^ () = ()
  negateV   = id

-- For 'Num' types:
--
-- instance AdditiveGroup n where {zeroV=0; (^+^) = (+); negateV = negate}

#define ScalarTypeCon(con,t) \
  instance con => AdditiveGroup (t) where {zeroV=0; (^+^) = (+); negateV = negate}

#define ScalarType(t) ScalarTypeCon((),t)

ScalarType(Int)
ScalarType(Integer)
ScalarType(Float)
ScalarType(Double)
ScalarType(CSChar)
ScalarType(CInt)
ScalarType(CShort)
ScalarType(CLong)
ScalarType(CLLong)
ScalarType(CIntMax)
ScalarType(CFloat)
ScalarType(CDouble)
ScalarTypeCon(Integral a,Ratio a)

instance (RealFloat v, AdditiveGroup v) => AdditiveGroup (Complex v) where
  zeroV   = zeroV :+ zeroV
  (^+^)   = (+)
  negateV = negate

-- Hm.  The 'RealFloat' constraint is unfortunate here.  It's due to a
-- questionable decision to place 'RealFloat' into the definition of the
-- 'Complex' /type/, rather than in functions and instances as needed.

instance (AdditiveGroup u,AdditiveGroup v) => AdditiveGroup (u,v) where
  zeroV             = (zeroV,zeroV)
  (u,v) ^+^ (u',v') = (u^+^u',v^+^v')
  negateV (u,v)     = (negateV u,negateV v)

instance (AdditiveGroup u,AdditiveGroup v,AdditiveGroup w)
    => AdditiveGroup (u,v,w) where
  zeroV                  = (zeroV,zeroV,zeroV)
  (u,v,w) ^+^ (u',v',w') = (u^+^u',v^+^v',w^+^w')
  negateV (u,v,w)        = (negateV u,negateV v,negateV w)

instance (AdditiveGroup u,AdditiveGroup v,AdditiveGroup w,AdditiveGroup x)
    => AdditiveGroup (u,v,w,x) where
  zeroV                       = (zeroV,zeroV,zeroV,zeroV)
  (u,v,w,x) ^+^ (u',v',w',x') = (u^+^u',v^+^v',w^+^w',x^+^x')
  negateV (u,v,w,x)           = (negateV u,negateV v,negateV w,negateV x)

-- Standard instance for an applicative functor applied to a vector space.
instance AdditiveGroup v => AdditiveGroup (a -> v) where
  zeroV   = pure   zeroV
  (^+^)   = liftA2 (^+^)
  negateV = fmap   negateV

-- Maybe is handled like the Maybe-of-Sum monoid
instance AdditiveGroup a => AdditiveGroup (Maybe a) where
  zeroV = Nothing
  Nothing ^+^ b'      = b'
  a' ^+^ Nothing      = a'
  Just a' ^+^ Just b' = Just (a' ^+^ b')
  negateV = fmap negateV

{-

Alexey Khudyakov wrote:

  I looked through vector-space package and found lawless instance. Namely Maybe's AdditiveGroup instance

  It's group so following relation is expected to hold. Otherwise it's not a group.
  > x ^+^ negateV x == zeroV

  Here is counterexample:

  > let x = Just 2 in x ^+^ negateV x == zeroV
  False

  I think it's not possible to sensibly define group instance for
  Maybe a at all.

I see that the problem here is in distinguishing 'Just zeroV' from
Nothing. I could fix the Just + Just line to use Nothing instead of Just
zeroV when a' ^+^ b' == zeroV, although doing so would require Eq a and
hence lose some generality. Even so, the abstraction leak would probably
show up elsewhere.

Hm.

-}

-- Memo tries
instance (HasTrie u, AdditiveGroup v) => AdditiveGroup (u :->: v) where
  zeroV   = pure   zeroV
  (^+^)   = liftA2 (^+^)
  negateV = fmap   negateV

-- | Monoid under group addition.  Alternative to the @Sum@ in
-- "Data.Monoid", which uses 'Num' instead of 'AdditiveGroup'.
newtype Sum a = Sum { getSum :: a }
  deriving (Eq, Ord, Read, Show, Bounded)

instance Functor Sum where
  fmap f (Sum a) = Sum (f a)

-- instance Applicative Sum where
--   pure a = Sum a
--   Sum f <*> Sum x = Sum (f x)

instance Applicative Sum where
  pure  = Sum
  (<*>) = inSum2 ($)

instance AdditiveGroup a => Semigroup (Sum a) where
  (<>) = liftA2 (^+^)

instance AdditiveGroup a => Monoid (Sum a) where
  mempty  = Sum zeroV
#if !(MIN_VERSION_base(4,11,0))
  mappend = (<>)
#endif

-- | Application a unary function inside a 'Sum'
inSum :: (a -> b) -> (Sum a -> Sum b)
inSum = getSum ~> Sum

-- | Application a binary function inside a 'Sum'
inSum2 :: (a -> b -> c) -> (Sum a -> Sum b -> Sum c)
inSum2 = getSum ~> inSum

instance AdditiveGroup a => AdditiveGroup (Sum a) where
  zeroV   = mempty
  (^+^)   = mappend
  negateV = inSum negateV

---- to go elsewhere

(~>) :: (a' -> a) -> (b -> b') -> ((a -> b) -> (a' -> b'))
(i ~> o) f = o . f . i

-- result :: (b -> b') -> ((a -> b) -> (a -> b'))
-- result = (.)

-- argument :: (a' -> a) -> ((a -> b) -> (a' -> b))
-- argument = flip (.)

-- g ~> f = result g . argument f

instance AdditiveGroup a => AdditiveGroup (Gnrx.Rec0 a s) where
  zeroV = Gnrx.K1 zeroV
  negateV (Gnrx.K1 v) = Gnrx.K1 $ negateV v
  Gnrx.K1 v ^+^ Gnrx.K1 w = Gnrx.K1 $ v ^+^ w
  Gnrx.K1 v ^-^ Gnrx.K1 w = Gnrx.K1 $ v ^-^ w
instance AdditiveGroup (f p) => AdditiveGroup (Gnrx.M1 i c f p) where
  zeroV = Gnrx.M1 zeroV
  negateV (Gnrx.M1 v) = Gnrx.M1 $ negateV v
  Gnrx.M1 v ^+^ Gnrx.M1 w = Gnrx.M1 $ v ^+^ w
  Gnrx.M1 v ^-^ Gnrx.M1 w = Gnrx.M1 $ v ^-^ w
instance (AdditiveGroup (f p), AdditiveGroup (g p)) => AdditiveGroup ((f :*: g) p) where
  zeroV = zeroV :*: zeroV
  negateV (x:*:y) = negateV x :*: negateV y
  (x:*:y) ^+^ (ξ:*:υ) = (x^+^ξ) :*: (y^+^υ)
  (x:*:y) ^-^ (ξ:*:υ) = (x^-^ξ) :*: (y^-^υ)


text:
module KalmanAxioms
where

import Numeric.LinearAlgebra
import Unsafe.Coerce

type Float_t = Double
type SDM = Matrix Float_t
type SSM = Matrix Float_t

float_zero :: Float_t
float_zero = 0

float_one :: Float_t
float_one = 1

float_opp :: Float_t -> Float_t
float_opp = negate

float_plus :: Float_t -> Float_t -> Float_t
float_plus = (+)

float_minus :: Float_t -> Float_t -> Float_t
float_minus = (-)

float_times :: Float_t -> Float_t -> Float_t
float_times = (*)

float_div :: Float_t -> Float_t -> Float_t
float_div = (/)

float_inv :: Float_t -> Float_t
float_inv = (1 /)

float_eq_dec :: Float_t -> Float_t -> Bool
float_eq_dec = (==)

float_lt :: a -> a' -> Bool
float_lt x y = (unsafeCoerce x :: Float_t) < (unsafeCoerce y :: Float_t)

cholesky_DC :: a -> Int -> sdm -> sdm'
cholesky_DC _ _ a = unsafeCoerce (chol (trustSym (unsafeCoerce a :: SDM)))

dense_sparse_mul_to_sparse :: a -> Int -> sdm -> ssm -> ssm'
dense_sparse_mul_to_sparse _ _ a b = unsafeCoerce ((unsafeCoerce a :: SDM) <> (unsafeCoerce b :: SSM))

sparse_dense_mul :: a -> Int -> ssm -> sdm -> sdm'
sparse_dense_mul _ _ a b = unsafeCoerce ((unsafeCoerce a :: SSM) <> (unsafeCoerce b :: SDM))

dense_sparse_mul :: a -> Int -> sdm -> ssm -> sdm'
dense_sparse_mul _ _ a b = unsafeCoerce ((unsafeCoerce a :: SDM) <> (unsafeCoerce b :: SSM))

solveR_upper :: a -> Int -> sdm -> sdm -> sdm'
solveR_upper _ _ a b = unsafeCoerce ((unsafeCoerce a :: SDM) <\> (unsafeCoerce b :: SDM))

solveR_lower :: a -> Int -> sdm -> sdm -> sdm'
solveR_lower _ _ a b = unsafeCoerce ((unsafeCoerce a :: SDM) <\> (unsafeCoerce b :: SDM))

dense_get :: a -> Int -> Int -> sdm -> Int -> Int -> Float_t
dense_get _ _ _ m i j = (unsafeCoerce m :: SDM) `atIndex` (i, j)

dense_mul :: a -> Int -> Int -> Int -> sdm -> sdm' -> sdm''
dense_mul _ _ _ _ a b = unsafeCoerce ((unsafeCoerce a :: SDM) <> (unsafeCoerce b :: SDM))

dense_fill :: a -> Int -> Int -> (Int -> Int -> Float_t) -> sdm
dense_fill _ m n f = unsafeCoerce (fromLists [ [f i j | j <- [0..n - 1]] | i <- [0..m-1]] :: SDM)

dense_elementwise_op :: a -> Int -> Int -> (Float_t -> Float_t -> Float_t) -> SDM -> SDM -> SDM
dense_elementwise_op me m n f a b = dense_fill me m n (\i j -> f (a `atIndex` (i, j)) (b `atIndex` (i, j)))

sparse_get :: a -> Int -> Int -> SSM -> Int -> Int -> Float_t
sparse_get _ _ _ m i j = m `atIndex` (i, j)

sparse_mul :: a -> Int -> Int -> Int -> SSM -> SSM -> SSM
sparse_mul _ _ _ _ a b = a <> b

sparse_fill :: a ->  Int -> Int -> (Int -> Int -> Float_t) -> SSM
sparse_fill _ m n f = fromLists [ [f i j | j <- [0..n - 1]] | i <- [0..m-1]]

sparse_elementwise_op :: a -> Int -> Int -> (Float_t -> Float_t -> Float_t) -> SSM -> SSM -> SSM
sparse_elementwise_op me m n f a b = sparse_fill me m n (\i j -> f (a `atIndex` (i, j)) (b `atIndex` (i, j)))


text:
module Numeric.SpecFunctions ( module X ) where

import Numeric.SpecFunctions.Internal as X

max_stars_repo_name: runeksvendsen/hs-gauge

max_issues_repo_name: runeksvendsen/hs-gauge

max_forks_repo_name: runeksvendsen/hs-gauge


text:
{-# LANGUAGE FlexibleInstances #-}
import Criterion.Main
import PCA.QModels
import PCA.Distribution
import Numeric.LinearAlgebra.Data (Vector, Matrix, fromRows, toRows)
import qualified Numeric.LinearAlgebra.HMatrix as H

import Control.DeepSeq

instance NFData (Distr Double Double) where
  rnf a = seq a ()

instance NFData (Distr (Vector Double) (Vector Double)) where
  rnf a = seq a ()

instance NFData (Distr (Vector Double) (Matrix Double)) where
  rnf a = seq a ()

instance NFData (Distr (Matrix Double) (Matrix Double)) where
  rnf a = seq a ()

prepareData n d = do
  trainT <- H.rand n d
  let tau = initQtau
      mu = initQmu d
      alpha = initQalpha d
      w = initQW d
      sigmaX = calcSigmaX tau w
  return (trainT, tau, mu, alpha, w, sigmaX)

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

main =
    defaultMain
        [ env (prepareData 100 3) $
          \small ->
               bgroup
                   "Small 100x3"
                   [ bench "Fool" $
                     whnf
                         (\(trainT,tau,mu,alpha,w,sigmaX) ->
                               foolCalcMX trainT tau sigmaX w mu)
                         small
                   , bench "Matrix" $
                     whnf
                         (\(trainT,tau,mu,alpha,w,sigmaX) ->
                               calcMX trainT tau sigmaX w mu)
                         small]
        , env (prepareData 10000 30) $
          \big ->
               bgroup
                   "Big 10000x30"
                   [ bench "Fool" $
                     whnf
                         (\(trainT,tau,mu,alpha,w,sigmaX) ->
                               foolCalcMX trainT tau sigmaX w mu)
                         big
                   , bench "Matrix" $
                     whnf
                         (\(trainT,tau,mu,alpha,w,sigmaX) ->
                               calcMX trainT tau sigmaX w mu)
                         big]]

max_stars_repo_name: DbIHbKA/vbpca

max_issues_repo_name: DbIHbKA/vbpca

max_forks_repo_name: DbIHbKA/vbpca


text:
--
-- Image : image data structure
--

module CNN.Image (
  Image
, Plain
, Class
, Trainer
, classNumToVec
) where

import qualified Data.Map as M
import Numeric.LinearAlgebra.Data

type Plain = Matrix R    -- 2D: X x Y pixels
type Image = [Plain]     -- n dimension

type Class = Vector R    -- teacher vector

type Trainer = (Image, Class)

-- FUNCTIONS

classNumToVec :: Int -> Int -> Class
classNumToVec n c = fromList ls
  where
    ls = take n (replicate c 0.0 ++ [1.0] ++ repeat 0.0)

max_stars_repo_name: eijian/deeplearning

max_issues_repo_name: eijian/deeplearning

max_forks_repo_name: eijian/deeplearning


text:
{-# LANGUAGE GeneralizedNewtypeDeriving #-}

module Valuation
    ( module Valuation
    , throwError
    ) where

import Control.Applicative
import Control.Monad.Except
import Control.Monad.Identity
import Control.Monad.Reader
import Statistics.Distribution
import Statistics.Distribution.Normal

type Year = Int

newtype Scenario a = Scenario
    { scenario :: a
    }

newtype Value s a = Value
    { runValue :: ReaderT (Scenario s) (ExceptT String Identity) a
    } deriving (Monad, Applicative, Alternative, Functor, MonadReader (Scenario s), MonadError String)

instance Num a => Num (Value s a) where
    negate = fmap negate
    (+) = liftM2 (+)
    (*) = liftM2 (*)
    fromInteger = pure . fromInteger
    abs = fmap abs
    signum = fmap signum

instance Fractional a => Fractional (Value s a) where
    (/) = liftM2 (/)
    fromRational = pure . fromRational
    recip x = 1 / x

type TimeValue s a = Year -> Value s a

runValuation :: s -> Value s a -> Either String a
runValuation s v = runIdentity (runExceptT (runReaderT (runValue v) (Scenario s)))

runV :: (Show a, RealFrac a) => TimeValue () a -> [Year] -> IO ()
runV v (y:ys) = do
    putStr $ show y
    putStr "  "
    case runValuation () $ v y of
        Left s -> print s
        Right result -> print result
    runV v ys
runV _ [] = return ()

discount :: TimeValue s Double -> Year -> TimeValue s Double -> TimeValue s Double
discount r t0 v t
    | t == t0 = v t
    | otherwise = discount r (t0 + 1) v t / (1 + r (t0 + 1))

npv :: TimeValue s Double -> TimeValue s Double -> TimeValue s Double
npv r v t = v t + (npv r v (t + 1) / (1 + r (t + 1)) <|> return 0)

linear :: Fractional a => TimeValue s a -> Year -> Year -> a -> TimeValue s a
linear v t0 t1 target t
    | t0 < t && t <= t1 = do
        base <- v t0
        return $ base + (target - base) * fromIntegral (t - t0) / fromIntegral (t1 - t0)
    | t1 < t = return target

pvAnnuity :: Double -> Double -> Int -> Double
pvAnnuity r pmt n = pmt * (1 - (1 + r) ^^ (-n)) / r

blackScholes :: Double -> Double -> Double -> Double -> Double -> Double
blackScholes sigma r t s k =
    let d1 = 1 / (sigma * sqrt t) * (log (s / k) + (r + sigma ** t / 2) * t)
        d2 = d1 - sigma * sqrt t
    in s * cumulative standard d1 - cumulative standard d2 * k * exp (-r * t)

max_stars_repo_name: sboehler/haskell-valuation

max_issues_repo_name: sboehler/haskell-valuation

max_forks_repo_name: sboehler/haskell-valuation


text:
{-# LANGUAGE ForeignFunctionInterface #-}
-----------------------------------------------------------------------------
{- |
Module      : Numeric.GSL.Fourier
Copyright   :  (c) Alberto Ruiz 2006
License     :  GPL-style

Maintainer  :  Alberto Ruiz (aruiz at um dot es)
Stability   :  provisional
Portability :  uses ffi

Fourier Transform.

<http://www.gnu.org/software/gsl/manual/html_node/Fast-Fourier-Transforms.html#Fast-Fourier-Transforms>

-}
-----------------------------------------------------------------------------
module Numeric.GSL.Fourier (
    fft,
    ifft
) where

import Data.Packed.Internal
import Data.Complex
import Foreign.C.Types
import System.IO.Unsafe (unsafePerformIO)

genfft code v = unsafePerformIO $ do
    r <- createVector (dim v)
    app2 (c_fft code) vec v vec r "fft"
    return r

foreign import ccall unsafe "gsl-aux.h fft" c_fft ::  CInt -> TCVCV

{- | Fast 1D Fourier transform of a 'Vector' @(@'Complex' 'Double'@)@ using /gsl_fft_complex_forward/. It uses the same scaling conventions as GNU Octave.

@> fft ('fromList' [1,2,3,4])
vector (4) [10.0 :+ 0.0,(-2.0) :+ 2.0,(-2.0) :+ 0.0,(-2.0) :+ (-2.0)]@

-}
fft :: Vector (Complex Double) -> Vector (Complex Double)
fft = genfft 0

-- | The inverse of 'fft', using /gsl_fft_complex_inverse/.
ifft :: Vector (Complex Double) -> Vector (Complex Double)
ifft = genfft 1

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
-----------------------------------------------------------------------------
-- |
-- Module      :  Graphics.Rendering.Plot.Render.Plot
-- Copyright   :  (c) A. V. H. McPhail 2010
-- License     :  BSD3
--
-- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
-- Stability   :  provisional
-- Portability :  portable
--
-- Rendering 'Figure's
--
-----------------------------------------------------------------------------

module Graphics.Rendering.Plot.Render.Plot (
                                       -- * Rendering
                                       renderPlots
                                       ) where

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

--import Data.Either

--import Data.Packed.Vector
--import Numeric.LinearAlgebra.Linear

--import Data.Word

--import Data.Maybe

--import Data.Colour.SRGB
--import Data.Colour.Names

import qualified Data.Array.IArray as A

import Data.Colour.Names

import qualified Graphics.Rendering.Cairo as C
--import qualified Graphics.Rendering.Pango as P

--import Control.Monad.Reader
import Control.Monad.State
--import Control.Monad.Trans

import Graphics.Rendering.Plot.Types
--import Graphics.Rendering.Plot.Defaults

--import Graphics.Rendering.Plot.Figure.Text

import Graphics.Rendering.Plot.Render.Types

import Graphics.Rendering.Plot.Render.Text
import Graphics.Rendering.Plot.Render.Plot.Axis
import Graphics.Rendering.Plot.Render.Plot.Data
import Graphics.Rendering.Plot.Render.Plot.Legend
import Graphics.Rendering.Plot.Render.Plot.Annotation

--import qualified Text.Printf as Printf

--import Prelude hiding(min,max)
--import qualified Prelude(max)

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

bbPlot :: Int -> Int -> (Int,Int) -> Render ()
bbPlot r c (px,py) = modify (\(BoundingBox x y w h) ->
  let w' = w/(fromIntegral c)
      h' = h/(fromIntegral r)
  in (BoundingBox
      (x+w'*((fromIntegral py)-1))
      (y+h'*((fromIntegral px)-1))
      w' h'))

renderPlots :: Plots -> Render ()
renderPlots d = do
  let ((x,y),(x',y')) = A.bounds d
      rows = x'-x+1
      cols = y'-y+1
  bb <- get
  mapM_ (\(i,e) -> do
    case e of
      Nothing -> return ()
      Just e' -> do
        bbPlot rows cols i
        renderPlot e'
        put bb) (A.assocs d)

renderPlot :: PlotData -> Render ()
renderPlot (Plot b c p hd r a bc sd d l an) = do
  tx <- bbCentreWidth
  ty <- bbTopHeight
  (_,th) <- renderText hd Centre TTop tx ty
  when (th /= 0) $ bbLowerTop (th+textPad)
  legend <- renderLegend l d
  (axes,padding) <- renderAxes p r a
  renderBorder b
  cairo C.save
  clipBoundary
  when (c /= white) (do
    cairo $ do
      setColour c
      C.paint)
  renderData r bc sd d
  renderAnnotations r an
  cairo C.restore
  cairo C.save
  legend padding
  cairo C.restore
  cairo C.save
  axes
  cairo C.restore

renderBorder :: Border -> Render ()
renderBorder False = return ()
renderBorder True  = do
  (BoundingBox x y w h) <- get
  cairo $ do
    C.setLineWidth 0.5
    C.rectangle (x+0.5) (y+0.5) w h
    C.stroke

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

max_stars_repo_name: JackTheEngineer/plot


text:
module HScheme.Parser where

import HScheme.Data

import Text.ParserCombinators.Parsec as P hiding ( spaces )
import Data.Char ( toUpper )
import Data.Ratio ( (%) )
import Data.Complex ( Complex((:+)) )
import Numeric ( readHex, readOct )
import Data.List ( foldl' )
import Control.Monad.Except
import qualified Data.Vector as V

symbol :: Parser Char
symbol = oneOf "!$%&|*+/:<=>?@^_~"

spaces :: Parser ()
spaces = skipMany1 space

parseBool :: Parser Value
parseBool = try $ do
    char '#'
    (char 't' >> return (Bool True)) <|> (char 'f' >> return (Bool False))

parseString :: Parser Value
parseString = do
    char '"'
    x <- many (noneOf "\"" <|> escapeChar)
    char '"'
    return $ String x
    where
        escapeChar :: Parser Char
        escapeChar = do
            oneOf "nrt\\" >>= \case
                'n' -> return '\n'
                'r' -> return '\r'
                't' -> return '\t'
                '\\'-> return '\\'

parseChar :: Parser Value
parseChar = do
    try $ string "#\\"
    chr <- (string "newline" <|> string "space") <|>
        do
            x <- anyChar
            notFollowedBy alphaNum
            return [x]
    return $ Character $ case chr of
        "\\" -> '\n'
        "newline" -> '\n'
        "space" -> ' '
        _ -> head chr

parseAtom :: Parser Value
parseAtom = do
    first <- letter <|> symbol
    rest <- many (letter <|> digit <|> symbol)
    let atom = first:rest
    return $ case atom of
        "#t" -> Bool True
        "#f" -> Bool False
        _    -> Atom atom

{-
    Parsing numbers
-}

parseDecimal :: Parser Value
parseDecimal = do
    sign <- option "" $ string "-"
    many1 digit >>= (return . Number . read . (sign <>))

parseDecimal' :: Parser Value
parseDecimal' = try $ string "#d" >> many1 digit >>= (return . Number . read)

-- | Parse hexadecimal number of form #x...
parseHex :: Parser Value
parseHex = try $ string "#x" >> many1 hexDigit >>= (return . Number . hexToDigit)
    where
        hexToDigit x = fst $ head $ readHex x

-- | Parse octadecimal number of form #o...
parseOct :: Parser Value
parseOct = try $ string "#o" >> many1 octDigit >>= (return . Number . octToDigit)
    where
        octToDigit x = fst $ head $ readOct x

-- | Parse binary of form #b...
parseBin :: Parser Value
parseBin = try $ string "#b" >> many1 (oneOf "10") >>= (return . Number . binToDigit)
    where
        binToDigit = foldl' binDigitToDigit 0
            where
                binDigitToDigit acc = \case
                    '1' -> 2 * acc + 1
                    '0' -> 2 * acc + 0

parseNumber :: Parser Value
parseNumber = do
    parseDecimal
    <|> parseDecimal'
    <|> parseHex
    <|> parseOct
    <|> parseBin

parseFloat :: Parser Value
parseFloat = try $ do
    sign <- option "" $ string "-"
    x <- many1 digit
    y <- (:) <$> char '.' <*> many1 digit
    return $ Float $ read $ sign <> x <> y

parseRatio :: Parser Value
parseRatio = try $ do
    (Number x) <- parseDecimal
    char '/'
    (Number y) <- parseDecimal
    return $ Ratio (x % y)

parseComplex :: Parser Value
parseComplex = try $ do
    real <- parseFloat <|> parseNumber
    char '+'
    imiginary <- parseFloat <|> parseNumber
    char 'i'
    return $ Complex (toDouble real :+ toDouble imiginary)
    where
        toDouble (Float f) = realToFrac f
        toDouble (Number n) = realToFrac n

parseExpr :: Parser Value
parseExpr = parseAtom
    <|> parseString
    <|> parseChar
    <|> parseRatio
    <|> parseComplex
    <|> parseFloat
    <|> parseNumber
    <|> parseBool
    <|> parseQuoted
    <|> parseVector
    <|> parseList

{-
    Recursive Parsers
-}

parseList :: Parser Value
parseList = do
    char '('
    x <- try parseSpacesList <|> parseDottedList
    char ')'
    return x

parseSpacesList :: Parser Value
parseSpacesList = List <$> sepBy parseExpr spaces

parseDottedList :: Parser Value
parseDottedList = do
    head <- endBy parseExpr spaces
    tail <- char '.' >> spaces >> parseExpr
    return $ DottedList head tail

parseQuoted :: Parser Value
parseQuoted = do
    char '\''
    x <- parseExpr
    return $ List [Atom "quote", x]

parseVector :: Parser Value
parseVector = do
    try $ string "#("
    values <- sepBy parseExpr spaces
    char ')'
    return $ Vector (V.fromList values)

readExpr :: String -> ThrowsErr Value
readExpr input = case parse parseExpr "" input of
    Left err -> throwError $ Parser err
    Right val -> return val

max_stars_repo_name: dvdvgt/hScheme

max_issues_repo_name: dvdvgt/hScheme

max_forks_repo_name: dvdvgt/hScheme


text:
module Main where

import System.Environment (getArgs)
import qualified Statistics.Test.Types as Stats
import Jakway.Blackjack.Random.Checks
import Jakway.Blackjack.Random.Options
import Jakway.Blackjack.Util (die)     -- ^ die is implemented here for compatibility with older GHC versions
import Jakway.Blackjack.AI

main :: IO ()
main = do
    args <- getArgs
    conf <- getConfig args

    -- | TODO: make the RNG max range an option
    let rngMaxRange = toInteger $ 2 ^ 25

    let distrib = distribution conf
        pval = pvalue conf
        n = sampleSize conf
        r = rngSampleSize conf
        dealerAI = BasicDealer
        playerAIs = [BasicPlayer]

    if (n <= 104) then die "Samples must be >104." else return ()

    case distrib of EvenDistribution -> testDeckEvenDistribution pval n dealerAI playerAIs >>= printResult
                    RNGDistribution -> testRNGDistribution pval n r rngMaxRange dealerAI playerAIs >>= printResult
                    _ -> die "Fatal error, unknown test!"

printResult :: Stats.TestResult -> IO ()
printResult (Stats.Significant) = putStrLn $ "The test result is: Significant.  There is evidence to conclude that the deck IS NOT an unbiased source of randomness."
printResult (Stats.NotSignificant) = putStrLn $ "The test result is: Not Significant.  There is insufficient evidence to conclude that the deck is not an unbiased source of randomness."

max_stars_repo_name: tjakway/haskell-blackjack

max_issues_repo_name: tjakway/haskell-blackjack

max_forks_repo_name: tjakway/haskell-blackjack