text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
-- This module uses the test-framework-quickcheck2 package.
module Main where

import Control.Monad
import qualified Data.Vector.Unboxed as V
import qualified Data.Vector.Storable as VS
import Data.Complex

import Test.Framework (defaultMain, testGroup)
import Test.Framework.Providers.QuickCheck2 (testProperty)
import Test.QuickCheck

import qualified Numeric.FFT.Vector.Invertible as I
import qualified Numeric.FFT.Vector.Invertible.Multi as IM
import qualified Numeric.FFT.Vector.Unitary as U
import qualified Numeric.FFT.Vector.Unitary.Multi as UM
import Numeric.FFT.Vector.Plan

main = defaultMain
            -- NB: There's no explicit tests for the Unnormalized package.
            -- However, its Planners are implicitly used by the other modules,
            -- so it's covered in the below tests.
            [ testGroup "invertibility"
              [ testProperty "I.dft" $ prop_invert I.dft I.idft
              , testProperty "I.dftR2C" $ prop_invert I.dftR2C I.dftC2R
              , testProperty "I.dct1" $ prop_invert I.dct1 I.idct1
              , testProperty "I.dct2" $ prop_invert I.dct2 I.idct2
              , testProperty "I.dct3" $ prop_invert I.dct3 I.idct3
              , testProperty "I.dct4" $ prop_invert I.dct4 I.idct4
              , testProperty "I.dst1" $ prop_invert I.dst1 I.idst1
              , testProperty "I.dst2" $ prop_invert I.dst2 I.idst2
              , testProperty "I.dst3" $ prop_invert I.dst3 I.idst3
              , testProperty "I.dst4" $ prop_invert I.dst4 I.idst4
              , testProperty "U.dft" $ prop_invert U.dft U.idft
              , testProperty "U.dftR2C" $ prop_invert U.dftR2C U.dftC2R
              , testProperty "U.dct2" $ prop_invert U.dct2 U.idct2
              ]
            , testGroup "orthogonality"
              [ testProperty "U.dft" $ prop_orthog U.dft
              , testProperty "U.idft" $ prop_orthog U.idft
              , testProperty "U.dftR2C" $ prop_orthog U.dftR2C
              , testProperty "U.dftC2R" $ prop_orthog U.dftR2C
              , testProperty "U.dct2" $ prop_orthog U.dct2
              , testProperty "U.idct2" $ prop_orthog U.idct2
              , testProperty "U.dct4" $ prop_orthog U.dct4
              ]
            , testGroup "invertibility ND"
              [ testProperty "IM.dft" $ prop_invertND IM.dft IM.idft
              , testProperty "IM.dftR2C" $ prop_invertND IM.dftR2C IM.dftC2R
              , testProperty "UM.dft" $ prop_invertND UM.dft UM.idft
              , testProperty "UM.dftR2C" $ prop_invertND UM.dftR2C UM.dftC2R
              ]
            , testGroup "orthogonality"
              [ testProperty "UM.dft" $ prop_orthogND UM.dft
              , testProperty "UM.idft" $ prop_orthogND UM.idft
              , testProperty "UM.dftR2C" $ prop_orthogND UM.dftR2C
              , testProperty "UM.dftC2R" $ prop_orthogND UM.dftR2C
              ]
            ]

-------------------
-- An instance of Arbitrary that probably belongs in another package.

instance (V.Unbox a, Arbitrary a) => Arbitrary (V.Vector a) where
    arbitrary = V.fromList `fmap` arbitrary

-------------------------
-- Support functions to compare Doubles for (near) equality.

class Num a => Mag a where
    mag :: a -> Double

instance Mag Double where
    mag = abs

instance Mag (Complex Double) where
    mag = magnitude

-- Robustly test whether two Doubles are nearly identical.
close :: Mag a => a -> a -> Bool
close x y = tol > mag (x-y) / max 1 (mag x + mag y)
  where
    tol = 1e-10

withinTol :: (Mag a, V.Unbox a) => V.Vector a -> V.Vector a -> Bool
withinTol a b
    | V.length a /= V.length b = False
    | otherwise = V.and $ V.zipWith close a b

---------------------
-- The actual properties

-- Test whether the inverse actually inverts the forward transform.
prop_invert f g a = let
                        p1 = plan f (V.length a)
                        p2 = plan g (V.length a)
                    in (V.length a > 1) ==> withinTol a $ execute p2 $ execute p1 a

-- Test whether the transform preserves the L2 (sum-of-squares) norm.
prop_orthog f a = let
                    p1 = plan f (V.length a)
                  in (V.length a > 1) ==> close (norm2 a) (norm2 $ execute p1 a)

data DimsAndValues a = DimsAndValues (VS.Vector Int) (V.Vector a)
  deriving (Show)

instance (Arbitrary a, V.Unbox a) => Arbitrary (DimsAndValues a) where
  arbitrary = do
    dims <- liftM (VS.fromList . map getPositive) arbitrary `suchThatMap` maybeReduceSize
    values <- V.replicateM (VS.product dims) arbitrary
    return (DimsAndValues dims values)
    where
      -- We use this to prevent test cases from growing too big
      maybeReduceSize ds =
        if VS.product ds < 1000 then Just ds else maybeReduceSize (VS.init ds)

prop_invertND f g (DimsAndValues ds a) = let
                        p1 = planND f ds
                        p2 = planND g ds
                    in (V.length a > 1) ==> withinTol a $ execute p2 $ execute p1 a

prop_orthogND f (DimsAndValues ds a) = let
                    p1 = planND f ds
                  in (V.length a > 1) ==> close (norm2 a) (norm2 $ execute p1 a)

norm2 a = sqrt $ V.sum $ V.map (\x -> x*x) $ V.map mag a

max_stars_repo_name: TravisWhitaker/vector-fftw

max_issues_repo_name: TravisWhitaker/vector-fftw

max_forks_repo_name: TravisWhitaker/vector-fftw


