text:
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE StandaloneDeriving    #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE DeriveAnyClass #-}
{-|
Module      : Grenade.Core.Pooling
Description : Max Pooling layer for 2D and 3D images
Copyright   : (c) Huw Campbell, 2016-2017
License     : BSD2
Stability   : experimental
-}
module Grenade.Layers.Pooling (
    Pooling (..)
  ) where

import           Data.Maybe
import           Data.Proxy
import           Data.Serialize
import           Data.Singletons.TypeLits hiding (natVal)
import           GHC.TypeLits
import           Control.DeepSeq (NFData)
import           GHC.Generics    (Generic)

import           Grenade.Core
import           Grenade.Layers.Internal.Pooling

import           Numeric.LinearAlgebra.Static as LAS hiding ((|||), build, toRows)

-- | A pooling layer for a neural network.
--
--   Does a max pooling, looking over a kernel similarly to the convolution network, but returning
--   maxarg only. This layer is often used to provide minor amounts of translational invariance.
--
--   The kernel size dictates which input and output sizes will "fit". Fitting the equation:
--   `out = (in - kernel) / stride + 1` for both dimensions.
--
data Pooling :: Nat -> Nat -> Nat -> Nat -> * where
  Pooling :: Pooling kernelRows kernelColumns strideRows strideColumns
  deriving (NFData, Generic)

instance Show (Pooling k k' s s') where
  show Pooling = "Pooling"

instance UpdateLayer (Pooling kernelRows kernelColumns strideRows strideColumns) where
  type Gradient (Pooling kernelRows kernelColumns strideRows strideColumns) = ()
  runUpdate _ Pooling _ = Pooling

instance RandomLayer (Pooling k k' s s') where
  createRandomWith _ _ = return Pooling

instance Serialize (Pooling kernelRows kernelColumns strideRows strideColumns) where
  put _ = return ()
  get = return Pooling

-- | A two dimentional image can be pooled.
instance ( KnownNat kernelRows
         , KnownNat kernelColumns
         , KnownNat strideRows
         , KnownNat strideColumns
         , KnownNat inputRows
         , KnownNat inputColumns
         , KnownNat outputRows
         , KnownNat outputColumns
         , ((outputRows - 1) * strideRows) ~ (inputRows - kernelRows)
         , ((outputColumns - 1) * strideColumns) ~ (inputColumns - kernelColumns)
         ) => Layer (Pooling kernelRows kernelColumns strideRows strideColumns) ('D2 inputRows inputColumns) ('D2 outputRows outputColumns) where
  type Tape (Pooling kernelRows kernelColumns strideRows strideColumns) ('D2 inputRows inputColumns) ('D2 outputRows outputColumns) = S ('D2 inputRows inputColumns)
  runForwards Pooling (S2D input) =
    let height = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        width  = fromIntegral $ natVal (Proxy :: Proxy inputColumns)
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelColumns)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideColumns)
        ex = extract input
        r  = poolForward 1 height width kx ky sx sy ex
        rs = fromJust . create $ r
    in  (S2D input, S2D rs)
  runBackwards Pooling (S2D input) (S2D dEdy) =
    let height = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        width  = fromIntegral $ natVal (Proxy :: Proxy inputColumns)
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelColumns)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideColumns)
        ex = extract input
        eo = extract dEdy
        vs = poolBackward 1 height width kx ky sx sy ex eo
    in  ((), S2D . fromJust . create $ vs)

-- | A three dimensional image can be pooled on each layer.
instance ( KnownNat kernelRows
         , KnownNat kernelColumns
         , KnownNat strideRows
         , KnownNat strideColumns
         , KnownNat inputRows
         , KnownNat inputColumns
         , KnownNat outputRows
         , KnownNat outputColumns
         , KnownNat channels
         , KnownNat (outputRows * channels)
         , ((outputRows - 1) * strideRows) ~ (inputRows - kernelRows)
         , ((outputColumns - 1) * strideColumns) ~ (inputColumns - kernelColumns)
         ) => Layer (Pooling kernelRows kernelColumns strideRows strideColumns) ('D3 inputRows inputColumns channels) ('D3 outputRows outputColumns channels) where
  type Tape (Pooling kernelRows kernelColumns strideRows strideColumns) ('D3 inputRows inputColumns channels) ('D3 outputRows outputColumns channels) = S ('D3 inputRows inputColumns channels)
  runForwards Pooling (S3D input) =
    let ix = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        iy = fromIntegral $ natVal (Proxy :: Proxy inputColumns)
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelColumns)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideColumns)
        ch = fromIntegral $ natVal (Proxy :: Proxy channels)
        ex = extract input
        r  = poolForward ch ix iy kx ky sx sy ex
        rs = fromJust . create $ r
    in  (S3D input, S3D rs)
  runBackwards Pooling (S3D input) (S3D dEdy) =
    let ix = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        iy = fromIntegral $ natVal (Proxy :: Proxy inputColumns)
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelColumns)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideColumns)
        ch = fromIntegral $ natVal (Proxy :: Proxy channels)
        ex = extract input
        eo = extract dEdy
        vs = poolBackward ch ix iy kx ky sx sy ex eo
    in  ((), S3D . fromJust . create $ vs)

-------------------- GNum instance --------------------

instance GNum (Pooling k k' s s') where
  _ |* _ = Pooling
  _ |+ _ = Pooling
  gFromRational _ = Pooling

max_stars_repo_name: koenigmaximilian/grenade

max_issues_repo_name: koenigmaximilian/grenade

max_forks_repo_name: koenigmaximilian/grenade


text:
module LabelMatrix where

-- import Data.Function    ((&))
import Foreign.Storable (Storable)

import           Data.Map              (Map)
-- import qualified Data.Map              as M
import qualified Data.Vector  as V
-- import qualified Data.Vector.Storable  as VS
import           Numeric.LinearAlgebra hiding (fromList)

{-
Matrices with labeled rows and columns.

Construction functions handle uniqueness & correct length of labels.

Provides similar convenience functions to pandas DataFrames. Contrary to the latter, all
elements in the matrix contain only one type.
-}

--bla = ((12><12) [0..] :: Matrix I) #> fromList [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0] :: Vector I

data LabelMatrix a r c =
    LabelMatrix
        { matrix :: Matrix a
        , rows   :: V.Vector r
        , cols   :: V.Vector c
        } deriving (Show)

fromMatrix :: Matrix a -> [r] -> [c] -> LabelMatrix a r c
fromMatrix = undefined

fromList :: Storable a => [a] -> [r] -> [c] -> LabelMatrix a r c
fromList content rowLabels columnLabels =
    let
        m = length rowLabels
        n = length columnLabels
        mat = (m><n) content
    in
        LabelMatrix mat (V.fromList rowLabels) (V.fromList columnLabels)

{-
v = m ¿ [0] -- extract column(s)

flatten v -- makes it into a vector, especially nice if the matrix is nx1 anyway.
-}

groupBy :: LabelMatrix a r c -> r -> Map a (LabelMatrix a r c)
groupBy = groupByRow

groupByRow :: LabelMatrix a r c -> r -> Map a (LabelMatrix a r c)
groupByRow = undefined

groupByColumn :: LabelMatrix a r c -> c -> Map a (LabelMatrix a r c)
groupByColumn = undefined

ext: hs


text:
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE UndecidableInstances #-}
module Math.HiddenMarkovModel.Private where

import qualified Math.HiddenMarkovModel.Distribution as Distr
import qualified Math.HiddenMarkovModel.CSV as HMMCSV
import Math.HiddenMarkovModel.Distribution (State(State))

import qualified Numeric.LinearAlgebra.Algorithms as Algo
import qualified Numeric.LinearAlgebra.Util as LinAlg
import qualified Numeric.Container as NC
import qualified Data.Packed.Development as Dev
import qualified Data.Packed.Matrix as Matrix
import qualified Data.Packed.Vector as Vector
import Numeric.Container ((<>))
import Data.Packed.Matrix (Matrix)
import Data.Packed.Vector (Vector)

import Control.DeepSeq (NFData, rnf)
import Foreign.Storable (Storable)

import qualified Data.NonEmpty.Class as NonEmptyC
import qualified Data.NonEmpty as NonEmpty
import qualified Data.Semigroup as Sg
import qualified Data.List as List
import Data.Traversable (Traversable, mapAccumL)
import Data.Tuple.HT (mapPair, mapFst, mapSnd, swap)

{- |
A Hidden Markov model consists of a number of (hidden) states
and a set of emissions.
There is a vector for the initial probability of each state
and a matrix containing the probability for switching
from one state to another one.
The 'distribution' field points to probability distributions
that associate every state with emissions of different probability.
Famous distribution instances are discrete and Gaussian distributions.
See "Math.HiddenMarkovModel.Distribution" for details.

The transition matrix is transposed
with respect to popular HMM descriptions.
But I think this is the natural orientation, because this way
you can write \"transition matrix times probability column vector\".

The type has two type parameters,
although the one for the distribution would be enough.
However, replacing @prob@ by @Distr.Probability distr@
would prohibit the derived Show and Read instances.
-}
data T distr prob =
   Cons {
      initial :: Vector prob,
      transition :: Matrix prob,
      distribution :: distr
   }
   deriving (Show, Read)

instance
   (NFData distr, NFData prob, Storable prob) =>
      NFData (T distr prob) where
   rnf hmm = rnf (initial hmm, transition hmm, distribution hmm)

