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

import           Data.Array.Repa      as R
import           Data.Array.Unboxed   as AU
import           Data.Complex
import           Data.List            as L
import           Data.Vector.Storable as VS
import           DFT.Plan
import           Filter.Utils
import           Types
import           Utils.Array
import           Utils.Parallel

data PowerMethodNormalizationOption
  = PowerMethodPatchNorm (VS.Vector (Complex Double))
  | PowerMethodGlobal
  | PowerMethodConnection [[(Int, Int)]]
  | PowerMethodConnectionMax [[(Int, Int)]]

-- {-# INLINE makePatchNormFilter #-}
-- makePatchNormFilter ::
--      DFTPlan
--   -> Int
--   -> Int
--   -> Bool
--   -> Int
--   -> IO (DFTPlan, PowerMethodNormalizationOption)
-- makePatchNormFilter plan cols rows True n = do
--   let arr =
--         makeFilter2D . fromFunction (Z :. cols :. rows) $ \(Z :. i :. j) ->
--           if (sqrt . fromIntegral $ (i - div cols 2) ^ 2 + (j - div rows 2) ^ 2) <=
--              fromIntegral n
--             then 1
--             else 0
--   lock <- getFFTWLock
--   (p, vec) <-
--     dft1dGPlan lock plan [cols, rows] [0, 1] . VS.convert . toUnboxed . computeS $
--     arr
--   (newPlan, _) <- idft1dGPlan lock p [cols, rows] [0, 1] $ vec
--   return (newPlan, PowerMethodPatchNorm vec)
-- makePatchNormFilter plan _ _ False _ = return (plan, PowerMethodGlobal)

-- {-# INLINE patchNorm #-}
-- patchNorm ::
--      (R.Source s (Complex Double))
--   => DFTPlan
--   -> VS.Vector (Complex Double)
--   -> R.Array s DIM4 (Complex Double)
--   -> IO (R.Array U DIM4 (Complex Double))
-- patchNorm plan filterF arr = do
--   let (Z :. numThetaFreq :. numScaleFreq :. cols :. rows) = extent arr
--   arrSum <- (R.sumP . rotate4D . rotate4D . R.map magnitude $ arr) >>= R.sumP
--   arrSumF <-
--     dftExecute plan (DFTPlanID DFT1DG [cols, rows] [0, 1]) .
--     VS.convert . toUnboxed . computeS . R.map (:+ 0) $
--     arrSum
--   arrPatchSum <-
--     fmap (fromUnboxed (Z :. cols :. rows) . VS.convert) .
--     dftExecute plan (DFTPlanID IDFT1DG [cols, rows] [0, 1]) $
--     VS.zipWith (*) arrSumF filterF
--   computeP . R.traverse2 arr arrPatchSum const $ \f fNorm idx@(Z :. _ :. _ :. i :. j) ->
--     if fNorm (Z :. i :. j) == 0
--       then 0
--       else f idx / fNorm (Z :. i :. j)

-- {-# INLINE powerMethodNormalization #-}
-- powerMethodNormalization ::
--      (R.Source s (Complex Double))
--   => DFTPlan
--   -> PowerMethodNormalizationOption
--   -> R.Array s DIM4 (Complex Double)
--   -> IO (R.Array U DIM4 (Complex Double))
-- powerMethodNormalization plan (PowerMethodPatchNorm filterF) arr =
--   patchNorm plan filterF arr
-- powerMethodNormalization _ PowerMethodGlobal arr = do
--   s <- R.sumAllP . rotate4D . rotate4D . R.map magnitude $ arr
--   return . computeS . R.map (/ (s :+ 0)) $ arr

{-# INLINE powerMethodNormalization #-}
powerMethodNormalization ::
     (R.Source s (Complex Double))
  => PowerMethodNormalizationOption
  -> R.Array s DIM4 (Complex Double)
  -> IO (R.Array U DIM4 (Complex Double))
powerMethodNormalization PowerMethodGlobal arr = do
  -- s <- R.sumAllP . rotate4D . rotate4D . R.map (\x -> (magnitude x) ^ 2) $ arr
  s <- R.foldAllP max 0 . rotate4D . rotate4D . R.map (\x -> magnitude x) $ arr
  return . computeS . R.map (/ (s :+ 0)) $ arr
powerMethodNormalization (PowerMethodConnection xss) arr = do
  let (Z :. _ :. _ :. cols :. rows) = extent arr
      arrNorm =
        fromListUnboxed (Z :. cols :. rows) . AU.elems $
        (AU.accumArray (+) 0 ((0, 0), (cols - 1, rows - 1)) .
         L.concat .
         parMap
           rdeepseq
           (\xs ->
              let s =
                    L.maximum .
                    L.map
                      (\(i, j) ->
                         R.foldAllS max 0 .
                         R.slice (R.map (\x -> magnitude x) arr) $
                         (Z :. All :. All :. i :. j)) $
                    xs
               in L.map
                    (\x ->
                       ( x
                       , if s == 0
                           then 0
                           else 1 / s))
                    xs) $
         xss :: UArray (Int, Int) Double)
  return . computeS . R.traverse2 arr arrNorm const $ \f fNorm idx@(Z :. _ :. _ :. i :. j) ->
    f idx * ((fNorm (Z :. i :. j)) :+ 0)
powerMethodNormalization (PowerMethodConnectionMax xss) arr = do
  let (Z :. numThetaFreq :. numScaleFreq :. cols :. rows) = extent arr
      arrMax =
        fromListUnboxed (Z :. cols :. rows) . AU.elems $
        (AU.accumArray (+) 0 ((0, 0), (cols - 1, rows - 1)) .
         L.concat .
         parMap
           rdeepseq
           (\xs ->
              let s =
                    sqrt .
                    L.maximum .
                    L.map
                      (\(i, j) ->
                         R.sumAllS .
                         R.slice (R.map (\x -> (magnitude x) ^ 2) arr) $
                         (Z :. All :. All :. i :. j)) $
                    xs
               in L.map (\x -> (x, s)) xs) $
         xss :: UArray (Int, Int) Double)
      arrSum =
        fromListUnboxed (Z :. cols :. rows) . AU.elems $
        (AU.accumArray (+) 0 ((0, 0), (cols - 1, rows - 1)) .
         L.concat .
         parMap
           rdeepseq
           (\xs ->
              L.map
                (\(i, j) ->
                   ( (i, j)
                   , sqrt .
                     R.sumAllS . R.slice (R.map (\x -> (magnitude x) ^ 2) arr) $
                     (Z :. All :. All :. i :. j))) $
              xs) $
         xss :: UArray (Int, Int) Double)
  return . computeS . R.traverse3 arr arrMax arrSum (\s1 s2 s3 -> s1) $ \f fMax fSum idx@(Z :. k :. l :. i :. j) ->
    if fSum (Z :. i :. j) > 0.9 * fMax (Z :. i :. j) && k == div numThetaFreq 2 && l == div numScaleFreq 2
      then 1 -- (f idx) / (fSum (Z :. i :. j) :+ 0)
      else 0

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 ScopedTypeVariables #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE RankNTypes #-}

module Hw7 where

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

import Control.Monad.State
import Control.Monad (forM_)
import System.Directory (doesDirectoryExist, listDirectory, Permissions(..),
                         getModificationTime, getPermissions, getFileSize, getDirectoryContents)
import System.FilePath ((</>), takeExtension)
import Control.Exception (handle, SomeException(..))

import Data.Time (UTCTime(..))
import Data.Time.Calendar(toGregorian)
import Data.Time.Clock(utctDay,UTCTime)
import Data.Time.Clock.POSIX(posixSecondsToUTCTime)

epochToUTC :: Integral a => a -> UTCTime
epochToUTC = posixSecondsToUTCTime . fromIntegral

epoch = epochToUTC 0

match :: String -> String -> Bool
match s p
  | length s < l = False
  | otherwise = (take l s) == p || match (drop 1 s) p
  where l = length p

class Applicative m => Predicate m where
   (>?) :: (Ord a) => m a -> a -> m Bool
   p >? a = (> a) <$> p
   infixl 4 >?

   (<?) :: (Ord a) => m a -> a -> m Bool
   p <? a = (< a) <$> p
   infixl 4 <?

   (==?) :: (Ord a) => m a -> a -> m Bool
   p ==? a = (== a) <$> p
   infixl 4 ==?

   (!=?) :: (Ord a) => m a -> a -> m Bool
   p !=? a = (/= a) <$> p
   infixl 4 !=?

   (~?) :: m String -> String -> m Bool
   p ~? s = (`match` s) <$> p
   infixl 4 ~?

   (&&?) :: m Bool -> m Bool -> m Bool
   p1 &&? p2 = (&&) <$> p1 <*> p2
   infixl 3 &&?

   (||?) :: m Bool -> m Bool -> m Bool
   p1 ||? p2 = (||) <$> p1 <*> p2
   infixl 3 ||?

   notP :: m Bool -> m Bool
   notP = fmap (not)

data FileInfo = FileInfo {
                    filePath :: FilePath,
                    filePermissions :: Maybe Permissions,
                    fileSize :: Maybe Integer,
                    fileLastModified :: Maybe UTCTime
                }

newtype FilePredicate a = FileP { runFileP :: FileInfo -> a } deriving (Functor,Applicative)

instance Predicate FilePredicate

pathP :: FilePredicate String
pathP = FileP filePath

extP :: FilePredicate String
extP = takeExtension <$> pathP

sizeP :: FilePredicate Integer
sizeP = maybe (-1) id <$> FileP fileSize

getPermissionP :: (Permissions -> Bool) -> FilePredicate Bool
getPermissionP f = maybe False f <$> FileP filePermissions

searchP:: FilePredicate Bool
searchP = getPermissionP searchable

readP :: FilePredicate Bool
readP = getPermissionP readable

writeP :: FilePredicate Bool
writeP = getPermissionP writable

execP :: FilePredicate Bool
execP = getPermissionP executable

timeP :: FilePredicate UTCTime
timeP = maybe epoch id <$> FileP fileLastModified

maybeIO :: IO a -> IO (Maybe a)
maybeIO act = handle (\(SomeException _) -> return Nothing) (fmap Just act)

getInfo :: FilePath -> IO FileInfo
getInfo path = do
  perms <- maybeIO (getPermissions path)
  size <- maybeIO (getFileSize path)
  modified <- maybeIO (getModificationTime path)
  return (FileInfo path perms size modified)

data FileAction = Done | Skip | Continue

type FileState s =  StateT s IO FileAction

find :: forall s. (FileInfo -> FileState s) -> FilePath -> FileState s

find getState path = do
  names <- lift $ listDirectory path
  iterate names
  where iterate :: [FilePath] -> FileState s
        iterate (name: names) = do
            let path' = path </> name
            let info = getInfo path'
            state <- getState =<< lift info
            isDirectory <- lift $ doesDirectoryExist path'

            case state of
                Done -> return state
                Continue
                    | isDirectory -> do
                        next <- find getState path'
                        case next of
                            Done -> return next
                            _ -> iterate names
                    | otherwise -> iterate names
                Skip -> iterate names

        iterate _ = return Continue

main :: IO ()

main = do
--    let downloads = "C:\\Users\\tzhao\\Downloads"
   let downloads = "/Users/nikolapeevski/Downloads"

   let yearP = ((\(x,_,_) -> x) . toGregorian . utctDay) <$> timeP
   let recurseP = yearP >? 2018

   let filterP = (extP ==? ".pdf" &&? sizeP >? 2^22) ||? execP

   let iter :: FileInfo -> FileState [(FilePath, Integer, UTCTime)]
       iter = \info -> do
            let eval = \p -> runFileP p info

            s <- get
            let r = (eval pathP, eval sizeP, eval timeP)
            let s' = if eval filterP then r : s else s
            put s'

            let a = if length s' >= 10 then Done
                    else if eval recurseP then Continue else Skip
            return a

   results <- execStateT (find iter downloads) []

   forM_ results $ \(p,n,t) ->
                            handle (\(SomeException _) -> return ())
                          $ putStrLn
                          $ (take 100 p) ++ "\t"
                              ++ show (round $ fromInteger n / 2^10) ++ " KB\t"
                              ++ show (utctDay t)

ext: hs

max_stars_repo_name: NikolaPeevski/Haskell-stuff

max_issues_repo_name: NikolaPeevski/Haskell-stuff

max_forks_repo_name: NikolaPeevski/Haskell-stuff


text:
module Statistics.Classification.ConfusionMatrix.Multi (
    MultiConfusionMatrix(..)
  , multiClassCount
  , multiTotalCount
  , truePositiveCount, falsePositiveCount, falseNegativeCount, trueNegativeCount
  , confusionMatrix
  , PosNeg
  , posVal, negVal, posVals, negVals
  , trueVal, falseVal
  , isPos, isNeg
  ) where

import qualified Data.Map.Strict as Map
import qualified Data.Set as Set

import Statistics.Classification.Types (ClassificationResult(..))

data MultiConfusionMatrix a = MultiConfusionMatrix {
    confCounts :: Map.Map (ClassificationResult a) Int
  , confClasses :: Set.Set a
  } deriving (Eq, Show, Read)

instance Ord a => Semigroup (MultiConfusionMatrix a) where
  (MultiConfusionMatrix mcm0 mcc0) <> (MultiConfusionMatrix mcm1 mcc1) =
    MultiConfusionMatrix (Map.unionWith (+) mcm0 mcm1) (Set.union mcc0 mcc1)

instance Ord a => Monoid (MultiConfusionMatrix a) where
  mempty = MultiConfusionMatrix Map.empty Set.empty
  mappend = (<>)

multiClassCount :: Ord a => MultiConfusionMatrix a -> ClassificationResult a -> Int
multiClassCount (MultiConfusionMatrix conf _) r = Map.findWithDefault 0 r conf

truePositiveCount,falsePositiveCount,falseNegativeCount,trueNegativeCount :: Ord a => MultiConfusionMatrix a -> PosNeg a -> Int
truePositiveCount (MultiConfusionMatrix conf _) posneg = sum $
  Map.filterWithKey (\(ClassificationResult a p) _-> isPos posneg a && isPos posneg p) conf
falsePositiveCount (MultiConfusionMatrix conf _) posneg = sum $
  Map.filterWithKey (\(ClassificationResult a p) _-> isNeg posneg a && isPos posneg p) conf
falseNegativeCount (MultiConfusionMatrix conf _) posneg = sum $
  Map.filterWithKey (\(ClassificationResult a p) _-> isPos posneg a && isNeg posneg p) conf
trueNegativeCount (MultiConfusionMatrix conf _) posneg = sum $
  Map.filterWithKey (\(ClassificationResult a p) _-> isNeg posneg a && isNeg posneg p) conf

multiTotalCount :: MultiConfusionMatrix a -> Int
multiTotalCount = sum . confCounts

confusionMatrix :: Ord a => [ClassificationResult a] -> MultiConfusionMatrix a
confusionMatrix = foldMap go
  where
    go k = MultiConfusionMatrix {
        confCounts = Map.singleton k 1
      , confClasses = Set.fromList [classifiedActual k, classifiedPredicted k]
      }

data PosNeg a = PosVal (Set.Set a) | NegVal (Set.Set a)
  deriving (Eq, Show, Read)

posVal,negVal :: a -> PosNeg a
posVal = PosVal . Set.singleton
negVal = NegVal . Set.singleton

posVals,negVals :: Ord a => [a] -> PosNeg a
posVals = PosVal . Set.fromList
negVals = NegVal . Set.fromList

trueVal,falseVal :: PosNeg Bool
trueVal = posVal True
falseVal = posVal False

isPos,isNeg :: Ord a => PosNeg a -> a -> Bool
isPos (PosVal pos) v = Set.member v pos
isPos (NegVal neg) v = Set.notMember v neg
isNeg (PosVal pos) v = Set.notMember v pos
isNeg (NegVal neg) v = Set.member v neg

max_stars_repo_name: tsbattman/rochs

max_issues_repo_name: tsbattman/rochs

max_forks_repo_name: tsbattman/rochs


text:
{-# LANGUAGE NoMonomorphismRestriction #-}
{-# LANGUAGE PartialTypeSignatures #-}
{-# LANGUAGE ScopedTypeVariables #-}

{-# OPTIONS_GHC -fno-warn-partial-type-signatures #-}

import qualified Prelude

import qualified Data.Complex as Complex

import Feldspar.Run
import Feldspar.Data.Vector
import Feldspar.Data.Buffered

import qualified Test.QuickCheck as QC
import qualified Test.QuickCheck.Monadic as QC

import Test.Tasty
import Test.Tasty.HUnit
import Test.Tasty.QuickCheck

import qualified Tut1_HelloWorld            as Tut1
import qualified Tut2_ExpressionsAndTypes   as Tut2
import qualified Tut3_Vectors               as Tut3
import qualified Tut4_MemoryManagement      as Tut4
import qualified Tut5_Matrices              as Tut5
import qualified Tut6_Testing               as Tut6
import qualified Tut7_ImperativeProgramming as Tut7
import qualified Tut8_SequentialVectors     as Tut8
import qualified Concurrent
import DFT
import FFT

almostEq a b
    =          Complex.magnitude d Prelude.< 1e-7
    Prelude.&& Complex.phase d     Prelude.< 1e-7
  where
    d = abs (a-b)

a ~= b = Prelude.and $ Prelude.zipWith almostEq a b

wrapStore :: (Syntax a, Finite (vec a), MonadComp m) =>
    (Store a -> vec a -> m b) -> vec a -> m b
wrapStore f v = do
    st <- newStore $ length v
    f st v

fftS u = wrapStore (flip fft u)  :: DManifest (Complex Double) -> _
ifftS  = wrapStore (flip ifft 1) :: DManifest (Complex Double) -> _

prop_fft_dft dft' fft' = QC.monadicIO $ do
    n   :: Int              <- QC.pick $ QC.choose (2,5)
    inp :: [Complex Double] <- QC.pick $ QC.vector (2 Prelude.^ n)
    outd <- QC.run $ dft' inp
    outf <- QC.run $ fft' inp
    QC.assert (outd ~= outf)

prop_inverse f fi = QC.monadicIO $ do
    n   :: Int              <- QC.pick $ QC.choose (2,5)
    inp :: [Complex Double] <- QC.pick $ QC.vector (2 Prelude.^ n)
    out1 <- QC.run $ f inp
    out2 <- QC.run $ fi out1
    QC.assert (inp ~= out2)

prop_fib fb1 fb2 = QC.monadicIO $ do
    n   <- QC.pick $ QC.choose (0,40)
    fs1 <- QC.run $ fb1 n
    fs2 <- QC.run $ fb2 n
    QC.assert (fs1 Prelude.== fs2)

main =
    marshalledM (return . dft)  $ \dft'  ->
    marshalledM (return . idft) $ \idft' ->
    marshalledM (fftS 1)        $ \fft1  ->
    marshalledM (fftS 2)        $ \fft2  ->
    marshalledM ifftS           $ \ifft' ->

    marshalled (return . Tut8.fibSeq) $ \fb1 ->
    marshalled (\n -> return $ Pull n Tut2.fib) $ \fb2 ->

      defaultMain $ testGroup "tests"
        [ testCase "Tut1"       Tut1.testAll
        , testCase "Tut2"       Tut2.testAll
        , testCase "Tut3"       Tut3.testAll
        , testCase "Tut4"       Tut4.testAll
        , testCase "Tut5"       Tut5.testAll
        , testCase "Tut6"       Tut6.testAll
        , testCase "Tut7"       Tut7.testAll
        , testCase "Tut8"       Tut8.testAll
        , testCase "Concurrent" Concurrent.testAll
        , testProperty "fft1_dft" $ prop_fft_dft dft' fft1
        , testProperty "fft2_dft" $ prop_fft_dft dft' fft2
        , testProperty "dft_idft" $ prop_inverse dft' idft'
        , testProperty "fft_ifft" $ prop_inverse fft1 ifft'
        , testProperty "fib"      $ prop_fib fb1 fb2
        ]

  where
    marshalledM = marshalled' def def {externalFlagsPost = ["-lm"]}

max_stars_repo_name: Abhiroop/mu-feldspar

max_issues_repo_name: Abhiroop/mu-feldspar

max_forks_repo_name: Abhiroop/mu-feldspar


text:
--Quantum circuits
--Copyright Laurence Emms 2018

module QuantumCircuit (bell) where

import Data.Complex
import Register
import Instruction

bell :: (Floating a, Show a, Ord a) => Circuit a
bell = let parameters@[Left a, Left b] = [Left (MetaQubitRegister "a"), Left (MetaQubitRegister "b")]
           instructions = [(Hadamard a),
                           (CNot a b)]
       in Circuit "BELL" parameters instructions

clear :: (Floating a, Show a, Ord a) => Circuit a
clear = let parameters@[Left q, Right scratch] = [Left (MetaQubitRegister "q"), Right (MetaRegister "scratch")]
            instructions = [(MeasureOut q scratch),
                            (JumpUnless "CLEAREND" scratch),
                            (PauliX q),
                            (Label "CLEAREND")]
        in Circuit "CLEAR" parameters instructions

euler :: (Floating a, Show a, Ord a) => Complex a -> Complex a -> Complex a -> Circuit a
euler alpha beta gamma = let parameters@[Left q] = [Left (MetaQubitRegister "q")]
                             instructions = [(RX (ComplexConstant alpha) q),
                                             (RY (ComplexConstant alpha) q),
                                             (RZ (ComplexConstant alpha) q)]
                         in Circuit "EULER" parameters instructions

max_stars_repo_name: WhatTheFunctional/Hasquil

max_issues_repo_name: WhatTheFunctional/Hasquil

max_forks_repo_name: WhatTheFunctional/Hasquil


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE RankNTypes #-}
module Statistics.BBVI.Examples
  ( genMixture
  , mixtureFit
  , genNormal
  , normalFit
  , genDirichlet
  , dirichletFit
  , genMixedMem
  , mixedMemFit
  )
where

import           Control.Monad.ST
import           Data.Propagator
import qualified Data.Vector                   as V
import           Numeric.AD                     ( grad'
                                                , auto
                                                , diff
                                                )
import           System.Random.MWC              ( create
                                                , uniform
                                                , initialize
                                                , GenST
                                                )
import qualified System.Random.MWC.Distributions
                                               as MWCD
import           Statistics.BBVI

--- helper functions
logSum :: (Floating c, Ord c) => V.Vector c -> V.Vector c -> c
logSum v1 = log . minimumProb . V.sum . V.map exp . V.zipWith (+) v1

minimumProb :: (Floating a, Ord a) => a -> a
minimumProb = max 1e-308

--- useful global params
-- for fit
globalMaxStep :: Int
globalMaxStep = 1000
globalDelta :: Double
globalDelta = 1e-16 -- 0.00001 --

-- for generating data
nDimension :: Int
nDimension = 1
nTime :: Int
nTime = 1000
nState :: Int
nState = 3 -- 10

-----------------
--  mixture
------------------
-- | log joint of a mixture model, with categorical marginalized out
mixtureJoint :: (Floating a, Ord a) => a -> V.Vector a -> a -> V.Vector a -> a
mixtureJoint std theta ob betas = logSum
  (V.map log theta)
  (V.map (\mu -> diffableNormalLogProb mu std ob) betas)

-- | transform a log joint to a gradient propagator, for use with
-- models with a mixture-like plate structure. Updates the mixture
-- components with the reparameterization gradient (Kucukelbir et al
-- 2017).
updateReparam
  :: (Dist b SampleVector, Differentiable c Double) -- Dist a SampleVector,
  => Int
  -> (  forall e
      . (Floating e, Ord e)
     => V.Vector e
     -> e
     -> V.Vector e
     -> e
     )
  -> GenST s
  -> V.Vector (V.Vector Double)
  -> DistInvariant b
  -> DistInvariant c
  -> DistCell b
  -> V.Vector (V.Vector (DistCell c))
  -> ST
       s
       (DistCell b, V.Vector (V.Vector (DistCell c)))
updateReparam nSamp jointF gen obss thetaG betaG thetaN betasN = do
  -- obss        <- V.replicateM nObs (resample xss gen)
  thetaSamp   <- V.replicateM nSamp (resample theta gen)
  betaSamples <- V.replicateM nSamp (epsilon ((betass V.! 0) V.! 0) gen)
  let
    (joints, gradJoints) = V.unzip $ V.zipWith
      (\e th ->
        V.foldl1'
            (\(x1, y1) (x2, y2) -> (x1 + x2, V.zipWith (V.zipWith (+)) y1 y2))
          $ V.map
              (\obs ->
                let (joints', gradss) = V.unzip $ V.zipWith
                      (\ob betas -> grad' (jointF (V.map auto th) (auto ob))
                                          (V.map (\d -> transform d e) betas)
                      )
                      obs
                      betass
                in  (V.sum joints', gradss)
              )
              (obss :: V.Vector (V.Vector Double))
      )
      betaSamples
      thetaSamp
  return
    $ ( gradientScore thetaG thetaN (nFactor, joints, thetaSamp)
      , V.imap
        (\dim v -> V.imap
          (\c d -> gradientReparam
            betaG
            d
            (nFactor, V.map ((V.! c) . (V.! dim)) gradJoints, betaSamples)
          )
          v
        )
        betasN
      )
 where
  nFactor = fromIntegral $ V.length obss
  theta   = dist thetaN
  betass  = V.map (V.map dist) betasN

-- | transform a log joint to a gradient propagator, for use with
-- models with a mixture-like plate structure. Updates the mixture
-- components with the score gradient (Ranganath et al 2014)
updateScore
  :: (Dist a1 c1, Dist a2 c2)
  => Int
  -> (c1 -> a3 -> V.Vector c2 -> Double)
  -> GenST s
  -> V.Vector (V.Vector a3)
  -> DistInvariant a1
  -> DistInvariant a2
  -> DistCell a1
  -> V.Vector (V.Vector (DistCell a2))
  -> ST
       s
       ( DistCell a1
       , V.Vector (V.Vector (DistCell a2))
       )
updateScore nSamp jointF gen obss thetaG betaG thetaN betasN = do
  -- obss        <- V.replicateM nObs (resample xs gen)
  thetaSamp   <- V.replicateM nSamp (resample theta gen)
  betaSamples <- V.replicateM nSamp
                              (V.mapM (V.mapM (\b -> resample b gen)) betass)
  let joints = V.zipWith
        (\theta' betass' ->
          V.foldl1' (V.zipWith (+))
            . V.map (\obs -> V.zipWith (jointF theta') obs betass')
            $ obss
        )
        thetaSamp
        betaSamples
  return
    $ ( gradientScore thetaG thetaN (nFactor, V.map V.sum joints, thetaSamp)
      , V.imap
        (\dim v -> V.imap
          (\c d -> gradientScore
            betaG
            d
            ( nFactor
            , V.map (V.! dim) joints
            , V.map ((V.! c) . (V.! dim)) betaSamples
            )
          )
          v
        )
        betasN
      )
 where
  -- xs     = dist xsN
  nFactor = fromIntegral $ V.length obss
  theta   = dist thetaN
  betass  = V.map (V.map dist) betasN

-- | generate data for testing mixture model
genMixture :: ST s (V.Vector (V.Vector Double))
genMixture = do
  gen <- create
  let theta' = MWCD.categorical (V.replicate nState 1.0) gen
  let std    = 1.0
  let
    mixtures = V.generate
      nState
      (\k ->
        (MWCD.normal
          ((fromIntegral k - (fromIntegral (nState - 1) * 0.5)) * 5)
          std
          gen
        )
      )
  xs <- V.replicateM nTime (V.replicateM nDimension . (mixtures V.!) =<< theta')
  return xs

-- | fit a mixture model, given data
mixtureFit
  :: V.Vector (V.Vector Double) -> (Dirichlet, V.Vector (V.Vector NormalDist))
mixtureFit xs = runST $ do
  genG <- create
  gen1 <- initialize =<< V.replicateM 256 (uniform genG)
  let priorTheta = dirichlet (V.replicate nState 1.0)
  let priorBeta = normalDistr 0 (5 * (fromIntegral nState - 1) :: Double)
  let nSamp      = 10
  let localStep  = (20 :: Int)

  -- initialize distribution cells (both invariant and variants)
  let thetaGrad = DistInvariant (fromIntegral $ V.length xs)
                                priorTheta
                                (rhoKuc defaultKucP) --
  qTheta <- cellWith $ mergeGeneric globalMaxStep globalDelta
  write qTheta $ defaultDistCell (dirichlet (V.replicate nState 1.0))
  let betaGrad = DistInvariant (fromIntegral $ V.length xs)
                               priorBeta
                               (rhoKuc defaultKucP) --
  qBetas <- cellWith $ mergeGenericss globalMaxStep globalDelta
  write qBetas =<< V.replicateM
    nDimension
    (V.generateM
      nState
      (\_i -> do
        mu <- resample priorBeta gen1
        return $ defaultDistCell (normalDistr mu (1.0 :: Double))
      )
    )

  -- attach gradient propagators to cells
  stepTogether
    (updateReparam nSamp (mixtureJoint 1.0) gen1 xs thetaGrad betaGrad)
    qTheta
    qBetas
  -- stepSeparate localStep
  --              globalDelta
  --              (updateScore nSamp (mixtureJoint 1.0) gen1 xs thetaGrad betaGrad)
  --              qTheta
  --              qBetas

  -- pull content (run network to quiesence)
  thetaF <- unsafeContent qTheta
  betaF  <- unsafeContent qBetas
  let betaDists = V.map (V.map dist) betaF
  return (dist thetaF, betaDists)

--------------------------------------
-- fit simple normal distribution mean
--------------------------------------
-- | generate normally distributed data for testing
genNormal :: ST s (V.Vector Double)
genNormal = do
  gen <- create
  xs  <- V.replicateM 1000 (resample (normalDistr (5.0 :: Double) 3.0) gen)
  return xs

-- | log joint of normal distribution
jointNormal :: Floating a => a -> a -> a -> a
jointNormal std ob mu = diffableNormalLogProb mu std ob

-- | fit a gaussian using reparam gradient
normalFit :: Dist (Obs a) SampleDouble => V.Vector a -> (NormalDist, Time)
normalFit xs = runST $ do
  genG <- create
  gen1 <- initialize =<< V.replicateM 256 (uniform genG)
  let nSamp = 100
  xProp <- known $ defaultObs xs
  q     <- cellWith $ mergeGeneric globalMaxStep globalDelta
  write q (defaultDistCell $ normalDistr 0.0 2.0)
  let qInvar = DistInvariant (fromIntegral $ V.length xs)
                             (normalDistr 0.0 2.0)
                             (rhoKuc defaultKucP)

  (\qP xP -> watch qP $ \q' -> with xP $ \xs' ->
      singleUpdateReparam nSamp (jointNormal 1.0) gen1 qInvar xs' q' >>= write q
    )
    q
    xProp
  q' <- unsafeContent q
  return (dist q', time q')

-- | transform a log joint to a gradient propagator for a single
-- distribution cell, using reparameterization graidient. TODO: more
-- work on generalizing these types of functions!
singleUpdateReparam
  :: (Dist b SampleDouble, Differentiable c Double)
  => Int
  -> (forall  e . (Floating e, Ord e) => e -> e -> e)
  -> GenST s
  -> DistInvariant c
  -> DistCell b
  -> DistCell c
  -> ST s (DistCell c)
singleUpdateReparam nSamp joint gen qG xsN qN = do
  obs     <- V.replicateM nSamp (resample xs gen)
  samples <- V.replicateM nSamp (epsilon q gen)
  let gradJoint = V.map
        (\mu -> V.sum $ V.map (\ob -> diff (joint (auto ob)) mu) obs)
        (V.map (transform q) samples)
  return $ gradientReparam qG qN ((fromIntegral nSamp), gradJoint, samples)
 where
  xs = dist xsN
  q  = dist qN

---------------------------
-- simple dirichlet example
---------------------------
-- | generate data for testing the fit to a dirichlet
genDirichlet :: ST s (V.Vector Int)
genDirichlet = do
  gen <- create
  xs  <- V.replicateM
    1000
    (   resample (dirichlet (V.fromList [10.0, 20.0])) gen
    >>= \cat -> MWCD.categorical (cat :: V.Vector Double) gen
    )
  return (xs :: V.Vector Int)

-- | git a dirichlet to samples form a categorical
dirichletFit :: V.Vector Int -> (Dirichlet, Time)
dirichletFit xs = runST $ do
  genG <- create
  gen1 <- initialize =<< V.replicateM 256 (uniform genG)
  let priorTheta = dirichlet (V.fromList [1.0, 1.0])
  let nSamp      = 100
  qTheta <- cellWith $ mergeGeneric globalMaxStep globalDelta
  let qInvar = DistInvariant 1 priorTheta (rhoKuc defaultKucP)
  write qTheta (defaultDistCell priorTheta)
  (\tP -> watch tP $ \theta' -> do
      upTh <- dirichletGradProp nSamp gen1 xs qInvar theta'
      write tP upTh
    )
    qTheta
  thetaF <- unsafeContent qTheta
  return (dist thetaF, time thetaF)

dirichletGradProp
  :: Dist a (V.Vector Double)
  => Int
  -> GenST s
  -> V.Vector Int
  -> DistInvariant a
  -> DistCell a
  -> ST s (DistCell a)
dirichletGradProp nSamp gen xs qInvar diriQ =
  V.replicateM nSamp (resample diri gen) >>= \samples -> return $ gradientScore
    qInvar
    diriQ
    (1, V.map (\z -> V.sum $ V.map (\i -> log (z V.! i)) xs) samples, samples)
  where diri = dist diriQ

----------------------
-- mixed membership
----------------------

-- | transform a log joint to a gradient propagator, for use with
-- models with a mixed-membership-like plate structure. Updates the
-- mixture components with the reparameterization gradient (Kucukelbir
-- et al 2017). TODO : generalize all these joint transformer
-- functions better around the idea of plates
updateMembershipReparam
  :: (Dist b SampleVector, Differentiable c Double) -- Dist a SampleVector,
  => Int
  -> (  forall e
      . (Floating e, Ord e)
     => V.Vector e
     -> e
     -> V.Vector e
     -> e
     )
  -> GenST s
  -> V.Vector (V.Vector Double)
  -> DistInvariant b
  -> DistInvariant c
  -> V.Vector (DistCell b)
  -> V.Vector (V.Vector (DistCell c))
  -> ST
       s
       ( V.Vector (DistCell b)
       , V.Vector (V.Vector (DistCell c))
       )
updateMembershipReparam nSamp jointF gen obss thetaG betaG thetaN betasN = do
  thetaSamples <- V.replicateM nSamp (V.mapM (\th -> resample th gen) thetas)
  betaSamples  <- V.replicateM nSamp (epsilon ((betass V.! 0) V.! 0) gen)
  let (joints, gradJoints) = V.unzip $ V.zipWith
        (\e ths -> V.unzip $ V.zipWith
          (\obs th ->
            let (joints', gradss) = V.unzip $ V.zipWith
                  (\ob betas -> grad' (jointF (V.map auto th) (auto ob))
                                      (V.map (\d -> transform d e) betas)
                  )
                  obs
                  betass
            in  (V.sum joints', gradss)
          )
          (obss :: V.Vector (V.Vector Double))
          ths
        )
        betaSamples
        thetaSamples
  return
    $ ( V.imap
        (\i thN -> gradientScore
          thetaG
          thN
          (nFactor, V.map (V.! i) joints, (V.map (V.! i) thetaSamples))
        )
        thetaN
      , V.imap
        (\dim v -> V.imap
          (\c d -> gradientReparam
            betaG
            d
            ( nFactor
            , V.map (V.sum . V.map ((V.! c) . (V.! dim))) gradJoints
            , betaSamples
            )
          )
          v
        )
        betasN
      )
 where
  nFactor = fromIntegral $ V.length obss
  thetas  = V.map dist thetaN
  betass  = V.map (V.map dist) betasN

-- | transform a log joint to a gradient propagator, for use with
-- models with a mixed membership-like plate structure. Updates the
-- mixture components with the score gradient (Ranganath et al
-- 2014). TODO: generalize these log joint to gradient propagator
-- transformer functions
updateMembershipScore
  :: (Dist a1 c1, Dist a2 c2)
  => Int
  -> (c1 -> a3 -> V.Vector c2 -> Double)
  -> GenST s
  -> V.Vector (V.Vector a3)
  -> DistInvariant a1
  -> DistInvariant a2
  -> V.Vector (DistCell a1)
  -> V.Vector (V.Vector (DistCell a2))
  -> ST
       s
       ( V.Vector (DistCell a1)
       , V.Vector (V.Vector (DistCell a2))
       )
updateMembershipScore nSamp jointF gen obss thetaG betaG thetasN betasN = do
  -- obss        <- V.replicateM nObs (resample xs gen)
  thetaSamples <- V.replicateM nSamp
                               (V.mapM (\theta -> resample theta gen) thetas)
  betaSamples <- V.replicateM nSamp
                              (V.mapM (V.mapM (\b -> resample b gen)) betass)
  let joints = V.zipWith
        (\thetas' betass' ->
          V.zipWith (\obs theta -> V.zipWith (jointF theta) obs betass') obss
            $ thetas'
        )
        thetaSamples
        betaSamples
  return
    $ ( V.imap
        (\i thetaN -> gradientScore
          thetaG
          thetaN
          (nFactor, V.map (V.sum . (V.! i)) joints, V.map (V.! i) thetaSamples)
        )
        thetasN
      , V.imap
        (\dim v -> V.imap
          (\c d -> gradientScore
            betaG
            d
            ( nFactor
            , V.map (V.sum . V.map (V.! dim)) joints
            , V.map ((V.! c) . (V.! dim)) betaSamples
            )
          )
          v
        )
        betasN
      )
 where
  -- xs     = dist xsN
  nFactor = fromIntegral $ V.length obss
  thetas  = V.map dist thetasN
  betass  = V.map (V.map dist) betasN

-- | generate data for testing a Gaussian mixed membership model
genMixedMem :: ST s (V.Vector (V.Vector (Double))) -- Maybe Double
genMixedMem = do
  gen <- create
  -- initialize dirichlet that favors sparse categoricals
  let diriParam  = 0.01
  let thetaParam = MWCD.dirichlet (V.replicate nState diriParam) gen
  thetasTrue <- V.replicateM nTime thetaParam
  let thetas' = V.map (\theta' -> MWCD.categorical theta' gen) thetasTrue
  let betaStd = 0.001
  let betas' = V.generate
        nState
        (\k -> V.replicate
          nDimension
          (MWCD.normal
            ((fromIntegral k - (fromIntegral (nState - 1) * 0.5)) * 5)
            betaStd
            gen
          )
        )
  xs <- V.generateM
    nTime
    (\day -> V.generateM
      nDimension
      (\loc -> (thetas' V.! day) >>= \z -> ((betas' V.! z) V.! loc)) -- Just <$>
    )
  return xs

-- | fit a Gaussian mixed membership model to data
mixedMemFit
  :: V.Vector (V.Vector Double)
  -> (V.Vector Dirichlet, V.Vector (V.Vector NormalDist))
mixedMemFit xs = runST $ do
  genG <- create
  gen1 <- initialize =<< V.replicateM 256 (uniform genG)
  let priorTheta = dirichlet (V.replicate nState 1.0)
  let priorBeta = normalDistr 0 (5 * (fromIntegral nState - 1) :: Double)
  let nSamp      = 10
  -- let localStep  = 20
  let betaGrad = DistInvariant (fromIntegral $ V.length xs)
                               priorBeta
                               (rhoKuc defaultKucP)
  qBetas <- cellWith $ mergeGenericss globalMaxStep globalDelta
  write qBetas =<< V.generateM
    nDimension
    (\_i -> V.replicateM
      nState
      (do
        mu <- resample priorBeta gen1
        return $ defaultDistCell (normalDistr mu 1.0)
      )
    )
  qThetas <- cellWith $ mergeGenerics globalMaxStep globalDelta
  write qThetas $ V.replicate nTime (defaultDistCell priorTheta)
  let thetaGrad = DistInvariant (fromIntegral $ V.length xs)
                                priorTheta
                                (rhoKuc defaultKucP)
  stepTogether
    (updateMembershipReparam nSamp (mixtureJoint 1.0) gen1 xs thetaGrad betaGrad
    )
    qThetas
    qBetas
  -- stepSeparate localStep
  --              globalDelta
  --              (updateMembershipScore nSamp (mixtureJoint 1.0) gen1 xs thetaGrad betaGrad)
  --              qThetas
  --              qBetas
  thetaF <- unsafeContent qThetas
  betaF  <- unsafeContent qBetas
  return (V.map dist thetaF, V.map (V.map dist) betaF)

max_stars_repo_name: massma/propagator-bbvi

max_issues_repo_name: massma/propagator-bbvi

max_forks_repo_name: massma/propagator-bbvi


text:
module Eris.Predict.KNNbased where

import Data.HashMap.Strict as Map
import Data.Maybe (isNothing, isJust, fromJust)
import Control.Monad (guard)
import qualified Numeric.LinearAlgebra as NL

import Eris.Meta.DataStructure
import Eris.Compute.Similarity

-- | Example:
-- records2eCount :  get ecount
-- pairWiseSimilarity : get similarity Matrix
-- neighborOf : get neighbor of a entity  A
-- f1 : filter neighbors of A which related to entity of the other type .
-- c1 : retrive relation between two type of entity.
-- t : compute knnbased recommendataion.
t :: NeighborRank -> NeighborSimilarity-> Maybe Rank
t nr ns
    | length ns /= length nr = Nothing
    | otherwise = Just $ knnbased  ns nr
    where knnbased ns nr =  let
                              v1 = NL.vector $ fmap snd ns
                              v2 = NL.vector $ fmap snd nr
                            in  v1 NL.<.> v2 / NL.norm_1 v1

f1 :: PID -> ECount -> NeighborSimilarity -> NeighborSimilarity
f1 pid ecount nsim = do
    r@(eid,_) <- nsim
    guard $ validRecord ecount eid pid
    return r

c1 :: PID -> ECount -> NeighborSimilarity -> NeighborRank
c1 pid ecount nsim =do
                    (eid, _) <- nsim
                    let justRank = getRank ecount eid pid
                    guard $ isJust justRank
                    return (eid, fromJust justRank)

-- | Given an ID of a entity of certain type.
--  Based on the similairty matrix of this type of entity.
--  collect its neighbors and sort by Order.
--  return a list of tuple [(EID, Rank)] , rank value is necessary in case any further filter or verification process.
neighborOf :: Order -> CID -> SimilarityMatrix -> NeighborSimilarity
neighborOf ord cus matrix
  | isNothing (Map.lookup cus matrix) = []
  | otherwise = let
      rlist = Map.toList matrix
      neighborsRank = collectRank rlist cus
      in sortByRank neighborsRank ord
  where
        collectRank :: [(EID, Map.HashMap EID Rank)] -> EID -> NeighborSimilarity
        collectRank ((eid,edict) : exs) tid
            | eid == tid = Map.toList edict
            | otherwise = (eid,Map.lookupDefault 0 tid edict): collectRank exs tid
        sortByRank :: NeighborSimilarity -> Order -> NeighborSimilarity
        sortByRank nr o = let
            ascOrder = quickSortRank nr
            in if o == Desc
               then
                 reverse ascOrder
                else
                  ascOrder
        quickSortRank :: NeighborSimilarity -> NeighborSimilarity
        quickSortRank [] = []
        quickSortRank ((eid,erank):ers) =
            let smallRanks = quickSortRank [e | e <- ers , snd e <= erank]
                biggerRanks = quickSortRank [ e | e<- ers, snd e > erank]
            in smallRanks ++ [(eid,erank)] ++ biggerRanks

-- | Auxiliary Functions
validRecord :: ECount -> CID -> PID -> Bool
validRecord ec c p=
        let cContainsP = do
                          pdict <- Map.lookup c ec
                          Map.lookup p pdict
        in isJust cContainsP

getRank :: ECount -> CID -> PID -> Maybe Rank
getRank ec c p = do
    pdict <- Map.lookup c ec
    Map.lookup p pdict

max_stars_repo_name: emmettng/eris

max_issues_repo_name: emmettng/eris

max_forks_repo_name: emmettng/eris


text:
{- |
Description :  Some handy (mathematical, string manipulation, and so on) tools
Copyright   :  (c) Dominik Schrempf 2017
License     :  GPLv3

Maintainer  :  dominik.schrempf@gmail.com
Stability   :  unstable
Portability :  non-portable (not tested)

Please see function definitions and documentations.

* Changelog

-}

module Tools
  ( allValues
  , harmonic
  , matrixSeparateSymSkew
  , matrixSetDiagToZero
  , nearlyEq
  , left
  , right
  ) where

import Numeric.LinearAlgebra

-- | Get all values of a bounded enumerated type.
allValues :: (Bounded a, Enum a) => [a]
allValues = [minBound..]

-- | Calculate the nth harmonic number.
harmonic :: Int -> Double
harmonic 1 = 1.0
harmonic n = 1.0 / fromIntegral n + harmonic (n-1)

-- | Separate a square matrix into a symmetric and a skew-symmetric matrix.
matrixSeparateSymSkew :: Matrix R -> (Matrix R, Matrix R)
matrixSeparateSymSkew m = (mSym, mSkew)
  where trM = tr m
        mSym  = scale 0.5 $ m + trM
        mSkew = scale 0.5 $ m - trM

-- | Set the diagonal entries of a matrix to zero.
matrixSetDiagToZero :: Matrix R -> Matrix R
matrixSetDiagToZero m = m - diag (takeDiag m)

-- | Test for equality with tolerance (needed because of machine precision).
nearlyEq :: Double -> Double -> Double -> Bool
nearlyEq tol a b = tol > abs (a-b)

-- Functions that fill a string 's' to a given width 'n' by adding a pad
-- character 'c' (c) to align right.
fillDiff :: Int -> String -> Int
fillDiff width entry =
  if l >= width then 0 else width - l
  where l = length entry

fillLeft :: Char -> Int -> String -> String
fillLeft c n s = s ++ replicate (fillDiff n s) c

fillRight :: Char -> Int -> String -> String
fillRight c n s = replicate (fillDiff n s) c ++ s

-- | Fill a string to a given width by adding spaces. Align left.
left :: Int -> String -> String
left = fillLeft ' '

-- | Fill a string to a given width by adding spaces. Align right.
right :: Int -> String -> String
right = fillRight ' '

max_stars_repo_name: dschrempf/bmm-simulate

max_issues_repo_name: pomo-dev/bmm-simulate

max_forks_repo_name: pomo-dev/bmm-simulate


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ForeignFunctionInterface #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}

module Main where

import Foreign (Ptr, alloca, poke, peek, Storable)
import Foreign.C.Types (CInt(..), CDouble(..))
import Foreign.Ptr (Ptr)
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Devel
import System.IO.Unsafe (unsafePerformIO)
import Text.Printf (printf)

--------------------------------------------------
-- * 0) Simple.
--
-- Import straight from C's stdio (IE not any custom C file).

foreign import ccall "exp" myExp :: Double -> Double

--------------------------------------------------
-- * 1) Hello.
--
-- Demonstrates IO from Fortran with no parameter passing.

foreign import ccall "hello_" hello :: IO ()

--------------------------------------------------
-- * 2) FPow.
--
-- Demonstrates Fortran parameter passing.

-- this is a pure function that wraps the Fortran function. The Fortran function
-- expects pointers, and this C wrapper takes values, and handles passing on the
-- references to the Fortran function.
foreign import ccall "fpow" fpow :: Double -> Double -> Double

-- This function is the unwrapped Fortran function, FFI'd into C, and it expects
-- pointers, not values.
foreign import ccall "fpow_" fpow_ :: Ptr Double -> Ptr Double -> IO Double

fpowViaRefs :: Double -> Double -> IO Double
fpowViaRefs b e =
  alloca $ \base -> -- get a pointer
  alloca $ \exponent -> do
    poke base b -- set its value
    poke exponent e
    fpow_ base exponent -- note: pointers don't escape, just the final value

--------------------------------------------------
-- * 3) Transpose.
--
-- Demonstrates interfacing between C and HMatrix. C and HMatrix both use
-- row-major arrays.

foreign import ccall "transpose"
  transpose_ :: Double ::> Double ::> IO CInt

transpose :: Matrix Double -> Matrix Double
transpose m = let
  nrow = rows m
  ncol = cols m
  in unsafePerformIO $ do
  -- allocate some memory for the output matrix
  m' <- createMatrix RowMajor ncol nrow
  -- apply function to memory locations
  () <- (m #! m') transpose_ #| "transpose_"
  pure m'

--------------------------------------------------
-- * 4) Mat * Scalar multiplication.
--
-- Demonstrates interfacing between Fortran and HMatrix. Fortran uses
-- colum-major arrays, so, HMatrix can correct that with `fmat`.

foreign import ccall "scalarmul"
  scalarmul :: Double ::> (Double -> IO CInt)

scalarMul :: Matrix Double -> Double -> Matrix Double
scalarMul m x = unsafePerformIO $ do
  let m' = fmat m -- copy array into Fortran's col-major ordering
  () <- (m' #! x) scalarmul #| "scalarMul"
  -- TODO: m' mutated. this is ugly, is there a better way?
  pure  m'

--------------------------------------------------
-- * Main.

main :: IO ()
main = do
  putStrLn "0_simple: `exp` comes from C's stdio:"
  putStrLn $ printf "myExp: e^3=%3.3f" (myExp 3)

  putStrLn "\n1_hello: Fortran can print to stdout on its own:"
  hello

  putStrLn "\n2_fpow: Passing data into Fortran:"
  putStrLn $ printf "pure fpow: %3.3f" (fpow 2 12)

  putStrLn "\n2_fpow: Passing pointers into Fortran:"
  putStrLn . printf "impure fpow: %3.3f" =<< fpowViaRefs 2 13

  putStrLn "\n3_transpose: Connecting C and HMatrix:"
  let m = (10><10) [0 :: Double ..]
      -- testing slices to make the problem more interesting, since HMatrix may
      -- still retain the larger array, but records the dimensions of the sliced
      -- array separately.
      sliceM = subMatrix (3, 3) (5, 7) m
  print sliceM
  print $ transpose sliceM

  putStrLn "\n4_scalarmul: Connecting Fortran and HMatrix:"
  let m2 = (10><10) [0 :: Double ..]
      sliceM2 = subMatrix (3, 3) (5, 7) m2
  print $ scalarMul sliceM2 1000

--------------------------------------------------
-- * HMatrix.Devel helpers

infixr 1 #
(#) :: TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
a # b = apply a b
{-# INLINE (#) #-}

-- | Apply vectors and matrices to a c function
--
-- Usage Ex1:
--   (param1 # param2 # param3 #! param4) c_function #| "c_function"
(#!) :: (TransArray c, TransArray c1) => c1 -> c -> Trans c1 (Trans c (IO r)) -> IO r
a #! b = a # b # id
{-# INLINE (#!) #-}

-- This idea taken from the HMatrix library. TransArray instances, together with
-- `apply`, will set multiple parameters to pass into C.
infixr 5 :>, ::>

-- | Type for vector params passed to C
type (:>)  t r = CInt -- size of vector
              -> Ptr t -- pointer to vector
              -> r -- a continuation, allowing composition of params

-- | Type for matrix params passed to C
type (::>) t r =  CInt -- slice height (nrows)
               -> CInt -- slice width (ncols)
               -> CInt -- matrix height
               -> CInt -- matrix width
               -> Ptr t -- array pointer
               -> r -- a continuation, allowing composition of params

instance TransArray Double where
  type TransRaw Double b = Double -> b
  type Trans Double b = Double -> b
  apply x f g = f (g x)
  applyRaw x f g = f (g x)

instance TransArray Int where
  type TransRaw Int b = Int -> b
  type Trans Int b = Int -> b
  apply x f g = f (g x)
  applyRaw x f g = f (g x)

max_stars_repo_name: freckletonj/haskell-fortran

max_issues_repo_name: freckletonj/haskell-fortran

max_forks_repo_name: freckletonj/haskell-fortran


text:
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ViewPatterns #-}
{-# LANGUAGE MultiWayIf #-}

module Data.AdaptiveCoordinateDescent
  ( adaptiveCoordinateDescent )
  where

import Control.Monad.State.Strict
import Data.AdaptiveCoordinateDescent.Internal
import Data.Foldable
import Data.List ( sortBy )
import Data.Ord
import Data.Traversable
import Numeric.LinearAlgebra hiding ( toList )
import Pipes hiding ( for )
import qualified Pipes.Prelude as P
import System.Random.MWC

-- | Performs adaptive coordinate descent.
adaptiveCoordinateDescent :: (MonadIO m, Traversable f, Eq (f Double))
                          => (f Double -> m Double)  -- ^ Cost function. Smaller cost = better model.
                          -> f Double                -- ^ Initial parameters from which to start optimization.
                          -> Double                  -- ^ Tolerance on how small changes are allowed.
                                                     --   If no changes larger
                                                     --   than this value are
                                                     --   found, then this
                                                     --   descent function
                                                     --   returns. You probably
                                                     --   want something small
                                                     --   like 0.0000001
                          -> Int                     -- ^ How many best N candidates to keep for PCA phase.
                                                     --   Good values are
                                                     --   probably around the
                                                     --   size of number of
                                                     --   your parameters.
                                                     --   Maybe more if you
                                                     --   only have like 2
                                                     --   parameters. Probably.
                                                     --   We don't know.
                          -> Producer (f Double, Double) m (f Double, Double)
adaptiveCoordinateDescent evaluate initial_params tolerance how_many_to_keep = do
  score <- lift $ evaluate initial_params
  yield (initial_params, score)

  loop_it initial_params initial_params score (toColumns $ ident $ length $ toList initial_params) [] (replicate (length $ toList initial_params) 1.0) []
 where
  k_succ = 2.0 :: Double
  k_unsucc = 0.5 :: Double

  loop_it original_params params score principal_components last_n_params step_sizes new_step_sizes | length last_n_params < 5 = do
    rng <- liftIO createSystemRandom
    items <- liftIO $ P.toListM (makePerturbed params rng 5)
    scored_items <- lift $ for (initial_params:items) $ \item -> (,) item <$> evaluate item
    loop_it original_params params score principal_components scored_items step_sizes new_step_sizes

  loop_it original_params params score (principal_component:principal_components) last_n_params (step_size':step_sizes) new_step_sizes = do
    let step_size = max tolerance step_size'
        scaled_pcomponent = cmap (*step_size) principal_component
        top_candidate    = params `plus` fromVector scaled_pcomponent initial_params
        bottom_candidate = params `minus` fromVector scaled_pcomponent initial_params

    top_score    <- lift $ evaluate top_candidate
    bottom_score <- lift $ evaluate bottom_candidate

    yield (top_candidate, top_score)
    yield (bottom_candidate, bottom_score)

    if | top_score < score
         -> loop_it original_params top_candidate top_score principal_components ((top_candidate, top_score):last_n_params) step_sizes (step_size*k_succ:new_step_sizes)
       | bottom_score < score
         -> loop_it original_params bottom_candidate bottom_score principal_components ((bottom_candidate, bottom_score):last_n_params) step_sizes (step_size*k_succ:new_step_sizes)
       | otherwise
         -> loop_it original_params params score principal_components last_n_params step_sizes (step_size*k_unsucc:new_step_sizes)

  loop_it original_params params score [] _last_n_params _step_sizes new_step_sizes
    | original_params == params &&
      all (\x -> x <= tolerance) new_step_sizes = yield (params, score) >> return (params, score)
  loop_it _original_params params score [] last_n_params _step_sizes new_step_sizes = do
    let num_coordinates = length (toList params)
        sorted_items = sortBy (comparing snd) last_n_params
        picked_items = take how_many_to_keep sorted_items

        items_mat = normalize $
                    ( (length picked_items><num_coordinates)
                      (mconcat $ fmap (toList . fst) picked_items) )

    let (_evec, mat2) = rightSV items_mat
    loop_it params params score (toColumns mat2) picked_items (reverse new_step_sizes) []

  loop_it _ _ _ _ _ _ _ = error "adaptiveCoordinateDescent: impossible case."

max_stars_repo_name: Noeda/adaptivecoordinatedescent

max_issues_repo_name: Noeda/adaptivecoordinatedescent

max_forks_repo_name: Noeda/adaptivecoordinatedescent


text:
{-# LANGUAGE ExistentialQuantification #-}
module LispTypes where

import           Control.Monad.Except
import           Text.ParserCombinators.Parsec
import           Data.Complex
import           Data.Array
import           Data.IORef
import           Text.Printf
import           System.IO

type IOThrowsError = ExceptT LispError IO

type Env = IORef [(String, IORef LispVal)]

data Unpacker = forall a. Eq a => AnyUnpacker (LispVal -> ThrowsError a)

data LispVal = Atom String
            | List [LispVal]
            | DottedList [LispVal] LispVal
            | Number Integer
            | String String
            | Bool Bool
            | Character Char
            | Float Double
            | Ratio Rational
            | Complex (Complex Double)
            | Vector (Array Int LispVal)
            | PrimitiveFunc ([LispVal] -> ThrowsError LispVal)
            | Func { params :: [String], vararg :: (Maybe String), body :: [LispVal], closure :: Env }
            | IOFunc ([LispVal] -> IOThrowsError LispVal)
            | Port Handle

instance Eq LispVal where
  (Atom x         ) == (Atom y         ) = x == y
  (List x         ) == (List y         ) = x == y
  (DottedList x x') == (DottedList y y') = x == y && x' == y'
  (Number    x    ) == (Number    y    ) = x == y
  (String    x    ) == (String    y    ) = x == y
  (Bool      x    ) == (Bool      y    ) = x == y
  (Character x    ) == (Character y    ) = x == y
  (Float     x    ) == (Float     y    ) = x == y
  (Ratio     x    ) == (Ratio     y    ) = x == y
  (Complex   x    ) == (Complex   y    ) = x == y
  (Vector    x    ) == (Vector    y    ) = x == y

data LispError = NumArgs Integer [LispVal]
                | TypeMismatch String LispVal
                | Parser ParseError
                | BadSpecialForm String LispVal
                | NotFunction String String
                | UnboundVar String String
                | Default String

unwordsVal :: [LispVal] -> String
unwordsVal = unwords . map showVal

showVal :: LispVal -> String
showVal (String contents     ) = "\"" ++ contents ++ "\""
showVal (Atom   name         ) = name
showVal (Number contents     ) = show contents
showVal (Bool   True         ) = "#t"
showVal (Bool   False        ) = "#f"
showVal (List   contents     ) = "(" ++ unwordsVal contents ++ ")"
showVal (DottedList head tail) = "(" ++ unwordsVal head ++ showVal tail ++ ")"
showVal (PrimitiveFunc _     ) = "<primitive>"
showVal (Func { params = args, vararg = varargs, body = body, closure = env })
  = "(lambda ("
    ++ unwords (map show args)
    ++ (case varargs of
         Nothing  -> ""
         Just arg -> " . " ++ arg
       )
    ++ ") ...)"
showVal (Port   _) = "<IO Port>"
showVal (IOFunc _) = "<IO primitive>"

showError :: LispError -> String
showError (UnboundVar     message  varname) = message ++ ": " ++ varname
showError (BadSpecialForm message  form   ) = message ++ ": " ++ show form
showError (NotFunction    message  func   ) = message ++ ": " ++ show func
showError (NumArgs        expected found  ) = printf
  "Expected %i args; found %i values %s"
  expected
  (length found)
  (unwordsVal found)
showError (TypeMismatch expected found) =
  "Invalid type: expected " ++ expected ++ ", found " ++ show found
showError (Parser  parseErr) = "Parse error at " ++ show parseErr
showError (Default s       ) = s

instance Show LispError where
  show = showError

instance Show LispVal where
  show = showVal

type ThrowsError = Either LispError


text:
module Types
(Data (..),
 plus,
 divide,
 multiply,
 generateBasedOf,
 isEmpty,
 hadamardV,
 hadamardM,
 sigV,
 sigV'
) where

import Numeric.LinearAlgebra.HMatrix

-- Construtor do Data
data Data = Data {
                wHidden :: Matrix R,
                bHidden :: Vector Double,
                aHidden :: Vector Double,
                zetaHidden :: Vector Double,
                wOutput :: Matrix R,
                bOutput :: Vector Double,
                aOutput :: Vector Double,
                zetaOutput :: Vector Double
                } deriving (Eq, Show)

-- Efetua soma entre os respectivos vetores e matrizes de ambos os Datas
plus :: Data -> Data -> Data
plus (Data wH1 bH1 aH1 zH1 wO1 bO1 aO1 zO1) (Data wH2 bH2 aH2 zH2 wO2 bO2 a02 z02) =
    (Data (add wH1 wH2) (add bH1 bH2) aH1 zH1 (add wO1 wO2) (add bO1 bO2) aO1 zO1)

-- Efetua divisao por esacalar entre os respectivos vetores e matrizes do Data
divide :: Data -> Double -> Data
divide (Data wH bH aH zH wO bO aO zO) constant = (Data (scale (1/constant) wH) (scale (1/constant) bH) aH zH (scale (1/constant) wO) (scale (1/constant) bO) aO zO)

-- Efetua multiplicacao por scalar entre os respectivos vetores e matrizes do Data
multiply :: Data -> Double -> Data
multiply (Data wH bH aH zH wO bO aO zO) constant = (Data (scale constant wH) (scale constant bH) aH zH (scale constant wO) (scale constant bO) aO zO)

-- Efetua hadamard entre dois vetores
hadamardV :: (Num a) => [a] -> [a] -> [a]
hadamardV list1 list2 = zipWith (*) list1 list2

-- Efetua hadamard entre duas matrizes
hadamardM :: (Num a) => [[a]] -> [[a]] -> [[a]]
hadamardM matrix1 matrix2 = zipWith hadamardV matrix1 matrix2

-- Computa a funcao sigmoid daquele valor
sig :: Double -> Double
sig x = 1.0 - (1.0 + exp(-x))

-- Computa a derivada da funcao sigmoid daquele valor
sig' :: Double -> Double
sig' x = sig x * (1 - sig x)

-- Computa a funcao sigmoid para cada elemento de uma lista
sigV :: [Double] -> [Double]
sigV l = [sig x | x <- l]

-- Computa a derivada da funcao sigmoid para cada elemento de uma lista
sigV' :: [Double] -> [Double]
sigV' l = [sig' x | x <- l]

-- Cria um data repleto de valores 0.0 com as mesmas dimensoes do modelo passado por parametro
generateBasedOf :: Data -> Data
generateBasedOf (Data wH bH aH zH wO bO aO zO) = let whRows = fst $ size wH
                                                     whColumns = snd $ size wH
                                                     bhLen = size bH
                                                     ahLen = size aH
                                                     zhLen = size zH
                                                     woRows = fst $ size wO
                                                     woColumns = snd $ size wO
                                                     boLen = size bO
                                                     aoLen = size aH
                                                     zoLen = size zH

                                                     newWH = (whRows><whColumns)(replicate (whRows * whColumns) 0.0) :: Matrix R
                                                     newBH = fromList $ replicate bhLen 0.0
                                                     newAH = fromList $ replicate ahLen 0.0
                                                     newZH = fromList $ replicate zhLen 0.0
                                                     newWO = (woRows><woColumns)(replicate (woRows * woColumns) 0.0) :: Matrix R
                                                     newBO = fromList $ replicate boLen 0.0
                                                     newAO = fromList $ replicate aoLen 0.0
                                                     newZO = fromList $ replicate zoLen 0.0
                                                 in (Data newWH newBH newAH newZH newWO newBO newAO newZO)

-- Checa se um Data tem algum elemento nulo
isEmpty :: Data -> Bool
isEmpty (Data wH bH aH zH wO bO aO zO) =
    (null (toList $ flatten wH)) || (null $ toList bH) || (null (toList $ flatten wO)) || (null $ toList bO)


text:
module Scene.Material.Skylike where

import Numeric.LinearAlgebra

import Scene.Base
import Math.Color
import Math.Ray
import Scene.Monad

import Control.Applicative
import Control.Monad.Random

skylike :: Color -> Color -> Material
skylike near far (Ray p normal) (Ray q incidence) = do
    let t = abs $ (normal <.> incidence)/(norm_2 normal * norm_2 incidence)
    return $! weigh [(t,near),(1-t,far)]

max_stars_repo_name: danrocag/raytracer

max_issues_repo_name: danrocag/raytracer

max_forks_repo_name: danrocag/raytracer


text:
module IOUtils where

import System.IO (FilePath, readFile, writeFile)
import Data.Complex (Complex)

parseLines :: [String] -> [Double]
parseLines ls = map read ls :: [Double]

readData :: FilePath -> IO [Double]
readData filePath = do
    content <- readFile filePath
    return $ parseLines $ lines content

writeData :: FilePath -> [Double] -> IO ()
writeData filePath d = writeFile filePath $ unlines $ map show d

printData :: Show a => [a] -> IO ()
printData d = putStrLn $ unlines $ map (++ ", ") $ map show d

max_stars_repo_name: vined/fourier

max_issues_repo_name: vined/fourier

max_forks_repo_name: vined/fourier


text:
{-|
Module      : kMeans
Description : k-Means
Copyright   : (c) Fabrício Olivetti, 2017
License     : GPL-3
Maintainer  : fabricio.olivetti@gmail.com

k-Means Clustering algorithm
-}

module KMeansPar where

import Control.Parallel.Strategies
import Data.List
import Data.Ord
import Data.List.Split (chunksOf)
--import qualified Numeric.LinearAlgebra as H

import Vector
import Dados

-- |'parseFile' parses a space separated file
-- to a list of lists of Double
parseFile :: String -> [[Double]]
parseFile file = map parseLine (lines file)
  where
    parseLine l = map toDouble (words l)
    toDouble  w = read w :: Double

firstClusters idxs points = [ points !! i | i <- idxs ]

kmeans :: Int -> ChunksOf [[Double]] -> [[Double]]
       -> [[Double]]
kmeans it points clusters
  | it == 0               = clusters
  | clusters' == clusters = clusters
  | otherwise             = kmeans (it-1) points clusters'
  where
    clusters' = emStep points clusters

emStep :: ChunksOf [[Double]] -> [[Double]] -> [[Double]]
emStep points clusters = map (\(idx, (xs,n)) -> xs ./ n)
                       $ mapReduceByKey (\pi -> (assign pi clusters, (pi,1.0))) soma points
  where
    soma (xs1, y1) (xs2, y2) = (xs1 .+. xs2, y1+y2)

assign :: [Double] -> [[Double]] -> Integer
assign point clusters = closestTo point
  where
    closestTo p = argmin $ map (euclid p) clusters
    euclid pi ci = sum $ map (^2) $ zipWith (-) pi ci
    argmin xs = fst $ head $ sortByValue $ zip [0..] xs

max_stars_repo_path: 09 - Agrupamento de Dados/KMeansPar.hs

max_stars_repo_name: douggribeiro/bigdata

max_issues_repo_path: 09 - Agrupamento de Dados/KMeansPar.hs

max_issues_repo_name: douggribeiro/bigdata

max_forks_repo_path: 09 - Agrupamento de Dados/KMeansPar.hs

max_forks_repo_name: douggribeiro/bigdata


text:
{-# LANGUAGE ScopedTypeVariables, BangPatterns, RecordWildCards #-}

{-| Functions used internally by the SDR.Filter module. Most of these are
    not actually used but exist for benchmarking purposes to determine the
    fastest filter implementation.
-}
module SDR.FilterInternal where

import           Control.Monad.Primitive
import           Control.Monad
import           Foreign.C.Types
import           Foreign.Ptr
import           Data.Coerce
import           Data.Complex
import           Foreign.Marshal.Array
import           Foreign.Marshal.Alloc
import           Foreign.Storable

import qualified Data.Vector.Generic               as VG
import qualified Data.Vector.Generic.Mutable       as VGM
import qualified Data.Vector.Storable              as VS
import qualified Data.Vector.Storable.Mutable      as VSM
import qualified Data.Vector.Fusion.Bundle         as VFB

import           SDR.VectorUtils
import           SDR.Util

{-# INLINE filterHighLevel #-}
filterHighLevel :: (PrimMonad m, Functor m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) => v b -> Int -> v a -> vm (PrimState m) a -> m ()
filterHighLevel coeffs num inBuf outBuf = fill (VFB.generate num dotProd) outBuf
    where
    dotProd offset = VG.sum $ VG.zipWith mult (VG.unsafeDrop offset inBuf) coeffs

{-# INLINE filterImperative1 #-}
filterImperative1 :: (PrimMonad m, Functor m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) => v b -> Int -> v a -> vm (PrimState m) a -> m ()
filterImperative1 coeffs num inBuf outBuf = go 0
    where
    go offset
        | offset < num = do
            let res = dotProd offset
            VGM.unsafeWrite outBuf offset res
            go $ offset + 1
        | otherwise    = return ()
    dotProd offset = VG.sum $ VG.zipWith mult (VG.unsafeDrop offset inBuf) coeffs

{-# INLINE filterImperative2 #-}
filterImperative2 :: (PrimMonad m, Functor m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) => v b -> Int -> v a -> vm (PrimState m) a -> m ()
filterImperative2 coeffs num inBuf outBuf = go 0
    where
    go offset
        | offset < num = do
            let res = dotProd (VG.unsafeDrop offset inBuf)
            VGM.unsafeWrite outBuf offset res
            go $ offset + 1
        | otherwise    = return ()
    dotProd buf = go 0 0
        where
        go !accum j
            | j < VG.length coeffs = go (VG.unsafeIndex buf j `mult` VG.unsafeIndex coeffs j  + accum) (j + 1)
            | otherwise            = accum

type FilterCRR = CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()
type FilterRR  = VS.Vector Float -> Int -> VS.Vector Float -> VS.MVector RealWorld Float -> IO ()
type FilterRC  = VS.Vector Float -> Int -> VS.Vector (Complex Float) -> VS.MVector RealWorld (Complex Float) -> IO ()

filterFFIR :: FilterCRR -> FilterRR
filterFFIR func coeffs num inBuf outBuf =
    VS.unsafeWith (coerce coeffs) $ \cPtr ->
        VS.unsafeWith (coerce inBuf) $ \iPtr ->
            VSM.unsafeWith (coerce outBuf) $ \oPtr ->
                func (fromIntegral num) (fromIntegral $ VG.length coeffs) cPtr iPtr oPtr

filterFFIC :: FilterCRR -> FilterRC
filterFFIC func coeffs num inBuf outBuf =
    VS.unsafeWith (coerce coeffs) $ \cPtr ->
        VS.unsafeWith (coerce inBuf) $ \iPtr ->
            VSM.unsafeWith (coerce outBuf) $ \oPtr ->
                func (fromIntegral num) (fromIntegral $ VG.length coeffs) cPtr iPtr oPtr

foreign import ccall unsafe "filterRR"
    filterRR_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

filterCRR :: FilterRR
filterCRR = filterFFIR filterRR_c

foreign import ccall unsafe "filterRC"
    filterRC_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

filterCRC :: FilterRC
filterCRC = filterFFIC filterRC_c

foreign import ccall unsafe "filterSSERR"
    filterSSERR_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

filterCSSERR :: FilterRR
filterCSSERR = filterFFIR filterSSERR_c

foreign import ccall unsafe "filterSSESymmetricRR"
    filterSSESymmetricRR_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

filterCSSESymmetricRR :: FilterRR
filterCSSESymmetricRR = filterFFIR filterSSESymmetricRR_c

foreign import ccall unsafe "filterSSERC"
    filterSSERC_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

filterCSSERC :: FilterRC
filterCSSERC = filterFFIC filterSSERC_c

foreign import ccall unsafe "filterSSERC2"
    filterSSERC2_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

filterCSSERC2 :: FilterRC
filterCSSERC2 = filterFFIC filterSSERC2_c

foreign import ccall unsafe "filterAVXRR"
    filterAVXRR_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

filterCAVXRR :: FilterRR
filterCAVXRR = filterFFIR filterAVXRR_c

foreign import ccall unsafe "filterAVXSymmetricRR"
    filterAVXSymmetricRR_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

filterCAVXSymmetricRR :: FilterRR
filterCAVXSymmetricRR = filterFFIR filterAVXSymmetricRR_c

foreign import ccall unsafe "filterAVXRC"
    filterAVXRC_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

filterCAVXRC :: FilterRC
filterCAVXRC = filterFFIC filterAVXRC_c

foreign import ccall unsafe "filterAVXRC2"
    filterAVXRC2_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

filterCAVXRC2 :: FilterRC
filterCAVXRC2 = filterFFIC filterAVXRC2_c

foreign import ccall unsafe "filterSSESymmetricRC"
    filterSSESymmetricRC_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

filterCSSESymmetricRC :: FilterRC
filterCSSESymmetricRC = filterFFIC filterSSESymmetricRC_c

foreign import ccall unsafe "filterAVXSymmetricRC"
    filterAVXSymmetricRC_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

filterCAVXSymmetricRC :: FilterRC
filterCAVXSymmetricRC = filterFFIC filterAVXSymmetricRC_c

-- Decimation

{-# INLINE decimateHighLevel #-}
decimateHighLevel :: (PrimMonad m, Functor m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) => Int -> v b -> Int -> v a -> vm (PrimState m) a -> m ()
decimateHighLevel factor coeffs num inBuf outBuf = fill x outBuf
    where
    x = VFB.map dotProd (VFB.iterateN num (+ factor) 0)
    dotProd offset = VG.sum $ VG.zipWith mult (VG.unsafeDrop offset inBuf) coeffs

type DecimateCRR = CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()
type DecimateRR  = Int -> VS.Vector Float -> Int -> VS.Vector Float -> VS.MVector RealWorld Float -> IO ()
type DecimateRC  = Int -> VS.Vector Float -> Int -> VS.Vector (Complex Float) -> VS.MVector RealWorld (Complex Float) -> IO ()

decimateFFIR :: DecimateCRR -> DecimateRR
decimateFFIR func factor coeffs num inBuf outBuf =
    VS.unsafeWith (coerce coeffs) $ \cPtr ->
        VS.unsafeWith (coerce inBuf) $ \iPtr ->
            VSM.unsafeWith (coerce outBuf) $ \oPtr ->
                func (fromIntegral num) (fromIntegral factor) (fromIntegral $ VG.length coeffs) cPtr iPtr oPtr

decimateFFIC :: DecimateCRR -> DecimateRC
decimateFFIC func factor coeffs num inBuf outBuf =
    VS.unsafeWith (coerce coeffs) $ \cPtr ->
        VS.unsafeWith (coerce inBuf) $ \iPtr ->
            VSM.unsafeWith (coerce outBuf) $ \oPtr ->
                func (fromIntegral num) (fromIntegral factor) (fromIntegral $ VG.length coeffs) cPtr iPtr oPtr

foreign import ccall unsafe "decimateRR"
    decimateCRR_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

decimateCRR :: DecimateRR
decimateCRR = decimateFFIR decimateCRR_c

foreign import ccall unsafe "decimateRC"
    decimateCRC_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

decimateCRC :: DecimateRC
decimateCRC = decimateFFIC decimateCRC_c

foreign import ccall unsafe "decimateSSERR"
    decimateSSERR_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

decimateCSSERR :: DecimateRR
decimateCSSERR = decimateFFIR decimateSSERR_c

foreign import ccall unsafe "decimateSSERC"
    decimateSSERC_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

decimateCSSERC :: DecimateRC
decimateCSSERC = decimateFFIC decimateSSERC_c

foreign import ccall unsafe "decimateSSERC2"
    decimateSSERC2_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

decimateCSSERC2 :: DecimateRC
decimateCSSERC2 = decimateFFIC decimateSSERC2_c

foreign import ccall unsafe "decimateAVXRR"
    decimateAVXRR_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

decimateCAVXRR :: DecimateRR
decimateCAVXRR = decimateFFIR decimateAVXRR_c

foreign import ccall unsafe "decimateAVXRC"
    decimateAVXRC_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

decimateCAVXRC :: DecimateRC
decimateCAVXRC = decimateFFIC decimateAVXRC_c

foreign import ccall unsafe "decimateAVXRC2"
    decimateAVXRC2_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

decimateCAVXRC2 :: DecimateRC
decimateCAVXRC2 = decimateFFIC decimateAVXRC2_c

foreign import ccall unsafe "decimateSSESymmetricRR"
    decimateSSESymmetricRR_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

decimateCSSESymmetricRR :: DecimateRR
decimateCSSESymmetricRR = decimateFFIR decimateSSESymmetricRR_c

foreign import ccall unsafe "decimateAVXSymmetricRR"
    decimateAVXSymmetricRR_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

decimateCAVXSymmetricRR :: DecimateRR
decimateCAVXSymmetricRR = decimateFFIR decimateAVXSymmetricRR_c

foreign import ccall unsafe "decimateSSESymmetricRC"
    decimateSSESymmetricRC_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

decimateCSSESymmetricRC :: DecimateRC
decimateCSSESymmetricRC = decimateFFIC decimateSSESymmetricRC_c

foreign import ccall unsafe "decimateAVXSymmetricRC"
    decimateAVXSymmetricRC_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

decimateCAVXSymmetricRC :: DecimateRC
decimateCAVXSymmetricRC = decimateFFIC decimateAVXSymmetricRC_c

-- Resampling
{-# INLINE resampleHighLevel #-}
resampleHighLevel :: (PrimMonad m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) => Int -> Int -> v b -> Int -> Int -> v a -> vm (PrimState m) a -> m Int
resampleHighLevel interpolation decimation coeffs filterOffset count inBuf outBuf = fill 0 filterOffset 0
    where
    fill i filterOffset inputOffset
        | i < count = do
            let dp = dotProd filterOffset inputOffset
            VGM.unsafeWrite outBuf i dp
            let (q, r)        = divMod (decimation - filterOffset - 1) interpolation
                inputOffset'  = inputOffset + q + 1
                filterOffset' = interpolation - 1 - r
            filterOffset' `seq` inputOffset' `seq` fill (i + 1) filterOffset' inputOffset'
        | otherwise = return filterOffset
    dotProd filterOffset offset = VG.sum $ VG.zipWith mult (VG.unsafeDrop offset inBuf) (stride interpolation (VG.unsafeDrop filterOffset coeffs))

foreign import ccall unsafe "resampleRR"
    resample_c :: CInt -> CInt -> CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

resampleCRR :: Int -> Int -> Int -> Int -> VS.Vector Float -> VS.Vector Float -> VS.MVector RealWorld Float -> IO ()
resampleCRR num interpolation decimation offset coeffs inBuf outBuf =
    VS.unsafeWith (coerce coeffs) $ \cPtr ->
        VS.unsafeWith (coerce inBuf) $ \iPtr ->
            VSM.unsafeWith (coerce outBuf) $ \oPtr ->
                resample_c (fromIntegral num) (fromIntegral $ VG.length coeffs) (fromIntegral interpolation) (fromIntegral decimation) (fromIntegral offset) cPtr iPtr oPtr

pad :: a -> Int -> [a] -> [a]
pad with num list = list ++ replicate (num - length list) with

strideList :: Int -> [a] -> [a]
strideList s xs = go 0 xs
    where
    go _ []     = []
    go 0 (x:xs) = x : go (s-1) xs
    go n (x:xs) = go (n - 1) xs

roundUp :: Int -> Int -> Int
roundUp num div = ((num + div - 1) `quot` div) * div

data Coeffs = Coeffs {
    numCoeffs  :: Int,
    numGroups  :: Int,
    increments :: [Int],
    groups     :: [[Float]]
}

prepareCoeffs :: Int -> Int -> Int -> [Float] -> Coeffs
prepareCoeffs n interpolation decimation coeffs = Coeffs {..}
    where
    numCoeffs   = maximum $ map (length . snd) dats
    numGroups   = length groups
    increments  = map fst dats

    groups      :: [[Float]]
    groups      = map (pad 0 (roundUp numCoeffs n)) $ map snd dats

    dats :: [(Int, [Float])]
    dats = func 0
        where

        func' 0      = []
        func' x      = func x

        func :: Int -> [(Int, [Float])]
        func offset = (increment, strideList interpolation $ drop offset coeffs) : func' offset'
            where
            (q, r)    = divMod (decimation - offset - 1) interpolation
            increment = q + 1
            offset'   = interpolation - 1 - r

resampleFFIR :: (Ptr CFloat -> Ptr CFloat -> IO CInt) -> VS.Vector Float -> VSM.MVector RealWorld Float -> IO Int
resampleFFIR func inBuf outBuf = liftM fromIntegral $
    VS.unsafeWith (coerce inBuf) $ \iPtr ->
        VSM.unsafeWith (coerce outBuf) $ \oPtr ->
            func iPtr oPtr

resampleFFIC :: (Ptr CFloat -> Ptr CFloat -> IO CInt) -> VS.Vector (Complex Float) -> VSM.MVector RealWorld (Complex Float) -> IO Int
resampleFFIC func inBuf outBuf = liftM fromIntegral $
    VS.unsafeWith (coerce inBuf) $ \iPtr ->
        VSM.unsafeWith (coerce outBuf) $ \oPtr ->
            func iPtr oPtr

type ResampleR = CInt -> CInt -> CInt -> CInt -> Ptr CInt -> Ptr (Ptr CFloat) -> Ptr CFloat -> Ptr CFloat -> IO CInt

mkResampler :: ResampleR -> Int -> Int -> Int -> [Float] -> IO (Int -> Int -> VS.Vector Float -> VS.MVector RealWorld Float -> IO Int)
mkResampler func n interpolation decimation coeffs = do
    groupsP     <- mapM newArray $ map (map realToFrac) groups
    groupsPP    <- newArray groupsP
    incrementsP <- newArray $ map fromIntegral increments
    return $ \offset num -> resampleFFIR $ func (fromIntegral num) (fromIntegral numCoeffs) (fromIntegral offset) (fromIntegral numGroups) incrementsP groupsPP
    where
    Coeffs {..} = prepareCoeffs n interpolation decimation coeffs

type ResampleRR = Int -> Int -> [Float] -> IO (Int -> Int -> VS.Vector Float -> VS.MVector RealWorld Float -> IO Int)

foreign import ccall unsafe "resample2RR"
    resample2_c :: CInt -> CInt -> CInt -> CInt -> Ptr CInt -> Ptr (Ptr CFloat) -> Ptr CFloat -> Ptr CFloat -> IO CInt

resampleCRR2 :: ResampleRR
resampleCRR2 = mkResampler resample2_c 1

foreign import ccall unsafe "resampleSSERR"
    resampleCSSERR_c :: CInt -> CInt -> CInt -> CInt -> Ptr CInt -> Ptr (Ptr CFloat) -> Ptr CFloat -> Ptr CFloat -> IO CInt

resampleCSSERR :: ResampleRR
resampleCSSERR = mkResampler resampleCSSERR_c 4

foreign import ccall unsafe "resampleAVXRR"
    resampleAVXRR_c :: CInt -> CInt -> CInt -> CInt -> Ptr CInt -> Ptr (Ptr CFloat) -> Ptr CFloat -> Ptr CFloat -> IO CInt

resampleCAVXRR :: ResampleRR
resampleCAVXRR = mkResampler resampleAVXRR_c 8

type ResampleRC = Int -> Int -> [Float] -> IO (Int -> Int -> VS.Vector (Complex Float) -> VS.MVector RealWorld (Complex Float) -> IO Int)

mkResamplerC :: ResampleR -> Int -> Int -> Int -> [Float] -> IO (Int -> Int -> VS.Vector (Complex Float) -> VS.MVector RealWorld (Complex Float) -> IO Int)
mkResamplerC func n interpolation decimation coeffs = do
    groupsP     <- mapM newArray $ map (map realToFrac) groups
    groupsPP    <- newArray groupsP
    incrementsP <- newArray $ map fromIntegral increments
    return $ \offset num -> resampleFFIC $ func (fromIntegral num) (fromIntegral numCoeffs) (fromIntegral offset) (fromIntegral numGroups) incrementsP groupsPP
    where
    Coeffs {..} = prepareCoeffs n interpolation decimation coeffs

foreign import ccall unsafe "resample2RC"
    resample2RC_c :: CInt -> CInt -> CInt -> CInt -> Ptr CInt -> Ptr (Ptr CFloat) -> Ptr CFloat -> Ptr CFloat -> IO CInt

resampleCRC :: ResampleRC
resampleCRC = mkResamplerC resample2RC_c 1

foreign import ccall unsafe "resampleSSERC"
    resampleCSSERC_c :: CInt -> CInt -> CInt -> CInt -> Ptr CInt -> Ptr (Ptr CFloat) -> Ptr CFloat -> Ptr CFloat -> IO CInt

resampleCSSERC :: ResampleRC
resampleCSSERC = mkResamplerC resampleCSSERC_c 4

foreign import ccall unsafe "resampleAVXRC"
    resampleAVXRC_c :: CInt -> CInt -> CInt -> CInt -> Ptr CInt -> Ptr (Ptr CFloat) -> Ptr CFloat -> Ptr CFloat -> IO CInt

resampleCAVXRC :: ResampleRC
resampleCAVXRC = mkResamplerC resampleAVXRC_c 8

{-
 - Cross buffer
-}

{-# INLINE decimateCrossHighLevel #-}
decimateCrossHighLevel :: (PrimMonad m, Functor m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) => Int -> v b -> Int -> v a -> v a -> vm (PrimState m) a -> m ()
decimateCrossHighLevel factor coeffs num lastBuf nextBuf outBuf = fill x outBuf
    where
    x = VFB.map dotProd (VFB.iterateN num (+ factor) 0)
    dotProd i = VG.sum $ VG.zipWith mult (VG.unsafeDrop i lastBuf VG.++ nextBuf) coeffs

{-# INLINE filterCrossHighLevel #-}
filterCrossHighLevel :: (PrimMonad m, Functor m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) => v b -> Int -> v a -> v a -> vm (PrimState m) a -> m ()
filterCrossHighLevel coeffs num lastBuf nextBuf outBuf = fill (VFB.generate num dotProd) outBuf
    where
    dotProd i = VG.sum $ VG.zipWith mult (VG.unsafeDrop i lastBuf VG.++ nextBuf) coeffs

{-# INLINE resampleCrossHighLevel #-}
resampleCrossHighLevel :: (PrimMonad m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) => Int -> Int -> v b -> Int -> Int -> v a -> v a -> vm (PrimState m) a -> m Int
resampleCrossHighLevel interpolation decimation coeffs filterOffset count lastBuf nextBuf outBuf = fill 0 filterOffset 0
    where
    fill i filterOffset inputOffset
        | i < count = do
            let dp = dotProd filterOffset inputOffset
            VGM.unsafeWrite outBuf i dp
            let (q, r)        = divMod (decimation - filterOffset - 1) interpolation
                inputOffset'  = inputOffset + q + 1
                filterOffset' = interpolation - 1 - r
            filterOffset' `seq` inputOffset' `seq` fill (i + 1) filterOffset' inputOffset'
        | otherwise = return filterOffset
    dotProd filterOffset i = VG.sum $ VG.zipWith mult (VG.unsafeDrop i lastBuf VG.++ nextBuf) (stride interpolation (VG.unsafeDrop filterOffset coeffs))

foreign import ccall unsafe "dcBlocker"
    c_dcBlocker :: CInt -> CFloat -> CFloat -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

dcBlocker :: Int -> Float -> Float -> VS.Vector Float -> VS.MVector RealWorld Float -> IO (Float, Float)
dcBlocker num lastSample lastOutput inBuf outBuf =
    alloca $ \fsp ->
        alloca $ \fop ->
            VS.unsafeWith (coerce inBuf) $ \iPtr ->
                VSM.unsafeWith (coerce outBuf) $ \oPtr -> do
                    c_dcBlocker (fromIntegral num) (realToFrac lastSample) (realToFrac lastOutput) fsp fop iPtr oPtr
                    r1 <- peek fsp
                    r2 <- peek fop
                    return (realToFrac r1, realToFrac r2)

max_stars_repo_name: adamwalker/sdr

max_issues_repo_name: peixian/sdr

max_forks_repo_name: adamwalker/sdr