text:
-- |
-- Module      : BattleHack.Utilities.Vector
-- Description :
-- Copyright   : (c) Jonatan H Sundqvist, 2015
-- License     : MIT
-- Maintainer  : Jonatan H Sundqvist
-- Stability   : experimental|stable
-- Portability : POSIX (not sure)
--

-- Created September 14 2015

-- TODO | - Is it silly to use Complex Double for everything because of Cairo (yes it is; fixed) (✓)
--        - Many of these functions should probably be moved to Southpaw

-- SPEC | -
--        -

--------------------------------------------------------------------------------------------------------------------------------------------
-- GHC Pragmas
--------------------------------------------------------------------------------------------------------------------------------------------

--------------------------------------------------------------------------------------------------------------------------------------------
-- API
--------------------------------------------------------------------------------------------------------------------------------------------
module BattleHack.Utilities.Vector where

--------------------------------------------------------------------------------------------------------------------------------------------
-- We'll need these
--------------------------------------------------------------------------------------------------------------------------------------------
import Data.Complex

import BattleHack.Types
import BattleHack.Lenses

--------------------------------------------------------------------------------------------------------------------------------------------
-- Functions
--------------------------------------------------------------------------------------------------------------------------------------------
-- |
dotwise :: (a -> a -> a) -> Complex a -> Complex a -> Complex a
dotwise f (x:+y) (x':+y') = f x x':+f y y'

-- |
dotmap :: (a -> a) -> Complex a -> Complex a
dotmap f (x:+y) = f x:+f y

-- |
toscalar :: (a -> a -> a) -> Complex a -> a
toscalar f (x:+y) = f x y

--  | Converts a 2-tuple to a vector
tovector :: (a, a) -> Complex a
tovector =  uncurry (:+)

-- |
vectorise :: RealFloat f => (f -> f -> a) -> Complex f -> a
vectorise f (x:+y) = f x y


text:
{-|
Module: Numeric.Morpheus.Statistics
Description: Statistics Functions
Copyright: (c) Alexander Ignatyev, 2017
License: BSD-3
Stability: experimental
Portability: POSIX
-}

module Numeric.Morpheus.Statistics
(
  mean
  , stddev_m
  , stddev
  , columnMean
  , rowMean
  , columnStddev_m
  , rowStddev_m
)

where

import Numeric.Morpheus.Utils(morpheusLayout)

import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Devel
import System.IO.Unsafe(unsafePerformIO)
import Foreign.C.Types
import Foreign.Ptr(Ptr)
import Foreign.Storable

{- morpheus_mean -}
foreign import ccall unsafe "morpheus_mean"
  c_morpheus_mean :: CInt -> Ptr R -> IO R

-- | Calculates mean (average).
mean :: Vector R -> R
mean x = unsafePerformIO $ do
    apply x id c_morpheus_mean

{- morpheus_stddev_m -}
foreign import ccall unsafe "morpheus_stddev"
  c_morpheus_stddev_m :: R -> CInt -> Ptr R -> IO R

-- | Calculates sample standard deviation using given mean value.
stddev_m :: R -> Vector R -> R
stddev_m m x = unsafePerformIO $ do
    apply x id (c_morpheus_stddev_m m)

{- morpheus_stddev -}
foreign import ccall unsafe "morpheus_stddev"
  c_morpheus_stddev :: CInt -> Ptr R -> IO R

-- | Calculates sample standard deviation.
stddev :: Vector R -> R
stddev x = unsafePerformIO $ do
    apply x id c_morpheus_stddev

{- morpheus_column_mean -}
foreign import ccall unsafe "morpheus_column_mean"
  c_morpheus_column_mean :: CInt -> CInt -> CInt -> Ptr R -> Ptr R -> IO ()

call_morpheus_column_mean :: CInt -> CInt -> CInt -> CInt -> Ptr R
                          -> CInt -> Ptr R
                          -> IO ()
call_morpheus_column_mean rows cols xRow xCol mat _ means = do
  let layout = morpheusLayout xCol cols
  c_morpheus_column_mean layout rows cols mat means

-- | Calculates mean (average) for every column.
columnMean :: Matrix R -> Vector R
columnMean m = unsafePerformIO $ do
    v <- createVector (cols m)
    apply m (apply v id) call_morpheus_column_mean
    return v

{- morpheus_row_mean -}
foreign import ccall unsafe "morpheus_row_mean"
  c_morpheus_row_mean :: CInt -> CInt -> CInt -> Ptr R -> Ptr R -> IO ()

call_morpheus_row_mean :: CInt -> CInt -> CInt -> CInt -> Ptr R
                       -> CInt -> Ptr R
                       -> IO ()
call_morpheus_row_mean rows cols xRow xCol mat _ means = do
  let layout = morpheusLayout xCol cols
  c_morpheus_row_mean layout rows cols mat means

-- | Calculates mean (average) for every row.
rowMean :: Matrix R -> Vector R
rowMean m = unsafePerformIO $ do
    v <- createVector (rows m)
    apply m (apply v id) call_morpheus_row_mean
    return v

{- morpheus_column_stddev_m -}
foreign import ccall unsafe "morpheus_column_stddev_m"
  c_morpheus_column_stddev_m :: CInt -> CInt -> CInt -> Ptr R -> Ptr R -> Ptr R -> IO ()

call_morpheus_column_stddev_m :: CInt -> Ptr R
                              -> CInt -> CInt -> CInt -> CInt -> Ptr R
                              -> CInt -> Ptr R
                              -> IO ()
call_morpheus_column_stddev_m _ means rows cols xRow xCol mat _ stddevs = do
  let layout = morpheusLayout xCol cols
  c_morpheus_column_stddev_m layout rows cols means mat stddevs

-- | Calculates sample standard deviation using given mean value for every column.
columnStddev_m :: Vector R -> Matrix R -> Vector R
columnStddev_m means m = unsafePerformIO $ do
    v <- createVector (cols m)
    apply means (apply m (apply v id)) call_morpheus_column_stddev_m
    return v

{- morpheus_row_stddev_m -}
foreign import ccall unsafe "morpheus_row_stddev_m"
  c_morpheus_row_stddev_m :: CInt -> CInt -> CInt -> Ptr R -> Ptr R -> Ptr R -> IO ()

call_morpheus_row_stddev_m :: CInt -> Ptr R
                              -> CInt -> CInt -> CInt -> CInt -> Ptr R
                              -> CInt -> Ptr R
                              -> IO ()
call_morpheus_row_stddev_m _ means rows cols xRow xCol mat _ stddevs = do
  let layout = morpheusLayout xCol cols
  c_morpheus_row_stddev_m layout rows cols means mat stddevs

-- | Calculates sample standard deviation using given mean value for every row.
rowStddev_m :: Vector R -> Matrix R -> Vector R
rowStddev_m means m = unsafePerformIO $ do
    v <- createVector (rows m)
    apply means (apply m (apply v id)) call_morpheus_row_stddev_m
    return v

ext: hs

max_stars_repo_name: Alexander-Ignatyev/morpheus

max_issues_repo_name: aligusnet/morpheus

max_forks_repo_name: aligusnet/morpheus


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

module Numerical.HBLAS.BLAS.FFI  where

import Foreign.Ptr
import Foreign()
import Foreign.C.Types
import Data.Complex

-- /*Set the number of threads on runtime.*/
--foreign import ccall unsafe "openblas_set_num_threads" openblas_set_num_threads_unsafe :: CInt -> IO ()

--foreign import ccall unsafe "goto_set_num_threads" goto_set_num_threads_unsafe :: CInt -> IO ()

{- | For All of the BlAS FFI operations,

-}
{-

typedef enum CBLAS_ORDER     {CblasRowMajor=101, CblasColMajor=102} CBLAS_ORDER;
typedef enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113, CblasConjNoTrans=114} CBLAS_TRANSPOSE;
typedef enum CBLAS_UPLO      {CblasUpper=121, CblasLower=122} CBLAS_UPLO;
typedef enum CBLAS_DIAG      {CblasNonUnit=131, CblasUnit=132} CBLAS_DIAG;
typedef enum CBLAS_SIDE      {CblasLeft=141, CblasRight=142} CBLAS_SIDE;

-}

newtype CBLAS_INDEX = CBIndex CSize
        deriving (Eq,Show)
newtype CBLAS_ORDERT = CBOInt CInt
    deriving (Eq,Show)
data BLASOrder = BLASRowMajor | BLASColMajor
    deriving (Eq,Show)

encodeOrder :: BLASOrder -> CBLAS_ORDERT
encodeOrder BLASRowMajor = CBOInt 101
encodeOrder BLASColMajor = CBOInt 102

newtype CBLAS_TRANSPOSET = CBLAS_TransposeT{ unCBLAS_TransposeT :: CInt } deriving (Eq, Show)

data BLAS_Transpose = BlasNoTranspose | BlasTranspose | BlasConjTranspose | BlasConjNoTranspose

encodeTranspose :: BLAS_Transpose -> CBLAS_TRANSPOSET
encodeTranspose  BlasNoTranspose = CBLAS_TransposeT 111
encodeTranspose  BlasTranspose = CBLAS_TransposeT 112
encodeTranspose  BlasConjTranspose =  CBLAS_TransposeT 113
encodeTranspose  BlasConjNoTranspose = CBLAS_TransposeT 114

newtype CBLAS_UPLOT = CBlasUPLO CInt
    deriving (Eq,Show)
data BLASUplo = BUpper | BLower
    deriving (Eq,Show)

encodeUPLO :: BLASUplo -> CBLAS_UPLOT
encodeUPLO BUpper = CBlasUPLO 121
encodeUPLO BLower = CBlasUPLO 122

newtype CBLAS_DIAGT = CBLAS_DiagT CUChar
    deriving (Show,Eq)

data BlasDiag = BlasNonUnit   | BlasUnit
    deriving (Eq,Show )

encodeDiag :: BlasDiag -> CBLAS_DIAGT
encodeDiag BlasNonUnit = CBLAS_DiagT 131
encodeDiag BlasUnit = CBLAS_DiagT 132

newtype CBLAS_SIDET = CBLAS_SideT { unCBLAS_SideT :: CUChar }
    deriving (Eq, Show)
data BlasSide = BlasLeft | BlasRight
    deriving (Eq,Show)
encodeSide :: BlasSide -> CBLAS_SIDET
encodeSide BlasLeft = CBLAS_SideT 141
encodeSide BlasRight = CBLAS_SideT 142

