text:
{-# LANGUAGE ForeignFunctionInterface, GeneralizedNewtypeDeriving #-}

module Numerical.HBLAS.BLAS.FFI.Level3  where

import Foreign.Ptr
import Foreign()
import Foreign.C.Types
import Data.Complex
import Numerical.HBLAS.BLAS.FFI

--------------------------------------------------------------------------------
------------------------------ | BLAS LEVEL 3 ROUTINES
--------------------------------------------------------------------------------
-----------------------  |  Level 3 ops are faster than Levels 1 or 2
--------------------------------------------------------------------------------

--void cblas_sgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, const blasint M, const blasint N, const blasint K,
--         const float alpha, const float *A, const blasint lda, const float *B, const blasint ldb, const float beta, float *C, const blasint ldc);

--void cblas_dgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, const blasint M, const blasint N, const blasint K,
--         const double alpha, const double *A, const blasint lda, const double *B, const blasint ldb, const double beta, double *C, const blasint ldc);
--void cblas_cgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, const blasint M, const blasint N, const blasint K,
--         const float *alpha, const float *A, const blasint lda, const float *B, const blasint ldb, const float *beta, float *C, const blasint ldc);
--void cblas_zgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, const blasint M, const blasint N, const blasint K,
--         const double *alpha, const double *A, const blasint lda, const double *B, const blasint ldb, const double *beta, double *C, const blasint ldc);

-- |  Matrix mult for general dense matrices
type GemmFunFFI scale el = CBLAS_ORDERT ->   CBLAS_TRANSPOSET -> CBLAS_TRANSPOSET->
        CInt -> CInt -> CInt -> {- scal A * B -} scale  -> {- Matrix A-} Ptr el  -> CInt -> {- B -}  Ptr el -> CInt->
            scale -> {- C -}  Ptr el -> CInt -> IO ()

{- C := alpha*op( A )*op( B ) + beta*C ,  -}

-- matrix mult!
foreign import ccall unsafe "cblas_sgemm"
    cblas_sgemm_unsafe :: GemmFunFFI Float Float

foreign import ccall unsafe "cblas_dgemm"
    cblas_dgemm_unsafe :: GemmFunFFI Double Double

foreign import ccall unsafe "cblas_cgemm"
    cblas_cgemm_unsafe :: GemmFunFFI (Ptr(Complex Float)) (Complex Float)

foreign import ccall unsafe "cblas_zgemm"
    cblas_zgemm_unsafe :: GemmFunFFI (Ptr (Complex Double)) (Complex Double)

-- safe ffi variant for large inputs
foreign import ccall "cblas_sgemm"
    cblas_sgemm_safe :: GemmFunFFI Float Float

foreign import ccall "cblas_dgemm"
    cblas_dgemm_safe :: GemmFunFFI Double Double

foreign import ccall "cblas_cgemm"
    cblas_cgemm_safe :: GemmFunFFI (Ptr(Complex Float)) (Complex Float)

foreign import ccall "cblas_zgemm"
    cblas_zgemm_safe :: GemmFunFFI (Ptr (Complex Double)) (Complex Double)

-----------------------------------------
----- |  Matrix mult for Symmetric Matrices
-----------------------------------------

--void cblas_ssymm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side, const enum CBLAS_UPLO Uplo, const blasint M, const blasint N,
--                 const float alpha, const float *A, const blasint lda, const float *B, const blasint ldb, const float beta, float *C, const blasint ldc);
--void cblas_dsymm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side, const enum CBLAS_UPLO Uplo, const blasint M, const blasint N,
--                 const double alpha, const double *A, const blasint lda, const double *B, const blasint ldb, const double beta, double *C, const blasint ldc);
--void cblas_csymm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side, const enum CBLAS_UPLO Uplo, const blasint M, const blasint N,
--                 const float *alpha, const float *A, const blasint lda, const float *B, const blasint ldb, const float *beta, float *C, const blasint ldc);
--void cblas_zsymm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side, const enum CBLAS_UPLO Uplo, const blasint M, const blasint N,
--                 const double *alpha, const double *A, const blasint lda, const double *B, const blasint ldb, const double *beta, double *C, const blasint ldc);

type SymmFunFFI scale el = CBLAS_ORDERT -> CBLAS_SIDET -> CBLAS_UPLOT ->
     CInt->CInt -> scale -> Ptr el -> CInt -> Ptr el -> CInt -> scale ->Ptr el -> CInt -> IO ()

foreign import ccall unsafe "cblas_ssymm"
    cblas_ssymm_unsafe :: SymmFunFFI Float Float

foreign import ccall unsafe "cblas_dsymm"
    cblas_dsymm_unsafe :: SymmFunFFI Double Double

foreign import ccall unsafe "cblas_csymm"
    cblas_csymm_unsafe :: SymmFunFFI (Ptr (Complex Float )) (Complex Float)

foreign import ccall unsafe "cblas_zsymm"
    cblas_zsymm_unsafe :: SymmFunFFI (Ptr (Complex Double)) (Complex Double)

-- safe ffi variant,
foreign import ccall  "cblas_ssymm"
    cblas_ssymm_safe :: SymmFunFFI Float Float

foreign import ccall  "cblas_dsymm"
    cblas_dsymm_safe :: SymmFunFFI Double Double

foreign import ccall  "cblas_csymm"
    cblas_csymm_safe :: SymmFunFFI (Ptr (Complex Float )) (Complex Float)

foreign import ccall  "cblas_zsymm"
    cblas_zsymm_safe :: SymmFunFFI (Ptr (Complex Double)) (Complex Double)

-----------------------------------
--- |  symmetric rank k  matrix update, C := alpha*A*A' + beta*C
--- or C = alpha*A'*A + beta*C
------------------------------------

--void cblas_ssyrk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE Trans,
--         const blasint N, const blasint K, const float alpha, const float *A, const blasint lda, const float beta, float *C, const blasint ldc);
--void cblas_dsyrk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE Trans,
--         const blasint N, const blasint K, const double alpha, const double *A, const blasint lda, const double beta, double *C, const blasint ldc);
--void cblas_csyrk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE Trans,
--         const blasint N, const blasint K, const float *alpha, const float *A, const blasint lda, const float *beta, float *C, const blasint ldc);
--void cblas_zsyrk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE Trans,
--         const blasint N, const blasint K, const double *alpha, const double *A, const blasint lda, const double *beta, double *C, const blasint ldc);

type SyrkFunFFI scale el = CBLAS_ORDERT -> CBLAS_UPLOT -> CBLAS_TRANSPOSET ->
     CInt -> CInt  -> scale -> Ptr el -> CInt -> scale -> Ptr el -> CInt -> IO ()
foreign import ccall unsafe "cblas_ssyrk"
    cblas_ssyrk_unsafe :: SyrkFunFFI Float Float
foreign import ccall unsafe "cblas_dsyrk"
    cblas_dsyrk_unsafe :: SyrkFunFFI Double Double
foreign import ccall unsafe "cblas_csyrk"
    cblas_csyrk_unsafe :: SyrkFunFFI (Ptr(Complex Float)) (Complex Float)
foreign import ccall unsafe "cblas_zsyrk"
    cblas_zsyrk_unsafe :: SyrkFunFFI (Ptr(Complex Double)) (Complex Double)

foreign import ccall safe "cblas_ssyrk"
    cblas_ssyrk_safe :: SyrkFunFFI Float Float
foreign import ccall safe "cblas_dsyrk"
    cblas_dsyrk_safe :: SyrkFunFFI Double Double
foreign import ccall safe "cblas_csyrk"
    cblas_csyrk_safe :: SyrkFunFFI (Ptr(Complex Float)) (Complex Float)
foreign import ccall safe "cblas_zsyrk"
    cblas_zsyrk_safe :: SyrkFunFFI (Ptr(Complex Double)) (Complex Double)
----------------------
----- | Symmetric Rank 2k matrix update, C= alpha* A*B' + alpha* B*A' + beta * C
----- or C= alpha* A'*B + alpha* B'*A + beta * C
-------------------

--void cblas_ssyr2k(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE Trans,
--          const blasint N, const blasint K, const float alpha, const float *A, const blasint lda, const float *B, const blasint ldb, const float beta, float *C, const blasint ldc);
--void cblas_dsyr2k(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE Trans,
--          const blasint N, const blasint K, const double alpha, const double *A, const blasint lda, const double *B, const blasint ldb, const double beta, double *C, const blasint ldc);
--void cblas_csyr2k(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE Trans,
--          const blasint N, const blasint K, const float *alpha, const float *A, const blasint lda, const float *B, const blasint ldb, const float *beta, float *C, const blasint ldc);
--void cblas_zsyr2k(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE Trans,
         --const blasint N, const blasint K, const double *alpha, const double *A, const blasint lda, const double *B, const blasint ldb, const double *beta, double *C, const blasint ldc);

type Syr2kFunFFI scale el = CBLAS_ORDERT -> CBLAS_UPLOT -> CBLAS_TRANSPOSET  ->
     CInt->CInt -> scale -> Ptr el -> CInt -> Ptr el -> CInt ->
     scale ->Ptr el -> CInt -> IO ()

foreign  import ccall unsafe "cblas_ssyr2k"
    cblas_ssyr2k_unsafe :: Syr2kFunFFI Float Float
foreign import ccall unsafe "cblas_dsyr2k"
    cblas_dsyr2k_unsafe :: Syr2kFunFFI Double Double
foreign  import ccall unsafe "cblas_csyr2k"
    cblas_csyr2k_unsafe :: Syr2kFunFFI (Ptr (Complex Float)) (Complex Float)
foreign  import ccall unsafe "cblas_zsyr2k"
    cblas_zsyr2k_unsafe :: Syr2kFunFFI (Ptr (Complex Double)) (Complex Double)

foreign  import ccall safe "cblas_ssyr2k"
    cblas_ssyr2k_safe :: Syr2kFunFFI Float Float
foreign import ccall safe "cblas_dsyr2k"
    cblas_dsyr2k_safe :: Syr2kFunFFI Double Double
foreign  import ccall safe "cblas_csyr2k"
    cblas_csyr2k_safe :: Syr2kFunFFI (Ptr (Complex Float)) (Complex Float)
foreign  import ccall safe "cblas_zsyr2k"
    cblas_zsyr2k_safe :: Syr2kFunFFI (Ptr (Complex Double)) (Complex Double)

-------------------------------
--------  |  matrix matrix product for triangular matrices
------------------------------

type TrmmFunFFI scale el = CBLAS_ORDERT -> CBLAS_SIDET -> CBLAS_UPLOT -> CBLAS_TRANSPOSET -> CBLAS_DIAGT ->
     CInt -> CInt -> scale -> Ptr el -> CInt -> Ptr el -> CInt -> IO ()
foreign  import ccall unsafe "cblas_strmm"
    cblas_strmm_unsafe :: TrmmFunFFI Float Float
foreign  import ccall unsafe "cblas_dtrmm"
    cblas_dtrmm_unsafe :: TrmmFunFFI Double Double
foreign  import ccall unsafe "cblas_ctrmm"
    cblas_ctrmm_unsafe :: TrmmFunFFI (Ptr (Complex Float )) (Complex Float)
foreign  import ccall unsafe "cblas_ztrmm"
    cblas_ztrmm_unsafe :: TrmmFunFFI (Ptr (Complex Double )) (Complex Double)

foreign  import ccall safe "cblas_strmm"
    cblas_strmm_safe :: TrmmFunFFI Float Float
foreign  import ccall safe "cblas_dtrmm"
    cblas_dtrmm_safe :: TrmmFunFFI Double Double
foreign  import ccall safe "cblas_ctrmm"
    cblas_ctrmm_safe :: TrmmFunFFI (Ptr (Complex Float )) (Complex Float)
foreign  import ccall safe "cblas_ztrmm"
    cblas_ztrmm_safe :: TrmmFunFFI (Ptr (Complex Double )) (Complex Double)
--void cblas_strmm(  enum CBLAS_ORDER Order,   enum CBLAS_SIDE Side,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,
--                   enum CBLAS_DIAG Diag,   CInt M,   CInt N,   Float alpha,   Float *A,   CInt lda, Float *B,   CInt ldb);
--void cblas_dtrmm(  enum CBLAS_ORDER Order,   enum CBLAS_SIDE Side,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,
--                   enum CBLAS_DIAG Diag,   CInt M,   CInt N,   Double alpha,   Double *A,   CInt lda, Double *B,   CInt ldb);
--void cblas_ctrmm(  enum CBLAS_ORDER Order,   enum CBLAS_SIDE Side,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,
--                   enum CBLAS_DIAG Diag,   CInt M,   CInt N,   Float *alpha,   Float *A,   CInt lda, Float *B,   CInt ldb);
--void cblas_ztrmm(  enum CBLAS_ORDER Order,   enum CBLAS_SIDE Side,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,
--                   enum CBLAS_DIAG Diag,   CInt M,   CInt N,   Double *alpha,   Double *A,   CInt lda, Double *B,   CInt ldb);

------------------------
--  |  triangular solvers
-----------------------

--
--TRSM solves  op(A)*X = alpha*B or  X*op(A) = alpha*B
--op(A) is one of op(A) = A, or op(A) = A', or op(A) = conjg(A').
-- A is a unit, or non-unit, upper or lower triangular matrix
----
type TrsmFunFFI scale el = CBLAS_ORDERT -> CBLAS_SIDET -> CBLAS_UPLOT -> CBLAS_TRANSPOSET -> CBLAS_DIAGT ->
     CInt -> CInt -> scale -> Ptr el -> CInt -> Ptr el -> CInt -> IO ()
foreign  import ccall unsafe "cblas_strsm"
    cblas_strsm_unsafe :: TrsmFunFFI Float Float
