text:
{-# Language ConstraintKinds            #-}
{-# Language FlexibleContexts           #-}
{-# Language FlexibleInstances          #-}
{-# Language GeneralizedNewtypeDeriving #-}
{-# Language TypeFamilies               #-}

module Geometry.Base
  ( -- * Types
    -- ** Point represenations
    Cmp
  , XY
  -- ** Directed values
  , Direction (..)
  -- *** Direction isomorphisms
  , deg, asDeg
  , rad, asRad
  , turns, asTurns
  -- * Points in affine space
  , Metric (..), Pnt
  , Affine (..), roundUp, asAffine
  -- ** Predicates
  , isOrthogonal, isCollinear, isOpposite, isZero, opposite
  -- ** Vector and point operations
  , dot, det, cross, norm, distance, (<-?->)
  , angle, normalize, columns, azimuth
  -- * Linear transformations
  , TMatrix, transformOrientation
  , Trans (..)
  , transformCmp, transformXY, transformAt
  -- ** transformations
  , translate', translate, superpose, at, at'
  , rotate, rotateAt, rotateAt', reflect, reflectAt
  , along, along', on
  , scale, scaleX, scaleY, scaleAt', scaleAt
  , scaleXAt', scaleXAt, scaleYAt', scaleYAt, scaleFig
  -- * Manifolds and curves
  , Manifold (..), Bounding (..)
  , (->@), (->@?), (@->), (@->?)
  , start, end, paramL, projectL
  , Curve (..), PointLocation (..), ClosedCurve(..)
  , tangentLine
  -- * Figures
  , Figure (..), Box, pointBox
  -- ** Figures' size and bounding box corners.
  , figureHeight, figureWidth
  , corner, left, right, lower, upper
  -- * Fuzzy equality
  , AlmostEq
  -- ** Fuzzy inequalities
  , (~<=), (~>=), (~=)
  , (#)
  )
where

--import Data.Fixed (mod')
import Data.Complex
import Data.List.Extra (minimumOn)
import Data.Semigroup hiding (First)
import Data.Maybe
import Numeric.MathFunctions.Comparison

------------------------------------------------------------
-- | Type alias for a complex number.
type Cmp = Complex Double

-- | Type alias for a coordinate representation.
type XY = (Double, Double)

-- | Constraint for a point in a metric affine space
type Pnt a = (Affine a, Trans a, Metric a)

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

-- | Flipped application used for right chaining of transformations.
infixl 5 #

(#) :: a -> (a -> c) -> c
(#) = flip ($)

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

infix 4 ~=, ~<=, ~>=

-- | Type class for values for which equality could be stated with some known tolerance
class AlmostEq a where
  -- | The approximate equality operator.
  (~=) :: a -> a -> Bool

instance AlmostEq Int where a ~= b = a == b

instance AlmostEq Double where
  a ~= b = a == b
           || abs (a - b) < 1e-10 * abs(a + b)
           || abs (a - b) < 1e-10

instance AlmostEq (Complex Double) where
  (ax:+ay) ~= (bx:+by) = ax ~= bx && ay ~= by

instance (AlmostEq a, AlmostEq b) => AlmostEq (a, b) where
  (a1, b1) ~= (a2, b2) = a1 ~= a2 && b1 ~= b2

instance AlmostEq a => AlmostEq [a] where
  as ~= bs = and $ zipWith (~=) as bs

instance AlmostEq a => AlmostEq (Maybe a) where
  Just a ~= Just b = a ~= b
  Nothing ~= Nothing = True
  _ ~= _ = False

-- | The less or almost equal relation.
(~<=) :: (AlmostEq a, Ord a) => a -> a -> Bool
a ~<= b = a ~= b || a < b

-- | The greater or almost equal relation.
(~>=) :: (AlmostEq a, Ord a) => a -> a -> Bool
a ~>= b = a ~= b || a > b

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

-- | The class for points in a metric space
class Metric a where
  {-# MINIMAL dist | dist2 #-}
  -- | The distance between two points, or a norm of the difference between two vectors.
  dist :: a -> a -> Double
  dist x y = sqrt (dist2 x y)

  -- | The squared distance between two points, or a norm of the difference between two vectors.
  dist2 :: a -> a -> Double
  dist2 x y = dist x y ** 2

instance Metric Double where
  dist x y = abs (x - y)
  dist2 x y = (x - y)**2

instance Metric Cmp where
  dist x y = magnitude (x - y)

instance Metric XY where
  dist2 (ax, ay) (bx, by) = (ax-bx)**2 + (ay-by)**2

instance Metric a => Metric (Maybe a) where
  dist x y = fromMaybe 0 (dist <$> x <*> y)
  dist2 x y = fromMaybe 0 (dist2 <$> x <*> y)

------------------------------------------------------------
mod' :: (Real a) => a -> a -> a
{-# INLINE mod' #-}
mod' n d = n - (fromInteger $ floor ((toRational n) / (toRational d))) * d

mod2pi :: Double -> Double
{-# INLINE mod2pi #-}
mod2pi n = n - (fromInteger $ floor (n / (2*pi))) * 2 * pi

mod360 :: Double -> Double
{-# INLINE mod360 #-}
mod360 n = n - (fromInteger $ floor (n / 360)) * 360

{- | The representation of a directed value isomorphic either to an angle
 or to a complex number or coordinates.
-}
newtype Direction = Direction Double
  deriving (Num, Fractional, Floating, Real, Enum)

-- | Constructs a directed value from an angle given in degrees. Invers of `deg`.
asDeg :: Double -> Direction
asDeg = Direction . mod360

-- | Returns a representation of a directed value as an angle in degrees. Invers of `asDeg`.
deg :: Direction -> Double
deg (Direction a) = mod360 a

-- | Constructs a directed value from an angle given in radians. Invers of `rad`.
asRad :: Double -> Direction
asRad r = asDeg (180 * r / pi)

-- | Returns a representation of a directed value as an angle in radians. Invers of `asDeg`.
rad :: Direction -> Double
rad (Direction a) = mod2pi (a / 180 * pi)

-- | Constructs a directed value from a number of turns. Invers of `turns`.
asTurns :: Double -> Direction
asTurns a = asDeg $ a*360

-- | Returns a representation of a directed value as a number of turns. Invers of `asTurns`.
turns :: Direction -> Double
turns (Direction a) = mod360 a / 360

instance Show Direction where
  show a = show (round $ deg a) <> "°"

instance Metric Direction where
  dist x y = abs (rad x - rad y)
  dist2 x y = (rad x - rad y)**2

instance AlmostEq Direction where  a ~= b = rad a ~= rad b

instance Eq Direction  where  a == b = a ~= b

instance Ord Direction where  a <= b = a == b || rad a < rad b

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

-- | The representation of a linear transformation matrix.
newtype TMatrix = TMatrix ((Double, Double, Double), (Double, Double, Double))
  deriving Show

instance Semigroup TMatrix where
  m1 <> m2 = TMatrix ((b1*c2+a1*a2, b1*d2+a1*b2, b1*y2+a1*x2+x1),
                      (c2*d1+a2*c1, d1*d2+b2*c1, d1*y2+y1+c1*x2))
    where
      TMatrix ((a1, b1, x1), (c1, d1, y1)) = m1
      TMatrix ((a2, b2, x2), (c2, d2, y2)) = m2

instance Monoid TMatrix where
  mempty = TMatrix ((1,0,0),(0,1,0))

inv :: TMatrix -> TMatrix
inv (TMatrix ((a, b, x), (c, d, y)))
  | d == 0 = error "Transformation is not invertible!"
  | otherwise = res
  where dt = a*d-b*c
        res = TMatrix ( (d/dt, -b/dt, (b*y-d*x)/dt)
                      , (-c/dt, a/dt, (c*x-a*y)/dt))

-- | Class for objects which could be transformed linearly.
class Trans a where
  {-# MINIMAL transform #-}
  -- | The general linear transformation.
  transform :: TMatrix -> a -> a

instance {-# OVERLAPPABLE #-} (Functor f, Trans a) => Trans (f a) where
  transform  = fmap . transform

instance {-# OVERLAPPING #-} Trans Cmp where
  transform = transformCmp

instance {-# OVERLAPPING #-} Trans XY where
  transform  = transformXY

instance Trans Direction where
  transform t  = asCmp . transformCmp t . cmp

-- | Returns @1@ if transformation preserves curve orientation, and @-1@ otherwise.
transformOrientation :: TMatrix -> Double
transformOrientation (TMatrix ((a,b,_),(c,d,_))) =
  signum $ det ((a,b),(c,d))

-- | The rotation matrix.
rotateT :: Double -> TMatrix
rotateT a = TMatrix ((cos a, -(sin a), 0), (sin a, cos a, 0))

-- | The reflection matrix.
reflectT :: Double -> TMatrix
reflectT a = TMatrix ((cos (2*a), sin (2*a), 0), (sin (2*a), -(cos (2*a)), 0))

-- | The translation matrix.
translateT :: XY -> TMatrix
translateT (x, y) = TMatrix ((1, 0, x), (0, 1, y))

-- | The scaling matrix.
scaleT :: Double -> Double -> TMatrix
scaleT x y = TMatrix ((x, 0, 0), (0, y, 0))

-- | The transformation of a complex number. Used in class implementations.
transformCmp :: TMatrix -> Cmp -> Cmp
transformCmp t = cmp . transformXY t . xy

-- | The transformation of coordinates. Used in class implementations.
transformXY :: TMatrix -> XY -> XY
transformXY (TMatrix ((a11, a12, sx), (a21, a22, sy))) (x, y) =
    (a12*y + a11*x + sx, a22*y + a21*x + sy)

-- | The transformation with shifted origin.
transformAt :: (Trans a, Affine p) => p -> (a -> a) -> a -> a
transformAt p t = translate' p . t . translate' (opposite p)

-- | The translation of an object.
translate' :: (Trans a, Affine p) => p -> a -> a
translate' = transform . translateT . xy

-- | The translation leading to a superposition of tho points.
superpose :: (Trans a, Affine p1, Affine p2) => p1 -> p2 -> a -> a
superpose p1 p2 = translate' (cmp p2 - cmp p1)

-- | The isotropic scaling of an object.
scale :: Trans a => Double -> a -> a
scale s = transform (scaleT s s)

-- | The scaling of an object in x-direction.
scaleX :: Trans a => Double -> a -> a
scaleX s = transform (scaleT s 1)

-- | The scaling of an object in y-direction.
scaleY :: Trans a => Double -> a -> a
scaleY s = transform (scaleT 1 s)

-- | The isotropic scaling of an object against a given point.
scaleAt' :: (Trans a, Affine p) => p -> Double -> a -> a
scaleAt' p s = transformAt p (scale s)

-- | The isotropic scaling of an object in x-direction against a given point.
scaleXAt' :: (Trans a, Affine p) => p -> Double -> a -> a
scaleXAt' p s = transformAt p (scaleX s)

-- | The isotropic scaling of an object in y-direction against a given point.
scaleYAt' :: (Trans a, Affine p) => p -> Double -> a -> a
scaleYAt' p s = transformAt p (scaleY s)

-- | The rotation of an object against the origin.
rotate :: Trans a => Direction -> a -> a
rotate = transform . rotateT . rad

-- | The rotation of an object against a given point.
rotateAt' :: (Trans a, Affine p) => p -> Direction -> a -> a
rotateAt' p a = transformAt p (rotate a)

-- | The reflection of an object against the direction passing through the origin.
reflect :: Trans a => Direction -> a -> a
reflect d = transform $ reflectT $ rad d

-- | The reflection of an object against the affine manifold.
reflectAt :: (Trans a, Manifold m, Affine m) => m -> a -> a
reflectAt l = transformAt (start l) (reflect (angle l))

-- | Moves an object along given vector.
translate :: Trans f => XY -> (f -> f)
translate = translate'

-- | Scales  an object simmetrically (isotropically) against a given point.
scaleAt :: Trans f => XY -> Double -> (f -> f)
scaleAt = scaleAt'

-- | Scales  an object along x-axis against a given point.
scaleXAt :: Trans f => XY -> Double -> (f -> f)
scaleXAt = scaleXAt'

-- | Scales  an object along y-axis against a given point.
scaleYAt :: Trans f => XY -> Double -> (f -> f)
scaleYAt = scaleYAt'

-- | Scales  an object simmetrically (isotropically) against a given point.
scaleFig :: (Figure f, Trans f) => Double -> (f -> f)
scaleFig s f = scaleAt' (refPoint f) s f

-- | Rotates  an object  against a given point.
rotateAt :: Trans f => XY -> Direction -> (f -> f)
rotateAt = rotateAt'

-- | Generalized version of  `at` transformer.
at' :: (Affine p, Figure f) => p -> (f -> f)
at' p fig = superpose (refPoint fig) p fig

-- | Moves an object so that it's `refPoint` coinsides with a given one.
at :: Figure f => XY -> (f -> f)
at = at'

-- | Generalized version of  `along` transformer.
along' :: (Figure f, Affine v, Affine f) => v -> (f -> f)
along' v l = rotateAt' (refPoint l) (angle v - angle l) l

{- | Rotates the figure which is `Affine` instance against it's `refPoint` so that it's
 refference angle (given by `angle`) councides with a given one.

> let a = aSegment # at (1,0) # along 30 #: "a"
>     t = aTriangle # along' a
>     s = aSquare # at (2,0) # along' t
> in a <+> t <+> s

<< figs/along.svg >>
-}
along :: (Figure f, Affine f) => Double -> (f -> f)
along = along' . asDeg

-- | Locates an affine object on a given curve at
-- given parameter and aligns it along a tangent to a curve at this point.
--
-- > let c = aCircle
-- > in c <+>
-- >    aPoint # on c 0.1 <+>
-- >    aSegment # scale 0.5 # on c 0.3 <+>
-- >    aSquare # scale 0.5 # on c 0.6 <+>
-- >    aTriangle # scale 0.5 # on c 0.9
--
-- << figs/on.svg>>
--
on :: (Figure f, Affine f, Curve c) => c -> Double -> (f -> f)
on c x = along' (tangent c x) . at' (c @-> x)

------------------------------------------------------------
{- | Class representing points in 2d affine space.
 Instances may be points, coordinates, complex numbers, 2d-vectors.
-}
class Affine a where
  {-# MINIMAL (asCmp, cmp) | (asXY, xy) #-}

  -- | Constructs an affine point from a complex number.
  asCmp :: Cmp -> a
  asCmp = asXY . xy

  -- | Returns a representation of an afine point as a complex number.
  cmp :: a -> Cmp
  cmp p = let (x, y) = xy p in x :+ y

  -- | Constructs an affine point from a coordinates.
  asXY :: XY -> a
  asXY = asAffine

  -- | Returns a representation of an afine point as a coordinate.
  xy :: a -> XY
  xy p = let x :+ y = cmp p in (x, y)

  -- | Returns x-coordinate of an afine point.
  getX :: a -> Double
  getX = fst . xy

  -- | Returns y-coordinate of an afine point.
  getY :: a -> Double
  getY = snd . xy

instance Affine Cmp where
  cmp = id
  asCmp = id

instance Affine XY where
  xy = id
  asXY = id

instance Affine Direction where
  cmp = cis . rad
  asCmp = asRad . phase

instance Affine a => Affine (Maybe a) where
  cmp = maybe 0 cmp
  asCmp = pure . asCmp

{- | A general homomorphism between affine instances (via `Cmp`).

>>> asAffine (1 :: Cmp) :: XY
(1.0,0.0)

>>> asAffine ((2,3) :: XY) :: Segment
Segment (0.0 :+ 0.0, 2.0 :+ 3.0)

>>> asAffine (aSegment # rotate 30) :: Direction
30
-}
asAffine :: (Affine a, Affine b) => a -> b
asAffine = asCmp . cmp

-- | Returns `True` if two points represent orthogonal vectors.
isOrthogonal :: (Affine a, Affine b) => a -> b -> Bool
isOrthogonal a b = not (isZero a || isZero b)
                   && normalize a `dot` normalize b ~= 0

-- | Returns `True` if two points represent collinear vectors.
isCollinear :: (Affine a, Affine b) => a -> b -> Bool
isCollinear a b = not (isZero a || isZero b)
                  && normalize a `cross` normalize b ~= 0

-- | Returns `True` if two points represent collinear opposite vectors.
isOpposite :: (Affine a, Affine b) => a -> b -> Bool
isOpposite a b = cmp a + cmp b ~= 0

-- | Returns `True` if the vector is trivial or a point is equal to the origin.
isZero :: Affine a => a -> Bool
isZero a = cmp a ~= 0

-- | Returns the opposite vector to a given one.
opposite :: Affine a => a -> a
opposite = asCmp . negate . cmp

-- | Returns the deriection given by two points.
azimuth :: (Affine a, Affine b) => a -> b -> Direction
azimuth p1 p2 = asCmp (cmp p2 - cmp p1)

-- | The dot product of two vectors (points).
dot :: (Affine a, Affine b) => a -> b -> Double
dot a b = let (xa, ya) = xy a
              (xb, yb) = xy b
              in xa*xb + ya*yb

-- | Determinant of a matrix, composed of two vector rows.
det :: (Affine a, Affine b) => (a, b) -> Double
det (a, b) = let (xa, ya) = xy a
                 (xb, yb) = xy b
             in xa*yb - ya*xb

-- | Z-component of a cross product of two vectors.
cross :: (Affine a, Affine b) => a -> b -> Double
cross a b = det (a, b)

-- | The norm of a vector, or a dist from a point to the origin.
norm :: Affine a => a -> Double
norm = magnitude . cmp

-- | The distance between two affine points.
distance :: (Affine a, Metric a, Affine b, Metric b) => a -> b -> Double
distance p1 p2 = cmp p1 `dist` cmp p2

-- | Operator form for `distance` function.
(<-?->):: (Affine a, Metric a, Affine b, Metric b) => a -> b -> Double
(<-?->) = distance

-- | The normalized vector, or a projection of a point on a unit circle.
normalize :: Affine a => a -> a
normalize v
  | cmp v == 0 = v
  | otherwise = asCmp . cis . phase . cmp $ v

roundUp :: Affine a => Double -> a -> a
roundUp d = asXY . (\(x,y) -> (rounding x, rounding y)) . xy
  where rounding x = fromIntegral (ceiling (x /d) :: Integer) * d

-- | The direction of the vector or a point.
angle :: Affine a => a -> Direction
angle = asAffine

-- | The matrix, composed of two vector columns.
columns :: Affine a => (a, a) -> (a, a)
columns (a, b) = ( asXY (getX a, getX b)
                 , asXY (getY a, getY b))

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

-- | Type representing range of the internal parameter @x@ of the manifold.
data Bounding = Unbound -- ^ x ∈ (−∞, ∞)
              | Semibound -- ^ x ∈ [0, ∞)
              | Bound -- ^ x ∈ [0, 1]
              deriving Show

-- | Class representing 1-dimensional manifolds, e.g. lines, curves of angles, parameterized by real value.
class (Metric (Domain m), Affine (Domain m)) => Manifold m where
  {-# MINIMAL param, project, (isContaining | (paramMaybe, projectMaybe)) #-}

  {- | The domain in which the manifold is embedded.
Here are some instances:

> type instance Domain Angle = Direction
> type instance Domain Circle = Cmp
> type instance Domain Polyline = Cmp
> type instance Domain (Plot a) = a
> type instance Domain Line = Cmp
-}
  type Domain m

  -- | Returns a point on a manifold for a given parameter.
  param :: m -> Double -> Domain m

  -- | Returns a parameter, corresponding to a normal point projection on a manifold.
  project :: m -> Domain m -> Double

  -- | Returns `True` if point belongs to the manifold.
  isContaining :: (Metric p, Affine p) => m -> p -> Bool
  isContaining c p =
    case projectMaybe c (asAffine p) >>= paramMaybe c of
      Just p' -> p' `distance` p ~= 0
      Nothing -> False

  {- | Returns a point on a manifold for given parameter, or Nothing
     if parameter runs out of the range defined for a manifold. -}
  paramMaybe :: m -> Double -> Maybe (Domain m)
  paramMaybe m x = if m `isContaining` p then Just p else Nothing
    where p = param m x

  -- | Returns a parameter for a point on a manifold, or nothing
  -- if parameter could not be found.
  projectMaybe :: m -> Domain m -> Maybe Double
  projectMaybe m p = if inBounds x then Just x else Nothing
    where x = project m p
          inBounds y = case bounds m of
            Unbound -> True
            Semibound -> y >= 0
            Bound -> y >= 0 && y <= 1

  {- | Returns bounds for a parameter of a manifold.
    @[a, b]@ for finite, @[a]@ for semibound manifolds (@a@ is lower bound),
    @[]@ for unbound manifolds
    -}
  bounds :: m -> Bounding
  bounds _ = Bound

  -- | The internal length unit of the curve,  which maps the parameter to the length of the curve.
  unit :: m -> Double
  unit _ = 1

  {- | The dist between a manifold and a given point.

>>> aCircle # distTo (point (2,0))
1.0
>>> aCircle # distTo (point (2,0))
1.0
-}
  distanceTo :: (Metric p, Affine p) => p -> m -> Double
  distanceTo p m = p `dist` p'
    where p' = case p ->@? m of
                 Just x -> asAffine $ m @-> x
                 Nothing -> minimumOn (dist p) (asAffine . param m <$> ends)
          ends = case bounds m of
            Unbound -> []
            Semibound -> [0]
            Bound -> [0,1]

instance Manifold m => Manifold (Maybe m) where
  type Domain (Maybe m) = Domain m
  param = maybe (const (asCmp 0)) param
  project = maybe (const 0) project
  isContaining = maybe (const False) isContaining
  paramMaybe m x = m >>= (`paramMaybe` x)
  projectMaybe m x = m >>= (`projectMaybe` x)
  bounds = maybe Bound bounds
  unit = maybe 1 unit

infix 8 @->
-- | Operator form of the `param` function.
(@->) ::  Manifold m => m -> Double -> Domain m
(@->) = param

infix 8 @->?
-- | Operator form of the `paramMaybe` function.
(@->?) :: Manifold m => m -> Double -> Maybe (Domain m)
(@->?) = paramMaybe

infix 8 ->@
-- | Generalized operator form for `project` with flipped arguments.
(->@) :: (Affine a, Manifold m) => a -> m -> Double
(->@) = flip project . asAffine

infix 8 ->@?
-- | Operator for `projectMaybe` with flipped arguments.
(->@?) :: (Affine a, Manifold m) => a -> m -> Maybe Double
(->@?) = flip projectMaybe . asAffine

{- | Point on the curve, parameterized by length.

>>> aCircle # paramL 1 # rad . asCmp
1.0
>>> (aTriangle # paramL 1) == (aTriangle # vertex 1)
True
>>> (aTriangle # paramL 2) == (aTriangle # vertex 2)
True
-}
paramL :: Manifold m => Double -> m -> Domain m
paramL l m = m @-> (l / unit m)

{- | Projection on a manifold, parameterized by length.

>>> aCircle # projectL (asRad 1)
1.0
>>> aTriangle # projectL (aTriangle # vertex 1)
1.0
>>> aTriangle # projectL (aTriangle # vertex 2)
2.0
-}
projectL :: (Affine p, Manifold m) => p -> m -> Double
projectL p m = unit m * (p ->@ m)

{- | The starting point on the manifold.

>>> start aCircle
1.0:+0.0
-}
start :: Manifold m => m -> Domain m
start m = param m 0

-- | Returns the end point of the segment.
end :: Manifold m => m -> Domain m
end m = param m 1

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

-- | Class representing a curve as a  continuous locally smooth manifold in affine space.
class (Figure c, Trans c, Manifold c) => Curve c  where
  {-# MINIMAL normal | tangent #-}

  -- | The tangent direction for a given parameter on the curve.
  tangent :: c -> Double -> Direction
  tangent c t = normal c t - 90

  -- | The normal direction for a given parameter on the curve.
  normal :: c -> Double -> Direction
  normal c t = tangent c t + 90

tangentLine :: (Curve c, Figure l, Affine l) => c -> Double -> l
tangentLine c x = asCmp 1 # along' (tangent c x) # at' (c @-> x)

instance Curve c => Curve (Maybe c) where
  tangent = maybe (const 0) tangent
  normal = maybe (const 0) normal

-- |  Class representing a closed region
class Curve c => ClosedCurve c where
  {-# MINIMAL location | isEnclosing #-}

  {- | Returns the location of a point with respect to the closed curve area.
>>> location aCircle <$> ([(0,0),(1,0),(2,0)] :: [XY])
[Inside, OnCurve, Outside]
    -}
  location :: (Metric p, Affine p) => c -> p -> PointLocation
  location c p | isContaining c p = OnCurve
               | isEnclosing c p = Inside
               | otherwise = Outside

  {- | Returns `True` if point belongs to the region
>>> aCircle `isEnclosing` <$> ([(0,0),(1,0),(2,0)] :: [XY])
[True, False, False]
    -}
  isEnclosing :: (Metric p, Affine p) => c -> p -> Bool
  isEnclosing c p = location c p == Inside

instance ClosedCurve c => ClosedCurve (Maybe c) where
  location = maybe (const Outside) location
  isEnclosing = maybe (const False) isEnclosing

-- | The type representing the relation /belongs to/ between an affine point
-- | and a curve.
data PointLocation = Inside
                   | Outside
                   | OnCurve
  deriving (Show, Eq)

------------------------------------------------------------
-- | The type representing lower-left and upper-right corners of `Figure`'s bounding box.
type Box = ((Min Double, Min Double), (Max Double, Max Double))

pointBox :: Affine p => p -> Box
pointBox p = ((Min x, Min y), (Max x, Max y))
    where (x, y) = xy p

-- | Class representing the interface for a figure on a chart
class Trans a => Figure a where
  {-# MINIMAL isTrivial, box #-}

  -- | A rectangular bounding box enclosing a figure.
  box :: a -> Box

  -- | Returns `True` is figure is trivial in a certain sence.
  isTrivial :: a -> Bool

  -- | Returns `True` is figure is not trivial in a certain sence.
  isNontrivial :: a -> Bool
  isNontrivial = not . isTrivial

  -- | Returns a refference point using for superposing and figure locataion.
  refPoint :: a -> Cmp
  refPoint = left . lower . corner

instance Bounded Double where
  minBound = -1/0
  maxBound = 1/0

instance
  {-# OVERLAPPABLE #-}
  (Functor t, Foldable t, Figure f) => Figure (t f) where
  box = foldMap box
  isTrivial = null

-- | Returns the total width of the figure.
figureWidth :: Figure f => f -> Double
figureWidth f = let ((Min xmin, _), (Max xmax, _)) = box f
          in abs $ xmax - xmin

-- | Returns the total height of the figure.
figureHeight :: Figure f => f -> Double
figureHeight f = let ((_, Min ymin), (_, Max ymax)) = box f
           in abs $ ymax - ymin

-- | Returns four corners of the figure.
corner :: Figure f => f -> ((Cmp, Cmp), (Cmp, Cmp))
corner p = ((x1:+y2, x2:+y2),(x1:+y1, x2:+y1))
    where ((Min x1,Min y1),(Max x2, Max y2)) = box p

-- | Selectors for corners.
lower :: (a,b) -> b
lower = snd

-- |
right :: (a,b) -> b
right = snd

-- |
upper :: (a,b) -> a
upper = fst

-- |
left :: (a,b) -> a
left  = fst

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

max_stars_repo_name: MethaHardworker/geometry

max_issues_repo_name: MethaHardworker/geometry

max_forks_repo_name: MethaHardworker/geometry


text:
-----------------------------------------------------------------------------
-- |
-- Module     : Numeric.LinearAlgebra.Matrix.ST
-- Copyright  : Copyright (c) , Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- Mutable matrices in the ST monad.

module Numeric.LinearAlgebra.Matrix.ST (
    -- * Mutable matrices
    STMatrix,
    IOMatrix,
    create,

    -- * Read-only matrices
    RMatrix(..),

    -- * Creating new matrices
    new_,
    new,

    -- * Copying matrices
    newCopy,
    copyTo,

    -- * Matrix views
    withSlice,
    withTakeRows,
    withDropRows,
    withSplitRowsAt,
    withTakeCols,
    withDropCols,
    withSplitColsAt,
    withSliceM,
    withTakeRowsM,
    withDropRowsM,
    withSplitRowsAtM,
    withTakeColsM,
    withDropColsM,
    withSplitColsAtM,

    -- * Matrix rows and columns
    rowTo,
    unsafeRowTo,
    setRow,
    unsafeSetRow,
    withCol,

    withColM,
    withColsM,

    swapRows,
    unsafeSwapRows,
    swapCols,
    unsafeSwapCols,

    -- * Matrix diagonals
    diagTo,
    setDiag,

    -- * Reading and writing matrix elements
    read,
    write,
    modify,
    getIndices,
    getElems,
    getElems',
    getAssocs,
    getAssocs',
    setElems,
    setAssocs,

    -- * List-like operations
    mapTo,
    zipWithTo,

    -- * Matrix math operations
    shiftDiagM_,
    shiftDiagWithScaleM_,
    addTo,
    subTo,
    scaleM_,
    addWithScaleM_,
    scaleRowsM_,
    scaleColsM_,
    negateTo,
    conjugateTo,

    -- * Linear algebra
    transTo,
    conjTransTo,
    rank1UpdateM_,

    -- ** Matrix-Vector multiplication
    mulVectorTo,
    mulVectorWithScaleTo,
    addMulVectorWithScalesM_,

    -- ** Matrix-Matrix multiplication
    mulMatrixTo,
    mulMatrixWithScaleTo,
    addMulMatrixWithScalesM_,

    -- * Conversions between mutable and immutable matrices
    freeze,

    -- * Vector views of matrices
    maybeWithVectorM,

    -- * Matrix views of vectors
    withFromVector,
    withFromCol,
    withFromRow,

    withFromVectorM,
    withFromColM,
    withFromRowM,

    -- * Unsafe operations
    unsafeCopyTo,
    unsafeAddWithScaleM_,

    ) where

import Prelude()
import Numeric.LinearAlgebra.Matrix.STBase

ext: hs

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
module Utils.Stats (
        pricesAtInterval
    ,   baseStats
    ,   pricePercChange
    ,   slope'
    ,   std'
    ,   mean
    ,   roundToNearest
) where

import qualified Data.List                   as L
import           Data.Maybe
import qualified Data.Vector.Unboxed         as U
import           Protolude
import           Statistics.LinearRegression
import           Types                       (BasicStats (..))

-- i wonder if this could be written with foldable fns
intervalPrices :: Num a  => Int -> [a] -> [[a]]
intervalPrices _ [] = []
intervalPrices i p = let (x, xs) = splitAt i p
                        in x : intervalPrices i xs

pricesAtInterval :: (Num a ) => Int -> [a] -> Maybe [a]
pricesAtInterval i p =
    let priceGroups = intervalPrices i $ reverse p
    in traverse head $ reverse priceGroups

mean :: [Double] -> Maybe Double
mean [] = Nothing
mean x  = Just $  sum x / fromIntegral (length x)

std' :: [Double] -> Maybe Double
std' [] = Nothing
std' x  =
    let maybeμ     = mean x
        length'    = fromIntegral (length x)
        maybeVar   = case maybeμ of
                        Nothing -> Nothing
                        Just μ  -> Just $ foldr(\c acc -> (μ - c)**2 + acc) 0 x
    in case maybeVar of
        Nothing  -> Nothing
        Just var -> Just $ sqrt (var / length')

slope' :: [Double] -> Maybe Double
slope' [] = Nothing
slope' y  =
    let length' = fromIntegral $ L.length y :: Double
        xs      = U.fromList [0..length']
        ys      = U.fromList y
        (_, beta) = linearRegression xs ys
    in Just beta

roundToNearest :: Int -> Double -> Double
roundToNearest n x = fromInteger (round (x * 10^n)) / 10^n

pricePercChange :: [Double] -> Maybe Double
pricePercChange [] = Nothing
pricePercChange x  =
    let maxPrice = maximum x -- maximum price for say past month
        currentAvgPrice' = mean $ reverse $ take 3 (reverse x)-- avg price for past 3 days
    in case currentAvgPrice' of
        Nothing              -> Nothing
        Just currentAvgPrice -> Just $ ((maxPrice - currentAvgPrice) / maxPrice) * 100

baseStats :: [Double] -> Maybe BasicStats
baseStats [] = Nothing
baseStats x  =
    let x'            = reverse x
        reverseTake y = reverse . take y -- i dont want to use drop
        direction     = roundToNearest 2 $ (* 1e5) $ fromMaybe 0 $ slope' (reverseTake 5 x')
        std''         = roundToNearest 2 $ (* 1e5) $ fromMaybe 0 $ std' (reverseTake 7 x')
        change'       = roundToNearest 2 $ fromMaybe 0  $ pricePercChange $ reverseTake 24 x'
    in Just $ BasicStats direction change' std''

max_stars_repo_name: epicallan/crypto-bot

max_issues_repo_name: epicallan/crypto-bot-hs

max_forks_repo_name: epicallan/crypto-bot-hs


text:
module Statistics.Quantile.Bench.Accuracy where

import qualified Data.Vector as V

import Statistics.Quantile.Types
import Statistics.Quantile.Util

import Statistics.Sample

import System.IO

err :: Double
    -> Double
    -> Double
err true estimate =
  let e = abs (true - estimate)
  in e * e

selectorAccuracy :: Stream IO
                 -> Quantile
                 -> Double
                 -> Selector IO
                 -> IO Double
selectorAccuracy src q true (Selector select) =
  err true <$> select q src

benchAccuracy :: Int
              -> FilePath
              -> Quantile
              -> Double
              -> Selector IO
              -> IO Deviation
benchAccuracy n fp q true s = do
  as <- V.replicateM n accuracy
  pure $ Deviation (mean as) (stdDev as)
  where accuracy = do
          h <- openFile fp ReadMode
          r <- selectorAccuracy (streamHandle h) q true s
          hClose h
          pure r

max_stars_repo_name: fractalcat/slides

max_issues_repo_name: fractalcat/slides

max_forks_repo_name: fractalcat/slides


text:
import Numeric.GSL.Polynomials
import Data.Complex

*Main> mapM_ print $ polySolve [0,2,-3,1]
(-5.421010862427522e-20) :+ 0.0
2.000000000000001 :+ 0.0
0.9999999999999996 :+ 0.0

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE TypeFamilies  #-}

-- | Provides spherical harmonic models of scalar-valued functions.
module Math.SphericalHarmonics
(
  SphericalHarmonicModel
, sphericalHarmonicModel
, scaledSphericalHarmonicModel
, evaluateModel
, evaluateModelCartesian
, evaluateModelGradient
, evaluateModelGradientCartesian
, evaluateModelGradientInLocalTangentPlane
)
where

import Data.Complex
import Data.VectorSpace hiding (magnitude)
import Math.SphericalHarmonics.AssociatedLegendre
import Numeric.AD

-- | Represents a spherical harmonic model of a scalar-valued function.
data SphericalHarmonicModel a = SphericalHarmonicModel [[(a, a)]]
  deriving (Functor)

-- | Creates a spherical harmonic model.
-- Result in an error if the length of the list is not a triangular number.
sphericalHarmonicModel :: (Fractional a) => [[(a, a)]] -- ^ A list of g and h coefficients for the model
                       -> SphericalHarmonicModel a -- ^ The spherical harmonic model
sphericalHarmonicModel cs | valid = SphericalHarmonicModel cs
                          | otherwise = error "The number of coefficients is not a triangular number."
  where
    valid = and $ zipWith (==) (fmap length cs) [1..length cs]

-- | Creates a spherical harmonic model, scaling coefficients for the supplied reference radius.
-- Result in an error if the length of the list is not a triangular number.
scaledSphericalHarmonicModel :: (Fractional a) => a -- ^ The reference radius
                       -> [[(a, a)]] -- ^ A list of g and h coefficients for the model
                       -> SphericalHarmonicModel a -- ^ The spherical harmonic model
scaledSphericalHarmonicModel r cs = sphericalHarmonicModel cs'
  where
    cs' = normalizeReferenceRadius r cs

instance(Fractional a, Eq a) => AdditiveGroup (SphericalHarmonicModel a) where
  zeroV = SphericalHarmonicModel [[(0,0)]]
  negateV = fmap negate
  (SphericalHarmonicModel m1) ^+^ (SphericalHarmonicModel m2) = SphericalHarmonicModel (combineCoefficients m1 m2)
    where
      combineCoefficients [] cs = cs
      combineCoefficients cs [] = cs
      combineCoefficients (c1:cs1) (c2:cs2) = zipWith addPairs c1 c2 : combineCoefficients cs1 cs2
      addPairs (g1, h1) (g2, h2) = (g1 + g2, h1 + h2)

instance (Fractional a, Eq a) => VectorSpace (SphericalHarmonicModel a) where
  type Scalar (SphericalHarmonicModel a) = a
  x *^ m = fmap (* x) m

normalizeReferenceRadius :: (Fractional a) => a -> [[(a, a)]] -> [[(a, a)]]
normalizeReferenceRadius r = zipWith (fmap . mapWholePair . transform) [0 :: Int ..]
  where
    transform n = (* (r ^ (2 + n)))

-- | Computes the scalar value of the spherical harmonic model at a specified spherical position.
evaluateModel :: (RealFloat a, Ord a) => SphericalHarmonicModel a -- ^ Spherical harmonic model
              -> a -- ^ Spherical radius
              -> a -- ^ Spherical colatitude (radian)
              -> a -- ^ Spherical longitude (radian)
              -> a -- ^ Model value
evaluateModel m r colat lon = evaluateModel' m r (cos colat) (cis lon)

-- | Computes the scalar value of the spherical harmonic model at a specified Cartesian position.
evaluateModelCartesian :: (RealFloat a, Ord a) => SphericalHarmonicModel a -- ^ Spherical harmonic model
                       -> a -- ^ X position
                       -> a -- ^ Y position
                       -> a -- ^ Z position
                       -> a -- ^ Model value
evaluateModelCartesian m x y z = evaluateModel' m r cosColat cisLon
  where
    r = sqrt $ (x*x) + (y*y) + (z*z)
    cosColat = z / r
    cisLon = normalize $ mkPolar x y

evaluateModel' :: (RealFloat a, Ord a) => SphericalHarmonicModel a
               -> a -- r
               -> a -- cosColat
               -> Complex a -- cisLon
               -> a
evaluateModel' (SphericalHarmonicModel cs) r cosColat cisLon = sum $ zipWith (*) (iterate (/ r) (recip r)) (zipWith evaluateDegree [0..] cs)
  where
    sines = 1 : iterate (* cisLon) cisLon
    evaluateDegree n cs' = sum $ zipWith3 evaluateOrder (fmap (schmidtSemiNormalizedAssociatedLegendreFunction n) [0..n]) cs' sines
    evaluateOrder p (g, h) cisMLon = ((g * realPart cisMLon) + (h * imagPart cisMLon)) * (p (cosColat))

-- | Computes the gradient of the scalar value of the spherical harmonic model, in spherical coordinates, at a specified location.
evaluateModelGradient :: (RealFloat a, Ord a) => SphericalHarmonicModel a -- ^ Spherical harmonic model
                      -> a -- ^ Spherical radius
                      -> a -- ^ Spherical colatitude (radian)
                      -> a -- ^ Spherical longitude (radian)
                      -> (a, a, a) -- ^ Radial, colatitudinal, and longitudinal components of gradient
evaluateModelGradient model r colat lon = makeTuple . fmap negate $ modelGrad [r, colat, lon]
  where
    modelGrad = grad (\[r', c', l'] -> evaluateModel (fmap auto model) r' c' l')

-- | Computes the gradient of the scalar value of the spherical harmonic model at a specified location, in Cartesian coordinates.
-- The result is expressed in right-handed coordinates centered at the origin of the sphere, with the positive Z-axis piercing the
-- north pole and the positive x-axis piercing the reference meridian.
evaluateModelGradientCartesian :: (RealFloat a, Ord a) => SphericalHarmonicModel a -- ^ Spherical harmonic model
                               -> a -- ^ X position
                               -> a -- ^ Y position
                               -> a -- ^ Z position
                               -> (a, a, a) -- X, Y, and Z components of gradient
evaluateModelGradientCartesian model x y z = makeTuple . fmap negate $ modelGrad [x, y, z]
  where
    modelGrad = grad (\[x', y', z'] -> evaluateModelCartesian (fmap auto model) x' y' z')

-- | Computes the gradient of the scalar value of the spherical harmonic model at a specified location, in Cartesian coordinates.
-- The result is expressed in a reference frame locally tangent to the sphere at the specified location.
evaluateModelGradientInLocalTangentPlane :: (RealFloat a, Ord a) => SphericalHarmonicModel a -- ^ Spherical harmonic model
                                         -> a -- ^ Spherical radius
                                         -> a -- ^ Spherical colatitude (radian)
                                         -> a -- ^ Spherical longitude (radian)
                                         -> (a, a, a) -- ^ East, North, and up components of gradient
evaluateModelGradientInLocalTangentPlane model r colat lon = (e, n, u)
  where
    (r', colat', lon') = evaluateModelGradient model r colat lon
    e = lon' / (r * sin colat)
    n = -colat' / r -- negated because the colatitude increase southward
    u = r'

normalize :: (RealFloat a) => Complex a -> Complex a
normalize r@(x :+ y) | isInfinite m' = 0
                     | otherwise = (x * m') :+ (y * m')
  where
    m' = recip . magnitude $ r

mapWholePair :: (a -> b) -> (a, a) -> (b, b)
mapWholePair f (a, b) = (f a, f b)

makeTuple :: [a] -> (a, a, a)
makeTuple [x, y, z] = (x, y, z)

max_stars_repo_name: dmcclean/igrf

max_issues_repo_name: dmcclean/igrf

max_forks_repo_name: dmcclean/igrf


text:
{- |
Description :  Boundary mutation model Markov process
Copyright   :  (c) Dominik Schrempf 2017
License     :  GPLv3

Maintainer  :  dominik.schrempf@gmail.com
Stability   :  unstable
Portability :  non-portable (not tested)

This module defines the transition rate matrix of the boundary mutation model.

-}

module BndModel
  ( BMModel(..)
  , MutModel
  , Heterozygosity
  , createBMM
  , createBMMNormalized
  , scaleTreeToBMM
  , getBMMInfoStr
  ) where

import           BndState
import           Data.Maybe
import           DNAModel              hiding (rateMatrix)
import           Numeric.LinearAlgebra
import           RateMatrix            hiding (State)
import           RTree
import           Tools

-- | The boundary mutation models uses an underlying mutation model.
type MutModel = DNAModel

-- | Define a heterozygosity to make function definitions clearer.
type Heterozygosity = Double

-- | A boundary mutation model is defined by its rate matrix. However, it is
-- convenient to keep the underlying mutation model, the population size the
-- heterozygosity and the stationary distribution at hand.
data BMModel = BMModel { bmmRateMatrix     :: RateMatrix
                       , bmmMutModel       :: MutModel
                       , bmmPopSize        :: PopSize
                       , bmmHeterozygosity :: Heterozygosity
                       , bmmStationaryDist :: StationaryDist }

-- Get the mutation rate from one allele to another.
mutRate :: MutModel -> Allele -> Allele -> Double
mutRate (DNAModel m _) a b = m ! i ! j
  where i = fromEnum a
        j = fromEnum b

-- Get state frequency for specific allele.
stationaryFreq :: StationaryDist -> Allele -> Double
stationaryFreq f a = f ! fromEnum a

-- Calculate the rate of frequency shifts.
moranCoef :: PopSize -> Int -> Double
moranCoef n i = iD * (nD - iD) / nD
  where iD = fromIntegral i
        nD = fromIntegral n

-- The transition rate from one boundary state to another.
rate :: MutModel -> State -> State -> Double
rate m s t
  | not $ connected s t = 0.0
  | otherwise           = rate' s t
  where rate' (Ply n i _ _) _ = moranCoef n i
        rate' (Bnd _ a) (Ply _ _ b c)
          | a == b    = mutRate m a c
          | a == c    = mutRate m a b
          | otherwise = error "Cannot compute rate between states."
        rate' _ _ = error "Cannot compute rate between states."

-- The transition rate from one state (index) to another.
rateById :: MutModel -> PopSize -> Int -> Int -> Double
rateById m n i j = rate m s t
  where s = idToState n i
        t = idToState n j

-- The build function (see below) has a weird way of assigning entries to
-- indices. The indices have to be the same data type as the entries. This is
-- just a helper function that changes the indices from Double to Int.
rateByDouble :: MutModel -> PopSize -> Double -> Double -> Double
rateByDouble m n x y = rateById m n (round x) (round y)

-- The rate matrix of the boundary mutation model. The dimension with four
-- alleles will be 4 + 6*(N-1).
rateMatrix :: MutModel -> PopSize -> RateMatrix
rateMatrix m n = setDiagonal $ build (s,s) (rateByDouble m n)
  where s = stateSpaceSize n

-- Normalize the rate matrix such that on average one event happens per unit
-- time.
normalizedRateMatrix :: MutModel -> PopSize -> RateMatrix
normalizedRateMatrix m n = normalizeRates f' m'
  where m' = rateMatrix m n
        -- f  = dnaModelSpecGetStateFreqVec (dnaModelSpec m)
        f' = stationaryDist m n

-- Calculate the heterozygosity at stationarity.
theta :: MutModel -> Heterozygosity
theta (DNAModel m s) = f <.> rDiagZero #> f
  where f         = dnaModelSpecGetStateFreqVec s
        e         = toExchMatrix m f
        (r, _)    = matrixSeparateSymSkew e
        -- The summation excludes the diagonal (a /= b).
        rDiagZero = matrixSetDiagToZero r

-- The normalization constant of the stationary distribution.
norm :: MutModel -> PopSize -> Double
norm m n = 1.0 + harmonic (n-1) * theta m

-- Get entries of the stationary measure (not normalized) for a boundary model state.
stationaryMeasEntry :: MutModel -> State -> Double
stationaryMeasEntry m s =
  let f   = dnaModelSpecGetStateFreqVec (dnaModelSpec m) in
  case s of
    (Bnd _ a)     -> f ! fromEnum a
    (Ply n i a b) -> piA * mAB / (fromIntegral n - fromIntegral i)
                                      + piB * mBA / fromIntegral i
      where piA = stationaryFreq f a
            mAB = mutRate m a b
            piB = stationaryFreq f b
            mBA = mutRate m b a

stationaryMeasEntryById :: MutModel -> PopSize -> Int -> Double
stationaryMeasEntryById m n i = stationaryMeasEntry m s
  where s = idToState n i

stationaryMeasEntryByDouble :: MutModel -> PopSize -> Double -> Double
stationaryMeasEntryByDouble m n x = stationaryMeasEntryById m n (round x)

-- Get the stationary distribution of a boundary mutation model.
stationaryDist :: MutModel -> PopSize -> StationaryDist
stationaryDist m n = scale (1.0/norm m n) $
  build s (stationaryMeasEntryByDouble m n)
  where s = stateSpaceSize n

-- Normalize the mutation coefficients such that the heterozygosity matches a
-- given level (see Eq. 12.14 in my thesis).
normalizeToTheta :: MutModel -> PopSize -> Heterozygosity -> MutModel
normalizeToTheta mo@(DNAModel m _) n h =
  mo { dnaRateMatrix = scale (h / (t * (1.0 - c * h))) m }
  where
    -- The heterozygosity of the boundary mutation model.
    t = theta mo
    -- The branch length multiplicative factor introduced by the coalescent.
    c = harmonic (n-1)

-- | Create a boundary mutation model using the minimal number of necessary
-- ingredients.
createBMM :: DNAModel -> PopSize -> Heterozygosity -> BMModel
createBMM m n h = BMModel rm m' n h f
  where m' = normalizeToTheta m n h
        rm = normalizedRateMatrix m' n
        f  = stationaryDist m n

-- | Create a boundary mutation model without providing a heterozygosity. This
-- means that the mutation model has to be normalized already.
createBMMNormalized :: DNAModel -> PopSize -> BMModel
createBMMNormalized m n = BMModel rm m n h f
  where rm = normalizedRateMatrix m n
        h  = theta m
        f  = stationaryDist m n

-- | The branch lengths of threes in the boundary mutation model are not
-- measured in average number of substitutions per site but in average number of
-- mutations or frequency shifts per site. The conversion factor is just the
-- square of the population size. This function converts the branch lengths of a
-- tree.
scaleTreeToBMM :: BMModel -> RTree a Double -> RTree a Double
scaleTreeToBMM bmm = fmap (* nSq)
  where nSq = fromIntegral (bmmPopSize bmm) ** 2

-- | Report the boundary mutation model specifications.
getBMMInfoStr :: BMModel
              -> Maybe Double
              -> Maybe [Double]
              -> String
getBMMInfoStr bmm ma mrs = unlines $
  [ "Population size: " ++ show (bmmPopSize bmm)
  , "Heterozygosity: " ++ show (bmmHeterozygosity bmm)
  , "Mutation model:"
  , getDNAModelInfoStr (bmmMutModel bmm)
  , "Gamma rate heterogeneity: " ++ show (isJust ma) ] ++
  gammaShape ++ gammaMeans
  where
    gammaShape = maybe [] (\a -> ["Shape parameter: " ++ show a]) ma
    gammaMeans = maybe [] (\rs -> ["This corresponds to uniformly distributed rates: " ++ show rs]) mrs

max_stars_repo_name: dschrempf/bmm-simulate

max_issues_repo_name: pomo-dev/bmm-simulate

max_forks_repo_name: pomo-dev/bmm-simulate


text:
{-# LANGUAGE ForeignFunctionInterface #-}
{-# LANGUAGE QuasiQuotes     #-}
{-# LANGUAGE TemplateHaskell #-}
module Main where

import           Codec.Compression.GZip       ( decompress )
import           Control.Monad
import qualified Data.Array.Repa              as R
import qualified Data.Array.Repa.Stencil      as Stencil
import qualified Data.Array.Repa.Stencil.Dim2 as Stencil
import qualified Data.Massiv.Array            as A
import qualified Data.ByteString              as BS
import qualified Data.ByteString.Lazy         as BL
import           Data.List                    as L
import qualified Numeric.LinearAlgebra        as LA
import           Prelude                      as P
import           System.Directory             ( getCurrentDirectory )
import           System.Environment
import           System.IO
import           System.Random                as Rand
import           Control.DeepSeq              as DS

import           Foreign.Ptr                  ( Ptr )
import           System.IO.Unsafe             ( unsafePerformIO )
import qualified Data.Vector.Storable as U    ( unsafeToForeignPtr0, unsafeFromForeignPtr0 )
import qualified Numeric.LinearAlgebra.Devel as U
import           Foreign                      ( mallocForeignPtrArray, withForeignPtr )

main :: IO ()
main = do
    args <- getArgs

    imgs <- loadData
    kernels <- initKernels
    case args of
      "--massiv":_  -> do
          data_ <- return . BL.toStrict =<< getData "train-images-idx3-ubyte.gz"
          let imgs_ :: A.Array A.U A.Ix3 Double
              imgs_ = A.fromLists' A.Par $
                  (getImage_ data_) <$> [0..49999]
              {-# INLINE imgs_ #-}
              kernels_ :: A.Array A.U A.Ix3 Double
              kernels_ = A.computeAs A.U $
                  A.uniformRangeArray (mkStdGen 42) (0, 1) (A.ParN 3) (A.Sz3 20 5 5)
                  -- TODO: convert to standard normal
                  -- distribution using Box-Mueller transform
          runMassiv imgs_ kernels_
      "--hmatrix":_ -> runHmatrix imgs kernels
      "--cublas":_  -> runCublas
      "--repa":_    -> runRepa imgs kernels
      "--cuda":_    -> runCuda
      _             -> putStrLn "Error"

runMassiv :: A.Array A.U A.Ix3 Double -> A.Array A.U A.Ix3 Double -> IO ()
runMassiv imgs kernels = do
--    let stencil = A.makeConvolutionStencilFromKernel . (ks !!)
--        {-# INLINE stencil #-}
--
--        base = A.computeAs A.U $ A.applyStencil A.noPadding (stenc 0) x
--        conv_ = ((\i -> (\s -> A.computeAs A.U $ A.applyStencil A.noPadding s i) <$> stencils) <$> is)
--        conv = (foldl1' (\p ch -> A.computeAs A.U $ A.append' 3 p ch) conv) :: A.Array A.U A.Ix3 Double
    imgs `deepseq` putStrLn "Done"

-- loadDataMassiv :: IO (A.Array A.U A.Ix3 Double)
-- loadDataMassiv = do
--     trainImgs <- getData "train-images-idx3-ubyte.gz"
--     let imgs = BL.toStrict trainImgs
--         imgArr :: A.Array A.U A.Ix3 Double
--         imgArr = foldl1 (\acc i -> A.computeAs A.U $ A.append' 3 acc i) [getImage_ n imgs | n <- [0..49999]]
--     return imgArr

getImage_ :: BS.ByteString -> Int -> [[Double]]
getImage_ imgs n = chunkList 28 [normalize $ BS.index imgs (16 + n*784 + s) | s <- [0..783]]

-- loadDataMassiv :: IO (A.Array A.U A.Ix3 Double)
-- loadDataMassiv = do
--     trainImgs <- getData "train-images-idx3-ubyte.gz"
--     let imgs = BL.toStrict trainImgs
--         imgArr :: A.Array A.U A.Ix3 Double
--         imgArr = foldl1 (\acc i -> A.computeAs A.U $ A.append' 3 acc i) [getImage_ n imgs | n <- [0..49999]]
--     return imgArr

-- getImage_ :: Int -> BS.ByteString -> A.Array A.U A.Ix3 Double
-- getImage_ n imgs = A.resize' (A.Sz (1 A.:> 28 A.:. 28)) $
--     A.fromList A.Par [normalize $ BS.index imgs (16 + n*784 + s) | s <- [0..783]]

-- time:   1.54s user
-- space: ~2GB?
runHmatrix :: [LA.Matrix Double] -> [LA.Matrix Double] -> IO ()
runHmatrix imgs kernels =
    let kcols = LA.fromRows $ P.map LA.flatten kernels
        icols = (LA.tr . im2col 5 5 1 1) <$> imgs
        res = P.map (kcols LA.<>) icols
     in res `deepseq` putStrLn $ "Done, size: " P.++ (show . length $ res)

    --let dCircKernel = doublyBlockedCirculant . zeroPad (28, 28) . flipKernel $ k
    --sequence_ $
    --    (LA.dispShort 10 10 3) <$> [dCircKernel LA.<> img | img <- imgs]

-- time:   47.93s user
-- space: ~5GB
runRepa :: [LA.Matrix Double] -> [LA.Matrix Double] -> IO ()
runRepa imgs kernels = print "Done"
--    let ks = hmatrix2repa <$> kernels
--        is = hmatrix2repa <$> imgs
--
--        imgMat = head is
--        kernel = [stencil2| 10 10 10 10 10
--                            10 10 10 10 10
--                            10 10 10 10 10
--                            10 10 10 10 10
--                            10 10 10 10 10 |] :: Stencil.Stencil DIM2 Double
--    res <- forM is $ \i ->
--        forM (replicate 20 kernel) $ \k ->
--            computeP ((mapStencil2 (BoundConst 0) k i) :: Array PC5 DIM2 Double)
--    print (res :: [[Array U DIM2 Double]])

runCublas = print ""

runCuda = print ""

--convolveP :: Array U DIM2 Double  -- ^ 28x28    Input image       (Unboxed Vector)
--          -> Array C DIM3 Double  -- ^ 20x5x5   Kernel tensor     (Cursored Array)
--          -> Array U DIM3 Double  -- ^ 20x24x24 Convolution layer (Unboxed Vector)

{- Kernel stuff -}
initKernels :: IO [LA.Matrix Double]
initKernels = replicateM 20 $ LA.randn 5 5

flipKernel :: LA.Matrix Double -> LA.Matrix Double
flipKernel = LA.flipud

doublyBlockedCirculant :: LA.Matrix Double  -- ^ input kernel
                       -> LA.Matrix Double
doublyBlockedCirculant m = LA.fromBlocks . circ $ circOfRows m

circOfRows :: LA.Matrix Double -> [LA.Matrix Double]
circOfRows m = [ LA.fromLists $ circ row | row <- rows ]
    where rows = reverse $ LA.toLists m

circ :: [a] -> [[a]]
circ xs = take (length xs) $ iterate (rotate 1) xs -- `(length xs)` is wrong here..

rotate :: Int -> [a] -> [a]
rotate n xs = drop diff xs <> take diff xs
    where diff = length xs - n

zeroPad :: (Int, Int) -> LA.Matrix Double -> LA.Matrix Double
zeroPad (r,c) m = LA.fromBlocks $
    [[LA.konst 0 (r-r', 1), LA.row $ replicate (c-c') 0],
     [m, LA.row $ replicate (c-c') 0]]
         where (r', c') = LA.size m

zeroPadTo :: Int -> Int -> [LA.Matrix Double] -> [LA.Matrix Double]
zeroPadTo r c kernels = [ LA.fromBlocks [[LA.konst 0 (r', 1), LA.row $ replicate r' 0], [k, LA.row $ replicate r' 0]] | k <- kernels]
    where (r', c') = (r-5, c-5)

getSizeZP :: Int -> Int -> (Int, Int)
getSizeZP m n = ((n-m+1)^2, n^2)

{- Load data -}
loadData :: IO [LA.Matrix Double]
loadData = do
    trainImgs <- getData "train-images-idx3-ubyte.gz"
    let imgs = BL.toStrict trainImgs
        trainData = chunkList 10 [getImage n imgs | n <- [0..49999]]
    return $ LA.fromColumns <$> trainData

getImage :: Int -> BS.ByteString -> LA.Vector Double
getImage n imgs = LA.fromList [normalize $ BS.index imgs (16 + n*784 + s) | s <- [0..783]]

getData :: FilePath -> IO BL.ByteString
getData path = do
    currentDir <- getCurrentDirectory
    fileData <- decompress <$> BL.readFile (currentDir P.<> "/mnist_dataset/" P.<> path)
    return fileData

chunkList :: Int -> [a] -> [[a]]
chunkList n xs = takeWhile (not.null) $ unfoldr (Just . splitAt n) xs

normalize :: (Integral a, Floating b) => a -> b
normalize x = fromIntegral x/255

hmatrix2repa :: LA.Matrix Double -> R.Array R.U R.DIM2 Double
hmatrix2repa mat = R.fromListUnboxed (R.Z R.:. (r::Int) R.:. (c::Int)) . LA.toList . LA.flatten $ mat
    where (r, c) = LA.size mat

foreign import ccall unsafe im2col_cpu :: Ptr Double -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> Ptr Double -> IO ()

im2col_c :: Int -> Int -> Int -> Int -> Int -> Int -> Int -> LA.Matrix Double -> LA.Matrix Double
im2col_c channels height width kernelRows kernelColumns strideRows strideColumns dataImg =
    let vec = LA.flatten dataImg
        rowOut = (height - kernelRows) `div` strideRows + 1
        colOut = (width - kernelColumns) `div` strideColumns + 1
        kernelSize = kernelRows * kernelColumns
        numberOfPatches = rowOut * colOut
    in unsafePerformIO $ do
        outPtr <- mallocForeignPtrArray (numberOfPatches * kernelSize * channels)
        let (inPtr, _) = U.unsafeToForeignPtr0 vec

        withForeignPtr inPtr $ \inPtr' ->
            withForeignPtr outPtr $ \outPtr' ->
                im2col_cpu inPtr' channels height width kernelRows kernelColumns strideRows strideColumns outPtr'

        let matVec = U.unsafeFromForeignPtr0 outPtr (numberOfPatches * kernelSize * channels)
        return $ U.matrixFromVector U.RowMajor numberOfPatches (kernelSize * channels) matVec

im2col :: Int -> Int -> Int -> Int -> LA.Matrix Double -> LA.Matrix Double
im2col kernelRows kernelColumns strideRows strideColumns dataIm =
  let channels = 1
      height = LA.rows dataIm
      width  = LA.cols dataIm
  in  im2col_c channels height width kernelRows kernelColumns strideRows strideColumns dataIm


text:
module Main where

import Criterion.Main
import Statistics.Sample.Histogram.Magnitude
import Data.Functor.Identity

main = defaultMain
  [ bgroup "resolution" [ bench "1"  . nf (foldHist 1) $ Identity (1 :: Double)
                        , bench "2"  . nf (foldHist 2) $ Identity (1 :: Double)
                        , bench "3"  . nf (foldHist 3) $ Identity (1 :: Double)
                        ]
  , bgroup "size"       [ bench "0"      $ nf (foldHist 2) ([] :: [Double])
                        , bench "0..9"   $ nf (foldHist 2) [0..9 :: Double]
                        , bench "0..99"  $ nf (foldHist 2) [0..99 :: Double]
                        ]
  ]

max_stars_repo_name: skedgeme/histogram-magnitude

max_issues_repo_name: skedgeme/histogram-magnitude

max_forks_repo_name: skedgeme/histogram-magnitude


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

import qualified ArrayFire       as A
import           Data.Int
import           Data.Word
import           Data.Complex
import           Data.Proxy
import           Foreign.C.Types
import           Test.Hspec

spec :: Spec
spec =
  describe "Signal spec" $ do
    it "Should do FFT in place" $ do
      A.fftInPlace (A.matrix @(Complex Double) (1,1) [[1 :+ 1]]) 10.2
        `shouldReturn` ()
    it "Should do FFT" $ do
      A.fft (A.matrix @(Complex Float) (1,1) [[1 :+ 1]]) 1 1
        `shouldBe` A.matrix @(Complex Float) (1,1) [[1 :+ 1]]

max_stars_repo_name: danielkroeni/arrayfire-haskell

max_issues_repo_name: danielkroeni/arrayfire-haskell

max_forks_repo_name: danielkroeni/arrayfire-haskell


text:
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}

module CircSpec where

import Prelude hiding ((.))
import Test.QuickCheck.Classes
import Test.QuickCheck hiding (output)
import Test.Hspec
import Test.QuickCheck.Checkers

import ConCat.Circ

import ConCat.Category hiding (it)
import ConCat.Rep
import ConCat.Free.VectorSpace hiding (it)
import ConCat.Circuit
import ConCat.Misc
import Data.Complex

instance (Arbitrary a) => Arbitrary (RLC a) where
  arbitrary = oneof [Res <$> arbitrary, Cap <$> arbitrary, Ind <$> arbitrary]

spec = do
  describe "LCirc built out of the Circuit Category" $ do
    it "Rep is coherent" $ {--do
      a <- (arbs @RLC ) 10
      let
        vcheck :: RLC -> IO ()
        vcheck r = do
          print $ "actual: " <> (show r)
          print $ "repr: " <> (show $ repr r)
          print $ "abst . repr " <> (show $ (abst @ RLC) . repr $ r)
      mapM (vcheck) a--}
      property $ ((\(x :: (RLC (Complex R))) -> (abst . repr $ x) == x))
    it "Vectorspace instance is coherent" $ {--do
      a <- (arbs @RLC ) 10
      let
        vcheck :: RLC -> IO ()
        vcheck r = do
          print $ "actual: " <> (show $ r)
          print $ "toVec: " <> (show $ (toV @Double @RLC) r)
          print $ "unVec . toVec " <> (show $ ((unV @Double @RLC) . (toV @Double @RLC) $ r))
      mapM vcheck a--}
      property $ ((\(x :: RLC (Complex R)) -> ((unV @Double @(RLC (Complex R))) . (toV @Double @(RLC (Complex R))) $ x) == x))
    it "can generate a simple resistor component" $ do
      let
        --r = namedC "r" @(RLC :> RLC)
        r0 :: CR
        r0 = (Res 0)
        r2 :: CR
        r2 = (Res 2)
        rn = addV r2 r0
        --brc = (fmapC @(:>) @Blackbox) r0
      --print r
      --print r0
      --print brc
      rn `shouldBe` r2

type CR = RLC (Complex R)

max_stars_repo_name: eco-energy/LCirc

max_issues_repo_name: eco-energy/LCirc

max_forks_repo_name: eco-energy/LCirc


text:
{-# LANGUAGE BangPatterns #-}
module FokkerPlanck.MonteCarlo
  ( -- runMonteCarloFourierCoefficients
    runMonteCarloFourierCoefficientsGPU
  , solveMonteCarloR2S1
  ) where

import           Array.UnboxedArray              as UA
import           Control.DeepSeq
import           Control.Monad                   as M
import           Control.Concurrent.Async
import           Control.Monad.Parallel          as MP
import           Data.Array.Accelerate           (constant, use)
import           Data.Array.Accelerate.LLVM.PTX
import           Data.Array.Repa                 as R
import           Data.Binary
import           Data.ByteString.Lazy            as BL
import           Data.Complex
import           Data.DList                      as DL
import           Data.Ix
import           Data.List                       as L
import           Data.Vector.Unboxed             as VU
import           FokkerPlanck.BrownianMotion
import           FokkerPlanck.FourierSeries
import           FokkerPlanck.FourierSeriesGPU
import           FokkerPlanck.Histogram
import           Foreign.CUDA.Driver             as CUDA
import           Statistics.Distribution
import           Statistics.Distribution.Laplace
import           Statistics.Distribution.Normal
import           Statistics.Distribution.Poisson
import           Statistics.Distribution.Uniform
import           System.Directory
import           System.IO                       as IO
import           System.Random.MWC
import           Text.Printf
import           Utils.List
import           Utils.Parallel
import           Utils.Time
import GHC.Float

{-# INLINE runBatch #-}
runBatch :: Int -> Int -> Int -> (Int -> a -> [GenIO] -> IO a) -> a -> IO a
runBatch !numGens !numTrails !batchSize monterCarloHistFunc !init = do
  let (!numBatch, !numLeftover) = divMod numTrails batchSize
  printf
    "%d trails and batch size is %d, %d batch in total.\n"
    numTrails
    batchSize
    (numBatch +
     if numLeftover > 0
       then 1
       else 0)
  gensList <- M.replicateM numBatch (M.replicateM numGens createSystemRandom)
  x <-
    M.foldM
      (\y gens -> do
         printCurrentTime ""
         monterCarloHistFunc batchSize y gens)
      init
      gensList
  if numLeftover > 0
    then do
      gens <- M.replicateM numGens createSystemRandom
      monterCarloHistFunc numLeftover x gens
    else return x

{-# INLINE computeHistogramFromMonteCarloParallel #-}
computeHistogramFromMonteCarloParallel ::
     (NFData hist, Binary hist, Show hist)
  => FilePath
  -> (GenIO -> IO particle)
  -> ([particle] -> IO hist)
  -> (hist -> hist -> hist)
  -> Int
  -> hist
  -> [GenIO]
  -> IO hist
computeHistogramFromMonteCarloParallel !filePath pointsGenerator histFunc addHist !n !initHist !gens = do
  xs <- MP.mapM (M.replicateM (div n . L.length $ gens) . pointsGenerator) gens
  let tmpFilePath = (filePath L.++ "_tmp")
  hist <- L.foldl' addHist initHist <$> M.mapM histFunc xs
  unless
    (L.null filePath)
    (do encodeFile tmpFilePath hist
        copyFile tmpFilePath filePath)
  return hist

{-# INLINE computeHistogramFromMonteCarloParallelSingleGPU #-}
computeHistogramFromMonteCarloParallelSingleGPU ::
     (NFData hist, Binary hist, Show hist)
  => FilePath
  -> (GenIO -> IO particle)
  -> ([particle] -> hist)
  -> (hist -> hist -> hist)
  -> Int
  -> Int
  -> hist
  -> [GenIO]
  -> IO hist
computeHistogramFromMonteCarloParallelSingleGPU !filePath pointsGenerator histFunc addHist !deviceID !n !initHist !gens = do
  xs <- MP.mapM (M.replicateM (div n . L.length $ gens) . pointsGenerator) gens
  let tmpFilePath = (filePath L.++ "_tmp")
      !hist = addHist initHist . histFunc . L.concat $ xs
  unless
    (L.null filePath)
    (do encodeFile tmpFilePath hist
        copyFile tmpFilePath filePath)
  return hist

{-# INLINE computeHistogramFromMonteCarloParallelMultipleGPU #-}
computeHistogramFromMonteCarloParallelMultipleGPU ::
     (NFData hist, Binary hist, Show hist, NFData particle)
  => FilePath
  -> (GenIO -> IO particle)
  -> (PTX -> [particle] -> hist)
  -> (hist -> hist -> hist)
  -> [PTX]
  -> Int
  -> hist
  -> [GenIO]
  -> IO hist
computeHistogramFromMonteCarloParallelMultipleGPU !filePath pointsGenerator histFunc addHist !ptxs !n !initHist !gens = do
  xs <-
    L.concat <$>
    mapConcurrently
      (M.replicateM (div n . L.length $ gens) . pointsGenerator)
      gens
  let tmpFilePath = filePath L.++ "_tmp"
      !hist =
        L.foldl' addHist initHist .
        parZipWith rdeepseq histFunc ptxs . divideListN (L.length ptxs) $
        xs
  unless
    (L.null filePath)
    (do encodeFile tmpFilePath hist
        copyFile tmpFilePath filePath)
  return hist

-- {-# INLINE runMonteCarloFourierCoefficients #-}
-- runMonteCarloFourierCoefficients ::
--      Int
--   -> Int
--   -> Int
--   -> Double
--   -> Double
--   -> Double
--   -> Double
--   -> [Double]
--   -> [Double]
--   -> [Double]
--   -> [Double]
--   -> Double
--   -> FilePath
--   -> Histogram (Complex Double)
--   -> IO (Histogram (Complex Double))
-- runMonteCarloFourierCoefficients !numGens !numTrails !batchSize !thetaSigma !scaleSigma !maxScale !tao !phiFreqs !rhoFreqs !thetaFreqs !rFreqs !deltaLog !filePath !initHist = do
--   when
--     (maxScale <= 1)
--     (error $
--      printf "runMonteCarloFourierCoefficients error: maxScale(%f) <= 1" maxScale)
--   let !thetaDist = normalDistrE 0 thetaSigma
--       !scaleDist = normalDistrE 0 scaleSigma
--       pointsGenerator = generatePath thetaDist scaleDist maxScale tao 1
--       histFunc =
--         computeFourierCoefficients
--           phiFreqs
--           rhoFreqs
--           thetaFreqs
--           rFreqs
--           (log maxScale)
--           deltaLog
--       monterCarloHistFunc =
--         computeHistogramFromMonteCarloParallel
--           filePath
--           pointsGenerator
--           histFunc
--           addHistogramUnsafe
--   hist <- runBatch numGens numTrails batchSize monterCarloHistFunc initHist
--   return hist

{-# INLINE runMonteCarloFourierCoefficientsSingleGPU #-}
runMonteCarloFourierCoefficientsSingleGPU ::
     Int
  -> Int
  -> Int
  -> Int
  -> Double
  -> Double
  -> Double
  -> Double
  -> [Double]
  -> [Double]
  -> [Double]
  -> [Double]
  -> Double -> Double
  -> FilePath
  -> Histogram (Complex Double)
  -> IO (Histogram (Complex Double))
runMonteCarloFourierCoefficientsSingleGPU !deviceID !numGens !numTrails !batchSize !thetaSigma !scaleSigma !maxScale !tao !phiFreqs !rhoFreqs !thetaFreqs !rFreqs !deltaLog !initScale !filePath !initHist = do
  when
    (maxScale <= 1)
    (error $
     printf "runMonteCarloFourierCoefficients error: maxScale(%f) <= 1" maxScale)
  gen <- createSystemRandom
  initialise []
  dev <- device deviceID
  ctx <- CUDA.create dev []
  ptx <- createTargetFromContext ctx
  let !thetaDist = normalDistrE 0 thetaSigma
      !scaleDist = normalDistrE 0 scaleSigma
      -- !scaleDist = uniformDistrE (- (log maxScale)) (log maxScale)
      !poissonDist = poissonE 0
      !freqArr =
        computeFrequencyArray
          (L.map double2Float phiFreqs)
          (L.map double2Float rhoFreqs)
          (L.map double2Float thetaFreqs)
          (L.map double2Float rFreqs)
      pointsGenerator =
        generatePath
          thetaDist
          scaleDist
          poissonDist
          (maxScale ^ 2)
          maxScale
          tao
          initScale
          1
      histFuncSingleGPU =
        computeFourierCoefficientsGPU
          phiFreqs
          rhoFreqs
          thetaFreqs
          rFreqs
          (use freqArr)
          (constant . double2Float $ maxScale)
          (constant . double2Float . log $ maxScale)
          ptx
      monterCarloHistFuncSingleGPU =
        computeHistogramFromMonteCarloParallelSingleGPU
          filePath
          pointsGenerator
          histFuncSingleGPU
          addHistogramUnsafe
          deviceID
  hist <-
    runBatch numGens numTrails batchSize monterCarloHistFuncSingleGPU initHist
  -- deepseq hist (destroy ctx)
  return hist

-- {-# INLINE runMonteCarloFourierCoefficientsMultipleGPU #-}
-- runMonteCarloFourierCoefficientsMultipleGPU ::
--      [Int]
--   -> Int
--   -> Int
--   -> Int
--   -> Double
--   -> Double
--   -> Double
--   -> Double
--   -> [Double]
--   -> [Double]
--   -> [Double]
--   -> [Double]
--   -> Double
--   -> Double
--   -> FilePath
--   -> Histogram (Complex Double)
--   -> IO (Histogram (Complex Double))
-- runMonteCarloFourierCoefficientsMultipleGPU !deviceIDs !numGens !numTrails !batchSize !thetaSigma !scaleSigma !maxScale !tao !phiFreqs !rhoFreqs !thetaFreqs !rFreqs !deltaLog !initScale !filePath !initHist = do
--   when
--     (maxScale <= 1)
--     (error $
--      printf "runMonteCarloFourierCoefficients error: maxScale(%f) <= 1" maxScale)
--   initialise []
--   devs <- M.mapM device deviceIDs
--   ctxs <- M.mapM (\dev -> CUDA.create dev []) devs
--   ptxs <- M.mapM createTargetFromContext ctxs
--   let !thetaDist = normalDistrE 0 thetaSigma
--       !scaleDist = normalDistrE 0 scaleSigma
--       !freqArr = computeFrequencyArray phiFreqs rhoFreqs thetaFreqs rFreqs
--       pointsGenerator = generatePath thetaDist scaleDist maxScale tao initScale
--       histFuncMultipleGPU =
--         computeFourierCoefficientsGPU
--           phiFreqs
--           rhoFreqs
--           thetaFreqs
--           rFreqs
--           (log maxScale)
--           deltaLog
--           (use freqArr)
--           (constant maxScale)
--           (constant (log maxScale))
--           (constant (deltaLog :+ 0))
--       monterCarloHistFuncMultipleGPU =
--         computeHistogramFromMonteCarloParallelMultipleGPU
--           filePath
--           pointsGenerator
--           histFuncMultipleGPU
--           addHistogramUnsafe
--           ptxs
--   hist <-
--     runBatch numGens numTrails batchSize monterCarloHistFuncMultipleGPU initHist
--   deepseq hist (M.mapM destroy ctxs)
--   return hist

{-# INLINE runMonteCarloFourierCoefficientsGPU #-}
runMonteCarloFourierCoefficientsGPU ::
     [Int]
  -> Int
  -> Int
  -> Int
  -> Double
  -> Double
  -> Double
  -> Double
  -> Double
  -> Double
  -> [Double]
  -> [Double]
  -> [Double]
  -> [Double]
  -> Double
  -> Double
  -> FilePath
  -> IO (Histogram (Complex Double))
runMonteCarloFourierCoefficientsGPU !deviceIDs !numGens !numTrails !batchSize !thetaSigma !scaleLambda !poissonLambda !sigma !tao !deltaT !phiFreqs !rhoFreqs !thetaFreqs !rFreqs !periodEnv !stdR2 !filePath = do
  when
    (periodEnv <= 1)
    (error $
     printf
       "runMonteCarloFourierCoefficients error: periodEnv(%f) <= 1"
       periodEnv)
  initialise []
  devs <- M.mapM device deviceIDs
  ctxs <- M.mapM (\dev -> CUDA.create dev []) devs
  ptxs <- M.mapM createTargetFromContext ctxs
  let !initHist =
        emptyHistogram
          [ L.length phiFreqs
          , L.length rhoFreqs
          , L.length thetaFreqs
          , L.length rFreqs
          ]
          0
      !thetaDist = normalDistrE 0 (thetaSigma * sqrt deltaT)
      !scaleDist = normalDistrE 0 (scaleLambda * sqrt deltaT)
      !poissonDist = poissonE (poissonLambda / deltaT)
      !freqArr =
        computeFrequencyArray
          (L.map double2Float phiFreqs)
          (L.map double2Float rhoFreqs)
          (L.map double2Float thetaFreqs)
          (L.map double2Float rFreqs)
      !maxRho =  sqrt periodEnv / 1 / sqrt 2
      !maxR = sqrt periodEnv / 1 / sqrt 2
      pointsGenerator =
        generatePath
          thetaDist
          scaleDist
          poissonDist
          maxRho
          maxR
          (tao / deltaT)
          deltaT
          stdR2
      -- pointsGenerator =
      --   generatePath'
      --     thetaSigma
      --     scaleDist
      --     poissonLambda
      --     maxRho
      --     maxR
      --     tao
      --     deltaT
      --     stdR2
      -- pointsGenerator =
      --   generatePath' thetaSigma scaleDist poissonLambda periodEnv (tao / deltaT) deltaT
      histFuncSingleGPU =
        computeFourierCoefficientsGPU
          phiFreqs
          rhoFreqs
          thetaFreqs
          rFreqs
          (use freqArr)
          (constant . double2Float $ sigma)
          (constant . double2Float $ log periodEnv)
      monterCarloHistFuncGPU =
        computeHistogramFromMonteCarloParallelMultipleGPU
          filePath
          pointsGenerator
          histFuncSingleGPU
          addHistogramUnsafe
          ptxs
  runBatch numGens numTrails batchSize monterCarloHistFuncGPU initHist

{-# INLINE countR2S1 #-}
countR2S1 ::
     GenIO
  -> Double
  -> (Int, Int)
  -> (Int, Int)
  -> Int
  -> [DList Particle]
  -> IO (Histogram Double)
countR2S1 randomGen thetaSigma xRange@(!xMin, !xMax) yRange@(!yMin, !yMax) !numOrientations !xs = do
  let !deltaTheta = 2 * pi / (fromIntegral numOrientations)
  ys <-
    fmap
      (L.filter
         (\((_, x, y), _) ->
            (inRange xRange x) && (inRange yRange y) && (x /= 0 || y /= 0)) .
       L.map
         (\(Particle phi rho theta' _ _) ->
            let !x = rho * cos phi
                !y = rho * sin phi
                !theta =
                  if theta' < 0
                    then theta' + 2 * pi
                    else theta'
             in ((floor $ theta / deltaTheta, round x, round y), 1)) .
       L.concat) .
    M.mapM
      (\particle@(Particle phi rho theta r _) -> do
         let delta = 1
             n = Prelude.floor $ r / delta
         zs <-
           M.mapM
             (\i -> do
                let thetaDist =
                      normalDistr
                        0
                        (thetaSigma * sqrt (delta * fromIntegral i / r))
                    (Particle a b c d _) =
                      FokkerPlanck.BrownianMotion.moveParticle
                        1
                        (Particle
                           phi
                           rho
                           theta
                           (Prelude.fromIntegral i * delta)
                           1)
                deltaThetaDiffusion <- genContVar thetaDist randomGen
                return (Particle a b (c `thetaPlus` deltaThetaDiffusion) d 1))
             [1 .. n - 1]
         return .
           L.concatMap
             (\(Particle phi' rho' theta' r' v) ->
                [ (Particle phi' rho' theta' r' v)
                , (Particle (-phi') rho' (-theta') r' v)
                ]) $
           (FokkerPlanck.BrownianMotion.moveParticle
              1
              (Particle phi rho theta r 1)) :
           zs) .
    DL.toList . DL.concat $
    xs
  return .
    Histogram
      [(yMax - yMin + 1), (xMax - xMin + 1), numOrientations]
      (L.length ys) .
    toUnboxedVector .
    UA.accum (+) 0 ((0, xMin, yMin), (numOrientations - 1, xMax, yMax)) $
    ys

{-# INLINE solveMonteCarloR2S1 #-}
solveMonteCarloR2S1 ::
     Int
  -> Int
  -> Int
  -> Int
  -> Int
  -> Int
  -> Double
  -> Double
  -> Double
  -> Double
  -> FilePath
  -> IO (R.Array U DIM3 Double)
solveMonteCarloR2S1 numGens numTrails batchSize xLen yLen numOrientations thetaSigma tao r initSpeed histFilePath = do
  gen <- createSystemRandom
  let !xShift = div xLen 2
      xRange' =
        if odd xLen
          then (-xShift, xShift)
          else (-xShift, xShift - 1)
      !yShift = div yLen 2
      yRange' =
        if odd yLen
          then (-yShift, yShift)
          else (-yShift, yShift - 1)
      thetaDist = normalDistrE 0 thetaSigma
      scaleDist = normalDistrE 0 0
      poissonDist = poissonE 0
      pointsGenerator =
        generatePath
          thetaDist
          scaleDist
          poissonDist
          (r^2)
          r -- (sqrt . fromIntegral $ xLen ^ 2 + yLen ^ 2)
          tao
          initSpeed
          1
      histFunc = countR2S1 gen thetaSigma xRange' yRange' numOrientations
      monterCarloHistFunc =
        computeHistogramFromMonteCarloParallel
          histFilePath
          pointsGenerator
          histFunc
          addHistogramUnsafe
  hist <-
    runBatch
      numGens
      numTrails
      batchSize
      monterCarloHistFunc
      (emptyHistogram [yLen, xLen, numOrientations] 0)
  return . getNormalizedHistogramArr $ hist

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 DataKinds #-}
{-# LANGUAGE NoMonomorphismRestriction #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE MultiParamTypeClasses #-}
module Engine.Matrix.NewMatrix where

import Foreign.Storable (Storable(..))
import qualified Data.Vector.Storable as V

import GHC.TypeLits
import Numeric.LinearAlgebra.Static

newtype Matrix w h t = Matrix (L w h)
  deriving (Num, Fractional, Floating)

newtype Vector l t = Vector (R l)
  deriving (Num, Fractional, Floating)

instance (Sized t (R l) V.Vector) => Sized t (Vector l t) V.Vector where
    konst = Vector . konst
    unwrap (Vector r) = unwrap r

(<&>) :: (KnownNat w1, KnownNat w2, KnownNat h) =>
    Matrix w1 h t -> Matrix w2 h t -> Matrix (w1+w2) h t
(<&>) (Matrix l1) (Matrix l2) = Matrix $ l1 —— l2

max_stars_repo_name: jaredloomis/Haskell-OpenGL


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

{-|
Module:      TextShow.Data.Ratio
Copyright:   (C) 2014-2017 Ryan Scott
License:     BSD-style (see the file LICENSE)
Maintainer:  Ryan Scott
Stability:   Provisional
Portability: GHC

'TextShow' instances for 'Ratio'.

/Since: 2/
-}
module TextShow.Data.Complex () where

import Data.Complex (Complex)

import TextShow.Classes (TextShow(..))
import TextShow.Data.Floating ()
import TextShow.TH.Internal (deriveTextShow1, makeShowbPrec)

-- | /Since: 2/
instance TextShow a => TextShow (Complex a) where
    {-# SPECIALIZE instance TextShow (Complex Float)  #-}
    {-# SPECIALIZE instance TextShow (Complex Double) #-}
    showbPrec = $(makeShowbPrec ''Complex)
    {-# INLINE showbPrec #-}

-- | /Since: 2/
$(deriveTextShow1 ''Complex)

max_stars_repo_name: janakg/duckling

max_issues_repo_name: janakg/duckling

max_forks_repo_name: janakg/duckling


text:
import System.Environment
import Data.Complex
import Paraiso

main = do
  args <- getArgs
  let arch = if "--cuda" `elem` args then
               CUDA 128 128
             else
               X86
  putStrLn $ compile arch code
    where
      code = do
         parallel 16384 $ do
           r <- allocate
           x <- allocate
           r =$ Rand 0.0 (4.0::Double)
           x =$ Rand 0.0 (1.0::Double)
           cuda $ do
             sequential 65536 $ do
               x =$ r * x * (1-x)
           output [r,x]


text:
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE DeriveGeneric #-}

module Lib
    ( genImage
    , ViewPort(..)
    ) where

import Codec.Picture
import Codec.Picture.Types
import Data.Fixed (mod')
import Data.Word (Word8)
import Data.Array
import Data.Hashable
import Data.Complex
import Control.Parallel.Strategies
import GHC.Generics

data ViewPort = ViewPort { _centerX :: Double
                         , _centerY :: Double
                         , _magnification :: Double
                         , _maxIterations :: Int
                         , _rotation :: Double
                         , _superScaling :: Int
                         } deriving (Eq, Show, Generic)

instance Hashable ViewPort

genImage :: Int -> Int -> ViewPort -> DynamicImage
genImage w h v = ImageRGB8 $ generateImage (renderImage v w h) w h

renderImage :: ViewPort -> Int -> Int -> (Int -> Int -> PixelRGB8)
renderImage (ViewPort cx cy mg mi _ s) w h =
    let aspectRatio = fromIntegral w / fromIntegral h :: Double
        (rw, rh) = if aspectRatio <= 4.0 / 3.0 then (4.0 / mg, 4.0 / mg / aspectRatio)
                                               else (3.0 / mg * aspectRatio, 3.0 / mg)
        f x y = let ps = [(x * s + ix, y * s + iy) | ix <- [0..s-1], iy <- [0..s-1]]
                    r x' = cx - (rw / 2.0) + ((fromIntegral x' + 0.5) / fromIntegral (w * s) :: Double) * rw
                    i y' = cy - (rh / 2.0) + ((fromIntegral y' + 0.5) / fromIntegral (h * s) :: Double) * rh
                    rs = fmap (\(x', y') -> (r x', i y')) ps
                    qs = fmap (\(r0, i0) -> color (fastRender (realToFrac r0 :+ realToFrac i0) mi)) rs
                    (rc, gc, bc) = foldr (\(PixelRGB8 r g b) (rc :: Int, gc :: Int, bc :: Int) -> (rc + fromIntegral r, gc + fromIntegral g, bc + fromIntegral b)) (0, 0, 0) qs
                 in PixelRGB8 (fromIntegral (rc `div` (s*s))) (fromIntegral (gc `div` (s*s))) (fromIntegral (bc `div` (s*s)))
    in f

mag2 :: (RealFloat a, Num a) => Complex a -> a
mag2 (r :+ i) = r*r + i*i

render :: (RealFloat a, Num a) => Complex a -> Int -> Maybe (Complex a, Int)
render c maxI = loop (0 :+ 0) 0
    where loop z i = if mag2 z <= 4 && i <= maxI
                        then loop (step z) (i + 1)
                        else if i < maxI
                            then Just (step $ step z, i)
                            else Nothing
          step z = z*z + c

fastRender :: (RealFloat a, Num a) => Complex a -> Int -> Maybe (Complex a, Int)
fastRender c maxI = fastLoop (0 :+ 0) 0
    where fastLoop z i = let zn = (step . step . step . step . step . step . step . step) z
                         in if mag2 zn <= 4 && (i + 8) <= maxI
                            then fastLoop zn (i + 8)
                            else slowLoop z i
          slowLoop z i = if mag2 z <= 4 && i <= maxI
                         then slowLoop (step z) (i + 1)
                         else if i < maxI
                              then Just (step $ step z, i)
                              else Nothing
          step z = z*z + c

color :: (RealFloat a, Num a) => Maybe (Complex a, Int) -> PixelRGB8
color Nothing = PixelRGB8 0 0 0
color (Just (z, i)) = let mu = ((fromIntegral (i + 1)) - (logBase 2 . log . mag2) z) / 80.0
                          frac = snd (properFraction mu)
                      in getColor $ if frac < 0 then frac + 1.0 else frac
{-
renderPoint :: Double     -- ^ The real coordinate
            -> Double     -- ^ The imaginary coordinate
            -> Int        -- ^ The max iterations
            -> PixelRGB8  -- ^ The resulting pixel
renderPoint r0 i0 mi = case loop (0, 0) 0 of
                       Left (p, iter) -> let (r2, i2) = step . step $ p
                                             mu = ((fromIntegral (iter + 1)) - (logBase 2 . log) (r2 * r2 + i2 * i2)) / 80.0
                                             frac = snd (properFraction mu)
                                         in getColor $ if frac < 0 then frac + 1.0 else frac
                       Right _ -> PixelRGB8 0 0 0
    where loop (r, i) iter = if r*r + i*i <= 4 && iter <= mi
                          then loop (step (r, i)) (iter + 1)
                          else if iter < mi
                               then Left ((r, i), iter)
                               else Right ()
          step (r, i) = (r*r - i*i + r0, 2*r*i + i0)
-}

palette :: (RealFloat a, Num a) => [(a, PixelRGB8)]
palette = [ (0.0,    PixelRGB8 0   7   100)
          , (0.16,   PixelRGB8 32  107 203)
          , (0.42,   PixelRGB8 237 255 255)
          , (0.6425, PixelRGB8 255 170 0)
          , (0.8575, PixelRGB8 0   2   0)
          , (1.0,    PixelRGB8 0   7   100)
          ]

getColor :: (RealFloat a, Num a) => a -> PixelRGB8
getColor f = let ((li, (PixelRGB8 lr lg lb)), (ri, (PixelRGB8 rr rg rb))) = pair palette
                 index = (f - li) / (ri - li)
                 index2 = index * index
                 index3 = index * index2
                 r = fromIntegral . floor $ (2.0 * (fromIntegral lr - fromIntegral rr) * index3) - (3.0 * (fromIntegral lr - fromIntegral rr) * index2) + fromIntegral lr
                 g = fromIntegral . floor $ (2.0 * (fromIntegral lg - fromIntegral rg) * index3) - (3.0 * (fromIntegral lg - fromIntegral rg) * index2) + fromIntegral lg
                 b = fromIntegral . floor $ (2.0 * (fromIntegral lb - fromIntegral rb) * index3) - (3.0 * (fromIntegral lb - fromIntegral rb) * index2) + fromIntegral lb
             in PixelRGB8 r g b
    where pair (c:cs) = if f >= fst (head cs) then pair cs else (c, head cs)

max_stars_repo_name: micahhahn/DeepBlue

max_issues_repo_name: micahhahn/DeepBlue

max_forks_repo_name: micahhahn/DeepBlue


text:
{-# LANGUAGE GADTs #-}

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

module Feldspar.Core.Middleend.FromTypeUtil where

import Data.Complex

import Feldspar.Core.Types
import Feldspar.Core.UntypedRepresentation hiding ( Lambda, UntypedFeldF(..)
                                                  , Size, ScalarType(..)
                                                  , Type(..), Signedness
                                                  , Op(..)
                                                  )
import qualified Feldspar.Core.UntypedRepresentation as Ut
import Feldspar.ValueInfo
import Feldspar.Range

untypeType :: TypeRep a -> Size a -> Ut.Type
untypeType UnitType _               = Ut.TupType []
untypeType BoolType _               = 1 Ut.:# Ut.BoolType
untypeType (IntType s n) _          = 1 Ut.:# Ut.IntType (convSign s) (convSize n)
untypeType FloatType _              = 1 Ut.:# Ut.FloatType
untypeType DoubleType _             = 1 Ut.:# Ut.DoubleType
untypeType (ComplexType t) _        = 1 Ut.:# Ut.ComplexType (untypeType t (defaultSize t))
untypeType (Tup2Type a b) (sa,sb)
  = Ut.TupType [untypeType a sa, untypeType b sb]
untypeType (Tup3Type a b c) (sa,sb,sc)
  = Ut.TupType [untypeType a sa, untypeType b sb, untypeType c sc]
untypeType (Tup4Type a b c d) (sa,sb,sc,sd)
  = Ut.TupType [ untypeType a sa, untypeType b sb, untypeType c sc
               , untypeType d sd
               ]
untypeType (Tup5Type a b c d e) (sa,sb,sc,sd,se)
  = Ut.TupType [ untypeType a sa, untypeType b sb, untypeType c sc
               , untypeType d sd, untypeType e se
               ]
untypeType (Tup6Type a b c d e f) (sa,sb,sc,sd,se,sf)
  = Ut.TupType [ untypeType a sa, untypeType b sb, untypeType c sc
               , untypeType d sd, untypeType e se, untypeType f sf
               ]
untypeType (Tup7Type a b c d e f g) (sa,sb,sc,sd,se,sf,sg)
  = Ut.TupType [ untypeType a sa, untypeType b sb, untypeType c sc
               , untypeType d sd, untypeType e se, untypeType f sf
               , untypeType g sg
               ]
untypeType (Tup8Type a b c d e f g h) (sa,sb,sc,sd,se,sf,sg,sh)
  = Ut.TupType [ untypeType a sa, untypeType b sb, untypeType c sc
               , untypeType d sd, untypeType e se, untypeType f sf
               , untypeType g sg, untypeType h sh
               ]
untypeType (Tup9Type a b c d e f g h i) (sa,sb,sc,sd,se,sf,sg,sh,si)
  = Ut.TupType [ untypeType a sa, untypeType b sb, untypeType c sc
               , untypeType d sd, untypeType e se, untypeType f sf
               , untypeType g sg, untypeType h sh, untypeType i si
               ]
untypeType (Tup10Type a b c d e f g h i j) (sa,sb,sc,sd,se,sf,sg,sh,si,sj)
  = Ut.TupType [ untypeType a sa, untypeType b sb, untypeType c sc
               , untypeType d sd, untypeType e se, untypeType f sf
               , untypeType g sg, untypeType h sh, untypeType i si
               , untypeType j sj
               ]
untypeType (Tup11Type a b c d e f g h i j k) (sa,sb,sc,sd,se,sf,sg,sh,si,sj,sk)
  = Ut.TupType [ untypeType a sa, untypeType b sb, untypeType c sc
               , untypeType d sd, untypeType e se, untypeType f sf
               , untypeType g sg, untypeType h sh, untypeType i si
               , untypeType j sj, untypeType k sk
               ]
untypeType (Tup12Type a b c d e f g h i j k l) (sa,sb,sc,sd,se,sf,sg,sh,si,sj,sk,sl)
  = Ut.TupType [ untypeType a sa, untypeType b sb, untypeType c sc
               , untypeType d sd, untypeType e se, untypeType f sf
               , untypeType g sg, untypeType h sh, untypeType i si
               , untypeType j sj, untypeType k sk, untypeType l sl
               ]
untypeType (Tup13Type a b c d e f g h i j k l m) (sa,sb,sc,sd,se,sf,sg,sh,si,sj,sk,sl,sm)
  = Ut.TupType [ untypeType a sa, untypeType b sb, untypeType c sc
               , untypeType d sd, untypeType e se, untypeType f sf
               , untypeType g sg, untypeType h sh, untypeType i si
               , untypeType j sj, untypeType k sk, untypeType l sl
               , untypeType m sm
               ]
untypeType (Tup14Type a b c d e f g h i j k l m n) (sa,sb,sc,sd,se,sf,sg,sh,si,sj,sk,sl,sm,sn)
  = Ut.TupType [ untypeType a sa, untypeType b sb, untypeType c sc
               , untypeType d sd, untypeType e se, untypeType f sf
               , untypeType g sg, untypeType h sh, untypeType i si
               , untypeType j sj, untypeType k sk, untypeType l sl
               , untypeType m sm, untypeType n sn
               ]
untypeType (Tup15Type a b c d e f g h i j k l m n o) (sa,sb,sc,sd,se,sf,sg,sh,si,sj,sk,sl,sm,sn,so)
  = Ut.TupType [ untypeType a sa, untypeType b sb, untypeType c sc
               , untypeType d sd, untypeType e se, untypeType f sf
               , untypeType g sg, untypeType h sh, untypeType i si
               , untypeType j sj, untypeType k sk, untypeType l sl
               , untypeType m sm, untypeType n sn, untypeType o so
               ]
untypeType (MutType a) sz           = Ut.MutType (untypeType a sz)
untypeType (RefType a) sz           = Ut.RefType (untypeType a sz)
untypeType (ArrayType a) (rs :> es) = Ut.ArrayType rs (untypeType a es)
untypeType (MArrType a) (rs :> es)  = Ut.MArrType rs (untypeType a es)
untypeType (ParType a) sz           = Ut.ParType (untypeType a sz)
untypeType (ElementsType a) (rs :> es) = Ut.ElementsType (untypeType a es)
untypeType (ConsType a b) (sa,sb)   = Ut.TupType $ untypeTup (ConsType a b) (sa,sb)
untypeType NilType        _         = Ut.TupType []
untypeType (TupleType t) sz         = untypeType t sz
untypeType (IVarType a) sz          = Ut.IVarType $ untypeType a sz
untypeType (FunType a b) (sa, sz)   = Ut.FunType (untypeType a sa) (untypeType b sz)
untypeType (FValType a) sz          = Ut.FValType (untypeType a sz)
untypeType typ _                    = error "untypeType: missing "

untypeTup :: TypeRep (RTuple a) -> Size (RTuple a) -> [Ut.Type]
untypeTup (ConsType a b) (sa,sb) = untypeType a sa : untypeTup b sb
untypeTup NilType        _       = []

convSign :: Signedness a -> Ut.Signedness
convSign U       = Unsigned
convSign S       = Signed

convSize :: BitWidth a -> Ut.Size
convSize N8      = S8
convSize N16     = S16
convSize N32     = S32
convSize N64     = S64
convSize NNative = S32

literal :: TypeRep a -> Size a -> a -> Lit
literal t@UnitType        sz a = literalConst t sz a
literal t@BoolType        sz a = literalConst t sz a
literal t@IntType{}       sz a = literalConst t sz a
literal t@FloatType       sz a = literalConst t sz a
literal t@DoubleType      sz a = literalConst t sz a
literal t@ComplexType{}   sz a = literalConst t sz a
literal t@ArrayType{}     sz a = literalConst t sz a
literal (Tup2Type ta tb) (sa,sb) (a,b)
    = LTup [literal ta sa a, literal tb sb b]

literal (Tup3Type ta tb tc) (sa,sb,sc) (a,b,c)
    = LTup [literal ta sa a, literal tb sb b, literal tc sc c]

literal (Tup4Type ta tb tc td) (sa,sb,sc,sd) (a,b,c,d)
    = LTup [ literal ta sa a, literal tb sb b, literal tc sc c
           , literal td sd d
           ]

literal (Tup5Type ta tb tc td te) (sa,sb,sc,sd,se) (a,b,c,d,e)
    = LTup [ literal ta sa a, literal tb sb b, literal tc sc c
           , literal td sd d, literal te se e
           ]

literal (Tup6Type ta tb tc td te tf) (sa,sb,sc,sd,se,sf) (a,b,c,d,e,f)
    = LTup [ literal ta sa a, literal tb sb b, literal tc sc c
           , literal td sd d, literal te se e, literal tf sf f
           ]

literal (Tup7Type ta tb tc td te tf tg) (sa,sb,sc,sd,se,sf,sg) (a,b,c,d,e,f,g)
    = LTup [ literal ta sa a, literal tb sb b, literal tc sc c
           , literal td sd d, literal te se e, literal tf sf f
           , literal tg sg g
           ]
literal t s a = error "Missing pattern: FromTyped.hs: literal"

literalConst :: TypeRep a -> Size a -> a -> Lit
literalConst UnitType        _  ()     = LTup []
literalConst BoolType        _  a      = LBool a
literalConst (IntType s n)   sz a      = LInt (convSign s) (convSize n) (toInteger a)
literalConst FloatType       _  a      = LFloat a
literalConst DoubleType      _  a      = LDouble a
literalConst (ArrayType t)   _  a      = LArray t' $ map (literalConst t (defaultSize t)) a
  where t' = untypeType t (defaultSize t)
literalConst (ComplexType t) _  (r:+i) = LComplex re ie
  where re = literalConst t (defaultSize t) r
        ie = literalConst t (defaultSize t) i

-- | Construct a ValueInfo from a TypeRep and a Size
toValueInfo :: TypeRep a -> Size a -> ValueInfo
toValueInfo UnitType          _             = VIProd []
-- FIXME: No range for boolean types yet.
toValueInfo BoolType          _             = VIBool    (Range 0 1)
toValueInfo (IntType U N8)    r             = VIWord8   r
toValueInfo (IntType S N8)    r             = VIInt8    r
toValueInfo (IntType U N16)   r             = VIWord16  r
toValueInfo (IntType S N16)   r             = VIInt16   r
toValueInfo (IntType U N32)   r             = VIWord32  r
toValueInfo (IntType S N32)   r             = VIInt32   r
toValueInfo (IntType U N64)   r             = VIWord64  r
toValueInfo (IntType S N64)   r             = VIInt64   r
toValueInfo (IntType U NNative) r           = VIWordN   r
toValueInfo (IntType S NNative) r           = VIIntN    r
-- FIXME: No range for FP types and ComplexType yet.
toValueInfo FloatType         _             = VIFloat
toValueInfo DoubleType        _             = VIDouble
toValueInfo (ComplexType _)   _             = VIProd []
toValueInfo (Tup2Type a b) (sa,sb)
  = VIProd [toValueInfo a sa, toValueInfo b sb]
toValueInfo (Tup3Type a b c) (sa,sb,sc)
  = VIProd [toValueInfo a sa, toValueInfo b sb, toValueInfo c sc]
toValueInfo (Tup4Type a b c d) (sa,sb,sc,sd)
  = VIProd [ toValueInfo a sa, toValueInfo b sb, toValueInfo c sc
           , toValueInfo d sd
           ]
toValueInfo (Tup5Type a b c d e) (sa,sb,sc,sd,se)
  = VIProd [ toValueInfo a sa, toValueInfo b sb, toValueInfo c sc
           , toValueInfo d sd, toValueInfo e se
           ]
toValueInfo (Tup6Type a b c d e f) (sa,sb,sc,sd,se,sf)
  = VIProd [ toValueInfo a sa, toValueInfo b sb, toValueInfo c sc
           , toValueInfo d sd, toValueInfo e se, toValueInfo f sf
           ]
toValueInfo (Tup7Type a b c d e f g) (sa,sb,sc,sd,se,sf,sg)
  = VIProd [ toValueInfo a sa, toValueInfo b sb, toValueInfo c sc
           , toValueInfo d sd, toValueInfo e se, toValueInfo f sf
           , toValueInfo g sg
           ]
toValueInfo (Tup8Type a b c d e f g h) (sa,sb,sc,sd,se,sf,sg,sh)
  = VIProd [ toValueInfo a sa, toValueInfo b sb, toValueInfo c sc
           , toValueInfo d sd, toValueInfo e se, toValueInfo f sf
           , toValueInfo g sg, toValueInfo h sh
           ]
toValueInfo (Tup9Type a b c d e f g h i) (sa,sb,sc,sd,se,sf,sg,sh,si)
  = VIProd [ toValueInfo a sa, toValueInfo b sb, toValueInfo c sc
           , toValueInfo d sd, toValueInfo e se, toValueInfo f sf
           , toValueInfo g sg, toValueInfo h sh, toValueInfo i si
            ]
toValueInfo (Tup10Type a b c d e f g h i j) (sa,sb,sc,sd,se,sf,sg,sh,si,sj)
  = VIProd [ toValueInfo a sa, toValueInfo b sb, toValueInfo c sc
           , toValueInfo d sd, toValueInfo e se, toValueInfo f sf
           , toValueInfo g sg, toValueInfo h sh, toValueInfo i si
           , toValueInfo j sj
           ]
toValueInfo (Tup11Type a b c d e f g h i j k) (sa,sb,sc,sd,se,sf,sg,sh,si,sj,sk)
  = VIProd [ toValueInfo a sa, toValueInfo b sb, toValueInfo c sc
           , toValueInfo d sd, toValueInfo e se, toValueInfo f sf
           , toValueInfo g sg, toValueInfo h sh, toValueInfo i si
           , toValueInfo j sj, toValueInfo k sk
           ]
toValueInfo (Tup12Type a b c d e f g h i j k l) (sa,sb,sc,sd,se,sf,sg,sh,si,sj,sk,sl)
  = VIProd [ toValueInfo a sa, toValueInfo b sb, toValueInfo c sc
           , toValueInfo d sd, toValueInfo e se, toValueInfo f sf
           , toValueInfo g sg, toValueInfo h sh, toValueInfo i si
           , toValueInfo j sj, toValueInfo k sk, toValueInfo l sl
           ]
toValueInfo (Tup13Type a b c d e f g h i j k l m) (sa,sb,sc,sd,se,sf,sg,sh,si,sj,sk,sl,sm)
  = VIProd [ toValueInfo a sa, toValueInfo b sb, toValueInfo c sc
           , toValueInfo d sd, toValueInfo e se, toValueInfo f sf
           , toValueInfo g sg, toValueInfo h sh, toValueInfo i si
           , toValueInfo j sj, toValueInfo k sk, toValueInfo l sl
           , toValueInfo m sm
           ]
toValueInfo (Tup14Type a b c d e f g h i j k l m n) (sa,sb,sc,sd,se,sf,sg,sh,si,sj,sk,sl,sm,sn)
  = VIProd [ toValueInfo a sa, toValueInfo b sb, toValueInfo c sc
           , toValueInfo d sd, toValueInfo e se, toValueInfo f sf
           , toValueInfo g sg, toValueInfo h sh, toValueInfo i si
           , toValueInfo j sj, toValueInfo k sk, toValueInfo l sl
           , toValueInfo m sm, toValueInfo n sn
           ]
toValueInfo (Tup15Type a b c d e f g h i j k l m n o) (sa,sb,sc,sd,se,sf,sg,sh,si,sj,sk,sl,sm,sn,so)
  = VIProd [ toValueInfo a sa, toValueInfo b sb, toValueInfo c sc
           , toValueInfo d sd, toValueInfo e se, toValueInfo f sf
           , toValueInfo g sg, toValueInfo h sh, toValueInfo i si
           , toValueInfo j sj, toValueInfo k sk, toValueInfo l sl
           , toValueInfo m sm, toValueInfo n sn, toValueInfo o so
           ]
toValueInfo (MutType a) sz                  = toValueInfo a sz
toValueInfo (RefType a) sz                  = toValueInfo a sz
toValueInfo (ArrayType a) (Range (WordN l) (WordN r) :> es)
  = VIProd [VIWord32 (Range l r), toValueInfo a es]
toValueInfo (MArrType a) (Range (WordN l) (WordN r) :> es)
  = VIProd [VIWord32 (Range l r), toValueInfo a es]
toValueInfo (ParType a) sz                  = toValueInfo a sz
toValueInfo (ElementsType a) (Range (WordN l) (WordN r) :> es)
  = VIProd [VIWord32 (Range l r), toValueInfo a es]
toValueInfo (ConsType a b) (sa,sb) = VIProd $ toValueInfo a sa : ss
  where VIProd ss = toValueInfo b sb
toValueInfo NilType _ = VIProd []
toValueInfo (TupleType t) sz = toValueInfo t sz
toValueInfo (IVarType a) sz                 = toValueInfo a sz
-- TODO: Maybe keep argument information for FunType.
toValueInfo (FunType a b) (sa, _)           = toValueInfo a sa
toValueInfo (FValType a) sz                 = toValueInfo a sz
toValueInfo typ _
  = error $ "toValueInfo: missing case for " ++ show typ

max_stars_repo_name: pjonsson/feldspar-language

max_issues_repo_name: pjonsson/feldspar-language

max_forks_repo_name: pjonsson/feldspar-language