--typedef enum CBLAS_ORDER     {CblasRowMajor=101, CblasColMajor=102} CBLAS_ORDER;
--typedef enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113, CblasConjNoTrans=114} CBLAS_TRANSPOSE;
--typedef enum CBLAS_UPLO      {CblasUpper=121, CblasLower=122} CBLAS_UPLO;
--typedef enum CBLAS_DIAG      {CblasNonUnit=131, CblasUnit=132} CBLAS_DIAG;
--typedef enum CBLAS_SIDE      {CblasLeft=141, CblasRight=142} CBLAS_SIDE;

--dot products
foreign import ccall unsafe "cblas_sdsdot" cblas_sdsdot_unsafe :: CInt -> Float -> Ptr Float -> CInt -> Ptr Float -> CInt -> IO Float
foreign import ccall unsafe "cblas_dsdot" cblas_dsdot_unsafe :: CInt -> Ptr Float -> CInt -> Ptr Float -> CInt -> IO Double
foreign import ccall unsafe "cblas_sdot" cblas_sdot_unsafe :: CInt -> Ptr Float -> CInt -> Ptr Float -> CInt -> IO Float
foreign import ccall unsafe "cblas_ddot" cblas_ddot_unsafe :: CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> IO Double
--Float  cblas_sdsdot(  CInt n,   Float alpha,   Float *x,   CInt incx,   Float *y,   CInt incy);
--Double cblas_dsdot (  CInt n,   Float *x,   CInt incx,   Float *y,   CInt incy);
--Float  cblas_sdot(  CInt n,   Float  *x,   CInt incx,   Float  *y,   CInt incy);
--Double cblas_ddot(  CInt n,   Double *x,   CInt incx,   Double *y,   CInt incy);

{-
not doing these right now, because requires handling return value as a complex number,
we can only handle pointers to complex numbers right now
-}
--openblas_complex_Float  cblas_cdotu(  CInt n,   Float  *x,   CInt incx,   Float  *y,   CInt incy);
--openblas_complex_Float  cblas_cdotc(  CInt n,   Float  *x,   CInt incx,   Float  *y,   CInt incy);
--openblas_complex_Double cblas_zdotu(  CInt n,   Double *x,   CInt incx,   Double *y,   CInt incy);
--openblas_complex_Double cblas_zdotc(  CInt n,   Double *x,   CInt incx,   Double *y,   CInt incy);

--- not sure what to do for these complex
--void  cblas_cdotu_sub(  CInt n,   Float  *x,   CInt incx,   Float  *y,   CInt incy, openblas_complex_Float  *ret);
--void  cblas_cdotc_sub(  CInt n,   Float  *x,   CInt incx,   Float  *y,   CInt incy, openblas_complex_Float  *ret);
--void  cblas_zdotu_sub(  CInt n,   Double *x,   CInt incx,   Double *y,   CInt incy, openblas_complex_Double *ret);
--void  cblas_zdotc_sub(  CInt n,   Double *x,   CInt incx,   Double *y,   CInt incy, openblas_complex_Double *ret);

---- absolute value
foreign import ccall unsafe "cblas_sasum" cblas_sasum_unsafe::
    CInt -> Ptr Float -> CInt -> IO Float
foreign import ccall unsafe "cblas_dasum" cblas_dasum_unsafe ::
    CInt -> Ptr Double -> CInt -> IO Double
foreign import ccall unsafe "cblas_scasum" cblas_casum_unsafe ::
    CInt -> Ptr (Complex Float)-> CInt -> IO Float
foreign import ccall unsafe "cblas_dzasum" cblas_zasum_unsafe ::
    CInt -> Ptr (Complex Double) -> CInt -> IO Double
--Float  cblas_sasum (  CInt n,   Float  *x,   CInt incx);
--Double cblas_dasum (  CInt n,   Double *x,   CInt incx);
--Float  cblas_scasum(  CInt n,   Float  *x,   CInt incx);
--Double cblas_dzasum(  CInt n,   Double *x,   CInt incx);

foreign import ccall unsafe "cblas_snrm2" cblas_snrm2_unsafe ::
    CInt -> Ptr Float -> CInt -> IO Float
foreign import ccall unsafe "cblas_dnrm2" cblas_dnrm2_unsafe ::
    CInt -> Ptr Double -> CInt -> IO Double
foreign import ccall unsafe "cblas_scnrm2" cblas_scnrm2_unsafe ::
    CInt -> Ptr (Complex Float)-> CInt -> IO Float
foreign import ccall unsafe "cblas_dznrm2" cblas_dznrm2_unsafe ::
    CInt -> Ptr (Complex Double) -> CInt -> IO Double

--Float  cblas_snrm2 (  CInt N,   Float  *X,   CInt incX);
--Double cblas_dnrm2 (  CInt N,   Double *X,   CInt incX);
--Float  cblas_scnrm2(  CInt N,   Float  *X,   CInt incX);
--Double cblas_dznrm2(  CInt N,   Double *X,   CInt incX);

foreign import ccall unsafe "cblas_isamax" cblas_isamax_unsafe ::
    CInt -> Ptr Float -> CInt -> IO CInt
foreign import ccall unsafe "cblas_idamax" cblas_idamax_unsafe ::
    CInt -> Ptr Float -> CInt -> IO CInt
foreign import ccall unsafe "cblas_icamax" cblas_icamax_unsafe ::
    CInt -> Ptr (Complex Float) -> CInt -> IO CInt
foreign import ccall unsafe "cblas_izamax" cblas_izamax_unsafe ::
    CInt -> Ptr (Complex Double) -> CInt -> IO CInt
--CBLAS_INDEX cblas_isamax(  CInt n,   Float  *x,   CInt incx);
--CBLAS_INDEX cblas_idamax(  CInt n,   Double *x,   CInt incx);
--CBLAS_INDEX cblas_icamax(  CInt n,   Float  *x,   CInt incx);
--CBLAS_INDEX cblas_izamax(  CInt n,   Double *x,   CInt incx);

--void cblas_saxpy(  CInt n,   Float alpha,   Float *x,   CInt incx, Float *y,   CInt incy);
--void cblas_daxpy(  CInt n,   Double alpha,   Double *x,   CInt incx, Double *y,   CInt incy);
--void cblas_caxpy(  CInt n,   Float *alpha,   Float *x,   CInt incx, Float *y,   CInt incy);
--void cblas_zaxpy(  CInt n,   Double *alpha,   Double *x,   CInt incx, Double *y,   CInt incy);

foreign import ccall unsafe "cblas_scopy" cblas_scopy_unsafe ::
    CInt -> Ptr Float -> CInt -> Ptr Float -> CInt -> IO ()
foreign import ccall unsafe "cblas_dcopy" cblas_dcopy_unsafe ::
    CInt -> Ptr Double-> CInt -> Ptr Double -> CInt -> IO ()
foreign import ccall unsafe "cblas_ccopy" cblas_ccopy_unsafe ::
    CInt -> Ptr (Complex Float) -> CInt -> Ptr (Complex Float) -> CInt -> IO ()
foreign import ccall unsafe "cblas_zcopy" cblas_zcopy_unsafe ::
    CInt -> Ptr (Complex Double)-> CInt -> Ptr (Complex Double) -> CInt -> IO ()

--void cblas_scopy(  CInt n,   Float *x,   CInt incx, Float *y,   CInt incy);
--void cblas_dcopy(  CInt n,   Double *x,   CInt incx, Double *y,   CInt incy);
--void cblas_ccopy(  CInt n,   Float *x,   CInt incx, Float *y,   CInt incy);
--void cblas_zcopy(  CInt n,   Double *x,   CInt incx, Double *y,   CInt incy);

foreign import ccall unsafe "cblas_sswap" cblas_sswap_unsafe ::
    CInt -> Ptr Float -> CInt -> Ptr Float -> CInt -> IO ()
foreign import ccall unsafe "cblas_dswap" cblas_dswap_unsafe ::
    CInt -> Ptr Double-> CInt -> Ptr Double -> CInt -> IO ()
foreign import ccall unsafe "cblas_cswap" cblas_cswap_unsafe ::
    CInt -> Ptr (Complex Float) -> CInt -> Ptr (Complex Float) -> CInt -> IO ()
foreign import ccall unsafe "cblas_zswap" cblas_zswap_unsafe ::
    CInt -> Ptr (Complex Double)-> CInt -> Ptr (Complex Double) -> CInt -> IO ()

--void cblas_sswap(  CInt n, Float *x,   CInt incx, Float *y,   CInt incy);
--void cblas_dswap(  CInt n, Double *x,   CInt incx, Double *y,   CInt incy);
--void cblas_cswap(  CInt n, Float *x,   CInt incx, Float *y,   CInt incy);
--void cblas_zswap(  CInt n, Double *x,   CInt incx, Double *y,   CInt incy);

--void cblas_srot(  CInt N, Float *X,   CInt incX, Float *Y,   CInt incY,   Float c,   Float s);
--void cblas_drot(  CInt N, Double *X,   CInt incX, Double *Y,   CInt incY,   Double c,   Double  s);

--void cblas_srotg(Float *a, Float *b, Float *c, Float *s);
--void cblas_drotg(Double *a, Double *b, Double *c, Double *s);

--void cblas_srotm(  CInt N, Float *X,   CInt incX, Float *Y,   CInt incY,   Float *P);
--void cblas_drotm(  CInt N, Double *X,   CInt incX, Double *Y,   CInt incY,   Double *P);

--void cblas_srotmg(Float *d1, Float *d2, Float *b1,   Float b2, Float *P);
--void cblas_drotmg(Double *d1, Double *d2, Double *b1,   Double b2, Double *P);

--void cblas_sscal(  CInt N,   Float alpha, Float *X,   CInt incX);
--void cblas_dscal(  CInt N,   Double alpha, Double *X,   CInt incX);
--void cblas_cscal(  CInt N,   Float *alpha, Float *X,   CInt incX);
--void cblas_zscal(  CInt N,   Double *alpha, Double *X,   CInt incX);
--void cblas_csscal(  CInt N,   Float alpha, Float *X,   CInt incX);
--void cblas_zdscal(  CInt N,   Double alpha, Double *X,   CInt incX);

---------------------
----- BLAS LEVEL 2
---------------------

{-
matrix vector product for general matrices
 perform one of the matrix-vector operations   y :=
      alpha*A*x + beta*y, or y := alpha*A'*x + beta*y,
-}

type GemvFunFFI sc el =
       CBLAS_ORDERT -> CBLAS_TRANSPOSET -> CInt -> CInt
    -> sc -> Ptr el -> CInt -> Ptr el -> CInt -> sc -> Ptr el -> CInt -> IO ()