text:
{-# LANGUAGE ConstraintKinds     #-}
{-# LANGUAGE DataKinds           #-}

{-# LANGUAGE DeriveTraversable   #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE RankNTypes          #-}
{-# LANGUAGE ScopedTypeVariables #-}

{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
module Test.Grenade.Recurrent.Layers.LSTM.Reference where

import           Data.Reflection
import           GHC.TypeLits                  (KnownNat)
import           Numeric.AD.Internal.Reverse   (Tape)
import           Numeric.AD.Mode.Reverse
import qualified Numeric.LinearAlgebra         as H
import qualified Numeric.LinearAlgebra.Static  as S

import           Grenade.Recurrent.Layers.LSTM (LSTMWeights (..))
import qualified Grenade.Recurrent.Layers.LSTM as LSTM
import           Grenade.Types

--
-- This module contains a set of list only versions of
-- an LSTM layer which can be used with the AD library.
--
-- Using this, we can check to make sure that our fast
-- back propagation implementation is correct.
--

-- | List only matrix deriving functor
newtype Matrix a = Matrix {
    matrixWeights         :: [[a]]
  } deriving (Functor, Foldable, Traversable, Eq, Show)

-- | List only vector deriving functor
newtype Vector a = Vector {
    vectorWeights         :: [a]
  } deriving (Functor, Foldable, Traversable, Eq, Show)

-- | List only LSTM weights
data RefLSTM a = RefLSTM
    { refLstmWf :: Matrix a -- Weight Forget     (W_f)
    , refLstmUf :: Matrix a -- Cell State Forget (U_f)
    , refLstmBf :: Vector a -- Bias Forget       (b_f)
    , refLstmWi :: Matrix a -- Weight Input      (W_i)
    , refLstmUi :: Matrix a -- Cell State Input  (U_i)
    , refLstmBi :: Vector a -- Bias Input        (b_i)
    , refLstmWo :: Matrix a -- Weight Output     (W_o)
    , refLstmUo :: Matrix a -- Cell State Output (U_o)
    , refLstmBo :: Vector a -- Bias Output       (b_o)
    , refLstmWc :: Matrix a -- Weight Cell       (W_c)
    , refLstmBc :: Vector a -- Bias Cell         (b_c)
    } deriving (Functor, Foldable, Traversable, Eq, Show)

lstmToReference :: (KnownNat a, KnownNat b) => LSTM.LSTMWeights a b -> RefLSTM RealNum
lstmToReference lw =
    RefLSTM
      { refLstmWf = Matrix . H.toLists . S.extract $ lstmWf lw -- Weight Forget     (W_f)
      , refLstmUf = Matrix . H.toLists . S.extract $ lstmUf lw -- Cell State Forget (U_f)
      , refLstmBf = Vector . H.toList  . S.extract $ lstmBf lw -- Bias Forget       (b_f)
      , refLstmWi = Matrix . H.toLists . S.extract $ lstmWi lw -- Weight Input      (W_i)
      , refLstmUi = Matrix . H.toLists . S.extract $ lstmUi lw -- Cell State Input  (U_i)
      , refLstmBi = Vector . H.toList  . S.extract $ lstmBi lw -- Bias Input        (b_i)
      , refLstmWo = Matrix . H.toLists . S.extract $ lstmWo lw -- Weight Output     (W_o)
      , refLstmUo = Matrix . H.toLists . S.extract $ lstmUo lw -- Cell State Output (U_o)
      , refLstmBo = Vector . H.toList  . S.extract $ lstmBo lw -- Bias Output       (b_o)
      , refLstmWc = Matrix . H.toLists . S.extract $ lstmWc lw -- Weight Cell       (W_c)
      , refLstmBc = Vector . H.toList  . S.extract $ lstmBc lw -- Bias Cell         (b_c)
      }

runLSTM :: Floating a => RefLSTM a -> Vector a -> Vector a -> (Vector a, Vector a)
runLSTM rl cell input =
    let -- Forget state vector
        f_t = sigmoid   $ refLstmBf rl #+ refLstmWf rl #> input #+ refLstmUf rl #> cell
        -- Input state vector
        i_t = sigmoid   $ refLstmBi rl #+ refLstmWi rl #> input #+ refLstmUi rl #> cell
        -- Output state vector
        o_t = sigmoid   $ refLstmBo rl #+ refLstmWo rl #> input #+ refLstmUo rl #> cell
        -- Cell input state vector
        c_x = fmap tanh $ refLstmBc rl #+ refLstmWc rl #> input
        -- Cell state
        c_t = f_t #* cell #+ i_t #* c_x
        -- Output (it's sometimes recommended to use tanh c_t)
        h_t = o_t #* c_t
    in (c_t, h_t)

runLSTMback :: forall a. Floating a => Vector a -> Vector a -> RefLSTM a -> RefLSTM a
runLSTMback cell input =
  grad f
    where
  f :: forall s. Reifies s Tape => RefLSTM (Reverse s a) -> Reverse s a
  f net =
    let cell'   = fmap auto cell
        input'  = fmap auto input
        (cells, forwarded) = runLSTM net cell' input'
    in  sum forwarded + sum cells

runLSTMbackOnInput :: forall a. Floating a => Vector a -> RefLSTM a -> Vector a -> Vector a
runLSTMbackOnInput cell net =
  grad f
    where
  f :: forall s. Reifies s Tape => Vector (Reverse s a) -> Reverse s a
  f input =
    let cell'   = fmap auto cell
        net'    = fmap auto net
        (cells, forwarded) = runLSTM net' cell' input
    in  sum forwarded + sum cells

runLSTMbackOnCell :: forall a. Floating a => Vector a -> RefLSTM a -> Vector a -> Vector a
runLSTMbackOnCell input net =
  grad f
    where
  f :: forall s. Reifies s Tape => Vector (Reverse s a) -> Reverse s a
  f cell =
    let input'  = fmap auto input
        net'    = fmap auto net
        (cells, forwarded) = runLSTM net' cell input'
    in  sum forwarded + sum cells

-- | Helper to multiply a matrix by a vector
matMult :: Num a => Matrix a -> Vector a -> Vector a
matMult (Matrix m) (Vector v) = Vector result
  where
    lrs = map length m
    l   = length v
    result = if all (== l) lrs
             then map (\r -> sum $ zipWith (*) r v) m
             else error $ "Matrix has rows of length " ++ show lrs ++
                          " but vector is of length " ++ show l

(#>) :: Num a => Matrix a -> Vector a -> Vector a
(#>) = matMult
infixr 8 #>

(#+) :: Num a => Vector a -> Vector a -> Vector a
(#+) (Vector as) (Vector bs) = Vector $ zipWith (+) as bs
infixl 6 #+

(#-) :: Num a => Vector a -> Vector a -> Vector a
(#-) (Vector as) (Vector bs) = Vector $ zipWith (-) as bs
infixl 6 #-

(#*) :: Num a => Vector a -> Vector a -> Vector a
(#*) (Vector as) (Vector bs) = Vector $ zipWith (*) as bs
infixl 7 #*

sigmoid :: (Functor f, Floating a) => f a -> f a
sigmoid xs = (\x -> 1 / (1 + exp (-x))) <$> xs

ext: hs

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE UndecidableInstances #-}
module LSI.FrequencyResponse where

import Data.Complex
import Data.Vector as V hiding (map, zip, foldr)
import Linear.V as V

import LSI.RationalFunction

computeFR :: RealFloat f
          => RationalFunction d f
          -> V d f
          -> Complex f
computeFR rf = (computeFR' rf) . V.toList .V.toVector
  where computeFR' f =
          case f of
            Monomial c es ->
              let fs = map toInteger . V.toList . V.toVector $ es in
                \vs -> {-# SCC monom #-}
                  let factors = map iexp $ zip zs fs
                      iexp (z, e) | e < 0 = 1/(z^(-e))
                      iexp (z, e) = z^e
                      zs = map (\w -> mkPolar 1.0 w) vs in
                    (c :+ 0.0) * (foldr (*) 1 factors)

            Add rf1 rf2 ->
              let fr1 = computeFR' rf1
                  fr2 = computeFR' rf2 in
                \vs -> (fr1 vs) + (fr2 vs)

            Mul rf1 rf2 ->
              let fr1 = computeFR' rf1
                  fr2 = computeFR' rf2 in
                \vs -> (fr1 vs) * (fr2 vs)

            Div rf1 rf2 ->
              let fr1 = computeFR' rf1
                  fr2 = computeFR' rf2 in
                \vs -> (fr1 vs) / (fr2 vs)

max_stars_repo_name: gdeest/lsi-systems

max_issues_repo_name: gdeest/lsi-systems

max_forks_repo_name: gdeest/lsi-systems


text:
import Test.Framework
import Test.Framework.Providers.QuickCheck2

import Vector
import STVector
import Matrix
import STMatrix
import Statistics

main :: IO ()
main = defaultMain tests
  where
    tests = [ tests_Vector
            , tests_STVector
            , tests_Matrix
            , tests_STMatrix
            , tests_Statistics
            ]

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
module Fractals (pMandelbrot, pJulia) where

import           Data.Complex

data PointResult = Convergent | Divergent Int (Complex Double)

-- Calculates a point in the Mandelbrot set
pMandelbrot :: Bool -> Double -> Int -> Double -> Complex Double -> Maybe Double
pMandelbrot aa pixelSize maxIterations power point =
   fPoint aa pixelSize maxIterations power (\d z -> z ** (d :+ 0) + point) point

-- Calculates a point in the Julia set
pJulia :: Bool -> Double -> Int -> Double -> Complex Double -> Complex Double -> Maybe Double
pJulia aa pixelSize maxIterations power c point =
   fPoint aa pixelSize maxIterations power (\d z -> z ** (d :+ 0) + c) point

fPoint
   :: Bool
   -> Double
   -> Int
   -> Double
   -> (Double -> Complex Double -> Complex Double)
   -> Complex Double
   -> Maybe Double
fPoint bAAEnabled pixelSize maxIterations power f z0 = if not bAAEnabled
   then normalize $ iterPoint 0 z0
   else
      (\(sum, count', isConvergent) ->
            if isConvergent then Nothing else Just $ (sum / (numberOfSubpoints - count'))
         )
         $ foldl
              (\(sum, countConv, isConvergent) y -> if countConv > maxConvergent
                 then (sum, countConv, True)
                 else case normalize . iterPoint 0 $ z0 + y of
                    Nothing -> (sum, countConv + 1, False)
                    Just z  -> (sum + z, countConv, False)
              )
              (0, 0, False)
              [ (-size') :+ (-size'), (-size') :+ 0, (-size') :+ size'
              , 0 :+ (-size'),        0 :+ 0,        0 :+ size'
              , size' :+ (-size'),    size' :+ 0,    size' :+ size'
              ]
 where
  maxMagnitude = 50
  -- Main iteration function
  iterPoint n z | n > maxIterations          = Convergent
                | magnitude z > maxMagnitude = Divergent n z
                | otherwise                  = iterPoint (n + 1) (f power z)
  -- Prevents the aliasing present in generators which color based purely on escape time
  normalize Convergent      = Nothing
  normalize (Divergent 0 _) = Just 0
  normalize (Divergent iterations final) =
     Just . (fromIntegral iterations + 1 -) . logBase power . logBase maxMagnitude . magnitude $ final
  size'             = pixelSize / 3
  numberOfSubpoints = 9
  maxConvergent     = 4

max_stars_repo_name: MaygeKyatt/Fractals

max_issues_repo_name: MaygeKyatt/Fractals

max_forks_repo_name: MaygeKyatt/Fractals


text:
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE GADTs #-}

-- {-# LANGUAGE TypeOperators #-}
-- {-# LANGUAGE DataKinds #-}
-- {-# LANGUAGE KindSignatures #-}
-- {-# LANGUAGE StandaloneDeriving #-}
-- {-# LANGUAGE FlexibleInstances #-}

{-# OPTIONS_GHC -XNoMonomorphismRestriction -O0 -Wtype-defaults -fexpose-all-unfoldings -ddump-simpl -dcore-lint -dsuppress-all -dno-suppress-type-applications -fprint-equality-relations -fplugin=Plugin.MatchPlugin #-}

module Test.PluginExample where

import           GHC.Float
import           Data.List

import           Deep.Expr
import           Deep.Delineate

import           Data.Char (ord)

import           GHC.Generics
import           Data.Void

import           Data.Bifunctor

import           Data.Complex

import           Data.Ratio

import           GHC.Real

import           Data.Proxy

import           CodeGen.C

-- import           Rewrite.Rewrite

-- import           Debug.Trace

-- deriving instance Generic (Ratio a) -- XXX: This instance should probably be in 'base'

default (Int, Double)

{-
data Three = First | Second | Third
  deriving (Generic)

instance GPURep Three

-- test :: Int
-- test =
--   case Second of
--     First -> 10
--     Second -> 100
--     Third -> 1000

-- TODO: Handle 'let's in plugin
test :: Complex Double
test =
  let p = (3.0 :+ 0, 2.0 :+ 0, 1.0 :+ 0) :: (Complex Double, Complex Double, Complex Double)
      {-# INLINE p #-}
  in
  internalize (externalize
    (case p of
       (x, y, re :+ im) ->
          im :+ 0
    ))
{-# NOINLINE test #-}

-- prod :: (Int, Int, (Int, Int))
-- prod = (3, 2, (1, 0))
-- {-# NOINLINE prod #-}

-- test :: (Int, Int)
-- test =
--   internalize (externalize
--     (case prod of
--        (x, y, (a, b)) -> (a, 5)))
-- {-# NOINLINE test #-}

main :: IO ()
main = return ()
-}

data Nat = Z | S Nat deriving (Generic, Show)

data Example = N Nat | B Bool deriving (Generic, Show)

eitherExample :: Int -> Either Int Bool
eitherExample x = Right False

example1 :: Int
example1 =
  internalize (externalize
    (case eitherExample 1 of
      Left x -> x * 2
      Right y -> fromEnum y))
{-# NOINLINE example1 #-}

data Example' = N' Int | B' Bool deriving (Show, Generic)

instance GPURep Example'

-- example2_ :: Int -> Example'
-- example2_ x = B' False
-- {-# NOINLINE example2_ #-}

-- example2 :: Example'
-- example2 =
--   internalize (externalize
--     (case example2_ 0 of
--       N' n -> N' (n+2)
--       B' b -> B' (not b)))
-- {-# NOINLINE example2 #-}

-- example3 :: (Bool, Int) -> Int
-- example3 p =
--   internalize (externalize
--     (case p of
--       (x, y) -> fromEnum x + y))

data Example4 = E1 Int | E2 Float | E3 Bool deriving (Show, Generic)

instance GPURep Example4

example4_ :: Int -> Example4
example4_ x = E2 23.0

-- example4 :: Int
-- example4 =
--   internalize (externalize
--     (case example4_ 0 of
--       E1 x -> 2
--       E2 y -> 4
--       E3 z -> 6))

data Example5 = A1 Float Float | A2 Int deriving (Show, Generic)

instance GPURep Example5

example5_ :: Int -> Example5
example5_ x = A1 2.3 7.5

-- example5 :: Float
-- example5 =
--   internalize (externalize
--     (case example5_ 0 of
--       A2 x -> fromIntegral x
--       A1 x y -> x + y))

-- | isEven

-- example6 :: Int -> Bool
-- example6 x =
--   (internalize (externalize
--     (case x == 0 of
--       True -> True
--       False ->
--         case x == 1 of
--           True -> False
--           False -> example6 (x - 2))))

data IntPair = IntPair Int Int deriving (Show, Generic)

instance GPURep IntPair where

example7_ :: Int -> IntPair
example7_ x = IntPair 1 2

-- example7 :: Int
-- example7 =
--   internalize (externalize
--     (case example7_ 0 of
--       IntPair x y ->
--         if x == 0
--           then y
--           else x))

example8 :: IntPair -> Int
example8 p =
  internalize (externalize
    (case p of
      IntPair x y -> y))

example9 :: IntPair -> Int
example9 p =
  internalize (externalize
    (case p of
      IntPair x y ->
        if x == 0
          then example8 p
          else example9 (IntPair (x-1) (x*y))))

-- floatTest :: Float
-- floatTest =
--   internalize (externalize
--     (case True of
--       False -> 3 * 2
--       True -> 3 * 3))

data ComplexPair where
  ComplexPair :: Complex Double
                 -> Complex Double
                 -> ComplexPair
  deriving (Generic, Show)

instance GPURep ComplexPair

sumComplexPair :: ComplexPair -> Complex Double
sumComplexPair t =
  internalize (externalize
    (case t of
      ComplexPair a b -> a + b))

realSum :: ComplexPair -> Double
realSum p =
  internalize (externalize
    (case p of
      ComplexPair (a :+ _) (b :+ _) ->
        a + b))

mandelbrot_nextZ :: (Complex Double, Complex Double) -> Complex Double
mandelbrot_nextZ t =
  internalize (externalize
    (case t of
      (c, z) -> (z*z) + c))
{-# NOINLINE mandelbrot_nextZ #-}

mandelbrot_helper :: (Int, Complex Double, Complex Double) -> Maybe Int
mandelbrot_helper t =
  internalize (externalize
    (case t of
      (iters, c, z) ->
        if iters == 50
          then Nothing
          else
            case z of
              real :+ imag ->
                if (real*real) + (imag*imag) > 4
                  then Just iters
                  else mandelbrot_helper (iters+1, c, mandelbrot_nextZ (c, z))))
{-# NOINLINE mandelbrot_helper #-}

mandelbrot_point :: Complex Double -> Maybe Int
mandelbrot_point c =
  internalize (externalize
    (mandelbrot_helper (0, c, 0)))

mandelbrot_pointE :: Complex Double -> GPUExp (Maybe Int)
mandelbrot_pointE c =
  externalize
    (mandelbrot_helper (0, c, 0))

data IntList = Nil | Cons Int IntList
  deriving (Generic, Show)

instance GPURep IntList

isEmpty :: IntList -> Bool
isEmpty t =
  internalize (externalize
    (case t of
      Nil -> True
      Cons x xs -> False))

-- total :: (Int, IntList) -> GPUExp Int
total = externalize (\p -> case p of
  (acc, Nil) -> acc
  (acc, Cons x xs) -> total (acc + x, xs))

intListLength_helper :: (Int, IntList) -> Int
intListLength_helper p =
  internalize (externalize
    (case p of
      (acc, t) ->
        case t of
          Nil -> acc
          Cons _ xs -> intListLength_helper (acc+1, xs)))
{-# NOINLINE intListLength_helper #-}

intListLength :: IntList -> Int
intListLength t = internalize (externalize (intListLength_helper (0, t)))

intListSum_helper :: (Int, IntList) -> Int
intListSum_helper p =
  internalize (externalize
    (case p of
       (acc, t) ->
         case t of
           Nil -> acc
           Cons x xs -> intListSum_helper (x+acc, xs)))
{-# NOINLINE intListSum_helper #-}

intListSum :: IntList -> Int
intListSum t =
  internalize (externalize (intListSum_helper (0, t)))

intListSumE :: IntList -> GPUExp Int
intListSumE t = externalize (intListSum_helper (0, t))

example_x :: Either Char Int
example_x = Left 'a'
{-# INLINE example_x #-}

example :: Int
example =
  internalize (externalize
    (case example_x of
      Left  c -> ord c
      Right i -> i))
{-# NOINLINE example #-}

data DoubleList = DNil | DCons Double DoubleList deriving (Generic)

instance GPURep DoubleList

doubleListSum_helper :: (Double, DoubleList) -> Double
doubleListSum_helper p =
  internalize (externalize
    (case p of
      (acc, t) ->
        case t of
          DNil -> acc
          DCons x xs -> doubleListSum_helper (x+acc, xs)))

doubleListSumE :: DoubleList -> GPUExp Double
doubleListSumE t = externalize (doubleListSum_helper (0, t))

-- data FinBitTree (n :: TL.Nat) where
--   FBTLeafZero :: FinBitTree n
--   FBTLeafOne  :: FinBitTree n
--   FBTNode     :: FinBitTree n -> FinBitTree m -> FinBitTree (n+m)

-- deriving instance Generic1 FinBitTree

-- nonterm_test0 :: Int -> Int
-- nonterm_test0 0 = 0
-- nonterm_test0 n = nonterm_test1 (n-1)

-- nonterm_test1 :: Int -> Int
-- nonterm_test1 0 = 1
-- nonterm_test1 n = nonterm_test0 (n-1)

-- nonterm_test0 :: Int -> Int
-- nonterm_test0 x =
--   if x == 0
--     then 0
--     else nonterm_test1 (x-1)
-- -- nonterm_test0 0 = 0
-- -- nonterm_test0 n = nonterm_test1 (n-1)

-- nonterm_test1 :: Int -> Int
-- nonterm_test1 x =
--   if x == 0
--     then 0
--     else nonterm_test0 (x-1)

-- nonterm_test1 0 = 1
-- nonterm_test1 n = nonterm_test0 (n-1)

-- nonterm_test :: Int -> Int
-- nonterm_test x =
--   internalize (externalize
--     (nonterm_test0 4))

main :: IO ()
main = do
  let intList = Cons 10 (Cons 100 (Cons 1000 (Cons 10000 Nil)))
      doubleList = DCons 3.2 (DCons 1 (DCons 2.5 DNil))
  -- print $ mandelbrot_point (1 :+ 0)
  putStrLn $ genProgram (mandelbrot_pointE (1 :+ 0))
  -- putStrLn (genProgram (doubleListSumE doubleList)) --putStrLn (genProgram cTest6)
  -- print (isEmpty Nil, isEmpty (Cons 1 Nil)
  --       ,intListLength intList
  --       ,intListSum intList
  --       ,example
  --       )
  -- print $ realSum (ComplexPair (2 :+ 100) (3 :+ 10000))
  -- putStrLn mandelbrotTestAscii

-- main :: IO ()
-- main =
--   -- print (example6 3, example6 4, example9 (IntPair 5 1))
--   print example1

mandelbrotTestAscii :: String
mandelbrotTestAscii =
  unlines
    (map go [0..mandelbrot_height-1])
  where
    go y = map (go2 y) [0..mandelbrot_width-1]

    go2 y x =
      case mandelbrot_point (mandelbrot_toCoord x y) of
        Just _ -> ' '
        Nothing -> '*'

mandelbrot_toCoord :: Int -> Int -> Complex Double
mandelbrot_toCoord x0 y0 =
    (mandelbrot_xMin + x * mandelbrot_xIncr) :+ (mandelbrot_yMin + y * mandelbrot_yIncr)
  where
    x, y :: Double
    x = fromIntegral x0
    y = fromIntegral y0

mandelbrot_xIncr :: Double
mandelbrot_xIncr = (mandelbrot_xMax - mandelbrot_xMin) / (fromIntegral mandelbrot_width - 1)

mandelbrot_yIncr :: Double
mandelbrot_yIncr = (mandelbrot_yMax - mandelbrot_yMin) / (fromIntegral mandelbrot_height - 1)

mandelbrot_xMin :: Double
mandelbrot_xMin = -2.5

mandelbrot_xMax :: Double
mandelbrot_xMax = 1

mandelbrot_yMin :: Double
mandelbrot_yMin = -1.5

mandelbrot_yMax :: Double
mandelbrot_yMax = 1

mandelbrot_width :: Int
mandelbrot_width = 200

mandelbrot_height :: Int
mandelbrot_height = 40

-- thExample2 :: Q Exp
-- thExample2 = do
--   exp <-
--     [| case (True, the @Int 7) of
--          (x, y) -> fromEnum x + y
--     |]
--   transformExpr exp

-- data Example' = N' Int | B' Bool deriving (Show, Generic)

-- thExample3 :: Q Exp
-- thExample3 = do
--   exp <-
--     [|
--       case B' False of
--         N' n -> N' (n+2)
--         B' b -> B' (not b)
--     |]
--   transformExpr exp


text:
------------------------------------------------------------
--
-- Nums.hs
-- Code sample accompanying topic 1.3.2 "Types and type signatures"
-- See README.md for details
--
-- Fundamentals of Practical Haskell Programming
-- By Richard Cook
--
------------------------------------------------------------

import Data.Complex
import Data.Ratio

n0 :: Int
n0 = 5

n1 :: Double
n1 = 5.0

n2 :: Complex Double
n2 = 2 :+ 3

n3 :: Ratio Int
n3 = 2 % 3

main :: IO ()
main = do
    print n0
    print n1
    print n2
    print n3

max_stars_repo_path: haskell/original/Fundamentals of Practical Haskell Programming [Video]/Volume_01_Section_03_Topic_02/Nums.hs


text:
{-# LANGUAGE GADTs #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# OPTIONS_GHC -Wall #-}
-- Show and Eq instances gives warnings.
{-# OPTIONS_GHC -Wno-orphans #-}

--
-- Copyright (c) 2009-2011, 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.Types
       ( module Feldspar.Core.Types
       ) where

import Data.Array.IO (IOArray)
import Data.Bits
import Data.Complex
import Data.Hash (Hashable(..), combine, hashInt, Hash)
import Data.Int
import Data.IORef
import Data.List
import Data.Typeable (Typeable)
import Data.Word
import qualified Control.Monad.Par as MonadPar
import GHC.TypeLits

import Feldspar.Core.Tuple (Tuply(..))
import Feldspar.Lattice
import Feldspar.Range
import Feldspar.Core.NestedTuples

--------------------------------------------------------------------------------
-- * Heterogenous lists
--------------------------------------------------------------------------------

-- | Heterogeneous list
data a :> b = a :> b
  deriving (Eq, Ord, Show)

infixr 5 :>

instance (Lattice a, Lattice b) => Lattice (a :> b)
  where
    bot = bot :> bot
    top = top :> top
    (a1:>a2) \/ (b1:>b2) = (a1 \/ b1) :> (a2 \/ b2)
    (a1:>a2) /\ (b1:>b2) = (a1 /\ b1) :> (a2 /\ b2)

--------------------------------------------------------------------------------
-- * Integers
--------------------------------------------------------------------------------

-- | Convenience alias for bounded integers
type BoundedInt a = (Bounded a, Integral a, FiniteBits a,
                     Integral (UnsignedRep a), FiniteBits (UnsignedRep a))

-- | Witness for a type level literal. Native is encoded as 0.
data BitWidth (n :: Nat) where
    N8      :: BitWidth 8
    N16     :: BitWidth 16
    N32     :: BitWidth 32
    N64     :: BitWidth 64

instance KnownNat a => Hashable (BitWidth a) where
  hash = hashInt . fromIntegral . natVal

-- | Witness for 'U' or 'S'
data Signedness (s :: Symbol) where
    U :: Signedness "U"
    S :: Signedness "S"

instance KnownSymbol a => Hashable (Signedness a) where
  hash = hash . symbolVal

signedness :: Signedness s -> String
signedness U = "Word"
signedness S = "Int"

-- | A generalization of unsigned and signed integers. The first parameter
-- represents the signedness and the second parameter the number of bits.
type family GenericInt s n where
  GenericInt "U" 8      = Word8
  GenericInt "S" 8      = Int8
  GenericInt "U" 16     = Word16
  GenericInt "S" 16     = Int16
  GenericInt "U" 32     = Word32
  GenericInt "S" 32     = Int32
  GenericInt "U" 64     = Word64
  GenericInt "S" 64     = Int64

-- FIXME: IntN/WordN are hardcoded, not based on Options for architecture.

-- | Signed "natural" representation on an architecture
type IntN = Int32

-- | Unsigned "natural" representation on an architecture
type WordN = Word32

type Length = WordN
type Index  = WordN

--------------------------------------------------------------------------------
-- * Mutable data
--------------------------------------------------------------------------------

-- TODO Make newtypes?

-- | Monad for manipulation of mutable data
type Mut = IO

-- | Mutable references
instance Show (IORef a)
  where
    show _ = "IORef"

-- | Mutable arrays
type MArr a = IOArray Integer a

instance Show (MArr a)
  where
    show _ = "MArr"

--------------------------------------------------------------------------------
-- * Par Monad
--------------------------------------------------------------------------------

-- | Monad for parallel constructs
type Par = MonadPar.Par

-- | Immutable references
type IV = MonadPar.IVar

instance Show (IV a)
  where
    show _ = "IVar"

--------------------------------------------------------------------------------
-- * Elements Language
--------------------------------------------------------------------------------

newtype Elements a = Elements { unE :: [(Index, a)] }

instance Show a => Show (Elements a)
  where
    show (Elements a) = "Elements " ++ show a

instance Eq a => Eq (Elements a)
  where
    _ == _ = False

--------------------------------------------------------------------------------
-- * Future values
--------------------------------------------------------------------------------

newtype FVal a = FVal {unFVal :: a}

instance Show (FVal a)
  where
    show _ = "future"

instance Eq a => Eq (FVal a)
  where
    (FVal a) == (FVal b) = a == b

--------------------------------------------------------------------------------
-- * Type representation
--------------------------------------------------------------------------------

-- | Representation of supported types
data TypeRep a
  where
    UnitType      :: TypeRep ()
    BoolType      :: TypeRep Bool
    IntType       :: ( KnownNat n, KnownSymbol s, BoundedInt (GenericInt s n)
                     , Size (GenericInt s n) ~ Range (GenericInt s n)
                     ) =>
                       Signedness s -> BitWidth n -> TypeRep (GenericInt s n)
    FloatType     :: TypeRep Float
    DoubleType    :: TypeRep Double
    ComplexType   :: RealFloat a => TypeRep a -> TypeRep (Complex a)
    ArrayType     :: TypeRep a -> TypeRep [a]
    Tup2Type      :: TypeRep (Tuple '[a, b]) -> TypeRep (a,b)
    Tup3Type      :: TypeRep (Tuple '[a, b, c]) -> TypeRep (a,b,c)
    Tup4Type      :: TypeRep (Tuple '[a, b, c, d]) -> TypeRep (a, b, c, d)
    Tup5Type      :: TypeRep (Tuple '[a, b, c, d, e]) -> TypeRep (a, b, c, d, e)
    Tup6Type      :: TypeRep (Tuple '[a, b, c, d, e, f])
                  -> TypeRep (a, b, c, d, e, f)
    Tup7Type      :: TypeRep (Tuple '[a, b, c, d, e, f, g])
                  -> TypeRep (a, b, c, d, e, f, g)
    Tup8Type      :: TypeRep (Tuple '[a, b, c, d, e, f, g, h])
                  -> TypeRep (a, b, c, d, e, f, g, h)
    Tup9Type      :: TypeRep (Tuple '[a, b, c, d, e, f, g, h, i])
                  -> TypeRep (a, b, c, d, e, f, g, h, i)
    Tup10Type     :: TypeRep (Tuple '[a, b, c, d, e, f, g, h, i, j])
                  -> TypeRep (a, b, c, d, e, f, g, h, i, j)
    Tup11Type     :: TypeRep (Tuple '[a, b, c, d, e, f, g, h, i, j, k])
                  -> TypeRep (a, b, c, d, e, f, g, h, i, j, k)
    Tup12Type     :: TypeRep (Tuple '[a, b, c, d, e, f, g, h, i, j, k, l])
                  -> TypeRep (a, b, c, d, e, f, g, h, i, j, k, l)
    Tup13Type     :: TypeRep (Tuple '[a, b, c, d, e, f, g, h, i, j, k, l, m])
                  -> TypeRep (a, b, c, d, e, f, g, h, i, j, k, l, m)
    Tup14Type     :: TypeRep (Tuple '[a, b, c, d, e, f, g, h, i, j, k, l, m, n])
                  -> TypeRep (a, b, c, d, e, f, g, h, i, j, k, l, m, n)
    Tup15Type     :: TypeRep (Tuple '[a, b, c, d, e, f, g, h, i, j, k, l, m, n, o])
                  -> TypeRep (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o)
    FunType       :: TypeRep a -> TypeRep b -> TypeRep (a -> b)
    MutType       :: TypeRep a -> TypeRep (Mut a)
    RefType       :: TypeRep a -> TypeRep (IORef a)
    MArrType      :: TypeRep a -> TypeRep (MArr a)
    ParType       :: TypeRep a -> TypeRep (Par a)
    ElementsType  :: TypeRep a -> TypeRep (Elements a)
    ConsType      :: TypeRep a -> TypeRep (Tuple b) -> TypeRep (Tuple (a ': b))
    NilType       :: TypeRep (Tuple '[])
    IVarType      :: TypeRep a -> TypeRep (IV a)
    FValType      :: TypeRep a -> TypeRep (FVal a)
      -- TODO `MArrType` Should have a target-specialized version. Or perhaps
      --      use a single type with a flag to distinguish between immutable and
      --      mutable arrays.

instance Show (TypeRep a)
  where
    show UnitType            = "()"
    show BoolType            = "Bool"
    show (IntType s n)       = signedness s ++ show (natVal n)
    show FloatType           = "Float"
    show DoubleType          = "Double"
    show (ComplexType t)     = "(Complex " ++ show t ++ ")"
    show (ArrayType t)       = "[" ++ show t ++ "]"
    show (Tup2Type t)        = show t
    show (Tup3Type t)        = show t
    show (Tup4Type t)        = show t
    show (Tup5Type t)        = show t
    show (Tup6Type t)        = show t
    show (Tup7Type t)        = show t
    show (Tup8Type t)        = show t
    show (Tup9Type t)        = show t
    show (Tup10Type t)       = show t
    show (Tup11Type t)       = show t
    show (Tup12Type t)       = show t
    show (Tup13Type t)       = show t
    show (Tup14Type t)       = show t
    show (Tup15Type t)       = show t
    show (FunType ta tb)                 = show ta ++ " -> " ++ show tb
    show (MutType ta)                    = unwords ["Mut", show ta]
    show (RefType ta)                    = unwords ["Ref", show ta]
    show (MArrType ta)                   = unwords ["MArr", show ta]
    show (ParType ta)                    = unwords ["Par", show ta]
    show (ElementsType ta)               = unwords ["Elements", show ta]
    show (ConsType ta tb)                = unwords [show ta, ":*", show tb]
    show NilType                         = "TNil"
    show (IVarType ta)                   = unwords ["IVar", show ta]
    show (FValType ta)                   = unwords ["FVal", show ta]

infixl 5 #
(#) :: Hashable a => Hash -> a -> Hash
h # x = h `combine` hash x

instance Hashable (TypeRep a) where
  hash UnitType              = hashInt 1
  hash BoolType              = hashInt 2
  hash (IntType sgn sz)      = hashInt 3 # sgn # sz
  hash FloatType             = hashInt 4
  hash DoubleType            = hashInt 5
  hash (ComplexType t)       = hashInt 6 # t
  hash (ArrayType t)         = hashInt 7 # t
  hash (Tup2Type t)          = hashInt 8 # t
  hash (Tup3Type t)          = hashInt 9 # t
  hash (Tup4Type t)          = hashInt 10 # t
  hash (Tup5Type t)          = hashInt 11 # t
  hash (Tup6Type t)          = hashInt 12 # t
  hash (Tup7Type t)          = hashInt 13 # t
  hash (Tup8Type t)          = hashInt 14 # t
  hash (Tup9Type t)          = hashInt 15 # t
  hash (Tup10Type t)         = hashInt 16 # t
  hash (Tup11Type t)         = hashInt 17 # t
  hash (Tup12Type t)         = hashInt 18 # t
  hash (Tup13Type t)         = hashInt 19 # t
  hash (Tup14Type t)         = hashInt 20 # t
  hash (Tup15Type t)         = hashInt 21 # t
  hash (FunType a b)         = hashInt 22 # a # b
  hash (MutType t)           = hashInt 23 # t
  hash (RefType t)           = hashInt 24 # t
  hash (MArrType t)          = hashInt 25 # t
  hash (ParType t)           = hashInt 26 # t
  hash (ElementsType t)      = hashInt 27 # t
  hash (ConsType a b)        = hashInt 28 # a # b
  hash NilType               = hashInt 29
  hash (IVarType t)          = hashInt 30 # t
  hash (FValType t)          = hashInt 31 # t

-- | The set of supported types
class (Eq a, Show a, Typeable a, Show (Size a), Lattice (Size a), TypeF a) => Type a
  where
    -- | Gives the type representation a value.
    typeRep  :: TypeRep a
    sizeOf   :: a -> Size a

instance Type ()      where typeRep = UnitType;          sizeOf = singletonRange
instance Type Bool    where typeRep = BoolType;          sizeOf = singletonRange
instance Type Word8   where typeRep = IntType U N8;      sizeOf = singletonRange
instance Type Int8    where typeRep = IntType S N8;      sizeOf = singletonRange
instance Type Word16  where typeRep = IntType U N16;     sizeOf = singletonRange
instance Type Int16   where typeRep = IntType S N16;     sizeOf = singletonRange
instance Type Word32  where typeRep = IntType U N32;     sizeOf = singletonRange
instance Type Int32   where typeRep = IntType S N32;     sizeOf = singletonRange
instance Type Word64  where typeRep = IntType U N64;     sizeOf = singletonRange
instance Type Int64   where typeRep = IntType S N64;     sizeOf = singletonRange
instance Type Float   where typeRep = FloatType;         sizeOf = singletonRange
instance Type Double  where typeRep = DoubleType;        sizeOf = singletonRange

instance Hashable a => Hashable (Complex a) where
  hash (re :+ im) = hash re `combine` hash im

instance (Type a, RealFloat a) => Type (Complex a)
  where
    typeRep             = ComplexType typeRep
    sizeOf _            = AnySize

instance Type a => Type [a]
  where
    typeRep       = ArrayType typeRep
    sizeOf as     = singletonRange (genericLength as) :> unions (map sizeOf as)

instance (Type a, Type b) => Type (a, b) where
  typeRep = Tup2Type typeRep
  sizeOf = sizeOf . unpack

instance (Type a, Type b, Type c) => Type (a, b, c) where
  typeRep = Tup3Type typeRep
  sizeOf = sizeOf . unpack

instance (Type a, Type b, Type c, Type d) => Type (a, b, c, d) where
  typeRep = Tup4Type typeRep
  sizeOf = sizeOf . unpack

instance (Type a, Type b, Type c, Type d, Type e) => Type (a, b, c, d, e) where
  typeRep = Tup5Type typeRep
  sizeOf = sizeOf . unpack

instance (Type a, Type b, Type c, Type d, Type e, Type f)
         => Type (a, b, c, d, e, f) where
  typeRep = Tup6Type typeRep
  sizeOf = sizeOf . unpack

instance (Type a, Type b, Type c, Type d, Type e, Type f, Type g)
         => Type (a, b, c, d, e, f, g) where
  typeRep = Tup7Type typeRep
  sizeOf = sizeOf . unpack

instance (Type a, Type b, Type c, Type d, Type e, Type f, Type g, Type h)
         => Type (a, b, c, d, e, f, g, h) where
  typeRep = Tup8Type typeRep
  sizeOf = sizeOf . unpack

instance (Type a, Type b, Type c, Type d, Type e, Type f, Type g, Type h,
          Type i)
         => Type (a, b, c, d, e, f, g, h, i) where
  typeRep = Tup9Type typeRep
  sizeOf = sizeOf . unpack

instance (Type a, Type b, Type c, Type d, Type e, Type f, Type g, Type h,
          Type i, Type j)
         => Type (a, b, c, d, e, f, g, h, i, j) where
  typeRep = Tup10Type typeRep
  sizeOf = sizeOf . unpack

instance (Type a, Type b, Type c, Type d, Type e, Type f, Type g, Type h,
          Type i, Type j, Type k)
         => Type (a, b, c, d, e, f, g, h, i, j, k) where
  typeRep = Tup11Type typeRep
  sizeOf = sizeOf . unpack

instance (Type a, Type b, Type c, Type d, Type e, Type f, Type g, Type h,
          Type i, Type j, Type k, Type l)
         => Type (a, b, c, d, e, f, g, h, i, j, k, l) where
  typeRep = Tup12Type typeRep
  sizeOf = sizeOf . unpack

instance (Type a, Type b, Type c, Type d, Type e, Type f, Type g, Type h,
          Type i, Type j, Type k, Type l, Type m)
         => Type (a, b, c, d, e, f, g, h, i, j, k, l, m) where
  typeRep = Tup13Type typeRep
  sizeOf = sizeOf . unpack

instance (Type a, Type b, Type c, Type d, Type e, Type f, Type g, Type h,
          Type i, Type j, Type k, Type l, Type m, Type n')
         => Type (a, b, c, d, e, f, g, h, i, j, k, l, m, n') where
  typeRep = Tup14Type typeRep
  sizeOf = sizeOf . unpack

instance (Type a, Type b, Type c, Type d, Type e, Type f, Type g, Type h,
          Type i, Type j, Type k, Type l, Type m, Type n', Type o)
         => Type (a, b, c, d, e, f, g, h, i, j, k, l, m, n', o) where
  typeRep = Tup15Type typeRep
  sizeOf = sizeOf . unpack

instance Type a => Type (IORef a)
  where
    typeRep = RefType typeRep

    sizeOf _ = universal

instance Type a => Type (MArr a)
  where
    typeRep = MArrType typeRep

    sizeOf _ = universal

instance Type a => Type (Elements a)
  where
    typeRep = ElementsType typeRep

    sizeOf _ = universal

instance (Type a, Typeable b, Type (Tuple b)) => Type (Tuple (a ': b)) where
    typeRep = ConsType typeRep typeRep
    sizeOf (x :* xs) = (sizeOf x, sizeOf xs)

instance Type (Tuple '[]) where
    typeRep = NilType
    sizeOf _ = universal

instance Type a => Type (IV a)
  where
    typeRep = IVarType typeRep

    sizeOf _ = universal

instance Type a => Type (FVal a)
  where
    typeRep = FValType typeRep

    sizeOf _ = universal

-- | Extend the class Type to higher order types
class (Typeable a, Show (Size a), Lattice (Size a)) => TypeF a where
  typeRepF :: TypeRep a
  sizeOfF :: a -> Size a

instance {-# OVERLAPPING #-} (TypeF a, TypeF b) => TypeF (a -> b) where
  typeRepF = FunType typeRepF typeRepF
  sizeOfF _ = universal

instance {-# OVERLAPPABLE #-} (Typeable a, Type a) => TypeF a where
  typeRepF = typeRep
  sizeOfF = sizeOf

--------------------------------------------------------------------------------
-- * Sized types
--------------------------------------------------------------------------------

data AnySize = AnySize
    deriving (Eq, Show, Ord)

anySizeFun :: AnySize -> AnySize
anySizeFun AnySize = AnySize

anySizeFun2 :: AnySize -> AnySize -> AnySize
anySizeFun2 AnySize AnySize = AnySize

instance Num AnySize
  where
    fromInteger _ = AnySize
    abs           = anySizeFun
    signum        = anySizeFun
    (+)           = anySizeFun2
    (-)           = anySizeFun2
    (*)           = anySizeFun2

instance Lattice AnySize
  where
    bot  = AnySize
    top  = AnySize
    (\/) = anySizeFun2
    (/\) = anySizeFun2

type family Size a where
  Size ()              = Range ()
  Size Bool            = Range Bool
  Size Word8           = Range Word8
  Size Int8            = Range Int8
  Size Word16          = Range Word16
  Size Int16           = Range Int16
  Size Word32          = Range Word32
  Size Int32           = Range Int32
  Size Word64          = Range Word64
  Size Int64           = Range Int64
  Size Float           = Range Float
  Size Double          = Range Double
  Size (Complex a)     = AnySize
  Size [a]             = Range Length :> Size a
  Size (a, b)          = Size (Tuple '[a, b])
  Size (a, b, c)       = Size (Tuple '[a, b, c])
  Size (a, b, c, d)    = Size (Tuple '[a, b, c, d])
  Size (a, b, c, d, e) = Size (Tuple '[a, b, c, d, e])
  Size (a, b, c, d, e, f) = Size (Tuple '[a, b, c, d, e, f])
  Size (a, b, c, d, e, f, g) = Size (Tuple '[a, b, c, d, e, f, g])
  Size (a, b, c, d, e, f, g, h) = Size (Tuple '[a, b, c, d, e, f, g, h])
  Size (a, b, c, d, e, f, g, h, i) = Size (Tuple '[a, b, c, d, e, f, g, h, i])
  Size (a, b, c, d, e, f, g, h, i, j)
    = Size (Tuple '[a, b, c, d, e, f, g, h, i, j])
  Size (a, b, c, d, e, f, g, h, i, j, k)
    = Size (Tuple '[a, b, c, d, e, f, g, h, i, j, k])
  Size (a, b, c, d, e, f, g, h, i, j, k, l)
    = Size (Tuple '[a, b, c, d, e, f, g, h, i, j, k, l])
  Size (a, b, c, d, e, f, g, h, i, j, k, l, m)
    = Size (Tuple '[a, b, c, d, e, f, g, h, i, j, k, l, m])
  Size (a, b, c, d, e, f, g, h, i, j, k, l, m, n)
    = Size (Tuple '[a, b, c, d, e, f, g, h, i, j, k, l, m, n])
  Size (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o)
    = Size (Tuple '[a, b, c, d, e, f, g, h, i, j, k, l, m, n, o])
  Size (a -> b)        = (Size a, Size b)
  Size (Mut a)         = Size a
  Size (IORef a)       = Size a
  Size (MArr a)        = Range Length :> Size a
  Size (Par a)         = Size a
  Size (Elements a)    = Range Length :> Size a
  Size (Tuple (a ': b)) = (Size a, Size (Tuple b))
  Size (Tuple '[])     = Size ()
  Size (IV a)          = Size a
  Size (FVal a)        = Size a

max_stars_repo_name: Feldspar/feldspar-language

max_issues_repo_name: Feldspar/feldspar-language

max_forks_repo_name: Feldspar/feldspar-language


text:
module Permutation where

import Rotation (MatrixRotation, identity)
import qualified Rotation as Rotation

import qualified Numeric.LinearAlgebra as LA

type N = Integer

toNumbers :: MatrixRotation -> [N]
toNumbers = map (go 0 . LA.toList) . LA.toColumns
  where
    go n [] = n
    go n (f:fs)
      | f /= 0 = n
      | otherwise = go (n + 1) fs

facs :: [N]
facs =
  let n = 9 * 6
  in map fac [n - 1, n - 2 .. 0]
  where
    fac 0 = 1
    fac n = n * fac (n - 1)

multiplierMappings :: [N -> N]
multiplierMappings = map mappingFor [0..]
  where
    mappingFor :: Integer -> Integer -> Integer
    mappingFor n m
      | n == m = 0
      | n < m = m
      | otherwise = m + 1

fromRotation :: MatrixRotation -> N
fromRotation = sum . zipWith3 go facs multiplierMappings . toNumbers
  where
    go fac multiplierMapping n = fac * multiplierMapping n

-- toRotation :: N -> MatrixRotation
toRotation = go facs
  where
    go [] _ = []
    go (f:fs) n
      | n >= f =
        let (d, r) = n `divMod` f
        in d : go fs r
      | otherwise = 0 : go fs n

max_stars_repo_name: runjak/hRubiks

max_issues_repo_name: runjak/hRubiks

max_forks_repo_name: runjak/hRubiks


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}

-- |
-- Module      :  Spiral.Spiral.SPL.Run
-- Copyright   :  (c) 2016-2017 Drexel University
-- License     :  BSD-style
-- Maintainer  :  mainland@drexel.edu

module Spiral.SPL.Run (
    toProgram,
    runSPL
  ) where

import Prelude hiding ((!!), read)

import Control.Monad (unless,
                      when)
import Data.Complex
import qualified Data.Vector as V
import Text.PrettyPrint.Mainland
import Text.PrettyPrint.Mainland.Class

import qualified Spiral.Array as A
import Spiral.Array (Computable,
                     DS,
                     MArray,
                     SArray,
                     Vector,
                     computeP,
                     forShapeP,
                     indexS,
                     write)
import Spiral.Array.Operators.Matrix
import Spiral.Array.Operators.Permute
import Spiral.Array.Repr.Complex
import Spiral.Array.Repr.Compute
import Spiral.Array.Repr.Concrete
import Spiral.Array.Repr.Slice
import Spiral.Array.Repr.Transform
import Spiral.Array.Shape
import Spiral.Exp
import Spiral.Monad
import Spiral.Program.Monad
import Spiral.Program.Syntax
import Spiral.SPL
import Spiral.Util.Trace

-- | Generate code for an SPL transform.
toProgram :: forall a m . (Num (Exp a), Typed a, MonadSpiral m)
          => String
          -> SPL (Exp a)
          -> m (Program a)
toProgram f a = do
    stms <- evalP $ runSPL a (fromGather x) >>= computeP y
    return $ Program f x y stms
   where
     x, y :: Vector C (Exp a)
     x = C (ix1 n) "X"
     y = C (ix1 m) "Y"

     m, n :: Int
     Z :. m :. n = extent a

infix 4 .<-.
(.<-.) :: forall r1 r2 m a . (MonadSpiral m, MArray r1 DIM1 a, Computable r2 DIM1 a)
       => Vector r1 a
       -> P m (Vector r2 a)
       -> P m ()
y .<-. k = do
    x <- k
    computeP y x

runSPL :: forall m a .
          ( MonadSpiral m
          , Typed a
          , Num (Exp a)
          )
       => SPL (Exp a)
       -> Vector T (Exp a)
       -> P m (Vector T (Exp a))
runSPL e@I{} x = do
    comment $ ppr e
    return x

runSPL e@(Pi p) x = do
    comment $ ppr e
    permuteP p x

runSPL e@(Diag v) x = do
    t <- gather x
    computeTransform e $ \y -> do
      comment $ ppr e
      shouldUnroll n >>= go t y
  where
    n :: Int
    n = V.length v

    go :: forall r1 r2 m . (MonadSpiral m, SArray r1 DIM1 (Exp a), MArray r2 DIM1 (Exp a))
       => Array r1 DIM1 (Exp a)
       -> Array r2 DIM1 (Exp a)
       -> Bool
       -> P m ()
    go x y True =
        forP 0 n $ \ei@(ConstE (IntC i)) ->
            write y (Z :. ei) (v V.! i * indexS x (Z :. ei))

    go x y _ = do
        d <- cacheArray (A.fromFunction (ix1 n) f)
        forShapeP (ix1 n) $ \eix ->
            write y eix (indexS d eix * indexS x eix)
      where
        f (Z :. i) = v V.! i

runSPL e@(KDiag n k) x = do
    t <- gather x
    computeTransform e $ \y -> do
      comment $ ppr e
      forShapeP (ix1 n) $ \eix ->
          write y eix (k * indexS t eix)

runSPL e@(Kron (I m) a) x = do
    t <- gather x
    computeTransform e $ \y -> do
      comment $ ppr e
      forP 0 m $ \i ->
        slice y (i*e_n) 1 n .<-. runSPL a (fromGather (slice t (i*e_n') 1 n'))
  where
    Z :. n :. n' = extent a

    e_n, e_n' :: Exp Int
    e_n  = fromIntegral n
    e_n' = fromIntegral n'

runSPL e@(Kron a (I n)) x = do
    t <- gather x
    computeTransform e $ \y -> do
      comment $ ppr e
      forP 0 n $ \i ->
        slice y i n m .<-. runSPL a (fromGather (slice t i n m'))
  where
    Z :. m :. m' = extent a

runSPL (Kron a b) x =
    runSPL ((I m ⊗ b) × (a ⊗ I n')) x
  where
    Z :. m :. _  = extent a
    Z :. _ :. n' = extent b

runSPL e@(DSum a b) x = do
    t <- gather x
    computeTransform e $ \y -> do
      comment $ ppr e
      slice y 0 1 m                .<-. runSPL a  (fromGather (slice t 0 1 m'))
      slice y (fromIntegral m) 1 n .<-. runSPL b  (fromGather (slice t (fromIntegral m') 1 n'))
  where
    Z :. m :. m' = extent a
    Z :. n :. n' = extent b

runSPL e@(Prod a b) x = do
    when (n' /= n) $
        faildoc $ "Incompatible matrix product:" <+> ppr n <+> char 'x' <+> ppr n' </> ppr e
    comment $ ppr e
    t <- runSPL b x
    runSPL a t
  where
    Z :. _m :.  n = extent a
    Z :. n' :. _p = extent b

runSPL (Re a) x = do
    t <- gather x
    computeTransform (Re a) $ \y -> do
      let t' :: Vector (CMPLX DS) (Exp (Complex a)) = toComplexArray t
      toComplexArray y .<-. runSPL a (fromGather t')

runSPL a x = do
    unless (shouldDefault a) $
      traceCg $ text "Falling back to default compilation path:" </> ppr a
    t <- gather x
    computeTransform a $ \y -> do
      comment $ ppr a
      mXvP (toMatrix a) t y
  where
    shouldDefault :: SPL (Exp a) -> Bool
    shouldDefault E{}  = True
    shouldDefault F2{} = True
    shouldDefault _    = False

computeTransform :: MonadSpiral m
                 => SPL a
                 -> (forall r m . (MArray r DIM1 a, MonadSpiral m) => Array r DIM1 a -> P m ())
                 -> P m (Array T DIM1 a)
computeTransform a k =
    return $ fromScatter $ fromCompute (ix1 m) k
  where
    Z :. m :. _n = extent a

max_stars_repo_name: mainland/hspiral

max_issues_repo_name: mainland/hspiral

max_forks_repo_name: mainland/hspiral


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
-- |
-- Module    : Statistics.Distribution.Gamma
-- Copyright : (c) 2009, 2011 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- The gamma distribution.  This is a continuous probability
-- distribution with two parameters, /k/ and ϑ. If /k/ is
-- integral, the distribution represents the sum of /k/ independent
-- exponentially distributed random variables, each of which has a
-- mean of ϑ.

module Statistics.Distribution.Gamma
    (
      GammaDistribution
    -- * Constructors
    , gammaDistr
    , gammaDistrE
    , improperGammaDistr
    , improperGammaDistrE
    -- * Accessors
    , gdShape
    , gdScale
    ) where

import Control.Applicative
import Data.Aeson           (FromJSON(..), ToJSON, Value(..), (.:))
import Data.Binary          (Binary(..))
import Data.Data            (Data, Typeable)
import GHC.Generics         (Generic)
import Numeric.MathFunctions.Constants (m_pos_inf, m_NaN, m_neg_inf)
import Numeric.SpecFunctions (incompleteGamma, invIncompleteGamma, logGamma, digamma)
import qualified System.Random.MWC.Distributions as MWC

import Statistics.Distribution.Poisson.Internal as Poisson
import qualified Statistics.Distribution as D
import Statistics.Internal

-- | The gamma distribution.
data GammaDistribution = GD {
      gdShape :: {-# UNPACK #-} !Double -- ^ Shape parameter, /k/.
    , gdScale :: {-# UNPACK #-} !Double -- ^ Scale parameter, ϑ.
    } deriving (Eq, Typeable, Data, Generic)

instance Show GammaDistribution where
  showsPrec i (GD k theta) = defaultShow2 "improperGammaDistr" k theta i
instance Read GammaDistribution where
  readPrec = defaultReadPrecM2 "improperGammaDistr" improperGammaDistrE

instance ToJSON GammaDistribution
instance FromJSON GammaDistribution where
  parseJSON (Object v) = do
    k     <- v .: "gdShape"
    theta <- v .: "gdScale"
    maybe (fail $ errMsgI k theta) return $ improperGammaDistrE k theta
  parseJSON _ = empty

instance Binary GammaDistribution where
  put (GD x y) = put x >> put y
  get = do
    k     <- get
    theta <- get
    maybe (fail $ errMsgI k theta) return $ improperGammaDistrE k theta

-- | Create gamma distribution. Both shape and scale parameters must
-- be positive.
gammaDistr :: Double            -- ^ Shape parameter. /k/
           -> Double            -- ^ Scale parameter, ϑ.
           -> GammaDistribution
gammaDistr k theta
  = maybe (error $ errMsg k theta) id $ gammaDistrE k theta

errMsg :: Double -> Double -> String
errMsg k theta
  =  "Statistics.Distribution.Gamma.gammaDistr: "
  ++ "k=" ++ show k
  ++ "theta=" ++ show theta
  ++ " but must be positive"

-- | Create gamma distribution. Both shape and scale parameters must
-- be positive.
gammaDistrE :: Double            -- ^ Shape parameter. /k/
            -> Double            -- ^ Scale parameter, ϑ.
            -> Maybe GammaDistribution
gammaDistrE k theta
  | k > 0 && theta > 0 = Just (GD k theta)
  | otherwise          = Nothing

-- | Create gamma distribution. Both shape and scale parameters must
-- be non-negative.
improperGammaDistr :: Double            -- ^ Shape parameter. /k/
                   -> Double            -- ^ Scale parameter, ϑ.
                   -> GammaDistribution
improperGammaDistr k theta
  = maybe (error $ errMsgI k theta) id $ improperGammaDistrE k theta

errMsgI :: Double -> Double -> String
errMsgI k theta
  =  "Statistics.Distribution.Gamma.gammaDistr: "
  ++ "k=" ++ show k
  ++ "theta=" ++ show theta
  ++ " but must be non-negative"

-- | Create gamma distribution. Both shape and scale parameters must
-- be non-negative.
improperGammaDistrE :: Double            -- ^ Shape parameter. /k/
                    -> Double            -- ^ Scale parameter, ϑ.
                    -> Maybe GammaDistribution
improperGammaDistrE k theta
  | k >= 0 && theta >= 0 = Just (GD k theta)
  | otherwise            = Nothing

instance D.Distribution GammaDistribution where
    cumulative = cumulative

instance D.ContDistr GammaDistribution where
    density    = density
    logDensity (GD k theta) x
      | x <= 0    = m_neg_inf
      | otherwise = log x * (k - 1) - (x / theta) - logGamma k - log theta * k
    quantile   = quantile

instance D.Variance GammaDistribution where
    variance (GD a l) = a * l * l

instance D.Mean GammaDistribution where
    mean (GD a l) = a * l

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

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

instance D.MaybeEntropy GammaDistribution where
  maybeEntropy (GD a l)
    | a > 0 && l > 0 =
      Just $
      a
      + log l
      + logGamma a
      + (1-a) * digamma a
    | otherwise = Nothing

instance D.ContGen GammaDistribution where
    genContVar (GD a l) = MWC.gamma a l

density :: GammaDistribution -> Double -> Double
density (GD a l) x
  | a < 0 || l <= 0   = m_NaN
  | x <= 0            = 0
  | a == 0            = if x == 0 then m_pos_inf else 0
  | x == 0            = if a < 1 then m_pos_inf else if a > 1 then 0 else 1/l
  | a < 1             = Poisson.probability (x/l) a * a / x
  | otherwise         = Poisson.probability (x/l) (a-1) / l

cumulative :: GammaDistribution -> Double -> Double
cumulative (GD k l) x
  | x <= 0    = 0
  | otherwise = incompleteGamma k (x/l)

quantile :: GammaDistribution -> Double -> Double
quantile (GD k l) p
  | p == 0         = 0
  | p == 1         = 1/0
  | p > 0 && p < 1 = l * invIncompleteGamma k p
  | otherwise      =
    error $ "Statistics.Distribution.Gamma.quantile: p must be in [0,1] range. Got: "++show p


text:
-- |
-- Module    : Statistics.Matrix
-- Copyright : 2011 Aleksey Khudyakov, 2014 Bryan O'Sullivan
-- License   : BSD3
--
-- Basic matrix operations.
--
-- There isn't a widely used matrix package for Haskell yet, so
-- we implement the necessary minimum here.

module Statistics.Matrix
    (
      Matrix(..)
    , Vector
    , fromList
    , fromVector
    , toVector
    , toList
    , dimension
    , center
    , multiply
    , multiplyV
    , transpose
    , power
    , norm
    , column
    , row
    , map
    , for
    , unsafeIndex
    , hasNaN
    , bounds
    , unsafeBounds
    ) where

import Prelude hiding (exponent, map, sum)
import Statistics.Function (for, square)
import Statistics.Matrix.Types
import Statistics.Sample.Internal (sum)
import qualified Data.Vector.Unboxed as U

-- | Convert from a row-major list.
fromList :: Int                 -- ^ Number of rows.
         -> Int                 -- ^ Number of columns.
         -> [Double]            -- ^ Flat list of values, in row-major order.
         -> Matrix
fromList r c = fromVector r c . U.fromList

-- | Convert from a row-major vector.
fromVector :: Int               -- ^ Number of rows.
           -> Int               -- ^ Number of columns.
           -> U.Vector Double   -- ^ Flat list of values, in row-major order.
           -> Matrix
fromVector r c v
  | r*c /= len = error "input size mismatch"
  | otherwise  = Matrix r c 0 v
  where len    = U.length v

-- | Convert to a row-major flat vector.
toVector :: Matrix -> U.Vector Double
toVector (Matrix _ _ _ v) = v

-- | Convert to a row-major flat list.
toList :: Matrix -> [Double]
toList = U.toList . toVector

-- | Return the dimensions of this matrix, as a (row,column) pair.
dimension :: Matrix -> (Int, Int)
dimension (Matrix r c _ _) = (r, c)

-- | Avoid overflow in the matrix.
avoidOverflow :: Matrix -> Matrix
avoidOverflow m@(Matrix r c e v)
  | center m > 1e140 = Matrix r c (e + 140) (U.map (* 1e-140) v)
  | otherwise        = m

-- | Matrix-matrix multiplication. Matrices must be of compatible
-- sizes (/note: not checked/).
multiply :: Matrix -> Matrix -> Matrix
multiply m1@(Matrix r1 _ e1 _) m2@(Matrix _ c2 e2 _) =
  Matrix r1 c2 (e1 + e2) $ U.generate (r1*c2) go
  where
    go t = sum $ U.zipWith (*) (row m1 i) (column m2 j)
      where (i,j) = t `quotRem` c2

-- | Matrix-vector multiplication.
multiplyV :: Matrix -> Vector -> Vector
multiplyV m v
  | cols m == c = U.generate (rows m) (sum . U.zipWith (*) v . row m)
  | otherwise   = error $ "matrix/vector unconformable " ++ show (cols m,c)
  where c = U.length v

-- | Raise matrix to /n/th power. Power must be positive
-- (/note: not checked).
power :: Matrix -> Int -> Matrix
power mat 1 = mat
power mat n = avoidOverflow res
  where
    mat2 = power mat (n `quot` 2)
    pow  = multiply mat2 mat2
    res | odd n     = multiply pow mat
        | otherwise = pow

-- | Element in the center of matrix (not corrected for exponent).
center :: Matrix -> Double
center mat@(Matrix r c _ _) =
    unsafeBounds U.unsafeIndex mat (r `quot` 2) (c `quot` 2)

-- | Calculate the Euclidean norm of a vector.
norm :: Vector -> Double
norm = sqrt . sum . U.map square

-- | Return the given column.
column :: Matrix -> Int -> Vector
column (Matrix r c _ v) i = U.backpermute v $ U.enumFromStepN i c r
{-# INLINE column #-}

-- | Return the given row.
row :: Matrix -> Int -> Vector
row (Matrix _ c _ v) i = U.slice (c*i) c v

unsafeIndex :: Matrix
            -> Int              -- ^ Row.
            -> Int              -- ^ Column.
            -> Double
unsafeIndex = unsafeBounds U.unsafeIndex

map :: (Double -> Double) -> Matrix -> Matrix
map f (Matrix r c e v) = Matrix r c e (U.map f v)

-- | Indicate whether any element of the matrix is @NaN@.
hasNaN :: Matrix -> Bool
hasNaN = U.any isNaN . toVector

-- | Given row and column numbers, calculate the offset into the flat
-- row-major vector.
bounds :: (Vector -> Int -> r) -> Matrix -> Int -> Int -> r
bounds k (Matrix rs cs _ v) r c
  | r < 0 || r >= rs = error "row out of bounds"
  | c < 0 || c >= cs = error "column out of bounds"
  | otherwise        = k v $! r * cs + c
{-# INLINE bounds #-}

-- | Given row and column numbers, calculate the offset into the flat
-- row-major vector, without checking.
unsafeBounds :: (Vector -> Int -> r) -> Matrix -> Int -> Int -> r
unsafeBounds k (Matrix _ cs _ v) r c = k v $! r * cs + c
{-# INLINE unsafeBounds #-}

transpose :: Matrix -> Matrix
transpose m@(Matrix r0 c0 e _) = Matrix c0 r0 e . U.generate (r0*c0) $ \i ->
  let (r,c) = i `quotRem` r0
  in unsafeIndex m c r

max_stars_repo_name: StefanHubner/statistics

max_issues_repo_name: StefanHubner/statistics

max_forks_repo_name: StefanHubner/statistics


text:
{-# OPTIONS_GHC -Wall #-}

import           Numeric.Sundials.ARKode.ODE
import           Numeric.LinearAlgebra

import           Data.List (intercalate)

import           Text.PrettyPrint.HughesPJClass

butcherTableauTex :: ButcherTable -> String
butcherTableauTex (ButcherTable m c b b2) =
  render $
  vcat [ text ("\n\\begin{array}{c|" ++ (concat $ replicate n "c") ++ "}")
       , us
       , text "\\hline"
       , text bs <+> text "\\\\"
       , text b2s <+> text "\\\\"
       , text "\\end{array}"
       ]
  where
    n = rows m
    rs = toLists m
    ss = map (\r -> intercalate " & " $ map show r) rs
    ts = zipWith (\i r -> show i ++ " & " ++ r) (toList c) ss
    us = vcat $ map (\r -> text r <+> text "\\\\") ts
    bs  = " & " ++ (intercalate " & " $ map show $ toList b)
    b2s = " & " ++ (intercalate " & " $ map show $ toList b2)

main :: IO ()
main = do

  let res = butcherTable (SDIRK_2_1_2 undefined)
  putStrLn $ show res
  putStrLn $ butcherTableauTex res

  let resA = butcherTable (KVAERNO_4_2_3 undefined)
  putStrLn $ show resA
  putStrLn $ butcherTableauTex resA

  let resB = butcherTable (SDIRK_5_3_4 undefined)
  putStrLn $ show resB
  putStrLn $ butcherTableauTex resB

  let resC = butcherTable (FEHLBERG_6_4_5 undefined)
  putStrLn $ show resC
  putStrLn $ butcherTableauTex resC

max_stars_repo_name: regnat/hmatrix-sundials

max_issues_repo_name: regnat/hmatrix-sundials

max_forks_repo_name: regnat/hmatrix-sundials


text:
-- By Christian Oliveros and Minmin Chen
module PingPong.Player.Stig.Kinematics where

import Control.Lens (view)
import Data.Bool (bool)
import Data.Geometry hiding (Vector, dot, head)
import Data.Maybe
import GHC.Float
import Numeric.LinearAlgebra hiding (Element)
import PingPong.Model
import PingPong.Player.Stig.General
import PingPong.Player.Stig.GeometryHelpers
import Prelude hiding ((<>))

-- Inverse and Forward Kinematics

-- | Calculate the forward kinematic matrix for a joint / link
calFwdMatrix :: Element -> Element -> [[Float]]
calFwdMatrix (Joint _ a) (Link _ l) =
  [ [cos a, - sin a, l * cos a],
    [sin a, cos a, l * sin a],
    [0, 0, 1]
  ]
calFwdMatrix _ _ = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]

data ArmKinematicPart
  = ArmBase
      { footParam :: R
      }
  | ArmJoint
      { translationParam :: R
      }
  | ArmBat
      { translationParam :: R
      }
  deriving (Show)

type ArmKinematic = [ArmKinematicPart]

-- | Goes from rough representation to a more refined one
fromFootArmToArmKinematic :: Float -> Arm -> ArmKinematic
fromFootArmToArmKinematic foot arm = ArmBase (float2Double foot) : processArm 0 arm
  where
    processArm :: R -> Arm -> ArmKinematic
    processArm t [Link _ t2] = [ArmBat $ t + float2Double t2]
    processArm t (Link _ t2 : as) = processArm (t + float2Double t2) as
    processArm t (Joint _ _ : as) = ArmJoint t : processArm 0 as

-- | Goes from rough representation to a more refined one.
-- As well, get the motion vector
getArmKinematicAndMotion :: Float -> Arm -> (ArmKinematic, Motion)
getArmKinematicAndMotion foot arm =
  (fromFootArmToArmKinematic foot simpleArm, getCurrentJoints simpleArm)
  where
    simpleArm = simplifyArm arm

-- | Converts a motion list to a vector of joint coordinates
motionToJointVector :: Motion -> Vector R
motionToJointVector m = fromList $ map float2Double m

-- | Converts vector of joint coordinates to a motion list
jointVectorToMotion :: Vector R -> Motion
jointVectorToMotion v = map double2Float $ toList v

-- | Translation in the X axis matrix
translateXTrans :: R -> Matrix R
translateXTrans t =
  fromLists [[1, 0, t], [0, 1, 0], [0, 0, 1]] :: Matrix R

-- | Translation in the X axis inverse matrix
translateXInvTrans :: R -> Matrix R
translateXInvTrans t =
  fromLists [[1, 0, - t], [0, 1, 0], [0, 0, 1]] :: Matrix R

-- | Rotation matrix
rotateTrans :: R -> Matrix R
rotateTrans a =
  fromLists [[c, - s, 0], [s, c, 0], [0, 0, 1]] :: Matrix R
  where
    c = cos a
    s = sin a

-- | Rotation inverse matrix
rotateInvTrans :: R -> Matrix R
rotateInvTrans a =
  fromLists [[c, s, 0], [- s, c, 0], [0, 0, 1]] :: Matrix R
  where
    c = cos a
    s = sin a

-- | Create an homogeneous 2D point
homogeneousPoint :: R -> R -> Vector R
homogeneousPoint x y = fromList [x, y, 1]

-- | Zero Point in Homogenous Coordinates
homogeneousZero :: Vector R
homogeneousZero = homogeneousPoint 0 0

-- | 2d Point to Homogeneous Coordinates
pointToHomogenousPoint :: Point 2 Float -> Vector R
pointToHomogenousPoint p = homogeneousPoint (float2Double $ view xCoord p) (float2Double $ view yCoord p)

-- | Create an homogeneous 2D vector
homogeneousVector :: R -> R -> Vector R
homogeneousVector x y = fromList [x, y, 0]

-- | X Vector
homogeneousVectorX :: Vector R
homogeneousVectorX = homogeneousVector 1 0

-- | Y Vector
homogeneousVectorY :: Vector R
homogeneousVectorY = homogeneousVector 0 1

-- | Homogeneous 2D Identity Matrix
homogeneousIdent :: Matrix R
homogeneousIdent = ident 3

data ArmKinematicPartMatrix
  = ArmBaseMatrix
      { footMatrix :: Matrix R
      }
  | ArmJointMatrix
      { translationMatrix :: Matrix R,
        jointMatrix :: Matrix R
      }
  | ArmBatMatrix
      { batMatrix :: Matrix R
      }
  deriving (Show)

type ArmKinematicMatrix = [ArmKinematicPartMatrix]

-- | Get Transform from Part Matrix
getTrans :: ArmKinematicPartMatrix -> Matrix R
getTrans (ArmBaseMatrix t) = t
getTrans (ArmJointMatrix t r) = t <> r
getTrans (ArmBatMatrix t) = t

-- | Get Inverse Transform from Part Matrix
getTransInv :: ArmKinematicPartMatrix -> Matrix R
getTransInv (ArmBaseMatrix t) = t
getTransInv (ArmJointMatrix t r) = r <> t
getTransInv (ArmBatMatrix t) = t

-- | Calculate the transforms for forward kinematics
applyForwardKinematicTrans ::
  ArmKinematic -> Vector R -> ArmKinematicMatrix
applyForwardKinematicTrans arm v = toTrans arm motion
  where
    motion = toList v :: [R]

    toTrans :: ArmKinematic -> [R] -> ArmKinematicMatrix
    toTrans (ArmBase t : as) js = (ArmBaseMatrix $ translateXTrans t <> rotateTrans (pi / 2)) : toTrans as js
    toTrans (ArmJoint t : as) (j : js) = ArmJointMatrix (translateXTrans t) (rotateTrans j) : toTrans as js
    toTrans [ArmBat t] [] = [ArmBatMatrix $ translateXTrans t]

-- | Calculate the inverse transforms for forward kinematics
applyForwardKinematicTransInv ::
  ArmKinematic -> Vector R -> ArmKinematicMatrix
applyForwardKinematicTransInv arm v = reverse $ toTrans arm motion
  where
    motion = toList v :: [R]

    toTrans :: ArmKinematic -> [R] -> ArmKinematicMatrix
    toTrans (ArmBase t : as) js = (ArmBaseMatrix $ rotateInvTrans (pi / 2) <> translateXInvTrans t) : toTrans as js
    toTrans (ArmJoint t : as) (j : js) = ArmJointMatrix (translateXInvTrans t) (rotateInvTrans j) : toTrans as js
    toTrans [ArmBat t] [] = [ArmBatMatrix $ translateXInvTrans t]

-- | Apply Forward Kinematic Transformations
applyForwardKinematicMatrixTrans :: ArmKinematicMatrix -> Matrix R
applyForwardKinematicMatrixTrans = foldr ((<>) . getTrans) homogeneousIdent

-- | Apply Forward Kinematic Inverse Transformations (the list must go from end effecto to base)
applyForwardKinematicMatrixTransInv :: ArmKinematicMatrix -> Matrix R
applyForwardKinematicMatrixTransInv = foldr ((<>) . getTransInv) homogeneousIdent

-- | Compress the transforms of a Forward Kinematic Arm to only the joints and the end effector
compressForwardKinematicsJointsTrans :: ArmKinematicMatrix -> ArmKinematicMatrix
compressForwardKinematicsJointsTrans a = compress a homogeneousIdent
  where
    compress :: ArmKinematicMatrix -> Matrix R -> ArmKinematicMatrix
    compress [ArmBatMatrix t] m = [ArmBatMatrix (m <> t)]
    compress (ArmJointMatrix t r : as) m = ArmJointMatrix (m <> t) r : compress as homogeneousIdent
    compress (a : as) m = compress as (m <> getTrans a)

-- | Compress the transforms of a Forward Kinematic Arm to only the joints and the end base
compressForwardKinematicsJointsInvTrans :: ArmKinematicMatrix -> ArmKinematicMatrix
compressForwardKinematicsJointsInvTrans a = compress a homogeneousIdent
  where
    compress :: ArmKinematicMatrix -> Matrix R -> ArmKinematicMatrix
    compress [ArmBaseMatrix t] m = [ArmBaseMatrix (t <> m)]
    compress (ArmJointMatrix t r : as) m = ArmJointMatrix (t <> m) r : compress as homogeneousIdent
    compress (a : as) m = compress as (getTransInv a <> m)

-- | Calculate the 2D Homogeneous Jacobian of an arm [J^T|0]^T .
-- fwdMatTrans : Kinematic Forward Transforms (From base to end effector).
-- fwdMatTransInv : Kinematic Forward Inverse Transforms (From end effector to base).
-- xLocal : Position of Tool on End Effector Local Coordinates
calculateJacobian ::
  ArmKinematicMatrix -> ArmKinematicMatrix -> Vector R -> Matrix R
calculateJacobian fwdMatTrans fwdMatTransInv xLocal = jH
  where
    -- Compressed form guarantees that the links were applied to the joints
    fwdTransCompressed = compressForwardKinematicsJointsTrans fwdMatTrans
    fwdTransInvCompressed = compressForwardKinematicsJointsInvTrans fwdMatTransInv

    -- Derivative of a revolute Joint
    revoluteDeriv = fromLists [[0, -1, 0], [1, 0, 0], [0, 0, 0]] :: Matrix R

    rollingTransInv :: ArmKinematicMatrix -> Matrix R -> [Matrix R]
    rollingTransInv [ArmBaseMatrix t] m = [t <> m]
    rollingTransInv (a : as) m = roll : rollingTransInv as roll
      where
        roll = getTransInv a <> m
    calculate ::
      ArmKinematicMatrix -> Matrix R -> [Matrix R] -> [Vector R]
    calculate [ArmBatMatrix _] _ [_] = [] -- We reach the bat, we no longer need to process
    calculate (j@(ArmJointMatrix _ _) : as) m (endToJ : mIs) =
      (baseToJ <> revoluteDeriv <> endToJ) #> xLocal : calculate as baseToJ mIs
      where
        baseToJ = m <> getTrans j

    -- Jacobian in this form: [J^T|0]^T there is a row of 0s at the end due to 2d homogeneous coordinates
    jH =
      fromColumns $
        calculate fwdTransCompressed homogeneousIdent (rollingTransInv fwdTransInvCompressed homogeneousIdent)

-- | Get the Jacobian from the Homogeneous Jacobian (Note this doesn't work for the inverse of the Jacobian)
getJacobianFromHomogeneousJacobian :: Matrix R -> Matrix R
getJacobianFromHomogeneousJacobian jH = jH ?? (DropLast 1, All)

-- | Newton Raphson Threshold
newtonRaphsonThreshold :: (Num r, Ord r, Fractional r) => r -> r -> r
newtonRaphsonThreshold = threshold 0.0000001

-- | Calculates a Newton-Raphson IK step (if singular it fails). Returns the new Joints and its error agains target
-- a : Arm Description.
-- q : Current Joints.
-- xLocal : Tool in End-Effector Local Coordinates
-- xTargetGlobal : Target in Global Coordinates
newtonRaphsonStep ::
  ArmKinematic ->
  Vector R ->
  Vector R ->
  Vector R ->
  Maybe (Vector R, R)
newtonRaphsonStep a q xLocal xTargetGlobal
  | singular = Nothing
  | otherwise = Just (qn, norm_2 eN)
  where
    -- Forward Transforms
    fwdT = applyForwardKinematicTrans a q
    fwdTI = applyForwardKinematicTransInv a q
    jH = calculateJacobian fwdT fwdTI xLocal

    -- Bat Global Position
    batGlobal = applyForwardKinematicMatrixTrans fwdT #> xLocal

    -- Error
    e = xTargetGlobal - batGlobal

    -- Drop Homogeneous part
    er = subVector 0 2 e
    j = getJacobianFromHomogeneousJacobian jH

    -- Calculate delta q
    dq = - pinv j #> er

    dqNorm = norm_2 dq

    qnUnormalized = q + dq
    -- Re normalizes angles for more accuracy
    qn = fromList $ map normalizeAngle $ toList qnUnormalized

    -- New Bat Position
    batN = applyForwardKinematicMatrixTrans (applyForwardKinematicTrans a qn) #> xLocal

    -- New Error
    eN = xTargetGlobal - batN

    -- If the move was singular
    singular = rank jH < 2 && newtonRaphsonThreshold 0 dqNorm == 0

-- | Maximum Newton Raphson Step
newtonRaphsonMaxStep :: Int
newtonRaphsonMaxStep = 10000

-- | Really Bad Step
newtonRaphsonBadStep :: forall a. Floating a => a
newtonRaphsonBadStep = 2

-- | When to do a random Restart Newton Raphson Step
newtonRaphsonMaxRandomRestartStep :: Int
newtonRaphsonMaxRandomRestartStep = round (fromIntegral newtonRaphsonMaxStep / 10)

-- | NewtonRaphson Loop Iteration
newtonRaphsonIKIter ::
  Int ->
  Int ->
  ArmKinematic ->
  (Vector R, Vector R) ->
  Vector R ->
  (Vector R, R) ->
  (Vector R, R)
newtonRaphsonIKIter i j a (xLocal, xTargetGlobal) q (qBest, eBest)
  | i >= newtonRaphsonMaxStep = (qBest, eBest)
  | newtonRaphsonThreshold 0 eBest == 0 = (qBest, eBest)
  | singular = newtonRaphsonIKIter (i + 1) 0 a (xLocal, xTargetGlobal) qR (qBest, eBest)
  | needReset = newtonRaphsonIKIter (i + 1) 0 a (xLocal, xTargetGlobal) qR (qBest, eBest)
  | eN < eBest = newtonRaphsonIKIter (i + 1) (j + 1) a (xLocal, xTargetGlobal) qN (qN, eN)
  | otherwise = newtonRaphsonIKIter (i + 1) (j + 1) a (xLocal, xTargetGlobal) qN (qBest, eBest)
  where
    -- Random Vector
    qSize = size q
    pseudoInt = round $ fromIntegral i ** (dot q qBest * bool (eBest + 1) (1 / eBest) (eBest < 1 && 0 < eBest))
    qR = pi / 2 * (2 * randomVector pseudoInt Uniform qSize - 1)

    -- Perform the step
    step = newtonRaphsonStep a q xLocal xTargetGlobal

    singular = case step of
      Nothing -> True
      Just _ -> False

    (qN, eN) = fromJust step

    needReset =
      ( eN >= ((eBest / (1 + fromIntegral j)) * newtonRaphsonBadStep)
          || (newtonRaphsonThreshold 0 (norm_2 (qN - q)) == 0)
      )
        && j >= newtonRaphsonMaxRandomRestartStep

-- | Newton Raphson IK Algorithm
newtonRaphsonIK ::
  ArmKinematic ->
  (Vector R, Vector R) ->
  Vector R ->
  (Vector R, R)
newtonRaphsonIK a (xLocal, xTargetGlobal) q = newtonRaphsonIKIter 0 0 a (xLocal, xTargetGlobal) q (q, eNorm)
  where
    -- Forward Transforms
    fwdT = applyForwardKinematicTrans a q

    -- Bat Global Position
    batGlobal = applyForwardKinematicMatrixTrans fwdT #> xLocal

    -- Error
    e = xTargetGlobal - batGlobal
    eNorm = norm_2 e

-- Calculates newton raphson step to approximate acos
newtonRaphsonAcosStep :: Int -> Int -> Double -> (Double, Double) -> Double -> (Double, Double)
newtonRaphsonAcosStep i j q (qB, eB) t
  | i == newtonRaphsonMaxStep = (qB, eB)
  | newtonRaphsonThreshold 0 eB == 0 = (qB, eB)
  | isSingular || needReset = newtonRaphsonAcosStep (i + 1) 0 qR (qB, eB) t
  | eN < eB = newtonRaphsonAcosStep (i + 1) (j + 1) qN (qN, eN) t
  | otherwise = newtonRaphsonAcosStep (i + 1) (j + 1) qN (qB, eB) t
  where
    -- Calculate derivative and function
    f = t - cos q
    fDeriv = sin q

    -- Check if we are in a singularity
    isSingular = newtonRaphsonThreshold 0 fDeriv == 0

    -- Step
    qN = normalizeAngle $ q - f / fDeriv

    -- New Error
    eN = t - cos qN

    -- Random Restarts
    needReset = eN >= eB * newtonRaphsonBadStep && j >= newtonRaphsonMaxRandomRestartStep

    pseudoInt = round $ fromIntegral i ** (q * qB * bool (eB + 1) (1 / eB) (eB < 1 && 0 < eB))
    randNum = atIndex (randomVector pseudoInt Uniform 1) 0
    qR = pi / 2 * (2 * randNum - 1)

-- | Approximate Acos with newton raphson
newtonRaphsonAcos :: Double -> Double -> (Double, Double)
newtonRaphsonAcos q t = newtonRaphsonAcosStep 0 0 qN (qN, eN) t
  where
    qN = normalizeAngle q
    eN = t - cos qN


text:
{-# LANGUAGE Haskell98 #-}
{-# LINE 1 "src/Data/Complex/Compat.hs" #-}

{-# LANGUAGE CPP, NoImplicitPrelude #-}
module Data.Complex.Compat (
  module Base
) where

import Data.Complex as Base

max_stars_repo_name: BinderDavid/fragnix

max_issues_repo_name: BinderDavid/fragnix

max_forks_repo_name: BinderDavid/fragnix


text:
module Binary where

import Nnaskell
import Numeric.LinearAlgebra.Data

binaryNN :: IO NN
binaryNN = loadNN "./data/binary.txt"

binaryTD :: [(Vector R, Vector R)]
binaryTD = zip (map binaryInput ns) (map binaryOutput ns)
    where ns = [0 .. 15]

binaryInput :: Int -> Vector R
binaryInput x = vector $ map fromIntegral [ x `mod` 2
                                          , (x `div` 2) `mod` 2
                                          , (x `div` 4) `mod` 2
                                          , (x `div` 8) `mod` 2]

binaryOutput :: Int -> Vector R
binaryOutput x = accum (vector $ replicate 16 0.0) (\a _ -> a) [(x, 1.0)]

max_stars_repo_name: tsoding/nnaskell

max_issues_repo_name: tsoding/nnaskell

max_forks_repo_name: tsoding/nnaskell


text:
module D03Spec (spec) where

import Test.Hspec
import D03
import Data.Complex

spec :: Spec
spec = do

  describe "utilities" $ do
    it "parses into claims" $ mkClaims "#1 @ 146,196: 19x14\n#2 @ 641,817: 27x28" `shouldBe` [ Claim 1 (146:+196) (164:+209), Claim 2 (641:+817) (667:+844) ]

  describe "intersections size" $ do
    it "simple example" $ claimsIntersectionSize "#1 @ 1,3: 4x4\n#2 @ 3,1: 4x4\n#3 @ 5,5: 2x2" `shouldBe` 4
    it "other example" $ claimsIntersectionSize "#1 @ 1,2: 3x3\n#2 @ 2,3: 3x4\n#3 @ 1,5: 3x3" `shouldBe` 8
    it "overlap example" $ claimsIntersectionSize "#1 @ 1,2: 3x3\n#2 @ 2,3: 3x4\n#3 @ 1,5: 3x3\n#4 @ 2,5: 2x3" `shouldBe` 10
    it "other overlap example" $ claimsIntersectionSize "#1 @ 1,2: 3x3\n#2 @ 2,3: 3x4\n#3 @ 1,5: 3x3\n#4 @ 2,5: 2x3\n#5 @ 2,2: 3x3" `shouldBe` 14
    it "submits answer" $ do input <- readFile "test/D03.txt"; claimsIntersectionSize input `shouldBe` 106501

  describe "isolated claim" $ do
    it "simple example" $ isolatedClaim "#1 @ 1,3: 4x4\n#2 @ 3,1: 4x4\n#3 @ 5,5: 2x2" `shouldBe` 3
    it "submits answer" $ do input <- readFile "test/D03.txt"; isolatedClaim input `shouldBe` 632


text:
-- file: Astro/Benchmarks.hs
-- Benchmarks for some of the functions

module Astro.Benchmarks
	where

import Astro.Stumpff
import Astro.Kepler
import Astro.Elements

import Numeric.GSL
import Numeric.LinearAlgebra
import Control.Monad

keplerPropagateBM :: Int -> IO ()
keplerPropagateBM nsteps = step 0 (rv,vv)
	where
	(rv, vv) 	= (3|>[1,0,0], 3|>[0,0.1,0])
	oe 		= posvelToOrbel rv vv
	h  		= 10*2*pi/(fromIntegral nsteps::Double)
	step n (r,v) = do
		putStrLn $ (show r)
		(when (n < nsteps) $ step (n+1) newS)
		where
		newS = keplerPropagate r v h (1e-15,100)

max_stars_repo_name: sageh/AstroHaskell

max_issues_repo_name: sageh/AstroHaskell

max_forks_repo_name: sageh/AstroHaskell