text:
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# Language QuasiQuotes #-}

module MLP.App (
   runInIO,
   run,
   train,
) where

import MLP.Types
    ( curEpoch,
      curError,
      dataSet,
      epochs,
      hiddenLayers,
      learningRate,
      minErrorRate,
      network,
      numInputs,
      numOutputs,
      params,
      runApp,
      topology,
      testSet,
      trainSet,
      App,
      MonadFileSystem(printText),
      Parameters,
      Pattern(Pattern),
      State(State) )
import MLP.Config (readConf, mkDataSet)
import MLP.Network (AllCon, ArrList, Learn(..), Net, MLP(..), NetInstance, arrLast)
import Numeric.LinearAlgebra.Static (R)
import Control.Lens (view, use, (+=), (^.), (.=))
import Control.Monad (when, forM_, forM)
import qualified Data.Vector.Storable as V (snoc, sum)
import Control.Monad.Writer.Class (tell)
import PyF (fmt)
import qualified Data.Text.Lazy as TL
import Control.Monad.Random.Class ( MonadRandom )
import Data.Singletons ( withSomeSing )
import Data.Singletons.Prelude.Tuple ()
import Data.Singletons.Prelude.List ()
import Fcf ()
import Fcf.Data.List ()
import Data.Singletons.TypeLits ( KnownNat, Nat, SNat(..) )
import Data.Constraint ( Dict(..) )
import Numeric.Natural ()
import Data.Maybe ( fromJust )
import MLP.NatList ( List(List), SNatList(..) )

type NetConstraints (i :: Nat) (hs :: [Nat]) (o :: Nat) =
   (AllCon KnownNat (i ': o ': hs), AllCon KnownNat (Topo (Net i hs o)), Learn (Net i hs o))