foreign import ccall unsafe "cblas_sgemv"
    cblas_sgemv_unsafe :: GemvFunFFI Float Float
foreign import ccall safe   "cblas_sgemv"
    cblas_sgemv_safe   :: GemvFunFFI Float Float

foreign import ccall unsafe "cblas_dgemv"
    cblas_dgemv_unsafe :: GemvFunFFI Double Double
foreign import ccall safe   "cblas_dgemv"
    cblas_dgemv_safe   :: GemvFunFFI Double Double

foreign import ccall unsafe "cblas_cgemv"
    cblas_cgemv_unsafe :: GemvFunFFI (Ptr (Complex Float)) (Complex Float)
foreign import ccall safe   "cblas_cgemv"
    cblas_cgemv_safe   :: GemvFunFFI (Ptr (Complex Float)) (Complex Float)

foreign import ccall unsafe "cblas_zgemv"
    cblas_zgemv_unsafe :: GemvFunFFI (Ptr (Complex Double)) (Complex Double)
foreign import ccall safe   "cblas_zgemv"
    cblas_zgemv_safe   :: GemvFunFFI (Ptr (Complex Double)) (Complex Double)

--void cblas_sgemv(  enum CBLAS_ORDER order,    enum CBLAS_TRANSPOSE trans,    CInt m,   CInt n,
--           Float alpha,   Float  *a,   CInt lda,    Float  *x,   CInt incx,    Float beta,  Float  *y,   CInt incy);
--void cblas_dgemv(  enum CBLAS_ORDER order,    enum CBLAS_TRANSPOSE trans,    CInt m,   CInt n,
--           Double alpha,   Double  *a,   CInt lda,    Double  *x,   CInt incx,    Double beta,  Double  *y,   CInt incy);
--void cblas_cgemv(  enum CBLAS_ORDER order,    enum CBLAS_TRANSPOSE trans,    CInt m,   CInt n,
--           Float *alpha,   Float  *a,   CInt lda,    Float  *x,   CInt incx,    Float *beta,  Float  *y,   CInt incy);
--void cblas_zgemv(  enum CBLAS_ORDER order,    enum CBLAS_TRANSPOSE trans,    CInt m,   CInt n,
--           Double *alpha,   Double  *a,   CInt lda,    Double  *x,   CInt incx,    Double *beta,  Double  *y,   CInt incy);

-- perform the rank 1 operation   A := alpha*x*y' + A,

type GerFunFFI el = CBLAS_ORDERT -> CInt -> CInt -> el -> Ptr el -> CInt -> Ptr el -> CInt -> Ptr el -> CInt -> IO ()

foreign import ccall unsafe "cblas_sger" cblas_sger_unsafe ::
        CBLAS_ORDERT -> CInt -> CInt -> Float -> Ptr Float -> CInt -> Ptr Float -> CInt -> Ptr Float -> CInt -> IO ()
foreign import ccall safe   "cblas_sger" cblas_sger_safe   ::
        CBLAS_ORDERT -> CInt -> CInt -> Float -> Ptr Float -> CInt -> Ptr Float -> CInt -> Ptr Float -> CInt -> IO ()
foreign import ccall unsafe "cblas_dger" cblas_dger_unsafe ::
        CBLAS_ORDERT -> CInt -> CInt -> Double -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> IO ()
foreign import ccall safe   "cblas_dger" cblas_dger_safe   ::
        CBLAS_ORDERT -> CInt -> CInt -> Double -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> IO ()

--void cblas_sger (  enum CBLAS_ORDER order,   CInt M,   CInt N,   Float   alpha,   Float  *X,   CInt incX,   Float  *Y,   CInt incY, Float  *A,   CInt lda);
--void cblas_dger (  enum CBLAS_ORDER order,   CInt M,   CInt N,   Double  alpha,   Double *X,   CInt incX,   Double *Y,   CInt incY, Double *A,   CInt lda);
--void cblas_cgeru(  enum CBLAS_ORDER order,   CInt M,   CInt N,   Float  *alpha,   Float  *X,   CInt incX,   Float  *Y,   CInt incY, Float  *A,   CInt lda);
--void cblas_cgerc(  enum CBLAS_ORDER order,   CInt M,   CInt N,   Float  *alpha,   Float  *X,   CInt incX,   Float  *Y,   CInt incY, Float  *A,   CInt lda);
--void cblas_zgeru(  enum CBLAS_ORDER order,   CInt M,   CInt N,   Double *alpha,   Double *X,   CInt incX,   Double *Y,   CInt incY, Double *A,   CInt lda);
--void cblas_zgerc(  enum CBLAS_ORDER order,   CInt M,   CInt N,   Double *alpha,   Double *X,   CInt incX,   Double *Y,   CInt incY, Double *A,   CInt lda);

--STRSV - solve one of the systems of equations   A*x = b, or A'*x = b, where A is a (non)unit upper(/lower) triangular matrix

--STRSV - solve one of the systems of equations   A*x = b, or A'*x = b, where A is a (non)unit upper(/lower) triangular matrix

type TrsvFunFFI el =
       CBLAS_ORDERT -> CBLAS_UPLOT -> CBLAS_TRANSPOSET -> CBLAS_DIAGT
    -> CInt -> Ptr el -> CInt -> Ptr el -> CInt -> IO ()

foreign import ccall unsafe "cblas_strsv"
  cblas_strsv_unsafe :: TrsvFunFFI Float
foreign import ccall safe   "cblas_strsv"
  cblas_strsv_safe   :: TrsvFunFFI Float

foreign import ccall unsafe "cblas_dtrsv"
  cblas_dtrsv_unsafe :: TrsvFunFFI Double
foreign import ccall safe   "cblas_dtrsv"
  cblas_dtrsv_safe   :: TrsvFunFFI Double

foreign import ccall unsafe "cblas_ctrsv"
  cblas_ctrsv_unsafe :: TrsvFunFFI (Complex Float)
foreign import ccall safe   "cblas_ctrsv"
  cblas_ctrsv_safe   :: TrsvFunFFI (Complex Float)

foreign import ccall unsafe "cblas_ztrsv"
  cblas_ztrsv_unsafe :: TrsvFunFFI (Complex Double)
foreign import ccall safe   "cblas_ztrsv"
  cblas_ztrsv_safe   :: TrsvFunFFI (Complex Double)

--void cblas_strsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,   CInt N,   Float *A,   CInt lda, Float *X,   CInt incX);
--void cblas_dtrsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,   CInt N,   Double *A,   CInt lda, Double *X,   CInt incX);
--void cblas_ctrsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,   CInt N,   Float *A,   CInt lda, Float *X,   CInt incX);
--void cblas_ztrsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,   CInt N,   Double *A,   CInt lda, Double *X,   CInt incX);

--void cblas_strmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,   CInt N,   Float *A,   CInt lda, Float *X,   CInt incX);
--void cblas_dtrmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,   CInt N,   Double *A,   CInt lda, Double *X,   CInt incX);
--void cblas_ctrmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,   CInt N,   Float *A,   CInt lda, Float *X,   CInt incX);
--void cblas_ztrmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,   CInt N,   Double *A,   CInt lda, Double *X,   CInt incX);

--void cblas_ssyr(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Float alpha,   Float *X,   CInt incX, Float *A,   CInt lda);
--void cblas_dsyr(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double alpha,   Double *X,   CInt incX, Double *A,   CInt lda);
--void cblas_cher(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Float alpha,   Float *X,   CInt incX, Float *A,   CInt lda);
--void cblas_zher(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double alpha,   Double *X,   CInt incX, Double *A,   CInt lda);

--void cblas_ssyr2(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,  CInt N,   Float alpha,   Float *X,
--                  CInt incX,   Float *Y,   CInt incY, Float *A,   CInt lda);
--void cblas_dsyr2(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double alpha,   Double *X,
--                  CInt incX,   Double *Y,   CInt incY, Double *A,   CInt lda);
--void cblas_cher2(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Float *alpha,   Float *X,   CInt incX,
--                  Float *Y,   CInt incY, Float *A,   CInt lda);
--void cblas_zher2(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double *alpha,   Double *X,   CInt incX,
--                  Double *Y,   CInt incY, Double *A,   CInt lda);

--void cblas_sgbmv(  enum CBLAS_ORDER order,   enum CBLAS_TRANSPOSE TransA,   CInt M,   CInt N,
--                   CInt KL,   CInt KU,   Float alpha,   Float *A,   CInt lda,   Float *X,   CInt incX,   Float beta, Float *Y,   CInt incY);
--void cblas_dgbmv(  enum CBLAS_ORDER order,   enum CBLAS_TRANSPOSE TransA,   CInt M,   CInt N,
--                   CInt KL,   CInt KU,   Double alpha,   Double *A,   CInt lda,   Double *X,   CInt incX,   Double beta, Double *Y,   CInt incY);
--void cblas_cgbmv(  enum CBLAS_ORDER order,   enum CBLAS_TRANSPOSE TransA,   CInt M,   CInt N,
--                   CInt KL,   CInt KU,   Float *alpha,   Float *A,   CInt lda,   Float *X,   CInt incX,   Float *beta, Float *Y,   CInt incY);
--void cblas_zgbmv(  enum CBLAS_ORDER order,   enum CBLAS_TRANSPOSE TransA,   CInt M,   CInt N,
--                   CInt KL,   CInt KU,   Double *alpha,   Double *A,   CInt lda,   Double *X,   CInt incX,   Double *beta, Double *Y,   CInt incY);

--void cblas_ssbmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   CInt K,   Float alpha,   Float *A,
--                   CInt lda,   Float *X,   CInt incX,   Float beta, Float *Y,   CInt incY);
--void cblas_dsbmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   CInt K,   Double alpha,   Double *A,
--                   CInt lda,   Double *X,   CInt incX,   Double beta, Double *Y,   CInt incY);

--void cblas_stbmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   CInt K,   Float *A,   CInt lda, Float *X,   CInt incX);
--void cblas_dtbmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   CInt K,   Double *A,   CInt lda, Double *X,   CInt incX);
--void cblas_ctbmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   CInt K,   Float *A,   CInt lda, Float *X,   CInt incX);
--void cblas_ztbmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   CInt K,   Double *A,   CInt lda, Double *X,   CInt incX);

