text:
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeFamilies #-}
-- |
-- Module     : Numeric.BLAS
-- Copyright  : Copyright (c) 2012 Aleksey Khudyakov <alexey.skladnoy@gmail.com>
-- License    : BSD3
-- Maintainer : Aleksey Khudyakov <alexey.skladnoy@gmail.com>
-- Stability  : experimental
--
-- BLAS operations for immutable vectors and matrices. Current
-- implementation tries hard to evaluate expression built using
-- functions from this module with minimal number of BLAS calls.
module Numeric.BLAS (
    -- * Type class based API
    LinSpace(..)
  , Mul(..)
  , trans
  , conj
    -- * Vector operations
  , dotProduct
  , hermitianProd
  , vectorNorm
  , absSum
  , absIndex
  ) where

import Control.Monad.ST

-- import Data.Complex
import Data.Vector.Generic (Mutable)

import Numeric.BLAS.Bindings (BLAS1,BLAS2,BLAS3,RealType,
                              Trans(..))
import Numeric.BLAS.Expression

-- Vector type classes
import           Data.Vector.Generic         (Vector)
import qualified Data.Vector.Generic         as G
-- Matrix type classes
import           Data.Matrix.Generic           (Transposed(..),Conjugated(..))
import qualified Data.Matrix.Generic         as Mat
-- Concrete vectors
import qualified Data.Vector.Storable         as S
import qualified Data.Vector.Storable.Strided as V
-- Concrete matrices
import           Data.Matrix.Dense     (Matrix)
import           Data.Matrix.Symmetric
  (SymmetricRaw,IsSymmetric,IsHermitian,Conjugate(..),NumberType,IsReal)

import qualified Numeric.BLAS.Mutable as M

import Numeric.BLAS.Mutable (MVectorBLAS)

----------------------------------------------------------------
-- Type class for addition and multiplication
----------------------------------------------------------------

-- | Addition and multiplication by scalar for vectors and matrices.
class LinSpace m a where
  (.+.) :: m a -> m a -> m a
  (.-.) :: m a -> m a -> m a
  ( *.) ::   a -> m a -> m a

