text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Data.Eigen.SparseMatrix (
    -- * SparseMatrix type
    SparseMatrix(..),
    SparseMatrixXf,
    SparseMatrixXd,
    SparseMatrixXcf,
    SparseMatrixXcd,
    -- * Matrix internal data
    values,
    innerIndices,
    outerStarts,
    innerNNZs,
    -- * Accessing matrix data
    cols,
    rows,
    coeff,
    (!),
    -- * Matrix conversions
    fromList,
    toList,
    fromVector,
    toVector,
    fromDenseList,
    toDenseList,
    fromMatrix,
    toMatrix,
    -- * Matrix properties
    norm,
    squaredNorm,
    blueNorm,
    block,
    nonZeros,
    innerSize,
    outerSize,
    -- * Basic matrix algebra
    add,
    sub,
    mul,
    -- * Matrix transformations
    pruned,
    scale,
    transpose,
    adjoint,
    -- * Matrix representation
    compress,
    uncompress,
    compressed,
    -- * Matrix serialization
    encode,
    decode,
    -- * Mutable matricies
    thaw,
    freeze,
    unsafeThaw,
    unsafeFreeze,
) where

import qualified Prelude as P
import Prelude hiding (map)
import qualified Data.List as L
import Data.Complex
import Data.Binary hiding (encode, decode)
import qualified Data.Binary as B
import Foreign.C.Types
import Foreign.C.String
import Foreign.Storable
import Foreign.Ptr
import Foreign.ForeignPtr
import Foreign.Marshal.Alloc
import Control.Monad
#if __GLASGOW_HASKELL__ >= 710
#else
import Control.Applicative
#endif
import qualified Data.Eigen.Matrix as M
import qualified Data.Eigen.Matrix.Mutable as MM
import qualified Data.Eigen.SparseMatrix.Mutable as SMM
import qualified Foreign.Concurrent as FC
import qualified Data.Eigen.Internal as I
import qualified Data.Vector.Storable as VS
import qualified Data.Vector.Storable.Mutable as VSM
import qualified Data.ByteString.Lazy as BSL

{-| 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>.
-}

data SparseMatrix a b where
    SparseMatrix :: I.Elem a b => !(ForeignPtr (I.CSparseMatrix a b)) -> SparseMatrix a b

-- | Alias for single precision sparse matrix
type SparseMatrixXf = SparseMatrix Float CFloat
-- | Alias for double precision sparse matrix
type SparseMatrixXd = SparseMatrix Double CDouble
-- | Alias for single previsiom sparse matrix of complex numbers
type SparseMatrixXcf = SparseMatrix (Complex Float) (I.CComplex CFloat)
-- | Alias for double prevision sparse matrix of complex numbers
type SparseMatrixXcd = SparseMatrix (Complex Double) (I.CComplex CDouble)

-- | Pretty prints the sparse matrix
instance (I.Elem a b, Show a) => Show (SparseMatrix a b) where
    show m = concat [
        "SparseMatrix ", show (rows m), "x", show (cols m),
        "\n", L.intercalate "\n" $ P.map (L.intercalate "\t" . P.map show) $ toDenseList m, "\n"]

-- | Basic sparse matrix math exposed through Num instance: @(*)@, @(+)@, @(-)@, `fromInteger`, `signum`, `abs`, `negate`
instance I.Elem a b => Num (SparseMatrix a b) where
    (*) = mul
    (+) = add
    (-) = sub
    fromInteger x = fromList 1 1 [(0,0,fromInteger x)]
    signum = _map signum
    abs = _map abs
    negate = _map negate

instance I.Elem a b => Binary (SparseMatrix a b) where
    put m = do
        put $ I.magicCode (undefined :: b)
        put $ rows m
        put $ cols m
        put $ toVector m

    get = do
        get >>= (`when` fail "wrong matrix type") . (/= I.magicCode (undefined :: b))
        fromVector <$> get <*> get <*> get

-- | Encode the sparse matrix as a lazy byte string
encode :: I.Elem a b => SparseMatrix a b -> BSL.ByteString
encode = B.encode

-- | Decode sparse matrix from the lazy byte string
decode :: I.Elem a b => BSL.ByteString -> SparseMatrix a b
decode = B.decode

-- | Stores the coefficient values of the non-zeros.
values :: I.Elem a b => SparseMatrix a b -> VS.Vector b
values = _getvec I.sparse_values

-- | Stores the row (resp. column) indices of the non-zeros.
innerIndices :: I.Elem a b => SparseMatrix a b -> VS.Vector CInt
innerIndices = _getvec I.sparse_innerIndices

-- | Stores for each column (resp. row) the index of the first non-zero in the previous two arrays.
outerStarts :: I.Elem a b => SparseMatrix a b -> VS.Vector CInt
outerStarts = _getvec I.sparse_outerStarts

-- | 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
innerNNZs :: I.Elem a b => SparseMatrix a b -> Maybe (VS.Vector CInt)
innerNNZs m
    | compressed m = Nothing
    | otherwise = Just $ _getvec I.sparse_innerNNZs m

-- | Number of rows for the sparse matrix
rows :: I.Elem a b => SparseMatrix a b -> Int
rows = _unop I.sparse_rows (return . I.cast)

-- | Number of columns for the sparse matrix
cols :: I.Elem a b => SparseMatrix a b -> Int
cols = _unop I.sparse_cols (return . I.cast)

-- | Matrix coefficient at given row and col
coeff :: I.Elem a b => Int -> Int -> SparseMatrix a b -> a
coeff row col (SparseMatrix fp) = I.performIO $ withForeignPtr fp $ \p -> alloca $ \pq -> do
    I.call $ I.sparse_coeff p (I.cast row) (I.cast col) pq
    I.cast <$> peek pq

-- | Matrix coefficient at given row and col
(!) :: I.Elem a b => SparseMatrix a b -> (Int, Int) -> 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 :: I.Elem a b => SparseMatrix a b -> a
norm = _unop I.sparse_norm (return . I.cast)

-- | 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 :: I.Elem a b => SparseMatrix a b -> a
squaredNorm = _unop I.sparse_squaredNorm (return . I.cast)

-- | 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 :: I.Elem a b => SparseMatrix a b -> a
blueNorm = _unop I.sparse_blueNorm (return . I.cast)

-- | Extract rectangular block from sparse matrix defined by startRow startCol blockRows blockCols
block :: I.Elem a b => Int -> Int -> Int -> Int -> SparseMatrix a b -> SparseMatrix a b
block row col rows cols = _unop (\p pq -> I.sparse_block p (I.cast row) (I.cast col) (I.cast rows) (I.cast cols) pq) _mk

-- | Number of non-zeros elements in the sparse matrix
nonZeros :: I.Elem a b => SparseMatrix a b -> Int
nonZeros = _unop I.sparse_nonZeros (return . I.cast)

-- | The matrix in the compressed format
compress :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b
compress = _unop I.sparse_makeCompressed _mk

-- | The matrix in the uncompressed mode
uncompress :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b
uncompress = _unop I.sparse_uncompress _mk

-- | Is this in compressed form?
compressed :: I.Elem a b => SparseMatrix a b -> Bool
compressed = _unop I.sparse_isCompressed (return . (/=0))

-- | Minor dimension with respect to the storage order
innerSize :: I.Elem a b => SparseMatrix a b -> Int
innerSize = _unop I.sparse_innerSize (return . I.cast)

-- | Major dimension with respect to the storage order
outerSize :: I.Elem a b => SparseMatrix a b -> Int
outerSize = _unop I.sparse_outerSize (return . I.cast)

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

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

-- | Transpose of the sparse matrix
transpose :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b
transpose = _unop I.sparse_transpose _mk

-- | Adjoint of the sparse matrix
adjoint :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b
adjoint = _unop I.sparse_adjoint _mk

-- | Adding two sparse matrices by adding the corresponding entries together. You can use @(+)@ function as well.
add :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b -> SparseMatrix a b
add = _binop I.sparse_add _mk

-- | Subtracting two sparse matrices by subtracting the corresponding entries together. You can use @(-)@ function as well.
sub :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b -> SparseMatrix a b
sub = _binop I.sparse_sub _mk

-- | Matrix multiplication. You can use @(*)@ function as well.
mul :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b -> SparseMatrix a b
mul = _binop I.sparse_mul _mk

-- | Construct sparse matrix of given size from the list of triplets (row, col, val)
fromList :: I.Elem a b => Int -> Int -> [(Int, Int, a)] -> SparseMatrix a b
fromList rows cols = fromVector rows cols . VS.fromList . P.map I.cast

-- | Construct sparse matrix of given size from the storable vector of triplets (row, col, val)
fromVector :: I.Elem a b => Int -> Int -> VS.Vector (I.CTriplet b) -> SparseMatrix a b
fromVector rows cols tris = I.performIO $ VS.unsafeWith tris $ \p -> alloca $ \pq -> do
    I.call $ I.sparse_fromList (I.cast rows) (I.cast cols) p (I.cast $ VS.length tris) pq
    peek pq >>= _mk

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

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

-- | Construct sparse matrix of two-dimensional list of values. Matrix dimensions will be detected automatically. Zero values will be compressed.
fromDenseList :: (I.Elem a b, Eq a) => [[a]] -> SparseMatrix a b
fromDenseList list = fromList rows cols $ do
    (row, vals) <- zip [0..] list
    (col, val) <- zip [0..] vals
    guard $ val /= 0
    return (row, col, val)
    where
        rows = length list
        cols = L.foldl' max 0 $ P.map length list

-- | Convert sparse matrix to (rows X cols) dense list of values
toDenseList :: I.Elem a b => SparseMatrix a b -> [[a]]
toDenseList m = [[coeff row col m | col <- [0 .. cols m - 1]] | row <- [0 .. rows m - 1]]

-- | Construct sparse matrix from dense matrix. Zero elements will be compressed
fromMatrix :: I.Elem a b => M.Matrix a b -> SparseMatrix a b
fromMatrix m1 = I.performIO $ alloca $ \pm0 ->
    M.unsafeWith m1 $ \vals rows cols -> do
        I.call $ I.sparse_fromMatrix vals rows cols pm0
        peek pm0 >>= _mk

-- | Construct dense matrix from sparse matrix
toMatrix :: I.Elem a b => SparseMatrix a b -> M.Matrix a b
toMatrix m1@(SparseMatrix fp) = I.performIO $ do
    m0 <- MM.new (rows m1) (cols m1)
    MM.unsafeWith m0 $ \vals rows cols ->
        withForeignPtr fp $ \pm1 ->
            I.call $ I.sparse_toMatrix pm1 vals rows cols
    M.unsafeFreeze m0

-- | Yield an immutable copy of the mutable matrix
freeze :: I.Elem a b => SMM.IOSparseMatrix a b -> IO (SparseMatrix a b)
freeze (SMM.IOSparseMatrix fp) = SparseMatrix <$> _clone fp

-- | Yield a mutable copy of the immutable matrix
thaw :: I.Elem a b => SparseMatrix a b -> IO (SMM.IOSparseMatrix a b)
thaw (SparseMatrix fp) = SMM.IOSparseMatrix <$> _clone fp