----------------
--- | solves  Ax=v where A is k+1 banded triangular matrix, and x and
----------------
--void cblas_stbsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   CInt K,   Float *A,   CInt lda, Float *X,   CInt incX);
--void cblas_dtbsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   CInt K,   Double *A,   CInt lda, Double *X,   CInt incX);
--void cblas_ctbsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   CInt K,   Float *A,   CInt lda, Float *X,   CInt incX);
--void cblas_ztbsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   CInt K,   Double *A,   CInt lda, Double *X,   CInt incX);

-------------------------------------------------------------------------
-- | matrix vector product Av, writes result into v, where A is a packed triangular nxn matrix
-------------------------------------------------------------------------
--void cblas_stpmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   Float *Ap, Float *X,   CInt incX);
--void cblas_dtpmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   Double *Ap, Double *X,   CInt incX);
--void cblas_ctpmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   Float *Ap, Float *X,   CInt incX);
--void cblas_ztpmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   Double *Ap, Double *X,   CInt incX);

--------------------------------------------------
---  | solve  Ax=v where A is a nxn packed triangular matrix, v vector input, writes the solution into x.
--------------------------------------------------
--void cblas_stpsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   Float *Ap, Float *X,   CInt incX);
--void cblas_dtpsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   Double *Ap, Double *X,   CInt incX);
--void cblas_ctpsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   Float *Ap, Float *X,   CInt incX);
--void cblas_ztpsv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   enum CBLAS_TRANSPOSE TransA,   enum CBLAS_DIAG Diag,
--                   CInt N,   Double *Ap, Double *X,   CInt incX);

----------------------------------
---- |  (unpacked) symmetric matrix vector product    x:=Av, writes result x into v
---------------------------------

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

--void cblas_ssymv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Float alpha,   Float *A,
--                   CInt lda,   Float *X,   CInt incX,   Float beta, Float *Y,   CInt incY);
--void cblas_dsymv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double alpha,   Double *A,
--                   CInt lda,   Double *X,   CInt incX,   Double beta, Double *Y,   CInt incY);

--------------------------------
---- | hermitian matrix vector product   x:=Av, writes result x into v
--------------------------------
--void cblas_chemv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Float *alpha,   Float *A,
--                   CInt lda,   Float *X,   CInt incX,   Float *beta, Float *Y,   CInt incY);
--void cblas_zhemv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double *alpha,   Double *A,
--                   CInt lda,   Double *X,   CInt incX,   Double *beta, Double *Y,   CInt incY);

---------------
--- | packed symmetric matrix * vector product  y:= alpha * Av  + beta * y
---------------

--void cblas_sspmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Float alpha,   Float *Ap,
--                   Float *X,   CInt incX,   Float beta, Float *Y,   CInt incY);
--void cblas_dspmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double alpha,   Double *Ap,
--                   Double *X,   CInt incX,   Double beta, Double *Y,   CInt incY);

--void cblas_sspr(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Float alpha,   Float *X,   CInt incX, Float *Ap);
--void cblas_dspr(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double alpha,   Double *X,   CInt incX, Double *Ap);

--void cblas_chpr(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Float alpha,   Float *X,   CInt incX, Float *A);
--void cblas_zhpr(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double alpha,   Double *X,  CInt incX, Double *A);

--void cblas_sspr2(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Float alpha,   Float *X,   CInt incX,   Float *Y,   CInt incY, Float *A);
--void cblas_dspr2(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double alpha,   Double *X,   CInt incX,   Double *Y,   CInt incY, Double *A);
--void cblas_chpr2(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Float *alpha,   Float *X,   CInt incX,   Float *Y,   CInt incY, Float *Ap);
--void cblas_zhpr2(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   Double *alpha,   Double *X,   CInt incX,   Double *Y,   CInt incY, Double *Ap);

--void cblas_chbmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   CInt K,
--           Float *alpha,   Float *A,   CInt lda,   Float *X,   CInt incX,   Float *beta, Float *Y,   CInt incY);
--void cblas_zhbmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,   CInt K,
--           Double *alpha,   Double *A,   CInt lda,   Double *X,   CInt incX,   Double *beta, Double *Y,   CInt incY);

--void cblas_chpmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,
--           Float *alpha,   Float *Ap,   Float *X,   CInt incX,   Float *beta, Float *Y,   CInt incY);
--void cblas_zhpmv(  enum CBLAS_ORDER order,   enum CBLAS_UPLO Uplo,   CInt N,
--           Double *alpha,   Double *Ap,   Double *X,   CInt incX,   Double *beta, Double *Y,   CInt incY);

{-

-}

--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
------------------------------ | 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 -> 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)

----------------------
----- | 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)) Float
foreign  import ccall unsafe "cblas_zsyr2k"
    cblas_zsyr2k_unsafe :: Syr2kFunFFI (Ptr (Complex Double)) 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 -> 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)

--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 -> Ptr el -> CInt -> IO ()
foreign  import ccall unsafe "cblas_strsm"
    cblas_strsm_unsafe :: TrmmFunFFI Float Float
foreign  import ccall unsafe "cblas_dtrsm"
    cblas_dtrsm_unsafe :: TrmmFunFFI Double Double
foreign  import ccall unsafe "cblas_ctrsm"
    cblas_ctrsm_unsafe :: TrmmFunFFI (Ptr (Complex Float )) (Complex Float)
foreign  import ccall unsafe "cblas_ztrsm"
    cblas_ztrsm_unsafe :: TrmmFunFFI (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 -> 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)

--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 HerkFun scale el = CBLAS_ORDERT -> CBLAS_SIDET-> CBLAS_TRANSPOSET ->
     CInt->CInt ->  scale -> Ptr el -> CInt -> Ptr el -> CInt ->scale ->Ptr el -> CInt -> IO ()

foreign  import ccall unsafe "cblas_cherk"
    cblas_cherk_unsafe :: HerkFun  Float  (Complex Float)
foreign  import ccall unsafe "cblas_zherk"
    cblas_zherk_unsafe :: HerkFun  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_SIDET -> 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)

--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: archblob/hblas

max_issues_repo_name: archblob/hblas

max_forks_repo_name: archblob/hblas