instance (LinSpaceM (Mutable m) a, Freeze m a, Num a) => LinSpace m a where
   x .+. y = eval $ Add () (Lit x) (Lit y)
   {-# INLINE (.+.) #-}
   x .-. y = eval $ Sub () (Lit x) (Lit y)
   {-# INLINE (.-.) #-}
   α  *. x = eval $ Scale () α (Lit x)
   {-# INLINE (*.) #-}

-- | Very overloaded operator for matrix and vector multiplication.
class Mul v u where
  -- | Result of multiplication,
  type MulRes v u :: *
  (.*.) :: v -> u -> MulRes v u

-- | Transpose vector or matrix.
trans :: mat a -> Transposed mat a
{-# INLINE trans #-}
trans = Transposed

-- | Conjugate transpose vector or matrix.
conj :: mat a -> Conjugated mat a
{-# INLINE conj #-}
conj  = Conjugated

infixl 6 .+.
infixl 7 .*., *.

----------------------------------------------------------------
-- BLAS 1
----------------------------------------------------------------

-- | Scalar product of vectors
dotProduct :: (BLAS1 a, Vector v a, MVectorBLAS (Mutable v))
           => v a -> v a -> a
{-# INLINE dotProduct #-}
dotProduct v u = runST $ do
  mv <- G.unsafeThaw v
  mu <- G.unsafeThaw u
  M.dotProduct mv mu

-- | Hermitian product of vectors. For real-valued vectors is same
--   as 'dotProduct'.
hermitianProd :: (BLAS1 a, Vector v a, MVectorBLAS (Mutable v))
              => v a -> v a -> a
{-# INLINE hermitianProd #-}
hermitianProd v u = runST $ do
  mv <- G.unsafeThaw v
  mu <- G.unsafeThaw u
  M.hermitianProd mv mu

-- | Euclidean norm of vector
vectorNorm :: (BLAS1 a, Vector v a, MVectorBLAS (Mutable v))
           => v a -> RealType a
{-# INLINE vectorNorm #-}
vectorNorm v
  = runST $ M.vectorNorm =<< G.unsafeThaw v

-- | Sum of absolute values of vector
absSum :: (BLAS1 a, Vector v a, MVectorBLAS (Mutable v))
       => v a -> RealType a
{-# INLINE absSum #-}
absSum v
  = runST $ M.absSum =<< G.unsafeThaw v

-- | Index of element with maximal absolute value
absIndex :: (BLAS1 a, Vector v a, MVectorBLAS (Mutable v))
         => v a -> Int
{-# INLINE absIndex #-}
absIndex v
  = runST $ M.absIndex =<< G.unsafeThaw v

----------------------------------------------------------------
-- Dot product
----------------------------------------------------------------

instance (BLAS1 a, a ~ a') => Mul (Transposed S.Vector a) (S.Vector a') where
  type MulRes (Transposed S.Vector a )
              (           S.Vector a')
             = a
  Transposed v .*. u = dotProduct v u
  {-# INLINE (.*.) #-}
instance (BLAS1 a, a ~ a') => Mul (Transposed V.Vector a) (V.Vector a') where
  type MulRes (Transposed V.Vector a )
              (           V.Vector a')
             = a
  Transposed v .*. u = dotProduct v u
  {-# INLINE (.*.) #-}

----------------------------------------------------------------
--  Vector x Vector => Matrix
----------------------------------------------------------------

instance (BLAS2 a, a ~ a') => Mul (S.Vector a) (Transposed S.Vector a') where
  type MulRes (           S.Vector a )
              (Transposed S.Vector a')
             = Matrix a
  v .*. Transposed u = eval $ VecT () (Lit v) (Lit u)
  {-# INLINE (.*.) #-}
instance (BLAS2 a, a ~ a') => Mul (S.Vector a) (Conjugated S.Vector a') where
  type MulRes (           S.Vector a )
              (Conjugated S.Vector a')
             = Matrix a
  v .*. Conjugated u = eval $ VecH () (Lit v) (Lit u)
  {-# INLINE (.*.) #-}

instance (BLAS2 a, a ~ a') => Mul (V.Vector a) (Transposed V.Vector a') where
  type MulRes (           V.Vector a )
              (Transposed V.Vector a')
             = Matrix a
  v .*. Transposed u = eval $ VecT () (Lit v) (Lit u)
  {-# INLINE (.*.) #-}
instance (BLAS2 a, a ~ a') => Mul (V.Vector a) (Conjugated V.Vector a') where
  type MulRes (           V.Vector a )
              (Conjugated V.Vector a')
             = Matrix a
  v .*. Conjugated u = eval $ VecH () (Lit v) (Lit u)
  {-# INLINE (.*.) #-}

----------------------------------------------------------------
-- Dense matrix x Vector
----------------------------------------------------------------

-- Strided
instance (BLAS2 a, a ~ a') => Mul (Matrix a) (V.Vector a') where
  type MulRes (Matrix   a )
              (V.Vector a')
             = V.Vector a
  m .*. v = eval $ MulMV () (Lit m) (Lit v)
  {-# INLINE (.*.) #-}
instance (BLAS2 a, a ~ a') => Mul (Transposed Matrix a) (V.Vector a') where
  type MulRes (Transposed Matrix a)
              (V.Vector a')
             = V.Vector a
  Transposed m .*. v = eval $ MulTMV () Trans (Lit m) (Lit v)
  {-# INLINE (.*.) #-}
instance (BLAS2 a, a ~ a') => Mul (Conjugated Matrix a) (V.Vector a') where
  type MulRes (Conjugated Matrix a)
              (V.Vector a')
             = V.Vector a
  Conjugated m .*. v = eval $ MulTMV () ConjTrans (Lit m) (Lit v)
  {-# INLINE (.*.) #-}

-- Storable
instance (BLAS2 a, a ~ a') => Mul (Matrix a) (S.Vector a') where
  type MulRes (Matrix   a )
              (S.Vector a')
             = S.Vector a
  m .*. v = eval $ MulMV () (Lit m) (Lit v)
  {-# INLINE (.*.) #-}
instance (BLAS2 a, a ~ a') => Mul (Transposed Matrix a) (S.Vector a') where
  type MulRes (Transposed Matrix a)
              (S.Vector a')
             = S.Vector a
  Transposed m .*. v = eval $ MulTMV () Trans (Lit m) (Lit v)
  {-# INLINE (.*.) #-}
instance (BLAS2 a, a ~ a') => Mul (Conjugated Matrix a) (S.Vector a') where
  type MulRes (Conjugated Matrix a)
              (S.Vector a')
             = S.Vector a
  Conjugated m .*. v = eval $ MulTMV () ConjTrans (Lit m) (Lit v)
  {-# INLINE (.*.) #-}

----------------------------------------------------------------
-- Symmetric matrix x Vector
----------------------------------------------------------------

instance (BLAS2 a, Conjugate a, a ~ a') => Mul (SymmetricRaw IsHermitian a) (S.Vector a') where
  type MulRes (SymmetricRaw IsHermitian a)
              (S.Vector a')
             = S.Vector a
  m .*. v = eval $ MulMV () (Lit m) (Lit v)
  {-# INLINE (.*.) #-}

instance (BLAS2 a, NumberType a ~ IsReal, a ~ a') => Mul (SymmetricRaw IsSymmetric a) (S.Vector a') where
  type MulRes (SymmetricRaw IsSymmetric a)
              (S.Vector a')
             = S.Vector a
  m .*. v = eval $ MulMV () (Lit m) (Lit v)
  {-# INLINE (.*.) #-}

----------------------------------------------------------------
-- Matrix x Matrix for dense matrices
----------------------------------------------------------------

instance (BLAS3 a, a ~ a') => Mul (Matrix a) (Matrix a') where
  type MulRes (Matrix a )
              (Matrix a')
             = Matrix a
  m .*. n = eval $ MulMM () NoTrans (Lit m) NoTrans (Lit n)
  {-# INLINE (.*.) #-}

instance (BLAS3 a, a ~ a') => Mul (Matrix a) (Transposed Matrix a') where
  type MulRes (           Matrix a )
              (Transposed Matrix a')
             = Matrix a
  m .*. Transposed n = eval $ MulMM () NoTrans (Lit m) Trans (Lit n)
  {-# INLINE (.*.) #-}

instance (BLAS3 a, a ~ a') => Mul (Matrix a) (Conjugated Matrix a') where
  type MulRes (           Matrix a )
              (Conjugated Matrix a')
             = Matrix a
  m .*. Conjugated n = eval $ MulMM () NoTrans (Lit m) ConjTrans (Lit n)
  {-# INLINE (.*.) #-}

instance (BLAS3 a, a ~ a') => Mul (Transposed Matrix a) (Matrix a') where
  type MulRes (Transposed Matrix a )
              (           Matrix a')
             = Matrix a
  Transposed m .*. n = eval $ MulMM () Trans (Lit m) NoTrans (Lit n)
  {-# INLINE (.*.) #-}

instance (BLAS3 a, a ~ a') => Mul (Transposed Matrix a) (Transposed Matrix a') where
  type MulRes (Transposed Matrix a )
              (Transposed Matrix a')
             = Matrix a
  Transposed m .*. Transposed n = eval $ MulMM () Trans (Lit m) Trans (Lit n)
  {-# INLINE (.*.) #-}

instance (BLAS3 a, a ~ a') => Mul (Transposed Matrix a) (Conjugated Matrix a') where
  type MulRes (Transposed Matrix a )
              (Conjugated Matrix a')
             = Matrix a
  Transposed m .*. Conjugated n = eval $ MulMM () Trans (Lit m) ConjTrans (Lit n)
  {-# INLINE (.*.) #-}

instance (BLAS3 a, a ~ a') => Mul (Conjugated Matrix a) (Matrix a') where
  type MulRes (Conjugated Matrix a )
              (           Matrix a')
             = Matrix a
  Conjugated m .*. n = eval $ MulMM () ConjTrans (Lit m) NoTrans (Lit n)
  {-# INLINE (.*.) #-}

instance (BLAS3 a, a ~ a') => Mul (Conjugated Matrix a) (Transposed Matrix a') where
  type MulRes (Conjugated Matrix a )
              (Transposed Matrix a')
             = Matrix a
  Conjugated m .*. Transposed n = eval $ MulMM () ConjTrans (Lit m) Trans (Lit n)
  {-# INLINE (.*.) #-}

instance (BLAS3 a, a ~ a') => Mul (Conjugated Matrix a) (Conjugated Matrix a') where
  type MulRes (Conjugated Matrix a )
              (Conjugated Matrix a')
             = Matrix a
  Conjugated m .*. Conjugated n = eval $ MulMM () ConjTrans (Lit m) ConjTrans (Lit n)
  {-# INLINE (.*.) #-}

----------------------------------------------------------------
-- Symmetric matrix x Dense matrix
----------------------------------------------------------------

instance (BLAS3 a, a ~ a') => Mul (Matrix a) (SymmetricRaw IsSymmetric a') where
  type MulRes (Matrix a)
              (SymmetricRaw IsSymmetric a')
             = Matrix a
  m .*. sym = eval $ MulSymMM () M.RightSide (Lit sym) (Lit m)
  {-# INLINE (.*.) #-}

instance (BLAS3 a, a ~ a') => Mul (SymmetricRaw IsSymmetric a') (Matrix a) where
  type MulRes (SymmetricRaw IsSymmetric a')
              (Matrix a)
             = Matrix a
  sym .*. m = eval $ MulSymMM () M.LeftSide (Lit sym) (Lit m)
  {-# INLINE (.*.) #-}

----------------------------------------------------------------
-- Symmetric matrix x Dense matrix
----------------------------------------------------------------

instance (BLAS3 a, Conjugate a, a ~ a') => Mul (Matrix a) (SymmetricRaw IsHermitian a') where
  type MulRes (Matrix a)
              (SymmetricRaw IsHermitian a')
             = Matrix a
  m .*. sym = eval $ MulHerMM () M.RightSide (Lit sym) (Lit m)
  {-# INLINE (.*.) #-}

instance (BLAS3 a, Conjugate a, a ~ a') => Mul (SymmetricRaw IsHermitian a') (Matrix a) where
  type MulRes (SymmetricRaw IsHermitian a')
              (Matrix a)
             = Matrix a
  sym .*. m = eval $ MulHerMM () M.LeftSide (Lit sym) (Lit m)
  {-# INLINE (.*.) #-}

max_stars_repo_name: Shimuuar/blas-lapack

max_issues_repo_name: Shimuuar/blas-lapack

max_forks_repo_name: Shimuuar/blas-lapack


text:
{-# LANGUAGE FlexibleContexts #-}

module Matrix.State
  (
    replaceMatrixElement
  , replaceVector
  , replaceElement
  , getElement
  , applyFMatrixElement
  , MatrixIndex
  , VectorIndex
  , MatrixMutation
  , VectorMutation
  , module Numeric.LinearAlgebra.HMatrix
  ) where

import Numeric.LinearAlgebra.HMatrix hiding (corr)

type MatrixIndex      = (MatrixRowIndex, VectorIndex)
type VectorIndex      = Int
type MatrixRowIndex   = Int
type MatrixMutation a = Matrix a -> Matrix a
type VectorMutation a = Vector a -> Vector a

getElement :: (Element a, Indexable (Vector a) a) => Matrix a -> MatrixIndex -> a
getElement m (i,j) = m ! i ! j

applyFMatrixElement :: (Element a, Indexable (Vector a) a)
  => (a -> a) -> MatrixIndex -> MatrixMutation a
applyFMatrixElement f i m = replaceMatrixElement i (f (getElement m i)) m

replaceMatrixElement :: (Element a) => MatrixIndex -> a -> MatrixMutation a
replaceMatrixElement (i,j) e m = replaceVector i v' m
  where
    v  = m ! i
    v' = replaceElement j e v

replaceVector :: (Element a) => MatrixRowIndex -> Vector a -> MatrixMutation a
replaceVector i v m = if not $ inBoundsMatrixRow i m
                      then m
                      else fromRows $ rowsAbove ++ v : rowsBelow
  where
    rowsAbove = toRows $ takeRows i       m
    rowsBelow = toRows $ dropRows (i + 1) m

replaceElement :: (Element a) => VectorIndex -> a -> VectorMutation a
replaceElement i e v = if not $ inBounds i v
                       then v
                       else fromList $ elemsBefore ++ e : elemsAfter
  where
    elemsBefore = take  i      $ lV
    elemsAfter  = drop (i + 1) $ lV
    lV         = toList v

inBounds :: (Element a) => VectorIndex -> Vector a -> Bool
inBounds i v = (i >= 0) && (i < (length $ toList v))

inBoundsMatrixRow :: MatrixRowIndex -> Matrix a -> Bool
inBoundsMatrixRow i m = (i >= 0) && (i < (rows m))

ext: hs

max_stars_repo_name: eklinkhammer/neural-algorithms

max_issues_repo_name: eklinkhammer/neural-algorithms

max_forks_repo_name: eklinkhammer/neural-algorithms


text:
import Data.Complex (Complex, cis)

rootsOfUnity :: (Enum a, Floating a) => a -> [Complex a]
rootsOfUnity n =
  [ cis (2 * pi * k / n)
  | k <- [0 .. n - 1] ]

main :: IO ()
main = mapM_ print $ rootsOfUnity 3

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module MyPhoto.Actions.Montage
    ( montage
    ) where

import           Control.Concurrent ( getNumCapabilities )
import           Control.Concurrent.Async ( mapConcurrently )
import           Control.Concurrent.MSem as MS
import           Control.Monad
import           GHC.Conc ( numCapabilities )
import           System.Console.GetOpt
import           System.IO.Temp
import           System.Directory
import           System.Exit
import           System.FilePath
import           System.Process
import           Graphics.Netpbm
import qualified Data.ByteString as B
import           Data.Maybe (maybe)
import qualified Statistics.Sample as S
import qualified Data.Vector as V (fromList)

import MyPhoto.Model
import MyPhoto.Utils

help :: String
help = let
    header = "Usage: montage OUTPUT_FILE"
  in usageInfo header options

data Options
  = Options
  { optOutputFile :: Maybe FilePath
  , optHelp       :: Bool
  } deriving Show

defaultOptions :: Options
defaultOptions
  = Options
  { optOutputFile = Nothing
  , optHelp       = False
  }

options :: [OptDescr (Options -> Options)]
options =
  [ Option ['o'] ["--output-file"]
      (OptArg (\ fp opts -> opts { optOutputFile = fp })
        "OUTPUT_FILE")
      "file to wriite the montage to (expected to be e.g. a PNG)"
  , Option ['h'] ["help"]
      (NoArg (\ opts -> opts { optHelp = True }))
      "print help"
  ]

getMyOpts :: [String] -> IO (Options, [String])
getMyOpts argv = case getOpt Permute options argv of
                   (o, n, []  ) -> return (foldl (flip id) defaultOptions o, n)
                   (_, _, errs) -> ioError (userError (concat errs ++ help))

montageImpl:: [String] -> [Img] -> PActionBody
montageImpl _    []            = return (Right [])
montageImpl args imgs@(img1:_) = do
  (opts, _) <- getMyOpts args

  if optHelp opts
  then return (Left help)
  else let
    (bn,ext) = splitExtensions img1
    outputFile = case optOutputFile opts of
          Just outputFile' -> outputFile'
          Nothing -> bn ++ "_" ++ show (length imgs) ++ "_MONTAGE.png"
    in do
      (_, _, _, pHandle) <- createProcess (proc "montage" (["-geometry", "100x100+2+2"] ++ imgs ++ [outputFile]))
      exitCode <- waitForProcess pHandle
      unless (exitCode == ExitSuccess) $
        fail ("Resize failed with " ++ show exitCode)
      return (Right imgs)

montage :: PrePAction
montage args = logSeparator "montage" <> PAction (montageImpl args)

max_stars_repo_name: maxhbr/myphoto

max_issues_repo_name: maxhbr/myphoto

max_forks_repo_name: maxhbr/myphoto


text:
{- HLINT ignore "Evaluate" -}

{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE DuplicateRecordFields #-}
{-# LANGUAGE NamedFieldPuns #-}
{-# LANGUAGE NumericUnderscores #-}
{-# LANGUAGE OverloadedLists #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}

-- |
-- Copyright: © 2021 IOHK
-- License: Apache-2.0
--
-- Provides functions for selecting coins for use as collateral from a UTxO
-- set.
--
-- See the documentation for 'performSelection' for more details.
--
module Cardano.Wallet.CoinSelection.Internal.Collateral
    (
    -- * Public API

      performSelection
    , PerformSelection
    , PerformSelectionOf
    , SelectionConstraints (..)
    , SelectionParams
    , SelectionParamsOf (..)
    , SelectionResult
    , SelectionResultOf (..)
    , selectionResultEmpty
    , SelectionCollateralError
    , SelectionCollateralErrorOf (..)
    , SearchSpaceLimit (..)
    , searchSpaceLimitDefault

    -- * Internal API

    -- ** Selecting collateral by giving priority to smallest values first
    , selectCollateralSmallest

    -- ** Selecting collateral by giving priority to largest values first
    , selectCollateralLargest

    -- ** Guarding search space size
    , SearchSpaceRequirement (..)
    , guardSearchSpaceSize

    -- ** Generating submaps
    , submaps

    -- ** Generating subsequences
    , subsequencesOfSize
    , numberOfSubsequencesOfSize

    -- ** Control flow
    , firstRight
    , takeUntil
    )
    where

import Cardano.Wallet.Primitive.Types.Coin
    ( Coin )
import Cardano.Wallet.Primitive.Types.Tx
    ( TxIn )
import Data.Function
    ( (&) )
import Data.IntCast
    ( intCast, intCastMaybe )
import Data.List.NonEmpty
    ( NonEmpty )
import Data.Map.Strict
    ( Map )
import Data.Maybe
    ( listToMaybe, mapMaybe )
import Data.Ord
    ( Down (..) )
import Data.Semigroup
    ( sconcat )
import Data.Set
    ( Set )
import GHC.Generics
    ( Generic )

import Prelude

import qualified Data.Foldable as F
import qualified Data.List as L
import qualified Data.Map.Strict as Map
import qualified Data.Set as Set
import qualified Math.Combinatorics.Exact.Binomial as MathExact
import qualified Numeric.SpecFunctions as MathFast

--------------------------------------------------------------------------------
-- Public API
--------------------------------------------------------------------------------

-- | The type of all functions that perform selections.
--
type PerformSelectionOf inputId =
    SelectionConstraints ->
    SelectionParamsOf inputId ->
    Either (SelectionCollateralErrorOf inputId) (SelectionResultOf inputId)

-- | The default type for 'PerformSelectionOf'.
--
type PerformSelection = PerformSelectionOf TxIn

-- | Specifies all constraints required for collateral selection.
--
-- Selection constraints:
--
--    - are dependent on the current set of protocol parameters.
--
--    - are not specific to a given selection.
--
--    - place limits on the selection algorithm, enabling it to produce
--      selections that are acceptable to the ledger.
--
data SelectionConstraints = SelectionConstraints
    { maximumSelectionSize
        :: Int
        -- ^ An upper bound on the number of unique coins that can be selected
        -- as collateral.
    , searchSpaceLimit
        :: SearchSpaceLimit
        -- ^ An upper bound on the search space size, to protect the wallet
        -- against computations that use excessive amounts of time or space.
    }
    deriving (Eq, Generic, Show)

-- | Specifies all parameters that are specific to a given selection.
--
data SelectionParamsOf inputId = SelectionParams
    { coinsAvailable
        :: Map inputId Coin
        -- ^ The set of all coins available for selection as collateral.
    , minimumSelectionAmount
        :: Coin
        -- ^ A lower bound on the sum of coins to be selected as collateral.
    }
    deriving (Eq, Generic, Show)

-- | The default type for 'SelectionParamsOf'.
--
type SelectionParams = SelectionParamsOf TxIn

-- | Specifies an upper bound on the search space size.
--
data SearchSpaceLimit
    = SearchSpaceLimit Int
    -- ^ Specifies an upper bound on the number of coin combinations that can
    -- be considered in any single step.
    | UnsafeNoSearchSpaceLimit
    -- ^ Specifies that there is no search space limit. This should only be
    -- used for testing purposes.
    deriving (Eq, Show)

-- | The default search space limit.
--
-- This constant is used by the test suite, so we can be reasonably confident
-- that performing selections with this limit will not use inordinate amounts
-- of time and space.
--
searchSpaceLimitDefault :: SearchSpaceLimit
searchSpaceLimitDefault = SearchSpaceLimit 1_000_000

-- | Represents a successful selection of collateral.
--
newtype SelectionResultOf inputId = SelectionResult
    { coinsSelected :: Map inputId Coin
        -- ^ The coins that were selected for collateral.
    }
    deriving (Eq, Generic, Show)

-- | The default type for 'SelectionResultOf'.
--
type SelectionResult = SelectionResultOf TxIn

-- | A completely empty result, with no inputs selected.
--
selectionResultEmpty :: SelectionResultOf inputId
selectionResultEmpty = SelectionResult
    { coinsSelected = Map.empty
    }

-- | Represents an unsuccessful attempt to select collateral.
--
data SelectionCollateralErrorOf inputId = SelectionCollateralError
    { largestCombinationAvailable :: Map inputId Coin
        -- ^ The largest combination of coins available.
    , minimumSelectionAmount :: Coin
        -- ^ A lower bound on the sum of coins to be selected as collateral.
    }
    deriving (Eq, Generic, Show)

-- | The default type for `SelectionCollateralErrorOf`.
--
type SelectionCollateralError = SelectionCollateralErrorOf TxIn

-- | Selects coins for collateral.
--
-- This function tries two strategies in the following order, picking the first
-- strategy that succeeds:
--
--    1. Attempt to select an amount of collateral that is as small as possible.
--    2. Attempt to select collateral from the largest coins available.
--
-- The first strategy, given unlimited computation time, will always produce an
-- optimal result: the smallest possible amount of collateral. However, if the
-- required search space is large, and if the 'searchSpaceLimit' parameter is
-- set to a value that's smaller than the required search space size, then this
-- strategy will fail without computing a result.
--
-- The second strategy sacrifices optimality and always produces a result if
-- one is available, by looking only at the very largest coins available. This
-- result can be computed very quickly, without using much search space.
--
-- The combination of these two strategies means that we can satisfy the
-- following properties:
--
-- If the attempt to select collateral succeeds:
--
--    >>> sum  coinsSelected ≥ minimumSelectionAmount
--    >>> size coinsSelected ≤ maximumSelectionSize
--    >>>      coinsSelected ⊆ coinsAvailable
--
-- If the attempt to select collateral fails:
--
--    >>> sum  largestCombinationAvailable < minimumSelectionAmount
--    >>> size largestCombinationAvailable ≤ maximumSelectionSize
--    >>>      largestCombinationAvailable ⊆ coinsAvailable
--
performSelection :: forall inputId. Ord inputId => PerformSelectionOf inputId
performSelection constraints =
    firstRight $ fmap ($ constraints)
        [ selectCollateralSmallest
        , selectCollateralLargest
        ]

--------------------------------------------------------------------------------
-- Internal API
--------------------------------------------------------------------------------

--------------------------------------------------------------------------------
-- Selecting collateral by giving priority to smallest values first
--------------------------------------------------------------------------------

-- | Attempts to select an amount of collateral that is as small as possible.
--
-- This function, given unlimited computation time, will always produce an
-- optimal result: the smallest possible amount of collateral. However, if the
-- required search space is large, and if the 'searchSpaceLimit' parameter is
-- set to a value that's smaller than the required search space size, then this
-- function will return without computing a result.
--
selectCollateralSmallest
    :: forall inputId. Ord inputId => PerformSelectionOf inputId
selectCollateralSmallest constraints params =
    case smallestValidCombination of
        Just coinsSelected ->
            Right SelectionResult {coinsSelected}
        Nothing ->
            Left SelectionCollateralError
                { largestCombinationAvailable = mempty
                , minimumSelectionAmount
                }
  where
    coinsToConsider :: [(inputId, Coin)]
    coinsToConsider = coinsAvailable
        & Map.toList
        & L.sortOn snd
        & takeUntil ((>= minimumSelectionAmount) . snd)

    numberOfCoinsToConsider :: Int
    numberOfCoinsToConsider = length coinsToConsider

    smallestValidCombination :: Maybe (Map inputId Coin)
    smallestValidCombination =
        listToMaybe $ L.sortOn F.fold validCombinations
      where
        validCombinations :: [Map inputId Coin]
        validCombinations =
            mapMaybe smallestValidCombinationOfSize
            [1 .. maximumSelectionSize]

    smallestValidCombinationOfSize :: Int -> Maybe (Map inputId Coin)
    smallestValidCombinationOfSize size =
        guardSearchSpaceSize searchSpaceRequirement searchSpaceLimit result
      where
        result :: Maybe (Map inputId Coin)
        result = coinsToConsider
            & (`subsequencesOfSize` size)
            & fmap (\ics -> (ics, F.foldMap snd ics))
            & L.sortOn snd
            & L.dropWhile ((< minimumSelectionAmount) . snd)
            & listToMaybe
            & fmap (Map.fromList . fst)

        searchSpaceRequirement :: SearchSpaceRequirement
        searchSpaceRequirement = maybe
            SearchSpaceRequirementUnknown
            SearchSpaceRequirement
            (numberOfCoinsToConsider `numberOfSubsequencesOfSize` size)

    SelectionConstraints
        { maximumSelectionSize
        , searchSpaceLimit
        } = constraints
    SelectionParams
        { coinsAvailable
        , minimumSelectionAmount
        } = params

--------------------------------------------------------------------------------
-- Selecting collateral by giving priority to largest values first
--------------------------------------------------------------------------------

-- | Selects collateral from the largest coins available.
--
-- This function sacrifices optimality and always produces a result if one is
-- available, by looking only at the very largest coins available.
--
-- This result can be computed very quickly, without using much search space.
--
selectCollateralLargest
    :: forall inputId. Ord inputId => PerformSelectionOf inputId
selectCollateralLargest constraints params =
    case smallestValidSubmapOfLargestCombinationAvailable of
        Just coinsSelected ->
            Right SelectionResult {coinsSelected}
        Nothing ->
            Left SelectionCollateralError
                { largestCombinationAvailable
                , minimumSelectionAmount
                }
  where
    largestCombinationAvailable :: Map inputId Coin
    largestCombinationAvailable =
        coinsAvailable
            & Map.toList
            & L.sortOn (Down . snd)
            & L.take maximumSelectionSize
            & Map.fromList

    smallestValidSubmapOfLargestCombinationAvailable :: Maybe (Map inputId Coin)
    smallestValidSubmapOfLargestCombinationAvailable =
        largestCombinationAvailable
            & submaps
            & Set.toList
            & fmap (\ics -> (ics, F.fold ics))
            & L.sortOn snd
            & L.dropWhile ((< minimumSelectionAmount) . snd)
            & fmap fst
            & listToMaybe

    SelectionConstraints
        { maximumSelectionSize
        } = constraints
    SelectionParams
        { coinsAvailable
        , minimumSelectionAmount
        } = params

--------------------------------------------------------------------------------
-- Generating submaps
--------------------------------------------------------------------------------

-- | Generates all submaps of a given map.
--
-- This function is analogous to 'Set.powerSet'.
--
-- For a map 'm' of size 'n', this function will generate all possible submaps,
-- including the empty map and the original map 'm'.
--
submaps :: forall a b. (Ord a, Ord b) => Map a b -> Set (Map a b)
submaps m = Set.map (Map.restrictKeys m) (Set.powerSet (Map.keysSet m))

--------------------------------------------------------------------------------
-- Guarding search space size
--------------------------------------------------------------------------------

data SearchSpaceRequirement
    = SearchSpaceRequirement Int
      -- ^ Indicates a known search space requirement.
    | SearchSpaceRequirementUnknown
      -- ^ Indicates that the search space requirement is unknown.

guardSearchSpaceSize
    :: SearchSpaceRequirement
    -- ^ The search space requirement
    -> SearchSpaceLimit
    -- ^ The search space limit
    -> Maybe a
    -- ^ A computation that potentially yields a value
    -> Maybe a
    -- ^ The guarded computation
guardSearchSpaceSize requirement limit =
    case requirement of
        -- When the search space requirement is unknown, err on the side of
        -- caution and avoid evaluating the computation, unless the caller
        -- has explicitly specified that there is no limit:
        SearchSpaceRequirementUnknown ->
            case limit of
                SearchSpaceLimit _       -> const Nothing
                UnsafeNoSearchSpaceLimit -> id
        -- When the search space requirement is known, only evaluate the
        -- computation if the requirement is not greater than the limit:
        SearchSpaceRequirement r ->
            case limit of
                SearchSpaceLimit l | l < r -> const Nothing
                SearchSpaceLimit _         -> id
                UnsafeNoSearchSpaceLimit   -> id

--------------------------------------------------------------------------------
-- Generating subsequences
--------------------------------------------------------------------------------

-- | Computes the number of subsequences generated by 'subsequencesOfSize'.
--
-- This function can be used to determine whether calling 'subsequencesOfSize'
-- would use an excessive amount of time and space, and if so, avoid calling
-- it.
--
-- Returns 'Nothing' if the result is larger than 'maxBound :: Int'.
--
numberOfSubsequencesOfSize
    :: Int
    -- ^ Indicates the size of the sequence.
    -> Int
    -- ^ Indicates the size of subsequences.
    -> Maybe Int
numberOfSubsequencesOfSize n k
    | k <  0 || n <   0      = Nothing
    | k == 0 || k ==  n      = Just 1
    | k == 1 || k == (n - 1) = Just n
    | resultOutOfBounds      = Nothing
    | otherwise              = intCastMaybe resultExact
  where
    resultExact :: Integer
    resultExact = MathExact.choose
        (intCast @Int @Integer n)
        (intCast @Int @Integer k)

    resultFast :: Integer
    resultFast = floor (MathFast.choose n k)

    resultOutOfBounds :: Bool
    resultOutOfBounds = False
        || resultFast  < 0
        || resultFast  > intCast @Int @Integer (maxBound @Int)
        || resultExact < 0
        || resultExact > intCast @Int @Integer (maxBound @Int)

-- | Generates all subsequences of size 'k' from a particular sequence.
--
-- Warning: this function can use an excessive amount of time and space.
--
-- To check how many results would be returned without actually generating
-- them, use the 'numberOfSubsequencesOfSize' function.
--
-- Properties:
--
--    >>> all (== k) (length <$> xs `subsequencesOfSize` k)
--
--    >>> length (xs `subsequencesOfSize` k) ==
--    >>>     length xs `numberOfSubsequencesOfSize` k
--
subsequencesOfSize
    :: [a]
    -- ^ The sequence from which to generate subsequences.
    -> Int
    -- ^ The size 'k' of subsequences to generate.
    -> [[a]]
    -- ^ All subsequences of size 'k'.
subsequencesOfSize xs k
    | k <= 0 =
        []
    | k > length xs =
        []
    | otherwise =
        case drop k (subsequencesBySize xs) of
            result : _ -> result
            [] -> []
  where
    subsequencesBySize [] = [[[]]]
    subsequencesBySize (y : ys) =
        zipWith
            (++)
            ([] : map (map (y :)) next)
            (next ++ [[]])
      where
        next = subsequencesBySize ys

--------------------------------------------------------------------------------
-- Control flow
--------------------------------------------------------------------------------

-- | Applies a sequence of functions to an argument until one succeeds.
--
-- This function iterates through the specified sequence from left to right,
-- applying each function to the given argument, and returning the very first
-- 'Right' result encountered, without evaluating the subsequent functions.
--
-- If none of the given functions produces a 'Right' result, then this function
-- returns the 'Left' result produced by the last function in the sequence.
--
firstRight :: NonEmpty (a -> Either e r) -> (a -> Either e r)
firstRight = sconcat

-- | Takes items from a list until a predicate becomes true.
--
-- The returned list is a prefix of the original list, and includes the very
-- first item that satisfies the predicate.
--
takeUntil :: (a -> Bool) -> [a] -> [a]
takeUntil p = foldr (\x ys -> x : if p x then [] else ys) []

max_stars_repo_name: runeksvendsen/cardano-wallet

max_issues_repo_name: cryptorites-scotti/cardano-wallet

max_forks_repo_name: cryptorites-scotti/cardano-wallet


text:
{-# LANGUAGE OverloadedLists #-}
module Tasks.Curve where
import qualified Data.Vector.Storable as VS
import           Numeric.LinearAlgebra (Matrix, Vector, flatten,(><))
import           Data.Monoid ((<>))

hermiteMatrix :: Matrix Float
hermiteMatrix = (4><4)
  [ 2,-2, 1, 1
  ,-3, 3,-2,-1
  , 0, 0, 1, 0
  , 1, 0, 0, 0 ]

curvePoints :: Matrix Float
curvePoints = (4><3)
  [ 20, 0, -10
  , 20, 0, -10
  , 100, 100, 0
  , 100, -100, 0]

countPosition :: Float -> Vector Float
countPosition t = flatten $ (1><4) [t^3,t^2,t^1,1] <> hermMat
  where
    hermMat = hermiteMatrix <> curvePoints

lowerPoints :: [Vector Float]
lowerPoints = map countPosition  [0,0.05..0.999]

curveVertices :: Vector Float
curveVertices = VS.concat . map (VS.++ vData) $ quad =<< [0..length lowerPoints - 2]
  where
    lowerTri i = [lVert $ i+1,lVert i,uVert i] :: [Vector Float]
    upperTri i = [lVert $ i+1, uVert i, uVert $ i+1] :: [Vector Float]
    quad i = lowerTri i ++ upperTri i
    lVert i = lowerPoints!!i
    uVert i = (lowerPoints!!i) + ([0,1,0] :: Vector Float)
    vData = [0,0,1,0,0,1,0,1,0,0,0,0,0,0] :: Vector Float

max_stars_repo_name: Antystenes/CPG

max_issues_repo_name: Antystenes/CPG

max_forks_repo_name: Antystenes/CPG


text:
{-# LANGUAGE ConstraintKinds        #-}
{-# LANGUAGE FlexibleContexts       #-}
{-# LANGUAGE FlexibleInstances      #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE GADTs                  #-}
{-# LANGUAGE InstanceSigs           #-}
{-# LANGUAGE MultiParamTypeClasses  #-}
{-# LANGUAGE QuantifiedConstraints  #-}
{-# LANGUAGE RankNTypes             #-}
{-# LANGUAGE TemplateHaskell        #-}
{-# LANGUAGE TypeOperators          #-}
{-# LANGUAGE UndecidableInstances   #-}

module Q.Stats.Arima where
import           Control.Monad.State
import           Data.Foldable
import           Data.Functor.Identity
import           Data.Random
import           Data.Random.Source
import           Data.RVar
import           Data.Time
import           Numeric.LinearAlgebra
import           Q.Stats.TimeSeries
import           System.Random.Mersenne.Pure64
import           Data.Random.Distribution
import           Data.Random.Distribution.Poisson
import           Data.Random.Distribution.T
import           Data.RVar
import           Statistics.Sample

data Ewma d = Ewma Double d

--ll :: (Ewma d) -> [DataPoint Double] -> (Double -> Double)
ll (Ewma lambda d) datapoints = mapM ll_ datapoints where
  ll_ :: DataPoint LocalTime Double -> State Double Double
  ll_ x@(DataPoint _ v) = do
    vart <- get
    let vart2 = lambda * vart + (1 - lambda) * v * v
    put vart2
    return $ logPdf d (sqrt (v  * v / vart))

--forecast :: (Distribution d Double) => (Ewma d) -> Int ->
forecast :: forall d. (Distribution d Double) => Ewma (d Double) -> StateT Double RVar Double
forecast (Ewma lambda d) = do
  y <- lift $ rvar d
  vart <- get
  let vart2 = lambda * vart + (1 - lambda) * y * y
  put vart2
  return (y * sqrt vart)

--forecastN :: Distribution d Double => Ewma (d Double) -> Int -> Double -> RVar ([Double], Double)
forecastN ewma var0 n =  sample $ runStateT (replicateM n (forecast ewma)) var0

max_stars_repo_name: ghais/lowq

max_issues_repo_name: ghais/lowq

max_forks_repo_name: ghais/lowq


text:
module Main where

import Control.Monad

import Criterion.Main

import Options.Applicative

import Statistics.Quantile.Bench
import Statistics.Quantile.Exact
import Statistics.Quantile.Approximate.Sampling
import Statistics.Quantile.Bench.Gen
import Statistics.Quantile.Types
import Statistics.Quantile.Util

import System.Environment

import System.IO

data Params = Params Bool
  deriving (Eq, Show)

main :: IO ()
main = do
  hSetBuffering stdin LineBuffering
  hSetBuffering stdout LineBuffering
  (fp, fh) <- openTempFile tmpDir "bench"
  void $ countOut fh stream
  hClose fh
  s <- lookupEnv "BENCH_SIZE"
  defaultMain $
    [ benchSelector [benchSize s] "median/exact" median external fp
    , benchBufferingSelector sqrtN [benchSize s] "median/sampling-sqrt(n)" median sampling fp
    , benchBufferingSelector sqrtN [benchSize s] "median/jackknife-sampling-sqrt(n)" median samplingJackknife fp
    ]
  where stream = streamHandle stdin

        benchSize = maybe 10000 read
--        stream = Stream $ gauss 3.14159 10000

max_stars_repo_name: fractalcat/slides

max_issues_repo_name: fractalcat/slides

max_forks_repo_name: fractalcat/slides


text:
{-# LANGUAGE DeriveFunctor                            #-}
{-# LANGUAGE DeriveGeneric                            #-}
{-# LANGUAGE FlexibleContexts                         #-}
{-# LANGUAGE GADTs                                    #-}
{-# LANGUAGE KindSignatures                           #-}
{-# LANGUAGE LambdaCase                               #-}
{-# LANGUAGE NoStarIsType                             #-}
{-# LANGUAGE RankNTypes                               #-}
{-# LANGUAGE RecordWildCards                          #-}
{-# LANGUAGE ScopedTypeVariables                      #-}
{-# LANGUAGE StandaloneDeriving                       #-}
{-# LANGUAGE TypeApplications                         #-}
{-# LANGUAGE TypeFamilyDependencies                   #-}
{-# LANGUAGE TypeInType                               #-}
{-# LANGUAGE TypeOperators                            #-}
{-# LANGUAGE UndecidableInstances                     #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}

module Numeric.Wavelet.Discrete (
    DWavelet(..)
  , DWT(..)
  , cdwt
  , dwt
  , idwt
  , haar
    -- DWD(..)
  -- , dwdApprox, dwdDetail
  -- , flattenDWD
  -- , denseDWD
  -- , haar
  -- , unHaar
  ) where

import           Data.Complex
import           Data.Finite
import           Data.Finite.Internal
import           Data.Kind
import           Data.Proxy
import           Data.Vector.Generic.Sized    (Vector)
import           GHC.Generics
import           GHC.TypeLits.Compare
import           GHC.TypeNats
import           Numeric.Wavelet.Internal.FFT
import qualified Data.Vector                  as UV
import qualified Data.Vector.Generic          as UVG
import qualified Data.Vector.Generic.Sized    as VG
import qualified Data.Vector.Sized            as V

data DWavelet v l a = DW
    { dwHiDecomp :: VG.Vector v l a      -- ^ High-pass decomposition filter
    , dwLoDecomp :: VG.Vector v l a      -- ^ Low-pass decomposition filter
    , dwHiRecon  :: VG.Vector v l a      -- ^ High-pass reconstruction filter
    , dwLoRecon  :: VG.Vector v l a      -- ^ Low-pass reconstruction filter
    }
  deriving (Show, Eq, Ord, Functor, Generic)

data DWT v n a = DWT
    { dwtApprox :: VG.Vector v n a
    , dwtDetail :: VG.Vector v n a
    }
  deriving (Show, Eq, Ord, Functor, Generic)

cdwt
    :: (UVG.Vector v (Complex a), KnownNat l, KnownNat n, 1 <= l, 1 <= n, FFTWReal a)
    => DWavelet  v l (Complex a)
    -> VG.Vector v n (Complex a)
    -> DWT v ((n + l - 1) `Div` 2) (Complex a)
cdwt DW{..} x = DWT{..}
  where
    dwtApprox = downsamp $ convolve x dwLoDecomp
    dwtDetail = downsamp $ convolve x dwHiDecomp

dwt
    :: (UVG.Vector v (Complex a), UVG.Vector v a, KnownNat l, KnownNat n, 1 <= l, 1 <= n, FFTWReal a)
    => DWavelet  v l a
    -> VG.Vector v n a
    -> DWT v ((n + l - 1) `Div` 2) a
dwt DW{..} x = DWT{..}
  where
    dwtApprox = downsamp $ rconvolve x dwLoDecomp
    dwtDetail = downsamp $ rconvolve x dwHiDecomp

idwt
    :: forall v l n a.
     ( UVG.Vector v (Complex a)
     , UVG.Vector v a
     , KnownNat l
     , KnownNat n
     , FFTWReal a
     , Div ((((Div ((n + l) - 1) 2 * 2) + l) - 1) - n) 2 <= ((((Div ((n + l) - 1) 2 * 2) + l) - 1) - n)
     )
    => DWavelet  v l a
    -> DWT v ((n + l - 1) `Div` 2) a
    -> VG.Vector v n a      -- todo: remove the choice
idwt DW{..} DWT{..} = VG.slice @_ @((((n+l-1)`Div`2)*2+l-1-n)`Div`2) @n @(((n+l-1)`Div`2)*2+l-1-n-((((n+l-1)`Div`2)*2+l-1-n)`Div`2)) Proxy $
                VG.zipWith (+) x y
  where
    x = rconvolve (upsamp dwtApprox) dwLoRecon
    y = rconvolve (upsamp dwtDetail) dwHiRecon

downsamp :: (UVG.Vector v a, KnownNat n) => VG.Vector v n a -> VG.Vector v (n `Div` 2) a
downsamp v = VG.generate $ \i -> v `VG.index` doubleFinite i
  where
    doubleFinite :: Finite (n `Div` 2) -> Finite n
    doubleFinite (Finite x) = Finite (x * 2 + 1)

upsamp :: (UVG.Vector v a, KnownNat n, Num a) => VG.Vector v n a -> VG.Vector v (n * 2) a
upsamp v = VG.generate $ \i -> case separateProduct @2 i of
    (j, k)
      | j == 0    -> v `VG.index` k
      | otherwise -> 0

-- upsamp :: forall v n a. (UVG.Vector v a, KnownNat n, Num a) => VG.Vector v (n `Div` 2) a -> VG.Vector v n a
-- upsamp v = VG.generate $ \(Finite i) -> case i `divMod` 2 of
--     (d, m)
--       | m == 0    -> v `VG.index` Finite d
--       | otherwise -> 0

    -- let (j, k) = case separateProduct @2 i
    -- in  if j == 0
    --       then

-- unhaarPass
--     :: (UVG.Vector v a, KnownNat n, Num a)
--     => Vector v (2 ^ n) a
--     -> Vector v (2 ^ n) a
--     -> Vector v (2 ^ (n + 1)) a
-- unhaarPass app det = VG.generate $ \i ->
--     let (j, k) = separateProduct @2 i
--         combiner
--           | j == 0    = (+)
--           | otherwise = (-)
--     in  combiner (app `VG.index` k) (det `VG.index` k)

haar :: (UVG.Vector v a, Floating a) => DWavelet v 2 a
haar = DW{..}
  where
    srtot = sqrt 2 / 2
    dwHiDecomp = VG.fromTuple (-srtot, srtot)
    dwLoDecomp = VG.fromTuple ( srtot, srtot)
    dwHiRecon  = VG.fromTuple ( srtot,-srtot)
    dwLoRecon  = VG.fromTuple ( srtot, srtot)

-- data DWD :: (Type -> Type) -> Nat -> Type -> Type where
--     DWDZ :: { _dwdDetailZ :: !a
--            , _dwdApprox  :: !a
--            } -> DWD v 0 a
--     DWDS :: { _dwdDetailS :: !(Vector  v (2 ^ (n + 1)) a)
--            , _dwdNext    :: !(DWD v n             a)
--            } -> DWD v (n + 1) a

-- deriving instance (Show (v a), Show a) => Show (DWD v n a)

-- dwdDetail :: UVG.Vector v a => DWD v n a -> Vector v (2 ^ n) a
-- dwdDetail = \case
--     DWDZ d _ -> VG.singleton d
--     DWDS d _ -> d

-- dwdApprox :: DWD v n a -> a
-- dwdApprox = \case
--     DWDZ _ a -> a
--     DWDS _ a -> dwdApprox a

-- testDec :: DWD UV.Vector 4 Double
-- testDec = DWDS (0 :: V.Vector 16 Double)
--         . DWDS (0 :: V.Vector  8 Double)
--         . DWDS (0 :: V.Vector  4 Double)
--         . DWDS (0 :: V.Vector  2 Double)
--         $ DWDZ 0 0

-- flattenDWD :: UVG.Vector v a => DWD v n a -> Vector v (2 ^ (n + 1)) a
-- flattenDWD = \case
--     DWDZ d  a  -> VG.fromTuple (a, d)
--     -- whoops this is O(n^2)
--     DWDS ds as -> flattenDWD as VG.++ ds

-- denseDWD
--     :: forall v n a. (KnownNat n, UVG.Vector v a)
--     => DWD v n a
--     -> V.Vector (n + 1) (Vector v (2 ^ n) a)
-- denseDWD = \case
--     DWDZ d  _  -> VG.singleton $ VG.singleton d
--     -- whoops this is O(n^3) or something
--     DWDS ds as -> (densify @_ @_ @(n - 1) <$> denseDWD as) `VG.snoc` ds

-- densify
--     :: (UVG.Vector v a, KnownNat n)
--     => Vector v (2 ^ n) a
--     -> Vector v (2 ^ (n + 1)) a
-- densify xs = VG.generate $ \i -> xs `VG.index` snd (separateProduct @2 i)

-- haar
--     :: forall v n a. (UVG.Vector v a, KnownNat n, Fractional a)
--     => Vector  v (2 ^ (n + 1)) a
--     -> DWD v n a
-- haar xs = case pNat (Proxy @n) of
--     PZ -> let x = xs `VG.index` 0
--               y = xs `VG.index` 1
--           in  DWDZ ((x - y) / 2) ((x + y) / 2)
--     PS -> let (a, d) = haarPass xs
--           in  DWDS d (haar a)

-- unHaar
--     :: (UVG.Vector v a, KnownNat n, Num a)
--     => DWD v n a
--     -> Vector v (2 ^ (n + 1)) a
-- unHaar = \case
--     DWDZ d  a  -> VG.fromTuple (a + d, a - d)
--     DWDS ds as -> unhaarPass (unHaar as) ds

-- haarPass
--     :: (UVG.Vector v a, KnownNat n, Fractional a)
--     => Vector v (2 ^ (n + 1)) a
--     -> (Vector v (2 ^ n) a, Vector v (2 ^ n) a)
-- haarPass xs = (app, det)
--   where
--     app = VG.generate $ \i -> ( xs `VG.index` combineProduct (0, i)
--                               + xs `VG.index` combineProduct (1, i)
--                               ) / 2
--     det = VG.generate $ \i -> ( xs `VG.index` combineProduct (0, i)
--                               - xs `VG.index` combineProduct (1, i)
--                               ) / 2

-- unhaarPass
--     :: (UVG.Vector v a, KnownNat n, Num a)
--     => Vector v (2 ^ n) a
--     -> Vector v (2 ^ n) a
--     -> Vector v (2 ^ (n + 1)) a
-- unhaarPass app det = VG.generate $ \i ->
--     let (j, k) = separateProduct @2 i
--         combiner
--           | j == 0    = (+)
--           | otherwise = (-)
--     in  combiner (app `VG.index` k) (det `VG.index` k)

max_stars_repo_name: mstksg/wavelets

max_issues_repo_name: mstksg/pure-wavelets

max_forks_repo_name: mstksg/pure-wavelets


text:
module Statistics
    ( runJSONFold
    , receivedCreatedF
    , module Statistics.Block
    , module Statistics.Chart
    , module Statistics.CSV
    , module Statistics.Focus
    , module Statistics.Graph
    , module Statistics.Histogram
    , module Statistics.MemPool
    , module Statistics.Report
    , module Statistics.Throughput
    , module Statistics.Tx
    ) where

import           Control.Foldl   (Fold (..))
import           Data.Map.Strict (Map)
import qualified Data.Map.Strict as M

import JSONLog
import Statistics.Block
import Statistics.Chart
import Statistics.CSV
import Statistics.Focus
import Statistics.Graph
import Statistics.Histogram
import Statistics.MemPool
import Statistics.Report
-- import Statistics.Relay
import Statistics.Throughput
import Statistics.Tx
import Types
import Universum
import Util.Pipes    (fold')

runJSONFold :: FilePath -> Fold IndexedJLTimedEvent a -> IO a
runJSONFold logDir fd = runParseLogs logDir $ fold' fd

receivedCreatedF :: Fold IndexedJLTimedEvent (Map TxHash (Maybe Timestamp))
receivedCreatedF = f <$> txFirstReceivedF <*> inBlockChainF
  where
    f :: Map TxHash Timestamp -> Map TxHash Timestamp -> Map TxHash (Maybe Timestamp)
    f rm cm = M.mapWithKey g rm
      where
        g :: TxHash -> Timestamp -> Maybe Timestamp
        g tx ts = maybe Nothing (\ts' -> Just $ ts' - ts) $ M.lookup tx cm

max_stars_repo_name: darkverg/cardano

max_issues_repo_name: darkverg/cardano

max_forks_repo_name: darkverg/cardano


text:
module Statistics.Dispersion where

import Statistics.Center

variance :: (Fractional a) => [a] -> a
variance vals = (sum $ zipWith (*) deviations deviations)/n
    where n = (fromIntegral $ length vals)
          mu = arithmeticMean vals
          deviations = map (\x -> x-mu) vals

stdev :: (Floating a) => [a] -> a
stdev vals = sqrt $ variance vals

max_stars_repo_name: hadrianhughes/Haskell

max_issues_repo_name: hadrianhughes/Haskell

max_forks_repo_name: hadrianhughes/Haskell


text:
{-# LANGUAGE FlexibleContexts #-}

module RandomUtil.RandomMatrix
  (
    selectRandomMatrixElement
  , randomMatrixIndex
  , randomizeRandomMatrixElement
  , addGaussianNoiseToMatrixElement
  , addGaussianNoiseToRandomMatrixElement
  , randomNMatrixElements
  , RandomMatrixMutation
  ) where

import System.Random
import Data.Random.Normal
import Matrix.State
import Numeric.LinearAlgebra.HMatrix hiding (corr)

import Util.Tuples

type RandomMatrixMutation g a = g -> Matrix a -> (Matrix a, g)
type RandomMatrixIndex g = (MatrixIndex, g)

randomMatrixIndex :: RandomGen g => g -> Matrix a -> RandomMatrixIndex g
randomMatrixIndex gen m = let (i, gen')  = randomR (0, (rows m) - 1) gen
                              (j, gen'') = randomR (0, (cols m) - 1) gen'
                           in ((i,j), gen'')

selectRandomMatrixElement :: (Element a, Indexable (Vector a) a, RandomGen g)
  => g -> Matrix a -> (a, g)
selectRandomMatrixElement gen m = let (i, gen') = randomMatrixIndex gen m
                                  in (getElement m i, gen')

randomizeRandomMatrixElement :: (Element a, Random a, RandomGen g)
  =>  (a,a) -> RandomMatrixMutation g a
randomizeRandomMatrixElement bounds gen m = let (i, gen')  = randomMatrixIndex gen m
                                                (e, gen'') = randomR bounds gen'
                                            in (replaceMatrixElement i e m, gen'')

type Mean = Double
type Sigma = Double
addGaussianNoiseToRandomMatrixElement :: RandomGen g =>  Mean -> Sigma -> RandomMatrixMutation g Double
addGaussianNoiseToRandomMatrixElement u s gen m =
  let (i,gen') = randomMatrixIndex gen m
  in addGaussianNoiseToMatrixElement i u s gen' m

addGaussianNoiseToMatrixElement :: RandomGen g
  => MatrixIndex -> Mean -> Sigma -> RandomMatrixMutation g Double
addGaussianNoiseToMatrixElement index u s gen m =
  let (sample, gen') = normal' (u,s) gen
  in (applyFMatrixElement (+ sample) index m, gen')

randomNMatrixElements :: (Element a, Random a, RandomGen g)
  => Int -> RandomMatrixMutation g a -> RandomMatrixMutation g a
randomNMatrixElements 0 _ gen m = (m, gen)
randomNMatrixElements n f gen m = let (m', gen') = f gen m
                                  in randomNMatrixElements (n-1) f gen' m'


text:
{-# LANGUAGE BangPatterns #-}

{- | Some extra facilities to work with 'Rand' monad and 'PureMT'
     random number generator.
-}

module Moo.GeneticAlgorithm.Random
    (
    -- * Random numbers from given range
      getRandomR
    , getRandom
    -- * Probability distributions
    , getNormal2
    , getNormal
    -- * Random samples and shuffles
    , randomSample
    , randomSampleIndices
    , shuffle
    -- * Building blocks
    , withProbability
    -- * Re-exports from random number generator packages
    , getBool, getInt, getWord, getInt64, getWord64, getDouble
    , runRandom, evalRandom, newPureMT
    , Rand, Random, PureMT
    ) where

import Control.Monad (liftM)
import Control.Monad.State.Strict (MonadState(..), State, runState, evalState)
import Data.Complex (Complex (..))
import System.Random (RandomGen, Random(..))
import System.Random.Mersenne.Pure64
import qualified System.Random.Shuffle as S
import qualified Data.Set as Set
import Data.Word (Word64)
import Data.Int (Int64)

-- | Monad for handling randomized computation. This uses the mersenne twister RNG.
type Rand = State PureMT

-- | Run a random computation using the generator @g@, returning the result
-- and the updated generator.
runRandom  :: Rand a -> PureMT -> (a, PureMT)
runRandom = runState

-- | Evaluate a random computation using the mersenne generator @g@.  Note that the
-- generator @g@ is not returned, so there's no way to recover the
-- updated version of @g@.
evalRandom :: Rand a -> PureMT -> a
evalRandom = evalState

-- | Yield a new boolean value from the generator.
getBool     :: Rand Bool
getBool     = fmap (<0) getInt

-- | Yield a new 'Int' value from the generator.
getInt      :: Rand Int
getInt      = state randomInt

-- | Yield a new 'Word' value from the generator.
getWord     :: Rand Word
getWord     = state randomWord

-- | Yield a new 'Int64' value from the generator.
getInt64    :: Rand Int64
getInt64    = state randomInt64

-- | Yield a new 53-bit precise 'Double' value from the generator.
getDouble :: Rand Double
getDouble = state randomDouble

-- | Yield a new 'Word64' value from the generator.
getWord64   :: Rand Word64
getWord64   = state randomWord64

-- | Yield a new randomly selected value of type @a@ in the range @(lo, hi)@.
-- See 'System.Random.randomR' for details.
getRandomR :: Random a => (a, a) -> Rand a
getRandomR range = state $ randomR range

-- | Yield a new randomly selected value of type @a@.
-- See 'System.Random.random' for details.
getRandom :: Random a => Rand a
getRandom = state random

-- | Yield two randomly selected values which follow standard
-- normal distribution.
getNormal2 :: Rand (Double, Double)
getNormal2 = do
  -- Box-Muller method
  u <- getDouble
  v <- getDouble
  let (c :+ s) = exp (0 :+ (2*pi*v))
  let r = sqrt $ (-2) * log u
  return (r*c, r*s)

-- | Yield one randomly selected value from standard normal distribution.
getNormal :: Rand Double
getNormal = fst `liftM` getNormal2

-- | Take at most n random elements from the list. Preserve order.
randomSample :: Int -> [a] -> Rand [a]
randomSample n xs =
  state $ \g -> select g n (length xs) xs []
  where
    select rng _ _ [] acc = (reverse acc, rng)
    select rng n m xs acc
        | n <= 0     = (reverse acc, rng)
        | otherwise  =
            let (k, rng') = randomR (0, m - n) rng
                (x:rest) = drop k xs
            in  select rng' (n-1) (m-k-1) rest (x:acc)

-- | Select @sampleSize@ numbers in the range from @0@ to @(populationSize-1)@.
-- The function works best when @sampleSize@ is much smaller than @populationSize@.
randomSampleIndices :: Int -> Int -> Rand [Int]
randomSampleIndices sampleSize populationSize =
    state $ \g ->
        let (sampleSet, g') = buildSampleSet g sampleSize Set.empty
        in  (Set.toList sampleSet, g')
  where
    buildSampleSet g n s
        | n <= 0 = (s, g)
        | otherwise =
            let (i, g') = randomR (0, populationSize-1) g
            in  if (i `Set.member` s)
                then buildSampleSet g' n s
                else buildSampleSet g' (n-1) (Set.insert i s)

-- | Randomly reorder the list.
shuffle :: [a] -> Rand [a]
shuffle xs = state $ randomShuffle xs (length xs)

-- | Given a sequence (e1,...en) to shuffle, its length, and a random
-- generator, compute the corresponding permutation of the input
-- sequence, return the permutation and the new state of the
-- random generator.
randomShuffle :: RandomGen gen => [a] -> Int -> gen -> ([a], gen)
randomShuffle elements len g =
    let (rs, g') = rseq len g
    in  (S.shuffle elements rs, g')
  where
  -- | The sequence (r1,...r[n-1]) of numbers such that r[i] is an
  -- independent sample from a uniform random distribution
  -- [0..n-i]
  rseq :: RandomGen gen => Int -> gen -> ([Int], gen)
  rseq n g = second lastGen . unzip $ rseq' (n - 1) g
      where
        rseq' :: RandomGen gen => Int -> gen -> [(Int, gen)]
        rseq' i gen
          | i <= 0    = []
          | otherwise = let (j, gen') = randomR (0, i) gen
                        in  (j, gen') : rseq' (i - 1) gen'
        -- apply a function on the second element of a pair
        second :: (b -> c) -> (a, b) -> (a, c)
        second f (x,y) = (x, f y)
        -- the last returned random number generator
        lastGen [] = g   -- didn't use the generator yet
        lastGen (lst:[]) = lst
        lastGen gens = lastGen (drop 1 gens)

-- |Modify value with probability @p@. Return the unchanged value with probability @1-p@.
withProbability :: Double -> (a -> Rand a) -> (a -> Rand a)
withProbability p modify x = do
  t <- getDouble
  if t < p
     then modify x
     else return x


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE ViewPatterns #-}

-- TODO: This module is in the wrong module (these are not linear programs)
{-|
Module      : Numeric.LinearProgramming.L1.Static
Description : Solve l1 programs with static vectors/matrices.
Copyright   : (c) Ryan Orendorff, 2020
License     : BSD3
Stability   : experimental
-}
module Numeric.LinearProgramming.L1.Static (
  l1Solve,
  l1SolveGT,
  l1SolveO,
  lInfSolveO,
  l1SolveU
)

where

import           Data.Maybe (fromJust)
import           GHC.TypeLits
import           Numeric.LinearAlgebra.Static
import qualified Numeric.LinearAlgebra as LA
import qualified Numeric.LinearProgramming.L1 as L1

-- Internal only. Could lead to runtime exceptions except in this case we trust
-- the output from the solvers is the correct dimensions.
mkL :: KnownNat n => LA.Vector ℝ -> R n
mkL = fromJust . create

-- | Solution in the L_1 norm, with L_1 regularization, of a linear system @Ax=b@.
--
-- @argmin_x  λ||x||_1 + ||Ax-b||_1@
l1Solve :: (KnownNat m, KnownNat n) =>
           Double -- ^ λ
        -> L m n  -- ^ A
        -> R m    -- ^ b
        -> R n    -- ^ x
l1Solve λ (extract -> a) (extract -> b) = mkL $ L1.l1Solve λ a b

-- | Solution in the L_1 norm, with L_1 regularization, of a system of linear inequalities @Ax>=b@.
--
-- @argmin_x  λ||x||_1 + ||step(b-Ax)||_1@
l1SolveGT :: (KnownNat m, KnownNat n) =>
             Double -- ^ λ
          -> L m n  -- ^ A
          -> R m    -- ^ b
          -> R n    -- ^ x
l1SolveGT λ (extract -> a) (extract -> b) = mkL $ L1.l1SolveGT λ a b

-- | L_1 solution of overconstrained system Ax=b.
--
-- @argmin_x ||Ax-b||_1@
l1SolveO :: (KnownNat m, KnownNat n) =>
            L m n -- ^ A
         -> R m   -- ^ b
         -> R n   -- ^ x
l1SolveO (extract -> a) (extract -> b) = mkL $ L1.l1SolveO a b

-- | L_inf solution of overconstrained system Ax=b.
--
-- @argmin_x ||Ax-b||_inf@
lInfSolveO :: (KnownNat m, KnownNat n) =>
              L m n -- ^ A
           -> R m   -- ^ b
           -> R n   -- ^ x
lInfSolveO (extract -> a) (extract -> b) = mkL $ L1.lInfSolveO a b

-- | L1 solution of underconstrained linear system Ax=b.
--
-- @argmin_x ||x||_1 such that Ax=b@
l1SolveU :: (KnownNat m, KnownNat n) =>
            L m n -- ^ A
         -> R m   -- ^ b
         -> R n   -- ^ x
l1SolveU (extract -> a) (extract -> b) = mkL $ L1.l1SolveU a b

max_stars_repo_name: ryanorendorff/convex

max_issues_repo_name: ryanorendorff/convex

max_forks_repo_name: ryanorendorff/convex


text:
module Lib
( Ising
, randModel
, runBatch
, runMC
, totalEnergy
, totalMagnetization
, flipRandSpin
, flipSpin
) where

import Model

import Data.Maybe
import Control.Monad.State
import System.Random (Random, random, randomR)
import Numeric.LinearAlgebra (accum, atIndex, toLists)

runBatch :: Show a => Int -> [Ising a] -> Ising ([Maybe [a]], IsingState)
runBatch steps prop = do
    props <- replicateM steps $ runMC prop
    state <- get
    return (props, state)

runMC :: [Ising a] -> Ising (Maybe [a])
runMC properties = do
    incrementStep
    state_i <- get
    dE <- flipRandSpin >>= getSpinEnergy >>= return . negate
    state_f <- get

    if dE <= 0
        then put $ state_f { nAccept = nAccept state_f + 1 }
    else do
        r <- randR (0,1)

        if r <= exp (-beta * dE)
            then put $ state_f { nAccept = nAccept state_f + 1 }
        else do
            put state_i            -- reset the state to before the spin flip
            r <- rand :: Ising Int -- increment RNG
            return ()

    state <- get
    if step state `mod` propFreq state == 0
        then liftM Just $ sequence properties
        else return Nothing

totalEnergy :: Ising Double
totalEnergy = do
    state <- get
    let n = dim state
    let indices = [(row, col) | row <- [0..n-1], col <- [0..n-1]]
    spinE <- liftM sum $ mapM getSpinEnergy indices
    let fieldE = sum $ map (* h state) (spins state)
    normalize $ spinE - fieldE

totalMagnetization :: Ising Double
totalMagnetization = do
        state <- get
        normalize $ sum (spins state)

getSpinEnergy :: (Int, Int) -> Ising Double
getSpinEnergy (row, col) = do
    state <- get
    neighbors <- getNeighbors (row, col)
    let spin = model state `atIndex` (row,col)
    return $ 2 * j state * spin * (sum neighbors)

getNeighbors :: (Int, Int) -> Ising [Spin]
getNeighbors (row, col) = do
    state <- get
    let n = dim state
    let m = model state
    let ns = map (atIndex m) [ ((row - 1) `mod` n, col)
                             , ((row + 1) `mod` n, col)
                             , (row, (col - 1) `mod` n)
                             , (row, (col + 1) `mod` n)
                             ]
    return ns
flipRandSpin :: Ising (Int, Int)
flipRandSpin = do
    state <- get
    let n = dim state
    row <- randR (0,n-1)
    col <- randR (0,n-1)
    flipSpin (row, col)
    return (row, col)

flipSpin :: (Int, Int) -> Ising ()
flipSpin (row, col) = do
    state <- get
    let spin = (model state) `atIndex` (row, col)
    let model' = accum (model state) (\a b -> a) [((row,col), -spin)]
    put $ state { model = model' }

-----
--
-- Helper functions
--
-----

beta = 1.0

normalize :: Fractional a => a -> Ising a
normalize x = do
    state <- get
    return $ x / (fromIntegral $ dim state)^2

spins :: IsingState -> [Spin]
spins state = concat . toLists $ model state

randR :: (Random a) => (a,a) -> Ising a
randR (a,b) = do
    state <- get
    let (r, rng') = randomR (a,b) $ rng state
    put $ state { rng = rng' }
    return r

rand :: (Random a) => Ising a
rand = do
    state <- get
    let (r, rng') = random $ rng state
    put $ state { rng = rng' }
    return r

incrementStep :: Ising Int
incrementStep = do
    state <- get
    let step' = step state + 1
    put $ state { step = step' }
    return step'


text:
{-# LANGUAGE DeriveDataTypeable, ScopedTypeVariables, NoMonomorphismRestriction, BangPatterns, FlexibleInstances, ViewPatterns #-}
module Math.Probably.RandIO where

import qualified Math.Probably.PDF as PDF
import Math.Probably.Sampler
import Math.Probably.MCMC hiding (traceIt)
import Math.Probably.StochFun
import Math.Probably.FoldingStats
import Math.Probably.Sampler
import qualified Control.Monad.State.Strict as S
import System.IO
import qualified Data.Packed.Matrix as L
import Control.Monad.Trans
import Control.Monad
import Data.IORef
import qualified Numeric.LinearAlgebra as L
import Math.Probably.NelderMead

import qualified Data.Vector.Storable as V
import Statistics.Autocorrelation

import Text.Printf
import Data.List
import Control.Monad.Par

import Debug.Trace

type RIO = S.StateT Seed IO

runRIO :: RIO a -> IO a
runRIO mx = do
  sd <- getSeedIO
  S.evalStateT mx sd

io :: IO a -> RIO a
io = liftIO

sample :: Sampler a -> RIO a
sample (Sam f) = do sd <- S.get
                    let (x,sd') = f sd
                    S.put sd'
                    return x

update :: IORef a -> (a->Sampler a) -> RIO ()
update rf sm = do
  x <- io $ readIORef rf
  newx <- sample $ sm x
  io $ writeIORef rf newx

runChainRIO :: Int -> (a->String) -> a -> (a->Sampler a) -> RIO [a]
runChainRIO n showit init sam = do
    seed <- S.get
    (nseed, xs) <- io $ go n seed init []
    S.put nseed
    return xs
     where go 0 s x xs = return (s, xs)
           go nn s x xs = do let (!xx, !ns) = unSam (sam x) s
                             when(nn `rem` chsz==0) $
                                 putStrLn $ show (((n-nn) `div` chsz)*2)++"%: " ++showit xx
                             go (nn-1) ns xx $ xx:xs
           chsz = n `div` 50

data AdaMetRunPars = AdaMetRunPars
     { nmTol :: Double,
       displayIt :: Maybe (L.Vector Double -> String),
       verboseNM :: Bool,
       amnsam :: Int,
       initw:: Int -> Double,
       minNM ::Int,
       maxNM ::Int }

defaultAM = AdaMetRunPars 0.5 Nothing False 1000 (const 0.02) 100 1000

traceit s x = trace (s++show x) x

getInitialCov :: Sampler (L.Vector Double)
              -> (L.Vector Double -> Double)
              -> RIO (Double, L.Vector Double, L.Matrix Double)
getInitialCov inisam posteriorL = go 500 where
  pdfv v = posteriorL v
  go n | n < 1
         = error $ "can't find suitable initial conditions. Try a narrower prior"
       | otherwise =  do
    initialv <- sample inisam
    let iniSimplex = genInitial (negate . pdfv) [] (const 0.02)  $ initialv
    if any (nanOrInf . snd) iniSimplex
       then go (n-1)
       else runNM n iniSimplex
  runNM n iniSimplex = do
    let finalSim = goNm (negate . pdfv) [] 0.01 100 1000 iniSimplex
        (maxPost,hess) = hessianFromSimplex (negate . pdfv) [] [] finalSim
    if any nanOrInf $ concat $ L.toLists hess
       then go (n-10)
       else getCov finalSim hess
  getCov finalSim hess = do
    let cor = case L.mbCholSH $ hess of
            Just _ -> L.inv hess
            Nothing -> case L.mbCholSH $ PDF.posdefify hess of
                         Just _ -> L.inv $ PDF.posdefify hess
                         Nothing -> L.diag $ L.mapVector (recip)
                                           $ L.takeDiag $ hess
        initV = fst $ head $ finalSim
    return (negate $ snd $ head $ finalSim,
            initV, cor)

laplaceApprox :: AdaMetRunPars -> PDF.PDF (L.Vector Double) -> [Int] -> [((Int, Int), Double)]
              -> L.Vector Double -> (L.Vector Double, Maybe (L.Matrix Double), Simplex)
laplaceApprox (AdaMetRunPars nmtol dispit verbnm nsam initw minNM maxNM) pdf isInt fixed init =
     let iniSim = genInitial (negate . pdf) isInt initw $ init
         finalSim = {-if verbnm then goNmVerbose (negate . pdf) isInt nmtol minNM maxNM iniSim
                              else -} goNm (negate . pdf) isInt nmtol minNM maxNM iniSim

         (maxPost,hess) = hessianFromSimplex (negate . pdf) isInt fixed finalSim
--     io $ print maxPost
         mbcor = case L.mbCholSH $ hess of
                   Just _ -> Just $ L.inv hess
                   Nothing -> case L.mbCholSH $ PDF.posdefify hess of
                                 Just _ -> Just $ L.inv $ PDF.posdefify hess
                                 Nothing -> Just $ L.diag $ L.mapVector (recip) $ L.takeDiag $ hess
         initV = fst $ head $ finalSim
     in (initV, mbcor, finalSim)

nmAdaMet :: AdaMetRunPars -> PDF.PDF (L.Vector Double) -> [Int] -> [((Int, Int), Double)]
            -> L.Vector Double -> RIO [L.Vector Double]
nmAdaMet (AdaMetRunPars nmtol dispit verbnm nsam initw minNM maxNM ) pdf isInt fixed init = do
     let iniSim = genInitial (negate . pdf) isInt initw $ init
     io $ print iniSim
     let finalSim =  goNm (negate . pdf) isInt nmtol minNM maxNM iniSim
     io $ print finalSim
     let (maxPost,hess) = hessianFromSimplex (negate . pdf) isInt fixed finalSim
--     io $ print maxPost
     io $ putStrLn "hessian"
     io $ print hess

     let mbcor = case L.mbCholSH hess of
                   Just _ -> Just $ L.inv hess
                   Nothing -> Nothing
     io $ putStrLn "maybe inverse hess"
     io $ print mbcor
     let initV = centroid finalSim
     let mbcorChol = mbcor >>= L.mbCholSH
     case (mbcor, mbcorChol) of
       (Just cor, Just _) ->  do --io $ putStrLn "chol cor"
                                 --io $ print $ L.inv cor
                                 let ampar = AMPar initV initV cor 2.4 (pdf initV) 0 0
                                 runAdaMetRIO nsam True ampar pdf
       _         -> do iniampar <- sample $ initialAdaMet 100 (const 5e-3) pdf initV
                       froampar <- runAndDiscard (nsam*2) (show . ampPar) iniampar $ adaMet False pdf
                       runAdaMetRIO (nsam*2) True froampar pdf

runAdaMetRIO :: Int -> Bool -> AMPar -> PDF.PDF (L.Vector Double) -> RIO [L.Vector Double]
runAdaMetRIO n freeze ampar pdf = do
    seed <- S.get
    (nseed, xs) <- io $ go n seed ampar []
    S.put nseed
    return xs
     where go 0 s amp vs = do print $ amp
                              return (s, reverse vs)
           go nn s amp vs = do let (!ampn, !ns) = unSam (adaMet freeze pdf amp) s
                               when(nn `rem` chsz==0) $
                                   putStrLn $ show (((n-nn) `div` chsz)*2)++"%: " ++showV (ampPar ampn) ++" LH="++printf "%.3g" (pdf (ampPar ampn)) ++ " accept="++acceptS ampn
                               go (nn-1) ns (ampn) $ (ampPar ampn):vs
           chsz = n `div` 50

runAdaMetRioESS :: Int -> Bool -> AMPar -> PDF.PDF (L.Vector Double) -> RIO [L.Vector Double]
runAdaMetRioESS want_ess freeze ampar pdf = do
    seed <- S.get
    (nseed, xs, amp) <- io $ go 200 seed ampar []
    S.put nseed
    return xs
     where go (0::Int) s amp vs = goChunks s amp vs
           go nn s amp vs = do let (!ampn, !ns) = unSam (adaMet freeze pdf amp) s
                               go (nn-1) ns (ampn) $ (ampPar ampn):vs

           goChunks s amp [] = do
              (nseed, xs, namp) <- go 200 s amp []
              goChunks nseed namp xs
           goChunks s amp xs = do
              let have_ess = min (realToFrac $ count_accept amp)
                                 $ calcESS  want_ess xs
                  drawn = length xs
              putStrLn $ show $ ampPar amp
              putStrLn $ "ESS="++show have_ess++" from "++show drawn
                       ++" drawn accept ratio="++acceptS amp
              if have_ess > realToFrac want_ess
                 then return (s, xs, amp)
                 else do
                   let need_ess = min (realToFrac $ want_ess `div` 5)
                                    $ max 1
                                    $ realToFrac want_ess - have_ess
                       samples_per_es = realToFrac drawn/have_ess
                       to_do = max 50 $ min 2000 $ samples_per_es * need_ess
                   putStrLn $ "now doing "++ show to_do
                   go (round to_do) s amp xs

runFixMetRioESS :: Double -> Int -> AMPar -> PDF.PDF (L.Vector Double) -> RIO [L.Vector Double]
runFixMetRioESS factor want_ess ampar pdf = do
    seed <- S.get
    (nseed, xs, amp) <- io $ go 200 seed ampar []
    S.put nseed
    return xs
     where go (0::Int) s amp vs = goChunks s amp vs
           go nn s amp vs = do let (!ampn, !ns) = unSam (fixedMet factor pdf amp) s
                               go (nn-1) ns (ampn) $ (ampPar ampn):vs

           goChunks s amp [] = do
              (nseed, xs, namp) <- go 200 s amp []
              goChunks nseed namp xs
           goChunks s amp xs = do
              let have_ess = min (realToFrac $ count_accept amp)
                             $ calcESS want_ess xs
                  drawn = length xs
              --putStrLn $ show $ ampPar amp
              putStrLn $ "ESS="++show have_ess++" from "++show drawn
                       ++" drawn accept ratio="++acceptS amp
              if have_ess > realToFrac want_ess
                 then return (s, xs, amp)
                 else do
                   let need_ess = min (realToFrac $ want_ess `div` 5)
                                    $ max 1
                                    $ realToFrac want_ess - have_ess
                       samples_per_es = realToFrac drawn/have_ess
                       to_do = max 50 $ min 2000 $ samples_per_es * need_ess * 1.2
                   putStrLn $ "now doing "++ show to_do
                   go (round to_do) s amp xs

runFixMetRioToFile :: Int -> Int -> String -> Double -> AMPar
                      -> PDF.PDF (L.Vector Double) -> RIO ()
runFixMetRioToFile samples thinn fileNm factor ampar pdf = do
    seed <- S.get
    h <- io $ openFile (fileNm) WriteMode
    nseed <- io $ go h samples seed ampar []
    S.put nseed
    return ()
     where go h (0::Int) s amp vs = return s
           go h nn s amp vs = do let (!ampn, !ns) = unSam (fixedMet factor pdf amp) s
                                 when (nn `rem` thinn == 0) $ do
                                     hPutStrLn h $ show $ L.toList $ ampPar ampn
                                 go h (nn-1) ns (ampn) $ (ampPar ampn):vs

runFixMetRio :: Double -> Int -> AMPar -> PDF.PDF (L.Vector Double) -> RIO [L.Vector Double]
runFixMetRio factor samples ampar pdf = do
    seed <- S.get
    (nseed, xs, amp) <- io $ go samples seed ampar []
    S.put nseed
    return xs
     where go (0::Int) s amp vs = return (s,vs, amp)
           go nn s amp vs = do let (!ampn, !ns) = unSam (fixedMet factor pdf amp) s
                               go (nn-1) ns (ampn) $ (ampPar ampn):vs

runFixMetRioBurn :: Double -> Int -> AMPar -> PDF.PDF (L.Vector Double) -> RIO AMPar
runFixMetRioBurn factor burn ampar pdf = do
    seed <- S.get
    (nseed, amp) <- io $ go burn seed ampar
    S.put nseed
    return amp
     where go (0::Int) s amp  = return (s,amp)
           go nn s amp  = do let (!ampn, !ns) = unSam (fixedMet factor pdf amp) s
                             go (nn-1) ns (ampn)

traceHead vs = trace ("head="++(show $ vs L.@> 0)) vs

most_es_per_s = 10

calcESS :: Int -> [L.Vector Double] -> Double
calcESS  want_ess mcmcOut
   | L.dim (head mcmcOut) <  want_ess * most_es_per_s * 2
      = calcESSprim mcmcOut
   | otherwise
      = let thinFactor = L.dim (head mcmcOut) `div` ( want_ess * most_es_per_s)
        in calcESSprim $ map (thinV $ thinFactor - 1) mcmcOut

thinV 0   =  id
thinV thin = V.ifilter $ \ix _ -> ix `mod` thin == 0

--thinV thin v = V.generate ((V.length v `div` (thin+1))+1) $ \i -> (V.!) v  (i*(thin))

calcESSprim ::  [L.Vector Double] -> Double
calcESSprim mcmcOut =
  let ndims = L.dim $ head mcmcOut
      len = realToFrac (length mcmcOut)
      acfs = mapP (\i->  V.sum $ V.takeWhile (>0.1) $ fst3 $ autocorrelation $ V.fromList $ map (L.@>i) mcmcOut) [0..ndims-1]
      ess = foldl1' min $ flip map acfs $ \acfsum-> (len/(1+2*acfsum)) -- realToFrac samples/(1+2*acfsum)
  in ess

mapP :: NFData b => (a-> b) -> [a] -> [b]
mapP f xs = runPar $
  forM xs (spawn . return . f) >>= mapM get

runAdaMetRIOtoFile :: Int -> Int -> String -> Bool -> AMPar -> PDF.PDF (L.Vector Double) -> RIO ()
runAdaMetRIOtoFile n thinn fileNm freeze ampar pdf = do
    seed <- S.get
    h <- io $ openFile (fileNm) WriteMode
    (nseed) <- io $ go h n seed ampar
    S.put nseed
    return ()
    --return xs
     where go h 0 s amp  = do print $ amp
                              return s --  return (s, reverse vs)
           go h nn s amp = do let (!ampn, !ns) = unSam (adaMet freeze pdf amp) s
                              when(nn `rem` chsz==0) $
                                   putStrLn $ show (((n-nn) `div` chsz)*2)++"%: " ++
                                              showV (ampPar ampn) ++" LH="++
                                              printf "%.3g" (pdf (ampPar ampn)) ++
                                              " accept="++acceptS ampn
                              when (nn `rem` thinn == 0) $ do
                                   hPutStrLn h $ show $ L.toList $ ampPar ampn
                              go h (nn-1) ns (ampn)
           chsz = n `div` 50

showV v = "<"++intercalate "," (map (printf "%.4g") $ L.toList v)++">"
acceptS ampar  | count ampar == 0 = "0/0"
               | otherwise = printf "%.3g" (rate::Double) ++ " ("++show yes++"/"++show total++")"where
   rate = realToFrac (yes) / realToFrac (total)
   yes = count_accept ampar
   total = count ampar

ampParRate ampar  | count ampar == 0 = 0
                  | otherwise = rate where
   rate = realToFrac (yes) / realToFrac (total)
   yes = count_accept ampar
   total = count ampar

runAdaMetRIOInterleaveInitial :: Int -> Bool -> L.Matrix Double -> AMPar -> PDF.PDF (L.Vector Double) -> RIO [L.Vector Double]
runAdaMetRIOInterleaveInitial n freeze cov ampar pdf = do
    seed <- S.get
    (nseed, xs) <- io $ go n seed ampar []
    S.put nseed
    return xs
     where go 0 s amp vs = do print $ amp
                              return (s, vs)
           go nn s amp vs = do let (!ampn, !ns) = unSam (adaMetInterleaveInitial freeze cov pdf amp) s
                               when(nn `rem` chsz==0) $
                                   putStrLn $ show (((n-nn) `div` chsz)*2)++"%: " ++show (ampPar ampn)-- ++" LH="++show (pdf (ampPar ampn))
                               go (nn-1) ns ampn $ (ampPar ampn):vs
           chsz = n `div` 50

runChainStat :: Int -> (Fold a b) -> a -> (a->Sampler a) -> RIO b
runChainStat n (F acc fldini k _) ini sam = do
    seed <- S.get
    let (nseed, stat) = go n seed fldini ini
    S.put nseed
    return stat
     where go 0 s y x = (s, k y)
           go nn s y x = let (xx, ns) = unSam (sam x) s
                         in go (nn-1) ns (acc y xx) xx

runAndWriteTo :: String -> Int -> (a -> String) -> a -> (a->Sampler a) -> RIO a
runAndWriteTo fnm n showit x sam = do
    h <- io $ openFile fnm WriteMode
    seed <- S.get
    (y,nseed) <- io $ runChainPrim h n showit seed x sam
    io $ hClose h
    S.put nseed
    return y

runChainPrim :: Handle -> Int -> (a -> String) -> Seed -> a -> (a->Sampler a) -> IO (a, Seed)
runChainPrim h n showit seed x sam = go n seed x
   where go 0 s x = return (x,s)
         go nn s x = do
            let (y, ns) = unSam (sam x) s
            hPutStrLn h $ showit y
            when(nn `rem` chsz==0) $
              putStrLn $ show (((n-nn) `div` chsz)*2)++"%: "++showit y
            go (nn-1) ns y
         chsz = n `div` 50

runOnce :: a -> (a->Sampler a) -> RIO a
runOnce init sam = runAndDiscardPrim 1 init sam

runAndDiscard :: Int -> (a->String) -> a -> (a->Sampler a) -> RIO a
--runAndDiscard = runAndDiscardPrim
runAndDiscard 0 _ x _ = return x
runAndDiscard n showit x sam = do
    go 50 x
 where chunksz = n `div` 50
       go 0 x = io (putStr "\n") >> return x
       go n x = do
          x' <- runAndDiscardPrim chunksz x sam
          io $ putStrLn $ show ((50-n)*2+2)++"%: "++showit x'
          go (n-1) x'

runMeanVar :: Floating b => Int -> (a->b) -> a -> (a->Sampler a) -> RIO (b,b)
runMeanVar n conv init sam = do
   xs <- sample $ runChainS n init sam
   return $ both meanF stdDevF `runStat` map conv xs

runAndDiscardPrim :: Int -> a -> (a->Sampler a) -> RIO a
runAndDiscardPrim n init sam = do
   seed <- S.get
   let (x,nseed) = runPrim seed n init sam
   S.put nseed
   return x

runPrim :: Seed -> Int -> a -> (a->Sampler a) -> (a,Seed)
runPrim s 0 x _ = (x,s)
runPrim s n x sam = let (x', s') = unSam (sam x) s
                    in runPrim s' (n-1) x' sam

metSample2P :: Show a => String -> (Double -> a -> a -> Sampler a)
               -> PDF.PDF (a,b) -> (Param a,b) -> Sampler (Param a)
metSample2P st prop p (par@(Param j t tt _ curw ada ini xi),y)
    = let accept xi pi pstar | notNanInf2 pi pstar =  min 1 $ exp (pstar - pi)
                             | otherwise = cond [(nanOrInf pi && nanOrInf pstar,
                                                        error $ "metropolisLnP "++st++" pi pi pstar :"++
                                                                show (pi,pstar)++"\n"++
                                                                show xi),
                                                 (nanOrInf pstar, -1), -- never accept
                                                 (nanOrInf pi, 2)] $ error "metropolisLn: the impossible happend"
      in  do
        let (nextw, nj, nt) = calcNextW ada curw j t tt
        u <-  unitSample
        xstar <- prop nextw ini xi
        let pstar = p (xstar,y)
        let pi = p (xi,y)
        return $ if u < accept par pi pstar
                      then Param (nj+1) (nt+1) (tt+1) pstar nextw ada ini xstar
                      else Param nj (nt+1) (tt+1) pi nextw ada ini xi
metSample2PC :: Show a => String -> (Double -> a -> a -> Sampler a)
               -> PDF.PDF (a,b) -> (Param a,b) -> Sampler (Param a)
metSample2PC st prop p (par@(Param j t tt lhi curw ada ini xi),y)
    = let accept xi pi pstar | notNanInf2 pi pstar =  min 1 $ exp (pstar - pi)
                             | otherwise = cond [(nanOrInf pi && nanOrInf pstar,
                                                        error $ "metropolisLnP "++st++" pi pi pstar :"++
                                                                show (pi,pstar)++"\n"++
                                                                show xi),
                                                 (nanOrInf pstar, -1), -- never accept
                                                 (nanOrInf pi, 2)] $ error "metropolisLn: the impossible happend"
      in  do
        let (nextw, nj, nt) = calcNextW ada curw j t tt
        u <-  unitSample
        xstar <- prop nextw ini xi
        let pstar = p (xstar,y)
        let pi = if notNanInf lhi then lhi else p (xi,y)
        return $ if u < accept par pi pstar
                      then Param (nj+1) (nt+1) (tt+1) pstar nextw ada ini xstar
                      else Param nj (nt+1) (tt+1) pi nextw ada ini xi

bestOfTwoCov inisam posterior = do
   set1@(postval1,_,_) <- getInitialCov inisam posterior
   set2@(postval2,_,_) <- getInitialCov inisam posterior
   if postval1>postval2 -- I THINK
      then return set1
      else return set2

mkAMPar init cov initp = AMPar init init cov 2.4 initp 10 5

initAdaMetFromCov nsam pdf initv retries cov = do
  --lift $ putStrLn $ "starting from existing cov; try number "++ show retries
  iniampar <- sample $ initialAdaMetFromCov nsam (pdf) initv
                                                 (PDF.posdefify $ cov)
  --lift $ print iniampar
  let rate = realToFrac (count_accept iniampar) / realToFrac nsam
  case () of
     _ | retries > 8 -> do lift $ putStrLn "initals ok."
                           return iniampar
     _ | rate > 0.5 -> initAdaMetFromCov nsam pdf initv (retries +1) $ L.scale 2 cov
     _ | rate > 0.40 -> initAdaMetFromCov nsam pdf initv (retries +1) $ L.scale 1.5 cov
     _ | rate < 0.025 ->  initAdaMetFromCov nsam pdf initv (retries +1) $ L.scale 0.1 cov
     _ | rate < 0.04 ->  initAdaMetFromCov nsam pdf initv (retries +1) $ L.scale 0.2 cov
     _ | rate < 0.12 ->  initAdaMetFromCov nsam pdf initv (retries +1) $ L.scale 0.3 cov
     _ | rate < 0.16 ->  initAdaMetFromCov nsam pdf initv (retries +1) $ L.scale 0.5 cov
     _ | rate < 0.20 ->  initAdaMetFromCov nsam pdf initv (retries +1) $ L.scale 0.8 cov
     _ | otherwise -> do lift $ putStrLn "initals ok."
                         return iniampar

max_stars_repo_name: glutamate/probably

max_issues_repo_name: glutamate/probably

max_forks_repo_name: glutamate/probably