text:
module Quantum where

--import Data.Kind
--import Control.Applicative
--import Data.Typeable
import Prelude hiding ((^))
--import Data.Singletons
--import Numeric.LinearAlgebra hiding (toInt) -- hmatrix library
import GHC.TypeLits

import Prelim hiding (One)
import Types
import Interface
import DeepEmbedding

import LinTrans2 hiding (cnot)

-- Signature
data QuantumSig exp = MkQubit
type Qubit = MkLType MkQubit

class HasMILL exp => HasQuantum exp where
  new :: Bool -> exp '[] Qubit
  meas :: exp γ Qubit -> exp γ (Lower Bool)
  unitary :: KnownType σ => Unitary σ -> exp γ σ -> exp γ σ

controlBy :: (HasQuantum exp,KnownType σ)
          => Unitary σ -> exp γ (Qubit ⊗ σ) -> exp γ (Qubit ⊗ σ)
controlBy u e = unitary (control u) e

-- -- first element is the control
cnot :: HasQuantum exp
     => exp γ (Qubit ⊗ Qubit) -> exp γ (Qubit ⊗ Qubit)
cnot = controlBy PauliX

-------------------------------------------------------
-- Deep embedding -------------------------------------
-------------------------------------------------------

type instance Effect Deep = DensityMonad
data instance LVal Deep Qubit = QId Int

data QuantumExp :: Sig where
  New     :: Bool -> QuantumExp '[] Qubit
  Meas    :: Deep γ Qubit -> QuantumExp γ (Lower Bool)
  Unitary :: KnownType σ => Unitary σ -> Deep γ σ -> QuantumExp γ σ

instance HasQuantum Deep where
  new = Dom . New
  meas = Dom . Meas
  unitary u = Dom . Unitary u

instance Domain QuantumExp where
  evalDomain (New b) _ = QId <$> newM b
  evalDomain (Meas e) ρ = do QId i <- eval e ρ
                             VPut <$> measM i
  evalDomain (Unitary u e) ρ = do v <- eval e ρ
                                  applyU u v
                                  return v

applyU :: forall σ. KnownType σ => Unitary σ -> LVal Deep σ -> DensityMonad ()
applyU u v = applyMatrix (interpU u) (valToQubits v)

class KnownType σ where
  numQubits :: Int
  valToQubits :: LVal Deep σ -> [Int]
instance KnownType Qubit where
  numQubits = 1
  valToQubits (QId i) = [i]
instance KnownType One where
  numQubits = 1
  valToQubits VUnit = []
instance (KnownType σ1,KnownType σ2) => KnownType (σ1 ⊗ σ2) where
  numQubits = numQubits @σ1 + numQubits @σ2
  valToQubits (VPair v1 v2) = valToQubits v1 ++ valToQubits v2
instance KnownType (Lower α) where
  numQubits = 0
  valToQubits (VPut _) = []

instance Show (Unitary σ) where
  show Identity = "I"
  show Hadamard = "H"
  show PauliX   = "X"
  show PauliY   = "Y"
  show PauliZ   = "Z"
  show (R m)    = "R " ++ show m
  show (Alt u0 u1) = "(" ++ show u0 ++ " ⊕ " ++ show u1 ++ ")"
  show (Transpose u) = show u ++ "†"