foreign  import ccall unsafe "cblas_dtrsm"
    cblas_dtrsm_unsafe :: TrsmFunFFI Double Double
foreign  import ccall unsafe "cblas_ctrsm"
    cblas_ctrsm_unsafe :: TrsmFunFFI (Ptr (Complex Float )) (Complex Float)
foreign  import ccall unsafe "cblas_ztrsm"
    cblas_ztrsm_unsafe :: TrsmFunFFI (Ptr (Complex Double )) (Complex Double)

foreign  import ccall safe "cblas_strsm"
    cblas_strsm_safe :: TrsmFunFFI Float Float
foreign  import ccall safe "cblas_dtrsm"
    cblas_dtrsm_safe :: TrsmFunFFI Double Double
foreign  import ccall safe "cblas_ctrsm"
    cblas_ctrsm_safe :: TrsmFunFFI (Ptr (Complex Float )) (Complex Float)
foreign  import ccall safe "cblas_ztrsm"
    cblas_ztrsm_safe :: TrsmFunFFI (Ptr (Complex Double )) (Complex Double)
--void cblas_strsm(  enum CBLAS_ORDER Order,   enum CBLAS_SIDE Side,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,
--                   enum CBLAS_DIAG Diag,   CInt M,   CInt N,   Float alpha,   Float *A,   CInt lda, Float *B,   CInt ldb);
--void cblas_dtrsm(  enum CBLAS_ORDER Order,   enum CBLAS_SIDE Side,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,
--                   enum CBLAS_DIAG Diag,   CInt M,   CInt N,   Double alpha,   Double *A,   CInt lda, Double *B,   CInt ldb);
--void cblas_ctrsm(  enum CBLAS_ORDER Order,   enum CBLAS_SIDE Side,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,
--                   enum CBLAS_DIAG Diag,   CInt M,   CInt N,   Float *alpha,   Float *A,   CInt lda, Float *B,   CInt ldb);
--void cblas_ztrsm(  enum CBLAS_ORDER Order,   enum CBLAS_SIDE Side,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,
--                   enum CBLAS_DIAG Diag,   CInt M,   CInt N,   Double *alpha,   Double *A,   CInt lda, Double *B,   CInt ldb);

-------------------------
-- | hermitian matrix mult
------------------------

type HemmFunFFI  el = CBLAS_ORDERT -> CBLAS_SIDET -> CBLAS_UPLOT ->
     CInt -> CInt -> Ptr el -> Ptr el -> CInt -> Ptr el -> CInt -> Ptr el -> Ptr el -> CInt -> IO ()

foreign  import ccall unsafe "cblas_chemm"
    cblas_chemm_unsafe :: HemmFunFFI (Complex Float)
foreign  import ccall unsafe "cblas_zhemm"
    cblas_zhemm_unsafe :: HemmFunFFI  (Complex Double)

foreign  import ccall safe "cblas_chemm"
    cblas_chemm_safe :: HemmFunFFI (Complex Float)
foreign  import ccall safe "cblas_zhemm"
    cblas_zhemm_safe :: HemmFunFFI  (Complex Double)

--void cblas_chemm(  enum CBLAS_ORDER Order,   enum CBLAS_SIDE Side,   enum CBLAS_UPLO Uplo,   CInt M,   CInt N,
--                   Float *alpha,   Float *A,   CInt lda,   Float *B,   CInt ldb,   Float *beta, Float *C,   CInt ldc);
--void cblas_zhemm(  enum CBLAS_ORDER Order,   enum CBLAS_SIDE Side,   enum CBLAS_UPLO Uplo,   CInt M,   CInt N,
--                   Double *alpha,   Double *A,   CInt lda,   Double *B,   CInt ldb,   Double *beta, Double *C,   CInt ldc);

type HerkFunFFI scale el = CBLAS_ORDERT -> CBLAS_UPLOT -> CBLAS_TRANSPOSET ->
     CInt -> CInt -> scale -> Ptr el -> CInt -> scale -> Ptr el -> CInt -> IO ()

foreign  import ccall unsafe "cblas_cherk"
    cblas_cherk_unsafe :: HerkFunFFI Float (Complex Float)
foreign  import ccall unsafe "cblas_zherk"
    cblas_zherk_unsafe :: HerkFunFFI Double (Complex Double)

foreign  import ccall safe "cblas_cherk"
    cblas_cherk_safe :: HerkFunFFI Float (Complex Float)
foreign  import ccall safe "cblas_zherk"
    cblas_zherk_safe :: HerkFunFFI Double (Complex Double)
--void cblas_cherk(  enum CBLAS_ORDER Order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE Trans,   CInt N,   CInt K,
--                   Float alpha,   Float *A,   CInt lda,   Float beta, Float *C,   CInt ldc);
--void cblas_zherk(  enum CBLAS_ORDER Order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE Trans,   CInt N,   CInt K,
--                   Double alpha,   Double *A,   CInt lda,   Double beta, Double *C,   CInt ldc);

type Her2kFunFFI scale el = CBLAS_ORDERT -> CBLAS_UPLOT -> CBLAS_TRANSPOSET ->
     CInt -> CInt -> Ptr el  -> Ptr el -> CInt -> Ptr el -> CInt -> scale ->Ptr el -> CInt -> IO ()

foreign  import ccall unsafe "cblas_cher2k"
    cblas_cher2k_unsafe :: Her2kFunFFI  Float  (Complex Float)
foreign  import ccall unsafe "cblas_zher2k"
    cblas_zher2k_unsafe :: Her2kFunFFI  Double  (Complex Double)

foreign  import ccall safe "cblas_cher2k"
    cblas_cher2k_safe :: Her2kFunFFI Float (Complex Float)
foreign  import ccall safe "cblas_zher2k"
    cblas_zher2k_safe :: Her2kFunFFI Double (Complex Double)
--void cblas_cher2k(  enum CBLAS_ORDER Order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE Trans,   CInt N,   CInt K,
--                    Float *alpha,   Float *A,   CInt lda,   Float *B,   CInt ldb,   Float beta, Float *C,   CInt ldc);
--void cblas_zher2k(  enum CBLAS_ORDER Order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE Trans,   CInt N,   CInt K,
--                    Double *alpha,   Double *A,   CInt lda,   Double *B,   CInt ldb,   Double beta, Double *C,   CInt ldc);

----void cblas_xerbla(CInt p, char *rout, char *form, ...);

max_stars_repo_name: schnecki/hblas

max_issues_repo_name: schnecki/hblas

max_forks_repo_name: schnecki/hblas


text:
module Numeric.FFT.Types
       ( VCD, MVCD, VVCD, VMVCD, VVVCD, VI
       , Direction (..), Plan (..), BaseTransform (..)
       ) where

import           Data.Complex
import qualified Data.Vector                 as V
import           Data.Vector.Unboxed
import qualified Data.Vector.Unboxed.Mutable as MV

-- | Some useful type synonyms.
type VCD = Vector (Complex Double)
type MVCD s = MV.MVector s (Complex Double)
type VVCD = V.Vector VCD
type VVVCD = V.Vector VVCD
type VMVCD a = V.Vector (MVCD a)
type VI = Vector Int

-- | Transform direction: 'Forward' is the normal FFT, 'Inverse' is
-- inverse FFT.
data Direction = Forward | Inverse deriving (Eq, Show)

-- | A FFT plan.  This depends only on the problem size and can be
-- pre-computed and reused to transform (and inverse transform) any
-- number of vectors of the given size.
data Plan = Plan { plDLInfo  :: V.Vector (Int, Int, VVVCD, VVVCD)
                   -- ^ Size information and diagonal matrix entries
                   -- for Danielson-Lanczos recursive decomposition of
                   -- problem size.
                 , plPermute :: Maybe VI
                   -- ^ Input vector permutation to use before base
                   -- transformation and recursive Danielson-Lanczos
                   -- composition.
                 , plBase    :: BaseTransform
                   -- ^ Base transformation used for each sub-vector
                   -- before performing recursive Danielson-Lanczos
                   -- steps to form the full FFT result.
                 } deriving (Eq, Show)

-- | A "base transform" used at the "bottom" of the recursive
-- Cooley-Tukey decomposition of the input problem size: either a
-- simple DFT, a special hard-coded small problem size case, or a
-- Rader prime-length FFT invocation.
data BaseTransform = SpecialBase { baseSize :: Int }
                     -- ^ Hard-coded small-size base transform.
                   | DFTBase { baseSize :: Int
                             , dftWsFwd :: VCD
                             , dftWsInv :: VCD }
                     -- ^ Simple DFT base transform, giving problem
                     -- size and powers of roots of unity needed for
                     -- transform.
                   | RaderBase { baseSize      :: Int
                               , raderOutPerm  :: VI
                               , raderBFwd     :: VCD
                               , raderBInv     :: VCD
                               , raderConvSize :: Int
                               , raderConvPlan :: Plan }
                     -- ^ Prime-length Rader FFT base transform,
                     -- giving problem size, output index permutation
                     -- (the input index permutation is folded into
                     -- the main input permutation of the full
                     -- transform), pre-transformed Rader b sequence
                     -- for forward and inverse problems, the (padded
                     -- or not) problem size for Rader sequence
                     -- convolution and a sub-plan (either of size
                     -- baseSize-1 or the next larger power of two)
                     -- for computing the Rader convolution.
                     deriving (Eq, Show)

ext: hs

max_stars_repo_name: ian-ross/arb-fft

max_issues_repo_name: ian-ross/arb-fft

max_forks_repo_name: ian-ross/arb-fft


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE OverloadedLists #-}

module Main where

import qualified Data.Vector.Unboxed         as U
import           Data.Vector.Unboxed         (Vector)
import qualified Data.Vector as V

-- DLA
import qualified Statistics.Matrix as M
import qualified Statistics.Matrix.Fast as MF
import qualified Statistics.Matrix.Fast.Algorithms as A

-- hmatrix
import qualified Numeric.LinearAlgebra as H

-- numhask
import qualified NumHask.Array as NH
import qualified NumHask.Prelude as NP

-- data.matrix from matrix
import qualified Data.Matrix as DMX

-- massiv
import qualified Data.Massiv.Array as MA

import qualified System.Random.MWC as Mwc

import qualified Weigh as W

n :: Int
n = 100

vectorGen :: IO (Vector Double)
vectorGen =  do
    gen <- Mwc.create
    Mwc.uniformVector gen (n*n)

matrixDLA :: IO M.Matrix
matrixDLA = do
    vec <- vectorGen
    return $ M.Matrix n n vec

matrixH :: IO (H.Matrix H.R)
matrixH = do
    vec <- vectorGen
    return $ (n H.>< n) $ U.toList $ vec

identH :: Int -> H.Matrix Double
identH = H.ident

identDMX :: Int -> DMX.Matrix Double
identDMX = DMX.identity

elemZero :: Double -> Double
elemZero = const 0

elemSqr :: Double -> Double
elemSqr x = x*x

mapH :: (Double -> Double) -> H.Matrix Double -> H.Matrix Double
mapH = H.cmap

main :: IO ()
main = do

    let

    vDLA <- vectorGen
    uDLA <- vectorGen

    let

    --
      subDLA = U.take n vDLA
      aDLA = M.Matrix n n vDLA
      bDLA = M.Matrix n n uDLA

    --
      vList = U.toList vDLA
      uList = U.toList uDLA

    --
      aH = (n H.>< n) vList
      bH = (n H.>< n) uList

      subH = H.fromList . take n $ vList
      vH = H.fromList vList

    --
      aNH = NP.fromList vList :: NH.Array V.Vector '[50, 50] Double
      bNH = NP.fromList uList :: NH.Array V.Vector '[50, 50] Double

      subNH = NP.fromList . take n $ vList :: NH.Array V.Vector '[50] Double
      vNH = NP.fromList vList :: NH.Array V.Vector '[2500] Double

    --
      aDMX = DMX.fromList n n vList
      bDMX = DMX.fromList n n uList

    --
      vMA = MA.fromList MA.Seq vList :: MA.Array MA.P MA.Ix1 Double
      aMA = MA.resize' (MA.Sz (n MA.:. n)) vMA :: MA.Array MA.P MA.Ix2 Double
      bMA = MA.resize' (MA.Sz (n MA.:. n)) $ MA.fromList MA.Seq uList :: MA.Array MA.P MA.Ix2 Double

    W.mainWith (do
               W.func "DLA - multiplication" (MF.multiply aDLA) bDLA
               W.func "DLA - qr factorization" A.qr aDLA
               W.func "DLA - transpose" M.transpose aDLA
               W.func "DLA - norm" MF.norm vDLA
               W.func "DLA - row" (M.row  aDLA) 0
               W.func "DLA - column" (M.column  aDLA) 0
               W.func "DLA - identity" M.ident n

               W.func "Hmatrix - multiplication" ((<>) aH) bH
               W.func "Hmatrix - qr factorization" H.qr aH
               W.func "Hmatrix - transpose" H.tr aH
               W.func "Hmatrix - norm" H.norm_2 vH
               W.func "Hmatrix - row" ((H.?) aH) [0]
               W.func "Hmatrix - column" ((H.¿) aH) [0]
               W.func "Hmatrix - identity" identH n

               W.func "NumHask Array - multiplication" (NH.mmult aNH) bNH
               W.func "NumHask Array - transpose" NH.transpose aNH
               W.func "NumHask Array - norm" (sqrt . (NP.<.> vNH)) vNH
               W.func "NumHask Array - row" (NH.row (NP.Proxy :: NP.Proxy 0)) aNH
               W.func "NumHask Array - column" (NH.col (NP.Proxy :: NP.Proxy 0)) aNH

               W.func "Massiv - norm" (sqrt . MA.foldlS (+) 0 . MA.zipWith (*) vMA) vMA
               W.func "Massiv - multiplication" ((MA.|*|) aMA) bMA
               W.func "Massiv - transpose" (MA.computeAs MA.P . MA.transpose) aMA
               W.func "Massiv - row" (MA.computeAs MA.P . (MA.!>) aMA) 0

               W.func "matrix - multiplication" (DMX.multStrassenMixed aDMX) bDMX
               W.func "matrix - transpose" DMX.transpose aDMX
               W.func "matrix - row" (DMX.getRow 1) aDMX
               W.func "matrix - column" (DMX.getCol 1) aDMX
               W.func "matrix - identity" identDMX n
               )