-- | Unsafe convert a mutable matrix to an immutable one without copying. The mutable matrix may not be used after this operation.
unsafeFreeze :: I.Elem a b => SMM.IOSparseMatrix a b -> IO (SparseMatrix a b)
unsafeFreeze (SMM.IOSparseMatrix 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 :: I.Elem a b => SparseMatrix a b -> IO (SMM.IOSparseMatrix a b)
unsafeThaw (SparseMatrix fp) = return $! SMM.IOSparseMatrix fp

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

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

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

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

_map :: I.Elem a b => (a -> a) -> SparseMatrix a b -> SparseMatrix a b
_map f m = fromVector (rows m) (cols m) . VS.map g . toVector $ m where
    g (I.CTriplet r c v) = I.CTriplet r c $ I.cast $ f $ I.cast v

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

max_stars_repo_name: osidorkin/haskell-eigen

max_issues_repo_name: osidorkin/haskell-eigen

max_forks_repo_name: osidorkin/haskell-eigen


text:
{-# LANGUAGE NoImplicitPrelude, TypeSynonymInstances, RankNTypes, StandaloneDeriving,
ScopedTypeVariables, TypeApplications, GADTs, TypeOperators, DeriveTraversable,
FlexibleInstances, AllowAmbiguousTypes, UndecidableInstances, TypeFamilies, LambdaCase #-}
module Lib
    where
import Numeric.LinearAlgebra
import Control.Category
import Prelude hiding ((.), id)
import Control.Arrow ((***))
someFunc :: IO ()
someFunc = putStrLn "someFunc"

{-
type V a = a -> Double
type D a = a -> Double
D (V a)
type V a  = ContT () Q a -- ContT b Q a = (a ->Q b) -> Q b. Has some way of dealing with a index. Some kind of
SelT -- takes a diag, unsummed?
-- (a -> Q b) -> Q a

newtype Dag a = Dag (a -> Q ())

class Dagger a where
    type Dual :: * -> * -- *?
    dag :: a -> Dual a
instance Dagger (D (Vec a)) where
    type Dual = Vec
    dag (D v) = join $ ContT ()
-- dagger needs to embed a search procedure?

instance Dagger (D (Vec Bool)) where
    type Dual = Vec Bool
    dag (D v) = join $ ContT ()
instance Dagger a b where
    dag :: a -> b
    dag' :: b -> a
instance Dagger k where
    type family Dual :: * -> *
    dag :: k a (Dual a) -- dag'?

-- also a kind of dot product
dag1 :: ((a -> Q ()) -> Q ()) -> (((a -> Q ()) -> Q ()) -> Q ()
dag1 f =

-- it's a kind of dot product.
lift :: f Double -> (a -> Double) -> Double
lift w f = fold (+) $ (mapWithIndex \k v -> (f k) * v) w

-- tensor is a_i b_i, mutiply out without summing
tensor :: f a ->  f a  -> f a
tensor = intersectWith (*)

type V' a = (a -> Q ()) -> Q ()
tensor :: V' a -> V' a -> V' a

dag :: (V' Bool -> Q ()) -> V' Bool
dag f = \g -> True False

dag :: (Bool -> Q ()) -> V' Bool
dag :: (Bool -> Q ()) -> (Bool -> Q ()) -> Q ()
dag f g = (f True) (g True) + (f False) (g False)

dag :: (Bool -> Q ()) -> Q Bool
dag f = (f True) >>= true + (f False) >>= false

dag :: Q Bool -> (Bool -> Q ())
dag xs b = lookup b xs
dag :: Eq a => Q a -> (a -> Q ()) -- maybe Ord
dag xs b = lookup b xs

-- If we mixed Hedges select sat solver with the Q monad?
-- (Bool -> Q Bool) -> Q Bool ? the unsummed trace? diag? Hedgest hing would come out weirder.

true :: () -> Q Bool
true _ = pure True
false :: () -> Q Bool
false _ = pure False

-- the seelect function.
-- (b -> m ()) -> m b

-}

type M = Matrix Double

{-

data BMatrix a b = BMatrix M M M M | Id -- A B C D

instance Category BMatrix where
    id = Id
    Id . x = x
    x . Id = x
    (BMatrix a b c d) . (BMatrix e f g h) = BMatrix <\> where
            q = d + e
            qlu = lupacked q
            qinv = lusolve qlu
            a' = a - (b <> (qinv c))
            b' =  b <> (qinv f)
            c' = g <> (qinv c)
            d' = h - g <> (qinv f)

par :: BMatrix a b -> BMatrix c d -> BMatrix (a,c) (b,d)
par (Matrix a b c d) (BMatrix e f g h) = BMatrix (diagBlock [a,e]) (diagBlock [b,f]) (diagBlock [c,g])  (diagBlock [d,h])
dup = BMatrix a (a,a)
dup = BMatrix (diagBlock [ident,ident]) (ident ||| ident) (ident === ident) (diagBlock [ident,ident])
id = BMatrix ident ident ident ident
avg =
-}
data QuadFun a b = QuadFun {mat :: M, vec :: Vector Double, c :: Double } deriving Show

--data QuadFun a b where
--   QuadFun :: (Bounded a, Bounded b, Enum a, Enum b) => M -> Vector Double -> Double -> QuadFun a b }

-- This is kron product.
instance (Enum a, Enum b, Bounded a) => Enum (a,b) where
  toEnum x =  let (a, b) = quotRem x (fromEnum (maxBound :: a)) in
                (toEnum a   , toEnum b)
  fromEnum (a,b) =  (fromEnum (maxBound :: a)) * (fromEnum a) + (fromEnum b)
  {-
instance (Enum a, Enum b, Bounded a) => Enum (Either a b) where
    toEnum x | x > (fromEnum (maxBound :: a)) = Left (toEnum x)
             | otherwise = Right $ toEnum (x - 1 - (fromEnum (maxBound :: a)))
    fromEnum (Left a) =  fromEnum a
    fromEnum (Right b) =  (fromEnum (maxBound :: a)) + 1 + (fromEnum b)
-}
-- deriving instance (Enum a, Enum b) => Enum (a,b)

count :: forall a. (Enum a, Bounded a) => Int
count = (fromEnum (maxBound @a)) - (fromEnum (minBound @a)) + 1

id ::forall a. (Enum a, Bounded a) => QuadFun a a
id = let s = 2 * (count @ a) in QuadFun (ident s) (konst 0 s) 0

--dup :: QuadFun a (a,a)
--dup  = QuadFun (ident

class (Bounded a, Enum a) => BEnum a where
instance (Bounded a, Enum a) => BEnum a where

t7 :: QuadFun () ()
t7 = QuadFun (konst 4 (2,2)) (konst 7 2) 3

t8 = (mat t7) <\> (vec t7)
t4 = par t7 t7
t9 = (mat t4) <\> (vec t4)

compose :: forall a b c. BEnum b => QuadFun b c -> QuadFun a b -> QuadFun a c
compose (QuadFun m' w' c') (QuadFun m w c) = QuadFun m'' w'' c'' where
   n = count @b
   k = rows m -- assuming square
   l = rows m'
   i = ident n
   q = konst 0 (k-n,n)  === i  -- I underneath zeros
   q' = -i === konst 0 (l-n,n)   -- -i above zeros
   m'' = fromBlocks [[m,    q, 0],
                     [tr q ,0, tr q'],
                     [0,    q', m']]
   w'' = vjoin [w, konst 0 n, w']
   c'' = c + c'

identOut :: forall a b. BEnum b => QuadFun a b -> M
identOut (QuadFun m w c) = konst 0 (k-n,n)  === i  where
    n = count @b
    k = rows m
    i = ident n
identIn :: forall a b. BEnum a => QuadFun a b -> M
identIn (QuadFun m w c) = -i === konst 0 (k-n,n)   where
    n = count @a
    k = rows m
    i = ident n
type a :+: b = Either a b
-- I kind of feel like our sign convention is flipped
par :: forall a b c d. (BEnum a, BEnum b, BEnum c, BEnum d) => QuadFun a c -> QuadFun b d -> QuadFun (a :+: b) (c :+: d)
par x@(QuadFun m w c) y@(QuadFun m' w' c') = QuadFun m'' w'' c'' where
    ia = ident (count @a)
    ib = ident (count @b)
    ia' = identIn x
    ib' = identIn y
    ia't = tr ia'
    ib't = tr ib'
    ic = - ident (count @c)
    id = - ident (count @d)
    ic' = identOut x
    id' = identOut y
    ic't = tr ic'
    id't = tr id'
    n = (count @a) + (count @b)
    n' = (count @c) + (count @d)
    --iab = fromBlocks [[0,0,ia,0], [0,0,0,ib], [ia, 0,0,0], [0,ib,0,0]]
    --iab' = fromBlocks [[tr (identIn x),0],  [0 , tr (identIn y)]]
    --zab = konst 0 (n,n)
    -- should be 2 ins + 2 langrange + 2 matr + 2lagarne + 2 outs = 10x10 block matrix
    m'' = fromBlocks [[0 ,0 , ia, 0, 0, 0, 0 ,0,0,0],
                    [0, 0,  0, ib, 0, 0, 0, 0,0,0 ] ,
                    [ia, 0, 0, 0,  ia't, 0,0,0,0, 0],
                    [0, ib, 0, 0,  0,ib't,0, 0,0 ,0],
                    [0, 0, ia',0,  m, 0 , ic',0,0,0],
                    [0, 0,  0, ib',0, m', 0,id',0,0],
                    [0, 0,  0, 0, ic't,0 ,0,0,ic,0],
                    [0, 0,  0, 0,  0,id't,0,0,0,id],
                    [0, 0,  0, 0,  0, 0,  ic,0,0,0],
                    [0, 0,  0, 0,  0, 0,  0,id,0,0]]
    w'' = vjoin [konst 0 (2*n), w,w', konst 0 (2*n')]
    c'' = c + c'

dup :: forall a. BEnum a => QuadFun a (a,a)
dup = QuadFun m 0 0 where -- 1 input, 2 lagrange mutipliers, and 2 outputs.
    ia = ident (count @a)
    m = fromBlocks [ [0, ia,ia, 0,0],
                     [ia, 0, 0, -ia,0],
                     [0, ia, 0, 0, -ia],
                     [0,-ia, 0, 0, 0 ],
                     [0, 0, -ia, 0,0]]

-- fst...?
-- for snd, we could just leave it alone
-- fuse?
-- swap.
data Void
fuse :: forall a. BEnum a => QuadFun (a,a) Void
fuse = QuadFun m 0 0 where -- 2 inputs, 1 lagrange multiplier.
    ia = ident (count @a)
    m = fromBlocks [[0,0,    ia],
                    [0 ,0 , -ia],
                    [ia, -ia, 0]]

-- The analog for this slicing for convex sets may be to slice the
-- dimensions into in and out dimensions. Then

data Cell a = Cell {phi :: Double, j :: Double, next :: a} deriving (Show, Traversable, Foldable, Functor) -- composition of cell is a statically sized vector. Derive applicative?

instance Applicative Cell where
    pure x = Cell 0 0 x
    (Cell phi1 j1 f) <*> (Cell phi2 j2 x) = Cell (phi1 + phi2) (j1 + j2) (f x) -- ? I doubt this makes any sense.
{-
data Cell a = Cell { phi :: a, j :: a}
type f :+: g = Product f g
gaussK :: Lens ((Cell :+: Cell) a) (Cell a)

gaussK :: -> Lens (Cell :*: f a) -- no Cell is very not this.

type Row = (Cell Cell Cell Cell Cell)

Lens (a, b, other) (b, other)

gaussK :: Lens

data Cell2 f a b = Cell2 {phi :: a , j :: a, next :: f b}
data Cell2 f a = Cell2 {phi :: f Double a , j ::f Double, next :: f a}

-- zipA :: f a -> f b -> f (a,b)
-- zipA x y = (,) <$> x <*> y
-- monProd

fmap2 = fmap . fmap
fmap4 = fmap2 . fmap2
fmap8 = fmap4 . fmap4

zipA x y = (,) <$> x <*> y
zip2 = zipA . (fmap zipA) -- hmm. Maybe we should be using Compose. we're going to need famp2. Compose will get us these instances for free.
zip4 = zip2 . (fmap2 zip2)
zip8 = zip4 . (fmap4 zip4)

parK :: Lens (f b) b -> Lens (g a) a -> Lens (f g a) a
parK :: forall b. Lens (f b) (f' b) -> Lens (g a) (g' a) -> Lens (f g a) (f' (g' a))
parK :: forall b. Lens (f g a) (f' g a) -> Lens (g a) (g' a) -> Lens (f g a) (f' (g' a))
parK = compose (fmap l2) l1

Do the y direction inductively

ydir :: Lens (Cell a, Cell a) (a,a)
ydir = \(Cell phi j x, Cell phi2 j2 y) -> ((x,y), \(x',y') -> (Cell phi j x, Cell phi2 j2 y)      )

ydir :: Lens (Cell a, Cell a) (a,a)
ydir = \(a,b) -> ((Cell phi j x) , f) = (gaussK a) in ((Cell phi j y) , g) = (gaussK b) in

Lens ( (X (), (X (), a) (X (), a)
ydir = \(r1, (r2, z) ->  zip8    )    -- = gaussK r1 in  = gaussK r2 in

Lens (f a) a -> Lens (f a) a -> Lens (f a, f a) a
Lens (Cell a) a -> Lens

Lens x y x' y'

-}

{-
can i just use regular lenses? I need to set in kind of a weird way though.
-}

data Lens a b = Lens (a -> (b, b -> a))
-- newtype Lens a b = Lens forall r. (a -> (b -> (b -> a) -> r) -> r) -- cpsify for speed? van Laarhoven?

-- SLens s a b = SLens (a -> ST s (b, b -> ST s a)) mutable update
-- MLens a b = MLens (a -> m (b , b -> m a)) -- monadic lens
-- KLens a b = KLens (a -> (b, k b a)) ---- this is sort of what Conal wrote.
comp :: Lens b c -> Lens a b -> Lens a c
comp (Lens f) (Lens g) = Lens $ \x -> let (y, g') = g x in
                                      let (z, f') = f y in
                                      (z, g' . f')

-- Gauss Seidel of the standard 1 -2 1 matrix
-- j1 is the effective source incliuding the influence from phi values up in the stack
-- Double Cell construction is rather wonky
-- gaussK :: Vec (S (S n)) Double -> Vec (S n)
-- we could also go for the lagrange mutiplier interface.
-- can also push the rho value needed for stability up and down.

-- -2 phi1 + 1 phi2 = ~j1
-- -1 phi1 + -2 phi2 + ...? = j2

{-
We are moving the lower diagonal to the right hand side as the splitting.
That ends up to claculating an effective j based on previous values.
We then triangular solve the upper diagonal, which we are able to find the new diagonal element in terms of the lower values.

-}
-- it's kind of weird that we mutate j on the way down but restore it on the way up.
-- But we do need a way to access j. (Phi (Phi a), J a) (Phi )

gaussK :: Lens (Cell (Cell a)) (Cell a) -- we need the context Cell
gaussK = Lens $ \case (Cell phi1 j1 (Cell phi2 j2 y)) -> (Cell phi2 (j2 - phi1) y , \case (Cell phi2' j2' z) -> let j1' = j1 - phi2' in  -- moving the triangular upsolve to the right hand side
                                                                                Cell (- j1' / 2) j1 (Cell phi2' j2 z))

-- interface in , internals, interface out. was the previous way of talking about it. But then the internals grow, which is fine
-- compose :: Lens in internal out -> Lens in' internal' out' -> Lens in (internal, out, in', internal') out
-- some ind of 2 category? Enriched? The morphisms have this internal structure.
-- this is more a containing relationship. The larger context can be converted into the smaller context.
--

g2 :: Lens (Cell (Cell (Cell a))) (Cell a)
g2 = gaussK `comp` gaussK
g4 = g2 `comp` g2
g8 = g4 `comp` g4
g16 = g8 `comp` g8
g32 = g16 `comp` g16

capZero :: Lens (Cell ()) () -- not even sure i really need this? does runGauss (f `comp` capZero) ~ runGauss f
capZero = Lens $ \case (Cell phi j _) -> ((), \_ -> Cell (- j / 2) j ())

-- runGauss ::  Lens a b -> a -> a -- removes the open ended nature of the thing.
runGauss ::  Lens a () -> a -> a -- this might make more sense as what you really want. This is some kind of cap operation.
runGauss (Lens f) x = let (y , trisolve)  = f x in trisolve y
startingVal :: Cell (Cell (Cell ()))
startingVal = (Cell 0 0 (Cell 0 1 (Cell 0 0 ())))

iters = iterate (runGauss (capZero `comp` g2)) startingVal
--iters' = iterate (runGauss g2) startingVal -- They are different.

sol = ((3><3) [-2,1,0,
               1,-2,1,
               0,1,-2])  <\> (vector [0,1,0])

-- can do inner block solves also (via an iterative method perhaps)
-- this is reminsecent of a multiscale attack.
parC :: Lens a a' -> Lens b b' -> Lens (a,b) (a',b')
parC (Lens f) (Lens g) = Lens $ \(x,y) -> let (x', f') = f x in
                                          let (y', g') = g y in
                                          ((x',y') , f' *** g')

-- profucntor optics paper
data FunList a b t = Done t | More a (FunList a b (b -> t))
-- contents and fill
-- contents (s -> a^n)
-- fill s -> b^n -> t
-- s -> exists n. (a^n, a^n -> s) is a traversal. Could use a Vec. But we already have a Vec. Cell is a Vec
-- s -> exists n, (Vec n a, Vec n a -> s)
-- maybe we do need an applicative. We sort of need to zip together two rows to start going 2d.

-- huh. A block wise schur is kind of the monoidal ppoduct here. No. monoidal product is pure dsum.
-- actually composition? is kind of what plays the game of block wise stuff.
{-

if we want this to be not the case,
comp :: Lens (a,a') a'

schur :: Lens a a' -> (a -> b) -> (b -> a) -> Lens b b' -> Lens (a,b) (a',b')
schur (Lens a) b c (Lens d) = Lens $ \(x,y) -> let (x', a') = a x in
                                               let (y', d') = d x in
                                               (       )
-}
-- Lens (f (g a)) (g a)

                                               -- we can change this all into a continuation form
--

{-
par :: forall a b c d. (BEnum a, BEnum b, BEnum c, BEnum d) => BMatrix a b ->
                         BMatrix c d -> BMatrix (a,c) (b,d)
par q@(QuadFun m w r) q'@(QuadFun m' w' r') = QuadFun q'' w'' r'' where
     a' = diagBlock [sliceA q, sliceA q']
     b' = diagBlock [sliceB q, sliceB q']
     c' = diagBlock [sliceC q, sliceC q']
     d' = diagBlock [sliceD q, sliceD q']
     q'' = fromBlocks [[a',b'], [c',d']]
     u' = vjoin [sliceU v, sliceU v']
     v' = vjoin [sliceV v, sliceV v']
     w'' = vjoin [u', v']
     r'' = r + r'

sliceA :: forall a b. (BEnum a, BEnum b) => QuadFun a b -> M
sliceA (QuadFun m v c) = subMatrix (0,0) (count @a, count @a)  m
sliceB :: forall a b. (BEnum a, BEnum b) => QuadFun a b -> M
sliceB (QuadFun m v c) = subMatrix (0,count @a) (count @a, count @b) m
sliceC :: forall a b. (BEnum a, BEnum b) => QuadFun a b -> M
sliceC (QuadFun m v c) = subMatrix (count @a,0) (count @b, count @a) m
sliceD :: forall a b. (BEnum a, BEnum b) => QuadFun a b -> M
sliceD (QuadFun m v c) = subMatrix (count @a,count @a) (count @b, count @b) m

sliceU :: forall a b. (BEnum a, BEnum b) => QuadFun a b -> Vector Double
sliceU (QuadFun m v c) = subVector 0 (count @a) v
sliceV :: forall a b. (BEnum a, BEnum b) => QuadFun a b -> Vector Double
sliceV (QuadFun m v c) = subVector (count @a) (count @b) v

n = count @a
m = count @b
q = rows m
p = q - m
a = subMatrix (0,0) (count @a, count @a)  m

a = subMatrix
--Other options
-- Keep a matrix with the end and the beginning are the

-- Keep Constraints too?

-- what in the world
-- stack build --ghc-options /usr/lib/libiconv.dylib

-- size
-- rows
-- cols

compose :: QuadFun b c -> QuadFun a b -> QuadFun a c
compose (QuadFun m w c) (QuadFun m' w' c') = QuadFun m'' w'' c'' where
    n = (count @b)
    corner =  fromBlocks [[ident n, 0],[0, -(ident n)]]
    corner' = fromBlocks [[0,0],[corner,0]]
    m'' = fromBlocks [[m ,corner'], [tr corner', m']]
    w'' = vjoin [w, const 0 n, w']
    c'' = c + c

par (QuadFun m w c) (QuadFun m' w' c') =

type a :+: b = Either a b
data QuadFun x = QuadFun

((a,b,c),
 (d,e,f),
 (g,h,i)) =

 data BMatrix' = BMatrix' M M M M M M M M M
 compose (BMatrix' a' b' c' d' e' f' g' h' i') (BMatrix' a b c d e f g h i) =
    BMatrix' a'' b'' c'' d'' e'' f'' g'' h'' i'' where
    a'' = a
    b'' = b ||| c ||| 0
    c'' = 0
    d'' = d === g === 0
    e'' = fromBlocks [[e,f,0], [h, i + a', b'], [0, d', e']]
    f'' = 0 === c' === f'
    g'' =  0
    h'' = 0 ||| g' ||| h'
    i'' = i'
compose' (BVec u' v' w') (BVec u v w) = BVec u (vjoin [v ,w + u', v']) w'

par :: QuadFun a b -> QuadFun c d -> QuadFun (Either a c) (Either b d)
par (BMatrix' a' b' c' d' e' f' g' h' i') (BMatrix' a b c d e f g h i) =
   BMatrix' a'' b'' c'' d'' e'' f'' g'' h'' i'' where
   a'' = diagBlock [a',a]
   b'' = diagBlock [b',b]
   c'' = diagBlock [c',c]
   d'' = diagBlock [d',d]
   e'' = diagBlock [e',e]
   f'' = diagBlock [f',f]
   g'' = diagBlock [g',g]
   h'' = diagBlock [h',h]
   i'' = diagBlock [i',i]

-- If we're just recording the entire program, we might as well just record it in a DSL
-- rather than directly building the matrix
-- data QuadFunDSL = Par QuadFun
-- FreeCat + (Lit HMatrix)

-- iterative matrix solution proximal method?
--
-- gauss seidel
type VBlock = [Vector Double]
-- (V, V)
-- Maybe schur solve works. But with possible splitting
-- Maybe unsolvable schur
-- could damp a little if the diagonal isn't dominant enough
-- (lusolve a (w - V <#| v', ) -- WHat I'm suggesting here is block Jacobi? No
-- No it isn't. Because we're passing back v' it is guass seidel.
-- in order to do block gauss seidel
-- type L = VBlock -> (VBlock, VBlock -> VBlock)

-- Writing in the monoidal categroy style makes parallelism more apparent.

-- [vjoin blocks] -> ([lesser blocks] , update [lesserblocks] -> [blocks]
-- \x:xs -> (a*x + xs, \b -> b )

-- storage in the lens?
-- In my AD lens, I had all the weights in the input. This was ungainly
-- Maybe a choice monad? Everyone could have access to a global state
-- forall s. (ActualLens s a, a -> (b, db -> da))
-- compsoe (ActualLens s a, a -> (b, db -> da))
-- par
-- ActualLens are functional references. Do they allow us to refiy sharing?
-- (ActualLens acc a,   )

-- alternative a -> State s (b, db -> State ds da)

-- Mutiple directions of composition. Horizontsl and vertical?
-- a 2-category?
-- ActualLens s a -> ActualLens s' s

-- Jules Hedges:
-- Lens sigma a b -> Lens sigm' b c -> Lens (sig,sig') a c
-- Lens :: Sigma -> (sigma,  ) -> Lens a b -- hides sigma, we'll never be able to get at it

-}

ext: hs


text:
{-# language
        BangPatterns
      , CPP
      , GeneralizedNewtypeDeriving
      , InstanceSigs
      , MagicHash
      , ScopedTypeVariables
      , StandaloneDeriving
      , UnboxedTuples
  #-}

#if __GLASGOW_HASKELL__ >= 800
{-# language TypeInType #-}
#endif

{-# OPTIONS_GHC
      -Weverything
      -fno-warn-unsafe
      -fno-warn-orphans
      -fno-warn-name-shadowing
      -fno-warn-missing-import-lists
      -fno-warn-implicit-prelude
      -O2
#-}

-- | Orphan instances for the 'Prim' typeclass.
module Data.Primitive.Instances () where

import Data.Complex (Complex(..))
#if !MIN_VERSION_primitive(0,7,0)
import Data.Functor.Const (Const(..))
#if MIN_VERSION_base(4,8,0)
import Data.Functor.Identity (Identity(..))
#endif
import qualified Data.Monoid as Monoid
import qualified Data.Semigroup as Semigroup
import Data.Ord (Down(..))
#endif
import Data.Primitive.Types (Prim(..), defaultSetOffAddr#, defaultSetByteArray#)
import GHC.Real (Ratio(..))
import Data.Word (Word64)
import GHC.Fingerprint (Fingerprint(..))
import GHC.Exts (State#, Int#, Addr#, MutableByteArray#, (+#), (*#))

instance Prim a => Prim (Complex a) where
  sizeOf# _ = 2# *# sizeOf# (undefined :: a)
  alignment# _ = alignment# (undefined :: a)
  indexByteArray# arr# i# =
    let x,y :: a
        x = indexByteArray# arr# (2# *# i#)
        y = indexByteArray# arr# (2# *# i# +# 1#)
    in x :+ y
  readByteArray# :: forall s a. (Prim a) => MutableByteArray# s -> Int# -> State# s -> (# State# s, Complex a #)
  readByteArray# arr# i# =
    \s0 -> case readByteArray# arr# (2# *# i#) s0 of
      (# s1#, x #) -> case readByteArray# arr# (2# *# i# +# 1#) s1# of
        (# s2#, y #) -> (# s2#, x :+ y #)
  writeByteArray# :: forall s a. (Prim a) => MutableByteArray# s -> Int# -> Complex a -> State# s -> State# s
  writeByteArray# arr# i# (a :+ b) =
    \s0 -> case writeByteArray# arr# (2# *# i#) a s0 of
      s1 -> case writeByteArray# arr# (2# *# i# +# 1#) b s1 of
        s2 -> s2
  setByteArray# = defaultSetByteArray#
  indexOffAddr# :: Addr# -> Int# -> Complex a
  indexOffAddr# addr# i# =
    let x,y :: a
        x = indexOffAddr# addr# (2# *# i#)
        y = indexOffAddr# addr# (2# *# i# +# 1#)
    in x :+ y
  readOffAddr# :: forall s a. (Prim a) => Addr# -> Int# -> State# s -> (# State# s, Complex a #)
  readOffAddr# addr# i# =
    \s0 -> case readOffAddr# addr# (2# *# i#) s0 of
      (# s1, x #) -> case readOffAddr# addr# (2# *# i# +# 1#) s1 of
        (# s2, y #) -> (# s2, x :+ y #)
  writeOffAddr# :: forall s a. (Prim a) => Addr# -> Int# -> Complex a -> State# s -> State# s
  writeOffAddr# addr# i# (a :+ b) =
    \s0 -> case writeOffAddr# addr# (2# *# i#) a s0 of
      s1 -> case writeOffAddr# addr# (2# *# i# +# 1#) b s1 of
        s2 -> s2
  setOffAddr# = defaultSetOffAddr#
  {-# INLINE sizeOf# #-}
  {-# INLINE alignment# #-}
  {-# INLINE indexByteArray# #-}
  {-# INLINE readByteArray# #-}
  {-# INLINE writeByteArray# #-}
  {-# INLINE setByteArray# #-}
  {-# INLINE indexOffAddr# #-}
  {-# INLINE readOffAddr# #-}
  {-# INLINE writeOffAddr# #-}
  {-# INLINE setOffAddr# #-}

instance (Integral a, Prim a) => Prim (Ratio a) where
  sizeOf# _ = 2# *# sizeOf# (undefined :: a)
  alignment# _ = alignment# (undefined :: a)
  indexByteArray# arr# i# =
    let x,y :: a
        x = indexByteArray# arr# (2# *# i#)
        y = indexByteArray# arr# (2# *# i# +# 1#)
    in x :% y
  readByteArray# :: forall s a. (Prim a) => MutableByteArray# s -> Int# -> State# s -> (# State# s, Ratio a #)
  readByteArray# arr# i# =
    \s0 -> case readByteArray# arr# (2# *# i#) s0 of
      (# s1#, x #) -> case readByteArray# arr# (2# *# i# +# 1#) s1# of
        (# s2#, y #) -> (# s2#, x :% y #)
  writeByteArray# :: forall s a. (Prim a) => MutableByteArray# s -> Int# -> Ratio a -> State# s -> State# s
  writeByteArray# arr# i# (a :% b) =
    \s0 -> case writeByteArray# arr# (2# *# i#) a s0 of
      s1 -> case writeByteArray# arr# (2# *# i# +# 1#) b s1 of
        s2 -> s2
  setByteArray# = defaultSetByteArray#
  indexOffAddr# :: Addr# -> Int# -> Ratio a
  indexOffAddr# addr# i# =
    let x,y :: a
        x = indexOffAddr# addr# (2# *# i#)
        y = indexOffAddr# addr# (2# *# i# +# 1#)
    in x :% y
  readOffAddr# :: forall s a. (Prim a) => Addr# -> Int# -> State# s -> (# State# s, Ratio a #)
  readOffAddr# addr# i# =
    \s0 -> case readOffAddr# addr# (2# *# i#) s0 of
      (# s1, x #) -> case readOffAddr# addr# (2# *# i# +# 1#) s1 of
        (# s2, y #) -> (# s2, x :% y #)
  writeOffAddr# :: forall s a. (Prim a) => Addr# -> Int# -> Ratio a -> State# s -> State# s
  writeOffAddr# addr# i# (a :% b) =
    \s0 -> case writeOffAddr# addr# (2# *# i#) a s0 of
      s1 -> case writeOffAddr# addr# (2# *# i# +# 1#) b s1 of
        s2 -> s2
  setOffAddr# = defaultSetOffAddr#
  {-# INLINE sizeOf# #-}
  {-# INLINE alignment# #-}
  {-# INLINE indexByteArray# #-}
  {-# INLINE readByteArray# #-}
  {-# INLINE writeByteArray# #-}
  {-# INLINE setByteArray# #-}
  {-# INLINE indexOffAddr# #-}
  {-# INLINE readOffAddr# #-}
  {-# INLINE writeOffAddr# #-}
  {-# INLINE setOffAddr# #-}

instance Prim Fingerprint where
  sizeOf# _ = 2# *# sizeOf# (undefined :: Word64)
  alignment# _ = alignment# (undefined :: Word64)
  indexByteArray# arr# i# =
    let x,y :: Word64
        x = indexByteArray# arr# (2# *# i#)
        y = indexByteArray# arr# (2# *# i# +# 1#)
    in Fingerprint x y
  readByteArray# :: forall s. MutableByteArray# s -> Int# -> State# s -> (# State# s, Fingerprint #)
  readByteArray# arr# i# =
    \s0 -> case readByteArray# arr# (2# *# i#) s0 of
      (# s1#, x #) -> case readByteArray# arr# (2# *# i# +# 1#) s1# of
        (# s2#, y #) -> (# s2#, Fingerprint x y #)
  writeByteArray# :: forall s. MutableByteArray# s -> Int# -> Fingerprint -> State# s -> State# s
  writeByteArray# arr# i# (Fingerprint a b) =
    \s0 -> case writeByteArray# arr# (2# *# i#) a s0 of
      s1 -> case writeByteArray# arr# (2# *# i# +# 1#) b s1 of
        s2 -> s2
  setByteArray# = defaultSetByteArray#
  indexOffAddr# :: Addr# -> Int# -> Fingerprint
  indexOffAddr# addr# i# =
    let x,y :: Word64
        x = indexOffAddr# addr# (2# *# i#)
        y = indexOffAddr# addr# (2# *# i# +# 1#)
    in Fingerprint x y
  readOffAddr# :: forall s. Addr# -> Int# -> State# s -> (# State# s, Fingerprint #)
  readOffAddr# addr# i# =
    \s0 -> case readOffAddr# addr# (2# *# i#) s0 of
      (# s1, x #) -> case readOffAddr# addr# (2# *# i# +# 1#) s1 of
        (# s2, y #) -> (# s2, Fingerprint x y #)
  writeOffAddr# :: forall s. Addr# -> Int# -> Fingerprint -> State# s -> State# s
  writeOffAddr# addr# i# (Fingerprint a b) =
    \s0 -> case writeOffAddr# addr# (2# *# i#) a s0 of
      s1 -> case writeOffAddr# addr# (2# *# i# +# 1#) b s1 of
        s2 -> s2
  setOffAddr# = defaultSetOffAddr#
  {-# INLINE sizeOf# #-}
  {-# INLINE alignment# #-}
  {-# INLINE indexByteArray# #-}
  {-# INLINE readByteArray# #-}
  {-# INLINE writeByteArray# #-}
  {-# INLINE setByteArray# #-}
  {-# INLINE indexOffAddr# #-}
  {-# INLINE readOffAddr# #-}
  {-# INLINE writeOffAddr# #-}
  {-# INLINE setOffAddr# #-}

#if !MIN_VERSION_primitive(0,7,0)

deriving instance Prim a => Prim (Down a)
#if MIN_VERSION_base(4,8,0)
deriving instance Prim a => Prim (Identity a)
deriving instance Prim a => Prim (Monoid.Dual a)
deriving instance Prim a => Prim (Monoid.Sum a)
deriving instance Prim a => Prim (Monoid.Product a)
#endif
#if MIN_VERSION_base(4,9,0)
deriving instance Prim a => Prim (Semigroup.First a)
deriving instance Prim a => Prim (Semigroup.Last a)
deriving instance Prim a => Prim (Semigroup.Min a)
deriving instance Prim a => Prim (Semigroup.Max a)
#endif
deriving instance Prim a => Prim (Const a b)

#endif

max_stars_repo_name: nprindle/prim-instances

max_issues_repo_name: nprindle/prim-instances

max_forks_repo_name: nprindle/prim-instances


text:
-- |
-- Module      : Copernicus
-- Description : Contains the core of the Copernicus project
-- Copyright   : (c) Jonatan H Sundqvist, 2015
-- License     : MIT
-- Maintainer  : Jonatan H Sundqvist
-- Stability   : experimental
-- Portability : POSIX (not sure)
--
-- Jonatan H Sundqvist
-- November 25 2014
--

-- TODO | - Model planetary motion and physics in 2D
--        - Cartoon earth with gravitational field
--        - Flexible event handling
--        - App typeclass (run, manage events, window properties, etc.)
--        - Lenses
--        - Options (eg. toggle grid) (cf. 'when')
--        - UI
--        - Use types to encode units (eg. SI, radians)
--        - Move to Cairo (branch?)
--        - Polymorphic types (not just Floats) (?)

-- SPEC | -
--        -

module Copernicus where

---------------------------------------------------------------------------------------------------
-- We'll need these
---------------------------------------------------------------------------------------------------
import Data.Complex
-- import Control.Monad (when)

import Graphics.Gloss.Geometry.Angle (degToRad, radToDeg, normaliseAngle)

---------------------------------------------------------------------------------------------------
-- Types
---------------------------------------------------------------------------------------------------
-- Make typeclass (mesh, body, bounding box, collision, etc) (?)
-- type Number a = Floating a -- Real number
type Vector = Complex -- TODO: Polymorphic (cf. related TODO item)
data Body f = Body (Vector f) (Vector f) (Vector f) deriving Show -- Add argument

---------------------------------------------------------------------------------------------------
-- Functions
---------------------------------------------------------------------------------------------------
clampAngle :: Float -> Float
clampAngle = radToDeg . normaliseAngle . degToRad

--
-- TODO: Make them change colour when bouncing (?)
animate :: (RealFloat f, Floating f) => f -> Body f -> Body f
animate t (Body p v a) = collide ground $ Body (parabola t p v a) (v + (t:+0)*a) a
	where ground = 0.0 -- 30+30/2-540/2 (previous hard-coded value)

--
parabola :: (RealFloat f, Floating f) => f -> Vector f -> Vector f -> Vector f -> Vector f
parabola t p v a = let (px:+py) = p
                       (vx:+vy) = v
                       (ax:+ay) = a
                   in (px + vx*t + 0.5*ax*t**2) :+ (py + vy*t + 0.5*ay*t**2)

-- collide
-- Very primitive for now
-- TODO: Use 'contains' function (Range -> Value -> Bool)
-- TODO: Don't hard-code bounds (left, right)
-- TODO: Take bounds of Body into account (don't hard-code that either)
collide :: (RealFloat f, Floating f) => f -> Body f -> Body f
collide gnd (Body (px:+py) (vx:+vy) a) = Body (px:+py) ((invertIf (\ _ -> (px <= left) || ( px >= right)) vx) :+ (invertIf (\ v -> (v < 0) && (py <= gnd)) vy)) a
	where invertIf p v | p v 	   = -v
	                   | otherwise =  v
	      (left, right) = (-5, 5) --(15-720/2, 720/2-30/2)

-- | ETA (estimated time of arrival)
-- TODO: Rename (eg. timeUntil, solveForT, etc)
-- TODO: Parabola type (eg. Parabola a v x)
-- eta :: Acceleration

-- Utilities (should eventually be moved to separate module or library ----------------------------
-- Python-style String formatting (eg. keyword interpolation, {0}, customisation, format specs.)
-- Parsec

---------------------------------------------------------------------------------------------------
-- Entry point
---------------------------------------------------------------------------------------------------

max_stars_repo_name: SwiftsNamesake/Copernicus

max_issues_repo_name: SwiftsNamesake/Copernicus

max_forks_repo_name: SwiftsNamesake/Copernicus


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE NoMonomorphismRestriction #-}
{-# OPTIONS_GHC -fno-warn-type-defaults #-}
{-# OPTIONS_GHC -fno-warn-type-defaults -fno-warn-orphans -Wall #-}

module Main where

import Algebra.Bridge.Singular
import Algebra.Field.Prime
import Algebra.Prelude.Core
import Cases
import Control.Monad (void)
import qualified Data.Text as T
import qualified Data.Vector.Unboxed as V
import Statistics.Resampling
import Statistics.Resampling.Bootstrap
import Statistics.Types
import qualified System.Random.MWC as Rand
import Prelude (read)

benchIdeal :: IsSingularPolynomial poly => Text -> Ideal poly -> IO Double
benchIdeal fun i =
  fmap ((/ 1000) . read . T.unpack) $
    singular $
      toProg fun i

toProg :: (SingularOrder (Arity poly) (MOrder poly), SingularCoeff (Coefficient poly), IsOrderedPolynomial poly) => Text -> Ideal poly -> Text
toProg fun i =
  let expr = funE fun [idealE' i]
   in prettySingular $ do
        directC "system(\"--ticks-per-sec\", 1000000)"
        void $ ringC "R" expr
        declOnlyC IdealT "G"
        directC "timer=1"
        directC "int t=rtimer"
        letC "G" expr
        directC "print(rtimer-t)"
        directC "exit"

analyse :: (IsSingularPolynomial poly) => String -> Text -> Ideal poly -> IO ()
analyse lab fun ideal = do
  gen <- Rand.create
  i2Gr <- V.replicateM 50 $ benchIdeal fun ideal
  res <- resample gen [Mean, StdDev] 1000 i2Gr
  let [Estimate mn (ConfInt lbmn ubmn _), Estimate dv _] =
        bootstrapBCA cl95 i2Gr res
  putStrLn lab
  mapM_
    (putStrLn . ('\t' :))
    [ "Mean:\t" ++ show mn ++ "(ms)"
    , "MeanLB:\t" ++ show lbmn
    , "MeanUB:\t" ++ show ubmn
    , "StdDev:\t" ++ show dv
    ]

runTestCases ::
  (IsMonomialOrder n o, KnownNat n) =>
  String ->
  Ideal (OrderedPolynomial Rational o n) ->
  IO ()
runTestCases lab i = do
  analyse (lab ++ " (Q, Lex, Sing(gr))") "groebner" $ fmap (changeOrder Lex) i
  analyse (lab ++ " (Q, Lex, Sing(sba))") "sba" $ fmap (changeOrder Lex) i
  analyse (lab ++ " (Q, Grevlex, Sing(gr))") "groebner" $ fmap (changeOrder Grevlex) i
  analyse (lab ++ " (Q, Grevlex, Sing(sba))") "sba" $ fmap (changeOrder Grevlex) i
  analyse (lab ++ " (F_65521, Lex, Sing(gr))") "groebner" $ fmap (mapCoeff ratToF . changeOrder Lex) i
  analyse (lab ++ " (F_65521, Lex, Sing(sba))") "sba" $ fmap (mapCoeff ratToF . changeOrder Lex) i
  analyse (lab ++ " (F_65521, Grevlex, Sing(gr))") "groebner" $ fmap (mapCoeff ratToF . changeOrder Grevlex) i
  analyse (lab ++ " (F_65521, Grevlex, Sing(sba))") "sba" $ fmap (mapCoeff ratToF . changeOrder Grevlex) i

main :: IO ()
main = do
  runTestCases "Cyclic-4" $ cyclic (sNat :: SNat 4)
  runTestCases "Cyclic-5" $ cyclic (sNat :: SNat 5)
  runTestCases "Cyclic-6" $ cyclic (sNat :: SNat 6)
  runTestCases "Katsura-5" $ katsura (sNat :: SNat 5)
  runTestCases "Katsura-6" $ katsura (sNat :: SNat 6)
  runTestCases "Katsura-7" $ katsura (sNat :: SNat 7)

ratToF :: Rational -> F 65521
ratToF = modRat'

max_stars_repo_name: hangingman/computational-algebra

max_issues_repo_name: hangingman/computational-algebra

max_forks_repo_name: hangingman/computational-algebra


text:
-------------------------------------------------------------------------------

  {-  LANGUAGE CPP #-}  -- specified in .cabal default-extensions

#define USE_NFDATA_SUPERCLASS 0

-------------------------------------------------------------------------------

  -- Later: I'm not so sure about this, actually; is the arithmetic
  -- on n actually not piling up thunks, without the bang-patterns?
  --
  -- It would be easy to get rid of the bang-patterns.
  -- The Complex instance is done as an example.
  -- If you do go with the case's, probably want -fno-warn-name-shadowing.
  {-# LANGUAGE BangPatterns #-}
  {-# OPTIONS_GHC -fno-warn-name-shadowing #-}

-------------------------------------------------------------------------------

-- |
-- Module      :  Control.DeepSeq.Bounded.NFDataN
-- Copyright   :  Andrew G. Seniuk 2014-2015
-- License     :  BSD-style (see the file LICENSE)
--
-- Maintainer  :  Andrew Seniuk <rasfar@gmail.com>
-- Stability   :  provisional
-- Portability :  portable
--
-- This module provides an overloaded function, 'deepseqn', for partially
-- (or fully) evaluating data structures to a given depth.
--

-------------------------------------------------------------------------------

  module Control.DeepSeq.Bounded.NFDataN
  (

      -- * Depth-bounded analogues of 'deepseq' and 'force'

      deepseqn
    , forcen

      -- * Class of things that can be evaluated to an arbitrary finite depth

#if JUST_ALIAS_GNFDATAN
    , rnfn
#else
    , NFDataN(..)
#endif

  )
  where

-------------------------------------------------------------------------------

--import qualified Data.Generics.Shape as S  -- not actually used

#if JUST_ALIAS_GNFDATAN
  import Generics.SOP ( Generic )
  import Control.DeepSeq.Bounded.Generic.GNFDataN ( grnfn )
#endif

#if USE_NFDATA_SUPERCLASS
  import Control.DeepSeq
#endif

--import Data.Data

  import Data.Int
  import Data.Word
  import Data.Ratio
  import Data.Complex
  import Data.Array
  import Data.Fixed
  import Data.Version

-------------------------------------------------------------------------------

#if JUST_ALIAS_GNFDATAN
  type NFDataN = Generic
#endif

#if JUST_ALIAS_GNFDATAN
  rnfn = grnfn
#endif

-------------------------------------------------------------------------------

-- infixr 0 $!!

-------------------------------------------------------------------------------

  -- | @'deepseqn' n x y@ evaluates @x@ to depth n, before returning @y@.
  --
  -- This is used when expression @x@ also appears in @y@, but mere
  -- evaluation of @y@ does not force the embedded @x@ subexpression
  -- as deeply as we wish.
  --
  -- The name 'deepseqn' is used to illustrate the relationship to 'seq':
  -- where 'seq' is shallow in the sense that it only evaluates the top
  -- level of its argument, @'deepseqn' n@ traverses (evaluates) the entire
  -- top @n@ levels of the data structure.
  --
  -- A typical use is to ensure any exceptions hidden within lazy
  -- fields of a data structure do not leak outside the scope of the
  -- exception handler; another is to force evaluation of a data structure in
  -- one thread, before passing it to another thread (preventing work moving
  -- to the wrong threads). Unlike <http://hackage.haskell.org/package/deepseq/docs/Control-DeepSeq.html DeepSeq>, potentially infinite values of coinductive
  -- data types are supported by principled bounding of deep evaluation.
  --
  -- It is also useful for diagnostic purposes when trying to understand
  -- and manipulate space\/time trade-offs in lazy code,
  -- and as a less indiscriminate substitute for 'deepseq'.
  --
  -- Furthermore, 'deepseqn' can sometimes be a better solution than
  -- using stict fields in your data structures, because the
  -- latter will behave strictly everywhere that its constructors
  -- are used, instead of just where its laziness is problematic.
  --
  -- There may be possible applications to the prevention of resource leaks
  -- in lazy streaming, but I'm not certain.
  --
  -- One of the special qualities of 'NFDataN' is shape oblivity:
  -- it doesn't care about any details of the shape of the term
  -- it's forcing, it only cares about stratifying levels of
  -- recursion depth.
  -- (I would say \"as contrasted with <http://hackage.haskell.org/package/deepseq-bounded-0.8.0.0/docs/Control-DeepSeq-Bounded-NFDataP.html NFDataP>\" but cannot, because
  -- <http://hackage.haskell.org/package/deepseq-bounded-0.8.0.0/docs/Control-DeepSeq-Bounded-NFDataP.html NFDataP>
  -- was extended to include
  -- <http://hackage.haskell.org/package/deepseq-bounded-0.8.0.0/docs/Control-DeepSeq-Bounded-NFDataN.html NFDataN>
  -- syntax\/capabilities, precisely to ammend this deficiency.)

  deepseqn :: NFDataN a => Int -> a -> b -> b
  deepseqn n a b = rnfn n a `seq` b
#if 1
  -- XXX Need to double-check that this makes sense; didn't think
  -- it through -- when b != a, things are not so simple.
  -- XXX Partially-applied; is that okay in GHC RULES?
  {-# RULES
    "deepseqn/composition"    forall n1 n2 x.  (.) (deepseqn n2) (deepseqn n1) x = deepseqn ( max n1 n2 ) x
      #-}
#endif

-------------------------------------------------------------------------------

#if 0
-- As per DeepSeq:
  -- | the deep analogue of '$!'.  In the expression @f $!! x@, @x@ is
  -- fully evaluated before the function @f@ is applied to it.
  ($!!) :: NFData a => (a -> b) -> a -> b
  f $!! x = x `deepseq` f x
#endif

-------------------------------------------------------------------------------

  -- | 'forcen' is a variant of 'deepseqn' that is useful in some circumstances.
  --
  -- > forcen n x = deepseqn n x x
  --
  -- @forcen n x@ evaluates @x@ to depth @n@, and then returns it.
  -- Recall that, in common with all Haskell functions, @forcen@ only
  -- performs evaluation when upstream demand actually occurs,
  -- so essentially it turns shallow evaluation into depth-@n@ evaluation.

  forcen :: NFDataN a => Int -> a -> a
  forcen n x = deepseqn n x x
#if 1
  {-# RULES
    "forcen/composition"    forall n1 n2 x.  (.) (forcen n2) (forcen n1) x = forcen ( max n1 n2 ) x
      #-}
#endif

-------------------------------------------------------------------------------

#if ! JUST_ALIAS_GNFDATAN

  -- | A class of types that can be evaluated to arbitrary depth.
#if USE_NFDATA_SUPERCLASS
  class NFData a => NFDataN a where
#else
  class NFDataN a where
#endif
#if 0
    -- | rnf should reduce its argument to normal form (that is, fully
    -- evaluate all sub-components), and then return '()'.
    --
    -- The default implementation of 'rnf' is
    --
    -- > rnf a = a `seq` ()
    --
    -- which may be convenient when defining instances for data types with
    -- no unevaluated fields (e.g. enumerations).
    rnf :: a -> ()
    rnf a = a `seq` ()
#endif
    rnfn :: Int -> a -> ()
#if USE_NFDATA_SUPERCLASS
    -- ?? not sure about this (whatever, we'll override it)
#if 1
    rnfn n a | n <= 0     = ()  -- case needed? (seems to be!)
             | otherwise  = rnf a
#else
    rnfn _ a = rnf a
#endif
#else
#if 0
#elif 1
    rnfn n x | n <= 0     = ()
             | otherwise  = x `seq` ()
#elif 0
    rnfn _ x = x `seq` ()
#elif 0
    rnfn _ _ = ()
#endif
#endif
#if 1
  {-# RULES
    "rnfn/composition"    forall n1 n2 x.  (.) (rnfn n2) (rnfn n1) x = rnfn ( max n1 n2 ) x
      #-}
--   "rnfp/composition"    forall p1 p2 x.  compose (rnfp p2) (rnfp p1) x = rnfp ( unionPat [p1, p2] ) x
--   "rnfp/composition"    forall p1 p2 x.  ( rnfp p2 . rnfp p1 ) x = rnfp ( unionPat [p1, p2] ) x
#endif

-------------------------------------------------------------------------------

  instance NFDataN Int
  instance NFDataN Word
  instance NFDataN Integer
  instance NFDataN Float
  instance NFDataN Double

  instance NFDataN Char
  instance NFDataN Bool
  instance NFDataN ()

  instance NFDataN Int8
  instance NFDataN Int16
  instance NFDataN Int32
  instance NFDataN Int64

  instance NFDataN Word8
  instance NFDataN Word16
  instance NFDataN Word32
  instance NFDataN Word64

  instance NFDataN (Fixed a)

  -- [Quoted from deepseq:]
  -- This instance is for convenience and consistency with 'seq'.
  -- This assumes that WHNF is equivalent to NF for functions.
  instance NFDataN (a -> b)

#if 0
  instance NFDataN a => NFDataN (IO a) where
    rnfn !n x = ()  -- XXX ignore if in IO (cludge easing an experiment...)
#endif

  -- not taken to be a level of depth
  instance (Integral a, NFDataN a) => NFDataN (Ratio a) where
    rnfn !n x = rnfn n (numerator x, denominator x)

  instance (RealFloat a, NFDataN a) => NFDataN (Complex a) where
#if 1
    rnfn n (x:+y) = case n of
     n | n <= 0  -> ()
     _           -> let n' = -1+n in
                    rnfn n' x `seq`
                    rnfn n' y `seq`
                    ()
#else
    rnfn !n (x:+y)
     | n <= 0    = ()
     | otherwise = let n' = -1+n in
                   rnfn n' x `seq`
                   rnfn n' y `seq`
                   ()
#endif

  instance NFDataN a => NFDataN (Maybe a) where
    rnfn _ Nothing  = ()
    rnfn !n (Just x)
     | n <= 0     = ()
     | otherwise  = rnfn (-1+n) x

  instance (NFDataN a, NFDataN b) => NFDataN (Either a b) where
    rnfn !n (Left x)
     | n <= 0     = ()
     | otherwise  = rnfn (-1+n) x
    rnfn !n (Right y)
     | n <= 0     = ()
     | otherwise  = rnfn (-1+n) y

  instance NFDataN Data.Version.Version where
    rnfn !n (Data.Version.Version branch tags)
     | n <= 0     = ()
     | otherwise  = let n' = -1+n in
                    rnfn n' branch `seq` rnfn n' tags

  instance NFDataN a => NFDataN [a] where
    rnfn _ [] = ()
    rnfn !n (x:xs)
     | n <= 0     = ()
     | otherwise  = let n' = -1+n in
                    rnfn n' x `seq` rnfn n' xs

  -- not taken to be a level of depth
  instance (Ix a, NFDataN a, NFDataN b) => NFDataN (Array a b) where
    rnfn !n x = rnfn n (bounds x, Data.Array.elems x)

  instance (NFDataN a, NFDataN b) => NFDataN (a,b) where
    rnfn !n (x,y)
     | n <= 0     = ()
     | otherwise  = let n' = -1+n in
                    rnfn n' x `seq` rnfn n' y

  instance (NFDataN a, NFDataN b, NFDataN c) => NFDataN (a,b,c) where
    rnfn !n (x,y,z)
     | n <= 0     = ()
     | otherwise  = let n' = -1+n in
                    rnfn n' x `seq` rnfn n' y `seq` rnfn n' z

  instance (NFDataN a, NFDataN b, NFDataN c, NFDataN d) => NFDataN (a,b,c,d) where
    rnfn !n (x1,x2,x3,x4)
     | n <= 0     = ()
     | otherwise  = let n' = -1+n in
                    rnfn n' x1 `seq`
                    rnfn n' x2 `seq`
                    rnfn n' x3 `seq`
                    rnfn n' x4

  instance (NFDataN a1, NFDataN a2, NFDataN a3, NFDataN a4, NFDataN a5) =>
           NFDataN (a1, a2, a3, a4, a5) where
    rnfn !n (x1,x2,x3,x4,x5)
     | n <= 0     = ()
     | otherwise  = let n' = -1+n in
                    rnfn n' x1 `seq`
                    rnfn n' x2 `seq`
                    rnfn n' x3 `seq`
                    rnfn n' x4 `seq`
                    rnfn n' x5

  instance (NFDataN a1, NFDataN a2, NFDataN a3, NFDataN a4, NFDataN a5, NFDataN a6) =>
           NFDataN (a1, a2, a3, a4, a5, a6) where
    rnfn !n (x1,x2,x3,x4,x5,x6)
     | n <= 0     = ()
     | otherwise  = let n' = -1+n in
                    rnfn n' x1 `seq`
                    rnfn n' x2 `seq`
                    rnfn n' x3 `seq`
                    rnfn n' x4 `seq`
                    rnfn n' x5 `seq`
                    rnfn n' x6

  instance (NFDataN a1, NFDataN a2, NFDataN a3, NFDataN a4, NFDataN a5, NFDataN a6, NFDataN a7) =>
           NFDataN (a1, a2, a3, a4, a5, a6, a7) where
    rnfn !n (x1,x2,x3,x4,x5,x6,x7)
     | n <= 0     = ()
     | otherwise  = let n' = -1+n in
                    rnfn n' x1 `seq`
                    rnfn n' x2 `seq`
                    rnfn n' x3 `seq`
                    rnfn n' x4 `seq`
                    rnfn n' x5 `seq`
                    rnfn n' x6 `seq`
                    rnfn n' x7

  instance (NFDataN a1, NFDataN a2, NFDataN a3, NFDataN a4, NFDataN a5, NFDataN a6, NFDataN a7, NFDataN a8) =>
           NFDataN (a1, a2, a3, a4, a5, a6, a7, a8) where
    rnfn !n (x1,x2,x3,x4,x5,x6,x7,x8)
     | n <= 0     = ()
     | otherwise  = let n' = -1+n in
                    rnfn n' x1 `seq`
                    rnfn n' x2 `seq`
                    rnfn n' x3 `seq`
                    rnfn n' x4 `seq`
                    rnfn n' x5 `seq`
                    rnfn n' x6 `seq`
                    rnfn n' x7 `seq`
                    rnfn n' x8

  instance (NFDataN a1, NFDataN a2, NFDataN a3, NFDataN a4, NFDataN a5, NFDataN a6, NFDataN a7, NFDataN a8, NFDataN a9) =>
           NFDataN (a1, a2, a3, a4, a5, a6, a7, a8, a9) where
    rnfn !n (x1,x2,x3,x4,x5,x6,x7,x8,x9)
     | n <= 0     = ()
     | otherwise  = let n' = -1+n in
                    rnfn n' x1 `seq`
                    rnfn n' x2 `seq`
                    rnfn n' x3 `seq`
                    rnfn n' x4 `seq`
                    rnfn n' x5 `seq`
                    rnfn n' x6 `seq`
                    rnfn n' x7 `seq`
                    rnfn n' x8 `seq`
                    rnfn n' x9

#endif

-------------------------------------------------------------------------------

max_stars_repo_name: phunehehe/deepseq-bounded

max_issues_repo_name: phunehehe/deepseq-bounded

max_forks_repo_name: phunehehe/deepseq-bounded


text:
{-# LANGUAGE CPP                   #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE UndecidableInstances  #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE OverloadedStrings     #-}
{-# LANGUAGE OverloadedLabels      #-}
{-# OPTIONS_GHC -Wno-incomplete-uni-patterns #-}

module Grenade.Layers.FullyConnected (
    FullyConnected (..)
  , FullyConnected' (..)
  , randomFullyConnected
  ) where

import           Control.DeepSeq
import           Control.Monad.Primitive        (PrimBase, PrimState)
import           GHC.Generics                   (Generic)
import           GHC.TypeLits
import           System.Random.MWC              hiding (create)
import           Data.List                      (foldl1')
import           Data.Either                    (fromRight)
import           Data.Proxy
import           Data.Serialize

import qualified Numeric.LinearAlgebra          as LA
import           Numeric.LinearAlgebra.Static

import           Grenade.Core
import           Grenade.Layers.Internal.Update
import           Grenade.Onnx
import           Grenade.Utils.LinearAlgebra
import           Grenade.Utils.ListStore
import           Grenade.Types

import           Lens.Micro

-- | A basic fully connected (or inner product) neural network layer.
data FullyConnected i o = FullyConnected
                        !(FullyConnected' i o)   -- Neuron weights
                        !(ListStore (FullyConnected' i o))   -- momentum store
                        deriving (Generic)

instance NFData (FullyConnected i o) where
  rnf (FullyConnected w store) = rnf w `seq` rnf store

data FullyConnected' i o = FullyConnected'
                         !(R o)   -- Bias
                         !(L o i) -- Activations
                        deriving (Generic)

instance NFData (FullyConnected' i o) where
  rnf (FullyConnected' b w) = rnf b `seq` rnf w

instance Show (FullyConnected i o) where
  show FullyConnected {} = "FullyConnected"

instance (KnownNat i, KnownNat o, KnownNat (i * o)) => UpdateLayer (FullyConnected i o) where
  type Gradient (FullyConnected i o) = (FullyConnected' i o)
  type MomentumStore (FullyConnected i o) = ListStore (FullyConnected' i o)
  runUpdate opt@OptSGD{} x@(FullyConnected (FullyConnected' oldBias oldActivations) store) (FullyConnected' biasGradient activationGradient) =
    let (FullyConnected' oldBiasMomentum oldMomentum) = getData opt x store
        VectorResultSGD newBias newBiasMomentum = descendVector opt (VectorValuesSGD oldBias biasGradient oldBiasMomentum)
        MatrixResultSGD newActivations newMomentum = descendMatrix opt (MatrixValuesSGD oldActivations activationGradient oldMomentum)
        newStore = setData opt x store (FullyConnected' newBiasMomentum newMomentum)
     in FullyConnected (FullyConnected' newBias newActivations) newStore
  runUpdate opt@OptAdam{} x@(FullyConnected (FullyConnected' oldBias oldActivations) store) (FullyConnected' biasGradient activationGradient) =
    let [FullyConnected' oldMBias oldMActivations, FullyConnected' oldVBias oldVActivations] = getData opt x store
        VectorResultAdam newBias newMBias newVBias = descendVector opt (VectorValuesAdam (getStep store) oldBias biasGradient oldMBias oldVBias)
        MatrixResultAdam newActivations newMActivations newVActivations = descendMatrix opt (MatrixValuesAdam (getStep store) oldActivations activationGradient oldMActivations oldVActivations)
        newStore = setData opt x store [FullyConnected' newMBias newMActivations, FullyConnected' newVBias newVActivations]
    in FullyConnected (FullyConnected' newBias newActivations) newStore

  reduceGradient grads = FullyConnected' (dvmap (/l) bs) (dmmap (/l) as)
    where
      FullyConnected' bs as = foldl1' (\(FullyConnected' bs as) (FullyConnected' bs' as') -> FullyConnected' (bs + bs') (as + as')) grads
      l = fromIntegral $ length grads :: RealNum

instance (KnownNat i, KnownNat o, KnownNat (i * o)) => LayerOptimizerData (FullyConnected i o) (Optimizer 'SGD) where
  type MomentumDataType (FullyConnected i o) (Optimizer 'SGD) = FullyConnected' i o
  getData opt x store = head $ getListStore opt x store
  setData opt x store = setListStore opt x store . return
  newData _ _ = FullyConnected' (konst 0) (konst 0)

instance (KnownNat i, KnownNat o, KnownNat (i * o)) => LayerOptimizerData (FullyConnected i o) (Optimizer 'Adam) where
  type MomentumDataType (FullyConnected i o) (Optimizer 'Adam) = FullyConnected' i o
  type MomentumExpOptResult (FullyConnected i o) (Optimizer 'Adam) = [FullyConnected' i o]
  getData = getListStore
  setData = setListStore
  newData _ _ = FullyConnected' (konst 0) (konst 0)

instance (KnownNat i, KnownNat o) => FoldableGradient (FullyConnected' i o) where
  mapGradient f (FullyConnected' bias activations) = FullyConnected' (dvmap f bias) (dmmap f activations)
  squaredSums (FullyConnected' bias activations) = [sumV . squareV $ bias, sumM . squareM $ activations]

instance (KnownNat i, KnownNat o, KnownNat (i * o)) => Layer (FullyConnected i o) ('D1 i) ('D1 o) where
  type Tape (FullyConnected i o) ('D1 i) ('D1 o) = R i
  -- Do a matrix vector multiplication and return the result.
  runForwards (FullyConnected (FullyConnected' wB wN) _) (S1D v) = (v, S1D (wB + wN #> v))

  -- Run a backpropogation step for a full connected layer.
  runBackwards (FullyConnected (FullyConnected' _ wN) _) x (S1D dEdy) =
          let wB'  = dEdy
              mm'  = dEdy `outer` x
              -- calcluate derivatives for next step
              dWs  = tr wN #> dEdy
          in  (FullyConnected' wB' mm', S1D dWs)

instance (KnownNat i, KnownNat o) => Serialize (FullyConnected i o) where
  put (FullyConnected w ms) = put w >> put ms
  get = FullyConnected <$> get <*> get

instance (KnownNat i, KnownNat o) => Serialize (FullyConnected' i o) where
  put (FullyConnected' b w) = do
    putListOf put . LA.toList . extract $ b
    putListOf put . LA.toList . LA.flatten . extract $ w
  get = do
      let f  = fromIntegral $ natVal (Proxy :: Proxy i)
      b     <- maybe (fail "Vector of incorrect size") return . create . LA.fromList =<< getListOf get
      k     <- maybe (fail "Vector of incorrect size") return . create . LA.reshape f . LA.fromList =<< getListOf get
      return $ FullyConnected' b k

instance (KnownNat i, KnownNat o, KnownNat (i*o)) => RandomLayer (FullyConnected i o) where
  createRandomWith = randomFullyConnected

randomFullyConnected :: forall m i o . (PrimBase m, KnownNat i, KnownNat o, KnownNat (i*o))
                     => WeightInitMethod -> Gen (PrimState m) -> m (FullyConnected i o)
randomFullyConnected m gen = do
  wN <- getRandomMatrix i o m gen
  wB <- getRandomVector i o m gen
  return $ FullyConnected (FullyConnected' wB wN) mkListStore
  where i = natVal (Proxy :: Proxy i)
        o = natVal (Proxy :: Proxy o)

instance OnnxOperator (FullyConnected i o) where
  onnxOpTypeNames _ = ["Gemm"]

instance (KnownNat i, KnownNat o) => OnnxLoadable (FullyConnected i o) where
  loadOnnxNode inits node = case (node ^. #input) of
    [_, b, c] -> do
      -- FIXME: Proper attribute checking
      -- node `doesNotHaveAttribute` "alpha"
      -- node `doesNotHaveAttribute` "transA"
      -- node `doesNotHaveAttribute` "transB"

      let beta = fromRight 1 (readFloatAttributeToRealNum "beta" node)
      loadedB <- readInitializerMatrix inits b
      loadedC <- readInitializerVector inits c

      return $ FullyConnected (FullyConnected' loadedC (dmmap (*beta) loadedB)) mkListStore
    _         -> onnxIncorrectNumberOfInputs

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-# LANGUAGE UndecidableInstances,
             FlexibleInstances,
             FlexibleContexts,
             TypeFamilies,
             ScopedTypeVariables #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Signal.Multichannel
-- Copyright   :  (c) Alexander Vivian Hugh McPhail 2010, 2014, 2015, 2016
-- License     :  BSD3
--
-- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
-- Stability   :  provisional
-- Portability :  uses Concurrency
--
-- Signal processing functions, multichannel datatype
--
-- link with '-threaded' and run with +RTS -Nn, where n is the number of CPUs
--
-----------------------------------------------------------------------------

--             IncoherentInstances,

module Numeric.Signal.Multichannel (
                       Multichannel,readMultichannel,writeMultichannel,
                       createMultichannel,
                       sampling_rate,precision,channels,samples,
                       detrended,filtered,
                       getChannel,getChannels,
                       toMatrix,
                       mapConcurrently,
                       detrend,filter,
                       slice,
                       histograms,
                       entropy_delta_phase,mi_phase
                ) where

-----------------------------------------------------------------------------

import qualified Numeric.Signal as S

--import Complex

import qualified Data.Array.IArray as I
import Data.Ix

--import Data.Word

import Control.Concurrent
--import Control.Concurrent.MVar

import System.IO.Unsafe(unsafePerformIO)

--import qualified Data.List as L

import Data.Binary
import Data.Maybe

import Foreign.Storable

import Numeric.LinearAlgebra hiding(range)

import qualified Data.Vector.Generic as GV

--import qualified Numeric.GSL.Fourier as F

import qualified Numeric.GSL.Histogram as H
import qualified Numeric.GSL.Histogram2D as H2

import qualified Numeric.Statistics.Information as SI

import Prelude hiding(filter)

import Control.Monad(replicateM)

{-
-------------------------------------------------------------------

instance (Binary a, Storable a) => Binary (Vector a) where
    put v = do
            let d = GV.length v
            put d
            mapM_ (\i -> put $ v @> i) [0..(d-1)]
    get = do
          d <- get
          xs <- replicateM d get
          return $ fromList xs

-------------------------------------------------------------------
-}

-----------------------------------------------------------------------------

-- | data type with multiple channels
data Multichannel a = MC {
                          _sampling_rate :: Int             -- ^ sampling rate
                          , _precision   :: Int             -- ^ bits of precision
                          , _channels    :: Int             -- ^ number of channels
                          , _length      :: Int             -- ^ length in samples
                          , _detrended   :: Bool            -- ^ was the data detrended?
                          , _filtered    :: Maybe (Int,Int) -- ^ if filtered the passband
                          , _data        :: I.Array Int (Vector a) -- ^ data
                         }

-----------------------------------------------------------------------------

instance Binary (Multichannel Double) where
    put (MC s p c l de f d) = do
                              put s
                              put p
                              put c
                              put l
                              put de
                              put f
                              put $! fmap convert d
        where convert v = let (mi,ma) = (minElement v,maxElement v)
                              v' = GV.map (\x -> round $ (x - mi)/(ma - mi) * (fromIntegral (maxBound :: Word64))) v
                          in (mi,ma,v' :: Vector Word64)

    get = do
          s <- get
          p <- get
          c <- get
          l <- get
          de <- get
          f <- get
          (d :: I.Array Int (Double,Double,Vector Word64)) <- get
          return $! (MC s p c l de f (seq d (fmap convert) d))
              where convert (mi,ma,v) = GV.map (\x -> ((fromIntegral x)) / (fromIntegral (maxBound :: Word64)) * (ma - mi) + mi) v

instance Binary (Multichannel Float) where
    put (MC s p c l de f d) = do
                              put s
                              put p
                              put c
                              put l
                              put de
                              put f
                              put $! fmap convert d
        where convert v = let (mi,ma) = (minElement v,maxElement v)
                              v' = GV.map (\x -> round $ (x - mi)/(ma - mi) * (fromIntegral (maxBound :: Word64))) v
                          in (mi,ma,v' :: Vector Word64)

    get = do
          s <- get
          p <- get
          c <- get
          l <- get
          de <- get
          f <- get
          (d :: I.Array Int (Float,Float,Vector Word32)) <- get
          return $! (MC s p c l de f (seq d (fmap convert) d))
              where convert (mi,ma,v) = GV.map (\x -> ((fromIntegral x)) / (fromIntegral (maxBound :: Word32)) * (ma - mi) + mi) v

instance Binary (Multichannel (Complex Double)) where
    put (MC s p c l de f d) = do
                              put s
                              put p
                              put c
                              put l
                              put de
                              put f
                              put $! fmap ((\(r,j) -> (convert r, convert j)) . fromComplex) d
        where convert v = let (mi,ma) = (minElement v,maxElement v)
                              v' = GV.map (\x -> round $ (x - mi)/(ma - mi) * (fromIntegral (maxBound :: Word64))) v
                          in (mi,ma,v' :: Vector Word64)

    get = do
          s <- get
          p <- get
          c <- get
          l <- get
          de <- get
          f <- get
          (d :: I.Array Int ((Double,Double,Vector Word64),(Double,Double,Vector Word64))) <- get
          return $! (MC s p c l de f (seq d (fmap (\(r,j) -> toComplex (convert r,convert j)) d)))
              where convert (mi,ma,v) = GV.map (\x -> ((fromIntegral x)) / (fromIntegral (maxBound :: Word64)) * (ma - mi) + mi) v

instance Binary (Multichannel (Complex Float)) where
    put (MC s p c l de f d) = do
                              put s
                              put p
                              put c
                              put l
                              put de
                              put f
                              put $! fmap ((\(r,j) -> (convert r, convert j)) . fromComplex) d
        where convert v = let (mi,ma) = (minElement v,maxElement v)
                              v' = GV.map (\x -> round $ (x - mi)/(ma - mi) * (fromIntegral (maxBound :: Word32))) v
                          in (mi,ma,v' :: Vector Word32)

    get = do
          s <- get
          p <- get
          c <- get
          l <- get
          de <- get
          f <- get
          (d :: I.Array Int ((Float,Float,Vector Word32),(Float,Float,Vector Word32))) <- get
          return $! (MC s p c l de f (seq d (fmap (\(r,j) -> toComplex (convert r,convert j)) d)))
              where convert (mi,ma,v) = GV.map (\x -> ((fromIntegral x)) / (fromIntegral (maxBound :: Word32)) * (ma - mi) + mi) v

-----------------------------------------------------------------------------

readMultichannel :: (Binary (Multichannel a)) => FilePath -> IO (Multichannel a)
readMultichannel = decodeFile

writeMultichannel :: (Binary (Multichannel a)) => FilePath -> Multichannel a -> IO ()
writeMultichannel = encodeFile

-----------------------------------------------------------------------------

-- | create a multichannel data type
createMultichannel :: Storable a
                   => Int               -- ^ sampling rate
                   -> Int               -- ^ bits of precision
                   -> [Vector a]        -- ^ data
                   -> Multichannel a    -- ^ datatype
createMultichannel s p d = let c = length d
                 in MC s p c (GV.length $ head d) False Nothing (I.listArray (1,c) d)

-- | the sampling rate
sampling_rate :: Multichannel a -> Int
sampling_rate = _sampling_rate

-- | the bits of precision
precision :: Multichannel a -> Int
precision = _precision

-- | the number of channels
channels :: Multichannel a -> Int
channels = _channels

-- | the length, in samples
samples :: Multichannel a -> Int
samples = _length

-- | extract one channel
getChannel :: Int -> Multichannel a -> Vector a
getChannel c d = (_data d) I.! c

-- | extract all channels
getChannels :: Multichannel a -> I.Array Int (Vector a)
getChannels d = _data d

-- | convert the data to a matrix with channels as rows
toMatrix :: Element a => Multichannel a -> Matrix a
toMatrix = fromRows . I.elems . _data

-- | was the data detrended?
detrended :: Multichannel a -> Bool
detrended = _detrended

-- | was the data filtered?
filtered :: Multichannel a -> Maybe (Int,Int)
filtered = _filtered

-----------------------------------------------------------------------------

-- | map a function executed concurrently
mapArrayConcurrently :: Ix i => (a -> b)    -- ^ function to map
                     -> I.Array i a         -- ^ input
                     -> (I.Array i b)     -- ^ output
mapArrayConcurrently f d = unsafePerformIO $ do
  let b = I.bounds d
  results <- replicateM (rangeSize b) newEmptyMVar
  mapM_ (forkIO . applyFunction f) $ zip results (I.assocs d)
  vectors <- mapM takeMVar results
  return $ I.array b vectors
    where applyFunction f' (m,(j,e)) = putMVar m (j,f' e)

{-
-- | map a function executed concurrently
mapListConcurrently :: (a -> b)             -- ^ function to map
                    -> [a]                  -- ^ input
                    -> [b]                  -- ^ output
mapListConcurrently f d = unsafePerformIO $ do
                                            results <- replicateM (length d) newEmptyMVar
                                            mapM_ (forkIO . applyFunction f) zip results d
                                            mapM takeMVar results
    where applyFunction f' (m,e) = putMVar m (f' e)
-}

-- | map a function executed concurrently
mapConcurrently :: Storable b
                => (Vector a -> Vector b)      -- ^ the function to be mapped
                -> Multichannel a              -- ^ input data
                -> Multichannel b              -- ^ output data
mapConcurrently f (MC sr p c _ de fi d) = let d' = mapArrayConcurrently f d
                                          in MC sr p c (GV.length $ d' I.! 1) de fi d'

-- | map a function
mapMC :: Storable b
      => (Vector a -> Vector b)                -- ^ the function to be mapped
      -> Multichannel a                        -- ^ input data
      -> Multichannel b                        -- ^ output data
mapMC f (MC sr p c _ de fi d) = let d' = fmap f d
                                in MC sr p c (GV.length $ d' I.! 1) de fi d'

-----------------------------------------------------------------------------

-- | detrend the data with a specified window size
detrend :: Int -> Multichannel Double -> Multichannel Double
detrend w m = let m' = mapConcurrently (S.detrend w) m
              in m' { _detrended = True }

-- | filter the data with the given passband
filter :: (S.Filterable a, Double ~ DoubleOf a) =>
         (Int,Int) -> Multichannel a -> Multichannel a
filter pb m = let m' = mapConcurrently (S.broadband_filter (_sampling_rate m) pb) m
              in m' { _filtered = Just pb }

-----------------------------------------------------------------------------

-- | extract a slice of the data
slice :: Storable a
      => Int                 -- ^ starting sample number
      -> Int                 -- ^ length
      -> Multichannel a
      -> Multichannel a
slice j w m = let m' = mapConcurrently (subVector j w) m
              in m' { _length = w }

-----------------------------------------------------------------------------

-- | calculate histograms
histograms :: (S.Filterable a, Double ~ DoubleOf a) =>
            I.Array Int (Vector a)
          -> Int -> (Double,Double)
          -> Int -> Int -> (Double,Double) -> (Double,Double) -- ^ bins and ranges
          -> (I.Array Int H.Histogram,I.Array (Int,Int) H2.Histogram2D)
histograms d' b (l,u) bx by (lx,ux) (ly,uy)
  = let d = fmap double d'
        (bl,bu) = I.bounds d
        br = ((bl,bl),(bu,bu))
        histarray = mapArrayConcurrently (H.fromLimits b (l,u)) d
        pairs = I.array br $ map (\(m,n) -> ((m,n),(d I.! m,d I.! n))) (range br)
        hist2array = mapArrayConcurrently (\(x,y) -> (H2.addVector (H2.emptyLimits bx by (lx,ux) (ly,uy)) x y)) pairs
    in (histarray,hist2array)

-----------------------------------------------------------------------------

-- | calculate the entropy of the phase difference between pairs of channels (fills upper half of matrix)
entropy_delta_phase :: (S.Filterable a, Double ~ DoubleOf a) =>
                Multichannel a      -- ^ input data
              -> Matrix Double
entropy_delta_phase m = let d = _data m
                            c = _channels m
                            b = ((1,1),(c,c))
                            r = I.range b
                            diff = I.listArray b (map (\j@(x,y) -> (j,if x <= y then Just (double $ (d I.! y)-(d I.! x)) else Nothing)) r) :: I.Array (Int,Int) ((Int,Int),Maybe (Vector Double))
                            h = mapArrayConcurrently (maybe Nothing (\di -> Just $ H.fromLimits 128 ((-2)*pi,2*pi) di)) (fmap snd diff)
                            ent = mapArrayConcurrently (\(j,difvec) -> case difvec of
                                                        Nothing -> 0 :: Double
                                                        Just da -> SI.entropy (fromJust (h I.! j)) da) diff
                        in fromArray2D ent

-----------------------------------------------------------------------------

-- | calculate the mutual information of the phase between pairs of channels (fills upper half of matrix)
mi_phase :: (S.Filterable a, Double ~ DoubleOf a) =>
           Multichannel a      -- ^ input data
         -> Matrix Double
mi_phase m = let d = _data m
                 (histarray,hist2array) = histograms d 128 (-pi,pi) 128 128 (-pi,pi) (-pi,pi)
                 indhist = I.listArray (I.bounds hist2array) (I.assocs hist2array)
                 mi = mapArrayConcurrently (doMI histarray (fmap double d)) indhist
             in fromArray2D mi
    where doMI histarray d ((x,y),h2)
              | x <= y     = SI.mutual_information h2 (histarray I.! x) (histarray I.! y) (d I.! x,d I.! y)
              | otherwise = 0

-----------------------------------------------------------------------------

max_stars_repo_name: amcphail/hsignal

max_issues_repo_name: amcphail/hsignal

max_forks_repo_name: amcphail/hsignal