text:
module Feldspar.Frontend where

import Prelude (Integral, Ord, RealFloat, RealFrac)
import qualified Prelude as P
import Prelude.EDSL

import Control.Monad.Identity
import Data.Bits (Bits, FiniteBits)
import qualified Data.Bits as Bits
import Data.Complex (Complex)
import Data.Int
import Data.List (genericLength)

import Language.Syntactic (Internal)
import Language.Syntactic.Functional
import qualified Language.Syntactic as Syntactic

import qualified Control.Monad.Operational.Higher as Oper
import Language.Embedded.Imperative (IxRange)
import qualified Language.Embedded.Imperative as Imp

import qualified Data.Inhabited as Inhabited
import Data.TypedStruct
import Feldspar.Primitive.Representation
import Feldspar.Representation
import Feldspar.Sugar ()

--------------------------------------------------------------------------------
-- * Pure expressions
--------------------------------------------------------------------------------

----------------------------------------
-- ** General constructs
----------------------------------------

-- | Force evaluation of a value and share the result. Note that due to common
-- sub-expression elimination, this function is rarely needed in practice.
share :: (Syntax a, Syntax b)
    => a         -- ^ Value to share
    -> (a -> b)  -- ^ Body in which to share the value
    -> b
share = shareTag ""
  -- Explicit sharing can be useful e.g. when the value to share contains a
  -- function or when the code motion algorithm for some reason doesn't work
  -- find opportunities for sharing.

-- | Explicit tagged sharing
shareTag :: (Syntax a, Syntax b)
    => String    -- ^ A tag (that may be empty). May be used by a back end to
                 --   generate a sensible variable name.
    -> a         -- ^ Value to share
    -> (a -> b)  -- ^ Body in which to share the value
    -> b
shareTag tag = sugarSymFeld (Let tag)

-- | For loop
forLoop :: Syntax st => Data Length -> st -> (Data Index -> st -> st) -> st
forLoop = sugarSymFeld ForLoop

-- | Conditional expression
cond :: Syntax a
    => Data Bool  -- ^ Condition
    -> a          -- ^ True branch
    -> a          -- ^ False branch
    -> a
cond = sugarSymFeld Cond

-- | Condition operator; use as follows:
--
-- @
-- cond1 `?` a $
-- cond2 `?` b $
-- cond3 `?` c $
--         default
-- @
(?) :: Syntax a
    => Data Bool  -- ^ Condition
    -> a          -- ^ True branch
    -> a          -- ^ False branch
    -> a
(?) = cond

infixl 1 ?

-- | Multi-way conditional expression
--
-- The first association @(a,b)@ in the list of cases for which @a@ is equal to
-- the scrutinee is selected, and the associated @b@ is returned as the result.
-- If no case matches, the default value is returned.
switch :: (Syntax a, Syntax b, PrimType (Internal a))
    => b        -- ^ Default result
    -> [(a,b)]  -- ^ Cases (match, result)
    -> a        -- ^ Scrutinee
    -> b        -- ^ Result
switch def [] _ = def
switch def cs s = P.foldr
    (\(c,a) b -> desugar c == desugar s ? a $ b)
    def
    cs

----------------------------------------
-- ** Literals
----------------------------------------

-- | Literal
value :: Syntax a => Internal a -> a
value = sugarSymFeld . Lit

false :: Data Bool
false = value False

true :: Data Bool
true = value True

instance Syntactic.Syntactic ()
  where
    type Domain ()   = FeldDomain
    type Internal () = Int32
    desugar () = unData 0
    sugar   _  = ()

-- | Example value
--
-- 'example' can be used similarly to 'undefined' in normal Haskell, i.e. to
-- create an expression whose value is irrelevant.
--
-- Note that it is generally not possible to use 'undefined' in Feldspar
-- expressions, as this will crash the compiler.
example :: Syntax a => a
example = value Inhabited.example

----------------------------------------
-- ** Primitive functions
----------------------------------------

instance (Bounded a, Type a) => Bounded (Data a)
  where
    minBound = value minBound
    maxBound = value maxBound

instance (Num a, PrimType a) => Num (Data a)
  where
    fromInteger = value . fromInteger
    (+)         = sugarSymFeld Add
    (-)         = sugarSymFeld Sub
    (*)         = sugarSymFeld Mul
    negate      = sugarSymFeld Neg
    abs         = sugarSymFeld Abs
    signum      = sugarSymFeld Sign

instance (Fractional a, PrimType a) => Fractional (Data a)
  where
    fromRational = value . fromRational
    (/) = sugarSymFeld FDiv

instance (Floating a, PrimType a) => Floating (Data a)
  where
    pi    = sugarSymFeld Pi
    exp   = sugarSymFeld Exp
    log   = sugarSymFeld Log
    sqrt  = sugarSymFeld Sqrt
    (**)  = sugarSymFeld Pow
    sin   = sugarSymFeld Sin
    cos   = sugarSymFeld Cos
    tan   = sugarSymFeld Tan
    asin  = sugarSymFeld Asin
    acos  = sugarSymFeld Acos
    atan  = sugarSymFeld Atan
    sinh  = sugarSymFeld Sinh
    cosh  = sugarSymFeld Cosh
    tanh  = sugarSymFeld Tanh
    asinh = sugarSymFeld Asinh
    acosh = sugarSymFeld Acosh
    atanh = sugarSymFeld Atanh

-- | Alias for 'pi'
π :: (Floating a, PrimType a) => Data a
π = pi

-- | Integer division truncated toward zero
quot :: (Integral a, PrimType a) => Data a -> Data a -> Data a
quot = sugarSymFeld Quot

-- | Integer remainder satisfying
--
-- > (x `quot` y)*y + (x `rem` y) == x
rem :: (Integral a, PrimType a) => Data a -> Data a -> Data a
rem = sugarSymFeld Rem

-- | Simultaneous @quot@ and @rem@
quotRem :: (Integral a, PrimType a) => Data a -> Data a -> (Data a, Data a)
quotRem a b = (q,r)
  where
    q = quot a b
    r = a - b * q

-- | Integer division truncated toward negative infinity
div :: (Integral a, PrimType a) => Data a -> Data a -> Data a
div = sugarSymFeld Div

-- | Integer modulus, satisfying
--
-- > (x `div` y)*y + (x `mod` y) == x
mod :: (Integral a, PrimType a) => Data a -> Data a -> Data a
mod = sugarSymFeld Mod

-- | Integer division assuming `unsafeBalancedDiv x y * y == x` (i.e. no
-- remainder)
--
-- The advantage of using 'unsafeBalancedDiv' over 'quot' or 'div' is that the
-- above assumption can be used for simplifying the expression.
unsafeBalancedDiv :: (Integral a, PrimType a) => Data a -> Data a -> Data a
unsafeBalancedDiv a b = guardValLabel
    InternalAssertion
    (rem a b == 0)
    "unsafeBalancedDiv: division not balanced"
    (sugarSymFeld DivBalanced a b)
  -- Note: We can't check that `result * b == a`, because `result * b` gets
  -- simplified to `a`.

-- | Construct a complex number
complex :: (Num a, PrimType a, PrimType (Complex a))
    => Data a  -- ^ Real part
    -> Data a  -- ^ Imaginary part
    -> Data (Complex a)
complex = sugarSymFeld Complex

-- | Construct a complex number
polar :: (Floating a, PrimType a, PrimType (Complex a))
    => Data a  -- ^ Magnitude
    -> Data a  -- ^ Phase
    -> Data (Complex a)
polar = sugarSymFeld Polar

-- | Extract the real part of a complex number
realPart :: (PrimType a, PrimType (Complex a)) => Data (Complex a) -> Data a
realPart = sugarSymFeld Real

-- | Extract the imaginary part of a complex number
imagPart :: (PrimType a, PrimType (Complex a)) => Data (Complex a) -> Data a
imagPart = sugarSymFeld Imag

-- | Extract the magnitude of a complex number's polar form
magnitude :: (RealFloat a, PrimType a, PrimType (Complex a)) =>
    Data (Complex a) -> Data a
magnitude = sugarSymFeld Magnitude

-- | Extract the phase of a complex number's polar form
phase :: (RealFloat a, PrimType a, PrimType (Complex a)) =>
    Data (Complex a) -> Data a
phase = sugarSymFeld Phase

-- | Complex conjugate
conjugate :: (RealFloat a, PrimType (Complex a)) =>
    Data (Complex a) -> Data (Complex a)
conjugate = sugarSymFeld Conjugate
  -- `RealFloat` could be replaced by `Num` here, but it seems more consistent
  -- to use `RealFloat` for all functions.

-- | Integral type casting
i2n :: (Integral i, Num n, PrimType i, PrimType n) => Data i -> Data n
i2n = sugarSymFeld I2N

-- | Cast integer to 'Bool'
i2b :: (Integral a, PrimType a) => Data a -> Data Bool
i2b = sugarSymFeld I2B

-- | Cast 'Bool' to integer
b2i :: (Integral a, PrimType a) => Data Bool -> Data a
b2i = sugarSymFeld B2I

-- | Round a floating-point number to an integer
round :: (RealFrac a, Num b, PrimType a, PrimType b) => Data a -> Data b
round = sugarSymFeld Round

-- | Boolean negation
not :: Data Bool -> Data Bool
not = sugarSymFeld Not

-- | Boolean conjunction
(&&) :: Data Bool -> Data Bool -> Data Bool
(&&) = sugarSymFeld And

infixr 3 &&

-- | Boolean disjunction
(||) :: Data Bool -> Data Bool -> Data Bool
(||) = sugarSymFeld Or

infixr 2 ||

-- | Equality
(==) :: PrimType a => Data a -> Data a -> Data Bool
(==) = sugarSymFeld Eq

-- | Inequality
(/=) :: PrimType a => Data a -> Data a -> Data Bool
a /= b = not (a==b)

-- | Less than
(<) :: (Ord a, PrimType a) => Data a -> Data a -> Data Bool
(<) = sugarSymFeld Lt

-- | Greater than
(>) :: (Ord a, PrimType a) => Data a -> Data a -> Data Bool
(>) = sugarSymFeld Gt

-- | Less than or equal
(<=) :: (Ord a, PrimType a) => Data a -> Data a -> Data Bool
(<=) = sugarSymFeld Le

-- | Greater than or equal
(>=) :: (Ord a, PrimType a) => Data a -> Data a -> Data Bool
(>=) = sugarSymFeld Ge

infix 4 ==, /=, <, >, <=, >=