max_stars_repo_name: lehins/fastest-matrices

max_issues_repo_name: lehins/fastest-matrices

max_forks_repo_name: lehins/fastest-matrices


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE DeriveGeneric #-}
#if MIN_VERSION_base(4,12,0)
{-# LANGUAGE DerivingVia #-}
#endif
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE StandaloneDeriving #-}

import Control.Monad ((>=>),forM)
import Control.Applicative (liftA2)
import Data.Complex
import Data.Either
import Data.Euclidean hiding (gcd)
import Data.Fixed
import Data.Functor.Const.Compat
import Data.Functor.Identity.Compat
import Data.HashMap.Strict
import Data.HashSet
import Data.Hashable
import Data.Int
import Data.IntMap.Strict
import Data.IntSet
import Data.Map
import Data.Monoid
import Data.Ord (Down(..))
import Data.Orphans ()
import Data.Proxy (Proxy(..))
import Data.Ratio
import Data.Semigroup.Compat
import Data.Semiring
import Data.Semiring.Generic
import Data.Sequence
import Data.Set
import Data.Word
import GHC.Generics
import GHC.Natural
import qualified Data.Foldable as F
import Prelude hiding (Num(..),(^))
import System.IO.Unsafe
import Test.QuickCheck hiding (Fixed)
import Test.QuickCheck.Gen (suchThat)
import Test.QuickCheck.Instances ()
import Test.Tasty (defaultMain, testGroup, TestTree)
import Text.Show.Functions ()
import qualified Data.Map as Map
import qualified Data.Sequence as Seq
import qualified GHC.Num as Num
import qualified Test.QuickCheck.Classes as QCC
import qualified Control.Exception as E

main :: IO ()
main = do
  QCC.lawsCheckMany namedTests
  F.fold pow_prop
  F.fold gcd_prop

forceError :: a -> IO (Either E.ErrorCall a)
forceError = E.try . E.evaluate

pow_prop :: [IO ()]
pow_prop = [quickCheck pow_prop1, quickCheck pow_prop2, quickCheck pow_prop3]

pow_prop1 :: Int -> Negative Int -> Property
pow_prop1 x y = ioProperty $ do
  p <- forceError (x ^ getNegative y)
  pure $ isLeft p

pow_prop2 :: Int -> Positive Int -> Property
pow_prop2 x y = ioProperty $ do
  p <- forceError (x ^ getPositive y)
  pure $ isRight p

pow_prop3 :: Int -> Property
pow_prop3 x = ioProperty $ do
  p <- forceError (x ^ 0)
  case p of
    Left e -> pure False
    Right x -> pure $ x == 1

gcd_prop :: [IO ()]
gcd_prop = [quickCheck gcd_prop1, quickCheck gcd_prop2]

gcd_prop1 :: Integer -> Integer -> Property
gcd_prop1 x y = Num.abs (fst (x `gcdExt` y)) === x `gcd` y

gcd_prop2 :: Integer -> Integer -> Property
gcd_prop2 x y = y /= 0 ==> (x * s) `mod` y === g `mod` y
  where
    (g, s) = x `gcdExt` y

type Laws = QCC.Laws

semiringLaws :: (Arbitrary a, Show a, Eq a, Semiring a) => Proxy a -> [Laws]
semiringLaws p = [QCC.semiringLaws p]

ringLaws :: (Arbitrary a, Show a, Eq a, Ring a) => Proxy a -> [Laws]
ringLaws p = [QCC.semiringLaws p, QCC.ringLaws p]

euclideanLaws :: (Arbitrary a, Show a, Eq a, Euclidean a) => Proxy a -> [Laws]
#if MIN_VERSION_quickcheck_classes(0,6,3)
euclideanLaws p = [QCC.gcdDomainLaws p, QCC.euclideanLaws p]
#else
euclideanLaws _ = []
#endif

namedTests :: [(String, [Laws])]
namedTests =
  [ ("Bool", semiringLaws pBool)
--  , ("Double", ringLaws pDouble) -- needs to be within some epsilon
--  , ("Float", ringLaws pFloat)   -- needs to be within some epsilon
--  , ("Complex", ringLaws pComplex) -- needs to be within some epsilon
  , ("Mod2", semiringLaws pMod2 ++ ringLaws pMod2 ++ euclideanLaws pMod2)
  , ("Int", semiringLaws pInt ++ ringLaws pInt ++ euclideanLaws pInt)
  , ("Word", semiringLaws pWord ++ ringLaws pWord ++ euclideanLaws pWord)
  , ("Natural", semiringLaws pNatural ++ euclideanLaws pNatural)
  , ("Rational", semiringLaws pRational ++ ringLaws pRational ++ euclideanLaws pRational)
  , ("()", semiringLaws pInt8 ++ ringLaws pUnit)
  , ("Maybe", semiringLaws pMaybe)
  , ("PosRatio", semiringLaws pPosRatio)
  , ("IO", semiringLaws pIO ++ ringLaws pIO)
  , ("Fixed", ringLaws pFixed)
  , ("Identity", ringLaws pIdentity)
  , ("Dual", ringLaws pDual)
  , ("(->)", ringLaws pFunction)
  , ("Down", ringLaws pDown)
  , ("Const", ringLaws pConst)
--  , ("IntMap", semiringLaws pIntMap) -- needs newtypes
  , ("Set", semiringLaws pSet)
--  , ("IntSet", semiringLaws pIntSet) -- needs newtypes
  , ("HashSet", semiringLaws pHashSet)
  , ("HashMap", semiringLaws pHashMap)
  , ("Map", semiringLaws pMap)
  , ("Predicate", semiringLaws pPredicate)
  , ("Equivalence", semiringLaws pEquivalence)
  , ("Op", semiringLaws pOp ++ ringLaws pOp)
  , ("Ap", semiringLaws pAp ++ ringLaws pAp)

  , ("IntSet Sum", semiringLaws pIntSetSum)
  , ("IntSet Product", semiringLaws pIntSetProduct)
  , ("IntSet Min", semiringLaws pIntSetMin)
  , ("IntSet Max", semiringLaws pIntSetMax)

  , ("IntMap Sum", semiringLaws pIntMapSum)
  , ("IntMap Product", semiringLaws pIntMapProduct)
  , ("IntMap Min", semiringLaws pIntMapMin)
  , ("IntMap Max", semiringLaws pIntMapMax)

  , ("AB", semiringLaws pAB)
  ]

#if !(MIN_VERSION_base(4,12,0))
newtype Ap f a = Ap { getAp :: f a }
  deriving (Eq, Functor, Applicative, Show)
instance (Semiring a, Applicative f) => Semiring (Ap f a) where
  zero = pure zero
  one = pure one
  plus = liftA2 plus
  times = liftA2 times
  fromNatural = pure . fromNatural
instance (Ring a, Applicative f) => Ring (Ap f a) where
  negate = fmap negate
#endif
deriving instance (Arbitrary (f a)) => Arbitrary (Ap f a)

newtype Predicate a = Predicate (a -> Bool)
  deriving (Eq, Show)
deriving instance Semiring (Predicate a)
deriving instance (Arbitrary a, CoArbitrary a) => Arbitrary (Predicate a)

newtype Equivalence a = Equivalence (a -> a -> Bool)
  deriving (Eq, Show)
deriving instance Semiring a => Semiring (Equivalence a)
deriving instance (Arbitrary a, CoArbitrary a) => Arbitrary (Equivalence a)

newtype Op a b = Op (b -> a)
  deriving (Eq, Show)
deriving instance Semiring a => Semiring (Op a b)
deriving instance Ring a => Ring (Op a b)
deriving instance (Arbitrary a, CoArbitrary a, CoArbitrary b) => Arbitrary (Op a b)

pAp = Proxy :: Proxy (Ap Identity Int)
pPredicate = Proxy :: Proxy (Predicate Int)
pEquivalence = Proxy :: Proxy (Equivalence Int)
pOp = Proxy :: Proxy (Op Int Int)

deriving instance Arbitrary a => Arbitrary (Down a)

-- | Rationals constrained to being positive
newtype PosRatio a = PosRatio { getPosRatio :: Ratio a }
  deriving (Eq, Show, Semiring)

instance (Arbitrary a, Integral a, Ord a) => Arbitrary (PosRatio a) where
  arbitrary =
    let nm = suchThat (arbitrary :: Gen a) (> 0)
        dm = suchThat (arbitrary :: Gen a) (> 0)
    in PosRatio <$> liftA2 (%) nm dm

instance (Arbitrary a, Eq b) => Eq (a -> b) where
  (==) f g =
    let x = generate arbitrary :: IO a
        x' = unsafePerformIO x
    in f x' == g x'

instance Show (Endo a) where
  show _ = "<<<Endo>>>"

deriving instance (Arbitrary a, Eq a) => Eq (Endo a)

-- no
instance Show (IO a) where
  show _ = "<<<IO action>>>"

-- no no
instance Eq a => Eq (IO a) where
  x == y = unsafePerformIO $ liftA2 (==) x y

-- this is fine
instance Arbitrary a => Arbitrary (IO a) where
  arbitrary = fmap pure (arbitrary :: Gen a)

deriving instance Hashable a => Hashable (Sum a)

deriving instance Arbitrary (IntSetOf a)
deriving instance Arbitrary v => Arbitrary (IntMapOf k v)

deriving instance Arbitrary Mod2

pBool = Proxy :: Proxy Bool
pDouble = Proxy :: Proxy Double
pFloat = Proxy :: Proxy Float
pInt = Proxy :: Proxy Int
pInt8 = Proxy :: Proxy Int8
pInt16 = Proxy :: Proxy Int16
pInt32 = Proxy :: Proxy Int32
pInt64 = Proxy :: Proxy Int64
pInteger = Proxy :: Proxy Integer
pNatural = Proxy :: Proxy Natural
pWord = Proxy :: Proxy Word
pWord8 = Proxy :: Proxy Word8
pWord16 = Proxy :: Proxy Word16
pWord32 = Proxy :: Proxy Word32
pWord64 = Proxy :: Proxy Word64
pUnit = Proxy :: Proxy ()
pMaybe = Proxy :: Proxy (Maybe Int)
pPosRatio = Proxy :: Proxy (PosRatio Int)
pIO = Proxy :: Proxy (IO Int)
pComplex = Proxy :: Proxy (Complex Double)
pFixed = Proxy :: Proxy (Fixed E0)
pMin = Proxy :: Proxy (Min Int)
pMax = Proxy :: Proxy (Max Int)
pIdentity = Proxy :: Proxy (Identity Int)
pDual = Proxy :: Proxy (Dual Int)
pEndo = Proxy :: Proxy (Endo (Sum Int))
pSum  = Proxy :: Proxy (Sum Int)
pProduct = Proxy :: Proxy (Product Int)
pDown = Proxy :: Proxy (Down Int)
pSeq = Proxy :: Proxy (Seq Int)
pSet = Proxy :: Proxy (Set (Sum Int))
pHashSet = Proxy :: Proxy (HashSet (Sum Int))
pFunction = Proxy :: Proxy (Int -> Int)
pMap = Proxy :: Proxy (Map (Sum Int) Int)
pHashMap = Proxy :: Proxy (HashMap (Sum Int) Int)
pConst = Proxy :: Proxy (Const Int Int)
pAlt = Proxy :: Proxy (Alt Maybe Int)
pRational = Proxy :: Proxy Rational
pMod2 = Proxy :: Proxy Mod2
pAB = Proxy :: Proxy AB

pIntSetSum = Proxy :: Proxy (IntSetOf (Sum Int))
pIntSetProduct = Proxy :: Proxy (IntSetOf (Product Int))
pIntSetMin = Proxy :: Proxy (IntSetOf (Min Int))
pIntSetMax = Proxy :: Proxy (IntSetOf (Max Int))

pIntMapSum = Proxy :: Proxy (IntMapOf (Sum Int) Int)
pIntMapProduct = Proxy :: Proxy (IntMapOf (Product Int) Int)
pIntMapMin = Proxy :: Proxy (IntMapOf (Min Int) Int)
pIntMapMax = Proxy :: Proxy (IntMapOf (Max Int) Int)

newtype A = A Integer
  deriving (Show, Eq, Semiring, Arbitrary)
newtype B = B Integer
  deriving (Show, Eq, Semiring, Arbitrary)

data AB = AB A B
  deriving (Show, Eq, Generic)
#if MIN_VERSION_base(4,12,0)
  deriving (Semiring) via (GenericSemiring AB)
#else
instance Semiring AB where
  zero = gzero
  one = gone
  plus = gplus
  times = gtimes
  fromNatural = gfromNatural
#endif

instance Arbitrary AB where
  arbitrary = AB <$> arbitrary <*> arbitrary
  shrink = genericShrink

max_stars_repo_name: hseg/semirings

max_issues_repo_name: hseg/semirings

max_forks_repo_name: hseg/semirings


text:
module Main where

import qualified Numeric.LinearAlgebra as LA
import qualified MachineLearning as ML
import qualified MachineLearning.Regression as MLR

main = do
  -- Step 1. Data loading
  m <- LA.loadMatrix "linear_regression/data.txt"
  let (x, y) = ML.splitToXY m

  -- Step 2. Feature Normalization
      muSigma = ML.meanStddev x
      xNorm = ML.featureNormalization muSigma x
      x1 = ML.addBiasDimension xNorm

  -- Step 3. Learning using 3 methods:
  -- NE: Normal Equation (exact solution)
  -- GD: Gradient Descent (basic iterative method)
  -- BFGS: BFGS (most advanced iterative method)
      zeroTheta = LA.konst 0 (LA.cols x1)
      thetaNE = MLR.normalEquation x1 y
      (thetaGD, optPathGD) = MLR.minimize (MLR.GradientDescent 0.01) MLR.LeastSquares 0.0001 5000 MLR.RegNone x1 y zeroTheta
      (thetaBFGS, optPathBFGS) = MLR.minimize (MLR.BFGS2 0.1 0.1)    MLR.LeastSquares 0.0001 1500 MLR.RegNone x1 y zeroTheta

  -- Step 4. Prediction
      xPredict = LA.matrix 2 [1650, 3]
      xPredict1 = ML.addBiasDimension $ ML.featureNormalization muSigma xPredict
      yPredictNE = MLR.hypothesis MLR.LeastSquares xPredict1 thetaNE
      yPredictGD = MLR.hypothesis MLR.LeastSquares xPredict1 thetaGD
      yPredictBFGS = MLR.hypothesis MLR.LeastSquares xPredict1 thetaBFGS

  -- Step 5. Printing results
  putStrLn "\n=== Linear Regression Sample Application ===\n"
  putStrLn "Feature normalization, mu:"
  LA.disp 4 $ fst muSigma
  putStrLn "Feature normalization, sigma:"
  LA.disp 4 $ snd muSigma
  putStrLn "Normilized features:"
  LA.dispShort 10 10 4 $ x LA.||| xNorm

  putStrLn ""
  putStrLn $ "Theta (Normal Equation):  " ++ (show thetaNE)
  putStrLn $ "Theta (Gradient Descent): " ++ (show thetaGD)
  putStrLn $ "Theta (BFGS):             " ++ (show thetaBFGS)

  putStrLn ""
  putStrLn $ "Prediction (Normal Equation):  " ++ (show yPredictNE)
  putStrLn $ "Prediction (Gradient Descent): " ++ (show yPredictGD)
  putStrLn $ "Prediction (BFGS):             " ++ (show yPredictBFGS)

  putStrLn ""
  putStrLn "Optimization Path (Gradient Descent):"
  LA.dispShort 20 10 4 optPathGD
  putStrLn "Optimization Path (BFGS2):"
  LA.disp 3 optPathBFGS

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
{-# LANGUAGE FlexibleContexts, UndecidableInstances #-}

module School.Unit.UnitParams
( UnitParams(..) ) where

import Numeric.LinearAlgebra (Container, Matrix, Vector)

data UnitParams a =
   AffineParams { affineWeights :: Matrix a
                , affineBias :: Vector a
                }
 | EmptyParams deriving (Show)

instance (Container Vector a, Eq a, Num a) => Eq (UnitParams a) where
  AffineParams { affineBias = b1, affineWeights = w1 } == AffineParams { affineBias = b2, affineWeights = w2 } =
    b1 == b2 && w1 == w2
  EmptyParams == EmptyParams = True
  _ == _ = False

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


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

import Lib
import qualified Data.Text.IO as T
import qualified MachineLearning.Types as T
import qualified Numeric.LinearAlgebra as LA
import qualified MachineLearning as ML
import qualified MachineLearning.Classification.Binary as BC

main = do
    -- Step 1. Data loading.
    designMatrix <- LA.loadMatrix "data.txt"
    let (x, y) = ML.splitToXY designMatrix
    -- Step 2. Feature normalization (skipped - we don't need feature normalization if we use BFGS2).
    -- See Linear Regression sample app for dedails.

    -- Step 3. Feature mapping.
        x1 = ML.addBiasDimension $ ML.mapFeatures 6 x

    -- Step 4. Learning
        zeroTheta = LA.konst 0 (LA.cols x1)
        (theta, _) = BC.learn (BC.BFGS2 0.1 0.1) 0.000002 1500 (BC.L2 1) x1 y zeroTheta

    -- Step 5. Prediction and checking accuracy
        accuracy = calcAccuracy x1 y theta

    -- Step 6. Printing results.
    putStrLn "\n=== Logistic Regression Sample Application ===\n"

    putStrLn ""
    putStrLn $ "Theta:             " ++ (show theta)

    putStrLn ""
    putStrLn $ "Accuracy on training set data (%): " ++ show (accuracy*100)

    putStrLn "This is me"
    someFunc'
    T.putStrLn "someFunc"
    putStrLn $ anotherFunc 5


text:
module Fractals (  iteratingFractal, mandelbrot, burningShip, tricorn, julia, flipForJulia, Formula, Parameters ) where

import Data.Complex ( Complex(..), conjugate, imagPart, magnitude, realPart )
import Data.Bool ( bool )
import Data.Tuple ( swap )

type Formula = Complex Float -> Complex Float -> Complex Float -- e.g z' = z^2 + c
type Parameters = (Complex Float, Complex Float) -- (z,c)

maxIter = 100

flipForJulia :: Parameters -> Bool -> Parameters
flipForJulia =  bool <*> swap

mandelbrot :: Formula
mandelbrot z c = z ^ 2 + c

burningShip :: Formula
burningShip z c = (abs (realPart z) :+ abs (imagPart z)) ^ 2 + c

tricorn :: Formula
tricorn z c = conjugate z ^ 2 + c

-- i think julia sets can be just our regular sets (mandelbrot for the classic ones) given a const c in reversed c z paremeter order.
-- for now let it be just a regular formula. for which we just reverse z and c
julia :: Formula
julia z _ = z ^ 2 + ((-0.835) :+ (-0.2321))

-- returns intensity in correspondence to a point on the complex plane
iteratingFractal :: Formula -> Parameters -> Int -> Int
iteratingFractal formula (z, c) iter
  | iter > maxIter = 0
  | otherwise = let z' = formula z c in
      if magnitude z' > 2 then iter
      else iteratingFractal formula (z', c) (iter + 1)


text:
module NeuralNet.NetSpec (netSpec) where

import Test.Hspec
import Numeric.LinearAlgebra
import NeuralNet.Net
import NeuralNet.Activation
import NeuralNet.Layer
import NeuralNet.Example
import Control.Exception (evaluate)
import Data.List

matrixSize :: Matrix Double -> (Int, Int)
matrixSize m = (rows m, cols m)

allUnique :: Matrix Double -> Bool
allUnique m = length (nub list) == length list
  where list = concat (toLists m)

createdLayer :: Activation -> (Int, Int) -> NeuronLayer -> Bool
createdLayer a s l = activationEq && wSizeEq && bCorrect && allUnique m
  where
    m = layerW l
    bCorrect = layerB l == col (replicate (fst s) 0)
    activationEq = layerActivation l == a
    wSizeEq = matrixSize m == s

netSpec :: SpecWith ()
netSpec =
  describe "NeuralNet.Net" $ do
    describe "initNN" $ do
      it "returns correctly for small definition" $
        let
          nn = initNN [1..] (8, [LayerDefinition ReLU 5, LayerDefinition Tanh 1, LayerDefinition Sigmoid 2])
          layers = nnLayers nn
        in do
          layers!!0 `shouldSatisfy` createdLayer ReLU (5, 8)
          layers!!1 `shouldSatisfy` createdLayer Tanh (1, 5)
          layers!!2 `shouldSatisfy` createdLayer Sigmoid (2, 1)

      it "returns correctly for single layer definition" $
        let
          nn = initNN [1..] (6, [LayerDefinition ReLU 5])
          layers = nnLayers nn
        in
          layers!!0 `shouldSatisfy` createdLayer ReLU (5, 6)

      it "throws exception for no layers" $
        evaluate (initNN [1..] (6, [])) `shouldThrow` anyException

      it "throws exception for no inputs" $
        evaluate (initNN [1..] (0, [LayerDefinition ReLU 5])) `shouldThrow` anyException

    describe "buildNNFromList" $ do
      it "creates log reg correctly" $
        let
          nn = buildNNFromList (5, [LayerDefinition ReLU 1]) [1, 2, 3, 4, 5, 6]
          layers = nnLayers nn
        in do
          length layers `shouldBe` 1
          layerW (layers!!0) `shouldBe` fromLists [[1, 2, 3, 4, 5]]
          layerB (layers!!0) `shouldBe` fromLists [[6]]

      it "creates multi output correctly" $
        let
          nn = buildNNFromList (3, [LayerDefinition ReLU 2]) [1, 2, 3, 4, 5, 6, 7, 8]
          layers = nnLayers nn
        in do
          length layers `shouldBe` 1
          layerW (layers!!0) `shouldBe` fromLists [[1, 2, 3], [4, 5, 6]]
          layerB (layers!!0) `shouldBe` fromLists [[7], [8]]

      it "creates multi layer correctly" $
        let
          nn = buildNNFromList (2, [LayerDefinition ReLU 2, LayerDefinition Tanh 1]) [1, 2, 3, 4, 5, 6, 7, 8, 9]
          layers = nnLayers nn
        in do
          length layers `shouldBe` 2
          layerW (layers!!0) `shouldBe` fromLists [[1, 2], [3, 4]]
          layerB (layers!!0) `shouldBe` fromLists [[5], [6]]
          layerW (layers!!1) `shouldBe` fromLists [[7, 8]]
          layerB (layers!!1) `shouldBe` fromLists [[9]]

      it "throw error for incorrect list size" $
        evaluate (buildNNFromList (5, [LayerDefinition ReLU 1]) [1, 2]) `shouldThrow` anyException

    describe "isExampleSetCompatibleWithNN" $ do
      it "calculates correctly for compatible" $
        let
          nn = buildNNFromList (3, [LayerDefinition ReLU 1]) [1, 2, 3, 4]
          examples = createExampleSet [([1, 2, 3], 0), ([4, 5, 6], 1)]
        in isExampleSetCompatibleWithNN examples nn `shouldBe` True

      it "calculates correctly for incompatible" $
        let
          nn = buildNNFromList (3, [LayerDefinition ReLU 1]) [1, 2, 3, 4]
          examples = createExampleSet [([1, 2], 0), ([4, 5], 1)]
        in isExampleSetCompatibleWithNN examples nn `shouldBe` False

    describe "nnForward" $ do
      -- Not working, can't figure out why. When change to shouldBe [] it does throw an error
--      it "throws errors when incorrect num inputs provided" $
--        let nn = buildNNFromList (3, [LayerDefinition ReLU 1]) [1, 2, 3, 4]
--        in evaluate (nnForward nn [1, 1]) `shouldThrow` anyErrorCall

      it "calculates correctly for logreg" $
        let nn = buildNNFromList (3, [LayerDefinition ReLU 1]) [1, 2, 3, 4]
        in nnForward nn [1, 1, 2] `shouldBe` [13.0]

      it "calculates correctly for multi layer" $
        let
          def = (3, [LayerDefinition Sigmoid 2, LayerDefinition ReLU 1])
          nums = [1, -1, -1
            , 1, 1, -1
            , 1, -1
            , 2, 3,
            0]
          nn = buildNNFromList def nums
        in
          nnForward nn [1, 1, 1] `shouldBe` [2.5]

    describe "nnForwardSet" $ do
      it "calculates correctly for logreg" $
        let
          nn = buildNNFromList (3, [LayerDefinition ReLU 1]) [1, 2, 3, 4]
          examples = createExampleSet [([1, 1, 2], 0)]
          resultAs = map forwardPropA (nnForwardSet nn examples)
        in resultAs `shouldBe` [fromLists [[1], [1], [2]], scalar 13.0]

      it "calculates correctly for multiple examples, multiple output" $
        let
          nn = buildNNFromList (3, [LayerDefinition ReLU 2, LayerDefinition ReLU 1]) [1, 4, 2, 5, 3, 6, 7, 8, 1, 1, 0]
          examples = createExampleSet [([1, 1, 2], 0), ([1, 1, 0], 0)]
          resultAs = map forwardPropA (nnForwardSet nn examples)
        in resultAs `shouldBe` [fromLists [[1, 1], [1, 1], [2, 0]], fromLists [[16, 12], [28, 16]], fromLists [[44, 28]]]

max_stars_repo_name: danielholmes/neural-net

max_issues_repo_name: danielholmes/neural-net

max_forks_repo_name: danielholmes/neural-net


text:
module Statistics.Information.Utils.KNN where

import qualified Data.KdTree.Static as KdTree
import Data.Matrix
import Data.Metric.Class
import Data.Metric.Vector.Real
import qualified Data.Vector as Vector
import Numeric.SpecFunctions
import Statistics.Information.Utils.Matrix

chebyshevDistance u v =
  distance (Chebyshev $ Vector.fromList u) (Chebyshev $ Vector.fromList v)

kdTree metric xs = KdTree.buildWithDist id metric xs

kdDistances tree metric xs k =
  map
    (\row -> let neighbor = KdTree.kNearest tree (k+1) row !! k in
      metric neighbor row)
    xs

knnDistances :: Real a => ([a] -> [a] -> a) -> [[a]] -> Int -> [a]
knnDistances metric xs k = kdDistances tree metric xs k where
  tree = kdTree metric xs

avgdigamma metric xs distances =
  mean . map (digamma . fromIntegral . length) $ radiusPoints
  where
    tree = kdTree metric (toLists xs)
    mean xs = sum xs / (fromIntegral $ length xs)
    radiusPoints = [KdTree.inRadius tree (d - 1e-15) x | (x, d) <- pds]
    pds = (map Vector.toList $ rows xs) `zip` distances

max_stars_repo_name: eligottlieb/Shannon

max_issues_repo_name: eligottlieb/Shannon

max_forks_repo_name: eligottlieb/Shannon


text:
{-# LANGUAGE GADTs #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ViewPatterns #-}

-- |
-- Module      : PSO.FromPython
-- Description : Communication with Python
-- Copyright   : (c) Tom Westerhout, 2018
-- License     : BSD3
-- Maintainer  : t.westerhout@student.ru.nl
-- Stability   : experimental
module PSO.FromPython
  ( fromPyComplex
  , fromPyList
  , readVector
  , readMatrix
  ) where

import Control.Applicative
import Control.Monad
import Debug.Trace
import Data.Attoparsec.Text
import Data.Complex
import Data.List (transpose)
import Data.Maybe
import Data.Text (Text)
import qualified Data.Text as T
import Data.Text.IO (hGetLine)
import qualified Data.Vector.Storable as V
import GHC.Float
import System.IO (FilePath, Handle, IOMode(..), withFile)

between :: Parser start -> Parser stop -> Parser a -> Parser a
between start stop p = do
  start
  x <- p
  stop
  return x

fromPyComplex :: Parser (Complex Double)
fromPyComplex = pBoth <|> pImag
  where
    pImag = do
      y <- double
      char 'j'
      return (0 :+ y)
    pBoth =
      between (char '(') (char ')') $ do
        x <- double
        y <- double
        char 'j'
        return (x :+ y)

fromPyList :: Parser a -> Parser [a]
fromPyList p = between (char '[') (char ']') (p `sepBy'` string ", ")

readVector :: Text -> Either String (V.Vector (Complex Float))
readVector txt =
  V.fromList <$> fmap (\(x :+ y) -> double2Float x :+ double2Float y) <$>
  parseOnly (fromPyList fromPyComplex <* endOfInput) txt

readMatrix :: Text -> Either String (V.Vector (Complex Float))
readMatrix txt =
  V.fromList <$> fmap (\(x :+ y) -> double2Float x :+ double2Float y) <$> concat <$>
  transpose <$>
    parseOnly (fromPyList (fromPyList fromPyComplex >>= \xs -> trace (show xs) (return xs)) ) txt

max_stars_repo_name: twesterhout/tcm-swarm

max_issues_repo_name: twesterhout/tcm-swarm

max_forks_repo_name: twesterhout/tcm-swarm


text:
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE TupleSections #-}
{-# LANGUAGE ImplicitParams #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Microbenchmark where

import Test.HUnit hiding (State)
import Test.Framework
import Test.Framework.Providers.HUnit

import           System.Random
import           Control.Monad.State
import           Control.Monad
import           Control.Lens

import GHC.Conc.Sync (setNumCapabilities)

import qualified Data.HashMap.Strict       as HM
import qualified Data.HashSet              as Set
import           Data.List
import qualified Data.Vector               as V
import           Data.IORef
import           Data.Time.Clock.POSIX
import           Data.Aeson                as AE
import           Data.ByteString.Lazy      as BS (writeFile)
import           Kvservice_Types
import           Kvstore.KVSTypes

import           Requests
import           TestSetup
import           ServiceConfig
import           Versions

import Options.Applicative
import Data.Semigroup ((<>))

import           Statistics.Sample (mean)

import           Debug.Trace

valueTemplate v = "value-" ++ show v
fieldTemplate f = "field-" ++ show f
keyTemplate k = "key-" ++ show k
tableTemplate t = "table-" ++ show t

data RangeGen = RangeGen Int Int StdGen
instance RandomGen RangeGen where
  next (RangeGen lo hi g) = let (i,g') = randomR (lo,hi) g in (i, RangeGen lo hi g')
  split (RangeGen lo hi g) = let (g1,g2) = split g in (RangeGen lo hi g1, RangeGen lo hi g2)

newtype LinearGen = LinearGen Int
instance RandomGen LinearGen where
  next (LinearGen i) = (i, LinearGen $ i+1)
  split (LinearGen i) = (LinearGen i, LinearGen i)

data BenchmarkState g = BenchmarkState {
                        _fieldCount :: Int,
                        _fieldSelection :: RangeGen,
                        _valueSizeGen :: RangeGen,
                        _tableCount :: Int,
                        _tableSelection :: RangeGen,
                        _keySelection :: g,
                        _operationSelection :: RangeGen,
                        _fieldCountSelection :: RangeGen,
                        _scanCountSelection :: RangeGen
                      }

makeLenses ''BenchmarkState
-- lens does not work with RankNTypes :(

createValue :: forall g. RandomGen g => StateT (BenchmarkState g) IO String
createValue = do
  bmState <- get
  let valSizeGenerator = view valueSizeGen bmState
      (size,valSizeGenerator') = next valSizeGenerator
  put $ over valueSizeGen (const valSizeGenerator') bmState
  (return . concatMap valueTemplate . take size) [1,2..]

createRequest :: forall g. RandomGen g => StateT (BenchmarkState g) IO KVRequest
createRequest = do
  bmState <- get
  let opSelector = view operationSelection bmState
      (op,opSelector') = next opSelector
      bmState' = over operationSelection (const opSelector') bmState

      tableSelector = view tableSelection bmState
      (tableId,tableSelector') = next tableSelector
      table = tableTemplate tableId
      bmState'' = over tableSelection (const tableSelector') bmState'

      keySelector = view keySelection bmState
      (keyId,keySelector') = next keySelector
      key = keyTemplate keyId
      bmState''' = over keySelection (const keySelector') bmState''
  put bmState'''
  case op of
        0 -> prepareINSERT table key <$> createINSERTEntry
        1 -> prepareUPDATE table key <$> createUPDATEEntry
        2 -> prepareREAD table key . Set.map fieldTemplate <$> getFields
        3 -> prepareSCAN table key <$> getScanCount <*> (Set.map fieldTemplate <$> getFields)
        4 -> return $ prepareDELETE table key
        _ -> error $ "No such operation: " ++ show op
    where
      getFieldsAndValues :: forall g. RandomGen g => [Int] -> StateT (BenchmarkState g) IO (HM.HashMap String String)
      getFieldsAndValues = fmap HM.fromList . mapM (\ i -> (,) <$> pure (fieldTemplate i) <*> createValue)
      createINSERTEntry :: forall g. RandomGen g => StateT (BenchmarkState g) IO (HM.HashMap String String)
      createINSERTEntry = getFieldsAndValues . flip take  [1,2..] . view fieldCount =<< get
      getFields :: forall g. RandomGen g => StateT (BenchmarkState g) IO (Set.HashSet Int)
      getFields = do
        s <- get
        let fieldCountGen = view fieldCountSelection s
            (fieldCount,fieldCountGen') = next fieldCountGen
            s' = over fieldCountSelection (const fieldCountGen') s
            fieldSel = view fieldSelection s'
            (fieldSel', fields) = foldl (\ (sel,l) _ -> let (f,sel') = next sel in (sel',l ++ [f]) )
                                       (fieldSel,[])
                                       $ take fieldCount [1,2..]
            s'' = over fieldSelection (const fieldSel') s'
        put s''
        return $ Set.fromList fields
      createUPDATEEntry :: forall g. RandomGen g => StateT (BenchmarkState g) IO (HM.HashMap String String)
      createUPDATEEntry = getFieldsAndValues . Set.toList =<< getFields
      getScanCount = do
        s <- get
        let scanCountSel = view scanCountSelection s
            (scanCount,scanCountSel') = next scanCountSel
        put $ over scanCountSelection (const scanCountSel') s
        return scanCount

workload :: forall g. RandomGen g => Int -> StateT (BenchmarkState g) IO (V.Vector KVRequest)
workload operationCount = V.fromList <$> mapM (const createRequest) [1..operationCount]

showState :: KVSState MockDB -> IO String
showState (KVSState cache dbRef _ _) = do
  db <- readIORef dbRef
  return $ "Cache:\n" ++ show cache ++ "\nDB:\n" ++ show db

loadDB keyCount = do
  s <- initState
  -- fill the db first
  (requests,_) <- runStateT (workload keyCount) $ BenchmarkState
                                            10 -- _fieldCount
                                            (RangeGen 0 10 $ mkStdGen 0) -- _fieldSelection
                                            (RangeGen 5 10 $ mkStdGen 0) -- _valueSizeGen
                                            1 -- _tableCount
                                            (RangeGen 1 1 $ mkStdGen 0) -- _tableSelection
                                            (LinearGen 1) -- _keySelection
                                            (RangeGen 0 0 $ mkStdGen 0) -- _operationSelection (INSERT only)
                                            (RangeGen 3 10 $ mkStdGen 0) -- _fieldCountSelection
                                            (RangeGen 5 10 $ mkStdGen 0) -- _scanCountSelection
  -- traceM "requests (INSERT):"
  -- mapM (\i -> traceM $ show i ++ "\n" ) requests
  (responses, s') <- flip runStateT s $ ?execRequests requests
  responses `seq` assertEqual "wrong number of responses" keyCount $ length responses
  return s'

  -- traceM "state after init:"
  -- traceM =<< showState s'
  -- traceM "done with insert."

reqBenchmarkState keyCount = BenchmarkState
                                10 -- _fieldCount
                                (RangeGen 0 10 $ mkStdGen 0) -- _fieldSelection
                                (RangeGen 5 10 $ mkStdGen 0) -- _valueSizeGen
                                1 -- _tableCount
                                (RangeGen 1 1 $ mkStdGen 0) -- _tableSelection
                                (RangeGen 1 keyCount $ mkStdGen 0) -- _keySelection
                                (RangeGen 1 3 $ mkStdGen 0) -- _operationSelection (no INSERT, no DELETE)
                                (RangeGen 3 10 $ mkStdGen 0) -- _fieldCountSelection
                                (RangeGen 5 10 $ mkStdGen 0) -- _scanCountSelection

currentTimeMillis = round . (* 1000) <$> getPOSIXTime

-- then run some requests
-- runRequests :: (?execRequests :: ExecReqFn)
--             => Int -> Int -> BenchmarkState RangeGen -> KVSState MockDB -> IO (KVSState MockDB, Integer)
runRequests operationCount keyCount bmState s = do
  (requests,_) <- runStateT (workload operationCount) bmState
  -- traceM $ "requests:"
  -- mapM (\i -> traceM $ show i ++ "\n" ) requests
  start <- currentTimeMillis
  (responses, s') <- requests `seq` flip runStateT s $ ?execRequests requests
  stop <- currentTimeMillis
  let execTime = stop - start
  -- traceM "???????????????????????????????????????"
  -- traceM $ "responses:"
  -- mapM (\i -> traceM $ show i ++ "\n" ) responses
  responses `seq` assertEqual "wrong number of responses" operationCount $ length responses
  return (s',execTime)

runSingleBatch :: (?execRequests :: ExecReqFn) => Int -> Int -> IO ()
runSingleBatch keyCount reqCount = do
  (_, execTime) <- runRequests reqCount keyCount (reqBenchmarkState keyCount) =<< loadDB keyCount
  traceM $ "exec time: " ++ show execTime
  return ()

runMultipleBatches :: (?execRequests :: ExecReqFn) => Int -> Int -> Int -> IO Double
runMultipleBatches keyCount reqCount batchCount = do
   s <- loadDB keyCount
   let bmState = reqBenchmarkState keyCount
   (_, execTimes) <- foldM (\ (st,execTimes) _ -> do
                                (st',execTime) <- runRequests reqCount keyCount bmState st
                                return $ st' `seq` (st',execTimes ++ [execTime]))
                           (s,[])
                           $ take batchCount [1,2..]
   let meanExecTime = mean $ V.fromList $ map fromIntegral execTimes
   -- traceM $ "mean execution time: " ++ show meanExecTime ++ " ms"
   return meanExecTime

runMultipleBatches_ :: (?execRequests :: ExecReqFn) => Int -> Int -> Int -> IO ()
runMultipleBatches_ keyCount reqCount batchCount = (\e -> traceM $ "mean execution time: " ++ show e ++ " ms") =<< runMultipleBatches keyCount reqCount batchCount

data ScalabilityResult = ScalabilityResult { version :: String
                                           , results :: [Double] }

scalability name numThreads = do
    results <- foldM (\res n -> do
                      _ <- setNumCapabilities n
                      r <- runMultipleBatches 2000 20 20
                      return $ res ++ [r])
                    []
                    $ take numThreads [1,2..]
    return (name,results)

buildSuite (BenchmarkConfig maxThreadCount) =
  [testCase "Done!" runAsSingleTest]
  where
    runAsSingleTest = do
      results <- forM versions $ \(version, name) -> let ?execRequests = version in scalability name maxThreadCount
      BS.writeFile "scalability.json" $ AE.encode $ HM.fromList results

data BenchmarkConfig = BenchmarkConfig { maxThreadCount :: Int }

benchmarkOptionsParser = BenchmarkConfig <$> (read <$> strOption
                                                         ( short 't'
                                                        <> long "num-threads"
                                                        <> metavar "NUM_THREADS"
                                                        <> help "Maximum number of threads."
                                                        <> showDefault
                                                        <> value "8" ))

max_stars_repo_name: sertel/kvstore

max_issues_repo_name: sertel/kvstore

max_forks_repo_name: sertel/kvstore


text:
-- example 3.19

import Statistics.Distribution
import Statistics.Distribution.Poisson

-- Let a half-hour be our unital time.
pd :: PoissonDistribution
pd = poisson 1.0

p :: Int -> Double
p = probability pd

main = do
  putStrLn $ "the patrol will miss during the 0th period: " ++ show (p 0)
  putStrLn $ "the patrol will miss during the 1th period: " ++ show (p 1)
  putStrLn $ "the patrol will miss during the 2th period: " ++ show (p 2)
  putStrLn $ "the patrol will miss at least once:         " ++ show (1 - p 0)

max_stars_repo_name: raydsameshima/Stat

max_issues_repo_name: raydsameshima/Stat

max_forks_repo_name: raydsameshima/Stat


text:
{-# LANGUAGE DeriveDataTypeable #-}

module BBrotConf where

import Data.Complex
import System.Console.CmdArgs

loCorner, hiCorner :: Complex Double
loCorner = (-2.0) :+ (-1.5)
hiCorner =   1.0  :+   1.5

data ColorPalette = Gray | Reddish | Flames
                    deriving (Show, Data, Typeable)

data Curve = Line | Square | Root
             deriving (Show, Data, Typeable)

data BBrotConf = Compute { seed       :: Maybe Int
                         , samples    :: Int
                         , minIters   :: Int
                         , maxIters   :: Int
                         , ocachepath :: Maybe String
                         , gridStep   :: Double
                         }
               | Render { xpixels     :: Int
                        , ypixels     :: Int
                        , icachepath  :: String
                        , isComplex   :: Bool
                        , dontRender  :: Bool
                        , imagepath   :: Maybe String
                        , palette     :: ColorPalette
                        , curve       :: Curve
                        }
               | ShowCells { gridStep :: Double
                           , maxIters :: Int
                           , animpath :: String
                           }
                 deriving (Show, Data, Typeable)

max_stars_repo_name: saffroy/buddhabrot

max_issues_repo_name: saffroy/buddhabrot

max_forks_repo_name: saffroy/buddhabrot


text:
module Fixtures where

import qualified Statistics.Matrix as M
import qualified Statistics.Matrix.Types as T

f :: Int -> Int -> Double
f a b = (c + d + 1) / (d + 1) where
  c = fromIntegral a
  d = fromIntegral b

g :: Int -> Int -> Double
g a b = (c + d + 1) / (c * d + 1) where
  c = fromIntegral a
  d = fromIntegral b

c :: [[Double]]
c = [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]

d :: [[Double]]
d = [[3.0, 2.5, 1.6], [5.5, 2.3, 6.9], [3.7, 9.1, 10.0]]

matA :: T.Matrix
matA = M.generate 3 3 f

matB :: T.Matrix
matB = M.generate 3 3 g

matC :: T.Matrix
matC = M.fromRowLists c

matD :: T.Matrix
matD = M.fromRowLists d

matId :: T.Matrix
matId = M.ident 3

max_stars_repo_name: Magalame/dh-core

max_issues_repo_name: Magalame/dh-core

max_forks_repo_name: Magalame/dh-core


text:
{-# LANGUAGE BangPatterns     #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE Strict #-}
module STC.Convolution where

import           Control.Arrow
import           Control.Monad              as M
import           Control.Monad.Parallel     as MP (bindM2, mapM)
import           Data.Array.IArray          as IA
import           Data.Array.Repa            as R
import           Data.Complex
import           Data.Ix
import           Data.List                  as L
import           Data.Vector.Storable       as VS
import           Data.Vector.Unboxed        as VU
import           Debug.Trace
import           DFT.Plan
import           Filter.Utils
import           FokkerPlanck.FourierSeries
import           Image.IO
import           STC.DFTArray
import           Types
import           Utils.Parallel
import FokkerPlanck.GPUKernel
import qualified Data.Array.Accelerate              as A
import qualified Data.Array.Accelerate.Data.Complex as A
import           Data.Array.Accelerate.LLVM.PTX
import Debug.Trace
import Utils.List
import Utils.BLAS

data Field
  = Source
  | Sink
  deriving (Read, Show)

{-# INLINE dftHarmonicsArray #-}
dftHarmonicsArray ::
     DFTPlan
  -> Int
  -> Double
  -> Int
  -> Double
  -> [Double]
  -> [Double]
  -> [Double]
  -> [Double]
  -> Double
  -> Double
  -> IO (IA.Array (Int, Int) (VS.Vector (Complex Double)))
dftHarmonicsArray !plan !numRows !deltaRow !numCols !deltaCol !phiFreqs !rhoFreqs !thetaFreqs !rFreqs !halfLogPeriod !cutoff = do
  let !arr =
        computeHarmonicsArray
          numRows
          deltaRow
          numCols
          deltaCol
          phiFreqs
          rhoFreqs
          thetaFreqs
          rFreqs
          halfLogPeriod
          cutoff
  fmap (listArray (bounds arr)) .
    dftExecuteBatchP plan (DFTPlanID DFT1DG [numCols, numRows] [0, 1]) .
    L.map
      (VS.convert .
       toUnboxed .
       computeUnboxedS . makeFilter2D . fromUnboxed (Z :. numCols :. numRows)) .
    IA.elems $
    arr

{-# INLINE dftHarmonicsArrayG #-}
dftHarmonicsArrayG ::
     DFTPlan
  -> Int
  -> Double
  -> Int
  -> Double
  -> [Double]
  -> [Double]
  -> [Double]
  -> [Double]
  -> Double
  -> Double -> VS.Vector (Complex Double)
  -> IO (IA.Array (Int, Int) (VS.Vector (Complex Double)))
dftHarmonicsArrayG !plan !numRows !deltaRow !numCols !deltaCol !phiFreqs !rhoFreqs !thetaFreqs !rFreqs !halfLogPeriod !cutoff !gaussian = do
  let !arr =
        computeHarmonicsArray
          numRows
          deltaRow
          numCols
          deltaCol
          phiFreqs
          rhoFreqs
          thetaFreqs
          rFreqs
          halfLogPeriod
          cutoff
  xs <-
    fmap (L.map (VS.zipWith (*) gaussian)) .
    dftExecuteBatchP plan (DFTPlanID DFT1DG [numCols, numRows] [0, 1]) .
    IA.elems $
    arr
  ys <- dftExecuteBatchP plan (DFTPlanID IDFT1DG [numCols, numRows] [0, 1]) xs
  fmap (listArray (bounds arr)) .
    dftExecuteBatchP plan (DFTPlanID DFT1DG [numCols, numRows] [0, 1]) .
    L.map
      (\vec ->
         let arr = fromUnboxed (Z :. numCols :. numRows) . VS.convert $ vec
         in VS.convert .
            toUnboxed . computeUnboxedS . makeFilter2D . R.traverse arr id $ \f idx@(Z :. i :. j) ->
              if (i == div numCols 2 && j == div numRows 2)
                 then 0
                 else f idx) $
    ys

{-# INLINE convolve #-}
convolve ::
     Field
  -> DFTPlan
  -> R.Array U DIM4 (Complex Double)
  -> IA.Array (Int, Int) (VS.Vector (Complex Double))
  -> DFTArray
  -> IO DFTArray
convolve !field !plan !coefficients !harmonicsArray !arr@(DFTArray rows cols thetaFreqs rFreqs vecs) = do
  dftVecs <- dftExecuteBatchP plan (DFTPlanID DFT1DG [cols, rows] [0, 1]) vecs
  let !initVec = VS.replicate (VS.length . L.head $ vecs) 0
      idx = (,) <$> (L.zip [0 ..] rFreqs) <*> (L.zip [0 ..] thetaFreqs)
  fmap (DFTArray rows cols thetaFreqs rFreqs) .
    dftExecuteBatchP plan (DFTPlanID IDFT1DG [cols, rows] [0, 1]) .
    parMap
      rdeepseq
      (\((!r, !rFreq), (!theta, !thetaFreq)) ->
         L.foldl'
           (\(!vec) (((!rho, !rhoFreq), (!phi, !phiFreq)), inputVec) ->
              VS.zipWith
                (+)
                vec
                (VS.map (* (coefficients R.! (Z :. r :. theta :. rho :. phi))) .
                 VS.zipWith
                   (*)
                   (getHarmonics harmonicsArray phiFreq rhoFreq thetaFreq rFreq) $
                 inputVec))
           initVec .
         L.zip idx $
         dftVecs) $
    idx

{-# INLINE convolve' #-}
convolve' ::
     Field
  -> DFTPlan
  -> R.Array U DIM4 (Complex Double)
  -> IA.Array (Int, Int) (VS.Vector (Complex Double))
  -> DFTArray
  -> IO DFTArray
convolve' !field !plan !coefficients !harmonicsArray !arr@(DFTArray rows cols thetaFreqs rFreqs vecs) = do
  dftVecs <- dftExecuteBatchP plan (DFTPlanID DFT1DG [cols, rows] [0, 1]) vecs
  let !initVec = VS.replicate (VS.length . L.head $ vecs) 0
      idxTheta = L.zip [0 ..] thetaFreqs
  fmap (DFTArray rows cols thetaFreqs rFreqs) .
    dftExecuteBatchP plan (DFTPlanID IDFT1DG [cols, rows] [0, 1]) .
    parMap
      rdeepseq
      (\(!theta, !thetaFreq) ->
         VS.map
           (\x ->
              case field of
                Source -> x
                Sink -> x * cis (pi * thetaFreq)) .
         L.foldl'
           (\vec1 (rho, rhoFreq) ->
              L.foldl'
                (\(!vec2) ((!phi, !phiFreq), inputVec) ->
                   VS.zipWith
                     (+)
                     vec2
                     (VS.map
                        (* (coefficients R.!
                            (Z :. (0 :: Int) :. theta :. rho :. phi))) .
                      VS.zipWith
                        (*)
                        (getHarmonics harmonicsArray phiFreq rhoFreq thetaFreq 0) $
                      inputVec))
                vec1 .
              L.zip idxTheta $
              dftVecs)
           initVec .
         L.zip [0 ..] $
         rFreqs) $
    idxTheta

{-# INLINE convolveG' #-}
convolveG' ::
     Field
  -> DFTPlan
  -> R.Array U DIM4 (Complex Double)
  -> IA.Array (Int, Int) (VS.Vector (Complex Double))
  -> VS.Vector (Complex Double)
  -> DFTArray
  -> IO DFTArray
convolveG' !field !plan !coefficients !harmonicsArray !gaussianFilter !arr@(DFTArray rows cols thetaFreqs rFreqs vecs) = do
  dftVecs <-
    L.map (VS.zipWith (*) gaussianFilter) <$>
    dftExecuteBatchP plan (DFTPlanID DFT1DG [cols, rows] [0, 1]) vecs
  let !initVec = VS.replicate (VS.length . L.head $ vecs) 0
      idxTheta = L.zip [0 ..] thetaFreqs
  fmap (DFTArray rows cols thetaFreqs rFreqs) .
    dftExecuteBatchP plan (DFTPlanID IDFT1DG [cols, rows] [0, 1]) .
    parMap
      rdeepseq
      (\(!theta, !thetaFreq) ->
         L.foldl'
           (\vec1 (rho, rhoFreq) ->
              L.foldl'
                (\(!vec2) ((!phi, !phiFreq), inputVec) ->
                   VS.zipWith
                     (+)
                     vec2
                     (VS.map
                        (* (case field of
                              Source ->
                                coefficients R.!
                                (Z :. (0 :: Int) :. theta :. rho :. phi)
                              Sink ->
                                (coefficients R.!
                                 (Z :. (0 :: Int) :. theta :. rho :. phi)) *
                                (cis $ (-(phiFreq + thetaFreq)) * pi))) .
                      VS.zipWith
                        (*)
                        (getHarmonics harmonicsArray phiFreq rhoFreq thetaFreq 0) $
                      inputVec))
                vec1 .
              L.zip idxTheta $
              dftVecs)
           initVec .
         L.zip [0 ..] $
         rFreqs) $
    idxTheta

-- {-# INLINE convolve'' #-}
-- convolve'' ::
--      Field
--   -> DFTPlan
--   -> R.Array U DIM4 (Complex Double)
--   -> IA.Array (Int, Int) (VS.Vector (Complex Double))
--   -> DFTArray
--   -> IO DFTArray
-- convolve'' !field !plan !coefficients !harmonicsArray !arr@(DFTArray rows cols phiFreqs rhoFreqs vecs) = do
--   let !initVec = VS.replicate (VS.length . L.head $ vecs) 0
--       idx = L.zip [0 ..] phiFreqs
--   return .
--     (DFTArray rows cols phiFreqs rhoFreqs) .
--     parMap
--       rdeepseq
--       (\(theta, thetaFreq) ->
--          L.foldl'
--            (\vec1 (rho, rhoFreq) ->
--               L.foldl'
--                 (\vec2 (phi, phiFreq) ->
--                    VS.zipWith (+) vec2 .
--                    VS.map
--                      (* (case field of
--                            Source ->
--                              coefficients R.!
--                              (Z :. (0 :: Int) :. theta :. rho :. phi)
--                            Sink ->
--                              (coefficients R.!
--                               (Z :. (0 :: Int) :. theta :. rho :. phi)) *
--                              (cis $ (-(thetaFreq + phiFreq)) * pi))) $
--                    (getHarmonics harmonicsArray phiFreq rhoFreq 0 0))
--                 vec1
--                 idx)
--            initVec .
--          L.zip [0 ..] $
--          rhoFreqs) $
--     idx

-- {-# INLINE dftHarmonicsArrayGPU #-}
-- dftHarmonicsArrayGPU ::
--      DFTPlan
--   -> Int
--   -> Double
--   -> Int
--   -> Double
--   -> [Double]
--   -> [Double]
--   -> [Double]
--   -> [Double]
--   -> Double
--   -> Double
--   -> IO (Acc (A.Array A.DIM4 (A.Complex Double)))
-- dftHarmonicsArrayGPU !plan !numRows !deltaRow !numCols !deltaCol !phiFreqs !rhoFreqs !thetaFreqs !rFreqs !halfLogPeriod !cutoff = do
--   let !arr =
--         computeHarmonicsArrayGPU
--           numRows
--           deltaRow
--           numCols
--           deltaCol
--           phiFreqs
--           rhoFreqs
--           thetaFreqs
--           rFreqs
--           halfLogPeriod
--           cutoff
--       rangeFunc2 xs ys =
--         (round (L.head xs - L.last ys), round (L.last xs - L.head ys))
--       (!tfLB, !tfUB) = rangeFunc2 phiFreqs thetaFreqs
--       (!rfLB, !rfUB) = rangeFunc2 rhoFreqs rFreqs
--   fmap
--     (A.use .
--      A.fromList
--        (A.Z A.:. (rfUB - rfLB + 1) A.:. (tfUB - tfLB + 1) A.:. numCols A.:.
--         numRows) .
--      VS.toList . VS.concat) .
--     dftExecuteBatchP plan (DFTPlanID DFT1DG [numCols, numRows] [0, 1]) .
--     L.map
--       (VS.convert .
--        toUnboxed .
--        computeUnboxedS . makeFilter2D . fromUnboxed (Z :. numCols :. numRows)) .
--     IA.elems $
    -- arr

-- {-# INLINE convolveGPU #-}
-- convolveGPU ::
--      PTX
--   -> Field
--   -> DFTPlan
--   -> Acc (A.Array A.DIM4 (A.Complex Double))
--   -> Acc (A.Array A.DIM4 (A.Complex Double))
--   -> Acc (A.Array A.DIM4 (A.Complex Double))
--   -> DFTArray
--   -> IO DFTArray
-- convolveGPU !ptx !field !plan !coefficients !coefficientsSink !harmonics !arr@(DFTArray rows cols thetaFreqs rFreqs vecs) = do
--   dftVecs <- dftExecuteBatchP plan (DFTPlanID DFT1DG [cols, rows] [0, 1]) vecs
--   let !numRFreq = L.length rFreqs
--       !numThetaFreq = L.length thetaFreqs
--       idx =
--         [ (r, theta)
--         | r <- [0 .. (numRFreq - 1)]
--         , theta <- [0 .. (numThetaFreq - 1)]
--         ]
--   case field of
--     Source ->
--       fmap (DFTArray rows cols thetaFreqs rFreqs) .
--       dftExecuteBatchP plan (DFTPlanID IDFT1DG [cols, rows] [0, 1]) .
--       L.map
--         (\(r, theta) ->
--            VS.fromList . A.toList $
--            (runNWith ptx (convolveKernel coefficients harmonics))
--              (A.fromList A.Z [theta])
--              (A.fromList A.Z [r])
--              (A.fromList
--                 (A.Z A.:. numRFreq A.:. numThetaFreq A.:. cols A.:. rows) .
--               VS.toList . VS.concat $
--               dftVecs)) $
--       idx
--     -- Sink ->
--     --   fmap (DFTArray rows cols thetaFreqs rFreqs) .
--     --   dftExecuteBatchP plan (DFTPlanID IDFT1DG [cols, rows] [0, 1]) .
--     --   L.map
--     --     (\(r, theta) ->
--     --        VS.fromList .
--     --        A.toList .
--     --        (runNWith ptx (convolveKernel coefficientsSink harmonics))
--     --          (A.fromList (Z :. (1 :: Int)) theta)
--     --          (A.fromList (Z :. (1 :: Int)) r)
--     --          (A.fromList (Z :. numRFreq :. numThetaFreq :. cols :. rows) .
--     --           VS.toList . VS.concat $
--     --           dftVecs)) $
--     --   idx

{-# INLINE convolveSingle #-}
convolveSingle ::
     Field
  -> R.Array U DIM4 (Complex Double)
  -> IA.Array (Int, Int) (R.Array U DIM2 (Complex Double))
  -> R.Array U DIM1 Double
  -> R.Array U DIM1 Double
  -> Double
  -> Double
  -> R.Array U DIM2 (Complex Double)
  -> VU.Vector (Complex Double)
convolveSingle !field !coefficients !harmonicsArray !thetaFreqs !rFreqs !x !y !input =
  let (Z :. cols :. rows) = extent (harmonicsArray IA.! (0, 0))
      rowCenter = div rows 2
      colCenter = div cols 2
      (Z :. _ :. _ :. numRFreq :. numThetaFreq) = extent coefficients
      -- harmonics =
      --   R.traverse2 thetaFreqs rFreqs (\_ _ -> extent coefficients) $ \fThetaFreq fRFreq (Z :. r :. theta :. rho :. phi) ->
      --     (harmonicsArray IA.!
      --      ( round (fRFreq (Z :. rho) - fRFreq (Z :. r))
      --      , round $ (fThetaFreq (Z :. phi) - (fThetaFreq (Z :. theta))))) R.!
      --     (Z :. (x + colCenter) :. (y + rowCenter))
      harmonics =
        R.traverse2 thetaFreqs rFreqs (\_ _ -> extent coefficients) $ \fPhiFreq fRhoFreq (Z :. r :. theta :. rho :. phi) ->
          let !tf = fPhiFreq (Z :. phi) - fPhiFreq (Z :. theta)
              !rf = fRhoFreq (Z :. rho) - fRhoFreq (Z :. r)
          in (x :+ y) ** (tf :+ 0) *
             ((x ^ 2 + y ^ 2) :+ 0) **
             (((-tf - 1) :+ 2 * pi * rf / (log $ (sqrt 2) * 512)) / 2)
      product =
        R.traverse2 (R.zipWith (*) coefficients harmonics) input const $ \fCoefficient fInput idx@(Z :. _ :. _ :. rho :. phi) ->
          fCoefficient idx * fInput (Z :. rho :. phi)
      arr =
        case field of
          Source -> product
          Sink ->
            error
              "convolveSingle: You should not use it to compute a sink field."
  in if x == 0 && y == 0
       then VU.replicate (numRFreq * numThetaFreq) 0
       else toUnboxed . sumS . sumS $ arr

-- {-# INLINE convolveSparse #-}
convolveSparse ::
     Field
  -> R.Array U DIM4 (Complex Double)
  -> IA.Array (Int, Int) (R.Array U DIM2 (Complex Double))
  -> R.Array U DIM1 Double
  -> R.Array U DIM1 Double
  -> Double
  -> [(Double, Double)]
  -> [R.Array U DIM2 (Complex Double)]
  -> [R.Array U DIM2 (Complex Double)]
convolveSparse !field !coefficients !harmonicsArray !thetaFreqs !rFreqs !cutoff !xs !inputs =
  L.map (fromUnboxed (extent . L.head $ inputs)) .
  parMap
    rdeepseq
    (\(x, y) ->
       L.foldl1' (VU.zipWith (+)) .
       L.zipWith
         (\(x', y') input ->
            let !x'' = x - x'
                !y'' = y - y'
                (!a, !b) =
                  if (x'' ^ 2 + y'' ^ 2) <= cutoff ^ 2
                    then (x'', y'')
                    else (0, 0)
            in convolveSingle
                 field
                 coefficients
                 harmonicsArray
                 thetaFreqs
                 rFreqs
                 a
                 b
                 input)
         xs $
       inputs) $
  xs

{-# INLINE convolveSingle' #-}
convolveSingle' ::
     Field
  -> R.Array U DIM4 (Complex Double)
  -> IA.Array (Int, Int) (R.Array U DIM2 (Complex Double))
  -> R.Array U DIM1 Double
  -> R.Array U DIM1 Double
  -> Double
  -> Double
  -> R.Array U DIM1 (Complex Double)
  -> VU.Vector (Complex Double)
convolveSingle' !field !coefficients !harmonicsArray !phiFreqs !rhoFreqs !x !y !input =
  let (Z :. cols :. rows) = extent (harmonicsArray IA.! (0, 0))
      rowCenter = div rows 2
      colCenter = div cols 2
      (Z :. _ :. _ :. numRhoFreq :. numPhiFreq) = extent coefficients
      -- harmonics =
      --   R.traverse2 phiFreqs rhoFreqs (\_ _ -> extent coefficients) $ \fPhiFreq fRhoFreq (Z :. _ :. theta :. rho :. phi) ->
      --     (harmonicsArray IA.!
      --      ( round (fRhoFreq (Z :. rho))
      --      , round (fPhiFreq (Z :. phi) - (fPhiFreq (Z :. theta))))) R.!
      --     (Z :. (x + colCenter) :. (y + rowCenter))
      harmonics =
        R.traverse2 phiFreqs rhoFreqs (\_ _ -> extent coefficients) $ \fPhiFreq fRhoFreq (Z :. _ :. theta :. rho :. phi) ->
          let !tf = fPhiFreq (Z :. phi) - fPhiFreq (Z :. theta)
              !rf = fRhoFreq (Z :. rho)
          in (x :+ y) ** (tf :+ 0) *
             ((x ^ 2 + y ^ 2) :+ 0) ** (((-tf - 1) :+ rf) / 2)
      product =
        R.traverse2 (R.zipWith (*) coefficients harmonics) input const $ \fCoefficient fInput idx@(Z :. _ :. _ :. _ :. phi) ->
          fCoefficient idx * fInput (Z :. phi)
      arr =
        case field of
          Source -> product
          Sink ->
            R.traverse2 product phiFreqs const $ \fProd fPhiFreq idx@(Z :. _ :. theta :. _ :. phi) ->
              fProd idx *
              (cis $ (-(fPhiFreq (Z :. phi) - fPhiFreq (Z :. theta))) * pi)
  in if x == 0 && y == 0
       then VU.replicate (numRhoFreq * numPhiFreq) 0
       else toUnboxed . sumS . sumS $ arr

convolveSparse' ::
     Field
  -> R.Array U DIM4 (Complex Double)
  -> IA.Array (Int, Int) (R.Array U DIM2 (Complex Double))
  -> R.Array U DIM1 Double
  -> R.Array U DIM1 Double
  -> Double
  -> [(Double, Double)]
  -> [R.Array U DIM1 (Complex Double)]
  -> [R.Array U DIM1 (Complex Double)]
convolveSparse' !field !coefficients !harmonicsArray !phiFreqs !rhoFreqs !cutoff !xs !inputs =
  L.map (fromUnboxed (extent . L.head $ inputs)) .
  parMap
    rdeepseq
    (\(x, y) ->
       L.foldl1' (VU.zipWith (+)) .
       L.zipWith
         (\(x', y') input ->
            let !x'' = x - x'
                !y'' = y - y'
                (!a, !b) =
                  if (x'' ^ 2 + y'' ^ 2) <= cutoff ^ 2
                    then (x'', y'')
                    else (0, 0)
            in convolveSingle'
                 field
                 coefficients
                 harmonicsArray
                 phiFreqs
                 rhoFreqs
                 a
                 b
                 input)
         xs $
       inputs) $
  xs

{-# INLINE gaussianFilter2D #-}
gaussianFilter2D :: DFTPlan -> Int -> Double -> IO (VS.Vector (Complex Double))
gaussianFilter2D plan numPoint stdG =
  dftExecute plan (DFTPlanID DFT1DG [numPoint, numPoint] [0, 1]) .
  VU.convert .
  toUnboxed . computeS . makeFilter2D . fromFunction (Z :. numPoint :. numPoint) $ \(Z :. i :. j) ->
    let x = fromIntegral $ i - div numPoint 2
        y = fromIntegral $ j - div numPoint 2
    in (exp (-(x ^ 2 + y ^ 2) / (2 * stdG ^ 2))) :+ 0

{-# INLINE convoluvePinhweelBasis #-}
convoluvePinhweelBasis ::
     R.Array U DIM4 (Complex Double)
  -> IA.Array (Int, Int) (VS.Vector (Complex Double))
  -> DFTArray
  -> DFTArray
convoluvePinhweelBasis coefficients harmonicsArray arr@(DFTArray xFreq yFreq thetaFreqs rFreqs vecs) =
 let initVec = VS.replicate (VS.length . L.head $ vecs) 0
     idxs = (,) <$> (L.zip [0 ..] rFreqs) <*> (L.zip [0 ..] thetaFreqs)
  in DFTArray xFreq yFreq thetaFreqs rFreqs .
     parMap
      rdeepseq
       (\((r, rFreq), (theta, thetaFreq)) ->
          L.foldl'
            (\vec (((rho, rhoFreq), (phi, phiFreq)), inputVec) ->
               VS.zipWith (+) vec .
               VS.map (* (coefficients R.! (Z :. r :. theta :. rho :. phi))) .
               VS.zipWith
                 (*)
                 (getHarmonics harmonicsArray phiFreq rhoFreq thetaFreq rFreq) $
               inputVec)
            initVec .
          L.zip idxs $
          vecs) $
     idxs

{-# INLINE convoluvePinhweelBasis' #-}
convoluvePinhweelBasis' ::
     R.Array U DIM4 (Complex Double)
  -> IA.Array (Int, Int) (VS.Vector (Complex Double))
  -> DFTArray
  -> DFTArray
convoluvePinhweelBasis' coefficients harmonicsArray arr@(DFTArray xFreq yFreq thetaFreqs rFreqs vecs) =
  let idxTheta = L.zip [0 ..] thetaFreqs
      initVec = VS.replicate (VS.length . L.head $ vecs) 0
  in DFTArray xFreq yFreq thetaFreqs rFreqs .
     parMap
       rdeepseq
       (\(!theta, !thetaFreq) ->
          L.foldl'
            (\vec1 (rho, rhoFreq) ->
               L.foldl'
                 (\(!vec2) ((!phi, !phiFreq), inputVec) ->
                    VS.zipWith
                      (+)
                      vec2
                      (VS.map
                         (* (coefficients R.!
                             (Z :. (0 :: Int) :. theta :. rho :. phi))) .
                       VS.zipWith
                         (*)
                         (getHarmonics
                            harmonicsArray
                            phiFreq
                            rhoFreq
                            thetaFreq
                            0) $
                       inputVec))
                 vec1 .
               L.zip idxTheta $
               vecs)
            initVec .
          L.zip [0 ..] $
          rFreqs) $
     idxTheta

{-# INLINE convolveFull'' #-}
convolveFull'' ::
     R.Array U DIM4 (Complex Double)
  -> IA.Array (Int, Int) (VS.Vector (Complex Double))
  -> DFTArray
  -> DFTArray
convolveFull'' !coefficients !harmonicsArray !arr@(DFTArray r2Freq _ thetaFreqs rFreqs vecs) =
  let idxTheta = L.zip [0 ..] thetaFreqs
      !initVec = VS.replicate (VS.length . L.head $ vecs) 0
  in DFTArray r2Freq r2Freq thetaFreqs rFreqs .
     parMap
       rdeepseq
       (\(!theta, !thetaFreq) ->
          L.foldl'
            (\vec1 (rho, rhoFreq) ->
               L.foldl'
                 (\(!vec2) ((!phi, !phiFreq), inputVec) ->
                    VS.zipWith
                      (+)
                      vec2
                      (VS.map
                         (* (coefficients R.!
                             (Z :. (0 :: Int) :. theta :. rho :. phi)))
                         (getHarmonics
                            harmonicsArray
                            phiFreq
                            rhoFreq
                            thetaFreq
                            0)))
                 vec1 .
               L.zip idxTheta $
               vecs)
            initVec .
          L.zip [0 ..] $
          rFreqs) $
     idxTheta

{-# INLINE convoluvePinhweelBasisTest #-}
convoluvePinhweelBasisTest ::
     VS.Vector (Complex Double)
  -> [VS.Vector (Complex Double)]
  -> [VS.Vector (Complex Double)]
  -> DFTArray
  -> IO DFTArray
convoluvePinhweelBasisTest coefficients harmonicsArray1 harmonicsArray2 arr@(DFTArray xFreq yFreq thetaFreqs rFreqs vecs) = do
  let vecs1 = parZipWith rdeepseq (VS.zipWith (*)) vecs harmonicsArray1
      m = L.length thetaFreqs * L.length rFreqs
      n = xFreq * yFreq
  vec2 <- gemmBLAS m n m coefficients . VS.concat $ vecs1
  let arr =
        fromUnboxed
          (Z :. (L.length rFreqs) :. (L.length thetaFreqs) :. xFreq :. yFreq) .
        VS.convert $
        vec2
      vecs2 =
        parMap
          rdeepseq
          (\(i, j) ->
             VS.convert . toUnboxed . computeS . R.slice arr $
             (Z :. i :. j :. All :. All))
          [ (i, j)
          | i <- [0 .. L.length rFreqs - 1]
          , j <- [0 .. L.length thetaFreqs - 1]
          ] :: [VS.Vector (Complex Double)]
      vecs3 = parZipWith rdeepseq (VS.zipWith (*)) vecs2 harmonicsArray2
  return . DFTArray xFreq yFreq thetaFreqs rFreqs $ vecs3

{-# INLINE convoluvePinhweelBasisTest' #-}
convoluvePinhweelBasisTest' ::
     R.Array U DIM4 (Complex Double)
  -> [VS.Vector (Complex Double)]
  -> [VS.Vector (Complex Double)]
  -> DFTArray
  -> IO DFTArray
convoluvePinhweelBasisTest' coefficients harmonicsArray1 harmonicsArray2 arr@(DFTArray xFreq yFreq thetaFreqs rFreqs vecs) = do
  let (Z :. numRFreq :. numThetaFreq :. numRhoFreq :. numPhiFreq) =
        extent coefficients
      coefVec = VU.convert . toUnboxed $ coefficients
      vecs1 =
        parZipWith rdeepseq (VS.zipWith (*)) harmonicsArray1 .
        L.concat . L.replicate numRhoFreq $
        vecs
      m = numRFreq * numThetaFreq
      n = xFreq * yFreq
      k = numRhoFreq * numPhiFreq
  print (m, n, k)
  vec2 <- gemmBLAS m n k coefVec . VS.concat $ vecs1
  let arr =
        fromUnboxed (Z :. numRFreq :. numThetaFreq :. xFreq :. yFreq) .
        VS.convert $
        vec2
      vecs2 =
        parMap
          rdeepseq
          (\(i, j) ->
             VS.convert . toUnboxed . computeS . R.slice arr $
             (Z :. i :. j :. All :. All))
          [(i, j) | i <- [0 .. numRFreq - 1], j <- [0 .. numThetaFreq - 1]] :: [VS.Vector (Complex Double)]
      vecs3 = parZipWith rdeepseq (VS.zipWith (*)) vecs2 harmonicsArray2
  return . DFTArray xFreq yFreq thetaFreqs rFreqs $ vecs3

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:
module ML.Unsupervised.GMM (
  fit,
  GMM(..),
  GMMComponent(..),
  FitOptions(..),
  defaultOptions,
  predict,
  empty
) where

import Numeric.LinearAlgebra
import qualified Data.Vector.Storable as V
import Statistics.Distribution
import qualified Data.List as L

import qualified ML.Unsupervised.Clustering.KMeans as KMeans

import Debug.Trace

data GMMComponent a = GMMComponent {
  gmmcPi :: Double,
  gmmcMean :: V.Vector a,
  gmmcSigma :: Matrix a
} deriving (Show)

data GMM a = GMM {
  mixtureComponents :: [GMMComponent a]
} deriving (Show)

empty :: GMM a
empty = GMM { mixtureComponents = [] }

normalPdf :: V.Vector Double -> Matrix Double -> V.Vector Double -> Double
normalPdf mu sigma x = (1 / sqrt (2 * pi * det sigma)) * exp (-0.5 * (x - mu) <.> ((inv sigma) #> (x - mu)))

vecsum (v:vs) = foldr (V.zipWith (+)) v vs
vecsum [] = error "Empty list in vecsum"

data FitOptions = FitOptions {
  foMaxIterations :: Int,
  foLikelihoodConvergenceRatio :: Double
}

defaultOptions = FitOptions {
  foMaxIterations = 100,
  foLikelihoodConvergenceRatio = 1.0e-12
}

likelihood :: GMM Double -> [V.Vector Double] -> Double
likelihood gmm points = sum $ fmap (logL gmm) points
  where
    logL gmm point = log $ sum $ fmap (logLFromComponent point) (mixtureComponents gmm)
    logLFromComponent point comp = (gmmcPi comp) * normalPdf (gmmcMean comp) (gmmcSigma comp) point

fit :: FitOptions -> Int -> [V.Vector Double] -> GMM Double
fit options nComponents points = fit' (GMM initialComponents) 0
  where
    fit' gmm itnum = let updatedGmm = update gmm in
                if (abs (likelihood updatedGmm points / likelihood gmm points - 1) < foLikelihoodConvergenceRatio options) || (itnum >= foMaxIterations options)
                  then gmm
                  else fit' updatedGmm (trace ("fit' " ++ show itnum) itnum + 1)

    calculateResponsibility gmm comp = zipWith (/) (fmap (respNumerator comp) points) (fmap (respDenominator gmm) points)

    respNumerator component point = (gmmcPi component) * normalPdf (gmmcMean component) (gmmcSigma component) point

    respDenominator gmm point = sum $ fmap (\x -> respNumerator x point) (mixtureComponents gmm)

    update gmm = GMM $ fmap (updateComponent gmm) (mixtureComponents gmm)

    updateComponent gmm comp = comp {gmmcPi = updatePi gmm comp, gmmcMean = updateMu gmm comp, gmmcSigma = updateSigma gmm comp }

    updatePi gmm comp = (sum $ calculateResponsibility gmm comp) / (fromIntegral $ length points)

    updateMu gmm comp = V.map (/ (sum $ calculateResponsibility gmm comp)) (vecsum $ zipWith (\x r -> V.map (* r) x) points (calculateResponsibility gmm comp))

    updateSigma gmm comp = (sum $ zipWith (zipR comp) (calculateResponsibility gmm comp) points) / scalar (sum $ calculateResponsibility gmm comp)

    zipR :: GMMComponent Double -> Double -> V.Vector Double -> Matrix Double
    zipR comp x p = scalar x * ((p - gmmcMean comp) `outer` (p - gmmcMean comp))

    dim = V.length $ head points

    initialComponents = let km = KMeans.fit [1..(nComponents)] points
                            predicted = fmap (\p -> (KMeans.predict km p, p)) points in
                          fmap (mkComponent predicted) (KMeans.centroids km)

    mkComponent predicted cent = let filteredPoints = filter (\x -> fst x == fst cent) predicted in
                         GMMComponent { gmmcPi = fromIntegral (length filteredPoints) / fromIntegral (length points),
                                        gmmcMean = snd cent,
                                        gmmcSigma = ident dim * scalar (sqrt ((foldr (\v x -> x + dist (snd cent) (snd v)) 0 filteredPoints) / fromIntegral (length points))) }

    dist :: V.Vector Double -> V.Vector Double -> Double
    dist v1 v2 = (v1 - v2) <.> (v1 - v2)

predict :: GMM Double -> V.Vector Double -> Int
predict gmm point = fst $ L.maximumBy (\x y -> snd x `compare` snd y) $ zip [0..] $ fmap (\comp -> (gmmcPi comp) * normalPdf (gmmcMean comp) (gmmcSigma comp) point) (mixtureComponents gmm)

max_stars_repo_name: asakul/small-ml

max_issues_repo_name: asakul/small-ml

max_forks_repo_name: asakul/small-ml


text:
module GDAL.Complex where

import qualified Data.Complex as C

newtype Complex a =
  Complex (C.Complex a)
  deriving (Eq)

instance Show a => Show (Complex a) where
  show (Complex c) = show c

instance Read a => Read (Complex a) where
  read = Complex . read

infix 6 :+
(:+) :: a -> a -> Complex a
a :+ b = Complex (a C.:+ b)

{-
c
Data.Complex.:+         Data.Complex.cis        Data.Complex.imagPart   Data.Complex.mkPolar    Data.Complex.polar
Data.Complex.Complex    Data.Complex.conjugate  Data.Complex.magnitude  Data.Complex.phase      Data.Complex.realPart
-}

max_stars_repo_name: peddie/bindings-gdal

max_issues_repo_name: peddie/bindings-gdal

max_forks_repo_name: peddie/bindings-gdal


text:
import AI.HNN.FF.Network
import Numeric.LinearAlgebra

samples :: Samples Double
samples = [
  (fromList [ 1, 1, 1, 1, 1
            , 0, 0, 0, 0, 1
            , 0, 0, 1, 1, 1
            , 0, 0, 0, 0, 1
            , 1, 1, 1, 1, 1 ], fromList [1]),   -- three

  (fromList [ 1, 1, 1, 1, 1
            , 0, 0, 0, 0, 1
            , 1, 1, 1, 1, 1
            , 0, 0, 0, 0, 1
            , 1, 1, 1, 1, 1 ], fromList [1]),   -- three

  (fromList [ 0, 1, 1, 1, 1
            , 0, 0, 0, 0, 1
            , 0, 0, 0, 1, 1
            , 0, 0, 0, 0, 1
            , 0, 1, 1, 1, 1 ], fromList [1]),   -- three

  (fromList [ 1, 1, 1, 1, 1
            , 1, 0, 0, 0, 1
            , 1, 0, 0, 0, 1
            , 1, 0, 0, 0, 1
            , 1, 1, 1, 1, 1 ], fromList [0]),   -- not a three

  (fromList [ 1, 1, 1, 1, 1
            , 1, 0, 0, 0, 0
            , 1, 0, 0, 0, 0
            , 1, 0, 0, 0, 0
            , 1, 0, 0, 0, 0 ], fromList [0]),   -- not a three

  (fromList [ 0, 1, 1, 1, 0
            , 0, 1, 0, 1, 0
            , 0, 1, 1, 1, 0
            , 0, 1, 0, 1, 0
            , 0, 1, 1, 1, 0 ], fromList [0]),   -- not a three

  (fromList [ 0, 0, 1, 0, 0
            , 0, 1, 1, 0, 0
            , 1, 0, 1, 0, 0
            , 0, 0, 1, 0, 0
            , 0, 0, 1, 0, 0 ], fromList [0]) ]  -- not a three

main :: IO ()
main = do
  n <- createNetwork 25 [250] 1
  let n' = trainNTimes 10000 0.5 tanh tanh' n samples
  mapM_ (print . output n' tanh . fst) samples
  putStrLn "-------------"
  print . output n' tanh $ testInput

  where testInput = fromList [ 0, 0, 1, 1, 1
                             , 0, 0, 0, 0, 1
                             , 0, 0, 1, 1, 1
                             , 0, 0, 0, 0, 1
                             , 0, 0, 1, 1, 1 ]

{-

OUTPUT:
fromList [0.9996325368507625]
fromList [0.9997784075859734]
fromList [0.9996165887689248]
fromList [-2.8107935971909852e-2]
fromList [7.001808876464477e-3]
fromList [2.54989546107178e-2]
fromList [5.286805464313172e-4]
-------------
fromList [0.9993713524712442]
-}

max_stars_repo_name: EditResearch/HsANN

max_issues_repo_name: EditResearch/HsANN

max_forks_repo_name: EditResearch/HsANN