netDict :: forall (i :: Nat) (o :: Nat) (hs :: [Nat]). (AllCon KnownNat (i ': o ': hs)) =>
           SNatList ('List hs) ->
           Dict (NetConstraints i hs o)
netDict =
   \case
      SNatNil -> Dict
      SNatCons _ SNatNil -> Dict
      SNatCons _ (SNatCons _ SNatNil) -> Dict
      SNatCons _ (SNatCons _ (SNatCons _ SNatNil)) -> Dict

runInIO :: FilePath -> IO ()
runInIO fl = do
   conf <- readConf fl
   let i  = conf ^. topology.numInputs
       hs = conf ^. topology.hiddenLayers
       o  = conf ^. topology.numOutputs

   withSomeSing i $ \(SNat :: SNat i) ->
      withSomeSing o $ \(SNat :: SNat o) ->
         withSomeSing (List (0:hs)) $
            \case
               (SNatCons _ snatlist) ->
                  case netDict @i @o snatlist of
                     (Dict :: Dict (NetConstraints i hs o)) -> run @i @hs @o @IO conf

run :: forall i hs o m. (NetConstraints i hs o, MonadFileSystem m, MonadRandom m) => Parameters -> m ()
run conf = do
   dt <- mkDataSet @i @o conf
   net <- create @(Net i hs o)

   let st = State dt conf net 0 0
       ((a, st'),log) = runApp train st

   mapM_ printText log

trainOnSet :: forall i o hs. NetConstraints i hs o => App i hs o Double
trainOnSet = do
   trainset <- use (dataSet.trainSet)
   lr <- use (params.learningRate)

   errorRates <- forM trainset $ \(Pattern i o) -> do
      net <- use network

      let outs   = netOut i net
          deltas = mkDeltas outs o net
          actOut = fromJust $ arrLast @o outs
          err = netError @(Net i hs o) actOut o

      network .= newNet lr i deltas outs net

      -- tell [fmt| \n ==> Pattern Error {err:.4f} \n |]
      -- tell [fmt| \n ==> Deltas {deltas:s} \n |]
      -- tell [fmt| \n ==> Expected Output {o:s} \n |]
      -- tell [fmt| \n ==> Actual Output {actOut:s} \n |]
      -- tell [fmt| \n Network Dump \n {net:s} \n|]

      return err

   let totalErr = sum errorRates / fromIntegral (length trainset)

   tell [[fmt| \n ==> Epoch Train Error {totalErr:f} \n |]]

   return totalErr

evaluateOnSet :: forall i o hs. NetConstraints i hs o => App i hs o Double
evaluateOnSet = do
   testset <- use (dataSet.testSet)
   net <- use network

   errorRates <- forM testset $ \(Pattern i o) -> do

      let outs   = netOut i net
          deltas = mkDeltas outs o net
          actOut = fromJust $ arrLast @o outs
          err = netError @(Net i hs o) actOut o

      return err

   let totalErr = sum errorRates / fromIntegral (length testset)
   curError .= totalErr

   tell [[fmt| \n ==> Epoch Test Error {totalErr:f} \n |]]

   return totalErr

loop :: App i hs o () -> App i hs o ()
loop trainAction = do
   epochTotal <- use (params.epochs)
   epochNow <- use curEpoch

   when (epochNow < epochTotal) $ do
      trainAction

      curError <- use curError
      minError <- use (params.minErrorRate)

      when (curError > minError) $ do
         loop trainAction

train :: NetConstraints i hs o => App i hs o ()
train = do
   tell ["Starting Training process...\n"]

   loop $ do
      epoch <- use curEpoch

      tell [[fmt| \n ==> Epoch {epoch:d} \n |]]

      trainOnSet
      evaluateOnSet

      curEpoch += 1

max_stars_repo_name: ovanr/Typed-MLP

max_issues_repo_name: ovanr/Typed-MLP

max_forks_repo_name: ovanr/Typed-MLP


text:
{-# LANGUAGE FlexibleContexts #-}
module STC.DFTFilter where

import           Control.Monad                as M
import           Control.Monad.IO.Class
import           Control.Monad.Parallel       as MP
import           Control.Monad.Trans.Resource
import           Data.Array.Repa              as R
import           Data.Binary
import           Data.ByteString.Lazy         as BL
import           Data.Complex
import           Data.Conduit
import           Data.Conduit.List            as C
import           Data.Int
import           Data.List                    as L
import           Data.Vector.Storable         as VS
import           Data.Vector.Unboxed          as VU
import           DFT.Plan
import           Filter.Utils
import           FokkerPlanck.Interpolation
import           System.IO
import           Utils.Array
import           Utils.Parallel               hiding ((.|))

{-# INLINE sourceFile #-}
sourceFile :: FilePath -> ConduitT () (VS.Vector (Complex Double)) (ResourceT IO) ()
sourceFile filePath = do
  h <- liftIO $ openBinaryFile filePath ReadMode
  sourceFunc h
  liftIO $ hClose h
  where
    sourceFunc handle = do
      flag <- liftIO $ hIsEOF handle
      if flag
        then sourceNull
        else do
          lenBS <- liftIO $ BL.hGet handle 8
          let len = fromIntegral (decode lenBS :: Int64) :: Int
          dataBS <- liftIO $ BL.hGet handle len
          yield . VS.fromList $ (decode dataBS :: [Complex Double])
          sourceFunc handle

{-# INLINE sinkFile #-}
sinkFile :: FilePath -> ConduitT ByteString Void (ResourceT IO) ()
sinkFile filePath = do
  h <- liftIO $ openBinaryFile filePath WriteMode
  sinkFunc h
  where
    sinkFunc handle = do
      x <- await
      case x of
        Nothing -> liftIO $ hClose handle
        Just y -> do
          let len = BL.length y
          liftIO $ BL.hPut handle (encode len)
          liftIO $ BL.hPut handle y
          sinkFunc handle

{-# INLINE interpolation #-}
interpolation ::
     (R.Source r Double)
  => DFTPlan
  -> (Double -> Double -> Double -> Double -> Int -> Int -> Complex Double)
  -> R.Array r DIM5 Double
  -> Int
  -> Int
  -> Double
  -> Double
  -> (VU.Vector Double)
  -> (VU.Vector Double)
  -> (VU.Vector Double)
  -> (VU.Vector Double)
  -> (Int, Int)
  -> IO ByteString
interpolation plan pinwheelFunc radialArr xLen yLen scaleFactor rMax thetaFreqs scaleFreqs theta0Freqs scale0Freqs (t, s) = do
  let tf = thetaFreqs VU.! t
      sf = scaleFreqs VU.! s
      pinwheelArr =
        traverse2
          (fromUnboxed (Z :. VU.length theta0Freqs) theta0Freqs)
          (fromUnboxed (Z :. VU.length scale0Freqs) scale0Freqs)
          (\(Z :. numTheta0Freq) (Z :. numScale0Freq) ->
             (Z :. numTheta0Freq :. numScale0Freq :. xLen :. yLen)) $ \ft0 fs0 (Z :. t0 :. s0 :. i :. j) ->
          pinwheelFunc
            (tf - ft0 (Z :. t0))
            (sf + fs0 (Z :. s0))
            rMax
            0
            (i - center xLen)
            (j - center yLen)
      filter =
        radialCubicInterpolation
          (R.slice radialArr (Z :. t :. s :. All :. All :. All))
          scaleFactor
          pinwheelArr
  filterF <-
    dftExecute
      plan
      (DFTPlanID
         DFT1DG
         [VU.length theta0Freqs, VU.length scale0Freqs, xLen, yLen]
         [2, 3]) .
    VU.convert . toUnboxed . computeS $
    filter
  return . encode . VS.toList $ filterF

{-# INLINE dftFilter2File #-}
dftFilter2File ::
     (R.Source r Double)
  => ParallelParams
  -> DFTPlan
  -> (Double -> Double -> Double -> Double -> Int -> Int -> Complex Double)
  -> R.Array r DIM5 Double
  -> Int
  -> Int
  -> Double
  -> Double
  -> (VU.Vector Double)
  -> (VU.Vector Double)
  -> (VU.Vector Double)
  -> (VU.Vector Double)
  -> FilePath
  -> IO ()
dftFilter2File parallelParams plan pinwheelFunc radialArr xLen yLen scaleFactor rMax thetaFreqs scaleFreqs theta0Freqs scale0Freqs filePath =
  runConduitRes $
  sourceList
    [ (t, s)
    | t <- [0 .. VU.length thetaFreqs]
    , s <- [0 .. VU.length scaleFreqs]
    ] .|
  parConduitIO
    parallelParams
    (interpolation
       plan
       pinwheelFunc
       radialArr
       xLen
       yLen
       scaleFactor
       rMax
       thetaFreqs
       scaleFreqs
       theta0Freqs
       scale0Freqs) .|
  sinkFile filePath

{-# INLINE convolutionHelper #-}
convolutionHelper ::
     DFTPlan
  -> Int
  -> Int
  -> Int
  -> Int
  -> VS.Vector (Complex Double)
  -> VS.Vector (Complex Double)
  -> IO (VU.Vector (Complex Double))
convolutionHelper plan numThetaFreq numScaleFreq xLen yLen vecF1 vecF2 =
  toUnboxed .
  sumS .
  sumS .
  rotate4D .
  rotate4D .
  fromUnboxed (Z :. numThetaFreq :. numScaleFreq :. xLen :. yLen) . VS.convert <$>
  dftExecute
    plan
    (DFTPlanID IDFT1DG [numThetaFreq, numScaleFreq, xLen, yLen] [2, 3])
    (VS.zipWith (*) vecF1 vecF2)

{-# INLINE convolveSlowly #-}
convolveSlowly ::
     ParallelParams
  -> FilePath
  -> DFTPlan
  -> R.Array U DIM4 (Complex Double)
  -> IO (R.Array U DIM4 (Complex Double))
convolveSlowly parallelParams filePath plan arr = do
  let (Z :. numThetaFreq :. numScaleFreq :. xLen :. yLen) = extent arr
  xs <-
    runConduitRes $
    sourceFile filePath .|
    parConduitIO
      parallelParams
      (convolutionHelper
         plan
         numThetaFreq
         numScaleFreq
         xLen
         yLen
         (VU.convert . toUnboxed $ arr)) .|
    C.consume
  return . fromUnboxed (extent arr) . VU.concat $ xs

{-# INLINE convolutionHelperSink #-}
convolutionHelperSink ::
     DFTPlan
  -> Int
  -> Int
  -> Int
  -> Int
  -> VU.Vector Double
  -> VS.Vector (Complex Double)
  -> ((Int, Int), VS.Vector (Complex Double))
  -> IO (VU.Vector (Complex Double))
convolutionHelperSink plan numThetaFreq numScaleFreq xLen yLen thetaFreqs vecF1 ((t, _), vecF2) = do
  let tf = thetaFreqs VU.! t
  arr <-
    fromUnboxed (Z :. numThetaFreq :. numScaleFreq :. xLen :. yLen) . VS.convert <$>
    dftExecute
      plan
      (DFTPlanID IDFT1DG [numThetaFreq, numScaleFreq, xLen, yLen] [2, 3])
      (VS.zipWith (*) vecF1 vecF2)
  return . toUnboxed . sumS . sumS . rotate4D . rotate4D . R.traverse arr id $ \f idx@(Z :. t0 :. _ :. i :. j) ->
    f idx * exp (0 :+ ((thetaFreqs VU.! t0) + tf) * pi)

{-# INLINE convolveSlowlySink #-}
convolveSlowlySink ::
     ParallelParams
  -> FilePath
  -> DFTPlan
  -> VU.Vector Double
  -> VU.Vector Double
  -> R.Array U DIM4 (Complex Double)
  -> IO (R.Array U DIM4 (Complex Double))
convolveSlowlySink parallelParams filePath plan thetaFreqs scaleFreqs arr = do
  let (Z :. numThetaFreq :. numScaleFreq :. xLen :. yLen) = extent arr
  xs <-
    runConduitRes $
    sourceFile filePath .|
    mergeSource
      (sourceList
         [ (t, s)
         | t <- [0 .. VU.length thetaFreqs]
         , s <- [0 .. VU.length scaleFreqs]
         ]) .|
    parConduitIO
      parallelParams
      (convolutionHelperSink
         plan
         numThetaFreq
         numScaleFreq
         xLen
         yLen
         thetaFreqs
         (VU.convert . toUnboxed $ arr)) .|
    C.consume
  return . fromUnboxed (extent arr) . VU.concat $ xs

ext: hs

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE NoImplicitPrelude #-}
{-# OPTIONS_GHC -Wno-type-defaults #-}
{-# OPTIONS_GHC -fno-warn-unused-imports #-}
{-# OPTIONS_GHC -fplugin Data.Singletons.TypeNats.Presburger #-}

module Algebra.Algorithms.Groebner.ZeroDimSpec where

import Algebra.Algorithms.Groebner
import Algebra.Algorithms.ZeroDim
import Algebra.Internal
import Algebra.Prelude.Core hiding ((%))
import Algebra.Ring.Ideal
import Algebra.Ring.Polynomial
import Algebra.Ring.Polynomial.Quotient
import Control.Monad
import Control.Monad.Random
import Data.Complex
import Data.Convertible (convert)
import qualified Data.Foldable as F
import qualified Data.Matrix as M
import Data.Maybe
import qualified Data.Sized as SV
import Data.Type.Natural.Lemma.Order
import qualified Data.Vector as V
import Numeric.Field.Fraction (Fraction, (%))
import Numeric.Natural
import Test.QuickCheck hiding (promote)
import qualified Test.QuickCheck as QC
import Test.Tasty
import Test.Tasty.ExpectedFailure
import Test.Tasty.HUnit
import Test.Tasty.QuickCheck
import Utils
import qualified Prelude as P

default (Fraction Integer, Natural)

asGenListOf :: Gen [a] -> a -> Gen [a]
asGenListOf = const

test_solveLinear :: TestTree
test_solveLinear =
  testGroup
    "solveLinear"
    [ testCase "solves data set correctly" $
        forM_ linSet $ \set ->
          solveLinear (M.fromLists $ inputMat set) (V.fromList $ inputVec set)
            @?= Just (V.fromList $ answer set)
    , testProperty "solves any solvable cases" $
        forAll (resize 10 arbitrary) $ \(MatrixCase ms) ->
          let mat = M.fromLists ms :: M.Matrix (Fraction Integer)
           in rank mat == M.ncols mat
                ==> forAll (vector (length $ head ms))
                $ \v ->
                  let ans = M.getCol 1 $ mat P.* M.colVector (V.fromList v)
                   in solveLinear mat ans == Just (V.fromList v)
    , ignoreTestBecause "Needs effective inputs for this" $
        testCase "cannot solve unsolvable cases" $
          pure ()
    ]

test_univPoly :: TestTree
test_univPoly =
  testGroup
    "univPoly"
    [ testProperty "produces monic generators of the elimination ideal" $
        withMaxSuccess 25 $
          QC.mapSize (const 4) $
            checkForTypeNat [2 .. 4] chk_univPoly
    ]

test_radical :: TestTree
test_radical =
  testGroup
    "radical"
    [ ignoreTestBecause "We can verify correctness by comparing with singular, but it's not quite smart way..." $
        testCase "really computes radical" $ pure ()
    ]

test_fglm :: TestTree
test_fglm =
  testGroup
    "fglm"
    [ testProperty "computes monomial basis" $
        withMaxSuccess 25 $
          mapSize (const 3) $
            checkForTypeNat [2 .. 4] $ \sdim ->
              case zeroOrSucc sdim of
                IsZero -> error "impossible"
                IsSucc k ->
                  withWitness (lneqSucc k) $
                    withWitness (lneqZero k) $
                      withKnownNat k $
                        forAll (zeroDimOf sdim) $ \(ZeroDimIdeal ideal) ->
                          let base =
                                reifyQuotient (mapIdeal (changeOrder Lex) ideal) $ \ii ->
                                  map quotRepr $ fromJust $ standardMonomials' ii
                           in stdReduced (snd $ fglm ideal) == stdReduced base
    , testProperty "computes lex base" $
        checkForTypeNat [2 .. 4] $ \sdim ->
          case zeroOrSucc sdim of
            IsZero -> error "impossible"
            IsSucc k -> withKnownNat k $
              withWitness (lneqSucc k) $
                withWitness (lneqZero k) $
                  forAll (zeroDimOf sdim) $ \(ZeroDimIdeal ideal) ->
                    stdReduced (fst $ fglm ideal)
                      == stdReduced (calcGroebnerBasisWith Lex ideal)
    , testProperty "returns lex base in descending order" $
        checkForTypeNat [2 .. 4] $ \sdim ->
          case zeroOrSucc sdim of
            IsZero -> error "impossible"
            IsSucc k -> withKnownNat k $
              case (lneqSucc k, lneqZero k) of
                (Witness, Witness) ->
                  forAll (zeroDimOf sdim) $ \(ZeroDimIdeal ideal) ->
                    isDescending (map leadingMonomial $ fst $ fglm ideal)
    ]

test_solve' :: TestTree
test_solve' =
  testGroup
    "solve'"
    [ testProperty "solves equation with admissible error" $
        withMaxSuccess 50 $
          mapSize (const 4) $ do
            checkForTypeNat [2 .. 4] $ chkIsApproximateZero 1e-5 (pure . solve' 1e-5)
    , testGroup "solves regression cases correctly" $
        map (approxZeroTestCase 1e-5 (pure . solve' 1e-5)) companionRegressions
    ]

test_solveViaCompanion :: TestTree
test_solveViaCompanion =
  testGroup
    "solveViaCompanion"
    [ testProperty "solves equation with admissible error" $
        withMaxSuccess 50 $
          mapSize (const 4) $
            checkForTypeNat [2 .. 4] $
              chkIsApproximateZero
                1e-5
                (pure . solveViaCompanion 1e-5)
    , testGroup "solves regression cases correctly" $
        map (approxZeroTestCase 1e-5 (pure . solveViaCompanion 1e-5)) companionRegressions
    ]

test_solveM :: TestTree
test_solveM =
  testGroup
    "solveM"
    [ testProperty "solves equation with admissible error" $
        withMaxSuccess 50 $
          mapSize (const 4) $
            checkForTypeNat [2 .. 4] $ chkIsApproximateZero 1e-9 solveM
    , testGroup "solves regressions correctly" $
        map (approxZeroTestCase 1e-9 solveM) companionRegressions
    ]

isDescending :: Ord a => [a] -> Bool
isDescending xs = and $ zipWith (>=) xs (drop 1 xs)

chkIsApproximateZero ::
  (KnownNat n) =>
  Double ->
  ( forall m.
    (0 < m, KnownNat m) =>
    Ideal (Polynomial (Fraction Integer) m) ->
    IO [Sized m (Complex Double)]
  ) ->
  SNat n ->
  Property
chkIsApproximateZero err solver sn =
  case zeroOrSucc sn of
    IsSucc _ ->
      withKnownNat sn $
        forAll (zeroDimOf sn) $ ioProperty . checkSolverApproxZero err solver
    IsZero -> error "chkIsApproximateZero must be called with non-zero typenats!"

companionRegressions ::
  [SolverTestCase]
companionRegressions =
  [ SolverTestCase @3 $
      let [x, y, z] = vars
       in ZeroDimIdeal $
            toIdeal
              [ 2 * x ^ 2 + (0.5 :: Fraction Integer) .*. x * y
              , - (0.5 :: Fraction Integer) .*. y ^ 2 + y * z
              , z ^ 2 - z
              ]
  , SolverTestCase @3 $
      let [x, y, z] = vars
       in ZeroDimIdeal $
            toIdeal
              [2 * x + 1, 2 * y ^ 2, - z ^ 2 - 1]
  , SolverTestCase @2 $
      let [x, y] = vars
       in ZeroDimIdeal $
            toIdeal
              [x ^ 2 - y, -2 * y ^ 2]
  ]

approxZeroTestCase ::
  Double ->
  ( forall n.
    (0 < n, KnownNat n) =>
    Ideal (Polynomial (Fraction Integer) n) ->
    IO [Sized n (Complex Double)]
  ) ->
  SolverTestCase ->
  TestTree
approxZeroTestCase err calc (SolverTestCase f@(ZeroDimIdeal i)) =
  testCase (show $ getIdeal f) $ do
    isZeroDimensional (F.toList i) @? "Not zero-dimensional!"
    checkSolverApproxZero err calc f
      @? "Solved correctly"

checkSolverApproxZero ::
  (0 < n, KnownNat n) =>
  Double ->
  (Ideal (Polynomial (Fraction Integer) n) -> IO [Sized n (Complex Double)]) ->
  ZeroDimIdeal n ->
  IO Bool
checkSolverApproxZero err solver (ZeroDimIdeal ideal) = do
  anss <- solver ideal
  let mul r d = convert r * d
  pure $
    all
      ( \as ->
          all ((< err) . magnitude . substWith mul as) $
            generators ideal
      )
      anss

data SolverTestCase where
  SolverTestCase ::
    (0 < n, KnownNat n) =>
    ZeroDimIdeal n ->
    SolverTestCase

chk_univPoly :: KnownNat n => SNat n -> Property
chk_univPoly sdim =
  forAll (zeroDimOf sdim) $ \(ZeroDimIdeal ideal) ->
    let ods = enumOrdinal sdim
     in conjoin $
          flip map ods $ \nth ->
            let gen = univPoly nth ideal
             in forAll (unaryPoly sdim nth) $ \f ->
                  (f `modPolynomial` [gen] == 0) == (f `isIdealMember` ideal)

rank :: (Ord r, Fractional r) => M.Matrix r -> Int
rank mat =
  let Just (u, _, _, _, _, _) = M.luDecomp' mat
   in V.foldr (\a acc -> if a /= 0 then acc + 1 else acc) (0 :: Int) $ M.getDiag u

data TestSet = TestSet
  { inputMat :: [[Fraction Integer]]
  , inputVec :: [Fraction Integer]
  , answer :: [Fraction Integer]
  }
  deriving (Show, Eq, Ord)

linSet :: [TestSet]
linSet =
  [ TestSet
      [ [1, 0, 0, 0, 0]
      , [0, (-2), (-2), (-2), (-2)]
      , [0, 0, 3 % 2, 0, (-1) % 2]
      , [0, 0, 0, 0, (-5) % 2]
      , [0, 1, 1, 1, 1]
      , [0, 0, (-2), 1, (-1)]
      ]
      [0, (-2), 19 % 5, 14 % 5, 1, 0]
      [0, (-81) % 25, 54 % 25, 16 % 5, (-28) % 25]
  ]

max_stars_repo_name: hangingman/computational-algebra

max_issues_repo_name: hangingman/computational-algebra

max_forks_repo_name: hangingman/computational-algebra


text:
module Main where

import qualified Numeric.LinearAlgebra as LA
import qualified MachineLearning.Types as T
import qualified MachineLearning as ML
import qualified MachineLearning.Classification.Binary as CB
import qualified MachineLearning.Classification.OneVsAll as OVA
import qualified MachineLearning.TerminalProgress as TP

processFeatures :: T.Matrix -> T.Matrix
processFeatures = ML.addBiasDimension . (ML.mapFeatures 2)

calcAccuracy :: T.Matrix -> T.Vector -> [T.Vector] -> Double
calcAccuracy x y thetas = OVA.calcAccuracy y yPredicted
  where yPredicted = OVA.predict x thetas

main = do
  putStrLn "\n=== Optical Recognition of Handwritten Digits Data Set ===\n"
  -- Step 1. Data loading.
  -- Step 1.1 Training Data loading.
  (x, y) <- pure ML.splitToXY <*> LA.loadMatrix "digits_classification/optdigits.tra"
  -- Step 1.1 Testing Data loading.
  (xTest, yTest) <- pure ML.splitToXY <*> LA.loadMatrix "digits_classification/optdigits.tes"

  -- Step 2. Outputs and features preprocessing.
  let numLabels = 10
      x1 = processFeatures x
      initialTheta = LA.konst 0 (LA.cols x1)
      initialThetas = replicate numLabels initialTheta
  -- Step 3. Learning.
  (thetas, optPath) <- TP.learnOneVsAllWithProgressBar (CB.learn (OVA.BFGS2 0.01 0.1) 0.001 30 (CB.L2 30) x1) y initialThetas 1

  -- Step 4. Prediction and checking accuracy
  let accuracyTrain = calcAccuracy x1 y thetas
      accuracyTest = calcAccuracy (processFeatures xTest) yTest thetas

  -- Step 5. Printing results.
  putStrLn $ "\nNumber of iterations to learn: " ++ show (map LA.rows optPath)

  putStrLn $ "Accuracy on train set (%): " ++ show (accuracyTrain*100)
  putStrLn $ "Accuracy on test set (%): " ++ show (accuracyTest*100)

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ScopedTypeVariables #-}

-- |
-- Module      :  Spiral.Util.Pretty.LaTeX
-- Copyright   :  (c) 2017-2020 Drexel University
-- License     :  BSD-style
-- Maintainer  :  mainland@drexel.edu

module Spiral.Util.Pretty.LaTeX (
    Pretty(..)
  ) where

import Data.Complex (Complex(..))
import Data.Complex.Cyclotomic
import Data.List (intersperse)
import qualified Data.Map as Map
import qualified Data.Matrix as M
import Data.Modular (Mod, unMod)
import Data.Monoid ((<>))
import Data.Ratio (Ratio, numerator, denominator)
import Data.Set (Set)
import qualified Data.Set as Set
import Data.Symbol (unintern)
import qualified Data.Text as T
import qualified Data.Vector as V
import GHC.TypeLits (KnownNat)
import Text.LaTeX
import Text.LaTeX.Base.Class (comm0,
                              comm1)
import Text.LaTeX.Packages.AMSFonts
import Text.LaTeX.Packages.AMSMath

import Spiral.Array (IArray,
                     Matrix,
                     (!))
import qualified Spiral.Array as A
import Spiral.Array.Shape
import Spiral.Exp
import Spiral.Permutation hiding (Inv)
import qualified Spiral.Permutation as P
import Spiral.SPL
import Spiral.Util.Name (Name(..))
import Spiral.Util.Pretty (Assoc(..),
                           Fixity(..),
                           HasFixity(..),
                           addPrec,
                           infixl_,
                           precOf)
import Spiral.Util.Uniq (Uniq(..))

class Pretty a where
  {-# MINIMAL pprPrec | ppr #-}
  ppr     :: a -> LaTeX
  pprPrec :: Int -> a -> LaTeX

  pprList     :: [a] -> LaTeX
  pprPrecList :: Int -> [a] -> LaTeX

  ppr       = pprPrec 0
  pprPrec _ = ppr

  pprPrecList _ = pprList
  pprList       = autoSquareBrackets . commasep . map ppr

parensIf :: Bool -> LaTeX -> LaTeX
parensIf False l = l
parensIf True l  = autoParens l

hsep :: LaTeX -> [LaTeX] -> LaTeX
hsep sep = mconcat . intersperse sep

commasep :: [LaTeX] -> LaTeX
commasep = hsep ","

mskip :: LaTeX -> LaTeX
mskip = comm1 "mskip"

nicefrac :: LaTeX -> LaTeX -> LaTeX
nicefrac e1 e2 = mempty ^: e1 <> (mskip "-2mu" <> "/" <> mskip "-1mu") !: e2
--nicefrac e1 e2 = e1 <> "/" <> e2

instance Pretty a => Pretty [a] where
  ppr     = pprList
  pprPrec = pprPrecList

instance Pretty a => Pretty (Set a) where
  ppr = autoBraces . commasep . map ppr . Set.toList

instance Pretty Bool where
  ppr True  = mathfrak "T"
  ppr False = mathfrak "F"

pprSignedIntegral :: (Integral a, Show a) => a -> LaTeX
pprSignedIntegral x
  | x < 0     = "-" <> pprSignedIntegral (-x)
  | otherwise = (fromString . show) x

pprRealFloat :: (RealFloat a, Show a) => a -> LaTeX
pprRealFloat x
    | isIntegral x = ppr (ceiling x :: Integer)
    | x < 0        = fromString $ showFloat (-x)
    | otherwise    = fromString $ showFloat x
  where
    isIntegral :: RealFrac a => a -> Bool
    isIntegral x = fromIntegral (ceiling x :: Integer) == x

instance Pretty Int where
  ppr = pprSignedIntegral

instance Pretty Integer where
  ppr = pprSignedIntegral

instance Pretty Float where
  ppr = pprRealFloat

instance Pretty Double where
  ppr = pprRealFloat

instance (Eq a, Num a, Pretty a) => Pretty (Ratio a) where
  ppr x
      | d == 1    = ppr n
      | otherwise = nicefrac (ppr n) (ppr d)
    where
      n = numerator x
      d = denominator x

-- | Pretty-print an imaginary number
pprIm :: (Eq a, Num a, Pretty a) => a -> LaTeX
pprIm 0    = mempty
pprIm 1    = "i"
pprIm (-1) = "-i"
pprIm i    = ppr i <> "i"

instance (Eq a, Num a, Pretty a) => Pretty (Complex a) where
  ppr (r :+ 0)    = ppr r
  ppr (0 :+ 1)    = "i"
  ppr (0 :+ (-1)) = "-i"
  ppr (0 :+ i)    = pprIm i
  ppr (r :+ i)    = case (T.unpack . render) (pprIm i) of
                      '-' : _ -> ppr r <> pprIm i
                      _ -> ppr r + pprIm i

instance Pretty Cyclotomic where
  pprPrec p (Cyclotomic n mp)
    | Map.null mp = "0"
    | null xs     = leadingTerm rat n ex
    | otherwise   = parensIf (p > addPrec) $
                    leadingTerm rat n ex <> mconcat (map (followingTerm n) xs)
    where
      ((ex,rat):xs) = Map.toList mp

      pprBaseExp :: Integer -> Integer -> LaTeX
      pprBaseExp n 1  = omega !: ppr n
      pprBaseExp n ex = omega !^ (ppr n, ppr ex)

      leadingTerm :: Rational -> Integer -> Integer -> LaTeX
      leadingTerm r _ 0 = ppr r
      leadingTerm r n ex
        | r == 1     = t
        | r == (-1)  = "-" <> t
        | r > 0      = ppr r <> t
        | r < 0      = "-" <> ppr (-r) <> t
        | otherwise  = mempty
        where
          t = pprBaseExp n ex

      followingTerm :: Integer -> (Integer, Rational) -> LaTeX
      followingTerm n (ex, r)
        | r == 1     = "+" <> t
        | r == (-1)  = "-" <> t
        | r > 0      = "+" <> ppr r <> t
        | r < 0      = "-" <> ppr (-r) <> t
        | otherwise  = mempty
        where
          t = pprBaseExp n ex

instance (Pretty i, KnownNat p) => Pretty (Mod i p) where
  ppr z = ppr (unMod z)

instance Pretty (Const a) where
  ppr (BoolC x)      = ppr x
  ppr (IntC x)       = ppr x
  ppr (IntegerC x)   = ppr x
  ppr (FloatC x)     = ppr x
  ppr (DoubleC x)    = ppr x
  ppr (RationalC x)  = ppr x
  ppr (ComplexC r i) = ppr (r :+ i)
  ppr (W _ 0 _)      = "1"
  ppr (W n 1 _)      = omega !: ppr n
  ppr (W n k _)      = omega !^ (ppr n, ppr k)
  ppr (CycC x)       = ppr x
  ppr (PiC k)        = ppr k <> pi_
  ppr (ModularC x)   = ppr x

instance Pretty Uniq where
  ppr (Uniq u) = ppr u

instance Pretty Name where
  ppr (Name sym Nothing)  = fromString (unintern sym)
  ppr (Name sym (Just u)) = fromString (unintern sym) !: ppr u

instance Pretty Var where
  ppr (Var n) = ppr n

instance (Pretty e, IArray r DIM2 e) => Pretty (Matrix r e) where
  ppr mat = bmatrix Nothing $ M.matrix m n f
    where
      Z :. m :. n = A.extent mat

      f :: (Int, Int) -> LaTeX
      f (i, j) = ppr (mat ! (i-1, j-1))

instance Pretty (Exp a) where
  pprPrec p (ConstE c) = pprPrec p c
  pprPrec p (VarE v)   = pprPrec p v

  pprPrec _ (UnopE Abs e) =
    autoBrackets "|" "|" (ppr e)

  pprPrec _ (UnopE Sqrt e) =
    tsqrt Nothing (ppr e)

  pprPrec p (UnopE op e) =
    unop p op e

  pprPrec p (BinopE op e1 e2) =
    infixop p op e1 e2

  pprPrec _ (IdxE ev eis) =
    ppr ev <> autoSquareBrackets (commasep (map ppr eis))

  pprPrec p (ComplexE er ei) =
    parensIf (p > addPrec) $
    ppr (er :+ ei)

  pprPrec _ (ReE e) =
    comm0 "Re" <> pprPrec 10 e

  pprPrec _ (ImE e) =
    comm0 "Im" <> pprPrec 10 e

  pprPrec p (BBinopE op e1 e2) =
    infixop p op e1 e2

  pprPrec p (IfE e1 e2 e3) =
    parensIf (p > 10) $
    operatorname (mathsf "if") <> ppr e1 <>
    operatorname (mathsf "then") <> ppr e2 <>
    operatorname (mathsf "else") <> ppr e3

instance Pretty Unop where
  ppr Neg    = "-"
  ppr Abs    = operatorname "abs"
  ppr Signum = operatorname "sgn"
  ppr Exp    = texp
  ppr Log    = tlog
  ppr Sqrt   = operatorname "sqrt"
  ppr Sin    = tsin
  ppr Cos    = tcos
  ppr Asin   = arcsin
  ppr Acos   = arccos
  ppr Atan   = arctan
  ppr Sinh   = tsinh
  ppr Cosh   = tcosh
  ppr Asinh  = operatorname "arsinh"
  ppr Acosh  = operatorname "arcosh"
  ppr Atanh  = operatorname "artanh"

instance Pretty Binop where
  ppr Add  = "+"
  ppr Sub  = "-"
  ppr Mul  = comm0 "cdot"
  ppr Quot = comm0 "quot"
  ppr Rem  = operatorname "rem"
  ppr Div  = operatorname "div"
  ppr Mod  = operatorname "mod"
  ppr FDiv = "/"

instance Pretty BBinop where
  ppr Eq = "="
  ppr Ne = comm0 "neq"
  ppr Lt = "<"
  ppr Le = comm0 "le"
  ppr Ge = comm0 "ge"
  ppr Gt = ">"

instance Pretty Permutation where
  ppr (L mn n)        =  "L" !^ (ppr n, ppr mn)
  ppr (J n)           = operatorname "J" !: ppr n
  ppr (CS n k)        = operatorname "S" !: ppr n <> autoParens (ppr k)
  ppr (CRT m n _ _)   = operatorname "CRT" !^ (ppr n, ppr m)
  ppr (CRT' m n _ _)  = operatorname "CRT'" !^ (ppr n, ppr m)
  ppr (Good m n _ _)  = operatorname "Good" !^ (ppr n, ppr m)
  ppr (Good' m n _ _) = operatorname "Good'" !^ (ppr n, ppr m)
  ppr (R p a)         = operatorname "R" !^ (ppr p, ppr a)
  ppr (P.Inv p _)     = ppr p ^: "-1"

data MatrixBinop = KOp
                 | DSOp
  deriving (Eq, Ord, Show)

instance HasFixity MatrixBinop where
  fixity KOp  = infixl_ 7
  fixity DSOp = infixl_ 6

instance Pretty MatrixBinop where
  ppr KOp  = comm0 "otimes"
  ppr DSOp = comm0 "oplus"

pprArgs :: Pretty a => [a] -> LaTeX
pprArgs = autoParens . commasep . map ppr

instance (Num e, Pretty e) => Pretty (SPL e) where
  pprPrec p m@E{}       = pprPrec p (toMatrix m)
  pprPrec _ (I n)       = operatorname "I" !: ppr n
  pprPrec _ (T e)       = pprPrec 10 e ^: comm0 "intercal"
  pprPrec _ (Inv e)     = pprPrec 10 e ^: "-1"
  pprPrec p (Pi pi)     = pprPrec p pi
  pprPrec _ (Rot alpha) = operatorname "R" !: ppr alpha
  pprPrec _ (Diag xs)   = operatorname "diag" <> pprArgs (V.toList xs)
  pprPrec p (KDiag _ e) = pprPrec p e

  pprPrec _ ((tl `Beside` tr) `Above` (bl `Beside` br)) =
    bmatrix Nothing $ M.matrix 2 2 f
      where
        f :: (Int, Int) -> LaTeX
        f (1, 1) = ppr tl
        f (1, 2) = ppr tr
        f (2, 1) = ppr bl
        f (2, 2) = ppr br
        f _      = error "can't happen"

  pprPrec _ (Above a b) =
    bmatrix Nothing $ M.matrix 2 1 f
      where
        f :: (Int, Int) -> LaTeX
        f (1, 1) = ppr a
        f (2, 1) = ppr b
        f _      = error "can't happen"

  pprPrec _ (Beside a b) =
    bmatrix Nothing $ M.matrix 1 2 f
      where
        f :: (Int, Int) -> LaTeX
        f (1, 1) = ppr a
        f (1, 2) = ppr b
        f _      = error "can't happen"

  pprPrec p (Kron a b)  = infixop p KOp a b
  pprPrec p (DSum a b)  = infixop p DSOp a b
  pprPrec _ (Prod a b)  = pprPrec 10 a <> pprPrec 10 b
  pprPrec _ (Circ xs)   = operatorname "circ" <> pprArgs (V.toList xs)
  pprPrec _ (Skew xs)   = operatorname "skew" <> pprArgs (V.toList xs)
  pprPrec _ (Toep xs)   = operatorname "toep" <> pprArgs (V.toList xs)
  pprPrec _ (Re a)      = comm0 "Re" <> autoParens (ppr a)
  pprPrec _ F2          = operatorname "F" !: "2"
  pprPrec _ (DFT n)     = operatorname "DFT" !: ppr n
  pprPrec _ (DFT' n)    = operatorname "DFT" !^ (ppr n, "-1")
  pprPrec _ (F n w)     = operatorname "DFT" !: ppr n <> autoParens (ppr w)
  pprPrec _ (F' n w)    = operatorname "DFT" !^ (ppr n, "-1") <> autoParens (ppr w)

unop :: (Pretty a, Pretty op, HasFixity op)
     => Int -- ^ precedence of context
     -> op  -- ^ operator
     -> a
     -> LaTeX
unop prec op x =
    parensIf (prec > precOf op) $
    ppr op <> pprPrec (precOf op) x

infixop :: (Pretty a, Pretty b, Pretty op, HasFixity op)
        => Int -- ^ precedence of context
        -> op  -- ^ operator
        -> a   -- ^ left argument
        -> b   -- ^ right argument
        -> LaTeX
infixop prec op l r =
    parensIf (prec > opPrec) $
    pprPrec leftPrec l <> ppr op <> pprPrec rightPrec r
  where
    leftPrec | opAssoc == RightAssoc = opPrec + 1
             | otherwise             = opPrec

    rightPrec | opAssoc == LeftAssoc = opPrec + 1
              | otherwise            = opPrec

    Fixity opAssoc opPrec = fixity op

max_stars_repo_name: mainland/hspiral

max_issues_repo_name: mainland/hspiral

max_forks_repo_name: mainland/hspiral


text:
module Statistics.Distribution.Empirical where

import File

import Data.List     (sort, partition)
import Data.Time.Clock
import Statistics.Distribution (Distribution, cumulative)

import qualified Data.IntMap as IM

newtype EmpiricalDistribution = EmpiricalDistribution (NominalDiffTime -> Double)

instance Distribution EmpiricalDistribution where
  cumulative (EmpiricalDistribution f) = f . realToFrac

cumulativeEmpiricalDistribution :: FileFormat -> [(Int, EmpiricalDistribution)]
cumulativeEmpiricalDistribution (FileFormat ls _ runs) =
  let es = IM.elems runs
      lData = [ map (sum . map fst)
                $ filter (\l -> length l == len)
                $ map ( reverse . dropWhile (not . snd) . reverse
                        . IM.elems . IM.filterWithKey (\k _ -> (k <= a))) es
              | (a, len) <- zip (map snd ls) [1..] ]
      fs = map toF lData
  in zip (map snd ls) fs

levelEmpiricalDistribution :: FileFormat -> [(Int, EmpiricalDistribution)]
levelEmpiricalDistribution ff =
  let (ls, lData) = unzip $ onlyValidSplits $ levelData ff
      fs = map toF lData
  in zip ls fs

toF :: [NominalDiffTime] -> EmpiricalDistribution
toF ls =
  let slData = sort ls
      tLen   = fromIntegral $ length slData
      f = \t -> let (lp, gp) = partition (<= t) slData
                in case (reverse lp, gp) of
                     ([], []) -> 0.0
                     ([], _ ) -> 0.0
                     (_,  []) -> 1.0
                     (y1:y0:_, y2:y3:_) ->
                       let m = mu (y1, y2) t
                           lLen = fromIntegral $ length lp
                           p = \o -> (lLen + o) / tLen
                       in hermiteInterpolate (p (-1), p 0, p 1, p 2) m
                     (l:_,  m:_ ) ->
                       let lLen = fromIntegral $ length lp
                       in linearInterpolate (l, m) (lLen / tLen, (lLen + 1) / tLen) t

  in EmpiricalDistribution f
  where
    linearInterpolate :: (NominalDiffTime, NominalDiffTime) -> (Double, Double) -> NominalDiffTime -> Double
    linearInterpolate domain (lRange, uRange) domainVal =
      let m    = mu domain domainVal
          notM = 1 - m
      in lRange * notM + uRange * m

    mu :: RealFrac a => (a, a) -> a -> Double
    mu (lDomain, uDomain) domainVal = realToFrac (domainVal - lDomain) / realToFrac (uDomain - lDomain)

    hermiteInterpolate :: (Double, Double, Double, Double) -> Double -> Double
    hermiteInterpolate (y0, y1, y2, y3) m =
      let m2 = m * m
          a0 = (-0.5) * y0 + 1.5 * y1 - 1.5 * y2 + 0.5 * y3
          a1 = y0 - 2.5 * y1 + 2 * y2 - 0.5 * y3
          a2 = (-0.5) * y0 + 0.5 * y2
          a3 = y1
      in a0 * m * m2 + a1 * m2 + a2 * m + a3

max_stars_repo_name: argiopetech/timer

max_issues_repo_name: argiopetech/timer

max_forks_repo_name: argiopetech/timer


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE TupleSections #-}
module Main where

import qualified Numeric.LinearAlgebra as Lina
import Numeric.LinearAlgebra hiding( Matrix, Vector )
import qualified Numeric.LinearAlgebra as Lina

import qualified Data.Csv as CSV
import Control.Monad.Except
import qualified Data.ByteString.Lazy as BS
import qualified Data.Vector as Vec
import Data.Char
import Foreign.C.Types( CInt )

type Matrix = Lina.Matrix Double
type Vector = Lina.Vector Double
type Label = CInt

type ErrT m a = ExceptT String m a

trainingDataFormat =
	CSV.defaultDecodeOptions

inputDataFormat =
	CSV.defaultDecodeOptions{
		CSV.decDelimiter = fromIntegral (ord ' ')
	}

-----------------------------------------------------------------
-- IO stuff:
-----------------------------------------------------------------

main :: IO ()
main =
	handleErrors $
	do
		mapM_
			(uncurry4 testWithData . uncurry testParamsFromLabels) $
			allPairs [3,5,7,8]
	where
		handleErrors x =
			do
				valOrErr <- runExceptT x
				either
					(\err -> putStrLn $ "ERROR: " ++ err)
					return
					valOrErr

testWithData :: FilePath -> FilePath -> Label -> Label -> ErrT IO ()
testWithData trainingFile1 trainingFile2 label1 label2 =
	do
		liftIO $ putStrLn $ concat $
			[ "----------------------------------------------\n"
			, "classifying to labels ", show [label1, label2]
			, " in files ", show [trainingFile1, trainingFile2]
			]
		[trainingSet1, trainingSet2] <-
			mapM (readData trainingDataFormat) $
			[ trainingFile1
			, trainingFile2
			]
		(inputLabels, inputData) <-
			prepareInputData (`elem` [fromIntegral label1, fromIntegral label2]) <$>
			readData inputDataFormat "resource/zip.test"
		let beta = calcBeta trainingSet1 trainingSet2
		let classified = classify (label1,label2) beta inputData
		let result =  calcClassificationQuality (cmap round $ inputLabels) classified
		liftIO $ putStrLn $
			descriptionString
				trainingSet1
				trainingSet2
				beta
				inputData
				classified
				result

readData :: CSV.DecodeOptions -> FilePath -> ErrT IO Matrix
readData fmtOpts path =
	(fromRows . Vec.toList) <$>
	(
		ExceptT $
		fmap (CSV.decodeWith fmtOpts CSV.NoHeader) $
		BS.readFile path
	)

-- (helpers: )
-----------------------------------------------------------------

prepareInputData :: (Double -> Bool) -> Matrix -> (Vector, Matrix)
prepareInputData filterRowsBy =
	(\rawInput ->
		(
			flatten $ (rawInput ?? (All, Take 1)),
			rawInput ?? (All, Drop 1)
		)
	)
	.
	fromLists
	.
	filter (filterRowsBy . head)
	.
	toLists

testParamsFromLabels x y =
	let
		[filePath1, filePath2] = map (("resource/train." ++) . show) [x,y]
	in
		(filePath1, filePath2, x, y)

descriptionString set1 set2 beta inputData classified result =
	unlines $
	[ concat $ ["set1 size:", show $ size set1]
	, concat $ ["set2 size:", show $ size set2]
	, concat $ ["beta size:", show $ size beta]
	, concat $ ["inputData size:", show $ size inputData]
	, concat $ ["result: --------------------"]
	, concat $ ["classification quality:", show $ result]
	]

-----------------------------------------------------------------
-- actual code:
-----------------------------------------------------------------

calcBeta :: Matrix -> Matrix -> Vector
calcBeta set1 set2 =
	flatten $ -- matrix to vector
	let
		x =
			prependOnes $
				set1
				===
				set2
		y =
			konst (-1) (count1,1)
			===
			konst 1 (count2,1)
		count1 = rows set1
		count2 = rows set2
		xTr_times_x_SAFE =
			let
				xTr_times_x = tr' x <> x
			in
				if det xTr_times_x > 0.01
				then xTr_times_x
				else
					xTr_times_x + (alpha * ident (rows xTr_times_x))
					where alpha = 0.01
	in
		inv xTr_times_x_SAFE <> tr' x <> y

calcClassificationQuality :: Lina.Vector Label -> Lina.Vector Label -> Double
calcClassificationQuality expected res =
	(/ fromIntegral (size res)) $
	sum $
	map (\x -> if x/=0 then 0 else 1) $
	toList $
	expected - res

classify :: (Label, Label) -> Vector -> Matrix -> Lina.Vector Label
classify (labelNeg, labelPos) beta input =
	cmap (assignLabels . sgn) $
		prependOnes input #> beta
	where
		sgn :: Double -> Double
		sgn x =
			let temp = signum x in
				if temp /= 0 then temp else 1
		assignLabels x
			| x < 0 = labelNeg
			| otherwise = labelPos

prependOnes m =
	konst 1 (rows m,1) ||| m

-----------------------------------------------------------------
-- utils:
-----------------------------------------------------------------

instance CSV.FromRecord Vector where
	parseRecord v =
		fmap fromList $ CSV.parseRecord v

-- return all pairs in a list
allPairs l =
	case l of
		(x:rest) ->
			map (x,) rest ++ allPairs rest
		_ -> []

uncurry4 f (a,b,c,d) = f a b c d


text:
{-# LANGUAGE ParallelListComp #-}
import Numeric.LinearAlgebra (linearSolve, toDense, (!), flatten)
import Data.Monoid ((<>), Sum(..))

rMesh n (ar, ac) (br, bc)
  | n < 2 = Nothing
  | any (\x -> x < 1 || x > n) [ar, ac, br, bc] = Nothing
  | otherwise = between a b <$> voltage
  where
    a = (ac - 1) + n*(ar - 1)
    b = (bc - 1) + n*(br - 1)

    between x y v = abs (v ! a - v ! b)

    voltage = flatten <$> linearSolve matrixG current

    matrixG = toDense $ concat [ element row col node
                               | row <- [1..n], col <- [1..n]
                               | node <- [0..] ]

    element row col node =
      let (Sum c, elements) =
            (Sum 1, [((node, node-n), -1)]) `when` (row > 1) <>
            (Sum 1, [((node, node+n), -1)]) `when` (row < n) <>
            (Sum 1, [((node, node-1), -1)]) `when` (col > 1) <>
            (Sum 1, [((node, node+1), -1)]) `when` (col < n)
      in [((node, node), c)] <> elements

    x `when` p = if p then x else mempty

    current  = toDense [ ((a, 0), -1) , ((b, 0),  1) , ((n^2-1, 0), 0) ]

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
{-# OPTIONS_GHC -fno-warn-orphans #-}
{-# LANGUAGE ScopedTypeVariables #-}
module Storage where

import qualified Data.ByteString.Lazy as BL
import Data.Csv
import qualified Data.Vector as V
import Data.Maybe
import GHC.Stack (HasCallStack)
import Data.Time.Clock
import Data.String.Utils (strip)

import Statistics.Sample

import Types

saveTiming :: Configuration -> RunInfo -> IO ()
saveTiming conf info =
  let scaleFactor = 1000
  in appendFile (storage conf)
    (show (round ((riTime info) * scaleFactor)::Int) ++ ", " ++ riCommand info ++ "\n")

-- | Load all runs; this may become expensive.
loadTiming :: Configuration -> IO [RunInfo]
loadTiming conf = do
  csvData <- BL.readFile (storage conf)
  case decode NoHeader csvData of
    Left err -> error (unlines [ "Could not decode storage."
                               , "  storage: "++storage conf
                               , "  error: "++err
                               ])
    Right v -> return $ V.toList ((\(time::NominalDiffTime, cmd :: String) -> RunInfo (strip cmd) time) <$> v)

instance FromField NominalDiffTime where
  parseField s = do
    count <- parseField s :: Parser Int
    pure (fromRational $ (toRational count) / 1000)

timesForCommand :: String -> [RunInfo] -> [NominalDiffTime]
timesForCommand cmd xs = mapMaybe maybeTime xs
  where
    maybeTime :: RunInfo -> Maybe NominalDiffTime
    maybeTime (RunInfo c t) | (strip c) == (strip cmd) = Just t
                            | otherwise = Nothing

statsForCommand :: HasCallStack => [NominalDiffTime] -> Stats
statsForCommand times =
  let realTimes :: [Double]
      realTimes = realToFrac <$> times
  in Stats
  { stCount = length times
  , stMax = maximum times
  , stMin = minimum times
  , stMean = mean $ V.fromList realTimes
  , stStdDev = stdDev $ V.fromList realTimes
  }

max_stars_repo_name: creswick/timeit

max_issues_repo_name: creswick/timeit

max_forks_repo_name: creswick/timeit


text:
{-# LANGUAGE OverloadedStrings #-}

module Statistics.Classification.ROC (
    ROC(..)
  , ROCPoint(..)
  , rocPoint
  , efficientROC
  , roc
  ) where

import Data.List (foldl')

import Data.Csv ((.:), (.=), DefaultOrdered(..), FromNamedRecord(..), ToNamedRecord(..), header, namedRecord)
import qualified Data.Vector as V

import Statistics.Classification.ConfusionMatrix (ConfusionMatrix, trueConfMatrix, truePosRate, falsePosRate)
import Statistics.Classification.Types (ClassificationResult(..), ClassificationScore(..), sortByScore)

data ROC = ROC {
    rocAUC :: !Double
  , rocCurve :: !(V.Vector ROCPoint)
  } deriving (Eq, Show, Read)

data ROCPoint = ROCPoint {
    rocTruePos :: {-# UNPACK #-}!Double
  , rocFalsePos :: {-# UNPACK #-}!Double
  , rocCutoff :: {-# UNPACK #-}!Double
  , rocConfusion :: !ConfusionMatrix
  } deriving (Eq, Show, Read)

rocPoint :: Double -> ConfusionMatrix -> ROCPoint
rocPoint cutoff cm = ROCPoint {
    rocTruePos = truePosRate cm
  , rocFalsePos = falsePosRate cm
  , rocCutoff = cutoff
  , rocConfusion = cm
  }

instance DefaultOrdered ROCPoint where
  headerOrder _ = header [
      "true_positive_rate"
    , "false_positive_rate"
    , "cutoff"
    ] <> headerOrder (undefined :: ConfusionMatrix)

instance ToNamedRecord ROCPoint where
  toNamedRecord rp = namedRecord [
      "true_positive_rate" .= rocTruePos rp
    , "false_positive_rate" .= rocFalsePos rp
    , "cutoff" .= rocCutoff rp
    ] <> toNamedRecord (rocConfusion rp)

instance FromNamedRecord ROCPoint where
  parseNamedRecord r = ROCPoint <$>
        r .: "true_positive_rate"
    <*> r .: "false_positive_rate"
    <*> r .: "cutoff"
    <*> parseNamedRecord r

data EffRocState = EffRocState {
    effRocCutoff :: {-# UNPACK #-}!Double
  , effRocTP :: {-# UNPACK #-}!Int
  , effRocFP :: {-# UNPACK #-}!Int
  , effRocTPprev :: {-# UNPACK #-}!Int
  , effRocFPprev :: {-# UNPACK #-}!Int
  , effRocOut :: ![ROCPoint]
  , effRocAUC :: {-# UNPACK #-}!Double
  } deriving (Eq, Show, Read)

emptyEffRocState :: EffRocState
emptyEffRocState = EffRocState inf 0 0 0 0 [] 0
  where inf = 1 / 0

trapezoidArea :: Int -> Int -> Int -> Int -> Double
trapezoidArea x1 x2 y1 y2 = base * height
  where
    base = fromIntegral . abs $ x2 - x1
    height = 0.5 * fromIntegral y2 + 0.5 * fromIntegral y1

addArea :: EffRocState -> EffRocState
addArea s0 = s0 { effRocAUC = effRocAUC s0 + area }
  where
    area = trapezoidArea (effRocFP s0) (effRocFPprev s0) (effRocTP s0) (effRocTPprev s0)

updRocState :: Int -> Int -> EffRocState -> ClassificationScore Bool -> EffRocState
updRocState p n s r
  | effRocCutoff s /= cutoff = nextExample . nextCutoff $ addArea s
  | otherwise = nextExample s
  where
    cutoff = classifiedScore r
    actual = classifiedExample r
    nextCutoff s0 = s0 { effRocCutoff = cutoff, effRocTPprev = effRocTP s0, effRocFPprev = effRocFP s0 }
    nextExample s0 =
      let s1 = if actual then s0 { effRocTP = effRocTP s0 + 1 } else s0 { effRocFP = effRocFP s0 + 1 }
          cm = trueConfMatrix (ClassificationResult p (effRocTP s1)) (ClassificationResult n (n - effRocFP s1))
      in s1 { effRocOut = rocPoint cutoff cm : effRocOut s1 }

-- O(n log n) ROC curve generation and AUC calculation
-- Algorithm 1 and 2 in
-- Fawcett, T., 2006. An Introduction to ROC Analysis. Pattern Recognition Letters 27 (861 - 874)
efficientROC :: Foldable f => Int -> Int -> f (ClassificationScore Bool) -> ROC
efficientROC p n = out . foldl' (updRocState p n) emptyEffRocState
  where
    finAUC effr
      | p == 0 || n == 0 = 0 -- degenerate case, no actual positive or negative examples.
      | otherwise = (effRocAUC effr + trapezoidArea n (effRocFPprev effr) p (effRocTPprev effr)) / fromIntegral (p * n)
    out s = ROC {
        rocAUC = finAUC s
      , rocCurve = V.fromList $ effRocOut s
      }

roc :: V.Vector (ClassificationScore Bool) -> ROC
roc v = efficientROC p n . V.reverse $ sortByScore v
  where
    p = V.length $ V.findIndices classifiedExample v
    n = V.length v - p

max_stars_repo_name: tsbattman/rochs

max_issues_repo_name: tsbattman/rochs

max_forks_repo_name: tsbattman/rochs


text:
{-# LANGUAGE BangPatterns , RankNTypes, GADTs, DataKinds #-}

{- | The 'Numerical.HBLAS.BLAS.Level3' module provides a fully general
yet type safe Level3 BLAS API.

When in doubt about the semantics of an operation,
consult your system's BLAS api documentation, or just read the documentation
for
<https://software.intel.com/sites/products/documentation/hpc/mkl/mklman/index.htm the Intel MKL BLAS distribution>

A few basic notes about how to invoke BLAS routines.

Many BLAS operations take one or more arguments of type 'Transpose'.
'Tranpose' has the following different constructors, which tell BLAS
routines what transformation to implicitly apply to an input matrix @mat@ with dimension @n x m@.

*  'NoTranspose' leaves the matrix @mat@ as is.

* 'Transpose' treats the @mat@ as being implicitly transposed, with dimension
    @m x n@. Entry @mat(i,j)@ being treated as actually being the entry
    @mat(j,i)@. For Real matrices this is also the matrix adjoint operation.
    ie @Tranpose(mat)(i,j)=mat(j,i)@

*  'ConjNoTranspose' will implicitly conjugate @mat@, which is a no op for Real ('Float' or 'Double') matrices, but for
'Complex Float' and 'Complex Double' matrices, a given matrix entry @mat(i,j)==x':+'y@
will be treated as actually being  @conjugate(mat)(i,j)=y':+'x@.

* 'ConjTranpose' will implicitly transpose and conjugate the input matrix.
ConjugateTranpose acts as matrix adjoint for both real and complex matrices.

The *gemm operations  work as follows (using 'sgemm' as an example):

* @'sgemm trLeft trRight alpha beta left right result'@, where @trLeft@ and @trRight@
are values of type 'Transpose' that respectively act on the matrices @left@ and @right@.

* the generalized matrix computation thusly formed can be viewed as being
@result = alpha * trLeft(left) * trRight(right) + beta * result@

the *gemv operations are akin to the *gemm operations, but with @right@ and @result@
being vectors rather than matrices.

the *trsv operations solve for @x@ in the equation @A x = y@ given @A@ and @y@.
The 'MatUpLo' argument determines if the matrix should be treated as upper or
lower triangular and 'MatDiag' determines if the triangular solver should treat
the diagonal of the matrix as being all 1's or not.  A general pattern of invocation
would be @'strsv' matuplo  tranposeMatA  matdiag  matrixA  xVector@.
A key detail to note is that the input vector is ALSO the result vector,
ie 'strsv' and friends updates the vector place.

-}

module Numerical.HBLAS.BLAS.Level3(
        dgemm
        ,sgemm
        ,cgemm
        ,zgemm

        ,chemm
        ,zhemm

        ,cherk
        ,zherk

        ,cher2k
        ,zher2k

        ,ssymm
        ,dsymm
        ,csymm
        ,zsymm

        ,ssyrk
        ,dsyrk
        ,csyrk
        ,zsyrk

        ,ssyr2k
        ,dsyr2k
        ,csyr2k
        ,zsyr2k

        ,strmm
        ,dtrmm
        ,ctrmm
        ,ztrmm

        ,strsm
        ,dtrsm
        ,ctrsm
        ,ztrsm
            ) where

import Numerical.HBLAS.UtilsFFI
import Numerical.HBLAS.BLAS.FFI.Level3
import Numerical.HBLAS.BLAS.Internal.Level3
import Control.Monad.Primitive
import Data.Complex

sgemm :: PrimMonad m=>  GemmFun Float  orient  (PrimState m) m
sgemm =  gemmAbstraction "sgemm"  cblas_sgemm_safe cblas_sgemm_unsafe (\x f -> f x )

dgemm :: PrimMonad m=>  GemmFun  Double orient  (PrimState m) m
dgemm = gemmAbstraction "dgemm"  cblas_dgemm_safe cblas_dgemm_unsafe  (\x f -> f x )

cgemm :: PrimMonad m=>  GemmFun (Complex Float) orient  (PrimState m) m
cgemm = gemmAbstraction "cgemm" cblas_cgemm_safe cblas_cgemm_unsafe  withRStorable_

zgemm :: PrimMonad m=>  GemmFun (Complex Double) orient  (PrimState m) m
zgemm = gemmAbstraction "zgemm"  cblas_zgemm_safe cblas_zgemm_unsafe withRStorable_

chemm :: PrimMonad m=>  HemmFun (Complex Float) orient (PrimState m) m
chemm = hemmAbstraction "chemm" cblas_chemm_safe cblas_chemm_unsafe withRStorable_

zhemm :: PrimMonad m=>  HemmFun (Complex Double) orient (PrimState m) m
zhemm = hemmAbstraction "zhemm" cblas_zhemm_safe cblas_zhemm_unsafe withRStorable_

cherk :: PrimMonad m=>  HerkFun Float (Complex Float) orient (PrimState m) m
cherk = herkAbstraction "cherk" cblas_cherk_safe cblas_cherk_unsafe (\x f -> f x)

zherk :: PrimMonad m=>  HerkFun Double (Complex Double) orient (PrimState m) m
zherk = herkAbstraction "zherk" cblas_zherk_safe cblas_zherk_unsafe (\x f -> f x)

cher2k :: PrimMonad m=>  Her2kFun Float (Complex Float) orient (PrimState m) m
cher2k = her2kAbstraction "cher2k" cblas_cher2k_safe cblas_cher2k_unsafe withRStorable_

zher2k :: PrimMonad m=>  Her2kFun Double (Complex Double) orient (PrimState m) m
zher2k = her2kAbstraction "zher2k" cblas_zher2k_safe cblas_zher2k_unsafe withRStorable_

ssymm :: PrimMonad m=>  SymmFun Float orient (PrimState m) m
ssymm = symmAbstraction "ssymm" cblas_ssymm_safe cblas_ssymm_unsafe (\x f -> f x)

dsymm :: PrimMonad m=>  SymmFun Double orient (PrimState m) m
dsymm = symmAbstraction "dsymm" cblas_dsymm_safe cblas_dsymm_unsafe (\x f -> f x)

csymm :: PrimMonad m=>  SymmFun (Complex Float) orient (PrimState m) m
csymm = symmAbstraction "csymm" cblas_csymm_safe cblas_csymm_unsafe withRStorable_

zsymm :: PrimMonad m=>  SymmFun (Complex Double) orient (PrimState m) m
zsymm = symmAbstraction "zsymm" cblas_zsymm_safe cblas_zsymm_unsafe withRStorable_

ssyrk :: PrimMonad m=>  SyrkFun Float orient (PrimState m) m
ssyrk = syrkAbstraction "ssyrk" cblas_ssyrk_safe cblas_ssyrk_unsafe (\x f -> f x)

dsyrk :: PrimMonad m=>  SyrkFun Double orient (PrimState m) m
dsyrk = syrkAbstraction "dsyrk" cblas_dsyrk_safe cblas_dsyrk_unsafe (\x f -> f x)

csyrk :: PrimMonad m=>  SyrkFun (Complex Float) orient (PrimState m) m
csyrk = syrkAbstraction "csyrk" cblas_csyrk_safe cblas_csyrk_unsafe withRStorable_

zsyrk :: PrimMonad m=>  SyrkFun (Complex Double) orient (PrimState m) m
zsyrk = syrkAbstraction "zsyrk" cblas_zsyrk_safe cblas_zsyrk_unsafe withRStorable_

ssyr2k :: PrimMonad m=>  Syr2kFun Float orient (PrimState m) m
ssyr2k = syr2kAbstraction "ssyr2k" cblas_ssyr2k_safe cblas_ssyr2k_unsafe (\x f -> f x)

dsyr2k :: PrimMonad m=>  Syr2kFun Double orient (PrimState m) m
dsyr2k = syr2kAbstraction "dsyr2k" cblas_dsyr2k_safe cblas_dsyr2k_unsafe (\x f -> f x)

csyr2k :: PrimMonad m=>  Syr2kFun (Complex Float) orient (PrimState m) m
csyr2k = syr2kAbstraction "csyr2k" cblas_csyr2k_safe cblas_csyr2k_unsafe withRStorable_

zsyr2k :: PrimMonad m=>  Syr2kFun (Complex Double) orient (PrimState m) m
zsyr2k = syr2kAbstraction "zsyr2k" cblas_zsyr2k_safe cblas_zsyr2k_unsafe withRStorable_

strmm :: PrimMonad m=> TrmmFun Float orient (PrimState m) m
strmm = trmmAbstraction "strmm" cblas_strmm_safe cblas_strmm_unsafe (\x f -> f x)

dtrmm :: PrimMonad m=> TrmmFun Double orient (PrimState m) m
dtrmm = trmmAbstraction "dtrmm" cblas_dtrmm_safe cblas_dtrmm_unsafe (\x f -> f x)

ctrmm :: PrimMonad m=> TrmmFun (Complex Float) orient (PrimState m) m
ctrmm = trmmAbstraction "ctrmm" cblas_ctrmm_safe cblas_ctrmm_unsafe withRStorable_

ztrmm :: PrimMonad m=> TrmmFun (Complex Double) orient (PrimState m) m
ztrmm = trmmAbstraction "ztrmm" cblas_ztrmm_safe cblas_ztrmm_unsafe withRStorable_

strsm :: PrimMonad m=> TrsmFun Float orient (PrimState m) m
strsm = trsmAbstraction "strsm" cblas_strsm_safe cblas_strsm_unsafe (\x f -> f x)

dtrsm :: PrimMonad m=> TrsmFun Double orient (PrimState m) m
dtrsm = trsmAbstraction "dtrsm" cblas_dtrsm_safe cblas_dtrsm_unsafe (\x f -> f x)

ctrsm :: PrimMonad m=> TrsmFun (Complex Float) orient (PrimState m) m
ctrsm = trsmAbstraction "ctrsm" cblas_ctrsm_safe cblas_ctrsm_unsafe withRStorable_

ztrsm :: PrimMonad m=> TrsmFun (Complex Double) orient (PrimState m) m
ztrsm = trsmAbstraction "ztrsm" cblas_ztrsm_safe cblas_ztrsm_unsafe withRStorable_

max_stars_repo_name: schnecki/hblas

max_issues_repo_name: schnecki/hblas

max_forks_repo_name: schnecki/hblas


text:
module Numeric.RBM where

import Control.Monad
import Control.Monad.IO.Class
import Math.Probable hiding (vector)
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Data
import Numeric.LinearAlgebra.HMatrix (norm_2)

data RBM = RBM
  { weights  :: !(Matrix Double)
  , weightsT :: !(Matrix Double)
  } deriving (Eq, Show)

pp :: RBM -> IO ()
pp (RBM w _) = disp 2 w

randMat :: Int -> Int -> RandT IO Double -> IO (Matrix Double)
randMat n p = mwc . fmap fromLists . replicateM n . listOf p

new :: Int -- ^ nb of visible units
    -> Int -- ^ nb of hidden units
    -> IO RBM
new nvisible nhidden = do
  w      <- randMat nvisible nhidden (normal 0 0.1)
  let w'  = asColumn (konst 0 nvisible) ||| w
      w'' = asRow (konst 0 (nhidden + 1)) === w'
  return $ RBM w'' (tr w'')

trainOnce :: Matrix Double -- ^ training data
          -> Double        -- ^ learning rate
          -> RBM
          -> IO RBM
trainOnce trainingData learningRate (RBM w wt) = do
  let dat            = asColumn (konst 1 n) ||| trainingData
      pos_hidden_act = dat <> w
      pos_hidden_p   = cmap logistic pos_hidden_act

  pos_hidden_st <- do
    randomM <- randMat n p (uniform 0 1)
    return $ step (pos_hidden_p - randomM)

  let pos_assocs = tr dat <> pos_hidden_p

      neg_visible_act = pos_hidden_st <> wt
      neg_visible_p   = cmap logistic neg_visible_act
      n'              = rows neg_visible_p
      neg_visible_p'  = asColumn (konst 1 n') ||| dropColumns 1 neg_visible_p
      neg_hidden_act  = neg_visible_p' <> w
      neg_hidden_p    = cmap logistic neg_hidden_act
      neg_assocs      = tr neg_visible_p' <> neg_hidden_p

      new_w           = (w +) . cmap (\x -> learningRate * x / fromIntegral n) $
                          pos_assocs - neg_assocs

      err             = norm_2 (dat - neg_visible_p')

  putStrLn $ "Error: " ++ show err
  return $ RBM new_w (tr new_w)

  where n = rows trainingData
        p = cols w

logistic :: Double -> Double
logistic x = 1 / (1 + exp (negate x))

trainN :: Matrix Double -- ^ training data
       -> Double        -- ^ learning rate
       -> Int           -- ^ epochs
       -> RBM
       -> IO RBM
trainN dat learningRate epochs = go epochs

  where go 0 rbm = return rbm
        go n rbm = do
          rbm' <- trainOnce dat learningRate rbm
          go (n - 1) rbm'

runVisible :: RBM
           -> Matrix Double      -- ^ sample rows of values for visible units
           -> IO (Matrix Double) -- ^ values of hidden units for the given visible units
runVisible (RBM w wt) inputData = do
  let dat        = asColumn (konst 1 n) ||| inputData
      hidden_act = dat <> w
      hidden_p   = cmap logistic hidden_act

  hidden_st <- do
    randomM <- randMat n p (uniform 0 1)
    return $ step (hidden_p - randomM)

  return $ dropColumns 1 hidden_st

  where n = rows inputData
        p = cols w

runHidden :: RBM
          -> Matrix Double      -- ^ sample rows of hidden units values
          -> IO (Matrix Double) -- ^ values of visible units for the given hidden ones
runHidden (RBM w wt) hiddenData = do
  let dat         = asColumn (konst 1 n) ||| hiddenData
      visible_act = dat <> wt
      visible_p   = cmap logistic visible_act

  visible_st <- do
    randomM <- randMat n p (uniform 0 1)
    return $ step (visible_p - randomM)

  return $ dropColumns 1 visible_st

  where n = rows hiddenData
        p = cols wt

daydream :: RBM
         -> Int                -- ^ number of samples we want to extract
         -> IO [Vector Double] -- ^ each row is a sample of the visible units
                               --   that the RBM has "daydreamed"
daydream (RBM w wt) nsamples = do
  sample0 <- fmap (fromList . (1:) . tail) . mwc $ listOf n (uniform 0 1)

  go sample0 nsamples

  where n = rows w
        p = cols w

        go :: Vector Double -> Int -> IO [Vector Double]
        go _v 0 = return []
        go v i = do
          let hidden_act = asRow v <> w
              hidden_p   = cmap logistic hidden_act

          hidden_st <- do
            randomV <- fmap vector . mwc $ listOf p (uniform 0 1)
            return $ step (head (toRows hidden_p) - randomV)

          let hidden_st' = asRow . vector . (1:) . tail . toList $ hidden_st
              visible_act = hidden_st' <> wt
              visible_p   = cmap logistic visible_act

          visible_st <- do
            randomV <- fmap vector . mwc $ listOf n (uniform 0 1)
            return $ step (head (toRows visible_p) - randomV)

          let res = vector . tail . toList $ visible_st
          fmap (res:) $ go visible_st (i - 1)

-- EXAMPLE

-- data given in:
-- http://blog.echen.me/2011/07/18/introduction-to-restricted-boltzmann-machines/

-- Harry Potter, Avatar, LOTR, Gladiator, Titanic, Glitter
samples :: Matrix Double
samples = fromLists
  [ [ 1, 1, 1, 0, 0, 0 ]
  , [ 1, 0, 1, 0, 0, 0 ]
  , [ 1, 1, 1, 0, 0, 0 ]
  , [ 0, 0, 1, 1, 1, 0 ]
  , [ 0, 0, 1, 1, 0, 0 ]
  , [ 0, 0, 1, 1, 1, 0 ]
  ]

new_user :: Matrix Double
new_user = fromLists [ [ 0, 0, 0, 1, 1, 0 ] ]

hiddens :: Matrix Double
hiddens = fromLists
  [ [0, 0]
  , [0, 1]
  , [1, 0]
  , [1, 1]
  ]

test :: Double -> Int -> IO ()
test rate k = do
  rbm <- new 6 2

  pp rbm
  trained <- trainN samples rate k rbm
  pp trained
  putStrLn "-------"
  putStrLn "VISIBLE"
  res <- runVisible trained new_user
  disp 2 res
  putStrLn "-------"
  putStrLn "HIDDEN"
  res' <- runHidden trained hiddens
  disp 2 res'

  putStrLn "-------"
  putStrLn "DAYDREAM"
  daydream trained 10 >>= mapM_ print

main :: IO ()
main = test 0.1 10000

max_stars_repo_name: alpmestan/rbm

max_issues_repo_name: alpmestan/rbm

max_forks_repo_name: alpmestan/rbm


text:
{-# LANGUAGE MultiParamTypeClasses #-}

module Q.Options.Bachelier (
    Bachelier(..)
  , euOption
  , eucall
  , euput
  , module Q.Options
) where
import           Data.Time                      ()
import           Q.Stochastic.Discretize        ()
import           Q.Stochastic.Process           ()
import           Q.Time                         ()
import           Statistics.Distribution        (cumulative, density)
import           Statistics.Distribution.Normal (standard)
import           Control.Monad.State
import           Data.Random                    (RVar, stdNormal)
import           Q.MonteCarlo
import           Q.Options
import           Q.Types

data Bachelier = Bachelier Forward Rate Vol deriving Show

-- | European option valuation with bachelier model.
euOption ::  Bachelier -> YearFrac -> OptionType -> Strike -> Valuation
euOption (Bachelier (Forward f) (Rate r) (Vol sigma)) (YearFrac t) cp (Strike k)
  = Valuation premium delta vega gamma where
    premium = Premium $ df * (q*(f - k)*n(q*d1) + sigma*sqrt(t)/sqrt2Pi * (exp(-0.5 *d1 * d1)))
    delta   = Delta   $ df * n (q * d1)
    vega    = Vega    $ df * (sqrt t) / sqrt2Pi * (exp (-0.5 * d1 * d1))
    gamma   = Gamma   $ (df/(sigma * (sqrt t)))*(recip sqrt2Pi)*(exp(-0.5 *d1 * d1))
    d1 = (f - k) / (sigma * sqrt(t))
    q = cpi cp
    sqrt2Pi = sqrt (2*pi)
    df =  exp $ (-r) * t
    n = cumulative standard

-- | see 'euOption'
euput b t =  euOption b t Put

-- | see 'euOption'
eucall b t = euOption b t Call

instance Model Bachelier Double where
  discountFactor (Bachelier _ r _) t1 t2 = return $ exp (scale dt r)
    where dt = t2 - t1

  evolve (Bachelier (Forward f) (Rate r) (Vol sigma)) (YearFrac t) = do
    (YearFrac t0, f0) <- get
    let dt = t - t0
    dW <- (lift stdNormal)::StateT (YearFrac, Double) RVar Double
    let ft = f0 * exp (r * dt) + sqrt(sigma*sigma/2*r * ((exp (2 * r * dt)) - 1)) * dW
    put (YearFrac t, ft)
    return ft

max_stars_repo_name: ghais/lowq

max_issues_repo_name: ghais/lowq

max_forks_repo_name: ghais/lowq


text:
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
module Statistics.Sample.Histogram.MagnitudeSpec where

import qualified Data.Vector.Unboxed as U
import           Data.Monoid
import           Data.Functor.Identity
import           Statistics.Sample.Histogram.Magnitude
import           Test.Hspec
import           Test.QuickCheck
import           Test.QuickCheck.Checkers
import           Test.QuickCheck.Classes

main :: IO ()
main = hspec spec

spec :: Spec
spec = parallel . describe "Histogram" $ do

  it "calculates a meaningful histogram" $ do
    let x = foldHist 1 [-3, 1, 4, 6, 9, 9.9, 4.032, 9.2]
        result = U.fromList [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 2, 0, 1, 0, 0, 3]
    histBuckets x `shouldBe` result

  describe "calculates magnitude" $ do
    it "is -1" $ property $ \(MagN x) -> histMagnitude (fromList 1 [x]) == -1
    it "is 0" $ property $ \(Mag0 x) -> histMagnitude (fromList 1 [x]) == 0
    it "is 1" $ property $ \(Mag1 x) -> histMagnitude (fromList 1 [x]) == 1
    it "is 2" $ property $ \(Mag2 x) -> histMagnitude (fromList 1 [x]) == 2

  it "blends equal resolutions and magnitudes" $ do
    let x = foldHist 1 [1, 4, 6, 9]
        y = foldHist 1 [1, 6]
        result = U.fromList [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 1, 0, 2, 0, 0, 1]
    histBuckets (x <> y) `shouldBe` result

  it "scales non equal magnitudes to the higher magnitude" $ do
    let x = foldHist 1 [1, 4, 6, 9]
        y = foldHist 1 [134, 666]
        result = U.fromList [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 1, 0, 0, 0, 0, 1, 0, 0, 0]
    histMagnitude (x <> y) `shouldBe` 2
    histBuckets (x <> y) `shouldBe` result

  it "blurs non equal resolutions to the lower resolution" $ do
    let x = foldHist 1 [1, 4, 6, 9]
        y = foldHist 2 [1, 6.6]
        result = U.fromList [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 1, 0, 2, 0, 0, 1]
    histResolution (x <> y) `shouldBe` 1
    histBuckets (x <> y) `shouldBe` result

  it "blurs and scales non equal resolution and magnitude" $ do
    let x = foldHist 2 [1, 4, 6, 9]
        y = foldHist 3 [-17, 1, 66]
        vec = (U.replicate 100 0)
        result = U.accum (+) vec [(82, 1)]
              <> U.accum (+) vec [(1, 2), (4, 1), (6, 1), (9, 1), (66, 1)]
    histBuckets (x <> y) `shouldBe` result

  describe "produces keys for given magnitudes and scales" $ do
    it "2-1" $
      keys (fromList 2 [1]) `shouldBe` U.fromList (reverse [0, -0.1 .. -9.9 ] <> [0, 0.1 .. 9.9 :: Double])
    it "2-10" $
      keys (fromList 2 [10]) `shouldBe` U.fromList ([-99, -98 .. 0] <> [0, 1 .. 99 :: Double])
    it "1-9" $
      keys (fromList 1 [9]) `shouldBe` U.fromList ([-9, -8 .. 0] <> [0, 1 .. 9 :: Double])

  it "allows insertion" $ do
    let x = foldHist 1 [1, 4, 6, 9]
        result = U.fromList [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1]
    histBuckets (insert (-3) x) `shouldBe` result

  testBatch $ monoid (undefined :: Histogram Double)

  it "is commutable" $
    property $ isCommut ((<>) :: Histogram Double -> Histogram Double -> Histogram Double)

  it "has decreasing resolution on mappend" $
    property $ \(x :: Histogram Double, y :: Histogram Double) ->
      histResolution (x <> y) == min (histResolution x) (histResolution y)

  it "has increasing magnitude on mappend" $
    property $ \(x :: Histogram Double, y :: Histogram Double) ->
      histMagnitude (x <> y) == max (histMagnitude x) (histMagnitude y)

  it "is equivalent to mappend and insert" $
    property $ \(xs :: [Double]) ->
      foldHist 2 xs == foldMap (foldHist 2 . Identity) xs

testBatch :: TestBatch -> Spec
testBatch (batchName, tests) = describe ("laws for: " ++ batchName) $
    foldr (>>) (return ()) (map (uncurry it) tests)

instance (RealFrac a, Floating a, Arbitrary a) => Arbitrary (Histogram a) where
  arbitrary = do
    ys <- arbitrary
    x <- choose (0, 3)
    pure $ fromList (fromIntegral (x :: Int)) (getNonEmpty ys :: [a])

instance EqProp a => EqProp (Histogram a) where
  x =-= y = eq x y

newtype MagN = MagN Double
  deriving (Show, Ord, Eq, Num)

instance Arbitrary MagN where
  arbitrary = MagN <$> choose (0.1, 0.9)

newtype Mag0 = Mag0 Double
  deriving (Show, Ord, Eq, Num)

instance Arbitrary Mag0 where
  arbitrary = Mag0 <$> choose (1, 9.9)

newtype Mag1 = Mag1 Double
  deriving (Show, Ord, Eq, Num)

instance Arbitrary Mag1 where
  arbitrary = Mag1 <$> choose (10, 99.9)

newtype Mag2 = Mag2 Double
  deriving (Show, Ord, Eq, Num)

instance Arbitrary Mag2 where
  arbitrary = Mag2 <$> choose (100, 999.9)

max_stars_repo_name: skedgeme/histogram-magnitude

max_issues_repo_name: skedgeme/histogram-magnitude

max_forks_repo_name: skedgeme/histogram-magnitude


text:
-- meanframe.hs
-- Compiles to executable.
-- We are passed a directory. We find any images in that directory,
-- average them, and save the average as average/average.png.
-- This tool is not used in the training or analysis of images.
import ImageToVector (loadImagesToVectors, vectorToImage)
import System.Environment (getArgs)
import Data.Maybe (catMaybes)
import System.Directory (getDirectoryContents)
import Data.Packed.Vector
import Numeric.LinearAlgebra.HMatrix (scale)
import Vision.Primitive.Shape (ix2)
import Vision.Image.Storage.DevIL
import Data.List (foldl')

main = do
  dir:_ <- getArgs
  let f = loadImagesToVectors 1280 720 . fmap ((dir++)"/"++)
  imgs <- fmap catMaybes $ f =<< getDirectoryContents dir
  let avgimg = vectorToImage (ix2 1280 720) $ avg imgs
  save PNG (dir ++ "/average/average.png") avgimg

avg :: [Vector Double] -> Vector Double
avg xs = scale (1 / fromIntegral num) sum
  where
  (num, sum) = foldl' f zero xs
  zero = (0, fromList . take (1280*720*3) . repeat $ 0)
  f (n,s) v = n `seq` s `seq` (n+1,v+s)


text:
module GenerateSet where

import Data.Word
import Data.List.Split
import Data.Scientific as Scientific
import Data.Complex
import Control.Monad
import Codec.Picture
import Data.List
import Data.Time.Clock
import Data.Time.Calendar
import Graphics.X11.Xlib
import System.Exit (exitWith, ExitCode(..))
import Control.Concurrent
import Data.Bits
import GHC.Conc (numCapabilities)

import Mand

generateSet mv = do
  step <- promptForImput mv "enter step"
  sr <- promptForImput mv "enter starting real"
  nr <- promptForImput mv "enter no of real points"
  si <- promptForImput mv "starting im"
  ni <- promptForImput mv "no of im"
  name <- promptForImput mv "Enter name:"

  let stepf = read step :: Double
  let srf = read sr :: Double -- starting real
  let nrf = read nr :: Double --number of real points
  let sif = read si :: Double -- starting imaginary
  let nif = read ni :: Double -- number of imaginary points

  d <- getCurrentTime
  desc <- promptForImput mv "description:"
  appendFile "mand/log.txt" $ "\n"++(take 19 $ show d)++", "++name++", "++desc

  --writePng name $ generateImage (\x y -> PixelRGB8 (fromIntegral 0) (fromIntegral 0) (fromIntegral $ mand $ ((read sr :: S) + (read (show x) :: S) * (read step :: S)) :+ ((read si :: S) + (read (show y) :: S) * (read step :: S)))) (read nr :: Int) (read ni :: Int)
  writePng ("mand/"++name) $ generateImage (\x y -> generatingFunction stepf srf sif (sif+ nif * stepf) x y ) (read nr :: Int) (read ni :: Int)
--writePng "mandv1" $ generateImage (\x y -> PixelRGB8 (fromIntegral 100) (fromIntegral 100) (fromIntegral 0)) 300 200
  print "completed"
  comment <- promptForImput mv "add comments:"
  appendFile "mand/log.txt" $ ", "++comment

generatingFunction stepf srf sif realend x y = normal stepf srf sif realend x y
  where
    m = general mand_iteration Nothing ((srf + (fromIntegral x) * stepf) :+ (realend - ((fromIntegral y) * stepf))) 255
    normal stepf srf sif realend x y = PixelRGB8 (fromIntegral 0)  (fromIntegral 0) (fromIntegral m)
    fullcolour stepf srf sif realend x y = PixelRGB8 p1 p2 p3
    max = 10
    scale x = round $ 255 * ((fromIntegral x) / (fromIntegral (max * max * max -1)))
    --p1 = scale $ fromIntegral $ max * max * (div m max* max)
    --p2 = scale $ fromIntegral $ max * max * (mod (div m max) max * max)
    --p3 = scale $ max * max * (fromIntegral $ mod m max)
    p1 = fromIntegral $ div (scaleUp 5000 m) 65536
    p2 = fromIntegral $ mod (div (scaleUp 5000 m) 255) 256
    p3 = fromIntegral $ mod (scaleUp 5000 m) 255

scaleUp :: Integral a => a -> a -> a
scaleUp max i = round $ 16777215 * fromIntegral i / fromIntegral max
    --m :: Int
    --m = mand $ Right $ fromFloatDigits (srf + (fromIntegral x) * stepf) :+ fromFloatDigits (sif + (fromIntegral y) * stepf)
    --m = mand $ ( (read :: String -> Scientific) (show (srf + (fromIntegral x) * stepf)) :=  (read :: String -> Scientific) (show (sif+ (fromIntegral y) * stepf))
    --m = if (mand $ (srf + (fromIntegral x) * stepf) :+ (sif + (fromIntegral y) * stepf) )== 255  then 0 else 255

--test = foldl (\acc x -> 'b':acc) []  [1..100]

colCorect1 :: Int -> Int
colCorect1 x = floor $ 256*(1 - exp(fromIntegral x/40))

colCorect2 :: Int -> Int
colCorect2 x = floor . (/255) . fromIntegral  $  x * x

promptForImput :: MVar Bool -> String ->  IO String
promptForImput mv message = do
  takeMVar mv
  print message
  input <- getLine
  putMVar mv True
  if input == "*" then promptForImput mv message else return input

max_stars_repo_name: AlwinHughes/fractal

max_issues_repo_name: AlwinHughes/fractal

max_forks_repo_name: AlwinHughes/fractal


text:
{-# LANGUAGE ExistentialQuantification #-}
{-# LANGUAGE ImpredicativeTypes #-}
{-# LANGUAGE RankNTypes #-}
{-# OPTIONS_GHC -fno-warn-incomplete-patterns #-}

module Data.Histogram.Util where

import Control.Foldl (fold)
import qualified Control.Foldl as L
import Data.Histogram (Histogram)
import qualified Data.Histogram as H
import Data.Histogram.Fill
import Data.Monoid
import qualified Data.Vector.Unboxed as V
import Statistics.Distribution
import qualified Statistics.Distribution.Combo as OC
import Statistics.Distribution.Exponential
import Statistics.Distribution.Gamma
import qualified Statistics.Distribution.Laplace as OL
import Statistics.Distribution.Normal
import Statistics.Distribution.Poisson

bin :: BinD
bin = binDn (-10.5) 1 10.5

hist :: BinD -> V.Vector Double -> H.Histogram BinD Double
hist b =
  fillBuilderVec (forceDouble -<< mkSimple b)

-- | make a percentile histogram
perc :: Histogram BinD Double -> Histogram BinD Double
perc h = H.map (\x -> x / H.sum h) h

-- reversing out mean and stdev
-- calculating mean using the trick that the Bin boundaries are constructed to be isomorphic to actual values
meanH :: Histogram BinD Double -> Double
meanH h =
  let t = H.sum h
  in H.bfoldl step begin h /
     t
  where begin = 0.0
        step b bv v = b + bv * v

meanSqH :: Histogram BinD Double -> Double
meanSqH h =
  let t = H.sum h
  in H.bfoldl step begin h /
     t
  where begin = 0.0
        step b bv v =
          b +
          (bv ^
           (2 :: Int)) *
          v

type Pdf = BinD -> [Double] -> V.Vector Double

pdfBin :: (Double -> Double) -> BinD -> V.Vector Double
pdfBin cdf h =
  V.map (\(x,y) -> cdf y - cdf x) binsV
  where binsV =
          binsList h :: V.Vector (Double,Double)

pdfBinNormal :: Pdf
pdfBinNormal b [] =
  pdfBin (cumulative $
          normalDistr 0.0 1.0)
         b
pdfBinNormal b [s] =
  pdfBin (cumulative $
          normalDistr 0.0 s')
         b
  where s' = max 1.0e-5 s
pdfBinNormal b (m:s:_) =
  pdfBin (cumulative $
          normalDistr m s')
         b
  where s' = max 1.0e-5 s

pdfBinLaplace :: Pdf
pdfBinLaplace b (r:l:sl:glue:_) =
  pdfBin (cumulative $
          OL.laplace r' l' sl' glue)
         b
  where l' = max 1.0e-5 l
        r' = max 1.0e-5 r
        sl' = max 1.0e-4 sl
pdfBinLaplace _ _ = mempty

pdfBinExponential :: Pdf
pdfBinExponential b (e:_) =
  pdfBin (cumulative $ exponential e') b
  where e' = max 1.0e-5 e
pdfBinExponential _ _ = mempty

pdfBinPoisson :: Pdf
pdfBinPoisson b (e:_) =
  pdfBin (cumulative $ poisson e') b
  where e' = max 1.0e-5 e
pdfBinPoisson _ _ = mempty

pdfBinGamma :: Pdf
pdfBinGamma b (m:k:v:_) =
  V.map (m *) $
  pdfBin (cumulative $
          gammaDistr k' v')
         b
  where k' = max 1.0e-5 k
        v' = max 1.0e-5 v
pdfBinGamma _ _ = mempty

pdfBinCombo :: Pdf
pdfBinCombo b (split:m:s:r:l:sl:_) =
  pdfBin (cumulative $
          OC.combo split' m s' r' l' sl')
         b
  where split' = max 1.0e-5 split
        s' = max 1.0e-5 s
        l' = max 1.0e-5 l
        r' = max 1.0e-5 r
        sl' = max 1.0e-4 sl
pdfBinCombo _ _ = mempty

err :: Histogram BinD Double -> Pdf -> [Double] -> Double
err h pdf p =
  do let pdfT = pdf (H.bins h) p
         pdfE = H.histData (perc h)
         errV =
           zipWith (-)
                   (V.toList pdfE)
                   (V.toList pdfT)
         err' = fold L.sum (fmap abs errV)
     err'


text:
{-# LANGUAGE BangPatterns #-}

module Data.Benchmarks.UnorderedCollections.Distributions
  ( makeRandomData
  , makeRandomVariateData
    -- * Workloads
  , insertWorkload
  , deleteWorkload
  , uniformLookupWorkload
  , exponentialLookupWorkload

  , loadOnly
  , loadAndUniformLookup
  , loadAndSkewedLookup
  , loadAndDeleteAll
  , loadAndDeleteSome
  , uniformlyMixed
  ) where

import qualified Control.Concurrent.Thread as Th
import           Control.DeepSeq
import           Control.Monad
import           Control.Monad.Reader
import           Control.Monad.Trans (liftIO)
import           Data.Benchmarks.UnorderedCollections.Types
import qualified Data.Vector as V
import qualified Data.Vector.Mutable as MV
import qualified Data.Vector.Unboxed as VU
import           Data.Vector (Vector)
import qualified Data.Vector.Algorithms.Shuffle as V
import           GHC.Conc (numCapabilities)
import           Statistics.Distribution
import           Statistics.Distribution.Exponential
import           System.Random.MWC

import           Criterion.Collection.Types

------------------------------------------------------------------------------
debug :: (MonadIO m) => String -> m ()
debug s = liftIO $ putStrLn s

------------------------------------------------------------------------------
makeRandomData :: (NFData k) =>
                  (GenIO -> IO k)
               -> Int
               -> WorkloadMonad (Vector (k,Int))
makeRandomData !genFunc !n = do
    rng <- getRNG
    debug $ "making " ++ show n ++ " data items"
    keys <- liftIO $ vreplicateM n rng genFunc
    let !v = keys `V.zip` vals
    let !_ = forceVector v
    debug $ "made " ++ show n ++ " data items"
    return $! v

  where
    vals      = V.enumFromN 0 n

------------------------------------------------------------------------------
makeRandomVariateData :: (Ord k, NFData k, Variate k) =>
                         Int
                      -> WorkloadMonad (Vector (k,Int))
makeRandomVariateData = makeRandomData uniform

------------------------------------------------------------------------------
insertWorkload :: (NFData k) => Vector (k,Int) -> Vector (Operation k)
insertWorkload = mapForce $ \(k,v) -> Insert k v

------------------------------------------------------------------------------
deleteWorkload :: (NFData k) => Vector (k,Int) -> Vector (Operation k)
deleteWorkload = mapForce $ \(k,_) -> Delete k

------------------------------------------------------------------------------
uniformLookupWorkload :: (NFData k) =>
                         Vector (k,Int)
                      -> Int
                      -> WorkloadMonad (Vector (Operation k))
uniformLookupWorkload !vec !ntimes = do
    rng <- getRNG
    debug $ "uniformLookupWorkload: generating " ++ show ntimes ++ " lookups"
    v <- liftIO $ vreplicateM ntimes rng f
    debug $ "uniformLookupWorkload: done"
    return v

  where
    !n = V.length vec
    f r = do
        idx <- pick
        let (k,_) = V.unsafeIndex vec idx
        return $ Lookup k
      where
        pick = uniformR (0,n-1) r

------------------------------------------------------------------------------
exponentialLookupWorkload :: (NFData k) =>
                             Double
                          -> Vector (k,Int)
                          -> Int
                          -> WorkloadMonad (Vector (Operation k))
exponentialLookupWorkload !lambda !vec !ntimes = do
    rng <- getRNG
    liftIO $ vreplicateM ntimes rng f
  where
    !dist = exponential lambda
    !n    = V.length vec
    !n1   = n-1
    !nd   = fromIntegral n

    f r = do
        x <- uniformR (0.1, 7.0) r
        let idx = max 0 . min n1 . round $ nd * density dist x
        let (k,_) = V.unsafeIndex vec idx
        return $! Lookup k

------------------------------------------------------------------------------
loadOnly :: (NFData k) =>
            (GenIO -> IO k)     -- ^ rng for keys
         -> WorkloadGenerator (Operation k)
loadOnly !genFunc !n = return $ Workload V.empty f
  where
    f _ = liftM insertWorkload $ makeRandomData genFunc n

------------------------------------------------------------------------------
loadAndUniformLookup :: (NFData k) =>
                        (GenIO -> IO k)  -- ^ rng for keys
                     -> WorkloadGenerator (Operation k)
loadAndUniformLookup !genFunc !n = do
    !vals <- makeRandomData genFunc n
    let !inserts = insertWorkload vals

    return $! Workload inserts $ uniformLookupWorkload vals

------------------------------------------------------------------------------
loadAndSkewedLookup :: (NFData k) =>
                       (GenIO -> IO k)  -- ^ rng for keys
                    -> WorkloadGenerator (Operation k)
loadAndSkewedLookup !genFunc !n = do
    !vals <- makeRandomData genFunc n
    let !inserts = insertWorkload vals
    return $! Workload inserts $ exponentialLookupWorkload 1.5 vals

------------------------------------------------------------------------------
loadAndDeleteAll :: (NFData k) =>
                    (GenIO -> IO k)     -- ^ key generator
                 -> WorkloadGenerator (Operation k)
loadAndDeleteAll !genFunc !n = do
    rng <- getRNG
    !vals <- makeRandomData genFunc n
    let !inserts = insertWorkload vals
    let !deletes = deleteWorkload $ V.shuffle rng vals

    return $ Workload inserts (const $ return deletes)

------------------------------------------------------------------------------
loadAndDeleteSome :: (NFData k) =>
                     (GenIO -> IO k)
                  -> WorkloadGenerator (Operation k)
loadAndDeleteSome !genFunc !n = do
    !vals <- makeRandomData genFunc n
    let !inserts = insertWorkload vals

    return $ Workload inserts $ f vals

  where
    f vals k = do
        rng <- getRNG
        return $ deleteWorkload $ V.take k $ V.shuffle rng vals

------------------------------------------------------------------------------
uniformlyMixed :: (NFData k) =>
                  (GenIO -> IO k)
               -> Double
               -> Double
               -> WorkloadGenerator (Operation k)
uniformlyMixed !genFunc !lookupPercentage !deletePercentage !n = do
    let !numLookups = ceiling (fromIntegral n * lookupPercentage)
    let !numDeletes = ceiling (fromIntegral n * deletePercentage)

    !vals <- makeRandomData genFunc n
    let !inserts = insertWorkload vals
    !lookups <- uniformLookupWorkload vals numLookups

    rng <- getRNG
    let !deletes = deleteWorkload $ V.take numDeletes $ V.shuffle rng vals
    let !out = V.shuffle rng $ V.concat [inserts, lookups, deletes]

    return $! Workload V.empty $ const $ return $ forceVector out

------------------------------------------------------------------------------
-- utilities
------------------------------------------------------------------------------
forceVector :: (NFData k) => Vector k -> Vector k
forceVector !vec = V.foldl' force () vec `seq` vec
  where
    force x v = x `deepseq` v `deepseq` ()

mapForce :: (NFData b) => (a -> b) -> Vector a -> Vector b
mapForce !f !vIn = let !vOut = V.map f vIn
                   in forceVector vOut

-- split a GenIO
splitGenIO :: GenIO -> IO GenIO
splitGenIO rng = VU.replicateM 256 (uniform rng) >>= initialize

-- vector replicateM is slow as dogshit.
vreplicateM :: Int -> GenIO -> (GenIO -> IO a) -> IO (Vector a)
vreplicateM n origRng act = do
    rngs <- replicateM numCapabilities (splitGenIO origRng)
    mv <- MV.new n
    let actions = map (f mv) (parts `zip` rngs)
    results <- liftM (map snd) $ mapM Th.forkIO actions
    _ <- sequence results
    V.unsafeFreeze mv

  where
    parts = partition (n-1) numCapabilities

    f mv ((low,high),rng) = do
        f' low
      where
        f' !idx | idx > high = return ()
                | otherwise = do
                                x <- act rng
                                MV.unsafeWrite mv idx x
                                f' (idx+1)

partition :: Int -> Int -> [(Int,Int)]
partition n k = ys `zip` xs
  where
    xs = map f [1..k]
    ys = 0:(map (+1) xs)
    f i = (i * n) `div` k

max_stars_repo_name: sanjaymsh/hashtables

max_issues_repo_name: sanjaymsh/hashtables

max_forks_repo_name: sanjaymsh/hashtables


text:
module VB where

import Numeric.SpecFunctions (digamma)
import Numeric.LinearAlgebra.HMatrix
       (Matrix, size, (!), fromRows, scalar, sumElements, toRows, reshape,
        flatten, tr, dot, (<>), ( #> ), flatten, asRow, matrix, fromList, Vector, cmap)
import Numeric.LinearAlgebra.Devel (zipVectorWith)

data LogType
    = Log  -- ^ log x
    | InvLog  -- ^ log (1 - x)

type ExpectationQZ = Matrix Double
type ExpectationQmu = (Matrix Double, Matrix Double)  -- ^ Matrix Log and Matrix InvLog
type QZ = Matrix Double
type Qmu = (Matrix Double, Matrix Double)  -- ^ Matrix A and matrix B

-- | Moments of logarithmically transformed random variables
expectationBeta :: Double  -- ^ alpha
                -> Double -- ^ beta
                -> LogType  -- ^ log x or log (1 - x)
                -> Double
expectationBeta alpha beta (Log) = digamma alpha -
    digamma (alpha + beta)
expectationBeta alpha beta (InvLog) = digamma beta -
    digamma (alpha + beta)

em :: Matrix Double
   -> Vector Double
   -> Vector Double
   -> Double
   -> Double
   -> Vector Double
em = go 0 (-1)
  where
    go lQthetaNew lQthetaOld matrixX vectorPi vectorAlpha a b
        | abs (lQthetaNew - lQthetaOld) <
              0.1 =
            vectorPi
        | otherwise =
            let (eQz,eQmu) = eStep matrixX initQZ initQmu
                vectorPiNew = mStep vectorPi vectorAlpha eQz
            in go
                   (lQTheta matrixX vectorPiNew vectorAlpha eQz eQmu)
                   lQthetaNew
                   matrixX
                   vectorPiNew
                   vectorAlpha
                   a
                   b
        where initQZ = fromRows (replicate nRow vectorPi)
              initQmu = (matrixAB a, matrixAB b)
              matrixAB v = matrix
                      nColumn
                      [v | _ <-
                              [1 .. nMixtures * nColumn]]
              (nRow,nColumn) = size matrixX

eStep :: Matrix Double -> QZ -> Qmu -> (ExpectationQZ, ExpectationQmu)
eStep matrixX qZ qMu = let eQmu = calculateEQmu qMu
                           eQz = calculateEQZ qZ
                           qZ' = calculateQZ matrixX eQmu qZ
                           eQz' = calculateEQZ qZ'
                           qMu' = calculateQmu matrixX eQz' qMu
                           eQmu' = calculateEQmu qMu'
    in go
           (lQ matrixX eQz' eQmu')
           (lQ matrixX eQz eQmu)
           qZ'
           eQz'
           qMu'
           eQmu'
    where go lQnew lQold qZn eQzn qMun eQmun
              | abs (lQnew - lQold) <
                    0.1 =
                  (eQzn, eQmun)
              | otherwise =
                  let qZn' = calculateQZ matrixX eQmun qZn
                      eQzn' = calculateEQZ qZn'
                      qMun' = calculateQmu matrixX eQzn' qMun
                      eQmun' = calculateEQmu qMun'
                  in go (lQ matrixX eQzn' eQmun') lQnew qZn' eQzn' qMun' eQmun'

lQ :: Matrix Double -> ExpectationQZ -> ExpectationQmu -> Double
lQ matrixX eQz (eQmuM,eQmuMi) = sum
        (map
             (\n ->
                   (eQz ! n) `dot`
                   flatten
                       (asRow (matrixX ! n) <>
                        teQmuM +
                        asRow
                             (1 -
                              (matrixX ! n)) <>
                        teQmuMi))
             [0 .. nRow - 1])
    where teQmuM = tr eQmuM
          teQmuMi = tr eQmuMi
          (nRow,_) = size matrixX

calculateEQZ :: QZ -> ExpectationQZ
calculateEQZ qZ = fromRows
        (map
             (\row ->
                   row /
                   scalar (sumElements row))
             (toRows qZ))

calculateQZ :: Matrix Double -> ExpectationQmu -> QZ -> QZ
calculateQZ matrixX (eQmuM,eQmuMi) qZ = qZ + matrixX <> tr eQmuM +
    (1 - matrixX) <>
    tr eQmuMi

calculateEQmu :: Qmu -> ExpectationQmu
calculateEQmu (aM, bM) = (reshape d qMuM, reshape d qMuMi)
  where
    aV = flatten aM
    bV = flatten bM
    qMuM = zipVectorWith (\a b -> expectationBeta a b Log) aV bV
    qMuMi = zipVectorWith (\a b -> expectationBeta a b InvLog) aV bV
    (k, d) = size aM

calculateQmu :: Matrix Double -> ExpectationQZ -> Qmu -> Qmu
calculateQmu matrixX eQz (aOld,bOld) = ( aOld - 1 + tr eQz <> matrixX
                                       , bOld - 1 +
                                         tr eQz <>
                                         (1 - matrixX))

lQTheta :: Matrix Double
        -> Vector Double
        -> Vector Double
        -> ExpectationQZ
        -> ExpectationQmu
        -> Double
lQTheta matrixX vectorPi vectorAlpha eQz eQmu = lQ matrixX eQz eQmu +
    vectorAlpha `dot`
    cmap log vectorPi

mStep :: Vector Double -> Vector Double -> ExpectationQZ -> Vector Double
mStep vectorPi vectorAlpha eQz = let vectorNewPi = fromList
                                             (map
                                                  (\k ->
                                                        -(vectorPi ! k) *
                                                         sumElements
                                                              (eqZt ! k) /
                                                         (vectorAlpha ! k - 1))
                                                  [0 .. nMixtures - 1])
    in norm vectorNewPi
    where eqZt = tr eQz
          norm v = v /
              scalar (sumElements v)

nMixtures :: Int
nMixtures = 10

max_stars_repo_name: DbIHbKA/mnistvb

max_issues_repo_name: DbIHbKA/mnistvb

max_forks_repo_name: DbIHbKA/mnistvb


text:
{-# LANGUAGE DeriveAnyClass  #-}
{-# LANGUAGE DeriveGeneric   #-}
{-# LANGUAGE RecursiveDo     #-}
{-# LANGUAGE TemplateHaskell #-}
module Main where

--import           Circuits.Network
import           Control.Lens
import           Control.Monad.Fix
import qualified Data.Vector                  as V
import qualified Data.Vector.Storable         as VS
import qualified Data.Vector.Storable.Mutable as VSM
import           GHC.Generics
import qualified Numeric.LinearAlgebra        as HMatrix
import qualified Numeric.LinearAlgebra.Data   as HMatrix
import           Text.Printf

import           Circuits.Analysis.DC
import           Circuits.Circuit
import           Circuits.Components.Linear

data TestCircuit = TestCircuit
  { _gnd :: Node
  , _n1  :: Node
  , _vs  :: VoltageSource
  , _r   :: Resistor
  }
  deriving (Show, Generic, SimulateDC, HasVariables)

c1 :: Circuit TestCircuit
c1 = buildCircuit $ do
  gnd <- ground
  n1 <- newNode
  vs <- newBranch
  return TestCircuit
    { _gnd = gnd
    , _n1  = n1
    , _vs  = VoltageSource (Independent 5) vs gnd n1
    , _r   = Resistor 100 gnd n1
    }

data RC = RC
  { _rcGnd :: Node
  , _rcN1  :: Node
  , _rcN2  :: Node
  , _rcVS  :: VoltageSource
  , _rcR   :: Resistor
  , _rcC   :: Capacitor
  }
  deriving (Show, Generic, SimulateDC, HasVariables)
makeLenses ''RC

rc :: Circuit RC
rc = buildCircuit $ do
  gnd <- ground
  n1 <- newNode
  n2 <- newNode
  vs <- newBranch
  return RC
    { _rcGnd = gnd
    , _rcN1  = n1
    , _rcN2  = n2
    , _rcVS  = VoltageSource (Independent 5) vs gnd n2
    , _rcR   = Resistor 100 gnd n1
    , _rcC   = Capacitor 0.01 n1 n2
    }

main :: IO ()
main = do
  let c2 = step SteadyState c1
  print c2
  let rcs = transient (Transient 0 4 0.1) rc
  forMOf_ (traverse . Control.Lens.to (over _2 $ view $ model . rcC . voltageAcross)) rcs $ uncurry $ printf "%.1f ~ %.2f\n"
  --mapM_ (\(t, c) -> printf "t = %.1f\n%s\n\n" t (show c)) rcs
  return ()

data Transient = Transient
  { _transientStartTime :: Double
  , _transientStopTime  :: Double
  , _transientTimeStep  :: Double
  }

transient :: (SimulateDC c, HasVariables c) => Transient -> Circuit c -> [(Double, Circuit c)]
transient (Transient start stop dt) c0 = go c0 start where
  go c t
    | t > stop = []
    | otherwise = (t, c) : go (step (Step t dt) c) (t + dt)

  {-
  sys <- newMSystem 1 1
  stamp sys (VoltageSource (Independent 5) (Var 1) Ground (Var 0))
  stamp sys (Resistor 100 Ground (Var 0) :: Resistor NodeRef NodeRef Double)
  stamp sys (Resistor 200 Ground (Var 0) :: Resistor NodeRef NodeRef Double)
  mv <- HMatrix.reshape 2 <$> VS.unsafeFreeze (_msystemCoeff sys)
  rhs <- HMatrix.asColumn <$> VS.unsafeFreeze (_msystemRhs sys)
  print mv
  print rhs
  print $ HMatrix.linearSolve mv rhs
-}
{-  let simple = empty "GND"
        & part "R1" .~ Just (Resistor 100 "N1" "GND")
        & part "C1" .~ Just (Capacitor 0.01 "N2" "N1")
        & part "Vs" .~ Just (VoltageSource 5 "GND" "N2")

  case compile simple of
    Left err -> putStrLn err
    Right c -> do
      let cs = take 100 $ iterate (step 0.01) c
      mapM_ print cs
-}

max_stars_repo_name: fatho/circuits

max_issues_repo_name: fatho/circuits

max_forks_repo_name: fatho/circuits


text:
module Assign_1_ExtraCredit where

import Data.Complex
-- see https://www.stackage.org/haddock/lts-8.24/base-4.9.1.0/Data-Complex.html

macid :: String
macid = "hutchm6"

{-
    - Assignment 1
    - Name: Mark Hutchison
    - Date: September 17th, 2019
 -}

{-
    Similar to the normal assignment, I will be using the 4 base equations that will be listed under this comment string. That being said, because the equations for Q R and Disc do not change, I will not be re showing the outputs for those above the functions.

    - cubicRealSolutions 1 (-6) 11 (-6) = should [1,2,3]
    - cubicRealSolutions 1 0 (-3) 0 = [0, sqrt 3, -sqrt 3]
    - cubicRealSolutions 1 (-3) 3 (-1) = should [1,1,1]
    - cubicRealSolutions 1 (-5) 8 (-4) = should [1,2,2]
-}
cubicQ :: Double -> Double -> Double -> Double
cubicQ a b c = ((3 * a * c) - (b ** 2)) / (9 * (a ** 2))

cubicR :: Double -> Double -> Double -> Double -> Double
cubicR a b c d = ((9 * a * b * c) - (27 * (a ** 2) * d) - (2 * (b ** 3))) / (54 * (a ** 3))

cubicDisc :: Double -> Double -> Double
cubicDisc q r = do
    if r == 0
        then q
        else ((q ** 3) + (r ** 2)) ** (1 / 3)

{-
    Here is the point where everything changes and we need to begin thinking in complex number systems.

    We need to be able to get the square root of a negative number, so the S and T function needed to change.
-}

cubicComplexS :: Double -> Double -> Complex Double
cubicComplexS q r = s where
    s = cubeRoot (r + sqrt (cubicDisc q r))

{-

-}

cubicComplexT :: Double -> Double -> Complex Double
cubicComplexT q r = t where
    t = cubeRoot (r - sqrt (cubicDisc q r))

cubicComplexSolutions :: Double -> Double -> Double -> Double -> [Complex Double]
cubicComplexSolutions a b c d = [x1, x2, x3] where
    q = cubicQ a b c
    r = cubicR a b c d
    s = cubicComplexS q r
    t = cubicComplexT q r
    x1 = (s + t - (b/(3*a)))
    x2 = ((-1) * ((s + t) / 2) - (b/(3*a)) + (((sqrt 3)/(2)) * (s - t)))
    x3 = ((-1) * ((s + t) / 2) - (b/(3*a)) - (((sqrt 3)/(2)) * (s - t)))

cubeRoot :: Double -> Double
cubeRoot n = do
    if n == 0
        then 0
        else
            if n > 0
                then n ** (1/3)
                else
                    if n < 0
                        then (-1) * ((-n) ** (1/3))
                        else n ** (1/3)

max_stars_repo_name: AstroCoding/haskell-cubic-solver

max_issues_repo_name: AstroCoding/haskell-cubic-solver

max_forks_repo_name: AstroCoding/haskell-cubic-solver