text:
{-# language FlexibleContexts, GeneralizedNewtypeDeriving, DeriveFunctor #-}
{-# language OverloadedStrings #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Control.Iterative
-- Copyright   :  (c) Marco Zocca 2017-2018
-- License     :  GPL-style (see the file LICENSE)
--
-- Maintainer  :  zocca marco gmail
-- Stability   :  experimental
-- Portability :  portable
--
-- Combinators and helper functions for iterative algorithms, with support for monitoring and exceptions.
--
-----------------------------------------------------------------------------
module Control.Iterative where

import Control.Applicative

import Control.Monad (when)
-- import Control.Monad.Trans.Reader
import Control.Monad.Reader
import Control.Monad.State.Strict
-- import Control.Monad.Trans.Writer.CPS
import Control.Monad.Trans.Class (MonadTrans(..), lift)
import qualified Control.Monad.Trans.State.Strict as MTS -- (runStateT)
import Control.Monad.Catch
import qualified Control.Monad.Log as L (MonadLog(..), WithSeverity(..), Severity(..), renderWithSeverity, LoggingT, runLoggingT, Handler, logMessage, logError, logDebug, logInfo, logNotice)

-- import Data.Bool (bool)
import Data.Char (toUpper)

import Data.Typeable
import qualified Control.Exception as E (Exception, Handler)

import Data.Foldable (foldrM)

-- import Control.Exception.Common
-- import Numeric.LinearAlgebra.Class
-- import Numeric.Eps

data ConvergenceStatus a = BufferNotReady
                         | Converging
                         | Converged a
                         | Diverging a a
                         | NotConverged
                           deriving (Eq, Show, Typeable)
instance (Typeable a, Show a) => Exception (ConvergenceStatus a)

-- FIXME I don't like this; logging shouldn't be in IO, `iterationView` is too general, etc.
data IterationConfig a b =
  IterConf { numIterationsMax :: Int     -- ^ Max.# of iterations
           , printDebugInfo :: Bool      -- ^ Print iteration info to stdout
           , iterationView :: a -> b     -- ^ Project state to a type `b`
           , printDebugIO :: b -> IO ()} -- ^ print function for type `b`
instance Show (IterationConfig a b) where
  show (IterConf n qd _ _) = unwords ["Max. # of iterations:",show n,", print debug information:", show qd]

-- | Iterative algorithms need configuration and logging; here we use a transformer stack or ReaderT on top of LoggingT (from `logging-effect`).
newtype App c msg m a =
  App {
    unApp :: ReaderT c (L.LoggingT msg m) a
      } deriving (Functor, Applicative, Alternative, Monad)

-- instance MonadTrans (App c msg) where

liftApp :: (c -> L.LoggingT msg m a) -> App c msg m a
liftApp = App . ReaderT

testApp :: (Monad m, MonadState b (t m), MonadTrans t, L.MonadLog msg (t m)) =>
     (b -> m b) -> (t1 -> b -> msg) -> t1 -> t m b
testApp fm logf c = do
  x <- get
  y <- lift $ fm x
  put y
  L.logMessage (logf c y)
  return y

-- | Configure and log a computation
runApp :: L.Handler m msg -> c -> App c msg m a -> m a
runApp logHandler config m = L.runLoggingT (runReaderT (unApp m) config) logHandler

runLogDebug :: L.MonadLog (L.WithSeverity msg) m => c -> App c msg m a -> m a
runLogDebug = runApp L.logDebug

-- -- | runApp, logging by printing to terminal
-- runDebug :: Show a => c -> App c a IO b -> IO b
-- runDebug = runApp print

-- | App1: LoggingT / ReaderT

-- newtype App1 c msg m a = App1 { unApp1 :: L.LoggingT msg (ReaderT c m) a } deriving (Functor, Applicative, Alternative, Monad)

-- liftApp1 logf f = App1 $ do
--   c <- ask
--   y <- lift . lift $ f c
--   L.logMessage (logf y)
--   return y

-- | App2: ReaderT / LoggingT / StateT

newtype App2 c msg s m a = App2 {
  unApp2 :: ReaderT c (L.LoggingT msg (StateT s m)) a
  } deriving (Functor, Applicative, Alternative, Monad)

runApp2 lh r x0 m = runStateT (L.runLoggingT (runReaderT (unApp2 m) r) lh) x0

-- app2 :: Monad m =>
--         (c -> a -> log)  -- ^ Builds a log entry from the configuration and the current state
--      -> (r -> t) -- ^ Builds
--      -> (r -> c)
--      -> (a -> t -> m a)
--      -> App2 r log a m a
app2 logf cconf lconf fm = App2 $ do
  x <- get
  cc <- asks cconf
  lc <- asks lconf
  y <- lift . lift . lift $ fm x cc
  L.logMessage (logf lc y)
  put y
  return y

-- logging test

-- -- logApp :: L.MonadLog (L.WithSeverity String) m => m ()
logApp = do
  L.logError "moo"
  L.logInfo "info"

asdf = L.runLoggingT logApp (putStrLn . withSeverity id)

-- -- renderWithSeverity k (L.WithSeverity u a)

bracketsUpp :: Show a => a -> String
bracketsUpp p = unwords ["[", map toUpper (show p), "]"]

withSeverity :: (t -> String) -> L.WithSeverity t -> String
withSeverity k (L.WithSeverity u a ) = unwords [bracketsUpp u, k a]

-- -- >>> renderWithSeverity id (WithSeverity Informational "Flux capacitor is functional")
-- -- [Informational] Flux capacitor is functional
-- renderWithSeverity
--   :: (a -> PP.Doc) -> (WithSeverity a -> PP.Doc)
-- renderWithSeverity k (WithSeverity u a) =
--   PP.brackets (PP.pretty u) PP.<+> PP.align (k a)

-- * Control primitives for bounded iteration with convergence check

-- -- | transform state until a condition is met
modifyUntil :: MonadState s m => (s -> Bool) -> (s -> s) -> m s
modifyUntil q f = modifyUntilM q (pure . f)

modifyUntilM :: MonadState s m => (s -> Bool) -> (s -> m s) -> m s
modifyUntilM q f = do
  x <- get
  y <- f x
  put y
  if q y then return y
         else modifyUntilM q f

modifyUntilM_ :: MonadState s m => (s -> Bool) -> (s -> m s) -> m s
modifyUntilM_ q f = do
  x <- get
  y <- f x
  if q y
    then pure y
    else do
        put y
        modifyUntilM_ q f

-- | modifyUntil with optional iteration logging to stdout
modifyUntil' :: L.MonadLog String m =>
   IterationConfig a b -> (a -> Bool) -> (a -> a) -> a -> m a
modifyUntil' config q f x0 = modifyUntilM' config q (pure . f) x0

modifyUntilM' :: L.MonadLog String m =>
   IterationConfig a b -> (a -> Bool) -> (a -> m a) -> a -> m a
modifyUntilM' config q f x0 = MTS.execStateT (go 0) x0 where
  pf = iterationView config
  go i = do
   x <- get
   y <- lift $ f x
   when (printDebugInfo config) $ do
     L.logMessage $ unwords ["Iteration", show i, "\n"]
   -- when (printDebugInfo config) $ liftIO $ do
   --   putStrLn $ unwords ["Iteration", show i, "\n"]
   --   printDebugIO config (pf y)
   put y
   if q y
     then return y
     else go (i + 1)

-- | `untilConvergedG0` is a special case of `untilConvergedG` that assesses convergence based on the L2 distance to a known solution `xKnown`
untilConvergedG0 :: (Normed v, MonadThrow m, L.MonadLog String m, Typeable (Magnitude v), Typeable s, Show s) =>
     String
     -> IterationConfig s v
     -> v                    -- ^ Known value
     -> (s -> s)
     -> s
     -> m s
untilConvergedG0 fname config xKnown f x0 =
  modifyInspectGuarded fname config norm2Diff nearZero qdiverg qfin f x0
   where
    qfin s = nearZero $ norm2 (xKnown ^-^ s)

-- | This function makes some default choices on the `modifyInspectGuarded` machinery: convergence is assessed using the squared L2 distance between consecutive states, and divergence is detected when this function is increasing between pairs of measurements.
untilConvergedG :: (Normed v, MonadThrow m, L.MonadLog String m, Typeable (Magnitude v), Typeable s, Show s) =>
        String
     -> IterationConfig s v
     -> (v -> Bool)
     -> (s -> s)
     -> s
     -> m s
untilConvergedG fname config =
  modifyInspectGuarded fname config norm2Diff nearZero qdiverg

-- | ", monadic version
untilConvergedGM ::
  (Normed v, MonadThrow m, L.MonadLog String m, Typeable (Magnitude v), Typeable s, Show s) =>
     String
     -> IterationConfig s v
     -> (v -> Bool)
     -> (s -> m s)
     -> s
     -> m s
untilConvergedGM fname config =
  modifyInspectGuardedM fname config norm2Diff nearZero qdiverg

-- | `modifyInspectGuarded` is a high-order abstraction of a numerical iterative process. It accumulates a rolling window of 3 states and compares a summary `q` of the latest 2 with that of the previous two in order to assess divergence (e.g. if `q latest2 > q prev2` then the function throws an exception and terminates). The process ends by either hitting an iteration budget or by relative convergence, whichever happens first. After the iterations stop, the function then assesses the final state with a predicate `qfinal` (e.g. for comparing the final state with a known one; if this is not available, the user can just supply `const True`)
modifyInspectGuarded ::
  (MonadThrow m, L.MonadLog String m, Typeable s, Typeable a, Show s, Show a) =>
        String              -- ^ Calling function name
     -> IterationConfig s v -- ^ Configuration
     -> ([v] -> a)          -- ^ State summary array projection
     -> (a -> Bool)         -- ^ Convergence criterion
     -> (a -> a -> Bool)    -- ^ Divergence criterion
     -> (v -> Bool)         -- ^ Final state acceptance criterion
     -> (s -> s)            -- ^ State evolution
     -> s                   -- ^ Initial state
     -> m s                 -- ^ Final state
modifyInspectGuarded fname config sf qc qd qfin f x0 =
  modifyInspectGuardedM fname config sf qc qd qfin (pure . f) x0

-- | ", monadic version
modifyInspectGuardedM ::
  (MonadThrow m, L.MonadLog String m, Typeable s, Show s, Typeable a, Show a) =>
     String
     -> IterationConfig s v
     -> ([v] -> a)
     -> (a -> Bool)
     -> (a -> a -> Bool)
     -> (v -> Bool)
     -> (s -> m s)
     -> s
     -> m s
modifyInspectGuardedM fname config sf qconverg qdiverg qfinal f x0
  | nitermax > 0 = MTS.execStateT (go 0 []) x0
  | otherwise = throwM (NonNegError fname nitermax)
  where
    lwindow = 3
    nitermax = numIterationsMax config
    pf = iterationView config
    checkConvergStatus y i ll
      | length ll < lwindow = BufferNotReady
      | qdiverg qi qt && not (qconverg qi) = Diverging qi qt
      | qconverg qi || qfinal (pf y) = Converged qi
      | i == nitermax - 1 = NotConverged
      | otherwise = Converging
      where llf = pf <$> ll
            qi = sf $ init llf         -- summary of latest 2 states
            qt = sf $ tail llf         -- "       "  previous 2 states
    go i ll = do
      x <- MTS.get
      y <- lift $ f x
      when (printDebugInfo config) $ do
        L.logMessage $ unwords ["Iteration", show i]
      -- when (printDebugInfo config) $ liftIO $ do
      --   putStrLn $ unwords ["Iteration", show i]
      --   printDebugIO config (pf y)
      case checkConvergStatus y i ll of
        BufferNotReady -> do
          MTS.put y
          let ll' = y : ll    -- cons current state to buffer
          go (i + 1) ll'
        Converged qi -> MTS.put y
        Diverging qi qt -> do
          MTS.put y
          throwM (DivergingE fname i qi qt)
        Converging -> do
          MTS.put y
          let ll' = init (y : ll) -- rolling state window
          go (i + 1) ll'
        NotConverged -> do
          MTS.put y
          throwM (NotConvergedE fname nitermax y)

-- | Some useful combinators

-- | Apply a function over a range of integer indices, zip the result with it and filter out the almost-zero entries
onRangeSparse :: Epsilon b => (Int -> b) -> [Int] -> [(Int, b)]
onRangeSparse f ixs = foldr ins [] ixs where
  ins x xr | isNz (f x) = (x, f x) : xr
           | otherwise = xr

-- | ", monadic version
onRangeSparseM :: (Epsilon b, Foldable t, Monad m) =>
     (a -> m b) -> t a -> m [(a, b)]
onRangeSparseM f ixs = unfoldZipM mf f ixs where
  mf x = isNz <$> f x

unfoldZipM0 :: (Foldable t, Monad m) =>
     (a -> Bool) -> (a -> b) -> t a -> m [(a, b)]
unfoldZipM0 q f = unfoldZipM (pure . q) (pure . f)

unfoldZipM :: (Foldable t, Monad m) =>
     (a -> m Bool) -> (a -> m b) -> t a -> m [(a, b)]
unfoldZipM q f ixs = foldrM insf [] ixs where
  insf x xr = do
    qx <- q x
    if qx
    then do
      y <- f x
      pure $ (x, y) : xr
    else pure xr

-- | A combinator I don't know how to call
combx :: Functor f => (a -> b) -> (t -> f a) -> t -> f b
combx g f x = g <$> f x

-- | Helpers

-- | Relative residual
relRes :: (Normed t, LinearVectorSpace t) =>
     MatrixType t -> t -> t -> Magnitude t
relRes aa b x = n / d where
  n = norm2 $ (aa #> x) ^-^ b
  d = norm2 b

-- meanl :: (Foldable t, Fractional a) => t a -> a
-- meanl xx = 1/fromIntegral (length xx) * sum xx

-- norm2l :: (Foldable t, Functor t, Floating a) => t a -> a
-- norm2l xx = sqrt $ sum (fmap (**2) xx)

-- | Squared difference of a 2-element list.
-- | NB: unsafe !
diffSqL :: Floating a => [a] -> a
diffSqL xx = (x1 - x2)**2 where [x1, x2] = [head xx, xx!!1]

-- | Relative tolerance :
-- relTol a b := ||a - b|| / (1 + min (||norm2 a||, ||norm2 b||))
relTol :: Normed v => v -> v -> Magnitude v
relTol a b = norm2 (a ^-^ b) / m where
  m = 1 + min (norm2 a) (norm2 b)

qdiverg :: Ord a => a -> a -> Bool
qdiverg = (>)

norm2Diff [s1, s0] = norm2 (s1 ^-^ s0)
norm2Diff _ = 1/0

-- test data

data S = S {unS1 :: Double, unS2 :: String} deriving (Eq, Show)
liftS1 f (S x i) = S (f x) i
s0 = S 1 "blah"
ic1 = IterConf 2 True unS1 print

-- playground

-- instance MonadThrow m => MonadThrow (WriterT w m) where
--   throwM = lift . throwM

-- -- | iter0 also accepts a configuration, e.g. for optional printing of debug info
-- -- iter0 :: MonadIO m =>
-- --      Int -> (s -> m s) -> (s -> String) -> IterationConfig s -> s -> m s
-- iter0 nmax f sf config x0 = flip runReaderT config $ MTS.execStateT (go (0 :: Int)) x0
--  where
--   go i = do
--     x <- get
--     c <- lift $ asks printDebugInfo  -- neat
--     y <- lift . lift $ f x           -- not neat
--     when c $ liftIO $ putStrLn $ sf y
--     put y
--     unless (i >= nmax) (go $ i + 1)

-- -- | iter1 prints output at every iteration until the loop terminates OR is interrupted by an exception, whichever happens first
-- -- iter1 :: (MonadThrow m, MonadIO m, Typeable t, Show t) =>
-- --      (t -> m t) -> (t -> String) -> (t -> Bool) -> (t -> Bool) -> t -> m t
-- iter1 f wf qe qx x0 = execStateT (go 0) x0 where
--  go i = do
--    x <- get
--    y <- lift $ f x
--    _ <- liftIO $ wf y
--    when (qx y) $ throwM (NotConvergedE "bla" (i+1)  y)
--    put y
--    unless (qe y) $ go (i + 1)

-- -- | iter2 concatenates output with WriterT but does NOT `tell` any output if an exception is raised before the end of the loop
-- iter2 :: (MonadThrow m, Monoid w, Typeable t, Show t) => (t -> m t)
--      -> (t -> w) -> (t -> Bool) -> (t -> Bool) -> t -> m (t, w)
-- iter2 f wf qe qx x0 = runWriterT $ execStateT (go 0) x0 where
--  go i = do
--    x <- get
--    y <- lift . lift $ f x
--    lift $ tell $ wf y
--    when (qx y) $ throwM (NotConvergedE "bla" (i+1)  y)
--    put y
--    unless (qe y) $ go (i + 1)

-- -- test :: IO (Int, [String])
-- test :: IO Int
-- -- test :: IO ()
-- test = do
--   (yt, w ) <- iter2 f wf qe qexc x0
--   putStrLn w
--   return yt
--   -- iter1 f wf qe qexc x0
--   where
--     f = pure . (+ 1)
--     wf v = unwords ["state =", show v]
--     qe = (== 5)
--     qexc = (== 3)
--     x0 = 0 :: Int

max_stars_repo_name: ocramz/iterative-algorithms

max_issues_repo_name: ocramz/iterative-algorithms

max_forks_repo_name: ocramz/iterative-algorithms


text:
{-# LANGUAGE TypeFamilies #-}
module Data.ExponentialFamily.Distribution.Binomial where

import Data.ExponentialFamily.Density
import Data.ExponentialFamily.ThetaEta

import Numeric.SpecFunctions (choose)

-- | Binomial distribution, where you conduct n trials with probability p.
data Binomial = Binomial { _n :: Int
                         , _p :: Double
                         } deriving (Show, Eq)

instance ProbDensity Binomial where
    -- | The probability of succeeding @success@ ∈ ℕ
    type Point Binomial = Int
    f (Binomial n p) success =
        (n `choose` success) * p^success * (1-p)^(n-success)

-- | (η,N) and (θ,N) for 'known' number of trials : N.
instance ThetaEta Binomial where
    type TEParam Binomial = (Double, Int)
    toθ (Binomial n p)  = Theta (log (p / (1-p)), n)
    fromθ (Theta (logit, n))
        = Binomial n (exp logit / (1+exp logit))
    toη (Binomial n p)  = Eta (fromIntegral n*p, n)
    fromη (Eta (np, n)) = Binomial n (np/fromIntegral n)

ext: hs

max_stars_repo_name: lesguillemets/HExpFams

max_issues_repo_name: lesguillemets/HExpFams

max_forks_repo_name: lesguillemets/HExpFams


text:
{-# LANGUAGE RecordWildCards #-}
import Data.Complex
import Control.Monad.Trans.Except
import Foreign.C.Types
import Data.Word
import System.IO
import Data.Monoid

import Control.Error.Util
import Data.Vector.Storable as VS hiding ((++))
import Data.Vector.Generic as VG hiding ((++))
import Pipes
import qualified Pipes.Prelude as P
import Foreign.Storable.Complex

import Options.Applicative

import SDR.Filter
import SDR.RTLSDRStream
import SDR.Util
import SDR.Demod
import SDR.Pulse
import SDR.Serialize
import SDR.ArgUtils
import SDR.CPUID

--The filter coefficients are stored in another module
import Coeffs

data Options = Options {
    frequency  :: Word32,
    input      :: Maybe FilePath,
    output     :: Maybe FilePath
}

optParser :: Parser Options
optParser = Options
          <$> option (fmap fromIntegral parseSize) (
                 long "frequency"
              <> short 'f'
              <> metavar "FREQUENCY"
              <> help "Frequency to tune to"
              )
          <*> optional (strOption (
                 long "input"
              <> short 'i'
              <> metavar "FILENAME"
              <> help "Input filename"
              ))
          <*> optional (strOption (
                 long "output"
              <> short 'o'
              <> metavar "FILENAME"
              <> help "Output filename"
              ))

opt :: ParserInfo Options
opt = info (helper <*> optParser) (fullDesc <> progDesc "Receive and demodulate broadcast FM radio" <> header "RTLSDR FM")

bufNum     = 1
bufLen     = 16384
samples    = fromIntegral (bufNum * bufLen) `quot` 2
decimation = 8
sqd        = samples `quot` decimation

{-
    Sampling frequency of the input is 1280 khz
    This is decimated by a factor of 8 and then demodulated
    Sampling frequency of demodulated signal is 160 khz
    Need audio output at 48 khz
    Resampling factor is 48/160 == 3/10
    FM pilot tone at 19khz (0.3958 * 48)
    Start audio filter cutoff at 15khz (0.3125 * 48)
-}

doIt Options{..} = do

    info <- lift getCPUInfo

    let rtlstream = do
            str <- sdrStream (defaultRTLSDRParams frequency 1280000) bufNum bufLen
            return $ str >-> P.map (interleavedIQUnsignedByteToFloatFast info)

    let fileStream fname = lift $ do
            h <- openFile fname ReadMode
            return $ fromHandle samples h
            --TODO: how do I ensure these handles get closed?

    inputSpectrum <- maybe rtlstream fileStream input

    let fileSink fname = do
            h <- openFile fname ReadMode
            return $ toHandle h

    sink <- lift $ maybe pulseAudioSink fileSink output

    deci <- lift $ fastDecimatorC info decimation coeffsRFDecim
    resp <- lift $ fastResamplerR info 3 10 coeffsAudioResampler
    filt <- lift $ fastFilterSymR info coeffsAudioFilter

    --Build the pipeline
    let pipeline :: Effect IO ()
        pipeline =   inputSpectrum
                 >-> firDecimator deci sqd
                 >-> fmDemod
                 >-> firResampler resp sqd
                 >-> firFilter filt sqd
                 >-> P.map (VG.map (* 0.2))
                 >-> sink

    --Run the pipeline
    lift $ runEffect pipeline

main = execParser opt >>= exceptT putStrLn return . doIt

max_stars_repo_name: adamwalker/sdr-apps

max_issues_repo_name: adamwalker/sdr-apps

max_forks_repo_name: adamwalker/sdr-apps


text:
module Pineappl.Quantum where

import           Data.Bifunctor
import           Data.Complex
import qualified Data.Map                      as M
import           Pineappl                       ( Prob(P)
                                                , LogProb(LP)
                                                )
import           Pineappl.Deterministic         ( runBDDist
                                                , bddist
                                                , hist
                                                , WrappedBDDist
                                                )

fromAmplitude
  :: (RealFloat a, Ord b)
  => WrappedBDDist (LogProb (Complex a)) b
  -> WrappedBDDist (Prob a) b
fromAmplitude =
  bddist
    . fmap (second (\(LP psi) -> P (magnitude psi ^ 2)))
    . M.toList
    . hist
    . runBDDist

max_stars_repo_name: mcwitt/pineappl

max_issues_repo_name: mcwitt/pineappl

max_forks_repo_name: mcwitt/pineappl


text:
import Criterion.Main
import qualified Data.Vector.Unboxed as U
import Numeric.SpecFunctions
import Numeric.Polynomial
import Text.Printf

-- Uniformly sample logGamma performance between 10^-6 to 10^6
benchmarkLogGamma logG =
  [ bench (printf "%.3g" x) $ nf logG x
  | x <- [ m * 10**n | n <- [ -8 .. 8 ]
                     , m <- [ 10**(i / tics) | i <- [0 .. tics-1] ]
         ]
  ]
  where tics = 3
{-# INLINE benchmarkLogGamma #-}

-- Power of polynomial to be evaluated (In other words length of coefficients vector)
coef_size :: [Int]
coef_size = [ 1,2,3,4,5,6,7,8,9
            , 10,    30
            , 100,   300
            , 1000,  3000
            , 10000, 30000
            ]
{-# INLINE coef_size #-}

-- Precalculated coefficients
coef_list :: [U.Vector Double]
coef_list = [ U.replicate n 1.2 | n <- coef_size]
{-# NOINLINE coef_list #-}

main :: IO ()
main = defaultMain
  [ bgroup "logGamma" $
    benchmarkLogGamma logGamma
  , bgroup "logGammaL" $
    benchmarkLogGamma logGammaL
  , bgroup "incompleteGamma" $
      [ bench (show p) $ nf (incompleteGamma p) p
      | p <- [ 0.1
             , 1,   3
             , 10,  30
             , 100, 300
             , 999, 1000
             ]
      ]
  , bgroup "factorial"
    [ bench (show n) $ nf factorial n
    | n <- [ 0, 1, 3, 6, 9, 11, 15
           , 20, 30, 40, 50, 60, 70, 80, 90, 100
           ]
    ]
  , bgroup "incompleteBeta"
    [ bench (show (p,q,x)) $ nf (incompleteBeta p q) x
    | (p,q,x) <- [ (10,      10,      0.5)
                 , (101,     101,     0.5)
                 , (1010,    1010,    0.5)
                 , (10100,   10100,   0.5)
                 , (100100,  100100,  0.5)
                 , (1001000, 1001000, 0.5)
                 , (10010000,10010000,0.5)
                 ]
    ]
  , bgroup "log1p"
      [ bench (show x) $ nf log1p x
      | x <- [ -0.9
             , -0.5
             , -0.1
             ,  0.1
             ,  0.5
             ,  1
             ,  10
             ,  100
             ]
      ]
  , bgroup "poly"
      $  [ bench ("vector_"++show (U.length coefs)) $ nf (\x -> evaluatePolynomial x coefs) (1 :: Double)
         | coefs <- coef_list ]
      ++ [ bench ("unpacked_"++show n) $ nf (\x -> evaluatePolynomialL x (map fromIntegral [1..n])) (1 :: Double)
         | n <- coef_size ]
  ]

max_stars_repo_name: chris-taylor/math-functions

max_issues_repo_name: chris-taylor/math-functions

max_forks_repo_name: chris-taylor/math-functions


text:
module Measurable.Measures where

import Measurable.Core
import Numeric.SpecFunctions (choose)
import Statistics.Distribution
import qualified Statistics.Distribution.Beta as Statistics
import qualified Statistics.Distribution.Binomial as Statistics
import qualified Statistics.Distribution.ChiSquared as Statistics
import qualified Statistics.Distribution.Gamma as Statistics
import qualified Statistics.Distribution.Exponential as Statistics
import qualified Statistics.Distribution.Normal as Statistics

standard :: Measure Double
standard = fromDensityFunction pdf where
  pdf = density Statistics.standard

normal :: Double -> Double -> Measure Double
normal m s = fromDensityFunction pdf where
  pdf = density $ Statistics.normalDistr m s

logNormal :: Double -> Double -> Measure Double
logNormal m s = fmap exp (normal m s)

exponential :: Double -> Measure Double
exponential r = fromDensityFunction pdf where
  pdf = density $ Statistics.exponential r

gamma :: Double -> Double -> Measure Double
gamma a b = fromDensityFunction pdf where
  pdf = density $ Statistics.gammaDistr a b

inverseGamma :: Double -> Double -> Measure Double
inverseGamma a b = fmap recip (gamma a b)

chiSquare :: Int -> Measure Double
chiSquare k = fromDensityFunction pdf where
  pdf = density $ Statistics.chiSquared k

beta :: Double -> Double -> Measure Double
beta a b = fromDensityFunction pdf where
  pdf = density $ Statistics.betaDistr a b

binomial :: Int -> Double -> Measure Int
binomial n p = fromMassFunction pmf [0..n] where
  pmf = binomialMass n p

bernoulli :: Double -> Measure Int
bernoulli = binomial 1

binomialMass :: Int -> Double -> Int -> Double
binomialMass n p k = bc * p ^ k * (1 - p) ^ (n - k) where
  bc = n `choose` k

max_stars_repo_name: jtobin/measurable

max_issues_repo_name: jtobin/measurable

max_forks_repo_name: jtobin/measurable


text:
{-# LANGUAGE TupleSections #-}
{-# LANGUAGE TemplateHaskell #-}
{-

Space in the middle: 9
Off-center: 7

-}

import Codec.Picture
import Codec.Picture.Gif
import Data.List.Split
import Data.Complex
import System.Environment
import System.FilePath
import qualified Data.ByteString as BS
import qualified Data.Bits.Bitwise as BA
import Options.Applicative
import Control.Monad
import Data.Monoid
import System.Exit
import Data.Char
import Data.FileEmbed
import Language.Haskell.TH
import Language.Haskell.TH.Syntax
import Text.Regex.Posix
import Numeric (readHex)

cFRAMES = 256          :: Integer
cCENTER = 9            :: Double
cLEN    = 33+cCENTER/2 :: Double  -- Lenght of one arm
cT0RATE = 1000000      :: Integer -- TODO: Need to calibrate this

data Color = R|G|B deriving (Show, Eq)

type Offset = Double
type Shift = Double
type Rotation = Double
type Speed = Double

type BitMaker = (Color -> Bool -> Complex Double -> Bool) -> [Bool]

-- Frame (0..1) to position (-cLen..cLen) to coordinate (-1..1 × -1..1)
type Transformation = Double -> Double -> Complex Double

projectCircular :: Offset -> Shift -> Rotation -> Transformation
projectCircular off shift rot =
    \ frame pos ->
        let rho = frame*2*pi
            cpos = pos :+ 0
        in ((cpos + displacement) * cis rho * cis rot_radians) / (radius :+ 0)
  where
    displacement = (shift :+ off)
    rot_radians = (-rot-3)/12*2*pi
    radius = sqrt ((cLEN + abs shift)^2 + off^2)

projectLinear :: Transformation
projectLinear frame pos = (2*frame -1) :+ pos/cLEN

data TransSpec
    = Circular Offset Shift Rotation
    | Linear
    deriving Show

project :: TransSpec -> Transformation
project (Circular off shift rot) = projectCircular off shift rot
project Linear                   = projectLinear

bitBuilder :: TransSpec -> BitMaker
bitBuilder spec f =
    [ f c a cpos
    | frame <- [0,1..cFRAMES-1]
    , c <- [R,G,B]
    , (a, pos) <- ((False,) <$> arm) ++ (((True,) . negate) <$> arm)
    , let cpos = trans (fromIntegral frame / fromIntegral cFRAMES) pos
    ]
  where
    trans = project spec
    arm = concatMap reverse $ chunksOf 8 $ map (cCENTER/2 +) [0..31]

getColor :: Color -> PixelRGB8 -> Bool
getColor R (PixelRGB8 r _ _ ) = r > 130
getColor G (PixelRGB8 _ g _ ) = g > 130
getColor B (PixelRGB8 _ _ b ) = b > 130

calibrationBits :: [Bool]
calibrationBits = map not $ foldr1 (zipWith (||)) $
    [ bitBuilder (Circular off 0 0) $ \c a z ->
        abs (phase z - off*pi/10) < pi/40 &&
       (c == oc) &&
       (if a then magnitude z > 0.75
             else (magnitude z < 0.75 && (magnitude z > 0.5 || off == 0)))
    | (oc, off) <- zip (cycle [R,G,B]) [-8..8] ]

dynImageToPackagedData builder dynImage =
    bitsToData $
    builder $ \col _ (x :+ y) ->
           let pix = pixelAt pixels
                  (round ((( x+1)/2) * w))
                  (round (((-y+1)/2) * h))
           in not $ getColor col pix
  where
    pixels = convertRGB8 dynImage

    w = fromIntegral (imageWidth pixels)
    h = fromIntegral (imageHeight pixels)

bitsToData = BS.pack . map BA.fromListLE . chunksOf 8

getImage :: BitMaker -> (FilePath, Double) -> IO [(BS.ByteString, Double)]
getImage builder source = case source of
        ("CALIBRATION",n) -> do
            putStrLn $ "Adding offset calibration image"
            pure [(bitsToData calibrationBits, n)]
        (filename, _) | takeExtension filename == ".gif" -> do
            gifData <- BS.readFile filename
            let gifImages = either error id $ decodeGifImages gifData
            let gifDelays = either error id $ getDelaysGifImages gifData
            putStrLn $ "Adding " ++ filename ++ " (" ++ show (length gifImages) ++ " frames)"

            pure $ zip (map (dynImageToPackagedData builder) gifImages)
                       (map ((/100) . fromIntegral) gifDelays)
        (filename, n) -> do
            putStrLn $ "Adding " ++ filename
            dynImage <- either error id <$> readImage filename
            let bits = dynImageToPackagedData builder dynImage
            pure [(bits, n)]

replace :: Int -> BS.ByteString -> BS.ByteString -> BS.ByteString
replace at replacement source
    | BS.length replacement > BS.length after
    = error "replaceIn: Replacement larger than remaining string"
    | otherwise = before <> replacement <> BS.drop (BS.length replacement) after
    where
        (before, after) = BS.splitAt at source

template :: BS.ByteString
template = $(embedFile "firmware/firmware.bin")

offsetInitialStep, offsetTiming, offsetImages :: Int
(offsetInitialStep, offsetTiming, offsetImages) = $(do
    qAddDependentFile "firmware/firmware.map"
    mapFile <- runIO $ readFile "firmware/firmware.map"
    let find name
            | [_,s] <- getAllTextSubmatches (mapFile =~ ("^C: +([0-9ABCDEF]{8}) +_" ++ name)) :: [String]
            , [(i,"")] <- readHex s
            = return i
            | otherwise
            = do reportError $ "Cound not find location of " ++ show name ++ " in firmware.map"
                 fail ""
    ois <- find "initial_step"
    ot  <- find "timing"
    oi  <- find "images"
    return $ TupE $ map LitE [integerL ois, integerL ot, integerL oi]
    )

work :: FilePath -> TransSpec -> Speed -> [(FilePath, Double)] -> IO ()
work output spec speed timed_sources = do
    let builder = bitBuilder spec
    timed_data <- concat <$> mapM (getImage builder) timed_sources
    let (imagesData, timings) = unzip timed_data
    let imageData = BS.concat imagesData

    when (length imagesData > 8) $ do
        putStrLn $ "Too many images"
        exitFailure

    let timingData    = BS.pack $ map fromIntegral $
            concatMap (\x -> [x`mod`256, x`div`256]) $
            map round $
            map (*256) $
            take 8 $ timings ++ repeat 0

    let initial_step_data = BS.pack $ map fromIntegral $
            concatMap (\x -> [x`mod`256, x`div`256]) $
            map round $
            [ speed / fromIntegral cFRAMES * fromIntegral cT0RATE ]

    putStrLn $ "Writing " ++ output
    BS.writeFile output $
        replace offsetImages imageData $
        replace offsetTiming timingData $
        replace offsetInitialStep initial_step_data $
        template

-- Argument handling

main :: IO ()
main = join . customExecParser (prefs showHelpOnError) $
  info (helper <*> parser)
  (  fullDesc
  <> header "YQ3008 image generator"
  -- <> progDesc "TODO"
  )
  where
    linearSpecParser :: Parser TransSpec
    linearSpecParser =
        flag' Linear
            (  long "linear"
            <> help "Scan the input image linearly, not circularly"
            )

    circularSpecParser :: Parser TransSpec
    circularSpecParser = Circular
        <$> option auto
            (  long "offset"
            <> metavar "OFFSET"
            <> help "vertical offset of the bar from the hub (in leds)"
            <> value 0
            )
        <*> option auto
            (  long "shift"
            <> metavar "SHIFT"
            <> help "horizontal shift of the bar from the hub (in leds)"
            <> value 0
            )
        <*> option auto
            (  long "rotation"
            <> metavar "ROTATION"
            <> help "position of the magnet [0..12]"
            <> value 12
            <> showDefault
            )

    parser :: Parser (IO ())
    parser = work
        <$> strOption
            (  long "output"
            <> short 'o'
            <> metavar "FILE"
            <> help "output file"
            )
        <*> (linearSpecParser <|> circularSpecParser)
        <*> option auto
            (  long "speed"
            <> metavar "SECONDS"
            <> help "how long one image scan is initially"
            <> value 0.5
            <> showDefault
            )
        <*> some (
            (,) <$> strArgument
                    (  metavar "IMAGE"
                    <> help "input file (can be CALIBRATION)"
                    )
                <*> argument auto
                    (  metavar "DURATION"
                    <> help "duration (in seconds)"
                    )
            )

max_stars_repo_name: DonniK/bSpokeLight

max_issues_repo_name: DonniK/bSpokeLight

max_forks_repo_name: DonniK/bSpokeLight


text:
module LogisticRegression where
import Util
import GradientDescent
import Data.List
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Util
import Numeric.GSL.Minimization

{- Trains multiple logistic regression classifiers and returns all
  the classifiers in a list of matrices [theta], where the i-th element
  corresponds to the classifier for label i -}
classifyOneVsAll :: Matrix Double -> Matrix Double -> Double -> Double -> Int -> ([Matrix Double],[Double])
classifyOneVsAll x y lambda precision max_iters =  (map (\ target -> trainClassifier x target lambda precision max_iters) targetMatrixList,targets)
  where
    (targetMatrix,targets) = mapOutputsVectorToTargetMatrix y
    targetMatrixList = map asColumn $ toColumns targetMatrix

{- Trains a logistic regression classifier, returning the learned theta -}
trainClassifier :: Matrix Double -> Matrix Double -> Double -> Double -> Int -> Matrix Double
trainClassifier x y lambda precision max_iters = asColumn $ fst $ minimizeVDGeneric VectorBFGS2 x y (zeros (cols x) 1) max_iters lambda precision

{- The cost function J
   Compute cost and gradient for logistic regression with regularization.
   Computes the cost of using theta as the parameter for regularized logistic regression and the
   gradient of the cost w.r.t. to the parameters.
   Returns double of cost and gradient -}
computeCost :: Matrix Double -> Matrix Double -> Matrix Double -> Double -> (Double, Matrix Double)
computeCost x y theta lambda = (cost, gradient)
    where
      m = fromIntegral $ rows y
      prediction = sigmoid $ x <> theta
      tempTheta = accum theta (*) [((0,0),0)]
      cost = - (1/m) * sumElements ( y * (log prediction) + (1 - y) * (log (1 - prediction)) ) + (lambda/(2*m)) * sumElements (tempTheta ^ 2)
      gradient = ctrans ((scalar (1/m)) * ( sumColumns  ((prediction - y) * x) ) ) + (scalar (lambda/m)) * tempTheta

{- Variants that only compute cost OR gradient (for min functions) -}
costOnly :: Matrix Double -> Matrix Double -> Matrix Double -> Double -> Double
costOnly x y theta lambda = - (1/m) * sumElements ( y * (log prediction) + (1 - y) * (log (1 - prediction)) ) + (lambda/(2*m)) * sumElements (tempTheta ^ 2)
  where
      m = fromIntegral $ rows y
      prediction = sigmoid $ x <> theta
      tempTheta = accum theta (*) [((0,0),0)]

gradientOnly :: Matrix Double -> Matrix Double -> Matrix Double -> Double -> Matrix Double
gradientOnly x y theta lambda = ctrans ((scalar (1/m)) * ( sumColumns  ((prediction - y) * x) ) ) + (scalar (lambda/m)) * tempTheta
  where
      m = fromIntegral $ rows y
      prediction = sigmoid $ x <> theta
      tempTheta = accum theta (*) [((0,0),0)]

-- Minimize with gradient descent
minimizeGD x y initial_theta alpha num_iters lambda = gradientDescent gradFunc initial_theta alpha num_iters
    where
      gradFunc t = snd (computeCost x y t lambda)

{-
Minimization methods:
ConjugateFR
ConjugatePR
VectorBFGS
VectorBFGS2
SteepestDescent
-}

minimizeVDGeneric method x y initial_theta num_iters lambda precision = minimizeVD method precision num_iters first_step tol costFunction gradientFunction (flatten initial_theta)
    where
      first_step = precision * 100
      tol = 0.1 -- see http://www.gnu.org/software/gsl/manual/html_node/Multimin-Algorithms-with-Derivatives.html
      costFunction thetaVect = costOnly x y (asColumn thetaVect) lambda
      gradientFunction thetaVect = flatten $ gradientOnly x y (asColumn thetaVect) lambda

-- Predict probability that the label is 0 or 1 using learned logistic
predict xi theta = sigmoid (xi <> theta)

-- Predict whether the label is 0 or 1 using learned logistic
predictThreshold :: Matrix Double -> Matrix Double -> Double -> Matrix Double
predictThreshold xi theta threshold = cond (predict xi theta) (scalar threshold) 0 1 1

-- Gauge the accuracy of a prediction
accuracy p y = mean $ toList $ flatten (cond p y 0 1 0) * 100

{- Predict the label for a trained one-vs-all classifier. The labels
   are in the range 1..K, where K = size(all_theta, 1).
   p = PREDICTONEVSALL(all_theta, X) will return a vector of predictions
   for each example in the matrix X. Note that X contains the examples in
   rows. all_theta is a matrix where the i-th row is a trained logistic
   regression theta vector for the i-th class. You should set p to a vector
   of values from 1..K (e.g., p = [1; 3; 1; 2] predicts classes 1, 3, 1, 2
   for 4 examples) -}
predictOneVsAll :: [Matrix Double] -> Matrix Double -> [Double] -> Matrix Double
predictOneVsAll all_theta x targetMapping = selectTarget
  where
    at_matrix = fromBlocks [all_theta]
    maxIndex xs = snd $ maximum $ zip xs targetMapping
    selectTarget = asColumn $ fromList $ map maxIndex $ toLists $ predict x at_matrix

max_stars_repo_name: thade/haskellml

max_issues_repo_name: thade/haskellml

max_forks_repo_name: thade/haskellml


text:
module ISRES where

import Numeric.LinearAlgebra
import Numeric.NLOPT

isres01 :: Int -> Either Result Solution
isres01 n =
  let objf x = x `dot` x + 22 -- define objective
      stop = ObjectiveRelativeTolerance 1e-12 :| [] -- define stopping criterion
      algorithm = ISRES objf [] [] (SeedValue 22) -- specify algorithm
      lowerbounds = fromList $ replicate n (-10) -- specify bounds
      upperbounds = fromList $ replicate n (10) -- specify bounds
      problem = GlobalProblem lowerbounds upperbounds stop algorithm
      x0 = fromList $ replicate n 5 -- specify initial guess
   in minimizeGlobal problem x0

isres02 :: Int -> Either Result Solution
isres02 n =
  let objf x = - norm_1 x
      stop = ObjectiveRelativeTolerance 1e-12 :| [] -- define stopping criterion
      algorithm = ISRES objf [] [eqCon] (SeedValue 22) -- specify algorithm
      lowerbounds = fromList $ replicate n 0 -- specify bounds
      upperbounds = fromList $ replicate n 10 -- specify bounds
      eqCon =
        EqualityConstraint
          { eqConstraintFunctions = Scalar (\x -> x ! 1 * (x ! 0 - 5)),
            {-               Scalar
                            ( \x ->
                                let x1 = subVector 0 (n `div` 2) x
                                    x2 = subVector ((n `div` 2) + 1) n x
                                 in x1 `dot` (x2 - upperbounds)
                            ),
             -}
            eqConstraintTolerance = 1e-12
          }

      problem = GlobalProblem lowerbounds upperbounds stop algorithm
      x0 = fromList $ replicate n 5 -- specify initial guess
   in minimizeGlobal problem x0

isres03 :: Int -> Either Result Solution
isres03 n =
  if n < 3
    then error "n must be at least 3"
    else
      let objf x = - norm_1 x
          stop = ObjectiveRelativeTolerance 1e-12 :| [] -- define stopping criterion
          algorithm = ISRES objf [] [eqCon] (SeedValue 22) -- specify algorithm
          lowerbounds = fromList $ replicate n 0 -- specify bounds
          upperbounds = fromList $ replicate n 10 -- specify bounds
          eqCon =
            EqualityConstraint
              { eqConstraintFunctions =
                  Vector
                    2
                    ( \x _ ->
                        fromList
                          [ x ! 1 * (x ! 0 - 10),
                            x ! 2 * (x ! 1 - (10 - x ! 0))
                          ]
                    ),
                eqConstraintTolerance = 1e-12
              }

          problem = GlobalProblem lowerbounds upperbounds stop algorithm
          x0 = fromList $ replicate n 0 -- specify initial guess
       in minimizeGlobal problem x0

{-
See https://www.mat.univie.ac.at/~neum/glopt/coconut/Benchmark/Library1_new_v1.html

Alkyl problem in AMPL format:
var x2 := 1.745, >= 0, <= 2;
var x3 := 1.2, >= 0, <= 1.6;
var x4 := 1.1, >= 0, <= 1.2;
var x5 := 3.048, >= 0, <= 5;
var x6 := 1.974, >= 0, <= 2;
var x7 := 0.893, >= 0.85, <= 0.93;
var x8 := 0.928, >= 0.9, <= 0.95;
var x9 := 8, >= 3, <= 12;
var x10 := 3.6, >= 1.2, <= 4;
var x11 := 1.45, >= 1.45, <= 1.62;
var x12 := 1, >= 0.99, <= 1.01010101010101;
var x13 := 1, >= 0.99, <= 1.01010101010101;
var x14 := 1, >= 0.9, <= 1.11111111111111;
var x15 := 1, >= 0.99, <= 1.01010101010101;

minimize obj:  - 6.3*x5*x8 + 5.04*x2 + 0.35*x3 + x4 + 3.36*x6;

subject to
e2:  - 0.819672131147541*x2 + x5 - 0.819672131147541*x6 = 0;
e3: 0.98*x4 - x7*(0.01*x5*x10 + x4) = 0;
e4:  - x2*x9 + 10*x3 + x6 = 0;
e5: x5*x12 - x2*(1.12 + 0.13167*x9 - 0.0067*x9*x9) = 0;
e6: x8*x13 - 0.01*(1.098*x9 - 0.038*x9*x9) - 0.325*x7 - 0.57425 = 0;
e7: x10*x14 + 22.2*x11 - 35.82 = 0;
e8: x11*x15 - 3*x8 + 1.33 = 0;
-}
alkyl :: [Double] -> Either Result Solution
alkyl x0 =
  let lb = fromList [0, 0, 0, 0, 0, 0.85, 0.9, 3, 1.2, 1.45, 0.99, 0.99, 0.9, 0.99]
      ub = fromList [2, 1.6, 1.2, 5, 2, 0.93, 0.95, 12, 4, 1.62, 1.01010101010101, 1.01010101010101, 1.11111111111111, 1.01010101010101]
      objf x = - 6.3 * x ! 3 * x ! 6 + 5.04 * x ! 0 + 0.35 * x ! 1 + x ! 2 + 3.36 * x ! 4
      stop = ObjectiveRelativeTolerance 1e-12 :| []
      eqCons =
        [ EqualityConstraint
            { eqConstraintFunctions = Scalar (\x -> -0.819672131147541 * x ! 0 + x ! 3 - 0.819672131147541 * x ! 4),
              eqConstraintTolerance = 1e-10
            },
          EqualityConstraint
            { eqConstraintFunctions = Scalar (\x -> 0.98 * x ! 2 - x ! 5 * (0.01 * x ! 3 * x ! 8 + x ! 2)),
              eqConstraintTolerance = 1e-10
            },
          EqualityConstraint
            { eqConstraintFunctions = Scalar (\x -> - x ! 0 * x ! 7 + 10 * x ! 1 + x ! 4),
              eqConstraintTolerance = 1e-10
            },
          EqualityConstraint
            { eqConstraintFunctions = Scalar (\x -> x ! 3 * x ! 10 - x ! 0 * (1.12 + 0.13167 * x ! 7 - 0.0067 * x ! 7 * x ! 7)),
              eqConstraintTolerance = 1e-10
            },
          EqualityConstraint
            { eqConstraintFunctions = Scalar (\x -> x ! 6 * x ! 11 - 0.01 * (1.098 * x ! 7 - 0.038 * x ! 7 * x ! 7) - 0.325 * x ! 5 - 0.57425),
              eqConstraintTolerance = 1e-10
            },
          EqualityConstraint
            { eqConstraintFunctions = Scalar (\x -> x ! 8 * x ! 12 + 22.2 * x ! 9 - 35.82),
              eqConstraintTolerance = 1e-10
            },
          EqualityConstraint
            { eqConstraintFunctions = Scalar (\x -> x ! 9 * x ! 13 - 3 * x ! 6 + 1.33),
              eqConstraintTolerance = 1e-10
            }
        ]
      algorithm = ISRES objf [] eqCons (SeedValue 22)
      problem = GlobalProblem lb ub stop algorithm
   in minimizeGlobal problem (fromList x0)


text:
module Statistics.Report
    ( reportTxFate
    ) where

import           Data.Map.Strict (Map)
import qualified Data.Map.Strict as M
import qualified Data.Set as S
import qualified Data.Text as T
import           System.IO (hPutStrLn)

import           Statistics.Block (TxFate (..))
import           Types
import           Universum

reportTxFate :: FilePath -> Map TxHash TxFate -> IO ([TxHash], [TxHash], [TxHash])
reportTxFate f m = withFile f WriteMode $ \h -> do
    let total = M.size m
        m'    = M.filter isFork m
        fork  = M.size m'
        m''   = M.filter isChain m
        chain = M.size m''
        m'''  = M.filter isLost m
        lost  = M.size m'''
        chains = fst <$> M.toList m''
        forks  = fst <$> M.toList m'
        losts  = fst <$> M.toList m''
    hPutStrLn h $ "total transactions: " ++ show total
    hPutStrLn h $ "in blockchain: " ++ show chain
    hPutStrLn h $ "in fork: " ++ show fork
    hPutStrLn h $ "lost: " ++ show lost
    hPutStrLn h ""
    hPutStrLn h "transactions in blockchain:"
    for_ (M.toList m'') $ \(tx, InBlockChain _ hash _) ->
        hPutStrLn h $ toString tx ++ ": " ++ toString (T.take 6 hash)
    hPutStrLn h ""
    hPutStrLn h "transactions in fork:"
    for_ (M.toList m') $ \(tx, InFork s) ->
        hPutStrLn h $ toString tx ++ ": " ++ toString (unwords $ map (T.take 6 . snd) $ S.toList s)
    hPutStrLn h ""
    hPutStrLn h "lost transactions:"
    for_ (M.toList m''') $ \(tx, InNoBlock) ->
        hPutStrLn h (toString tx)
    return (chains, forks, losts)

isFork :: TxFate -> Bool
isFork (InFork _) = True
isFork _          = False

isLost :: TxFate -> Bool
isLost InNoBlock = True
isLost _         = False

isChain :: TxFate -> Bool
isChain (InBlockChain{}) = True
isChain _                = False

max_stars_repo_name: jasper-perf/cardano-sl


text:
{-# LANGUAGE NamedFieldPuns, TemplateHaskell #-}

module School.Unit.Test.MultiNoulli
( multiNoulliTest ) where

import Conduit ((.|), await, liftIO, yield)
import Data.Default.Class (def)
import Numeric.LinearAlgebra ((><))
import School.TestUtils (addClasses, assertRight, diffCost, randomMatrix,
                         randomNNInts, singleClassOutput, testState, unitCorrect)
import School.Train.AppTrain (AppTrain)
import School.Train.ForwardPass (forwardPass)
import School.Types.FloatEq ((~=))
import School.Types.Slinky (Slinky(..), slinkySingleton)
import School.Unit.CostFunction (CostFunction(..))
import School.Unit.CostParams (CostParams(..))
import School.Unit.LogSoftMax (logSoftMax)
import School.Unit.MultiNoulli
import School.Unit.Unit (Unit(..))
import School.Unit.UnitActivation (UnitActivation(..))
import School.Unit.UnitGradient (UnitGradient(..))
import School.Unit.UnitParams (UnitParams(..))
import School.Utils.LinearAlgebra (compareDoubleMatrix)
import School.Utils.Tuple (trd3)
import Test.Tasty (TestTree)
import Test.Tasty.QuickCheck hiding ((><))
import Test.Tasty.TH
import Test.QuickCheck.Monadic (monadicIO, assert, pre)

eps :: Double
eps = 1e-5

prec :: Double
prec = 5e-2

multi :: CostFunction Double (AppTrain a)
multi = multiNoulli Nothing Nothing

prop_correct_classes :: (Positive Int) -> (Positive Int) -> Property
prop_correct_classes (Positive c) (Positive b) = monadicIO $ do
  classes <- liftIO $ randomNNInts (c - 1) b
  let out = singleClassOutput c classes
  let params = slinkySingleton $ BatchClassTarget classes
  let cost = computeCost multi out params
  assert $ cost ~= Right (-1)

prop_wrong_classes :: (Positive Int) -> (Positive Int) -> Property
prop_wrong_classes (Positive c) (Positive b) = monadicIO $ do
  pre $ c > 1
  classes <- liftIO $ randomNNInts (c - 1) b
  let wrong = map (\k -> mod (k + 1) c) classes
  let out = singleClassOutput c wrong
  let params = slinkySingleton $ BatchClassTarget classes
  let cost = computeCost multi out params
  assert $ cost == Right 0

prop_norm :: (Positive Int) -> (Positive Int) -> Property
prop_norm (Positive c) (Positive b) = monadicIO $ do
  classes <- liftIO $ randomNNInts (c - 1) b
  input <- liftIO $ BatchActivation <$> randomMatrix b c
  let out = apply logSoftMax EmptyParams input
  let params = slinkySingleton $ BatchClassTarget classes
  let cost = computeCost multi out params
  assertRight (>= (-1)) cost

prop_numerical_deriv :: (Positive Int) -> (Positive Int) -> Property
prop_numerical_deriv (Positive c) (Positive b) = monadicIO $ do
  pre $ b < 23 && c < 23
  classes <- liftIO $ slinkySingleton . BatchClassTarget <$> randomNNInts (c - 1) b
  input <- liftIO $ BatchActivation <$> randomMatrix b c
  let deriv = derivCost multi input classes
  let check = (b >< c) [ diffCost multi
                                  input
                                  eps
                                  (j, k)
                                  classes
                               | j <- [0..b-1], k <- [0..c-1] ]
  assertRight (\(BatchGradient g) -> compareDoubleMatrix prec g check)
              deriv

prop_correct_forward_pass :: (Positive Int) -> (Positive Int) -> Property
prop_correct_forward_pass (Positive c) (Positive b) = monadicIO $ do
  classes <- liftIO $ randomNNInts (c - 1) b
  let forward = forwardPass [unitCorrect classes] multi
  activation <- liftIO $ randomMatrix b c
  let input = ([BatchActivation $ addClasses classes activation], SNil)
  let pass = yield input
          .| forward
          .| await
  result <- testState pass def
  assertRight ((maybe False ((~= (-1)) . trd3)) . fst)
              result

multiNoulliTest :: TestTree
multiNoulliTest = $(testGroupGenerator)

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


text:
{-# LANGUAGE FlexibleContexts #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Statistics.Information
-- Copyright   :  (c) A. V. H. McPhail 2010, 2014
-- License     :  BSD3
--
-- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
-- Stability   :  provisional
-- Portability :  portable
--
-- Shannon entropy
--
-----------------------------------------------------------------------------

module Numeric.Statistics.Information (
                                   entropy
                                   , mutual_information
                                  ) where

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

import Numeric.Statistics.PDF

import qualified Numeric.LinearAlgebra as LA
--import Numeric.LinearAlgebra.Data hiding(Vector)

--import qualified Data.Vector as DV
import qualified Data.Vector.Storable as V

import Prelude hiding(map,zip)

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

zeroToOne x
    | x == 0.0  = 1.0
    | otherwise = x

logE = V.map (log . zeroToOne)

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

-- | the entropy \sum p_i l\ln{p_i} of a sequence
entropy :: PDF a Double
        => a                       -- ^ the underlying distribution
        -> LA.Vector Double           -- ^ the sequence
        -> Double                  -- ^ the entropy
entropy p x = let ps = probability p x
              in negate $ (LA.dot ps (logE ps))

-- | the mutual information \sum_x \sum_y p(x,y) \ln{\frac{p(x,y)}{p(x)p(y)}}
mutual_information :: (PDF a Double, PDF b (Double,Double))
                   => b                                          -- ^ the underlying distribution
                   -> a                                          -- ^ the first dimension distribution
                   -> a                                          -- ^ the second dimension distribution
                   -> (LA.Vector Double, LA.Vector Double)             -- ^ the sequence
                   -> Double         -- ^ the mutual information
mutual_information p px py (x,y) = let ps = probability p $ V.zipWith (,) x y
                                       xs = probability px x
                                       ys = probability py y
                                   in (LA.dot ps (logE ps - logE (xs*ys)))

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

max_stars_repo_name: amcphail/hstatistics

max_issues_repo_name: amcphail/hstatistics

max_forks_repo_name: amcphail/hstatistics


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

module Audio.FFT.Inplace where

import Audio.Wav.Sample
import Conduit
import Control.Arrow ((>>>))
import qualified Control.Foldl as L
import Control.Monad (forM_, (<$!>))
import Control.Monad.Loops (iterateUntilM)
import Data.Bits (Bits (bit, shiftR))
import Data.Complex (Complex ((:+)), magnitude)
import Data.Conduit.Utils (chunkedVector)
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Generic.Mutable as MG
import qualified Data.Vector.Unboxed.Mutable as MU
import Math.NumberTheory.Logarithms
import Streaming
import qualified Streaming as S
import qualified Streaming.Prelude as S

fftC ::
  (Monad m, G.Vector v (Complex Double), G.Vector v Double) =>
  Int ->
  ConduitM Double (v Double) m ()
fftC n = chunkedVector n .| mapC (G.map magnitude . simpleFFT . G.map (:+ 0))

fftS ::
  (Monad m, G.Vector v (Complex Double), G.Vector v Double, PrimMonad m) =>
  Int ->
  Stream (Of Sample) m () ->
  Stream (Of (v Double)) m ()
{-# INLINE fftS #-}
fftS n =
  S.chunksOf n
    >>> S.mapped (L.impurely S.foldM L.vectorM)
    >>> S.map (G.map magnitude . simpleFFT . G.map (:+ 0))

simpleFFT :: G.Vector v (Complex Double) => v (Complex Double) -> v (Complex Double)
{-# INLINEABLE simpleFFT #-}
simpleFFT = G.modify $ \inps -> do
  reverseBitM inps
  let !theta = 2 * pi / fromIntegral (MG.length inps)
  go (cos theta) (sin theta) inps
  pure ()
  where
    {-# INLINE go #-}
    go !c !s inps
      | n <= 1 = pure ()
      | otherwise = do
        let half = n `div` 2
            (lh, uh) = MG.splitAt half inps
            !dblCos = 2 * c * c - 1
            !dblSin = 2 * s * c
            !w = c :+ s
        go dblCos dblSin lh
        go dblCos dblSin uh
        forM_ [0 .. half - 1] $ \k -> do
          !ek <- MG.read lh k
          !ok <- MG.read uh k
          MG.write inps k $ ek + w ^ k * ok
          MG.write inps (half + k) $! ek + w ^ (half + k) * ok
      where
        !n = MG.length inps

reverseBit :: (G.Vector v a) => v a -> v a
{-# INLINE reverseBit #-}
reverseBit = G.modify reverseBitM

reverseBitM :: (MG.MVector v a, PrimMonad m) => v (PrimState m) a -> m ()
{-# INLINEABLE reverseBitM #-}
reverseBitM vec = do
  let !n = intLog2 $ MG.length vec
      !m = bit (n `shiftR` 1)
  table <- MU.replicate m 0
  MU.write table 0 0

  void $
    iterateUntilM
      (\(!pk, !pl) -> pl + 1 >= pk)
      ( \(!pk, !pl) -> do
          let !k = bit (pk - 1)
              !l = bit pl
          forM_ [0 .. l - 1] $ \j ->
            MU.write table (l + j) . (+ k) =<< MU.read table j
          pure (pk - 1, pl + 1)
      )
      (n, 0)

  if even n
    then forM_ [1 .. m - 1] $ \i -> forM_ [0 .. i - 1] $ \j -> do
      !ji <- (j +) <$!> MU.read table i
      !ij <- (i +) <$!> MU.read table j
      MG.swap vec ji ij
    else forM_ [1 .. m - 1] $ \i -> forM_ [0 .. i - 1] $ \j -> do
      !ji <- (j +) <$!> MU.read table i
      !ij <- (i +) <$!> MU.read table j
      MG.swap vec ji ij
      MG.swap vec (ji + m) (ij + m)

max_stars_repo_name: konn/genetic-twelve

max_issues_repo_name: konn/genetic-twelve

max_forks_repo_name: konn/genetic-twelve


text:
module LeastSquares where

import Numeric.LinearAlgebra

n = 20
x = linspace n (-3,7::Double)
y0 = 3 * x

main = do
        noise <- randn 1 n
        let y = (flatten noise) + y0
        let sampleMatrix = (asColumn x) ||| (konst 1 (n,1))
        let sol = (sampleMatrix <\> y)
        print $ "Best fit is y = " ++ show (sol ! 0) ++ " * x + " ++ (show (sol ! 1))


text:
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}

{-# LANGUAGE AllowAmbiguousTypes   #-}
{-# LANGUAGE CPP                   #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DeriveAnyClass        #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedStrings     #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE UndecidableInstances  #-}
{-|
Module      : Grenade.Layers.SamePadPooling
Description : Pooling layer that supports SAME_UPPER padding
Maintainer  : Theo Charalambous
License     : BSD2
Stability   : experimental
-}
module Grenade.Layers.SamePadPooling (
  -- * Layer Definition
    SamePadPooling (..)
  ) where

import           Control.DeepSeq
import           Data.Function                   ((&))
import           Data.Kind                       (Type)
import           Data.Maybe
import           Data.Proxy
import           Data.Serialize
import           Data.Singletons.TypeLits        hiding (natVal)
import           GHC.Generics
import           GHC.TypeLits
import           Numeric.LinearAlgebra.Static    (create, extract)

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

-- | A pooling layer for a neural network, for when auto_pad is SAME_UPPER or SAME_LOWER
--
--   Pads on the X and Y dimension of an image.
data SamePadPooling  :: Nat -> Nat -> Nat -> Nat -> Nat -> Nat -> Nat -> Nat -> Type where
  SamePadPooling :: SamePadPooling kernelRows kernelColumns strideRows strideColumns padLeft padTop padRight padBottom
  deriving (NFData, Generic)

instance Show (SamePadPooling kernelRows kernelColumns strideRows strideColumns padLeft padTop padRight padBottom) where
  show SamePadPooling = "SamePadPooling"

instance UpdateLayer (SamePadPooling kernelRows kernelColumns strideRows strideColumns padLeft padTop padRight padBottom) where
  type Gradient (SamePadPooling kernelRows kernelColumns strideRows strideColumns padLeft padTop padRight padBottom) = ()
  runUpdate _ x _ = x
  reduceGradient _ = ()

instance RandomLayer (SamePadPooling kernelRows kernelColumns strideRows strideColumns padLeft padTop padRight padBottom)  where
  createRandomWith _ _ = return SamePadPooling

instance Serialize (SamePadPooling kernelRows kernelColumns strideRows strideColumns padLeft padTop padRight padBottom) where
  put _ = return ()
  get = return SamePadPooling

-- | A two dimentional image can be padded.
instance ( KnownNat padLeft
         , KnownNat padTop
         , KnownNat padRight
         , KnownNat padBottom
         , KnownNat strideRows
         , KnownNat strideColumns
         , KnownNat kernelRows
         , KnownNat kernelColumns
         , KnownNat inputRows
         , KnownNat inputColumns
         , KnownNat outputRows
         , KnownNat outputColumns
         , (padLeft + padRight) ~ ((outputColumns - 1) * strideColumns + kernelColumns - inputColumns)
         , (padTop + padBottom) ~ ((outputRows - 1) * strideRows + kernelRows - inputRows)
         ) => Layer (SamePadPooling kernelRows kernelColumns strideRows strideColumns padLeft padTop padRight padBottom) ('D2 inputRows inputColumns) ('D2 outputRows outputColumns) where
  type Tape (SamePadPooling kernelRows kernelColumns strideRows strideColumns padLeft padTop padRight padBottom) ('D2 inputRows inputColumns) ('D2 outputRows outputColumns)  = ()
  runForwards SamePadPooling (S2D input) =
    let padl  = fromIntegral $ natVal (Proxy :: Proxy padLeft)
        padt  = fromIntegral $ natVal (Proxy :: Proxy padTop)
        padr  = fromIntegral $ natVal (Proxy :: Proxy padRight)
        padb  = fromIntegral $ natVal (Proxy :: Proxy padBottom)
        kr    = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        kc    = fromIntegral $ natVal (Proxy :: Proxy kernelColumns)
        sr    = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sc    = fromIntegral $ natVal (Proxy :: Proxy strideColumns)
        h     = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        w     = fromIntegral $ natVal (Proxy :: Proxy inputColumns)

        m     = extract input

        r     = validPadPoolForwards 1 h w kr kc sr sc padl padt padr padb m

    in  ((), S2D . fromJust . create $ r)
  runBackwards _ _ _ = error "backward pass for SamePadPooling not implemented"

-- | A two dimentional image can be padded.
instance ( KnownNat padLeft
         , KnownNat padTop
         , KnownNat padRight
         , KnownNat padBottom
         , KnownNat strideRows
         , KnownNat strideColumns
         , KnownNat kernelRows
         , KnownNat kernelColumns
         , KnownNat inputRows
         , KnownNat inputColumns
         , KnownNat outputRows
         , KnownNat outputColumns
         , KnownNat channels
         , KnownNat (inputRows * channels)
         , KnownNat (outputRows * channels)
         , (padLeft + padRight) ~ ((outputColumns - 1) * strideColumns + kernelColumns - inputColumns)
         , (padTop + padBottom) ~ ((outputRows - 1) * strideRows + kernelRows - inputRows)
         ) => Layer (SamePadPooling kernelRows kernelColumns strideRows strideColumns padLeft padTop padRight padBottom) ('D3 inputRows inputColumns channels) ('D3 outputRows outputColumns channels) where
  type Tape (SamePadPooling kernelRows kernelColumns strideRows strideColumns padLeft padTop padRight padBottom) ('D3 inputRows inputColumns channels) ('D3 outputRows outputColumns channels)  = ()
  runForwards SamePadPooling (S3D input) =
    let padl  = fromIntegral $ natVal (Proxy :: Proxy padLeft)
        padt  = fromIntegral $ natVal (Proxy :: Proxy padTop)
        padr  = fromIntegral $ natVal (Proxy :: Proxy padRight)
        padb  = fromIntegral $ natVal (Proxy :: Proxy padBottom)
        kr    = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        kc    = fromIntegral $ natVal (Proxy :: Proxy kernelColumns)
        sr    = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sc    = fromIntegral $ natVal (Proxy :: Proxy strideColumns)
        h     = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        w     = fromIntegral $ natVal (Proxy :: Proxy inputColumns)
        c     = fromIntegral $ natVal (Proxy :: Proxy channels)

        m     = extract input

        r     = validPadPoolForwards c h w kr kc sr sc padl padt padr padb m
    in  ((), S3D . fromJust . create $ r)

  runBackwards _ _ _ = error "backward pass for SamePadPooling not implemented"

instance OnnxOperator (SamePadPooling kernelRows kernelColumns strideRows strideColumns padLeft padTop padRight padBottom) where
  onnxOpTypeNames _ = ["MaxPool"]

instance ( KnownNat padLeft
         , KnownNat padTop
         , KnownNat padRight
         , KnownNat padBottom
         , KnownNat strideRows
         , KnownNat strideColumns
         , KnownNat kernelRows
         , KnownNat kernelColumns
         ) => OnnxLoadable (SamePadPooling kernelRows kernelColumns strideRows strideColumns padLeft padTop padRight padBottom) where

  loadOnnxNode _ node = do
    node & hasSupportedDilations

    (node `hasMatchingShape` "kernel_shape") kernelShape
    (node `hasMatchingShape` "strides")      strideShape

    -- todo: check that attribute is one of: SAME_UPPER or SAME_LOWER

    return SamePadPooling
      where
        kernelShape = [natVal (Proxy :: Proxy kernelRows), natVal (Proxy :: Proxy kernelColumns)]
        strideShape = [natVal (Proxy :: Proxy strideRows), natVal (Proxy :: Proxy strideColumns)]

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


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

import CV.ColourUtils
import CV.DFT
import CV.Drawing
import CV.Edges
import CV.Filters
import CV.Image
import CV.ImageOp
import CV.Pixelwise
import Control.Applicative hiding ((<**>))
import qualified CV.ImageMath as IM
import qualified CV.Transforms as T
import Data.Complex

fi = fromIntegral

main = do
    Just x <- loadImage "smallLena.jpg"
    let  d :: Pixelwise (Complex D32)
         d  = fromImage $ dft x
         mpper f = \(x,y) -> f (x,y) + (sin (7*fi x))
         d2' = remapImage mpper x
         d2  = dft d2'
         d2c = d2 <# circleOp (0:+0) (130,102) 13 Filled
                  <# circleOp (0:+0) (71,102)  13 Filled

         dftd1 = (toImage d) <# circleOp (0:+0) (101,101) 30 Filled
    saveImage "dftmanip.png" $ montage (2,3) 10
        [
         logarithmicCompression $ fst $ dftSplit $ dftd1
        ,stretchHistogram $ idft dftd1
        ,stretchHistogram $ d2'
        ,logarithmicCompression $ fst $ dftSplit $ d2
        ,stretchHistogram $ fst $ dftSplit $ d2c
        ,stretchHistogram $ idft $ d2c
        ]

max_stars_repo_name: maaleske/CV

max_issues_repo_name: maaleske/CV

max_forks_repo_name: maaleske/CV


text:
module SentimentAnalysisSpec where

import Numeric.LinearAlgebra as LA
import SentimentAnalysis
import Test.Hspec

theta :: Vector R
theta = LA.vector [-1.0, 1.0]

spec :: Spec
spec =
  describe "Hinge Loss Tests" $ do

    it "Hinge loss single" $ do
      let x = LA.vector [1, 2]
          y = 1
          theta0 = -0.2
          m = hingeLossSingle x y theta theta0
      fromInteger (round (m * 10)) / 10 `shouldBe` 0.2

    it "Hinge loss full" $ do
      let x = (2><2) [1.0, 2.0, 1.0, 2.0]
          y = LA.vector [1.0, 1.0]
          theta0 = -0.2
          m = hingeLossFull x y theta theta0
      fromInteger (round (m * 10)) / 10 `shouldBe` 0.2

    it "Perceptron single update" $ do
        let x = LA.vector [1, 2]
            y = 1
            theta0 = -1.5
        perceptronSingleStepUpdate x y theta theta0 `shouldBe` (LA.vector [0, 3], -0.5)

    it "Perceptron 1" $ do
        let x1 = (1><2) [1.0, 2.0]
            y = LA.vector [1]
            t = 1
        perceptron x1 y t [0] `shouldBe` (LA.vector [1, 2], 1)

    it "Perceptron 2" $ do
        let x = (2><2) [1.0, 2.0, -1.0, 0]
            y = LA.vector [1, 1]
            t = 1
        perceptron x y t [0, 1] `shouldBe` (LA.vector [0, 2], 2)

    it "Perceptron 3" $ do
        let x3 = (1><2) [1.0, 2.0]
            y = LA.vector [1]
            t = 2
        perceptron x3 y t [0] `shouldBe` (LA.vector [1, 2], 1)

    it "Perceptron 4" $ do
        let x = (2><2) [1.0, 2.0, -1.0, 0]
            y = LA.vector [1, 1]
            t = 2
        perceptron x y t [0,1] `shouldBe` (LA.vector [0, 2], 2)

    it "Average Perceptron 1" $ do
        let x = (1><2) [1.0, 2.0]
            y = LA.vector [1]
            t = 1
        averagePerceptron x y t [0] `shouldBe` (LA.vector [1, 2], 1)

    it "Average Perceptron 2" $ do
        let x = (2><2) [-1.0, 0, 1.0, 2.0]
            y = LA.vector [1, 1]
            t = 1
        averagePerceptron x y t [0, 1] `shouldBe` (LA.vector [-0.5, 1], 1.5)

    it "Average Perceptron 3" $ do
        let x = (1><2) [1.0, 2.0]
            y = LA.vector [1]
            t = 2
        averagePerceptron x y t [0] `shouldBe` (LA.vector [1, 2], 1)

    it "Average Perceptron 4" $ do
        let x = (2><2) [-1.0, 0, 1.0, 2.0] :: Matrix Double
            y = LA.vector [1, 1]
            t = 2
        averagePerceptron x y t [0, 1] `shouldBe` (LA.vector [-0.25, 1.5], 1.75)

    it "Pegasos single update 1" $ do
        let x4 = LA.vector [1, 2]
            y1 = 1
            lambda = 0.2
            eta = 0.1
            theta0 = -1.5
        pegasosSingleStepUpdate x4 y1 lambda eta theta theta0 `shouldBe` (LA.vector [-0.88, 1.18], -1.4)

    it "Pegasos single update 2" $ do
        let x = LA.vector [1, 1]
            y2 = 1
            lambda = 0.2
            eta = 0.1
            theta0 = 1
        pegasosSingleStepUpdate x y2 lambda eta theta theta0 `shouldBe` (LA.vector [-0.88, 1.08], 1.1)

    it "Pegasos single update 3" $ do
        let x = LA.vector [1, 2]
            y = 1
            lambda = 0.2
            eta = 0.1
            theta0 = -2
        pegasosSingleStepUpdate x y lambda eta theta theta0 `shouldBe` (LA.vector [-0.88, 1.18], -1.9)

    it "Pegasos 1" $ do
        let x = (1><2) [1, 2]
            y = LA.vector [1]
            t = 1
            lambda = 0.2
        pegasos x y t lambda [0] `shouldBe` (LA.vector [1, 2], 1)

    it "Pegasos 2" $ do
        let x = (2><2) [1, 1, 1, 1]
            y = LA.vector [1, 1]
            t = 1
            lambda = 1
        pegasos x y t lambda [0, 1] `shouldBe` (LA.vector [1 - 1 / sqrt 2, 1 - 1 / sqrt 2], 1)


text:
---------------------------------------------------------------------------
-- | Module    : Numeric.Optimization.Basic.LP
-- Copyright   : (c) 2016 Marco Zocca
-- License     : GPL-3
--
-- Maintainer  : zocca.marco gmail
-- Stability   : experimental
-- Portability : portable
--
-- This module implements a general interior-point method for LP, based on the logarithmic barrier penalization
--
-- * [1] Primal-dual interior point algorithms for linear programming, G. Tzallas-Regas
--
--------------------------------------------------------------------------
module Numeric.Optimization.Basic.LP where

import Control.Monad.State

import Data.Sparse -- (reexports Data.Sparse.Common from sparse-linear-algebra)

import Numeric.LinearAlgebra.Sparse

feasX :: (Foldable t, Ord a, Num a) => t a -> Bool
feasX = all (> 0)

linesearch amax x s dx ds = undefined where

-- * Karmarkar algorithm

karmIpLP :: SpMatrix Double -> SpVector Double -> KarmLPData
karmIpLP aa c = execState (untilConverged (\x -> c `dot` _xKarm x) (karmarkarLPstep aa c ar)) kInit
  where
   ar = alpha * r
   kInit = KarmLPData x0 y0 0
   (m, n) = dimSM aa
   n' = fromIntegral n
   r = 1 / sqrt (n' * (n'-1))
   alpha = (n' - 1)/ (3*n')
   x0 = recip (fromIntegral n) .* onesSV n
   y0 = x0

data KarmLPData =
  KarmLPData { _xKarm, _yKarm:: SpVector Double, _niter :: Int } deriving (Eq, Show)

-- | alphar = alpha * r
karmarkarLPstep ::
  SpMatrix Double -> SpVector Double -> Double -> KarmLPData -> KarmLPData
karmarkarLPstep aa c alphar (KarmLPData x y k) = KarmLPData xnew ynew (k + 1) where
  (m, n) = dim aa
  ddk = diagonalSM x
  pp = (aa #~# ddk) -=- svToSM (onesSV n)
  cdk = c <# ddk
  cp = cdk ^-^ (transposeSM pp #> ((pp ##^ pp) <\> (pp #> cdk)))
  cp0 = normalize 2 cp
  ynew = y ^-^ (alphar .* cp0)
  dy = ddk #> ynew
  xnew = recip (sum dy) .* dy
  -- z = c `dot` xnew

-- untilConverged :: MonadState a m => (a -> SpVector Double) -> (a -> a) -> m a
untilConverged fproj = modifyInspectN 200 (normDiffConverged fproj)

normDiffConverged f [x1,x2] = nearZero $ f x1 - f x2

-- normDiffConverged fp xx = nearZero $ normSq (foldrMap fp (^-^) (zeroSV 0) xx)

{- |
Interior point, log-penalized LP (PLP):

max_x <c, x> - rho B(x)

s.t. A x = b

where B(x) := sum_j( log(x_j) )  (i.e. non-negativity constraint)

1. Write Lagrangian for PLP L(x, rho, y)
2. Write KKT conditions (first order optimality of L(x, rho, y))
3. Write Newton equations for KKT system
-}

-- | Update primal `x`, dual `y` and slack `s` variables by applying Newton's method to the KKT system of the log-penalized LP.
-- This requires one linear solve and a line search for a stepsize `alpha` that preserves feasibility : (x, s) > 0.
-- NB : It assumes the dimensions of vectors are x :: Rn , y :: Rm , s :: Rn

-- pdIpLPstep ::
--   (SpVector Double ->
--    SpVector Double ->
--    SpVector Double ->
--    SpVector Double ->
--    Double) ->
--   SpMatrix Double ->
--   Double ->
--   (SpVector Double, SpVector Double, SpVector Double, Int) ->
--   (SpVector Double, SpVector Double, SpVector Double, Int)
-- pdIpLPstep chooseAlpha aa tau (x, y, s, k) = (xNew, yNew, sNew, k+1) where
--   gamma = (x `dot` s) / fromIntegral n  -- duality measure
--   rho = tau * gamma                     -- Lagrange multiplier of barrier penalty
--   ss = diagonalSM s
--   xx = diagonalSM x
--   arow1 = zn -||- transposeSM aa -||- eye n
--   arow2 = aa -||- zm             -||- zm
--   arow3 = ss -||- zm             -||- xx
--   amat = arow1 -=- (arow2 -=- arow3)
--   rhs = concatSV zvn $ concatSV zvm $ negated (xx #> s) ^+^ (rho .* onesSV m)
--   xysStep = amat <\> rhs
--   xStep = takeSV n xysStep            -- 0 .. n-1
--   yStep = takeSV m (dropSV n xysStep) -- n .. n+m-1
--   sStep = dropSV (m + n) xysStep      -- n+m .. 2n+m
--   xNew = x ^+^ (alpha .* xStep)
--   yNew = y ^+^ (alpha .* yStep)
--   sNew = s ^+^ (alpha .* sStep)
--   alpha = chooseAlpha x s xStep sStep
--   (m, n) = dim aa
--   zm = zeroSM m m
--   zn = zeroSM n n
--   zvn = zeroSV n
--   zvm = zeroSV m

max_stars_repo_name: ocramz/optimization-basic

max_issues_repo_name: ocramz/optimization-basic

max_forks_repo_name: ocramz/optimization-basic


text:
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE OverloadedStrings, ScopedTypeVariables, BangPatterns #-}
{-# LANGUAGE DeriveGeneric, DeriveAnyClass, RankNTypes #-}

module Optimizations(
   processVectorsBV
 , fusionOnListFoldl
 , fusionOnListFoldl'
 , fusionOnListFoldr
 , fusedOnListFoldl
 , fusedOnListFoldl'
 , fusionOnStream
 , fusedOnListFoldr
 , fusedOnStream
 , sumOnStream
 , sumOnChan
 , sumOnVector
 , sumOnUnagi
 , sumUsingTasks
             ) where

import Common
import Process

import Data.List as L
import Control.Concurrent.Async
import Control.Concurrent.MVar
import Control.Concurrent (forkIO)
import qualified Data.Vector.Unboxed as BV
import Control.DeepSeq
import Control.Exception.Safe (catch, catchAny, onException, finally, handleAny, bracket
                              , SomeException, throwIO, throw, Exception, MonadMask
                              , withException, displayException)
import Criterion.Main
import Criterion.Main.Options
import Criterion.Types
import Statistics.Types
import Data.Hashable
import System.Random
import qualified Data.Vector.Unboxed as BV
import qualified System.IO.Streams as S
import qualified System.IO.Streams.Combinators as S
import qualified System.IO.Streams.List as S
import qualified System.IO.Streams.Vector as S
import qualified Control.Concurrent.Chan.Unagi as UStream

-- --------------------------------
-- Bench lazy vs strict

-- process-begin
processVectorsBV :: Int -> Int -> Int -> Bool -> Int -> Int
processVectorsBV !seed !vectorDimension !numberOfVectors !useStrictAlgo !lastSum
  = case numberOfVectors of
      0 -> lastSum
      _ -> let bv = randomBV seed vectorDimension
               !h = case useStrictAlgo of
                     True -> slowComputationOnVector bv seed
                     False -> slowComputationOnVector2 bv seed
               !r = processVectorsBV (seed + 1) vectorDimension (numberOfVectors - 1) useStrictAlgo (lastSum + h)
           in r
-- process-end

-- -------------------------------
-- Bench stream fusion

-- fusion-list-begin
fusionOnListFoldl :: [Int] -> Int
fusionOnListFoldl l = L.foldl (+) 0 $ map (* 2) $ filter odd l

fusionOnListFoldl' :: [Int] -> Int
fusionOnListFoldl' l = L.foldl' (+) 0 $ map (* 2) $ filter odd l

fusionOnListFoldr :: [Int] -> Int
fusionOnListFoldr l = L.foldr (+) 0 $ map (* 2) $ filter odd l

fusedCalcOnFoldl :: Int -> Int -> Int
fusedCalcOnFoldl s x
  = case odd x of
      True -> s + (x * 2)
      False -> s
{-# INLINE fusedCalcOnFoldl #-}

fusedCalcOnFoldr :: Int -> Int -> Int
fusedCalcOnFoldr x s = fusedCalcOnFoldl s x
{-# INLINE fusedCalcOnFoldr #-}

fusedOnListFoldl' :: [Int] -> Int
fusedOnListFoldl' l = L.foldl' fusedCalcOnFoldl 0 l

fusedOnListFoldl :: [Int] -> Int
fusedOnListFoldl l = L.foldl fusedCalcOnFoldl 0 l

fusedOnListFoldr :: [Int] -> Int
fusedOnListFoldr l = L.foldr fusedCalcOnFoldr 0 l
-- fusion-list-end

-- fusion-stream-begin
fusionOnStream :: [Int] -> IO Int
fusionOnStream l = do
  inS1 <- S.fromList l
  inS2 <- S.filter odd inS1
  inS3 <- S.map (* 2) inS2
  S.fold (+) 0 inS3

fusedOnStream :: [Int] -> IO Int
fusedOnStream l = do
  inS1 <- S.fromList l
  S.fold fusedCalcOnFoldl 0 inS1
-- fusion-stream-end

-- -----------------------------------
-- Bench channels vs stream vs vector

-- sumOnStream-begin
sumOnStream :: BV.Vector Int -> IO Int
sumOnStream l = do
  inS1 <- S.fromVector l
  S.fold fusedCalcOnFoldl 0 inS1
-- sumOnStream-end

-- sumOnVector-begin
sumOnVector :: BV.Vector Int -> Int
sumOnVector l = BV.foldl' fusedCalcOnFoldl 0 l
-- sumOnVector-end

-- sumOnChan-begin
sumOnChan :: BV.Vector Int -> IO Int
sumOnChan l = do
  ch <- newChan

  (_, r) <- concurrently
              (do BV.mapM_ (writeChan ch) l
                  writeEOFChan ch)
              (sumOnChan' ch 0)
  return r
 where

   sumOnChan' ch !s = do
     mx <- readChan ch
     case mx of
       Nothing -> return s
       Just x -> sumOnChan' ch (fusedCalcOnFoldl s x)
-- sumOnChan-end

-- sumOnUnagi-begin
sumOnUnagi :: BV.Vector Int -> IO Int
sumOnUnagi l = do
  (inCh, outCh) <- UStream.newChan

  (_, r) <- concurrently
              (do BV.mapM_ (\i -> UStream.writeChan inCh (Just i)) l
                  UStream.writeChan inCh Nothing)
              (sumOnChan' outCh 0)
  return r
 where

   sumOnChan' ch !s = do
     mx <- UStream.readChan ch
     case mx of
       Nothing -> return s
       Just x -> sumOnChan' ch (fusedCalcOnFoldl s x)
-- sumOnUnagi-end

-- ---------------------------
-- Threads as Task

-- sumUsingTasks-begin
sumUsingTasks :: BV.Vector Int -> IO Int
sumUsingTasks l = do

  finalResult :: MVar Int <- newEmptyMVar

  currentState :: MVar (Int, Int) <- newEmptyMVar

  nrOfElements <- BV.foldM (\i x -> do forkIO $ sumTask x finalResult currentState
                                       return $ i + 1) 0 l

  putMVar currentState (0, nrOfElements)

  takeMVar finalResult

 where

   sumTask :: Int -> MVar Int -> MVar (Int, Int) -> IO ()
   sumTask x finalResult currentState = do
     (s1, leftSums) <- takeMVar currentState
     let !currentState2@(!s2, !leftSums2) = (fusedCalcOnFoldl s1 x, leftSums - 1)
     case leftSums2 of
       0 -> putMVar finalResult s2
       _ -> putMVar currentState currentState2

-- sumUsingTasks-end

max_stars_repo_name: massimo-zaniboni/threads-post

max_issues_repo_name: massimo-zaniboni/threads-post

max_forks_repo_name: massimo-zaniboni/threads-post


text:
import Debug.Trace
-- import Control.Monad.State.Lazy

import System.Environment (getArgs)
import System.IO (readFile)
import Data.Map.Strict (Map, (!), insert, elems, fromList, toList, findWithDefault, size, empty, member, findMin, findMax, singleton, filter)
import qualified Data.Map.Strict as M
--import qualified Data.Array as A
import Data.List (permutations)
import Data.List.Split (splitOn)
import Data.Complex (Complex((:+)), realPart, imagPart)

import UI.NCurses

type Instructions = Map Integer Integer

data ComputerState = Ready | Running | Blocked | Done
  deriving (Show, Eq)

data Computer = Computer {
  state :: ComputerState,
  memory :: Map Integer Integer,
  iptr :: Integer, -- instruction pointer
  base :: Integer, -- base offset
  input :: [Integer],
  output :: [Integer]}

computer0 = Computer { state = Ready, memory = empty, iptr = 0, base = 0, input = [], output = [] }

run :: Computer -> Computer -- (Integer, Integer) -> Instructions -> [Integer] -> [Integer]
run c
  | state c == Ready = run $ c {state = Running}
  | state c == Done = c
  | state c == Blocked = if input c == [] then c else run $ c {state = Running}
  | otherwise = -- running
    --traceShow (iptr c, base c, memory c) $
    case instr `mod` 100 of
      1 -> -- add
        run $ c {iptr = i+4, memory = insert (addr 3) (arg 1 + arg 2) instructions}
      2 -> -- multiply
        run $ c {iptr = i+4, memory = insert (addr 3) (arg 1 * arg 2) instructions}
      3 -> -- read input
        if (null $ input c)
        then c {state = Blocked}
        else run $ c {iptr = i+2, memory = insert (addr 1) (head $ input c) instructions, input = tail $ input c}
      4 -> -- output
        --traceShow (arg 1) $
        run $ c {iptr = i+2, output = output c ++ [arg 1]}
      5 -> -- jump-if-true
        run $ c {iptr = if arg 1 == 0 then i+3 else arg 2}
      6 -> -- jump-if-false
        run $ c {iptr = if arg 1 == 0 then arg 2 else i+3}
      7 -> -- less than
        run $ c {iptr = i+4, memory = insert (addr 3) (if arg 1 < arg 2 then 1 else 0) instructions}
      8 -> -- equals
        run $ c {iptr = i+4, memory = insert (addr 3) (if arg 1 == arg 2 then 1 else 0) instructions}
      9 -> -- set relative base
        run $ c {iptr = i+2, base = base c + arg 1}
      99 -> -- halt
        run $ c {state = Done}
      _ -> error "unknown opcode"
  where instructions = memory c
        i = iptr c
        instr = instructions!i
        ii x = findWithDefault 0 x instructions
        arg :: Integer -> Integer
        arg n = case (instr `mod` (100*10^n)) `div` (10*10^n) of
          0 -> ii $ ii $ i+n
          1 -> ii $ i+n
          2 -> ii ((ii $ i+n) + base c)
          _ -> error "bad argument mode"
        addr n = case (instr `mod` (100*10^n)) `div` (10*10^n) of
          0 -> ii $ i+n
          1 -> error "address in mode 1"
          2 -> --trace "address in mode 2" $
            ii (i+n) + base c
          _ -> error $ "address in unknown mode"

instance Ord a => Ord (Complex a) where
  compare a b | ar < br || ar == br && ai < bi = LT
              | ar > br || ar == br && ai > bi = GT
              | otherwise = EQ
    where (ar, ai) = toParts a
          (br, bi) = toParts b

toParts :: Complex a -> (a, a)
toParts c = (realPart c, imagPart c)

toIntParts :: RealFrac a => Complex a -> (Int, Int)
toIntParts c = (round $ realPart c, round $ imagPart c)

fromParts :: (a, a) -> Complex a
fromParts (a,b) = a :+ b

data Game = Game {
  m :: Map (Integer, Integer) Integer, -- map
  b :: [(Integer, Integer)], -- ball
  p :: [(Integer, Integer)]} -- paddle
  deriving (Show)

game0 = Game { m = empty, b = [], p = [] }
updateGame :: Game -> [Integer] -> Game
updateGame g [] = g
updateGame g (x:y:t:r) = updateGame (g' { m = insert (x,y) t $ m g}) r
  where g' | t == 4 = g { b = (x,y): b g }
           | t == 3 = g { p = (x,y): p g }
           | otherwise = g
updateGame g _ = error "bad input"

recommendMove :: Game -> Integer
recommendMove g =
  if null (p g) || null (b g) then 0
  else if x == x' && (y == y' || y - 1 == y') then 0
       else case length (b g) of
              0 -> 0
              1 -> signum $ (fst $ head $ b g) - x
              _ -> let [x2, x1] = map fst $ take 2 $ b g
                in signum $ x2 + signum (x2 - x1) - x
  where (x,y) = head $ p g
        (x',y') = head $ b g

display :: [Integer] -> IO ()
display drawCommands = runCurses $ do
    setEcho False
    w <- defaultWindow
    updateWindow w $ do
        draw drawCommands
    render
    waitFor w (\ev -> ev == EventCharacter 'q' || ev == EventCharacter 'Q')
      where draw [] = do
              moveCursor 0 10
              drawString "Done"
            draw (x:y:t:r) =
              if (x,y) == (-1, 0)
              then do
                moveCursor 0 0
                drawString $ show t
                draw r
              else do
                moveCursor (y+5) x
                drawString (case t of
                              0 -> " "
                              1 -> "W"
                              2 -> "B"
                              3 -> "_"
                              4 -> "*")
                draw r

waitFor :: Window -> (Event -> Bool) -> Curses ()
waitFor w p = loop where
    loop = do
        ev <- getEvent w Nothing
        case ev of
            Nothing -> loop
            Just ev' -> if p ev' then return () else loop

play :: Game -> Computer -> IO (Game)
play game c = runCurses $ do
  setEcho False
  w <- defaultWindow
  loop w game $ run c
    where loop :: Window -> Game -> Computer -> Curses(Game)
          loop w g c
            | state c == Done = do
                processDrawCommands w g' $ output c
                updateWindow w $ do
                  moveCursor 0 10
                  drawString "Press q to quit"
                render
                let qloop = do
                      ev <- getEvent w Nothing
                      case ev of
                        Nothing -> qloop
                        Just (EventCharacter 'q') -> return g'
                        Just (EventCharacter 'Q') -> return g'
                        _ -> qloop
                qloop
            | state c == Blocked = do -- waiting for input
                processDrawCommands w g' $ output c
                moveEvent <- getEvent w (Just 20)
                let joystick = case moveEvent of
                                 Just (EventSpecialKey KeyLeftArrow) -> -1
                                 Just (EventSpecialKey KeyRightArrow) -> 1
                                 Nothing -> recommendMove g'
                                 Just (EventCharacter 'j') -> -1
                                 Just (EventCharacter 'J') -> -1
                                 Just (EventCharacter 'k') -> 1
                                 Just (EventCharacter 'K') -> 1
                                 _ -> 0
                updateWindow w $ do
                  moveCursor 0 10
                  drawString $ case joystick of
                    -1 -> "move left"
                    0 -> "no move"
                    1 -> "move right"
                render
                loop w g' $ run $ c {state = Running, output = [], input = [joystick]}
            | otherwise = loop w game $ run c
              where g' = updateGame g $ output c
          processDrawCommands :: Window -> Game -> [Integer] -> Curses()
          processDrawCommands w g [] = render
          processDrawCommands w g (x:y:t:r) = do
            updateWindow w $ do
              if (x,y) == (-1, 0)
              then do
              moveCursor 0 0
              drawString $ show t
              else do
              moveCursor 1 0
              drawString $ show $ take 25 $ b g
              moveCursor 2 0
              drawString $ show $ take 25 $ p g
              moveCursor (y+5) x
              drawString (case t of
                             0 -> " "
                             1 -> "W"
                             2 -> "B"
                             3 -> "_"
                             4 -> "*"
                             _ -> error ("Unknown input " ++ show t)
                         )
            processDrawCommands w g r

main = do
  [instructionFile] <- getArgs
  instructionStrings <- readFile instructionFile
  let instructions = fromList . zip [0 ..] $ map read $ splitOn "," instructionStrings

  -- putStrLn "Part 1"
  -- let c = run $ computer0 {memory = instructions}
  -- -- putStrLn $ show ("state", state c, "drawCommands", output c)
  -- let game = updateGame game0 $ output c
  -- putStrLn $ show $ size $ M.filter (== 2) $ m game
  -- putStrLn "Part 1b"
  -- display $ output c
  putStrLn "Part 2"
  let instructions2 = insert 0 2 instructions -- play for free by setting memory address 0 to 2
  game <- play game0 $ computer0 {memory = instructions2}
  putStrLn $ show $ head $ p game
  -- let moves = moves >>= (play . run (0,0) instructions2)
  -- -- moves <- play $ run (0,0) instructions2 $ repeat 0
  -- moves' <- moves
  -- putStrLn $ show $ length moves'

  -- let moves = moves >>= (play . run (0,0) instructions2)
  -- moves' <- moves
  -- putStrLn $ show moves'
  -- moves' <- play $ run (0,0) instructions2 []
  -- putStrLn $ show ("moves", moves')
  -- let moves = run (0,0) instructions2 $ repeat 0
  -- putStrLn $ show ("moves", moves)
  -- display moves

  -- let out = run (0,0) instructions $ (startColors!origin) : processOutput
  --     (processOutput, finalp, finald, colors) = process (origin,up, startColors) $ out
  -- putStrLn
  --   $ unlines
  --   $ ["done Part 1"
  --     , show $ ("tiles painted: ", size colors)
  --     , show $ ("finalp", finalp, "finald", robot finald)
  --     ]
  -- putStr $ showColors finalp finald colors

  -- putStrLn "Part 2"
  -- let out = run (0,0) instructions $ (startColors2!origin) : processOutput
  --     (processOutput, finalp, finald, colors) = process (origin,up, startColors2) $ out
  -- putStrLn
  --   $ unlines
  --   $ ["done Part 2"
  --     , show $ ("tiles painted: ", size colors)
  --     , show $ ("finalp", finalp, "finald", robot finald)
  --     ]
  -- putStr $ showColors finalp finald colors