-- | Return the smallest of two values
min :: (Ord a, PrimType a) => Data a -> Data a -> Data a
min a b = a<=b ? a $ b
  -- There's no standard definition of min/max in C:
  -- <http://stackoverflow.com/questions/3437404/min-and-max-in-c>
  --
  -- There is `fmin`/`fminf` for floating-point numbers, but these are
  -- implemented essentially as above (except that they handle `NaN`
  -- specifically:
  -- <https://sourceware.org/git/?p=glibc.git;a=blob;f=math/s_fmin.c;hb=HEAD>

-- | Return the greatest of two values
max :: (Ord a, PrimType a) => Data a -> Data a -> Data a
max a b = a>=b ? a $ b

----------------------------------------
-- ** Bit manipulation
----------------------------------------

-- | Bit-wise \"and\"
(.&.) :: (Bits a, PrimType a) => Data a -> Data a -> Data a
(.&.) = sugarSymFeld BitAnd

-- | Bit-wise \"or\"
(.|.) :: (Bits a, PrimType a) => Data a -> Data a -> Data a
(.|.) = sugarSymFeld BitOr

-- | Bit-wise \"xor\"
xor :: (Bits a, PrimType a) => Data a -> Data a -> Data a
xor = sugarSymFeld BitXor

-- | Bit-wise \"xor\"
(⊕) :: (Bits a, PrimType a) => Data a -> Data a -> Data a
(⊕) = xor

-- | Bit-wise complement
complement :: (Bits a, PrimType a) => Data a -> Data a
complement = sugarSymFeld BitCompl

-- | Left shift
shiftL :: (Bits a, PrimType a)
    => Data a      -- ^ Value to shift
    -> Data Int32  -- ^ Shift amount (negative value gives right shift)
    -> Data a
shiftL = sugarSymFeld ShiftL

-- | Right shift
shiftR :: (Bits a, PrimType a)
    => Data a      -- ^ Value to shift
    -> Data Int32  -- ^ Shift amount (negative value gives left shift)
    -> Data a
shiftR = sugarSymFeld ShiftR

-- | Left shift
(.<<.) :: (Bits a, PrimType a)
    => Data a      -- ^ Value to shift
    -> Data Int32  -- ^ Shift amount (negative value gives right shift)
    -> Data a
(.<<.) = shiftL

-- | Right shift
(.>>.) :: (Bits a, PrimType a)
    => Data a      -- ^ Value to shift
    -> Data Int32  -- ^ Shift amount (negative value gives left shift)
    -> Data a
(.>>.) = shiftR

infixl 8 `shiftL`, `shiftR`, .<<., .>>.
infixl 7 .&.
infixl 6 `xor`
infixl 5 .|.

bitSize :: forall a . FiniteBits a => Data a -> Length
bitSize _ = P.fromIntegral $ Bits.finiteBitSize (a :: a)
  where
    a = P.error "finiteBitSize evaluates its argument"

-- | Set all bits to one
allOnes :: (Bits a, Num a, PrimType a) => Data a
allOnes = complement 0

-- | Set the @n@ lowest bits to one
oneBits :: (Bits a, Num a, PrimType a) => Data Int32 -> Data a
oneBits n = complement (allOnes .<<. n)

-- | Extract the @k@ lowest bits
lsbs :: (Bits a, Num a, PrimType a) => Data Int32 -> Data a -> Data a
lsbs k i = i .&. oneBits k

-- | Integer logarithm in base 2. Returns \(\lfloor log_2(x) \rfloor\).
-- Assumes \(x>0\).
ilog2 :: (FiniteBits a, Integral a, PrimType a) => Data a -> Data a
ilog2 a = guardValLabel InternalAssertion (a >= 1) "ilog2: argument < 1" $
    snd $ P.foldr (\ffi vr -> share vr (step ffi)) (a,0) ffis
  where
    step (ff,i) (v,r) =
        share (b2i (v > fromInteger ff) .<<. value i) $ \shift ->
          (v .>>. i2n shift, r .|. shift)

    -- [(0x1, 0), (0x3, 1), (0xF, 2), (0xFF, 3), (0xFFFF, 4), ...]
    ffis
        = (`P.zip` [0..])
        $ P.takeWhile (P.<= (2 P.^ (bitSize a `P.div` 2) - 1 :: Integer))
        $ P.map ((subtract 1) . (2 P.^) . (2 P.^))
        $ [(0::Integer)..]
  -- Based on this algorithm:
  -- <http://graphics.stanford.edu/~seander/bithacks.html#IntegerLog>

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

-- | Index into an array
arrIx :: Syntax a => IArr a -> Data Index -> a
arrIx arr i = resugar $ mapStruct ix $ unIArr arr
  where
    ix :: forall b . PrimType' b => Imp.IArr Index b -> Data b
    ix arr' = sugarSymFeldPrim
      (GuardVal InternalAssertion "arrIx: index out of bounds")
      (i < length arr)
      (sugarSymFeldPrim (ArrIx arr') (i + iarrOffset arr) :: Data b)

class Indexed a
  where
    type IndexedElem a

    -- | Indexing operator. If @a@ is 'Finite', it is assumed that
    -- @i < `length` a@ in any expression @a `!` i@.
    (!) :: a -> Data Index -> IndexedElem a

infixl 9 !

-- | Linear structures with a length. If the type is also 'Indexed', the length
-- is the successor of the maximal allowed index.
class Finite a
  where
    -- | The length of a finite structure
    length :: a -> Data Length

instance Finite (Arr a)  where length = arrLength
instance Finite (IArr a) where length = iarrLength

-- | Linear structures that can be sliced
class Slicable a
  where
    -- | Take a slice of a structure
    slice
        :: Data Index   -- ^ Start index
        -> Data Length  -- ^ Slice length
        -> a            -- ^ Structure to slice
        -> a

instance Syntax a => Indexed (IArr a)
  where
    type IndexedElem (IArr a) = a
    (!) = arrIx

instance Slicable (Arr a)
  where
    slice from len (Arr o l arr) = Arr o' l' arr
      where
        o' = guardValLabel InternalAssertion (from<=l) "invalid Arr slice" (o+from)
        l' = guardValLabel InternalAssertion (from+len<=l) "invalid Arr slice" len

instance Slicable (IArr a)
  where
    slice from len (IArr o l arr) = IArr o' l' arr
      where
        o' = guardValLabel InternalAssertion (from<=l) "invalid IArr slice" (o+from)
        l' = guardValLabel InternalAssertion (from+len<=l) "invalid IArr slice" len

----------------------------------------
-- ** Syntactic conversion
----------------------------------------

desugar :: Syntax a => a -> Data (Internal a)
desugar = Data . Syntactic.desugar

sugar :: Syntax a => Data (Internal a) -> a
sugar = Syntactic.sugar . unData

-- | Cast between two values that have the same syntactic representation
resugar :: (Syntax a, Syntax b, Internal a ~ Internal b) => a -> b
resugar = Syntactic.resugar

----------------------------------------
-- ** Assertions
----------------------------------------

-- | Guard a value by an assertion (with implicit label @`UserAssertion` ""@)
guardVal :: Syntax a
    => Data Bool  -- ^ Condition that is expected to be true
    -> String     -- ^ Error message
    -> a          -- ^ Value to attach the assertion to
    -> a
guardVal = guardValLabel $ UserAssertion ""

-- | Like 'guardVal' but with an explicit assertion label
guardValLabel :: Syntax a
    => AssertionLabel  -- ^ Assertion label
    -> Data Bool       -- ^ Condition that is expected to be true
    -> String          -- ^ Error message
    -> a               -- ^ Value to attach the assertion to
    -> a
guardValLabel c cond msg = sugarSymFeld (GuardVal c msg) cond

----------------------------------------
-- ** Unsafe operations
----------------------------------------

-- | Turn a 'Comp' computation into a pure value. For this to be safe, the
-- computation should be free of side effects and independent of its
-- environment.
unsafePerform :: Syntax a => Comp a -> a
unsafePerform = sugarSymFeld . UnsafePerform . fmap desugar

--------------------------------------------------------------------------------
-- * Programs with computational effects
--------------------------------------------------------------------------------

-- | Monads that support computational effects: mutable data structures and
-- control flow
class Monad m => MonadComp m
  where
    -- | Lift a 'Comp' computation
    liftComp :: Comp a -> m a
    -- | Conditional statement
    iff :: Data Bool -> m () -> m () -> m ()
    -- | For loop
    for :: (Integral n, PrimType n) =>
        IxRange (Data n) -> (Data n -> m ()) -> m ()
    -- | While loop
    while :: m (Data Bool) -> m () -> m ()

instance MonadComp Comp
  where
    liftComp        = id
    iff c t f       = Comp $ Imp.iff c (unComp t) (unComp f)
    for range body  = Comp $ Imp.for range (unComp . body)
    while cont body = Comp $ Imp.while (unComp cont) (unComp body)

----------------------------------------
-- ** References
----------------------------------------

-- | Create an uninitialized reference
newRef :: (Syntax a, MonadComp m) => m (Ref a)
newRef = newNamedRef "r"

-- | Create an uninitialized named reference
--
-- The provided base name may be appended with a unique identifier to avoid name
-- collisions.
newNamedRef :: (Syntax a, MonadComp m)
    => String  -- ^ Base name
    -> m (Ref a)
newNamedRef base = liftComp $ fmap Ref $
    mapStructA (const $ Comp $ Imp.newNamedRef base) typeRep

-- | Create an initialized named reference
initRef :: (Syntax a, MonadComp m) => a -> m (Ref a)
initRef = initNamedRef "r"

-- | Create an initialized reference
--
-- The provided base name may be appended with a unique identifier to avoid name
-- collisions.
initNamedRef :: (Syntax a, MonadComp m)
    => String  -- ^ Base name
    -> a       -- ^ Initial value
    -> m (Ref a)
initNamedRef base =
    liftComp . fmap Ref . mapStructA (Comp . Imp.initNamedRef base) . resugar

-- | Get the contents of a reference.
getRef :: (Syntax a, MonadComp m) => Ref a -> m a
getRef = liftComp . fmap resugar . mapStructA (Comp . Imp.getRef) . unRef

-- | Set the contents of a reference.
setRef :: (Syntax a, MonadComp m) => Ref a -> a -> m ()
setRef r
    = liftComp
    . sequence_
    . zipListStruct (\r' a' -> Comp $ Imp.setRef r' a') (unRef r)
    . resugar

-- | Modify the contents of reference.
modifyRef :: (Syntax a, MonadComp m) => Ref a -> (a -> a) -> m ()
modifyRef r f = setRef r . f =<< unsafeFreezeRef r

-- | Freeze the contents of reference (only safe if the reference is not updated
--   as long as the resulting value is alive).
unsafeFreezeRef :: (Syntax a, MonadComp m) => Ref a -> m a
unsafeFreezeRef
    = liftComp
    . fmap resugar
    . mapStructA (Comp . Imp.unsafeFreezeRef)
    . unRef

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

-- | Create an uninitialized array
newArr :: (Type (Internal a), MonadComp m) => Data Length -> m (Arr a)
newArr = newNamedArr "a"

-- | Create an uninitialized named array
--
-- The provided base name may be appended with a unique identifier to avoid name
-- collisions.
newNamedArr :: (Type (Internal a), MonadComp m)
    => String  -- ^ Base name
    -> Data Length
    -> m (Arr a)
newNamedArr base l = liftComp $ fmap (Arr 0 l) $
    mapStructA (const (Comp $ Imp.newNamedArr base l)) typeRep

-- | Create an array and initialize it with a constant list
constArr :: (PrimType (Internal a), MonadComp m)
    => [Internal a]  -- ^ Initial contents
    -> m (Arr a)
constArr = constNamedArr "a"

-- | Create a named array and initialize it with a constant list
--
-- The provided base name may be appended with a unique identifier to avoid name
-- collisions.
constNamedArr :: (PrimType (Internal a), MonadComp m)
    => String        -- ^ Base name
    -> [Internal a]  -- ^ Initial contents
    -> m (Arr a)
constNamedArr base as =
    liftComp $ fmap (Arr 0 len . Single) $ Comp $ Imp.constNamedArr base as
  where
    len = value $ genericLength as

-- | Get an element of an array
getArr :: (Syntax a, MonadComp m) => Arr a -> Data Index -> m a
getArr arr i = do
    assertLabel
      InternalAssertion
      (i < length arr)
      "getArr: index out of bounds"
    liftComp
      $ fmap resugar
      $ mapStructA (Comp . flip Imp.getArr (i + arrOffset arr))
      $ unArr arr

-- | Set an element of an array
setArr :: forall m a . (Syntax a, MonadComp m) =>
    Arr a -> Data Index -> a -> m ()
setArr arr i a = do
    assertLabel
      InternalAssertion
      (i < length arr)
      "setArr: index out of bounds"
    liftComp
      $ sequence_
      $ zipListStruct
          (\a' arr' -> Comp $ Imp.setArr arr' (i + arrOffset arr) a') rep
      $ unArr arr
  where
    rep = resugar a :: Struct PrimType' Data (Internal a)

-- | Copy the contents of an array to another array. The length of the
-- destination array must not be less than that of the source array.
--
-- In order to copy only a part of an array, use 'slice' before calling
-- 'copyArr'.
copyArr :: MonadComp m
    => Arr a  -- ^ Destination
    -> Arr a  -- ^ Source
    -> m ()
copyArr arr1 arr2 = do
    assertLabel
      InternalAssertion
      (length arr1 >= length arr2)
      "copyArr: destination too small"
    liftComp $ sequence_ $
      zipListStruct
        (\a1 a2 ->
            Comp $ Imp.copyArr
              (a1, arrOffset arr1)
              (a2, arrOffset arr2)
              (length arr2)
        )
        (unArr arr1)
        (unArr arr2)

-- | Freeze a mutable array to an immutable one. This involves copying the array
-- to a newly allocated one.
freezeArr :: (Type (Internal a), MonadComp m) => Arr a -> m (IArr a)
freezeArr arr = liftComp $ do
    arr2 <- newArr (length arr)
    copyArr arr2 arr
    unsafeFreezeArr arr2
  -- This is better than calling `freezeArr` from imperative-edsl, since that
  -- one copies without offset.

-- | A version of 'freezeArr' that slices the array from 0 to the given length
freezeSlice :: (Type (Internal a), MonadComp m) =>
    Data Length -> Arr a -> m (IArr a)
freezeSlice len = fmap (slice 0 len) . freezeArr

-- | Freeze a mutable array to an immutable one without making a copy. This is
-- generally only safe if the the mutable array is not updated as long as the
-- immutable array is alive.
unsafeFreezeArr :: MonadComp m => Arr a -> m (IArr a)
unsafeFreezeArr arr
    = liftComp
    $ fmap (IArr (arrOffset arr) (length arr))
    $ mapStructA (Comp . Imp.unsafeFreezeArr)
    $ unArr arr

-- | A version of 'unsafeFreezeArr' that slices the array from 0 to the given
-- length
unsafeFreezeSlice :: MonadComp m => Data Length -> Arr a -> m (IArr a)
unsafeFreezeSlice len = fmap (slice 0 len) . unsafeFreezeArr

-- | Thaw an immutable array to a mutable one. This involves copying the array
-- to a newly allocated one.
thawArr :: (Type (Internal a), MonadComp m) => IArr a -> m (Arr a)
thawArr arr = liftComp $ do
    arr2 <- unsafeThawArr arr
    arr3 <- newArr (length arr)
    copyArr arr3 arr2
    return arr3

-- | Thaw an immutable array to a mutable one without making a copy. This is
-- generally only safe if the the mutable array is not updated as long as the
-- immutable array is alive.
unsafeThawArr :: MonadComp m => IArr a -> m (Arr a)
unsafeThawArr arr
    = liftComp
    $ fmap (Arr (iarrOffset arr) (length arr))
    $ mapStructA (Comp . Imp.unsafeThawArr)
    $ unIArr arr

-- | Create an immutable array and initialize it with a constant list
constIArr :: (PrimType (Internal a), MonadComp m) =>
    [Internal a] -> m (IArr a)
constIArr = constArr >=> unsafeFreezeArr

----------------------------------------
-- ** Control-flow
----------------------------------------

-- | Conditional statement that returns an expression
ifE :: (Syntax a, MonadComp m)
    => Data Bool  -- ^ Condition
    -> m a        -- ^ True branch
    -> m a        -- ^ False branch
    -> m a
ifE c t f = do
    res <- newRef
    iff c (t >>= setRef res) (f >>= setRef res)
    unsafeFreezeRef res

-- | Break out from a loop
break :: MonadComp m => m ()
break = liftComp $ Comp Imp.break

-- | Assertion (with implicit label @`UserAssertion` ""@)
assert :: MonadComp m
    => Data Bool  -- ^ Expression that should be true
    -> String     -- ^ Message in case of failure
    -> m ()
assert = assertLabel $ UserAssertion ""

-- | Like 'assert' but tagged with an explicit assertion label
assertLabel :: MonadComp m
    => AssertionLabel  -- ^ Assertion label
    -> Data Bool       -- ^ Expression that should be true
    -> String          -- ^ Message in case of failure
    -> m ()
assertLabel c cond msg =
    liftComp $ Comp $ Oper.singleInj $ Assert c cond msg

----------------------------------------
-- ** Misc.
----------------------------------------

-- | Force evaluation of a value and share the result (monadic version of
-- 'share')
shareM :: (Syntax a, MonadComp m) => a -> m a
shareM = initRef >=> unsafeFreezeRef
  -- This function is more commonly needed than `share`, since code motion
  -- doesn't work across monadic binds.

max_stars_repo_name: Abhiroop/mu-feldspar

max_issues_repo_name: Abhiroop/mu-feldspar

max_forks_repo_name: Abhiroop/mu-feldspar


text:
module Main where

import Data.Monoid

import Core
import Render
import Codec.Picture
import Codec.Picture.Types
import Numeric.LinearAlgebra(vector)

settings :: RenderSettings
settings = RenderSettings {
    background = PixelRGBF 1 1 1,
    width = 800,
    height = 800,
    path = "out.png",
    topLeft = vector [1,1,0],
    topRight = vector [1,-1,0],
    bottomRight = vector [-1, -1,0],
    origin = vector [0,0,-1],
    antialiasing = 10}

scene
    = plane (diffuse (scaleColor 0.5)) (Ray (vec3 (-2) 0 0) (vec3 (-2) 0 0))
    <> sphere 0.9 (vec3 0 0 2) (diffuse (scaleColor 0.5))
    <> plane (skylike white blue) (Ray (vec3 10 0 0) (vec3 1 0 0))

main :: IO ()
main = do
    img <- render settings scene
    savePngImage (path settings) img

ext: hs

max_stars_repo_name: danrocag/raytracer

max_issues_repo_name: danrocag/raytracer

max_forks_repo_name: danrocag/raytracer


text:
{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, ScopedTypeVariables, RecordWildCards #-}

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

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

import           Foreign.Storable.Complex
import           Criterion.Main
import           Test.QuickCheck
import           Test.QuickCheck.Monadic

-- | A class for things that can be multiplied by a scalar.
class Mult a b where
    mult :: a -> b -> a

instance (Num a) => Mult a a where
    mult = (*)

instance (Num a) => Mult (Complex a) a where
    mult (x :+ y) z = (x * z) :+ (y * z)

-- | Fill a mutable vector from a monadic stream
{-# INLINE fill #-}
fill :: (PrimMonad m, Functor m, VGM.MVector vm a) => VFS.MStream m a -> vm (PrimState m) a -> m ()
fill str outBuf = void $ VFSM.foldM' put 0 str
    where
    put i x = do
        VGM.unsafeWrite outBuf i x
        return $ i + 1

{-# INLINE stride #-}
stride :: VG.Vector v a => Int -> v a -> v a
stride str inv = VG.unstream $ VFS.unfoldr func 0
    where
    len = VG.length inv
    func i | i >= len  = Nothing
           | otherwise = Just (VG.unsafeIndex inv i, i + str)

-- | The functions to be benchmarked

-- | Filters

filterHighLevel :: (PrimMonad m, Functor m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) => Int -> v b -> v a -> vm (PrimState m) a -> m ()
filterHighLevel num coeffs inBuf outBuf = fill (VFSM.generate num dotProd) outBuf
    where
    dotProd offset = VG.sum $ VG.zipWith mult (VG.unsafeDrop offset inBuf) coeffs

type FilterCRR = CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()
type FilterRR  = Int -> VS.Vector Float -> VS.Vector Float -> VS.MVector RealWorld Float -> IO ()
type FilterRC  = Int -> VS.Vector Float -> VS.Vector (Complex Float) -> VS.MVector RealWorld (Complex Float) -> IO ()

filterFFIR :: FilterCRR -> FilterRR
filterFFIR func num coeffs inBuf outBuf =
    VS.unsafeWith (unsafeCoerce coeffs) $ \cPtr ->
        VS.unsafeWith (unsafeCoerce inBuf) $ \iPtr ->
            VSM.unsafeWith (unsafeCoerce outBuf) $ \oPtr ->
                func (fromIntegral num) (fromIntegral $ VG.length coeffs) cPtr iPtr oPtr

filterFFIC :: FilterCRR -> FilterRC
filterFFIC func num coeffs inBuf outBuf =
    VS.unsafeWith (unsafeCoerce coeffs) $ \cPtr ->
        VS.unsafeWith (unsafeCoerce inBuf) $ \iPtr ->
            VSM.unsafeWith (unsafeCoerce outBuf) $ \oPtr ->
                func (fromIntegral num) (fromIntegral $ VG.length coeffs) cPtr iPtr oPtr

-- | Decimation

decimateHighLevel :: (PrimMonad m, Functor m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) => Int -> Int -> v b -> v a -> vm (PrimState m) a -> m ()
decimateHighLevel num factor coeffs inBuf outBuf = fill x outBuf
    where
    x = VFSM.map dotProd (VFSM.iterateN num (+ factor) 0)
    dotProd offset = VG.sum $ VG.zipWith mult (VG.unsafeDrop offset inBuf) coeffs

type DecimateCRR = CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()
type DecimateRR  = Int -> Int -> VS.Vector Float -> VS.Vector Float -> VS.MVector RealWorld Float -> IO ()
type DecimateRC  = Int -> Int -> VS.Vector Float -> VS.Vector (Complex Float) -> VS.MVector RealWorld (Complex Float) -> IO ()

decimateFFIR :: DecimateCRR -> DecimateRR
decimateFFIR func num factor coeffs inBuf outBuf =
    VS.unsafeWith (unsafeCoerce coeffs) $ \cPtr ->
        VS.unsafeWith (unsafeCoerce inBuf) $ \iPtr ->
            VSM.unsafeWith (unsafeCoerce outBuf) $ \oPtr ->
                func (fromIntegral num) (fromIntegral factor) (fromIntegral $ VG.length coeffs) cPtr iPtr oPtr

decimateFFIC :: DecimateCRR -> DecimateRC
decimateFFIC func num factor coeffs inBuf outBuf =
    VS.unsafeWith (unsafeCoerce coeffs) $ \cPtr ->
        VS.unsafeWith (unsafeCoerce inBuf) $ \iPtr ->
            VSM.unsafeWith (unsafeCoerce outBuf) $ \oPtr ->
                func (fromIntegral num) (fromIntegral factor) (fromIntegral $ VG.length coeffs) cPtr iPtr oPtr

-- | Rational downsampling
resampleHighLevel :: (PrimMonad m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) => Int -> Int -> Int -> Int -> v b -> v a -> vm (PrimState m) a -> m Int
resampleHighLevel count interpolation decimation filterOffset coeffs inBuf outBuf = fill 0 filterOffset 0
    where
    fill i filterOffset inputOffset
        | i < count = do
            let dp = dotProd filterOffset inputOffset
            VGM.unsafeWrite outBuf i dp
            let (q, r)        = quotRem (decimation - filterOffset - 1) interpolation
                inputOffset'  = inputOffset + q + 1
                filterOffset' = interpolation - 1 - r
            filterOffset' `seq` inputOffset' `seq` fill (i + 1) filterOffset' inputOffset'
        | otherwise = return filterOffset
    dotProd filterOffset offset = VG.sum $ VG.zipWith mult (VG.unsafeDrop offset inBuf) (stride interpolation (VG.unsafeDrop filterOffset coeffs))

pad :: a -> Int -> [a] -> [a]
pad with num list = list ++ replicate (num - length list) with

strideList :: Int -> [a] -> [a]
strideList s xs = go 0 xs
    where
    go _ []     = []
    go 0 (x:xs) = x : go (s-1) xs
    go n (x:xs) = go (n - 1) xs

roundUp :: Int -> Int -> Int
roundUp num div = ((num + div - 1) `quot` div) * div

data Coeffs = Coeffs {
    numCoeffs  :: Int,
    numGroups  :: Int,
    increments :: [Int],
    groups     :: [[Float]]
}

prepareCoeffs :: Int -> Int -> Int -> [Float] -> Coeffs
prepareCoeffs n interpolation decimation coeffs = Coeffs {..}
    where
    numCoeffs   = maximum $ map (length . snd) dats
    numGroups   = length groups
    increments  = map fst dats

    groups      :: [[Float]]
    groups      = map (pad 0 (roundUp numCoeffs n)) $ map snd dats

    dats :: [(Int, [Float])]
    dats = func 0
        where

        func' 0      = []
        func' x      = func x

        func :: Int -> [(Int, [Float])]
        func offset = (increment, strideList interpolation $ drop offset coeffs) : func' offset'
            where
            (q, r)    = quotRem (decimation - offset - 1) interpolation
            increment = q + 1
            offset'   = interpolation - 1 - r

resampleFFIR :: (Ptr CFloat -> Ptr CFloat -> IO ()) -> VS.Vector Float -> VSM.MVector RealWorld Float -> IO ()
resampleFFIR func inBuf outBuf =
    VS.unsafeWith (unsafeCoerce inBuf) $ \iPtr ->
        VS.unsafeWith (unsafeCoerce outBuf) $ \oPtr ->
            func iPtr oPtr

type ResampleR = CInt -> CInt -> CInt -> CInt -> Ptr CInt -> Ptr (Ptr CFloat) -> Ptr CFloat -> Ptr CFloat -> IO ()

mkResampler :: ResampleR -> Int -> Int -> Int -> [Float] -> IO (Int -> Int -> VS.Vector Float -> VS.MVector RealWorld Float -> IO ())
mkResampler func n interpolation decimation coeffs = do
    groupsP     <- mapM newArray $ map (map realToFrac) groups
    groupsPP    <- newArray groupsP
    incrementsP <- newArray $ map fromIntegral increments
    return $ \num offset -> resampleFFIR $ func (fromIntegral num) (fromIntegral numCoeffs) (fromIntegral offset) (fromIntegral numGroups) incrementsP groupsPP
    where
    Coeffs {..} = prepareCoeffs n interpolation decimation coeffs

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

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

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

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

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

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

    --Benchmarks
    defaultMain [
            bgroup "filter" [
                bgroup "real" [
                    bench "highLevel"   $ nfIO $ filterHighLevel        num coeffs inBuf outBuf
                ],
                bgroup "complex" [
                    bench "highLevel"   $ nfIO $ filterHighLevel        num coeffs  inBufComplex outBufComplex
                ]
            ],
            bgroup "decimate" [
                bgroup "real" [
                    bench "highLevel"   $ nfIO $ decimateHighLevel        (num `quot` decimation) decimation coeffs inBuf outBuf
                ],
                bgroup "complex" [
                    bench "highLevel"   $ nfIO $ decimateHighLevel      (num `quot` decimation) decimation coeffs  inBufComplex outBufComplex
                ]
            ],
            bgroup "resample" [
                bgroup "real" [
                    bench "highLevel"   $ nfIO $ resampleHighLevel      (num `quot` decimation) interpolation decimation 0 coeffs inBuf outBuf
                ],
                bgroup "complex" [
                    bench "highLevel"   $ nfIO $ resampleHighLevel      (num `quot` decimation) interpolation decimation 0 coeffs inBufComplex outBufComplex
                ]
            ]
        ]

theTest = quickCheck $ conjoin [propFiltersComplex]
    where
    sizes           = elements [1024, 2048, 4096, 8192, 16384, 32768, 65536]
    numCoeffs       = elements [32, 64, 128, 256, 512]
    factors         = elements [1, 2, 3, 4, 7, 9, 12, 15, 21]
    factors'        = [1, 2, 3, 4, 7, 9, 12, 15, 21]
    propFiltersReal = forAll sizes $ \size ->
                          forAll (vectorOf size (choose (-10, 10))) $ \inBuf ->
                              forAll numCoeffs $ \numCoeffs ->
                                  forAll (vectorOf numCoeffs (choose (-10, 10))) $ \coeffs ->
                                      testFiltersReal size numCoeffs coeffs inBuf
    testFiltersReal :: Int -> Int -> [Float] -> [Float] -> Property
    testFiltersReal size numCoeffs coeffs inBuf = monadicIO $ do
        let vCoeffsHalf = VS.fromList coeffs
            vCoeffs     = VS.fromList $ coeffs ++ reverse coeffs
            vInput      = VS.fromList inBuf
            num         = size - numCoeffs*2 + 1

        r1 <- run $ getResult num $ filterHighLevel       num vCoeffs     vInput

    propFiltersComplex = forAll sizes $ \size ->
                             forAll (vectorOf size (choose (-10, 10))) $ \inBufR ->
                                 forAll (vectorOf size (choose (-10, 10))) $ \inBufI ->
                                     forAll numCoeffs $ \numCoeffs ->
                                         forAll (vectorOf numCoeffs (choose (-10, 10))) $ \coeffs ->
        assert $ all (r1 `eqDelta`) [r1]
                                             testFiltersComplex size numCoeffs coeffs $ zipWith (:+) inBufR inBufI
    testFiltersComplex :: Int -> Int -> [Float] -> [Complex Float] -> Property
    testFiltersComplex size numCoeffs coeffs inBuf = monadicIO $ do
        let vCoeffsHalf = VS.fromList coeffs
            vCoeffs     = VS.fromList $ coeffs ++ reverse coeffs
            vInput      = VS.fromList inBuf
            num         = size - numCoeffs*2 + 1
            --vCoeffs2    = VG.fromList $ duplicate $ coeffs ++ reverse coeffs

        r1 <- run $ getResult num $ filterHighLevel       num vCoeffs     vInput

    propDecimationReal = forAll sizes $ \size ->
                             forAll (vectorOf size (choose (-10, 10))) $ \inBuf ->
                                 forAll numCoeffs $ \numCoeffs ->
                                     forAll (vectorOf numCoeffs (choose (-10, 10))) $ \coeffs ->
                                        forAll factors $ \factor ->
        assert $ all (r1 `eqDeltaC`) [r1]
                                             testDecimationReal size numCoeffs factor coeffs inBuf
    testDecimationReal :: Int -> Int -> Int -> [Float] -> [Float] -> Property
    testDecimationReal size numCoeffs factor coeffs inBuf = monadicIO $ do
        let vCoeffsHalf = VS.fromList coeffs
            vCoeffs     = VS.fromList $ coeffs ++ reverse coeffs
            vInput      = VS.fromList inBuf
            num         = (size - numCoeffs*2 + 1) `quot` factor

        r1 <- run $ getResult num $ decimateHighLevel       num factor vCoeffs     vInput

    propDecimationComplex = forAll sizes $ \size ->
                                forAll (vectorOf size (choose (-10, 10))) $ \inBufR ->
                                    forAll (vectorOf size (choose (-10, 10))) $ \inBufI ->
                                        forAll numCoeffs $ \numCoeffs ->
                                            forAll (vectorOf numCoeffs (choose (-10, 10))) $ \coeffs ->
                                                forAll factors $ \factor ->
        assert $ all (r1 `eqDelta`) [r1]
                                                    testDecimationComplex size numCoeffs factor coeffs $ zipWith (:+) inBufR inBufI
    testDecimationComplex :: Int -> Int -> Int -> [Float] -> [Complex Float] -> Property
    testDecimationComplex size numCoeffs factor coeffs inBuf = monadicIO $ do
        let vCoeffsHalf = VS.fromList coeffs
            vCoeffs     = VS.fromList $ coeffs ++ reverse coeffs
            vInput      = VS.fromList inBuf
            num         = (size - numCoeffs*2 + 1) `quot` factor
            --vCoeffs2    = VG.fromList $ duplicate $ coeffs ++ reverse coeffs

        r1 <- run $ getResult num $ decimateHighLevel       num factor vCoeffs     vInput

    propResamplingReal = forAll sizes $ \size ->
                             forAll (vectorOf size (choose (-10, 10))) $ \inBuf ->
                                 forAll numCoeffs $ \numCoeffs ->
                                     forAll (vectorOf numCoeffs (choose (-10, 10))) $ \coeffs ->
                                        forAll (elements $ tail factors') $ \decimation ->
                                            forAll (elements $ filter (< decimation) factors') $ \interpolation ->
        assert $ all (r1 `eqDeltaC`) [r1]
                                                 testResamplingReal size numCoeffs interpolation decimation coeffs inBuf
    testResamplingReal :: Int -> Int -> Int -> Int -> [Float] -> [Float] -> Property
    testResamplingReal size numCoeffs interpolation decimation coeffs inBuf = monadicIO $ do
        let vCoeffsHalf = VS.fromList coeffs
            vCoeffs     = VS.fromList $ coeffs ++ reverse coeffs
            vInput      = VS.fromList inBuf
            num         = (size - numCoeffs*2 + 1) `quot` decimation

        r1 <- run $ getResult num $ resampleHighLevel       num interpolation decimation 0 vCoeffs vInput

        assert $ all (r1 `eqDelta`) [r1]
    getResult :: (VSM.Storable a) => Int -> (VS.MVector RealWorld a -> IO b) -> IO [a]
    getResult size func = do
        outBuf <- VGM.new size
        func outBuf
        out :: VS.Vector a <- VG.freeze outBuf
        return $ VG.toList out
    eqDelta x y = all (uncurry eqDelta') $ zip x y
        where
        eqDelta' x y = abs (x - y) < 0.01
    eqDeltaC x y = all (uncurry eqDelta') $ zip x y
        where
        eqDelta' x y = magnitude (x - y) < 0.01
    duplicate :: [a] -> [a]
    duplicate = concatMap func
        where func x = [x, x]

main = theBench

max_stars_repo_name: adamwalker/sdr

max_issues_repo_name: peixian/sdr

max_forks_repo_name: adamwalker/sdr


text:
module HollowPinwheel where

import           Control.Monad                  as M
import qualified Data.Array.Accelerate          as A
import           Data.Array.Accelerate.LLVM.PTX
import           Data.Array.IArray              as IA
import           Data.Array.Repa                as R
import           Data.Complex
import           Data.List                      as L
import           Data.Vector.Storable           as VS
import           Data.Vector.Unboxed            as VU
import           FokkerPlanck.FourierSeries
import           FokkerPlanck.GreensFunction
import           FokkerPlanck.Histogram
import           Foreign.CUDA.Driver            as CUDA
import           FourierMethod.BlockMatrixAcc
import           FourierMethod.FourierSeries2D
import           Image.IO
import           Pinwheel.FourierSeries2D
import           STC
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Utils.Array
import           Utils.List
import           Utils.Time
import Graphics.Rendering.Chart.Easy
import Utils.Parallel
import Graphics.Rendering.Chart.Backend.Cairo
import Utils.SimpsonRule

main = do
  args@(deviceIDsStr:numPointsStr:numR2FreqStr:periodR2Str:phiFreqsStr:rhoFreqsStr:thetaFreqsStr:scaleFreqsStr:stdStr:radiusStr:sStr:idxStr:numBatchStr:aStr:bStr:deltaStr:radialFreqStr:numThreadStr:_) <-
    getArgs
  let deviceIDs = read deviceIDsStr :: [Int]
      numPoints = read numPointsStr :: Int
      numR2Freq = read numR2FreqStr :: Int
      periodR2 = read periodR2Str :: Double
      phiFreq = read phiFreqsStr :: Int
      phiFreqs = L.map fromIntegral [-phiFreq .. phiFreq]
      rhoFreq = read rhoFreqsStr :: Int
      rhoFreqs = L.map fromIntegral [-rhoFreq .. rhoFreq]
      thetaFreq = read thetaFreqsStr :: Int
      thetaFreqs = L.map fromIntegral [-thetaFreq .. thetaFreq]
      scaleFreq = read scaleFreqsStr :: Int
      scaleFreqs = L.map fromIntegral [-scaleFreq .. scaleFreq]
      std = read stdStr :: Double
      radius = read radiusStr :: Double
      s = read sStr :: Double
      idx = read idxStr :: (Int, Int)
      numBatch = read numBatchStr :: Int
      numThread = read numThreadStr :: Int
      folderPath = "output/test/HollowPinwheel"
  createDirectoryIfMissing True folderPath
  initialise []
  devs <- M.mapM device deviceIDs
  ctxs <- M.mapM (\dev -> CUDA.create dev []) devs
  ptxs <- M.mapM createTargetFromContext ctxs
  let harmonicsArray' =
        pinwheelFourierCoefficientsAnatical
          numR2Freq
          phiFreq
          rhoFreq
          thetaFreq
          scaleFreq
          (-s)
          periodR2
          (periodR2 * sqrt 2)
      hpf = idealLowPassFilter radius periodR2 numR2Freq
      hpfMat =
        A.use .
        A.fromList (A.Z A.:. (numR2Freq ^ 2) A.:. (1 :: Int)) . VS.toList $
        hpf
  plotImageRepaComplex (folderPath </> "Filter_Coef.png") .
    ImageRepa 8 .
    fromUnboxed (Z :. (1 :: Int) :. numR2Freq :. numR2Freq) . VS.convert $
    hpf
  plan <-
    makePlan
      folderPath
      emptyPlan
      numPoints
      numPoints
      numR2Freq
      (2 * thetaFreq + 1)
      (2 * scaleFreq + 1)
  harmonicsArray <- convolveFrequency plan numR2Freq hpf harmonicsArray'
  -- let harmonicsArray = centerHollow numR2Freq harmonicsArray'
  let pinwheelMat =
        A.use .
        A.fromList (A.Z A.:. (numR2Freq ^ 2) A.:. (1 :: Int)) . VS.toList $
        (harmonicsArray' IA.! idx)
      pinwheelHPMat =
        A.use .
        A.fromList (A.Z A.:. (numR2Freq ^ 2) A.:. (1 :: Int)) . VS.toList $
        (harmonicsArray IA.! idx)
  plotImageRepaComplex (folderPath </> "FilteredPinhweel_Coef.png") .
    ImageRepa 8 .
    fromUnboxed (Z :. (1 :: Int) :. numR2Freq :. numR2Freq) . VS.convert $
    (harmonicsArray IA.! idx)
  pinwheelR2 <-
    computeFourierSeriesR2StreamAcc
      ptxs
      numR2Freq
      numPoints
      1
      periodR2
      1
      numBatch
      pinwheelMat
  pinwheelHPR2 <-
    computeFourierSeriesR2StreamAcc
      ptxs
      numR2Freq
      numPoints
      1
      periodR2
      1
      numBatch
      pinwheelHPMat
  hpfR2' <-
    computeFourierSeriesR2StreamAcc
      ptxs
      numR2Freq
      numPoints
      1
      periodR2
      1
      numBatch
      hpfMat
  plotImageRepaComplex (folderPath </> "Pinwheel.png") .
    ImageRepa 8 .
    computeS . extend (Z :. (1 :: Int) :. All :. All) . sumS . rotate3D $
    pinwheelR2
  printf
    "pinwheel center: %f\n"
    (magnitude $
     (pinwheelR2 R.!
      (Z :. (0 :: Int) :. (div numPoints 2) :. (div numPoints 2 - 1))))
  plotImageRepaComplex (folderPath </> "Pinwheel_HighPass.png") .
    ImageRepa 8 .
    computeS . extend (Z :. (1 :: Int) :. All :. All) . sumS . rotate3D $
    pinwheelHPR2
  printf
    "lowpassed pinwheel center: %f\n"
    (magnitude $
     (pinwheelHPR2 R.!
      (Z :. (0 :: Int) :. (div numPoints 2) :. (div numPoints 2 - 1))))
  plotImageRepaComplex (folderPath </> "Pinwheel_Coef.png") .
    ImageRepa 8 .
    fromUnboxed (Z :. (1 :: Int) :. numR2Freq :. numR2Freq) . VS.convert $
    (harmonicsArray' IA.! idx)
  printf
    "filter center: %f\n"
    (magnitude $
     (hpfR2' R.! (Z :. (0 :: Int) :. (div numPoints 2) :. (div numPoints 2 - 1))))
  plotImageRepaComplex (folderPath </> "HighPass.png") .
    ImageRepa 8 .
    computeS . extend (Z :. (1 :: Int) :. All :. All) . sumS . rotate3D $
    hpfR2'
  plotImageRepaComplex (folderPath </> "PinwheelHollow.png") .
    ImageRepa 8 .
    computeS .
    extend (Z :. (1 :: Int) :. All :. All) .
    sumS . rotate3D . R.zipWith (-) pinwheelR2 $
    pinwheelHPR2

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
-- |
-- Module    : Statistics.Distribution.Geometric
-- Copyright : (c) 2009 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- The Geometric distribution. There are two variants of
-- distribution. First is the probability distribution of the number
-- of Bernoulli trials needed to get one success, supported on the set
-- [1,2..] ('GeometricDistribution'). Sometimes it's referred to as
-- the /shifted/ geometric distribution to distinguish from another
-- one.
--
-- Second variant is probability distribution of the number of
-- failures before first success, defined over the set [0,1..]
-- ('GeometricDistribution0').
module Statistics.Distribution.Geometric
    (
      GeometricDistribution
    , GeometricDistribution0
    -- * Constructors
    , geometric
    , geometricE
    , geometric0
    , geometric0E
    -- ** Accessors
    , gdSuccess
    , gdSuccess0
    ) where

import Control.Applicative
import Control.Monad       (liftM)
import Data.Aeson          (FromJSON(..), ToJSON, Value(..), (.:))
import Data.Binary         (Binary(..))
import Data.Data           (Data, Typeable)
import GHC.Generics        (Generic)
import Numeric.MathFunctions.Constants (m_pos_inf, m_neg_inf)
import qualified System.Random.MWC.Distributions as MWC

import qualified Statistics.Distribution as D
import Statistics.Internal

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

-- | Distribution over [1..]
newtype GeometricDistribution = GD {
      gdSuccess :: Double
    } deriving (Eq, Typeable, Data, Generic)

instance Show GeometricDistribution where
  showsPrec i (GD x) = defaultShow1 "geometric" x i
instance Read GeometricDistribution where
  readPrec = defaultReadPrecM1 "geometric" geometricE

instance ToJSON GeometricDistribution
instance FromJSON GeometricDistribution where
  parseJSON (Object v) = do
    x <- v .: "gdSuccess"
    maybe (fail $ errMsg x) return  $ geometricE x
  parseJSON _ = empty

instance Binary GeometricDistribution where
  put (GD x) = put x
  get = do
    x <- get
    maybe (fail $ errMsg x) return  $ geometricE x

instance D.Distribution GeometricDistribution where
    cumulative = cumulative

instance D.DiscreteDistr GeometricDistribution where
    probability (GD s) n
      | n < 1     = 0
      | otherwise = s * (1-s) ** (fromIntegral n - 1)
    logProbability (GD s) n
       | n < 1     = m_neg_inf
       | otherwise = log s + log (1-s) * (fromIntegral n - 1)

instance D.Mean GeometricDistribution where
    mean (GD s) = 1 / s

instance D.Variance GeometricDistribution where
    variance (GD s) = (1 - s) / (s * s)

instance D.MaybeMean GeometricDistribution where
    maybeMean = Just . D.mean

instance D.MaybeVariance GeometricDistribution where
    maybeStdDev   = Just . D.stdDev
    maybeVariance = Just . D.variance

instance D.Entropy GeometricDistribution where
  entropy (GD s)
    | s == 0 = m_pos_inf
    | s == 1 = 0
    | otherwise = negate $ (s * log s + (1-s) * log (1-s)) / s

instance D.MaybeEntropy GeometricDistribution where
  maybeEntropy = Just . D.entropy

instance D.DiscreteGen GeometricDistribution where
  genDiscreteVar (GD s) g = MWC.geometric1 s g

instance D.ContGen GeometricDistribution where
  genContVar d g = fromIntegral `liftM` D.genDiscreteVar d g

cumulative :: GeometricDistribution -> Double -> Double
cumulative (GD s) x
  | x < 1        = 0
  | isInfinite x = 1
  | isNaN      x = error "Statistics.Distribution.Geometric.cumulative: NaN input"
  | otherwise    = 1 - (1-s) ^ (floor x :: Int)

-- | Create geometric distribution.
geometric :: Double                -- ^ Success rate
          -> GeometricDistribution
geometric x = maybe (error $ errMsg x) id $ geometricE x

-- | Create geometric distribution.
geometricE :: Double                -- ^ Success rate
           -> Maybe GeometricDistribution
geometricE x
  | x >= 0 && x <= 1 = Just (GD x)
  | otherwise        = Nothing

errMsg :: Double -> String
errMsg x = "Statistics.Distribution.Geometric.geometric: probability must be in [0,1] range. Got " ++ show x

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

-- | Distribution over [0..]
newtype GeometricDistribution0 = GD0 {
      gdSuccess0 :: Double
    } deriving (Eq, Typeable, Data, Generic)

instance Show GeometricDistribution0 where
  showsPrec i (GD0 x) = defaultShow1 "geometric0" x i
instance Read GeometricDistribution0 where
  readPrec = defaultReadPrecM1 "geometric0" geometric0E

instance ToJSON GeometricDistribution0
instance FromJSON GeometricDistribution0 where
  parseJSON (Object v) = do
    x <- v .: "gdSuccess0"
    maybe (fail $ errMsg x) return  $ geometric0E x
  parseJSON _ = empty

instance Binary GeometricDistribution0 where
  put (GD0 x) = put x
  get = do
    x <- get
    maybe (fail $ errMsg x) return  $ geometric0E x

instance D.Distribution GeometricDistribution0 where
    cumulative (GD0 s) x = cumulative (GD s) (x + 1)

instance D.DiscreteDistr GeometricDistribution0 where
    probability    (GD0 s) n = D.probability    (GD s) (n + 1)
    logProbability (GD0 s) n = D.logProbability (GD s) (n + 1)

instance D.Mean GeometricDistribution0 where
    mean (GD0 s) = 1 / s - 1

instance D.Variance GeometricDistribution0 where
    variance (GD0 s) = D.variance (GD s)

instance D.MaybeMean GeometricDistribution0 where
    maybeMean = Just . D.mean

instance D.MaybeVariance GeometricDistribution0 where
    maybeStdDev   = Just . D.stdDev
    maybeVariance = Just . D.variance

instance D.Entropy GeometricDistribution0 where
  entropy (GD0 s) = D.entropy (GD s)

instance D.MaybeEntropy GeometricDistribution0 where
  maybeEntropy = Just . D.entropy

instance D.DiscreteGen GeometricDistribution0 where
  genDiscreteVar (GD0 s) g = MWC.geometric0 s g

instance D.ContGen GeometricDistribution0 where
  genContVar d g = fromIntegral `liftM` D.genDiscreteVar d g

-- | Create geometric distribution.
geometric0 :: Double                -- ^ Success rate
           -> GeometricDistribution0
geometric0 x = maybe (error $ errMsg0 x) id $ geometric0E x

-- | Create geometric distribution.
geometric0E :: Double                -- ^ Success rate
            -> Maybe GeometricDistribution0
geometric0E x
  | x >= 0 && x <= 1 = Just (GD0 x)
  | otherwise        = Nothing

errMsg0 :: Double -> String
errMsg0 x = "Statistics.Distribution.Geometric.geometric0: probability must be in [0,1] range. Got " ++ show x

max_stars_repo_name: intricate/statistics

max_issues_repo_name: intricate/statistics

max_forks_repo_name: intricate/statistics


text:
{-|
    Module      : Vocoder
    Description : Phase vocoder
    Copyright   : (c) Celina Pawlińska, 2020
                      Marek Materzok, 2021
    License     : BSD2

This module implements the phase vocoder algorithms.
The implementation is designed to be used directly or to be integrated
into some convenient abstraction (streaming or FRP).
-}
module Vocoder (
      Moduli,
      Phase,
      PhaseInc,
      Frame,
      Window,
      HopSize,
      Length,
      STFTFrame,
      FFTOutput,
      VocoderParams,
      vocoderParams,
      vocFrameLength,
      vocInputFrameLength,
      vocFreqFrameLength,
      vocHopSize,
      vocWindow,
      doFFT,
      doIFFT,
      analysisBlock,
      analysisStep,
      analysisStage,
      synthesisBlock,
      synthesisStep,
      synthesisStage,
      zeroPhase,
      volumeCoeff,
      frameFromComplex,
      frameToComplex,
      addFrames
    ) where

import Data.List
import Data.Complex
import Data.Fixed
import Data.Tuple
import Control.Arrow
import Numeric.FFT.Vector.Invertible as FFT
import Numeric.FFT.Vector.Plan as FFTp
import qualified Data.Vector.Storable as V

-- | Complex moduli of FFT frames. Represent signal amplitudes.
type Moduli = V.Vector Double

-- | Complex arguments of FFT frames. Represent signal phases.
type Phase = V.Vector Double

-- | Phase increments. Represent the deviation of the phase difference
-- between successive frames from the expected difference for the center
-- frequencies of the FFT bins.
type PhaseInc = V.Vector Double

-- | Time domain frame.
type Frame = V.Vector Double

-- | Sampled STFT window function.
type Window = Frame

-- | Offset between successive STFT frames, in samples.
type HopSize = Int

-- | Size in samples.
type Length = Int

-- | STFT processing unit.
type STFTFrame = (Moduli, PhaseInc)

-- | Frequency domain frame.
type FFTOutput = V.Vector (Complex Double)

-- | Type of FFT plans for real signals.
type FFTPlan = FFTp.Plan Double (Complex Double)

-- | Type of IFFT plans for real signals.
type IFFTPlan = FFTp.Plan (Complex Double) Double

-- | Configuration parameters for the phase vocoder algorithm.
data VocoderParams = VocoderParams{
    -- | FFT plan used in analysis stage.
    vocFFTPlan  :: FFTPlan,
    -- | FFT plan used in synthesis stage.
    vocIFFTPlan :: IFFTPlan,
    -- | STFT hop size.
    vocHopSize :: HopSize,
    -- | Window function used during analysis and synthesis.
    vocWindow :: Window
    -- TODO thread safety?
}

-- | FFT frequency frame length.
vocFreqFrameLength :: VocoderParams -> Length
vocFreqFrameLength par = planOutputSize $ vocFFTPlan par

-- | FFT frame length. Can be larger than `vocInputFrameLength` for zero-padding.
vocFrameLength :: VocoderParams -> Length
vocFrameLength par = planInputSize $ vocFFTPlan par

-- | STFT frame length.
vocInputFrameLength :: VocoderParams -> Length
vocInputFrameLength par = V.length $ vocWindow par

-- | Create a vocoder configuration.
vocoderParams :: Length -> HopSize -> Window -> VocoderParams
vocoderParams len hs wnd = VocoderParams (plan dftR2C len) (plan dftC2R len) hs wnd

-- | Apply a window function on a time domain frame.
applyWindow :: Window -> Frame -> Frame
applyWindow = V.zipWith (*)

-- | Change the vector indexing so that the sample at the middle has the number 0.
-- This is done so that the FFT of the window has zero phase, and therefore does not
-- introduce phase shifts in the signal.
rewind :: (V.Storable a) => V.Vector a -> V.Vector a
rewind vec = uncurry (V.++) $ swap $ V.splitAt (V.length vec `div` 2) vec

-- | Zero-pad the signal symmetrically from both sides.
addZeroPadding :: Length
    -> Frame
    -> Frame
addZeroPadding len v
    | diff < 0  = error $ "addZeroPadding: input is " ++ (show diff) ++ " samples longer than target length"
    | diff == 0 = v
    | otherwise = res
    where
    l = V.length v
    diff = len - l
    halfdiff = diff - (diff `div` 2)
    res = (V.++) ((V.++) (V.replicate halfdiff 0) v) (V.replicate (diff-halfdiff) 0)

-- | Perform FFT processing, which includes the actual FFT, rewinding, zero-paddding
-- and windowing.
doFFT :: VocoderParams -> Frame -> FFTOutput
doFFT par =
    FFT.execute (vocFFTPlan par) . rewind . addZeroPadding (vocFrameLength par) . applyWindow (vocWindow par)

-- | Perform analysis on a sequence of frames. This consists of FFT processing
-- and performing analysis on frequency domain frames.
analysisStage :: Traversable t => VocoderParams -> Phase -> t Frame -> (Phase, t STFTFrame)
analysisStage par ph = mapAccumL (analysisBlock par) ph

-- | Perform FFT transform and frequency-domain analysis.
analysisBlock :: VocoderParams -> Phase -> Frame -> (Phase, STFTFrame)
analysisBlock par prev_ph vec = analysisStep (vocHopSize par) (vocFrameLength par) prev_ph (doFFT par vec)

-- | Analyze a frequency domain frame. Phase from a previous frame must be supplied.
-- It returns the phase of the analyzed frame and the result.
analysisStep :: HopSize -> Length -> Phase -> FFTOutput -> (Phase, STFTFrame)
analysisStep h eN prev_ph vec =
    (ph,(mag,ph_inc))
    where
    (mag, ph) = frameFromComplex vec
    ph_inc = V.imap (calcPhaseInc eN h) $ V.zipWith (-) ph prev_ph

-- | Wraps an angle (in radians) to the range [-pi : pi].
wrap :: Double -> Double
wrap e = (e+pi) `mod'` (2*pi) - pi

calcPhaseInc :: Length -> HopSize -> Int -> Double -> Double
calcPhaseInc eN hop k ph_diff =
    (omega + wrap (ph_diff - omega)) / fromIntegral hop
    where
    omega = (2*pi*fromIntegral k*fromIntegral hop) / fromIntegral eN

-- | Perform synthesis on a sequence of frames. This consists of performing
-- synthesis and IFFT processing.
synthesisStage :: Traversable t => VocoderParams -> Phase -> t STFTFrame -> (Phase, t Frame)
synthesisStage par ph frs = mapAccumL (synthesisBlock par) ph frs

-- | Perform frequency-domain synthesis and IFFT transform.
synthesisBlock :: VocoderParams -> Phase -> STFTFrame -> (Phase, Frame)
synthesisBlock par ph fr = (id *** doIFFT par) $ synthesisStep (vocHopSize par) ph fr

-- | Synthesize a frequency domain frame. Phase from the previously synthesized frame
-- must be supplied. It returns the phase of the synthesized frame and the result.
synthesisStep :: HopSize -> Phase -> STFTFrame -> (Phase, FFTOutput)
synthesisStep hop ph (mag, ph_inc) =
    (new_ph, frameToComplex (mag, new_ph))
    where
    new_ph = V.zipWith (+) ph $ V.map (* fromIntegral hop) ph_inc

-- | Perform IFFT processing, which includes the actual IFFT, rewinding, removing padding
-- and windowing.
doIFFT :: VocoderParams -> FFTOutput -> Frame
doIFFT par =
    applyWindow (vocWindow par) . cutCenter (vocInputFrameLength par) . rewind . FFT.execute (vocIFFTPlan par)

-- | Cut the center of a time domain frame, discarding zero padding.
cutCenter :: (V.Storable a) => Length -> V.Vector a -> V.Vector a
cutCenter len vec = V.take len $ V.drop ((V.length vec - len) `div` 2) vec

-- | Zero phase for a given vocoder configuration.
-- Can be used to initialize the synthesis stage.
zeroPhase :: VocoderParams -> Phase
zeroPhase par = V.replicate (vocFreqFrameLength par) 0

-- | An amplitude change coefficient for the processing pipeline.
-- Can be used to ensure that the output has the same volume as the input.
volumeCoeff :: VocoderParams -> Double
volumeCoeff par = fromIntegral (vocHopSize par) / V.sum (V.map (**2) $ vocWindow par)

-- | Converts frame representation to complex numbers.
frameToComplex :: STFTFrame -> FFTOutput
frameToComplex = uncurry $ V.zipWith mkPolar

-- | Converts frame representation to magnitude and phase.
frameFromComplex :: FFTOutput -> STFTFrame
frameFromComplex = V.map magnitude &&& V.map phase

-- | Adds STFT frames.
addFrames :: STFTFrame -> STFTFrame -> STFTFrame
addFrames f1 f2 = frameFromComplex $ V.zipWith (+) (frameToComplex f1) (frameToComplex f2)

max_stars_repo_name: tilk/vocoder

max_issues_repo_name: tilk/vocoder

max_forks_repo_name: tilk/vocoder


text:
--Mihai Ciobanu
--Student number : s1519734
import Graphics.UI.GLUT
import Graphics.Rendering.OpenGL
import Control.Monad
import Data.Complex -- we are using complex numbers
--some functions included in Data.Complex: x :+ y = x +yi, magnitude z = sqrt ((Re z)^2 + (Im z)^2)).

 ----------------------Settings-------------------------------------------------------------------------------------------------------------------------------------------------------------
type Saturation = GLfloat
shsv = 0.8 :: Saturation -- controls the saturation of the colors used + is included in the interval [0,1]

type ValueH = GLfloat
valueH = 1 :: ValueH -- value that influences the outcome of the coloring + is included in the interval [0,1]

type Chroma = GLfloat
chroma = shsv * valueH -- the chroma is calculated as the product of the saturation and the value

diff = valueH - chroma -- it's used to smooth out the colors

type Smoothcolor = GLfloat -- HSV value that will be converted to RGB

maxiter = 1000 :: Int -- maximum number of iterations - the bigger it is, the clearer the image obtained, but the program will render the image a lot slower.
constant = 0.32 :+ 0.42 :: Complex GLfloat --constant that will directly influence the formula used for generating the image - we will obtain different fractals for different constants
bound = 20 :: GLfloat -- the constant that decides when a value "escapes".
nx = 1.5 :: GLfloat -- controls the zoom of the image as explained before isJulia

wSize = Size 700 700 -- controls the size of the window
width = 1000 -- controls the width of the image
height = 1000 -- controls the height of the image
wxh = (fromIntegral width) * (fromIntegral height) -- constant we'll be needing in our calculations.

--We need to adjust the width and height of certain fractals in order to center the image.
valuewl = 0 :: Int
valuewr = 0 :: Int
valuehu = 0 :: Int
valuehd = 0 :: Int

red1 = 0.3 :: GLfloat
green1 = 1.3 :: GLfloat
blue1 = 1 :: GLfloat
red2 = 1 :: GLfloat
green2 = 1 :: GLfloat
blue2 = 1 :: GLfloat
red3 = 1 :: GLfloat
green3 = 1 :: GLfloat
blue3 = 1 :: GLfloat
red4 = 1 :: GLfloat
green4 = 1 :: GLfloat
blue4 = 1 :: GLfloat
red5 = 1 :: GLfloat
green5 = 1 :: GLfloat
blue5 = 1 :: GLfloat
red6 = 1 :: GLfloat
green6 = 1 :: GLfloat
blue6 = 1 :: GLfloat -- these variables allow us a greater freedom in coloring our fractal - must be bigger than 0, preferrably in the interval [0..1]
--Selects the function used to do the iterations.
function :: Complex GLfloat -> Complex GLfloat
function = f1

--Different functions with their settings for different fractals (specified in brackets)
f1 :: Complex GLfloat -> Complex GLfloat
f1 z = z ^ 2
-- Method of coloring with exp (nx = 2) : -0.835 + 0.2321i (constant1), (-0.8) :+ 0.156 (constant2) , 0.285 :+ 0.01 (constant3)
-- Method of coloring with log (nx = 1) : -0.79 :+ 0.15 (constant 1), .28 + .008 (constant2), (-0.163) :+ 1.04 (constant3), (-0.12) :+ (-0.77) (constant 4) (-0.75) :+ (0.1) (constant5 - r1 0.3 g1 1.3 b1 1), 0.21 :+ 0.58 (constant6),
-- 0.15 :+ 0.58 (constant 7) 0.28 :+ 0.49 (constant 8), 0.32 :+ 0.42 (constant 9), (-0.57) :+ 0.48 (constant 10), (-0.86) :+ 0.21 (constant 11)

--Failed attempts (with examples) :
f2 :: Complex GLfloat -> Complex GLfloat
f2 z = z ^ 5 -- constant = 0.544 :+ 0 bound = 8 nx = 2 shsv = 0.8 valueH = 1 maxiter = 1000 width = 1000 height = 1000 value__ = 0 red 1 green 0.8 blue 0.5

f3 :: Complex GLfloat -> Complex GLfloat
f3 z = sqrt (sinh (z ^ 2)) -- constant 0.544 :+ 0 bound = 5 nx = 2 shsv = 0.8 valueH = 1 maxiter = 1000 height = 1000 value__ = 0 red 1 green 1 blue 1
--constant 0.065,0.122i
f4 :: Complex GLfloat -> Complex GLfloat
f4 z = (z ^ 2 / 2) + 0.26 * (z ^ 4)

--OBS : keep constant 0 :+ 0 and bound 20 and maxiter 50 for the next ones.
g5 :: Complex GLfloat -> Complex GLfloat
g5 z = z ^ 5 - 1

g5' :: Complex GLfloat -> Complex GLfloat
g5' z = 5 * (z ^ 4)

f5 :: Complex GLfloat -> Complex GLfloat
f5 z | magnitude z > 0.0001 = z - (g5 z) / (g5' z)
     | otherwise = (0 :+ 0)

g6 :: Complex GLfloat -> Complex GLfloat
g6 z = sin z

g6' :: Complex GLfloat -> Complex GLfloat
g6' z = cos z
--f6 creates two vertical "lines"
f6 :: Complex GLfloat -> Complex GLfloat
f6 z | magnitude (cos z) > 0.0001 = z - ((g6 z) / (g6' z))
     | otherwise = f1 z

f7 :: Complex GLfloat -> Complex GLfloat
f7 z | magnitude (z) > 0.00001 = (1 + 3 * z ^ 4) / (4 * z ^ 3)
     | otherwise = f1 z

f8 :: Complex GLfloat -> Complex GLfloat
f8 z | magnitude (cos z) > 0.005 = z - (sin z) / (cos z)
     | otherwise = f1 z

c1f9 = 0 :+ (-1) :: Complex GLfloat
c2f9 = 0.23 :+ 0.678 :: Complex GLfloat

f9 :: Complex GLfloat -> Complex GLfloat
f9 z = 1 + c1f9 * z ^ 5
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

----------------------OpenGL stuff----------------------------------------------------------------------------------------------------------------------------------------------------------
main :: IO ()
main = do
  (_progName, _args) <- getArgsAndInitialize
  _window <- createWindow "(Not Responding) (Not Responding) (Not Responding) Julia Set (Not Responding) (Not Responding)" --just a way to vent my frustration, ignore it
  keyboardMouseCallback $= Just keyboardMouse
  initialDisplayMode $= [DoubleBuffered]
  windowSize $= wSize
  reshapeCallback $= Just reshape
  displayCallback $= display
  actionOnWindowClose $= MainLoopReturns
  mainLoop

keyboardMouse :: KeyboardMouseCallback
keyboardMouse _key _state _modifiers _position = return ()

display = do
   clear [ColorBuffer]
   loadIdentity
   scale 0.001 0.001 (0.2::GLfloat) -- rescales the picture
   forM_ points $ \(x,y,z,t) -> -- plots the points and their color
    preservingMatrix $ do
      juliaa points
      swapBuffers
   flush

juliaa :: [(GLfloat,GLfloat,GLfloat,Color3 GLfloat)] -> IO ()
juliaa asdf = renderPrimitive Points $ do
            mapM_ drawing asdf
                where
                  drawing (x,y,z,c) = do
                  color c
                  vertex $ Vertex3 x y z

reshape :: ReshapeCallback
reshape size = do
  viewport $= (Position 0 0, size)
  postRedisplay Nothing

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

---------Functions used to generate the fractal-based image---------------------------------------------------------------------------------------------------------------------------------

-- This function decides if a point is going to be colored.
funcRec :: Complex GLfloat -> Int -> Smoothcolor -> (Bool,Smoothcolor)
funcRec z 0 sm = (False, 0)
funcRec z n sm = if (magnitude fz < bound) then funcRec (fz) (n-1) (sm) else (True, (fromIntegral (maxiter - n) - logBase 2 ((log (magnitude fz)) / (log bound))) / fromIntegral (maxiter)) where fz = function z + constant
--A smooth coloring method (will work for the majority of our functions)
-- funcRec z n sm = if ((magnitude z) < bound) then funcRec ((function z) + constant) (n-1) (sm + exp(-magnitude z)) else (True, (sm / fromIntegral (maxiter - n)))

--Calculates the mod 2 of a GLfloat - the function is needed for converting HSV to RGB in coloring.
modulo2 :: Smoothcolor -> Smoothcolor
modulo2 x | x < 2 = x
          | x > 2 = modulo2 (x-2)

--A small transformation to further smooth our coloring.
intermvalueX :: Smoothcolor -> GLfloat
intermvalueX sm = chroma * (1 - abs ((modulo2 sm) - 1))

--Converts HSV to RGB.
--We can have fun by making one of those lines black in order to see exactly which points are affected by the coloring.
coloring :: Smoothcolor -> Color3 GLfloat
coloring u | u == 0 = Color3 0 0 0
           | u >= 0 && u < 1 = Color3 ((valueH) * red1) (((intermvalueX u) + diff) * green1) ((diff) * blue1)
           | u >= 1 && u < 2 = Color3 (((intermvalueX u) + diff) * red2) ((valueH) * green2) ((diff) * blue2)
           | u >= 2 && u < 3 = Color3 ((diff) * red3) ((valueH) * green3) (((intermvalueX u) + diff) * blue3)
           | u >= 3 && u < 4 = Color3 ((diff) * red4) (((intermvalueX u) + diff) * green4) ((valueH) * blue4)
           | u >= 4 && u < 5 = Color3 (((intermvalueX u) + diff) * red5) ((diff) * green5) ((valueH) * blue5)
           | u >= 5 && u < 6 = Color3 ((valueH) * red6) ((diff) * green6) (((intermvalueX u) + diff) * blue6)

-- nx - bigger number => zoom out the picture, smaller => zoom in (I've done this before I discovered I could scale the picture in the OpenGL panel). A more correct name of the function would be isNotJulia, but we'll stick to this notation.
isJulia :: Complex GLfloat -> (Bool,Smoothcolor)
isJulia (x :+ y) = funcRec ((nx * x / (fromIntegral width)) :+ (nx * y / (fromIntegral height))) maxiter 0

-- Here we evaluate each point in the rectangle (-width, -height) (width, height) and color it if it's not in the Julia Set.
juliaP :: [(GLfloat,GLfloat,Smoothcolor)]
juliaP = [(x,y,snd j) | x <- [fromIntegral (-width + valuewl)..fromIntegral (width + valuewr)], y <- [fromIntegral (-height + valuehd)..fromIntegral (height + valuehu)], let j = isJulia (x :+ y), fst j]

--The points we have to plot and their color - "black" if the point is in the Julia Set, colored if it's not (if it escapes after a number of iterations)
points :: [(GLfloat,GLfloat,GLfloat,Color3 GLfloat)]
points =  [(x,y,0,(coloring (sm * 6))) | (x,y,sm) <- juliaP]
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

----------------------------References------------------------------------------------------------------------------------------------------------------------------------------------------
--https://en.wikipedia.org/wiki/HSL_and_HSV                                                   - for information on how HSV works and how to convert it to RGB;
--http://stackoverflow.com/questions/369438/smooth-spectrum-for-mandelbrot-set-rendering      - one smooth coloring algorithm for Julia Sets;
--http://hackage.haskell.org/package/GLUT-2.7.0.3/docs/Graphics-UI-GLUT-Callbacks-Window.html - information about OpenGL.
--https://wiki.haskell.org/OpenGLTutorial1#Lines                                              -
--https://wiki.haskell.org/OpenGLTutorial2                                                    - tutorial on OpenGL.
--http://www.juliasets.dk/Pictures_of_Julia_and_Mandelbrot_sets.htm                           - useful information about Julia Sets, Mandelbrot fractals and Newton fractals.
--http://www.karlsims.com/julia.html                                                          - interesting Julia Fractals

max_stars_repo_path: haskell_fractals/Mihai Ciobanu/Base.hs

max_stars_repo_name: MihaiAC/Other-Projects

max_issues_repo_path: haskell_fractals/Mihai Ciobanu/Base.hs

max_issues_repo_name: MihaiAC/Other-Projects

max_forks_repo_path: haskell_fractals/Mihai Ciobanu/Base.hs

max_forks_repo_name: MihaiAC/Other-Projects


text:
{-# LANGUAGE TypeFamilies #-}

import qualified Data.Foldable as F
import Test.HUnit
import Test.Framework
import Test.Framework.Providers.HUnit

import qualified Linear.V2 as V
import Numeric.AD
import Numeric.LinearAlgebra.HMatrix

import Data.SQP

import Debug.Trace

startState :: State2D
startState = vector [(-1), (-1)]

goalState :: State2D
goalState = vector [1, 1]

internalPathLength :: Int
internalPathLength = 10

costMatrix :: Int -> Matrix Double
costMatrix numInternalStates
  | numInternalStates < 1 =
      error "costMatrix can't be generated with fewer than 1 internal states!"
  | otherwise =
      let mainDiag = replicate (2 * numInternalStates) 2
          aboveDiag =
            map (\ix -> ((ix, ix + 2), (-1))) $
              take (2 * (numInternalStates - 1)) [0..]
          belowDiag =
            map (\ix -> ((ix + 2, ix), (-1))) $
              take (2 * (numInternalStates - 1)) [0..]
      in  diagBlock [diagl mainDiag +
            assoc (2 * numInternalStates, 2 * numInternalStates) 0
              (aboveDiag ++ belowDiag), matrix 1 [0]]

cost :: Vector Double -> Double
cost pathV =
  let (internalStates, _) = vectorToPath pathV
      states = startState : internalStates ++ [goalState]
      squaredDists = zipWith (\x x' -> let dx = x - x' in dot dx dx)
        states $ tail states
  in  0.5 * sum squaredDists

approxCost :: Vector Double -> (Matrix Double, Vector Double, Double)
approxCost _ = ( costMatrix internalPathLength
               , vjoin [(vjoin [negate startState, konst 0.0 (2 * (internalPathLength - 1))] +
                 vjoin [konst 0.0 (2 * (internalPathLength - 1)), negate goalState]), vector [0]]
               , 0.5 * (dot startState startState + dot goalState goalState) )

costTest :: Assertion
costTest =
  let pathV = fromList [(-10), 10]
      (costMat, costVec, costAffine) = approxCost pathV
      testApproxCost = 0.5 * (pathV `dot` (costMat #> pathV)) +
        dot pathV costVec + costAffine
  in  (testApproxCost, cost pathV) @?= (400, 400)

type State2D = Vector Double

pathToVector :: [State2D] -> Double -> Vector Double
pathToVector states stepSize = vjoin $ states ++ [vector [stepSize]]

vectorToPath :: Vector Double -> ([State2D], Double)
vectorToPath pathV = ( takesV (replicate internalPathLength 2) pathV
                     , pathV ! (2 * internalPathLength) )

obstacleRadius :: Double
obstacleRadius = 0.5

obstaclePosition :: Vector Double
obstaclePosition = vector [0.001, 0.001]

signedDistance :: State2D -> (Double, Vector Double)
signedDistance p =
  let dp = p - obstaclePosition
      dpNorm = norm_2 dp
      signedDist = dpNorm - obstacleRadius
  in  (signedDist, dp / scalar dpNorm)

ineqs :: Vector Double -> Vector Double
ineqs pathV =
  let (states, stepSize) = vectorToPath pathV
      collisionIneqs = map (negate . fst . signedDistance) states
      boundIneqs = negate stepSize
  in  fromList $ collisionIneqs ++ [boundIneqs]

approxIneqs :: Vector Double -> (Matrix Double, Vector Double)
approxIneqs pathV =
  let (states, _) = vectorToPath pathV
      (dists, normals) = unzip $ map signedDistance states
      assocList =
        concatMap (\(ix, n) -> [ ((ix, 2 * ix), negate $ n ! 0)
                               , ((ix, 2 * ix + 1), negate $ n ! 1)]) $
          zip [0..] normals
      ineqMat = diagBlock
        [ assoc (internalPathLength, 2 * internalPathLength) 0 assocList
        , matrix 1 [(-1)] ]
      affine =
        fromList $ (zipWith (+) (map negate dists) $ zipWith dot states normals)
          ++ [0]
  in  (ineqMat, affine)

unconcat :: [Int] -> [a] -> [[a]]
unconcat [] [] = []
unconcat _ [] = []
unconcat (n : ns) xs = let (xh, xt) = splitAt n xs
                       in  xh : unconcat ns xt

eqsL :: (Mode a, Scalar a ~ Double) => [a] -> [a]
eqsL pathL =
  let stepSize = last pathL
      internalStatesFlat = take (2 * internalPathLength) pathL
      internalStates = unconcat (repeat 2) internalStatesFlat
      states = (map auto $ toList startState) : internalStates ++ [map auto $ toList goalState]
      dist2 s s' = sum $ map (^ (2 :: Int)) $ zipWith (-) s s'
      squaredDists = zipWith dist2 states $ tail states
      stepSizeSqrd = stepSize * stepSize
  in  map (\d2 -> d2 - stepSizeSqrd) squaredDists

eqs :: Vector Double -> Vector Double
eqs pathV =
  vector $ eqsL $ toList pathV

approxEqs :: Vector Double -> (Matrix Double, Vector Double)
approxEqs pathV =
  let f = fromList $ eqsL $ toList pathV
      f' = fromLists $ jacobian eqsL $ toList pathV
  in  ( f', f - (f' #> pathV) )

planningTest :: Assertion
planningTest =
  let problem = Problem
                { _cost = cost
                , _approxCost = approxCost
                , _trueIneqs = ineqs
                , _approxAffineIneqs = approxIneqs
                , _eqs = eqs
                , _approxEqs = approxEqs
                , _numVariables = 2 * internalPathLength + 1
                , _numIneqs = internalPathLength + 1
                , _numEqs = internalPathLength + 1
                }
      initPathStates =
        map (\t -> startState + scalar (fromIntegral t / fromIntegral (internalPathLength + 1)) *
              (goalState - startState))
          [1..internalPathLength]
      initPathStatesPerturbed = map (+ vector [0.0, 0.0001]) initPathStates
      shortestPathLength = norm_2 $ goalState - startState
      initPathVec = pathToVector initPathStatesPerturbed $
       shortestPathLength / fromIntegral (internalPathLength + 1)
      (xResult, fResult) = optimize problem initPathVec
      (xTrue, fTrue) = ((-50), (-50)) -- whatever
  in  (xResult, fResult) @?= (xTrue, fTrue)

main = defaultMain tests
  where tests = [ -- testCase "cost-test" costTest
                -- , testCase "ineq-test" ineqsTest
                 testCase "planning-test" planningTest ]

max_stars_repo_name: giogadi/hs-trajopt

max_issues_repo_name: giogadi/hs-trajopt

max_forks_repo_name: giogadi/hs-trajopt


text:
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}

{-# LANGUAGE CPP                   #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DeriveAnyClass        #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE UndecidableInstances  #-}
{-# LANGUAGE OverloadedStrings     #-}
{-|
Module      : Grenade.Layers.Reshape
Description : Multipurpose reshaping layer
Copyright   : (c) Huw Campbell, 2016-2017
License     : BSD2
Stability   : experimental
-}
module Grenade.Layers.Reshape
  ( Reshape(..)
  ) where

import           Control.DeepSeq                (NFData (..))
import           Data.Serialize
import           GHC.Generics                   (Generic)
import           GHC.TypeLits
import           Numeric.LinearAlgebra.Data     as LA (flatten)
import           Numeric.LinearAlgebra.Static

import           Grenade.Core

import           Grenade.Onnx

-- | Reshape Layer
--
-- The Reshape layer can flatten any 2D or 3D image to 1D vector with the
-- same number of activations, as well as cast up from 1D to a 2D or 3D
-- shape.
--
-- Can also be used to turn a 3D image with only one channel into a 2D image
-- or vice versa.
data Reshape = Reshape
  deriving (Show,Generic,NFData)

instance UpdateLayer Reshape where
  type Gradient Reshape = ()
  runUpdate _ _ _ = Reshape
  reduceGradient _ = ()

instance RandomLayer Reshape where
  createRandomWith _ _ = return Reshape

instance (KnownNat a, KnownNat x, KnownNat y, a ~ (x * y)) => Layer Reshape ('D2 x y) ('D1 a) where
  type Tape Reshape ('D2 x y) ('D1 a) = ()
  runForwards _ (S2D y)   =  ((), fromJust' . fromStorable . flatten . extract $ y)
  runBackwards _ _ (S1D y) = ((), fromJust' . fromStorable . extract $ y)

instance (KnownNat a, KnownNat x, KnownNat y, KnownNat (x * z), KnownNat z, a ~ (x * y * z)) => Layer Reshape ('D3 x y z) ('D1 a) where
  type Tape Reshape ('D3 x y z) ('D1 a) = ()
  runForwards _ (S3D y)     = ((), fromJust' . fromStorable . flatten . extract $ y)
  runBackwards _ _ (S1D y) = ((), fromJust' . fromStorable . extract $ y)

instance (KnownNat y, KnownNat x, KnownNat z, z ~ 1) => Layer Reshape ('D3 x y z) ('D2 x y) where
  type Tape Reshape ('D3 x y z) ('D2 x y) = ()
  runForwards _ (S3D y)    = ((), S2D y)
  runBackwards _ _ (S2D y) = ((), S3D y)

instance (KnownNat y, KnownNat x, KnownNat z, z ~ 1) => Layer Reshape ('D2 x y) ('D3 x y z) where
  type Tape Reshape ('D2 x y) ('D3 x y z) = ()
  runForwards _ (S2D y)    = ((), S3D y)
  runBackwards _ _ (S3D y) = ((), S2D y)

instance (KnownNat a, KnownNat x, KnownNat y, a ~ (x * y)) => Layer Reshape ('D1 a) ('D2 x y) where
  type Tape Reshape ('D1 a) ('D2 x y) = ()
  runForwards _ (S1D y)   =  ((), fromJust' . fromStorable . extract $ y)
  runBackwards _ _ (S2D y) = ((), fromJust' . fromStorable . flatten . extract $ y)

instance (KnownNat a, KnownNat x, KnownNat y, KnownNat (x * z), KnownNat z, a ~ (x * y * z)) => Layer Reshape ('D1 a) ('D3 x y z) where
  type Tape Reshape ('D1 a) ('D3 x y z) = ()
  runForwards _ (S1D y)     = ((), fromJust' . fromStorable . extract $ y)
  runBackwards _ _ (S3D y)  = ((), fromJust' . fromStorable . flatten . extract $ y)

instance (KnownNat a, KnownNat b, KnownNat c, KnownNat w, KnownNat x, KnownNat y, KnownNat z, (a * b * c) ~ (w * x * y * z)) => Layer Reshape ('D3 a b c) ('D4 w x y z) where
  type Tape Reshape ('D3 a b c) ('D4 w x y z) = ()
  runForwards _ (S3D y)     = ((), fromJust' . fromStorable . flatten . extract $ y)
  runBackwards _ _ (S4D y)  = ((), fromJust' . fromStorable . flatten . extract $ y)

instance (KnownNat a, KnownNat b, KnownNat c, KnownNat d, KnownNat x, KnownNat y, KnownNat z, (a * b * c * d) ~ (x * y * z)) => Layer Reshape ('D4 a b c d) ('D3 x y z) where
  type Tape Reshape ('D4 a b c d) ('D3 x y z) = ()
  runForwards _ (S4D y)     = ((), fromJust' . fromStorable . flatten . extract $ y)
  runBackwards _ _ (S3D y)  = ((), fromJust' . fromStorable . flatten . extract $ y)

instance Serialize Reshape where
  put _ = return ()
  get = return Reshape

fromJust' :: Maybe x -> x
fromJust' (Just x) = x
fromJust' Nothing  = error "Reshape error: data shape couldn't be converted."

instance OnnxOperator Reshape where
  onnxOpTypeNames _ = ["Flatten", "Reshape"]

instance OnnxLoadableActivation Reshape where
  activationLayer = Reshape

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
module LispType
    (
    LispVal(..)
  , Array
    ) where
import Data.Ratio
import Data.Complex
import Data.Array

{-
台数的データ型の一例
LispVal型の変数が持つことのできる値の集合を定めている
選択肢のそれぞれ(コンストラクタ、と呼ばれ、| で区切られる)は、
、コンストラクタのタグとそのコンストラクタが持つことのできるデータの型を含む
この例では、LispValは次のどれか
  1. Atom - そのアトムの示す文字列を格納します。
  2. List - 他のLispValのリストを保持します(Haskellのリストは角括弧で表されます)。properリストとも呼ばれます。
  3. DottedList - Schemeの(a b . c)を表し、improperリストとも呼ばれます。これは最後以外全ての要素のリストを持ち、最後の要素を別に格納します。
  4. Number - Haskellの整数を保持します。
  5. String - Haskellの文字列を保持します。
  6. Bool - Haskellの真偽値を保持します。
-}
data LispVal =  Atom        String
              | String      String
              | Character   Char
              | Bool        Bool
              | Number      Integer
              | Float       Double
              | Complex    (Complex Double)
              | Ratio       Rational
              | Vector     (Array Int LispVal)
              | List       [LispVal]
              | DottedList [LispVal] LispVal

max_stars_repo_name: bokuo-okubo/scheme-of-bko

max_issues_repo_name: bokuo-okubo/scheme-of-bko

max_forks_repo_name: bokuo-okubo/scheme-of-bko


text:
{-# LANGUAGE RankNTypes #-}

module Observable.MCMC.AffineInvariantEnsemble where

import Control.Monad.Primitive
import Control.Monad.State.Strict
import Data.Vector.Unboxed (Vector)
import qualified Data.Vector.Unboxed as V
import Observable.Core
import Statistics.Distribution
import Statistics.Distribution.Normal

type Ensemble = Vector [Double]

-- | Generate a random value from a distribution having the property that
--   g(1/z) = z g(z).
symmetricVariate :: PrimMonad m => Observable m Double
symmetricVariate = do
  z <- unit
  return $ 0.5 * (z + 1) ^ 2

perturbParticle w0 w1 = do
  zs <- replicateM (length w0) symmetricVariate

V.zipWith (+) (V.map (* z) w0) (V.map (* (1 - z)) w1)

-- acceptanceRatio
--   :: [Double] -> [Double]
--   -> Double -> Double
--   -> ([Double] -> Double)
--   -> ([Double], Int)
-- acceptanceRatio target w0 w1 z zc =
--   let val      = target proposal - target w0 + (fromIntegral (length w0) - 1) * log z
--       proposal = zipWith (+) (map (*z) w0) (map (*(1-z)) w1)
--   in  if zc <= min 1 (exp val) then (proposal, 1) else (w0, 0)
--
-- -- | Calculate the acceptance ratio for a proposed move.
-- acceptRejectRatio
--   :: Target Double -> Double -> Vector Double -> Vector Double -> Double
-- acceptRejectRatio target e current proposed = exp . min 0 $
--     logObjective target proposed + log (isoGauss current proposed e)
--   - logObjective target current  - log (isoGauss proposed current e)

max_stars_repo_name: jtobin/deprecated-observable

max_issues_repo_name: jtobin/deprecated-observable

max_forks_repo_name: jtobin/deprecated-observable


text:
{-# LANGUAGE CPP #-}
module Control.Instances.Package.Haskell2010 () where
#ifdef MIN_VERSION_haskell2010
#if MIN_VERSION_haskell2010(1,1,2)
import Prelude ()
import Control.Monad ()
import Data.Array ()
import Data.Bits ()
import Data.Char ()
import Data.Complex ()
import Data.Int ()
import Data.Ix ()
import Data.List ()
import Data.Maybe ()
import Data.Ratio ()
import Data.Word ()
import Foreign ()
import Foreign.C ()
import Foreign.C.Error ()
import Foreign.C.String ()
import Foreign.C.Types ()
import Foreign.ForeignPtr ()
import Foreign.Marshal ()
import Foreign.Marshal.Alloc ()
import Foreign.Marshal.Array ()
import Foreign.Marshal.Error ()
import Foreign.Marshal.Utils ()
import Foreign.Ptr ()
import Foreign.StablePtr ()
import Foreign.Storable ()
import Numeric ()
import System.Environment ()
import System.Exit ()
import System.IO ()
import System.IO.Error ()
#elif MIN_VERSION_haskell2010(1,1,1)
import Data.Array ()
import Data.Char ()
import Data.Complex ()
import System.IO ()
import System.IO.Error ()
import Data.Ix ()
import Data.List ()
import Data.Maybe ()
import Control.Monad ()
import Data.Ratio ()
import System.Environment ()
import System.Exit ()
import Numeric ()
import Prelude ()
import Data.Int ()
import Data.Word ()
import Data.Bits ()
import Foreign ()
import Foreign.Ptr ()
import Foreign.ForeignPtr ()
import Foreign.StablePtr ()
import Foreign.Storable ()
import Foreign.C ()
import Foreign.C.Error ()
import Foreign.C.String ()
import Foreign.C.Types ()
import Foreign.Marshal ()
import Foreign.Marshal.Alloc ()
import Foreign.Marshal.Array ()
import Foreign.Marshal.Error ()
import Foreign.Marshal.Utils ()
#elif MIN_VERSION_haskell2010(1,1,0)
import Data.Array ()
import Data.Char ()
import Data.Complex ()
import System.IO ()
import System.IO.Error ()
import Data.Ix ()
import Data.List ()
import Data.Maybe ()
import Control.Monad ()
import Data.Ratio ()
import System.Environment ()
import System.Exit ()
import Numeric ()
import Prelude ()
import Data.Int ()
import Data.Word ()
import Data.Bits ()
import Foreign ()
import Foreign.Ptr ()
import Foreign.ForeignPtr ()
import Foreign.StablePtr ()
import Foreign.Storable ()
import Foreign.C ()
import Foreign.C.Error ()
import Foreign.C.String ()
import Foreign.C.Types ()
import Foreign.Marshal ()
import Foreign.Marshal.Alloc ()
import Foreign.Marshal.Array ()
import Foreign.Marshal.Error ()
import Foreign.Marshal.Utils ()
#elif MIN_VERSION_haskell2010(1,1,0)
import Data.Array ()
import Data.Char ()
import Data.Complex ()
import System.IO ()
import System.IO.Error ()
import Data.Ix ()
import Data.List ()
import Data.Maybe ()
import Control.Monad ()
import Data.Ratio ()
import System.Environment ()
import System.Exit ()
import Numeric ()
import Prelude ()
import Data.Int ()
import Data.Word ()
import Data.Bits ()
import Foreign ()
import Foreign.Ptr ()
import Foreign.ForeignPtr ()
import Foreign.StablePtr ()
import Foreign.Storable ()
import Foreign.C ()
import Foreign.C.Error ()
import Foreign.C.String ()
import Foreign.C.Types ()
import Foreign.Marshal ()
import Foreign.Marshal.Alloc ()
import Foreign.Marshal.Array ()
import Foreign.Marshal.Error ()
import Foreign.Marshal.Utils ()
#elif MIN_VERSION_haskell2010(1,0,0)
import Data.Array ()
import Data.Char ()
import Data.Complex ()
import System.IO ()
import System.IO.Error ()
import Data.Ix ()
import Data.List ()
import Data.Maybe ()
import Control.Monad ()
import Data.Ratio ()
import System.Environment ()
import System.Exit ()
import Numeric ()
import Prelude ()
import Data.Int ()
import Data.Word ()
import Data.Bits ()
import Foreign ()
import Foreign.Ptr ()
import Foreign.ForeignPtr ()
import Foreign.StablePtr ()
import Foreign.Storable ()
import Foreign.C ()
import Foreign.C.Error ()
import Foreign.C.String ()
import Foreign.C.Types ()
import Foreign.Marshal ()
import Foreign.Marshal.Alloc ()
import Foreign.Marshal.Array ()
import Foreign.Marshal.Error ()
import Foreign.Marshal.Utils ()
#endif
#endif

max_stars_repo_name: clintonmead/ghc-instances

max_issues_repo_name: clintonmead/ghc-instances

max_forks_repo_name: clintonmead/ghc-instances


text:
{-# LANGUAGE LambdaCase #-}
-- |
-- Module: Analysis.PerfModel
--
-- An experimental approach to obtaining strengths, which takes Elo (or
-- Elo-like) ratings as derived from a performance model. With a performance
-- model in hand, we can use https://stats.stackexchange.com/a/44142 to
-- obtain the victory probability of a 1500-rated racer given a field of
-- racers, and use that as a measure of strength.
--
-- According to the model, player performances follow a gamma distribution,
-- with the outcome 0 corresponding to an ideal performance. The shape
-- parameter of the distribution is fixed, while the rate parameter varies
-- across players. A higher rate parameter amounts to results that are both
-- more consistent and closer to the ideal. (Pinning both mean and standard
-- deviation of player performance to a single parameter is the main
-- simplifying assumption of the model.)
--
-- Note that the performance model winning probabilities for pairs of players
-- agree exactly with the conventional Elo formula if the shape parameter of
-- the gamma distribution is taken to be 1.
module Analysis.PerfModel
    ( perfModelStrength
    , perfModelTopStrength
    , example
    ) where

import Orbital (orbitalDistr, initialRating)
import qualified Util.Combinations as Util

import qualified Numeric.Integration.TanhSinh as Integration
import Statistics.Distribution

import Data.Profunctor
import qualified Control.Foldl as L
import qualified Data.IntMap.Strict as IntMap
import Data.IntMap.Strict (IntMap, (!))
import qualified Data.IntSet as IntSet
import Data.IntSet (IntSet)
import qualified Control.Foldl as L

-- | PDF for a race victory against a field of racers.
raceWinPDF
    :: Int       -- ^ Gamma shape.
    -> [Double]  -- ^ Ratings of the field of opponents.
    -> Double    -- ^ Rating of the racer.
    -> Double    -- ^ Performance.
    -> Double
raceWinPDF gsh rs r t = density (orbitalDistr gsh r) t
    * L.fold (lmap winVersus L.product) rs
    where
    winVersus r' = 1 - cumulative (orbitalDistr gsh r') t

-- | Integrating the race win PDF gives the likelihood of victory, whose
-- reciprocal we use as a strength metric.
perfModelStrength
    :: Int       -- ^ Gamma shape.
    -> [Double]  -- ^ Ratings of the field of opponents.
    -> Double
perfModelStrength gsh rs = 1 / integ
    where
    integ = Integration.result . Integration.absolute 1e-6 $
        Integration.nonNegative Integration.trap (raceWinPDF gsh rs initialRating)

-- TODO: positionPDF and positionPDFPre are still the same as in
-- Analysis.PerfModel.Reference.Precompute . Conside redefining them so that
-- topPDF can share code with them.

-- | Like 'raceWinPDF', but for an arbitrary position on the scoreboard.
positionPDF
    :: Int       -- ^ Gamma shape.
    -> Int       -- ^ Target position.
    -> [Double]  -- ^ Ratings of the field of opponents.
    -> Double    -- ^ Rating of the racer.
    -> Double    -- ^ Performance.
    -> Double
positionPDF gsh pos rs r t = density (orbitalDistr gsh r) t
    * L.fold L.sum (L.fold L.product . toFactors <$> combos)
    where
    n = length rs
    combos = Util.combinations n (pos-1) [1..n]
    loseVersus r' = cumulative (orbitalDistr gsh r') t
    -- winVersus r' = 1 - loseVersus r'
    lvs = IntMap.fromList (zip [1..n] (loseVersus <$> rs))
    wvs = (1 -) <$> lvs
    toFactors (above, below) = map (lvs !) above ++ map (wvs !) below

-- | A variant of 'positionPDF' that takes precomputed probabilities.
positionPDFPre :: Int                            -- ^ Target position.
               -> (IntMap Double, IntMap Double) -- ^ Win and loss probabilities
                                                 -- (given a probe laptime) for
                                                 -- for each racer on the
                                                 -- scoreboard.
               -> Double                         -- ^ Probability density of the
                                                 -- probe laptime.
               -> Double
positionPDFPre pos (lvs, wvs) p =
    p * L.fold L.sum (L.fold L.product . toFactors <$> combos)
    where
    n = IntMap.size lvs
    combos = Util.combinations n (pos-1) [1..n]
    toFactors (above, below) = map (lvs !) above ++ map (wvs !) below

-- | Like 'positionPDF', but for a top-N finish.
topPDF
    :: Int       -- ^ Gamma shape.
    -> Int       -- ^ Target position.
    -> [Double]  -- ^ Ratings of the field of opponents.
    -> Double    -- ^ Rating of the racer.
    -> Double    -- ^ Performance.
    -> Double
topPDF gsh pos rs r t = (probeDensity *) . L.fold L.sum $
    L.fold L.sum . partials <$> validPositions
    where
    validPositions = zipWith const [1..pos] (r : rs)
    loseVersus r' = cumulative (orbitalDistr gsh r') t
    lvs = IntMap.fromList (zip [1..] (loseVersus <$> rs))
    wvs = (1 -) <$> lvs
    probeDensity = density (orbitalDistr gsh r) t
    partials = Util.processCombsInt alg (length rs) . subtract 1
    alg = \case
        Util.LeafF Nothing -> 0
        Util.LeafF (Just rest) ->
            L.fold L.product ((wvs !) <$> IntSet.toList rest)
        Util.FlowerF a rest -> (lvs ! a)
            * L.fold L.product ((wvs !) <$> IntSet.toList rest)
        Util.BranchF a bs -> (lvs ! a) * L.fold L.sum bs

-- | Like 'perfModelStrength', but for a top-N finish. Note that the
-- computational cost grows quickly with the length of the list of ratings.
perfModelTopStrength
    :: Int       -- ^ Gamma shape.
    -> Int       -- ^ Target position.
    -> [Double]  -- ^ Ratings of the field of racers.
    -> Double
perfModelTopStrength gsh pos rs = 1 / integ
    where
    integ = Integration.result . Integration.absolute 1e-6 $
        Integration.nonNegative Integration.trap (topPDF gsh pos rs initialRating)

example :: [Double]
example = [2200,2100,1900,1870,1850,1600]
--example = [2200,2100,1900,1870,1850,1600,1590,1570,1510,1420,1370,1350,1225]

-- >$> perfModelStrength [2200,2100,1900,1870,1850]
--
-- >$> :set +s
--

-- >$> perfModelTopStrength 5 example

max_stars_repo_name: duplode/elo-zs

max_issues_repo_name: duplode/elo-zs

max_forks_repo_name: duplode/elo-zs