text:
{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE KindSignatures      #-}
{-# LANGUAGE RecordWildCards     #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications    #-}
{-# LANGUAGE TypeInType          #-}
{-# LANGUAGE TypeOperators       #-}

module Eigen.SparseMatrix
  ( -- * Types
    SparseMatrix(..)
  , SparseMatrixXf
  , SparseMatrixXd
  , SparseMatrixXcf
  , SparseMatrixXcd

    -- * Matrix internal data
  , elems
  , values
  , innerIndices
  , outerStarts
  , innerNNZs

    -- * Accessors
  , cols
  , rows
  , coeff
  , (!)
  , getRow
  , getCol
  , getRows
  , getCols
  , unsafeGetCol
  , unsafeGetRow
  , unsafeFromRows
  , unsafeFromCols

    -- * Matrix conversions
  , fromList
  , toList
  , fromVector
  , toVector
  , fromVectorSized
  , fromListSized
  , fromDenseList
  , toDenseList
  , fromMatrix
  , toMatrix

    -- * Matrix properties
  , norm
  , squaredNorm
  , blueNorm
  , block
  , unsafeBlock
  , nonZeros
  , innerSize
  , outerSize

    -- * Basic matrix algebra
  , add
  , sub
  , mul

    -- * Matrix tranformations
  , pruned
  , scale
  , transpose
  , adjoint
  , map
  , imap

    -- * Matrix representation
  , compress
  , uncompress
  , compressed

    -- * Matrix serialisation
  , encode
  , decode

    -- * Mutable matrices
  , thaw
  , freeze
  , unsafeThaw
  , unsafeFreeze
  ) where

import Control.Monad (when, guard)
import Control.Monad.Primitive (PrimMonad(..), unsafePrimToPrim)
import qualified Prelude
import Prelude hiding (read, map)
import Data.Binary (Binary(..))
import qualified Data.Binary as Binary
import qualified Data.ByteString.Lazy as BSL
import Data.Complex (Complex)
import Data.Kind (Type)
import qualified Data.Vector.Storable as VS
import qualified Data.Vector.Storable.Mutable as VSM
import Foreign.C.String (CString)
import Foreign.C.Types (CInt)
import Foreign.ForeignPtr (ForeignPtr, withForeignPtr, touchForeignPtr)
import Foreign.Marshal.Alloc (alloca)
import Foreign.Ptr (Ptr)
import Foreign.Storable (Storable(..))
import qualified Foreign.Concurrent as FC
import GHC.TypeLits (Nat, KnownNat, type (<=))
import Eigen.Internal
  ( Elem
  , Cast(..)
  , CSparseMatrix
  , CSparseMatrixPtr
  , natToInt
  , Row(..)
  , Col(..)
  , CTriplet(..)
  )
import qualified Data.List as List
import qualified Eigen.Internal as Internal
import qualified Eigen.Matrix as M
import qualified Eigen.Matrix.Mutable as MM
import qualified Eigen.SparseMatrix.Mutable as SMM

{-| A versatible sparse matrix representation.
SparseMatrix is the main sparse matrix representation of Eigen's sparse module.
It offers high performance and low memory usage.
It implements a more versatile variant of the widely-used Compressed Column (or Row) Storage scheme.
It consists of four compact arrays:
* `values`: stores the coefficient values of the non-zeros.
* `innerIndices`: stores the row (resp. column) indices of the non-zeros.
* `outerStarts`: stores for each column (resp. row) the index of the first non-zero in the previous two arrays.
* `innerNNZs`: stores the number of non-zeros of each column (resp. row). The word inner refers to an inner vector that is a column for a column-major matrix, or a row for a row-major matrix. The word outer refers to the other direction.
This storage scheme is better explained on an example. The following matrix
@
0   3   0   0   0
22  0   0   0   17
7   5   0   1   0
0   0   0   0   0
0   0   14  0   8
@
and one of its possible sparse, __column major__ representation:
@
values:         22  7   _   3   5   14  _   _   1   _   17  8
innerIndices:   1   2   _   0   2   4   _   _   2   _   1   4
outerStarts:    0   3   5   8   10  12
innerNNZs:      2   2   1   1   2
@
Currently the elements of a given inner vector are guaranteed to be always sorted by increasing inner indices.
The "\_" indicates available free space to quickly insert new elements. Assuming no reallocation is needed,
the insertion of a random element is therefore in @O(nnz_j)@ where @nnz_j@ is the number of nonzeros of the
respective inner vector. On the other hand, inserting elements with increasing inner indices in a given inner
vector is much more efficient since this only requires to increase the respective `innerNNZs` entry that is a @O(1)@ operation.
The case where no empty space is available is a special case, and is refered as the compressed mode.
It corresponds to the widely used Compressed Column (or Row) Storage schemes (CCS or CRS).
Any `SparseMatrix` can be turned to this form by calling the `compress` function.
In this case, one can remark that the `innerNNZs` array is redundant with `outerStarts` because we the equality:
@InnerNNZs[j] = OuterStarts[j+1]-OuterStarts[j]@. Therefore, in practice a call to `compress` frees this buffer.
The results of Eigen's operations always produces compressed sparse matrices.
On the other hand, the insertion of a new element into a `SparseMatrix` converts this later to the uncompressed mode.
For more infomration please see Eigen <http://eigen.tuxfamily.org/dox/classEigen_1_1SparseMatrix.html documentation page>.
-}
newtype SparseMatrix :: Nat -> Nat -> Type -> Type where
  SparseMatrix :: ForeignPtr (CSparseMatrix a) -> SparseMatrix n m a

instance forall n m a. (Elem a, Show a, KnownNat n, KnownNat m) => Show (SparseMatrix n m a) where
  show m = concat
    [ "SparseMatrix"
    , show (rows_ m)
    , "x"
    , show (cols_ m)
    , "\n"
    , List.intercalate "\n" $ Prelude.map (List.intercalate "\t" . Prelude.map show) $ toDenseList m
    , "\n"
    ]

instance forall n m a. (Elem a, KnownNat n, KnownNat m) => Binary (SparseMatrix n m a) where
  put mat = do
    put $ Internal.magicCode (undefined :: C a)
    put $ natToInt @n
    put $ natToInt @m
    put $ toVector mat

  get = do
    get >>= (`when` fail "wrong matrix type") . (/= Internal.magicCode (undefined :: C a))
    fromVector <$> get

-- | Encode the sparse matrix as a lazy bytestring
encode :: (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> BSL.ByteString
encode = Binary.encode

-- | Decode the sparse matrix from a lazy bytestring
decode :: (Elem a, KnownNat n, KnownNat m) => BSL.ByteString -> SparseMatrix n m a
decode = Binary.decode

-- | Alias for single precision sparse matrix
type SparseMatrixXf  n m = SparseMatrix n m Float
-- | Alias for double precision sparse matrix
type SparseMatrixXd  n m = SparseMatrix n m Double
-- | Alias for single previsiom sparse matrix of complex numbers
type SparseMatrixXcf n m = SparseMatrix n m (Complex Float)
-- | Alias for double prevision sparse matrix of complex numbers
type SparseMatrixXcd n m = SparseMatrix n m (Complex Double)

-- | Get the coefficient values of the non-zeros.
values :: Elem a => SparseMatrix n m a -> VS.Vector a
values = VS.map fromC . _getvec Internal.sparse_values

-- | Get the row indices of the non-zeros.
innerIndices :: Elem a => SparseMatrix n m a -> VS.Vector Int
innerIndices = VS.map fromC . _getvec Internal.sparse_innerIndices

-- | Gets for each column the index of the first non-zero in the previous two arrays.
outerStarts :: Elem a => SparseMatrix n m a -> VS.Vector Int
outerStarts = VS.map fromC . _getvec Internal.sparse_outerStarts

-- | Gets the number of non-zeros of each column.
-- The word inner refers to an inner vector that is a column for a column-major matrix, or a row for a row-major matrix.
-- The word outer refers to the other direction
innerNNZs :: Elem a => SparseMatrix n m a -> Maybe (VS.Vector Int)
innerNNZs m
  | compressed m = Nothing
  | otherwise    = Just $ VS.map fromC $ _getvec Internal.sparse_innerNNZs m

-- | Number of rows in the sparse matrix
rows :: forall n m a. (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> Row n
rows _ = Row @n

-- | Number of colums in the sparse matrix
cols :: forall n m a. (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> Col m
cols _ = Col @m

-- | Number of rows in the sparse matrix
rows_ :: forall n m a. (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> Int
rows_ _ = natToInt @n

-- | Number of colums in the sparse matrix
cols_ :: forall n m a. (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> Int
cols_ _ = natToInt @m

-- | Sparse matrix coefficient at the given row and column
coeff :: forall n m r c a. (Elem a, KnownNat n, KnownNat m, KnownNat r, KnownNat c, r <= n, c <= m)
  => Row r -> Col c -> SparseMatrix n m a -> a
coeff _ _ (SparseMatrix fp) =
  let !c_row = toC $! natToInt @r
      !c_col = toC $! natToInt @c
  in Internal.performIO $ withForeignPtr fp $ \p -> alloca $ \pq -> do
       Internal.call $ Internal.sparse_coeff p c_row c_col pq
       fromC <$> peek pq

-- | Matrix coefficient at the given row and column
(!) :: forall n m r c a. (Elem a, KnownNat n, KnownNat m, KnownNat r, KnownNat c, r <= n, c <= m)
  => SparseMatrix n m a -> (Row r, Col c) -> a
(!) m (row,col) = coeff row col m

{-| For vectors, the l2 norm, and for matrices the Frobenius norm.
    In both cases, it consists in the square root of the sum of the square of all the matrix entries.
    For vectors, this is also equals to the square root of the dot product of this with itself.
-}
norm :: Elem a => SparseMatrix n m a -> a
norm = _unop Internal.sparse_norm (pure . fromC)

-- | For vectors, the squared l2 norm, and for matrices the Frobenius norm. In both cases, it consists in the sum of the square of all the matrix entries. For vectors, this is also equals to the dot product of this with itself.
squaredNorm :: Elem a => SparseMatrix n m a -> a
squaredNorm = _unop Internal.sparse_squaredNorm (pure . fromC)

-- | The l2 norm of the matrix using the Blue's algorithm. A Portable Fortran Program to Find the Euclidean Norm of a Vector, ACM TOMS, Vol 4, Issue 1, 1978.
blueNorm :: Elem a => SparseMatrix n m a -> a
blueNorm = _unop Internal.sparse_blueNorm (pure . fromC)

unsafeBlock ::
     (Elem a)
  => Int
  -> Int
  -> Int
  -> Int
  -> SparseMatrix n m a
  -> SparseMatrix n1 m1 a
unsafeBlock !sr !sc !br !bc =
  let !c_startRow = toC $! sr
      !c_startCol = toC $! sc
      !c_rows     = toC $! br
      !c_cols     = toC $! bc
  in _unop (\p pq -> Internal.sparse_block p c_startRow c_startCol c_rows c_cols pq) _mk

-- | Extract a rectangular block from the sparse matrix, given a startRow, startCol, blockRows, blockCols
block :: forall sr sc br bc n m a.
     (Elem a, KnownNat sr, KnownNat sc, KnownNat br, KnownNat bc, KnownNat n, KnownNat m)
  => (sr <= n, sc <= m, br <= n, bc <= m)
  => Row sr
  -> Col sc
  -> Row br
  -> Col bc
  -> SparseMatrix n m a
  -> SparseMatrix br bc a
block _ _ _ _ =
  let !c_startRow = toC $! natToInt @sr
      !c_startCol = toC $! natToInt @sc
      !c_rows     = toC $! natToInt @br
      !c_cols     = toC $! natToInt @bc
  in _unop (\p pq -> Internal.sparse_block p c_startRow c_startCol c_rows c_cols pq) _mk

-- | Number of non-zeros elements in the sparse matrix
nonZeros :: Elem a => SparseMatrix n m a -> Int
nonZeros = _unop Internal.sparse_nonZeros (pure . fromC)

-- | Number of elements in the sparse matrix, including zeros
elems :: forall n m a. (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> Int
elems _ = (natToInt @n) * (natToInt @m)

-- | The matrix in the compressed format
compress :: Elem a => SparseMatrix n m a -> SparseMatrix n m a
compress = _unop Internal.sparse_makeCompressed _mk

-- | The matrix in the uncompressed format
uncompress :: Elem a => SparseMatrix n m a -> SparseMatrix n m a
uncompress = _unop Internal.sparse_uncompress _mk

-- | Is the matrix compressed?
compressed :: Elem a => SparseMatrix n m a -> Bool
compressed = _unop Internal.sparse_isCompressed (pure . (/=0))

-- | Minor dimension with respect to the storage order
innerSize :: Elem a => SparseMatrix n m a -> Int
innerSize = _unop Internal.sparse_innerSize (pure . fromC)

-- | Major dimension with respect to the storage order
outerSize :: Elem a => SparseMatrix n m a -> Int
outerSize = _unop Internal.sparse_outerSize (pure . fromC)

-- | Suppresses all nonzeros which are much smaller than the reference under the tolerance @epsilon@
pruned :: Elem a => a -> SparseMatrix n m a -> SparseMatrix n m a
pruned r = _unop (\p pq -> alloca $ \pr -> poke pr (toC r) >> Internal.sparse_prunedRef p pr pq) _mk

-- | Multiply matrix on a given scalar
scale :: Elem a => a -> SparseMatrix n m a -> SparseMatrix n m a
scale x = _unop (\p pq -> alloca $ \px -> poke px (toC x) >> Internal.sparse_scale p px pq) _mk

-- | Transpose of the sparse matrix
transpose :: Elem a => SparseMatrix n m a -> SparseMatrix m n a
transpose = _unop Internal.sparse_transpose _mk

-- | Adjoint of the sparse matrix
adjoint :: Elem a => SparseMatrix n m a -> SparseMatrix m n a
adjoint = _unop Internal.sparse_adjoint _mk

-- | Add two sparse matrices by adding the corresponding entries together.
add :: Elem a => SparseMatrix n m a -> SparseMatrix n m a -> SparseMatrix n m a
add = _binop Internal.sparse_add _mk

-- | Subtract two sparse matrices by subtracting the corresponding entries together.
sub :: Elem a => SparseMatrix n m a -> SparseMatrix n m a -> SparseMatrix n m a
sub = _binop Internal.sparse_sub _mk

-- | Matrix multiplication.
mul :: Elem a => SparseMatrix p q a -> SparseMatrix q r a -> SparseMatrix p r a
mul = _binop Internal.sparse_mul _mk

-- | Map a function over the 'SparseMatrix'.
map :: (Elem a, Elem b, KnownNat n, KnownNat m) => (a -> b) -> SparseMatrix n m a -> SparseMatrix n m b
map f m = fromVector . VS.map g . toVector $ m where
  g (CTriplet r c v) = CTriplet r c $ (toC . f . fromC) v

-- | Map an indexed function over the 'SparseMatrix'.
imap :: (Elem a, Elem b, KnownNat n, KnownNat m) => (Int -> Int -> a -> b) -> SparseMatrix n m a -> SparseMatrix n m b
imap f m = fromVector . VS.map g . toVector $ m where
  g (CTriplet r c v) =
    let !_r = fromC r
        !_c = fromC c
        !_v = fromC v
    in CTriplet r c $ toC $ f _r _c _v

-- | Construct asparse matrix of the given size from the storable vector of triplets (row, col, val)
fromVector :: forall n m a. (Elem a, KnownNat n, KnownNat m)
  => VS.Vector (CTriplet a)
  -> SparseMatrix n m a
fromVector = fromVectorSized (natToInt @n) (natToInt @m)

-- | Convert a sparse matrix to the list of triplets (row, col, val). Compressed elements will not be included.
toVector :: Elem a => SparseMatrix n m a -> VS.Vector (CTriplet a)
toVector m@(SparseMatrix fp) = Internal.performIO $ do
  let !size = nonZeros m
  tris <- VSM.new size
  withForeignPtr fp $ \p ->
    VSM.unsafeWith tris $ \q ->
      Internal.call $ Internal.sparse_toList p q (toC size)
  VS.unsafeFreeze tris

-- | Convert a sparse matrix to the list of triplets (row, col, val). Compressed elements will not be included.
toList :: Elem a => SparseMatrix n m a -> [(Int, Int, a)]
toList = Prelude.map fromC . VS.toList . toVector

-- | Construct a sparse matrix from a list of triples (row, val, col)
--
fromList :: (Elem a, KnownNat n, KnownNat m) => [(Int, Int, a)] -> SparseMatrix n m a
fromList = fromVector . VS.fromList . fmap toC

fromListSized :: (Elem a, KnownNat n, KnownNat m) => Int -> Int -> [(Int, Int, a)] -> SparseMatrix n m a
fromListSized _rows _cols = fromVectorSized _rows _cols . VS.fromList . Prelude.map Internal.toC

fromVectorSized :: forall n m a. (Elem a, KnownNat n, KnownNat m)
  => Int
  -> Int
  -> VS.Vector (CTriplet a)
  -> SparseMatrix n m a
fromVectorSized _rows _cols tris =
  let !c_rs = toC $! _rows
      !c_cs = toC $! _cols
      !len  = toC $! VS.length tris
  in Internal.performIO $ VS.unsafeWith tris $ \p ->
       alloca $ \pq -> do
         Internal.call $ Internal.sparse_fromList c_rs c_cs p len pq
         peek pq >>= _mk

-- | Construct asparse matrix of the given size from the storable vector of triplets (row, col, val)
--fromVector :: forall n m a. (Elem a, KnownNat n, KnownNat m)
--  => VS.Vector (CTriplet a)
--  -> SparseMatrix n m a
--fromVector tris =
--  let !c_rs = toC $! natToInt @n
--      !c_cs = toC $! natToInt @m
--      !len  = toC $! VS.length tris
--  in Internal.performIO $ VS.unsafeWith tris $ \p ->
--       alloca $ \pq -> do
--         Internal.call $ Internal.sparse_fromList c_rs c_cs p len pq
--         peek pq >>= _mk

-- | Convert a sparse matrix to a (n X m) dense list of values.
toDenseList :: forall n m a. (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> [[a]]
toDenseList mat = [[_unsafeCoeff row col mat | col <- [0 .. _unsafeCols mat - 1]] | row <- [0 .. _unsafeRows mat - 1]]

-- | Construct a sparsematrix from a two-dimensional list of values.
--   If the dimensions of the list do not match that of the list, 'Nothing' is returned.
--   Zero values will be compressed.
fromDenseList :: forall n m a. (Elem a, Eq a, KnownNat n, KnownNat m) => [[a]] -> Maybe (SparseMatrix n m a)
fromDenseList list =
  let _rows = List.length list
      _cols = List.foldl' max 0 $ List.map length list
  in if ((_rows /= (natToInt @n)) || (_cols /= (natToInt @m)))
    then Nothing
    else Just $ fromList $ do
      (row, vals) <- zip [0..] list
      (col, val) <- zip [0..] vals
      guard $ val /= 0
      return (row, col, val)

--  fromList $ do
-- (row, vals) <- zip [0..] list
-- (col, val) <- zip [0..] vals
-- guard $ val /= 0
-- return (row, col, val)
-- where
--   rows = List.length list
--   cols = List.foldl' max 0 $ List.map length list

-- | Construct a dense matrix from a sparse matrix
toMatrix :: (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> M.Matrix n m a
toMatrix (SparseMatrix fp) = Internal.performIO $ do
  m0 :: MM.IOMatrix n m a <- MM.new
  MM.unsafeWith m0 $ \_vals _rows _cols ->
    withForeignPtr fp $ \pm1 ->
      Internal.call $ Internal.sparse_toMatrix pm1 _vals _rows _cols
  M.unsafeFreeze m0

-- | Construct a sparse matrix from a dense matrix. zero-elements will be compressed.
fromMatrix :: (Elem a, KnownNat n, KnownNat m) => M.Matrix n m a -> SparseMatrix n m a
fromMatrix m1 = Internal.performIO $ alloca $ \pm0 ->
  M.unsafeWith m1 $ \_vals _rows _cols -> do
    Internal.call $ Internal.sparse_fromMatrix _vals _rows _cols pm0
    peek pm0 >>= _mk

-- | Yield an immutable copy of the mutable matrix
freeze :: (Elem a, PrimMonad p) => SMM.MSparseMatrix n m (PrimState p) a -> p (SparseMatrix n m a)
freeze (SMM.MSparseMatrix fp) = SparseMatrix <$> _clone fp

-- | Yield a mutable copy of the immutable matrix.
thaw :: (Elem a, PrimMonad p) => SparseMatrix n m a -> p (SMM.MSparseMatrix n m (PrimState p) a)
thaw (SparseMatrix fp) = SMM.MSparseMatrix <$> _clone fp

-- | Unsafely convert a mutable matrix to an immutable one without copying. The mutable matrix may not be used after this operation.
unsafeFreeze :: (Elem a, PrimMonad p) => SMM.MSparseMatrix n m (PrimState p) a -> p (SparseMatrix n m a)
unsafeFreeze (SMM.MSparseMatrix fp) = return $! SparseMatrix fp

-- | Unsafely convert an immutable matrix to a mutable one without copying. The immutable matrix may not be used after this operation.
unsafeThaw :: (Elem a, PrimMonad p) => SparseMatrix n m a -> p (SMM.MSparseMatrix n m (PrimState p) a)
unsafeThaw (SparseMatrix fp) = return $! SMM.MSparseMatrix fp

-- | Return a single row of the sparse matrix.
getRow :: forall n m r a. (Elem a, KnownNat n, KnownNat m, KnownNat r, r <= n, 1 <= n) => Row r -> SparseMatrix n m a -> SparseMatrix 1 m a
getRow row mat = block row (Col @0) (Row @1) (Col @m) mat

-- | Return a single column of the sparse matrix.
getCol :: forall n m c a. (Elem a, KnownNat n, KnownNat m, KnownNat c, c <= m, 1 <= m) => Col c -> SparseMatrix n m a -> SparseMatrix n 1 a
getCol col mat = block (Row @0) col (Row @n) (Col @1) mat

-- | Return all rows of the sparse matrix.
getRows :: forall n m a. (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> [SparseMatrix 1 m a]
getRows mat = fmap (flip unsafeGetRow mat) [0 .. (natToInt @n) - 1]

-- | Return all columns of a sparse matrix.
getCols :: forall n m a. (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> [SparseMatrix n 1 a]
getCols mat = fmap (flip unsafeGetCol mat) [0 .. (natToInt @m) - 1]

-- | Returns a single row of the sparse matrix. This is unsafe because it allows out-of-bounds access.
unsafeGetRow :: forall n m a. (Elem a, KnownNat n, KnownNat m) => Int -> SparseMatrix n m a -> SparseMatrix 1 m a
unsafeGetRow row mat = unsafeBlock row 0 1 (natToInt @m) mat

-- | Returns a single column of the sparse matrix. This is unsafe because it allows out-of-bounds access.
unsafeGetCol :: forall n m a. (Elem a, KnownNat n, KnownNat m) => Int -> SparseMatrix n m a -> SparseMatrix n 1 a
unsafeGetCol col mat = unsafeBlock 0 col (natToInt @n) 1 mat

-- | Get a sparse matrix from a list of columns.
--   This operation is unsafe because we cannot unify the length of the list with the type variable 'm'.
unsafeFromCols :: forall n m a. (Elem a, KnownNat n, KnownNat m) => [SparseMatrix n 1 a] -> SparseMatrix n m a
unsafeFromCols [] = fromListSized 0 0 []
unsafeFromCols _cols = fromListSized n m . concatMap toList . zipWith updateColIdx [0..] $ _cols
  where
    n = maximum . fmap _unsafeRows $ _cols
    m = length _cols

-- | Get a sparse matrix from a list of rows.
--   This operation is unsafe because we cannot unify the length of the list with the type variable 'n'.
unsafeFromRows :: forall n m a. (Elem a, KnownNat n, KnownNat m) => [SparseMatrix 1 m a] -> SparseMatrix n m a
unsafeFromRows [] = fromListSized 0 0 []
unsafeFromRows _rows = fromListSized n m . concatMap toList . zipWith updateRowIdx [0..] $ _rows
  where
    n = length _rows
    m = maximum . fmap _unsafeCols $ _rows

-- | Update the row indices of a sparse matrix.
updateRowIdx :: forall n m a. (KnownNat n, KnownNat m, Elem a) => Int -> SparseMatrix n m a -> SparseMatrix n m a
updateRowIdx row mat = fromListSized n m . fmap (\(_, j, v) -> (row, j, v)) . toList $ mat
  where
    n = row + 1
    m = _unsafeCols mat

-- | Update the column indices of a sparse matrix.
updateColIdx :: forall n m a. (KnownNat n, KnownNat m, Elem a) => Int -> SparseMatrix n m a -> SparseMatrix n m a
updateColIdx col mat = fromListSized n m . fmap (\(i, _, v) -> (i, col, v)) . toList $ mat
  where
    n = _unsafeRows mat
    m = col + 1

--fromRows :: forall n m a. (Elem a, KnownNat n, KnownNat m)
--  => [SparseMatrix 1 m a]
--  -> SparseMatrix n m a

_unop :: Storable b => (CSparseMatrixPtr a -> Ptr b -> IO CString) -> (b -> IO c) -> SparseMatrix n m a -> c
_unop f g (SparseMatrix fp) = Internal.performIO $
  withForeignPtr fp $ \p ->
    alloca $ \pq -> do
      Internal.call (f p pq)
      peek pq >>= g

_binop :: Storable b => (CSparseMatrixPtr a -> CSparseMatrixPtr a -> Ptr b -> IO CString) -> (b -> IO c) -> SparseMatrix n m a -> SparseMatrix n1 m1 a -> c
_binop f g (SparseMatrix fp1) (SparseMatrix fp2) = Internal.performIO $
  withForeignPtr fp1 $ \p1 ->
    withForeignPtr fp2 $ \p2 ->
      alloca $ \pq -> do
        Internal.call (f p1 p2 pq)
        peek pq >>= g

_getvec :: (Elem a, Storable b) => (Ptr (CSparseMatrix a) -> Ptr CInt -> Ptr (Ptr b) -> IO CString) -> SparseMatrix n m a -> VS.Vector b
_getvec f (SparseMatrix fm) = Internal.performIO $
  withForeignPtr fm $ \m ->
    alloca $ \ps ->
      alloca $ \pq -> do
        Internal.call $ f m ps pq
        s <- fromIntegral <$> peek ps
        q <- peek pq
        fr <- FC.newForeignPtr q $ touchForeignPtr fm
        pure $! VS.unsafeFromForeignPtr0 fr s

_clone :: (PrimMonad p, Elem a) => ForeignPtr (CSparseMatrix a) -> p (ForeignPtr (CSparseMatrix a))
_clone fp = unsafePrimToPrim $ withForeignPtr fp $ \p -> alloca $ \pq -> do
  Internal.call $ Internal.sparse_clone p pq
  q <- peek pq
  FC.newForeignPtr q $ Internal.call $ Internal.sparse_free q

_mk :: Elem a => Ptr (CSparseMatrix a) -> IO (SparseMatrix n m a)
_mk p = SparseMatrix <$> FC.newForeignPtr p (Internal.call $ Internal.sparse_free p)

-- | Number of rows in the sparse matrix
_unsafeRows :: forall n m a. (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> Int
{-# INLINE _unsafeRows #-}
_unsafeRows _ = natToInt @n

-- | Number of colums in the sparse matrix
_unsafeCols :: forall n m a. (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> Int
{-# INLINE _unsafeCols #-}
_unsafeCols _ = natToInt @m

_unsafeCoeff :: (Elem a, KnownNat n) => Int -> Int -> SparseMatrix n m a -> a
{-# INLINE _unsafeCoeff #-}
_unsafeCoeff !row !col (SparseMatrix fp) =
  let !c_row = toC row
      !c_col = toC col
  in Internal.performIO $ withForeignPtr fp $ \p -> alloca $ \pq -> do
       Internal.call $ Internal.sparse_coeff p c_row c_col pq
       fromC <$> peek pq

max_stars_repo_name: nilsalex/eigen

max_issues_repo_name: nilsalex/eigen

max_forks_repo_name: nilsalex/eigen


text:
{-# LANGUAGE BangPatterns #-}
module FourierSeries where

import           Control.Monad              as M
import           Data.Array.Repa            as R
import           Data.Complex
import           Data.List                  as L
import           Data.Vector.Storable       as VS
import           Data.Vector.Unboxed        as VU
import           FokkerPlanck.FourierSeries
import           Foreign.CUDA.BLAS          as BLAS
import           Foreign.CUDA.Driver        as CUDA
import           Image.IO
import           Numeric.LinearAlgebra      as NL hiding (All)
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Utils.BLAS
import           Utils.Distribution
import           Utils.Parallel
import           Utils.SimpsonRule
import           Utils.Time

main = do
  args@(gpuIDStr:sizeStr:deltaStr:periodStr:freqStr:angularFreqStr:radialFreqStr:stdStr:_) <-
    getArgs
  printCurrentTime . show $ args
  let gpuID = read gpuIDStr :: Int
      numPoints = read sizeStr :: Int
      delta = read deltaStr :: Double
      period = read periodStr :: Double
      freq = read freqStr :: Int
      numFreq = 2 * freq + 1
      angularFreq = read angularFreqStr :: Int
      radialFreq = read radialFreqStr :: Int
      std = read stdStr :: Double
      folderPath = "output/test/FourierSeries"
      len' = round $ (fromIntegral numPoints / delta :: Double)
      len =
        if odd len'
          then len'
          else len' + 1
      center = div len 2
      !numAngularFreq = 2 * angularFreq + 1
      !numRadialFreq = 2 * radialFreq + 1
      !numPinwheels = numAngularFreq * numRadialFreq
      pinwheel =
        fromFunction (Z :. len :. len :. numAngularFreq :. numRadialFreq) $ \(Z :. i :. j :. af :. rf) ->
          let x = delta * (fromIntegral $ i - center)
              y = delta * (fromIntegral $ j - center)
          in fourierMellin 0.5 (af - angularFreq) (rf - radialFreq) (x, y)
      rectangularHarmonics' =
        fromFunction (Z :. numFreq :. numFreq :. len :. len) $ \(Z :. xFreq :. yFreq :. i :. j) ->
          let !x = fromIntegral $ i - center
              !y = fromIntegral $ j - center
          in cis
               ((-2) * pi * delta *
                ((fromIntegral $ xFreq - freq) * x +
                 (fromIntegral $ yFreq - freq) * y) /
                period)
      simpsonWeights =
        computeUnboxedS $ computeWeightArrFromListOfShape [len, len] :: R.Array U DIM2 (Complex Double)
      inverseR2Harmonics' =
        fromFunction (Z :. numPoints :. numPoints :. numFreq :. numFreq) $ \(Z :. i :. j :. xFreq :. yFreq) ->
          cis
            (2 * pi *
             fromIntegral
               ((xFreq - freq) * (i - (div numPoints 2)) +
                (yFreq - freq) * (j - (div numPoints 2))) /
             period)
      !weight = (delta / 3) ^ 2 :+ 0
  createDirectoryIfMissing True folderPath
  createDirectoryIfMissing True (folderPath </> "Recons")
  createDirectoryIfMissing True (folderPath </> "Coefficients")
  print len
  !weightedPinwheels <-
    computeUnboxedP . R.traverse2 pinwheel simpsonWeights const $ \fP fS idx@(Z :. a :. b :. _ :. _) ->
      fP idx * fS (Z :. a :. b)
  !rectangularHarmonics <- computeUnboxedP rectangularHarmonics'
  !inverseR2Harmonics <- computeUnboxedP inverseR2Harmonics'
  -- printCurrentTime "CPU starts .."
  -- let !gaussianCoefficients =
  --       computeUnboxedS . R.fromFunction (Z :. numFreq :. numFreq) $ \(Z :. xFreq :. yFreq) ->
  --         weight *
  --         (gaussian2D
  --            (fromIntegral $ xFreq - freq)
  --            (fromIntegral $ yFreq - freq)
  --            std :+
  --          0)
  -- M.mapM_
  --   (\(af, rf) -> do
  --      let coefficients =
  --            sumS .
  --            sumS .
  --            R.traverse2
  --              rectangularHarmonics
  --              (R.slice weightedPinwheels (Z :. All :. All :. af :. rf))
  --              const $ \fR fP idx@(Z :. _ :. _ :. i :. j) ->
  --              fR idx * fP (Z :. i :. j)
  --          recon =
  --            sumS .
  --            sumS .
  --            R.traverse2
  --              inverseR2Harmonics
  --              (coefficients *^ gaussianCoefficients)
  --              const $ \fR fC idx@(Z :. _ :. _ :. i :. j) ->
  --              fR idx * fC (Z :. i :. j)
  --      plotImageRepaComplex
  --        (folderPath </> "Recon" </> printf "Recon%03d.png" (af * numRadialFreq + rf + 1)) .
  --        ImageRepa 8 . computeS . extend (Z :. (1 :: Int) :. All :. All) $
  --        recon)
  --   [(af, rf) | af <- [0 .. numAngularFreq - 1], rf <- [0 .. numRadialFreq - 1]]
  -- printCurrentTime "Done"
  -- coefficients <-
  --   (sumP . R.traverse2 rectangularHarmonics weightedPinwheel const $ \fR fP idx@(Z :. _ :. _ :. i :. j) ->
  --      fR idx * fP (Z :. i :. j)) >>=
  --   sumP
  --     gaussianCoefficients =
  --       computeUnboxedS . R.traverse coefficients id $ \f idx@(Z :. xFreq :. yFreq) ->
  --         weight * (f idx) *
  --         (gaussian2D
  --            (fromIntegral $ xFreq - freq)
  --            (fromIntegral $ yFreq - freq)
  --            std :+
  --          0)
  -- recon <-
  --   (sumP . R.traverse2 inverseR2Harmonics gaussianCoefficients const $ \fR fC idx@(Z :. _ :. _ :. i :. j) ->
  --      fR idx * fC (Z :. i :. j)) >>=
  --   sumP
  -- -- Repa
  -- let coefficients =
  --       sumS . sumS . R.traverse2 rectangularHarmonics weightedPinwheel const $ \fR fP idx@(Z :. _ :. _ :. i :. j) ->
  --         fR idx * fP (Z :. i :. j)
  --     gaussianCoefficients =
  --       computeUnboxedS . R.traverse coefficients id $ \f idx@(Z :. xFreq :. yFreq) ->
  --         weight * (f idx) *
  --         (gaussian2D
  --            (fromIntegral $ xFreq - freq)
  --            (fromIntegral $ yFreq - freq)
  --            std :+
  --          0)
  --     recon =
  --       sumS . sumS . R.traverse2 inverseR2Harmonics gaussianCoefficients const $ \fR fC idx@(Z :. _ :. _ :. i :. j) ->
  --         fR idx * fC (Z :. i :. j)
  -- -- HMatrix
  -- let rectangularHarmonicsMat =
  --       ((numFreq ^ 2) >< (len ^ 2)) . R.toList $ rectangularHarmonics
  --     inverseR2HarmonicsMat =
  --       ((numPoints ^ 2) >< (numFreq ^ 2)) . R.toList $ inverseR2Harmonics
  --     weightedPinwheelVec = NL.fromList . R.toList $ weightedPinwheel
  --     gaussianCoefficientsVec =
  --       NL.fromList . R.toList . R.fromFunction (Z :. numFreq :. numFreq) $ \(Z :. xFreq :. yFreq) ->
  --         weight *
  --         (gaussian2D
  --            (fromIntegral $ xFreq - freq)
  --            (fromIntegral $ yFreq - freq)
  --            std :+
  --          0)
  --     coefficientsVec = rectangularHarmonicsMat #> weightedPinwheelVec
  --     recon =
  --       fromListUnboxed (Z :. numPoints :. numPoints) . NL.toList $
  --       inverseR2HarmonicsMat #> (coefficientsVec * gaussianCoefficientsVec)
  -- --BLAS mxv
  -- let rectangularHarmonicsMat = VU.convert . toUnboxed $ rectangularHarmonics
  --     inverseR2HarmonicsMat = VU.convert . toUnboxed $ inverseR2Harmonics
  --     weightedPinwheelVec = VU.convert . toUnboxed $ weightedPinwheel
  --     gaussianCoefficientsVec =
  --       VU.convert .
  --       toUnboxed . computeS . R.fromFunction (Z :. numFreq :. numFreq) $ \(Z :. xFreq :. yFreq) ->
  --         weight *
  --         (gaussian2D
  --            (fromIntegral $ xFreq - freq)
  --            (fromIntegral $ yFreq - freq)
  --            std :+
  --          0)
  -- coefficientsVec <-
  --   mxv (numFreq ^ 2) (len ^ 2) rectangularHarmonicsMat weightedPinwheelVec
  -- reconVec <-
  --   mxv
  --     (numPoints ^ 2)
  --     (numFreq ^ 2)
  --     inverseR2HarmonicsMat
  --     (VS.zipWith (*) coefficientsVec gaussianCoefficientsVec)
  -- let recon = fromUnboxed (Z :. numPoints :. numPoints) . VS.convert $ reconVec
  --BLAS mxm
  -- let !rectangularHarmonicsMat = VU.convert . toUnboxed $ rectangularHarmonics
  --     !inverseR2HarmonicsMat = VU.convert . toUnboxed $ inverseR2Harmonics
  --     !weightedPinwheelVec = VU.convert . toUnboxed $ weightedPinwheels
  --     !gaussianCoefficients =
  --       computeUnboxedS . R.fromFunction (Z :. numFreq :. numFreq) $ \(Z :. xFreq :. yFreq) ->
  --         weight *
  --         (gaussian2D
  --            (fromIntegral $ xFreq - freq)
  --            (fromIntegral $ yFreq - freq)
  --            std :+
  --          0)
  -- printCurrentTime "Start.."
  -- coefficientsVec <-
  --   gemmBLAS
  --     (numFreq ^ 2)
  --     numPinwheels
  --     (len ^ 2)
  --     rectangularHarmonicsMat
  --     weightedPinwheelVec
  -- printCurrentTime "coefficients done"
  -- let coefficientsArr =
  --       fromUnboxed (Z :. numFreq :. numFreq :. numAngularFreq :. numRadialFreq) .
  --       VS.convert $
  --       coefficientsVec
  --     weightedCoefficientsArr =
  --       computeS .
  --       R.traverse2
  --         coefficientsArr
  --         gaussianCoefficients
  --         const $ \f1 f2 idx@(Z :. i :. j :. _ :. _) ->
  --         f1 idx * f2 (Z :. i :. j)
  --     weightedCoefficients = VS.convert . toUnboxed $ weightedCoefficientsArr
  -- M.mapM_
  --   (\(af, rf) ->
  --      plotImageRepaComplex
  --        (folderPath </> "Coefficients" </>
  --         printf
  --           "Coef(%d,%d).png"
  --           (af - div numAngularFreq 2)
  --           (rf - div numRadialFreq 2)) .
  --      ImageRepa 8 .
  --      computeS .
  --      extend (Z :. (1 :: Int) :. All :. All) . R.slice coefficientsArr $
  --      (Z :. All :. All :. af :. rf))
  --   [(af, rf) | af <- [0 .. numAngularFreq - 1], rf <- [0 .. numRadialFreq - 1]]
  -- printCurrentTime "Done"
  -- reconVec <-
  --   gemmBLAS
  --     (numPoints ^ 2)
  --     numPinwheels
  --     (numFreq ^ 2)
  --     inverseR2HarmonicsMat
  --     weightedCoefficients
  -- let !recon =
  --       fromUnboxed
  --         (Z :. numPoints :. numPoints :. numAngularFreq :. numRadialFreq) .
  --       VS.convert $
  --       reconVec
  -- M.mapM_
  --   (\(af, rf) ->
  --      plotImageRepaComplex
  --        (folderPath </> "Recons" </>
  --         printf
  --           "Recon(%d,%d).png"
  --           (af - div numAngularFreq 2)
  --           (rf - div numRadialFreq 2)) .
  --      ImageRepa 8 .
  --      computeS . extend (Z :. (1 :: Int) :. All :. All) . R.slice recon $
  --      (Z :. All :. All :. af :. rf))
  --   [(af, rf) | af <- [0 .. numAngularFreq - 1], rf <- [0 .. numRadialFreq - 1]]
  -- printCurrentTime "Done"
  --BLAS mxmGPU
  let !rectangularHarmonicsMat = VU.convert . toUnboxed $ rectangularHarmonics
      !inverseR2HarmonicsMat = VU.convert . toUnboxed $ inverseR2Harmonics
      !weightedPinwheelVec = VU.convert . toUnboxed $ weightedPinwheels
      !gaussianCoefficients =
        computeUnboxedS . R.fromFunction (Z :. numFreq :. numFreq) $ \(Z :. xFreq :. yFreq) ->
          weight *
          (gaussian2D
             (fromIntegral $ xFreq - freq)
             (fromIntegral $ yFreq - freq)
             std :+
           0)
  initialise []
  dev <- device gpuID
  ctx <- CUDA.create dev []
  handle <- BLAS.create
  printCurrentTime "Start.."
  coefficientsVec <-
    gemmCuBLAS11
      handle
      (numFreq ^ 2)
      numPinwheels
      (len ^ 2)
      rectangularHarmonicsMat
      weightedPinwheelVec
  let coefficientsArr =
        fromUnboxed (Z :. numFreq :. numFreq :. numAngularFreq :. numRadialFreq) .
        VS.convert $
        coefficientsVec
      weightedCoefficients =
        VU.convert .
        toUnboxed .
        computeS . R.traverse2 coefficientsArr gaussianCoefficients const $ \f1 f2 idx@(Z :. i :. j :. _ :. _) ->
          f1 idx * f2 (Z :. i :. j)
  printCurrentTime "coefficients done"
  reconVec <-
    gemmCuBLAS11
      handle
      (numPoints ^ 2)
      numPinwheels
      (numFreq ^ 2)
      inverseR2HarmonicsMat
      weightedCoefficients
  let !recon =
        fromUnboxed (Z :. numPoints :. numPoints :. numPinwheels) . VS.convert $
        reconVec
  M.mapM_
    (\i ->
       plotImageRepaComplex
         (folderPath </> "Recons" </> printf "Recon%03d.png" (i + 1)) .
       ImageRepa 8 .
       computeS . extend (Z :. (1 :: Int) :. All :. All) . R.slice recon $
       (Z :. All :. All :. i))
    [0 .. numPinwheels - 1]
  printCurrentTime "Done"
  M.mapM_
    (\(af, rf) ->
       plotImageRepaComplex
         (folderPath </> "Coefficients" </>
          printf
            "Coef(%d,%d).png"
            (af - div numAngularFreq 2)
            (rf - div numRadialFreq 2)) .
       ImageRepa 8 .
       computeS .
       extend (Z :. (1 :: Int) :. All :. All) . R.slice coefficientsArr $
       (Z :. All :. All :. af :. rf))
    [(af, rf) | af <- [0 .. numAngularFreq - 1], rf <- [0 .. numRadialFreq - 1]]
  -- plotImageRepaComplex (folderPath </> "Grating.png") .
  --   ImageRepa 8 .
  --   computeS .
  --   extend (Z :. (1 :: Int) :. All :. All) . R.slice rectangularHarmonics $
  --   (Z :. (0 :: Int) :. freq :. All :. All)
  -- plotImageRepaComplex (folderPath </> "Pinwheel.png") .
  --   ImageRepa 8 .
  --   computeS . fromFunction (Z :. (1 :: Int) :. numPoints :. numPoints) $ \(Z :. _ :. i :. j) ->
  --   fourierMellin
  --     angularFreq
  --     radialFreq
  --     ( (fromIntegral $ i - (div numPoints 2))
  --     , (fromIntegral $ j - (div numPoints 2)))

ext: hs

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE MagicHash, UnboxedTuples #-}

module Render (render, parRender) where

import Control.Parallel.Strategies

import Numeric.Vector
import Numeric.Scalar (fromScalar, scalar)
import Types
import Control.Monad.State
import Camera (rayFrom)
import Raycast (raycast)
import Shader (reflectance, sampleBRDF, probability)
import Scene (makeContext)
import VecUtil (unpack3, vec3Of)
import System.Random
import SceneRandom
import Lights

map2 :: (a -> b) -> [[a]] -> [[b]]
map2 f = map (map f)

vec23 :: Vec2d -> (Double, Double, Double)
vec23 v = case unpackV2# v of
            (# u, v #) -> (u, v, 1)

possibly :: (a -> b) -> b -> Maybe a -> b
possibly _ def Nothing = def
possibly f _ (Just x) = f x

shade :: Int -> State SceneContext Vec3d
shade bounceCount = get >>= (\ctx -> if bounceCount > s_getBounces (ss_getScene ctx) then return (vec3Of 0) else do
  let hit = ss_getHit ctx
  let brdf = rh_getShader $ hit
  col <- reflectance brdf
  let scene = ss_getScene ctx
  let lights = s_getLights scene
  let norm = rh_getNormal hit

  mReflectRay <- sampleBRDF brdf
  (reflectColor, reflectProbability) <- case mReflectRay of
    Nothing -> return (0, 0)
    Just reflectRay -> do
      prob <- probability brdf (r_getDir reflectRay)
      case raycast scene reflectRay of
        Nothing -> return $ (s_getSkyColor scene, )prob
        Just hit -> do
          let (v, s) = runState (shade (bounceCount + 1)) $ makeContext scene (ss_getGen ctx) hit
          put s
          return (v, prob) -- 1 / (2 * PI)

  let newRayPos = (rh_getPos hit) + (norm * (vec3Of 0.01))

  lighting <- sequence $ map (sampleLight newRayPos norm) lights
  let lighting' = filter ((>0) . fst) lighting
  lightProbabilities <- sequence $ map (probability brdf . snd) lighting'
  let lightColors = map fst lighting'

  let probabilitySum = reflectProbability + sum lightProbabilities

  let adjustedReflectColor = reflectColor * vec3Of reflectProbability
  let adjustedIncomingLighting = sum $ zipWith (*) lightColors (map vec3Of lightProbabilities)

  let incomingLight = adjustedReflectColor + adjustedIncomingLighting

  let lit = col * incomingLight

  reflectColor <- reflectance brdf
  return lit
  )

samplePixel :: Scene -> (Int, Int) -> Int -> Vec3d
samplePixel s (x, y) sample = let w = s_getWidth s
                                  h = s_getHeight s
                                  cam = s_getCamera s
                                  skyCol = s_getSkyColor s
                                  -- Random seed is an arbitrary formula based on pixel location to form distinct results
                                  gen = mkStdGen $ (x + y * w + x * (x + y) + (s_getSeed s)) * sample
                              in  possibly (fst . runState (shade 0) . makeContext s gen) (skyCol / (2 * pi)) $
                                  raycast s $
                                  rayFrom cam $
                                  (fromIntegral x / fromIntegral w * 2.0 - 1.0, (fromIntegral h / fromIntegral w) - fromIntegral y / fromIntegral w * 2.0)

renderPixel :: Scene -> (Int, Int) -> Vec3d
renderPixel s p = let numSamples = s_getSamples s
                      samples = map (samplePixel s p)  [1..numSamples]
                  -- Multiply by 2pi for monte-carlo integration because the volume of the integration bounds is 2pi
                  -- (integration bounds form a hemisphere)
                  in  sum samples / fromIntegral numSamples * vec3Of (2 * pi)

render :: Scene -> [[(Double, Double, Double)]]
render s = let w = s_getWidth s
               h = s_getHeight s
               cam = s_getCamera s
               skyCol = s_getSkyColor s
               seed = 1
               gen = mkStdGen seed
          in  map2 unpack3 $
              map2 (renderPixel s) $
              [ [ (x, y) | x <- [1..w] ] | y <- [1..h] ]

-- Found how to evaluate the list in parallel from
-- https://stackoverflow.com/a/5606176
parRender :: Scene -> [[(Double, Double, Double)]]
parRender s = (render s) `using` parList rdeepseq


text:
{-# LANGUAGE ScopedTypeVariables #-}

import           Control.Monad.Primitive
import           Control.Monad
import           Foreign.C.Types
import           Foreign.Ptr
import           Unsafe.Coerce
import           Data.Complex

import qualified Data.Vector.Generic               as VG
import qualified Data.Vector.Generic.Mutable       as VGM
import qualified Data.Vector.Storable              as VS
import qualified Data.Vector.Storable.Mutable      as VSM
import qualified Data.Vector.Fusion.Bundle         as VFS
import qualified Data.Vector.Fusion.Stream.Monadic as VFSM

import           Foreign.Storable.Complex
import           Criterion.Main

import           SDR.FilterInternal
import           SDR.Util
import           SDR.CPUID

theBench :: IO ()
theBench = do
    --Setup
    let size       =  16384
        numCoeffs  =  128
        num        =  size - numCoeffs + 1
        decimation =  4
        interpolation = 3
        numCoeffsDiv2  =  64

        coeffsList :: [Float]
        coeffsList = take numCoeffs [0 ..]
        coeffs    :: VS.Vector Float
        coeffs    =  VG.fromList $ take numCoeffs [0 ..]
        coeffsSym    :: VS.Vector Float
        coeffsSym    =  VG.fromList $ take numCoeffsDiv2 [0 ..]
        inBuf     :: VS.Vector Float
        inBuf     =  VG.fromList $ take size [0 ..]
        inBufComplex :: VS.Vector (Complex Float)
        inBufComplex =  VG.fromList $ take size $ do
            i <- [0..]
            return $ i :+ i
        inBufRTLSDR  :: VS.Vector CUChar
        inBufRTLSDR  =  VG.fromList $ take size [0 ..]
        inBufBladeRF :: VS.Vector CShort
        inBufBladeRF =  VG.fromList $ take size [0 ..]

        numConv   = 16386
        inBufConv :: VS.Vector CUChar
        inBufConv = VG.fromList $ take size $ concat $ repeat [0 .. 255]

        duplicate :: [a] -> [a]
        duplicate = concatMap func
            where func x = [x, x]

        coeffs2 :: VS.Vector Float
        coeffs2 =  VG.fromList $ duplicate $ take numCoeffs [0 ..]

    outBuf        :: VS.MVector RealWorld Float <- VGM.new size
    outBufComplex :: VS.MVector RealWorld (Complex Float) <- VGM.new size

    info <- getCPUInfo

    let hasFeatures :: [(CPUInfo -> Bool, a)] -> [a]
        hasFeatures = map snd . filter (($ info) . fst)

    resampler3 <- resampleCRR2   interpolation decimation coeffsList
    resampler4 <- resampleCSSERR interpolation decimation coeffsList
    resampler5 <- resampleCAVXRR interpolation decimation coeffsList

    resampler3C <- resampleCRC    interpolation decimation coeffsList
    resampler4C <- resampleCSSERC interpolation decimation coeffsList
    resampler5C <- resampleCAVXRC interpolation decimation coeffsList

    --Benchmarks
    defaultMain [
            bgroup "filter" [
                bgroup "real" $ hasFeatures [
                    (const True, bench "highLevel"   $ nfIO $ filterHighLevel          coeffs    num inBuf outBuf),
                    (const True, bench "imperative1" $ nfIO $ filterImperative1        coeffs    num inBuf outBuf),
                    (const True, bench "imperative2" $ nfIO $ filterImperative2        coeffs    num inBuf outBuf),
                    (const True, bench "c"           $ nfIO $ filterCRR                coeffs    num inBuf outBuf),
                    (hasSSE42,   bench "cSSE"        $ nfIO $ filterCSSERR             coeffs    num inBuf outBuf),
                    (hasSSE42,   bench "cSSESym"     $ nfIO $ filterCSSESymmetricRR    coeffsSym num inBuf outBuf),
                    (hasAVX,     bench "cAVX"        $ nfIO $ filterCAVXRR             coeffs    num inBuf outBuf),
                    (hasAVX,     bench "cAVXSym"     $ nfIO $ filterCAVXSymmetricRR    coeffsSym num inBuf outBuf)
                ],
                bgroup "complex" $ hasFeatures [
                    (const True, bench "highLevel"   $ nfIO $ filterHighLevel          coeffs    num inBufComplex outBufComplex),
                    (const True, bench "c"           $ nfIO $ filterCRC                coeffs    num inBufComplex outBufComplex),
                    (hasSSE42,   bench "cSSE"        $ nfIO $ filterCSSERC             coeffs2   num inBufComplex outBufComplex),
                    (hasSSE42,   bench "cSSE2"       $ nfIO $ filterCSSERC2            coeffs    num inBufComplex outBufComplex),
                    (hasSSE42,   bench "cSSESym"     $ nfIO $ filterCSSESymmetricRC    coeffsSym num inBufComplex outBufComplex),
                    (hasAVX,     bench "cAVX"        $ nfIO $ filterCAVXRC             coeffs2   num inBufComplex outBufComplex),
                    (hasAVX,     bench "cAVX2"       $ nfIO $ filterCAVXRC2            coeffs    num inBufComplex outBufComplex),
                    (hasAVX,     bench "cAVXSym"     $ nfIO $ filterCAVXSymmetricRC    coeffsSym num inBufComplex outBufComplex)
                ]
            ],
            bgroup "decimate" [
                bgroup "real" $ hasFeatures [
                    (const True, bench "highLevel"   $ nfIO $ decimateHighLevel        decimation coeffs    (num `quot` decimation) inBuf outBuf),
                    (const True, bench "c"           $ nfIO $ decimateCRR              decimation coeffs    (num `quot` decimation) inBuf outBuf),
                    (hasSSE42,   bench "cSSE"        $ nfIO $ decimateCSSERR           decimation coeffs    (num `quot` decimation) inBuf outBuf),
                    (hasSSE42,   bench "cSSESym"     $ nfIO $ decimateCSSESymmetricRR  decimation coeffsSym (num `quot` decimation) inBuf outBuf),
                    (hasAVX,     bench "cAVX"        $ nfIO $ decimateCAVXRR           decimation coeffs    (num `quot` decimation) inBuf outBuf),
                    (hasAVX,     bench "cAVXSym"     $ nfIO $ decimateCAVXSymmetricRR  decimation coeffsSym (num `quot` decimation) inBuf outBuf)
                ],
                bgroup "complex" $ hasFeatures [
                    (const True, bench "highLevel"   $ nfIO $ decimateHighLevel        decimation coeffs    (num `quot` decimation) inBufComplex outBufComplex),
                    (const True, bench "c"           $ nfIO $ decimateCRC              decimation coeffs    (num `quot` decimation) inBufComplex outBufComplex),
                    (hasSSE42,   bench "cSSE"        $ nfIO $ decimateCSSERC           decimation coeffs2   (num `quot` decimation) inBufComplex outBufComplex),
                    (hasSSE42,   bench "cSSE2"       $ nfIO $ decimateCSSERC2          decimation coeffs    (num `quot` decimation) inBufComplex outBufComplex),
                    (hasSSE42,   bench "cSSESym"     $ nfIO $ decimateCSSESymmetricRC  decimation coeffsSym (num `quot` decimation) inBufComplex outBufComplex),
                    (hasAVX,     bench "cAVX"        $ nfIO $ decimateCAVXRC           decimation coeffs2   (num `quot` decimation) inBufComplex outBufComplex),
                    (hasAVX,     bench "cAVX2"       $ nfIO $ decimateCAVXRC2          decimation coeffs    (num `quot` decimation) inBufComplex outBufComplex),
                    (hasAVX,     bench "cAVXSym"     $ nfIO $ decimateCAVXSymmetricRC  decimation coeffsSym (num `quot` decimation) inBufComplex outBufComplex)
                ]
            ],
            bgroup "resample" [
                bgroup "real" $ hasFeatures [
                    (const True, bench "highLevel"   $ nfIO $ resampleHighLevel        interpolation decimation coeffs 0 (num `quot` decimation) inBuf outBuf),
                    (const True, bench "c"           $ nfIO $ resampleCRR              (num `quot` decimation) interpolation decimation 0 coeffs inBuf outBuf),
                    (const True, bench "c2"          $ nfIO $ resampler3               (num `quot` decimation) 0 inBuf outBuf),
                    (hasSSE42,   bench "cSSE"        $ nfIO $ resampler4               (num `quot` decimation) 0 inBuf outBuf),
                    (hasAVX,     bench "cAVX"        $ nfIO $ resampler5               (num `quot` decimation) 0 inBuf outBuf)
                ],
                bgroup "complex" $ hasFeatures [
                    (const True, bench "highLevel"   $ nfIO $ resampleHighLevel        interpolation decimation coeffs 0 (num `quot` decimation) inBufComplex outBufComplex),
                    (const True, bench "c"           $ nfIO $ resampler3C               (num `quot` decimation) 0 inBufComplex outBufComplex),
                    (const True, bench "SSE"         $ nfIO $ resampler4C               (num `quot` decimation) 0 inBufComplex outBufComplex),
                    (const True, bench "AVX"         $ nfIO $ resampler5C               (num `quot` decimation) 0 inBufComplex outBufComplex)
                ]
            ],
            bgroup "scaling" $ hasFeatures [
                (const True, bench "c"               $ nfIO $ scaleC    0.3 inBuf outBuf),
                (hasSSE42,   bench "cSSE"            $ nfIO $ scaleCSSE 0.3 inBuf outBuf),
                (hasAVX,     bench "cAVX"            $ nfIO $ scaleCAVX 0.3 inBuf outBuf)
            ],
            bgroup "conversion" [
                bgroup "RTLSDR" $ hasFeatures [
                    (const True, bench "h"    $ nf (interleavedIQUnsigned256ToFloat :: VS.Vector CUChar -> VS.Vector (Complex Float)) inBufRTLSDR),
                    (const True, bench "c"    $ nf interleavedIQUnsignedByteToFloat    inBufRTLSDR),
                    (hasSSE42,   bench "cSSE" $ nf interleavedIQUnsignedByteToFloatSSE inBufRTLSDR),
                    (hasAVX2,    bench "cAVX" $ nf interleavedIQUnsignedByteToFloatAVX inBufRTLSDR)
                ],
                bgroup "BladeRF" $ hasFeatures [
                    (const True, bench "h"    $ nf (interleavedIQSigned2048ToFloat :: VS.Vector CShort -> VS.Vector (Complex Float))  inBufBladeRF),
                    (const True, bench "c"    $ nf interleavedIQSignedWordToFloat    inBufBladeRF),
                    (hasSSE42,   bench "cSSE" $ nf interleavedIQSignedWordToFloatSSE inBufBladeRF),
                    (hasAVX2,    bench "cAVX" $ nf interleavedIQSignedWordToFloatAVX inBufBladeRF)
                ]
            ]
        ]

main = theBench

max_stars_repo_name: adamwalker/sdr

max_issues_repo_name: peixian/sdr

max_forks_repo_name: adamwalker/sdr


text:
{-# LANGUAGE FlexibleInstances #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}

module Data.Configurator.Instances () where

import Data.Configurator.Types.Internal
import Data.Complex (Complex)
import Data.Fixed (Fixed, HasResolution)
import Data.Int (Int8, Int16, Int32, Int64)
import Data.Text.Encoding (encodeUtf8)
import Data.Ratio (Ratio, denominator, numerator)
import Data.Word (Word8, Word16, Word32, Word64)
import Foreign.C.Types (CDouble, CFloat)
import qualified Data.ByteString as B
import qualified Data.ByteString.Lazy as LB
import qualified Data.Text as T
import qualified Data.Text.Lazy as L

instance Configured Value where
    convert = Just

instance Configured Bool where
    convert (Bool v) = Just v
    convert _        = Nothing

convertNumberToNum :: (Num a) => Value -> Maybe a
convertNumberToNum (Number r)
    | denominator r == 1 = Just $ fromInteger $ numerator r
convertNumberToNum _ = Nothing

instance Configured Int where
    convert = convertNumberToNum

instance Configured Integer where
    convert = convertNumberToNum

instance Configured Int8 where
    convert = convertNumberToNum

instance Configured Int16 where
    convert = convertNumberToNum

instance Configured Int32 where
    convert = convertNumberToNum

instance Configured Int64 where
    convert = convertNumberToNum

instance Configured Word where
    convert = convertNumberToNum

instance Configured Word8 where
    convert = convertNumberToNum

instance Configured Word16 where
    convert = convertNumberToNum

instance Configured Word32 where
    convert = convertNumberToNum

instance Configured Word64 where
    convert = convertNumberToNum

convertNumberToFractional :: (Fractional a) => Value -> Maybe a
convertNumberToFractional (Number r) = Just $ fromRational r
convertNumberToFractional _ = Nothing

instance Configured Double where
    convert = convertNumberToFractional

instance Configured Float where
    convert = convertNumberToFractional

instance Configured CDouble where
    convert = convertNumberToFractional

instance Configured CFloat where
    convert = convertNumberToFractional

instance Integral a => Configured (Ratio a) where
    convert = convertNumberToFractional

instance RealFloat a => Configured (Complex a) where
    convert = convertNumberToFractional

instance HasResolution a => Configured (Fixed a) where
    convert = convertNumberToFractional

instance Configured T.Text where
    convert (String v) = Just v
    convert _          = Nothing

instance Configured Char where
    convert (String txt) | T.length txt == 1 = Just $ T.head txt
    convert _ = Nothing

    convertList = fmap T.unpack . convert

instance Configured L.Text where
    convert = fmap L.fromStrict . convert

instance Configured B.ByteString where
    convert = fmap encodeUtf8 . convert

instance Configured LB.ByteString where
    convert = fmap (LB.fromChunks . (:[])) . convert

instance (Configured a, Configured b) => Configured (a,b) where
    convert (List [a,b]) = (,) <$> convert a <*> convert b
    convert _            = Nothing

instance (Configured a, Configured b, Configured c) => Configured (a,b,c) where
    convert (List [a,b,c]) = (,,) <$> convert a <*> convert b <*> convert c
    convert _              = Nothing

instance (Configured a, Configured b, Configured c, Configured d)
    => Configured (a,b,c,d) where
    convert (List [a,b,c,d]) = (,,,) <$> convert a <*> convert b <*> convert c
                                     <*> convert d
    convert _                = Nothing

max_stars_repo_name: unisonweb/configurator

max_issues_repo_name: unisonweb/configurator

max_forks_repo_name: unisonweb/configurator


text:
module Ranks
  ( initWeights,
    ilsrPairwise,
  )
where

import Data.List ((!!), iterate)
import Import hiding ((<>), Vector, fromList)
import Numeric.LinearAlgebra

mean :: Vector R -> R
mean v = sumElements v / fromIntegral (size v)

expTransform :: Vector R -> Vector R
expTransform weights = weights' / scalar (mean weights')
  where
    weights' = exp (weights - scalar (mean weights))

logTransform :: Vector R -> Vector R
logTransform weights = params - scalar (mean params)
  where
    params = log weights

initWeights :: Int -> Vector R
initWeights = konst 1

iterLsr :: Vector R -> Vector R
iterLsr = expTransform

initChain :: Int -> R -> Matrix R
initChain n alpha = scalar alpha * konst 1 (n, n)

lsrPairwise :: Matrix R -> R -> Vector R -> Vector R
lsrPairwise bm alpha weights = logTransform (statdist chain)
  where
    n = size weights
    denom = asRow weights
    chain' = initChain n alpha + tr' bm / (denom + tr' denom)
    chain = chain' - (diagl . fmap sumElements . toRows) chain'

ilsrPairwise :: Matrix R -> R -> Vector R
ilsrPairwise bm alpha = iterate (lsr . iterLsr) (lsr weights) !! 100
  where
    weights = initWeights (rows bm)
    lsr = lsrPairwise bm alpha

statdist :: Matrix R -> Vector R
statdist chain = res / scalar (mean res)
  where
    (_, u, _, _) = lu (tr' chain)
    left = u ?? (DropLast 1, DropLast 1)
    right = negate (u ?? (DropLast 1, TakeLast 1))
    res = vjoin [flatten (triSolve Upper left right), 1.0]

max_stars_repo_name: jsza/isaacranks

max_issues_repo_name: jsza/isaacranks

max_forks_repo_name: jsza/isaacranks


text:
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE RoleAnnotations #-}
{-# LANGUAGE EmptyDataDecls #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE Trustworthy #-}
{-# LANGUAGE Rank2Types #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE CPP #-}
#if __GLASGOW_HASKELL__ >= 800
{-# LANGUAGE UndecidableSuperClasses #-}
#endif
#if __GLASGOW_HASKELL__ >= 708 && __GLASGOW_HASKELL__ < 710
{-# LANGUAGE NullaryTypeClasses #-}
#endif
-----------------------------------------------------------------------------
-- |
-- Module      :  Data.Constraint
-- Copyright   :  (C) 2011-2015 Edward Kmett,
-- License     :  BSD-style (see the file LICENSE)
--
-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
-- Stability   :  experimental
-- Portability :  non-portable
--
-- @ConstraintKinds@ made type classes into types of a new kind, @Constraint@.
--
-- @
-- 'Eq' :: * -> 'Constraint'
-- 'Ord' :: * -> 'Constraint'
-- 'Monad' :: (* -> *) -> 'Constraint'
-- @
--
-- The need for this extension was first publicized in the paper
--
-- <http://research.microsoft.com/pubs/67439/gmap3.pdf Scrap your boilerplate with class: extensible generic functions>
--
-- by Ralf Lämmel and Simon Peyton Jones in 2005, which shoehorned all the
-- things they needed into a custom 'Sat' typeclass.
--
-- With @ConstraintKinds@ we can put into code a lot of tools for manipulating
-- these new types without such awkward workarounds.
----------------------------------------------------------------------------
module Data.Constraint
  (
  -- * The Kind of Constraints
    Constraint
  -- * Dictionary
  , Dict(Dict)
  , withDict
  -- * Entailment
  , (:-)(Sub)
  , (\\)
  , weaken1, weaken2, contract
  , strengthen1, strengthen2
  , (&&&), (***)
  , trans, refl
  , Bottom(no)
  , top, bottom
  -- * Dict is fully faithful
  , mapDict
  , unmapDict
  -- * Reflection
  , Class(..)
  , (:=>)(..)
  ) where
import Control.Applicative
import Control.Category
import Control.DeepSeq
import Control.Monad
import Data.Complex
import Data.Ratio
import Data.Semigroup
import Data.Data
import qualified GHC.Exts as Exts (Any)
import GHC.Exts (Constraint)
import Data.Bits (Bits)
import Data.Functor.Identity (Identity)
#if MIN_VERSION_base(4,8,0)
import Numeric.Natural (Natural)
#endif
#if !MIN_VERSION_base(4,8,0)
import Data.Word (Word)
#endif

-- | Values of type @'Dict' p@ capture a dictionary for a constraint of type @p@.
--
-- e.g.
--
-- @
-- 'Dict' :: 'Dict' ('Eq' 'Int')
-- @
--
-- captures a dictionary that proves we have an:
--
-- @
-- instance 'Eq' 'Int
-- @
--
-- Pattern matching on the 'Dict' constructor will bring this instance into scope.
--
data Dict :: Constraint -> * where
  Dict :: a => Dict a
  deriving Typeable

instance (Typeable p, p) => Data (Dict p) where
  gfoldl _ z Dict = z Dict
  toConstr _ = dictConstr
  gunfold _ z c = case constrIndex c of
    1 -> z Dict
    _ -> error "gunfold"
  dataTypeOf _ = dictDataType

dictConstr :: Constr
dictConstr = mkConstr dictDataType "Dict" [] Prefix

dictDataType :: DataType
dictDataType = mkDataType "Data.Constraint.Dict" [dictConstr]

deriving instance Eq (Dict a)
deriving instance Ord (Dict a)
deriving instance Show (Dict a)

instance NFData (Dict c) where
  rnf Dict = ()

-- | From a 'Dict', takes a value in an environment where the instance
-- witnessed by the 'Dict' is in scope, and evaluates it.
--
-- Essentially a deconstruction of a 'Dict' into its continuation-style
-- form.
--
withDict :: Dict a -> (a => r) -> r
withDict d r = case d of
                 Dict -> r

infixr 9 :-

-- | This is the type of entailment.
--
-- @a ':-' b@ is read as @a@ \"entails\" @b@.
--
-- With this we can actually build a category for 'Constraint' resolution.
--
-- e.g.
--
-- Because @'Eq' a@ is a superclass of @'Ord' a@, we can show that @'Ord' a@
-- entails @'Eq' a@.
--
-- Because @instance 'Ord' a => 'Ord' [a]@ exists, we can show that @'Ord' a@
-- entails @'Ord' [a]@ as well.
--
-- This relationship is captured in the ':-' entailment type here.
--
-- Since @p ':-' p@ and entailment composes, ':-' forms the arrows of a
-- 'Category' of constraints. However, 'Category' only became sufficiently
-- general to support this instance in GHC 7.8, so prior to 7.8 this instance
-- is unavailable.
--
-- But due to the coherence of instance resolution in Haskell, this 'Category'
-- has some very interesting properties. Notably, in the absence of
-- @IncoherentInstances@, this category is \"thin\", which is to say that
-- between any two objects (constraints) there is at most one distinguishable
-- arrow.
--
-- This means that for instance, even though there are two ways to derive
-- @'Ord' a ':-' 'Eq' [a]@, the answers from these two paths _must_ by
-- construction be equal. This is a property that Haskell offers that is
-- pretty much unique in the space of languages with things they call \"type
-- classes\".
--
-- What are the two ways?
--
-- Well, we can go from @'Ord' a ':-' 'Eq' a@ via the
-- superclass relationship, and then from @'Eq' a ':-' 'Eq' [a]@ via the
-- instance, or we can go from @'Ord' a ':-' 'Ord' [a]@ via the instance
-- then from @'Ord' [a] ':-' 'Eq' [a]@ through the superclass relationship
-- and this diagram by definition must \"commute\".
--
-- Diagrammatically,
--
-- >                    Ord a
-- >                ins /     \ cls
-- >                   v       v
-- >             Ord [a]     Eq a
-- >                cls \     / ins
-- >                     v   v
-- >                    Eq [a]
--
-- This safety net ensures that pretty much anything you can write with this
-- library is sensible and can't break any assumptions on the behalf of
-- library authors.
newtype a :- b = Sub (a => Dict b)
  deriving Typeable

type role (:-) nominal nominal

-- TODO: _proper_ Data for @(p ':-' q)@ requires @(:-)@ to be cartesian _closed_.
--
-- This is admissable, but not present by default

-- constraint should be instance (Typeable p, Typeable q, p |- q) => Data (p :- q)
instance (Typeable p, Typeable q, p, q) => Data (p :- q) where
  gfoldl _ z (Sub Dict) = z (Sub Dict)
  toConstr _ = subConstr
  gunfold _ z c = case constrIndex c of
    1 -> z (Sub Dict)
    _ -> error "gunfold"
  dataTypeOf _ = subDataType

subConstr :: Constr
subConstr = mkConstr dictDataType "Sub" [] Prefix

subDataType :: DataType
subDataType = mkDataType "Data.Constraint.:-" [subConstr]

-- | Possible since GHC 7.8, when 'Category' was made polykinded.
instance Category (:-) where
  id  = refl
  (.) = trans

-- | Assumes 'IncoherentInstances' doesn't exist.
instance Eq (a :- b) where
  _ == _ = True

-- | Assumes 'IncoherentInstances' doesn't exist.
instance Ord (a :- b) where
  compare _ _ = EQ

instance Show (a :- b) where
  showsPrec d _ = showParen (d > 10) $ showString "Sub Dict"

instance a => NFData (a :- b) where
  rnf (Sub Dict) = ()

infixl 1 \\ -- required comment

-- | Given that @a :- b@, derive something that needs a context @b@, using the context @a@
(\\) :: a => (b => r) -> (a :- b) -> r
r \\ Sub Dict = r

--------------------------------------------------------------------------------
-- Constraints form a Category
--------------------------------------------------------------------------------

-- | Transitivity of entailment
--
-- If we view @(':-')@ as a Constraint-indexed category, then this is @('.')@
trans :: (b :- c) -> (a :- b) -> a :- c
trans f g = Sub $ Dict \\ f \\ g

-- | Reflexivity of entailment
--
-- If we view @(':-')@ as a Constraint-indexed category, then this is 'id'
refl :: a :- a
refl = Sub Dict

--------------------------------------------------------------------------------
-- (,) is a Bifunctor
--------------------------------------------------------------------------------

-- | due to the hack for the kind of @(,)@ in the current version of GHC we can't actually
-- make instances for @(,) :: Constraint -> Constraint -> Constraint@, but @(,)@ is a
-- bifunctor on the category of constraints. This lets us map over both sides.
(***) :: (a :- b) -> (c :- d) -> (a, c) :- (b, d)
f *** g = Sub $ Dict \\ f \\ g

--------------------------------------------------------------------------------
-- Constraints are Cartesian
--------------------------------------------------------------------------------

-- | Weakening a constraint product
--
-- The category of constraints is Cartesian. We can forget information.
weaken1 :: (a, b) :- a
weaken1 = Sub Dict

-- | Weakening a constraint product
--
-- The category of constraints is Cartesian. We can forget information.
weaken2 :: (a, b) :- b
weaken2 = Sub Dict

strengthen1 :: Dict b -> a :- c -> a :- (b,c)
strengthen1 d e = unmapDict (const d) &&& e

strengthen2 :: Dict b -> a :- c -> a :- (c,b)
strengthen2 d e = e &&& unmapDict (const d)

-- | Contracting a constraint / diagonal morphism
--
-- The category of constraints is Cartesian. We can reuse information.
contract :: a :- (a, a)
contract = Sub Dict

-- | Constraint product
--
-- > trans weaken1 (f &&& g) = f
-- > trans weaken2 (f &&& g) = g
(&&&) :: (a :- b) -> (a :- c) -> a :- (b, c)
f &&& g = Sub $ Dict \\ f \\ g

--------------------------------------------------------------------------------
-- Initial and terminal morphisms
--------------------------------------------------------------------------------

-- | Every constraint implies truth
--
-- These are the terminal arrows of the category, and @()@ is the terminal object.
--
-- Given any constraint there is a unique entailment of the @()@ constraint from that constraint.
top :: a :- ()
top = Sub Dict

-- | 'Any' inhabits every kind, including 'Constraint' but is uninhabited, making it impossible to define an instance.
class Exts.Any => Bottom where
  no :: a

-- |
-- This demonstrates the law of classical logic <http://en.wikipedia.org/wiki/Principle_of_explosion "ex falso quodlibet">
bottom :: Bottom :- a
bottom = Sub no

--------------------------------------------------------------------------------
-- Dict is fully faithful
--------------------------------------------------------------------------------

-- | Apply an entailment to a dictionary.
--
-- From a category theoretic perspective 'Dict' is a functor that maps from the category
-- of constraints (with arrows in ':-') to the category Hask of Haskell data types.
mapDict :: (a :- b) -> Dict a -> Dict b
mapDict p Dict = case p of Sub q -> q

-- |
-- This functor is fully faithful, which is to say that given any function you can write
-- @Dict a -> Dict b@ there also exists an entailment @a :- b@ in the category of constraints
-- that you can build.
unmapDict :: (Dict a -> Dict b) -> a :- b
unmapDict f = Sub (f Dict)

type role Dict nominal

--------------------------------------------------------------------------------
-- Reflection
--------------------------------------------------------------------------------

-- | Reify the relationship between a class and its superclass constraints as a class
--
-- Given a definition such as
--
-- @
-- class Foo a => Bar a
-- @
--
-- you can capture the relationship between 'Bar a' and its superclass 'Foo a' with
--
-- @
-- instance 'Class' (Foo a) (Bar a) where 'cls' = 'Sub' 'Dict'
-- @
--
-- Now the user can use 'cls :: Bar a :- Foo a'
class Class b h | h -> b where
  cls :: h :- b

infixr 9 :=>
-- | Reify the relationship between an instance head and its body as a class
--
-- Given a definition such as
--
-- @
-- instance Foo a => Foo [a]
-- @
--
-- you can capture the relationship between the instance head and its body with
--
-- @
-- instance Foo a ':=>' Foo [a] where 'ins' = 'Sub' 'Dict'
-- @
class b :=> h | h -> b where
  ins :: b :- h

-- Bootstrapping

instance Class () (Class b a) where cls = Sub Dict
instance Class () (b :=> a) where cls = Sub Dict

instance Class b a => () :=> Class b a where ins = Sub Dict
instance (b :=> a) => () :=> (b :=> a) where ins = Sub Dict

instance Class () () where cls = Sub Dict
instance () :=> () where ins = Sub Dict

-- Local, Prelude, Applicative, C.M.I and Data.Monoid instances

-- Eq
instance Class () (Eq a) where cls = Sub Dict
instance () :=> Eq () where ins = Sub Dict
instance () :=> Eq Int where ins = Sub Dict
instance () :=> Eq Bool where ins = Sub Dict
instance () :=> Eq Integer where ins = Sub Dict
instance () :=> Eq Float where ins = Sub Dict
instance () :=> Eq Double where ins = Sub Dict
instance Eq a :=> Eq [a] where ins = Sub Dict
instance Eq a :=> Eq (Maybe a) where ins = Sub Dict
instance Eq a :=> Eq (Complex a) where ins = Sub Dict
instance Eq a :=> Eq (Ratio a) where ins = Sub Dict
instance (Eq a, Eq b) :=> Eq (a, b) where ins = Sub Dict
instance (Eq a, Eq b) :=> Eq (Either a b) where ins = Sub Dict
instance () :=> Eq (Dict a) where ins = Sub Dict
instance () :=> Eq (a :- b) where ins = Sub Dict
instance () :=> Eq Word where ins = Sub Dict
instance Eq a :=> Eq (Identity a) where ins = Sub Dict
#if MIN_VERSION_base(4,8,0)
instance Eq a :=> Eq (Const a b) where ins = Sub Dict
instance () :=> Eq Natural where ins = Sub Dict
#endif

-- Ord
instance Class (Eq a) (Ord a) where cls = Sub Dict
instance () :=> Ord () where ins = Sub Dict
instance () :=> Ord Bool where ins = Sub Dict
instance () :=> Ord Int where ins = Sub Dict
instance ():=> Ord Integer where ins = Sub Dict
instance () :=> Ord Float where ins = Sub Dict
instance ():=> Ord Double where ins = Sub Dict
instance () :=> Ord Char where ins = Sub Dict
instance Ord a :=> Ord (Maybe a) where ins = Sub Dict
instance Ord a :=> Ord [a] where ins = Sub Dict
instance (Ord a, Ord b) :=> Ord (a, b) where ins = Sub Dict
instance (Ord a, Ord b) :=> Ord (Either a b) where ins = Sub Dict
instance Integral a :=> Ord (Ratio a) where ins = Sub Dict
instance () :=> Ord (Dict a) where ins = Sub Dict
instance () :=> Ord (a :- b) where ins = Sub Dict
instance () :=> Ord Word where ins = Sub Dict
instance Ord a :=> Ord (Identity a) where ins = Sub Dict
#if MIN_VERSION_base(4,8,0)
instance Ord a :=> Ord (Const a b) where ins = Sub Dict
instance () :=> Ord Natural where ins = Sub Dict
#endif

-- Show
instance Class () (Show a) where cls = Sub Dict
instance () :=> Show () where ins = Sub Dict
instance () :=> Show Bool where ins = Sub Dict
instance () :=> Show Ordering where ins = Sub Dict
instance () :=> Show Char where ins = Sub Dict
instance () :=> Show Int where ins = Sub Dict
instance Show a :=> Show (Complex a) where ins = Sub Dict
instance Show a :=> Show [a] where ins = Sub Dict
instance Show a :=> Show (Maybe a) where ins = Sub Dict
instance (Show a, Show b) :=> Show (a, b) where ins = Sub Dict
instance (Show a, Show b) :=> Show (Either a b) where ins = Sub Dict
instance (Integral a, Show a) :=> Show (Ratio a) where ins = Sub Dict
instance () :=> Show (Dict a) where ins = Sub Dict
instance () :=> Show (a :- b) where ins = Sub Dict
instance () :=> Show Word where ins = Sub Dict
instance Show a :=> Show (Identity a) where ins = Sub Dict
#if MIN_VERSION_base(4,8,0)
instance Show a :=> Show (Const a b) where ins = Sub Dict
instance () :=> Show Natural where ins = Sub Dict
#endif

-- Read
instance Class () (Read a) where cls = Sub Dict
instance () :=> Read () where ins = Sub Dict
instance () :=> Read Bool where ins = Sub Dict
instance () :=> Read Ordering where ins = Sub Dict
instance () :=> Read Char where ins = Sub Dict
instance () :=> Read Int where ins = Sub Dict
instance Read a :=> Read (Complex a) where ins = Sub Dict
instance Read a :=> Read [a] where ins = Sub Dict
instance Read a :=> Read (Maybe a) where ins = Sub Dict
instance (Read a, Read b) :=> Read (a, b) where ins = Sub Dict
instance (Read a, Read b) :=> Read (Either a b) where ins = Sub Dict
instance (Integral a, Read a) :=> Read (Ratio a) where ins = Sub Dict
instance () :=> Read Word where ins = Sub Dict
instance Read a :=> Read (Identity a) where ins = Sub Dict
#if MIN_VERSION_base(4,8,0)
instance Read a :=> Read (Const a b) where ins = Sub Dict
instance () :=> Read Natural where ins = Sub Dict
#endif

-- Enum
instance Class () (Enum a) where cls = Sub Dict
instance () :=> Enum () where ins = Sub Dict
instance () :=> Enum Bool where ins = Sub Dict
instance () :=> Enum Ordering where ins = Sub Dict
instance () :=> Enum Char where ins = Sub Dict
instance () :=> Enum Int where ins = Sub Dict
instance () :=> Enum Integer where ins = Sub Dict
instance () :=> Enum Float where ins = Sub Dict
instance () :=> Enum Double where ins = Sub Dict
instance Integral a :=> Enum (Ratio a) where ins = Sub Dict
instance () :=> Enum Word where ins = Sub Dict
#if MIN_VERSION_base(4,9,0)
instance Enum a :=> Enum (Identity a) where ins = Sub Dict
instance Enum a :=> Enum (Const a b) where ins = Sub Dict
#endif
#if MIN_VERSION_base(4,8,0)
instance () :=> Enum Natural where ins = Sub Dict
#endif

-- Bounded
instance Class () (Bounded a) where cls = Sub Dict
instance () :=> Bounded () where ins = Sub Dict
instance () :=> Bounded Ordering where ins = Sub Dict
instance () :=> Bounded Bool where ins = Sub Dict
instance () :=> Bounded Int where ins = Sub Dict
instance () :=> Bounded Char where ins = Sub Dict
instance (Bounded a, Bounded b) :=> Bounded (a,b) where ins = Sub Dict
instance () :=> Bounded Word where ins = Sub Dict
#if MIN_VERSION_base(4,9,0)
instance Bounded a :=> Bounded (Identity a) where ins = Sub Dict
instance Bounded a :=> Bounded (Const a b) where ins = Sub Dict
#endif

-- Num
instance Class () (Num a) where cls = Sub Dict
instance () :=> Num Int where ins = Sub Dict
instance () :=> Num Integer where ins = Sub Dict
instance () :=> Num Float where ins = Sub Dict
instance () :=> Num Double where ins = Sub Dict
instance RealFloat a :=> Num (Complex a) where ins = Sub Dict
instance Integral a :=> Num (Ratio a) where ins = Sub Dict
instance () :=> Num Word where ins = Sub Dict
#if MIN_VERSION_base(4,9,0)
instance Num a :=> Num (Identity a) where ins = Sub Dict
instance Num a :=> Num (Const a b) where ins = Sub Dict
#endif
#if MIN_VERSION_base(4,8,0)
instance () :=> Num Natural where ins = Sub Dict
#endif

-- Real
instance Class (Num a, Ord a) (Real a) where cls = Sub Dict
instance () :=> Real Int where ins = Sub Dict
instance () :=> Real Integer where ins = Sub Dict
instance () :=> Real Float where ins = Sub Dict
instance () :=> Real Double where ins = Sub Dict
instance Integral a :=> Real (Ratio a) where ins = Sub Dict
instance () :=> Real Word where ins = Sub Dict
#if MIN_VERSION_base(4,9,0)
instance Real a :=> Real (Identity a) where ins = Sub Dict
instance Real a :=> Real (Const a b) where ins = Sub Dict
#endif
#if MIN_VERSION_base(4,8,0)
instance () :=> Real Natural where ins = Sub Dict
#endif

-- Integral
instance Class (Real a, Enum a) (Integral a) where cls = Sub Dict
instance () :=> Integral Int where ins = Sub Dict
instance () :=> Integral Integer where ins = Sub Dict
instance () :=> Integral Word where ins = Sub Dict
#if MIN_VERSION_base(4,9,0)
instance Integral a :=> Integral (Identity a) where ins = Sub Dict
instance Integral a :=> Integral (Const a b) where ins = Sub Dict
#endif
#if MIN_VERSION_base(4,8,0)
instance () :=> Integral Natural where ins = Sub Dict
#endif

-- Bits
instance Class (Eq a) (Bits a) where cls = Sub Dict
instance () :=> Bits Bool where ins = Sub Dict
instance () :=> Bits Int where ins = Sub Dict
instance () :=> Bits Integer where ins = Sub Dict
instance () :=> Bits Word where ins = Sub Dict
#if MIN_VERSION_base(4,9,0)
instance Bits a :=> Bits (Identity a) where ins = Sub Dict
instance Bits a :=> Bits (Const a b) where ins = Sub Dict
#endif
#if MIN_VERSION_base(4,8,0)
instance () :=> Bits Natural where ins = Sub Dict
#endif

-- Fractional
instance Class (Num a) (Fractional a) where cls = Sub Dict
instance () :=> Fractional Float where ins = Sub Dict
instance () :=> Fractional Double where ins = Sub Dict
instance RealFloat a :=> Fractional (Complex a) where ins = Sub Dict
instance Integral a :=> Fractional (Ratio a) where ins = Sub Dict
#if MIN_VERSION_base(4,9,0)
instance Fractional a :=> Fractional (Identity a) where ins = Sub Dict
instance Fractional a :=> Fractional (Const a b) where ins = Sub Dict
#endif

-- Floating
instance Class (Fractional a) (Floating a) where cls = Sub Dict
instance () :=> Floating Float where ins = Sub Dict
instance () :=> Floating Double where ins = Sub Dict
instance RealFloat a :=> Floating (Complex a) where ins = Sub Dict
#if MIN_VERSION_base(4,9,0)
instance Floating a :=> Floating (Identity a) where ins = Sub Dict
instance Floating a :=> Floating (Const a b) where ins = Sub Dict
#endif

-- RealFrac
instance Class (Real a, Fractional a) (RealFrac a) where cls = Sub Dict
instance () :=> RealFrac Float where ins = Sub Dict
instance () :=> RealFrac Double where ins = Sub Dict
instance Integral a :=> RealFrac (Ratio a) where ins = Sub Dict
#if MIN_VERSION_base(4,9,0)
instance RealFrac a :=> RealFrac (Identity a) where ins = Sub Dict
instance RealFrac a :=> RealFrac (Const a b) where ins = Sub Dict
#endif

-- RealFloat
instance Class (RealFrac a, Floating a) (RealFloat a) where cls = Sub Dict
instance () :=> RealFloat Float where ins = Sub Dict
instance () :=> RealFloat Double where ins = Sub Dict
#if MIN_VERSION_base(4,9,0)
instance RealFloat a :=> RealFloat (Identity a) where ins = Sub Dict
instance RealFloat a :=> RealFloat (Const a b) where ins = Sub Dict
#endif

-- Semigroup
instance Class () (Semigroup a) where cls = Sub Dict
instance () :=> Semigroup () where ins = Sub Dict
instance () :=> Semigroup Ordering where ins = Sub Dict
instance () :=> Semigroup [a] where ins = Sub Dict
instance Semigroup a :=> Semigroup (Maybe a) where ins = Sub Dict
instance (Semigroup a, Semigroup b) :=> Semigroup (a, b) where ins = Sub Dict
instance Semigroup a :=> Semigroup (Const a b) where ins = Sub Dict
#if MIN_VERSION_base(4,9,0)
instance Semigroup a :=> Semigroup (Identity a) where ins = Sub Dict
#endif
#if MIN_VERSION_base(4,10,0)
instance Semigroup a :=> Semigroup (IO a) where ins = Sub Dict
#endif

-- Monoid
#if MIN_VERSION_base(4,11,0)
instance Class (Semigroup a) (Monoid a) where cls = Sub Dict
#else
instance Class () (Monoid a) where cls = Sub Dict
#endif
instance () :=> Monoid () where ins = Sub Dict
instance () :=> Monoid Ordering where ins = Sub Dict
instance () :=> Monoid [a] where ins = Sub Dict
instance Monoid a :=> Monoid (Maybe a) where ins = Sub Dict
instance (Monoid a, Monoid b) :=> Monoid (a, b) where ins = Sub Dict
instance Monoid a :=> Monoid (Const a b) where ins = Sub Dict
#if MIN_VERSION_base(4,9,0)
instance Monoid a :=> Monoid (Identity a) where ins = Sub Dict
instance Monoid a :=> Monoid (IO a) where ins = Sub Dict
#endif

-- Functor
instance Class () (Functor f) where cls = Sub Dict
instance () :=> Functor [] where ins = Sub Dict
instance () :=> Functor Maybe where ins = Sub Dict
instance () :=> Functor (Either a) where ins = Sub Dict
instance () :=> Functor ((->) a) where ins = Sub Dict
instance () :=> Functor ((,) a) where ins = Sub Dict
instance () :=> Functor IO where ins = Sub Dict
instance Monad m :=> Functor (WrappedMonad m) where ins = Sub Dict
instance () :=> Functor Identity where ins = Sub Dict
instance () :=> Functor (Const a) where ins = Sub Dict

-- Applicative
instance Class (Functor f) (Applicative f) where cls = Sub Dict
instance () :=> Applicative [] where ins = Sub Dict
instance () :=> Applicative Maybe where ins = Sub Dict
instance () :=> Applicative (Either a) where ins = Sub Dict
instance () :=> Applicative ((->)a) where ins = Sub Dict
instance () :=> Applicative IO where ins = Sub Dict
instance Monoid a :=> Applicative ((,)a) where ins = Sub Dict
instance Monoid a :=> Applicative (Const a) where ins = Sub Dict
instance Monad m :=> Applicative (WrappedMonad m) where ins = Sub Dict

-- Alternative
instance Class (Applicative f) (Alternative f) where cls = Sub Dict
instance () :=> Alternative [] where ins = Sub Dict
instance () :=> Alternative Maybe where ins = Sub Dict
instance MonadPlus m :=> Alternative (WrappedMonad m) where ins = Sub Dict

-- Monad
#if MIN_VERSION_base(4,8,0)
instance Class (Applicative f) (Monad f) where cls = Sub Dict
#else
instance Class () (Monad f) where cls = Sub Dict
#endif
instance () :=> Monad [] where ins = Sub Dict
instance () :=> Monad ((->) a) where ins = Sub Dict
instance () :=> Monad (Either a) where ins = Sub Dict
instance () :=> Monad IO where ins = Sub Dict
instance () :=> Monad Identity where ins = Sub Dict

-- MonadPlus
#if MIN_VERSION_base(4,8,0)
instance Class (Monad f, Alternative f) (MonadPlus f) where cls = Sub Dict
#else
instance Class (Monad f) (MonadPlus f) where cls = Sub Dict
#endif
instance () :=> MonadPlus [] where ins = Sub Dict
instance () :=> MonadPlus Maybe where ins = Sub Dict

--------------------------------------------------------------------------------
-- UndecidableInstances
--------------------------------------------------------------------------------

instance a :=> Enum (Dict a) where ins = Sub Dict
instance a => Enum (Dict a) where
  toEnum _ = Dict
  fromEnum Dict = 0

instance a :=> Bounded (Dict a) where ins = Sub Dict
instance a => Bounded (Dict a) where
  minBound = Dict
  maxBound = Dict

instance a :=> Read (Dict a) where ins = Sub Dict
deriving instance a => Read (Dict a)

instance () :=> Semigroup (Dict a) where ins = Sub Dict
instance Semigroup (Dict a) where
  Dict <> Dict = Dict

instance a :=> Monoid (Dict a) where ins = Sub Dict
instance a => Monoid (Dict a) where
#if !(MIN_VERSION_base(4,11,0))
  mappend = (<>)
#endif
  mempty = Dict

max_stars_repo_name: phadej/constraints

max_issues_repo_name: phadej/constraints

max_forks_repo_name: phadej/constraints


text:
{-# LANGUAGE ConstraintKinds, DataKinds, FlexibleContexts        #-}
{-# LANGUAGE MultiParamTypeClasses, NoImplicitPrelude, PolyKinds #-}
{-# LANGUAGE QuasiQuotes, TemplateHaskell                        #-}
{-# OPTIONS_GHC -fno-warn-type-defaults #-}
module Main (module Algebra.Algorithms.Groebner, module Algebra.Prelude
            , module Main
            ) where
import           Algebra.Algorithms.Groebner
import           Algebra.Algorithms.ZeroDim
import           Algebra.Ring.Ideal
import           Algebra.Ring.Polynomial.Quotient
import           Algebra.Prelude
import           Control.Monad.Random             hiding (fromList)
import           Data.Complex
import           Data.Convertible
import           Data.List                        (find, nub, partition, sortBy)
import qualified Data.Matrix                      as M
import qualified Data.Sized               as SV
import qualified Data.Vector                      as V
import           Debug.Trace
import qualified Numeric.LinearAlgebra            as LA
import qualified Prelude                          as P

tr :: Show a => a -> a
tr a = trace (show a) a

x, y, z :: Polynomial (Fraction Integer) 3
[x, y, z] = vars

seed :: Polynomial (Fraction Integer) 3
seed = -412742019532366985 * x -7641395389638504101 * y + 4362835172800530323 * z

seedMat :: LA.Matrix Double
seedMat = LA.fromLists $ map (map toDouble) $ reifyQuotient eqn02 $ \pxy -> matrixRep (modIdeal' pxy seed)

toDouble :: Fractional a => Fraction Integer -> a
toDouble rat = fromIntegral (numerator rat) P./ fromIntegral (denominator rat)

fromRight :: Either t t1 -> t1
fromRight (Right a) = a
fromRight _ = error "fromRight"

printLvl :: Show a => Int -> a -> IO ()
printLvl lvl = putStrLn . unlines . map (replicate lvl '\t' ++) . lines . show

eqn01 :: Ideal (Polynomial (Fraction Integer) 3)
eqn01 = toIdeal [x^2 - 2*x*z + 5, x*y^2+y*z+1, 3*y^2 - 8*x*z]

eqn02 :: Ideal (Polynomial (Fraction Integer) 3)
eqn02 =
  toIdeal [x^2 + 2*y^2 - y - 2*z
          ,x^2 - 8*y^2 + 10*z - 1
          ,x^2 - 7*y*z
          ]

eqn03 :: Ideal (Polynomial (Fraction Integer) 3)
eqn03 = toIdeal [x^2 + y^2 + z^2 - 2*x
                ,x^3 - y*z - x
                ,x - y + 2*z
                ]

jdeal :: Ideal (Polynomial (Fraction Integer) 3)
jdeal = toIdeal [x*y + z - x*z, x^2 - z, 2*x^3 - x^2 * y * z - 1]

vs :: [V.Vector (Fraction Integer)]
vs = reifyQuotient eqn03 $ \pxy -> map (vectorRep . modIdeal' pxy) [var 0 ^ i | i <- [0..6::Natural]]

mat :: M.Matrix (Fraction Integer)
mat = fromCols $ take 4 vs

fromCols :: [V.Vector a] -> M.Matrix a
fromCols = foldr1 (M.<|>) . map M.colVector

findUnivar :: (CoeffRing r, IsMonomialOrder n ord, KnownNat n)
           => OrderedPolynomial r ord n -> Maybe (Ordinal n)
findUnivar poly =
  let os = enumOrdinal (sArity' poly)
      ms = map snd $ getTerms poly
  in find (\a -> all (`isPowerOf` (leadingMonomial (var a `asTypeOf` poly))) ms) os

toCoeffList :: (CoeffRing r,  KnownNat n, IsMonomialOrder n ord) => Ordinal n -> OrderedPolynomial r ord n -> [r]
toCoeffList on f =
  let v = var on  `asTypeOf` f
  in [ coeff (leadingMonomial $ v ^ i) f | i <- [0.. fromIntegral (totalDegree' f)]]

showSols :: (KnownNat n, IsMonomialOrder n order, Convertible a Double)
         => Double -> Ideal (OrderedPolynomial a order n) -> [Sized n1 (Complex Double)] -> IO ()
showSols err eqn sols = do
  let (rs, is) = partition (all ((<err).P.abs.imagPart)) $ map SV.toList sols
      subs a b c = generators $
                   mapIdeal (magnitude . substWith (*) (SV.unsafeFromList' [a, b, c]) . mapCoeff toComplex)
                            eqn
      showCase [a,b,c] = print (a, b, c) >> putStr "\terror: ">> print (maximum $ subs a b c)
  putStrLn $ "- " ++ show (length rs) ++ " real solution(s):"
  mapM_ showCase $ sortBy (comparing $ map magnitude) rs
  putStrLn $ "- " ++ show (length is) ++ " imaginary solution(s):"
  mapM_  showCase $ sortBy (comparing $ map magnitude) is
  let errs = concatMap (\ [a,b,c] -> subs a b c) $ rs ++ is
  putStrLn $ "- maximum error: " ++ show (maximum errs)
  putStrLn $ "- minimum error: " ++ show (minimum errs)
  putStrLn $ "- average error: " ++ show (sum errs P./ fromIntegral (length errs))

main :: IO ()
main = do
  putStrLn "---- solving equation system"
  let err = 1e-10
  putStrLn "< naive method"
  showSols err eqn01 $ solve' err eqn01
  putStrLn "\n< randomized method"
  showSols err eqn01 =<< evalRandIO (solveM eqn01)
  putStrLn "\n< companion characteristics"
  showSols err eqn01 $ solveViaCompanion err eqn01
  putStrLn "\n< univariate spanning"
  showSols err eqn01 $ solve' err eqn01

  putStrLn "\n\n---- exercise 8"
  putStrLn "< Solving 1-6"
  putStrLn "< Naive Method: "
  showSols err eqn02 $ nub $ solve' err eqn02
  putStrLn "\n< new method"
  showSols err eqn02 =<< evalRandIO (solveM eqn02)

  putStrLn "\n< Solving 1-7"
  putStrLn "< Naive Method: "
  showSols err eqn03 $ nub $ solve' err eqn03
  putStrLn "\n< new method"
  showSols err eqn03 =<< evalRandIO (solveM eqn03)
  putStrLn "\n\n---- FGLM Algorithm"
  print $ fglm jdeal
  print $ calcGroebnerBasisWith Lex jdeal
  print $ univPoly 0 jdeal
  print $ univPoly 1 jdeal
  print $ univPoly 2 jdeal
  return ()

substIdeal :: IsMonomialOrder 3 order
           => [OrderedPolynomial (Fraction Integer) Grevlex 3]
           -> Ideal (OrderedPolynomial (Fraction Integer) order 3)
           -> Ideal (OrderedPolynomial (Fraction Integer) Grevlex 3)
substIdeal = mapIdeal . substWith (.*.) . SV.unsafeFromList'

toComplex :: Convertible r Double => r -> Complex Double
toComplex = (:+ 0) . convert

max_stars_repo_name: hangingman/computational-algebra

max_issues_repo_name: hangingman/computational-algebra

max_forks_repo_name: hangingman/computational-algebra


text:
-- Modulo destinado a leitura e escrita de arquivos.

module InputOutput
(getTraining,
 getTest,
 readIn,
 writeIn,
 isEmptyFile,
 definitiveAnswer,
 printEpoch
 ) where

import Data.List
import Data.Char
import Data.Maybe
import System.Directory
import System.Random
import Types
import Text.Printf
import Numeric.LinearAlgebra.HMatrix

type Image = [Double]
type Sample = (Int, Image)

-- Transforma uma String em uma lista
strToArr::String->[Double]
strToArr line = [(read x::Double)| x<-(words line)]

-- Transforma uma String em uma matrix
strToMatrix::String->[[Double]]
strToMatrix mt = [(strToArr x)| x<-(lines mt)]

-- Transforma uma lista em uma String
arrToStr::[Double]->String
arrToStr [] = ""
arrToStr (h:t) = (show h::String) ++ " " ++ (arrToStr t)

-- Transforma uma matriz em uma String
matrixToStr::[[Double]]->String
matrixToStr [] = ""
matrixToStr (h:t) = (arrToStr h) ++ "\n" ++ (matrixToStr t)

-- Dado o seu tamanho, gera uma lista
-- com valores randomicos
listRandom::Int->IO [Double]
listRandom 0 = return []
listRandom sz = do
    head <- randomIO::IO Double
    tail <- listRandom (sz-1)
    return (head:tail)

-- Dado o seu tamanho, gera uma matrix
-- com valores randomicos
matrixRandom::Int->Int->IO [[Double]]
matrixRandom 0 y = return []
matrixRandom x y = do
    head <- listRandom y
    tail <- matrixRandom (x-1) y
    return (head:tail)

-- Dado o seu caminho, le o arquivo
-- e o transforma em uma lista
getArchiveL::String->Int->IO [Double]
getArchiveL path sz = do
    emp <- isEmptyFile path
    if not emp then
        do
            archive <- (readFile path)
            return $ strToArr archive
    else
        do
            ret <- listRandom sz
            return ret

-- Dado o seu caminho, le o arquivo
-- e o transforma em uma matriz
getArchiveM::String->Int->Int->IO [[Double]]
getArchiveM path x y = do
    emp <- isEmptyFile path
    if not emp then
        do
            archive <- readFile path
            return $ strToMatrix archive
    else
        do
            ret <- matrixRandom x y
            return ret

-- Le os vetores e matrizes presentes nos arquivos para montar e retornar um Data
readIn::IO Data
readIn = do
    wH <- getArchiveM "Data/Weight_hidden.txt" 30 784
    bH <- getArchiveL "Data/Biases_hidden.txt" 30
    aH <- getArchiveL "Data/Activation_hidden.txt" 30
    zH <- getArchiveL "Data/Zeta_hidden.txt" 30

    wO <- getArchiveM "Data/Weight_Output.txt" 10 30
    bO <- getArchiveL "Data/Biases_Output.txt" 10
    aO <- getArchiveL "Data/Activation_Output.txt" 10
    zO <- getArchiveL "Data/Zeta_Output.txt" 10

    return $ (Data (fromLists wH) (fromList bH) (fromList aH) (fromList zH) (fromLists wO) (fromList bO) (fromList aO) (fromList zO))

-- Escreve o Data nos arquivos
writeIn::Data->IO()
writeIn elem = do
    writeFile "NedeReural_Haskell/Data/Weight_hidden.txt" (matrixToStr (toLists $ wHidden elem))
    writeFile "NedeReural_Haskell/Data/Biases_hidden.txt" (arrToStr (toList $ bHidden elem))
    writeFile "NedeReural_Haskell/Data/Activation_hidden.txt" (arrToStr (toList $ aHidden elem))
    writeFile "NedeReural_Haskell/Data/Zeta_hidden.txt" (arrToStr (toList $ zetaHidden elem))

    writeFile "NedeReural_Haskell/Data/Weight_Output.txt" (matrixToStr (toLists $ wOutput elem))
    writeFile "NedeReural_Haskell/Data/Biases_Output.txt" (arrToStr (toList $ bOutput elem))
    writeFile "NedeReural_Haskell/Data/Activation_Output.txt" (arrToStr (toList $ aOutput elem))
    writeFile "NedeReural_Haskell/Data/Zeta_Output.txt" (arrToStr (toList $ zetaOutput elem))

-- Verifica se o arquivo esta vazio
isEmptyFile::String->IO Bool
isEmptyFile path = do
    elem <- readFile path

    if (length elem) == 0 then return True
    else return False

-- Pega o numero que define a resposta do teste/treino
pickNumber::FilePath->Int
pickNumber "" = -1
pickNumber (h:t)
                | h == '-' = digitToInt $ t!!0
                | otherwise = (pickNumber t)

-- Constroi uma imagem
makeImage::FilePath->IO Image
makeImage path = do
    elem <- readFile path
    let ret = (strToArr elem)
    return ret

-- Constroi um Sample
makeSample::FilePath->String->IO Sample
makeSample archive path = do
    let number = (pickNumber archive)
    mt <- (makeImage $ path++archive)
    return (number, mt)

-- Cria uma lista de Sample dado os arquivos e o local
listSample::[FilePath]->String->IO [Sample]
listSample [] path = return []
listSample (h:t) path = do
                head <- makeSample h path
                tail <- listSample t path
                return $ head:tail

-- Retorna o TestSet da rede
getTest::IO [Sample]
getTest = do
    let caminho = "NedeReural_Haskell/Tests" -- Diretorio de testes
    elems <- (getDirectoryContents caminho) -- Colocar diretorio
    ret <- listSample (filter (/= ".") $ filter (/= "..") elems) caminho
    return ret

-- Retorna o TraninigSet da rede
getTraining::IO [Sample]
getTraining = do
    let caminho = "NedeReural_Haskell/Training" -- Diretorio de treino
    elems <- (getDirectoryContents caminho) -- Colocar diretorio
    ret <- listSample (filter (/= ".") $ filter (/= "..") elems) caminho
    return ret

-- Imprime na tela a taxa de acerto de uma epoca de teste
printEpoch :: Int -> Int -> Int -> IO()
printEpoch epochIndex correctCnt total = putStrLn $ "EPOCH #" ++ (show epochIndex) ++ " - " ++ (show correctCnt) ++ " / " ++ (show total)

-- Recebe a lista de ativacao como parametro
-- Retorna uma tupla: o melhor sigmoid e seu respectivo indice
getBestSigmoid :: [Double] -> (Double, Int)
getBestSigmoid activationValues = let first = maximum (activationValues)
                                      second = fromJust $ (elemIndex (maximum (activationValues)) (activationValues))
                                  in (first, second)

-- Recebe a tupla (sigmoid, indice)
-- Retorna a formatacao em string que representa
-- o numero e a probabilidade de ser esse numero
toPercentage :: (Double, Int) -> [Double] -> String
toPercentage sigmoid activationValues =
    "[" ++ (printf "%d" (snd sigmoid + 1)) ++ " - " ++ (printf "%.2f" ((fst sigmoid) * 100 / sum (activationValues))) ++ "%]"

-- Computa a resposta definitiva a partir do
-- array de valores de ativacao
definitiveAnswer :: Data -> IO String
definitiveAnswer activationValues = return ("Resposta definitiva: " ++ toPercentage (getBestSigmoid (toList $ aOutput activationValues)) (toList $ aOutput activationValues))


text:
module CNC.IntegrationTests where

import CNC.Declarative
import CNC.HCode

import Data.Complex
import Control.Monad

flake_side = do
  [p1,p2,p3,p4,p5] <- declarePoints 5
  xSize p1 p2 1
  len p2 p3 1
  len p3 p4 1
  xAngle p2 p3 (pi / 3)
  xAngle p3 p4 (- pi / 3)
  xSize p4 p5 1

renderPolygon :: Int -> Path -> HCode ()
renderPolygon n path = render_side 0 path n
  where
    render_side p0 path 0 = return ()
    render_side p0 path k = do
      let d = deltaPath path
      renderPath (posToExpr p0) path
      render_side (p0 + d) (rotate (- 2*pi / fromIntegral n) path) (k-1)

test1 = do putStrLn "a snowflake"
           print $ (figure flake_side)
           putHCode $ renderPolygon 3 (figure flake_side)

integrationTests = do
  putStrLn "Integration tests"
  test1

max_stars_repo_name: akamaus/gcodec

max_issues_repo_name: akamaus/gcodec

max_forks_repo_name: akamaus/gcodec