text:
{-# LANGUAGE DefaultSignatures          #-}
{-# LANGUAGE DeriveGeneric              #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE StandaloneDeriving         #-}

{-|
Module: Data.Semiring
Description: Haskell semirings
License: MIT
Maintainer: mail@doisinkidney.com
Stability: experimental
-}

module Data.Semiring
  ( Semiring(..)
  , Add(..)
  , Mul(..)
  , mul
  , add
  ) where

import           Data.Coerce           (coerce)

import           Data.Functor.Identity (Identity (..))

import           Data.Monoid
import           Data.Semigroup        (Max (..), Min (..))

import           Data.Complex          (Complex)
import           Data.Fixed            (Fixed, HasResolution)
import           Data.Ratio            (Ratio)
import           Numeric.Natural       (Natural)

import           Data.Set              (Set)
import qualified Data.Set              as Set

import           Data.Int              (Int16, Int32, Int64, Int8)
import           Data.Word             (Word16, Word32, Word64, Word8)
import           Foreign.C.Types       (CChar, CClock, CDouble, CFloat, CInt,
                                        CIntMax, CIntPtr, CLLong, CLong,
                                        CPtrdiff, CSChar, CSUSeconds, CShort,
                                        CSigAtomic, CSize, CTime, CUChar, CUInt,
                                        CUIntMax, CUIntPtr, CULLong, CULong,
                                        CUSeconds, CUShort, CWchar)

import           Foreign.Ptr           (IntPtr, WordPtr)
import           System.Posix.Types    (CCc, CDev, CGid, CIno, CMode, CNlink,
                                        COff, CPid, CRLim, CSpeed, CSsize,
                                        CTcflag, CUid, Fd)

import           GHC.Generics          (Generic, Generic1)

import           Test.QuickCheck       (Arbitrary)

-- | A <https://en.wikipedia.org/wiki/Semiring Semiring> is like the
-- the combination of two 'Data.Monoid.Monoid's. The first
-- is called '<+>'; it has the identity element 'zero', and it is
-- commutative. The second is called '<.>'; it has identity element 'one',
-- and it must distribute over '<+>'.
--
-- = Laws
-- == Normal 'Precursor.Algebra.Monoid.Monoid' laws
-- * @(a '<+>' b) '<+>' c = a '<+>' (b '<+>' c)@
-- * @'zero' '<+>' a = a '<+>' 'zero' = a@
-- * @(a '<.>' b) '<.>' c = a '<.>' (b '<.>' c)@
-- * @'one' '<.>' a = a '<.>' 'one' = a@
--
-- == Commutativity of '<+>'
-- * @a '<+>' b = b '<+>' a@
--
-- == Distribution of '<.>' over '<+>'
-- * @a '<.>' (b '<+>' c) = (a '<.>' b) '<+>' (a '<.>' c)@
-- * @(a '<+>' b) '<.>' c = (a '<.>' c) '<+>' (b '<.>' c)@
--
-- Another useful law, annihilation, may be deduced from the axioms
-- above:
--
-- * @'zero' '<.>' a = a '<.>' 'zero' = 'zero'@
class Semiring a where
  -- | The identity of '<+>'.
  zero :: a
  -- | The identity of '<.>'.
  one :: a
  -- | An associative binary operation, which distributes over '<+>'.
  infixl 7 <.>
  (<.>) :: a -> a -> a
  -- | An associative, commutative binary operation.
  infixl 6 <+>
  (<+>) :: a -> a -> a

  default zero :: Num a => a
  default one :: Num a => a
  default (<+>) :: Num a => a -> a -> a
  default (<.>) :: Num a => a -> a -> a

  zero = 0
  one = 1
  (<+>) = (+)
  (<.>) = (*)

------------------------------------------------------------------------
-- Instances
------------------------------------------------------------------------
instance Semiring Bool where
  one = True
  zero = False
  (<+>) = (||)
  (<.>) = (&&)

instance Semiring () where
  one = ()
  zero = ()
  _ <+> _ = ()
  _ <.> _ = ()

cartProd :: (Ord a, Monoid a) => Set a -> Set a -> Set a
cartProd xs ys =
  Set.foldl' (\a x ->
                Set.foldl' (flip (Set.insert . mappend x)) a ys)
  Set.empty xs

-- | The 'Set' 'Semiring' is 'Data.Set.union' for '<+>', and a Cartesian
-- product for '<.>'.
instance (Ord a, Monoid a) => Semiring (Set a) where
  (<.>) = cartProd
  (<+>) = Set.union
  zero = Set.empty
  one = Set.singleton mempty

------------------------------------------------------------------------
-- Addition and multiplication newtypes
------------------------------------------------------------------------

type WrapBinary f a = (a -> a -> a) -> f a -> f a -> f a

-- | Monoid under '<+>'. Analogous to 'Data.Monoid.Sum', but uses the
-- 'Semiring' constraint, rather than 'Num'.
newtype Add a = Add
  { getAdd :: a
  } deriving (Eq, Ord, Read, Show, Enum, Bounded, Generic, Generic1, Num
             ,Arbitrary)

-- | Monoid under '<.>'. Analogous to 'Data.Monoid.Product', but uses the
-- 'Semiring' constraint, rather than 'Num'.
newtype Mul a = Mul
  { getMul :: a
  } deriving (Eq, Ord, Read, Show, Enum, Bounded, Generic, Generic1, Num
             ,Arbitrary)

add :: (Semiring a, Foldable f) => f a -> a
add = getAdd . foldMap Add

mul :: (Semiring a, Foldable f) => f a -> a
mul = getMul . foldMap Mul

instance Functor Add where fmap = coerce

instance Functor Mul where fmap = coerce

instance Foldable Add where
  foldr   =
    (coerce :: ((a -> b -> c) -> (b -> a -> c))
            -> (a -> b -> c)
            -> (b -> Add a -> c)) flip
  foldl   = coerce
  foldMap = coerce
  length  = const 1

instance Foldable Mul where
  foldr   =
    (coerce :: ((a -> b -> c) -> (b -> a -> c))
            -> (a -> b -> c)
            -> (b -> Mul a -> c)) flip
  foldl   = coerce
  foldMap = coerce
  length  = const 1

instance Applicative Add where
  pure = coerce
  (<*>) =
    (coerce :: ((a -> b) -> a -> b)
            -> (Add (a -> b) -> Add a -> Add b)) ($)

instance Applicative Mul where
  pure = coerce
  (<*>) =
    (coerce :: ((a -> b) -> a -> b)
            -> (Mul (a -> b) -> Mul a -> Mul b)) ($)

instance Monad Add where
  (>>=) = flip coerce

instance Monad Mul where
  (>>=) = flip coerce

instance Semiring a => Monoid (Add a) where
  mempty = Add zero
  mappend = (coerce :: WrapBinary Add a) (<+>)

instance Semiring a => Monoid (Mul a) where
  mempty = Mul one
  mappend = (coerce :: WrapBinary Mul a) (<.>)

instance Semiring a => Semiring (Add a) where
  zero = Add zero
  one = Add one
  (<+>) = (coerce :: WrapBinary Add a) (<+>)
  (<.>) = (coerce :: WrapBinary Add a) (<.>)

instance Semiring a => Semiring (Mul a) where
  zero = Mul zero
  one = Mul one
  (<+>) = (coerce :: WrapBinary Mul a) (<+>)
  (<.>) = (coerce :: WrapBinary Mul a) (<.>)

------------------------------------------------------------------------
-- Ord wrappers
------------------------------------------------------------------------

-- | The 'Semiring' for 'Max' uses the 'max' operation for '<+>', and
-- normal '+' for '<.>'.
instance (Ord a, Bounded a, Semiring a) => Semiring (Max a) where
  (<+>) = mappend
  zero = mempty
  (<.>) = (coerce :: WrapBinary Max a) (<+>)
  one = Max zero

-- | The 'Semiring' for 'Min' uses the 'min' operation for '<+>', and
-- normal '+' for '<.>'.
instance (Ord a, Bounded a, Semiring a) => Semiring (Min a) where
  (<+>) = mappend
  zero = mempty
  (<.>) = (coerce :: WrapBinary Min a) (<+>)
  one = Min zero

------------------------------------------------------------------------
-- (->) instance
------------------------------------------------------------------------

-- | The @(->)@ instance is analogous to the one for 'Monoid'.
instance Semiring b => Semiring (a -> b) where
  zero = const zero
  one = const one
  (f <+> g) x = f x <+> g x
  (f <.> g) x = f x <.> g x

------------------------------------------------------------------------
-- Endo instance
------------------------------------------------------------------------

-- | The 'Endo' semiring uses function composition for '<.>', and
-- pointwise 'mappend' for '<+>'. The underlying 'Monoid' needs to be
-- commutative.
instance Monoid a => Semiring (Endo a) where
  (<.>) = mappend
  one = mempty
  Endo f <+> Endo g = Endo (\x -> f x `mappend` g x)
  zero = Endo (const mempty)

------------------------------------------------------------------------
-- Boring instances
------------------------------------------------------------------------
instance Semiring Float
instance Semiring Double
instance Semiring Int
instance Semiring Int8
instance Semiring Int16
instance Semiring Int32
instance Semiring Int64
instance Semiring Integer
instance Semiring Word
instance Semiring Word8
instance Semiring Word16
instance Semiring Word32
instance Semiring Word64
instance Semiring CUIntMax
instance Semiring CIntMax
instance Semiring CUIntPtr
instance Semiring CIntPtr
instance Semiring CSUSeconds
instance Semiring CUSeconds
instance Semiring CTime
instance Semiring CClock
instance Semiring CSigAtomic
instance Semiring CWchar
instance Semiring CSize
instance Semiring CPtrdiff
instance Semiring CDouble
instance Semiring CFloat
instance Semiring CULLong
instance Semiring CLLong
instance Semiring CULong
instance Semiring CLong
instance Semiring CUInt
instance Semiring CInt
instance Semiring CUShort
instance Semiring CShort
instance Semiring CUChar
instance Semiring CSChar
instance Semiring CChar
instance Semiring IntPtr
instance Semiring WordPtr
instance Semiring Fd
instance Semiring CRLim
instance Semiring CTcflag
instance Semiring CSpeed
instance Semiring CCc
instance Semiring CUid
instance Semiring CNlink
instance Semiring CGid
instance Semiring CSsize
instance Semiring CPid
instance Semiring COff
instance Semiring CMode
instance Semiring CIno
instance Semiring CDev
instance Semiring Natural
instance Integral a => Semiring (Ratio a)
deriving instance Semiring a => Semiring (Product a)
deriving instance Semiring a => Semiring (Sum a)
instance RealFloat a => Semiring (Complex a)
instance HasResolution a => Semiring (Fixed a)
deriving instance Semiring a => Semiring (Identity a)

------------------------------------------------------------------------
-- Very boring instances
------------------------------------------------------------------------
instance (Semiring a, Semiring b) => Semiring (a,b) where
        zero = (zero, zero)
        (a1,b1) <+> (a2,b2) =
                (a1 <+> a2, b1 <+> b2)
        one = (one, one)
        (a1,b1) <.> (a2,b2) =
                (a1 <.> a2, b1 <.> b2)

instance (Semiring a, Semiring b, Semiring c) => Semiring (a,b,c) where
        zero = (zero, zero, zero)
        (a1,b1,c1) <+> (a2,b2,c2) =
                (a1 <+> a2, b1 <+> b2, c1 <+> c2)
        one = (one, one, one)
        (a1,b1,c1) <.> (a2,b2,c2) =
                (a1 <.> a2, b1 <.> b2, c1 <.> c2)

instance (Semiring a, Semiring b, Semiring c, Semiring d) => Semiring (a,b,c,d) where
        zero = (zero, zero, zero, zero)
        (a1,b1,c1,d1) <+> (a2,b2,c2,d2) =
                (a1 <+> a2, b1 <+> b2,
                 c1 <+> c2, d1 <+> d2)
        one = (one, one, one, one)
        (a1,b1,c1,d1) <.> (a2,b2,c2,d2) =
                (a1 <.> a2, b1 <.> b2, c1 <.> c2, d1 <.> d2)

instance (Semiring a, Semiring b, Semiring c, Semiring d, Semiring e) =>
                Semiring (a,b,c,d,e) where
        zero = (zero, zero, zero, zero, zero)
        (a1,b1,c1,d1,e1) <+> (a2,b2,c2,d2,e2) =
                (a1 <+> a2, b1 <+> b2, c1 <+> c2,
                 d1 <+> d2, e1 <+> e2)
        one = (one, one, one, one, one)
        (a1,b1,c1,d1,e1) <.> (a2,b2,c2,d2,e2) =
                (a1 <.> a2, b1 <.> b2, c1 <.> c2,
                 d1 <.> d2, e1 <.> e2)

max_stars_repo_name: oisdk/Playgrounds

max_issues_repo_name: oisdk/Playgrounds

max_forks_repo_name: oisdk/Playgrounds


text:
{-# OPTIONS_GHC -Wall #-}

{-|

  Module      : AutoBench.Internal.Types
  Description : Datatypes and associated helper functions\/defaults.
  Copyright   : (c) 2018 Martin Handley
  License     : BSD-style
  Maintainer  : martin.handley@nottingham.ac.uk
  Stability   : Experimental
  Portability : GHC

  Datatypes used internally throughout AutoBench's implementation and any
  associated helper functions\/defaults.

-}

{-
   ----------------------------------------------------------------------------
   <TO-DO>:
   ----------------------------------------------------------------------------
   - 'DataOpts' Discover setting;
   - Make AnalOpts in TestSuite a maybe type? In case users don't want to
     analyse right away;
   - 'UserInputs' PP doesn't wrap;
   - 'TestSuite's PP isn't alphabetical by test program name;
-}

module AutoBench.Internal.Types
  (

  -- * Re-exports
    module AutoBench.Types
  -- * User inputs
  -- ** Test data options
  , toHRange               -- Convert @Gen l s u :: DataOpts@ to a Haskell range.
  -- ** Internal representation of user inputs
  , UserInputs(..)         -- A data structure maintained by the system to classify user inputs.
  , initUserInputs         -- Initialise a 'UserInputs' data structure.
  -- * Benchmarking
  , BenchReport(..)        -- A report to summarise the benchmarking phase of testing.
  , Coord                  -- (Input Size, Runtime) results as coordinates for unary test programs.
  , Coord3                 -- (Input Size, Input Size, Runtime) results as coordinates for binary test programs.
  , DataSize(..)           -- The size of unary and binary test data.
  , SimpleReport(..)       -- A simplified version of Criterion's 'Report'. See 'Criterion.Types.Report'.
  -- * Test results
  , TestReport(..)         -- A report to summarise the system's testing phase.
  -- ** QuickBench
  , QuickReport(..)        -- A report to summarise the QuickBench phase of testing.
  -- * Statistical analysis
  , AnalysisReport(..)     -- A report to summarise the system's analysis phase.
  , CVStats(..)            -- Fitting statistics calculated for regression models per each iteration of cross-validation.
  , Improvement            -- An efficiency improvement is an ordering between two test programs and a rating
                           -- 0 <= d <= 1 that corresponds to the percentage of test cases that support the ordering.
  , Exp                    -- Expressions with 'Double' literals.
  , LinearCandidate(..)    -- The details of a regression model necessary to fit it to a given dataset.
  , LinearFit(..)          -- A regression model's fitting statistics and helper functions: predicting y-coordinates, pretty printing.
  , SimpleResults(..)      -- Simple statistical analysis results for each test program.
  , numPredictors          -- Number of predictors for each type of model.
  , simpleReportToCoord    -- Convert a 'SimpleReport' to a (input size(s), runtime) coordinate, i.e., 'Coord' or 'Coord3'.
  , simpleReportsToCoords  -- Convert a list of 'SimpleReport' to a (input size(s), runtime) coordinate, i.e., 'Coord' or 'Coord3'.
  -- ** QuickBench
  , QuickAnalysis(..)      -- A report to summarise the system's analysis phase for QuickBenching.
  , QuickResults(..)       -- Simple statistical analysis results for each test program for QuickBenching.
  -- * Errors
  -- ** System errors
  , SystemError(..)        -- System errors.
  -- ** Input errors
  , InputError(..)         -- User input errors.

  ) where

import           Control.Exception.Base (Exception)
import           Criterion.Types        (OutlierEffect(..))
import           Data.Either            (partitionEithers)
import           Numeric.LinearAlgebra  (Vector)

import qualified AutoBench.Internal.Expr as E
import           AutoBench.Types  -- Re-export.

import AutoBench.Internal.AbstractSyntax
  ( HsType
  , Id
  , ModuleElem(..)
  , TypeString
  )

-- * User inputs

-- ** Test suites

-- | Convert @Gen l s u :: DataOpts@ to a Haskell range.
toHRange :: DataOpts -> [Int]
toHRange Manual{}    = []
toHRange (Gen l s u) = [l, (l + s) .. u]

-- ** Internal representation of user inputs

-- | While user inputs are being analysed by the system, a 'UserInputs' data
-- structure is maintained. The purpose of this data structure is to classify
-- user inputs according to the properties they satisfy. For example, when the
-- system first interprets a user input file, all of its definitions are added
-- to the '_allElems' list. This list is then processed to determine which
-- definitions have function types that are syntactically compatible with the
-- requirements of the system (see 'AutoBench.Internal.StaticChecks').
-- Definitions that are compatible are added to the '_validElems' list, and
-- those that aren't are added to the '_invalidElems' list. Elements in the
-- '_validElems' list are then classified according to, for example, whether
-- they are nullary, unary, or binary functions. This check process continues
-- until all user  inputs are classified according to the list headers below.
-- Note that both static ('AutoBench.Internal.StaticChecks') and dynamic
-- ('AutoBench.Internal.DynamicChecks') checks are required to classify user
-- inputs.
--
-- Notice that each /invalid/ definitions has one or more input errors
-- associated with it.
--
-- After the system has processed all user inputs, users can review this data
-- structure to see how the system has classified their inputs, and if any
-- input errors have been generated.
data UserInputs =
  UserInputs
   {
     _allElems           :: [(ModuleElem, Maybe TypeString)]         -- ^ All definitions in a user input file.
   , _invalidElems       :: [(ModuleElem, Maybe TypeString)]         -- ^ Syntactically invalid definitions (see 'AutoBench.Internal.AbstractSyntax').
   , _validElems         :: [(Id, HsType)]                           -- ^ Syntactically valid definitions (see 'AutoBench.Internal.AbstractSyntax').
   , _nullaryFuns        :: [(Id, HsType)]                           -- ^ Nullary functions.
   , _unaryFuns          :: [(Id, HsType)]                           -- ^ Unary functions.
   , _binaryFuns         :: [(Id, HsType)]                           -- ^ Binary functions.
   , _arbFuns            :: [(Id, HsType)]                           -- ^ Unary/binary functions whose input types are members of the Arbitrary type class.
   , _benchFuns          :: [(Id, HsType)]                           -- ^ Unary/binary functions whose input types are members of the NFData type class.
   , _nfFuns             :: [(Id, HsType)]                           -- ^ Unary/binary functions whose result types are members of the NFData type class.
   , _invalidData        :: [(Id, HsType, [InputError])]             -- ^ Invalid user-specified test data.
   , _unaryData          :: [(Id, HsType, [Int])]                    -- ^ Valid user-specified test data for unary functions /with size information/.
   , _binaryData         :: [(Id, HsType, [(Int, Int)])]             -- ^ Valid user-specified test data for binary functions /with size information/.
   , _invalidTestSuites  :: [(Id, [InputError])]                     -- ^ Invalid test suites.
   , _testSuites         :: [(Id, TestSuite)]                        -- ^ Valid test suites.
   }

-- | Initialise a 'UserInputs' data structure by specifying the '_allElems'
-- list.
initUserInputs :: [(ModuleElem, Maybe TypeString)] -> UserInputs
initUserInputs xs =
  UserInputs
    {
      _allElems          = xs
    , _invalidElems      = []
    , _validElems        = []
    , _nullaryFuns       = []
    , _unaryFuns         = []
    , _binaryFuns        = []
    , _arbFuns           = []
    , _benchFuns         = []
    , _nfFuns            = []
    , _invalidData       = []
    , _unaryData         = []
    , _binaryData        = []
    , _invalidTestSuites = []
    , _testSuites        = []
    }

-- * Benchmarking

-- | (Input Size, Runtime) results as coordinates for unary test programs.
type Coord = (Double, Double)

-- | (Input Size, Input Size, Runtime) results as coordinates for binary test
-- programs.
type Coord3 = (Double, Double, Double)

-- | The size of unary and binary test data.
data DataSize =
    SizeUn Int       -- ^ The size of unary test data.
  | SizeBin Int Int  -- ^ The size of binary test data.
    deriving (Ord, Eq)

-- | A report to summarise the benchmarking phase of testing.
data BenchReport =
  BenchReport
    {
      _reports   :: [[SimpleReport]]  -- ^ Individual reports for each test case, per test program.
    , _baselines :: [SimpleReport]    -- ^ Baseline measurements (will be empty if '_baseline' is set to @False@).
    }

-- | A simplified version of Criterion's 'Report' datatype, see
-- 'Criterion.Types.Report'.
data SimpleReport =
  SimpleReport
   {
     _name       :: Id              -- ^ Name of test program.
   , _size       :: DataSize        -- ^ Size of test data.
   , _samples    :: Int             -- ^ Number of samples used to calculate statistics below.
   , _runtime    :: Double          -- ^ Estimate runtime.
   , _stdDev     :: Double          -- ^ Estimate standard deviation.
   , _outVarEff  :: OutlierEffect   -- ^ Outlier effect.
   , _outVarFrac :: Double          -- ^ Outlier effect as a percentage.
   }

-- * Test results

-- | A report to summarise the system's testing phase.
data TestReport =
  TestReport
    {
      _tProgs    :: [String]          -- ^ Names of all test programs.
    , _tDataOpts :: DataOpts          -- ^ Which test data options were used.
    , _tNf       :: Bool              -- ^ Whether test cases were evaluated to normal form.
    , _tGhcFlags :: [String]          -- ^ Flags used when compiling the benchmarking file.
    , _eql       :: Bool              -- ^ Whether test programs are semantically equal according to QuickCheck testing.
    , _br        :: BenchReport       -- ^ Benchmarking report.
    }

-- ** QuickBench

-- | A report to summarise the QuickBench testing phase.
data QuickReport =
  QuickReport
    {
      _qName     :: Id                        -- ^ Name of test program.
    , _qRuntimes :: Either [Coord] [Coord3]   -- ^ [(Input size(s), mean runtime)].
    } deriving Show

-- * Statistical analysis

-- | A report to summarise the system's analysis phase.
data AnalysisReport =
  AnalysisReport
    {
      _anlys :: [SimpleResults]       -- ^ Simple statistical analysis results per test program.
    , _imps  :: [Improvement]         -- ^ Improvement results.
    , _blAn  :: Maybe SimpleResults   -- ^ Analysis of baseline measurements, if applicable.
    }

-- | Simple statistical analysis results for each test program.
data SimpleResults =
  SimpleResults
   {
     _srIdt           :: Id                           -- ^ Name of test program.
   , _srRaws          :: Either [Coord] [Coord3]      -- ^ Raw input size/runtime results.
   , _srStdDev        :: Double                       -- ^ Standard deviation of all runtime results.
   , _srAvgOutVarEff  :: OutlierEffect                -- ^ Average outlier effect.
   , _srAvgPutVarFrac :: Double                       -- ^ Average outlier effect as a percentage.
   , _srFits          :: [LinearFit]                  -- ^ Fitting statistics for each candidate model.
   }

-- | An efficiency improvement is an ordering between two test programs and a
-- rating 0 <= d <= 1 that corresponds to the percentage of test cases that
-- support the ordering.
type Improvement = (Id, Ordering, Id, Double)

-- | Fitting statistics calculated for regression models per each iteration of
-- cross-validation. Cumulative fitting statistics are produced by combining
-- 'CVStats' from all iterations, for example, PMSE and PMAE. See 'Stats'.
data CVStats =
  CVStats
   {
     _cv_mse    :: Double   -- ^ Mean squared error.
   , _cv_mae    :: Double   -- ^ Mean absolute error.
   , _cv_ss_tot :: Double   -- ^ Total sum of squares.
   , _cv_ss_res :: Double   -- ^ Residual sum of squares.
   } deriving Eq

-- | Expressions with 'Double' literals.
type Exp = E.Expr Double

-- | Each 'LinearType' gives rise to a 'LinearCandidate' that is then fitted to
-- a given data set generating a 'LinearFit'. Unlike a 'LinearType', which
-- just describes a particular regression model, a 'LinearCandidate'
-- encompasses the required information to fit a model to a given data set.
-- For example, it includes '_fxs' to transforms the raw x-coordinates
-- of the dataset before fitting the model, and '_fyhat' which can be used to
-- generate y-coordinates predicted by the model once it has been fit.
--
-- For example, if fitting a 'Log b 1' model, '_fxs' will transform each
-- x-coordinate in the data set to log_b(x) before fitting. Then the linear
-- relationship between the /resulting/ xy-coordinates corresponds to a
-- logarithmic relationship between the /initial/ xy-coordinates.
--
-- When the coefficients of a model are determined by regression analysis,
-- the corresponding 'LinearCandidate' gives rise to a 'LinearFit'.
data LinearCandidate =
  LinearCandidate
   {
     _lct    :: LinearType                                         -- ^ The model.
   , _fxs    :: Vector Double -> Vector Double                     -- ^ A function to transform x-coords before fitting.
   , _fex    :: Vector Double -> Exp                               -- ^ A function to generate the model's equation as an 'Exp'.
   , _fyhat  :: Vector Double -> Vector Double -> Vector Double    -- ^ A function to generate model's predicted y-coords.
   }

-- | When a 'LinearCandidate' is fitted to a given data set and its coefficients
-- determined by regression analysis, a 'LinearFit' is generated.
-- A 'LinearFit' primarily includes the fitting statistics ('Stats') of the
-- model to be used to compare it against other models also fitted to the same
-- data set. In order to be able to plot a 'LinearFit' on a results graph
-- as a line of best fit, the '_yhat' function generates y-coordinates
-- predicted by the model for a given set of x-coordinates. To pretty print the
-- 'LinearFit', the '_ex' function generates an 'Exp' expression for the model's
-- equation that has a pretty printing function.
data LinearFit =
  LinearFit
   {
     _lft  :: LinearType                       -- ^ The model.
   , _cfs  :: Vector Double                    -- ^ The coefficients of the model.
   , _ex   :: Exp                              -- ^ The model's equation as an 'Exp'.
   , _yhat :: Vector Double -> Vector Double   -- ^ A function to generate the model's predicted y-coords for a given set of x-coords.
   , _sts  :: Stats                            -- ^ Fitting statistics.
   }

-- | Number of predictors for each type of model.
numPredictors :: LinearType -> Int
numPredictors (Poly      k) = k + 1
numPredictors (Log     _ k) = k + 1
numPredictors (PolyLog _ k) = k + 1
numPredictors Exp{}         = 2

-- | Convert a list of 'SimpleReport's to a list of (input size(s), runtime)
-- coordinates, i.e., a list 'Coord's or 'Coord3's. The name of each
-- simple report is verified against the given test program identifier.
simpleReportsToCoords :: Id -> [SimpleReport] -> Either [Coord] [Coord3]
simpleReportsToCoords idt srs = case (cs, cs3) of
  ([], _) -> Right cs3
  (_, []) -> Left cs
  _       -> Left [] -- Shouldn't happen.
  where
    srs' = filter (\sr -> _name sr == idt) srs
    (cs, cs3) = partitionEithers (fmap simpleReportToCoord srs')

-- | Convert a 'SimpleReport' to a (input size(s), runtime) coordinate,
-- i.e., 'Coord' or 'Coord3'.
simpleReportToCoord :: SimpleReport -> Either Coord Coord3
simpleReportToCoord sr = case _size sr of
  SizeUn n      -> Left  (fromIntegral n, _runtime sr)
  SizeBin n1 n2 -> Right (fromIntegral n1, fromIntegral n2, _runtime sr)

-- ** QuickBench

-- | A report to summarise the system's analysis phase for QuickBenching.
data QuickAnalysis =
  QuickAnalysis
    {
      _qAnlys :: [QuickResults]    -- ^ Quick results per test program.
    , _qImps  :: [Improvement]     -- ^ Improvement results.
    }

-- | Simple statistical analysis results for each test program for QuickBenching.
data QuickResults =
  QuickResults
   {
     _qrIdt   :: Id                           -- ^ Name of test program.
   , _qrRaws  :: Either [Coord] [Coord3]      -- ^ Raw input size/runtime results.
   , _qrFits  :: [LinearFit]                  -- ^ Fitting statistics for each candidate model.
   }

-- * Errors

-- | Errors raised by the system due to implementation failures. These can be
-- generated at any time but are usually used to report unexpected IO results.
-- For example, when dynamically checking user inputs (see
-- 'AutoBench.Internal.UserInputChecks'), system errors are used to relay
-- 'InterpreterError's thrown by functions in the hint package in cases
-- where the system didn't expect errors to result.
data SystemError = InternalErr String

-- Note: needed for the 'Exception' instance.
instance Show SystemError where
  show (InternalErr s) = "Internal error: " ++ s ++ "\n** Please report on GitHub **"

instance Exception SystemError

-- ** Input errors

-- | Input errors are generated by the system while analysing user input
-- files. Examples input errors include erroneous test options, invalid test
-- data, and test programs with missing Arbitrary/NFData instances.
--
-- In general, the system always attempts to continue with its execution for as
-- long as possible. Therefore, unless a critical error is encountered, such as
-- a filepath or file access error, it will collate all non-critical input
-- errors. These will then be summarised after the user input file has been
-- fully analysed.
data InputError =
    FilePathErr   String    -- ^ Invalid filepath.
  | FileErr       String    -- ^ File access error.
  | TestSuiteErr  String    -- ^ Invalid test suite.
  | DataOptsErr   String    -- ^ Invalid data options.
  | AnalOptsErr   String    -- ^ Invalid statistical analysis options.
  | TypeErr       String    -- ^ Invalid type signature.
  | InstanceErr   String    -- ^ One or more missing instances
  | TestReportErr String    -- ^ Invalid test report.
  | QuickOptsErr  String    -- ^ Invalid quick options. See AutoBench.QuickBench.
  | QuickBenchErr String    -- ^ QuickBench error. See AutoBench.QuickBench.

-- Note: needed for the 'Exception' instance.
instance Show InputError where
  show (FilePathErr   s) = "File path error: "        ++ s
  show (FileErr       s) = "File error: "             ++ s
  show (TestSuiteErr  s) = "Test suite error: "       ++ s
  show (DataOptsErr   s) = "Test data error: "        ++ s
  show (AnalOptsErr   s) = "Analysis options error: " ++ s
  show (TypeErr       s) = "Type error: "             ++ s
  show (InstanceErr   s) = "Instance error: "         ++ s
  show (TestReportErr s) = "Test report error: "      ++ s
  show (QuickOptsErr  s) = "Quick options error: "    ++ s
  show (QuickBenchErr s) = "QuickBench test error: "  ++ s

instance Exception InputError

ext: hs

max_stars_repo_name: recursion-ninja/AutoBench

max_issues_repo_name: recursion-ninja/AutoBench

max_forks_repo_name: recursion-ninja/AutoBench


text:
{-# LANGUAGE CPP                  #-}
{-# LANGUAGE FlexibleContexts     #-}
{-# LANGUAGE FlexibleInstances    #-}
{-# LANGUAGE ScopedTypeVariables  #-}
{-# LANGUAGE UndecidableInstances #-}

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

-----------------------------------------------------------------------------
{- |
Module      :  Numeric.LinearAlgebra.Tests.Instances
Copyright   :  (c) Alberto Ruiz 2008
License     :  BSD3
Maintainer  :  Alberto Ruiz
Stability   :  provisional

Arbitrary instances for vectors, matrices.

-}

module Test.Numeric.LinearAlgebra.Tests.Instances(
    Sq(..),     rSq,cSq,
    Sq2WC(..),  rSq2WC,cSq2WC,
    Rot(..),    rRot,cRot,
                rHer,cHer,
    WC(..),     rWC,cWC,
    SqWC(..),   rSqWC, cSqWC, rSymWC, cSymWC,
    PosDef(..), rPosDef, cPosDef,
    Consistent(..), rConsist, cConsist,
    RM,CM, rM,cM,
    FM,ZM, fM,zM
) where

import           System.Random

import           Control.Monad                 (replicateM)
import           Numeric.LinearAlgebra.HMatrix hiding (vector)
import           Test.QuickCheck               (Arbitrary, arbitrary, choose, shrink,
                                                sized, vector)

import           Data.Proxy                    (Proxy (..))
import           GHC.TypeLits
import qualified Numeric.LinearAlgebra.Static  as Static
#if MIN_VERSION_base(4,11,0)
import           Prelude                       hiding ((<>))
#endif

shrinkListElementwise :: (Arbitrary a) => [a] -> [[a]]
shrinkListElementwise []     = []
shrinkListElementwise (x:xs) = [ y:xs | y  <- shrink x                 ]
                            ++ [ x:ys | ys <- shrinkListElementwise xs ]

shrinkPair :: (Arbitrary a, Arbitrary b) => (a,b) -> [(a,b)]
shrinkPair (a,b) = [ (a,x) | x <- shrink b ] ++ [ (x,b) | x <- shrink a ]

chooseDim = sized $ \m -> choose (1,max 1 m)

instance (Field a, Arbitrary a) => Arbitrary (Vector a) where
    arbitrary = do m <- chooseDim
                   l <- vector m
                   return $ fromList l
    -- shrink any one of the components
    shrink = map fromList . shrinkListElementwise . toList

instance KnownNat n => Arbitrary (Static.R n) where
    arbitrary = do
      l <- vector n
      return (Static.fromList l)

      where
        n :: Int
        n = fromIntegral (natVal (Proxy :: Proxy n))

    shrink _v = []

instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where
    arbitrary = do
        m <- chooseDim
        n <- chooseDim
        l <- vector (m*n)
        return $ (m><n) l

    -- shrink any one of the components
    shrink a = map (rows a >< cols a)
               . shrinkListElementwise
               . concat . toLists
                     $ a

instance (KnownNat n, KnownNat m) => Arbitrary (Static.L m n) where
    arbitrary = do
      l <- vector (m * n)
      return (Static.fromList l)

      where
        m :: Int
        m = fromIntegral (natVal (Proxy :: Proxy m))

        n :: Int
        n = fromIntegral (natVal (Proxy :: Proxy n))

    shrink _mat = []

-- a square matrix
newtype (Sq a) = Sq (Matrix a) deriving Show
instance (Element a, Arbitrary a) => Arbitrary (Sq a) where
    arbitrary = do
        n <- chooseDim
        l <- vector (n*n)
        return $ Sq $ (n><n) l

    shrink (Sq a) = [ Sq b | b <- shrink a ]

-- a pair of square matrices
newtype (Sq2WC a) = Sq2WC (Matrix a, Matrix a) deriving Show
instance (ArbitraryField a, Numeric a) => Arbitrary (Sq2WC a) where
    arbitrary = do
        n <- chooseDim
        l <- vector (n*n)
        r <- vector (n*n)
        l' <- makeWC $ (n><n) l
        r' <- makeWC $ (n><n) r
        return $ Sq2WC (l', r')
        where
            makeWC m = do
              let (u,_,v) = svd m
                  n = rows m
              sv' <- replicateM n (choose (1,100))
              let s = diag (fromList sv')
              return $ u <> real s <> tr v

-- a unitary matrix
newtype (Rot a) = Rot (Matrix a) deriving Show
instance (Field a, Arbitrary a) => Arbitrary (Rot a) where
    arbitrary = do
        Sq m <- arbitrary
        let (q,_) = qr m
        return (Rot q)

-- a complex hermitian or real symmetric matrix
instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Herm a) where
    arbitrary = do
        Sq m <- arbitrary
        let m' = m/2
        return $ sym m'

class (Field a, Arbitrary a, Element (RealOf a), Random (RealOf a)) => ArbitraryField a
instance ArbitraryField Float
instance ArbitraryField (Complex Float)

-- a well-conditioned general matrix (the singular values are between 1 and 100)
newtype (WC a) = WC (Matrix a) deriving Show
instance (Numeric a, ArbitraryField a) => Arbitrary (WC a) where
    arbitrary = do
        m <- arbitrary
        let (u,_,v) = svd m
            r = rows m
            c = cols m
            n = min r c
        sv' <- replicateM n (choose (1,100))
        let s = diagRect 0 (fromList sv') r c
        return $ WC (u <> real s <> tr v)

-- a well-conditioned square matrix (the singular values are between 1 and 100)
newtype (SqWC a) = SqWC (Matrix a) deriving Show
instance (ArbitraryField a, Numeric a) => Arbitrary (SqWC a) where
    arbitrary = do
        Sq m <- arbitrary
        let (u,_,v) = svd m
            n = rows m
        sv' <- replicateM n (choose (1,100))
        let s = diag (fromList sv')
        return $ SqWC (u <> real s <> tr v)

-- a positive definite square matrix (the eigenvalues are between 0 and 100)
newtype (PosDef a) = PosDef (Matrix a) deriving Show
instance (Numeric a, ArbitraryField a, Num (Vector a))
    => Arbitrary (PosDef a) where
    arbitrary = do
        m <- arbitrary
        let (_,v) = eigSH m
            n = rows (unSym m)
        l <- replicateM n (choose (0,100))
        let s = diag (fromList l)
            p = v <> real s <> tr v
        return $ PosDef (0.5 * p + 0.5 * tr p)

-- a pair of matrices that can be multiplied
newtype (Consistent a) = Consistent (Matrix a, Matrix a) deriving Show
instance (Field a, Arbitrary a) => Arbitrary (Consistent a) where
    arbitrary = do
        n <- chooseDim
        k <- chooseDim
        m <- chooseDim
        la <- vector (n*k)
        lb <- vector (k*m)
        return $ Consistent ((n><k) la, (k><m) lb)

    shrink (Consistent (x,y)) = [ Consistent (u,v) | (u,v) <- shrinkPair (x,y) ]

type RM = Matrix Float
type CM = Matrix (Complex Float)
type FM = Matrix Float
type ZM = Matrix (Complex Float)

rM m = m :: RM
cM m = m :: CM
fM m = m :: FM
zM m = m :: ZM

rHer m = unSym m :: RM
cHer m = unSym m :: CM

rRot (Rot m) = m :: RM
cRot (Rot m) = m :: CM

rSq  (Sq m)  = m :: RM
cSq  (Sq m)  = m :: CM

rSq2WC (Sq2WC (a, b)) = (a, b) :: (RM, RM)
cSq2WC (Sq2WC (a, b)) = (a, b) :: (CM, CM)

rWC (WC m) = m :: RM
cWC (WC m) = m :: CM

rSqWC (SqWC m) = m :: RM
cSqWC (SqWC m) = m :: CM

rSymWC (SqWC m) = sym m :: Herm R
cSymWC (SqWC m) = sym m :: Herm C

rPosDef (PosDef m) = m :: RM
cPosDef (PosDef m) = m :: CM

rConsist (Consistent (a,b)) = (a,b::RM)
cConsist (Consistent (a,b)) = (a,b::CM)

max_stars_repo_name: schnecki/hmatrix-float

max_issues_repo_name: schnecki/hmatrix-float

max_forks_repo_name: schnecki/hmatrix-float


text:
{-# LANGUAGE BangPatterns #-}

module School.Types.Decoding
( binToDouble
, binToInt
, getDouble
, getInt
, binToMatrixDouble
, binToListInt
, runGet
) where

import Control.Monad (replicateM)
import Data.ByteString (ByteString, length, splitAt, unpack)
import Data.Int (Int32)
import Data.Serialize.Get (Get, getInt32be, runGet)
import Data.Serialize.IEEE754 (getFloat64be)
import Numeric.LinearAlgebra ((><), Matrix, I, R)
import Prelude hiding (length, splitAt)
import School.Types.DataType (DataType(..))
import School.Types.Error (Error)

getDouble :: Get R
getDouble = getFloat64be

binToDouble :: ByteString -> Either Error R
binToDouble = runGet getDouble

binToMatrixDouble :: DataType
                  -> Int
                  -> Int
                  -> (ByteString -> Either Error
                                           (Matrix R))
binToMatrixDouble DBL64B nRows nCols =
  runGet (getDoubleMatrixDouble nRows nCols)
binToMatrixDouble INT08B nRows nCols =
   Right
 . (nRows >< nCols)
 . map (fromIntegral . fromEnum)
 . unpack
binToMatrixDouble dType _ _ = const . Left $
  "Decoding to Matrix Double undefined for " ++ show dType

getDoubleMatrixDouble :: Int
                      -> Int
                      -> Get (Matrix Double)
getDoubleMatrixDouble nRows nCols = do
  let nElements = nRows * nCols
  list <- replicateM nElements getDouble
  return $ (nRows >< nCols) list

getInt :: Get Int32
getInt = getInt32be

binToInt :: ByteString -> Either Error I
binToInt = fmap fromIntegral . runGet getInt

binToListInt :: DataType
             -> (ByteString -> Either Error [Int])
binToListInt DBL64B = const . Left $
  "Reject conversion from floating point DBL64B to integral"
binToListInt INT32B = loop [] where
  loop acc bytes = let len = length bytes in
    if len < 4
      then return acc
      else do
        let !(next, current) = splitAt (len - 4) bytes
        int <- fromIntegral <$> binToInt current
        loop (int:acc) next
binToListInt INT08B =
    Right
  . map (fromIntegral . fromEnum)
  . unpack

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


text:
import Data.List (transpose)
import Data.Complex

type Matrix a = [[a]]

main :: IO ()
main =
    mapM_ (\a -> do
        putStrLn "\nMatrix:"
        mapM_ print a
        putStrLn "Conjugate Transpose:"
        mapM_ print (conjTranspose a)
        putStrLn $ "Hermitian? " ++ show (isHermitianMatrix a)
        putStrLn $ "Normal? " ++ show (isNormalMatrix a)
        putStrLn $ "Unitary? " ++ show (isUnitaryMatrix a))
        ([[[3,       2:+1],
           [2:+(-1), 1   ]],

          [[1, 1, 0],
           [0, 1, 1],
           [1, 0, 1]],

          [[sqrt 2/2:+0, sqrt 2/2:+0,     0   ],
           [0:+sqrt 2/2, 0:+ (-sqrt 2/2), 0   ],
           [0,           0,               0:+1]]] :: [Matrix (Complex Double)])

isHermitianMatrix, isNormalMatrix, isUnitaryMatrix :: RealFloat a => Matrix (Complex a) -> Bool
isHermitianMatrix a = a `approxEqualMatrix` conjTranspose a
isNormalMatrix a = (a `mmul` conjTranspose a) `approxEqualMatrix` (conjTranspose a `mmul` a)
isUnitaryMatrix a = (a `mmul` conjTranspose a) `approxEqualMatrix` ident (length a)

approxEqualMatrix :: (Fractional a, Ord a) => Matrix (Complex a) -> Matrix (Complex a) -> Bool
approxEqualMatrix a b = length a == length b && length (head a) == length (head b) &&
                        and (zipWith approxEqualComplex (concat a) (concat b))
    where approxEqualComplex (rx :+ ix) (ry :+ iy) = abs (rx - ry) < eps && abs (ix - iy) < eps
          eps = 1e-14

mmul :: Num a => Matrix a -> Matrix a -> Matrix a
mmul a b = [[sum (zipWith (*) row column) | column <- transpose b] | row <- a]

ident :: Num a => Int -> Matrix a
ident size = [[fromIntegral $ div a b * div b a | a <- [1..size]] | b <- [1..size]]

conjTranspose :: Num a => Matrix (Complex a) -> Matrix (Complex a)
conjTranspose = map (map conjugate) . transpose

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE FlexibleContexts  #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MonoLocalBinds    #-}
{-# LANGUAGE TypeOperators     #-}
module Utils.Array where

import           Control.Monad        (when)
import           Data.Array.Repa      as R
import           Data.Binary
import           Data.ByteString.Lazy as BL
import           Data.Complex
import           Data.Int
import           Data.List            as L
import           Data.Vector.Unboxed  as VU
import           System.IO            (IOMode (..), withBinaryFile)
import           Text.Printf
import           Types

{-# INLINE rotate3D #-}
rotate3D :: (R.Source s e) => R.Array s DIM3 e -> R.Array D DIM3 e
rotate3D arr =
  let (Z :. k :. i :. j) = extent arr
   in R.backpermute
        (Z :. i :. j :. k)
        (\(Z :. a :. b :. c) -> (Z :. c :. a :. b))
        arr

{-# INLINE rotate3DR #-}
rotate3DR :: (R.Source s e) => R.Array s DIM3 e -> R.Array D DIM3 e
rotate3DR arr =
  let (Z :. i :. j :. k) = extent arr
   in R.backpermute
        (Z :. k :. i :. j)
        (\(Z :. a :. b :. c) -> (Z :. b :. c :. a))
        arr

{-# INLINE rotate4D #-}
rotate4D :: (R.Source s e) => R.Array s DIM4 e -> R.Array D DIM4 e
rotate4D arr =
  let (Z :. k :. l :. i :. j) = extent arr
   in R.backpermute
        (Z :. l :. i :. j :. k)
        (\(Z :. a :. b :. c :. d) -> (Z :. d :. a :. b :. c))
        arr

{-# INLINE rotate4D' #-}
rotate4D' :: (R.Source s e) => R.Array s DIM4 e -> R.Array D DIM4 e
rotate4D' arr =
  let (Z :. k :. l :. i :. j) = extent arr
   in R.backpermute
        (Z :. j :. k :. l :. i)
        (\(Z :. a :. b :. c :. d) -> (Z :. b :. c :. d :. a))
        arr

{-# INLINE rotate4D2 #-}
rotate4D2 :: (R.Source s e) => R.Array s DIM4 e -> R.Array D DIM4 e
rotate4D2 arr =
  let (Z :. k :. l :. i :. j) = extent arr
  in R.backpermute
       (Z :. i :. j :. k :. l)
       (\(Z :. a :. b :. c :. d) -> (Z :. c :. d :. a :. b))
       arr

{-# INLINE rotateR2Z1T0Array #-}
rotateR2Z1T0Array :: (R.Source s e) => R.Array s DIM4 e -> R.Array D DIM4 e
rotateR2Z1T0Array arr =
  let (Z :. k :. l :. i :. j) = extent arr
   in R.backpermute
        (Z :. k :. i :. j :. l)
        (\(Z :. a :. b :. c :. d) -> (Z :. a :. d :. b :. c))
        arr

{-# INLINE rotateR2Z2T0S0Array #-}
rotateR2Z2T0S0Array :: (R.Source s e) => R.Array s DIM6 e -> R.Array D DIM6 e
rotateR2Z2T0S0Array arr =
  let (Z :. k :. l :. m :. n :. i :. j) = extent arr
  in R.backpermute
       (Z :. k :. l :. i :. j :. m :. n)
       (\(Z :. a :. b :. c :. d :. e :. f) -> (Z :. a :. b :. e :. f :. c :. d))
       arr

{-# INLINE rotateR2Z2T0S0ArrayDeconv #-}
rotateR2Z2T0S0ArrayDeconv :: (R.Source s e) => R.Array s DIM6 e -> R.Array D DIM6 e
rotateR2Z2T0S0ArrayDeconv arr =
  let (Z :. k :. l :. m :. n :. i :. j) = extent arr
  in R.backpermute
       (Z :. m :. n :. i :. j :. k :. l)
       (\(Z :. a :. b :. c :. d :. e :. f) -> (Z :. e :. f :. a :. b :. c :. d))
       arr

{-# INLINE writeRepaArray #-}
writeRepaArray ::
     (Shape sh, Binary e, R.Source s e) => FilePath -> R.Array s sh e -> IO ()
writeRepaArray filePath arr =
  withBinaryFile filePath WriteMode $ \h -> do
    let shapeBS = encode . listOfShape . extent $ arr
        elemBS = encode . R.toList $ arr
    BL.hPut h . encode . BL.length $ shapeBS -- ByteString lenght is 64 bits = 8 bytes
    BL.hPut h shapeBS
    BL.hPut h . encode . BL.length $ elemBS  -- ByteString lenght is 64 bits = 8 bytes
    BL.hPut h elemBS

{-# INLINE readRepaArray #-}
readRepaArray ::
     (R.Source U e, Binary e, Unbox e, Shape sh)
  => FilePath
  -> IO (R.Array U sh e)
readRepaArray filePath =
  withBinaryFile filePath ReadMode $ \h -> do
    shapeLen <- decode <$> BL.hGet h 8 :: IO Int64
    sh <- (shapeOfList . decode) <$> BL.hGet h (fromIntegral shapeLen)
    elemLen <- decode <$> BL.hGet h 8 :: IO Int64
    elem <- decode <$> BL.hGet h (fromIntegral elemLen)
    when
      (size sh /= L.length elem)
      (error
         (printf
            "readRepaArray: array shape size (%d) /= element size (%d)\n"
            (size sh)
            (L.length elem)))
    return $ fromListUnboxed sh elem

{-# INLINE pad #-}
pad ::
     (Source s e, Shape sh)
  => [Int]
  -> e
  -> R.Array s sh e
  -> R.Array D sh e
pad newDims padVal arr
  | L.length newDims /= L.length oldDimList =
    error $
    printf
      "pad: Dimensions are different. (%d vs. %d)"
      (L.length oldDimList)
      (L.length newDims)
  | L.all (== 0) diff = delay arr
  | otherwise =
    backpermuteDft
      (fromFunction (shapeOfList dimList) (const padVal))
      (\sh' ->
         let idx = L.zipWith (-) (listOfShape sh') diff
         in if L.or (L.zipWith (\i j -> i < 0 || (i >= j)) idx oldDimList)
              then Nothing
              else Just $ shapeOfList idx)
      arr
  where
    !oldDimList = listOfShape . extent $ arr
    !dimList = L.zipWith max newDims oldDimList
    !diff =
      L.zipWith
        (\a b ->
           if a - b <= 0
             then 0
             else if a - b == 1
                    then 1
                    else div (a - b) 2)
        newDims
        oldDimList

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


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

{-| Fast FFTs using FFTW -}
module SDR.FFT (
    -- * Windows
    hamming,
    hanning,
    blackman,

    -- * FFTs
    fftw',
    fftw,
    fftwReal',
    fftwReal,
    fftwParallel
    ) where

import           Control.Monad                as CM
import           Foreign.Storable
import           Foreign.Storable.Complex
import           Foreign.C.Types
import           Data.Complex
import           Foreign.ForeignPtr
import           Control.Concurrent           hiding (yield)
import qualified Data.Map                     as Map
import           Data.Coerce

import qualified Data.Vector.Generic          as VG
import qualified Data.Vector.Storable         as VS
import qualified Data.Vector.Storable.Mutable as VSM

import           Pipes
import           Numeric.FFTW

import           SDR.FilterDesign
import           SDR.VectorUtils

mallocForeignBufferAligned :: forall a. Storable a => Int -> IO (ForeignPtr a)
mallocForeignBufferAligned elems = do
    ptr <- fftwMalloc $ fromIntegral $ elems * sizeOf (undefined :: a)
    newForeignPtr fftwFreePtr ptr

-- | Creates a function that performs a complex to complex DFT.
fftw' :: (VG.Vector v (Complex Double))
     => Int -- ^ The size of the input and output buffers
     -> IO (v (Complex Double) -> IO (VS.Vector (Complex Double)))
fftw' samples = do
    ina <- mallocForeignBufferAligned samples
    out <- mallocForeignBufferAligned samples

    plan <- withForeignPtr ina $ \ip ->
        withForeignPtr out $ \op ->
            planDFT1d samples ip op Forward fftwEstimate

    return $ \inv' -> do
        out <- mallocForeignBufferAligned samples
        ina <- mallocForeignBufferAligned samples
        let inv = VSM.unsafeFromForeignPtr0 ina samples

        copyInto inv inv'

        let (fp, offset, length) = VSM.unsafeToForeignPtr inv

        withForeignPtr (coerce fp) $ \fpp ->
            withForeignPtr (coerce out) $ \op ->
                executeDFT plan fpp op

        return $ VS.unsafeFromForeignPtr0 out samples

-- | Creates a Pipe that performs a complex to complex DFT.
fftw :: (VG.Vector v (Complex Double))
     => Int -- ^ The size of the input and output buffers
     -> IO (Pipe (v (Complex Double)) (VS.Vector (Complex Double)) IO ())
fftw samples = do
    func <- fftw' samples
    return $ for cat $ \dat -> lift (func dat) >>= yield

-- | Creates a function that performs a real to complex DFT.
fftwReal' :: (VG.Vector v Double)
         => Int -- ^ The size of the input Vector
         -> IO (v Double -> IO (VS.Vector (Complex Double)))
fftwReal' samples = do
    --Allocate in and out buffers that wont be used because there doesnt seem to be a way to create a plan without them
    ina <- mallocForeignBufferAligned samples
    out <- mallocForeignBufferAligned samples

    plan <- withForeignPtr ina $ \ip ->
        withForeignPtr out $ \op ->
            planDFTR2C1d samples ip op fftwEstimate

    return $ \inv' -> do
        out <- mallocForeignBufferAligned ((samples `quot` 2) + 1)
        ina <- mallocForeignBufferAligned samples
        let inv = VSM.unsafeFromForeignPtr0 ina samples

        copyInto inv inv'
        let (fp, offset, length) = VSM.unsafeToForeignPtr inv

        withForeignPtr (coerce fp)  $ \fpp ->
            withForeignPtr (coerce out) $ \op ->
                executeDFTR2C plan fpp op

        return $ VS.unsafeFromForeignPtr0 out samples

-- | Creates a pipe that performs a real to complex DFT.
fftwReal :: (VG.Vector v Double)
         => Int -- ^ The size of the input Vector
         -> IO (Pipe (v Double) (VS.Vector (Complex Double)) IO ())
fftwReal samples = do
    func <- fftwReal' samples
    return $ for cat $ \dat -> lift (func dat) >>= yield

{-| Creates a pipe that uses multiple threads to perform complex to complex DFTs in
    a pipelined fashion. Each time a buffer is consumed, it is given to
    a pool of threads to perform the DFT. Then, if a thread has finished
    performing a previous DFT, the result is yielded.
-}
fftwParallel :: (VG.Vector v (Complex Double))
             => Int -- ^ The number of threads to use
             -> Int -- ^ The size of the input Vector
             -> IO (Pipe (v (Complex Double)) (VS.Vector (Complex Double)) IO ())
fftwParallel threads samples = do
    --plan the DFT
    ina <- mallocForeignBufferAligned samples
    out <- mallocForeignBufferAligned samples

    plan <- withForeignPtr ina $ \ip ->
        withForeignPtr out $ \op ->
            planDFT1d samples ip op Forward fftwEstimate

    --setup the channels and worker threads
    inChan <- newChan
    outMap <- newMVar Map.empty

    CM.replicateM threads $ forkIO $ forever $ do
        (idx, res) <- readChan inChan

        out <- mallocForeignBufferAligned samples
        ina <- mallocForeignBufferAligned samples
        let inv = VSM.unsafeFromForeignPtr0 ina samples

        copyInto inv res

        let (fp, offset, length) = VSM.unsafeToForeignPtr inv

        withForeignPtr (coerce fp) $ \fpp ->
            withForeignPtr (coerce out) $ \op ->
                executeDFT plan fpp op

        theMap <- takeMVar outMap
        putMVar outMap $ Map.insert idx (VS.unsafeFromForeignPtr0 out samples) theMap

    --build the pipe
    let pipe nextIn nextOut = do
            dat <- await
            lift $ writeChan inChan (nextIn, dat)

            theMap <- lift $ takeMVar outMap
            case Map.lookup nextOut theMap of
                Nothing  -> do
                    lift $ putMVar outMap theMap
                    pipe (nextIn + 1) nextOut
                Just dat -> do
                    lift $ putMVar outMap $ Map.delete nextOut theMap
                    yield dat
                    pipe (nextIn + 1) (nextOut + 1)

    return $ pipe 0 0

max_stars_repo_name: adamwalker/sdr

max_issues_repo_name: peixian/sdr

max_forks_repo_name: adamwalker/sdr


text:
{-# LANGUAGE CPP                 #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE PolyKinds           #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies        #-}
{-# LANGUAGE TypeOperators       #-}

{-|
Module      : Grenade.Core.WeightInitialization
Description : Defines the Weight Initialization methods of Grenade.
Copyright   : (c) Manuel Schneckenreither, 2018
License     : BSD2
Stability   : experimental

This module defines the weight initialization methods.

-}

module Grenade.Core.WeightInitialization
    ( getRandomVector
    , getRandomMatrix
    , WeightInitMethod (..)
    ) where

import           Control.Monad
import           Control.Monad.Primitive         (PrimBase, PrimState)
import           Data.Proxy
import           Data.Singletons.TypeLits
import           GHC.TypeLits                    hiding (natVal)
import           Numeric.LinearAlgebra.Static
import           System.Random.MWC
import           System.Random.MWC.Distributions

-- ^ Weight initialization method.
data WeightInitMethod
  = UniformInit -- ^ W_l,i ~ U(-1/sqrt(n_l),1/sqrt(n_l))                   where n_l is the number of nodes in layer l
  | Xavier      -- ^ W_l,i ~ U(-sqrt (6/n_l+n_{l+1}),sqrt (6/n_l+n_{l+1})) where n_l is the number of nodes in layer l
  | HeEtAl      -- ^ W_l,i ~ N(0,sqrt(2/n_l))                              where n_l is the number of nodes in layer l

-- | Get a random vector initialized according to the specified method.
getRandomVector ::
     forall m n. (PrimBase m, KnownNat n)
  => Integer
  -> Integer
  -> WeightInitMethod
  -> Gen (PrimState m)
  -> m (R n)
getRandomVector i o method gen = do
  unifRands <- vector <$> replicateM n (uniformR (-1, 1) gen)
  gaussRands <- vector <$> replicateM n (realToFrac <$> standard gen)
  return $
    case method of
      UniformInit -> (1 / sqrt (fromIntegral i)) * unifRands
      Xavier      -> (sqrt 6 / sqrt (fromIntegral i + fromIntegral o)) * unifRands
      HeEtAl      -> sqrt (2 / fromIntegral i) * gaussRands
  where
    n = fromIntegral $ natVal (Proxy :: Proxy n)

-- | Get a matrix with weights initialized according to the specified method.
getRandomMatrix ::
     forall m r n. (PrimBase m, KnownNat r, KnownNat n, KnownNat (n * r))
  => Integer
  -> Integer
  -> WeightInitMethod
  -> Gen (PrimState m)
  -> m (L r n)
getRandomMatrix i o method gen = do
  unifRands <- matrix <$> replicateM nr (uniformR (-1, 1) gen)
  gaussRands <- matrix <$> replicateM nr (realToFrac <$> standard gen)
  return $
    case method of
      UniformInit -> (1 / sqrt (fromIntegral i)) * unifRands
      Xavier      -> (sqrt 6 / sqrt (fromIntegral i + fromIntegral o)) * unifRands
      HeEtAl      -> sqrt (2 / fromIntegral i) * gaussRands
  where
    nr = fromIntegral $ natVal (Proxy :: Proxy (n * r))

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
module Types where

import Data.IORef
import Data.Ratio
import Data.Complex
import Data.Maybe (isJust)

import System.IO (Handle)
import Control.Monad.Error
import Text.ParserCombinators.Parsec (ParseError)

-- AST of the implemented scheme
data LispVal = Atom String
             | List [LispVal]
             | DottedList [LispVal] LispVal
             | Number Integer
             | Float Double
             | Complex (Complex Double)
             | Rational Rational
             | String String
             | Bool Bool
             | Character Char
             | PrimitiveFunc ([LispVal] -> ThrowsError LispVal)
             | IOFunc ([LispVal] -> IOThrowsError LispVal)
             | Port Handle
             | Func { params :: [String]
                    , vararg :: Maybe String
                    , body   :: [LispVal]
                    , closure:: Env
                    }

showVal :: LispVal -> String
showVal (Atom s)         = s
showVal (List l)         = "(" ++ showListVal l ++ ")"
showVal (DottedList h t) = "(" ++ showListVal h ++ " . " ++  showVal t ++ ")"
showVal (Number i)       = show i
showVal (Float d)        = show d
showVal (Rational r)     = show (numerator r) ++ "/" ++ show (denominator r)
showVal (String s)       = "\"" ++ s ++ "\""
showVal (Bool True)      = "#t"
showVal (Bool False)     = "#f"
showVal (Character c)    = [c]
showVal (Complex (real :+ img)) = show real ++ " + " ++ show img ++ "i"
showVal (PrimitiveFunc _)= "<primitive function>"
showVal (IOFunc _) = "<IO Primitive>"
showVal (Port _) = "<IO Port>"
showVal (Func args varargs _ _) =
  "(lambda (" ++ unwords (map show args) ++
  (case varargs of
     Nothing  -> ""
     Just arg -> " . " ++ arg) ++ ") ...)"

showListVal :: [LispVal] -> String
showListVal = unwords . map showVal

instance Show LispVal where show = showVal

-- Type to support errors during parsing/interpretation
data LispError = NumArgs Integer [LispVal]
               | TypeMismatch String LispVal
               | Parser ParseError
               | BadSpecialForm String LispVal
               | NotFunction String String
               | UnboundVar String String
               | Default String

showError :: LispError -> String
showError (UnboundVar message varname)  = message ++ ": " ++ varname
showError (BadSpecialForm message form) = message ++ ": " ++ show form
showError (NotFunction message func)    = message ++ ": " ++ show func
showError (Parser parseErr)             = "Parse error at " ++ show parseErr
showError (Default message)             = message
showError (TypeMismatch expected found) = "Invalid type: expected " ++ expected
                                       ++ ", found " ++ show found
showError (NumArgs expected found)      = "Expected " ++ show expected ++
                                          " args; found values " ++ unwordsList found
  where unwordsList = unwords . map show

instance Show LispError where show = showError

instance Error LispError where
  noMsg  = Default "An error has occured"
  strMsg = Default

type ThrowsError = Either LispError

trapError :: (MonadError e m, Show e) => m String -> m String
trapError action = catchError action (return . show)

extractValue :: ThrowsError a -> a
extractValue (Right val) = val

-- LispError raised to the IO Monad
type IOThrowsError = ErrorT LispError IO

liftThrows :: ThrowsError a -> IOThrowsError a
liftThrows (Left err)  = throwError err
liftThrows (Right val) = return val

runIOThrows :: IOThrowsError String -> IO String
runIOThrows action = liftM extractValue (runErrorT $ trapError action)

-- Type to hold REPL environment
type Env = IORef [(String, IORef LispVal)]

nullEnv :: IO Env
nullEnv = newIORef []

isBound :: Env -> String -> IO Bool
isBound envRef var = liftM (isJust . lookup var) (readIORef envRef)

max_stars_repo_name: kimsnj/writeAscheme

max_issues_repo_name: kimsnj/writeAscheme

max_forks_repo_name: kimsnj/writeAscheme


text:
import qualified Data.Vector as U
import Data.Complex

main = print . U.sum $ U.replicate (1000000000 :: Int) (1 :+ 1 ::Complex Double)

max_stars_repo_name: alexkalderimis/vector

max_issues_repo_name: alexkalderimis/vector

max_forks_repo_name: alexkalderimis/vector


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

{-|
Module      : Data.Picture
Description : picture manipulation functions
Copyright   : (c) Mahdi Dibaiee, 2016
License     : GPL-3
Maintainer  : mdibaiee@aol.com
Stability   : experimental
Portability : POSIX
-}
module Data.Picture ( Picture
                    -- * Manipulation functions
                    , grayscale
                    , fade
                    , rotate
                    , contrast
                    , brightness
                    , gamma
                    , invert
                    , compress
                    , embed
                    , resize
                    , Data.Picture.scale
                    -- * Converting between Image and Picture
                    , fromImage
                    , toImage
                    -- * IO operations
                    , readPicture
                    , writePicturePng
                    )
  where

    import Codec.Picture
    import Numeric.LinearAlgebra
    import qualified Data.Vector.Storable as V
    import System.IO
    import Data.Maybe
    import Debug.Trace
    import Data.List (zipWith4)

    -- | (R, G, B, A) color channels
    type Picture = (Matrix Double, Matrix Double, Matrix Double, Matrix Double)

    -- | Converts a JuicyPixel 'Image PixelRGBA8' to 'Picture'
    fromImage :: Image PixelRGBA8 -> Picture
    fromImage Image { imageWidth = w, imageHeight = h, imageData = vec } =
     let [r, g, b, a] = map (reshape w . V.fromList . reverse) (snd $ V.foldl' gp (0, [[],[],[],[]]) (V.map fromIntegral vec))
     in (r, g, b, a)
     where
       gp acc x =
         case acc of
           (0, [r, g, b, a]) -> (1, [x:r, g, b, a])
           (1, [r, g, b, a]) -> (2, [r, x:g, b, a])
           (2, [r, g, b, a]) -> (3, [r, g, x:b, a])
           (3, [r, g, b, a]) -> (0, [r, g, b, x:a])

    -- | Converts a 'Picture' to JuicyPixel 'Image PixelRGBA8'
    toImage :: Picture -> Image PixelRGBA8
    toImage (r, g, b, a) =
      let (fr, fg, fb, fa) = (toList $ flatten r, toList $ flatten g, toList $ flatten b, toList $ flatten a)
          img = V.map (fromIntegral . floor) . V.concat $ zipWith4 (\a b c d -> vector [a, b, c, d]) fr fg fb fa
      in Image { imageWidth = cols r, imageHeight = rows r, imageData = img }

    -- | Reads a 'Picture' from specified path
    readPicture :: FilePath -> IO (Either String Picture)
    readPicture path = do
      img <- readImage path
      return $ case img of
        Left err -> Left err
        Right im -> Right $ fromImage (convertRGBA8 im)

    -- | Write the specified 'Picture' to a PNG file
    writePicturePng :: FilePath -> Picture -> IO ()
    writePicturePng path pic = writePng path (toImage pic)

    -- | Turn the 'Picture' grayscale
    grayscale :: Picture -> Picture
    grayscale (r, g, b, a) =
      let (fr, fg, fb) = (flatten r, flatten g, flatten b)
          mean = reshape (cols r) $ V.map (/ 3) (fr + fg + fb)
      in (mean, mean, mean, a)

    -- | Fade the 'Picture' by a number between 0 and 1
    fade :: Double -> Picture -> Picture
    fade opacity (r, g, b, a) = (r, g, b, f a)
      where
        f = cmap (*opacity)

    -- | Set contrast level of 'Picture', a number between -255 and 255
    contrast :: Double -> Picture -> Picture
    contrast level (r, g, b, a) = (f r, f g, f b, a)
      where
        cfactor = (259 * (level + 255)) / (255 * (259 - level))
        f = cmap (\x -> pixelBound $ cfactor * (x - 128) + 128)

    -- | Set brightness level of 'Picture', a number between -255 and 255
    brightness :: Double -> Picture -> Picture
    brightness level (r, g, b, a) = (f r, f g, f b, a)
      where
        f = cmap (pixelBound . (+level))

    -- | Set gamma level of 'Picture'
    gamma :: Int -> Picture -> Picture
    gamma level (r, g, b, a) = (f r, f g, f b, a)
      where
        f = cmap (\x -> pixelBound $ 255 * (x / 255) ^ level)

    -- | Inverts the 'Picture'
    invert :: Picture -> Picture
    invert (r, g, b, a) = (f r, f g, f b, a)
      where
        f = cmap (`subtract` 255)

    -- | Rotate 'Picture' for the specified degrees, around the specified origin.
    -- If the origin is `Nothing`, rotates around the center
    rotate :: Double -> Maybe (Int, Int) -> Picture -> Picture
    rotate deg orig (r, g, b, a) = (f r, f g, f b, f a)
      where
        -- rotation in radians
        rad = deg * pi / 180

        -- rotation matrix (clockwise)
        rm = fromLists [[cos rad, sin rad],
                        [negate $ sin rad, cos rad]]

        -- origin of rotation
        (originX, originY) = if isJust orig then fromJust orig else (cols r `div` 2, rows r `div` 2)

        -- all index pairs
        indices = [vector [fromIntegral x - fromIntegral originX, fromIntegral y - fromIntegral originY] | y <- [0..rows r - 1], x <- [0..cols r - 1]]

        -- rotate them using rotation matrix
        rotatedIndices :: [[Int]] = map (\r -> toList . V.map (fromIntegral . round) $ rm #> r) indices

        -- move them back to the origin
        movedIndices = map (\[x, y] -> [x + originX, y + originY]) rotatedIndices

        f m = reshape (cols m) $ fromList $ map (\[x, y] -> if y < 0 || y >= rows r || x < 0 || x >= cols r then 0 else m `atIndex` (y, x)) movedIndices

    -- | Compress the image using SVD
    -- note: this is not size compression, it's just a k-rank approximation of the image
    compress :: Int -> Picture -> Picture
    compress rate (r, g, b, a) = (f r, f g, f b, a)
      where
        k = cols r - rate
        f m =
          let (u, s, v) = svd m
              si = diagRect 0 s (rows m) (cols m)
              (mu, ms, mv) = (u ?? (All, Take k), si ?? (Take k, Take k), (tr v) ?? (Take k, All))
          in mu <> ms <> mv

    -- | Embed a 'Picture' into another one, in the specified position-}
    embed :: Picture -> (Int, Int) -> Picture -> Picture
    embed (br, bg, bb, ba) (x, y) (lr, lg, lb, la) = (f br lmr, f bg lmg, f bb lmb, maxAlpha)
      where
        (lmr, lmg, lmb, lma) = (fit lr, fit lg, fit lb, fit la)
        scaledAlpha = cmap (/255) lma

        fit m =
          let distance = y * cols br
              total = rows br * cols br
              xPush = (rows m><x) (repeat 0)
              xPast = (rows m><(cols br - x - cols m)) (repeat 0)
              positioned = xPush ||| m ||| xPast
              flat = toList $ flatten positioned
          in (rows br><cols br) $ replicate distance 0 ++ flat ++ repeat 0

        f b lm = (b * (cmap (1-) scaledAlpha)) + (lm * scaledAlpha)
        maxAlpha = (rows ba><cols ba) $ zipWith max (toList . flatten $ ba) (toList . flatten . fit $ la)

    -- | Resize an image using nearest-neighbor interpolation
    resize :: (Int, Int) -> Picture -> Picture
    resize (sWidth, sHeight) (r, g, b, a) = (f r, f g, f b, f a)
      where
        initial = vector [0..fromIntegral sWidth * fromIntegral sHeight - 1]
        (width, height) = (cols r, rows r)
        factor = 2 ^ 16
        (xRatio, yRatio) = (width * factor `div` sWidth + 1, height * factor `div` sHeight + 1)
        f m = reshape sWidth $ V.map replace initial
          where
            v = flatten m
            replace index =
              let (x, y) = (floor index `mod` sWidth, floor index `div` sWidth)
                  (px, py) = (x * xRatio `div` factor, y * yRatio `div` factor)
              in v ! (py * width + px)

    -- | Scale an image using the resize function
    scale :: Double -> Picture -> Picture
    scale 1 p = p
    scale s (r, g, b, a) = resize (floor $ s * width, floor $ s * height) (r, g, b, a)
      where
        (width, height) = (fromIntegral $ cols r, fromIntegral $ rows r)

    bound (l, u) x = max l $ min u x
    pixelBound = bound (0, 255)

max_stars_repo_name: mdibaiee/picedit

max_issues_repo_name: mdibaiee/picedit

max_forks_repo_name: mdibaiee/picedit


text:
module NeuralA where

import Numeric.LinearAlgebra

x = vector [-5.0,-4.9..5.0]
y1 = cmap sigmoid x

sigmoid :: Double -> Double
sigmoid x = 1.0 / (1.0 + exp(-x))

step_func :: Double -> Double
step_func x
  | x > 0.0 = 1.0
  | otherwise = 0.0

relu :: Double -> Double
relu x
  | x > 0.0 = x
  | otherwise = 0.0

max_stars_repo_name: chupaaaaaaan/nn-with-haskell

max_issues_repo_name: chupaaaaaaan/nn-with-haskell

max_forks_repo_name: chupaaaaaaan/nn-with-haskell


text:
{-# LANGUAGE AllowAmbiguousTypes                      #-}
{-# LANGUAGE ConstraintKinds                          #-}
{-# LANGUAGE DataKinds                                #-}
{-# LANGUAGE DefaultSignatures                        #-}
{-# LANGUAGE FlexibleContexts                         #-}
{-# LANGUAGE FlexibleInstances                        #-}
{-# LANGUAGE GADTs                                    #-}
{-# LANGUAGE RankNTypes                               #-}
{-# LANGUAGE ScopedTypeVariables                      #-}
{-# LANGUAGE TypeApplications                         #-}
{-# LANGUAGE TypeFamilies                             #-}
{-# LANGUAGE TypeOperators                            #-}
{-# LANGUAGE UndecidableInstances                     #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}

module Backprop.Learn.Initialize (
    Initialize(..)
  , gInitialize
  , initializeNormal
  , initializeSingle
  -- * Reshape
  , reshapeR
  , reshapeLRows
  , reshapeLCols
  ) where

import           Control.Monad.Primitive
import           Data.Complex
import           Data.Proxy
import           Data.Type.Equality
import           Data.Type.Tuple
import           Data.Vinyl
import           GHC.TypeLits.Compare
import           GHC.TypeNats
import           Generics.OneLiner
import           Numeric.LinearAlgebra.Static.Vector
import           Statistics.Distribution
import           Statistics.Distribution.Normal
import qualified Data.Vector.Generic                 as VG
import qualified Data.Vector.Generic.Sized           as SVG
import qualified Numeric.LinearAlgebra.Static        as H
import qualified System.Random.MWC                   as MWC

-- | Class for types that are basically a bunch of 'Double's, which can be
-- initialized with a given identical and independent distribution.
class Initialize p where
    initialize
        :: (ContGen d, PrimMonad m)
        => d
        -> MWC.Gen (PrimState m)
        -> m p

    default initialize
        :: (ADTRecord p, Constraints p Initialize, ContGen d, PrimMonad m)
        => d
        -> MWC.Gen (PrimState m)
        -> m p
    initialize = gInitialize

-- | 'initialize' for any instance of 'Generic'.
gInitialize
    :: (ADTRecord p, Constraints p Initialize, ContGen d, PrimMonad m)
    => d
    -> MWC.Gen (PrimState m)
    -> m p
gInitialize d g = createA' @Initialize (initialize d g)

-- | Helper over 'inititialize' for a gaussian distribution centered around
-- zero.
initializeNormal
    :: (Initialize p, PrimMonad m)
    => Double                               -- ^ standard deviation
    -> MWC.Gen (PrimState m)
    -> m p
initializeNormal = initialize . normalDistr 0

-- | 'initialize' definition if @p@ is a single number.
initializeSingle
    :: (ContGen d, PrimMonad m, Fractional p)
    => d
    -> MWC.Gen (PrimState m)
    -> m p
initializeSingle d = fmap realToFrac . genContVar d

instance Initialize Double where
    initialize = initializeSingle
instance Initialize Float where
    initialize = initializeSingle

-- | Initializes real and imaginary components identically
instance Initialize a => Initialize (Complex a) where

instance Initialize T0
instance Initialize a => Initialize (TF a)
instance (Initialize a, Initialize b) => Initialize (a :# b)

instance RPureConstrained Initialize as => Initialize (T as) where
    initialize d g = rtraverse (fmap TF)
                   $ rpureConstrained @Initialize (initialize d g)

-- instance (Initialize a, ListC (Initialize <$> as), Known Length as) => Initialize (NETup (a ':| as)) where
--     initialize d g = NET <$> initialize d g
--                          <*> initialize d g

instance Initialize ()
instance (Initialize a, Initialize b) => Initialize (a, b)
instance (Initialize a, Initialize b, Initialize c) => Initialize (a, b, c)
instance (Initialize a, Initialize b, Initialize c, Initialize d) => Initialize (a, b, c, d)
instance (Initialize a, Initialize b, Initialize c, Initialize d, Initialize e) => Initialize (a, b, c, d, e)

instance (VG.Vector v a, KnownNat n, Initialize a) => Initialize (SVG.Vector v n a) where
    initialize d = SVG.replicateM . initialize d

instance KnownNat n => Initialize (H.R n) where
    initialize d = fmap vecR . initialize d
instance KnownNat n => Initialize (H.C n) where
    initialize d = fmap vecC . initialize d

instance (KnownNat n, KnownNat m) => Initialize (H.L n m) where
    initialize d = fmap vecL . initialize d
instance (KnownNat n, KnownNat m) => Initialize (H.M n m) where
    initialize d = fmap vecM . initialize d

-- | Reshape a vector to have a different amount of items  If the matrix is
-- grown, new weights are initialized according to the given distribution.
reshapeR
    :: forall i j d m. (ContGen d, PrimMonad m, KnownNat i, KnownNat j)
    => d
    -> MWC.Gen (PrimState m)
    -> H.R i
    -> m (H.R j)
reshapeR d g x = case Proxy @j %<=? Proxy @i of
    LE  Refl      -> pure . vecR . SVG.take @_ @j @(i - j) . rVec $ x
    NLE Refl Refl -> (x H.#) <$> initialize @(H.R (j - i)) d g

-- | Reshape a matrix to have a different amount of rows  If the matrix
-- is grown, new weights are initialized according to the given
-- distribution.
reshapeLRows
    :: forall i j n d m. (ContGen d, PrimMonad m, KnownNat n, KnownNat i, KnownNat j)
    => d
    -> MWC.Gen (PrimState m)
    -> H.L i n
    -> m (H.L j n)
reshapeLRows d g x = case Proxy @j %<=? Proxy @i of
    LE Refl       -> pure . rowsL . SVG.take @_ @j @(i - j) . lRows $ x
    NLE Refl Refl -> (x H.===) <$> initialize @(H.L (j - i) n) d g

-- | Reshape a matrix to have a different amount of columns.  If the matrix
-- is grown, new weights are initialized according to the given
-- distribution.
reshapeLCols
    :: forall i j n d m. (ContGen d, PrimMonad m, KnownNat n, KnownNat i, KnownNat j)
    => d
    -> MWC.Gen (PrimState m)
    -> H.L n i
    -> m (H.L n j)
reshapeLCols d g x = case Proxy @j %<=? Proxy @i of
    LE Refl       -> pure . colsL . SVG.take @_ @j @(i - j) . lCols $ x
    NLE Refl Refl -> (x H.|||) <$> initialize @(H.L n (j - i)) d g

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE DeriveTraversable #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE RankNTypes #-}
#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 702
{-# LANGUAGE Trustworthy #-}
{-# LANGUAGE DeriveGeneric #-}
#endif
{-# LANGUAGE DeriveDataTypeable #-}
#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ < 708
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE ScopedTypeVariables #-}
#endif
-----------------------------------------------------------------------------
-- |
-- License     :  BSD-style (see the file LICENSE)
-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
-- Stability   :  provisional
-- Portability :  portable
--
-- Operations on affine spaces.
-----------------------------------------------------------------------------
module Linear.Affine where

import Control.Applicative
import Control.Lens
import Data.Complex (Complex)
import Data.Data
import Data.Distributive
import Data.Foldable as Foldable
import Data.Functor.Bind
import Data.Functor.Rep as Rep
import Data.HashMap.Lazy (HashMap)
import Data.Hashable
import Data.IntMap (IntMap)
import Data.Ix
import Data.Map (Map)
import Data.Vector (Vector)
import Foreign.Storable
#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 702
import GHC.Generics (Generic)
#endif
#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 706
import GHC.Generics (Generic1)
#endif
import Linear.Epsilon
import Linear.Metric
import Linear.Plucker
import Linear.Quaternion
import Linear.V
import Linear.V0
import Linear.V1
import Linear.V2
import Linear.V3
import Linear.V4
import Linear.Vector

#ifdef HLINT
{-# ANN module "HLint: ignore Unused LANGUAGE pragma" #-}
#endif

-- | An affine space is roughly a vector space in which we have
-- forgotten or at least pretend to have forgotten the origin.
--
-- > a .+^ (b .-. a)  =  b@
-- > (a .+^ u) .+^ v  =  a .+^ (u ^+^ v)@
-- > (a .-. b) ^+^ v  =  (a .+^ v) .-. q@
class Additive (Diff p) => Affine p where
  type Diff p :: * -> *

  infixl 6 .-.
  -- | Get the difference between two points as a vector offset.
  (.-.) :: Num a => p a -> p a -> Diff p a

  infixl 6 .+^
  -- | Add a vector offset to a point.
  (.+^) :: Num a => p a -> Diff p a -> p a

  infixl 6 .-^
  -- | Subtract a vector offset from a point.
  (.-^) :: Num a => p a -> Diff p a -> p a
  p .-^ v = p .+^ negated v
  {-# INLINE (.-^) #-}

-- | Compute the quadrance of the difference (the square of the distance)
qdA :: (Affine p, Foldable (Diff p), Num a) => p a -> p a -> a
qdA a b = Foldable.sum (fmap (join (*)) (a .-. b))
{-# INLINE qdA #-}

-- | Distance between two points in an affine space
distanceA :: (Floating a, Foldable (Diff p), Affine p) => p a -> p a -> a
distanceA a b = sqrt (qdA a b)
{-# INLINE distanceA #-}

#define ADDITIVEC(CTX,T) instance CTX => Affine T where type Diff T = T ; \
  (.-.) = (^-^) ; {-# INLINE (.-.) #-} ; (.+^) = (^+^) ; {-# INLINE (.+^) #-} ; \
  (.-^) = (^-^) ; {-# INLINE (.-^) #-}
#define ADDITIVE(T) ADDITIVEC((), T)

ADDITIVE([])
ADDITIVE(Complex)
ADDITIVE(ZipList)
ADDITIVE(Maybe)
ADDITIVE(IntMap)
ADDITIVE(Identity)
ADDITIVE(Vector)
ADDITIVE(V0)
ADDITIVE(V1)
ADDITIVE(V2)
ADDITIVE(V3)
ADDITIVE(V4)
ADDITIVE(Plucker)
ADDITIVE(Quaternion)
ADDITIVE(((->) b))
ADDITIVEC(Ord k, (Map k))
ADDITIVEC((Eq k, Hashable k), (HashMap k))
ADDITIVEC(Dim n, (V n))

-- | A handy wrapper to help distinguish points from vectors at the
-- type level
newtype Point f a = P (f a)
  deriving ( Eq, Ord, Show, Read, Monad, Functor, Applicative, Foldable
           , Traversable, Apply, Additive, Metric
           , Fractional , Num, Ix, Storable, Epsilon
           , Hashable
#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 702
           , Generic
#endif
#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 706
           , Generic1
#endif
#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 708
           , Typeable, Data
#endif
           )

#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ < 708
instance forall f. Typeable1 f => Typeable1 (Point f) where
  typeOf1 _ = mkTyConApp (mkTyCon3 "linear" "Linear.Affine" "Point") [] `mkAppTy`
              typeOf1 (undefined :: f a)

deriving instance (Data (f a), Typeable1 f, Typeable a) => Data (Point f a)
#endif

lensP :: Lens' (Point g a) (g a)
lensP afb (P a) = P <$> afb a
{-# INLINE lensP #-}

_Point :: Iso' (Point f a) (f a)
_Point = iso (\(P a) -> a) P
{-# INLINE _Point #-}

instance (t ~ Point g b) => Rewrapped (Point f a) t
instance Wrapped (Point f a) where
  type Unwrapped (Point f a) = f a
  _Wrapped' = _Point
  {-# INLINE _Wrapped' #-}

instance Bind f => Bind (Point f) where
  join (P m) = P $ join $ fmap (\(P m')->m') m

instance Distributive f => Distributive (Point f) where
  distribute = P . collect (\(P p) -> p)

instance Representable f => Representable (Point f) where
  type Rep (Point f) = Rep f
  tabulate f = P (tabulate f)
  {-# INLINE tabulate #-}
  index (P xs) = Rep.index xs
  {-# INLINE index #-}

type instance Index (Point f a) = Index (f a)
type instance IxValue (Point f a) = IxValue (f a)

instance Ixed (f a) => Ixed (Point f a) where
  ix l = lensP . ix l
  {-# INLINE ix #-}

instance Traversable f => Each (Point f a) (Point f b) a b where
  each = traverse
  {-# INLINE each #-}

instance R1 f => R1 (Point f) where
  _x = lensP . _x
  {-# INLINE _x #-}

instance R2 f => R2 (Point f) where
  _y = lensP . _y
  {-# INLINE _y #-}
  _xy = lensP . _xy
  {-# INLINE _xy #-}

instance R3 f => R3 (Point f) where
  _z = lensP . _z
  {-# INLINE _z #-}
  _xyz = lensP . _xyz
  {-# INLINE _xyz #-}

instance R4 f => R4 (Point f) where
  _w = lensP . _w
  {-# INLINE _w #-}
  _xyzw = lensP . _xyzw
  {-# INLINE _xyzw #-}

instance Additive f => Affine (Point f) where
  type Diff (Point f) = f
  P x .-. P y = x ^-^ y
  {-# INLINE (.-.) #-}
  P x .+^ v = P (x ^+^ v)
  {-# INLINE (.+^) #-}
  P x .-^ v = P (x ^-^ v)
  {-# INLINE (.-^) #-}

-- | Vector spaces have origins.
origin :: (Additive f, Num a) => Point f a
origin = P zero

-- | An isomorphism between points and vectors, given a reference
--   point.
relative :: (Additive f, Num a) => Point f a -> Iso' (Point f a) (f a)
relative p0 = iso (.-. p0) (p0 .+^)
{-# INLINE relative #-}

max_stars_repo_name: abbradar/linear

max_issues_repo_name: abbradar/linear

max_forks_repo_name: abbradar/linear


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

module Language.Hakaru.Distribution where

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

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

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

bern :: Double -> Dist Bool
bern p = Dist {logDensity = (\ (Discrete x) -> log (if x then p else 1 - p)),
               distSample = (\ g -> do t <- MWC.uniformR (0,1) g
                                       return $ Discrete (t <= p))}

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

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

marsaglia :: (MWC.Variate a, Ord a, Floating a, PrimMonad m) => PRNG m -> m (a, a)
marsaglia g = do -- "Marsaglia polar method"
  x <- MWC.uniformR (-1,1) g
  y <- MWC.uniformR (-1,1) g
  let s = x * x + y * y
      q = sqrt ((-2) * log s / s)
  if 1 >= s && s > 0 then return (x * q, y * q) else marsaglia g

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

chooseIndex :: (PrimMonad m) => [Double] -> PRNG m -> m Int
chooseIndex probs g = do
  p <- MWC.uniform g
  return $ fromMaybe (error ("chooseIndex: failure p=" ++ show p))
           (findIndex (p <=) (scanl1 (+) probs))

normal_rng :: (Real a, Floating a, MWC.Variate a, PrimMonad m) =>
              a -> a -> PRNG m -> m a
normal_rng mu sd g | sd > 0 = do (x, _) <- marsaglia g
                                 return (mu + sd * x)
normal_rng _ _ _ = error "normal: invalid parameters"

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

normal :: Double -> Double -> Dist Double
normal mu sd = Dist {logDensity = normalLogDensity mu sd . fromLebesgue,
                     distSample = (\g -> liftM Lebesgue $ normal_rng mu sd g)}

categoricalLogDensity :: (Eq b, Floating a) => [(b, a)] -> b -> a
categoricalLogDensity list x = log $ fromMaybe 0 (lookup x list)

categoricalSample :: (Num b, Ord b, PrimMonad m, MWC.Variate b) =>
    [(t,b)] -> PRNG m -> m t
categoricalSample list g = do
  let total = sum $ map snd list
  p <- MWC.uniformR (0, total) g
  let sumList = scanl1 (\acc (a, b) -> (a, b + snd(acc))) list
      elem' = fst $ head $ filter (\(_,p0) -> p <= p0) sumList
  return elem'

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

lnFact :: Int -> Double
lnFact = logFactorial

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

         make_poisson :: (PrimMonad m) => PRNG m -> m Int
         make_poisson g = do u <- MWC.uniformR (-0.5,0.5) g
                             v <- MWC.uniformR (0,1) g
                             let us = 0.5 - abs u
                                 k = floor $ (2*a / us + b)*u + lambda + 0.43
                             case () of
                               () | us >= 0.07 && v <= vr -> return k
                               () | k < 0 -> make_poisson g
                               () | us <= 0.013 && v > us -> make_poisson g
                               () | accept_region us v k -> return k
                               _  -> make_poisson g

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

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

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

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

gamma :: Double -> Double -> Dist Double
gamma shape scl = Dist {logDensity = gammaLogDensity shape scl . fromLebesgue,
                        distSample = (\g -> liftM Lebesgue $ gamma_rng shape scl g)}

beta_rng :: (PrimMonad m) => Double -> Double -> PRNG m -> m Double
beta_rng a b g | a <= 1.0 && b <= 1.0 = do
                 u <- MWC.uniformR (0.0, 1.0) g
                 v <- MWC.uniformR (0.0, 1.0) g
                 let x = u ** (recip a)
                     y = v ** (recip b)
                 case (x+y) <= 1.0 of
                   True -> return $ x / (x + y)
                   False -> beta_rng a b g
beta_rng a b g = do ga <- gamma_rng a 1 g
                    gb <- gamma_rng b 1 g
                    return $ ga / (ga + gb)

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

beta :: Double -> Double -> Dist Double
beta a b = Dist {logDensity = betaLogDensity a b . fromLebesgue,
                 distSample = (\g -> liftM Lebesgue $ beta_rng a b g)}

laplace_rng :: (PrimMonad m) => Double -> Double -> PRNG m -> m Double
laplace_rng mu sd g = MWC.uniformR (0.0, 1.0) g >>= sample
   where sample u = return $ case u < 0.5 of
                               True  -> mu + sd * log (u + u)
                               False -> mu - sd * log (2.0 - u - u)

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

laplace :: Double -> Double -> Dist Double
laplace mu sd = Dist {logDensity = laplaceLogDensity mu sd . fromLebesgue,
                      distSample = (\g -> liftM Lebesgue $ laplace_rng mu sd g)}

-- Consider having dirichlet return Vector
-- Note: This is actually symmetric dirichlet
dirichlet_rng :: (PrimMonad m) => Int ->  Double -> PRNG m -> m [Double]
dirichlet_rng n' a g' = liftM normalize $ gammas g' n'
  where gammas _ 0 = return ([], 0)
        gammas g n = do (xs, total) <- gammas g (n-1)
                        x <- gamma_rng a 1 g
                        return ((x : xs), x+total)
        normalize (b, total) = map (/ total) b

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

max_stars_repo_name: SamuelSchlesinger/hakaru

max_issues_repo_name: SamuelSchlesinger/hakaru

max_forks_repo_name: SamuelSchlesinger/hakaru


text:
{-#LANGUAGE BangPatterns #-}

module Tests.Correlation
    ( tests ) where

import Test.Framework
import Test.Framework.Providers.QuickCheck2
import Test.Framework.Providers.HUnit
import Test.HUnit (Assertion, (@=?))
import qualified Data.Vector as V
import Statistics.Correlation.Kendall

tests :: Test
tests = testGroup "Correlation"
    [ testProperty "Kendall test -- general" testKendall
    , testCase "Kendall test -- special cases" testKendallSpecial
    ]

testKendall :: [(Double, Double)] -> Bool
testKendall xy | isNaN r1 = isNaN r2
               | otherwise = r1 == r2
  where
    r1 = kendallBruteForce xy
    r2 = kendall $ V.fromList xy

testKendallSpecial :: Assertion
testKendallSpecial = ys @=? map (kendall.V.fromList) xs
  where
    (xs, ys) = unzip testData
    testData :: [([(Double, Double)], Double)]
    testData = [ ( [(1,1), (2,2), (3,1), (1,5), (2,2)], -0.375 )
               , ( [(1,3), (1,3), (1,3), (3,2), (3,5)], 0)
               ]

kendallBruteForce :: [(Double, Double)] -> Double
kendallBruteForce xy = (n_c - n_d) / sqrt ((n_0 - n_1) * (n_0 - n_2))
  where
    allPairs = f xy
    (n_c, n_d, n_1, n_2) = foldl g (0,0,0,0) allPairs
    n_0 = fromIntegral.length $ allPairs
    g (!nc, !nd, !n1, !n2) ((x1, y1), (x2, y2))
      | (x2 - x1) * (y2 - y1) > 0 = (nc+1, nd, n1, n2)
      | (x2 - x1) * (y2 - y1) < 0 = (nc, nd+1, n1, n2)
      | otherwise = if x1 == x2
                       then if y1 == y2
                               then (nc, nd, n1+1, n2+1)
                               else (nc, nd, n1+1, n2)
                       else (nc, nd, n1, n2+1)
    f (x:xs) = zip (repeat x) xs ++ f xs
    f _ = []

max_stars_repo_name: StefanHubner/statistics

max_issues_repo_name: StefanHubner/statistics

max_forks_repo_name: StefanHubner/statistics