text:
{-# LANGUAGE BangPatterns, CPP, MultiParamTypeClasses, TypeFamilies, FlexibleContexts #-}
#if __GLASGOW_HASKELL__ >= 707
{-# LANGUAGE DeriveDataTypeable, StandaloneDeriving #-}
#endif
{-# OPTIONS_HADDOCK hide #-}

-- |
-- Module      : Data.Vector.Unboxed.Base
-- Copyright   : (c) Roman Leshchinskiy 2009-2010
-- License     : BSD-style
--
-- Maintainer  : Roman Leshchinskiy <rl@cse.unsw.edu.au>
-- Stability   : experimental
-- Portability : non-portable
--
-- Adaptive unboxed vectors: basic implementation
--

module Data.Vector.Unboxed.Base (
  MVector(..), IOVector, STVector, Vector(..), Unbox
) where

import qualified Data.Vector.Generic         as G
import qualified Data.Vector.Generic.Mutable as M

import qualified Data.Vector.Primitive as P

import Control.DeepSeq ( NFData(rnf) )

import Control.Monad.Primitive
import Control.Monad ( liftM )

import Data.Word ( Word8, Word16, Word32, Word64 )
import Data.Int  ( Int8, Int16, Int32, Int64 )
import Data.Complex

#if !MIN_VERSION_base(4,8,0)
import Data.Word ( Word )
#endif

#if __GLASGOW_HASKELL__ >= 707
import Data.Typeable ( Typeable )
#else
import Data.Typeable ( Typeable1(..), Typeable2(..), mkTyConApp,
#if MIN_VERSION_base(4,4,0)
                       mkTyCon3
#else
                       mkTyCon
#endif
                     )
#endif

import Data.Data     ( Data(..) )

-- Data.Vector.Internal.Check is unused
#define NOT_VECTOR_MODULE
#include "vector.h"

data family MVector s a
data family Vector    a

type IOVector = MVector RealWorld
type STVector s = MVector s

type instance G.Mutable Vector = MVector

class (G.Vector Vector a, M.MVector MVector a) => Unbox a

instance NFData (Vector a) where rnf !_ = ()
instance NFData (MVector s a) where rnf !_ = ()

-- -----------------
-- Data and Typeable
-- -----------------
#if __GLASGOW_HASKELL__ >= 707
deriving instance Typeable Vector
deriving instance Typeable MVector
#else
#if MIN_VERSION_base(4,4,0)
vectorTyCon = mkTyCon3 "vector"
#else
vectorTyCon m s = mkTyCon $ m ++ "." ++ s
#endif

instance Typeable1 Vector where
  typeOf1 _ = mkTyConApp (vectorTyCon "Data.Vector.Unboxed" "Vector") []

instance Typeable2 MVector where
  typeOf2 _ = mkTyConApp (vectorTyCon "Data.Vector.Unboxed.Mutable" "MVector") []
#endif

instance (Data a, Unbox a) => Data (Vector a) where
  gfoldl       = G.gfoldl
  toConstr _   = error "toConstr"
  gunfold _ _  = error "gunfold"
  dataTypeOf _ = G.mkType "Data.Vector.Unboxed.Vector"
  dataCast1    = G.dataCast

-- ----
-- Unit
-- ----

newtype instance MVector s () = MV_Unit Int
newtype instance Vector    () = V_Unit Int

instance Unbox ()

instance M.MVector MVector () where
  {-# INLINE basicLength #-}
  {-# INLINE basicUnsafeSlice #-}
  {-# INLINE basicOverlaps #-}
  {-# INLINE basicUnsafeNew #-}
  {-# INLINE basicInitialize #-}
  {-# INLINE basicUnsafeRead #-}
  {-# INLINE basicUnsafeWrite #-}
  {-# INLINE basicClear #-}
  {-# INLINE basicSet #-}
  {-# INLINE basicUnsafeCopy #-}
  {-# INLINE basicUnsafeGrow #-}

  basicLength (MV_Unit n) = n

  basicUnsafeSlice _ m (MV_Unit _) = MV_Unit m

  basicOverlaps _ _ = False

  basicUnsafeNew n = return (MV_Unit n)

  -- Nothing to initialize
  basicInitialize _ = return ()

  basicUnsafeRead (MV_Unit _) _ = return ()

  basicUnsafeWrite (MV_Unit _) _ () = return ()

  basicClear _ = return ()

  basicSet (MV_Unit _) () = return ()

  basicUnsafeCopy (MV_Unit _) (MV_Unit _) = return ()

  basicUnsafeGrow (MV_Unit n) m = return $ MV_Unit (n+m)

instance G.Vector Vector () where
  {-# INLINE basicUnsafeFreeze #-}
  basicUnsafeFreeze (MV_Unit n) = return $ V_Unit n

  {-# INLINE basicUnsafeThaw #-}
  basicUnsafeThaw (V_Unit n) = return $ MV_Unit n

  {-# INLINE basicLength #-}
  basicLength (V_Unit n) = n

  {-# INLINE basicUnsafeSlice #-}
  basicUnsafeSlice _ m (V_Unit _) = V_Unit m

  {-# INLINE basicUnsafeIndexM #-}
  basicUnsafeIndexM (V_Unit _) _ = return ()

  {-# INLINE basicUnsafeCopy #-}
  basicUnsafeCopy (MV_Unit _) (V_Unit _) = return ()

  {-# INLINE elemseq #-}
  elemseq _ = seq

-- ---------------
-- Primitive types
-- ---------------

#define primMVector(ty,con)                                             \
instance M.MVector MVector ty where {                                   \
  {-# INLINE basicLength #-}                                            \
; {-# INLINE basicUnsafeSlice #-}                                       \
; {-# INLINE basicOverlaps #-}                                          \
; {-# INLINE basicUnsafeNew #-}                                         \
; {-# INLINE basicInitialize #-}                                        \
; {-# INLINE basicUnsafeReplicate #-}                                   \
; {-# INLINE basicUnsafeRead #-}                                        \
; {-# INLINE basicUnsafeWrite #-}                                       \
; {-# INLINE basicClear #-}                                             \
; {-# INLINE basicSet #-}                                               \
; {-# INLINE basicUnsafeCopy #-}                                        \
; {-# INLINE basicUnsafeGrow #-}                                        \
; basicLength (con v) = M.basicLength v                                 \
; basicUnsafeSlice i n (con v) = con $ M.basicUnsafeSlice i n v         \
; basicOverlaps (con v1) (con v2) = M.basicOverlaps v1 v2               \
; basicUnsafeNew n = con `liftM` M.basicUnsafeNew n                     \
; basicInitialize (con v) = M.basicInitialize v                         \
; basicUnsafeReplicate n x = con `liftM` M.basicUnsafeReplicate n x     \
; basicUnsafeRead (con v) i = M.basicUnsafeRead v i                     \
; basicUnsafeWrite (con v) i x = M.basicUnsafeWrite v i x               \
; basicClear (con v) = M.basicClear v                                   \
; basicSet (con v) x = M.basicSet v x                                   \
; basicUnsafeCopy (con v1) (con v2) = M.basicUnsafeCopy v1 v2           \
; basicUnsafeMove (con v1) (con v2) = M.basicUnsafeMove v1 v2           \
; basicUnsafeGrow (con v) n = con `liftM` M.basicUnsafeGrow v n }

#define primVector(ty,con,mcon)                                         \
instance G.Vector Vector ty where {                                     \
  {-# INLINE basicUnsafeFreeze #-}                                      \
; {-# INLINE basicUnsafeThaw #-}                                        \
; {-# INLINE basicLength #-}                                            \
; {-# INLINE basicUnsafeSlice #-}                                       \
; {-# INLINE basicUnsafeIndexM #-}                                      \
; {-# INLINE elemseq #-}                                                \
; basicUnsafeFreeze (mcon v) = con `liftM` G.basicUnsafeFreeze v        \
; basicUnsafeThaw (con v) = mcon `liftM` G.basicUnsafeThaw v            \
; basicLength (con v) = G.basicLength v                                 \
; basicUnsafeSlice i n (con v) = con $ G.basicUnsafeSlice i n v         \
; basicUnsafeIndexM (con v) i = G.basicUnsafeIndexM v i                 \
; basicUnsafeCopy (mcon mv) (con v) = G.basicUnsafeCopy mv v            \
; elemseq _ = seq }

newtype instance MVector s Int = MV_Int (P.MVector s Int)
newtype instance Vector    Int = V_Int  (P.Vector    Int)
instance Unbox Int
primMVector(Int, MV_Int)
primVector(Int, V_Int, MV_Int)

newtype instance MVector s Int8 = MV_Int8 (P.MVector s Int8)
newtype instance Vector    Int8 = V_Int8  (P.Vector    Int8)
instance Unbox Int8
primMVector(Int8, MV_Int8)
primVector(Int8, V_Int8, MV_Int8)

newtype instance MVector s Int16 = MV_Int16 (P.MVector s Int16)
newtype instance Vector    Int16 = V_Int16  (P.Vector    Int16)
instance Unbox Int16
primMVector(Int16, MV_Int16)
primVector(Int16, V_Int16, MV_Int16)

newtype instance MVector s Int32 = MV_Int32 (P.MVector s Int32)
newtype instance Vector    Int32 = V_Int32  (P.Vector    Int32)
instance Unbox Int32
primMVector(Int32, MV_Int32)
primVector(Int32, V_Int32, MV_Int32)

newtype instance MVector s Int64 = MV_Int64 (P.MVector s Int64)
newtype instance Vector    Int64 = V_Int64  (P.Vector    Int64)
instance Unbox Int64
primMVector(Int64, MV_Int64)
primVector(Int64, V_Int64, MV_Int64)

newtype instance MVector s Word = MV_Word (P.MVector s Word)
newtype instance Vector    Word = V_Word  (P.Vector    Word)
instance Unbox Word
primMVector(Word, MV_Word)
primVector(Word, V_Word, MV_Word)

newtype instance MVector s Word8 = MV_Word8 (P.MVector s Word8)
newtype instance Vector    Word8 = V_Word8  (P.Vector    Word8)
instance Unbox Word8
primMVector(Word8, MV_Word8)
primVector(Word8, V_Word8, MV_Word8)

newtype instance MVector s Word16 = MV_Word16 (P.MVector s Word16)
newtype instance Vector    Word16 = V_Word16  (P.Vector    Word16)
instance Unbox Word16
primMVector(Word16, MV_Word16)
primVector(Word16, V_Word16, MV_Word16)

newtype instance MVector s Word32 = MV_Word32 (P.MVector s Word32)
newtype instance Vector    Word32 = V_Word32  (P.Vector    Word32)
instance Unbox Word32
primMVector(Word32, MV_Word32)
primVector(Word32, V_Word32, MV_Word32)

newtype instance MVector s Word64 = MV_Word64 (P.MVector s Word64)
newtype instance Vector    Word64 = V_Word64  (P.Vector    Word64)
instance Unbox Word64
primMVector(Word64, MV_Word64)
primVector(Word64, V_Word64, MV_Word64)

newtype instance MVector s Float = MV_Float (P.MVector s Float)
newtype instance Vector    Float = V_Float  (P.Vector    Float)
instance Unbox Float
primMVector(Float, MV_Float)
primVector(Float, V_Float, MV_Float)

newtype instance MVector s Double = MV_Double (P.MVector s Double)
newtype instance Vector    Double = V_Double  (P.Vector    Double)
instance Unbox Double
primMVector(Double, MV_Double)
primVector(Double, V_Double, MV_Double)

newtype instance MVector s Char = MV_Char (P.MVector s Char)
newtype instance Vector    Char = V_Char  (P.Vector    Char)
instance Unbox Char
primMVector(Char, MV_Char)
primVector(Char, V_Char, MV_Char)

-- ----
-- Bool
-- ----

fromBool :: Bool -> Word8
{-# INLINE fromBool #-}
fromBool True = 1
fromBool False = 0

toBool :: Word8 -> Bool
{-# INLINE toBool #-}
toBool 0 = False
toBool _ = True

newtype instance MVector s Bool = MV_Bool (P.MVector s Word8)
newtype instance Vector    Bool = V_Bool  (P.Vector    Word8)

instance Unbox Bool

instance M.MVector MVector Bool where
  {-# INLINE basicLength #-}
  {-# INLINE basicUnsafeSlice #-}
  {-# INLINE basicOverlaps #-}
  {-# INLINE basicUnsafeNew #-}
  {-# INLINE basicInitialize #-}
  {-# INLINE basicUnsafeReplicate #-}
  {-# INLINE basicUnsafeRead #-}
  {-# INLINE basicUnsafeWrite #-}
  {-# INLINE basicClear #-}
  {-# INLINE basicSet #-}
  {-# INLINE basicUnsafeCopy #-}
  {-# INLINE basicUnsafeGrow #-}
  basicLength (MV_Bool v) = M.basicLength v
  basicUnsafeSlice i n (MV_Bool v) = MV_Bool $ M.basicUnsafeSlice i n v
  basicOverlaps (MV_Bool v1) (MV_Bool v2) = M.basicOverlaps v1 v2
  basicUnsafeNew n = MV_Bool `liftM` M.basicUnsafeNew n
  basicInitialize (MV_Bool v) = M.basicInitialize v
  basicUnsafeReplicate n x = MV_Bool `liftM` M.basicUnsafeReplicate n (fromBool x)
  basicUnsafeRead (MV_Bool v) i = toBool `liftM` M.basicUnsafeRead v i
  basicUnsafeWrite (MV_Bool v) i x = M.basicUnsafeWrite v i (fromBool x)
  basicClear (MV_Bool v) = M.basicClear v
  basicSet (MV_Bool v) x = M.basicSet v (fromBool x)
  basicUnsafeCopy (MV_Bool v1) (MV_Bool v2) = M.basicUnsafeCopy v1 v2
  basicUnsafeMove (MV_Bool v1) (MV_Bool v2) = M.basicUnsafeMove v1 v2
  basicUnsafeGrow (MV_Bool v) n = MV_Bool `liftM` M.basicUnsafeGrow v n

instance G.Vector Vector Bool where
  {-# INLINE basicUnsafeFreeze #-}
  {-# INLINE basicUnsafeThaw #-}
  {-# INLINE basicLength #-}
  {-# INLINE basicUnsafeSlice #-}
  {-# INLINE basicUnsafeIndexM #-}
  {-# INLINE elemseq #-}
  basicUnsafeFreeze (MV_Bool v) = V_Bool `liftM` G.basicUnsafeFreeze v
  basicUnsafeThaw (V_Bool v) = MV_Bool `liftM` G.basicUnsafeThaw v
  basicLength (V_Bool v) = G.basicLength v
  basicUnsafeSlice i n (V_Bool v) = V_Bool $ G.basicUnsafeSlice i n v
  basicUnsafeIndexM (V_Bool v) i = toBool `liftM` G.basicUnsafeIndexM v i
  basicUnsafeCopy (MV_Bool mv) (V_Bool v) = G.basicUnsafeCopy mv v
  elemseq _ = seq

-- -------
-- Complex
-- -------

newtype instance MVector s (Complex a) = MV_Complex (MVector s (a,a))
newtype instance Vector    (Complex a) = V_Complex  (Vector    (a,a))

instance (RealFloat a, Unbox a) => Unbox (Complex a)

instance (RealFloat a, Unbox a) => M.MVector MVector (Complex a) where
  {-# INLINE basicLength #-}
  {-# INLINE basicUnsafeSlice #-}
  {-# INLINE basicOverlaps #-}
  {-# INLINE basicUnsafeNew #-}
  {-# INLINE basicInitialize #-}
  {-# INLINE basicUnsafeReplicate #-}
  {-# INLINE basicUnsafeRead #-}
  {-# INLINE basicUnsafeWrite #-}
  {-# INLINE basicClear #-}
  {-# INLINE basicSet #-}
  {-# INLINE basicUnsafeCopy #-}
  {-# INLINE basicUnsafeGrow #-}
  basicLength (MV_Complex v) = M.basicLength v
  basicUnsafeSlice i n (MV_Complex v) = MV_Complex $ M.basicUnsafeSlice i n v
  basicOverlaps (MV_Complex v1) (MV_Complex v2) = M.basicOverlaps v1 v2
  basicUnsafeNew n = MV_Complex `liftM` M.basicUnsafeNew n
  basicInitialize (MV_Complex v) = M.basicInitialize v
  basicUnsafeReplicate n (x :+ y) = MV_Complex `liftM` M.basicUnsafeReplicate n (x,y)
  basicUnsafeRead (MV_Complex v) i = uncurry (:+) `liftM` M.basicUnsafeRead v i
  basicUnsafeWrite (MV_Complex v) i (x :+ y) = M.basicUnsafeWrite v i (x,y)
  basicClear (MV_Complex v) = M.basicClear v
  basicSet (MV_Complex v) (x :+ y) = M.basicSet v (x,y)
  basicUnsafeCopy (MV_Complex v1) (MV_Complex v2) = M.basicUnsafeCopy v1 v2
  basicUnsafeMove (MV_Complex v1) (MV_Complex v2) = M.basicUnsafeMove v1 v2
  basicUnsafeGrow (MV_Complex v) n = MV_Complex `liftM` M.basicUnsafeGrow v n

instance (RealFloat a, Unbox a) => G.Vector Vector (Complex a) where
  {-# INLINE basicUnsafeFreeze #-}
  {-# INLINE basicUnsafeThaw #-}
  {-# INLINE basicLength #-}
  {-# INLINE basicUnsafeSlice #-}
  {-# INLINE basicUnsafeIndexM #-}
  {-# INLINE elemseq #-}
  basicUnsafeFreeze (MV_Complex v) = V_Complex `liftM` G.basicUnsafeFreeze v
  basicUnsafeThaw (V_Complex v) = MV_Complex `liftM` G.basicUnsafeThaw v
  basicLength (V_Complex v) = G.basicLength v
  basicUnsafeSlice i n (V_Complex v) = V_Complex $ G.basicUnsafeSlice i n v
  basicUnsafeIndexM (V_Complex v) i
                = uncurry (:+) `liftM` G.basicUnsafeIndexM v i
  basicUnsafeCopy (MV_Complex mv) (V_Complex v)
                = G.basicUnsafeCopy mv v
  elemseq _ (x :+ y) z = G.elemseq (undefined :: Vector a) x
                       $ G.elemseq (undefined :: Vector a) y z

-- ------
-- Tuples
-- ------

#define DEFINE_INSTANCES
#include "unbox-tuple-instances"


text:
{-# LANGUAGE NoMonomorphismRestriction, Rank2Types #-}
module STVector (
    tests_STVector,
    mutatesToVector,
    readOnlyVector,
    ) where

import Control.Monad
import Control.Monad.ST
import Data.AEq
import Data.Complex( magnitude )
import Debug.Trace
import Foreign( Storable )
import Test.Framework
import Test.Framework.Providers.QuickCheck2
import Test.QuickCheck hiding ( vector )
import qualified Test.QuickCheck as QC

import Numeric.LinearAlgebra.Vector( Vector, STVector )
import qualified Numeric.LinearAlgebra.Vector as V

import Test.QuickCheck.LinearAlgebra( TestElem(..), Dim(..), Index(..),
    Assocs(..), VectorPair(..), VectorTriple(..) )
import qualified Test.QuickCheck.LinearAlgebra as Test
import Test.QuickCheck.Property(rejected)

import Typed

--- this could be a bad idea ..
instance Testable () where
  property _ = property rejected

tests_STVector = testGroup "STVector"
    [ testPropertyI "new_" prop_new_
    , testPropertyI "new" prop_new
    , testPropertyI "newCopy" prop_newCopy
    , testPropertyI "copyTo" prop_copyTo
    , testPropertyDZ "swap" prop_swap prop_swap
    , testPropertyI "read" prop_read
    , testPropertyI "write" prop_write
    , testPropertyI "modify" prop_modify
    , testPropertyI "getElems" prop_getElems
    , testPropertyI "getElems'" prop_getElems'
    , testPropertyI "getAssocs" prop_getAssocs
    , testPropertyI "getAssocs'" prop_getAssocs'
    , testPropertyI "setElems" prop_setElems
    , testPropertyI "setAssocs" prop_setAssocs
    , testPropertyI "mapTo" prop_mapTo
    , testPropertyI "zipWithTo" prop_zipWithTo
    , testPropertyDZ "getSumAbs" prop_getSumAbs prop_getSumAbs
    , testPropertyDZ "getNorm2" prop_getNorm2 prop_getNorm2
    , testPropertyDZ "getDot" prop_getDot prop_getDot
    , testPropertyDZ "kroneckerTo" prop_kroneckerTo prop_kroneckerTo
    , testPropertyDZ "addTo" prop_addTo prop_addTo
    , testPropertyDZ "subTo" prop_subTo prop_subTo
    , testPropertyDZ "mulTo" prop_mulTo prop_mulTo
    , testPropertyDZ "negateTo" prop_negateTo prop_negateTo
    , testPropertyDZ "conjugateTo" prop_conjugateTo prop_conjugateTo
    , testPropertyDZ "absTo" prop_absTo prop_absTo
    , testPropertyDZ "signumTo" prop_signumTo prop_signumTo
    , testPropertyDZ "divTo" prop_divTo prop_divTo
    , testPropertyDZ "recipTo" prop_recipTo prop_recipTo
    , testPropertyDZ "sqrtTo" prop_sqrtTo prop_sqrtTo
    , testPropertyDZ "expTo" prop_expTo prop_expTo
    , testPropertyDZ "logTo" prop_logTo prop_logTo
    , testPropertyDZ "powTo" prop_powTo prop_powTo
    , testPropertyDZ "sinTo" prop_sinTo prop_sinTo
    , testPropertyDZ "cosTo" prop_cosTo prop_cosTo
    , testPropertyDZ "tanTo" prop_tanTo prop_tanTo
    , testPropertyDZ "asinTo" prop_asinTo prop_asinTo
    , testPropertyDZ "acosTo" prop_acosTo prop_acosTo
    , testPropertyDZ "atanTo" prop_atanTo prop_atanTo
    , testPropertyDZ "sinhTo" prop_sinhTo prop_sinhTo
    , testPropertyDZ "coshTo" prop_coshTo prop_coshTo
    , testPropertyDZ "tanhTo" prop_tanhTo prop_tanhTo
    , testPropertyDZ "asinhTo" prop_asinhTo prop_asinhTo
    , testPropertyDZ "acoshTo" prop_acoshTo prop_acoshTo
    , testPropertyDZ "atanhTo" prop_atanhTo prop_atanhTo
    ]

{-
  getSumAbs:

    Double: [Failed]
*** Failed! Falsifiable (after 34 tests and 1 shrink):
fromList [48.00809219086489,-28.969358622663314,63.04503159639685,5.714612973491514,-5.312024428882535,60.56484370204321,20.339282130713237,-25.609691263383986,11.776052643355653,-72.22738635520177]
LHS
(used seed -5655945031033394469)

    Complex Double: [Failed]
*** Failed! Falsifiable (after 46 tests):
fromList [(-126.63684010249744) :+ (-25.732604105949363),48.30271669448071 :+ (-39.257907001851365),(-742.581416038894) :+ (-40.94772884905545),(-21.033310146590043) :+ 11.040319311574454,(-32.47272728995846) :+ (-75.38889146104881),34.61215019691849 :+ 63.02796474654257,25.495435934718248 :+ 14.86496960897311,285.1569736733574 :+ 20.592531068434027,(-75.03887610996705) :+ 62.568157533798285,13.719037128011735 :+ 49.186886152517744,(-23.33031638959207) :+ (-41.802953657765805),46.325514491999364 :+ (-40.3054052139413),(-97.8445174356864) :+ 62.27686352261467,(-210.08994152230045) :+ 23.4837419831381,(-47.12015960420664) :+ 9.502863461401958]
LHS

    Double: [Failed]
*** Failed! Falsifiable (after 61 tests):
VectorPair (fromList [24.074823693680273,166.2322847046224,-2.8103216076602617,6.544349669641704,157.92295537758915,-52.792112116283036,-56.334929255514766,12.079209274206134,-51.61912866052746,163.39401573569245,-49.34434955256682,-17.020978906318945,116.41247221011058]) (fromList [91.62365332716433,-46.23902424675873,151.71325376854423,34.663441501788725,-2.5921683575991397,-48.720845208466315,51.34009695413845,41.16706632522662,-52.1781413047511,18.0562293933656,-198.82058422020668,-4.213635327473868,40.59717128981193])
LHS
LHS
(used seed -6941761855370166739)

    Complex Double: [OK, passed 100 tests]

-}

prop_new_ t (Dim n) =
    (V.dim $ typed t $ V.create $ V.new_ n) === n

prop_new t (Dim n) e =
    (V.create $ V.new n e) === (typed t $ V.constant n e)

prop_newCopy t x =
    (V.create $ V.newCopy x) === x
  where
    _ = typed t x

prop_copyTo t (VectorPair x y) = runST $
    x `readOnlyVector` \mx ->
    y `mutatesToVector` x $ \my ->
        V.copyTo my mx
  where
    _ = typed t x

prop_swap t (VectorPair x y) = runST $
    x `mutatesToVector` y $ \mx ->
    y `mutatesToVector` x $ \my ->
        V.swap mx my
  where
    _ = typed t x

prop_read t (Index n i) =
    forAll (typed t `fmap` Test.vector n) $ \x -> runST $
        x `readOnlyVector` \mx -> do
            e <- V.read mx i
            return $ e === V.at x i

prop_write t (Index n i) e =
    forAll (Test.vector n) $ \x -> runST $
        x `mutatesToVector` (x `V.update` [(i,e)]) $ \mx -> do
            V.write mx i e
  where
    _ = e == t

prop_modify t (Index n i) (Blind f) =
    forAll (Test.vector n) $ \x -> runST $
        x `mutatesToVector`
            (typed t $ x `V.update` [(i, f $ V.at x i)]) $ \mx ->
                V.modify mx i f

prop_getElems t x =
    forAll arbitrary $ \(Blind f) -> runST $
        x `mutatesToVector` (V.map f x) $ \mx -> do
            es <- V.getElems mx
            V.mapTo mx f mx
            return $ es === V.elems (V.map f x)
  where
    _ = typed t x

prop_getElems' t x =
    forAll arbitrary $ \(Blind f) -> runST $
        x `mutatesToVector` (V.map f x) $ \mx -> do
            es <- V.getElems' mx
            V.mapTo mx f mx
            return $ es === V.elems x
  where
    _ = typed t x

prop_getAssocs t x =
    forAll arbitrary $ \(Blind f) -> runST $
        x `mutatesToVector` (V.map f x) $ \mx -> do
            ies <- V.getAssocs mx
            V.mapTo mx f mx
            return $ ies === V.assocs (V.map f x)
  where
    _ = typed t x

prop_getAssocs' t x =
    forAll arbitrary $ \(Blind f) -> runST $
        x `mutatesToVector` (V.map f x) $ \mx -> do
            ies <- V.getAssocs' mx
            V.mapTo mx f mx
            return $ ies === V.assocs x
  where
    _ = typed t x

prop_setElems t x =
    forAll (QC.vector n) $ \es -> runST $
        x `mutatesToVector` (V.fromList n es) $ \mx ->
            V.setElems mx es
  where
    n = V.dim x
    _ = typed t x

prop_setAssocs t (Assocs n ies) =
    forAll (Test.vector n) $ \x -> runST $
        x `mutatesToVector` (typed t $ V.update x ies) $ \mx ->
            V.setAssocs mx ies

prop_mapTo t (Blind f) = binaryProp t
    (\x -> V.map f x)
    (\dst mx -> V.mapTo dst f mx)

prop_zipWithTo t (Blind f) = ternaryProp t
    (\x y -> V.zipWith f x y)
    (\dst mx my -> V.zipWithTo dst f mx my)

prop_getSumAbs t x = runST $
    x `readOnlyVector` \mx -> do
        s <- V.getSumAbs mx
        return $ s === V.sumAbs x
  where
    _ = typed t x

prop_getNorm2 t x = runST $
    x `readOnlyVector` \mx -> do
        n <- V.getNorm2 mx
        return $ n === V.norm2 x
  where
    _ = typed t x

prop_getDot t (VectorPair x y) = runST $
    x `readOnlyVector` \mx ->
    y `readOnlyVector` \my -> do
        d <- V.getDot mx my
        return $ d === V.dot x y
  where
    _ = typed t x

prop_kroneckerTo t = sized $ \s -> resize (s `div` 2) $
    forAll arbitrary $ \x ->
    forAll arbitrary $ \y ->
    forAll (Test.vector $ V.dim x * V.dim y) $ \z -> runST $
        x `readOnlyVector` \mx ->
        y `readOnlyVector` \my ->
        z `mutatesToVector` (typed t $ V.kronecker x y) $ \mz ->
            V.kroneckerTo mz mx my

prop_addTo t = ternaryProp t V.add V.addTo
prop_subTo t = ternaryProp t V.sub V.subTo
prop_mulTo t = ternaryProp t V.mul V.mulTo
prop_negateTo t = binaryProp t V.negate V.negateTo
prop_conjugateTo t = binaryProp t V.conjugate V.conjugateTo
prop_absTo t = binaryProp t V.abs V.absTo
prop_signumTo t = binaryProp t V.signum V.signumTo

prop_divTo t = ternaryProp t V.div V.divTo
prop_recipTo t = binaryProp t V.recip V.recipTo

prop_sqrtTo t = binaryProp t V.sqrt V.sqrtTo
prop_expTo t = binaryProp t V.exp V.expTo
prop_logTo t = binaryProp t V.log V.logTo
prop_powTo t = ternaryProp t V.pow V.powTo
prop_sinTo t = binaryProp t V.sin V.sinTo
prop_cosTo t = binaryProp t V.cos V.cosTo
prop_tanTo t = binaryProp t V.tan V.tanTo
prop_asinTo t = binaryProp t V.asin V.asinTo
prop_acosTo t = binaryProp t V.acos V.acosTo
prop_atanTo t = binaryProp t V.atan V.atanTo
prop_sinhTo t = binaryProp t V.sinh V.sinhTo
prop_coshTo t = binaryProp t V.cosh V.coshTo
prop_tanhTo t = binaryProp t V.tanh V.tanhTo
prop_asinhTo t = binaryProp t V.asinh V.asinhTo
prop_acoshTo t = binaryProp t V.acosh V.acoshTo
prop_atanhTo t = binaryProp t V.atanh V.atanhTo

binaryProp :: (AEq e, Arbitrary e, Show e, Storable e, Testable a)
            => e
            -> (Vector e -> Vector e)
            -> (forall s . STVector s e -> STVector s e -> ST s a)
            -> Property
binaryProp t imm_f f = let
    prop2 t (VectorPair x y) = runST $ let _ = typed t x in
        x `readOnlyVector` \mx ->
        y `mutatesToVector` (imm_f x) $ \my -> f my mx
      where _ = typed t x
    in label "" prop2

ternaryProp :: (AEq e, Arbitrary e, Show e, Storable e, Testable a)
            => e
            -> (Vector e -> Vector e -> Vector e)
            -> (forall s . STVector s e -> STVector s e -> STVector s e -> ST s a)
            -> Property
ternaryProp t imm_f f = let
    prop3 t (VectorTriple x y z) = runST $
        x `readOnlyVector` \mx ->
        y `readOnlyVector` \my ->
        z `mutatesToVector` (imm_f x y) $ \mz -> f mz mx my
      where _ = typed t x
    in label "" prop3

readOnlyVector :: (Storable e, AEq e, Testable prop, Show e)
               => Vector e
               -> (STVector s e -> ST s prop)
               -> ST s Property
readOnlyVector x = mutatesToVector x x

mutatesToVector :: (Storable e, AEq e, Testable prop, Show e)
          => Vector e
          -> Vector e
          -> (STVector s e -> ST s prop)
          -> ST s Property
mutatesToVector x x_new f = do
    mx <- V.newCopy x
    prop <- f mx
    x' <- V.freeze mx
    return $ prop .&. (x' === x_new)

ext: hs

max_stars_repo_name: cartazio/hs-cblas

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
{-# LANGUAGE AllowAmbiguousTypes                      #-}
{-# LANGUAGE ConstraintKinds                          #-}
{-# LANGUAGE DataKinds                                #-}
{-# LANGUAGE DefaultSignatures                        #-}
{-# LANGUAGE FlexibleContexts                         #-}
{-# LANGUAGE FlexibleInstances                        #-}
{-# LANGUAGE GADTs                                    #-}
{-# LANGUAGE RankNTypes                               #-}
{-# LANGUAGE ScopedTypeVariables                      #-}
{-# LANGUAGE TypeApplications                         #-}
{-# LANGUAGE TypeFamilies                             #-}
{-# LANGUAGE TypeOperators                            #-}
{-# LANGUAGE UndecidableInstances                     #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}

module Backprop.Learn.Initialize (
    Initialize(..)
  , gInitialize
  , initializeNormal
  , initializeSingle
  ) where

import           Control.Monad.Primitive
import           Data.Complex
import           Data.Type.Length
import           Data.Type.NonEmpty
import           Data.Type.Tuple
import           GHC.TypeNats
import           Generics.OneLiner
import           Numeric.LinearAlgebra.Static.Vector
import           Statistics.Distribution
import           Statistics.Distribution.Normal
import           Type.Class.Known
import           Type.Family.List
import qualified Data.Vector.Generic                 as VG
import qualified Data.Vector.Generic.Sized           as SVG
import qualified Numeric.LinearAlgebra.Static        as H
import qualified System.Random.MWC                   as MWC

-- | Class for types that are basically a bunch of 'Double's, which can be
-- initialized with a given identical and independent distribution.
class Initialize p where
    initialize
        :: (ContGen d, PrimMonad m)
        => d
        -> MWC.Gen (PrimState m)
        -> m p

    default initialize
        :: (ADTRecord p, Constraints p Initialize, ContGen d, PrimMonad m)
        => d
        -> MWC.Gen (PrimState m)
        -> m p
    initialize = gInitialize

-- | 'initialize' for any instance of 'Generic'.
gInitialize
    :: (ADTRecord p, Constraints p Initialize, ContGen d, PrimMonad m)
    => d
    -> MWC.Gen (PrimState m)
    -> m p
gInitialize d g = createA' @Initialize (initialize d g)

initializeNormal
    :: (Initialize p, PrimMonad m)
    => Double                               -- ^ standard deviation
    -> MWC.Gen (PrimState m)
    -> m p
initializeNormal = initialize . normalDistr 0

-- | 'initialize' definition if @p@ is a single number.
initializeSingle
    :: (ContGen d, PrimMonad m, Fractional p)
    => d
    -> MWC.Gen (PrimState m)
    -> m p
initializeSingle d = fmap realToFrac . genContVar d

instance Initialize Double where
    initialize = initializeSingle
instance Initialize Float where
    initialize = initializeSingle

-- | Initializes real and imaginary components identically
instance Initialize a => Initialize (Complex a) where

instance Initialize T0
instance (Initialize a, Initialize b) => Initialize (T2 a b)
instance (Initialize a, Initialize b, Initialize c) => Initialize (T3 a b c)

instance (ListC (Initialize <$> as), Known Length as) => Initialize (T as) where
    initialize d g = constTA @Initialize (initialize d g) known

instance (Initialize a, ListC (Initialize <$> as), Known Length as) => Initialize (NETup (a ':| as)) where
    initialize d g = NET <$> initialize d g
                         <*> initialize d g

instance Initialize ()
instance (Initialize a, Initialize b) => Initialize (a, b)
instance (Initialize a, Initialize b, Initialize c) => Initialize (a, b, c)

instance (VG.Vector v a, KnownNat n, Initialize a) => Initialize (SVG.Vector v n a) where
    initialize d = SVG.replicateM . initialize d

instance KnownNat n => Initialize (H.R n) where
    initialize d = fmap vecR . initialize d
instance KnownNat n => Initialize (H.C n) where
    initialize d = fmap vecC . initialize d

instance (KnownNat n, KnownNat m) => Initialize (H.L n m) where
    initialize d = fmap vecL . initialize d
instance (KnownNat n, KnownNat m) => Initialize (H.M n m) where
    initialize d = fmap vecM . initialize d

constTA
    :: forall c as f. (ListC (c <$> as), Applicative f)
    => (forall a. c a => f a)
    -> Length as
    -> f (T as)
constTA x = go
  where
    go :: forall bs. ListC (c <$> bs) => Length bs -> f (T bs)
    go LZ     = pure TNil
    go (LS l) = (:&) <$> x <*> go l

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
{-# LANGUAGE TypeApplications #-}
module ArrayFire.StatisticsSpec where

import ArrayFire    hiding (not)

import Data.Complex
import Test.Hspec

spec :: Spec
spec =
  describe "Statistics spec" $ do
    it "Should find the mean" $ do
      mean (vector @Double 10 [1..]) 0
        `shouldBe`
           5.5
    it "Should find the weighted-mean" $ do
      meanWeighted (vector @Double 10 [1..]) (vector @Double 10 [1..]) 0
        `shouldBe`
           7.0
    it "Should find the variance" $ do
      var (vector @Double 8 [1..8]) False 0
        `shouldBe`
           6.0
    it "Should find the weighted variance" $ do
      varWeighted (vector @Double 8 [1..]) (vector @Double 8 (repeat 1)) 0
        `shouldBe`
           5.25
    it "Should find the standard deviation" $ do
      stdev (vector @Double 10 (cycle [1,-1])) 0
        `shouldBe`
           1.0
    it "Should find the covariance" $ do
      cov (vector @Double 10 (repeat 1)) (vector @Double 10 (repeat 1)) False
        `shouldBe`
           0.0
    it "Should find the median" $ do
      median (vector @Double 10 [1..]) 0
        `shouldBe`
           5.5
    it "Should find the mean of all elements across all dimensions" $ do
      fst (meanAll (matrix @Double (2,2) [[10,10],[10,10]]))
        `shouldBe`
           10
    it "Should find the weighted mean of all elements across all dimensions" $ do
      fst (meanAllWeighted (matrix @Double (2,2) [[10,10],[10,10]]) (matrix @Double (2,2) [[10,10],[10,10]]))
        `shouldBe`
           10
    it "Should find the variance of all elements across all dimensions" $ do
      fst (varAll (vector @Double 10 (repeat 10)) False)
        `shouldBe`
           0
    it "Should find the weighted variance of all elements across all dimensions" $ do
      fst (varAllWeighted (vector @Double 10 (repeat 10)) (vector @Double 10 (repeat 10)))
        `shouldBe`
           0
    it "Should find the stdev of all elements across all dimensions" $ do
      fst (stdevAll (vector @Double 10 (repeat 10)))
        `shouldBe`
           0
    it "Should find the median of all elements across all dimensions" $ do
      fst (medianAll (vector @Double 10 [1..]))
        `shouldBe`
           5.5
    it "Should find the correlation coefficient" $ do
      fst (corrCoef (vector @Int 10 [1..] ) ( vector @Int 10 [10,9..] ))
        `shouldBe`
           (-1.0)
    it "Should find the top k elements" $ do
      let (vals,indexes) = topk ( vector @Double 10 [1..] ) 3 TopKDefault
      vals `shouldBe` vector @Double 3 [10,9,8]
      indexes `shouldBe` vector @Double 3 [9,8,7]

max_stars_repo_name: danielkroeni/arrayfire-haskell

max_issues_repo_name: danielkroeni/arrayfire-haskell

max_forks_repo_name: danielkroeni/arrayfire-haskell


text:
{-# LANGUAGE GADTs #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}

module Control.Quantum
    ( Amplitude
    , WaveFunction
    , Quantum
    , quantize
    , MonadQuantum
    , runQuantum
    , runQuantum'
    , measure
    , measurements
    )
    where

import qualified Data.Map  as M
import           Data.Map (Map)
import           Data.Monoid
import           Control.Monad
import           Control.Applicative
import           Control.Monad.Random as R
import           Data.Complex

nearZero :: Double -> Bool
nearZero x = abs x <= 1e-12

type Amplitude = Complex Double

-- | The Wave function is the distribution of amplitudes over the
-- possible states of the underlying data type
type WaveFunction a = [(a, Amplitude)]

class (Monad m) => MonadQuantum m where
    quantize  :: (Ord a) => WaveFunction a -> m a
    quantize' ::            WaveFunction a -> m a
    condition :: Bool -> m ()

liftQ  :: (Ord b, MonadQuantum m) => (a -> b) -> m a -> m b
liftQ f q1 = q1 >>= always . f

liftQ2 :: (Ord c, MonadQuantum m) => (a -> b -> c) -> m a -> m b -> m c
liftQ2 f p1 p2 = do x1 <- p1
                    x2 <- p2
                    always (f x1 x2)

liftQ3 :: (Ord d, MonadQuantum m) => (a -> b -> c -> d) -> m a -> m b -> m c -> m d
liftQ3 f p1 p2 p3 = do x1 <- p1
                       x2 <- p2
                       x3 <- p3
                       always (f x1 x2 x3)

always a = quantize [(a, 1.0:+0.0)]
always' a = quantize' [(a, 1.0:+0.0)]

-- | Describes the quantized version of an arbitrary data type.
-- The underlying data type defines the basis that the state is measured in.
data Quantum a where
    Quantum :: Ord a => Map (Maybe a) Amplitude -> Quantum a
    QuantumAny :: [(Maybe a, Amplitude)]  -> Quantum a

noState  :: (Ord a) => Quantum a
noState  = Quantum (M.singleton Nothing (1.0:+0.0))

noState' :: Quantum a
noState' = QuantumAny [(Nothing, 1.0 :+ 0.0)]

deriving instance (Show a) => Show (Quantum a)

instance Functor Quantum where
    fmap = liftM

instance Applicative Quantum where
    pure  = return
    (<*>) = ap

instance Monad Quantum where
    return  = always'

    m >>= f = if unitary then next
                         else error "Non-unitary transformation applied to quantum state!"
        where
            unitary = nearZero $ l2norm (map snd $ toList' next) - 1.0
            next = collect [multAmpl q (go a) | (a, q) <- toList' m]
            go = maybe noState' f

multAmpl :: Amplitude -> Quantum a -> Quantum a
multAmpl q (Quantum x) = Quantum $ M.map (* conjugate q) x
multAmpl q (QuantumAny x) = QuantumAny [ (a, q * r) | (a, r) <- x ]

toList' :: Quantum a -> [(Maybe a, Amplitude)]
toList' (Quantum x) = M.toList x
toList' (QuantumAny x) = x

toList :: (Ord a) => Quantum a -> [(Maybe a, Amplitude)]
toList (Quantum x) = M.toList x
toList (QuantumAny x) = merge x

collect :: [Quantum a] -> Quantum a
collect [ ]        = QuantumAny []
collect [x]        = x
collect (Quantum x:t) = case collect t of
                        Quantum y -> Quantum (M.unionWith (+) x y)
                        QuantumAny y -> Quantum (M.unionWith (+) x (M.fromList y))
collect (QuantumAny x:t) = case collect t of
                        Quantum y -> Quantum (M.unionWith (+) (M.fromList x) y)
                        QuantumAny y -> QuantumAny (x ++ y)

merge :: (Ord a) => WaveFunction a -> WaveFunction a
merge = M.toList . M.fromListWith (+)

instance MonadQuantum Quantum where
    quantize    = Quantum . M.fromListWith (+) . normalize
    quantize'   = QuantumAny . normalize
    condition test = if test then always () else noState

instance (Ord a, Monoid a) => Monoid (Quantum a) where
    mempty  = always mempty
    mappend = liftQ2 mappend

normalize :: WaveFunction a -> [(Maybe a, Amplitude)]
normalize xs = map (\(a, q) -> (Just a, q / total)) xs
    where
        total = sum $ zipWith (*) (map conjugate ampl) ampl
        ampl  = map snd xs

-- | Remove all impossible states from the quantum state and renormalize
collapse :: [(Maybe a, Amplitude)] -> WaveFunction a
collapse xs = [ (x, q / (norm:+0)) | (Just x, q) <- xs ]
    where
        norm = l2norm [ q | (Just x, q) <- xs ]

l2norm :: [Amplitude] -> Double
l2norm xs = sqrt $ sum $ map (\x -> (magnitude x)**2) xs

-- | Extracts the wave function from a quantum state.
runQuantum :: (Ord a) => Quantum a -> WaveFunction a
runQuantum = collapse . toList

runQuantum' :: Quantum a -> WaveFunction a
runQuantum' = collapse . toList'

-- | Measures a quantum state. Returns a single realization of the underlying data type
-- with probability equal to the squared magnitude of the amplitude of that realization.
-- Note that if this Module would be backed by a real quantum processor, this would be
-- the only valid way to extract information from the Monad.
measure :: (Ord a) => Quantum a -> IO a
measure q = do
  x <- evalRandIO $ R.fromList $ map (\(a, b) -> (a, toRational b)) $ measurements q
  return x

-- | Converts a quantum state into a probability state through measurement.
-- Equivalent to performing repeated measurements on equally prepared quantum states
-- and noting the frequency of each possible realization.
measurements :: (Ord a) => Quantum a -> [(a, Double)]
measurements = f . runQuantum
    where f xs = [(x, (magnitude q)**2) | (x, q) <- xs]

max_stars_repo_name: ibab/haskell-quantum

max_issues_repo_name: ibab/haskell-quantum

max_forks_repo_name: ibab/haskell-quantum


text:
import Statistics.Distribution -- (mean, stdDev)
import Statistics.Distribution.Gamma

gd :: GammaDistribution
gd = gammaDistr a b
       where
         a = 3.1 -- shape parameter k
         b = 2.0 -- scale parameter theta

m = mean gd
s = stdDev gd

main = do
  putStrLn $ "The mean is: " ++ show m ++ " and the standard deviation: " ++ show s
  let k = (22.5 - m)/s
  putStrLn $ "k-value is: " ++ show k
  putStrLn $ "P( Y>22.5)= " ++ show (complCumulative gd 22.5)
  let tchebycheff = 1/k^2
  putStrLn $ "Using Tchebysheff's theorem, P( |Y-m|> " ++ show k ++ "*s ) <= 1/k^2 = " ++ show tchebycheff

max_stars_repo_name: raydsameshima/Stat

max_issues_repo_name: raydsameshima/Stat

max_forks_repo_name: raydsameshima/Stat


text:
{-# language TypeFamilies #-}
{-# language FlexibleInstances #-}
{-# language FlexibleContexts #-}

module Geometry.Circle
  ( Circle, Circular (..)
  , trivialCircle, mkCircle
  ) where

import Data.Complex
import Data.Semigroup
import Geometry.Base

-- | Class for circle and decorated circle
class (Trans c, Manifold c, Curve c, ClosedCurve c, Figure c) =>
  Circular c where
  {-# MINIMAL toCircle, asCircle #-}
  toCircle :: c -> Circle
  asCircle :: Circle -> c

  -- | Center of the circle.
  center :: c -> Cmp
  center = center . toCircle

  -- | Radius of the circle
  radius :: c -> Double
  radius = radius . toCircle

  -- | The angle of the starting point.
  phaseShift :: c -> Direction
  phaseShift = phaseShift . toCircle

  -- | The orientation of the circle (positive -- CCW).
  orientation :: c -> Direction
  orientation = orientation . toCircle

-- | Represents a circle with given center, radius vector and tangent direction.
data Circle = Circle Cmp Cmp Double
  deriving Show

instance Circular Circle where
  toCircle = id
  asCircle = id
  center (Circle c _ _) = c
  radius (Circle _ r _) = norm r
  orientation (Circle _ _ o) = asDeg (signum o)
  phaseShift (Circle _ r _) = angle r

instance Circular c => Circular (Maybe c) where
  toCircle = maybe trivialCircle toCircle
  asCircle = Just . asCircle

-- | The trivial circle with zero radius.
trivialCircle :: Circle
trivialCircle = mkCircle 0 0

-- | The constructor for a circle with given center and radius.
-- The radius-vector has zero angle, and circle is CCW-oriented.
mkCircle :: Double -> Cmp -> Circle
mkCircle r c = Circle c (r:+0) 1

instance Eq Circle where
  c1 == c2 = radius c1 ~= radius c2 &&
             center c1 ~= center c2 &&
             phaseShift c1 ~= phaseShift c2

instance Trans Circle where
  transform t (Circle c r o) = Circle c' r' o'
    where c' = transformCmp t c
          r' = transformCmp t (c + r) - c'
          o' = transformOrientation t * o

instance Manifold Circle where
  type Domain Circle = Cmp
  param c t = center c + mkPolar (radius c) (rad x)
    where
      ph = turns $ phaseShift c
      x = asTurns $ ph + t * deg (orientation c)

  project c p = turns (x - ph)
    where
      ph = phaseShift c
      x = orientation c * azimuth (center c) p

  isContaining c p = distance p (center c) ~= radius c
  unit _ = 2 * pi

  distanceTo p c = abs (distance p (center c) - radius c)

instance Curve Circle where
  normal c t = azimuth (center c) (c @-> t)
  tangent c t = normal c t + 90 * orientation c

instance ClosedCurve Circle where
  location c p = res
    where res | r' ~= radius c = OnCurve
              | r' < radius c  = Inside
              | r' > radius c  = Outside
          r' = distance (cmp p) (center c)

instance Figure Circle where
  isTrivial c = radius c <= 0
  refPoint = center
  box cir = ((Min x1, Min y1), (Max x2, Max y2))
    where c = center cir
          r = radius cir
          x1:+y1 = c-(r:+r)
          x2:+y2 = c+(r:+r)

max_stars_repo_name: MethaHardworker/geometry

max_issues_repo_name: MethaHardworker/geometry

max_forks_repo_name: MethaHardworker/geometry


text:
{-# LANGUAGE TypeApplications #-}

module MinkowskiScalarWave where

import qualified Data.Vector.Unboxed as U
import Data.VectorSpace
import qualified Numeric.LinearAlgebra as L
import SSSFN

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

-- | Scalar wave in double-null coordinates:
-- > (- dt^2 + dx^2) u = 0
-- > t = v + u   dt = dv + du
-- > x = v - u   dx = dv - du
-- > (- (dv + du)^2 + (dv - du)^2) u = 0
-- > - 4 dv du u = 0
dalembert :: (Eq a, L.Field a, U.Unbox a) => Op a
dalembert = -4 *^ (derivV * derivU)

incomingU :: forall a. (L.Container L.Vector a, Fractional a, Ord a) => Op a
incomingU = omap isneg (bndU @a)
  where
    isneg a = fromIntegral (fromEnum (a < 0))

incomingV :: forall a. (L.Container L.Vector a, Fractional a, Ord a) => Op a
incomingV = omap isneg (bndV @a)
  where
    isneg a = fromIntegral (fromEnum (a < 0))

originUV :: forall a. (L.Container L.Vector a, Fractional a, Ord a) => Op a
originUV = ozipWith (*) incomingU incomingV

potential :: forall a. (L.Container L.Vector a, Floating a) => Grid a
potential = gzipWith (\u v -> 0) (coordU @a) (coordV @a)

incoming :: forall a. (L.Container L.Vector a, Floating a) => Grid a
incoming = gsample bcond

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

op :: (L.Field a, Ord a, U.Unbox a) => Op a
op = intr * dalembert + bndu * incomingU + bndv * incomingV + orig * originUV
  where
    bndu = incomingU - originUV
    bndv = incomingV - originUV
    orig = originUV
    intr = 1 - (bndu + bndv + orig)

rhs :: (Floating a, L.Numeric a, Ord a) => Grid a
rhs =
  intr #> potential
    ^+^ bndu #> incoming
    ^+^ bndv #> incoming
    ^+^ orig #> incoming
  where
    bndu = incomingU - originUV
    bndv = incomingV - originUV
    orig = originUV
    intr = 1 - (bndu + bndv + orig)

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

bcond :: Floating a => (a, a) -> a
bcond (u, v) = cos (pi / 4 * t) * cos (pi / 4 * x)
  where
    t = v + u
    x = v - u

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

main :: IO ()
main =
  do
    putStrLn "1+1D Scalar Field in the Minkowski spacetime in double-Null coordinates"
    -- putStrLn $ "intr " ++ show (intr :: Op Double)
    -- putStrLn $ "bndu " ++ show (bndu :: Op Double)
    -- putStrLn $ "bndv " ++ show (bndv :: Op Double)
    -- putStrLn $ "orig " ++ show (orig :: Op Double)
    -- putStrLn $ "dal " ++ show (omap approx (dalembert @Double))
    -- putStrLn $ "op " ++ show (omap approx (op @Double))
    -- let inc = incoming @Double
    -- putStrLn $ "inc " ++ show (gmap approx inc)
    -- putStrLn $ "rhs " ++ show (gmap approx (rhs @Double))
    -- let ires = dalembert #> inc
    -- putStrLn $ "ires " ++ show (gmap approx ires)
    -- putStrLn $ "|ires| " ++ show (gmaxabs ires)
    -- let jres = op #> inc ^-^ rhs
    -- putStrLn $ "jres " ++ show (gmap approx jres)
    -- putStrLn $ "|jres| " ++ show (gmaxabs jres)

    let sol = solve op (rhs @Double)
    putStrLn $ "sol " ++ show (gmap approx sol)
    let res = op #> sol ^-^ rhs
    putStrLn $ "res " ++ show (gmap approx res)
    putStrLn $ "|res| " ++ show (gmaxabs res)
    let err = sol ^-^ incoming
    putStrLn $ "err " ++ show (gmap approx err)
    putStrLn $ "|err| " ++ show (gmaxabs err)
    putStrLn "Done."
  where
    approx :: RealFrac a => a -> a
    approx x = fromInteger (round (1.0e+4 * x)) / 1.0e+4
-- intr = 1 - (incomingU + incomingV - originUV)
-- bndu = incomingU - originUV
-- bndv = incomingV - originUV
-- orig = originUV

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

-- quantum measurements: the quantum system measures the state of the
-- classical system

-- Beatrice Bonga: conserved current from varying boundary term from
-- lagrangian <https://arxiv.org/abs/1911.04514>

max_stars_repo_name: eschnett/sssfn

max_issues_repo_name: eschnett/sssfn

max_forks_repo_name: eschnett/sssfn


text:
module Core(
    RenderSettings(..),
    module Scene.Base,
    module Scene.Object.Plane,
    module Scene.Material.Diffuse,
    module Scene.Object.SimpleSphere,
    module Math.Color,
    module Math.Ray,
    module Scene.Material.Skylike
    ) where

import Scene.Base(Trace, traceRay, Scene, Material)
--import Scene.Object.SimpleSphere(sphere)
import Scene.Object.Plane
import Scene.Object.SimpleSphere
import Scene.Material.Skylike
import Scene.Material.Diffuse
import Math.Color
import Math.Ray
import Numeric.LinearAlgebra
import Scene.Object.Plane
import GHC.Generics (Generic)
import Control.DeepSeq

data RenderSettings  = RenderSettings {
    background :: Color,
    width :: Int,
    height :: Int,
    topLeft :: Vec3,
    topRight :: Vec3,
    bottomRight :: Vec3,
    origin :: Vec3,
    path :: FilePath,
    antialiasing :: Int}

max_stars_repo_name: danrocag/raytracer

max_issues_repo_name: danrocag/raytracer

max_forks_repo_name: danrocag/raytracer


text:
module Nnaskell ( NN(..)
                , randomNN
                , activateNN
                , cost
                , optimizeCost
                , loadNN
                , saveNN
                ) where

import Data.Function
import Data.List
import Control.Monad
import System.Random
import Numeric.LinearAlgebra.Data
import Numeric.LinearAlgebra.HMatrix

data NN = NN { nnWs :: [Matrix R]
             , nnBs :: [Vector R]
             } deriving (Show, Read)

randomVec :: Int -> IO (Vector R)
randomVec n = vector <$> (replicateM n $ randomRIO (-1.0, 1.0))

randomMatrix :: Int -> Int -> IO (Matrix R)
randomMatrix n m = matrix m <$> (replicateM ((*) n m) $ randomRIO (-1.0, 1.0))

randomLayer :: Int -> Int -> IO (Matrix R, Vector R)
randomLayer n m = liftM2 (\ws bs -> (ws, bs)) (randomMatrix m n) (randomVec m)

randomNN :: [Int] -> IO NN
randomNN ns = do (ws, bs) <- unzip <$> (sequence $ map (uncurry randomLayer) $ pairScan ns)
                 return $ NN { nnWs = ws
                             , nnBs = bs
                             }

sigmoid :: R -> R
sigmoid t = 1 / (1 + exp (-1 * t))

sigmoidVec :: Vector R -> Vector R
sigmoidVec = fromList . map sigmoid . toList

activateLayer :: Vector R -> (Matrix R, Vector R) -> Vector R
activateLayer as (ws, bs) = sigmoidVec (ws #> as + bs)

activateNN :: NN -> Vector R -> Vector R
activateNN nn as = foldl activateLayer as $ zip ws bs
    where ws = nnWs nn
          bs = nnBs nn

pairScan :: [a] -> [(a, a)]
pairScan (x:y:rest) = (x, y) : pairScan (y:rest)
pairScan _ = []

cost :: [(Vector R, Vector R)] -> NN -> R
cost td nn = (sum $ map inputCost td) / n
    where inputCost (input, output) = sumElements $ cmap (** 2) (activateNN nn input - output)
          n = fromIntegral $ length td

stepBias :: NN -> Int -> R -> NN
stepBias nn updateIdx s = nn { nnBs = map (\(v, idx) -> if idx <= updateIdx && updateIdx < idx + size v
                                                        then accum v (+) [(updateIdx - idx, s)]
                                                        else v)
                                      $ zip bs
                                      $ scanl (\idx v -> idx + size v) 0 bs
                             }
    where bs = nnBs nn

stepWeight :: NN -> Int -> R -> NN
stepWeight nn updateIdx s = nn { nnWs = map (\(mx, idx) -> let (n, m) = size mx
                                                               effIdx = updateIdx - idx
                                                           in if idx <= updateIdx && updateIdx < idx + (n * m)
                                                              then accum mx (+) [((effIdx `div` m, effIdx `mod` m), s)]
                                                              else mx)
                                        $ zip ws
                                        $ scanl (\idx m -> idx + uncurry (*) (size m)) 0 ws
                               }
    where ws = nnWs nn

class Domain d where
    countArgs :: d -> Int
    stepArg :: d -> Int -> R -> d

instance Domain NN where
    countArgs nn = bsCount + wsCount
        where bsCount = sum $ map (length . toList) $ nnBs nn
              wsCount = sum $ map length $ concatMap toLists $ nnWs nn
    stepArg nn idx s = if idx < bsCount
                       then stepBias nn idx s
                       else stepWeight nn (idx - bsCount) s
        where bsCount = sum $ map (length . toList) $ nnBs nn

optimizeCost :: Domain d => (d -> R) -> d -> [d]
optimizeCost cost d = scanl optimizeStep d $ cycle [0 .. n - 1]
    where n = countArgs d
          optimizeStep d idx = minimumBy (compare `on` cost) $ map (stepArg d idx) [-step, 0, step]
          step = 0.1

loadNN :: FilePath -> IO NN
loadNN filePath = read <$> readFile filePath

saveNN :: FilePath -> NN -> IO ()
saveNN filePath nn = writeFile filePath $ show nn

max_stars_repo_name: tsoding/nnaskell

max_issues_repo_name: tsoding/nnaskell

max_forks_repo_name: tsoding/nnaskell


text:
{-# LANGUAGE ScopedTypeVariables, RankNTypes #-}

module Dna.Graph (
readConnection,readAllConnection,
allPossibleNode,graphFromConnection,
Graph,euclidianCycle,
unbalanceNodes,
canBeBalance, balance, initOverlaGraph
)
where

import Data.List
import Data.Tuple
import Data.Functor
import Data.STRef
import Control.Monad
import Data.Maybe
import Control.Applicative
import Numeric.LinearAlgebra (rows, (!), cols, toLists, (?), (¿), toList)
import Numeric.LinearAlgebra.Data (Matrix, I)
import Numeric.LinearAlgebra.Devel
import Control.Monad.ST

import qualified Data.ByteString.Lazy.Char8 as C
import qualified Data.Set as S

data Graph a = GraphInternal {
    node :: S.Set a,
    connection :: Matrix I,
    addCon :: Maybe (a,a)
}

instance (Ord a, Show a) => Show (Graph a) where
    show graph = fromMaybe matrixShow (do
                    (from,to) <- addCon graph
                    return $ matrixShow ++ "\n With Link: " ++ show from  ++ "->" ++ show to)
        where
            linkCount row col = connection graph ! row ! col
            isLink row col =   linkCount row col >= 1
            matrixShow = intercalate "\n"
                [link | row <- [0 .. rows (connection graph) - 1],
                        any (isLink row) [0 .. cols (connection graph) - 1],
                        link <- [show (S.elemAt row (node graph)) ++ " -> " ++
                                 intercalate "," (concatMap (\col -> replicate (fromIntegral $ linkCount row col) $ show  $ S.elemAt col (node graph))
                                                    (filter (isLink row) [0 .. cols (connection graph) - 1]))]]

initOverlaGraph :: Ord a =>  [a] -> (forall s. STMatrix s I -> (a -> Int) -> ST s ()) -> Graph a
initOverlaGraph nodes f = GraphInternal setNode (runSTMatrix $ do
        m <- newMatrix 0 lenNodes lenNodes
        f m (`S.findIndex` setNode)
        return m
    ) Nothing
    where
        setNode = S.fromList nodes
        lenNodes  = length setNode

readConnection :: (Read a) => String -> (a, [a])
readConnection info = case lex info of
    [(node, left)] -> (read node, readListOfNode (snd $ head $ lex left))

readListOfNode :: (Read a) => String -> [a]
readListOfNode nodes = case lex nodes of
    [("", "")] -> []
    [(node, left)] -> read node : readListOfNode (snd $ head $ lex left)

readAllConnection :: (Read a) => String -> [(a, [a])]
readAllConnection = map readConnection . lines

allPossibleNode :: [(a,[a])] -> [a]
allPossibleNode = concatMap $ uncurry (:)

graphFromConnection ::Ord a =>  [(a,[a])] -> Graph a
graphFromConnection connections = GraphInternal
        nodes
        (runSTMatrix $ do
            m <- newMatrix 0 nodesLen nodesLen
            forM_ connections $ \ (from,connectee) ->
                forM_ connectee $ \ to ->
                    modifyMatrix m (findIndex from) (findIndex to) (1+)
            return m)
        Nothing
    where
        nodes = S.fromList $ allPossibleNode connections
        findIndex = flip S.findIndex nodes
        nodesLen = S.size nodes

euclidianCycle :: Ord a => Graph a -> Maybe [(a,a)]
euclidianCycle g = do
    (GraphInternal nodes m newCon) <- balance g
    conn <- runST $ do
        traverseM <- thawMatrix m
        path <- findEuclidianCycleAt nodes traverseM (Just 0)
        findNext nodes path traverseM
    return $ maybe conn (\(from,to) -> init $ uncurry (++) (swap (break ( (to == ).fst ) conn))) newCon
  where findNext nodes path m = do
            newPath <- findNextEuclidianPath nodes m path
            if isJust newPath
                then findNext nodes newPath m
                else return path

findNextEuclidianPath :: Ord a => S.Set a -> STMatrix s I -> Maybe [(a,a)] -> ST s (Maybe [(a,a)])
findNextEuclidianPath nodes mat path = do
    inMat <- freezeMatrix mat
    let newPath = do
            foundPath <- path
            let (left,right) = break (\(from,to) -> sum (map ti $ toList (inMat ! S.findIndex from nodes)) > 0) foundPath
            if null right
                then Nothing
                else Just (right ++ left)
    let maybeNext = fmap ((`S.findIndex` nodes) . fst . head) newPath
    maybeNextPath <- findEuclidianCycleAt nodes mat maybeNext
    return $ liftM2 (++) maybeNextPath newPath

findEuclidianCycleAt :: Ord a => S.Set a -> STMatrix s I -> Maybe Int -> ST s (Maybe [(a,a)])
findEuclidianCycleAt nodes mat startNode = if isJust startNode
                                                then findNextNod startNode
                                                else return Nothing
    where
        findNextNod n =  do
            next <- findNextEuclidianAt nodes mat (fromJust n)
            let cycleEnd = liftM2 (==) next startNode
            case (cycleEnd, next) of
                (Just True, Just to) -> return (Just [(S.elemAt (fromJust n) nodes,S.elemAt to nodes)])
                (Just False, to) -> do
                    list <- findNextNod to
                    return $ liftM2 (:) (Just (S.elemAt (fromJust n) nodes,S.elemAt (fromJust to) nodes)) list
                _ -> return Nothing

findNextEuclidianAt :: Ord a => S.Set a -> STMatrix s I -> Int -> ST s (Maybe Int)
findNextEuclidianAt nodes mat n = do
    con <- mapM (readMatrix mat n) [0..(S.size nodes - 1)]
    let next = find ((>0).snd) (zip [0..] con)
    when (isJust next) $ modifyMatrix mat n (fst . fromJust $ next) (subtract 1)
    return $ fst <$> next

mapFst :: (a -> b) -> [(a,c)] -> [(b,c)]
mapFst f = map (uncurry $ (,) . f)

canBeBalance :: Ord a => Graph a -> Bool
canBeBalance g = case unbal of
    [(_,-1), (_,1)] -> True
    [(_,1), (_,-1)] -> True
    [] -> True
    _ -> False
  where unbal = unbalanceNodes g

balance :: Ord a => Graph a -> Maybe (Graph a)
balance g@(GraphInternal node connection addCon) = case unbal of
    [(to,-1), (from,1)] -> return (GraphInternal node (fst (mutable (correctingMatrix (S.findIndex from node) (S.findIndex to node)) connection)) (Just (from,to)))
    [(from,1), (to,-1)] -> return (GraphInternal node (fst (mutable (correctingMatrix (S.findIndex from node) (S.findIndex to node)) connection)) (Just (from,to)))
    [] -> return g
    _ -> Nothing
  where
    unbal = unbalanceNodes g
    correctingMatrix f t _ m = modifyMatrix m f t (+1)

unbalanceNodes :: Ord a => Graph a -> [(a,Int)]
unbalanceNodes (GraphInternal nodes m _) =filter ((0 /= ) . snd) (mapFst (`S.elemAt` nodes) (map (\idx -> (idx,eulerianNode idx m)) [0..(S.size nodes - 1)]))

eulerianNode :: Int -> Matrix I -> Int
eulerianNode n m = ti (sum (concat colConnection)) - ti (sum (concat rowConnection))
    where
        colConnection = toLists $ m ¿ [n]
        rowConnection = toLists $ m ? [n]

eulerianNodeST :: Int -> STMatrix s I -> ST s Bool
eulerianNodeST n m = do
    colm <- extractMatrix m AllRows (Col n)
    rowm <- extractMatrix m (Row n) AllCols
    return $ (sum . concat . toLists) colm == (sum . concat . toLists) rowm

max_stars_repo_name: mathk/NautilusDna

max_issues_repo_name: mathk/NautilusDna

max_forks_repo_name: mathk/NautilusDna


text:
module Scheme.Parser
  ( readExpr
  ) where

import Control.Monad (liftM)
import Control.Monad.Error
import Data.Array
import Data.Char (toLower)
import Data.Complex (Complex(..))
import Data.Ratio (Rational, (%))
import Numeric (readHex, readOct)
import Scheme.Data
import System.Environment
import Text.ParserCombinators.Parsec hiding (spaces)

readExpr :: String -> ThrowsError LispVal
readExpr input =
  case parse parseExpr "lisp" input of
    Left err -> throwError $ Parser err
    Right val -> return val

--
-- LispVal Parsers
--
parseExpr :: Parser LispVal
parseExpr =
  parseAtom <|> parseString <|> parseNumber <|> parseQuoted <|> parseQuasiQuoted <|> do
    char '('
    x <- try parseList <|> parseDottedList
    char ')'
    return x

parseQuasiQuoted :: Parser LispVal
parseQuasiQuoted = do
  char '`'
  x <- parseExpr
  return $ List [Atom "quasiquote", x]

parseUnQuote :: Parser LispVal
parseUnQuote = do
  char ','
  x <- parseExpr
  return $ List [Atom "unquote", x]

parseAtom :: Parser LispVal
parseAtom = do
  first <- letter <|> symbol
  rest <- many (letter <|> digit <|> symbol)
  (return . Atom) (first : rest)

parseVector :: Parser LispVal
parseVector = do
  arrayValues <- sepBy parseExpr spaces
  return $ Vector (listArray (0, (length arrayValues - 1)) arrayValues)

parseNumber :: Parser LispVal
parseNumber = parsePlainNumber <|> parseRadixNumber

parsePlainNumber :: Parser LispVal
parsePlainNumber = many1 digit >>= return . Number . read

parseRadixNumber :: Parser LispVal
parseRadixNumber =
  char '#' >> (parseDecimal <|> parseBinary <|> parseOctal <|> parseHex)

parseDecimal :: Parser LispVal
parseDecimal = do
  char 'd'
  n <- many1 digit
  (return . Number . read) n

parseBinary :: Parser LispVal
parseBinary = do
  char 'b'
  n <- many $ oneOf "01"
  (return . Number . bin2int) n

parseOctal :: Parser LispVal
parseOctal = do
  char 'o'
  n <- many $ oneOf "01234567"
  (return . Number . (readWith readOct)) n

parseHex :: Parser LispVal
parseHex = do
  char 'x'
  n <- many $ oneOf "0123456789abcdefABCDEF"
  (return . Number . (readWith readHex)) n

parseRatio :: Parser LispVal
parseRatio = do
  num <- fmap read $ many1 digit
  char '/'
  denom <- fmap read $ many1 digit
  (return . Ratio) (num % denom)

parseFloat :: Parser LispVal
parseFloat = do
  whole <- many1 digit
  char '.'
  decimal <- many1 digit
  return $ Float (read (whole ++ "." ++ decimal))

parseComplex :: Parser LispVal
parseComplex = do
  r <- fmap toDouble (try parseFloat <|> parsePlainNumber)
  char '+'
  i <- fmap toDouble (try parseFloat <|> parsePlainNumber)
  char 'i'
  (return . Complex) (r :+ i)
  where
    toDouble (Float x) = x
    toDouble (Number x) = fromIntegral x

parseString :: Parser LispVal
parseString = do
  char '"'
  s <- many (escapedChars <|> (noneOf ['\\', '"']))
  char '"'
  (return . String) s

parseChar :: Parser LispVal
parseChar = do
  string "#\\"
  s <- many1 letter
  return $
    case (map toLower s) of
      "space" -> Char ' '
      "newline" -> Char '\n'
      [x] -> Char x

parseBool :: Parser LispVal
parseBool = do
  char '#'
  c <- oneOf "tf"
  return $
    case c of
      't' -> Bool True
      'f' -> Bool False

parseList :: Parser LispVal
parseList = liftM List $ sepBy parseExpr spaces

parseDottedList :: Parser LispVal
parseDottedList = do
  head <- endBy parseExpr spaces
  tail <- char '.' >> spaces >> parseExpr
  return $ DottedList head tail

parseQuoted :: Parser LispVal
parseQuoted = do
  char '\''
  x <- parseExpr
  return $ List [Atom "quote", x]

--
-- Helpers
--
escapedChars :: Parser Char
escapedChars = do
  char '\\'
  c <- oneOf ['\\', '"', 'n', 'r', 't']
  return $
    case c of
      '\\' -> c
      '"' -> c
      'n' -> '\n'
      'r' -> '\r'
      't' -> '\t'

symbol :: Parser Char
symbol = oneOf "!$%&|*+-/:<=>?@^_~"

spaces :: Parser ()
spaces = skipMany1 space

bin2int :: String -> Integer
bin2int s = sum $ map (\(i, x) -> i * (2 ^ x)) $ zip [0 ..] $ map p (reverse s)
  where
    p '0' = 0
    p '1' = 1

readWith f s = fst $ f s !! 0

max_stars_repo_name: underscorenico/Scheme

max_issues_repo_name: underscorenico/Scheme

max_forks_repo_name: underscorenico/Scheme


text:
-- $ ghc --make -O -rtsopts -threaded parallel.hs
-- $ ./parallel 3000 +RTS -N4 -s -A200M

import System.Environment(getArgs)
import Numeric.LinearAlgebra
import Control.Parallel.Strategies
import System.Time

inParallel = parMap rwhnf id

-- matrix product decomposed into p parallel subtasks
parMul p x y = fromBlocks [ inParallel ( map (x <>) ys ) ]
    where [ys] = toBlocksEvery (rows y) (cols y `div` p) y

main = do
    n <- (read . head) `fmap` getArgs
    let m = ident n :: Matrix Double
    time $ print $ maxElement $ takeDiag $ m <> m
    time $ print $ maxElement $ takeDiag $ parMul 2 m m
    time $ print $ maxElement $ takeDiag $ parMul 4 m m
    time $ print $ maxElement $ takeDiag $ parMul 8 m m

time act = do
    t0 <- getClockTime
    act
    t1 <- getClockTime
    print $ tdSec $ normalizeTimeDiff $ diffClockTimes t1 t0

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
{-|
Module: Numeric.Morpheus.MatrixReduce
Description: Matrix reduce by column/row operations.
Copyright: (c) Alexander Ignatyev, 2017
License: BSD-3
Stability: experimental
Portability: POSIX
-}

module Numeric.Morpheus.MatrixReduce
(
  columnPredicate
  , rowPredicate
  , columnSum
  , rowSum
  , columnMaxIndex
  , columnMinIndex
  , rowMaxIndex
  , rowMinIndex
)

where

import Numeric.Morpheus.Utils(morpheusLayout)

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

type Predicate = R -> R -> R
foreign import ccall "wrapper" mkPredicate :: Predicate -> IO (FunPtr Predicate)

{- morpheus_column_predicate -}
foreign import ccall safe "morpheus_column_predicate"
  c_morpheus_column_predicate :: FunPtr Predicate -> CInt -> CInt -> CInt
                              -> Ptr R -> Ptr R -> IO ()

call_morpheus_column_predicate :: FunPtr Predicate
                               -> CInt -> CInt -> CInt -> CInt -> Ptr R
                               -> CInt -> Ptr R
                               -> IO ()
call_morpheus_column_predicate f rows cols xRow xCol matPtr _ vecPtr = do
  let layout = morpheusLayout xCol cols
  c_morpheus_column_predicate f layout rows cols matPtr vecPtr

-- | Scan every column of the given matrix.
-- Predicate takes and accumulator and next value of the column, returns new accumulator.
-- Returns accumulator values for every column.
columnPredicate :: (R -> R -> R) -> Matrix R -> Vector R
columnPredicate f m = unsafePerformIO $ do
  v <- createVector (cols m)
  fpred <- mkPredicate f
  apply m (apply v id) (call_morpheus_column_predicate fpred)
  return v

{- morpheus_row_predicate -}
foreign import ccall safe "morpheus_row_predicate"
  c_morpheus_row_predicate :: FunPtr Predicate -> CInt -> CInt -> CInt
                           -> Ptr R -> Ptr R -> IO ()

call_morpheus_row_predicate :: FunPtr Predicate
                            -> CInt -> CInt -> CInt -> CInt -> Ptr R
                            -> CInt -> Ptr R
                            -> IO ()
call_morpheus_row_predicate f rows cols xRow xCol matPtr _ vecPtr = do
  let layout = morpheusLayout xCol cols
  c_morpheus_row_predicate f layout rows cols matPtr vecPtr

-- | Scan every row of the given matrix.
-- Predicate takes and accumulator and next value of the row, returns new accumulator.
-- Returns accumulator values for every row.
rowPredicate :: (R -> R -> R) -> Matrix R -> Vector R
rowPredicate f m = unsafePerformIO $ do
  v <- createVector (rows m)
  fpred <- mkPredicate f
  apply m (apply v id) (call_morpheus_row_predicate fpred)
  return v

{- morpheus_column_sum -}
foreign import ccall unsafe "morpheus_column_sum"
  c_morpheus_column_sum :: CInt -> CInt -> CInt -> Ptr Double -> Ptr Double -> IO ()

call_morpheus_column_sum :: CInt -> CInt -> CInt -> CInt -> Ptr Double
                         -> CInt -> Ptr Double
                         -> IO ()
call_morpheus_column_sum rows cols xRow xCol matPtr _ vecPtr = do
  let layout = morpheusLayout xCol cols
  c_morpheus_column_sum layout rows cols matPtr vecPtr

-- | Calculates sums of elements of every column of the given matrix
columnSum :: Matrix Double -> Vector Double
columnSum m = unsafePerformIO $ do
    v <- createVector (cols m)
    apply m (apply v id) call_morpheus_column_sum
    return v

{- morpheus_row_sum -}
foreign import ccall unsafe "morpheus_row_sum"
  c_morpheus_row_sum :: CInt -> CInt -> CInt -> Ptr Double -> Ptr Double -> IO ()

call_morpheus_row_sum :: CInt -> CInt -> CInt -> CInt -> Ptr Double
                      -> CInt -> Ptr Double
                      -> IO ()
call_morpheus_row_sum rows cols xRow xCol matPtr _ vecPtr = do
  let layout = morpheusLayout xCol cols
  c_morpheus_row_sum layout rows cols matPtr vecPtr

-- | Calculates sums of elements of every row of the given matrix
rowSum :: Matrix Double -> Vector Double
rowSum m = unsafePerformIO $ do
    v <- createVector (rows m)
    apply m (apply v id) call_morpheus_row_sum
    return v

{- morpheus_column_max_index -}
foreign import ccall unsafe "morpheus_column_max_index"
  c_morpheus_column_max_index :: CInt -> CInt -> CInt -> Ptr Double -> Ptr Double -> Ptr CInt -> IO ()

call_morpheus_column_max_index :: CInt -> CInt -> CInt -> CInt -> Ptr Double
                                  -> CInt -> Ptr Double
                                  -> CInt -> Ptr CInt
                                  -> IO ()
call_morpheus_column_max_index rows cols xRow xCol matPtr _ vecPtr _ idxPtr = do
  let layout = morpheusLayout xCol cols
  c_morpheus_column_max_index layout rows cols matPtr vecPtr idxPtr

-- | Finds maximum values and their indices of every column of the given matrix
columnMaxIndex :: Matrix Double -> (Vector R, Vector I)
columnMaxIndex m = unsafePerformIO $ do
    v <- createVector (cols m)
    i <- createVector (cols m)
    apply m (apply v (apply i id)) call_morpheus_column_max_index
    return (v, i)

{- morpheus_column_min_index -}
foreign import ccall unsafe "morpheus_column_min_index"
  c_morpheus_column_min_index :: CInt -> CInt -> CInt -> Ptr Double -> Ptr Double -> Ptr CInt -> IO ()

call_morpheus_column_min_index :: CInt -> CInt -> CInt -> CInt -> Ptr Double
                                  -> CInt -> Ptr Double
                                  -> CInt -> Ptr CInt
                                  -> IO ()
call_morpheus_column_min_index rows cols xRow xCol matPtr _ vecPtr _ idxPtr = do
  let layout = morpheusLayout xCol cols
  c_morpheus_column_min_index layout rows cols matPtr vecPtr idxPtr

-- | Finds minimum values and their indices of every column of the given matrix
columnMinIndex :: Matrix Double -> (Vector R, Vector I)
columnMinIndex m = unsafePerformIO $ do
    v <- createVector (cols m)
    i <- createVector (cols m)
    apply m (apply v (apply i id)) call_morpheus_column_min_index
    return (v, i)

{- morpheus_row_max_index -}
foreign import ccall unsafe "morpheus_row_max_index"
  c_morpheus_row_max_index :: CInt -> CInt -> CInt -> Ptr Double -> Ptr Double -> Ptr CInt -> IO ()

call_morpheus_row_max_index :: CInt -> CInt -> CInt -> CInt -> Ptr Double
                               -> CInt -> Ptr Double
                               -> CInt -> Ptr CInt
                               -> IO ()
call_morpheus_row_max_index rows cols xRow xCol matPtr _ vecPtr _ idxPtr = do
  let layout = morpheusLayout xCol cols
  c_morpheus_row_max_index layout rows cols matPtr vecPtr idxPtr

-- | Finds maximum values and their indices of every row of the given matrix
rowMaxIndex :: Matrix Double -> (Vector R, Vector I)
rowMaxIndex m = unsafePerformIO $ do
    v <- createVector (rows m)
    i <- createVector (rows m)
    apply m (apply v (apply i id)) call_morpheus_row_max_index
    return (v, i)

{- morpheus_row_min_index -}
foreign import ccall unsafe "morpheus_row_min_index"
  c_morpheus_row_min_index :: CInt -> CInt -> CInt -> Ptr Double -> Ptr Double -> Ptr CInt -> IO ()

call_morpheus_row_min_index :: CInt -> CInt -> CInt -> CInt -> Ptr Double
                               -> CInt -> Ptr Double
                               -> CInt -> Ptr CInt
                               -> IO ()
call_morpheus_row_min_index rows cols xRow xCol matPtr _ vecPtr _ idxPtr = do
  let layout = morpheusLayout xCol cols
  c_morpheus_row_min_index layout rows cols matPtr vecPtr idxPtr

-- | Finds minimum values and their indices of every row of the given matrix
rowMinIndex :: Matrix Double -> (Vector R, Vector I)
rowMinIndex m = unsafePerformIO $ do
    v <- createVector (rows m)
    i <- createVector (rows m)
    apply m (apply v (apply i id)) call_morpheus_row_min_index
    return (v, i)

max_stars_repo_name: Alexander-Ignatyev/morpheus

max_issues_repo_name: aligusnet/morpheus

max_forks_repo_name: aligusnet/morpheus


text:
{-# LANGUAGE TemplateHaskell #-}
module Turtle (main) where
import Data.Complex
import Ros.Node
import Ros.Topic (cons, repeatM)
import Ros.TopicUtil (tee, filterBy, everyNew, interruptible, gate, share)
import Ros.Turtlesim.Pose
import Ros.Turtlesim.Velocity
import Ros.Logging
import System.IO (hFlush, stdout)

-- A type synonym for a 2D point.
type Point = Complex Float

-- A Topic of user-supplied waypoint trajectories.
getTraj :: Topic IO [Point]
getTraj = repeatM (do putStr "Enter waypoints: " >> hFlush stdout
                      $(logInfo "Waiting for new traj")
                      (map (uncurry (:+)) . read) `fmap` getLine)

wrapAngle theta
  | theta < 0    = theta + 2 * pi
  | theta > 2*pi = theta - 2 * pi
  | otherwise    = theta

angleDiff x y = wrapAngle (x' + pi - y') - pi
  where x' = if x < 0 then x + 2*pi else x
        y' = if y < 0 then y + 2*pi else y

-- Produce a unit value every time a goal is reached.
arrivalTrigger :: Topic IO Point -> Topic IO Pose -> Topic IO ()
arrivalTrigger goals poses = fmap (const ()) $
                             filterBy (fmap arrived goals) (fmap p2v poses)
  where arrived goal pose = magnitude (goal - pose) < 1.5
        p2v (Pose x y _ _ _) = x :+ y

-- Navigate to a goal given a current pose estimate.
navigate :: (Point, Pose) -> Velocity
navigate (goal, pos) = Velocity (min 2 (magnitude v)) angVel
  where v        = goal - (x pos :+ y pos)
        thetaErr = (angleDiff (phase v) (theta pos)) * (180 / pi)
        angVel   = signum thetaErr * (min 2 (abs thetaErr))

main = runNode "HaskellBTurtle" $
       do enableLogging (Just Warn)
          poses <- subscribe "/turtle1/pose"
          (t1,t2) <- liftIO . tee . interruptible $ getTraj
          let goals = gate t1 (cons () (arrivalTrigger t2 poses))
          advertise "/turtle1/command_velocity" $
                    fmap navigate (everyNew goals poses)

max_stars_repo_name: rgleichman/roshask

max_issues_repo_name: rgleichman/roshask

max_forks_repo_name: rgleichman/roshask