emission ::
   (Distr.EmissionProb distr,
    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
   T distr prob -> emission -> Vector prob
emission  =  Distr.emissionProb . distribution

forward ::
   (Distr.EmissionProb distr,
    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,
    Traversable f) =>
   T distr prob -> NonEmpty.T f emission -> prob
forward hmm = NC.sumElements . NonEmpty.last . alpha hmm

alpha ::
   (Distr.EmissionProb distr,
    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,
    Traversable f) =>
   T distr prob ->
   NonEmpty.T f emission -> NonEmpty.T f (Vector prob)
alpha hmm (NonEmpty.Cons x xs) =
   NonEmpty.scanl
      (\alphai xi -> NC.mul (emission hmm xi) (transition hmm <> alphai))
      (NC.mul (emission hmm x) (initial hmm))
      xs

backward ::
   (Distr.EmissionProb distr,
    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,
    Traversable f) =>
   T distr prob -> NonEmpty.T f emission -> prob
backward hmm (NonEmpty.Cons x xs) =
   NC.sumElements $
   NC.mul (initial hmm) $
   NC.mul (emission hmm x) $
   NonEmpty.head $ beta hmm xs

beta ::
   (Distr.EmissionProb distr,
    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,
    Traversable f) =>
   T distr prob ->
   f emission -> NonEmpty.T f (Vector prob)
beta hmm =
   NonEmpty.scanr
      (\xi betai -> NC.mul (emission hmm xi) betai <> transition hmm)
      (NC.constant 1 (NC.dim $ initial hmm))

alphaBeta ::
   (Distr.EmissionProb distr,
    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,
    Traversable f) =>
   T distr prob ->
   NonEmpty.T f emission ->
   (prob, NonEmpty.T f (Vector prob), NonEmpty.T f (Vector prob))
alphaBeta hmm xs =
   let alphas = alpha hmm xs
       betas = beta hmm $ NonEmpty.tail xs
       recipLikelihood = recip $ NC.sumElements $ NonEmpty.last alphas
   in  (recipLikelihood, alphas, betas)

xiFromAlphaBeta ::
   (Distr.EmissionProb distr,
    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
   T distr prob -> prob ->
   NonEmpty.T [] emission ->
   NonEmpty.T [] (Vector prob) ->
   NonEmpty.T [] (Vector prob) ->
   [Matrix prob]
xiFromAlphaBeta hmm recipLikelihood xs alphas betas =
   zipWith3
      (\x alpha0 beta1 ->
         NC.scale recipLikelihood $
         NC.mul
            (NC.outer (NC.mul (emission hmm x) beta1) alpha0)
            (transition hmm))
      (NonEmpty.tail xs)
      (NonEmpty.init alphas)
      (NonEmpty.tail betas)

zetaFromXi ::
   (NC.Product prob) => T distr prob -> [Matrix prob] -> [Vector prob]
zetaFromXi hmm xis =
   map (NC.constant 1 (Matrix.rows $ transition hmm) <>) xis

zetaFromAlphaBeta ::
   (NC.Container Vector prob) =>
   prob ->
   NonEmpty.T [] (Vector prob) ->
   NonEmpty.T [] (Vector prob) ->
   NonEmpty.T [] (Vector prob)
zetaFromAlphaBeta recipLikelihood alphas betas =
   fmap (NC.scale recipLikelihood) $
   NonEmptyC.zipWith NC.mul alphas betas

{- |
In constrast to Math.HiddenMarkovModel.reveal
this does not normalize the vector.
This is slightly simpler but for long sequences
the product of probabilities might be smaller
than the smallest representable number.
-}
reveal ::
   (Distr.EmissionProb distr,
    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,
    Traversable f) =>
   T distr prob -> NonEmpty.T f emission -> NonEmpty.T f State
reveal hmm (NonEmpty.Cons x xs) =
   fmap State $
   uncurry (NonEmpty.scanr Dev.at') $
   mapFst NC.maxIndex $
   mapAccumL
      (\alphai xi ->
         swap $ mapSnd (NC.mul (emission hmm xi)) $
         matrixMaxMul (transition hmm) alphai)
      (NC.mul (emission hmm x) (initial hmm)) xs

matrixMaxMul ::
   (NC.Container Vector a) =>
   Matrix a -> Vector a -> (Vector Int, Vector a)
matrixMaxMul m v =
   mapPair (Vector.fromList, Vector.fromList) $ unzip $
   map ((\x -> (NC.maxIndex x, NC.maxElement x)) . NC.mul v) $
   Matrix.toRows m

{- |
A trained model is a temporary form of a Hidden Markov model
that we need during the training on multiple training sequences.
It allows to collect knowledge over many sequences with 'mergeTrained',
even with mixed supervised and unsupervised training.
You finish the training by converting the trained model
back to a plain modul using 'finishTraining'.

You can create a trained model in three ways:

* supervised training using an emission sequence with associated states,

* unsupervised training using an emission sequence and an existing Hidden Markov Model,

* derive it from state sequence patterns, cf. "Math.HiddenMarkovModel.Pattern".
-}
data Trained distr prob =
   Trained {
      trainedInitial :: Vector prob,
      trainedTransition :: Matrix prob,
      trainedDistribution :: distr
   }
   deriving (Show, Read)

instance
   (NFData distr, NFData prob, Storable prob) =>
      NFData (Trained distr prob) where
   rnf hmm =
      rnf (trainedInitial hmm, trainedTransition hmm, trainedDistribution hmm)

sumTransitions ::
   (NC.Container Vector e) =>
   T distr e -> [Matrix e] -> Matrix e
sumTransitions hmm =
   List.foldl' NC.add (NC.konst 0 $ LinAlg.size $ transition hmm)

{- |
Baum-Welch algorithm
-}
trainUnsupervised ::
   (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr,
    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
   T distr prob -> NonEmpty.T [] emission -> Trained tdistr prob
trainUnsupervised hmm xs =
   let (recipLikelihood, alphas, betas) = alphaBeta hmm xs
       zetas = zetaFromAlphaBeta recipLikelihood alphas betas

   in  Trained {
          trainedInitial = NonEmpty.head zetas,
          trainedTransition =
             sumTransitions hmm $
             xiFromAlphaBeta hmm recipLikelihood xs alphas betas,
          trainedDistribution =
             Distr.accumulateEmissions $ map (zip (NonEmpty.flatten xs)) $
             List.transpose $ map Vector.toList $ NonEmpty.flatten zetas
       }

mergeTrained ::
   (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr,
    Distr.Probability distr ~ prob) =>
   Trained tdistr prob -> Trained tdistr prob -> Trained tdistr prob
mergeTrained hmm0 hmm1 =
   Trained {
      trainedInitial = NC.add (trainedInitial hmm0) (trainedInitial hmm1),
      trainedTransition =
         NC.add (trainedTransition hmm0) (trainedTransition hmm1),
      trainedDistribution =
         Distr.combine
            (trainedDistribution hmm0) (trainedDistribution hmm1)
   }

instance
   (Distr.Estimate tdistr, Distr.Distribution tdistr ~ distr,
    Distr.Probability distr ~ prob) =>
      Sg.Semigroup (Trained tdistr prob) where
   (<>) = mergeTrained

toCells ::
   (Distr.CSV distr, Algo.Field prob, Show prob) =>
   T distr prob -> [[String]]
toCells hmm =
   (HMMCSV.cellsFromVector $ initial hmm) :
   (HMMCSV.cellsFromMatrix $ transition hmm) ++
   [] :
   (Distr.toCells $ distribution hmm)

parseCSV ::
   (Distr.CSV distr, Algo.Field prob, Read prob) =>
   HMMCSV.CSVParser (T distr prob)
parseCSV = do
   v <- HMMCSV.parseNonEmptyVectorCells
   m <- HMMCSV.parseSquareMatrixCells $ Vector.dim v
   HMMCSV.skipEmptyRow
   distr <- Distr.parseCells $ Vector.dim v
   return $ Cons {
      initial = v,
      transition = m,
      distribution = distr
   }

max_stars_repo_name: rybern/hmm-hmatrix

max_issues_repo_name: rybern/hmm-hmatrix

max_forks_repo_name: rybern/hmm-hmatrix


text:
module Polestar.PrettyPrint
  (prettyPrintTypeP
  ,prettyPrintType
  ,prettyPrintCanonicalTypeP
  ,prettyPrintICanonicalTypeP
  ,prettyPrintTermP
  ,prettyPrintTerm
  ) where
import Polestar.Type
import Polestar.Parse (NameBinding(..))
import Data.List (intersperse)
import Data.Complex

varNames :: [NameBinding] -> [String]
varNames = foldr (\x ys -> case x of
                             NVarBind n -> n:ys
                             _ -> ys) []

tyVarNames :: [NameBinding] -> [String]
tyVarNames = foldr (\x ys -> case x of
                               NTyVarBind n -> n:ys
                               _ -> ys) []

rename :: [String] -> String -> String
rename ctx name | name `notElem` ctx = name
                | otherwise = loop 1
  where loop :: Int -> String
        loop i = case name ++ show i of
          m | m `notElem` ctx -> m
          _ -> loop (i + 1)

showPrimType :: PrimType -> ShowS
showPrimType p = case p of
  PTyZero -> showString "Zero"
  PTyNat -> showString "Nat"
  PTyInt -> showString "Int"
  PTyNNReal -> showString "NNReal"
  PTyReal -> showString "Real"
  PTyImaginary -> showString "Imaginary"
  PTyComplex -> showString "Complex"
  PTyTrue -> showString "True"
  PTyFalse -> showString "False"
  PTyBool -> showString "Bool"
  PTyUnit -> showString "Unit"

-- <3>: <primitive type> | <type variable> | <0>
-- <2>: <3> -> <2>
-- <1>: <2> & .. & <2>
-- <0>: forall a <: <1>. <0> | <1>
prettyPrintTypeP :: Int -> [NameBinding] -> Type -> ShowS
prettyPrintTypeP p ctx t = case t of
  TyPrim p -> showPrimType p
  TyArr s t -> showParen (p > 2) $ prettyPrintTypeP 3 ctx s . showString " -> " . prettyPrintTypeP 2 (NAnonymousBind : ctx) t
  TyRef i | i < length ctx -> case ctx !! i of
              NTyVarBind n -> showString n
              n -> showString "<invalid type variable reference #" . shows i . showString ": " . shows n . showChar '>'
          | otherwise -> showString "<invalid type variable reference #" . shows i . showString ", index out of range>"
  TyAll (Id name) Nothing t -> showParen (p > 0) $ showString "forall " . showString name' . showString ". " . prettyPrintTypeP 0 (NTyVarBind name' : ctx) t
    where name' = rename (tyVarNames ctx) name
  TyAll (Id name) (Just b) t -> showParen (p > 0) $ showString "forall " . showString name' . showString "<:" . prettyPrintTypeP 1 ctx b . showString ". " . prettyPrintTypeP 0 (NTyVarBind name' : ctx) t
    where name' = rename (tyVarNames ctx) name
  TyInter tys -> showParen (p > 1) $ foldr (.) id $ intersperse (showString " & ") $ map (prettyPrintTypeP 2 ctx) tys
  TyTuple components -> showParen True $ foldr (.) id $ intersperse (showString ", ") $ map (prettyPrintTypeP 0 ctx) components

prettyPrintType :: Type -> String
prettyPrintType t = prettyPrintTypeP 0 [] t ""

prettyPrintCanonicalTypeP :: Int -> [NameBinding] -> CanonicalType -> ShowS
prettyPrintCanonicalTypeP p ctx [] = showString "<untyped>"
prettyPrintCanonicalTypeP p ctx [t] = prettyPrintICanonicalTypeP p ctx t
prettyPrintCanonicalTypeP p ctx ts = showParen (p > 1) $ foldl1 (\x y -> x . showString " & " . y) $ map (prettyPrintICanonicalTypeP 2 ctx) ts

prettyPrintICanonicalTypeP :: Int -> [NameBinding] -> ICanonicalType -> ShowS
prettyPrintICanonicalTypeP p ctx t = case t of
  CTyPrim p -> showPrimType p
  CTyArr s t -> showParen (p > 2) $ prettyPrintCanonicalTypeP 3 ctx s . showString " -> " . prettyPrintICanonicalTypeP 2 (NAnonymousBind : ctx) t
  CTyRef i | i < length ctx -> case ctx !! i of
               NTyVarBind n -> showString n
               n -> showString "<invalid type variable reference #" . shows i . showString ": " . shows n . showChar '>'
           | otherwise -> showString "<invalid type variable reference #" . shows i . showString ", index out of range>"
  CTyAll (Id name) [] t -> showParen (p > 0) $ showString "forall " . showString name' . showString ". " . prettyPrintICanonicalTypeP 0 (NTyVarBind name' : ctx) t
    where name' = rename (tyVarNames ctx) name
  CTyAll (Id name) b t -> showParen (p > 0) $ showString "forall " . showString name' . showString "<:" . prettyPrintCanonicalTypeP 1 ctx b . showString ". " . prettyPrintICanonicalTypeP 0 (NTyVarBind name' : ctx) t
    where name' = rename (tyVarNames ctx) name
  CTyTuple components -> showParen True $ foldr (.) id $ intersperse (showString ", ") $ map (prettyPrintICanonicalTypeP 0 ctx) components

-- <2>: <primitive> | <variable> | (<0>)
-- <1>: <1> <2> | <1> [ty] | <1> as <type>
-- <0>: \x:t. <0> | ?a. <0> | if <0> then <0> else <0>
prettyPrintTermP :: Int -> [NameBinding] -> Term -> ShowS
prettyPrintTermP p ctx t = case t of
  TmPrim (PVZero) -> showChar '0'
  TmPrim (PVInt x) -> shows x
  TmPrim (PVReal x) -> shows x
  TmPrim (PVImaginary x) -> shows x . showChar 'i'
  TmPrim (PVComplex (x :+ y)) -> showParen True $ shows x . showChar '+' . shows y . showChar 'i'
  TmPrim (PVBool x) -> if x then showString "true" else showString "false"
  TmPrim PVUnit -> showString "unit"
  TmPrim (PVBuiltin f) -> showString $ case f of
    BNegate -> "negate"
    BLogicalNot -> "not"
    BNatToInt -> "natToInt"
    BNatToNNReal -> "natToNNReal"
    BIntToNat -> "intToNat"
    BIntToReal -> "intToReal"
    BNNRealToReal -> "nnrealToReal"
    BRealToComplex -> "realToComplex"
    BMkImaginary -> "mkImaginary"
    BImaginaryToComplex -> "imaginaryToComplex"
    BRealPart -> "realPart"
    BImagPart -> "imagPart"
    BAbs -> "abs"
    BSqrt -> "sqrt"
    BExp -> "exp"
    BExpm1 -> "expm1"
    BLog -> "log"
    BLog1p -> "log1p"
    BSin -> "sin"
    BCos -> "cos"
    BTan -> "tan"
    BSinh -> "sinh"
    BCosh -> "cosh"
    BTanh -> "tanh"
    BAsin -> "asin"
    BAcos -> "acos"
    BAtan -> "atan"
    BAsinh -> "asinh"
    BAcosh -> "acosh"
    BAtanh -> "atanh"
    BAdd -> "add"
    BSub -> "sub"
    BMul -> "mul"
    BDiv -> "div"
    BPow -> "pow"
    BTSubNat -> "tSub"
    BLt -> "lt"
    BLe -> "le"
    BEqual -> "equal"
    BMax -> "max"
    BMin -> "min"
    BIntDiv -> "intDiv"
    BIntMod -> "mod"
    BGcd -> "gcd"
    BLcm -> "lcm"
    BLogicalAnd -> "and"
    BLogicalOr -> "or"
    BIterate -> "iterate"
  TmAbs (Id name) ty body -> showParen (p > 0) $ showChar '\\' . showString name' . showChar ':' . prettyPrintTypeP 1 ctx ty . showString ". " . prettyPrintTermP 0 (NVarBind name' : ctx) body
    where name' = rename (varNames ctx) name
  TmTyAbs (Id name) Nothing body -> showParen (p > 0) $ showString "?" . showString name' . showString ". " . prettyPrintTermP 0 (NTyVarBind name' : ctx) body
    where name' = rename (tyVarNames ctx) name
  TmTyAbs (Id name) (Just b) body -> showParen (p > 0) $ showString "?" . showString name' . showString "<:" . prettyPrintTypeP 1 ctx b . showString ". " . prettyPrintTermP 0 (NTyVarBind name' : ctx) body
    where name' = rename (tyVarNames ctx) name
  TmRef i | i < length ctx -> case ctx !! i of
              NVarBind n -> showString n
              n -> showString "<invalid variable reference #" . shows i . showString ": " . shows n . showChar '>'
          | otherwise -> showString "<invalid variable reference #" . shows i . showString ", index out of range>"
  TmApp u v -> showParen (p > 1) $ prettyPrintTermP 1 ctx u . showChar ' ' . prettyPrintTermP 2 ctx v
  TmTyApp u t -> showParen (p > 1) $ prettyPrintTermP 1 ctx u . showString " [" . prettyPrintTypeP 0 ctx t . showChar ']'
  TmLet (Id name) def body -> showParen (p > 0) $ showString "let " . showString name' . showString " = " . prettyPrintTermP 0 ctx def . showString " in " . prettyPrintTermP 0 (NVarBind name' : ctx) body
    where name' = rename (varNames ctx) name
  TmIf cond then_ else_ -> showParen (p > 0) $ showString "if " . prettyPrintTermP 0 ctx cond . showString " then " . prettyPrintTermP 0 ctx then_ . showString " else " . prettyPrintTermP 0 ctx else_
  TmCoerce x t -> showParen (p > 1) $ prettyPrintTermP 1 ctx x . showString " as " . prettyPrintTypeP 1 ctx t
  TmAlt (Id name) candidates body -> showParen (p > 0) $ showString "for " . showString name' . showString " in " . foldr (.) id (intersperse (showString ", ") $ map (prettyPrintTypeP 1 ctx) candidates) . showString ". " . prettyPrintTermP 0 (NTyVarBind name' : ctx) body
    where name' = rename (tyVarNames ctx) name
  TmTuple components -> showParen True $ foldr (.) id $ intersperse (showString ", ") $ map (prettyPrintTermP 0 ctx) components
  TmProj tuple i -> showParen (p > 1) $ prettyPrintTermP 1 ctx tuple . showString " @" . shows i
  TmCoherentTuple components -> showString "<coherent tuple: " . (foldr (.) id $ intersperse (showString ", ") $ map (prettyPrintTermP 0 ctx) components) . showString ">"