text:
{-# LANGUAGE GADTs #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}

--
-- Copyright (c) 2009-2011, ERICSSON AB
-- All rights reserved.
--
-- Redistribution and use in source and binary forms, with or without
-- modification, are permitted provided that the following conditions are met:
--
--     * Redistributions of source code must retain the above copyright notice,
--       this list of conditions and the following disclaimer.
--     * Redistributions in binary form must reproduce the above copyright
--       notice, this list of conditions and the following disclaimer in the
--       documentation and/or other materials provided with the distribution.
--     * Neither the name of the ERICSSON AB nor the names of its contributors
--       may be used to endorse or promote products derived from this software
--       without specific prior written permission.
--
-- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-- AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
-- DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
-- FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
-- DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
-- SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
-- CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
-- OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
-- OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
--

module Feldspar.Core.Types where

import Data.Array.IO
import Data.Bits
import Data.Complex
import Data.Int
import Data.IORef
import Data.List
import Data.Typeable (Typeable, Typeable1)
import Data.Word
import Test.QuickCheck
import qualified Control.Monad.Par as MonadPar

import Data.Patch
import Data.Proxy

import Feldspar.Lattice
import Feldspar.Range

--------------------------------------------------------------------------------
-- * Heterogenous lists
--------------------------------------------------------------------------------

-- | Heterogeneous list
data a :> b = a :> b
  deriving (Eq, Ord, Show)

infixr 5 :>

instance (Lattice a, Lattice b) => Lattice (a :> b)
  where
    bot = bot :> bot
    top = top :> top
    (a1:>a2) \/ (b1:>b2) = (a1 \/ b1) :> (a2 \/ b2)
    (a1:>a2) /\ (b1:>b2) = (a1 /\ b1) :> (a2 /\ b2)

--------------------------------------------------------------------------------
-- * Integers
--------------------------------------------------------------------------------

-- | Target-dependent unsigned integers
newtype WordN = WordN Word32
  deriving
    ( Eq, Ord, Num, Enum, Ix, Real, Integral, Bits, Bounded, Typeable
    , Arbitrary )

type instance UnsignedRep WordN = Word32

-- | Target-dependent signed integers
newtype IntN = IntN Int32
  deriving
    ( Eq, Ord, Num, Enum, Ix, Real, Integral, Bits, Bounded, Typeable
    , Arbitrary )

type instance UnsignedRep IntN = Word32

instance Show WordN
  where
    show (WordN a) = show a

instance Show IntN
  where
    show (IntN a) = show a

-- | Type representation of 8 bits
data N8

-- | Type representation of 16 bits
data N16

-- | Type representation of 32 bits
data N32

-- | Type representation of 64 bits
data N64

-- | Type representation of the native number of bits on the target
data NNative

-- | Witness for 'N8', 'N16', 'N32', 'N64' or 'NNative'
data BitWidth n
  where
    N8      :: BitWidth N8
    N16     :: BitWidth N16
    N32     :: BitWidth N32
    N64     :: BitWidth N64
    NNative :: BitWidth NNative

bitWidth :: BitWidth n -> String
bitWidth N8      = "8"
bitWidth N16     = "16"
bitWidth N32     = "32"
bitWidth N64     = "64"
bitWidth NNative = "N"

-- | Type representation of \"unsigned\"
data U

-- | Type representation of \"signed\"
data S

-- | Witness for 'U' or 'S'
data Signedness s
  where
    U :: Signedness U
    S :: Signedness S

signedness :: Signedness s -> String
signedness U = "Word"
signedness S = "Int"

-- | A generalization of unsigned and signed integers. The first parameter
-- represents the signedness and the sectond parameter the number of bits.
type family GenericInt s n
type instance GenericInt U N8      = Word8
type instance GenericInt S N8      = Int8
type instance GenericInt U N16     = Word16
type instance GenericInt S N16     = Int16
type instance GenericInt U N32     = Word32
type instance GenericInt S N32     = Int32
type instance GenericInt U N64     = Word64
type instance GenericInt S N64     = Int64
type instance GenericInt U NNative = WordN
type instance GenericInt S NNative = IntN

type family WidthOf a
type instance WidthOf Word8  = N8
type instance WidthOf Int8   = N8
type instance WidthOf Word16 = N16
type instance WidthOf Int16  = N16
type instance WidthOf Word32 = N32
type instance WidthOf Int32  = N32
type instance WidthOf Word64 = N64
type instance WidthOf Int64  = N64
type instance WidthOf WordN  = NNative
type instance WidthOf IntN   = NNative

type family SignOf a
type instance SignOf Word8  = U
type instance SignOf Int8   = S
type instance SignOf Word16 = U
type instance SignOf Int16  = S
type instance SignOf Word32 = U
type instance SignOf Int32  = S
type instance SignOf Word64 = U
type instance SignOf Int64  = S
type instance SignOf WordN  = U
type instance SignOf IntN   = S

fromWordN :: BitWidth n -> WordN -> GenericInt U n
fromWordN N8      = fromInteger . toInteger
fromWordN N16     = fromInteger . toInteger
fromWordN N32     = fromInteger . toInteger
fromWordN N64     = fromInteger . toInteger
fromWordN NNative = id
  -- TODO Check that the number fits

fromIntN :: BitWidth n -> IntN -> GenericInt S n
fromIntN N8      = fromInteger . toInteger
fromIntN N16     = fromInteger . toInteger
fromIntN N32     = fromInteger . toInteger
fromIntN N64     = fromInteger . toInteger
fromIntN NNative = id
  -- TODO Check that the number fits

genericLen :: BitWidth n -> [a] -> GenericInt U n
genericLen N8      = genericLength
genericLen N16     = genericLength
genericLen N32     = genericLength
genericLen N64     = genericLength
genericLen NNative = genericLength

type Length = WordN
type Index  = WordN

--------------------------------------------------------------------------------
-- * Arrays
--------------------------------------------------------------------------------

-- | Array whose length is represented by an @n@-bit word
data TargetArr n a = TargetArr (GenericInt U n) [a]

--------------------------------------------------------------------------------
-- * Monadic Types
--------------------------------------------------------------------------------

-- | This class is used to allow constructs to be abstract in the monad
class MonadType m
  where
    voidTypeRep :: TypeRep (m ())
  -- TODO Since the `Mut` monad is already abstract, this class is probably only
  --      needed if we want to be able to use monadic constructs with different
  --      monad types.

--------------------------------------------------------------------------------
-- * Mutable data
--------------------------------------------------------------------------------

-- TODO Make newtypes?

-- | Monad for manipulation of mutable data
type Mut = IO

-- | Mutable references
instance Show (IORef a)
  where
    show _ = "IORef"

-- | Mutable arrays
type MArr a = IOArray Index a

instance Show (MArr a)
  where
    show _ = "MArr"

instance MonadType Mut
  where
    voidTypeRep = MutType UnitType

--------------------------------------------------------------------------------
-- * Par Monad
--------------------------------------------------------------------------------

-- | Monad for parallel constructs
type Par = MonadPar.Par

deriving instance Typeable1 Par

-- | Immutable references
type IV = MonadPar.IVar

deriving instance Typeable1 IV

instance Show (IV a)
  where
    show _ = "IVar"

instance Eq (IV a)
  where
    _ == _ = False

instance MonadType Par
  where
    voidTypeRep = ParType UnitType

--------------------------------------------------------------------------------
-- * Future values
--------------------------------------------------------------------------------

newtype FVal a = FVal {unFVal :: a}

deriving instance Typeable1 FVal

instance Show (FVal a)
  where
    show _ = "future"

instance Eq a => Eq (FVal a)
  where
    (FVal a) == (FVal b) = a == b

--------------------------------------------------------------------------------
-- * Type representation
--------------------------------------------------------------------------------

-- | Representation of supported types
data TypeRep a
  where
    UnitType      :: TypeRep ()
    BoolType      :: TypeRep Bool
    IntType       :: ( BoundedInt (GenericInt s n)
                     , Size (GenericInt s n) ~ Range (GenericInt s n)
                     ) =>
                       Signedness s -> BitWidth n -> TypeRep (GenericInt s n)
    FloatType     :: TypeRep Float
    DoubleType    :: TypeRep Double
    ComplexType   :: RealFloat a => TypeRep a -> TypeRep (Complex a)
    ArrayType     :: TypeRep a -> TypeRep [a]
    TargetArrType :: BitWidth n -> TypeRep a -> TypeRep (TargetArr n a)
    Tup2Type      :: TypeRep a -> TypeRep b -> TypeRep (a,b)
    Tup3Type      :: TypeRep a -> TypeRep b -> TypeRep c -> TypeRep (a,b,c)
    Tup4Type      :: TypeRep a -> TypeRep b -> TypeRep c -> TypeRep d -> TypeRep (a,b,c,d)
    Tup5Type      :: TypeRep a -> TypeRep b -> TypeRep c -> TypeRep d -> TypeRep e -> TypeRep (a,b,c,d,e)
    Tup6Type      :: TypeRep a -> TypeRep b -> TypeRep c -> TypeRep d -> TypeRep e -> TypeRep f -> TypeRep (a,b,c,d,e,f)
    Tup7Type      :: TypeRep a -> TypeRep b -> TypeRep c -> TypeRep d -> TypeRep e -> TypeRep f -> TypeRep g -> TypeRep (a,b,c,d,e,f,g)
    FunType       :: TypeRep a -> TypeRep b -> TypeRep (a -> b)
    MutType       :: TypeRep a -> TypeRep (Mut a)
    RefType       :: TypeRep a -> TypeRep (IORef a)
    MArrType      :: TypeRep a -> TypeRep (MArr a)
    ParType       :: TypeRep a -> TypeRep (Par a)
    IVarType      :: TypeRep a -> TypeRep (IV a)
    FValType      :: TypeRep a -> TypeRep (FVal a)
      -- TODO `MArrType` Should have a target-specialized version. Or perhaps
      --      use a single type with a flag to distinguish between immutable and
      --      mutable arrays.

instance Show (TypeRep a)
  where
    show UnitType            = "()"
    show BoolType            = "Bool"
    show (IntType s n)       = signedness s ++ bitWidth n
    show FloatType           = "Float"
    show DoubleType          = "Double"
    show (ComplexType t)     = "(Complex " ++ show t ++ ")"
    show (ArrayType t)       = "[" ++ show t ++ "]"
    show (TargetArrType _ t) = "[" ++ show t ++ "]"
    show (Tup2Type ta tb)                = showTup [show ta, show tb]
    show (Tup3Type ta tb tc)             = showTup [show ta, show tb, show tc]
    show (Tup4Type ta tb tc td)          = showTup [show ta, show tb, show tc, show td]
    show (Tup5Type ta tb tc td te)       = showTup [show ta, show tb, show tc, show td, show te]
    show (Tup6Type ta tb tc td te tf)    = showTup [show ta, show tb, show tc, show td, show te, show tf]
    show (Tup7Type ta tb tc td te tf tg) = showTup [show ta, show tb, show tc, show td, show te, show tf, show tg]
    show (FunType ta tb)                 = show ta ++ " -> " ++ show tb
    show (MutType ta)                    = unwords ["Mut", show ta]
    show (RefType ta)                    = unwords ["Ref", show ta]
    show (MArrType ta)                   = unwords ["MArr", show ta]
    show (ParType ta)                    = unwords ["Par", show ta]
    show (IVarType ta)                   = unwords ["IVar", show ta]
    show (FValType ta)                   = unwords ["FVal", show ta]

argType :: TypeRep (a -> b) -> TypeRep a
argType (FunType ta _) = ta

resType :: TypeRep (a -> b) -> TypeRep b
resType (FunType _ tb) = tb

-- | Type equality witness
data TypeEq a b
  where
    TypeEq :: TypeEq a a

defaultSize :: TypeRep a -> Size a
defaultSize UnitType = universal
defaultSize BoolType = universal
defaultSize (IntType _ _) = universal
defaultSize FloatType = universal
defaultSize DoubleType = universal
defaultSize (ComplexType _) = universal
defaultSize (ArrayType t) = universal :> defaultSize t
--defaultSize (TargetArrType n t) = universal :> defaultSize t -- TODO
defaultSize (Tup2Type ta tb) =  ( defaultSize ta
                                , defaultSize tb
                                )
defaultSize (Tup3Type ta tb tc) = ( defaultSize ta
                                  , defaultSize tb
                                  , defaultSize tc
                                  )
defaultSize (Tup4Type ta tb tc td) = ( defaultSize ta
                                     , defaultSize tb
                                     , defaultSize tc
                                     , defaultSize td
                                     )
defaultSize (Tup5Type ta tb tc td te) = ( defaultSize ta
                                        , defaultSize tb
                                        , defaultSize tc
                                        , defaultSize td
                                        , defaultSize te
                                        )
defaultSize (Tup6Type ta tb tc td te tf) = ( defaultSize ta
                                           , defaultSize tb
                                           , defaultSize tc
                                           , defaultSize td
                                           , defaultSize te
                                           , defaultSize tf
                                           )
defaultSize (Tup7Type ta tb tc td te tf tg) = ( defaultSize ta
                                              , defaultSize tb
                                              , defaultSize tc
                                              , defaultSize td
                                              , defaultSize te
                                              , defaultSize tf
                                              , defaultSize tg
                                              )
defaultSize (FunType ta tb) = (defaultSize ta, defaultSize tb)
defaultSize (MutType ta) = defaultSize ta
defaultSize (RefType ta) = defaultSize ta
defaultSize (MArrType ta) = universal :> defaultSize ta
defaultSize (ParType ta) = defaultSize ta
defaultSize (IVarType ta) = defaultSize ta
defaultSize (FValType ta) = defaultSize ta

-- | Type equality on 'Signedness'
signEq :: Signedness s1 -> Signedness s2 -> Maybe (TypeEq s1 s2)
signEq U U = Just TypeEq
signEq S S = Just TypeEq
signEq _ _ = Nothing

-- | Type equality on 'BitWidth'
widthEq :: BitWidth n1 -> BitWidth n2 -> Maybe (TypeEq n1 n2)
widthEq N8      N8      = Just TypeEq
widthEq N16     N16     = Just TypeEq
widthEq N32     N32     = Just TypeEq
widthEq N64     N64     = Just TypeEq
widthEq NNative NNative = Just TypeEq
widthEq _ _ = Nothing

-- | Type equality on 'TypeRep'
typeEq :: TypeRep a -> TypeRep b -> Maybe (TypeEq a b)
typeEq UnitType UnitType = Just TypeEq
typeEq BoolType BoolType = Just TypeEq
typeEq (IntType s1 n1) (IntType s2 n2) = do
    TypeEq <- signEq s1 s2
    TypeEq <- widthEq n1 n2
    return TypeEq
typeEq FloatType FloatType = Just TypeEq
typeEq DoubleType DoubleType = Just TypeEq
typeEq (ComplexType t1) (ComplexType t2) = do
    TypeEq <- typeEq t1 t2
    return TypeEq
typeEq (ArrayType t1) (ArrayType t2) = do
    TypeEq <- typeEq t1 t2
    return TypeEq
typeEq (TargetArrType n1 t1) (TargetArrType n2 t2) = do
    TypeEq <- widthEq n1 n2
    TypeEq <- typeEq t1 t2
    return TypeEq
typeEq (Tup2Type a1 b1) (Tup2Type a2 b2) = do
    TypeEq <- typeEq a1 a2
    TypeEq <- typeEq b1 b2
    return TypeEq
typeEq (Tup3Type a1 b1 c1) (Tup3Type a2 b2 c2) = do
    TypeEq <- typeEq a1 a2
    TypeEq <- typeEq b1 b2
    TypeEq <- typeEq c1 c2
    return TypeEq
typeEq (Tup4Type a1 b1 c1 d1) (Tup4Type a2 b2 c2 d2) = do
    TypeEq <- typeEq a1 a2
    TypeEq <- typeEq b1 b2
    TypeEq <- typeEq c1 c2
    TypeEq <- typeEq d1 d2
    return TypeEq
typeEq (Tup5Type a1 b1 c1 d1 e1) (Tup5Type a2 b2 c2 d2 e2) = do
    TypeEq <- typeEq a1 a2
    TypeEq <- typeEq b1 b2
    TypeEq <- typeEq c1 c2
    TypeEq <- typeEq d1 d2
    TypeEq <- typeEq e1 e2
    return TypeEq
typeEq (Tup6Type a1 b1 c1 d1 e1 f1) (Tup6Type a2 b2 c2 d2 e2 f2) = do
    TypeEq <- typeEq a1 a2
    TypeEq <- typeEq b1 b2
    TypeEq <- typeEq c1 c2
    TypeEq <- typeEq d1 d2
    TypeEq <- typeEq e1 e2
    TypeEq <- typeEq f1 f2
    return TypeEq
typeEq (Tup7Type a1 b1 c1 d1 e1 f1 g1) (Tup7Type a2 b2 c2 d2 e2 f2 g2) = do
    TypeEq <- typeEq a1 a2
    TypeEq <- typeEq b1 b2
    TypeEq <- typeEq c1 c2
    TypeEq <- typeEq d1 d2
    TypeEq <- typeEq e1 e2
    TypeEq <- typeEq f1 f2
    TypeEq <- typeEq g1 g2
    return TypeEq
typeEq (FunType a1 b1) (FunType a2 b2) = do
    TypeEq <- typeEq a1 a2
    TypeEq <- typeEq b1 b2
    return TypeEq
typeEq (MutType t1) (MutType t2) = do
    TypeEq <- typeEq t1 t2
    return TypeEq
typeEq (RefType t1) (RefType t2) = do
    TypeEq <- typeEq t1 t2
    return TypeEq
typeEq (MArrType a1) (MArrType a2) = do
    TypeEq <- typeEq a1 a2
    return TypeEq
typeEq (ParType t1) (ParType t2) = do
    TypeEq <- typeEq t1 t2
    return TypeEq
typeEq (IVarType t1) (IVarType t2) = do
    TypeEq <- typeEq t1 t2
    return TypeEq
typeEq (FValType t1) (FValType t2) = do
    TypeEq <- typeEq t1 t2
    return TypeEq

typeEq _ _ = Nothing

showTup :: [String] -> String
showTup as =  "(" ++ intercalate "," as ++ ")"

type family TargetType n a
type instance TargetType n ()              = ()
type instance TargetType n Bool            = Bool
type instance TargetType n Word8           = Word8
type instance TargetType n Int8            = Int8
type instance TargetType n Word16          = Word16
type instance TargetType n Int16           = Int16
type instance TargetType n Word32          = Word32
type instance TargetType n Int32           = Int32
type instance TargetType n Word64          = Word64
type instance TargetType n Int64           = Int64
type instance TargetType n WordN           = GenericInt U n
type instance TargetType n IntN            = GenericInt S n
type instance TargetType n Float           = Float
type instance TargetType n Double          = Double
type instance TargetType n (Complex a)     = Complex (TargetType n a)
type instance TargetType n [a]             = TargetArr n (TargetType n a)
type instance TargetType n (a,b)           = (TargetType n a, TargetType n b)
type instance TargetType n (a,b,c)         = (TargetType n a, TargetType n b, TargetType n c)
type instance TargetType n (a,b,c,d)       = (TargetType n a, TargetType n b, TargetType n c, TargetType n d)
type instance TargetType n (a,b,c,d,e)     = (TargetType n a, TargetType n b, TargetType n c, TargetType n d, TargetType n e)
type instance TargetType n (a,b,c,d,e,f)   = (TargetType n a, TargetType n b, TargetType n c, TargetType n d, TargetType n e, TargetType n f)
type instance TargetType n (a,b,c,d,e,f,g) = (TargetType n a, TargetType n b, TargetType n c, TargetType n d, TargetType n e, TargetType n f, TargetType n g)
type instance TargetType n (IORef a)       = IORef (TargetType n a)
type instance TargetType n (MArr a)        = MArr (TargetType n a)
type instance TargetType n (IV a)          = IV (TargetType n a)
type instance TargetType n (FVal a)        = FVal (TargetType n a)

-- | The set of supported types
class (Eq a, Show a, Typeable a, Show (Size a), Lattice (Size a)) => Type a
  where
    -- | Gives the type representation a value.
    typeRep  :: TypeRep a
    sizeOf   :: a -> Size a
    toTarget :: BitWidth n -> a -> TargetType n a

instance Type ()      where typeRep = UnitType;          sizeOf _ = AnySize;      toTarget _ = id
instance Type Bool    where typeRep = BoolType;          sizeOf _ = AnySize;      toTarget _ = id
instance Type Word8   where typeRep = IntType U N8;      sizeOf = singletonRange; toTarget _ = id
instance Type Int8    where typeRep = IntType S N8;      sizeOf = singletonRange; toTarget _ = id
instance Type Word16  where typeRep = IntType U N16;     sizeOf = singletonRange; toTarget _ = id
instance Type Int16   where typeRep = IntType S N16;     sizeOf = singletonRange; toTarget _ = id
instance Type Word32  where typeRep = IntType U N32;     sizeOf = singletonRange; toTarget _ = id
instance Type Int32   where typeRep = IntType S N32;     sizeOf = singletonRange; toTarget _ = id
instance Type Word64  where typeRep = IntType U N64;     sizeOf = singletonRange; toTarget _ = id
instance Type Int64   where typeRep = IntType S N64;     sizeOf = singletonRange; toTarget _ = id
instance Type WordN   where typeRep = IntType U NNative; sizeOf = singletonRange; toTarget = fromWordN
instance Type IntN    where typeRep = IntType S NNative; sizeOf = singletonRange; toTarget = fromIntN
instance Type Float   where typeRep = FloatType;         sizeOf _ = AnySize;      toTarget _ = id
instance Type Double  where typeRep = DoubleType;        sizeOf _ = AnySize;      toTarget _ = id

instance (Type a, RealFloat a) => Type (Complex a)
  where
    typeRep             = ComplexType typeRep
    sizeOf _            = AnySize
    toTarget _ (_ :+ _) = error "TODO" -- toTarget n r :+ toTarget n i

instance Type a => Type [a]
  where
    typeRep       = ArrayType typeRep
    sizeOf as     = singletonRange (genericLength as) :> unions (map sizeOf as)
    toTarget n as = TargetArr (genericLen n as) $ map (toTarget n) as

instance (Type a, Type b) => Type (a,b)
  where
    typeRep = Tup2Type typeRep typeRep

    sizeOf (a,b) =
        ( sizeOf a
        , sizeOf b
        )

    toTarget n (a,b) =
        ( toTarget n a
        , toTarget n b
        )

instance (Type a, Type b, Type c) => Type (a,b,c)
  where
    typeRep = Tup3Type typeRep typeRep typeRep

    sizeOf (a,b,c) =
        ( sizeOf a
        , sizeOf b
        , sizeOf c
        )

    toTarget n (a,b,c) =
        ( toTarget n a
        , toTarget n b
        , toTarget n c
        )

instance (Type a, Type b, Type c, Type d) => Type (a,b,c,d)
  where
    typeRep = Tup4Type typeRep typeRep typeRep typeRep

    sizeOf (a,b,c,d) =
        ( sizeOf a
        , sizeOf b
        , sizeOf c
        , sizeOf d
        )

    toTarget n (a,b,c,d) =
        ( toTarget n a
        , toTarget n b
        , toTarget n c
        , toTarget n d
        )

instance (Type a, Type b, Type c, Type d, Type e) => Type (a,b,c,d,e)
  where
    typeRep = Tup5Type typeRep typeRep typeRep typeRep typeRep

    sizeOf (a,b,c,d,e) =
        ( sizeOf a
        , sizeOf b
        , sizeOf c
        , sizeOf d
        , sizeOf e
        )

    toTarget n (a,b,c,d,e) =
        ( toTarget n a
        , toTarget n b
        , toTarget n c
        , toTarget n d
        , toTarget n e
        )

instance (Type a, Type b, Type c, Type d, Type e, Type f) => Type (a,b,c,d,e,f)
  where
    typeRep = Tup6Type typeRep typeRep typeRep typeRep typeRep typeRep

    sizeOf (a,b,c,d,e,f) =
        ( sizeOf a
        , sizeOf b
        , sizeOf c
        , sizeOf d
        , sizeOf e
        , sizeOf f
        )

    toTarget n (a,b,c,d,e,f) =
        ( toTarget n a
        , toTarget n b
        , toTarget n c
        , toTarget n d
        , toTarget n e
        , toTarget n f
        )

instance (Type a, Type b, Type c, Type d, Type e, Type f, Type g) => Type (a,b,c,d,e,f,g)
  where
    typeRep = Tup7Type typeRep typeRep typeRep typeRep typeRep typeRep typeRep

    sizeOf (a,b,c,d,e,f,g) =
        ( sizeOf a
        , sizeOf b
        , sizeOf c
        , sizeOf d
        , sizeOf e
        , sizeOf f
        , sizeOf g
        )

    toTarget n (a,b,c,d,e,f,g) =
        ( toTarget n a
        , toTarget n b
        , toTarget n c
        , toTarget n d
        , toTarget n e
        , toTarget n f
        , toTarget n g
        )

instance (Type a, Show (IORef a)) => Type (IORef a)
  where
    typeRep = RefType typeRep

    sizeOf _ = universal

    toTarget = error "toTarget: IORef"  -- TODO Requires IO

instance Type a => Type (MArr a)
  where
    typeRep = MArrType typeRep

    sizeOf _ = universal

    toTarget = error "toTarget: MArr"  -- TODO Requires IO

instance Type a => Type (IV a)
  where
    typeRep = IVarType typeRep

    sizeOf _ = universal

    toTarget = error "toTarget: IVar" -- TODO Requires IO

instance Type a => Type (FVal a)
  where
    typeRep = FValType typeRep

    sizeOf _ = universal

    toTarget = error "toTarget: FVal" -- TODO

typeRepByProxy :: Type a => Proxy a -> TypeRep a
typeRepByProxy _ = typeRep

--------------------------------------------------------------------------------
-- * Sized types
--------------------------------------------------------------------------------

data AnySize = AnySize
    deriving (Eq, Show, Ord)

anySizeFun :: AnySize -> AnySize
anySizeFun AnySize = AnySize

anySizeFun2 :: AnySize -> AnySize -> AnySize
anySizeFun2 AnySize AnySize = AnySize

instance Num AnySize
  where
    fromInteger _ = AnySize
    abs           = anySizeFun
    signum        = anySizeFun
    (+)           = anySizeFun2
    (-)           = anySizeFun2
    (*)           = anySizeFun2

instance Lattice AnySize
  where
    bot  = AnySize
    top  = AnySize
    (\/) = anySizeFun2
    (/\) = anySizeFun2

type family Size a

type instance Size ()              = AnySize
type instance Size Bool            = AnySize
type instance Size Word8           = Range Word8
type instance Size Int8            = Range Int8
type instance Size Word16          = Range Word16
type instance Size Int16           = Range Int16
type instance Size Word32          = Range Word32
type instance Size Int32           = Range Int32
type instance Size Word64          = Range Word64
type instance Size Int64           = Range Int64
type instance Size WordN           = Range WordN
type instance Size IntN            = Range IntN
type instance Size Float           = AnySize
type instance Size Double          = AnySize
type instance Size (Complex a)     = AnySize
type instance Size [a]             = Range Length :> Size a
type instance Size (TargetArr n a) = Range (GenericInt U n) :> Size a
type instance Size (a,b)           = (Size a, Size b)
type instance Size (a,b,c)         = (Size a, Size b, Size c)
type instance Size (a,b,c,d)       = (Size a, Size b, Size c, Size d)
type instance Size (a,b,c,d,e)     = (Size a, Size b, Size c, Size d, Size e)
type instance Size (a,b,c,d,e,f)   = (Size a, Size b, Size c, Size d, Size e, Size f)
type instance Size (a,b,c,d,e,f,g) = (Size a, Size b, Size c, Size d, Size e, Size f, Size g)
type instance Size (a -> b)        = (Size a, Size b)
type instance Size (Mut a)         = Size a
type instance Size (IORef a)       = Size a
type instance Size (MArr a)        = Range Length :> Size a
type instance Size (Par a)         = Size a
type instance Size (IV a)          = Size a
type instance Size (FVal a)        = Size a

-- | A generalization of 'Range' that serves two purposes: (1) Adding an extra
-- 'Universal' constructor to support unbounded types ('Range' can only
-- represent bounded ranges), and (2) pack a 'BoundedInt' constraint with the
-- 'RangeSet' constructor. This is what allows 'sizeToRange' to be defined as a
-- total function with 'Type' as the only constraint.
data RangeSet a
  where
    RangeSet  :: BoundedInt a => Range a -> RangeSet a
    Universal :: RangeSet a

-- | Cast a 'Size' to a 'RangeSet'
sizeToRange :: forall a . Type a => Size a -> RangeSet a
sizeToRange sz = case typeRep :: TypeRep a of
    IntType _ _ -> RangeSet sz
    _           -> Universal

tIntN :: Patch IntN IntN
tIntN = id

tWordN :: Patch WordN WordN
tWordN = id

tIndex :: Patch Index Index
tIndex  = id

tLength :: Patch Length Length
tLength = id

tArr :: Patch a a -> Patch [a] [a]
tArr _ = id

max_stars_repo_name: kffaxen/feldspar-language

max_issues_repo_name: kffaxen/feldspar-language

max_forks_repo_name: kffaxen/feldspar-language


text:
{-# LANGUAGE RecordWildCards #-}

module Numeric.Solver where

import           Control.Exception
import           Data.Function
import           Data.List
import qualified Data.Vector.Storable       as Storable
import           Numeric.GSL.ODE
import           Numeric.LinearAlgebra.Data
import           Numeric.SSystem
import           Text.Parse.ODEBench

type Objective = SSystem Double -> Double

logLikelihood :: Matrix Double -- ^ Observed values
              -> Matrix Double -- ^ Standard deviations
              -> Matrix Double -- ^ Simulated values
              -> Double
logLikelihood o d s =
  (-0.5) * Storable.sum
  (Storable.zipWith3 ll
   (flatten o)
   (flatten d)
   (flatten s))
  where ll obs stddev sim = ((sim - obs) / stddev) ** 2

data Experiment = Experiment
  { times           :: Storable.Vector Double
  , initalVals      :: [Double]
  , observedVals    :: Matrix Double
  , observedStddevs :: Matrix Double }

experiments :: Problem -> [Experiment]
experiments = map toExpr . groupBy ((==) `on` experimentNum) . samples where
  toExpr = Experiment <$> Storable.fromList . map time
                      <*> Storable.toList . values . head
                      <*> fromRows . map values
                      <*> fromRows . map stdDevs

solve :: SSystem Double -> [Double] -> Storable.Vector Double -> Matrix Double
solve = odeSolve . const . toEqns'

likelihood :: [Experiment] -> Objective
likelihood exprs syst = foldl' (+) 0 (map f exprs) where
  f Experiment {..} =
    logLikelihood
      observedVals
      observedStddevs
      (solve syst initalVals times)

type SafeObjective = SSystem Double -> IO (Either SomeException Double)

safeLikelihood :: [Experiment] -> SafeObjective
safeLikelihood expr syst = (try . evaluate) (likelihood expr syst)

-- | >>> testSolve exampleSystem [0.7,0.12,0.14,0.16,0.18] 6 (0,0.25)
-- (6><5)
--  [                0.7,               0.12,               0.14,               0.16,                0.18
--  , 1.1167718442535226, 0.5288004914151249, 0.9744853978475282, 1.0888238837276851, 0.40157122204470663
--  , 0.9667747407434063,  0.827179819885076, 0.9938518679839954, 1.0925434854072784,  0.8240710601831301
--  , 0.8394130067100113, 0.8627637078683361, 0.9980826732970266, 0.9920051282448126,  0.9522071825140638
--  , 0.7765858895457368, 0.8274811279700081,  0.999411143907833,  0.943541476235036,  0.9572057285536645
--  , 0.7499481576778694, 0.7905920430606405, 0.9998291350909083, 0.9281457408781921,  0.9424146353794791 ]
testSolve :: SSystem Double -> [Double] -> Int -> (Double,Double) -> Matrix Double
testSolve s t n i = solve s t (linspace n i)

max_stars_repo_name: oisdk/SSystemOpt

max_issues_repo_name: oisdk/SSystemOpt

max_forks_repo_name: oisdk/SSystemOpt


text:
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedLabels #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE UndecidableInstances #-}

module Day12 where

import Control.Monad.State (MonadState, State, execState)
import Data.ByteString qualified as B
import Data.Complex (Complex((:+)), cis)
import Data.Functor (($>))
import Data.Text qualified as T
import Data.Text.Encoding qualified as TE
import Data.Void (Void)
import Flow ((.>))
import Optics
    ( Optic, A_Setter, Is, view, noPrefixFieldLabels, makeFieldLabelsWith )
import Optics.State (use)
import Optics.State.Operators ((%=))
import Text.Megaparsec qualified as P
import Text.Megaparsec.Char qualified as PC
import Text.Megaparsec.Char.Lexer qualified as PCL
-- import Text.Pretty.Simple (pPrint)

type Parser = P.Parsec Void T.Text

data NavAction
    = MoveNorth
    | MoveSouth
    | MoveEast
    | MoveWest
    | TurnLeft
    | TurnRight
    | MoveForward
    deriving (Bounded, Enum, Eq, Ord, Read, Show)

data NavInstr = MkNavInstr
    { navAction :: !NavAction
    , navValue :: !Int
    } deriving (Eq, Ord, Read, Show)

data ShipState = MkShipState
    { shipPhase :: !Double
    , shipPosition :: !(Complex Double)
    , wayPosition :: !(Complex Double)
    } deriving (Eq, Read, Show)

makeFieldLabelsWith noPrefixFieldLabels ''ShipState

main :: IO ()
main = do
    file <- readFileUtf8 "day-12/input.txt"
    case P.parse (parseNavInstrs <* P.eof) "day 12 input" file of
        Left err -> putStrLn $ P.errorBundlePretty err
        Right navInstrs -> do
            print $ part1 navInstrs
            print $ part2 navInstrs

part2 :: [NavInstr] -> Double
part2 = pilotShip2 .> view #shipPosition .> manhattanNorm

part1 :: [NavInstr] -> Double
part1 = pilotShip1 .> view #shipPosition .> manhattanNorm

pilotShip2 :: [NavInstr] -> ShipState
pilotShip2 instrs = execState (applyInstrs2 instrs) $
    MkShipState
        { shipPhase = 0
        , shipPosition = pure 0
        , wayPosition = 10 :+ 1
        }

pilotShip1 :: [NavInstr] -> ShipState
pilotShip1 instrs = execState (applyInstrs1 instrs) $
    MkShipState
        { shipPhase = 0
        , shipPosition = pure 0
        , wayPosition = pure 0
        }

applyInstrs1 :: [NavInstr] -> State ShipState ()
applyInstrs1 = foldr (\instr st -> applyInstr1 instr *> st) (pure ())

applyInstr1 :: NavInstr -> State ShipState ()
applyInstr1 (MkNavInstr action intVal) = case action of
    MoveNorth -> #shipPosition += (0 :+ val)
    MoveSouth -> #shipPosition -= (0 :+ val)
    MoveEast -> #shipPosition += (val :+ 0)
    MoveWest -> #shipPosition -= (val :+ 0)
    TurnLeft -> #shipPhase %= (+ val) .> (`modulo` 360)
    TurnRight -> #shipPhase %= subtract val .> (`modulo` 360)
    MoveForward -> do
        phase <- use #shipPhase
        #shipPosition += fmap (* val) (cisDegrees phase)
  where
    val = fromIntegral intVal

applyInstrs2 :: [NavInstr] -> State ShipState ()
applyInstrs2 = foldr (\instr st -> applyInstr2 instr *> st) (pure ())

applyInstr2 :: NavInstr -> State ShipState ()
applyInstr2 (MkNavInstr action intVal) = case action of
    MoveNorth -> #wayPosition += (0 :+ val)
    MoveSouth -> #wayPosition -= (0 :+ val)
    MoveEast -> #wayPosition += (val :+ 0)
    MoveWest -> #wayPosition -= (val :+ 0)
    TurnLeft -> #wayPosition *= cisDegrees val
    TurnRight -> #wayPosition *= cisDegrees (-val)
    MoveForward -> do
        wayPos <- use #wayPosition
        #shipPosition += fmap (* val) wayPos
  where
    val = fromIntegral intVal

cisDegrees :: Double -> Complex Double
cisDegrees degRaw
    | deg == 0 = 1 :+ 0
    | deg == 90 = 0 :+ 1
    | deg == 180 = (-1) :+ 0
    | deg == 270 = 0 :+ (-1)
    | otherwise = cis rad
  where
    deg = degRaw `modulo` 360
    rad = deg * pi / 180

modulo :: Double -> Double -> Double
modulo x m = x - floorDbl (x / m) * m
  where
    floorDbl = floor .> fromIntegral @Int

manhattanDist :: RealFloat a => Complex a -> Complex a -> a
manhattanDist u v = sum $ fmap abs $ v - u

manhattanNorm :: RealFloat a => Complex a -> a
manhattanNorm = manhattanDist (pure 0)

parseNavInstrs :: Parser [NavInstr]
parseNavInstrs = parseNavInstr `P.sepEndBy` PC.space

parseNavInstr :: Parser NavInstr
parseNavInstr = do
    action <- parseNavAction
    val <- PC.space *> PCL.decimal
    pure $ MkNavInstr action val

parseNavAction :: Parser NavAction
parseNavAction = P.choice
    [ PC.char 'N' $> MoveNorth
    , PC.char 'S' $> MoveSouth
    , PC.char 'E' $> MoveEast
    , PC.char 'W' $> MoveWest
    , PC.char 'L' $> TurnLeft
    , PC.char 'R' $> TurnRight
    , PC.char 'F' $> MoveForward
    ]

readFileUtf8 :: FilePath -> IO T.Text
readFileUtf8 path = TE.decodeUtf8 <$> B.readFile path

(+=)
    :: (Is k A_Setter, MonadState s m, Num a)
    => Optic k is s s a a
    -> a
    -> m ()
optic += x = optic %= (+ x)

(-=)
    :: (Is k A_Setter, MonadState s m, Num a)
    => Optic k is s s a a
    -> a
    -> m ()
optic -= x = optic %= subtract x

(*=)
    :: (Is k A_Setter, MonadState s m, Num a)
    => Optic k is s s a a
    -> a
    -> m ()
optic *= x = optic %= (* x)


text:
--
-- PoolLayer : pooling layer
--

module CNN.PoolLayer (
  poolMax
, depoolMax
, reversePooling
) where

import Numeric.LinearAlgebra

import CNN.Algebra
import CNN.Image
import CNN.LayerType

type Pix = (Double, Double)

{- |
poolMax

  In : pooling size (x = y)
       image

  OUT: updated image
       position of max values

>>> let im = [fromLists [[1.0,2.0,3.0,4.0],[8.0,7.0,6.0,5.0],[1.0,3.0,5.0,7.0],[2.0,4.0,6.0,8.0]], fromLists [[1.0,2.0,3.0,4.0],[2.0,2.0,3.0,4.0],[3.0,3.0,3.0,4.0],[4.0,4.0,4.0,4.0]]]
>>> poolMax 2 im
[[(2><2)
 [ 8.0, 6.0
 , 4.0, 8.0 ],(2><2)
 [ 2.0, 4.0
 , 4.0, 4.0 ]],[(2><2)
 [ 2.0, 2.0
 , 3.0, 3.0 ],(2><2)
 [ 1.0, 1.0
 , 2.0, 1.0 ]]]

-}

poolMax :: Int -> Image -> [Image]
poolMax s im = [os, is]
  where
    pl = head im
    x  = cols pl `div` s
    y  = rows pl `div` s
    ps = [(i*s, j*s) | i <- [0..(x-1)], j <- [0..(y-1)]]
    (os, is) = unzip $ map (toPlain x y . maxPix s ps) im

toPlain :: Int -> Int -> [Pix] -> (Plain, Plain)
toPlain x y pls = ((x><y) op, (x><y) ip)
  where
    (op, ip) = unzip pls

maxPix :: Int -> [(Int, Int)] -> Plain -> [Pix]
maxPix s ps is = map (max' . toPix) ps
  where
    toPix :: (Int, Int) -> [Pix]
    toPix p = zip (concat $ toLists $ subMatrix p (s, s) is) [0.0..]

max' :: [Pix] -> Pix
max' [] = error "empty list!"
max' [x] = x
max' (x:xs) = maximum' x (max' xs)

maximum' :: Pix -> Pix -> Pix
maximum' a@(v1, _) b@(v2, _) = if v1 < v2 then b else a

-- back prop

{- |
depoolMax

>>> let im = [fromLists [[0.0,2.0],[3.0,1.0]]]
>>> let dl = [fromLists [[0.1,0.2],[0.3,0.4]]]
>>> depoolMax 2 im dl
([(4><4)
 [ 0.1, 0.0, 0.0, 0.0
 , 0.0, 0.0, 0.2, 0.0
 , 0.0, 0.0, 0.0, 0.4
 , 0.0, 0.3, 0.0, 0.0 ]],Nothing)

-}

depoolMax :: Int -> Image -> Delta -> (Delta, Maybe Layer)
depoolMax s im d = (zipWith (concatWith (expand s 0.0)) im d, Nothing)
  where
    concatWith :: (Double -> Double -> Matrix R) ->  Matrix R -> Matrix R
               -> Matrix R
    concatWith f i d = fromBlocks $ zipWith (zipWith f) is ds
      where
        is = toLists i
        ds = toLists d

{- |
expand

  IN : size of pooling
       filling value
       positions
       delta values

>>> expand 2 0.0 0 3
(2><2)
 [ 3.0, 0.0
 , 0.0, 0.0 ]
>>> expand 2 0.0 2 4
(2><2)
 [ 0.0, 0.0
 , 4.0, 0.0 ]
>>> expand 3 0.0 5 3
(3><3)
 [ 0.0, 0.0, 0.0
 , 0.0, 0.0, 3.0
 , 0.0, 0.0, 0.0 ]
>>> expand 3 0.0 8 3
(3><3)
 [ 0.0, 0.0, 0.0
 , 0.0, 0.0, 0.0
 , 0.0, 0.0, 3.0 ]

-}

expand :: Int -> Double -> Double -> Double -> Matrix R
expand s r p d = (s><s) $ (replicate p1 r ++ [d] ++ replicate p2 r)
  where
    p1 = truncate p
    p2 = max 0 (s * s - p1 - 1)

-- reverse

reversePooling :: Int -> Layer
reversePooling = MaxPoolLayer

max_stars_repo_name: eijian/deeplearning

max_issues_repo_name: eijian/deeplearning

max_forks_repo_name: eijian/deeplearning


text:
{-# LANGUAGE FlexibleContexts                         #-}
{-# LANGUAGE GADTs                                    #-}
{-# LANGUAGE TypeInType                               #-}
{-# LANGUAGE TypeOperators                            #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}

module AOC.Common.Numeric (
    fft
  , ifft
  , convolve
  , rconvolve
  , zconvolve
  , FFT.FFTWReal
  ) where

import           Data.Complex
import           GHC.TypeNats
import qualified Data.Array.CArray         as CA
import qualified Data.Array.IArray         as IA
import qualified Data.Ix                   as Ix
import qualified Data.Vector.Generic       as VG
import qualified Data.Vector.Generic.Sized as SVG
import qualified Foreign.Storable          as FS
import qualified Math.FFT                  as FFT
import qualified Math.FFT.Base             as FFT

fft :: (FFT.FFTWReal a, VG.Vector v (Complex a))
    => SVG.Vector v n (Complex a)
    -> SVG.Vector v n (Complex a)
fft = SVG.withVectorUnsafe $
        fromCA
      . FFT.dft
      . toCA

ifft
    :: (FFT.FFTWReal a, VG.Vector v (Complex a))
    => SVG.Vector v n (Complex a)
    -> SVG.Vector v n (Complex a)
ifft = SVG.withVectorUnsafe $
        fromCA
      . FFT.idft
      . toCA

fromCA
    :: (FS.Storable a, VG.Vector v (Complex a))
    => CA.CArray Int (Complex a)
    -> v (Complex a)
fromCA v = VG.generate (Ix.rangeSize (IA.bounds v)) (v IA.!)

toCA
    :: (FS.Storable a, VG.Vector v (Complex a))
    => v (Complex a)
    -> CA.CArray Int (Complex a)
toCA v = IA.listArray (0, VG.length v - 1) (VG.toList v)

-- | FFT-based convolution
convolve
    :: ( VG.Vector v (Complex a)
       , KnownNat n, 1 <= n
       , KnownNat m, 1 <= m
       , FFT.FFTWReal a
       )
    => SVG.Vector v n (Complex a)
    -> SVG.Vector v m (Complex a)
    -> SVG.Vector v (n + m - 1) (Complex a)
convolve x y = ifft $ fft x' * fft y'
  where
    x' = x SVG.++ 0
    y' = y SVG.++ 0

-- | FFT-based real-valued convolution
rconvolve
    :: ( VG.Vector v (Complex a)
       , VG.Vector v a
       , KnownNat n, 1 <= n
       , KnownNat m, 1 <= m
       , FFT.FFTWReal a
       )
    => SVG.Vector v n a
    -> SVG.Vector v m a
    -> SVG.Vector v (n + m - 1) a
rconvolve x y = SVG.map realPart $ convolve (SVG.map (:+ 0) x) (SVG.map (:+ 0) y)

-- | FFT-based integral convolution
zconvolve
    :: ( VG.Vector v (Complex Double)
       , VG.Vector v Double
       , VG.Vector v a
       , KnownNat n, 1 <= n
       , KnownNat m, 1 <= m
       , Integral a
       )
    => SVG.Vector v n a
    -> SVG.Vector v m a
    -> SVG.Vector v (n + m - 1) a
zconvolve x y = SVG.map (round @Double) $
    rconvolve (SVG.map fromIntegral x) (SVG.map fromIntegral y)