text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE NoStarIsType #-}
{-# OPTIONS_GHC -Wno-redundant-constraints #-}
{-# OPTIONS_GHC -fno-warn-incomplete-uni-patterns #-}

-- | Arrays with a fixed shape, with a HMatrix representation.
module NumHask.Array.HMatrix
  ( -- $usage
    Array (..),

    -- * Representation
    --
    -- With no functor instance, we instead supply the representable API
    index,
    tabulate,

    -- * Conversion
    shape,
    toDynamic,
    toFixed,
    fromFixed,

    -- * Operators
    reshape,
    transpose,
    diag,
    ident,
    singleton,
    selects,
    selectsExcept,
    folds,
    concatenate,
    insert,
    append,
    reorder,
    expand,
    slice,
    squeeze,

    -- * Scalar
    --
    -- Scalar specialisations
    Scalar,
    fromScalar,
    toScalar,

    -- * Vector
    --
    -- Vector specialisations.
    Vector,

    -- * Matrix
    --
    -- Matrix specialisations.
    Matrix,
    col,
    row,
    safeCol,
    safeRow,
    mmult,
  )
where

import qualified Data.Vector as V
import GHC.Exts (IsList (..))
import GHC.TypeLits
import qualified NumHask.Array.Dynamic as D
import qualified NumHask.Array.Fixed as F
import NumHask.Array.Shape
import NumHask.Prelude as P
import qualified Numeric.LinearAlgebra as H
import qualified Numeric.LinearAlgebra.Devel as H
import qualified Prelude
import Data.Proxy (Proxy(..))

-- $setup
-- >>> :set -XDataKinds
-- >>> :set -XOverloadedLists
-- >>> :set -XTypeFamilies
-- >>> :set -XFlexibleContexts
-- >>> import GHC.TypeLits
-- >>> import GHC.Int
-- >>> let s = [1] :: Array ('[] :: [Nat]) Int64 -- scalar
-- >>> let v = [1,2,3] :: Array '[3] Int64       -- vector
-- >>> let m = [0..11] :: Array '[3,4] Int64     -- matrix
-- >>> let a = [1..24] :: Array '[2,3,4] Int64

-- | a multidimensional array with a type-level shape
--
-- >>> let a = [1..24] :: Array '[2,3,4] Int64
-- >>> a
-- [[[1, 2, 3, 4],
--   [5, 6, 7, 8],
--   [9, 10, 11, 12]],
--  [[13, 14, 15, 16],
--   [17, 18, 19, 20],
--   [21, 22, 23, 24]]]
--
-- >>> [1,2,3] :: Array '[2,2] Int64
-- [[*** Exception: NumHaskException {errorMessage = "shape mismatch"}
newtype Array s a = Array {unArray :: H.Matrix a}
  deriving (Show, Generic)

instance
  ( Additive a,
    HasShape s,
    H.Container H.Vector a,
    Num a
  ) =>
  Additive (Array s a)
  where
  (+) (Array x1) (Array x2) = Array $ H.add x1 x2

  zero = Array $ H.konst zero (n, m)
    where
      s = shapeVal (toShape @s)
      [n, m] = s

instance
  ( Multiplicative a,
    HasShape s,
    H.Container H.Vector a,
    Num (H.Vector a),
    Num a
  ) =>
  Multiplicative (Array s a)
  where
  (*) (Array x1) (Array x2) = Array $ H.liftMatrix2 (Prelude.*) x1 x2

  one = Array $ H.konst one (n, m)
    where
      s = shapeVal (toShape @s)
      [n, m] = s

-- (<.>) (Array a) (Array b) = H.sumElements $ H.liftMatrix2 (Prelude.*) a b

instance
  ( Subtractive a,
    HasShape s,
    H.Container H.Vector a,
    Num a
  ) =>
  Subtractive (Array s a)
  where
  negate (Array x1) = Array $ H.cmap negate x1

instance
  (HasShape s, Multiplicative a, H.Container H.Vector a, Num a) =>
  MultiplicativeAction (Array s a) a
  where
  (.*) s (Array r) = Array $ H.cmap (s *) r
  {-# INLINE (.*) #-}

  (*.) (Array r) s = Array $ H.cmap (*s) r
  {-# INLINE (*.) #-}

instance
  (HasShape s, Additive a, H.Container H.Vector a, Num a) =>
  AdditiveAction (Array s a) a
  where
  (.+) s (Array r) = Array $ H.cmap (s +) r
  {-# INLINE (.+) #-}

  (+.) (Array r) s = Array $ H.cmap (+s) r
  {-# INLINE (+.) #-}

instance
  (HasShape s, Subtractive a, H.Container H.Vector a, Num a) =>
  SubtractiveAction (Array s a) a
  where
  (.-) s (Array r) = Array $ H.cmap (s -) r
  {-# INLINE (.-) #-}

  (-.) (Array r) s = Array $ H.cmap (\x -> x - s) r
  {-# INLINE (-.) #-}

instance
  (HasShape s, Divisive a, H.Container H.Vector a, Num a) =>
  DivisiveAction (Array s a) a
  where
  (./) s (Array r) = Array $ H.cmap (s /) r
  {-# INLINE (./) #-}

  (/.) (Array r) s = Array $ H.cmap (/ s) r
  {-# INLINE (/.) #-}

-- | from flat list
instance
  ( HasShape s,
    H.Element a
  ) =>
  IsList (Array s a)
  where
  type Item (Array s a) = a

  fromList l =
    bool
      (throw (NumHaskException "shape mismatch"))
      (Array $ H.reshape n $ H.fromList l)
      ((length l == 1 && null s) || (length l == size s))
    where
      s = shapeVal (toShape @s)
      n = Prelude.last s

  toList (Array v) = H.toList $ H.flatten v

-- | Get shape of an Array as a value.
--
-- >>> shape a
-- [2,3,4]
shape :: forall a s. (HasShape s) => Array s a -> [Int]
shape _ = shapeVal $ toShape @s
{-# INLINE shape #-}

-- | Convert to a dynamic array.
toDynamic :: (HasShape s, H.Element a) => Array s a -> D.Array a
toDynamic a@(Array h) = D.fromFlatList (shape a) (mconcat $ H.toLists h)

-- | Convert to a fixed array.
toFixed :: (HasShape s, H.Element a) => Array s a -> F.Array s a
toFixed (Array h) = fromList (mconcat $ H.toLists h)

-- | Convert from a fixed array.
fromFixed :: (HasShape s, H.Element a) => F.Array s a -> Array s a
fromFixed a = fromList (P.toList a)

-- | with no fmap, we supply the representable API
index ::
  forall s a.
  ( HasShape s,
    H.Element a,
    H.Container H.Vector a
  ) =>
  Array s a ->
  [Int] ->
  a
index (Array v) i = H.flatten v `H.atIndex` flatten s i
  where
    s = shapeVal (toShape @s)

-- | tabulate an array with a generating function
--
-- >>> tabulate [2,3,4] ((1+) . flatten [2,3,4]) == a
-- True
tabulate ::
  forall s a.
  ( HasShape s,
    H.Element a
  ) =>
  ([Int] -> a) ->
  Array s a
tabulate f =
  fromList (V.toList $ V.generate (size s) (f . shapen s))
  where
    s = shapeVal (toShape @s)

-- | Reshape an array (with the same number of elements).
--
-- >>> reshape a :: Array '[4,3,2] Int64
-- [[[1, 2],
--   [3, 4],
--   [5, 6]],
--  [[7, 8],
--   [9, 10],
--   [11, 12]],
--  [[13, 14],
--   [15, 16],
--   [17, 18]],
--  [[19, 20],
--   [21, 22],
--   [23, 24]]]
reshape ::
  forall a s s'.
  ( Size s ~ Size s',
    HasShape s,
    HasShape s',
    H.Container H.Vector a
  ) =>
  Array s a ->
  Array s' a
reshape a = tabulate (index a . shapen s . flatten s')
  where
    s = shapeVal (toShape @s)
    s' = shapeVal (toShape @s')

-- | Reverse indices eg transposes the element A/ijk/ to A/kji/.
--
-- >>> index (transpose a) [1,0,0] == index a [0,0,1]
-- True
transpose :: forall a s. (H.Element a, H.Container H.Vector a, HasShape s, HasShape (Reverse s)) => Array s a -> Array (Reverse s) a
transpose a = tabulate (index a . reverse)

-- | The identity array.
--
-- >>> ident :: Array '[3,2] Int64
-- [[1, 0],
--  [0, 1],
--  [0, 0]]
ident :: forall a s. (H.Element a, H.Container H.Vector a, HasShape s, Additive a, Multiplicative a) => Array s a
ident = tabulate (bool zero one . isDiag)
  where
    isDiag [] = True
    isDiag [_] = True
    isDiag [x, y] = x == y
    isDiag (x : y : xs) = x == y && isDiag (y : xs)

-- | Extract the diagonal of an array.
--
-- >>> diag (ident :: Array '[3,2] Int64)
-- [1, 1]
diag ::
  forall a s.
  ( HasShape s,
    HasShape '[Minimum s],
    H.Element a,
    H.Container H.Vector a
  ) =>
  Array s a ->
  Array '[Minimum s] a
diag a = tabulate go
  where
    go [] = throw (NumHaskException "Rank Underflow")
    go (s' : _) = index a (replicate (length ds) s')
    ds = shapeVal (toShape @s)

-- | Create an array composed of a single value.
--
-- >>> singleton one :: Array '[3,2] Int64
-- [[1, 1],
--  [1, 1],
--  [1, 1]]
singleton :: (H.Element a, H.Container H.Vector a, HasShape s) => a -> Array s a
singleton a = tabulate (const a)

-- | Select an array along dimensions.
--
-- >>> let s = selects (Proxy :: Proxy '[0,1]) [1,1] a
-- >>> :t s
-- s :: Array '[4] Int64
--
-- >>> s
-- [17, 18, 19, 20]
selects ::
  forall ds s s' a.
  ( HasShape s,
    HasShape ds,
    HasShape s',
    s' ~ DropIndexes s ds,
    H.Element a,
    H.Container H.Vector a
  ) =>
  Proxy ds ->
  [Int] ->
  Array s a ->
  Array s' a
selects _ i a = tabulate go
  where
    go s = index a (addIndexes s ds i)
    ds = shapeVal (toShape @ds)

-- | Select an index /except/ along specified dimensions.
--
-- >>> let s = selectsExcept (Proxy :: Proxy '[2]) [1,1] a
-- >>> :t s
-- s :: Array '[4] Int64
--
-- >>> s
-- [17, 18, 19, 20]
selectsExcept ::
  forall ds s s' a.
  ( HasShape s,
    HasShape ds,
    HasShape s',
    s' ~ TakeIndexes s ds,
    H.Element a,
    H.Container H.Vector a
  ) =>
  Proxy ds ->
  [Int] ->
  Array s a ->
  Array s' a
selectsExcept _ i a = tabulate go
  where
    go s = index a (addIndexes i ds s)
    ds = shapeVal (toShape @ds)

-- | Fold along specified dimensions.
--
-- >>> folds sum (Proxy :: Proxy '[1]) a
-- [68, 100, 132]
folds ::
  forall ds st si so a b.
  ( HasShape st,
    HasShape ds,
    HasShape si,
    HasShape so,
    si ~ DropIndexes st ds,
    so ~ TakeIndexes st ds,
    H.Element a,
    H.Container H.Vector a,
    H.Element b,
    H.Container H.Vector b
  ) =>
  (Array si a -> b) ->
  Proxy ds ->
  Array st a ->
  Array so b
folds f d a = tabulate go
  where
    go s = f (selects d s a)

-- | Concatenate along a dimension.
--
-- >>> :t concatenate (Proxy :: Proxy 1) a a
-- concatenate (Proxy :: Proxy 1) a a :: Array '[2, 6, 4] Int64
concatenate ::
  forall a s0 s1 d s.
  ( CheckConcatenate d s0 s1 s,
    Concatenate d s0 s1 ~ s,
    HasShape s0,
    HasShape s1,
    HasShape s,
    KnownNat d,
    H.Element a,
    H.Container H.Vector a
  ) =>
  Proxy d ->
  Array s0 a ->
  Array s1 a ->
  Array s a
concatenate _ s0 s1 = tabulate go
  where
    go s =
      bool
        (index s0 s)
        ( index
            s1
            ( addIndex
                (dropIndex s d)
                d
                ((s !! d) - (ds0 !! d))
            )
        )
        ((s !! d) >= (ds0 !! d))
    ds0 = shapeVal (toShape @s0)
    d = fromIntegral $ natVal @d Proxy

-- | Insert along a dimension at a position.
--
-- >>> insert (Proxy :: Proxy 2) (Proxy :: Proxy 0) a ([100..105])
-- [[[100, 1, 2, 3, 4],
--   [101, 5, 6, 7, 8],
--   [102, 9, 10, 11, 12]],
--  [[103, 13, 14, 15, 16],
--   [104, 17, 18, 19, 20],
--   [105, 21, 22, 23, 24]]]
insert ::
  forall a s s' d i.
  ( DropIndex s d ~ s',
    CheckInsert d i s,
    KnownNat i,
    KnownNat d,
    HasShape s,
    HasShape s',
    HasShape (Insert d s),
    H.Element a,
    H.Container H.Vector a
  ) =>
  Proxy d ->
  Proxy i ->
  Array s a ->
  Array s' a ->
  Array (Insert d s) a
insert _ _ a b = tabulate go
  where
    go s
      | s !! d == i = index b (dropIndex s d)
      | s !! d < i = index a s
      | otherwise = index a (decAt d s)
    d = fromIntegral $ natVal @d Proxy
    i = fromIntegral $ natVal @i Proxy

-- | Insert along a dimension at the end.
--
-- >>>  :t append (Proxy :: Proxy 0) a
-- append (Proxy :: Proxy 0) a
--   :: Array '[3, 4] Int64 -> Array '[3, 3, 4] Int64
append ::
  forall a d s s'.
  ( DropIndex s d ~ s',
    CheckInsert d (Dimension s d - 1) s,
    KnownNat (Dimension s d - 1),
    KnownNat d,
    HasShape s,
    HasShape s',
    HasShape (Insert d s),
    H.Element a,
    H.Container H.Vector a
  ) =>
  Proxy d ->
  Array s a ->
  Array s' a ->
  Array (Insert d s) a
append d = insert d (Proxy :: Proxy (Dimension s d - 1))

-- | Change the order of dimensions.
--
-- >>> let r = reorder (Proxy :: Proxy '[2,0,1]) a
-- >>> :t r
-- r :: Array '[4, 2, 3] Int64
reorder ::
  forall a ds s.
  ( HasShape ds,
    HasShape s,
    HasShape (Reorder s ds),
    CheckReorder ds s,
    H.Element a,
    H.Container H.Vector a
  ) =>
  Proxy ds ->
  Array s a ->
  Array (Reorder s ds) a
reorder _ a = tabulate go
  where
    go s = index a (addIndexes [] ds s)
    ds = shapeVal (toShape @ds)

-- | Product two arrays using the supplied binary function.
--
-- For context, if the function is multiply, and the arrays are tensors,
-- then this can be interpreted as a tensor product.
--
-- https://en.wikipedia.org/wiki/Tensor_product
--
-- The concept of a tensor product is a dense crossroad, and a complete treatment is elsewhere.  To quote:
--
-- ... the tensor product can be extended to other categories of mathematical objects in addition to vector spaces, such as to matrices, tensors, algebras, topological vector spaces, and modules. In each such case the tensor product is characterized by a similar universal property: it is the freest bilinear operation. The general concept of a "tensor product" is captured by monoidal categories; that is, the class of all things that have a tensor product is a monoidal category.
--
-- >>> expand (*) v v
-- [[1, 2, 3],
--  [2, 4, 6],
--  [3, 6, 9]]
expand ::
  forall s s' a b c.
  ( HasShape s,
    HasShape s',
    HasShape ((++) s s'),
    H.Element a,
    H.Container H.Vector a,
    H.Element b,
    H.Container H.Vector b,
    H.Element c
  ) =>
  (a -> b -> c) ->
  Array s a ->
  Array s' b ->
  Array ((++) s s') c
expand f a b = tabulate (\i -> f (index a (take r i)) (index b (drop r i)))
  where
    r = rank (shape a)

-- | Select elements along positions in every dimension.
--
-- >>> let s = slice (Proxy :: Proxy '[[0,1],[0,2],[1,2]]) a
-- >>> :t s
-- s :: Array '[2, 2, 2] Int64
--
-- >>> s
-- [[[2, 3],
--   [10, 11]],
--  [[14, 15],
--   [22, 23]]]
--
-- >>> let s = squeeze $ slice (Proxy :: Proxy '[ '[0], '[0], '[0]]) a
-- >>> :t s
-- s :: Array '[] Int64
--
-- >>> s
-- 1
slice ::
  forall (pss :: [[Nat]]) s s' a.
  ( HasShape s,
    HasShape s',
    KnownNatss pss,
    KnownNat (Rank pss),
    s' ~ Ranks pss,
    H.Element a,
    H.Container H.Vector a
  ) =>
  Proxy pss ->
  Array s a ->
  Array s' a
slice pss a = tabulate go
  where
    go s = index a (zipWith (!!) pss' s)
    pss' = natValss pss

-- | Remove single dimensions.
--
-- >>> let a = [1..24] :: Array '[2,1,3,4,1] Int64
-- >>> a
-- [[[[[1],
--     [2],
--     [3],
--     [4]],
--    [[5],
--     [6],
--     [7],
--     [8]],
--    [[9],
--     [10],
--     [11],
--     [12]]]],
--  [[[[13],
--     [14],
--     [15],
--     [16]],
--    [[17],
--     [18],
--     [19],
--     [20]],
--    [[21],
--     [22],
--     [23],
--     [24]]]]]
-- >>> squeeze a
-- [[[1, 2, 3, 4],
--   [5, 6, 7, 8],
--   [9, 10, 11, 12]],
--  [[13, 14, 15, 16],
--   [17, 18, 19, 20],
--   [21, 22, 23, 24]]]
--
-- >>> squeeze ([1] :: Array '[1,1] Double)
-- 1.0
squeeze ::
  forall s t a.
  (t ~ Squeeze s) =>
  Array s a ->
  Array t a
squeeze (Array x) = Array x

-- $scalar
-- Scalar specialisations

-- | <https://en.wikipedia.org/wiki/Scalarr_(mathematics) Wiki Scalar>
--
-- An Array '[] a despite being a Scalar is never-the-less a one-element vector under the hood. Unification of representation is unexplored.
type Scalar a = Array ('[] :: [Nat]) a

-- | Unwrapping scalars is probably a performance bottleneck.
--
-- >>> let s = [3] :: Array ('[] :: [Nat]) Int64
-- >>> fromScalar s
-- 3
fromScalar :: (H.Element a, H.Container H.Vector a, HasShape ('[] :: [Nat])) => Array ('[] :: [Nat]) a -> a
fromScalar a = index a ([] :: [Int])

-- | Convert a number to a scalar.
--
-- >>> :t toScalar 2
-- toScalar 2 :: Num a => Array '[] a
toScalar :: (H.Element a, H.Container H.Vector a, HasShape ('[] :: [Nat])) => a -> Array ('[] :: [Nat]) a
toScalar a = fromList [a]

-- | <https://en.wikipedia.org/wiki/Vector_(mathematics_and_physics) Wiki Vector>
type Vector s a = Array '[s] a

-- | <https://en.wikipedia.org/wiki/Matrix_(mathematics) Wiki Matrix>
type Matrix m n a = Array '[m, n] a

instance
  ( Multiplicative a,
    P.Distributive a,
    Subtractive a,
    H.Numeric a,
    KnownNat m,
    HasShape '[m, m],
    H.Element a,
    H.Container H.Vector a
  ) =>
  Multiplicative (Matrix m m a)
  where
  (*) = mmult

  one = ident

-- | Extract specialised to a matrix.
--
-- >>> row 1 m
-- [4, 5, 6, 7]
row :: forall m n a. (H.Element a, H.Container H.Vector a, KnownNat m, KnownNat n, HasShape '[m, n]) => Int -> Matrix m n a -> Vector n a
row i (Array a) = fromList $ H.toList $ H.subVector (i * n) n (H.flatten a)
  where
    n = fromIntegral $ natVal @n Proxy

-- | Row extraction checked at type level.
--
-- >>> safeRow (Proxy :: Proxy 1) m
-- [4, 5, 6, 7]
--
-- >>> safeRow (Proxy :: Proxy 3) m
-- ...
-- ... index outside range
-- ...
safeRow :: forall m n a j. (H.Element a, H.Container H.Vector a, 'True ~ CheckIndex j m, KnownNat j, KnownNat m, KnownNat n, HasShape '[m, n]) => Proxy j -> Matrix m n a -> Vector n a
safeRow _j (Array a) = fromList $ H.toList $ H.subVector (j * n) n (H.flatten a)
  where
    n = fromIntegral $ natVal @n Proxy
    j = fromIntegral $ natVal @j Proxy

-- | Extract specialised to a matrix.
--
-- >>> col 1 m
-- [1, 5, 9]
col :: forall m n a. (H.Element a, H.Container H.Vector a, KnownNat m, KnownNat n, HasShape '[m, n]) => Int -> Matrix m n a -> Vector n a
col i (Array a) = Array $ H.takeColumns i a

-- | Column extraction checked at type level.
--
-- >>> safeCol (Proxy :: Proxy 1) m
-- [1, 5, 9]
--
-- >>> safeCol (Proxy :: Proxy 4) m
-- ...
-- ... index outside range
-- ...
safeCol :: forall m n a j. (H.Element a, H.Container H.Vector a, 'True ~ CheckIndex j n, KnownNat j, KnownNat m, KnownNat n, HasShape '[m, n]) => Proxy j -> Matrix m n a -> Vector n a
safeCol _j (Array a) = Array $ H.takeColumns j a
  where
    j = fromIntegral $ natVal @j Proxy

-- | Matrix multiplication.
--
-- This is dot sum (*) specialised to matrices
--
-- >>> let a = [1, 2, 3, 4] :: Array '[2, 2] Int64
-- >>> let b = [5, 6, 7, 8] :: Array '[2, 2] Int64
-- >>> a
-- [[1, 2],
--  [3, 4]]
--
-- >>> b
-- [[5, 6],
--  [7, 8]]
--
-- >>> mmult a b
-- [[19, 22],
--  [43, 50]]
mmult ::
  forall m n k a.
  ( KnownNat k,
    KnownNat m,
    KnownNat n,
    HasShape [m, n],
    Ring a,
    H.Numeric a
  ) =>
  Array [m, k] a ->
  Array [k, n] a ->
  Array [m, n] a
mmult (Array x) (Array y) = Array $ x H.<> y


text:
{-# LANGUAGE TemplateHaskell, GeneralizedNewtypeDeriving, Rank2Types, NoMonomorphismRestriction #-}
import Test.QuickCheck
import Test.QuickCheck.Gen.Unsafe
import Data.List
import Data.Int
import Data.Word
import Data.Version
import System.Exit
import Data.Complex
import Text.ParserCombinators.ReadP (readP_to_S)

newtype Path a = Path [a] deriving (Show, Functor)

instance Arbitrary a => Arbitrary (Path a) where
  arbitrary = do
    x <- arbitrary
    fmap Path (pathFrom 100 x)
    where
      pathFrom n x =
        fmap (x:) $
        case shrink x of
          [] -> return []
          _ | n == 0 -> return []
          ys -> oneof [pathFrom (n-1) y | y <- ys]

  shrink (Path xs) = map Path [ ys | ys <- inits xs, length ys > 0 && length ys < length xs ]

path :: (a -> Bool) -> Path a -> Bool
path p (Path xs) = all p xs

somePath :: (a -> Bool) -> Path a -> Property
somePath p = expectFailure . withMaxSuccess 1000 . path (not . p)

newtype Extremal a = Extremal { getExtremal :: a } deriving (Show, Eq, Ord, Num, Enum, Real, Integral)

instance (Arbitrary a, Bounded a) => Arbitrary (Extremal a) where
  arbitrary =
    fmap Extremal $
    frequency
      [(1, return minBound),
       (1, return maxBound),
       (8, arbitrary)]
  shrink (Extremal x) = map Extremal (shrink x)

smallProp :: Integral a => Path a -> Bool
smallProp = path (\x -> (x >= -100 || -100 `asTypeOf` x >= 0) && x <= 100)

largeProp :: Integral a => Path a -> Property
largeProp = somePath (\x -> x < -1000000 || x > 1000000)

prop_int :: Path Int -> Bool
prop_int = smallProp

prop_int32 :: Path Int32 -> Property
prop_int32 = largeProp

prop_word :: Path Word -> Property
prop_word = largeProp

prop_word32 :: Path Word32 -> Property
prop_word32 = largeProp

prop_integer :: Path Integer -> Bool
prop_integer = smallProp

prop_small :: Path (Small Int) -> Bool
prop_small = smallProp

prop_large :: Path (Large Int) -> Property
prop_large = largeProp

prop_smallWord :: Path (Small Word) -> Bool
prop_smallWord = smallProp

prop_largeWord :: Path (Large Word) -> Property
prop_largeWord = largeProp

data Choice a b = Choice a b deriving Show
instance (Arbitrary a, Arbitrary b) => Arbitrary (Choice a b) where
  arbitrary = do
    Capture eval <- capture
    return (Choice (eval arbitrary) (eval arbitrary))

idemProp :: (Eq a, Arbitrary a, Arbitrary b) => (b -> a) -> Choice a b -> Bool
idemProp f (Choice x y) = x == f y

prop_fixed_length :: Arbitrary a => Path (Fixed a) -> Bool
prop_fixed_length (Path xs) = length xs == 1

prop_fixed_idem = idemProp getFixed
prop_blind_idem = idemProp getBlind

prop_ordered_list = path (\(Ordered xs) -> sort xs == xs)
prop_nonempty_list = path (\(NonEmpty xs) -> not (null xs))

pathInt, somePathInt ::
  (Arbitrary (f (Extremal Int)), Show (f (Extremal Int)),
   Arbitrary (f Integer), Show (f Integer),
   Arbitrary (f (Extremal Int8)), Show (f (Extremal Int8)),
   Arbitrary (f (Extremal Int16)), Show (f (Extremal Int16)),
   Arbitrary (f (Extremal Int32)), Show (f (Extremal Int32)),
   Arbitrary (f (Extremal Int64)), Show (f (Extremal Int64)),
   Arbitrary (f (Extremal Word)), Show (f (Extremal Word)),
   Arbitrary (f (Extremal Word8)), Show (f (Extremal Word8)),
   Arbitrary (f (Extremal Word16)), Show (f (Extremal Word16)),
   Arbitrary (f (Extremal Word32)), Show (f (Extremal Word32)),
   Arbitrary (f (Extremal Word64)), Show (f (Extremal Word64))) =>
  Bool -> (forall a. f a -> a) -> (forall a. Integral a => a -> Bool) -> Property
pathInt word f p =
  conjoin
    [counterexample "Int" (path ((p :: Int -> Bool) . getExtremal . f)),
     counterexample "Integer" (path ((p :: Integer -> Bool) . f)),
     counterexample "Int8" (path ((p :: Int8 -> Bool) . getExtremal . f)),
     counterexample "Int16" (path ((p :: Int16 -> Bool) . getExtremal . f)),
     counterexample "Int32" (path ((p :: Int32 -> Bool) . getExtremal . f)),
     counterexample "Int64" (path ((p :: Int64 -> Bool) . getExtremal . f)),
     counterexample "Word" (not word .||. path ((p :: Word -> Bool) . getExtremal . f)),
     counterexample "Word8" (not word .||. path ((p :: Word8 -> Bool) . getExtremal . f)),
     counterexample "Word16" (not word .||. path ((p :: Word16 -> Bool) . getExtremal . f)),
     counterexample "Word32" (not word .||. path ((p :: Word32 -> Bool) . getExtremal . f)),
     counterexample "Word64" (not word .||. path ((p :: Word64 -> Bool) . getExtremal . f))]
somePathInt word f p = expectFailure (pathInt word f (not . p))

prop_positive = pathInt True getPositive (> 0)
prop_positive_bound = somePathInt True getPositive (== 1)

prop_nonzero = pathInt True getNonZero (/= 0)
prop_nonzero_bound_1 = somePathInt True getNonZero (== 1)
prop_nonzero_bound_2 = somePathInt True getNonZero (== -1)

prop_nonnegative = pathInt True getNonNegative (>= 0)
prop_nonnegative_bound = somePathInt True getNonNegative (== 0)

prop_negative = pathInt False getNegative (< 0)
prop_negative_bound = somePathInt False getNegative (== -1)

prop_nonpositive = pathInt True getNonPositive (<= 0)
prop_nonpositive_bound = somePathInt True getNonPositive (== 0)

reachesBound :: (Bounded a, Integral a, Arbitrary a) =>
  a -> Property
reachesBound x = withMaxSuccess 1000 (expectFailure (x < 3 * (maxBound `div` 4)))

prop_reachesBound_Int8 = reachesBound :: Int8 -> Property
prop_reachesBound_Int16 = reachesBound :: Int16 -> Property
prop_reachesBound_Int32 = reachesBound :: Int32 -> Property
prop_reachesBound_Int64 = reachesBound :: Int64 -> Property
prop_reachesBound_Word8 = reachesBound :: Word8 -> Property
prop_reachesBound_Word16 = reachesBound :: Word16 -> Property
prop_reachesBound_Word32 = reachesBound :: Word32 -> Property
prop_reachesBound_Word64 = reachesBound :: Word64 -> Property

-- Shrinking should not loop.
noShrinkingLoop :: (Eq a, Arbitrary a) => Path a -> Bool
noShrinkingLoop (Path (x:xs)) = x `notElem` xs

prop_no_shrinking_loop_Unit = noShrinkingLoop :: Path () -> Bool
prop_no_shrinking_loop_Bool = noShrinkingLoop :: Path Bool -> Bool
prop_no_shrinking_loop_Ordering = noShrinkingLoop :: Path Ordering -> Bool
prop_no_shrinking_loop_Maybe = noShrinkingLoop :: Path (Maybe Int) -> Bool
prop_no_shrinking_loop_Either = noShrinkingLoop :: Path (Either Int Int) -> Bool
prop_no_shrinking_loop_List = noShrinkingLoop :: Path [Int] -> Bool
prop_no_shrinking_loop_Ratio = noShrinkingLoop :: Path Rational -> Bool
prop_no_shrinking_loop_Complex = noShrinkingLoop :: Path (Complex Double) -> Bool
prop_no_shrinking_loop_Fixed = noShrinkingLoop :: Path (Fixed Int) -> Bool
prop_no_shrinking_loop_Pair = noShrinkingLoop :: Path (Int, Int) -> Bool
prop_no_shrinking_loop_Triple = noShrinkingLoop :: Path (Int, Int, Int) -> Bool
prop_no_shrinking_loop_Integer = noShrinkingLoop :: Path Integer -> Bool
prop_no_shrinking_loop_Int = noShrinkingLoop :: Path Int -> Bool
prop_no_shrinking_loop_Int8 = noShrinkingLoop :: Path Int8 -> Bool
prop_no_shrinking_loop_Int16 = noShrinkingLoop :: Path Int16 -> Bool
prop_no_shrinking_loop_Int32 = noShrinkingLoop :: Path Int32 -> Bool
prop_no_shrinking_loop_Int64 = noShrinkingLoop :: Path Int64 -> Bool
prop_no_shrinking_loop_Word = noShrinkingLoop :: Path Word -> Bool
prop_no_shrinking_loop_Word8 = noShrinkingLoop :: Path Word8 -> Bool
prop_no_shrinking_loop_Word16 = noShrinkingLoop :: Path Word16 -> Bool
prop_no_shrinking_loop_Word32 = noShrinkingLoop :: Path Word32 -> Bool
prop_no_shrinking_loop_Word64 = noShrinkingLoop :: Path Word64 -> Bool
prop_no_shrinking_loop_Char = noShrinkingLoop :: Path Char -> Bool
prop_no_shrinking_loop_Float = noShrinkingLoop :: Path Float -> Bool
prop_no_shrinking_loop_Double = noShrinkingLoop :: Path Double -> Bool
prop_no_shrinking_loop_Version = noShrinkingLoop :: Path Version -> Bool
prop_no_shrinking_loop_ExitCode = noShrinkingLoop :: Path ExitCode -> Bool

-- Check that shrinking a Double always produces a shrinking candidate.
prop_shrink_candidate_double :: Property
prop_shrink_candidate_double =
  forAllShrink gen shrink $ \x ->
    x > 0 ==>
    not (null (shrink x))
  where
    gen :: Gen Double
    gen = oneof [arbitrary, fmap fromInteger arbitrary]

-- Bad shrink: infinite list
--
-- remove unexpectedFailure in prop_B1, shrinking should not loop forever.
data B1 = B1 Int deriving (Eq, Show)

instance Arbitrary B1 where
    arbitrary = fmap B1 arbitrary
    shrink x = x : shrink x

prop_B1 :: B1 -> Property
prop_B1 (B1 n) = expectFailure $ n === n + 1

return []
main = do True <- $forAllProperties (quickCheckWithResult stdArgs { maxShrinks = 10000 }); return ()

ext: hs

max_stars_repo_name: isovector/quickcheck

max_issues_repo_name: isovector/quickcheck

max_forks_repo_name: isovector/quickcheck


text:
{-# LANGUAGE BangPatterns #-}

import Gauge.Main (defaultMain, bench, whnf)

import Data.List (unfoldr)
import Data.Maybe (fromJust)
import Data.Primitive.Contiguous (PrimArray, fromListN)
import Statistics.Array.Types (AscList)
import System.Random (random, mkStdGen)

import qualified Statistics.Array as Stats
import qualified Statistics.Array.Types as Asc

main :: IO ()
main = defaultMain
  [ bench "sort-16" $ whnf (fromJust . Asc.fromArray) input16
  , bench "sort-128" $ whnf (fromJust . Asc.fromArray) input128
  , bench "sort-1024" $ whnf (fromJust . Asc.fromArray) input1024
  , bench "sort-65536" $ whnf (fromJust . Asc.fromArray) input65536
  , bench "list-derivative" $ whnf Stats.listDerivative asc1024
  , bench "median-of-absolute-deviations" $ whnf Stats.mad asc1024
  ]
  where
  input16, input128, input1024, input65536 :: PrimArray Int
  !input16 = fromListN 16 $ take 16 . drop 0 $ inputInf
  !input128 = fromListN 128 $ take 128 . drop 16 $ inputInf
  !input1024 = fromListN 1024 $ take 1024 . drop (16+128) $ inputInf
  !input65536 = fromListN 65536 $ take 65536 . drop (16+128+1024) $ inputInf
  asc1024 :: AscList PrimArray Int
  asc1024 = fromJust $ Asc.fromArray input1024

inputInf :: [Int]
inputInf = unfoldr (Just . random) seed
  where
  seed = mkStdGen 1234567

max_stars_repo_name: andrewthad/array-statistics

max_issues_repo_name: andrewthad/array-statistics

max_forks_repo_name: andrewthad/array-statistics


text:
{-# OPTIONS_GHC -Wno-orphans #-}

module Data.Function.FastMemo.Instances () where

import Data.Complex (Complex)
import Data.Function.FastMemo.Char ()
import Data.Function.FastMemo.Class (Memoizable)
import Data.Function.FastMemo.Int ()
import Data.Functor.Identity (Identity)
import Data.Proxy (Proxy)
import Data.Version (Version)
import Data.Void (Void)

instance Memoizable Bool

instance Memoizable a => Memoizable (Maybe a)

instance (Memoizable a, Memoizable b) => Memoizable (Either a b)

instance Memoizable Void

instance Memoizable ()

instance Memoizable (Proxy a)

instance Memoizable Ordering

instance Memoizable a => Memoizable (Complex a)

instance Memoizable Version

instance Memoizable a => Memoizable (Identity a)

instance (Memoizable a, Memoizable b, Memoizable c) => Memoizable (a, b, c)

instance (Memoizable a, Memoizable b, Memoizable c, Memoizable d) => Memoizable (a, b, c, d)

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

instance (Memoizable a, Memoizable b, Memoizable c, Memoizable d, Memoizable e, Memoizable f) => Memoizable (a, b, c, d, e, f)

instance (Memoizable a, Memoizable b, Memoizable c, Memoizable d, Memoizable e, Memoizable f, Memoizable g) => Memoizable (a, b, c, d, e, f, g)

max_stars_repo_name: davidspies/fastmemo

max_issues_repo_name: davidspies/fastmemo

max_forks_repo_name: davidspies/fastmemo


text:
-----------------------------------------------------------------------------
-- |
-- Module     : Numeric.LinearAlgebra.Packed.Cholesky
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- Cholesky factorizations of symmetric (Hermitian) positive-definite
-- matrices.
--

module Numeric.LinearAlgebra.Packed.Cholesky (
    -- * Immutable interface
    cholFactor,
    cholSolveVector,
    cholSolveMatrix,

    -- * Mutable interface
    cholFactorM,
    cholSolveVectorM_,
    cholSolveMatrixM_,
    ) where

import Control.Monad( when )
import Control.Monad.ST( ST, runST, unsafeIOToST )
import Text.Printf( printf )

import qualified Foreign.LAPACK as LAPACK

import Numeric.LinearAlgebra.Types

import Numeric.LinearAlgebra.Packed.Base( Packed, RPacked, STPacked )
import qualified Numeric.LinearAlgebra.Packed.Base as P

import Numeric.LinearAlgebra.Matrix( Matrix, STMatrix )
import qualified Numeric.LinearAlgebra.Matrix as M

import Numeric.LinearAlgebra.Vector( Vector, STVector )
import qualified Numeric.LinearAlgebra.Vector as V

-- | @cholFactor a@ tries to compute the Cholesky
-- factorization of @a@.  If @a@ is positive-definite then the routine
-- returns @Right@ with the factorization.  If the leading minor of order @i@
-- is not positive-definite then the routine returns @Left i@.
cholFactor :: (LAPACK e)
                 => Herm Packed e
                 -> Either Int (Chol Packed e)
cholFactor (Herm uplo a) = runST $ do
    ma <- P.newCopy a
    cholFactorM (Herm uplo ma)
        >>= either (return . Left) (\(Chol uplo' ma') -> do
                a' <- P.unsafeFreeze ma'
                return $ Right (Chol uplo' a')
                )

-- | @cholSolveVector a x@ returns @a \\ x@.
cholSolveVector :: (LAPACK e)
                      => Chol Packed e
                      -> Vector e
                      -> Vector e
cholSolveVector a x = V.create $ do
    x' <- V.newCopy x
    cholSolveVectorM_ a x'
    return x'

-- | @cholSolveMatrix a b@ returns @a \\ b@.
cholSolveMatrix :: (LAPACK e)
                      => Chol Packed e
                      -> Matrix e
                      -> Matrix e
cholSolveMatrix a c = M.create $ do
    c' <- M.newCopy c
    cholSolveMatrixM_ a c'
    return c'

-- | @cholFactorM a@ tries to compute the Cholesky
-- factorization of @a@ in place.  If @a@ is positive-definite then the
-- routine returns @Right@ with the factorization, stored in the same
-- memory as @a@.  If the leading minor of order @i@ is not
-- positive-definite then the routine returns @Left i@.
-- In either case, the original storage of @a@ is destroyed.
cholFactorM :: (LAPACK e)
             => Herm (STPacked s) e
             -> ST s (Either Int (Chol (STPacked s) e))
cholFactorM (Herm uplo a) = do
    n <- P.getDim a
    unsafeIOToST $
        P.unsafeWith a $ \pa -> do
            info <- LAPACK.pptrf uplo n pa
            return $ if info > 0 then Left info
                                 else Right (Chol uplo a)

-- | @cholSolveVectorM_ a x@ sets @x := a \\ x@.
cholSolveVectorM_ :: (LAPACK e, RPacked p)
                  => Chol p e
                  -> STVector s e
                  -> ST s ()
cholSolveVectorM_ a x =
    M.withFromColM x $ \x' ->
        cholSolveMatrixM_ a x'

-- | @cholSolveMatrixM_ a b@ sets @b := a \\ b@.
cholSolveMatrixM_ :: (LAPACK e, RPacked p)
                  => Chol p e
                  -> STMatrix s e
                  -> ST s ()
cholSolveMatrixM_ (Chol uplo a) b = do
    na <- P.getDim a
    (mb,nb) <- M.getDim b
    let (n,nrhs) = (mb,nb)

    when ((not . and) [ na == n
                      , (mb,nb) == (n,nrhs)
                      ]) $ error $
        printf ("cholSolveMatrixM_"
                ++ " (Chol _ <packed matrix with dim %d>)"
                ++ " <matrix with dim (%d,%d)>"
                ++ ": dimension mismatch")
               na mb nb

    unsafeIOToST $
        P.unsafeWith a $ \pa ->
        M.unsafeWith b $ \pb ldb ->
            LAPACK.pptrs uplo n nrhs pa pb ldb

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
{-# LANGUAGE BangPatterns #-}
--{-# XRankNTypes #-}

module RprGD2 (OptimisationType (Rapid , FromData ,WithDepth)
    , mmlRprGdEvl,mml2DsGD -- ,mmlRprTest
    , mml2DsGen, RprMix2.aBar,mapP --  , stdTest
    --,mmlRprTest
    ) where

import Numeric.LinearAlgebra
import Numeric.GSL.Statistics (stddev,lag1auto)
import qualified Data.ByteString.Lazy.Char8 as L (writeFile, pack, unpack, intercalate,appendFile, append,concat)
import System.Random
import Control.Concurrent
import Control.DeepSeq
import Foreign.Storable
import Control.Parallel (par, pseq)
import qualified Control.Parallel.Strategies as St-- (parMap)
import System.IO.Unsafe
import Data.Maybe (isJust,fromMaybe,maybe,listToMaybe,fromJust)
import Data.Time
import Data.Char (isDigit)
--import Text.Printf
import Data.Ord (comparing)
import Control.Monad (liftM,liftM2)
import Data.List (transpose, foldr1, sort, foldl' ,nub, nubBy,sortBy,  minimumBy, partition)
import qualified Data.Map  as M (Map,empty,insertWith' ,mapWithKey,filterWithKey,toList, fromList)
--import IO
--------------------------------------------------------------------------------------------
import HMDIFPrelude (bsReadChr,bs2Int')
import RprMix2 -- (vLength)
import LINreg2
import ResultTypes
import ProcessFile (writeToFile,log2File,writeOrAppend,myTime)
import LMisc hiding (force)
import ListStats (nCr,normal,binomial ,myTake,-- corr,
                  mean, adjustedMatrix,cdf, meanBy, stDivBy, -- stDiv,
                  adjustedMatrixN, adjustedMatrixPN,stNorm,mean,stDivM)
---------------------------------------------------------------------------------------------
----                 optimisation type
data OptimisationType = Rapid | FromData | WithDepth Int deriving (Show, Eq)-- | WithNum Int
---------------------------------------------------------------------------------------------
---- Paralell definitions
concatPMap ::  (a -> [b]) -> [a] -> [b]
concatPMap f = concat . St.parMap (St.parList St.rseq) f

-- parallel map and parallel concatMap
--mapP  = St.parMap (St.dot St.rseq St.rpar)

mapP' :: (a -> b) -> [a] -> [b]
mapP' f = foldl' (\bs a -> let fa = (f a) `St.using`  St.rseq in
                                fa : bs ) []
--
mapP :: (a -> b) -> [a] -> [b]
mapP  f = St.parMap St.rseq f

concatPMap' f   = foldl'Rnf  (\xs x -> xs ++ [St.using (f x) (St.dot St.rseq St.rpar) ] ) []

--- strict fold vector
-- foldVector :: Storable a => (a -> b -> b) -> b -> Vector a -> b
foldVector' :: (NFData b, Storable a) => (b -> a -> b) -> b -> Vector a -> b
foldVector' f a vs = foldVector (\b g x -> g (seqit $ f x b)) id vs a
    where seqit a = rnf a `pseq` a
--------------------------------------------------------------------------------------------------
myMinimum :: (NFData a, Ord a) => [a] -> a
myMinimum   (x:xs) = foldl'Rnf  min x xs
myMinimum     _    = error "error: myMinimum empty list"
--
--mmlRprGD :: Vector Double -> Double -> Double -> Matrix Double -> Mmodel -> (Mmodel, Double)
{-# NOINLINE mmlRprGD #-}
mmlRprGD g ys n cts nSLP dss (ml, model) = mmlRprGD_aux g ys n cts nSLP dss  model

mmlRprGD_aux :: (RandomGen g) => g
     -> Vector Double
     -> Double
     -> [Int]
     -> Maybe Bool
     -> Matrix Double
     -> Mmodel
     -> IO OModel
mmlRprGD_aux g ys n cts nSLP dss model =  do
    return $ Omod (nAbaMod model taU) oldMsgLn
        where
            ------- compute values from the input model, eg the weighs-----------
            nAbaMod  (MkModel a b c x _ )  k  =  (MkModel a b c x k)
            ----
            fL             = (+ 1) . fromIntegral . length
            aA'            = fL cts
            -------------------------------------------------------
            siG            =  sigD  model
            siGE           =  sigE  model
            taU            =  tau   model
            ---------------------------------------------------------
            xXa            =  invS nSLP (dss <> (aBar model)) ---
            abar           = aBar model

            aA             =  fromIntegral . (\ a -> if a == 2 then 1 else a `div` 2)  $ dim  abar
            errFun  a b    =  (r_n model a b)
            rN             =   zipVectorWith  errFun xXa ys -- (r_n model)
            lnR            =   dim rN   -- get the magnitude of rN
            wData :: Matrix Double
            wData          =  diagRect 0 rN lnR lnR
            wError :: Matrix Double
            wError         =  diagRect 0 (mapVector (1 -) rN) lnR lnR
            ---
            xtWx            =  ((trans dss) <> wData) <> dss
            ------------------------------------------------------------
            ---
            nAlpha         =  foldVector (+) 0 rN
            nAlpA1         =  nAlpha - aA - 1 --
            --
            ysSubXa        =  zipVectorWith (-) ys  xXa -- the data minus the predictions
            ySubWySub ww   =  ysSubXa <.> (ww <> ysSubXa) -- y
            --
            dehAbarSq      = (abar <.> abar)
            aBarSqTau1     =   dehAbarSq / (2 * taU^2)
            ----------------------------------------------------------------------------------------
            ------------------- calculate the message length of the input model --------------------
            ----------------------------------------------------------------------------------------
            !oldMsgLn     =   l1 `par` l2 `par` l3 `par` l4 `pseq`   (l1 + l2 + l3 + l4)
            ----------------------------------log 0.01 *------------------------------------------------------
            restrPr       = ((aA + 1)* log 2)/ 2
            restrPrior     = maybe restrPr (\_ -> restrPr) nSLP
            l1             =   (aA + 2.0) * logg taU + aBarSqTau1 - restrPrior  + 0.5 * ((n-2) * log (2 * pi))
            ----------------------------------- log 2 * (aA + 1)/2 ---------------------------------
            l2              =   nAlpA1 * logg siG + (1/(2 * siG^2)) * (ySubWySub wData)
            l3              =  (n - nAlpha) * logg siGE + (1/(2 * siGE^2)) * (ySubWySub wError)  + 0.5 * logg (2 * nAlpha)
            mlLogs          = 0.5 * log (det xtWx) + 0.5 * logg (2*(n - nAlpha)) +  0.5 * log (pi * (aA + 3))
            l4              =   mlLogs - 1.319
----------------------------------------------------------------------------------------------------
dummyMod       = Omod (MkModel 0 0 0 (fromList [1]) 0) 100000000
vunStdNorm' vs = vunStdNorm (vmean vs) (vstDiv vs) vs
{----------------------------------------------------------------------------------------------------
mmlRprGdEvl: the mmlRPR  evaluation function. Takes:
n  - the range in years of roughness data (which corresponds to the number of data points
      when dealing with simulation data)
ms - the number of chains that are joined together - the maintenance strategy.
     Note that the mml2DS uses the mmlRPR on adjaent sections joined together. Hence this
     parameter says how many sections are to be joined for calls to mmlRprGdEvl. When this
     functions is called outside of the mml2DS, this valus is always 1 (i.e we do not make any
     prior assumptions about maintenances on the chain)
ys  - a lsit of the data values for this chain
pp  - the interval in years between maintenenace interventions
gen - a random generator to initialize the process
nSlopes - determines whether to exclude negative slopes: if true, negative slopes are excluded
           otherwise they are not
rprOnly - toggles the applicaion of the MMLRPR function only. The likelihood of maintenance
          is not applied to each interventions
----------------------------------------------------------------------------------------------------}
mmlRprGdEvl :: (RandomGen g) => Int ->
    Int -> -- error search depth limit
    [Vector Double]  -> -- Either (Vector Double) [Double] -> -- the data in standardized for
    Double ->
    g ->
    Maybe Bool -> -- negative slopes (temporarily used for turning on and off the mixture model)
    Int ->
    IO (Maybe ( (OModel,(Vector Double,[Int]) ) , [Double])) --ms
mmlRprGdEvl n dpt ys pp gen nSlopes mtd = -- mN std
    mmlRprAuxGD n dpt m p (filter ((== mLen) . dim) ys) pp gen nSlopes mtd True -- ms -- mN std
        where
            p     =  1 / pp
            --- don need this claculation. we can filter out irelevant values
            m     =  n `div` mtd
            mLen  = maybe 0 (\_ -> maximum $ map dim ys) (listToMaybe ys)
            ---------------------------------------------------------------------------------------------------
            unStd :: (OModel,(Vector Double,[Int])) -> (OModel,(Vector Double,[Int]))
            unStd  (Omod (MkModel c b e ks d)  ml,(vs,ys)) =  (Omod (MkModel c b e (vunStdNorm' ks) d)  ml,(vunStdNorm' vs,ys))
            unSd ms =
                case ms of
                    Nothing -> Nothing
                    Just (a,b) -> Just (unStd a, b)

------------
mmlRprAuxGD n _ m p ysv pp gen nSlopes mtd toOpt = do
    return minMod
    where
        minMod = liftM applyMixture $
                    maybe Nothing (\_ -> Just $ minimumBy (comparing (fst . fst))  rprs1) (listToMaybe rprs1)
        --
        result  rprss =  (snd (minimumBy (comparing (fst . snd)) rprss), map fst rprss)
        f         =  fromIntegral
        (g1,_)    =  split gen
        abr' ::  Matrix Double -> [(Double, (Vector Double, (Double, Vector Double))) ]
        abr'      = (: []) . minimumBy (comparing fst)
                    . zipWith (\f a -> f a) (map (\ys -> mmlLinearRegLP ys 1 1) setLen)
                    . replicate len -- (length ysv)
            where
                len    = length ysv
                setLen = replicate len $ join ysv
        ----------------------------------------------------------------------------
        initM b s =  initModel gen s b -
        initLz (ml, (prd, (s,ab))) =  (ml, (prd , initM ab s))
        --
        mod :: Matrix Double -> [(Double, (Vector Double, Mmodel))]
        mod  mxx  =  map   initLz $ abr'  mxx
        ----------------------------------------------------------------------------
        nN = f . sum $ map dim ysv -- n
        nNs = map ((\n -> [1 .. n]) . f . dim) ysv
        fLn = (+ 1) . f . length
        msgLm (Omod _ l) = l
        appMML g1 nM ct  mmX (ml,(prd , inMd))  = unsafePerformIO $ mmlRprGD g1 (join  ysv) nM ct nSlopes mmX (ml, inMd)
        iTr (mX, cts) = (map (appMML gen nN cts  mX) (mod mX), (mX, cts))
        cMML  = mapP (getPredictions nSlopes) . mapSnd . iTr
        nps cmm  =   realignedNpieces cmm   mtd
        --------------------------------------------------------------------------------------------
        ----------- applying the mixture model to the discovered intervention, only ----------------
        rprs1  =   [ (kx, cm1) | cm1 <- [0 .. m],  kx <- ( concatPMap (mapSnd . iTr) . (nps cm1)) nNs  ]  `St.using` (St.parList St.rseq)
        applyMixture ( (Omod md mln , (mx , cts)), cm) =
            let   (!nMd , !mxL)  = unsafePerformIO $ nextModelR' gen cts nSlopes mx (join ysv) md
                  ncma    = n
                  nncma   =  ncma `nCr` cm
                  ncr     =   (binomial n cm p) / f nncma -- else 1
                  lncr    =   log ncr
                  newOMod = Omod nMd ((lncr + mln) + mxL)
                  prds    = getPredictions nSlopes (newOMod, (mx, cts))
            in    (prds , [lncr]  )
        --------------   end applying the Mixture Model at the end ----------------------------------
        getPredictions nSlp (mm, (xConfig, cts))  =   (mm , (invS nSlp  prds , cts))
            where
               prds     =   xConfig <> (aBar (oMod mm))
        --- retain the domain of the function
        retDom :: (a -> [b]) -> a -> [(a , b)]
        retDom  f  a  =   map (\x -> (a, x)) ( f a)
        ---
        nAbaMod  (MkModel a b c _ x)  k  =  (MkModel a b c k x)

{----------------------------------------------------------------------------------------------------
  The mml2DsGD functions applies the MMLRPR to alogorithm to a list of data
  points (i.e. the readings or chainages) and after
  sectioning them off mml2DsGD years distance gen [[Double]]
--}
--- returning results
mml2DsGen :: (RandomGen g) => g -> TdsParms ->
    OptimisationType -> -- the optimisation type (Rapid , FromData ,WithDepth)
    SecOrSchm ->
    IO SecOrSchm
mml2DsGen g tprms optType ssm =
    case ssm of
        -- analysis results for sections
        Left  sect -> do
            let (a, b , c ,d) = (sCondition sect, sScannerCode sect, sScannerSorrg sect, sCWay sect)
            let mkSect = Section a b c d (sChains sect)
            liftM (Left . mkSect . omd2Results) $ mml2DsGD g tprms optType [xs | xs <- sChains sect, length (cData xs) > 1]
        -- results for schemes
        Right schm -> do
            nScmRecs <- mapM (mkSchemeRecRes g tprms optType) (schmRecs schm)
            return $ Right $ Scheme (schemeName schm) (schmCondition schm) nScmRecs
    where
        mkSchemeRecRes :: (RandomGen g) => g ->
            TdsParms ->
            OptimisationType ->
            SchemeRec ->
            IO SchemeRec
        mkSchemeRecRes g tp opT src = liftM (mkSchmRec . omd2Results) $ mml2DsGD g tp opT chn
            where
                mkSchmRec = SchemeRec scd srd srnm srg stc sec chn
                --
                scd = scmScannerCode src
                srd = roadClassCode src
                srnm = roadName  src
                srg  = scmScannerSorrg  src
                --- carriage way
                stc = startChain   src
                sec = endChain  src
                chn = scmChains src

--}
-- returning a OMod, etc
mml2DsGD :: (RandomGen g) => g ->
    TdsParms -> -- the parameters for the mmlTds
    OptimisationType -> -- the optimisation type (Rapid , FromData ,WithDepth)
    [Chain] -> --- a list of all the data for each chain in a section or scheme
    IO ([(OModel, ([Double], [Int]))], [Int])
mml2DsGD  gen tdPrs optType chxss = do
    ---------------------------------------------
    scrollLog ("Identifying candidate groups for section with length: " ++ show len) -- (max 3 (len `div` 15))
    !mks1  <-   (joinPairsV dist optType 2 len 25 xss)
    ------
    let total = length mks1
    let msgP1 = "Fitting progression rate and identifying outliers for section length: " ++ show len ++". "
    let msgP2 = show total ++": candidate groups found."
    onlycuts <- liftM (maybe 0 (fromJust . bs2Int') . listToMaybe . head) $ bsReadChr  ',' "./mmlTdsFiles/onlycuts.txt"
    -- print ("negslopes is : " ++ show nSlopes)
    if len > 0 then
        if null mks1 then  do
            let noGrupMsg = "NO Groups for section length: " ++ show len
            logLog ("-------------------"++noGrupMsg)
            return  dummy
        else   if onlycuts == 1 then do
            logLog ("only first cuts printed")
            return dummy
        else do
            --
            scrollLog  (msgP1 ++ msgP2)
            let mks2  = [(a,b,c) | ((a,b),c) <- zip (zip mks1 (rGens gen)) [1 ..]]
            let applyRpr = mapP (mml2DsGD_aux1 total minTrend scrollLog)
            let appMin  = liftM (minimumBy compML'')  . sequence .  applyRpr  -- ($!)
            appMin mks2
    else do
        logLog ("Empty section list sent to MML2DS: " ++ show len ++". limit: "++ show lenLim)
        return dummy
    where
        ---- projections --
        xss                  =  map cData chxss
        (!dpt, !pp , !dist)  =  (searchDept tdPrs, rngInYrs tdPrs, rngInDist tdPrs)
        (nSlopes,lim,minAnl) =  (negSlopes tdPrs, maxChn2Aln tdPrs, minChn2Aln tdPrs)
        minTrend             =  minPts2Trd tdPrs
        loggers              =  logFun tdPrs
        scrollLog            =  maybe (\_ -> return ())  (\_ -> loggers !! 0) (listToMaybe loggers)
        logLog               =  maybe (\_ -> return ())  (\_ -> loggers !! 1) (listToMaybe loggers)
        maxChain             =  (foldl' (\n xs -> max (length xs) n) 0 xss)
        lenLim               =  lim * maxChain
        len                  =  length xss
        minLen               =  minAnl * maxChain
        -------------------------------- end of projections ------------------
        dummy                =  ([(dummyMod , ([0] ,[0]))] ,[])
        pm                   =  1 / dist
        compML''             =  comparing  (sumBy (oModMl . fst). fst)
        f                    =  fromIntegral
        {--------------------------- logging the output from the meddage lengths --------------------
        printCutsInfo :: OptimisationType -> ([(OModel, ([Double], [Int]))], [Int]  ) -> IO ()
        printCutsInfo opTp (ks,cs) = do
            let cutsMsg = "\n cuts: "++ printElm cs
            let mlgMsg  = "; \t message length: " ++ (show $ sumBy (oModMl . fst) ks)
            let name = case opTp of
                            Rapid         ->  "fromData" -- "priorOnly"
                            FromData      ->  "fromData"
                            WithDepth _   ->   "fixedDepth"
            let file = "./SimulationFiles/"++name++"groups.txt"
            -- let header = "\n-----------\t cuts and message lengths ------------\n"
            let output = L.pack (cutsMsg ++ mlgMsg)
            writeOrAppend file output
        ---------------------------------------------------------------------------------------}
        mml2DsGD_aux1 :: RandomGen g =>  Int ->
                           Int ->
                         (String -> IO()) ->
                         ( ([([Vector Double], Int)],[Int]), g, Int)  ->
                         IO ([(OModel, ([Double], [Int]))], [Int])
        mml2DsGD_aux1 total minToTrend write ((pgroups,pgcuts),gen,num)  =  do
            write ("Calculating mmlrpr on group: "++ show num ++ " of " ++ show total)
            let max1  = sumBy snd pgroups -- sum pgcut
            ----
            let  calc = [ (result,pgcuts) |
                            let !rprsCalcs =  unsafePerformIO . sequence $ mapP (calcRpR max1 gen minToTrend) pgroups -- do the mlrpr on each grou
                            , all isJust rprsCalcs -- .  filter isJust
                            , let result   = (map fromJust . filter isJust) rprsCalcs]
            if null calc then do
                write ("Invalid mmlrpr group for group "++ show num)
                return dummy

            else do
                retMinComp calc
                where
                    --
                    retMinComp =  return . minimumBy  compML''
                    --
                    calcRpR :: RandomGen g => Int -> g -> Int ->
                               ([Vector Double], Int) -> IO (Maybe (OModel,([Double],[Int]))) -- , [Double]) -
                    calcRpR ns g mtd (vys, ms) = do
                        !rprs <- mmlRprGdEvl ms dpt vys pp g nSlopes mtd
                        case rprs of
                            Just ( (Omod r1 r2, prs),nn) -> do
                                       let smm = ncr' + r2
                                       let str = ("\nncr' is: "++ show ncr' ++" mLen is: "++ show r2 ++". Their sum is: " ++ show smm)
                                       -- L.appendFile "tempLog1.txt" (L.pack str)
                                       return $ Just (Omod r1 smm , (toList (fst prs), snd prs) ) --
                            Nothing ->   return $ Nothing
                        where
                            ncr'  =   (binomial ns ms pm) / f (ns `nCr` ms)
-----------------------------------------------------------------------------------------------------}
------------
--------------------------------------------------------------------------------------------------
findGrpsFromData :: Double -> -- the distance of intervention points
    Int              ->  -- lenght of the input list
    [[Double]] ->           -- ^ the data to aling
    IO [([([Vector Double], Int)], [Int])]
findGrpsFromData dist inLen   =
    return . findCandidateSet .  uncurry calculateFixedGroups1 . findInterventions_aux
    where
        findCandidateSet xs
            | length xs > 100 = maintainBy' avgRelCorr (Just 45) xs
            | otherwise       = (take 10 .  sortBy (comparing avgRelCorr )) xs
        --
        avgRelCorr    =  (meanBy (relativeCorr . fst)) . fst
        -- find possible interventions from the data
        findInterventions_aux ::  [[Double]] ->    ([[Double]] , [[Int]])
        findInterventions_aux  xys =   (xys ,  (nub . findGrps1 0.7) xys)
            where
                findGrps1 tolerance xs =  [ ys | let stDiv = stddev . fromList
                                                , ys <-  findGrps tolerance 38 mmlCorr xs
                                          ] `St.using` (St.parList St.rseq)
        -- calculate fixed groups
        calculateFixedGroups1 :: [[Double]] ->  [[Int]]  -> [([([Vector Double], Int)], [Int])]
        calculateFixedGroups1  xys  =  map (unzip . foldr mkLst [] . filter (not . null . fst) .  mkIntervals vs )
            where
                vs              = (map fromList xys)
                mkLst (ks,k) ps =  ((ks, sumBy dim ks) , k) :  ps
        -- relative correlaton between adjacent chains
        relativeCorr :: [Vector Double] -> Double
        relativeCorr xs
            | length xs < 2 =  4 -- bias against groups with only one chian
            | otherwise     =  mcD xs
            where
                corrDist ps qs =  abs (nCorrV ps qs - 1)
                mcD            =  meanBy (exp . uncurry corrDist) . zipp

----------------------------------------------------------------------------------------------------
--- an infinite list of random generators
rGens :: RandomGen g =>  g -> [g]
rGens g = g : (g1 : rGens g2)
    where (g1, g2) = split g

------------------------------------------------------------------------------------------------------------------------
joinPairsV :: Double -> -- the distance of maintenence interventions
            OptimisationType -> -- the optimisation type (Rapid , FromData ,WithDepth)
            Int ->              -- ^ the minimum number of chains to aligh
            --Int ->            -- ^ maximum
            Int ->              -- ^ length of the input list
            Int ->              -- ^ the maximum number of chains to align
            [[Double]] ->  -- ^ the data to aling
            IO [([([Vector Double], Int)], [Int])]
            --[([(Vector Double, Int, Double )],[Int])]
            -- (joinPairsV pm optType minLen lenLim 3 len lim) xss
joinPairsV dist optType malg xysLen lim xys -- minL maxL
            | null  xys           =  return []
            | xysLen <  4         =  return joinFew -- ' pm --xys
            --- check the length of the combined list
            | lenData < 10        =  return joinFew
            | otherwise           =
                case optType of
                    Rapid         ->  findGrpsFromData dist xysLen xys
                    FromData      ->  findGrpsFromData dist xysLen xys
                    WithDepth d   ->  calculateFixedGroups vxys [(repeat d)]
            where
                lenData           = sumBy length xys
                --
                vxys              =  map fromList xys
                ------------------
                vs                =  join vxys
                ms                =  dim vs
                joinFew           =  [([([vs],dim vs)],[])]

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

-- given a list of possible cuts and the list of chains, we generate the a list
-- of possible maintenance groups defeind by the cuts
calculateFixedGroups :: [Vector Double] ->  [[Int]]  -> IO [([([Vector Double], Int)], [Int])]
calculateFixedGroups  vs  = return . map (unzip . foldr mkLst [] . filter (not . null . fst) . mkIntervals vs)
    where
       mkLst (ks,k) ps =  ((ks, sumBy dim ks) , k) :  ps
       -- return [ (mkIntervals vs cs , cs)  | cs <- ks]

max_stars_repo_name: rawlep/MML

max_issues_repo_name: rawlep/MML

max_forks_repo_name: rawlep/MML


text:
module LibSpec where

import           Debug.Trace
import           Lib                        (sigmoid)
import           Numeric.LinearAlgebra      (( #> ))
import           Numeric.LinearAlgebra.Data (Matrix, Vector, cond, tr, vector,
                                             (><))
import           Test.Hspec                 hiding (Spec)
import           Test.Hspec                 as Hspec
import           Test.HUnit                 (Assertion, assertBool)

spec :: Hspec.Spec
spec = do
  describe "sigmoid" $
    it "Should be correct." $ do
      (sigmoid ((3 >< 2) [1.0 ..]) :: Matrix Double) `shouldBe`
        (((3 >< 2)
            [ 0.7310585786300049
            , 0.8807970779778823
            , 0.9525741268224334
            , 0.9820137900379085
            , 0.9933071490757153
            , 0.9975273768433653
            ]) :: Matrix Double)
  describe "cond behaviour" $
    it "Should be correct." $ do
      (traceShowId (cond ((3 >< 2) [1.0 ..] :: Matrix Double) 3 0 1 0)) `shouldBe`
        ((3 >< 2) [0, 0, 1, 0, 0, 0] :: Matrix Double)
  describe "m.v behaviour" $
    it "Should be correct." $ do
      let x = (3 >< 2) [1.0 ..] :: Matrix Double
          theta = vector [10.0, 5.0] :: Vector Double
      (x #> theta) `shouldBe` (vector [20.0, 50.0, 80.0] :: Vector Double)

max_stars_repo_name: krisajenkins/sketch

max_issues_repo_name: krisajenkins/sketch

max_forks_repo_name: krisajenkins/sketch


text:
module Main
    ( main
    ) where

-------------------------------------------------------------------------------
import           Criterion
import           Criterion.Main
import           Data.Foldable
-------------------------------------------------------------------------------
import           Statistics.RollingAverage
-------------------------------------------------------------------------------

main :: IO ()
main = defaultMain [
  bgroup "ravgAdd" $ flip map [1, 100, 1000] $ \n ->
    bgroup (show n) [
        bench "Int" $ whnf (avgList n) ([1..fromIntegral n] :: [Int])
      , bench "Integer" $ whnf (avgList n) ([1..fromIntegral n] :: [Integer])
      ]
  , bgroup "ravg" $ flip map [1, 100, 1000] $ \n ->
      bgroup (show n) [
          bench "Int" $ whnf ravg' (avgList n ([1..fromIntegral n] :: [Int]))
        , bench "Integer" $ whnf ravg' (avgList n ([1..fromIntegral n] :: [Integer]))
      ]
  ]

-------------------------------------------------------------------------------
avgList :: Num a => Int -> [a] -> RollingAvg a
avgList lim = foldl' ravgAdd (mkRavg lim)

-------------------------------------------------------------------------------
ravg' :: (Integral a) => RollingAvg a -> Double
ravg' = ravg

max_stars_repo_name: Soostone/rolling-average

max_issues_repo_name: Soostone/rolling-average

max_forks_repo_name: Soostone/rolling-average


text:
--
-- Pool : image pool
--

module Pool (
  PoolClass
, Pool
, getImages
, getImagesRandomly
, nSample
, nClass
, initSamplePool
, initFilePool
, loadImage
) where

import           Control.Exception
import           Control.Monad
import           Data.List
import qualified Data.Map               as Map
import           Data.Maybe
import           Debug.Trace
import           Numeric.LinearAlgebra
import           System.Directory
import qualified System.IO.Strict       as SIO
import qualified System.Random.Mersenne as MT

import CNN.Image

class PoolClass p where
  {-
  getImages
    IN : pool
         batch size
         offset
  -}
  getImages :: p -> Int -> Int -> IO [Trainer]
  getImagesRandomly :: p -> Int -> Int -> IO [Trainer]
  nSample :: p -> Int
  nClass :: p -> Int

type ImageFile = (Int, String)
data Pool =
  MemPool {
    m :: Map.Map Int Trainer
  , c :: Int
  } |
  FilePool {
    imageFiles :: [[ImageFile]]
  , numclass :: Int
  , classvec :: [Class]
  }

instance PoolClass Pool where
  getImages p@(MemPool m _) bt ofs = do
    let
      s = nSample p
      o = ofs `mod` s
      mx0 = o + bt - 1
      mx2 = mx0 - s
      mx1 = if mx2 < 0 then mx0 else s - 1
      im0 = mapMaybe (`Map.lookup` m) [o..mx1]
      im1 = mapMaybe (`Map.lookup` m) [0..mx2]
    return (im0 ++ im1)
  getImages (FilePool flist nc cv) bt ofs = do
    let flist' = selectPerClass flist bt ofs
    ts <- forM flist' $ \(c, fn) -> do
      im <- loadImage fn
      return $ Just (im, cv !! c)
    let trainers = catMaybes ts
    return trainers

  getImagesRandomly p@(MemPool m _) bt s = do
    r0 <- MT.randomIO :: IO Double
    let
      r = if r0 == 1.0 then 0.0 else r0
      s = nSample p
      ofs = floor (r * (fromIntegral s))
    getImages p bt ofs
  getImagesRandomly p@(FilePool flist nc _) bt s = do
    r0 <- MT.randomIO :: IO Double
    let
      r = if r0 == 1.0 then 0.0 else r0
      s = nSample p `div` nc
      ofs = floor (r * (fromIntegral s))
    getImages p bt ofs

  nSample (MemPool m _) = Map.size m
  nSample (FilePool flist _ _) = count flist
    where
      count :: [[ImageFile]] -> Int
      count [] = 0
      count (x:xs) = length x + count xs

  nClass (MemPool _ c)  = c
  nClass (FilePool _ c _) = c

---- IMAGE POOL ON MEMORY ----
------------------------------
{-
initSamplePool
  IN : channel
       image size
       output size (#class)
       probablity
       n samples

  OUT: memory pool
-}

initSamplePool :: Int -> (Int, Int) -> Int -> Double -> Int -> IO Pool
initSamplePool c (sx, sy) o p n = do
  s0 <- forM [0..(n-1)] $ \i -> do
    let
      cl = i `mod` o  -- class of this image

    -- Image data
    s1 <- forM [1..c] $ \j -> do
      w <-forM [0..(sy-1)] $ \y -> do
        let
          p' = if y `div` st == cl then p else 1-p
        forM [1..sx] $ \x -> pixel p'
      return $ fromLists w
    -- Trainer data
    e1 <- forM [0..(o-1)] $ \j ->
      return $ if j == cl then 1.0 else 0.0
    return (s1, fromList e1)
  return (MemPool (Map.fromList $ zip [0..] s0) c)
  where
    st = sy `div` o
    pixel :: Double -> IO Double
    pixel p = do
      v <- MT.randomIO :: IO Double
      let v' = if v < p then 0.5 else 0.0
      return v'

---- IMAGE POOL FROM FILES ----
-------------------------------
{- |
  initFilePool
  IN : pool dir path
       number of class
  OUT:
-}

initFilePool :: String -> Int -> IO Pool
initFilePool path nc = do
  let classv = map (classNumToVec nc) [0..(nc-1)]
  flist <- forM [0..(nc-1)] $ \i -> do
    let dir = path ++ "/" ++ (show i) ++ "/"
    --files <- listDirectory dir
    files <- getDirectoryContents dir
    let files' = map (makeImageFile i dir) $ filter (isSuffixOf ".dat") files
    return files'
  return (FilePool flist nc classv)
  where
    makeImageFile :: Int -> String -> String -> ImageFile
    makeImageFile i dir fn = (i, dir ++ fn)

loadImage :: String -> IO Image
loadImage fn = do
  res <- try $ (SIO.run $ SIO.readFile fn) :: IO (Either SomeException String)
  return $ case res of
    Left  e -> error ("read error:" ++ (show e))
    Right s -> read s :: Image

--
-- PRIVATE FUNCTIONS
--

selectPerClass :: [[ImageFile]] -> Int -> Int -> [ImageFile]
selectPerClass [] _ _ = []
selectPerClass (x:xs) bt ofs = fl ++ (selectPerClass xs bt ofs)
  where
    len = length x
    ofs' = ofs `mod` len
    ed = ofs' + bt - 1
    ed' = if ed > len then (ed - len) `mod` len else ed
    range = if ed' >= ofs'
      then [ofs' .. ed']
      else [0..ed'] ++ [ofs'..(len - 1)]
    fl = map (x !!) range

max_stars_repo_name: eijian/deeplearning

max_issues_repo_name: eijian/deeplearning

max_forks_repo_name: eijian/deeplearning


text:
{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}
module LispValTest (lispValTests) where

import Test.Tasty
import Test.Tasty.HUnit
import Test.Tasty.SmallCheck as SC
import Test.Tasty.QuickCheck as QC
import Test.SmallCheck.Series

import Data.Char (isAscii)
import Data.Complex (Complex(..))
import Data.List (isPrefixOf)
import Data.IORef (newIORef)
import Data.Ratio (numerator, denominator)
import qualified Data.Vector as V
import qualified Data.HashMap.Strict as Map
import Control.Monad (liftM2, (>=>))
import Control.Monad.Trans (lift)
import System.IO (stdout)

import Val
import EvaluationMonad (unsafeEMtoIO, liftIO)
import System.IO.Unsafe (unsafePerformIO) -- be careful! for arbitrary instance
import Types (extractValue) -- these functions can probably be removed entirely.
import Primitives
import Primitives.WriteLib (writeSharedSH, ushowString)
import Types.Port (pStdout)

instance Serial IO Val where
    series = cons1 (Symbol . pack)
                \/ (series >>= lift . unsafeEMtoIO . makeMutableList)
                \/ (makeImmutableList <$> series)
                \/ cons1 Number
                \/ cons1 (IString . pack)
                \/ (series >>= lift . fmap String . newIORef . pack)
                \/ cons1 Char
                \/ cons1 Bool
                \/ cons0 Undefined
                \/ cons0 Nil
                \/ cons0 EOF
                \/ (series >>= lift . fmap Vector . V.thaw . V.fromList)
                -- \/ cons3 Primitive
                -- Primitive function types left off for now
                -- \/ cons5 Closure
                -- Env is probably difficult/impossible to give Serial instance
                -- \/ cons1 Port
                -- Handle is probably difficult/impossible to give Serial instance

instance Serial IO LispErr where
    series = cons2 NumArgs
                \/ cons2 TypeMismatch
                -- Parser ParseError left off for now
                \/ cons1 BadSpecialForm
                \/ cons2 NotFunction
                \/ cons2 (\str name -> UnboundVar str (pack name))
                \/ cons1 (EvaluateDuringInit . pack)
                \/ cons1 Default
                \/ cons0 Quit

instance Monad m => Serial m Arity where
    series = cons1 Exactly
          \/ cons1 AtLeast
          \/ cons2 Between

instance Monad m => Serial m Number where
    series = cons1 Real \/ cons1 Complex

instance Monad m => Serial m RealNumber where
    series = cons1 Bignum \/ cons1 Ratnum \/ cons1 Flonum

instance Arbitrary Val where
    arbitrary = sized lispval'
      where lispval' 0 = oneof simpleCons
            lispval' n = oneof $ simpleCons ++
                            [ do num <- choose (0, n)
                                 unsafePerformIO
                                    . unsafeEMtoIO
                                    . makeMutableList
                                  <$> vectorOf num subval
                            , do num <- choose (0, n)
                                 unsafePerformIO
                                   . fmap Vector . V.thaw . V.fromList
                                  <$> vectorOf num subval
            -- TODO: there's quite a conundrum here. I'm not _super_
            -- concerned about being able to test immutable data. After
            -- all, all the implementations should mirror the mutable versions
            -- exactly anyway. However there's always the chance that I do
            -- something really dumb.
            -- However, it's an invariant that once we start making immutable
            -- data, everything that that data contains must also be immutable.
            -- How to encode into the Arbitrary instance? Perhaps we should
            -- choose mutable/immutable with probability favoring mutable
            -- and then we need a helper that only generates immutable data.
                            ]
              where subval = lispval' $ n `div` 2
            simpleCons = [ Symbol . pack <$> arbitrary
                         , Number <$> arbitrary
                         , IString . pack <$> arbitrary
                         , Char <$> arbitrary
                         , Bool <$> arbitrary
                         ]

instance Arbitrary LispErr where
    arbitrary = oneof [ do n <- choose (0, 3)
                           NumArgs <$> arbitrary <*> vectorOf n arbitrary
                      , liftM2 TypeMismatch arbitrary arbitrary
                      , BadSpecialForm <$> arbitrary
                      , liftM2 NotFunction arbitrary arbitrary
                      , liftM2 UnboundVar arbitrary (pack <$> arbitrary)
                      , EvaluateDuringInit . pack <$> arbitrary
                      , Default <$> arbitrary
                      , return Quit
                      ]

instance Arbitrary Arity where
    arbitrary = oneof
      [ Exactly . QC.getPositive <$> arbitrary
      , AtLeast . QC.getPositive <$> arbitrary
      , mkBetween <$> arbitrary <*> arbitrary
      ]
      where
        mkBetween lo hi = Between (QC.getPositive lo) (QC.getPositive hi)

instance Arbitrary Number where
    arbitrary = oneof [Real <$> arbitrary, Complex <$> arbitrary]

instance Arbitrary RealNumber where
    arbitrary = oneof
      [ Bignum <$> arbitrary
      , Ratnum <$> arbitrary
      , Flonum <$> arbitrary
      ]

lispValTests :: TestTree
lispValTests = testGroup "Val" [unitTests, propTests]

unitTests :: TestTree
unitTests = testGroup "Unit tests"
    [ testShowPort
    , testShowPrimitives
    , testShowFunctions
    ]

propTests :: TestTree
propTests = testGroup "Property tests" [qcTests, scTests]

qcTests :: TestTree
qcTests = testGroup "(QuickCheck)"
    [ prop_TerminationErrors
    , prop_ShowVal
    ]

scTests :: TestTree
scTests = testGroup "(SmallCheck)"
    [ prop_ShowErrPrefix
    ]

-- UNIT TESTS
testShowPort :: TestTree
testShowPort = testCase "Ports show correctly" $
    (show :: Val -> String) (Port pStdout) @?= "#<port>"

testShowPrimitives :: TestTree
testShowPrimitives = testCase "Primitives show correctly" $
    mapM_
        (\(key, func) ->
            let str = case func of
                  Primitive arity _ _ -> "#<procedure "
                                           ++ symbolAsString key ++ " "
                                           ++ displayArity arity
                                           ++ ">"
                  PrimMacro {} -> "#<macro " ++ symbolAsString key ++ ">"

            in show func @?= str)
        -- Test only some of the primitives list as eventually it will be quite large
        . take 50 $ Map.toList primitives
  where
    displayArity (Exactly n) = "(" ++ simple n ++ ")"
    displayArity (Between lo hi) = "(" ++ between lo hi ++ ")"
    displayArity (AtLeast 0) = "_"
    displayArity (AtLeast n) = "(" ++ simple n ++ " . _)"

    simple 0 = ""
    simple 1 = "_"
    simple n = "_ " ++ simple (n - 1)

    between 0 0  = ""
    between 0 1  = "[_]"
    between 0 hi = "[_] " ++ between 0 (hi - 1)
    between n hi = "_ " ++ between (n-1) (hi - 1)

testShowFunctions :: TestTree
testShowFunctions = testCase "Functions show correctly" $
    do let emptyEnv = []
       show (Closure [] Nothing [] emptyEnv (Just "testFunc"))
         @?= "#<procedure testFunc ()>"
       show (Closure ["x"] Nothing [] emptyEnv (Just "testFunc"))
         @?= "#<procedure testFunc (x)>"
       show (Closure [] Nothing [] emptyEnv Nothing)
         @?= "#<procedure lambda ()>"
       show (Closure ["x"] Nothing [] emptyEnv Nothing)
         @?= "#<procedure lambda (x)>"
       show (Closure [] (Just "xs") [] emptyEnv (Just "testFunc"))
         @?= "#<procedure testFunc xs>"
       show (Closure ["x"] (Just "xs") [] emptyEnv (Just "testFunc"))
         @?= "#<procedure testFunc (x . xs)>"
       show (Closure ["x"] (Just "xs") [] emptyEnv Nothing)
         @?= "#<procedure lambda (x . xs)>"
       show (Closure ["x", "y", "z"] (Just "others") [] emptyEnv Nothing)
         @?= "#<procedure lambda (x y z . others)>"

-- PROPERTY TESTS
prop_TerminationErrors :: TestTree
prop_TerminationErrors = QC.testProperty "Only Quit is a termination error" $
    withMaxSuccess 50 $
        \input -> isTerminationError input == case input of
            Quit -> True
            _    -> False

prop_ShowVal :: TestTree
prop_ShowVal = QC.testProperty "Showing a Val produces correct string" $
    \input -> QC.ioProperty $
        (==) <$> unsafeEMtoIO (referenceShowV input) <*> writeSharedSH input
  where referenceShowV = freezeList >=> referenceShow
        -- im not sure that this is right if s has non-ascii characters
        -- TODO [r7rs]: it's explicitly wrong, but so is 'write' atm.
        referenceShow (FNotList (Symbol s)) = pure $ unpack s
        referenceShow (FNotList EOF)        = pure   "#<eof>"
        referenceShow (FNotList (Number n)) = pure $ showNumber n
        referenceShow (FNotList (IString s)) = pure $ ushowString (unpack s) ""
        referenceShow (FNotList (Char c)) = pure $ case c of
            ' '  -> "#\\space"
            '\t' -> "#\\tab"
            '\n' -> "#\\newline"
            '\r' -> "#\\carriage-return"
            _    -> "#\\" ++ [c]
        referenceShow (FNotList (Bool b)) = pure $
            if b then "#t" else "#f"
        referenceShow (FNotList (Vector v)) =
            (\es -> "#(" ++ unwords es ++ ")")
            <$> (liftIO (V.freeze v) >>= mapM referenceShowV . V.toList)

        referenceShow (FList vs) = (\es -> "(" ++ unwords es ++ ")")
            <$> mapM referenceShowV vs
        referenceShow (FDottedList vs v) =
            (\es e -> "(" ++ unwords es ++ ". " ++ e ++ ")")
            <$> mapM referenceShowV vs <*> referenceShowV v
        --    "(" ++ unwords (map show vs) ++ " . " ++ show v ++ ")"

        showNumber (Real (Bignum i)) = show i
        showNumber (Real (Ratnum r))
          | denominator r == 1 = show $ numerator r
          | otherwise = show (numerator r)
                      ++ "/"
                      ++ show (denominator r)
        showNumber (Real (Flonum f))
          | isInfinite f = sign : "inf.0"
          | isNaN f      = "+nan.0"
          | otherwise    = show f
          where sign = if f < 0 then '-' else '+'
        showNumber (Complex (r :+ i))
          | isExactZero $ Real i = showNumber (Real r)
          | isExactZero $ Real r = showNumber (Real i) ++ "i"

          -- minus sign will be given since i is negative
          | i < 0 || isNegativeZero i = showNumber (Real r)
                                      ++ showNumber (Real i)
                                      ++ "i"
          | otherwise = showNumber (Real r)
                      ++ (if isInfinite i || isNaN i then "" else "+")
                      ++ showNumber (Real i)
                      ++ "i"

prop_ShowErrPrefix :: TestTree
prop_ShowErrPrefix = SC.testProperty "Showing LispErr is prefixed with 'Error:'" $
    changeDepth (const 3) $
        \input -> "Error: " `isPrefixOf` show (input :: LispErr)

max_stars_repo_name: JKTKops/Thcheme

max_issues_repo_name: JKTKops/Thcheme

max_forks_repo_name: JKTKops/Thcheme


text:
module School.Train.GradientDescent
( gradientDescent ) where

import Conduit ((.|), ConduitM, mapMC, runConduit, sinkNull, takeWhileC)
import Control.Monad.State.Lazy (get)
import Data.Maybe (isJust)
import Numeric.LinearAlgebra (Matrix)
import School.Train.AppTrain (AppTrain, runAppTrain)
import School.Train.GradientDescentPass (gradientDescentPass)
import School.Train.IterationHandler (IterationHandler(..))
import School.Train.UpdateParams (UpdateParams)
import School.Train.StoppingCondition (StoppingCondition(..))
import School.Train.TrainState (TrainState)
import School.Unit.CostFunction (CostFunction(..))
import School.Unit.Unit (Unit)
import School.Unit.UnitBackward (BackwardStack)

stopping :: StoppingCondition a
         -> ConduitM (BackwardStack a)
                     (Maybe ())
                     (AppTrain a)
                     ()
stopping condition = mapMC $ \_ -> do
  state <- get
  if runCondition condition state
    then return Nothing
    else return (Just ())

gradientDescent :: ConduitM () (Matrix a) (AppTrain a) ()
                -> [Unit a]
                -> CostFunction a (AppTrain a)
                -> UpdateParams a
                -> StoppingCondition a
                -> IterationHandler a
                -> TrainState a
                -> IO (Either String
                              (TrainState a))
gradientDescent source
                units
                cost
                update
                condition
                handler
                initState = do
  let pass = setupCost cost source
          .| gradientDescentPass units cost update
          .| runHandler handler
          .| stopping condition
          .| takeWhileC isJust
          .| sinkNull
  result <- runAppTrain initState . runConduit $ pass
  return $ snd <$> result

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


text:
module Examples where

import Control.Lens hiding (para)
import Prelude hiding (id, (.))
import System.Random
import Control.Monad

import Numeric.LinearAlgebra.Array
import Numeric.LinearAlgebra.Array.Util

import CategoricDefinitions
import Autodiff.GAD
import Autodiff.D
import Ops
import Para
import TrainUtils
import TensorUtils
import OnesLike

{-
f (p, q) a = (p + q*a)
-}
linRegFn :: (Additive a, Num a) => DType ((a, a), a) a
linRegFn = addC . (id `x` mulC) . assocL

-- returns (inp, out) pair we want to learn
sampleData :: IO (Double, Double)
sampleData = do
    x <- randomIO :: IO Double
    return (x, 3 + 7*x)

run :: IO ()
run = do
    initialParams <- randomIO :: IO (Double, Double)
    let initialLearner = Learner initialParams (Para linRegFn) sgd
        sampler = zip3 [0..] (repeat sampleData) (repeat sqDiff)

    finalLearner <- foldM trainStepWithCost initialLearner (take 10000 sampler)
    putStrLn $ "Starting network parameters:\n" ++ show initialParams
    putStrLn $ "Final network parameters:\n" ++ show (finalLearner ^. p)

    return ()

sampleDataTensor :: IO (Tensor, Tensor)
sampleDataTensor = do
    d <- randomTensor [5, 3] "bf"
    return (d, 3 + 7 * d)

-- multiply two arrays you get and sum the "b" axis
l :: _ => ParaDType (Tensor, Tensor) Tensor Tensor
l = Para $ sumAxes "b" . linRegFn

run1 :: IO _
run1 = do
    p1 <- randomTensor [3, 2] "fo"
    p2 <- randomTensor [2] "o"
    let initialLearner = Learner (p1, p2) l sgd
        sampler = zip3 [0..] (repeat sampleDataTensor) (repeat (sumAxes "b" . sqDiff))

    finalLearner <- foldM trainStepWithCost initialLearner (take 1 sampler)
    putStrLn $ "Starting network parameters:\n" ++ arrShow (p1, p2)
    putStrLn $ "Final network parameters:\n" ++ arrShow (finalLearner ^. p)

    return finalLearner


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

module Lib
    ( someFunc
    ) where

import Numeric.GSL.Fitting
import Numeric.LinearAlgebra
import Numeric.AD
import Numeric.AD.Internal.Reverse
import Data.Reflection

someFunc :: IO ()
someFunc = putStrLn "someFunc"

type FitRes = ([(Double, Double)], Matrix Double)
type FitData = [([Double],([Double],Double))]

-- type ManualJacob = ([Double] -> [Double] -> [[Double]])

data JacobType = AutoJacob | ManualJacob ([Double] -> [Double] -> [[Double]])

data FitParams = FitParams { jacob :: JacobType
                           , guess :: [Double]
                           , iter :: Int
                           , absTol :: Double
                           , relTol :: Double
                           }

defaultParams = FitParams { jacob = AutoJacob,
                            guess = repeat 1,
                            iter = 1000,
                            absTol = 1E-20,
                            relTol = 1E-20 }

     -- ([Double] -> [Double] -> [Double]) model
     -- ([Double] -> [Double] -> [[Double]]) model der
-- Only Doubles since hmatrix requires it
fit :: [Double]
     -> [Double]
     -> [Double]
     -> (forall a . Floating a => [a]->[a]->[a])
     -> JacobType
     -> [Double]->  Maybe FitRes

fit xs ys sigma model jacobian guess
    | xeqy && ( seqx || seqone)  = Just $ fitModelScaled 1E-20 1E-20 500 (model, sanitizedJacob jacobian) fitdata guess
    | otherwise  = Nothing
    where
    lx = length xs
    ly = length ys
    ls = length sigma
    xeqy = lx == ly
    seqx = ls == lx
    seqone = ls == 1
    fitdata = unsafeFormatData xs ys sigma
    sanitizedJacob jacob = case jacob of AutoJacob -> mkJac model
                                         ManualJacob f -> f

{-
fit :: [Double]
     -> [Double]
     -> [Double]
     -> (forall a . Floating a => [a]->[a]->[a])
     -> ([Double] -> [Double] -> [[Double]])
     -> [Double]->  Maybe FitRes

fit xs ys sigma model jacobian guess
    | xeqy && ( seqx || seqone)  = Just $ fitModelScaled 1E-20 1E-20 500 (model, mkJac model) fitdata guess
    | otherwise  = Nothing
    where
    lx = length xs
    ly = length ys
    ls = length sigma
    xeqy = lx == ly
    seqx = ls == lx
    seqone = ls == 1
    fitdata = unsafeFormatData xs ys sigma
-}

{-
fit :: [Double]
     -> [Double]
     -> [Double]
     -> ([Double] -> [Double] -> [Double])
     -> ([Double] -> [Double] -> [[Double]])
     -> [Double]->  Maybe FitRes

fit xs ys sigma model jacobian guess
    | xeqy && ( seqx || seqone)  = Just $ fitModelScaled 1E-20 1E-20 500 (model, jacobian) fitdata guess
    | otherwise  = Nothing
    where
    lx = length xs
    ly = length ys
    ls = length sigma
    xeqy = lx == ly
    seqx = ls == lx
    seqone = ls == 1
    fitdata = unsafeFormatData xs ys sigma -}

unsafeFormatData :: [Double] -> [Double] -> [Double] -> FitData
unsafeFormatData xs ys sigma
    |  length sigma == length xs = zip (return <$> xs) $ zip (return <$> ys) sigma
    |  length sigma == 1 = zip (return <$> xs) $ zip (return <$> ys) (repeat (head sigma))

--jac :: (Num a, Traversable f, Functor g) =>
  --     (forall s. Reifies s Tape => f (Reverse s a) -> [Reverse s a] -> g (Reverse s a))
    --   -> f a
      -- -> g (f a)

jac2 :: (Num a, Traversable f1, Traversable f2, Functor g) =>
       (forall s. Reifies s Tape => f1 (Reverse s a) -> f2 (Reverse s a) -> g (Reverse s a))
       -> f1 a
       -> f2 a

       -> g (f1 a)

jac2 f p x = jacobian ((flip f) (fmap auto x)) p

jac3 :: (Floating a, Traversable f1, Traversable f2, Functor g) =>
       (forall s. Reifies s Tape => f1 (Reverse s a) -> f2 (Reverse s a) -> g (Reverse s a))
       -> f1 a
       -> f2 a

       -> g (f1 a)

jac3 f p x = jacobian ((flip f) (fmap auto x)) p

mkJac :: (Num a) => (forall s. Reifies s Tape => [Reverse s a] -> [Reverse s a] -> [Reverse s a])
       -> [a]
       -> [a]

       -> [[a]]

mkJac f p x = jacobian ((flip f) (fmap auto x)) p

-- to fit data
expModel [a,lambda,b] [t] = [a * exp (-lambda * t)+b]

expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]]

-- to generate data
expModel2 [a,lambda,b] t = a * exp (-lambda * t)+b

expModelDer2 [a,lambda,b] t = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]]
-- to go from unbox to box:
-- a = (\p -> \x -> return (expModel2 p (x !! 0)))::Floating a => [a] -> [a] -> [a]

xs = [a | a <- [0..100]]
sigma = [0.1]
ys = (map (expModel2 [5,0.1,1]) xs)

max_stars_repo_name: Magalame/Spinell

max_issues_repo_name: Magalame/Spinell

max_forks_repo_name: Magalame/Spinell


text:
-- Programatica Front-End Commands, level 0
module Pfe0Cmds where
import Prelude hiding (print,putStr,putStrLn,catch,readFile)
import Data.List((\\))
import Control.Monad(unless,when)
import Data.Maybe(fromJust)

import PFE0

import HsName(ModuleName(..),isMainModule,sameModuleName)
--import HsConstants(main_mod)
import HsLexerPass1(rmSpace,Pos(line)) -- for lines of code metrics

import PfeParse
import PrettyPrint
import MUtils
import SimpleGraphs(reverseGraph,reachable)
import DirUtils(optCreateDirectory)
import AbstractIO
import Statistics

pfe0 parseModule args = runPFE0Cmds () pfe0Cmds parseModule args

--runPFE0Cmds ext cmds = runPFE0 ext $ doCmd (cmds, projectStatus)
runPFE0Cmds ext cmds = runCmds (runPFE0 ext) cmds

pfe0Cmds = projman++graphqueries

addHelpCmd prg cmds0 = cmds where cmds=cmds0++helpCmd prg cmds

helpCmd prg cmds =
  [("help",(noArgs (putStrLn (usage prg cmds)),"list available commands"))]

projman =
  -- Project management (a project is a set of source files)
  [("new"     , (qfileArgs  new, "create a new project containing zero or more files")),
   ("add"     , (qfileArgs  add, "add files to the project")),
   ("remove"  , (qfileArgs  remove, "remove files from the project")),
   ("prune"   , (moduleArgs prune,"remove unreachable modules from the project")),
   ("files"   , (noArgs     files,"list files in the project")),
   ("options" , (noArgs     options,"show options in effect"))]

graphqueries =
  -- Module graph queries
  [("modules" , (moduleArgs modules,"show a topologically sorted list of modules")),
   ("graph"   , (graphArgs graph,"show module dependecy (sub)graph")),
--   ("dirgraph", (moduleArgs dirgraph,"show source directory dependecy (sub)graph")),
--   ("dotdirgraph", (moduleArgs dotdirgraph,"dot format source directory dependecy (sub)graph")),
   ("unused"  , (moduleArgs unused,"list unimported and unreachable modules")),
   ("file"    , (moduleArg  file,"which file is the module in")),
   ("module"  , (fileArg    module',"which module does the file contain")),
   ("loc"     , (moduleArgs locModules, "number of lines of code")),
   ("sizemetrics" , (noArgs sizeMetrics, "number of lines per module metrics")),
   ("locmetrics" , (noArgs locMetrics, "number of lines of code per module metrics")),
   ("importmetrics",(noArgs importMetrics, "number of imports per module metrics")),
   ("exportmetrics",(noArgs exportMetrics, "number of exports (importers) per module metrics"))]

--- Project management ---------------------------------------------------------
new    quiet args = do newProject; addPaths quiet args
add    quiet args = addPaths quiet args
remove quiet args = removePaths quiet args

files = putStr . unlines =<< allFiles

options = print =<< parserFlags

qfileArgs f = f #@ kwOption "-quiet" <@ filenames

--- Module graph queries -------------------------------------------------------
modules ms = do g <- getSortedSubGraph (just ms)
	        let sccs = map (map (fst.snd)) g
	        putStrLn "Topologically sorted strongly connected components:"
	        putStr (unlines (map (unwords.map pp) sccs))

graphArgs = moduleArgs' graphOpts
  where
    graphOpts = (,,) #@ kwOption "-rev" <@ kwOption "-dot" <@ kwOption "-dir"

graph (rev,dot,dir) =
    pput . out_conv_dirs @@ getSubGraph . just
  where
    out_conv_dirs  = if dir then out_conv . srcdirGraph  else out_conv . map snd
    out_conv = out . conv
    out   = if dot then dotFormat    else makeFormat
    conv  = if rev then reverseGraph else id

-- Dot format functions contributed by Claus Reinke:

dotFormat g =
   "digraph ModuleGraph {size=\"7.5,10.5\";ratio=fill;" $$
   vcat [q x<>"->"<>braces (fsep [q d<>";"|d<-xs])|(x,xs)<-g] $$
   "}"
  where
    q = doubleQuotes

makeFormat g = vcat [x<>":"<+>fsep xs|(x,xs)<-g]

--dirgraph    ms = pput . makeFormat =<< srcdirGraph ms
--dotdirgraph ms = pput . dotFormat =<< srcdirGraph ms

srcdirGraph g = dg
  where moddirs = [(m,dirname path)|(path,(m,_))<-g]
	moddir m = fromJust (lookup m moddirs)
	dg = mapSnd usort $
	     collectByFst
	       [(moddir m,dir')|(_,(m,ms))<-g,dir'<-usort (map moddir ms)]

dirname = current . reverse . dropWhile (/='/') . reverse
  where current "" = "./"
	current path = path

file m =
    do g <- getCurrentModuleGraph
       pput $ vcat [m<>":"<+>f|(f,(m',_))<-g,m'==m]

module' f =
    do g <- getCurrentModuleGraph
       pput $ vcat [f<>":"<+>m|(f',(m,_))<-g,f'==f]

getUnused ms =
  do fg <- getCurrentModuleGraph
     let g = map snd fg
         unimported = map fst g \\ [m|(m,_:_)<-reverseGraph g]
	 roots = if null ms then mains else ms
	   where mains = [m|(m,_)<-g,isMainModule m]
	 r = reachable g roots
	 unreached = map fst g \\ r
     when (null roots) $
       fail "No root modules given, no Main module in the project"
     return (fg,roots,unimported,unreached)

mainModules =
  do fg <- getCurrentModuleGraph
     return [m|(f,(m,_))<-fg,isMainModule m]

unused ms =
  do (_,roots,us,unr) <- getUnused ms
     if null us -- unlikely...
       then putStrLn "All modules are imported somewhere"
       else pput $ "The following modules are never imported:" $$
	           nest 4 (fsep us)
     unless (null unr) $
       do pput $
	    sep [ppi "The following modules are unreachable from",
		 nest 2 (fsep roots<>":")]
	  pput $ nest 4 (fsep unr)

prune ms =
  do (fg,_,_,unreached) <- getUnused ms
     removePaths False [f|(f,(m,_))<-fg,m `elem` unreached]

sizeMetrics =
    pput . ppStatistics "number of lines" "module"
      =<< mapM (forkM (fst.snd) (size #. readFile . fst))
      =<< getCurrentModuleGraph
  where
    size = length . lines

locMetrics =
    pput . ppStatistics "number of lines of code" "module"
      =<< mapM (forkM (fst.snd) (loc . snd #. lex0SourceFile.fst))
      =<< getCurrentModuleGraph

--forkM f g x = (,) (f x) # g x
forkM f g x =
  do y <- g x
      -- seq to process one file at a time in locMetrics & sizeMetrics
     seq y (return (f x,y))

loc = length . squeezeDups . map (line.fst.snd) . rmSpace

locModules = mapM_ locModule
locModule = locFile @@ findFile
locFile = pput . loc . snd @@ lex0SourceFile

importMetrics = graphMetrics id "number of imports"
exportMetrics = graphMetrics reverseGraph "number of importers"

graphMetrics f lbl =
    do g <- mapSnd length . f . map snd # getCurrentModuleGraph
       pput $ ppStatistics lbl "module" g

max_stars_repo_name: bubba/HaRe

max_issues_repo_name: bubba/HaRe

max_forks_repo_name: bubba/HaRe


text:
{-# LANGUAGE OverloadedStrings, RecordWildCards, NamedFieldPuns #-}
module Numeric.LP where

import           Control.Monad
import           Data.Array.IArray           (Array, (//))
import qualified Data.Array.IArray           as A
import           Data.Maybe
import qualified Data.Text.Lazy              as T
import qualified Data.Text.Lazy.Builder      as TB
import qualified Formatting                  as F
import           Numeric.LinearAlgebra       (I, Matrix, Vector, atIndex, cols,
                                              flatten, ident, minIndex, rows,
                                              subMatrix, subVector, (!), (><),
                                              (?), (|||), (¿))
import qualified Numeric.LinearAlgebra       as LA
import           Numeric.LinearAlgebra.Devel (atM', modifyMatrix, runSTMatrix,
                                              thawMatrix)
import           Text.UnicodeTable

data LP = LP {
    a        :: Matrix Double,
    lhs, rhs :: Vector I}

lp2Table :: Int -> Int -> LP -> Table
lp2Table nPrec nColW LP{..} = mkTable ta cws hs vs
    where rN = rows a
          cN = cols a
          taElem r c | r>0 && c>0 = F.bprint (F.fixed nPrec) $ a `atIndex` (r-1, c-1)
                     | r==0 && c>1 = TB.fromString $ "x"++show (rhs!(c-1))
                     | r==1 && c==0 =  TB.fromString "z"
                     | r>1 && c==0 = TB.fromString $ "x"++show (lhs!(r-1))
                     | otherwise = TB.fromString ""
          ta = A.array ((0,0), (rN, cN)) [((r,c), taElem r c) | r <- [0..rN], c <- [0..cN]]
          cws = A.listArray (0,cN) (replicate (cN+1) nColW)
          hs = konstArray0 (rN+2) NoLine // [(0, SingleLine), (1, DoubleLine), (rN+1, SingleLine)]
          vs = konstArray0 (cN+2) NoLine // [(0, SingleLine), (1, DoubleLine), (cN+1, SingleLine)]

konstArray0 :: Int -> e -> Array Int e
konstArray0 n x = A.listArray (0,n-1) (replicate n x)

instance Show LP where
    show = T.unpack . TB.toLazyText . fromTable . lp2Table 4 15

transfromToEquationalForm :: Matrix Double -> Matrix Double
transfromToEquationalForm a = a ||| (ident (rows a))

choosePivot :: LP -> Maybe (Int, Int)
choosePivot LP {..} =
    let c = subVector 1 (cols a-1) (flatten $ a ? [0])
        b = subVector 1 (cols a-1) (flatten $ a ¿ [0])
        a' = subMatrix (1,1) (rows a-1, cols a-1) a
    in case LA.find (>0) c of
            s:_ ->
                let as = flatten $ a' ¿ [s] -- column s as a vector
                in Just (1 + minIndex (b / as), 1+s)
            _ -> Nothing

-- see <https://en.wikipedia.org/wiki/Simplex_algorithm>
-- Let's solve the following problem:
--      maximize $c^t x$ subject to
--      Ax ≤ b
--      ∀i, xᵢ ≥ 0
--      x,c ∈ ℝⁿ, b ∈ ℝ^p, A ∈ ℝ^pxn

ex1 :: LP
ex1 = LP {a = (3><3)[0,2,-4, 2,-6,1, 8,3,-4],
          lhs = LA.fromList [0,1,4],
          rhs = LA.fromList [0,2,3]
        }

simplexStep :: LP -> Int -> Int -> LP
simplexStep (LP {..}) sr sc =
    let theta = atM' a sr sc -- save the pivot element
        pc = flatten $ a ¿ [sc] -- save the whole pivot column
        a' = runSTMatrix
            (do
                a'm <- thawMatrix a
                -- replace each row -- except the pivot row -- by that linear combination
                -- of itself and the pivot row which makes its pivot column entry zero.
                forM_ [r | r <- [0..(rows a - 1)], r /= sr] $ \r -> do
                    let m = (-1) * (a `atIndex` (r,sc)) / theta
                    forM_ [0..(cols a - 1)] $ \i -> do
                        modifyMatrix a'm r i (\x -> x + m * (a `atIndex` (sr,i)))
                -- divide the pivot row by the negative of the pivot
                forM_ [0..(cols a - 1)] $ \i -> do
                        modifyMatrix a'm sr i (\x -> x/(-theta))
                -- replace the pivot by the negative of its saved reciprocal
                modifyMatrix a'm sr sc (\_ -> 1/theta)
                -- replace the other elements of the pivot column by its saved values divided
                -- by the saved pivot element
                forM_ [r | r <- [0..(rows a - 1)], r /= sr] $ \j -> do
                        modifyMatrix a'm j sc (\_ -> pc!j/theta)
                return a'm)
    in LP a' lhs rhs

main :: IO ()
main = do
    putStrLn (show ex1)
    let (LP {..}) = ex1
    let (sr, sc) = fromJust $ choosePivot ex1
    putStrLn ("Pivot row/column: "++(show sr)++", "++(show sc))
    let theta = atM' a sr sc
    putStrLn ("Pivot element: "++(show theta))

max_stars_repo_name: tkemps/int-prog

max_issues_repo_name: tkemps/int-prog

max_forks_repo_name: tkemps/int-prog