prettyPrintTerm :: Term -> String
prettyPrintTerm t = prettyPrintTermP 0 [] t ""

max_stars_repo_name: minoki/polestar

max_issues_repo_name: minoki/polestar

max_forks_repo_name: minoki/polestar


text:
-- file: Astro/Kepler.hs
-- Various forms of Kepler's equation with efficient solvers

module Astro.Kepler ( eKepler
		    , gKepler
		    , solveEKepler
		    , solveGKepler
		    , keplerPropagate
		    )
	where

import Astro.Stumpff
import Astro.Misc

import Complex
import Debug.Trace
import Numeric.GSL
import Numeric.LinearAlgebra
import Text.Printf (printf)

-- Shorthands and such
type CD = Complex Double
type VD = Vector Double

-- Elliptic Kepler's equation
eKepler :: (Floating a) => a -> a -> a
eKepler e ea = ea - e * sin ea

-- General Variable Formulation of Kepler's eq
gKepler :: Double -> Double -> Double -> Double -> Double -> Double
gKepler r0 eta0 zeta0 beta x =
	r0 * x + eta0 * stumpff_G 2 beta x + zeta0 * stumpff_G 3 beta x

-- Solve Elliptic Kepler's equation iteratively
solveEKepler :: (Floating a) => a -> a -> a
solveEKepler e m =
	iterate 0 m m
	where
	--iterate n ea m | trace ("n,E: " ++ show n ++ " " ++ show ea) False = undefined
	iterate n ea m
		| ea' == ea || n > maxit = ea'
		| otherwise = iterate (n+1) ea' m
		where
		    fd = 1.0 - e * cos ea
		    fdd = e * sin ea
		    f = eKepler e ea - m
		    ea' = ea - f / sqrt (fd^2 - f * fdd)
	maxit = 10

-- Solve General Kepler's eq iteratively
--solveGKepler :: Double -> Double -> Double -> Double -> Double -> Double -> Double
solveGKepler r0 eta0 zeta0 beta x0 t (tol,maxIt) =
	if abs lastErr > tol
		then head sol `debug`
			("solveGKepler: iter. incomp., error: "
				++(show lastErr))
		else head sol
	where
		(sol,path) = root Hybrids tol maxIt f [x0]
		lastErr = (last.last) $ toLists path
		f vx = [gKepler r0 eta0 zeta0 beta (head vx) - t]

-- Solve r, v given r0, v0, gravitational parameter, delta_t
--keplerPropagate ::  VD -> VD -> Double -> (Double, Int) -> (VD,VD)
keplerPropagate dt (tol,maxIt) rv vv =
	(scale f rv + scale g vv , scale fd rv + scale gd vv)
	where (r0, v0) = (norm rv, norm vv)
	      tau = dt -- working in units where G(m_1+m_2) = 1
	      eta0 = rv <.> vv
	      beta = 2/r0 - v0^2
	      zeta0 = 1 - beta * r0
	      x = solveGKepler r0 eta0 zeta0 beta 0.0 tau (tol,maxIt)
	      r = r0 + eta0 * stumpff_G 1 beta x + zeta0 * stumpff_G 2 beta x
	      f = 1 - stumpff_G 2 beta x / r0
	      g = tau - stumpff_G 3 beta x
	      fd = - stumpff_G 1 beta x / (r * r0)
	      gd = 1 - stumpff_G 2 beta x / r

max_stars_repo_name: sageh/AstroHaskell

max_issues_repo_name: sageh/AstroHaskell

max_forks_repo_name: sageh/AstroHaskell


text:
module Main where

import Numeric.NLOPT
-- import Numeric.LinearAlgebra (dot)
import Numeric.LinearAlgebra (fromList, toList)

main :: IO ()
main = do
    let -- objf x = x `dot` x + 22
        -- objf x = sum (toList x)
        objf x = let [x1, x2, x3] = toList x
                 in 3 * x1 * x1 + 2 * x1 * x2 + x1 * x3
                    + 2.5 * x2 * x2 + 2 * x2 * x3 + 2 * x3 * x3
                    - 8 * x1 - 3 * x2 - 3 * x3

        stop = ObjectiveRelativeTolerance 1e-9 :| []
        -- algorithm = SBPLX objf [] Nothing
        algorithm = NELDERMEAD objf [] Nothing
        subproblem = LocalProblem 3 stop algorithm
        x0 = fromList [0, 0, 0]

        solU = minimizeLocal subproblem x0

    print solU

    -- define constraint function:
    let -- constraintf x = sum (toList x) - 1.0
        -- constraintf x = let [x1, x2] = toList x
        --                 in x1 * x1 + x2 * x2 - 2
        c1f x = let [x1, _, x3] = toList x
                in x1 + x3 - 3
        c2f x = let [_, x2, x3] = toList x
                in x2 + x3

    -- define constraint object to pass to the algorithm:
    let -- constraint = EqualityConstraint (Scalar constraintf) 1e-6
        c1 = EqualityConstraint (Scalar c1f) 1e-6
        c2 = EqualityConstraint (Scalar c2f) 1e-6
        -- problem = AugLagProblem [constraint] [] (AUGLAG_EQ_LOCAL subproblem)
        problem = AugLagProblem [c1, c2] [] (AUGLAG_EQ_LOCAL subproblem)

        -- x = (2, -1, 1)
        solC = minimizeAugLag problem x0

    print solC

    putStrLn "-- Done."

max_stars_repo_name: cbpark/hmatrix-nlopt-test

max_issues_repo_name: cbpark/hmatrix-nlopt-test

max_forks_repo_name: cbpark/hmatrix-nlopt-test


text:
module AI.Learning.Perceptron where

import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Util

perceptronPredict :: Vector Double -> Matrix Double -> Vector Double
perceptronPredict weights x = step (x <> weights)

perceptronCost :: Vector Double -> Vector Double -> Vector Double -> Double
perceptronCost y yhat = undefined

----------------------
-- Gradient Descent --
----------------------

gradientDescent :: (Vector Double -> Vector Double) -- f
                -> Vector Double                    -- x0
                -> Double                           -- alpha
                -> Double                           -- tol
                -> Vector Double
gradientDescent g x0 alpha tol = go x0 (fun x0)
    where
        go x x' = if converged x x'
                    then x'
                    else go x' (fun x')
        converged a b = norm b / norm a - 1 < tol
        fun x = gradientDescentStep g alpha x

gradientDescentStep :: (Vector Double -> Vector Double) -> Double -> Vector Double -> Vector Double
gradientDescentStep g alpha x = x - alpha `scale` g x

max_stars_repo_name: cagix/aima-haskell


text:
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DeriveAnyClass        #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-|
Module      : Grenade.Layers.Logit
Description : Exponential linear unit layer
Copyright   : (c) Huw Campbell, 2016-2017
License     : BSD2
Stability   : experimental
-}
module Grenade.Layers.Elu (
    Elu (..)
  ) where

import           Control.DeepSeq              (NFData)
import           Data.Serialize
import           GHC.Generics                 (Generic)

import           GHC.TypeLits
import           Grenade.Core

import qualified Numeric.LinearAlgebra.Static as LAS

-- | An exponential linear unit.
--   A layer which can act between any shape of the same dimension, acting as a
--   diode on every neuron individually.
data Elu = Elu
  deriving (Generic, NFData, Show)

instance UpdateLayer Elu where
  type Gradient Elu = ()
  runUpdate _ _ _ = Elu

instance RandomLayer Elu where
  createRandomWith _ _ = return Elu

instance Serialize Elu where
  put _ = return ()
  get = return Elu