-- instance Show (QuantumLExp lang g τ) where
--   show (New b)  = "New(" ++ show b ++ ")"
--   show (Meas q) = "Meas(" ++ show q ++ ")"
--   show (Unitary u e) = "Unitary (" ++ show u ++ ") " ++ show e
-- --  show (ControlBy _ e e') = show e ++ "`ControlBy`" ++ show e'

-- Quantum Data

-- Add more?
data Unitary (σ :: LType) where
  Identity  :: KnownType σ => Unitary σ
  Hadamard  :: Unitary Qubit
  PauliX    :: Unitary Qubit -- (NOT)
  PauliY    :: Unitary Qubit
  PauliZ    :: Unitary Qubit
  R         :: Int -> Unitary Qubit
  Alt       :: KnownType σ => Unitary σ -> Unitary σ -> Unitary (Qubit ⊗ σ)
  Transpose :: Unitary σ -> Unitary σ

control :: KnownType σ  => Unitary σ -> Unitary (Qubit ⊗ σ)
control = Alt Identity

{-
type KnownQubits σ = SingI (NumQubits σ)

unitarySing :: forall σ. Unitary σ -> Sing (NumQubits σ)
unitarySing Identity = sing
unitarySing Hadamard = one
unitarySing PauliX   = one
unitarySing PauliY   = one
unitarySing PauliZ   = one
unitarySing (R _)    = one
unitarySing (Alt u0 _)    = SS $ unitarySing u0
unitarySing (Transpose u) = unitarySing u
-}

-- instance HasQuantumEffect ('Sig DensityMonad sigs) where
-- --  type QUnitary ('Sig DensityMonad sigs) = Density
--type QUnitary σ = Squared (NumQubits σ)

interpU :: forall σ. KnownType σ => Unitary σ -> Matrix
interpU Identity = ident $ numQubits @σ
interpU Hadamard = hadamard
interpU PauliX   = pauliX
interpU PauliY   = pauliY
interpU PauliZ   = pauliZ
interpU (R m)    = undefined -- rotation about the z axis by 2πi/2^m?
interpU (Alt (u0 :: Unitary σ') u1) =
    (newD False `kron` interpU u0) + (newD True `kron` interpU u1)
interpU (Transpose u) = transpose $ interpU u

type family   NumQubits (τ :: LType) :: Nat
type instance NumQubits One            = 0
type instance NumQubits Qubit          = 1
type instance NumQubits (τ1 ⊗ τ2) = NumQubits τ1 + NumQubits τ2

{-
instance HasQuantumDom lang => Domain QuantumDom (lang :: Lang sig) where

  evalDomain _ (New b)   = do
    i <- newQubit @sig b
    return $ vqubit i
  evalDomain ρ (Meas e)  = do
    VQubit i <- toDomain @QuantumDom <$> eval' ρ e
    b <- measQubit @sig i
    return $ vput b
  evalDomain ρ (Unitary u e) = do
    v  <- eval' ρ e
    qs <- valToQubits @sig v
    applyU @sig u qs
    return v

type Qubits (τ :: LType sig) = [QId]

valToQubits :: forall sig (lang :: Lang sig) τ.
              HasQuantumDom lang => LVal lang τ -> SigEffect sig (Qubits τ)
valToQubits v = case toDomain' @QuantumDom v of
    Just (VQubit i) -> return [i]
    Nothing -> case toDomain' @OneDom v of
      Just VUnit -> return []
      Nothing -> case toDomain' @TensorDom v of
        Just (VPair v1 v2) -> liftA2 (++) (valToQubits v1) (valToQubits v2)
        Nothing -> return []
--        case fromLVal' proxyLower v of
--          Just (VPut _) -> return []
--          Nothing       -> error "Cannot extract qubits from the given value"

type family   NumQubits (τ :: LType sig) :: Nat
type instance NumQubits ('LType _ 'OneSig)            = 'Z
type instance NumQubits ('LType _ 'QubitSig)          = 'S 'Z
type instance NumQubits ('LType _ ('TensorSig τ1 τ2)) = NumQubits τ1 `Plus` NumQubits τ2
-}

--type instance NumQubits ('LType _ ('PlusSig τ1 τ2))   = NumQubits τ1 `Plus` NumQubits τ2
type instance NumQubits (Lower _)      = 0

-- EXAMPLES

-- Flip -------------------------------------------

qflip :: Lin Deep Bool
qflip = suspend $ meas (unitary Hadamard (new False))

-- ----------------------------------------------------
-- -- Teleportation -----------------------------------
-- ----------------------------------------------------

plus_minus :: HasQuantum exp => Bool -> Lift exp Qubit
plus_minus b = suspend $ unitary Hadamard $ new b

share :: HasQuantum exp => Lift exp (Qubit ⊸ Qubit ⊗ Qubit)
share = suspend . λ $ \q -> cnot (q ⊗ new False)

bell00 :: HasQuantum exp => Lift exp (Qubit ⊗ Qubit)
bell00 = suspend $ force (plus_minus False) `letin` \a ->
                   force share ^ a

alice :: HasQuantum exp => Lift exp (Qubit ⊸ Qubit ⊸ Lower (Bool,Bool))
alice = suspend . λ $ \q -> λ $ \a ->
            cnot (q ⊗ a) `letPair` \(q,a) ->
            meas (unitary Hadamard q) >! \x ->
            meas a >! \y ->
            put (x,y)

bob :: HasQuantum exp => Bool -> Bool -> Lift exp (Qubit ⊸ Qubit)
bob x y = suspend . λ $ \b ->
    if y then unitary PauliX b else b `letin` \b ->
    if x then unitary PauliZ b else b

teleport :: HasQuantum exp => Lift exp (Qubit ⊸ Qubit)
teleport = suspend . λ $ \q ->
    force bell00 `letPair` \(a,b) ->
    force alice ^ q ^ a >! \(x,y) ->
    force (bob x y) ^ b

teleport0 :: Lin Deep Bool
teleport0 = suspend . meas $ force teleport ^ new False

main' :: Lin Deep Bool
main' = suspend $ force bell00 `letPair` \(q1,q2) ->
                 meas q1 >! \x ->
                 meas q2 >! \y ->
                 put (x == y)

main :: IO ()
main = print $ run main'

-- Dependent fourier transform

type family (σ :: LType) ⊗⊗ (n :: Nat) :: LType where
    σ ⊗⊗ n = CompareOrd (CmpNat n 1)
                        One                -- n = 0
                        σ                  -- n = 1
                        (σ ⊗ (σ ⊗⊗ (n-1))) -- n >= 2
--  σ ⊗⊗ Z = One
--  σ ⊗⊗ S Z = σ
--  σ ⊗⊗ (S (S n)) = σ ⊗ (σ ⊗⊗ S n)

{- TODO: fix
rotations :: forall (m :: Nat) (n :: Nat) exp. (HasQuantum exp)
          => Sing m -> Sing n -> Lift exp (Qubit ⊗⊗ S n ⊸ Qubit ⊗⊗ S n)
rotations _ SZ      = idL
rotations _ (SS SZ) = idL
rotations m (SS n'@(SS _)) = suspend . λpair $ \(c,qs) -> qs `letPair` \(q,qs') ->
    force (rotations m n') ^ (c ⊗ qs') `letPair` \(c,qs') ->
    controlBy (R $ 2 + toInt m - toInt n') (c ⊗ q) `letPair` \(c,q) ->
    c ⊗ (q ⊗ qs')

fourier :: forall n exp. (HasQuantum exp)
        => Sing n -> Lift exp (Qubit ⊗⊗ n ⊸ Qubit ⊗⊗ n)
fourier SZ = idL
fourier (SS SZ) = suspend . λ $ unitary Hadamard
fourier (SS n'@(SS _)) = Suspend . λpair $ \(q,w) ->
        force (fourier n') ^ w `letin` \w ->
        force (rotations (SS n') n') ^ (q ⊗ w)
-}

max_stars_repo_name: jpaykin/LNLHaskell

max_issues_repo_name: jpaykin/LNLHaskell

max_forks_repo_name: jpaykin/LNLHaskell


text:
{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE DeriveFunctor       #-}
{-# LANGUAGE ScopedTypeVariables #-}
-- | Chi-Squared test for uniformity.
module Uniformity (testUniformity) where

import Data.List                    (intercalate)
import Data.List                    (foldl')
import Numeric                      (showFFloat)
import Numeric.SpecFunctions        (incompleteGamma)
import Test.Framework.Providers.API (Test, TestName)

import qualified Data.Map as Map

import MiniQC as QC

-- | \( \lim_{n\to\infty} \mathrm{Pr}(V \le v) = \ldots \)
chiDist
    :: Int     -- ^ k, categories
    -> Double  -- ^ v, value
    -> Double
chiDist k x = incompleteGamma (0.5 * v) (0.5 * x) where
  v = fromIntegral (k - 1)

-- | When the distribution is uniform,
--
-- \[
-- \frac{1}{n} \sum_{s = 1}^k \frac{Y_s^2}{p_s} - n
-- \]
--
-- simplifies to
--
-- \[
-- \frac{k}{n} \sum_{s=1}^k Y_s^2 - n
-- \]
--
-- when \(p_s = \frac{1}{k} \), i.e. \(k\) is the number of buckets.
--
calculateV :: Int -> Map.Map k Int -> Double
calculateV k data_ = chiDist k v
  where
    v          = fromIntegral k * fromIntegral sumY2 / fromIntegral n - fromIntegral n
    V2 n sumY2 = foldl' sumF (V2 0 0) (Map.elems data_) where
        sumF (V2 m m2) x = V2 (m + x) (m2 + x * x)

-- Strict pair of 'Int's, used as an accumulator.
data V2 = V2 !Int !Int

countStream :: Ord a => Stream a -> Int -> Map.Map a Int
countStream = go Map.empty where
    go !acc s n
        | n <= 0    = acc
        | otherwise = case s of
            x :> xs -> go (Map.insertWith (+) x 1 acc) xs (pred n)

testUniformityRaw :: forall a. (Ord a, Show a) => Int -> Stream a -> Either String Double
testUniformityRaw k s
    | Map.size m > k = Left $ "Got more elements (" ++ show (Map.size m, take 5 $ Map.keys m) ++ " than expected (" ++ show k ++ ")"
    | p > 0.999999   = Left $
        "Too impropabable p-value: " ++ show p ++ "\n" ++ table
        [ [ show x, showFFloat (Just 3) (fromIntegral y / fromIntegral n :: Double) "" ]
        | (x, y) <- take 20 $ Map.toList m
        ]
    | otherwise      = Right p
  where
    -- each bucket to have roughly 128 elements
    n :: Int
    n = k * 128

    -- buckets from the stream
    m :: Map.Map a Int
    m = countStream s n

    -- calculate chi-squared value
    p :: Double
    p = calculateV k m

testUniformityQC :: (Ord a, Show a) => Int -> Stream a -> QC.Property
testUniformityQC k s = case testUniformityRaw k s of
    Left err -> QC.counterexample err False
    Right _  -> QC.property True

-- | Test that generator produces values uniformly.
--
-- The size is scaled to be at least 20.
--
testUniformity
    :: forall a b. (Ord b, Show b)
    => TestName
    -> QC.Gen a  -- ^ Generator to test
    -> (a -> b)    -- ^ Partitioning function
    -> Int         -- ^ Number of partittions
    -> Test
testUniformity name gen f k = QC.testMiniProperty name
    $ QC.forAllBlind (streamGen gen)
    $ testUniformityQC k . fmap f

-------------------------------------------------------------------------------
-- Infinite stream
-------------------------------------------------------------------------------

data Stream a = a :> Stream a deriving (Functor)
infixr 5 :>

streamGen :: QC.Gen a -> QC.Gen (Stream a)
streamGen g = gs where
    gs = do
        x <- g
        xs <- gs
        return (x :> xs)

-------------------------------------------------------------------------------
-- Table
-------------------------------------------------------------------------------

table :: [[String]] -> String
table cells = unlines rows
  where
    cols      :: Int
    rowWidths :: [Int]
    rows      :: [String]

    (cols, rowWidths, rows) = foldr go (0, repeat 0, []) cells

    go :: [String] -> (Int, [Int], [String]) -> (Int, [Int], [String])
    go xs (c, w, yss) =
        ( max c (length xs)
        , zipWith max w (map length xs ++ repeat 0)
        , intercalate "   " (take cols (zipWith fill xs rowWidths))
          : yss
        )

    fill :: String -> Int -> String
    fill s n = s ++ replicate (n - length s) ' '

ext: hs

max_stars_repo_name: Mikolaj/splitmix

max_issues_repo_name: Mikolaj/splitmix

max_forks_repo_name: Mikolaj/splitmix


text:
{-# LANGUAGE FlexibleContexts, Rank2Types, NoMonomorphismRestriction,
             StandaloneDeriving #-}

module Numeric.MaxEnt.Linear where

import Control.Applicative

import Data.List (transpose)
import qualified Data.Vector.Storable as S

import Numeric.MaxEnt.ConjugateGradient (minimize, dot)
import Numeric.Optimization.Algorithms.HagerZhang05 (Result, Statistics)
import Numeric.AD

multMV :: (Num a) => [[a]] -> [a] -> [a]
multMV mat vec = map (\row -> dot row vec) mat

probs :: (Floating a) => [[a]] -> [a] -> [a]
probs matrix ls = result where
    norm = partitionFunc matrix ls
    result = map (\x -> exp x / norm ) $ (transpose matrix) `multMV` ls

partitionFunc :: (Floating a) => [[a]] -> [a] -> a
partitionFunc matrix ws = sum . map exp . multMV (transpose matrix) $ ws

-- This is almost the sam as the objectiveFunc

objectiveFunc :: (Floating a) => [[a]] -> [a] -> [a] -> a
objectiveFunc as moments ls = log $ partitionFunc as ls - dot ls moments

data LinearConstraints = LC
  { unLC :: forall a. (Floating a) => ([[a]], [a]) }

-- These instances default the underlying numeric type of `LC` to `Double`,
-- which may be problematic for some usages.
deriving instance Eq LinearConstraints
deriving instance Show LinearConstraints

-- | This is for the linear case Ax = b
--   @x@ in this situation is the vector of probablities.
--
--  Consider the 1 dimensional circular convolution using hmatrix.
--
--  >>> import Numeric.LinearAlgebra
--  >>> fromLists [[0.68, 0.22, 0.1], [0.1, 0.68, 0.22], [0.22, 0.1, 0.68]] <> fromLists [[0.2], [0.5], [0.3]]
--  (3><1) [0.276, 0.426, 0.298]
--
--   Now if we were given just the convolution and the output, we can use 'linear' to infer the input.
--
--   >>> linear 3.0e-17 $ LC ([[0.68, 0.22, 0.1], [0.1, 0.68, 0.22], [0.22, 0.1, 0.68]], [0.276, 0.426, 0.298])
--   Right (fromList [0.2000000000000001,0.49999999999999983,0.30000000000000004])
--
--   I fell compelled to point out that we could also just invert the original
--   convolution matrix. Supposedly using maxent can reduce errors from noise if
--   the convolution matrix is not properly estimated.
linear :: Double
       -- ^ Tolerance for the numerical solver
       -> LinearConstraints
       -- ^ The matrix A and column vector b
       -> Either (Result, Statistics) (S.Vector Double)
       -- ^ Either a description of what went wrong or the probability
       --   distribution
linear tolerance constraints  =
    let (matrix, output) = unLC constraints
        obj = objectiveFunc matrix output
        n = length output
    in (S.fromList . probs matrix . S.toList) <$> minimize tolerance n obj

--------------------------------------------------------------------------------
-- I updated everything below to work with the new types, but it's not clear to
-- me what it's for.  -- EP
--------------------------------------------------------------------------------

linear' :: (Floating a, Ord a)
        => LinearConstraints
        -- ^ The matrix A and column vector b
        -> [[a]]
        -- ^ Either a description of what went wrong or the probability
        --   distribution
linear' constraints =
    let (matrix, output) = unLC constraints
        obj = objectiveFunc matrix output
        guess = 1 : replicate (length output - 1) 0
    in map (probs matrix) . gradientDescent obj $ guess

linear'' :: (Floating a, Ord a)
         => LinearConstraints
         -- ^ The matrix A and column vector b
         -> [[a]]
         -- ^ Either a description of what went wrong or the probability
         --   distribution
linear'' constraints =
    let (matrix, output) = unLC constraints
        obj = objectiveFunc matrix output
        guess = 1 : replicate (length output - 1) 0
    in map (probs matrix) . conjugateGradientDescent obj $ guess

--test1 = LC ( [ [0.892532,0.003851,0.063870,0.001593,0.038155]
--             , [0.237713,0.111149,0.326964,0.271535,0.052639]
--             , [0.133708,0.788233,0.051543,0.003976,0.022539]
--             , [0.238064,0.263171,0.112279,0.270452,0.116034]
--             , [0.844155,0.011312,0.001470,0.001826,0.141237]
--             ]
--           ,
--             [0.246323,0.235600,0.071699,0.211339,0.238439]
--           )

max_stars_repo_name: jfischoff/maxent

max_issues_repo_name: jfischoff/maxent

max_forks_repo_name: jfischoff/maxent


text:
{-# LANGUAGE ConstraintKinds     #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TemplateHaskell     #-}
{-# LANGUAGE TypeOperators       #-}

{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
module Test.Grenade.Layers.Internal.Pooling where

import           Grenade.Layers.Internal.Pooling

import           Control.Monad

import           Numeric.LinearAlgebra                  hiding (konst,
                                                         uniformSample, (===))

import           Hedgehog
import qualified Hedgehog.Gen                           as Gen
import qualified Hedgehog.Range                         as Range

import qualified Test.Grenade.Layers.Internal.Reference as Reference
import           Test.Hedgehog.Compat

import           Grenade.Types

prop_poolForwards_poolBackwards_behaves_as_reference =
  let ok extent kernel = [stride | stride <- [1..extent], (extent - kernel) `mod` stride == 0]
      output extent kernel stride = (extent - kernel) `div` stride + 1
  in  property $ do
        height   <- forAll $ choose 2 100
        width    <- forAll $ choose 2 100
        kernel_h <- forAll $ choose 1 (height - 1)
        kernel_w <- forAll $ choose 1 (width - 1)
        stride_h <- forAll $ Gen.element (ok height kernel_h)
        stride_w <- forAll $ Gen.element (ok width kernel_w)
        input    <- forAll $ (height >< width) <$> Gen.list (Range.singleton $ height * width) (Gen.realFloat $ Range.linearFracFrom 0 (-100) 100)

        let outFast       = poolForward 1 height width kernel_h kernel_w stride_h stride_w input
        let retFast       = poolBackward 1 height width kernel_h kernel_w stride_h stride_w input outFast

        let outReference  = Reference.poolForward kernel_h kernel_w stride_h stride_w (output height kernel_h stride_h) (output width kernel_w stride_w) input
        let retReference  = Reference.poolBackward kernel_h kernel_w stride_h stride_w  input outReference

        outFast === outReference
        retFast === retReference

prop_same_pad_pool_behaves_as_reference_when_zero_pad =
  let output extent kernel_dim stride = (extent - kernel_dim) `div` stride + 1
      kernel i s = let x = ceiling ((fromIntegral i :: RealNum) / (fromIntegral s :: RealNum)) in i - (x - 1) * s
  in  property $ do
        height   <- forAll $ choose 2 100
        width    <- forAll $ choose 2 100

        stride_h <- forAll $ choose 1 (height - 1)
        stride_w <- forAll $ choose 1 (width - 1)

        let kernel_h = kernel height stride_h
            kernel_w = kernel width stride_w

        input    <- forAll $ (height >< width) <$> Gen.list (Range.singleton $ height * width) (Gen.realFloat $ Range.linearFracFrom 0 (-100) 100)

        guard $ output height kernel_h stride_h == (ceiling $ (fromIntegral height :: RealNum) / (fromIntegral stride_h :: RealNum))
        guard $ output width kernel_w stride_w  == (ceiling $ (fromIntegral width :: RealNum)  / (fromIntegral stride_w :: RealNum))

        let outFast       = validPadPoolForwards 1 height width kernel_h kernel_w stride_h stride_w 0 0 0 0 input
        let outReference  = poolForward 1 height width kernel_h kernel_w stride_h stride_w input

        assert $ norm_Inf (outFast - outReference) < 0.000001

prop_same_pad_pool_behaves_correctly_at_edges = withTests 1 $ property $ do
  let input           = (2 >< 2) [-0.01, -0.04, -0.02, -0.03]
      expected_output = (2 >< 2) [-0.01, -0.03, -0.02, -0.03]

      out = validPadPoolForwards 1 2 2 2 2 1 1 0 0 1 1 input

  assert $ norm_Inf (out - expected_output) < 0.000001

prop_same_pad_pool_behaves_correctly_at_edges_three_channels = withTests 1 $ property $ do
  let input           = (6 >< 2) [ 0.7, -0.9, -1.4, -0.1, -0.3, 0.5, 0.1, 0.2, -1.1, -0.7, 0.5, -0.7]
      expected_output = (6 >< 2) [ 0.7, -0.1, -0.1, -0.1, 0.5, 0.5, 0.2, 0.2, 0.5, -0.7, 0.5, -0.7]

      out = validPadPoolForwards 3 2 2 2 2 1 1 0 0 1 1 input

  assert $ norm_Inf (out - expected_output) < 0.000001

tests :: IO Bool
tests = checkParallel $$(discover)

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
module Quantum
  ( QRegister,
    QGate,
    C,
    qbitsNum,
    dimensionsNum,
    qbit,
    qregister,
    qbitL,
    qregisterL,
    makeGateC,
    makeGate,
    (|>),
    (|+>),
    measure,
  )
where

import Data.Complex
import Data.List (sort)
import Matrix
import System.Random (randomRIO)
import Utils

type C = Complex Double

data QRegister = QRegister
  { qbitsNum :: Int,
    states :: Matrix C
  }

instance Show QRegister where
  show (QRegister n xm) = show xm

data QGate = QGate
  { dimensionsNum :: Int,
    operator :: Matrix C
  }

instance Show QGate where
  show (QGate n xm) = show xm

class Entanglable a where
  entangle :: a -> a -> a

instance Entanglable QRegister where
  entangle (QRegister n1 xm) (QRegister n2 ym) = QRegister n zm
    where
      n = n1 + n2
      zm = kron xm ym

instance Entanglable QGate where
  entangle (QGate n1 xm) (QGate n2 ym) = QGate n zm
    where
      n = n1 + n2
      zm = kron xm ym

qbitL :: Int -> [C]
qbitL 0 = [1 :+ 0, 0 :+ 0]
qbitL 1 = [0 :+ 0, 1 :+ 0]
qbitL _ = error "Expecting |0> or |1>"

qbit :: Int -> QRegister
qbit = QRegister 1 . vectorV . qbitL

applyGate :: QGate -> QRegister -> QRegister
applyGate (QGate n1 xm) (QRegister n2 ym)
  | n1 /= n2 = error "Dimensions of gate and number of qbits in register disagree"
  | otherwise = QRegister n2 (mult xm ym)

qregister :: [Int] -> QRegister
qregister = foldl1 entangle . map qbit

qregisterL :: [Int] -> [C]
qregisterL = concat . toList . states . qregister

makeGateC :: Matrix C -> QGate
makeGateC xm
  | r /= c = error "Gate matrix must be square"
  | not . powerOfTwo $ r = error "Gate matrix dimensions must be power of two"
  | otherwise = QGate n xm
  where
    r = rowsNum xm
    c = colsNum xm
    n = log2 r
    powerOfTwo n = odd n || (elem n . take (log2 n + 1) . map (2 ^) $ [1 ..])
    log2 = fromIntegral . truncate . logBase 2.0 . fromIntegral

makeGate :: Matrix Double -> QGate
makeGate = makeGateC . fmap (:+ 0.0)

infixl 7 |>

(|>) :: QRegister -> QGate -> QRegister
(|>) = flip applyGate

infixl 8 |+>

(|+>) :: Entanglable a => a -> a -> a
(|+>) = entangle

measure :: QRegister -> IO Int
measure (QRegister n xm) = do
  rnd <- randomRIO (0, 1) :: IO Double
  let state = length . takeWhile (< rnd) $ probs
  return (state - 1)
  where
    probs = sort . (\xs -> 0.0 : xs ++ [1]) . scanl1 (+) . map ((^ 2) . magnitude) . concat . toList $ xm

max_stars_repo_name: meownoid/Deutsch-Jozsa

max_issues_repo_name: meownoid/Deutsch-Jozsa

max_forks_repo_name: meownoid/Deutsch-Jozsa


text:
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables   #-}
module Bio.ChromVAR.Background (getBackgroundPeaks) where

import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Storable as S
import qualified Data.Vector as V
import qualified Data.Matrix.Unboxed as MU
import qualified Data.Vector.Unboxed.Mutable as UM
import Control.Monad
import Data.Ord
import Data.List
import Control.Arrow
import qualified Data.Matrix.Static.Dense as D
import qualified Data.Matrix.Static.Generic as D
import Data.Matrix.Dynamic (Dynamic(..), matrix)

import Statistics.Sample
import System.Random.MWC
import System.Random.MWC.Distributions (categorical)

import Bio.ChromVAR.Utils

type PeakGroup = ( MU.Matrix Double   -- ^ Group to group distance
                 , V.Vector (U.Vector Int)  -- ^ Peaks in each group
                 , U.Vector Int       -- ^ Peak-group membership
                 )

getBackgroundPeaks :: Int -> U.Vector (Double, Double) -> GenIO -> IO [U.Vector Int]
getBackgroundPeaks n xs gen =
    let pg = mkPeakGroup xs
    in replicateM n $ getBackgroundPeak gen pg

getBackgroundPeak :: GenIO -> PeakGroup -> IO (U.Vector Int)
getBackgroundPeak gen (weightMat, bins, membership) = U.generateM n $ \i -> do
    let ws = weightMat `MU.takeRow` (membership U.! i)
    grp <- categorical ws gen
    let peaks = bins V.! grp
    idx <- uniformR (0, U.length peaks - 1) gen
    return $ peaks U.! idx
  where
    n = U.length membership
{-# INLINE getBackgroundPeak #-}

-- | Make peak groups.
mkPeakGroup :: U.Vector (Double, Double) -> PeakGroup
mkPeakGroup raw = (weights, V.fromList groups, membership)
  where
    transformed = case matrix (map (\(x,y) -> [x,y]) $ U.toList raw) of
        Dynamic mat@(D.Matrix _) -> U.fromList $
           map ((\[x,y] -> (x,y)) . S.toList) $ D.toRows $ whiten Cholesky mat
    weights = MU.generate (U.length points, U.length points) $ \(i,j) ->
        weight (points U.! i) (points U.! j)
      where
        points = U.fromList $ flip map groups $ \is -> mean *** mean $
            U.unzip $ U.map (transformed U.!) is
    membership = U.create $ do
        v <- UM.new $ U.length transformed
        forM_ (zip [0..] groups) $ \(x, is) -> U.forM_ is $ \i ->
            UM.unsafeWrite v i x
        return v
    groups = go [] 0 0 $ sortBy (comparing fst) $ zip (U.toList transformed) [0..]
      where
        go acc i j (((x,y), idx) : rest)
            | null acc || i' == i || j' == j = go (idx : acc) i' j' rest
            | i' > i || j' > j = U.fromList acc : go [idx] i' j' rest
            | otherwise = error "Impossible"
          where
            i' = truncate $ (x - x_min) / x_step :: Int
            j' = truncate $ (y - y_min) / y_step :: Int
        go _ _ _ _ = []
    (xs, ys) = U.unzip transformed
    x_min = U.minimum xs
    x_max = U.maximum xs
    x_step = (x_max - x_min) / n
    y_min = U.minimum ys
    y_max = U.maximum ys
    y_step = (y_max - y_min) / n
    n = 50
{-# INLINE mkPeakGroup #-}

max_stars_repo_name: kaizhang/ChromVAR

max_issues_repo_name: kaizhang/ChromVAR

max_forks_repo_name: kaizhang/ChromVAR


text:
module Mandelbrot where

import Data.Complex
import Data.Maybe
import Data.Word
import qualified Data.Array.Repa as R
import Data.Array.Repa hiding (map,zip,foldr)

mandelbrot :: Complex Double -> Complex Double -> Complex Double
mandelbrot c z = z*z + c

isMandelbrot :: Int -> Complex Double -> Maybe Int
isMandelbrot n c = case (n' < n) of
                     True -> Just n'
                     False -> Nothing
  where
    m = iterate (mandelbrot c) 0
    n' = length . take n . takeWhile (\z -> magnitude z < 2) $ m

isMandelbrot2 :: Int -> Complex Double -> Int
isMandelbrot2 n c = case (n' < n) of
                      True ->   n'
                      False -> (-1)
  where
    m = iterate (mandelbrot c) 0
    n' = length. take n . takeWhile (\z -> magnitude z <2) $ m

array :: Double -> Complex Double -> Int -> Array D DIM2 (Complex Double)
array size centre n = fromFunction (Z:.n:.n) posToVal
  where
    dx = size / (fromIntegral n)
    posToVal (Z:.i:.j) = (((fromIntegral i)*dx - size/2) :+ ((fromIntegral j)*dx - size/2)) + centre

genMandelbrot :: Int -> Int -> Double -> Complex Double -> Array D DIM2 (Word8,Word8,Word8)
genMandelbrot res n size centre = cs
  where
    zs = array size centre res
    cs = R.map (toColor2 . isMandelbrot2 n) zs
    --cs = R.map toColor ms

toColor2 :: Int -> (Word8,Word8,Word8)
toColor2 n = case (n==(-1)) of
               True  -> (0,0,0)
               False -> (255 - (fromIntegral n),fromIntegral n,50)

toColor :: Maybe Int -> (Word8,Word8,Word8)
toColor (Just n) = (255-(fromIntegral n),(fromIntegral n),(fromIntegral n))
toColor Nothing  = (0,0,0)

toBw2 :: Int -> (Word8, Word8, Word8)
toBw2 n = case (n==(-1)) of
            True  -> (0,0,0)
            False ->(255-n', 255-n', 255-n')
  where n' = fromIntegral n

toBw :: Maybe Int -> (Word8,Word8,Word8)
toBw (Just n') = let n = fromIntegral n' in (255-n,255-n,255-n)
toBw Nothing = (0,0,0)

max_stars_repo_name: ahgibbons/HaskMandelbrot

max_issues_repo_name: ahgibbons/HaskMandelbrot

max_forks_repo_name: ahgibbons/HaskMandelbrot


text:
{-# LANGUAGE DefaultSignatures #-}
{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE UndecidableInstances #-}

{-| This module defines `Distributive` for a fixed `Monad`
    instead of an arbitrary `Functor`.
-}

module Data.Distributive.Monadic (
  -- * Type class
    DistributiveM(..)

  -- ** distributem and collectm in terms of each other
  , distributemByCollectm
  , collectmByDistributem

  -- *  Utilities defined using `DistributiveM`
  , wrapFoldM

  -- ** ConstT
  , ConstT(..)

  -- *  MonadTrans instances
  , distributemTrans
  , collectmTrans

  -- *  Generics
  , gdistributem
  , gcollectm
  , GDistributiveM(..)
  ) where

import Control.Comonad (Comonad(..))
import Control.Comonad.Cofree (Cofree(..))
import Control.Comonad.Trans.Cofree (CofreeT(..))
import qualified Control.Comonad.Trans.Cofree as Trans (CofreeF(..))
import Control.Comonad.Trans.Coiter (CoiterT(..))
import Control.Comonad.Trans.Env (EnvT(..), ask, lowerEnvT)
import Control.Comonad.Trans.Store (StoreT(..), pos)
import Control.Comonad.Trans.Traced (TracedT(..))
import Control.Foldl (FoldM(..))
import Control.Foldl.Utils (extractM, innerJoinFoldM, stepFoldM, unwrapFoldM)
import Control.Monad (join, liftM2)
import Control.Monad.Free (Free(..))
import Control.Monad.Free.Church (F(..), fromF, toF)
import Control.Monad.Trans.Accum (AccumT(..))
import Control.Monad.Trans.Class (MonadTrans(..))
import Control.Monad.Trans.Const
import Control.Monad.Trans.Cont (ContT(..))
import Control.Monad.Trans.Except (ExceptT(..))
import Control.Monad.Trans.Free (FreeT(..))
import qualified Control.Monad.Trans.Free as Trans (FreeF(..))
import Control.Monad.Trans.Free.Church (FT(..), fromFT, toFT)
import Control.Monad.Trans.Identity (IdentityT(..))
import Control.Monad.Trans.Iter (IterT(..))
import Control.Monad.Trans.Maybe (MaybeT(..))
import Control.Monad.Trans.RWS (RWST(..))
import qualified Control.Monad.Trans.RWS.Strict as Strict (RWST(..))
import Control.Monad.Trans.Reader (ReaderT(..))
import Control.Monad.Trans.State (StateT(..))
import qualified Control.Monad.Trans.State.Strict as Strict (StateT(..))
import Control.Monad.Trans.Writer (WriterT(..))
import qualified Control.Monad.Trans.Writer.Strict as Strict (WriterT(..))
import Data.Bifunctor (Bifunctor(..))
import Data.Complex (Complex)
import Data.Distributive (Distributive(..))
import Data.Functor.Compose (Compose(..))
import Data.Functor.Identity (Identity(..))
import Data.Proxy (Proxy(..))
import Data.Semigroup (Dual, Product, Sum)
import GHC.Generics
  ( (:*:)(..)
  , (:.:)(..)
  , Generic1(..)
  , M1(..)
  , Par1(..)
  , Rec1(..)
  , U1(..)
  )

-- | This is the categorical dual of `Traversable`, for a `Monad` fixed by
-- our choice of @f@.
--
-- If you define one of `distributem`, `collectm` you may define the other
-- using `collectmByDistributem`, `distributemByCollectm`, respectively.
--
-- All `MonadTrans`formers that result in a `Monad` are instances of
-- `DistributiveM`, see `collectmTrans`.
class Monad m => DistributiveM m f | f -> m where
  -- | The dual of `sequence`, for a fixed `Monad`
  --
  -- @
  --  `distributem` = `collectm` `id`
  -- @
  --
  distributem :: m (f a) -> f (m a)
  default distributem :: (Generic1 f, GDistributiveM m (Rep1 f)) => m (f a) -> f (m a)
  distributem = gdistributem

  -- | Map then `distributem`:
  --
  -- @
  --  `collectm` f = `distributem` . `fmap` f
  -- @
  --
  collectm :: (a -> f b) -> m a -> f (m b)
  default collectm :: (Generic1 f, GDistributiveM m (Rep1 f)) => (a -> f b) -> m a -> f (m b)
  collectm = gcollectm

instance Monad m => DistributiveM m (FoldM m a) where
  distributem f = FoldM (fmap unwrapFoldM . flip stepFoldM) f (return . extractM)

  collectm = collectmByDistributem

-- | If there were no functional dependencies, this would be the "other"
-- instance for `Identity`:
--
-- @
--  instance `Monad` m => `DistributiveM` m `Identity` where
--    `distributem` :: m (`Identity` a) -> `Identity` (m a)
--    `distributem` = `Identity` . `fmap` `runIdentity`
-- @
--
-- But it'd violate the functional dependency since `Identity` does not
-- determine @m@.
--
-- By ignoring the chosen `Monad`, we can recreate the instance.
--
-- See `wrappedCofreeFToProd` for an example.
instance Monad m => DistributiveM m (ConstT m) where
  distributem = ConstT . fmap runConstT
  collectm f = ConstT . fmap (runConstT . f)

-- distributem and collectm in terms of each other

-- | `distributem` in terms of `collectm`
distributemByCollectm :: DistributiveM m f => m (f a) -> f (m a)
distributemByCollectm = collectm id

-- | `collectm` in terms of `distributem`
collectmByDistributem :: DistributiveM m f => (a -> f b) -> m a -> f (m b)
collectmByDistributem f = distributem . fmap f

-- Utilities defined using `DistributiveM`

-- | Implement `wrap` using `distributem` and `innerJoinFoldM`
wrapFoldM :: Monad m => m (FoldM m a b) -> FoldM m a b
wrapFoldM = innerJoinFoldM . distributem

-- GHC.Generics instances

instance (DistributiveM m f, DistributiveM m g) =>
         DistributiveM m (f :*: g)

instance (DistributiveM m f, DistributiveM m g, Functor f) =>
         DistributiveM m (f :.: g)

instance (DistributiveM m f, DistributiveM m g, Functor f) =>
         DistributiveM m (Compose f g)

instance DistributiveM Identity Par1

instance DistributiveM m f => DistributiveM m (M1 i c f)

instance DistributiveM m f => DistributiveM m (Rec1 f)

instance DistributiveM Identity U1

-- Lifted instances

instance DistributiveM Identity Identity

instance DistributiveM Identity Complex

instance DistributiveM Identity Dual

instance DistributiveM Identity Sum

instance DistributiveM Identity Product

instance DistributiveM Identity Proxy

-- MonadTrans instances

instance Monad m => DistributiveM m (ReaderT a m) where
  distributem = distributemTrans -- ReaderT . distribute . fmap runReaderT
  collectm = collectmTrans -- f xs = ReaderT $ \x -> flip (runReaderT . f) x <$> xs

instance (Monoid w, Monad m) => DistributiveM m (AccumT w m) where
  distributem = distributemTrans -- AccumT . fmap (>>= fmap (first return)) . collect runAccumT
  collectm = collectmTrans -- f xs = AccumT $ \x -> xs >>= fmap (first return) . flip (runAccumT . f) x

instance Monad m => DistributiveM m (ContT r m) where
  distributem = distributemTrans -- xs = ContT $ \f -> xs >>= ($ f . return) . runContT
  collectm = collectmTrans -- f xs = ContT $ \g -> xs >>= \y -> (runContT . f) y (g . return)

instance Monad m => DistributiveM m (ExceptT e m) where
  distributem = distributemTrans -- ExceptT . (>>= fmap (fmap return) . runExceptT)
  collectm = collectmTrans -- f xs = ExceptT $ xs >>= fmap (fmap return) . runExceptT . f

instance Monad m => DistributiveM m (IdentityT m) where
  distributem = distributemTrans
  collectm = collectmTrans

instance Monad m => DistributiveM m (MaybeT m) where
  distributem = distributemTrans
  collectm = collectmTrans

instance (Monoid w, Monad m) => DistributiveM m (RWST r w s m) where
  distributem = distributemTrans
  collectm = collectmTrans

instance (Monoid w, Monad m) => DistributiveM m (Strict.RWST r w s m) where
  distributem = distributemTrans
  collectm = collectmTrans

instance Monad m => DistributiveM m (StateT s m) where
  distributem = distributemTrans
  collectm = collectmTrans

instance Monad m => DistributiveM m (Strict.StateT s m) where
  distributem = distributemTrans
  collectm = collectmTrans

instance (Monoid w, Monad m) => DistributiveM m (WriterT w m) where
  distributem = distributemTrans
  collectm = collectmTrans

instance (Monoid w, Monad m) => DistributiveM m (Strict.WriterT w m) where
  distributem = distributemTrans
  collectm = collectmTrans

-- | `distributem` for any `MonadTrans` that's a `Monad` when its base functor is.
distributemTrans :: (Monad m, MonadTrans t, Monad (t m)) => m (t m a) -> t m (m a)
distributemTrans = fmap return . join . lift

-- | `collectm` for any `MonadTrans` that's a `Monad` when its base functor is.
--
-- We `lift` the `m` then @(`>>=`)@ the function to the result.
-- Finally, we `return` inside of the transformer.
collectmTrans :: (Monad m, MonadTrans t, Monad (t m)) => (a -> t m b) -> m a -> t m (m b)
collectmTrans f = fmap return . (>>= f) . lift

-- Comonad transformers

instance (Comonad m, Monad m) => DistributiveM m (EnvT e m) where
  distributem = liftM2 EnvT (ask . extract) (fmap lowerEnvT)
  collectm = collectmByDistributem

instance (Comonad m, Monad m) => DistributiveM m (StoreT s m) where
  distributem = liftM2 StoreT (fmap (\(~(StoreT acc _)) -> distribute acc)) (pos . extract)
  collectm = collectmByDistributem

instance (Monoid w, Monad m) => DistributiveM m (TracedT w m) where
  distributem = TracedT . fmap (distribute . runTracedT)
  collectm f = TracedT . (>>= runTracedT . fmap return . f)

-- Free

instance (Comonad m, Monad m) => DistributiveM m (Cofree m) where
  distributem xs =
    fmap extract xs :< (xs >>= (\(~(_ :< ys)) -> fmap return <$> ys))
  collectm = collectmByDistributem

-- | A wrapped representation of `CofreeF` without the outermost
-- `Comonad`ic layer
newtype WrappedCofreeF f w a = WrappedCofreeF
  { runWrappedCofreeF :: Trans.CofreeF f a (CofreeT f w a)
  }

instance (Functor f, Functor w) => Functor (WrappedCofreeF f w) where
  fmap f = WrappedCofreeF . bimap f (fmap f) . runWrappedCofreeF

-- | Convert `CofreeT` to our wrapped representation
toWrappedCofreeF :: Functor w => CofreeT f w a -> w (WrappedCofreeF f w a)
toWrappedCofreeF = fmap WrappedCofreeF . runCofreeT

-- | Convert `WrappedCofreeF` to a product, which has a suitable automatically
-- derived `DistributiveM` instance.
{-# INLINE wrappedCofreeFToProd #-}
wrappedCofreeFToProd ::
     WrappedCofreeF f w a -> (ConstT f :*: Compose (IdentityT f) (CofreeT f w)) a
wrappedCofreeFToProd ~(WrappedCofreeF (x Trans.:< xs)) =
  ConstT x :*: Compose (IdentityT xs)

-- | Convert `WrappedCofreeF` from a product
{-# INLINE wrappedCofreeFFromProd #-}
wrappedCofreeFFromProd ::
     (ConstT f :*: Compose (IdentityT f) (CofreeT f w)) a -> WrappedCofreeF f w a
wrappedCofreeFFromProd ~(ConstT x :*: Compose (IdentityT xs)) =
  WrappedCofreeF $ x Trans.:< xs

instance (DistributiveM m f, Functor f) =>
         DistributiveM m (WrappedCofreeF m f) where
  distributem = wrappedCofreeFFromProd . collectm wrappedCofreeFToProd
  collectm f = wrappedCofreeFFromProd . collectm (wrappedCofreeFToProd . f)

instance (DistributiveM m f, Functor f) => DistributiveM m (CofreeT m f) where
  distributem =
    CofreeT .
    fmap runWrappedCofreeF . getCompose . collectm (Compose . toWrappedCofreeF)
  collectm f =
    CofreeT .
    fmap runWrappedCofreeF .
    getCompose . collectm (Compose . toWrappedCofreeF . f)

instance Monad m => DistributiveM m (CoiterT m) where
  distributem = CoiterT . (>>= fmap (bimap return (fmap return)) . runCoiterT)
  collectm f =
    CoiterT . (>>= fmap (bimap return (fmap return)) . runCoiterT . f)

-- | Deconstruct a `Free`
{-# INLINE free #-}
free :: (a -> b) -> (f (Free f a) -> b) -> Free f a -> b
free f _ (Pure x) = f x
free _ g ~(Free x) = g x

instance Monad m => DistributiveM m (Free m) where
  distributem = Free . fmap (fmap return)
  collectm f =
    Free . (>>= free (return . return . return) (fmap return <$>) . f)

instance Monad m => DistributiveM m (F m) where
  distributem = toF . collectm fromF
  collectm = collectmByDistributem

-- | A wrapped representation of `FreeF` without the outermose `Monad`ic layer
newtype WrappedFreeF f m a = WrappedFreeF
  { runWrappedFreeF :: Trans.FreeF f a (FreeT f m a)
  }

instance (Functor f, Monad m) => Functor (WrappedFreeF f m) where
  fmap f = WrappedFreeF . bimap f (fmap f) . runWrappedFreeF

-- | Convert `FreeT` to our wrapped representation
toWrappedFreeF :: Functor m => FreeT f m a -> m (WrappedFreeF f m a)
toWrappedFreeF = fmap WrappedFreeF . runFreeT

-- | Deconstruct a `FreeF`
{-# INLINE freeF #-}
freeF :: (a -> c) -> (f b -> c) -> Trans.FreeF f a b -> c
freeF f _ (Trans.Pure x) = f x
freeF _ g ~(Trans.Free x) = g x

instance (Monad m, Monad n) => DistributiveM m (WrappedFreeF m n) where
  distributem =
    WrappedFreeF .
    Trans.Free .
    (>>= freeF (return . return . return) (fmap return <$>) . runWrappedFreeF)
  collectm f =
    WrappedFreeF .
    Trans.Free .
    (>>= freeF (return . return . return) (fmap return <$>) .
         runWrappedFreeF . f)

instance (DistributiveM m n, Monad n) => DistributiveM m (FreeT m n) where
  distributem =
    FreeT .
    fmap runWrappedFreeF . getCompose . collectm (Compose . toWrappedFreeF)
  collectm f =
    FreeT .
    fmap runWrappedFreeF . getCompose . collectm (Compose . toWrappedFreeF . f)

instance (DistributiveM m n, Monad n) => DistributiveM m (FT m n) where
  distributem = toFT . collectm fromFT
  collectm = collectmByDistributem

instance Monad m => DistributiveM m (IterT m) where
  distributem = IterT . (>>= fmap (bimap return (fmap return)) . runIterT)
  collectm f = IterT . (>>= fmap (bimap return (fmap return)) . runIterT . f)

-- Generics

-- | A default implementation of `distributem` for `Generic1` types
gdistributem :: (GDistributiveM m (Rep1 f), Generic1 f) => m (f a) -> f (m a)
gdistributem = to1 . gcollectm' from1

-- | A default implementation of `collectm` for `Generic1` types
gcollectm :: (GDistributiveM m (Rep1 f), Generic1 f) => (a -> f b) -> m a -> f (m b)
gcollectm f = to1 . gcollectm' (from1 . f)

-- | The typeclass implementing `gdistributem` (as `gdistributem'`)
-- for `Generic1` types
class Monad m => GDistributiveM (m :: * -> *) (f :: * -> *) | f -> m where
  gdistributem' :: m (f a) -> f (m a)
  gdistributem' = gcollectm' id

  gcollectm' :: (a -> f b) -> m a -> f (m b)

instance GDistributiveM m f => GDistributiveM m (M1 i c f) where
  gcollectm' f = M1 . gcollectm' (unM1 . f)

instance GDistributiveM Identity U1 where
  gcollectm' _ _ = U1

instance GDistributiveM Identity Par1 where
  gcollectm' f = Par1 . fmap (unPar1 . f)

instance (GDistributiveM m f, GDistributiveM m g) =>
         GDistributiveM m (f :*: g) where
  gcollectm' f xs = gcollectm' leftF xs :*: gcollectm' rightF xs
    where
      leftF = (\(~(ys :*: _)) -> ys) . f
      rightF = (\(~(_ :*: zs)) -> zs) . f

instance DistributiveM m f => GDistributiveM m (Rec1 f) where
  gcollectm' f xs = Rec1 $ collectm (unRec1 . f) xs

instance (DistributiveM m f, GDistributiveM m g, Functor f) =>
         GDistributiveM m (f :.: g) where
  gcollectm' f xs = Comp1 $ gdistributem' <$> collectm (unComp1 . f) xs

max_stars_repo_name: michaeljklein/algebraic-foldl

max_issues_repo_name: michaeljklein/algebraic-foldl

max_forks_repo_name: michaeljklein/algebraic-foldl


text:
{-|
Module: MachineLearning.LeastSquaresModel
Description: Least Squares Model
Copyright: (c) Alexander Ignatyev, 2016
License: BSD-3
Stability: experimental
Portability: POSIX

-}

module MachineLearning.LeastSquaresModel
(
  LeastSquaresModel(..)
)

where

import qualified Numeric.LinearAlgebra as LA
import Numeric.LinearAlgebra((<>), (#>), (<.>))
import qualified Numeric.LinearAlgebra.Data as LAD
import qualified Data.Vector.Storable as V

import qualified MachineLearning.Regularization as R

import MachineLearning.Model

data LeastSquaresModel = LeastSquares

instance Model LeastSquaresModel where
  hypothesis LeastSquares x theta = x #> theta

  cost LeastSquares lambda x y theta =
    let m = x #> theta - y
        nExamples = fromIntegral $ LA.rows x
        regTerm = R.costReg lambda theta
    in (LA.sumElements (m * m) * 0.5 + regTerm) / nExamples

  gradient LeastSquares lambda x y theta = ((LA.tr x) #> (x #> theta - y) + regTerm) / nExamples
    where nExamples = fromIntegral $ LAD.rows x
          regTerm = R.gradientReg lambda theta

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
-- | This module provides the core functionality of the ansigraph package:
--   terminal-based graphing for vectors and matrices of real and complex numbers.
--
--   This is implemented via a 'Graphable' type class.
--
--   __Ansigraph is intended to be used in on of two ways:__
--
--   * __By importing "System.Console.Ansigraph"__.
--   This provides all the functionality we typically use, including the FlexibleInstances
--   extension which makes it easier to use graphing functions by allowing instances like
--   'Graphable [Double]'.
--
--
--   * __By directly importing "System.Console.Ansigraph.Core"__, which does not activate
--   FlexibleInstances but includes everything else provided by the other module. This just means
--   you must use one of a few newtype wrappers, namely: 'Graph', 'PosGraph', 'CGraph',
--   'Mat', 'CMat'. They are also available from the standard module.
module System.Console.Ansigraph.Core (

  -- * Core Functionality
  -- ** The Graphable class
    Graphable (..)
  , graph
  , animateWith
  , animate
  , transientAnim
  , transientAnimWith

  -- *** Graphing options
  , GraphSettings (..)

  -- **** Default options
  , graphDefaults
  , blue, pink, white, red, green
  , noColoring

  -- *** ANSI data
  -- **** Basic types from ANSI package
  , Color (..)
  , ColorIntensity (..)

-- **** Custom composite data types
  , AnsiColor (..)
  , Coloring (..)

-- *** ANSI helpers
  , mkColoring
  , fromFG
  , fromBG
  , realColors
  , imagColors
  , colorSets
  , invert
  , interpAnsiColor
  , setColor
  , clear
  , clearLn
  , applyColoring
  , colorStr
  , colorStrLn
  , boldStr
  , boldStrLn

  -- * Graphable wrapper types
  , Graph (..)
  , CGraph (..)
  , PosGraph (..)
  , Mat (..)
  , CMat (..)

  -- * Graphing
  -- *** Horizontal vector graphing (IO actions)
  , displayPV
  , displayRV
  , displayCV

  -- *** Horizontal rendering logic (producing strings)
  , renderPV
  , renderRV
  , renderCV

  -- *** Matrix graphing
  , displayMat
  , displayCMat
  , matShow

-- *** Simple (non-ANSI) graphing for strictly-positive data
  , posGraph
  , posAnim

-- *** For clearing
  , clearBack

) where

import System.Console.Ansigraph.Internal.Core
import System.Console.Ansigraph.Internal.Horizontal
import System.Console.Ansigraph.Internal.Matrix

import System.Console.ANSI
import Control.Concurrent     (threadDelay)
import Control.Monad          (replicateM_)
import Data.Complex           (Complex)
import Control.Monad.IO.Class (MonadIO, liftIO)

-- | Things that ansigraph knows how to render at the terminal are instances of this class.
--
--   In general, when ANSI codes are involved, a 'graphWith' method should fush stdout when
--   finished, and whenever codes are invoked to i.e. change terminal colors. This is easily
--   handled by defining it in terms of 'colorStr' and 'colorStrLn'.
--
--   The 'graphHeight' function specifies how many vertical lines a graph occupies and is
--   needed for animations to work properly
class Graphable a where

  -- | Render a graph to standard output.
  graphWith :: MonadIO m => GraphSettings -> a -> m ()

  -- | The number of vertical lines a graph occupies.
  graphHeight :: a -> Int

-- | Invokes the 'Graphable' type class method 'graphWith' with the
--   default 'GraphSettings' record, 'graphDefaults'.
graph :: MonadIO m => Graphable a => a -> m ()
graph = graphWith graphDefaults

---- IO / ANSI helpers ----

-- | Clear the last @n@ lines of terminal text. Used to make graph animations. Rexported as
--   a handy convenience for other uses.
clearBack :: MonadIO m => Int -> m ()
clearBack n = do
  putStr' "\r"  -- return cursor to horizontal position 0
  replicateM_ n (liftIO $ cursorUpLine 1 >> clearLine)

-- | For some number of frames per second, return the corresponding time delta in microseconds.
deltaFromFPS :: Int -> Int
deltaFromFPS fps = 1000000 `div` fps

---- Animation ----

clearGraph :: MonadIO m => Graphable a => a -> m ()
clearGraph = clearBack . graphHeight

animationFrame :: MonadIO m => Graphable a => GraphSettings -> a -> m ()
animationFrame s x = do
  graphWith s x
  liftIO . threadDelay . deltaFromFPS . framerate $ s
  clearGraph x

-- | Any list of a 'Graphable' type can be made into an animation, by
--   'graph'ing each element with a time delay and screen-clear after each.
--   'GraphSettings' are used to determine the time delta and any coloring/scaling options.
animateWith :: MonadIO m => Graphable a => GraphSettings -> [a] -> m ()
animateWith _ []       = return ()
animateWith s [x]      = graphWith s x
animateWith s (x:y:zs) = animationFrame s x >> animateWith s (y:zs)

-- | Perform 'animateWith' using default options. Equivalent to 'graph'ing each member
--   of the supplied list with a short delay and screen-clear after each.
animate :: MonadIO m => Graphable a => [a] -> m ()
animate = animateWith graphDefaults

-- | Like 'animateWith', only it does not leave the final frame of the animation visible.
transientAnimWith :: MonadIO m => Graphable a => GraphSettings -> [a] -> m ()
transientAnimWith = mapM_ . animationFrame

-- | Like 'animate', only it does not leave the final frame of the animation visible.
transientAnim :: (MonadIO m, Graphable a) => [a] -> m ()
transientAnim = transientAnimWith graphDefaults

---- Wrappers to avoid needing FlexibleInstances ----

-- | Wrapper type for graph of a real vector/function.
newtype Graph = Graph { unGraph :: [Double] }

-- | Wrapper type for graph of a complex vector/function.
newtype CGraph = CGraph { unCGraph :: [Complex Double] }

-- | Wrapper type for graph of a non-negative real vector/function.
newtype PosGraph = PosGraph { unPosGraph :: [Double] }

-- | Wrapper type for graph of a real two-index vector/two-argument function.
newtype Mat = Mat { unMat :: [[Double]] }

-- | Wrapper type for graph of a complex two-index vector/two-argument function.
newtype CMat = CMat { unCMat :: [[Complex Double]] }

instance Graphable Graph where
  graphWith s = displayRV s . unGraph
  graphHeight _ = 2

instance Graphable CGraph where
  graphWith s = displayCV s . unCGraph
  graphHeight _ = 4

instance Graphable PosGraph where
  graphWith s = displayPV s . unPosGraph
  graphHeight _ = 1

instance Graphable Mat where
  graphWith s = displayMat s . unMat
  graphHeight = length . unMat

instance Graphable CMat where
  graphWith s = displayCMat s . unCMat
  graphHeight = length . unCMat

---- helpers for graphing/animating strictly-positive real functions ----

-- | Display a graph of the supplied (non-negative) real vector.
posGraph :: MonadIO m => [Double] -> m ()
posGraph = graph . PosGraph

-- | Display an animation of the supplied list of (non-negative) real vectors.
posAnim :: MonadIO m => [[Double]] -> m ()
posAnim = animate . map PosGraph

max_stars_repo_name: fieldstrength/ansigraph

max_issues_repo_name: fieldstrength/ansigraph

max_forks_repo_name: fieldstrength/ansigraph


text:
{-# OPTIONS_GHC  -fno-warn-unused-binds -fno-warn-unused-matches -fno-warn-name-shadowing -fno-warn-missing-signatures #-}
{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses, UndecidableInstances, FlexibleContexts, TypeSynonymInstances #-}

---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------
-- |
-- | Module : Fequency Moments
-- | Creator: Xiao Ling
-- | Created: 12/17/2015
-- |
---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------

module FeqMoments where

import Prelude hiding (replicate)
import Control.Monad.Random.Class
import Control.Monad.Random
import Control.Monad.State

import Data.Conduit
import Data.Foldable (toList)
import qualified Data.Sequence as S
import qualified Data.Conduit.List as Cl
import Data.Sequence (Seq,(|>),update,empty)

import Core
import Statistics

{-----------------------------------------------------------------------------
  Types
------------------------------------------------------------------------------}

{-----------------------------------------------------------------------------
  Approximate Median
------------------------------------------------------------------------------}


text:
module Polymorphic where

import Data.Complex
import Data.Ratio
import Prelude hiding (product)

stringLength :: String -> Int
stringLength x = case x of
  [] -> 0
  ['a',_] -> 1
  _:xs -> 1 + stringLength xs

integerListLength :: [Integer] -> Int
integerListLength n = case n of
  [] -> 0
  _:xs -> 1 + integerListLength xs

-- | A polymorphic length function
polymorphicLength :: [a] -> Int
polymorphicLength list = case list of
  []     -> 0
  _ : xs -> 1 + polymorphicLength xs

-- | A polymorphic reverse function
reverseOf :: [a] -> [a]
reverseOf n = case n of
  [] -> []
  x:xs -> reverseOf xs ++ [x]

-- | A polymorphic isPalindrome function
isPalindrome :: (Eq a) => [a] -> Bool
isPalindrome list = list == reverseOf list

-- | A polymorphic list equality function
listEqual :: (Eq a) => [a] -> [a] -> Bool
listEqual list1 list2 = case (list1,list2) of
  ([],[] )-> True
  ([],_) -> False
  (_,[]) -> False
  (x:xs,y:ys) -> (x==y) && (xs `listEqual` ys)

isMonotonicallyIncreasing :: (Ord a) => [a] -> Bool
isMonotonicallyIncreasing list = case list of
  [] -> True
  i1:i2:is
    |i1<=i2 -> isMonotonicallyIncreasing (i2:is)
    |otherwise -> False
  [_] -> True

normaliseVector :: (Floating a) => [a] -> [a]
normaliseVector vector = divideByScalar vector (norm vector)

divideByScalar :: (Fractional a) =>[a] -> a -> [a]
divideByScalar vector' scalar = case vector' of
   [] -> []
   f: fs -> (f / scalar): divideByScalar fs scalar

norm :: (Floating a) => [a] -> a
norm vector' = sqrt(sumSqr vector')
  where
  sumSqr :: (Num a)=> [a] -> a
  sumSqr vector'' = case vector'' of
    []    -> 0
    f: fs -> f*f + sumSqr fs

reallyPolymorphicLength :: (Integral b) => [a] -> b
reallyPolymorphicLength list = case list of
  []    -> 0
  _: xs -> 1 + reallyPolymorphicLength xs

printLargestAsString :: (Ord a, Show a) => [a] -> String
printLargestAsString list =
  "The largest element in the list is " ++ show (maximum list)

rationalZero :: Ratio Integer
floatZero    :: Float
doubleZero   :: Double
complexZero  :: Complex Double

rationalZero = 0
floatZero    = 0.0
doubleZero   = 0.0
complexZero  = mkPolar 0.0 0.0

product :: (Num a, Eq a) => a -> a -> a
product x y = case y of
  0->0
  1->x
  _->x+product x(y-1)


text:
{-
 -  <DeclList> ::= { (<FunDecl> | <ValDecl>) }
 -
 -  <FunDecl> ::= 'fun' <Ident> <Ident> '=' <Exp>
 -
 -  <ValDecl> ::= 'val' <Ident> '=' <Exp>
 -
 -  <Expr> ::= <Comp>
 -           | 'if' <Expr> 'then' <Expr> 'else' <Expr>
 -           | 'let' <DeclList> 'in' <Expr> 'end'
 -           | 'fn' <Ident> '=>' <Expr>
 -
 -  <Comp> ::= <Plus> { ('>' | '=' | '<') <Plus> }
 -
 -  <Plus> ::= <Mult> { ('+' | '-') <Mult> }
 -
 -  <Mult> ::= <App>  { ('*' | '/') <App> }
 -
 -  <App>  ::= <Fact> { <Fact> }
 -
 -  <Fact> ::= '(' <Expr> ')'
 -          |  <Integer>
 -          |  <Identifier>
 -}

 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
 {-# LANGUAGE RankNTypes #-}

 {-# LANGUAGE FlexibleContexts #-}

 module Hw8 () where

 -- import Hw1
 -- import Hw2
 -- import Hw3
 -- import Hw4
 -- import Hw5
 -- import Hw6
 -- import Hw7
 import Data.Time.Clock
 import Data.Complex

 import Text.Parsec
 import Text.Parsec.Expr
 import Text.Parsec.Prim (unexpected)
 import Text.Parsec.Language (emptyDef)
 import Data.Either
 import qualified Text.Parsec.Token as Token

 -- Abstract Syntax Tree ----------------------------------------------------------------------------------

 -- function/variable declaration
 data Decl = Fun String String Exp -- fun f x = e;
           | Val String Exp        -- val x = e;

 -- expressions
 data Exp = Lt Exp Exp     -- e1 < e2
         | Gt Exp Exp     -- e1 > e2
         | Eq Exp Exp     -- e1 = e2
         | Plus Exp Exp   -- e1 + e2
         | Minus Exp Exp  -- e1 - e2
         | Times Exp Exp  -- e1 * e2
         | Div Exp Exp    -- e1 div e2
         | Var String     -- x
         | If Exp Exp Exp -- if e0 then e1 else e2
         | Fn String Exp  -- fn x => e
         | Let [Decl] Exp -- let val x = e0; fun f = e1; in e2 end
         | App Exp Exp    -- e1 e2
         | Const Integer  -- n

 instance Show Decl where
   show (Fun f x e) = "fun " ++ f ++ " " ++ x ++ " = " ++ show e
   show (Val x e) = "val " ++ x ++ " = " ++ show e

 instance Show Exp where
   show (Const x) = show x
   show (Plus t1 t2) = showBop t1 "+" t2
   show (Times t1 t2) = showBop t1 "*" t2
   show (Minus t1 t2) = showBop t1 "-" t2
   show (Div t1 t2) = showBop t1 "/" t2
   show (Lt t1 t2) = showBop t1 "<" t2
   show (Gt t1 t2) = showBop t1 ">" t2
   show (Eq t1 t2) = showBop t1 "=" t2
   show (If t0 t1 t2) = "if " ++ show t0 ++ " then " ++ show t1 ++ " else " ++ show t2
   show (Var s) = s
   show (App t1 t2) = "(" ++ show t1 ++ " " ++ show t2 ++ ")"
   show (Fn x t) = "(fn " ++ x ++ " => " ++ show t ++ ")"
   show (Let decls e) = "(let " ++ show decls ++ " in " ++ show e ++ " end)"

 showBop t1 op t2 = "(" ++ show t1 ++ " " ++ op ++ " " ++ show t2 ++ ")"

 -- Parser definitions ----------------------------------------------------------------------------------

 type Parser = Parsec [Char] ()

 -- Runner

 -- print AST or error
 run :: Show a => Parser a -> String -> IO()
 run p input = putStrLn $ case parse p "" input of
       Left error -> show error
       Right x    -> show x

 -- return Either AST or parse error
 runParser :: Parser a -> String -> Either ParseError a
 runParser p input = parse p "" input

 -- Lexer

 lexeme p = do { x <- p; skipMany space <?> ""; return x }  -- get rid of trailing spaces

 symbol name = lexeme $ string name                         -- parse string and then remove spaces

 lexer :: Token.TokenParser ()
 lexer = Token.makeTokenParser style
   where style = emptyDef {
           Token.reservedOpNames = ["+","*","-","=", ">", "<", "in"],
           Token.reservedNames = ["let", "in", "end"]
           }

 parens :: Parser a -> Parser a
 parens = Token.parens lexer

 ident :: Parser String                                     -- raw identifier
 ident = (lexeme $ many1 letter) <?> "identifier"

 integer :: Parser Integer                                  -- parse integer and return numeric value
 integer = (read <$> (lexeme $ many1 digit)) <?> "integer"

 keyword :: String -> Parser ()                             -- make keyword parser with backtrack
 keyword name = (try $ symbol name) >> return ()

 if_ = keyword "if"                                         -- keyword parsers
 then_ = keyword "then"
 else_ = keyword "else"
 fn_ = keyword "fn"
 fun_ = keyword "fun"
 let_ = keyword "let"
 in_ = keyword "in"
 end_ = keyword "end"
 val_ = keyword "val"
 assign_ = keyword "="
 gt_ = keyword ">"
 lt_ = keyword "<"
 eq_ = keyword "="
 pl_ = keyword "+"
 mn_ = keyword "-"
 tm_ = keyword "*"
 dv_ = keyword "/"
 fatArrow_ = keyword "=>"
 sCol_ = keyword ";"

 identifier :: Parser String                                -- identifier that won't collide with keywords
 identifier = try $ do
     name <- ident
     if name `elem` reserved
       then unexpected ("reserved word " ++ name) -- throw parse error
       else return name
   where
     reserved = ["if", "then", "else", "fn", "fun", "let", "in", "end", "val"]

 -- Parser for the (phrase) grammar

 prog :: Parser [Decl]
 prog = skipMany space >> declList                              -- skip leading spaces

 declList :: Parser [Decl]                                      -- parse a list of function or variable declarations
 declList = many $ valDecl <|> funDecl

 funDecl :: Parser Decl
 funDecl = do
   fun_
   name <- identifier
   param <- identifier
   assign_
   t1 <- parseExpression
   return (Fun name param t1)

 valDecl :: Parser Decl
 valDecl = do
   val_
   name <- identifier
   assign_
   t1 <- parseExpression
   return (Val name t1)

 parseExpression :: Parser Exp
 parseExpression =  parseIf <|> parseLet <|> expr

 parseApp :: Parser Exp
 parseApp = ((factor `chainl1` (return App)) <?> "parseApp expr")

 expr = expr' `chainl1` compOp
 expr' = term `chainl1` addOp
 term    = parseApp `chainl1` mulOp
 factor  =  parseVar <|> parseInt <|> parseFn <|> parens expr

 compOp = do {gt_ ; return gtExpr}
         <|> do {lt_ ; return ltExpr}
         <|> do {eq_ ; return eqExpr}

 mulOp   =   do{ tm_ ; return tmExpr   }
         <|> do{ dv_ ; return dvExpr }

 addOp   =   do{ pl_ ; return plExpr }
         <|> do{ mn_ ; return mlExpr }

 gtExpr :: Exp -> Exp -> Exp
 gtExpr = Gt

 ltExpr :: Exp -> Exp -> Exp
 ltExpr = Lt

 eqExpr :: Exp -> Exp -> Exp
 eqExpr = Eq

 tmExpr :: Exp -> Exp -> Exp
 tmExpr = Times

 dvExpr :: Exp -> Exp -> Exp
 dvExpr = Div

 plExpr :: Exp -> Exp -> Exp
 plExpr = Plus

 mlExpr :: Exp -> Exp -> Exp
 mlExpr = Minus

 parseIf :: Parser Exp
 parseIf = do
   if_
   cond <- parseExpression
   then_
   true <- parseExpression
   else_
   false <- parseExpression
   return (If cond true false)

 parseFn :: Parser Exp
 parseFn = do
   fn_
   name <- identifier
   fatArrow_
   exp <- parseExpression
   return (Fn name exp)

 -- parseLet :: Parser Exp
 -- parseLet = do
 --   let_
 --   decl <- parseExpression
 --   in_
 --   exp <- parseExpression
 --   end_
 --   return (Let decl exp)

 parseVar :: Parser Exp
 parseVar = Var `fmap` identifier

 parseInt :: Parser Exp
 parseInt = Const `fmap` integer

 parseLet :: Parser Exp
 parseLet = do
   let_
   decls <- declList
   in_
   exp <- parseExpression
   end_
   return (Let decls exp)
 -- Test code ----------------------------------------------------------------------------------

--  main :: IO ()
--  main = do

--          let fact = "fun fact x = if x < 1 then 1 else x * fact (x-1)\n"
--          let it = "val it = let val y = 10 val z = 15 in fact y end"
--          let fnFn = "fun f x = fn z => y"

--          let d = (fact ++ it)
--          let test =  "val x = x * 5 / 5 - 5\n"
--          let test2 =  "fun fact x = 5\n"
--          let test3 = "if x > 1 then 1 else -1"
--          let test4 = "fun fact x = y"
--          let bla = (test ++ test2 ++ test3)
--          -- run prog bla
--          run prog d

 -- main = do
     -- Hw2
     -- let n = 2^5
     -- let s1 = map (\x -> sin(20*pi*x) + sin(40*pi*x)/2) $ range 0 1 n
     -- -- print $ split [1..64]
     -- -- print(rd 3 s1)
     -- start <- getCurrentTime
     -- let dft1 = map (\x -> x/n) $ absolute $ dft s1
     -- print(rd 2 dft1)
     -- end <- getCurrentTime
     -- print (diffUTCTime end start)
     -- start2 <- getCurrentTime
     -- let fft1 = map (\x -> x/n) $ absolute $ fft s1
     -- print(rd 2 fft1)
     -- end2 <- getCurrentTime
     -- print (diffUTCTime end2 start2)
     -- let v1 = Vec [1.0,2.0,3.0]
     -- let v2 = Vec [2,3,4]
     -- let v3 = Vec [-10,0,10]

     -- Hw3
     -- print $ v1 == v2
     -- print $ show v3
     -- print $ v1 + v2
     -- print $ v1 - v2
     -- print $ v1 * v2
     -- print $ v1 / v3
     -- print $ negate v1
     -- print $ signum v3
     -- print $ abs v3
     -- print $ foldr (*) 1 v1
     -- print $ sin $ v1 * (pi / 2)
     -- print $ v1 + (pure' $ sqrt 2)
     -- print $ realV v1
     -- print $ pure' 'a'

     -- let v1 = Vec [1.123,2.213,3.321]
     -- let v2 = Vec [2,3,4]
     -- let v3 = Vec [-10,0,10]
     -- print $ v1 + v2
     -- print $ v1 - v2
     -- print $ v1 * v2
     -- print $ v1 / v2
     -- print $ negate v1
     -- print $ signum v3
     -- print $ abs v3
     -- print $ v1 + 10
     -- print $ v2 + 1.2
     -- print $ v1 + (pure' $ sqrt 2)
     -- print $ realV v1
     -- print $ imagV v1
     -- print $ realV v1 + imagV v2
     -- print $ sin $ v1 * (pi / 2)
     -- print $ sum v1

     -- Hw4
     -- let n = fromIntegral 2^8
     -- let s1 = realV $ Vec $ fmap (\x -> sin(20*pi*x) + sin(40*pi*x)/2) $ range 0 1 n
     -- print $ fmap (\z -> conjugate z) s1
     -- print s1
     -- start <- getCurrentTime
     -- let dft1 = fmap (/n) $ absolute $ dft s1
     -- print(rd 2 dft1)
     -- let n = fromIntegral 2^8
     -- let s1 = fmap (\x -> sin(20*pi*x) + sin(40*pi*x)/2) $ Vec $ range 0 1 n
     -- print(rd 3 s1)
     -- print s1
     -- let dft1 = fmap (/n) $ absolute $ dft s1
     -- print(rd 2 dft1)
     -- print(rd 2 $ fmap (/n) $ absolute $ idft $ dft s1)
     -- print(rd 2 $ fmap (/n) $ absolute $ ifft $fft s1)

     -- end <- getCurrentTime
     -- print (diffUTCTime end start)
     -- start2 <- getCurrentTime
     -- let fft1 = fmap (/n) $ absolute $ fft s1
     -- print(rd 2 fft1)
     -- end2 <- getCurrentTime
     -- print (diffUTCTime end2 start2)

     -- let v1 = Vec [1,2,3]
     -- print $ v1 + v2
     -- print $ rd 2 v1
     -- let n = fromIntegral 2^8
     -- let s1 = fmap (\x -> sin(20*pi*x) + sin(40*pi*x)/2) $ Vec $ range 0 1 n
     -- print(rd 3 s1)
     -- print(rd 3 $ fmap (\(r:+_) -> r) $ low_pass' 15 $ realV s1)
     -- print(rd 3 $ fmap (\(r:+_) -> r) $ low_pass 15 $ realV s1)
     -- let a = [2,4,1,11,9]
     -- let b = [3,1,5,0,2]
     -- print $ maxList a
     -- print $ maxList' (<) a
     -- print $ rd 1 $ makeNoise 12345 10 0.0 1.0
     -- let n' = 2^8
     -- let noise = makeNoise 0 n' (-1.0) 1.0
     -- let n = fromIntegral n'
     -- let s = fmap (\x -> sin(20*pi*x) + sin(40*pi*x)/2) $ Vec $ range 0 1 n
     -- let s1 = s + noise
     -- -- noisy signal
     -- print(rd 2 s1)
     -- -- noisy frequency
     -- print (rd 2 $ fmap (/n) $ absolute $ dft $ realV s1)
     -- -- filtered signal (dft)
     -- print (rd 3 $ fmap (\(r:+_) -> r) $ low_pass' 15 $ realV s1)
     -- -- filtered signal (fft)
     -- print (rd 2 $ fmap (\(r:+_) -> r) $ low_pass 15 $ realV s1)

max_stars_repo_name: NikolaPeevski/Haskell-stuff

max_issues_repo_name: NikolaPeevski/Haskell-stuff

max_forks_repo_name: NikolaPeevski/Haskell-stuff


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
-- |
-- Module    : Statistics.Distribution.Uniform
-- Copyright : (c) 2011 Aleksey Khudyakov
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Variate distributed uniformly in the interval.
module Statistics.Distribution.Uniform
    (
      UniformDistribution
    -- * Constructors
    , uniformDistr
    , uniformDistrE
    -- ** Accessors
    , uniformA
    , uniformB
    ) where

import Control.Applicative
import Data.Data           (Data, Typeable)
import GHC.Generics        (Generic)
import qualified System.Random.MWC       as MWC

import qualified Statistics.Distribution as D
import Statistics.Internal

-- | Uniform distribution from A to B
data UniformDistribution = UniformDistribution {
      uniformA :: {-# UNPACK #-} !Double -- ^ Low boundary of distribution
    , uniformB :: {-# UNPACK #-} !Double -- ^ Upper boundary of distribution
    } deriving (Eq, Typeable, Data, Generic)

instance Show UniformDistribution where
  showsPrec i (UniformDistribution a b) = defaultShow2 "uniformDistr" a b i
instance Read UniformDistribution where
  readPrec = defaultReadPrecM2 "uniformDistr" uniformDistrE

-- | Create uniform distribution.
uniformDistr :: Double -> Double -> UniformDistribution
uniformDistr a b = maybe (error errMsg) id $ uniformDistrE a b

-- | Create uniform distribution.
uniformDistrE :: Double -> Double -> Maybe UniformDistribution
uniformDistrE a b
  | b < a     = Just $ UniformDistribution b a
  | a < b     = Just $ UniformDistribution a b
  | otherwise = Nothing
-- NOTE: failure is in default branch to guard against NaNs.

errMsg :: String
errMsg = "Statistics.Distribution.Uniform.uniform: wrong parameters"

instance D.Distribution UniformDistribution where
  cumulative (UniformDistribution a b) x
    | x < a     = 0
    | x > b     = 1
    | otherwise = (x - a) / (b - a)

instance D.ContDistr UniformDistribution where
  density (UniformDistribution a b) x
    | x < a     = 0
    | x > b     = 0
    | otherwise = 1 / (b - a)
  quantile (UniformDistribution a b) p
    | p >= 0 && p <= 1 = a + (b - a) * p
    | otherwise        =
      error $ "Statistics.Distribution.Uniform.quantile: p must be in [0,1] range. Got: "++show p
  complQuantile (UniformDistribution a b) p
    | p >= 0 && p <= 1 = b + (a - b) * p
    | otherwise        =
      error $ "Statistics.Distribution.Uniform.complQuantile: p must be in [0,1] range. Got: "++show p

instance D.Mean UniformDistribution where
  mean (UniformDistribution a b) = 0.5 * (a + b)

instance D.Variance UniformDistribution where
  -- NOTE: 1/sqrt 12 is not constant folded (#4101) so it's written as
  --       numerical constant. (Also FIXME!)
  stdDev   (UniformDistribution a b) = 0.2886751345948129 * (b - a)
  variance (UniformDistribution a b) = d * d / 12 where d = b - a

instance D.MaybeMean UniformDistribution where
    maybeMean = Just . D.mean

instance D.MaybeVariance UniformDistribution where
    maybeStdDev   = Just . D.stdDev

instance D.Entropy UniformDistribution where
  entropy (UniformDistribution a b) = log (b - a)

instance D.MaybeEntropy UniformDistribution where
  maybeEntropy = Just . D.entropy

instance D.ContGen UniformDistribution where
    genContVar (UniformDistribution a b) gen = MWC.uniformR (a,b) gen

max_stars_repo_name: vmchale/statistics

max_issues_repo_name: vmchale/statistics

max_forks_repo_name: vmchale/statistics


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

module Parser where

import           Control.Monad.Except
import           Data.Array
import           Data.Char
import           Data.Complex
import           Data.Ratio
import           Lib
import           Numeric
import           System.Environment
import           Text.Printf
import           Data.IORef
import           Text.ParserCombinators.Parsec
                                         hiding ( spaces )
import           LispTypes

trapError
  :: forall  a (m :: * -> *) . (MonadError a m, Show a) => m String -> m String
trapError action = action `catchError` (return . show)

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

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

spaces :: Parser ()
spaces = skipMany1 space

parseString :: Parser LispVal
parseString = do
  char '"'
  x <- many (escapedChars <|> noneOf "\"")
  char '"'
  return $ String x

parseBool :: Parser LispVal
parseBool = do
  char '#'
  x <- oneOf "tf"
  return $ case x of
    't' -> Bool True
    'f' -> Bool False

escapedChars :: Parser Char
escapedChars = do
  char '\\'
  x <- oneOf "\\\"nrt"
  return $ case x of
    '\\' -> x
    '\"' -> '"'
    'n'  -> '\n'
    'r'  -> '\r'
    't'  -> '\t'

parseAtom :: Parser LispVal
parseAtom = do
  first <- letter <|> symbol
  rest  <- many (letter <|> symbol <|> digit)
  return $ Atom (first : rest)

parseChar :: Parser LispVal
parseChar = do
  string "#\\"
  value <- many1 alphaNum
  if length value == 1
    then return $ Character (head value)
    else case (map toLower value) of
      "space"   -> return $ Character ' '
      "newline" -> return $ Character '\n'
      _         -> fail "Invalid character literal"

parseDecimal1 :: Parser LispVal
parseDecimal1 = many1 digit >>= return . Number . read

parseDecimal2 :: Parser LispVal
parseDecimal2 = do
  try $ string "#d"
  parseDecimal1

parseOct :: Parser LispVal
parseOct = do
  try $ string "#o"
  x <- many1 octDigit
  return $ Number (oct2dig x)
  where oct2dig y = fst $ readOct y !! 0

parseHex :: Parser LispVal
parseHex = do
  try $ string "#h"
  x <- many1 hexDigit
  return $ Number (hex2dig x)
  where hex2dig y = fst $ readHex y !! 0

parseBin :: Parser LispVal
parseBin = do
  try $ string "#b"
  x <- many1 digit
  return $ Number (bin2dig x)
 where
  bin2dig = bin2dig' 0
  bin2dig' n "" = n
  bin2dig' n (x : xs) =
    let current = (n * 2) + (if x == '1' then 1 else 0) in bin2dig' current xs

parseInteger :: Parser LispVal
parseInteger =
  parseDecimal1 <|> parseDecimal2 <|> parseOct <|> parseHex <|> parseBin

parseFloat :: Parser LispVal
parseFloat = do
  ds <- many1 digit
  char '.'
  decs <- many1 digit
  return $ Float $ toFloat (ds ++ "." ++ decs)
  where toFloat xs = fst $ head $ readFloat xs

parseRatio :: Parser LispVal
parseRatio = do
  numerator <- many1 digit
  char '/'
  denominator <- many1 digit
  return $ Ratio (read numerator % read denominator)

toDouble :: LispVal -> Double
toDouble (Float  x) = x
toDouble (Number x) = fromIntegral x

parseComplex :: Parser LispVal
parseComplex = do
  r <- try parseFloat <|> parseInteger
  char '+'
  i <- try parseFloat <|> parseInteger
  char 'i'
  return $ Complex (toDouble r :+ toDouble i)

parseNumber :: Parser LispVal
parseNumber =
  try parseComplex <|> try parseFloat <|> try parseRatio <|> parseInteger

parseList :: Parser LispVal
parseList = do
  values <- sepBy parseExpr spaces
  return $ List values

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

parseLists :: Parser LispVal
parseLists = do
  char '('
  value <- try parseList <|> parseDottedList
  char ')'
  return value

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

parseQuasiquote :: Parser LispVal
parseQuasiquote = do
  char '`'
  exp <- parseExpr
  return $ List [Atom "quasiquote", exp]

parseUnquote :: Parser LispVal
parseUnquote = do
  char ','
  exp <- parseExpr
  return $ List [Atom "unquote", exp]

parseVector :: Parser LispVal
parseVector = do
  string "#("
  values <- sepBy parseExpr spaces
  char ')'
  return $ Vector (listArray (0, (length values) - 1) values)

parseQuotes :: Parser LispVal
parseQuotes = parseQuoted <|> parseQuasiquote <|> parseUnquote

parseExpr :: Parser LispVal
parseExpr =
  parseAtom
    <|> parseString
    <|> try parseBool
    <|> try parseVector
    <|> parseNumber
    <|> parseQuotes
    <|> parseLists

readOrThrow :: Parser a -> String -> ThrowsError a
readOrThrow parser input = case parse parser "lisp" input of
  Left err -> throwError $ Parser err
  Right val -> return val

readExpr :: String -> ThrowsError LispVal
readExpr = readOrThrow parseExpr

readExprList = readOrThrow (endBy parseExpr spaces)

parseIt :: String -> Maybe LispVal
parseIt s = case parse parseExpr "lisp" s of
  Left  _ -> Nothing
  Right s -> Just s


text:
{-# LANGUAGE NoMonomorphismRestriction, ImpredicativeTypes ,RankNTypes, TypeOperators #-}
import Kalman
import qualified Data.Vector as V
import Vectorization.Instances
import Local
import Linear.V1
import Data.Distributive
import qualified Data.Foldable as F
import Linear.V2
import Data.Complex
import CartesianProduct
import System.Random.MWC.Distributions
import Vectorization
import Foreign.Storable
import Foreign.Cholesky
import Space.Class
import Control.Applicative
import Control.Monad
import System.Random.MWC
import Graphics.Gnuplot.Simple
import Statistics.Autocorrelation
import Statistics.Sample
import Data.Functor.Compose
import Foreign.FStorable
import Foreign.FStorable

type Model = V1 :|: V1 :|: V1 :|: V1

predSin :: V1 Double -> Model Double -> Model Double
predSin  dt s@(theta :|: omega :|: amp :|: accel)   = s |+|  ((omega*dt )  :|: 0 :|: (dt*accel ) :|: 0 )

projSin :: Model Double -> V1 Double
projSin (theta :|: omega :|: amp :|: accel) = amp * (fmap sin theta )

nx :: Model (Model Double)
nx = diag $ fmap (*1e-3)((1/3) :|: 1 :|: 0.1 :|: 0.1)

getMean :: Functor f => (forall b . g b -> V1 b) -> f (g b , g (g b)) -> f b
getMean s = fmap ( unV1 . s . fst )

getCov :: (Num b ,Functor f )=> (forall b .  g b -> V1 b) -> f (g b , g (g b)) -> f b
getCov s = fmap ( abs .unV1 . s .  unV1 . s . snd )

sigmaBound :: Num b => (forall b . g b -> V1 b) -> [(g b , g (g b))] -> [[b]]
sigmaBound s f =  [zipWith (+) (getMean s f) (getCov s f), getMean s  f,zipWith (-) (getMean s f) (getCov s f)]

sigmas :: Num b => (forall b . g b -> V1 b) -> [[(g b , g (g b))]] -> [[b]]
sigmas s = concat . map (sigmaBound s)

allSigmas :: Num b => [forall b . g b -> V1 b] -> [[(g b , g (g b))]] -> [[[b]]]
allSigmas s d =  map (flip sigmas $ d) s

unV1 (V1 x) = x
execute = do
   seed <- create
   let n = 500
   v <- replicateM n (normal 0 0.25 seed)
   let noisedData = testeSet v
       cleanData = testeSet (replicate n 0)
   let ny :: (V1 (V1 Double))
       ny = 1
       xo :: Model Double
       xo = 0.1 :|: 0.1 :|: 1 :|: 1e-3
       alpha = 1e-2
       iteration  r = full (predSin 1) nx projSin ny r
       backwardIteration = smoothing (predSin 1 ) nx
       forwards =   scanl (flip iteration)  (xo,nx) (fmap V1 $ testeSet v)
       backwards =  reverse $ scanl1 (flip backwardIteration)  $ reverse forwards
       project = fmap (unV1 . projSin . fst )
       getAmp = fmap ( unV1 . amp . fst )
       getAmpCov = fmap ( unV1 .  amp . unV1 . amp . snd)
       get x = fmap ( unV1 . x . fst )
   -- print projections
   plotLists [] [testeSet (replicate n 0) ,testeSet v,project backwards  ]
   -- print all parameters
   mapM_ (plotLists []) (allSigmas [amp ,accel,omega] [backwards])
   let ampCov = V.fromList noiseSignalBackward
       noiseSignalBackward = zipWith (-)  cleanData  (project backwards)
       noiseSignalForward = zipWith (-)  cleanData (project forwards)
   -- Print residuals
   plotLists [] [noiseSignalForward , noiseSignalBackward ,v]
   print $ mean ampCov
   print $ stdDev ampCov
   print $ mean $ V.fromList $  zipWith (-) (getAmp forwards) (getAmp backwards)
   print $ stdDev $ V.fromList $  zipWith (-) (getAmp forwards) (getAmp backwards)

amp :: Model a -> V1 a
amp ( x :|: _ :|: amp :|: _ ) = amp
accel ( x :|: _ :|: amp :|: accel ) = accel
omega ( x :|: omega  :|: amp :|: accel ) = omega
theta ( theta :|: _ :|: amp :|: accel ) = theta

matI :: Num a => Int -> [[a]]
matI n = [ [fromIntegral $ fromEnum $ i == j | i <- [1..n]] | j <- [1..n]]

diag = liftIdent
liftIdent ::(R f,Num a,Storable a,Vectorize f,Functor f,Applicative f)=> f a -> f ( f a)
liftIdent x = fmap (liftA2 (*) x) (fromLists $ matI n )
    where n = dimM x

testeSet v =  zipWith (+)  (map (function.fromIntegral) [1..length v]) v where
  function x =  (1+x/500) * sin(16 * x/500 * pi)

max_stars_repo_name: massudaw/mtk

max_issues_repo_name: massudaw/mtk

max_forks_repo_name: massudaw/mtk


text:
{-# LANGUAGE FlexibleInstances,
             DeriveDataTypeable #-}

module XAlg.Languages.Arithmetic.Number where

import XAlg.Foundation.Core

import Prelude hiding               (Float, Int)
import qualified  Data.Complex as C (Complex)
import Data.Complex          hiding (Complex)
import Data.Data

type Z = Integer
type R = Double
type C = C.Complex Double

deriving instance Ord (C.Complex Double)

data Number = Real R
            | Int Z
            | Complex C
            deriving (Eq, Ord, Data)

instance Show Number where
  show (Real r) = show r
  show (Int z) = show z
  show (Complex c) = show c

instance Z <: R where
  incl = fromInteger

instance Z <: C where
  incl = fromInteger

instance R <: C where
  incl = (0 :+)

compat :: (Number,Number) -> (Number,Number)
compat p@(Int n,Int m)          = p
compat   (Int n,r@(Real _))     = (Real (incl n),r)
compat   (Int n,c@(Complex _))  = (Complex (incl n),c)
compat   (Real n,Int m)         = (Real n, Real (incl m))
compat p@(Real _,Real _)        = p
compat   (Real n,c@(Complex _)) = (Complex (incl n),c)
compat p@(Complex _,Complex _)  = p
compat   (Complex n, Int m)     = (Complex n, Complex (incl m))
compat   (Complex n, Real m)    = (Complex n, Complex (incl m))

caseFun :: (Z -> Z) -> (R -> R) -> (C -> C) -> Number -> Number
caseFun zf rf cf (Real n)    = Real $ rf n
caseFun zf rf cf (Int n)     = Int $ zf n
caseFun zf rf cf (Complex n) = Complex $ cf n

caseFun2 :: (Z -> Z -> Z) -> (R -> R -> R) -> (C -> C -> C) -> Number -> Number -> Number
caseFun2 zf rf cf (Real n)    (Real m)    = Real $ rf n m
caseFun2 zf rf cf (Real n)    (Int m)     = Real $ rf n (incl m)
caseFun2 zf rf cf (Real n)    (Complex m) = Complex $ cf (incl n) m
caseFun2 zf rf cf (Int n)     (Real m)    = Real $ rf (incl n) m
caseFun2 zf rf cf (Int n)     (Int m)     = Int $ zf n m
caseFun2 zf rf cf (Int n)     (Complex m) = Complex $ cf (incl n) m
caseFun2 zf rf cf (Complex n) (Real m)    = Complex $ cf n (incl m)
caseFun2 zf rf cf (Complex n) (Int m)     = Complex $ cf n (incl m)
caseFun2 zf rf cf (Complex n) (Complex m) = Complex $ cf n m

instance Num Number where
  fromInteger = Int
  (+) = caseFun2 (+) (+) (+)
  (-) = caseFun2 (-) (-) (-)
  (*) = caseFun2 (*) (*) (*)
  abs = caseFun abs abs abs
  signum = caseFun signum signum signum

instance Fractional Number where
  (/) (Int x) (Int y) = Real $ fromInteger x / fromInteger y
  (/) x y = uncurry (/) . compat $ (x,y)
  recip (Int x) = Real . recip . fromInteger $ x
  recip x = caseFun undefined recip recip x
  fromRational = Real . fromRational

instance Floating Number where
  pi = Real pi
  exp (Int z) = Real $ exp $ fromInteger z
  exp (Real r) = Real $ exp r
  exp (Complex c) = Complex $ exp c
  log (Int z) = Real $ log $ fromInteger z
  log (Real r) = Real $ log r
  log (Complex c) = Complex $ log c
  sin (Int z) = Real $ sin $ fromInteger z
  sin (Real r) = Real $ sin r
  sin (Complex c) = Complex $ sin c
  cos (Int z) = Real $ cos $ fromInteger z
  cos (Real r) = Real $ cos r
  cos (Complex c) = Complex $ cos c
  asin (Int z) = Real $ asin $ fromInteger z
  asin (Real r) = Real $ asin r
  asin (Complex c) = Complex $ asin c
  acos (Int z) = Real $ acos $ fromInteger z
  acos (Real r) = Real $ acos r
  acos (Complex c) = Complex $ acos c
  atan (Int z) = Real $ atan $ fromInteger z
  atan (Real r) = Real $ atan r
  atan (Complex c) = Complex $ atan c
  sinh (Int z) = Real $ sinh $ fromInteger z
  sinh (Real r) = Real $ sinh r
  sinh (Complex c) = Complex $ sinh c
  cosh (Int z) = Real $ cosh $ fromInteger z
  cosh (Real r) = Real $ cosh r
  cosh (Complex c) = Complex $ cosh c
  asinh (Int z) = Real $ asinh $ fromInteger z
  asinh (Real r) = Real $ asinh r
  asinh (Complex c) = Complex $ asinh c
  acosh (Int z) = Real $ acosh $ fromInteger z
  acosh (Real r) = Real $ acosh r
  acosh (Complex c) = Complex $ acosh c
  atanh (Int z) = Real $ atanh $ fromInteger z
  atanh (Real r) = Real $ atanh r
  atanh (Complex c) = Complex $ atanh c

max_stars_repo_name: BlackBrane/xalg

max_issues_repo_name: BlackBrane/xalg

max_forks_repo_name: BlackBrane/xalg


text:
module MachineLearning.Optimization.MinibatchGradientDescentTest
(
  tests
)

where

import Test.Framework (testGroup)
import Test.Framework.Providers.HUnit
import Test.HUnit
import Test.HUnit.Approx
import Test.HUnit.Plus

import qualified Data.Vector.Storable as V
import qualified Numeric.LinearAlgebra as LA

import MachineLearning.Types (Vector)
import qualified MachineLearning as ML
import MachineLearning.Regularization (Regularization(..))
import MachineLearning.LeastSquaresModel (LeastSquaresModel(..))
import MachineLearning.Optimization.MinibatchGradientDescent

import MachineLearning.DataSets (dataset1)

(x, y) = ML.splitToXY dataset1

muSigma = ML.meanStddev x
xNorm = ML.featureNormalization muSigma x
x1 = ML.addBiasDimension xNorm
initialTheta :: Vector
initialTheta = LA.konst 0 (LA.cols x1)
lsExpectedTheta = LA.vector [325010.120,113889.649,5234.404]
eps = 1e-3

isInDescendingOrder :: [Double] -> Bool
isInDescendingOrder lst = and . snd . unzip $ scanl (\(prev, _) current -> (current, prev >= current)) (1/0, True) lst

testMinibatchGradientDescent model expectedTheta = do
  let (theta, optPath) = minibatchGradientDescent 0 16 0.01 model eps 5000 RegNone x1 y initialTheta
      js = V.toList $ (LA.toColumns optPath) !! 1
  assertVector "theta" 0.01 expectedTheta theta
  assertBool "non-increasing errors" $ isInDescendingOrder js

tests = [testGroup "minibatchGradientDescent" [
            testCase "leastSquares" $ testMinibatchGradientDescent LeastSquares lsExpectedTheta
            ]
        ]


text:
{-|
Module: MachineLearning.NeuralNetwork.LogisticLoss
Description: Logistic Loss.
Copyright: (c) Alexander Ignatyev, 2017
License: BSD-3
Stability: experimental
Portability: POSIX

Logistic Loss.
-}

module MachineLearning.NeuralNetwork.LogisticLoss
(
  scores
  , gradient
  , loss
)

where

import qualified Numeric.LinearAlgebra as LA
import MachineLearning.Types (R, Matrix)
import qualified MachineLearning.LogisticModel as LM

scores :: Matrix -> Matrix
scores = LM.sigmoid

gradient :: Matrix -> Matrix -> Matrix
gradient scores y = scores - y

-- Logistic Loss function
loss :: Matrix -> Matrix -> R
loss x y = (LA.sumElements $ (-y) * log(tau + x) - (1-y) * log ((1+tau)-x))/m
  where tau = 1e-7
        m = fromIntegral $ LA.rows x