instance ( KnownNat i) => Layer Elu ('D1 i) ('D1 i) where
  type Tape Elu ('D1 i) ('D1 i) = LAS.R i

  runForwards _ (S1D y) = (y, S1D (elu y))
    where
      elu = LAS.dvmap (\a -> if a <= 0 then exp a - 1 else a)
  runBackwards _ y (S1D dEdy) = ((), S1D (elu' y * dEdy))
    where
      elu' = LAS.dvmap (\a -> if a <= 0 then exp a else 1)

instance (KnownNat i, KnownNat j) => Layer Elu ('D2 i j) ('D2 i j) where
  type Tape Elu ('D2 i j) ('D2 i j) = S ('D2 i j)

  runForwards _ (S2D y) = (S2D y, S2D (elu y))
    where
      elu = LAS.dmmap (\a -> if a <= 0 then exp a - 1 else a)
  runBackwards _ (S2D y) (S2D dEdy) = ((), S2D (elu' y * dEdy))
    where
      elu' = LAS.dmmap (\a -> if a <= 0 then exp a else 1)

instance (KnownNat i, KnownNat j, KnownNat k) => Layer Elu ('D3 i j k) ('D3 i j k) where

  type Tape Elu ('D3 i j k) ('D3 i j k) = S ('D3 i j k)

  runForwards _ (S3D y) = (S3D y, S3D (elu y))
    where
      elu = LAS.dmmap (\a -> if a <= 0 then exp a - 1 else a)
  runBackwards _ (S3D y) (S3D dEdy) = ((), S3D (elu' y * dEdy))
    where
      elu' = LAS.dmmap (\a -> if a <= 0 then exp a else 1)

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

instance GNum Elu where
  _ |* Elu = Elu
  _ |+ Elu = Elu
  gFromRational _ = Elu


text:
module Main (main) where

import           Control.Monad
import           System.IO
import           Text.Printf (printf)

import           Data.Number.LogFloat (LogFloat)
import qualified Data.Number.LogFloat as LF
import           Statistics.Distribution
import           Statistics.Distribution.Binomial

-- | Compute the probability of seeing more than @k@ blocks in @n@ slots
moreThanK ::
     Double  -- ^ Active slot coefficient
  -> Int     -- ^ Security parameter (@k@)
  -> Int     -- ^ Number of slots
  -> LogFloat
moreThanK f k n =
    LF.sum [LF.logToLogFloat $ logProbability d i | i <- [k + 1 .. n]]
  where
    d :: BinomialDistribution
    d = binomial n f

defaultS ::
     Double  -- ^ Active slot coefficient
  -> Int     -- ^ Security parameter (@k@)
  -> Int
defaultS f k = floor (fromIntegral k / f) `div` 4

main :: IO ()
main = do
    forM_ [s .. 4 * s] $ \n -> do
      putStrLn $ show n ++ "\t" ++ showLogFloat (moreThanK f k n)
      hFlush stdout
  where
    f = 0.05
    k = 2160
    s = defaultS f k

{-------------------------------------------------------------------------------
  LogFloat util
-------------------------------------------------------------------------------}

showLogFloat :: LogFloat -> String
showLogFloat lf = printf "%6.4f * 10 ^ %d" m e
  where
    (m, e) = logFloatToScientific lf

logFloatToScientific :: LogFloat -> (Double, Int)
logFloatToScientific lf = (m, e)
  where
    l :: Double
    l = LF.logFromLogFloat lf

    e :: Int
    e = floor $ l / log10

    m :: Double
    m = exp $ l - log10 * fromIntegral e

    log10 :: Double
    log10 = log 10

max_stars_repo_name: RyanGlScott/ouroboros-network

max_issues_repo_name: RyanGlScott/ouroboros-network

max_forks_repo_name: RyanGlScott/ouroboros-network


text:
module FFT.Orig (
    fft,
    dft,
) where
import Data.Complex
import FFT.Samples
import System.Environment
import Control.Parallel
import Strategies

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

-- Discrete Fourier Transform -- O(n^2)
dft :: [Complex Float] -> [Complex Float]
dft xs = [ sum [ xs!!j * tw n (j*k) | j <- [0..n']] | k <- [0..n']]
  where
    n = length xs
    n' = n-1

-- Fast Fourier Transform

-- In case you are wondering, this is the Decimation in Frequency (DIF)
-- radix 2 Cooley-Tukey FFT

fft :: [Complex Float] -> [Complex Float]
fft [a] = [a]
fft as = interleave ls rs
  where
    (cs,ds) = bflyS as
    ls = fft cs
    rs = fft ds

interleave [] bs = bs
interleave (a:as) bs = a : interleave bs as

bflyS :: [Complex Float] -> ([Complex Float], [Complex Float])
bflyS as = (los,rts)
  where
    (ls,rs) = halve as
    los = zipWith (+) ls rs
    ros = zipWith (-) ls rs
    rts = zipWith (*) ros [tw (length as) i | i <- [0..(length ros) - 1]]

-- split the input into two halves
halve as = splitAt n' as
  where
    n' = div (length as + 1) 2

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 TypeFamilies #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Signal.Internal
-- Copyright   :  (c) Alexander Vivian Hugh McPhail 2010, 2014, 2015, 2016
-- License     :  BSD3
--
-- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
-- Stability   :  provisional
-- Portability :  uses FFI
--
-- low-level interface
--
-----------------------------------------------------------------------------

module Numeric.Signal.Internal (
                Convolvable(..),
                Filterable(..),
                freqz,
                pwelch,
                hilbert
                ) where

import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Devel

--import Numeric.LinearAlgebra.Linear

import qualified Numeric.GSL.Fourier as F
import Foreign

--import Data.Complex
import Foreign.C.Types

import Prelude hiding(filter)

import System.IO.Unsafe(unsafePerformIO)

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

infixr 1 #
a # b = applyRaw a b
{-# INLINE (#) #-}

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

type PD = Ptr Double
type PC = Ptr (Complex Double)
type PF = Ptr Float

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

class Convolvable a where
    -- | convolve two containers, output is the size of the second argument, no zero padding
    convolve :: a -> a -> a

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

class (Storable a, Container Vector a, Num (Vector a)
      , Convert a, Floating (Vector a), RealElement a
      , Num a)
       => Filterable a where
    -- | convert from Vector Double
    fromDouble :: Vector Double -> Vector a
--       b ~ ComplexOf a, Container Vector b, Convert b) => Filterable a where
    -- | filter a signal
    filter_ :: Vector a -- ^ zero coefficients
            -> Vector a -- ^ pole coefficients
            -> Vector a -- ^ input signal
            -> Vector a -- ^ output signal
    -- | coefficients of a Hamming window
    hamming_ :: Int           -- ^ length
           -> Vector a -- ^ the Hamming coeffficents
    -- | the complex power : real $ v * (conj v)
    complex_power_ :: Vector (Complex Double) -- ^ input
                  -> Vector a                       -- ^ output
     -- | resample, take one sample every n samples in the original
    downsample_ :: Int -> Vector a -> Vector a
    -- | the difference between consecutive elements of a vector
    deriv_ :: Vector a -> Vector a
    -- | unwrap the phase of signal (input expected to be within (-pi,pi)
    unwrap_ :: Vector a -> Vector a
    -- | evaluate a real coefficient polynomial for complex arguments
    polyEval_ :: Vector a           -- ^ the real coefficients
         -> Vector (Complex Double) -- ^ the points at which to be evaluated
         -> Vector (Complex Double) -- ^ the values
    -- | the cross covariance of two signals
    cross_covariance_ :: Int            -- ^ maximum delay
                     -> Vector a       -- ^ time series
                     -> Vector a       -- ^ time series
                     -> (a,a,Vector a) -- ^ (sd_x,sd_y,cross_covariance)
    -- | the cumulative sum of a signal
    cumulative_sum_ :: Vector a    -- ^ time series
                    -> Vector a    -- ^ result

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

instance Convolvable (Vector Double) where
    convolve x y = fst $ fromComplex $ F.ifft $ (F.fft (complex x) * F.fft (complex y))
--    convolve = convolve_vector_double

convolve_vector_double c a = unsafePerformIO $ do
                             r <- createVector (size a)
                             (c # a # r # id) signal_vector_double_convolve #| "signalDoubleConvolve"
                             return r

foreign import ccall "signal-aux.h vector_double_convolve" signal_vector_double_convolve :: CInt -> PD -> CInt -> PD -> CInt -> PD -> IO CInt

instance Convolvable (Vector Float) where
    convolve x y = single $ fst $ fromComplex $ F.ifft $ (F.fft (complex $ double x) * F.fft (complex $ double y))
--    convolve = convolve_vector_double

convolve_vector_float c a = unsafePerformIO $ do
                             r <- createVector (size a)
                             (c # a # r # id ) signal_vector_float_convolve #| "signalFloatConvolve"
                             return r

foreign import ccall "signal-aux.h vector_float_convolve" signal_vector_float_convolve :: CInt -> PF -> CInt -> PF -> CInt -> PF -> IO CInt

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

instance Convolvable (Vector (Complex Double)) where
    convolve x y = F.ifft $ (F.fft x * F.fft y)
--    convolve = convolve_vector_complex

convolve_vector_complex c a = unsafePerformIO $ do
                              r <- createVector (size a)
                              (c # a # r # id) signal_vector_complex_convolve #| "signalComplexConvolve"
                              return r

foreign import ccall "signal-aux.h vector_complex_convolve" signal_vector_complex_convolve :: CInt -> PC -> CInt -> PC -> CInt -> PC -> IO CInt

instance Convolvable (Vector (Complex Float)) where
    convolve x y = single $ F.ifft $ (F.fft (double x) * F.fft (double y))

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

instance Filterable Double where
    fromDouble = id
    filter_ = filterD
    hamming_ = hammingD
    complex_power_ = complex_powerD
    downsample_ = downsampleD
    deriv_ = derivD
    unwrap_ = unwrapD
    polyEval_ = polyEval
    cross_covariance_ = crossCovarianceD
    cumulative_sum_ = cumSumD

instance Filterable Float where
    fromDouble = single
    filter_ = filterF
    hamming_ = hammingF
    complex_power_ = complex_powerF
    downsample_ = downsampleF
    deriv_ = derivF
    unwrap_ = unwrapF
    polyEval_ c = polyEval (double c)
    cross_covariance_ = crossCovarianceF
    cumulative_sum_ = cumSumF

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

-- | filters the signal
filterD :: Vector Double -- ^ zero coefficients
       -> Vector Double -- ^ pole coefficients
       -> Vector Double -- ^ input signal
       -> Vector Double -- ^ output signal
filterD l k v = unsafePerformIO $ do
               r <- createVector (size v)
               (l # k # v # r # id) signal_filter_double #| "signalFilter"
               return r

foreign import ccall "signal-aux.h filter_double" signal_filter_double :: CInt -> PD -> CInt -> PD -> CInt -> PD -> CInt -> PD -> IO CInt

-- | filters the signal
filterF :: Vector Float -- ^ zero coefficients
       -> Vector Float -- ^ pole coefficients
       -> Vector Float -- ^ input signal
       -> Vector Float -- ^ output signal
filterF l k v = unsafePerformIO $ do
               r <- createVector (size v)
               (l # k # v # r # id) signal_filter_float #| "signalFilter"
               return r

foreign import ccall "signal-aux.h filter_float" signal_filter_float :: CInt -> PF -> CInt -> PF -> CInt -> PF -> CInt -> PF -> IO CInt

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

-- | Hilbert transform with original vector as real value, transformed as imaginary
hilbert :: Vector Double -> Vector (Complex Double)
hilbert v = unsafePerformIO $ do
            let r = complex v
            -- could use (complex v) to make a complex vector in haskell rather than C
            (r # id) signal_hilbert #| "hilbert"
            return r

foreign import ccall "signal-aux.h hilbert" signal_hilbert :: CInt -> PC -> IO CInt

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

-- | Welch (1967) power spectrum density using periodogram/FFT method
pwelch :: Int            -- ^ window size (multiple of 2)
       -> Vector Double  -- ^ input signal
       -> Vector Double  -- ^ power density
pwelch w v = unsafePerformIO $ do
             let r = konst 0.0 ((w `div` 2) + 1)
             (complex v # r # id) (signal_pwelch $ fromIntegral w) #| "pwelch"
             return r

foreign import ccall "signal-aux.h pwelch" signal_pwelch :: CInt -> CInt -> PC -> CInt -> PD -> IO CInt

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

-- | coefficients of a Hamming window
hammingD :: Int           -- ^ length
        -> Vector Double -- ^ the Hamming coeffficents
hammingD l
    | l == 1          = konst 1.0 1
    | otherwise       = unsafePerformIO $ do
                        r <- createVector l
                        (r # id) signal_hamming_double #| "Hamming"
                        return r

foreign import ccall "signal-aux.h hamming_double" signal_hamming_double :: CInt -> PD -> IO CInt

-- | coefficients of a Hamming window
hammingF :: Int           -- ^ length
        -> Vector Float -- ^ the Hamming coeffficents
hammingF l
    | l == 1          = konst 1.0 1
    | otherwise       = unsafePerformIO $ do
                        r <- createVector l
                        (r # id) signal_hamming_float #| "Hamming"
                        return r

foreign import ccall "signal-aux.h hamming_float" signal_hamming_float :: CInt -> PF -> IO CInt

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

-- | determine the frequency response of a filter
{-freqz :: (Filterable a, Storable a, Container Vector a, Convert a, RealElement a,
         DoubleOf a ~ DoubleOf (RealOf b), RealElement c, c ~ DoubleOf a, c ~ DoubleOf (RealOf b),
         b ~ Complex a, b ~ ComplexOf a, Convert b, Container Vector b,
         Container Vector c, Convert c, b ~ ComplexOf a, b ~ ComplexOf c)
        ⇒ Vector a     -- ^ zero coefficients
      -> Vector a       -- ^ pole coefficients
      -> Vector a       -- ^ points (between 0 and 2*pi)
      -> Vector a       -- ^ response
-}
freqz :: (Filterable a, Complex Double ~ ComplexOf (DoubleOf a)
        ,Filterable (DoubleOf a)) =>
        Vector a       -- ^ zero coefficients
      -> Vector a       -- ^ pole coefficients
      -> Vector a       -- ^ points (between 0 and 2*pi)
      -> Vector a       -- ^ response
freqz b a w = let k = max (size b) (size a)
                  hb = polyEval_ (postpad b k) (exp (scale (0 :+ 1) ((complex $ double w))))
                  ha = polyEval_ (postpad a k) (exp (scale (0 :+ 1) ((complex $ double w))))
              in complex_power_ (hb / ha)

postpad v n = let d = size v
              in if d < n then vjoin [v,(konst 0.0 (n-d))]
              else v

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

-- | evaluate a real coefficient polynomial for complex arguments
polyEval :: Vector Double           -- ^ the real coefficients
         -> Vector (Complex Double) -- ^ the points at which to be evaluated
         -> Vector (Complex Double) -- ^ the values
polyEval c z = unsafePerformIO $ do
               r <- createVector (size z)
               (c # z # r # id) signal_real_poly_complex_eval #| "polyEval"
               return r

foreign import ccall "signal-aux.h real_poly_complex_eval" signal_real_poly_complex_eval :: CInt -> PD -> CInt -> PC -> CInt -> PC -> IO CInt

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

-- | the complex power : real $ v * (conj v)
complex_powerD :: Vector (Complex Double) -- ^ input
              -> Vector Double           -- ^ output
complex_powerD v = unsafePerformIO $ do
                  r <- createVector (size v)
                  (v # r # id) signal_complex_power_double #| "complex_power"
                  return r

foreign import ccall "signal-aux.h complex_power_double" signal_complex_power_double :: CInt -> PC -> CInt -> PD -> IO CInt

-- | the complex power : real $ v * (conj v)
complex_powerF :: Vector (Complex Double) -- ^ input
              -> Vector Float             -- ^ output
complex_powerF v = unsafePerformIO $ do
                  r <- createVector (size v)
                  (v # r # id) signal_complex_power_float #| "complex_power"
                  return r

foreign import ccall "signal-aux.h complex_power_float" signal_complex_power_float :: CInt -> PC -> CInt -> PF -> IO CInt

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

-- | resample, take one sample every n samples in the original
downsampleD :: Int -> Vector Double -> Vector Double
downsampleD n v = unsafePerformIO $ do
               r <- createVector (size v `div` n)
               (v # r # id) (signal_downsample_double $ fromIntegral n) #| "downsample"
               return r

foreign import ccall "signal-aux.h downsample_double" signal_downsample_double :: CInt -> CInt -> PD -> CInt -> PD -> IO CInt

-- | resample, take one sample every n samples in the original
downsampleF :: Int -> Vector Float -> Vector Float
downsampleF n v = unsafePerformIO $ do
               r <- createVector (size v `div` n)
               (v # r # id) (signal_downsample_float $ fromIntegral n) #| "downsample"
               return r

foreign import ccall "signal-aux.h downsample_float" signal_downsample_float :: CInt -> CInt -> PF -> CInt -> PF -> IO CInt

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

-- | the difference between consecutive elements of a vector
derivD :: Vector Double -> Vector Double
derivD v = unsafePerformIO $ do
          r <- createVector (size v - 1)
          (v # r # id) (signal_diff_double) #| "diff"
          return r

foreign import ccall "signal-aux.h vector_diff_double" signal_diff_double :: CInt -> PD -> CInt -> PD -> IO CInt

-- | the difference between consecutive elements of a vector
derivF :: Vector Float -> Vector Float
derivF v = unsafePerformIO $ do
          r <- createVector (size v - 1)
          (v # r # id) (signal_diff_float) #| "diff"
          return r

foreign import ccall "signal-aux.h vector_diff_float" signal_diff_float :: CInt -> PF -> CInt -> PF -> IO CInt

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

-- | unwrap the phase of signal (input expected to be within (-pi,pi)
unwrapD :: Vector Double -> Vector Double
unwrapD v = unsafePerformIO $ do
           r <- createVector $ size v
           (v # r # id) signal_unwrap_double #| "unwrap"
           return r

foreign import ccall "signal-aux.h unwrap_double" signal_unwrap_double :: CInt -> PD -> CInt -> PD -> IO CInt

-- | unwrap the phase of signal (input expected to be within (-pi,pi)
unwrapF :: Vector Float -> Vector Float
unwrapF v = unsafePerformIO $ do
           r <- createVector $ size v
           (v # r # id) signal_unwrap_float #| "unwrap"
           return r

foreign import ccall "signal-aux.h unwrap_float" signal_unwrap_float :: CInt -> PF -> CInt -> PF -> IO CInt

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

-- | compute the cross covariance of two signals
crossCovarianceD :: Int -> Vector Double -> Vector Double -> (Double,Double,Vector Double)
crossCovarianceD l x y = unsafePerformIO $ do
                           r <- createVector (2*l)
                           alloca $ \sx ->
                               alloca $ \sy -> do
                                 (x # y # r # id) (signal_cross_covariance_double (fromIntegral l) sx sy) #| "cross_covariance"
                                 sx' <- peek sx
                                 sy' <- peek sy
                                 return (sx',sy',r)

foreign import ccall "signal-aux.h cross_covariance_double"
        signal_cross_covariance_double :: CInt -> PD -> PD -> CInt -> PD -> CInt
                                       -> PD -> CInt -> PD -> IO CInt

-- | compute the cross covariance of two signals
crossCovarianceF :: Int -> Vector Float -> Vector Float -> (Float,Float,Vector Float)
crossCovarianceF l x y = unsafePerformIO $ do
                           r <- createVector (2*l)
                           alloca $ \sx ->
                               alloca $ \sy -> do
                                 (x # y # r # id) (signal_cross_covariance_float (fromIntegral l) sx sy) #| "cross_covariance"
                                 sx' <- peek sx
                                 sy' <- peek sy
                                 return (sx',sy',r)

foreign import ccall "signal-aux.h cross_covariance_float"
        signal_cross_covariance_float :: CInt -> PF -> PF -> CInt -> PF -> CInt
                                       -> PF -> CInt -> PF -> IO CInt

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

cumSumD :: Vector Double -> Vector Double
cumSumD v = unsafePerformIO $ do
              r <- createVector (size v)
              (v # r # id) signal_cum_sum_double #| "cumSumD"
              return r

cumSumF :: Vector Float -> Vector Float
cumSumF v = unsafePerformIO $ do
              r <- createVector (size v)
              (v # r # id) signal_cum_sum_float #| "cumSumF"
              return r

foreign import ccall "signal-aux.h cum_sum_double"
        signal_cum_sum_double :: CInt -> PD -> CInt -> PD -> IO CInt
foreign import ccall "signal-aux.h cum_sum_float"
        signal_cum_sum_float :: CInt -> PF -> CInt -> PF -> IO CInt

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

max_stars_repo_name: amcphail/hsignal

max_issues_repo_name: amcphail/hsignal

max_forks_repo_name: amcphail/hsignal


text:
{-# LANGUAGE OverloadedStrings #-}

module Main where

import Data.Csv
import Numeric.LinearAlgebra
import qualified Data.ByteString.Lazy as BL
import qualified Data.Vector as V

import Lib

-- Use Iris dataset to test functionality

data Iris = Iris
  { sepalLength :: !Double,
    sepalWidth  :: !Double,
    petalLength :: !Double,
    petalWidth  :: !Double,
    className   :: !String
  }

encodeClassName :: Iris -> [Double]
encodeClassName i = case className i of
  "Iris-setosa"     -> [1,0,0]
  "Iris-versicolor" -> [0,1,0]
  "Iris-virginica"  -> [0,0,1]

instance FromNamedRecord Iris where
  parseNamedRecord r = Iris <$> r .: "sepalLength"
    <*> r .: "sepalWidth"
    <*> r .: "petalLength"
    <*> r .: "petalWidth"
    <*> r .: "className"

irisToInput :: Iris -> [Double]
irisToInput iris = map ($ iris) [sepalLength,sepalWidth,petalLength,petalWidth]

-- Hyperparameters
regularizationTerm=0.01
learningRate=0.02
layerSizes=[4,5,3]

-- Define Model
initializerScale=0.0005
initialBiases = biasesWithInitializer initializerScale layerSizes
initialWeights = weightsWithInitializer initializerScale layerSizes
initialParameters = zip initialWeights initialBiases
activationFunctions = [sigmoid,sigmoid]
layers = zipWith getLayer initialParameters activationFunctions

-- TODO: Write function to get validation set performance

main = do
  irisTrainingData <- BL.readFile "data/training.csv"
  case decodeByName irisTrainingData of
    Left err -> putStrLn err
    Right (_,vTraining) -> do
      let
        irisTrainingList = map irisToInput $ V.toList vTraining
        targetOutputList = map encodeClassName $ V.toList vTraining
        irisTrainingInputMatrix = fromLists irisTrainingList
        irisTrainingTargetOutputMatrix = fromLists targetOutputList
        updatedParameters = runEpochs 3000 learningRate initializerScale irisTrainingInputMatrix irisTrainingTargetOutputMatrix layers
        updatedLayers = zipWith getLayer updatedParameters activationFunctions
      irisTestData <- BL.readFile "data/test.csv"
      case decodeByName irisTestData of
        Left err -> putStrLn err
        Right (_,vTest) -> do
          let
            irisTestList = map irisToInput $ V.toList vTest
            targetOutputList = map encodeClassName $ V.toList vTest
            irisTestInputMatrix = fromLists irisTestList
            irisTestTargetOutputMatrix = fromLists targetOutputList
            testActivations = forwardPropogate irisTestInputMatrix updatedLayers
            outputs = (postActivation . last) testActivations :: Matrix Double
            targetOutputRows = toRows irisTestTargetOutputMatrix
            outputRows = toRows outputs
            targetArgmaxList = map maxIndex targetOutputRows
            outputArgmaxList = map maxIndex outputRows
            zippedArgmaxLists = zip targetArgmaxList outputArgmaxList
            total = length zippedArgmaxLists
            correct = length [(a,b) | (a,b) <- zippedArgmaxLists,a == b]
          print total
          print correct
          print outputs
          print irisTestTargetOutputMatrix

max_stars_repo_name: dmurphs/hsnet

max_issues_repo_name: dmurphs/hsnet

max_forks_repo_name: dmurphs/hsnet


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE UndecidableInstances #-}
{-# OPTIONS_GHC -Wno-orphans #-}
{- Module      : Text.Loquate.Instances
   Description : Default Loquate Instances
-}
module Text.Loquate.Instances
  (
  ) where

--import qualified Data.ByteString as BS
--import qualified Data.ByteString.Lazy as LBS
import qualified Control.Exception.Base as E
import           Data.Array.Unboxed (Array, IArray, Ix, UArray, elems)
import           Data.Complex (Complex(..))
import           Data.Fixed
import           Data.Foldable
import           Data.Int (Int16, Int32, Int64, Int8)
import           Data.List.NonEmpty (NonEmpty)
import           Data.Ratio (Ratio, denominator, numerator)
import qualified Data.Text as T
import qualified Data.Text.Lazy as LT
import           Data.Text.Prettyprint.Doc (Pretty(..))
import           Data.Typeable
import           Data.Version (Version, showVersion)
import           Data.Void (Void)
import           Data.Word (Word16, Word32, Word64, Word8)
import           GHC.Stack
import           Numeric.Natural
import           Text.Loquate.Class
import           Text.Loquate.Doc

-- misc types

instance Lang l => Loquate l () where loq _ = pretty
instance Lang l => Loquate l Bool where loq _ = pretty
instance Lang l => Loquate l Version where loq _ = pretty . showVersion
instance Lang l => Loquate l Void where loq _ = pretty

-- textual types

instance Lang l => Loquate l Char where loq _ = pretty
instance {-# OVERLAPS #-} Lang l => Loquate l String where loq _ = pretty
-- instance Loquate l BS.ByteString where loq _ = pretty
-- instance Loquate l LBS.ByteString where loq _ = pretty
instance Lang l => Loquate l T.Text where loq _ = pretty
instance Lang l => Loquate l LT.Text where loq _ = pretty

-- numeric types

instance Lang l => Loquate l Double where loq _ = pretty
instance Lang l => Loquate l Float where loq _ = pretty
instance (Lang l, Typeable t, HasResolution t) => Loquate l (Fixed t) where loq _ = viaShow
-- | Integral types should loquate via Integer
instance Lang l => Loquate l Integer where loq _ = pretty
instance (Lang l, Loquate l Integer) => Loquate l Natural where loq l = loq l . toInteger

instance (Lang l, Loquate l Integer) => Loquate l Int where loq l = loq l . toInteger
instance (Lang l, Loquate l Integer) => Loquate l Int8 where loq l = loq l . toInteger
instance (Lang l, Loquate l Integer) => Loquate l Int16 where loq l = loq l . toInteger
instance (Lang l, Loquate l Integer) => Loquate l Int32 where loq l = loq l . toInteger
instance (Lang l, Loquate l Integer) => Loquate l Int64 where loq l = loq l . toInteger

instance (Lang l, Loquate l Integer) => Loquate l Word where loq l = loq l . toInteger
instance (Lang l, Loquate l Integer) => Loquate l Word8 where loq l = loq l . toInteger
instance (Lang l, Loquate l Integer) => Loquate l Word16 where loq l = loq l . toInteger
instance (Lang l, Loquate l Integer) => Loquate l Word32 where loq l = loq l . toInteger
instance (Lang l, Loquate l Integer) => Loquate l Word64 where loq l = loq l . toInteger

-- polymorphic numeric types

instance Loquate l t => Loquate l (Complex t) where
  loq l (a :+ b) = loq l a <+> "+" <+> loq l b <> "i"

instance Loquate l t => Loquate l (Ratio t) where
  loq l x = loq l (numerator x) <> "/" <> loq l (denominator x)

-- polymorphic containers in dependencies of loquacious

instance Loquate l t => Loquate l (Maybe t) where
  loq _ Nothing  = mempty
  loq l (Just x) = loq l x

instance (Loquate l t, Loquate l s) => Loquate l (t, s) where
  loq l (a, b) = align $ tupled [loq l a, loq l b]

instance (Loquate l t, Loquate l s, Loquate l r) => Loquate l (t, s, r) where
  loq l (a, b, c) = align $ tupled [loq l a, loq l b, loq l c]

instance (Loquate l t) => Loquate l [t] where
  loq l xs = align . list $ loq l <$> xs

instance (Loquate l t) => Loquate l (NonEmpty t) where
  loq l xs = align . list $ loq l <$> toList xs

instance (Typeable i, Loquate l t) => Loquate l (Array i t) where
  loq l xs = align . list $ loq l <$> toList xs

instance (Loquate l t, IArray UArray t, Ix i, Typeable i) => Loquate l (UArray i t) where
  loq l xs = align . list $ loq l <$> elems xs

instance (Loquate l t, Loquate l s) => Loquate l (Either t s) where
  loq l (Left x)  = loq l x
  loq l (Right x) = loq l x

-- showing types rather than values

instance Lang l => Loquate l TypeRep where loq _ = viaShow

instance (Lang l, Typeable a) => Loquate l (Proxy a) where
  loq l _ = loq l $ typeRep (Proxy :: Proxy (Proxy a))

instance (Lang l, Typeable a, Typeable b) => Loquate l (a -> b) where
  loq l f = loq l $ typeOf f

instance (Lang l, Typeable a, Typeable b) => Loquate l (a :~: b) where
  loq l _ = loq l (typeRep (Proxy :: Proxy a)) <+> ":~:" <+> loq l (typeRep (Proxy :: Proxy b))

instance (Lang l, Typeable a, Typeable b) => Loquate l (a :~~: b) where
  loq l _ = loq l (typeRep (Proxy :: Proxy a)) <+> ":~~:" <+> loq l (typeRep (Proxy :: Proxy b))

-- call stacks!

instance Lang l => Loquate l CallStack where
  loq l = withFrozenCallStack ( align . vsep . fmap loqStack . getCallStack )
    where loqStack (func, srcLoc) = fromString func <+> "←" <+> loq l srcLoc

instance Lang l => Loquate l SrcLoc where
  loq l SrcLoc{..} = locFile <+> parens locPackage
    where
      locFile = fromString srcLocFile <> colon <> loq l srcLocStartLine <> colon <> loq l srcLocStartCol
      locPackage = fromString srcLocPackage <> colon <> fromString srcLocModule

-- Exception types

instance Lang l => Loquate l E.AllocationLimitExceeded where loq _ e = viaShow e
instance Lang l => Loquate l E.ArithException where loq _ e = viaShow e
instance Lang l => Loquate l E.ArrayException where loq _ e = viaShow e
instance Lang l => Loquate l E.AssertionFailed where loq _ e = viaShow e
instance Lang l => Loquate l E.AsyncException where loq _ e = viaShow e
instance Lang l => Loquate l E.BlockedIndefinitelyOnMVar where loq _ e = viaShow e
instance Lang l => Loquate l E.BlockedIndefinitelyOnSTM where loq _ e = viaShow e
instance Lang l => Loquate l E.CompactionFailed where loq _ e = viaShow e
instance Lang l => Loquate l E.Deadlock where loq _ e = viaShow e
instance Lang l => Loquate l E.FixIOException where loq _ e = viaShow e
instance Lang l => Loquate l E.IOException where loq _ e = viaShow e
instance Lang l => Loquate l E.NestedAtomically where loq _ e = viaShow e
instance Lang l => Loquate l E.NoMethodError where loq _ e = viaShow e
instance Lang l => Loquate l E.NonTermination where loq _ e = viaShow e
instance Lang l => Loquate l E.PatternMatchFail where loq _ e = viaShow e
instance Lang l => Loquate l E.RecConError where loq _ e = viaShow e
instance Lang l => Loquate l E.RecSelError where loq _ e = viaShow e
instance Lang l => Loquate l E.RecUpdError where loq _ e = viaShow e
instance Lang l => Loquate l E.SomeAsyncException where loq _ e = viaShow e
instance Lang l => Loquate l E.SomeException where loq _ e = viaShow e
instance Lang l => Loquate l E.TypeError where loq _ e = viaShow e

max_stars_repo_name: imuli/loquate

max_issues_repo_name: imuli/loquate

max_forks_repo_name: imuli/loquate


text:
module Approx where
import Data.Complex

approxDbl :: Int -> Double -> Double
approxDbl n x = fromInteger (round $ x * (10^n)) / (10.0^^n)

approx :: Int -> Complex Double -> Complex Double
approx n z = approxDbl n (realPart z) :+ approxDbl n (imagPart z)

max_stars_repo_name: stla/elliptic

max_issues_repo_name: stla/elliptic

max_forks_repo_name: stla/elliptic


text:
{-# LANGUAGE RecordWildCards #-}

module Regression.Logistic
  ( computeCost
  , gradientDescent
  , minimizeBFGS2
  ) where

import Numeric.GSL.Minimization (MinimizeMethodD(VectorBFGS2), minimizeVD)
import Numeric.LinearAlgebra (( #> ), (<.>), add)
import Numeric.LinearAlgebra.Data
  ( Matrix
  , R
  , Vector
  , (?)
  , atIndex
  , cmap
  , cols
  , flatten
  , fromList
  , scalar
  , size
  , toList
  , tr'
  )
import Common (MinimizationOpts(..), sigmoid, sigmoidMatrix)

computeCost :: Matrix R -> Vector R -> Vector R -> R
computeCost x y theta =
  let s = sigmoidMatrix (x #> theta)
      m = size y
   in ((y <.> log (s)) + ((1 - y) <.> log (1 - s))) / fromIntegral (-m)

gradientDescent :: Matrix R -> Vector R -> Vector R -> R -> Int -> R -> Vector R
gradientDescent _ _ _ alpha numIters regFactor
  | regFactor < 0 = error "Regularization factor cannot be < 0"
  | alpha < 0 = error "Learning factor cannot be < 0"
  | numIters < 0 = error "Number of iterations cannot be < 0"
gradientDescent x y theta alpha depth regFactor = go 0 theta
  where
    alpha' = scalar alpha
    --
    m = fromIntegral $ size y
    --
    go k accum
      | k >= depth = accum
      | otherwise =
        let delta = computeDelta accum
         in go (k + 1) (accum - (alpha' * delta))
    --
    penalizedTheta = fromList $ 0 : (tail $ ((* regFactor) <$> toList theta))
    --
    computeDelta th =
      let delta = (tr' x) #> ((sigmoidMatrix (x #> th)) - y)
       in (delta + penalizedTheta) / scalar m

minimizeBFGS2 ::
     Matrix R
  -> Vector R
  -> Vector R
  -> Int
  -> MinimizationOpts
  -> R
  -> Vector R
minimizeBFGS2 _ _ _ numIters _ regFactor
  | regFactor < 0 = error "Regularization factor cannot be < 0"
  | numIters < 0 = error "Number of iterations cannot be < 0"
minimizeBFGS2 x y theta numIters MinimizationOpts {..} regFactor =
  let (params, _) =
        minimizeVD
          VectorBFGS2
          precision
          numIters
          sizeOfFirstTrialStep
          tolerance
          (computeCost x y)
          computeGradients
          theta
   in params
  where
    m = fromIntegral $ size y
    --
    xTr = tr' x
    --
    penalizedTheta = fromList $ 0 : (tail $ ((* regFactor) <$> toList theta))
    --
    computeGradients :: Vector R -> Vector R
    computeGradients th =
      let delta = (tr' x) #> ((sigmoidMatrix (x #> th)) - y)
       in (delta + penalizedTheta) / scalar m

max_stars_repo_name: sebashack/mlToolBox

max_issues_repo_name: sebashack/mlToolBox

max_forks_repo_name: sebashack/mlToolBox


text:
-- |
-- Module    : Statistics.Resampling.Bootstrap
-- Copyright : (c) 2009, 2011 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- The bootstrap method for statistical inference.

module Statistics.Resampling.Bootstrap
    ( bootstrapBCA
    -- * References
    -- $references
    ) where

import           Data.Vector.Generic ((!))
import qualified Data.Vector.Unboxed as U

import Statistics.Distribution (cumulative, quantile)
import Statistics.Distribution.Normal
import Statistics.Resampling (Bootstrap(..), jackknife)
import Statistics.Sample (mean)
import Statistics.Types (Sample, CL, Estimate, ConfInt, estimateFromInterval,
                         estimateFromErr, CL, significanceLevel)
import qualified Statistics.Resampling as R

data T = {-# UNPACK #-} !Double :< {-# UNPACK #-} !Double
infixl 2 :<

-- | Bias-corrected accelerated (BCA) bootstrap. This adjusts for both
--   bias and skewness in the resampled distribution.
--
--   BCA algorithm is described in ch. 5 of Davison, Hinkley "Confidence
--   intervals" in section 5.3 "Percentile method"
bootstrapBCA
  :: CL Double       -- ^ Confidence level
  -> Sample          -- ^ Full data sample
  -> [(R.Estimator, Bootstrap U.Vector Double)]
  -- ^ Estimates obtained from resampled data and estimator used for
  --   this.
  -> [Estimate ConfInt Double]
bootstrapBCA confidenceLevel sample resampledData
  = map e resampledData
  where
    e (est, Bootstrap pt resample)
      | U.length sample == 1 || isInfinite bias =
          estimateFromErr      pt (0,0) confidenceLevel
      | otherwise =
          estimateFromInterval pt (resample ! lo, resample ! hi) confidenceLevel
      where
        -- Quantile estimates for given CL
        lo    = max (cumn a1) 0
          where a1 = bias + b1 / (1 - accel * b1)
                b1 = bias + z1
        hi    = min (cumn a2) (ni - 1)
          where a2 = bias + b2 / (1 - accel * b2)
                b2 = bias - z1
        -- Number of resamples
        ni    = U.length resample
        n     = fromIntegral ni
        -- Corrections
        z1    = quantile standard (significanceLevel confidenceLevel / 2)
        cumn  = round . (*n) . cumulative standard
        bias  = quantile standard (probN / n)
          where probN = fromIntegral . U.length . U.filter (<pt) $ resample
        accel = sumCubes / (6 * (sumSquares ** 1.5))
          where (sumSquares :< sumCubes) = U.foldl' f (0 :< 0) jack
                f (s :< c) j = s + d2 :< c + d2 * d
                    where d  = jackMean - j
                          d2 = d * d
                jackMean     = mean jack
        jack  = jackknife est sample

-- $references
--
-- * Davison, A.C; Hinkley, D.V. (1997) Bootstrap methods and their
--   application. <http://statwww.epfl.ch/davison/BMA/>

max_stars_repo_name: runeksvendsen/hs-gauge

max_issues_repo_name: runeksvendsen/hs-gauge

max_forks_repo_name: runeksvendsen/hs-gauge


text:
{-# LANGUAGE TemplateHaskell, OverloadedLists, DeriveGeneric, DeriveAnyClass, Strict #-}

module Data.PhysicsData where

import           Control.Lens      hiding (indices)
import           Number.Quaternion (T)
import           Numeric.LinearAlgebra (Matrix, Vector,(><), scale)

import           Utils.Quaternions
import Control.DeepSeq
import GHC.Generics
import Control.Arrow

import Physics.Collision

data PhysicsData = PhysicsData {
  _acc      :: Vector Float,
  _speed    :: Vector Float,
  _inertiaD :: Matrix Float,
  _massInv  :: Float,
  _inertiaInv :: Matrix Float,
  _angularS :: Vector Float,
  _torque   :: Vector Float } deriving Show

makeLenses ''PhysicsData

defaultInertiaDiag :: Float -> Matrix Float
defaultInertiaDiag m =
  (3><3)[ m, 0, 0
        , 0, m, 0
        , 0, 0, m ]

setRectInertia x y z m = set inertiaD . scale (m/12) $
  (3><3) [(x^2+z^2),         0,         0
         ,        0, (x^2+z^2),         0
         ,        0,         0, (x^2+y^2)]

carPhysics =
  setRectInertia 1 1 10 &&& startPhysics >>> uncurry (.)

startPhysics :: Float -> T Float -> PhysicsData
startPhysics mInv quat =
  PhysicsData { _acc        = [0,0,0]
              , _speed      = [0,0,0]
              , _inertiaD   = defaultInertiaDiag (1/mInv)
              , _massInv    = mInv
              , _inertiaInv = quatToLInv (defaultInertiaDiag (1/mInv)) quat
              , _angularS   = [0,0,0]
              , _torque     = [0,0,0] }

max_stars_repo_name: Antystenes/CPG

max_issues_repo_name: Antystenes/CPG

max_forks_repo_name: Antystenes/CPG


text:
-- (C) Copyright Chris Banks 2011-2012

-- This file is part of The Continuous Pi-calculus Workbench (BondCalculusWB).

--     BondCalculusWB is free software: you can redistribute it and/or modify
--     it under the terms of the GNU General Public License as published by
--     the Free Software Foundation, either version 3 of the License, or
--     (at your option) any later version.

--     BondCalculusWB is distributed in the hope that it will be useful,
--     but WITHOUT ANY WARRANTY; without even the implied warranty of
--     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
--     GNU General Public License for more details.

--     You should have received a copy of the GNU General Public License
--     along with BondCalculusWB.  If not, see <http://www.gnu.org/licenses/>.

module CPi.Plot
    (plotTimeSeries,
     plotTimeSeriesD,
     plotTimeSeriesFiltered,
     plotTimeSeriesFilteredD,
     plotTimeSeriesToFile,
     plotTimeSeriesToFileFiltered,
     phasePlot2,
     phasePlot2ToFile,
     phasePlot2D
    ) where

import Graphics.Rendering.Chart
import Graphics.Rendering.Chart.Gtk
import Graphics.Rendering.Chart.Backend.Cairo
import Data.Colour.Names
import Data.Colour.SRGB
import Data.Colour
-- import Data.Accessor
import Data.Default.Class
import Control.Lens
import qualified Control.Exception as X
import qualified Numeric.LinearAlgebra as LA
import Control.Monad
import Graphics.UI.Gtk.Misc.DrawingArea
import qualified Data.List as L

import BioCalcLib.Plot
import BioCalcLib.Lib
import CPi.Lib

-- Takes data from the ODE solver and plots them
plotTimeSeries :: LA.Vector Double -> LA.Matrix Double -> [Species] -> IO ()
plotTimeSeries ts soln ss
    = plot
      (LA.toList ts)
      (zip (map pretty ss) (map LA.toList (LA.toColumns soln)))

plotTimeSeriesD :: DrawingArea -> LA.Vector Double -> LA.Matrix Double -> [Species] -> IO ()
plotTimeSeriesD da ts soln ss
    = plotD
      da
      (LA.toList ts)
      (zip (map pretty ss) (map LA.toList (LA.toColumns soln)))

-- Plots the data to a PDF file
plotTimeSeriesToFile :: LA.Vector Double -> LA.Matrix Double -> [Species] -> String -> IO ()
plotTimeSeriesToFile ts soln ss
    = plotToFile
      (LA.toList ts)
      (zip (map pretty ss) (map LA.toList (LA.toColumns soln)))

-- Only plots selected species
plotTimeSeriesFiltered :: LA.Vector Double -> LA.Matrix Double -> [Species] -> [Species]
                       -> IO ()
plotTimeSeriesFiltered ts soln ss ss'
    = plot
      (LA.toList ts)
      (filter (\(s,_)-> s `elem` map specName ss')
       (zip (map specName ss) (map LA.toList (LA.toColumns soln))))

plotTimeSeriesFilteredD :: DrawingArea -> LA.Vector Double -> LA.Matrix Double -> [Species] -> [Species] -> IO ()
plotTimeSeriesFilteredD drawingArea ts soln ss ss'
    = plotD
      drawingArea
      (LA.toList ts)
      (filter (\(s,_)-> s `elem` map specName ss')
       (zip (map specName ss) (map LA.toList (LA.toColumns soln))))

-- Only plots selected species to a PDF file
plotTimeSeriesToFileFiltered :: LA.Vector Double -> LA.Matrix Double -> [Species] -> [Species]
                             -> String -> IO ()
plotTimeSeriesToFileFiltered ts soln ss ss'
    = plotToFile
      (LA.toList ts)
      (filter (\(s,_)-> s `elem` map specName ss')
       (zip (map pretty ss) (map LA.toList (LA.toColumns soln))))

---------------
-- Phase plots:
---------------

-- a plot of two dimensions:
phasePlot2 :: LA.Vector Double
           -> LA.Matrix Double
           -> [Species]
           -> (Species,Species)
           -> IO ()
phasePlot2 ts soln ss ss' = plotPhase $ filterPhasePlot ts soln ss ss'

filterPhasePlot :: LA.Vector Double
                -> LA.Matrix Double
                -> [Species]
                -> (Species,Species)
                -> [(String, [Double])]
filterPhasePlot _ soln ss ss'
    = filter (\(s,_) -> s == specName (fst ss') || s == specName (snd ss'))
             (zip (map specName ss) (map LA.toList $ LA.toColumns soln))

phasePlot2ToFile :: LA.Vector Double
                 -> LA.Matrix Double
                 -> [Species]
                 -> (Species, Species)
                 -> String
                 -> IO ()
phasePlot2ToFile ts soln ss ss'
    = plotPhaseToFile (filterPhasePlot ts soln ss ss')

phasePlot2D :: DrawingArea
            -> LA.Vector Double
            -> LA.Matrix Double
            -> [Species]
            -> (Species,Species)
            -> IO ()
phasePlot2D da ts soln ss ss' = plotPhaseD da $ filterPhasePlot ts soln ss ss'

max_stars_repo_name: twright/bondwb

max_issues_repo_name: twright/bondwb

max_forks_repo_name: twright/bondwb


text:
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DeriveAnyClass        #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE Strict                #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-|
Module      : Grenade.Layers.LeakyTanh
Description : Hyperbolic tangent nonlinear layer
Copyright   : (c) Huw Campbell, 2016-2017
License     : BSD2
Stability   : experimental
-}
module Grenade.Layers.LeakyTanh
  ( LeakyTanh(..)
  , SpecLeakyTanh (..)
  , specLeakyTanh1D
  , specLeakyTanh2D
  , specLeakyTanh3D
  , specLeakyTanh
  , leakyTanhLayer
  ) where

import           Control.DeepSeq                (NFData (..))
import           Data.Constraint                (Dict (..))
import           Data.Proxy
import           Data.Reflection                (reifyNat)
import           Data.Serialize
import           Data.Singletons
import           GHC.Generics                   (Generic)
import           GHC.TypeLits
import qualified Numeric.LinearAlgebra.Static   as LAS
import           Unsafe.Coerce                  (unsafeCoerce)

import           Grenade.Core
import           Grenade.Dynamic
import           Grenade.Dynamic.Internal.Build
import           Grenade.Types
import           Grenade.Utils.Conversion
import           Grenade.Utils.Vector

-- | A LeakyTanh layer. A layer which can act between any shape of the same dimension, performing a tanh function. The maximum value is given in Promille, i.e. 995 is 0.995, and is positive and negative. I.e. LeakyTanh(x) = min 0.995 (max (-0.995) x).
data LeakyTanh (v :: Nat) = LeakyTanh
  deriving (Generic,NFData,Show)

instance UpdateLayer (LeakyTanh maxVal) where
  type Gradient (LeakyTanh maxVal) = ()
  runUpdate _ l@LeakyTanh{} _ = l

instance RandomLayer (LeakyTanh maxVal) where
  createRandomWith _ _ = return LeakyTanh

instance Serialize (LeakyTanh maxVal) where
  put _ = return ()
  get = return LeakyTanh

instance (a ~ b, SingI a, KnownNat maxVal) => Layer (LeakyTanh maxVal) a b where
  type Tape (LeakyTanh maxVal) a b = S a
  runForwards _ (S1DV v) = (S1DV v, S1DV $ mapVector (tanhMax maxVal) v)
    where
      maxVal = (/ 1000) $ fromIntegral $ max 0 $ min 1000 $ natVal (Proxy :: Proxy maxVal)
  runForwards _ (S2DV v) = (S2DV v, S2DV $ mapVector (tanhMax maxVal) v)
    where
      maxVal = (/ 1000) $ fromIntegral $ max 0 $ min 1000 $ natVal (Proxy :: Proxy maxVal)
  runForwards _ (S1D a) = (S1D a, S1D $ LAS.dvmap (tanhMax maxVal) a)
    where
      maxVal = (/ 1000) $ fromIntegral $ max 0 $ min 1000 $ natVal (Proxy :: Proxy maxVal)
  runForwards _ (S2D a) = (S2D a, S2D $ LAS.dmmap (tanhMax maxVal) a)
    where
      maxVal = (/ 1000) $ fromIntegral $ max 0 $ min 1000 $ natVal (Proxy :: Proxy maxVal)
  runForwards _ (S3D a) = (S3D a, S3D $ LAS.dmmap (tanhMax maxVal) a)
    where
      maxVal = (/ 1000) $ fromIntegral $ max 0 $ min 1000 $ natVal (Proxy :: Proxy maxVal)
  runBackwards _ (S1DV v) (S1DV gs) = ((), S1DV $ zipWithVector (\t g -> tanhMax' maxVal t * g) v gs)
    where
      maxVal = (/ 1000) $ fromIntegral $ max 0 $ min 1000 $ natVal (Proxy :: Proxy maxVal)
  runBackwards _ (S2DV v) (S2DV gs) = ((), S2DV $ zipWithVector (\t g -> tanhMax' maxVal t * g) v gs)
    where
      maxVal = (/ 1000) $ fromIntegral $ max 0 $ min 1000 $ natVal (Proxy :: Proxy maxVal)
  runBackwards _ (S1D a) (S1D g) = ((), S1D $ LAS.dvmap (tanhMax' maxVal) a * g)
    where
      maxVal = (/ 1000) $ fromIntegral $ max 0 $ min 1000 $ natVal (Proxy :: Proxy maxVal)
  runBackwards _ (S2D a) (S2D g) = ((), S2D $ LAS.dmmap (tanhMax' maxVal) (tanh' a) * g)
    where
      maxVal = (/ 1000) $ fromIntegral $ max 0 $ min 1000 $ natVal (Proxy :: Proxy maxVal)
  runBackwards _ (S3D a) (S3D g) = ((), S3D $ LAS.dmmap (tanhMax' maxVal) (tanh' a) * g)
    where
      maxVal = (/ 1000) $ fromIntegral $ max 0 $ min 1000 $ natVal (Proxy :: Proxy maxVal)
  runBackwards l x y = runBackwards l x (toLayerShape x y)

tanhMax :: (Ord a, Fractional a) => a -> a -> a
tanhMax m v
  | v > m = leaky * (v - m) + m
  | v < -m = leaky * (v - m) - m
  | otherwise = v
  where
    leaky = 0.05

tanhMax' :: (Ord a, Fractional a, Floating a) => a -> a -> a
tanhMax' m v
  | v > m = leaky
  | v < -m = -leaky
  | otherwise = max 0.005 (tanh' v)
  where
    leaky = 0.05

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

-------------------- DynamicNetwork instance --------------------

instance (KnownNat maxVal) => FromDynamicLayer (LeakyTanh maxVal) where
  fromDynamicLayer inp _ LeakyTanh = SpecNetLayer $ SpecLeakyTanh maxVal (tripleFromSomeShape inp)
    where maxVal = (/1000) $ fromIntegral $ max 0 $ min 1000 $ natVal (Proxy :: Proxy maxVal)

instance ToDynamicLayer SpecLeakyTanh where
  toDynamicLayer _ _ (SpecLeakyTanh maxVal (rows, cols, depth)) =
     reifyNat rows $ \(_ :: (KnownNat rows) => Proxy rows) ->
     reifyNat cols $ \(_ :: (KnownNat cols) => Proxy cols) ->
     reifyNat depth $ \(_ :: (KnownNat depth) => Proxy depth) ->
     reifyNat (round $ 1000 * maxVal) $ \(_ :: (KnownNat maxVal) => Proxy maxVal) ->
     case (rows, cols, depth) of
         (_, 1, 1) -> case (unsafeCoerce (Dict :: Dict()) :: Dict ()) of
           Dict -> return $ SpecLayer (LeakyTanh :: LeakyTanh maxVal) (sing :: Sing ('D1 rows)) (sing :: Sing ('D1 rows))
         (_, _, 1) -> case (unsafeCoerce (Dict :: Dict()) :: Dict ()) of
           Dict -> return $ SpecLayer (LeakyTanh :: LeakyTanh maxVal) (sing :: Sing ('D2 rows cols)) (sing :: Sing ('D2 rows cols))
         _         -> case (unsafeCoerce (Dict :: Dict()) :: Dict (KnownNat (rows GHC.TypeLits.* depth))) of
           Dict -> return $ SpecLayer (LeakyTanh :: LeakyTanh maxVal) (sing :: Sing ('D3 rows cols depth)) (sing :: Sing ('D3 rows cols depth))

-- | Create a specification for a LeakyTanh layer.
specLeakyTanh1D :: RealNum -> Integer -> SpecNet
specLeakyTanh1D maxVal i = specLeakyTanh3D maxVal (i, 1, 1)

-- | Create a specification for a LeakyTanh layer.
specLeakyTanh2D :: RealNum -> (Integer, Integer) -> SpecNet
specLeakyTanh2D maxVal (i, j) = specLeakyTanh3D maxVal (i, j, 1)

-- | Create a specification for a LeakyTanh layer.
specLeakyTanh3D :: RealNum -> (Integer, Integer, Integer) -> SpecNet
specLeakyTanh3D maxVal = SpecNetLayer . SpecLeakyTanh maxVal

-- | Create a specification for a LeakyTanh layer.
specLeakyTanh :: RealNum -> (Integer, Integer, Integer) -> SpecNet
specLeakyTanh maxVal = SpecNetLayer . SpecLeakyTanh maxVal

-- | Add a LeakyTanh layer to your build.
leakyTanhLayer :: RealNum -> BuildM ()
leakyTanhLayer maxVal
  | maxVal <= 0 || maxVal > 1 = error "The maxVal parameter of tanhMaxLayer has to be in (0,1]"
  | otherwise = buildGetLastLayerOut >>= buildAddSpec . SpecNetLayer . SpecLeakyTanh maxVal

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

instance GNum (LeakyTanh maxVal) where
  _ |* _ = LeakyTanh
  _ |+ _ = LeakyTanh
  sumG _ = LeakyTanh

max_stars_repo_name: schnecki/grenade

max_issues_repo_name: schnecki/grenade

max_forks_repo_name: schnecki/grenade


text:
{-# LANGUAGE ForeignFunctionInterface #-}
{-|
Module      : Grenade.Layers.Internal.Add
Description : Fast addition functions that call efficient C function
Maintainer  : Theo Charalambous
License     : BSD2
Stability   : experimental
-}
module Grenade.Layers.Internal.Add (
  addPerChannel
) where

import qualified Data.Vector.Storable        as U (unsafeFromForeignPtr0,
                                                   unsafeToForeignPtr0)

import           Foreign                     (mallocForeignPtrArray, withForeignPtr)
import           Foreign.Ptr                 (Ptr)
import           Numeric.LinearAlgebra       (Matrix, Vector, flatten)
import qualified Numeric.LinearAlgebra.Devel as U
import           System.IO.Unsafe            (unsafePerformIO)

import           Grenade.Types

-- | Add the nth element to a vector to every pixel in the nth channel of a matrix.
--   It assumes that the size of the vector and the number of channels in the matrix
--   are both equal to channels.
addPerChannel :: Int -> Int -> Int -> Matrix RealNum -> Vector RealNum -> Matrix RealNum
addPerChannel channels rows cols m b
  = let outMatSize      = rows * cols * channels
        vec             = flatten m
    in unsafePerformIO $ do
      outPtr        <- mallocForeignPtrArray outMatSize
      let (inPtr, _) = U.unsafeToForeignPtr0 vec
          (bPtr, _)  = U.unsafeToForeignPtr0 b

      withForeignPtr inPtr $ \inPtr' ->
        withForeignPtr bPtr $ \bPtr' ->
          withForeignPtr outPtr $ \outPtr' ->
            add_per_channel_cpu inPtr' channels rows cols bPtr' outPtr'

      let matVec = U.unsafeFromForeignPtr0 outPtr outMatSize
      return (U.matrixFromVector U.RowMajor (rows * channels) cols matVec)
{-# INLINE addPerChannel #-}

foreign import ccall unsafe
    add_per_channel_cpu
      :: Ptr RealNum -> Int -> Int -> Int -> Ptr RealNum -> Ptr RealNum -> IO ()

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-# LANGUAGE Arrows, ExistentialQuantification, ScopedTypeVariables #-}

module SignalLevel.SpecAnalysis20 where

import Data.Complex
import Data.Maybe (catMaybes,listToMaybe)
import qualified Data.Map as M
import Control.Arrow.ArrowP
import Control.DeepSeq (NFData)
import Control.SF.SF (SF(..))
import FRP.UISF (SEvent,asyncVT,Time,DeltaT)
import FRP.UISF.Asynchrony (Automaton(..))
import FRP.UISF.Graphics.Color (rgbE)

import Euterpea hiding (SF)
import HSoM
import HSoM.MUI (asyncV)
import HSoM.Examples.FFT (fftA,quantize,presentFFT)

import SignalLevel.SigFun19

dft :: RealFloat a => [Complex a] -> [Complex a]
dft xs =
  let lenI = length xs
      lenR = fromIntegral lenI
      lenC = lenR :+ 0
  in [let i = -2*pi * fromIntegral k / lenR
      in (1/lenC) * sum [(xs !! n) * exp (0 :+ i * fromIntegral n)
                          | n <- [0,1..lenI-1]]
       | k <- [0,1..lenI-1]]

-- num = # samples
-- n   = term

mkTerm :: Int -> Double -> [Complex Double]
mkTerm num n = let f = 2*pi / fromIntegral num
  in [sin (n*f*fromIntegral i) / n :+ 0
       | i <- [0,1..num-1]]

mkxa,mkxb,mkxc :: Int -> [Complex Double]
mkxa num = mkTerm num 1
mkxb num = zipWith (+) (mkxa num) (mkTerm num 3)
mkxc num = zipWith (+) (mkxb num) (mkTerm num 5)

mkPulse :: Int -> [Complex Double]
mkPulse n = 100 : take (n-1) (repeat 0)

x1 num = let f = 2*pi*pi / fromIntegral num
  in map (:+ 0) [sin (f * fromIntegral i)
                  | i <- [0,1..num-1]]

mkPolars :: [Complex Double] -> [Complex Double]
mkPolars = map ((\(m,p) -> m :+ p) . polar)

---------------
-- Exercises --
---------------

{- Exercise 20.1©:Write a Haskell function idft that implements the inverse DFT
 - as captured in Equation 20.3. Test your code by applying idft to one of the
 - signals used earlier in this section. In other words, show empirically that,
 - up to round-off errors, idft (dft xs) == xs. -}

idft :: RealFloat a => [Complex a] -> [Complex a]
idft xs = let lenI = length xs
              lenR = fromIntegral lenI
          in [let i = 2*pi* fromIntegral n / lenR
              in sum [(xs !! k) * exp (0 :+ i * fromIntegral k)
                     | k <- [0..lenI-1]]
             | n <- [0..lenI-1]]

{- Exercise 20.2: Use dft to analyze some of the signals generated using signal
 - functions defined in Chapter 19. -}

fft1 :: forall p . Clock p => SigFun p () (SEvent (M.Map Double Double))
fft1 = s1 >>> (fftA 128 256) >>> arr (fmap $ presentFFT $ rate (undefined :: p))

{- Exercise 20.3:© Define a function mkSqWave :: Int -> Int -> [Complex Double]
 - such that mkSqWave num n is the sum of the first n terms of the Fourier
 - Series of a square wave, having num samples in the result. -}

mkSqWave :: Int -> Int -> [Complex Double]
mkSqWave num n =
  foldr (zipWith (+) . (\i -> mkTerm i))
        (replicate num 0)
        (filter odd [1..(2*n-1)])
  where mkTerm x = let f = 2 * pi / fromIntegral num
                       t = fromIntegral x
                   in [sin (t*f*fromIntegral i)/t:+0
                      | i <- [0,1..num-1]]

{- Exercise 20.4: -}

{- a. © Prove mathematically that x and x^ are inverses.

  Let X^ and X be the NxN matrices defined by

  X^(k,n) = (1/N) * w^(kn)
  X(n,k)    = (1/N) * w^(-(kn))

  where w = e^(-j2π/N)
  and (a,b) corresponds to the entry in the ath row and bth column of the matrix

  Prove X^ and X are inverses
  Proof:

    It suffices to prove (X^)X = I

    Thus I show
      if p=q then X^X(p,q) = 1
      otherwise   X^X(p,q) = 0

    X^X(p,q) = X^(p,*)·X(*,q) = (1/N) ∑w^(pn)*w^(-(n*q))

    where X^(p,*) is the pth row of X^ , X(*,q) is the qth column of X and
    ∑ is a sum from 0 to N-1

    if p = q then
      X^X(p,q) = (1/N) ∑ w^(pn)/w^(nq) = (1/N) ∑ 1 =  (1/N)*N = 1

    if p ≠ q then
      X^X(p,q) = (1/N) ∑ w^(n(p-q))
               = (1 - (e^(-j*2π(p-q)/N))^N)/(1 - e^(-j*2π(p-q)/N))
               = (1 - (e^(-j*2π(p-q)))/(1 - e^(-j*2π(p-q)/N))
               = (1 - 1/(cos(2π(p-q)) + jsin(2π(p-q))))/(1 - e^(-j*2π(p-q)/N))
               = (1 - 1)/(1 - e^(-j*2π(p-q)/N))
               = 0

    Thus X^X = I, proving X^ and X are inverses.
-}

{-
e1 =
  \i -> (1/lenC)*sum[(cxs !! n)*exp(0 :+ i*fromIntegral n) | n <- [0,1..lenI-1]]

f = \exp -> let lenI = length cxs
            in let lenR = fromIntegral lenI
               in let lenC = lenR :+ 0
                  in exp

-> fdft = (f (concatMap
               (\k -> let i = -2*pi*fromIntegral k/lenR in e1 i)
               [0,1..lenI-1]
             )
          )

fix exp = \x -> let x = exp x in x
-}

{- b. Also prove, using equational reasoning, that dft and idft are
 - inverses. (..., you may assume that Haskell numeric types
 - obey the standard axioms of real arithmetic.)

idft(dft(cxs))  = cxs

idft(dft(cxs))
-> idft (let lenI = length cxs
             lenR = fromIntegral lenI
             lenC = lenR :+ 0
         in [let i = -2*pi*fromIntegral k/lenR
             in (1/lenC)*sum[(cxs !! n)*exp (0 :+ i*fromIntegral n)
                             | n <-  [0,1..lenI-1]]
              |k <- [0,1..lenI-1]] )

-> idft (let lenI = length cxs
             lenR = fromIntegral lenI
             lenC = lenR :+ 0
         in [let i = -2*pi*fromIntegral k/lenR
             in e1 i | k <- [0,1..lenI-1]] )

  -> concatMap
       (\k -> let i = -2*pi*fromIntegral k/lenR in e1 i)
       [0,1..lenI-1]

-> idft (f (concatMap
             (\k -> let i = -2*pi*fromIntegral k/lenR in e1 i)
             [0,1..lenI-1]
           )
        )

-> let lenI2 = length fdft
       lenR2 = fromIntegral lenI2
   in [let i2 = 2*pi* fromIntegral n2 / lenR2
       in sum [(fdft !! k) * exp (0 :+ i2 * fromIntegral k) | k <- [0..lenI2-1]]
        | n2 <- [0..lenI2-1]]

-> let lenI2 = length fdft
       lenR2 = fromIntegral lenI2
   in [let i2 = 2*pi* fromIntegral n2 / lenR2
       in sum [(fdft !! k) * exp (0 :+ i2 * fromIntegral k) | k <- [0..lenI2-1]]
        | n2 <- [0..lenI2-1]]
-}

fftEx :: UISF () ()
fftEx =
  proc _ -> do
    f <- withDisplay $ hSlider (1,2000) 440 -< ()
    d <- clockedSFToUISF (0.1) fftDataSig -< f
    let (ds,ts) = unzip d
    _ <- histogram (makeLayout (Fixed {fixedSize = 500})
                               (Fixed {fixedSize = 150}))
                                -< listToMaybe (catMaybes (snd $ unzip ds))
    _ <- realtimeGraph (makeLayout (Fixed {fixedSize = 500})
                                   (Fixed {fixedSize = 150}))
                                   1 Black -< (zip (fst $ unzip ds) ts)
    m  <- withDisplay $ hSlider (1,2000) 20 -< ()
    d1 <- clockedSFToUISF (0.1) fftDataSig -< m
    let (ds1,ts1) = unzip d1
    _  <- histogram (makeLayout (Fixed {fixedSize = 500})
                               (Fixed {fixedSize = 150}))
                                -< listToMaybe (catMaybes (snd $ unzip ds1))
    _ <- realtimeGraph (makeLayout (Fixed {fixedSize = 500})
                                   (Fixed {fixedSize = 150}))
                                   1 Black -< (zip (fst $ unzip ds1) ts1)
    d3 <- clockedSFToUISF (0.1) fftDataSig -< (f+m)
    let (ds3,ts3) = unzip d3
    _  <- histogram (makeLayout (Fixed {fixedSize = 500})
                               (Fixed {fixedSize = 150}))
                                -< listToMaybe (catMaybes (snd $ unzip ds3))
    outA -< ()
    where
      fftDataSig :: SigFun CtrRate Double (Double, SEvent [Double])
      fftDataSig = proc f -> do
        s <- osc (tableSinesN 4096 [1]) 0 -< f
        fftData <- fftA 128 256 -< s
        outA -< (s,fftData)

t0 = runMUI (defaultMUIParams
               {uiSize       = (500,900),
                uiBackground = rgbE 93 10 13,
                uiTitle      = "FFT Test"})
            fftEx

clockedSFToUISF :: forall a b c . (NFData b,Clock c) =>
                   DeltaT -> SigFun c a b -> UISF a [(b,Time)]
clockedSFToUISF buffer ~(ArrowP sf) = let r = rate (undefined :: c)
  in asyncV r buffer (toAutomaton sf)

toAutomaton :: forall a b . SF a b -> Automaton (->) a b
toAutomaton ~(SF f) = Automaton $ \a -> let (b,sf) = f a in (b,toAutomaton sf)

{- Exercise 20.5: Modify the program in 20.5 (fftEx) in the following ways:
 -
 - 1. Add a second slider, and use it to control the frequency of a second
 - oscillator.©
 - 2. Let s1 and s2 be the names of the signals whose frequencies are controlled
 - by the first and second sliders, respectively. Instead of displaying the FFT
 - of just s1, try a variety of combinations of s1 and s2, such as
 - s1 + s2
 - s1 - s2
 - s1 * s2
 - 1/s1 + 1/s2
 - s1/s2
 - Comment on the results.
 - 3. Use s2 to control the frequency of s1 (as was done with vibrato in Ch.19).
 - Plot the fft of s1 and comment on the result.
 - 4. Instead of using osc to generate a pure sine wave, try using other
 - oscillators and/or table generators to create more complex tones, and plot
 - their FFTs. Comment on the results. -}

----------
-- Util --
----------
printComplexL :: [Complex Double] -> IO ()
printComplexL xs =
  let f (i, rl :+ im) =
        do putStr (spaces (3 - length (show i)))
           putStr (show i ++ ":    ("          )
           putStr (niceNum rl ++ ", "         )
           putStr (niceNum im ++ ") \n"        )
   in mapM_ f (zip [0..length xs - 1] xs)

niceNum :: Double -> String
niceNum d =
  let d' = fromIntegral (round (1e10*d))/1e10
      (dec,fra) = break (== '.') (show d')
      (fra',exp) = break (== 'e') fra
  in spaces (3 - length dec) ++ dec ++ take 11 fra'
     ++ exp ++ spaces (12 - length fra' - length exp)

spaces :: Int -> String
spaces n = take n (repeat ' ')

max_stars_repo_name: REBradley/music-lab

max_issues_repo_name: REBradley/music-lab

max_forks_repo_name: REBradley/music-lab